Remarks

• • Our Cox model generalizes the framework of Badescu et al. (Reference Badescu, Lin and Tang2016) by modeling claim arrivals at the individual policyholder level, using covariates $\boldsymbol{x}_i$ to capture heterogeneity across risk profiles. This allows the claim intensity to reflect both individual characteristics and broader environmental variation such as seasonality. Moreover, in contrast to their assumption of an Erlang-distributed intensity over $[d_l, d_{l+1})$ , we model intensity as a covariate-dependent constant within each period. This choice simplifies estimation and leads to a Poisson mixture distribution, which is appropriate for overdispersed claim count data. While the Erlang assumption yields a Pascal mixture, it adds complexity without a clear improvement in model performance. Moreover, with appropriate segmentation of the timeline, the constant-intensity assumption remains a reasonable approximation.

• • A special case of our framework arises when modeling claim arrivals across lines of business, with index i denoting each line. If the regression function $\lambda^{(j)}(\!\cdot\!)$ is constant, the model reduces to the contemporaneously conditionally independent multivariate HMM of MacDonald and Zucchini (Reference MacDonald and Zucchini2016). This HMM structure would capture dependence between lines of business through shared hidden states – an important feature given the emphasis in the literature on modeling such dependencies in reserving. Alternatively, i can represent claim types within a single line of business, allowing the model to account for heterogeneity in mark distributions across types.

• • While our model assumes predetermined time intervals $[d_{l-1}, d_l)$ , it would be ideal to select these intervals in a data-driven manner to better capture changes in intensity over time. Nonetheless, using fixed intervals – such as weekly, monthly, or quarterly – offers a practical and sufficiently flexible solution in most applications.

4.2.1. The multinomial model

With the aim of joint modeling of frequency and delay, we deviate from the continuous-time framework to a discrete-time framework by aggregating the events towards a coarser time scale. To this end, we assume that claims can be reported with a maximum delay of D periods. This means that if a claim occurs in period t, it will be reported in periods $t+d$ ( $d=0,1,\dots,D$ ). Note that a period t corresponds to the time interval $[d_{t-1},d_t)$ . We also, for simplicity, assume that the predetermined time points $d_t$ $(t=0,\dots,T)$ are equidistant, so that $d_0=0$ , $d_{t+1}-d_{t}=1$ and T represent the last observed period before the valuation date $\tau$ . That is, we assume that we consider periods of the same length, where each period represents one unit of time.

We denote the probability that a policyholder’s claim with risk attributes $\boldsymbol{x}_i$ that occurred in period t is reported in period $t+d$ by $p_{t}(d;\,\boldsymbol{x}_i)$ . It is easy to show that the discretely observed claim arrival process, $\{N_{i,1},N_{i,2},\dots,N_{i,T} \}$ , the discretely observed reported claim process, $\{N_{i,1}^{\text{r}},\dots,N_{i,T}^{\text{r}} \}$ , and the discretely observed IBNR claim process, $\{N_{i,1}^{\text{ibnr}},\dots,N_{i,T}^{\text{ibnr}}\}$ , fall under the class of Poisson-HMMs, with the same hidden parameter process $\{C_1,C_2,\dots\}$ as $N(t;\,\boldsymbol{x}_i)$ , and with state-dependent Poisson intensities given by

\begin{align*}\tilde{\lambda}_t^{(j)}(\boldsymbol{x}_i)=e_{i,t}\lambda^{(j)}(\boldsymbol{x}_i), \quad \mu_t^{(j)}(\boldsymbol{x}_i) = e_{i,t}\lambda^{(j)}(\boldsymbol{x}_i)p_{i,t}^{\text{r}}, \quad \nu_t^{(j)}(\boldsymbol{x}_i) =e_{i,t} \lambda^{(j)}(\boldsymbol{x}_i)(1-p_{i,t}^{\text{r}}),\end{align*}

respectively, where $p_{i,t}^{\text{r}} = \sum_{d=0}^{\min\{D,T-t\}}p_{t}(d;\,\boldsymbol{x}_i)$ is the probability that a claim from policyholder i that occurred in period t is reported before the reserve valuation date.

Let $Z_{i,t,d}=\sum_{k=1}^\infty 1\left\{T_{i,k} \in [t-1,t), T_{i,k}+U_{i,k} \in [t+d-1, t+d)\right\}$ denote the number of claims from policyholder i that occurred in period t and are reported in period $t+d$ . Note that $N_{i,t}^{\text{r}} = \sum_{d=0}^{\min\{D,T-t\}}Z_{i,t,d}$ . Our observed data are now given by

\begin{align*}\mathcal{O}_D = \{n_{i,t}^{\text{r}},z_{i,t,d}\mid i=1,\dots,m;\,t = 1,\dots,T;\,d=0,\dots,D;\, t+d\le T \},\end{align*}

where $n_{i,t}^{\text{r}}$ and $z_{i,t,d}$ are the realizations of $N_{i,t}^{\text{r}}$ and $Z_{i,t,d}$ , respectively.

Similar to Verbelen et al. (Reference Verbelen, Antonio, Claeskens and Crevecoeur2022), we make the assumption that

\begin{align*}(Z_{i,t,0}, \dots, Z_{i,t,\min\{D,T-t\}}) \,\big|\, N_{i,t}^{\text{r}} \sim \text{Multinomial}\left(N_{i,t}^{\text{r}}, \left(\frac{p_t(d;\,\boldsymbol{x}_i)}{p_{i,t}^{\text{r}}}\right)_{d=0}^{\min\{D,T-t\}}\right).\end{align*}

Note that for $t \le T-D$ , all claims occurred in period t is fully observed, that is, $N_{i,t}^{\text{r}} = N_{i,t}$ , and the counts $(Z_{i,t,0},\dots,Z_{i,t,D})\mid N_{i,t}^{\text{r}}\sim \text{Multinomial}(N_{i,t}^{\text{r}},\boldsymbol{p}_t(\boldsymbol{x}_i)),$ where $\boldsymbol{p}_t(\boldsymbol{x}_i)=(p_t(0;\,\boldsymbol{x}_i),\dots,p_t(D;\,\boldsymbol{x}_i))$ .

The likelihood of the observed data is then given by

(4.2) \begin{equation} \begin{split} \mathcal{L}^{(2)}(\boldsymbol{\Phi}_2\mid \mathcal{O}_D) & \propto P\left(\boldsymbol{N}_1^{({\text{r}},1:T)}=\boldsymbol{n}_1^{(\text{r},1:T)},\dots,\boldsymbol{N}_m^{(\text{r},1:T)}=\boldsymbol{n}_m^{(\text{r},1:T)}\right) \\ &\times \prod_{i=1}^m \prod_{t=1}^T \frac{n_{i,t}^{\text{r}}!}{\prod_{d=0}^{\min\{D,T-t\}}z_{i,t,d}!} \prod_{d=0}^{\min\{D,T-t\}} \left(\frac{p_{t}(d;\,\boldsymbol{x}_i)}{p_{i,t}^{\text{r}}}\right)^{z_{i,t,d}}, \end{split}\end{equation}

where $\boldsymbol{\Phi}_2$ represents the vector of the model’s parameters, including the initial distribution $\boldsymbol{\pi}_1$ and the transition matrix $\boldsymbol{\Gamma}$ of the HMM, the regression coefficients $\boldsymbol{\theta}_j$ for the $\lambda^{(j)}$ ’s, and the parameters related to the modeling of $\boldsymbol{p}_t(\boldsymbol{x}_i)$ , denoted by $\boldsymbol{\delta}$ . An equivalent way to write the multinomial component of the likelihood above is to decompose it into a series of conditional binomial likelihoods. That is, we can write $\mathcal{L}^{(2)}$ as:

(4.3) \begin{equation} \begin{split} \mathcal{L}^{(2)}(\boldsymbol{\Phi}_2\mid \mathcal{O}_D) & \propto P\left(\boldsymbol{N}_1^{({\text{r}},1:T)}=\boldsymbol{n}_1^{(\text{r},1:T)},\dots,\boldsymbol{N}_m^{(\text{r},1:T)}=\boldsymbol{n}_m^{(\text{r},1:T)}\right) \\ &\times \prod_{i=1}^m \prod_{t=1}^T \prod_{d=1}^{\min\{D,T-t\}} \binom{\sum_{j=0}^{d}z_{i,t,j}}{z_{i,t,d}} q_t(d;\,\boldsymbol{x}_i)^{z_{i,t,d}}(1-q_t(d;\,\boldsymbol{x}_i))^{\sum_{j=0}^{d-1}z_{i,t,d}}, \end{split}\end{equation}

