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("https://www.math.ttu.edu/~atrindad/tsdata/FTSdatasets/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
q     = ncol(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="")


# Scatter plots and correlations
################################
par(mfrow=c(1,3))
plot(rates[,1],rates[,2],xlab="Growth rate - UK",ylab="Growth rate (Canada)")
abline(lm(rates[,2]~rates[,1])$coef,col=2)
title(paste("Correlation = ",round(cor(rates[,1],rates[,2]),3),sep=""))
plot(rates[,1],rates[,3],xlab="Growth rate - UK",ylab="Growth rate (US)")
abline(lm(rates[,3]~rates[,1])$coef,col=2)
title(paste("Correlation = ",round(cor(rates[,1],rates[,3]),3),sep=""))
plot(rates[,2],rates[,3],xlab="Growth rate - Canada",ylab="Growth rate (US)")
abline(lm(rates[,3]~rates[,2])$coef,col=2)
title(paste("Correlation = ",round(cor(rates[,2],rates[,3]),3),sep=""))


# Fitting VAR(2) via OLS (step-by-step)
#######################################
p     = 2
n     = nrow(rates)
y     = rates[(p+1):n,]
y     = as.matrix(y)
const = rep(1,n-p)
lag1  = rates[p:(n-1),]
lag2  = rates[(p-1):(n-2),]
X     = cbind(const,lag1,lag2)
X     = as.matrix(X)
XtX   = t(X)%*%X
iXtX  = solve(XtX)
Xty   = t(X)%*%y
bhat  = iXtX%*%Xty
B     = c(bhat[,1],bhat[,2],bhat[,3])
A     = y - X%*%bhat
Smle  = t(A)%*%A/(n-p)

Cov  = kronecker(Smle,iXtX)
seB  = sqrt(diag(Cov))
para = cbind(B,seB,round(2*(1-pnorm(abs(B)/seB)),3))
colnames(para) = c("beta","std error","p-value")

para

Smle

# Impulse response function
#######################
B1 = t(bhat[2:(q+1),])
B2 = t(bhat[(q+2):(p*q+1),])

h = 11
Psi = array(0,c(h,q,q))
Psi[1,,] = diag(1,q)
Psi[2,,] = B1
for (i in 3:h)
  Psi[i,,] = B1%*%Psi[i-1,,]+B2%*%Psi[i-2,,]
 
names = c("United Kingdom","Canada","United States")
par(mfrow=c(1,q))
for (j in 1:q){
    plot(0:(h-1),Psi[,1,j],ylim=range(Psi),type="b",col=1:3,main=paste("Shock in ",names[j],sep=""),
    ylab="Response",xlab="horizon")
    lines(0:(h-1),Psi[,2,j],col=2,type="b")
    lines(0:(h-1),Psi[,3,j],col=3,type="b")
    abline(h=0,lty=2)
}
legend("topright",legend=names,col=1:3,lty=1)    
    
    


# Fitting VAR(2) via Ruey Tsay's MTS package
############################################
#https://cran.r-project.org/web/packages/MTS/MTS.pdf
install.packages("MTS")
library(MTS)

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

# 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)

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

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



# Bayesian VAR 
##############
# Kuschnig, N., & Vashold, L. (2021). BVAR: Bayesian Vector Autoregressions 
#  with Hierarchical Prior Selection in R. Journal of Statistical Software, 100(14), 1-27. 
# https://www.jstatsoft.org/article/view/v100i14
# https://cran.r-project.org/web/packages/BVAR/BVAR.pdf

set.seed(42) 
install.packages("BVAR")
library("BVAR")
x <- fred_qd[1:243, c("GDPC1", "PCECC96", "GPDIC1",  "HOANBS", "GDPCTPI", "FEDFUNDS")]
x <- fred_transform(x, codes = c(4, 4, 4, 4, 4, 1))

par(mfrow=c(2,3))
for (i in 1:6)
  ts.plot(x[,i])
  
mn <- bv_minnesota(lambda = bv_lambda(mode = 0.2, sd = 0.4, min = 0.0001, max = 5), alpha = bv_alpha(mode = 2), var = 1e07)
soc <- bv_soc(mode = 1, sd = 1, min = 1e-04, max = 50)
sur <- bv_sur(mode = 1, sd = 1, min = 1e-04, max = 50)
priors <- bv_priors(hyper = "auto", mn = mn, soc = soc, sur = sur)


mh <- bv_metropolis(scale_hess = c(0.05, 0.0001, 0.0001), adjust_acc = TRUE, acc_lower = 0.25, acc_upper = 0.45)


run <- bvar(x, lags = 5, n_draw = 50000, n_burn = 25000, n_thin = 1, priors = priors, mh = mh, verbose = TRUE)



print(run)

plot(run)

par(mfrow=c(3,3))

plot(run, type = "dens", vars_response = "GDPC1", vars_impulse = "GDPC1-lag1")


fitted(run, type = "mean")


plot(residuals(run, type = "mean"), vars = c("GDPC1", "PCECC96"))



opt_irf <- bv_irf(horizon = 16, identification = TRUE)
irf(run) <- irf(run, opt_irf, conf_bands = c(0.05, 0.16)) 
plot(irf(run), area = TRUE, vars_impulse = c("GDPC1", "FEDFUNDS"), vars_response = c(1:2, 6))


predict(run) <- predict(run, horizon = 16, conf_bands = c(0.05, 0.16))
plot(predict(run), area = TRUE, t_back = 32, vars = c("GDPC1", "GDPCTPI", "FEDFUNDS"))