## Vector autoregression—simulation, estimation, and inference in Stata

\(\newcommand{\epsb}{{\boldsymbol{\epsilon}}}

\newcommand{\mub}{{\boldsymbol{\mu}}}

\newcommand{\thetab}{{\boldsymbol{\theta}}}

\newcommand{\Thetab}{{\boldsymbol{\Theta}}}

\newcommand{\etab}{{\boldsymbol{\eta}}}

\newcommand{\Sigmab}{{\boldsymbol{\Sigma}}}

\newcommand{\Phib}{{\boldsymbol{\Phi}}}

\newcommand{\Phat}{\hat{{\bf P}}}\)Vector autoregression (VAR) is a useful tool for analyzing the dynamics of multiple time series. VAR expresses a vector of observed variables as a function of its own lags.

**Simulation**

Let’s begin by simulating a bivariate VAR(2) process using the following specification,

\[

\begin{bmatrix} y_{1,t}\\ y_{2,t}

\end{bmatrix}

= \mub + {\bf A}_1 \begin{bmatrix} y_{1,t-1}\\ y_{2,t-1}

\end{bmatrix} + {\bf A}_2 \begin{bmatrix} y_{1,t-2}\\ y_{2,t-2}

\end{bmatrix} + \epsb_t

\]

where \(y_{1,t}\) and \(y_{2,t}\) are the observed series at time \(t\), \(\mub\) is a \(2 \times 1\) vector of intercepts, \({\bf A}_1\) and \({\bf A}_2\) are \(2\times 2\) parameter matrices, and \(\epsb_t\) is a \(2\times 1\) vector of innovations that is uncorrelated over time. I assume a \(N({\bf 0},\Sigmab)\) distribution for the innovations \(\epsb_t\), where \(\Sigmab\) is a \(2\times 2\) covariance matrix.

I set my sample size to 1,100 and generate variables to hold the observed series and innovations.