where

\begin{align*}q_t(d;\,\boldsymbol{x}_i) = \frac{p_t(d;\,\boldsymbol{x}_i)}{\sum_{j=0}^{d}p_t(j;\,\boldsymbol{x}_i)}\end{align*}

represents the conditional probability that a claim incurred at period t by policyholder with risk attributes $\boldsymbol{x}_i$ is reported in period $t+d$ , given that the claim is reported by period $t+d$ . Thus, we break down the multinomial likelihood by considering the binomial likelihood of observing each count given the cumulative counts and adjusting the success probabilities accordingly. Note that $p_t(d;\,\boldsymbol{x}_i)$ can be obtained iteratively from $q_t(d;\,\boldsymbol{x}_i)$ by:

\begin{align*} & p_t(D;\,\boldsymbol{x}_i) = q_t(D;\,\boldsymbol{x}_i),\quad p_t(d;\,\boldsymbol{x}_i) = q_t(d;\,\boldsymbol{x}_i)\times \left(1-\sum_{j=d+1}^D p_t(j;\,\boldsymbol{x}_i)\right) \quad (d=D-1,\dots,1), \\ & p_t(0;\,\boldsymbol{x}_i) = 1 - \sum_{j=1}^D p_t(j;\,\boldsymbol{x}_i).\end{align*}

Remark

Note that t used in a discrete-time setting differs from t in $ N_i(t) = \sum_{j=1}^\infty \unicode{x1D7D9}\{T_{ij} \lt t\} $ . In the continuous-time setting, t represents calendar time and $ N_i(t) $ counts the number of claims of the policyholder i that occurred before t. In the discrete-time setting, t represents calendar time discretized into periods (e.g., months), where $ t = 1 $ corresponds to the first period, $ t = 2 $ to the second, and so on. For simplicity, we use the same notation t in both settings, but its meaning is clear from the context.

4.2.2. The Dirichlet-multinomial model

The probabilities of reporting delay often exhibit temporal variability, which can be characterized as either explainable or unexplainable. Explainable variations in reporting delay are systematic variations that arise from differences in the period of occurrence of claims and/or the risk attributes associated with policyholders. For example, when working with daily data, claims occurring on weekends may exhibit different reporting patterns compared to those occurring on weekdays. Similarly, the month in which a claim occurs can also influence reporting delay probabilities. In contrast, unexplainable variability refers to unsystematic fluctuations in reporting delay that are not captured by our covariates. In cases where data exhibit high random variability in reporting delay, both continuous-time and multinomial models may produce interval estimates for IBNR claim counts that are significantly off from the actual IBNR claim count, as we will observe in our data analysis. Although an exact cause for the high fluctuations cannot be precisely identified without additional data collection, one possible explanation is the presence of unobserved covariates that have direct and interaction effects on the delay probabilities.

To this end, we add randomness to the reporting delay probabilities $p_{t}(d;\,\boldsymbol{x}_i)$ ’s by assuming that $\boldsymbol{p}_t(\boldsymbol{x}_i) = (p_{t}(0;\,\boldsymbol{x}_i),\dots,p_{t}(D;\,\boldsymbol{x}_i))$ follows a Dirichlet distribution with parameters $\boldsymbol{\eta}_{i,t} = (\eta_{t,0}(\boldsymbol{x}_i),\dots,\eta_{t,D}(\boldsymbol{x}_i))$ . This choice is motivated by the fact that the Dirichlet distribution serves as a conjugate prior for the multinomial distribution, which significantly simplifies the mathematical treatment and estimation of our model.

Parameters $\boldsymbol{\eta}_{i,t} = (\eta_{t,0}(\boldsymbol{x}_i), \dots, \eta_{t,D}(\boldsymbol{x}_i))$ can be specified in various ways. One approach is to model $\boldsymbol{\eta}_{i,t}$ using a regression with fixed effects, where the parameters are expressed as functions of observed covariates such as policyholder characteristics or temporal factors. Alternatively, the Dirichlet assumption can be incorporated as a random effect by assuming that, for each period t, the delay probability vector $\boldsymbol{p}_t$ follows a Dirichlet distribution with global parameters $\boldsymbol{\eta}_t$ . In this paper, we adopt the fixed-effects regression structure for $\boldsymbol{\eta}_{i,t}$ .

We make use of the following proposition which can also be found in Lawless (Reference Lawless1994):

Proposition 3. If $\boldsymbol{p}_t(\boldsymbol{x}_i) = (p_{t}(0;\,\boldsymbol{x}_i),\dots,p_{t}(D;\,\boldsymbol{x}_i))$ follows a Dirichlet distribution with parameters $\boldsymbol{\eta}_{i,t} = (\eta_{t,0}(\boldsymbol{x}_i),\dots,\eta_{t,D}(\boldsymbol{x}_i))$ , then the conditional probabilities $q_t(d;\,\boldsymbol{x}_i)= \frac{p_t(d;\,\boldsymbol{x}_i)}{\sum_{j=0}^{d}p_t(j;\,\boldsymbol{x}_i)}$ $(d=1,\dots,D)$ are independently distributed beta random variables with

\begin{align*}q_t(d;\,\boldsymbol{x}_i) \sim \text{Beta}(\eta_{t,d}(\boldsymbol{x}_i),\sum_{j=0}^{d-1}\eta_{t,j}(\boldsymbol{x}_i)).\end{align*}

Thus, the binomial likelihoods components in $\mathcal{L}^{(2)}$ would be replaced with beta-binomial likelihoods components, and the likelihood for our model with the Dirichlet assumption becomes:

(4.4) \begin{equation} \begin{split} \mathcal{L}^{(3)}(\boldsymbol{\Phi}_3\mid \mathcal{O}_D) & \propto P\left(\boldsymbol{N}_1^{({\text{r}},1:T)}=\boldsymbol{n}_1^{(\text{r},1:T)},\dots,\boldsymbol{N}_m^{(\text{r},1:T)}=\boldsymbol{n}_m^{(\text{r},1:T)}\right) \\ &\times \prod_{i=1}^m \prod_{t=1}^T \prod_{d=1}^{\min\{D,T-t\}} \binom{\sum_{j=0}^{d}z_{i,t,j}}{z_{i,t,d}}\frac{B(z_{i,t,d}+\eta_{t,d}(\boldsymbol{x}_i),\sum_{j=0}^{d-1}(z_{i,t,d}+\eta_{t,j}(\boldsymbol{x}_i)))}{B(\eta_{t,d}(\boldsymbol{x}_i),\sum_{j=0}^{d-1}\eta_{t,j}(\boldsymbol{x}_i))}, \end{split}\end{equation}

where B is the beta function, and $\boldsymbol{\Phi}_3$ represents the vector of the model’s parameters, including: the initial distribution $\boldsymbol{\pi}_1$ and the transition matrix $\boldsymbol{\Gamma}$ of the HMM, the regression coefficients $\boldsymbol{\theta}_j$ for the $\lambda^{(j)}$ ’s, and the parameters related to the modeling of the $\boldsymbol{\eta}$ ’s, denoted by $\boldsymbol{\delta}$ .

Remark

With $ g=1 $ (single state), we have that $ N_{i,t}^{\text{r}} $ follows a Poisson distribution, and so, this Dirichlet assumption can be seen as an extension of the model of Verbelen et al. (Reference Verbelen, Antonio, Claeskens and Crevecoeur2022). We provide the EM algorithm for this model in the next section.

The Forward-Backward Algorithm

Before sketching the EM algorithm, we first define the forward and backward probabilities as follows:

\begin{align*} & \alpha_{tj} \,:\!=\, P(\boldsymbol{N}_{1}^{(r,1:t)}=\boldsymbol{n}_{1}^{(r,1:t)},\dots,\boldsymbol{N}_{m}^{(r,1:t)}=\boldsymbol{n}_{m}^{(r,1:t)}, C_t=j), \\ & \beta_{tj} \,:\!=\, P(\boldsymbol{N}_{1}^{(r,t+1:T)}=\boldsymbol{n}_{1}^{(r,t+1:T)},\dots,\boldsymbol{N}_{m}^{(r,t+1:T)}=\boldsymbol{n}_{m}^{(r,t+1:T)}\mid C_t=j),\end{align*}

respectively. These probabilities are used to iteratively update the hidden state probabilities during the EM algorithm, and they have the following properties (see MacDonald and Zucchini, Reference MacDonald and Zucchini2016 for details):

