### 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.