###########################################################################
# Charles B. Moss                                                         #
# Food and Resource Economics Department                                  #
# University of Florida                                                   #
# Email: cbmoss@ufl.edu   Websites: http://ricardo.ifas.ufl.edu           #
#                                   http://www.charlesbmoss.com:8080      #
###########################################################################
# Introduction to Bayesian Econometrics                                   #
# Lecture 8 - Simulation and Empirical Loss Functions                     #
###########################################################################
# Estimation of a linear equation using a Metropolis Hastings Monte Carlo #
# drawing system.                                                         #
# Load the 1960 through 1969 KLEM capital share data                      #

library(rgl)

dta <- read.csv("Lect05-2013.csv")

# dta[,1] - Share of cost spent on capital                                #
# dta[,2] - Logarithm of agricultural output                              #
# dta[,3] - Logarithm of capital input price                              #
# dta[,4] - Logarithm of labor input price                                #
# dta[,5] - Logarithm of energy input price                               #
#                                                                         #
# Setup the X matrix (independent variables)                              #

xx <- cbind(matrix(1,nrow=nrow(dta),ncol=1),as.matrix(dta[,2:5]))

###########################################################################
# Prior distributions                                                     #
#  Mean of the beta vector                                                #

bprior <- rbind(mean(dta[,1]),0,0,0,0)
# msep   <- (1/nrow(dta))*t(dta[,1]-xx%*%bprior)%*%(dta[,1]-xx%*%bprior)
# vprior <- 2.0*msep[1,1]*solve(t(xx)%*%xx)
# vprior <- diag(diag(vprior),5,5)
vprior <- 1.50*diag(1,5,5)

like <- function(b) {
    ll <- 0
    for (j in 1:10) {
      ll <- ll - t(dta[j,1]-xx[j,]%*%b)%*%(dta[j,1]-xx[j,]%*%b)  }
    return (exp(ll/2)) }

###########################################################################
# bdraw - Draw from prior distribution                                    #
# nprob - New probability (value of the likelihood function)              #
# res   - Result matrix                                                   #

mcmc <- cbind(0.00,0.05,0.10,0.25,0.50,0.75)

for (m in mcmc) {
  bdraw <- bprior + chol(vprior)%*%rnorm(5,mean=0,sd=1)
  nprob <- like(bdraw)
  while (nprob == 0) {
    bdraw <- bprior + chol(vprior)%*%rnorm(5,mean=0,sd=1)
    nprob <- like(bdraw)}
  res   <- cbind(1,nprob,t(bdraw),nprob[1,1]*t(bdraw))
  best  <- nprob[1,1]*bdraw

  for (i in 2:45000) {
    odraw <- bdraw
    oprob <- nprob
    bdraw <- bprior + chol(vprior)%*%rnorm(5,mean=0,sd=1)
    nprob <- like(bdraw)
    lrat <- nprob/oprob
###########################################################################
#  Metropolis Hastings Monte Carlo section                                #
    if (lrat <= m || is.nan(lrat)) {                                      
      while (nprob == 0 || udraw==0) {
        udraw <- runif(1)
        tdraw <- udraw*odraw + (1-udraw)*bdraw 
        nprob <- like(tdraw) }
      bdraw <- tdraw
     }
###########################################################################
    res <- rbind(res,cbind(i,nprob,t(bdraw),nprob[1,1]*t(bdraw)))
    best <- best + nprob[1,1]*bdraw
    }

  best <- best / sum(res[,2])
  vhat <- matrix(0,nrow=5,ncol=5)
  for (i in 1:nrow(res)) {
    vhat <- vhat + res[i,2]/sum(res[,2])*
         t(res[i,cbind(3:7)]-t(best))%*%(res[i,cbind(3:7)]-t(best)) }
  if (m==0.0) tres <- res[,2] else tres <- cbind(tres,res[,2])
  if (m==0.0) fres <- rbind(m,best,as.matrix(sqrt(diag(vhat))),sum(res[,2]),sum(res[,2])/nrow(res)) else
       fres <- cbind(fres,rbind(m,best,as.matrix(sqrt(diag(vhat))),sum(res[,2]),sum(res[,2])/nrow(res)))
}

plot(res[,cbind(4,5)],xlab="ln(Output)",
     ylab="ln(Capital Price)")

plot3d(res[,cbind(4,5,6)],xlab="ln(Output)",ylab="ln(Capital Price}",
       zlab="ln(Labor Price)")

print(fres)

bols <- solve(t(xx)%*%xx)%*%(t(xx)%*%dta[,1])
mse  <- 1/(nrow(dta)-5)*t(dta[,1]-xx%*%bols)%*%(dta[,1]-xx%*%bols)
vols <- mse[1,1]*solve(t(xx)%*%xx)
print(cbind(bols,sqrt(diag(vols))))