\begin{align*} & P(C_t=j\mid \boldsymbol{N}_{1}^{(r,1:T)}=\boldsymbol{n}_{1}^{(r,1:T)},\dots,\boldsymbol{N}_{m}^{(r,1:T)}=\boldsymbol{n}_{m}^{(r,1:T)}) = \frac{\alpha_{tj}\beta_{tj}}{L_T}, \\ & P(C_{t-1}=j,C_t=k\mid \boldsymbol{N}_{1}^{(r,1:T)}=\boldsymbol{n}_{1}^{(r,1:T)},\dots,\boldsymbol{N}_{m}^{(r,1:T)}=\boldsymbol{n}_{m}^{(r,1:T)}) = \frac{\alpha_{t-1,j}\gamma_{jk}\mathbb{P}_k(\boldsymbol{n}_t^r)\beta_{tk}}{L_T},\end{align*}

where $\mathbb{P}_k(\boldsymbol{n}_t^r)=\prod_{i=1}^m \boldsymbol{P}(N_{i,t}^r=n_{i,t}^r\mid C_t=k)$ .

They are computed using the following recursive relations:

\begin{align*} & \boldsymbol{\alpha}_1 = \boldsymbol{\pi}_1\boldsymbol{P}_1(\boldsymbol{n}_1^r),\quad \boldsymbol{\alpha}_t=\boldsymbol{\alpha}_{t-1}\boldsymbol{\Gamma}\boldsymbol{P}_t(\boldsymbol{n}_t^r),\quad t=2,\dots,T,\\ & \boldsymbol{\beta}_T = \mathbf{1}, \quad \boldsymbol{\beta}_t=\boldsymbol{\Gamma}\boldsymbol{P}_{t+1}(\boldsymbol{n}_{t+1}^r) \boldsymbol{\beta}_{t+1},\quad t=T-1,T-2,\dots,1.\end{align*}

We now proceed with the details of the EM algorithm.

Complete Data Likelihood

The EM algorithm iteratively estimates the parameters $\boldsymbol{\Phi}_2$ , where the number of hidden states, g, is treated as preset. The algorithm consists of two main steps: the E-step and the M-step, which will be demonstrated in the following subsections. We first introduce the following notations:

• • $\boldsymbol{c}^{(T)}=(c_1,\dots,c_T)$ : the unobserved states of the HMM,

• • $u_{tj} \,:\!=\, \unicode{x1D7D9}\{C_t = j\}$ , and

• • $v_{tjk} \,:\!=\, \unicode{x1D7D9}\{C_{t-1}=j, C_{t}=k\}$ .

We let our complete data consist of all claims that occurred before our reserve valuation date (both reported and unreported), as well as the realizations of the states of the HMM. Thus, the complete data is given by:

\begin{align*} \mathcal{C} &= \{n_{i,t}, z_{i,t,d}, c_t \mid i=1, \dots, m;\, \, t=1, \dots, T;\, \, d=0, \dots, D\}.\end{align*}

Note that the complete data assumes that we know the total number of claims that occurred in each period for each policyholder, as well as the periods in which these claims are reported. The likelihood of the complete data is then given by

\begin{equation*}\begin{split} \mathcal{L}^{(2)}_{\text{C}}(\boldsymbol{\Phi}_2\mid \mathcal{C}) & =\pi_{1,c_1}\prod_{t=2}^T\gamma_{c_{t-1},c_t}\prod_{t=1}^T \prod_{i=1}^m \mathbb{P}_{c_t}(n_{i,t}) \\ & \times \prod_{i=1}^m \prod_{t=1}^T \prod_{d=1}^{D} \binom{\sum_{j=0}^{d}z_{i,t,j}}{z_{i,t,d}} q_t(d;\,\boldsymbol{x}_i)^{z_{i,t,d}}(1-q_t(d;\,\boldsymbol{x}_i))^{\sum_{j=0}^{d-1}z_{i,t,d}},\end{split}\end{equation*}

where $\mathbb{P}_{c_t}(n_{i,t}) = P(N_{i,t}=n_{i,t}\mid C_t = c_t)$ .

It is straightforward to see that we can write the log-likelihood of the complete data as

\begin{align*} l_{\text{C}}^{(2)}(\boldsymbol{\Phi}_2\mid \mathcal{C}) & = \sum_{j=1}^g u_{1j}\log \pi_{1j} +\sum_{j=1}^g\sum_{k=1}^g\left(\sum_{t=2}^Tv_{tjk}\right)\log \gamma_{jk} + \sum_{j=1}^g\sum_{t=1}^T \sum_{i=1}^m u_{tj}\log \mathbb{P}_{j}(n_{i,t}) \\ & + \sum_{i=1}^m \sum_{t=1}^T \sum_{d=1}^D \left[\log \binom{\sum_{j=0}^{d}z_{i,t,j}}{z_{i,t,d}} + z_{i,t,d}\log q_t(d;\,\boldsymbol{x}_i) + \left(\sum_{j=0}^{d-1}z_{i,t,d}\right)\log(1-q_t(d;\,\boldsymbol{x}_i))\right],\end{align*}

where $l_{\text{C}}^{(2)}(\boldsymbol{\Phi}_2\mid \mathcal{C})\,:\!=\,\log(\mathcal{L}_{\text{C}}^{(2)}(\boldsymbol{\Phi}_2\mid \mathcal{C}))$ . We now outline the E-step and M-step of the EM algorithm that are applied to maximize the complete-data log-likelihood.

E-step

In the k-th E-step, we compute the conditional expectation of the complete-data log-likelihood given the observed data and the current estimator $\boldsymbol{\Phi}_2^{(k-1)}$ of the model parameters $\boldsymbol{\Phi}_2$ . Note that the unobserved components of our complete data consist of the hidden states, $c^{(T)}$ , as well as the IBNR claims, $\{n_{i,t}^{\text{ibnr}}, z_{i,t,d} \mid i=1,\dots,m;\, t=1,\dots,T;\, t+d \gt T\}$ . The expectation of these components given the observed data and $\boldsymbol{\Phi}_2^{(k-1)}$ should thus be derived. The conditional expectation is given by

(5.2) \begin{equation} \begin{split} \mathbb{E}[l_{\text{C}}^{(2)}(\boldsymbol{\Phi}_2 \mid \mathcal{C}) \mid \boldsymbol{\Phi}_2^{(k-1)}, \mathcal{O}_D ] & = \sum_{j=1}^g \hat{u}_{1j}^{(k)} \log \pi_{1j} + \sum_{j=1}^g \sum_{k=1}^g \left(\sum_{t=2}^T \hat{v}_{tjk}^{(k)}\right) \log \gamma_{jk} \\ & + \sum_{j=1}^g \sum_{t=1}^T \sum_{i=1}^m \hat{u}_{tj}^{(k)} \left[\hat{n}_{i,t,j}^{(k)} \log(\lambda^{(j)}(\boldsymbol{x}_i)) - \lambda^{(j)}(\boldsymbol{x}_i)\right] \\ & + \sum_{i=1}^m \sum_{t=1}^T \sum_{d=1}^D \left[\hat{z}_{i,t,d}^{(k)} \log q_t(d;\, \boldsymbol{x}_i) + \left(\sum_{j=0}^{d-1} \hat{z}_{i,t,d}^{(k)}\right) \log(1-q_t(d;\, \boldsymbol{x}_i))\right] \\ & + \text{constant}, \end{split}\end{equation}

where the expressions for $\hat{u}_{tj}^{(k)}$ and $\hat{v}_{tjl}^{(k)}$ are given by

(5.3) \begin{equation} \begin{split} \hat{u}_{tj}^{(k)} & = \mathbb{E}\left(u_{tj} \mid \boldsymbol{n}^{(r,1:T)}, \boldsymbol{\Phi}_2^{(k-1)}\right) = \mathbb{P}\left(C_t = j \mid \boldsymbol{n}^{(r,1:T)}, \boldsymbol{\Phi}_2^{(k-1)}\right) = \frac{\alpha_{tj}^{(k-1)} \beta_{tj}^{(k-1)}}{L_T^{(k-1)}}, \end{split}\end{equation}

and

(5.4) \begin{equation} \begin{split} \hat{v}_{tjl}^{(k)} & = \mathbb{E}\left(v_{tjl} \mid \boldsymbol{n}^{(r,1:T)}, \boldsymbol{\Phi}_2^{(k-1)}\right)= \mathbb{P}\left(C_{t-1} = j, C_t = l \mid \boldsymbol{n}^{(r,1:T)}, \boldsymbol{\Phi}_2^{(k-1)}\right) \\ & = \frac{\alpha_{t-1,j}^{(k-1)} \gamma_{jl}^{(k-1)} \prod_{i=1}^{m} P(N_{i,t}^{\text{r}} = n_{i,t}^{\text{r}} \mid C_t = l, \boldsymbol{\Phi}_2^{(k-1)}) \beta_{tl}^{(k-1)}}{L_T^{(k-1)}} \end{split}\end{equation}

