2. A log-normal chain ladder model with calendar year effects

For $i\in \{1,\dots, m\}$ and $k\in \{1,\dots, n\}$ let $C_{i,k}\gt0$ denote the cumulative claims of accident year i after development year k (payments or incurred amounts). Then $C_{i,k}$ is the cumulative claims amount of the accident year i after the calendar year $i+k-1$ . We assume that the random variables $C_{i,k}$ are observable up to the calendar year p (i.e., for $i+k\le p+1$ ) for a p with $\max (m,n)\le p \lt m+n-1$ . Let ${A}\, :\!= \,\{1,\dots,m\}\times\{1,\dots,n\}$ denote the set of all pairs of accident and development years and let $\mathcal{O}\, :\!= \,\{(i,k)\in\mathcal{A}\,|\, i+k\le p+1\}$ be the subset containing the indices of the observable $C_{i,k}$ . The set of indices corresponding to calendar year q is denoted by $\mathcal{C}(q)\, :\!= \,\{(i,k)\in \mathcal{A}\,|\, i+k=q+1\}$ . The purpose of a reserving method is to predict the cumulative amounts $C_{i,k}$ with $(i,k)\in \mathcal{P}\, :\!= \,\mathcal{A}\setminus\mathcal{O}$ . We define

\begin{align*}I(k)\, :\!= \,\min\{p-k+1,m\}\quad\text{and}\quad K(i)\, :\!= \,\min\{p-i+1,n\} \end{align*}

such that $C_{1,k},\dots,C_{I(k),k}$ are the observable cumulative claims with development year k and $C_{i,1},\dots,C_{i,K(i)}$ are the observable cumulative claims from accident year i. The LDFs

\begin{align*}R_{i,k}\, :\!= \,\frac{C_{i,k}}{C_{i,k-1}} \end{align*}

are defined for $(i,k)\in \mathcal{A}^*\, :\!= \,\{(i,k)\in \mathcal{A}\,|\, i\ge 2\}$ and observable for $(i,k)\in \mathcal{O}^*\, :\!= \,\mathcal{O}\cap \mathcal{A}^*$ . The developments of calendar year q are described by the LDFs $R_{i,k}$ with indices $(i,k)\in \mathcal{C}^*(q)\, :\!= \,\mathcal{C}(q)\cap \mathcal{A}^*$ . This notation is illustrated in Figure 1.

Figure 1.

Illustration of the notation.

For simplicity, we assume that the random variables $C_{1,1},\dots,C_{m,1}$ are deterministic. Alternatively, the entire analysis could be performed conditionally, given $C_{1,1},\dots,C_{m,1}$ . Let

\begin{align*}\mathscr{B}_{i,k}\, :\!= \,\sigma(C_{i,1},\dots,C_{i,k}) \end{align*}

be the sigma algebra generated by $C_{i,1},\dots,C_{i,k}$ and let $\operatorname{P}^{X|\mathscr{B}}$ denote the conditional distribution of a random variable X given a $\sigma$ -algebra $\mathscr{B}$ .

The following two model assumptions define the LMCL model (see Riegel, Reference Riegel2026), which will be used as basis for our study. By design, the LMCL model adheres closely to Mack’s chain ladder assumptions. In particular, its variance assumptions are considerably more in line with Mack’s framework than those of Hertig’s model. Below, we will extend the model assumptions to consider calendar year effects.

Log-normal Mack chain ladder (LMCL) model.

1. LMCL1: The accident years are independent.

2. LMCL2: There are parameters $\mu_2,\dots,\mu_{n}$ and $\sigma_2,\dots, \sigma_{n}\gt0$ such that

   \begin{align*} \operatorname{P}^{R_{i,k}|\mathscr{B}_{i,k-1}} = \operatorname{LogNormal}\!(\mu_{k}-\sigma^2_k/(2C_{i,k-1}),\sigma^2_k/C_{i,k-1}). \end{align*}

With $m_k\, :\!= \,\exp\!(\mu_k)$ and $s_k\, :\!= \,m_k\sigma_k$ assumption LMCL2 implies

(1) \begin{equation}\operatorname{E}\!(R_{i,k}|\mathscr{B}_{i,k-1})=\exp\!(\mu_k)=m_k \end{equation}

and

(2) \begin{equation}\operatorname{Var}\!(R_{i,k}|\mathscr{B}_{i,k-1})=\exp\!(\mu_k)^2\cdot\left[\exp\!\left(\frac{\sigma_k^2}{C_{i,k-1}}\right)-1\right]\approx \exp\!(\mu_k)^2\cdot \frac{\sigma_k^2}{C_{i,k-1}}=\frac{s_k^2}{C_{i,k-1}}, \end{equation}

where we have used the linear approximation $\exp\!(x)-1\approx \exp\!(0)+\exp'(0)x-1=x$ in (2). Consequently, the LMCL model satisfies the chain ladder assumptions CL1 and CL2 and satisfies CL3 approximately.

Riegel (Reference Riegel2026) proposes an iteration for computing maximum likelihood estimators of the parameters $\mu_2,\dots, \mu_n$ and $\sigma_2,\dots,\sigma_n$ , and uses them to introduce a new reserving method, referred to as the the LMCL method. The paper also discusses the accuracy of the underlying approximation and presents several numerical examples demonstrating that the LMCL method is closely aligned with the chain ladder method.

In the following, we introduce the LMCLcy model, which extends the LMCL model by incorporating calendar year effects. For $(i,k)\in \mathcal{A}^*$ , we model the calendar year effects with log-normally distributed random variables $F_{i,k}$ , which scale the LDFs $R_{i,j}$ . The moments

\begin{align*}f_{i,k}\, :\!= \,\operatorname{E}\!(F_{i,k})\gt0\quad \text{and}\quad t_{i,k}^2\, :\!= \,\operatorname{Var}\!(F_{i,k})\gt0 \end{align*}

are provided as prior information. $f_{i,k}$ represents the average impact of a calendar year effect in the market, whereas $t_{i,k}$ reflects the degree to which this effect varies across portfolios. A small value of $t_{i,k}$ indicates that all portfolios are affected in a similar way, whereas a large value suggests substantial heterogeneity.

