# A Bayesian Skew-heavy-tailed modelling for loss reserving

William L. Leão1 and Viviana G. R. Lobo1111Corresponding author: Viviana G R
Lobo, Departamento de Métodos Estatísticos, Instituto de Matemática,
Universidade Federal do Rio de Janeiro, Av. Athos da Silveira Ramos, Centro de
Tecnologia, Bloco C, CEP 21941-909.
E-mail<EMAIL_ADDRESS>Homepage:
https://sites.google.com/site/dme/viviana
1Instituto de Matemática, Departamento de Métodos Estatísticos
Universidade Federal do Rio de Janeiro, Brazil

###### Abstract

This paper focuses on modelling loss reserving to pay outstanding claims. As
the amount liable on any given claim is not known until settlement, we propose
a flexible model via heavy-tailed and skewed distributions to deal with
outstanding liabilities. The inference relies on Markov chain Monte Carlo via
Gibbs sampler with adaptive Metropolis algorithm steps allowing for fast
computations and providing efficient algorithms. An illustrative example
emulates a typical dataset based on a runoff triangle and investigates the
properties of the proposed models. Also, a case study is considered and shows
that the proposed model outperforms the usual loss reserving models well
established in the literature in the presence of skewness and heavy tails.

Keywords: bayesian inference; heavy-tailed; skewness; loss reserving; adaptive
MCMC

##  1 Introduction

###  1.1 Background

An insurance policy guarantees the payment provided by an insurer to the
policyholder if a loss occurs conditionally on a premium payment. However, in
some cases, insurance companies often do not settle outstanding claims
immediately because some claims are reported with some delay, which means that
claims should be paid with a time delay of months, years, or decades. In this
context, is crucial for the insurance company be the ability to predict as
accurately as possible the amount necessary to compose the reserve for claims
that incurred but not yet reported (IBNR, which means the liabilities that
have arisen but are reported after the accounting date) or settled to maintain
their solvency.

Traditionally, some of the existing statistical approaches for dealing with
the problem, such as the chain ladder method (Taylor, (1984)),
Bornhuetter–Ferguson algorithms (Bornhuetter and Ferguson, (1972)) and Mack
method (Mack, (1993)) to compute loss reserves are not able to provide
satisfactory performance, because is based on deterministic methods that
return a pointwise estimate and do not take into account stochastic
uncertainties involved. Also, the usual techniques are not flexible to capture
irregular and extreme claims in the dataset.

Several papers have proposed models based on stochastic techniques or
modifications of the usual chain-ladder method for loss reserving. For
example, Verrall, (2000) investigates the relationship between chain-ladder
and some stochastic models. Other authors, Kremer, (1982), Renshaw and
Verrall, (1998), Verrall, (1996), Haberman and Renshaw, (1996) consider the
use of log-linear models, such as the analysis of variance via ANOVA-type and
ANCOVA-type modelling. Verrall, (1989), de Jong and Zehnwirth, (1983),
Verrall, (1990), and Ntzoufras and Dellaportas, (2002) consider flexibility of
the models based on state-space. For more details about loss reserving
stochastic methods, see Ntzoufras and Dellaportas, (2002) and England and
Verrall, (2002).

A usual assumption to models loss reserving takes into account that reserves
have Gaussian errors (see Verrall, (1989), Ntzoufras and Dellaportas, (2002),
and among others shown previously). However, these processes often have common
characteristics such as non-Gaussianity, outlying claims, and non-constant
variance. In this way, if these behaviour are not incorporated into the model
assumption, the predictions for loss reserving might be poor. Relaxing the
assumption, some proposals consider modelling reserves via scale mixture
forms. Chan et al., (2008) proposed a generalized-t error distribution via a
scale mixture of uniforms (SMN) to take into account the presence of potential
outliers. Choy et al., (2016) consider scale mixtures (Gaussian and Uniform)
based on heavy-tailed distributions including Cauchy, Student-t, EP, logistic,
and Laplace in a univariate loss reserving data modelling. They consider
modelling the amount of claims paid to the insureds of an insured product
during the 18 delayed years and conclude that the log-ANCOVA-t model affords
better performance in the heavy-tailed class proposed. Shi et al., (2012) and
Goudarzi and Zokaei, (2021) consider the same approach by Verrall, (1989) and
Choy et al., (2016) however in a multivariate loss reserve setting. The
authors include an additional parameter associated with the calendar year. The
calendar effects are captured by a random effect term, via a random-walk model
or an autoregressive model. Recently, Usman and Chan, (2022) propose adopting
the conditional autoregressive range model to the mean function of the loss
reserve distribution to explain the persistence of the development year
effect.

Although these models allow fatter tails than Normal, they are not able to
capture skewed claims. In this manner, Azzalini, (1985) and Azzalini and
Valle, (1996) propose a class of distributions based on skew-normal (SNLM).
Other authors, such as Arellano-Valle et al., (2005) and Arellano-Valle et
al., (2007) propose skew-normal linear mixed models (SNLMM). Pigeon et al.,
(2013) consider a family of multivariate skew-symmetric (MSS) based on the
skew-normal distribution to model individual reserve claims. Although this
model allows skewness, they do not take into account heavy-tail distributions.
To allow simultaneously both skewness and thick tails in the modelling,
Fernández and Steel, (1998) propose a family of univariate skew-t
distribution. Another family of skew-t distribution is proposed by Kim, (2001)
extending the idea in Azzalini, (1985) where the skew-normal is a special
case. Furthermore, Wang and Genton, (2006) introduce the multivariate skew-
slash model as an alternative to skew-t and skew-normal models.

Choy et al., (2016) address the importance of improving their proposed model
considering modelling skewness. However, few propose to accommodate both
skewness and high kurtosis have been proposed yet in the loss reserving
literature. Vilela, (2013) performs the estimation of the log-ANOVA and
dynamic ANOVA models using the normal, Student’s t, Skewed Normal and Skewed
Student’s t distributions. Vieira, (2019) considers a class of loss-reserving
models based on the quantile regression approach developed for modelling the
tails of the claims considering the asymmetric Laplace distributions and Usman
and Chan, (2022) consider the Generalized beta type II distribution to capture
the characteristics of loss reserving distribution including skewness,
peakness and long right tails. We propose flexible models based on a family of
Skew-heavy-tailed distributions that deal with outstanding liabilities. Our
class of models include Student-t, Slash and Variance-Gamma (also known as
generalized Laplace, Madan and Seneta, (1990)) based on a scale mixture of
skew-Normal distribution. This proposal gives a suitable characterization for
reserve claim data. The inference procedure follows the Bayesian framework
that has the advantage of taking into account the uncertainty associated with
the estimates through the posterior distribution. Also, we can access the loss
reserving of the predictive distribution for outstanding liabilities and the
prior knowledge might be elicited over past experience from an insurance
company.

From a full Bayesian view, the posterior distribution is commonly not
available from their analytical form, and numerical integration and simulation
are considered, in particular, Markov chain Monte Carlo (MCMC) methods
(Gamerman and Lopes,, 2006) are used in this work. Several actuarial papers
for loss reserving such as, Choy et al., (2016), Goudarzi and Zokaei, (2021),
and Usman and Chan, (2022) obtain the posterior distribution via MCMC through
open sources software (OpenBUGS, WinBUGS and RStan Lunn et al., (2000), Lunn
et al., (2009), Stan Development Team, (2023)). Although BUGS and Stan are
user-friendly Bayesian software, it could be very restrictive if considering
additional complexity in the modelling process and/or computationally
exhaustive. For this, we address the posterior distribution via adaptive MCMC
(Andrieu and Thoms, (2008)) jointly with the Gibbs sampler techniques,
allowing for fast and smart feasible computations. For the Skew-Variance-Gamma
distribution, we regard the technique proposed by Philippe, (1997) to generate
truncated Gamma distributions.

Briefly, our proposed method for loss reserving allows accommodating skewness,
heavy tails, and outlying claims behaviour in the dataset considering Skew-
heavy-tailed models and taking into account Bayesian inference procedure
through Gibbs sampler with adaptive Metropolis algorithm steps MCMC techniques
with computational efficiency that we can obtain through the development of
the proposed method in the free software R, with the use of the Rcpp library
(Eddelbuettel,, 2013; Eddelbuettel and Balamuta,, 2018) that provides easy
integration between the programming language C++ and R, and through the use of
the linear algebra library Armadillo for language C++ (Sanderson and Curtin,,
2016, 2018) and loops implemented efficiently, we were able to save
considerable computational time.

###  1.2 Outline of the paper

The remaining of the paper is organised as follows. Section 2 describes the
class of Skew-heavy-tailed distributions. In subsection 2.2, the Bayesian
procedure for loss reserving is introduced and algorithms are presented (see
subsection 2.3). Section 3 includes some simulation studies and examples that
have already been presented in Choy et al., (2016) in order to illustrate the
performance of our proposed models. Concluding remarks are given in Section 4.
Some aspects related to properties of the proposed model, the posterior
distribution, model comparison, and real datasets, are presented in Appendices
A, B, C and D, respectively.

##  2  Skew-heavy-tailed modelling

In this section, we propose new approaches modelling based on a family of
Skew-heavy-tailed models for outstanding claims predictions. We define a scale
mixture of the skew-normal family and study some of its properties.

The insurance data follows a structure typically presented by a run-off
triangle (see Table 1), showing the development of the claims over time for
each period. Let a random variable $Y_{ij}$ the amount that the insurance
company paid with a delay of $j-1$ years (development year) for accidents that
occurred in year $i$ (policy year), $i=1,\ldots,n$, $j=1,\ldots,n-i+1$. After
the business has been running for $t$ years, the data available has the form
$Y_{\Delta}:=\left\\{Y_{ij}:j\leq n-i+1,i\leq n\right\\}$ and $Y_{ij}>0$,
$\forall i$ and $\forall j$ (Verrall, (1989), Kremer, (1982)), usually called
run-off triangle. For example, $y_{2,2}$ means that the event occurred in year
2 and was reported in the same year, whereas $y_{2,6}$ means that the event
was settled after four years before the policy year. The main interest lies on
predicting the unobserved claim amounts, which means the lower triangle (light
grey cells, $Y_{\nabla}:=\left\\{Y_{ij}:j>n-i+1\right\\}$) should be
estimated.

Table 1: Run-off triangle: structure of claim amount dataset. policy year |  lag year $j$
---|---
$i$ | 1 | 2 | $\ldots$ | $n-1$ | $n$
1 | $y_{1,1}$ | $y_{1,2}$ | $\ldots$ | $y_{1,n-1}$ | $y_{1,n}$
2 | $y_{2,1}$ | $y_{2,2}$ | $\ldots$ | $y_{2,n-1}$ |
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ |
$n-1$ | $y_{n-1,1}$ | $y_{n-1,2}$ |
$n$ | $y_{n,1}$ |

###  2.1 A class of Skew-heavy-tailed distributions

Let $Z_{ij}=log(Y_{ij})$ be the log claim amount paid to the policyholder from
the policy-year $i$ with $n-1$ years of delay, and $j\leq n-i+1$, $i\leq n$.
The stochastic representation of the variable for loss reserving is given by
the formulation as a mixture of scales as follows

$Z_{ij}=\mu_{ij}+\sigma\frac{\varepsilon_{ij}}{\sqrt{\lambda_{ij}}},$ (1)

where $\mu_{ij}$ represents the mean of the log claim amount, originated in
year $i$ and was paid with a delay of $n-1$ years. The variable
$\varepsilon_{ij}$ is a zero mean Skew-Normal distribution, with $\sigma^{2}$
as a scale parameter, and $\kappa$ skewness parameter. The component
$\lambda_{ij}$ is a positive mixing variable and it is responsible for
controlling the tail of the distribution and is independent of
$\varepsilon_{ij}$. The p.d.f. of $Z_{ij}$ is said to have a Skew-Mixture-
Normal representation with $\mu_{ij}\in\Re$, $\sigma^{2}\in\Re^{+}$ and
$\kappa\in\Re$ if it can be expressed as

$\displaystyle f_{Z}(z\mid\mu,\sigma,\kappa)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathcal{S}\mathcal{N}(z;\mu,\sigma^{2},\lambda,\kappa)p(\lambda)d\lambda,$
(2) $\displaystyle=$ $\displaystyle
2\int_{0}^{\infty}\phi(z\mid\mu,\lambda^{-1}\sigma^{2})\Phi\left(\kappa\frac{z-\mu}{\sigma}\right)p(\lambda)d\lambda,$

where $\phi$ follows a Gaussian distribution and $\Phi$ follows a cumulative
distribution of the standard Normal distribution, respectively. Conditional on
$\lambda_{ij}$, the variable $Z_{ij}$ has to mean function $\mu_{ij}$ that is
specified for parameters that depend on policy-year and lag-year effects,
variance $\lambda_{ij}^{-1}\sigma^{2}$ and skewness
$\lambda_{ij}^{-1/2}\kappa$. Gaussian behaviour is only assumed given the
mixing variable $\lambda$. The representation is given by
$Z_{ij}\mid\mu_{ij},\sigma^{2},\lambda_{ij},\kappa\sim\mathcal{S}\mathcal{N}(\mu_{ij},\lambda^{-1}_{ij}\sigma^{2},\lambda_{ij}^{-1/2}\kappa)$.
Notice that when $\lambda_{ij}\neq 1$, $Z_{ij}$ has a heterogeneous variance,
and if $\lambda_{ij}$ is integrated out the resulting is non-Gaussian
distribution. On the other hand, if $\lambda_{ij}=1$, the resulting model is
the usual Skew-Normal (Azzalini and Valle, (1996)). For $\kappa=0$ and
conditional on $\lambda_{ij}$, the Skew-Normal distribution reduces to the
Normal distribution, that is, the loss reserving follows
$Z_{ij}\mid\mu_{ij},\sigma^{2},\lambda_{ij}\sim\mathcal{N}(\mu_{ij},\lambda_{ij}^{-1}\sigma^{2})$.

To keep the model feasible for the inference procedure, we define the loss
reserving modelling through a hierarchical representation of distributions
according to the stochastic representation proposed by Henze, (1986). Let
$T_{1}$ and $T_{2}$ standard Normal variables. If $\varepsilon$ follows a
standard Skew-Normal distribution with density given by equation (2) with
$\lambda=1$, then $\varepsilon$ is able to be written as a location-scale
mixture of Gaussian distribution given by

$\varepsilon=\rho|T_{1}|+\sqrt{(1-\rho^{2})}T_{2},$ (3)

where $\rho=\frac{\kappa}{\sqrt{1+\kappa^{2}}}$, $|T_{1}|$ denotes the
absolute value of $T_{1}$, and $T_{1},T_{2}$ are independent
$\mathcal{N}(0,1)$. As we can see, $\varepsilon$ is written as a mixture of
variables from $T$ that follows a standard Gaussian distribution.

Thus, equation (1) can be written as

$Z_{ij}=\mu_{ij}+\sigma\lambda_{ij}^{-1/2}\left(\rho|T_{1;i,j}|+\sqrt{(1-\rho^{2})}T_{2;i,j}\right).$
(4)

Therefore,

