Near-To-Optimal Harvesting Strategiesfor a Stochastic Multicohort Fishery

(1)

SNF-REPORT NO. 57/1999

Near-To-Optimal Harvesting Strategies for a Stochastic Multicohort Fishery

by

Erling Moxnes

SNF-project No. 5774

«Flerbestandsforvaltning under usikkerhet»

This project is financed by the Research Council of Norway

Centre for Fisheries Economics Report No. 56

FOUNDATION FOR RESEARCH IN ECONOMICS AND BUSINESS ADMINISTRATION BERGEN, NOVEMBER 1999

Ytterligere eksemplarfremstilling uten avtale og i strid med åndsverkloven er straffbart og kan medføre erstatningsansvar.

(2)

ISBN 82-491-0010-7 ISSN 0803-4036

(3)

FOREWORD

This is a report in a series dealing with optimal management of multispecies fisheries under uncertainty. The case for all reports is the Barents Sea, where cod and capelin are the two main species being harvested. In this report we investigate near-optimal harvesting strategies for a multicohort representation of cod, given natural variability and measurement error. The results are compared to previous results from yield-per- recruit analysis and from stochastic optimization in aggregate surplus growth models.

The project is financed by the Research Council of Norway, under the program Marine Resource Management.

(4)

TABLE OF CONTENTS

FOREWORD I

1. INTRODUCTION 1

2. A MULTICOHORT MODEL 3

3. PARAMETER ESTIMATES 7

3.1. Recruitment 7

3.2. Weight (and length) 8

3.3. Mortality, fishing selectivity and environmental variation 19

3.4. Summary of model parameters 29

3.5. Aggregate model 30

4. STOCHASTIC OPTIMIZATION IN POLICY SPACE, SOPS 33

5. FROM YPR TO A SIMPLE FULL MODEL STRATEGY 36

6. NEAR-TO-OPTIMAL FISHING STRATEGIES 41

6.1. Non-linear one-dimensional policies 41

6.2. Linear multi-dimensional policies 42

6.3. Non-linear two-dimensional policies 44

6.4. Non-linear three-dimensional policies 45

7. COMPARISONS TO AGGREGATE MODELS 47

8. MEASUREMENT ERROR 50

9. SENSITIVITY TO NATURAL MORTALITY 54

10. CONCLUSIONS 55

REFERENCES 57

APPENDIX 59

(5)

1. INTRODUCTION

In this report we search harvesting strategies for a multicohort fish resource. The case is cod in the Barents Sea. According to control theory, optimal fishing strategies should in general make use of information about the size of all age classes, see e.g. Mendelssohn (1978). Except for in highly simplified cases, such strategies are not easily found. Hence, various simplifications are resorted to.

Of greatest practical importance to date is the yield-per-recruit method. This method is used to find fishing mortalities which produce a maximal (Fmax) or an otherwise desirable yield per recruit (e.g. F0.1). The method takes a certain fishing pattern for granted, the fishing mortality does not reflect the current age distribution, it does not explicitly consider environmental disturbances, it does not reflect stock effects on weights and on catch per unit effort, and it does not reflect harvesting costs and fish prices. Nor does the method take account of the stock effect on recruitment. However, the latter effect is normally accounted for by the use of a separate stock-recruitment model, indicating upper limits for fishing mortality (Fmed or Fhigh).

Here we use a numerical method termed “stochastic optimization in policy space”

(SOPS) to find simplified harvesting strategies for a complex multicohort model. This method has two major advantages over other methods: It is intuitively appealing, and it can be used to find near-to-optimal solutions to highly complex problems. In short the method works as follows. We search for a fishing strategy, i.e. a function which relates measurements of biomass to next year’s quotas. To begin with we make a guess about this function. Then we simulate a model of the fish population and of the fleet economics to see how well the policy function does. Since natural variation (stochasticity) is included in the model, repeated simulations are made to see how well the policy does on average (Monte Carlo simulations). Then a minor change is made in the policy function, and again simulations are used to evaluate how well the new policy does. This process goes on until no further improvements are made. The method represents only a minor deviation from the ordinary use of simulation models. It is related to an array of new methods currently applied to complex decision problems like for instance chess, see e.g.

Bertsekas and Tsitsiklis (1996). Due to the method’s ability to tackle complexity, the most important restricting assumptions¹ made when using the yield-per-recruit method are relaxed one by one, except the one about a given fishing pattern. For the complex

1 Similar to other one-species cohort models and policy analyses, we do not consider a vast array of other complicating assumptions like species interactions and spatial distribution. While these effects could be tackled by the method, they are left for further research.

(6)

problem at hand, the method will only provide approximate solutions. However, through a progression of policy improvements we are able to indicate rapidly diminishing returns to further refinement.

The approximate solutions we provide might seem inferior to exact results obtained by analytical methods or stochastic dynamic programming, see e.g. Mendelssohn (1978), Spulber (1983), Spulber (1985), and Naqib and Stollery (1982). However, one must keep in mind that existing exact results rely on highly simplifying assumptions about the problem at hand. While these results provide basic insight about the studied models, they may be of limited value for practical decision making since the problem of blending and extrapolating the basic results to practical situations is very difficult. Non-linearities and feedback underlay this difficulty.

