This is page 383 [600212]

This is page 383
Printer: Opaque this
11
Vector Autoregressive Models for
Multivariate Time Series
11.1 Introduction
The vector autoregression (VAR) model is one of the most successful, flexi-
ble, and easy to use models for the analysis of multivariate time series. It is
a natural extension of the univariate autoregressive model to dynamic mul-
tivariate time series. The VAR model has proven to be especially useful fordescribing the dynamic behavior of economic and financial time series and
for forecasting. It often provides superior forecasts to those from univari-
ate time series models and elaborate theory-based simultaneous equations
models. Forecasts from VAR models are quite flexible because they can be
made conditional on the potential future paths of speci fied variables in the
model.
In addition to data description and forecasting, the VAR model is also
used for structural inference and policy a nalysis. In structural analysis, cer-
tain assumptions about the causal structure of the data under investiga-tion are imposed, and the resulting causal impacts of unexpected shocks or
innovations to speci fied variables on the variables in the model are summa-
rized. These causal impacts are usua lly summarized with impulse response
functions and forecast error variance decompositions.
This chapter focuses on the analysis of covariance stationary multivari-
ate time series using VAR models. The following chapter describes the
analysis of nonstationary multivariate time series using VAR models thatincorporate cointegration relationships.

384 11. Vector Autoregressive Models for Multivariate Time Series
This chapter is organized as follows. Section 11.2 describes speci fication,
estimation and inference in VAR models and introduces the S+FinMetrics
function VAR. Section 11.3 covers forecasting from VAR model. The discus-
sion covers traditional forecasting algorithms as well as simulation-basedforecasting algorithms that can impose certain types of conditioning infor-mation. Section 11.4 summarizes the types of structural analysis typically
performed using VAR models. These analyses include Granger-causality
tests, the computation of impulse resp onse functions, and forecast error
variance decompositions. Section 11.5 gives an extended example of VARmodeling. The chapter concludes wit h a brief discussion of Bayesian VAR
models.
This chapter provides a relatively non-technical survey of VAR models.
VAR models in economics were made popular by Sims (1980). The de finitive
technical reference for VAR models is L¨ utkepohl (1991), and updated sur-
veys of VAR techniques are given in Watson (1994) and L¨ utkepohl (1999)
and Waggoner and Zha (1999). Applications of VAR models to financial
data are given in Hamilton (1994), Campbell, Lo and MacKinlay (1997),
Cuthbertson (1996), Mills (1999) and Tsay (2001).
11.2 The Stationary Vector Autoregression Model
LetYt=(y1t,y2t,…,y nt)0denote an ( n×1) vector of time series variables.
The basic p-lagvector autoregressive (VAR(p)) model has the form
Yt=c+Π1Yt−1+Π2Yt−2+···+ΠpYt−p+εt,t=1,…,T (11.1)
whereΠiare (n×n)c o efficient matrices and εtis an (n×1) unobservable
zero mean white noise vector process (se rially uncorrelated or independent)
with time invariant covariance matrix Σ. For example, a bivariate VAR(2)
model equation by equation has the form
µy1t
y2t¶
=µc1
c2¶
+µπ1
11π112
π1
21π122¶µy1t−1
y2t−1¶
(11.2)
+µπ2
11π212
π221π222¶µy1t−2
y2t−2¶
+µε1t
ε2t¶
(11.3)
or
y1t=c1+π1
11y1t−1+π1
12y2t−1+π2
11y1t−2+π2
12y2t−2+ε1t
y2t=c2+π1
21y1t−1+π1
22y2t−1+π2
21y1t−1+π2
22y2t−1+ε2t
wherecov(ε1t,ε2s)=σ12fort=s; 0 otherwise. Notice that each equation
has the same regressors — lagged values of y1tandy2t. Hence, the VAR( p)
model is just a seemingly unrelated regression (SUR) model with lagged
variables and deterministic terms as common regressors.

11.2 The Stationary Vecto r Autoregression Model 385
In lag operator notation, the VAR( p) is written as
Π(L)Yt=c+εt
whereΠ(L)=In−Π1L−…−ΠpLp.T h eV A R ( p) is stable if the roots of
det (In−Π1z−···−Πpzp)=0
lie outside the complex unit circle (have modulus greater than one), or,
equivalently, if the eigenvalues of the companion matrix
F=
Π1Π2···Πn
In0 ··· 0
0…0…
00 I n0

have modulus less than one. Assuming that the process has been initialized
in the in finite past, then a stable VAR( p) process is stationary and ergodic
with time invariant means, variances, and autocovariances.
IfYtin (11.1) is covariance stationary, then the unconditional mean is
given by
µ=(In−Π1−···−Πp)−1c
The mean-adjusted form of the VAR( p)i st h e n
Yt−µ=Π1(Yt−1−µ)+Π2(Yt−2−µ)+···+Πp(Yt−p−µ)+εt
The basic VAR( p) model may be too restrictive to represent su fficiently
the main characteristics of the data. In particular, other deterministic terms
such as a linear time trend or seasonal dummy variables may be required
to represent the data properly. Additionally, stochastic exogenous variablesmay be required as well. The general form of the VAR( p)m o d e lw i t hd e –
terministic terms and exogenous variables is given by
Y
t=Π1Yt−1+Π2Yt−2+···+ΠpYt−p+ΦDt+GXt+εt (11.4)
where Dtrepresents an ( l×1) matrix of deterministic components, Xt
represents an ( m×1) matrix of exogenous variables, and ΦandGare
parameter matrices.
Example 64 Simulating a stationary VAR(1) model using S-PLUS
A stationary VAR model may be easily simulated in S-PLUS using the
S+FinMetrics function simulate.VAR . The commands to simulate T=
250 observations from a bivariate VAR(1) model
y1t=−0.7+0.7y1t−1+0.2y2t−1+ε1t
y2t=1.3+0.2y1t−1+0.7y2t−1+ε2t

386 11. Vector Autoregressive Models for Multivariate Time Series
with
Π1=µ0.70.2
0.20.7¶
,c=µ−0.7
1.3¶
,µ=µ1
5¶
,Σ=µ10.5
0.51¶
and normally distributed errors are
> pi1 = matrix(c(0.7,0.2,0.2,0.7),2,2)
> mu.vec = c(1,5)> c.vec = as.vector((diag(2)-pi1)%*%mu.vec)
> cov.mat = matrix(c(1,0.5,0.5,1),2,2)
> var1.mod = list(const=c.vec,ar=pi1,Sigma=cov.mat)> set.seed(301)
> y.var = simulate.VAR(var1.mod,n=250,
+ y0=t(as.matrix(mu.vec)))> dimnames(y.var) = list(NULL,c("y1","y2"))
The simulated data are shown in Figure 11.1. The VAR is stationary since
the eigenvalues of Π
1are less than one:
> eigen(pi1,only.values=T)
$values:[1] 0.9 0.5
$vectors:
NULL
Notice that the intercept values are quite di fferent from the mean values of
y
1andy2:
> c.vec
[1] -0.7 1.3> colMeans(y.var)
y1 y2
0.8037 4.751
11.2.1 Estimation
Consider the basic VAR( p) model (11.1). Assume that the VAR( p)m o d e l
is covariance stationary, and there are no restrictions on the parameters oft h em o d e l .I nS U Rn o t a t i o n ,e a c he q u a t i o ni nt h eV A R ( p)m a yb ew r i t t e n
as
y
i=Zπi+ei,i=1,…,n
where yiis a (T×1) vector of observations on the ithequation, Zis
a(T×k)m a t r i xw i t h tthrow given by Z0
t=( 1,Y0
t−1,…,Y0
t−p),k=
np+1,πiis a (k×1) vector of parameters and eiis a (T×1) error with
covariance matrix σ2
iIT.Since the VAR( p) is in the form of a SUR model

11.2 The Stationary Vecto r Autoregression Model 387
timey1,y2
0 50 100 150 200 250- 4 – 2 02468 1 0y1
y2
FIGURE 11.1. Simulated stationary VAR(1) model.
where each equation has the same explanatory variables, each equation may
be estimated separately by ordinary least squares without losing e fficiency
relative to generalized least squares. Let ˆΠ=[ˆπ1,…,ˆπn]d e n o t et h e( k×n)
matrix of least squares coe fficients for the nequations.
Letvec(ˆΠ) denote the operator that stacks the columns of the ( n×k)
matrix ˆΠi n t oal o n g( nk×1) vector. That is,
vec(ˆΠ)=
ˆπ1
…
ˆπn

