Healthcare workers: diminishing burnout symptoms through self-care

1.   Introduction

Credibility theory is one of the most important and applicable methods of property and casualty ratemaking, particularly in the context of automobile insurance. It is called a "cornerstone" in the field of actuarial science, and has now been widely used in many insurance companies. Numerous authors devote themselves to developing tools to evaluate the existing credibility models, and then designing improved Bonus-Malus Systems (BMS) or No-Claim Discount (NCD). Contributions to the related literature include, among others, references ^[1,2,3,4].

As we all know, an insurance company insures many kinds of risks, and it is one of the main tasks of actuaries to analyze individual risks quantitatively, predict future claims and calculate equitable premium for each policyholder. As a rule, ratemaking process in non-life actuarial science always consists of two separate steps. In the first step, the policyholders are grouped into several classes on the basis of their observable characteristics, and those people in the same class are collected the same premium because they are considered to have similar expected losses. Examples of such risk characteristics are the age, gender and occupation of the policyholders, the type and use of their vehicles, the place where they reside, etc.

However, many important factors, such as respectfulness towards the law, reaction speed, aggressiveness behind the wheel and knowledge of the highway code, cannot be taken into account at this stage because they are usually impossible to measure and incorporate in a price list. Therefore, some heterogeneity still remains within the risk classes. Credibility models are used to distribute premiums fairly among a heterogeneous group of policyholders during the second step. The unobserved heterogeneity is taken into account by the introduction of a latent random variable $\Theta$, with Gamma, LogNormal or Inverse Gaussian distribution. Making use of Bayesian methodology, the base premium determined in the first step will be adjusted by discounts or surcharges depending on past claims of the policyholder. As result of this experience rating process, BMS or NCD is established to motivate people to drive more safely by rewarding good drivers and penalizing bad drivers.

Traditionally, actuaries assume that the random effect $\Theta$ in the credibility model does not vary with time. One of major problems in this so-called time-independent (or static) heterogeneous model is that the age of claims has not been taken into account. In other words, the classical models conflict with the intuition that the predictability based on a period of the policyholder's history should decrease with time. To overcome this defect, more and more actuaries have been paying attention to the evolutionary (or dynamic) credibility models, in which the unknown characteristics are represented by time-varying random effects $\{\Theta_t, t \geq 1\}$, see for instance ^{[5,6,7,8,9,10]} and the references therein. However, one drawback of these models is that the closed-form expressions for the predictive distribution and predictive premium are difficult to derive.

On the other hand, the dynamics on insurance data can also be obtained from endogenous approaches. As a typical representation among these methods, Gourieroux and Jasiak ^[11] has applied a heterogeneous integer-valued time series to describe the serial dependence in claim counts process of the policyholders. The authors show that their model is more flexible in computing predictive premiums by deriving the closed-form expressions of the Bayesian premium. They also find that the premium depends on not only the number of claims but also the timing of claim arrivals.

In the model of Gourieroux and Jasiak ^[11], the famous integer-valued autoregressive (INAR) process plays a central role. The classical INAR(1) process proposed by Al-Osh and Alzaid ^[12] can be written as

where $\phi \in [0, 1)$, and the so-called thinning operator "$\circ$" is defined by

in which $\{B_{t, k}, t = 2, 3, \ldots, k = 1, 2, \ldots\}$ is an array of independent and identically distributed (i.i.d.) Bernoulli random variables with mean $\phi$, and independent of the sequence $\{\epsilon_2, \epsilon_3, \ldots\}$.

Model (1.1) implies that the number of claims in period $t$ consists of two parts. The first part is a sum of Bernoulli variables indexed by the claims occurred in the preceding period. This term captures claim trends of the policyholders by providing a causality relationship and temporal dependence among the claim numbers of different periods. It can be seen that each of the claims occurred in the previous period could contribute one claim to the claim count of the next period, with probability $\phi$. This could be explained by the incentives to careful driving in automobile insurance, i.e., car accidents may modify a driver's perception of danger behind the wheel, such that a policyholder who reports claims in previous period would like to lower his risk level during next period, in consideration of safety and BMS. Furthermore, the other part $\epsilon_t$ is a number of claims without relationship with those occurred in the past, representing an innovation or fluctuations.