In these expressions, the forward probabilities, backward probabilities, and observed likelihood are obtained using the model parameters $\boldsymbol{\Phi}_2 = \boldsymbol{\Phi}_2^{(k-1)}$ from the $(k-1)$ -th iteration. The terms $\hat{u}_{tj}^{(k)}$ and $\hat{v}_{tjl}^{(k)}$ represent the updated probabilities of being in state j at time t and the transition probability from state j to state l at time t, respectively, based on the current model parameter estimates $\boldsymbol{\Phi}_2^{(k-1)}$ and the observed data.

The expressions for $\hat{n}_{i,t,j}^{(k)}$ and $\hat{z}_{i,t,d}^{(k)}$ are given by

\begin{align*} & \hat{n}_{i,t,j}^{(k)} \,:\!=\, \mathbb{E}(N_{i,t} \mid C_t = j, \mathcal{O}_D, \boldsymbol{\Phi}_2^{(k-1)}) = n_{i,t}^{\text{r}} + \unicode{x1D7D9}\{t \gt T-D\} e_{i,t} \lambda_j^{(k-1)}(\boldsymbol{x}_i) \sum_{d=T-t+1}^D p_{t}^{(k-1)}(d;\, \boldsymbol{x}_i), \\ & \hat{z}_{i,t,d}^{(k)} \,:\!=\, \mathbb{E}(Z_{i,t,d} \mid \mathcal{O}_D, \boldsymbol{\Phi}_2^{(k-1)}) = \begin{cases} z_{i,t,d}, & t \le T-D \\ \left(\sum_{j=1}^g \pi_{tj}^{(k-1)} e_{i,t} \lambda_j^{(k-1)}(\boldsymbol{x}_i) \right) p_{t}^{(k-1)}(d;\, \boldsymbol{x}_i), & t \gt T-D \end{cases}\end{align*}

where $p_t^{(k-1)}(d;\, \boldsymbol{x}_i)$ and $\lambda_j^{(k-1)}(\boldsymbol{x}_i)$ are the estimates of $p_t(d;\, \boldsymbol{x}_i)$ and $\lambda^{(j)}(\boldsymbol{x}_i)$ , respectively, computed using the parameter estimates obtained from the M-step in the $(k-1)$ -th iteration. The term $\hat{n}_{i,t,j}^{(k)}$ represents the updated expectation of the total number of claims that occurred in period t for policyholder i, given that the state of the HMM is j, while the term $\hat{z}_{i,t,d}^{(k)}$ represents the updated expectation of the portion of these claims that are reported in period $t+d$ . Note that all claims that occurred at $t \le T-D$ are reported and are part of the observed data $\mathcal{O}_D$ , and thus, the expected values of $N_{i,t}$ and $Z_{i,t,d}$ ( $t \le T-D;\, d=0, \dots, D$ ) given the observed data are equal to $n_{i,t}$ and $z_{i,t,d}$ , respectively.

By maximizing this expected complete-data log-likelihood, we can iteratively update the model parameters in the M-step of the algorithm, leading to improved parameter estimates that better fit the observed data.

M-step

The M-step in the EM algorithm involves maximizing the expected complete-data log-likelihood (Equation (5.2)) with respect to the model parameters $\boldsymbol{\Phi}_2$ , subject to the constraints

\begin{align*} & \sum_{j=1}^g\pi_{1j}=1, \text{ and } \sum_{l=1}^g\gamma_{jl}=1,\quad j=1,\dots,g.\end{align*}

For this step, we can obtain the updated estimates for $\boldsymbol{\pi}_1$ and $\boldsymbol{\Gamma}$ separately by maximizing the first two terms of Equation (5.2), respectively. The optimal values for $\pi_{1j}^{(k)}$ and $\gamma_{jl}^{(k)}$ at the kth iteration are given by

(5.5) \begin{align}\pi_{1j}^{(k)} = \hat{u}_{1j}^{(k)},\quad j=1,\dots,g, \end{align}

(5.6) \begin{align} \gamma_{jl}^{(k)} = \frac{\sum_{t=2}^T \hat{v}_{tjl}^{(k)}}{\sum_{l=1}^g\sum_{t=2}^T \hat{v}_{tjl}^{(k)}},\quad j,l=1,\dots,g. \end{align}

As for the third term:

\begin{align*}\sum_{t=1}^T \sum_{i=1}^m \hat{u}_{tj}^{(k)}\left[\hat{n}_{i,t,j}^{(k)}\log(\lambda^{(j)}(\boldsymbol{x}_i)) - \lambda^{(j)}(\boldsymbol{x}_i)\right],\end{align*}

it is easy to see that for each j, this is a weighted Poisson log-likelihood, with the exception that $\hat{n}_{i,t,j}^{(k)}$ can be non-integer. Yet, we can still obtain the updated regression parameters for each $\lambda^{(j)}$ using quasi-likelihood methods available through standard packages in R. Finally, for the term:

\begin{align*}\sum_{i=1}^m \sum_{t=1}^T \left[\hat{z}_{i,t,d}^{(k)}\log q_t(d;\,\boldsymbol{x}_i) + \left(\sum_{j=0}^{d-1}\hat{z}_{i,t,d}^{(k)}\right)\log(1-q_t(d;\,\boldsymbol{x}_i))\right],\end{align*}

it is easy to see that for each d, this is a binomial log-likelihood, again with the exception that $\hat{z}_{i,t,d}^{(k)}$ can be non-integer. For each d, we model $q_t(d;\,\boldsymbol{x}_i)$ as a regression function of the policy-holder risk attributes, $\boldsymbol{x}_i$ , and time-dependent covariates (e.g. the month of the period t). The binomial log-likelihood can then be optimized using quasi-likelihood methods available through standard packages in R.

The model parameters are then given by $\boldsymbol{\Phi}_2 = \{\boldsymbol{\pi}_1,\boldsymbol{\Gamma},\boldsymbol{\theta}_1,\dots,\boldsymbol{\theta}_g,\boldsymbol{\delta}_1,\dots,\boldsymbol{\delta}_D\}$ , where $\boldsymbol{\theta}_j$ is the vector of regression parameters for $\lambda^{(j)}(\boldsymbol{x})$ , and $\boldsymbol{\delta}_d$ is the vector of regression parameters for $q_t(d;\,\boldsymbol{x}).$

Convergence

We continue iterating the E-step and M-step until the relative distance between the estimates from consecutive iterations falls below a predetermined threshold. The relative distance is defined as

\begin{align*} d(\boldsymbol{\Phi}_2^{(k-1)},\boldsymbol{\Phi}_2^{(k)}) & = \sum_{j=1}^g \left| \frac{\pi_{1j}^{(k-1)}-\pi_{1j}^{(k)}}{\pi_{1j}^{(k-1)}} \right| + \sum_{j=1}^g \sum_{l=1}^g\left|\frac{\gamma_{jl}^{(k-1)}-\gamma_{jl}^{(k)}}{\gamma_{jl}^{(k-1)}} \right| + \sum_{j=1}^g \sum_{s=1}^{||\boldsymbol{\theta}_j||_0}\left|\frac{\theta_{js}^{(k-1)}-\theta_{js}^{(k)}}{\theta_{js}^{(k-1)}}\right|\\&+\sum_{d=1}^D \sum_{s=1}^{||\boldsymbol{\delta}_j||_0}\left|\frac{\delta_{js}^{(k-1)}-\delta_{js}^{(k)}}{\delta_{js}^{(k-1)}}\right|.\end{align*}

The EM algorithm converges when the relative distance is below the threshold.

Initialization of Model Parameters

To effectively apply the EM algorithm, appropriate initialization of model parameters is crucial. Note that we have complete data (all claims are reported) for periods $t \le T - D$ , so we can fit the model for such data to obtain parameter estimates as initialization

Model Selection

Finally, we need to decide on the number of states g. Typically, a larger g would provide a better fit, but it could also lead to overfitting. We usually select the model with the lowest chosen information criterion (e.g., AIC or BIC). For fast and efficient computation, we can employ a backward fitting strategy similar to Badescu et al. (Reference Badescu, Chen, Lin and Tang2019). We start by fitting the model with a large g, and then iteratively delete a state until the chosen information criterion stops decreasing. In this strategy, we use the parameter estimates obtained from fitting with g states to initialize the algorithm for $ g-1 $ states, with the initial values for $ \boldsymbol{\pi}_1 $ and $ \boldsymbol{\Gamma} $ being normalized estimates excluding the deleted state. This approach should make the estimation process very efficient.