Under standard assumptions regarding the behavior of stationary and er-
godic VAR models (see Hamilton (1994) or L¨ utkepohl (1991)) vec(ˆΠ)i s
consistent and asymptotically normally distributed with asymptotic covari-
ance matrix
[avar(vec(ˆΠ)) =ˆΣ⊗(Z0Z)−1
where
ˆΣ=1
T−kTX
t=1ˆεtˆε0
t
andˆεt=Yt−ˆΠ0Ztis the multivariate least squares residual from (11.1) at
timet.

388 11. Vector Autoregressive Models for Multivariate Time Series
11.2.2 Inference on Coe fficients
Theithelement of vec(ˆΠ), ˆπi, is asymptotically normally distributed with
standard error given by the square root of ithdiagonal element of ˆΣ⊗(Z0Z)−1.
Hence, asymptotically valid t-tests on individual coe fficients may be con-
structed in the usual way. More general linear hypotheses of the form
R·vec(Π)=rinvolving coe fficients across di fferent equations of the VAR
may be tested using the Wald statistic
Wald =(R·vec(ˆΠ)−r)0n
Rh
[avar(vec(ˆΠ))i
R0o−1
(R·vec(ˆΠ)−r) (11.5)
Under the null, (11.5) has a limiting χ2(q) distribution where q=rank (R)
gives the number of linear restrictions.
11.2.3 Lag Length Selection
The lag length for the VAR( p) model may be determined using model
selection criteria. The general approach is to fitV A R (p) models with orders
p=0,…,pmax and choose the value of pwhich minimizes some model
selection criteria. Model selection criteria for VAR( p)m o d e l sh a v et h ef o r m
IC(p)=l n |˜Σ(p)|+cT·ϕ(n,p)
where ˜Σ(p)=T−1PT
t=1ˆεtˆε0
tis the residual covariance matrix without a de-
grees of freedom correction from a VAR( p)m o d e l , cTis a sequence indexed
b yt h es a m p l es i z e T,a n dϕ(n,p) is a penalty function which penalizes
large VAR( p) models. The three most common information criteria are the
Akaike (AIC), Schwarz-Bayesian (BIC) and Hannan-Quinn (HQ):
AIC(p)=l n |˜Σ(p)|+2
Tpn2
BIC (p)=l n |˜Σ(p)|+lnT
Tpn2
HQ(p)=l n |˜Σ(p)|+2l nl nT
Tpn2
The AIC criterion asymptotically overestimates the order with positive
probability, whereas the BIC and HQ criteria estimate the order consis-tently under fairly general conditions if the true order pis less than or
equal to p
max. For more information on the use of model selection criteria
in VAR models see L¨ utkepohl (1991) chapter four.
11.2.4 Estimating VAR Models Using the S+FinMetrics
Function VAR
The S+FinMetrics function VARis designed to fit and analyze VAR models
as described in the previous section. VARproduces an object of class “ VAR”

11.2 The Stationary Vecto r Autoregression Model 389
for which there are print ,summary, plot and predict m e t h o d sa sw e l l
as extractor functions coefficients ,residuals ,fitted and vcov .T h e
calling syntax of VARis a bit complicated because it is designed to handle
multivariate data in matrices, data frames as well as “ timeSeries ”o b j e c t s .
The use of VARis illustrated with the following example.
Example 65 Bivariate VAR model for exchange rates
This example considers a bivariate VAR model for Yt=(∆st,fpt)0,
wherestis the logarithm of the monthly spot exchange rate between the US
and Canada, fpt=ft−st=iUS
t−iCA
tis the forward premium or interest
rate differential, and ftis the natural logarithm of the 30-day forward
exchange rate. The data over the 20 year period March 1976 through June
1996 is in the S+FinMetrics “timeSeries ”lexrates.dat . The data for
the VAR model are computed as
> dspot = diff(lexrates.dat[,"USCNS"])
> fp = lexrates.dat[,"USCNF"]-lexrates.dat[,"USCNS"]
> uscn.ts = seriesMerge(dspot,fp)> colIds(uscn.ts) = c("dspot","fp")> uscn.ts@title = "US/CN Exchange Rate Data"
> par(mfrow=c(2,1))
> plot(uscn.ts[,"dspot"],main="1st difference of US/CA spot+ exchange rate")> plot(uscn.ts[,"fp"],main="US/CN interest rate
+ differential")
Figure 11.2 illustrates the monthly return ∆s
tand the forward premium
fptover the period March 1976 through June 1996. Both series appear to be
I(0) (which can be con firmed using the S+FinMetrics functions unitroot
orstationaryTest )w i t h∆stmuch more volatile than fpt.fptalso ap-
pears to be heteroskedastic.
Specifying and Estimating the VAR(p) Model
To estimate a VAR(1) model for Ytuse
> var1.fit = VAR(cbind(dspot,fp)~ar(1),data=uscn.ts)Note that the VAR model is speci fied using an S-PLUS formula, with the
multivariate response on the left hand side of the ~operator and the built-
in AR term specifying the lag length of the model on the right hand side.The optional data argument accepts a data frame or “ timeSeries ”o b –
ject with variable names matching those used in specifying the formula.
If the data are in a “ timeSeries ” object or in an unattached data frame
(“timeSeries ” objects cannot be attached) then the data argument must
be used. If the data are in a matrix then the data argument may be omit-
ted. For example,

390 11. Vector Autoregressive Models for Multivariate Time Series
1st difference of US/CN spot exchange rate
1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996-0.06 -0.02 0.02
US/CN interest rate differential
1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996-0.005 -0.001 0.003
FIGURE 11.2. US/CN forward premium and spot rate.
> uscn.mat = as.matrix(seriesData(uscn.ts))
> var2.fit = VAR(uscn.mat~ar(1))
I ft h ed a t aa r ei na“ timeSeries ” object then the start andendoptions
may be used to specify the estimation sample. For example, to estimate theVAR(1) over the sub-period January 1980 through January 1990
> var3.fit = VAR(cbind(dspot,fp)~ar(1), data=uscn.ts,
+ start="Jan 1980", end="Jan 1990", in.format="%m %Y")
m a yb eu s e d .T h eu s eo f in.format=\%m %Y" sets the format for the date
strings speci fied in the start andendoptions to match the input format
of the dates in the positions slot of uscn.ts .
The VAR model may be estimated with the lag length pdetermined using
as p e c i fied information criterion. For example, to estimate the VAR for the
exchange rate data with pset by minimizing the BIC with a maximum lag
p
max=4u s e
> var4.fit = VAR(uscn.ts,max.ar=4, criterion="BIC")
> var4.fit$info
ar(1) ar(2) ar(3) ar(4)
BIC -4028 -4013 -3994 -3973
When a formula is not speci fied and only a data frame, “ timeSeries ”o r
matrix is supplied that contains the variables for the VAR model, VARfits

11.2 The Stationary Vecto r Autoregression Model 391
all VAR( p) models with lag lengths pless than or equal to the value given
tomax.ar , and the lag length is determined as the one which minimizes
the information criterion speci fied by the criterion option. The default
criterion is BICbut other valid choices are logL ,AICandHQ.I nt h ec o m –
putation of the information criteria, a common sample based on max.ar
is used. Once the lag length is determined, the VAR is re-estimated us-
ing the appropriate sample. In the above example, the BIC values were
computed using the sample based on max.ar=4 andp= 1 minimizes BIC.
The VAR(1) model was automatically re-estimated using the sample size
appropriate for p=1 .
Print and Summary Methods
The function VARproduces an object of class “ VAR” with the following
components.
> class(var1.fit)
[1] "VAR"
> names(var1.fit)
[1] "R" "coef" "fitted" "residuals"[5] "Sigma" "df.resid" "rank" "call"[9] "ar.order" "n.na" "terms" "Y0"
To see the estimated coe fficients of the model use the print method:
> var1.fit
Call:
VAR(formula = cbind(dspot, fp) ~ar(1), data = uscn.ts)
Coefficients:
dspot fp
(Intercept) -0.0036 -0.0003
dspot.lag1 -0.1254 0.0079
fp.lag1 -1.4833 0.7938
Std. Errors of Residuals:
dspot fp
0.0137 0.0009
Information Criteria:
logL AIC BIC HQ2058 -4104 -4083 -4096
total residual
Degree of freedom: 243 240Time period: from Apr 1976 to Jun 1996

392 11. Vector Autoregressive Models for Multivariate Time Series
The first column under the label “ Coefficients :” gives the estimated
coefficients for the ∆stequation, and the second column gives the estimated
coefficients for the fptequation:
∆st=−0.0036−0.1254 ·∆st−1−1.4833 ·fpt−1
fpt=−0.0003 + 0.0079 ·∆st−1+0.7938 ·fpt−1
Since uscn.ts is a “ timeSeries ” object, the estimation time period is also
displayed.
The summary method gives more detailed information about the fitted
VAR:
> summary(var1.fit)
Call:
VAR(formula = cbind(dspot, fp) ~ar(1), data = uscn.ts)
Coefficients:
dspot fp
(Intercept) -0.0036 -0.0003
(std.err) 0.0012 0.0001
(t.stat) -2.9234 -3.2885
dspot.lag1 -0.1254 0.0079
(std.err) 0.0637 0.0042
(t.stat) -1.9700 1.8867
fp.lag1 -1.4833 0.7938
(std.err) 0.5980 0.0395
(t.stat) -2.4805 20.1049
Regression Diagnostics:
dspot fp
R-squared 0.0365 0.6275
Adj. R-squared 0.0285 0.6244
Resid. Scale 0.0137 0.0009
Information Criteria:
logL AIC BIC HQ2058 -4104 -4083 -4096
total residual
Degree of freedom: 243 240Time period: from Apr 1976 to Jun 1996
In addition to the coe fficient standard errors and t-statistics, summary also
displaysR
2measures for each equation (which are valid because each equa-

11.2 The Stationary Vecto r Autoregression Model 393
tion is estimated by least squares). The summary output shows that the
coefficients on ∆st−1andfpt−1in both equations are statistically signi fi-
cant at the 10% level and that the fitf o rt h e fptequation is much better
than the fitf o rt h e ∆stequation.
As an aside, note that the S+FinMetrics function OLSmay also be used
to estimate each equation in a VAR model. For example, one way to com-
pute the equation for ∆stusing OLSis
> dspot.fit = OLS(dspot~ar(1)+tslag(fp),data=uscn.ts)
> dspot.fit
Call:
OLS(formula = dspot ~ar(1) + tslag(fp), data = uscn.ts)
Coefficients:
(Intercept) tslag(fp) lag1-0.0036 -1.4833 -0.1254
Degrees of freedom: 243 total; 240 residual
Time period: from Apr 1976 to Jun 1996Residual standard error: 0.01373
Graphical Diagnostics
Theplot method for “ VAR” objects may be used to graphically evaluate the
fitted VAR. By default, the plot method produces a menu of plot options:
> plot(var1.fit)
Make a plot selection (or 0 to exit):
1: plot: All
2: plot: Response and Fitted Values
3: plot: Residuals4: plot: Normal QQplot of Residuals5: plot: ACF of Residuals
6: plot: PACF of Residuals
7: plot: ACF of Squared Residuals8: plot: PACF of Squared ResidualsSelection:
Alternatively, plot.VAR may be called directly. The function plot.VAR
has arguments
> args(plot.VAR)
function(x, ask = T, which.plots = NULL, hgrid = F, vgrid= F, …)

394 11. Vector Autoregressive Models for Multivariate Time Series
-0.06 -0.04 -0.02 0.00 0.02
1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996dspot-0.004 0.000 0.004fpResponse and Fitted Values
FIGURE 11.3. Response and fitted values from VAR(1) model for US/CN ex-
change rate data.
To create all seven plots without using the menu, set ask=F .T oc r e a t e
the Residuals plot without using the menu, set which.plot=2 . The optional
arguments hgrid andvgrid control printing of horizontal and vertical grid
lines on the plots.
Figures 11.3 and 11.4 give the Respon se and Fitted Values and Residuals
plots for the VAR(1) fit to the exchange rate data. The equation for fptfits
much better than the equation for ∆st. The residuals for both equations
look fairly random, but the residuals for the fptequation appear to be
heteroskedastic. The qq-plot (not shown) indicates that the residuals for
the∆stequation are highly non-normal.
Extractor Functions
The residuals and fitted values for each equation of the VAR may be ex-
tracted using the generic extractor functions residuals andfitted :
> var1.resid = resid(var1.fit)
> var1.fitted = fitted(var.fit)> var1.resid[1:3,]
Positions dspot fp
Apr 1976 0.0044324 -0.00084150May 1976 0.0024350 -0.00026493Jun 1976 0.0004157 0.00002435

11.2 The Stationary Vecto r Autoregression Model 395-0.06 -0.04 -0.02 0.00 0.02
1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996dspot-0.004 -0.002 0.000 0.002fpResiduals versus Time
FIGURE 11.4. Residuals from VAR(1) model fit to US/CN exchange rate data.
Notice that since the data are in a “ timeSeries ” object, the extracted
residuals and fitted values are also “ timeSeries ” objects.
The coefficients of the VAR model may be extracted using the generic
coef function:
> coef(var1.fit)
dspot fp
(Intercept) -0.003595149 -0.0002670108
dspot.lag1 -0.125397056 0.0079292865
fp.lag1 -1.483324622 0.7937959055
Notice that coef produces the (3 ×2) matrix ˆΠwhose columns give the
estimated coe ffic i e n t sf o re a c he q u a t i o ni nt h eV A R ( 1 ) .
To test stability of the VAR, extract the matrix Π1and compute its
eigenvalues
> PI1 = t(coef(var1.fit)[2:3,])
> abs(eigen(PI1,only.values=T)$values)[1] 0.7808 0.1124
Since the modulus of the two eigenvalues of Π
1are less than 1, the VAR(1)
is stable.

396 11. Vector Autoregressive Models for Multivariate Time Series
Testing Linear Hypotheses
Now, consider testing the hypothesis that Π1=0(i.e.,Yt−1does not help
to explain Yt) using the Wald statistic (11.5). In terms of the columns of
vec(Π) the restrictions are π1=(c1,0,0)0andπ2=(c2,0,0) and may be
expressed as Rvec(Π)=rwith
R=
010000
001000
000010000001
,r=
0
0
00

The Wald statistic is easily constructed as follows
> R = matrix(c(0,1,0,0,0,0,
+ 0,0,1,0,0,0,+ 0,0,0,0,1,0,+ 0,0,0,0,0,1),
+ 4,6,byrow=T)
> vecPi = as.vector(var1.fit$coef)> avar = R%*%vcov(var1.fit)%*%t(R)> wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi)
> wald
[,1]
[1,] 417.1> 1-pchisq(wald,4)
[1] 0
Since the p-value for the Wald statistic based on the χ
2(4) distribution
is essentially zero, the hypothesis that Π1=0should be rejected at any
reasonable signi ficance level.
11.3 Forecasting
Forecasting is one of the main objectives of multivariate time series analysis.
Forecasting from a VAR model is similar to forecasting from a univariateAR model and the following gives a brief description.
11.3.1 Traditional Forecasting Algorithm
Consider first the problem of forecasting future values of Ytwhen the
parameters Πof the VAR( p) process are assumed to be known and there
are no deterministic terms or exogenous variables. The best linear predictor,in terms of minimum mean squared error (MSE), of Y
t+1or 1-step forecast
based on information available at time Tis
YT+1 |T=c+Π1YT+···+ΠpYT−p+1

11.3 Forecasting 397
Forecasts for longer horizons h(h-step forecasts) may be obtained using
thechain-rule of forecasting as
YT+h|T=c+Π1YT+h−1|T+···+ΠpYT+h−p|T
where YT+j|T=YT+jforj≤0. Theh-step forecast errors may be ex-
pressed as
YT+h−YT+h|T=h−1X
s=0ΨsεT+h−s
where the matrices Ψsare determined by recursive substitution
Ψs=p−1X
j=1Ψs−jΠj (11.6)
withΨ0=InandΠj=0f o rj>p .1The forecasts are unbiased since all of
the forecast errors have expectation zero and the MSE matrix for Yt+h|T
is
Σ(h)=MSE¡
YT+h−YT+h|T¢
=h−1X
s=0ΨsΣΨ0
s (11.7)
Now consider forecasting YT+hwhen the parameters of the VAR( p)
process are estimated using multivariate least squares. The best linear pre-
dictor of YT+his now
ˆYT+h|T=ˆΠ1ˆYT+h−1|T+···+ˆΠpˆYT+h−p|T (11.8)
where ˆΠjare the estimated parameter matrices. The h-step forecast error
is now
YT+h−ˆYT+h|T=h−1X
s=0ΨsεT+h−s+³
YT+h−ˆYT+h|T´
(11.9)
and the term³
YT+h−ˆYT+h|T´
captures the part of the forecast error due
to estimating the parameters of the VAR. The MSE matrix of the h-step
forecast is then
ˆΣ(h)=Σ(h)+MSE³
YT+h−ˆYT+h|T´
1The S+FinMetrics fucntion VAR.ar2ma computes the Ψsmatrices given the Πjma-
trices using (11.6).

398 11. Vector Autoregressive Models for Multivariate Time Series
In practice, the second term MSE³
YT+h−ˆYT+h|T´
is often ignored and
ˆΣ(h) is computed using (11.7) as
ˆΣ(h)=h−1X
s=0ˆΨsˆΣˆΨ0
s (11.10)
with ˆΨs=Ps
j=1ˆΨs−jˆΠj.L ¨utkepohl (1991, chapter 3) gives an approxi-
mation to MSE³
YT+h−ˆYT+h|T´
which may be interpreted as a finite
sample correction to (11.10).
Asymptotic (1 −α)·100% con fidence intervals for the individual elements
ofˆYT+h|Tare then computed as
£
ˆyk,T+h|T−c1−α/2ˆσk(h),ˆyk,T+h|T+c1−α/2ˆσk(h)¤
wherec1−α/2is the (1−α/2) quantile of the standard normal distribution
and ˆσk(h) denotes the square root of the diagonal element of ˆΣ(h).
Example 66 Forecasting exchange rates from a bivariate VAR
Consider computing h-step forecasts, h=1,…, 12, along with estimated
forecast standard errors from the bivariate VAR(1) model for exchangerates. Forecasts and forecast standard errors from the fitted VAR may be
computed using the generic S-PLUS predict method
> uscn.pred = predict(var1.fit,n.predict=12)
The predict function recognizes var1.fit as a “ VAR” object, and calls the
appropriate method function predict.VAR . Alternatively, predict.VAR
can be applied directly on an object inheriting from class “ VAR”. See the
online help for explanations of the arguments to predict.VAR .
The output of predict.VAR is an object of class “ forecast ”f o rw h i c h
there are print ,summary andplot methods. To see just the forecasts, the
print method will su ffice:
> uscn.pred
Predicted Values:
dspot fp
1-step-ahead -0.0027 -0.00052-step-ahead -0.0026 -0.00063-step-ahead -0.0023 -0.0008
4-step-ahead -0.0021 -0.0009
5-step-ahead -0.0020 -0.00106-step-ahead -0.0018 -0.00117-step-ahead -0.0017 -0.0011

11.3 Forecasting 399
8-step-ahead -0.0017 -0.0012
9-step-ahead -0.0016 -0.0012
10-step-ahead -0.0016 -0.0013
11-step-ahead -0.0015 -0.001312-step-ahead -0.0015 -0.0013
The forecasts and their standard errors can be shown using summary :
> summary(uscn.pred)
Predicted Values with Standard Errors:
dspot fp
1-step-ahead -0.0027 -0.0005
(std.err) 0.0137 0.0009
2-step-ahead -0.0026 -0.0006
(std.err) 0.0139 0.0012
…
12-step-ahead -0.0015 -0.0013
(std.err) 0.0140 0.0015
L¨utkepohl’s finite sample correction to the forecast standard errors com-
puted from asymptotic theory may be obtained by using the optional ar-
gument fs.correction=T in the call to predict.VAR .
The forecasts can also be plotted together with the original data using
the generic plot function as follows:
> plot(uscn.pred,uscn.ts,n.old=12)
where the n.old optional argument speci fie st h en u m b e ro fo b s e r v a t i o n st o
plot from uscn.ts .I fn.old is not speci fied, all the observations in uscn.ts
will be plotted together with uscn.pred . Figure 11.5 shows the forecasts
produced from the VAR(1) fit to the US/CN exchange rate data
2.A tt h e
beginning of the forecast horizon the spot return is below its estimated
mean value, and the forward premium is above its mean values. The spotreturn forecasts start o ffnegative and grow slowly toward the mean, and the
forward premium forecasts decline sharply toward the mean. The forecast
standard errors for both sets of forecasts, however, are fairly large.
2Notice that the dates associated with the forecasts are not shown. This is the result
of “timeDate ” objects not having a well de fined frequency from which to extrapolate
dates.

400 11. Vector Autoregressive Models for Multivariate Time Series
-0.01 0.0 0.01
235 240 245 250 255dspot
-0.002 -0.001 0.0235 240 245 250 255
fp
indexvalues
FIGURE 11.5. Predicted values from VAR(1) model fitt oU S / C Ne x c h a n g er a t e
data.
11.3.2 Simulation-Based Forecasting
T h ep r e v i o u ss u b s e c t i o ns h o w e dh o wt og e nerate multivariate forecasts from
afitted VAR model, using the chain-rule of forecasting (11.8). Since the
multivariate forecast errors (11.9) are asymptotically normally distributed
with covariance matrix (11.10), the forecasts of Yt+hcan be simulated
by generating multivariate normal random variables with mean zero andcovariance matrix (11.10). These simu lation-based forecasts can be ob-
tained by setting the optional argument method to"mc" in the call to
predict.VAR .
When method="mc" , the multivariate normal random variables are ac-
tually generated as a vector of standard normal random variables scaled
by the Cholesky factor of the covariance matrix (11.10). Instead of using
standard normal random variables, one could also use the standardizedresiduals from the fitted VAR model. Simulation-based forecasts based on
this approach are obtained by settin g the optional argument method to
"bootstrap" in the call to predict.VAR .
Example 67 Simulation-based forecasts of exchange rate data from bivari-
ate VAR
Theh-step forecasts ( h=1,…, 12) for∆s
t+handfpt+husing the Monte
Carlo simulation method are

11.3 Forecasting 401
> uscn.pred.MC = predict(var1.fit,n.predict=12,method="mc")
> summary(uscn.pred.MC)
Predicted Values with Standard Errors:
dspot fp
1-step-ahead -0.0032 -0.0005
(std.err) 0.0133 0.0009
2-step-ahead -0.0026 -0.0006
(std.err) 0.0133 0.0012
…
12-step-ahead -0.0013 -0.0013
(std.err) 0.0139 0.0015
The Monte Carlo forecasts and forecast standard errors for fpt+hare almost
identical to those computed using the chain-rule of forecasting. The MonteCarlo forecasts for ∆s
t+hare slightly di fferent and the forecast standard
errors are slightly larger than the corresponding values computed from the
chain-rule.
Theh-step forecasts computed from the bootstrap simulation method
are
> uscn.pred.boot = predict(var1.fit,n.predict=12,
+ method="bootstrap")> summary(uscn.pred.boot)
Predicted Values with Standard Errors:
dspot fp
1-step-ahead -0.0020 -0.0005
(std.err) 0.0138 0.0009
2-step-ahead -0.0023 -0.0007
(std.err) 0.0140 0.0012
…
12-step-ahead -0.0023 -0.0013
(std.err) 0.0145 0.0015
As with the Monte Carlo forecasts, the bootstrap forecasts and forecast
standard errors for fp
t+hare almost identical to those computed using the
chain-rule of forecasting. The bootstrap forecasts for ∆st+hare slightly
different from the chain-rule and Monte Carlo forecasts. In particular, the
bootstrap forecast standard errors are larger than corresponding values
from the chain-rule and Monte Carlo methods.
The simulation-based forecasts described above are di fferent from the
traditional simulation-based approach taken in VAR literature, e.g., see

402 11. Vector Autoregressive Models for Multivariate Time Series
Runkle (1987). The traditional approach is implemented using the following
procedure:
1. Obtain VAR coe fficient estimates Πand residuals εt.
2. Simulate the fitted VAR model by Monte Carlo simulation or by
bootstrapping the fitted residuals ˆεt.
3. Obtain new estimates of Πand forecasts of Yt+hbased on the sim-
ulated data.
The above procedure is repeated many times to obtain simulation-based
forecasts as well as their con fidence intervals. To illustrate this approach,
generate 12-step ahead forecasts from the fitted VAR object var1.fit by
Monte Carlo simulation using the S+FinMetrics function simulate.VAR
as follows:
> set.seed(10)
> n.pred=12
> n.sim=100
> sim.pred = array(0,c(n.sim, n.pred, 2))> y0 = seriesData(var1.fit$Y0)> for (i in 1:n.sim) {
+ dat = simulate.VAR(var1.fit,n=243)
+ dat = rbind(y0,dat)+ mod = VAR(dat~ar(1))+ sim.pred[i,,] = predict(mod,n.pred)$values
+}
The simulation-based forecasts are o btained by averaging the simulated
forecasts:
> colMeans(sim.pred)
[,1] [,2]
[1,] -0.0017917 -0.0012316
[2,] -0.0017546 -0.0012508[3,] -0.0017035 -0.0012643[4,] -0.0016800 -0.0012741
[5,] -0.0016587 -0.0012814
[6,] -0.0016441 -0.0012866[7,] -0.0016332 -0.0012904[8,] -0.0016253 -0.0012932
[9,] -0.0016195 -0.0012953
[10,] -0.0016153 -0.0012967[11,] -0.0016122 -0.0012978[12,] -0.0016099 -0.0012986

11.3 Forecasting 403
Comparing these forecasts with those in uscn.pred computed earlier, one
can see that for the first few forecasts, these simulated forecasts are slightly
different from the asymptotic forecasts. However, at larger steps, they ap-
proach the long run stable values of the asymptotic forecasts.
Conditional Forecasting
The forecasts algorithms considered up to now are unconditional multivari-
ate forecasts. However, sometimes it is desirable to obtain forecasts of some
variables in the system conditional on some knowledge of the future pathof other variables in the system. For example, when forecasting multivari-ate macroeconomic variables using quarterly data from a VAR model, it
may happen that some of the future values of certain variables in the VAR
model are known, because data on these variables are released earlier thandata on the other variables. By incorporating the knowledge of the futurepath of certain variables, in principle it should be possible to obtain more
reliable forecasts of the other variables in the system. Another use of con-
ditional forecasting is the generation of forecasts conditional on di fferent
“policy” scenarios. These scenario-based conditional forecasts allow one toanswer the question: if something hap pens to some variables in the system
in the future, how will it a ffect forecasts of other variables in the future?
S+FinMetrics provides a generic function cpredict for computing con-
ditional forecasts, which has a method cpredict.VAR for “ VAR”o b j e c t s .
The algorithms in cpredict.VAR are based on the conditional forecasting
algorithms described in Waggoner and Zha (1999). Waggoner and Zha clas-
sify conditional information into “hard” conditions and “soft conditions”.The hard conditions restrict the future values of certain variables at fixed
values, while the soft conditions restr ict the future values of certain vari-
ables in speci fied ranges. The arguments taken by cpredict.VAR are:
> args(cpredict.VAR)
function(object, n.predict = 1, newdata = NULL, olddata = NULL,
method = "mc", unbiased = T, variables.conditioned =NULL, steps.conditioned = NULL, upper = NULL, lower =NULL, middle = NULL, seed = 100, n.sim = 1000)
Like most predict methods in S-PLUS ,t h e first argument must be a fitted
model object, while the second argument, n.predict , speci fies the number
of steps to predict ahead. The arguments newdata andolddata can usually
be safely ignored, unless exogenous variables were used in fitting the model.
With classical forecasts that ignore the uncertainty in coe fficient esti-
mates, hard conditional forecasts can be obtained in closed form as shown
by Doan, Litterman and Sims (1984), and Waggoner and Zha (1999). Toobtain hard conditional forecasts, the argument middle is used to specify
fixed values of certain variables at certain steps. For example, to fixt h e

404 11. Vector Autoregressive Models for Multivariate Time Series
1-step ahead forecast of dspot invar1.fit at -0.005 and generate other
predictions for 2-step ahead forecasts, use the following command:
> cpredict(var1.fit, n.predict=2, middle=-0.005,
+ variables="dspot", steps=1)
Predicted Values:
dspot fp
1-step-ahead -0.0050 -0.00052-step-ahead -0.0023 -0.0007
In the call to cpredict , the optional argument variables is used to specify
the restricted variables, and steps to specify the restricted steps.
To specify a soft condition, the optional arguments upper and lower
are used to specify the upper bound and lower bound, respectively, of asoft condition. Since closed form results are not available for soft condi-
tional forecasts, either Monte Carlo simulation or bootstrap methods are
used to obtain the actual forecasts. T he simulations follow a similar proce-
dure implemented in the function predict.VAR , except that a reject/accept
method to sample from the distribution conditional on the soft conditions
is used. For example, to restrict the range of the first 2-step ahead forecasts
ofdspot to be (−0.004,−0.001) use:
> cpredict(var1.fit, n.predict=2, lower=c(-0.004, -0.004),
+ upper=c(-0.001, -0.001), variables="dspot",+ steps=c(1,2))
Predicted Values:
dspot fp
1-step-ahead -0.0027 -0.0003
2-step-ahead -0.0029 -0.0005
11.4 Structural Analysis
The general VAR( p) model has many parameters, and they may be di fficult
to interpret due to complex interactions and feedback between the variablesin the model. As a result, the dynamic properties of a VAR( p)a r eo f t e n
summarized using various types of structural analysis .T h et h r e em a i nt y p e s
of structural analysis summaries are (1) Granger causality tests ;( 2 ) impulse
response functions ;a n d( 3 ) forecast error variance decompositions .T h e
following sections give brief descr iptions of these summary measures.

11.4 Structural Analysis 405
11.4.1 Granger Causality
One of the main uses of VAR models is forecasting. The structure of the
VAR model provides information about a variable’s or a group of variables’
forecasting ability for other variables. The following intuitive notion of a
variable’s forecasting ability is due to Granger (1969). If a variable, orgroup of variables, y
1is found to be helpful for predicting another variable,
or group of variables, y2theny1is said to Granger-cause y2;o t h e r w i s ei t
is said to fail to Granger-cause y2. Formally, y1fails to Granger-cause y2
if for all s>0 the MSE of a forecast of y2,t+sbased on ( y2,t,y2,t−1,…)i s
the same as the MSE of a forecast of y2,t+sbased on ( y2,t,y2,t−1,…)a n d
(y1,t,y1,t−1,…). Clearly, the notion of Granger causality does not imply
true causality. It only implies forecasting ability.
Bivariate VAR Models
In a bivariate VAR( p)m o d e lf o r Yt=(y1t,y2t)0,y2fails to Granger-cause
y1if all of the pVAR coe fficient matrices Π1,…,Πpare lower triangular.
That is, the VAR( p) model has the form
µy1t
y2t¶
=µc1
c2¶
+µπ1
11 0
π1
21π122¶µy1t−1
y2t−1¶
+···
+µπp
11 0
πp
21πp22¶µy1t−p
y2t−p¶
+µε1t
ε2t¶
so that all of the coe fficients on lagged values of y2are zero in the equation
fory1. Similarly, y1fails to Granger-cause y2if all of the coe fficients on
lagged values of y1are zero in the equation for y2.T h eplinear coe fficient
restrictions implied by Granger non-causality may be tested using the Waldstatistic (11.5). Notice that if y
2fails to Granger-cause y1andy1fails
to Granger-cause y2,t h e nt h eV A Rc o e fficient matrices Π1,…,Πpare
diagonal.
General VAR Models
Testing for Granger non-causality in general nvariable VAR( p)m o d e l s
follows the same logic used for bivariate models. For example, consider aVAR(p)m o d e lw i t h n=3a n d Y
t=(y1t,y2t,y3t)0.I nt h i sm o d e l , y2does
not Granger-cause y1if all of the coe fficients on lagged values of y2are zero
in the equation for y1. Similarly, y3does not Granger-cause y1if all of the
coefficients on lagged values of y3are zero in the equation for y1.T h e s e
simple linear restrictions may be tested using the Wald statistic (11.5). The
reader is encouraged to consult L¨ utkepohl (1991) or Hamilton (1994) for
more details and examples.
Example 68 Testing for Granger causality in bivariate VAR(2) model for
exchange rates

406 11. Vector Autoregressive Models for Multivariate Time Series
Consider testing for Granger causality in a bivariate VAR(2) model for
Yt=(∆st,fpt)0. Using the notation of (11.2), fptdoes not Granger cause
∆stifπ1
12=0a n d π2
12= 0. Similarly, ∆stdoes not Granger cause fptif
π1
21=0a n d π2
21= 0. These hypotheses are easily tested using the Wald
statistic (11.5). The restriction matrix Rfor the hypothesis that fptdoes
not Granger cause ∆stis
R=µ0010000000
0000100000¶
and the matrix for the hypothesis that ∆stdoes not Granger cause fptis
R=µ0000001000
0000000010¶
The S-PLUS commands to compute and evaluate these Granger causality
Wald statistics are
> var2.fit = VAR(cbind(dspot,fp)~ar(2),data=uscn.ts)
> # H0: fp does not Granger cause dspot
> R = matrix(c(0,0,1,0,0,0,0,0,0,0,
+ 0,0,0,0,1,0,0,0,0,0),+ 2,10,byrow=T)> vecPi = as.vector(coef(var2.fit))
> avar = R%*%vcov(var2.fit)%*%t(R)
> wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi)> wald
[,1]
[1,] 8.468844
> 1-pchisq(wald,2)[1] 0.01448818
> R = matrix(c(0,0,0,0,0,0,1,0,0,0,
+ 0,0,0,0,0,0,0,0,1,0),+ 2,10,byrow=T)> vecPi = as.vector(coef(var2.fit))
> avar = R%*%vcov(var2.fit)%*%t(R)
> wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi)> wald
[,1]
[1,] 6.157
> 1-pchisq(wald,2)[1] 0.04604
Thep-values for the Wald tests indicate a fairly strong rejection of the null
thatfp
tdoes not Granger cause ∆stbut only a weak rejection of the null
that∆stdoes not Granger cause fpt. Hence, lagged values of fptappear

