####################################################################################################
#
#  Example 7.3 (pages 251-2) – Figure 7.2 (page 252)  - Random walk Metropolis-Hastings algorithm 
#
#
#  Data presented by Ratkowski (1983) on the temporal evolution of the dry weight of onion bulbs.  
#  Gelfand, Dey and Chang (1992) consider 2 possible non-linear models in the form 
#
#                        y[t]~N(f(theta,t),sigma2) 
#
#   where
#
#   Model 1 - Logistic model   :   f(theta,t) = theta[1]/(1+theta[2]*theta[3]*t)
#   Model 2 - Gompertz model   :   f(theta,t) = theta[1] + exp(theta[2]*theta[3]*t}
#
#   with theta[2]>0 and 0 < theta[3] < 1.
#
#          psi[1] = theta[1]
#          psi[2] = log(theta[2])
#          psi[3] = log(theta[3]/(1-theta[3]))
#
#  Non-informative prior: p(psi,sigma2) = 1/sigma2
#
##############################################################################################################
#
# Author : Hedibert Freitas Lopes
#          Graduate School of Business
#          University of Chicago
#          5807 South Woodlawn Avenue
#          Chicago, Illinois, 60637
#          Email : hlopes@ChicagoGSB.edu
#
###############################################################################################################
like.logistic = function(theta){
  -sum(dnorm(y,theta[1]/(1+exp(theta[2])*(exp(theta[3])/(1+exp(theta[3])))^t),exp(theta[4]/2),log=T))
}
post.logistic = function(theta){
  sum(dnorm(y,theta[1]/(1+exp(theta[2])*(exp(theta[3])/(1+exp(theta[3])))^t),exp(theta[4]/2),log=T))
}
like.gompertz = function(theta){
  -sum(dnorm(y,theta[1] + exp(exp(theta[2])*(exp(theta[3])/(1+exp(theta[3])))^t),exp(theta[4]/2),log=T))
}
quant025 = function(x){
  quantile(x,0.025)
}
quant975 = function(x){
  quantile(x,0.975)
}

# Dataset
# -------
n = 15
t = 1:n
y = c(16.08,33.83,65.80,97.20,191.55,326.20,386.87,520.53,590.03,651.92,724.93,699.56,689.96,637.56,717.41)

theta = c(700,36,0.5946)
xs    = seq(1,15,length=100)
f     = theta[1]/(1+theta[2]*theta[3]^xs)
theta = c(theta,mean((y-theta[1]/(1+theta[2]*theta[3]^t))^2))
theta0= c(theta[1],log(theta[2]),log(theta[3]/(1-theta[3])),log(theta[4]))
mle   = nlm(like.logistic,theta0)$estimate
theta.mle = c(mle[1],exp(mle[2]),exp(mle[3])/(1+exp(mle[3])),exp(mle[4]))
round(cbind(theta,theta.mle),2)

#  theta   theta.mle
#  700.00     702.87
#   36.00      84.99
#    0.59       0.50
# 2316.85     595.33

par(mfrow=c(1,1))
plot(t,y,xlab="time",ylab="weight",axes=F)
axis(1)
axis(2)
lines(xs,mle[1]/(1+exp(mle[2])*(exp(mle[3])/(1+exp(mle[3])))^xs),lty=2)

# Random Walk Metropolis (pilot run)
# ---------------------------------
M0     = 1000
M      = 10000
sd     = c(10,0.1,0.05,0.01)
theta  = mle
thetas = theta
for (i in 1:(M0+M-1)){
  theta1 = rnorm(4,theta,sd)
  accept = min(0,post.logistic(theta1)-post.logistic(theta))
  if (log(runif(1))<accept){
      theta = theta1
  }
  thetas = rbind(thetas,t(theta))
}

# Random Walk Metropolis
# ----------------------
M0     = 1000
M      = 1000
sd     = c(10,0.1,0.05,0.01)
V      = var(thetas)
cV = t(chol(V))
theta  = mle
thetas = theta
for (i in 1:(M0+M-1)){
  #theta1 = rnorm(4,theta,sd)
  theta1 = theta + cV%*%rnorm(4)
  accept = min(0,post.logistic(theta1)-post.logistic(theta))
  if (log(runif(1))<accept){
      theta = theta1
  }
  thetas = rbind(thetas,t(theta))
}
thetas[,2] = exp(thetas[,2])
thetas[,3] = exp(thetas[,3])/(1+exp(thetas[,3]))
thetas[,4] = exp(thetas[,4])

# MCMC output
# -----------
names = c("theta1","theta2","theta3","sigma2")
par(mfrow=c(4,2))
for (i in 1:4){
  ts.plot(thetas[(M0+1):(M0+M),i],xlab="iterations",ylab="",main=names[i])
  acf(thetas[(M0+1):(M0+M),i],xlab="lag",ylab="",main=names[i])
}

# Posterior summary and curve fitting
# -----------------------------------
fs = NULL
for (i in (M0+1):(M0+M-1))
  fs = rbind(fs,thetas[i,1]/(1+thetas[i,2]*thetas[i,3]^xs))
f025 = apply(fs,2,quant025)
f050 = apply(fs,2,median)
f975 = apply(fs,2,quant975)

par(mfrow=c(2,2))
for (i in 1:3)
  hist(thetas[(M0+1):(M0+M),i],prob=T,xlab="",ylab="",main=names[i])
plot(t,y,xlab="time",ylab="weight",axes=F)
axis(1)
axis(2)
lines(xs,f025)
lines(xs,f050)
lines(xs,f975)