Transmission dynamics and stability of fractional order derivative model for COVID-19 epidemic with optimal control analysis

1.   Introduction

The insurance company is a common financial institution in our real life. Its profit mainly comes from two aspects: premium income and investment income, and the risks it needs to face mainly include: compensation risk and investment risk. In the past decade, more and more scholars have focused on building appropriate risk models $($r.m.$)$ to describe the various situations that insurance companies may face^[1,2,3]. At first, the r.m. studied by the researchers was a classical risk model that only considered a company's claims as a negative jump. For example, Zhang et al.^[4] studied a new method to estimate the Gerber-Shiu discount penalty function $($p.f.$)$ under the classical r.m., and Peng et al.^[5] studied a r.m. of dividend payment with perturbations. But in reality, an insurance company's random returns should also be taken into account. To better fit the actual situation, Boucheire et al.^[6] first proposed the two-sided jumps r.m., which is used to extend the r.m. of a single jump. Here, it is considered that the company's revenue is random, which is also a random variable, then the random revenue is a non-negative jump, and a negative jump is a claim. Since then, this model has been paid much attention by many researchers. For example, E.C.K. Cheung^[7] studied a renewal model with continuous expenses and bidirectional jumps, where the amplitude of the jumps and the time intervals of arrival time are random. From this, E.C.K. Cheung obtained the updated equation of the discounted penalty funtion (e.d.p.f.) with defects. Zhang^[8] considered the problem of e.d.p.f. for a two-sided jumps r.m. with dividend payout and obtained some explicit expressions. Wang et al. ^[9] considered the investment r.m. under the bilateral jump and tried to obtain the maximum surplus through the appropriate investment proportion. Xu et al.^[10] studied the problem of ruin probability under bilateral jumps with random observations. For more research on two-side jump r.m., we can refer to references ^{[11,12,13,14,15]}.

Subsequently, some scholars put forward the dividend barrier strategy, that is, they set a threshold value $b > 0$, and pay dividends to shareholders when the company's earnings are greater than $b$. The strategy was first proposed by De Finetti. Then, Gerber et al.^[16] studied the threshold dividend strategy, and Yin et al.^[17] and Cossette et al.^[18] put forward the horizontal barrier strategy. The multi-tier dividend strategy can be learned from Xie and Zou^[19]. To make the r.m. more realistic, some scholars have added dividend barriers to the study of bilateral jump risk models. For example, Bo et al. ^[20] studied the Lévy model with bilateral jumps under the dividend barrier strategy and Chen et al.^[21] studied the dividend payment and the reward and e.d.p.f. of the dividend strategy with a threshold under the compound Poisson (c.p.) model. The integral differential equation (IDE) is derived under the boundary conditions and the approximate solution (a.s.) is approximated by the sinc numerical method. When studying the c.p. model with proportional investment, Chen and Ou^[21] added the dividend with threshold value. Inspired by the above research, we propose a bilateral jumps model with a threshold strategy under random observation.

We introduce our work in the following parts. In the second section, we construct the two-sided jumps risk model with investment interest under random observation, and the observational intervals obey a same exponential distribution. In the third section, we obtain the IDE of the expected discounted dividend payment (e.d.d.p.) function. To solve this equation, in the fourth section, we introduce an excellent numerical method to the solution of the IDEs and get the upper boundary of the error between the a.s. and the real solution. This numerical method is called the sinc numerical method. In the last section, we give some numerical examples to explore the effects of the included parameters on the e.d.d.p..

2.   The model

According to the previous research on the bilateral jump r.m., we define

