Random-walk dynamics of exploited fish populations

Abstract

Introduction

Fluctuations in marine fish landings have been a concern for many years (Botsford et al., 1997; Pauly et al., 2002; Halley and Stergiou, 2005; Niwa, 2006a). The variation in supply (the stock size) is a central issue when fishery managers decide how best to regulate fishing effort. A fishery scientist's standard approach to trying to ensure sustainable harvesting is to identify the location of the stock-recruitment curve in relation to the replacement line defining the recruitment level needed to replace the spawning stock in future (Ricker, 1954; Beverton and Holt, 1957; Ulltang, 1980; Caddy and Sharp, 1986; Cook et al., 1997; Hjermann et al., 2004). The place where the two lines intersect is the equilibrium point to which the population is attracted. The resulting theory on sustainably harvesting fish stocks would be mathematically convenient for resource managers, but it fails abysmally to fit historical data on abundance (Shepherd and Cushing, 1990; Dixon et al., 1999). The fitted curves have been plagued by the extreme scatter of data points, and there is always doubt about the precise level of sustainability (Schnute and Richards, 2001; Longhurst, 2006). The term “equilibrium” may even be an empty concept. Despite substantial effort directed towards elucidating the population dynamics of exploited fish stocks, declines and recoveries in abundance remain an ecological enigma (Myers et al., 1995; Hutchings, 2000; Sibly et al., 2005). Characterizing temporal patterns in population size is a major challenge for ecological time-series analysis, especially attempting to show evidence or suggestions of the source of the variability.

Material and methods

Time-series data

Spawning-stock biomass (SSB) data are generally derived from sequential population analysis, which can lead to systematic error in the abundance computed. Because time-series results are the product of a quite complex and potentially biased estimation process (results for successive years are inferred on the basis of landed catches), they may be subject to change when data are tuned in various ways. Nevertheless, it is generally considered by scientists that there is good correspondence between time-series data and actual levels of population abundance.

Here, I extract information on the temporal variability in SSB from time-series data on commercial fish stocks in the North Atlantic. The data are derived from the 2005 working group reports of the International Council for the Exploration of the Sea (ICES), prepared to support advice on living resources and their harvesting provided by the ICES Advisory Committee on Fishery Management (ACFM). The information can be found at <http://www.ices.dk/advice/icesadvice.asp>. I select from the database those stocks for which the annual time-series of SSB encompassed at least 30 y, in all 27 marine fish stocks addressed by ICES. The stocks and the distributions analysed are listed in Table 1, but for completeness here, the species I investigated were cod (Gadus morhua), Greenland halibut (Reinhardtius hippoglossoides), haddock (Melanogrammus aeglefinus), herring (Clupea harengus), plaice (Pleuronectes platessa), saithe (Pollachius virens), sole (Solea vulgaris), sprat (Sprattus sprattus), and whiting (Merlangius merlangus). The ICES fishing areas, subareas, and divisions for which these stocks are analysed are outlined on the ICES website (http://www.ices.dk/aboutus/icesareas.asp).

Population fluctuations

Here, I analyse how the biomass of exploited fish populations ebb and flow in the time-series data across all 27 stocks. Let S(t)be the SSB at time t. Temporal variation is considered by examining the natural logarithm of the ratio of successive abundances, i.e. the population growth rate

1

which approximates the relative changes [S(t + 1) – S(t)]/S(t), an annual increment of Sper population size, leading to a rate equation:

2

providing the limit of small changes in S.

Statistical tests rely on the distributions of abundance fluctuations r(t) and their serial correlations over time, in applying a test fit to a Gaussian model, a test of more time, more variation (Lawton, 1988; Inchausti and Halley, 2001, 2002), and a test of draw-downs and -ups (Johansen and Sornette, 1998, 2001) to examine time-series of population fluctuations. One would expect some stocks to grow more slowly than others, resulting in their having less resilience to population decline. To compare different stocks, population variability is standardized by dividing r(t) by its standard deviation σ_r.