Section 2 presents a multicohort model of the cod fisheries of the Barents Sea. In section 3 we estimate a recruitment function for cod, we study the relationship between the amount of cod and individual weights, and we use catch-age analysis to estimate fishing selectivities, natural variation, and natural mortality. A summary of all parameters is given in section 3.4. Time-series from simulations with the resulting model are used to estimate an aggregate model for cod. Section 4 presents the optimization method. In section 5 we simplify the model to such an extent that it produces the yield-per-recruit policy for this problem. Using the proportional harvesting strategy implied by the yield- per-recruit method, we examine how this strategy changes as the model is made more and more complex. Using the full model, we search for more correct harvesting policies in section 6. It turns out that only marginal improvements result as non-linearities are introduced and as the age distribution is considered. In section 7 policies for the cohort model are compared to policies for aggregate models. The difference between the policies is highly sensitive to assumptions about elastic prices and costs. A comparison to earlier results for cohort models is also made. Section 8 presents near-to-optimal harvesting strategies when error is introduced in biomass assessments. In section 9 the sensitivity of the harvesting policy to uncertainty in the estimate of natural mortality is investigated. Section 10 concludes.

(7)

2. A MULTICOHORT MODEL

We want to maximise the expected net present value V for the cod fishery in the Barents Sea:

V E p p H H c c c e

e e c e

t

t t

t t t

¦

f ^U

^{

⁽ ⁰ ¹ ⁾ ⁽ ⁰ ⁽ ¹ ⁰⁾⁽ ^{) )}^D

^}

0 0

2 0 (1)

The discount factor is denoted by U, etis the applied fishing effort, and e₀ denotes constant fishing capacity. Price of fish is given as a linear function of harvest with parameters p0 and p1. Unit variable costs equal c0 at zero effort, and they equal c1 when effort equals capacity e0. Increasing marginal costs are ensured by assuming D>0. The per unit leasing cost of capacity is c2. We explicitly avoid maximizing a social welfare function for the fishing nations. Most of the harvest is exported and the domestic prices reflect export prices.

The Barents Sea cod fishery is managed by yearly quotas. Therefore we focus on a strategy for quota setting. The quota or harvest H_t will be found as a function of total biomass or biomass measures for different age classes. Capacity e0 will also be considered a policy variable when appropriate. The selectivity of the fishing gear will be kept constant and is not considered a policy variable in this study.

In the following we use capital letters to denote fish in biomass terms (million tonnes), while lower case letters are used to denote numbers (billion fish). Also note that the choice of quota or harvest as the decision variable has implication for the way in which the model is formulated. If fishing mortality had been in focus, a more traditional approach would have been chosen.

Total effort e_t is derived from harvest using a standard instantaneous catch per unit effort relationship, h ( /e X X₀)^E, where X X₀ is the biomass for which effort is defined equal to harvest:

e X X

X

X H

t X

t t t

0 0

1

0

1 1

{( ) ^E ( ) ^E}/ ( E) (2)

The expression is found by solving the catch per unit effort relationship for e and by integrating over X from X_t H_t to X_t, see Clark (1985). Ideally, there should have

(8)

been a stochastic variable in this equation since the effort needed to catch a given quota is likely to vary from year to year. However, a test shows that such an extra random variable is of little importance for the optimal harvesting policy. The catch per unit effort relationship should also be expected to be related to the selectivity of the gear. This is no problem since we will keep the selectivity constant. However, in future studies it is important to reconsider the catch per unit effort relationship if selectivity is allowed to vary.

A cohort model is used to describe the biology of cod. Several papers have shown that age distributions matter, for the case of stochastic models see Mendelssohn (1978) and Spulber (1983). These papers however, limit themselves to rather simple models since according to Mendelssohn: “The large increase in analytic complexity caused by the addition of even the simplest interaction term is cause for both consternation and chal- lenge.” Since our method allows for greater model complexity, we introduce “interaction” terms to capture vital feedback mechanisms, e.g. recruitment and weight.

The resource dynamics are given by the following equations:

x₃_,_t S_t₃exp(r₀ r S₁ _t₃ r J₂ _t v_{r t}_, )r_s (3) x_i 1_,_t 1 x_{i t}_, exp(m v_i _{i t}_, )h_{i t}_, exp(m v_i _{i t}_, / ) 2 i=3,4,...,13 (4) x₁₅_,_t₁ x₁₄_,_texp(m v_{14 14}_,_t)h₁₄_,_t exp(m v_{14 14}_,_t / )2

x₁₅_,_texp(m v_{15 15}_,_t)h₁₅_,_texp(m v_{15 15}_,_t / )2 (5)

where x_3,_t represents recruitment of three year old cod. S_t₃ is the biomass of the spawning stock at the appropriate point in time, Y_t is a measure of cannibalistic cod juveniles, weighted according to suitability, v_{r t}_, ~ N( ,0V represents random recruit-_r) ment variability and r_s is the maximum recruitment. Yearclass harvest is denoted by h_{i t}_, , m_i is the natural mortality for yearclass i, and v_{i t}_, ~ N( ,1V represents random_m) variations in natural mortality. Even though one might expect natural mortalities for age classes to be influenced by some of the same environmental forces, we disregard this possibility here and assume independence. Suitability matrices (stomach content analysis) indicate that there is a certain cannibalism on three year old cod. We ignore this direct relationship since the bulk of cannibalism is supposed to be captured by the recruitment function.

