Stock Recruitment Relationship

In general, the relationship between spawning stock \(B_S\) and recruitment \(R\) is highly variable owing to intrinsic variability in factors governing early life history survival and to measurement error in the estimates of recruitment and the spawning biomass that generated it. The stock-recruitment relationship ultimately defines the sustainable yield 8 curve and its expected variability assuming that the stochastic processes of growth, maturation, and natural mortality are density-independent and stationary throughout the time horizon. Quinn and Deriso (1999) provide a useful discussion of stock-recruitment models, renewal processes, and sustainable yield. Note that the assumed stock-recruitment relationship does not affect the initial population abundance at the beginning of the time horizon (see Initial Population Abundance).

A total of twenty one stochastic recruitment models are available for population projection in the AGEPRO software. Thirteen of the recruitment models are functionally dependent on \(B_S\) while eight do not depend on spawning biomass. Five of the recruitment models have time-dependent parameters, twelve are time-invariant, and four may include time as a predictor, or not. The user is responsible for the choice and parameterization of the recruitment models. A description of each of the recruitment models follows. Important: note that the absolute units for recruitment \(R\) are numbers of age-1 fish, while for spawning biomass \(B_S\), the absolute units are kilograms of spawning biomass in each of the recruitment models below.

Model 1. Markov Matrix

A Markov matrix approach to modeling recruitment may be useful when there is uncertainty about the functional form of the stock-recruitment relationship. A Markov matrix contains transition probabilities that define the probability of obtaining a given level of recruitment given that \(B_S\) was within a defined interval range. In particular, the distribution of recruitment is assumed to follow a multinomial distribution conditioned on the spawning biomass interval or spawning state of the stock. The Markov matrix model depends on spawning biomass and is time-invariant.

An empirical approach to estimate a Markov matrix uses stock-recruitment data to determine the parameters of a multinomial distribution for each spawning biomass state. In this case, matrix elements can be empirically determined by counting the number of times that a recruitment observation interval lies within a given spawning biomass state, defined by an interval of spawning biomass, and normalizing over all spawning states. To do this, assume that there are \(K\) recruitment values and \(J\) spawning biomass states. The spawning biomass states are defined by disjoint intervals on the spawning biomass axis

\[ I_1 = {\bigl[0,B_{s,1}\bigr)}\ and\ for\ j=1, ..., J-2\ \ I_j = {\bigl[B_{S,j-1},B_{S,J}\bigr)}\ and\ I_j = {\bigl[B_{S,J-1}, \infty\bigr)} \tag{14}\label{eq:14} \]

where the spawning biomass values \(B_{S,j}\) are the input endpoints of the disjoint intervals of categories of spawning biomass (e.g., high, medium, low). Note that the spawning biomass intervals are defined by the cut points \(B_{S,j}\). The conditional probability of realizing the \(k^{th}\) recruitment value given that observed spawning biomass \(B_{S,Observed}\) is in the \(j^{th}\) interval is \(P{j,k}\). Here \(P_{j,k}\) is the element in the \(j^{th}\) row and \(k^{th}\) column of the Markov matrix \(\underline{P}={\bigl[ P_{j,k}\bigr]}\) of conditional recruitment probabilities with elements

\[ P_{j,k} = Pr{\bigl(R_k | B_{S,Observed} \in I_j\bigr)} \tag{15}\label{eq:15} \]

These individual conditional probabilities can be estimated by the computing the number of points in the stock recruitment data set that fall within a selected recruitment \([O_{k-1},O_k]\) range conditioned on the spawning biomass interval \(I_j\). If \(x_{j,k}\) represents the number of stock-recruitment observations in cell \(I_j \times O_k\) and there is at least one observation in spawning state \(j\), then the empirical maximum likelihood estimate of \(P_{j,k}\) is

\[ \Pr(r=O_k\ |\ B_S \in I_j)= \frac{x_{j,k}}{ \sum\limits_{k}x_{j,k} } \tag{16}\label{eq:16} \] Here \(P_{j,k} \ge 0\) and \(\sum\limits_{k=1}P_{j,k}=1\)

Up to 25 recruitment values and up to 10 spawning biomass states can be used in the Markov matrix model. For each spawning biomass interval, the user needs to specify the conditional probabilities of realizing the expected recruitment level, e.g., the \(P_{j,k}\). Given the conditional probabilities \(P_{j,k}\), simulated values of \(\hat{R}\) are generated by randomly sampling the conditional distribution \(\hat{R}(t) = Pr(R=O_k\ |\ B_S(t) \in I_j)\) through time.

Model 2. Empirical Recruits Per Spawning Biomass Distribution

