Optimal Surplus-Dependent Reinsurance under Regime-Switching in a Brownian Risk Model

Abstract

In this paper, we consider a company that wishes to determine the optimal reinsurance strategy minimising the total expected discounted amount of capital injections needed to prevent the ruin. The company’s surplus process is assumed to follow a Brownian motion with drift, and the reinsurance price is modelled by a continuous-time Markov chain with two states. The presence of regime-switching substantially complicates the optimal reinsurance problem, as the surplus-independent strategies turn out to be suboptimal. We develop a recursive approach that allows to represent a solution to the corresponding Hamilton–Jacobi–Bellman (HJB) equation and the corresponding reinsurance strategy as the unique limits of the sequence of solutions to ordinary differential equations and their first- and second-order derivatives. Via Ito’s formula, we prove the constructed function to be the value function. Two examples illustrate the recursive procedure along with a numerical approach yielding the direct solution to the HJB equation.

1. Introduction

Writing red numbers is generally considered a bad sign for the financial health of a (insurance) company. An old but also highly criticised concept to measure company’s riskiness is the ruin probability, i.e., the probability that the company’s surplus will become negative in finite time. There is a vast of literature on ruin probability in different settings and under various assumptions, see, for instance, Rolski et al. (1999); Asmussen and Albrecher (2010) and further reference therein.

As ruin probabilities do not take into account the time and the severity of ruin, the related concept of capital injections incorporating both features has been suggested in Pafumi (1998) in the discussion of Gerber and Shiu (1998). The risk is measured by the expected discounted amount of capital injections needed to keep the surplus non-negative. If the discounting rate, or rather the preference rate of the insurer, is positive, then the amount of capital injections is minimised if one injects just as much as it is necessary to shift the surplus to zero (but not above) and to inject just when the surplus becomes negative (but not before, anticipating a possible ruin), see, for instance, Eisenberg and Schmidli (2009).

A well-established way to reduce the risk of an insurance portfolio is to buy reinsurance. Finding the optimal or fair reinsurance in different settings is a popular and widely investigated topic in insurance mathematics, see, for instance, Azcue and Muler (2005); Ben Salah and Garrido (2018); Brachetta and Ceci (2020) and the references therein. However, the reinsurance premia are usually higher than the premia of the first insurer. Otherwise, it would create an arbitrage opportunity for the first insurer, who could transfer the entire risk to the reinsurer (i.e., the amount of the necessary capital injection is zero) and still receive a risk-free gain in from of the remaining premium payments. If we consider a model including both capital injections and reinsurance, the capital injection process will naturally depend on the chosen reinsurance strategy. That is, one can control the capital injections—representing the company’s riskiness—by reinsurance. In this context, the problem of finding a reinsurance strategy leading to a minimal possible value of expected discounted capital injections has been solved in Eisenberg and Schmidli (2009). There, the optimal reinsurance strategy is given by a constant, meaning that the insurance company is choosing a retention level once and forever. This result has been obtained under the assumption that the parameters describing the evolution of the insurer’s and reinsurer’s surplus never change. However, the reality offers a contrary picture. The state of economy has an enormous impact on the insurance/reinsurance companies, adding an exogeneously given source of uncertainty.

In financial literature, regime-switching models have become very popular because they take into account possible macroeconomic changes. Originally proposed by Hamilton to model stock returns, this class of models has been adopted also in insurance mathematics, see, e.g., Asmussen (1989); Bargès et al. (2013); Bäuerle and Kötter (2007); Gajek and Rudź (2020); Jiang and Pistorius (2012). In this connection, one should not forget the models containing hidden information. Reinsurance companies deciding over the price of their reinsurance products have to take into account the competition on the market and the consequences of adverse selection, see, for instance, Chiappori et al. (2006) and references therein.

In the present paper, we model the surplus of the first insurer by a Brownian motion with drift. The insurer is obligated to inject capital if the surplus becomes negative and is allowed to buy proportional reinsurance. In order to account for the macroeconomic changes that are assumed to happen in circles, we allow the price of the reinsurance—represented through a safety loading—to depend on the current regime of the economy. A continuous-time Markov chain with two states describes the length of the regimes and the switching intensity from one state into another. We target to find a reinsurance strategy that minimises the value of expected discounted capital injections where the discounting rate is a positive regime-independent constant. If the discounting rate would be assumed to be negative in one of the states, it might become optimal to inject capital even if the surplus is still positive (see, for instance, in Eisenberg and Krühner (2018)), which would substantially complicate the problem. For the same reason, we do not incorporate hidden information or moral hazard into this model

We solve the optimisation problem via the Hamilton–Jacobi–Bellman (HJB) equation, which is in this case a system of equations. Differently than it has been done in the one-regime case, we cannot guess the optimal strategy and prove the corresponding return function to solve the HJB equation. Instead, we solve the system of HJB equations recursively. First of all, the system of HJB equations is rewritten as a system of ordinary differential equations. Then, we assume that the value function, say in the second regime, is an exponential function and solve the corresponding ordinary differential equation for the first regime. The obtained solution is inserted into the ordinary differential equation for the second regime. Proceeding in this way, we obtain a monotone uniformly converging sequence of solutions, whose limit functions solve the original HJB equation. Here, it is of crucial importance to choose correctly the exponential of the starting function in the recursion. We present an equation system providing the only correct choice of the starting function.

The aim of the present paper is to develop an algorithm for finding a candidate for the value function. Like in the case with just one regime, the HJB equation is rewritten as a differential equation with boundary conditions. Here, we are facing a boundary value problem, i.e., the conditions are specified at different boundaries, with one boundary being even infinity. Therefore, using Volterra-type representations and comparison theorems, we prove the existence and uniqueness of a solution to the HJB equation. Ito’s formula allows one to show that the constructed solution is indeed the value function. We show that the optimal reinsurance strategies are increasing in one regime and decreasing in the other, depending on the parameters. This fact reflects the dependence of the strategies on the reinsurance prices along with the switching intensities. For instance, being in a regime with a low reinsurance price and a relative high switching intensity into a state with a high reinsurance price would produce a decreasing proportion of the self-insured risk.

As we do not get a closed form expressions for the value function and for the optimal strategies, we give a numerical illustration of both the algorithm and the value function. Here, we follow the approach of Auzinger et al. (2019).

The remainder of the paper is structured as follows. In Section 2, we give a mathematical formulation of the problem and present the Hamilton–Jacobi–Bellman equation. In Section 3, we shortly revise the case of constant controls and prove that those cannot be optimal except for the case when it is optimal to buy no reinsurance. In Section 4 and Section 5 we recursively construct a function solving the HJB equation and prove it to be the value function. Finally, we explore the problem numerically in Section 6 and conclude in Section 7.

2. The Model

In the following, we give a mathematical formulation of the problem and present the heuristically derived Hamilton–Jacobi–Bellman (HJB) equation. We are acting on a probability space $(\Omega ,\mathcal{F},\mathbb{P})$.

In the classical risk model, the surplus process of an insurance company is given by

$$X_t=x+ct-\sum _{i=1}^{N_t}{Z}_i,$$

where $\left\{N_t\right\}$ is a Poisson process with intensity $\lambda$ and the claim sizes $Z_i$ are iid. with $\mathbb{E}\left[Z_1\right]=\mu$ and $\mathbb{E}\left[Z_1^2\right]={\mu }_2$ and independent of $\left\{N_t\right\}$. Furthermore, x denotes the initial capital and $c>0$ the premium rate. For further details on the classical risk model, see, e.g., Chapter 5.1 in Schmidli (2017).

The insurer can buy proportional reinsurance at a retention level $0\le b\le 1$, i.e., for a claim $Z_i$, the cedent pays $b{Z}_i$, and the reinsurer pays the remaining claim $(1-b)Z_i$. Assume the expected value principle for the calculation of the insurance and reinsurance premia with safety loadings $\eta >0$ and $\theta >0$, respectively, where $\eta <\theta$ (see Chapter 1.10 in Schmidli (2017)) transforms the surplus of the insurer, denoted now by $X^b$ to

$$\begin{array}{c}\hfill X^b\left(t\right)=x+c\left(b\right)t-\sum _{i=1}^{N_t}b{Z}_i.\end{array}$$

(1)

The new premium rate depends on the retention level and is given by $c\left(b\right)=\left(b\right(1+\theta )-(\theta -\eta \left)\right)\lambda \mu$, being old premia reduced for the premia paid to the reinsurer (see, e.g., Chapter 5.7 in Schmidli (2017)).

Usually, optimisation problems can be tackled more easily if the surplus is given by a Brownian motion. Therefore, diffusion approximation of the classical risk model is a popular concept in optimisation problems in insurance. A diffusion approximation to (1) by adopting a dynamic reinsurance strategy $B=\left\{b_t\right\}$, that is, the retention level $b_t$ changes continuously in time, is given by

$$\begin{array}{c}\hfill X_t^B=x+\theta {\int }_0^t{b}_s\mathrm{d}s-\lambda \mu (\theta -\eta )t+\sqrt{\lambda {\mu }_2}{\int }_0^t{b}_s\mathrm{d}{W}_s,\end{array}$$

(2)

dfsuch that the first two moments of (1) and (2) remain the same; see, for instance, Appendix D.3 in Schmidli (2008) for details. In addition to buying reinsurance, the insurance company has to inject capital in order to keep the surplus non-negative. The process describing the accumulated capital injections up to time t under a reinsurance strategy B will be denoted by $Y^B=\left\{Y_t^B\right\}$. The surplus process under a reinsurance strategy B and capital injections Y is given by

$$\begin{array}{c}\hfill X_t^{B,Y}=X_t^B+Y_t.\end{array}$$

(3)

Further, we introduce a continuous-time Markov chain $J=\left\{J_t\right\}$ with state space $\mathcal{S}=\{1,2\}$. We assume that J and W are independent, and that J has a strongly irreducible generator $Q={\left[q_{ij}\right]}_{2\times 2}$, where $q_{ij}=-q_{ii}$ for $i\ne j$, and we consider the filtration $\left\{{\mathcal{F}}_t\right\}$, generated by the pair $(W,J)$. That is, the economy can be in two different regimes, and accordingly the parameters in (2) are no longer constant, but depend on the state. In order to emphasise the dependence on the reinsurance price, we let the safety coefficient of the reinsurer depend on the current regime by letting all other variables unchanged. Thus, instead of (2), we now consider the process

$$\begin{array}{cc}\hfill X_t^B& =x+{\int }_0^t{\theta }_{J_s}{b}_s-\lambda \mu ({\theta }_{J_s}-\eta )\mathrm{d}s+{\int }_0^t\sqrt{\lambda {\mu }_2}{b}_s\mathrm{d}{W}_s,\hfill \end{array}$$

The set of admissible reinsurance strategies will be denoted by $\mathcal{B}$ and is formally defined as

$$\mathcal{B}=\{B=\left\{b_t\right\},b_t\in [0,1],b_t{\mathcal{F}}_t-\mathrm{adapted}\}.$$