where $u_0$ represents the company's initial surplus on the account and $u_0$ is greater than zero. In addition, $\{U(t)\}_{t\ge 0}$ stands for the surplus process, while $c$ represents the premium rate paid by the insured, so obviously $c > 0$. Here the two stochastic processes $S_{1t} = \sum\limits_{i = 1}^{M_1(t)}{Y_i}$ and $S_{2t} = \sum\limits_{i = 1}^{M_2(t)}{Z_i}$, are both c.p. processes, representing the total claims and returns until time $t$, respectively, and $M_{1}(t)$ and $M_{2}(t)$ are homogeneous Poisson processes with parameters $\lambda_1 > 0$ and $\lambda_2 > 0$. The claim size is determined by the cumulative distribution function (c.d.f.) $F_Y(\cdot)$ and the probability density function (p.d.f.) ${f_Y}(\cdot)$ of independent and identically distributed (i.i.d.) positive random variables (r.v.s) $\{Y_i\}_{i = 1}^{\infty}$. The random return is given by the c.d.f. $F_Z(\cdot)$ and the p.d.f. ${f_Z}(\cdot)$ of the positive r.v.s $\{Z_i\}_{i = 1}^{\infty}$. Define $M_1(t) = \sup\{ j:S_{11} + S_{12} +\cdots + S_{1j} \le t\}$ and $M_2(t) = \sup\{ j:S_{21} + S_{22} +\cdots + S_{2j} \le t\}$, where inter-claim times $\{S_{1j}\}_{j = 1}^{\infty}$ and inter-return times $\{S_{2j}\}_{j = 1}^{\infty}$ follow the exponential distribution of intensity $\lambda_1$ and $\lambda_2$, respectively.

In reality, to protect the interests of the manager and the insured, the manager needs to have a reasonable plan for the surplus funds. Under normal circumstances, insurance companies generally take a portfolio of risk and risk-free investments for surplus funds^[22]. As investment income becomes a larger share of insurance company's total revenue, we need to take into account investment ratio factors. Therefore, suppose that the manager uses part of the surplus funds for risk-free investment and the other part for risk investment. In that way, risk-free investment $\{R_t\}_{t\geq0}$ satisfies

where $r$ is the interest rate on a risk-free asset, so obviously $r$ should be greater than zero. Risk asset $\{Q_t\}_{t\geq0}$ is defined as

where $\{W_t, t\geq0\}$ is a standard Brownian motion, and $\sigma$ and $a$ represent the volatility and expected rate of return of risk assets, respectively, both of which are greater than zero. So the risk asset process $\{{Q}_t\}_{t\geq0}$ satisfies

Let $q \in (0, 1)$ represent the proportion of the insurance company's surplus invested in risky assets, and then $1-q$ represents the proportion invested in risk-free assets. So $U(t)$ satisfies

where $\xi = (a+\frac{1}{2}\sigma^2)q +(1 - q)r$, $U(t-)$ is the left limit of $U(t)$ at $t$, and the loading condition to ensure that the formula holds is $c+\lambda_2{E[Z_1]} > \lambda_1{E[Y_1]}$.

We consider the dividend problems of the above model under the dividend strategy: when $U(t)$ is greater than threshold $b$, dividends are paid consecutively in $\alpha$, where $\alpha$ is constant and greater than zero; when $U(t)$ is greater than zero and less than $b > 0$, no dividends are paid; and when $U(t)$ is less than zero, bankruptcy occurs at this time $($but, in practice, the state of this moment may not be observed and therefore is still meaningful in the short term$)$. Combined with Eq (2.5), the surplus process with threshold $b$ is represented by $\{U_b(t), t\geq 0\}$, and $\{U_b(t), t\geq 0\}$ satisfies

where $V(Q, R, q, t) = (1 - q)\frac{{\rm{d}}{R_t}}{{R_t}} + q\frac{{\rm{d}}{Q_t}}{{Q_t}}$.

Let the cumulative dividend paid until the ruin time $t$ be $\mathscr{D}(t)$, and $T_b = \inf\{t:U_b(t) \leq 0 \}$ is the ruin time. The present value of accumulated dividends before the ruin time $T_b$ is $\mathscr{D}_{u, b}$, so

where $\delta$ is the interest force and is greater than zero, and $I(\cdot)$ stands for the indicator function. According to the above definition, it is not difficult to derive $0 < \mathscr{D}_{u, b} < \frac{\alpha}{\delta}$, which provides convenience for the subsequent derivation of the boundary of the IDEs. For $u \in \mathbb{R}$, the expectation of $\mathscr{D}_{u, b}$ is represented by