For some stocks, the distribution of recruits per spawner may be independent of the number of spawners over the range of observed data. The recruitment per spawning biomass \((R/B_s)\) model randomly generates recruitment under the assumption that the distribution of the \(R/B_s\) ratio is stationary and independent of stock size. The empirical recruits per spawning biomass distribution model depends on spawning biomass and is time-invariant.

To describe this nonparametric approach, let \(S_t\) be the \(R/B_s\) ratio for the tth stock recruitment data point assuming age-1 recruitment

\[ S-t = \frac{R(t)}{B_s(t-1)} \tag{17}\label{eq:17} \]

and let \(R_S\) be the thS element in the ordered set of \(S_t\). The empirical probability density function for \(R_S\), denoted as \(g(R_S)\), assigns a probability of \(1/T\) to each value of \(R/B_S\) where \(T\) is the number of stock-recruitment data points. Let \(G(R_S)\) denote the cumulative distribution function (cdf) for \(R_S\). Set the values of \(G\) at the minimum and maximum observed \(R_S\) to be \(G(R_{min}) = 0\) and \(G(R_{max}) = 0\) so that the cdf of \(R_S\) can be written as

\[ G(R_S) = \frac{s-1}{T-1} \tag{18}\label{eq:18} \]

Random values of \(R_S\) can be generated by applying the probability integral transform to the empirical cdf. To do this, let \(U\) be a uniformly distributed random variable on the interval [0,1]. The value of \(R_S\) corresponding to a randomly chosen value of \(U\) is determined by applying the inverse function of the cdf \(G(R_S)\). In particular, if \(U\) is an integer multiple of \(1/(T-1)\) so that \(U=s/(T-1)\) then \(R_S = G^{-1}(U)\). Otherwise \(R_S\) can be obtained by linear interpolation when \(U\) is not a multiple of \(1/(T-1)\).

In particular, if \((s-1)/(T-1)<U<s/(T-1)\), then \(\hat{R_S}\) is interpolated between \(R_S\) and \(R_{S+1}\) as

\[ U = \left(\frac{\dfrac{s}{T-1} - \dfrac{s-1}{T-1}}{R_{s+1}-R_S} \right) (\widehat{R_S} - R_S) + \frac{s-1}{T-1} \tag{19}\label{eq:19} \]

Solving for \(\widehat{R_S}\) as a function of \(U\) yields

\[ \widehat{R_S} - (T - 1)(R_{S+1} - R_S)\left(U - \frac{s-1}{T-1}\right) + R_{S}\ \tag{20}\label{eq:20} \]

Where the interpolation index s is determined as the greatest integer in \(1+U(T-1)\). Given the value of \(\widehat{R_S}\), recruitment is generated as

\[ R(t) = N_1(t) = B_S(t-1) \cdot \widehat{R_S} \tag{21}\label{eq:21} \]

The AGEPRO program can generate stochastic recruitments under model 2 based on thousands of input stock-recruitment data points (i.e., the stock-recruitment data array size is defined as a long int variable in the C language and is user specified with the input variable MaxRecObs, see Table 3, keyword RECRUIT) conditioned on available computer memory resources.

Model 3. Empirical Recruitment Distribution

Another simple model for generating recruitment is to draw randomly from the observed set of recruitments \(\underline{R}= \bigl\{ R(1), R(2), ..., R(T) \bigr\}\) .This may be a useful approach when the recruitment has randomly fluctuated about its mean and appears to beindependent of spawning biomass over the observed range of data. In this case, the recruitment distribution may be modeled as a multinomial random variable where the probability of randomly choosing a particular recruitment is \(1 / T\) given \(T\) observed recruitments. The empirical recruitment distribution model does not depend on spawning biomass and is time-invariant.

In this model, realized recruitment \(\hat{R}\) is simulated from the empirical recruitment distribution as

\[ Pr\left(\widehat{R}=R(t)\right) = \frac{1}{T}, for\ t \in \{1,2,...,T\} \tag{22}\label{eq:22} \] The empirical recruitment distribution approach is nonparametric and assumes that future recruitment is totally independent of spawning stock biomass. When current levels of \(B_S\) are near the midrange of historical values this assumption seems reasonable. However, if contemporary \(S_B\) values are near the bottom of the range, then this approach could be overly optimistic, for it assumes that all historicallyobserved recruitment levels are possible, regardless of \(B_S\). The AGEPRO program can generate stochastic recruitments under model 3 based on thousands of input recruitment data points. Note that the empirical recruitment distribution model can be used to make deterministic projections by specifying a single observed recruitment.

Model 4. Two-Stage Empricial Recruits Per Spwaning Biomass