We are interested in the minimal value of expected discounted capital injections by starting in state i with initial capital x over all admissible strategies, i.e., we minimise

$$V^B(i,x):={\mathbb{E}}_{i,x}\left[{\int }_0^{\infty }{e}^{-\delta t}\mathrm{d}{Y}_t^B\right],(i,x)\in \{1,2\}\times [0,\infty ).$$

Here, and in the following, we use the common notation $\mathbb{E}[.|X_0=x,J_0=i]={\mathbb{E}}_{i,x}[.]$.

Our target is to find an admissible reinsurance strategy $B^{\ast }$ such that the value function

$$\begin{array}{c}\hfill V(i,x):=\underset{B\in \mathcal{B}}{inf}{V}^B(i,x).\end{array}$$

can be written as the return function corresponding to the strategy $B^{\ast }$, i.e., $V(i,x)=V^{B^{\ast }}(i,x)$.

According to the theory of stochastic control, we expect the value function V is to solve the HJB equation (see Schmidli (2008) or Jiang and Pistorius (2012) for a model with Markov-switching)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \underset{b\in [0,1]}{inf}\frac{\lambda {\mu }_2{b}^2}{2}{V}^{″}(i,x)+\lambda \mu \left({\theta }_{i}b-{\theta }_i+\eta \right)V^{\prime }(i,x)-(\delta -q_{ii})V(i,x)-q_{ii}V(j,x)=0,\end{array}\end{array}$$

(4)

The boundary condition $V^{\prime }(i,0)=-1$ arises from the requirement of smooth fit (${\mathcal{C}}^1$-fit) at zero. As we do not allow the surplus to become negative, it is clear that the value function for $x<0$ fulfils $V(i,x)=-x+V(i,0)$, i.e., we immediately inject as much capital as it is needed to shift the surplus process to zero, meaning $V^{\prime }(i,x)=-1$. The second boundary condition $\underset{x\to \infty }{lim}V(i,x)=0$ originates from the fact that a Brownian motion with a positive drift converges to infinity almost surely, i.e., for $x\to \infty$ the amount of expected discounted capital injections converges to 0, see, for instance, Rolski et al. (1999).

The HJB equation can be formally derived as the infinitesimal version of the dynamic programming principle, upon assuming that V has the regularity needed to apply Ito’s formula for Markov-modulated diffusion processes (as in the proof of Lemma 1 below); we refer to Chapter 2 in Schmidli (2008) for a textbook treatment.

It is clear that $b^{\ast }(i,x)=-\frac{\mu {\theta }_i{V}^{\prime }(i,x)}{{\mu }_2{V}^{″}(i,x)}\wedge 1$. If $b^{\ast }<1$, the HJB becomes for $i,j\in \{1,2\}$$i\ne j$

$$\begin{array}{c}\hfill -\frac{{\mu }^2{\theta }_i^2{V}^{\prime }{(i,x)}^2}{2{\mu }_2{V}^{″}(i,x)}-\lambda \mu \left({\theta }_i-\eta \right)V^{\prime }(i,x)-(\delta -q_{ii})V(i,x)-q_{ii}V(j,x)=0.\end{array}$$

(5)

Technically, HJB Equation (5) is a system of 2 ordinary differential equations, coupled through the transition rates of the underlying Markov chain. It is a hard task to explicitly solve these equations and show that the solutions are decreasing and convex functions of the initial capital. Therefore, we use a recursive method to obtain the value function as a limit. However, first we look at the constant strategies and investigate why none of those can be optimal in the case of more than one regime.

3. Constant Strategies

It is known (see, for instance, Eisenberg and Schmidli (2009)) that in the one-regime case the optimal reinsurance strategy is given by a constant. In this section, we show that in the two-regime case a constant strategy (other than “no reinsurance at all”) cannot be optimal.

Let $b_1,b_2\in [0,1)$, then $\widehat{B}:=\left\{b_{J_t}\right\}$ is an admissible reinsurance strategy. Further, we let $\widehat{Y}$ and $\widehat{X}$ denote the capital injection process triggered by the reinsurance strategy $\widehat{B}$ and the surplus process under $\widehat{B}$ and after capital injections, respectively. Thus, for the return function $\widehat{V}(i,x)$, corresponding to $\widehat{B}$, it holds

$$\widehat{V}(i,x)={\mathbb{E}}_{i,x}\left[{\int }_0^{\infty }{e}^{-\delta t}\mathrm{d}{\widehat{Y}}_t\right].$$

Lemma 1.

If $\widehat{u}$ is a solution to the system of ODEs, $i\in \{1,2\}$, $j\ne i$

$$\begin{array}{c}\hfill \frac{\lambda {\mu }_2{b}_i^2}{2}{\widehat{u}}^{″}(i,x)+\lambda \mu ({\theta }_i{b}_i-{\theta }_i+\eta ){\widehat{u}}^{\prime }(i,x)-(\delta -q_{ii})\widehat{u}(i,x)-q_{ii}\widehat{u}(j,x)=0,\end{array}$$

(6)

with boundary conditions ${\widehat{u}}^{\prime }(i,0)=-1$, $\underset{x\to \infty }{lim}\widehat{u}(i,x)=0$, then $\widehat{u}(i,x)=\widehat{V}(i,x)$.

Proof. 

First, we look at the Equation (6). A general solution to (6) fulfilling $\underset{x\to \infty }{lim}\widehat{u}(i,x)=0$ is given by

$$\begin{array}{c}\hfill \widehat{u}(i,x)=C_{i1}{e}^{A_{1}x}+C_{i2}{e}^{A_{2}x},i\in \{1,2\},\end{array}$$

(7)

where $A_1,A_2<0$. The coefficients are uniquely given by

$$\begin{array}{cc}\hfill & C_{11}\left(A_1^2\frac{\lambda {\mu }_2{b}_1^2}{2}+A_1\lambda \mu ({\theta }_1{b}_1-{\theta }_1+\eta )-(\delta -q_{11})\right)-q_{11}{C}_{21}=0,\hfill \\ \hfill & C_{12}\left(A_2^2\frac{\lambda {\mu }_2{b}_1^2}{2}+A_2\lambda \mu ({\theta }_1{b}_1-{\theta }_1+\eta )-(\delta -q_{11})\right)-q_{11}{C}_{22}=0,\hfill \\ \hfill & C_{21}\left(A_1^2\frac{\lambda {\mu }_2{b}_2^2}{2}+A_1\lambda \mu ({\theta }_2{b}_2-{\theta }_2+\eta )-(\delta -q_{22})\right)-q_{22}{C}_{11}=0,\hfill \\ \hfill & C_{22}\left(A_2^2\frac{\lambda {\mu }_2{b}_2^2}{2}+A_2\lambda \mu ({\theta }_2{b}_2-{\theta }_2+\eta )-(\delta -q_{22})\right)-q_{22}{C}_{12}=0,\hfill \\ \hfill & C_{11}{A}_1+C_{12}{A}_2=-1,C_{21}{A}_1+C_{22}{A}_2=-1.\hfill \end{array}$$

(8)

Now, arguing like in Shreve et al. (1984), we show that $\widehat{u}=\widehat{V}$. Using a generalised form of Itô’s formula, like it has been done, for instance, in Jiang and Pistorius (2012), we get

$$\begin{array}{cc}\hfill & e^{-\delta t}\widehat{u}(J_t,{\widehat{X}}_t)=\widehat{u}(J_0,{\widehat{X}}_0)+{\int }_0^t{e}^{-\delta s}{\widehat{u}}^{\prime }(J_s,{\widehat{X}}_s)\mathrm{d}{W}_s+M_t\hfill \\ & +{\int }_0^t{e}^{-\delta s}\{\frac{\lambda {\mu }_2{b}_s^2}{2}{\widehat{u}}^{″}(J_s,{\widehat{X}}_s)+\lambda \mu ({\theta }_{J_s}{b}_s-{\theta }_{J_s}+\eta ){\widehat{u}}^{\prime }(J_s,{\widehat{X}}_s)\hfill \\ & -(\delta -q_{J_s,J_s})\widehat{u}(J_s,{\widehat{X}}_s)-q_{J_s,J_s}\widehat{u}\left(11_{[J_s=1]}+1,{\widehat{X}}_s\right)\}\mathrm{d}s\hfill \\ & +{\int }_0^t{e}^{-\delta s}{\widehat{u}}^{\prime }(J_s,{\widehat{X}}_s)\mathrm{d}{\widehat{Y}}_s\hfill \end{array}$$

(9)

where M is a local martingale associated with the regime switching mechanism. That is, M is given by

$$\begin{array}{c}\hfill M_t={\int }_{[0,t]\times [0,2]}\widehat{u}({\widehat{X}}_{s-},j)-\widehat{u}({\widehat{X}}_{s-},J_{s-})\tilde{\pi }(\mathrm{d}s,\mathrm{d}j),\end{array}$$

where $\tilde{\pi }=\pi -\nu$ is a compensated random measure as defined in Jacod and Shiryaev (2003). It holds

$$\begin{array}{cc}\hfill & \pi (\mathrm{d}t,\mathrm{d}j)=\sum _{s\ge 0}11_{[\Delta J_s\left(\omega \right)\ne 0]}11_{(\mathrm{d}t,\mathrm{d}j)}(s,J_s\left(\omega \right)),\hfill \\ & \nu (\mathrm{d}t,\mathrm{d}j)=11_{[\Delta J_t\left(\omega \right)\ne 0]}{q}_{J_{t-},j}\mathcal{P}\left(\mathrm{d}j\right)\mathrm{d}t,\hfill \end{array}$$

where $\mathcal{P}$ is the counting measure on $\{1,2\}$ and $11_{(\mathrm{d}t,\mathrm{d}j)}(s,J_s\left(\omega \right))$ is the Dirac measure at the point $(s,J_s\left(\omega \right))$.

Note that M is bounded because

$$|\widehat{u}(j,{\widehat{X}}_{s-})-\widehat{u}(J_{s-},{\widehat{X}}_{s-})|\le \underset{i\in \{1,2\}}{max}\widehat{u}(i,{\widehat{X}}_{s-}),$$

As $\widehat{u}$ is bounded, we can conclude that M is a martingale with expectation zero. Because ${\widehat{u}}^{\prime }$ is bounded we can conclude that also the stochastic integral is a martingale with expectation zero. Further, as $\widehat{u}$ solves Equation (6) with $u^{\prime }(i,0)=-1$, building expectations on the both sides in (9) yields

$$\begin{array}{c}\hfill \mathbb{E}\left[e^{-\delta t}\widehat{u}(J_t,{\widehat{X}}_t)\right]=\widehat{u}(i,x)-\mathbb{E}\left[{\int }_0^t{e}^{-\delta s}\mathrm{d}{\widehat{Y}}_s\right].\end{array}$$