$\displaystyle Z_{ij}\mid T_{ij},\lambda_{ij}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}\left(\mu_{ij}+\rho
T_{ij};\sigma^{2}\lambda_{ij}^{-1}(1-\rho^{2})\right)$ $\displaystyle
T_{ij}\mid\lambda_{ij}$ $\displaystyle\sim$
$\displaystyle\mathcal{N}\mathcal{T}\left(0,\sigma^{2}\lambda_{ij}^{-1};(0,\infty)\right)$
(5) $\displaystyle\lambda_{ij}$ $\displaystyle\sim$ $\displaystyle
p(\cdot\mid\nu),$

where $T_{ij}=\lambda_{ij}^{-1/2}|T_{1;i,j}|$ and $\nu$ controls the tails of
the distribution. Then, to generate a Gaussian distribution, we proceed
through three steps, that is, we generate first from the mixture distribution
$\lambda_{ij}$, after from the conditional $T_{ij}\mid\lambda_{ij}$ generate a
Truncate Gaussian distribution and next from the conditional distribution
$Z_{ij}\mid T_{ij},\lambda_{ij}$ using the stochastic representation given in
equation (4). The p.d.f. of $Z_{ij}$ is rewritten as

$\displaystyle f_{Z}(z\mid\cdot)$
$\displaystyle=\int_{-\infty}^{\infty}\int_{0}^{\infty}f(z,\lambda,t\mid\cdot)dtd\lambda$
$\displaystyle=\int_{-\infty}^{\infty}\int_{0}^{\infty}f(z|\lambda,t,\cdot)f(t\mid\lambda,\cdot)p(\lambda)dtd\lambda$
$\displaystyle=\int_{-\infty}^{\infty}\int_{0}^{\infty}\mathcal{N}\left(z\mid\mu+t\rho,\sigma^{2}\lambda^{-1}(1-\rho^{2})\right)\mathcal{NT}\left(z\mid
0,\sigma^{2}\lambda^{-1};(0,\infty)\right)p(\lambda\mid\nu)dtd\lambda.$ (6)

Some properties of the Skew-Heavy-Tailed class The resulting mean and variance
of the distribution $Z_{ij}$, defined in equation (1) is given by

$E(Z_{ij}\mid\cdot)=\mu_{ij}+\sigma\rho\sqrt{\frac{2}{\pi}}E\left(\lambda_{ij}^{-1/2}\right)$
(7)

$Var(Z_{ij}\mid\cdot)=\sigma^{2}\left(E\left(\lambda_{ij}^{-1}\right)-\rho^{2}\left(\frac{2}{\pi}\right)E^{2}\left(\lambda_{ij}^{-1/2}\right)\right).$
(8)

Notice that if $\lambda_{ij}=1$, then
$E(Z_{ij})=\mu_{ij}+\sqrt{\frac{2}{\pi}}\sigma\rho$ and
$Var(Z_{ij})=\sigma^{2}\left(1-\frac{2}{\pi}\rho^{2}\right)$ follows an usual
Skew-Normal distribution, it implies that the distribution does not have
ticker tails than Skew-Normal. Also, if $\lambda_{ij}=1$ and $\kappa=0$, then
$\rho=0$ and we return a Gaussian distribution. As follows, we addressed the
prior specification of the mixing variable $\lambda_{ij}$ and the parameter
associated with skewness $\kappa$.