The two-stage recruits per spawning biomass model is a direct generalization of the \(R/B_S\) model where the spawning stock of the population is categorized into “low” and “high” states. The two-stage empirical recruits per spawning biomass distribution model depends on spawning biomass and is time-invariant.

In this model, there is an \(R/B_S\) distribution for the low spawning biomass state and an \(R/B_S\) distribution for the high spawning biomass state. Let \(G_{Low}\) be the cdf and let \(T_Low\) be the number of \(R/B_S\) values for the low \(B_S\) state. Similarly, let \(G_{High}\) be the cdf and let \(T_{High}\) be the number of \(R/B_S\) values for the high \(B_S\) state. Further, let \(B^*_S\) B denote the cutoff level of \(B_S\) such that, if \(B_S > B^*_S\), then \(B_S\) falls in the high state. Conversely if \(B_S <B^*_S\), then \(B_S\) falls in the low state. Recruitment is stochastically generated from \(G_{Low}\) or \(G_{High}\) using equations (\(\ref{eq:20}\)) and (\(\ref{eq:21}\)) dependent on the \(B_S\) state. The AGEPRO program can generate stochastic recruitments under model 4 based on thousands of input stock recruitment data points per \(B_S\) state.

Model 5. Beverton-Holt Curve with Lognormal Error

The Beverton-Holt curve (Beverton and Holt 1957) with lognormal errors is a parametric model of recruitment generation where survival to recruitment age is density dependent and subject to stochastic variation. The Beverton-Holt curve with lognormal error model depends on spawning biomass and is time-invariant.

The Beverton-Holt curve with lognormal error generates recruitment as

\[ \begin{split} \hat{r}(t) &= \frac{\alpha \cdot b_s (t-1)}{\beta + b_s (t-1)} \cdot e^w\ \\ & where\ w ~ \sim N \bigl(0, \sigma_{w}^2\bigr), \widehat{R}(t)= c_R \cdot \widehat{r}(t), and\ B_S(t) = c_{B}\cdot b_S(t) \\ \end{split} \tag{23}\label{eq:23} \]

The stock-recruitment parameters \(\alpha\), \(\beta\), and \(\sigma_w^2\) and the conversion coefficients for recruitment \(c_R\) and spawning stock biomass \(c_B\) are specified by the user. Here it is assumed that the parameter estimates for the Beverton-Holt curve have been estimated in relative units of recruitment \(r(t)\) and spawning biomass \(b_s(t)\), which are converted to absolute values using the conversion coefficients. Note that the absolute value of recruitment is expressed as numbers of fish, while for spawning biomass, the absolute value is expressed as kilograms of \(B_S\). For example, if the stock-recruitment curve was estimated with stock-recruitment data that were measured in millions of fish and thousands of metric tons of \(B_S\), then \(c_R=10^6\) and \(c_R=10^6\). It may be important to estimate the parameters of the stock-recruitment curve in relative units to reduce the potential effects of roundoff error on parameter estimates. It is important to note that the expected value of the lognormal error term is not unity but is \(exp\left(\dfrac{1}{2}\sigma_w^2\right)\). Therefore, in order to generate a recruitment model that has a lognormal error term that is equal to 1, one needs to multiply the parameter \(\alpha\) by \(exp\left(-\dfrac{1}{2}\sigma_w^2\right)\). This bias correction applies when the lognormal error used to fit the Beverton-Holt curve has a log-scale error term \(w\) with zero mean.

The Beverton-Holt curve is often reparameterized in a modified form with parameters for steepness \(h\) , unfished recruitment \(R_0\), and unfished spawning biomass \(B_0\). The modified Beverton-Holt curve produces \(h=R_0\) recruits when \(B_S = 0.2B_0\) and has the form

\[ \hat{R} = \frac{4hR_0B_S}{B_0(1-h)+B_S(5h-1)} \tag{24}\label{eq:24} \]

The parameters \(\alpha\) and \(\beta\) can be expressed as functions of the parameters of the modified Beverton-Holt curve as

\[ \alpha = \frac{4hR_0}{5h-1} = 4B_0\frac{h}{\left(\dfrac{B_0}{R_0}\right)(5h-1)} \tag{25}\label{eq:25} \] and \[ \beta = \frac{B_0(1-h)}{(5h-1)} - \frac{\alpha \left(\dfrac{B_0}{R_0}\right)(h^{-1}-1)}{4} \tag{26}\label{eq:26} \] Thus, parameter estimates for the modified curve can be used to determine the Beverton-Holt parameters for the AGEPRO program.

Model 6. Ricker Curve with Lognormal Error