Test 1: distribution of abundance fluctuations

To characterize the population process, the probability distribution of population growth rate is plotted. As the distribution of abundance S(t) of a stock through time is often approximately lognormal (MacArthur, 1960; Preston, 1962; Sugihara, 1980; Engen et al., 1998; Sæther et al., 2000; Bjørnstad and Grenfell, 2001; Halley and Inchausti, 2002; Niwa, 2006b), one would expect the distribution of r(t) to be normal (Gaussian): the ratio of two numbers, S(t + 1)/S(t), each drawn from a lognormal distribution, is also lognormally distributed. A rank-ordering technique is used to plot cumulative distribution (MacArthur, 1960; Johansen and Sornette, 2001; Halley and Inchausti, 2002; Niwa, 2006b), rank-ordering statistics and cumulative distribution plots being equivalent. Let L + 1 records be observed in a time-series for a stock (i.e. a time-series of L + 1 y). The cumulative number k_ of r ≤ 0 (the number of records in which growth rates are less than or equal to r_k_), i.e. the rank in ascending order, r₁ < r₂ < … ≤ 0, and the complementary cumulative number k₊of r ≥ 0, i.e. the rank in descending order, r₁ > r₂ > … > 0, are normalized by L. I then plot

3

showing the cumulative distributions of standardized growth rates x₋ ≤ 0 and the complementary cumulative distributions of x₊ ≥ 0. If the Gaussian model is a valid description of reality, the rank-ordering plots of standardized data [Equation (3)] for the 27 fish stocks in the North Atlantic collapse onto each other, falling onto the standard normal (complementary) cumulative distribution functions (CDFs):

4

for x ≥ 0 or x ≤ 0, where the Gaussian error function is defined by erf(x) = (2/√π) ∫ ₀^xexp( − z²)dz.

Test 2: time evolution of standard deviation

I here focus attention on the dynamics of the population process, i.e. the serial correlations of abundance fluctuations over time. To do this, I study the population time-series to probe the nature of the stochastic process underlying them, by measuring the time dependence of successive abundance changes.

It is well known that the variance of population abundance increases with the length of a time-series (Pimm and Redfearn, 1988; Keitt and Stanley, 1998; Inchausti and Halley, 2001, 2002), and that is related to the “more time more variation” effect, as noted by Lawton (1988). Therefore, it is necessary to analyse the variability in exploited population abundance, examining the breadth of population fluctuations for an arbitrary time lag. The breadth of population fluctuations for time lag Δt( = 1, 2, …),

5

( < … > denotes the sample mean for each stock) corresponds to the standard deviation of the distribution function for the Δt-year growth rate:

6

For a random-walk process in which successive changes are uncorrelated, the distribution width increases as a power-law function of time lag Δt, i.e.

7

where H = 0.5. A series of accumulations of random, uncorrelated numbers over time is a random walk (the running sum of a Gaussian white noise is a Brownian motion). Complex dynamics generate serial correlations over time, so that the standard deviation increases more rapidly (0.5 < H ≤ 1, implying persistence correlations) or more slowly (0 ≤ H < 0.5, implying antipersistence) than for a random walk: movement in one direction tends to be followed by movement in the same (H > 0.5) or opposite (H < 0.5) direction.

I test for dependence in the succession of abundance records from year to year by examining how the standardized statistic

8

increases with time lag Δt, where denotes the average across all 27 stocks.

Test 3: persistent moves of populations

Another test of whether successive abundance changes are independent (i.e. uncorrelated) is to quantify the duration and the magnitude of a cumulative decrease or increase (a draw-down or -up) of population abundance over time. A draw-down is defined as the logarithmic loss from a peak (local maximum) to the subsequent valley (local minimum). For instance, if the signs (up or down) of successive changes in population trajectory over 12 consecutive years are ++−−+−+−−−−+, the first draw-down lasts 2 y, the second draw-down lasts 1 y, and the third lasts 4 y. Under the hypothesis of uncorrelated consecutive population moves with probability p_–of observing a population decline, the probability density function (PDF) observing a cumulative decrease from the last local maximum to the next local minimum of a given duration, τyears, is