Modelling the mixing $\lambda_{ij}$ A first and natural approach would be to
assume that
$\lambda_{ij}\sim\mathcal{G}\left(\frac{\nu}{2},\frac{\nu}{2}\right)$. This
implies that $Z_{ij}$ is a Skew-Student-t (or just Skew-t) distribution with
$\nu>0$ degrees of freedom, denoted by
$Z_{ij}\mid\mu_{ij},\sigma^{2},\nu,\kappa\sim\mathcal{S}\mathcal{S}\mathcal{T}(\mu_{ij},\sigma^{2},\nu,\kappa)$
and the Gamma distribution ($\mathcal{G}(\cdot,\cdot)$). Therefore, the
expected value and variance given in equations (7) and (8), unconditional on
the mixing parameter are obtained by computing
$E(\lambda_{ij}^{-1/2})=\left[{\Gamma(\frac{\nu-1}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{1/2}$
and $E\left(\lambda_{ij}^{-1}\right)=\frac{\nu}{2}/(\frac{\nu}{2}-1)$, with
$\frac{1}{\lambda_{ij}}\sim\mathcal{I}\mathcal{G}(\frac{\nu}{2},\frac{\nu}{2})$.

Alternatively, if we assume a Skew-Slash model, then $\lambda_{ij}$ has Beta
distribution ($\mathcal{B}e(\cdot,\cdot)$), with mean is $\frac{\nu}{\nu+1}$,
and variance is $\frac{\nu}{(\nu+1)^{2}(\nu+2)}$, hence
$\lambda_{ij}\sim\mathcal{B}{e}(\nu,1)$, and
$Z_{ij}\mid\mu_{ij},\sigma^{2},\nu,\kappa\sim\mathcal{S}\mathcal{S}(\mu_{ij},\sigma^{2},\nu,\kappa)$.
This implies that the prior marginal distribution of $\lambda_{ij}$ is
concentrated around one for bigger values of $\nu$, and as $\nu$ decreases,
the distribution becomes more spread out. Analogously,
$E(\lambda_{ij}^{-1/2})=\frac{\nu}{\nu-1/2}$ and
$E\left(\lambda_{ij}^{-1}\right)=\frac{\nu}{\nu-1}$, respectively.

Another distribution that belongs to the skew-heavy-tailed family assumes that
$\lambda_{ij}\sim\mathcal{I}\mathcal{G}(\frac{\nu}{2},\frac{\nu}{2})$, where
$\nu$ is responsible to capture ticker tails, then
$Z_{ij}\mid\mu_{ij},\sigma^{2},\nu,\kappa\sim\mathcal{S}\mathcal{V}\mathcal{G}(\mu_{ij},\sigma^{2},\nu,\kappa)$.
The mean and variance of Skew-Variance-Gamma distribution follow the equations
that have been mentioned previously considering
$E(\lambda_{ij}^{-1/2})=\left[{\Gamma(\frac{\nu+1}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{-1/2}$
and $E\left(\lambda_{ij}^{-1}\right)=1$. For more details about the
distributions and their properties see Appendix A.

|
---|---
(a) symmetric-heavy-tail distributions | (b) skew-heavy-tail distributions
Figure 1: Illustration for the class of heavy-tailed distribution compared
with Normal distribution (solid line): (a) symmetric distributions and (b)
asymmetric distributions, Student-t (dashed line), Slash (dotted line), and
Variance-Gamma (dotdash line).

Figure 1 illustrates the behaviour of heavy tails distribution when compared
to the Gaussian ($\mu=0$, $\sigma=1$). As we can see, Panel (a) presents the
symmetric distributions: the Student-t ($\mu=0$, $\sigma=1$, $\nu=3$), the
Slash ($\mu=0$, $\sigma=1$, $\nu=3$) and the Variance-Gamma ($\mu=0$,
$\sigma=1$, $\nu=1$). In particular, for the Student-t distribution and the
Slash distribution, as $\nu\rightarrow\infty$ we retrieve Normal tails, while
smaller values of $\nu$ will induce thicker tails. Moreover, for the Variance-
Gamma, parameters $\mu$, $\sigma$, and $\nu=3$ are responsible for, location,
spread, and excess kurtosis (peakedness), respectively. Notice that, the tails
of the Variance-Gamma distribution decrease slower than the Normal
distribution and for $\frac{1}{\nu}\rightarrow 0$ implies Normal tails because
kurtosis is $3(1+\frac{1}{\nu})\approx 3$. Otherwise, larger values of $\nu$
suggest peakedness. Panel (b) of Figure 1 shows an example of skewness in the
distributions with $\kappa=5$.

Modelling the mean $\mu_{ij}$ The mean function can be modelled assuming a
dynamic linear structure where the parameters depend on each other in a time-
recursive way. The idea is to take into account the interaction between
accident and development year in the mean. The dynamic linear mean can be
written as

$\mu_{ij}=\mu+\alpha_{i}+\beta_{ij}+\gamma_{t=i+j-1},$ (9)

for $i=2,\ldots,n$, $j=2,\ldots,n-i+1$ and $\mu$ is a general mean between $i$
and $j$, $\alpha_{i}$, $\beta_{ij}$ and $\gamma_{t}$ represent a factor which
reflects expected changes because of the accident year, a factor the
interaction among accident year and development year, and calendar year
effects, respectively. Notice that the dynamic variable $\beta_{ij}$ allows
the development year pattern to change with accident years, which means, they
capture the accident-development year interaction effects with hardly any
parameters. Following a Bayesian formulation, $\alpha_{i}=\alpha_{i-1}+u_{i}$,
$u_{i}\sim N(0,\sigma^{2}_{\alpha})$, $\beta_{ij}=\beta_{i-1,j}+v_{i}$,
$v_{i}\sim N(0,\sigma^{2}_{\beta})$ and $\gamma_{t}$ as a random walk,
$\gamma_{t}=\gamma_{t-1}+\omega_{t}$, $\omega_{t}\sim
N(0,\sigma^{2}_{\gamma})$. Also, for parameter identification set the
constraints as $\alpha_{1}=\beta_{1j}=\gamma_{1}=0$. More details see
Ntzoufras and Dellaportas, (2002), Goudarzi and Zokaei, (2021), Shi et al.,
(2012) and Usman and Chan, (2022).

###  2.2 Posterior distribution and inference procedure

Assume that for each period where $i\leq n$ and $j\leq n-i+1$ the loss
reserving observations are available. Let
$\mathbf{z}=(z_{1}^{\prime},\ldots,z_{n}^{\prime})$ the observational vector,
$\mathbf{T}=(T_{1}^{\prime},\ldots,T_{n}^{\prime})$, and mixing vector
$\boldsymbol{\lambda}=(\lambda_{1}^{\prime},\ldots,\lambda_{n}^{\prime})$, the
vectors that contains the information about $z_{ij}$, $T_{ij}$, and
$\lambda_{ij}$, respectively. The distribution obtained from equation (1) is
given by hierarchical structure as follows

$\displaystyle Z_{ij}\mid T_{ij},\lambda_{ij}$
$\displaystyle\sim\mathcal{N}\left(\mu+\alpha_{i}+\beta_{ij}+\gamma_{t}+\rho
T_{ij},\lambda_{ij}^{-1}\sigma^{2}\left(1-\rho^{2}\right)\right)$
$\displaystyle\alpha_{i}$
$\displaystyle\sim\mathcal{N}\left(\alpha_{i-1},\sigma_{\alpha}^{2}\right)$
$\displaystyle\beta_{ij}$
$\displaystyle\sim\mathcal{N}\left(\beta_{i-1,j},\sigma_{\beta}^{2}\right)$
$\displaystyle\gamma_{t}$
$\displaystyle\sim\mathcal{N}\left(\gamma_{t-1},\sigma_{\gamma}^{2}\right)$
$\displaystyle T_{ij}\mid\lambda_{ij}$
$\displaystyle\sim\mathcal{NT}_{[0,\infty)}\left(0,\sigma^{2}\lambda_{ij}^{-1}\right)$
$\displaystyle\lambda_{ij}$ $\displaystyle\sim p(\lambda_{ij}\mid\nu),$ (10)

where $z_{i}^{\prime}=(z_{i1},\ldots,z_{in})$,
$T_{i}^{\prime}=(T_{i1},\ldots,T_{in})$ and
$\lambda_{i}^{\prime}=(\lambda_{i1},\ldots,\lambda_{in})$. As mentioned
previously, we follow the Bayesian approach to make inference, predictions,
access the goodness-of-fit, and model comparison that are obtained from the
joint posterior distribution of the parameters. In particular, we address the
hierarchical structure of our proposed model in equation (2.2) in the proposed
estimation algorithm to generate from the joint posterior and to make
predictions.

The model description is complete after assigning a prior distribution for the
static parameter vector
$\boldsymbol{\theta}=\left(\mu,\rho,\sigma^{2},\sigma_{\alpha}^{2},\sigma_{\beta}^{2},\sigma_{\gamma}^{2},\nu\right)$.
We assign vague independent priors to the parameters in $\boldsymbol{\theta}$.
A sensitivity analysis of priors is considered for these choices (see Section
3.2). In particular, we assume $\mu\sim\mathcal{N}(0,s_{\mu}^{2}=100)$,
$(1+\rho)/2\sim\mathcal{B}e(c=1,d=1)$,
$\sigma^{2}\sim\mathcal{G}(a_{\sigma},b_{\sigma})$,
$\sigma_{\alpha}^{2}\sim\mathcal{G}(a_{\alpha},b_{\alpha})$,
$\sigma_{\beta}^{2}\sim\mathcal{G}(a_{\beta},b_{\beta})$ e
$\sigma_{\gamma}^{2}\sim\mathcal{G}(a_{\gamma},b_{\gamma})$, with $a=0.001$
and $b=0.001$. For mixing parameter $\nu$ associated with the Student-t, the
Slash, and the Variance-Gamma distribution, we assign a prior
$\mathcal{G}(a_{t}=12,b_{t}=0.8)$, $\mathcal{G}(a_{s}=0.2,b_{s}=0.05)$,
$\mathcal{G}(a_{vg}=12,b_{vg}=0.8)$, respectively. Notice that larger values
of $\nu$ lead to approximate skew-Normal, while smaller values of $\nu$
suggest very thick tails, with $\kappa>0$ and hence $-1\leq\rho\leq 1$.

Let
$\boldsymbol{\Theta}=(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\theta})$
be the parameter vector associated with the skew-heavy-tailed model such that,
$\boldsymbol{\alpha}=(\alpha_{2},\ldots,\alpha_{n})^{\prime}$,
$\boldsymbol{\beta}=(\beta_{2}^{\prime},\ldots,\beta_{n}^{\prime})$, with
$\beta_{i}^{\prime}=(\beta_{i1},\ldots,\beta_{in})$, and
$\boldsymbol{\gamma}=(\gamma_{2},\ldots,\gamma_{n})$ the dynamical parameters
for $i=2,\ldots,n$ and $j=2,\ldots,n$ with
$\alpha_{1}=\beta^{\prime}_{1}=\gamma_{1}=0$. Following Bayes’ theorem, the
posterior distribution of $\boldsymbol{\Theta}$, $\boldsymbol{\lambda}$, ${\bf
T}$, given the observed data, $\mathbf{z}$, is proportional to

$\displaystyle
p(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\theta},\boldsymbol{\lambda},{\bf
T}\mid\mathbf{z})$ $\displaystyle\propto$
$\displaystyle\left(\prod_{i=1}^{n}\prod_{j=1}^{n-i+1}f_{\mathcal{N}}(z_{ij}\mid\boldsymbol{\theta},\lambda_{ij},T_{ij})f_{\mathcal{N}\mathcal{T}}(T_{ij}\mid\lambda_{ij})p(\lambda_{ij}\mid\nu)\right)$
(11) $\displaystyle\times$
$\displaystyle\left(\prod_{i=2}^{n}p(\alpha_{i}\mid\alpha_{i-1})\prod_{i=2}^{n}\prod_{j=2}^{n-i+1}p(\beta_{ij}\mid\beta_{i-1,j})\prod_{t=2}^{n}p(\gamma_{t}\mid\gamma_{t-1})\right)\times\pi(\boldsymbol{\theta})$
$\displaystyle\propto$
$\displaystyle\prod_{i=1}^{n}\prod_{j=1}^{n-i+1}\left\\{\left[\frac{\lambda_{ij}}{\sigma^{2}\left(1-\rho^{2}\right)}\right]^{1/2}\exp\left\\{-\frac{\lambda_{ij}}{2}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{\sigma^{2}\left(1-\rho^{2}\right)}\right]\right\\}\right\\}$
$\displaystyle\times$
$\displaystyle\prod_{i=1}^{n}\prod_{j=1}^{n-i+1}\left\\{\left[\frac{\lambda_{ij}}{\sigma^{2}}\right]^{1/2}\exp\left\\{-\frac{\lambda_{ij}}{2\sigma^{2}}T_{ij}^{2}\right\\}\boldsymbol{I}_{(0,\infty)}(T_{ij})p\left(\lambda_{ij}|\nu\right)\right\\}$
$\displaystyle\times$
$\displaystyle\left(\sigma_{\alpha}^{2}\right)^{-(n-1)/2}\exp\left\\{-\frac{1}{2\sigma_{\alpha}^{2}}\sum_{i=2}^{n}\left(\alpha_{i}-\alpha_{i-1}\right)^{2}\right\\}$
$\displaystyle\times$
$\displaystyle\left(\sigma_{\beta}^{2}\right)^{-(n-2)(n-1)/4}\exp\left\\{-\frac{1}{2\sigma_{\beta}^{2}}\sum_{i=2}^{n-1}\sum_{j=2}^{n-i+1}\left(\beta_{ij}-\beta_{i-1,j}\right)^{2}\right\\}$
$\displaystyle\times$
$\displaystyle\left(\sigma_{\gamma}^{2}\right)^{-(n-1)/2}\exp\left\\{-\frac{1}{2\sigma_{\gamma}^{2}}\sum_{t=2}^{n}\left(\gamma_{t}-\gamma_{t-1}\right)^{2}\right\\}$
$\displaystyle\times$
$\displaystyle\exp\left\\{-\frac{\mu^{2}}{2s_{\mu}^{2}}\right\\}\left(\sigma^{2}\right)^{-\left(a_{\sigma}+1\right)}\exp\left\\{-\frac{b_{\sigma}}{\sigma^{2}}\right\\}(1+\rho)^{c-1}(1-\rho)^{d-1}\boldsymbol{I}_{(-1,1)}(\rho)\nu^{a_{l}-1}\exp\left\\{-b_{l}\nu\right\\}$
$\displaystyle\times$
$\displaystyle\left(\sigma_{\alpha}^{2}\right)^{-(a_{\alpha}+1)}\exp\left\\{-\frac{b_{\alpha}}{\sigma_{\alpha}^{2}}\right\\}\left(\sigma_{\beta}^{2}\right)^{-(a_{\beta}+1)}\exp\left\\{-\frac{b_{\beta}}{\sigma_{\beta}^{2}}\right\\}\left(\sigma_{\gamma}^{2}\right)^{-\left(a_{\gamma}+1\right)}\exp\left\\{-\frac{b_{\gamma}}{\sigma_{\gamma}^{2}}\right\\},$

where $f_{\mathcal{N}}(\cdot)$ denotes the density of the Gaussian
distribution, $f_{\mathcal{N}\mathcal{T}}(\cdot)$ the density of the Truncated
Gaussian distribution, $\pi(\cdot)$ stands for the prior distribution of the
static parameters in $\boldsymbol{\theta}$, and $\boldsymbol{I}_{A}(\cdot)$ is
a characteristic function of the subset $A$ and $l\in\\{t,s,vg\\}$. See that
the resulting posterior distribution does not have an analytical solution, and
we appeal to Markov chain Monte Carlo techniques (MCMC) (Gamerman and Lopes,,
2006) to obtain samples from this distribution. The inference procedure is
detailed in Appendix B.

###  2.3 Computational issues

Conditional on the latent variables $\mathbf{T}$ and $\boldsymbol{\lambda}$
the Gaussianity is preserved and hence simple to build an algorithm to
estimate the model in Section 2.1. Below, we describe the sampling scheme from
the posterior full conditional distributions through the Gibbs algorithm of
the proposed model.

1 Initialize $\boldsymbol{\theta}^{(0)}$, $\boldsymbol{\alpha}^{(0)}$,
$\boldsymbol{\beta}^{(0)}$, $\boldsymbol{\gamma}^{(0)}$, $\mathbf{T}^{(0)}$
and $\boldsymbol{\lambda}^{(0)}$

2 for  _$l=1,\dots,M$_ do

3 $\boldsymbol{\alpha}^{(l)}\sim
p\left(\boldsymbol{\alpha}\mid\boldsymbol{\alpha}^{(l-1)},\boldsymbol{\beta}^{(l-1)},\boldsymbol{\gamma}^{(l-1)},\boldsymbol{\theta}^{(l-1)},\mathbf{T}^{(l-1)},\boldsymbol{\lambda}^{(l-1)},\mathbf{z}\right)$

4 $\boldsymbol{\beta}^{(l)}\sim
p\left(\boldsymbol{\beta}\mid\boldsymbol{\alpha}^{(l)},\boldsymbol{\beta}^{(l-1)},\boldsymbol{\gamma}^{(l-1)},\boldsymbol{\theta}^{(l-1)},\mathbf{T}^{(l-1)},\boldsymbol{\lambda}^{(l-1)},\mathbf{z}\right)$

5 $\boldsymbol{\gamma}^{(l)}\sim
p\left(\boldsymbol{\gamma}\mid\boldsymbol{\alpha}^{(l)},\boldsymbol{\beta}^{(l)},\boldsymbol{\gamma}^{(l-1)},\boldsymbol{\theta}^{(l-1)},\mathbf{T}^{(l-1)},\boldsymbol{\lambda}^{(l-1)},\mathbf{z}\right)$

6 $\boldsymbol{\theta}^{(l)}\sim
p\left(\boldsymbol{\theta}\mid\boldsymbol{\alpha}^{(l)},\boldsymbol{\beta}^{(l)},\boldsymbol{\gamma}^{(l)},\boldsymbol{\theta}^{(l-1)},\mathbf{T}^{(l-1)},\boldsymbol{\lambda}^{(l-1)},\mathbf{z}\right)$

7 $\mathbf{T}^{(l)}\sim
p\left(\mathbf{T}\mid\boldsymbol{\alpha}^{(l)},\boldsymbol{\beta}^{(l)},\boldsymbol{\gamma}^{(l)},\boldsymbol{\theta}^{(l)},\mathbf{T}^{(l-1)},\boldsymbol{\lambda}^{(l-1)},\mathbf{z}\right)$

8 $\boldsymbol{\lambda}^{(l)}\sim
p\left(\boldsymbol{\lambda}\mid\boldsymbol{\alpha}^{(l)},\boldsymbol{\beta}^{(l)},\boldsymbol{\gamma}^{(l)},\boldsymbol{\theta}^{(l)},\mathbf{T}^{(l)},\boldsymbol{\lambda}^{(l-1)},\mathbf{z}\right)$

9

10 end for

Algorithm 1 Gibbs algorithm with adaptive Metropolis step.

For the parameter vector $\boldsymbol{\theta}$ whose complete conditional
distribution is unknown ($\nu$ and $\rho$), we consider an additional step
based on the adaptive Metropolis-Hastings algorithm (Andrieu and Thoms,
(2008)). We employ the same procedure for the mixing parameter
$\boldsymbol{\lambda}$.

More details about the complete conditional distributions of the parameters,
latent effects and mixing variables are described in Appendix B. The algorithm
was coded in R using RStudio Version 2022.12.0-353 (Team, (2021)). Due to the
high autocorrelation of the samples generated from the latent effects and the
unknown parameters obtained through the adaptive method, the algorithm will
require many iterations to generate samples from the posterior distribution
$p(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\theta},\boldsymbol{\lambda},{\bf
T}\mid\mathbf{z})$. For this, we use the Rcpp library to build improved loops
of complete conditional distributions by performing part of the code
implementation with the C++ programming language and with the open source C++
library for statistical computing, armadillo-11.4.3, available at the website
https://arma.sourceforge.net/.

###  2.4 Predictions and Model Comparison

In order to check the predictive accuracy of the proposed model, consider a
vector $(\mathbf{z}^{obs},\mathbf{z}^{pred})$, with $(n+h)$ dimensional where
$\mathbf{z}^{obs}$ and $\mathbf{z}^{pred}$ representing, respectively, the
$n$-th dimensional observed and the $h$-th dimensional out-of-sample values of
$\mathbf{z}$. Here, the vector $\mathbf{z}^{obs}$ represents the information
available in the upper triangle unless the observations are left out in the
inference procedure, represented by $\mathbf{z}^{pred}$. The idea is to verify
if the proposed model is able to provide good performance to predict the lower
triangle.

To obtain samples from the posterior predictive distribution
$p(\mathbf{z}^{pred}\mid\mathbf{z}^{obs})$, we appeal to composition sampling
where assuming that
$\boldsymbol{\theta}^{(k)},\boldsymbol{\alpha}^{(k)},\boldsymbol{\beta}^{(k)},\boldsymbol{\gamma}^{(k)},\boldsymbol{\lambda}^{(k)},\mathbf{T}^{(k)}$
represent the $k$-th sample from the joint posterior distribution
$p(\boldsymbol{\theta},\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\lambda}^{obs},\mathbf{T}^{obs}\mid\mathbf{z}^{obs})$.
Thus, samples from $p(\mathbf{z}^{pred}\mid\mathbf{z}^{obs})$ are obtained by
computing from

  * (i)

$\boldsymbol{\lambda}^{pred}\mid\boldsymbol{\lambda}^{(k)},\nu^{(k)}$ and

  * (ii)

$\mathbf{T}^{pred}\mid\mathbf{T}^{(k)},\boldsymbol{\lambda}^{pred}$ and

  * (iii)

$\mathbf{z}^{pred}\mid\mathbf{z}^{obs},\boldsymbol{\lambda}^{pred},\mathbf{T}^{pred},\boldsymbol{\theta}^{(k)},\boldsymbol{\alpha}^{(k)},\boldsymbol{\beta}^{(k)},\boldsymbol{\gamma}^{(k)}$.

Notice that (i) represents the distribution responsible for accommodating the
tails of the distribution as seen previously in Section 2.1. The distributions
in (ii) and (iii) are the Truncated Gaussian and the Gaussian, respectively.
The same procedure is considered to obtain the lower triangle.

To model comparison, we consider measures based on scoring rules and the root
mean square prediction error (RMSPE) assessing the predictive performance of
the competing models. Scoring rules provide summary measures for the
evaluation of probabilistic forecasts, by assigning a numerical score based on
the predictive distribution and on the event or value that materializes
(Gneiting and Raftery,, 2007). In this case, we consider the Interval Score
(IS), the Width of Credible Interval (CWI), and the Continuous Ranked
Probability Score (CRPS). A brief description is presented in Appendix C.

###  2.5 Bayesian residual analysis

Consider $Z_{ij}$ observations at policy year $i$ with development $j$ as
defined in equation (2.2) such that $Z_{ij}\mid
T_{ij},\lambda_{ij},\boldsymbol{\Theta}\sim\mathcal{N}\left(\mu_{ij}+\rho
T_{ij},\lambda_{ij}^{-1}\sigma^{2}\left(1-\rho^{2}\right)\right)$. Then the
standardized Bayesian residual for the mixture model is

$r_{ij}=\frac{Z_{ij}-E\left(Z_{ij}\mid\boldsymbol{\Theta},\lambda_{ij},T_{ij}\right)}{\sqrt{Var\left(Z_{ij}\mid\boldsymbol{\Theta},\lambda_{ij},T_{ij}\right)}},$
(12)

with $i\leq n$ and $j\leq n-i+1$ . If the errors have Gaussian distribution,
then approximately 95% of the individual residuals are expected to be in the
interval $[-2,2]$. Otherwise, there is some evidence that this observation
could be a potential outlier.

##  3 Applications

This section presents two applications to the discussion about loss reserving
in insurance. A realistic simulated dataset that emulates a real behaviour
present in loss-reserving data to evaluate the performance of our proposed
methodology; and a case study that considers the claim amounts paid to the
insureds of an insurance company from 1978 to 1995. These data have been
previously analysed by Chan et al., (2008) and Choy et al., (2016). The idea
of the case study is to illustrate the flexibility of our proposed model which
is able to account jointly skewness and heavy tails when compared to well-
known models used in the literature of loss reserving. Codes with the
simulated study and inference procedure are available at
https://github.com/williamleao/SHTLR.

We fit the models described in Table 2, which are particular cases of the
general model proposed in Section 2. To verify how our proposed models are
able to predict throughout the lower triangle in the case study, we consider a
cross-validation predictive evaluation where a part of the upper triangle was
used for training, and the remaining was left out of the inference procedure
for validation.

Table 2: Competing models fitted to data applications: Normal ($\mathcal{N}$), Student-t ($\mathcal{ST}$), Slash ($\mathcal{S}$), Variance-Gamma ($\mathcal{VG}$), Skew-Normal ($\mathcal{SN}$), Skew-Student-t ($\mathcal{SST}$), Skew-Slash ($\mathcal{SS}$), Skew-Variance-Gamma ($\mathcal{SVG}$).  class | model | $\lambda$ | $T$
---|---|---|---
| $\mathcal{N}$ | $1$ | -
non-Skew | $\mathcal{ST}$ | $\mathcal{G}(\nu/2,\nu/2)$ | -
| $\mathcal{S}$ | $\mathcal{B}e(\nu,1)$ | -
| $\mathcal{VG}$ | $\mathcal{I}\mathcal{G}(\nu/2,\nu/2)$ | -
| $\mathcal{SN}$ | 1 | $\mathcal{N}\mathcal{T}_{[0,\infty)}(0,\sigma^{2})$
Skew | $\mathcal{SST}$ | $\mathcal{G}(\nu/2,\nu/2)$ | $\mathcal{N}\mathcal{T}_{[0,\infty)}(0,\sigma^{2}\lambda^{-1})$
| $\mathcal{SS}$ | $\mathcal{B}e(\nu,1)$ | $\mathcal{N}\mathcal{T}_{[0,\infty)}(0,\sigma^{2}\lambda^{-1})$
| $\mathcal{SVG}$ | $\mathcal{I}\mathcal{G}(\nu/2,\nu/2)$ | $\mathcal{N}\mathcal{T}_{[0,\infty)}(0,\sigma^{2}\lambda^{-1})$

###  3.1 A simulated Skew-Heavy-tailed dataset

This section presents a simulated study to investigate the convergence of the
parameters and predictive performance of the competing models in identifying
skewness and heavy tails. For this, we consider simulating an artificial run-
off triangle with $16\times 16$ (totalling $n=136$ observations in the upper
triangle) from a Skew Student-t distribution, assuming that the mean function
$\mu_{ij}=\mu+\alpha_{i}+\beta_{ij}+\gamma_{t}$, with $\mu=9$,
$\alpha_{i}\mid\alpha_{i-1}\sim\mathcal{N}(0,\sigma^{2}_{\alpha}=0.13)$,
$\beta_{ij}\mid\beta_{i-1j}\sim\mathcal{N}(0,\sigma^{2}_{\beta}=0.05)$,
$\gamma_{t}\mid\gamma_{t-1}\sim\mathcal{N}(0,\sigma^{2}_{\gamma}=0.13)$,
$i=2,\ldots,n$, $j=2,\ldots,n-i+1$, $t=i+j-1$, and parameters related on
variability, kurtosis and asymmetry as $\nu=3$, $\rho=-0.89$ and
$\sigma^{2}=0.14$, respectively, which were defined based on the posterior
point estimates obtained for the parameters in the dataset from Choy et al.,
(2016) . We are interested in predicting the total reserve (light grey cells
as shown in Table 1, which means $Y_{\nabla}$). The true total loss reserving
given by the lower triangle is 33,047.39 currency units. Priors on parameters
are the same as shown in Section 2.2. The burn-in and lag for spacing of the
chain were selected so that the effective sample size was around 5,000
samples. To model comparison, we consider the use of scoring rules and
evaluate parameters responsible for kurtosis and asymmetry and compare four
space-state competing models: $\mathcal{N}$, $\mathcal{ST}$, $\mathcal{SN}$
and $\mathcal{SST}$, that last representing the data-generating model. Also,
the usual chain ladder method was considered to project the loss reserving as
a naive technique.

Figure 2 shows the behaviour of the simulated dataset in the log scale (see
panel (a)) indicating skewness and ticker tails. See that the behaviour of the
claims is very irregular and contains some extreme outliers. Panels (b) and
(c) of Figure 2 show the trend of effects development year and accident year
in terms of the empirical mean, respectively. As we can see, in panel (b) the
average of these trends decreases over development years. To the effect of the
average accident year (panel (c)), we notice a volatility trend over the
years.

Posterior summaries are addressed for the competing models. The parameters to
be estimated are the dynamic coefficients ($\alpha_{i}$, $\beta_{ij}$,
$\gamma_{t})$), the variance parameters ($\sigma^{2}_{\alpha}$,
$\sigma^{2}_{\beta}$, $\sigma^{2}_{\gamma}$, $\sigma^{2}$), the function of
the skew parameter $\rho$, the mixing parameter $\nu$ and the latent mixing
variables ($\lambda_{ij}$, $T_{ij}$), $i=1,\dots,n$ and $j=1,\ldots,n$. Figure
3 presents the posterior distribution for the dynamic parameters $\alpha_{i}$,
$\beta_{ij}$ and $\gamma_{t}$ compared with the true values (red dot),
indicating that the estimates are well recovered by the Skew-t model. We have
already expected this result since the data are generated from a Skew-t
process. The first panel (a) shows a decreasing trend over accident years with
low volatility, whereas, in panel (c), the state parameter $\beta_{ij}$ shows
a trend of increase of the development year during the first two years and
then in the third, the development year decrease abruptly over accident years.
Notice that the third lag-year effect is low due to the fact the value of
observations decreases over accident years. After, they turn to increase and
maintain less volatility over the years. For the calendar year parameter, see
an increasing trend over time and is relatively stable over in the first years
(see panel (b)).