. clear all . set seed 2016 . local T = 1100 . set obs `T' number of observations (_N) was 0, now 1,100 . gen time = _n . tsset time time variable: time, 1 to 1100 delta: 1 unit . generate y1 = . (1,100 missing values generated) . generate y2 = . (1,100 missing values generated) . generate eps1 = . (1,100 missing values generated) . generate eps2 = . (1,100 missing values generated)

In lines 1–6, I set the seed for the random-number generator, set my sample size to 1,100, and generate a time variable, **time**. In the remaining lines, I generate variables **y1**, **y2**, **eps1**, and **eps2** to hold the observed series and innovations.

**Setting parameter values**

I choose the parameter values for the VAR(2) model as follows:

\[

\mub = \begin{bmatrix} 0.1 \\ 0.4 \end{bmatrix}, \quad {\bf A}_1 =

\begin{bmatrix} 0.6 & -0.3 \\ 0.4 & 0.2 \end{bmatrix}, \quad {\bf A}_2 =

\begin{bmatrix} 0.2 & 0.3 \\ -0.1 & 0.1 \end{bmatrix}, \quad \Sigmab =

\begin{bmatrix} 1 & 0.5 \\ 0.5 & 1 \end{bmatrix}

\]

. mata: ------------------------------------------------- mata (type end to exit) ----- : mu = (0.1\0.4) : A1 = (0.6,-0.3\0.4,0.2) : A2 = (0.2,0.3\-0.1,0.1) : Sigma = (1,0.5\0.5,1) : end -------------------------------------------------------------------------------

In Mata, I create matrices **mu**, **A1**, **A2**, and **Sigma** to hold the parameter values. Before generating my data, I check whether these values correspond to a stable VAR(2) process. Let

\[

{\bf F} = \begin{bmatrix} {\bf A}_1 & {\bf A_2} \\ {\bf I}_2 & {\bf 0}

\end{bmatrix}

\]

denote a \(4\times 4\) matrix where \({\bf I}_2\) is a \(2\times 2\) identity matrix and \({\bf 0}\) is a \(2\times 2\) matrix of zeros. The VAR(2) process is stable if the modulus of all eigenvalues of \({\bf F}\) is less than 1. The code below computes the eigenvalues.

. mata: ------------------------------------------------- mata (type end to exit) ----- : K = p = 2 // K = number of variables; p = number of lags : F = J(K*p,K*p,0) : F[1..2,1..2] = A1 : F[1..2,3..4] = A2 : F[3..4,1..2] = I(K) : X = L = . : eigensystem(F,X,L) : L' 1 +----------------------------+ 1 | .858715598 | 2 | -.217760515 + .32727213i | 3 | -.217760515 - .32727213i | 4 | .376805431 | +----------------------------+ : end -------------------------------------------------------------------------------

I construct the matrix **F** as defined above and use the function **eigensystem()** to compute its eigenvalues. The matrix **X** holds the eigenvectors, and **L** holds the eigenvalues. All eigenvalues in **L** are less than 1. The modulus of the second and third complex eigenvalues is \(\sqrt{r^2 + c^2} = 0.3931\), where \(r\) is the real part and \(c\) is the complex part. After checking the stability condition, I generate draws for the VAR(2) model.

**Drawing innovations from a multivariate normal**

I draw random normal values from \(N({\bf 0},\Sigmab)\) and assign them to Stata variables **eps1** and **eps2**.

. mata: ------------------------------------------------- mata (type end to exit) ----- : T = strtoreal(st_local("T")) : u = rnormal(T,2,0,1)*cholesky(Sigma) : epsmat = . : st_view(epsmat,.,"eps1 eps2") : epsmat[1..T,.] = u : end -------------------------------------------------------------------------------

I assign the sample size, defined in Stata as a local macro **T**, to a Mata numeric variable. This simplifies changing the sample size in the future; I only need to do this once at the beginning. In Mata, I use two functions: **st_local()** and **strtoreal()** to assign the sample size. The first function obtains strings from Stata macros, and the second function converts them into a real value.

The second line draws a \(1100 \times 2\) matrix of normal errors from a \(N({\bf 0},\Sigmab)\) distribution. I use the **st_view()** function to assign the draws to the Stata variables **eps1** and **eps2**. This function creates a matrix that is a view on the current Stata dataset. I create a null matrix **epsmat** and use **st_view()** to modify **epsmat** based on the values of the Stata variables **eps1** and **eps2**. Finally, I assign this matrix to hold the draws stored in **u**, effectively populating the Stata variables **eps1** and **eps2** with the random draws.

**Generating the observed series**

Following Lütkepohl (2005, 708), I generate the first two observations so that their correlation structure is the same as the rest of the sample. I assume a bivariate normal distribution with mean equal to the unconditional mean \(\thetab = ({\bf I}_K – {\bf A}_1 – {\bf A}_2)^{-1}\mub\). The covariance matrix of the first two observations of the two series is

\[

\text{vec}(\Sigmab_y) = ({\bf I}_{16} – {\bf F} \otimes {\bf F})^{-1}\text{vec}(\Sigmab_{\epsilon})

\]

where \(\text{vec}()\) is an operator that stacks matrix columns, \({\bf I}_{16}\) is a \(16\times 16\) identity matrix, and \(\Sigmab_{\epsilon} = \begin{bmatrix} \Sigmab & {\bf 0}\\ {\bf 0} & {\bf 0} \end{bmatrix}\) is a \(4\times 4\) matrix. The first two observations are generated as

\[

\begin{bmatrix} {\bf y}_0 \\ {\bf y}_{-1} \end{bmatrix} = {\bf Q} \etab + \Thetab

\]

where **Q** is a \(4\times 4\) matrix such that \({\bf QQ}’ = \Sigmab_y\), \(\etab\) is a \(4\times 1\) vector of standard normal innovations and \(\Thetab = \begin{bmatrix} \thetab \\ \thetab \end{bmatrix}\) is a \(4\times 1\) vector of means.

The following code generates the first two observations and assigns the values to the Stata variables **y1** and **y2**.

. mata: ------------------------------------------------- mata (type end to exit) ----- : Sigma_e = J(K*p,K*p,0) : Sigma_e[1..K,1..K] = Sigma : Sigma_y = luinv(I((K*p)^2)-F#F)*vec(Sigma_e) : Sigma_y = rowshape(Sigma_y,K*p)' : theta = luinv(I(K)-A1-A2)*mu : Q = cholesky(Sigma_y)*rnormal(K*p,1,0,1) : data = . : st_view(data,.,"y1 y2") : data[1..p,.] = ((Q[3..4],Q[1..2]):+mu)' : end -------------------------------------------------------------------------------

After generating the first two observations, I can generate the rest of the series from a VAR(2) process in Stata as follows:

. forvalues i=3/`T' { 2. qui { 3. replace y1 = 0.1 + 0.6*l.y1 - 0.3*l.y2 + 0.2*l2.y1 + 0.3*l2.y2 + eps1 in ` > i' 4. replace y2 = 0.4 + 0.4*l.y1 + 0.2*l.y2 - 0.1*l2.y1 + 0.1*l2.y2 + eps2 in ` > i' 5. } 6. } . drop in 1/100 (100 observations deleted)

I added the **quietly** statement to suppress the output generated after the **replace** command. Finally, I drop the first 100 observations as burn-in to mitigate the effect of initial values.

**Estimation**

I use the **var** command to fit a VAR(2) model.

. var y1 y2 Vector autoregression Sample: 103 - 1100 Number of obs = 998 Log likelihood = -2693.949 AIC = 5.418735 FPE = .7733536 HQIC = 5.43742 Det(Sigma_ml) = .7580097 SBIC = 5.467891 Equation Parms RMSE R-sq chi2 P>chi2 ---------------------------------------------------------------- y1 5 1.14546 0.5261 1108.039 0.0000 y2 5 .865602 0.4794 919.1433 0.0000 ---------------------------------------------------------------- ------------------------------------------------------------------------------ | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- y1 | y1 | L1. | .5510793 .0324494 16.98 0.000 .4874797 .614679 L2. | .2749983 .0367192 7.49 0.000 .20303 .3469667 | y2 | L1. | -.3080881 .046611 -6.61 0.000 -.3994439 -.2167323 L2. | .2551285 .0425803 5.99 0.000 .1716727 .3385844 | _cons | .1285357 .0496933 2.59 0.010 .0311387 .2259327 -------------+---------------------------------------------------------------- y2 | y1 | L1. | .3890191 .0245214 15.86 0.000 .340958 .4370801 L2. | -.0190324 .027748 -0.69 0.493 -.0734175 .0353527 | y2 | L1. | .1944531 .035223 5.52 0.000 .1254172 .263489 L2. | .0459445 .0321771 1.43 0.153 -.0171215 .1090106 | _cons | .4603854 .0375523 12.26 0.000 .3867843 .5339865 ------------------------------------------------------------------------------

By default, **var** estimates the parameters of a VAR model with two lags. The parameter estimates are significant and are similar to the true values used to generate the bivariate series.

**Inference: Impulse response functions**

Impulse response functions (IRF) are useful to analyze the response of endogenous variables in the VAR model due to an exogenous impulse to one of the innovations. For example, in a bivariate VAR of inflation and interest rate, IRFs can trace out the response of interest rates over time due to exogenous shocks to the inflation equation.

Let’s consider the bivariate model that I used earlier. Suppose I want to estimate the effect of a unit shock in \(\epsb_t\) to the endogenous variables in my system. I do this by converting the VAR(2) process into an MA(\(\infty\)) process

as

\[

{\bf y}_t = \sum_{p=0}^\infty \Phib_p \mub + \sum_{p=0}^\infty \Phib_p \epsb_{t-p}

\]

where \(\Phib_p\) is the MA(\(\infty\)) coefficient matrices for the *p*th lag. Following Lütkepohl (2005, 22–23), the MA coefficient matrices are related to the AR matrices as follows:

\begin{align*}

\Phib_0 &= {\bf I}_2 \\

\Phib_1 &= {\bf A}_1 \\

\Phib_2 &= \Phib_1 {\bf A}_1 + {\bf A}_2 = {\bf A}_1^2 + {\bf A}_2\\

&\vdots \\

\Phib_i &= \Phib_{i-1} {\bf A}_1 + \Phib_{i-2}{\bf A}_2

\end{align*}

More formally, the future \(h\) responses of the *i*th variable due to a unit shock in the *j*th equation at time \(t\) are given by

\[

\frac{\partial y_{i,t+h}}{\partial \epsilon_{j,t}} = \{\Phib_h\}_{i,j}

\]

The responses are the collection of elements in the *i*th row and *j*th column of the MA(\(\infty\)) coefficient matrices. For the VAR(2) model, the first few responses using the estimated AR parameters are as follows:

\begin{align*}

\begin{bmatrix} y_{1,0} \\ y_{2,0} \end{bmatrix} &=

\Phib_0 = {\bf I}_2 =

\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix} \\

\begin{bmatrix} y_{1,1} \\ y_{2,1} \end{bmatrix} &=

\hat{\Phib}_1 = \hat{{\bf A}}_1 =

\begin{bmatrix} 0.5510 & -0.3081 \\ 0.3890 & 0.1945 \end{bmatrix} \\

\begin{bmatrix} y_{1,2} \\ y_{2,2} \end{bmatrix} &= \hat{\Phib}_2 =

\hat{{\bf A}}_1^2 + \hat{{\bf A}}_2 =

\begin{bmatrix} 0.4588 & 0.0254 \\ 0.2710 & -0.0361 \end{bmatrix}

\end{align*}

Future responses of \({\bf y}_t\) for \(t>2\) is computed using similar recursions. The response of the first variable due to an impulse in the first equation is the vector \((1,0.5510,0.4588,\dots)\). I obtain the impulse responses using the **irf create** command as follows:

. irf create firstirf, set(myirf) (file myirf.irf created) (file myirf.irf now active) (file myirf.irf updated)

This command estimates the IRFs and other statistics in **firstirf** and saves the file as **myirf.irf**. The **set()** option sets the filename **myirf.irf** as the active file. I can list the table of responses of **y1** due to an impulse in the same equation by typing

. irf table irf, impulse(y1) response(y1) noci Results from firstirf +--------------------+ | | (1) | | step | irf | |--------+-----------| |0 | 1 | |1 | .551079 | |2 | .458835 | |3 | .42016 | |4 | .353356 | |5 | .305343 | |6 | .263868 | |7 | .227355 | |8 | .196142 | +--------------------+ (1) irfname = firstirf, impulse = y1, and response = y1

The default horizon is 8, and I specify the **noci** option to suppress the confidence limits. Notice that the first few responses are similar to the ones I computed earlier. IRFs are better analyzed using a graph. I can visualize this in a graph along with the 95% confidence bands by typing

. irf graph irf, impulse(y1) response(y1)

A unit impulse to the first equation increases **y1** by a unit

contemporaneously. The response of **y1** slowly declines over time to its long-run level.

**Orthogonalized impulse–response functions**

In the previous section, I showed the response of **y1** due to a unit impulse on the same equation while holding the other impulse constant. The variance–covariance matrix \(\Sigmab\), however, implies a strong positive correlation between the two equations. I list the contents of the estimated variance–covariance matrix by typing

. matrix list e(Sigma) symmetric e(Sigma)[2,2] y1 y2 y1 1.3055041 y2 .4639629 .74551376

The estimated covariance of the two equations is positive. This implies that I cannot assume the other impulse will remain constant. An impulse to the **y2** equation has a contemporaneous effect on **y1**, and vice versa.

Orthogonalized impulse-response functions (OIRF) address this by decomposing the estimated variance–covariance matrix \(\hat{\Sigmab}\) into a lower triangular matrix. This type of decomposition isolates the contemporaneous response of **y1** arising solely because of an impulse in the same equation. Nonetheless, the impulse on the first equation will still contemporaneously affect **y2**. For example, if **y1** is inflation and **y2** is interest rate, this decomposition implies that shock to inflation affects both inflation and interest rates. However, a shock to the interest rate equation only affects interest rates.

To estimate the OIRFs, let \({\bf P}\) denote the Cholesky decomposition of \(\Sigmab\) such that \({\bf P} {\bf P}’=\Sigmab\). Let \({\bf u}_t\) denote a \(2\times 1\) vector such that \({\bf P u}_t = \epsb_t\), which also implies \({\bf

u}_t = {\bf P}^{-1} \epsb_t\). The errors in \({\bf u}_t\) are uncorrelated by construction because \(E({\bf u}_t{\bf u}_t’) = {\bf P}^{-1}E(\epsb_t \epsb_t’) {\bf P’}^{-1}={\bf I}_2\). This allows us to interpret the OIRFs as a one standard-deviation impulse to \({\bf u}_t\).

I rewrite the MA(\(\infty\)) representation from earlier in terms of the \({\bf u}_t\) vector as

\[

{\bf y}_t = \sum_{p=0}^\infty \Phib_p \mub + \sum_{p=0}^\infty \Phib_p {\bf P u}_{t-p}

\]

The OIRFs are simply

\[

\frac{\partial y_{i,t+h}}{\partial u_{j,t}} = \{\Phib_h {\bf

P}\}_{i,j}

\]

the product of the MA coefficient matrices and the lower triangular matrix **P**.

I obtain the estimate of \(\Phat\) by typing

. matrix Sigma_hat = e(Sigma) . matrix P_hat = cholesky(Sigma_hat) . matrix list P_hat P_hat[2,2] y1 y2 y1 1.1425866 0 y2 .40606367 .76198823

Using this matrix, I compute the first few responses as follows:

\begin{align*}

\begin{bmatrix} y_{1,0} \\ y_{2,0} \end{bmatrix} &=

\Phib_0 \Phat = {\bf I}_2 \Phat =

\begin{bmatrix} 1.1426 & 0\\ 0.4061 & 0.7620 \end{bmatrix} \\

\begin{bmatrix} y_{1,1} \\ y_{2,1} \end{bmatrix} &=

\hat{\Phib}_1 \Phat = \hat{{\bf A}}_1 \Phat =

\begin{bmatrix} 0.5046 & -0.2348 \\ 0.5234 & 0.1482 \end{bmatrix} \\

\begin{bmatrix} y_{1,2} \\ y_{2,2} \end{bmatrix} &= \hat{\Phib}_2 \Phat

= (\hat{{\bf A}}_1^2 + \hat{{\bf A}}_2) \Phat =

\begin{bmatrix} 0.5346 & 0.0194 \\ 0.2950 & -0.0275 \end{bmatrix}

\end{align*}

I list all the OIRFs in a table and plot the response of **y1** due to an impulse in **y1**.

. irf table oirf, noci Results from firstirf +--------------------------------------------------------+ | | (1) | (2) | (3) | (4) | | step | oirf | oirf | oirf | oirf | |--------+-----------+-----------+-----------+-----------| |0 | 1.14259 | .406064 | 0 | .761988 | |1 | .504552 | .523448 | -.23476 | .148171 | |2 | .534588 | .294977 | .019384 | -.027504 | |3 | .476019 | .279771 | -.0076 | .013468 | |4 | .398398 | .242961 | -.010024 | -.00197 | |5 | .346978 | .206023 | -.003571 | -.003519 | |6 | .299284 | .178623 | -.004143 | -.001973 | |7 | .257878 | .154023 | -.003555 | -.002089 | |8 | .222533 | .13278 | -.002958 | -.001801 | +--------------------------------------------------------+ (1) irfname = firstirf, impulse = y1, and response = y1 (2) irfname = firstirf, impulse = y1, and response = y2 (3) irfname = firstirf, impulse = y2, and response = y1 (4) irfname = firstirf, impulse = y2, and response = y2

**irf table oirf** requests the OIRFs. Notice that the estimates in the first three rows are the same as the ones computed above.

. irf graph oirf, impulse(y1) response(y1)

The graph corresponds to the response of **y1** as a result of a one standard-deviation impulse to the same equation.

**Conclusion**

In this post, I showed how to simulate data from a stable VAR(2) model. I estimated the parameters of this model using the **var** command. I showed how to estimate IRFs and OIRFs. The latter obtains the responses using the lower triangular decomposition of the covariance matrix.

**Reference**

Lütkepohl, H. 2005. *New Introduction to Multiple Time Series Analysis*. New York: Springer.