# VARexamples.ssc script file for examples used in VAR chapter # # author: Eric Zivot # created: December 4, 2001 # revised: March 16, 2005 # revised: July 6, 2007 # #################################################################### # # simulate VAR models # # simulate using loop nobs = 250 set.seed(301) e.var = rmvnorm(nobs,sd=c(1,1),rho=0.5) e.var = t(e.var) 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) y.var = matrix(0,2,nobs) y.var[,1] = mu.vec for (i in 2:nobs) { y.var[,i] = c.vec+pi1%*%y.var[,i-1]+e.var[,i] } y.var = t(y.var) dimnames(y.var) = list(NULL,c("y1","y2")) eigen(pi1,only.values=T) c.vec colMeans(y.var) tsplot(y.var,lty=1:2,ylab="y1, y2") legend(200,9,legend=colIds(y.var),lty=1:2) # simulate using simulate.VAR args(simulate.VAR) 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")) eigen(pi1,only.values=T) c.vec colMeans(y.var) tsplot(y.var,lty=1:2,ylab="y1, y2") legend(200,9,legend=colIds(y.var),lty=1:2) # # Example: bivariate VAR for spot and forward exchange rates # create data for VAR model of y(t) = (Ds(t), f(t)-s(t))' # motivation: can forward premium (f(t)-s(t)) be used to # predict spot exchange rates and vice-versa? # 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/CN spot exchange rate") plot(uscn.ts[,"fp"],main="US/CN interest rate differential") par(mfrow=c(1,1)) # # fit VAR with 1 lag - this lag is selected by BIC # var1.fit = VAR(cbind(dspot,fp)~ar(1),data=uscn.ts) args(VAR.formula) # # fit VAR with data in a matrix # uscn.mat = as.matrix(seriesData(uscn.ts)) var2.fit = VAR(uscn.mat~ar(1)) args(VAR.default) # # fit VAR over subsample January 1980 - January 1990 # var3.fit = VAR(cbind(dspot,fp)~ar(1),data=uscn.ts, start="Jan 1980",end="Jan 1990",in.format="%m %Y") var3.fit # # fit VAR using automatic lag selection # Note: default is BIC # Note: should the VAR object save the info values so that these may be plotted or # investigated etc? # var4.fit = VAR(uscn.ts,max.ar=4,criterion="BIC") var4.fit\$info # # show VAR object information # class(var1.fit) names(var1.fit) var1.fit summary(var1.fit) # # compute residual standard errors - make sure to normalize them since Sigma is unnormalized # sqrt(diag(var1.fit\$Sigma/var1.fit\$df.resid)) # # use OLS to compute VAR equations # dspot.fit = OLS(dspot~ar(1)+tslag(fp),data=uscn.ts) dspot.fit # # Illustrate plot method # args(plot.VAR) plot(var1.fit,ask=F) plot(var1.fit,which.plot=2) # extract residuals and fitted values var1.resid = resid(var1.fit) var1.fitted = fitted(var1.fit) var1.resid[1:3,] # extract coefficients coef(var1.fit) # test stability of VAR PI1 = t(coef(var1.fit)[2:3,]) abs(eigen(PI1,only.values=T)\$values) # test H0: PI1 = 0 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-pchisq(wald,4) # # 2. forecasting from bivariate VAR(1) for exchange rates # # compute standard error using asymptotic approximations uscn.pred = predict(var1.fit,n.predict=12) uscn.pred summary(uscn.pred) # plot forecasts with actual data plot(uscn.pred,uscn.ts,n.old=12) plot(uscn.pred,uscn.ts) # compute forecasts using Monte Carlo and bootstrap simulations uscn.pred.MC = predict(var1.fit,n.predict=12,method="mc") summary(uscn.pred.MC) uscn.pred.boot = predict(var1.fit,n.predict=12,method="bootstrap") summary(uscn.pred.boot) # illustrate traditional simulation-based forecasts 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 } # conditional predictions args(cpredict.VAR) # hard conditional forecasts cpredict(var1.fit, n.predict=2, middle=-0.005, variables="dspot", steps=1) # soft conditional forecasts cpredict(var1.fit, n.predict=2, lower=c(-0.004, -0.004), upper=c(-0.001, -0.001), variables="dspot", steps=c(1,2)) # # 3. Structural Analysis # # # 3.1 Granger causality tests in VAR(2) # var2.fit = VAR(cbind(dspot,fp)~ar(2),data=uscn.ts) var2.fit # H0: fp does not Granger cause dspot # R = matrix(c(0,0,1,0,0,0,0,0,0,0, # 0,0,0,0,0,1,0,0,0,0), # 2,10,byrow=T) # typo fixed on feb 22, 2005 # error reported by Frank Hansen 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-pchisq(wald,2) # H0: dspot does not Granger cause fp 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-pchisq(wald,2) # # 3.2 impulse response analysis from VAR(1) # # order dspot first # default ordering is based on ordering of variables in VAR fit args(impRes) uscn.irf = impRes(var1.fit,period=12) uscn.irf # compute asymptotic standard errors uscn.irf = impRes(var1.fit,period=12, std.err="asymptotic") # plot with standard errors plot(uscn.irf) uscn.irf = impRes(var1.fit,period=12, std.err="asymptotic",plot=T) # mention print and summary methods summary(uscn.irf) # order fp first: use order= option impRes(var1.fit,period=12,std.err="asymptotic", order=c("fp","dspot"),plot=T) # # check residual correlation - ordering shouldn't matter # Sigma = var1.fit\$Sigma sd.vals = sqrt(diag(Sigma)) cor.mat = Sigma/outer(sd.vals,sd.vals) cor.mat # # 3.3 compute forecast error variance decompositions # uscn.fevd = fevDec(var1.fit,period=12, std.err="asymptotic") uscn.fevd # plot with standard errors uscn.fevd = fevDec(var1.fit,period=12, std.err="asymptotic",plot=T) summary(uscn.fevd) # compute using different ordering uscn.fevd2 = fevDec(var1.fit,period=12, std.err="asymptotic",order=c("fp","dspot"),plot=T) # # 4. Bong-Soo Lee VAR for real stock returns, inflation, real output and real interest rates # # plot data smpl = (positions(varex.ts) >= timeDate("1/1/1947") & positions(varex.ts) < timeDate("1/1/1988")) par(mfrow=c(2,2)) plot(varex.ts[smpl,"MARKET.REAL"],main="Real Return on Market") plot(varex.ts[smpl,"RF.REAL"],main="Real Interest Rate") plot(varex.ts[smpl,"IPG"],main="Industrial Production Growth") plot(varex.ts[smpl,"INF"],main="Inflation") par(mfrow=c(1,1)) varex.acf = acf(varex.ts[smpl,]) par(mfrow=c(1,1)) # estimate VAR(6) as in Lee's paper var.out = VAR(cbind(MARKET.REAL,RF.REAL,IPG,INF)~AR(6), data=varex.ts,start="Jan 1947",end="Dec 1987", in.format="%m %Y") # let BIC choose the lag length # There appears to be an error in VAR.default. The below doesn't work with start # and end options specified. The following error is created: # Problem in VAR.formula(data = list(structure..: NAs are not allo # wed when na.rm=F. # Use traceback() to see the call stack varBIC.out = VAR(varex.ts,max.ar=6,criterion="BIC", start="Jan 1947",end="Dec 1987", in.format="%m %Y") # BIC picks one lag varAIC.out = VAR(varex.ts,max.ar=6,criterion="AIC", start="Jan 1947",end="Dec 1987", in.format="%m %Y") # AIC picks two lags varHQ.out = VAR(varex.ts,max.ar=6,criterion="HQ", start="Jan 1947",end="Dec 1987", in.format="%m %Y") varHQ.out # HQ picks two lags # plot all info criteria together tsplot(t(rbind(varBIC.out\$info,varAIC.out\$info,varHQ.out\$info)), main="Information Criteria") # look at some diagnostic plots: actual vs. fitted, ACF of residuals, Normal QQ plot(varAIC.out, ask=F) # # Pair-wise Granger non-causality tests # # 4 eqns, 2 lags, 4 variables and constant # mkt doesn't granger-cause T-bill bigR = matrix(0,2,36) bigR[1,11]=bigR[2,15]=1 vecPi = as.vector(coef(varAIC.out)) avar = bigR%*%vcov(varAIC.out)%*%t(bigR) wald = t(bigR%*%vecPi)%*%solve(avar)%*%(bigR%*%vecPi) wald 1-pchisq(wald,2) # mkt doesn't granger-cause inflation bigR = matrix(0,2,36) bigR[1,20]=bigR[2,24]=1 vecPi = as.vector(coef(varAIC.out)) avar = bigR%*%vcov(varAIC.out)%*%t(bigR) wald = t(bigR%*%vecPi)%*%solve(avar)%*%(bigR%*%vecPi) wald 1-pchisq(wald,2) # mkt doesn't granger-cause IPG bigR = matrix(0,2,36) bigR[1,29]=bigR[2,33]=1 vecPi = as.vector(coef(varAIC.out)) avar = bigR%*%vcov(varAIC.out)%*%t(bigR) wald = t(bigR%*%vecPi)%*%solve(avar)%*%(bigR%*%vecPi) as.numeric(wald) 1-pchisq(wald,2) # T-bill doesn't granger-cause mkt R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # T-bill doesn't granger cause inflation R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # T-bill doesn't granger-cause IPG R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # inflation doesn't granger-cause mkt R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # inflation doesn't granger-cause T-bill R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # inflation doesn't granger-cause IPG R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # IPG doesn't granger-cause mkt R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # IPG doesn't granger-cause T-bill R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # IPG doesn't granger-cause inflation R = matrix(c(0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,1,0,0,0,0), 2,10,byrow=T) vecPi = as.vector(coef(varAIC.out)) avar = R%*%vcov(varAIC.out)%*%t(R) wald = t(R%*%vecPi)%*%solve(avar)%*%(R%*%vecPi) wald 1-pchisq(wald,2) # compute impulse response functions and forecast error variance decompositions varAIC.irf = impRes(varAIC.out,period=24, order=c("MARKET.REAL","RF.REAL","IPG","INF"), std.err="asymptotic",plot=T) varAIC.fevd = fevDec(varAIC.out,period=24, order=c("MARKET.REAL","RF.REAL","IPG","INF"), std.err="asymptotic",plot=T) # print out 24 month FEVD varAIC.fevd\$values[,,24] varAIC.fevd\$std.err[,,24] summary(varAIC.fevd) # only care about 24th element in array # compute residual correlation matrix sd.vals = sqrt(diag(varAIC.out\$Sigma)) cor.mat = varAIC.out\$Sigma/outer(sd.vals,sd.vals) cor.mat # compute LM test for zero correlation values? # compute and plot 24 month forecasts var.pred = predict(varAIC.out,n.predict=24) plot(var.pred,varex.ts,n.old=12)