|  |
---|---|---
(a) log-claims distribution | (b) development year | (c) accident year
Figure 2: Simulated dataset: (a) histogram of the log claims, (b) the empirical mean of the development lag year trend, and (c) the empirical mean of the accident year. |
---|---
(a) | (b)
(c)
Figure 3: Simulated dataset: posterior summaries for the dynamic mean effects $\alpha_{i}$, $\gamma_{t}$ and $\beta_{ij}$ under the Skew-t model. The box-plot represent the distribution of the parameters and the red dots represent the true values generated from the Skew-t process. |
---|---
(a) $\sigma^{2}/\lambda_{ij}$ | (b) $\sigma^{2}(1-\rho^{2})/\lambda_{ij}$
Figure 4: Simulated dataset: posterior median with 95% posterior credible
intervals for the variance estimated structure assuming (a) the non-skew
models $\sigma^{2}/\lambda_{ij}$ and (b) the skew models
$\sigma^{2}(1-\rho^{2})/\lambda_{ij}$. The red line represents the variance of
the Gaussian class.

In the most complete model given by equation (2.2), the mean of $Z_{ij}$
conditional on $\lambda_{ij}$ and $T_{ij}$ is
$\mu_{ij}=\mu+\alpha_{i}+\beta_{ij}+\gamma_{t}+\rho T_{ij}$ and the variance
is $\sigma^{2}_{ij}=\sigma^{2}(1-\rho^{2})/\lambda_{ij}$, from which we can
deduce the variance structures over development and accident years. Figure 4
illustrates the variance $\sigma^{2}_{ij}$ estimated for the Gaussian, the
Student-t, the Skew-Normal, and the Skew-Student-t models. Note that the
variance is non-constant with some peaks over the development periods. As we
can see, the class of the Student-t models can recover peaks of variability
and potential outliers. Although the Skew-Normal model is able to capture
changes in the skews of the distribution, the model does not take into account
volatility over the development periods, which means $\lambda_{ij}=1$ for all
$i=1,\ldots,n$ and $j=1,\ldots,n$ (see panel (b)). Notice that if we assume
$\lambda_{ij}=1$ and $\rho=0$ the posterior variance for an usual dynamical
Gaussian model is considered (see red line in panel (a)). The variance of the
non-skew class is higher than the skew class (see axis $y$ in Figure 4). Table
3 shows the posterior summary for the Skew-t model. See that, the median
posterior of $\rho$ is equal to -0.877 (with true value -0.890) indicating the
Skew-t model captures the presence of skewness in the data. Panels of Figure 5
present the 95% credible intervals posterior summaries of the mean $\mu_{ij}$
for the competing models compared with the true values (black dots) indicating
that the proposed models recover the true mean over the period. Although the
Gaussian model covers the true values, the model has a higher range of 95%
posterior credible intervals for $\mu_{ij}$ compared to the others.

For prediction purposes, we assess the amount of loss reserve for all
competing models. Figure 6 plots the medians and 95% credible intervals for
the predictions of loss reserve for the first calendar year. The one-year
ahead reserve could be important in short-term financial decisions for the
insurance company, as they are the reserves for settling that are expected to
occur soon. The class of the Skew models provides interval widths that are
smaller in the majority of first calendar year data when compared to the non-
skewed models. Figure 7 shows the predictive posterior total loss reserving
distribution for all competing models compared with the true total reserve and
chain ladder reserve. The uncertainty is worse for the non-skew class, mainly
for the Gaussian model. Besides that, the chain ladder method underestimates
the loss reserving producing an amount smaller than would be expected of
22,835.28 (dashed red line). Remember that the true total reserve is
33,047.39. In fact, the chain ladder technique does not include any calendar
year effects and is unsuitable for volatile data. In the actuarial context,
this would be an aggravating factor due to the fact that the insurance company
does not save enough to constitute an efficient reserve for soon payments. For
the dynamical linear models, the total reserve posterior median and 95%
credible interval are $Y_{\nabla\mathcal{N}}=47,993.08(18,354.42;155,179.7)$,
$Y_{\nabla\mathcal{T}}=35,407.88(11,387.48;181,571.9)$,
$Y_{\nabla\mathcal{SN}}=33,339.87(13,486.54;103,837.2)$,
$Y_{\nabla\mathcal{ST}}=34,339.96(12,106.10;146,763.2)$. Table 4 shows model
comparison based on RMSPE and scoring rules criteria indicating that the
Skew-t model provides better performance when compared to others. This
suggests that taking into account heterogeneous variance and skewness is
essential in order to produce good predictions in loss reserving context,
especially for extreme values.

---
(a) non-Skew models
(b) Skew models
Figure 5: Simulated dataset: posterior claims with 95% credible interval and true values (black dotted) for the dynamical mean $\mu_{ij}$, (a) the non-skew models and (b) the skew models.  |
---|---
(a) Normal | (b) Student-t
|
(c) Skew-Normal | (d) Skew-t
Figure 6: Simulated dataset: Posterior predictive medians (dots) versus the
true observed loss reserving for the first calendar year with vertical lines
representing the 95% credible predictive intervals. The diagonal line
indicates $y=x$. (a) the Normal, (b) the Student-t, (c) the Skew-Normal, and
(d) the Skew-t models.
---
Figure 7: Simulated dataset: predictive posterior of the total loss reserving for all competing models. The dashed grey line represents the true total loss reserve and the dashed red line represents the chain ladder reserve (CL). Table 3: Simulated dataset: Posterior summary for the Skew-t model compared with the true values. Parameters | True value | Median | S.E. | IC95%
---|---|---|---|---
$\mu$ | 9.00 | 8.815 | 0.284 | (8.328 , 9.478)
$\rho$ | -0.89 | -0.877 | 0.211 | (-0.996 , -0.187)
$\nu$ | 3.00 | 4.579 | 1.458 | (2.606 , 8.234)
$\sigma^{2}$ | 0.14 | 0.162 | 0.080 | (0.061 , 0.366)
$\sigma_{\alpha}^{2}$ | 0.13 | 0.249 | 0.139 | (0.105 , 0.638)
$\sigma_{\beta}^{2}$ | 0.05 | 0.060 | 0.021 | (0.030 , 0.112)
$\sigma_{\gamma}^{2}$ | 0.13 | 0.094 | 0.061 | (0.033 , 0.258)
Table 4: Simulated dataset: Model comparison based on the Root Mean Square Prediction Error (RMSPE), the Interval Score (IS), the Width of Credibility Interval (WCI), and the Continuous Ranking Predictive Score criteria for the predicted total loss reserve (the lower triangle) under all fitted models. Competing | Measures
---|---
Models | average IS | average WCI | RMSPE | average CRPS
$\mathcal{N}$ | 5.307 | 5.307 | 16.957 | 0.545
$\mathcal{ST}$ | 4.452 | 4.452 | 14.497 | 0.437
$\mathcal{SN}$ | 4.732 | 4.732 | 15.408 | 0.489
$\mathcal{SST}$ | 4.421 | 4.421 | 14.354 | 0.434

###  3.2 Case study: a real dataset from Chan et al. (2008)

We investigate the same dataset considered by Chan et al., (2008) and Choy et
al., (2016). The data refers to the claim amounts paid to the insureds of an
insurance company during the period from 1978 to 1995 (see Appendix D, Table
D.2). There are $n=171$ observations in the upper triangle. Given a policy
year, the amount of claims paid follows an increasing trend in the first 4 to
6 lag years and then a decreasing trend thereafter (Chan et al., (2008)). For
prediction evaluation, we consider a cross-validation proposal where $13\times
13$ observations were used for the inference procedure, being that the last
five diagonals (calendar year) we left out for validation (for more details,
see Appendix D \- highlighted in bold). Now, our data is composed of $n=91$
observations (upper triangle, $Y_{\Delta}$) and the remaining representing the
lower triangle (total reserve, $Y_{\nabla}$ highlighted in dark and light
grey).

The mean, median, standard deviation, skewness, and excess kurtosis of the
log-transformed data log($Y_{ij}$) are 8.214 8.450 1.079 -1.759 and 3.767,
respectively. The behaviour of the claims is irregular and contains potential
outliers. Panel (a) of Figure 8 presents the distribution of the log claims
suggesting that skewness and heavy tails may be accommodated in the modelling
process.

|  |
---|---|---
(a) log-claims data | (b) residual empirical median | (c) q-q plot of residual
Figure 8: Case study: (a) the histogram of log-claims distribution (b) the
residual median of the Gaussian model and (c) the Normal q-q plot of residuals
for the dataset.

Initially, we fit a Gaussian model, which means $\lambda_{ij}=1$, for all
$i,j=1,\ldots,n$ and $\kappa=0$. Panels (b) and (c) of Figure 8 show the
residual empirical median and the theoretical quantiles for the Gaussian
distribution versus the empirical quantiles based on the residuals from the
Gaussian model suggesting that the model with constant variance and without
asymmetry might be unsuitable to predict the loss reserving. Note that the
empirical excess kurtosis and skewness of the residuals are 3.767 and -1.759,
respectively. These features corroborate that we need to consider a model that
is able to capture this behaviour. Given the results from Figure 8, we fit all
competing models described in Table 2. The parameters to be estimated are
related to dynamical mean $\mu_{ij}$ that is
$(\mu,\alpha_{i},\beta_{ij},\gamma_{t})$, the variance parameters
$(\sigma^{2},\sigma^{2}_{\alpha},\sigma^{2}_{\beta},\sigma^{2}_{\gamma})$, the
mixing parameters $(\nu,\rho)$, and the latent variables
$(\lambda_{ij},T_{ij})$.

Sensitivity to the prior parameter specification is an important step in the
inference procedure. An unsuitable specification of the prior parameter can
lead a misleading results. For this, we investigate the performance of the
parameters based on a sensitivity analysis by considering different
specifications of the prior distributions. The proposed priors are presented
in Table 5 and the choice of the prior distributions are detailed in Gelman,
(2006), Goudarzi and Zokaei, (2021), and Pérez-Rodríguez et al., (2018) for
the Skew-t model. The prior sensitivity analysis for the Skew-t model is shown
in Table 6 through the posterior median, 95% credible interval and the
predicted posterior total reserve. See that there are no significant
differences for the specifications of $\mu$, $\rho$, $\sigma^{2}$,
$\sigma_{\alpha}^{2}$, $\sigma_{\beta}^{2}$ and $\sigma_{\gamma}^{2}$
distributions. Although there are significant differences in the choice of the
$\nu$ parameter ($\mathcal{G}(12,0.8)$, $Exp(0.3)$, $\mathcal{U}(2,40)$ and
Jeffreys prior distribution proposed by Fonseca et al., (2008)), there is no
impact on the total reserve for several scenarios. The following results shown
in this study are based on scenario 2.

Table 5: Case study: Prior distributions scenario for the Skew-t model. | Prior distributions
---|---
scenario | $\mu$ | $\rho$ | $\nu$ |  $\sigma^{2},\sigma_{\alpha}^{2},\sigma_{\beta}^{2}$ and $\sigma_{\gamma}^{2}$
1 | $\mathcal{N}(0,100)$ | $\mathcal{B}e(1,1)$ | $\mathcal{G}(12,0.8)$ | $\mathcal{G}(0.01,0.01)$
2 | $\mathcal{N}(0,100)$ | $\mathcal{B}e(1,1)$ | $\mathcal{G}(12,0.8)$ | $\mathcal{G}(0.001,0.001)$
3 | $\mathcal{N}(0,100)$ | $\mathcal{B}(0.5,0.5)$ | $\mathcal{G}(12,0.8)$ | $\mathcal{G}(0.001,0.001)$
4 | $\mathcal{N}(0,100)$ | $\mathcal{B}e(1,1)$ | $\mathcal{U}(2,40)$ | $\mathcal{G}(0.001,0.001)$
5 | $\mathcal{N}(0,100)$ | $\mathcal{B}e(1,1)$ | $Exp(0.3)$ | $\mathcal{G}(0.001,0.001)$
6 | $\mathcal{N}(0,100)$ | $\mathcal{B}e(1,1)$ | Jeffreys | $\mathcal{G}(0.001,0.001)$
Table 6: Case study: Prior sensitivity analysis for the Skew-t model.

