# The effects of random and seasonal environmental fluctuations on optimal
harvesting and stocking

Alexandru Hening  Department of Mathematics
Texas A&M University
Mailstop 3368
College Station, TX 77843-3368
United States<EMAIL_ADDRESS>,  Ky Quan Tran  Department of Mathematics and
Statistics
The State University of New York in Korea
119 SongdoMoonhwa-Ro, Yeonsu-Gu
Yeonsu-Gu, Incheon 21985
Korea<EMAIL_ADDRESS>and  Sergiu C. Ungureanu  Department of
Economics
City, University of London
Northampton Square
London EC1V 0HB
United Kingdom<EMAIL_ADDRESS>
###### Abstract.

We analyze the harvesting and stocking of a population that is affected by
random and seasonal environmental fluctuations. The main novelty comes from
having three layers of environmental fluctuations. The first layer is due to
the environment switching at random times between different environmental
states. This is similar to having sudden environmental changes or
catastrophes. The second layer is due to seasonal variation, where there is a
significant change in the dynamics between seasons. Finally, the third layer
is due to the constant presence of environmental stochasticity—between the
seasonal or random regime switches, the species is affected by fluctuations
which can be modelled by white noise. This framework is more realistic because
it can capture both significant random and deterministic environmental shifts
as well as small and frequent fluctuations in abiotic factors. Our framework
also allows for the price or cost of harvesting to change deterministically
and stochastically, something that is more realistic from an economic point of
view.

The combined effects of seasonal and random fluctuations make it impossible to
find the optimal harvesting-stocking strategy analytically. We get around this
roadblock by developing rigorous numerical approximations and proving that
they converge to the optimal harvesting-stocking strategy. We apply our
methods to multiple population models and explore how prices, or costs, and
environmental fluctuations influence the optimal harvesting-stocking strategy.
We show that in many situations the optimal way of harvesting and stocking is
not of threshold type.

###### Key words and phrases:

Harvesting; stochastic environment; density-dependent price; controlled
diffusion; switching environment; seasonality

###### 2010 Mathematics Subject Classification:

92D25, 60J70, 60J60

###### Contents

  1. 1 Introduction
  2. 2 Model and Results
    1. 2.1 Numerical Scheme
  3. 3 Extensions
    1. 3.1 Variable effort harvesting-stocking strategies
      1. 3.1.1 Constructing the approximation scheme.
    2. 3.2 Uncertain price functions
      1. 3.2.1 Constructing the approximation scheme.
    3. 3.3 The combined effects of seasonality and Markovian switching
      1. 3.3.1 Constructing the approximation scheme.
  4. 4 Numerical Examples
    1. 4.1 Introducing switching in a Verhulst-Pearl system
    2. 4.2 Analysis of the effect of the control cost
    3. 4.3 Large and small transition rates with switching
    4. 4.4 Switching in the cost and price functions
    5. 4.5 Interlude with intuition and conjecture
    6. 4.6 Variable effort harvesting examples
    7. 4.7 Stochastic pricing examples
    8. 4.8 Seasonal and periodic variability
  5. 5 Discussion
  6. A Transition Probabilities
    1. A.1 The formulation from Section 2
    2. A.2 Variable effort harvesting-stocking strategies
    3. A.3 Uncertain price functions
    4. A.4 The combined effects of seasonality and Markovian switching
  7. B Continuous–time interpolation
  8. C Convergence
  9. D Numerical Experiments
    1. D.1 Varying the cost dependency
    2. D.2 The Gompertz model of population growth
    3. D.3 The Nisbet-Gurney model of population growth

##  1\. Introduction

A fundamental problem in conservation ecology is finding the optimal strategy
for harvesting a species. This problem is important because excessive
harvesting can drive species extinct while under-harvesting ensures the loss
of valuable resources. If one assumes that the dynamics and harvesting happen
in continuous time, there has been significant progress in finding the optimal
harvesting strategies which maximize the total discounted or asymptotic
harvest yield—see the work by Abakuks & Prajneshu (1981); Lungu & Øksendal
(1997); Alvarez E. & Shepp (1998); Hening, Nguyen, Ungureanu & Wong (2019);
Alvarez E & Hening (2020); Cohen et al. (2021). These studies have shown that
in a very general setting the optimal strategy is of threshold or bang-bang
type: there exists a threshold $w>0$ such that whenever the population is
under the threshold there is no harvesting while when the population is above
the threshold one harvests at the maximal (possibly infinite) rate. In this
paper we look at whether this result is true in more general and realistic
models.

The mathematical framework we will use is the one of stochastic differential
equations with switching (SSDE). A SSDE has a discrete component that keeps
track of the environment, and which changes at random times. In a fixed
environmental state the system is modelled by a stochastic differential
equation. This way one can capture the more realistic behaviour of two types
of environmental fluctuations:

  * •

major environmental shifts (daily or seasonal changes, catastrophes),

  * •

abiotic fluctuations within each environment.

We focus on the most natural setting, when the switching rates are constant.
In this case, the process spends an exponentially distributed time in each
environmental state and then switches to a different state. Environmental
switches have been shown to fundamentally alter the fate of ecological
communities by reversing competitive exclusion into coexistence or having
other unexpected results (Hening & Strickler, 2019; Hening & Nguyen, 2020;
Benaïm & Lobry, 2016; Bourquin, 2021; Hening et al., 2021). We explore how
switching impacts the harvesting and stocking of a species.

Even though the general properties for SSDE have been studied thoroughly (Yin
& Zhu, 2009; Zhu & Yin, 2009; Nguyen et al., 2017), there are few results
regarding the persistence or harvesting of ecological systems modelled by SSDE
(Song et al., 2011; Tran & Yin, 2015, 2017; Bao & Shao, 2016; Song & Zhu,
2016; Hening & Li, 2020). We fill this gap by providing an analysis
complemented by an in depth look at some specific illuminating examples. In
particular, we look at the logistic equation which has been used extensively
in fisheries and other harvesting settings (Clark, 2010; Alvarez E. & Shepp,
1998).

The SSDE framework is generalized even further by including deterministic
seasonal variation, which will make the various coefficients depend explicitly
on time in a periodic fashion. There are few studies which look at the
interaction of harvesting and seasonal variability. Some focus on very
specific models or look only at the purely deterministic setting (Cromer,
1988; Fan & Wang, 1998; Brauer & Sànchez, 2003; Xu et al., 2005; Braverman &
Mamdani, 2008; Bohner & Streipert, 2016). The current paper provides important
generalizations to these previous results as we can analyze a very wide range
of models.

In addition to modelling the ecological dynamics, one also has to have a
robust way of modelling the economics. The price of the harvested species can
depend on the population size, the state of the environment and also
explicitly on time. Our framework also includes a realistic cost that is
incurred through stocking or harvesting. This cost is due to fees associated
with harvesting-stocking policies, state constraints, certain taxes, or
incentives that the manager must follow or can receive.

The price and cost functions depend on time both directly, and indirectly
through the population state and the environment switching state.
Economically, the indirect time dependence is important since it captures the
growth-delayed-return trade-off, which is a typical feature of harvesting
models. The direct time dependence adds another important layer of realism,
since it can accommodate exogenous price changes. These can come from varying
commodity prices (external supply and demand shocks), and varying costs of
inputs of production.
The harvested quantity can indirectly change the price in the market by
affecting the total supply. Because our price and cost functions depend
directly on the harvested quantities, we can accommodate for this endogenous
effect on the total supply. The numerical methods can simulate the price
changes with a suitable choice of price and cost functions. Because the market
side is not modeled explicitly, this can only be done ad hoc, for each
particular application. A discussion of the importance of market variables in
fisheries management can be found in Sylvia (1994); Pooley (1987); Asche et
al. (2015).

The main novelties of our work are the following:

  1. (1)

We formulate the harvesting-stocking problem for a species that is influenced
by three types of environmental fluctuations: the first due to major regime
shifts, the second due to constant abiotic changes, and the third due to
seasonality.

  2. (2)

The price for harvesting or stocking is realistic and depends on the
population size as well as the state of the environment. Furthermore, stocking
and harvesting incur a cost as well.

  3. (3)

We develop rigorous numerical approximation schemes based on the Markov chain
approximation method.

  4. (4)

We discover interesting new phenomena by analyzing in depth some important
examples.

The rest of the paper is organized as follows. In Section 2 we describe our
model and the main results. In Section 3 we discuss several extensions of the
proposed model. Particular examples are explored using the newly developed
numerical schemes in Section 4. The discussion of our results is in Section 5.
Finally, all the technical proofs appear in the appendices.

##  2\. Model and Results

Assume we have a probability space $(\Omega,\mathcal{F},{\mathbb{P}})$
satisfying the usual conditions. Denote by $X(t)$ the size of the population
at time $t\geq 0$. For a natural population, without harvesting or stocking,
the dynamics is given by

(2.1)
$dX(t)=b\big{(}X(t),\alpha(t)\big{)}dt+\sigma\big{(}X(t),\alpha(t)\big{)}dw(t),$

where $w(\cdot)$ is a standard Brownian motion, $\alpha(t)$ is an irreducible
continuous-time Markov chain taking values in
$\mathcal{M}=\\{1,\dots,m_{0}\\}$, and
$b,\sigma:\mathbb{R}_{+}\times{\mathcal{M}}\to{\mathbb{R}}$ are smooth enough
functions. Furthermore, we assume that the Brownian motion $w(\cdot)$ and the
Markov chain $\alpha(\cdot)$ are independent and denote by $q_{ij}$ the
transition rates of $\alpha(t)$. This means that

(2.2)
$\begin{array}[]{ll}&\displaystyle{\mathbb{P}}\\{\alpha(t+\Delta)=j~{}|~{}\alpha(t)=i,X(s),\alpha(s),s\leq
t\\}=q_{ij}\Delta+o(\Delta),\text{ if }i\neq j.\\\
&\displaystyle{\mathbb{P}}\\{\alpha(t+\Delta)=i~{}|~{}\alpha(t)=i,X(s),\alpha(s),s\leq
t\\}=1+q_{ii}\Delta+o(\Delta),\end{array}$

where $q_{ii}:=-\sum_{j\neq i}q_{ij}$. As usual, we assume that
$b(0,\alpha)=\sigma(0,\alpha)=0$ since if the population goes extinct, it
should not be able to get resurrected without external intervention (like a
repopulation/stocking event). If $X(t_{0})=0$ for some $t_{0}\geq 0$, then
$X(t)=0$ for any $t\geq t_{0}$. Thus, $X(t)\in{\mathbb{R}}_{+}$ for any $t\geq
0$.

Biological interpretation: If at time $t$ the state of the environment is
$\alpha(t)=k$ the dynamics is governed by the SDE

$dX(t)=b\big{(}X(t),k\big{)}dt+\sigma\big{(}X(t),k\big{)}dw(t)$

with nonlinear drift term $\mu(x,k)$ and diffusion term $\sigma(x,k)$.
However, the environment can switch to a different state $j$ with the
transition rate $q_{kj}$. In the small time $\Delta$ the probability of
transition to state $j$ is approximately $q_{kj}\Delta$.

Note that since the transition matrix $Q=(q_{ij})_{m_{0}\times m_{0}}$ is
independent of the population of the species, the jump times of $\alpha(t)$
will be exponentially distributed and independent of the process $X$.

Let $U(t)$ be the harvesting-stocking rate of the population. This means that
in the small time-interval $(t,t+\Delta t)$ one harvests $U(t)\Delta t$ if
$U(t)>0$ and one stocks $U(t)\Delta t$ if $U(t)<0$.

The dynamics of the population $X(t)$ that includes harvesting and stocking
becomes

(2.3)
$X(t)=x+\int\limits_{0}^{t}\big{(}b(X(s),\alpha(s))-U(s)\big{)}ds+\int\limits_{0}^{t}\sigma\big{(}X(s),\alpha(s)\big{)}dw(s).$

In order to have a well defined system we also need to impose some initial
conditions:

(2.4) $X(0)=x\in{\mathbb{R}}_{+},\quad\alpha(0)=\alpha\in{\mathcal{M}}.$

Note that $\max\\{U(t),0\\}$ and $\max\\{-U(t),0\\}$ are the magnitudes of the
harvesting and the stocking rate, respectively. We suppose that $U(t)$ takes
values in a compact subset $\mathcal{U}\subset{\mathbb{R}}$ and
$0\in\mathcal{U}$. This means that the maximal harvesting and stocking rates
are bounded and that one can always choose to not interfere with the
population, that is, pick $U(t)=0$ at certain times $t\geq 0$.

Let $\mathcal{A}_{x,\alpha}$ denote the collection of all admissible controls
with initial value $(x,\alpha)\in{\mathbb{R}}_{+}\times{\mathcal{M}}$. The
collection $\mathcal{A}_{x,\alpha}$ of admissible controls is defined to
contain the elements $U(\cdot)$ such that

  * •

$U(t)$ is $\mathcal{F}(t)$-adapted,

  * •

$U(t)\in\mathcal{U}$ for any $t\geq 0$,

  * •

$X(t)\geq 0$ for any $t\geq 0$.

Let $P:\mathbb{R}_{+}\times{\mathcal{M}}\to[0,\infty)$ be the price of the
species per unit. The price depends on the population size and the
environmental state. The accumulated income can be written as

$\int_{0}^{\infty}e^{-\delta s}P\big{(}X(s),\alpha(s)\big{)}\cdot U(s)ds,$

where $\delta>0$ is the time discount rate. In addition, we also suppose that
the harvesting-stocking is costly, something that is described by the cost
function
$C:{\mathbb{R}}_{+}\times{\mathcal{M}}\times\mathcal{U}\to[0,\infty)$. Thus,
for a control strategy $U(\cdot)\in\mathcal{A}_{x,\alpha}$, we define the
performance function as

(2.5)
$J\big{(}x,\alpha,U(\cdot)\big{)}:={\mathbb{E}}_{x,\alpha}\int_{0}^{\infty}e^{-\delta
s}\Big{[}P(X(s),\alpha(s))\cdot
U(s)-C\big{(}X(s),\alpha(s),U(s)\big{)}\Big{]}ds,$

where ${\mathbb{E}}_{x,\alpha}$ denotes the expectation with respect to the
probability law when the process $(X(t),\alpha(t))$ starts with initial
condition $(x,\alpha)$. The goal is to maximize $J(x,\alpha,U(\cdot))$ and
find an optimal strategy $U^{*}(\cdot)$ such that

(2.6)
$J\big{(}x,\alpha,U^{*}(\cdot)\big{)}=V(x,\alpha):=\sup\limits_{U(\cdot)\in\mathcal{A}_{x,\alpha}}J\big{(}x,\alpha,U(\cdot)\big{)}.$

For notational simplicity, we collect the price function $P(\cdot)$ and the
cost function $C(\cdot)$ into the price-cost function
$p:{\mathbb{R}}_{+}\times{\mathcal{M}}\times\mathcal{U}\to{\mathbb{R}}$ given
by

$p(x,\alpha,u)=P(x,\alpha)\cdot u-C\big{(}x,\alpha,u\big{)}.$

As a result, the performance function can be written as

(2.7) $J(x,\alpha,U(\cdot))={\mathbb{E}}_{x,\alpha}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X(s),\alpha(s),U(s)\big{)}ds.$

The following standing assumptions are made throughout the paper.

######  Assumption 2.1.

  * (a)

The functions $b(\cdot,\alpha)$ and $\sigma(\cdot,\alpha)$ are locally
Lipschitz continuous for each $\alpha\in{\mathcal{M}}$. Moreover, for any
initial condition $(x,\alpha)\in{\mathbb{R}}_{+}\times{\mathcal{M}}$, the
uncontrolled system (2.1) has a unique global solution.

  * (b)

The function $p(\cdot,\alpha,\cdot)$ is bounded and continuous for each
$\alpha\in{\mathcal{M}}$.

######  Remark 2.2.

We recall that in recent papers by Hening, Tran, Phan & Yin (2019); Hening &
Tran (2020) the authors assume that the stocking price is higher than the
harvesting price. In the current paper, however, the stocking and harvesting
have the same price which is a function of the current population size and the
environmental regime. We take into account the control cost by introducing the
price-cost function $p(\cdot)$. As a result our setting is more general than
those considered by Hening, Tran, Phan & Yin (2019); Hening & Tran (2020). The
cost function $C(\cdot)$ allows us to take into account taxes, state
constraints, and incentives for sustainability that arise in practical
settings. We will discuss several extensions and related formulations in
Section 3. In the first extension, we look at the setting where harvesting in
a short time $dt$ looks like $U(t)X(t)dt$, and $U(t)$ lives in a compact set
including $0$—note that the harvesting rate in time can be unbounded in this
setting. In the second extension, we look at the setting where the price is
time-dependent, uncertain and its evolution is described by a stochastic
differential equation. In the third extension, we consider the combined
effects of random and seasonal (periodic) environmental fluctuations.

A common approach in treating stochastic control problems is to characterize
the value function as a viscosity solution of a Hamilton-Jacobi-Bellman
equation. Formally, the associated equation of the underlying problem is given
by

(2.8)
$\begin{array}[]{ll}\max\limits_{u\in\mathcal{U}}\left[V^{\prime}(x,\alpha)\big{(}b(x,\alpha)-u\big{)}+\dfrac{1}{2}\sigma^{2}(x,\alpha)V^{\prime\prime}(x,\alpha)\right.\\\
\left.+\sum\limits_{k=1}^{m_{0}}q_{\alpha k}V(x,k)+p(x,\alpha,u)-\delta
V(x,\alpha)\right]=0,\end{array}$

