########################################################################
# 
# AN INTRODUCTION TO R
# A BVAR(1) EXAMPLE BASED ON BRAZILIAN MACROECONOMIC DATA
#
########################################################################
rm(list=ls())
riw = function(v,V){
	q  = nrow(V)
  iV = solve(V)
	U  = chol(iV)
	T  = matrix(0,q,q)
	for (i in 2:q)
    for (j in 1:(i-1))
      T[i,j] = rnorm(1)
  for (i in 1:q)
    T[i,i] = sqrt(rgamma(1,(v-i+1)/2,1/2))
	C = solve(t(T)%*%U)
	return(C%*%t(C))
}

pdf(file="BVAR1-braziliandata.pdf",width=20,height=15)

##############################
# reading the external data
##############################
names = c("Inflation","Interest Rate","Exchange Rate","Output Gap")
data  = matrix(scan("braziliandata.txt"),byrow=TRUE,ncol=6)
y     = data[,3:6]
q     = ncol(y)
n     = nrow(y)
date  = data[,1]+data[,2]/4

#####################################
# ploting the (demeaned) time series
#####################################
par(mfrow=c(2,2))
for (i in 1:q)
  plot(date,y[,i],main=names[i],xlab="Quarter",ylab="",pch=16,type="b")

##################################################
# Univariate OLS autoregressions
##################################################
x = y[1:(n-1),]

# Matrix derivation
P = solve(t(x)%*%x)%*%t(x)
B = matrix(0,q,q)
for (i in 1:q)
  B[i,] = P%*%y[2:n,i]

# Using R function lm
B = matrix(0,q,q)
for (i in 1:4)
  B[i,] = lm(y[2:n,i]~x-1)$coef
Bhat = B

# Error terms
e = matrix(0,n-1,q)
for (i in 1:4)
  e[,i] = lm(y[2:n,i]~x-1)$res
Shat = t(e)%*%e/(n-1)


##################################################
# Bayesian VAR(1) - unknown B and knwon S
##################################################

# Prior hyperparameters
b0  = matrix(0,q,q)
B0  = matrix(1000,q,q)

# Other summaries
y1  = y[2:n,]
X   = y[1:(n-1),]
XtX = t(X)%*%X
Xty = t(X)%*%y1
iB0 = 1/B0
iB0b0 = iB0*b0

# SIMPLE MONTE CARLO SAMPLING
set.seed(132154)
M        = 5000
draws.b1 = array(0,c(M,q,q))
V        = array(0,c(q,q,q))
m        = matrix(0,q,q)
tcV      = array(0,c(q,q,q))
for (i in 1:q){
  V[i,,] = solve(XtX/Shat[i,i]+diag(iB0[i,]))
  m[i,]  = V[i,,]%*%(Xty[,i]/Shat[i,i]+iB0b0[i,])
  tcV[i,,]    = t(chol(V[i,,]))
}
for (iter in 1:M)
  for (i in 1:q)
    draws.b1[iter,i,] = m[i,] + tcV[i,,]%*%rnorm(q)

# Posterior distributions (marginals)
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
  	  breaks = seq(min(draws.b1[,i,j]),max(draws.b1[,i,j]),length=20)
    hist(draws.b1[,i,j],prob=TRUE,xlab="",main="",breaks=breaks,col=grey(0.9))
    title(paste("B",i,j,sep=""))
  }

##################################################
# Bayesian VAR(1) - unknown B and S
##################################################

# Prior hyperparameters
b0  = matrix(0,q,q)
B0  = matrix(1000,q,q)
V0  = diag(0.001,q)
v0  = 0.001

# Other summaries
y1  = y[2:n,]
X   = y[1:(n-1),]
XtX = t(X)%*%X
Xty = t(X)%*%y1
v1  = v0 + (n-1)
iB0 = 1/B0
iB0b0 = iB0*b0
# Initial value
b  = Bhat

# Gibbs sampler
set.seed(132154)
burnin  = 2000
niter   = burnin + M
draws.b = array(0,c(niter,q,q))
draws.s = array(0,c(niter,q,q))
e       = matrix(0,n-1,q)
for (iter in 1:niter){
  print(iter)
  # Sampling the variance-covariance matrix
  e = y1-X%*%t(b)
  V1 = V0+t(e)%*%e
  s  = riw(v1,V1)
  # Sampling the autoregressive coefficients
  for (i in 1:q){
  	  V     = solve(XtX/s[i,i]+diag(iB0[i,]))
    m     = V%*%(Xty[,i]/s[i,i]+iB0b0[i,])
    b[i,] = m + t(chol(V))%*%rnorm(q)
  }
  # Storing the draws
  draws.b[iter,,] = b
  draws.s[iter,,] = s
}

# Trace plots
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
    ts.plot(draws.b[,i,j],xlab="Iterations",ylab="",main="")
    title(paste("B",i,j,sep=""))
  }
      
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
    ts.plot(draws.s[,i,j],xlab="Iterations",ylab="",main="")
    title(paste("S",i,j,sep=""))
  }
  