By the bounded convergence theorem, we can interchange limit $t\to \infty$ and expectations and get $\widehat{u}(i,x)=\mathbb{E}\left[{\int }_0^{\infty }{e}^{-\delta s}\mathrm{d}{\widehat{Y}}_s\right]$.    □

Remark 1.

Note that it holds $C_{ij}\ne 0$ for all $i,j\in \{1,2\}$. If, for example, $C_{11}=0$ then we have from (8) that it must also hold $C_{21}=0$ and simultaneously

$$\begin{array}{cc}\hfill & C_{22}=-\frac{1}{A_2}=C_{12},\hfill \\ & C_{22}\left(A_2^2\frac{\lambda {\mu }_2{b}_1^2}{2}+A_2\lambda \mu ({\theta }_1{b}_1-{\theta }_1+\eta )-(\delta -q_{11})-q_{11}\right)=0,\hfill \\ & C_{22}\left(A_2^2\frac{\lambda {\mu }_2{b}_2^2}{2}+A_2\lambda \mu ({\theta }_2{b}_2-{\theta }_2+\eta )-(\delta -q_{22})-q_{22}\right)=0,\hfill \end{array}$$

which leads to a contradiction.

Lemma 2.

The return function corresponding to a strategy $B=\left\{b_{J_s}\right\}$ with $b_i\in [0,1)$, $i=1,2$, does not solve HJB Equation (4).

Proof. 

We proof this lemma by contradiction. Assume that there is a strategy $B=\left\{b_{J_s}\right\}$ with $b_i\in [0,1)$, $i=1,2$, so that the corresponding return function, $\widehat{V}$ solves the HJB Equation (4), i.e.,

$$\begin{array}{c}\hfill 0=\underset{b\in [0,1]}{inf}\left[\frac{\lambda {\mu }_2{b}^2}{2}{\widehat{V}}^{″}(i,x)+\lambda \mu ({\theta }_{i}b-{\theta }_i+\eta ){\widehat{V}}^{\prime }(i,x)-(\delta -q_{ii})\widehat{V}(i,x)-q_{ii}\widehat{V}(j,x)\right].\end{array}$$

Subtracting Equation (6) from the latter, one finds

$$\begin{array}{cc}\hfill 0& =\underset{b\in [0,1]}{inf}\left[\frac{\lambda {\mu }_2{b}^2}{2}{\widehat{V}}^{″}(i,x)+\lambda \mu ({\theta }_{i}b-{\theta }_i+\eta ){\widehat{V}}^{\prime }(i,x)-(\delta -q_{ii})\widehat{V}(i,x)-q_{ii}\widehat{V}(j,x)\right.\hfill \\ \hfill & -\left.\left(\frac{\lambda {\mu }_2{\widehat{b}}_i^2}{2}{\widehat{V}}^{″}(i,x)+\lambda \mu ({\theta }_i{\widehat{b}}_i-{\theta }_i+\eta ){\widehat{V}}^{\prime }(i,x)-(\delta -q_{ii})\widehat{V}(i,x)-q_{ii}\widehat{V}(j,x)\right)\right],\hfill \end{array}$$

which is equivalent to

$$\begin{array}{c}\hfill \frac{\lambda {\mu }_2{\widehat{b}}^2}{2}{\widehat{V}}^{″}(i,x)+\lambda \mu {\theta }_i{\widehat{b}}_i{\widehat{V}}^{\prime }(i,x)=\underset{b\in [0,1]}{inf}\left[\frac{\lambda {\mu }_2{b}^2}{2}{\widehat{V}}^{″}(i,x)+\lambda \mu {\theta }_i{b}_i{\widehat{V}}^{\prime }(i,x)\right].\end{array}$$

Then, it should hold that $b_i=-\frac{\mu {\theta }_i}{{\mu }_2}\frac{{\widehat{V}}^{\prime }(i,x)}{{\widehat{V}}^{″}(i,x)}\wedge 1$. As we assumed ${\widehat{b}}_i<1$, for $i=1,2$, the expression $-\frac{\mu {\theta }_i}{{\mu }_2}\frac{{\widehat{V}}^{\prime }(i,x)}{{\widehat{V}}^{″}(i,x)}$ should not depend on x. Keeping in mind the boundary conditions $\underset{x\to \infty }{lim}\widehat{V}(i,x)=0$ and $V^{\prime }(i,0)=-1$, we get

$$\begin{array}{c}\hfill {\widehat{V}}^{\prime }(i,x)=-e^{-\frac{\mu {\theta }_i}{{\mu }_2{\widehat{b}}_i}x},\end{array}$$

which contradicts (7).    □

4. Recursion

In the following, we establish an algorithm that allows us to calculate the value function as a limit of a sequence of twice continuously differentiable, decreasing and convex functions. For simplicity, we let

$$\begin{array}{cc}\hfill & {\Delta }_i:=\frac{\lambda {\mu }^2{\theta }_i^2}{2{\mu }_2}+\delta -q_i,B_i:=\frac{{\Delta }_i}{\lambda \mu ({\theta }_i-\eta )},{\tilde{B}}_i:=\frac{{\Delta }_i+q_i}{\lambda \mu ({\theta }_i-\eta )}.\hfill \end{array}$$

(10)

We will see that the behaviour of the optimal reinsurance strategy will depend on the relations between $B_1$, ${\tilde{B}}_1$, $B_2$ and ${\tilde{B}}_2$. As there are many possibilities to arrange the above quantities, we consider just one path, omitting considering the case of no reinsurance, in order not to complicate the explanations. However, the algorithm proposed below can be applied to any combination of parameters.

Assumption 1.

W.l.o.g. we assume that $B_1>{\tilde{B}}_1>B_2>{\tilde{B}}_2>max\{\frac{\mu {\theta }_2}{{\mu }_2},\frac{\mu {\theta }_1}{{\mu }_2}\}$, which is equivalent to

$$\frac{1}{B_1}<\frac{1}{{\tilde{B}}_1}<\frac{1}{B_2}<\frac{1}{{\tilde{B}}_2}<max\left\{\frac{{\mu }_2}{\mu {\theta }_2},\frac{{\mu }_2}{\mu {\theta }_1}\right\}.$$

(11)

In the case of just one regime, the problem could be solved by conjecturing that the optimal strategy is constant and the corresponding return function is an exponential function. This allowed us to verify easily that the solution, say v, to Differential Equation (5) with $q_{ii}=0$ was strictly increasing, convex and fulfilled $-\frac{\mu {\theta }_i{v}^{\prime }}{{\mu }_2{v}^{″}}<1$ or, if the case maybe, the optimal strategy was not to buy reinsurance at all. In our case of two regimes, the situation changes as we have seen in Section 3. The return functions corresponding to the constant strategies do not solve the HJB equation in general.

As it is impossible to guess the optimal strategy and to subsequently check whether the return function corresponding to this strategy is the value function, we slightly change the solving procedure. At first, we look at the HJB equation in the form of Differential Equation (5), such has been done, for instance, in Eisenberg and Schmidli (2009). The next, very technical step is to solve the obtained differential equation and to check whether the solution, say f, fulfils $-\frac{\mu {\theta }_i{f}^{\prime }(i,x)}{{\mu }_2{f}^{″}(i,x)}\in (0,1)$. Then, we show that the gained solution f is indeed the return function corresponding to the reinsurance strategy $-\frac{\mu {\theta }_i{f}^{\prime }(i,x)}{{\mu }_2{f}^{″}(i,x)}\in (0,1)$. Thus, in this way we find an admissible strategy whose return function solves HJB Equation (4). Verification theorem proves this return function to be the value function.

In the following, we describe the steps of an algorithm allowing to get a strictly decreasing and convex solution to Differential Equation (5) under Assumption (11). The procedure consists in choosing a starting function, fixing, say $i=1$, and replacing the unknown function $V(2,x)$ in (5) by the chosen starting function. Using the method of Højgaard and Taksar (1998), we show the existence and uniqueness of a solution. In the next step, now it holds $i=2$, the unknown function $V(1,x)$ in (5) is replaced by the function obtained in the first step. Letting the number of steps go to infinity, we get a solution to (5). We will see that the starting value of the recursion plays a crucial role in obtaining a solution with the desired properties: convexity and monotonicity. Therefore, we start by explaining how to choose the starting function in Step 0.

4.1. Step 0

The solutions to the differential equations

$$-\frac{\lambda {\mu }^2{\theta }_i^2{f}^{\prime }{\left(x\right)}^2}{{\mu }_2{f}^{″}\left(x\right)}-\lambda \mu \left({\theta }_i-\eta \right)f^{\prime }\left(x\right)-\delta f\left(x\right)=0.$$

(12)

with boundary conditions $\underset{x\to \infty }{lim}f\left(x\right)=0$ and $f^{\prime }\left(0\right)=-1$ are well known and given by $\frac{1}{{\tilde{B}}_i}{e}^{-{\tilde{B}}_{i}x}$. Note that due to Assumption (11) it holds that $\frac{\mu {\theta }_i}{{\mu }_2{\tilde{B}}_i}<1$.

The optimal strategy in the case of two regimes is not constant, see Section 3. However, we conjecture that the value function in the case of two regimes fulfils $\underset{x\to \infty }{lim}\frac{-V^{\prime }(1,x)}{V^{″}(1,x)}=\underset{x\to \infty }{lim}\frac{-V^{\prime }(2,x)}{V^{″}(2,x)}\in [1/{\tilde{B}}_1,1/{\tilde{B}}_2]$, i.e., the ratio of the first and second derivatives converges to the same value does not matter the initial regime state. One can see it as a sort of averaging of the optimal strategies from the one-regime cases. This means, for instance, that if in one-regime case the optimal reinsurance level was low in the first state and high in the second, in the two-regime case the optimal level in the first state will go up and go down in the second.

Mathematically, the above explanations are reflected in the starting function of our algorithm:

$$W_0\left(x\right):=\frac{1}{\Lambda }{e}^{-\Lambda x},$$

(13)

where $\Lambda$ fulfils

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \lambda \mu ({\theta }_1-\eta )\left\{\frac{B_1}{\Lambda }-1\right\}+\frac{q_{11}}{\Lambda }{e}^{\alpha }=0,\hfill \\ & \lambda \mu ({\theta }_2-\eta )\left\{\frac{B_2}{\Lambda }-1\right\}+\frac{q_{22}}{\Lambda }{e}^{-\alpha }=0.\hfill \end{array}\end{array}$$

(14)

It means that $\Lambda$ and $\alpha$ are uniquely given by

$$\begin{array}{cc}\hfill & \Lambda =\frac{B_1+B_2-\sqrt{D}}{2}\mathrm{and}\alpha =ln\left(\frac{B_1-B_2+\sqrt{D}}{-2q_{11}/\lambda \mu ({\theta }_1-\eta )}\right),\hfill \\ & D:={(B_1-B_2)}^2+\frac{4q_{11}{q}_{22}}{{\lambda }^2{\mu }^2({\theta }_1-\eta )({\theta }_2-\eta )}.\hfill \end{array}$$