9

where p₊( = 1 − p_–) is the probability of observing a positive abundance variation, and the term p₊/p₋ensures normalization of P(τ). The probability that a draw-down lasts τconsecutive years (i.e. that a draw-down is a run of length τyears) is then an exponential function with the decay constant b = −  ln p–. The mean duration of < τ > consecutive drops in a random-walk population,

10

is calculated to give 〈τ〉 = /(1–e^−b) years for a long duration, L, of a time-series.

A draw-down (magnitude) is the cumulative decline from the last local maximum Ŝ (year t) to the next local minimum S̆(year t + τ) of abundance in the natural logarithm:

11

(τ ≥ 1). A draw-up is defined as the change (in the log-transformed value) between a local minimum and the subsequent local maximum, i.e. a draw-up is the event that follows the draw-down, and vice versa. With this definition of draw-downs and the null hypothesis that consecutive log-abundance changes are uncorrelated, the PDF to observe a draw-down of given magnitude uis

12

where p(r) is a PDF of annual rates rof population growth, giving p₋ = 1 − p₊ = ∫ _−∞⁰p(r)dr; the term p₊/ p₋ensures normalization of P(u). The Dirac delta function δ(u − ∑_{j = 1}^τr_j) ensures that the sum over τ on the right-hand side of Equation (12) covers all possible run durations τ of consecutive drops r₁, r₂, …, r_τthat sum up to a given u. The tail of the distribution of draw-downs is an exponential for general distributions p(r) when their tails decay no slower than an exponential (Johansen and Sornette, 2001). It is expected that the magnitude distribution of draw-downs (or -ups) for random-walk populations results in the exponential fall-off of a Poisson distribution, as well as a persistence-time distribution. Equation (12) is solved for an independent, identically distributed (iid) Gaussian variate fluctuating around zero, i.e. p(r) is the Gaussian PDF with mean zero and standard deviation σ_r. Replacing the delta function with a Gaussian PDF (sharply peaking at the origin with width 0.05σ_r), I compare the solution with the empirical distributions of draw-downs and -ups.

Results

Test for fit to a Gaussian model

Figure 1a groups the cumulative probability distributions of standardized growth rates [Equation (3)] of 27 fish stocks in the North Atlantic (semi-logarithmic scale plot). The negative and the positive parts are the cumulative and complementary cumulative distributions (cumulative and complementary cumulative numbers are divided by L). The solid lines show the standard normal (complementary) CDFs [Equation (4)]. The values of sample means < r > and standard deviations σ_rfor individual stocks are given in Table 1; the mean rate of population growth is at the level of noise, < r > ≈0.

Figure 1.

Open in new tabDownload slide

Fitting population time-series to a Gaussian model (27 stocks in the North Atlantic, listed in Table 1; data from ICES Advice, ACFM Documents, May and October 2005). (a) Relative number of times a given level of observed population growth rate is plotted against standardized rates r/σ_rfor each individual stock (semi-logarithmic scale plot). The solid lines show the standard normal (complementary) CDFs [Equation (4)]. (b) Error plot, reading ln[y/P_CDF(x)] of data points on (a). (c) Probability density distribution of population growth rates across all 27 stocks (semi-logarithmic scale plot). Standardized data shown in (a) are aggregated, and the solid line shows the standard normal PDF.

Figure 1b is the error plot, i.e. log-differences between relative (complementary) cumulative number y_±(≡ k_±/L) and standard normal (complementary) CDF [Equation (4)],

13

are plotted against standardized growth rates x_±(≡ r_k±/σ_r). The central part (95%) of the distribution (i.e. |r| < 2σ_r) is tolerably fitted by the Gaussian model (i.e. |log-differences| < ln 2≈0.7).

