On the Benefits of Including Age-Structure in Harvest Control Rules

Abstract

This paper explores the benefits of including age structure in the control rule (HCR) when decision makers regard their (age-structured) models as approximations. We find that introducing age structure into the HCR reduces both the volatility of the spawning biomass and the yield. Although the benefits are lower at a fairly imprecise level, there are still major advantages for the actual precision with which the case study is assessed. Moreover, we find that when age-structure is included in the HCR the relative ranking of different policies in terms of variance in biomass and yield does not differ. These results are shown both theoretically and numerically by applying the model to the Southern Hake fishery.

Appendix

1.1 HCR and stock dynamics in the three age-class model

The steps used for the three age model are equivalent to the ones described in “Appendix”. We explicitly define the matrices that solve for the optimal HCR and the stock dynamics. By solving the following system of equations given by the Riccati equation, we obtain the solution stated in Sect. 3.

The HCR is obtained by solving the following Riccati equation

$$\begin{aligned} \mathbf P&= \left[ \begin{array}{cccc} P_{11} &{}P_{12} &{} P_{13} &{} P_{14} \\ P_{12} &{}P_{22} &{} P_{23} &{} P_{24} \\ P_{13} &{}P_{23} &{} P_{33} &{} P_{34} \\ P_{14} &{}P_{24} &{} P_{34} &{} P_{44} \end{array}\right] \\&= \left[ \begin{array}{ll} \beta P_{11} - \frac{(beta^2 P_{14}^2 s2^2)}{(\beta P_{44} s2^2 - 1)} &{} \beta ( P_{13} + P_{12} \rho ) - \frac{(\beta ^2 P_{14} s2^2 ( P_{34} + P_{24} \rho ))}{(\beta P_{44} s2^2 - 1)} \\ \beta ( P_{13} \!+\! P_{12} \rho ) \!-\! \frac{(\beta ^2 P_{14} s2^2 ( P_{34} \!+\! P_{24} rho))}{(\beta P_{44} s2^2 \!-\! 1)}&{} \beta P_{33} \!-\! hxx11 {\uplambda }\!+\! \beta P_{22} \rho ^2 \!+\! 2 \beta P_{23} \rho \!-\! \frac{(\beta ^2 s2^2 ( P_{34} + P_{24} rho)^2)}{(\beta P_{44} s2^2 - 1)}\\ -\frac{(\beta P_{14})}{(\beta P_{44} s2^2 - 1)}&{} - hxx12 {\uplambda }- \frac{(\beta ( P_{34} + P_{24} \rho ))}{(\beta P_{44} s2^2 - 1)}\\ 0 &{}-hxx13 {\uplambda }\end{array}\right. \\&\qquad \qquad \left. \begin{array}{ll} -\frac{(\beta P_{14})}{(\beta P_{44} s2^2 - 1)} &{} 0\\ - hxx12 {\uplambda }- \frac{(\beta ( P_{34} + P_{24} \rho ))}{(\beta P_{44} s2^2 - 1)} &{} -hxx13 {\uplambda }\\ - hxx22 {\uplambda }- \frac{(\beta P_{44})}{(\beta P_{44} s2^2 - 1} &{} -hxx23 {\uplambda }\\ -hxx23 {\uplambda }&{}-hxx33 {\uplambda }\end{array}\right] \end{aligned}$$

We solve the equation by using a symbolic Matlab code. Given $\mathbf P$ we compute $G=-\left( \mathbf Q + \mathbf B ^T \mathbf P \mathbf B \right) ^{-1} \mathbf B ^T \mathbf P \mathbf A$ to obtain the optimal harvesting control rule. Thus

$$\begin{aligned} \Delta F =\mathbf G \mathbf x&= \left[ 0 \displaystyle \frac{\beta {\uplambda }p_2 (\mu _{x_1} e^{x_{1,tar}} \mu _{x_2} e^{x_{2,tar}} + \mu _{z} e^{z_{tar}} \mu _{x_2} e^{x_{2,tar}} \rho )}{\beta {\uplambda }\mu ^2_{x_2} e^{2x_{2,tar}} p_2^2 + 1} \displaystyle \frac{\beta {\uplambda }\mu ^2_{x_2} e^{2x_{2,tar}} p_2}{\beta {\uplambda }\mu ^2_{x_2} e^{2x_{2,tar}} p_2^2 + 1} 0 \right] \\&\times \left[ \begin{array}{c} 1 \\ z_t \\ \Delta x_{1,t}\\ \Delta x_{2,t} \end{array}\right] \end{aligned}$$

Given $G$, we introduce the stochastic shocks to generate the optimal trajectories by using $\mathbf x '= \mathbf A \mathbf x + \mathbf B \mathbf G \mathbf x + \varepsilon '=\mathbf D \mathbf x + \varepsilon '$ where

$$\begin{aligned} \mathbf D = \left[ \begin{array}{cccc} 1 &{} 0 &{} 0 &{} 0\\ 0 &{} \rho &{} 0 &{} 0\\ 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} -\displaystyle \frac{\beta {\uplambda }p_2 (\mu _{x_1} e^{x_{1,tar}} \mu _{x_2} e^{x_{2,tar}} + \mu _{z} e^{z_{tar}} \mu _{x_2} e^{x_{2,tar}} \rho )}{\beta {\uplambda }\mu ^2_{x_2} e^{2x_{2,tar}} p_2^2 + 1} &{}\displaystyle \frac{1}{\beta {\uplambda }\mu ^2_{x_2} e^{2x_{2,tar}} p_2^2 + 1} &{} 0 \end{array}\right] \end{aligned}$$