For calendar years q in which no exceptional market-wide effects on the LDFs are observed, the intention is to set $f_{i,k}=1$ for all $(i,k)\in \mathcal{C}^*(q)$ . Note that $f_{i,k}=1$ for all $(i,k)\in \mathcal{A}^*$ does not imply the absence of calendar year effects; for instance, constant calendar year inflation does not distort the LDFs in a triangle. The factors $F_{i,k}$ within a fixed calendar year q may differ in their moments; however, since they are intended to model a common calendar year effect, they are assumed to be comonotonic. Calendar year factors belonging to different calendar years are assumed to be independent.

Remark 3. In the case study in the online supplementary material, we estimate the means $f_{i,k}$ from the market-wide LDFs provided in Table 1. If LDFs are available for all (or at least many) portfolios in the market, it may also be possible to estimate the variances $t_{i,k}^2$ directly from the data. However, in practice, a pragmatic approach is typically required to specify reasonable values for $t_{i,k}$ .

In the case study, we set $f_{i,k}=1$ for all normal calendar years, that is, all years which are not systematically affected by COVID-19. In such a situation, one possible approach is to choose a parameter $\rho \ge 0$ and define

\begin{align*} t_{i,k}\, :\!= \,\rho\cdot |f_{i,k}-1|. \end{align*}

The choice $\rho=0$ corresponds to the assumption that the calendar year effects observed in the market fully apply to the individual triangle under consideration. As $\rho$ increases, less credibility is assigned to the market prior along calendar year diagonals. Note that this approach permits calendar year effects in the individual triangle only in those years, where market-wide effects were observed, that is, for calendar years q for which $f_{i,k}\neq 1$ for some $(i,k)\in \mathcal{C}^*(q)$ . Strategies for selecting the parameter $\rho$ are discussed in the online supplementary material.

An alternative method for selecting the parameters $f_{i,k}$ and $t_{i,k}^2$ , which is more closely aligned with the approach taken in Riegel (Reference Riegel2023), is provided in Appendix A.

Remark 4. At first glance, one might expect to obtain a non-informative prior as $t_{i,k}\to \infty$ . However, this is not the case because we assume log-normal distributions for the $F_{i,k}$ . As $t_{i,k}$ increases, the distribution of $F_{i,k}$ becomes increasingly skewed, with a pronounced mode that approaches zero. This behavior can lead to unexpected effects. As long as

\begin{align*} \frac{t_{i,k}}{f_{i,k}}=\operatorname{CV}\!(F_{i,k})\le 0.2, \end{align*}

the log-normal density remains close in shape to a normal density with the same first two moments, and no issues arise. In contrast, cases with large coefficients of variation (e.g., $t_{i,k}/f_{i,k}\ge 0.5$ ) require careful handling.

Remark 5. In the model assumptions, it is convenient to use the moments

\begin{align*}\varphi_{i,k}\, :\!= \,\operatorname{E}\!(\!\log(F_{i,k}))\quad \text{and}\quad \tau^2_{i,k}\, :\!= \,\operatorname{Var}\!(\!\log(F_{i,k})) \end{align*}

of the logarithmized calendar year factors $\log(F_{i,k})$ – that is, the parameters of the log-normal distributions – instead of the moments $f_{i,k}$ and $t_{i,k}^2$ . The parameters $\varphi_{i,k}$ and $\tau_{i,k}^2$ can be calculated from $f_{i,k}$ and $t_{i,k}^2$ with the formulas

\begin{align*} \tau_{i,k}^2=\log\left(t_{i,k}^2/f_{i,k}^2+1\right)\quad\text{and}\quad \varphi_{i,k}=\log(f_{i,k})-\tau_{i,k}^2/2. \end{align*}

Vice versa, we have

\begin{align*}f_{i,k}=\exp\!(\varphi_{i,k}+\tau_{i,k}^2/2)\quad\text{and}\quad t_{i,k}^2=f_{i,k}^2\cdot (\!\exp\!(\tau_{i,k}^2)-1). \end{align*}

Let $\mathscr{F}_q$ denote the sigma algebra generated by $\{F_{i,k}\,|\,(i,k)\in \mathcal{C}^*(q)\}$ . Straightforward generalization of the model assumptions LMCL1 and LMCL2 leads to the following model.

Log-normal Mack chain ladder model with calendar year effects (LMCLcy model).

1. LMCLcy1: The calendar year factors of different calendar years are independent, that is, the $\sigma$ -algebras $\mathscr{F}_2,\dots,\mathscr{F}_{m+n-1}$ are independent. For every calendar year $q=2,\dots,m+n-1$ the vector

   \begin{align*}(F_{i,k})_{(i,k)\in\mathcal{C}^*(q)} \end{align*}
   is comonotonic. Conditionally, given $\mathscr{F}\, :\!= \,\sigma(\mathscr{F}_2,\dots,\mathscr{F}_{m+n-1})$ , the accident years are independent, that is, the vectors
   \begin{align*} \{(C_{i,1},\dots, C_{i,n})\, |\, i=1,\dots,m\} \end{align*}
   are conditionally independent.

2. LMCLcy2: There are known parameters $\varphi_{i,k}, \tau_{i,k}$ , $(i,k)\in \mathcal{A}^*$ , such that

   \begin{align*} F_{i,k}\sim \operatorname{LogNormal}\!(\varphi_{i,k}, \tau^2_{i,k}) \end{align*}
   (prior distribution of the calendar year factors).

3. LMCLcy3: There are parameters $\mu_2,\dots,\mu_{n}$ and $\sigma_2,\dots, \sigma_{n}\gt0$ such that

   \begin{align*} \operatorname{P}^{R_{i,k}|\sigma(\mathscr{B}_{i,k-1}, \mathscr{F})} = \operatorname{LogNormal}\!(\mu_{k}-\sigma^2_k/(2C_{i,k-1}) + \log(F_{i,k}),\sigma^2_k/C_{i,k-1}) \end{align*}
   or equivalently
   \begin{align*} \operatorname{P}^{R_{i,k}/F_{i,k}|\sigma(\mathscr{B}_{i,k-1}, \mathscr{F})} = \operatorname{LogNormal}\!(\mu_{k}-\sigma^2_k/(2C_{i,k-1}) ,\sigma^2_k/C_{i,k-1}). \end{align*}

For the remainder of this paper, we work in this LMCLcy model.

Remark 6. As a consequence of the model assumptions LMCLcy1 and LMCLcy2, there are independent standard normally distributed random variables $X_q$ , $q\in \{2,\dots,m+n-1\}$ , such that