The Ricker curve (Ricker 1954) with lognormal error is a parametric model of recruitment generation where survival to recruitment age is density dependent and subject to stochastic variation. The Ricker curve with lognormal error model depends on spawning biomass and is time invariant.

The Ricker curve with lognormal error generates recruitment as

\[ \begin{split} &\hat{r}(t) = \alpha \cdot b_S(t-1) \cdot e^{-\beta \cdot b_S(t-1)}\ \cdot e^w \\ & where\ w \sim N(0,\sigma_w^2), \hat{R}(t) = c_R \cdot \hat{r}(t),\ and\ B_S(t) = c_B \cdot b_S (t) \end{split} \tag{27}\label{eq:27} \]

The stock-recruitment parameters \(\alpha\), \(\beta\), and \(\sigma_w^2\) and the conversion coefficients for recruitment \(C_R\) and spawning stock biomass \(c_B\) are specified by the user. Here it is assumed that the parameter estimates for the Beverton-Holt curve have been estimated in relative units of recruitment \(r(t)\) and spawning biomass \(b_s(t)\) which are converted to absolute values using the conversion coefficients. It is important to note that the expected value of the lognormal error term is not unity but is \(exp\left(\dfrac{1}{2}\sigma^2_w\right)\). To generate a recruitment model that has a lognormal error term that is equal to 1, premultiply the parameter \(\alpha\) by \(exp\left(-\dfrac{1}{2}\sigma^2_w\right)\); this mean correction applies when the lognormal error used to fit the Ricker curve has a log-scale error term \(w\) with zero mean.

Model 7. Shepherd Curve with Lognormal Error

The Shepherd curve (Shepherd 1982) with lognormal error is a parametric model of recruitment generation where survival to recruitment age is density dependent and subject to stochastic variation. The Shepherd curve with lognormal error model depends on spawning biomass and is time-invariant.

The Shepherd curve with lognormal error generates recruitment as

\[ \begin{split} &\hat{r}(t) = \dfrac{\alpha \cdot b_S(t-1)}{ 1 + \left( \dfrac{b_s(t-1)}{k} \right)^{\beta} } \cdot e^w \\ &\\ & where\ w \sim N(0,\sigma_w^2), \hat{R}(t) = c_R \cdot \hat{r}(t),\ and\ B_S(t) = c_B \cdot b_S (t) \end{split} \tag{28}\label{eq:28} \]

The stock-recruitment parameters \(\alpha\), \(\beta\), \(k\), and \(\sigma^2_w\) and the conversion coefficients for recruitment \(c_R\) and spawning stock biomass \(c_B\) are specified by the user. Here it is assumed that the parameter estimates for the Beverton-Holt curve have been estimated in relative units of recruitment \(r(t)\) and spawning biomass \(b_S(t)\) which are converted to absolute values using the conversion coefficients. It is important to note that the expected value of the lognormal error term is not unity but is \(exp\left(\dfrac{1}{2}\sigma^2_w\right)\). To generate a recruitment model that has a lognormal error term that is equal to 1, premultiply the parameter \(\alpha\) by \(exp\left(-\dfrac{1}{2}\sigma^2_w\right)\); this mean correction applies when the lognormal error used to fit the Ricker curve has a log-scale error term w with zero mean.

Model 8. Lognormal Distribution

The lognormal distribution provides a parametric model for stochastic recruitment generation. The lognormal distribution model does not depend on spawning biomass and is time-invariant.

The lognormal distribution generates recruitment as

\[ \begin{split} & \hat{r}(t) = e^w \\ & \\ & where\ w \sim N (\mu_{log(r)},\sigma^2_{log(r)})\ and \ \hat{R}(t) = c_R \cdot \hat{r}(t) \end{split} \tag{29}\label{eq:29} \]

The lognormal distribution parameters \(\mu_{log(r)}\) and \(\sigma^2_{log(r)}\) as well as the conversion coefficient for recruitment \(C_R\) are specified by the user. It is assumed that the parameters of the lognormal distribution have been estimated in relative units \(r(t)\) and then converted to absolute values with the conversion coefficients.

Model 9. [DEPRECATED] Time-Varying Empirical Recruitment Distribution

This recruitment model has been deprecated. The model for a time-varying empirical recruitment distribution can be fully implemented using model 3.

Model 10. Beverton-Holt Curve with Autocorrelated Lognormal Error

The Beverton-Holt curve with autocorrelated lognormal errors is a parametric model of recruitment generation where survival to recruitment age is density dependent and subject to serially-correlated stochastic variation. The Beverton-Holt curve with autocorrelated lognormal error model depends on spawning biomass and is time-dependent.

