Last lecture:

• Cointegration: concepts and testing

• Matrix representation of AR: VAR

This lecture:

• Identify the shocks from a VAR system

• Stationarity and Stability of VAR

• Examples

• \(n\) entries in the column vector: \[ \mathbf {e} ={\begin{bmatrix}e_{1}\\\vdots \\e_{n}\end{bmatrix}}\] are random variables, each with finite variance, then the covariance matrix \(\Sigma\) is the matrix whose \((i,j)\) entry is the covariance

\[\sigma _{ij}=\mbox{Cov} (e_{i},e_{j})=\mathbb{E} [(e_{i}-\mathbb{E}[e_{i}])(e_{j}-\mathbb{E}[e _{j}])]\]

• As \(\mathbb{E} [(e_{1}-\mathbb{E}[e _{1}])(e_{1}-\mathbb{E}[e _{1}])] = \sigma_{11}\), \(\mathbb{E} [(e_{1}-\mathbb{E}[e _{1}])(e_{n}-\mathbb{E}[e _{n}])]=\sigma_{1n}\)

and the covariance matrix is

\[\Sigma ={\begin{bmatrix}\sigma_{11}&\cdots &\sigma_{1n}\\\\\vdots &\ddots &\vdots \\\\\sigma_{n1}&\cdots &\sigma_{nn}\end{bmatrix}}.\]

• The pair \(\{e_{1,t}, e_{2,t} \}\) are multivariate white noise if each individual series is white noise and the CCF for all \(k\neq 0\) \[\gamma_{e_1, e_2}(k)=\mbox{Cov}(e_{1,t}, e_{2,t+k}) = 0,\,\, \gamma_{e_1,e_2}(0)=\Sigma.\]

library(mvtnorm)
covmat = matrix(c(1,0.5,0.5,1) , nr=2)
e = rmvnorm(1000, sigma = covmat)
cov(e)

##           [,1]      [,2]
## [1,] 1.0060044 0.5006652
## [2,] 0.5006652 0.9768115

• Recall the VAR(1) form: \[\underset{(n \times1)}{\mathbf{Y}_{t}}=\underset{(n\times n)}{\mathbf{\Phi}}\underset{(n\times1)}{\mathbf{Y}_{t-1}}+\underset{(n\times n)}{\mathbf{C}} \underset{(n\times1)}{\mathbf{w}_t}\] Note \[\underset{(n \times n)}{\mathbf{C}} \underset{(n\times1)}{\mathbf{w}_t} = \underset{(n\times1)}{\mathbf{e}_t} , \, \mathbf{C} \mathbb{E}[\mathbf{w}_t \mathbf{w}_t^{T}] \mathbf{C}^{T} = \underset{(n\times n)}{\mathbf{C} \mathbf{C}^{T}} = \mbox{Var}(\mathbf{e}_t) = \underset{(n \times n)}{ \Sigma}\] where \(\mathbf{w}_t\) is the “standardized” multivariate white noise \(\mbox{Var}(\mathbf{w}_t) = \mathbf{I}\).

• \(\Sigma\) is a real symmetric square matrix.

• The reaction of the system following a shock.

• The IRF is the path that \(Y\) follows if it is kicked by a single unit shock, i.e., \(Y_{t−j}= 0\), \(Y_t = 1\), \(Y_{t+j} = 0\).

• Causes and effects: the cause of one unit shock leads to its effect on the whole system.

• For example, the effective change of GNP to a shock in M1. VAR interprets the result as the effect on GNP of monetary policy that includes total supply of money, interest rates, exchange rates, equities and bonds, etc.

Y = rep(0, 50); Y[25] = 1; C = rep(0,50) # two variables
for (i in 25:50){ Y[i+1] = 0.8* Y[i]
  C[i+1] = - 0.9 * Y[i]} 
ts.plot(C, ylim=c(-1,1), col="red"); lines(Y, col="blue")

Y = rep(0, 50); Y[25] = 1; C = rep(0,50) # two variables
for (i in 25:50){ Y[i+1] = 0.8* Y[i]
  C[i+1] = -0.9 * Y[i] - 0.9*C[i] } 
ts.plot(C, ylim=c(-1,1), col="red"); lines(Y, col="blue")

Y = rep(0, 50); Y[25] = 1; C = rep(0,50) # two variables
for (i in 25:50){ Y[i+1] =  0.8* Y[i]
  C[i+1] = 0.9 * Y[i] + 0.9*C[i] } 
