Bayesian Estimation of Renewal Function Based on Pareto-distributed Inter-arrival Times via an MCMC Algorithm

undefined

1. Introduction

In the literature, we often see the application of a single-parameter Pareto (SPP) distribution as a model to describe income and wealth. In insurance applications, heavy-tailed distributions of the SPP sort are essential tools for modeling extreme losses, especially for more risky types of insurance, such as medical malpractice insurance. In the field of actuarial science, the SPP distribution is combined with another distribution to model the portion of data characterized by large claims with low frequencies, i.e., high-severity risks. For the other portion of data characterized by small claims with high frequencies, i.e., low-severity risks, other distributions such as the exponential or the log-normal are used. Such models are called composite models. Deng and Aminzadeh (2020) developed a Bayesian predictive Weibull-Pareto composite model, and Aminzadeh and Deng (2019) developed a Bayesian predictive inverse gamma–Pareto composite model. Recent work in the field of network traffic has considered an SPP distribution as the interarrival distribution. Fras et al. (2010) considered the distribution for modeling self-similar network traffic. Harrison et al. (2014) used a network simulator to generate statistics on arrival times of traffic and concluded that an SPP distribution is more appropriate than the conventional exponential distribution in modeling interarrival times on a bursty network. The SPP distribution has also been used for insurance models: Hanafy (2020) used the distribution for interarrivals of claims on data from Egyptian fire insurance companies.

The Lomax distribution (also known as a Pareto Type II distribution) is a heavy-tailed probability distribution used in business, economics, actuarial science, queueing theory, and Internet traffic modeling. It is essentially the SPP distribution, but its support begins at zero and it has a probability density function (PDF) given as

$$f(x)=\frac{\alpha}{m}\big(1+\frac{x}{m}\big)^{-(\alpha+1)}, \hskip .1in x>0, \ \ \alpha, m >0$$

The above PDF can also be written as

$$f(x)=\frac{\alpha}{m}\big(\frac{x}{m}\big)^{-(\alpha+1)}, \hskip .1in x>m, \ \ \alpha, m> 0.$$

which is the standard form of the PDF for the SPP and which is also used by Ramsay (2006). Garsva et al. (2014) considered a Pareto Type II distribution as a model for the packet interarrival time distribution in an academic computer network, and they concluded via a Kolmogorov-Smirnov test that it was an appropriate choice.

In this paper, we analyze flood and tornado insurance damage data sets (provided in the Emergency Events Database (EM-DAT); see appendix) in the United States from the year 2000 to the year 2020. We show that both data sets fit the SPP distribution well. The interarrivals between two insured damages are defined as the difference between the start date of an occurrence and the start date of the subsequent occurrence.

Our purpose in this paper is to use SPP as the interarrival distribution of a renewal process and provide the ML- and MCMC-based Bayes estimators of the renewal function. Martino (2018) provides a thorough review of MCMC methods and compares them with numerical methods. In the current paper, we use the Metropolis-Hastings algorithm—see Casella and Berger (2004)—for the MCMC method to estimate the renewal function.

In Section 2, we define the renewal function and give computational formulas for the renewal function based on Ramsay (2006). These formulas are used in the Mathematica code for simulations and estimation of the renewal function. Section 3 provides the ML estimation for the parameters of the SPP distribution. Section 4 provides the derivation of two conditional posterior distributions that are based on shifted-exponential and gamma priors. In Section 5 we outline the steps taken in the Mathematica code to use the MCMC algorithm and discuss a summary of the simulation results. Section 6 provides numerical examples to illustrate the computations, and Section 7 concludes. The Mathematica code used for the simulations and the estimation of the renewal function in this paper is available from the authors.

2. Renewal function of the SPP distribution

Let $\{N(t), t\geq 0\}$ be a counting process, where $N(t)$ denotes the number of events that occur by time $t.$ Let $X_n$ denote the time between the $(n-1)$st and $n$th events of this process, $n\geq1.$ If $\{X_1,X_2,...\}$ is a sequence of independent identically distributed (iid) random variables, then the process $\{N(t), t\geq 0\}$ is called a renewal process.

Renewal theory has many applications in stochastic modeling, and as a result, statistical inference with respect to the renewal function is of interest in many applied fields. For example, in reliability theory, a renewal process can be used to model successive repairs of a failed machine. In actuarial science, a renewal process can be used to model the successive occurrences of risks. Hardy (2006) provides many examples for risk measures with actuarial applications, such as in property and casualty (P&C) insurance, life insurance, and the banking industry. Necir, Rassoul, and Meraghni (2010) provide a semiparametric estimate of the renewal function for heavy-tailed risks. Woo (2011) provides two-sided bounds for the renewal function considering a variety of reliability classifications.

In this paper, we consider ML- and MCMC-based Bayes estimators of the renewal function. Actuaries should find the results useful in applications where events (risks) occur according to a counting process and the renewal random variable has an SPP distribution. One of the measures associated with a renewal process is the renewal function, $M(t)=E[N(t)].$ In some applications of renewal processes, one is interested in estimating the number of renewals by time $t.$ Provided that the parameters associated with the interarrival distribution are known or estimated, the computation can be accomplished by using the integral equation

$$\begin{align} M(t)&=E[N(t)]=\int_0^\infty E[N(t)|X_1=x]f(x)dx\\ &=F(t)+\int_0^\infty M(t-x)f(x)dx, \end{align} \tag{1}$$