(15)

Remark 2.

It is a straightforward calculation to show $\Lambda \in ({\tilde{B}}_2,{\tilde{B}}_1)$ and $\alpha >0$ using the definition of $B_i$ and Assumption (11).

Note that due to our assumption, it holds that $-\frac{\mu {\theta }_i{W}_0^{\prime }\left(x\right)}{{\mu }_2{W}_0^{″}\left(x\right)}<1$.

We will see by establishing the algorithm below that Equation (14) for $\Lambda$ will be crucial in order to obtain a well-defined solution to (5). The elaborated mathematical meaning and explanation of the Equation (14) will be demonstrated in the following steps.

4.2. Step 1

Assuming (11), we start investigating Differential Equation (5) and substitute the term $-q_{11}V(2,x)$ by $-q_{11}{W}_0\left(x\right)$ defined in (13), i.e., we look now at the differential equation

$$-\frac{\lambda {\mu }^2{\theta }_1^2{f}^{\prime }{\left(x\right)}^2}{{\mu }_2{f}^{″}\left(x\right)}-\lambda \mu \left({\theta }_1-\eta \right)f^{\prime }\left(x\right)-(\delta -q_{11})f\left(x\right)-q_{11}{W}_0\left(x\right)=0.$$

(16)

Although we know the function $W_0$, differential Equation (16) still cannot be solved in a way that we could easily prove the solution to be strictly decreasing and convex. Therefore, we use the following technique, introduced in Højgaard and Taksar (1998).

We assume that there is a strictly increasing, bijective on ${\mathbb{R}}_{+}$ function g such that the derivative of the solution f to (16) fulfils $f^{\prime }\left(g\left(x\right)\right)=-e^{-x}$. Then, it holds

$$f^{″}\left(g\left(x\right)\right)=\frac{e^{-x}}{g^{\prime }\left(x\right)}\mathrm{and}{f}^{‴}\left(g\left(x\right)\right)=-\frac{e^{-x}}{g^{\prime }{\left(x\right)}^2}-\frac{e^{-x}}{g^{\prime }{\left(x\right)}^3}.$$

Differentiating (16) yields

$$-\frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}\left(2f^{\prime }\left(x\right)-\frac{f^{\prime }{\left(x\right)}^2{f}^{‴}\left(x\right)}{{\left(f^{″}\left(x\right)\right)}^2}\right)-\left({\theta }_1-\eta \right)f^{″}\left(x\right)-(\delta -q_{11})f^{\prime }\left(x\right)-q_{11}{W}_0^{\prime }\left(x\right)=0.$$

Changing the variable to $g\left(x\right)$ leads to a new differential equation for the function g

$$\begin{array}{c}\hfill \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}\left(2e^{-x}+\frac{-\frac{e^{-3x}}{{\left(g^{\prime }\left(x\right)\right)}^2}\left(1+\frac{g^{″}\left(x\right)}{g^{\prime }\left(x\right)}\right)}{e^{-2x}/{\left(g^{\prime }\left(x\right)\right)}^2}\right)-\lambda \mu \left({\theta }_1-\eta \right)\frac{e^{-x}}{g^{\prime }\left(x\right)}+(\delta -q_{11})e^{-x}\\ \hfill -q_{11}{W}_0^{\prime }\left(g\left(x\right)\right)=0,\end{array}$$

which can be further simplified by multiplying by $e^{-x}{g}^{\prime }\left(x\right)$ and inserting the definition of $B_1$:

$$\begin{array}{c}\hfill \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{g}^{″}\left(x\right)=\lambda \mu ({\theta }_1-\eta )\left\{g^{\prime }\left(x\right)B_1-1\right\}+q_{11}{g}^{\prime }\left(x\right)e^{-\Lambda g\left(x\right)+x}.\end{array}$$

(17)

As the function g should be bijective, we will prove the existence and uniqueness of a solution to (17) on ${\mathbb{R}}_{+}$ with the boundary conditions guaranteeing $g\left({\mathbb{R}}_{+}\right)={\mathbb{R}}_{+}$ and $g^{\prime }>0$. In particular, the term $e^{-\Lambda g\left(x\right)+x}$ determines the unique condition yielding $g^{\prime }>0$ and $\underset{x\to \infty }{lim}{g}^{\prime }\left(x\right)\ne \infty$, namely, $\underset{x\to \infty }{lim}{g}^{\prime }\left(x\right)=1/\Lambda$.

In order to guide the reader through the auxiliary results below, we provide a roadmap identifying the key findings of Step 1.

Note that when investigating (17) we are looking at a boundary value problem. To show the existence and uniqueness of a solution, we will translate the boundary value problem into an initial value problem, i.e., we shift the condition $g^{\prime }\left(x\right)\to 1/\Lambda$ as $x\to \infty$ to $x=0$ by using Volterra type representation for (17).

• First, we show that if (17) has a solution, say g, with the boundary values $g\left(0\right)=0$ and $g^{\prime }\left(n\right)=1/\Lambda$, for some $n\in \mathbb{N}$, then $g^{\prime }\left(0\right)\in \left(1/{\tilde{B}}_1,1/\Lambda \right)$.
• In the second step, we show that there is a unique solution ${\xi }_n$ to (17) with the boundary conditions ${\xi }_n\left(0\right)=0$ and ${\xi }_n^{\prime }\left(n\right)=1/\Lambda$.
• We prove the existence of a solution $g_1$ to (17) with $g\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_1\left(x\right)=1/\Lambda$.
• It holds $g_1^{\prime }\left(0\right)\in (1/{\tilde{B}}_1,1/\Lambda )$ and $g_1^{″}\left(x\right)>0$.
• The inverse function $h_1$ of $g_1$ fulfils $h_1^{\prime }\in (\Lambda ,{\tilde{B}}_1)$ and $h_1^{″}\left(x\right)<0$.
• $h_1\left(x\right)\in (\Lambda x,{\tilde{B}}_{1}x)$ for all $x>0$ and $\underset{x\to \infty }{lim}\left(h_1\left(x\right)-\Lambda x\right)=\alpha$ with $\alpha$ given in (15).

Lemma 3.

If there is a solution g to (17) with the boundary conditions $g\left(0\right)=0$, $g^{\prime }\left(n\right)=1/\Lambda$, for some $n\ge 1$, then it holds that

$$g^{\prime }\left(0\right)\in \left(1/{\tilde{B}}_1,1/\Lambda \right).$$

Proof. 

We prove the claim by contradiction. Let g be a solution to (17) with the boundary conditions $g\left(0\right)=0$, $g^{\prime }\left(n\right)=1/\Lambda$.