Under-/Overflow Problems

As noted by MacDonald and Zucchini (Reference MacDonald and Zucchini2016), the computation of forward and backward probabilities, along with the subsequent calculations in the E- and M- steps, is susceptible to under- or overflow issues. MacDonald and Zucchini (Reference MacDonald and Zucchini2016) recommend scaling the forward and backward probability computations by taking the logarithms of these quantities and provides code solutions for this purpose. However, for our model, additional considerations are needed.

The computation of forward and backward probabilities involves matrix multiplication of the matrices $\boldsymbol{P}_t(\boldsymbol{n}_t^r)$ . Unlike normal HMMs fitted to single time series, where these matrices are diagonal matrices containing the probability of the count process obtaining its observed value given the different states, our model’s probabilities are the product of the observed count process for each individual. For our model, this matrix is defined as

\begin{align*} \boldsymbol{P}_t(\boldsymbol{n}_t^r) = \text{diag}\left\{\prod_{i=1}^m \boldsymbol{P}(N_{i,t}^r=n_{i,t}^r\mid C_t=1),\dots,\prod_{i=1}^m \boldsymbol{P}(N_{i,t}^r=n_{i,t}^r\mid C_t=g)\right\}.\end{align*}

It becomes evident that the elements of the diagonal matrix are extremely small, potentially causing underflow problems in computing forward and backward probabilities, as well as subsequent calculations, even after scaling. To address this issue, we implement the LogSumExp algorithm (e.g., see Blanchard et al., Reference Blanchard, Higham and Higham2021), which is effective in mitigating numerical instability arising from the small probabilities.

Complete Data Likelihood

We now aim to find the parameter estimates that maximize $\mathcal{L}^{(3)}$ using the EM algorithm. As we illustrated, the Dirichlet-multinomial model differs from the multinomial model by the assumption that the probability vector $\boldsymbol{p}_{t}(\boldsymbol{x}_i)$ is random. To this end, we let our complete data be given by

\begin{align*} \tilde{\mathcal{C}} &= \{n_{i,t}, z_{i,t,d}, p_{i,t,d}, c_t \mid t=1,\dots,T;\, d=0,\dots,D\},\end{align*}

where $p_{i,t,d}$ is the realization of $p_d(t;\, \boldsymbol{x}_i)$ . That is, we now assume we know the actual reporting probability vector $\boldsymbol{p}_t(\boldsymbol{x}_i)$ . The likelihood of the complete data is then given by:

\begin{equation*}\begin{split} \mathcal{L}^{(3)}_{\text{C}}(\boldsymbol{\Phi}_3\mid \tilde{\mathcal{C}}) & =\pi_{1,c_1}\prod_{t=2}^T\gamma_{c_{t-1},c_t}\prod_{t=1}^T \prod_{i=1}^m \mathbb{P}_{c_t}(n_{i,t}) \\ & \times \prod_{i=1}^m \prod_{t=1}^T \left(\underbrace{\frac{n_{i,t}!}{\prod_{d=0}^D z_{i,t,d}!}\prod_{d=1}^{D} p_{i,t,d}^{z_{i,t,d}}}_{\text{Multinomial Distribution}} \times \underbrace{\frac{1}{B(\boldsymbol{\eta}_{i,t})}\prod_{d=0}^D p_{i,t,d}^{\eta_{t,d}(\boldsymbol{x}_i)-1}}_{\text{Dirichlet Distribution}}\right),\end{split}\end{equation*}

where $\mathbb{P}_{c_t}(n_{i,t}) = P(N_{i,t}=n_{i,t}\mid C_t = c_t)$ . Note that the “Multinomial Component” in the complete data likelihood is constant since it does not depend on the model’s parameters. The log-likelihood of the complete data is then given by

\begin{align*} l_{\text{C}}^{(3)}(\boldsymbol{\Phi}_3\mid \tilde{\mathcal{C}}) & = \sum_{j=1}^g u_{1j}\log \pi_{1j} +\sum_{j=1}^g\sum_{k=1}^g\left(\sum_{t=2}^Tv_{tjk}\right)\log \gamma_{jk} + \sum_{j=1}^g\sum_{t=1}^T \sum_{i=1}^m u_{tj}\log \mathbb{P}_{j}(n_{i,t}) \\ & + \sum_{i=1}^m \sum_{t=1}^T -\log B(\boldsymbol{\eta}_{i,t}) \sum_{d=1}^D (\eta_{t,d}(\boldsymbol{x}_i)-1)\log p_{i,t,d}+ \text{constant},\end{align*}

where $l_{\text{C}}^{(3)}(\boldsymbol{\Phi}_2\mid \mathcal{C})\,:\!=\,\log(\mathcal{L}_{\text{C}}^{(3)}(\boldsymbol{\Phi}_2\mid \mathcal{C}))$ .

E-Step

In the k-th iteration of the E-step, we take the expectation of $l_{\text{C}}^{(3)}(\boldsymbol{\Phi}_3\mid \tilde{\mathcal{C}})$ with respect to $\boldsymbol{\Phi}_3^{(k-1)}$ and $\mathcal{O}_D$ . We have:

(5.7) \begin{equation} \begin{split} \mathbb{E}\left[l_{\text{C}}^{(3)}(\boldsymbol{\Phi}_3 \mid \tilde{\mathcal{C}}) \mid \boldsymbol{\Phi}_3^{(k-1)}, \mathcal{O}_D \right] &= \sum_{j=1}^g \hat{u}_{1j}^{(k)} \log \pi_{1j} + \sum_{j=1}^g \sum_{l=1}^g \left(\sum_{t=2}^T \hat{v}_{tjl}^{(k)}\right) \log \gamma_{jl} \\ & \quad + \sum_{j=1}^g \sum_{t=1}^T \sum_{i=1}^m \hat{u}_{tj}^{(k)} \left[\hat{n}_{i,t,j}^{(k)} \log (\lambda^{(j)}(\boldsymbol{x}_i)) - \lambda^{(j)}(\boldsymbol{x}_i) \right] \\ & \quad + \sum_{i=1}^m \sum_{t=1}^T - \log B(\boldsymbol{\eta}_{i,t}) \sum_{d=1}^D (\eta_{t,d}(\boldsymbol{x}_i) - 1) \widehat{\log p_{i,t,d}}^{(k)} \\ & \quad + \text{constant}, \end{split}\end{equation}

where $\hat{u}_{tj}^{(k)}$ and $\hat{v}_{tjl}^{(k)}$ are given by Equations (5.3) and (5.4), respectively, and

\begin{align*}\hat{n}_{i,t,j}^{(k)} &\,:\!=\, \mathbb{E}\left[N_{i,t} \mid C_t = j, \mathcal{O}_D, \boldsymbol{\Phi}_3^{(k-1)}\right] \\&= n_{i,t}^{\text{r}} + \unicode{x1D7D9}\{t \gt T - D\} e_{i,t} \lambda_j^{(k-1)}(\boldsymbol{x}_i) \mathbb{E}\left[\sum_{d=T-t+1}^D p_{t}(d;\, \boldsymbol{x}_i) \mid C_t = j, \mathcal{O}_D, \boldsymbol{\Phi}_3^{(k-1)}\right], \\\widehat{\log p_{i,t,d}}^{(k)} &\,:\!=\, \mathbb{E}\left[\log p_{t}(d;\, \boldsymbol{x}_i) \mid \mathcal{O}_D, \boldsymbol{\Phi}_3^{(k-1)}\right],\end{align*}

where $\lambda_j^{(k-1)}(\boldsymbol{x}_i)$ is the estimate of $\lambda^{(j)}(\boldsymbol{x}_i)$ computed using the parameter estimates obtained from the M-step in the $(k-1)$ -th iteration. We do not have closed-form analytical expressions for $\hat{n}_{i,t,j}^{(k)}$ ( $t \gt T-D$ ) and $\widehat{\log p_{i,t,d}}^{(k)}$ , and thus, numerical methods are needed to compute these expectations. We can obtain Monte Carlo estimates of these expectations using samples from the conditional distribution $\boldsymbol{p}_t(\boldsymbol{x}_i) = (p_t(0;\, \boldsymbol{x}_i), \dots, p_t(D;\, \boldsymbol{x}_i))$ , using the parameters $\boldsymbol{\Phi}_3$ . This, however, adds an extra layer of challenge because there is no guarantee that the likelihood will increase with each iteration, unlike in the case of the normal EM algorithm, as we are replacing the actual expectation with a numerical approximation.