for $(x,\alpha)\in\mathbb{R}_{+}\times{\mathcal{M}}$. One can try to adopt
this approach by following the techniques used by Zhu (2011). However, as
pointed out by Zhu (2011), under regime switching the associated Hamilton-
Jacobi-Bellman equation is complicated and a closed-form solution is virtually
impossible to obtain. Therefore, in order to treat the underlying formulation
and extensions in Section 3, we will approximate the original problem by a
sequence of discrete-time optimal control problems of approximating Markov
chains.

######  Remark 2.3.

Although our work is motivated by the presence of a Markovian switching
environment, the numerical schemes we developed can also be used to solve
harvesting problems for diffusion models without switching. In particular, one
can simply take ${\mathcal{M}}=\\{1\\}$ and put $q_{11}=0$.

###  2.1. Numerical Scheme

We provide a numerical approach to gain information about the value function
and the optimal harvesting-stocking strategy. Using the Markov chain
approximation method developed by Kushner (1990) and Kushner & Dupuis (1992),
we construct a controlled Markov chain in discrete time that approximates the
controlled stochastic processes. As pointed out by Kushner & Dupuis (1992) a
probabilistic approach using the Markov chain approximation method for
controlled diffusions has some important advantages. First, the Markov chain
approximation method allows one to use physical insights derived from the
dynamics of the controlled diffusion in obtaining a suitable approximation
scheme. Second, the Markov chain approximation method does not require
significant regularity of the controlled processes (2.3) nor does it rely on
the uniqueness properties of the associated HJB equation (2.8).

Let $h>0$ be a discretization parameter. Define
$S_{h}:=\\{kh:k\in\mathbb{Z}_{\geq 0}\\}$ and let
$\\{(X^{h}_{n},\alpha^{h}_{n}):n\in\mathbb{Z}_{\geq 0}\\}$ be a discrete-time
controlled Markov chain with state space $S_{h}\times{\mathcal{M}}$. At each
discrete-time step $n$, the magnitude of the harvesting-stocking component
$U^{h}_{n}$ must be specified. Let $U^{h}=\\{U^{h}_{n}\\}$ be a sequence of
controls. We denote by $q^{h}\left((x,k),(y,l)|u\right)$ the transition
probability from state $(x,k)$ to another state $(y,l)$ under the control $u$.
Denote $\mathcal{F}^{h}_{n}=\sigma\\{X^{h}_{m},\alpha^{h}_{m},U^{h}_{m},m\leq
n\\}$.

The sequence $U^{h}$ is said to be admissible if it satisfies the following
conditions:

  * (a)

$U^{h}_{n}$ is
$\sigma\\{X^{h}_{0},\dots,X^{h}_{n},\alpha^{h}_{0},\dots,\alpha^{h}_{n},U^{h}_{0},\dots,U^{h}_{n-1}\\}-\text{adapted}$
and $U^{h}_{n}\in\mathcal{U}$ for each $n$;

  * (b)

For any $(x,k)\in S_{h}\times{\mathcal{M}}$, we have

$\begin{split}{\mathbb{P}}\left\\{\left(X^{h}_{n+1},\alpha^{h}_{n+1}\right)=(x,k)|\mathcal{F}^{h}_{n}\right\\}&={\mathbb{P}}\left\\{\left(X^{h}_{n+1},\alpha^{h}_{n+1}\right)=(x,k)|X^{h}_{n},\alpha^{h}_{n},U^{h}_{n}\right\\}\\\
&=q^{h}\left((X^{h}_{n},\alpha^{h}_{n}),(x,k)|U^{h}_{n}\right);\end{split}$

  * (c)

$X^{h}_{n}\in S_{h}$ for all $n\in\mathbb{Z}_{\geq 0}$.

The class of all admissible control sequences $U^{h}$ for initial state
$(x,\alpha)$ will be denoted by $\mathcal{A}^{h}_{x,\alpha}$.

For each $(x,\alpha,u)\in S_{h}\times{\mathcal{M}}\times\mathcal{U}$, we
define a family of interpolation intervals $\Delta t^{h}(x,\alpha,u)$. The
values of $\Delta t^{h}(x,\alpha,u)$ will be specified later. Then we define

$t^{h}_{0}=0,\quad\Delta t^{h}_{m}=\Delta
t^{h}(X^{h}_{m},\alpha^{h}_{m},U^{h}_{m}),\quad
t^{h}_{n}=\sum\limits_{m=0}^{n-1}\Delta t^{h}_{m}.$

For $(x,\alpha)\in S_{h}\times{\mathcal{M}}$ and
$U^{h}\in\mathcal{A}^{h}_{x,\alpha}$, the performance function and the value
function for the controlled Markov chain is defined as

(2.9) $J^{h}(x,\alpha,U^{h})={\mathbb{E}}\sum_{m=0}^{\infty}e^{-\delta
t_{m}^{h}}p(X^{h}_{m},\alpha^{h}_{m},U^{h}_{m})\Delta t_{m}^{h},\quad
V^{h}(x,\alpha)=\sup\limits_{U^{h}\in\mathcal{A}^{h}_{x,\alpha}}J^{h}(x,\alpha,U^{h}).$

The corresponding dynamic programming equation for the discrete approximation
is given by

$\begin{split}V^{h}(x,\alpha)=\max_{u\in\mathcal{U}}\bigg{[}e^{-\delta\Delta
t^{h}(x,\alpha,u)}\sum\limits_{(y,\beta)\in
S_{h}\times{\mathcal{M}}}V^{h}(y,\beta)q^{h}\big{(}(x,\alpha),(y,\beta)|u\big{)}+p(x,\alpha,u)\Delta
t^{h}(x,\alpha,u)\bigg{]}.\end{split}$

We will construct the transition probabilities and interpolation intervals
such that the Markov chain $\left\\{(X^{h}_{n},\alpha^{h}_{n})\right\\}$
approximates the process $\left\\{(X(\cdot),\alpha(\cdot)\right\\}$ well, in
the sense that they are locally consistent. Then the similarity between (2.7)
and (2.9) suggests that the values $V^{h}(x,\alpha)$ and $V(x,\alpha)$ will be
close for small $h$, and this will turn out to be the case. Solving the
optimal harvesting problem for the chain
$\left\\{(X^{h}_{n},\alpha^{h}_{n})\right\\}$, we obtain an approximating
optimal value and an approximating optimal strategy for the continuous-time
process $\big{(}X(\cdot),\alpha(\cdot)\big{)}$. The main convergence result is
given below.

######  Theorem 2.4.

Suppose Assumptions 2.1 holds. Then for any
$(x,\alpha)\in{\mathbb{R}}_{+}\times{\mathcal{M}}$, $V^{h}(x,\alpha)\to
V(x,\alpha)$ as $h\to 0$. Thus, for sufficiently small $h$, a near-optimal
harvesting-stocking strategy of the controlled Markov chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ is also a near-optimal harvesting-stocking
policy of $\big{(}X(\cdot),\alpha(\cdot)\big{)}$.

######  Remark 2.5.

From a control theoretic point of view, an admissible control
$U(\cdot)\in\mathcal{A}_{x,\alpha}$ is called $\varepsilon$-optimal with
respect to $(x,\alpha)$ if $J(x,\alpha,U(\cdot))\geq V(x,\alpha)-\varepsilon$.
A family of admissible controls
$U_{\varepsilon}(\cdot)\in\mathcal{A}_{x,\alpha}$ parameterized by
$\varepsilon>0$ is called near-optimal with respect to $(x,\alpha)$ if
$J(x,\alpha,U_{\varepsilon}(\cdot))\geq V(x,\alpha)-\psi(\varepsilon)$ for
sufficiently small $\varepsilon$, where $\psi(\cdot)$ is a function of
$\varepsilon$ satisfying $\psi(\varepsilon)\to 0$ as $\varepsilon\to 0$. The
$\varepsilon$-optimal controls and near-optimal controls for the chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ can be defined similarly. By virtue of
Theorem 2.4 (see also Theorem C.5), as $h$ approaches $0$, a near-optimal
harvesting-stocking strategy of the controlled Markov chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ will provide a near-optimal harvesting-
strategy for$\big{(}X(\cdot),\alpha(\cdot)\big{)}$. As a result, the original
problem can be reduced to solving for a near-optimal harvesting-stocking
policy of the Markov chain $\\{(X^{h}_{n},\alpha^{h}_{n})\\}$.

##  3\. Extensions

###  3.1. Variable effort harvesting-stocking strategies

In this subsection, we analyze an extension of our framework for which the
harvesting-stocking rates are not forced to be bounded as in (2.3)—see Alvarez
& Shepp (1998, Section 3) and work by I. Kharroubi (2019). This allows us to
take into account not only the impact of the harvesting-stocking effort, but
also the abundance of the species itself. It is natural to assume that the
harvesting or stocking rate is proportional to the population size (Hening,
Nguyen, Ungureanu & Wong, 2019). These types of harvesting strategies, where a
constant fraction of the population is harvested, are called constant effort
harvesting strategies and are widely used in modeling fisheries. We suppose
that $\mathcal{U}$ is a compact subset of $\mathbb{R}$ and that
$0\in\mathcal{U}$. Our harvesting-stocking strategy will be such that in the
small time $dt$ we harvest or stock the amount $U(t)X(t)\,dt$ where
$U(t)\in\mathcal{U}$. We call these strategies variable effort harvesting-
stocking strategies. The dynamics that includes stocking and harvesting is
given by

(3.1)
$X(t)=x+\int\limits_{0}^{t}\big{(}b(X(s),\alpha(s))-U(s)X(s)\big{)}ds+\int\limits_{0}^{t}\sigma\big{(}X(s),\alpha(s)\big{)}dw(s),\quad\alpha(0)=\alpha\in{\mathcal{M}}.$

Let $\mathcal{A}_{x,\alpha}$ denote the collection of all admissible controls
with initial value $(x,\alpha)\in{\mathbb{R}}_{+}\times{\mathcal{M}}$, where a
strategy $U(\cdot)$ will be in $\mathcal{A}_{x,\alpha}$ if $U(t)$ is
$\mathcal{F}(t)$-adapted, $U(t)\in\mathcal{U}$ for any $t\geq 0$, and
$X(t)\geq 0$ for any $t\geq 0$. The corresponding performance function and
value function are given by

$\displaystyle p(x,\alpha,u)$ $\displaystyle=P(x,\alpha)\cdot x\cdot
u-C\big{(}x,\alpha,u\big{)},$ $\displaystyle\overline{J}(x,\alpha,U(\cdot)):$
$\displaystyle={\mathbb{E}}_{x,\alpha}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X(s),\alpha(s),U(s)\big{)}ds,$ $\displaystyle\overline{V}(x,\alpha)$
$\displaystyle=\sup\limits_{U(\cdot)\in\mathcal{A}_{x,\alpha}}\overline{J}\big{(}x,\alpha,U(\cdot)\big{)}.$

Compared with the stocking-harvesting model (2.3), in this setting, one can
harvest with higher rates when $X(t)$ is large, while we can only stock at
lower rates when $X(t)$ is small. It is also clear that the stocking rate when
$X(t)=0$ is $U(t)X(t)=0$ so that $V(0,\alpha)=0$ for any
$\alpha\in{\mathcal{M}}$. This implies that one cannot restock a species once
it is extinct.

####  3.1.1. Constructing the approximation scheme.

We use the same method as in the preceding section to construct a discrete-
time controlled Markov chain
$\\{(X^{h}_{n},\alpha^{h}_{n}):n\in\mathbb{Z}_{\geq 0}\\}$ with state space
$S_{h}\times{\mathcal{M}}=\\{kh:k\in\mathbb{Z}_{\geq
0}\\}\times{\mathcal{M}}$. At each step $n$, we need to specify a control
$U^{h}_{n}\in\mathcal{U}$. We need to define transition probabilities
$q^{h}((x,k),(y,l)|u)$ and interpolation intervals $\Delta t^{h}(x,\alpha,u)$
so that the controlled Markov chain $\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ is
locally consistent with respect to the controlled diffusion (3.1). For
$(x,\alpha,u)\in S_{h}\times{\mathcal{M}}\times\mathcal{U}$ and the family of
the interpolation intervals $\Delta t^{h}(x,\alpha,u)$, we define

$t^{h}_{0}=0,\quad\Delta t^{h}_{m}=\Delta
t^{h}(X^{h}_{m},\alpha^{h}_{m},U^{h}_{m}),\quad
t^{h}_{n}=\sum\limits_{m=0}^{n-1}\Delta t^{h}_{m}.$

For $(x,\alpha)\in S_{h}\times{\mathcal{M}}$, the class
$\mathcal{A}^{h}_{x,\alpha}$ of all admissible control sequences $U^{h}$ for
initial state $(x,\alpha)$ can be defined as before. Then the performance
function for the controlled Markov chain is defined as

$\overline{J}^{h}(x,\alpha,U^{h})={\mathbb{E}}\sum_{m=0}^{\infty}e^{-\delta
t_{m}^{h}}p(X^{h}_{m},\alpha^{h}_{m},U^{h}_{m})\Delta t_{m}^{h}.$

The value function of the controlled Markov chain is

$\overline{V}^{h}(x,\alpha)=\sup\limits_{U^{h}\in\mathcal{A}^{h}_{x,\alpha}}\overline{J}^{h}(x,\alpha,U^{h}).$

The convergence result given in Theorem 2.4 is also valid in this case, if
$\overline{V}^{h}(x,\alpha)$ and $\overline{V}(x,\alpha)$ are substituted for
$V^{h}(x,\alpha)$ and $V(x,\alpha)$, respectively.

###  3.2. Uncertain price functions

Since the harvesting-stocking problem we investigate is on an infinite time
horizon, it is natural to look at prices that can change in time. In this
section we consider a more realistic scenario with time dependent price
functions (Alvarez & Koskela, 2007; I. Kharroubi, 2019). Working with the
population model given by (2.3), we suppose that the per-unit price-cost
function is $p_{0}(x,\alpha,u)+\Phi(t)$ for
$(t,x,\alpha,u)\in{\mathbb{R}}^{2}_{+}\times{\mathcal{M}}\times\mathcal{U},$
where $p_{0}(\cdot)$ is a deterministic function and $\Phi(\cdot)$ satisfies
the following stochastic differential equation

(3.2)
$d\Phi(t)=b_{0}\big{(}\Phi(t),\alpha(t)\big{)}dt+\sigma_{0}\big{(}\Phi(t),\alpha(t)\big{)}dw_{0}(t),\quad\Phi(0)=\phi\in\mathbb{R}_{+}.$

We assume that $w_{0}(\cdot)$ is a standard Brownian motion in $\mathbb{R}$
which is independent of the Brownian motion $w(\cdot)$ that drives the
dynamics of the species. We also suppose that $b_{0}(\cdot,\alpha)$ and
$\sigma_{0}(\cdot,\alpha)$ are locally Lipschitz continuous function for each
$\alpha\in{\mathcal{M}}$ and are such that (3.2) has a unique solution
$\Phi(\cdot)$, with $\Phi(t)\geq 0,t\geq 0$ and $\sup\limits_{t\geq
0}{\mathbb{E}}|\Phi(t)|<\infty$. We note that the price depends on the current
population size, the environmental regime, time (through (3.2)), and is also
influenced by randomness. Note that if $\sigma_{0}(x,\alpha)\equiv 0$,
$\Phi(\cdot)$ is simply time dependent and no longer uncertain. In this
setting the performance and the value functions are given by

(3.3)
$\begin{split}\widehat{J}(x,\alpha,U(\cdot)):&={\mathbb{E}}_{x,\alpha}\int_{0}^{\infty}e^{-\delta
s}\big{(}p_{0}(X(s),\alpha(s),U(s))+\Phi(s)\big{)}U(s)ds,\\\
\widehat{V}(x,\alpha)&=\sup\limits_{U(\cdot)\in\mathcal{A}_{x,\alpha}}\widehat{J}\big{(}x,\alpha,U(\cdot)\big{)}.\end{split}$

####  3.2.1. Constructing the approximation scheme.

In order to treat the uncertain price function given by (3.2), we combine
$\big{(}X(\cdot),\alpha(\cdot)\big{)}$ and $\Phi(\cdot)$ into one controlled
process $\big{(}\Phi(\cdot),X(\cdot),\alpha(\cdot)\big{)}$ with initial
condition $(\Phi(0),x,\alpha)=(\phi,x,\alpha)$. The set of admissible
strategies $\mathcal{A}_{x,\alpha}$ is defined as in Section 2 and does not
depend on $\Phi(0)$. With the price function above, the performance function
becomes

$\widehat{J}\big{(}\phi,x,\alpha,U(\cdot)\big{)}:={\mathbb{E}}_{\phi,x,\alpha}\int_{0}^{\infty}e^{-\delta
s}\big{(}p_{0}(X(s),\alpha(s),U(s))+\Phi(s)\big{)}U(s)ds,$

where ${\mathbb{E}}_{\phi,x,\alpha}$ denotes the expectation with respect to
the probability law when the process $(\Phi(t),X(t),\alpha(t))$ starts with
initial condition $(\phi,x,\alpha)$ and the value function is

$\widehat{V}(\phi,x,\alpha)=\sup\limits_{U(\cdot)\in\mathcal{A}_{x,\alpha}}\widehat{J}\big{(}\phi,x,\alpha,U(\cdot)\big{)}.$

To approximate the controlled process
$\big{(}\Phi(\cdot),X(\cdot),\alpha(\cdot)\big{)}$, we construct a discrete-
time controlled Markov chain
$\left\\{(\Phi^{h}_{n},X^{h}_{n},\alpha^{h}_{n}):n\in\mathbb{Z}_{\geq
0}\right\\}$ with state space $\widehat{S}_{h}\times{\mathcal{M}}$, where

