###############################################################################
#
# Modeling Petrobras' log-prices, y[t]
#
# y[t] ~ N(alpha + beta*y[t-1],sig2)
#
# Prior: alpha ~ N(a,A), beta ~ N(b,B) e sig2 ~ IG(c,d)
#
# Time period: 1/2/01 - 12/31/13 (3267 obs.)
#
##############################################################################
rm(list=ls())
dig = function(x,a,b){dgamma(1/x,a,b)/x^2}

pdf(file="workedexample-posterior-via-Gibbssampler-AR1-model.pdf",width=8,height=8)

data      = read.csv("petrobras.csv",header=TRUE)
n         = nrow(data)
y         = log(data[,3])

x = y[1:(n-1)]
y = y[2:n]
n = n-1
X = cbind(1,x)

par(mfrow=c(1,1))
plot(x,y,xlab="x: log-price at t-1",ylab="log price at t",pch=16)

# MLE
coefhat = solve(t(X)%*%X)%*%t(X)%*%y
s2hat   = mean((y-X%*%coefhat)^2)

# Prior hyperparameters
a    = 0
A    = 100
b    = 0
B    = 100
c    = 3
d    = (c-1)*(0.01)

########################################################################
#
# Gibbs sampler:
#
# alpha | beta,sig2,X,y   ~ Normal
# beta  | alpha,sig2,X,y  ~ Normal
# sig2  | alpha,beta,X,y  ~ Inverse-Gamma
#
########################################################################
set.seed(12356)
alpha  = a
beta   = b
burnin = 5000
M      = 5000
lag    = 1
niter  = burnin+M*lag
pars  = matrix(0,niter,3)
for (iter in 1:niter){
  p1    = c+(n-1)/2
  p2    = d+sum((y-alpha-beta*x)^2)/2
  sig2  = 1/rgamma(1,p1,p2)

  v     = 1/(1/A+(n-1)/sig2)
  m     = v*(a/A+sum(y-beta*x)/sig2)
  alpha = rnorm(1,m,sqrt(v))

  v     = 1/(1/B+sum(x^2)/sig2)
  m     = v*(b/B+sum((y-alpha)*x)/sig2)
  beta  = rnorm(1,m,sqrt(v))
  pars[iter,] = c(alpha,beta,sig2)
}

pars  = pars[seq(burnin+1,niter,by=lag),]
pars1 = pars

# Checking Gibbs sampler output
par(mfrow=c(3,3))
plot(pars[,1],xlab="iterations",ylab="",main=expression(alpha),type="l")
acf(pars[,1],main="")
hist(pars[,1],xlab="",prob=TRUE,main="");box()
plot(pars[,2],xlab="iterations",ylab="",main=expression(beta),type="l")
acf(pars[,2],main="")
hist(pars[,2],xlab="",prob=TRUE,main="");box()
plot(pars[,3],xlab="iterations",ylab="",main=expression(sigma^2),type="l")
acf(pars[,3],main="")
hist(pars[,3],xlab="",prob=TRUE,main="");box()

par(mfrow=c(1,1))
plot(pars[,1],pars[,2],xlab=expression(alpha),ylab=expression(beta),pch=16)

########################################################################
#
# Gibbs sampler:
#
# alpha,beta | sig2,X,y   ~ Normal
# sig2  | alpha,beta,X,y  ~ Inverse-Gamma
#
########################################################################
m    = matrix(c(a,b),2,1)
C    = diag(c(A,B))
iC   = solve(C)
iCm  = iC%*%m
XtX  = t(X)%*%X
Xty  = t(X)%*%y

# Gibbs sampler
set.seed(12356)
alpha  = a
beta   = b
pars  = matrix(0,niter,3)
for (iter in 1:niter){
  p1    = c+(n-1)/2
  p2    = d+sum((y-alpha-beta*x)^2)/2
  sig2  = 1/rgamma(1,p1,p2)
  C1 = solve(iC + XtX/sig2)
  m1 = C1%*%(iCm + Xty/sig2)
  draw = m1 + t(chol(C1))%*%rnorm(2)
  alpha = draw[1]
  beta  = draw[2]
  pars[iter,] = c(alpha,beta,sig2)
}

pars  = pars[seq(burnin+1,niter,by=lag),]
pars2 = pars

# Checking Gibbs sampler output
par(mfrow=c(3,3))
plot(pars[,1],xlab="iterations",ylab="",main=expression(alpha),type="l")
acf(pars[,1],main="")
hist(pars[,1],xlab="",prob=TRUE,main="");box()
plot(pars[,2],xlab="iterations",ylab="",main=expression(beta),type="l")
acf(pars[,2],main="")
hist(pars[,2],xlab="",prob=TRUE,main="");box()
plot(pars[,3],xlab="iterations",ylab="",main=expression(sigma^2),type="l")
acf(pars[,3],main="")
hist(pars[,3],xlab="",prob=TRUE,main="");box()

par(mfrow=c(1,1))
plot(pars[,1],pars[,2],xlab=expression(alpha),ylab=expression(beta),pch=16)


########################################################################
#
# Gibbs sampler: centering the dependent variable
#
# alpha | beta,sig2,X,y   ~ Normal
# beta  | alpha,sig2,X,y  ~ Normal
# sig2  | alpha,beta,X,y  ~ Inverse-Gamma
#
########################################################################
set.seed(12356)
xbar   = mean(x)
x1     = x-xbar
a      = a-b*xbar
A      = A+B*xbar^2
alpha  = a
beta   = b
burnin = 5000
M      = 5000
lag    = 1
niter  = burnin+M*lag
pars  = matrix(0,niter,3)
for (iter in 1:niter){
  p1    = c+(n-1)/2
  p2    = d+sum((y-alpha-beta*x1)^2)/2
  sig2  = 1/rgamma(1,p1,p2)

  v     = 1/(1/A+(n-1)/sig2)
  m     = v*(a/A+sum(y-beta*x1)/sig2)
  alpha = rnorm(1,m,sqrt(v))

  v     = 1/(1/B+sum(x1^2)/sig2)
  m     = v*(b/B+sum((y-alpha)*x1)/sig2)
  beta  = rnorm(1,m,sqrt(v))
  pars[iter,] = c(alpha,beta,sig2)
}

pars     = pars[seq(burnin+1,niter,by=lag),]
pars[,1] = pars[,1]-pars[,2]*xbar
pars3    = pars

# Checking Gibbs sampler output
par(mfrow=c(3,3))
plot(pars[,1],xlab="iterations",ylab="",main=expression(alpha),type="l")
acf(pars[,1],main="")
hist(pars[,1],xlab="",prob=TRUE,main="");box()
plot(pars[,2],xlab="iterations",ylab="",main=expression(beta),type="l")
acf(pars[,2],main="")
hist(pars[,2],xlab="",prob=TRUE,main="");box()
plot(pars[,3],xlab="iterations",ylab="",main=expression(sigma^2),type="l")
acf(pars[,3],main="")
hist(pars[,3],xlab="",prob=TRUE,main="");box()

par(mfrow=c(1,1))
plot(pars[,1],pars[,2],xlab=expression(alpha),ylab=expression(beta),pch=16)


dev.off()