11.4 Structural Analysis 407
to be useful for forecasting future values of ∆stand lagged values of ∆st
appear to be useful for forecasting future values of fpt.
11.4.2 Impulse Response Functions
Any covariance stationary VAR( p) process has a Wold representation of
the form
Yt=µ+εt+Ψ1εt−1+Ψ2εt−2+··· (11.11)
where the ( n×n) moving average matrices Ψsare determined recursively
using (11.6). It is tempting to interpret the ( i,j)-th element, ψs
ij,o ft h e
matrixΨsas the dynamic multiplier or impulse response
∂yi,t+s
∂εj,t=∂yi,t
∂εj,t−s=ψs
ij,i , j =1,…,n
However, this interpretation is only possible if var(εt)=Σis a diagonal
matrix so that the elements of εtare uncorrelated. One way to make the
errors uncorrelated is to follow Sims (1980) and estimate the triangular
structural VAR(p)m o d e l
y1t=c1+γ0
11Yt−1+···+γ0
1pYt−p+η1t (11.12)
y2t=c1+β21y1t+γ0
21Yt−1+···+γ0
2pYt−p+η2t
y3t=c1+β31y1t+β32y2t+γ0
31Yt−1+···+γ0
3pYt−p+η3t
…
ynt=c1+βn1y1t+···+βn,n−1yn−1,t+γ0
n1Yt−1+···+γ0
npYt−p+ηnt
In matrix form, the triangular structural VAR( p)m o d e li s
BYt=c+Γ1Yt−1+Γ2Yt−2+···+ΓpYt−p+ηt (11.13)
where
B=
10 ··· 0
−β21 10 0
…………
−βn1−βn2··· 1
(11.14)
is a lower triangular matrix with 10salong the diagonal. The algebra of
least squares will ensure that the est imated covariance matrix of the error
vectorηtis diagonal. The uncorrelated/orthogonal errors ηtare referred
to as structural errors.
The triangular structural model (11.12) imposes the recursive causal or-
dering
y1→y2→···→yn (11.15)

