\section{Introduction}
A fundamental problem in population biology has been to find conditions for when interacting species coexist or go extinct. Since the dynamics of interacting populations is invariably influenced by the random fluctuations of the environment, realistic mathematical models need to take into account the joint effects of biotic interactions and environmental stochasticity. A successful way of analysing the persistence and extinction of interacting species has been to look at Markov processes, in either discrete or continuous time, and describe their asymptotic properties. There has been a recent resurgence in stochastic population dynamics, and significant progress has been made for stochastic differential equations (\cite{SBA11, HN18}), piecewise deterministic Markov processes (\cite{BL16, HS19, HN19}), stochastic difference equations (\cite{BS19}), and general Markov processes (\cite{B18}). The first focus of this paper is to present new results for persistence and extinction in the setting of stochastic difference equations when there are one or two interacting species. These results significantly generalize the work by \cite{CE89, E89} which only treated competition models and had no extinction results as well as the more recent work by \cite{BS19} which only looks at compact state spaces. We are able to give explicit conditions for extinction and persistence in the setting of competitive or predator-prey Ricker equations with random coefficients, adding to the previously known results \cite{E89, VH97, FH02, SBA11}. Our results involve computing the invasion rates (\cite{T78, C82, E84, CE89}) of each species into the random equilibrium of the other species. We show that if both invasion rates are strictly positive then there is coexistence. If instead, one invasion rate is positive and one is negative, the species with the positive invasion rate persists while the one with the negative invasion rate goes extinct. If there is coexistence we prove that, under natural conditions, the populations converge to a unique invariant probability measure. If there is extinction we show that, with probability one, one or both species go extinct exponentially fast. The general theory for the setting with $n>2$ interacting species will appear in future work by the author and his collaborators (\cite{CHN19}).

Once criteria for persistence and extinction are established, our focus shifts towards a key problem from conservation ecology: what is the optimal strategy for harvesting species? This is a delicate issue as overharvesting can sometimes lead to extinction while underharvesting can mean the loss of precious economic resources. In continuous time models, recent studies have been able to find the optimal harvesting strategy, which maximizes either the discounted total yield or the asymptotic yield under very general assumptions if the ecosystem has only one species (\cite{Alvarez98, HNUW18, AH19}). For multiple species the theory is less developed. Nevertheless, partial results exist (\cite{LO97, ALO16, Ky15, Ky16, Ky18, HT19}).

Quite often harvesting models are intrinsically discrete in time. For example, if one looks at the management of fisheries, most models (\cite{GH89, HW92, C10, HL19}) assume that the population in a given year can be described by a single continuous variable and that without harvesting the population levels in successive years are related by
\[
x_{n+1}=F(x_n)
\]
where $F$ is the so-called recruitment function or the reproduction function. Most discrete time harvesting results ignore random environmental fluctuations and their effects on the availability of food, competition rates, growth and death rates, strength of predation and other key factors. Ignoring environmental stochasticity can create significant problems, in some cases making the models unrealistic (\cite{May78}) and hard to fit to data (\cite{L73}).
A series of key studies where environmental fluctuations are included was done by \cite{R78, R79}. Reed looked at the setting where there is one species whose dynamics in the absence of harvesting is given by
\[
X_{n+1}=Z_n F(X_n)
\]
where $(Z_n)_{n\in\mathbb{Z}_+}$ is a sequence of i.i.d. random variables. We extend Reed's analysis in two ways. First, we study the more general stochastic difference equation
\[
X_{n+1}=F(X_n, \xi_{n+1})
\]
where $(\xi_n)_{n\in\mathbb{Z}_+}$ is an i.i.d sequence. Second, we are able to analyze systems of two interacting species. To our knowledge these are the first results in discrete time that study the harvesting of multiple species.
\subsection*{Single species ecosystem.}
We are able to give exact conditions under which harvesting leads to persistence or to extinction. In particular, we show that if there is only one species present, then the criteria for persistence only involve the harvesting rate of the population at $0$. We are able to find the maximal harvesting rate which does not lead to extinction. If the species $Y_t$ undergoing harvesting persists we prove it converges in law to a random variable $Y_\infty$ and, if the fraction of the population that gets harvested is given by the strategy $h(y)$, we show that the long run average and the expected long term harvest both converge to $\mathbb{E} h(Y_\infty)$. In many applications one is interested in seeing how the environmental fluctuations change the long term yield. We show that in most settings the environmental fluctuations are detrimental and lower the harvesting yield. Only in special cases, can we have that the maximum deterministic sustainable yield is equal to the steady state harvest yield of the stochastic system.

An interesting corollary of our results is that threshold harvesting strategies (also called constant-escapement strategies), where one does not harvest anything below a threshold and harvests everything above that threshold, do not influence the persistence of species as long as the threshold at which one starts harvesting is strictly positive. We showcase two examples where environmental fluctuations are not detrimental for threshold harvesting: 1) The threshold $w$ at which we harvest is self-sustaining, i.e., if at the start of the year we are at level $w$ the fluctuations of the environment can not push the population's size under $w$. In this setting the expected value of the long term yield in the stochastic model equals the yield from the equivalent deterministic model. The downside is that the variance of the yield is higher due to the environmental fluctuations. 2) The threshold $w$ is not self sustaining and the maximum yield of the dynamics happens at a self-sustaining threshold $\overline x<w$. In this case the expected yield of the constant escapement strategy is strictly greater than the yield of the same strategy in the deterministic system.

When looking at constant effort harvesting strategies, where a constant proportion of the return is captured every year, we show that even though the deterministic model might say that we harvest at a sustainable rate, the environmental fluctuations might lead to extinction.

We are able to say more in the setting of the Ricker model. We give conditions under which we can get the same maximal yield in the deterministic and stochastic settings. This includes giving information about the threshold for which the yield is maximized. We also find the optimal harvesting strategy if we restrict ourselves to proportional harvest strategies.
\subsection*{Two interacting species.}
We analyze a system of two interacting species that can be exploited through harvesting. We show that threshold harvesting strategies do influence the persistence criteria, unlike in the single species setting. In order to be able to compute things explicitly we focus on Ricker, also called discrete time Lotka-Volterra, models and assume that the harvesting strategy is of proportional type, where we harvest a fraction $q\in [0,1]$ of the first species and a fraction $r\in [0,1]$ of the second species.

The first studied model is a predator-prey system where species $1$ is the prey and species $2$ the predator. We give analytical expressions for when one has the persistence of both species, the persistence of the prey and the extinction of the predator as well as the extinction of both species. These expressions tell us exactly for which rates $q,r$ we get one of the three scenarios above. Which strategy, among all proportional harvesting strategies, maximizes the expected long term harvesting yield? We find by using both analytical results and numerical simulations that it is never optimal to harvest both the predator and the prey. Either we drive the predator extinct and we harvest the prey or we do not harvest the prey at all and we harvest the predator.

The second model we look at consists of an ecosystem where the two species compete with each other for resources. We  show that, depending on the inter and intracompetition coefficients of the system, one can have two different regimes each having three regions which depend on the harvesting rates $q, r$:
\begin{itemize}
  \item [a)] I) Persistence of species $1$ and extinction of species $2$; II) Extinction of species $1$ and persistence of species $2$; III) Coexistence
  \item [b)] I) Persistence of species $1$ and extinction of species $2$; II) Extinction of species $1$ and persistence of species $2$; III) Bistability
\end{itemize}
We show that harvesting can facilitate coexistence in certain cases. When species $1$ is dominant and drives species $2$ extinct in the absence of harvesting, it is possible to harvest species $1$ and ensure the persistence of both species.

Finally, we look at the optimal harvesting strategies for the competitive system. Combining analytical proofs and numerical simulations we see that, in contrast to the predator-prey setting, it can be optimal, depending on the inter and intra competition rates, to harvest one or both of the species.

\section{Stochastic population dynamics}
We start by describing the stochastic population models we will be working with. To include the effects of random environmental fluctuations, ecologists and mathematicians often use stochastic difference equations of the form
\begin{equation}\label{e:SDE}
X^i_{t+1} = X_t^i f_i(\mathbf{X}_t, \xi_{t+1}).
\end{equation}
Here the vector $\mathbf{X}_t:=(X_t^1,\dots,X_t^n)\in \mathcal{S}\subset \mathbb{R}_+^n$ records the abundances of the $n$ populations at time $t\in\mathbb{Z}_+$ and $\xi_{t+1}$ is a random variable that describes the environmental conditions between time $t$ and $t+1$. The subset $\mathcal{S}$ will denote the state space of the dynamics. It will either be a compact subset of $\mathbb{R}_+^n$ or all of $\mathbb{R}_+^n$.
The coexistence set is the subset $\mathcal{S}_+=\{\mathbf{x}\in \mathcal{S}~|~x_i>0, i=1,\dots n\}$ of the state space where no species is extinct.
The real function $f_i(\mathbf{X}_t,\xi_{t+1})$ captures the fitness of the $i$-th population at time $t$ and depends both on the population sizes and the environmental state. Models of this type can capture complex short-term life histories and include predation, cannibalism, competition, and seasonal variations.

We have to differentiate between the setting where the dynamics is bounded, and the process enters and remains in a compact set, and the case when the dynamics is unbounded.  We will make the following assumptions throughout the paper:
\begin{itemize}
\item[\textbf{(A1)}] $\xi_1,\dots,\xi_n,\dots$ is a sequence of i.i.d. random variables taking values in a Polish space $E$.
\item[\textbf{(A2)}] For each $i$ the fitness function $f_i(\mathbf{x},\xi)$ is continuous in $\mathbf{x}$ on $\mathcal{S}$, measurable in $(\mathbf{x},\xi)$ and strictly positive.
\end{itemize}

Assumptions (A1) and (A2) ensure that the process $\mathbf{X}_t$ is a Feller process that lives on $\mathcal{S}_+$, i.e. $\mathbf{X}_t\in\mathcal{S}_+, t\in \mathbb{Z}_+$ whenever $\mathbf{X}_0\in\mathcal{S}_+$. One has to make extra assumptions (see (A3) or (A4) in Appendix \ref{s:a1}) in order to ensure the process does not blow up or fluctuate too abruptly between $0$ and $\infty$. We note that most ecological models will satisfy these assumptions. For more details see the work by \cite{BS19, CHN19}.
\begin{rmk}\label{r:1}
Suppose the dynamics is given by the more general model of the type
\begin{equation}\label{e:SDE2}
X^i_{t+1} = F_i(\mathbf{X}_t, \xi_{t+1}).
\end{equation}
Note that \eqref{e:SDE2} reduces to \eqref{e:SDE} if $F_i$ is $C^1$ and $F_i(\mathbf{x})=0$ whenever $x_i=0$.This means that $F_i$ is a nice, sufficiently smooth, vector field whick takes the value $0$ if species $i$ is extinct -- this is a natural assumption as there is no reason the population should be able to come back from extinction. Under these assumptions we can see that \eqref{e:SDE} is satisfied by setting
\begin{equation*}
   f_i(\mathbf{x},\xi) = \begin{cases}
              \frac{F_i(\mathbf{x},\xi)}{x_i} & \text{if } x_i>0,\\
               \frac{\partial F_i(\mathbf{x},\xi)}{\partial x_i} & \text{if } x_i = 0.
          \end{cases}
\end{equation*}
\end{rmk}

We will sometimes compare the stochastic model \eqref{e:SDE} with its averaged deterministic counterpart
\begin{equation}\label{e:det}
\mathbf{x}_{t+1}^i=x_t^i \overline f_i(\mathbf{x}_t)
\end{equation}
where $\overline f_i(\mathbf{x}) := \mathbb{E} f_i(\mathbf{x},\xi_1)$. Note that
\[
\mathbb{E}[X^i_{t+1}~|~\mathbf{X}_t=\mathbf{x}] = x^i\mathbb{E} f_i(\mathbf{x},\xi_1) = x^i  \overline f_i(\mathbf{x}),
\]
so that \eqref{e:det} is the average of \eqref{e:SDE} in this sense.

