A Cox model for gradually disappearing events

Abstract

Innovations in medicine provide us longer and healthier life, leading lower mortality. Sooner rather than later, much greater longevity would be possible for us due to artificial intelligence advances in health care. Similarly, Advanced Driver Assistance Systems (ADAS) in highly automated vehicles may reduce or even eventually eliminate accidents by perceiving dangerous situations, which would minimize the number of accidents and lead to fewer loss claims for insurance companies. To model the survivor function capturing greater longevity as well as the number of claims reflecting less accidents in the long run, in this paper, we study a Cox process whose intensity process is piecewise-constant and decreasing. We derive its ultimate distributional properties, such as the Laplace transform of intensity integral process, the probability generating function of point process, their associated moments and cumulants, and the probability of no more claims for a given time point. In general, this simple model may be applicable in many other areas for modeling the evolution of gradually disappearing events, such as corporate defaults, dividend payments, trade arrivals, employment of a certain job type (e.g., typists) in the labor market, and release of particles. In particular, we discuss some potential applications to insurance.

1. Introduction

Technological advances, such as medical innovations, have delivered significant improvements in survival and quality of life for us. Many diseases and epidemics were fatal to human being throughout a long history, and recently become curable. A notable example is Malaria. Dr. Tu Youyou managed to develop drugs based on artemisinin from Chinese traditional herbal medicines, which have led to the survival and improved health of millions of people. In 2015, she got the Nobel Prize in Physiology or Medicine for her discoveries concerning a novel therapy against Malaria. In addition, the National Cancer Institute recently found sudden reductions in mortality rates for prostate cancer, which is likely due to effective treatments, screening methods for early diagnosis and public health programs [18]. Moreover, although the recent coronavirus (COVID-19) pandemic is still ongoing currently, it is possible that the deaths caused by the COVID-19 would be gradually disappearing in the long run with more effective government interventions, vaccine development, and further enhancement of public health systems.

Similarly, the Advanced Driver Assistance Systems (ADAS) in highly automated vehicles may reduce or even ultimately eliminate accidents by perceiving dangerous situations. The development of highly automated vehicles and the ADAS is ongoing even though it is a challenge. Highly automated vehicles will minimize the number of accidents leading to a much lower number of claims for insurance companies ultimately. In addition, we can observe that the annual fatality rate per 100 million miles that motor vehicle traveled in the U.S. provided by the National Highway Traffic Safety Administration (NHTSA) [33] has declined gradually from 1920 to 1930 as plotted in Figure 1. The 1920's average rate is 18.52, and the 2010's average rate is only 1.13, where we can see that it declines substantially in a long run. Of course, other positive factors, such as improved infrastructure developments and safety regulations, may all together contribute to this decline. The longevity of ourselves has significantly increased in the last decade. It is certain for us to live longer and healthier on average in the long run due to artificial intelligence (AI) advances in health care. AI applications in the automotive industry will also change the risk landscape of insurance industry. They will help avoid accidents and saving lives, and consequently insurers will have fewer claims.

Figure 1. Motor vehicle traffic fatality rates, 1921–2018 (NHTSA).

In the actuarial literature, there is a plethora of papers aiming at modeling improvements in life expectancy due to systematic factors, notably, the Lee–Carter model [31] and its extensions which are mortality projection models on time-series mortality. Beside, alternative continuous-time stochastic models for modeling mortality can be found in for example [5,12,26,27,32,36]. The aim of this paper is to develop a specific continuous-time model which is based on doubly stochastic Poisson processes or Cox processes [9]. It could be used as a model component for the long-term survivor function capturing the greater longevity within a general competing-risks framework similarly as the literature of survival analysis for a particular failure type of interest, see for example [19,22,29, Sect. 8]. More precisely, it is a Cox Process with piecewise-constant decreasing intensity, since its underlying intensity process is piecewise-constant, and it jumps downward with random sizes to lower (but still positive) levels at random times as illustrated in Figure 2. Apparently, the occurrence timing of each innovation is uncertain in nature, so jump times in the underlying intensity process is random in our model, which is different from a simpler inhomogeneous Poisson process with a piecewise-constant decreasing intensity function. More generally, our model may be applicable to modeling many events of other types which are gradually disappearing. Recently, alternative intensity functions for Cox processes have been specified, see [1,2,14,25,37] and just to name a few.

Figure 2. Piecewise-constant decreasing intensity process $\lambda _t$.

We first obtain the Laplace transform of ultimate intensity integral in an analytic form. For insurance modeling based on a Cox process in general, the integral of the claim intensity process is crucial, as it is linked to the probability generating function (PGF) of claim numbers. More broadly, the integral of an underlying stochastic process has been played an important role to many applications in both finance and insurance. For example, the integral of the interest rate process is used to price zero-coupon bonds [11,23,24,38], and the integral of the stock price process is also used for pricing Asian-type options [4,8,15,17,34,35]. The integral of the default intensity process is also required to obtain the default probability in reduced-form credit risk models [16,30]. The integrated hazard rate is required to derive the survivor function [13,26] in life insurance. In non-life insurance, Jarrow [28] proposed pricing formulas for catastrophe bonds, where the integral of the intensity process is used to derive the probability of no catastrophic event. Under certain assumptions, other key distributional properties, such as the PGF of point process, their associated moments and cumulants, and the probability of no more claims for a given time point, are derived in analytic forms (with some additional assumptions), which are all important for model applications. We then apply our results to calculate the survival probability in life insurance and reinsurance premium in non-life insurance.

This paper is structured as follows. In Section 2, we define our model of the Cox Process with a piecewise-constant decreasing intensity, and analyze its theoretical distributional properties. In Section 3, we apply our results to calculate the survival probability for life insurance and stop-loss reinsurance premium for non-life insurance, respectively. Section 4 makes a brief conclusion.

2. A Cox process with piecewise-constant decreasing intensity

In this section, we explain how to construct a Cox process with a piecewise-constant decreasing intensity for modeling gradually disappearing events (e.g., insurance claims of a certain type) in general, and then obtain its key distributional properties.