It should be emphasized that the surplus can be observed randomly in this paper. In practice, however, the executive director of an insurance company randomly reviews the balance of the company's books to determine whether dividends are being paid or whether it is ruined (e.g., ^[23,24,25]). Suppose $\{T_j\}_{j = 0}^{\infty}$ is a series of discrete time points of the moments of observing surplus, where $T_j$ is the $j$th observation time. In addition, we stipulate that $T_0 = 0$ and $T_{j^*}$ is the time when the company goes to ruin, where $j^* = {\rm{inf}}\{j\ge 1: M(j) \leq 0 \}$. Suppose $\{S_j\}_{j = 0}^{\infty}$ is an i.i.d sequence, where $S_j = T_j-T_{j-1}$ is the $j$th observation interval and $S_j$ are positive r.v.s, which are subject to an exponential distribution of intensity $\gamma > 0$. Suppose $\{{Y_i}\} _{i = 1}^\infty$, $\{{Z_i}\} _{i = 1}^\infty$, $\{M_1(t)\}_{t\geq0}$, $\{M_2(t)\}_{t\geq0}$, $\{W_t, t\geq0\}$, and $\{S_j\}_{j = 0}^{\infty}$ are independent of each other. Let the surplus level of the $j$th observation be $M(j) = U(T_j)$, and combine (2.5) to obtain

3.   IDEs of $V(u; b)$

In this section, our work is to give the IDEs of e.d.d.p. $V(u; b)$. Before we begin, we need to discuss the range of values of $u$, considering a time interval $(0, dt]$. If a claim occurred before observation, it is possible that $U_b(t) < 0$ was not observed. Therefore, the range of values of $u$ extends to the entire field of real numbers. In addition, it is not difficult to find that for different initial surplus $u$, $V(u; b)$ behaves differently. For convenience, let us set

Here are the following conclusions.

Theorem 3.1. For $u\in\left(-\infty, 0 \right]$, $V_1(u; b)$ satisfies

For $u\in\left(0, b \right]$, $V_2(u; b)$ satisfies

and for $u\in(b, \infty)$, $V_3(u; b)$ satisfies

The following boundary conditions are satisfied

Proof. Consider an infinitesimal interval (0, ${\rm{d}}t$], and discuss whether claims and benefits occur or not. The cumulative distribution function of $Y_i$ and $Z_i$ is continuous. For $u\in\left(-\infty, 0 \right]$,

and for $u\in\left(0, b \right]$,

amd for $u\in(b, \infty)$,

where

According to the Itô formula, we get

where

and $o({\rm{d}}t)$ stands for the infinitesimal of higher order ${\rm{d}}t$.

Substitute Eqs (3.11)–(3.14) into Eqs (3.8)–(3.10), respectively. Divide both sides of the equation by ${\rm{d}}t$ and let ${\rm{d}}t$ approach zero infinitely. According to the properties of higher order infinitesimals and some careful calculation, we can get the $IDEs$ (3.1)–(3.3).

With further analysis, if the initial surplus $U_{b} < 0$, the ruin occurs immediately, at which time no dividend is paid; then $T_{b} = 0$. If $0 < U_{b} < b$, then the ruin did not occur and the dividend is always paid at rate $\alpha$. If $U_{b} > b$, then the shares are always paid at rates $\alpha-c$, so $T_{b} = \infty$. □

Remark 3.1. Referring to the analysis of Albrecher ^[26], we can also find that $V(u; b)$ is not differentiable when $u = 0$ in general. Similarly, to fully describe the solution of Theorem 3.1, we also use $V_1(0-;b) = V_2(0+;b)$ and $V_1(b-; b) = V_2(b+; b)$, and the boundary conditions (3.4) and (3.5).

4.   Sinc asymptotic analysis

The sinc numerical method was proposed by James H. Wilkinson in the 1950s and developed by Frank Stenger in the 1990s. Frank Stenger summarized his work results in ^[27], which caused a great response in various fields (e.g., ^[28,29]). The real solutions to Eqs (3.1)–(3.3) are theoretically difficult to obtain. Therefore, we changed the angle, tried to obtain the a.s. by a numerical method, and then carried out an error analysis. Nowadays, the commonly used numerical methods for solving integral differential equations include the RK-Fehlberg method, the sinc method, the Runge-Kutta method, the Adams method, and so on. The sinc numerical method has high accuracy and good convergence when the sampling interval is small enough, which makes it perform well in high-precision numerical results. At the same time, the sinc method has an adaptive sampling interval. When the sampling interval is small, the sinc method can accurately reflect the details of the original function, to achieve high-precision numerical calculation. When the sampling interval is large, the sinc method can effectively smooth the function and avoid the ringing effect ^[30] in the interpolation process. Therefore, we also use this numerical method here.