\begin{align*} \log(F_{i,k})=\varphi_{i,k}+\tau_{i,k}X_{i+k-1} \end{align*}

for all $(i,k)\in \mathcal{A}^*$ . We refer to the $X_q$ as residual calendar year effects. Note that the $\sigma$ -algebra $\mathscr{F}_q$ is generated by $X_q$ .

4. Prediction error

In this section, we derive an estimator for the MSEP of the LMCLcy method. As mentioned above, a more robust version of the method will be introduced in Section 5. There, we will also derive formulas for the MSEP of this robustified method based on the formulas that are presented in the current section.

Let $w\, :\!= \,(w_1,\dots,w_m)\in [0,1]^m$ and

\begin{align*}R_w\, :\!= \,\sum_{i=1}^m w_i\big(C_{i,n}-C_{i,K(i)}\big), \quad\text{and}\quad \widehat{R}_w^{\,(\infty)}\, :\!= \,\sum_{i=1}^m w_i \big(\widehat{C}_{i,n}^{\,(\infty)}-C_{i,K(i)}\big). \end{align*}

If $w_{i}=1$ for all $i\in \{1,\dots,m\}$ , then $R_w$ is the total reserve for all accident years. If there is an $i\in \{1,\dots,m\}$ such that $w_i=1$ and $w_j=0$ for all $j\neq i$ , then $R_w$ is the reserve for the accident year i. Let

\begin{align*}\mathscr{D}\, :\!= \,\sigma\{C_{i,k}\,|\,(i,k)\in\mathcal{O}\} \end{align*}

be the $\sigma$ -algebra generated by the observable $C_{i,k}$ . We are interested in the conditional MSEP

\begin{align*} \operatorname{msep}_{R_w|\mathscr{D}}( \widehat{R}_w^{\,(\infty)})\;:\!=\; \operatorname{E}\!\left( (R_w -\widehat{R}_w^{\,(\infty)})^2 \,|\, \mathscr{D}\right) = \operatorname{Var}\!(R_w\,|\, \mathscr{D}) + (\widehat{R}_w^{\,(\infty)}-\operatorname{E}\!(R_w\,|\, \mathscr{D}))^2 \end{align*}

with the process variance $ \operatorname{Var}\!(R_w\,|\, \mathscr{D}) $ and the squared estimation error $(\widehat{R}_w^{\,(\infty)}-\operatorname{E}\!(R_w\,|\, \mathscr{D}))^2$ .

Estimator 1 (Process Variance). We estimate the process variance $\operatorname{Var}\!(R_w\,|\, \mathscr{D})$ by

\begin{align*} \widehat{\operatorname{Var}} (R_w\,|\, \mathscr{D}) & \, :\!= \, \sum_{i=I(n)+1}^{m} w_i^2 (\widehat{C}_{i,n}^{\,(\infty)})^2\sum_{k=K(i)+1}^{n}(\!\exp\!((\widehat{\sigma}^{\,(\infty)}_{k})^2/\widehat{C}^{\,(\infty)}_{i,k-1}+\tau_{i,k}^2)-1). \\[3pt] & \phantom{aaaaa} + 2\sum_{i=I(n)+1}^{m-1} \sum_{j=i+1}^{m} w_iw_j \cdot\widehat{C}^{\,(\infty)}_{j,n}\cdot\widehat{C}^{\,(\infty)}_{i,n}\cdot\left(\!\exp\left(\sum_{k=K(i)+1}^n \tau_{j,k-j+i} \tau_{i,k}\right)-1\right).\end{align*}

This estimator is derived in Appendix B. Note that the parameters $\tau_{i,k}$ are treated as known prior information.

In order to estimate the squared estimation error, we calculate the observed information matrix, defined as the negative Hessian of the log-likelihood function (see Oakes, Reference Oakes1999). For this purpose, let $\widehat{a}_q\, :\!= \,a_q(\widehat{\Theta}^{\,(\infty)})$ , $\widehat{b}_q\, :\!= \,b_q(\widehat{\Theta}^{\,(\infty)})$ , $\widehat{c}_{i,k}\, :\!= \,c_{i,k}(\widehat{\Theta}^{\,(\infty)})$ , and $\,\widehat{\!d}_{i,k}\, :\!= \,d_{i,k}(\widehat{\Theta}^{\,(\infty)})$ . Let

\begin{align*}\pi_2\colon \mathbb{N}\times\mathbb{N}\to \mathbb{N}, \quad (i,k)\mapsto k \end{align*}

denote the projection on the second coordinate. We define

\begin{align*}\mathcal{K}^*(q)\, :\!= \,\pi_2(\mathcal{C}^*(q))=\{k\in \{2,\dots,n\} \,|\,q-m+1 \le k\le q \}. \end{align*}

For $q\in\{2,\dots,p\}$ and $l\in\{2,\dots,n\}$ with $l\in \mathcal{K}^*(q)$ , we define

\begin{align*} \widehat{\alpha}^{\,(\mu)}_{q,l} &\, :\!= \,0, \\[3pt] \widehat{\alpha}^{\,(\sigma)}_{q,l} & \, :\!= \, \widehat{a}_q^2\cdot \frac{2\tau_{q-l+1,l}^2}{(\widehat{\sigma}_l^{\,(\infty)})^3/C_{q-l+1,l-1}},\\[3pt] \widehat{\beta}^{\,(\mu)}_{q,l} & \, :\!= \, -\widehat{a}_q \frac{\tau_{q-l+1,l}^2}{(\widehat{\sigma}_l^{\,(\infty)})^2/C_{q-l+1,l-1}} , \\[3pt] \widehat{\beta}^{\,(\sigma)}_{q,l} & \, :\!= \, \widehat{\alpha}^{\,(\sigma)}_{q,l}\sum_{(i,k)\in \mathcal{C}^*(q)} \left(\frac{ T_{i,k}-( \widehat{\mu}_k^{\,(\infty)} + \varphi_{i,k})}{(\widehat{\sigma}_k^{\,(\infty)})^2/C_{i,k-1}}+\frac{1}{2}\right)\!\cdot\tau_{i,k} - 2\widehat{a}_q\cdot \frac{T_{q-l+1,l}-(\widehat{\mu}_l^{\,(\infty)}+\varphi_{q-l+1,l})}{(\widehat{\sigma}_l^{\,(\infty)})^3/C_{q-l+1,l-1}}\cdot\tau_{q-l+1,l},\end{align*}

