I suggest to check the following links which are very useful for you.

https://bookdown.org/singh_pratap_tejendra/intro_time_series_r/multivariate-ts-analysis.html

https://stackoverflow.com/questions/1714280/multivariate-time-series-modelling-in-r

You can check this site

https://faculty.chicagobooth.edu/ruey.tsay/teaching/ama/

The tsDyn package offers basic functionalities for the estimation of Vector Autoregressive models with possible regime-switching behavior and or co-integration. There is also a dedicated vars package. If you care about modeling volatility spillovers as well, then rmgarch allows you to estimate VAR structures along with multivariate autoregressive conditional heteroscedastic behavior.

From the details you provide it is not immediately clear what question you are trying to solve and what a sensible approach is to obtain a fair result. Since you refer to your other covariates as "independent", it seems that you initial hunch is to not model the multivariate behavior but simply focus on the temporal dynamics of the key variable of interest, along with possible dependencies on other covariates. That can be achieved using basic lm or AR functionalities that are available in the standard stats package.

You mention that the assumptions of OLS are met; note that the assumptions behind LLN and CLT in a time-series context differ from the usual context in which the I.I.D. assumption is imposed. In particular, when series are dependent, for example in the temporal or spatial dimension, then, for example, the LLN relies on stable behavior of the underlying process and the CLT relies on suitable forms of weak dependence. This in turn may be implied by reasonably well fitting models under certain ranges of possible parameter values. Often it is the case that assumptions are stricter for more complex models, and many dynamic models will rely on more general likelihood or moment approaches.

Finally, you did not mention the frequency of the data. Given that you have a total of 13 variables, you may, if the frequency of the data is low, need to confine your analysis to a subset of key variables of interest if you interest is in modeling the entire dynamical system. Otherwise you may not have enough observations to estimate a multivariate time series structure. In particular, as the the number of lags increases, the number of required parameters will sharply increase when you are considering many variables.

Good luck with you analysis, I hope this helps.

Hi Bo Andree, thank you for your explanations. They are very useful.

Hi, Antoine

The tsDyn package would be very helpful in your study. It provides the tools to estimate time series behavior, in particular non-linear.

Bo Andree , I collected data on some Burundian economic indicators from World Bank official web site. The literature on the drivers of economic growth shows some factors that affect economic growth. My aim is to test those factors the way they affect Burundian economic growth. So, collected data covers the period between 1980 to 2017 which means I have 38 observations.

Since those data are related to time, I would like to know how to make a multivariate time series.

Niyungeko Antoine , 38 observations is a very limited information set. In this case it is best not to use complex procedures and instead focus only on simple methods that may still be enough to learn interesting properties of the data. For example, you might want to tabulate the data in useful ways and compare trends along with a simple linear regression. This may help support some theories that you may already have based on the literature you mention. Either way, it will be important to be modest in your conclusions as typically you will need more data to do any substantive time series analysis.

If you insist on multivariate modeling, you can look into Bayesian Multivariate Autoregressive Time Series Models (or penalized estimators if you want to stay within frequentist frameworks). These procedures may have better small sample properties. In this case it may be important to provide some simulation evidence relevant to your exact data dimensions and parameter complexity. This is a big hike in complexity in contrast to my other suggestion so you may want to explore your data a bit more before deciding to invest in this. Good luck.

Hi Serei, apparently tsdyn is for univariate analysis according to to the literature.

Bo Andree ,thank you so much for your advice.

OLS can only be used if all the variables are stationary at level I.e I(0)

Toun, according to you, when all variables are stationary I can us OLS?Right?

toun, I do not agree because using OLS leads to not take account the fact that data have autoregressive part and moving average part.

I doubt it would be possible to take care of an autoregressive part and a moving average part with only 38 observations of one variable and 12 independent variables. And that all these 13 variables have a stationary behaviour. Look in the literature for what others do in that case.

tsDyn includes both univariate and multivariate procedures. The description given on CRAN ( https://cran.r-project.org/web/packages/tsDyn/index.html) is as follows -

Implements nonlinear autoregressive (AR) time series models. For univariate series, a non-parametric approach is available through additive nonlinear AR. Parametric modeling and testing for regime switching dynamics is available when the transition is either direct (TAR: threshold AR) or smooth (STAR: smooth transition AR, LSTAR). For multivariate series, one can estimate a range of TVAR or threshold cointegration TVECM models with two or three regimes. Tests can be conducted for TVAR as well as for TVECM (Hansen and Seo 2002 and Seo 2006).

Version:10-1.1