$\widehat{S}_{h}:=\left\\{(k_{1}h,k_{2}h)^{\prime}\in\mathbb{R}^{2}:k_{i}\in\mathbb{Z}_{\geq
0}\right\\}.$

Note that we add $\\{\Phi^{h}_{n}\\}$ to approximate the evolution of the
price component $\Phi(\cdot)$. Therefore, the computations are more involved
than those in the preceding settings. At each step $n$, we need to specify a
control $U^{h}_{n}\in\mathcal{U}$. We need to define transition probabilities
$q^{h}\big{(}(\phi,x,k),(\psi,y,l)|u\big{)}$ and interpolation intervals
$\Delta t^{h}(\phi,x,k,u)$ so that the controlled Markov chain
$\left\\{(\Phi^{h}_{n},X^{h}_{n},\alpha^{h}_{n})\right\\}$ is locally
consistent with respect to the controlled diffusion (2.3) and (3.2).

For $(\phi,x,\alpha)\in\widehat{S}_{h}\times{\mathcal{M}}$, the class
$\mathcal{A}^{h}_{x,\alpha}$ of all admissible control sequences $U^{h}$ for
initial state $(\phi,x,\alpha)$ can be defined as before. Then the performance
function for the controlled Markov chain is defined as

$\widehat{J}^{h}(\phi,x,\alpha,U^{h})={\mathbb{E}}_{\phi,x,\alpha}\sum_{m=0}^{\infty}e^{-\delta
t_{m}^{h}}\big{(}p_{0}(X^{h}_{m},\alpha^{h}_{m},U^{h}_{m})+\Phi^{h}_{n}\big{)}U^{h}_{m}\Delta
t_{m}^{h}.$

The value function of the controlled Markov chain is

$\widehat{V}^{h}(\phi,x,\alpha)=\sup\limits_{U^{h}\in\mathcal{A}_{x,\alpha}}\widehat{J}^{h}(\phi,x,\alpha,U^{h}).$

As a natural analogue to Theorem 2.4 we get the following convergence result.

######  Theorem 3.1.

Suppose Assumption 2.1 holds, $b_{0}(\cdot,\alpha)$ and
$\sigma_{0}(\cdot,\alpha)$ are locally Lipschitz continuous function for any
fixed $\alpha\in{\mathcal{M}}$, and that (3.2) has a unique solution
$\Phi(\cdot)$ for which $\Phi(t)\geq 0$ and $\sup\limits_{t\geq
0}{\mathbb{E}}|\Phi(t)|<\infty$. Then for any
$(x,\alpha)\in{\mathbb{R}}_{+}\times{\mathcal{M}}$ and $\phi=\Phi(0)$ one has

$\widehat{V}^{h}(\phi,x,\alpha)\to\widehat{V}(\phi,x,\alpha)$

as $h\to 0$. Furthermore, for sufficiently small $h$, a near-optimal
harvesting-stocking strategy of the controlled Markov chain
$\left\\{(\Phi^{h}_{n},X^{h}_{n},\alpha^{h}_{n})\right\\}$ is also a near-
optimal harvesting-stocking policy of
$\big{(}\Phi(\cdot),X(\cdot),\alpha(\cdot)\big{)}$.

###  3.3. The combined effects of seasonality and Markovian switching

In this section, we focus on another extension of (2.3) in which the
population model is periodic. This is very natural if one considers that
seasonal effects are periodic and strongly influence the dynamics of species.

There have been multiple papers treating periodic environments in ecology
(Cushing, 1977, 1980; Henson & Cushing, 1997; Rinaldi et al., 1993). Some of
these have shown how periodic forcing can create interesting new phenomena. In
White & Hastings (2020) the authors present a synthesis on the important role
of seasonality in ecology. They explain how our knowledge on seasonal dynamics
is limited both empirically and theoretically. Few studies have looked at the
joint effects of periodic and random fluctuations. The current section is a
first step in the direction of finding the optimal harvesting-stocking
strategies when both random and seasonal effects are taken into account.

If we include seasonality the dynamics is given by

(3.4)
$X(t)=x+\int\limits_{0}^{t}\big{(}b(s,X(s),\alpha(s))-U(s)\big{)}ds+\int\limits_{0}^{t}\sigma\big{(}s,X(s),\alpha(s)\big{)}dw(s),$

where $b(\cdot,x,\alpha)$, $\sigma(\cdot,x,\alpha)$ are periodic with period
$T>0$. We also suppose that the price-cost function $p(t,x,\alpha,u)$ is
periodic with period $T$ for each $(x,\alpha,u)$. The performance function
becomes

(3.5)
$\widetilde{J}(x,\alpha,U(\cdot)):={\mathbb{E}}_{x,\alpha}\int_{0}^{\infty}e^{-\delta
s}p\big{(}s,X(s),\alpha(s),U(s)\big{)}ds,$

and the value function is defined in the standard way. Since the functions
$b(\cdot)$, $\sigma(\cdot)$, and $p(\cdot)$ are time-dependent and the time
horizon of the control problem is infinite, the approach we used for the
previous formulations no longer works. In order to treat this case, we
introduce a deterministic function $\Gamma(\cdot)$ to capture the time and
convert it into the interval $[0,T)$. Thus, we will study the total expected
discounted value starting from any time $\gamma\in[0,T)$. Specifically, we
consider a generalization of (3.4)–(3.5) given by

(3.6) $\begin{split}&\\!\\!\\!\displaystyle
X(t)=x+\int\limits_{\gamma}^{t}\big{(}b(\Gamma(s),X(s),\alpha(s))-U(s)\big{)}ds+\int\limits_{\gamma}^{t}\sigma(\Gamma(s),X(s),\alpha(s))dw(s),\\\
&\\!\\!\\!\displaystyle\Gamma(t)=t-mT,\quad\text{ where }m\in\mathbb{Z}_{\geq
0}\text{ such that }t\in[mT,mT+T).\end{split}$

We consider the combined process
$\big{(}\Gamma(\cdot),X(\cdot),\alpha(\cdot)\big{)}$. The performance function
and the value function are given by

(3.7)
$\begin{split}\widetilde{J}(\gamma,x,\alpha,U(\cdot)):&={\mathbb{E}}_{\gamma,x,\alpha}\int_{\gamma}^{\infty}e^{-\delta
s}p\big{(}\Gamma(s),X(s),\alpha(s),U(s)\big{)}ds,\\\
\widetilde{V}(\gamma,x,\alpha)&=\sup\limits_{U(\cdot)\in\mathcal{A}_{\gamma,x,\alpha}}\widetilde{J}\big{(}\gamma,x,\alpha,U(\cdot)\big{)},\end{split}$

where ${\mathbb{E}}_{\gamma,x,\alpha}$ denotes the expectation with respect to
the probability law of the process $\big{(}\Gamma(t),X(t),\alpha(t)\big{)}$
having initial conditions
$\big{(}\Gamma(\gamma),X(\gamma),\alpha(\gamma)\big{)}=(\gamma,x,\alpha)$.

####  3.3.1. Constructing the approximation scheme.

In this case, we need to use two positive parameters $h_{1}$ and $h_{2}$,
where $T$ is a multiple of $h_{1}$. Let $h=(h_{1},h_{2})$. We construct a
discrete-time controlled Markov chain
$\left\\{(\Gamma^{h}_{n},X^{h}_{n},\alpha^{h}_{n}):n\in\mathbb{Z}_{\geq
0}\right\\}$ with state space $\widetilde{S}_{h}\times{\mathcal{M}}$, where

$\widetilde{S}_{h}:=\left\\{(\gamma,x)=(k_{1}h_{1},k_{2}h_{2})^{\prime}\in\mathbb{R}^{2}:k_{i}\in\mathbb{Z}_{\geq
0},k_{1}\leq T/h_{1}\right\\}.$

Note that we add $\\{\Gamma^{h}_{n}\\}$ to approximate the time
$\Gamma(\cdot)$. At each step $n$, we need to specify a control
$U^{h}_{n}\in\mathcal{U}$. We need to define transition probabilities
$q^{h}\big{(}(\gamma,x,k),(\lambda,y,l)|u\big{)}$ and interpolation intervals
$\Delta t(\gamma,x,k,u)$ so that the controlled Markov chain
$\left\\{(\Gamma^{h}_{n},X^{h}_{n},\alpha^{h}_{n})\right\\}$ is locally
consistent with respect to the controlled diffusion (3.6).

For $(\gamma,x,\alpha)\in\widetilde{S}_{h}\times{\mathcal{M}}$, the class
$\mathcal{A}^{h}_{\gamma,x,\alpha}$ of all admissible control sequences
$U^{h}$ for initial state $(\gamma,x,\alpha)$ can be defined as before. Then
the performance function for the controlled Markov chain is defined as

$\widetilde{J}^{h}(\gamma,x,\alpha,U^{h})={\mathbb{E}}\sum_{m=0}^{\infty}e^{-\delta
t_{m}^{h}}p(\Gamma^{h}_{n},X^{h}_{m},\alpha^{h}_{m},U^{h}_{m})\Delta
t_{m}^{h}.$

The value function of the controlled Markov chain is

$\widetilde{V}^{h}(\gamma,x,\alpha)=\sup\limits_{U^{h}\in\mathcal{A}_{\gamma,x,\alpha}}\widetilde{J}^{h}(\gamma,x,\alpha,U^{h}).$

It should be noted from the periodicity of the problem that for any
$(\gamma,x,\alpha)\in[0,T)\times\mathbb{R}_{+}\times{\mathcal{M}}$ one has

$\widetilde{V}(\gamma,x,\alpha)=e^{-\delta
T}\widetilde{V}(\gamma+T,x,\alpha),\quad\widetilde{V}^{h}(\gamma,x,\alpha)=e^{-\delta
T}\widetilde{V}^{h}(\gamma+T,x,\alpha).$

This property will allow us to work with the compact time interval $[0,T]$
instead of the entire infinite horizon $[0,\infty).$ The convergence result is
given below.

######  Theorem 3.2.

Suppose the following assumptions hold:

  1. (1)

The functions $b(\cdot,\cdot,\alpha)$ and $\sigma(\cdot,\cdot,\alpha)$ are
locally Lipschitz continuous for each $\alpha\in{\mathcal{M}}$.

  2. (2)

For any initial condition
$(\gamma,x,\alpha)\in[0,T)\times{\mathbb{R}}_{+}\times{\mathcal{M}}$, the
uncontrolled system (3.4) has a unique global solution.

  3. (3)

The price $p(\cdot,\cdot,\alpha,\cdot)$ is bounded and continuous for all
$\alpha\in{\mathcal{M}}$.

  4. (4)

For each $(x,\alpha)\in\mathbb{R}_{+}\times{\mathcal{M}}$ and
$u\in\mathcal{U}$, the functions $b(\cdot,x,\alpha)$,
$\sigma(\cdot,x,\alpha)$, $p(\cdot,x,\alpha,u)$ are periodic with period
$T>0$.

Then for any
$(\gamma,x,\alpha)\in[0,T)\times{\mathbb{R}}_{+}\times{\mathcal{M}}$, we have

$\widetilde{V}^{h}(\gamma,x,\alpha)\to\widetilde{V}(\gamma,x,\alpha)$

as $h\to 0$. Furthermore, for sufficiently small $h$, a near-optimal
harvesting-stocking strategy of the controlled Markov chain
$\left\\{(\Gamma^{h}_{n},X^{h}_{n},\alpha^{h}_{n})\right\\}$ is also a near-
optimal harvesting-stocking policy of
$\big{(}\Gamma(\cdot),X(\cdot),\alpha(\cdot)\big{)}$.

##  4\. Numerical Examples

We explore some numerical examples, using various models and assumptions.
Throughout this section, we suppose the time discount rate is $\delta=0.02$.
Let $B\in(0,\infty)$ be an upper bound introduced for computational purposes.
If not specified, we will take $B=4$ in our examples. We assume that the
environment only switches between two states, so ${\mathcal{M}}=\\{1,2\\}$. We
exhibit various price-cost formulations to explore for features, and possibly
new phenomena, in the various considered models.

###  4.1. Introducing switching in a Verhulst-Pearl system

A model that is extensively used in biology is the Verhulst-Pearl one. The
dynamics that includes switching is given by

$dX(t)=X(t)\big{(}\mu(\alpha(t))-\kappa(\alpha(t))X(t)\big{)}\,dt+\sigma\big{(}\alpha(t)\big{)}X(t)dw(t).$

Here $\mu(\alpha)$ is the growth rate in environment $\alpha$,
$\kappa(\alpha)$ is the intraspecific competition rate in environment
$\alpha$, and $\sigma^{2}(\alpha)$ is the variance of the per-capita
environmental fluctuations. Suppose there are $m_{0}$ environmental states in
environment $\alpha(\cdot)$. One can fix an environment $\alpha$ and look at
the dynamics of the population in that environment

$dX^{\alpha}(t)=X^{\alpha}(t)\big{(}\mu(\alpha)-\kappa(\alpha)X^{\alpha}(t)\big{)}\,dt+\sigma(\alpha)X^{\alpha}(t)dw(t).$

The stochastic growth rate of the species (Chesson & Ellner, 1989; Chesson,
2000; Schreiber et al., 2011; Hening & Nguyen, 2018) in environment $\alpha$
is given by

$r(\alpha)=\mu(\alpha)-\frac{\sigma^{2}(\alpha)}{2}.$

Following Evans et al. (2015), this stochastic growth rate completely
determines the long term behavior of the population in each fixed environment.
If $r(\alpha)>0$ we get the convergence to a unique stationary distribution on
$(0,\infty)$. If $r(\alpha)<0$ the population goes extinct in environment
$\alpha$ and with probability $1$ we have

$\lim_{t\to\infty}\frac{\ln X^{\alpha}(t)}{t}=r(\alpha)<0.$

If $r(\alpha)=0$, by Evans et al. (2015) the population process is null
recurrent in environment $\alpha$. This means the population does not go
extinct but also does not converge to a stationary distribution; it keeps
fluctuating between large and small values. Let $(\nu_{1},\dots,\nu_{m_{0}})$
be the stationary distribution of the Markov chain $\alpha(t)$. If

$r=\sum_{k=1}^{m_{0}}\nu_{k}\left(\mu(k)-\frac{\sigma^{2}(k)}{2}\right)>0,$

we get by Hening & Li (2020) that the population converges, if there is
switching, to a stationary distribution.

We consider a Verhulst population model with
$\mu(\alpha)=4-\alpha,\;\kappa(\alpha)=2,\;\sigma(\alpha)=1,\;(x,\alpha)\in{\mathbb{R}}_{+}\times\\{1,2\\}.$
As a result the stochastic growth rate in environment $\alpha$ is

$r(\alpha)=(4-\alpha)-\frac{1}{2},$

which implies $r(1)>0$ and $r(2)>0$, so if the environment would be fixed to
either of the $\alpha$ values, the species would converge to a unique
stationary distribution on ${\mathbb{R}}_{+}$. Suppose that the generator $Q$
of the Markov chain $\alpha(\cdot)$ is given by $q_{11}=-0.1$, $q_{12}=0.1$,
$q_{21}=0.1$, $q_{22}=-0.1$. This implies that the stationary distribution of
$\alpha(t)$ is $(\nu_{1},\nu_{2})=(0.5,0.5)$. Henceforth, let the set of
controls be $\mathcal{U}=\\{u:u=k/500,k\in\mathbb{Z},-1000\leq k\leq 1500\\}$
where unspecified.

In a first experiment, we compare this model to a baseline model, where
$\mu(\alpha)=2.5$, which is in between the two values for the switching
environments. To isolate the effect of switching, in our first example we keep
the price and cost functions simple: $P(x,\alpha)=1$ and $C(x,\alpha,u)=0$.

We have presented the formulation theoretically in Section 2, and we describe
a detailed procedure in Section 2.1 and Appendix A. In the following numerical
formulations, we will apply similar procedures, described in the same
sections. For an admissible strategy $U(\cdot)$ we have

$J\big{(}x,\alpha,U(\cdot)\big{)}={\mathbb{E}}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X(s),\alpha(s),U(s)\big{)}ds.$

Based on the algorithm constructed above and in the Appendixes, we carry out
the computation by using the methods in Kushner & Dupuis (1992). We take the
initial control $U_{0}(x,\alpha)\equiv 0$ and set the initial values
$V^{h}_{0}(x,\alpha)\equiv 0$. We outline how to find the sequence of values
of $V_{n}^{h}(\cdot)$ as follows. At each level $x=h,2h,\dots,B$,
$\alpha\in{\mathcal{M}}$, and control $u\in\mathcal{U}$, we compute

$V_{n+1}^{h}(x,\alpha\,|\,u)=e^{-\delta\Delta
t^{h}(x,\alpha,u)}\sum\limits_{(y,\beta)\in
S_{h}\times{\mathcal{M}}}V^{h}_{n}(y,\beta)q^{h}\big{(}(x,\alpha),(y,\beta)|u\big{)}+p(x,\alpha,u)\Delta
t^{h}(x,\alpha,u).$

Working with the compact set state space $[0,B]$, we use reflection if $x=B$;
that is, $V^{h}_{n}(B+h,\beta)=V^{h}_{n}(B,\beta)$ for any
$\beta\in{\mathcal{M}}$. Then we choose the control $U^{h}_{n+1}(x,\alpha)$
and record an improved value $V^{h}_{n+1}(x,\alpha)$ by