Let us first start by a brief review for the Cox process [9] in general. If $\lambda _{t}$ is the intensity process of a Cox point process $N_t$, it is well known that the PGF of $N_{t}$ is given by

(2.1) \begin{equation} \mathbb{E}\left[ \theta ^{N_t}\mid \lambda_0 \right] = \mathbb{E}\left[ \exp \left.\left( - (1-\theta) \int_0^{t} \lambda_s \,\mathrm{d}s \right) \right| \lambda_0 \right],\quad \theta \in [0,1], \end{equation}

which suggests that the problem of finding the distribution of point process $N_{t}$ is equivalent to the problem of finding the distribution of the integral of intensity $\lambda _t$. By convention, we denote the intensity integral process by

$$\Lambda_t :=\int_0^{t}\lambda_s\,\mathrm{d}s.$$

In the large family of Cox processes, we can consider various candidates of non-negative stochastic processes for $\lambda _t$, which provide us a great flexibility for modeling event arrivals in practice. For more details about Cox processes, see for example [3,13], and the books by [6,7,10,20,21].

To construct a Cox process with piecewise-constant decreasing intensity, we introduce a positive continuous-time stochastic process $X_t$ as a state process for intensity process $\lambda _t$. That is, $\lambda _t$ is a deterministic function for a given state of $X_t$. We denote this function by $h(\cdot )$, that is, $\lambda _t=h(X_t)$, and further assume that it is a strictly increasing function of $u$ on the positive real line. For example, $h(u)$ could be a nonlinear function, such as

(2.2) \begin{equation} h(u)=u^{1+c},\quad u\in \mathbb{R}^{+},\enspace c> 0, \end{equation}

where the constant $c$ can be considered as a measure of the nonlinear amplification effect of states to the intensity, for example, the impact of vaccine advance to the fatality rate of COVID-19. In addition, we assume that state process $X_t$ is piecewise-constant and stochastically decreasing. More precisely, $X_t$ is constant until a random downward jump (i.e., random drop) potentially occurs at time point $T_i$ (for any given $i=1,2,\ldots$): it jumps to a new (lower) level $Y_i$ if $Y_i< X_{T_{i_{-}}}$, or, it just stays (the same) at $X_{T_{i_{-}}}$ if $Y_i>X_{T_{i_{-}}}$, where jump sizes $\{Y_i\}_{i=1,2,\ldots }$ are independent identical distributed (i.i.d.) with the cumulative distribution function (CDF) $G(y)$ and density function $g(y),y \in \mathbb {R}^{+}$. The interarrival times of these jumps at $\{T_i\}_{i=1,2,\ldots }$ are also i.i.d. with the density function $p(t)$, and they are independent of jump sizes $\{Y_i\}_{i=1,2,...}$. Through the functional transformation $h(\cdot )$, the resulting intensity process $\lambda _t=h(X_t)$, therefore, is still piecewise-constant and decreasing as visualized in Figure 2. Note that these arrival times $\{T_i\}_{i=1,2,\ldots }$ are labeled for the jump times in the intensity process $\lambda _t$ rather than the jump times generated from the point process $N_t$. The evolution of technological developments and breakthroughs, such as AI algorithms which have been developed for health care and automotive industry, could be aggregately modeled by this state process $X_t$ as a proxy. The impacts of these technological advances to the health care and automotive industries are measured by the nonlinear function of $h(u)$. The intensity process $\lambda _t$ remains constant until a breakthrough, represented by a drop to a new level as the result of this breakthrough.

2.1. Laplace transform of ultimate intensity integral

For convenience, the Laplace transform of intensity integral at the ultimate time $t\rightarrow \infty$ conditional on $X_0=x>0$ is denoted by

$$\phi(x) :=\mathbb{E}\left[ e^{{-}v \Lambda_{\infty }}\mid X_0=x \right],\quad v \geq 0.$$

Proposition 2.1. The Laplace transform of ultimate intensity integral $\Lambda _{\infty }$ conditional on $X_0=x>0$ satisfies

(2.3) \begin{equation} \phi(x)=\hat{p}\left( v h(x) \right) \left[\int_{0}^{x}\phi(y)g(y)\,\mathrm{d}y + \phi(x) \bar{G}(x) \right], \end{equation}

where

$$\hat{p}(v) := \int_{0}^{\infty }e^{{-}v u} p (u)\,\mathrm{d}u,\quad \bar{G}(u) := \int_{u}^{\infty } g(y)\,\mathrm{d}y.$$

Proof. By construction, the interarrival times of these jumps at $\{T_i\}_{i=1,2,\ldots }$ are independent of jump sizes $\{Y_i\}_{i=1,2,\ldots }$, therefore, $T_1$ and $X_{T_1}$ are independent, and we have

\begin{align*} \phi(x) & = \mathbb{E}\left[ \exp\left({-}v \int_0^{T_1}h(X_s)\,\mathrm{d}s \right) \times \exp\left({-}v \left.\int_{T_1}^{\infty}h(X_s)\,\mathrm{d}s \right) \right| X_0=x \right] \\ & = \mathbb{E}\left[ e^{{-}v h(x) T_1} \times \exp\left({-}v \left.\int_{T_1}^{\infty}h(X_s)\,\mathrm{d}s \right)\right| X_0=x \right] \\ & = \mathbb{E}\left[ \mathbb{E}\left[ e^{{-}v h(x) T_1} \times \exp\left({-}v \left.\int_{T_1}^{\infty}h(X_s)\,\mathrm{d}s \right) \right| X_0=x, T_1 \right] \right]\\ & = \mathbb{E}\left[ e^{{-}v h(x) T_1} \times \mathbb{E} \left[ \exp\left({-}v \left.\int_{T_1}^{\infty}h(X_s)\,\mathrm{d}s \right) \right| X_0=x, T_1 \right] \right] \\ & = \mathbb{E}\left[ e^{{-}v h(x) T_1} \times \phi(X_{T_1}) \right] \\ & = \hat{p}\left(v h(x)\right) \times \mathbb{E}\left[ \phi(X_{T_1}) \right]. \end{align*}