Figure 1c reflects the probability distribution of population growth rate. After binning and averaging the distribution of standardized data (shown in Figure 1a) across all 27 stocks (bin size is 0.25, in units of standard deviation), a quadratic (i.e. parabolic) behaviour on a semi-logarithmic scale is obtained. The aggregated data are well described by the standard Gaussian profile (solid line), the normal PDF with zero mean and unit variance. This result implies that each distribution shown in Figure 1a is symmetrical around r = 0, suggesting that there is no important overall trend in population abundance. Moreover, the fact that mixing distributions with means of zero and identical standard deviations (unity) results in the Gaussian with mean zero and the same standard deviation (unity) is evidence of Gaussian behaviour of each population time-series.

Test for “more time, more variation”

Figure 2(double logarithmic scale plot) shows the width of the distribution functions for the Δt-year growth rate (Δt = 1, 2, …, 20 y); the distribution for time lag Δt = 1 y corresponds to Figure 1c. The standardized statistic [Equation (8)] increases with time lag Δtas in Equation (7), with H = 0.516 ± 0.013 (the adjusted R², i.e. the coefficient of determination, is calculated, giving R² = 0.9878); the solid line shows the fit of the linear regression using log-transformed data. This result implies that successive changes in the logarithms of population abundance are not correlated, i.e. the population process is a random process. The mechanism by which the distribution width 〈r_Δt²〉^1/2evolves is nothing but the consequence of the central limit theorem; the population growth rate r(t) is an iid Gaussian variate. The stochastic process S(t) corresponds to a geometric random walk.

Figure 2.

Open in new tabDownload slide

Characterization of abundance moves in the population time-series (data from ICES–ACFM Documents of 2005), in a double logarithmic scale plot. Dots show the standard deviations of Δt-year growth rates (r_Δt/σ_r) across all 27 stocks [Equation 8], and the solid line is the fit to Equation (7).

Test for draw-downs and draw-ups

Figure 3a shows the complementary CDFs of the duration of persistent decreases (triangles) and increases ( + signs) for the population time-series across all 27 stocks. The linear plot on a semi-logarithmic scale represents the exponential form of a Poisson distribution for uncorrelated data. Providing that the draw-downs and -ups are symmetrical, the fit with

14

(the number of events that the persistent duration is ≥ τ y) is done using the natural logarithm, depicted by the solid line; the linear regression gives b = 0.62 ± 0.03 (adjusted R² = 0.9725). The result that the decay constant bis close to ln 2 = 0.693 y⁻¹is evidence against a “biased coin tossing” feature, suggesting that p₊ = p₋ = 1/2 for symmetrical distributions of population growth rates. The mean duration [Equation (10)] of consecutive drops (or rises) in a random-walk population is < τ > = 2 y. For real data aggregated across 27 stocks, the mean duration is = 2.6 y for draw-downs and = 2.1 y for draw-ups.

Figure 3.

Open in new tabDownload slide

Persistent decreases and increases in the population time-series (data from ICES–ACFM Documents of 2005). (a) Complementary CDFs of persistence time (semi-logarithmic scale plot), triangles indicating draw-downs and plus signs draw-ups, and the line depicting the exponential fit to Equation (14), if the draw-downs and -ups are symmetrical. (b) The number of times a given level of draw-downs (triangles) and draw-ups (plus signs) are observed, i.e. complementary CDFs of persistent moves (absolute values) in units of standard deviation, |u|/σ_r, in a semi-logarithmic scale plot. Here, the solid line shows the symmetrical fit to the exponential CDF [Equation (15)], and the dashed line is the numerical integration of Equation (17) with Equation (12).

Figure 3b shows the number of times a given level of draw-downs (triangles) and -ups ( + signs) of SSB in the absolute values has occurred in the 27 stocks; the abscissa reads the standardized draw-downs (or -ups), |u|/σ_r. There is an apparent symmetry between draw-downs and -ups. The number of draw-downs (or -ups) such that their magnitude is ≥ |u| is well fitted by an exponential law (solid line):

15

