########################################################################
# 
# BIVARIATE VAR(1): A BAYESIAN APPROACH
#
########################################################################
rm(list=ls())
riw2 = function(v,V){
  iV = solve(V)
	U  = chol(iV)
	T  = matrix(c(sqrt(rgamma(1,v/2,1/2)),rnorm(1),0,sqrt(rgamma(1,(v-1)/2,1/2))),2,2)
	C = solve(t(T)%*%U)
	return(C%*%t(C))
}

# Simulating 1000 observations from a bivariate VAR(1) 
set.seed(13579)
n  = 360
q  = 2
p  = 1
B  = matrix(c(0.85,0.1,0.0,0.95),q,q,byrow=TRUE)
#S  = matrix(c(0.26,0.03,0.03,0.09),q,q,byrow=TRUE)
S  = matrix(c(1,0.2,0.2,1),q,q,byrow=TRUE)
e  = matrix(rnorm(q*2*n),2*n,q)%*%chol(S)
y  = matrix(0,2*n,q)
for (t in 2:(2*n)){
  y[t,1] = B[1,1]*y[t-1,1]+B[1,2]*y[t-1,2]+e[t,1]
  y[t,2] = B[2,1]*y[t-1,1]+B[2,2]*y[t-1,2]+e[t,2]
}
y = y[(n+1):(2*n),]

pdf(file="bvar1-example-data.pdf",width=15,height=10)
par(mfrow=c(1,2))
ts.plot(y[,1],ylim=range(y),xlab="Time",ylab="",main="Time Series")
lines(y[,2],col=2)
plot(y,xlab="Y1",ylab="Y2",pch=16)
dev.off()

# OLS
eq1   = lm(y[2:n,1]~y[1:(n-1),1]+y[1:(n-1),2]-1)
eq2   = lm(y[2:n,2]~y[1:(n-1),1]+y[1:(n-1),2]-1)
b.ols = rbind(eq1$coef,eq2$coef)
s.ols = var(cbind(eq1$res,eq2$res))

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

# Other summaries
v1  = v0 + (n-1)
X   = y[1:(n-1),]
XtX = t(X)%*%X
Xty = t(X)%*%y[2:n,]

# Initial value
b  = B

# Gibbs sampler
set.seed(132154)
niter   = 10000
draws.b = array(0,c(niter,q,q))
draws.s = array(0,c(niter,q,q))
for (iter in 1:niter){
  # Sampling the variance-covariance matrix
  e1 = y[2:n,1]-b[1,1]*y[1:(n-1),1]-b[1,2]*y[1:(n-1),2]
  e2 = y[2:n,2]-b[2,1]*y[1:(n-1),1]-b[2,2]*y[1:(n-1),2]
  e  = cbind(e1,e2)
  V1 = V0+t(e)%*%e
  s  = riw2(v1,V1)

  # Sampling the coefficients of the 1st equation
  V     = solve(XtX/s[1,1]+diag(1/B0[1,]))
  m     = V%*%(Xty[,1]/s[1,1]+b0[1,]/B0[1,])
  b[1,] = m + t(chol(V))%*%rnorm(2)

  # Sampling the coefficients of the 2nd equation
  V     = solve(XtX/s[2,2]+diag(1/B0[2,]))
  m     = V%*%(Xty[,2]/s[2,2]+b0[2,]/B0[2,])
  b[2,] = m + t(chol(V))%*%rnorm(2)
  
  # Storing the draws
  draws.b[iter,,] = b
  draws.s[iter,,] = s
}

# MCMC and posterior summaries
ind = seq(1,niter,by=10)
pdf(file="bvar1-example-coef.pdf",width=15,height=10)
par(mfrow=c(2,2*q))
for (i in 1:q)
  for (j in 1:q){
    plot(ind,draws.b[ind,i,j],xlab="Iterations",ylab="",main="",ylim=range(draws.b[,i,j],B[i,j]),type="l")
    title(paste("B",i,j,sep=""))
  }
for (i in 1:q)
  for (j in 1:q){
  	  xxx = seq(min(draws.b[,i,j]),max(draws.b[,i,j]),length=30)
    hist(draws.b[,i,j],main="",prob=TRUE,xlab="",breaks=xxx)
    abline(v=B[i,j],col=2,lwd=2)  
    abline(v=b.ols[i,j],col=3,lwd=2)
    if (i==j){
      title(paste("Pr(B",i,j,">1)=",round(mean(draws.b[,i,i]>1),5),sep=""))
    }
  }
dev.off()

pdf(file="bvar1-example-var.pdf",width=15,height=10)
par(mfrow=c(2,q*(q+1)/2))
for (i in 1:q)
  for (j in 1:i){
    plot(ind,draws.s[ind,i,j],xlab="Iterations",ylab="",main="",ylim=range(draws.s[,i,j],S[i,j]),type="l")
    title(paste("S",i,j,sep=""))
  }
for (i in 1:q)
  for (j in 1:i){
  	  xxx = seq(min(draws.s[,i,j],S[i,j]),max(draws.s[,i,j],S[i,j]),length=30)
    hist(draws.s[,i,j],main="",prob=TRUE,xlab="",breaks=xxx)
    abline(v=S[i,j],col=2,lwd=2)  
    abline(v=s.ols[i,j],col=3,lwd=2)
  }
dev.off()