For example, if $f(x,\xi)=\xi u(x)$ and $\xi_1$ is a random variable with expectation $\mathbb{E} \xi_1=1$ then
\[
\mathbb{E}[X^i_{t+1}~|~\mathbf{X}_t=\mathbf{x}] = x^i u(\mathbf{x}).
\]
This is the setting used by \cite{R78}.
\subsection{Stochastic persistence.}
We define the extinction set, where at least one species is extinct, by
\[
\mathcal{S}_0:=\mathcal{S}\setminus \mathcal{S}_+=\{\mathbf{x}\in\mathcal{S}~:~\min_i x_i=0\}.
\]
The transition operator $P: \mathcal{B}\to \mathcal{B}$ of the process $\mathbf{X}$ is an operator which acts on Borel functions $\mathcal{B}:=\{h:\mathcal{S}\to\mathbb{R}~|~h~\text{Borel}\}$ as
\[
Ph(\mathbf{x})=\mathbb{E}_\mathbf{x}[h(\mathbf{X}(1))]:=\mathbb{E}[h(\mathbf{X}(1))~|~\mathbf{X}(0)=\mathbf{x}],  ~\mathbf{x}\in\mathcal{S}.
\]
The operator $P$ acts by duality on Borel probability measures $\mu$ by $\mu\to \mu P$ where $\mu P$ is the probability measure given by
\[
\int_{\mathcal{S}}h(\mathbf{x}) (\mu P)(d\mathbf{x}):=\int_{\mathcal{S}}P h(\mathbf{x}) \mu(d\mathbf{x})
\]
for all $h\in C(\mathcal{S})$.
A Borel probability measure $\mu$ on $\mathcal{S}$ is called an \textit{invariant probability measure} if
\[
\mu P = \mu
\]
where $P$ is the transition operator of the Markov process $\mathbf{X}_t$. An invariant probability measure or stationary distribution is a way of describing a `random equilibrium'. If the process starts with $\mathbf{X}_0$ having an initial distribution given by the invariant probability measure $\mu$ then the distribution of $\mathbf{X}_t$ is $\mu$ for all $t\in\mathbb{Z}_+$. In a sense this is the random analogue of a fixed point of a deterministic dynamical system. It turns out that a key concept is the \textit{realized per-capita growth rate} of species $i$ when introduced in the community described by an invariant probability measure $\mu$
\begin{equation}\label{e:r}
r_i(\mu) =\int_{\mathbb{R}_+^n}\mathbb{E}[\ln f_i(\mathbf{x},\xi_1)]\,\mu(d\mathbf{x}) = \int r_i(\mathbf{x})\mu(d\mathbf{x})
\end{equation}
where
\[
r_i(\mathbf{x}) = \mathbb{E} [\ln f_i(\mathbf{x},\xi_1)]
\]
is the mean per-capita growth rate of species $i$ at population state $\mathbf{x}$. This quantity tells us whether species $i$ tends to increase or decrease when introduced at an infinitesimally small density into the the subcommunity described by $\mu$. If the $i$th species is among the ones supported by the subcommunity given by $\mu$, i.e., $i$ lies in the support of $\mu$, then this species is in a sense `at equilibrium' and one can prove that
\begin{equation}\label{e:equ}
r_i(\mu)=0.
\end{equation}
The only directions $i$ in which $r_i(\mu)$ can be non-zero are those which are not supported by $\mu$.

One can show that the invariant probability measures living on the extinction set $\mathcal{S}_0$ fully describe the long term behavior of the system. In a sense, if any such invariant probability measure is a \textit{repeller} which pushes the process away from the boundary in at least one direction, then the system persists. Let $\Conv(\mathcal{M})$ denote the set of all invariant probability measures supported on $\mathcal{S}_0$. In order to have the convergence of the process to a unique stationary distribution one needs some irreducibility conditions which keep the process from being too degenerate \citep{CHN19, MT, B18}. The following theorem characterizes the coexistence of the ecosystem.
\begin{thm}\label{t:pers1_disc}
Suppose that for all $\mu\in\Conv(\mathcal{M})$ we have
\begin{equation}\label{e:pers}
\max_{i} r_i(\mu)>0.
\end{equation}
Then the system is almost surely stochastically persistent and stochastically persistent in probability. Under additional irreducibility conditions, there exists a unique invariant probability measure $\pi$ on $\mathcal{S}_+$ and as $t\to\infty$ the distribution of $\mathbf{X}_t$ converges in total variation to $\pi$ whenever $\mathbf{X}(0)=\mathbf{x}\in\mathcal{S}_+$. Furthermore, if $w:\mathcal{S}_+\to\mathbb{R}$ is bounded
then $$\lim_{t\to\infty}\mathbb{E} w(\mathbf{X}_t) = \int_{\mathcal{S}_+} w(\mathbf{x})\,\pi(d\mathbf{x}).$$
\end{thm}
A sketch of the proof of this result appears in Appendix \ref{s:a1}.
\subsection{Classification of two species dynamics}\label{s:2d}
Sometimes one is not only interested in persistence and coexistence, but also in conditions which lead to extinction. Extinction results are more delicate and require a technical analysis. Some extinction results appeared in work by \cite{HN18, BS19}. The sharpest and most complete results for stochastic difference equations and species feedbacks will appear in \cite{CHN19}. We restrict our discussion to a system with two species.
In this setting \eqref{e:SDE} becomes
\begin{equation}\label{e:2d_sys}
\begin{aligned}
X^1_{t+1}&=X^1_t f_1(X^1_t, X^2_t,\xi_{t+1}),\\
X^2_{t+1}&=X^2_t f_2(X^1_t, X^2_t,\xi_{t+1}).
\end{aligned}
\end{equation}
The exact assumptions and technical results are found in Appendix \ref{s:a2d}. We can classify the dynamics as follows. We first look at the Dirac delta measure $\delta_0$ at the origin $(0,0)$
\[
r_i(\delta_0) = \mathbb{E}[\ln f_i(0,\xi_1)], i=1,2.
\]
If $r_i(\delta_0)>0$ then species $i$ survives on its own and converges to a unique invariant probability measure $\mu_i$ supported on $\mathcal{S}_+^i := \{\mathbf{x}\in\mathcal{S}~|~x_i\neq 0, x_j=0, i\neq j\}$. The realized per-capita growth rates can be computed as
  \[
  r_i(\mu_j)=\int_{(0,\infty)}\mathbb{E}[\ln f_i(x,\xi_1)]\mu_j(dx).
  \]
\begin{enumerate}[label=(\roman*)]
  \item Suppose $r_1(\delta_0)>0, r_2(\delta_0)>0$.
  \begin{itemize}
    \item If $r_1(\mu_2)>0$ and $r_2(\mu_1)>0$ we have coexistence and convergence of the distribution of $\mathbf{X}_t$ to the unique invariant probability measure $\pi$ on $\mathcal{S}_+$.
    \item If $r_1(\mu_2)>0$ and $r_2(\mu_1)<0$ we have the persistence of $X^1$ and extinction of $X^2$.
    \item If $r_1(\mu_2)<0$ and $r_2(\mu_1)>0$ we have the persistence of $X^2$ and extinction of $X^1$.
    \item If $r_1(\mu_2)<0$ and $r_2(\mu_1)<0$ we have that for any $\mathbf{X}_0=\mathbf{x}\in\mathcal{S}_+$
    \[
    p_{\mathbf{x},1}+ p_{\mathbf{x},2}=1,
    \]
    where $p_{\mathbf{x},j}>0$ is the probability that species $j$ persists and species $i\neq j$ goes extinct.
  \end{itemize}
  \item Suppose $r_1(\delta_0)>0, r_2(\delta_0)<0$. Then species $1$ survives on its own and converges to its unique invariant probability measure $\mu_1$ on $\mathcal{S}^1_+$.
  \begin{itemize}
    \item If $r_2(\mu_1)>0$ we have the persistence of both species and convergence of the distribution of $\mathbf{X}_t$ to the unique invariant probability measure $\pi$ on $\mathcal{S}_+$.
    \item If $r_2(\mu_1)<0$ we have the persistence of $X^1$ and the extinction of $X^2$.
  \end{itemize}
\item Suppose $r_1(\delta_0)<0, r_2(\delta_0)<0$. Then both species go extinct with probability one.
\end{enumerate}
We note that our results are significantly more general than those from \cite{E89}. In \cite{E89} the author only gives conditions for coexistence, and does not treat the possibility of the extinction of one or both species.

\begin{example}\label{ex:1}
The simplest case is when the noise is multiplicative, that is
\begin{equation}\label{e:2d_syss}
\begin{aligned}
X^1_{t+1}&=X^1_t Z_{t+1}^1f_1(X^1_t, X^2_t)\\
X^2_{t+1}&=X^2_t Z_{t+1}^2 f_2(X^1_t, X^2_t),
\end{aligned}
\end{equation}
where $Z^1_{1}, Z^1_{2},\dots$ is an i.i.d sequence of random variables and $Z^2_{1}, Z^2_{2}, \dots$ is an independent sequence of i.i.d random variables.
In this case for $i=1,2$ we have
\begin{equation}\label{e:2d_pers}
\begin{aligned}
r_i(\delta_0) &= \mathbb{E}[\ln (Z_{t+1}^i f_i(0))]\\
&= \mathbb{E} \ln Z_1 + \ln f_i(0).
\end{aligned}
\end{equation}
The growth rates at $0$ in the stochastic model differ from the growth rates at $0$ of the deterministic model only by the term $\mathbb{E} \ln Z_1$.
\end{example}

\subsection{Harvesting}\label{s:harvesting}
We next describe how the harvesting effects are taken into account. We assume that the harvesting takes place during a short harvest season. The size of the population at the beginning of the harvest season in year $t$ will be denoted by $\mathbf{Y}_t$ and will be called \textit{return in year $t$}. If we assume the harvest season is short so that growth and natural mortality can be neglected during the harvesting and that the harvesting strategy is \textit{stationary}, i.e., the size of the harvest in any year depends only on the size of the population return $\mathbf{Y}$ in that year we can write
\begin{equation}\label{e:harv1}
X^i_t = Y_t^i - h_i(\mathbf{Y}_t) = u_i(\mathbf{Y}_t)
\end{equation}
where $X^i_t$ is the escapement of the $i$th population from the harvest and $h_i(\mathbf{Y}_t)$ is the amount of species $i$ that is harvested at time $t$. The function $u_i$ is called the \textit{escapement function} and measures how much is left after harvesting. Note that, since we cannot harvest a negative amount or more than the total population size we will always have
\[
0\leq h_i(\mathbf{y})\leq y_i.
\]
 Set $\mathbf{u}(\mathbf{y}):=(u_1(\mathbf{y}),\dots,u_n(\mathbf{y}))$. Once the harvesting is done, the population evolves according to \eqref{e:SDE} so that the size of the return in year $t+1$ is related to the escapement in year $t$ via