where N_uis the total number of draw-downs (or -ups), given by N_u = L_total/2 < τ > , with L_total = 1137, the number of aggregated data points (L_total = ∑_{i = 1}²⁷L_i, where L_i + 1 denotes the length of the population time-series of each stock, i = 1, 2, … , 27). In a random-walk population, it is expected that N_u = 284, which is not too far from the exact number N_Δ = 209 of draw-downs and N₊ = 207 of draw-ups. The mean of persistent moves (standardized data) across all 27 stocks is for draw-downs and for draw-ups. Fitting the time-series data to Equation (15) yields N_u = 213 and u₀ = 1.82 ± 0.07 (adjusted R² = 0.9730), where the draw-downs and -ups are assumed to be symmetrical for the fit. The mean of consecutive drops in the North Atlantic stocks is given by

16

(the mean draw-up is the same magnitude).

Equation (12), numerically integrated, asymptotically yields an exponential distribution. The dashed line in Figure 3b shows the cumulative number of draw-downs (or -ups) for iid Gaussian variations,

17

with the mean draw-downs (or -ups) |〈u〉 |/σ_r = 1.56. One sees an exponential fall-off for a Gaussian process to be compatible with the empirical results.

Discussion

Here, I have studied the nature of time-series of exploited North Atlantic fish stocks. What is remarkable about the result shown in Figure 1is that the distribution is, to a good approximation, symmetrical, so that as many populations are increasing in abundance as are decreasing; such a result would be expected for a random process. Both the analysis of runs of successive log-abundance changes (Figure 2) and the persistence statistics (Figure 3) consistently support the statistical independence of series of successive abundance changes; correlations are absent in the exploited population process. Upward and downward steps of the population process are determined by a “fair coin toss”. The exploited population trajectory in logarithms of abundance is a series of the accumulation of random, uncorrelated numbers over time, leading to random-walk dynamics: the running sum of an iid Gaussian variate r(t) is a random walk of ln S(t). This is natural for a stochastic population process, if one considers the annual change of population abundance to be the result of many independent “shocks”, i.e. replenishing (spawning and recruitment of progeny) and mortality from natural (predation and starvation) and human (fishery) causes; the fluctuations r(t) should be Gaussian.

One might suppose that stocks have been under different and varying exploitation regimes and that such stocks exhibit longer-term trends, e.g. that data (on successive abundance changes) in early years of a population's time-series (perhaps a period of general draw-down in a stock) are correlated in time. The empirical results from the ICES–ACFM time-series data are clear evidence against dependence in the successions of annual rates of population growth.

Moreover, one might suppose that each stock has different distributions of r(t) associated with different exploitation regimes in a population time-series: the volatility, σ_r, of population growth rates would increase with increasing mortality rate, because each year's SSB would, at greater rates of mortality, be based on the contributions of fewer cohorts. Therefore, the distribution of a time-series would exhibit fat-tailed behaviour compared with that of the Gaussian, attributed to mixing of Gaussian distributions with different standard deviations σ_rat various regimes of exploitation (Allen et al., 2001). The results shown in Figure 1are clear evidence against there being crucial differences in distribution, i.e. that population data records {r(t)} in a time-series are identically Gaussian distributed.

If sample means are the same ( < r > = 0) but standard deviations vary among Gaussian distributions, contributing to a mixture, the result from mixing Gaussians is exponentially distributed. Combining population time-series data (unstandardized) across 27 stocks, it is clear that the distribution of population growth rates is not Gaussian, but exponential, i.e. it has a much broader tail than a Gaussian (Figure 4). Therefore, the resulting distribution depends not only on the shape of the tails of each series, but also on the widths of the individual Gaussian distributions. Figure 4plots on a semi-logarithmic scale the CDF of population growth rates with rank ordering: the cumulative number of r ≤ 0 (negative tail) and the complementary cumulative number of r ≥ 0 (positive tail) are fitted to the Laplace (double-sided exponential) CDF:

18

