##################################################################################################
# A FUNCTION TO GENERATE RANDOM MULTIVARIATE GAUSSIANS
##################################################################################################

rmultinorm <- function(num.vals, mu.vec, vcmat, tol = 1e-08)  {
    k <- ncol(vcmat)
    if(length(mu.vec)!=k)
        stop(paste("rmultinorm error: rmultinorm: mu.vec vector wrong length:",length(mu.vec)))
    if(max(abs(vcmat - t(vcmat))) > tol)
        stop("rmultinorm error: variance-covariance matrix not symmetric")
    vs <- svd(vcmat)
    vcsqrt <- t(vs$v %*% (t(vs$u) * sqrt(vs$d)))
    ans.mat <- sweep(matrix(rnorm(num.vals*k), nrow = num.vals) %*% vcsqrt,2,mu.vec,"+")
    dimnames(ans.mat) <- list(NULL, dimnames(vcmat)[[2]])
    return(ans.mat)
}

##################################################################################################
# A FUNCTION TO DEAL WITH MACHINE ROUNDING/TRUNCATION ISSUES
##################################################################################################

make.symmetric <- function(in.mat)  {
         for (i in 1:nrow(in.mat))
                 for (j in 1:ncol(in.mat))
                         if (i > j) in.mat[i,j] <- in.mat[j,i]
         return(in.mat)
}

##################################################################################################
# SETUP DATA, HYPERPARAMETERS AND START VALUES
##################################################################################################

library(msm)                            # FOR THE TRUNCATED NORMAL GENERATOR: rtnorm
library(foreign)                        # TO IMPORT THE STATA FILE
library(survival)                       # FOR THE survreg FUNCTION
library(LearnBayes)			# FOR THE rigamma FUNCTION
norr.df <- read.table("http://jgill.wustl.edu/data/norr.dat",header=TRUE)

# Current sentencing policy   		ds9395p			PRIOR DIRECTION
# Past execution rate          		ep4089n			+
# Political culture                   	polcul			+
# Current opinion                 	d8892r2			+
# Citizen ideology                	id8892m2		+
# Murder rate                         	murder90		+
# Catholic                          	catholic		-
# Black                               	black90			-
# Urban                                	urban90			-
# Past laws                       	law3689			+
# Past opinion                 		favdth36		+

##################################################################################################
# RUN A STANDARD OLS AND A TOBIT MODEL TO REPLICATE NORRANDER, PAGE 788-9 
##################################################################################################

# MODEL B1
  norr.ols <- lm(ds9395p ~ ep4089n+polcul+d8892r2,data=norr.df,na.action="na.fail")
  norr.tob <- survreg(Surv(ds9395p,ds9395p>0,type='left') ~ 
	                 ep4089n+polcul+d8892r2, dist='gaussian',data=norr.df)
# MODEL B2
  norr.ols <- lm(ds9395p ~ ep4089n+polcul+d8892r2+id8892m2,data=norr.df,na.action="na.fail")
  norr.ols <- lm(ds9395p ~ ep4089n+polcul+id8892m2,data=norr.df,na.action="na.fail")
# MODEL B3
  norr.ols <- lm(ds9395p ~ ep4089n+polcul+d8892r2+murder90,data=norr.df,na.action="na.fail")
# MODEL B4
  norr.ols <- lm(d8892r2 ~ catholic+black90+urban90+law3689+ep4089n,
	      data=norr.df,na.action="na.fail")
# MODEL 5
  norr.ols <- lm(d8892r2 ~ catholic+black90+urban90+law3689+murder90+id8892m2,
	      data=norr.df,na.action="na.fail")
# MODEL 6
  norr.ols <- lm(ds9395p ~ ep4089n+polcul+d8892r2+favdth36, data=norr.df,na.action="na.fail")

##################################################################################################
# SETUP DATA STRUCTURES FOR GIBBS
##################################################################################################

attach(norr.df)
Y <- c(ds9395p);  num.zeros <- 15
X <- cbind(rep(1,nrow(norr.df)), ep4089n, polcul, d8892r2, id8892m2, murder90)
detach(norr.df)
dimnames(X)[[2]] <- c("Constant","Past Rates","Political Culture","Current Opinion","Ideology",
	"Murder Rate")

##################################################################################################
# RUN BIVARIATE REGRESSIONS TO GET PRIOR MEANS FOR BETAS
##################################################################################################