Note that, given the initial state level $X_0=x>0$ and the realization of the first jump size $Y_1=y$, the state process at the first-jump time, $X_{T_1}$, stays at the same level as $X_{T_1^{-}}$ if $y>X_{T_1^{-}}$, or moves down to a new level $y$ if $y< X_{T_1^{-}}$, so we have

$$\mathbb{E}\left[ \phi(X_{T_1}) \right]= \int_{0}^{x}\phi(y)g(y)\,\mathrm{d}y + \int_x^{\infty }\phi(x)g(y)\,\mathrm{d}y = \int_{0}^{x}\phi(y)g(y)\,\mathrm{d}y + \phi(x) \bar{G}(x).$$

Corollary 2.1. If jump sizes follow a uniform distribution on $[0,1]$, then, the Laplace transform of $\Lambda _{\infty }$ conditional on $X_0=x \in (0,1)$ is given by

(2.4) \begin{equation} \phi(x)=\hat{p}\left( v h(x) \right) \left[\int_{0}^{x}\phi(y)\,\mathrm{d}y + \left( 1-x \right) \phi(x) \right]. \end{equation}

Proof. Set $g(y) =1$ in (2.3) with the condition $X_0=x \in (0,1)$, the result follows immediately.

Corollary 2.2. If jump sizes follow a uniform distribution on $[0,1]$ and interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p(t) =\nu \,e^{-\nu t}$, then the Laplace transform of $\Lambda _{\infty }$ conditional on $X_0=x \in (0,1)$ is given by

(2.5) \begin{equation} \phi(x)=\exp \left( - \int_{0}^{x}\frac{ v h^{\prime}(y)}{\nu y+v h(y)}\mathrm{d}y \right). \end{equation}

Proof. Set $p(t) =\nu e^{-\nu t}$ in (2.4), then we have

$$\phi(x)=\frac{\nu }{\nu +v h(x)} \left[\int_{0}^{x}\phi(y)\,\mathrm{d}y+\left( 1-x \right) \phi(x)\right],$$

or

$$\nu \int_{0}^{x}\phi(y)\,\mathrm{d}y=\left[ \nu x +v h(x)\right] \phi(x).$$

Differentiate it with respect to $x$ both sides, then we have

$$\left[ \nu x +v h(x) \right] \phi^{\prime}(x)+v h^{\prime}(x)\phi(x)=0,$$

or

$$\phi^{\prime}(x)={-}\frac{v h^{\prime }(x)}{\nu x +v h(x)}\phi(x).$$

Given the initial condition $\phi (0)=1$, we can solve this ODE by

$$\ln \phi(x) - \ln\phi(0)={-}\int_{0}^{x}\frac{v h^{\prime}(y)}{\nu y+v h(y)}\mathrm{d}y,$$

and the result follows.

By further specifying the functional form for $h(x)$, then, nicely, we can find the exact distribution of $\Lambda _{\infty }$.

Corollary 2.3. If jump sizes follow a uniform distribution on $[0,1]$, interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p\left ( t\right ) =\nu \,e^{-\nu t}$ and $h(u)=u^{1+c},\ c> 0$, then the Laplace transform of $\Lambda _{\infty }$ conditional on $X_0=x \in (0,1)$ is given by

(2.6) \begin{equation} \phi(x)=\left(\frac{\frac{\nu }{x^{c}}}{\frac{\nu }{x^{c}}+v }\right)^{{(c+1)}/{c}}, \end{equation}

which implies

(2.7) \begin{equation} \Lambda_{\infty}\mid X_0=x \sim\textsf{Gamma}\left(\frac{c+1}{c},\frac{\nu }{x^{c}}\right). \end{equation}

Proof. Set $h(u)=u^{1+c},\ c>0$, then we have $h^{\prime }(u)= (c+1)u^{c}$. Hence, from (2.5), we have

\begin{align*} \phi(x) & = \exp \left( - (c+1) \int_{0}^{x}\frac{ v y^{c}}{\nu y+v y^{1+c}}\mathrm{d}y\right) \\ & = \exp \left( -\frac{c+1}{c} \int_{0}^{x}\frac{cy^{c-1}}{\frac{\nu }{v }+y^{c}}\mathrm{d}y \right) \\ & = \exp \left( -\frac{c+1}{c} \int_{0}^{x^{c} }\frac{1}{\frac{\nu }{v }+u}\mathrm{d}u \right) \\ & = \exp\left( - \frac{c+1}{c}\ln \left( \frac{\frac{\nu }{v }+x^{c}}{\frac{\nu }{v }} \right) \right)\\ & = \left( \frac{\frac{\nu }{x^{c}}}{\frac{\nu }{x^{c}}+v }\right)^{{(c+1)}/{c}}, \end{align*}

which is the Laplace transform of a gamma distribution with the shape parameter ${(c+1)}/{c}$ and rate parameter ${\nu }/{x^{c}}$.

2.2. Probability generating function of event numbers

Let us find the expression for the PGF of $N_\infty$.

Corollary 2.4. If jump sizes follow a uniform distribution on $[0,1]$ and interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p\left ( t\right ) =\nu e^{-\nu t}$, then the PGF of $N_{\infty }$ conditional on $X_0=x \in (0,1)$ is given by

(2.8) \begin{equation} \mathbb{E}\left[ \theta ^{N_{\infty }}\mid X_0=x \right] = \exp \left( - \int_{0}^{x}\frac{( 1-\theta )h^{\prime }(y)}{\nu y+\left( 1-\theta \right) h(y)}\mathrm{d}y\right). \end{equation}

Proof. Based on the relationship between the PGF of $N_{\infty }$ and the Laplace transform of $\Lambda _{\infty }$ as given by in (2.1), we can set $v=1-\theta$ in (2.5) and hence obtain (2.8) immediately.

Corollary 2.5. If jump sizes follow a uniform distribution on $[0,1]$, interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p\left ( t\right ) =\nu e^{-\nu t}$ and $h(u)=u^{1+c}, c>0$, then the PGF of $N_{\infty }$ conditional on $X_0=x \in (0,1)$ is given by