ts.plot(C, ylim=c(0,3), col="red"); lines(Y, col="blue")

Y = rep(0, 50); Y[25] = 1; C = rep(0,50) # two variables
for (i in 25:50){ Y[i+1] =  0.8* Y[i]
  C[i+1] = 0.9 * Y[i] - 0.9*C[i] } 
ts.plot(C, ylim=c(-0.3,1), col="red"); lines(Y, col="blue")

Y = rep(0, 50); Y[25] = 1; C = rep(0,50) # two variables
for (i in 25:50){ Y[i+1] =  0.8* Y[i]
  C[i+1] = - 0.9 * Y[i] + 0.9*C[i] } 
ts.plot(C, ylim=c(-3,1), col="red"); lines(Y, col="blue")

Y = rep(0, 50); Y[25] = 1; C = rep(0,50) # two variables
for (i in 25:50){ Y[i+1] =  0.8* Y[i] + 0.1*C[i]
  C[i+1] = - 0.9 * Y[i] + 0.9*C[i] } 
ts.plot(C, ylim=c(-4,4), col="red"); lines(Y, col="blue")

• Understand how does a shock affect each component of the system.

• IRF for all component of a VAR(1) \[\underset{(n \times1)}{\mathbf{Y}_{t}}=\underset{(n\times n)}{\mathbf{\Phi}}\underset{(n\times1)}{\mathbf{Y}_{t-1}}+\underset{(n\times n)}{\mathbf{C}} \underset{(n\times1)}{\mathbf{w}_t}= \underset{(n\times n)}{\mathbf{\Phi}}\underset{(n\times1)}{\mathbf{Y}_{t-1}}+\underset{(n\times1)}{\mathbf{e}_t}\] can be written as \[ \begin{align} (\mathbf{I} - \mathbf{\Phi} \mathbb{L}) \mathbf{Y}_{t} = & \Phi(\mathbb{L}) \mathbf{Y}_{t} = \mathbf{e}_t \\ \mathbf{Y}_{t} = & \Phi(\mathbb{L})^{-1} \mathbf{e}_t = \Phi(\mathbb{L})^{-1}\mathbf{C} \mathbf{w}_t. \end{align}\]

• ISSUE: Several invertible \(\mathbf{C}\) can be “planted” here. It leads to an undetermined linear combination \(\mathbf{C} \mathbf{w}_t\) for the shocks.

• Orthogonalization: construct the orthogonal shocks that are uncorrelated.

• Construct the matrix \(\mathbf{C}\) via the Cholesky decomposition: \[\Sigma = \mathbf{C} \mathbf{C}^{T} \] where \(\mathbf{C}\) is a triangular matrix.

• Orthogonal shock is importnat. If the M1 shock is correlated with the GNP shock, then the observed response of GNP to money can also come from a investment shock (to GNP) that happened at the same time (e.g. the central bank sees the GNP shock and accommodates it by using monetary policy).

A = as.matrix(data.frame(c(3,4,3),c(4,8,6),c(3,6,9)))
A.chol = chol(A); A.chol

##      c.3..4..3. c.4..8..6. c.3..6..9.
## [1,]   1.732051   2.309401   1.732051
## [2,]   0.000000   1.632993   1.224745
## [3,]   0.000000   0.000000   2.121320

t(A.chol) %*% A.chol

##            c.3..4..3. c.4..8..6. c.3..6..9.
## c.3..4..3.          3          4          3
## c.4..8..6.          4          8          6
## c.3..6..9.          3          6          9

Cholesky decompsition satisfies:

• the errors are orthogonal so that the instantaneous response of one variable to the other shock is zero.

• For IRF, \[\mathbf{Y}_{t} = \Phi(\mathbb{L})^{-1}\mathbf{C} \mathbf{w}_t = \Phi(\mathbb{L})^{-1} \begin{pmatrix} C_{11} & 0 & 0 \\ C_{21} & C_{22} & 0 \\ \cdots & \cdots & \cdots\\ \end{pmatrix} \begin{pmatrix} w_{1} \\ w_{2} \\ \cdots \\ \end{pmatrix}\] we have one shock \(w_1\) only affects one variable through \(C_{11}\). So this contemporaneous shock through \(C_{11}\) does not affect the other variables in the equation.