To generalize the application of model (1.1), Bermúdez et al. ^[13] adopts the bivariate INAR(1) process to price an automobile insurance contract with two types of coverage, in which the dependence structures arising from different sources of the claim numbers are taken into account. Afterwards, Bermúdez and Karlis ^[14] extends the results in the aforementioned paper to a multivariate setting using the multivariate INAR(1) process, and Zhang et al. ^[15] provides a flexible framework to fit the claim count data sets that contain a large number of zeros by modelling the innovations of multivariate INAR(1) process with zero-inflated distributions. Recently, Hu and Yao ^[16] proposes a modification of a combined INAR($p$) process to address higher-order dependence among the number of claims in an insurance portfolio.

In practice, the real-world data sets sometimes exhibit nonlinear structure, such as the so-called piecewise phenomenon, which means that the autoregressive parameter $\phi$ in model (1.1) may be affected by some factors, saying for example, the value of $N_{t-1}$. To capture this attribute, Monteiro et al. ^[17] proposes the integer-valued self-exciting threshold autoregressive process of order one with two regimes (referred to as SETINAR(2, 1) for short). In this paper, it is our goal to extend the classical SETINAR(2, 1) process to accommodate the unobserved heterogeneity and then apply the modified SETINAR(2, 1) process to credibility model for ratemaking problem. We will show that our approach provides a superior pattern to the predictive premiums calculation when compared to the outcomes of several alternative models.

The layout of this paper is organized as follows. In Sections 2, we introduce the proposed model and derive the Bayesian premium formula. In Section 3, some numerical examples are conducted to illustrate our methods and results. In section 4, we present an empirical analysis of insurance claims count data. Section 5 concludes this paper.

2.   Modelling and results

For automobile insurance, retemaking should be calculated based on both claim frequency and severity. However, it is believed that the number of claims can better reveal the inherent risk of a driver. In this paper, we follow most of the actuarial literature and focus only on the frequency part.

Let us consider one policyholder from an insurance portfolio, suppose $T$ to be the number of periods elapsed since this policy has been issued, and $N_{t}, t = 1, 2, \ldots, T$, denotes the number of claims reported by this policyholder during period $t$ of insurance. Our task is to draw prediction on the premium for the subsequent period $T+1$ of this policyholder. To this end, we assume that:

A1. The count variables $N_{1}, \ldots, N_{T}, N_{T+1}$ follow the SETINAR(2, 1) process:

where

● $I_{1, t} = I\{N_{t-1} \leq r \}, \; \; I_{2, t} = 1-I_{1, t} = I\{N_{t-1} > r\}$, in which $r$ is the so-called threshold variable;

● the thinning operator "$\circ$" is defined as

in which