\begin{equation}\label{e:harv2}
Y^i_{t+1} = X_t^i f_i(\mathbf{X}_t, \xi_{t+1}).
\end{equation}
Combining \eqref{e:harv1} and \eqref{e:harv2} we get
\begin{equation}\label{e:harv}
Y^i_{t+1} = u_i(\mathbf{Y}_t)  f_i(\mathbf{u}(\mathbf{Y}_t), \xi_{t+1}).
\end{equation}
In order to be able to analyze the process $\mathbf{Y}_t$ we have to make sure that it can be written in the \textit{Kolmogorov} form \eqref{e:SDE}. In order to get this we assume that
\begin{asp}\label{a:1}
For all $i=1,\dots,n$ the following properties hold
\begin{itemize}
\item[(a)] The function $u_i$ is strictly positive on $\mathcal{S}_+$, with
\[
u_i(\mathbf{y}) \leq y_i.
\]
\item [(b)] The function $u_i$ is continuous on $\mathcal{S}$ and continuously differentiable at $y_i=0$.
\end{itemize}
\end{asp}
\begin{rem}\label{r:2}
Note that Assumption \ref{a:1} implies that $u_i(\mathbf{y})=0$ if $y_i=0$ and
\[
\frac{\partial u_i}{\partial y_i}(\mathbf{y})=\lim_{y_i\to 0} \frac{u_i(\mathbf{y})}{y_i}\leq 1
\]
if $\mathbf{y}\in\mathcal{S}$ with $y_i=0$.
\end{rem}

\subsection{Persistence with harvesting.} Since overharvesting can lead to extinction, we want to find sufficient conditions which ensure the process $\mathbf{Y}_t$ converges to a unique invariant probability measure on $\mathcal{S}$.  Note that we need to put \eqref{e:harv} into the form \eqref{e:SDE}. For this, using Remark \ref{r:1}, let
\begin{equation}\label{e:g}
   g_i(\mathbf{x},\xi) := \begin{cases}
              \frac{u_i(\mathbf{x})}{x_i} f_i(\mathbf{u}(\mathbf{x}), \xi)& \text{if } x_i>0,\\
               \left(\frac{\partial u_i}{\partial x_i}(\mathbf{x})\right)f_i(\mathbf{u}(\mathbf{x}), \xi) & \text{if } x_i = 0.
          \end{cases}
\end{equation}
We can write \eqref{e:harv} as
\begin{equation}\label{e:harv3}
Y^i_{t+1} = Y^i_t  g_i(\mathbf{Y}_t, \xi_{t+1}).
\end{equation}
In order to use Theorem \ref{t:pers1_disc} we have to make sure that conditions A1)-A4) are satisfied, that the process $\mathbf{Y}_t$ is $\phi$-irreducible and that \eqref{e:pers} holds.
If $\mu$ is an invariant probability measure of $\mathbf{Y}_t$ living on the extinction set $\mathcal{S}_+$ the realized per capita growth rates will be given by
\begin{equation}\label{e:r2}
r_i(\mu) =\int_{ \mathcal{S}}\mathbb{E}[\ln g_i(\mathbf{x},\xi_1)]\,\mu(d\mathbf{x}).
\end{equation}
Specifically, if we look at the Dirac mass at $0$ we get
\begin{equation}\label{e:nd_pers3}
\begin{aligned}
r_i(\delta_0)&=\int_{\mathcal{S}_0}\mathbb{E}[\ln g_i(\mathbf{x},\xi_1)]\,\delta_0(d\mathbf{x})\\
&= \mathbb{E}[\ln g_i(0,\xi_1)] \\
&= \mathbb{E}\left[\ln \left(\frac{\partial u_i}{\partial x_i}(0)f_i(0, \xi)\right)\right].
\end{aligned}
\end{equation}
\begin{example}
If the noise is multiplicative and we are in the setting of Example \ref{ex:1}, i.e.,
\begin{equation}\label{e:nd_sys}
\begin{aligned}
Y^i_{t+1}&=Y^i_t Z_{t+1}^ig_i(\mathbf{Y}_t)
\end{aligned}
\end{equation}
then
\begin{equation}\label{e:nd_pers2}
\begin{aligned}
r_i(\delta_0) &= \mathbb{E}[\ln (Z_{t+1}^i g_i(0))]\\
&= \mathbb{E} \ln Z_1 +\ln \left(\frac{\partial u_i}{\partial x_i}(0)f_i(0)\right)\\
&= \mathbb{E} \ln Z_1 +\ln \left(\frac{\partial u_i}{\partial x_i}(0)\right) + \ln f_i(0).
\end{aligned}
\end{equation}
\textit{\textbf{Biological interpretation:} The noise can either help or inhibit the survival of the $i$th species, according to whether $\mathbb{E} \ln Z_1<0$ or $\mathbb{E} \ln Z_1>0$. Since $\frac{\partial u_i}{\partial x_i}(0)\leq 1$ we always have $\ln \left(\frac{\partial u_i}{\partial x_i}(0)\right)\leq 0$ so that, as expected, harvesting is always detrimental to the survival of each individual species. The minimal escapement rate under which one can still have persistence is such that
\[
\frac{\partial u^{min}_i}{\partial x_i}(0) > e^{- (\mathbb{E} \ln Z_1 + \ln f_i(0))} = f_i(0)  e^{- \mathbb{E} \ln Z_1}.
\]
Since $h_i(\mathbf{y})=1-u_i(\mathbf{y})$ this implies that the maximal harvesting rate satisfies
\[
\frac{\partial h^{max}_i}{\partial x_i}(0) < 1- e^{- (\mathbb{E} \ln Z_1 + \ln f_i(0))}=1- f_i(0)  e^{- \mathbb{E} \ln Z_1}.
\]
}
\end{example}

\section{Single species harvesting}\label{s:1s}
This section explores the setting when there is only one species in the ecosystem. The results can be seen as an extension of the results from \cite{R78}. Our results show that environmental fluctuations are usually detrimental to the optimal harvesting yield. Actually, only under very special conditions, is it possible for the stochastic dynamics to have the same maximal expected long term yield as the related deterministic dynamics. Even in that case, the nonzero variance of the stochastic long term yield, makes it more risky than its deterministic analogue.

We can show that in some special cases of constant-escapement strategies, it is possible for the stochastic expected long term yield to be higher than the deterministic yield.

One can see from \eqref{e:harv3} that the dynamics of the return will be given by
\begin{equation}\label{e:harv4}
Y_{t+1} = Y_t  g(Y_t, \xi_{t+1})
\end{equation}
for
\begin{equation*}
   g(x) := \begin{cases}
              \frac{u(x)}{x} f(u(x), \xi)& \text{if } x>0,\\
               \left(\frac{\partial u}{\partial x}(0)\right)f(u(0), \xi) & \text{if } x = 0.
          \end{cases}
\end{equation*}
We will work under the assumption that without harvesting we have
\begin{equation}\label{e:r_single_pers}
\mathbb{E}\left[\ln \left(f(0, \xi_1)\right)\right]>0
\end{equation}
so that the species persists.
Suppose the assumptions of one of the Theorems \ref{t:pers_single1}, \ref{t:pers_single2}, \ref{t:pers_single3} or \ref{t:pers1_disc} hold. Then, in order to have persistence we need
\begin{equation}\label{e:r_single2}
r_1(\delta_0) =\int_{\partial \mathbb{R}_+}\mathbb{E}[\ln g(x,\xi_1)]\,\delta_0(dx) = \mathbb{E}\left[\ln \left(\left(\frac{\partial u}{\partial x}(0)\right)f(0, \xi_1)\right)\right]>0
\end{equation}
where $\delta_0$ is the point mass at $0$ and we made use of \eqref{e:r2} and Assumption \ref{a:1}.
We can express this result as
\begin{equation}\label{e:r_single}
\frac{\partial u}{\partial x}(0)> e^{-  \mathbb{E}\ln f(0, \xi_1)}.
\end{equation}
Let us next compute the expected long term harvest yield. If the assumptions of Theorem \ref{t:pers1_disc} are satisfied we will have
\[
Y_{t}\to Y_\infty
\]
in distribution as $t\to\infty$. Here $Y_\infty$ is a random variable whose distribution equal to the the invariant probability measure $\pi_u$. In many models, and for well behaved functions $h$ one can show by Theorem \ref{t:pers1_disc} that $\mathbb{E} h(Y_\infty)$ exists and is finite. As a result, we have that with probability one
\[
\lim_{T\to \infty} \frac{\sum_{t=0}^Th(Y_t)}{T} =  \mathbb{E} h(Y_\infty).
\]
This tells us that the \textit{long-run average harvest yield} converges to a steady yield $\mathbb{E} h(Y_\infty)$. Furthermore, we can also see that \textit{expected yield} also converges to the same quantity
\[
\lim_{T\to \infty} \mathbb{E} h(Y_T) = \mathbb{E} h(Y_\infty) = \int_{(0,\infty)}h(x)\,\pi_u(dx).
\]
From now on we will call $\mathbb{E} h(Y_\infty)$ the \textit{expected steady-state yield}.
In general it is not possible to find $\mathbb{E} h(Y_\infty)$. However, in certain instances we can exploit the fact that, at stationarity, the realized per-capita growth rates in the directions supported by the measure $\pi_u$ are all zero \citep{CHN19}. In other words
\begin{equation}\label{e:r0}
0=r_1(\pi) = \int_{(0,\infty)}\mathbb{E}\left[\ln\left( \frac{u(x)}{x} f(u(x), \xi_1)\right)\right]\,\pi_u(dx).
\end{equation}
\subsection{Stochastic versus deterministic harvesting.}\label{s:stocdet}
Let us compare the stochastic dynamics \eqref{e:harv4} with its deterministic average
\begin{equation}\label{e:det1}
x_{t+1}=x_t \overline f(x_t)
\end{equation}
where $\overline f(x) := \mathbb{E} f(x,\xi_1)$ and $\overline F(x) := x \overline f(x)$. If $h$ is any stationary harvesting strategy the deterministic equilibrium return $y$ satisfies
\[
y=\overline F(u(y))
\]
and the equilibrium yield is
\[
h(y)=y-u(y) = \overline F(u(y))-u(y)=\overline G (u(y))
\]
where $\overline G(x):=\overline F(x)-x$. The \textit{deterministic maximum sustainable yield} (DMSY) is obtained by keeping the escapement $u(y)$ at the level $x_1$ at which $\overline G$ attains its maximum, i.e., at the point $x_1$ where
\[
0=\overline G'(x_1)=\overline F'(x_1)-1=\overline f(x_1) + x_1 \overline f'(x_1)-1.
\]
The DMSY $M_{\det}$ will be
\[
M_{\det} = \overline G(x_1)=x_1\overline f(x_1)-x_1.
\]

\begin{thm}\label{t:maximal_yield}
The expected value of the steady state harvest yield $\mathbb{E} h(Y_\infty)$ of any stationary harvesting policy $h$ of the model \eqref{e:harv4} is always dominated by the maximum deterministic sustainable yield of the equivalent deterministic model \eqref{e:det1},
\[
\mathbb{E} h(Y_\infty) \leq M_{\det}.
\]
 The only way to achieve an equality in the above is when the following conditions are satisfied:
\begin{itemize}
  \item [1)] The unharvested dynamics $X_{t+1}=X_tf(X_t,\xi_{t+1})$ is able to go to a level greater or equal to $x_1$.
  \item [2)] The harvesting policy is bang-bang with threshold $x_1$, that is
  \begin{equation*}
   h^*(y) := \begin{cases}
             y-x_1 & \text{if } y>x_1,\\
               0 & \text{if } y\leq x_1.
          \end{cases}
\end{equation*}
  \item [3)] The level $x_1$ is self-sustaining i.e. the stochastic effects do not make the population ever go below $x_1$ once it reaches this level.