and for $l\not\in \mathcal{K}^*(q)$

\begin{align*}\widehat{\alpha}^{\,(\mu)}_{q,l} \, :\!= \,\widehat{\alpha}^{\,(\sigma)}_{q,l} \, :\!= \,\widehat{\beta}^{\,(\mu)}_{q,l}\, :\!= \,\widehat{\beta}^{\,(\sigma)}_{q,l}\, :\!= \,0. \end{align*}

Moreover, let

\begin{align*}\widehat{\gamma}^{\,(\mu)}_{i,k,l} \, :\!= \,\tau_{i,k}^2\widehat{\alpha}^{\,(\mu)}_{i+k-1,l},\quad \widehat{\gamma}^{\,(\sigma)}_{i,k,l} \, :\!= \,\tau_{i,k}^2\widehat{\alpha}^{\,(\sigma)}_{i+k-1,l},\quad\widehat{\delta}^{\,(\mu)}_{i,k,l} \, :\!= \,\tau_{i,k}\widehat{\beta}^{\,(\mu)}_{i+k-1,l},\quad\widehat{\delta}^{\,(\sigma)}_{i,k,l} \, :\!= \,\tau_{i,k}\widehat{\beta}^{\,(\sigma)}_{i+k-1,l}. \end{align*}

Let

\begin{align*}\delta_{k,l}\, :\!= \,\begin{cases} 1, & \mbox{if } k=l \\[3pt] 0, & \mbox{otherwise} \end{cases} \end{align*}

denote the Kronecker symbol. We define $\,\widehat{\!\mathcal{I}}^{\,\,(\infty)}\, :\!= \,(\,\widehat{\iota}^{\,\,(\infty)}_{i,k})_{i,k}\in \mathbb{R}^{2(n-1)\times 2(n-1)}$ by

\begin{align*} \widehat{\iota}^{\,\,(\infty)}_{k-1,l-1} \, :\!= \, & \delta_{k,l}\cdot\left( \sum_{i=1}^{I(k)}\frac{C_{i,k-1}}{(\widehat{\sigma}_k^{\,(\infty)})^2} \right) +\sum_{i=1}^{I(k)}\frac{\widehat{\delta}^{\,(\mu)}_{i,k,l}}{(\widehat{\sigma}_k^{\,(\infty)})^2/C_{i,k-1}} , \\[3pt] \widehat{\iota}^{\,\,(\infty)}_{k-1,n+l-2} \, :\!= \, & \widehat{\iota}^{\,\,(\infty)}_{n+l-2,k-1}\, :\!= \, 2\delta_{k,l}\cdot\left(\sum_{i=1}^{I(k)}\frac{-\,\widehat{\!d}_{i,k} + T_{i,k}-\widehat{\mu}^{\,(\infty)}_k}{(\widehat{\sigma}_k^{\,(\infty)})^3/C_{i,k-1}}\right) + \sum_{i=1}^{I(k)}\frac{\widehat{\delta}^{\,(\sigma)}_{i,k,l}}{(\widehat{\sigma}_k^{\,(\infty)})^2/C_{i,k-1}}, \\[3pt] \widehat{\iota}^{\,\,(\infty)}_{n+k-2,n+l-2} \, :\!= \, & \delta_{k,l}\cdot\sum_{i=1}^{I(k)}\left(-\frac{1}{(\widehat{\sigma}_k^{\,(\infty)})^2}+3 \frac{\widehat{c}_{i,k} +\,\widehat{\!d}^{\,2}_{i,k} -2 \,\widehat{\!d}_{i,k}(T_{i,k}-\widehat{\mu}^{\,(\infty)}_k) + (T_{i,k}-\widehat{\mu}^{\,(\infty)}_k)^2 }{(\widehat{\sigma}_k^{\,(\infty)})^4/C_{i,k-1}} +\frac{1}{4C_{i,k-1}}\right) \\[3pt] & \phantom {a} -\sum_{i=1}^{I(k)} \frac{\widehat{\gamma}^{\,(\sigma)}_{i,k,l} +2\,\widehat{\!d}_{i,k} \widehat{\delta}_{i,k,l}^{\,(\sigma)}-2\widehat{\delta}_{i,k,l}^{\,(\sigma)}(T_{i,k}-\widehat{\mu}^{\,(\infty)}_k ) }{(\widehat{\sigma}_k^{\,(\infty)})^3/C_{i,k-1}}\end{align*}

for $k,l=2,\dots,n$ .

The proof of Theorem 1 is based on a result from Oakes (Reference Oakes1999) and is found in Appendix C. We use $(\,\widehat{\!\mathcal{I}}^{\,\,(\infty)})^{-1}$ as an approximation for the covariance matrix of the estimator $\widehat{\Theta}^{\,(\infty)}$ .

In order to obtain an approximation for the squared estimation error, we define

\begin{align*}\Phi_w\colon \mathbb{R}^{n-1}\times (0,\infty)^{n-1}\to \mathbb{R} \end{align*}

by

\begin{align*}\Phi_w(\mu_2,\dots,\mu_n,\sigma_2,\dots,\sigma_n)\, :\!= \,\sum_{i=I(n)+1}^n w_i\cdot C_{i,K(i)} \prod_{k=K(i)+1}^n f_{i, k}\cdot \exp\!(\mu_{k}). \end{align*}

Then we have

\begin{align*}\Phi_w(\widehat{\Theta}^{\,(\infty)})= \sum_{i=I(n)+1}^n w_i\cdot\widehat{C}^{\,(\infty)}_{i,n} \end{align*}

and the Jacobian

\begin{align*}J_{\Phi_w}(\widehat{\Theta}^{\,(\infty)})^t=\left(\begin{array}{c} \sum_{i=I(2)+1}^n w_i\cdot \widehat{C}^{\,(\infty)}_{i,n} \\[3pt] \vdots \\[3pt] \sum_{i=I(n)+1}^n w_i\cdot \widehat{C}^{\,(\infty)}_{i,n} \\[1ex] 0 \\[3pt] \vdots \\[3pt] 0\end{array}\right). \end{align*}

Using a result of Cramér (see Witting and Müller-Funk, Reference Witting and Müller-Funk1995, Section 5.3.1), we obtain the following estimator.