where $f$ and $F$ are, respectively, the PDF and the cumulative distribution function (CDF) of the interarrival random variable $X.$ For example: (a) for the uniform(0,1) interarrival distribution, $M(t)=e^t-1, 0\leq t\leq 1,$ and (b) for the exponential interarrival distribution with the PDF $f(x)=\frac{1}{\lambda}e^{-\frac{x}{\lambda}},$ where $\lambda$ is the scale parameter, $M(t)=\lambda t, t > 0,$ the reason being that when the interarrival times have an exponential distribution, $N(t)\sim \text{Poisson}(\lambda t).$ However, the derivation of $M(t)$ in a closed form for many interarrival distributions is not straightforward. Let $S_n=\sum_{i=1}^nX_i$ denote the time of $n$th renewal. It can be shown that

$$M(t)=\sum_{n=1}^\infty F_{S_n}(t),$$

$$\partial M(t)/\partial t=M^{'}(t)=\sum_{n=1}^\infty f_{S_n}(t),\tag{2}$$

where $f_{S_n}$ and $F_{S_n},$ respectively, are the PDF and the CDF of $S_n.$ Therefore, to find $M(t)$ in a closed form, not only must we identify the distribution of $S_n,$ but we must also evaluate $\sum_{n=1}^\infty F_{S_n}(t),$ which may not be trivial.

Ross (2007) supplies us with a method to approximate the renewal function. The method computes $M_2(t), M_3(t),$ …, $M_n(t)$ recursively to approximate the renewal function $M(t)$ using

$$M_r(t)= \\{{\sum_{k=1}^{r-1}(1+M_{r-k}(t))E(X^ke^{-\delta_r X})(\delta_r^k/k!)+E(e^{-\delta_r X})}\over {1-E(e^{-\delta_r X})}},$$ $$M_1(t)={{E(e^{-\delta_1 X})}\over {1-E(e^{-\delta_1 X})}},$$

where $\delta_r=r/t$ is used for the computation of $M_r, r=1,2,....$ Therefore, $M_1(t), M_2(t), M_3(t),..., M_n(t)$ are found consecutively to approximate the renewal function $M(t).$ Ross indicates that for a large $r,$ $M(t)$ is approximated by $M_r(t)$ quite well. The method is useful when $E(X^ke^{-\delta_r X})$ and $E(e^{-\delta_r X})$ have explicit expressions for an interarrival distribution. For example, if the interarrival distribution is gamma (2,1) (shape parameter = 2, and scale parameter = 1), then

$$E[X^k e^{-\frac{r X}{t}}]=t^2(\frac{r+t}{t})^{-k}\frac{\Gamma(2+k)}{(r+t)^2}, \ \ \ \ \frac{r}{t}> -1,$$

Then using (2) it can be shown that the exact expression for the renewal function is

$$M(t)=.5t -.25(1-e^{-2t}), t>0.$$

For $t=10,$ $M(10)=4.75$ is the exact value, where its approximate value based on $r=50$ is calculated in Ross (2007) as 4.75. As another example, for the interarrival random variable with PDF $f(x)=1.4e^{-2x}+.3e^{-x}$ and CDF $F(x)=1-.3e^{-x}-.7^{-2x},$ $M(10)$ is approximated as 15.5089 (see Ross 2007). For some interarrival distributions, $M(t)$ cannot be derived in a closed form. For example, consider the inverse Gaussian (IG) distribution and let $X_i \sim \text{IG}(\phi,\lambda);$ it can be shown that $S_n\sim \text{IG}(\phi/n, n^2\lambda),$ but $\sum_{n=1}^\infty F_{S_n}(t)$ cannot be written in a closed form. See Aminzadeh (2011), where the author provides a Bayes estimator of the renewal function using IG distributed interarrivals. For $X \sim \text{IG}(\phi,\lambda),$ with given values for $r$ and $t$ and unknown values for $\phi$ and $\lambda,$ Mathematica cannot provide a closed form for

$$\int_0^\infty x^k e^{- (r/t ) x}f(x; \phi,\lambda)dx.$$

However, for example, with $\phi=2,\lambda=3, t=20,$ and $r=5,$ for the above expected value we get

$$2^k 3^{1/4+k/2}5^{1/4+k/2} e^{1.5} \text{BesselK}[.5-k,\sqrt{15}/2]/\sqrt{\pi}$$

via Mathematica with a long integration time, and the expression is quite complicated. SPP$(m,\alpha),$ which is considered in the current paper, is an even more difficult case, as

$$\int_m^\infty x^k e^{- (r/t ) x}f(x; m,\alpha)dx$$

cannot be obtained in a closed formula, even for selected values of $m,\alpha, t,$ and $r.$

For the SPP distribution, the approximation Ross provides is not practical. The reason is that $E(X^ke^{-\delta_r X})$ cannot be written in a closed form and therefore an approximate expression for the expected values must be derived or evaluated numerically, which would depend on $k.$ Even if a "good"expression for the expected value is found, different values of $k$ would produce different expressions/approximations for the expected value. In addition, the recursive formula for $M_r(t)$ uses previous approximations for the expected value. Therefore even for a small value of $r,$ $E(X^ke^{-\delta_r X})$ in the above sum has to be approximated for $k=1,2,...(r-1).$ In addition to the $r-1$ approximation errors, one has to take into account the accumulated approximation errors for the expected value that are involved in the computation of $M_{r-k}(t).$ Thus, due to the many accumulated approximation errors involved for approximating $E(X^ke^{-\delta_r X})$ and $M_r(t),$ Ross’s method is not considered in this paper.

The PDF of the SPP$(m,\alpha)$ random variable $X,$ with parameters $\alpha$ and $m,$ is given by

$$\begin{align} f_{X}(x;\alpha,m)&= \alpha x^{-(\alpha+1)}m^{\alpha},\\ &\quad x>m,\ \ \text{and}\ \ \\ &\quad \alpha, m>0. \end{align} \tag{3}$$

Consider a renewal process $\{N(t), t\geq 0\}$ and let $X_n,$ the time between the $(n-1)$st and $n$th renewal, be distributed as SPP$(\alpha, m),$ $n =1,2,3,....$ Suppose from this process a random sample $X_1,...,X_H$ is collected. The random sample will be used to find the ML and an approximate Bayes estimate for $M(t)$ via an MCMC algorithm.

Ramsay (2006) applied Laplace transforms and generalized exponential integrals to develop the following exact expressions for PDF $f_{S_n}$ and CDF $F_{S_n}$ of $S_n=\sum_{i=1}^n X_i$ under the assumption that $X_i, i=1,2,....$ are iid from the SPP distribution and that $\alpha$ is a positive integer. Recall that $S_n$ is the time of $n$th renewal. The following formulas are from Ramsay (2006):

$$f_{S_n}(s)=\frac{1}{nm}\int_0^\infty e^{-(1+\frac{s}{nm})v}\phi_{\alpha,n}(\frac{v}{n})dv,$$

where,

$$\begin{align} \phi_{\alpha,n}(v)&=(-1)^{n+1}\alpha^n \sum_{r=0}^{\lfloor (n-1)/2 \rfloor}(-\pi^2)^r \\ &\quad \binom{n} {2r+1}(G_{\alpha+1}(v))^{n-2r-1}\big(\frac{v^\alpha}{\alpha!}\big)^{2r+1}, \end{align}$$

$$\begin{align} G_\alpha(v)&=\frac{v^{\alpha-1}}{(\alpha-1)!}\big[\gamma +\ln(v)-\sum_{r=1}^{\alpha-1}\frac{1}{r}\big]\\ &\quad +\sum_{r=0, r\ne \alpha-1}^\infty \frac{v^r}{(r-\alpha+1)r!}, \\ \gamma&=.5772156649, \end{align}$$

and

$$\begin{align} F_{S_n}(t)&=\int_0^t f_{S_n}(s)ds \\ &=\int_0^\infty \frac{1}{v}(1-e^{-\frac{tv}{m}})e^{-nv}\phi_{\alpha,n}(v)dv. \end{align} \tag{4}$$

The approach for computing $F_{S_n}(t)$ in Ramsay is to use a series expansion and then use

$$\begin{align} e^{-nv}\phi_{\alpha,n}(v)&=(-1)^{n+1}\alpha^n\sum_{r=0}^{\lfloor (n-1)/2 \rfloor}(-\pi^2)^r \\ &\quad \binom{n} {2r+1}(e^{-v}G_{\alpha+1}(v))^{n-2r-1}\\ &\quad \big(\frac{e^{-v}v^\alpha}{\alpha!}\big)^{2r+1} \end{align}$$

in (4) to compute $F_{S_n}(t).$ The above formulas are built into the Mathematica code of the current paper, and the NIntegrate function of Mathematica is used to compute (4) based on input values $m,\alpha,n,$ and $t.$

Although the results given in Ramsay assume that $\alpha$ is a positive integer, simulation studies confirm that for both ML and Bayes estimates of the renewal function, once the ML and Bayes estimates of the parameters $\alpha$ and $m$ are obtained, the estimated values of $\alpha$ can be rounded to the nearest integer so that $F_{S_n}(t)$ and as a result an estimate for $M(t)$ can be obtained. The reason that in Ramsay $\alpha$ is assumed to be a positive integer is that to derive (4), for complex $z$ and Re$(z)>0,$ it is shown that the Laplace transform of $f_n(t)$ is given by $f^{*}_n(z)=(f^{*}(z))^n,$ where $f^{*}(z)=\alpha e^{mz}E_{\alpha+1}(mz),$ and $E_{\alpha}(z)=\int_1^\infty \frac{e^{-zt}}{t^{\alpha}} dt,$ for $\alpha=1,2,3,...,$ is the generalized exponential integral. And then $f_n^{*}(z)$ is inverted to get the Bromwich integral for $f_n(t),$ and consequently the CDF in (4) is derived. Simulations confirm that by rounding the estimated values of $\alpha,$ the average squared error (ASE) of the ML and Bayes estimates of the renewal function is still quite small and that the MCMC-based Bayesian method outperforms the ML method with regard to accuracy when the sample size $H$ is small. For large sample sizes, the ML and Bayes methods have comparable accuracies.

3. Maximum likelihood estimation of M(t)

Based on a random sample $X_{1},...,X_{H}$ from SPP$(\alpha, m),$ and using the likelihood function

$$\begin{align} L(\alpha, m)&= \alpha^H(\prod_{i=1}^H X_i)^{-(\alpha+1)}m^{H\alpha}, \\ &\quad X_i>m, \hskip.1in i=1,2,...,H, \end{align} \tag{5}$$

it can be shown that the ML estimates of $\alpha$ and $m$ are

$$\hat{m}=X_{(1)}, \hskip .1in \hat{\alpha}=\frac{H}{T-H \ln \big(X_{(1)}\big)},\hskip .1in T=\sum_{i=1}^H \ln(X_i),$$

and that $(X_{(1)}, T)$ is a minimal sufficient statistic for $(m, \alpha).$ $X_{(1)}$ is the first-order statistic. The goal is to estimate

$$M(t)=\int_0^t M^{'}(s)ds=\sum_{n=1}^ \infty F_{S_n}(t). \tag{6}$$

Therefore, for each value of $n$ in (6), and given estimates of $(\alpha,m)$ (based on a random sample $X_1,...,X_H),$ we can compute $F_{S_n}(t)$ via (4). Computations confirm that for selected values $\alpha, m,$ and $t,$ the CDF $F_n(t)$ in (4) is a decreasing function of $n.$ The Mathematica code in a while loop searches for an appropriate upper limit $k$ for the sum in (6). The loop starts with $k=1$ and continues until $F_k(t)-F_{k-1}(t)>0,$ and, using the value of $k,$ computes $\sum_{n=1}^{k-1} F_n(t)$ as an estimate for $M(t).$

It is also noted that through an estimate of $M(\tau),$ where $\tau$ is a truncated time, the excess life of a renewal at time $\tau$ can be estimated. Let $Y(\tau)$ denote the excess life at $\tau$ and $\mu=E[X]$ be the mean of interarrival random variable $X.$ We can use the well-known formula

$$M(\tau)=\frac{\tau}{\mu}-1+\frac{E[Y(\tau)]}{\mu},$$

to find the ML estimate of $E[Y(\tau)]$ via the ML estimate $\hat{\mu}=\frac{\hat{\alpha}\hat{m}}{\hat{\alpha}-1}$ and the ML estimate of $M(\tau).$ The same formula above can be used to get an approximate (using simple substitution) Bayes estimate of $E[Y(\tau)],$ once Bayes estimates of $\alpha, m,$ and $M(\tau)$ are found through the MCMC algorithm given in Section 5.

4. Bayesian estimation of M(t)

In the literature, we find many instances where authors consider Bayesian inference for the parameters of the SPP distribution. Sun et al. (2021) adopt the Jeffreys prior to obtain Bayes estimates for $m$ and $\alpha.$ Kim, Kang, and Lee (2009) also provide Bayes estimates via the non-informative Jeffreys prior for SPP. Noor and Aslam (2012) provide non-informative as well as informative Bayes estimates of SPP parameters; however, in their paper, gamma priors are used for both $m$ and $\alpha.$

The characteristics of the shifted exponential distribution with regard to its support $(a, \infty)$ are similar to the support of SPP. This is the main reason that in the current paper we use the shifted exponential prior as opposed to another distribution with support $(0,\infty).$ The PDF of the shifted exponential prior is given as

$$\rho(m|a,b)=\frac{1}{b}e^{-(\frac{m-a}{b})}, \ \ m > a, \ \ a, b >0, \tag{7}$$

where $a$ (location parameter) and $b$ (scale parameter) are hyperparameters associated with the prior distribution.

For the conditional prior of $\alpha|m,$ gamma$(c,d)$ is selected, with the PDF

$$\rho(\alpha|m)\propto \alpha^{c-1}e^{-\frac{\alpha}{d}},\ \ c ,d >0,\tag{8}$$

where $c$ (shape parameter) and $d$ (scale parameter) are hyperparameters associated with the conditional prior distribution. As mentioned earlier, although computation of the CDF $F_{S_n}(t)$ requires that $\alpha$ be an integer (as the formulas in (4) use only positive integers for $\alpha),$ using a gamma distribution for the conditional prior does not cause an issue with regard to the accuracy of the ML and Bayes estimates of the renewal function. For the current paper, in the Mathematica code for simulations, once a Bayes estimate for $\alpha$ is obtained, it is rounded to the nearest integer for computational purposes. From (7) and (8), we get the joint prior as

$$\rho(\alpha,m)\propto \alpha^{c-1}e^{-(\frac{\alpha}{d}+\frac{m}{b})},\ \ m>a.\tag{10}$$

Now, using the likelihood function for a sample of size $H,$ the joint posterior is given by

$$\pi(\alpha, m |\underline{x})\propto \alpha^{H+c-1}e^{-\alpha(\frac{1}{d}+T)}m^{H\alpha} e^{-\frac{m}{b}}, \ \ m>a. \tag{11}$$

To apply an MCMC algorithm, (11) is used to find the conditional posteriors as

$$\pi(\alpha |\underline{x}, m)\propto \alpha^{H+c-1}e^{-\alpha(\frac{1}{d}+T)}m^{H\alpha}, \tag{12}$$

and

$$\pi(m |\underline{x},\alpha)\propto e^{-\frac{m}{b}}m^{(H\alpha+1)-1}, \ \ m>a. \tag{13}$$

Simon and Adler (2021) provide a thorough reference on MCMC algorithms and compare several MCMC methods for estimating the parameters of the Pareto/negative binomial distribution (NBD) model and their performance. The Handbook of Markov Chain Monte Carlo (Brooks et al. 2011) covers MCMC foundations, methodology, and algorithms, as well as MCMC’s use in a variety of practical applications in areas such as educational research, astrophysics, brain imaging, ecology, and sociology. Many authors have discussed the MCMC method via conditional posterior distributions and the Metropolis algorithm. In particular, the description and implementation of the Metropolis algorithm can be found in Casella and Berger (2004), among other books and articles.

Ji, Aminzadeh, and Deng (2020) propose an MCMC Bayesian approach using the Metropolis algorithm and conditional posteriors to obtain estimates of the parameters as well as an approximate predictive density via simulation, which can be used to compute life expectancy and other measures of interest in a Bayesian framework. Aminzadeh and Deng (2021) employ an MCMC method based on the Metropolis algorithm and conditional posterior distributions to estimate ruin probability based on nonhomogeneous Poisson process claim arrivals and inverse Gaussian–distributed claim aggregates. The MCMC algorithm we use in this paper is based on the same logic as the aforementioned articles.

Considering the conditional posterior distributions in (12) and (13), we can see that the kernel of the PDF in (12) does not belong to the gamma$(H+c, \frac{d}{1+dT})$ distribution, due to the extra term $m^{H\alpha}.$ Therefore, we use gamma$(P_1, P_2)$ as a proposal distribution for the PDF in (12) within the MCMC algorithm. $P_1$ and $P_2$ are, respectively, the shape and scale parameters of the proposal distribution gamma$(P_1,P_2).$ For selected "true"values of $\alpha$ and $m,$ in simulation studies, values of $P_1$ and $P_2$ are selected in such a way that $P_1P_2$ is close enough to $\alpha$ while $P_2$ is very small. This is because the expected value and variance of the proposal gamma$(P_1,P_2)$ distribution are, respectively, $P_1P_2$ and $P_1P_2^2.$

As for the conditional posterior (13), we can see that it resembles the kernel of the generalized gamma$(g,\beta,\xi,\mu)$ distribution with the PDF

$$f(y|g,\beta,\xi,\mu)=\xi e^{-(\frac{y-\mu}{\beta})^\xi}(\frac{y-\mu}{\beta})^{g \xi-1}, y>\mu,$$

where $\beta=b, \mu=a, \xi=1,$ and $g=H\alpha+1.$ Note that $\mu$ is the location parameter and $\beta$ is the scale parameter.

5. Simulation

For selected values of $(\alpha, m, t),$ the Mathematica code (available from the authors upon request) first computes the "true"value of the renewal function $M(t)$ via (4). The code then generates random samples $X_1,...,X_H$ from the SPP$(\alpha,m)$ distribution using a built-in function in Mathematica. For the generated sample $i, (i=1,2,....,N=100),$ the following steps are taken:

1. ML estimates of $\alpha$ and $m$ are computed via the formulas in Section 3.

2. For each value of $n$ in the While loop, ML estimates from Step 1 are used to compute the CDF in (4). The while loop starts with $k=1$ and computes $F_{S_1}(t),F_{S_2 }(t), F_{S_3}(t),...$ via (4). Then for each value of $k (k=1,2,3...),$ $\sum_{n=1}^k F_{S_n}(t)$ is calculated. Computations confirm that for given values of $\alpha, m,$ and $t,$ the CDF $F_{S_n}(t)$ in (4) is a decreasing function of $n.$ The Mathematica code in the while loop searches for an appropriate upper limit $k$ for the sum in (6). Computations continue until $F_k(t)-F_{k-1}(t)>0.$ And, using the value of $k,$ the loop ends and ${M_{i,ML}(t)}= \sum_{n=1}^{k-1} F_{S_n}(t)$ is reported as the ML estimate of the renewal function based on sample $i.$ It is noted that in this step, the ML estimates of $\alpha$ and $m,$ which are based on a generated sample of size $H$ from the SPP, are used in (4) to find $F_{S_n}(t)$ for $n=1,2,3,...$ and that $S_n$ is a symbolic notation in the CDF $F_{S_n}(t).$ The computation of (4) does not depend on a value of $S_n.$

3. The mean and ASE of estimates $M_{i,ML}(t), (i=1,2,...,N)$ of the renewal function based on $N$ simulated samples are respectively computed as

$$\begin{align} \overline{\hat{M}}_{ML}(t)&=\frac{\sum_{i=1}^N M_{i,ML}(t)}{N}, \\ \xi_{ML}(\hat{M}(t))&=\frac{\sum_{i=1}^N(M_{i,ML}(t)-M(t))^2}{N}. \end{align}$$

Note that the term MSE (mean squared error) refers to the expected value of an estimator under the squared-error loss function. In this paper, the term ASE refers to $\sum_{i=1}^N (\hat{\theta}_i-\theta)^2/N,$ where $\hat{\theta}_i$ is an estimate of the "true"parameter $\theta$ based on a generated sample $i,$ $(i=1,2,...N).$

To compute the MCMC-based Bayes estimate of the renewal function, we use the conditional posterior distributions in (12) and (13). Let $f_1(\alpha |m_0)=\alpha^{H+c-1}e^{-\alpha(\frac{1}{d}+T)}m_0^{H\alpha}$ denote the conditional posterior PDF (RHS in (12)) of $\alpha$ given $m=m_0,$ and $h_1(\alpha)$ is the PDF of the proposal distribution gamma $(P_1,P_2)$ for $\alpha.$ As mentioned earlier, in simulation studies, for selected "true"values of $\alpha$ and $m,$ the values of $P_1$ and $P_2$ are selected in such a way that $P_1P_2$ is close enough to $\alpha$ while $P_2$ is very small. Also, the hyperparameters $c$ and $d$ are selected such that $c d\approx \alpha,$ while $d$ is small. Again, this is due to the expected value and variance of the gamma distribution. For example, if we select $\alpha=3$ in simulation, then one possibility would be $c = 300, d = .01, P_1 = 10,$ and $P_2 = .3.$

It is noted that the expected value and the variance of the generalized gamma$(H\alpha+1,b,1,a),$ which is used as the conditional posterior of $m|\alpha,\underline{x},$ are, respectively,

$$E[m|\alpha]=a+Hb\alpha, \hskip .1in \text{Var}(m|\alpha)=b^2(H\alpha+1),$$

Therefore, the hyperparameters $a$ and $b$ are selected so that $m \approx a+ Hb\alpha$. We choose a value for $a$ so that it is close to $m(a \lt m)$ and then find the corresponding value for $b=\frac{m-a}{H\alpha}$. This way the Bayes estimates of $m$ will be stable with a higher precision. For example, if the selected values are $m=2.2$, $H=30$, and $\alpha=3.0$, then one possibility is $a=2.1$ and $b=.0011$. In practice, given a sample from the SPP distribution, one does not have information on the "true"values of $\alpha$ and $m$. If $P_1$ and $P_2$ are selected such that their product is close to the ML estimate of $\alpha$, while $P_2$ is small, the Bayes estimates of $\alpha$ and $m$ would be stable with a small ASE, and as a result the renewal function can be estimated with a higher accuracy.

For example, if the ML estimates are $\hat{m}=5.15, \hat{\alpha}=2.15\approx 2.0=\alpha^{*},$ and $H=30,$ then we can choose hyperparameters as

$$\begin{align} c=500&, d=.001, P_1=10, \\ P_2=.5&, a=2, b=\frac{\hat{m}-a}{H\alpha^{*}}=.0525 \end{align}$$

The same approach is used in Aminzadeh and Deng (2018) and Deng and Aminzadeh (2020) for the selection of hyperparameter values. Using ML estimates through a random sample helps to identify appropriate values for hyperparameters that are consistent with the sample information. Selection of hyperparameters using ML estimates is considered in the simulation studies (see Section 5 for a summary of results) as well as in the numerical example in Section 6.

4. The MCMC algorithm in this step is based on the Metropolis algorithm (see Casella and Berger 2004). The idea is to use a proposal distribution that is close enough to a conditional posterior distribution such as (12), which is not a recognizable distribution, and implement the Metropolis algorithm through a very large number of iterations (we used $n_s=$ 20,000) so that the Bayes estimates converge.

To get the MCMC-based estimate of the renewal function, for each generated sample $i, (i=1,...., N)$ from SPP$(\alpha, m),$ the following steps are taken:

(a) Generate an initial value $\alpha_0$ from gamma$(P_1,P_2)$ and set the first entry of the array $\alpha [1]=\alpha_0.$ The method for choosing the shape parameter $P_1$ and the scale paremeter $P_2$ is outlined in the step 3.

(b) Generate a random observation $W_1$ from the generalized gamma $(H\alpha_0+1, b, 1, a),$ set $m_0=W_1,$ and let the first entry of the array $m[j], (j=1,2,...,ns=20,000)$ be $m [1]=m_0,$ where $ns$ is the number of iterations based on the MCMC method.

(c) For $j=2,..., ns,$ repeat the following steps:

(i) Generate $U_1$ from uniform(0,1).

(ii) Generate $V_1$ from gamma$(P_1,P_2).$

(iii) Let $p=\text{Min}[\frac{f_1(V_1| m[j-1])h_1(\alpha[j-1])}{f_1(\alpha[j-1]|m[j-1])h_1(V_1)},1]$

(iv) If $U_1< p, \alpha[j]=V_1,$ otherwise $\alpha [j]=\alpha [j-1],$ set $\alpha_0=\alpha [j]$

(v) Generate $W_1$ from the generalized gamma $(H\alpha_0+1, b, 1, a),$ set $m[j]=W_1,$ and $m_0=m[j]$

(vi) Let $j=j+1$ and go to (i).

The first initial $\tau$ = 5,000 estimates of $\alpha$ and $m$ are discarded as burned-in values, and Bayes estimates of the parameters based on sample $i (i=1,2,...,N)$ are found as

$$\hat{m_i}_B=\frac{\sum_{j=\tau+1}^{ns} m[j]}{ns-\tau}, \ \ \hat{\alpha_i}_{B}=\frac{\sum_{j=\tau+1}^{ns} \alpha[j] }{ns-\tau}.$$

5. Given $\hat{m_i}_B$ and $\hat{\alpha_i}_{B},$ the code uses (4) and the same process as in step 2 above (except this time Bayes estimates of the parameters $\alpha$ and $m$ are used in computations) to compute an approximate Bayes estimate $M_{(i,B)}(t)$ of the renewal function. The mean of $N$ Bayes estimates for the renewal function and ASE are computed

$$\begin{align} \overline{M}_{Bayes}(t)&=\frac{\sum_{i=1}^N M_{(i,B)}(t)}{N}, \\ \xi_{B}(\hat{M}(t))&=\frac{\sum_{i=1}^N(M_{(B,i)}(t)-M(t))^2}{N} \end{align}$$

The following tables provide a summary of the simulations. To compare the accuracy of the ML method versus the Bayesian method, the same sample $i (i=1,2,...,N)$ is used to compute both $M_{(i,B)}(t)$ and $M_{i,ML}(t).$ The line over a parameter estimate such as $\overline{\hat{m}_{ML}}$ denotes the average of estimated values in simulations.

Tables 1 and 3 reveal that as sample size $H$ increases, the ASE of the ML and Bayes estimates for $m, \alpha,$ and $M(t)$ decrease. Also, we can see that for a small sample size $(H<100),$ the Bayes estimator of $M(t)$ outperforms the ML, as the ASE of the Bayes estimator is much smaller. Furthermore, for large values of $H,$ the accuracies of the ML and Bayes estimators are comparable.

For the same generated samples used in Tables 1 and 3, we test the goodness of fit of the generated samples using Mathematica. Four distributions are considered: SPP, gamma, exponential, and inverse Gaussian. Given selected values of parameters $(m, \alpha),$ 100 samples of size $H$ are generated from the SPP, and the averages of the 100 p-values under the null hypotheses that the samples are from those distributions are listed in Tables 2 and 4. Tables 2 and 4 reveal that, in fact, the SPP is an appropriate distribution for the generated samples, as the p-values under the SPP distribution are large.

Table 5 shows a summary of the simulation studies when ML estimates are used to choose hyperparameter values (see step 3 in Section 5). The conclusion is the same as for Tables 1 and 3. That is, the Bayes estimator is more accurate than the ML estimator when the sample size $H$ is small. For larger samples, they have a similar precision.

To confirm convergence of MCMC-based estimators for $m$ and $\alpha,$ Figures 1 through 4, which are based on a generated sample of size $H$ from the SPP$(m=4.5, \alpha=3),$ respectively provide an MCMC trace plot for the Bayes estimate of $\alpha$ ("true"value 3.0), a histogram for estimated values of $\alpha,$ an MCMC trace plot for the Bayes estimate of $m$ ("true"value 4.5), and a histogram for estimated values of $m$ in $n_s -\tau=15,000$ iterations.

6. Numerical examples

EM-DAT provides worldwide data on dates and impacts of natural events from 1900 to the present. As an illustration of computations, we consider flood and tornado insurance damage in the United States from the year 2000 to the year 2020. The detailed flood and tornado data sets and examples of interarrival calculations are given in the appendix. Interarrivals between two insured damages are defined as the difference between the end date of an occurrence and the start date of the subsequent occurrence. To confirm that, in fact, both data sets fit SPP$(\alpha, m),$ we use the goodness-of-fit test from Gulati and Shapiro (2008), given below. A review of goodness-of-fit tests for the SPP distribution can be found in Chu, Dickin, and Nadarajah (2019). The data sets are sorted for computational convenience to find the ML estimate of $m.$ Using notation from Chu, Dickin, and Nadarajah (2019), let

$$\Lambda_0=\Lambda_1^2+\Lambda_2^2$$

where

$$\begin{align} \Lambda_1&=\sqrt{\frac{12}{H-1}}(\overline{U}-.5), \\ \Lambda_2&=\sqrt{\frac{5}{4(H+2)(H-1)(H-2)}}\\ &\quad (H-2+6H\overline{U}-12\sum_{i=1}^{H-1}\frac{i U_i}{H-1}), \end{align}$$

$$\begin{align} \overline{U}&=\sum_{i=1}^{H-1}\frac{U_i}{H-1}, \ \ Y_i^{*}=\sum_{j=1}^iY_i, \\ U_i&=\frac{Y_i^{*}}{Y_n^{*}},\ \ Y_i=(H-i+1)[X_{i:H}-X_{i-1,H}], \\ X_{0:H}&=0. \end{align}$$

Under the null hypothesis $H_0:$ the sample is from SPP$(\alpha,m),$ the test statistic is $\Lambda_0\sim \chi^2_2,$ and $H_0$ is rejected if $\Lambda_0 > \chi^2_{\nu, 2},$ where $\nu$ is a significance level. The following Mathematica code is used to find the p-value associated with a calculated value of $\Lambda_0.$

x: is the list for the sample observations

H is the sample size

Array[Y, H];

Array[U, H];

Array[YS, H];

Y[1] = H*x[[1]];

YS[1] = Y[1];

j = 2;

While[j < H + 1,

Y[j]= (H - j + 1) (x[[j]] - x[[j - 1]]);

YS[j] = Sum[Y[i], {i, 1, j}];

j++;

];

i = 1;

While[i < H + 1,

U[i] = YS[i]/YS[H];

i++;

];

Ubar = Sum[U[i]/(H - 1), {i, 1, H - 1}];

Λ₁ = Sqrt[12/(H - 1)] (Ubar - .5);

Λ₂ = Sqrt[5/(4 (H + 2) (H - 1) (H - 2))] 
  (H - 2 + 6 HUbar - 12 Sum[(i U[i])/(H - 1), 
  {i, 1, H - 1}]);

Λ₀ = Λ₁² + Λ₂²

P-value = 1 - CDF[ChiSquareDistribution[2], Λ₀] 

Flood data

The following denote $X_i, i=1,2,...,23$ values: 1.1, 1.2, 1.266, 2.066, 2.233, 2.433, 3.633, 5.6, 5.766, 6.5, 6.833, 7.5, 7.6, 8.666, 9.4, 10.733, 13.5, 18.9, 22.433, 23.628, 25.1, 30.996.

For the flood data, we get $\Lambda_0= .0004$ with p-value = .999. Therefore, there is no significant evidence against the SPP distribution. We have $H=23, \hat{m}_{ML}=1.1, \hat{\alpha}_{ML} =.5959\approx 1.$ As mentioned earlier, we can use ML estimates to choose hyperparameters $a, b, c,$ and $d.$ We let $a=.57, c=100, d=.01, P_1=1,000,$ and $P2=.001,$ and $b$ is calculated as $0.00024712.$ Using the above information, the Mathematica code based on $t=30$ provides

$$\hat{\alpha}_B=.99054, \ \ \hat{m}_B=1.490,\ \ \hat{M}_B(30)=7.28, \ \ \hat{M}_{ML}(30)=9.16.$$

Table 6 lists the p-values for goodness of fit of the flood data under four null hypotheses that data are from the four distributions (SPP, gamma, exponential, and inverse Gaussian). We can see that there is no strong evidence against SPP and exponential; however, the p-value for SPP is the largest among the four p-values.

Tornado data

The following denote $X_i, i=1,2,...,50$ values:

0.033, .1, 0.133, 0.166, 0.266, 0.3, 0.333, 0.366, 0.4, 0.4, 0.466, 0.566, 0.633, 0.766, 0.833, 0.9, 0.966, 1, 1.0333, 1.133, 1.2333, 1.266, 1.3, 1.466, 1.866, 2.066, 3.533, 4.233, 4.266, 5.133, 5.266, 5.333, 6.066, 7.366, 7.7, 8.066, 8.466, 8.466, 8.533, 8.9, 8.9, 9.533, 9.8, 10.8, 10.933, 13.066, 14.1, 16.23, 16.566, 24.133

For the tornado data, we get $\Lambda_0=.002$ with p-value = .9989. Therefore, there is no significant evidence against the SPP distribution. We have $H=50, \hat{m}_{ML}=.033$, and $\hat{\alpha}_{ML} =.2427\approx1$. As mentioned earlier, we can use ML estimates to choose the hyperparameters $a, b, c,$ and $d$. We let $a=.03, c=100, d=.01, P_1=1,000,$ and $P2=.001$, and $b$ is calculated as $0.0386637$. Using the above information, the Mathematica code based on $t=2$ provides

$$\begin{align} \hat{\alpha}_B=0.874733&, \ \ \hat{m}_B=0.0410618,\\ \hat{M}_B(20)=54.1&, \ \ \hat{M}_{ML}(20)=55.6 \end{align}$$

Table 7 lists the p-values for goodness of fit of the tornado data under four null hypotheses that data are from the four distributions (SPP, gamma, exponential, and inverse Gaussian). We can see that there is no strong evidence against the SPP distribution, as the p-values for the other three distributions are very small.

It is noted that $M(t)$ estimates for the tornado data are larger than for the flood data. This is because the average interarrival for the flood data is 9.36 and for the tornado data, 4.90. Therefore, for a specified $t$, the number of floods by time $t$ would be smaller than the number of tornados.

We also note that for both of the above data sets, the ML and Bayes estimates for $\alpha$ are less than 1. For computation of estimates of $M(t)$, as outlined in Section 2 using Ramsay (2006), we need $\alpha$ to be a positive integer. Therefore, in the Mathematica code both estimates are rounded up to 1 to obtain estimates for the renewal function. It is also noted that the expected value and the variance of the SPP random variable are defined, respectively, if $\alpha>1$ and $\alpha>2$. Hence, even though the first and the second moments are not defined when $0<\alpha<1,$ the delta moment for SPP($m,\alpha)$ is defined as

$$E[X^\delta]=\frac{\alpha m^\delta}{\alpha-\delta}, \hskip .1in 0 < \delta < \alpha.$$

A generated sample from SPP

The following data are generated from the SPP distribution using the built-in program in Mathematica, and the following denote $X_i, i=1,2,...,20$ values:

4.53749, 4.62575, 4.63825, 4.72415, 4.84672, 4.8955, 5.15921, 5.23818, 5.248, 5.39497, 5.54926, 5.62348, 5.69356, 5.72578, 5.87845, 6.42667, 6.51168, 8.86315, 9.56578, 16.6426.

For the above sample, we get $\Lambda_0=.0967$ with p-value = .952. Therefore, as expected, the interarrivals have an SPP distribution. The p-values for goodness of fit of gamma, exponential, and inverse Gaussian, respectively, are .0003, 3.8x$10^{-8}$, and .006.

$H=20, T=35.6201, \hat{m}_{Ml}=4.53749,$ and $\hat{\alpha}_{Ml}=3.7226\approx 4$. Based on ML estimates, we choose appropriate values for the hyperparameters, as follows:

$$\begin{align} a&=4.4, b=\frac{\ \hat{m}_{Ml}-a}{H \hat{\alpha}_{Ml}}=.00184, c=400, \\ d&=.01,, P_1=1,000, P_2=.004. \end{align}$$

Based on the above information, the Mathematica code computes the ML- and MCMC-based estimates for the renewal function at times $t=30$ and $t=50$ as

$$\begin{align} \hat{M}_B(30)=20.75, \ \ \hat{M}_{ML}(30)=20.70, \\ \hat{M}_B(50)=30.14, \ \ \hat{M}_{ML}(50)=30.11. \end{align}$$

Table 9 shows data extracted from FEMA.gov (https://www.fema.gov/openfema-data-page/fima-nfip-redacted-claims-v1) for hurricanes hitting New Orleans from September 29, 2021, through January 18, 2022. Table 10 shows the p-values for goodness of fit of the SPP, gamma, exponential, and inverse Gaussian distributions for the third column, interarrivals (days).

Table 10 reveals that whereas all four distributions can be used for the data, the p-value for the SPP distribution is the largest.

7. Summary

SPP distribution is considered as an inter-arrival distribution of a renewal process. Computational formulas from Ramsay (2006), for the sum of iid random variables from the SPP distribution, are used to provide ML- and MCMC-based Bayes estimates for the renewal function. The Bayes estimates of parameters $(m, \alpha)$ are based on shifted exponential and gamma priors. Two conditional posterior distributions $\pi(\alpha |\underline{x}, m)$ and $\pi(m |\underline{x},\alpha)$ are derived. Since $\pi(\alpha |\underline{x}, m)$ is not a recognized PDF, an MCMC method based on the Metropolis algorithm is employed. The parameter estimates are then used to compute the ML estimate and an approximate Bayes estimate for the renewal function. Simulation studies confirm that if selected values for hyperparameters of the prior distributions of $\alpha$ and $m$ are consistent with selected "true"values of $\alpha$ and $m$ in simulations, then the Bayes estimator for the renewal function would have a much smaller ASE and it outperforms the MLE with regard to accuracy for small samples. Also, for large samples, simulations confirm that the accuracies of the ML and Bayes estimators are comparable. In practice, only one sample is available from the SPP distribution, and choosing hyperparameters that are consistent with the ML estimates of $\alpha$ and $m,$ as outlined in Section 5, would provide more stable Bayes estimates for the renewal function, as demonstrated in Table 5. For illustration purposes of the computations, three data sets are considered: two data sets from EM-DAT and one data set from FEMA.gov for hurricanes hitting New Orleans.

---