\end{itemize}
\end{thm}
\begin{proof}
Let $G(x,\xi):= x f(x,\xi)-x$ and note that $\mathbb{E} G(x,\xi_1)= \overline G(x)$. For the stochastic model if we use the harvesting policy $h$, the long run average yield is
\begin{equation}\label{e:hY}
\begin{split}
  \mathbb{E} h(Y_\infty) &=  \mathbb{E} Y_\infty - \mathbb{E}[u(Y_\infty)] \\
   &= \mathbb{E} u(Y_\infty) f(u(Y_\infty),\xi_1) - \mathbb{E}[u(Y_\infty)] \\
   &= \mathbb{E} G(u(Y_\infty)),\xi_1)\\
   &= \int \int G(u(y), \xi) \mathbb{P}(Y_\infty\in dy, \xi_1\in d\xi)\\
   &= \int \int G(u(y), \xi) \mathbb{P}(Y_\infty\in dy)\mathbb{P}(\xi_1\in d\xi)\\
   &= \int \left(\int G(u(y), \xi)\mathbb{P}(\xi_1\in d\xi)\right) \mathbb{P}(Y_\infty\in dy)\\
   &= \mathbb{E} \overline G (u(Y_\infty))
\end{split}
\end{equation}
where we used the fact that $u(Y_\infty)$ and $u(Y_\infty) f(u(Y_\infty),\xi_1)$ have the same distribution and $Y_\infty$ is independent of $\xi_1$. Since $\overline G$ attains its maximum at $x_1$ we have
\begin{equation}\label{e:hYin}
 \mathbb{E} h(Y_\infty) = \mathbb{E} \overline G (u(Y_\infty))\leq \overline G(x_1)= M_{\det}.
\end{equation}
In order to have equality in \eqref{e:hYin} we need the law of $u(Y_\infty)$ to be the point mass $\delta_{x_1}$ at $x_1$. This means that with probability 1
\[
Y_\infty - h(Y_\infty) = x_1.
\]
One can achieve this if:
\begin{itemize}
\item [1)] the population can get to a level that is equal or greater to $x_1$,
\item [2)] one uses the \textit{bang-bang}, also called \textit{constant escapement} or \textit{threshold}, harvest policy at the level $x_1$
\begin{equation*}
   h^*(y) := \begin{cases}
             y-x_1 & \text{if } y>x_1,\\
               0 & \text{if } y\leq x_1,
          \end{cases}
\end{equation*}
and
\item [3)] once the population reaches the level $x_1$, it never decreases to a lower level, that is, if $X_t=x_1$ then
\[
Y_{t+1} = X_1f(X_1,\xi_1)=x_1f(x_1,\xi_1) \geq x_1
\]
with probability $1$.
\end{itemize}
 The last property is equivalent to having $\mathbb{P}(f(x_1,\xi_1)\geq 1)=1$ -- if this is true, we say that the level $x_1$ is \textit{self-sustaining}.
\end{proof}
\textit{\textbf{Biological interpretation:} In general, it is not possible to have the same maximal yield in the stochastic setting as in the deterministic setting. Due to environmental fluctuations the expected long term yield of any harvesting strategy $h$ will be dominated by the deterministic maximum sustainable yield. The only case when the maximal yields in the stochastic and deterministic setting are equal, is when one uses a constant escapement strategy with threshold $x_1$ (which maximizes the deterministic MSY), the stochastic dynamics can reach levels greater or equal to $x_1$ and $x_1$ and then never goes below $x_1$ due to environmental fluctuations. These very specific conditions will not usually hold. As such, for most situations we cannot expect to get the same optimal harvesting yields in the deterministic and stochastic settings.
}

We note that threshold (or bang-bang) harvesting strategies do not influence the persistence criterion in the one dimensional case. The unharvested system
\begin{equation}\label{e:SDE_1d}
X_{t+1} = X_t f(X_t, \xi_{t+1})
\end{equation}
has
\[
r^X_1(\delta_0) = \mathbb{E} \ln f(0,\xi_1) >0.
\]
If one adds harvesting, then
\[
Y_{t+1} = Y_t  g(Y_t, \xi_{t+1})
\]
where
\begin{equation*}
   g(x) := \begin{cases}
              \frac{u_w(x)}{x} f(u_w(x), \xi)& \text{if } x>0,\\
               \left(\frac{\partial u_w}{\partial x}(0)\right)f(u_w(0), \xi) & \text{if } x = 0.
          \end{cases}
\end{equation*}
The bang-bang strategy
\begin{equation*}
   h_w(y) := \begin{cases}
             y-w & \text{if } y>w,\\
               0 & \text{if } y\leq w,
          \end{cases}
\end{equation*}
with $w>0$ also has
\[
r^Y_1(\delta_0) = \mathbb{E} \ln g(0,\xi_1)= \mathbb{E} \ln \left(\frac{\partial u_w}{\partial x}(0)\right)f(u(0), \xi) = \mathbb{E} \ln f(0,\xi_1)=r^X_1(\delta_0).
\]
\textit{\textbf{Biological interpretation:} This implies that a constant escapement strategy with a threshold $w>0$ does not change the per-capita growth rate and thus does not interfere with persistence. This is one reason why bang-bang harvesting strategies are robust and make sense when there is only one species present. This is not the case anymore when there are multiple species present.}

\subsection{Constant effort harvesting}
Quite often in fisheries a constant effort harvesting method is used. These strategies are such that the same fixed proportion of the return is captured every year. In other words, for some fixed $\theta \in (0,1)$ we have
\[
h_\theta(x)= \theta x
\]
and
\[
u_\theta(x)=(1-\theta)x.
\]
The persistence criteria \eqref{e:r_single} becomes
\[
\theta < \theta_{\max} := 1 - e^{-\mathbb{E}\ln f(0,\xi_1)}
\]
where $\theta_{\max}$ is the \textit{maximum sustainable rate of exploitation}.
Let us compare this with the deterministic system
\[
x_{t+1}=x_t \overline f(x_t) =\overline F(x_t)
\]
where $\overline f(x) = \mathbb{E} f(x,\xi_1)$. In this setting the maximum sustainable rate of exploitation $\theta_{\det}$ is given by
\[
\theta_{\det} = 1 - \frac{1}{\overline f(0)}=1- e^{-\ln \mathbb{E} f(0,\xi_1) }.
\]
Since the logarithm is a concave function, Jensen's inequality implies that
\[
 \mathbb{E} \ln f(0,\xi_1) \leq \ln \mathbb{E} f(0,\xi_1),
\]
with equality if and only if $f(0,\xi_1)$ is constant with probability one. As a result
\[
\theta_{\max}<\theta_{\det},
\]
which was shown by \cite{R78} in a simpler model.

\textit{\textbf{Biological interpretation:} This inequality is important because it says that if one uses a deterministic model, then the rate of exploitation might seem to be sustainable. Nevertheless, using this strategy will lead to extinction due to the presence of random environmental fluctuations.}

It is well known that in the setting of \eqref{e:det1} the deterministic maximum sustainable yield (DMSY) is achieved when the rate of exploitation is
\[
\theta_{DMSY} = 1 - \frac{1}{\overline f(x_1)}
\]
for $x_1$ satisfying $$\overline f'(x_1)=1.$$ Under environmental conditions which are large enough we can have
\[
e^{-\ln \mathbb{E} f(0,\xi_1) } > \frac{1}{\overline f(x_1)}.
\]
\textit{\textbf{Biological interpretation:} If the environmental fluctuations are significant, one has $\theta_{DMSY}>\theta_{\max}$. This shows that if large environmental fluctuations are possible and we harvest the population according to the deterministic MSY rate of exploitation we will drive it to extinction.}
\begin{thm}\label{t:CE}
If the deterministic averaged system \eqref{e:det1} has a strictly concave $\overline F$, and the dynamics \eqref{e:SDE_1d} is not purely deterministic, then the asymptotic expected yield of any constant effort harvest strategy is strictly lower than the deterministic yield of that harvesting policy.
\end{thm}
\begin{proof}
For a constant-effort policy $h_\theta(y)=\theta y$ the asymptotic expected yield is given by
\[
\mathbb{E} h_\theta(Y_\infty)=\theta \mathbb{E} Y_\infty.
\]
Set $u_\theta(y)=(1-\theta)y$. Using that $ u_\theta(Y_\infty) f(u_\theta(Y_\infty),\xi_1)$ and $u_\theta(Y_\infty)$ have the same distribution an argument similar to the one from \eqref{e:hY} shows that
\[
\mathbb{E} Y_\infty = \mathbb{E} Y_\infty f(u_\theta(Y_\infty),\xi_1) = \frac{1}{1-\theta}\mathbb{E} \overline F(u_\theta(Y_\infty)).
\]
If the function
\[
\overline F (y) = y \mathbb{E} f(y,\xi_1)
\]
is strictly concave then by Jensen's inequality
\begin{equation}\label{e:ineq}
 \frac{1}{1-\theta}\mathbb{E} \overline F(u_\theta(Y_\infty)) \leq  \frac{1}{1-\theta} \overline F(\mathbb{E} u_\theta(Y_\infty)).
\end{equation}
One can have equality in \eqref{e:ineq} only if $Y_\infty$ is with probability one a constant random variable. This implies that
\[
\mathbb{E} Y_\infty \leq \overline F ((1-\theta) \mathbb{E} Y_\infty).
\]
As we know, in the deterministic model, using the same policy with harvest rate $\theta$, the equilibrium return $\widehat y_\theta$ satisfies
\[
\widehat y_\theta = \overline F((1-\theta)\widehat y_\theta).
\]
This together with the strict concavity of $\overline F$ implies that
\[
\mathbb{E} Y_\infty \leq \widehat y_\theta.
\]
Equality can only hold if $Y_\infty$ is with probability one a constant, which means the dynamics is deterministic.
\end{proof}

\subsection{Bang-bang threshold harvesting.} Bang-bang or constant-escapement harvesting strategies are important and are used in many theoretical models as well as in actual harvesting situations, like fisheries. These policies have been shown to be optimal in many instances both for the continuous (\cite{LO97,Alvarez98, HNUW18, AH19}) and discrete time (\cite{R78, R79}) settings. Constant escapement strategies turned out to be optimal for maximizing discounted yield, asymptotic yield, as well as discounted economic revenue under many different conditions. In discrete time the work by \cite{R79} implies that a bang-bang policy maximizes the expected discounted net revenue in a discrete time stochastic model. It has not been shown in discrete time, to our knowledge, in generality that the expected steady state yield is maximized under a bang-bang strategy. However, both heuristic arguments and analytical results in specific cases hint that these strategies are probably the ones that will in general be optimal. In addition, these are the strategies that are most widely used in fisheries, where the escapement is controlled. We will explore how well these bang-bang strategies do in the stochastic harvesting setting \eqref{e:harv4} in comparison to the deterministic setting \eqref{e:det1}.

Suppose we harvest according to the bang-bang strategy
\begin{equation}\label{e:bang}
   h_w(y) := \begin{cases}
             y-w & \text{if } y>w,\\
               0 & \text{if } y\leq w,
          \end{cases}
\end{equation}
with $w>0$.
Let $r_s$ be the maximum self-sustaining level
\[
r_s = \max \{x~|~\mathbb{P}(f(x,\xi_1)\geq 1)=1\}.
\]
We have to differentiate between two cases:
\begin{itemize}
  \item[1)] The level $w$ is self-sustaining, i.e.
  \[
  \mathbb{P}(f(w,\xi_1)\geq 1)=1.
  \]
  \item[2)] The level $w$ is not self-sustaining.
