rm(list=ls())

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

data = read.table("http://faculty.chicagobooth.edu/ruey.tsay/teaching/bs41202/sp2013/q-gdp-ukcaus.txt",header=TRUE)
n = nrow(data)
date = data[,1]+data[,2]/12
gdps    = log(data[,3:5])
rates    = gdps[2:126,] - gdps[1:125,]
rates    = 100*rates

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

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...){
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(0, 1, 0, 1))
    r <- abs(cor(x, y))
    txt <- format(c(r, 0.123456789), digits = digits)[1]
    txt <- paste0(prefix, txt)
    if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
    text(0.5, 0.5, txt, cex = cex.cor * r)
}
pairs(rates,labels=c("United Kingdom","Canada","United States"),
         lower.panel=panel.smooth,upper.panel=panel.cor)

# Fitting VAR(2) via OLS (step-by-step)
################################
y         = rates[3:125,]
y         = as.matrix(y)
const = rep(1,123)
lag1   = rates[2:124,]
lag2  = rates[1:123,]
X      = cbind(const,lag1,lag2)
X      = as.matrix(X)
XtX  = t(X)%*%X
iXtX = solve(XtX)
Xty   = t(X)%*%y
bhat = iXtX%*%Xty
bhat
beta = c(bhat[,1],bhat[,2],bhat[,3])

A = y - X%*%bhat
Shat = t(A)%*%A/(123-(3+1)*2-1)
Shat

Cov = kronecker(Shat,iXtX)
se = sqrt(diag(Cov))
para = cbind(beta,se,round(2*(1-pnorm(abs(beta)/se)),3))
para

S.mle = t(A)%*%A/(123-2)


# Fitting VAR(2) via Ruey's MTS package
##################################
install.packages("MTS")
library(MTS)

# Maximum likelihood fit
mle.fit = VAR(rates,2)

# Bayesian fit via Minnesota prior
C = 0.1*diag(7)       ### lambda = 0.1
V0 = diag(3)           ### V0 = I_3
bayes.fit = BVAR(rates,p=2,C,V0)


# VAR order selection
##################
rates1 = rates/100
order = VARorder(rates1)

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


# Multivariate Ljung-Box Q Statistics
residuals = mle.fit$res

mq(residuals,adj=18)


# Forecasting
forec = VARpred(mle.fit,8)
colMeans(rates)
sqrt(apply(rates,2,var))

# Forecast Error Variance Decomposition
fit1 = refVAR(mle.fit)
FEVD = FEVdec(fit1$Phi, fit1$Phi0,fit1$Sigma, lag = 9)
vd.uk = matrix(FEVD$OmegaR[1,],10,3,byrow=TRUE)
vd.can = matrix(FEVD$OmegaR[2,],10,3,byrow=TRUE)
vd.us = matrix(FEVD$OmegaR[3,],10,3,byrow=TRUE)

par(mfrow=c(1,3))
ts.plot(vd.uk,main="UK",col=1:3)
ts.plot(vd.can,main="CAN",col=1:3)
ts.plot(vd.us,main="US",col=1:3)
legend("topright",legend=c("UK","CAN","US"),lty=1,col=1:3)

# Impulse response functions
irf.orig = VARirf(mle.fit$Phi,mle.fit$Sigma,lag=12,orth=F)                 # Original innovations

irf.orth = VARirf(mle.fit$Phi,mle.fit$Sigma,lag=12)                             # Orthogonal innovations

irf.uk  = matrix(irf.orig$irf[1,],13,3,byrow=TRUE)
irf.can = matrix(irf.orig$irf[2,],13,3,byrow=TRUE)
irf.us  = matrix(irf.orig$irf[3,],13,3,byrow=TRUE)
irf1.uk  = matrix(irf.orth$orthirf[1,],13,3,byrow=TRUE)
irf1.can = matrix(irf.orth$orthirf[2,],13,3,byrow=TRUE)
irf1.us  = matrix(irf.orth$orthirf[3,],13,3,byrow=TRUE)

par(mfrow=c(2,3))
ts.plot(irf.uk,main="UK",col=1:3)
ts.plot(irf.can,main="CAN",col=1:3)
ts.plot(irf.us,main="US",col=1:3)
legend("topright",legend=c("UK","CAN","US"),lty=1,col=1:3)
ts.plot(irf1.uk,main="UK",col=1:3)
ts.plot(irf1.can,main="CAN",col=1:3)
ts.plot(irf1.us,main="US",col=1:3)
legend("topright",legend=c("UK","CAN","US"),lty=1,col=1:3)