beta0 <- rep(NA,ncol(X)); beta0[1] <- 1
for (i in 2:ncol(X))   beta0[i] <- summary(lm(Y~X[,i]))$coef[2,1]

##################################################################################################
# SET PARAMATERS FOR GIBBS
##################################################################################################

gamma0 <- 300; gamma1 <- 100 					# HYPERPARAMETERS
B0 <- 0.02 							# SCALING ON SIGMA^0
mc.start <- 1; mc.stop <- 50000					# MCMC CONTROL
norr.tob <- survreg(Surv(Y,Y>0,type='left') ~ X[,-1], dist='gaussian',data=norr.df)
B.mat  <- matrix(summary(norr.tob)$coef,nrow=1)			# STARTING POINTS
B.mat[1,] <- jitter(B.mat[1,],20)				# ALTER VALUES FOR GR DIAGNOSTIC
s.sq   <- matrix(rep(10,6), nrow=1)				# PRIOR on SIGMA^2 OF BETA
y.star <- Y							# LATENT DATA START

##################################################################################################
# GIBBS SAMPLER
##################################################################################################

for (i in mc.start:mc.stop)  {
    for (j in 1:length(Y)) {
        if (Y[j] == 0)  y.star[j] <- rtnorm(1,X[j,]%*%B.mat[i,],s.sq[i,],lower=-Inf,upper=0)
    }
    delta <- t(y.star - X %*% B.mat[i,]) %*% (y.star - X %*% B.mat[i,])
    s.temp <- rep(Inf,ncol(s.sq)); while(!is.finite(sum(s.temp))) 
    		       	s.temp <- rigamma(ncol(s.sq), (gamma0+length(Y))/2, (gamma1+delta)/2 )
    s.sq <- rbind(s.sq,s.temp) 
    B.hat <- solve( B0 + t(X)%*%X ) %*% (B0*beta0 + t(X)%*%y.star)
    B.var <-  make.symmetric( solve( s.sq[(i+1)]^(-1)*B0 + s.sq[(i+1)]^(-1)*t(X)%*%X ) )
    B.mat <- rbind( B.mat, rmultinorm(1,B.hat,B.var,tol=1e-06) )
    if (i %% 100 == 0)  print(paste("iteration:",i))
}

##################################################################################################
# RESULTS
##################################################################################################

B.mat6 <- B.mat[40001:50000,]	# TAKE ONLY THE LAST 10,000 DRAWS OF 50,000
save.image()

write.table(B.mat2,"Book.Bayes.III/Example.Tobit/mcmc.mat.sep.15.2011",sep=" ")
write.table(B.mat3,"Book.Bayes.III/Example.Tobit/mcmc.mat3.oct.3.2011",sep=" ")
write.table(B.mat4,"Book.Bayes.III/Example.Tobit/mcmc.mat4.oct.3.2011",sep=" ")
write.table(B.mat5,"Book.Bayes.III/Example.Tobit/mcmc.mat5.oct.3.2011",sep=" ")
write.table(B.mat6,"Book.Bayes.III/Example.Tobit/mcmc.mat6.oct.3.2011",sep=" ")
B.mat2 <- read.table("Book.Bayes.III/Example.Tobit/mcmc.mat.sep.15.2011",sep=" ")
B.mat2.mcmc <- mcmc(B.mat2[1:10000,])
B.mat3.mcmc <- mcmc(B.mat3[1:10000,])
B.mat4.mcmc <- mcmc(B.mat4[1:10000,])
B.mat5.mcmc <- mcmc(B.mat5[1:10000,])
B.mat6.mcmc <- mcmc(B.mat6[1:10000,])
dimnames(B.mat6.mcmc) <- dimnames(B.mat5.mcmc) <- dimnames(B.mat4.mcmc) <- dimnames(B.mat3.mcmc) <- dimnames(B.mat2.mcmc)

B.list <- mcmc.list(B.mat2.mcmc,B.mat3.mcmc,B.mat4.mcmc,B.mat5.mcmc,B.mat6.mcmc)
B.list <- mcmc.list("B.mat2.mcmc"=B.mat2.mcmc, "B.mat3.mcmc"=B.mat3.mcmc, "B.mat4.mcmc"=B.mat4.mcmc,
	       "B.mat5.mcmc"=B.mat5.mcmc, "B.mat6.mcmc"=B.mat6.mcmc)