i.e. the number of records of growth rates (absolute values) larger than or equal to |r|. The prefactor N_rreads the total number of records (downward or upward steps); σ_Lis the standard deviation of the growth-rate distribution of the mixture across 27 stocks. The symmetrical fit (dashed line) by the linear regression using the log-transformed data yields N_r = 553 and σ_L = 0.2609 ± 0.0013 (adjusted R² = 0.9744). There is no clear difference between downward and upward steps in the central body of the distribution—the population fluctuations r(t) in each time-series are symmetrical around r = 0. The exact numbers of downward and upward steps are N_{r ≤ 0} = 603 and N_{r ≥ 0} = 534, respectively. The average across all 27 stocks is ⁠, and the empirical standard deviation σ_L = 0.2517. Stock (species and area) differences contribute to differences in the amplitudes of population fluctuations, resulting in an exponential law in population growth-rate distribution when the time-series (unstandardized data) are combined across 27 stocks, and the statistical distribution of time-series values r(t) at various exploitation regimes is the same (i.e. Gaussian) for each stock.

Figure 4.

Open in new tabDownload slide

Mixing Gaussian distributions yields an exponential. The number of times a given level of population growth rate is observed is plotted on a semi-logarithmic scale against unstandardized rates in the combined time-series across the 27 stocks (data from ICES–ACFM Documents of 2005). The dashed line shows the symmetrical fit to Equation (18).

The variance 〈r_Δt²〉 of population abundance grows at a rate proportional to the length of the time-series. The finding of the time-series analysis implies that the population process is non-stationary. It has become accepted (cf. Halley, 2005, and references therein) that population variability increases with census length (i.e. the number of years included in the calculation). Such populations are not regulated towards a stable equilibrium point, nor do they have bounded variabilities. The usual theories on harvesting are under a hypothesis of population regulation (e.g. logistic-type or Ricker, 1954, models) with the replacement line; they imply the existence of an equilibrium point to which the population is attracted, or a “basin of attraction” wherein population abundance is constrained to lie (Beverton and Holt, 1957; Turchin, 1999). There is now some doubt whether populations are persistent in growing (or declining) on a path to equilibrium, which is the usual assumption in fisheries management (Steele and Henderson, 1984). The dynamics of exploited stocks result in a random walk down the population path.

I have here examined the variance of population fluctuations and the persistent duration and magnitude of population moves over time on North Atlantic commercial fish stocks. There is yet another approach to the problem of stationarity and population regulation. It may be valuable to examine frequency variations (the variance per unit frequency interval, in units of y⁻¹) of exploited population time-series, i.e. the spectra of population fluctuations, to investigate whether exploited populations are moving non-stationarily (low-frequency variability has greater amplitude than high-frequency variability, corresponding to unregulated non-equilibrium dynamics) or whether the population process S(t) contains equal-amplitude components of the spectrum at all frequencies (the population spectrum is white noise, corresponding to the extreme case of tight regulation).

In exploited fish populations, time-series analysis of landings has provided empirical evidence against stationarity (Halley and Stergiou, 2005). Besides the random nature of escapement from a fishery (i.e. SSB time-series), it is necessary to examine the statistical properties of the manner in which landed catches from fishing a stock change over time. Because the future harvest from a fish stock is uncertain, the relationship between the escapement of catchable adults from a fishery in one year and landing streams in future is of great importance.

In summary, it is a rational approximation that the annual population growth rate r(t) is an iid Gaussian variate. The exploited population process S(t) is modelled by geometric Brownian motion as the limit of a geometric random walk. The population dynamics are governed by the following stochastic differential equation [equivalent to Equation (2)]:

19

where dB(t) is the Gaussian white-noise term of normal distribution with mean zero and variance dt. B(t) is the standard Brownian motion, i.e. by integrating a white-noise process (the continuous analogue of taking cumulative sums), Brownian motion is obtained.

Acknowledgements

I thank two anonymous reviewers for their constructive comments on two draft versions of the manuscript, and editor Andrew Payne for his encouragement.

References