408 11. Vector Autoregressive Models for Multivariate Time Series
The ordering (11.15) means that the contemporaneous values of the vari-
ables to the left of the arrow →affect the contemporaneous values of the
variables to the right of the arrow but not vice-versa. These contempora-
neous effects are captured by the coe fficientsβijin (11.12). For example,
the ordering y1→y2→y3imposes the restrictions: y1taffectsy2tandy3t
buty2tandy3tdo not affecty1;y2taffectsy3tbuty3tdoes not a ffecty2t.
Similarly, the ordering y2→y3→y1imposes the restrictions: y2taffects
y3tandy1tbuty3tandy1tdo not affecty2;y3taffectsy1tbuty1tdoes not
affecty3t.F o raV A R ( p)w i t hnvariables there are n! possible recursive
causal orderings. Which ordering to use in practice depends on the context
and whether prior theory can be used to justify a particular ordering. Re-
sults from alternative orderings can always be compared to determine thesensitivity of results to the imposed ordering.
Once a recursive ordering has been established, the Wold representation
ofY
tbased on the orthogonal errors ηtis given by
Yt=µ+Θ0ηt+Θ1ηt−1+Θ2ηt−2+··· (11.16)
whereΘ0=B−1is a lower triangular matrix. The impulse responses to
the orthogonal shocks ηjtare
∂yi,t+s
∂ηj,t=∂yi,t
∂ηj,t−s=θs
ij,i , j =1,…,n ;s>0 (11.17)
whereθs
ijis the (i,j)t he l e m e n to f Θs. A plot of θs
ijagainstsis called the
orthogonal impulse response function (IRF) of yiwith respect to ηj.W i t h
nvariables there are n2possible impulse response functions.
In practice, the orthogonal IRF (11.17) based on the triangular VAR( p)
(11.12) may be computed directly from the parameters of the non triangularVAR(p) (11.1) as follows. First, decompose the residual covariance matrix
Σas
Σ=ADA
0
where Ais an invertible lower triangular matrix with 10salong the diagonal
andDis a diagonal matrix with positive diagonal elements. Next, de fine
the structural errors as
ηt=A−1εt
These structural errors are orthogonal by construction since var(ηt)=
A−1ΣA−10=A−1ADA0A−10=D. Finally, re-express the Wold represen-
tation (11.11) as
Yt=µ+AA−1εt+Ψ1AA−1εt−1+Ψ2AA−1εt−2+···
=µ+Θ0ηt+Θ1ηt−1+Θ2ηt−2+···
whereΘj=ΨjA. Notice that the structural Bmatrix in (11.13) is equal
toA−1.