library(Quandl); library(xts); library(lubridate); library(forecast)
library(lattice);library(zoo)
BE.econIndicators = c("GDP Per Capita" = "WWDI/BEL_NY_GDP_PCAP_KN",
                  "GDP Per Capita Growth" = "WWDI/BEL_NY_GDP_PCAP_KD_ZG",
                  "Real Interest Rate" = "WWDI/BEL_FR_INR_RINR",
                  "Exchange Rate" = "WWDI/BEL_PX_REX_REER",
                  "Inflation" = "WWDI/BEL_FP_CPI_TOTL_ZG",
                  "Labor Force Part Rate" = "WWDI/BEL_SL_TLF_ACTI_ZS")
BE = Quandl(BE.econIndicators, type="xts", start_date="1962-01-01",
            end_date="2017-01-01")
colnames(BE)= names(BE.econIndicators)

v2 = data.frame(BE$"Inflation", BE$"GDP Growth Rate")
acf(v2)

• Consider VAR(1) \[\left(\begin{array}{c} Y_{1,t}\\ Y_{2,t} \end{array}\right)= \left(\begin{array}{c} c_{1}\\ c_{2} \end{array}\right) + \left(\begin{array}{cc} \phi_{11} & \phi_{12}\\ \phi_{21} & \phi_{22} \end{array}\right)\left(\begin{array}{c} Y_{1,t-1}\\ Y_{2,t-1} \end{array}\right)+\left(\begin{array}{c} \varepsilon_{1,t}\\ \varepsilon_{2,t} \end{array}\right)\]

• This is equivalent to \[ \begin{align} Y_{1,t}=& c_{1}+\phi_{11}Y_{1,t-1}+\phi_{12}Y_{2,t-1}+\varepsilon_{1,t}\\ Y_{2,t}=& c_{2}+\phi_{21}Y_{1,t-1}+\phi_{22}Y_{2,t-1}+\varepsilon_{2,t} \end{align} \] where \(\mbox{Cov}(\varepsilon_{1,t},\varepsilon_{2,s})=\sigma_{12}\) for \(t=s\), otherwise \(\mbox{Cov}(\varepsilon_{1,t},\varepsilon_{2,s})=0\).

library(vars); library(forecast); library(tseries)
BE.VAR = VAR(v2, p=1); coef(BE.VAR)

## $Inflation
##                       Estimate Std. Error     t value     Pr(>|t|)
## Inflation.l1        0.80847642 0.07175584 11.26704723 1.891446e-15
## GDP.Growth.Rate.l1  0.32459649 0.09860020  3.29204700 1.810241e-03
## const              -0.01058949 0.38817258 -0.02728037 9.783426e-01
## 
## $GDP.Growth.Rate
##                       Estimate Std. Error    t value    Pr(>|t|)
## Inflation.l1       -0.06776988 0.09594399 -0.7063484 0.483186082
## GDP.Growth.Rate.l1  0.29807659 0.13183730  2.2609428 0.028058695
## const               1.74568679 0.51902153  3.3634189 0.001467306

var.resid = resid(BE.VAR)
acf(var.resid)

plot(var.resid)

ts.plot(var.resid[,1])
lines((var.resid[,2]), col="red")

• The system can identify the shocks \[ \begin{align} \mbox{Inflation}_{t}=& -0.01 + 0.32 \mbox{Growth}_{t-1}+ 0.8 \mbox{Inflation}_{t-1}+\varepsilon_{1,t} \\ \mbox{Growth}_{t}=& 1.75 + 0.3 \mbox{Growth}_{t-1} +\varepsilon_{2,t} \end{align} \]

• The coefficient matrix \(\mathbf{\Phi}\) is a upper triangular matrix.

• The shocks of the growth will give impact to inflation, while inflation cannot obstruct the growth in this empirical study.

irf.inf = irf(BE.VAR, impulse = "Inflation", 
                 n.ahead = 20, ci = 0.95)
plot(irf.inf)

irf.growth = irf(BE.VAR, impulse = "GDP.Growth.Rate", 
                 n.ahead = 20, ci = 0.95)
plot(irf.growth)

• Recall the stability and stationarity requirement for AR(1) \(\lim_{j\rightarrow\infty}\phi^{j}=0\) when \(|\phi|<1.\)

• Similarly for \(n\)-Vector AR(1) \[\lim_{j\rightarrow\infty} \det \mathbf{\Phi}^{j}=0\] where the determinant is \(\det \mathbf{\Phi} = \prod_{i=1}^{n} \lambda_i\) for eigenvalues \(\underline{\lambda}=(\lambda_1, \dots, \lambda_n)\) such that \(\mathbf{\Phi} x = \lambda x\) where \[|\Phi-\lambda I|=(\lambda _{1}-\lambda )(\lambda _{2}-\lambda )\cdots (\lambda _{n}-\lambda ).\]

