#AGE-STRUCTURED SEQUENTIAL POPULATION ANALYSIS
#This code is a modified version of that used in 
# Millar and Meyer (2000) "Bayesian state-space modeling of 
# age structured data: fitting a model is just the beginning. CJFAS 57:43-50.
#
#The data are from Myers and Cadigan (1995a)

model cod;
{
######Vague priors on first row and column#######
for(y in 1:(Y+1)) {LogRecruits[y]~dnorm(10,.01);
                  LogN[1,y]<-LogRecruits[y]; }
for(a in 2:(A+1)){LogInitials[a-1]~dnorm(10,.01);
                  LogN[a,1]<-LogInitials[a-1]; }

########Distribution for data, I[a,y]######
for (a in 1:A) {
  for (y in 1:Y) { MeanLogI[a,y]<-Logq[a]+(LogN[a,y]+11*LogN[a+1,y+1])/12;
                   Prec[a,y]<-PrecLogI/a;
                   I[a,y]~dlnorm(MeanLogI[a,y],Prec[a,y]); } }

#####Prior on mortalities, hierarchical######NB: log(0.2)=-1.61
for(a in 1:A){ MeanLogAgeM[a]~dnorm(-1.61,1);}
for(y in 1:Y){ MeanLogDeltaM[y]~dnorm(0,1);}
PrecLogM~dgamma(1,1);
for(a in 1:A) {
  for(y in 1:Y) { MeanLogM[a,y]<-MeanLogAgeM[a]+MeanLogDeltaM[y];
                  M[a,y]~dlnorm(MeanLogM[a,y],PrecLogM)I(0,10.0); } }

#####Cohort equation#####
for(a in 2:(A+1)) {
  for(y in 2:(Y+1)) {
    LogN[a,y]<-log(max(exp(LogN[a-1,y-1])*exp(-M[a-1,y-1])-C[a-1,y-1]*exp(-0.5*M[a-1,y-1]),10)); } }


######Prior on q's######
MeanLogq<--10; PrecLogq<-1;
for (a in 1:A){ Logq[a]~dnorm(MeanLogq,PrecLogq); }

######Prior on q's, quadratic######
#gamma~dnorm(-8,1)I(-10,-6); 
#alpha~dunif(0.0001,0.02); 
#beta~dnorm(5,0.25)I(1,11);
#for (a in 1:A){ MeanLogq[a]<-gamma-alpha*(a-beta)*(a-beta);
#                Logq[a]~dnorm(MeanLogq[a],16); }

######Prior on precision of log(I)######
PrecLogI~dgamma(10,1);

######Vector of 1983 numbers for monitoring#####
for(a in 1:A) {N1983s[a]<-exp(LogN[a,6])/1000; }
}

list(
A=10, Y=16,
I = structure(.Data = c(5.39, 1.94, 2.48, 5.12, 5.87, 12.22, 10.79, 7.27, 
 4.77, 2.04, 3.93, 8.98, 10.93, 3.35, 1.78, 0.6, 11.51, 11.78, 3.83, 
 2.74, 5.92, 10.62, 15.23, 12.35, 20.7, 4.03, 3.2, 8.3, 12.95, 13.97, 
 2.3, 0.83, 13.95, 16.79, 13.23, 3.26, 3.83, 10.83, 11.34, 10.01, 31.29, 
 13.23, 5.29, 6.2, 8.61, 9, 2.72, 0.34, 5.51, 10.53, 13.31, 9.67, 2.79, 
 3.87, 9.59, 7.28, 21.28, 11.61, 10.57, 6.52, 5.64, 3.31, 1.42, 0.22, 
 1.62, 2.27, 4.99, 8.79, 5.82, 2.43, 2.3, 4.24, 10.14, 4.38, 10.13, 8.23,
 3.9, 1.1, 0.35, 0.04, 0.63, 0.92, 1.19, 3.66, 5.31, 5.33, 1.37, 0.92, 
 5.26, 2.67, 2.58, 4.84, 3.98, 0.5, 0.04, 0.01, 0.47, 0.31, 0.37, 0.74, 
 2.59, 2.93, 2.09, 0.78, 1.37, 1.38, 1.55, 1.62, 1.68, 0.35, 0.02, 0.001, 
 0.33, 0.26, 0.23, 0.23, 0.57, 1.42, 1.3, 0.67, 0.58, 0.34, 0.79, 0.98, 
 0.55, 0.16, 0.01, 0.001, 0.12, 0.19, 0.11, 0.1, 0.16, 0.36, 0.54, 0.41, 
 0.68, 0.17, 0.15, 0.43, 0.23, 0.04, 0.001, 0.001, 0.09, 0.06, 0.16, 0.11,
 0.09, 0.14, 0.28, 0.15, 0.42, 0.19, 0.11, 0.16, 0.12, 0.02, 0.01, 0.001), 
 .Dim = c(10,16)), 
C = structure(.Data = c(1323, 1152, 2554, 2185, 1702, 2585, 
 782, 650, 831, 2329, 2779, 1696, 7693, 3111, 430, 940, 17556, 12361, 
 12025, 7172, 31286, 13616, 14871, 14824, 15219, 9217, 14651, 17639, 
 40557, 31654, 3860, 4993, 39206, 37439, 28814, 13191, 19003, 42602, 
 31760, 36614, 44168, 32340, 20184, 21150, 36410, 53805, 14535, 3343, 
 20319, 29202, 30016, 24800, 14297, 19028, 38624, 33922, 45869, 49061, 
 47917, 25212, 22695, 29553, 12211, 1940, 7711, 10982, 18017, 22014, 
 25435, 12044, 12503, 28006, 26025, 28469, 45725, 38708, 16390, 9064, 
 4526, 700, 3078, 3460, 4830, 11848, 16930, 14701, 7246, 7050, 14722, 
 19505, 18608, 28499, 17940, 6164, 1372, 147, 1530, 1300, 1217, 3175, 
 11936, 8934, 8910, 3836, 3104, 5818, 9026, 8696, 9156, 4745, 376, 21, 
 1083, 757, 520, 779, 1923, 6341, 4227, 5162, 2000, 1346, 4337, 3640, 
 2865, 1696, 199, 0.1, 437, 560, 232, 309, 338, 1018, 2536, 2905, 1977, 
 676, 774, 1695, 1084, 641, 104, 0.1, 324, 299, 285, 320, 246, 338, 597, 
 1935, 1675, 1264, 788, 816, 581, 338, 27, 0.1), .Dim = c(10,16))
)


list(
LogRecruits=c(13,12,12,13,13,13,13,13,12,12,12,12,12,11,11,11,11),
LogInitials=c(13,12,12,11,10,9,9,8,8,7),
PrecLogI= 10)