The VPA model reconstructs the historical population based on an expanded catch at age matrix and an index of abundance. From an estimate of the fishing mortality in the terminal year, historical abundance and fishing mortality are back-calculated. Natural mortality is a required input parameter. For a biomass-based index, weight at age is also required. Maturity at age is also required if spawning biomass is to be calculated. The dynamics equations are based on VPA-2BOX (Porch 2018), although 2-stock mixing and diffusion processes are not modeled here. The maximum age in the model is a plus-group accumulator age.
Dynamics equations
In this model, fishing mortality \(F_{a,t}\) of age \(a\) in year \(t\) is generally solved from the observed catch of the corresponding age and year \(C_{a,t}\) and the estimated abundance of the same cohort in the following year \(N_{a+1,t+1}\): \[F_{a,t} = \dfrac{Z_{a,t}C_{a,t}}{[\exp(Z_{a,t}) - 1]N_{a+1,t+1}},\] where \(Z_{a,t} = F_{a,t} + M_a\), and \(M\) is natural mortality. Equation 1 is Baranov’s equation with the substitution: \(N_{a,t}=N_{a+1,t+1}\exp(Z_{a,t})\). Newton’s method is used to numerically solve Equation 1.
With \(F_{a,t}\), the abundance of the corresponding age and year is solved also from Baranov’s equation: \[N_{a,t} = \dfrac{C_{a,t}Z_{a,t}}{F_{a,t}[1 - \exp(-Z_{a,t})]}.\]
There are two exceptions to this algorithm. First, for the terminal year \(t=T\), F is estimated in the model, followed by the calculations for the terminal year abundances with Equation 2. This starts the back-calculations for \(F\) and \(N\) for the preceding years with equations 1 and 2, respectively. The terminal F-at-age vector \(F_{a,T}\) can be free parameters with age or constrained as a logistic or double-normal (dome) function with age.
Second, the fishing mortality rates of the plus-group age \(a=A\) and the preceding age \(a=A-1\) are subject to additional constraints and are solved separately. First, note that the abundance of the plus-group in year \(t+1\) must satisfy the equation \[ N_{A,t+1} = N_{A,t}\exp(-Z_{A,t}) + N_{A-1,t}\exp(-Z_{A-1,t}),\] with \[\begin{align} N_{A,t} &= \dfrac{C_{A,t}Z_{A,t}}{F_{A,t}[1 - \exp(-Z_{A,t})]}\\ N_{A-1,t} &= \dfrac{C_{A-1,t}Z_{A-1,t}}{F_{A-1,t}[1 - \exp(-Z_{A-1,t})]} \end{align}.\]
An additional constraint is needed for the plus-group \(F\), where \[F_{A,t} = \alpha F_{A-1,t},\] where \(\alpha\) is a numeric constant either fixed (often to 1) or estimated in the model.
From \(N_{A,t+1}\), \(C_{A,t}\), \(C_{A-1,t}\), \(M_A\), and \(M_{A-1}\), we can solve \(F_{A-1,t}\) and \(F_{A,t}\) as the values that satisfy equations 3-6. The corresponding abundances \(N_{A-1,t}\) and \(N_{A,t}\) are then obtained with equation 2.
Fishery selectivity \(v_{a,t}\) can derived from \(F_{a,t}\): \[v_{a,t} = \dfrac{F_{a,t}}{\max \limits_a F_{a,t}}\]
Likelihoods
The VPA is tuned to an index via maximum likelihood. The predicted index \(s\) is \[ \hat{I}_{s,t} = \hat{q_s}\sum_a N_{a,t} \tilde{v}_{a,s} w_a\] and \[ \hat{I}_{s,t} = \hat{q_s}\sum_a N_{a,t} \tilde{v}_{a,s}\] for a biomass-based and abundance-based index, respectively. Index selectivity \(\tilde{v}_{a,s}\) is fixed in the model.
The likelihood of the observed indices is
\[L^I = \sum_s \lambda_s \sum_t \left(-\log(\sigma_{s,t}) - 0.5\left[\dfrac{\log(I_{s,t}) - \log(\hat{I}_{s,t})}{\sigma_{s,t}}\right]^2\right),\] with \(\lambda_s\) as the weighting coefficient for index \(s\).
Additional constraints
In any VPA model, estimates of F and N of the youngest age classes are highly uncertain. To stabilize estimates, random walk penalties can be applied to the recruitment and vulnerability estimates in the most recent years: \[ \begin{align} \textrm{pen}_N &= \sum_{i=T-\tilde{t}}^T \left(-\log(\sigma_{\tilde{a}}) - 0.5 \left[\dfrac{\log(N_{\tilde{a},i}) - \log(N_{\tilde{a},i-1})}{\sigma_{\tilde{a}}}\right]^2\right) \\ \textrm{pen}_v &= \sum_{i=T-\tilde{t}}^T \left(-\log(\sigma_v) - 0.5 \left[\dfrac{\log(v_{a,i}) - \log(v_{v,i-1})}{\sigma_v}\right]^2\right) \end{align},\] where \(\tilde{a}\) is the smallest age in the VPA.
Note that these penalties can significantly alter model results, and generate spurious trends in F and N especially if there are significantly shifts in recent vulnerability and recruitment. The implications of incorporating these random walk penalties need to be evaluated on a case-by-case basis.
Reference points
From the VPA output, reference points can be calculated. Estimates of recruitment (the youngest age class \(\tilde{a}\)) and spawning biomass are used to estimate a stock-recruitment relationship to obtain MSY and unfished reference points.
References
Porch, C.E. 2018. VPA-2BOX 4.01 User Guide. NOAA Tech. Memo. NMFS-SEFSC-726. 67 pp.