$U^{h}_{n+1}(x,\alpha)=\hbox{argmax}_{u\in\mathcal{U}}V_{n+1}^{h}(x,\alpha\,|\,u),\quad
V^{h}_{n+1}(x,\alpha)=V^{h}_{n+1}(x,\alpha\,|\,U^{h}_{n+1}(x,\alpha)).$

The iterations stop as soon as the increment
$V^{h}_{n+1}(\cdot)-V^{h}_{n}(\cdot)$ reaches some tolerance level. We set the
error tolerance to be $10^{-8}$.

(a)  (b)

Figure 1. Value function (left) and optimal harvesting-stocking rate (right)
for a model with switching affecting $\mu(\alpha)=4-\alpha$, compared with a
baseline model with no switching and $\mu(\alpha)=2.5.$

The numerical result is shown in Figure 1. The value function in the left
panel is increasing in both the switching and the baseline models, as
expected. Moreover, it is concave, which implies that the marginal value of
one small unit of population decreases. This is expected because population
growth becomes less favourable as the population size increases—a result of
intra-species competition.

Because the transition probability rates were chosen to be small, the value
function of the baseline model is in-between the value functions for the two
$\alpha$ states of the switching model. When the transition probability rates
increase, that will not always be the case, because the value function in one
state is determined, in part, by the value function in the other state.

Another expected result is that the optimal harvesting rate in the second
panel of Figure 1 is bang-bang—it switches between extremes at a certain
population threshold. This should be the case, because the cost of harvesting-
stocking is 0, and therefore independent of the value. This was shown
theoretically for a model without switching in Hening, Nguyen, Ungureanu &
Wong (2019). In the second panel, we also see that the optimal stocking-
harvesting transition happens at a lower population level in the growth-
unfavorable state $\alpha=2$, when compared with the intermediate growth
baseline and the growth-favorable state $\alpha=1$. Intuitively, lower growth
prospects imply optimal extraction should start at a lower population level.

Numerical experiments where $\sigma(\alpha)$ and $\kappa(\alpha)$ are also
allowed to depend on $\alpha$ show results with similar features, so we omit
them for brevity.

###  4.2. Analysis of the effect of the control cost

We expect to find all or nothing harvesting-stocking when the cost is not
convex in the rate, even when there is environmental switching. Using the same
set-up as above, but now with a convex cost function $C(x,\alpha,u)=u^{2}/2$,
we find that the optimal control is not bang-bang anymore (Figure 2). Further
numerical experiments (Figures 9, 10, Appendix D) support the following
conjecture, which says that, with cost functions that are not convex, the
optimal strategies are bang-bang, with thresholds that may depend on the
states $\alpha$.

######  Conjecture 4.1.

Suppose we have one species that evolves according to (2.3) and suppose
Assumption 2.1 holds, the price $P$ only varies with the environment, the cost
does not depend on the population size, $C(x,\alpha,u)=C(\alpha,u)$ and the
cost is not convex in the harvesting rate $u$. Furthermore, assume that in the
absence of harvesting the species persists, i.e.,

$\lambda(\mu_{0})=\sum_{k=1}^{m_{0}}\nu_{k}\left(\mu(k)-\frac{\sigma^{2}(k)}{2}\right)>0.$

One can construct the optimal harvesting strategy $U^{*}$ as follows. There
exists a threshold $0\leq u^{*}(k)\leq\infty$ for the population such that, if
$X(0-)=x$, then

(4.1) $\displaystyle U^{*}(t)$
$\displaystyle=\begin{cases}\inf\;\mathcal{U},&\mbox{\text{ if }
$t>0,\;\alpha(t)=k,\;X(t)<u^{*}(k)$},\\\ \sup\;\mathcal{U},&\mbox{\text{ if }
$t>0,\;\alpha(t)=k,\;X(t)\geq u^{*}(k)$}.\end{cases}$

(a)  (b)

Figure 2. Value function (left) and optimal harvesting-stocking rate (right)
for a model with switching affecting $\mu(\alpha)=4-\alpha$, and convex cost
$C(x,\alpha,u)=u^{2}/2$, compared with a baseline model with no switching,
$\mu(\alpha)=2.5$, and the same cost function.

From the right panel of Figure 2, we note that the growth favorable state,
$\alpha=1$, is more conducive to expensive harvesting and stocking, which is
intuitive and expected. The value functions are again concave and increasing
in population, as expected. This will be common across all numerical examples
considered. Henceforth we may omit the value function from figures.

###  4.3. Large and small transition rates with switching

In the left panel of Figure 2, the apparent overlap between the value function
of the baseline model and the switching model in the state $\alpha=1$ is
coincidental. Numerical experiments show that increasing the switching
probability rates pulls the value function graphs for $\alpha=1,2$ toward each
other. If the switching rates were very small, the green baseline value
function would be sandwiched between the red and blue value functions. As the
switching rates become large, the value functions of the two switching states
will overlap in the limit.111The value function in a switching state is
influenced by the value function in the other state, since there is a
probability of transition. Figure 3 shows this process for the same
specification as above, when $q_{12}=q_{21}\in\\{0,0.01,0.1,1,10,\infty\\}$.
The following remark confirms the numerical experiments that suggest
convergence.

(a)  (b)

Figure 3. The left and right panels show the same value function graph, but on
different vertical scales. The model is discussed in Sections 4.2, 4.3. The
green solid curves in both panels corresponds to a model without switching and
$\mu=2.5$, or to a model with $\mu=4-\alpha$ and infinitely large switching
rates. The blue long-dashed curves in the left panel show models with no
switching and $\mu=3$, $\mu=2$ respectively, or equivalently a model with
$\mu=4-\alpha$ and infinitely slow switching rates. The red dashed curves in
the left panel are for a model with $\mu=4-\alpha$ and $q_{12}=q_{21}=0.01$.
The black dotted curves in both panels correspond to the same model, but
$q_{12}=q_{21}=0.1$. The teal dot-dashed curves in the right panel correspond
to $q_{12}=q_{21}=1$, and the orange dashed curves in the right panel to
$q_{12}=q_{21}=10$.

######  Remark 4.2.

If we divide the switching rates by $\varepsilon$ and let
$\varepsilon\downarrow 0$ we have a slow-fast dynamics where the switching
happens on a much shorter time scale. One can show that the dynamics without
harvesting (2.1) will converge (Hening & Li, 2020) to

(4.2) $d\bar{\xi}(t)=\bar{b}(\bar{\xi}(t))dt+\bar{\sigma}(\xi(t))dw(t),$

where

(4.3) $\bar{b}(\mathbf{x})=\sum_{k=1}^{m_{0}}b(k)\nu_{k}.$

and

(4.4) $\bar{\sigma}=\sqrt{\sum_{k=1}^{m_{0}}\sigma^{2}(k)\nu_{k}}.$

Here $(\nu_{1},\dots,\nu_{m_{0}})$ is the stationary distribution of the
Markov chain $\alpha(t)$. We conjecture that the optimal stocking-harvesting
strategies $U^{\varepsilon,*}$ of the system (2.3) converge as
$\varepsilon\downarrow 0$ to the optimal harvesting-stocking strategy
$\bar{U}^{*}$ of the system

(4.5)
$\bar{X}(t)=x+\int\limits_{0}^{t}\big{(}\bar{b}(\bar{X}(s),\alpha(s))-\bar{U}(s)\big{)}ds+\int\limits_{0}^{t}\sigma(\bar{X}(s),\alpha(s))dw(s).$

Further numerical experiments, which we omit, show that if the switching rates
are such that they favour one state over the other in the stationary
distribution, the value function of a simplified system with only the favoured
state will be closer to the value function with switching, than the value
function of a simplified system with only the unfavoured state. That is
because the environment is mostly in the favored state, and its conditions are
therefore more relevant to the valuation.

###  4.4. Switching in the cost and price functions

We have shown how the dependence of the cost function on the harvesting rate
determines whether harvesting is bang-bang. It turns out that non-trivial
price functions may also smooth out the optimal harvesting rate. We have
looked at three typical functional forms for a per-unit price-cost function
$p(x,\alpha,u)=p_{0}(u)$, two of which can capture monopolistic competition,
i.e., imperfect competition. The first is a simple constant price function,
$p_{0}(u)=p,\;p>0,$ which is appropriate when the market influence of the
extraction operation is low. The second has a piecewise linear decreasing form
given by

(4.6)
$p_{0}(u)=\begin{cases}\overline{p},&\kappa_{1}-\kappa_{2}u\geq\overline{p},\\\
(\kappa_{1}-\kappa_{2}u),\quad 0<&\kappa_{1}-\kappa_{2}u<\overline{p},\\\
0,&\text{else}.\end{cases}$

The third has a constant price elasticity of demand form, and is given by

(4.7)
$p_{0}(u)=\begin{cases}\overline{p},&\kappa_{1}|\kappa_{2}+u|^{1/\epsilon}\geq\overline{p},\\\
\kappa_{1}|\kappa_{2}+u|^{1/\epsilon},\quad
0\leq&\kappa_{1}|\kappa_{2}+u|^{1/\epsilon}<\overline{p}.\\\ \end{cases}$

The functions represent instantaneous market demand, and obey common sense
assumptions. In particular, the price decreases as the supplied quantity
increases, and the price is always positive. Both forms are typical choices in
modeling market behavior. The economic model implied by the price dependence
is the following: the instantaneous harvested rate $u$ is a quantity placed in
a common market, where the total supplied by all sources is what determines
the instantaneous common final price $u$, and hence the returns.
The constants $\kappa_{1},\kappa_{2}$ are positive. They capture, among other
things, the contribution to total harvested quantities from different sources.
In the third formulation, $\epsilon$ is negative and represents the price
elasticity. The constant $\overline{p}$ is added to account for unbounded
negative flow rates (e.g., unrealistically large stocking events) in the
economic demand model. Large negative values of $u$ could make the overall
market supply negative—something evidently not possible—so we ignore such
events for practical applications. That is, we want to keep prices below
$\overline{p}$ in our examples.

Figure 4 shows, as an example, the harvesting strategy corresponding to two
price functions, also known as demand functions, based on the functional forms
discussed. The left panel corresponds to linear demand, $p(u)=1-u/4$, and the
right panel to $p(u)=(1+u/3)^{-1}$. The second form is for a good with unit
elasticity.222The (price) elasticity of a demand function is a measure of the
responsiveness of the quantity demanded to price changes, measured in an unit
free manner, defined as $\epsilon=\frac{du}{dp}\cdot\frac{p}{u}$, and widely
used in economics. $\epsilon=-1$ is defined as “unit elasticity.” The
constants were chosen so that price will not be negative in the harvesting
range, but are otherwise arbitrary. Different choices of slopes and elasticity
values lead to similar qualitative conclusions—with no cost, when $up(u)$ is
concave, we have a harvesting strategy that is continuous in the population
level.

(a)  (b)

Figure 4. The left and right panels show the optimal harvesting functions for
the model in Section 4.2, for two price functions, $p(u)=1-u/4$ on the left,
and $p(u)=(1+u/3)^{-1}$ on the right. There is no dependency on $\alpha$ here,
and $\mu=2.5$. The other parameters are as before.

###  4.5. Interlude with intuition and conjecture

In the previous examples, the interpretation of the control value was
“population units extracted per unit of time.” Based on marginal cost-benefit
analysis and intuition, we suspect that a consequence of such a cost
formulation, together with others that ensure that the value functions is
increasing and concave, is that the optimal harvesting strategies are
monotonically increasing in the population level (see also Figures 1, 2 and
4).

Our intuition is as follows: Let’s take for granted that the value function is
increasing and concave in the population variable. Consider the cost-benefit
of seeding an extra unit of population, or harvesting one less. If the
population increases, the sale value of one future unit of population is the
same, but future growth is less favourable at higher population levels,
because of intra-species competition. Taking into account both, we can say
that we have a decreasing marginal benefit of having an extra unit of
population, as population increases. We can also see it in the concave shape
of the value function. Its rate of increase per unit of population decreases
as population increases. Moreover, with constant effort harvesting
strategies—that is, cost dependent on the absolute value of stocking or
harvesting—marginal cost per unit does not depend on the population level. The
difference between marginal benefit and marginal cost is also monotonic
decreasing, which means the net marginal benefit of one unit of population is
decreasing in the population. This implies that marginally increasing the
population level would make stocking less favourable or harvesting more
favourable, so we have a monotonic increasing harvesting strategy.

We have not found a way to restrict our general model to provide a clear
conjecture statement, mainly because our model only allows for finite rates of
harvesting and stocking—which we consider more realistic. However, there are
many suggestive results in the literature that are supportive of our thinking,
if we relax this assumption and allow possibly infinite stocking-harvesting
rates. Song et al. (2011) show in Theorem 4.4 that the value function is
continuous for a model similar to our baseline model, with switching. Alvarez
E. & Shepp (1998); Alvarez (2000), working with a Verhulst-Pearl model and a
general diffusion model respectively, with no switching, constant price-cost
$p$, and a possibly infinite rate of harvesting, show that the value function
is continuous and increasing. Hening, Tran, Phan & Yin (2019), working with a
multi-dimensional population system with harvesting and stocking, no
switching, state-dependent price-cost $p$, and a possibly infinite rate of
harvesting, show in Proposition 2.5 that the value function is increasing and
Lipschitz continuous. It can be easily shown that Proposition 2.5 from Hening,
Tran, Phan & Yin (2019) can be generalised to a model with switching.

###  4.6. Variable effort harvesting examples

Monotonic strategies are evidently simpler to interpret, to implement, and to
parametrize. But we want to know what happens to a model if control is borne
differently. What if the cost depends on the intensity of extraction, rather
than on the extracted quantity itself? Numerical experiments, e.g., the one in
Figure 5, show that the optimal harvesting strategy may not be monotonic in
such cases; for example when the cost of control applies to the fractional
harvesting or seeding time rate. That is because the cost incurred through
harvesting or seeding per unit of population, per unit of time, can now
decrease as the population increases—even with an increasing cost function!
With a possibly decreasing marginal cost of extraction, we may have a non-
monotonic net marginal benefit.

This makes the v-shape in the right panel possible. If the population stock is
very low, while growth is favourable, the cost per unit seeded is not
favourable at all. As the population grows, seeding costs per unit go down, so
the optimal strategy reflects that. After further increases, costs stay low,
but growth prospects diminish as well. Hence optimal control slowly shifts
from seeding to extraction. The adaptation of our general model for such a
type of control was described in Section 4.6. The discussion of the numerical
method is detailed in Appendix A.2.

(a)  (b)

Figure 5. The left and right panels show the value function and optimal
harvesting rate for a model with variable effort harvesting, and switching
that affects the population growth rate. $\mu(\alpha)=2-\alpha/2$,
$\kappa(\alpha)=1/2$, and the cost function is $u^{2}$. The other parameters
are as in Section 4.1. For simplicity, we show only the value function and the
harvesting strategy for $\alpha=1$. The other state has similar values and
controls.

###  4.7. Stochastic pricing examples

A realistic model of population management may want to take into account that
the value of the resource varies over time. We have already considered price
dependence on gross quantity extracted. But price can also change over time
because of external factors. For example, usually economists model the price
of a commodity as a stochastic variable (Osborne, 1959; Miltersen, 2003). If
this variability is large enough, we would like to understand what effects it
has on our optimal population management strategy. The following is not the
standard model of price stochasticity used in finance, but it is convenient
for our set-up.

Consider the price-cost function

$p(t,x,\alpha,u)=\left(p_{0}(\alpha,u)+\Phi(t)\right)u-C(t,x,\alpha,u),\quad(x,\alpha)\in\mathbb{R}_{+}\times\mathcal{M},t\in\mathbb{R}^{+},$

where $p_{0}(\cdot)=1$, and $C(t,x,\alpha,u)=u^{2}/2$. If we ignore random
perturbations, the evolution of $\Phi(\cdot)$ is given by

$d\Phi(t)=\Phi(t)\big{(}0.4-\Phi(t)\big{)}dt,\quad\Phi(0)\in[0,0.4).$

It can be seen that $\Phi(t)\in(0,0.4)$ for any $t\geq 0$ and $\Phi(t)\to 0.4$
as $t\to\infty$. Now we suppose that the evolution of $\Phi(t)$ is subject to
a white noise $w_{0}(\cdot)$ and is given by

$d\Phi(t)=\Phi(t)\big{(}0.4-\Phi(t)\big{)}dt+0.5\Phi(t)\big{(}0.4-\Phi(t)\big{)}dw_{0}(t).$

We also have $\Phi(t)\in[0,0.4)$ with probability one.

(a)  (b)

Figure 6. Harvesting rate in regime $\alpha=1$, 3D plot on the left, contour
plot on the right. Stochastic price given by $p_{0}(\alpha,u)+\Phi(t)$, cost
by $C(u)=u^{2}/2$, and the other parameters as in the baseline model.
Discussion in Section 4.7.

The harvesting policies for the two regimes share the same shape, so Figure 6
shows the harvesting rate as a function of population size $x$ and the
observation of the price $\Phi(\cdot)$, for one of the regimes, $\alpha=1$,
for brevity. The value function has the usual properties, and we omit its
graph. It is increasing and concave in population for all price cross-sections
in the range. Moreover, fixing the population variable, we see that higher
prices are more conducive to both harvesting and seeding, as expected.

###  4.8. Seasonal and periodic variability

Finally, in this subsection, we consider applications of our method to models
with periodic parameters. Here, periodicity can be thought of as representing
seasonal changes in the environment. In this example, we consider a simple
adaptation of the Verhulst model, where growth is additively affected by a
sine wave. There is little justification for this functional form, except to
point out that our sine wave would be one of the two likely lowest order terms
in the Fourier expansion of any periodic dependency we may conceive.