(2.9) \begin{equation} \mathbb{E}\left[ \theta ^{N_{\infty }}\mid X_0=x \right]= \left( \frac{1- \frac{x^{c}}{\nu +x^{c}} }{1- \frac{x^{c}}{\nu +x^{c}}\theta }\right)^{{(c+1)}/{c}}, \end{equation}

which implies

(2.10) \begin{equation} N_{\infty }\mid X_0=x \sim \textsf{Negative Binomial}\left( \frac{c+1}{c}, \frac{x^{c}}{\nu +x^{c}} \right). \end{equation}

Proof. Set $v =1-\theta$ in (2.8), then we have

$$\mathbb{E}\left[ \theta ^{N_{\infty }}\mid X_0=x \right]= \left( \frac{\frac{\nu }{x^{c}}}{\frac{\nu }{x^{c}}+ 1-\theta }\right)^{{(c+1)}/{c}}= \left( \frac{1- \frac{x^{c}}{\nu +x^{c}} }{1- \frac{x^{c}}{\nu +x^{c}}\theta }\right)^{{(c+1)}/{c}},$$

which is the PGF of a negative binomial distribution with parameters ${(c+1)}/{c}$ and ${x^{c}}/{(\nu +x^{c})}$.

2.3. Cumulants and moments of intensity integral

Let us find the expressions for the cumulants and the first two moments of $\Lambda _{\infty }$.

Theorem 2.1. The $m$th cumulant of $\Lambda _{\infty }$ conditional on $X_0=x \in (0,1)$ is given by

(2.11) \begin{equation} \kappa_m = m! \int_{0}^{x} \frac{h^{\prime }(y)} {\nu y} \left( \frac{h(y)}{\nu y} \right)^{m-1} \mathrm{d}y,\quad m=1,2,\ldots, \end{equation}

Proof. By taking logarithm for the Laplace transform (2.5), we obtain the cumulant generating function

\begin{align*} \ln \phi(x) & = \ln \mathbb{E}\left[e^{{-}v \Lambda_{\infty }}\mid X_0=x \right] \nonumber \\ & ={-} \int_{0}^{x}\frac{ v h^{\prime}(y)}{\nu y+ v h(y)}\mathrm{d}y \nonumber \\ & ={-}\int_{0}^{x} \frac{ \frac{h^{\prime }(y)} {\nu y} v }{ 1+ \frac{h(y)}{\nu y} v} \mathrm{d}y \nonumber\\ & ={-}\int_{0}^{x} \frac{h^{\prime }(y)} {\nu y} v \sum_{k=0}^{\infty} \left( -\frac{h(y)}{\nu y} v \right)^{k} \mathrm{d}y \nonumber\\ & = \sum_{k=0}^{\infty} \left[ \int_{0}^{x} \frac{h^{\prime }(y)} {\nu y} ({-}1)^{k+1} \left( \frac{h(y)}{\nu y} \right)^{k} \mathrm{d}y \right] v^{k+1}\quad m=k+1 \nonumber\\ & = \sum_{m=1}^{\infty} \left[ \int_{0}^{x} \frac{h^{\prime }(y)} {\nu y} ({-}1)^{m} \left(\frac{h(y)}{\nu y} \right)^{m-1} \mathrm{d}y \right] v^{m} \nonumber\\ & = \sum_{m=1}^{\infty} ({-}1)^{m} \frac{\kappa_m}{m!} v^{m}. \end{align*}

Hence, we have

$$\int_{0}^{x} \frac{h^{\prime }(y)} {\nu y} ({-}1)^{m} \left( \frac{h(y)}{\nu y} \right)^{m-1} \mathrm{d}y = ({-}1)^{m}\frac{\kappa_m}{m!},$$

and the cumulants (2.11) for any $m=1,2,\ldots$.

The expression for any moment of $\Lambda _{\infty }$ given $X_0=x$ can easily be obtained using its cumulants (2.11). For example, we provide the mean and variance as below.

Corollary 2.6. The mean and variance of $\Lambda _{\infty }$ conditional on $X_0=x \in (0,1)$ are, respectively, given by

(2.12) \begin{align} \mathbb{E}\left[ \Lambda_{\infty }\mid X_0=x \right] & = \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y, \end{align}

(2.13) \begin{align} \mathrm{Var}\left[ \Lambda_{\infty }\mid X_0=x \right] & = 2\int_{0}^{x} \frac{h^{\prime }(y) h(y)}{\nu ^{2}y^{2}} \mathrm{d}y. \end{align}

By specifying the form of $h(u),$ we may obtain the mean and variance explicitly.

Corollary 2.7. If $h(u)=u^{1+c},\ c>0$, then the mean and variance of $\Lambda _{\infty }$ conditional on $X_0=x \in (0,1)$ are respectively given by

\begin{align*} \mathbb{E}\left[ \Lambda_{\infty }\mid X_0=x \right] & = \frac{1}{\nu} \left( \frac{c+1}{c}\right) x^{c}, \\ \mathrm{Var}\left[ \Lambda_{\infty }\mid X_0=x \right] & = \frac{1}{\nu ^{2}}\left( \frac{c+1}{c}\right) x^{2c}. \end{align*}

Proof. The results follow immediately from (2.7) in Corollary 2.3. Alternatively, set $h(u)=u^{1+c}$ and $h^{\prime }(u)= (c+1)u^{c}$ in (2.12) and (2.13), respectively, then the results follow.

2.4. Moments of event numbers

Let us find the expressions for the first two moments of $N_{\infty }$.

Corollary 2.8. The mean and variance of $N_{\infty }$ conditional on $X_0=x \in (0,1)$ are respectively given by

\begin{align*} \mathbb{E}\left[ N_{\infty }\mid X_0=x \right] & = \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y, \\ \mathrm{Var}\left[ N_{\infty }\mid X_0=x \right] & = 2\int_{0}^{x} \frac{h^{\prime }(y) h(y)}{\nu ^{2}y^{2}} \mathrm{d}y + \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y. \end{align*}

Proof. Since

$$\mathbb{E}\left[ N_{\infty} \mid X_0=x \right] = \left.\frac{\mathrm{d}}{\mathrm{d} \theta }\mathbb{E} \left[\theta^{N_{\infty}}\mid X_0=x \right]\right|_{\theta =1}= \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y,$$