Remark

The computation of the forward probabilities, backward probabilities, the likelihood, and consequently, $\hat{u}_{tj}$ and $\hat{v}_{tjk}$ in the E-step differs slightly for the two models considered. The differences arise when computing the probability $ P(N_{i,t}^{\text{r}} = n_{i,t}^{\text{r}} \mid \mathcal{O}_D, \boldsymbol{\Phi}^{(k-1)}) $ . As previously discussed, this probability depends on the reporting delay. For the discrete-time Multinomial model, the estimation of the probability vector $ \boldsymbol{p}_t(\boldsymbol{x}_i) $ changes after each iteration of the algorithm, and so we use $ \boldsymbol{p}_t(\boldsymbol{x}_i) $ estimated from the $ k-1 $ -th iteration in our computation in the k -th iteration. As for the Dirichlet-multinomial model, $ \boldsymbol{p}_t(\boldsymbol{x}_i) $ is a Dirichlet random vector, and thus we compute

\begin{align*} & P(N_{i,t}^{\text{r}} = n_{i,t}^{\text{r}} \mid \mathcal{O}_D, \boldsymbol{\Phi}_3^{(k-1)}) = \\ & \int_{\boldsymbol{p}} P(N_{i,t}^{\text{r}} = n_{i,t}^{\text{r}} \mid \boldsymbol{p}_t(\boldsymbol{x}_i) = \boldsymbol{p}, \mathcal{O}_D, \boldsymbol{\Phi}_3^{(k-1)}) \times P(\boldsymbol{p}_t(\boldsymbol{x}_i) = \boldsymbol{p} \mid \mathcal{O}_D, \boldsymbol{\Phi}_3^{(k-1)}) \, d\boldsymbol{p}.\end{align*}

Distribution of $\boldsymbol{p}_t \mid \mathcal{O}_D$

In the case where $ t \le T - D $ , we have:

\begin{align*}\boldsymbol{p}_t(\boldsymbol{x}_i) \mid \mathcal{O}_D \sim \text{Dirichlet}(\eta_{t,0}(\boldsymbol{x}_i) + z_{i,t,0}, \dots, \eta_{t,D}(\boldsymbol{x}_i) + z_{i,t,D}),\end{align*}

since the Dirichlet distribution is a conjugate prior of the multinomial distribution. Challenges arise in the case where $ t \gt T - D $ , as we do not observe $ z_{i,t,d} $ for $ d \gt T - t $ , and thus we lose the nice conjugate property. Therefore, we need to derive the conditional distribution of $\boldsymbol{p}_t$ for $ t \gt T - D $ .

The probability of $\boldsymbol{p}_t(\boldsymbol{x}_i)$ equal $\boldsymbol{p}=(p_0,\dots,p_D)$ given the observed data is given by

\begin{align*} P(\boldsymbol{p}_t(\boldsymbol{x}_i) = \boldsymbol{p}\mid \mathcal{O}_D) & \propto P(\mathcal{O}_D\mid \boldsymbol{p}_t(\boldsymbol{x}_i)=\boldsymbol{p}) \times P(\boldsymbol{p}_t(\boldsymbol{x}_i) = \boldsymbol{p}),\end{align*}

where

\begin{align*} & P(\mathcal{O}_D\mid \boldsymbol{p}_t(\boldsymbol{x}_i)=\boldsymbol{p}) \propto P(N_{i,t}^{\text{r}}=n_{i,t}^{\text{r}}\mid \boldsymbol{p}_t(\boldsymbol{x}_i)=\boldsymbol{p})\times\\&P(Z_{i,t,0}=z_{i,t,0},\dots,Z_{i,t,T-t}=z_{i,t,T-t}\mid N_{i,t}^{\text{r}}=n_{i,t}^{\text{r}},\boldsymbol{p}_t(\boldsymbol{x}_i)=\boldsymbol{p}).\end{align*}

Thus,

\begin{align*} & P(\boldsymbol{p}_t(\boldsymbol{x}_i) = \boldsymbol{p}\mid \mathcal{O}_D) \propto P(\boldsymbol{p}_t(\boldsymbol{x}_i) = \boldsymbol{p}) \times P(N_{i,t}^{\text{r}}=n_{i,t}^{\text{r}}\mid \boldsymbol{p}_t(\boldsymbol{x}_i)=\boldsymbol{p}) \times \\ & P(Z_{i,t,0}=z_{i,t,0},\dots,Z_{i,t,T-t}=z_{i,t,T-t}\mid N_{i,t}^{\text{r}}=n_{i,t}^{\text{r}},\boldsymbol{p}_t(\boldsymbol{x}_i)=\boldsymbol{p}) \\ & \propto \prod_{d=0}^D (p_d)^{\eta_{t,d}(\boldsymbol{x}_i)-1} \times \left(\sum_{j=1}^g \pi_{tj}\frac{e^{-\lambda^{(j)}(\boldsymbol{x}_i)p^{\text{r}}}(\lambda^{(j)}(\boldsymbol{x}_i)p^{\text{r}})^{n_{i,t}^{\text{r}}}}{n_{i,t}^{\text{r}}!}\right)\times\frac{n_{i,t}^{\text{r}}!}{\prod_{d=0}^{T-t}z_{i,t,d}!} \prod_{d=0}^{T-t}\left(\frac{p_d}{p^{\text{r}}}\right)^{z_{i,t,d}} \\ & \propto \underbrace{\left(\prod_{d=0}^{T-t} (p_d)^{z_{i,t,d}+\eta_{t,d}(\boldsymbol{x}_i)-1}\right) \times \left( \prod_{d=T-t+1}^D (p_d)^{\eta_{t,d}(\boldsymbol{x}_i)-1}\right)}_{h} \times \underbrace{\sum_{j=1}^g \pi_{tj}\left(\lambda^{(j)}(\boldsymbol{x}_i)\right)^{n_{i,t}^{\text{r}}}e^{-\lambda^{(j)}(\boldsymbol{x}_i)p^{\text{r}}}}_{g},\end{align*}

where $p^{{\text{r}}}=\sum_{d=0}^{T-t}p_d$ . Note that the conditional distribution in this case is not a standard distribution, and thus, we cannot sample from it directly.

Rejection Sampling

Samples from the conditional distribution of $\boldsymbol{p}_t(\boldsymbol{x}_i)$ can be selected by multivariate rejection sampling. We can write $P(\boldsymbol{p}_t(\boldsymbol{x}_i)=p\mid \mathcal{O}_D)=ah(\boldsymbol{p})g(\boldsymbol{p}).$ It is straightforward to see that h is the density function of Dirichlet distribution with parameters $(z_{i,t,0}+\eta_{t,0}(\boldsymbol{x}_i),\dots,z_{i,t,T-t}+\eta_{t,T-t}(\boldsymbol{x}_i),\eta_{t,T-t+1}(\boldsymbol{x}_i),\dots,\eta_{t,D}(\boldsymbol{x}_i)),$ and thus, we can easily sample from it. We perform the rejection sampling as follows:

1. • Step 1: sample $\boldsymbol{p}$ from h, and independently, sample u from a Uniform(0,1) distribution.

2. • Step 2: If $u \le \frac{g(\boldsymbol{p})}{\sup_{\boldsymbol{p}}\{g(\boldsymbol{p})\}},$ then accept $\boldsymbol{p}$ . If not, go to Step 1.

Note that $\sup_{\boldsymbol{p}}\{g(\boldsymbol{p})\}\le \sum_{j=1}^g \pi_{tj}(\lambda^{(j)}(\boldsymbol{x}_i))^{n_{i,t}^{\text{r}}}$ . By obtaining a sample for $\boldsymbol{p}_t(\boldsymbol{x}_i)$ given $\mathcal{O}_D$ and $\boldsymbol{\Phi}_3^{(k-1)},$ we can estimate $\hat{n}_{i,t,j}^{(k)}$ and $\widehat{\log p_{i,t,d}}^{(k)}$ using Monte Carlo methods.

Remark: The formulas presented above in the derivation of the conditional distribution of $\boldsymbol{p}_t$ implicitly assume an exposure of $e_{t,i}=1$ . To account for exposure, simply multiply the Poisson intensity by the exposure factor.

M-Step