Estimator 2 (Estimation error). If the observed information matrix $\,\widehat{\!\mathcal{I}}^{\,\,(\infty)}$ is invertible, then we estimate the squared estimation error $(\widehat{R}_w^{\,(\infty)}-\operatorname{E}\!(R_w\,|\, \mathscr{D}))^2$ by

\begin{align*} J_{\Phi_w}(\widehat{\Theta}^{\,(\infty)}) \cdot (\,\widehat{\!\mathcal{I}}^{\,\,(\infty)})^{-1}\cdot J_{\Phi_w}(\widehat{\Theta}^{\,(\infty)})^t. \end{align*}

Note that in Estimator 2, we neglect the bias and actually approximate the squared estimation error by the asymptotic variance of $\Phi_w(\widehat{\Theta}^{\,(\infty)})$ . The estimator makes use of the parameters $\varphi_{i,k}$ and $\tau_{i,k}$ , which are assumed to be known a priori.

5. Robustification of the method

As mentioned in Section 3, the predictions of the LMCLcy model can be instable since the variance parameters $\sigma_2,\dots,\sigma_n$ are being estimated for all development years individually. In particular, changing the prior information on the calendar year effects can change the relativities between the estimators of the $\sigma_k$ , which can lead to unexpected effects. Therefore, it is desirable to stabilize the estimation of the $\sigma_k$ and to reduce the number of parameters that are estimated in the EM algorithm.

There are various ways to tackle this issue. For instance, one could use regression or apply a variance function in each iteration of the EM algorithm (see Riegel, Reference Riegel2023, Reference Riegel2025). However, as detailed in Riegel (Reference Riegel2025), there are decisive advantages in fixing the relativities between the estimators of $\sigma_k$ before applying the EM algorithm. We therefore follow the approach from Section 6 in Riegel (Reference Riegel2025) and assumexs

\begin{align*}\sigma_k=\eta_k\sigma_0 \end{align*}

for all k, where $\eta=(\eta_2,\dots,\eta_n)$ is a known vector with $\eta_k\gt0$ and $\sigma_0\gt0$ is an unknown parameter. As we will see in Theorem 2 below, this assumption allows to prove convergence of the resulting EM algorithm under relatively weak conditions.

We suggest to choose the vector $\eta$ – before starting the EM algorithm – based on the maximum likelihood estimators $\widehat{\mu}_k^{\operatorname{LMCL}}$ and $\widehat{\sigma}_k^{\operatorname{LMCL}}$ of the LMCL model. Setting $c_{i,k}\equiv 0$ and $d_{i,k}\equiv 0$ , we can obtain these estimators from the iteration for parameter estimation in Section 3 as

\begin{align*}\widehat{\mu}_k^{\operatorname{LMCL}}\, :\!= \,\widehat{\mu}_k^{(1)}=\lim_{\kappa\to \infty}\widehat{\mu}^{(0|\kappa)} \quad \text{and}\quad \widehat{\sigma}_k^{\operatorname{LMCL}}\, :\!= \,\widehat{\sigma}_k^{(1)}=\lim_{\kappa\to \infty}\widehat{\sigma}^{(0|\kappa)}. \end{align*}

In Riegel (Reference Riegel2026), the same estimators are denoted by $\widehat{\mu}_k^{\,(\infty)}$ and $\widehat{\sigma}_k^{\,(\infty)}$ .

One possibility is to directly use $\eta_k\, :\!= \,\widehat{\sigma}_k^{\operatorname{LMCL}}$ for all $k=2,\dots,n$ . We will use this approach in our case study in the online supplementary material. This option has the advantage that one obtains the same predictions as with the LMCL method for $f_{i,k}=1$ and $t_{i,k}\to 0$ .

If the estimated expected development factors $\exp\!(\widehat{\mu}_k^{\operatorname{LMCL}})$ are all greater than 1, one can alternatively choose a monotonically increasing function $V\colon (0,\infty) \to (0,\infty)$ , for instance, a power function $V(x)=x^p$ with $p\ge 0$ , and define

\begin{align*}\eta_k^2\, :\!= \,\frac{V(\!\exp\!(\widehat{\mu}_k^{\operatorname{LMCL}})-1)}{\exp\!(\widehat{\mu}_k^{\operatorname{LMCL}})^2}. \end{align*}

The motivation for this approach is as follows: let $S_{i,k}\, :\!= \,C_{i,k}-C_{i,k-1}$ for $k\ge 2$ denote the incremental claims. Then

\begin{align*}\operatorname{Var}\left(\left.\frac{S_{i,k}}{C_{i,k-1}}\right|\mathscr{B}_{i,k-1}\right) & =\operatorname{Var}\left(\left.\frac{C_{i,k}}{C_{i,k-1}}\right|\mathscr{B}_{i,k-1}\right)= \operatorname{Var}\!(R_{i,k}|\mathscr{B}_{i,k-1})\approx \frac{\exp\!(\widehat{\mu}_k^{\operatorname{LMCL}})^2\cdot(\widehat{\sigma}_k^{\operatorname{LMCL}})^2}{C_{i,k-1}} \\[3pt]&\approx \sigma_0^2\cdot \frac{\exp\!(\widehat{\mu}_k^{\operatorname{LMCL}})^2\cdot\eta_k^2}{C_{i,k-1}} = \frac{\sigma_0^2}{C_{i,k-1}} \cdot V(\!\exp\!(\widehat{\mu}_k^{\operatorname{LMCL}})-1)\\[3pt]&\approx \frac{\sigma_0^2}{C_{i,k-1}} \cdot V\left(\!\operatorname{E}\!\left(\left.\frac{S_{i,k}}{C_{i,k-1}}\right|\mathscr{B}_{i,k-1}\right)\right),\end{align*}

that is, the meaning of V is – to some extent – similar to a variance function in a generalized linear model.

We will now describe the changes to the iteration for parameter estimation in this setup. Below, we also adjust the estimators for the MSEP to account for the given relativities $\eta_k$ . Consider the mapping

\begin{align*} \iota\colon \quad \mathbb{R}^{n-1}\times (0,\infty) & \to \mathbb{R}^{n-1}\times (0,\infty)^{n-1}, \\[3pt] (\mu,\sigma_0) & \mapsto (\mu,\eta_2\sigma_0,\dots, \eta_n\sigma_0).\end{align*}