\end{itemize}
\begin{prop}
If the threshold level $w$ is self-sustaining then the expected value of the long term yield $\mathbb{E} h_w(Y_\infty)$ is equal to the deterministic yield of the same strategy $\overline G (w)$. The variance of the yield $h_w(Y_\infty)$ is given by
\[
\sigma^2(h_w(Y_\infty)):=w^2 \mathbb{E}[f^2(w,\xi_1)- \overline f(w)^2].
\]
\end{prop}
\begin{proof}
Suppose $w$ is self-sustaining. Then
\[
\mathbb{E} h_w(Y_\infty)=\mathbb{E}(Y_\infty-w) = \mathbb{E} w f(w,\xi_1)-w=\overline F(w)-w=\overline G(w).
\]
For the variance of the yield we get
\begin{equation*}
\begin{split}
  \sigma^2(h_w(Y_\infty))&= \mathbb{E}[h_w(Y_\infty)^2] - (\mathbb{E} h(Y_\infty))^2\\
  &= \mathbb{E} [(wf(w,\xi_1)-w)^2] - (\overline F(w)-w)^2\\
  &=\mathbb{E} w^2 f^2(w,\xi_1) -2 w^2 \mathbb{E} f(w,\xi_1) + w^2 -\overline F^2(w) + 2\overline F(w) w - w^2\\
  &=w^2 \mathbb{E}[f^2(w,\xi_1)-(\mathbb{E} f(w,\xi_1))^2].
\end{split}
\end{equation*}
This completes the proof.
\end{proof}
\textit{\textbf{Biological interpretation:} Self-sustaining thresholds $w$, where the fluctuations never push the population's size under this threshold, make it possible to have the same yield in the stochastic and deterministic settings. However, the environmental fluctuations make the variance of the yield increase.
}
\begin{prop}
Suppose the following properties hold:
\begin{itemize}
  \item $w$ is not self-sustaining
  \item all the levels $x\in [0,r_s]$ are self-sustaining
  \item $\overline G(x)=x\overline F(x) - x$ is unimodal with its maximum at $\overline x$
  \item $\overline x$ is self-sustaining.
\end{itemize}
The expected steady state yield of the harvesting strategy $h_w$ is strictly greater than the deterministic nominal yield of the same strategy
\[
\mathbb{E} h_w(Y_\infty)>\overline G(w).
\]
\end{prop}
\begin{proof}
If $h_w$ is given by \eqref{e:bang} for some $w>0$ then
\begin{equation}\label{e:bang2}
   u_w(y) := \begin{cases}
             w & \text{if } y>w,\\
               y & \text{if } y\leq w.
          \end{cases}
\end{equation}
Note that $Y_\infty$ will be supported by a subset of $[r_s,\infty)$. By assumption $w>r_s>\overline x$, so that
\[
\overline G(y) > \overline G(w), y\in [r_s,w].
\]
This implies that with probability one
\[
\overline G(Y_\infty) \mathbf{1}\{Y_\infty\in [r_s,w]\}>\overline G(w)\mathbf{1}\{Y_\infty\in [r_s,w]\}.
\]
Using that the function $\overline G$ is nonincreasing on $[r_s,w]$ together with \eqref{e:bang2} and the last inequality we see that
\begin{equation*}
\begin{split}
\mathbb{E} h_w(Y_\infty) &= \mathbb{E} \overline G(u_w(Y_\infty))\\
&=\mathbb{E} [\overline G(Y_\infty) \mathbf{1}\{Y_\infty\in [r_s,w]\} ] + \mathbb{E}[\overline G(w) \mathbf{1}\{Y_\infty> w\}]\\
&>  \mathbb{E} [\overline G(w) \mathbf{1}\{Y_\infty\in [r_s,w]\} ] + \mathbb{E}[\overline G(w) \mathbf{1}\{Y_\infty> w\}]\\
&= \overline G(w).
\end{split}
\end{equation*}
with $w>0$.
\end{proof}
\textit{\textbf{Biological interpretation:} Suppose one picks a harvesting threshold $w$ which is not self-sustaining while the maximum yield of the deterministic dynamics happens at a threshold $\overline x<w$ which is self-sustaining. Then the expected yield of the constant escapement strategy with threshold $w$ for the stochastic dynamics is strictly greater than the expected yield of the same strategy in the deterministic system. The environmental fluctuations will push the population size into the region $(\overline x, w)$ and in this region, the function $\overline G$, which measures the size of the deterministic harvest, is strictly decreasing. This makes it more favorable to go below $w$, something which is not possible in the deterministic dynamics.
}

\section{The Ricker model: single species}
In this section we will provide an in depth analysis of the \textit{Ricker model}.
Its dynamics is given by the functional response
\[
f(x,\xi) = e^{\rho-\alpha x}.
\]
Here the randomness comes from $\xi:=(\rho,\alpha)$. The quantity $\rho_t$ is the fluctuating growth rate and $\alpha_t$ is the competition rate. We assume that $\rho_1, \rho_2,\dots$ are i.i.d random variables on $\mathbb{R}$, and $\alpha_1,\dots$ are independent i.i.d random variables supported on $\mathbb{R}$.
In this setting, one can see that without harvesting
\[
r^X(\delta_0) = \mathbb{E} \rho_1
\]
while with harvesting strategy $h(y)$ (or escapement strategy $u(y)$)
\[
r^Y(\delta_0)= \ln\left(\frac{\partial u}{\partial x}(0) h\right) + \mathbb{E} \rho_1.
\]
The maximal harvesting rate at $0$ which does not lead to extinction is
\[
\frac{\partial h}{\partial x}(0)<1- e^{-  \mathbb{E}\rho_1}.
\]
\subsection{Maximum sustainable yield.} We want to see when we can apply the results of Theorem \ref{t:maximal_yield}. Suppose that $\rho_1$ is such that $\mathbb{E} e^{\rho_1} = K_1>0$ and assume for simplicity that $\alpha_1>0$ is a constant. Then
\[
\overline f (x) = \mathbb{E} e^{\rho_1-\alpha_1 x} =  K_1 e^{-\alpha_1 x},
\]
$\overline F(x) = x \overline f (x)= xK_1 e^{-\alpha_1 x}$ and $\overline G(x) = x K_1 e^{-\alpha_1 x} - x$. By the analysis from Section \ref{s:stocdet} the deterministic maximum yield is achieved at the point $x_1$ where
\[
\overline G'(x_1)=\overline f(x_1) +x_1 \overline f'(x_1) -1 = K_1e^{-\alpha_1 x_1}-\alpha_1 x_1  K_1e^{-\alpha_1 x_1} - 1=0.
\]
Define the function
\[
q(x) = K_1 e^{-\alpha_1 x} - \alpha_1 K_1x e^{-\alpha_1 x} - 1, x\in \mathbb{R}_+.
\]
\begin{lm}
If $q(0)<0$ then the equation $q(x)=0$ has no solutions on $(0,\infty)$. If instead $q(0)>0$ then the equation $q(x)=0$ has exactly one solution $x_1>0$.
\end{lm}
\begin{proof}
Note that
\[
q'(x)= \alpha_1K_1 e^{-\alpha_1 x}(-2+\alpha_1x)
\]
and
\[
q''(x)=-\alpha_1^2 K_1  e^{-\alpha_1 x} (-2+\alpha_1 x) + \alpha_1^2 K_1 e^{-\alpha_1 x}.
\]
This shows that starting from $x=0$, the function $q$ decreases to its minimum at $x=\frac{2}{\alpha_1}$ and then increases from there on forever. However, once $q$ goes below zero it will never go above zero again. This happens because of the above properties and the fact that
\[
\lim_{x\to\infty} q(x) =\lim_{x\to\infty} (K_1 e^{-\alpha_1 x} - \alpha_1 K_1x e^{-\alpha_1 x} - 1)=-1.
\]
This implies that if $q(0)<0$ there are no solutions to $q(x)=0$. If we assume $q(0)>0$ we get in combination with $\lim_{x\to \infty }<0$ by the Intermediate Value Theorem that there exists a solution to $q(x)=0$. It is also clear by the properties of $q(x)$ that there exists exactly one solution to $q(x)=0$ and the solution has to lie in the interval $\left(0, \frac{2}{\alpha_1}\right)$.
\end{proof}
In order to be able to achieve this yield in the stochastic setting, according to Theorem \ref{t:maximal_yield} we need to ensure that $x_1$ is self-sustainable. This boils down to
\[
 \mathbb{P}\left(f(x_1,\xi)\geq 1\right) =\mathbb{P}\left(e^{\rho_1-\alpha_1 x_1}\geq 1\right)=1,
\]
or
\[
\mathbb{P}\left(\rho_1\geq \alpha_1 x_1 \right)=1.
\]
Since $x_1\in \left(0, \frac{2}{\alpha_1}\right)$ we see that if $\rho_1\geq 2$ with probability one, then the self-sustaining harvesting policy given by
\begin{equation*}
   h^*(y) := \begin{cases}
             y-x_1 & \text{if } y>x_1,\\
               0 & \text{if } y\leq x_1,
          \end{cases}
\end{equation*}
where $x_1$ is the unique solution to $q(x)=0$, maximizes the expected long term yield and makes it equal to the deterministic maximal sustainable yield. The value of the optimal expected long term yield will be
\[
\mathbb{E} h^*(Y_\infty) = \overline G(x_1) = x_1(K_1 e^{-\alpha_1 x_1}-1).
\]
\subsection{Maximal constant effort policy.} Suppose we use a constant effort policy $h(x)=\theta x$ for some $\theta\in(0,1)$ and that both $\rho_1$ and $\alpha_1$ are random.
The condition for persistence (see Theorem \ref{t:pers1_disc} and Theorem \ref{t:pers_single3}) is given by
\[
\mathbb{E}\rho_1 + \ln(1-\theta)> 0.
\]
This forces that $\theta\in \left(0,1-e^{-\mathbb{E}\rho_1}\right)$. Assume this condition holds so that $Y_t$ converges to a stationary distribution $\pi_\theta$. Then \eqref{e:r0} becomes
\[
0 = \ln(1-\theta) +\mathbb{E}\rho_1 - (1-\theta) \mathbb{E}\alpha_1 \int_{(0,\infty)}x\,\pi_\theta(dx).
\]
We can use this to show that the long run expected yield is given by
\[
H(\theta):=\mathbb{E} h(Y_\infty) = \int_{(0,\infty)}h(x)\,\pi_\theta(dx) = \frac{\theta (\mathbb{E}\rho_1 +\ln(1-\theta))}{(1-\theta) \mathbb{E}\alpha_1}.
\]

 \begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2in]{Hrho1alpha1234.pdf}
 		\caption{Graph of the long run average yield $H(\cdot)$ as a function of the harvesting rate $\theta$ when $\mathbb{E} \rho_1=1$ and $\mathbb{E}\alpha_1=1, 2, 3, 4$.} \label{fig1}\end{center}
 \end{figure}

 The Intermediate Value Theorem shows there is a solution $\theta^*\in\left(0,1-e^{-\mathbb{E}\rho_1}\right)$ to