The Beverton-Holt curve with autocorrelated lognormal error generates recruitment as

\[ \begin{split} \hat{r}(t) &= \frac{\alpha \cdot b_s (t-1)}{\beta + b_s (t-1)} \cdot e^{\varepsilon_t} \\ & where\ \varepsilon_t = \phi\varepsilon_{t-1} + w_t \ \ with\ w_t \sim N(0,\sigma^2_w), \\ & \hat{R}(t) = c_r \cdot \hat{r}(t), and\ B_s(t) = c_B \cdot b_S(t) \end{split} \tag{30}\label{eq:30} \]

The stock-recruitment parameters \(\alpha\), \(\beta\), \(\phi\), \(\varepsilon_0\), and \(\sigma^2_w\) and the conversion coefficients for recruitment \(c_R\) and spawning stock biomass \(c_B\) are specified by the user. The parameter \(\varepsilon_0\) is the log-scale residual for the stock-recruitment fit in the time prior to the projection. If this value is not known, the default is to set \(\varepsilon_0 = 0\).

Model 11. Ricker Curve with Autocorrelated Lognormal Error

The Ricker curve with autocorrelated lognormal error is a parametric model of recruitment generation where survival to recruitment age is density dependent and subject to serially correlated stochastic variation. The Ricker curve with autocorrelated lognormal error model depends on spawning biomass and is time-dependent.

The Ricker curve with autocorrelated lognormal error generates recruitment as

\[ \begin{split} \hat{r}(t) &= \alpha \cdot b_S(t-1) \cdot e^{-\beta\cdot b_S(t-1)} \cdot e^{\varepsilon_t} \\ & where\ \varepsilon_t = \phi\varepsilon_{t-1} + w_t \ \ with\ w_t \sim N(0,\sigma^2_w), \\ & \hat{R}(t) = c_r \cdot \hat{r}(t), and\ B_s(t) = c_B \cdot b_S(t) \end{split} \tag{31}\label{eq:31} \] The stock-recruitment parameters \(\alpha\), \(\beta\), \(\phi\), \(\varepsilon_0\), and \(\sigma^2_w\) and the conversion coefficients for recruitment \(c_R\) and spawning stock biomass \(c_B\) are specified by the user. The parameter \(\varepsilon_0\) is the log-scale residual for the stock-recruitment fit in the time prior to the projection. If this value is not known, the default is to set \(\varepsilon_0 = 0\).

Model 12. Shepherd Curve with Autocorrelated Lognormal Error

The Shepherd curve with autocorrelated lognormal error is a parametric model of recruitment generation where survival to recruitment age is density dependent and subject to serially-correlated stochastic variation. The Shepherd curve with autocorrelated lognormal error model depends on spawning biomass and is time-dependent.

The Shepherd curve with autocorrelated lognormal error generates recruitment as

\[ \begin{split} & \hat{r}(t) = \dfrac{\alpha \cdot b_S(t-1)}{ 1 + \left( \dfrac{b_s(t-1)}{k} \right)^{\beta} } \cdot e^{\varepsilon_t} \\ & \\ & where\ \varepsilon_t = \phi\varepsilon_{t-1} + w_t \ \ with\ w_t \sim N(0,\sigma^2_w), \\ & \hat{R}(t) = c_r \cdot \hat{r}(t), and\ B_s(t) = c_B \cdot b_S(t) \end{split} \tag{32}\label{eq:32} \] The stock-recruitment parameters \(\alpha\), \(\beta\), \(k\), \(\phi\), \(\varepsilon_0\), and \(\sigma^2_w\) and the conversion coefficients for recruitment \(c_R\) and spawning stock biomass \(c_B\) are specified by the user. The parameter \(\varepsilon_0\) is the log-scale residual for the stock-recruitment fit in the time prior to the projection. If this value is not known, the default is to set \(\varepsilon_0 = 0\).

Model 13. Autocorrelated Lognormal Distribution

The autocorrelated lognormal distribution provides a parametric model for stochastic recruitment generation with serial correlation. The autocorrelated lognormal distribution model does not depend on spawning biomass and is time-dependent.

The autocorrelated lognormal distribution is

\[ \begin{split} & n_r(t) = e^{s_t}\\ & \\ & where\ \varepsilon_t = \phi\varepsilon_{t-1} + w_t \ with \ w_t \sim N \left(\mu_{log(r)} \cdot \sigma^2_{log(r)}\right) \\ & and\ R(t) = c_R \cdot n_r(t) \end{split} \tag{33}\label{eq:33} \]