Instead of $\Theta=(\mu_2,\dots,\mu_n,\sigma_2,\dots,\sigma_n)$ , we estimate the parameter vector $\Psi=(\mu_2,\dots,\mu_n,\sigma_0)$ using the EM algorithm. Since there is no risk of misinterpretation, we drop the $\iota$ if we apply the functions from Section 3 to $\iota(\Psi)$ . For instance, we write $a_q(\Psi)$ instead of $a_q(\iota(\Psi))$ . Given an estimate $\widetilde{\Psi}^{(\nu)}=(\widetilde{\mu}_2^{(\nu)}, \dots,\widetilde{\mu}_n^{(\nu)},\widetilde{\sigma}_0^{(\nu)})$ , we obtain the following expected log-likelihood in the E step:

\begin{align*} Q(\Psi|\widetilde{\Psi}^{(\nu)})& = \sum_{q=2}^{p} \left( - \frac{a_q(\widetilde{\Psi}^{(\nu)}) + b_q(\widetilde{\Psi}^{(\nu)})^2}{2} \right)- \sum_{(i,k)\in \mathcal{O}^*}\log(\eta_k\sigma_0) \\[3pt] & \phantom{aa} - \sum_{(i,k)\in \mathcal{O}^*}\frac{c_{i,k}(\widetilde{\Psi}^{(\nu)})+d_{i,k}(\widetilde{\Psi}^{(\nu)})^2 -2d_{i,k}(\widetilde{\Psi}^{(\nu)})(T_{i,k}-\mu_k) + (T_{i,k}-\mu_k)^2}{2\eta_k^2\sigma_0^2/C_{i,k-1}} \\[3pt] & \phantom{aa} -\frac{1}{2} \sum_{(i,k)\in \mathcal{O}^*} \left(T_{i,k}-\mu_k- d_{i,k}(\widetilde{\Psi}^{(\nu)}) \right)-\frac{1}{8} \sum_{(i,k)\in \mathcal{O}^*} \eta_k^2\sigma_0^2/C_{i,k-1} +c\end{align*}

with a constant c. Let $|M|$ denote the cardinality of a set M. In order to maximize $Q(\cdot\,|\widetilde{\Psi}^{(\nu)})$ in the modified M step, we have to solve the following equations for $\Psi$ to obtain $\widetilde{\Psi}^{\,(\nu+1)}$ :

\begin{align*} 0 & = \frac{\partial}{\partial \mu_k}Q(\Psi|\widetilde{\Psi}^{(\nu)}) =\sum_{i=1}^{I(k)}\left(-\frac{d_{i,k}(\widetilde{\Psi}^{(\nu)})+\mu_k-T_{i,k}}{\eta_k^2\sigma_0^2/C_{i,k-1}}+\frac{1}{2}\right), \\[3pt] 0 & = \frac{\partial}{\partial \sigma_0}Q(\Psi|\widetilde{\Psi}^{(\nu)}) =\sum_{(i,k)\in \mathcal{O}^*}\left(-\frac{1}{\sigma_0}+ \frac{C_{i,k-1}[c_{i,k}(\widetilde{\Psi}^{(\nu)})+( T_{i,k}-\mu_k-d_{i,k}(\widetilde{\Psi}^{(\nu)}))^2] }{\eta_k^2\sigma_0^3}-\frac{\eta_k^2\sigma_0}{4C_{i,k-1}}\right)\\[3pt] &= -\frac{|\mathcal{O}^*|}{\sigma_0} +\frac{1}{\sigma_0^3}\cdot\sum_{(i,k)\in \mathcal{O}^*} \frac{C_{i,k-1}[c_{i,k}(\widetilde{\Psi}^{(\nu)})+( T_{i,k}-\mu_k-d_{i,k}(\widetilde{\Psi}^{(\nu)}))^2] }{\eta_k^2} -\sigma_0 \cdot\sum_{(i,k)\in \mathcal{O}^*}\frac{\eta_k^2}{4C_{i,k-1}} .\end{align*}

The equations $\frac{\partial}{\partial \mu_k}Q(\Psi|\widetilde{\Psi}^{(\nu)})=0$ yield

(3) \begin{equation}\mu_k=\frac{I(k)\cdot \eta_k^2\sigma_0^2}{2\sum_{i=1}^{I(k)} C_{i,k-1}}+\frac{\sum_{i=1}^{I(k)} C_{i,k-1}(T_{i,k}-d_{i,k}(\widetilde{\Psi}^{(\nu)}))}{\sum_{i=1}^{I(k)} C_{i,k-1}}.\end{equation}

The equation $\frac{\partial}{\partial \sigma_0}Q(\Psi|\widetilde{\Psi}^{(\nu)})=0$ yields

\begin{align*}0 = \left(\sum_{(i,k)\in \mathcal{O}^*} \frac{\eta_k^2}{4C_{i,k-1}}\right)\sigma_0^4+|\mathcal{O}^*|\cdot\sigma_0^2 - \left(\sum_{(i,k)\in \mathcal{O}^*} \frac{C_{i,k-1}}{\eta_k^2}\Big[c_{i,k}(\widetilde{\Psi}^{(\nu)})+[ T_{i,k}-\mu_k- d_{i,k}(\widetilde{\Psi}^{(\nu)})]^2\Big]\right), \end{align*}

that is,

\begin{align*}\sigma_k^2=\frac{-|\mathcal{O}^*|+\sqrt{|\mathcal{O}^*|^2+\sum\limits_{(j,l)\in \mathcal{O}^*}\frac{\eta_l^2}{ C_{i,l-1}} \sum\limits_{(i,k)\in \mathcal{O}^*} \kern-1ex \frac{C_{i,k-1}}{\eta_k^2}\Big[c_{i,k}(\widetilde{\Psi}^{(\nu)})+[ T_{i,k}-\mu_k- d_{i,k}(\widetilde{\Psi}^{(\nu)})]^2\Big]}}{\frac{1}{2} \sum\limits_{(i,k)\in \mathcal{O}^*}\frac{\eta_k^2}{ C_{i,k-1}}} . \end{align*}

This leads to the following iteration.

Iteration for parameter estimation with fixed relativities. Let

