######################################################################
# Econometrics III - Time series
#
# Bayesian inference via Gibbs sampler for the Gaussian CAPM model
#
######################################################################
rm(list=ls())

bayes.capm = function(y,x,a0,A0,b0,B0,c0,d0,M0,M,thin){
  alpha = a0
  beta  = b0
  sig2  = d0/(c0+1)
  niter = M0+M*thin
  draws = matrix(0,niter,3)
  for (iter in 1:niter){
    # Sampling from p(alpha|beta,sig2,y)
    var = 1/(1/A0+n/sig2)
    mean = var*(a0/A0+sum(y-beta*x)/sig2)
    alpha = rnorm(1,mean,sqrt(var))
    # Sampling from p(beta|alpha,sig2,y)
    var = 1/(1/B0+sum(x^2)/sig2)
    mean = var*(b0/B0+sum((y-alpha)*x)/sig2)
    beta = rnorm(1,mean,sqrt(var))
    # Sampling from p(sig2|alpha,beta,y)
    z = y-alpha-beta*x
    c1 = c0 + n/2
    d1 = d0 + sum(z^2)/2
    sig2 = 1/rgamma(1,c1,d1)
    # Storing the Gibbs output for posterior inference
    draws[iter,] = c(alpha,beta,sqrt(sig2))
  }
  ind = seq(M0+1,niter,by=thin)
  draws = draws[ind,]
  return(list(alphas=draws[,1],betas=draws[,2],sigmas=draws[,3]))
}

data = read.table("http://faculty.chicagobooth.edu/ruey.tsay/teaching/fts3/m-fac9003.txt",header=TRUE)
companies = c("Alcoa","A.G.Edwards","Caterpillar","Ford Motor","FedEx","General Motors","Hewlett-Packard",
              "Kimberly-Clark","Mellon Financial","New York Times","Procter & Gamble","Chicago Tribune",
              "Texas Instrument","S&P 500 Index")
n = nrow(data)
p = ncol(data)

ind = seq(1,n,by=12)
year = 1990:2003
par(mfrow=c(3,5))
for (i in 1:p){
  plot(data[,i],main=companies[i],axes=FALSE,xlab="Month",ylab="",type="l")
  axis(2);box();axis(1,at=ind,lab=year)
}             

# Let us first see what happens to a single time series
# We pick Texas Instrument

i = 13 # Texas Instrument
y = data[,i]
x = data[,14]

par(mfrow=c(1,2))
plot(y,main="Excess returns",axes=FALSE,xlab="Month",ylab="",type="l")
axis(2);box();axis(1,at=ind,lab=year)
lines(x,col=2)
legend("topleft",legend=companies[c(i,14)],col=1:2,lty=1)

plot(x,y,xlab=companies[14],ylab=companies[i])
abline(0,1,lty=2)
abline(lm(y~x),col=2)


# Bayesian inference via Gibbs sampler
# Priors and hyperparameters
# alpha~N(a0,A0), beta~N(b0,B0) and sig2~IG(c0,d0)
a0=0.0
A0=1.0
b0=1.0
B0=1.0
c0=0.1
d0=0.1

# Gibbs sampler set up
M0=1000
M=10000
thin=1

run = bayes.capm(y,x,a0,A0,b0,B0,c0,d0,M0,M,thin)

# Some plots of the marginal posteriors
par(mfrow=c(2,3))
ts.plot(run$alphas,xlab="iterations",ylab="",main=expression(alpha))
ts.plot(run$betas,xlab="iterations",ylab="",main=expression(beta))
ts.plot(run$sigmas,xlab="iterations",ylab="",main=expression(sigma^2))
hist(run$alphas,xlab="",prob=TRUE,main="")
hist(run$betas,xlab="",prob=TRUE,main="")
hist(run$sigmas,xlab="",prob=TRUE,main="")

pairs(cbind(run$alphas,run$betas,run$sigmas),labels=c("alpha","beta","sigma"))


# Now let us take a look at all of them (indivually, of course!)
summary = matrix(0,13,4)
for (i in 1:13){
  fit = lm(data[,i]~x)
  summary[i,] = c(fit$coef,summary(fit)$sigma,summary(fit)$r.squared)
}
rownames(summary)=companies[1:13]
colnames(summary)=c("alpha","beta","sigma","R2")
ord = order(summary[,2],decreasing=TRUE)
round(summary[ord,],3)

alphas = matrix(0,13,M)
betas  = matrix(0,13,M)
sigmas = matrix(0,13,M)
for (i in 1:13){
  run = bayes.capm(data[,i],data[,14],a0,A0,b0,B0,c0,d0,M0,M,thin)
  alphas[i,] = run$alphas
  betas[i,]  = run$betas
  sigmas[i,] = run$sigmas
}

par(mfrow=c(1,3),mar=c(5,9,4,2))
boxplot.matrix(t(alphas[ord,]),outline=FALSE,xlab=expression(alpha),ylab="",axes=FALSE,horizontal=TRUE)
axis(1);axis(2,at=1:13,lab=companies[ord[1:13]],las=2)
abline(v=0,col=2,lwd=3)
for (i in 1:13)
  points(summary[ord[i],1],i,col=4,pch=16)

boxplot.matrix(t(betas[ord,]),outline=FALSE,xlab=expression(beta),ylab="",axes=FALSE,horizontal=TRUE)
axis(1);axis(2,at=1:13,lab=companies[ord[1:13]],las=2)
abline(v=1,col=2,lwd=3)
for (i in 1:13)
  points(summary[ord[i],2],i,col=4,pch=16)

boxplot.matrix(t(sigmas[ord,]),outline=FALSE,xlab=expression(sigma),ylab="",axes=FALSE,horizontal=TRUE)
axis(1);axis(2,at=1:13,lab=companies[ord[1:13]],las=2)
abline(v=1,col=2,lwd=3)
for (i in 1:13)
  points(summary[ord[i],3],i,col=4,pch=16)
legend("topright",legend=c("MLE","Bayes"),col=c(1,4),pch=16)

med.beta  = apply(betas,1,median)
med.sigma = apply(sigmas,1,median)

par(mfrow=c(1,1))
plot(med.beta,med.sigma,xlab=expression(beta),ylab=expression(sigma),col=0)
text(med.beta,med.sigma,names(data)[1:13],col=c(1:4,6))
abline(v=1,lty=2)