Guidelines for the stochastic production model in continuous time (SPiCT)

Main assumptions and input data for SPiCT

• Catch data should be representative of both landings and bycatch. It is also possible to use landings only, but then the interpretation of the results changes. If available, seasonal catches should be used as input. Catches are assumed to be taken over a time interval (e.g. years or quarters), thus the associated time vector in SPiCT inp$timeC should reflect the beginning of each catch interval (e.g. 2002.25 and 2002.75 for the second and fourth quarter catches, respectively). Additionally, the vector inp$dtc should reflect the length of each time interval (e.g. 1 for annual and 0.25 for quarterly catches, respectively).

• Stock size indices should be in terms of biomass (not numbers) and representative of the part of the stock vulnerable to the commercial fleets, the so called exploitable stock biomass (ESB). In many cases, the gear selectivity of the commercial and scientific fleets do not coincide and thus the stock size indices have to be corrected to exclude individuals that are not represented in the commercial fleets.

• Biomass indices are assumed to be snapshots at given points in time. Therefore, the timing of survey indices inp$timeI has to be given as decimal years reflecting the timing of the survey (e.g. 1995.5 for the middle of the year). The timing of the survey will be matched to the closest model time which is dependent on inp$dteuler (see below). Commercial CPUE indices should be associated with the midpoint of the interval of the corresponding catches, i.e. middle of the year if they are based on yearly aggregated catches and effort.

• If input data spans long periods and there is ecological evidence for potential changes in productivity, it is possible to fit productivity regime shifts or long-term gradual changes in productivity with SPiCT (Mildenberger et al. 2019).

• The Euler discretisation has to be finer than the data, i.e. for yearly data the Euler time step inp$dteuler has to be smaller than 1. A value of 1 changes the model to a discrete time production model with different model assumptions. The default value is 1/16 year, which is sufficient in most cases.

Checklist for the acceptance of a SPiCT assessment

ICES category 3 stocks can be managed using the official advice rule 3.1.1 and 3.1.2 in ICES (2017), which requires an accepted SPiCT assessment. Following points should be considered for the acceptance of a SPiCT assessment.

1. The assessment converged (fit$opt$convergence equals \(0\)).

2. All variance parameters of the model parameters are finite (all(is.finite(fit$sd)) should be TRUE).

3. No violation of model assumptions based on one-step-ahead residuals (bias, auto-correlation, normality). This means, that p-values are insignificant (\(> 0.05\)), indicated by green titles in the graphs of plotspict.diagnostic(fit). Slight violations of these assumptions do not necessarily invalidate model results.

4. Consistent patterns in the retrospective analysis (fit <- retro(fit)). This means that there is no tendency of consistent under- or overestimation of the relative fishing mortality (\(\frac{F}{F_{MSY}}\)) and relative biomass (\(\frac{B}{B_{MSY}}\)) in successive assessment. The retrospective trajectories of those two quantities should be inside the credible intervals of the base run.

5. Realistic production curve. The shape of the production curve should not be too skewed (\(\frac{B_{MSY}}{K}\) should be between 0.1 and 0.9). Low values of \(\frac{B_{MSY}}{K}\) allow for an infinite population growth rate (calc.bmsyk(fit)).

6. High assessment uncertainty can indicate a lack of contrast in the input data or violation of the ecological model assumptions. The main variance parameters (logsdb, logsdc, logsdi, logsdf) should not be unrealistically high. Credible intervals for \(\frac{B}{B_{MSY}}\) and \(\frac{F}{F_{MSY}}\) should not span more than 1 order of magnitude (calc.om(fit)).

7. Initial values do not influence the parameter estimates (fit <- check.ini(fit)). The estimates should be the same for all initial values (fit$check.ini$resmat). Runs which did not converge should not be considered in this regard.

Optional model and assessment tuning

Disclaimer: An uncertain assessment is okay; the model should not be tweaked and parameters should not be fixed unless there is sufficient information and evidence to make such assumptions. Rather than avoiding uncertainty, it should be accounted for by means of stochastic harvest control rules (ICES 2019).

• Increase iterations of optimisation. The special error code false convergence (8) could indicate that the optimisation exceeded the maximum number of iterations max.iter. inp$optimiser.control = list(iter.max = 1e5, eval.max = 1e5)

• If the catch time series is longer than the survey(s), shortening the catch time series to cover only the period where there is an available biomass index may help model convergence. inp = shorten.inp(inp, 2005, 2018)