1.2 Proof of Proposition 1

From the expression for the optimal LQ HCR, if $z=0$, then the optimal fishing mortality and stock level is the same as in the heuristic case, with $B_{tar}=B$ and $F=F_{tar}$. However, if we take into account the recruitment structure, $z>0$, then $\Delta F=F-F_{tar}$. Therefore, the two optimal policies are no longer equivalent. Additionally assume that $\hat{\mu }_z$ is a low value, $\hat{\mu }_z \rightarrow 0$, then the optimal LQ HCR sets a higher fishing mortality, meaning that heuristic HCR underestimate fishing mortality. Then, If $\hat{\mu }_z \rightarrow 0$ then $\Delta F= F-F_{tar}= \displaystyle \frac{\Theta }{p\Theta +1} \left( \frac{\hat{\mu }_x^1}{\hat{\mu }_x^2 }\right) z$. $\square$

1.3 Proof of Proposition 2

Using the standard formula for calculating the variance of a random variable, the Biomass volatility is given by:

$$\begin{aligned} Var(B)&= \frac{1}{8\sigma ^3}\left[ 4\sigma ^2\text {sinh}(2\sigma ) \left[ \hat{\mu }^2_{x^1}+\hat{\mu }^2_{z}\right] + 16\hat{\mu }_{x^1}\hat{\mu }_{z}\sigma \text {sinh}(\sigma )^2\right. \\&\left. +\, \frac{2\hat{\mu }^2_{x^2}\sigma \text {sinh}(2\sigma \theta _1)\text {sinh}(2\sigma \theta _2)}{\theta _1\theta _2} +\frac{16\hat{\mu }_{x^2}\hat{\mu }_{z}\text {sinh}(\sigma \theta _1)\text {sinh}(\sigma \theta _2)\text {sinh}(\sigma )}{\theta _1\theta _2} \right. \\&\left. +\, \frac{8\hat{\mu }_{x^1}\hat{\mu }_{x^2}\sigma \text {sinh}(\sigma \theta _2)(e^{2\sigma + 2\sigma \theta _1} - 1)}{\theta _2(\theta _1 + 1)e^{\sigma + \sigma \theta _1}}\right] -E(B) \end{aligned}$$

where $\text {sinh}(a)=\frac{e^a-e^{-a}}{2}$ is the hyperbolic sine function. Now consider a HCR that implements $F_{tar}$ for all possible states, such as a constant catch rule. In our model this is equivalent to calculating the solution of the optimal HCR when ${\uplambda }=0$. In that case $\theta _1=0$ and $\theta _2=1$. Then

$$\begin{aligned} Var(B|{\uplambda }=0)&= \frac{1}{8\sigma ^3}\left[ 4\sigma ^2\text {sinh}(2\sigma ) \left[ \hat{\mu }^2_{x^1}+\hat{\mu }^2_{z}\right] + 16\hat{\mu }_{x^1}\hat{\mu }_{z}\sigma \text {sinh}(\sigma )^2 \right. \\&\left. +\, \frac{8\hat{\mu }_{x^1}\hat{\mu }_{x^2}\sigma \text {sinh}(\sigma )(e^{2\sigma } - 1)}{e^{\sigma }} +\lim _{\theta _1\rightarrow 0}\left( \frac{2\hat{\mu }^2_{x^2}\sigma \text {sinh}(2\sigma \theta _1)\text {sinh}(2\sigma )}{\theta _1} \right. \right. \\&\left. +\,\left. \frac{16\hat{\mu }_{x^2}\hat{\mu }_{z}\text {sinh}(\sigma \theta _1)\text {sinh}(\sigma )^2}{\theta _1}\right) \right] -E(B) \end{aligned}$$