\begin{equation}\label{e:H'}
0=H'(\theta^*)=\frac{\mathbb{E}\rho_1-\theta^* + \ln(1-\theta^*)}{ \mathbb{E}\alpha_1(1-\theta^*)^2}.
\end{equation}
Since the function $p(x)=\ln(1-x)-x$ is strictly decreasing on $(0,1)$ we also get that the solution $x^*$ is unique.
Taking another derivative, evaluating at $x^*$ and using \eqref{e:H'} we get
\[
H''(\theta^*)= \frac{2\mathbb{E}\rho_1-\theta^*+ 2\ln(1-\theta^*)-2}{\mathbb{E}\alpha_1(1-\theta^*)^3} = \frac{\theta^*-2}{\mathbb{E}\alpha_1(1-\theta^*)^3}<0.
\]
This implies that $\theta^*$ is a global maximum of $H(\theta)$ on $\left[0,1-e^{-\mathbb{E}\rho_1}\right]$. The maximal expected constant effort harvesting yield will be
\[
H(\theta^*) =  \frac{(\theta^*)^2}{(1-\theta^*) \mathbb{E}\alpha_1}.
\]

\section{Harvesting of two interacting species}
In this section we analyze the situation when there are two interacting species that can be harvested.
The system is modelled in the absence of harvesting by
\begin{equation*}
\begin{aligned}
X^1_{t+1}&=X^1_{t} f_1(X^1_t, X^2_t,\xi_{t+1}),\\
X^2_{t+1}&=X^1_{t} f_2(X^1_t, X^2_t,\xi_{t+1}).
\end{aligned}
\end{equation*}
As the theory from Appendix \ref{s:a2d} shows, one needs to first look at the quantities
\[
r_i(\delta_0) = \mathbb{E}[\ln f_i(0,\xi_1)], i=1,2.
\]
If $r_i(\delta_0)>0$ then species $i$ survives on its own and converges to a unique invariant probability measure $\mu_i$ supported on $(0,\infty)$.
Suppose $r_1(\delta_0)>0, r_2(\delta_0)>0$. The realized per-capita growth rates can be computed via
  \[
  r_i(\mu_j)=\int_{(0,\infty)}\mathbb{E}[\ln f_i(x,\xi_1)]\mu_j(dx).
  \]
If $r_1(\mu_2)>0$ and $r_2(\mu_1)>0$ by Theorem \ref{t:pers1_disc} we have the convergence to a unique stationary distribution $\pi$ supported on $\mathcal{S}_+$.
\subsection{Two species with harvesting.} Assume next that we harvest according to the strategies $h_1(x_1,x_2)$ and $h_2(x_1, x_2)$. Using \eqref{e:g} and \eqref{e:harv3} the dynamics becomes
\begin{equation}\label{e:harv2d}
Y^i_{t+1} = Y^i_t  g_i(\mathbf{Y}_t, \xi_{t+1}).
\end{equation}
where for $i=1,2$
\begin{equation}\label{e:g2}
   g_i(y_1,y_2,\xi) := \begin{cases}
              \frac{u_i(y_1,y_2)}{y_i} f_i(u_1(y_1,y_2), u_2(y_1,y_2), \xi)& \text{if } y_i>0,\\
               \left(\frac{\partial u_i}{\partial x_i}(y_1,y_2)\right)f_i(u_1(y_1,y_2), u_2(y_1,y_2), \xi) & \text{if } y_i = 0.
          \end{cases}
\end{equation}
Species $Y^i$ persists on its own with harvesting if
\begin{equation}\label{e:2d_pers2}
\begin{aligned}
r^Y_i(\delta_0) &= \mathbb{E}[\ln g_i(0,0,\xi_1)]\\
&= \mathbb{E}\ln\left[\left(\frac{\partial u_i}{\partial x_i}(0,0)\right)f_i(u_1(0,0), u_2(0,0), \xi_1)\right]\\
&= \ln \left(\frac{\partial u_i}{\partial x_i}(0,0)\right) + \mathbb{E}\ln f_i(0,0, \xi_1)\\
&>0.
\end{aligned}
\end{equation}
or equivalently
\begin{equation}\label{e:r_two}
\frac{\partial u_i}{\partial x_i}(0,0)> e^{-  \mathbb{E}\ln f_i(0,0, \xi_1)}.
\end{equation}
At this point there are three possibilities one may want to look at:
\begin{itemize}
  \item [1)] $\mathbb{E}\ln f_i(0,0, \xi_1)>0, i=1,2$ and $\frac{\partial u_i}{\partial x_i}(0,0)> e^{-  \mathbb{E}\ln f_i(0,0, \xi_1)}, i=1,2$ so that both harvested species persist on their own and have unique invariant probability measures $ \mu_1$ and $ \mu_2$ on the two positive axes. This describes the harvesting of a competitive system.
  \item [2)]$\mathbb{E}\ln f_i(0,0, \xi_1)>0, i=1,2$, $\frac{\partial u_1}{\partial x_1}(0,0)> e^{-  \mathbb{E}\ln f_1(0,0, \xi_1)}$ and $\frac{\partial u_2}{\partial x_2}(0,0)< e^{-  \mathbb{E}\ln f_2(0,0,\xi_1)}$. In this case, there are two species which compete with each other, both species persist on their own when there is no harvesting, and species $1$ persists with harvesting on its own, while species $2$ goes extinct if it is on its own and gets harvested.
  \item [3)] $\mathbb{E}\ln f_1(0,0, \xi_1)>0$, $\mathbb{E}\ln f_2(0,0, \xi_1)<0$, and $\frac{\partial u_1}{\partial x_1}(0,0)> e^{-  \mathbb{E}\ln f_1(0,0, \xi_1)}$. In this setting, species $1$ is a prey that persists on its own both with harvesting and without harvesting while species $2$ is a predator that can not persist on its own.
\end{itemize}
\begin{example}

Assume we work with constant threshold harvesting strategies, so that we harvest species $1$ according to
\begin{equation*}
   h(y_1,y_2) := \begin{cases}
             y_1-w & \text{if } y_1>w,\\
               0 & \text{if } y_1\leq w
          \end{cases}
\end{equation*}
where $w>0$. We will suppose species $2$ does not get harvested.
As we have seen in Section \ref{s:stocdet}, constant escapement harvesting strategies do not influence the persistence of a single species (as long as the threshold is strictly positive). However, we can show that they do change the persistence criteria if there are two interacting species.

Suppose species $1$ persists one its own: $\mathbb{E}\ln f_1(0,0, \xi_1)>0$. Without harvesting it converges to a stationary distribution $\tilde \mu_1$ while with harvesting it converges to a different stationary distribution $ \mu^{h}_1$.
Without harvesting we have
\[
r_2(\mu_1) = \int_{(0,\infty)}\mathbb{E}[\ln f_2(x,0,\xi_1)]\tilde \mu_1(dx)
\]
while with harvesting
\begin{equation}\label{e:2d_1}
\begin{aligned}
r_2(\mu_1^{h}) &= \int_{(0,\infty)}\mathbb{E}[\ln g_2(x,0,\xi_1)]\mu_1^{h}(dx) \\
&=  \int_{(0,\infty)}\mathbb{E}[\ln f_2(u_1(x_1,0),0,\xi_1)]\mu_1^{h}(dx_1)\\
&=\int_{(0,w)}\mathbb{E}[\ln f_2(x_1,0,\xi_1)]\mu_1^{h}(dx_1) + \int_{(w,\infty)}\mathbb{E}[\ln f_2(w,0,\xi_1)]\mu_1^{h}(dx_1).
\end{aligned}
\end{equation}
We see that in general, since $\mu_1^{h}\neq \tilde \mu_1$, we will have $r_2(\mu_1)\neq r_2(\mu_1^{h})$.  Therefore the persistence criteria are influenced by the threshold policies.
\end{example}

\subsection{Two dimensional  Lotka--Volterra predator-prey model.}
Suppose that we have model with a predator and a prey, that get harvested proportionally at rates $q,r\in (0,1)$. This means that $u_1(y_1,y_2)=(1-q)y_1$ and $u_2(y_1,y_2)=(1-r)y_2$. Using this in \eqref{e:harv2d} we get
\begin{eqnarray*}
Y^1_{t+1} &=& Y^1_t \exp \left(\rho_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)Y^1_t-a_{t+1}(1-r)Y^2_t\right)\\
Y^2_{t+1}&=& Y^2_t \exp \left(-d_{t+1}+\ln(1-r)-c_{t+1}(1-r)Y^2_t+b_{t+1}(1-q)Y^1_t\right)
\end{eqnarray*}
where the random coefficients have the following interpretations: $d_{t+1}>0$ is the predator's death rate, $a_{t+1}>0$ is the predator's attack rate on the prey, $b_{t+1}>0$ is the predator's conversion rate of prey, $c_{t+1}>0$ is the predator's intraspecific competition rate. Let $d_r:=\mathbb{E}[d_1]-\ln(1-r)>0, \rho_q:=\mathbb{E}[\rho_1]+\ln(1-q)$ and assume $(\rho_t)_{t\in \mathbb{Z}_+}, (\alpha_t)_{t\in \mathbb{Z}_+}, (a_t)_{t\in \mathbb{Z}_+}, (d_t)_{t\in \mathbb{Z}_+}, (c_t)_{t\in \mathbb{Z}_+}$ and $(b_t)_{t\in \mathbb{Z}_+}$ form independent sequences of i.i.d random variables. We assume for simplicity that the different random variables have compact support and are absolutely continuous with respect to Lebesgue measure. Then one can show by \cite{H87, BS19} that there is $K>0$ such that the process $\mathbf{Y}_t$ eventually enters and then stays forever in the compact set $\mathcal{K}=[0,K]^2$. We assume that the boundaries $\{0\}\times (0,K) $ and $(0,K)\times \{0\}$ are accessible for all $q,r\in(0,1)$. This can be easily checked for example if $\rho_1$ has $(0,L)$ for some $L>0$ in its support.

 As long as $\rho_q>0$, by our previous results, there exists a unique stationary distribution $\mu^q_1$ on $(0,\infty)\times \{0\}$ and
$$
 \int x_1\,\mu^q_1(dx_1) =\frac{\rho_q}{(1-q)\mathbb{E} \alpha_1} = \frac{\mathbb{E}[\rho_1]+\ln(1-q)}{(1-q)\mathbb{E} \alpha_1}.
$$
This can be used to get
$$
r_2(\mu^q_1) = -d_r +\mathbb{E}[b_1] \frac{\rho_q}{\mathbb{E} \alpha_1}= -(\mathbb{E}[d_1]-\ln(1-r)) +\mathbb{E}[b_1] \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1}.
$$
If $r_2(\mu_1^q)<0$ the predator $Y^2$ will go extinct with probability one. This means that, for a given harvesting rate $q\in \left(0, 1 - e^{-\mathbb{E}\rho_1}\right)$ of the prey, the maximal harvesting rate of the predator that does not lead to its extinction is
\[
r_{\max} = 1- \exp\left(\mathbb{E} d_1 - \mathbb{E}[b_1] \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1} \right).
\]

As long as $r_2(\mu_1^q)>0$, or equivalently
\[
\frac{d_r}{\rho_q}<  \frac{\mathbb{E} b_1}{\mathbb{E} \alpha_1},
\]
we get the existence of a unique invariant probability measure $\mu^{q,r}_{12}$ supported on a subset of $(0,\infty)^2$. Putting all the conditions together we get the following classification of the harvested dynamics:
\begin{itemize}
  \item If
   \begin{equation*}
\begin{split}
0&<q< 1 - e^{-\mathbb{E}\rho_1}\\
0&<r<r_{\max}=1- \exp\left(\mathbb{E} d_1 - \mathbb{E}[b_1] \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1} \right)
\end{split}
\end{equation*}
then the two species coexist and there is a unique invariant probability measure when
  \item If
   \begin{equation*}
\begin{split}
0&<q< 1 - e^{-\mathbb{E}\rho_1}\\
1&>r\geq r_{\max}=1- \exp\left(\mathbb{E} d_1 - \mathbb{E}[b_1] \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1} \right)
\end{split}
\end{equation*}
then the prey persists and the predator goes extinct with probability $1$.
  \item If $$1>q\geq  1 - e^{-\mathbb{E}\rho_1}$$
then both the prey and the predator go extinct with probability $1$.
\end{itemize}
We depict one example of the three possible regions in Figure \ref{fig2}.
\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{regions.pdf}
 		\caption{The regions of the harvesting rates $q,r$ for which both species persist, for which just the prey persists and for which both species go extinct. The parameters are $\mathbb{E} a_1=\mathbb{E} d_1= \mathbb{E} \alpha_1 = \mathbb{E} \rho_1 =1, \mathbb{E} b_1=2, \mathbb{E} c_1=1.5$.} \label{fig2}\end{center}
 \end{figure}