We now maximize Equation (5.7) from the E-step after replacing $\hat{n}_{i,t,j}^{(k)}$ ( $t \gt T-D$ ) and $\widehat{\log p_{i,t,d}}^{(k)}$ with their Monte Carlo estimates. Again, the updated estimates of $\boldsymbol{\pi}_1$ and $\boldsymbol{\Gamma}$ are given by Equations (5.5) and (5.6), respectively. Moreover, similar to the Multinomial model case, given each j, it is straightforward to see that the third term is a weighted Poisson log-likelihood, and thus can be fit using standard packages in R. Finally, for the fourth and final term, it is straightforward to see that this is a Dirichlet log-likelihood. As previously illustrated, there are several flexible ways to model $\boldsymbol{\eta}$ . For our data analysis, we model $\boldsymbol{\eta}$ as a regression function of $\boldsymbol{x}_i$ and time-dependent covariates. We obtain the parameter estimates for this Dirichlet regression using the DirichletReg package in R (Maier, Reference Maier2014), which rescales the probability estimates from the E-step to ensure that the constraints are satisfied.

Convergence Criterion, Model Initialization, and Model Selection

We deviate from the standard EM algorithm and use the MCEM algorithm, as we compute the expectations in the E-step using Monte Carlo methods. In the MCEM algorithm, it is not guaranteed that the likelihood will increase after each iteration, unlike the standard EM algorithm (see Booth and Hobert, Reference Booth and Hobert1999). Thus, special care should be taken. Many methods have been proposed for approximating the EM algorithm (see Ruth, Reference Ruth2024 for a review). From a practical point of view, we use the same convergence criterion as we have for the previous model (the relative distance criterion), as this would have a negligible effect on the parameter estimates even though the likelihood can fluctuate.

For model initialization, we fit the model to our complete data up to period $ T - D $ (as we did in the Multinomial case) and use the parameter estimates as the initialization for the algorithm. Finally, we choose the number of states g using a backward strategy as explained before.

Data Description

Our dataset contains information on Physical Damage policies issued by the company between January 2009 and December 2017. During these years, a total of 293,709 one-year policies were issued, resulting in 106,187 reported claims. For each policy, the dataset includes the policy start and end dates, as well as various vehicle and policyholder attributes. These attributes include the age of the driver (18–90), the age of the car (0–26), the class of the car (A, B, or C), the fuel type (Gasoline or Diesel), whether the contract is a renewal or issued to a new policyholder, and the region (five regions). For each claim, the dataset records the day on which the claim occurred, the day on which the claim was reported, and the progression of claim payments. In this study, we focus on estimating the IBNR claim count, which requires us to concentrate on the time of occurrence and reporting delay of each claim. Therefore, we will not consider the payments in our analysis.

Reporting Delay

From the 106,197 claims reported between January 2009 and December 2017, only two were reported more than two years after the occurrence date. This suggests that nearly all claims occurring before 2016 have been reported. Table 1 displays the distribution of reporting delays in months for these claims. A delay of zero months indicates a claim was reported in the same month it occurred, while a delay of one month means it was reported the following month, and so on. The majority of claims (81.71%) were reported in the same month they occurred. Over 99% of claims were reported within five months, and only 0.1% were reported after ten months. This indicates that most claims are reported promptly after occurrence, with only a small fraction experiencing delays beyond ten months.

Table 1. Reporting delay distribution for claims that occurred before 2016.

Frequency

Figure 1 illustrates the trends in exposure (left) and claim count per exposure (right) from January 2009 to December 2015. The exposure, which represents the number of policies at risk, exhibited a steady increase from January 2009 until December 2012. Following this period, it decreased consistently until January 2015, after which it began to rise again. In contrast, the claim count per exposure displayed significant fluctuations over the same period. It reached a minimum of 0.023 claims per exposure in December 2014 and peaked at 0.045 claims per exposure in May 2013. This variation in claim frequency highlights the need for a sophisticated model capable of capturing such dynamic changes. The HMM employed in our model is well-suited for this purpose, as it can account for the underlying state changes that drive fluctuations in claim frequency over time.

Figure 1. Monthly exposure (left) and claim count per exposure (right) between January 2009 and December 2015.

Specifications for the Fitted Models

• • The Continuous-time Model (CM): As previously explained, our approach begins by fitting the distribution for the reporting delay, which is then utilized to fit our model, and consequently estimate the IBNR claim count. A notable challenge in this process is the right truncation of the reporting delay, which complicates the fitting of the distribution. We provide the continuous-time model with an advantage over the discrete-time models as we fit the reporting delay using all claims that occurred before the reserve valuation date $\tau$ , including those that have not yet been reported. This overcomes the issue of right truncation, allowing for a more accurate fit.

While this method is not feasible in practice – since unreported claims cannot be known beforehand – it serves to highlight the strengths of the discrete-time models, which we will demonstrate in the results section. In a practical scenario, one could use a “back censoring” technique, where a time threshold within the training period is identified. Claims occurring before this threshold would have been reported by the valuation date, allowing the model to be fit without concern for right truncation.

In our analysis, the log-logistic proportional hazard regression performed well in fitting the reporting delay data, as evidenced by the Cox-Snell residuals plot in Figure 2. The covariates used in the regression include the the vehicle and policyholder risk attributes introduced in the data description (Section 7.1), the month in which the claim occurred and the day of occurrence.

Figure 2. Cox-Snell residual plot for assessing the fit of the reporting delay model

• • The Multinomial Model (MM): In Section 5.1, we outlined the EM algorithm required to fit the multinomial model, where we model the binomial probabilities $ q_t(d;\,\boldsymbol{x}_i) $ for each $ d \in \{1,\dots,D\} $ as a regression function of the vehicle and policyholder’s risk attributes $\boldsymbol{x}_i$ and time-dependent covariates. The vehicle and policyholder’s covariates are the vehicle and policyholder risk attributes introduced in the data description (Section 7.1), while the time-dependent covariates include the month of occurrence and the weekday on which the period ends, both of which were found to be significant in our regression analysis. Given that $ q_t(d;\,\boldsymbol{x}_i) $ is small for $ d \gt 1 $ , we use a complementary log-log link function to model it. For $ d = 1 $ , we use a standard log link function, as is typical for binomial GLMs. We fit the binomial component of the likelihood using the standard glm function in R.

• • The Dirichlet-Multinomial Model (DM): The fitting process for the Dirichlet-Multinomial model is thoroughly detailed in Section 5.2. We model the Dirichlet prarameters $\boldsymbol{\eta}$ as a regression function of the same covariates used for $q_t(d;\,\boldsymbol{x}_i)$ in the Multinomial Model.

• • The Chain Ladder (CL): We estimate the IBNR claim count using the CL method on data aggregated monthly, serving as the baseline for comparing our models. As with the discrete-time models, we assume that claims are reported no later than 9 months after the month of occurrence.

Remarks

• • For the continuous-time model, we do not assume a maximum reporting delay, so by following the steps in Section 6, we can estimate the full number of IBNR claims, unlike the discrete-time models. However, for a fair comparison with the discrete-time models, we use the estimated parameters after fitting the continuous-time model to estimate only the IBNR claims with a maximum delay period of 9 months after the month of claim occurrence. This is achieved by adjusting Step 5 so that the probability factor $\left(1 - \int_{d_{l-1}}^{d_l} F_{U\mid t,\boldsymbol{x}_i}(\tau-t) \, dt\right)$ is replaced with:

\begin{align*}\int_{d_{l-1}}^{d_l} P(U+t \in [\tau, t_{\text{max}}]\mid t,\boldsymbol{x}_i) \, dt,\end{align*}

• where t is the day of claim occurrence, U is the reporting delay random variable, and $t_{\text{max}}$ is the date by which the claim should be reported if we assume a maximum delay of 9 months after the month of claim occurrence. For example, if the claim occurrence date is in January, then $t_{\text{max}}$ will be the end of October for the same year. This adjustment ensures that only the first 9 months of reporting delay are considered in the estimation, aligning the continuous-time model’s output with the assumptions made in the discrete-time models.

• • The covariates used for the state-dependent intensities, $\lambda^{(j)}(\!\cdot\!)$ , are the vehicle and policyholder risk attributes introduced in the data description (Section 7.1). More specifically, these attributes include the age of the driver, the age of the car, the car class, whether the contract is a renewal or issued to a new policyholder, and the policyholder’s region.

Performance Comparison of the Chain Ladder Method and HMM-Based Models

The classical CL method performs worse than our proposed HMM-based models. The CL method shows a mean and median APE of around 16%, with high variability in the error, as indicated by both the boxplot and the standard deviation values presented in Table 2. This highlights the limitations of the CL method compared to more advanced approaches.

Performance of the Proposed Models