11.4 Structural Analysis 409
Computing the Orthogonal Impulse Response Function Using the
S+FinMetrics Function impRes
The orthogonal impulse response functi on (11.17) from a triangular struc-
tural VAR model (11.13) may be computed using the S+FinMetrics func-
tion impRes . The function impRes has arguments
> args(impRes)
function(x, period = NULL, std.err = "none", plot = F,
unbiased = T, order = NULL, …)
where xis an object of class “ VAR”a n d period speci fies the number of
responses to compute. By default, no standard errors for the responsesare computed. To compute asymptotic standard errors for the responses,
specify std.err="asymptotic" .T oc r e a t eap a n e lp l o to fa l lt h er e s p o n s e
functions, specify plot=T . The default recursive causal ordering is based on
t h eo r d e r i n go ft h ev a r i a b l e si n Y
twhen the VAR model is fit. The optional
argument order m a yb eu s e dt os p e c i f yad i fferent recursive causal order-
ing for the computation of the impulse responses. The argument order ac-
cepts a character vector of variable names whose order de fines the recursive
causal ordering. The output of impRes is an object of class “ impDecomp ”f o r
which there are print ,summary andplot methods. The following example
illustrates the use of impRes .
Example 69 IRF from VAR(1) for exchange rates
Consider again the VAR(1) model for Yt=(∆st,fpt)0.F o rt h ei m p u l s e
response analysis, the initial orderin g of the variables imposes the assump-
tion that structural shocks to fpthave no contemporaneous e ffect on∆st
but structural shocks to ∆stdo have a contemporaneous e ffect onfpt.T o
compute the four impulse response functions
∂∆st+h
∂η1t,∂∆st+h
∂η2t,∂fpt+h
∂η1t,∂fpt+h
∂η2t
forh=1,…, 12 we use S+FinMetrics function impRes .T h e first twelve
impulse responses from the VAR(1) model for exchange rates are computedusing
> uscn.irf = impRes(var1.fit, period=12, std.err="asymptotic")
The print method shows the impulse response values without standard
errors:
> uscn.irf
Impulse Response Function:
(with responses in rows, and innovations in columns)

410 11. Vector Autoregressive Models for Multivariate Time Series
, , lag.0
dspot fp
dspot 0.0136 0.0000
fp 0.0000 0.0009
, , lag.1
dspot fp
dspot -0.0018 -0.0013
fp 0.0001 0.0007
, , lag.2
dspot fp
dspot 0.0000 -0.0009
fp 0.0001 0.0006
…
, , lag.11
dspot fp
dspot 0.0000 -0.0001
fp 0.0000 0.0001
The summary method will display the responses with standard errors and
t-statistics. The plot method will produce a four panel Trellis graphics
plot of the impulse responses
> plot(uscn.irf)
A plot of the impulse responses can also be created in the initial call to
impRes by using the optional argument plot=T .
Figure 11.6 shows the impulse respon se functions along with asymptotic
standard errors. The top row shows the responses of ∆stto the structural
shocks, and the bottom row shows the responses of fptto the structural
shocks. In response to the first structural shock, η1t,∆stinitially increases
but then drops quickly to zero after 2 months. Similarly, fptinitially in-
creases, reaches its peak response in 2 months and then gradually dropsoffto zero after about a year. In response to the second shock, η
2t,b y
assumption ∆sthas no initial response. At one month, a sharp drop occurs
in∆stfollowed by a gradual return to zero after about a year. In contrast,
fptinitially increases and then gradually drops to zero after about a year.
The orthogonal impulse responses in Figure 11.6 are based on the recur-
sive causal ordering ∆st→fpt.I tm u s ta l w a y sb ek e p ti nm i n dt h a tt h i s
ordering identi fies the orthogonal structural shocks η1tandη2t.I ft h eo r –
dering is reversed, then a di fferent set of structural shocks will be identi fied,
and these may give very di fferent impulse response functions. To compute

11.4 Structural Analysis 4110.0 0.00010 0.00020
02468 1 0Inno.: dspotResp.: fp
0.0 0.0004 0.0008Inno.: fpResp.: fp0.0 0.005 0.010Inno.: dspotResp.: dspot
-0.0015 -0.0005 0.002468 1 0
Inno.: fpResp.: dspot
StepsImpulse ResponseOrthogonal Impulse Response Function
FIGURE 11.6. Impulse response function from VAR(1) model fitt oU S / C Ne x –
change rate data with ∆stordered first.
the orthogonal impulse responses using the alternative ordering fpt→∆st
specify order=c("fp","dspot") in the call to impRes :
> uscn.irf2 = impRes(var1.fit,period=12,std.err="asymptotic",
+ order=c("fp","dspot"),plot=T)
These impulse responses are presented in Figure 11.7 and are almost iden-
tical to those computed using the ordering ∆st→fpt. The reason for this
response is that the reduced form VAR residuals ˆ ε1tand ˆε2tare almost
uncorrelated. To see this, the residua l correlation matrix may be computed
using
> sd.vals = sqrt(diag(var1.fit$Sigma))
> cor.mat = var1.fit$Sigma/outer(sd.vals,sd.vals)
> cor.mat
dspot fp
dspot 1.000000 0.033048
fp 0.033048 1.000000
Because of the near orthogonality in the reduced form VAR errors, the
error in the ∆stequation may be interpreted as an orthogonal shock to the
exchange rate and the error in the fptequation may be interpreted as an
orthogonal shock to the forward premium.