Suppose next that the two species persist. If we set $\overline Y^1 := \int_{(0,\infty)^2} x_1\mu^{q,r}_{12}(du), \overline Y^2 := \int_{(0,\infty)^2} x_2\mu^{q,r}_{12}(du)$, since by \eqref{e:equ} the per-capita growth rates at stationarity are zero, we get
\begin{eqnarray*}
r_1(\mu^{q,r}_{12})&=& 0 = \rho_q-\mathbb{E}[\alpha_{1}](1-q)\overline Y^1-\mathbb{E}[a_{1}](1-r)\overline Y^2\\
r_2(\mu^{q,r}_{12}) &=&0=-d_r-\mathbb{E}[c_{1}](1-r)\overline Y^2+(1-q)\mathbb{E}[b_{1}]\overline Y^1.
\end{eqnarray*}
Solving this linear system we get the unique solution
\begin{equation}\label{e:Y}
\begin{split}
\overline Y^1&= \frac{\mathbb{E}[c_1]\rho_q+\mathbb{E}[a_1]d_r}{1-q}\frac{1}{\mathbb{E}[\alpha_1]\mathbb{E}[c_1]+\mathbb{E}[a_1]\mathbb{E}[b_1]}\\
\overline Y^2&=\frac{\mathbb{E}[b_1]\rho_q-\mathbb{E}[\alpha_1]d_r}{1-r}\frac{1}{\mathbb{E}[\alpha_1]\mathbb{E}[c_1]+\mathbb{E}[a_1]\mathbb{E}[b_1]}
\end{split}
\end{equation}

On the other hand, if $r_2(\mu_1^q)<0$ we get by the results from Appendix \ref{s:a2d} that $\lim_{t\to\infty}\frac{\ln Y_t^2}{t}=r_2(\mu_1^q)<0$ so that the predator goes extinct exponentially fast.
\begin{prop}\label{p:1}
If $r_2(\mu_1^q)<0$ and $r_1(\delta_0)=\rho_q>0$ the prey will persist and $\mathbb{E} Y_t^1\to \frac{\rho_q}{(1-q)\mathbb{E}[\alpha_1]}$ as $t\to\infty$.
\end{prop}
\begin{proof}Pick $\varepsilon>0$ small. Since $Y^2_t<\varepsilon$ for $t>T$ large enough, we have with high probability
\begin{equation}\label{e:conv}
\begin{split}
   Y^1_t e^{ \rho_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)Y^1_t-a_{t+1}(1-r)\varepsilon} &<Y^1_{t+1} \\
   &= Y^1_t e^{ \rho_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)Y^1_t-a_{t+1}(1-r)Y^2_t}  \\
     &<Y^1_t e^{\rho_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)Y^1_t}.
     \end{split}
\end{equation}
Define the processes $(Y_t^\varepsilon)_{t\in\mathbb{Z}_+}$ and $(\tilde Y_t)_{t\in\mathbb{Z}_+}$ via
\[
Y_{t+1}^\varepsilon= Y^\varepsilon_t e^{ \rho_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)Y^\varepsilon_t-a_{t+1}(1-r)\varepsilon}
\]
and
\[
\tilde Y^1_{t+1}=\tilde Y^1_t e^{\rho_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)\tilde Y^1_t}.
\]
By the results from Appendix \ref{s:a1} it is easy to see that
\[
\lim_{t\to\infty} \mathbb{E} \tilde Y^1_{t+1}= \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1}
\]
and
\[
\lim_{t\to\infty} \mathbb{E} Y^\varepsilon_{t+1}= \frac{\mathbb{E} \rho_1 + \ln(1-q)-\varepsilon}{\mathbb{E} \alpha_1}.
\]
By \eqref{e:conv} we get
\[
\frac{\mathbb{E} \rho_1 + \ln(1-q)-\varepsilon}{\mathbb{E} \alpha_1}=\lim_{t\to\infty} \mathbb{E} \tilde Y^\varepsilon_{t+1}\leq \limsup_{t\to\infty} \mathbb{E} Y_t^1 \leq \lim_{t\to\infty} \mathbb{E} \tilde Y^1_{t+1}= \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1}
\]
and
\[
\frac{\mathbb{E} \rho_1 + \ln(1-q)-\varepsilon}{\mathbb{E} \alpha_1}=\lim_{t\to\infty} \mathbb{E} \tilde Y^\varepsilon_{t+1}\leq \liminf_{t\to\infty} \mathbb{E} Y_t^1 \leq \lim_{t\to\infty} \mathbb{E} \tilde Y^1_{t+1}= \frac{\mathbb{E} \rho_1 + \ln(1-q)}{\mathbb{E} \alpha_1}.
\]
Letting $\varepsilon\downarrow 0$ yields
\[
\lim_{t\to\infty}\mathbb{E} Y_t^1= \frac{\rho_q}{(1-q)\mathbb{E}[\alpha_1]}
\]
which finishes the proof.
\end{proof}