• Assume for the moment that $g^{\prime }\left(0\right)\le 1/{\tilde{B}}_1$. Then,

  $$\begin{array}{c}\hfill \left\{\begin{array}{cc}\frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{g}^{″}\left(0\right)=\lambda \mu \left({\theta }_1-\eta \right)(g^{\prime }\left(0\right){\tilde{B}}_1-1)<0\hfill & :\mathrm{if}{g}^{\prime }\left(0\right)<1/{\tilde{B}}_1,\hfill \\ \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{g}^{‴}\left(0\right)=q_{11}{g}^{\prime }\left(0\right)(-g^{\prime }\left(0\right)\Lambda +1)<0\hfill & :\mathrm{if}{g}^{\prime }\left(0\right)=1/{\tilde{B}}_1,\hfill \end{array}\right.\end{array}$$

  meaning that $g^{\prime }\left(x\right)$ stays positive but below $1/{\tilde{B}}_1$ in an environment of 0. As ${\tilde{B}}_1>\Lambda$, the function $e^{-\Lambda g\left(x\right)+x}$ is increasing in an $\epsilon$-environment of 0, i.e., $e^{-\Lambda g\left(x\right)+x}>1$, which means that

  $$\begin{array}{cc}\hfill \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{g}^{″}\left(x\right)& =(g^{\prime }\left(x\right)B_1-1)\lambda \mu \left({\theta }_1-\eta \right)+q_{11}{g}^{\prime }\left(x\right)e^{-\Lambda g\left(x\right)+x}\hfill \\ & <(g^{\prime }\left(x\right){\tilde{B}}_1-1)\lambda \mu \left({\theta }_1-\eta \right)<0.\hfill \end{array}$$
  Thus, $g^{″}$ will stay negative and $g^{\prime }$ will never arrive at $1/\Lambda >1/{\tilde{B}}_1$.
• On the other hand, if $g^{\prime }\left(0\right)\ge 1/\Lambda >1/{\tilde{B}}_1$, then in a similar way one concludes that $g^{\prime }$ stays above $1/\Lambda$ for all $x\in (0,n]$, contradicting $g^{\prime }\left(n\right)=1/\Lambda$.
• Thus, we can conclude that $g^{\prime }\left(0\right)\in \left(1/{\tilde{B}}_1,1/\Lambda \right)$.

   □

Lemma 4.

For every $n\ge 1$ there is a unique solution ${\xi }_n\left(x\right)$ to (17) on $[0,n]$ fulfilling ${\xi }_n\left(0\right)=0$, ${\xi }_n^{\prime }\left(n\right)=1/\Lambda$.

Proof. 

As the proof is very technical, we postpone it to Appendix A.    □

Lemma 5.

Let ${\xi }_n$ be the unique solution to (17) with the boundary conditions ${\xi }_n\left(0\right)=0$ and ${\xi }_n^{\prime }\left(n\right)=1/\Lambda$. Then, ${\xi }_n^{″}\left(x\right)>0$ on $[0,\infty )$.

Proof. 

See Appendix A.    □

Proposition 1.

There exists a unique solution $g_1$ to (17) with the boundary conditions $g_1\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_1^{\prime }\left(x\right)=1/\Lambda$, $g_1^{\prime }\in (1/{\tilde{B}}_1,1/\Lambda )$ and $g_1^{″}>0$ on $(0,\infty )$.

Proof. 

See Appendix A.    □

Remark 3.

Proposition 1 implies $-\Lambda g_1\left(x\right)+x>0$ for all $x>0$. Moreover, due to Equation (14), it holds that

$$\underset{x\to \infty }{lim}(-\Lambda g_1\left(x\right)+x)=\alpha .$$

(18)

where α is given in (15).

Note that the definition of g yields $g^{\prime }\left(x\right)=-\frac{f^{\prime }\left(g\left(x\right)\right)}{f^{″}\left(g\left(x\right)\right)}$. The boundary conditions imply $\underset{x\to \infty }{lim}g\left(x\right)=\infty$ and $\underset{x\to \infty }{lim}{g}^{″}\left(x\right)=0$. Thus, letting $x\to \infty$ in (17) and using (18), we get the first equation in (14). This provides the first idea and the meaning of the choice of $W_0$.

Corollary 1.

There is a strictly increasing and concave inverse function of $g_1$ on ${\mathbb{R}}_{+}$: $g_1^{-1}\left(x\right)=:h_1\left(x\right)$. Further, it holds

• $h_1$ fulfils $h_1^{\prime }\left(x\right)>0$, $h_1^{\prime }\in (\Lambda ,{\tilde{B}}_1)$, $\underset{x\to \infty }{lim}{h}_1^{\prime }\left(x\right)=\Lambda$ and $h_1^{″}\left(x\right)<0$.
• $h_1^{\prime }\left(x\right)-\Lambda =\frac{1}{{\tilde{g}}_1^{\prime }\left(h_1\left(x\right)\right)}-\Lambda >0$, i.e., $h_1\left(x\right)-\Lambda x$ is strictly increasing with $h_1\left(x\right)>\Lambda x$ for $x>0$.

Proof. 

The function $g_1$ fulfils $g_1\left({\mathbb{R}}_{+}\right)={\mathbb{R}}_{+}$ and $g_1^{\prime }\left(x\right)>0$ for all $x\in {\mathbb{R}}_{+}$, i.e., $g_1$ is a bijective function, which implies the existence of an inverse function $h_1$. All other properties follow from the properties of $g_1$.    □

We can now let

$$W_1\left(x\right)={\int }_x^{\infty }{e}^{-h_1\left(y\right)}\mathrm{d}y,$$

(19)

i.e., $W_1^{\prime }\left(x\right)=-e^{-h_1\left(x\right)}$. Note that $W_1$ is well defined due to Corollary 1 and solves Differential Equation (17) with the boundary conditions $W_1^{\prime }\left(0\right)=-1$ and $\underset{x\to \infty }{lim}{W}_1\left(x\right)=0$. In particular, due to (11) it holds $-\frac{\mu {\theta }_1{W}_1^{\prime }}{{\mu }_2{W}_1^{″}}<1$.

In the following second step, we construct in a similar way a function $g_2$.

4.3. Step 2

In the second step, we add the term $-q_{22}{g}^{\prime }\left(x\right)W_1^{\prime }\left(g\left(x\right)\right)e^x$ to Differential Equation (12), i.e., we are looking at the differential equation

$$\begin{array}{cc}\hfill \frac{\lambda {\mu }^2{\theta }_2^2}{2{\mu }_2}{g}^{″}\left(x\right)& =\lambda \mu \left({\theta }_2-\eta \right)(g^{\prime }\left(x\right)B_2-1)-q_{22}{g}^{\prime }\left(x\right)W_1^{\prime }\left(g\left(x\right)\right)e^x\hfill \\ & =\lambda \mu \left({\theta }_2-\eta \right)(g^{\prime }\left(x\right)B_2-1)+q_{22}{g}^{\prime }\left(x\right)e^{-h_1\left(g\left(x\right)\right)+x}.\hfill \end{array}$$

(20)

Note that $h_1\left(x\right)\in C^{\infty }$, which implies Lipschitz-continuity on compacts. The existence of a solution $g_2$ with the boundary conditions $g_2\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_2\left(x\right)=1/\Lambda$ can be shown similar to Step 1.

The main findings of Step 2 are as follows:

• There is a unique solution $g_2$ to (20) with the boundary conditions $g_2\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_2\left(x\right)=1/\Lambda$.
• It holds that $g_2^{\prime }\left(0\right)\in (1/\Lambda ,1/{\tilde{B}}_2)$ and $g_2^{″}\left(x\right)<0$.
• The inverse function $h_2$ of $g_2$ fulfils $h_2^{\prime }\in ({\tilde{B}}_2,\Lambda )$ and $h_2^{″}\left(x\right)>0$.
• $h_2\left(x\right)\in ({\tilde{B}}_{2}x,\Lambda x)$ for all $x>0$.

In the following, we prove only the results that cannot be easily transferred from Step 1.

Lemma 6.

If there is a solution $g_2$ to Differential Equation (20) with the boundary conditions $g_2\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_2^{\prime }\left(x\right)=1/\Lambda$ then $g_2^{\prime }\left(0\right)\in (1/\Lambda ,1/{\tilde{B}}_2)$.

Proof. 

See Appendix A.    □

Lemma 7.

Let $g_2$ be the unique solution to Differential Equation (20) with the boundary conditions $g_2\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_2^{\prime }\left(x\right)=1/\Lambda$. Then, $g_2^{″}\left(x\right)<0$ for all $x\in (0,\infty )$.

Proof. 

Lemma 6 yields $g_2^{\prime }\left(0\right)\in (1/\Lambda ,1/{\tilde{B}}_2)$. It means (see (A1)) that $g_2^{″}\left(0\right)<0$. Let $\widehat{x}:=inf\{x>0:g_2^{″}\left(x\right)=0\}$, then $g_2^{\prime }\left(\widehat{x}\right)\in (1/\Lambda ,1/{\tilde{B}}_2)$ because if $g_2^{\prime }\left(x\right)=1/\Lambda$ and $g_2^{″}\left(x\right)<0$ then Lemma 6 gives $\underset{x\to \infty }{lim}{g}_2^{\prime }\left(x\right)\ne 1/\Lambda$. Further, we also have

$$g_2^{‴}\left(\widehat{x}\right)=q_{22}{g}_2^{\prime }\left(\widehat{x}\right)\left(-g_2^{\prime }\left(\widehat{x}\right)h_1^{\prime }\left(g_2\left(\widehat{x}\right)\right)+1\right)e^{-h_1\left(g_2\left(\widehat{x}\right)\right)+\widehat{x}}>0$$

because $h_1^{\prime }>\Lambda$ due to Corollary 1. Then, $g_2^{″}$ becomes positive, i.e., $g_2^{\prime }$ becomes increasing and stays increasing for $g_2^{\prime }>1/\Lambda$, i.e., bounded away from $1/\Lambda$, which yields a contradiction.    □

Corollary 2.

Let $h_2\left(x\right)$ be the inverse function of $g_2\left(x\right)$. Then, $h_2^{\prime }\in ({\tilde{B}}_2,\Lambda )$, $h_2^{″}>0$, $\underset{x\to \infty }{lim}{h}_2^{\prime }\left(x\right)=\Lambda$, $h_2\left(x\right)\in ({\tilde{B}}_{2}x,\Lambda x)$.

Proof. 

The proof is a direct consequence of Lemma 7.    □

Remark 4.

• Let

  $$\beta :=\underset{x\to \infty }{lim}(-h_1\left(g_2\left(x\right)\right)+x).$$

  (21)
  Then, due to Equation (14) it holds $\beta =-\alpha$.
• Furthermore, it follows easily

  $$\begin{array}{cc}\hfill \beta =\underset{x\to \infty }{lim}(-h_1\left(g_2\left(x\right)\right)+x)& =\underset{x\to \infty }{lim}\left(-h_1\left(g_2\left(h_2\left(x\right)\right)\right)+h_2\left(x\right)\right)\hfill \\ & =\underset{x\to \infty }{lim}\left(-h_1\left(x\right)+\Lambda x-\Lambda x+h_2\left(x\right)\right),\hfill \end{array}$$
  and using (18) we get $\underset{x\to \infty }{lim}(-h_1\left(x\right)+\Lambda x)=\underset{x\to \infty }{lim}(-x+\Lambda g_1\left(x\right))=-\alpha$.
  Therefore, we conclude

  $$\underset{x\to \infty }{lim}\left(-\Lambda x+h_2\left(x\right)\right)=\beta +\alpha =0$$

  (22)
  as $h_2\left(x\right)\le \Lambda x$, see Corollary 2

Remark 4 explains the second equation in (14), if we let $x\to \infty$ in Differential Equation (20).

4.4. Step 2m+1

In this step, we are searching for the function $h_{2m+1}$ as the inverse of the solution $g_{2m+1}$ to the differential equation ${\mathcal{U}}_m\left(g\right)=0$, where

$$\begin{array}{c}\hfill {\mathcal{U}}_m\left(g\right):=\frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{g}^{″}\left(x\right)-\lambda \mu \left({\theta }_1-\eta \right)(g^{\prime }\left(x\right)B_1-1)-q_{11}{g}^{\prime }\left(x\right)e^{-h_{2m}\left(g\left(x\right)\right)+x}.\end{array}$$

(23)

The existence of a solution $g_{2m+1}$ can be proven similarly to Step 1. The boundary conditions are $g_{2m+1}\left(0\right)=0$ and $\underset{x\to \infty }{lim}{g}_{2m+1}^{\prime }\left(x\right)=1/\Lambda$.

Our main target is to show that the obtained sequences of functions $\left(g_{2m+1}\right)$, $\left(g_{2m+1}^{\prime }\right)$ and $\left(h_{2m+1}\right)$ are monotone. We carry out the proof by induction using as the induction step $h_2\left(x\right)<\Lambda x$ on $(0,\infty )$, see Corollary 2 in Step 2.

The main findings of Step $2m+1$ are summarised in the following remark.

Remark 5.

Similar to Step 1, we get for $g_{2m+1}$ and its inverse function $h_{2m+1}$ that

• $g_{2m+1}\left(0\right)=h_{2m+1}\left(0\right)=0$;
• $g_{2m+1}^{\prime }\in (1/{\tilde{B}}_1,1/\Lambda )$, $h_{2m+1}^{\prime }\in (\Lambda ,{\tilde{B}}_1)$;
• $g_{2m+1}^{″}\left(x\right)>0$ and $h_{2m+1}^{″}\left(x\right)<0$; and
• $\underset{x\to \infty }{lim}{g}_{2m+1}^{\prime }\left(x\right)=1/\Lambda$, $\underset{x\to \infty }{lim}{h}_{2m+1}^{\prime }\left(x\right)=\Lambda$.

In Lemma 8 we show

• $g_{2m+1}^{\prime }>g_{2m-1}^{\prime }$ on ${\mathbb{R}}_{+}$, $g_{2m+1}>g_{2m-1}$ and $h_{2m+1}<h_{2m-1}$ on $(0,\infty )$.

Lemma 8.

Assume that the functions $h_{2k}$ obtained in Steps $2k$, $0\le k\le m$, fulfil

1.

$h_{2k}\left({\mathbb{R}}_{+}\right)={\mathbb{R}}_{+}$, $h_{2k}\left(0\right)=0$, $h_{2k}^{\prime }\in ({\tilde{B}}_2,\Lambda )$, $h_{2k}\left(x\right)\le \Lambda x$, $\underset{x\to \infty }{lim}{h}_{2k}=\Lambda$ and $h_{2k}^{″}\left(x\right)>0$, and

2.

$h_{2k}\left(x\right)<h_{2k-2}\left(x\right)$ for $x>0$.

Then, $g_{2m+1}^{\prime }>g_{2m-1}^{\prime }$ on ${\mathbb{R}}_{+}$, $g_{2m+1}>g_{2m-1}$ and $h_{2m+1}<h_{2m-1}$ on $(0,\infty )$.

Proof. 

See Appendix A.    □

Remark 6.

Similar to Step 2, (22) we conclude by induction from $\underset{x\to \infty }{lim}\left(-h_{2m}\left(g_{2m+1}\left(x\right)\right)+x\right)=\alpha$ that

$$\underset{x\to \infty }{lim}\left(-\Lambda x+h_{2m+1}\left(x\right)\right)=m\beta +(m+1)\alpha =\alpha$$

(24)

with α given in (15).

4.5. Step $2m+2$

In Step $2m+2$, we are searching for the function $h_{2m+2}$ as the inverse of the solution $g_{2m+2}$ to the differential equation ${\mathcal{G}}_m\left(g_{2m+2}\right)=0$, where

$$\begin{array}{c}\hfill {\mathcal{G}}_m\left(g\right)=\frac{\lambda {\mu }^2{\theta }_2^2}{2{\mu }_2}{g}^{″}\left(x\right)-\lambda \mu \left({\theta }_2-\eta \right)(g^{\prime }\left(x\right)\Lambda -1)-q_{22}{g}^{\prime }\left(x\right)e^{-h_{2m+1}\left(g\left(x\right)\right)+x}.\end{array}$$

The main findings of this Step are summarised below.

Remark 7.

Similar to Step $2m+1$, we get for $g_{2m+2}$ and its inverse function $h_{2m+2}$:

• $g_{2m+2}\left(0\right)=h_{2m+2}\left(0\right)=0$;
• $g_{2m+2}^{\prime }\in (1/\Lambda ,1/{\tilde{B}}_2)$, $h_{2m+2}^{\prime }\in ({\tilde{B}}_2,\Lambda )$;
• $g_{2m+2}^{″}\left(x\right)<0$ and $h_{2m+2}^{″}\left(x\right)>0$;
• $\underset{x\to \infty }{lim}{g}_{2m+2}^{\prime }\left(x\right)=1/\Lambda$, $\underset{x\to \infty }{lim}{h}_{2m+2}^{\prime }\left(x\right)=\Lambda$.
• $g_{2m+1}^{\prime }>g_{2m-1}^{\prime }$ on ${\mathbb{R}}_{+}$, $g_{2m+1}>g_{2m-1}$ and $h_{2m+1}<h_{2m-1}$ on $(0,\infty )$.

Remark 8.

Similar to Step $2m+1$, (24) we conclude by induction from $\underset{x\to \infty }{lim}(-h_{2m+1}\left(g_{2m+2}\left(x\right)\right)+x)=\beta =-\alpha$ that

$$\underset{x\to \infty }{lim}\left(-\Lambda x+h_{2m+2}\left(x\right)\right)=(m+1)\beta +(m+1)\alpha =0.$$

(25)

Note that the choice of $\Lambda$ and $\alpha$ given by (14) is the only choice leading to $\beta =-\alpha$. In this way, one makes sure that it holds $\underset{x\to \infty }{lim}\left(-\Lambda x+h_{2m+2}\left(x\right)\right)=0$ also in the $2m+2$-th step, implying in this way $0\le h_{2m+2}\left(x\right)\le \Lambda x$. A different choice of $\Lambda$ and $\alpha$ would eliminate the upper boundary for $h_{2m+2}$, destroying the well-definiteness of the limiting function $\underset{m\to \infty }{lim}{h}_{2m+2}$.

5. The Value Function

In this section, we first construct a candidate for the value function by letting $m\to \infty$ for the sequences $\left(g_{2m+1}\right)$, $\left(g_{2m+2}\right)$, $\left(g_{2m+1}^{\prime }\right)$, $\left(g_{2m+2}^{\prime }\right)$, $\left(h_{2m+1}\right)$ and $\left(h_{2m+2}\right)$. Then, we prove the candidate to be the value function via a verification theorem.

We know from Remarks 5 and 7 that the sequences $\left(g_{2m+1}\right)$, $\left(g_{2m+2}\right)$, $\left(g_{2m+1}^{\prime }\right)$, $\left(g_{2m+2}^{\prime }\right)$, $\left(h_{2m+1}\right)$ and $\left(h_{2m+2}\right)$ are monotone and thus pointwise convergent. In the following lemma, we show that the convergence is uniform on compacts.

Lemma 9.

The sequences $\left(g_{2m+1}\right)$, $\left(g_{2m}\right)$, $\left(g_{2m+1}^{\prime }\right)$, $\left(g_{2m}^{\prime }\right)$, $\left(g_{2m+1}^{″}\right)$, $\left(g_{2m}^{″}\right)$, $h_{2m}$, $h_{2m+1}$ converge uniformly on compacts to ${\mathfrak{g}}_{\mathfrak{1}}$, ${\mathfrak{g}}_{\mathfrak{2}}$, ${\left({\mathfrak{g}}_{\mathfrak{1}}\right)}^{\prime }$, ${\left({\mathfrak{g}}_{\mathfrak{2}}\right)}^{\prime }$, ${\left({\mathfrak{g}}_{\mathfrak{1}}\right)}^{″}$, ${\left({\mathfrak{g}}_{\mathfrak{2}}\right)}^{″}$, ${\mathfrak{h}}_{\mathfrak{2}}={\mathfrak{g}}_{\mathfrak{2}}^{-1}$, ${\mathfrak{h}}_{\mathfrak{1}}={\mathfrak{g}}_{\mathfrak{1}}^{-1}$, respectively.

Proof. 

Note first that Lemma 8 and Remark 7 yield the monotonicity of the sequences $\left(g_{2m+1}\right)$, $\left(g_{2m}\right)$, $\left(g_{2m+1}^{\prime }\right)$, $\left(g_{2m}^{\prime }\right)$, $h_{2m}$ and $h_{2m+1}$. Therefore, these sequences converge pointwise implying the pointwise convergence of $\left(g_{2m+1}^{″}\right)$ and $\left(g_{2m}^{″}\right)$.

In the following, we show that $\left(g_{2m+1}\right)$, $\left(g_{2m+1}^{\prime }\right)$, $\left(g_{2m+1}^{″}\right)$ and $h_{2m+1}$ converge uniformly on compacts.

Assume $\left(g_{2m+1}\right)$, $\left(g_{2m+1}^{\prime }\right)$, $\left(g_{2m+1}^{″}\right)$ and $h_{2m+1}$ converge pointwise to ${\mathfrak{g}}_{\mathfrak{1}}$, w, u and $\chi$, respectively. Note that it holds by definition of $g_{2m+1}$: ${\mathcal{U}}_m\left(g_{2m+1}\right)=0$ with ${\mathcal{U}}_m$ defined in (23). As $g_{2m+1}^{″}>0$ and $g_{2m+1}^{\prime }>0$, see Step $2m+1$, it holds

$$\begin{array}{c}\hfill \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{g}_{2m+1}^{″}\left(x\right)<\lambda \mu \left({\theta }_1-\eta \right)(g_{2m+1}^{″}\left(x\right)B_1-1).\end{array}$$

Integrating both sides of the above inequality, yields

$$\begin{array}{cc}\hfill \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}\left(g_{2m+1}^{\prime }\left(x\right)-g_{2m+1}^{\prime }\left(0\right)\right)& ={\int }_0^x{g}_{2m+1}^{″}\left(y\right)\mathrm{d}y\hfill \\ & <\frac{2{\mu }_2\lambda \mu \left({\theta }_1-\eta \right)}{\lambda {\mu }^2{\theta }_1^2}{\int }_0^x(g_{2m+1}^{\prime }\left(y\right)B_1-1)\mathrm{d}y\hfill \\ & =\frac{2{\mu }_2\lambda \mu \left({\theta }_1-\eta \right)}{\lambda {\mu }^2{\theta }_1^2}\left(g_{2m+1}\left(x\right)B_1-x\right),\hfill \end{array}$$

