State-space modeling of stock-recruitment data

WinBUGS code from Su and Peterman (2012) "Performance of a Bayesian state-space model of semelparous species for stock-recruitment data subject to measurement error "

WinBUGS program for the state-space Ricker model

model

{

    # Observation equation

    for (t in 1:T)

    {

        E[t] ~ dlnorm(mu.logE[t], prec.E); mu.logE[t] <- log(S[t]);                        

    }

    # Prior distributions for the state variables

    for (t in 1:k)

    {

        S[t] <- S1[t]; S1[t] ~ dlnorm(0.0, 1.0E-3)   

    }

    for (t in 1:(T - k)) {

        mu[t] <-  b[1] + log(S[t]) - b[2] * S[t];

        R[t] ~ dlnorm(mu[t], prec.R)I(C[t], 1.0E+9); # R and C represent Rt+k and Ct+k, respectively

        S[t + k] <- R[t] - C[t];  

    }

    # Prior distributions for model parameters

    # a bivariate normal prior for a_r and beta

    b[1:2] ~ dmnorm(mu.b[], prec[,])

    a_r <- b[1]; alpha <- exp(b[1]);

    beta <- b[2]

    # Priors for standard deviations

    tau ~ dunif(0, 1.0E+9); tau2 <- tau * tau; prec.R <- 1/tau2

    sig ~ dunif(0, 1.0E+9); sig2 <- sig * sig; prec.E <- 1/sig2

    # Management-oriented quantities

    SMSY <- hMSY/beta

    hMSY <- 0.5*a_r - 0.07*a_r *a_r

}

Note: This program uses a bivariate normal prior distribution for a_r and b. This instructs WinBUGS to update a_r and b in a block as was done for our Bayesian state-space model (SSR). However, a_r and b are not restricted to be positive as in the case for a_r and b in SSR. Nevertheless, the results from the WinBUGS program can be compared to those from SSR by discarding those draws with non-positive a_r and b.

Here is an example data set generated from the conditions: k = 2, alpha = 4, beta = 1, sig2 = 0.01, tau2 = 0.003, rho = 0.25, l = 0, T = 28, and the “variable” harvest situation:

Data: list(mu.b = c(1, 1), prec= structure(.Data= c(1.0E-03, 0, 0, 1.0E-06), .Dim=c(2, 2)), T = 28, k = 2, E = c(1.062, 1.179, 1.385, 1.213, 1.176, 1.385, 1.141, 0.944, 1.087, 1.119, 0.966, 0.918, 0.793, 0.775, 0.645, 0.661, 0.548, 0.703, 0.807, 0.792, 0.9, 1.007, 1.222, 1.269, 1.001, 1.108, 1.032, 1.099), C = c(0.165, 0.23, 0.272, 0.34, 0.385, 0.409, 0.508, 0.596, 0.609, 0.669, 0.674, 0.779, 0.708, 0.864, 0.611, 0.658, 0.631, 0.583, 0.548, 0.534, 0.556, 0.513, 0.332, 0.319, 0.235, 0.2))

Inits:

list(b=c(6, 2), tau = 2, sig = 3, R=c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), S1 =c(1, 1))

 list(b=c(3, 0.5), tau = 1, sig = 1, R=c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2), S1 =c(2, 2))