####################################################################################
#
#  Source: Description of Dataset found on the Wiley web site for 
#  Econometric Analysis of Panel Data, by Badi H. Baltagi. 
#  http://www.wiley.co.uk/wileychi/baltagi/datasets.html
#
#  Source: Grunfeld (1958)
#  Bayesian analysis: Zellner (1971), Section 8.5, Table 8.1, page 244-245.
#
#
#  Description:  Panel Data, 10 U.S. firms over 20 years, 1935 1954.
#  Variables:  (1) FN = Firm Number.
#              (2) YR = Year.
#              (3) I = Annual real gross investment.
#              (4) F = Real value of the firm (shares outstanding).
#              (5) C = Real value of the capital stock.
#
#  Firms:  General Electric
#          Westinghouse
#          U.S. Steel
#          Diamond Match
#          Atlantic Refining
#          Union Oil
#          Goodyear
#          General Motors
#          Chrysler
#          IBM
#
####################################################################################
# Author : Hedibert Freitas Lopes
#          Booth School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoBooth.edu
#
####################################################################################
rm(list=ls())
M       = 10
T       = 20
K       = 3
years   = 1935:1954
data    = matrix(scan("Grunfeld.txt"),T*M,5,byrow=T)
names   = c("Annual real gross investment","Real value of the firm (shares outstanding)","Real value of the capital stock")
companies = c("General Electric","Westinghouse","U.S. Steel","Diamond Match","Atlantic Refining","Union Oil",
"Goodyear","General Motors","Chrysler","IBM")
par(mfrow=c(1,3))
for (k in 1:K){
  plot(years,data[data[,1]==1,2+k],type="l",xlab="Years",ylab=names[k],ylim=range(data[,2+k]))
  for (m in 2:M)
    lines(years,data[data[,1]==m,2+k],col=m)
}

y = data[,3]  
X = cbind(1,data[,4:5])
coefs = matrix(0,M,K)
s2s   = rep(0,M)
ses   = matrix(0,M,K)
es    = matrix(0,M,T)
pdf(file="SUR-graph0.pdf",width=20,height=10)
par(mfrow=c(2,5))
for (m in 1:M){
	x1   = X[data[,1]==m,]
	y1   = y[data[,1]==m]
	V    = solve(t(x1)%*%x1)
  coef = V%*%t(x1)%*%y1 	
  fit  = x1%*%coef
  e    = y1 - fit
  s2   = sum(e^2)/(T-3)
  se   = sqrt(s2*diag(V))
	olsm = lm(y[data[,1]==m]~X[data[,1]==m,])
  plot(y[data[,1]==m],fit,xlim=range(y[data[,1]==m],fit),ylim=range(y[data[,1]==m],fit),
       xlab="Real gross",ylab="Fitted real gross",pch=16,main="")
  title(paste(companies[m],"\n s2=",round(s2,2),sep=""))
  abline(0,1,col=2)
  coefs[m,] = coef
  s2s[m]   = s2 
  ses[m,]  = se
  es[m,]   = e
}
dev.off()
cbind(round(coefs[,1],2),round(coefs[,2],4),round(coefs[,3],3))
cbind(round(ses[,1],2),round(ses[,2],4),round(ses[,3],4))

# Seemingly unrelated regression, SUR 
# -----------------------------------
yy = t(matrix(data[,3],M,T,byrow=TRUE))
y1 = matrix(yy,M*T,1)
X1 = NULL
xx = array(0,c(T,M,K*M))
for (t in 1:T){
  x = matrix(0,M,K*M)
  for (m in 1:M)
    x[m,((m-1)*K+1):(m*K)] = X[(m-1)*T+t,]
  xx[t,,] = x
  X1 = rbind(X1,x)
}  

# Simple joint OLS
olsj = solve(t(X1)%*%(X1))%*%t(X1)%*%y1
s2j = sum((y1 - X1%*%olsj)^2)/(M*T-M*K)

# Prior hyperparameters
b0    = rep(0,M*K)
V0    = diag(10000,M*K)
iV0   = solve(V0)
iV0b0 = iV0%*%b0
nu0   = M+3
Phi0  = diag(0.1,M)
iPhi0 = solve(Phi0)
nu1   = nu0+T

# Initial value
beta = matrix(t(coefs),M*K)
Phi  = diag(1/s2s)

# MCMC scheme
set.seed(123456)
M0     = 10000
LAG    = 10
M1     = 1000
niter  = M0+M1*LAG
betas  = matrix(0,niter,M*K)
Sigmas = array(0,c(niter,M,M))
for (iter in 1:niter){
	print(iter)
  # Sampling Phi
  Phi = iPhi0
  for (t in 1:T){
    e = yy[t,]-xx[t,,]%*%beta
    Phi = Phi + e%*%t(e)  
  }
  tcS   = t(chol(solve(Phi)))
  draws = matrix(rnorm(nu1*M),nu1,M)
  Phi   = tcS%*%t(draws)%*%draws%*%t(tcS)
  
  # Sampling beta
  iV = iV0
  iVb = iV0b0
  for (t in 1:T){
    iV = iV + t(xx[t,,])%*%Phi%*%xx[t,,]
    iVb = iVb + t(xx[t,,])%*%Phi%*%yy[t,]
  }
  V = solve(iV)
  b = V%*%iVb
  beta = b + t(chol(V))%*%rnorm(M*K)
  
  betas[iter,] = beta
  Sigmas[iter,,] = solve(Phi)
}

ind = seq(M0+1,niter,by=LAG)
coef1 = matrix(t(coefs),M*K,1)
pdf(file="SUR-graph1.pdf",width=20,height=10)
par(mfrow=c(K,M))
l = 0
for (k in 1:K)
  for (m in 1:M){
  	  l = l + 1
    ts.plot(betas[ind,l],ylim=range(betas[ind,l],coef1[l]),xlab="",ylab="",main=paste(companies[m],"\n BETA",k,sep=""))
    abline(h=coef1[l],col=2)
  }
dev.off()

  
Sigma = diag(s2s)
pdf(file="SUR-graph2.pdf",width=20,height=10)
par(mfrow=c(2,5))
for (m in 1:M){
  ts.plot(sqrt(Sigmas[ind,m,m]),ylim=range(sqrt(Sigmas[ind,m,m]),sqrt(Sigma[m,m])),xlab="",ylab="")
  title(companies[m])
  abline(h=sqrt(Sigma[m,m]),col=2)
}
dev.off()


pdf(file="SUR-graph3.pdf",width=15,height=10)
par(mfrow=c(6,8))
for (i in 1:(M-1))
  for (j in (i+1):M){
    hist(Sigmas[ind,i,j]/sqrt(Sigmas[ind,i,i]*Sigmas[ind,j,j]),prob=TRUE,xlab="",ylab="",main="")
    abline(v=0,col=2)
    title(paste(companies[i],"\n",companies[j]))
  }
dev.off()