The lognormal distribution parameters \(\mu_{log(r)} \cdot \sigma^2_{log(r)}\),\(\phi\), \(\varepsilon_0\) as well as the conversion coefficient for recruitment \(C_R\) are specified by the user. It is assumed that the parameters of the lognormal distribution have been estimated in relative units \(r(t)\) and then converted to absolute values with the conversion coefficient.

Model 14. Empirical Cumulative Distribution Function of Recruitment

The empirical cumulative distribution function of recruitment can be used to randomly generates recruitment under the assumption that the distribution of \(R\) is stationary and independent of stock size. The empirical cumulative distribution function of recruitment model does not depend on spawning biomass and is time-invariant.

To describe this nonparametric approach, let \(R_S\) denote the \(S^{th}\) element in the ordered set of observed recruitment values. The empirical probability density function for \(R_S\), denoted as \(g(R_s)\), assigns a probability of \(1/T\) to each value of \(R(t)\) where \(T\) is the number of stock-recruitment data points. Let \(G(R_S)\) denote the cumulative distribution function (cdf) for \(R_S\). Set the values of \(G\) at the minimum and maximum observed \(R_S\) to be \(G(R_{min})=0\) and \(G(R_{max})=0\) so that the cdf of \(R_S\) can be written as

\[ G(R_S) = \frac{s-1}{T-1} \tag{34}\label{eq:34} \]

Random values of \(R_S\) can be generated by applying the probability integral transform to the empirical cdf. To do this, let \(U\) be a uniformly distributed random variable on the interval [0,1]. The value of \(R_S\) corresponding to a randomly chosen value of \(U\) is determined by applying the inverse function of the cdf \(G(R_S)\). In particular, if \(U\) is an integer multiple of \(1/(T-1)\) so that \(U = s/(T-1)\) then \(\widehat{R_s} = G^{-1}(U)\). Otherwise \(\widehat{R_s}\) can be obtained by linear interpolation when \(U\) is not a multiple of \(1/(T-1)\).

In particular, if \((s-1) / (T-1) < U < s/ (T-1)\), then \(\widehat{R_s}\) is interpolated between \(R_S\) and \(R_{S+1}\) as

\[ U = \left(\frac{\dfrac{s}{T-1}-\dfrac{s-1}{T-1}}{R_{S+1} - R_S} \right)\left(\widehat{R_S}-R_S\right) + \frac{s-1}{T-1} \tag{35}\label{eq:35} \]

Solving for \(\widehat{R_s}\) as a function of \(U\) yields

\[ \widehat{R_S} = (T-1)(R_{S+1}-R_S)\left(U-\frac{s-1}{T-1}\right)+R_S \tag{36}\label{eq:36} \]

where the interpolation index \(s\) is determined as the greatest integer in \(1+U(T-1)\). Given the value of \(R_S\), recruitment is set to be

\[ \widehat{R}(t) = \widehat{R_S} \tag{37}\label{eq:37} \] The AGEPRO program can generate stochastic recruitments under model 14 based on thousands of input recruitment data points.

Model 15. Two Stage Empirical Cumulative Distribution Function of Recruitment

The two-stage empirical cumulative distribution function of recruitment model is an extension of Model 14 where the spawning stock of the population is categorized into low and high states. The two-stage empirical cumulative distribution function of recruitment model depends on spawning biomass and is time-invariant.

In this model, there is an empirical recruitment distribution \(R_{Low}\) for the low spawning biomass state and an empirical recruitment distribution \(R_{High}\) for the high spawning biomass state. Let \(G_{Low}\) be the cdf and let \(T_{Low}\) be the number of \(R\) values for the low \(B_S\) state. Similarly, let \(G_{High}\) be the cdf and let \(T_{High}\) be the number of \(R\) values for the high \(B_S\) state. Further, let \(B_S^*\) denote the cutoff level of \(B_S\) such that, if \(B_S > B_S^*\), then \(B_S\) falls in the high state. Conversely if \(B_S < B_S^*\), then \(B_S\) falls in the low state. Recruitment is stochastically generated from \(G_{Low}\) or \(G_{High}\) using equations (\(\ref{eq:36}\)) and (\(\ref{eq:37}\)) dependent on the \(B_S\) state. The AGEPRO program can generate stochastic recruitments under model 15 based on thousands of input stock-recruitment data points per \(B_S\) state.

Model 16. Linear Recruits Per Spawning Biomass Predictor with Normal Error