Suppose we want to find the optimal harvesting strategy. Since the profit from harvesting prey or predators might be different we let $\beta> 0$ represent the relative value of the predator compared to the prey.
The problem then becomes maximizing the function
$$
H(q,r):=\lim_{T\to \infty} \mathbb{E}[h_1( Y^1_T)+\beta h_2(Y^2_T)] = \lim_{T\to \infty} \mathbb{E}[q Y^1_T+\beta r Y^2_T]=\lim_{T\to \infty} \frac{\sum_{n=1}^T q Y^1_n+r\beta Y^2_n}{T} = q \overline Y^1 + r\beta \overline Y^2,
$$
for $q,r\in[0,1]^2$.  Using the expressions for $\overline Y^1$, $\overline Y^2$ from \eqref{e:Y} together with Proposition \ref{p:1}, with the understanding that we set $\overline Y^1=\frac{\rho_q}{(1-q)\mathbb{E}[\alpha_1]}, \overline Y^2=0$ if the prey persists and the predator goes extinct, and the domain regions identified above we get
$$
H(q,r) = \left\{
        \begin{array}{ll}
            \left(\frac{\mathbb{E}[c_1]q\rho_q+\mathbb{E}[a_1]qd_r}{1-q} + \frac{\beta\mathbb{E}[b_1]r\rho_q-\beta\mathbb{E}[\alpha_1]rd_r}{1-r}\right)\frac{1}{\mathbb{E}[\alpha_1]\mathbb{E}[c_1]+\mathbb{E}[a_1]\mathbb{E}[b_1]} & \quad  r<r_{\max}, 0<q< 1 - e^{-\mathbb{E}\rho_1} \\
           \frac{q\rho_q}{(1-q)\mathbb{E}[\alpha_1]} & \quad r\geq r_{\max}, 0<q< 1 - e^{-\mathbb{E}\rho_1}\\
           0&\quad q\geq 1 - e^{-\mathbb{E}\rho_1}.
        \end{array}
    \right.
$$

\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{Hgraph.pdf}
 		\caption{The graph of $H(q,r)$. The parameters are $\beta = \mathbb{E} a_1=\mathbb{E} d_1= \mathbb{E} \alpha_1 = \mathbb{E} \rho_1 =1, \mathbb{E} b_1=2, \mathbb{E} c_1=1.5$.} \label{fig3}\end{center}
 \end{figure}
\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{Hgraph2.pdf}
 		\caption{The graph of $H(q,r)$. The parameters are $\beta = 5, \mathbb{E} a_1=\mathbb{E} d_1= \mathbb{E} \alpha_1 = \mathbb{E} \rho_1 =1, \mathbb{E} b_1=2, \mathbb{E} c_1=1.5$.} \label{fig4}\end{center}
 \end{figure}

 \textit{\textbf{Biological interpretation:} It is never optimal to harvest both the predator and the prey. If the relative price $\beta$ of the predator compared to the prey is low, it is always optimal to harvest the predator to extinction (see Figure \ref{fig3}). This then lets the prey population increase, and one gains by harvesting the prey. If instead the relative price $\beta$ is high, it is optimal to never harvest the prey (see Figure \ref{fig4}). This leads to an increase in the predator population, which then increases the harvesting yield of the predators. }

\subsection{Two dimensional Lotka-Volterra competition mode.}
We look at a two-species discrete Lotka-Volterra competition model when the two species get harvested proportionally at rates $q,r\in (0,1)$. The harvested dynamics is given by
\begin{eqnarray*}
Y^1_{t+1} &=& Y^1_t \exp \left(\rho^1_{t+1}+\ln(1-q)-\alpha_{t+1}(1-q)Y^1_t-a_{t+1}(1-r)Y^2_t\right)\\
Y^2_{t+1}&=& Y^2_t \exp \left(\rho^2_{t+1}+\ln(1-r)-c_{t+1}(1-r)Y^2_t-b_{t+1}(1-q)Y^1_t\right)
\end{eqnarray*}
We set $\overline\rho^1_q:=\mathbb{E}[\rho^1_1]+\ln(1-q)$ and $\overline\rho^2_q:=\mathbb{E}[\rho^2_1]+\ln(1-r)>0$. We assume the different random coefficients are independent and form sequences of i.i.d random variables. Furthermore, we make the same assumptions that were made in the predator-prey system. These ensure that the state space is compact and that the boundaries are accessible.

 As long as $\overline\rho^1_q, \overline\rho^2_r>0$, by our previous results, there exists a unique stationary distribution $\mu^q_1$ (respectively $\mu^r_2$) on $(0,\infty)\times \{0\}$ (respectively $\{0\}\times (0,\infty) $) and
$$
 \int x_1\,\mu^q_1(dx_1) =\frac{\overline\rho^1_q}{(1-q)\mathbb{E} \alpha_1} = \frac{\mathbb{E}[\rho^1_1]+\ln(1-q)}{(1-q)\mathbb{E} \alpha_1},
$$
$$
 \int x_2\,\mu^r_2(dx_2) =\frac{\overline\rho^2_r}{(1-r)\mathbb{E} c_1} = \frac{\mathbb{E}[\rho^2_1]+\ln(1-r)}{(1-r)\mathbb{E} c_1}.
$$
One can then compute the per-capita growth rates
$$
r_1(\mu^r_2) = \overline \rho_q^1 - \mathbb{E} a_1 (1-r) \int x_2\,\mu^r_2(dx_2) = (\mathbb{E}[\rho^1_1]+\ln(1-q)) -\mathbb{E} a_1 \frac{\mathbb{E} \rho^2_1 + \ln(1-r)}{\mathbb{E} c_1},
$$
and
$$
r_2(\mu^q_1) = \overline \rho_r^2 - \mathbb{E} b_1 (1-q) \int x_1\,\mu^q_1(dx_1) = (\mathbb{E}[\rho^2_1]+\ln(1-r)) -\mathbb{E} b_1 \frac{\mathbb{E} \rho^1_1 + \ln(1-q)}{\mathbb{E} \alpha_1}.
$$

We get the following classification of the dynamics:
\begin{itemize}
  \item If
   \begin{equation*}
\begin{split}
0&<q< 1 - e^{-\mathbb{E}\rho^1_1}\\
0&<r< 1 - e^{-\mathbb{E}\rho^2_1}\\
 \frac{\mathbb{E} a_1}{\mathbb{E} c_1} &<  \frac{\mathbb{E} \rho^1_1 +\ln(1-q)}{\mathbb{E} \rho^2_1 + \ln(1-r)} < \frac{\mathbb{E} \alpha_1}{\mathbb{E} b_1}
\end{split}
\end{equation*}
then the two species coexist and the process converges to its unique invariant probability measure.
  \item If
  \begin{equation*}
\begin{split}
0&<q< 1 - e^{-\mathbb{E}\rho^1_1}\\
0&<r< 1 - e^{-\mathbb{E}\rho^2_1}\\
 \frac{\mathbb{E} a_1}{\mathbb{E} c_1} &<  \frac{\mathbb{E} \rho^1_1 +\ln(1-q)}{\mathbb{E} \rho^2_1 + \ln(1-r)}\\
  \frac{\mathbb{E} b_1}{\mathbb{E} \alpha_1} &>  \frac{\mathbb{E} \rho^2_1 +\ln(1-r)}{\mathbb{E} \rho^1_1 + \ln(1-q)}
\end{split}
\end{equation*}
then species 1 persists and species 2 goes extinct with probability $1$.
 \item If
   \begin{equation*}
\begin{split}
0&<q< 1 - e^{-\mathbb{E}\rho^1_1}\\
0&<r< 1 - e^{-\mathbb{E}\rho^2_1}\\
 \frac{\mathbb{E} a_1}{\mathbb{E} c_1} &>  \frac{\mathbb{E} \rho^1_1 +\ln(1-q)}{\mathbb{E} \rho^2_1 + \ln(1-r)}\\
  \frac{\mathbb{E} b_1}{\mathbb{E} \alpha_1} &<  \frac{\mathbb{E} \rho^2_1 +\ln(1-r)}{\mathbb{E} \rho^1_1 + \ln(1-q)}
\end{split}
\end{equation*}
then species 2 persists and species 1 goes extinct with probability $1$.
 \item If $(Y_0^1, Y_0^2)=(x,y)\in (0,\infty)^2$ then let $p_{x,y} = \mathbb{P}(Y_t^1\to \mu_1^q, Y_t^2\to 0~|~(Y_0^1, Y_0^2)=(x,y))$. If
\begin{equation*}
\begin{split}
0&<q< 1 - e^{-\mathbb{E}\rho^1_1}\\
0&<r< 1 - e^{-\mathbb{E}\rho^2_1}\\
 \frac{\mathbb{E} a_1}{\mathbb{E} c_1} &>  \frac{\mathbb{E} \rho^1_1 +\ln(1-q)}{\mathbb{E} \rho^2_1 + \ln(1-r)}\\
  \frac{\mathbb{E} b_1}{\mathbb{E} \alpha_1} &>  \frac{\mathbb{E} \rho^2_1 +\ln(1-r)}{\mathbb{E} \rho^1_1 + \ln(1-q)}
\end{split}
\end{equation*}
then we have bistability, that is, $p_{x,y}\in (0,1)$ and $1-p_{x,y}= \mathbb{P}(Y_t^2\to \mu_2^r, Y_t^1\to 0~|~(Y_0^1, Y_0^2)=(x,y))$.
\end{itemize}
Note that for a given set of coefficients one cannot have all the 4 regions if we vary $q$ and $r$. There are two possibilities, each having three regions. One is to have coexistence, the persistence of species 1 and extinction of species 2, or the extinction of species 1 and the persistence of species 2 (see Figure \ref{fig5}). The other possibility is to have bistability, the persistence of species 1 and extinction of species 2, or the extinction of species 1 and the persistence of species 2 (see Figure \ref{fig6}).

\textit{\textbf{Biological interpretation:} Harvesting can facilitate the coexistence of the two species. For example, suppose that
\[
 \frac{\mathbb{E} a_1}{\mathbb{E} c_1} < \frac{\mathbb{E} \alpha_1}{\mathbb{E} b_1}<  \frac{\mathbb{E} \rho^1_1 }{\mathbb{E} \rho^2_1 }.
\]
This implies that if there is no harvesting species $1$ persists and species $2$ goes extinct. It is clear that there exists $q\in (0,1)$ such that
\[
 \frac{\mathbb{E} a_1}{\mathbb{E} c_1}<  \frac{\mathbb{E} \rho^1_1 +\ln(1-q)}{\mathbb{E} \rho^2_1  } < \frac{\mathbb{E} \alpha_1}{\mathbb{E} b_1}
\]
which leads to coexistence. If one species has a competitive advantage so that without harvesting it drives the other competitor extinct, one can harvest this dominant species and get coexistence.
}
\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{nobistability.pdf}
 		\caption{Figure showing the regions of the harvesting rates $q,r$ for which both species persist, for which just the species 1 persists, and the region for which just species 2 persists. The parameters are $\mathbb{E} a_1=\mathbb{E} b_1=1, \mathbb{E} c_1=\mathbb{E} \alpha_1=1.5, \mathbb{E} \rho^1_1 =1, \mathbb{E} \rho^2_1 =1.5$.} \label{fig5}\end{center}
 \end{figure}

\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{nocoexistence.pdf}
 		\caption{Figure showing the regions of the harvesting rates $q,r$ for which there is bistability, for which just the species 1 persists, and the region for which just species 2 persists. The parameters are $\mathbb{E} a_1=1.5, \mathbb{E} b_1=2, \mathbb{E} c_1=1, \mathbb{E} \alpha_1=1, \mathbb{E} \rho^1_1 =1, \mathbb{E} \rho^2_1 =1.5$.} \label{fig6}\end{center}
 \end{figure}
If there is coexistence and the system converges to an invariant probability measure $\mu_{12}^{q, r}$ on $(0,\infty)^2$ we see by \eqref{e:equ} that the per-capita growth rates at stationarity are zero. This shows that
\begin{eqnarray*}
r_1(\mu_{12})=0&=& \overline \rho_q^1-\mathbb{E}[\alpha_{1}](1-q)\overline Y^1-\mathbb{E}[a_{1}](1-r)\overline Y^2\\
r_2(\mu_{12})=0&=&  \overline \rho_r^2-\mathbb{E}[c_{1}](1-r)\overline Y^2-(1-q)\mathbb{E}[b_{1}]\overline Y^1,
\end{eqnarray*}
where as before $\overline Y^1 = \int_{(0,\infty)^2} x_1\mu^{q,r}_{12}(du), \overline Y^2= \int_{(0,\infty)^2} x_2\mu^{q,r}_{12}(du)$ are the expected values of the two species at stationarity.
Solving this linear system yields the unique solution
\begin{eqnarray*}
\overline Y^1&=& \frac{\mathbb{E}[c_1]\overline \rho^1_q-\mathbb{E}[a_1]\overline\rho^2_r}{1-q}\frac{1}{\mathbb{E}[\alpha_1]\mathbb{E}[c_1]-\mathbb{E}[a_1]\mathbb{E}[b_1]}\\
\overline Y^2&=&\frac{\mathbb{E}[\alpha_1]\overline \rho^2_r-\mathbb{E}[b_1]\overline \rho^1_q}{1-r}\frac{1}{\mathbb{E}[\alpha_1]\mathbb{E}[c_1]-\mathbb{E}[a_1]\mathbb{E}[b_1]}.
\end{eqnarray*}
One can prove the following analogue of Proposition \ref{p:1}.
\begin{prop}\label{p:2}
If $r_2(\mu^q_1)<0$ and $r_1(\mu^r_2)>0$ then species 1 will persist and $\mathbb{E} Y_t^1\to \frac{\overline\rho^1_q}{(1-q)\mathbb{E}[\alpha_1]}$ as $t\to\infty$.
If $r_1(\mu^r_2) <0$ and $r_2(\mu^q_1)>0$ then species 2 will persist and $\mathbb{E} Y_t^2\to \frac{\overline \rho^2_r}{(1-r)\mathbb{E}[c_1]}$ as $t\to\infty$.
\end{prop}
Suppose the coefficients are such that the coexistence of the two species is possible. We are interested in maximizing the function
$$
H(q,r)=\lim_{T\to \infty} \mathbb{E}[q Y^1_T+\beta rY^2_T] =\lim_{T\to \infty} \frac{\sum_{n=1}^T q Y^1_n+r\beta Y^2_n}{T} = q \overline Y^1 + r\beta \overline Y^2,
$$
where $\beta> 0$ represents the relative value of species 2 compared to species 1. Using the expressions for $\overline Y^1$, $\overline Y^2$, together with Proposition \ref{p:2} and the domain regions identified above we get
$$
H(q,r) = \left\{
        \begin{array}{ll}
            \left(q\frac{\mathbb{E}[c_1]\overline \rho^1_q-\mathbb{E}[a_1]\overline\rho^2_r}{1-q} + \beta r\frac{\mathbb{E}[\alpha_1]\overline \rho^2_r-\mathbb{E}[b_1]\overline \rho^1_q}{1-r}\right)\frac{1}{\mathbb{E}[\alpha_1]\mathbb{E}[c_1]-\mathbb{E}[a_1]\mathbb{E}[b_1]} & \quad  r_2(\mu^q_1)>0, r_1(\mu^r_2)>0 \\
          q \frac{\overline\rho^1_q}{(1-q)\mathbb{E}[\alpha_1]} & \quad r_2(\mu^q_1)<0, r_1(\mu^r_2)>0\\
           \beta r\frac{\overline\rho^2_r}{(1-r)\mathbb{E}[c_1]} & \quad r_2(\mu^q_1)>0, r_1(\mu^r_2)<0\\
           0&\quad q\geq 1 - e^{-\mathbb{E}\rho^1_1}, r\geq 1 - e^{-\mathbb{E}\rho^2_1}.
        \end{array}
    \right.
$$

\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{Hcompetition.pdf}
 		\caption{The graph of $H(q,r)$. The parameters are $\beta=1, \mathbb{E} a_1=\mathbb{E} b_1=1, \mathbb{E} c_1=\mathbb{E} \alpha_1=1.5, \mathbb{E} \rho^1_1 =1, \mathbb{E} \rho^2_1 =1.5$} \label{fig7}\end{center}
 \end{figure}

\begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{Hcompetition2.pdf}
 		\caption{The graph of $H(q,r)$. The parameters are $\beta=1, \mathbb{E} a_1=\mathbb{E} b_1=1, \mathbb{E} c_1=\mathbb{E} \alpha_1=1.5, \mathbb{E} \rho^1_1 =1.4, \mathbb{E} \rho^2_1 =1.5$} \label{fig8}\end{center}
 \end{figure}

 \begin{figure}[h!tb]
 	\begin{center}
 		\includegraphics[height=2.8in]{Hcompetition3.pdf}
 		\caption{The graph of $H(q,r)$. The parameters are $\beta=1, \mathbb{E} a_1=\mathbb{E} b_1=1.4, \mathbb{E} c_1=\mathbb{E} \alpha_1=1.5, \mathbb{E} \rho^1_1 =1.4, \mathbb{E} \rho^2_1 =1.5$} \label{fig9}\end{center}
 \end{figure}
\textit{\textbf{Biological interpretation:} Depending on the interaction coefficients, growth rates, and the relative value of the species there are three possible scenarios for the optimal harvesting strategy. In one case we harvest species $1$ to extinction and maximize the yield from harvesting species $2$. In other instances it is best to harvest species $2$ to extinction and maximize the harvest from species $1$. The third instance is the one of coexistence: the optimal harvesting strategy is to keep both species alive. In Figure \ref{fig7} we can see, that since the growth rate of species $2$ is greater than that of species $1$, while the other coefficients are identical, it is optimal to harvest species $1$ to extinction and to get a higher harvesting yield from species $2$. In the example from Figure \ref{fig8}, when the species are similar to each other, it is optimal to keep both species alive. However, once we increase the competition, it becomes optimal to drive one species extinct through harvesting (see Figure \ref{fig9}). These examples show that there is a delicate balance one has to take into account when looking for the optimal harvesting strategies. The intra- and interspecific competition rates, growth rates, and the prices of the species turn out to play key roles.
}

\textbf{Acknowledgements:} The author thanks Dang Nguyen and Sergiu Ungureanu for helpful discussions related to the paper.
\bibliographystyle{agsm}