#Douglas P. Swain, Ian D. Jonsen, James E. Simon, and Ransom A. Myers. 2009. Assessing threats to stage-structured state–space models: mortality trends in skate populations. Ecological Applications 19:1347–1364.
#Appendix D. Example model code.
#1. Z-model for the eastern Scotian Shelf, allowing decadal variation in Z

model;
{
for (s in 1:3) {
      q[s] ~ dlnorm(0, 100)  # Informative Priors equivalent to q=1 (1/trawlable_units) with sd=0.1
      for (d in 1:3) {
            z[s,d] ~ dunif(0, 3)    # Priors for instantaneous mortality rates with decadal variation
            } 
      sdpro[s] ~ dunif(0, 5) # Process error variance
      isigma2[s] <- pow(sdpro[s], -2)
      }

# Priors for observation error variance
sdobs[1]~dunif(0.497,2)
sdobs[2]~dunif(0.586,2)
sdobs[3]~dunif(0.513,2)

for(i in 1:3){
    itau2[i] <- pow(sdobs[i], -2)
      }

# Priors (informative) for transistion probabilities
theta[1] ~ dbeta(20, 80)
theta[2] ~ dbeta(25, 75)

# Prior for recruitment rate
recrate ~ dlnorm(1.178655, 4)I(0.01,1000)

# Informative priors for inital abundance (based on trawlable abundance in 1000s)
# sd=20
# Nad[1] is adult abundance in 1965
N69[1] ~ dlnorm(0, 6.25)
N69[2] ~ dlnorm(-1.203973, 6.25)

for (i in 1:5) {Nad[i] ~ dlnorm(1.252763, 11.11111)}
# Process equation
for (s in 1:3) {
      for (t in 1:10 ) {
            st[s,t] <- exp(-z[s,1])
            st[s,t+10] <- exp(-z[s,2])
            st[s,t+20] <- exp(-z[s,3])
            }
      for(t in 31:35){
            st[s,t] <- exp(-z[s,3])
            }
     }

nmod[1,1] <- (Nad[1] * (recrate / 2) + N69[1] * (1 - theta[1]))* st[1,1]
nmod[2,1] <- (N69[1] * theta[1] + N69[2] * (1 - theta[2]))* st[2,1]
nmod[3,1] <- (N69[2] * theta[2] + Nad[5])* st[3,1]

for (s in 1:3) {
      logn[s,1] <- log(nmod[s,1])
      N[s,1] ~ dlnorm(logn[s,1], isigma2[s])
      }

nmod[1,2] <- (Nad[2] * (recrate / 2) + N[1,1]  * (1 - theta[1]))* st[1,2]
nmod[2,2] <- (N[1,1]  * theta[1] + N[2,1] * (1 - theta[2]))* st[2,2]
nmod[3,2] <- (N[2,1]  * theta[2] + N[3,1])* st[3,2]

for (s in 1:3) {
      logn[s,2] <- log(nmod[s,2])
      N[s,2] ~ dlnorm(logn[s,2], isigma2[s])
      }

for (t in 3:5) {
      nmod[1,t] <- (Nad[t] * (recrate / 2) + N[1,t-1]  * (1 - theta[1]))* st[1,t]
      nmod[2,t] <- (N[1,t-1]  * theta[1] + N[2,t-1]  * (1 - theta[2]))* st[2,t]
      nmod[3,t] <- (N[2,t-1]  * theta[2] + N[3,t-1])* st[3,t]

      for (s in 1:3) {
            logn[s,t] <- log(nmod[s,t])
            N[s,t] ~ dlnorm(logn[s,t], isigma2[s])
           }
      }

for (t in 6:35 ) {
      nmod[1,t] <- (N[3,t-5] * (recrate / 2) + N[1,t-1]  * (1 - theta[1]))* st[1,t]
      nmod[2,t] <- (N[1,t-1]  * theta[1] + N[2,t-1]  * (1 - theta[2]))* st[2,t]
      nmod[3,t] <- (N[2,t-1]  * theta[2] + N[3,t-1])* st[3,t]

      for (s in 1:3) {
            logn[s,t] <- log(nmod[s,t])
            N[s,t] ~ dlnorm(logn[s,t], isigma2[s])
            }
      }

# Observation equation
for (t in 1:35 ) {
            for (s in 1:3) {
            Im[s,t] <- q[s] * N[s,t]
            logIm[s,t] <- log(Im[s,t])
            Iobs[s,t] ~ dlnorm(logIm[s,t], itau2[s])
            #resid[s,t] <- Iobs[s,t] - Im[s,t]
            }
      }
} # End model