4.1.   Approximate solution of $V(u; b)$

Since the domain of $u$ is the entire real axis, in order to construct approximations on $\mathbb{R}$, we consider conformal mappings. According to Algorithm 1.5.18 of Stenger^[27], we define an injective mapping from $\mathbb{R} \to \mathbb{R}$

where $z \in \mathbb{R}$. Define the grid point ${z_k}$ of sinc as

where $k \in \mathbb{Z}$, $h > 0$. Based on the sinc method, the basis function of $z\in \mit\Gamma$ on the interval $(-\infty, \infty)$ is given by the following composite function

Following the steps of the sinc method, we arrange Eqs (3.1)–(3.3) into the following integral differential

By the nature of convolution, Eq (4.2) is rewritten as

According to formulas (3.4) and (3.5), and Definition 1.5.2 in reference^[27], we have

where $\zeta(u) = {\rm{e}}^{\phi(u)} = {\rm{e}}^u$, when $t_1\rightarrow {-\infty}$, $t_2\rightarrow {\infty}$. Set

and then $W(u) \in {\mathscr{L}_{\tilde \alpha, \tilde \beta }(\delta)}$, where $\mathscr{L}_{\tilde{\alpha}, \tilde{\beta}}(\delta)$ is the function space for the sinc approximation over the finite interval $(\tilde{\alpha}, \tilde{\beta})$ (p. 72 in ^[27]).

When $u \to -\infty$, $u \to\infty$

Substituting (4.5)–(4.7) into (4.3), by simple calculation, we have

where $\mu_0(u) = \frac{(qu\delta)^2}{2}$, $\mu_1(u) = \xi u+c-\alpha{{{I}}_{(u > b)}}$, $\mu_2(u) = -(\delta+\lambda_1+\lambda_2),$

When $h > 0$, define the sinc grid point as

Then consulting reference ^[27], according to Theorem 1.5.13, Theorem 1.5.14, and Theorem 1.5.20, we can get

where $A$ and $B$ are resemble diagonal matrices $\Lambda$, with $A_{ij}$ and $B_{ij}$ denoting the elements at $(i, j)$ in $A$ and $B$, respectively. The approximate value of $W({u_j})$ is expressed by $U_j$.

Substituting (4.13)–(4.15) into Eq (4.8), replacing the integral term on the right side of Eq (4.8) with Eqs (4.13)–(4.15), and replacing $u$ with $u_k$ for $k = n_2, \cdots, n_1$, where $u_k$ is the sinc grid point, we have

where

Substituting (4.17)–(4.19) into Eq (4.16), we have