• The VAR is stable then \(\mbox{det}(\Phi(Z)- \lambda \mathbf{I})=0\) has all roots within the complex unit circle. (E.g. Cointegration may make two non-stationary process stable.)

VAR stability

1. Let \(\mathbf{I} - \mathbf{\Phi}\mathbb{L} = \Phi(\mathbb{L})\), \(\det (\Phi (z)) \neq 0\) has all roots outside the complex unit circle

2. \(\mbox{det}(\mathbf{\Phi}- \lambda \mathbf{I}) = 0\) has all roots within the complex unit circle.

• 1 and 2 are equivalent.

plot(stability(BE.VAR)) 

Roots of AR characteristic polynomial (if \(|\lambda_1|, |\lambda_2|< 1\), then stable):

roots(BE.VAR)

## [1] 0.760952 0.345601

Null hypothesis: There is no serial correlated residul.

serial.test(BE.VAR)

## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object BE.VAR
## Chi-squared = 57.758, df = 60, p-value = 0.5581

• Use insights from AR(1) to study VAR(1). As in AR(1): \[ Y_{t+h} = \phi^{h} Y_t + \phi^{h-1}\varepsilon_{t+1} + \cdots + \phi \varepsilon_{t+h-1} +\varepsilon_{t+h}.\]

• Similarly in VAR(1), there is \[\mathbf{Y}_{t+h}=\mathbf{\Phi}^{h}\mathbf{Y}_{t}+\mathbf{\Phi}^{h-1}\mathbf{C}\mathbf{w}_{t+1}+\cdots+\mathbf{\Phi} \mathbf{C} \mathbf{w}_{t+h-1}+ \mathbf{C}\mathbf{w}_{t+h}\]

• Forecast: \(\mathbb{E}[\mathbf{Y}_{t+h} | \mathcal{F}_t] = \mathbf{\Phi}^h \mathbf{Y}_t\).

• Forecast error: \[\mathbf{Y}_{t+1} - \mathbb{E}[\mathbf{Y}_{t+1} | \mathcal{F}_t] = \mathbf{C} \mathbf{w}_{t+1}.\]

• Forecast error variance: \[ \mbox{Var}(\mathbf{Y}_{t+1} | \mathcal{F}_t) = \mathbf{C} \mathbf{C}^{T}.\]

• \(2\)-period ahead forecast error: \[\mathbf{Y}_{t+2} - \mathbb{E}[\mathbf{Y}_{t+2} | \mathcal{F}_t] = \mathbf{C} \mathbf{w}_{t+2} + \mathbf{\Phi}\mathbf{C} \mathbf{w}_{t+1}.\]

• \(2\)-period ahead forecast error variance: \[ \mbox{Var}(\mathbf{Y}_{t+2} | \mathcal{F}_t) = \mathbf{C} \mathbf{C}^{T} + \mathbf{\Phi} \mathbf{C} \mathbf{C}^{T} \mathbf{\Phi}^{T}.\]

• \(h\)-period ahead forecast error variance: \[\mbox{Var}(\mathbf{Y}_{t+h} | \mathcal{F}_t) = \sum_{j=0}^{h-1} \mathbf{\Phi}^{j} \mathbf{C} \mathbf{C}^{T} (\mathbf{\Phi}^{j})^{T}.\]

• A recursive form: \[\begin{align} \mathbb{E}[\mathbf{Y}_{t+h} | \mathcal{F}_t] = & \mathbf{\Phi}\mathbb{E}[\mathbf{Y}_{t+h-1} | \mathcal{F}_t], \\ \mbox{Var}(\mathbf{Y}_{t+h} | \mathcal{F}_t) =& \mathbf{C} \mathbf{C}^{T} + \mathbf{\Phi} \mbox{Var}(\mathbf{Y}_{t+h - 1} | \mathcal{F}_t) \mathbf{\Phi}^{T} . \end{align}\]

• The recursive form simplifies the program (we will see its usefulness in a recursive filter scheme).

BE.pred = predict(BE.VAR, n.ahead = 5)
plot(BE.pred)

• A single shock can affect the whole system.

• Special matrix structure allows us to orthogonalize the shocks and understand their impluse responses.

• Stationarity and Stability in the matrix structure.

• An empirical example.