The linear recruits per spawning biomass predictor with normal error is a parametric model to simulate random values of recruits per spawning biomass \(\dfrac{R}{B_S}\) and realized recruitment values. The predictors in the linear model \(X_p(t)\) can be any continuous variable and may typically be survey indices of cohort abundance or environmental covariates that are correlated with recruitment strength. Input values of each predictor are required for each time period. If a value of a predictor is missing or not known for one or more periods, the missing values can be imputed using appropriate measures of central tendency, e.g., mean or median values. Similarly, if this model has zero probability in a given time period (e.g., is not a member of the set of probable models), then dummy values can be input for each predictor. For each time period and simulation, a random value of \(\dfrac{R}{B_S}\) is generated using the linear model

\[ \dfrac{R}{B_S} = \beta_0 + \sum^{N_P}_{p=1}\beta_p\ \cdot X_p(t) + \varepsilon \tag{38}\label{eq:38} \]

where \(N_p\) is the number of predictors, \(\beta_0\) is the intercept, \(\beta_p\) is the linear coefficient of the \(p^{th}\) predictor and \(\varepsilon\) is a normal distribution with zero mean and constant variance \(\sigma^2\). It is possible negative values of \(\dfrac{R}{B_S}\) to be generated using this formulation; such values are excluded from the set of simulated values of \(\dfrac{R}{B_S}\) from equation (\(\ref{eq:35}\)) by testing if \(\dfrac{R}{B_S} < 0\) repeating the random sampling until an feasible positive value of \(\dfrac{R}{B_S}\) is obtained. This model randomly generates \(\dfrac{R}{B_S}\) values under the assumption that the linear predictor of the \(\dfrac{R}{B_S}\) ratio is stationary and independent of stock size. Random values of \(\dfrac{R}{B_S}\) are multiplied by realized spawning biomass to generate recruitment in each time period. The linear recruits per spawning biomass predictor with normal error depends on spawning biomass and is time-invariant unless time is used as a predictor.

Model 17. Loglinear Recruits Per Spawning Biomass Predictor with Lognormal Error

The loglinear recruits per spawning biomass predictor with lognormal error is a parametric model to simulate random values of recruits per spawning biomass \(\dfrac{R}{B_S}\) and associated random recruitments. Predictors for the loglinear model \(X_p(t)\) can be any continuous variable and could include survey indices of cohort abundance or environmental covariates that are correlated with recruitment strength. Input values of each predictor are required for each time period. If a value of a predictor is missing or not known for one or more periods, the missing values can be imputed using appropriate measures of central tendency, e.g., mean or median values. If this model has zero probability in a given time period (e.g., is not a member of the set of probable models), then dummy values can be input for each predictor. For each time period and simulation, a random value of the natural logarithm of \(\dfrac{R}{B_S}\) is generated using the loglinear model

\[ \log \left(\dfrac{R}{B_s}\right) = \beta_0 + \sum\limits_{p=1}^{N_p} \beta_p\ \cdot X_p(t) + \varepsilon \tag{39}\label{eq:39} \]

where \(N_P\) is the number of predictors, \(\beta_0\) is the intercept, \(\beta_p\) is the linear coefficient of the \(p^{th}\) predictor and \(\varepsilon\) is a normal distribution with constant variance \(\sigma^2\) and mean equal to \(-0.5\sigma^2\). In this case, the mean of \(\varepsilon\) implies that the expected value of the lognormal error term is unity. This model generates positive random values of \(\dfrac{R}{B_S}\) under the assumption that the linear predictor of the \(\dfrac{R}{B_S}\) ratio is stationary and independent of stock size. Simulated values of \(\dfrac{R}{B_S}\) are multiplied by realized spawning biomass to generate recruitment in each time period. The loglinear recruits per spawning biomass predictor with lognormal error depends on spawning biomass and is time-invariant unless time is used as a predictor.

Model 18. Linear Recruitment Predictor with Normal Error

The linear recruitment predictor with normal error is a parametric model to simulate representative random values of recruitment. The predictors in the linear model \(X_p(t)\) can be any continuous variable and could represent survey indices of cohort abundance or environmental covariates correlated with recruitment strength. Input values of each predictor are required for each time period. If a value of a predictor is missing or not known for one or more periods, the missing values can be imputed using appropriate measures of central tendency, e.g., mean or median values. Similarly, if this model has zero probability in a given time period (e.g., is not a member of the set of probable models), then dummy values can be input for each predictor. For each time period and simulation, a random value of R is generated using the linear model

\[ \begin{split} & n_r(t) = \beta_0+\sum\limits_{p=1}^{N_p}\beta_p\ \cdot X_p(t) + \varepsilon\\ & with \ R(t) = c_r \cdot n_r(t)\\ \end{split} \tag{40}\label{eq:40} \]