which means that the sequence $\left(g_{2m+1}^{″}\right)$ is dominated by an integrable function. By Lebesgue’s convergence theorem ${\int }_0^x{g}_{2m+1}^{″}\left(y\right)\mathrm{d}y$ converges pointwise to ${\int }_0^{x}u\left(y\right)\mathrm{d}y$. Recall that ${\int }_0^{x}u\left(y\right)\mathrm{d}y$ is a continuous function of x and because of the uniqueness of the pointwise limit $\left(g_{2m+1}^{\prime }\right)$ converges pointwise to $w={\int }_0^{x}u\left(y\right)\mathrm{d}y$. That is, as $\left(g_{2m+1}^{\prime }\right)$ is a decreasing sequence, Dini’s theorem yields the uniform convergence of $\left(g_{2m+1}^{\prime }\right)$ to w on compacts.

With the same arguments we get that $\left(g_{2m+1}\right)$ converges uniformly to ${\mathfrak{g}}_{\mathfrak{1}}$ and it holds $w={\mathfrak{g}}_{\mathfrak{1}}^{\prime }$ on compacts.

In a similar way, one can conclude that $\left(g_{2m}\right)$ and $\left(g_{2m}^{\prime }\right)$ converge uniformly on compacts to ${\mathfrak{g}}_{\mathfrak{2}}$ and ${\mathfrak{g}}_{\mathfrak{2}}^{\prime }$, respectively. So that we can conclude, compare, for instance, (de Souza and Silva 2001, pp. 60, 297), that the sequence of the inverse functions $\left(h_{2m}\right)$ converges uniformly on compacts to ${\mathfrak{h}}_{\mathfrak{2}}$, the inverse of ${\mathfrak{g}}_{\mathfrak{2}}$.

As a consequence of Differential Equation ${\mathcal{U}}_m\left(g_{2m+1}\right)=0$, also the sequence $\left(g_{2m+1}^{″}\right)$ converges uniformly to ${\mathfrak{g}}_{\mathfrak{1}}^{″}$ on compacts.    □

Lemma 10.

The limiting functions ${\mathfrak{g}}_{\mathfrak{1}}$, ${\mathfrak{g}}_{\mathfrak{2}}$, ${\mathfrak{h}}_{\mathfrak{1}}$ and ${\mathfrak{h}}_{\mathfrak{2}}$ fulfil on $(0,\infty )$

