###########################################################################
# 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 7 - Simulation and Empirical Loss Functions                     #
###########################################################################
# Starting with a know sample - a normal variable with mean of 5 and      #
# variance of 2.5.                                                        #
#                                                                         #

rdraw <- rnorm(500,mean=5, sd=sqrt(2.5))

###########################################################################
# Next, we want to set up a likelihood function based on a draw from the  #
# prior                                                                   #

like <- function(a) {
   ll <- 0
   for (i in 1:ndraw) {
     ll <- ll + (x[i]-a)^2 }
   return (exp(-ll/(2*10))) }

###########################################################################
# Next, let's draw 10 values from the prior distribution and compute the  #
# likelihood function                                                     #


x <- rdraw[1:10]
ndraw <- 10

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- cbind(10,cmean,bmean,bmean/cmean)

# 20 values                                                               #

x <- rdraw[1:20]
ndraw <- 20

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(20,cmean,bmean,bmean/cmean))

# 30 values                                                               #

x <- rdraw[1:30]
ndraw <- 30

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(30,cmean,bmean,bmean/cmean))

# 40 values                                                               #

x <- rdraw[1:40]
ndraw <- 40

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(40,cmean,bmean,bmean/cmean))

# 50 values                                                               #

x <- rdraw[1:50]
ndraw <- 50

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(50,cmean,bmean,bmean/cmean))

# 75 values                                                               #

x <- rdraw[1:75]
ndraw <- 75

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(75,cmean,bmean,bmean/cmean))

# 100 values                                                              #

x <- rdraw[1:100]
ndraw <- 100

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(100,cmean,bmean,bmean/cmean))

# 150 values                                                              #

x <- rdraw[1:150]
ndraw <- 150

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(150,cmean,bmean,bmean/cmean))

# 200 values                                                              #

x <- rdraw[1:200]
ndraw <- 200

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(200,cmean,bmean,bmean/cmean))

# 300 values                                                              #

x <- rdraw[1:300]
ndraw <- 300

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(300,cmean,bmean,bmean/cmean))

# 400 values                                                              #

x <- rdraw[1:400]
ndraw <- 400

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(400,cmean,bmean,bmean/cmean))

# 500 values                                                              #

x <- rdraw[1:500]
ndraw <- 500

for (s in 1:2000) {
   a <- rnorm(1,mean=2.5,sd=1.0)
   if (s==1) res <- cbind(a,like(a),a*like(a)) else 
               res <- rbind(res,cbind(a,like(a),a*like(a)))
   }

bmean <- sum(res[,3])/sum(res[,2])
cmean <- mean(x)
rres <- rbind(rres,cbind(500,cmean,bmean,bmean/cmean))


print(rres)