412 11. Vector Autoregressive Models for Multivariate Time Series
-0.0020 -0.0005 0.001002468 1 0Inno.: fpResp.: dspot
0.0 0.005 0.010Inno.: dspotResp.: dspot0.0 0.0004 0.0008Inno.: fpResp.: fp
0.0 0.00005 0.0001502468 1 0
Inno.: dspotResp.: fp
StepsImpulse ResponseOrthogonal Impulse Response Function
FIGURE 11.7. Impulse response function from VAR(1) model fitt oU S / C Ne x –
change rate with fptordered first.
11.4.3 Forecast Error Variance Decompositions
The forecast error variance decomposition (FEVD) answers the question:
what portion of the variance of the forecast error in predicting yi,T+his
due to the structural shock ηj? Using the orthogonal shocks ηttheh-step
ahead forecast error vector, with known VAR coe fficients, may be expressed
as
YT+h−YT+h|T=h−1X
s=0ΘsηT+h−s
For a particular variable yi,T+h, this forecast error has the form
yi,T+h−yi,T+h|T=h−1X
s=0θs
i1η1,T+h−s+···+h−1X
s=0θs
inηn,T+h−s
Since the structural errors are orthogonal, the variance of the h-step fore-
cast error is
var(yi,T+h−yi,T+h|T)=σ2
η1h−1X
s=0(θs
i1)2+···+σ2
ηnh−1X
s=0(θs
in)2

11.4 Structural Analysis 413
whereσ2
ηj=var(ηjt). The portion of var(yi,T+h−yi,T+h|T)d u et os h o c k
ηjis then
FEVD i,j(h)=σ2
ηjPh−1
s=0¡
θs
ij¢2
σ2η1Ph−1
s=0(θs
i1)2+···+σ2ηnPh−1
s=0(θs
in)2,i , j =1,…,n
(11.18)
In a VAR with nvariables there will be n2FEVD i,j(h) values. It must be
kept in mind that the FEVD in (11.18) depends on the recursive causal or-
dering used to identify the structural shocks ηtand is not unique. Di fferent
causal orderings will produce di fferent FEVD values.
Computing the FEVD Using the S+FinMetrics Function fevDec
Once a VAR model has been fit, the S+FinMetrics function fevDec may
be used to compute the orthogonal FEVD. The function fevDec has argu-
ments
> args(fevDec)
function(x, period = NULL, std.err = "none", plot = F,
unbiased = F, order = NULL, …)
where xis an object of class “ VAR”a n d period speci fies the number of
responses to compute. By default, no standard errors for the responses
are computed and no plot is created . To compute asymptotic standard
errors for the responses, specify std.err="asymptotic" and to plot the
decompositions, specify plot=T . The default recursive causal ordering is
based on the ordering of the variables in Ytwhen the VAR model is fit. The
optional argument order may be used to specify a di fferent recursive causal
ordering for the computation of the FEVD. The argument order accepts
a text string vector of variable names whose order de fines the recursive
causal ordering. The output of fevDec is an object of class “ impDecomp ”
for which there are print ,summary andplot methods. The use of fevDec
is illustrated with the following example.
Example 70 FEVD from VAR(1) for exchange rates
The orthogonal FEVD of the forecast errors from the VAR(1) model
fit to the US/CN exchange rate data using the recursive causal ordering
∆st→fptis computed using
> uscn.fevd = fevDec(var1.fit,period=12,
+ std.err="asymptotic")> uscn.fevd
Forecast Error Variance Decomposition:
(with responses in rows, and innovations in columns)

414 11. Vector Autoregressive Models for Multivariate Time Series
, , 1-step-ahead
dspot fp
dspot 1.0000 0.0000
fp 0.0011 0.9989
, , 2-step-ahead
dspot fp
dspot 0.9907 0.0093
fp 0.0136 0.9864
…
, , 12-step-ahead
dspot fp
dspot 0.9800 0.0200
fp 0.0184 0.9816
The summary method adds standard errors to the above output if they are
computed in the call to fevDec .T h e plot method produces a four panel
Trellis graphics plot of the decompositions:
> plot(uscn.fevd)The FEVDs in Figure 11.8 show that most of the variance of the forecast
errors for ∆s
t+sat all horizons sis due to the orthogonal ∆stinnovations.
Similarly, most of the variance of the forecast errors for fpt+sis due to the
orthogonal fptinnovations.
The FEVDs using the alternativ e recursive causal ordering fpt→∆st
are computed using
> uscn.fevd2 = fevDec(var1.fit,period=12,
+ std.err="asymptotic",order=c("fp","dspot"),plot=T)
and are illustrated in Figure 11.9. Since the residual covariance matrix is
almost diagonal (see analysis of IRF above), the FEVDs computed usingthe alternative ordering are almost identical to those computed with the
initial ordering.
11.5 An Extended Example
In this example the causal relations and dynamic interactions among monthly
real stock returns, real interest rates , real industrial production growth and
the in flation rate is investigated using a VAR model. The analysis is similar
to that of Lee (1992). The variables are in the S+FinMetrics “timeSeries ”
object varex.ts
> colIds(varex.ts)

11.5 An Extended Example 4150.0 0.01 0.02 0.03 0.04
2468 1 0 1 2Inno.: dspotResp.: fp
0.96 0.97 0.98 0.99 1.00Inno.: fpResp.: fp0.97 0.98 0.99 1.00Inno.: dspotResp.: dspot
0.0 0.01 0.02 0.032468 1 0 1 2
Inno.: fpResp.: dspot
Forecast StepsProportion of Variance ContributionForecast Error Variance Decomposition
FIGURE 11.8. Orthogonal FEVDs computed from VAR(1) model fitt oU S / C N
exchange rate data using the re cursive causal ordering with ∆stfirst.
[1] "MARKET.REAL" "RF.REAL" "INF" "IPG"
Details about the data are in the documentation slot of varex.ts
> varex.ts@documentation
To be comparable to the results in Lee (1992), the analysis is conducted
over the postwar period January 1947 through December 1987
> smpl = (positions(varex.ts) >= timeDate("1/1/1947") &
+ positions(varex.ts) < timeDate("1/1/1988"))
The data over this period is displayed in Figure 11.10. All variables appear
to beI(0), but the real T-bill rate and the in flation rate appear to be highly
persistent.
To begin the analysis, autocorrelations and cross correlations at leads
and lags are computed using
> varex.acf = acf(varex.ts[smpl,])
and are illustrated in Figure 11.11. The real return on the market shows a
significant positive first lag autocorrelation, and in flation appears to lead
the real market return with a negative sign. The real T-bill rate is highly
positively autocorrelated, and in flation appears to lead the real T-bill rate
strongly with a negative sign. In flation is also highly positively autocorre-
lated and, interestingly, the real T-bill rate appears to lead in flation with

416 11. Vector Autoregressive Models for Multivariate Time Series
0.0 0.01 0.02 0.03 0.04
2468 1 0 1 2Inno.: fpResp.: dspot
0.96 0.97 0.98 0.99 1.00Inno.: dspotResp.: dspot0.975 0.985 0.995Inno.: fpResp.: fp
0.0 0.010 0.0202468 1 0 1 2
Inno.: dspotResp.: fp
Forecast StepsProportion of Variance ContributionForecast Error Variance Decomposition
FIGURE 11.9. Orthogonal FEVDs from VAR(1) model fit to US/CN exchange
rate data using recursive causal ordering with fptfirst.
a positive sign. Finally, industrial pro duction growth is slightly positively
autocorrelated, and the real market return appears to lead industrial pro-
d u c t i o ng r o w t hw i t hap o s i t i v es i g n .
The VAR( p)m o d e li s fit with the lag length selected by minimizing the
AIC and a maximum lag length of 6 months:
> varAIC.fit = VAR(varex.ts,max.ar=6,criterion="AIC",
+ start="Jan 1947",end="Dec 1987",
+ in.format="%m %Y")
The lag length selected by minimizing AIC is p=2 :
> varAIC.fit$info
ar(1) ar(2) ar(3) ar(4) ar(5) ar(6)
AIC -14832 -14863 -14853 -14861 -14855 -14862> varAIC.fit$ar.order
[1] 2
The results of the VAR(2) fita r e
> summary(varAIC.out)
Call:
VAR(data = varex.ts, start = "Jan 1947", end = "Dec 1987",

11.5 An Extended Example 417
Real Return on Market
1930 1950 1970 1990-0.3 -0.1 0.1 0.3Real Interest Rate
1930 1950 1970 1990-0.015 -0.005 0.005
Industrial Production Growth
1930 1950 1970 1990-0.10 0.00 0.10Inflation
1930 1950 1970 1990-0.02 0.02 0.05
FIGURE 11.10. Monthly data on stock returns, interest rates, output growth and
inflation.
max.ar = 6, criterion = "AIC", in.format = "%m %Y")
Coefficients:
MARKET.REAL RF.REAL INF IPG
(Intercept) 0.0074 0.0002 0.0010 0.0019
(std.err) 0.0023 0.0001 0.0002 0.0007
(t.stat) 3.1490 4.6400 4.6669 2.5819
MARKET.REAL.lag1 0.2450 0.0001 0.0072 0.0280
(std.err) 0.0470 0.0011 0.0042 0.0146
(t.stat) 5.2082 0.0483 1.7092 1.9148
RF.REAL.lag1 0.8146 0.8790 0.5538 0.3772
(std.err) 2.0648 0.0470 0.1854 0.6419
(t.stat) 0.3945 18.6861 2.9867 0.5877
INF.lag1 -1.5020 -0.0710 0.4616 -0.0722
(std.err) 0.4932 0.0112 0.0443 0.1533
(t.stat) -3.0451 -6.3147 10.4227 -0.4710
MARKET.REAL RF.REAL INF IPG
IPG.lag1 -0.0003 0.0031 -0.0143 0.3454

418 11. Vector Autoregressive Models for Multivariate Time Series
MARKET.REAL ACF
0 5 10 15 200.0 0.2 0.4 0.6 0.8 1.0 MARKET.REAL and RF.REAL
0 5 10 15 20-0.05 0.0 0.05 MARKET.REAL and INF
0 5 10 15 20-0.20 -0.15 -0.10 -0.05 0.0 0.05 0.10 MARKET.REAL and IPG
0 5 10 15 20-0.15 -0.10 -0.05 0.0 0.05
RF.REAL and MARKET.REALACF
-20 -15 -10 -5 0-0.05 0.0 0.05 0.10 RF.REAL
0 5 10 15 200.0 0.2 0.4 0.6 0.8 1.0 RF.REAL and INF
0 5 10 15 20-0.3 -0.2 -0.1 0.0 0.1 RF.REAL and IPG
0 5 10 15 20-0.05 0.0 0.05 0.10
INF and MARKET.REALACF
-20 -15 -10 -5 0-0.15 -0.10 -0.05 0.0 0.05 INF and RF.REAL
-20 -15 -10 -5 0-0.15 -0.10 -0.05 0.0 0.05 0.10 INF
0 5 10 15 200.0 0.2 0.4 0.6 0.8 1.0 INF and IPG
0 5 10 15 20-0.05 0.0 0.05
IPG and MARKET.REAL
LagACF
-20 -15 -10 -5 0-0.10 -0.05 0.0 0.05 0.10 0.15 0.20 IPG and RF.REAL
Lag-20 -15 -10 -5 0-0.05 0.0 0.05 IPG and INF
Lag-20 -15 -10 -5 0-0.20 -0.15 -0.10 -0.05 0.0 0.05 IPG
Lag0 5 10 15 20-0.2 0.0 0.2 0.4 0.6 0.8 1.0
FIGURE 11.11. Autocorrelations and cross correlations at leads and lags of data
in VAR model.
(std.err) 0.1452 0.0033 0.0130 0.0452
(t.stat) -0.0018 0.9252 -1.0993 7.6501
MARKET.REAL.lag2 -0.0500 0.0022 -0.0066 0.0395
(std.err) 0.0466 0.0011 0.0042 0.0145
(t.stat) -1.0727 2.0592 -1.5816 2.7276
RF.REAL.lag2 -0.3481 0.0393 -0.5855 -0.3289
(std.err) 1.9845 0.0452 0.1782 0.6169
(t.stat) -0.1754 0.8699 -3.2859 -0.5331
INF.lag2 -0.0602 0.0079 0.2476 -0.0370
(std.err) 0.5305 0.0121 0.0476 0.1649
(t.stat) -0.1135 0.6517 5.1964 -0.2245
MARKET.REAL RF.REAL INF IPG
IPG.lag2 -0.1919 0.0028 0.0154 0.0941
(std.err) 0.1443 0.0033 0.0130 0.0449
(t.stat) -1.3297 0.8432 1.1868 2.0968
Regression Diagnostics:
MARKET.REAL RF.REAL INF IPG