Suppose $b(t,x,\alpha)=x\left(4-\alpha+\sin(2\pi t)-2x\right)$,
$\sigma(t,x,\alpha)=x$, $P(t,x,\alpha,u)=1$, and $C(t,x,\alpha,u)=u^{2}/2$,
where $(t,x,\alpha,u)\in\mathbb{R}^{2}_{+}\times\\{1,2\\}\times\mathcal{U}$.
To help with numerical estimation, we reduce the population range, the control
range and the sampling density. Now
$\mathcal{U}=\\{u:u=k/20,k\in\mathbb{Z},-20\leq k\leq 40\\}$. The time step
$h_{1}$ is $T/4000$. The other parameters are as before. Note that
$b(\cdot,x,\alpha)$, $\sigma(\cdot,x,\alpha)$, and $p(\cdot,x,\alpha,u)$ are
periodic with period $T=1$.

(a)  (b)

Figure 7. Harvesting rate in regime $\alpha=1$, 3D plot on the left, contour
plot on the right. This corresponds to an adapted Verhulst model, with drift
affected by a sine wave, $b(t,x,a)=x(4-\alpha-\sin(2\pi t)-2x)$,
$\sigma(t,x,\alpha)=x$, constant pricing and quadratic cost, $C(u)=u^{2}/2$.
Discussion in subsection 4.8.

Figure 7 shows the estimation results for the harvesting rate as a function of
population size $x$ and the time $t$. As we can see the periodicity of the
growth rate is apparent in the graphs. It is interesting that the maximal and
minimal effects on the harvesting rate do not match the maximum and minimum of
the sine wave at $t=1/4$ and $t=3/4$. Not only is the effect of harvesting out
of phase with the sine wave, the effects on harvesting and stocking seem to be
out of phase with each other. At this point, we are not sure if this
observation is an artefact of the limited sampling in the numerical
estimation, or a robust result. We omit the value function and the harvesting
regime for $\alpha=2$ for brevity.

(a)  (b)

Figure 8. Harvesting rate in regime $\alpha=1$, 3D plot on the left, contour
plot on the right. The parameters correspond to the baseline Verhulst model,
$b(t,x,a)=x(4-\alpha-2x)$, $\sigma(t,x,\alpha)=x$, price is constant, and the
quadratic cost is affected by a sine wave, $C(t,u)=(1+\sin(2\pi t))u^{2}/2$.
$h_{1}=T/16000$, and $\mathcal{U}\in[-1,-1]$. Discussion in subsection 4.8.

Seasonality can affect other parameters as well. In the following example, we
have seasonality affecting costs, but not the growth of the population. Figure
8 shows what are intuitive results. Harvesting when the cost is high is low,
and vice-versa. Again, we omit the value function and the harvesting regime
for $\alpha=2$ for brevity.

##  5\. Discussion

The objective of this paper was to develop a theoretically sound grounding for
numerical methods, that would allow us to analyze realistic, complex,
stochastic population models. Furthermore, we wanted to use these methods to
explore the effects of adding features or extensions to classic models like
Verhulst-Pearl, Gompertz and Nisbet. Among the extensions of interest that we
explored are Markovian environmental switching, periodicity, and non-trivial
harvesting and stocking cost and price dependencies. These are very important
in any realistic model of population management.

We showed that the optimal harvesting and stocking strategy is not always of
bang-bang type—that is, all or nothing. The crucial factor seems to be whether
the dependence of the cost function is convex in the harvested time rate or
not. A similar observation applies to price dependency. Regarding
environmental switching, we determined that, in the limit of switching rates
favouring one state versus another, the value function and harvesting rate
converge as expected to the ones of the favoured state. Moreover, in the limit
of fast switching, the value function and harvesting rate converge to the one
representing an “average environment,” as is expected from the slow-fast
dynamics analysis of stochastic systems without harvesting (Yin & Zhu, 2009;
Hening & Li, 2020; Du et al., 2021).

In general, we have monotonicity of the value function in the population size,
but also of the harvesting rate. The harvesting rate is increasing in the
population size, which is intuitively expected, as a higher population should
be more suitable for harvesting, everything else being equal. However, this
observation breaks down when the cost of harvesting depends on the proportion
harvested from the available population per unit of time, as in the so-called
variable effort strategies.

Seasonal variation in the model parameters also has implications that are
expected or intuitive. The maximum harvesting or stocking values correspond to
times of maximal fertility.333Observe that we have not considered models with
a distinction between young and old. All individuals are the same at all
points in time, and have the same harvested value. In the slow interval of the
season, both harvesting and stocking are reduced.

Allowing the harvesting strategy to depend on a stochastic price variable
leads to a variable strategy where harvesting is favoured by high prices, and
stocking is favoured by low prices, as expected.

Overall, we think that proper numerical techniques are essential to test model
validity by numerical experiments, especially when the complexity of the model
does not favour finding explicit solutions. Moreover, such numerical methods
are good at finding points of interest in parameter space, and for formulating
conjectures.

Acknowledgments. A. Hening is supported by the NSF through the grant DMS
1853463. K. Q. Tran is supported by the National Research Foundation of Korea
grant funded by the Korea Government (MIST) NRF-2021R1F1A1062361.

## References

  * (1)
  * Abakuks & Prajneshu (1981) Abakuks, A. & Prajneshu (1981), ‘An optimal harvesting policy for a logistic model in a randomly varying environment’, Math. Biosci. 55(3-4), 169–177.
  * Alvarez E & Hening (2020) Alvarez E, L. H. & Hening, A. (2020), ‘Optimal sustainable harvesting of populations in random environments’, Stochastic Processes and their Applications .
  * Alvarez E. & Shepp (1998) Alvarez E., L. H. R. & Shepp, L. A. (1998), ‘Optimal harvesting of stochastically fluctuating populations’, J. Math. Biol. 37(2), 155–177.
  * Alvarez (2000) Alvarez, L. H. (2000), ‘Singular stochastic control in the presence of a state-dependent yield structure’, Stochastic processes and their applications 86(2), 323–343.
  * Alvarez & Koskela (2007) Alvarez, L. H. R. & Koskela, E. (2007), ‘Optimal harvesting under resource stock and price uncertainty’, Journal of Economic Dynamics and Control 31, 2461–2485.
  * Alvarez & Shepp (1998) Alvarez, L. H. R. & Shepp, L. A. (1998), ‘Optimal harvesting of stochastically fluctuating populations’, J. Math. Biol 37, 155–177.
  * Asche et al. (2015) Asche, F., Chen, Y. & Smith, M. D. (2015), ‘Economic incentives to target species and fish size: prices and fine-scale product attributes in norwegian fisheries’, ICES Journal of Marine Science 72(3), 733–740.
  * Bao & Shao (2016) Bao, J. & Shao, J. (2016), ‘Permanence and extinction of regime-switching predator-prey models’, SIAM Journal on Mathematical Analysis 48(1), 725–739.
  * Benaïm & Lobry (2016) Benaïm, M. & Lobry, C. (2016), ‘Lotka Volterra in fluctuating environment or “how switching between beneficial environments can make survival harder”’, Ann. Appl. Probab. .  to appear.
  * Bohner & Streipert (2016) Bohner, M. & Streipert, S. (2016), ‘Optimal harvesting policy for the beverton–holt model’, Mathematical Biosciences & Engineering 13(4), 673.
  * Bourquin (2021) Bourquin, A. (2021), ‘Persistence in randomly switched lotka-volterra food chains’, arXiv preprint arXiv:2109.03003 .
  * Brauer & Sànchez (2003) Brauer, F. & Sànchez, D. A. (2003), ‘Periodic environments and periodic harvesting’, Natural Resource Modeling 16(3), 233–244.
  * Braverman & Mamdani (2008) Braverman, E. & Mamdani, R. (2008), ‘Continuous versus pulse harvesting for population models in constant and variable environment’, Journal of Mathematical Biology 57(3), 413–434.
  * Chesson (2000) Chesson, P. (2000), ‘General theory of competitive coexistence in spatially-varying environments’, Theoretical Population Biology 58(3), 211–237.
  * Chesson & Ellner (1989) Chesson, P. L. & Ellner, S. (1989), ‘Invasibility and stochastic boundedness in monotonic competition models’, Journal of Mathematical Biology 27(2), 117–138.
  * Clark (2010) Clark, C. W. (2010), Mathematical Bioeconomics, Pure and Applied Mathematics (Hoboken), third edn, John Wiley & Sons, Inc., Hoboken, NJ.  The mathematics of conservation.
  * Cohen et al. (2021) Cohen, A., Hening, A. & Sun, C. (2021), ‘Optimal ergodic harvesting under ambiguity’, arXiv preprint arXiv:2104.10626 .
  * Cromer (1988) Cromer, T. (1988), ‘Harvesting in a seasonal environment’, Mathematical and Computer Modelling 10(6), 445–450.
  * Cushing (1977) Cushing, J. M. (1977), ‘Periodic time-dependent predator-prey systems’, SIAM Journal on Applied Mathematics 32(1), 82–95.
  * Cushing (1980) Cushing, J. M. (1980), ‘Two species competition in a periodic environment’, Journal of Mathematical Biology 10(4), 385–400.
  * Du et al. (2021) Du, N. H., Hening, A., Nguyen, D. H. & Yin, G. (2021), ‘Dynamical systems under random perturbations with fast switching and slow diffusion: Hyperbolic equilibria and stable limit cycles’, Journal of Differential Equations 293, 313–358.
  * Evans et al. (2015) Evans, S. N., Hening, A. & Schreiber, S. J. (2015), ‘Protected polymorphisms and evolutionary stability of patch-selection strategies in stochastic environments’, J. Math. Biol. 71(2), 325–359.
  * Fan & Wang (1998) Fan, M. & Wang, K. (1998), ‘Optimal harvesting policy for single population with periodic coefficients’, Mathematical biosciences 152(2), 165–178.
  * Hening & Li (2020) Hening, A. & Li, Y. (2020), ‘Stationary distributions of persistent ecological systems’, arXiv preprint arXiv:2003.04398 .
  * Hening & Nguyen (2018) Hening, A. & Nguyen, D. (2018), ‘Coexistence and extinction for stochastic Kolmogorov systems’, Annals of Applied Probability 28(3), 1893–1942.
  * Hening & Nguyen (2020) Hening, A. & Nguyen, D. H. (2020), ‘The competitive exclusion principle in stochastic environments’, Journal of Mathematical Biology 80, 1323––1351.
  * Hening et al. (2021) Hening, A., Nguyen, D. H., Nguyen, N. & Watts, H. (2021), ‘Random switching in an ecosystem with two prey and one predator’, arxiv .
  * Hening, Nguyen, Ungureanu & Wong (2019) Hening, A., Nguyen, D. H., Ungureanu, S. C. & Wong, T. K. (2019), ‘Asymptotic harvesting of populations in random environments’, Journal of mathematical biology 78(1-2), 293–329.
  * Hening & Strickler (2019) Hening, A. & Strickler, E. (2019), ‘On a predator-prey system with random switching that never converges to its equilibrium’, SIAM Journal on Mathematical Analysis 51(5), 3625–3640.
  * Hening, Tran, Phan & Yin (2019) Hening, A., Tran, K., Phan, T. & Yin, G. (2019), ‘Harvesting of interacting stochastic populations’, Journal of mathematical biology 79(2), 533–570.
  * Hening & Tran (2020) Hening, A. & Tran, K. Q. (2020), ‘Harvesting and seeding of stochastic populations: analysis and numerical approximation’, Journal of Mathematical Biology 81, 65–112.
  * Henson & Cushing (1997) Henson, S. M. & Cushing, J. M. (1997), ‘The effect of periodic habitat fluctuations on a nonlinear insect population model’, Journal of Mathematical Biology 36(2), 201–226.
  * I. Kharroubi (2019) I. Kharroubi, T. Lim, V. L. V. (2019), ‘Optimal exploitation of a resource with stochastic population dynamics and delayed renewal’, Journal of Mathematical Analysis and Applications 477(1), 627–656.
  * Kushner (1984) Kushner, H. J. (1984), Approximation and Weak Convergence Methods for Random Processes, with Applications to Stochastic Systems Theory, MIT Press.
  * Kushner (1990) Kushner, H. J. (1990), ‘Numerical methods for stochastic control problems in continuous time’, SIAM J. Control Optim 28(5), 999–1048.
  * Kushner & Dupuis (1992) Kushner, H. J. & Dupuis, P. G. (1992), Numerical methods for stochastic control problems in continuous time, Springer-Verlag.
  * Lungu & Øksendal (1997) Lungu, E. M. & Øksendal, B. (1997), ‘Optimal harvesting from a population in a stochastic crowded environment’, Math. Biosci. 145(1), 47–75.
  * Miltersen (2003) Miltersen, K. R. (2003), ‘Commodity price modelling that matches current observables: A new approach’.
  * Nguyen et al. (2017) Nguyen, D. H., Yin, G. & Zhu, C. (2017), ‘Certain properties related to well posedness of switching diffusions’, Stochastic Process. Appl. .
  * Osborne (1959) Osborne, M. F. (1959), ‘Brownian motion in the stock market’, Operations research 7(2), 145–173.
  * Pooley (1987) Pooley, S. G. (1987), ‘Demand considerations in fisheries management–hawaii’s market for bottom fish’, Tropical snappers and groupers: biology and fisheries management pp. 605–38.
  * Rinaldi et al. (1993) Rinaldi, S., Muratori, S. & Kuznetsov, Y. (1993), ‘Multiple attractors, catastrophes and chaos in seasonally perturbed predator-prey communities’, Bulletin of mathematical Biology 55(1), 15–35.
  * Schreiber et al. (2011) Schreiber, S. J., Benaïm, M. & Atchadé, K. A. S. (2011), ‘Persistence in fluctuating environments’, J. Math. Biol. 62(5), 655–683.
  * Song et al. (2006) Song, Q. S., Yin, G. & Zhang, Z. (2006), ‘Numerical methods for controlled regime-switching diffusions and regime-switching jump diffusions’, Automatica 2(7), 1147–1157.
  * Song et al. (2011) Song, Q., Stockbridge, R. H. & Zhu, C. (2011), ‘On optimal harvesting problems in random environments’, SIAM journal on control and optimization 49(2), 859–889.
  * Song & Zhu (2016) Song, Q. & Zhu, C. (2016), ‘On singular control problems with state constraints and regime-switching: A viscosity solution approach’, Automatica 70, 66–73.
  * Sylvia (1994) Sylvia, G. (1994), ‘Market information and fisheries management: A multiple-objective analysis’, North American Journal of Fisheries Management 14(2), 278–290.
  * Tran & Yin (2015) Tran, K. & Yin, G. (2015), ‘Optimal harvesting strategies for stochastic competitive Lotka–Volterra ecosystems’, Automatica 55, 236–246.
  * Tran & Yin (2017) Tran, K. & Yin, G. (2017), ‘Optimal harvesting strategies for stochastic ecosystems’, IET Control Theory & Applications 11(15), 2521–2530.
  * White & Hastings (2020) White, E. R. & Hastings, A. (2020), ‘Seasonality in ecology: Progress and prospects in theory’, Ecological Complexity 44, 100867.
  * Winsor (1932) Winsor, C. P. (1932), ‘The gompertz curve as a growth curve’, Proceedings of the National Academy of Sciences of the United States of America 18(1), 1.
  * Xu et al. (2005) Xu, C., Boyce, M. S. & Daley, D. J. (2005), ‘Harvesting in seasonal environments’, Journal of Mathematical Biology 50(6), 663–682.
  * Yin & Zhu (2009) Yin, G. G. & Zhu, C. (2009), Hybrid switching diffusions: properties and applications, Vol. 63, Springer Science & Business Media.
  * Yin et al. (2003) Yin, G., Zhang, Q. & Badowski, G. (2003), ‘Discrete-time singularly perturbed markov chains: Aggregation, occupation measures, and switching diffusion limit’, Advances in Applied Probability 35, 449–476.
  * Zeide (1993) Zeide, B. (1993), ‘Analysis of growth equations’, Forest science 39(3), 594–616.
  * Zhu (2011) Zhu, C. (2011), ‘Optimal control of the risk process in a regime-switching environment’, Automatica 47, 1570–1579.
  * Zhu & Yin (2009) Zhu, C. & Yin, G. (2009), ‘On competitive lotka–volterra model in random environments’, Journal of Mathematical Analysis and Applications 357(1), 154–170.

##  Appendix A Transition Probabilities

###  A.1. The formulation from Section 2

We first look at the details we need for the setting from Section 2. With the
notation defined in Section 2.1, let $(x,\alpha,u)\in
S_{h}\times{\mathcal{M}}\times\mathcal{U}$ and denote by
${\mathbb{E}}^{h,u}_{x,\alpha,n}$, ${\mathbb{C}ov}^{h,u}_{x,\alpha,n}$ the
conditional expectation and covariance given by

$\\{X_{m}^{h},\alpha^{h}_{m},U_{m}^{h},m\leq
n,X_{n}^{h}=x,\alpha^{h}_{n}=\alpha,U^{h}_{n}=u\\},$

respectively. Define $\Delta X^{h}_{n}=X^{h}_{n+1}-X^{h}_{n}$. Our objective
in this subsection is to define transition probabilities
$q^{h}((x,k),(y,l)|u)$ so that the controlled Markov chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ is locally consistent with respect to the
controlled diffusion (2.3). By this we mean that the following conditions
hold:

(A.1) $\begin{array}[]{ll}&\displaystyle{\mathbb{E}}^{h,u}_{x,k,n}\Delta
X_{n}^{h}=({b}(x,k)-u)\Delta t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u)),\\\
&\displaystyle{\mathbb{V}ar}^{h,u}_{x,k,n}\Delta
X_{n}^{h}=\sigma^{2}(x,k,u)\Delta t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u)),\\\
&\displaystyle{\mathbb{P}}^{h,u}_{x,k,n}(\alpha^{h}_{n+1}=l)=q_{kl}\Delta
t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u))\quad\text{for}\quad l\neq k,\\\
&\displaystyle{\mathbb{P}}^{h,u}_{x,k,n}(\alpha^{h}_{n+1}=k)=1+q_{kk}\Delta
t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u)),\\\ &\displaystyle\sup\limits_{n,\
\omega}|\Delta X_{n}^{h}|\to 0\quad\text{as}\quad h\to 0.\end{array}$

Using the procedure developed by Kushner (1990), for $(x,\alpha)\in
S_{h}\times{\mathcal{M}}$ and $u\in\mathcal{U}$, we define

(A.2) $\begin{array}[]{ll}&\displaystyle
Q_{h}(x,k,u)=\sigma^{2}(x,k)+h|b(x,k)-u|-h^{2}q_{kk}+h,\\\ &\displaystyle
q^{h}\left((x,k),(x+h,k)|u\right)=\dfrac{\sigma^{2}(x,k)/2+\big{(}b(x,k)-u\big{)}^{+}h}{Q_{h}(x,k,u)},\\\
&\displaystyle
q^{h}\left((x,k),(x-h,k)|u\right)=\dfrac{\sigma^{2}(x,k)/2+\left(b(x,k)-u\right)^{-}h}{Q_{h}(x,k,u)},\\\
&\displaystyle
q^{h}\left((x,k),(x,l)|u\right)=\dfrac{h^{2}q_{kl}}{Q_{h}(x,k,u)}\quad\text{
for }\quad k\neq l,\\\ &\displaystyle
q^{h}\left((x,k),(x,k)|u\right)=\dfrac{h}{Q_{h}(x,k,u)},\quad\Delta
t^{h}(x,k,u)=\dfrac{h^{2}}{Q_{h}(x,k,u)},\end{array}$

where for a real number $r$, $r^{+}=\max\\{r,0\\}$, $r^{-}=-\min\\{0,r\\}$.
Set $q^{h}\left((x,k),(y,l)|u\right)=0$ for all unlisted values of $(y,l)\in
S_{h}\times{\mathcal{M}}$. Note that $\sup_{x,k,u}\Delta t^{h}(x,k,u)\to 0$ as
$h\to 0$. Using the above transition probabilities, we can check that the
locally consistent conditions of $\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ are
satisfied.

######  Lemma A.1.

The Markov chain $\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ with transition
probabilities $\\{q^{h}(\cdot)\\}$ defined in (A.2) satisfies the local
consistence in (A.1).

###  A.2. Variable effort harvesting-stocking strategies

For $(x,\alpha,u)\in S_{h}\times{\mathcal{M}}\times\mathcal{U}$, let
${\mathbb{E}}^{h,u}_{x,\alpha,n}$, ${\mathbb{C}ov}^{h,u}_{x,\alpha,n}$ denote
the conditional expectation and covariance given by

$\\{X_{m}^{h},\alpha^{h}_{m},U_{m}^{h},m\leq
n,X_{n}^{h}=x,\alpha^{h}_{n}=\alpha,U^{h}_{n}=u\\},$

respectively. Define $\Delta X^{h}_{n}=X^{h}_{n+1}-X^{h}_{n}$. In order to
approximate the process $(X(\cdot),\alpha(\cdot))$ given in (3.1), the
controlled Markov chain $\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ must be locally
consistent with respect to $(X(\cdot),\alpha(\cdot))$ in the sense that the
following conditions hold

(A.3) $\begin{array}[]{ll}&\displaystyle{\mathbb{E}}^{h,u}_{x,k,n}\Delta
X_{n}^{h}=({b}(x,k)-ux)\Delta t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u)),\\\
&\displaystyle{\mathbb{V}ar}^{h,u}_{x,k,n}\Delta
X_{n}^{h}=\sigma^{2}(x,k,u)\Delta t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u)),\\\
&\displaystyle{\mathbb{P}}^{h,u}_{x,k,n}(\alpha^{h}_{n+1}=l)=q_{kl}\Delta
t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u))\quad\text{for}\,l\neq k,\\\
&\displaystyle{\mathbb{P}}^{h,u}_{x,k,n}(\alpha^{h}_{n+1}=k)=1+q_{kk}\Delta
t^{h}(x,k,u)+o(\Delta t^{h}(x,k,u)),\\\ &\displaystyle\sup\limits_{n,\
\omega}|\Delta X_{n}^{h}|\to 0\quad\text{as}\quad h\to 0.\end{array}$

To this end, we define the transition probabilities $q^{h}((x,k),(y,l)|u)$ as
follows. For $(x,k)\in S_{h}\times{\mathcal{M}}$ and $u\in\mathcal{U}$, let

(A.4) $\begin{array}[]{ll}&\displaystyle
Q_{h}(x,k,u)=\sigma^{2}(x,k)+h|b(x,k)-ux|-h^{2}q_{kk}+h,\\\ &\displaystyle
q^{h}\left((x,k),(x+h,k)|u\right)=\dfrac{\sigma^{2}(x,k)/2+\big{(}b(x,k)-ux)^{+}h}{Q_{h}(x,k,u)},\\\
&\displaystyle
q^{h}\left((x,k),(x-h,k)|u\right)=\dfrac{\sigma^{2}(x,k)/2+\left(b(x,k)-ux\right)^{-}h}{Q_{h}(x,k,u)},\\\
&\displaystyle
q^{h}\left((x,k),(x,l)|u\right)=\dfrac{h^{2}q_{kl}}{Q_{h}(x,k,u)}\quad\text{
for }\quad k\neq l,\\\ &\displaystyle
q^{h}\left((x,k),(x,k)|u\right)=\dfrac{h}{Q_{h}(x,k,u)},\quad\Delta
t^{h}(x,k,u)=\dfrac{h^{2}}{Q_{h}(x,k,u)}.\end{array}$

Set $q^{h}\left((x,k),(y,l)|u\right)=0$ for all unlisted values of $(y,l)\in
S_{h}\times{\mathcal{M}}$.

###  A.3. Uncertain price functions

Let $(\phi,x,\alpha,u)\in\widehat{S}_{h}\times{\mathcal{M}}\times\mathcal{U}$
and denote by ${\mathbb{E}}^{h,u}_{\phi,x,\alpha,n}$,
${\mathbb{C}ov}^{h,u}_{\phi,x,\alpha,n}$ the conditional expectation and
covariance given by

$\\{\Phi^{h}_{m},X_{m}^{h},\alpha^{h}_{m},U_{m}^{h},m\leq
n,\Phi^{h}_{n}=\phi,X_{n}^{h}=x,\alpha^{h}_{n}=\alpha,U^{h}_{n}=u\\},$

respectively. Define $\Delta X^{h}_{n}=X^{h}_{n+1}-X^{h}_{n}$ and
$\Delta\Phi^{h}_{n}=\Phi^{h}_{n+1}-\Phi^{h}_{n}$. In order to approximate
$(\Phi(\cdot),X(\cdot),\alpha(\cdot))$ given by (2.3)-and-(3.2), the
controlled Markov chain $\\{(\Phi^{h}_{n},X^{h}_{n},\alpha^{h}_{n})\\}$ must
be locally consistent with respect to $(\Phi(\cdot),X(\cdot),\alpha(\cdot))$
in the sense that the following conditions hold:

(A.5) $\begin{array}[]{ll}&\displaystyle{\mathbb{E}}^{h,u}_{\phi,x,k,n}\Delta
X_{n}^{h}=({b}(x,k)-u)\Delta t^{h}(\phi,x,k,u)+o(\Delta t^{h}(\phi,x,k,u)),\\\
&\displaystyle{\mathbb{V}ar}^{h,u}_{\phi,x,k,n}\Delta
X_{n}^{h}=\sigma^{2}(x,k,u)\Delta t^{h}(\phi,x,k,u)+o(\Delta
t^{h}(\phi,x,k,u)),\\\
&\displaystyle{\mathbb{E}}^{h,u}_{\phi,x,k,n}\Delta\Phi_{n}^{h}={b}_{0}(x,k)\Delta
t^{h}(\phi,x,k,u)+o(\Delta t^{h}(\phi,x,k,u)),\\\
&\displaystyle{\mathbb{V}ar}^{h,u}_{\phi,x,k,n}\Delta\Phi_{n}^{h}=\sigma^{2}_{0}(x,k)\Delta
t^{h}(\phi,x,k,u)+o(\Delta t^{h}(\phi,x,k,u)),\\\
&\displaystyle{\mathbb{P}}^{h,u}_{\phi,x,k,n}(\alpha^{h}_{n+1}=l)=q_{kl}\Delta
t^{h}(\phi,x,k,u)+o(\Delta t^{h}(\phi,x,k,u))\quad\text{for}\quad l\neq k,\\\
&\displaystyle{\mathbb{P}}^{h,u}_{\phi,x,k,n}(\alpha^{h}_{n+1}=k)=1+q_{kk}\Delta
t^{h}(\phi,x,k,u)+o(\Delta t^{h}(\phi,x,k,u)),\\\
&\displaystyle\sup\limits_{n,\ \omega}\big{(}|\Delta
X_{n}^{h}|+|\Delta\Phi^{h}_{n}|\big{)}\to 0\quad\text{as}\quad h\to
0.\end{array}$

To this end, we define the transition probabilities
$q^{h}((\phi,x,k),(\psi,y,l)|u)$ as follows. For
$(\phi,x,k)\in\widehat{S}_{h}\times{\mathcal{M}}$ and $u\in\mathcal{U}$, let

(A.6) $\begin{array}[]{ll}&\displaystyle
Q_{h}(\phi,x,k,u)=\sigma^{2}(x,k)+h|b(x,k)-u|+\sigma_{0}^{2}(x,k)+h|b_{0}(x,k)|-h^{2}q_{kk}+h,\\\
&\displaystyle
q^{h}\left((\phi,x,k),(\phi,x+h,k)|u\right)=\dfrac{\sigma^{2}(x,k)/2+\big{(}b(x,k)-u\big{)}^{+}h}{Q_{h}(\phi,x,k,u)},\\\
&\displaystyle
q^{h}\left((\phi,x,k),(\phi,x-h,k)|u\right)=\dfrac{\sigma^{2}(x,k)/2+\left(b(x,k)-u\right)^{-}h}{Q_{h}(\phi,x,k,u)},\\\
&\displaystyle
q^{h}\left((\phi,x,k),(\phi,x,l)|u\right)=\dfrac{h^{2}q_{kl}}{Q_{h}(\phi,x,k,u)}\quad\text{for}\quad
k\neq l,\\\ &\displaystyle
q^{h}\left((\phi,x,k),(\phi,x,k)|u\right)=\dfrac{h}{Q_{h}(\phi,x,k,u)},\quad\Delta
t^{h}(\phi,x,k,u)=\dfrac{h^{2}}{Q_{h}(\phi,x,k,u)},\\\ &\displaystyle
q^{h}\left((\phi,x,k),(\phi+h,x,k)|u\right)=\dfrac{\sigma^{2}_{0}(x,k)/2+hb^{+}_{0}(x,k)}{Q_{h}(\phi,x,k,u)},\\\
&\displaystyle
q^{h}\left((\phi,x,k),(\phi-h,x,k)|u\right)=\dfrac{\sigma_{0}^{2}(x,k)/2+hb^{-}_{0}(x,k)}{Q_{h}(\phi,x,k,u)}.\end{array}$

Set $q^{h}\left((\phi,x,k),(\psi,y,l)|u\right)=0$ for all unlisted values of
$(\psi,y,l)\in\widehat{S}_{h}\times{\mathcal{M}}$.

###  A.4. The combined effects of seasonality and Markovian switching

Recall that
$\widetilde{S}_{h}:=\\{(\gamma,x)=(k_{1}h_{1},k_{2}h_{2})^{\prime}\in\mathbb{R}^{2}:k_{i}\in\mathbb{Z}_{\geq
0},k_{1}\leq T/h_{1}\\}.$ Let
$(\gamma,x,\alpha,u)\in\widetilde{S}_{h}\times{\mathcal{M}}\times\mathcal{U}$
and denote by ${\mathbb{E}}^{h,u}_{\gamma,x,\alpha,n}$,
${\mathbb{C}ov}^{h,u}_{\gamma,x,\alpha,n}$ the conditional expectation and
covariance given by

$\\{\Gamma^{h}_{m},X_{m}^{h},\alpha^{h}_{m},U_{m}^{h},m\leq
n,\Gamma^{h}_{n}=\gamma,X_{n}^{h}=x,\alpha^{h}_{n}=\alpha,U^{h}_{n}=u\\},$

respectively. Define $\Delta X^{h}_{n}=X^{h}_{n+1}-X^{h}_{n}$ and
$\Delta\Gamma^{h}_{n}=\Gamma^{h}_{n+1}-\Gamma^{h}_{n}$. In order to
approximate $(\Gamma(\cdot),X(\cdot),\alpha(\cdot))$ given by (3.4), the
controlled Markov chain $\\{(\Gamma^{h}_{n},X^{h}_{n},\alpha^{h}_{n})\\}$ must
be locally consistent with respect to $(\Gamma(\cdot),X(\cdot),\alpha(\cdot))$
in the sense that the following conditions hold:

(A.7)
$\begin{array}[]{ll}&\displaystyle{\mathbb{E}}^{h,u}_{\gamma,x,k,n}\Delta
X_{n}^{h}=({b}(\gamma,x,k)-u)\Delta t^{h}(\gamma,x,k,u)+o(\Delta
t^{h}(\gamma,x,k,u)),\\\
&\displaystyle{\mathbb{V}ar}^{h,u}_{\gamma,x,k,n}\Delta
X_{n}^{h}=\sigma^{2}(\gamma,x,k,u)\Delta t^{h}(\gamma,x,k,u)+o(\Delta
t^{h}(\gamma,x,k,u)),\\\ &\displaystyle\Delta\Gamma_{n}^{h}=\Delta
t^{h}(\gamma,x,k,u),\quad\text{if}\quad\gamma+\Delta t^{h}(\phi,x,k,u)<T,\\\
&\displaystyle\Delta\Gamma_{n}^{h}=\gamma+\Delta
t^{h}(\gamma,x,k,u)-T\quad\text{if}\quad\gamma+\Delta t^{h}(\gamma,x,k,u)\geq
T,\\\
&\displaystyle{\mathbb{P}}^{h,u}_{\gamma,x,k,n}(\alpha^{h}_{n+1}=l)=q_{kl}\Delta
t^{h}(\gamma,x,k,u)+o(\Delta t^{h}(\gamma,x,k,u))\quad\text{for }\,l\neq k,\\\
&\displaystyle{\mathbb{P}}^{h,u}_{\gamma,x,k,n}(\alpha^{h}_{n+1}=k)=1+q_{kk}\Delta
t^{h}(\gamma,x,k,u)+o(\Delta t^{h}(\gamma,x,k,u)),\\\
&\displaystyle\sup\limits_{n,\ \omega}\big{(}|\Delta
X_{n}^{h}|+\Delta\Gamma^{h}_{n}\big{)}\to 0\quad\text{as}\quad h\to
0.\end{array}$

To this end, we define the transition probabilities
$q^{h}((\gamma,x,k),(\lambda,y,l)|u)$ as follows. Let
$(\gamma,x,k)\in\widetilde{S}_{h}\times{\mathcal{M}}$ and $u\in\mathcal{U}$.
If $\gamma+h_{1}=T$, $\gamma+h_{1}$ in the following definition is understood
as $0$. Let

(A.8) $\begin{array}[]{ll}&\\!\\!\\!\displaystyle\Delta
t^{h}(\gamma,x,k,u)=h_{1},\\\ &\\!\\!\\!\displaystyle
q^{h}\left((\gamma,x,k),(\gamma+h_{1},x+h_{2},k)|u\right)=\dfrac{\Big{(}\sigma^{2}(\gamma,x,k)/2+\big{(}b(\gamma,x,k)-u\big{)}^{+}h_{2}\Big{)}h_{1}}{h_{2}^{2}},\\\
&\\!\\!\\!\displaystyle
q^{h}\left((\gamma,x,k),(\gamma+h_{1},x-h_{2},k)|u\right)=\dfrac{\Big{(}\sigma^{2}(\gamma,x,k)/2+\big{(}b(\gamma,x,k)-u\big{)}^{-}h_{2}\Big{)}h_{1}}{h_{2}^{2}},\\\
&\\!\\!\\!\displaystyle
q^{h}\big{(}(\gamma,x,k),(\gamma+h_{1},x,l)|u\big{)}=h_{1}q_{kl},\quad\text{for}\quad
l\neq k,\\\ &\\!\\!\\!\displaystyle
q^{h}\big{(}(\gamma,x,k),(\gamma+h_{1},x,k)|u\big{)}=1-q^{h}\big{(}(\gamma,x,k),(\gamma+h_{1},x+h_{2},k)|u\big{)}\\\
&\\!\\!\\!\displaystyle\qquad-q^{h}\big{(}(\gamma,x,k),(\gamma+h_{1},x-h_{2},k)|u\big{)}-\sum\limits_{l\neq
k}q^{h}\big{(}(\gamma,x,k),(\gamma+h_{1},x,l)|u\big{)}.\end{array}$

Set $q^{h}\left((\gamma,x,k),(\lambda,y,l)|u\right)=0$ for all unlisted values
of $(\lambda,y,l)\in\widetilde{S}_{h}\times{\mathcal{M}}$.

