#install.packages("MCMCpack")
library(MCMCpack)

##read data
data1=read.table(file='D:\\CFE09_tutorial\\R_codes\\BF\\usedcar.txt')
##name data
y =matrix(data1[,1])
y=y/1000   # (one thousand dollar for one unit)
x1=matrix(data1[,2])
x2=matrix(data1[,3])
x3=matrix(data1[,4])
x4=matrix(data1[,5])
x5=matrix(data1[,6])


##########################################
##  Model selection by Bayes Factor


svar<- var(y)
## MCMC regress with burn-in 2000
## "Laplace" in which case the Laplace approximation (see Kass and Raftery, 1995) is used,
## and "Chib95" in which case the method of Chib (1995) is used.
model1<-MCMCregress(y~x1+x2, burnin=2000,
b0=c(0, 0, 0),
B0=c(1e-2, 1e-2, 1e-2), c0=5, d0=svar*3, marginal.likelihood="Chib95", mcmc=10000, verbose=2000)


##MCMC regress with burn-in 2000
model2<-MCMCregress(y~x1+x2+x3+x4, burnin=2000,
b0=c(0, 0, 0, 0, 0),
B0=c(1e-2, 1e-2, 1.0e-2, 1.0e-2,1.0e-2), c0=5, d0=svar*3, marginal.likelihood="Chib95", mcmc=10000, verbose=2000)


##MCMC regress with burn-in 2000
model3<-MCMCregress(y~x1+x2+x3+x4+x5, burnin=2000,
b0=c(0, 0, 0, 0, 0, 0),
B0=c(1e-2, 1e-2, 1e-2, 1e-2, 1e-2, 1e-2), c0=5, d0=svar*3, marginal.likelihood="Chib95", mcmc=10000, verbose=2000)


BF <- BayesFactor(model1, model2, model3)
print(BF)



#The matrix of Bayes Factors is:
#       model1 model2 model3
#model1  1.000  0.589   3.29
#model2  1.697  1.000   5.59
#model3  0.304  0.179   1.00

#The matrix of the natural log Bayes Factors is:
#       model1 model2 model3
#model1  0.000 -0.529   1.19
#model2  0.529  0.000   1.72
#model3 -1.192 -1.720   0.00