Due to the choice of total harvest as the decision variable, it is most practical to use age- class harvests and not fishing mortalities in these equations. To facilitate this, we have

(9)

made use of Pope’s approximation, i.e. harvest is assumed to take place in the middle of the year. This approximation is thought to yield good results for cod². Equation 5 shows that the survivors of age class 15 are re-entered into this age class. This is not a perfect way to represent fish older than 15 years of age since fish weight stays constant. With normal fishing activity, however, there are very few fish in the upper age classes such that this approximation should be of little concern.

The spawning stock biomass is given by ogives o_i, age class weights w_{i t}_, , and age class numbers x_{i t}_, .

S_t o w x_i _{i t} _{i t}

¦

i¹⁵₃ ^, ^, ⁽⁶⁾

The total biomass of harvestable fish (3 years and older) is

X_t w x_{i t} _{i t}

¦

i ^, ^, 3 15

(7)

Juveniles J_t represent a weighted average of biomass in lower age classes. The weights s_i reflect suitability of pre-recruitment cod for these age groups.

J_t s w x_i _{i t} _{i t}

¦

i ^, ^, 4

15

(8)

The harvest from each age class is derived from the total harvest H_t, the policy variable.

h H q x

q x w

i t

t i i t

i i t i t i

,

, ,

¦

¹²³ ⁽⁹⁾

Here q_i represents the selectivity of the fishing gear. One can easily see that the indivi- dual harvests in biomass terms ( h w_{i t}_, _{i t}_, ) sum to H_t. However, when the selectivity var- ies over age classes, Equation 7 does not ensure that harvests will be less than popula- tion numbers in each age class. This problem is easily seen in the case that H_t B_t. In this case all age classes should have been harvested completely. The problem is caused by discretization in time. In a continuous world the q_i’s could stay constant while the population numbers would gradually decrease. In turn the declining populations numbers would serve to limit harvests to what is available. Fortunately this is only a problem

2 Personal communication with Bjarte Bogstad at the Institute of Marine Research, IMR, Bergen.

(10)

when H_t is close to B_t. Since any reasonable harvesting policy will keep a good dis- tance, the weakness of the formulation usually presents no problem³.

Based on observed patterns seen in VPA data, harvesting selectivities are given by a logistic function:

q_i e^u^{i t}^, / {1 q_h /i ^q^e} i=3,4,...,15 (10) For older age classes q_i tends towards 1.0, q_hdenotes the age at which q_i equals 0.5, and the exponent q_e influences the steepness of the function. Selectivities are also influenced by natural variation, u_{q t}_, ~ N( ,0 V ._q)

Nearly all model studies we have come across ignore the effect of intraspecies competi- tion in terms of the effect of own stock biomass on own weight. One exception is Ault and Olson (1996). We assume that the weight of each age class is given by a reference weight for this age class w_{i ,0} times a weight index w_t:

w_{i t}_, w w_i_,₀ _t (11)

where the weight index depends on the cod biomass and a random variable v_{w t}_, ~ N( ,0V :_w)