Parameters scenario Reserve $\mu$ $\rho$ $\sigma^{2}$ $\nu$
$\sigma_{\alpha}^{2}$ $\sigma_{\beta}^{2}$ $\sigma_{\gamma}^{2}$ 1 630165.9
8.891 -0.960 0.278 11.187 0.113 0.061 0.131 (8.229,9.624) (-0.998,-0.337)
(0.111,0.545) (5.175,21.119) (0.039,0.365) (0.028,0.126) (0.053,0.389) 2
621282.4 8.921 -0.959 0.282 11.150 0.112 0.059 0.129 (8.206,9.679)
(-0.999,-0.455) (0.112,0.556) (4.958,21.064) (0.038,0.342) (0.027,0.122)
(0.050,0.391) 3 613902.6 8.955 -0.989 0.323 11.462 0.103 0.058 0.127
(8.246,9.743) (-1.000,-0.794) (0.140,0.607) (5.207,21.401) (0.035,0.324)
(0.027,0.120) (0.049,0.388) 4 591693.9 8.773 -0.896 0.106 2.947 0.091 0.058
0.118 (8.132,9.566) (-0.996,-0.216) (0.034,0.353) (2.035,11.827) (0.030,0.289)
(0.027,0.117) (0.046,0.374) 5 586052.7 8.734 -0.875 0.091 2.655 0.087 0.058
0.114 (8.122,9.557) (-0.994,-0.209) (0.031,0.264) (2.025,5.890) (0.029,0.279)
(0.028,0.116) (0.044,0.348) 6 593969.0 8.743 -0.876 0.088 2.558 0.086 0.058
0.115 (8.119,9.590) (-0.995,-0.181) (0.031,0.273) (2.021,5.636) (0.027,0.275)
(0.027,0.121) (0.042,0.346)

Posterior summaries for all competing models are provided by Figures 9, 10,
and 11. Panels in Figure 9 show the posterior median and the limits of the 95%
credible interval (in grey) compared to the observed values. The skew class
exhibits better performance, effectively capturing the impact of atypical
observations in the log claims when compared to the non-skew class.

Figure 10 presents the behaviour of the estimated standard deviation for the
proposed skew class. It is possible to identify four points (highlighted in
dashed red line), $Y_{1,10}=333$, $Y_{2,10}=156$, $Y_{1,11}=199$, and
$Y_{2,11}=35$, common to all panels, where models with heavy tails recognize
these observations as potential outliers. The Skew-Normal model inflates the
variability for all observations to adjust the standard deviation, while
models that take into account heavy tails are able to identify atypical
observations and precisely increase the necessary variability of a given log-
claim. Figure 11 presents the box-plots of the effects of accident year
$\alpha_{i}$, calendar year $\gamma_{t}$ and the interaction between accident
and development years $\beta_{ij}$. The dynamics of the three effects show
significant behaviour, it can be seen that the values of $\alpha_{i}$’s show
an increasing trend, the values of $\gamma_{t}$’s show a decreasing trend and
the values of $\beta_{ij}$’s show increasing, decreasing and even almost
constant variations for a given year of development.

|  |  |
---|---|---|---
(a) Normal | (b) Student-t | (c) Slash | (d) Variance Gamma
|  |  |
(e) Skew Normal | (f) Skew-t | (g) Skew Slash | (h) Skew Variance Gamma
Figure 9: Case study: Posterior fit value medians (red line) with 95% credible
interval (grey) versus the true observed (black dots) $Y_{\Delta}$. (a) the
Normal model, (b) the Student-t model, (c) the Slash model, (d) the Variance-
Gamma model, (e) the Skew-Normal model, (f) the Skew-t model, (g) the Skew
Slash model and (h) the Skew Variance-Gamma model.

Figure 10: Case study: Posterior smoothed mean of observational standard
deviation $\lambda_{ij}$ with Skew Student t distribution $(\mathcal{SVG})$
(top panel), Skew Slash distribution $(\mathcal{SS})$ (center panel) and Skew
Variance Gamma distribution $(\mathcal{SVG})$ (bottom panel). The vertical red
dotted lines indicates the points $Y_{1,10}=333$, $Y_{2,10}=156$,
$Y_{1,11}=199$ and $Y_{2,11}=35$. The horizontal dotted line indicates the
observational standard deviation estimate of the Skew Normal model.

Figure 11: Case study: box-plot for the posterior distribution of the accident
year ($\alpha_{i}$), calendar year effect ($\gamma_{t}$) and $\beta_{ij}$
(development year) for the Skew-t model.

In order to check the presence of outliers in the dataset, we consider
computing the Bayesian standardized residuals (see Section 2.5) for all
competing models. The panels in Figure 12 display the residuals, indicating
that there are no outliers for the skew class, whereas outliers are evident in
the non-skew models, being unsuitable to accommodate this behaviour of the
data.

|  |  |
---|---|---|---
(a) Normal | (b) Student-t | (c) Slash | (d) Variance Gamma
|  |  |
(e) Skew Normal | (f) Skew-t | (g) Skew Slash | (h) Skew Variance Gamma
Figure 12: Case study: the standardized residuals for all competing models.

For predictive purposes, panels of Figures 13 and 14 show the summary of the
posterior predictive distribution obtained for the calendar years (1991 and
1995). In general, the class of skew models provided interval widths smaller
them the non-skewed models, resulting in the prediction of outstanding claim
amounts with higher accuracy. The RMSPE measure, the average IS, the CRPS, and
the CWI are shown in Table 7, where as expected, the Gaussian model is the one
with the worst predictive performance. As mentioned previously, this model is
not adequate to accommodate irregular claims and extreme outliers. Although we
are taking into account models that capture thicker tails than Normal, the
skewness is not allowed. Notice that models that regard this characteristic
perform better than the non-skewed models (see ranking in Table 7). The
measures indicate that the Skew-t has the better performance.

Table 7: Case study: Model comparison based on the Width of Credibility Interval (WCI),Root Mean Square Prediction Error (RMSPE) and the Continuous Ranking Predictive Score criteria for the predicted total loss reserve (the lower triangle) under all fitted models. Competing | Measures | Ranking
---|---|---
Models | average IS | average WCI | RMSPE | average CRPS
$\mathcal{N}$ | 5.584 | 4.515 | 1.513 | 0.618 | 8
$\mathcal{ST}$ | 5.405 | 4.137 | 1.404 | 0.578 | 5
$\mathcal{S}$ | 5.527 | 4.354 | 1.480 | 0.595 | 7
$\mathcal{VG}$ | 5.451 | 4.228 | 1.431 | 0.587 | 6
$\mathcal{SN}$ | 5.289 | 4.019 | 1.378 | 0.575 | 3
$\mathcal{SST}$ | 5.232 | 3.743 | 1.329 | 0.564 | 1
$\mathcal{SS}$ | 5.316 | 4.012 | 1.612 | 0.617 | 4
$\mathcal{SVG}$ | 5.284 | 3.856 | 1.351 | 0.575 | 2

Table 8 presents the posterior summaries derived from the Skew-t model for
predicting claim payments, as per comparison measures. To aid in diagnosing
the convergence of the chains, we employ the convergence diagnostic (CD) based
on the range of -1.96 to +1.96, frequently linked with the 95% confidence
interval for a standard normal distribution. We can observe that both chains
appear to display convergence behaviour. The parameter responsible for
skewness, $\rho$ is significant and equals -0.909 indicating that the model
can capture the asymmetry in the dataset. Besides that, the posterior mean for
$\nu$ equals 11.614 shows that the proposed model presents a heavier tail than
the Gaussian model, indicating a relatively higher probability of observing
extreme values.

Table 8: Case study: Posterior summary for the Skew-t model. Parameter | Mean | Std. Dev. | IC(95%) | CD
---|---|---|---|---
$\mu$ | 8.926 | 0.379 | $\left(8.206,9.679\right)$ | -0.606
$\rho$ | -0.909 | 0.161 | $\left(-0.999,-0.455\right)$ | 0.035
$\sigma^{2}$ | 0.294 | 0.114 | $\left(0.112,0.556\right)$ | -0.807
$\nu$ | 11.614 | 4.159 | $\left(4.958,21.064\right)$ | -0.578
$\sigma_{\alpha}^{2}$ | 0.131 | 0.085 | $\left(0.038,0.342\right)$ | 0.776
$\sigma_{\beta}^{2}$ | 0.063 | 0.025 | $\left(0.027,0.122\right)$ | -0.393
$\sigma_{\gamma}^{2}$ | 0.151 | 0.090 | $\left(0.050,0.391\right)$ | -0.110
|  |  |
---|---|---|---
(a) Normal | (b) Student-t | (c) Slash | (d) Variance Gamma
|  |  |
(e) Skew Normal | (f) Skew-t | (g) Skew Slash | (h) Skew Variance Gamma
Figure 13: Case study: Posterior predictive medians (dots) versus the true observed loss reserving for the calendar year 1991 with vertical lines representing the 95% credible predictive intervals. The diagonal line indicates $y=x$. (a) the Normal model, (b) the Student-t model, (c) the Slash model, (d) the Variance-Gamma model, (e) the Skew-Normal model, (f) the Skew-t model, (g) the Skew Slash model and (h) the Skew Variance-Gamma model.  |  |  |
---|---|---|---
(a) Normal | (b) Student-t | (c) Slash | (d) Variance Gamma
|  |  |
(e) Skew Normal | (f) Skew-t | (g) Skew Slash | (h) Skew Variance Gamma
Figure 14: Case study: Posterior predictive medians (dots) versus the true
observed loss reserving for the calendar year 1995 with vertical lines
representing the 95% credible predictive intervals. The diagonal line
indicates $y=x$. (a) the Normal model, (b) the Student-t model, (c) the Slash
model, (d) the Variance-Gamma model, (e) the Skew-Normal model, (f) the Skew-t
model, (g) the Skew Slash model and (h) the Skew Variance-Gamma model.

Figure 15 displays box-plots and the posterior predictive distribution of
total loss reserves for non-skewed and skewed models, in comparison to the
chain ladder method. The plots indicate that the predictive distribution of
total loss reserves for the non-skewed models is wider than for the skewed
models. It is worth noting that the validation process excluded the last five
years (1991-1995). Note that during this validation, a portion of the true
reserve for the triangle $Y_{\nabla}$ is known, with a value of 191,274.00.
This finding suggests that the chain ladder method underestimates the total
loss reserve, providing a total reserve of 123,776.90. Table 9 presents the
predictive posterior total reserve for quantiles 20%, 35%, 50%, 65%, and 80%
for all competing models. It is observed that higher quantiles (above 65%)
result in some overestimation in loss reserving, whereas lower quantiles
(below 50%) exhibit higher accuracy when compared to the chain ladder
technique.

Figure 15: Case study: predictive posterior total loss reserving for all competing models. The dashed red line represents the chain ladder loss reserve. Table 9: Case study: Predictive posterior quantiles of 20%, 35%, 50%, 65% and 80% of the total reserve for all competing models. Models | Predictive Posterior Quantiles
---|---
$20\%$ | $35\%$ | $50\%$ | $65\%$ | $80\%$
$\mathcal{N}$ | 415,528.3 | 625,046.3 | 887,277.5 | 1,281,023.6 | 2,095,537
$\mathcal{ST}$ | 374,646.8 | 534,794.5 | 735,967.6 | 1,018,925.9 | 1,609,357
$\mathcal{S}$ | 393,488.7 | 571,226.8 | 829,155.1 | 1,236,204.8 | 2,009,094
$\mathcal{VG}$ | 384,843.8 | 562,636.8 | 773,439.6 | 1,105,735.7 | 1,751,480
$\mathcal{SN}$ | 354,978.1 | 501,559.2 | 675,223.2 | 927,993.2 | 1,370,806
$\mathcal{SST}$ | 343,518.2 | 479,568.9 | 632,237.3 | 831,590.2 | 1,212,584
$\mathcal{SS}$ | 363,925.6 | 495,431.3 | 647,159.1 | 859,286.6 | 1,218,528
$\mathcal{SVG}$ | 346,078.8 | 475,412.1 | 638,935.6 | 872,867.0 | 1,288658

##  4 Conclusions

We have introduced a versatile framework for loss reserving that builds upon
the well-established space-state loss reserving modelling. Our framework
accommodates heavy tails and skewness in the run-off triangle structure
through a scale mixture of skew-normal distributions. The dynamic evolution in
the mean model proposed in equation (2.2) incorporates an additional term to
capture changes over calendar years, as previously suggested by Shi et al.,
(2012). Furthermore, the dynamic evolution in the variance model is capable of
accommodating various regimes of variability over time being a crucial aspect
in the context of loss-reserving modelling. As inference is performed under
the Bayesian paradigm using MCMC techniques allowing fast computations for
inference and prediction. It takes advantage of the conditionally skew-normal
distributions obtained when we condition the distribution of $Z_{ij}$ and
$T_{ij}$ on the mixing latent variables.

We performed a simulation study to investigate the ability of the proposed
model to capture asymmetry and fat tails. Our simulated example indicates that
the correct model is selected. The non-skew models have worse performance when
the data does have the presence of irregular claims and potential outliers and
the chain ladder technique fails to produce estimates of loss reserving
leading to underestimation. Moreover, it seems that the predictive scoring
rules used to compare the models are adequate measures of good predictive
performance. Also, we illustrate the performance of our proposed model through
a real dataset presented by Chan et al., (2008) and Choy et al., (2016) as a
case study. Our study showed that the skewed class has a better performance in
predicting the loss reserving when compared with non-skewed models.

We conclude that considering the class of the scale mixture of skew-normal
distribution produces robust inferences and predictions for the total reserve.
Therefore, an extension of this work involves contemplating a multivariate
loss reserving model, as in Goudarzi and Zokaei, (2021), but now featuring a
class of skew-heavy-tailed distributions.

##  Appendix A Proofs of properties of the proposed model

In this appendix, we prove the results shown in Section 2 and provide de
kurtosis and skewness of the proposed model. For this, consider the model
given by equation (1). We are interested in centers on the first four moments
to compute skewness and kurtosis. Remember that the variable
$\varepsilon_{ij}$ is a zero mean Skew-Normal distribution, with $\sigma^{2}$
and $\kappa$ skewness parameter, with
$\rho=\frac{\kappa}{\sqrt{1+\kappa^{2}}}$. Then,

$E(\varepsilon_{ij})=\rho\sqrt{\frac{2}{\pi}},\hskip
5.69046ptVar(\varepsilon_{ij})=E\left[\left(\varepsilon_{ij}-E(\varepsilon_{ij})\right)^{2}\right]=1-\frac{2}{\pi}\rho^{2}$
$E(\varepsilon_{ij}^{2})=1,\hskip
8.5359ptE(\varepsilon_{ij}^{3})=\sqrt{\frac{2}{\pi}}\left(3\rho-\rho^{3}\right),\hskip
8.5359ptE(\varepsilon_{ij}^{4})=3.$

The mean and variance of $Z_{ij}$ as mentioned in Section 2 is given by

$E(Z_{ij})=E\left(\mu_{ij}+\sigma\frac{\varepsilon_{ij}}{\lambda_{ij}}\right)=\mu_{ij}+\sigma\rho\sqrt{\frac{2}{\pi}}E\left(\lambda_{ij}^{-1/2}\right),$

