#State-space surplus production model for South Albacore tuna
#############################################################

model Tuna;
{
######Distribution of biomass index#####
for (i in 1:N) { 
  Imean[i] <- log(Q*P[i])
  I[i] ~ dlnorm(Imean[i],itau2) }

######Distribution of P=B/K series######
Pmean[1] <- 0;
P[1] ~ dlnorm(Pmean[1],isigma2)

for (i in 2:N) {
    Pmean[i] <- log(max(P[i-1] + r*P[i-1]*(1-P[i-1]) - C[i-1]/K,0.01));
    P[i] ~ dlnorm(Pmean[i],isigma2) }

###Lognormal prior on K giving 10th and 90th percentiles of 80 and 300 (1000's t)###
eta1 <- 5.042905; v1 <- 0.26593; prec1 <- 1/v1; #(prec1=3.7604);
K ~ dlnorm(eta1,prec1)  #Baseline: Mean of K is 176.95

###Lognormal prior on r (obtained from meta analysis)###
eta2 <- -1.38; v2 <- 0.26010; prec2 <- 1/v2; #(prec2=3.8447);
r~dlnorm(eta2,prec2)  #Baseline: Mean of r is 0.28652

#######Prior on precision of process error#####
alpha1 <- 3.7855;  beta1 <- 0.010223;
isigma2 ~ dgamma(alpha1,beta1); 
sigma2<-1/isigma2;  #Baseline: Mean of sigma2 is 0.0036701

#####Prior on precision of observation error#####
alpha2 <- 1.7086; beta2 <- 0.0086139;
itau2 ~ dgamma(alpha2,beta2); 
tau2<-1/itau2;

#####Reference prior on Q#####
Q ~ dgamma(0.001,0.001)
q <- Q/K;

#####Define maximum sustainable yield#####
MSY <- r*K/4;

#######Logical nodes used in sensitivity analysis#######
###Parameter of interest
D[1] <- MSY;
###Nodes used for sensitivity of MSY to lognormal priors on K and r
D[2] <- log(K)/v1; D[3] <- log(r)/v2
###Nodes used for sensitivity of MSY to inverse gamma priors on variances
D[4] <- isigma2*(-alpha1); D[5] <- itau2*(-alpha2)
###Nodes used for posterior sensitivity of K and r to their priors
D[6] <- K/exp(eta1+0.5*v1); D[7] <- r/exp(eta2+0.5*v2)
###Nodes used for posterior sensitivity of variances to their priors
D[8]<-sigma2*(alpha1-1)/alpha1; D[9]<-tau2*(alpha2-1)/alpha2
}

#####DATA#####
list( N=23,
C=c(15.9, 25.7, 28.5, 23.7, 25.0, 33.3, 28.2, 19.7, 17.5, 19.3, 21.6, 23.1,
              22.5, 22.5, 23.6, 29.1, 14.4, 13.2, 28.4, 34.6, 37.5, 25.9, 25.3),
I=c(61.89, 78.98, 55.59, 44.61, 56.89, 38.27, 33.84, 36.13, 41.95, 36.63, 36.33, 38.82,
      34.32, 37.64, 34.01, 32.16, 26.88, 36.61, 30.07, 30.75, 23.36, 22.36, 21.91)
)
#####INITS#####
list(P=c(0.99,0.98,0.96,0.94,0.92,0.90,0.88,0.86,0.84,0.82,0.80,0.78,0.76,0.74,
0.72,0.70,0.68,0.66,0.64,0.62,0.60,0.58,0.56),r=0.1, K=200 Q=50,isigma2=100,itau2=100)