Applying l’Hopital Rule, we have

$$\begin{aligned} Var(B|{\uplambda }=0)&= \frac{1}{8\sigma ^3}\left[ 4\sigma ^2\text {sinh}(2\sigma ) \left[ \hat{\mu }^2_{x^1}+\hat{\mu }^2_{z}\right] + 16\hat{\mu }_{x^1}\hat{\mu }_{z}\sigma \text {sinh}(\sigma )^2 \right. \\&\left. +\, \frac{8\hat{\mu }_{x^1}\hat{\mu }_{x^2}\sigma \text {sinh}(\sigma )(e^{2\sigma } - 1)}{e^{\sigma }} +2\hat{\mu }^2_{x^2}\sigma \text {cosh}(0)\text {sinh}(2\sigma ) \right. \\&\left. +\, 16\hat{\mu }_{x^2}\hat{\mu }_{z}\text {cosh}(0)\text {sinh}(\sigma )^2\right] -E(B) \end{aligned}$$

Then

$$\begin{aligned} Var(B|{\uplambda }\!=\!0)\!-\!Var(B)&= \frac{1}{8\sigma ^3}\left[ 8\hat{\mu }_{x^1}\hat{\mu }_{x^2}\sigma \!\left( \!\frac{ \text {sinh}(\sigma )(e^{2\sigma } \!-\! 1)}{e^{\sigma }} \!-\!\frac{ \text {sinh}(\sigma \theta _2)(e^{2\sigma + 2\sigma \theta _1} \!-\! 1)}{\theta _2(\theta _1 + 1)e^{\sigma + \sigma \theta _1}}\!\right) \right. \\&\left. +\,2\hat{\mu }^2_{x^2}\sigma \left( \text {cosh}(0)\text {sinh}(2\sigma ) - \frac{ \text {sinh}(2\sigma \theta _1)\text {sinh}(2\sigma \theta _2)}{\theta _1\theta _2} \right) \right. \\&\left. +\,16\hat{\mu }_{x^2}\hat{\mu }_{z}\left( \text {cosh}(0)\text {sinh}(\sigma )^2\!-\! \frac{\text {sinh}(\sigma \theta _1)\text {sinh}(\sigma \theta _2)\text {sinh}(\sigma )}{\theta _1\theta _2} \!\right) \!\right] \end{aligned}$$

Given that $\forall x \in [-1,1]$, $\text {sinh}(\sigma ) - \frac{\text {sinh}(\sigma x)}{x}$ and $\frac{ (e^{2\sigma } - 1)}{e^{\sigma }} >\frac{ (e^{2\sigma + 2\sigma \theta _1} - 1)}{(\theta _1 + 1)e^{\sigma + \sigma \theta _1}}$ (see Fig. 6)

$$\begin{aligned} Var(B|{\uplambda }=0)-Var(B)&> \frac{1}{8\sigma ^3}\left\{ 2\hat{\mu }^2_{x^2}\sigma \text {sinh}(2\sigma ) \left[ \text {cosh}(0) - \text {sinh}(2\sigma ) \right] \right. \\&\left. +\,16\hat{\mu }_{x^2}\hat{\mu }_{z} \text {sinh}(\sigma )^2 \left[ \text {cosh}(0)- \text {sinh}(\sigma ) \right] \right\} \end{aligned}$$

Note that the limiting variance is positive if the hyperbolic cosine function in sinh1 (1) zero is 1. Then $Var(B|{\uplambda }=0)-Var(B)>0$ if $\sigma \le \displaystyle \frac{\text {sinh}^{-1}(1)}{2}=0.4407$ $\square$

Fig. 6

Hyperbolic sin function properties. The left hand side shows that $\forall x \in [-1,1]$, $\text {sinh}(\sigma ) - \frac{\text {sinh}(\sigma x)}{x}>0$. The right hand side shows that $\forall \theta _1 \in [-1,0]$, $\frac{ (e^{2\sigma } - 1)}{e^{\sigma }} -\frac{ (e^{2\sigma + 2\sigma \theta _1} - 1)}{(\theta _1 + 1)e^{\sigma + \sigma \theta _1}}>0.$