here \(N_p\) is the number of predictors, \(\beta_0\) is the intercept, \(\beta_p\) is the linear coefficient of the \(p^{th}\) predictor and \(\varepsilon\) is a normal distribution with zero mean and constant variance \(\sigma^2\) and the conversion coefficients for recruitment is \(c_R\). It is possible that negative values of \(R\) can be generated using this formulation; such values are excluded from the set of simulated values of \(R\) from equation (\(\ref{eq:37}\)) by testing if \(R\) repeating the random sampling until an feasible positive value of \(R\) is obtained. This model randomly generates \(R\) values under the assumption that the linear predictor of \(R\) is stationary and independent of stock size. The linear recruitment predictor with normal error does not depend on spawning biomass and is time-invariant unless time is used as a predictor.

Model 19. Loglinear Recruitment Predictor with Lognormal Error

The loglinear recruitment predictor with lognormal error is a parametric model to simulate random values of recruitment R. Predictors for the loglinear model \(X_p(t)\) can be any continuous variable such as survey indices of cohort abundance or environmental covariates that are correlated with recruitment strength. Input values of each predictor are required for each time period. If a value of a predictor is missing or not known for one or more periods, the missing values can be imputed using appropriate measures of central tendency, e.g., mean or median values. If this model has zero probability in a given time period (e.g., is not a member of the set of probable models), then dummy values can be input for each predictor. For each time period and simulation, a random value of the natural logarithm of \(R\) is generated using the loglinear model

\[ \begin{split} & \log{(n_r(t))} = \beta_0+\sum\limits_{p=1}^{N_p}\beta_p\ \cdot X_p(t) + \varepsilon\\ & with \ R(t) = c_r \cdot n_r(t)\\ \end{split} \tag{41}\label{eq:41} \] here \(N_p\) is the number of predictors, \(\beta_0\) is the intercept, \(\beta_p\) is the linear coefficient of the \(p^{th}\) predictor and \(\varepsilon\) is a normal distribution with constant variance \(\sigma^2\) and and mean equal to \(-0.5\sigma^2\), and the conversion coefficients for recruitment is \(c_R\). In this case, the mean of \(\varepsilon\) implies that the expected value if the lognormal error term is unity. This model generates positive random values of \(R\) under the assumption that the linear predictor of \(R\) is stationary and independent of stock size. The loglinear recruitment predictor with lognormal error does not depend on spawning biomass and is time-invariant unless time is used as a predictor.

Model 20. Fixed Recruitment

The fixed recruitment predictor applies a specified value of recruitment for each year of the time horizon. The vector of input recruitment values in relative units is \(\underline{n}_r = \Big[ n_r(1), n_r(2, ..., n_r(Y)\Big]\) for projections years \(1\) through \(Y\). The fixed recruitment model predicts recruitment as

\[ R(t) = c_r \cdot n_r(t) \tag{42}\label{eq:42} \]

where the conversion coefficient for input recruitment to absolute numbers is \(c_R\). The fixed recruitment model does not depend on spawning biomass and is time-invariant.

Model 21. Empirical Cumulative Distribution Function of Recruitment with Linear Decline to Zero

The empirical cumulative distribution function of recruitment with linear decline to zero model is an extension of the empirical cumulative distribution function of recruitment (Model 14) in which recruitment strength declines when the spawning stock biomass falls below a threshold \(B^*_S\). The decline in recruitment randomly generated from the empirical cdf of the observed recruitment is proportional to the ratio of current spawning stock biomass to the threshold \(\dfrac{B_S}{B^*_S}\) when \(B_S < B^*_S\). In particular, predicted recruitment values are randomly generated using equation (\(\ref{eq:37}\)) given the input time series of observed recruitment. That is, if the current spawning biomass is at or above the threshold with \(B_S \ge B^*_S\) then the predicted recruitment is

\[ R = c_R \cdot \left[ (T-1)(R_{S+1} - R_S)\left( U-\dfrac{s-1}{T-1}\right)+R_S\ \right] \tag{43}\label{eq:43} \]

where the conversion coefficient for input recruitment to absolute numbers is \(c_R\).

Otherwise, if the current spawning biomass falls below the threshold with \(B_S < B^*_S\) then the predicted recruitment is reduced to be

\[ R = c_R \cdot \frac{B_S}{B^*_S} \left[\ (T-1)(R_{S+1} - R_S)\left( U-\dfrac{s-1}{T-1}\right)+R_S\ \right] \tag{44}\label{eq:44} \]

where the conversion coefficient for input recruitment to absolute numbers is \(c_R\). The empirical cumulative distribution function of recruitment with linear decline to zero model depends on spawning biomass and is time-invariant.