• Adjust the prior for the parameter logn determining the shape of the production curve.

  • Alternative priors for logn could be based on meta studies, e.g. for all species pooled based on Thorson et al. (2012): inp$priors$logn <- c(log(1.478), 0.6, 1)

  • Tighter Schaefer prior for logn (or any other prior e.g. from meta studies): inp$priors$logn <- c(log(2), 0.5, 1)

  • Fixing n to resemble the Schaefer production model (or the meta study, alternatively): inp$ini$logn <- log(2); inp$phases$logn <- -1

Be aware that these modifications might cause smaller credible bands.

• Use a prior for the initial depletion level based on available information, FOR EXAMPLE:

  • If evidence or expert knowledge allows to infer that there was low or no exploitation before the beginning of the available data: initial depletion level could be assumed to be close to the carrying capacity (e.g. inp$priors$logbkfrac <- c(log(0.8),0.5,1))

  • If evidence or expert knowledge allows to infer that there was high exploitation before the beginning of the available data: initial depletion level could be assumed to be a certain fraction of the carrying capacity (e.g. inp$priors$logbkfrac <- c(log(0.2),0.5,1))

• If information on the level of uncertainty in the biomass index and/or the uncertainty of the catch is available prior distribution could be used for the observation error term of the indices logsdi and catches logsdc, respectively. This requires to remove the priors for the ratios of process to observation errors logalpha and logbeta.

Management based on an accepted SPiCT assessment

SPiCT estimates MSY based reference levels, which can be used to calculate quantities relevant for fisheries management, such as the Total Allowable Catch (TAC). The TAC can be calculated based on fishing at \(F_{MSY}\) or based more elaborate harvest control rules (HCR) which include additional components, e.g. a linearly reduction in F if the biomass is below a certain reference level (e.g., \(\text{MSY} B_{\text{trigger}}\)) as the hockey-stick MSY rule recommended by ICES (2017).

Mildenberger et al. (2022) demonstrated the importance of the quantification and consideration of assessment uncertainty into the calculation of the catch advice. The authors introduced a modification of the hockey-stick MSY rule that reduces catch advice as a function of uncertainty by using a percentile lower than 50 (i.e., median) for relative biomass \(\frac{B}{B_{\text{MSY}}}\) and catch \(C_{y+1}\) and higher than 50 for relative fishing mortality \(\frac{F}{F_{\text{MSY}}}\).

Based on this work, ICES currently recommends calculating TAC based on an MSY hockey-stick rule with a biomass trigger reference point at \(\text{MSY} B_{\text{trigger}} = 0.5 B_{\text{MSY}}\), a biomass limit reference point at \(B_{\text{lim}} = 0.3 B_{\text{MSY}}\), and the 35th percentile of the predicted catch distribution (ICES 2025). Under this rule, the target fishing mortality \(F^{\tau}\) equals \(F_{\text{MSY}}\) when \(B\) is above \(\text{MSY}B_{\text{trigger}}\); when \(B\) is between 0 and \(\text{MSY}B_{\text{trigger}}\), \(F^{\tau}\) is reduced linearly with \(B\); and when \(B\) falls below \(B_{\text{lim}}\), \(F^{\tau}\) is set to zero. This can be expressed as

\[\begin{equation*} F^{\tau}_y = \begin{cases} F_{\text{MSY}} \min\!\left(1, \dfrac{B_y}{\text{MSY}B_{\text{trigger}}}\right), & \text{if } B_y \ge B_{\text{lim}},\\[6pt] 0, & \text{if } B_y < B_{\text{lim}}. \end{cases} \end{equation*}\]

Given the target fishing mortality, the advised catch \(A_y\) for year \(y\) is taken as the 35th percentile of the predicted catch distribution:

\[\begin{equation*} A_y = \Phi^{-1}_{C_{\text{pred},y} \mid F^{\tau}_y}(0.35), \end{equation*}\]

where \(\Phi^{-1}\) denotes the inverse cumulative distribution function, so that \(\Phi^{-1}_{C_{\text{pred},y} \mid F^{\tau}_y}(0.35)\) is the catch corresponding to the 35th percentile of the predicted catch distribution.

This HCR is one of the 8 default management scenarios of the manage() function in spict but can als be run exclusively with:

fit <- manage(fit, scenarios = "ices").

This function is actually only a wrapper function and calls internally the add.man.scenario() function that makes the three key components of the ICES (2025) HCR explicit:

fit <- add.man.scenario(fit, "ICES (2025)", fractiles = list(catch = 0.35), breakpointB = 0.5, limitB = 0.3).

Another way to calculate the TAC based on the ICES (2025) HCR directly is using the get.TAC() function:

get.TAC(fit, fractiles = list(catch = 0.35), breakpointB = 0.5, limitB = 0.3).

It is important to note that informative priors (i.e. priors with smaller standard deviations) can influence parameter estimates and credible intervals in SPiCT assessments, and thereby affect the resulting TAC.