and

$\displaystyle Var(Z_{ij})$ $\displaystyle=$ $\displaystyle
E\left[\left(Z_{ij}-E(Z_{ij})\right)^{2}\right]$ $\displaystyle=$
$\displaystyle\sigma^{2}\left(E\left(\lambda_{ij}^{-1}\varepsilon_{ij}^{2}\right)-E^{2}\left(\lambda_{ij}^{-1/2}\varepsilon_{ij}\right)\right)$
$\displaystyle=$
$\displaystyle\sigma^{2}\left(E\left(\lambda_{ij}^{-1}\right)E\left(\varepsilon_{ij}^{2}\right)-E^{2}\left(\lambda_{ij}^{-1/2}\varepsilon_{ij}\right)\right)$
$\displaystyle=$
$\displaystyle\sigma^{2}\left(E\left(\lambda_{ij}^{-1}\right)-E^{2}\left(\lambda_{ij}^{-1/2}\right)E^{2}\left(\varepsilon_{ij}\right)\right)$
$\displaystyle=$
$\displaystyle\sigma^{2}\left(E\left(\lambda_{ij}^{-1}\right)-\rho^{2}\left(\frac{2}{\pi}\right)E^{2}\left(\lambda_{ij}^{-1/2}\right)\right).$

The expression of the skewness, unconditional on the mixing parameter, is
given by

$\displaystyle Skew\left(Z_{ij}\right)$ $\displaystyle=$
$\displaystyle\frac{E\left[\left(Z_{ij}-E(Z_{ij})\right)^{3}\right]}{\left[Var(Z_{ij})\right]^{3/2}},$

where the third central moment
$\eta_{3}=E\left[\left(Z_{ij}-E(Z_{ij})\right)^{3}\right]$ is