##  Appendix B Continuous–time interpolation

We will present the convergence analysis for the formulation in Section 2. The
other formulas can be handled in a similar way. Our procedure and methods are
similar to those in Kushner (1990), Kushner & Dupuis (1992), Song et al.
(2006). The convergence result is based on a continuous-time interpolation of
the controlled Markov chain, which will be constructed to be piecewise
constant on the time interval $[t^{h}_{n},t^{h}_{n+1}),n\geq 0$. To this end,
we define $n^{h}(t)=\max\\{n:t^{h}_{n}\leq t\\},t\geq 0.$ The piecewise
constant interpolation of $\\{(X^{h}_{n},\alpha^{h}_{n},U^{h}_{n})\\}$,
denoted by $\big{(}X^{h}(t),\alpha^{h}(t),U^{h}(t)\big{)}$ is naturally
defined as

(B.1) $\begin{array}[]{ll}&\displaystyle
X^{h}(t)=X^{h}_{n^{h}(t)},\quad\alpha^{h}(t)=\alpha^{h}_{n^{h}(t)},\quad
U^{h}(t)=U^{h}_{n^{h}(t)},\quad t\geq 0.\end{array}$

Define $\mathcal{F}^{h}(t)=\sigma\\{X^{h}(s),\alpha^{h}(s),U^{h}(s):s\leq
t\\}=\mathcal{F}^{h}_{n^{h}(t)}$. Also define

$M^{h}(0)=0,\quad M^{h}(t)=\sum\limits_{m=0}^{n^{h}(t)-1}(\Delta
X^{h}_{m}-{\mathbb{E}}^{h}_{m}\Delta X_{m})\quad\text{for}\quad t\geq 0.$

It is obvious that

(B.2) $X^{h}(t)=x+\sum\limits_{m=0}^{n^{h}(t)-1}{\mathbb{E}}^{h}_{m}\Delta
X^{h}_{m}+M^{h}(t).$

Recall that $\Delta
t^{h}_{m}=h^{2}/Q_{h}(X^{h}_{m},\alpha^{h}_{m},U^{h}_{m})$. It follows that

(B.3)
$\begin{array}[]{ll}\sum\limits_{m=0}^{n^{h}(t)-1}{\mathbb{E}}^{h}_{m}\Delta
X^{h}_{m}&\\!\\!\\!\displaystyle=\sum\limits_{m=0}^{n^{h}(t)-1}\left[b(X^{h}_{m},\alpha^{h}_{m})+U^{h}_{m}\right]\Delta
t^{h}_{m}\\\
&\\!\\!\\!\displaystyle=\int_{0}^{t}\left[b(X^{h}(s),\alpha^{h}(s))+U^{h}(s)\right]ds-\int_{t^{h}_{n^{h}(t)}}^{t}\left[b(X^{h}(s),\alpha^{h}(s))+U^{h}(s)\right]ds\\\
&\\!\\!\\!\displaystyle=\int_{0}^{t}\left[b(X^{h}(s),\alpha^{h}(s))+U^{h}(s)\right]ds+\varepsilon^{h}_{1}(t),\end{array}$

with $\\{\varepsilon_{1}^{h}(\cdot)\\}$ being an $\mathcal{F}^{h}(t)$-adapted
process satisfying

$\lim\limits_{h\to
0}\sup\limits_{t\in[0,T_{0}]}{\mathbb{E}}|\varepsilon_{1}^{h}(t)|=0\quad\text{for
any }0<T_{0}<\infty.$

For simplicity, we suppose that $\inf\limits_{(x,k)}1/|\sigma(x,k)|>0$ (if
this is not the case, we can use the trick from (Kushner & Dupuis 1992,
p.288-289)). Define $w^{h}(\cdot)$ by

(B.4)
$w^{h}(t)=\sum\limits_{m=0}^{n^{h}(t)-1}\big{[}1/\sigma(X^{h}_{m},\alpha^{h}_{m})\big{]}(\Delta
X^{h}_{m}-{\mathbb{E}}^{h}_{m}\Delta X^{h}_{m}).$

Then we can write

(B.5)
$M^{h}(t)=\int_{0}^{t}\sigma(X^{h}(s),\alpha^{h}(s))dw^{h}(s)+\varepsilon_{2}^{h}(t),$

with $\\{\varepsilon_{2}^{h}(\cdot)\\}$ being an $\mathcal{F}^{h}(t)$-adapted
process satisfying

$\lim\limits_{h\to
0}\sup\limits_{t\in[0,T_{0}]}{\mathbb{E}}|\varepsilon_{2}^{h}(t)|=0\quad\text{for
any }0<T_{0}<\infty.$

Using (B.3) and (B.5), we can write (B.2) as

(B.6)
$X^{h}(t)=x+\int_{0}^{t}\left[b(X^{h}(s),\alpha^{h}(s))+U^{h}(s)\right]ds+\int_{0}^{t}\sigma(X^{h}(s),\alpha^{h}(s))dw^{h}(s)+\varepsilon^{h}(t),$

where $\varepsilon^{h}(\cdot)$ is an $\mathcal{F}^{h}(t)$-adapted process
satisfying

$\lim\limits_{h\to
0}\sup\limits_{t\in[0,T_{0}]}{\mathbb{E}}|\varepsilon^{h}(t)|=0\quad\text{for
any }0<T_{0}<\infty.$

The performance function from (2.9) can be rewritten as

(B.7) $J^{h}(x,\alpha,U^{h}(\cdot))={\mathbb{E}}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X^{h}(s),\alpha^{h}(s),U^{h}(s)\big{)}ds.$

##  Appendix C Convergence

The convergence of the algorithms is established via the weak convergence
method. To proceed, let $D[0,\infty)$ denote the space of functions that are
right continuous and have left-hand limits endowed with the Skorokhod
topology. All the weak analysis will be on this space or its $k$-fold products
$D^{k}[0,\infty)$ for appropriate $k$. We follow (Kushner & Dupuis 1992,
Section 4.6) in order to introduce relaxed control representations, which we
need in order to prove the weak convergence.

######  Definition C.1.

Let $\mathcal{B}(\mathcal{U}\times[0,\infty))$ be the $\sigma$-algebra of
Borel subsets of $\mathcal{U}\times[0,\infty)$. An admissible relaxed control,
which we will call a relaxed control, $m(\cdot)$ is a measure on
$\mathcal{B}(\mathcal{U}\times[0,\infty))$ such that

$m(\mathcal{U}\times[0,t])=t\quad\text{for all}\quad t\geq 0.$

Given a relaxed control $m(\cdot)$, there is a probability measure
$m_{t}(\cdot)$ defined on the $\sigma$-algebra $\mathcal{B}(\mathcal{U})$ such
that $m(dudt)=m_{t}(du)dt$. Let $\mathcal{R}(\mathcal{U}\times[0,\infty))$
denote the set of all relaxed controls on $\mathcal{U}\times[0,\infty)$.

With the given probability space, we say that $m(\cdot)$ is an admissible
relaxed (stochastic) control if (i) for each fixed $t\geq 0$, $m(t,\cdot)$ is
a random variable taking values in $\mathcal{R}(\mathcal{U}\times[0,\infty))$,
and for each fixed $\omega$, $m(\cdot,\omega)$ is a deterministic relaxed
control; (ii) the function defined by $m(A\times[0,t])$ is
$\mathcal{F}(t)$-adapted for any $A\in\mathcal{B}(\mathcal{U})$. As a result,
with probability one, there is a measure $m_{t}(\cdot,\omega)$ on the Borel
$\sigma$-algebra $\mathcal{B}(\mathcal{U})$ such that $m(dcdt)=m_{t}(dc)dt$.

######  Remark C.2.

For a sequence of controls $U^{h}=\\{U^{h}_{n}:n\in\mathbb{Z}_{\geq 0}\\}$, we
define a sequence of relaxed control equivalence as follows. First, we set
$m_{t}^{h}(du)=\delta_{U^{h}(t)}(du)$ for $t\geq 0$, where
$\delta_{U^{h}(t)}(\cdot)$ is the probability measure concentrated at
$U^{h}(t)$. Then $m^{h}(\cdot)$ is defined by $m^{h}(dudt)=m_{t}(du)dt$; that
is,

$m^{h}(B\times[0,t])=\int_{0}^{t}\Big{(}\int_{B}\delta_{U^{h}(s)}(du)\Big{)}ds,\quad
B\in\mathcal{B}({\mathcal{U}})\quad\text{and}\quad t\geq 0.$

Recall that $\mathcal{R}(\mathcal{U}\times[0,\infty))$ is the space of all
relaxed controls on $\mathcal{U}\times[0,\infty)$. Then
$\mathcal{R}(\mathcal{U}\times[0,\infty))$ can be metrized using the Prokhorov
metric in the usual way as in (Kushner & Dupuis 1992, pages 263–264). With the
Prokhorov metric, $\mathcal{R}(\mathcal{U}\times[0,\infty))$ is a compact
space. It follows that any sequence of relaxed controls has a convergent
subsequence. Moreover, a sequence $(\eta_{n})_{n\in\mathbb{N}}$ with
$\eta_{n}\in\mathcal{R}(\mathcal{U}\times[0,\infty))$ converges to
$\eta\in\mathcal{R}(\mathcal{U}\times[0,\infty))$ if and only if for any
continuous functions with compact support $\Psi(\cdot)$ on
$\mathcal{U}\times[0,\infty)$ one has

$\int_{\mathcal{U}\times[0,\infty)}\Psi(u,s)\eta_{n}(du,ds)\to\int_{\mathcal{U}\times[0,\infty)}\Psi(u,s)\eta(du,ds)$

as $n\to\infty$. Note that for a sequence of ordinary controls
$U^{h}=\\{U^{h}_{n}:n\in\mathbb{Z}_{\geq 0}\\}$, the associated relaxed
control ${m}^{h}(dcdt)$ belongs to $\mathcal{R}(\mathcal{U}\times[0,\infty))$.
Note also that the limits of the “relaxed control representations” of the
ordinary controls might not be ordinary controls, but only relaxed controls.

With the notion of relaxed control given above, we can write (B.6) and (B.7)
as

(C.1)
$X^{h}(t)=x+\int_{0}^{t}\left[b(X^{h}(s),\alpha^{h}(s))+U(s)\right]m^{h}_{s}(du)ds+\int_{0}^{t}\sigma(X^{h}(s),\alpha^{h}(s))dw^{h}(s)+\varepsilon^{h}(t),$
(C.2) $J^{h}(x,\alpha,m^{h}(\cdot))={\mathbb{E}}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X^{h}(s),\alpha^{h}(s),u\big{)}m_{s}^{h}(du)ds.$

The value function defined in (2.6) can be rewritten as

$V(x,\alpha)=\sup\\{J(x,\alpha,m(\cdot)):m(\cdot)\,\,\,\text{ is an admissible
relaxed control}\\},$

where

$J(x,\alpha,m(\cdot)):={\mathbb{E}}_{x,\alpha}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X(s),\alpha(s),u\big{)}m_{s}(du)ds.$

######  Lemma C.3.

The process $\\{\alpha^{h}(\cdot)\\}$ converges weakly to $\alpha(\cdot)$,
which is a Markov chain with generator $Q=(q_{kl})$.

###### Proof.

The proof is similar to that of (Yin et al. 2003, Theorem 3.1) and is
therefore omitted. $\qquad\Box$

######  Theorem C.4.

Suppose Assumption 2.1 holds. Let the chain $\\{(X^{h}_{n},\alpha^{h}_{n})\\}$
be constructed using the transition probabilities defined in (A.2),
$\big{(}X^{h}(\cdot),\alpha^{h}(\cdot),w^{h}(\cdot)\big{)}$ be the continuous-
time interpolation defined in (B.1) and (B.4), $\\{U^{h}_{n}\\}$ be an
admissible strategy and $m^{h}(\cdot)$ be the relaxed control representation
of $\\{U^{h}_{n}\\}$. Then the following assertions hold.

  * (a)

The family of processes
$H^{h}(\cdot)=\big{(}{X}^{h}(\cdot),{\alpha}^{h}(\cdot),m^{h}(\cdot),{w}^{h}(\cdot)\big{)}$
is tight. As a result, it has a weakly convergent subsequence with limit
$H(\cdot)=\big{(}{X}(\cdot),{\alpha}(\cdot),m(\cdot),{w}(\cdot)\big{)}.$

  * (b)

Let $\mathcal{F}(t)$ be the $\sigma$-algebra generated by $\big{\\{}H(s):s\leq
t\big{\\}}$. Then $w(\cdot)$ is a standard $\mathcal{F}(t)$ adapted Brownian
motion, $m(\cdot)$ is an admissible control, and

(C.3)
$X(t)=x+\int_{0}^{t}\left[b(X(s),\alpha(s))+u\right]m_{s}(du)ds+\int_{0}^{t}\sigma(X(s),\alpha(s))dw(s),\quad
t\geq 0.$

###### Proof.

(a) We use the tightness criteria in (Kushner 1984, p. 47). Specifically, a
sufficient condition for tightness of a sequence of processes
$\zeta^{h}(\cdot)$ with paths in $D^{k}[0,\infty)$ is that for any
$T_{0},\rho\in(0,\infty)$,

$\begin{array}[]{ll}&\\!\\!\\!\displaystyle{\mathbb{E}}_{t}^{h}\big{|}\zeta^{h}(t+s)-\zeta^{h}(t)\big{|}^{2}\leq{\mathbb{E}}^{h}_{t}\gamma(h,\rho)\quad\text{for
all}\quad s\in[0,\rho],\quad t\leq T_{0},\\\
&\\!\\!\\!\displaystyle\lim\limits_{\rho\to 0}\limsup\limits_{h\to
0}{\mathbb{E}}\gamma(h,\rho)=0.\end{array}$

The tightness of $\\{\alpha^{h}(\cdot)\\}$ is obvious by the preceding lemma.
The process $\\{m^{h}(\cdot)\\}$ is tight since its range space is compact. It
is standard to show the tightness of $\\{w^{h}(\cdot)\\}$ and
$\\{X^{h}(\cdot)\\}$ – see Song et al. (2006) for details. As a result, a
subsequence of
$H^{h}(\cdot)=\big{(}X^{h}(\cdot),\alpha^{h}(\cdot),m^{h}(\cdot),w^{h}(\cdot)\big{)}$
converges weakly to the limit
$H(\cdot)=\big{(}X(\cdot),\alpha(\cdot),m(\cdot),w(\cdot)\big{)}$.

(b) For the rest of the proof, we assume the probability space is chosen as
required by Skorokhod representation. Thus, with a slight abuse of notation,
we assume that $H^{h}(\cdot)$ converges to the limit $H(\cdot)$ with
probability one via Skorokhod representation.

To characterize ${w}(\cdot)$, let $\widetilde{k},\widetilde{j}$ be arbitrary
positive integers. Pick $t>0$, $\rho>0$ and $\\{t_{k}:k\leq\widetilde{k}\\}$
such that $t_{k}\leq t\leq t+\rho$ for each $k$. Let $\phi_{j}(\cdot)$ be
real-valued continuous functions that are compactly supported on
$\mathcal{U}\times[0,\infty)$ for any $j\leq\widetilde{j}$. Define
$(\phi_{j},m)_{t}:=\int_{0}^{t}\int_{\mathcal{U}}\phi_{j}(u,s)m(duds)$.

Let $\Psi(\cdot)$ be a real-valued and continuous function of its arguments
with compact support. By the definition of $w^{h}(\cdot)$ in (B.4),
$w^{h}(\cdot)$ is an $\mathcal{F}^{h}(t)$-martingale. Thus, we have

(C.4)
${\mathbb{E}}\Psi\big{(}X^{h}(t_{k}),\alpha^{h}(t_{k}),w^{h}(t_{k}),(\phi_{j},m^{h})_{t_{k}},j\leq\widetilde{j},k\leq\widetilde{k}\big{)}\big{[}{w}^{h}(t+\rho)-{w}^{h}(t)\big{]}=0,$

and

(C.5)
${\mathbb{E}}\Psi\big{(}X^{h}(t_{k}),\alpha^{h}(t_{k}),w^{h}(t_{k}),(\phi_{j},m^{h})_{t_{k}},j\leq\widetilde{j},k\leq\widetilde{k}\big{)}\big{[}\big{(}{w}^{h}(t+\rho)\big{)}^{2}-\big{(}{w}^{h}(t)\big{)}^{2}-\rho-{\varepsilon}^{h}(\rho)\big{]}=0.$

By using the Skorokhod representation and the dominated convergence theorem,
letting $h\to 0$ in (C.4), we obtain

(C.6)
${\mathbb{E}}\Psi\big{(}X(t_{k}),\alpha(t_{k}),w(t_{k}),(\phi_{j},m)_{t_{k}},j\leq\widetilde{j},k\leq\widetilde{k}\big{)}\big{[}{w}(t+\rho)-{w}(t)\big{]}=0.$

Since ${w}(\cdot)$ has continuous paths with probability one, (C.6) implies
that ${w}(\cdot)$ is a continuous ${\mathcal{F}}(\cdot)$-martingale. Moreover,
(C.5) gives us that

(C.7)
${\mathbb{E}}\Psi\big{(}X(t_{k}),\alpha(t_{k}),w(t_{k}),(\phi_{j},m)_{t_{k}},j\leq\widetilde{j},k\leq\widetilde{k}\big{)}\big{[}\big{(}{w}(t+\rho)\big{)}^{2}-\big{(}{w}(t)\big{)}^{2}-\rho\big{]}=0.$

