Surplus production model

Dynamics equations

The surplus production model uses the Fletcher (1978) formulation. The biomass \(B_t\) in year \(t\) is \[B_t = B_{t-1} + P_{t-1} - C_{t-1},\] where \(C_t\) is the observed catch and \(P_t\) is the surplus production given by: \[P_t = \gamma \times MSY \times \left(\dfrac{B_t}{K}-\left[\dfrac{B_t}{K}\right]^n\right), \] where \(K\) is the carrying capacity, \(MSY\) is the maximum sustainable yield, and \(n\) is the parameter that controls shape of the production curve, and \(\gamma\) is \[\gamma = \dfrac{1}{n-1}n^{n/(n-1)}.\]

The fishing mortality \(F_t\) is conditioned on the observed catch \(C_t\). If fishing occurs as a pulse fishery, then \[\hat{F}_t = \dfrac{C_t}{\hat{B}_t},\] where \(F_t\) is effectively an exploitation rate. Alternatively, smaller time steps are used in the model to approximate continuous production and fishing. Given the biomass in the start of the year and assuming a constant fishing mortality over the time steps within a year, the fishing mortality that produces the observed annual catch is solved iteratively. By default, four sub-annual time steps are used to solve for \(F_t\).

Estimated and derived parameters

The model is configured with \(F_{MSY}\) and \(MSY\) as leading parameters to be estimated.

The biomass \(B_{MSY}\) at \(MSY\) is \[B_{MSY} = \dfrac{MSY}{F_{MSY}},\] the carrying capacity \(K\) is \[K = n^{1/(n-1)} \times B_{MSY} ,\] and the intrinsic rate of population increase \(r\) is \[ r = n \times F_{MSY}.\]

Schaefer model

The production parameter \(n\) is typically fixed, and the Schaefer model with a symmetric productive curve (\(B_{MSY}/K = 0.5\)) is produced when \(n = 2\) (default).

Fox model (SP_Fox)

The Fox model is the limiting case of the Fletcher parameterization as \(n \rightarrow 1\), where

\[K = e \times B_{MSY}\] \[r = F_{MSY}\] \[P_t = -e \times MSY \times \dfrac{B_t}{K} \times \log\left(\dfrac{B_t}{K}\right) \] where \(e\) is Euler’s number.

State-space version (SP_SS)

In the state-state version, annual biomass deviates are estimated as random effects. Similar to Meyer and Millar (1999), the biomass \(B_t\) in year \(t\) is \[B_t = (B_{t-1} + P_{t-1} - C_{t-1})\exp(\delta_t - 0.5 \tau^2),\] where \(\delta_t \sim N(0, \tau^2)\) are biomass deviations in lognormal space and \(\tau\) is the standard deviation of the biomass deviations.

Likelihood

By conditioning the model on observed catch, the model is fitted to the indices of abundance \(I_{s,t}\). The likelihood of the observed indices, using assuming a lognormal distribution, is \[L^I = \sum_s \lambda_s \sum_t \left(-\log(\sigma_{s,t}) - 0.5\left[\dfrac{\log(I_{s,t}) - \log(\hat{q}_s \hat{B}_t)}{\sigma_{s,t}}\right]^2\right)\] where \(lamba_s\) is the weighting coefficient of index \(s\), \(q_s\) is the cathabilitiy coefficient, and the circumflex denotes an estimate.

Log-biomass deviations \(\hat{\delta}_t\) are typically estimated as penalized parameters towards the likelihood, with the penalty: \[L^{\delta} = \sum_t \left(-\log(\tau) - 0.5\left[\dfrac{\delta_t}{\tau}\right]^2\right).\]

Prior for r

To generate the prior for the intrinsic rate of increase, natural mortality \(M\) and steepness \(h\) are sampled from a distribution. Natural mortality is modelled to be age-invariant and is sampled from a lognormal distribution.

With either a Beverton-Holt or Ricker stock-recruit relationship, steepness is sampled from a transformed beta or transformed lognormal distribution, respectively.

For each pair of sampled M and h values, the corresponding value of \(r\) is obtained by solving a modified Euler-Lotka equation: \[\alpha \sum_a l_a f_a \exp(-r \times a) = 1.\]

Equation 15 is modified to include the \(\alpha\) term from the stock-recruit relationship (Stanley et al. 2009). In this way, the recruits-per-spawner at low stock sizes, i.e., as spawning biomass approaches zero, is considered for calculating \(r\).

The numbers-per-recruit at age \(a\) is \[ l_a = \begin{cases} 1 & a = 1\\ l_{a-1} \exp(-M_{a-1}) & a = 2, \ldots, A-1\\ \dfrac{l_{a-1} \exp(-M_{a-1})}{1 - \exp(-M_a)} & a = A \\ \end{cases}. \]

Fecundity at age \(f_a\) is \[m_a = w_a\left[1 + \exp\left(-\log(19) \dfrac{a - \tilde{a}_{50}}{\tilde{a}_{95} - \tilde{a}_{50}}\right)\right]^{-1}, \] where \(\tilde{a}_{50}\) and \(\tilde{a}_{95}\) are the ages of 50% and 95% maturity, respectively.

Weight-at-age \(w_a\) is \[w_a = W_{\infty}(1 - \exp[K\{a-a_0\}])^b.\]

The recruits per spawner at the origin of the stock-recruit relationship \(\alpha\) is \[\alpha = \dfrac{4h}{(1-h)\phi_0}, \] or \[\alpha = \dfrac{(5h)^{1.25}}{\phi_0},\] for a Beverton-Holt and Ricker stock-recruit relationship, respectively. Unfished recruits-per-spawner \(\phi_0\) is \[\phi_0 = \sum_a l_a f_a.\]

A normal distribution is assumed for the prior (a penalty to the likelihood) with the mean (\(\mu_r\)) and standard deviation (\(\sigma_r\)) calculated from the values of \(r\) calculated using the above procedure:

\[L^r = - 0.5\left(\dfrac{\hat{r} - \mu_r}{\sigma_r}\right)^2,\]

References

Fletcher, R.I. 1978. On the restructuring of the Pella-Tomlinson system. Fishery Bulletin 76:515-521.

Meyer, R., and Millar, R.B. 1999. BUGS in Bayesian stock assessments. Canadian Journal of Fisheries and Aquatic Science 56:1078-1086.

Stanley, R.D., M. McAllister, P. Starr and N. Olsen. 2009. Stock assessment for bocaccio (Sebastes paucispinis) in British Columbia waters. DFO Can. Sci. Advis. Sec. Res. Doc. 2009/055. xiv + 200 p.