are arrays of i.i.d. Bernoulli random variables with mean $\phi_j \in (0, 1)$, respectively. Furthermore, {$\{B^{(1)}_{t, k}, k = 1, 2, \ldots\}$ and $\{B^{(2)}_{t, k}, k = 1, 2, \ldots\}$ are supposed to be independent of $\epsilon_t$ and $\{N_{t-1}, \ldots, N_1\}$ for any fixed $t$.

A2. Given the unobservable heterogeneity $\Theta = \theta$, $N_{1}$ follows the Poisson distribution with mean $\lambda \theta$, i.e.,

where $\lambda = e^{\mathit{\boldsymbol{\beta}}\mathit{\boldsymbol{x}}'_{ 1}}$, in which $\mathit{\boldsymbol{x}}_{ 1} = (1, x_{ 1, 1}, x_{ 1, 2}, \ldots, x_{ 1, q})$ is the observable risk characteristics of the policyholder in the first period, $\mathit{\boldsymbol{\beta}} = (\beta_0, \beta_1, \ldots, \beta_q)$ collects the regression coefficients and $\mathit{\boldsymbol{x}}_{ 1}'$ denotes the transpose of $\mathit{\boldsymbol{x}}_{ 1}$.

A3. Given the unobservable heterogeneity $\Theta = \theta$, $\{\epsilon_{t}, t = 2, 3, \ldots\}$ is a sequence of independent non-negative integer-valued random variables confirming to the Poisson distribution with mean $\eta_t \theta$, where $\eta_t = e^{\mathit{\boldsymbol{\omega}}\mathit{\boldsymbol{x}}'_t}$ with $\mathit{\boldsymbol{x}}_{ t} = (1, x_{ t, 1}, x_{ t, 2}, \ldots, x_{ t, q})$ representing the characteristics of the policyholder in the period $t$ and $\mathit{\boldsymbol{\omega}} = (\omega_0, \omega_1, \ldots, \omega_q)$ denoting the regression coefficients. Furthermore, $\epsilon_{t}$ and $N_{1}, \ldots, N_{t-1}$ are independent conditional on $\Theta$ and $\mathit{\boldsymbol{x}}_{ t}, \ldots, \mathit{\boldsymbol{x}}_{ 1}$ for any fixed $t$.

A4. The unobservable heterogeneity $\Theta$ follows Gamma distribution whose density function is

in which the shape parameter and scale parameter are identical so that the premium for this policyholder in the first period results from a priori rating system and equals to $\lambda = e^{\mathit{\boldsymbol{\beta}}\mathit{\boldsymbol{x}}'_1}$.

Remark 2.1. The model defined by (2.1) implies that the autoregressive parameter is a function of $N_{t-1}$, i.e., the effect of $N_{t-1}$ on $N_t$ varies with different values of $N_{t-1}$. In the insurance practice, this could be explained that policyholders may adjust their efforts to prevent losses according to their experience with past claims. Taking the NCD system of automobile insurance in China for example, if a policyholder has reported more than one claim, his premium will be raised up to 2 times of the base premium. Therefore, this policyholder might become more risk-averse and may have a greater of initiative to reduce the claim numbers in the future, so that it is reasonable for us to believe the constraint condition of ${\phi_1} > \phi_2$. By contrast, more claims mean higher risks, i.e., the claims of period $t-1$ could produce another claim in period $t$ with a higher probability, making us could acknowledge that ${\phi_1} < \phi_2$. In brief, it is very necessary to consider the SETINAR(2, 1) process with two different autoregressive coefficients for actuarial applications from a practical view.

Remark 2.2. The distributional assumptions A2–A4 on top of SETINAR(2, 1) process are made in accordance with the Poisson-Gamma model that have been widely applied for modeling claim counts when the portfolio is heterogeneous to construct the optimal experience rating system. Under these assumptions, we are able to derive the closed-form expression for Bayesian pricing formula and illustrate more effectively our main results regarding the threshold structure of SETINAR(2, 1) process. However, the limitations of these assumptions and the potential impacts they may have on the usefulness and applicability of our proposed model are also very obvious. For better and more flexible application in practice, our proposed model can be generalized to some alternatives to the Poisson-Gamma assumption, the computations become however more complicated, because numerical methods to perform integration are usually required. Perhaps credibility premium, the linear estimator (based on past claims data) that best approximates the Bayesian premium, can be considered to overcome these difficulties in these cases.

From Denuit et al. ^[3], we know that given the observations of claim counts for $T$ periods $N_{1} = n_{ 1}, \ldots, N_{T} = n_{ T}$, the Bayesian premium of the policyholder for period $T+1$ with the quadratic loss function can be obtained as

As shown in (2.4), the pure premium depends on the conditional distribution of the heterogeneity given the claim history. Therefore, we obtain the analytical expression of the conditional distribution of the heterogeneity in the following theorem.

Theorem 2.1. Under the assumptions A1–A4, we have that:

(1) For $T = 1$, given $N_{1} = n_{ 1}$, the heterogeneity $\Theta$ follows Gamma distribution with conditional density function

(2) For any $T \geq 2$, given $N_{1} = n_{ 1}, \ldots, N_{T} = n_{ T}$, the conditional density function of $\Theta$ is given by

in which

and $\gamma(\theta; \tilde{\alpha}_{ 1}, \tilde{\alpha}_{ 2})$ is the density function of Gamma distribution with shape parameter $\tilde{\alpha}_{ 1}$ and scale parameter $\tilde{\alpha}_{ 2}$.

Proof. For $T = 1$, because $N_{1}|\Theta = \theta \sim \mathrm{P}\left(\lambda\theta\right)$, we know that the joint distribution of $N_{1}$ and $\Theta$ is

then, we have

which implies that (2.5) holds.

Similarly, for any $T \geq 2$, we have

in which

Noting that

is the probability mass function of the sum of two discrete random variables, we can rewrite it as the following convolutions:

By applying (2.11), (2.13) and (2.14) into (2.12), we can obtain (2.6). □

One can see that the conditional distribution of the heterogeneity $\Theta$ is a mixture of Gamma distributions with parameters that depend on the claim history. From (2.4) and the above theorem, it is easy to derive the analytical formulas of the pure premium for period $T+1$ as follows.

Theorem 2.2. If the assumptions A1–A4 hold, then we have:

(1) When $T = 0$, no claim history is available at the beginning of the contract, and the base premium determined in the first step of ratemaking process is collected at the time of signing the insurance contract, i.e.,

(2) When $T = 1$, the Bayesian premium of the policyholder for period $T = 2$ is

where the predicted heterogeneity $\hat{\Theta}_{1}$ is given by

(3) When $T \geq 2$, the {Bayesian premium} of the policyholder for period $T+1$ is

where the predicted heterogeneity $\hat{\Theta}_{T}$ is given by

in which $\tilde{\alpha}_{ 1}(z_{ T}, \ldots, z_{ 2}, n_{ T}, \ldots, n_{ 1})$, $\tilde{\alpha}_{ 2}$ and $\pi(z_{ T}, \ldots, z_{ 2}, n_{ T}, \ldots, n_{ 1})$ are accordingly defined by the Eqs (2.8)–(2.10).

3.   Numerical illustration

In this section, we present some simulation results to illustrate the premium updating for one given policyholder. To this end, we fix all the values of the parameters, compute and compare the Bayesian premiums for different models.

We mainly focus on the following three different models for the number of claims.

Model 1: Credibility model based on SETINAR(2, 1) process with Gamma heterogeneity, abbreviated as CM-SETINAR(2, 1). In this model. we assume that $\lambda = 0.4286$, $\alpha = 9$ and $\eta_t = \eta = 0.3$ for $t = 2, 3, \ldots.$ From Remark 1, we know that the policyholders who report more than one claim would be penalized by some premium surcharges. Therefore, we take $r = 1$ to illustrate the results. Furthermore, in order to compare different situations, we consider the following two scenarios for $\phi_1$ and $\phi_2$:

(a) $\phi_1 = 0.3$, $\phi_2 = 0.2$;

(b) $\phi_1 = 0.3$, $\phi_2 = 0.4$.

Model 2: Credibility model based on INAR(1) process with Gamma heterogeneity, abbreviated as CM-INAR(1). When $\phi_1 = \phi_2 = \phi$, CM-SETINAR(2, 1) reduces to the case studied in Gourieroux and Jasiak ^[11]. For the other parameters, we assume that $\lambda = 0.4286$, $\phi = 0.3$, $\alpha = 9$ and $\eta_t = \eta = 0.3$ for $t = 2, 3, \ldots.$

Model 3: Credibility model based on mixed Poisson distribution with Gamma heterogeneity, abbreviated as CM-MP. When $\phi_1 = \phi_2 = 0$, Model 1 becomes the standard case, i.e., negative binomial distribution. Furthermore, let $\lambda = 0.4286$ and $\alpha = 9$.

Remark 3.1. In Gourieroux and Jasiak ^[11], the number of claims is fitted by the following model:

where conditional on the unobserved heterogeneity $\Theta = \theta$, the variables $\{\epsilon_t, t = 2, 3, \ldots\}$ are independent with distribution $P(\eta \theta)$, so that $\{N_t, t = 1, 2, \ldots\}$ has Poisson marginal that is distributed by $P(\frac{\eta}{1-\phi} \theta)$, to make the process itself satisfy the stationarity. Therefore, we also set $\lambda = \frac{\eta}{1-\phi}$ in our model for the convenience of comparison. Additionally, the mixed Poisson distribution with Gamma heterogeneity implies

which results in the standard negative binomial model.

Remark 3.2. Note that we fix $\eta_t = \eta$ for $t = 2, 3, \ldots$ in CM-SETINAR(2, 1) and CM-INAR(1), so that we can use the same setting of the parameters with Gourieroux and Jasiak ^[11] to calculate the premiums, and show the effect of the threshold on the results. In practice, $\eta_t$ could vary because the risk characteristics (e.g., age, coverage, etc.) of the policyholder may be different for every period. In this case, the claims history $\{N_t, t = 1, 2, \ldots\}$ may be not stationary any more. As stated in Asmussen ^[18], the transient distributions in a Bonus-Malus system are quite far from the stationary ones usually, i.e., the rate of convergence to stationarity may be slow in comparison to the typical sojourn time of a customer in the portfolio. We accept this view in this paper and relax the condition of stationarity. In the next section (Empirical Study), we will illustrate that our proposed model is able to fit the real data better than some existing models.

As an illustration, we consider several different claim histories, and provide a summary of the evolution of the corresponding premiums over the first $4$ periods in Table 1.

From Table 1, it is easy to see that the predicted premiums of CM-SETINAR(2, 1) and CM-INAR(1) are determined by the explicit claim frequency history, and not just by the total number of claims as that in the case of CM-MP. For example, we compare the predicted premiums based on claim histories $(0, 2)$ and $(2, 0)$, respectively. For CM-MP, $P_3$ has the same value 0.4783 because the total numbers of claims in the first 2 periods are both $2$. However, for CM-SETINAR(2, 1) and CM-INAR(1), the two claim histories are distinguishable and the claim arrival (i.e., in the first or in the second year) matters. Therefore, the predicted premiums are different. Furthermore, $P_3$ for claim history $(0, 2)$ is higher than $P_3$ for claim history $(2, 0)$ because the former is the premium after a claim-free period and a second period with $2$ claims, while the later is the premium after a first period with $2$ claims and a free-claim period. This result shows that the most recent claims are more predictive than older ones.

Let us now consider the claim history $(1, 2, 0)$ to present the difference of predicted premiums between CM-SETINAR(2, 1) (including Scenario (a) and (b)) and CM-INAR(1).

Ⅰ. Initially, all the models have the same premium 0.4286, because all of them calculate $P_1$ according to $\lambda$.

Ⅱ. If there is Ⅰ claim filed during the first period, $P_2$ for CM-SETINAR(2, 1) (a), CM-SETINAR(2, 1) (b) and CM-INAR(1) also have the same value 0.6182, because $N_1$ does not exceed the threshold value $r = 1$ at this moment, implying that both CM-SETINAR(2, 1) (a) and CM-SETINAR(2, 1) (b) predict $N_2$ according to $N_2 = \phi_1 \circ N_{1}+\epsilon_2,$ which is identical to that of CM-INAR(1).

Ⅲ. As time passes by, when Ⅱ claims are reported in the next period, and then a claim-free period comes, the predicted premiums for CM-SETINAR(2, 1) (a), CM-SETINAR(2, 1) (b) and CM-INAR(1) become quite different.

(1) The values of $P_3$ for CM-SETINAR(2, 1) (a) and CM-INAR(1) are 0.7479 and 0.9479, respectively. Moreover, it can be calculated that the surcharges for these two models accordingly are $(0.7479 - 0.6182)/0.6182 = 20.98%$ and $(0.9479 - 0.6182)/0.6182 = 53.33%$, implying that the policyholder for CM-SETINAR(2, 1) (a) will get milder penalty than the individual for CM-INAR(1), which is in line with the situation $\phi_1 > \phi_2$ (see Remark 1). Correspondingly, after a claim-free period, the policyholder for CM-SETINAR(2, 1) (a) will receive lower reward than the driver for CM-INAR(1), since the discounts for these two models are $(0.7479 - 0.3374)/0.7479 = 54.89%$ and $(0.9479 - 0.3374)/0.9479 = 64.41%$, respectively.

(2) By contrast, for CM-SETINAR(2, 1) (b) with $\phi_1 < \phi_2$, the policyholder will get more severe penalty or receives higher reward than driver that is described in CM-INAR(1).

To sum up, compared to CM-INAR(1), CM-SETINAR(2, 1) (a) brings less variation to the premiums, while CM-SETINAR(2, 1) (b) causes more variation to the premiums. We can also get these trends visually from Figure 1. The same conclusion can also be made from the results of the other claim histories, revealing that the threshold structure has significant implications on the dynamics of the predictive premiums.

4.   Empirical study

In Section 3, we fix the values of the parameters for different models in a special way to numerically illustrate the calculation of the Bayesian premiums for one policyholder. In this section, we implement all the models presented in this paper with a sample of insurance data, obtain the statistical estimates of the parameters and compare the three models in a more scientific way.

Suppose that the insurance portfolio consists of $I$ policies, denote the unobservable heterogeneity of the $i$th policyholder by $\Theta_i$, $i = 1, \ldots, I$ and define $N_{i, t}$ to be the claim count of the $i$th policyholder in period $t$, $t = 1, \ldots, T_i$. Then, for each policyholder in the portfolio, the related variables can be expressed as $(\Theta_i, N_{i, 1}, \ldots, N_{i, T_i})$, in which the time length $T_i$ may differ among the policyholders in the case when the data are unbalanced. Moreover, the sequences $\{(\Theta_i, N_{i, 1}, \ldots, N_{i, T_i}), i = 1, 2, \ldots, I\}$ are supposed to be independent at the portfolio level.

Rewrite the parameters to be estimated as $\mathit{\boldsymbol{\Phi}} = (\mathit{\boldsymbol{\beta}}, \mathit{\boldsymbol{\omega}}, \alpha, \phi_1, \phi_2)$ in our proposed model, and let $\mathit{\boldsymbol{n}}_i = (n_{i, 1}, \ldots, n_{i, { T}_i})$ denote the observations of $\mathit{\boldsymbol{N}}_i = (N_{i, 1}, \ldots, N_{i, { T}_i})$, then the likelihood for the $i$th policyholder could be shown as

in which $\pi(z_{ i, T_i}, \ldots, z_{ i, 2}, n_{ i, T_i}, \ldots, n_{ i, 1})$ is defined by (2.10). Then, we can obtain the estimators $\hat{\mathit{\boldsymbol{\Phi}}} = (\hat{\mathit{\boldsymbol{\beta}}}, \hat{\mathit{\boldsymbol{\omega}}}, \hat{\alpha}, \hat{\phi}_1, \hat{\phi}_2)$ by maximizing the following total log-likelihood function of the observed data derived by summing up the logarithm of (4.1) over all policyholders:

In the following part, we apply our proposed model to a panel data set for demonstrating our methods and results. The data are drawn from Wisconsin Local Government Property Insurance Fund (LGPIF), and have been used previously in Frees et al. ^[19], Quan and Valdez ^[20], Oh et al. ^[21,22] and Chen et al. ^[23]. The LGPIF provides property insurance such as buildings insurance, vehicles insurance and equipment insurance, for various governmental entities, mostly including counties, cities, towns, villages, school districts, fire departments and other miscellaneous entities. Detailed information on the LGPIF data can be found in Frees et al. ^[19]. There are six types of insurance coverage in the data, and here we will consider only the collision coverage for new vehicles (i.e., vehicles of current model year, or 1$\sim$2 years prior to the current model year), which provides coverage for the impact of a vehicle with an object, impact of vehicle with an attached vehicle, or the overturn of a vehicle. We first remove the policyholders who have zero coverage and then take the policyholders that have complete claim histories over the policy years from 2006 to 2010 as the training dataset for simplicity. Therefore, a longitudinal data with 199 governmental entities is used for our analysis to calibrate the model. We describe some summary statistics for the numbers of claims and the risk characteristics in Tables 2 and 3, respectively. As can be seen from the results, there is a clear increasing trend in the frequency, implying that the sequence of claim numbers for a policyholder could be nonstationary during these years.

Table 4 displays the effects of risk characteristics on the numbers of the claims, showing that theses variables are significant. For one thing, it seems that the county entities have the highest accident rate, while the town entities have the lowest risk. For another, we separate all the policyholders into two classes according to the median of the coverage, so that CoverageCN1 means the level with coverage $\in (0, 0.2051]$, and CoverageCN2 indicates the level with coverage $\in (0.2051, +\infty)$. As we can see, the drivers with higher coverage have more accidents, probably because of the moral hazard, i.e., the coverage of insurance provides incentives for reckless driving behavior, and make the policyholders have less risk aversion towards the loss.

In Frees et al. ^[19], the authors fit several commonly employed count models which include Poisson, zero-inflated Poisson, zero-one-inflated Poisson, negative binomial, zero-inflated negative binomial and zero-one-inflated negative binomial to the LGPIF dataset, and use Chi-square goodness-of-fit statistics to compare these different models. Their results show that negative binomial is significantly better than other models apart from zero-one-inflated negative binomial. The authors prefer and recommend the simpler model, so negative binomial is selected for further investigation. Follow this related study, it is reasonable for us to take CM-MP as the benchmark, then compare the fitting effect and predictive performance of CM-SETINAR(2, 1) and CM-INAR(1) with CM-MP, to illustrate our main discoveries.

We now apply the three models aforementioned, i.e., CM-SETINAR(2, 1), CM-INAR(1) and CM-MP, to the LGPIF data. The parameters are estimated via the MLE (maximum likelihood estimation) method based on the log-likelihood function defined by (4.2), and the estimates obtained using optim function in R are presented in Table 5. In the meantime, we also exhibit two goodness-of-fit measures, including the log-likelihood function calculated at the estimated parameters (LogLik) and the Akaike information criterion (AIC), to evaluate different models which we focus on. It can be seen that both of the criteria support our proposed model, because it has the largest value of LogLik and the smallest value of AIC.

It should be noted that we have to determine the value of the threshold parameter $r$ in our proposed CM-SETINAR(2, 1) before estimating the other parameters. To this end, we can search the integer that maximizes the LogLik criterion function, in an appropriate subset of positive integers $[1, 2\ldots, R]$. For the LGPIF data, we may set $R = 14$ because a maximum of 15 claims has been filed by the policyholders. Therefore, the algorithm to get the unknown threshold $r$ can be summarized in the following two steps:

Step 1. For each fixed $r \in [1, 2, \ldots, 14]$, find $\hat{\mathit{\boldsymbol{\Phi}}}^{(r)}$ such that

in which $l(\mathit{\boldsymbol{\Phi}})$ is defined by (4.2).

Step 2. The threshold parameter is estimated by searching over all candidates, i.e.,

In practice, we also can search the threshold parameter by minimizing the values of AIC or optimize other criterion functions. Figures 2 and 3 show different $r$ and the corresponding values of LogLik and AIC, respectively. It is found that $r = 6$ is the best choice. Accordingly, we can then obtain the estimates of the other parameters, as shown in Table 5.

At last, we compare the predictive performance of CM-SETINAR(2, 1), CM-INAR(1) and CM-MP by performing an out-of-sample experiment. To this end, the observations in year 2011 with 1098 policyholders are considered. Of these samples, only 174 policyholders are identical to the governmental entities in the training dataset. Therefore, a longitudinal data with the 174 samples is reserved as the testing dataset for model validation. Since the dataset used is balanced now, we denote $N_{j, T+1}$ to be the claim number in the testing period of the $j$th policyholder, $j = 1, \ldots, J$, where $J$ is the size of the testing dataset and equals to 174 in our analysis. Correspondingly, the observation of $N_{j, T+1}$ is denoted by $n_{ j, T+1}$ In insurance claims modeling, it is more important to check the overall distribution for all policyholders rather than prediction of the claim frequency for each policyholder, which could be more useful for premium calculation, risk management and so forth. As a result, we adopt three measures to assess the predictive performance for each model, i.e., the log-likelihood function (LogLik), the mean squared error of prediction (MSPE) and the mean absolute error of prediction (MAPE), which are accordingly defined by

From the results displayed in Table 6, it is clear that our proposed model, CM-SETINAR(2, 1), has the best predictive performance with the largest LogLik among other two models. Furthermore, the smallest MSPE and MAPE also show that CM-SETINAR(2, 1) outperforms other two models, which is consistent with the fitting result in Table 5.

5.   Conclusions

In this paper, we apply the integer-valued self-exciting threshold autoregressive time series with Gamma distribution for the heterogeneity to the credibility theory and propose a credibility model based on heterogeneous SETINAR(2, 1) process. We derive the explicit Bayesian premium formula, and consider the ratemaking problem of pricing the automobile insurance contract. Numerical examples show that our proposed model can account for the claim history and the behavior of the policyholders more flexibly. A real data analysis demonstrates that our approach provides a superior pattern to the Bayesian premium calculation when compared to the outcomes of two existing models, in view of goodness of fit and predictive ability.

Finally, the study of this paper might be extended in several directions. First, some other distributions, such as inverse Gaussian, LogNormal, mixed Gamma, etc., could be taken into account for the unobserved heterogeneity. Besides, the dynamic heterogeneity is also a very interesting topic for further research. Second, higher-order processes could be adopted to make the credibility models more practical. Third, multivariate integer-valued time series models could be considered to deal with the contracts with different types of coverage.

Use of AI tools declaration

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

Acknowledgments

The authors would like to thank the editor and the referees for their constructive and pertinent comments and suggestions.

Conflict of interest

The authors declare that they have no competing interests.