Vector Autoregressive modeling

Seasonally adjusted GDP - UK, US, Canada

1 Seasonally adjusted GDP - UK, US, Canada

Quartelry growth rates, in percentages, of real gross domestic product (GDP) of United Kingdom, Canada, and United States from the second quarter of 1980 to the second quarter of 2011. The data were seasonally adjusted and downloaded from the database of Federal Reserve Bank of St. Louis. The GDP were in millions of local currency, and the growth rate denotes the differenced series of log DGP.

Source: Example 2.3, page 51 of Tsay (2014) Multivariate Time Series Analysis with R and Financial Applications. John Wiley. Hoboken, NJ.

rm(list=ls())
data  = read.table("https://hedibert.org/wp-content/uploads/2025/12/Example-2.3-page51-Tsay2014.txt",header=TRUE)
data = data[,c(1,2,3,5,4)]
n     = nrow(data)
date  = data[,1]+data[,2]/12
gdps  = log(data[,3:5])
rates = gdps[2:n,] - gdps[1:(n-1),]
rates = 100*rates

1.1 Exploratory data analysis

par(mfrow=c(3,3),mar=c(4,4,2,2))
plot(date,gdps[,1],type="l",xlab="Quarter",ylab="Log real DGP",main="United Kingdom")
plot(date,gdps[,2],type="l",xlab="Quarter",ylab="Log real DGP",main="United States")
plot(date,gdps[,3],type="l",xlab="Quarter",ylab="Log real DGP",main="Canada")
plot(date[2:n],rates[,1],type="l",xlab="Quarter",ylab="Growth rate")
plot(date[2:n],rates[,2],type="l",xlab="Quarter",ylab="Growth rate")
plot(date[2:n],rates[,3],type="l",xlab="Quarter",ylab="Growth rate")
acf(rates[,1],main="")
acf(rates[,2],main="")
acf(rates[,3],main="")

1.2 Standardizing the rates

rates = scale(rates)
par(mfrow=c(1,1))
plot(date[2:n],rates[,1],type="l",xlab="",ylab="Growth rate (standardized)",ylim=range(rates),lwd=2)
lines(date[2:n],rates[,2],col=2,lwd=2)
lines(date[2:n],rates[,3],col=3,lwd=2)
abline(h=0,lty=2)
legend("topright",legend=c("UK","USA","Canada"),col=1:3,lwd=2,bty="n")

1.3 Sample correlation matrix

cor(rates)

##           uk        us        ca
## uk 1.0000000 0.5074002 0.4093792
## us 0.5074002 1.0000000 0.6452283
## ca 0.4093792 0.6452283 1.0000000

2 Fitting VAR(3) via OLS (step-by-step)

