####### rjags code to fit the Lognormal model to the income data #######

# unbuffer the output

# install rjags if you haven't already done so

library( rjags ) 

# change directory to the place where the income data file is

y <- scan( 'U.S.-family-income-2009.txt' )

income.data <- list( y = y, n = length( y ) )

income.lognormal.initial.values <- 
  list( mu = 0.0, log.tau = 0.0, y.new = mean( y ) )

A <- 2

income.lognormal.run.1 <- jags.model( file = 
  'ams-206-income-lognormal-analysis-rjags-model.txt', 
  data = income.data, inits = income.lognormal.initial.values, 
  n.chains = A )

# choose your own random number seed, different from
# my choice of 42, in the next line of code

set.seed( 42 )

n.burnin <- 1000

system.time( 

  update( income.lognormal.run.1, n.iter = n.burnin )

)

n.monitor <- 100000

system.time( 

  income.lognormal.run.1.results <- jags.samples( income.lognormal.run.1, 
    variable.names = c( 'mu', 'sigma.squared', 'theta', 'y.new' ), 
    n.iter = n.monitor )

)

mu.star <- income.lognormal.run.1.results$mu

sigma.squared.star <- income.lognormal.run.1.results$sigma.squared

theta.star <- income.lognormal.run.1.results$theta

y.new.star <- income.lognormal.run.1.results$y.new

n.monitor <- A * n.monitor

########################################################################

par( mfrow = c( 2, 2 ) )

plot( 1:n.monitor, mu.star, type = 'l', xlab = 'Iteration Number',
  ylab = 'mu' )

plot( density( mu.star ), xlab = 'mu', main = '' )

acf( as.ts( mu.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

pacf( as.ts( mu.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

c( mean( mu.star ), sd( mu.star ), quantile( mu.star, 
  probs = c( 0.025, 0.975 ) ) )

print( rho.1.hat.mu.star <- cor( mu.star[ 1:( n.monitor - 1 ) ],
  mu.star[ 2:n.monitor ] ) )

print( se.mu.star.mean <- ( sd( mu.star ) / sqrt( n.monitor ) ) *
  sqrt( ( 1 + rho.1.hat.mu.star ) / ( 1 - rho.1.hat.mu.star ) ) )

par( mfrow = c( 1, 1 ) )

########################################################################

par( mfrow = c( 2, 2 ) )

plot( 1:n.monitor, sigma.squared.star, type = 'l', 
  xlab = 'Iteration Number', ylab = 'sigma squared' )

plot( density( sigma.squared.star ), xlab = 'sigma.squared', main = '' )

acf( as.ts( sigma.squared.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

pacf( as.ts( sigma.squared.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

c( mean( sigma.squared.star ), sd( sigma.squared.star ), 
  quantile( sigma.squared.star, probs = c( 0.025, 0.975 ) ) )

print( rho.1.hat.sigma.squared.star <- cor( 
  sigma.squared.star[ 1:( n.monitor - 1 ) ],
  sigma.squared.star[ 2:n.monitor ] ) )

print( se.sigma.squared.star.mean <- ( sd( sigma.squared.star ) / 
  sqrt( n.monitor ) ) * sqrt( ( 1 + rho.1.hat.sigma.squared.star ) / 
  ( 1 - rho.1.hat.sigma.squared.star ) ) )

par( mfrow = c( 1, 1 ) )

########################################################################

par( mfrow = c( 2, 2 ) )

plot( 1:n.monitor, theta.star, type = 'l', xlab = 'Iteration Number',
  ylab = 'theta' )

plot( density( theta.star ), xlab = 'theta', main = '' )

acf( as.ts( theta.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

pacf( as.ts( theta.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

c( mean( theta.star ), sd( theta.star ), quantile( theta.star, 
  probs = c( 0.025, 0.975 ) ) )

print( rho.1.hat.theta.star <- cor( theta.star[ 1:( n.monitor - 1 ) ],
  theta.star[ 2:n.monitor ] ) )

print( se.theta.star.mean <- ( sd( theta.star ) / sqrt( n.monitor ) ) *
  sqrt( ( 1 + rho.1.hat.theta.star ) / ( 1 - rho.1.hat.theta.star ) ) )

par( mfrow = c( 1, 1 ) )

########################################################################

par( mfrow = c( 2, 2 ) )

plot( 1:n.monitor, y.new.star, type = 'l', xlab = 'Iteration Number',
  ylab = 'y.new' )

plot( density( y.new.star ), xlab = 'y.new', main = '' )

acf( as.ts( y.new.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

pacf( as.ts( y.new.star[ 1:n.monitor ], start = 1, 
  end = n.monitor, frequency = 1 ), lag.max = 10, main = '' )

par( mfrow = c( 1, 1 ) )

quantile( income.data$y, 
  probs = c( 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99 ) )

quantile( y.new.star, 
  probs = c( 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99 ) )

qqplot( income.data$y, y.new.star, xlab = 'income data',
  ylab = 'Predictive distribution from Lognormal model' )

abline( 0, 1, lwd = 2, col = 'red' )

c( mean( y.new.star ), sd( y.new.star ), quantile( y.new.star, 
  probs = c( 0.025, 0.975 ) ) )

print( rho.1.hat.y.new.star <- cor( y.new.star[ 1:( n.monitor - 1 ) ],
  y.new.star[ 2:n.monitor ] ) )

print( se.y.new.star.mean <- ( sd( y.new.star ) / sqrt( n.monitor ) ) *
  sqrt( ( 1 + rho.1.hat.y.new.star ) / ( 1 - rho.1.hat.y.new.star ) ) )

c( mean( income.data$y ), sd( income.data$y ) )

c( mean( y.new.star ), sd( y.new.star ) )

########################################################################

system.time( 

  lognormal.dic.summary <- dic.samples( income.lognormal.run.1, 
    n.iter = 100000, type = 'pD' )

)

print( lognormal.mean.deviance <- sum( lognormal.dic.summary$deviance ) )

print( lognormal.complexity <- sum( lognormal.dic.summary$penalty ) )

print( lognormal.dic <- lognormal.mean.deviance + lognormal.complexity )

########################################################################