Multiplying Eq (4.20) by $\frac{{h^2}}{[\phi^{'}({u_k})]^2}$, we have

Since

formula (4.21) can be turned into

Set ${I^{(m)}} = [\delta _{kj}^{(m)}]_{(n_2+n_1+1)\times(n_2+n_1+1)}$, and $m = -1, 0, 1, 2$. We rewrite Eq (4.22) as

where

So solving Eq (4.23), we get the expression of the approximate solution (a.s) of (4.5):

The meanings of the symbols mentioned in the above process are shown in Table 1.

4.2.   Error analysis

In the previous subsection, we obtained an inexact solution $($e.s.$)$ of the IDEs by using the sinc method. Therefore, in this section, we need to analyze the discrepancy between the a.s.s and the actual solutions. According to references ^[27,31], we find an upper bound of the error. Moreover, in reality, $u$ is non-negative. Therefore, in this subsection, our discussion takes place under the condition $u > 0$. Multiply $\frac{1}{\mu_{0}(u)}$ by both sides of Eq (4.8), and we set

so we have

where $\tilde{\mu_1}(u) = \frac{\mu_1(u)}{\mu_0(u)}$, $\tilde{\mu_1}(u) = \frac{\mu_2(u)}{\mu_0(u)}$.

Assumption 4.1. Let $\tilde{\mu_1}(u) / {\zeta^{'}}$, $1/(\zeta^{'})^{'}$, and $\tilde{\mu_2}(u)/(\zeta^{'})^2$ be elements of $\mathscr{W}^{\infty}(\mathscr{D})$, and we are given that $G/(\zeta^{'})^2 \in \mathscr{L}_{\hat{\alpha}(\mathscr{D})}$ and Eq (4.25) possess a single solution $W \in \mathscr{L}_{\hat{\alpha}}(\mathscr{D})$.

In the above assumption, $\mathscr{W}^{\infty}(\mathscr{D})$ represents the family of all functions of $W(u)$ that are analytically and uniformly bounded by $\mathscr{D}$, and $\mathscr{L}_{\hat{\alpha}}(\mathscr{D}) = \mathscr{L}_{\hat{\alpha}, \hat{\alpha}}(\mathscr{D})$.

Theorem 4.2. If the aforementioned assumption is true, $W$ represents the e.s. of Eq (4.25), $\tilde{W}$ represents the a.s. of Eq (4.24), and $\boldsymbol{U} = (U_{-n_2}, \ \cdots, \ U_{n_1})^{T}$ represents the e.s. of Eq (4.23). So there is a constant $\tilde{c} > 0$, and different from $N$, such that

Proof. Let

By using the triangle inequality, it is easily obtained that

Based on Theorem 4.4 in ^[31], there is a constant $c^* > 0$, and different from $N$, that according to Assumption 3.1, $W \in \mathscr{L}_{\hat{\alpha}(\mathscr{D})}$, and we have

For inequality (4.28), $|\mathscr{O}_{N}(u)-\tilde{W}(u)|$ fulfills

Similar to Theorem 3.8 in ^[31], if $u\in \Gamma$, then $\sum_{k\in\ \mathbb{Z}}|C(j, h)\circ\zeta(u)|^2 = 1$, and we can obtain

where $\textbf{W} = (W_{-n_2}, \ \cdots, \ W_{n_1})^{T}$ and $c^{**} > 0$ that is not dependent on $N$. Let us take $\tilde{c} = \rm{max}\{c^{*}, c^{**}\}$, and therefore, inequality (4.25) is obtained by formulas (4.27) − (4.31). □

Through formulas (4.4), (4.24), and (4.25), we get

5.   Numerical example

In this subsection, we provide specific numerical examples to demonstrate the effectiveness of the sinc method, and study the effects of investment ratio $q$ and fluctuation parameter $\sigma$ on the expected discounted dividend payout under exponential and lognormal distributions, respectively.

5.1.   The exponential distribution

All numerical examples in this section are assumed to be obtained under

and

Then,

and

Formulas (4.8) and (4.11) are converted to

and

Next, we examine how parameters $q$ and $\sigma$ affect $V(u; b)$. If not specified, the following example parameters are set as follows: $\delta = 0.06, \ \tilde{\alpha} = \frac{\pi}{4}, \ \tilde{\beta} = \frac{\pi}{4},$ $a = 0.6, \ c = 0.3, \ r = 0.05, \ \alpha = 0.2$, $d = \frac{\pi}{4}, \ N = 15, \ \lambda_1 = 1, \ \lambda_2 = 2, \ \eta_1 = 3, \ \eta_2 = 1$.

Example 5.1. The effect of the investment ratio $q$ on the e.d.d.p. is considered in the case of the exponential distribution of claims and returns. Set parameter $\sigma = 0.2$. As depicted in Figure 1, it becomes evident that as the proportion of surplus invested in risk assets increases, the corresponding fluctuation of $V(u; b)$ also increases. The value of $V(u; b)$ when $q$ changes is presented in Table 2 partially.

Example 5.2. The effect of volatility parameter $\sigma$ on the e.d.d.p. is considered in the case of the exponential distribution of claims and returns. Set parameter $q = 0.2$. As depicted in Figure 2, the greater the change of parameter $\sigma$, the greater the fluctuation of the curve corresponding to $V(u; b)$. Partial data is presented in Table 3.

As can be seen from Examples 5.1 and 5.2, the impact of two factors on the e.d.d.p. is considered: the proportion of risk investment $q$ and the volatility of risk assets $\sigma$. First, when a company invests a higher proportion of its surplus in risky assets, the dividend payout is higher, but also more volatile, while the dividend payout is more stable when the investment ratio is lower. This means high risk, high reward, danger, and opportunity. In addition, if the proportion of risk investment is fixed, choosing investment products with more volatile risk assets will bring higher profits, but also bear higher risks. On the contrary, they will earn lower profits and take lower risks. This is in line with reality.

5.2.   The lognormal distribution

In this section, it is assumed that $f_{Y}(y)$ and $f_{Z}(z)$ obey a lognormal distribution of parameter $(\eta_{3}, \ 2v_{1}^2)$ and $(\eta_{4}, \ 2v_{2}^2)$, respectively, where $\eta_{3} = \ln y$ and $\eta_{4} = \ln z$, and $2v_{1}^2$ and $2v_{2}^2$ represent the variance, so that $f_{Y}(y)$ and $f_{Z}(z)$ are defined as

Then,

and

Therefore, the formulas (4.8) and (4.11) can be rewritten as:

and

The next example is given in real condition: $\delta = 0.06, \ \tilde{\alpha} = \frac{\pi}{4}, \ \tilde{\beta} = \frac{\pi}{4}, \ a = 0.5, \ c = 0.4, \ r = 0.06, \ \alpha = 0.1, \ d = \frac{\pi}{4}, \ N = 10, \ $ $\lambda_1 = 1$, $\lambda_2 = 2$, $\eta_3 = 3$, $\eta_4 = 1$, $v_2 = 0.03$, $v_1 = 0.03$.

Example 5.3. In the case of a lognormal distribution of claims and returns, let us discuss the effect of investment ratio $q$ on $V(u; b)$. Set parameter $\sigma = 0.2$. It is not difficult to see from Figure 3 that when a company invests more surplus into risk assets, the growth of its expected discounted dividend curve experiences significant fluctuations. Partial data is presented in Table 4.

Example 5.4. In the case of a lognormal distribution of claims and returns, let us discuss the effect of investment ratio $\sigma$ on $V(u; b)$. Set parameter $q = 0.2$. It is not difficult to see from Figure 4 that when the company chooses a product investment with greater risk fluctuation, the growth of its expected discounted dividend curve exhibits substantial variability. Partial data is presented in Table 5.

From Examples 5.3 and 5.4, it can be seen that parameters $q$ and $\sigma$ have different effects on the e.d.d.p. $V(u; b)$ under a lognormal distribution of claims and returns. Other parameters being equal, the expected discounted dividend payout curve fluctuates more when a company invests a larger proportion of its earnings or invests in risky products with a higher freezing rate. It should be noted that when the claim amount and income follow the lognormal distribution, $V(u; b)$ shows a higher sensitivity to the above parameter changes.

6.   Conclusions

We explore a model with two-sided jumps, incorporating random observations and a dividend barrier strategy. By referring to the existing relevant literature, we find that the existing research is the classic model with a dividend strategy or the two-sided jump risk model. We want to know the situation of the dividend barrier strategy under double risk. According to this idea, through the literature review, we find that the model has very important practical significance. At the same time, we find that no scholars have introduced random observation into this model, but this is exactly what is for random observation in real life. In the process of research, we also find that there is no closed solution to the integral differential equation of this model after introducing random observation. To solve this problem, we obtained an a.s. by the sinc numerical method and analyzed the upper limit of the error. Perhaps one day in the future, we will have a better way to find the e.s. to this model.

Author contributions

Chunwei Wang: Methodology, supervision, resources, funding acquisition, writing-review & editing; Shaohua Li: Methodology, software, visualization, writing-original draft; Shujing Wang: Software, visualization; Jiaen Xu: Methodology, validation. All authors have read and agreed to the published version of the manuscript.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research was supported by the National Natural Science Foundation of China (No. 71801085). The authors would like to thank the referees for their valuable comments and suggestions.

Conflict of interest

All authors declare no conflicts of interest in this paper.