and

\begin{align*} \mathbb{E}\left[ N_{\infty}\left( N_{\infty}-1\right) \mid X_0=x \right] & = \left.\frac{\mathrm{d}^{2}}{\mathrm{d}\theta^{2}}\mathbb{E} \left[\theta^{N_{\infty}} \mid X_0=x \right]\right|_{\theta =1} \\ & =\left(\int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y \right)^{2} + 2\int_{0}^{x} \frac{h^{\prime }(y) h(y)}{\nu ^{2}y^{2}}\mathrm{d}y, \end{align*}

the variance of $N_{\infty }$ can be obtained by

\begin{align*} \mathrm{Var}[N_{\infty} \mid X_0=x] & = \mathbb{E}\left[ N_{\infty}\left( N_{\infty}-1\right) \mid X_0=x \right] + \mathbb{E}\left[ N_{\infty} \mid X_0=x \right] - \left(\mathbb{E}\left[ N_{\infty} \mid X_0=x\right] \right)^{2} \\ & = \left( \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y \right)^{2} + 2\int_{0}^{x} \frac{h^{\prime }(y) h(y)}{\nu ^{2}y^{2}} \mathrm{d}y + \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y - \left( \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y \right)^{2} \\ & = 2\int_{0}^{x} \frac{h^{\prime }(y) h(y)}{\nu ^{2}y^{2}} \mathrm{d}y + \int_{0}^{x} \frac{h^{\prime }(y)}{\nu y} \mathrm{d}y. \end{align*}

Corollary 2.9. If $h(u)=u^{1+c}, c>0$, then the mean and variance of $N_{\infty }$ conditional on $X_0=x \in (0,1)$ are respectively given by

\begin{align*} \mathbb{E}\left[ N_{\infty }\mid X_0=x \right] & = \frac{1}{\nu}\left( \frac{c+1}{c}\right) x^{c}, \\ \mathrm{Var}\left[ N_{\infty }\mid X_0=x \right] & = \frac{1}{\nu ^{2}}\left( \frac{c+1}{c}\right) x^{2c} + \frac{1}{\nu }\left( \frac{c+1}{c}\right) x^{c}. \end{align*}

Proof. The results follow immediately from (2.10) in Corollary 2.5.

For example, by setting the initial state $x=0.9$ and varying the parameter $c$ for amplification effect, the conditional means and variances of $N_\infty$ for $\nu =0.5,1$ are plotted in Figure 3.

Figure 3. Conditional mean and variance of $N_\infty$ against the parameter $c$ for amplification effect.

2.5. Probability of last event

As the intensity process is decreasing and eventually approaching to zero, the resulting points (i.e., events) are gradually disappearing, and an interesting problem is the probability of no more event (e.g., insurance claim of a certain type) beyond a given time point $t$, that is,

(2.14) \begin{equation} \Pr \left\{ T^{*} < t \right\} = \Pr \left\{ N_{[t,\infty)} =0 \right\} = \mathbb{E}\left[ \exp \left(-\int_{t}^{\infty }h(X_{s})\,\mathrm{d}s\right)\right], \end{equation}

where $T^{*}$ is the time point at which the last event occurs. In order to find (2.14), let us first derive the CDF of $X_t$ in Lemma 2.1.

Lemma 2.1. If interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p\left ( t\right ) =\nu \,e^{-\nu t}$, then the CDF of $X_{t}$ conditional on $X_0=x>0$ is given by

(2.15) \begin{equation} \Pr \left\{ X_{t}<\varsigma \right\} = 1-e^{-\nu t G(\varsigma )} \textbf{1}_{\left\{ \varsigma < x \right\}}. \end{equation}

Proof. By construction, $X_{t}$ is piecewise-constant and decreasing, which can be expressed by

$$X_{t}=\min\left\{ x,Y_{1},Y_{2},\ldots, Y_{M_{t}}\right\},$$

where $M_{t}$ is a homogeneous Poisson process of constant rate $\nu$. Given $M_{t}=m$, the CDF of the minimum of $m$ i.i.d. random variables $\{Y_i\}_{i=1,\ldots,m}$ with the CDF $G(\varsigma )$ is $\left [1-\left (1-G(\varsigma )\right )^{m} \right ]$, that is,

$$\Pr \left\{\min\left\{Y_{1},Y_{2},\ldots, Y_{M_{t}}\right\} < \varsigma \mid M_{t}=m\right\} = 1-\left( 1-G(\varsigma)\right)^{m}.$$

Then, we have

\begin{align*} \Pr \left\{\min\left\{Y_{1},Y_{2},\ldots, Y_{M_{t}}\right\} <\varsigma\right\} & = \sum_{m=0}^{\infty}\left[ 1-\left( 1-G\left( \varsigma \right) \right)^{m} \right] \frac{e^{-\nu t}\left( \nu t\right) ^{m}}{m!} \\ & = 1 - \sum_{m=0}^{\infty}\left( 1-G\left( \varsigma \right) \right)^{m} \frac{e^{-\nu t}\left( \nu t\right) ^{m}}{m!}, \end{align*}

where

$$\sum_{m=0}^{\infty}\left(1-G\left( \varsigma \right) \right)^{m} \frac{e^{-\nu t}\left( \nu t\right) ^{m}}{m!} = e^{ -\nu t G\left( \varsigma \right) }.$$

Note that,

\begin{align*} \Pr \left\{ X_{t}<\varsigma\right\} & = \Pr \left\{\min \left\{\min\left\{ Y_{1},Y_{2},\ldots, Y_{M_{t}}\right\}, x\right\} <\varsigma \right\} \\ & = 1- \Pr \left\{ \min \left\{ \min\left\{ Y_{1},Y_{2},\ldots, Y_{M_{t}}\right\}, x\right\} >\varsigma \right\} \\ & = 1- \Pr \left\{ \min\left\{ Y_{1},Y_{2},\ldots, Y_{M_{t}}\right\} > \varsigma \right\} \times \textbf{1}_{\left\{ x > \varsigma \right\}} \\ & = 1- e^{-\nu t G(\varsigma )} \times \textbf{1}_{\left\{ \varsigma < x \right\}}, \end{align*}

