# change directory to wherever you stored the file
#   lognormal-model.txt

# setwd( "C:/Teaching/AMS-206-132/Winter-2017" )

library( rjags )

income.data <- list( y = c( 
0.403, 0.606, 0.623, 0.984, 1.265, 2.182, 2.31, 2.343, 2.427, 
2.472, 2.726, 2.799, 3.014, 3.156, 3.494, 3.526, 3.697, 4.029, 
4.113, 4.49, 4.717, 4.836, 5.046, 5.521, 5.821, 5.936, 6.017, 
6.037, 6.173, 6.182, 6.536, 7.056, 7.27, 8.386, 8.422, 8.565, 
8.803, 8.807, 9.1, 9.117, 9.538, 9.588, 9.697, 9.708, 9.922, 
10.073, 10.224, 10.488, 10.496, 10.753, 11.371, 11.428, 11.627, 
11.72, 11.911, 11.948, 11.973, 12.014, 12.104, 12.184, 12.23, 
12.368, 12.486, 12.487, 12.491, 12.562, 13.151, 13.18, 13.203, 
13.373, 13.589, 13.592, 13.667, 13.96, 14.061, 14.383, 14.658, 
14.858, 15.055, 15.255, 15.394, 15.603, 15.924, 15.969, 16.298, 
16.579, 16.58, 17.036, 17.142, 17.226, 17.696, 17.826, 17.849, 
17.957, 18.276, 18.836, 18.846, 18.872, 18.927, 19.092, 19.118, 
19.351, 19.622, 19.883, 19.888, 19.915, 19.946, 20.071, 20.094, 
20.142, 20.248, 20.557, 20.819, 20.898, 20.94, 21.021, 21.052, 
21.068, 21.154, 21.178, 21.209, 21.378, 21.513, 21.544, 21.596, 
21.643, 21.857, 22.031, 22.24, 22.243, 22.494, 22.56, 22.871, 
22.895, 22.948, 23.281, 23.387, 23.75, 24.048, 24.05, 24.258, 
24.26, 24.39, 24.441, 24.676, 24.792, 24.879, 24.9, 25.215, 25.302, 
25.302, 25.361, 25.373, 25.55, 25.6, 25.679, 25.681, 25.684, 
25.724, 25.83, 26.182, 26.202, 26.485, 26.643, 26.666, 26.981, 
27.06, 27.11, 27.143, 27.154, 27.193, 27.302, 27.33, 27.372, 
27.374, 27.661, 27.904, 28.059, 28.314, 28.43, 28.565, 28.611, 
29.014, 29.037, 29.045, 29.085, 29.102, 29.755, 29.92, 29.955, 
29.976, 30.051, 30.055, 30.348, 30.616, 30.79, 30.803, 30.909, 
30.912, 30.915, 31.105, 31.275, 31.367, 31.449, 31.656, 31.89, 
31.915, 31.927, 31.97, 31.979, 32.044, 32.142, 32.344, 32.52, 
32.72, 32.72, 32.9, 32.998, 33.272, 33.274, 33.482, 33.558, 33.892, 
34.002, 34.035, 34.083, 34.103, 34.164, 34.421, 34.455, 34.498, 
34.962, 35.031, 35.152, 35.247, 35.393, 35.527, 35.665, 35.786, 
35.955, 36.05, 36.091, 36.094, 36.151, 36.188, 36.193, 36.296, 
36.669, 36.759, 36.78, 36.863, 36.998, 37.245, 37.451, 37.577, 
37.615, 37.71, 37.764, 37.808, 37.841, 38.12, 38.978, 39.197, 
39.419, 39.517, 39.527, 39.587, 39.749, 39.779, 39.802, 39.93, 
40.188, 40.353, 40.383, 40.428, 40.564, 40.676, 40.687, 40.83, 
40.929, 41.01, 41.234, 41.405, 41.446, 41.716, 41.727, 42.078, 
42.33, 42.436, 42.723, 42.741, 42.913, 42.942, 43.034, 43.126, 
43.127, 43.299, 43.474, 43.49, 43.555, 43.823, 43.92, 44.009, 
44.069, 44.218, 44.681, 44.746, 44.755, 45.082, 45.117, 45.171, 
45.712, 45.779, 45.978, 46.04, 46.235, 46.287, 46.381, 46.485, 
46.685, 46.897, 47.154, 47.303, 47.349, 47.597, 47.706, 48.018, 
48.06, 48.136, 48.155, 48.577, 48.967, 49.009, 49.023, 49.13, 
49.14, 49.275, 49.359, 49.418, 49.694, 49.728, 49.803, 49.806, 
49.836, 49.871, 50.068, 50.086, 50.2, 50.587, 50.657, 50.692, 
50.822, 50.822, 50.936, 50.967, 51.034, 51.083, 51.084, 51.158, 
51.373, 51.623, 51.652, 51.833, 51.869, 52.08, 52.206, 52.443, 
52.519, 52.56, 52.604, 52.675, 52.728, 52.99, 53.12, 53.142, 
53.357, 53.463, 53.519, 53.916, 53.958, 54.109, 54.404, 54.474, 
54.561, 54.969, 54.987, 55.056, 55.144, 55.233, 55.246, 55.628, 
55.843, 56.022, 56.047, 56.241, 56.364, 56.701, 57.04, 57.206, 
57.233, 57.634, 57.796, 57.851, 57.971, 58.169, 58.218, 58.401, 
58.484, 58.858, 58.87, 59.004, 59.144, 59.33, 59.586, 59.628, 
59.675, 59.931, 59.955, 60.142, 60.228, 60.23, 60.281, 60.326, 
60.829, 60.872, 60.902, 61.103, 61.267, 61.498, 61.743, 61.943, 
62.079, 62.398, 62.433, 62.484, 63.194, 63.535, 63.872, 63.98, 
64.162, 64.242, 64.623, 64.869, 64.874, 64.973, 65.075, 65.17, 
65.193, 65.27, 65.272, 65.425, 65.504, 65.76, 65.786, 65.858, 
65.991, 66.139, 66.39, 66.442, 66.483, 66.607, 66.653, 66.872, 
67.274, 67.874, 67.905, 68.207, 68.406, 68.591, 68.646, 68.708, 
68.937, 69.064, 69.116, 69.215, 69.413, 69.506, 69.637, 69.717, 
69.771, 70.083, 70.096, 70.234, 70.441, 70.576, 70.603, 70.754, 
70.97, 71.018, 71.179, 71.264, 71.301, 71.304, 71.319, 71.464, 
71.495, 71.787, 72.447, 72.611, 72.894, 73.149, 73.227, 73.439, 
73.821, 73.948, 74.319, 74.479, 74.58, 74.694, 75.091, 75.113, 
75.319, 75.403, 75.916, 76.315, 76.559, 76.566, 76.667, 76.707, 
76.732, 76.769, 76.828, 76.843, 77.027, 77.064, 77.607, 77.681, 
77.76, 78.243, 78.262, 78.376, 78.72, 78.739, 78.812, 78.901, 
78.901, 79.232, 79.578, 79.741, 80.048, 80.079, 80.08, 80.245, 
80.388, 80.402, 80.404, 80.542, 81.288, 81.357, 81.455, 81.648, 
82.652, 82.672, 82.686, 83.092, 83.1, 83.318, 83.336, 83.36, 
83.569, 83.649, 83.675, 83.682, 83.822, 84.007, 84.014, 84.714, 
85.829, 85.906, 86.272, 86.323, 86.795, 87.092, 87.431, 87.63, 
87.898, 88.267, 88.44, 88.907, 88.948, 89.132, 89.178, 89.297, 
89.347, 89.448, 90.168, 90.326, 90.72, 90.879, 91.463, 91.595, 
91.825, 92.005, 92.173, 92.221, 92.356, 92.474, 92.972, 93.06, 
93.448, 93.484, 93.602, 93.66, 93.782, 93.969, 94.141, 94.309, 
94.789, 95.033, 95.127, 95.187, 95.207, 95.506, 95.948, 96.175, 
96.4, 97.238, 97.256, 97.741, 97.882, 98.035, 98.107, 98.215, 
98.23, 98.585, 98.988, 99.021, 99.432, 99.564, 99.692, 100.303, 
100.417, 100.763, 101.986, 102.66, 103.8, 104.474, 105.238, 105.353, 
105.459, 105.536, 105.654, 105.817, 106.096, 106.351, 106.46, 
106.726, 107.068, 107.182, 107.33, 107.63, 107.817, 107.975, 
109.148, 109.637, 109.716, 109.81, 112.025, 112.096, 112.735, 
113.019, 113.199, 113.513, 113.856, 113.959, 114.109, 114.573, 
114.595, 114.95, 115.12, 115.312, 116.659, 116.666, 116.722, 
116.944, 117.475, 118.151, 118.489, 118.676, 118.852, 118.958, 
120.502, 121.004, 121.294, 121.624, 122.098, 122.178, 122.438, 
122.889, 123.647, 123.805, 124.233, 125.181, 125.266, 126.595, 
126.671, 126.787, 128.209, 128.384, 128.743, 129.137, 129.653, 
129.703, 129.803, 130.168, 130.243, 130.538, 130.565, 131.358, 
131.587, 131.888, 132.201, 132.287, 132.305, 132.494, 132.842, 
133.182, 134.049, 134.378, 135.706, 136.275, 137.43, 137.532, 
137.763, 137.81, 138.102, 138.129, 138.164, 138.704, 139.153, 
139.681, 139.786, 140.216, 140.301, 140.436, 141.131, 141.207, 
141.512, 141.961, 142.386, 142.663, 142.831, 143.135, 143.643, 
144.102, 144.922, 145.102, 145.41, 145.435, 145.526, 145.788, 
145.833, 148.819, 149.299, 149.338, 149.492, 152.745, 153.903, 
153.999, 154.544, 155.61, 156.059, 156.181, 157.377, 159.989, 
160.486, 161.313, 161.937, 162.872, 163.361, 164.642, 165.107, 
167.107, 167.913, 168.174, 168.509, 168.984, 169.393, 170.868, 
170.987, 171.3, 171.582, 171.812, 172.09, 173.63, 175.291, 176.55, 
177.837, 179.94, 180.784, 181.63, 182.424, 182.632, 189.035, 
189.051, 190.28, 193.699, 193.876, 195.126, 198.048, 198.311, 
199.291, 201.772, 202.172, 202.478, 203.964, 206.532, 210.209, 
211.767, 212.361, 212.828, 216.334, 221.223, 224.637, 227.766, 
228.366, 230.538, 232.842, 234.736, 238.839, 240.615, 241.071, 
243.32, 249.356, 280.831, 300.903, 363.072, 420.303, 434.876, 
446.507, 448.495, 600.195, 607, 624.36, 690.056, 792.878, 827.447, 
844.563, 848.537, 866.2, 897.563, 955.887 ), n = 842 )

lognormal.inits <- list( mu = 0, log.tau = 0, y.new = 100 )

set.seed( 1 )

n.parallel.chains <- 2

lognormal.run.1 <- jags.model( file = "lognormal-model.txt", 
  data = income.data, inits = lognormal.inits, 
  n.chains = n.parallel.chains )

n.burnin <- 1000

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

n.monitor <- 50000

system.time(

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

)

n.monitor <- n.monitor * n.parallel.chains

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

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

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

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

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 ) ) )

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 ) ) )

plot( 1:n.monitor, theta.star, type = 'l', xlab = 'Iteration thetamber',
  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 ) ) )

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 = '' )

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 ) )

par( mfrow = c( 1, 1 ) )

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

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

c( mean( income.data$y ), sd( income.data$y ), 
  quantile( income.data$y, probs = c( 0.025, 0.975 ) ) )

system.time(

  print( lognormal.dic <- dic.samples( lognormal.run.1, n.iter = 50000, 
    type = 'pD' ) )

)