w_t B_t^M^/(¹^M⁾e^v^{w t}^, w_s (12)

To avoid simultaneous equations, we replace cod biomass X_t by an approximation based on standard weights, B_t

¦

_i¹⁵³w x_i^,⁰ _{i t}^, . Weight is assumed to stay below an upper limit w_s in case cod biomass becomes very low. Using one common weight index implies that we ignore possible differences between age classes. By ignoring time lags, we simplify the model in that we do not introduce new state variables.

3 If measurement errors are introduced the problem reappears since measurement errors could lead to very high harvests relative to total biomasses. Later in this report, a formulation which increases low q_i’s as H_t approaches B_t will be used to alleviate the problem. Another solution is to use effort rather than harvest as the decision variable. At low stocks an overoptimistic effort would then lead to limited harvests due to decreasing catch per unit effort. However, we keep harvest as the decision variable, since this is the instrument the government has chosen to use in terms of quotas.

(11)

3. PARAMETER ESTIMATES

The model is meant to be representative of the Barents Sea cod fisheries. In the following we make an effort to find realistic parameters for the model.

3.1. Recruitment

Recruitment is described by Equation 3. Taking logarithms we get the following linear form:

ln(x₃_,_t)ln(S_t₃) r₀ r S₁ _t₃r J₂ _t v_{r t}_, (13) Data are shown in Table A1 of the appendix. The OLS-results are shown in Table 1.

The upper part of the table shows results for case where the juvenile biomass is based on data about suitability derived from stomach content analysis. Next we exclude some of the predating age classes to see if we get sharper results. We end up by, and report, the case where only two of the dominant age classes are included in the definition of the juveniles, J_t^*.

When using weights based on suitability matrix data, all parameters come out with the expected signs. The F-value signals a highly significant model. The effect of cannibalism, r2, is not significant at the 5 percent level. When focusing on the age classes being 3 and 4 years older than the recruits, the effect of cannibalism turns out to be highly significant. The higher parameter value is largely explained by the lower average of J_t^* than of J_t. A corresponding reduction in the average of J_t (multiplying by 0.48) leads to a parameter estimate of -0.54 for r2. Hence, the estimates of the two models are quite close. We will use the definition of juveniles based on suitability data, since the suitability data represent an informed prior in a Bayesian sense.

The similarity of the r1 and r2 estimates do not imply that cannibalism is just as important as the spawning stock biomass. A model without the variable J_t produces an R² of 0.28 as compared to 0.31 when J_t is included. The moderate importance of cannibalism in this model might also reflect weaknesses of the model. If we had modelled predation of juveniles on all age classes including the pre-recruitment zero, first, and second year cohorts, the effect might have turned out differently. It is beyond the scope of this study to do this.

(12)

Table 1: Recruitment model.

J_t ^R² Std.Error Residual F-value model Observations

0.31 0.64 10.5 46

Parameter Estimate Std Error t ratio p-value

r0 0.85 0.24 3.6 0.001

r1 -0.70 0.22 -3.1 0.003

r2 -0.25 0.18 -1.4 0.16

J_t^* ^R² Std.Error Residual F-value model Observations

0.40 0.59 11.0 46

Parameter Estimate Std Error t ratio p-value

r0 0.96 0.19 5.2 0.000

r1 -0.63 0.20 -3.2 0.003

r2 -0.77 0.25 -3.1 0.004

We also ignore an observed autocorrelation of the residuals. The Durbin-Watson observator is 1.04. A more detailed study of autocorrelations show the following significant effects. There are positive correlations with the first and the sixth lag, and negative correlations with the third and the fourth. In other words, there are indications of fluctuations with a period of around 5 years. We will not model this cyclical tendency in this study, however, we acknowledge its existence for further studies.

Also note that recruitment peaks at a spawning stock of 1.4 million tonnes; an obser- vation which is not sensitive to the value of J_t. This is a rather high spawning stock and this finding is likely to be important for ensuing optimization results.

3.2. Weight (and length)

Next we attempt to estimate the effect of cod biomass on the weight index in Equation 12. Due to limited amounts of weight data for cod, we also make use of somewhat longer time-series for cod length. The results will be compared. In a yearclass model, relative weight growth can be expressed as follows:

w w

g f

i t i t

i

1, 1 ,

( ) (14)

Here g_i is the average weight growth for age class i over time. For now, we do not specify which variables influence relative weight growth, we only indicate that growth is determined by a function f_i( ) for each age class i. Our first task is to see if the relative growth rates of the different age classes are correlated. If they are, we could simplify the model greatly by replacing all f_i( ) ’s by a common function f ( ) .

(13)

Weight data for this test are shown in Table A2, corresponding length data are shown in Table A3. Typically, there are fewer observations of younger and older fish than of medium aged fish. Thus measurements are less reliable for high and low age classes than for the medium ones. Also note that older data are likely to be less reliable than newer data. For this reason we concentrate on data from 1973 to the present, a period for which we also have data on capelin abundance, see Table A5 for time-series on biomasses of cod and capelin and of temperature. The older data will be used in a supplementary test. Relative growth rates for weight and length⁴ are shown in Figures 1 and 2. Tables 2 and 3 show correlations between age classes for respectively relative weight and relative length growth.

Figures and tables give a clear message about correlation between age classes for both relative weight and length growth. The figures show that age groups 3 to 6 largely follow the same pattern for both weight and length. Concerning correlations, all p-values are less than 0.02 for weight and 0.06 for length within these age groups; and usually the p-values are much smaller. Age groups 1 and 2 are also largely correlated with age groups 3 to 6. In most comparisons, age group 7 is not significantly correlated with the rest. A possible explanation is poor quality of data. However, a significant correlation between age groups 6 and 7 for weight growth, indicates that the nature of growth changes in higher age classes. Since we do not have data to investigate this possibility further, and since the age classes below 7 have by far the largest population numbers, we conclude that relative growth in weight and length could be assumed to be determined by the same factors for all age classes.

Relative weight growth

0 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4

1985 1987 1989 1991 1993 1995

1 2 3 4 5 6 7

Figure 1: Relative weight growth for cod in different age classes.

4 Growth rates for length for the period 1974 to 1988 are based on somewhat refined growth data in

Jørgensen (1992), that is, compared to the ones in our Table A3.

(14)

Relative length growth

0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80

1974 1977 1980 1983 1986 1989 1992 1995

2 3 4 5 6 7

Figure 2: Relative length growth for cod in different age classes.

A test for correlations does not reveal differences with respect to the size of variations in relative growth between age classes. Table 4 shows standard deviations for the time- series of relative growth for weight and length. There seems to be a tendency towards decreasing variation in relative weight growth with age, while variation increases with age for length growth. However, the differences are small and we ignore them⁵.

Table 2: Correlations between relative weight growth for different age classes.

2 3 4 5 6 7

1 U 0.82 0.66 0.76 0.59 0.57 0.41

p-value 0.004 0.038 0.011 0.070 0.087 0.246

N 10 10 10 10 10 10

2 U 0.76 0.78 0.76 0.69 0.35

p-value 0.011 0.008 0.010 0.028 0.319

N 10 10 10 10 10

3 U 0.88 0.91 0.74 0.55

p-value 0.000 0.000 0.006 0.064

N 12 12 12 12

4 U 0.89 0.67 0.31

p-value 0.000 0.018 0.333

N 12 12 12

5 U 0.81 0.52

p-value 0.001 0.086

N 12 12

6 U 0.75

p-value 0.005

N 12

5 The effect of age on variation could be implemented by a modification of Equation 14:

w w

g f

i t i t

i

ai

1 1

, ,

( ( ) )

D , where Dai would decline with age.

(15)

Table 3: Correlations between relative length growth for different age classes.

3 4 5 6 7

2 U 0.40 0.63 0.66 0.45 0.36

p-value 0.087 0.004 0.002 0.054 0.133

N 19 19 19 19 19

3 U 0.73 0.75 0.40 0.16

p-value 0.000 0.000 0.061 0.497

N 23 23 23 20

4 U 0.83 0.57 0.24

p-value 0.000 0.005 0.318

N 23 23 20

5 U 0.70 0.20

p-value 0.000 0.389

N 23 20

6 U 0.16

p-value 0.505

N 20

Table 4: Standard deviations for relative growth in weight and length over age classes.

2 3 4 5 6 7

Weight 0.50 0.44 0.44 0.43 0.36 0.36

Length 0.16 0.21 0.25 0.26 0.31 0.45

Finally we note that relative growth in weight and relative growth in length is highly correlated within each of the age classes 2 to 5, see Table 5. Thus, the longer time-series for length data could be used to support findings from the shorter time-series for weight.

Table 5: Correlations between relative length and weight growth for different age classes

2 3 4 5 6 7

U 0.91 0.80 0.94 0.89 0.49 0.34

p-value 0.000 0.002 0.000 0.000 0.102 0.284

N 10 12 12 12 12 12

Now we turn to the dynamics of weight growth. From Equation 14 we see that weight in age class i+1 is influenced by two factors. First there is the weight that is inherited from the lower age class, w_{i t}_, . It takes time before a change in a lower age class reaches higher age classes. The delay is longer, the higher age class one considers. Second there is the exogenous factor f ( ) influencing relative growth in each and every age class.

This effect is direct with no lag. The combined effect should be expected to be delayed, however, with a time lag that is shorter than the one caused by the first factor in isola- tion. The time lag should be expected to increase with age.

Before we turn to the data, note that the model in Equation 14 is not entirely consistent with much literature on weight growth, particularly the literature on recovery or catch-

(16)

up after starvation. Blasco et al. (1992) and Russel and Wootton (1992) find that growth increases above normal after periods of starvation such that in the long run weight is not reduced by the starvation incident. If this is correct for Arctic cod, we should expect to see that the exogenous factor dominates in the recovery phase, i.e.

higher age classes should gain weight in phase with lower age classes.

However, the literature is not unanimous on this point. Karlsen et al. (1995) report a weight loss also in the longer-term. Since cod seems to grow faster in aquacultures than in the Barents Sea⁶, it could be that wild cod is more or less permanently starved.

Hence, after a period of severe starvation, more ample food supplies in the ensuing period might still not be sufficient for the weight to catch-up. If this effect dominates, we should expect delayed weight increase in higher weight classes compared to lower age classes.

A first investigation of possible delays can be made by inspecting Figures 3 and 4. For both weight and length we see that variations in upper age classes (higher weights and longer lengths) are delayed compared to variations in lower age classes (indicated by dashed lines connecting troughs and peaks). If higher age classes were rapidly catching up after periods with starvation, we should have seen a simultaneous recovery of all age classes. Alternatively, the delay could be argued to be caused by a delayed food recovery (prey) for higher age classes. However, this is not very likely since relative growth rates are highly correlated between the different age classes, see above. The patterns that emerge in Figures 3 and 4, are consistent with patterns produced by a simulation model based on Equation 14. Figure 5 shows simulation results when g_i=10.0 for all i, and f_i( ) is the same sine wave varying from 0.6 to 1.4 for all age classes with a period of 6 years. Note that variation diminishes at high ages. In fact for the highest age class, the fluctuations disappear completely in the simulation model. The inherited and the direct effects come out of phase with each other and tend to cancel each other rather than strengthening each other as they do for lower age classes. Data in Figures 3 and 4 seem to show the same tendency.

Thus far we have found that relative growth rates are similar over most age classes, and that weight and length in higher age classes should be expected to be delayed relative to the factors that influence growth. We have not yet discussed which factors contribute to growth. In our one-species model, concern is with the effect of the stock size of cod

6 Growth rates around twice that of arctic cod has been obtained in aquacultures. Correction for different temperatures has not been made. Source: personal communication with Jens Kristian Holm, Austevoll.

(17)

Cod weight

0.00 1.00 2.00 3.00 4.00 5.00 6.00

1985 1987 1989 1991 1993 1995 1997

Figure 3: Observed weights for age classes 3 to 8.

Cod length

0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0

1973 1976 1979 1982 1985 1988 1991 1994 1997

Figure 4: Observed lengths for age classes 2 to 8.

Simulation

0 10 20 30 40 50 60 70

0 2 4 6 8 10 12 14 16 18 20

Figure 5: Simulated lengths using Equation 14

(18)

X_t. To enable more precise estimates we will also include temperature T_t and a key prey for cod, capelin L_t.

Since we at this stage seek a simple model of weight, we do not estimate Equation 14 for each age class. Rather we formulate a direct relationship between explanatory variables and weight while allowing for a time delay. This relationship is estimated for different age classes for both weight and length. Below we will use these estimates to formulate the weight equation to be used when searching for optimal policies. The weight equation is:

w_{i t}_,₁ w X_{i t}^M_,¹ _t^M²L T e e^M_t³ _t^M⁴ ^M⁵ ^H^{i t}^, (15)

Table 6: Results of estimation using weight data (p-values for parameters in parentheses).

age M1 M2 M3 M4 M5 R² F p-value Std.err.

2 0.40 -0.91 -0.04 1.38 -3.31 0.96 32 0.001 0.14

(0.337) (0.009) (0.744) (0.096) (0.102)

3 0.12 -0.57 0.18 1.35 -2.88 0.94 29.1 0.000 0.14

(0.419) (0.001) (0.011) (0.044) (0.013)

4 0.24 -0.19 0.15 1.21 -2.03 0.88 13.37 0.002 0.17

(0.218) (0.245) (0.025) (0.110) (0.095)

5 0.54 -0.15 0.08 1.11 -1.59 0.84 8.96 0.007 0.15

(0.017) (0.287) (0.115) (0.129) (0.177)

6 0.68 -0.11 0.06 0.97 -1.30 0.88 12.38 0.003 0.11

(0.002) (0.272) (0.078) (0.083) (0.166)

7 0.80 -0.05 0.03 0.79 -1.02 0.72 4.539 0.040 0.11

(0.010) (0.608) (0.329) (0.171) (0.344)

A multiplicative form is chosen to capture the non-linear, saturating effect of capelin on cod weight⁷. The time delay is captured by M1. From a linearized version of the model the lag time can be calculated as W 1/ ln(M1) . The parameters M², M³, and M⁴ repre- sent elasticities, i.e. they measure of the relative change in w_{i t}_,₁ per unit of relative change in the explanatory variables⁸. To find the long-term elasticities for the explanatory variables, each elasticity should be divided by⁹: (1-M¹). Logarithms are taken when estimating. Results for the case with weights are shown in Table 6, while the results for length (using the same equation) are shown in Table 7. The corresponding parameter estimates for age classes 3 to 6 are also shown in Figures 6 and 7.

7 The multiplicative model produces much more precise estimates of M1 than a corresponding linear one. A weakness of the nonlinear model is that zero capelin implies zero cod weight. This is not realistic since cod benefits from other preys as well. However, even though there are some very low observations of capelin in the data, the nonlinear model does not produce much different results than a linear one for capelin.

8 Note that the elasticity for cod biomass is written differently than in the model, Equation 12. This is because in Equation 12 it is a weighted measure of population numbers that enter the equation, while in the estimation we apply a biomass measure.

9 This is easily seen when one sets w_{i t}_,₁ w_{i t}_, .

(19)

Table 7: Results of estimation using length data (p-values for parameters in parentheses).

age M1 M2 M3 M4 M5 R² F p-value Std.err.

2 0.39 -0.114 0.028 -0.01 1.96 0.68 10.23 0.000 0.08

(0.255) (0.075) (0.245) (0.929) (0.068)

3 0.27 -0.131 0.041 0.12 2.42 0.78 16.71 0.000 0.06

(0.073) (0.000) (0.001) (0.266) (0.000)

4 0.46 -0.090 0.030 0.29 1.58 0.83 23.75 0.000 0.04

(0.000) (0.000) (0.000) (0.001) (0.000)

5 0.63 -0.053 0.024 0.17 1.23 0.74 13.16 0.000 0.04

(0.000) (0.016) (0.001) (0.033) (0.024)

6 0.75 -0.034 0.018 0.13 0.84 0.69 10.357 0.000 0.04

(0.000) (0.061) (0.003) (0.055) (0.174)

7 0.43 -0.035 0.005 0.07 2.35 0.25 1.566 0.224 0.06

(0.050) (0.201) (0.611) (0.498) (0.016)

The length model produces highly significant estimates for most parameters, the major exception is age class 7. All parameters with limited p-values have the expected signs. F- values peak for age class 4 and declines to both sides. This is what one should expect from the simulations in Figure 5 since the weight variation tends to peak at medium age.

Important to note here is that the lack of effect for age class 7 should not be taken as a proof of lacking importance of cod, capelin, and temperature for this age class. A more likely explanation is the balancing effects of inherited and the direct weight growth.

A Shapiro-Wilk test shows that we cannot discard the hypothesis of normally distributed residuals. A test of partial autocorrelations of up to 16 lags for the 6 time-series, show that only 4 out of 96 cases come out as significantly autocorrelated. This is slightly less than what should be expected in case of no autocorrelation, using a significance level of 5 percent.

For the weight regressions, the signs are also as expected, F and p-values for the entire models are quite good, while most parameters are not significant at the 5 percent level.

A comparison of the obtained parameter values for the weight and length models show considerable consistency however. We use this consistency to argue that if the obtained estimates for the length model are reliable, the estimates of the weight model are also reliable. The higher p-values for the weight model are likely to stem from the shorter time-series. Residuals are accepted as normally distributed and there are no signs of autocorrelation.

We note the following interesting details seen in the figures. First, the delay time increases with age for both weight and length, as postulated above. The delay time is slightly longer for length than for weight. This is consistent with the fact that the

(20)

intervals between age classes increases for weight, while they are nearly constant for length, see Figures 3 and 4. The increasing intervals imply that the direct effect of the explanatory variables get relatively more weight compared to the inherited effect from lower age classes.

Elasticities and time lag for weight

0.00 1.00 2.00 3.00

3 4 5 6

0 1 2 3 4

Cod Capelin Temp.

Delay Delay

Elasticities

Figure 6: Results of weight estimation for the age classes 3 to 6, absolute values of long-term elasticities and time delay (years).

Elasticities and time lag for length

0.00 0.25 0.50 0.75 1.00

3 4 5 6

0 1 2 3 4

Cod Capelin Temp.

Delay Delay

Elasticities

Figure 7: Results of length estimation for the age classes 3 to 6, absolute values of long-term elasticities and time delay (years).

Second, the long-term elasticities for weight in age classes 3 to 6 are approximately 3 times as strong as the ones for length (notice the different scales on the elasticity axes).

This is consistent with the fact that weight depends on volume, which is related to the cube of length. To the extent that proportions are maintained during growth, the factor 3 should be expected. This is easily seen by raising equation 15 to the power of 3.

(21)

Third, the effects are quite strong. To take age class 4 as an example: If the stock of cod increases from 1 to 2 million tonnes, the long-term effect on weight is a reduction of 16 percent. An increase in capelin from 1 to 6 million tonnes leads to a weight increase of 43 percent. Finally, a temperature increase from 5 to 6 ^oC, leads to a weight increase of 34 percent.

The statistics are quite good, particularly if each regression is seen as an independent experiment. However, it is important to keep in mind that we have only repeated observations of effects for one and the same set of explanatory variables. On the positive side, the smooth development in parameters and statistical indicators over age classes indicate that the measurements of weight and length are quite reliable¹⁰. On the negative side, our results could be influenced by autocorrelation in weights, in lengths, and in explanatory variables. Hence, while we find precise estimates, the exact parameters could change somewhat as new data become available.

The more incomplete and probably less reliable data for the early years in Table A3, can be used to assess the ability of the model to predict the situation in the 1950s when the cod biomass was considerably higher than in the period used for the estimation. The average cod stock for the period 1953 to 1958 was 3.42 million tonnes, as compared to 1.69 million tonnes for the period from 1973 to 1996. Average lengths for the earlier period were respectively 2.2, 5.6, 5.5, and 5.4 percent lower than in the later period for the age classes 4 to 7. Using the elasticity M2 / (1M1) for the long-term effect of cod on length, the model predicts lengths to be reduced by respectively 11.1, 9.6, 9.1, and 4.2 percent for the same age classes. Averaging over the four age classes, lengths are measured to be 4.7 percent lower in the earlier period, while our models predict them to be 8.5 percent lower. Temperatures are not important since they were nearly the same in the two periods. The average stock of capelin was 3.4 million tonnes in the later period.

What it was during the 1950s is not known. According to the estimated model, likely average long-term effects on lengths vary from -4.2 to 3.3 percent for capelin stocks in the range from 1 to 8 million tonnes. Adding the effect of cod to the likely range for the effect of capelin, we find a likely range for the reduction in length spanning from 5.2 to 12.7. The entire range denotes a stronger effect than what was measured, 4.7 percent.

The above test indicates that the sensitivity to cod is overestimated. To find an estimate of the degree of overestimation, we perform the following test. First we repeat the estimation of Equation 15 over the period from 1973 to 1996, however, excluding the explanatory variable capelin. Next, the same estimation is performed over the available

10 This conclusion requires that measurement errors are not strongly correlated.

(22)

data points from 1953 to 1996 (n=37). Since there are no older data for age classes two and three, they are not test. Table 8 shows the estimated elasticities for cod for the two tests. F-values for the entire models indicate that it is only the short time-span model for age class 7 that does not produce a p-value below 0.05. Except for age class 7, all F- values are lower than for the corresponding model with capelin, Table 7.

Table 8: Length elasticities with respect to cod (p-values for parameters in parentheses).

age M2 short F short M2 long F long ratio M2’s:

long/short

ratio M2’s:

cap./no cap.

4 -0.079 11.7 -0.042 7.7 0.53 1.14

(0.01) (0.06)

5 -0.044 7.7 -0.034 6.3 0.79 1.13

(0.10) (0.09)

6 -0.025 6.3 -0.015 5.2 0.62 1.36

(0.25) (0.34)

7 -0.033 2.1 -0.027 2.9 0.83 1.06

(0.22) (0.14)

Only half the M2 parameters are significant at the 10 percent level. However, the ratio between M2’s for models using long and short time-series data, do not vary much. In all cases the longer time-series imply lower sensitivity to cod by a factor varying from 0.53 to 0.83. On average for the four age classes, cod elasticities for the longer time-series are 31 percent lower. Comparing elasticities for the model with short time-series to the parameters for the model with capelin in Table 7, shows that cod elasticities in the model with capelin data are on average 17 percent higher. The cod elasticities to be used in the following analysis should reflect both the influence of capelin and of the longer time- series data. Hence we take the cod elasticities in Tables 6 and 7 as our starting point, and we adjust them downwards by 31 percent to reflect the older data.

The elasticity M to be used in the static relationship in Equation 12, is set equal to the average long-term weight elasticity, M2 / (1M1), for age classes 4 to 7, which equals 0 30. , times the adjustment from older data, 0.69. Hence M=-0.2. The long-term elasticity is chosen because this is the elasticity of importance for the optimal solution. By simulating equations 11 and 12 from 1985 to 1996, we estimate the standard deviation of the error term for age classes 4 to 7. Standard deviations vary between 0.21 and 0.47.

The standard deviation for v_{w t}_, is set equal to the average, i.e. Vw=0.34. To be correct, the index should have been multiplied by a constant, however, this constant comes out close to 1.0 in the simulation.

Using B_t (the sum of population numbers multiplied by average weights for each age class) rather than X_t, does not seem to be a problem. Regressions using population numbers rather than biomass for cod lead to quite similar results. Hence, it is the under-

(23)

lying variation in population numbers and not the smaller variation in weight that matters for the biomass measure.

3.3. Mortality, fishing selectivity and environmental variation

When it comes to estimation, mortality, fishing selectivity and environmental variation are not independent measures. Hence, all these measures should be estimated simul- taneously. We make an effort in this direction using a variant of “catch-age” analysis, see e.g. Deriso et al. (1985) and Megrey (1989). Using this technique, parameters determining mortality, fishing selectivity, and natural variation are found along with parameters denoting initial stock sizes by year class and parameters denoting each year’s recruitment.

We will use a maximum likelihood method¹¹. Harvests from the individual age classes and from different years are assumed to be lognormally and independently distributed.

The probabilistic model for each observations is:

> @

p y

y

y h

i t

i t i t

( _, | ) exp ln( ) ln( )

,

, ,

II ®

¯

½¾

¿ 1

2

1 2

2 2

2

SV V i=3,4,...,15, t=1,2,..,52 (16)

where y_{i t}_, is measured harvest, h_{i t}_, is predicted harvest, i denotes age class, t denotes the year, and II { , , , , , }m q q_h _e V r x_t _i,1 is a vector of model parameters to be estimated, where m denotes natural mortality, q_h and q_e determine fishing selectivities, V is a constant standard error, r_t is yearly recruitment, and x_{i ,1} represent initial conditions (where r1= x_{3 1}_, ).

The probability of all the data y is given by the product of the individual probabilities.

Taking the logarithm of this product, we get the log-likelihood of the sample:

> @

l y_{i t} y h

i t

i t i t

i t

( | )yII ^{13 52} ¹₂ln(² ²)

¦ ¦

ln( _, )₂¹

¦ ¦

ln( _, )ln( )_,

3 15

1 52

2

2 3

15

1 52

SV V (17)

Maximizing l( | )yII for II yields the maximum likelihood estimates of the parameters. To find predicted values of h_{i t}_, we simulated the equations 3, 4, 5, and 9. Hence we adapt the catch-age analysis to our choice of model. This implies that we consider total yearly

11 A least square method has also been tested with almost identical results. This indicates that it is not very important whether the criterion measures relative error or absolute error. Using absolute errors, only abundant and younger age classes will be of any importance for the criterion.

(24)

harvests as an exogenous variable, and that the model predicts harvests for each age class. The equations are repeated below with the proper adaptations to the estimation problem:

,

, ,

h h q x

i t q x

t i i t i i t

¦

¹⁵i ³ ⁽¹⁸⁾

x_3,_t r_t (19)

x_i 1_,_t 1 x_{i t}_, exp( m) h exp(_{i t}_, m/ ) i=3,...,132 (20) x₁₅_,_t₁ x₁₄_,_texp( m) h₁₄_,_texp(m/ )2 ^x14_,t exp( ^m) ^h exp(12_,t ^m/ )2 (21) The simulations start out with the initial values x_{i ,1}. Harvests h_{i t}_, are predicted based on stock sizes x_{i t}_, , selectivities q_i and total yearly harvests in numbers h_t

¦

_i¹⁵³y_{i t}^, ^{. This} equation differs from Equation 9 in that total yearly harvests are measured in numbers rather than in biomass. Consistent with this, weights do not enter in the denominator. By inspection we see that

h_{i t}, h

i t

¦

¹⁵³ . The reason for doing this is that the explanatory or exogenous variable total yearly harvest is measured in numbers and not in biomass, which was the case for the decision variable in the optimization model. If total yearly harvested biomass H_t (based on weights and harvests in numbers) had been used as the explanatory variable in the estimation, errors in predicted stocks x_{i t}_, would cause total yearly harvests in numbers to deviate from the observed numbers h_t. The estimates we obtain will also be the correct ones for the optimization model.

To simplify the estimation, we estimate only one constant mortality rate m over all age classes. Similarly, we operate with a fixed variance for the total prediction error over age classes and years. Fishing selectivities depend on two parameters q_h and q_e:

q_i 1/ {1 q_h /i ^q^e} i=3,4,...,15 (22) The simulation is carried out in a straightforward manner, where no attempt is made to update stock estimates and harvest estimates based on recent measurements of year class harvests (ballistic simulation). This method is ideal when there is only measurement error (this point is easily seen from the Kalman filter equations). In our case there is also process error (in mortality and in fishing selectivity). Thus the model is not ideal. To get an indication of its goodness (bias in estimates), we will produce synthetic data by a