and the result (2.11) follows.

Based on Lemma 2.1, let us now derive the PGF of $T^{*}$ in Theorem 2.2.

Theorem 2.2. If jump sizes follow a uniform distribution on $[0,1]$ and interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p\left ( t\right ) =\nu \,e^{-\nu t}$, then the CDF of $T^{*}$ conditional on $X_0=x \in (0,1)$ is given by

(2.16) \begin{equation} \Pr \left\{ T^{*} < t \right\}=\nu t\int_{0}^{x} \exp \left( -\int_{0}^{\varsigma } \frac{h^{\prime }(y)}{\nu y+h(y)}\mathrm{d}y\right) e^{-\nu t\varsigma}\,\mathrm{d}\varsigma. \end{equation}

Proof. Set $v =1$ in (2.5), then we have

$$\mathbb{E}\left[ \exp \left( -\int_{t}^{\infty }h(X_s)\,\mathrm{d}s\right) \mid X_{t} \right] = \exp \left(-\int_{0}^{X_{t}}\frac{h^{\prime }(y)}{\nu y+h(y)}\mathrm{d}y\right).$$

Hence,

$$\Pr \left\{ T^{*} < t \right\}=\mathbb{E}\left[ \exp \left(-\int_{t}^{\infty }h(X_{s})\,\mathrm{d}s\right)\right] = \mathbb{E}\left[ \exp \left( -\int_{0}^{X_{t}}\frac{h^{\prime }(y)}{\nu y+h(y)}\mathrm{d}y\right) \right].$$

Based on the CDF (2.11) for $X_t$, we have the density of $X_t$ as

$$\Pr \left\{ X_{t} \in \mathrm{d} \varsigma \right\} = \nu t \,e^{-\nu t G(\varsigma )}g(\varsigma ) \textbf{1}_{\left\{ \varsigma < x \right\}} \mathrm{d} \varsigma,$$

and then,

$$\Pr \left\{ T^{*} < t\right\}=\nu t \int_{0}^{x} \exp \left(-\int_{0}^{\varsigma } \frac{h^{\prime }(y)}{\nu y+h(y)}\mathrm{d}y\right) e^{-\nu t G(\varsigma )}g(\varsigma)\,\mathrm{d}\varsigma.$$

With a uniform distribution on the interval $[0,1]$ for jump sizes, (2.12) follows immediately.

Corollary 2.10. If jump sizes follow a uniform distribution on $[0,1]$, interarrival times of jumps follow an exponential distribution with parameter $\nu >0$, that is, $p\left ( t\right ) =\nu \,e^{-\nu t}$ and $h(u)=u^{1+c},\ c>0$, then the CDF of $T^{*}$ conditional on $X_0=x \in (0,1)$ is given by

(2.17) \begin{equation} \Pr \left\{ T^{*} < t \right\} =\nu^{2+\frac{1}{c}} t \int_{0}^{x} \frac{1}{e^{\nu t \varsigma } \left( \nu + \varsigma^{c}\right)^{{(c+1)}/{c}}} \mathrm{d}\varsigma. \end{equation}

Proof. Based on Theorem 2.2 and the additional assumption of $h(u)=u^{1+c},\ c>0$, we have

\begin{align*} \Pr \left\{ T^{*} < t \right\} & = \nu t \int_{0}^{x} \exp \left( -\frac{c+1}{c}\int_{0}^{\varsigma } \frac{c y^{c-1}}{\nu +y^{c}}\mathrm{d}y\right) e^{-\nu t \varsigma}\,\mathrm{d}\varsigma \\ & = \nu t \int_{0}^{x} \exp \left( -\frac{c+1}{c} \ln\left( \frac{\nu + \varsigma^{c}}{\nu} \right) \right) e^{-\nu t \varsigma }\,\mathrm{d}\varsigma \\ & = \nu t \int_{0}^{x} \left( \frac{\nu}{\nu + \varsigma^{c}} \right)^{{(c+1)}/{c}} e^{-\nu t \varsigma}\,\mathrm{d}\varsigma\\ & = \nu^{2+\frac{1}{c}} t \int_{0}^{x} \frac{1}{e^{\nu t \varsigma } \left( \nu + \varsigma^{c}\right)^{{(c+1)}/{c}} } \mathrm{d}\varsigma. \end{align*}

We can calculate the probability of no event (e.g., insurance claim of a certain type) beyond time $t$ as given in (2.17) via numerical integration. For example, by setting parameters $\nu =2,\ c=1.2$ with the initial state $x=0.9$, the associated probabilities are reported in Table 1 and plotted in Figure 4. As time $t$ is getting larger, the intensity level is getting smaller and smaller. The intensity integral increases only marginally and events are still generated, but it will eventually arrive at the time point of no event anymore. Hence, the larger the time $t$ is, the higher the probability of no more event beyond time $t$ is.

Table 1. Probability of no event beyond time $t$.

Figure 4. Probability of no event beyond time $t$.

3. Applications in insurance

In this section, we discuss some potential applications of our model as a component to calculate the survival probability in life insurance and reinsurance premium in non-life insurance, respectively.

3.1. Life insurance

Medical innovation has made us live longer, healthier, and more prosperous lives. Additionally, due to artificial intelligence advances in health care, a much greater longevity would be possible in the long term as $t \rightarrow \infty$ or even within a couple of decades. Hence, we may apply our results to calculate survival probabilities which are the key inputs for life insurance. In general, by setting $v =1$ in (2.3) and solving the equation, we can obtain the ultimate survival probability (or survivor function) due to a particular event-type of interest as

(3.1) \begin{equation} \Pr \left\{ \tau^{*}_x = \infty \right\} = \mathbb{E}\left[ e^{-\Lambda_{\infty }}\mid X_0=x \right] = \hat{p}\left( h(x) \right) \times\left[ \int_{0}^{x}\phi(y)g(y)\,\mathrm{d}y + \phi(x) \bar{G}(x) \right], \end{equation}

where $\tau ^{*}$ is the first jump-arrival time in the point process $N_t$, that is,