Our proposed models – the continuous-time model (CM), the multinomial model (MM), and the Dirichlet-mMultinomial model (DM) – perform much better than the CL method, demonstrating their effectiveness in estimating IBNR claim counts. The discrete-time models (MM and DM) generally outperform the CM, as evidenced by their lower median and mean APE and lower variability (see Figure 3) across all values of g (2, 3, and 4). The CM model shows a median absolute error close to that of the MM and DM models, but the mean absolute error is significantly larger due to its high variability, as reflected in the boxplot.

In general, the CM model tends to produce good estimates of IBNR claim counts when the integral

\begin{align*}\int_{d_{l-1}}^{d_l} F_{U\mid t,\boldsymbol{x}_i}(\tau-t) \, dt,\end{align*}

representing the probability that a claim occurring within the period $(d_{l-1}, d_l]$ is reported before the valuation date $\tau$ , aligns closely with the actual percentage of such claims. However, when this integral deviates significantly from the actual percentage, the CM model’s estimates deteriorate. For example, at the end of January 2016, the average value of the integral

\begin{align*}\int_{d_{T-1}}^{d_T} F_{U\mid t,\boldsymbol{x}_i}(\tau-t) \, dt,\end{align*}

where $d_T =\tau$ (which represents the probability that a claim that occurred in January 2016 was reported before the end of the month), across policyholders was 0.81095, while the actual percentage was 0.7465753 (i.e., 75% of such claims were reported by the end of the month) – a difference that diverged from the usual pattern. This discrepancy, which was not captured by our model for reporting delay despite using all available data, significantly impacted the IBNR claim count estimate; the underestimated value of

\begin{align*}1-\int_{d_{T-1}}^{d_T} F_{U\mid t,\boldsymbol{x}_i}(\tau-t) \, dt\end{align*}

led to an underestimated IBNR claim count.

Effect of Increasing States on Model Performance

As the number of states g increases, the performance of the discrete-time models (MM and DM) improves, suggesting that these models benefit from more states in capturing the underlying dynamics of the data. However, this trend does not hold for the continuous-time model (CM). The CM model’s performance remains inconsistent due to its sensitivity to the fitting of reporting delays and the two-step maximization of the likelihood process, as explained in the previous analysis. Therefore, while adding states improves the accuracy of discrete-time models, it does not necessarily translate into better performance for the continuous-time model.

BIC for the Fitted Models

We report the BIC for the fitted models at the first reserve valuation date. The BIC values for the Continuous-time, Multinomial, and Dirichlet-Multinomial models with $g = 2, 3, 4$ are shown in Table 3.

Table 3. BIC values for three fitted HMM-based models with $g=2,3,4$ at the first reserve valuation date.

The results indicate that the BIC is lowest for $g = 4$ across all three models, suggesting that $g = 4$ provides the best fit. However, the BIC values for $g = 3$ and $g = 4$ are very close, indicating that the improvement in fit from $g = 3$ to $g = 4$ is relatively small.

Effect of Incorporating the HMM

Table 2 also shows summary statistics for the APE in the case $ g = 1 $ , which translates to the Poisson process. It is clear that the performance of both the discrete-time models improves by incorporating the HMM. For both models, the mean APE decreases significantly as the number of states increases. This indicates that incorporating the HMM structure, with its ability to account for different state transitions, provides better accuracy in the IBNR claim count estimates. The improvement is seen not only in the mean APE but also in the reduction of the standard deviation of the APE, reflecting the model’s increased stability and precision with more states.

In contrast, the continuous-time CM model shows minimal change when incorporating the HMM. The APE remains largely consistent, indicating that the continuous-time model is less affected by the number of states. The same reasoning we outlined in the previous subsection, when talking about the effect of increasing states on model performance, applies here.

Effect of Including Covariates

To evaluate the effect of including the vehicle and policyholder’s risk characteristics, we fitted our models without incorporating any covariates in the framework. The results indicate that, across all numbers of hidden states g, the APE ranged between 16 and 17% for all three models (the errors are all very similar so we do not include a table). This is comparable to the results obtained using the CL method. These findings highlight the significant role of covariates in improving model performance, as their inclusion enables the models to account for variations in claim counts arising from individual-level characteristics, which cannot be captured in a purely aggregate framework.

Comparison with the IPW Method (Vanegas et al., Reference Vanegas, Badescu and Lin2024)

To provide additional context for the performance of our proposed models, we compared them to the inverse probability weighting (IPW) method proposed in Vanegas et al. (Reference Vanegas, Badescu and Lin2024). This method was implemented using the same reporting delay regression structure as in our continuous-time model. The mean APE for the IPW method was 15%, only 2% less than the APE obtained by the CL.

It is important to note that our proposed models fundamentally differ from the IPW method in their approach. While the IPW focuses on estimating development factors through reporting delay regressions, our models incorporate the temporal dynamics of the claim arrival process. This temporal component enables our models to capture additional patterns and dependencies in the data, resulting in superior performance relative to the IPW method.

Comparison with Chain Ladder under Yearly Aggregation

Figure 4 reveals a clear seasonal pattern in the IBNR claim counts. Although Tee and Käärik (Reference Tee and Käärik2022) argue that CL estimation tends to be more accurate under finer (e.g., monthly) aggregation, it is important to also consider coarser aggregations when seasonal effects in reporting structure are present. In practice, practitioners may prefer yearly aggregation to better capture full seasonal cycles in reporting structures.

To assess the impact of aggregation choice on model performance, we compare our models against the CL method applied to yearly aggregated data. Since the CL method with yearly aggregation does not allow restriction to claims reported within a maximum delay of 9 months, we estimate the total IBNR claim count for all models to ensure a fair comparison. For our discrete-time models, which assume a maximum delay of $D = 9$ , we estimate the long-delay component (i.e., for delays $d \gt 9$ ) using the empirical average:

\begin{align*}\hat{z}_{t,d} = \frac{1}{|\mathcal{L}_{d}|} \sum_{l \in \mathcal{L}_{d}} z_{l,d}, \quad \text{where } \mathcal{L}_{d} = \{l \mid l + d \leq T\},\end{align*}

and $z_{l,d}$ is the number of claims occurring at time l and reported with delay d. This correction is applied to all censored entries with $d \gt 9$ and $t + d \gt T$ , up to the maximum observed reporting delay in the data, allowing for a complete IBNR claim count estimate.

We compute IBNR claims counts at 24 monthly valuation dates spanning January 2014 to December 2015. Given that the maximum reporting delay in the dataset is two years and data are available up to December 2017, we are able to compare these estimates with the actual observed IBNR claim counts. The mean APE for IBNR claim counts across the 24 valuation dates - calculated for the CL method with monthly aggregation, the CL method with yearly aggregation, and the three proposed HMM-based models with 4 states - is shown in Table 4.

Table 4. Mean APE for IBNR claim count estimates across models.

Figure 4. Error bar plots showing the 95% confidence intervals for simulated IBNR claim count estimates using the DM and MM models with $ g = 2, 3, 4 $ . The red points indicate when the actual IBNR claim count falls outside the interval, and the blue points indicate when it falls within the interval

The results confirm that yearly aggregation improves CL performance compared to monthly aggregation, likely due to better alignment with seasonal patterns in reporting delays. Notably, the performance of CL-Yearly is comparable to our continuous-time model. Nevertheless, both of our discrete-time HMM-based models outperform the CL benchmark under both aggregation schemes.

Analysis of Confidence Intervals for DM and MM Models

The error bar plots in Figure 4 illustrate the 95% confidence intervals for the simulated IBNR claim count estimates using the DM and MM models with $ g = 2, 3, 4 $ , while the points represent the actual IBNR claim counts. The red points indicate instances where the actual IBNR claim count falls outside the interval, and the blue points indicate instances where it falls within the interval. It is evident that the DM model’s confidence intervals are generally wider and more frequently contain the actual IBNR claim counts compared to the MM model. Specifically, for the DM model, the confidence intervals include the actual IBNR claim count in 28, 29, and 28 out of the 36 valuation dates for $ g = 2, 3, 4 $ , respectively. In contrast, the MM model’s intervals contain the actual IBNR claim counts only 17, 20, and 21 times for $ g = 2, 3, 4 $ , respectively. The increased width of the confidence intervals in the DM model reflects the introduction of variability in the reporting probability vector through the Dirichlet assumption. This added variability leads to a more appropriate distribution of the IBNR claim counts, capturing the inherent uncertainty and providing a more realistic estimate. This highlights the importance of considering such variability for better risk management, as underestimating or overestimating reserves can have significant financial implications.