• ${\mathfrak{g}}_i\left(0\right)={\mathfrak{h}}_i\left(0\right)=0$;
• ${\mathfrak{g}}_{\mathfrak{1}}^{\prime }\in (1/{\tilde{B}}_1,1/\Lambda )$, ${\mathfrak{h}}_{\mathfrak{1}}^{\prime }\in (\Lambda ,{\tilde{B}}_1)$;
• ${\mathfrak{g}}_{\mathfrak{2}}^{\prime }\in (1/\Lambda ,1/{\tilde{B}}_2)$, ${\mathfrak{h}}_{\mathfrak{2}}^{\prime }\in ({\tilde{B}}_2,\Lambda )$;
• ${\mathfrak{g}}_{\mathfrak{1}}^{″}\left(x\right)>0$ and ${\mathfrak{g}}_{\mathfrak{2}}^{″}<0$;
• ${\mathfrak{h}}_{\mathfrak{1}}^{″}\left(x\right)>0$ and ${\mathfrak{h}}_{\mathfrak{2}}^{″}\left(x\right)<0$.
• ${\mathfrak{g}}_{\mathfrak{1}}$, ${\mathfrak{g}}_{\mathfrak{2}}$, ${\mathfrak{h}}_{\mathfrak{1}}$ and ${\mathfrak{h}}_{\mathfrak{2}}$ fulfil on $(0,\infty )$

  $$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{\lambda {\mu }^2{\theta }_1^2}{2{\mu }_2}{\mathfrak{g}}_{\mathfrak{1}}^{″}\left(x\right)=\lambda \mu \left({\theta }_1-\eta \right)({\mathfrak{g}}_{\mathfrak{1}}^{\prime }\left(x\right)B_1-1)+q_{11}{\mathfrak{g}}_{\mathfrak{1}}^{\prime }\left(x\right)e^{-{\mathfrak{h}}_{\mathfrak{2}}\left({\mathfrak{g}}_{\mathfrak{1}}\left(x\right)\right)+x},\hfill \\ & \frac{\lambda {\mu }^2{\theta }_2^2}{2{\mu }_2}{\mathfrak{g}}_{\mathfrak{2}}^{″}\left(x\right)=\lambda \mu \left({\theta }_2-\eta \right)({\mathfrak{g}}_{\mathfrak{2}}^{\prime }\left(x\right)B_2-1)+q_{22}{\mathfrak{g}}_{\mathfrak{2}}^{\prime }\left(x\right)e^{-{\mathfrak{h}}_{\mathfrak{1}}\left({\mathfrak{g}}_{\mathfrak{2}}\left(x\right)\right)+x}.\hfill \end{array}\end{array}$$

  (26)

Proof. 

From Lemma 9, one gets immediately the above inequalities with ≥ and ≤ instead of > and < along with Differential Equation (26).

The strict inequalities follow easily using the methods presented in Step 1.    □

Lemma 11.

For $i\in \{1,2\}$ it holds $\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{i}}^{\prime }\left(x\right)=1/\Lambda$ and $\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{i}}^{″}\left(x\right)=0$.

Proof. 

Note that ${\mathfrak{g}}_{\mathfrak{1}}^{\prime }$ is increasing and ${\mathfrak{g}}_{\mathfrak{2}}^{\prime }$ is decreasing with ${\mathfrak{g}}_{\mathfrak{1}}^{\prime }\le 1/\Lambda$, ${\mathfrak{g}}_{\mathfrak{2}}^{\prime }\ge 1/\Lambda$, ${\mathfrak{h}}_{\mathfrak{1}}^{\prime }\ge \Lambda$ and ${\mathfrak{h}}_{\mathfrak{2}}^{\prime }\le \Lambda$. It means

$$-{\mathfrak{h}}_{\mathfrak{2}}^{\prime }\left({\mathfrak{g}}_{\mathfrak{1}}\left(x\right)\right){\mathfrak{g}}_{\mathfrak{1}}^{\prime }\left(x\right)+1\ge 0and-{\mathfrak{h}}_{\mathfrak{1}}^{\prime }\left({\mathfrak{g}}_{\mathfrak{2}}\left(x\right)\right){\mathfrak{g}}_{\mathfrak{2}}^{\prime }\left(x\right)+1\le 0,$$

meaning that $e^{-{\mathfrak{h}}_{\mathfrak{2}}\left({\mathfrak{g}}_{\mathfrak{1}}\left(x\right)\right)+x}$ is increasing and $e^{-{\mathfrak{h}}_{\mathfrak{1}}\left({\mathfrak{g}}_{\mathfrak{2}}\left(x\right)\right)+x}$ is decreasing. If $\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{1}}^{\prime }\left(x\right)<1/\Lambda$, then $\underset{x\to \infty }{lim}{e}^{-{\mathfrak{h}}_{\mathfrak{2}}\left({\mathfrak{g}}_{\mathfrak{1}}\left(x\right)\right)+x}=\infty$ contradicting ${\mathfrak{g}}_{\mathfrak{1}}^{″}\ge 0$. If $\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{2}}^{\prime }\left(x\right)>1/\Lambda$, then $\underset{x\to \infty }{lim}{e}^{-{\mathfrak{h}}_{\mathfrak{2}}\left({\mathfrak{g}}_{\mathfrak{2}}\left(x\right)\right)+x}=0$ leads to the contradiction $\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{2}}^{″}\left(x\right)>0$.

$\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{i}}^{″}\left(x\right)=0$ is a direct consequence from the above.    □

Corollary 3.

It holds ${\int }_x^{\infty }{e}^{-{\mathfrak{h}}_{\mathfrak{i}}\left(y\right)}\mathrm{d}y<\infty$, $i\in \{1,2\}$.

Proof. 

As ${\mathfrak{h}}_{\mathfrak{1}}\left(0\right)=0={\mathfrak{h}}_{\mathfrak{2}}\left(0\right)$ and ${\mathfrak{h}}_{\mathfrak{1}}^{\prime }\in (\Lambda ,{\tilde{B}}_1)$ and ${\mathfrak{h}}_{\mathfrak{2}}^{\prime }\in ({\tilde{B}}_2,\Lambda )$ we conclude ${\mathfrak{h}}_{\mathfrak{1}}\in (\Lambda x,{\tilde{B}}_{1}x)$ and ${\mathfrak{h}}_{\mathfrak{2}}\in ({\tilde{B}}_{2}x,\Lambda x)$. Therefore,

$$\begin{array}{c}\hfill {\int }_x^{\infty }{e}^{-{\mathfrak{h}}_{\mathfrak{i}}\left(y\right)}\mathrm{d}y<{\int }_x^{\infty }{e}^{-{\tilde{B}}_{2}y}\mathrm{d}y<\infty .\end{array}$$

   □

Definition 1.

We let for $i\in \{1,2\}$

$$\tilde{V}(i,x):={\int }_x^{\infty }{e}^{-{\mathfrak{h}}_{\mathfrak{i}}\left(y\right)}\mathrm{d}y$$

(27)

and

$$b^{\ast }(i,x):=\frac{\mu {\theta }_i}{{\mu }_2{\mathfrak{h}}_{\mathfrak{i}}^{\prime }\left(x\right)}.$$

(28)

Lemma 12.

The function $\tilde{V}$ and $b^{\ast }$ defined in (27) and (28), respectively, fulfil

• ${\tilde{V}}^{\prime }(i,x)=-e^{-{\mathfrak{h}}_{\mathfrak{i}}\left(x\right)}$, $b^{\ast }(i,x)\in (0,1)$.
• $\tilde{V}(i,x)$ solves the system of differential equations for $i,j\in \{1,2\}$ with $i\ne j$

  $$\begin{array}{c}\hfill \frac{\lambda {\mu }_2{b}^{\ast }{(i,x)}^2}{2}{\tilde{V}}^{″}(i,x)+\lambda \mu ({\theta }_i{b}^{\ast }(i,x)-{\theta }_i+\eta ){\tilde{V}}^{\prime }(i,x)-(\delta -q_{ii})\tilde{V}(i,x)\\ \hfill -q_{ii}\tilde{V}(j,x)=0\end{array}$$

  with the boundary conditions ${\tilde{V}}^{\prime }(i,0)=-1$ and $\underset{x\to \infty }{lim}\tilde{V}(i,x)=0$.
• $\tilde{V}(i,x)$ is the return function corresponding to the strategy $b^{\ast }(i,x)$.

Proof. 

• It follows directly from (27), (28), (11), ${\mathfrak{h}}_{\mathfrak{2}}\in ({\tilde{B}}_2,\Lambda )$ and ${\mathfrak{h}}_{\mathfrak{1}}\in (\Lambda ,{\tilde{B}}_1)$.
• The functions ${\mathfrak{g}}_{\mathfrak{1}}$, ${\mathfrak{g}}_{\mathfrak{2}}$, ${\mathfrak{h}}_{\mathfrak{1}}$ and ${\mathfrak{h}}_{\mathfrak{2}}$ solve the system of equations (26) with boundary conditions ${\mathfrak{g}}_{\mathfrak{1}}\left(0\right)={\mathfrak{g}}_{\mathfrak{2}}\left(0\right)=0$, $\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{1}}^{\prime }\left(x\right)=\underset{x\to \infty }{lim}{\mathfrak{g}}_{\mathfrak{1}}^{\prime }\left(x\right)=1/\Lambda$ and ${\mathfrak{g}}_{\mathfrak{1}}\left({\mathfrak{h}}_{\mathfrak{1}}\left(x\right)\right)=x$, ${\mathfrak{g}}_{\mathfrak{2}}\left({\mathfrak{h}}_{\mathfrak{2}}\left(x\right)\right)=x$.
  It holds ${\tilde{V}}^{\prime }(1,{\mathfrak{g}}_{\mathfrak{i}}\left(x\right))=-e^{-x}$, i.e.,

  $$\begin{array}{cc}\hfill & {\tilde{V}}^{″}(i,{\mathfrak{g}}_{\mathfrak{i}}\left(x\right))=\frac{e^{-x}}{{\mathfrak{g}}_{\mathfrak{i}}^{\prime }\left(x\right)}and{\tilde{V}}^{‴}(i,{\mathfrak{g}}_{\mathfrak{i}}\left(x\right))=-\frac{e^{-x}}{{\mathfrak{g}}_{\mathfrak{i}}{\left(x\right)}^2}-\frac{e^{-x}{\mathfrak{g}}_{\mathfrak{i}}^{″}\left(x\right)}{{\mathfrak{g}}_{\mathfrak{i}}^{\prime }{\left(x\right)}^3}.\hfill \end{array}$$
  Substituting ${\mathfrak{g}}_{\mathfrak{i}}$ by $\tilde{V}(i,x)$ in (26) yields the desired result.
• Similar to Shreve et al. (1984) and Section 3, one gets that

  $$\tilde{V}(i,x)=\mathbb{E}\left[{\int }_0^{\infty }{e}^{-\delta t}\mathrm{d}{Y}_t^{b^{\ast }}\right],$$

  where $\left\{Y_t^{B^{\ast }}\right\}$ describes the capital injection process corresponding to the strategy $B^{\ast }=\left\{b^{\ast }\right\}$ with $b^{\ast }$ defined in (28).

   □

Proposition 2.

The function $\tilde{V}(i,x)$ defined in (27) is strictly decreasing, convex, fulfils $V^{\prime }(i,0)=-1$, $\underset{x\to \infty }{lim}V(i,x)=0$ and solves HJB Equation (4).

Proof. 

The proof follows easily from Lemma 12.    □

Theorem 1 (Verification Theorem).

The strategy $B^{\ast }=\left\{b^{\ast }\right\}$ with

