From cbd01b3da911f64ae858648c2f8cdc9d4a5d2b11 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 07:03:39 +1000 Subject: [PATCH 1/9] Add Maximum Sustainable Yield and Fishery Collapse lecture (#755) New lecture lectures/msy_fishery.md, slotted under Nonlinear Dynamics after solow. Builds the deterministic Schaefer MSY model (logistic growth, 45-degree diagrams, sustainable yield, MSY via grid search and optional calculus), then adds environmental randomness to show why a constant MSY quota sits on a knife-edge and collapses, versus a self-correcting constant-effort policy. Includes the Peruvian anchovy and Atlantic cod collapses, plus exercises on precautionary quotas and large-r logistic chaos. Adds schaefer1954, gordon1954, may1976, clark1990 to quant-econ.bib. Closes #755 Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/_static/quant-econ.bib | 40 ++ lectures/_toc.yml | 1 + lectures/msy_fishery.md | 852 ++++++++++++++++++++++++++++++++ 3 files changed, 893 insertions(+) create mode 100644 lectures/msy_fishery.md diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index 61f867e24..b20cef55b 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -2917,3 +2917,43 @@ @incollection{sargent1989least year={1989}, publisher={Cambridge University Press} } + +@article{schaefer1954, + title={Some aspects of the dynamics of populations important to the management of the commercial marine fisheries}, + author={Schaefer, Milner B.}, + journal={Bulletin of the Inter-American Tropical Tuna Commission}, + volume={1}, + number={2}, + pages={23--56}, + year={1954} +} + +@article{gordon1954, + title={The economic theory of a common-property resource: the fishery}, + author={Gordon, H. Scott}, + journal={Journal of Political Economy}, + volume={62}, + number={2}, + pages={124--142}, + year={1954}, + publisher={University of Chicago Press} +} + +@article{may1976, + title={Simple mathematical models with very complicated dynamics}, + author={May, Robert M.}, + journal={Nature}, + volume={261}, + number={5560}, + pages={459--467}, + year={1976}, + publisher={Nature Publishing Group} +} + +@book{clark1990, + title={Mathematical Bioeconomics: The Optimal Management of Renewable Resources}, + author={Clark, Colin W.}, + edition={2nd}, + year={1990}, + publisher={Wiley} +} diff --git a/lectures/_toc.yml b/lectures/_toc.yml index c9f39454f..eb17a8968 100644 --- a/lectures/_toc.yml +++ b/lectures/_toc.yml @@ -47,6 +47,7 @@ parts: chapters: - file: scalar_dynam - file: solow + - file: msy_fishery - file: cobweb - file: olg - file: commod_price diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md new file mode 100644 index 000000000..5b334bec8 --- /dev/null +++ b/lectures/msy_fishery.md @@ -0,0 +1,852 @@ +--- +jupytext: + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.17.2 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Maximum Sustainable Yield and Fishery Collapse + +## Overview + +In this lecture we study one of the central ideas in the management of renewable +resources: the **maximum sustainable yield** (MSY). + +The MSY is the largest catch a fish stock can support year after year without +running itself down. + +We build the idea up with almost no mathematics. + +The model is just a year-to-year bookkeeping rule, and the best catch can be +found by reading the peak off a graph. + +A short optional section shows how calculus recovers the same answers in closed +form, for readers who want it. + +But the MSY also has a darker side. + +After the Second World War, fisheries managers around the world adopted MSY (and +close relatives) as a target. + +Several great fisheries then collapsed --- among them the Peruvian anchovy in the +1970s and the Atlantic cod off Newfoundland in 1992. + +A major reason is that the simple MSY model ignores **randomness and risk**. + +The ocean is a noisy place, and a policy that looks safe in a deterministic model +can be dangerously fragile once we allow the environment to fluctuate. + +So in the second half of the lecture we add random shocks to the model and watch a +fishery collapse. + +The MSY framework is due to {cite}`schaefer1954` and {cite}`gordon1954`; a +classic textbook treatment is {cite}`clark1990`. + +Let's start with some standard imports. + +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +``` + +We will use a standard parameterization throughout. + +```{code-cell} ipython3 +r = 0.5 # intrinsic growth rate (per year) +K = 1000.0 # carrying capacity (tonnes) +q = 0.01 # catchability coefficient +``` + +## The model + +Let $x_t$ be the **stock biomass** --- the total weight of fish, say in tonnes --- +at the start of year $t$. + +We build the model in two steps: first how the stock grows on its own, then what +fishing does to it. + +### Growth without fishing + +Left alone, the population adds new fish each year according to **logistic +growth**, + +$$ + G(x) = r\,x\left(1 - \frac{x}{K}\right), +$$ (eq:logistic) + +with **intrinsic growth rate** $r > 0$ and **carrying capacity** $K > 0$. + +Growth is small when the stock is small (few spawners) and small again when the +stock is near $K$ (crowding, limited food), and is largest in between. + +With no fishing, next year's stock is just this year's stock plus that growth: + +$$ + x_{t+1} = x_t + G(x_t). +$$ + +Let's encode the growth function and the one-year update rule. + +```{code-cell} ipython3 +def G(x): + "Logistic growth added over one year." + return r * x * (1 - x / K) + +def update(x, E): + "Next year's stock given current stock x and fishing effort E." + return x + G(x) - q * E * x +``` + +(The `update` function already includes fishing through the effort $E$, which we +introduce below; setting $E = 0$ gives the unfished dynamics.) + +A clean way to see where the dynamics lead is a **45-degree diagram**: we plot +next year's stock $x_{t+1}$ against this year's stock $x_t$. + +Wherever the curve crosses the $45^\circ$ line we have $x_{t+1} = x_t$ --- a +**steady state**, a stock that exactly reproduces itself. + +We can then trace the dynamics by "staircasing": from a starting stock go *up* to +the curve (that gives next year's stock), *across* to the $45^\circ$ line (that +becomes this year's stock), and repeat. + +The next function draws such a diagram. + +```{code-cell} ipython3 +def plot_45(ax, E, x0, x_max, steady_state, ss_label, map_label, n_years=30): + "Draw a 45-degree (cobweb) diagram for the yearly stock update at effort E." + grid = np.linspace(0, x_max, 400) + ax.plot(grid, update(grid, E), color='C0', lw=2.5, label=map_label) + ax.plot(grid, grid, color='0.6', lw=1, ls='--', + label=r'$45^\circ$ line $x_{t+1}=x_t$') + # cobweb staircase starting from x0 + x = x0 + cx, cy = [x], [0.0] + for _ in range(n_years): + y = update(x, E) + cx += [x, y] + cy += [y, y] + x = y + ax.plot(cx, cy, color='black', lw=1, alpha=0.9) + ax.plot([steady_state], [steady_state], 'o', color='black', ms=8, zorder=5) + ax.annotate(ss_label, xy=(steady_state, steady_state), + xytext=(0.55 * x_max, 0.28 * x_max), fontsize=11, + arrowprops=dict(arrowstyle='->', color='black', lw=1)) + ax.set_xlabel('stock this year $x_t$ (tonnes)') + ax.set_ylabel('stock next year $x_{t+1}$ (tonnes)') + ax.set_xlim(0, x_max) + ax.set_ylim(0, x_max) + ax.set_aspect('equal') + ax.legend(loc='upper left', frameon=False, fontsize=9) + ax.spines[['top', 'right']].set_visible(False) +``` + +Here is the unfished case. + +```{code-cell} ipython3 +fig, ax = plt.subplots(figsize=(5.5, 5.5)) +plot_45(ax, E=0, x0=80, x_max=1100, + steady_state=K, ss_label=r'$x=K$ (carrying capacity)', + map_label=r'$x_{t+1}=x_t+G(x_t)$') +ax.set_title('Stock dynamics without fishing') +plt.tight_layout() +plt.show() +``` + +Starting from any positive stock, the staircase climbs to the carrying capacity +$x = K$. + +The unfished population fills up the environment and stays there. + +(The empty ocean $x = 0$ is also a steady state, but an unstable one --- any small +stock grows away from it.) + +### Adding fishing + +Now let a fishing fleet remove a catch each year. + +Following {cite}`schaefer1954`, the catch is proportional to fishing **effort** +$E$ (e.g. boat-days) and to the stock available to be caught: + +$$ + h_t = q\,E\,x_t, +$$ (eq:harvest) + +where $q > 0$ is the **catchability coefficient**. + +Subtracting the catch, next year's stock becomes + +$$ + x_{t+1} + \;=\; x_t + \;+\; \underbrace{r\,x_t\!\left(1-\frac{x_t}{K}\right)}_{\text{growth}} + \;-\; \underbrace{qE\,x_t}_{\text{catch}}. +$$ (eq:update) + +That is the whole model --- a rule you could run forward in a spreadsheet. + +On the 45-degree diagram, fishing pulls the update curve **down** by the harvest +term, so it now crosses the $45^\circ$ line at a **lower** steady state. + +A fished population settles below its unfished level. + +```{code-cell} ipython3 +E_demo = 20.0 # an illustrative fixed effort level + +def x_star(E): + "Sustainable (steady-state) stock at effort E." + return K * (1 - q * E / r) + +fig, ax = plt.subplots(figsize=(5.5, 5.5)) +plot_45(ax, E=E_demo, x0=80, x_max=1100, + steady_state=x_star(E_demo), ss_label=r'$x^*(E)$', + map_label=r'$x_{t+1}=x_t+G(x_t)-qEx_t$') +ax.set_title(f'Stock dynamics at a fixed effort $E={E_demo:.0f}$') +plt.tight_layout() +plt.show() +``` + +The natural question is now: which effort level $E$ gives the largest catch we can +take year after year? + +## When is a catch sustainable? + +Suppose the fleet applies a constant effort $E$ every year. + +The catch is **sustainable** if the stock holds steady from one year to the next +--- it neither grows without limit nor dwindles away. + +In symbols, the stock repeats itself: + +$$ + x_{t+1} = x_t. +$$ + +Looking at {eq}`eq:update`, this happens exactly when growth replaces the catch, +$G(x) = qEx$. + +Writing it out and cancelling the empty-ocean case $x = 0$: + +$$ + r\left(1 - \frac{x}{K}\right) = qE + \quad\Longrightarrow\quad + x^*(E) = K\left(1 - \frac{qE}{r}\right). +$$ (eq:xstar) + +So **each effort level $E$ pins down one steady-state stock** $x^*(E)$ --- +provided $E < r/q$. + +Any more effort than that drives the stock to zero. + +The catch this steady state delivers, year after year, is + +$$ + Y(E) = q\,E\,x^*(E). +$$ (eq:yield) + +```{code-cell} ipython3 +def sustainable_yield(E): + "Catch that the steady state delivers each year at effort E." + return q * E * x_star(E) +``` + +Notice that finding this needed only algebra --- no calculus. + +## The maximum sustainable yield + +Different effort levels give different sustainable catches $Y(E)$. + +The **maximum sustainable yield** is simply the largest of them: + +$$ + \text{MSY} \;=\; \max_{0 \le E < r/q} \; Y(E). +$$ + +We can find it without any calculus: evaluate $Y(E)$ on a fine grid of effort +levels and pick the biggest. + +```{code-cell} ipython3 +Es = np.linspace(0, r / q, 100_001) # efforts from 0 up to the wipe-out level r/q +Ys = sustainable_yield(Es) +best = np.argmax(Ys) + +E_msy = Es[best] +MSY = Ys[best] +x_msy = x_star(E_msy) + +print(f"Best effort E_MSY = {E_msy:.2f}") +print(f"Sustainable stock x* = {x_msy:.1f} tonnes (= K/2)") +print(f"Max sustainable yield MSY = {MSY:.1f} tonnes/year (= rK/4)") +``` + +The grid search lands on a sustainable stock of $x^* = 500 = K/2$ and a maximum +catch of $\text{MSY} = 125 = rK/4$. + +Those round numbers are no accident --- the optional calculus section below shows +where they come from. + +## Visualizing the MSY + +### The sustainable-yield curve + +As effort climbs from $0$ toward $r/q$, the sustainable stock $x^*(E)$ slides from +$K$ down to $0$, and the catch it supports traces out the growth curve $G(x)$. + +So the hump below is the full menu of sustainable (stock, catch) pairs --- and its +peak is the MSY. + +```{code-cell} ipython3 +x = np.linspace(0, K, 400) + +fig, ax = plt.subplots(figsize=(8, 5)) +ax.plot(x, G(x), lw=2.5, color='C0', + label=r'sustainable yield $Y=G(x)=rx(1-x/K)$') + +ax.plot([x_msy], [MSY], 'o', color='black', ms=9, zorder=5) +ax.vlines(x_msy, 0, MSY, ls='--', color='black', lw=1) +ax.hlines(MSY, 0, x_msy, ls='--', color='black', lw=1) +ax.annotate(f'MSY = rK/4 = {MSY:.0f}', + xy=(x_msy, MSY), xytext=(x_msy + 60, MSY - 18), + fontsize=11, color='black') +ax.annotate(r'$x_{MSY}=K/2$', xy=(x_msy, 0), xytext=(x_msy + 10, 6), + fontsize=11, color='black') + +ax.set_xlabel('stock biomass $x$ (tonnes)') +ax.set_ylabel('growth / sustainable yield (tonnes/year)') +ax.set_title('Maximum sustainable yield: the peak of the growth curve') +ax.set_xlim(0, K) +ax.set_ylim(0, MSY * 1.25) +ax.legend(loc='upper right', frameon=False) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +### Growth versus harvest + +Now plot the growth curve $G(x)$ together with the harvest lines $h = qEx$ for +several effort levels. + +Each crossing of a harvest line with the growth curve is a sustainable steady +state $x^*(E)$, and the catch there is the height of the line. + +The MSY effort line meets the growth curve exactly at its peak, $x = K/2$. + +```{code-cell} ipython3 +fig, ax = plt.subplots(figsize=(8, 5)) +ax.plot(x, G(x), lw=2.5, color='C0', label='growth $G(x)$') + +efforts = [0.5 * E_msy, E_msy, 1.5 * E_msy] +labels = [r'$E < E_{MSY}$ (underfishing)', + r'$E = E_{MSY}$', + r'$E > E_{MSY}$ (overfishing)'] +colors = ['C2', 'C3', 'C1'] + +for E, lab, c in zip(efforts, labels, colors): + ax.plot(x, q * E * x, lw=1.8, color=c, label=lab) + ax.plot([x_star(E)], [q * E * x_star(E)], 'o', color=c, ms=7, zorder=5) + +ax.set_xlabel('stock biomass $x$ (tonnes)') +ax.set_ylabel('rate (tonnes/year)') +ax.set_title('Sustainable steady states: growth $G(x)$ vs. harvest $qEx$') +ax.set_xlim(0, K) +ax.set_ylim(0, MSY * 1.4) +ax.legend(loc='upper left', frameon=False, fontsize=10) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +Notice that the under- and over-fishing lines cross the growth curve at the +**same height** (same catch) but at different stock levels. + +Only $E_{MSY}$ reaches the very top. + +Increasing effort past $E_{MSY}$ buys a *smaller* sustained catch from a *more +depleted* stock --- biologically and economically the worst of both worlds. + +### The yield-effort curve + +Plotting the sustainable catch $Y(E)$ directly against effort gives the classic +dome-shaped Schaefer curve that fisheries managers use. + +```{code-cell} ipython3 +E_grid = np.linspace(0, r / q, 400) +Y_grid = sustainable_yield(E_grid) + +fig, ax = plt.subplots(figsize=(8, 5)) +ax.plot(E_grid, Y_grid, lw=2.5, color='C0', label=r'$Y(E)=qEK\,(1-qE/r)$') +ax.plot([E_msy], [MSY], 'o', color='black', ms=9, zorder=5) +ax.vlines(E_msy, 0, MSY, ls='--', color='black', lw=1) +ax.hlines(MSY, 0, E_msy, ls='--', color='black', lw=1) +ax.annotate(f'MSY = {MSY:.0f}', xy=(E_msy, MSY), xytext=(E_msy + 1, MSY - 18), + fontsize=11, color='black') +ax.annotate(r'$E_{MSY}=r/2q$', xy=(E_msy, 0), xytext=(E_msy + 1, 6), + fontsize=11, color='black') + +ax.fill_between(E_grid, 0, Y_grid, where=(E_grid > E_msy), color='C1', alpha=0.15) +ax.text(E_msy * 1.45, MSY * 0.35, 'overfishing\n(yield falls)', + color='C1', fontsize=10, ha='center') + +ax.set_xlabel('fishing effort $E$') +ax.set_ylabel('sustainable yield $Y$ (tonnes/year)') +ax.set_title('Yield-effort curve (Schaefer model)') +ax.set_xlim(0, r / q) +ax.set_ylim(0, MSY * 1.25) +ax.legend(loc='upper right', frameon=False) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +### Does the stock settle there? + +Finally, let's run the yearly recursion {eq}`eq:update` forward from several +starting stocks, under the MSY effort. + +```{code-cell} ipython3 +def simulate(x0, E, years=40): + "Run the deterministic yearly stock recursion forward from x0." + x = np.empty(years + 1) + x[0] = x0 + for t in range(years): + x[t + 1] = update(x[t], E) + return np.arange(years + 1), x + +fig, ax = plt.subplots(figsize=(8, 5)) +for x0 in [100, 300, 700, 950]: + t, xt = simulate(x0, E_msy) + ax.plot(t, xt, 'o-', ms=3, lw=1.2, label=f'$x_0={x0}$') + +ax.axhline(x_msy, ls='--', color='black', lw=1.5) +ax.annotate(r'$x^*=K/2$', xy=(0, x_msy), xytext=(1, x_msy + 25), color='black') +ax.set_xlabel('year $t$') +ax.set_ylabel('stock biomass $x_t$ (tonnes)') +ax.set_title(r'Year-by-year stock path under MSY effort $E=E_{MSY}$') +ax.set_xlim(0, 40) +ax.legend(frameon=False) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +Under the MSY effort, every path climbs or falls toward $x^* = K/2$ and stays. + +So the catch settles down to the MSY. + +```{note} +The smooth, steady convergence above relies on $r$ being modest. + +For large $r$ the yearly model can overshoot and oscillate --- even behave +chaotically {cite}`may1976` --- but that is an artifact of the one-year time step, +not the biology. + +The continuous-time version in the next section converges smoothly for *any* +$r > 0$, which is one reason it is the textbook standard. +``` + +## The same result with calculus + +```{note} +This section is optional. + +Everything above used only algebra and a graph; readers who have met calculus can +use it to pin down the round numbers exactly. +``` + +Shrink the time step to zero and the yearly recursion {eq}`eq:update` becomes the +differential equation + +$$ + \dot{x} = G(x) - qEx. +$$ + +The sustainability condition is unchanged --- a steady state still needs +$G(x) = qEx$ --- so $x^*(E) = K(1 - qE/r)$ and $Y(E) = qEK(1 - qE/r)$ exactly as +before. + +Now $Y$ is a smooth, concave parabola in $E$, and its peak is where the slope +vanishes: + +$$ + Y'(E) = qK\left(1 - \frac{2qE}{r}\right) = 0 + \quad\Longrightarrow\quad + E_{MSY} = \frac{r}{2q}, +$$ + +which gives $x^* = K/2$ and $\text{MSY} = rK/4$ --- exactly the values the grid +search found. + +Equivalently, since the sustainable catch equals growth, $Y = G(x^*)$, the MSY is +just the top of the growth curve, where $G'(x) = 0$ at $x = K/2$. + +In words: **fish the stock down to where it grows fastest, and no further.** + +```{note} +Readers who have seen the Solow-Swan growth model may recognize the structure. + +A single control --- here fishing effort, there the savings rate --- picks one +point on a one-parameter family of steady states. + +The sustainable flow --- here the yield, there consumption --- is hump-shaped in +that control, and the optimum sits at an interior peak. + +Over-fishing is then the exact analogue of over-saving (dynamic inefficiency): +more of the control, less of the flow. +``` + +## Randomness, risk and collapse + +So far the ocean has been a perfectly predictable place. + +Each year the stock grows by exactly $G(x_t)$. + +Real fish populations are not like this. + +Recruitment depends on water temperature, currents, food supply, predators and +disease --- all of which vary from year to year, sometimes wildly. + +Let's add this randomness to the model. + +A simple way is to multiply each year's growth by a random **environmental shock** +$\xi_t$ with mean one: + +$$ + x_{t+1} = x_t + \xi_t \, G(x_t) - h_t, +$$ (eq:stochastic) + +where $h_t$ is the catch taken in year $t$. + +We take $\xi_t$ to be lognormal with mean one, controlled by a volatility +parameter $\sigma$. + +```{code-cell} ipython3 +def shocks(σ, n, rng): + "Return n lognormal shocks with mean one and volatility σ." + return np.exp(σ * rng.standard_normal(n) - σ**2 / 2) +``` + +Now we face a management choice that did not matter in the deterministic world but +matters enormously here. + +How should we set the catch $h_t$? + +We compare two policies, both designed to take the MSY *on average*: + +* **constant effort**: set $h_t = q E_{MSY} x_t$, so the catch scales with the + current stock, and +* **constant quota**: set $h_t = \text{MSY}$ every year, regardless of the stock. + +The constant quota is the more literal reading of "take the maximum sustainable +yield each year". + +It is also closer to how catch limits have often been set in practice. + +Let's simulate both, starting every fishery from the MSY stock $x^* = K/2$. + +```{code-cell} ipython3 +def simulate_stochastic(policy, σ, years, rng, x0=x_msy): + "Simulate the stochastic fishery under a given harvest policy." + x = np.empty(years + 1) + x[0] = x0 + ξ = shocks(σ, years, rng) + for t in range(years): + growth = ξ[t] * G(x[t]) + if policy == 'effort': + harvest = q * E_msy * x[t] + else: # constant quota + harvest = MSY + x[t + 1] = max(0.0, x[t] + growth - harvest) + return x +``` + +Here are a few sample paths under each policy, with a moderate amount of +environmental noise. + +```{code-cell} ipython3 +σ = 0.15 +years = 100 +rng = np.random.default_rng(seed=1) + +fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True) + +for ax, policy, title in zip( + axes, ['effort', 'quota'], + ['constant effort $h_t = qE_{MSY}x_t$', + 'constant quota $h_t = \\mathrm{MSY}$']): + for k in range(8): + path = simulate_stochastic(policy, σ, years, rng) + ax.plot(path, lw=1, alpha=0.8) + ax.axhline(x_msy, ls='--', color='black', lw=1) + ax.set_title(title) + ax.set_xlabel('year $t$') + ax.set_xlim(0, years) + ax.set_ylim(0, K) + ax.spines[['top', 'right']].set_visible(False) + +axes[0].set_ylabel('stock biomass $x_t$ (tonnes)') +axes[0].annotate(r'$x^*=K/2$', xy=(years, x_msy), xytext=(-20, 8), + textcoords='offset points', color='black') +fig.suptitle(f'Two MSY policies under environmental noise ($\\sigma={σ}$)') +plt.tight_layout() +plt.show() +``` + +The difference is stark. + +Under **constant effort**, the paths jiggle around $x^* = K/2$ but always recover: +a bad year cuts the catch (because the catch is proportional to the stock), giving +the population room to bounce back. + +The policy is **self-correcting**. + +Under **constant quota**, several paths slide down and hit zero --- the fishery +**collapses**. + +The quota keeps demanding the full MSY even after a run of bad years has thinned +the stock, and eventually the catch exceeds what the depleted population can +replace. + +Once that happens there is no way back. + +### Why the quota is a knife-edge + +The collapse is not bad luck --- it is built into the constant-quota policy. + +Set the catch to the constant $h = \text{MSY} = rK/4$ in the deterministic model +and look for steady states, $G(x) = rK/4$: + +$$ + r\,x\left(1 - \frac{x}{K}\right) = \frac{rK}{4} + \quad\Longrightarrow\quad + \left(x - \frac{K}{2}\right)^2 = 0. +$$ + +The two steady states have merged into a single **tangent** point at $x = K/2$. + +This point is *semi-stable*: stable from above but unstable from below. + +If the stock ever drifts below $K/2$, then $G(x) < rK/4$, so the constant quota +takes more than the stock can grow, and the stock falls further --- a one-way +ratchet down to zero. + +Fishing exactly at the MSY with a fixed quota leaves **no margin for error**. + +The deterministic model hides this fragility because nothing ever pushes the stock +off its perch. + +Add even a little randomness and the knife-edge reveals itself. + +### How often does the fishery collapse? + +Let's quantify the risk by simulating many fisheries under each policy and +counting how often the stock collapses within 100 years. + +```{code-cell} ipython3 +def collapse_fraction(policy, σ, years=100, n_paths=2000, seed=0): + "Fraction of simulated fisheries that collapse within the horizon." + rng = np.random.default_rng(seed) + collapses = 0 + for _ in range(n_paths): + path = simulate_stochastic(policy, σ, years, rng) + if path[-1] < 1.0: + collapses += 1 + return collapses / n_paths + +for policy in ['effort', 'quota']: + frac = collapse_fraction(policy, σ=0.15) + print(f"{policy:>8s} policy: collapse within 100 years = {frac:.1%}") +``` + +With this level of noise the constant-effort fishery essentially never collapses, +while the constant-quota fishery collapses much of the time. + +The two policies deliver the same average catch in calm conditions, yet one is +robust and the other is fragile. + +This is the heart of the matter: **the MSY is a deterministic, steady-state +concept, and steering straight at it ignores the risk that randomness creates.** + +## Historical collapses + +These are not just theoretical worries. + +The **Peruvian anchovy** fishery was the largest in the world in the early 1970s. + +Management was guided by sustainable-yield calculations, but those calculations +left out the **El Niño** warming events that periodically disrupt the cold, +nutrient-rich currents the anchovy depend on. + +When a strong El Niño struck in 1972--73, the stock --- already fished hard --- +collapsed, taking tens of thousands of jobs with it. + +The **Atlantic cod** off Newfoundland tells a similar story. + +Catch limits set with too little margin, combined with overoptimistic stock +estimates, allowed fishing to continue as the population fell. + +In 1992 the stock had dropped by more than 90% and Canada closed the fishery +entirely; it has still not fully recovered. + +In both cases the deterministic logic of "take the maximum sustainable yield" +proved dangerously fragile once the real, noisy ocean was taken into account. + +The lesson modern fisheries science has drawn is to manage **precautionarily** --- +to aim below the MSY, leaving a buffer against the bad years that the simple model +pretends do not exist. + +## Exercises + +```{exercise} +:label: msy_ex1 + +The constant-quota policy above demanded the *full* MSY every year, and we saw +that this is a knife-edge. + +A natural fix is to be **precautionary**: take a fixed quota that is only a +fraction $\alpha \in (0, 1]$ of the MSY. + +For each $\alpha$ in a grid from $0.5$ to $1.0$, estimate the probability that the +fishery collapses within 100 years under the constant quota $h = \alpha \cdot +\text{MSY}$, using $\sigma = 0.15$. + +Plot the collapse probability against $\alpha$. + +How much does pulling the quota back below the MSY reduce the risk? +``` + +```{solution-start} msy_ex1 +:class: dropdown +``` + +We adapt `simulate_stochastic` to take an arbitrary fixed quota. + +```{code-cell} ipython3 +def collapse_fraction_quota(quota, σ=0.15, years=100, n_paths=2000, seed=0): + "Collapse probability under a fixed quota." + rng = np.random.default_rng(seed) + collapses = 0 + for _ in range(n_paths): + x = x_msy + ξ = shocks(σ, years, rng) + for t in range(years): + x = max(0.0, x + ξ[t] * G(x) - quota) + if x <= 0.0: + break + if x < 1.0: + collapses += 1 + return collapses / n_paths + +alphas = np.linspace(0.5, 1.0, 11) +probs = [collapse_fraction_quota(α * MSY) for α in alphas] + +fig, ax = plt.subplots() +ax.plot(alphas, probs, 'o-', lw=2) +ax.set_xlabel(r'quota as a fraction $\alpha$ of MSY') +ax.set_ylabel('collapse probability within 100 years') +ax.set_title('Precaution pays: lower quotas are far safer') +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +The collapse probability falls steeply as the quota is pulled back below the MSY. + +A modest precautionary margin --- taking, say, 80% of the MSY --- buys a large +reduction in the risk of collapse, at the cost of only a small reduction in the +average catch. + +This is exactly the trade-off that precautionary fisheries management is built +around. + +```{solution-end} +``` + +```{exercise} +:label: msy_ex2 + +The note after the convergence plot warned that for large $r$ the *yearly* +recursion can misbehave, even though the underlying biology is benign. + +Investigate this for the unfished model $x_{t+1} = x_t + r x_t(1 - x_t/K)$. + +Rescale by writing $u_t = x_t / K$, so that + +$$ + u_{t+1} = u_t + r\,u_t(1 - u_t). +$$ + +Simulate this map from $u_0 = 0.3$ for $r = 1.5$, $r = 2.2$ and $r = 2.7$, and plot +the three paths. + +Describe what happens as $r$ increases. +``` + +```{solution-start} msy_ex2 +:class: dropdown +``` + +```{code-cell} ipython3 +def logistic_path(r_val, u0=0.3, years=40): + u = np.empty(years + 1) + u[0] = u0 + for t in range(years): + u[t + 1] = u[t] + r_val * u[t] * (1 - u[t]) + return u + +fig, ax = plt.subplots(figsize=(8, 5)) +for r_val in [1.5, 2.2, 2.7]: + ax.plot(logistic_path(r_val), 'o-', ms=3, lw=1, label=f'$r={r_val}$') +ax.axhline(1.0, ls='--', color='black', lw=1) +ax.set_xlabel('year $t$') +ax.set_ylabel(r'scaled stock $u_t = x_t / K$') +ax.set_title('The yearly logistic map as $r$ grows') +ax.legend(frameon=False) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +For $r = 1.5$ the path converges smoothly to the carrying capacity $u = 1$. + +For $r = 2.2$ it overshoots and settles into a steady oscillation (a two-year +cycle). + +For $r = 2.7$ the path never settles --- it bounces around erratically, the onset +of the chaotic behavior analyzed by {cite}`may1976`. + +The biology is unchanged; it is the large discrete time step that manufactures +these wild dynamics, which is why the smooth continuous-time version is the +textbook standard. + +```{solution-end} +``` + +## Summary + +The key formulas of the deterministic Schaefer model are collected below. + +| Quantity | Formula | Value | +|---|---|---| +| Stock at MSY | $x_{MSY} = K/2$ | 500 t | +| Maximum sustainable yield | $\text{MSY} = rK/4$ | 125 t/yr | +| Effort at MSY | $E_{MSY} = r/2q$ | 25 | + +The MSY is a deterministic, single-species, steady-state concept. + +It ignores economic costs and prices, age and size structure, environmental +randomness, and interactions with other species. + +As we saw, the omission of randomness is not a minor technicality. + +A constant quota at the MSY sits on a knife-edge, and in a noisy ocean it leads +the fishery to collapse. + +Real fisheries management therefore treats the MSY as a reference point --- often +a cautionary upper bound --- rather than a precise target. From 7c0af033606fc9f7d81cd1e7111b76cc269bd4fe Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 09:56:34 +1000 Subject: [PATCH 2/9] Revise MSY lecture: rename effort to e, yield to y*(e), add early graphic - Rename effort symbol E -> e and sustainable yield Y(E) -> y*(e) throughout prose, math, and code, matching the x*(e) notation. - Introduce y*(e) graphically with a growth-vs-harvest figure right where it is defined, marking x*(e) and y*(e) at effort e=20. - Shrink the two 45-degree diagrams by 10% (figsize 5.5 -> 4.95). Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 224 ++++++++++++++++++++++++---------------- 1 file changed, 134 insertions(+), 90 deletions(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index 5b334bec8..f2b3518f2 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -98,13 +98,13 @@ def G(x): "Logistic growth added over one year." return r * x * (1 - x / K) -def update(x, E): - "Next year's stock given current stock x and fishing effort E." - return x + G(x) - q * E * x +def update(x, e): + "Next year's stock given current stock x and fishing effort e." + return x + G(x) - q * e * x ``` -(The `update` function already includes fishing through the effort $E$, which we -introduce below; setting $E = 0$ gives the unfished dynamics.) +(The `update` function already includes fishing through the effort $e$, which we +introduce below; setting $e = 0$ gives the unfished dynamics.) A clean way to see where the dynamics lead is a **45-degree diagram**: we plot next year's stock $x_{t+1}$ against this year's stock $x_t$. @@ -119,17 +119,17 @@ becomes this year's stock), and repeat. The next function draws such a diagram. ```{code-cell} ipython3 -def plot_45(ax, E, x0, x_max, steady_state, ss_label, map_label, n_years=30): - "Draw a 45-degree (cobweb) diagram for the yearly stock update at effort E." +def plot_45(ax, e, x0, x_max, steady_state, ss_label, map_label, n_years=30): + "Draw a 45-degree (cobweb) diagram for the yearly stock update at effort e." grid = np.linspace(0, x_max, 400) - ax.plot(grid, update(grid, E), color='C0', lw=2.5, label=map_label) + ax.plot(grid, update(grid, e), color='C0', lw=2.5, label=map_label) ax.plot(grid, grid, color='0.6', lw=1, ls='--', label=r'$45^\circ$ line $x_{t+1}=x_t$') # cobweb staircase starting from x0 x = x0 cx, cy = [x], [0.0] for _ in range(n_years): - y = update(x, E) + y = update(x, e) cx += [x, y] cy += [y, y] x = y @@ -150,8 +150,8 @@ def plot_45(ax, E, x0, x_max, steady_state, ss_label, map_label, n_years=30): Here is the unfished case. ```{code-cell} ipython3 -fig, ax = plt.subplots(figsize=(5.5, 5.5)) -plot_45(ax, E=0, x0=80, x_max=1100, +fig, ax = plt.subplots(figsize=(4.95, 4.95)) +plot_45(ax, e=0, x0=80, x_max=1100, steady_state=K, ss_label=r'$x=K$ (carrying capacity)', map_label=r'$x_{t+1}=x_t+G(x_t)$') ax.set_title('Stock dynamics without fishing') @@ -172,10 +172,10 @@ stock grows away from it.) Now let a fishing fleet remove a catch each year. Following {cite}`schaefer1954`, the catch is proportional to fishing **effort** -$E$ (e.g. boat-days) and to the stock available to be caught: +$e$ (e.g. boat-days) and to the stock available to be caught: $$ - h_t = q\,E\,x_t, + h_t = q\,e\,x_t, $$ (eq:harvest) where $q > 0$ is the **catchability coefficient**. @@ -186,7 +186,7 @@ $$ x_{t+1} \;=\; x_t \;+\; \underbrace{r\,x_t\!\left(1-\frac{x_t}{K}\right)}_{\text{growth}} - \;-\; \underbrace{qE\,x_t}_{\text{catch}}. + \;-\; \underbrace{qe\,x_t}_{\text{catch}}. $$ (eq:update) That is the whole model --- a rule you could run forward in a spreadsheet. @@ -197,27 +197,27 @@ term, so it now crosses the $45^\circ$ line at a **lower** steady state. A fished population settles below its unfished level. ```{code-cell} ipython3 -E_demo = 20.0 # an illustrative fixed effort level +e_demo = 20.0 # an illustrative fixed effort level -def x_star(E): - "Sustainable (steady-state) stock at effort E." - return K * (1 - q * E / r) +def x_star(e): + "Sustainable (steady-state) stock at effort e." + return K * (1 - q * e / r) -fig, ax = plt.subplots(figsize=(5.5, 5.5)) -plot_45(ax, E=E_demo, x0=80, x_max=1100, - steady_state=x_star(E_demo), ss_label=r'$x^*(E)$', - map_label=r'$x_{t+1}=x_t+G(x_t)-qEx_t$') -ax.set_title(f'Stock dynamics at a fixed effort $E={E_demo:.0f}$') +fig, ax = plt.subplots(figsize=(4.95, 4.95)) +plot_45(ax, e=e_demo, x0=80, x_max=1100, + steady_state=x_star(e_demo), ss_label=r'$x^*(e)$', + map_label=r'$x_{t+1}=x_t+G(x_t)-qex_t$') +ax.set_title(f'Stock dynamics at a fixed effort $e={e_demo:.0f}$') plt.tight_layout() plt.show() ``` -The natural question is now: which effort level $E$ gives the largest catch we can +The natural question is now: which effort level $e$ gives the largest catch we can take year after year? ## When is a catch sustainable? -Suppose the fleet applies a constant effort $E$ every year. +Suppose the fleet applies a constant effort $e$ every year. The catch is **sustainable** if the stock holds steady from one year to the next --- it neither grows without limit nor dwindles away. @@ -229,58 +229,100 @@ $$ $$ Looking at {eq}`eq:update`, this happens exactly when growth replaces the catch, -$G(x) = qEx$. +$G(x) = qex$. Writing it out and cancelling the empty-ocean case $x = 0$: $$ - r\left(1 - \frac{x}{K}\right) = qE + r\left(1 - \frac{x}{K}\right) = qe \quad\Longrightarrow\quad - x^*(E) = K\left(1 - \frac{qE}{r}\right). + x^*(e) = K\left(1 - \frac{qe}{r}\right). $$ (eq:xstar) -So **each effort level $E$ pins down one steady-state stock** $x^*(E)$ --- -provided $E < r/q$. +So **each effort level $e$ pins down one steady-state stock** $x^*(e)$ --- +provided $e < r/q$. Any more effort than that drives the stock to zero. -The catch this steady state delivers, year after year, is +The catch this steady state delivers, year after year, is the **sustainable +yield** $$ - Y(E) = q\,E\,x^*(E). + y^*(e) = q\,e\,x^*(e). $$ (eq:yield) ```{code-cell} ipython3 -def sustainable_yield(E): - "Catch that the steady state delivers each year at effort E." - return q * E * x_star(E) +def sustainable_yield(e): + "Sustainable catch delivered each year at effort e." + return q * e * x_star(e) ``` -Notice that finding this needed only algebra --- no calculus. +We can read both of these quantities --- the steady-state stock $x^*(e)$ and the +sustainable catch $y^*(e)$ --- straight off a picture. + +Plot the growth curve $G(x)$ together with the harvest line $q e x$ for a single +effort level. + +The two cross where growth exactly replaces the catch: that crossing sits at the +steady-state stock $x^*(e)$. + +The height of the line there is the sustainable catch $y^*(e)$. + +```{code-cell} ipython3 +x = np.linspace(0, K, 400) + +fig, ax = plt.subplots(figsize=(8, 5)) +ax.plot(x, G(x), lw=2.5, color='C0', label=r'growth $G(x)$') +ax.plot(x, q * e_demo * x, lw=2, color='C3', label=r'harvest $q e x$') + +xs = x_star(e_demo) +ys = sustainable_yield(e_demo) +ax.plot([xs], [ys], 'o', color='black', ms=8, zorder=5) +ax.vlines(xs, 0, ys, ls='--', color='black', lw=1) +ax.hlines(ys, 0, xs, ls='--', color='black', lw=1) +ax.annotate(r'$x^*(e)$', xy=(xs, 0), xytext=(xs + 12, 8), fontsize=12) +ax.annotate(r'$y^*(e)$', xy=(0, ys), xytext=(10, ys + 5), fontsize=12) + +ax.set_xlabel('stock biomass $x$ (tonnes)') +ax.set_ylabel('rate (tonnes/year)') +ax.set_title(f'Sustainable stock and catch at effort $e={e_demo:.0f}$') +ax.set_xlim(0, K) +ax.set_ylim(0, G(K / 2) * 1.4) +ax.legend(loc='upper right', frameon=False) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` + +As we change $e$, the harvest line pivots about the origin, the crossing slides +along the growth curve, and the sustainable catch $y^*(e)$ is always the height of +that crossing. + +Notice that finding all this needed only algebra --- no calculus. ## The maximum sustainable yield -Different effort levels give different sustainable catches $Y(E)$. +Different effort levels give different sustainable catches $y^*(e)$. The **maximum sustainable yield** is simply the largest of them: $$ - \text{MSY} \;=\; \max_{0 \le E < r/q} \; Y(E). + \text{MSY} \;=\; \max_{0 \le e < r/q} \; y^*(e). $$ -We can find it without any calculus: evaluate $Y(E)$ on a fine grid of effort +We can find it without any calculus: evaluate $y^*(e)$ on a fine grid of effort levels and pick the biggest. ```{code-cell} ipython3 -Es = np.linspace(0, r / q, 100_001) # efforts from 0 up to the wipe-out level r/q -Ys = sustainable_yield(Es) -best = np.argmax(Ys) +e_vals = np.linspace(0, r / q, 100_001) # efforts from 0 up to the wipe-out level r/q +y_vals = sustainable_yield(e_vals) +best = np.argmax(y_vals) -E_msy = Es[best] -MSY = Ys[best] -x_msy = x_star(E_msy) +e_msy = e_vals[best] +MSY = y_vals[best] +x_msy = x_star(e_msy) -print(f"Best effort E_MSY = {E_msy:.2f}") +print(f"Best effort e_MSY = {e_msy:.2f}") print(f"Sustainable stock x* = {x_msy:.1f} tonnes (= K/2)") print(f"Max sustainable yield MSY = {MSY:.1f} tonnes/year (= rK/4)") ``` @@ -295,7 +337,7 @@ where they come from. ### The sustainable-yield curve -As effort climbs from $0$ toward $r/q$, the sustainable stock $x^*(E)$ slides from +As effort climbs from $0$ toward $r/q$, the sustainable stock $x^*(e)$ slides from $K$ down to $0$, and the catch it supports traces out the growth curve $G(x)$. So the hump below is the full menu of sustainable (stock, catch) pairs --- and its @@ -306,7 +348,7 @@ x = np.linspace(0, K, 400) fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(x, G(x), lw=2.5, color='C0', - label=r'sustainable yield $Y=G(x)=rx(1-x/K)$') + label=r'sustainable yield $y^*=G(x)=rx(1-x/K)$') ax.plot([x_msy], [MSY], 'o', color='black', ms=9, zorder=5) ax.vlines(x_msy, 0, MSY, ls='--', color='black', lw=1) @@ -330,11 +372,13 @@ plt.show() ### Growth versus harvest -Now plot the growth curve $G(x)$ together with the harvest lines $h = qEx$ for -several effort levels. +We just saw how a single harvest line picks out one sustainable point. + +To compare effort levels, plot the growth curve $G(x)$ together with the harvest +lines $q e x$ for several values of $e$. -Each crossing of a harvest line with the growth curve is a sustainable steady -state $x^*(E)$, and the catch there is the height of the line. +Each crossing is a sustainable steady state $x^*(e)$, and the catch there is the +height of the line. The MSY effort line meets the growth curve exactly at its peak, $x = K/2$. @@ -342,19 +386,19 @@ The MSY effort line meets the growth curve exactly at its peak, $x = K/2$. fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(x, G(x), lw=2.5, color='C0', label='growth $G(x)$') -efforts = [0.5 * E_msy, E_msy, 1.5 * E_msy] -labels = [r'$E < E_{MSY}$ (underfishing)', - r'$E = E_{MSY}$', - r'$E > E_{MSY}$ (overfishing)'] +efforts = [0.5 * e_msy, e_msy, 1.5 * e_msy] +labels = [r'$e < e_{MSY}$ (underfishing)', + r'$e = e_{MSY}$', + r'$e > e_{MSY}$ (overfishing)'] colors = ['C2', 'C3', 'C1'] -for E, lab, c in zip(efforts, labels, colors): - ax.plot(x, q * E * x, lw=1.8, color=c, label=lab) - ax.plot([x_star(E)], [q * E * x_star(E)], 'o', color=c, ms=7, zorder=5) +for e, lab, c in zip(efforts, labels, colors): + ax.plot(x, q * e * x, lw=1.8, color=c, label=lab) + ax.plot([x_star(e)], [q * e * x_star(e)], 'o', color=c, ms=7, zorder=5) ax.set_xlabel('stock biomass $x$ (tonnes)') ax.set_ylabel('rate (tonnes/year)') -ax.set_title('Sustainable steady states: growth $G(x)$ vs. harvest $qEx$') +ax.set_title('Sustainable steady states: growth $G(x)$ vs. harvest $qex$') ax.set_xlim(0, K) ax.set_ylim(0, MSY * 1.4) ax.legend(loc='upper left', frameon=False, fontsize=10) @@ -366,36 +410,36 @@ plt.show() Notice that the under- and over-fishing lines cross the growth curve at the **same height** (same catch) but at different stock levels. -Only $E_{MSY}$ reaches the very top. +Only $e_{MSY}$ reaches the very top. -Increasing effort past $E_{MSY}$ buys a *smaller* sustained catch from a *more +Increasing effort past $e_{MSY}$ buys a *smaller* sustained catch from a *more depleted* stock --- biologically and economically the worst of both worlds. ### The yield-effort curve -Plotting the sustainable catch $Y(E)$ directly against effort gives the classic +Plotting the sustainable catch $y^*(e)$ directly against effort gives the classic dome-shaped Schaefer curve that fisheries managers use. ```{code-cell} ipython3 -E_grid = np.linspace(0, r / q, 400) -Y_grid = sustainable_yield(E_grid) +e_grid = np.linspace(0, r / q, 400) +y_grid = sustainable_yield(e_grid) fig, ax = plt.subplots(figsize=(8, 5)) -ax.plot(E_grid, Y_grid, lw=2.5, color='C0', label=r'$Y(E)=qEK\,(1-qE/r)$') -ax.plot([E_msy], [MSY], 'o', color='black', ms=9, zorder=5) -ax.vlines(E_msy, 0, MSY, ls='--', color='black', lw=1) -ax.hlines(MSY, 0, E_msy, ls='--', color='black', lw=1) -ax.annotate(f'MSY = {MSY:.0f}', xy=(E_msy, MSY), xytext=(E_msy + 1, MSY - 18), +ax.plot(e_grid, y_grid, lw=2.5, color='C0', label=r'$y^*(e)=qeK\,(1-qe/r)$') +ax.plot([e_msy], [MSY], 'o', color='black', ms=9, zorder=5) +ax.vlines(e_msy, 0, MSY, ls='--', color='black', lw=1) +ax.hlines(MSY, 0, e_msy, ls='--', color='black', lw=1) +ax.annotate(f'MSY = {MSY:.0f}', xy=(e_msy, MSY), xytext=(e_msy + 1, MSY - 18), fontsize=11, color='black') -ax.annotate(r'$E_{MSY}=r/2q$', xy=(E_msy, 0), xytext=(E_msy + 1, 6), +ax.annotate(r'$e_{MSY}=r/2q$', xy=(e_msy, 0), xytext=(e_msy + 1, 6), fontsize=11, color='black') -ax.fill_between(E_grid, 0, Y_grid, where=(E_grid > E_msy), color='C1', alpha=0.15) -ax.text(E_msy * 1.45, MSY * 0.35, 'overfishing\n(yield falls)', +ax.fill_between(e_grid, 0, y_grid, where=(e_grid > e_msy), color='C1', alpha=0.15) +ax.text(e_msy * 1.45, MSY * 0.35, 'overfishing\n(yield falls)', color='C1', fontsize=10, ha='center') -ax.set_xlabel('fishing effort $E$') -ax.set_ylabel('sustainable yield $Y$ (tonnes/year)') +ax.set_xlabel('fishing effort $e$') +ax.set_ylabel('sustainable yield $y^*$ (tonnes/year)') ax.set_title('Yield-effort curve (Schaefer model)') ax.set_xlim(0, r / q) ax.set_ylim(0, MSY * 1.25) @@ -411,24 +455,24 @@ Finally, let's run the yearly recursion {eq}`eq:update` forward from several starting stocks, under the MSY effort. ```{code-cell} ipython3 -def simulate(x0, E, years=40): +def simulate(x0, e, years=40): "Run the deterministic yearly stock recursion forward from x0." x = np.empty(years + 1) x[0] = x0 for t in range(years): - x[t + 1] = update(x[t], E) + x[t + 1] = update(x[t], e) return np.arange(years + 1), x fig, ax = plt.subplots(figsize=(8, 5)) for x0 in [100, 300, 700, 950]: - t, xt = simulate(x0, E_msy) + t, xt = simulate(x0, e_msy) ax.plot(t, xt, 'o-', ms=3, lw=1.2, label=f'$x_0={x0}$') ax.axhline(x_msy, ls='--', color='black', lw=1.5) ax.annotate(r'$x^*=K/2$', xy=(0, x_msy), xytext=(1, x_msy + 25), color='black') ax.set_xlabel('year $t$') ax.set_ylabel('stock biomass $x_t$ (tonnes)') -ax.set_title(r'Year-by-year stock path under MSY effort $E=E_{MSY}$') +ax.set_title(r'Year-by-year stock path under MSY effort $e=e_{MSY}$') ax.set_xlim(0, 40) ax.legend(frameon=False) ax.spines[['top', 'right']].set_visible(False) @@ -464,27 +508,27 @@ Shrink the time step to zero and the yearly recursion {eq}`eq:update` becomes th differential equation $$ - \dot{x} = G(x) - qEx. + \dot{x} = G(x) - qex. $$ The sustainability condition is unchanged --- a steady state still needs -$G(x) = qEx$ --- so $x^*(E) = K(1 - qE/r)$ and $Y(E) = qEK(1 - qE/r)$ exactly as +$G(x) = qex$ --- so $x^*(e) = K(1 - qe/r)$ and $y^*(e) = qeK(1 - qe/r)$ exactly as before. -Now $Y$ is a smooth, concave parabola in $E$, and its peak is where the slope +Now $y^*$ is a smooth, concave parabola in $e$, and its peak is where the slope vanishes: $$ - Y'(E) = qK\left(1 - \frac{2qE}{r}\right) = 0 + \frac{d y^*}{d e} = qK\left(1 - \frac{2qe}{r}\right) = 0 \quad\Longrightarrow\quad - E_{MSY} = \frac{r}{2q}, + e_{MSY} = \frac{r}{2q}, $$ which gives $x^* = K/2$ and $\text{MSY} = rK/4$ --- exactly the values the grid search found. -Equivalently, since the sustainable catch equals growth, $Y = G(x^*)$, the MSY is -just the top of the growth curve, where $G'(x) = 0$ at $x = K/2$. +Equivalently, since the sustainable catch equals growth, $y^* = G(x^*)$, the MSY +is just the top of the growth curve, where $G'(x) = 0$ at $x = K/2$. In words: **fish the stock down to where it grows fastest, and no further.** @@ -539,7 +583,7 @@ How should we set the catch $h_t$? We compare two policies, both designed to take the MSY *on average*: -* **constant effort**: set $h_t = q E_{MSY} x_t$, so the catch scales with the +* **constant effort**: set $h_t = q\, e_{MSY}\, x_t$, so the catch scales with the current stock, and * **constant quota**: set $h_t = \text{MSY}$ every year, regardless of the stock. @@ -559,7 +603,7 @@ def simulate_stochastic(policy, σ, years, rng, x0=x_msy): for t in range(years): growth = ξ[t] * G(x[t]) if policy == 'effort': - harvest = q * E_msy * x[t] + harvest = q * e_msy * x[t] else: # constant quota harvest = MSY x[t + 1] = max(0.0, x[t] + growth - harvest) @@ -578,7 +622,7 @@ fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True) for ax, policy, title in zip( axes, ['effort', 'quota'], - ['constant effort $h_t = qE_{MSY}x_t$', + ['constant effort $h_t = q e_{MSY} x_t$', 'constant quota $h_t = \\mathrm{MSY}$']): for k in range(8): path = simulate_stochastic(policy, σ, years, rng) @@ -836,7 +880,7 @@ The key formulas of the deterministic Schaefer model are collected below. |---|---|---| | Stock at MSY | $x_{MSY} = K/2$ | 500 t | | Maximum sustainable yield | $\text{MSY} = rK/4$ | 125 t/yr | -| Effort at MSY | $E_{MSY} = r/2q$ | 25 | +| Effort at MSY | $e_{MSY} = r/2q$ | 25 | The MSY is a deterministic, single-species, steady-state concept. From aaf84fc51c5470cb37350395661c5d9787926456 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 10:54:00 +1000 Subject: [PATCH 3/9] Implement TODOs in revised first half of MSY lecture - Add a figure showing the shape of the logistic growth curve G(x); move the G definition above it. - Refactor plot_45 to take the one-year map as a function argument, so the no-fishing 45-degree diagram works before fishing/effort and the update() function are introduced; update both call sites to pass lambdas. - Add the two-effort 45-degree diagram (staircases omitted) showing x*(e) falling with effort. - Fix the multi-effort growth-vs-harvest figure: use explicit efforts [12.5, 25, 37.5] and G(K/2) for the y-limit (e_msy/MSY are not yet defined there), and fix the label typo. - Restore the grid-search computation of e_msy, MSY, x_msy in the MSY section. - Drop the under/over-fishing framing from the yield-effort curve per review note, since under randomness no constant policy is fully safe. Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 342 +++++++++++++++++++--------------------- 1 file changed, 163 insertions(+), 179 deletions(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index f2b3518f2..66d1b315d 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -18,18 +18,14 @@ kernelspec: In this lecture we study one of the central ideas in the management of renewable resources: the **maximum sustainable yield** (MSY). -The MSY is the largest catch a fish stock can support year after year without +The MSY is the largest catch a fishery can support year after year without running itself down. -We build the idea up with almost no mathematics. +We begin with a relatively elementary discrete time treatment. -The model is just a year-to-year bookkeeping rule, and the best catch can be -found by reading the peak off a graph. +We then discuss a continuous time formulation that conveys the same ideas via calculus. -A short optional section shows how calculus recovers the same answers in closed -form, for readers who want it. - -But the MSY also has a darker side. +Next, we turn to problems associated with MSY-based fishing policy. After the Second World War, fisheries managers around the world adopted MSY (and close relatives) as a target. @@ -37,13 +33,15 @@ close relatives) as a target. Several great fisheries then collapsed --- among them the Peruvian anchovy in the 1970s and the Atlantic cod off Newfoundland in 1992. -A major reason is that the simple MSY model ignores **randomness and risk**. +A major reason is that the simple MSY model ignores randomness and risk. + +One issue is that fish stocks are difficult to track, so policies may be based +on incorrect measurements. -The ocean is a noisy place, and a policy that looks safe in a deterministic model -can be dangerously fragile once we allow the environment to fluctuate. +Another is that ocean environments are complex and nonstationary: a policy that looks safe in a deterministic model +can be dangerously fragile once we admit randomness. -So in the second half of the lecture we add random shocks to the model and watch a -fishery collapse. +To illustrate these ideas, we add random shocks to the model and see how collapse can easily occur. The MSY framework is due to {cite}`schaefer1954` and {cite}`gordon1954`; a classic textbook treatment is {cite}`clark1990`. @@ -55,7 +53,7 @@ import numpy as np import matplotlib.pyplot as plt ``` -We will use a standard parameterization throughout. +We will use the following parameterization throughout. ```{code-cell} ipython3 r = 0.5 # intrinsic growth rate (per year) @@ -73,38 +71,47 @@ fishing does to it. ### Growth without fishing -Left alone, the population adds new fish each year according to **logistic -growth**, +In the model, the population adds new fish each year according the **logistic +growth** function $$ - G(x) = r\,x\left(1 - \frac{x}{K}\right), + G(x) = r\,x\left(1 - \frac{x}{K}\right) $$ (eq:logistic) -with **intrinsic growth rate** $r > 0$ and **carrying capacity** $K > 0$. - -Growth is small when the stock is small (few spawners) and small again when the -stock is near $K$ (crowding, limited food), and is largest in between. - -With no fishing, next year's stock is just this year's stock plus that growth: - -$$ - x_{t+1} = x_t + G(x_t). -$$ +Here $r > 0$ is called **intrinsic growth rate** and $K > 0$ is called the **carrying capacity**. -Let's encode the growth function and the one-year update rule. +Let's encode the growth function in Python. ```{code-cell} ipython3 def G(x): - "Logistic growth added over one year." + "Logistic growth over one year." return r * x * (1 - x / K) +``` -def update(x, e): - "Next year's stock given current stock x and fishing effort e." - return x + G(x) - q * e * x +As the next figure shows, growth is small when the stock is small (few spawners) +and small again when the stock is near $K$ (crowding, limited food), and is +largest in between. + +```{code-cell} ipython3 +x_grid = np.linspace(0, K, 400) + +fig, ax = plt.subplots(figsize=(8, 5)) +ax.plot(x_grid, G(x_grid), lw=2.5, color='C0') +ax.set_xlabel('stock biomass $x$ (tonnes)') +ax.set_ylabel('annual growth $G(x)$ (tonnes/year)') +ax.set_title('Logistic growth') +ax.set_xlim(0, K) +ax.set_ylim(0, G(K / 2) * 1.15) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() ``` -(The `update` function already includes fishing through the effort $e$, which we -introduce below; setting $e = 0$ gives the unfished dynamics.) +With no fishing, next year's stock is just this year's stock plus current growth: + +$$ + x_{t+1} = x_t + G(x_t). +$$ A clean way to see where the dynamics lead is a **45-degree diagram**: we plot next year's stock $x_{t+1}$ against this year's stock $x_t$. @@ -118,18 +125,21 @@ becomes this year's stock), and repeat. The next function draws such a diagram. +It takes the one-year update rule as a function argument `update_fn`, since at +this point we have not yet introduced fishing. + ```{code-cell} ipython3 -def plot_45(ax, e, x0, x_max, steady_state, ss_label, map_label, n_years=30): - "Draw a 45-degree (cobweb) diagram for the yearly stock update at effort e." +def plot_45(ax, update_fn, x0, x_max, steady_state, ss_label, map_label, n_years=30): + "Draw a 45-degree (cobweb) diagram for a one-year stock update rule." grid = np.linspace(0, x_max, 400) - ax.plot(grid, update(grid, e), color='C0', lw=2.5, label=map_label) + ax.plot(grid, update_fn(grid), color='C0', lw=2.5, label=map_label) ax.plot(grid, grid, color='0.6', lw=1, ls='--', label=r'$45^\circ$ line $x_{t+1}=x_t$') # cobweb staircase starting from x0 x = x0 cx, cy = [x], [0.0] for _ in range(n_years): - y = update(x, e) + y = update_fn(x) cx += [x, y] cy += [y, y] x = y @@ -151,7 +161,7 @@ Here is the unfished case. ```{code-cell} ipython3 fig, ax = plt.subplots(figsize=(4.95, 4.95)) -plot_45(ax, e=0, x0=80, x_max=1100, +plot_45(ax, lambda x: x + G(x), x0=80, x_max=1100, steady_state=K, ss_label=r'$x=K$ (carrying capacity)', map_label=r'$x_{t+1}=x_t+G(x_t)$') ax.set_title('Stock dynamics without fishing') @@ -167,6 +177,7 @@ The unfished population fills up the environment and stays there. (The empty ocean $x = 0$ is also a steady state, but an unstable one --- any small stock grows away from it.) + ### Adding fishing Now let a fishing fleet remove a catch each year. @@ -189,12 +200,16 @@ $$ \;-\; \underbrace{qe\,x_t}_{\text{catch}}. $$ (eq:update) -That is the whole model --- a rule you could run forward in a spreadsheet. -On the 45-degree diagram, fishing pulls the update curve **down** by the harvest -term, so it now crosses the $45^\circ$ line at a **lower** steady state. +Here's Python code to implement this update rule -A fished population settles below its unfished level. +```{code-cell} ipython3 +def update(x, e): + "Next year's stock given current stock x and fishing effort e." + return x + G(x) - q * e * x +``` + +Here's the 45 degree diagram, now with the fishing term included. ```{code-cell} ipython3 e_demo = 20.0 # an illustrative fixed effort level @@ -204,7 +219,7 @@ def x_star(e): return K * (1 - q * e / r) fig, ax = plt.subplots(figsize=(4.95, 4.95)) -plot_45(ax, e=e_demo, x0=80, x_max=1100, +plot_45(ax, lambda x: update(x, e_demo), x0=80, x_max=1100, steady_state=x_star(e_demo), ss_label=r'$x^*(e)$', map_label=r'$x_{t+1}=x_t+G(x_t)-qex_t$') ax.set_title(f'Stock dynamics at a fixed effort $e={e_demo:.0f}$') @@ -212,26 +227,51 @@ plt.tight_layout() plt.show() ``` -The natural question is now: which effort level $e$ gives the largest catch we can -take year after year? +On the 45-degree diagram, fishing pulls the update curve **down** by the harvest +term, so it now crosses the $45^\circ$ line at a **lower** steady state. + +Since this steady state is a function of $e$ now, we denote it by $x^*(e)$. -## When is a catch sustainable? +Here's the dynamics for two different levels of $e$, with the staircases omitted. + +```{code-cell} ipython3 +grid = np.linspace(0, 1100, 400) -Suppose the fleet applies a constant effort $e$ every year. +fig, ax = plt.subplots(figsize=(4.95, 4.95)) +ax.plot(grid, grid, color='0.6', lw=1, ls='--', label=r'$45^\circ$ line') +for e, c in zip((10.0, 30.0), ('C0', 'C3')): + ax.plot(grid, update(grid, e), lw=2.5, color=c, label=f'$e={e:.0f}$') + xs = x_star(e) + ax.plot([xs], [xs], 'o', color=c, ms=8, zorder=5) + +ax.set_xlabel('stock this year $x_t$ (tonnes)') +ax.set_ylabel('stock next year $x_{t+1}$ (tonnes)') +ax.set_title('Stock dynamics at two effort levels') +ax.set_xlim(0, 1100) +ax.set_ylim(0, 1100) +ax.set_aspect('equal') +ax.legend(loc='upper left', frameon=False, fontsize=9) +ax.spines[['top', 'right']].set_visible(False) +plt.tight_layout() +plt.show() +``` -The catch is **sustainable** if the stock holds steady from one year to the next ---- it neither grows without limit nor dwindles away. +Not surprisingly, the steady state $x^*(e)$ is decreasing in fishing effort. -In symbols, the stock repeats itself: -$$ - x_{t+1} = x_t. -$$ + +## When is a catch sustainable? + +Suppose, as above, that the fleet applies a constant effort $e$ every year. + +Let's look for the associated steady state. + +At any steady state, we must have $x_{t+1} = x_t$. Looking at {eq}`eq:update`, this happens exactly when growth replaces the catch, $G(x) = qex$. -Writing it out and cancelling the empty-ocean case $x = 0$: +Writing it out and ignoring the empty-ocean case $x = 0$, we get $$ r\left(1 - \frac{x}{K}\right) = qe @@ -239,16 +279,13 @@ $$ x^*(e) = K\left(1 - \frac{qe}{r}\right). $$ (eq:xstar) -So **each effort level $e$ pins down one steady-state stock** $x^*(e)$ --- -provided $e < r/q$. -Any more effort than that drives the stock to zero. +Provided $e < r/q$, we see that each effort level $e$ pins down one steady-state stock $x^*(e)$. -The catch this steady state delivers, year after year, is the **sustainable -yield** +The catch this steady state delivers, year after year, is the **sustainable yield**, defined as $$ - y^*(e) = q\,e\,x^*(e). + y^*(e) := q\,e\,x^*(e). $$ (eq:yield) ```{code-cell} ipython3 @@ -257,11 +294,8 @@ def sustainable_yield(e): return q * e * x_star(e) ``` -We can read both of these quantities --- the steady-state stock $x^*(e)$ and the -sustainable catch $y^*(e)$ --- straight off a picture. - -Plot the growth curve $G(x)$ together with the harvest line $q e x$ for a single -effort level. +To visualize the sustainable yield, we plot the growth curve $G(x)$ together +with the harvest line $q e x$ for a single effort level $e$. The two cross where growth exactly replaces the catch: that crossing sits at the steady-state stock $x^*(e)$. @@ -294,102 +328,17 @@ plt.tight_layout() plt.show() ``` -As we change $e$, the harvest line pivots about the origin, the crossing slides -along the growth curve, and the sustainable catch $y^*(e)$ is always the height of -that crossing. - -Notice that finding all this needed only algebra --- no calculus. - -## The maximum sustainable yield - Different effort levels give different sustainable catches $y^*(e)$. -The **maximum sustainable yield** is simply the largest of them: - -$$ - \text{MSY} \;=\; \max_{0 \le e < r/q} \; y^*(e). -$$ - -We can find it without any calculus: evaluate $y^*(e)$ on a fine grid of effort -levels and pick the biggest. - -```{code-cell} ipython3 -e_vals = np.linspace(0, r / q, 100_001) # efforts from 0 up to the wipe-out level r/q -y_vals = sustainable_yield(e_vals) -best = np.argmax(y_vals) - -e_msy = e_vals[best] -MSY = y_vals[best] -x_msy = x_star(e_msy) - -print(f"Best effort e_MSY = {e_msy:.2f}") -print(f"Sustainable stock x* = {x_msy:.1f} tonnes (= K/2)") -print(f"Max sustainable yield MSY = {MSY:.1f} tonnes/year (= rK/4)") -``` - -The grid search lands on a sustainable stock of $x^* = 500 = K/2$ and a maximum -catch of $\text{MSY} = 125 = rK/4$. - -Those round numbers are no accident --- the optional calculus section below shows -where they come from. - -## Visualizing the MSY - -### The sustainable-yield curve - -As effort climbs from $0$ toward $r/q$, the sustainable stock $x^*(e)$ slides from -$K$ down to $0$, and the catch it supports traces out the growth curve $G(x)$. - -So the hump below is the full menu of sustainable (stock, catch) pairs --- and its -peak is the MSY. - -```{code-cell} ipython3 -x = np.linspace(0, K, 400) - -fig, ax = plt.subplots(figsize=(8, 5)) -ax.plot(x, G(x), lw=2.5, color='C0', - label=r'sustainable yield $y^*=G(x)=rx(1-x/K)$') - -ax.plot([x_msy], [MSY], 'o', color='black', ms=9, zorder=5) -ax.vlines(x_msy, 0, MSY, ls='--', color='black', lw=1) -ax.hlines(MSY, 0, x_msy, ls='--', color='black', lw=1) -ax.annotate(f'MSY = rK/4 = {MSY:.0f}', - xy=(x_msy, MSY), xytext=(x_msy + 60, MSY - 18), - fontsize=11, color='black') -ax.annotate(r'$x_{MSY}=K/2$', xy=(x_msy, 0), xytext=(x_msy + 10, 6), - fontsize=11, color='black') - -ax.set_xlabel('stock biomass $x$ (tonnes)') -ax.set_ylabel('growth / sustainable yield (tonnes/year)') -ax.set_title('Maximum sustainable yield: the peak of the growth curve') -ax.set_xlim(0, K) -ax.set_ylim(0, MSY * 1.25) -ax.legend(loc='upper right', frameon=False) -ax.spines[['top', 'right']].set_visible(False) -plt.tight_layout() -plt.show() -``` - -### Growth versus harvest - -We just saw how a single harvest line picks out one sustainable point. - -To compare effort levels, plot the growth curve $G(x)$ together with the harvest +In the next figure we plot the growth curve $G(x)$ together with the harvest lines $q e x$ for several values of $e$. -Each crossing is a sustainable steady state $x^*(e)$, and the catch there is the -height of the line. - -The MSY effort line meets the growth curve exactly at its peak, $x = K/2$. - ```{code-cell} ipython3 fig, ax = plt.subplots(figsize=(8, 5)) ax.plot(x, G(x), lw=2.5, color='C0', label='growth $G(x)$') -efforts = [0.5 * e_msy, e_msy, 1.5 * e_msy] -labels = [r'$e < e_{MSY}$ (underfishing)', - r'$e = e_{MSY}$', - r'$e > e_{MSY}$ (overfishing)'] +efforts = [12.5, 25.0, 37.5] +labels = [r'low $e$', r'moderate $e$', r'high $e$'] colors = ['C2', 'C3', 'C1'] for e, lab, c in zip(efforts, labels, colors): @@ -398,27 +347,58 @@ for e, lab, c in zip(efforts, labels, colors): ax.set_xlabel('stock biomass $x$ (tonnes)') ax.set_ylabel('rate (tonnes/year)') -ax.set_title('Sustainable steady states: growth $G(x)$ vs. harvest $qex$') +ax.set_title('Steady states: growth $G(x)$ vs. harvest $qex$') ax.set_xlim(0, K) -ax.set_ylim(0, MSY * 1.4) +ax.set_ylim(0, G(K / 2) * 1.4) ax.legend(loc='upper left', frameon=False, fontsize=10) ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` -Notice that the under- and over-fishing lines cross the growth curve at the -**same height** (same catch) but at different stock levels. +Each dot marks a steady state, and its height is the sustainable catch at that +effort. + +As effort rises, the steady-state stock falls, while the catch rises and then +falls --- exactly the trade-off that the maximum sustainable yield captures. + + + +## The maximum sustainable yield + +The **maximum sustainable yield** is defined by -Only $e_{MSY}$ reaches the very top. +$$ + \text{MSY} \;=\; \max_{0 \le e < r/q} \; y^*(e). +$$ + +This is the largest steady state catch attainable, assuming a constant effort +rate $e$. -Increasing effort past $e_{MSY}$ buys a *smaller* sustained catch from a *more -depleted* stock --- biologically and economically the worst of both worlds. +We can compute it numerically by evaluating $y^*(e)$ on a fine grid of effort +levels and picking the largest. + +```{code-cell} ipython3 +e_search = np.linspace(0, r / q, 100_001) +i = np.argmax(sustainable_yield(e_search)) + +e_msy = e_search[i] +MSY = sustainable_yield(e_msy) +x_msy = x_star(e_msy) + +print(f"effort at MSY e_MSY = {e_msy:.2f}") +print(f"stock at MSY x* = {x_msy:.1f} tonnes (= K/2)") +print(f"maximum yield MSY = {MSY:.1f} tonnes/year (= rK/4)") +``` + +The search lands on the round numbers $x^* = K/2$ and $\text{MSY} = rK/4$; the +optional calculus section below shows why. ### The yield-effort curve -Plotting the sustainable catch $y^*(e)$ directly against effort gives the classic -dome-shaped Schaefer curve that fisheries managers use. +To visualize the MSY, we plot the sustainable catch $y^*(e)$ against effort $e$. + +This gives the classic dome-shaped Schaefer curve that fisheries managers use. ```{code-cell} ipython3 e_grid = np.linspace(0, r / q, 400) @@ -434,10 +414,6 @@ ax.annotate(f'MSY = {MSY:.0f}', xy=(e_msy, MSY), xytext=(e_msy + 1, MSY - 18), ax.annotate(r'$e_{MSY}=r/2q$', xy=(e_msy, 0), xytext=(e_msy + 1, 6), fontsize=11, color='black') -ax.fill_between(e_grid, 0, y_grid, where=(e_grid > e_msy), color='C1', alpha=0.15) -ax.text(e_msy * 1.45, MSY * 0.35, 'overfishing\n(yield falls)', - color='C1', fontsize=10, ha='center') - ax.set_xlabel('fishing effort $e$') ax.set_ylabel('sustainable yield $y^*$ (tonnes/year)') ax.set_title('Yield-effort curve (Schaefer model)') @@ -449,6 +425,17 @@ plt.tight_layout() plt.show() ``` +The sustainable yield rises with effort up to the peak at $e_{MSY}$ and then +falls. + +Pushing effort beyond $e_{MSY}$ is doubly costly: the catch is smaller *and* the +stock left in the ocean is smaller. + +It is tempting to call any effort below $e_{MSY}$ "safe", but we resist that +label: once randomness enters the picture, as it does in the next section, no +constant policy is entirely safe. + + ### Does the stock settle there? Finally, let's run the yearly recursion {eq}`eq:update` forward from several @@ -497,14 +484,16 @@ $r > 0$, which is one reason it is the textbook standard. ## The same result with calculus -```{note} -This section is optional. -Everything above used only algebra and a graph; readers who have met calculus can -use it to pin down the round numbers exactly. -``` +The following section is optional. + +Our discussion above used only algebra and figures. -Shrink the time step to zero and the yearly recursion {eq}`eq:update` becomes the +For readers who are familiar with calculus, we now provide a continuous time version of the model. + +The continuous time version is frequently used and has some advantages (and some disadvantages). + +If we shrink the time step to zero, the recursion {eq}`eq:update` becomes the differential equation $$ @@ -524,13 +513,8 @@ $$ e_{MSY} = \frac{r}{2q}, $$ -which gives $x^* = K/2$ and $\text{MSY} = rK/4$ --- exactly the values the grid -search found. - -Equivalently, since the sustainable catch equals growth, $y^* = G(x^*)$, the MSY -is just the top of the growth curve, where $G'(x) = 0$ at $x = K/2$. +which gives $x^* = K/2$ and $\text{MSY} = rK/4$. -In words: **fish the stock down to where it grows fastest, and no further.** ```{note} Readers who have seen the Solow-Swan growth model may recognize the structure. From cd0f547d56d5a739affa92b47fe5a0fa15dfab63 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 11:49:31 +1000 Subject: [PATCH 4/9] Align MSY lecture with QuantEcon style guide; fix two typos Following manual/styleguide (figures.md, math.md, references.md): - Remove all embedded matplotlib titles (ax.set_title / fig.suptitle) and add mystnb figure captions + names to every figure cell. - Keep the figure box: drop all ax.spines[...].set_visible(False). - Standardize data-line widths to lw=2; drop figsize from standard single-panel charts (keep square figsize on the cobweb diagrams and the stacked stochastic panel as justified exceptions). - Label the two policy panels with in-axes text rather than titles. - Use {cite:t} for in-sentence citations (Schaefer/Gordon/Clark, Following Schaefer, analyzed by May); keep {cite} for parenthetical. - Fix typos: "according the" -> "according to the"; "called intrinsic" -> "called the intrinsic". Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 178 +++++++++++++++++++++++++--------------- 1 file changed, 110 insertions(+), 68 deletions(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index 66d1b315d..1f9594213 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -43,8 +43,8 @@ can be dangerously fragile once we admit randomness. To illustrate these ideas, we add random shocks to the model and see how collapse can easily occur. -The MSY framework is due to {cite}`schaefer1954` and {cite}`gordon1954`; a -classic textbook treatment is {cite}`clark1990`. +The MSY framework is due to {cite:t}`schaefer1954` and {cite:t}`gordon1954`; a +classic textbook treatment is {cite:t}`clark1990`. Let's start with some standard imports. @@ -71,14 +71,14 @@ fishing does to it. ### Growth without fishing -In the model, the population adds new fish each year according the **logistic +In the model, the population adds new fish each year according to the **logistic growth** function $$ G(x) = r\,x\left(1 - \frac{x}{K}\right) $$ (eq:logistic) -Here $r > 0$ is called **intrinsic growth rate** and $K > 0$ is called the **carrying capacity**. +Here $r > 0$ is called the **intrinsic growth rate** and $K > 0$ is called the **carrying capacity**. Let's encode the growth function in Python. @@ -93,16 +93,20 @@ and small again when the stock is near $K$ (crowding, limited food), and is largest in between. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Logistic growth + name: fig:logistic-growth +--- x_grid = np.linspace(0, K, 400) -fig, ax = plt.subplots(figsize=(8, 5)) -ax.plot(x_grid, G(x_grid), lw=2.5, color='C0') -ax.set_xlabel('stock biomass $x$ (tonnes)') -ax.set_ylabel('annual growth $G(x)$ (tonnes/year)') -ax.set_title('Logistic growth') +fig, ax = plt.subplots() +ax.plot(x_grid, G(x_grid), lw=2, color='C0') +ax.set_xlabel('stock biomass $x$') +ax.set_ylabel('annual growth $G(x)$') ax.set_xlim(0, K) ax.set_ylim(0, G(K / 2) * 1.15) -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -132,9 +136,9 @@ this point we have not yet introduced fishing. def plot_45(ax, update_fn, x0, x_max, steady_state, ss_label, map_label, n_years=30): "Draw a 45-degree (cobweb) diagram for a one-year stock update rule." grid = np.linspace(0, x_max, 400) - ax.plot(grid, update_fn(grid), color='C0', lw=2.5, label=map_label) + ax.plot(grid, update_fn(grid), color='C0', lw=2, label=map_label) ax.plot(grid, grid, color='0.6', lw=1, ls='--', - label=r'$45^\circ$ line $x_{t+1}=x_t$') + label=r'$45^\circ$ line') # cobweb staircase starting from x0 x = x0 cx, cy = [x], [0.0] @@ -144,27 +148,31 @@ def plot_45(ax, update_fn, x0, x_max, steady_state, ss_label, map_label, n_years cy += [y, y] x = y ax.plot(cx, cy, color='black', lw=1, alpha=0.9) - ax.plot([steady_state], [steady_state], 'o', color='black', ms=8, zorder=5) + ax.plot([steady_state], [steady_state], 'o', color='black', ms=6, zorder=5) ax.annotate(ss_label, xy=(steady_state, steady_state), xytext=(0.55 * x_max, 0.28 * x_max), fontsize=11, arrowprops=dict(arrowstyle='->', color='black', lw=1)) - ax.set_xlabel('stock this year $x_t$ (tonnes)') - ax.set_ylabel('stock next year $x_{t+1}$ (tonnes)') + ax.set_xlabel('stock this year $x_t$') + ax.set_ylabel('stock next year $x_{t+1}$') ax.set_xlim(0, x_max) ax.set_ylim(0, x_max) ax.set_aspect('equal') ax.legend(loc='upper left', frameon=False, fontsize=9) - ax.spines[['top', 'right']].set_visible(False) ``` Here is the unfished case. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Stock dynamics without fishing + name: fig:cobweb-unfished +--- fig, ax = plt.subplots(figsize=(4.95, 4.95)) plot_45(ax, lambda x: x + G(x), x0=80, x_max=1100, steady_state=K, ss_label=r'$x=K$ (carrying capacity)', map_label=r'$x_{t+1}=x_t+G(x_t)$') -ax.set_title('Stock dynamics without fishing') plt.tight_layout() plt.show() ``` @@ -182,7 +190,7 @@ stock grows away from it.) Now let a fishing fleet remove a catch each year. -Following {cite}`schaefer1954`, the catch is proportional to fishing **effort** +Following {cite:t}`schaefer1954`, the catch is proportional to fishing **effort** $e$ (e.g. boat-days) and to the stock available to be caught: $$ @@ -212,6 +220,12 @@ def update(x, e): Here's the 45 degree diagram, now with the fishing term included. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Stock dynamics at a fixed effort + name: fig:cobweb-fished +--- e_demo = 20.0 # an illustrative fixed effort level def x_star(e): @@ -222,7 +236,6 @@ fig, ax = plt.subplots(figsize=(4.95, 4.95)) plot_45(ax, lambda x: update(x, e_demo), x0=80, x_max=1100, steady_state=x_star(e_demo), ss_label=r'$x^*(e)$', map_label=r'$x_{t+1}=x_t+G(x_t)-qex_t$') -ax.set_title(f'Stock dynamics at a fixed effort $e={e_demo:.0f}$') plt.tight_layout() plt.show() ``` @@ -235,23 +248,27 @@ Since this steady state is a function of $e$ now, we denote it by $x^*(e)$. Here's the dynamics for two different levels of $e$, with the staircases omitted. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Stock dynamics at two effort levels + name: fig:cobweb-two-efforts +--- grid = np.linspace(0, 1100, 400) fig, ax = plt.subplots(figsize=(4.95, 4.95)) ax.plot(grid, grid, color='0.6', lw=1, ls='--', label=r'$45^\circ$ line') for e, c in zip((10.0, 30.0), ('C0', 'C3')): - ax.plot(grid, update(grid, e), lw=2.5, color=c, label=f'$e={e:.0f}$') + ax.plot(grid, update(grid, e), lw=2, color=c, label=f'$e={e:.0f}$') xs = x_star(e) ax.plot([xs], [xs], 'o', color=c, ms=8, zorder=5) -ax.set_xlabel('stock this year $x_t$ (tonnes)') -ax.set_ylabel('stock next year $x_{t+1}$ (tonnes)') -ax.set_title('Stock dynamics at two effort levels') +ax.set_xlabel('stock this year $x_t$') +ax.set_ylabel('stock next year $x_{t+1}$') ax.set_xlim(0, 1100) ax.set_ylim(0, 1100) ax.set_aspect('equal') ax.legend(loc='upper left', frameon=False, fontsize=9) -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -260,7 +277,7 @@ Not surprisingly, the steady state $x^*(e)$ is decreasing in fishing effort. -## When is a catch sustainable? +## Sustainable yield Suppose, as above, that the fleet applies a constant effort $e$ every year. @@ -303,10 +320,16 @@ steady-state stock $x^*(e)$. The height of the line there is the sustainable catch $y^*(e)$. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Sustainable stock and catch + name: fig:sustainable-point +--- x = np.linspace(0, K, 400) -fig, ax = plt.subplots(figsize=(8, 5)) -ax.plot(x, G(x), lw=2.5, color='C0', label=r'growth $G(x)$') +fig, ax = plt.subplots() +ax.plot(x, G(x), lw=2, color='C0', label=r'growth $G(x)$') ax.plot(x, q * e_demo * x, lw=2, color='C3', label=r'harvest $q e x$') xs = x_star(e_demo) @@ -317,13 +340,11 @@ ax.hlines(ys, 0, xs, ls='--', color='black', lw=1) ax.annotate(r'$x^*(e)$', xy=(xs, 0), xytext=(xs + 12, 8), fontsize=12) ax.annotate(r'$y^*(e)$', xy=(0, ys), xytext=(10, ys + 5), fontsize=12) -ax.set_xlabel('stock biomass $x$ (tonnes)') -ax.set_ylabel('rate (tonnes/year)') -ax.set_title(f'Sustainable stock and catch at effort $e={e_demo:.0f}$') +ax.set_xlabel('stock biomass $x$') +ax.set_ylabel('catch') ax.set_xlim(0, K) ax.set_ylim(0, G(K / 2) * 1.4) -ax.legend(loc='upper right', frameon=False) -ax.spines[['top', 'right']].set_visible(False) +ax.legend(loc='upper left', frameon=False) plt.tight_layout() plt.show() ``` @@ -334,24 +355,28 @@ In the next figure we plot the growth curve $G(x)$ together with the harvest lines $q e x$ for several values of $e$. ```{code-cell} ipython3 -fig, ax = plt.subplots(figsize=(8, 5)) -ax.plot(x, G(x), lw=2.5, color='C0', label='growth $G(x)$') +--- +mystnb: + figure: + caption: Steady states at several effort levels + name: fig:steady-states +--- +fig, ax = plt.subplots() +ax.plot(x, G(x), lw=2, color='C0', label='growth $G(x)$') efforts = [12.5, 25.0, 37.5] labels = [r'low $e$', r'moderate $e$', r'high $e$'] colors = ['C2', 'C3', 'C1'] for e, lab, c in zip(efforts, labels, colors): - ax.plot(x, q * e * x, lw=1.8, color=c, label=lab) + ax.plot(x, q * e * x, lw=2, color=c, label=lab) ax.plot([x_star(e)], [q * e * x_star(e)], 'o', color=c, ms=7, zorder=5) -ax.set_xlabel('stock biomass $x$ (tonnes)') -ax.set_ylabel('rate (tonnes/year)') -ax.set_title('Steady states: growth $G(x)$ vs. harvest $qex$') +ax.set_xlabel('stock $x$') +ax.set_ylabel('catch') ax.set_xlim(0, K) ax.set_ylim(0, G(K / 2) * 1.4) ax.legend(loc='upper left', frameon=False, fontsize=10) -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -401,26 +426,30 @@ To visualize the MSY, we plot the sustainable catch $y^*(e)$ against effort $e$. This gives the classic dome-shaped Schaefer curve that fisheries managers use. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Yield-effort curve + name: fig:yield-effort +--- e_grid = np.linspace(0, r / q, 400) y_grid = sustainable_yield(e_grid) -fig, ax = plt.subplots(figsize=(8, 5)) -ax.plot(e_grid, y_grid, lw=2.5, color='C0', label=r'$y^*(e)=qeK\,(1-qe/r)$') +fig, ax = plt.subplots() +ax.plot(e_grid, y_grid, lw=2, color='C0', label=r'$y^*(e)=qeK\,(1-qe/r)$') ax.plot([e_msy], [MSY], 'o', color='black', ms=9, zorder=5) ax.vlines(e_msy, 0, MSY, ls='--', color='black', lw=1) ax.hlines(MSY, 0, e_msy, ls='--', color='black', lw=1) ax.annotate(f'MSY = {MSY:.0f}', xy=(e_msy, MSY), xytext=(e_msy + 1, MSY - 18), fontsize=11, color='black') -ax.annotate(r'$e_{MSY}=r/2q$', xy=(e_msy, 0), xytext=(e_msy + 1, 6), +ax.annotate(r'$e_{MSY}$', xy=(e_msy, 0), xytext=(e_msy + 1, 6), fontsize=11, color='black') ax.set_xlabel('fishing effort $e$') -ax.set_ylabel('sustainable yield $y^*$ (tonnes/year)') -ax.set_title('Yield-effort curve (Schaefer model)') +ax.set_ylabel('sustainable yield $y^*(e)$') ax.set_xlim(0, r / q) ax.set_ylim(0, MSY * 1.25) ax.legend(loc='upper right', frameon=False) -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -442,6 +471,12 @@ Finally, let's run the yearly recursion {eq}`eq:update` forward from several starting stocks, under the MSY effort. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Stock paths under MSY effort + name: fig:msy-convergence +--- def simulate(x0, e, years=40): "Run the deterministic yearly stock recursion forward from x0." x = np.empty(years + 1) @@ -450,19 +485,17 @@ def simulate(x0, e, years=40): x[t + 1] = update(x[t], e) return np.arange(years + 1), x -fig, ax = plt.subplots(figsize=(8, 5)) +fig, ax = plt.subplots() for x0 in [100, 300, 700, 950]: t, xt = simulate(x0, e_msy) - ax.plot(t, xt, 'o-', ms=3, lw=1.2, label=f'$x_0={x0}$') + ax.plot(t, xt, 'o-', ms=3, lw=2, label=f'$x_0={x0}$') -ax.axhline(x_msy, ls='--', color='black', lw=1.5) +ax.axhline(x_msy, ls='--', color='black', lw=1) ax.annotate(r'$x^*=K/2$', xy=(0, x_msy), xytext=(1, x_msy + 25), color='black') ax.set_xlabel('year $t$') -ax.set_ylabel('stock biomass $x_t$ (tonnes)') -ax.set_title(r'Year-by-year stock path under MSY effort $e=e_{MSY}$') +ax.set_ylabel('stock biomass $x_t$') ax.set_xlim(0, 40) ax.legend(frameon=False) -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -598,30 +631,32 @@ Here are a few sample paths under each policy, with a moderate amount of environmental noise. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: "Two MSY policies under noise: constant effort (top) and constant quota (bottom)" + name: fig:msy-policies +--- σ = 0.15 years = 100 rng = np.random.default_rng(seed=1) -fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True) +fig, axes = plt.subplots(2, 1, figsize=(6, 8), sharey=True) -for ax, policy, title in zip( +for ax, policy, label in zip( axes, ['effort', 'quota'], ['constant effort $h_t = q e_{MSY} x_t$', 'constant quota $h_t = \\mathrm{MSY}$']): for k in range(8): path = simulate_stochastic(policy, σ, years, rng) - ax.plot(path, lw=1, alpha=0.8) + ax.plot(path, lw=2, alpha=0.8) ax.axhline(x_msy, ls='--', color='black', lw=1) - ax.set_title(title) + ax.text(0.5, 0.92, label, transform=ax.transAxes, ha='center', va='top') ax.set_xlabel('year $t$') + ax.set_ylabel('stock biomass $x_t$') ax.set_xlim(0, years) ax.set_ylim(0, K) - ax.spines[['top', 'right']].set_visible(False) -axes[0].set_ylabel('stock biomass $x_t$ (tonnes)') -axes[0].annotate(r'$x^*=K/2$', xy=(years, x_msy), xytext=(-20, 8), - textcoords='offset points', color='black') -fig.suptitle(f'Two MSY policies under environmental noise ($\\sigma={σ}$)') plt.tight_layout() plt.show() ``` @@ -641,7 +676,6 @@ The quota keeps demanding the full MSY even after a run of bad years has thinned the stock, and eventually the catch exceeds what the depleted population can replace. -Once that happens there is no way back. ### Why the quota is a knife-edge @@ -756,6 +790,12 @@ How much does pulling the quota back below the MSY reduce the risk? We adapt `simulate_stochastic` to take an arbitrary fixed quota. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Collapse probability versus quota + name: fig:precautionary-quota +--- def collapse_fraction_quota(quota, σ=0.15, years=100, n_paths=2000, seed=0): "Collapse probability under a fixed quota." rng = np.random.default_rng(seed) @@ -778,8 +818,6 @@ fig, ax = plt.subplots() ax.plot(alphas, probs, 'o-', lw=2) ax.set_xlabel(r'quota as a fraction $\alpha$ of MSY') ax.set_ylabel('collapse probability within 100 years') -ax.set_title('Precaution pays: lower quotas are far safer') -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -821,6 +859,12 @@ Describe what happens as $r$ increases. ``` ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: The yearly logistic map for several values of r + name: fig:logistic-map +--- def logistic_path(r_val, u0=0.3, years=40): u = np.empty(years + 1) u[0] = u0 @@ -828,15 +872,13 @@ def logistic_path(r_val, u0=0.3, years=40): u[t + 1] = u[t] + r_val * u[t] * (1 - u[t]) return u -fig, ax = plt.subplots(figsize=(8, 5)) +fig, ax = plt.subplots() for r_val in [1.5, 2.2, 2.7]: - ax.plot(logistic_path(r_val), 'o-', ms=3, lw=1, label=f'$r={r_val}$') + ax.plot(logistic_path(r_val), 'o-', ms=3, lw=2, label=f'$r={r_val}$') ax.axhline(1.0, ls='--', color='black', lw=1) ax.set_xlabel('year $t$') ax.set_ylabel(r'scaled stock $u_t = x_t / K$') -ax.set_title('The yearly logistic map as $r$ grows') ax.legend(frameon=False) -ax.spines[['top', 'right']].set_visible(False) plt.tight_layout() plt.show() ``` @@ -847,7 +889,7 @@ For $r = 2.2$ it overshoots and settles into a steady oscillation (a two-year cycle). For $r = 2.7$ the path never settles --- it bounces around erratically, the onset -of the chaotic behavior analyzed by {cite}`may1976`. +of the chaotic behavior analyzed by {cite:t}`may1976`. The biology is unchanged; it is the large discrete time step that manufactures these wild dynamics, which is why the smooth continuous-time version is the From 7613242beed387a37b9c2a4bd661bea64b16ea5e Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 11:53:57 +1000 Subject: [PATCH 5/9] Reflow msy-policies figure caption Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index 1f9594213..801526d43 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -634,7 +634,8 @@ environmental noise. --- mystnb: figure: - caption: "Two MSY policies under noise: constant effort (top) and constant quota (bottom)" + caption: 'Two MSY policies under noise: constant effort (top) and constant quota + (bottom)' name: fig:msy-policies --- σ = 0.15 From 5b14da4c792f20066ce7fa9d1b6355979dab7fc4 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 12:09:23 +1000 Subject: [PATCH 6/9] Fix LaTeX PDF build: drop figure floats inside solution blocks The mystnb figure captions on the two code-cells inside the {solution-start}/{solution-end} dropdowns produced LaTeX figure floats nested in a directive, triggering "! LaTeX Error: Not in outer par mode." Per the admonitions style guide (use image, not figure floats, inside exercise/solution directives for PDF builds), remove the caption metadata from those two solution cells; the top-level captioned figures are unaffected. Verified with a local pdflatex build (24-page PDF, no error). Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index 801526d43..e5fb05f0e 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -791,12 +791,6 @@ How much does pulling the quota back below the MSY reduce the risk? We adapt `simulate_stochastic` to take an arbitrary fixed quota. ```{code-cell} ipython3 ---- -mystnb: - figure: - caption: Collapse probability versus quota - name: fig:precautionary-quota ---- def collapse_fraction_quota(quota, σ=0.15, years=100, n_paths=2000, seed=0): "Collapse probability under a fixed quota." rng = np.random.default_rng(seed) @@ -860,12 +854,6 @@ Describe what happens as $r$ increases. ``` ```{code-cell} ipython3 ---- -mystnb: - figure: - caption: The yearly logistic map for several values of r - name: fig:logistic-map ---- def logistic_path(r_val, u0=0.3, years=40): u = np.empty(years + 1) u[0] = u0 From 659f66f1267e2f1d39fa4f21adf548586df8ce83 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 15:29:57 +1000 Subject: [PATCH 7/9] Restructure MSY lecture around quotas; fix yield-effort build break Reframe as "Fishery Management: Quotas and MSY": compute the MSY after introducing the yield-effort dome, add Dynamics and Adding randomness subsections, list x/r/K under the growth equation, and credit the "Maths from the past" discussion. Fix the reorder regression: the yield-effort figure referenced MSY in set_ylim before MSY is computed (NameError); use y_grid.max() instead. Also fix two typos ("sevaral" -> "several"; missing question mark). Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 178 +++++++++++++++++++++------------------- 1 file changed, 95 insertions(+), 83 deletions(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index e5fb05f0e..555f09071 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -11,41 +11,40 @@ kernelspec: name: python3 --- -# Maximum Sustainable Yield and Fishery Collapse +# Fishery Management: Quotas and MSY ## Overview -In this lecture we study one of the central ideas in the management of renewable -resources: the **maximum sustainable yield** (MSY). +In this lecture we study fishery management through quotas and **maximum sustainable yield** (MSY). -The MSY is the largest catch a fishery can support year after year without -running itself down. +The MSY is defined as the largest catch a fishery can support year after year without running itself down. -We begin with a relatively elementary discrete time treatment. +After the Second World War, fisheries managers around the world adopted MSY and +imposed quotas and restrictions based on MSY. -We then discuss a continuous time formulation that conveys the same ideas via calculus. +In some cases these management efforts were successful. -Next, we turn to problems associated with MSY-based fishing policy. +However, this hasn't always been the case. -After the Second World War, fisheries managers around the world adopted MSY (and -close relatives) as a target. +For example, several major fisheries collapsed under MSY-based management, +including the Peruvian anchovy in the 1970s and the Atlantic cod off +Newfoundland in 1992. -Several great fisheries then collapsed --- among them the Peruvian anchovy in the -1970s and the Atlantic cod off Newfoundland in 1992. +In this lecture our main task is to explain how MSY is computed. -A major reason is that the simple MSY model ignores randomness and risk. +We begin with a relatively elementary treatment in discrete time. -One issue is that fish stocks are difficult to track, so policies may be based -on incorrect measurements. +Next we discuss a continuous time formulation that conveys the same ideas via calculus. -Another is that ocean environments are complex and nonstationary: a policy that looks safe in a deterministic model -can be dangerously fragile once we admit randomness. +Finally, we turn to problems associated with MSY and quotas derived from MSY. To illustrate these ideas, we add random shocks to the model and see how collapse can easily occur. The MSY framework is due to {cite:t}`schaefer1954` and {cite:t}`gordon1954`; a classic textbook treatment is {cite:t}`clark1990`. +This lecture was partly inspired by [a discussion of MSY](https://maths-from-the-past.org/2025/05/12/oversimplified-models-fishery-collapses-and-msy/) posted on [Maths from the past](https://maths-from-the-past.org). + Let's start with some standard imports. ```{code-cell} ipython3 @@ -63,9 +62,6 @@ q = 0.01 # catchability coefficient ## The model -Let $x_t$ be the **stock biomass** --- the total weight of fish, say in tonnes --- -at the start of year $t$. - We build the model in two steps: first how the stock grows on its own, then what fishing does to it. @@ -78,7 +74,11 @@ $$ G(x) = r\,x\left(1 - \frac{x}{K}\right) $$ (eq:logistic) -Here $r > 0$ is called the **intrinsic growth rate** and $K > 0$ is called the **carrying capacity**. +Here + +* $x$ is the **stock biomass** --- the total weight of fish in tonnes, +* $r > 0$ is called the **intrinsic growth rate**, and +* $K > 0$ is called the **carrying capacity**. Let's encode the growth function in Python. @@ -117,7 +117,7 @@ $$ x_{t+1} = x_t + G(x_t). $$ -A clean way to see where the dynamics lead is a **45-degree diagram**: we plot +One way to see where the dynamics lead is via a 45-degree diagram, which plots next year's stock $x_{t+1}$ against this year's stock $x_t$. Wherever the curve crosses the $45^\circ$ line we have $x_{t+1} = x_t$ --- a @@ -177,8 +177,8 @@ plt.tight_layout() plt.show() ``` -Starting from any positive stock, the staircase climbs to the carrying capacity -$x = K$. +Starting from any positive stock and following the staircase, we see that $x_t$ +climbs to the steady state carrying capacity $x = K$. The unfished population fills up the environment and stays there. @@ -194,10 +194,10 @@ Following {cite:t}`schaefer1954`, the catch is proportional to fishing **effort* $e$ (e.g. boat-days) and to the stock available to be caught: $$ - h_t = q\,e\,x_t, + h_t = q\,e\,x_t $$ (eq:harvest) -where $q > 0$ is the **catchability coefficient**. +Here $q > 0$ is the **catchability coefficient**. Subtracting the catch, next year's stock becomes @@ -240,8 +240,8 @@ plt.tight_layout() plt.show() ``` -On the 45-degree diagram, fishing pulls the update curve **down** by the harvest -term, so it now crosses the $45^\circ$ line at a **lower** steady state. +On the 45-degree diagram, fishing pulls the update curve down by the harvest +term, so it now crosses the $45^\circ$ line at a lower steady state. Since this steady state is a function of $e$ now, we denote it by $x^*(e)$. @@ -273,15 +273,15 @@ plt.tight_layout() plt.show() ``` -Not surprisingly, the steady state $x^*(e)$ is decreasing in fishing effort. - +Not surprisingly, the steady state $x^*(e)$ is decreasing in fishing effort +(more fishing leads to less fish). ## Sustainable yield Suppose, as above, that the fleet applies a constant effort $e$ every year. -Let's look for the associated steady state. +Let's look for the associated steady state and associated catch. At any steady state, we must have $x_{t+1} = x_t$. @@ -299,7 +299,7 @@ $$ (eq:xstar) Provided $e < r/q$, we see that each effort level $e$ pins down one steady-state stock $x^*(e)$. -The catch this steady state delivers, year after year, is the **sustainable yield**, defined as +The catch this steady state delivers is called the **sustainable yield**, and defined as $$ y^*(e) := q\,e\,x^*(e). @@ -388,37 +388,6 @@ As effort rises, the steady-state stock falls, while the catch rises and then falls --- exactly the trade-off that the maximum sustainable yield captures. - -## The maximum sustainable yield - -The **maximum sustainable yield** is defined by - -$$ - \text{MSY} \;=\; \max_{0 \le e < r/q} \; y^*(e). -$$ - -This is the largest steady state catch attainable, assuming a constant effort -rate $e$. - -We can compute it numerically by evaluating $y^*(e)$ on a fine grid of effort -levels and picking the largest. - -```{code-cell} ipython3 -e_search = np.linspace(0, r / q, 100_001) -i = np.argmax(sustainable_yield(e_search)) - -e_msy = e_search[i] -MSY = sustainable_yield(e_msy) -x_msy = x_star(e_msy) - -print(f"effort at MSY e_MSY = {e_msy:.2f}") -print(f"stock at MSY x* = {x_msy:.1f} tonnes (= K/2)") -print(f"maximum yield MSY = {MSY:.1f} tonnes/year (= rK/4)") -``` - -The search lands on the round numbers $x^* = K/2$ and $\text{MSY} = rK/4$; the -optional calculus section below shows why. - ### The yield-effort curve To visualize the MSY, we plot the sustainable catch $y^*(e)$ against effort $e$. @@ -437,37 +406,62 @@ y_grid = sustainable_yield(e_grid) fig, ax = plt.subplots() ax.plot(e_grid, y_grid, lw=2, color='C0', label=r'$y^*(e)=qeK\,(1-qe/r)$') -ax.plot([e_msy], [MSY], 'o', color='black', ms=9, zorder=5) -ax.vlines(e_msy, 0, MSY, ls='--', color='black', lw=1) -ax.hlines(MSY, 0, e_msy, ls='--', color='black', lw=1) -ax.annotate(f'MSY = {MSY:.0f}', xy=(e_msy, MSY), xytext=(e_msy + 1, MSY - 18), - fontsize=11, color='black') -ax.annotate(r'$e_{MSY}$', xy=(e_msy, 0), xytext=(e_msy + 1, 6), - fontsize=11, color='black') - ax.set_xlabel('fishing effort $e$') ax.set_ylabel('sustainable yield $y^*(e)$') ax.set_xlim(0, r / q) -ax.set_ylim(0, MSY * 1.25) +ax.set_ylim(0, y_grid.max() * 1.25) ax.legend(loc='upper right', frameon=False) plt.tight_layout() plt.show() ``` -The sustainable yield rises with effort up to the peak at $e_{MSY}$ and then -falls. +The sustainable yield rises with effort and then falls. + +The maximum height of this curve is called the **maximum sustainable yield** (MSY). -Pushing effort beyond $e_{MSY}$ is doubly costly: the catch is smaller *and* the -stock left in the ocean is smaller. +It is the largest steady state catch attainable, assuming a constant effort rate $e$. -It is tempting to call any effort below $e_{MSY}$ "safe", but we resist that -label: once randomness enters the picture, as it does in the next section, no -constant policy is entirely safe. +### Calculating the maximum sustainable yield -### Does the stock settle there? +The maximum sustainable yield can be defined more mathematically as -Finally, let's run the yearly recursion {eq}`eq:update` forward from several +$$ + \text{MSY} \;=\; \max_{0 \le e < r/q} \; y^*(e). +$$ + +The effort level that solves this maximization problem is represented by $e_{MSY}$. + +We can compute it either by calculus or numerically, by evaluating $y^*(e)$ on a fine grid of effort +levels and picking the largest. + +Below we look at the solution using calculus. + +For now we will use the numerical approach: + +```{code-cell} ipython3 +e_search = np.linspace(0, r / q, 100_001) +i = np.argmax(sustainable_yield(e_search)) + +e_msy = e_search[i] +MSY = sustainable_yield(e_msy) +x_msy = x_star(e_msy) + +print(f"effort at MSY e_MSY = {e_msy:.2f}") +print(f"stock at MSY x* = {x_msy:.1f} tonnes (= K/2)") +print(f"maximum yield MSY = {MSY:.1f} tonnes/year (= rK/4)") +``` + +The search lands on the round numbers $x^* = K/2$ and $\text{MSY} = rK/4$; the +optional calculus section below shows why. + + +### Dynamics + +What happens to biomass when we set $e = e_{MSY}$ from some arbitrary initial +$x_0$? + +Below we run the yearly recursion {eq}`eq:update` forward from several starting stocks, under the MSY effort. ```{code-cell} ipython3 @@ -500,9 +494,9 @@ plt.tight_layout() plt.show() ``` -Under the MSY effort, every path climbs or falls toward $x^* = K/2$ and stays. +We see that, under the MSY effort, every path climbs or falls toward $x^* = K/2$ and stays. -So the catch settles down to the MSY. +As a result, the catch settles down to the MSY. ```{note} The smooth, steady convergence above relies on $r$ being modest. @@ -564,6 +558,24 @@ more of the control, less of the flow. ## Randomness, risk and collapse +As mentioned in the introduction to this lecture, MSY-based policies have not +always been successful in sustaining fish populations. + +A major reason is that the simple MSY model ignores randomness and risk. + +This randomness and risk is associated with several features of fisheries. + +One feature is that fish stocks are difficult to track, so policies may be based +on incorrect measurements. + +Another is that ocean environments are complex and nonstationary: a policy that looks safe in a deterministic model +can be dangerously fragile once we admit randomness. + +In this section we investigate how one form of randomness --- stochastic growth +--- can impact fishery dynamics under MSY-based management. + +### Adding randomness + So far the ocean has been a perfectly predictable place. Each year the stock grows by exactly $G(x_t)$. From ab22a82887c1884982eafcc5c376048c45e8b7a4 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 16:38:47 +1000 Subject: [PATCH 8/9] Add MSY success story: Pacific Coast lingcod recovery Insert a new section before "Randomness, risk and collapse" so the lecture shows MSY-based management working before discussing its failures. Uses RAM Legacy Stock Assessment Database series for U.S. Pacific Coast lingcod (B/B_MSY and F/F_MSY, 1928-2009): sustained overfishing drives biomass to ~0.27 B_MSY, then an MSY-anchored rebuilding plan (overfished 1999, rebuilt 2005) pulls F below F_MSY and biomass recovers and overshoots -- the model's negative feedback in real data. Includes honest caveats that tee up the randomness section. Vendor the small series as datasets/lingcod_msy_recovery.csv (read locally, no build-time download); add pandas import and the ricard2012 bib entry. Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/_static/quant-econ.bib | 11 ++ lectures/datasets/lingcod_msy_recovery.csv | 83 +++++++++++++++ lectures/msy_fishery.md | 113 +++++++++++++++++++++ 3 files changed, 207 insertions(+) create mode 100644 lectures/datasets/lingcod_msy_recovery.csv diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index b20cef55b..87f9d61d0 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -2957,3 +2957,14 @@ @book{clark1990 year={1990}, publisher={Wiley} } + +@article{ricard2012, + title={Evaluating the knowledge base and status of commercially exploited marine species with the {RAM} Legacy Stock Assessment Database}, + author={Ricard, Daniel and Minto, Coilin and Jensen, Olaf P. and Baum, Julia K.}, + journal={Fish and Fisheries}, + volume={13}, + number={4}, + pages={380--398}, + year={2012}, + publisher={Wiley} +} diff --git a/lectures/datasets/lingcod_msy_recovery.csv b/lectures/datasets/lingcod_msy_recovery.csv new file mode 100644 index 000000000..ca84787c7 --- /dev/null +++ b/lectures/datasets/lingcod_msy_recovery.csv @@ -0,0 +1,83 @@ +year,B_over_Bmsy,F_over_Fmsy +1928,3.859289617,0.12815534 +1929,3.842213115,0.177669903 +1930,3.825136612,0.198058252 +1931,3.790983607,0.19223301 +1932,3.756830601,0.142718447 +1933,3.756830601,0.224271845 +1934,3.705601093,0.14368932 +1935,3.705601093,0.170873786 +1936,3.671448087,0.133009709 +1937,3.671448087,0.173786408 +1938,3.654371585,0.146601942 +1939,3.637295082,0.127184466 +1940,3.637295082,0.145631068 +1941,3.637295082,0.117475728 +1942,3.637295082,0.067961165 +1943,3.654371585,0.131067961 +1944,3.654371585,0.132038835 +1945,3.654371585,0.130097087 +1946,3.637295082,0.216504854 +1947,3.603142077,0.430097087 +1948,3.50068306,0.462135922 +1949,3.398224044,0.413592233 +1950,3.31284153,0.45631068 +1951,3.227459016,0.453398058 +1952,3.142076503,0.350485437 +1953,3.090846995,0.246601942 +1954,3.073770492,0.290291262 +1955,3.056693989,0.308737864 +1956,3.039617486,0.355339806 +1957,3.005464481,0.526213592 +1958,2.920081967,0.539805825 +1959,2.851775956,0.466990291 +1960,2.800546448,0.417475728 +1961,2.74931694,0.449514563 +1962,2.698087432,0.386407767 +1963,2.663934426,0.391262136 +1964,2.629781421,0.340776699 +1965,2.612704918,0.413592233 +1966,2.578551913,0.518446602 +1967,2.527322404,0.57961165 +1968,2.476092896,0.612621359 +1969,2.407786885,0.575728155 +1970,2.339480874,0.850485437 +1971,2.237021858,1.077669903 +1972,2.117486339,1.563106796 +1973,1.929644809,1.766990291 +1974,1.741803279,1.951456311 +1975,1.572745902,1.951456311 +1976,1.449795082,2.077669903 +1977,1.321721311,1.485436893 +1978,1.270491803,1.941747573 +1979,1.166325137,2.718446602 +1980,1.00068306,3.378640777 +1981,0.840163934,3.378640777 +1982,0.712090164,3.504854369 +1983,0.594262295,2.67961165 +1984,0.544740437,2.805825243 +1985,0.488387978,3.951456311 +1986,0.403005464,3.854368932 +1987,0.37568306,4.019417476 +1988,0.37397541,4.310679612 +1989,0.358606557,4.834951456 +1990,0.3125,4.553398058 +1991,0.285177596,4.32038835 +1992,0.268101093,4.174757282 +1993,0.256147541,3.485436893 +1994,0.266393443,2.27184466 +1995,0.315915301,2.019417476 +1996,0.365437158,2.038834951 +1997,0.396174863,1.815533981 +1998,0.425204918,1.300970874 +1999,0.479849727,1.398058252 +2000,0.539617486,0.63592233 +2001,0.650614754,0.49223301 +2002,0.800887978,1.029126214 +2003,0.988729508,1.291262136 +2004,1.243169399,0.195145631 +2005,1.656420765,0.261165049 +2006,2.083333333,0.249514563 +2007,2.510245902,0.144660194 +2008,2.885928962,0.093203883 +2009,3.193306011, diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index 555f09071..a8c207f5f 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -49,6 +49,7 @@ Let's start with some standard imports. ```{code-cell} ipython3 import numpy as np +import pandas as pd import matplotlib.pyplot as plt ``` @@ -556,6 +557,118 @@ Over-fishing is then the exact analogue of over-saving (dynamic inefficiency): more of the control, less of the flow. ``` +## A success story: Pacific Coast lingcod + +Before turning to the risks, it is worth seeing the MSY framework succeed. + +A clean example is the lingcod (*Ophiodon elongatus*) fishery off the U.S. +Pacific Coast. + +Lingcod is managed using MSY-based reference points: a target biomass and a +target fishing pressure that together define the maximum sustainable yield. + +To follow the fishery we use two dimensionless ratios. + +The first is $B / B_{MSY}$, the stock biomass relative to the biomass $B_{MSY}$ +that supports the MSY. + +(In our model this reference stock is $x^* = K/2$.) + +The second is $F / F_{MSY}$, fishing pressure relative to the pressure +$F_{MSY}$ that achieves the MSY. + +(In our model this is the MSY effort $e_{MSY}$.) + +Each ratio is measured against its own MSY reference point, so the value $1$ is +simultaneously the target and the limit. + +The data come from the RAM Legacy Stock Assessment Database {cite:t}`ricard2012`. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Pacific Coast lingcod recovery + name: fig:lingcod +--- +lingcod = pd.read_csv('datasets/lingcod_msy_recovery.csv') + +fig, ax = plt.subplots() + +ax.axhline(1.0, color='black', lw=1, ls='--') +ax.text(1930, 1.05, r'MSY reference ($=1$)', fontsize=9, va='bottom') + +ymax = max(lingcod['B_over_Bmsy'].max(), lingcod['F_over_Fmsy'].max()) * 1.06 + +# shade the years when fishing pressure exceeded the MSY level +over = lingcod['F_over_Fmsy'] > 1 +ax.fill_between(lingcod['year'], 0, ymax, where=over, + color='C3', alpha=0.08) + +ax.plot(lingcod['year'], lingcod['B_over_Bmsy'], color='C0', lw=2, + label=r'$B / B_{MSY}$ (biomass)') +ax.plot(lingcod['year'], lingcod['F_over_Fmsy'], color='C3', lw=2, + label=r'$F / F_{MSY}$ (fishing pressure)') + +ax.axvline(1999, color='grey', lw=0.8, ls=':') +ax.axvline(2005, color='grey', lw=0.8, ls=':') +ax.text(2000, ymax * 0.97, 'overfished (1999)', ha='left', va='top', + fontsize=8.5, color='dimgrey') +ax.text(2000, ymax * 0.88, 'rebuilt (2005)', ha='left', va='top', + fontsize=8.5, color='dimgrey') + +ax.set_xlabel('year') +ax.set_ylabel('ratio to MSY reference point') +ax.set_xlim(lingcod['year'].min(), lingcod['year'].max()) +ax.set_ylim(0, ymax) +ax.legend(loc='upper center', frameon=False) +plt.tight_layout() +plt.show() +``` + +The story unfolds in three acts. + +Through the 1960s the stock is lightly exploited: $F$ sits well below $F_{MSY}$ +and biomass drifts down slowly from a high level. + +From the early 1970s fishing pressure climbs above $F_{MSY}$ and stays there for +nearly three decades. + +Biomass falls steadily, reaching a trough of about $0.27\,B_{MSY}$ in 1993 --- +close to collapse. + +In 1999 the stock was formally declared *overfished*, and a rebuilding plan cut +catches sharply. + +Fishing pressure dropped well below $F_{MSY}$, and biomass climbed back through +$B_{MSY}$ by 2004. + +The stock was declared fully *rebuilt* in 2005, ahead of schedule, and then +overshot its target. + +This is exactly the negative feedback built into the model. + +Push fishing above the MSY level and the stock falls; pull it below and the +stock recovers toward $B_{MSY}$. + +Used as a steering target, the MSY reference point did its job. + +A few caveats keep the story honest --- and set up the next section. + +First, attribution is never clean: a strong 1999 year-class also helped, so the +cut in fishing was necessary but recruitment luck set the speed of recovery. + +Second, $B_{MSY}$ and $F_{MSY}$ are *estimated* quantities, revised at each +assessment; the famous MSY failures are often failures of estimating the +reference point rather than of the control rule itself. + +Third, the recovery rode partly on coast-wide measures aimed at the whole +groundfish community, a reminder that single-species MSY targets sit inside a +multispecies system. + +The next section takes these caveats seriously, and asks what randomness does to +MSY-based management. + ## Randomness, risk and collapse As mentioned in the introduction to this lecture, MSY-based policies have not From da50ea185dee420633db8e9d5449f09c86c32ec3 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 2 Jun 2026 17:41:45 +1000 Subject: [PATCH 9/9] Use default color cycle for data series; minor prose/tense fixes Per QuantEcon convention, let matplotlib pick data-series colors (C0/C1 standard pair) instead of hard-coding C0/C3; reserve explicit color for reference elements only (black 45-degree line/cobweb/markers/guides, grey milestone lines). Markers now follow their line via get_color(). Lingcod overfishing band switched from red to neutral grey. Also fold in the user's prose edits (overview rewording, G(x) definition, hide-input on plot_45, simplified update equation, tightened lingcod narrative) and fix a past/present tense slip in that narrative. Co-Authored-By: Claude Opus 4.8 (1M context) --- lectures/msy_fishery.md | 131 +++++++++++++++++++--------------------- 1 file changed, 61 insertions(+), 70 deletions(-) diff --git a/lectures/msy_fishery.md b/lectures/msy_fishery.md index a8c207f5f..a08194498 100644 --- a/lectures/msy_fishery.md +++ b/lectures/msy_fishery.md @@ -24,13 +24,15 @@ imposed quotas and restrictions based on MSY. In some cases these management efforts were successful. -However, this hasn't always been the case. +In other cases outcomes were disappointing. -For example, several major fisheries collapsed under MSY-based management, +In fact, several major fisheries collapsed under MSY-based management, including the Peruvian anchovy in the 1970s and the Atlantic cod off Newfoundland in 1992. -In this lecture our main task is to explain how MSY is computed. +In this lecture we provide an introduction to MSY-based management of fisheries. + +Our main task is to explain how MSY is computed. We begin with a relatively elementary treatment in discrete time. @@ -77,7 +79,8 @@ $$ (eq:logistic) Here -* $x$ is the **stock biomass** --- the total weight of fish in tonnes, +* $x$ is the **current biomass** --- the total weight of fish in tonnes, +* $G(x)$ is **annual growth** in biomass --- the difference between current and next year's biomass, * $r > 0$ is called the **intrinsic growth rate**, and * $K > 0$ is called the **carrying capacity**. @@ -103,7 +106,7 @@ mystnb: x_grid = np.linspace(0, K, 400) fig, ax = plt.subplots() -ax.plot(x_grid, G(x_grid), lw=2, color='C0') +ax.plot(x_grid, G(x_grid), lw=2) ax.set_xlabel('stock biomass $x$') ax.set_ylabel('annual growth $G(x)$') ax.set_xlim(0, K) @@ -121,30 +124,34 @@ $$ One way to see where the dynamics lead is via a 45-degree diagram, which plots next year's stock $x_{t+1}$ against this year's stock $x_t$. -Wherever the curve crosses the $45^\circ$ line we have $x_{t+1} = x_t$ --- a -**steady state**, a stock that exactly reproduces itself. +Wherever the curve crosses the $45^\circ$ line we have $x_{t+1} = x_t$. + +The corresponding value of $x$ obeys $x = x + G(x)$. + +Such an $x$ is called a **steady state**: a stock that exactly reproduces itself. -We can then trace the dynamics by "staircasing": from a starting stock go *up* to +We can trace the dynamics of the model by "staircasing": from a starting stock go *up* to the curve (that gives next year's stock), *across* to the $45^\circ$ line (that becomes this year's stock), and repeat. -The next function draws such a diagram. - -It takes the one-year update rule as a function argument `update_fn`, since at -this point we have not yet introduced fishing. +We start with some plotting code. ```{code-cell} ipython3 -def plot_45(ax, update_fn, x0, x_max, steady_state, ss_label, map_label, n_years=30): - "Draw a 45-degree (cobweb) diagram for a one-year stock update rule." +:tags: [hide-input] + +def plot_45( + ax, f, x0, x_max, steady_state, ss_label, map_label, n_years=30 + ): + "Draw a 45-degree (cobweb) diagram for a function f." grid = np.linspace(0, x_max, 400) - ax.plot(grid, update_fn(grid), color='C0', lw=2, label=map_label) + ax.plot(grid, f(grid), lw=2, label=map_label) ax.plot(grid, grid, color='0.6', lw=1, ls='--', label=r'$45^\circ$ line') # cobweb staircase starting from x0 x = x0 cx, cy = [x], [0.0] for _ in range(n_years): - y = update_fn(x) + y = f(x) cx += [x, y] cy += [y, y] x = y @@ -189,7 +196,7 @@ stock grows away from it.) ### Adding fishing -Now let a fishing fleet remove a catch each year. +Now let a fishing fleet remove some catch quantity $h_t$ each year. Following {cite:t}`schaefer1954`, the catch is proportional to fishing **effort** $e$ (e.g. boat-days) and to the stock available to be caught: @@ -203,10 +210,7 @@ Here $q > 0$ is the **catchability coefficient**. Subtracting the catch, next year's stock becomes $$ - x_{t+1} - \;=\; x_t - \;+\; \underbrace{r\,x_t\!\left(1-\frac{x_t}{K}\right)}_{\text{growth}} - \;-\; \underbrace{qe\,x_t}_{\text{catch}}. + x_{t+1} = x_t + G(x_t) - qe\,x_t. $$ (eq:update) @@ -259,10 +263,10 @@ grid = np.linspace(0, 1100, 400) fig, ax = plt.subplots(figsize=(4.95, 4.95)) ax.plot(grid, grid, color='0.6', lw=1, ls='--', label=r'$45^\circ$ line') -for e, c in zip((10.0, 30.0), ('C0', 'C3')): - ax.plot(grid, update(grid, e), lw=2, color=c, label=f'$e={e:.0f}$') +for e in (10.0, 30.0): + line, = ax.plot(grid, update(grid, e), lw=2, label=f'$e={e:.0f}$') xs = x_star(e) - ax.plot([xs], [xs], 'o', color=c, ms=8, zorder=5) + ax.plot([xs], [xs], 'o', color=line.get_color(), ms=8, zorder=5) ax.set_xlabel('stock this year $x_t$') ax.set_ylabel('stock next year $x_{t+1}$') @@ -330,8 +334,8 @@ mystnb: x = np.linspace(0, K, 400) fig, ax = plt.subplots() -ax.plot(x, G(x), lw=2, color='C0', label=r'growth $G(x)$') -ax.plot(x, q * e_demo * x, lw=2, color='C3', label=r'harvest $q e x$') +ax.plot(x, G(x), lw=2, label=r'growth $G(x)$') +ax.plot(x, q * e_demo * x, lw=2, label=r'harvest $q e x$') xs = x_star(e_demo) ys = sustainable_yield(e_demo) @@ -363,15 +367,14 @@ mystnb: name: fig:steady-states --- fig, ax = plt.subplots() -ax.plot(x, G(x), lw=2, color='C0', label='growth $G(x)$') +ax.plot(x, G(x), lw=2, label='growth $G(x)$') efforts = [12.5, 25.0, 37.5] labels = [r'low $e$', r'moderate $e$', r'high $e$'] -colors = ['C2', 'C3', 'C1'] -for e, lab, c in zip(efforts, labels, colors): - ax.plot(x, q * e * x, lw=2, color=c, label=lab) - ax.plot([x_star(e)], [q * e * x_star(e)], 'o', color=c, ms=7, zorder=5) +for e, lab in zip(efforts, labels): + line, = ax.plot(x, q * e * x, lw=2, label=lab) + ax.plot([x_star(e)], [q * e * x_star(e)], 'o', color=line.get_color(), ms=7, zorder=5) ax.set_xlabel('stock $x$') ax.set_ylabel('catch') @@ -406,7 +409,7 @@ e_grid = np.linspace(0, r / q, 400) y_grid = sustainable_yield(e_grid) fig, ax = plt.subplots() -ax.plot(e_grid, y_grid, lw=2, color='C0', label=r'$y^*(e)=qeK\,(1-qe/r)$') +ax.plot(e_grid, y_grid, lw=2, label=r'$y^*(e)=qeK\,(1-qe/r)$') ax.set_xlabel('fishing effort $e$') ax.set_ylabel('sustainable yield $y^*(e)$') ax.set_xlim(0, r / q) @@ -564,10 +567,12 @@ Before turning to the risks, it is worth seeing the MSY framework succeed. A clean example is the lingcod (*Ophiodon elongatus*) fishery off the U.S. Pacific Coast. -Lingcod is managed using MSY-based reference points: a target biomass and a -target fishing pressure that together define the maximum sustainable yield. +Lingcod is managed using MSY-based analysis. -To follow the fishery we use two dimensionless ratios. +This analysis leads to a target biomass and a target fishing pressure that +together define the maximum sustainable yield. + +To follow the fishery we use two ratios. The first is $B / B_{MSY}$, the stock biomass relative to the biomass $B_{MSY}$ that supports the MSY. @@ -579,9 +584,6 @@ $F_{MSY}$ that achieves the MSY. (In our model this is the MSY effort $e_{MSY}$.) -Each ratio is measured against its own MSY reference point, so the value $1$ is -simultaneously the target and the limit. - The data come from the RAM Legacy Stock Assessment Database {cite:t}`ricard2012`. ```{code-cell} ipython3 @@ -603,71 +605,60 @@ ymax = max(lingcod['B_over_Bmsy'].max(), lingcod['F_over_Fmsy'].max()) * 1.06 # shade the years when fishing pressure exceeded the MSY level over = lingcod['F_over_Fmsy'] > 1 ax.fill_between(lingcod['year'], 0, ymax, where=over, - color='C3', alpha=0.08) + color='0.5', alpha=0.12) -ax.plot(lingcod['year'], lingcod['B_over_Bmsy'], color='C0', lw=2, +ax.plot(lingcod['year'], lingcod['B_over_Bmsy'], lw=2, label=r'$B / B_{MSY}$ (biomass)') -ax.plot(lingcod['year'], lingcod['F_over_Fmsy'], color='C3', lw=2, +ax.plot(lingcod['year'], lingcod['F_over_Fmsy'], lw=2, label=r'$F / F_{MSY}$ (fishing pressure)') ax.axvline(1999, color='grey', lw=0.8, ls=':') ax.axvline(2005, color='grey', lw=0.8, ls=':') -ax.text(2000, ymax * 0.97, 'overfished (1999)', ha='left', va='top', - fontsize=8.5, color='dimgrey') -ax.text(2000, ymax * 0.88, 'rebuilt (2005)', ha='left', va='top', - fontsize=8.5, color='dimgrey') ax.set_xlabel('year') ax.set_ylabel('ratio to MSY reference point') ax.set_xlim(lingcod['year'].min(), lingcod['year'].max()) ax.set_ylim(0, ymax) -ax.legend(loc='upper center', frameon=False) +ax.legend(loc='upper left', frameon=False) plt.tight_layout() plt.show() ``` -The story unfolds in three acts. -Through the 1960s the stock is lightly exploited: $F$ sits well below $F_{MSY}$ -and biomass drifts down slowly from a high level. +Through the 1960s the stock was lightly exploited: $F$ sat well below $F_{MSY}$ +and biomass drifted down slowly from a high level. -From the early 1970s fishing pressure climbs above $F_{MSY}$ and stays there for +From the early 1970s fishing pressure climbed above $F_{MSY}$ and stayed there for nearly three decades. -Biomass falls steadily, reaching a trough of about $0.27\,B_{MSY}$ in 1993 --- +Biomass fell steadily, reaching a trough of about $0.27\,B_{MSY}$ in 1993 --- close to collapse. -In 1999 the stock was formally declared *overfished*, and a rebuilding plan cut -catches sharply. +In 1999 the stock was formally declared *overfished*, and a rebuilding plan +based around MSY cut catches sharply. Fishing pressure dropped well below $F_{MSY}$, and biomass climbed back through $B_{MSY}$ by 2004. -The stock was declared fully *rebuilt* in 2005, ahead of schedule, and then +The stock was declared fully rebuilt in 2005, ahead of schedule, and then overshot its target. -This is exactly the negative feedback built into the model. - -Push fishing above the MSY level and the stock falls; pull it below and the -stock recovers toward $B_{MSY}$. +As suggested by the model, fishing above the MSY level has led to falling stock, while +keeping it below has led to recovery. -Used as a steering target, the MSY reference point did its job. -A few caveats keep the story honest --- and set up the next section. +Of course, the real world is not as clean as the model and random factors +also influenced outcomes. -First, attribution is never clean: a strong 1999 year-class also helped, so the -cut in fishing was necessary but recruitment luck set the speed of recovery. +For example, a strong 1999 year-class also helped, so the +cut in fishing was important but recruitment luck set the speed of recovery. -Second, $B_{MSY}$ and $F_{MSY}$ are *estimated* quantities, revised at each -assessment; the famous MSY failures are often failures of estimating the -reference point rather than of the control rule itself. +In addition, the recovery was aided by coast-wide measures aimed at the whole +groundfish community. -Third, the recovery rode partly on coast-wide measures aimed at the whole -groundfish community, a reminder that single-species MSY targets sit inside a -multispecies system. +The next section expands our model somewhat, with the aim of introducing more +realistic dynamics. -The next section takes these caveats seriously, and asks what randomness does to -MSY-based management. ## Randomness, risk and collapse