(Embedded for self-containment. John may also drag-drop the converted markdown onto this issue.)
\documentclass[12pt]{article}
\usepackage{amsmath, amssymb, amsthm}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{booktabs}
\usepackage{enumitem}
\usepackage{xcolor}
\usepackage{natbib}
\geometry{margin=1.2in}
\hypersetup{
colorlinks=true,
linkcolor=blue!60!black,
citecolor=blue!60!black,
urlcolor=blue!60!black
}
\newtheorem{remark}{Remark}
\newtheorem{example}{Example}
\newtheorem{definition}{Definition}
\title{\textbf{Bayesian Surplus Production Models for Fisheries Stock
Assessment}\\[6pt]
\large{Notes Toward a QuantEcon Lecture}\\[4pt]
\normalsize{Third in the sequence: Unemployment Dynamics $\to$
Population State Space $\to$ Fisheries}}
\author{John Stachurski}
\date{\today}
\begin{document}
\maketitle
\begin{abstract}
\noindent
These notes describe a QuantEcon lecture on Bayesian estimation of a
state space surplus production model (SSPM) for fisheries stock
assessment. The lecture is the third in a sequence: the first
(unemployment dynamics) introduced Bayesian estimation of nonlinear
time series with directly observed states; the second (animal
population dynamics) introduced latent state inference via a logistic
growth state space model. This lecture adds two elements that
complete the transition to a research-grade fisheries model: a catch
term in the biomass dynamics, and a catchability-linked observation
equation connecting CPUE survey indices to latent biomass. The
theoretical foundation is Meyer and Millar (1999) and Punt et al.\
(2003). The implementation uses NumPyro/JAX. The application
dataset is Pacific Coast lingcod from the RAM Legacy Stock Assessment
Database. The lecture ends by connecting the posterior output
directly to a POMDP/CVaR policy layer, which is the entry point to
the GRIPS fisheries management research programme.
\end{abstract}
\tableofcontents
\bigskip
%=============================================================
\section{Context and Motivation}
%=============================================================
\subsection{Position in the lecture sequence}
This is the third lecture in a Bayesian methods sequence designed to
build, step by step, toward a research-grade Bayesian fisheries stock
assessment tool in JAX. The sequence is:
\begin{enumerate}[itemsep=4pt]
\item \textbf{Nonlinear unemployment dynamics.} Bayesian estimation
of a $\sinh$-type nonlinear AR model for the US unemployment
rate. States are directly observed. Teaches: nonlinear
likelihoods, prior design, NUTS in NumPyro, posterior
interpretation, comparison with Metropolis--Hastings.
\item \textbf{Animal population dynamics.} Bayesian estimation of
a logistic growth state space model for the Yellowstone wolf
population. Introduces latent states, the Kalman filter as a
linear Gaussian benchmark, Neal's funnel and non-centred
parameterisation, the signal-to-noise identification problem,
and posterior predictive checks. The model is a
single-species surplus production model without harvesting.
\item \textbf{Fisheries surplus production model (this lecture).}
Adds a catch term to the population dynamics and a
catchability-linked observation equation. Applies the full
model to a real stock assessment dataset. Connects the
posterior output to a management policy layer.
\end{enumerate}
A student who has worked through lectures one and two is
computationally and conceptually prepared for this lecture. The
NumPyro model for the SSPM differs from the animal population model
by approximately ten lines of code. The new intellectual content is:
(a) the economics and ecology of harvesting and surplus production;
(b) the additional identification problems introduced by catchability
and multiple CPUE indices; and (c) the connection between the
posterior and downstream decision-making under uncertainty.
\subsection{Why a JAX-native implementation?}
The standard tools for Bayesian surplus production models are JAGS
(used by Meyer and Millar 1999) and the R package JABBA (Winker et
al.\ 2018), which wraps JAGS with fisheries-specific scaffolding.
These tools are adequate for producing point estimates and confidence
intervals, but they have three limitations that matter for the
research programme motivating this lecture.
First, JAGS and JABBA do not run on GPU and are not differentiable
end-to-end. A NumPyro/JAX implementation is differentiable through
the posterior computation, which opens the door to gradient-based
policy optimisation.
Second, the posterior samples from a JAX implementation are JAX
arrays. This means the posterior distribution over current biomass
$B_T$ — the belief state in a POMDP formulation — can be passed
directly into a JAX-based value function computation without any data
serialisation. In JABBA, the analogous operation requires reading
JAGS output from disk and reimporting it into a separate optimisation
environment.
Third, a clean NumPyro implementation is pedagogically transparent.
JABBA has approximately 2000 lines of R code covering diagnostics,
plotting, and multiple model variants; the Meyer-Millar model itself
is perhaps 30 lines of BUGS code embedded within it. A standalone
NumPyro implementation is 100--150 lines and exposes every modelling
assumption explicitly.
\subsection{The gap in the literature}
As of mid-2025, there is no published, general-purpose, JAX-native
Bayesian surplus production model. The closest things are:
\begin{itemize}[itemsep=4pt]
\item JABBA (Winker et al.\ 2018): JAGS/R, not differentiable
\item \texttt{MSEtool} (Huynh et al.\ 2020): R/TMB, frequentist
\item Various Stan implementations in the grey literature: not
GPU-native, not JAX
\item FIMS (NOAA, 2023--): C++/TMB, not Python
\end{itemize}
A clean NumPyro implementation, validated against published
assessments and written as a QuantEcon lecture, would be a genuine
contribution to the open-source fisheries stock assessment ecosystem.
%=============================================================
\section{The Schaefer Surplus Production Model}
%=============================================================
\subsection{Biological background}
A surplus production model (SPM) describes the dynamics of total fish
biomass $B_t$ at the stock level, without age or size structure. The
foundational concept is \emph{surplus production}: the net biomass
produced by the stock above what is needed to maintain its current
size. In the absence of fishing, surplus production goes into
biomass growth. Fishing removes biomass; the stock is sustainable
when catch does not exceed surplus production.
The Schaefer (1954) model specifies surplus production as a logistic
function of current biomass:
\begin{equation}
P(B) \;=\; rB\!\left(1 - \frac{B}{K}\right),
\label{eq:surplus}
\end{equation}
where $r > 0$ is the intrinsic growth rate and $K > 0$ is the
carrying capacity (the unfished equilibrium biomass). Surplus
production is zero at $B = 0$ and $B = K$, and is maximised at $B =
K/2$, the biomass that supports maximum sustainable yield (MSY).
The MSY harvest level is:
\begin{equation}
\mathrm{MSY} \;=\; P(K/2) \;=\; \frac{rK}{4},
\label{eq:msy}
\end{equation}
and the corresponding fishing mortality rate is $F_{\mathrm{MSY}} =
r/2$. These reference points appear in the posterior predictive
analysis and in the policy layer.
\subsection{The Pella-Tomlinson generalisation}
The Schaefer model imposes that MSY occurs at exactly $B = K/2$.
Empirically this is often too restrictive. The Pella-Tomlinson (1969)
model generalises equation~\eqref{eq:surplus} to:
\begin{equation}
P(B) \;=\; \frac{r}{m} B\!\left[1 - \left(\frac{B}{K}\right)^{m-1}\right],
\label{eq:pt}
\end{equation}
where $m > 0$ is a shape parameter. The Schaefer model is the
special case $m = 2$. For $m < 2$, the MSY-maximising biomass is
below $K/2$ (left-skewed surplus production); for $m > 2$, it is
above $K/2$.
For the main lecture the Schaefer model ($m = 2$) is the right
starting point. The Pella-Tomlinson generalisation is a natural
extension exercise: add one parameter ($m$, or equivalently the
ratio $B_{\mathrm{MSY}}/K$), assign it a prior, and re-estimate.
JABBA focuses heavily on the Pella-Tomlinson model because many
tuna and billfish stocks have $m \ne 2$; for the lingcod
application the Schaefer model is adequate.
%=============================================================
\section{The State Space Surplus Production Model}
%=============================================================
\subsection{Overview}
The state space surplus production model (SSPM) adds two layers of
realism to the deterministic Schaefer model: stochastic biomass
dynamics (process error) and noisy observation of biomass through
CPUE indices (observation error). This is exactly the structure of
the animal population model from lecture two, with the addition of
the catch term.
\subsection{Notation}
\begin{center}
\begin{tabular}{lll}
\toprule
Symbol & Meaning & Units \\
\midrule
$B_t$ & True (latent) biomass at time $t$ & $10^3$ tonnes \\
$C_t$ & Observed catch at time $t$ & $10^3$ tonnes \\
$I_t$ & CPUE index (survey or commercial) & kg per unit effort \\
$r$ & Intrinsic growth rate & per year \\
$K$ & Carrying capacity & $10^3$ tonnes \\
$q$ & Catchability coefficient & dimensionless \\
$\sigma_p$ & Process noise (log scale) & dimensionless \\
$\sigma_o$ & Observation noise (log scale) & dimensionless \\
$B_0/K$ & Initial depletion (relative to $K$) & dimensionless \\
\bottomrule
\end{tabular}
\end{center}
Catch data $C_{1:T}$ are treated as \textbf{known without error} —
this is the standard assumption in SPMs and is reasonable when
logbook data are available. The latent variable is $B_{1:T}$.
\subsection{Process equation}
\begin{equation}
\log B_t \;=\; \log\!\left[
B_{t-1} + r B_{t-1}\!\left(1 - \frac{B_{t-1}}{K}\right)
- C_{t-1}
\right]
+ \eta_t,
\qquad \eta_t \sim \mathcal{N}(0,\, \sigma_p^2).
\label{eq:process}
\end{equation}
Several points deserve comment.
\paragraph{Log-space noise.} The process noise enters on the log
scale, making the shocks multiplicative on the biomass scale. This
is the Meyer-Millar parameterisation and ensures $B_t > 0$ almost
surely. An alternative is to put the noise on the arithmetic scale
inside the brackets; the log-space version is more tractable and more
standard.
\paragraph{Catch timing.} Equation~\eqref{eq:process} assumes catch
$C_{t-1}$ is removed at the start of period $t-1$ before growth
occurs (a ``pulse'' harvesting assumption). An alternative is to
remove catch at the end of the period (after growth). For annual
models the distinction matters little; for seasonal models it can be
important. The lecture should make the convention explicit.
\paragraph{Constraint $B_t > C_{t-1}$.} The argument of the outer
logarithm must be positive, which requires:
\[
B_{t-1} + rB_{t-1}(1 - B_{t-1}/K) > C_{t-1}.
\]
In practice this is enforced by the prior: if the posterior puts
significant mass on biomass trajectories that dip below the catch
level, it signals a misspecified model or implausible priors. Some
implementations add a small positive floor (e.g.\ $0.001 \cdot K$) to
prevent numerical issues. This should be discussed in the lecture.
\subsection{Observation equation}
\begin{equation}
\log I_t \;=\; \log q + \log B_t + \epsilon_t,
\qquad \epsilon_t \sim \mathcal{N}(0,\, \sigma_o^2).
\label{eq:observation}
\end{equation}
The CPUE index $I_t$ is assumed proportional to biomass, with
proportionality constant $q$ (catchability) and log-normal
multiplicative observation error. The parameter $q$ is not directly
interpretable in isolation; it absorbs units and the relationship
between fishing effort and encounter rate.
\paragraph{Analytical profile-out of $q$.} Because equation
\eqref{eq:observation} is linear in $\log q$, the catchability
parameter can be concentrated out of the likelihood analytically:
\begin{equation}
\hat q \;=\; \exp\!\left(
\frac{1}{T}\sum_{t=1}^T (\log I_t - \log B_t)
\right).
\label{eq:qhat}
\end{equation}
This is the Meyer-Millar trick: replace $q$ as a free parameter with
its analytic conditional MLE, reducing the parameter space by one
dimension. In a Bayesian context this corresponds to using an
improper flat prior on $\log q$ and analytically marginalising.
The lecture should present both approaches: the full Bayesian version
with a proper prior on $\log q$, and the profiled version. The
latter is faster and avoids the $q$--$K$ identification problem
described in Section~\ref{sec:identification}.
\subsection{Initial condition}
The initial biomass $B_1$ is parameterised as a fraction of the
carrying capacity:
\begin{equation}
B_1 \;=\; \phi K, \qquad \phi \in (0,1].
\label{eq:initial}
\end{equation}
Assigning a prior to $\phi$ rather than to $B_1$ directly makes the
prior scale-invariant: the same prior on $\phi$ applies regardless of
whether $K$ is measured in tonnes or kilograms. The standard prior
is $\phi \sim \mathrm{Uniform}(0.2, 1.0)$ or
$\phi \sim \mathrm{Beta}(2, 2)$ restricted to $(0.2, 1.0)$, reflecting
the belief that the stock was not heavily depleted before the data
period began. This prior is often consequential — it should be
treated as a sensitivity analysis parameter.
\subsection{Full model in one place}
Collecting equations \eqref{eq:process} and \eqref{eq:observation}
with the initial condition \eqref{eq:initial}:
\begin{align}
B_1 &= \phi K \label{eq:full1} \\[4pt]
\log B_t &= \log\!\bigl[B_{t-1} + rB_{t-1}(1 - B_{t-1}/K)
- C_{t-1}\bigr] + \eta_t,
\quad \eta_t \sim \mathcal{N}(0, \sigma_p^2),
\quad t = 2,\ldots,T
\label{eq:full2} \\[4pt]
\log I_t &= \log q + \log B_t + \epsilon_t,
\quad \epsilon_t \sim \mathcal{N}(0, \sigma_o^2),
\quad t = 1,\ldots,T.
\label{eq:full3}
\end{align}
Parameters: $\theta = (r, K, q, \sigma_p, \sigma_o, \phi)$.
Latent states: $B_{2:T}$ (given $B_1 = \phi K$).
%=============================================================
\section{Identification Problems}
\label{sec:identification}
%=============================================================
The SSPM has three overlapping identification challenges. Making
these visible through the joint posterior is one of the most
important contributions of the lecture — practitioners using JABBA or
SS3 often do not appreciate how weakly identified the parameters are,
because those tools report only point estimates with asymptotic
standard errors.
\subsection{The $r$--$K$ ridge}
The equilibrium biomass is $B^* \approx K$ (unfished) and $B^* =
K/2$ at MSY. The speed of convergence to these equilibria depends on
$r$. Over a short time series, a stock that is far from equilibrium
makes $r$ and $K$ observationally very similar for different
combinations: high $r$ combined with high $K$ can generate the same
trajectory over 30 years as moderate $r$ combined with moderate $K$.
This appears as a ridge in the joint posterior of $(r, K)$, exactly
analogous to the $(\beta, \lambda)$ ridge in the unemployment lecture
and the $(r, K)$ ridge in the population dynamics lecture.
The ridge is most severe when the stock has not been pushed to
extreme depletion and has not recovered close to $K$ within the data
period. Reference points such as MSY $= rK/4$ may be well-identified
even when $r$ and $K$ individually are not, because MSY lies along
the ridge. The lecture should make this distinction explicit: show
the $(r, K)$ scatter plot, show that MSY $= rK/4$ has a much tighter
posterior, and explain why.
\subsection{The $q$--$K$ problem}
Catchability $q$ and carrying capacity $K$ are identified by the
\emph{level} of the CPUE series relative to the model prediction.
If the CPUE index is high, it could be because the stock is large
(high $K$) or because the gear is very efficient (high $q$). Only
the product $qK$ is identified by the level; separating them requires
either an informative prior on one or absolute abundance data (e.g.\
acoustic surveys that give $\hat B_t$ directly without the $q$
multiplier).
This is the reason for profiling out $q$ analytically (Section 3.4):
doing so removes one dimension of the non-identification. In the
full Bayesian version, an informative prior on $q$ (from tagging
studies or hydroacoustic calibrations) is needed to identify $K$.
The lecture should show the joint posterior of $(q, K)$ in the full
version to illustrate the problem, then show that profiling resolves
it.
\subsection{The signal-to-noise ratio}
The ratio $q_{\sigma} = \sigma_p^2 / \sigma_o^2$ determines how
much of the variation in the CPUE index is attributed to genuine
biomass fluctuations versus observation noise. This is the same
identification problem as in the population dynamics lecture,
compounded here by the $r$--$K$ and $q$--$K$ issues. A large
$\sigma_p$ gives a wiggly latent biomass trajectory that closely
tracks the CPUE; a small $\sigma_p$ gives a smooth trajectory
governed mainly by the deterministic Schaefer dynamics.
The standard finding in the literature (Punt et al.\ 2003, Winker et
al.\ 2018) is that the posterior of $\sigma_p$ is often diffuse
and the data cannot distinguish high process noise from high
observation noise. This is not a failure of the model; it is the
correct probabilistic summary of what the data contain. Managers
should be told about this uncertainty rather than given a false sense
of precision.
%=============================================================
\section{Priors}
%=============================================================
\subsection{Canonical choices}
The following priors are adapted from Meyer and Millar (1999) and
Winker et al.\ (2018), converted to NumPyro parameterisation. All
are stated in the units of the lingcod application (biomass in $10^3$
tonnes, time in years).
\begin{align}
r &\;\sim\; \mathrm{HalfNormal}(0.5)
\tag{P1}\label{prior:r} \\[4pt]
K &\;\sim\; \mathrm{LogNormal}(\mu_K,\, 0.5^2)
\tag{P2}\label{prior:K} \\[4pt]
\log q &\;\sim\; \mathcal{N}(0,\, 1.0^2) \quad
\text{[or profiled out]}
\tag{P3}\label{prior:q} \\[4pt]
\sigma_p &\;\sim\; \mathrm{HalfNormal}(0.2)
\tag{P4}\label{prior:sigmap} \\[4pt]
\sigma_o &\;\sim\; \mathrm{HalfNormal}(0.3)
\tag{P5}\label{prior:sigmao} \\[4pt]
\phi &\;\sim\; \mathrm{Beta}(2,\, 2) \cdot
\mathbf{1}[\phi \in (0.2,\, 1.0)]
\tag{P6}\label{prior:phi}
\end{align}
\paragraph{Setting $\mu_K$.} The prior mean of $\log K$ should be
set to the log of a rough prior estimate of unfished biomass. For
lingcod on the US Pacific Coast, published stock assessments suggest
$K$ in the range 20,000--60,000 tonnes; set $\mu_K = \log(30{,}000)$
with the log-Normal scale of 0.5 giving a 95\% prior interval of
roughly $(11{,}000,\, 81{,}000)$ tonnes. The lecture should
demonstrate prior sensitivity by re-running with $\mu_K = \log(15{,}000)$
and $\mu_K = \log(60{,}000)$.
\paragraph{Prior on $r$.} For groundfish such as lingcod, life-history
theory suggests $r \in (0.1, 0.5)$ per year. The
$\mathrm{HalfNormal}(0.5)$ prior is diffuse over this range. An
informative alternative is $r \sim \mathrm{LogNormal}(\log 0.2, 0.3^2)$,
which encodes the belief that $r \approx 0.2$ with moderate
uncertainty. FishBase life-history parameters (via the
\texttt{rfishbase} R package or direct lookup) provide justification
for these choices.
\paragraph{Prior on $\phi$.} The initial depletion prior is among
the most consequential choices in the model, particularly when the
time series begins after substantial exploitation has already
occurred. A $\mathrm{Beta}(2,2)$ prior on $(0.2, 1.0)$ places most
mass around $\phi = 0.6$, consistent with a moderately exploited
stock. If historical catch records suggest heavy exploitation before
the data period, a prior concentrated at low $\phi$ values (e.g.\
$\phi \sim \mathrm{Beta}(1.5, 4)$ restricted to $(0.05, 0.6)$) is
more appropriate. The lecture should treat this as an explicit
sensitivity parameter.
%=============================================================
\section{Data: Pacific Lingcod}
%=============================================================
\subsection{Why lingcod?}
Pacific Coast lingcod (\textit{Ophiodon elongatus}) is an ideal
teaching dataset for three reasons.
First, it is already used in the QuantEcon MSY lecture (building the
RAM Legacy data pipeline, computing MSY reference points from surplus
production fits). The student has already cleaned the data and
computed point estimates. This lecture re-uses that infrastructure
and adds the Bayesian uncertainty layer.
Second, lingcod has a well-documented stock history: heavy
exploitation in the 1970s--80s brought the stock to low levels,
followed by strong recovery after management reform in the early
2000s. The latent biomass trajectory therefore has clear shape —
it is not a flat or slowly trending series — making the posterior
over $B_{1:T}$ interesting and interpretable.
Third, the Pacific Fishery Management Council (PFMC) publishes
regular stock assessments using SS3, including biomass estimates with
uncertainty. These provide an external benchmark: the Bayesian SSPM
posterior should be broadly consistent with the SS3 assessment,
modulo the different model structure.
\subsection{Data sources}
\begin{center}
\begin{tabular}{llp{7cm}}
\toprule
Series & Source & Notes \\
\midrule
Catch ($C_t$, annual, 1940--2020) &
RAM Legacy Database v4.3 &
Series ID: \texttt{LINGPACcatch}. Units: $10^3$ tonnes.
Already processed in the MSY lecture. \\[6pt]
CPUE index ($I_t$, annual, 1980--2020) &
RAM Legacy Database v4.3 &
Series ID: \texttt{LINGPACcpue}. Several CPUE series are
available; trawl survey biomass index preferred over commercial
CPUE (see below). \\[6pt]
Biomass index (alternative) &
PFMC stock assessment 2021 &
Published biomass estimates with CVs; can be used in place of
raw CPUE to validate the model. \\
\bottomrule
\end{tabular}
\end{center}
\paragraph{RAM Legacy Database access.}
The RAM Legacy Database (v4.3) is available from
\url{https://www.ramlegacy.org}. Download the \texttt{.zip} archive,
which contains \texttt{.RData} files loadable in R or — more
conveniently for this lecture — \texttt{.csv} exports of the main
tables. The relevant columns are \texttt{stockid},
\texttt{year}, \texttt{catch\_landings}, and the CPUE or biomass
index series. The full data pipeline from the MSY lecture should be
reused here; only the output changes (from point estimates to
posterior samples).
\paragraph{Choosing between CPUE series.}
The RAM Legacy Database may include multiple CPUE series for lingcod:
commercial trawl CPUE, commercial hook-and-line CPUE, and trawl
survey indices. For the SSPM, a trawl survey index is preferable to
commercial CPUE because: (a) survey effort is standardised and the
catchability $q$ is more stable over time; (b) commercial CPUE is
confounded with changes in fishing technology and targeting
behaviour. The lecture should discuss this choice and run the model
on both series as a sensitivity analysis.
\subsection{Data preparation}
The following steps should be documented explicitly in the lecture
notebook:
\begin{enumerate}[itemsep=4pt]
\item \textbf{Align catch and index series.} Catch typically runs
from the mid-20th century; the survey index may only be
available from the 1980s or later. The model can be run on
the full catch series with observations missing for the
pre-survey years: the latent state $B_t$ is propagated
forward by the process equation for years without an index
observation, without contributing a likelihood term. This
is handled naturally in NumPyro by looping over only the
observed index years.
\item \textbf{Check for catch exceeding reasonable biomass.} If
$C_t > B_{t-1}$ for any $t$ under prior-plausible biomass
values, the constraint in equation~\eqref{eq:process} will
be violated. Inspect the early catch history; implausibly
high catches may indicate data errors or unreported catch.
\item \textbf{Log-transform the index.} The observation equation
is in log-space; store $\log I_t$ as the working variable.
\item \textbf{Scale.} Work in units of $10^3$ tonnes to keep
the numerical scale of $K$ moderate (order 10--100 rather
than 10,000--100,000). Equivalently, normalise the CPUE
index to have mean 1.0 over the data period and absorb the
scale into $q$.
\end{enumerate}
%=============================================================
\section{Implementation Notes}
%=============================================================
\subsection{Model structure}
The NumPyro model follows the same architecture as the logistic
growth model from lecture two. The main differences are:
\begin{enumerate}[itemsep=4pt]
\item The transition equation subtracts catch $C_{t-1}$ before
taking the logarithm, requiring careful handling of the
positivity constraint.
\item The observation equation includes the catchability parameter
$q$ (or its profiled-out version).
\item Catch $C_{1:T}$ enters as observed data (a fixed array), not
as a stochastic variable.
\item The initial state is $B_1 = \phi K$ rather than a freely
estimated $x_0$.
\end{enumerate}
\subsection{Non-centred parameterisation}
As in the population dynamics lecture, the centred parameterisation
samples $B_t$ directly and suffers from Neal's funnel when $\sigma_p$
is small. The non-centred version samples standardised innovations
$\tilde\eta_t \sim \mathcal{N}(0,1)$ and reconstructs:
\[
\log B_t \;=\; \log\!\bigl[B_{t-1} + rB_{t-1}(1 - B_{t-1}/K)
- C_{t-1}\bigr]
+ \sigma_p \tilde\eta_t.
\]
The non-centred parameterisation should be the default in the lecture
code. The centred version can be shown first for clarity, with the
divergences used to motivate the switch.
\subsection{Numerical stability}
The argument of the $\log$ in the process equation,
\[
A_t \;=\; B_{t-1} + rB_{t-1}(1 - B_{t-1}/K) - C_{t-1},
\]
can become negative during NUTS exploration if the sampler visits
implausible parameter combinations. This causes numerical failures.
Two strategies:
\begin{itemize}[itemsep=4pt]
\item \textbf{Soft floor:} Replace $A_t$ with $\max(A_t,\,
\epsilon K)$ for a small $\epsilon > 0$ (e.g.\ $0.001$).
This is a minor model modification but prevents crashes.
\item \textbf{Informative prior:} Tighten the prior on $r$ and
$K$ to reduce the probability of visiting pathological
parameter combinations. Better in principle but may require
domain knowledge to implement correctly.
\end{itemize}
The lecture should use the soft floor with a comment explaining the
trade-off.
\subsection{Profiling out $q$}
The profiled log-likelihood replaces $\log q$ with
$\hat{\log q} = \frac{1}{T_{\mathrm{obs}}} \sum_{t \in
\mathcal{T}} (\log I_t - \log B_t)$, where $\mathcal{T}$ is the set
of years with observed CPUE. In the NumPyro model this is computed
as a deterministic quantity inside the model function:
\begin{verbatim}
log_q_hat = jnp.mean(log_I_obs - log_B[obs_idx])
resid = log_I_obs - (log_q_hat + log_B[obs_idx])
numpyro.sample("I_obs", dist.Normal(0.0, sigma_o),
obs=resid)
\end{verbatim}
This reduces the parameter count by one and eliminates the $q$--$K$
ridge. It should be presented as a standard technique (Meyer and
Millar 1999, Section 3) rather than a trick.
\subsection{Multiple CPUE indices}
Some stocks have multiple simultaneous CPUE series from different
fleets or survey programmes. The observation equation generalises to:
\begin{equation}
\log I_{it} \;=\; \log q_i + \log B_t + \epsilon_{it},
\qquad \epsilon_{it} \sim \mathcal{N}(0,\, \sigma_{oi}^2),
\label{eq:multi_cpue}
\end{equation}
for index $i = 1,\ldots,n_I$. Each index has its own catchability
$q_i$ and observation noise $\sigma_{oi}$, but all share the same
latent biomass trajectory $B_t$. In the profiled version, each $q_i$
is profiled out separately. This is the main extension in Winker
et al.\ (2018) and is straightforward to implement in NumPyro.
For lingcod, if both commercial CPUE and survey indices are available,
the multi-index version can be fit and compared to the single-index
version via LOO-CV. This is a good exercise in model comparison.
%=============================================================
\section{Posterior Analysis}
%=============================================================
\subsection{Parameter posteriors}
The primary quantities of interest are:
\begin{itemize}[itemsep=4pt]
\item Marginal posteriors of $r$, $K$, $\sigma_p$, $\sigma_o$, $\phi$
\item Joint posterior of $(r, K)$: visualise the ridge and the
implied posterior of MSY $= rK/4$
\item Posterior of MSY and $F_{\mathrm{MSY}} = r/2$: these are the
management reference points
\item Posterior of $B_{\mathrm{MSY}} = K/2$: the target biomass
\end{itemize}
\subsection{Latent biomass trajectory}
The posterior of $B_{1:T}$ gives a probability band around the
inferred biomass history. For lingcod this should show:
\begin{itemize}
\item Declining biomass through the 1970s--80s
\item Continued low biomass through the 1990s
\item Recovery post-2000 following management reform
\item Posterior uncertainty that widens in data-sparse years
\end{itemize}
Compare the posterior median trajectory to the SS3 assessment biomass
estimates (available from the PFMC assessment document). Systematic
discrepancies reveal the structural differences between an age-structured
model (SS3) and a biomass-aggregated SPM.
\subsection{Current stock status}
The posterior of $B_T / K$ and $B_T / B_{\mathrm{MSY}}$ gives a
full probability distribution over current stock status. Regulatory
thresholds in the US Pacific groundfish fishery are:
\begin{itemize}
\item Overfished if $B_T < 0.25 K$ (PFMC limit reference point)
\item Below MSY target if $B_T < B_{\mathrm{MSY}} = 0.5 K$
\end{itemize}
Compute $P(B_T < 0.25K \mid \text{data})$ and $P(B_T < B_{\mathrm{MSY}}
\mid \text{data})$ directly from the posterior samples. These are
the quantities that enter risk-based harvest decisions and motivate
the CVaR policy layer.
\subsection{Posterior predictive checks}
Generate replicated CPUE series $\tilde I_{1:T}^{(s)}$ from the
posterior predictive distribution. Compare to observed $I_{1:T}$ on:
\begin{itemize}
\item Distribution of residuals $\log I_t - \log \hat I_t$ (should
be approximately Gaussian with mean zero)
\item Autocorrelation of residuals (a positive autocorrelation
signals missing dynamics, such as recruitment variation or
climate-driven productivity shifts)
\item Maximum drawdown: the largest single-year decline in the CPUE
index (tests whether the model can replicate the worst
observed decline)
\end{itemize}
%=============================================================
\section{Diagnostics}
%=============================================================
The full diagnostics checklist from lecture two applies. The
following issues are specific to the SSPM.
\paragraph{Divergences.} Expect divergences near $\sigma_p \approx 0$
(funnel) and near $B_t \approx C_{t-1}$ (constraint boundary). The
non-centred parameterisation resolves the former; the soft floor
resolves the latter.
\paragraph{Slow mixing in $(r, K)$.} The ridge geometry in the
$(r, K)$ posterior slows NUTS mixing even though NUTS is much better
than Metropolis--Hastings. Monitor the ESS for $r$ and $K$
separately; if ESS is below 400 per chain, increase the number of
warmup steps or reduce the step size target acceptance rate from 0.8
to 0.9.
\paragraph{$\hat R$ for latent states.} With $T = 40$ years of data,
the latent state vector $B_{1:T}$ has 40 components. Check $\hat R$
for all of them (ArviZ does this automatically); pay particular
attention to the early years $B_1,\ldots,B_{10}$ where the prior
dominates and mixing can be slow.
\paragraph{Energy diagnostic.} ArviZ's \texttt{plot\_energy()} is
particularly useful for this model. A wide gap between the marginal
energy and the energy transition distribution indicates poor geometry
exploration, typically caused by the non-centred parameterisation
not being applied consistently.
%=============================================================
\section{The Policy Layer: Connection to POMDP and CVaR}
%=============================================================
\subsection{The belief state}
In a POMDP formulation of the fisheries management problem, the
manager observes a noisy CPUE index at each period and must choose a
harvest level without knowing the true biomass. The \emph{belief
state} is the manager's posterior distribution over current biomass
given all past data:
\[
\mu_t \;=\; p(B_t \mid I_{1:t}, C_{1:t-1}, \theta).
\]
In the SSPM framework, this is exactly the filtering distribution
computed by the Bayesian state space model. The terminal posterior
$p(B_T \mid \text{data})$ is the belief state at the end of the data
period, and it is the starting point for forward-looking policy
analysis.
The connection between the SSPM and the POMDP is therefore direct:
run NUTS on the SSPM, extract the posterior samples of $B_T$, and
use those samples as the empirical belief state distribution. No
additional modelling is required. The JAX implementation makes this
particularly clean: the posterior samples are already JAX arrays and
can be passed directly into a JAX-based value function computation.
\subsection{One-period ahead forecast}
Given $M$ posterior samples $\{(r^{(m)}, K^{(m)}, \sigma_p^{(m)},
B_T^{(m)})\}_{m=1}^M$ and a proposed harvest $C_T$, the predictive
distribution of $B_{T+1}$ is:
\begin{equation}
\log B_{T+1}^{(m)} \;\sim\;
\mathcal{N}\!\left(
\log\!\left[B_T^{(m)} + r^{(m)} B_T^{(m)}\!\left(1 - \frac{B_T^{(m)}}{K^{(m)}}\right)
- C_T\right],\;
(\sigma_p^{(m)})^2
\right).
\end{equation}
Draw one sample of $\log B_{T+1}^{(m)}$ for each parameter draw.
The resulting empirical distribution $\{B_{T+1}^{(m)}\}$ is the
one-period-ahead predictive distribution of biomass under harvest $C_T$.
\subsection{CVaR stock collapse risk}
Define stock collapse as the event $\{B_{T+1} < B_{\mathrm{lim}}\}$
where $B_{\mathrm{lim}}$ is a regulatory limit reference point (e.g.\
$0.1K$ or $0.25K$ in the US Pacific groundfish system). The
probability of collapse under harvest $C_T$ is:
\begin{equation}
p_{\mathrm{collapse}}(C_T)
\;=\;
P(B_{T+1} < B_{\mathrm{lim}} \mid \text{data}, C_T)
\;\approx\;
\frac{1}{M}\sum_{m=1}^M \mathbf{1}[B_{T+1}^{(m)} < B_{\mathrm{lim}}].
\end{equation}
The CVaR at level $\alpha$ of the biomass loss is:
\begin{equation}
\mathrm{CVaR}_\alpha(C_T)
\;=\;
\mathbb{E}\!\left[B_{\mathrm{lim}} - B_{T+1} \mid B_{T+1} < q_\alpha,\,
\text{data},\, C_T\right],
\end{equation}
where $q_\alpha$ is the $\alpha$-quantile of $B_{T+1}$ under harvest
$C_T$. This is the expected shortfall of biomass in the worst
$\alpha$ fraction of scenarios.
A \emph{CVaR-constrained} harvest policy solves:
\begin{equation}
\max_{C_T \geq 0} \; \mathbb{E}[C_T]
\quad \text{subject to} \quad
\mathrm{CVaR}_\alpha(C_T) \;\leq\; \delta,
\end{equation}
for a tolerance $\delta$. Since $B_{T+1}^{(m)}$ is differentiable
in $C_T$ (through the deterministic part of the transition), this
optimisation can be solved with JAX gradient methods. The lecture
should present this as an outlook and pointer to the research
programme, not as a solved exercise.
\subsection{The multi-player extension}
The research programme at GRIPS (Pacific saury / NPFC context)
extends this single-player framework to a multi-player game: several
nations share a common stock, each observing the same noisy CPUE
index but choosing harvests independently. The POMDP becomes a
\emph{partially observable stochastic game} (POSG). The belief state
remains the same; the policy layer must account for strategic
interaction. This extension is far beyond the scope of the lecture
but should be mentioned as the eventual destination.
%=============================================================
\section{Suggested Lecture Structure}
%=============================================================
\begin{enumerate}[itemsep=6pt]
\item \textbf{Recap and motivation} (10 min). Recall the logistic
growth state space model from lecture two. Show that adding
a catch term gives the Schaefer surplus production model.
Motivate the application: why do fisheries managers need
biomass estimates with uncertainty?
\item \textbf{The Schaefer model} (15 min). Surplus production,
MSY, $F_{\mathrm{MSY}}$, and $B_{\mathrm{MSY}}$. Plot the
surplus production curve for representative $(r, K)$ values.
Connect to the MSY lecture.
\item \textbf{The state space model} (15 min). Process equation,
observation equation, initial condition. Discuss catch-as-known
assumption. Discuss the positivity constraint and the soft
floor.
\item \textbf{Identification problems} (15 min). The $r$--$K$
ridge, the $q$--$K$ problem, the signal-to-noise ratio.
Show joint prior samples to build intuition before showing
posterior. Introduce the profile-likelihood trick.
\item \textbf{Lingcod data and priors} (10 min). Data pipeline
(already familiar from MSY lecture). Prior choices and
justification from life-history theory. Prior predictive
simulation.
\item \textbf{NumPyro implementation} (20 min). Walk through the
model code. Run NUTS. Show divergences in the centred
parameterisation. Switch to non-centred. Show diagnostics.
\item \textbf{Posterior analysis} (20 min). Marginal posteriors.
Biomass trajectory with credible bands. Current stock
status probabilities. Comparison with SS3 assessment.
Posterior predictive checks.
\item \textbf{The policy layer} (10 min). Show the one-period
predictive distribution of $B_{T+1}$ as a function of $C_T$.
Plot $p_{\mathrm{collapse}}(C_T)$. Sketch the CVaR
formulation. Describe the multi-player extension.
\end{enumerate}
%=============================================================
\section{Notes for a Continuation Session}
%=============================================================
This section records open questions and design choices that have not
yet been resolved, for the benefit of a future working session.
\subsection{Resolved choices}
\begin{itemize}[itemsep=4pt]
\item Dataset: Pacific lingcod from RAM Legacy Database v4.3.
\item Model: Schaefer SSPM (Meyer-Millar 1999 parameterisation).
\item Implementation: NumPyro/JAX with NUTS and non-centred
parameterisation.
\item Catchability: start with analytical profile-out; add full
Bayesian $q$ as an extension.
\item Validation benchmark: PFMC 2021 stock assessment for lingcod.
\end{itemize}
\subsection{Open choices}
\begin{enumerate}[itemsep=6pt]
\item \textbf{Which CPUE series to use.} The RAM Legacy Database
may contain multiple CPUE series for lingcod (commercial
trawl, hook-and-line, survey biomass index). The right
choice should be investigated before writing the lecture
notebook. The PFMC assessment document describes each series.
\item \textbf{Catch series start date.} RAM Legacy catch data for
lingcod begins around 1940. The CPUE index likely starts
later (1980s or later). A decision is needed on whether to
run the process equation from 1940 (long pre-observation
period, driven entirely by catch and prior) or from the first
CPUE observation year. Running from 1940 provides a longer
depletion history but requires the prior on $\phi$ to apply
to 1940 abundance, not the start of the CPUE series.
\item \textbf{Pella-Tomlinson extension.} Whether to include the
shape parameter $m$ in the main model or as an extension
exercise. The extra parameter worsens identification but
may matter for posterior predictive fit.
\item \textbf{lax.scan vs Python loop.} The model should use
\texttt{lax.scan} for the biomass recursion for performance.
The exact NumPyro idiom for non-centred scan needs to be
tested: the interaction between \texttt{lax.scan} and
\texttt{numpyro.sample} inside the scan body has some
version-specific quirks (check NumPyro changelog for
version $\geq$ 0.13).
\item \textbf{Missing CPUE years.} If there are years with no
CPUE observation (survey gaps), the likelihood contribution
for those years should be omitted. The NumPyro implementation
needs to handle a ragged observation set gracefully — either
via masking or by only looping over observed years.
\item \textbf{Comparison with a simple MLE.} It may be worth
including a section showing maximum likelihood estimation of
the same model (optimising the log-likelihood via Optax) as
a baseline, before introducing NUTS. This mirrors the
lecture-one structure (start with MLE, then go Bayesian) and
helps students who are still building intuition for what NUTS
adds.
\item \textbf{Depth of policy layer.} The current plan mentions
CVaR as an outlook section only. A decision is needed on
whether to include executable code for the one-period CVaR
computation or leave it as a mathematical sketch. Including
the code would make the lecture substantially more valuable
as a bridge to the research programme.
\end{enumerate}
%=============================================================
\bibliographystyle{plainnat}
\begin{thebibliography}{99}
\bibitem[Hilborn and Walters(1992)]{HilbornWalters1992}
Hilborn, R.\ and Walters, C.J.\ (1992).
\newblock \textit{Quantitative Fisheries Stock Assessment: Choice, Dynamics
and Uncertainty}.
\newblock Chapman and Hall, New York.
\bibitem[Huynh et al.(2020)]{Huynh2020}
Huynh, Q.C., Carruthers, T.R., and Punt, A.E.\ (2020).
\newblock \texttt{MSEtool}: Stock assessment and management strategy evaluation
toolkit.
\newblock R package. \url{https://github.com/Blue-Matter/MSEtool}.
\bibitem[Meyer and Millar(1999)]{MeyerMillar1999}
Meyer, R.\ and Millar, R.B.\ (1999).
\newblock BUGS in Bayesian stock assessments.
\newblock \textit{Canadian Journal of Fisheries and Aquatic Sciences},
56(6):1078--1087.
\bibitem[Pella and Tomlinson(1969)]{PellaTomlinson1969}
Pella, J.J.\ and Tomlinson, P.K.\ (1969).
\newblock A generalised stock production model.
\newblock \textit{Bulletin of the Inter-American Tropical Tuna Commission},
13(3):419--496.
\bibitem[Punt et al.(2003)]{Punt2003}
Punt, A.E., Smith, A.D.M., and Cui, G.\ (2003).
\newblock In pursuit of a robust harvest rule for Pacific sardine
(\textit{Sardinops sagax}).
\newblock \textit{Fisheries Oceanography}, 12(4--5):379--392.
\bibitem[Schaefer(1954)]{Schaefer1954}
Schaefer, M.B.\ (1954).
\newblock Some aspects of the dynamics of populations important to the
management of commercial marine fisheries.
\newblock \textit{Bulletin of the Inter-American Tropical Tuna Commission},
1(2):27--56.
\bibitem[Stock and Watson(1998)]{StockWatson1998}
Stock, J.H.\ and Watson, M.W.\ (1998).
\newblock Median unbiased estimation of coefficient variance in a time-varying
parameter model.
\newblock \textit{Journal of the American Statistical Association},
93(441):349--358.
\bibitem[Winker et al.(2018)]{Winker2018}
Winker, H., Carvalho, F., and Kapur, M.\ (2018).
\newblock JABBA: Just Another Bayesian Biomass Assessment.
\newblock \textit{Fisheries Research}, 204:275--288.
\bibitem[Methot and Wetzel(2013)]{MethotWetzel2013}
Methot, R.D.\ and Wetzel, C.R.\ (2013).
\newblock Stock Synthesis: A biological and statistical framework for fish
stock assessment and fishery management.
\newblock \textit{Fisheries Research}, 142:86--99.
\bibitem[Rockafellar and Uryasev(2000)]{RockafellarUryasev2000}
Rockafellar, R.T.\ and Uryasev, S.\ (2000).
\newblock Optimization of conditional value-at-risk.
\newblock \textit{Journal of Risk}, 2(3):21--42.
\bibitem[Papaspiliopoulos et al.(2007)]{Papa2007}
Papaspiliopoulos, O., Roberts, G.O., and Sk{\"o}ld, M.\ (2007).
\newblock A general framework for the parametrization of hierarchical models.
\newblock \textit{Statistical Science}, 22(1):59--73.
\end{thebibliography}
\end{document}
Goal
Add a new lecture,
fisheries_sspm.md— Bayesian Estimation of the Schaefer Surplus Production Model: Methodology and Validation — teaching Bayesian estimation of a state-space surplus production model (SSPM) with NumPyro/NUTS, and validating the whole estimation pipeline on synthetic data before any real-world data is touched.This is the first of two fisheries lectures (see the companion application lecture #920, the lingcod assessment). The split is deliberate: this lecture answers "does the method work, and what are its failure modes?"; the companion answers "what does it say about a real stock?". The methodology lecture can be written and fully executed on synthetic data with no external dataset, so it is not blocked on the data task that gates the application lecture.
Position in the sequence
Fourth in the Bayesian state-space arc: unemployment dynamics (#910) → population SSM (#911) → SSPM methodology (this lecture) → SSPM application (lingcod). A fifth lecture on harvest policy (POMDP / risk-sensitive control) is eventual future work, not in scope here.
The model is the stochastic, estimated counterpart of the deterministic Schaefer model taught in the existing intro-series MSY lecture,
msy_fishery.md. There,r, K, qare known and the question is how to fish; here they are unknown and the question is how to infer them — and how uncertain the resulting reference points (MSY = rK/4,B_MSY = K/2,F_MSY = r/2) really are.Origin
Drafted from
fisheries_sspm_lecture.tex(third in John's "Notes Toward a QuantEcon Lecture" series; full source embedded below). The core model — the Meyer–Millar (1999) state-space Schaefer model — is sound and standard. The.texis a spec, not a draft, and bundles an estimation lecture and a policy/POMDP/CVaR layer; per discussion the policy material is deferred and the estimation material is split across two lectures.Review findings (correctness checklist — fix when drafting)
The
.texcontains several technical errors that must be corrected in the lecture:m=2as written. Eq. (3),P(B) = (r/m)·B·[1 − (B/K)^{m-1}], gives(r/2)·B·(1 − B/K)atm=2— half the Schaefer eq. (1)rB(1 − B/K). Either the prefactor is wrong orris being silently rescaled; state the reparameterisation explicitly. (Pella–Tomlinson is an extension/exercise only, so this can be deferred — but fix it if included.)C_{t-1}(growth-then-removal). Make prose and equation agree.q= analytically marginalising over a flat prior onlog q" is false. Profiling plugs in the conditional MLEq̂; marginalising integratesqout. For this Gaussian-in-log qmodel the two differ by a factor∝ σ_o/√T, which does depend onσ_o— one of the weakly-identified parameters this lecture is about. Present profiling honestly as an approximation (it also costs one degree of freedom:Tresiduals constrained to sum to zero, fed toNormal(0, σ_o)), not as exact marginalisation. (q-profiling is an optional extension.)q–Kproblem" overstates. It removesqas a parameter, but by discarding the level information;Kis then identified only by trajectory shape. Don't claim the identifiability is recovered..texwrites the signal-to-noise ratio asq_σ = σ_p²/σ_o², reusingq(catchability). Rename (e.g.ρ), consistent with the fix already agreed in [population_ssm] New lecture: Bayesian state-space models (linear filtering to nonlinear population dynamics) #911.numpyro.contrib.control_flow.scan, not rawjax.lax.scan(same fix as [population_ssm] New lecture: Bayesian state-space models (linear filtering to nonlinear population dynamics) #911). Use the non-centred parameterisation from the start.r ~ HalfNormal(0.5)has its mode at 0 — prefer the LogNormal alternative the.texalready mentions.Agreed design
θ, simulateB_{1:T}andI_{1:T}, run the full NUTS pipeline; check marginal posteriors cover the truth, the latent-biomass band contains the true path, and derived reference points (MSY,B_MSY,F_MSY) recover. Design touch: use the real lingcod catch seriesC_t(known anyway) but simulate the index from knownθ, so the experiment matches the real fit in scale/timing.log I = log q + log B + εis linear-Gaussian in(log q, log B), so recoveringqis trivial and the sampler must nail it exactly. Quick plumbing check before the hard nonlinear fit.r–Kridge appears even when the truth is known (proves it's data geometry, not misspecification);σ_p, σ_oand show the posterior ofσ_p²/σ_o²is diffuse — a controlled demonstration of the signal-to-noise non-identification, rather than an assertion.B_{2:T}withB_1 = φK; soft floormax(A_t, εK)on the log argument, with the trade-off explained.Section plan
r, K, q; here we estimate them." (brief)B_{1:T}.r–Kridge,q–Klevel confounding, signal-to-noise (ρ).numpyro.contrib.control_flow.scan; soft floor.σ_p²/σ_o².R̂/ESS/divergences; centred-vs-non-centred funnel;plot_energy.Closes on: "we've earned the right to trust this method and we understand its failure modes — next we apply it to a real stock."
Optional/exercises (keep core lean): Pella–Tomlinson
m; profiling-outq; multiple CPUE indices; an MLE-via-Optax warm-up (mirrors the unemployment lecture's MLE-then-Bayes arc).Cross-referencing & style guide (note for eventual alignment)
lecture-python.myst): standard{ref}/{doc}tokalman.md,ar1_bayes.md,bayes_nonconj.md, and the plannedpopulation_ssm.md([population_ssm] New lecture: Bayesian state-space models (linear filtering to nonlinear population dynamics) #911) /unemployment_bayes.md([unemployment_bayes] New lecture: Bayesian estimation of nonlinear unemployment dynamics #910).lecture-python-intro: plain full URL,https://intro.quantecon.org/msy_fishery.html(the repos use full-URL links across series — there is no intersphinx).style-guiderepo is the successor to the olderQuantEcon.manual(which documents the cross-referencing conventions). The manual/style-guide conventions — including cross-referencing setup — should eventually be applied across this whole Bayesian sequence; flagging here so we bring everything into line in a later pass.Tasks
C_t+ simulated index); coverage of 90% CIsq-recovery unit testr–Kridge under known truth (short vs informative series); diffuseσ_p²/σ_o²lectures/fisheries_sspm.mdin MyST (GPU admonition include; numpyro/jax/arviz imports), with all correctness fixes abovelectures/_toc.ymlnear the Bayesian lectureskalman.md,ar1_bayes.md,bayes_nonconj.md,population_ssm.md,unemployment_bayes.md, and (full URL) the MSY lectureGating: conceptually follows #910/#911 (shared machinery), though it can be drafted in parallel since it is self-contained on synthetic data.
Original draft source —
fisheries_sspm_lecture.tex(Embedded for self-containment. John may also drag-drop the converted markdown onto this issue.)