In cboettig/multiple_uncertainty: Multiple Uncertainty
Abstract
We explore the role of multiple forms of uncertainty on optimal management policies
Introduction
-
Influence of Reed 1979 in our understanding of management under uncertainty
-
Importance of multiple uncertainty
-
The surprising sparse history of approaches (at least in optimal fisheries)
-
Context of partial observability in other areas, e.g. POMDP (cite Chades papers, Springborn paper)
Model and Methods
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE, cache = TRUE)
library("dplyr")
library("tidyr")
library("ggplot2")
library("multidplyr")
library("multipleuncertainty")
Model formulation
Our model formulation closely parallels @Sethi2005, which extends
the classic formulation of @Reed1979 to consider implementation and
measurement uncertainty. Population growth of the stock, $x_t$,
is determined by a stochastic stock $z_t^g$ to the growth function
$G(s_{t-1})$ of the previous year's escaped stock $s_{t-1}$,
$$ x_t = z_t^g G(s_{t-1}) $$
where escapement of a given year is the true stock minus what is harvested,
$$s_t = x_t - h_t$$
Measured stock, $m$ is determined by a random shock $z_t^m$ to the true stock size $x_t$,
$$ m_t = z_t^m x_t $$
The harvest size is determined by a stochastic stock $z_t^i$ to the
quota $q_t$ set by management policy, or the true stock size
(whichever is smaller)
$$ h_t = \min(x_t, z_t^i q_t) $$
For discount rate $\alpha$ and utility function $U(m_t, q_t(m_t))$,
the fisheries management problem can be written as the choice
(policy) of quotas ${ q_t }$ for each year that maximize the
expected total reward:
$$ \max_{{ q_t} \geq 0} \mathbb{E} \left[ \sum_0^{\infty} \alpha^t U \left(m_t, q_t(m_t) \right) \right] $$
At each year $t$ the manager measures the stock (with error) $m_t$,
and must then set a quota $q_t$. The optimal solution can thus be
found by dynamic orogramming after defining a Bellman recusion
equation in terms of the observed states $m_t$ and quotas $q_t$:
$$ V_t(m_t) = \max_{q_t} \left[ U_t(m_t, q_t) + \alpha V_{t+1}(m_t) \right] \label{bellman} $$
The trick is simply to define the appropriate utility for each
possible (observed) state $m_t$ and action (quota) $q_t$, and the
appropriate transition probability $P(m_{t+1} | m_t, q_t)$ for each
possible quota $q_t$ and state $m_t$. In this model, the transition
probability will depend on the growth model and all three sources
of uncertainty.
The expected utility is simpler to construct. The utility $U_t$
of choosing a quota $q_t$ having measured state $m_t$ is given by
the rewards derived from realizing a harvest $h_t$, integrated over
the probability that a harvest $h_t$ would be realized by a measurement
$m_t$ and quota $q_t$:
$$ U(m_t, q_t) = \int_x P(x_t | m_t) \int_h P(h_t | q_t) R(x_t,h_t) dh_t dx_t \label{utility} $$
Computing the probability $P(h|q)$ is straight forward, as it follows
immediately from the distribution from which we choose shocks $h_t
= z_t^i q_t$. However, the probability $P(x_t | m_t)$ is more
subtle[^1], as we have thus far only defined the inverse relation,
$P(m_t | x_t)$ by defining the shock $z_m^t$ in the expression $m_t
= z_m^t x_t$. The conditional probability can be reversed with the
help of Bayes Law:
$$P(x_t | m_t ) = \frac{P(m_t | x_t) P(x_t)}{\int P(m_t | x_t) P(x_t) dx_t,}$$
assuming a simple but naive choice of prior belief $P(x_t)$, a
uniform prior. While it would be preferable to define a prior that
was conditional on the previous measurements $m_{{t-1}\dots m_0}$,
this would forefit the Markovian assumption need to make the problem
amenable to solution. Such a calculation is not ammenable to a stochastic dynamic programming
approach, in which transition probabilities must be defined only
in terms of the current state and would thus have to consider all
possible previous states $m_{t-1}$, as well as current states $m_t$,
in determining the optimal quota $q_t$. The increased dimension
problem is beyond our scope. The uniform prior neatly sidesteps
this.[^2] Though @Sethi2005 gloss over any detail about how $P(x_t | m_t)$
is constructed, confirms that they have
taken exactly this approach.
[^1]: Introducing measurement uncertainty into Markov Decision Processes
is the focus of Partially Observed Markov Decision Processes (POMDP) research,
an active area in artificial intelligence community. Solutions do not use dynamic
programming algorithm, relying instead on a variety of approximation algorithms
which the mathematical biology community would recognize more as Expectation-Maximization (EM)
algorithms for Hidden Markov Models. These models permit active or online learning
about the process, rather than the single optimal policy functions considered here.
Their analysis is beyond the scope of this paper.
[^2]: It is also possible to demonstrate that the qualitative
patterns shown here are not sensitive to the precise choice of the
resulting probability distribution for $P( x_t | m_t)$, see
results.
We compare these results to what we would find using the respective
inverse probability function (inverse uniform, inverse log-normal),
which compared to the approach discussed above, puts relatively
less weight on stock sizes larger than the observed measurement and
more weight on smaller values. Note that this only impacts
calculations involving measurement error (not implementation error).
Numerical implementation
We solve the optimization problem using Stochastic Dynamic Programming, by first
discritizing the problem described above and then determining the resluting utility
matrix, $\mathcal{U}{mq}$ and transition probability between measured states for
each possible choice of quota, $\mathcal{T}{m_t, m_{t+1}, q_t}$. Given this matrix
and tensor, the optimization can be solved through the standard policy iteration algorithm.
We provide complete R code and Matlab/Octave code implementations for readers interested
in the details or in running the algorithm and exploring on their own.
We discritize the state space
(stock $x_t$, observed stock $m_t$, harvest $h_t$, quota $q_t$) onto
identical uniform grids between values between 0 and 200 at 0.5 unit step size.
(In general each dimension could have a different discritization; for instance,
allowing the measured stock sizes and/or chosen quotas to limited to a coarser
grid than the true dynamics, which could reflect real-world limitations and
not merely numerical necessity.)
The Supplemental material compares our results across a range of grid
lengths and step sizes to confirm this choice is sufficient to minimize
numerical errors introduced by the discritization.
Note that @Sethi2005 perform a similar discritization but present only cubic-spline interpolations
of the resulting policies. While this can supress numerical artifacts of the
discritization, it can also merely blurr them or introduce additional artifacts.
We present the policies as determined by discrete optimization without smoothing
by splines, making interpretation more straight-forward.
In the discrete problem, the Utility function $U(m_t, q_t)$ can be written as a
matrix $\mathcal{U}{mq}$, of size length(quota grid) by length(measurement grid). To do so, we
begin by defining the matrix $\mathcal{R}{xh}$ for the reward function $R(x_t, h_t)$ for harvest
of $h_t$ at stock size $x_t$ for each of their respective grid point values.
We can then define a matrix $\mathcal{P}{hq}$ as the probability of a harvest $h_t$ given a
quota set at $q_t$ from the assumed probability distribution $P(h_t | q_t)$.
Likewise given the probability of a true stock size of $x_t$ given a measurement of $m_t$,
$P(x_t | m_t)$ we define the probability matrix $\mathcal{P}{xm}$. Then the matrix products
$$\mathcal{U}{mq} = \mathcal{P}{xm} \mathcal{R}{xh} \mathcal{P}{hq}$$
are the discrete equivalent of integrating over $h$ and $x$ in Eqn~\eqref{utility},
giving us the utility matrix $\mathcal{U}_{mq}$.
The transition probability tensor is likewise defined in the space of the
manager (measured stock, quota), which must also be computed from the
true state of the system (stock, harvest).
Results
Figure 1 compares the optimal policies under the assumption of
uniform noise and lognormal noise, using a logistic growth model
with $r = 1$, $K = 100$, on a grid of $x \in [0,200]$ with interval
$\Delta = 0.5$, all following @Sethi2005. As in
@Sethi2005, the discounting rate is set at 5%, the reward function
places a fixed price per unit of fish harvested, with no cost to
harvesting effort. (The appendices explore the robustness of these
results to variation in each of these choices.)
fig3 <- function(noise){
grid <- seq(0, 200, by=0.5)
## Scale noise comparably for log-normal vs uniform
if(noise == "lognormal"){
lo <- 0.0577 # 0.1 / sqrt(3)
hi <- 0.2887 # 0.5 / sqrt(3)
} else {
lo <- 0.1
hi <- 0.5
}
small <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise)
growth <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise)
measure <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise)
implement <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise)
large <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = hi, sigma_m = hi, sigma_i = hi, noise_dist = noise)
df <- data.frame(y_grid = grid, small = small, growth = growth,
measure = measure, implement = implement, large = large) %>%
tidyr::gather(scenario, value, -y_grid)
}
cluster <- create_cluster(2)
df <-
data.frame(noise = c("uniform", "lognormal")) %>%
dplyr::group_by(noise) %>%
multidplyr::partition(noise, cluster = cluster) %>%
dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) +
geom_line() +
facet_wrap(~ noise) +
xlab("Measured Stock Size") +
ylab("Proposed Escapement") +
coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) +
theme_bw()
The left panel of Figure 1 is thus almost directly comparable @Sethi2005's Figure 3,
Our plot differs in three regards: First, we plot
the optimal policy for measured stock values that exceed the carrying capacity
(measured stock sizes up to 150), instead of cutting off the policy at carrying
capacity (stock size of 100). In the absence of measurement error the policy
above carrying capacity is of little interest, since such stock sizes are rarely
realized. Once we incorporate measurment error, observations of values in this
region are not at all unlikely (even if the true stock size is consistently below
carrying capacity), and thus some useful insight can be found by examining these
policies. Second, as described in the methods, we have ommitted the spline smoothing
step to show the policy that results directly from the discrete optimization, and
thus reveals the numerical discontinuities that result. Note that the impact of
discritization is more acute in the case of uniform noise distribution, which itself is
not smooth and may be split in-between grid boundaries. Finally, we have added the additional
curve showing the case when all noise sources are set to the high variance levels to permit
further comparison of the interactions between noise sources. The right panel shows
identical conditions but for the replacement of the uniform noise assumption
with lognormal noise of equivalent variance.
In Figure 1, all sources of noise are present in each plot at "low" levels, while varying
which noise source is present at a "high" level. While this permits direct comparison
to @Sethi2005, it is often instructive to consider each form of uncertainty individually
while removing the other sources of uncertainty. Such a plot is shown in Figure 2.
Observe that the deterministic and growth-only solutions now recover the
constant-escapement rule exactly, as expected.
fig3 <- function(noise){
grid <- seq(0, 200, by=0.5)
## Scale noise comparably for log-normal vs uniform
if(noise == "lognormal"){
hi <- 0.2887 # 0.5 / sqrt(3)
} else {
hi <- 0.5
}
lo <- 0.0
small <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise)
growth <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise)
measure <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise)
implement <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise)
df <- data.frame(y_grid = grid, none = small, growth = growth,
measure = measure, implement = implement) %>%
tidyr::gather(scenario, value, -y_grid)
}
df <-
data.frame(noise = c("uniform", "lognormal")) %>%
dplyr::group_by(noise) %>%
dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) +
geom_line() +
facet_wrap(~ noise) +
xlab("Measured Stock Size") +
ylab("Proposed Escapement") +
coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) +
theme_bw()
The most salient result seen in Figures 1 & 2 is that both implementation
error and measurement error can both result in substantial deviations
from the classic "Constant Escapment" result of @Reed1979. Such a result
also stands in contrast to the conclusions of @Sethi2005, which emphasize
that only measurement error can be responsible for significant deviations.
Understanding how the optimal strategy differs from the classic constant escapement
solution is also of interest. In particular, it is interesting to ask if or when accounting
for these additional uncertainties will make the resulting policy more conservative or less
conservative than the policy based on growth noise alone. Recall that @Reed1979 proves that
accounting for uncertainty in the growth process can only make the solution more conservative,
but not less: $S >= D$, where $S$ is the constant escapement level of the stochastic model and
$D$ that of the deterministic model. Moreover, Reed showed that for sufficiently small noise,
$S=D$.
The impacts of measurement and implementation uncertainties are more complex: sometimes
more conservative, sometimes more agressive. In general, very large measured stock
sizes are treated more conservatively after accounting for measurement error than they
would be if that error were ignored. For stock sizes larger than carrying capacity,
this is intuitive: such values are much more likely to be the result of measurement error.
Note that the same conservative pattern holds true for implementation error, even when
the measured stock sizes are known exactly (Figure 2, implementation error policy). In
this case, the inituition is different but equally straight forward: because implementation
error is multiplicative, large harvests come with much greater variance than small harvests.
For instance, under the uniform noise assumption a quota set at 100 is just as likely to harvest
150 as it is to harvest 100 units.
A similar intuition operates at the other extreme. The combination of population growth and
the harvest policy push the stock towards the optimal escapement level. Therefore, Measurements of
stock sizes significantly below that target escapement are most likely to result from measurement
error (when present) underestimating the stock size. Consequently, the quota set for small stock
sizes is higher (thus the escapement is lower) under the assumption of measurment error than without it.
The intepretation for implementation error is more subtle, and indeed the effect is less extreme. Just
as for very high stock sizes the proposed quotas were very large and very uncertain, for low stock sizes
the quotas are very small. In fact, for the deterministic case, quotas are zero below constant escapement
level. If we add only implementation noise (Figure 2), below the deterministic constant escapement level
there is no reason to deviate from the deterministic policy: a quota of zero gets a harvest of exactly zero.
As the stock sizes get just above this, harvesting begins at small amounts. If only implementation error
is present, it makes sense to harvest slightly more (setting escapement slightly lower) than in the deterministic
case because the value (dock-side revenue) from under-harvesting a little bit is less than from over-harvesting
due to the discounting and concave nature of the stock-recruitment curve. (Knocking the population down a little too low
means you derive the immediate benefit of the higher harvest, at the cost of the stock next year being lower than it would
be a the optimal harvest. Fishing too little imposes an immediate cost of less harvest than optimal, and a slight
benefit of higher-than-expected stock the following year. However, future benefits are discounted, and growth
rate at higher population densities slows).
Implementation error can also interact with other forms of uncertainty to result in a more aggressive
harvest even for stock values measured below the deterministic constant escapement. By comparing the
results of all uncertainties set to "small" noise in Figure 1 to implementation uncertainty being large
(while others are small but non-zero), we see that non-zero quotas begin even before they would under
the same amount of measurement error but smaller implementation error (i.e. the more agressive quota
seen for measured stock sizes below the deterministic escapement is not driven by measurement uncertainty alone.)
At small stock sizes, the possibility of bumper growth (from the stochastic recruitment process) and the possibility of understimated
stock sizes make it worthwhile to set a more agressive quota when implementation is very uncertain than when it is
relatively precise.
sigma_g = 0.1
sigma_m = 0.5
sigma_i = 0.1
f = "logistic"
noise_dist = "uniform"
delta <- 0.05
grid <- seq(0,250,length=401)
S <- multiple_uncertainty(f = f,
x_grid = grid,
sigma_g = sigma_g,
sigma_m = sigma_m,
sigma_i = sigma_i,
delta = delta,
noise_dist = noise_dist)
sims <- function(i)
scenario(S, f = f,
x_grid = grid,
x0= 75,
Tmax = 500,
sigma_g = sigma_g,
sigma_m = sigma_m,
sigma_i = sigma_i,
noise_dist = noise_dist)
sim <- purrr::map_df(1:100, sims)
sim %>%
mutate(s = x - h) %>%
select(t, x, y, s) %>%
tidyr::gather(variable, value, -t) %>%
ggplot() + geom_density(aes(value, fill=variable), alpha=0.2)
While such arguments help us understand the edge behavior, only the precise optimal policy calculation
can deterime exactly where the policies cross -- that is, where increased uncertainties results in a more
conservative or more agressive fishing quota. Under uniform noise, we see a stark difference between only
implementation uncertainty and only measurement uncertainty in the range of measured stock sizes that is of greatest
interest (e.g. most likely to be observed): As measurement uncertainty increases, the policy is consistently more agressive (higher quotas, lower escapement)
across this entire range, while as implementation uncertainty increases the policy becomes more conservative (lower quotas, higher escapement).
This is by no means a generic result. Under lognormal noise, the results of measurement vs implementation
uncertainty are far more similar. Policy under implementation error is always more conservative than an equivalent magnitude
of measurement error, all else equal, but measurement error crosses the deterministic policy far sooner.
If we consider much smaller growth rates, the pattern is dramatically different:
fig3 <- function(noise){
grid <- seq(0, 200, by=0.5)
## Scale noise comparably for log-normal vs uniform
if(noise == "lognormal"){
lo <- 0 #0.0577 # 0.1 / sqrt(3)
hi <- 0.2887 # 0.5 / sqrt(3)
} else {
lo <- 0.0
hi <- 0.5
}
model <- function(x, h) logistic(x, h, r = 0.1)
small <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist)
growth <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist)
measure <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise_dist)
implement <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise_dist)
df <- data.frame(y_grid = grid, small = small, growth = growth,
measure = measure, implement = implement) %>%
tidyr::gather(scenario, value, -y_grid)
}
df <-
data.frame(noise = c("uniform", "lognormal")) %>%
dplyr::group_by(noise) %>%
dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) +
geom_line() +
facet_wrap(~ noise) +
xlab("Measured Stock Size") +
ylab("Proposed Escapement") +
coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) +
theme_bw()
If we use a compensating instead of overcompensating recruitment function such as Beverton-Holt, the pattern is also quite different.
fig3 <- function(noise){
grid <- seq(0, 200, by=0.5)
## Scale noise comparably for log-normal vs uniform
if(noise == "lognormal"){
lo <- 0 #0.0577 # 0.1 / sqrt(3)
hi <- 0.2887 # 0.5 / sqrt(3)
} else {
lo <- 0.0
hi <- 0.5
}
model <- "bevertonholt"
small <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist)
growth <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist)
measure <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise_dist)
implement <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise_dist)
df <- data.frame(y_grid = grid, small = small, growth = growth,
measure = measure, implement = implement) %>%
tidyr::gather(scenario, value, -y_grid)
}
df <-
data.frame(noise = c("uniform", "lognormal")) %>%
dplyr::group_by(noise) %>%
dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) +
geom_line() +
facet_wrap(~ noise) +
xlab("Measured Stock Size") +
ylab("Proposed Escapement") +
coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) +
theme_bw()
The magnitude of the noise makes some difference as well, with larger magnitudes driving earlier crossing points.
Discussion
Review of general lit
Consistent with Sethi et al, we observe:
- Small uncertainties lead to only small deviations from classical expectations based on the deterministic model or growth-noise only model.
- Large measurement uncertainty can lead to a significant deviation from "constant escapement"
In contrast to Sethi et al, we find that:
-
These observations are sensitive to way uncertainty is modelled, rather than being generally robust as claimed.
-
Large implementation uncertainty can also result in substantial deviations from the the "constant escapement" rule.
This is most pronounced in more realisitic models of uncertainty, such as log-normally distributed noise, but visible in uniform noise as well. Small deviations from constant escapement are visible in the implementation error results of @Sethi2005, but are understated. This contrasts to the @Sethi claim that
... implementation sources of uncertainty is high, a constant-escapement policy is qualitatively appropriate since the slopes of these policies beyond the shutdown point are virtually flat.
Substantial deviations from the optimal policy can be found for non-convex growth functions. In particular, when there is small or negligible cost to harvest, the optimal quota can signficantly exceed the measured stock size, as the true stock size and true harvest can differ substantially.
cboettig/multiple_uncertainty documentation built on May 13, 2019, 2:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Abstract
We explore the role of multiple forms of uncertainty on optimal management policies
Introduction
-
Influence of Reed 1979 in our understanding of management under uncertainty
-
Importance of multiple uncertainty
-
The surprising sparse history of approaches (at least in optimal fisheries)
-
Context of partial observability in other areas, e.g. POMDP (cite Chades papers, Springborn paper)
Model and Methods
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE, cache = TRUE)
library("dplyr") library("tidyr") library("ggplot2") library("multidplyr") library("multipleuncertainty")
Model formulation
Our model formulation closely parallels @Sethi2005, which extends the classic formulation of @Reed1979 to consider implementation and measurement uncertainty. Population growth of the stock, $x_t$, is determined by a stochastic stock $z_t^g$ to the growth function $G(s_{t-1})$ of the previous year's escaped stock $s_{t-1}$,
$$ x_t = z_t^g G(s_{t-1}) $$
where escapement of a given year is the true stock minus what is harvested,
$$s_t = x_t - h_t$$
Measured stock, $m$ is determined by a random shock $z_t^m$ to the true stock size $x_t$,
$$ m_t = z_t^m x_t $$
The harvest size is determined by a stochastic stock $z_t^i$ to the quota $q_t$ set by management policy, or the true stock size (whichever is smaller)
$$ h_t = \min(x_t, z_t^i q_t) $$
For discount rate $\alpha$ and utility function $U(m_t, q_t(m_t))$, the fisheries management problem can be written as the choice (policy) of quotas ${ q_t }$ for each year that maximize the expected total reward:
$$ \max_{{ q_t} \geq 0} \mathbb{E} \left[ \sum_0^{\infty} \alpha^t U \left(m_t, q_t(m_t) \right) \right] $$
At each year $t$ the manager measures the stock (with error) $m_t$, and must then set a quota $q_t$. The optimal solution can thus be found by dynamic orogramming after defining a Bellman recusion equation in terms of the observed states $m_t$ and quotas $q_t$:
$$ V_t(m_t) = \max_{q_t} \left[ U_t(m_t, q_t) + \alpha V_{t+1}(m_t) \right] \label{bellman} $$
The trick is simply to define the appropriate utility for each possible (observed) state $m_t$ and action (quota) $q_t$, and the appropriate transition probability $P(m_{t+1} | m_t, q_t)$ for each possible quota $q_t$ and state $m_t$. In this model, the transition probability will depend on the growth model and all three sources of uncertainty.
The expected utility is simpler to construct. The utility $U_t$ of choosing a quota $q_t$ having measured state $m_t$ is given by the rewards derived from realizing a harvest $h_t$, integrated over the probability that a harvest $h_t$ would be realized by a measurement $m_t$ and quota $q_t$:
$$ U(m_t, q_t) = \int_x P(x_t | m_t) \int_h P(h_t | q_t) R(x_t,h_t) dh_t dx_t \label{utility} $$
Computing the probability $P(h|q)$ is straight forward, as it follows immediately from the distribution from which we choose shocks $h_t = z_t^i q_t$. However, the probability $P(x_t | m_t)$ is more subtle[^1], as we have thus far only defined the inverse relation, $P(m_t | x_t)$ by defining the shock $z_m^t$ in the expression $m_t = z_m^t x_t$. The conditional probability can be reversed with the help of Bayes Law:
$$P(x_t | m_t ) = \frac{P(m_t | x_t) P(x_t)}{\int P(m_t | x_t) P(x_t) dx_t,}$$
assuming a simple but naive choice of prior belief $P(x_t)$, a uniform prior. While it would be preferable to define a prior that was conditional on the previous measurements $m_{{t-1}\dots m_0}$, this would forefit the Markovian assumption need to make the problem amenable to solution. Such a calculation is not ammenable to a stochastic dynamic programming approach, in which transition probabilities must be defined only in terms of the current state and would thus have to consider all possible previous states $m_{t-1}$, as well as current states $m_t$, in determining the optimal quota $q_t$. The increased dimension problem is beyond our scope. The uniform prior neatly sidesteps this.[^2] Though @Sethi2005 gloss over any detail about how $P(x_t | m_t)$ is constructed, confirms that they have taken exactly this approach.
[^1]: Introducing measurement uncertainty into Markov Decision Processes is the focus of Partially Observed Markov Decision Processes (POMDP) research, an active area in artificial intelligence community. Solutions do not use dynamic programming algorithm, relying instead on a variety of approximation algorithms which the mathematical biology community would recognize more as Expectation-Maximization (EM) algorithms for Hidden Markov Models. These models permit active or online learning about the process, rather than the single optimal policy functions considered here. Their analysis is beyond the scope of this paper.
[^2]: It is also possible to demonstrate that the qualitative patterns shown here are not sensitive to the precise choice of the resulting probability distribution for $P( x_t | m_t)$, see results. We compare these results to what we would find using the respective inverse probability function (inverse uniform, inverse log-normal), which compared to the approach discussed above, puts relatively less weight on stock sizes larger than the observed measurement and more weight on smaller values. Note that this only impacts calculations involving measurement error (not implementation error).
Numerical implementation
We solve the optimization problem using Stochastic Dynamic Programming, by first discritizing the problem described above and then determining the resluting utility matrix, $\mathcal{U}{mq}$ and transition probability between measured states for each possible choice of quota, $\mathcal{T}{m_t, m_{t+1}, q_t}$. Given this matrix and tensor, the optimization can be solved through the standard policy iteration algorithm. We provide complete R code and Matlab/Octave code implementations for readers interested in the details or in running the algorithm and exploring on their own.
We discritize the state space (stock $x_t$, observed stock $m_t$, harvest $h_t$, quota $q_t$) onto identical uniform grids between values between 0 and 200 at 0.5 unit step size. (In general each dimension could have a different discritization; for instance, allowing the measured stock sizes and/or chosen quotas to limited to a coarser grid than the true dynamics, which could reflect real-world limitations and not merely numerical necessity.)
The Supplemental material compares our results across a range of grid lengths and step sizes to confirm this choice is sufficient to minimize numerical errors introduced by the discritization.
Note that @Sethi2005 perform a similar discritization but present only cubic-spline interpolations of the resulting policies. While this can supress numerical artifacts of the discritization, it can also merely blurr them or introduce additional artifacts. We present the policies as determined by discrete optimization without smoothing by splines, making interpretation more straight-forward.
In the discrete problem, the Utility function $U(m_t, q_t)$ can be written as a matrix $\mathcal{U}{mq}$, of size length(quota grid) by length(measurement grid). To do so, we begin by defining the matrix $\mathcal{R}{xh}$ for the reward function $R(x_t, h_t)$ for harvest of $h_t$ at stock size $x_t$ for each of their respective grid point values. We can then define a matrix $\mathcal{P}{hq}$ as the probability of a harvest $h_t$ given a quota set at $q_t$ from the assumed probability distribution $P(h_t | q_t)$. Likewise given the probability of a true stock size of $x_t$ given a measurement of $m_t$, $P(x_t | m_t)$ we define the probability matrix $\mathcal{P}{xm}$. Then the matrix products
$$\mathcal{U}{mq} = \mathcal{P}{xm} \mathcal{R}{xh} \mathcal{P}{hq}$$
are the discrete equivalent of integrating over $h$ and $x$ in Eqn~\eqref{utility}, giving us the utility matrix $\mathcal{U}_{mq}$.
The transition probability tensor is likewise defined in the space of the manager (measured stock, quota), which must also be computed from the true state of the system (stock, harvest).
Results
Figure 1 compares the optimal policies under the assumption of uniform noise and lognormal noise, using a logistic growth model with $r = 1$, $K = 100$, on a grid of $x \in [0,200]$ with interval $\Delta = 0.5$, all following @Sethi2005. As in @Sethi2005, the discounting rate is set at 5%, the reward function places a fixed price per unit of fish harvested, with no cost to harvesting effort. (The appendices explore the robustness of these results to variation in each of these choices.)
fig3 <- function(noise){ grid <- seq(0, 200, by=0.5) ## Scale noise comparably for log-normal vs uniform if(noise == "lognormal"){ lo <- 0.0577 # 0.1 / sqrt(3) hi <- 0.2887 # 0.5 / sqrt(3) } else { lo <- 0.1 hi <- 0.5 } small <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise) growth <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise) measure <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise) implement <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise) large <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = hi, sigma_m = hi, sigma_i = hi, noise_dist = noise) df <- data.frame(y_grid = grid, small = small, growth = growth, measure = measure, implement = implement, large = large) %>% tidyr::gather(scenario, value, -y_grid) }
cluster <- create_cluster(2) df <- data.frame(noise = c("uniform", "lognormal")) %>% dplyr::group_by(noise) %>% multidplyr::partition(noise, cluster = cluster) %>% dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) + geom_line() + facet_wrap(~ noise) + xlab("Measured Stock Size") + ylab("Proposed Escapement") + coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) + theme_bw()
The left panel of Figure 1 is thus almost directly comparable @Sethi2005's Figure 3, Our plot differs in three regards: First, we plot the optimal policy for measured stock values that exceed the carrying capacity (measured stock sizes up to 150), instead of cutting off the policy at carrying capacity (stock size of 100). In the absence of measurement error the policy above carrying capacity is of little interest, since such stock sizes are rarely realized. Once we incorporate measurment error, observations of values in this region are not at all unlikely (even if the true stock size is consistently below carrying capacity), and thus some useful insight can be found by examining these policies. Second, as described in the methods, we have ommitted the spline smoothing step to show the policy that results directly from the discrete optimization, and thus reveals the numerical discontinuities that result. Note that the impact of discritization is more acute in the case of uniform noise distribution, which itself is not smooth and may be split in-between grid boundaries. Finally, we have added the additional curve showing the case when all noise sources are set to the high variance levels to permit further comparison of the interactions between noise sources. The right panel shows identical conditions but for the replacement of the uniform noise assumption with lognormal noise of equivalent variance.
In Figure 1, all sources of noise are present in each plot at "low" levels, while varying which noise source is present at a "high" level. While this permits direct comparison to @Sethi2005, it is often instructive to consider each form of uncertainty individually while removing the other sources of uncertainty. Such a plot is shown in Figure 2. Observe that the deterministic and growth-only solutions now recover the constant-escapement rule exactly, as expected.
fig3 <- function(noise){ grid <- seq(0, 200, by=0.5) ## Scale noise comparably for log-normal vs uniform if(noise == "lognormal"){ hi <- 0.2887 # 0.5 / sqrt(3) } else { hi <- 0.5 } lo <- 0.0 small <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise) growth <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise) measure <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise) implement <- multiple_uncertainty(f = logistic, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise) df <- data.frame(y_grid = grid, none = small, growth = growth, measure = measure, implement = implement) %>% tidyr::gather(scenario, value, -y_grid) } df <- data.frame(noise = c("uniform", "lognormal")) %>% dplyr::group_by(noise) %>% dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) + geom_line() + facet_wrap(~ noise) + xlab("Measured Stock Size") + ylab("Proposed Escapement") + coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) + theme_bw()
The most salient result seen in Figures 1 & 2 is that both implementation error and measurement error can both result in substantial deviations from the classic "Constant Escapment" result of @Reed1979. Such a result also stands in contrast to the conclusions of @Sethi2005, which emphasize that only measurement error can be responsible for significant deviations.
Understanding how the optimal strategy differs from the classic constant escapement solution is also of interest. In particular, it is interesting to ask if or when accounting for these additional uncertainties will make the resulting policy more conservative or less conservative than the policy based on growth noise alone. Recall that @Reed1979 proves that accounting for uncertainty in the growth process can only make the solution more conservative, but not less: $S >= D$, where $S$ is the constant escapement level of the stochastic model and $D$ that of the deterministic model. Moreover, Reed showed that for sufficiently small noise, $S=D$.
The impacts of measurement and implementation uncertainties are more complex: sometimes more conservative, sometimes more agressive. In general, very large measured stock sizes are treated more conservatively after accounting for measurement error than they would be if that error were ignored. For stock sizes larger than carrying capacity, this is intuitive: such values are much more likely to be the result of measurement error.
Note that the same conservative pattern holds true for implementation error, even when the measured stock sizes are known exactly (Figure 2, implementation error policy). In this case, the inituition is different but equally straight forward: because implementation error is multiplicative, large harvests come with much greater variance than small harvests. For instance, under the uniform noise assumption a quota set at 100 is just as likely to harvest 150 as it is to harvest 100 units.
A similar intuition operates at the other extreme. The combination of population growth and the harvest policy push the stock towards the optimal escapement level. Therefore, Measurements of stock sizes significantly below that target escapement are most likely to result from measurement error (when present) underestimating the stock size. Consequently, the quota set for small stock sizes is higher (thus the escapement is lower) under the assumption of measurment error than without it.
The intepretation for implementation error is more subtle, and indeed the effect is less extreme. Just as for very high stock sizes the proposed quotas were very large and very uncertain, for low stock sizes the quotas are very small. In fact, for the deterministic case, quotas are zero below constant escapement level. If we add only implementation noise (Figure 2), below the deterministic constant escapement level there is no reason to deviate from the deterministic policy: a quota of zero gets a harvest of exactly zero. As the stock sizes get just above this, harvesting begins at small amounts. If only implementation error is present, it makes sense to harvest slightly more (setting escapement slightly lower) than in the deterministic case because the value (dock-side revenue) from under-harvesting a little bit is less than from over-harvesting due to the discounting and concave nature of the stock-recruitment curve. (Knocking the population down a little too low means you derive the immediate benefit of the higher harvest, at the cost of the stock next year being lower than it would be a the optimal harvest. Fishing too little imposes an immediate cost of less harvest than optimal, and a slight benefit of higher-than-expected stock the following year. However, future benefits are discounted, and growth rate at higher population densities slows).
Implementation error can also interact with other forms of uncertainty to result in a more aggressive
harvest even for stock values measured below the deterministic constant escapement. By comparing the
results of all uncertainties set to "small" noise in Figure 1 to implementation uncertainty being large
(while others are small but non-zero), we see that non-zero quotas begin even before they would under
the same amount of measurement error but smaller implementation error (i.e. the more agressive quota
seen for measured stock sizes below the deterministic escapement is not driven by measurement uncertainty alone.)
At small stock sizes, the possibility of bumper growth (from the stochastic recruitment process) and the possibility of understimated
stock sizes make it worthwhile to set a more agressive quota when implementation is very uncertain than when it is
relatively precise.
sigma_g = 0.1 sigma_m = 0.5 sigma_i = 0.1 f = "logistic" noise_dist = "uniform" delta <- 0.05 grid <- seq(0,250,length=401) S <- multiple_uncertainty(f = f, x_grid = grid, sigma_g = sigma_g, sigma_m = sigma_m, sigma_i = sigma_i, delta = delta, noise_dist = noise_dist)
sims <- function(i) scenario(S, f = f, x_grid = grid, x0= 75, Tmax = 500, sigma_g = sigma_g, sigma_m = sigma_m, sigma_i = sigma_i, noise_dist = noise_dist) sim <- purrr::map_df(1:100, sims) sim %>% mutate(s = x - h) %>% select(t, x, y, s) %>% tidyr::gather(variable, value, -t) %>% ggplot() + geom_density(aes(value, fill=variable), alpha=0.2)
While such arguments help us understand the edge behavior, only the precise optimal policy calculation can deterime exactly where the policies cross -- that is, where increased uncertainties results in a more conservative or more agressive fishing quota. Under uniform noise, we see a stark difference between only implementation uncertainty and only measurement uncertainty in the range of measured stock sizes that is of greatest interest (e.g. most likely to be observed): As measurement uncertainty increases, the policy is consistently more agressive (higher quotas, lower escapement) across this entire range, while as implementation uncertainty increases the policy becomes more conservative (lower quotas, higher escapement).
This is by no means a generic result. Under lognormal noise, the results of measurement vs implementation uncertainty are far more similar. Policy under implementation error is always more conservative than an equivalent magnitude of measurement error, all else equal, but measurement error crosses the deterministic policy far sooner.
If we consider much smaller growth rates, the pattern is dramatically different:
fig3 <- function(noise){ grid <- seq(0, 200, by=0.5) ## Scale noise comparably for log-normal vs uniform if(noise == "lognormal"){ lo <- 0 #0.0577 # 0.1 / sqrt(3) hi <- 0.2887 # 0.5 / sqrt(3) } else { lo <- 0.0 hi <- 0.5 } model <- function(x, h) logistic(x, h, r = 0.1) small <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist) growth <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist) measure <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise_dist) implement <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise_dist) df <- data.frame(y_grid = grid, small = small, growth = growth, measure = measure, implement = implement) %>% tidyr::gather(scenario, value, -y_grid) } df <- data.frame(noise = c("uniform", "lognormal")) %>% dplyr::group_by(noise) %>% dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) + geom_line() + facet_wrap(~ noise) + xlab("Measured Stock Size") + ylab("Proposed Escapement") + coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) + theme_bw()
If we use a compensating instead of overcompensating recruitment function such as Beverton-Holt, the pattern is also quite different.
fig3 <- function(noise){ grid <- seq(0, 200, by=0.5) ## Scale noise comparably for log-normal vs uniform if(noise == "lognormal"){ lo <- 0 #0.0577 # 0.1 / sqrt(3) hi <- 0.2887 # 0.5 / sqrt(3) } else { lo <- 0.0 hi <- 0.5 } model <- "bevertonholt" small <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist) growth <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = hi, sigma_m = lo, sigma_i = lo, noise_dist = noise_dist) measure <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = hi, sigma_i = lo, noise_dist = noise_dist) implement <- multiple_uncertainty(f = model, x_grid = grid, sigma_g = lo, sigma_m = lo, sigma_i = hi, noise_dist = noise_dist) df <- data.frame(y_grid = grid, small = small, growth = growth, measure = measure, implement = implement) %>% tidyr::gather(scenario, value, -y_grid) } df <- data.frame(noise = c("uniform", "lognormal")) %>% dplyr::group_by(noise) %>% dplyr::do(fig3(.$noise))
df %>% ggplot(aes(x = y_grid, y = value, col = scenario)) + geom_line() + facet_wrap(~ noise) + xlab("Measured Stock Size") + ylab("Proposed Escapement") + coord_cartesian(xlim = c(0, 150), ylim = c(0,100)) + theme_bw()
The magnitude of the noise makes some difference as well, with larger magnitudes driving earlier crossing points.
Discussion
Review of general lit
Consistent with Sethi et al, we observe:
- Small uncertainties lead to only small deviations from classical expectations based on the deterministic model or growth-noise only model.
- Large measurement uncertainty can lead to a significant deviation from "constant escapement"
In contrast to Sethi et al, we find that:
-
These observations are sensitive to way uncertainty is modelled, rather than being generally robust as claimed.
-
Large implementation uncertainty can also result in substantial deviations from the the "constant escapement" rule.
This is most pronounced in more realisitic models of uncertainty, such as log-normally distributed noise, but visible in uniform noise as well. Small deviations from constant escapement are visible in the implementation error results of @Sethi2005, but are understated. This contrasts to the @Sethi claim that
... implementation sources of uncertainty is high, a constant-escapement policy is qualitatively appropriate since the slopes of these policies beyond the shutdown point are virtually flat.
Substantial deviations from the optimal policy can be found for non-convex growth functions. In particular, when there is small or negligible cost to harvest, the optimal quota can signficantly exceed the measured stock size, as the true stock size and true harvest can differ substantially.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.