$\displaystyle\eta_{3}=$ $\displaystyle
E\left[\left(\mu_{ij}+\sigma\varepsilon_{ij}\lambda_{ij}^{-1/2}-\mu_{ij}-\sigma\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{3}\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{3}E\left[\left(\varepsilon_{ij}\lambda_{ij}^{-1/2}-\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{3}\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{3}E\left[\left(\varepsilon_{ij}\lambda_{ij}^{-1/2}-\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{2}\left(\varepsilon_{ij}\lambda_{ij}^{-1/2}-\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{3}\left[E(\varepsilon_{ij}^{3})E(\lambda_{ij}^{-3/2})-3\rho\sqrt{\frac{2}{\pi}}E(\varepsilon_{ij}^{2})E(\lambda_{ij}^{-1})E\left[E(\lambda_{ij}^{-1/2})\right]+\frac{4}{\pi}\rho^{2}E(\varepsilon_{ij})E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{2}\right]\right.$
$\displaystyle+$
$\displaystyle\left.\frac{2}{\pi}\rho^{2}E(\varepsilon_{ij})E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{2}\right]-\left(\frac{2}{\pi}\right)^{3/2}\rho^{3}E\left[E(\lambda_{ij}^{-1/2})^{3}\right]\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{3}\left[\sqrt{\frac{2}{\pi}}\left(3\rho-\rho^{3}\right)E(\lambda_{ij}^{-3/2})-3\rho\sqrt{\frac{2}{\pi}}E(\lambda_{ij}^{-1})E\left[E(\lambda_{ij}^{-1/2})\right]+\frac{4}{\pi}\rho^{3}\sqrt{\frac{2}{\pi}}E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{2}\right]\right.$
$\displaystyle+$
$\displaystyle\left.\left(\frac{2}{\pi}\right){3/2}\rho^{3}E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{2}\right]-\left(\frac{2}{\pi}\right)^{3/2}\rho^{3}E\left[E(\lambda_{ij}^{-1/2})^{3}\right]\right].$

Notice that $\lambda_{ij}\rightarrow 1$, then
$\eta_{3}=\sigma^{3}\left[\sqrt{\frac{2}{\pi}}\rho^{3}\left(\frac{4-\pi}{\pi}\right)\right]$
represents a third central moment of a Skew-Normal.

The expression of the kurtosis, unconditional on the mixing parameter, is
given by

$Kurt(Z_{ij})=\frac{E\left[\left(Z_{ij}-E(Z_{ij})\right)^{4}\right]}{\left[Var(Z_{ij})\right]^{2}},$

where the fourth central moment
$\eta_{4}=E\left[\left(Z_{ij}-E(Z_{ij})\right)^{4}\right]$ is

$\displaystyle\eta_{4}=$ $\displaystyle
E\left[\left(\mu_{ij}+\sigma\varepsilon_{ij}\lambda_{ij}^{-1/2}-\mu_{ij}-\sigma\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{4}\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{4}E\left[\left(\varepsilon_{ij}\lambda_{ij}^{-1/2}-\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{4}\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{4}E\left[\left(\varepsilon_{ij}\lambda_{ij}^{-1/2}-\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{2}\left(\varepsilon_{ij}\lambda_{ij}^{-1/2}-\rho\sqrt{\frac{2}{\pi}}E\left(\lambda^{-1/2}\right)\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{4}\left[E(\varepsilon_{ij}^{4})E(\lambda_{ij}^{-2})-4\rho\sqrt{\frac{2}{\pi}}E(\varepsilon_{ij}^{3})E(\lambda_{ij}^{-3/2})E\left[E(\lambda_{ij}^{-1/2})\right]+\frac{12}{\pi}\rho^{2}E(\varepsilon_{ij}^{2})E(\lambda_{ij}^{-1})E\left[E(\lambda_{ij}^{-1/2})\right]\right.$
$\displaystyle-$
$\displaystyle\left.4\rho^{3}\left(\frac{2}{\pi}\right)^{3/2}E(\varepsilon_{ij})E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{3}\right]+\rho^{4}\frac{4}{\pi^{2}}E\left[E(\lambda_{ij}^{-1/2})^{4}\right]\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{4}\left[3E(\lambda_{ij}^{-2})-4\rho\sqrt{\frac{2}{\pi}}\sqrt{\frac{2}{\pi}}(3\rho-\rho^{3})E(\lambda_{ij}^{-3/2})E\left[E(\lambda_{ij}^{-1/2})\right]+\frac{12}{\pi}\rho^{2}E(\lambda_{ij}^{-1})E\left[E(\lambda_{ij}^{-1/2})\right]\right.$
$\displaystyle-$
$\displaystyle\left.4\rho^{3}\left(\frac{2}{\pi}\right)^{3/2}\rho\sqrt{\frac{2}{\pi}}E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{3}\right]-\rho^{4}\left(\frac{4}{\pi^{2}}\right)E\left[E(\lambda_{ij}^{-1/2})^{4}\right]\right]$
$\displaystyle=$
$\displaystyle\thinspace\sigma^{4}\left[3E(\lambda_{ij}^{-2})+\frac{12}{\pi}\rho^{2}E(\lambda_{ij}^{-1})E\left[E(\lambda_{ij}^{-1/2})\right]-\frac{16}{\pi^{2}}\rho^{4}E(\lambda_{ij}^{-1/2})E\left[E(\lambda_{ij}^{-1/2})^{3}\right]\right.$
$\displaystyle-$
$\displaystyle\left.\left(\frac{24}{\pi}\rho^{2}-\frac{8}{\pi}\rho^{4}\right)E(\lambda_{ij}^{-3/2})E\left[E(\lambda_{ij}^{-1/2})\right]+\frac{4}{\pi^{2}}\rho^{4}E\left[E(\lambda_{ij}^{-1/2})^{4}\right]\right].$

Notice that $\lambda_{ij}\rightarrow 1$, then
$\eta_{4}=3-\frac{12}{\pi}\rho^{2}+\left(\frac{8\pi-12}{\pi^{2}}\right)\rho^{4}$
represents a fourth central moment a Skew-Normal. The skewness and kurtosis
for the Skew-Gaussian model are

$Skew(Z_{ij})=\frac{\eta_{3}}{\eta_{2}^{3/2}}=\frac{\sqrt{2}(4-\pi)\rho^{3}}{(\pi-2\rho^{2})^{3/2}}$

and

$Kurt(Z_{ij})=\frac{\eta_{4}}{\eta_{2}^{2}}=3+\frac{8(\pi-3)\rho^{4}}{(\pi-2\rho^{2})^{2}}.$

For the class of heavy tails the skewness and kurtosis can be computed
considering the expected value for $\lambda_{ij}$ shown in Table A.1.

Table A.1:  Expected value for $\lambda_{ij}$ considering skew heavy tail models. | $\mathcal{SST}$ | $\mathcal{SS}$ | $\mathcal{SVG}$
---|---|---|---
$E(\lambda_{ij}^{-1/2})$ | $\left[{\Gamma(\frac{\nu-1}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{1/2}$ | $\frac{\nu}{\nu-1/2}$ | $\left[{\Gamma(\frac{\nu+1}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{-1/2}$
$E(\lambda_{ij}^{-3/2})$ | $\left[{\Gamma(\frac{\nu-3}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{3/2}$ | $\frac{\nu}{\nu-3/2}$ | $\left[{\Gamma(\frac{\nu+3}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{-3/2}$
$E(\lambda_{ij}^{-1})$ | $\frac{\nu}{2}/(\frac{\nu}{2}-1)$ | $\frac{\nu}{\nu-1}$ | $1$
$E(\lambda_{ij}^{-2})$ | $\left[{\Gamma(\frac{\nu-4}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{2}$ | $\frac{\nu}{\nu-2}$ | $\left[{\Gamma(\frac{\nu+4}{2})}/{\Gamma(\frac{\nu}{2})}\right]\left(\frac{\nu}{2}\right)^{-2}$

##  Appendix B Posterior computations

The prior distributions considered for the parameters, the complete
conditional distributions, and proposal densities considered in the MCMC
algorithm are detailed as follows.

###  B.1 Full conditional of $\mu$

Under the specified Normal prior distribution, the full conditional of $\mu$
is given by

$\displaystyle p(\mu\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{\mu^{2}}{2s_{\mu}^{2}}-\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}\right]\right\\}$

that is,

$\displaystyle\mu\mid\cdot\sim\mathcal{N}\left(m,s^{2}\right)$

where
$m=s^{2}\left[\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}\left(z_{ij}-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)}{\sigma^{2}\left(1-\rho^{2}\right)}\right]$ and
$s^{2}=\left[\frac{1}{\sigma^{2}\left(1-\rho^{2}\right)}\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\lambda_{ij}+\frac{1}{s_{\mu}^{2}}\right]^{-1}$.

###  B.2 Full conditional of $\sigma^{2}$

Under the specified Inverse Gamma prior distribution, the full conditional of
$\sigma^{2}$ is given by

$\displaystyle p(\sigma^{2}\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\left(\frac{1}{\sigma^{2}}\right)^{a_{\sigma}+1+\frac{n(n+1)}{2}}\exp\left\\{-\frac{b_{\sigma}}{\sigma^{2}}-\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}+z_{ij}^{2}\right]\right\\}$

that is,

$\displaystyle\sigma^{2}\mid\cdot\sim\mathcal{IG}\left(a,b\right)$

where $a=a_{\sigma}+\frac{n(n+1)}{2}$ and
$b=b_{\sigma}+\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\lambda_{ij}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{i+j-1}-\rho
Z_{ij}\right)^{2}}{\left(1-\rho^{2}\right)}+T_{ij}^{2}\right]$.

###  B.3 Full conditional of $\rho$

Under the specified beta prior distribution, the full conditional of $\rho$ is
given by

$\displaystyle p(\rho\mid\cdot)$ $\displaystyle\propto$
$\displaystyle(1+\rho)^{c-1-\frac{n(n+1)}{4}}(1-\rho)^{d-1-\frac{n(n+1)}{4}}$
$\displaystyle\times$
$\displaystyle\exp\left\\{-\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}\right]\right\\}\boldsymbol{I}_{(-1,1)}(\rho).$

Note that $p(\rho\mid\cdot)$ does not have a known distribution, so we
obtained samples using the Adaptive Metropolis-Hastings algorithm. Using the
transformation of $\rho$ defined as
$\zeta=\frac{1}{2}\log\left(\frac{1+\rho}{1-\rho}\right)=\mbox{atanh}(\rho)$,
so that $\zeta$ has support in $\mathbb{R}$. Then, the density function of
$\zeta$ is given by

$\displaystyle p\left(\zeta\mid\cdot\right)$ $\displaystyle\propto$
$\displaystyle p(\rho(\zeta)\mid\cdot)\times\mbox{sech}^{2}(\zeta).$

For this, we considered the global adaptive scaling seen in (Andrieu and
Thoms, (2008)), assuming that the actual value of $\zeta$ is $\zeta_{i}$ and
then we will update $\zeta$ so the in following iteration we have
$\zeta_{i+1}$. We follow the steps below:

  1. 1.

Initialise $\zeta_{0}$, $\mu_{0}$, $\kappa_{0}$ and $\sigma_{0}^{2}$.

  2. 2.

At iteration $i+1$, given $\zeta_{i}$, $\mu_{i}$, $\sigma_{i}^{2}$ and
$\kappa_{i}$, sample
$\zeta_{i+1}\sim\mathcal{N}\left(\zeta_{i},\kappa_{i}\sigma_{i}^{2}\right)$
and accept it with probability $\alpha\left(\zeta_{i},\zeta_{i+1}\right)$,
otherwise $\zeta_{i+1}=\zeta_{i}$ and update:

$\displaystyle\gamma_{i+1}$ $\displaystyle=$ $\displaystyle(i+1)^{-a}$
$\displaystyle\log(\kappa_{i+1})$ $\displaystyle=$
$\displaystyle\log(\kappa_{i})+\gamma_{i+1}\left[\alpha\left(\zeta_{i},\zeta_{i+1}\right)-\alpha_{*}\right]$
$\displaystyle\mu_{i+1}$ $\displaystyle=$
$\displaystyle\mu_{i}+\gamma_{i+1}(\zeta_{i+1}-\mu_{i})$
$\displaystyle\sigma_{i+1}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2}+\gamma_{i+1}\left[(\zeta_{i+1}-\mu)^{2}-\sigma_{i}^{2}\right]$

Once $\zeta_{i+1}$ has been obtained, compute $\rho=\mbox{tanh}(\zeta_{i+1})$.
We set $\mu_{0}=0$, $\sigma_{0}^{2}=1$, $\kappa_{0}=2.38$, $a=0.8$ and
$\alpha_{*}=0.234$.

###  B.4 Full conditional of $\sigma_{\alpha}^{2}$

Under the specified Inverse Gamma prior distribution, the full conditional of
$\sigma_{\alpha}^{2}$ is given by

$\displaystyle
p(\sigma_{\alpha}^{2}\mid\cdot)\propto\left(\frac{1}{\sigma_{\alpha}^{2}}\right)^{a_{\alpha}+\frac{n-1}{2}}\exp\left\\{-\frac{b_{\alpha}}{\sigma_{\alpha}^{2}}-\frac{1}{2\sigma_{\alpha}^{2}}\sum_{i=2}^{n}\left(\alpha_{i}-\alpha_{i-1}\right)^{2}\right\\},$

that is,

$\displaystyle\sigma_{\alpha}^{2}\mid\cdot\sim\mathcal{IG}\left(a,b\right),$

where $a=a_{\alpha}+\frac{n-1}{2}$ and
$b=b_{\alpha}+\frac{1}{2}\sum_{i=2}^{n}\left(\alpha_{i}-\alpha_{i-1}\right)^{2}$.

###  B.5 Full conditional of $\sigma_{\beta}^{2}$

Under the specified Inverse Gamma prior distribution, the full conditional of
$\sigma_{\beta}^{2}$ is given by

$\displaystyle
p(\sigma_{\beta}^{2}\mid\dots)\propto\left(\frac{1}{\sigma_{\beta}^{2}}\right)^{a_{\beta}+1+\frac{n(n-1)}{4}}\exp\left\\{-\frac{b_{\beta}}{\sigma_{\beta}^{2}}-\frac{1}{2\sigma_{\beta}^{2}}\sum_{i=2}^{n}\sum_{j=1}^{n-i+1}\left(\beta_{ij}-\beta_{i-1,j}\right)^{2}\right\\},$

that is,

$\displaystyle\sigma_{\beta}^{2}\mid\dots\sim\mathcal{IG}\left(a,b\right),$

where $a=a_{\beta}+\frac{n(n-1)}{4}$ and
$b=b_{\beta}+\frac{1}{2}\sum_{i=2}^{n}\sum_{j=1}^{n-i+1}\left(\beta_{ij}-\beta_{i-1,j}\right)^{2}$.

###  B.6 Full conditional of $\sigma_{\gamma}^{2}$

Under the specified Inverse Gamma prior distribution, the full conditional of
$\sigma_{\gamma}^{2}$ is given by

$\displaystyle
p(\sigma_{\gamma}^{2}\mid\cdot)\propto\left(\frac{1}{\sigma_{\gamma}^{2}}\right)^{a_{\gamma}+1+\frac{n-1}{2}}\exp\left\\{-\frac{b_{\gamma}}{\sigma_{\gamma}^{2}}-\frac{1}{2\sigma_{\gamma}^{2}}\sum_{t=2}^{n}\left(\gamma_{t}-\gamma_{t-1}\right)^{2}\right\\},$

that is,

$\displaystyle\sigma_{\gamma}^{2}\mid\cdot\sim\mathcal{IG}\left(a,b\right),$

where $a=a_{\gamma}+\frac{n-1}{2}$ and
$b=b_{\gamma}+\frac{1}{2}\sum_{t=2}^{n}\left(\gamma_{t}-\gamma_{t-1}\right)^{2}$.

###  B.7 Full conditional of $\lambda_{ij}$ and $\nu$

  * •

Student t case:

As $\lambda_{ij}\sim\mathcal{G}(\nu/2,\nu/2)$, the full conditional of
$\lambda_{ij}$ is given by

$\displaystyle p(\lambda_{ij}\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\lambda_{ij}^{\frac{\nu}{2}}\exp\left\\{-\left(\frac{\nu}{2}+\frac{1}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}+T_{ij}^{2}\right]\right)\lambda_{ij}\right\\},$

that is,

$\displaystyle\lambda_{ij}\mid\cdot\sim\mathcal{G}(a,b),$

where $a=\frac{\nu}{2}+1$ and
$b=\frac{\nu}{2}+\frac{1}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}+T_{ij}^{2}\right]$.

Under the specified Gamma prior distribution, the full conditional of $\nu$ is
given by

$\displaystyle p(\nu\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\nu^{a_{t}-1}\left[\frac{\left(\nu/2\right)^{\frac{\nu}{2}}}{\Gamma(\nu/2)}\right]^{n(n+1)/2}\left(\prod_{i=1}^{n}\prod_{j=1}^{n-i+1}\lambda_{ij}\right)^{\frac{\nu}{2}}\exp\left\\{-\left[b_{t}+\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\lambda_{ij}\right]\nu\right\\},$

Note that $p(\nu\mid\cdot)$ does not have a known distribution, so we obtained
samples using the Adaptive Metropolis-Hastings algorithm. Using the
transformation of $\nu$ defined as $\xi=\log\left(\nu\right)$, so that $\xi$
has support in $\mathbb{R}$. Then the density function of $\zeta$ is given by

$\displaystyle p\left(\xi\mid\cdot\right)$ $\displaystyle\propto$
$\displaystyle p(\nu(\xi)\mid\cdot)\times\exp(\xi)$

Assume that the actual value of $\xi$ is $\xi_{i}$ and then we will update
$\xi$ so the in following iteration we have $\xi_{i+1}$ as:

  1. 1.

Initialise $\xi_{0}$, $\mu_{0}$, $\kappa_{0}$ and $\sigma_{0}^{2}$.

  2. 2.

At iteration $i+1$, given $\xi_{i}$, $\mu_{i}$, $\sigma_{i}^{2}$ and
$\kappa_{i}$, sample
$\xi_{i+1}\sim\mathcal{N}\left(\xi_{i},\kappa_{i}\sigma_{i}^{2}\right)$ and
accept it with probability $\alpha\left(\xi_{i},\xi_{i+1}\right)$, otherwise
$\xi_{i+1}=\xi_{i}$ and update:

$\displaystyle\gamma_{i+1}$ $\displaystyle=$ $\displaystyle(i+1)^{-a}$
$\displaystyle\log(\kappa_{i+1})$ $\displaystyle=$
$\displaystyle\log(\kappa_{i})+\gamma_{i+1}\left[\alpha\left(\xi_{i},\xi_{i+1}\right)-\alpha_{*}\right]$
$\displaystyle\mu_{i+1}$ $\displaystyle=$
$\displaystyle\mu_{i}+\gamma_{i+1}(\xi_{i+1}-\mu_{i})$
$\displaystyle\sigma_{i+1}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2}+\gamma_{i+1}\left[(\xi_{i+1}-\mu)^{2}-\sigma_{i}^{2}\right]$

Once $\xi_{i+1}$ has been obtained, compute $\nu=\exp(\xi_{i+1})$. We set
$\mu_{0}=0$, $\sigma_{0}^{2}=1$, $\kappa_{0}=2.38$, $a=0.8$ and
$\alpha_{*}=0.234$.

  * •

Slash case:

As $\lambda_{ij}\sim\mathcal{B}e(\nu,1)$, the full conditional of
$\lambda_{ij}$ is given by

$\displaystyle p(\lambda_{ij}\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\lambda_{ij}^{\nu}\exp\left\\{-\left(\frac{\nu}{2}+\frac{1}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}+T_{ij}^{2}\right]\right)\lambda_{ij}\right\\}\boldsymbol{I}_{(0,1)}(\lambda_{ij})$

that is,

$\displaystyle\lambda_{ij}\mid\cdot\sim\mathcal{G}_{(0,1)}(a,b)$

where $\mathcal{G}_{(0,1)}(\cdot,\cdot)$ is the right truncated gamma
distribution, $a=\nu+1$ and
$b=\frac{\nu}{2}+\frac{1}{2\sigma^{2}}\left[\frac{\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{1-\rho^{2}}+T_{ij}^{2}\right]$.

Under the specified gamma prior distribution, the full conditional of $\nu$ is
given by

$\displaystyle p(\nu\mid\dots)$ $\displaystyle\propto$
$\displaystyle\nu^{a_{s}+\frac{n(n+1)}{2}-1}\exp\left\\{-\left(b_{s}-\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\log(\lambda_{ij})\right)\nu\right\\}\boldsymbol{I}_{(1,\infty)}(\nu)$

that is,

$\displaystyle\nu\mid\dots\sim\mathcal{G}_{(1,\infty)}(a,b)$

where $\mathcal{G}_{(1,\infty)}(\cdot,\cdot)$ is the right truncated Gamma
distribution, $a=a_{s}+\frac{n(n+1)}{2}$ and
$b=b_{s}-\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\log(\lambda_{ij})$.

  * •

Variance Gamma case:

As $\lambda_{ij}\sim\mathcal{IG}(\nu/2,\nu/2)$, the full conditional of
$\lambda_{ij}$ is given by

$\displaystyle p(\lambda_{ij}\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\lambda_{ij}^{-\frac{\nu}{2}}\exp\left\\{-\frac{1}{2}\left[\lambda_{ij}\left(\frac{T_{ij}^{2}}{\sigma^{2}}+\frac{\left(y_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{\sigma^{2}\left(1-\rho^{2}\right)}\right)+\frac{\nu}{\lambda_{ij}}\right]\right\\}$

that is,

$\displaystyle\lambda_{ij}\mid\cdot\sim\mathcal{GIG}\left(\omega,\chi,\psi\right)$

where $\mathcal{GIG}(\cdot,\cdot,\cdot)$ is the generalized Inverse Gaussian
distribution, $\omega=1-\frac{\nu}{2}$, $\chi=\nu$ and
$\psi=\frac{T_{ij}^{2}}{\sigma^{2}}+\frac{\left(y_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{i+j-1}-\rho
T_{ij}\right)^{2}}{\sigma^{2}\left(1-\rho^{2}\right)}$.

Under the specified gamma prior distribution, the full conditional of $\nu$ is
given by

$\displaystyle p(\nu\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\nu^{a_{vg}-1}\left[\frac{\left(\nu/2\right)^{\frac{\nu}{2}}}{\Gamma(\nu/2)}\right]^{n(n+1)/2}\left(\prod_{i=1}^{n}\prod_{j=1}^{n-i+1}\lambda_{ij}\right)^{-\frac{\nu}{2}}\exp\left\\{-\left[b_{vg}+\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}\frac{1}{\lambda_{ij}}\right]\nu\right\\}$

Note that $p(\nu\mid\dots)$ does not have a known distribution, the Adaptive
Metropolis-Hastings is considered. Under a transformation of $\nu$ defined as
$\xi=\log\left(\nu\right)$, so that $\xi$ has support in $\mathbb{R}$. Then
the density function of $\xi$ is given by

$\displaystyle p\left(\xi|\cdot\right)$ $\displaystyle\propto$ $\displaystyle
p(\nu(\xi)\mid\dots)\cdot\exp(\xi)$

We follow the steps below:

  1. 1.

Initialise $\xi_{0}$, $\mu_{0}$, $\kappa_{0}$ and $\sigma_{0}^{2}$.

  2. 2.

At iteration $i+1$, given $\xi_{i}$, $\mu_{i}$, $\sigma_{i}^{2}$ and
$\kappa_{i}$, sample
$\xi_{i+1}\sim\mathcal{N}\left(\xi_{i},\kappa_{i}\sigma_{i}^{2}\right)$ and
accept it with probability $\alpha\left(\xi_{i},\xi_{i+1}\right)$, otherwise
$\xi_{i+1}=\xi_{i}$ and update:

$\displaystyle\gamma_{i+1}$ $\displaystyle=$ $\displaystyle(i+1)^{-a}$
$\displaystyle\log(\kappa_{i+1})$ $\displaystyle=$
$\displaystyle\log(\kappa_{i})+\gamma_{i+1}\left[\alpha\left(\xi_{i},\xi_{i+1}\right)-\alpha_{*}\right]$
$\displaystyle\mu_{i+1}$ $\displaystyle=$
$\displaystyle\mu_{i}+\gamma_{i+1}(\xi_{i+1}-\mu_{i})$
$\displaystyle\sigma_{i+1}^{2}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2}+\gamma_{i+1}\left[(\xi_{i+1}-\mu)^{2}-\sigma_{i}^{2}\right]$

Once $\xi_{i+1}$ has been obtained, compute $\nu=\exp(\xi_{i+1})$. We set
$\mu_{0}=0$, $\sigma_{0}^{2}=1$, $\kappa_{0}=2.38$, $a=0.8$ and
$\alpha_{*}=0.234$.

###  B.8 Full conditional of $T_{ij}$

The full conditional of $T_{ij}$ is given by

$\displaystyle p\left(T_{ij}\mid\cdot\right)$ $\displaystyle\propto$
$\displaystyle\exp\left\\{-\lambda_{ij}\frac{\left(T_{ij}-\rho(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t})\right)^{2}}{\sigma^{2}(1-\rho^{2})}\right\\}\boldsymbol{I}_{0,\infty}(T_{ij})$

that is,

$\displaystyle
T_{ij}\mid\cdot\sim\mathcal{TN}_{[0,\infty)}\left(m,s^{2}\right)$

where $m=\rho(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t})$ and
$s^{2}=\frac{\sigma^{2}(1-\rho^{2})}{\lambda_{ij}}$.

###  B.9 Full conditional of $\alpha_{i}$

Under the initial condition that $\alpha_{1}=0$, the full conditional of
$\alpha_{i}$ is given by

$\displaystyle p(\alpha_{i}\mid\dots)$ $\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{1}{2}\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{\sigma^{2}\left(1-\rho^{2}\right)}\right\\}$
$\displaystyle\times$
$\displaystyle\exp\left\\{-\frac{1}{2\sigma_{\alpha}^{2}}\left(\alpha_{i}-\alpha_{i-1}\right)^{2}\right\\}\exp\left\\{-\frac{1}{2\sigma_{\alpha}^{2}}\left(\alpha_{i+1}-\alpha_{i}\right)^{2}\boldsymbol{I}_{(i+1<n)}\right\\}$

that is,

$\displaystyle\alpha_{i}\mid\dots\sim\mathcal{N}(m,s^{2})$

where
$m=s^{2}\left[\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}\left(z_{ij}-\mu-\beta_{ij}-\gamma_{i+j-1}-\rho
T_{ij}\right)}{\sigma^{2}\left(1-\rho^{2}\right)}+\frac{\left(\alpha_{i-1}+\alpha_{i+1}\boldsymbol{I}_{(i+1<n)}\right)}{\sigma_{\alpha}^{2}}\right]$
and
$s^{2}=\left[\sum_{j=1}^{n-i+1}\frac{\lambda_{ij}}{\sigma^{2}\left(1-\rho^{2}\right)}+\frac{1+\boldsymbol{I}_{(i+1<n)}}{\sigma_{\alpha}^{2}}\right]^{-1}$.

###  B.10 Full conditional of $\beta_{ij}$

Under the initial conditions that $\beta_{i1}=0$ and $\beta_{1j}=0$ for all i
and j, the full conditional of $\beta_{ij}$ is given by

$\displaystyle p(\beta_{ij}\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{1}{2}\frac{\lambda_{ij}\left(z_{ij}-\mu-\alpha_{i}-\beta_{ij}-\gamma_{t}-\rho
T_{ij}\right)^{2}}{\sigma^{2}\left(1-\rho^{2}\right)}\right\\}$
$\displaystyle\times$
$\displaystyle\exp\left\\{-\frac{1}{2\sigma_{\beta}^{2}}\left(\beta_{ij}-\beta_{i-1,j}\right)^{2}\right\\}\exp\left\\{-\frac{1}{2\sigma_{\beta}^{2}}\left(\beta_{i+1,j}-\beta_{ij}\right)^{2}\boldsymbol{I}_{(i+1<n)}\right\\}$

that is,

$\displaystyle\beta_{ij}\mid\cdot\sim\mathcal{N}(m,s^{2})$

where
$m=s^{2}\left[\frac{\lambda_{ij}\left(z_{ij}-\mu-\alpha_{i}-\gamma_{t}-\rho
T_{ij}\right)}{\sigma^{2}\left(1-\rho^{2}\right)}+\frac{\left(\beta_{i-1,j}+\beta_{i+1,j}\boldsymbol{I}_{(i+1<n)}\right)}{s_{\beta}^{2}}\right]$
and
$s^{2}=\left[\frac{\lambda_{ij}}{\sigma^{2}\left(1-\rho^{2}\right)}+\frac{1+\boldsymbol{I}_{(i+1<n)}}{\sigma_{\beta}^{2}}\right]^{-1}$.

###  B.11 Full conditional of $\gamma_{t}$

Under the initial condition that $\gamma_{1}=0$, the full conditional of
$\gamma_{t}$ is given by

$\displaystyle p(\gamma_{t}\mid\cdot)$ $\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{1}{2}\sum_{j=1}^{t}\frac{\lambda_{t-j+1,j}\left(z_{t-j+1,j}-\mu-\alpha_{t-j+1}-\beta_{t-j+1,j}-\gamma_{t}-\rho
T_{t-j+1,k}\right)^{2}}{\sigma^{2}\left(1-\rho^{2}\right)}\right\\}$
$\displaystyle\times$
$\displaystyle\exp\left\\{-\frac{1}{2\sigma_{\gamma}^{2}}\left(\gamma_{t}-\gamma_{t-1}\right)^{2}\right\\}\exp\left\\{-\frac{1}{2\sigma_{\gamma}^{2}}\left(\gamma_{t+1}-\gamma_{t}\right)^{2}\boldsymbol{I}_{(t+1<n)}\right\\}$

that is,

$\displaystyle\gamma_{t}\mid\cdot\sim\mathcal{N}(m,s^{2})$

where
$m=s^{2}\left[\sum_{j=1}^{t}\frac{\lambda_{t-j+1,j}\left(z_{t-j+1,j}-\mu-\alpha_{t-j+1}-\beta_{t-j+1,j}-\rho
T_{t-j+1,j}\right)}{\sigma^{2}\left(1-\rho^{2}\right)}+\frac{\left(\gamma_{t-1}+\gamma_{t+1}\boldsymbol{I}_{(t+1<n)}\right)}{\sigma_{\gamma}^{2}}\right]$
and
$s^{2}=\left[\sum_{j=1}^{t}\frac{\lambda_{t-j+1,j}}{\sigma^{2}\left(1-\rho^{2}\right)}+\frac{1+\boldsymbol{I}_{(t+1<n)}}{\sigma_{\gamma}^{2}}\right]^{-1}$.

##  Appendix C Model Comparison Criteria

We regard the root mean square prediction error and scoring rules in a
Bayesian context as measures for comparing models based on their posterior
predictive distribution. In particular, we consider the Interval Score, the
Width of Credibility Interval and the Continuous Ranked Probability Score.

###  C.1 Root Mean Square Prediction Error

The root mean squared prediction error measures the expected squared distance
between the predictive value and the true observed value (validation
observation).

$RMSPE(z,z^{pred})=\sqrt{\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\frac{1}{n_{0}}\sum_{j=1}^{n_{0}}\left(z_{ij}-z_{ij}^{pred}\right)^{2}},$
(C.3)

where $z$ is the validation observation, $z^{pred}$ the predictive value and
$n_{0}$ the length of $z$.

###  C.2 Interval Score

Interval forecast is a crucial special case of quantile prediction (Gneiting
and Raftery,, 2007). It compares the predictive credibility interval with the
true observed value (validation observation), and it considers the uncertainty
in the predictions such that the model is penalised if an interval is too
narrow and misses the true value. The Interval Score (IS) is given by

${\rm
IS}(u,l;z)=(u-l)+\frac{2}{\psi}(l-z)I_{\left[z<l\right]}+\frac{2}{\psi}(z-u)I_{\left[z>u\right]},$
(C.4)

where $l$ and $u$ represent for the forecaster quoted $\frac{\psi}{2}$ and
$1-\frac{\psi}{2}$ quantiles based on the predictive distribution and $z$ is
the validation observation. If $\psi=0.05$ the resulting interval has 95%
credibility.

###  C.3 Width of Credibility Interval

The Width of Credibility Interval (WCI) is defined based on credible intervals
and is given by

$WCI(u,l)=(u-l),$ (C.5)

where $l$ and $u$ represent the lower and upper intervals, respectively shown
in the Interval Score measure.

###  C.4 Continuous Ranked Probability Score

This measure generalizes the absolute error which it reduce if $F$ is a
deterministic forecast, that is, a point measure. Thus, the CRPS provides a
direct way to compare deterministic and probabilistic forecasts Gneiting and
Raftery, (2007). The Continuous Ranked Probability Score is given by

$CRPS(F,z)=E_{F}|z^{pred}-{\tilde{z}}^{pred}|-\frac{1}{2}E_{F}|z^{pred}-z|,$
(C.6)

where $z^{pred}$ and $\tilde{z}^{pred}$ are independent replicas of the
predictive distribution, $z$ the validation observation and $F$ the predictive
distribution of the model.

##  Appendix D Run-off triangle data

Table D.2: Dataset1: The amount of claims paid to the insureds of an insurance
product during the period from 1978 to 1995.

accident development year $j$ year $i$ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
17 18 1978 3323 8332 9572 10172 7631 3855 3252 4433 2188 333 199 692 311 0 405
293 76 14 1979 3785 10342 8330 7849 2839 3577 1404 1721 1065 156 35 259 250
420 6 1 0 1980 4677 9989 8745 10228 8572 5787 3855 1445 1612 626 1172 589 438
473 370 31 1981 5288 8089 12839 11828 7560 6383 4118 3016 1575 1985 2645 266
38 45 115 1982 2294 9869 10242 13808 8775 5419 2424 1597 4149 1296 917 295 428
359 1983 3600 7514 8247 9327 8584 4245 4096 3216 2014 593 1188 691 368 1984
3642 7394 9838 9733 6377 4884 11920 4188 4492 1760 944 921 1985 2463 5033 6980
7722 6702 7834 5579 3622 1300 3069 1370 1986 2267 5959 6175 7051 8102 6339
6978 4396 3107 903 1987 2009 3700 5298 6885 6477 7570 5855 5751 3871 1988 1860
5282 3640 7538 5157 5766 6862 2572 1989 2331 3517 5310 6066 10149 9265 5262
1990 2314 4487 4112 7000 11163 10057 1-19 1991 2607 3952 8228 7895 9317 1992
2595 5403 6579 15546 1993 3155 4974 7961 1994 2626 5704 1995 2827

## References

  * Andrieu and Thoms, (2008) Andrieu, C. and Thoms, J. (2008).  A tutorial on adaptive MCMC.  Statistics and Computing, 18, issue 4:343 – 373.  Publisher: Springer Netherlands.
  * Arellano-Valle et al., (2007) Arellano-Valle, R., Bolfarine, H., and Lachos, V. (2007).  Bayesian inference for skew-normal linear mixed models.  Journal of Applied Statistics, 34(6):663–682.
  * Arellano-Valle et al., (2005) Arellano-Valle, R. B., Bolfarine, H., and Lachos, V. H. (2005).  Skew-normal linear mixed models.  Journal of Data Science, 3(4):415–438.
  * Azzalini, (1985) Azzalini, A. (1985).  A class of distributions which includes the normal ones.  Scandinavian Journal of Statistics, 12(2):171–178.
  * Azzalini and Valle, (1996) Azzalini, A. and Valle, A. D. (1996).  The multivariate skew-normal distribution.  Biometrika, 83(4):715–726.
  * Bornhuetter and Ferguson, (1972) Bornhuetter, R. L. and Ferguson, R. E. (1972).  The actuary and IBNR.  (59).
  * Chan et al., (2008) Chan, J. S., Boris Choy, S., and Makov, U. E. (2008).  Robust bayesian analysis of loss reserves data using the generalized-t distribution.  ASTIN Bulletin, 38(1):207–230.
  * Choy et al., (2016) Choy, S. B., Chan, J. S., and Makov, U. E. (2016).  Robust bayesian analysis of loss reserving data using scale mixtures distributions.  Journal of Applied Statistics, 43(3):396–411.
  * de Jong and Zehnwirth, (1983) de Jong, P. and Zehnwirth, B. (1983).  Claims reserving, state-space models and the kalman filter.  Journal of the Institute of Actuaries, 110(1):157–181.
  * Eddelbuettel, (2013) Eddelbuettel, D. (2013).  Seamless R and C++ Integration with Rcpp.  Springer, New York.  ISBN 978-1-4614-6867-7.
  * Eddelbuettel and Balamuta, (2018) Eddelbuettel, D. and Balamuta, J. J. (2018).  Extending extitR with extitC++: A Brief Introduction to extitRcpp.  The American Statistician, 72(1):28–36.
  * England and Verrall, (2002) England, P. and Verrall, R. (2002).  Stochastic claims reserving in general insurance.  British Actuarial Journal, 8(3):443–518.
  * Fernández and Steel, (1998) Fernández, C. and Steel, M. F. J. (1998).  On bayesian modeling of fat tails and skewness.  Journal of the American Statistical Association, 93(441):359–371.
  * Fonseca et al., (2008) Fonseca, T. C. O., Ferreira, M. A. R., and Migon, H. S. (2008).  Objective Bayesian analysis for the Student-t regression model.  Biometrika, 95(2):325–333.
  * Gamerman and Lopes, (2006) Gamerman, D. and Lopes, H. (2006).  Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference.  Texts in Statistical Science. Taylor & Francis.
  * Gelman, (2006) Gelman, A. (2006).  Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).  Bayesian Analysis, 1(3):515 – 534.
  * Gneiting and Raftery, (2007) Gneiting, T. and Raftery, A. E. (2007).  Strictly proper scoring rules, prediction and estimation.  Journal of the American Statistical Association, 102(477):360–378.
  * Goudarzi and Zokaei, (2021) Goudarzi, M. and Zokaei, M. (2021).  Bayesian modeling of multivariate loss reserving data based on scale mixtures of multivariate normal distributions: estimation and case influence diagnostics.  Communications in Statistics - Theory and Methods, 50(21):4934–4962.
  * Haberman and Renshaw, (1996) Haberman, S. and Renshaw, A. E. (1996).  Generalized linear models and actuarial science.  Journal of the Royal Statistical Society: Series D (The Statistician), 45(4):407–436.
  * Henze, (1986) Henze, N. (1986).  A probabilistic representation of the ’skew-normal’ distribution.  Scandinavian Journal of Statistics, 13:271–275.
  * Kim, (2001) Kim, H.-J. (2001).  On a skew-t distribution.  Communications for Statistical Applications and Methods, 8:867–873.
  * Kremer, (1982) Kremer, E. (1982).  Ibnr-claims and the two-way model of anova.  Scandinavian Actuarial Journal, 1982:47–55.
  * Lunn et al., (2009) Lunn, D., Spiegelhalter, D., Thomas, A., and Best, N. (2009).  The BUGS project: Evolution, critique and future directions (with discussion).  Stat Med, 28:3049–3082.
  * Lunn et al., (2000) Lunn, D. J., Thomas, A., Best, N., and Spiegelhalter, D. (2000).  Winbugs a bayesian modelling framework: Concepts, structure, and extensibility.  Statistics and Computing, 10(4):325–337.
  * Mack, (1993) Mack, T. (1993).  Distribution-free calculation of the standard error of chain ladder reserve estimates.  ASTIN Bulletin, 23(2):213–225.
  * Madan and Seneta, (1990) Madan, D. B. and Seneta, E. (1990).  The variance gamma (v.g.) model for share market returns.  The Journal of Business, 63(4):511–524.
  * Ntzoufras and Dellaportas, (2002) Ntzoufras, I. and Dellaportas, P. (2002).  Bayesian modelling of outstanding liabilities incorporating claim count uncertainty.  North American Actuarial Journal, 6(1):113–125.
  * Philippe, (1997) Philippe, A. (1997).  Simulation of right- and left-truncated gamma distribution by mixtures.  Statistics and Computing, 7:173–181.
  * Pigeon et al., (2013) Pigeon, M., Antonio, K., and Denuit, M. (2013).  Individual loss reserving with the multivariate skew normal framework.  ASTIN Bulletin, 43(3):399–428.
  * Pérez-Rodríguez et al., (2018) Pérez-Rodríguez, P., Acosta-Pech, R., Pérez-Elizalde, S., Cruz, C. V., Espinosa, J. S., and Crossa, J. (2018).  A Bayesian Genomic Regression Model with Skew Normal Random Errors.  G3 Genes—Genomes—Genetics, 8(5):1771–1785.
  * Renshaw and Verrall, (1998) Renshaw, A. and Verrall, R. (1998).  A stochastic model underlying the chain-ladder technique.  British Actuarial Journal, 4(4):903–923.
  * Sanderson and Curtin, (2016) Sanderson, C. and Curtin, R. (2016).  Armadillo: a template-based c++ library for linear algebra.  Journal of Open Source Software, 1(2):26.
  * Sanderson and Curtin, (2018) Sanderson, C. and Curtin, R. (2018).  A user-friendly hybrid sparse matrix class in c$+$$+$.  In Mathematical Software – ICMS 2018, pages 422–430. Springer International Publishing.
  * Shi et al., (2012) Shi, P. A., Basu, S., and Glenn G. Meyers, F. (2012).  A bayesian log-normal model for multivariate loss reserving.  North American Actuarial Journal, 16(1):29–51.
  * Stan Development Team, (2023) Stan Development Team (2023).  RStan: the R interface to Stan.  R package version 2.21.8.
  * Taylor, (1984) Taylor, G. (1984).  Claim reserving in non-life insurance.  North-Holland.
  * Team, (2021) Team, R. C. (2021).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Usman and Chan, (2022) Usman, F. and Chan, J. S. (2022).  New loss reserve models with persistence effects to forecast trapezoidal losses in run-off triangles.  ASTIN Bulletin: The Journal of the IAA, 52(3):877–920.
  * Verrall, (1990) Verrall, R. (1990).  Bayes and empirical bayes estimation for the chain ladder model.  ASTIN Bulletin, 20(2):217–243.
  * Verrall, (1996) Verrall, R. (1996).  Claims reserving and generalised additive models.  Insurance: Mathematics and Economics, 19(1):31–43.
  * Verrall, (2000) Verrall, R. (2000).  An investigation into stochastic claims reserving models and the chain-ladder technique.  Insurance: Mathematics and Economics, 26(1):91–99.
  * Verrall, (1989) Verrall, R. J. (1989).  A state space representation of the chain ladder linear model.  Journal of the Institute of Actuaries, 116(3):589–609.
  * Vieira, (2019) Vieira, F. V. G. (2019).  Bayesian quantile analysis for loss reserving models.  Master’s thesis, Universidade Federal do Rio de Janeiro - UFRJ.
  * Vilela, (2013) Vilela, J. M. (2013).  Modelos para provisão IBNR: uma abordagem bayesiana.  Master’s thesis, Universidade Federal do Rio de Janeiro - UFRJ.
  * Wang and Genton, (2006) Wang, J. and Genton, M. G. (2006).  The multivariate skew-slash distribution.  Journal of Statistical Planning and Inference, 136(1):209–220.