# VARexamples.ssc			script file for examples used in VAR chapter
#
# author: Eric Zivot
# created: December 4, 2001
# revised: May 11, 2002
#
####################################################################

#
# 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)
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(var.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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var2.fit))
avar = R%*%vcov(var2.fit)%*%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(var.out,n.predict=24)
plot(var.pred,varex.ts,n.old=12)