########### Get the data ready for analysis ###########

library(AER)
data("RecreationDemand")
dat<-RecreationDemand
dat$incomeC <- dat$income - mean(dat$income)
hist(dat$trips, main = "", 
     xlab = "Number of recreational boating trips")

############### load packages ###############
library(runjags)
library(kableExtra)


############### analyses ###############

###### Poisson ######
Poisson_Model <- "model{
   ## Likelihood ##
   for (i in 1:N){
    Y[i] ~ dpois(lambda[i])
    log(lambda[i]) <- beta0 + beta1*X[i]
   }
   
   ## priors for coefficients
    beta0 ~ dnorm(0, 1/10000)
    beta1 ~ dnorm(0, 1/10000)  
    
   ## exponentiate the paramters 
   exp_beta0 <- exp(beta0)
   exp_beta1 <- exp(beta1)
}"

inits <- list(list(beta0 = rnorm(1, 0, 0.1), 
              beta1 = rnorm(1, 0, 0.1)),
              list(beta0 = 1, beta1 = 1),
              list(beta0 = -1, beta1 = -1))

data <- list(Y = dat$trips, X = dat$incomeC, 
             N = nrow(dat))

Pois_Est <- run.jags(model = Poisson_Model, 
  monitor = c("beta0", "beta1", "exp_beta0", 
  "exp_beta1"), data = data, n.chains = 3, 
  inits = inits, method = "simple", adapt = 1000, 
  burnin = 3000, sample = 10000)


res11<-cbind(round(Pois_Est$HPD[,c(1,3)],2), 
  round(Pois_Est$summary$statistics[,1:2],2), 
  round(Pois_Est$psrf$psrf[,1],4))
colnames(res11) <- c("Lower 95", "Upper 95", 
  "Mean", "SD", "psrf")
kable(res11, caption = "Poisson runjags  
  Output", "simple")
par(mfrow = c(1, 2)) 
plot(Pois_Est, plot.type = "trace")


###### ZIP ######
ZIP_model <- "model{
 ## likelihood
 for (i in 1:N){
  Y[i] ~ dpois(lambda[i])
  lambda[i] <- W[i]*mu[i]
  W[i] ~ dbern(1-p[i])
  log(mu[i]) <- beta0 + beta1*X[i]
  logit(p[i]) <- gamma0 + gamma1*B[i]
 }
 
 ## prior 
  beta0 ~ dnorm(0, 1/10000)
  beta1 ~ dnorm(0, 1/10000) 
  gamma0 ~ dnorm(0, 1/10000)
  gamma1 ~ dnorm(0, 1/10000) 
  
  ## exponentiate the paramters 
  exp_beta0 <- exp(beta0)
  exp_beta1 <- exp(beta1)
  exp_gamma0 <- exp(gamma0)
  exp_gamma1 <- exp(gamma1)
}"

W <- dat$trips
W[dat$trips>0] <- 1

inits <- list(list(beta0 = rnorm(1, 0, 0.1), 
              beta1 = rnorm(1, 0, 0.1), 
              gamma0 = rnorm(1, 0, 0.1), 
              gamma1 = rnorm(1, 0, 0.1),  
  	      W = W), list(beta0 = 1, beta1 = 1, 
	      gamma0 = 1, gamma1 = 1, W = W),
              list(beta0 = -1, beta1 = -1,  
  	      gamma0 = 1, 
              gamma1 = 1, W = W))

data <- list(Y = dat$trips, X = dat$incomeC, 
             B = dat$incomeC, N = nrow(dat))

ZIP_Est <- run.jags(model = ZIP_model, monitor = 
   c("beta0", "beta1", "gamma0", "gamma1",  
   "exp_beta0", "exp_beta1", "exp_gamma0",
   "exp_gamma1"), data = data, n.chains = 3,
    inits = inits, method = "simple", adapt = 1000, 
    burnin = 3000, sample = 10000,
    keep.jags.files = T, tempdir = T)
res2<-cbind(round(ZIP_Est$HPD[,c(1,3)],2), 
  round(ZIP_Est$summary$statistics[,1:2],2), 
  round(ZIP_Est$psrf$psrf[,1],4))
colnames(res2) <- c("Lower 95", "Upper 95", "Mean", 
  "SD", "psrf")
res22<-round(exp(res2[,c(1, 2, 3)]),2)

kable(res22, caption = "ZIP runjags Exponentiated 
  Output", "simple")
par(mfrow = c(1, 2)) 
plot(ZIP_Est, plot.type = "trace")


###### Hurdle ######
hurdle_model <- "model{
 ## likelihood
 C <- 10000
  for (i in 1:N){
  zeros[i] ~ dpois(-ll[i] + C)
  truncPois[i] <- Y[i]*log(mu[i]) - mu[i] -  
  	(log(1-exp(-mu[i])) + logfact(Y[i]))
                  
  l1[i] <- (1-z[i])*log(p[i])
  l2[i] <- z[i]*(log(1-p[i]) + truncPois[i])
  ll[i] <- l1[i] + l2[i]
  
  log(mu[i]) <- beta0 + beta1*X[i]
  logit(p[i]) <- gamma0 + gamma1*B[i]
 }
 
 ## prior 
  beta0 ~ dnorm(0, 1/10000)
  beta1 ~ dnorm(0, 1/10000) 
  gamma0 ~ dnorm(0, 1/10000)
  gamma1 ~ dnorm(0, 1/10000) 
  
 ## exponentiate the paramters 
  exp_beta0 <- exp(beta0)
  exp_beta1 <- exp(beta1)
  exp_gamma0 <- exp(gamma0)
  exp_gamma1 <- exp(gamma1)
}"

z<-dat$trips
z[dat$trips > 0] <- 1

inits <- list(list(beta0 = rnorm(1, 0, 0.1),  
  		beta1 = rnorm(1, 0, 0.1), 
		gamma0 = rnorm(1, 0, 0.1),  
  		gamma1 = rnorm(1, 0, 0.1)), 
  		list(beta0 = 1, beta1 = 1,  
  		gamma0 = 1, gamma1 = 1), 
  		list(beta0 = -1, beta1 = -1,  
  		gamma0 = 1,  
  		gamma1 = 1))

data <- list(Y = dat$trips, X = dat$incomeC,   
 		B = dat$incomeC,   
  		N = nrow(dat), z = z,  
  		zeros = rep(0, nrow(dat)))

hurdle_Est <- run.jags(model = hurdle_model,  
  		monitor = c("beta0", "beta1", 
  		 "gamma0", "gamm a1", "exp_beta0",  
  		"exp_beta1", "exp_gamma0", "exp_gamma1"),  
  		data = data, n.chains = 3,  
  		inits = inits, method = "simple",  
  		adapt = 1000,  burnin = 3000,  
  		sample = 10000,  
  		keep.jags.files = T, tempdir = T)

res33<-cbind(round(hurdle_Est$HPD[,c(1,3)],2), 
          round(hurdle_Est$summary$statistics[,1:2],2), 
          round(hurdle_Est$psrf$psrf[,1],4)) 
  	  colnames(res33) <- c("Mean", "SD",  
  	  "Lower 95", "Upper 95","psrf")
res33<-round(exp(res33[,c(1, 2, 3)]),2)
kable(res33, caption = "Hurdle runjags  
  Exponentiated Output", "simple")
plot(hurdle_Est, plot.type = "trace")