\begin{align*}\widehat{m}_k\, :\!= \,\frac{\sum_{i=1}^{I(k)} C_{i,k}}{\sum_{i=1}^{I(k)} C_{i,k-1}}, \quad\widetilde{\mu}_k^{(0)}\, :\!= \, \log(\widehat{m}_k), \end{align*}

and

\begin{align*}(\widetilde{\sigma}^{(0)}_0)^2\, :\!= \,\frac{1}{|\mathcal{O}^*|}\sum_{(i,k)\in \mathcal{O}^*}\frac{C_{i,k-1}}{\eta_k^2\cdot \widehat{m}_k^2}\left(\frac{C_{i,k}}{C_{i,k-1}}-\widehat{m}_k \right)^2\!. \end{align*}

Let $\widetilde{\Psi}^{(0)}\, :\!= \,(\widetilde{\mu}_2^{(0)},\dots,\widetilde{\mu}_n^{(0)}, \widetilde{\sigma}^{(0)}_0)$ .

Given $\widetilde{\Psi}^{(\nu)}$ , the step of the EM algorithm to calculate $\widetilde{\Psi}^{\,(\nu+1)}$ is itself an iteration: let $\widetilde{\Psi}^{\,(\nu|0)}\, :\!= \,\widetilde{\Psi}^{(\nu)}$ and iterate

\begin{align*}\widetilde{\mu}_k^{(\nu|\kappa+1)}\, :\!= \,\frac{I(k)\cdot(\eta_k\widetilde{\sigma}_0^{(\nu|\kappa)})^2}{2\sum_{i=1}^{I(k)} C_{i,k-1}}+\frac{\sum_{i=1}^{I(k)} C_{i,k-1}(T_{i,k}-d_{i,k}(\widetilde{\Psi}^{(\nu)}))}{\sum_{i=1}^{I(k)} C_{i,k-1}} \end{align*}

and

\begin{align*} (\widetilde{\sigma}^{\,(\nu|\kappa+1)}_0)^2 :\!= \,\frac{-|\mathcal{O}^*|+\sqrt{|\mathcal{O}^*|^2+\sum\limits_{(j,l)\in \mathcal{O}^*}\frac{\eta_l^2}{ C_{j,l-1}}\sum\limits_{(i,k)\in \mathcal{O}^*} \frac{C_{i,k-1}}{\eta_k^2}\Big[c_{i,k}(\widetilde{\Psi}^{(\nu)})+[ T_{i,k}-\widetilde{\mu}_k^{(\nu|\kappa)}- d_{i,k}(\widetilde{\Psi}^{(\nu)})]^2\Big]}}{\frac{1}{2} \sum\limits_{(i,k)\in \mathcal{O}^*}\frac{\eta_k^2}{ C_{i,k-1}}} .\end{align*}

If this iteration converges, let $\widetilde{\Psi}^{\,(\nu+1)}\, :\!= \,\lim_{\kappa\to \infty}\widetilde{\Psi}^{\,(\nu|\kappa)}$ and continue with the next EM step.

If the algorithm converges, we define

\begin{align*}\widetilde{\mu}_k^{\,(\infty)}\, :\!= \,\lim_{\nu\to\infty}\widetilde{\mu}_k^{(\nu)},\quad \widetilde{\sigma}^{\,(\infty)}_0\, :\!= \,\lim_{\nu\to\infty}\widetilde{\sigma}^{\,(\nu)}_0, \quad \widetilde{\Psi}^{\,(\infty)}\, :\!= \,\lim_{\nu\to\infty}\widetilde{\Psi}^{\,(\nu)}, \end{align*}

and for $k=2,\dots,n$

\begin{align*}\widetilde{\sigma}_k^{\,(\infty)}\, :\!= \,\eta_k\cdot \widetilde{\sigma}^{\,(\infty)}_0. \end{align*}

Moreover, we define the posterior calendar year factors analogously to Section 3 by

\begin{align*}\widetilde{F}_{i,k}^{\,(\infty)}\, :\!= \,\operatorname{E}\!(F_{i,k}\,|\,T;\widetilde{\Psi}^{\,(\infty)})=\operatorname{E}\!(\!\exp\!(G_{i,k})\,|\,T;\widetilde{\Psi}^{\,(\infty)})=\exp\big[d_{i,k}(\widetilde{\Psi}^{\,(\infty)})+c_{i,k}(\widetilde{\Psi}^{\,(\infty)})/2\big] \end{align*}

for $(i,k)\in \mathcal{O}^*$ .

The predictions $\,\widetilde{\!C}_{i,k}^{\,(\infty)}$ can then be calculated analogously to the predictions $\widehat{C}_{i,k}^{\,(\infty)}$ in Section 3:

LMCLcy method with fixed relativities. For $(i,k)\in \mathcal{P}$ , we predict the future loss developments $C_{i,k}$ by

\begin{align*}\,\widetilde{\!C}_{i,k}^{\,(\infty)} \, :\!= \,C_{i,K(i)}\cdot \exp\!(\widetilde{\mu}^{\,(\infty)}_{K(i)+1})\cdot f_{i,K(i)+1}\cdot\ldots\cdot\exp\!(\widetilde{\mu}^{\,(\infty)}_k) \cdot f_{i,k}. \end{align*}

We define $\,\widetilde{\!R}_w^{\,\,(\infty)}$ analogously to $\widehat{R}_w^{\,(\infty)}$ above:

\begin{align*}\,\widetilde{\!R}_w^{\,(\infty)}\, :\!= \,\sum_{i=1}^m w_i \big(\,\widetilde{\!C}_{i,n}^{\,(\infty)}-C_{i,K(i)}\big). \end{align*}

We can directly use the formulas from Estimator 1 in Section 4 to estimate the process variance.

Estimator 3 (Process Variance). In case of fixed relativities, we estimate the process variance $\operatorname{Var}\!(R_w\,|\, \mathscr{D})$ by

\begin{align*} \widehat{\operatorname{Var}} (R_w\,|\, \mathscr{D}) & \, :\!= \, \sum_{i=I(n)+1}^{m} w_i^2 (\,\widetilde{\!C}_{i,n}^{\,(\infty)})^2\sum_{k=K(i)+1}^{n}(\!\exp\!((\widetilde{\sigma}^{\,(\infty)}_{k})^2/\,\widetilde{\!C}^{\,(\infty)}_{i,k-1}+\tau_{i,k}^2)-1). \\[3pt] & \phantom{aaaaa} + 2\sum_{i=I(n)+1}^{m-1} \sum_{j=i+1}^{m} w_iw_j \cdot\,\widetilde{\!C}^{\,(\infty)}_{j,n}\cdot\,\widetilde{\!C}^{\,(\infty)}_{i,n}\cdot\left(\!\exp\left(\sum_{k=K(i)+1}^n \tau_{j,k-j+i} \tau_{i,k}\right)-1\right)\!.\end{align*}

