#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter6: Section6.2.1
#
# Bayesian analysis of the ordered probit model
#
# Author : Tomohiro Ando
#
################################################################################
library(MCMCpack)


#Note: the primary heading in "Ratingdata" should be deleted when you read the dataset

D <- matrix(scan("Ratingdata.txt"),ncol=4,byrow=T)
D <- as.data.frame(D)
names(D) <- c("ROA","logSales","CF-Sales ratio","Ratings")
n <- nrow(D)


#Model 1 
#posterior <- MCMCoprobit(Ratings~., data = D, thin=5, burnin = 1000, mcmc = 1000, b0=0, B0 = 0.0001, mcmc.method="AC") 


plot(posterior)
summary(posterior)


Marloglike <- 0

Y <- matrix(0,n,4)
for(i in 1:n){Y[i,D[i,4]] <-1}

for(j in 1:nrow(posterior)){

X <- as.matrix(cbind(1,D[,1:3]))
Beta <- as.vector(posterior[j,1:4])
O <- X%*%Beta

G  <- as.vector(posterior[j,-(1:4)])
p1 <- pnorm(0-O)-0
p2 <- pnorm(G[1]-O)-pnorm(0-O)
p3 <- pnorm(G[2]-O)-pnorm(G[1]-O)
p4 <- 1-pnorm(G[2]-O)

P <- cbind(p1,p2,p3,p4)

Loglike <- sum( log(P^Y) )

Marloglike <- Marloglike + 1/Loglike

}

print(Marloglike)



#Model 2 

posterior <- MCMCoprobit(Ratings~ROA+logSales, data = D, thin=5, burnin = 1000, mcmc = 1000, b0=0, B0 = 0.0001, mcmc.method="AC") 


plot(posterior)
summary(posterior)

Marloglike <- 0

Y <- matrix(0,n,4)
for(i in 1:n){Y[i,D[i,4]] <-1}

for(j in 1:nrow(posterior)){

X <- as.matrix(cbind(1,D[,1:2]))
Beta <- as.vector(posterior[j,1:3])
O <- X%*%Beta

G  <- as.vector(posterior[j,-(1:3)])
p1 <- pnorm(0-O)-0
p2 <- pnorm(G[1]-O)-pnorm(0-O)
p3 <- pnorm(G[2]-O)-pnorm(G[1]-O)
p4 <- 1-pnorm(G[2]-O)

P <- cbind(p1,p2,p3,p4)

Loglike <- sum( log(P^Y) )

Marloglike <- Marloglike + 1/Loglike

}

print(Marloglike)