Thus, the quadratic variation of $w(t)$ is $t$, which implies that $w(\cdot)$
is a standard $\mathcal{F}(t)$ adapted Brownian motion.

By the convergence with probability one via Skorokhod representation, we have

${\mathbb{E}}\left|\int_{0}^{t}\int_{\mathcal{U}}\left[b(X^{h}(s),\alpha^{h}(s))+u\right]m_{s}^{h}(du)ds-\int_{0}^{t}\int_{\mathcal{U}}\left[b(X(s),\alpha(s))+u\right]m_{s}^{h}(du)ds\right|\to
0$

uniformly in $t$ as $h\to 0$.

Also, by the weak convergence of $\\{m^{h}(\cdot)\\}$, for any bounded and
continuous function $\phi(\cdot)$ with compact support,
$(\phi,m^{h})_{\infty}\to(\phi,m)_{\infty}$; see also Remark C.2. The weak
convergence and the Skorokhod representation imply that

$\int_{0}^{t}\int_{\mathcal{U}}\big{[}b(X(s),\alpha(s))+u\big{]}m_{s}^{h}(du)ds-\int_{0}^{t}\int_{\mathcal{U}}\big{[}b(X(s),\alpha(s))+u\big{]}m_{s}(du)ds\to
0$

uniformly in $t$ on any bounded interval with probability one.

For each positive constant $\rho$ and process ${\nu}(\cdot)$, define the
piecewise constant process ${\nu}^{\rho}(\cdot)$ by
${\nu}^{\rho}(t)={\nu}(k\rho)$ for
$t\in[k\rho,k\rho+\rho),k\in\mathbb{Z}_{\geq 0}$. Then, by the tightness of
$({X}^{h}(\cdot),{\alpha}^{h}(\cdot))$, (C.1) can be rewritten as

${X}^{h}(t)=x+\int_{0}^{t}\int_{\mathcal{U}}\big{[}b({X}^{h}(s),{\alpha}^{h}(s))+u\big{]}m_{s}^{h}(du)ds+\int_{0}^{t}\sigma({X}^{h,\rho}(s),{\alpha}^{h,\rho}(s))d{w}^{h}(s)+{\varepsilon}^{h,\rho}(t),$

where $\lim\limits_{\rho\to 0}\limsup\limits_{h\to
0}{\mathbb{E}}|{\varepsilon}^{h,\rho}(t)|=0.$ Noting that the processes
${X}^{h,\rho}(\cdot)$ and ${\alpha}^{h,\rho}(\cdot)$ take constant values on
the intervals $[n\rho,n\rho+\rho)$, we have

$\int_{0}^{t}\sigma({X}^{h,\rho}(s),{\alpha}^{h,\rho}(s))d{w}^{h}(s)\to\int_{0}^{t}\sigma({X}^{\rho}(s),{\alpha}^{\rho}(s))d{w}(s)\quad\text{
as }\quad h\to 0.$

The integrals above are well defined with probability one since they can be
written as finite sums. Combining the last results, we have

${X}(t)=x+\int_{0}^{t}\int_{\mathcal{U}}\big{[}b({X}(s),{\alpha}(s))+u\big{]}m_{s}(du)ds+\int_{0}^{t}\sigma({X}^{\rho}(s),{\alpha}^{\rho}(t))d{w}(s)+{\varepsilon}^{\rho}(t),$

where $\lim\limits_{\rho\to 0}E|{\varepsilon}^{\rho}(t)|=0.$ Taking the limit
as $\rho\to 0$ finishes the proof. $\qquad\Box$

######  Theorem C.5.

Suppose Assumption 2.1 holds. Let $V^{h}(x,\alpha)$ and $V(x,\alpha)$ be the
value functions defined in (2.6) and (2.9). Then $V^{h}(x,\alpha)\to
V(x,\alpha)$ as $h\to 0$.

###### Proof.

The proof is motivated by that of Theorem 7 in Song et al. (2006). Let
$U^{h}(\cdot)$ be an admissible strategy for the chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ and $m^{h}(\cdot)$ be the corresponding
relaxed control representation. Without loss of generality (passing to an
additional subsequence if needed), we assume that
$\big{(}{X}^{{h}}(\cdot),{\alpha}^{{h}}(\cdot),{w}^{{h}}(\cdot),{m}^{{h}}(\cdot)\big{)}$
converges weakly to
$\big{(}{X}(\cdot),{\alpha}(\cdot),{w}(\cdot),{m}(\cdot)\big{)}$. We show that
as $h\to 0$ we have

(C.8) $J^{h}(x,\alpha,U^{h}(\cdot))\to J(x,\alpha,m(\cdot)).$

From (C.2) one has

(C.9)
$\begin{array}[]{ll}J^{h}(x,\alpha,U^{h}(\cdot))&\\!\\!\\!\displaystyle={\mathbb{E}}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X^{h}(s),\alpha^{h}(s),u\big{)}m_{s}^{h}(du)ds.\end{array}$

By the weak convergence and the Skorokhod representation, as $h\to 0$,

$J^{h}(x,\alpha,U^{h}(\cdot))\to{\mathbb{E}}\int_{0}^{\infty}e^{-\delta
s}p\big{(}X(s),\alpha(s),u)\big{)}m_{s}(du)ds.$

This yields that $J^{h}(x,\alpha,U^{h}(\cdot))\to J(x,\alpha,m(\cdot))$ as
$h\to 0$.

Next, we prove that

(C.10) $\limsup\limits_{h\to 0}V^{h}(x,\alpha)\leq V(x,\alpha).$

For any small positive constant $\varepsilon$, let $\widetilde{U}^{h}(\cdot)$
be an $\varepsilon$-optimal harvesting strategy for the chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$; that is,

$V^{h}(x,\alpha)=\sup\limits_{U^{h}(\cdot)}J^{h}(x,\alpha,U^{h}(\cdot))\leq
J^{h}(x,\alpha,\widetilde{U}^{h}(\cdot))+\varepsilon.$

Choose a subsequence $\\{\widetilde{h}\\}$ of $\\{h\\}$ such that

(C.11) $\limsup\limits_{{h}\to
0}V^{{h}}(x,\alpha)=\lim\limits_{\widetilde{h}\to
0}V^{\widetilde{h}}(x,\alpha)\leq\limsup\limits_{\widetilde{h}\to
0}J^{\widetilde{h}}(x,\alpha,{\widetilde{U}}^{\widetilde{h}}(\cdot))+\varepsilon.$

Let $\widetilde{m}^{\widetilde{h}}(\cdot)$ be the relaxed control
representation of $\widetilde{U}^{\widetilde{h}}(\cdot)$. Without loss of
generality (passing to an additional subsequence if needed), we may assume
that
$\big{(}{X}^{\widetilde{h}}(\cdot),{\alpha}^{\widetilde{h}}(\cdot),{w}^{\widetilde{h}}(\cdot),{m}^{\widetilde{h}}(\cdot)\big{)}$
converges weakly to
$\big{(}{X}(\cdot),{\alpha}(\cdot),{w}(\cdot),{m}(\cdot)\big{)}$. It follows
from our claim in the beginning of the proof that

(C.12) $\lim\limits_{\widetilde{h}\to
0}J^{\widetilde{h}}(x,\alpha,{\widetilde{U}}^{\widetilde{h}}(\cdot))=\lim\limits_{\widetilde{h}\to
0}J^{\widetilde{h}}(x,\alpha,{\widetilde{m}}^{\widetilde{h}}(\cdot))=J(x,\alpha,m(\cdot))\leq
V(x,\alpha),$

where $J(x,\alpha,m(\cdot))\leq V(x,\alpha)$ by the definition of
$V(x,\alpha)$. Since $\varepsilon$ is arbitrarily small, (C.10) follows from
(C.11) and (C.12).

To prove the reverse inequality $\liminf\limits_{h}V^{h}(x,\alpha)\geq
V(x,\alpha)$, for any small positive constant $\varepsilon$, we choose a
particular $\varepsilon$-optimal strategy $\overline{m}(\cdot)$ for
(2.3)-(2.4) such that the approximation can be applied to the chain
$\\{(X^{h}_{n},\alpha^{h}_{n})\\}$ and the associated cost compared with
$V^{h}(x,\alpha)$. By the chattering lemma (see for instance (Kushner 1990,
Theorem 3.1)), for any given $\varepsilon>0$, there is a constant $\lambda>0$
and an ordinary control $\overline{U}^{\varepsilon}(\cdot)$ for (2.3)-(2.4)
with the following properties:

  1. (a)

$\overline{U}^{\varepsilon}(\cdot)$ takes only finitely many values (denoted
by $\mathcal{U}_{\varepsilon}$ the set of all such values);

  2. (b)

$\overline{U}^{\varepsilon}(\cdot)$ is constant on the intervals
$[k\lambda,k\lambda+\lambda)$ for $k\in\mathbb{Z}_{\geq 0};$

  3. (c)

with $\overline{m}^{\varepsilon}(\cdot)$ denoting the relaxed control
representation of $\overline{U}^{\varepsilon}(\cdot)$, we have that
$(\overline{X}^{\varepsilon}(\cdot),\overline{\alpha}^{\varepsilon}(\cdot),\overline{w}^{\varepsilon}(\cdot),\overline{m}^{\varepsilon}(\cdot))$
converges weakly to
$(\overline{X}(\cdot),\overline{\alpha}(\cdot),\overline{w}(\cdot),\overline{m}(\cdot))$
as $\varepsilon\to 0$;

  4. (d)

$J(x,\alpha,\overline{m}^{\varepsilon}(\cdot))\geq V(x,\alpha)-\varepsilon$.

For $\varepsilon>0$ and the corresponding $\lambda$ in the chattering lemma,
consider an optimal control problem for (2.3) subject to (2.4), but where the
controls are constants over the interval $[k\lambda,k\lambda+\lambda)$ for
$k\in\mathbb{Z}_{\geq 0}$ and take values in $\mathcal{U}_{\varepsilon}$ (the
set of control values of $\overline{U}^{\varepsilon}(\cdot)$). This
corresponds to controlling the discrete-time Markov process that is obtained
by sampling $X(\cdot)$ and $\alpha(\cdot)$ at times $k\lambda$ for
$k\in\mathbb{Z}_{\geq 0}$. Let $\widehat{U}^{\varepsilon}(\cdot)$ denote the
$\varepsilon$-optimal control, $\widehat{m}^{\varepsilon}(\cdot)$ denote the
relaxed control representation, and let $\widehat{X}^{\varepsilon}(\cdot)$
denote the associated state process. Since $\widehat{m}^{\varepsilon}(\cdot)$
is $\varepsilon$-optimal in the chosen class of controls, we have

$J(x,\alpha,\widehat{m}^{\varepsilon}(\cdot))\geq
J(x,\alpha,\overline{m}^{\varepsilon}(\cdot))-\varepsilon\geq
V(x,\alpha)-2\varepsilon.$

We next approximate $\widehat{U}^{\varepsilon}(\cdot)$ by a suitable function
of $w(\cdot)$ and $\alpha(\cdot)$. Using the same method as in Song et al.
(2006), we can approximate $\widehat{U}^{\varepsilon}(\cdot)$ by the ordinary
control $U^{\varepsilon,\theta}(\cdot)$ with the corresponding relaxed control
$m^{\varepsilon,\theta}(\cdot)$ and the state process
$X^{\varepsilon,\theta}(\cdot)$ such that

$m^{\varepsilon,\theta}(\cdot)\to\widehat{m}^{\varepsilon}(\cdot)$

as $\theta\to 0$ and

$J(x,\alpha,m^{\varepsilon,\theta}(\cdot))\geq
J(x,\alpha,\widehat{m}^{\varepsilon}(\cdot))-\varepsilon\geq
V(x,\alpha)-3\varepsilon.$

Then a sequence of ordinary controls $\\{\overline{U}^{h}_{n}\\}$ for the
chain $\\{(X^{h}(\cdot),\alpha^{h}(\cdot))\\}$ can be constructed with the
relaxed control representation $\\{\overline{m}^{h}_{n}\\}$ such that as $h\to
0$, the
$(X^{h}(\cdot),\alpha^{h}(\cdot),\overline{m}^{h}(\cdot),w^{h}(\cdot))$
converges weakly to
$(X^{\varepsilon,\theta}(\cdot),\alpha(\cdot),m^{\varepsilon,\theta}(\cdot),w(\cdot))$.
By the optimality of $V^{h}(x,\alpha)$ and the weak convergence above, we have
as $h\to 0$,

$V^{h}(x,\alpha)\geq J(x,\alpha,\overline{m}^{h}(\cdot))\to
J(x,\alpha,m^{\varepsilon,\theta}(\cdot)).$

It follows that $V^{h}(x,\alpha)\geq V(x,\alpha)-4\varepsilon$ for
sufficiently small $h$. Since any subsequence of $H^{h}(\cdot)$ has a
subsequence that converges weakly and $\varepsilon$ is arbitrary, we have
$\liminf\limits_{h}V^{h}(x,\alpha)\geq V(x,\alpha)$. The conclusion follows.
$\qquad\Box$

##  Appendix D Numerical Experiments

###  D.1. Varying the cost dependency

We want to see what effect different specifications of the cost function have
on the shape of the optimal harvesting rate, and in particular whether it is
bang-bang. We suspect that the convexity of the cost function leads to bang-
bang (all or nothing) optimal harvesting. In Figure 9, we have as an example a
cost functions of the form $C(u)=\sqrt{|u|})$. The rest of the parameters are
kept the same as in Subsection 4.1. This example has concave costs, but there
is a point of convexity at 0. Experiments with other partly concave cost
functions show a similar pattern. Piecewise linear costs like $C(u)=|u|$, or
$C(u)=\ln(1+|u|)$, lead to optimal controls that are step functions. However,
when we use a purely concave cost function, like $C(u)=\ln(1+u/3)$, seen in
Figure 10, we again obtain bang-bang optimal control. Further experiments
confirm the observation.

(a)  (b)

Figure 9. Value function (left) and optimal harvesting-stocking rate (right)
for a model with switching affecting $\mu(\alpha)=4-\alpha$, and a cost
function $C(u)=\sqrt{|u|}$. Other parameters described in Subsection 4.1.

(a)  (b)

Figure 10. Value function (left) and optimal harvesting-stocking rate (right)
for a model with switching affecting $\mu(\alpha)=4-\alpha$, and a cost
function $C(u)=\ln(1+u/3)$. Other parameters described in Subsection 4.1.

###  D.2. The Gompertz model of population growth

In this example, the dynamics of the population size without harvesting is
given by a Gompertz model (Winsor 1932, Zeide 1993) of the form

$dX(t)=\big{[}b(X(t),\alpha(t))-U(t)X(t)\big{]}dt+\sigma(X(t),\alpha(t))dw(t),$

where

$\displaystyle
b(x,\alpha)=(4-\alpha)x\ln\dfrac{2}{x},\quad\sigma(x,\alpha)=x,$
$\displaystyle\mathcal{U}=\\{u:u=k/500,k\in\mathbb{Z},-1000\leq k\leq
1500\\},\quad(x,\alpha)\in\mathbb{R}_{+}\times\\{1,2\\}.$

The generator $Q$ of the Markov chain $\alpha(\cdot)$ is given by

$q_{11}=-0.1,\quad q_{12}=0.1,\quad q_{21}=0.1,\quad q_{22}=-0.1.$

(a)  (b)

Figure 11. Value function (left) and optimal harvesting-stocking rate (right)
for a Gompertz model with absolute harvesting, with switching affectin
$b(x,\alpha)=(4-\alpha)x\ln\frac{2}{x}$, constant price $P(\cdot)=1$, and a
cost function $C(\cdot)=u^{2}/2$. Other parameters described in Subsection
4.1.

Figure 11 shows the value function and the optimal stocking-harvesting rate as
a function of population size $X(t)$ and the environmental state $\alpha$. In
the Gompertz model, the deterministic rate of growth near extinction goes to
$\infty$, unlike in the logistic model where it is linear. Comparing these
results with the ones in Figure 2, we can also see that low population values
in the Gompertz model are much less unfavorable, both in terms of future value
and in terms of the benefit of extraction.

###  D.3. The Nisbet-Gurney model of population growth

In this model, the evolution of the population size without harvesting is
given by a switched Nisbet-Gurney model; thus,

$dX(t)=\big{[}b(X(t),\alpha(t))-U(t)\big{]}dt+\sigma(X(t),\alpha(t))dw(t),$

where

$\displaystyle b(x,\alpha)=(4-\alpha)xe^{-x}-x,\quad\sigma(x,\alpha)=x,$
$\displaystyle\mathcal{U}=\\{u:u=k/500,k\in\mathbb{Z},-1000\leq k\leq
1500\\},\quad(x,\alpha)\in\mathbb{R}_{+}\times\\{1,2\\}.$

The generator $Q$ of the Markov chain $\alpha(\cdot)$ is given by

$q_{11}=-0.1,\quad q_{12}=0.1,\quad q_{21}=0.1,\quad q_{22}=-0.1.$

(a)  (b)

Figure 12. Value function (left) and optimal harvesting-stocking rate (right)
for a Nisbet-Gurney model with absolute harvesting, with switching affecting
the growth rate $b(x,\alpha)=(4-\alpha)xe^{-x}-x$, constant price
$P(\cdot)=1$, and a cost function $C(\cdot)=u^{2}/2$. Other parameters
described in Subsection 4.1.

Figure 12 shows a numerical estimation of this model. The value function has
the usual features, being increasing and concave. The harvesting rate is
monotonic, which is not a surprise considering the cost choice and our
discussion in Section 4.6. Again, the control in state $\alpha=1$ shows higher
harvesting and seeding, which is consistent with this state being more
favourable for growth.