For the squared estimation error, the formulas from Section 4 have to be adjusted. Let $J_\iota(\Psi)$ denote the Jacobian matrix of $\iota$ at $\Psi$ , that is,

We calculate $\widetilde{\mathcal{I}}^{\,(\infty)}$ using the definition of $\,\widehat{\!\mathcal{I}}^{\,\,(\infty)}$ from Section 4 but replacing $\widehat{\Theta}^{\,(\infty)}$ by $\iota(\widetilde{\Psi}^{\,(\infty)})$ . Then $\widetilde{\mathcal{I}}^{\,(\infty)}$ is the observed information matrix at $\iota(\widetilde{\Psi}^{\,\,(\infty)})$ (see Theorem 1). The observed information matrix $\,\,\,\,\widetilde{\!\!\mathcal{J}}^{\,\,(\infty)}$ at $\widetilde{\Psi}^{\,(\infty)}$ can be calculated as

\begin{align*}\,\,\,\,\widetilde{\!\!\mathcal{J}}^{\,\,(\infty)}=J_\iota(\widetilde{\Psi}^{\,(\infty)})^t\widetilde{\mathcal{I}}^{\,(\infty)}J_\iota(\widetilde{\Psi}^{\,(\infty)}). \end{align*}

Let $\Phi^\iota_w\, :\!= \,\Phi_w\circ \iota$ . Then, $\Phi^\iota_w\colon \mathbb{R}^{n-1}\times (0,\infty)\to \mathbb{R}$ is given by

\begin{align*}\Phi^\iota_w(\Psi)\, :\!= \,\sum_{i=I(n)+1}^n w_i\cdot C_{i,K(i)} \prod_{k=K(i)+1}^n f_{i, k}\cdot \exp\!(\mu_{k}). \end{align*}

We have

\begin{align*}\Phi^\iota_w(\widetilde{\Psi}^{\,(\infty)})= \sum_{i=I(n)+1}^n w_i\cdot\,\widetilde{\!C}^{\,(\infty)}_{i,n} \end{align*}

and

\begin{align*}J_{\Phi^\iota_w}(\widetilde{\Psi}^{\,(\infty)})^t=\left(\begin{array}{c} \sum_{i=I(2)+1}^n w_i\cdot \,\widetilde{\!C}^{\,(\infty)}_{i,n} \\[3pt] \vdots \\[3pt] \sum_{i=I(n)+1}^n w_i\cdot \,\widetilde{\!C}^{\,(\infty)}_{i,n} \\[1ex] 0\end{array}\right). \end{align*}

Estimator 4 (Estimation error). If the information matrix $\,\,\,\,\widetilde{\!\!\mathcal{J}}^{\,\,(\infty)}$ is invertible, then we estimate the squared estimation error $(\,\widetilde{\!R}_w^{\,(\infty)}-\operatorname{E}\!(R_w\,|\, \mathscr{D}))^2$ by

\begin{align*} J_{\Phi^\iota_w}(\widetilde{\Psi}^{\,(\infty)}) \cdot (\,\,\widetilde{\!\!\mathcal{J}}^{\,\,(\infty)})^{-1}\cdot J_{\Phi^\iota_w}(\widetilde{\Psi}^{\,(\infty)})^t. \end{align*}

To study the existence of MLEs and the convergence of the EM algorithm, the following definition is needed.

Definition. We say that the observed triangle is singular if there exist $\mu_2,\dots,\mu_n$ and $x_2,\dots,x_p$ such that

\begin{align*} T_{i,k}=\mu_k+\varphi_{i,k}+\tau_{i,k}x_{i+k-1} \end{align*}

for all $(i,k)\in \mathcal{O}^*$ .

Theorem 2. Let the observed triangle be non-singular. Then the log-likelihood $\widetilde{\ell}_T\, :\!= \,\ell_T\circ \iota$ of the model with fixed relativities has a global maximum in $\mathbb{R}^{n-1}\times (0,\infty)$ . If there is a unique MLE

\begin{align*} \widetilde{\Psi}^{\operatorname{MLE}}=(\tilde{\mu}_2^{\operatorname{MLE}},\dots,\tilde{\mu}_n^{\operatorname{MLE}},\widetilde{\sigma}_0^{\operatorname{MLE}}) \end{align*}

and

\begin{align*} \frac{I(k)\cdot \eta_k^2(\widetilde{\sigma}_0^{\operatorname{MLE}})^2}{2\sum_{i=1}^{I(k)}C_{i,k-1}}\lt1 \end{align*}

for $k=2,\dots,n$ , then there exists a neighborhood U of $\widetilde{\Psi}^{\operatorname{MLE}}$ , such that the EM iteration with fixed relativities converges to the MLE for all initial parameter vectors in U.

A proof of Theorem 2 is provided in Appendix D.

6. Conclusions

This paper introduces the LMCLcy model, an extension of the LMCL model, that enables the incorporation of external market information on calendar year effects within a principled framework. Calendar year effects are represented by multiplicative log-normal factors with known prior moments, and model parameters are estimated using an EM algorithm. Closed-form expressions for prediction uncertainty are derived, and a robustified variant of the method is proposed to enhance stability. For this robustified version, convergence of the EM algorithm is established under mild regularity conditions on the loss development triangle.

Simulation results presented in the online supplementary material indicate that, for triangles with low volatility, the LMCLcy method – with reasonable choices of $\rho$ – provides an effective framework for incorporating prior market information on calendar year effects. For triangles with medium to high volatility, the LMCLcy method is robust and does not exhibit convergence issues in the EM algorithm. However, in these settings, it does not improve predictive performance relative to the classical chain ladder method. For extremely volatile triangles, the LMCLcy method may become unstable, in which case the classical chain ladder method is the preferrable. Overall, the proposed framework provides a practical bridge between market-level information and portfolio-specific reserving.

It should be noted that the LMCLcy method is specifically designed for long-tail triangles; its applicability to short-tail segments appears limited.