1. Introduction
In this paper, the surplus process of an insurance company is described by the perturbed compound Poisson model
where $u\geq 0$ is the initial surplus, and $c>0$ is the constant premium rate per time. The aggregate claims process $S_t=\sum _{n=1}^{N_t} X_n$ is a compound Poisson process, where $\{N_t\}_{t\geq 0}$ is a Poisson claim number process with constant arrival rate $\lambda >0$, and the individual claim sizes $\{X_n\}_{n\geq 1}$ form a sequence of positive valued i.i.d. random variables with common probability density function $f_X$ and mean $\mu _X=\int _0^\infty x f_X(x)\,dx$. Finally, $\{B_t\}_{t\geq 0}$, independent of $\{S_t\}_{t\geq 0}$, is a standard Brownian motion starting from zero, and $\sigma >0$ is a volatility parameter.
The above perturbed compound Poisson model was first proposed by Gerber [Reference Gerber11], and it has been studied by many authors; see, for example, Gerber and Landry [Reference Gerber and Landry12], Wang [Reference Wang26], Tsai [Reference Tsai23,Reference Tsai24], Tsai and Willmot [Reference Tsai and Willmot25], Chiu and Yin [Reference Chiu and Yin7] and references therein. In this paper, we consider the model (1.1) modified by a constant barrier dividend strategy. Given a finite barrier of level $b>0$, we assume that whenever the surplus process reaches level $b$, dividends are paid off continuously such that the surplus stays at level $b$ until it becomes less than $b$. Let $U^b_t$ denote the modified surplus process under the above barrier dividend strategy, and let $\tau _b=\inf \{t\geq 0: U_t^b\leq 0 \}$ be the time of ruin. The present value of total dividends paid before the ruin time $\tau _b$ is given by
where $\delta >0$ is the force of interest for valuation, and $D(t)$ is the aggregate dividends paid by time $t$. Given the initial surplus level $0\leq u\leq b$, the expected present value of total dividend payments before ruin is defined by
Furthermore, we are interested in the expected discounted penalty function associated with model $U^b$, which is defined by
where $w(\cdot )$ is a nonnegative penalty function of the deficit at ruin. Note that ruin can be caused either by a claim or oscillation due to the Brownian motion, then we have the following decomposition of $\phi (u;b)$,
where
are respectively the Laplace transform of the ruin time when ruin is caused by oscillation, and the expected discounted penalty function when ruin is due to a claim.
Dividend problem was first proposed by de Finetti [Reference de Finetti10] in a binomial model, and since then it has been studied in a number of papers. For example, see Lin et al. [Reference Lin, Willmot and Drekic17], Dickson and Waters [Reference Dickson and Waters8], Li and Garrido [Reference Li and Garrido16], Albrecher et al. [Reference Albrecher, Hartinger and Tichy2] and Li [Reference Li15] to name a few. We know that a common assumption in these papers is that the probability characteristics of the model $U^b$ are known, so that analytic approach such as differential equation, Laplace transform, renewal theory can be used to study the functions $V(u;b)$, $\phi _d(u;b)$ and $\phi _c(u;b)$. However, usually the quantities such as the Poisson intensity $\lambda$, the diffusion parameter $\sigma$ and the claim size density $f_X$ are all unknown. Hence, it is very interesting to consider the statistical estimation of the dividend problems given that the random samples on claim number, individual claim sizes and the surplus levels are available. When the probability characteristics of the surplus process are unknown, some estimation problems have been considered under the model without dividend payments. For example, Zhang and Yang [Reference Zhang and Yang31,Reference Zhang and Yang32] proposed some nonparametric estimators of the ruin probability by some Fourier-based methods. The same estimation problem has also been considered by Cai et al. [Reference Cai, Chen and You3] and Cai and You [Reference Cai and You4] by numerical Laplace inversion transform, respectively. The expected discounted penalty functions are estimated by Shimizu [Reference Shimizu19,Reference Shimizu20] and Shimizu and Zhang [Reference Shimizu and Zhang21] under different risk models.
Recently, Xie and Zhang [Reference Xie and Zhang27] study the statistical estimation of dividend problems in the classical compound Poisson model. In their paper, the Fourier-cosine (COS) method was first used to estimate the expected present value of dividend payments before ruin, then the expected discounted penalty function was estimated by using the dividends-penalty identity. We remark that the COS method was first proposed by Fang and Oosterlee [Reference Fang and Oosterlee9] to price European options, and since then, it has been widely used in the field of financial mathematics. In insurance risk theory, Chan et al. [Reference Chau, Yam and Yang5,Reference Chau, Yam and Yang6] applied the COS method to approximate the ruin probability and the expected discounted penalty function, respectively; Zhang used the COS method to approximate the density function of the time to ruin; Yang et al. [Reference Yang, Su and Zhang28] applied a two-dimensional COS method to estimate the discounted density function of the deficit at ruin; Lee et al. [Reference Lee, Li, Liu, Shi and Yam14] studied the finite time ruin probabilities by the COS method. On the estimation of the expected discounted penalty function with dividend payments, we also would like to point out the work Strini and Thonhauser [Reference Strini and Thonhauser22]. In their paper, ruin problems in a renewal risk model with a general surplus dependent premium rate are studied. Note that both barrier and band type strategies can be expressed by some special settings of their premium rate function. Different from Xie and Zhang [Reference Xie and Zhang27], they first study the expected discounted penalty function by numerical solution of a partial-integro-differential equation, and then consider the statistical estimation through the estimated model parameters.
In this paper, we shall follow the approach in Xie and Zhang [Reference Xie and Zhang27] to estimate the functions $V(u;b), \phi _d(u;b)$ and $\phi _c(u;b)$ under the model $U^b$. Note that different from Xie and Zhang [Reference Xie and Zhang27], in this paper, we add an independent Brownian motion to describe the uncertainty of the surplus flow of the insurance company. This additional noise puts an obstacle for constructing the estimators of the interesting functions from the random samples on the individual claim sizes and claim number. As is shown in Section 4, the information on dividend payments has to be used to construct the estimators. When studying the consistency properties of our estimators, we need to study the uniform convergence properties of some functions related to the Laplace exponent of the process $U^\infty _t$, and this is more difficult to deal with than the compound Poisson model in Xie and Zhang [Reference Xie and Zhang27]. The remainder of this paper is organized as follows. In Section 2, we present some preliminaries on $V(u;b)$, $\phi _{c}(u;b)$ and $\phi _{d}(u;b)$. In Section 3, we show how to approximate $V(u;b)$, $\phi _{c}(u;b)$ and $\phi _{d}(u;b)$ by the COS method. In Section 4, we show how to estimate those functions. The consistency properties for the estimators are studied in Section 5. Finally, in Section 6, we present some numerical results to illustrate the performance of our method.
2. Some preliminaries
In this section, we recall some known results on the expected present value of dividend payments before ruin and the expected discounted penalty functions. Throughout this paper, for any function $f$ defined on the positive real line, we denote its Fourier transform and Laplace transform by
where $s$ is such that the above integrals are well defined.
Note that $U^\infty _t$ is a spectrally negative Lévy process, and its Laplace exponent is defined by
For each $q\geq 0$, let
be the right inverse of $\psi _U$. For the special case $q=\delta$, we shall put $\rho =\rho _\delta$ for simplicity.
For each $q, x\geq 0$, let $W_q(x)$ denote the $q$-scale function associated with the process $U^\infty _t$, which is a strictly increasing and continuous function with Laplace transform given by
We can extend $W_q$ to the whole real line by setting $W_q(x)=0$ for $x<0$. The $q$-scale function and its various extensions have many applications in the fields related to the spectrally negative Lévy processes. For a thorough introduction to the $q$-scale function and some of its applications, we refer the interested readers to the survey Kuznetsov et al. [Reference Kuznetsov, Kyprianou and Rivero13].
In our applications, it is useful to introduce the following auxiliary function,
It follows from Zhang and Cui [Reference Zhang and Cui30] that $h(x)$ is the probability density function of the random variable $U^\infty _0-U^\infty _{\mathbf {e}_\delta }$, where $\mathbf {e}_\delta$, independent of $U^\infty$, denotes an exponential random variable with rate $\delta >0$.
For the expected present value of total dividend payments before ruin, it follows from Albrecher and Gerber [Reference Albrecher and Gerber1] that
Then, by formula (2.9), in Xie and Zhang [Reference Xie and Zhang27], we obtain
where
Let $\tau _\infty =\inf \{t\geq 0: U^\infty _t\leq 0\}$ denote the ruin time of the model $U^\infty$, and accordingly define
to be the Laplace transform of ruin time when ruin is caused by the Brownian motion, and the expected discounted penalty function when ruin is due to a claim. For the expected discounted penalty functions $\phi _d(u;b)$ and $\phi _c(u;b)$, they can be expressed via the dividends-penalty identities as follows,
For example, see Lin et al. [Reference Lin, Willmot and Drekic17].
3. The COS approximation
In this section, we shall use the COS method to approximate $V(u;b)$, $\phi _d(u;b)$ and $\phi _c(u;b)$. Recall that for an integrable function $f$ with a finite support $[a_1, a_2]$, it has the following cosine series expansion,
where $\mathop {{\sum }'}$ denotes a summation with its first term weighted by a half, and the cosine coefficients are given by
where $\text {Re}(\cdot )$ means taking real part and $i=\sqrt {-1}$ is the imaginary unit.
If $f$ is an integrable function supported on the positive real line $[0, \infty )$, we can expand $f$ on a closed interval $[0, a]$, and the COS method suggests that, for large enough $a$, the COS coefficients can be approximated as follows,
Hence, we have
where $K$ is a large integer applied to truncate the infinite series.
The following lemma is proved in Xie and Zhang [Reference Xie and Zhang27], which gives the approximation error for the COS method.
Lemma 1. For real-valued integrable function $f$ supported on $[0, \infty )$, suppose that $|f'(0+)|<\infty$, $|f'(a)|<\infty$, and $\int _0^\infty |f''(y)|\,dy<\infty$. Then, for some positive constants $c_1$ and $c_2$, we have
Remark 1. Lemma 1 shows that the uniform approximation error of the COS method in the interval $[0, a]$ depends on the parameters $a$ and $K$. Note that the parameter $a$ is in fact an integration domain truncation parameter in the COS method, which is used in the approximation of the COS coefficients by the Fourier transform $\mathcal {F}f$. The parameter $K$ is used to truncate the COS series expansion formula. Usually, large $K$ will result in good approximation. The upper bound in (3.5) consists of two terms. The first error term $aK^{-1}$ means that larger $K$ can yields better approximation, but large $a$ may slow down the convergence rate. Usually, we should choose an appropriate truncation parameter $a$ to effectively capture the support of the objective function $f$. In the literature, it is usually selected according to some cumulant-based methods (see [Reference Fang and Oosterlee9]). The second error term $Ka^{-1}\int _a^\infty |f(y)|\,dy$ means that larger $a$ can yield better approximation, since larger $a$ usually can result in more accurate approximation of the COS coefficients. However, this error term is increasing w.r.t. $K$, since larger $K$ means more COS coefficients have to be approximated. Usually, the term $a^{-1}\int _a^\infty |f(y)|\,dy$ can decay very fast w.r.t. $a$, especially when the tail of the function $f$ converges to zero at the exponential rate. Overall, the first error can dominate the second error. Hence, for a fixed parameter $a$, the approximation error in (3.5) can achieve order $O(K^{-1})$.
Remark 2. Usually, the decay rate of the COS coefficients depends heavily on the smoothness of the objective function $f$. If $f$ is infinitely times differentiable, the series $\{A_{f,k}\}$ will show exponential convergence; otherwise, it will yield algebraic convergence. See, for example, Fang and Oosterlee [Reference Fang and Oosterlee9]. In our paper, note that the function $f$ that we approximate usually has support $[0, \infty )$, then using integration by parts we can find that
which yields that $A_{f,k}=O(k^{-2})$ for fixed $a$. Furthermore, if $f$ is $2n$ times differentiable and $f^{(2n)}$ is integrable, and $f^{(j)}(0+)=0$ for $j=1,3,\ldots, 2n-1$, then repeatedly using integration by parts we can prove that
Hence, if the derivatives $f^{(1)}, f^{(3)}, \ldots, f^{(2n-1)}$ all have exponential decay rate, we can first choose large $a$ to ignore the term $a^{-1}\sum _{j=1}^{n} ({a}/{k})^{2j} f^{(2j-1)}(a)$, then for such $a$, $A_{f,k}=O(k^{-2n})$, which can further lead to faster convergence for the COS approximation.
Remark 3. Although faster convergence can be achieved under the additional conditions in Remark 2, these conditions are usually not satisfied by the functions in this paper. For example, the objective functions to be approximated in this paper may take the form of exponential functions or their finite mixture, and we can only obtain algebraic convergence. However, in this case, spectral filters can be applied to accelerate the convergence rate according to Ruijter et al. [Reference Ruijter, Versteegh and Oosterlee18]. Recall that filtering is carried out in Fourier space and the idea is to pre-multiply the expansion coefficients by a decreasing function. It should be noted that filtering does not add any significant computational costs, and usually can improve the algebraic convergence. In Section 6, we shall illustrate more details on filtering.
In order to use the COS method to approximate $V(u;b)$, we should first approximate $h_+$ and $g_+$. It follows from Xie and Zhang [Reference Xie and Zhang27] that
and
Since $W_\delta (0)=0$ as $\sigma >0$, we have from (2.3) that $h_+(0)= {\delta } / {\psi _U'(\rho )}$, and this yields
and
Now using the closed-form Fourier transforms $\mathcal {F}h_+(s)$ and $\mathcal {F}g_+(s)$, the functions $h_+$ and $g_+$ can be approximated by the COS method as follows,
and
where the COS coefficients are given by
As for the expected present value of dividend payments before ruin, using formula (2.5) we can approximate it as follows,
where we take $a>b$.
Remark 4. Suppose that $|h'(0+)|<\infty,\ |h''(0+)|<\infty$, $|h'(a)|<\infty$, $|h''(a)|<\infty$ and
Furthermore, let $C$ be a positive generic constant that may vary at different steps. Then by Lemma 1 we can obtain
where
Next, we consider how to use the COS method to approximate the expected discounted penalty functions $\phi _d(u;b)$ and $\phi _c(u;b)$. To this end, we use the dividends-penalty identities given in (2.6) and (2.7). It is easily seen that we should use the COS method to approximate $\phi _{d}(u;\infty ), \phi _c(u;\infty )$ and their first-order derivatives.
For $\phi _d(u;\infty )$, by Zhang [Reference Zhang29], we know that its Laplace transform and Fourier transform are given by
Then, the Laplace transform of the derivative $\phi _d'(u;\infty )$ is given by
from which we obtain
For $\phi _c(u;\infty )$, its Laplace transform and Fourier transform are given by
where $\omega (u)=\int _u^\infty w(x-u)f_X(x)\,dx$, $u\geq 0$. Furthermore, since $\phi _c(0;\infty )=0$, we have
and
Now, using the COS method, we have
where, for $k=0,1,\ldots$, the COS coefficients are given by
Finally, by the dividends-penalty identities (2.6) and (2.7), we obtain
and
Remark 5. Again, the approximating error can be obtained by Lemma 1. Suppose all the conditions in Remark 4 hold true. Furthermore, suppose that $|\phi _d'(0+;\infty )|<\infty,\ |\phi _d''(0+;\infty )|<\infty$, $|\phi _d'(a;\infty )|<\infty$, $|\phi _d''(a;\infty )|<\infty$ and
then we can obtain from Lemma 1 and Remark 4 that
where $\mathcal {E}_h(a,K)$ is defined in (3.15), and
Similarly, suppose that $|\phi _c'(0+;\infty )|<\infty,\ |\phi _c''(0+;\infty )|<\infty$, $|\phi _c'(a;\infty )|<\infty$, $|\phi _c''(a;\infty )|<\infty$ and
Then, we can obtain
where
4. The estimation method
In this section, we show how to estimate $V(u;b)$, $\phi _d(u;b)$ and $\phi _c(u;b)$ from a random sample on the surplus process. We assume that the premium rate $c$ is known, but the Poisson intensity $\lambda$, the claim size density function $f_X$ and the diffusion parameter $\sigma$ are all unknown.
Suppose that the surplus process can be observed in the long time interval $[0, T]$, and the following random sample on the claim number and individual claim sizes during $[0, T]$ is available,
where $n_T$ denotes the number of claims received up to time $T$ and it is assumed to be strictly positive w.l.o.g. In addition, suppose that the surplus process $U^b_t$, the aggregate claims process $S_t$ and the aggregate dividends process $D(t)$ can be observed at a sequence of discrete time points so that the following sample is available
where $\Delta >0$ is a fixed sampling step satisfying $n\Delta =T$.
First, we use formula (3.13) to propose an estimator for the expected present value of dividend payments before ruin. It is easily seen that we need to estimate the following quantities,
Since $\rho$ is the positive root of equation $\psi _U(s)=\delta$, we should estimate the Laplace exponent $\psi _U(s)$ by formula (2.1). We estimate the Poisson intensity by $\hat {\lambda }={n_T}/{T}$, and estimate the Laplace transform $\mathcal {L}f_X(s)$ by $\widehat {\mathcal {L}f}_X(s) =({1}/{n_T})\sum _{j=1}^{n_T} e^{-sX_j}$. It remains to estimate the diffusion parameter $\sigma$. Since $U^b_t=U^\infty _t-D(t)$, then we have
For $k=1,\ldots, n$, set
which are available since the premium rate $c$ is known. Now, we estimate $\sigma ^2$ by
It is easily seen that $\hat {\lambda }, \widehat {\mathcal {L}f}_X(s)$ and $\hat {\sigma }^2$ are all unbiased estimators. Now using these estimators, we can estimate the Laplace exponent $\psi _U(s)$ and its derivative $\psi _U'(s)$ as follows,
Since $\rho$ is the positive root of equation $\psi _U(s)=\delta$, we define its estimator, denoted by $\hat {\rho }$, to be the positive root of the following estimating equation
Then, we can estimate $\psi _U'(\rho )$ by $\hat {\psi }_U'(\hat {\rho })$.
In order to estimate $\tilde {h}_+$ and $\tilde {g}_+$, we should first estimate the Fourier transforms ${\mathcal {F}h}_+(s)$ and $\mathcal {F}g_+(s)$. It follows from (3.7) and (3.9) that they can be estimated by the following estimators
Then, using formulas (3.10) and (3.11), we propose the following estimators for $\tilde {h}_+$ and $\tilde {g}_+$,
and
where the COS coefficients are given by
Now, using formula (3.13), we obtain the following estimator of $V(u;b)$,
Next, we consider how to estimate the expected discounted penalty function $\phi _d(u;b)$. To this end, we shall use the dividends-penalty identity (2.6). Note that the Fourier transform of $\phi _d(u;\infty )$ and its derivative can be estimated by
So we can estimate $\phi _d(u;\infty )$ and $\phi _d'(u;\infty )$ as follows,
where, for $k=0,1,\ldots,$
By the dividends-penalty identity (2.6), we propose the following estimator for $\phi _d(u;b)$,
Finally, we estimate the expected discounted penalty function $\phi _c(u;b)$. It is easily seen that
from which we obtain the following estimators for $\mathcal {L}\omega (\rho )$ and $\mathcal {F}\omega (s)$,
Now the Fourier transforms $\mathcal {F}\phi _c(s;\infty )$ and $\mathcal {F}\phi _c'(s;\infty )$ can be estimated by
Then $\phi _c(u;\infty )$ and $\phi _c'(u;\infty )$ are estimated as follows,
where, for $k=0,1,\ldots,$
By the dividends-penalty identity (2.7), we propose the following estimator for $\phi _c(u;b)$,
5. Consistency properties
In this section, we derive the consistency properties of our estimators when the observation interval $[0, T]$ is very large. First, it is easily seen that
Here, the notation $O_p(1)$ denotes a sequence that is bounded in probability; or more generally in our paper, for a given sequence of $\{R_{T,K,a}\}$, $O_p(R_{T,K,a})$ denotes a sequence that is bounded in probability at rate $R_{T,K,a}$. The following result on the estimator $\hat {\rho }$ is also well known. See, for example, Zhang [Reference Zhang29].
Lemma 2. Suppose that $c>\lambda EX$ and $EX^2<\infty$, then for each $\delta >0$, we have $\hat {\rho }-\rho =O_p(T^{-1/2})$.
Suppose the conditions in Lemma 2 hold true. For the estimator $\hat {\psi }_U'(\hat {\rho })$, note that
Lemma 2 and $\hat {\sigma }^2-\sigma ^2=O_p(T^{-1/2})$ imply that $(\hat {\sigma }^2\hat {\rho }-\sigma ^2\rho )=O_p(T^{-1/2})$. By Proposition 4 in Xie and Zhang [Reference Xie and Zhang27], we obtain
Hence, it follows from (5.1) that
Recall that $\rho$ is the root of equation $\psi _U(s)=\delta$, then we have
where $l(s)= ({\sigma ^2}/{2})(-is+\rho )+c+\lambda \int _0^\infty e^{isx}\int _0^x e^{-(\rho +is)y}\,dyf_X(x)\,dx$. Note that
Similarly, since $\hat {\rho }$ is the root of equation $\hat {\psi }_U(s)=\delta$, we have
where $\hat {l}(s)= ({\hat {\sigma }^2}/{2})(-is+\hat {\rho })+c +({1}/{T})\sum _{j=1}^{n_T} e^{is X_j} \int _0^{X_j} e^{-(\hat {\rho }+is)y}\,dy$. Furthermore, note that
which together with (5.6) gives
The following lemma shows the uniform convergence of ${1}/{(\hat {\psi }_U(-is)-\delta )}$ and ${s}/{(\hat {\psi }_U(-is)-\delta )}$, which plays an important role in studying the convergence of our estimators.
Lemma 3. Suppose that $c>\lambda EX$, $EX^2<\infty$ and $a=o(K)$. Then, we have
and
Proof. First, we study the uniform convergence of ${1}/{(\hat {\psi }_U(-is)-\delta )}$. By (5.3) and (5.5), we have for any real number $s$,
and by (5.6) and (5.8) we have for any real number $s$,
Recall the condition $c>\lambda EX$. In the remainder of this proof, suppose that $c -({1}/{T})\sum _{j=1}^{n_T} X_j>0$, since $({1}/{T})\sum _{j=1}^{n_T} X_j$ converges to $\lambda EX$ a.s.
Now using (5.12) and (5.13) and the following result
we can obtain
The convergence rate $\hat {\sigma }^2-\sigma ^2=O_p(T^{-1/2})$ implies that
It follows from Lemma 3 in Xie and Zhang [Reference Xie and Zhang27] that, for $a=o(K)$,
leading to
which together with (5.15), (5.16) yields (5.10).
Next, we prove (5.11). By (5.3) and (5.6), we have
which gives
It follows from (5.4) that
which together with (5.8) and the convergence rate $\hat {\sigma }^2-\sigma ^2=O_p(T^{-1/2})$ gives
By (5.5), (5.8) and (5.17), we obtain
In order to derive the estimation error of $\hat {V}(u;b)$, we should first consider the estimation errors of $\hat {h}_+$ and $\hat {g}_+$. We have the following results.
Proposition 1. Suppose that $c>\lambda EX$, $EX^2<\infty$ and $a=o(K)$. Then, we have
and
Proof. We only prove (5.22), since (5.23) can be proved similarly by using (5.11). First, it is easily seen that
where $\mathcal {S}=\{k\pi /a: k=0,1,\ldots, K\}$. For each $s\in \mathcal {S}$, we have
where
By Lemma 2 and (5.2), we can easily prove that
By Lemma 3, we have
Hence, we have
which together with (5.24) completes the proof.
By the convergence rate given in (5.2), Lemma 2 and Proposition 1, we can easily obtain the following result.
Proposition 2. Suppose that $c>\lambda EX$, $EX^2<\infty$ and $a=o(K)$. Then, we have
Next, we derive the convergence rate for the estimators of the expected discounted penalty functions. It follows from formulas (4.11) and (4.16) that we should first study the consistency properties of $\hat {\phi }_d(u;\infty )$ and $\hat {\phi }_d'(u;\infty )$.
Proposition 3. Suppose that $c>\lambda EX$, $EX^2<\infty$ and $a=o(K)$. Then, we have
Proof. First, we prove (5.26). By formulas (3.22) and (4.9), we can obtain
By (3.19) and (4.8), we obtain
Furthermore, by Lemma 3, we have
Similarly, we have
Hence, (5.29) gives
which together with (5.28) yields (5.26).
Next, we prove (5.27). Again, we can obtain
where
Furthermore, we have
By Lemma 3 and (5.17) and using the same arguments in (5.30) and (5.31), we can obtain
and
Hence, we obtain
Next, we derive the convergence rates for $\hat {\phi }_c(u;\infty )$ and $\hat {\phi }_c'(u;\infty )$.
Proposition 4. Suppose that $c>\lambda EX$, $EX^4<\infty$, $a=o(K)$, and
Then, we have
Proof. First, we prove (5.33). It is easily obtained that
where
Furthermore, by Xie and Zhang [Reference Xie and Zhang27], we have
Then, by Lemma 3, we can obtain
which together with (5.35) gives (5.33). The proof of (5.34) is similar, and we omit the detailed arguments.
Now by combining the results in Propositions 2–4, we can obtain the convergence rates of $\hat {\phi }_d(u;b)$ and $\hat {\phi }_c(u;b)$.
Proposition 5. Suppose that $c>\lambda EX$, $EX^2<\infty$ and $a=o(K)$. Then, we have
Proposition 6. Suppose that $c>\lambda EX$, $EX^4<\infty$, $a=o(K)$, and
Then, we have
Remark 6. Propositions 5 and 6 show that the convergence rate depend the parameters $a, K$ and $T$. Here, the parameters $a$ and $K$ appear in the COS approximation step, while the parameter $T$ stands for sample size in the statical estimation step. Usually, the parameter $a$ need not to be large enough, and we can use some cumulant-based methods for selection (see Section 6).
For the series truncation parameter $K$, although a larger value can lead to better COS approximation, it also implies that more COS coefficients have to be estimated in the statistical estimation step. Hence, the statistical estimation error is increasing w.r.t. $K$. We remark that the condition $a=o(K)$ is just used for theoretical derivation (see also [Reference Xie and Zhang27]), which means that the parameter $a$ should be dominated by the series truncation parameter $K$. This condition is satisfied in our paper, since when using the cumulant-based method to choose $a$, we usually have $a\leq 50$, but when choosing the truncation parameter $K$, we usually take $K=2^8, 2^9, 2^{10},$ etc., larger than $a$. Although we observe the error propagation w.r.t $K$, we can enlarge the sample size (i.e., take large $T$) to get satisfactory estimators. Anyway, if we simply look at the estimation of $\tilde {\phi }_d(u;b)$ and $\tilde {\phi }_c(u;b)$ (i.e., for fixed $K$ and $a$), Propositions 5 and 6 show that we can nearly obtain $O_p(T^{-{1}/{2}})$ convergence rate.
6. Numerical results
In this section, we present some numerical results to illustrate the performance of our method. All computations are performed in MATLAB on a desktop computer with Intel(R) Core(TM) i7-6700 CPU, 3.4 GHz and a RAM of 8 GB. In this part, we set $c=4$, $\lambda =3$, $\sigma =1$, $\delta =0.1$ and $b=30$. We consider the following three claim size density functions,
1. Exp(1): $f_X(x)=e^{-x},\ x>0$;
2. Erlang(2,2): $f_X(x)=4xe^{-2x},\ x>0$;
3. ${\rm Gamma}(\frac {3}{2},\frac {3}{2})$: $f_X(x)=\frac {3\sqrt {6x}e^{-{3x}/{2}}}{2\sqrt {\pi }},\ x>0$.
The choice of the truncation parameter $a$ is selected according to the following rule (see [Reference Fang and Oosterlee9]):
where
and $f$ is a probability density function. As in Xie and Zhang [Reference Xie and Zhang27], we set $L=10$. When computing $\tilde {h}_{+}$ and $\tilde {g}_{+}$, we take $\mathcal {F}f=\mathcal {F}{h}_{+}$; when computing $\tilde {\phi }_{d}$ and $\tilde {\phi }_{d}'$, we take $\mathcal {F}f=\mathcal {F}{\phi }_{d}$; when computing $\tilde {\phi }_{c}$ and $\tilde {\phi }_{c}'$, we take $\mathcal {F}f=\mathcal {F}{\phi }_{c}$.
For the expected present value of total dividend payments before ruin, we easily obtain the explicit formula by (2.4) when $f_X$ is exponential or Erlang(2,2). When $f_X$ is ${\rm Gamma}(\frac {3}{2},\frac {3}{2})$, in order to provide a benchmark, we compute the reference values by the COS method with truncation parameter $K=2^{13}$. For the expected discounted penalty function, we set $w \equiv 1$, then $\phi _d(u;b)$ is the Laplace transform of the time of ruin by oscillation, and $\phi _c(u;b)$ is the Laplace transform of the time of ruin due to a claim. When $f_X$ is exponential or Erlang(2,2), we can use Laplace inversion transform to compute the explicit formulas by (3.16) and (3.19) for $\phi _d(u;b)$ and $\phi _c(u;b)$, respectively. When $f_X$ is ${\rm Gamma}(\frac {3}{2},\frac {3}{2})$, we also compute the reference values by the COS method with truncation parameter $K=2^{13}$.
First, we study the approximation performance of the COS method. In order to measure the accuracy of our method, we consider average absolute errors for $\tilde {V}(u;b)$, $\tilde {\phi }_d(u;b)$ and $\tilde {\phi }_c(u;b)$, respectively, which are defined by
and
where $V^{{\rm ref}}(u;b)$, $\phi _d^{{\rm ref}}(u;b)$ and $\phi _c^{{\rm ref}}(u;b)$ are the reference values. We take $\mathcal {U}=\{0.1, 0.2, \ldots, 30\}$, $\mathcal {U}'=\{0.1, 0.2,\ldots, 5\}$ and $\mathcal {U}''=\{0.1, 0.2,\ldots, 15\}$, since $\phi _d(u;b)$ is very small for $u>5$, and $\phi _c(u;b)$ is very close to zero for $u>15$. The average absolute errors are reported in Table 1, where we consider the truncation parameter $K=2^p$ for $p=6, 7, 8, 9, 10$. All results in the following tables are rounded at the sixth decimal place. In Table 1, it can be easily seen that the absolute errors decrease with $p$ (or equivalently $K$) for $\tilde {V}(u;b)$, $\tilde {\phi }_d(u;b)$ and $\tilde {\phi }_c(u;b)$, which means that the large truncation parameter can reduce the error in our examples.
Next, we consider the effectiveness of our statistical estimation method, based on a random sample on the individual claim sizes $\{X_1, \ldots, X_{n_T} \}$ and claim number $n_T$. For the parameter $a$, we consider the following empirical version of the cumulant-based method,
where
When computing $\hat {h}_+$ and $\hat {g}_+$, we take $\widehat {\mathcal {F}f}=\widehat {\mathcal {F}h_+}$; when computing $\hat {\phi }_d$ and $\hat {\phi }_{d}'$, we take $\widehat {\mathcal {F}f}=\widehat {\mathcal {F}\phi _d}$; when computing $\hat {\phi }_c$ and $\hat {\phi }_{c}'$, we take $\widehat {\mathcal {F}f}=\widehat {\mathcal {F}\phi _c}$. Again, we set $L=10$.
Set the truncation parameter $K=2^{10}$ and the inter-observation interval $\Delta =1$. We perform $200$ experiments and consider average absolute errors for $\hat {V}(u;b)$, $\hat {\phi }_d(u;b)$ and $\hat {\phi }_c(u;b)$, respectively, which are defined by
and
where $\hat {V}_j(u;b)$, $\hat {\phi }_{d,j}(u;b)$ and $\hat {\phi }_{c,j}(u;b)$ denote the $j$th experiment values of $\hat {V}(u;b)$, $\hat {\phi }_d(u;b)$ and $\hat {\phi }_c(u;b)$, respectively.
In Table 2, we list the empirical average absolute errors of the estimator for the three functions. For the observation interval $[0, T]$, we take $T=1000\times 2^q$ for $q=0$, 1, 2, 3, 4, 5. It is obvious that the estimation errors are also decreasing with $q$ (or equivalently the time $T$), that is, the approximation results become closer to reference values, which means that our estimator performs better when the sample sizes become large.
In order to show the stability of our method, we plot $200$ estimated curves (green curves) and the true curves (bold blue curves) on the same picture. For the expected present value of total dividends before ruin with the exponential claim size density function, we plot the curves of $\hat {V}_j(u;b)$ as functions of $u$ when $T=1000\times 2^1$, $1000\times 2^3$ and $1000\times 2^5$in Figure 1. It is obvious that all the lines increase with $u$ for all the three models, which is consistent with the actual situation that the expected present value of total dividends before ruin increases with the initial surplus. Besides, as $T$ increases, the estimator $\hat {V}(u;b)$ tends to be stable and converges to ${V}(u;b)$ .
The functions $\phi _d(u;b)$, $\phi _c(u;b)$ and $\phi (u;b)$ are also investigated for Erlang(2,2) claim size density function. In Figure 2, we find that $\phi _d(u;b)$ is a decreasing function of the initial surplus $u$, which means that ruin is more likely to happen caused by oscillation when $u$ is small. At the same time, we can also observe that as $T$ increases, the estimator $\hat {\phi }_d(u;b)$ tends to be stable and converges to $\phi _d(u;b)$. As for $\hat {\phi }_c(u;b)$ and $\hat {\phi }(u;b)$, we plot the corresponding curves in Figures 3 and 4. Again, we can observe that the estimation becomes better as $T$ becomes larger.
Finally, we investigate the effectiveness of our method when further using spectral filters to accelerate the decay rate of the COS coefficients. Recall that a real and symmetric function $\Gamma (s)$ is a filter of order $\upsilon$ if: (i) $\Gamma (0)=1$, $\Gamma ^{(l)}(0)=0$, $1\leq l\leq \upsilon -1$; (ii) $\Gamma (s)=0$ for $|s|\geq 1$; (iii) $\Gamma (s)\in C^{\upsilon -1}$, $s\in \mathbb {R}$, where in particular $\Gamma ^{(l)}(\pm 1)=0$ for $0\leq l\leq \upsilon -1$. For a function $f$, its filter-COS approximation is simply defined by
In the following experiments, we select an exponential filter $\Gamma (s)=\exp (-\nu s^{\upsilon })$ with $\upsilon =6$ and $\nu =-\log \epsilon$, where $\epsilon$ represents the machine epsilon so that $\Gamma (1)=\epsilon \approx 0$. We take $\epsilon =10^{-6}$. In Table 3, we list the absolute errors of the three functions by using the filter-COS method. Similarly, we can find that the calculation errors decrease with the truncation parameter $K$. Compared with Table 1, we can find that the errors calculated by the filter-COS method are much smaller than the errors obtained without filter. In particular, we observed that for the $\tilde {V}(u;b)$ and $\tilde {\phi }_c(u;b)$ functions, the errors are significantly improved with the addition of filtering. We also list the empirical errors of the estimator for the three functions by the filter-COS method in Table 4, and almost all errors obtained by using the filter-COS method are also reduced compared with Table 2. In summary, spectral filters can accelerate the convergence rate.
Acknowledgments
The authors would like to thank the anonymous referees for useful suggestions. The research of Zhimin Zhang was supported by the National Natural Science Foundation of China [grant numbers: 11871121, 12171405] and the Fundamental Research Funds for the Central Universities [grant numbers 2021CDSKXYJG012, 2020CDJSK02ZH03].