The VAR(3) Gaussian model is fit for \(y_t=(y_{t1},y_{t2},y_{t3})'\) as follows: \[ y_t = B_1 y_{t-1} + B_2 y_{t-2} + B_3 y_{t-3} + \varepsilon_t \qquad \varepsilon \ \ iid \ \ N(0,\Sigma), \] with \(B_1\), \(B_2\) and \(B_3\), the \((3 \times 3)\) autoregressive coefficient matrices and \(\Sigma\) the residual covariance matrix. It is ease to rewrite the above VAR(3) as a standard Gaussian multivariate regression model: \[ y_t = B x_t + \varepsilon_t \] where \(B = (B_1,B_2,B_3)\) is a \((3 \times 9)\) matrix, while \(x_t=(y_{t-1}',y_{t-2}',y_{t-3}')'\) is a \(9 \times 1\) vector. Therefore, \[ \widehat{B} = (X'X)^{-1}X'Y \] where \(Y=(y_4,y_5,\ldots,y_n)'\) is \((n-3) \times 3\) and \(X=(x_4,x_5,\ldots,x_n)'\) is \((n-3) \times 9\).

n     = n - 1
p     = 3
y     = rates[(p+1):n,]
y     = as.matrix(y)
lag1  = rates[p:(n-1),]
lag2  = rates[(p-1):(n-2),]
lag3  = rates[(p-2):(n-3),]
X     = cbind(lag1,lag2,lag3)
X     = as.matrix(X)
XtX   = t(X)%*%X
iXtX  = solve(XtX)
Xty   = t(X)%*%y
bhat  = iXtX%*%Xty
beta  = c(bhat[,1],bhat[,2],bhat[,3])
A     = y - X%*%bhat
Shat  = t(A)%*%A/(123-(3+1)*2-1)
Cov   = kronecker(Shat,iXtX)
se    = sqrt(diag(Cov))
para  = cbind(beta,se,round(2*(1-pnorm(abs(beta)/se)),3))

2.1 Estimates

# B1
round(t(bhat),4)[,1:3]

##        uk     us     ca
## uk 0.4308 0.0587 0.1274
## us 0.4380 0.2405 0.1747
## ca 0.2938 0.4887 0.3231

# B2
round(t(bhat),4)[,4:6]

##         uk      us      ca
## uk  0.0008 -0.0439  0.1482
## us -0.2100  0.1710 -0.0946
## ca -0.1665  0.0228 -0.1915

# B3
round(t(bhat),4)[,7:9]

##        uk      us      ca
## uk 0.1010  0.1866 -0.2989
## us 0.0184 -0.0306 -0.0968
## ca 0.0422 -0.0858  0.0840

Shat

##            uk        us         ca
## uk 0.52718715 0.1543285 0.07821377
## us 0.15432847 0.5851116 0.22866668
## ca 0.07821377 0.2286667 0.49518778

3 Fitting VAR(3) - MTS package

#install.packages("MTS")
library("MTS")

mle.fit = VAR(rates,3)

## Constant term: 
## Estimates:  0.02894381 -0.004408488 -0.02165809 
## Std.Error:  0.06639475 0.07000534 0.06437023 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]   [,2]  [,3]
## [1,] 0.430 0.0576 0.128
## [2,] 0.438 0.2407 0.175
## [3,] 0.294 0.4895 0.322
## standard error 
##        [,1]   [,2]  [,3]
## [1,] 0.0906 0.0988 0.107
## [2,] 0.0955 0.1042 0.113
## [3,] 0.0878 0.0958 0.104
## AR( 2 )-matrix 
##           [,1]    [,2]    [,3]
## [1,]  0.000512 -0.0454  0.1487
## [2,] -0.209965  0.1713 -0.0947
## [3,] -0.166295  0.0240 -0.1919
## standard error 
##        [,1]  [,2]  [,3]
## [1,] 0.1020 0.106 0.110
## [2,] 0.1075 0.112 0.116
## [3,] 0.0988 0.103 0.106
## AR( 3 )-matrix 
##        [,1]    [,2]    [,3]
## [1,] 0.1011  0.1867 -0.2982
## [2,] 0.0184 -0.0306 -0.0969
## [3,] 0.0421 -0.0859  0.0835
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0932 0.0992 0.0945
## [2,] 0.0983 0.1046 0.0996
## [3,] 0.0904 0.0961 0.0916
##   
## Residuals cov-mtx: 
##           [,1]      [,2]      [,3]
## [1,] 0.4917830 0.1443357 0.0737094
## [2,] 0.1443357 0.5467243 0.2135770
## [3,] 0.0737094 0.2135770 0.4622492
##   
## det(SSE) =  0.09379619 
## AIC =  -1.934631 
## BIC =  -1.323715 
## HQ  =  -1.686448

3.1 VAR order selection

order = VARorder(rates)

## selected order: aic =  2 
## selected order: bic =  1 
## selected order: hq =  1 
## Summary table:  
##        p     AIC     BIC      HQ     M(p) p-value
##  [1,]  0 -1.6743 -1.6743 -1.6743   0.0000  0.0000
##  [2,]  1 -2.6013 -2.3976 -2.5185 115.1329  0.0000
##  [3,]  2 -2.6825 -2.2752 -2.5171  23.5389  0.0051
##  [4,]  3 -2.6418 -2.0309 -2.3936  10.4864  0.3126
##  [5,]  4 -2.6154 -1.8008 -2.2845  11.5767  0.2382
##  [6,]  5 -2.5001 -1.4819 -2.0864   2.7406  0.9737
##  [7,]  6 -2.4294 -1.2075 -1.9330   6.7822  0.6598
##  [8,]  7 -2.3362 -0.9107 -1.7571   4.5469  0.8719
##  [9,]  8 -2.4752 -0.8461 -1.8134  24.4833  0.0036
## [10,]  9 -2.4079 -0.5751 -1.6633   6.4007  0.6992
## [11,] 10 -2.3176 -0.2812 -1.4903   4.3226  0.8889
## [12,] 11 -2.3219 -0.0818 -1.4119  11.4922  0.2435
## [13,] 12 -2.3365  0.1072 -1.3437  11.8168  0.2238
## [14,] 13 -2.3901  0.2572 -1.3146  14.1266  0.1179

order

## $aic
##  [1] -1.674261 -2.601264 -2.682517 -2.641831 -2.615361 -2.500058 -2.429379
##  [8] -2.336182 -2.475226 -2.407881 -2.317578 -2.321865 -2.336480 -2.390056
## 
## $aicor
## [1] 2
## 
## $bic
##  [1] -1.67426063 -2.39762551 -2.27523971 -2.03091557 -1.80080654 -1.48186544
##  [7] -1.20754762 -0.91071230 -0.84611752 -0.57513390 -0.28119237 -0.08184084
## [13]  0.10718259  0.25724589
## 
## $bicor
## [1] 1
## 
## $hq
##  [1] -1.674261 -2.518536 -2.517062 -2.393649 -2.284450 -2.086420 -1.933014
##  [8] -1.757089 -1.813405 -1.663333 -1.490302 -1.411862 -1.343749 -1.314597
## 
## $hqor
## [1] 1
## 
## $Mstat
##  [1] 115.132873  23.538916  10.486416  11.576661   2.740611   6.782171
##  [7]   4.546892  24.483290   6.400690   4.322613  11.492247  11.816829
## [13]  14.126634
## 
## $Mpv
##  [1] 0.000000000 0.005092992 0.312559422 0.238240339 0.973697721 0.659786679
##  [7] 0.871885558 0.003599212 0.699241723 0.888925635 0.243469753 0.223833719
## [13] 0.117891395

par(mfrow=c(1,1))
yrange = range(order$aic,order$bic,order$hq)
plot(0:13,order$aic,xlab="Order",ylab="Value",type="b",pch=16,ylim=yrange)
lines(0:13,order$bic,pch=16,col=2,type="b")
lines(0:13,order$hq,pch=16,col=4,type="b")
legend("topleft",legend=c("AIC","BIC","HQ"),col=c(1,2,4),lty=1,pch=16)

3.2 Forecast Error Variance Decomposition

H      = 20
fit1   = refVAR(mle.fit)

## Constant term: 
## Estimates:  0 0 0 
## Std.Error:  0 0 0 
## AR coefficient matrix 
## AR( 1 )-matrix 
##       [,1]  [,2]  [,3]
## [1,] 0.446 0.000 0.141
## [2,] 0.431 0.231 0.150
## [3,] 0.306 0.472 0.327
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0824 0.0000 0.0822
## [2,] 0.0933 0.0999 0.1076
## [3,] 0.0848 0.0915 0.0934
## AR( 2 )-matrix 
##        [,1]  [,2]   [,3]
## [1,]  0.000 0.000  0.137
## [2,] -0.196 0.139  0.000
## [3,] -0.138 0.000 -0.175
## standard error 
##        [,1]  [,2]   [,3]
## [1,] 0.0000 0.000 0.1023
## [2,] 0.0956 0.102 0.0000
## [3,] 0.0851 0.000 0.0801
## AR( 3 )-matrix 
##        [,1]  [,2]   [,3]
## [1,] 0.0856 0.193 -0.305
## [2,] 0.0000 0.000 -0.145
## [3,] 0.0000 0.000  0.000
## standard error 
##        [,1]   [,2]   [,3]
## [1,] 0.0805 0.0974 0.0911
## [2,] 0.0000 0.0000 0.0819
## [3,] 0.0000 0.0000 0.0000
##   
## Residuals cov-mtx: 
##            [,1]      [,2]       [,3]
## [1,] 0.49497222 0.1438478 0.07301169
## [2,] 0.14384784 0.5516213 0.21667644
## [3,] 0.07301169 0.2166764 0.46840981
##   
## det(SSE) =  0.09657337 
## AIC =  -2.065452 
## BIC =  -1.680802 
## HQ  =  -1.909189

FEVD   = FEVdec(fit1$Phi, fit1$Phi0,fit1$Sigma,lag=H-1)

## Order of the ARMA mdoel:  
## [1] 3 0
## Standard deviation of forecast error:  
##           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
## [1,] 0.7035426 0.7822320 0.8320078 0.8934633 0.9108739 0.9147865 0.9189064
## [2,] 0.7427121 0.8584609 0.9084963 0.9265413 0.9373641 0.9389351 0.9394480
## [3,] 0.6844047 0.9002062 0.9616725 0.9797238 0.9917759 0.9970366 0.9976170
##           [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
## [1,] 0.9207173 0.9209199 0.9210735 0.9211810 0.9211937 0.9211967 0.9212009
## [2,] 0.9398933 0.9400134 0.9400395 0.9400487 0.9400527 0.9400536 0.9400536
## [3,] 0.9979197 0.9981924 0.9982483 0.9982567 0.9982626 0.9982644 0.9982647
##          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]
## [1,] 0.9212015 0.9212016 0.9212018 0.9212018 0.9212018 0.9212019
## [2,] 0.9400537 0.9400537 0.9400538 0.9400538 0.9400538 0.9400538
## [3,] 0.9982648 0.9982648 0.9982648 0.9982649 0.9982649 0.9982649
## Forecast-Error-Variance Decomposition 
## Forecast horizon:  1 
##            [,1]      [,2]      [,3]
## [1,] 1.00000000 0.0000000 0.0000000
## [2,] 0.07578528 0.9242147 0.0000000
## [3,] 0.02299206 0.1599804 0.8170276
## Forecast horizon:  2 
##           [,1]        [,2]       [,3]
## [1,] 0.9850854 0.002442194 0.01247240
## [2,] 0.2388760 0.749430152 0.01169389
## [3,] 0.1607360 0.316564870 0.52269909
## Forecast horizon:  3 
##           [,1]       [,2]       [,3]
## [1,] 0.9341113 0.02127101 0.04461774
## [2,] 0.2479349 0.73192133 0.02014375
## [3,] 0.2202995 0.32081438 0.45888612
## Forecast horizon:  4 
##           [,1]       [,2]       [,3]
## [1,] 0.8865944 0.06485824 0.04854732
## [2,] 0.2618325 0.71879139 0.01937616
## [3,] 0.2289797 0.32620771 0.44481256
## Forecast horizon:  5 
##           [,1]       [,2]       [,3]
## [1,] 0.8834373 0.06480959 0.05175313
## [2,] 0.2652673 0.70889471 0.02583795
## [3,] 0.2352403 0.32923218 0.43552750
## Forecast horizon:  6 
##           [,1]       [,2]       [,3]
## [1,] 0.8819903 0.06666172 0.05134802
## [2,] 0.2660172 0.70652763 0.02745517
## [3,] 0.2377681 0.32696405 0.43526788
## Forecast horizon:  7 
##           [,1]       [,2]       [,3]
## [1,] 0.8793437 0.06782211 0.05283419
## [2,] 0.2659328 0.70583382 0.02823339
## [3,] 0.2379812 0.32658598 0.43543281
## Forecast horizon:  8 
##           [,1]       [,2]       [,3]
## [1,] 0.8778561 0.06765323 0.05449068
## [2,] 0.2658154 0.70519073 0.02899392
## [3,] 0.2379757 0.32640368 0.43562061
## Forecast horizon:  9 
##           [,1]       [,2]       [,3]
## [1,] 0.8777303 0.06762619 0.05464353
## [2,] 0.2657478 0.70508186 0.02917036
## [3,] 0.2379173 0.32623937 0.43584332
## Forecast horizon:  10 
##           [,1]       [,2]       [,3]
## [1,] 0.8775803 0.06762038 0.05479927
## [2,] 0.2657369 0.70507965 0.02918345
## [3,] 0.2378906 0.32625742 0.43585193
## Forecast horizon:  11 
##           [,1]       [,2]       [,3]
## [1,] 0.8774524 0.06760636 0.05494121
## [2,] 0.2657318 0.70507129 0.02919694
## [3,] 0.2378878 0.32626153 0.43585068
## Forecast horizon:  12 
##           [,1]       [,2]       [,3]
## [1,] 0.8774341 0.06761196 0.05495393
## [2,] 0.2657296 0.70507219 0.02919820
## [3,] 0.2378852 0.32626036 0.43585439
## Forecast horizon:  13 
##           [,1]       [,2]       [,3]
## [1,] 0.8774312 0.06761164 0.05495712
## [2,] 0.2657296 0.70507183 0.02919860
## [3,] 0.2378845 0.32626235 0.43585318
## Forecast horizon:  14 
##           [,1]       [,2]       [,3]
## [1,] 0.8774267 0.06761107 0.05496221
## [2,] 0.2657296 0.70507178 0.02919864
## [3,] 0.2378844 0.32626236 0.43585324
## Forecast horizon:  15 
##           [,1]      [,2]       [,3]
## [1,] 0.8774261 0.0676114 0.05496254
## [2,] 0.2657296 0.7050717 0.02919864
## [3,] 0.2378845 0.3262623 0.43585319
## Forecast horizon:  16 
##           [,1]       [,2]       [,3]
## [1,] 0.8774261 0.06761139 0.05496253
## [2,] 0.2657296 0.70507165 0.02919873
## [3,] 0.2378845 0.32626232 0.43585316
## Forecast horizon:  17 
##           [,1]       [,2]       [,3]
## [1,] 0.8774260 0.06761138 0.05496265
## [2,] 0.2657296 0.70507164 0.02919873
## [3,] 0.2378845 0.32626230 0.43585318
## Forecast horizon:  18 
##           [,1]       [,2]       [,3]
## [1,] 0.8774260 0.06761137 0.05496267
## [2,] 0.2657297 0.70507161 0.02919873
## [3,] 0.2378845 0.32626231 0.43585315
## Forecast horizon:  19 
##           [,1]       [,2]       [,3]
## [1,] 0.8774260 0.06761137 0.05496266
## [2,] 0.2657297 0.70507161 0.02919873
## [3,] 0.2378845 0.32626231 0.43585314
## Forecast horizon:  20 
##           [,1]       [,2]       [,3]
## [1,] 0.8774260 0.06761138 0.05496267
## [2,] 0.2657297 0.70507161 0.02919873
## [3,] 0.2378845 0.32626231 0.43585314

vd.uk  = matrix(FEVD$OmegaR[1,],H,3,byrow=TRUE)
vd.can = matrix(FEVD$OmegaR[2,],H,3,byrow=TRUE)
vd.us  = matrix(FEVD$OmegaR[3,],H,3,byrow=TRUE)

par(mfrow=c(1,3))
ts.plot(vd.uk,main="UK",col=1:3,type="b",lwd=2)
ts.plot(vd.can,main="CAN",col=1:3,type="b",lwd=2)
ts.plot(vd.us,main="US",col=1:3,type="b",lwd=2)
legend("topright",legend=c("UK","CAN","US"),lwd=2,col=1:3,bty="n")

4 Impulse response functions

Bhat = matrix(para[,1],3,3*p,byrow=TRUE)
mB1  = Bhat[,1:3]
mB2  = Bhat[,4:6]
mB3  = Bhat[,7:9]

H        = 20
Psi      = array(0,c(H,p,p))
Psi[1,,] = diag(1,p)
Psi[2,,] = mB1%*%Psi[1,,]
Psi[3,,] = mB1%*%Psi[2,,]+mB2%*%Psi[1,,]
Psi[4,,] = mB1%*%Psi[3,,]+mB2%*%Psi[2,,]+mB3%*%Psi[1,,]
for (s in 5:H)
  Psi[s,,] = mB1%*%Psi[s-1,,]+mB2%*%Psi[s-2,,]+mB3%*%Psi[s-3,,]
  
names = c("UK","USA","CAD")
par(mfrow=c(p,p),mar=c(2,2,2,2))
for (i in 1:p)
  for (j in 1:p){
    plot(Psi[2:H,i,j],type="b",lwd=2,ylim=range(Psi[2:H,,]))
    title(paste("Response of ",names[j]," to shock in ",names[i],sep=""))
    abline(h=0,lty=3)
  }

4.1 IRF with confidence bands

H    = 20
R    = 1000
Psi  = array(0,c(H,p,p))
Psis = array(0,c(R,H,p,p))
for (r in 1:R){
  draw = matrix(t(chol(Cov))%*%rnorm(9*p),3,3*p,byrow=TRUE)
  dB1  = mB1+draw[,1:3]
  dB2  = mB2+draw[,4:6]
  dB3  = mB3+draw[,7:9]
  Psi[1,,] = diag(1,p)
  Psi[2,,] = dB1%*%Psi[1,,]
  Psi[3,,] = dB1%*%Psi[2,,]+dB2%*%Psi[1,,]
  Psi[4,,] = dB1%*%Psi[3,,]+dB2%*%Psi[2,,]+dB3%*%Psi[1,,]
  for (s in 5:H)
    Psi[s,,] = dB1%*%Psi[s-1,,]+dB2%*%Psi[s-2,,]+dB3%*%Psi[s-3,,]
  Psis[r,,,] = Psi
}

quants = array(0,c(p,p,H,3))
for (i in 1:p)
  for (j in 1:p)
    quants[i,j,,] = t(apply(Psis[,,i,j],2,quantile,c(0.025,0.5,0.975)))
   
par(mfrow=c(p,p),mar=c(2,2,2,2))
for (i in 1:p)
  for (j in 1:p){
    boxplot(Psis[,2:H,i,j],outline=FALSE,ylim=range(Psis[,2:H,,]))
    abline(h=0,lty=2)
    title(paste("Response of ",names[j]," to shock in ",names[i],sep=""))
}

par(mfrow=c(p,p),mar=c(2,2,2,2))
for (i in 1:p)
  for (j in 1:p){
    ts.plot(quants[i,j,2:H,],lty=c(3,1,3),ylim=range(Psis[,2:H,,]))
    abline(h=0,lty=2)
    title(paste("Response of ",names[j]," to shock in ",names[i],sep=""))
}