11.5 An Extended Example 419
R-squared 0.1031 0.9299 0.4109 0.2037
Adj. R-squared 0.0882 0.9287 0.4011 0.1905
Resid. Scale 0.0334 0.0008 0.0030 0.0104
Information Criteria:
logL AIC BIC HQ
7503 -14935 -14784 -14875
total residual
Degree of freedom: 489 480
Time period: from Mar 1947 to Nov 1987
The signs of the statistically signi ficant coefficient estimates corroborate
the informal analysis of the multivariate autocorrelations and cross lag au-tocorrelations. In particular, the real market return is positively related toits own lag but negatively related to the first lag of in flation. The real T-bill
rate is positively related to its own lag, negatively related to the first lag of
inflation, and positively related to the first lag of the real market return. In-
dustrial production growth is positively related to its own lag and positivelyrelated to the first two lags of the real market return. Judging from the
coefficients it appears that in flation Granger causes the real market return
and the real T-bill rate, the real T-bill rate Granger causes in flation, and
the real market return Granger causes the real T-bill rate and industrialproduction growth. These observations are con firmed with formal tests for
Granger non-causality. For example, t he Wald statistic for testing the null
hypothesis that the real market return does not Granger-cause industrialproduction growth is
> bigR = matrix(0,2,36)
> bigR[1,29]=bigR[2,33]=1> vecPi = as.vector(coef(varAIC.fit))
> avar = bigR%*%vcov(varAIC.fit)%*%t(bigR)
> wald = t(bigR%*%vecPi)%*%solve(avar)%*%(bigR%*%vecPi)> as.numeric(wald)[1] 13.82
> 1-pchisq(wald,2)
[1] 0.0009969
The 24-period IRF using the recursive causal ordering MARKET.REAL →
RF.REAL →IPG→INFis computed using
> varAIC.irf = impRes(varAIC.fit,period=24,
+ order=c("MARKET.REAL","RF.REAL","IPG","INF"),+ std.err="asymptotic",plot=T)
and is illustrated in Figure 11.12. The responses of MARKET.REAL to unex-
pected orthogonal shocks to the other variables are given in the first row of

420 11. Vector Autoregressive Models for Multivariate Time Series
-0.0003 0.0 0.00020 5 10 15 20Inno.: MARKET.REALResp.: INF
0.0 0.0004Inno.: RF.REALResp.: INF
-0.0002 0.0 0.0002
0 5 10 15 20Inno.: IPGResp.: INF
0.0 0.0015 0.0030Inno.: INFResp.: INF0.0 0.0010 0.0025Inno.: MARKET.REALResp.: IPG
-0.0002 0.0004Inno.: RF.REALResp.: IPG
0.0 0.004 0.008Inno.: IPGResp.: IPG
-0.0010 -0.0004 0.0002Inno.: INFResp.: IPG-1.5*10^-4 0Inno.: MARKET.REALResp.: RF.REAL
0.0002 0.0006Inno.: RF.REALResp.: RF.REAL
-2*10^-5 6*10^-5Inno.: IPGResp.: RF.REAL
-0.0004 -0.0001Inno.: INFResp.: RF.REAL0.0 0.01 0.03Inno.: MARKET.REALResp.: MARKET.REAL
-0.0015 0.0 0.00150 5 10 15 20
Inno.: RF.REALResp.: MARKET.REAL
-0.003 -0.001 0.001Inno.: IPGResp.: MARKET.REAL
-0.006 -0.003 0.00 5 10 15 20
Inno.: INFResp.: MARKET.REAL
StepsImpulse ResponseOrthogonal Impulse Response Function
FIGURE 11.12. IRF using the recursive causal ordering MARKET.REAL →RF.REAL
→IPG→INF.
thefigure. Most notable is the strong negative response of MARKET.REAL to
an unexpected increase in in flation. Notice that it takes about ten months
for the e ffect of the shock to dissipate. The responses of RF.REAL to the
orthogonal shocks is given in the second row of the figure. RF.REAL also
reacts negatively to an in flation shock and the e ffe c to ft h es h o c ki sf e l tf o r
about two years. The responses of IPGto the orthogonal shocks is given
in the third row of the figure. Industrial production growth responds posi-
tively to an unexpected shock to MARKET.REAL and negatively to shocks to
RF.REAL andINF.T h e s ee ffects, however, are generally short term. Finally,
the fourth row gives the responses of INFto the orthogonal shocks. In fla-
tion responds positively to a shock to the real T-bill rate, but this e ffect is
short-lived.
The 24 month FEVD computed using t he recursive causal ordering as
speci fied by MARKET.REAL →RF.REAL →IPG→INF,
> varAIC.fevd = fevDec(varAIC.out,period=24,
> order=c("MARKET.REAL","RF.REAL","IPG","INF"),> std.err="asymptotic",plot=T)

11.5 An Extended Example 421-0.002 0.004 0.008
51 0 1 5 2 0Inno.: MARKET.REALResp.: INF
0.01 0.03 0.05Inno.: RF.REALResp.: INF
-0.002 0.004 0.010
51 0 1 5 2 0Inno.: IPGResp.: INF
0.94 0.98Inno.: INFResp.: INF0.0 0.04Inno.: MARKET.REALResp.: IPG
-0.002 0.004 0.008Inno.: RF.REALResp.: IPG
0.90 0.94 0.98Inno.: IPGResp.: IPG
0.0 0.010 0.025Inno.: INFResp.: IPG0.01 0.03Inno.: MARKET.REALResp.: RF.REAL
0.5 0.7 0.9Inno.: RF.REALResp.: RF.REAL
-0.005 0.005 0.015Inno.: IPGResp.: RF.REAL
0.0 0.2 0.4Inno.: INFResp.: RF.REAL0.94 0.98Inno.: MARKET.REALResp.: MARKET.REAL
-0.0005 0.001051 0 1 5 2 0
Inno.: RF.REALResp.: MARKET.REAL
-0.002 0.004 0.010Inno.: IPGResp.: MARKET.REAL
0.0 0.02 0.04 0.0651 0 1 5 2 0
Inno.: INFResp.: MARKET.REAL
Forecast StepsProportion of Variance ContributionForecast Error Variance Decomposition
FIGURE 11.13. FEVDs using the recursive causal ordering MARKET.REAL →
RF.REAL →IPG→INF.
is illustrated in Figure 11.13. The first row gives the variance decom-
positions for MARKET.REAL and shows that most of the variance of the
forecast errors is due to own shocks. The second row gives the decomposi-tions for RF.REAL . At short horizons, most of the variance is attributable
to own shocks but at long horizons in flation shocks account for almost
half the variance. The third row gives the variance decompositions for IPG.
Most of the variance is due to own shocks and a small fraction is due toMARKET.REAL shocks. Finally, the fourth row shows that the forecast error
variance of INFis due mostly to its own shocks.
The IRFs and FEVDs computed above depend on the imposed recursive
causal ordering. However, in this example, the ordering of the variables willhave little e ffect on the IRFs and FEVDs because the errors in the reduced
form VAR are nearly uncorrelated:
> sd.vals = sqrt(diag(varAIC.out$Sigma))
> cor.mat = varAIC.out$Sigma/outer(sd.vals,sd.vals)> cor.mat
MARKET.REAL RF.REAL INF IPG
MARKET.REAL 1.00000 -0.16855 -0.04518 0.03916
RF.REAL -0.16855 1.00000 0.13046 0.03318
INF -0.04518 0.13046 1.00000 0.04732IPG 0.03916 0.03318 0.04732 1.00000

422 11. Vector Autoregressive Models for Multivariate Time Series
11.6 Bayesian Vector Autoregression
VAR models with many variables and long lags contain many parameters.
Unrestricted estimation of these models reqires lots of data and often theestimated parameters are not very precise, the results are hard to inter-pret, and forecasts may appear more precise than they really are because
standard error bands do not account for parameter uncertainty. The esti-
mates and forecasts can be improved if one has prior information aboutthe structure of the model or the possible values of the parameters or func-tions of the parameters. In a classical framework, it is di fficult to incorpo-
rate non-sample information into the estimation. Nonsample information
is easy to incorporate in a Bayesian framework. A Bayesian frameworkalso naturally incorporates parameter uncertainty into common measures
of precision. This section brie fly describes the Bayesian VAR modeling
tools in S+FinMetrics and illustrates these tools with an example. Details
of underlying Bayesian methods are given in Sims and Zha (1998) and Zha(1998).
11.6.1 An Example of a Bayesian VAR Model
S+FinMetrics comes with a “ timeSeries ”o b j e c t policy.dat ,w h i c hc o n –
tains six U.S. macroeconomic variables:
> colIds(policy.dat)
[1] "CP" "M2" "FFR" "GDP" "CPI" "U"
which represent IMF’s index of world commodity prices, M2 money stock,
federal funds rate, real GDP, consumer price index for urban consumers,
and civilian unemployment rate. The data set contains monthly observa-
tions from January 1959 to March 1998. Tao Zha and his co-authors haveanalyzed this data set in a number of papers, for example see Zha (1998).To use the same time period as in Zh a (1998), create a subset of the data:
> zpolicy.dat = policy.dat[1:264,]
> zpolicy.mat = as.matrix(seriesData(zpolicy.dat))
which contains monthly observations from January 1959 to December
1980.
Estimating a Bayesian VAR Model
To estimate a Bayesian vector autoregression model, use the S+FinMetrics
function BVAR . For macroeconomic modeling, it is usually found that many
trending macroeconomic variables have a unit root, and in some cases,
they may also have a cointegrating rel ationship (as described in the next
chapter). To incorporate these types of prior beliefs into the model, usetheunit.root.dummy and coint.dummy optional arguments to the BVAR

