###########################################################################
# 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("Lect08-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:nrow(dta)) {
      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:60000) {
    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)))
}

print(fres)