$$b^{\ast }(i,x)=\frac{\mu {\theta }_i}{{\mu }_2{\mathfrak{h}}_{\mathfrak{i}}^{\prime }\left(x\right)}<1$$

is the optimal strategy, and the corresponding return function $\tilde{V}$, given in (27), is the value function.

Proof. 

Let $B=\left\{b_t\right\}$ be an arbitrary admissible strategy and $X^B$ the surplus process under B and after the capital injections. Following the steps from lemma 1 we get

$$\begin{array}{cc}\hfill & e^{-\delta \left(t\right)}\tilde{V}(J_t,X_t^B)=\tilde{V}(J_0,X_0^B)+{\int }_0^t{e}^{-\delta s}{\tilde{V}}^{\prime }(J_s,X_s^B)\mathrm{d}{W}_s+M_t\hfill \\ & +{\int }_0^t{e}^{-\delta s}\{\frac{\lambda {\mu }_2{b}_s^2}{2}{\tilde{V}}^{″}(J_s,X_s^B)+\lambda \mu ({\theta }_{J_s}{b}_s-{\theta }_{J_s}+\eta ){\tilde{V}}^{\prime }(J_s,X_s^B)\hfill \\ & -(\delta -q_{J_s,J_s})\tilde{V}(J_s,X_s^B)-q_{J_s,J_s}\tilde{V}\left(11_{[J_s=1]}+1,X_s^B\right)\}\mathrm{d}s\hfill \\ & +{\int }_0^t{e}^{-\delta s}{\tilde{V}}^{\prime }(J_s,X_s^B)\mathrm{d}{Y}_s^B\hfill \end{array}$$

(29)

where M is again a martingale with expectation 0 as M is bounded

$$|\tilde{V}(X_{s-}^B,j)-\tilde{V}(X_{s-}^B,J_{s-})|\le \underset{i\in \{1,2\}}{max}{V}^1(X_{s-}^B,i),$$

where $V^1$ is the return function corresponding to the strategy “no reinsurance”, i.e., $b\equiv 1$. Because ${\tilde{V}}^{\prime }$ is bounded we can conclude that also the stochastic integral is a martingale with expectation zero. Further, as $\tilde{V}$ solves the HJB equation and is convex with $V^{\prime }(i,0)=-1$, it follows that

$$\begin{array}{c}\hfill \frac{\lambda {\mu }_2{b}_s^2}{2}{\tilde{V}}^{″}(J_s,X_s^B)+\lambda \mu ({\theta }_{J_s}{b}_s-{\theta }_{J_s}+\eta ){\tilde{V}}^{\prime }(J_s,X_s^B)-(\delta -q_{J_s,J_s})\tilde{V}(J_s,X_s^B)\\ \hfill -q_{J_s,J_s}\tilde{V}\left(11_{[J_s=1]}+1,X_s^B\right)\ge 0\end{array}$$

and ${\tilde{V}}^{\prime }(i,x)\ge -1$. Thus, building expectations on the both sides in (29) yields

$$\begin{array}{c}\hfill \mathbb{E}\left[e^{-\delta t}\tilde{V}(J_t,X_t^B)\right]\ge \tilde{V}(i,x)-\mathbb{E}\left[{\int }_0^t{e}^{-\delta s}\mathrm{d}{Y}_s^B\right].\end{array}$$

By the bounded convergence theorem, we can interchange limit $t\to \infty$ and expectations and get $\tilde{V}(i,x)\le \mathbb{E}\left[{\int }_0^{\infty }{e}^{-\delta s}\mathrm{d}{Y}_s^B\right]$. For the strategy $B^{\ast }=\left\{b^{\ast }(J_s,X_s)\right\}$, we get the equality.    □

6. Numerical Illustrations

All numerical computations were performed on Matlab R2020b using the library bvpsuite2.01.

The package bvpsuite2.0 has been developed at the Institute for Analysis and Scientific Computing, Vienna University of Technology, and can be used—among other applications—for the numerical solution of boundary value problems in ordinary differential equations on semi-infinite intervals. The library uses collocation for the numerical solution of the underlying boundary value problems, which is a piecewise polynomial function which satisfies the given ODE at a finite number of nodes (collocation points). This approach shows advantageous convergence properties compared to other direct higher order methods.2

The subsequent example has an illustrative scope and aims solely at providing numerical evidence of the convergence of the recursive algorithm developed in the previous theoretical sections. In order to provide a clear numerical illustration of the results, the following choice of parameters turns out to be suitable: $\mu =1$, $\lambda =1$, ${\mu }_2=4$, $\eta =0.3$, ${\theta }_1=0.33$, ${\theta }_2=0.8$, $q_{11}=-0.6$ and $q_{22}=-0.4$.

6.1. Illustration of the Recursive Procedure

We start with illustrating the recursive procedure described in Section 4; that is, we consider the functions $h_{2k}^{\prime }$ and $h_{2k+1}^{\prime }$ and their convergence behaviour.

The fast convergence of each $h_i^{\prime }\left(x\right)\to \Lambda$ for $x\to \infty$ results in a very badly conditioned differential equation (in this example, $\Lambda \sim 0.295522$). As a consequence, we had to truncate the solution interval and set the boundary conditions at $x=500$, i.e., $h_i^{\prime }\left(500\right)=\Lambda$. The short horizon leads the solver to overshoot the solution at the beginning and to compensate later on, creating a characteristic initial “hump” for $h_{2k+1}^{\prime }$, which is evened out more and more with each iteration. As a matter of fact, the larger k is, the more $h_{2k+1}^{\prime }$ is converging towards a concave shape. On the other hand, the convergence of $h_{2k}^{\prime }$ is faster, as for 6 iterations the functions are already very close to each other, see Figure 1.

6.2. Solving the HJB Directly

Differently than in the section above, the HJB Equation (4) is now solved directly using again bvpsuite2.0. The optimal reinsurance strategy $b^{\ast }(i,x)$ and the ratio $-\frac{V^{″}(i,x)}{V^{\prime }(i,x)}$ corresponding to the limit of the sequences $\left(h_{2m+1}^{\prime }\right)$ (for $i=1$) and $\left(h_{2m}^{\prime }\right)$ (for $i=2$) are illustrated in Figure 2.

In Figure 2, left picture, one sees that in both regimes the optimal reinsurance strategy is non-constant with respect to the surplus level. The red line representing the optimal strategy $b^{\ast }(1,x)$ is decreasing and the strategy $b^{\ast }(2,x)$, in blue, is increasing. The reason for this behaviour is the relation ${\tilde{B}}_1>\Lambda >{\tilde{B}}_2$, where ${\tilde{B}}_1=1.787083$ and ${\tilde{B}}_2=0.24$. The optimal strategies for the one-regime cases are

$$\frac{\mu {\theta }_1}{{\mu }_2{\tilde{B}}_1}=0.046165<\frac{\mu {\theta }_1}{{\mu }_2\Lambda }=0.350947$$

for the regime 1 and

$$\frac{\mu {\theta }_2}{{\mu }_2{\tilde{B}}_2}=0.83>\frac{\mu {\theta }_1}{{\mu }_2\Lambda }=0.676769$$

for regime 2, respectively. Thus, due to the possibility to change into a regime with a different reinsurance price, the optimal strategy changes.

In Figure 2, right picture, we see that the ratio $-\frac{V^{″}(i,x)}{V^{\prime }(i,x)}$ converges for both regimes to the level $0.295522=\Lambda$, the value towards which the sequences $h_{2k}^{\prime }$ and $h_{2k+1}^{\prime }$ converge as can be seen in Section 6.1.

Using bvpsuite optimisation, each step of the iteration takes between 30 and 40 s to compute on a single 3GHz core.

These findings confirm the validity of the iterative algorithm developed in Section 4 and illustrated numerically in Section 6.1.

7. Conclusions

In this paper, we study the problem faced by an insurance company that aims at finding the optimal proportional reinsurance strategy minimising the expected discounted capital injections. We assume that the cost of entering the proportional reinsurance contract depends on the current busyness cycle of a two-state economy, and we model this by letting the safety loading of the reinsurance be modulated by a continuous-time Markov chain. This leads to an optimal reinsurance problem under regime switching. In order to simplify our explanations, we assume a certain relation between the crucial parameters of the two regimes. Considering all possible combinations would be space-consuming with just a marginal additional value.

Differently to Eisenberg and Schmidli (2009)—where no regime switching has been considered—we find that the optimal reinsurance cannot be independent of the current value of the surplus process, but should instead be given as a feedback strategy $b^{\ast }$, also depending on the current regime. However, due to the complex nature of the resulting HJB equation, determining an explicit expression for $b^{\ast }$ turns out to be a challenging task. For this reason, we develop a recursive algorithm that hinges on the construction of two sequences of functions converging uniformly to a classical solution to the HJB equation and simultaneously providing the optimal strategies for both regimes. The obtained optimal strategies are monotone with respect to the surplus level and converge for both regimes to the same explicitly calculated constant as the surplus goes to infinity. The algorithm is illustrated by a numerical example, where one can also see that the convergence to the solution of the HJB equation is quite fast.

The recursive scheme represents the main contribution of this paper and might also be applied (with necessary adjustments) to other optimisation problems containing regime-switching. For this reason, we retrace here the main steps and ideas of the algorithm.

The differential equation for the value function is first translated into a differential equation for an auxiliary function, transforming the derivative of the value function into an exponential function, using the method of Højgaard and Taksar (1998). In order to get a solution to the system of equations, say Equations (a) and (b), we solve the differential Equation (a) assuming that the solution to Equation (b) is given by an exponential function $e^{-\Lambda }$, which is the starting function of our algorithm. Then, we solve Equation (b) by inserting the solution to Equation (a) from the previous step. Proceeding in this manner, we obtain two sequences of uniformly converging functions whose limiting functions solve the original HJB equation system.

Note that we are facing a boundary value problem, i.e., the boundary conditions on the value function and its derivative, $\underset{x\to \infty }{lim}V(i,x)=0$ and $V^{\prime }(i,0)=-1$ are given at different boundaries, with one boundary being infinity. Therefore, the usual Picard–Lindelöf approach does not work. Instead, we use Volterra form representations and comparison theorems to show the existence and uniqueness of a solution with the desired properties.

One of the crucial points in the above considerations is the starting function of the algorithm. It turns out that there is a uniquely given constant, $\Lambda$, allowing to obtain the desired properties of the limiting functions.

In particular, we show that the derivatives of the auxiliary functions lie in suitable intervals $({\tilde{B}}_2,\Lambda )$ or $(\Lambda ,{\tilde{B}}_1)$, depending on the differential equation we are looking at.

Thus, we are able to show that the value function is twice continuously differentiable and the optimal strategy has a monotone character and converges for $x\to \infty$ to an explicitly calculated value.

It would be interesting to implement the considerations from Chiappori et al. (2006) to extend the presented model by hidden information, for instance, by introducing a hidden Markov chain governing the reinsurance price over the parameter $\theta$. This topic will be one of the directions of our future research.