$$\tau^{*}_x :=\inf\left\{ t>0: N_{t}=1\mid N_{0}=0,\ X_0=x \right\},$$

which is modeling the first event-arrival time of a particular type. More specifically, if jump sizes follow a uniform distribution on $[0,1]$, interarrival times of jumps follow an exponential distribution with its parameter $\nu$, that is, $p\left ( t\right ) =\nu \,e^{-\nu t}$ and $h(u)=u^{1+c},\ c>0$, then the survival probability $\Pr \left \{ \tau ^{*}_x = \infty \right \}$ can be immediately obtained by setting $v=1$ in (2.6). For example, by setting the initial state $x=0.9$ and varying the parameter $c$ for amplification effect, the survival probabilities for $\nu =0.5,1$, are plotted in Figure 5, and the associated detailed numerical output is reported in Table 2.

Figure 5. Survival probability against the parameter $c$ for amplification effect.

Table 2. Survival probability against the parameter $c$ for amplification effect.

Note that $\tau ^{*}_x$ is considered as a defective (improper) random variable with a point mass at time $t=\infty$. The survival probability (3.1) is corresponding to a particular type of gradually disappearing events (e.g., deaths caused by traffic accidents or deaths caused by the COVID-19) rather than overall events. Essentially, it is a marginal probability function, also well known as the subdistribution for a particular failure type of interest in the literature of survival analysis, see for example [19,22,29, Sect. 8]. The Cox process with piecewise-constant decreasing intensity as introduced in Section 2 provides us a key model component for modeling these types of gradually disappearing events rather than for all events within a general competing-risks framework. Of course, each life is not immortal in the reality, and we can simply add another model component acting as the force for all rest risks (e.g., the mortality force is increasing with respect to the age, year, or many other factors). This can be done by introducing another random time $\tau ^{*}_{b}$ acting as a competing risk that will eventually produce a death event at time

$$\tau^{*} := \min\left\{\tau^{*}_x,\tau^{*}_{b}\right\},$$

where $\tau ^{*}_{b}$ may be assumed to follow a simple Poisson process or more generally another Cox process with time-varying intensity $b(t)> 0$, and $b(t)$ could depends on years, ages, or many other factors. Therefore, the overall intensity is

$$\lambda_t + b(t),$$

where $b(t)$ is the baseline intensity for $\tau ^{*}$, and $\lambda _t$ is the type-specific intensity. Trivially, the ultimate overall survival probability (due to overall events) is zero, that is, $\Pr \left \{ \tau ^{*} = \infty \right \} = 0$. Therefore, here we are mainly interested in the nontrivial result of ultimate (marginal) survival probability (3.1) caused by events of a particular type, $\tau ^{*}_x$, driven by technological advances in a very long run. In fact, this research motivation is similar as the classical problem of ultimate ruin probability which are extensively studied in actuarial mathematics. In this paper, we are also dealing with the ultimate probability but within a different context considering the risk landscape change in the insurance industry.

3.2. Non-life insurance

The Advanced Driver Assistance Systems (ADAS) in highly automated vehicles may reduce or even eventually eliminate accidents by perceiving dangerous situations. Nearly perfect automated vehicles could be deployed as $t\rightarrow \infty$ or even within a couple of decades, which will minimize the number of accidents leading to fewer and fewer loss claims. Hence, we apply our results to the stop-loss reinsurance contract for a portfolio purely concentrated on the motor insurance business in the long run. Standard non-life insurance contracts are typically underwritten in the short term. However, they are often automatically renewed, and it is worth studying the long-term risk landscape in insurance similarly as the classical ruin problem in risk theory.

Let $\{Z_{i}\}_{i=1,2,\ldots }$ be the claim amounts, which are assumed to be i.i.d. with the CDF $H(z), z>0$. The total loss excess over the retention limit $K>0$ up to the time of $t\rightarrow \infty$ is given by

$$\left( C_{\infty }-K \right)^{+} :=\max\left\{ C_{\infty }-K,0\right\},$$

where $C_{\infty }=\sum _{i=1}^{N_{\infty }}Z_{i}$, and $N_{\infty }$ is the ultimate number of claims. Therefore, the stop-loss reinsurance premium is given by

$$\mathbb{E}\left[ \left( C_{\infty }-K \right)^{+} \mid X_0=x \right],$$

which is similar to a perpetual call option in finance.

Corollary 3.1. If jump sizes follow a uniform distribution on $[0,1]$, the interarrival times of jumps follow an exponential distribution with parameter $\nu$, that is, $p(t) =\nu \,e^{-\nu t}$ and $h(u)=u^{1+c},\ c>0$, and the claim sizes follow an Erlang distribution $\textsf {Erlang}(\varphi,\beta )$, then the expected ultimate total loss conditional on $X_0=x \in (0,1)$ is given by

$$\mathbb{E}\left[ C_{\infty} \mid X_0=x \right] = \frac{x^{c}}{\nu}\left( \frac{c+1}{c}\right) \frac{\varphi }{\beta },$$

and the stop-loss reinsurance premium with the retention limit $K$ is given by

(3.2) \begin{align} & \mathbb{E}\left[ \left( C_{\infty} - K\right)^{+} \mid X_0=x \right] \nonumber\\ & \quad =\left( \frac{x^{c}}{\nu +x^{c}} \right)^{{(c+1)}/{c}}\sum_{n=1}^{\infty } \frac{\Gamma(n+\frac{c+1}{c})}{n! \Gamma(n)}\left( \frac{\nu}{\nu +x^{c}} \right)^{n} \nonumber\\ & \qquad \times\left[\frac{n\varphi }{\beta } \left( 1-\frac{\gamma(n\varphi+1,\beta K)}{(n\varphi)!} \right) -K\left( 1-\frac{\gamma(n\varphi,\beta K)}{(n\varphi-1)!} \right)\right], \end{align}

where $\gamma (\cdot,\cdot )$ is the lower incomplete gamma function.

Proof. We assume that the claim sizes follow an Erlang distribution denoted by $\textsf {Erlang}(\varphi,\beta )$, that is,