# Discarding the burnin draws
draws.b = draws.b[(burnin+1):niter,,]
draws.s = draws.s[(burnin+1):niter,,]

# Computing posterior means  
Bpost = matrix(0,q,q)
Spost = matrix(0,q,q)
for (i in 1:q)
  for (j in 1:q){
  	  Bpost[i,j] = mean(draws.b[,i,j])
  	  Spost[i,j] = mean(draws.s[,i,j])
  	}

# Posterior distributions (marginals)
par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
  	  breaks = seq(min(draws.b[,i,j]),max(draws.b[,i,j]),length=20)
    hist(draws.b[,i,j],prob=TRUE,xlab="",main="",breaks=breaks,col=grey(0.9))
    abline(v=Bhat[i,j],col=2,lwd=2)
    abline(v=Bpost[i,j],col=4,lwd=2)
    title(paste("B",i,j,sep=""))
  }

par(mfrow=c(q,q))
for (i in 1:q)
  for (j in 1:q){
  	  breaks = seq(min(draws.s[,i,j]),max(draws.s[,i,j]),length=20)
    hist(draws.s[,i,j],prob=TRUE,xlab="",main="",breaks=breaks,col=grey(0.9))
    abline(v=Shat[i,j],col=2,lwd=2)
    abline(v=Spost[i,j],col=4,lwd=2)
    title(paste("S",i,j,sep=""))
  }

par(mfrow=c(q,q))
for (i in 1:q){
  for (j in 1:q){
  	  if (i<=j){
  	    breaks = seq(min(draws.s[,i,j]),max(draws.s[,i,j]),length=20)
      hist(draws.s[,i,j],prob=TRUE,xlab="",main="",breaks=breaks,col=grey(0.9))
      title(paste("S",i,j,sep=""))
  	  	}else{
  	  	  corr = draws.s[,i,j]/sqrt(draws.s[,i,i]*draws.s[,j,j])
  	    breaks = seq(min(corr),max(corr),length=20)
      hist(corr,prob=TRUE,xlab="",main="",xlim=c(-1,1),breaks=breaks,col=grey(0.9))
      title(paste("R",i,j,sep=""))
      abline(v=0,lwd=2)
    }
  }
}

################################
#
# FORECASTING WITH A BVAR(1)
#
################################
h = 11

# Forecasting assuming that B=Bhat and S=Shat
yf0 = matrix(0,h+1,q)
yf0[1,] = y[n,]
for (i in 2:(h+1))
    yf0[i,] = Bhat%*%yf0[i-1,]
yf0 = yf0[2:(h+1),]

# Forecasting based on posterior of B (for S=Shat)
yf1 = array(0,c(M,h+1,q))
yf1[,1,] = matrix(y[n,],M,q,byrow=TRUE)
tcShat = t(chol(Shat))
for (i in 2:(h+1)){
  for (j in 1:M){
    yf1[j,i,] = draws.b1[j,,]%*%yf1[j,i-1,] + tcShat%*%rnorm(q)
  }
}
yf1 = yf1[,2:(h+1),]
quants1 = array(0,c(q,h,3))
for (i in 1:q)
  quants1[i,,]= t(apply(yf1[,,i],2,quantile,c(0.025,0.5,0.975)))
  
  
# Forecasting based on posterior of (B,S)
yf2      = array(0,c(M,h+1,q))
yf2[,1,] = matrix(y[n,],M,q,byrow=TRUE)
for (i in 2:(h+1)){
  for (j in 1:M){
    yf2[j,i,] = draws.b[j,,]%*%yf2[j,i-1,] + t(chol(draws.s[j,,]))%*%rnorm(q)
  }
}
yf2 = yf2[,2:(h+1),]
quants2 = array(0,c(q,h,3))
for (i in 1:q)
  quants2[i,,]= t(apply(yf2[,,i],2,quantile,c(0.025,0.5,0.975)))
  
  
ind = seq((date[n]+0.25),(date[n]+h/4),by=0.25)
par(mfrow=c(2,2))
for (i in 1:q){
  plot(date,y[,i],ylim=range(y[,i],quants1[i,,],quants2[i,,]),xlim=c(date[1],date[n]+h/4),
       main=names[i],xlab="Quarter",ylab="",pch=16,type="b")
  lines(ind,yf0[,i],type="b",pch=15)
  lines(ind,quants1[i,,2],col=2,type="b",pch=16)
  lines(ind,quants1[i,,1],type="b",pch=16,col=2)
  lines(ind,quants1[i,,3],type="b",pch=16,col=2)
  lines(ind,quants2[i,,2],col=4,type="b",pch=17)
  lines(ind,quants2[i,,1],type="b",pch=17,col=4)
  lines(ind,quants2[i,,3],type="b",pch=17,col=4)
  abline(v=data[n],lty=3)
  if (i==1){
  	  legend(2005,6,legend=c("(Bhat,Shat)","(B|Shat)","(B,S)"),col=c(1,2,4),lty=c(1,1,1),pch=15:17,bty="n",cex=2)
  	}
}
dev.off()