11.6 Bayesian Vector Autoregression 423
function, which add some dummy observ ations to the beginning of the data
to reflect these beliefs:
> zpolicy.bar13 = BVAR(zpolicy.mat~ar(13), unit.root=T,
+ coint=T)> class(zpolicy.bar13)[1] "BVAR"
The returned object is of class “ BVAR ”, which inherits from “ VAR”, so many
method functions for “ VAR” objects work similarly for “ BVAR ”o b j e c t s ,s u c h
as the extractor functions, impulse r esponse functions, and forecast error
variance decomposition functions.
The Bayesian VAR models are controlled through a set of hyper param-
eters, which can be speci fied using the optional argument control ,w h i c h
is usually a list returned by the function BVAR.control . For example, the
tightness of the belief in the unit root prior and cointegration prior is spec-
ified by mu5andmu6, respectively. To see what default values are used for
these hyper parameters, use
> args(BVAR.control)
function(L0 = 0.9, L1 = 0.1, L2 = 1, L3 = 1, L4 = 0.05,
m u 5=5 ,m u 6=5 )
For the meanings of these hyper parameters, see the online help file for
BVAR.control .
Adding Exogenous Variables to the Model
Other exogenous variables can be added to the estimation formula, just
as for OLSand VARfunctions. The BVAR function and related functions
will automatically take that into consideration and return the coe fficient
estimates for those variables.
Unconditional Forecasts
To forecast from a fitted Bayesian VAR model, use the generic predict
function, which automatically calls the method function predict.BVAR for
an object inheriting from class “ BVAR ”. For example, to compute 12-step
ahead forecasts use
> zpolicy.bpred = predict(zpolicy.bar13,n.predict=12)
> class(zpolicy.bpred)[1] "forecast"> names(zpolicy.bpred)
[1] "values" "std.err" "draws"
> zpolicy.bpred
Predicted Values:

424 11. Vector Autoregressive Models for Multivariate Time Series
4.5 4.6 4.7 4.8CP
4.3 4.4 4.5 4.6245 250 255 260 265 270 275
CPI
0.10 0.12 0.14 0.16 0.18 0.20FFR8.42 8.43 8.44 8.45 8.46 8.47
245 250 255 260 265 270 275GDP
7.25 7.30 7.35 7.40 7.45M2
0.06 0.07 0.08 0.09 0.10
245 250 255 260 265 270 275U
indexvalues
FIGURE 11.14. Forecasts from Bayesian VAR model.
CP M2 FFR GDP CPI U
1-step-ahead 4.6354 7.3794 0.1964 8.4561 4.4714 0.0725
2-step-ahead 4.6257 7.3808 0.1930 8.4546 4.4842 0.0732
3-step-ahead 4.6247 7.3834 0.1823 8.4505 4.4960 0.07464-step-ahead 4.6310 7.3876 0.1670 8.4458 4.5065 0.07635-step-ahead 4.6409 7.3931 0.1515 8.4414 4.5160 0.0785
6-step-ahead 4.6503 7.3998 0.1384 8.4394 4.5244 0.0810
7-step-ahead 4.6561 7.4075 0.1309 8.4390 4.5321 0.08338-step-ahead 4.6552 7.4159 0.1307 8.4403 4.5397 0.08529-step-ahead 4.6496 7.4242 0.1362 8.4428 4.5475 0.0867
10-step-ahead 4.6415 7.4323 0.1451 8.4453 4.5561 0.0879
11-step-ahead 4.6321 7.4402 0.1546 8.4473 4.5655 0.088912-step-ahead 4.6232 7.4476 0.1618 8.4482 4.5753 0.0899
The forecasts can also be plotted along with the original data using
> plot(zpolicy.bpred, zpolicy.mat, n.old=20)
The resulting plot is shown in Figure 11.14. The Bayesian forecasts usually
have wider error bands than classical forecasts, because they take into ac-
count the uncertainty in the coe fficient estimates. To ignore the uncertainty
in coefficient estimates, one can call the classical VAR predict method
function, predict.VAR , directly instead of the generic predict function.

11.6 Bayesian Vector Autoregression 425
The forecasts from Bayesian VAR models are of class “ forecast ”, and
are computed using Monte Carlo inte gration. By default, 1000 simulation
draws are used. To change the number of simulation draws and random
seed, specify the n.sim and seed optional arguments, respectively. For
forecasts from Bayesian VAR models, there is one more component in thereturned object: draws , which contains all the simulated forecasts. This
can be used to assess other statisti cal properties of the forecasts.
11.6.2 Conditional Forecasts
As mentioned earlier, conditional forecasts from classical VAR models ig-
nore the uncertainty in estimated coe fficients. In contrast, conditional fore-
casts from Bayesian VAR models take into account the uncertainty asso-
ciated with estimated coe fficients. To perform conditional forecasts from a
fitted Bayesian VAR model, use the generic cpredict function. For exam-
ple, if it is known that FFRin January 1981 is between 0 .185 and 0 .195, one
can incorporate this (soft condition) information into the forecasts using:
> zpolicy.spred = cpredict(zpolicy.bar13, 12, steps=1,
+ variables="FFR", upper=0.195, lower=0.185)> zpolicy.spred
Predicted Values:
CP M2 FFR GDP CPI U
1-step-ahead 4.6374 7.3797 0.1910 8.4554 4.4714 0.0729
2-step-ahead 4.6286 7.3816 0.1855 8.4540 4.4840 0.07363-step-ahead 4.6279 7.3850 0.1743 8.4498 4.4954 0.07524-step-ahead 4.6349 7.3899 0.1587 8.4452 4.5057 0.0768
5-step-ahead 4.6447 7.3960 0.1443 8.4414 4.5149 0.0791
6-step-ahead 4.6525 7.4033 0.1324 8.4406 4.5231 0.08147-step-ahead 4.6549 7.4114 0.1270 8.4412 4.5307 0.08358-step-ahead 4.6523 7.4201 0.1283 8.4428 4.5383 0.0851
9-step-ahead 4.6453 7.4284 0.1349 8.4457 4.5461 0.0864
10-step-ahead 4.6389 7.4365 0.1432 8.4482 4.5547 0.087611-step-ahead 4.6317 7.4444 0.1516 8.4501 4.5641 0.088512-step-ahead 4.6264 7.4519 0.1572 8.4511 4.5741 0.0896
For conditional forecasts with soft conditions, a Monte Carlo integra-
tion with acceptance/rejection method is used. By default, 1000 simulationdraws are used. However, it may occur that only a small number of draws
satisfy the soft conditions if the intervals are very tight. To see how many
draws satis fied the soft conditions and thus are used for inference, simply
check the dimension of the draws component of the returned object (see
t h eo n – l i n eh e l p file for forecast.object for details):

426 11. Vector Autoregressive Models for Multivariate Time Series
> dim(zpolicy.spred$draws)
[1] 372 72
In this case, only 372 out of 1000 simulation draws satis fied the con-
ditions. To continue simulating from the posterior moment distribution,
use the same command as before, with seed set to the current value of
.Random.seed :
> zpolicy.spred2 = cpredict(zpolicy.bar13, 12, steps=1,
+ variables="FFR", upper=0.195, lower=0.185, seed=.Random.seed)> dim(zpolicy.spred2$draws)[1] 389 72
Note that the draws in zpolicy.spred2 can be combined with the draws
inzpolicy.spred to obtain an updated and more accurate estimate of
conditional forecasts.
To ignore the coe fficient uncertainty for the conditional forecasts, call
the classical method function cpredict.VAR directly on a fitted Bayesian
VAR object. The technique introduced above can also be used for classical
prediction with soft conditions.
11.7 References
[1]C a m p b e l l ,J .A .L oa n dC .M a c K i n l a y (1997). The Economet-
rics of Financial Markets . Princeton University Press, New Jersey.
[2]Culbertson, K . (1996). Quantitative Financial Economics: Stocks,
Bonds and Foreign Exchange . John Wiley and Sons, Chichester.
[3]Doan, T. A., Litterman, R. B., and Sims, C. A. (1984). “Fore-
casting and Conditional Projection Using Realistic Prior Distribu-tions”, Econometric Reviews , 3, 1-100.
[4]Granger, C.W.J. (1969). “Investigating Causal Relations by Econo-
metric Models and Cross Spectral Methods,” Econometrica , 37, 424-
438.
[5]Hamilton, J.D. (1994). Time Series Analysis . Princeton University
Press, Princeton.
[6]Lee, B.-S. (1992). “Causal Relations Among Stock Returns, Interest
Rates, Real Activity, and In flation,” Journal of Finance , 47, 1591-
1603.
[7]Lutkepohl, H. (1991). Introduction to Multiple Time Series Anal-
ysis. Springer-Verlag, Berlin.

11.7 References 427
[8]Lutkepohl, H. (1999). “Vector Autoregressi ons,” unpublished manuscript,
Institut f¨ ur Statistik und ¨Okonometrie, Humboldt-Universit¨ at zu Berlin.
[9]Mills, T.C. (1999). T h eE c o n o m e t r i cM o d e l i n go fF i n a n c i a lT i m e
Series, Second Edition . Cambridge University Press, Cambridge.
[10]Runkle, D. E. (1987). “Vector Autoregressions and Reality,” Jour-
nal of Business and Economic Statistics , 5 (4), 437-442.
[11]Sims, C.A. (1980). “Macroeconomics and Reality,” Econometrica ,
48, 1-48.
[12]Sims, C. A., and Zha, T. (1998). “Bayesian Methods for Dynamic
Multivariate Models”, International Economic Review , 39 (4), 949-
968
[13]Stock, J.H. and M.W. Watson (2001). “Vector Autoregressions,”
Journal of Economic Perspectives , 15, 101-115.
[14]Tsay, R. (2001). Analysis of Financial Time Series . John Wiley &
Sons. New York.
[15]Watson, M. (1994). “Vector Autoregressions and Cointegration,” in
Handbook of Econometrics, Volume IV . R.F. Engle and D. McFadden
(eds.). Elsevier Science Ltd., Amsterdam.
[16]Waggoner, D. F., and Zha, T. (1999). “Conditional Forecasts in
Dynamic Multivariate Models,” Review of Economics and Statistics ,
81 (4), 639-651.
[17]Zha, T. (1998). “Dynamic Multivariate Model for Use in Formulating
Policy”, Economic Review , Federal Reserve Bank of Atlanta, First
Quarter, 1998.