$$h(z)=\frac{\beta^{\varphi }z^{\varphi -1}e^{-\beta z}}{(\varphi-1)!},\quad z,\beta >0,\enspace \varphi \geq 1,$$

where $\beta$ is the rate parameter and $\varphi$ is the shape parameter. Note that,

$$\sum_{i=1}^{N_{\infty}}Z_{i} \mid N_{\infty}\sim\textsf{Erlang}(N_{\infty}\varphi,\beta),$$

so, $C_{\infty }$ follows an Erlang-mixture distribution with the density function

$$\sum_{n=1}^{\infty }\Pr \left\{ N_{\infty }=n \mid X_0=x \right\}\times \frac{\beta ^{n\varphi }c^{n\varphi -1}e^{-\beta c}}{\left( n\varphi -1\right) !}.$$

Then, we have

$$\mathbb{E}\left[ C_{\infty} \mid X_0=x \right]= \mathbb{E}\left[ N_{\infty } \mid X_0=x \right] \times \mathbb{E}\left[ Z\right]= \frac{1}{\nu}\left( \frac{c+1}{c}\right) x^{c} \times\frac{\varphi }{\beta },$$

and

\begin{align*} & \mathbb{E}\left[ \left( C_{\infty} - K\right)^{+} \mid X_0=x \right] \\ & \quad =\mathbb{E}\left[ \left( C_{\infty }-K\right) \textbf{1}\{ C_{\infty }\geq K\} \mid X_0=x \right] \\ & \quad =\sum_{n=1}^{\infty }\Pr \left\{ N_{\infty }=n \mid X_0=x\right\} \times \left[\int_{K}^{\infty } c \frac{\beta ^{n\varphi }c^{n\varphi -1} e^{-\beta c}}{\left( n\varphi -1\right) !} \mathrm{d}c - K\int_{K}^{\infty } \frac{\beta ^{n\varphi }c^{n\varphi -1}e^{-\beta c}}{\left( n\varphi -1\right) !}\mathrm{d}c\right]\\ & \quad =\sum_{n=1}^{\infty }\Pr \left\{ N_{\infty }=n \mid X_0=x \right\} \times \left[\frac{n\varphi }{\beta }\int_{K}^{\infty } \frac{\beta ^{n\varphi+1}c^{n\varphi }e^{-\beta c}}{\left( n\varphi \right) !}\mathrm{d}c - K \int_{K}^{\infty }\frac{\beta ^{n\varphi }c^{n\varphi -1}e^{-\beta c}}{\left( n\varphi -1\right) !}\mathrm{d}c\right]\\ & \quad =\sum_{n=1}^{\infty }\Pr \left\{ N_{\infty }=n \mid X_0=x \right \} \times \left[\frac{n\varphi }{\beta } \left( 1-\frac{\gamma(n\varphi+1,\beta K)}{(n\varphi)!} \right) - K \left(1-\frac{\gamma(n\varphi,\beta K)}{(n\varphi-1)!} \right)\right] \\ & \quad =\sum_{n=1}^{\infty }\frac{\Gamma(n+\frac{c+1}{c})}{n! \Gamma(n)} \left(\frac{x^{c}}{\nu +x^{c}} \right)^{{(c+1)}/{c}}\left(\frac{\nu}{\nu +x^{c}} \right)^{n}\\ & \qquad \times\left[\frac{n\varphi }{\beta } \left( 1-\frac{\gamma(n\varphi+1,\beta K)}{(n\varphi)!} \right) - K \left(1-\frac{\gamma(n\varphi,\beta K)}{(n\varphi-1)!} \right)\right], \end{align*}

since

$$\int_{x}^{\infty} \frac{\beta ^{\varphi }z^{\varphi -1}e^{-\beta z}}{(\varphi-1)!} \mathrm{d}z = 1-\frac{\gamma(\varphi,\beta x)}{(\varphi-1)!},$$

and

$$N_{\infty} \mid X_0=x \sim \textsf{Negative Binomial}\left( \frac{c+1}{c}, \frac{x^{c}}{\nu +x^{c}} \right),$$

we have the probability mass function (PMF) of $N_{\infty }$ conditional on $X_0=x$ explicitly as

$$\Pr\left\{ N_{\infty }=n \mid X_0=x \right\}=\frac{\Gamma(n+\frac{c+1}{c})}{n! \Gamma(n)} \left(\frac{x^{c}}{\nu +x^{c}} \right)^{{(c+1)}/{c}} \left(\frac{\nu}{\nu +x^{c}} \right)^{n},\quad n=0,1,\ldots.$$

For example, by setting parameters $c=1.2,\varphi =1,\beta =1$ with the initial state $x=0.9$, we can calculate the expected total loss and stop-loss reinsurance premiums for $\nu =0.5,1$, which are plotted in Figures 6 and 7, respectively. The associated detailed numerical output is reported in Table 3.

Figure 6. Expected total loss against the parameter $c$ for amplification effect.

Figure 7. Stop-loss reinsurance premium.

Table 3. Stop-loss reinsurance premium.

4. Conclusion

A Cox process with piecewise-constant decreasing intensity has been introduced in this paper. The intensity process is a stochastic piecewise-constant decreasing function of time. Subject to the arrival of a breakthrough event, the intensity process might, or might not, be reduced to a new minimum level. The long-term properties are of our focus. We derive its ultimate distributional properties, such as the Laplace transform of the intensity integral process, the probability generating function of the point process, their associated moments and cumulants, and the probability of no more claims for a given time point. Using our model as a component within a general competing-risks framework, these results may be potentially applied to calculate survival probability for life insurance and reinsurance premium for non-life insurance in a very long run. Similarly, this model may be applicable in many other areas for modeling the evolution of gradually disappearing events, such as corporate defaults in credit risk modeling, trade arrivals in market microstructure, dividend payments by a stock, employment of a certain job type (e.g., typists) in labor market, and release of particles, as long as the underlying intensity process of the associated event arrivals is piecewise-constant and decreasing. Similarly as the classical ruin problem in insurance, the distributional properties and associated applications for a given finite-time horizon are much more difficult to be analytically studied, which are proposed for further research.