\documentclass[x11names]{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
  \usepackage{amsmath}
  \usepackage[round]{natbib}
  \usepackage[english]{babel}
  \usepackage[scaled=0.9]{helvet}
  \usepackage[sc]{mathpazo}
  \usepackage{framed,booktabs}
  \usepackage[noae]{Sweave}

  %\VignetteIndexEntry{Credibility theory}
  %\VignettePackage{actuar}

  \title{Credibility theory features of \pkg{actuar}}
  \author{Christophe Dutang \\ Université du Mans \\[3ex]
    Vincent Goulet \\ Université Laval \\[3ex]
    Xavier Milhaud \\ Université Claude Bernard Lyon 1 \\[3ex]
    Mathieu Pigeon \\ Université du Québec à Montréal}
  \date{}

  %% Colors
  \usepackage{xcolor}
  \definecolor{link}{rgb}{0,0.4,0.6}             % internal links
  \definecolor{url}{rgb}{0.6,0,0}                % external links
  \definecolor{citation}{rgb}{0,0.5,0}           % citations
  \definecolor{codebg}{named}{LightYellow1}      % R code background

  %% Hyperlinks
  \usepackage{hyperref}
  \hypersetup{%
    pdfauthor={Vincent Goulet},
    colorlinks = {true},
    linktocpage = {true},
    urlcolor = {url},
    linkcolor = {link},
    citecolor = {citation},
    pdfpagemode = {UseOutlines},
    pdfstartview = {Fit},
    bookmarksopen = {true},
    bookmarksnumbered = {true},
    bookmarksdepth = {subsubsection}}

  %% Sweave environments redefined to have a color background (using
  %% environment snugshade* of package framed).
  \DefineVerbatimEnvironment{Sinput}{Verbatim}{}
  \DefineVerbatimEnvironment{Soutput}{Verbatim}{}
  \fvset{listparameters={\setlength{\topsep}{0pt}}}
  \renewenvironment{Schunk}{%
    \setlength{\topsep}{0pt}
    \colorlet{shadecolor}{codebg}
    \begin{snugshade*}}%
    {\end{snugshade*}}

  %% Some new commands. Commands \E and \VAR redefined to fit within
  %% the author's writing habits...
  \newcommand{\E}[1]{\mathrm{E}[ #1 ]}
  \newcommand{\proglang}[1]{\textsf{#1}}
  \newcommand{\pkg}[1]{\textbf{#1}}
  \newcommand{\code}[1]{\texttt{#1}}
  \newcommand{\pt}{{\scriptscriptstyle \Sigma}}

  \bibliographystyle{plainnat}

<<echo=FALSE>>=
library(actuar)
options(width = 60, digits = 4)
@

\begin{document}

\maketitle

\section{Introduction}
\label{sec:introduction}

Credibility models are actuarial tools to distribute premiums fairly
among a heterogeneous group of policyholders (henceforth called
\emph{entities}). More generally, they can be seen as prediction methods
applicable in any setting where repeated measures are made for
subjects with different risk levels.

The credibility theory features of \pkg{actuar} consist of matrix
\code{hachemeister} containing the famous data set of
\cite{Hachemeister_75} and function \code{cm} to fit hierarchical
(including Bühlmann, Bühlmann-Straub) and regression credibility
models. Furthermore, function \code{simul} can simulate portfolios of
data satisfying the assumptions of the aforementioned credibility
models; see the \code{"simulation"} vignette for details.


\section{Hachemeister data set}
\label{sec:hachemeister}

The data set of \cite{Hachemeister_75} consists of private passenger
bodily injury insurance average claim amounts, and the corresponding
number of claims, for five U.S.\ states over 12 quarters between July
1970 and June 1973. The data set is included in the package in the
form of a matrix with 5 rows and 25 columns. The first column contains
a state index, columns 2--13 contain the claim averages and columns
14--25 contain the claim numbers:
<<echo=TRUE>>=
data(hachemeister)
hachemeister
@

\section{Hierarchical credibility model}
\label{sec:hierarchical}

The linear model fitting function of \proglang{R} is named \code{lm}.
Since credibility models are very close in many respects to linear
models, and since the credibility model fitting function of
\pkg{actuar} borrows much of its interface from \code{lm}, we named
the credibility function \code{cm}.

Function \code{cm} acts as a unified interface for all credibility
models supported by the package. Currently, these are the
unidimensional models of \cite{Buhlmann_69} and \cite{BS_70}, the
hierarchical model of \cite{Jewell_75} (of which the first two are
special cases) and the regression model of \cite{Hachemeister_75},
optionally with the intercept at the barycenter of time
\citep[Section~8.4]{Buhlmann_Gisler}. The modular design of \code{cm}
makes it easy to add new models if desired.

This subsection concentrates on usage of \code{cm} for hierarchical
models.

There are some variations in the formulas of the hierarchical model in
the literature. We compute the credibility premiums as given in
\cite{BJ_87} or \cite{Buhlmann_Gisler}. We support three types of
estimators of the between variance structure parameters: the unbiased
estimators of \cite{Buhlmann_Gisler} (the default), the slightly
different version of \cite{Ohlsson} and the iterative
pseudo-estimators as found in \cite{LivreVert} or \cite{Goulet_JAP}.

Consider an insurance portfolio where contracts are classified into
cohorts. In our terminology, this is a two-level hierarchical
classification structure. The observations are claim amounts
$S_{ijt}$, where index $i = 1, \dots, I$ identifies the cohort, index
$j = 1, \dots, J_i$ identifies the contract within the cohort and
index $t = 1, \dots, n_{ij}$ identifies the period (usually a year).
To each data point corresponds a weight --- or volume -- $w_{ijt}$.
Then, the best linear prediction for the next period outcome of a
contract based on ratios $X_{ijt} = S_{ijt}/w_{ijt}$ is
\begin{equation}
  \label{eq:hierarchical:premiums}
  \begin{split}
    \hat{\pi}_{ij} &= z_{ij} X_{ijw} + (1 - z_{ij}) \hat{\pi}_i \\
    \hat{\pi}_i    &= z_i X_{izw} + (1 - z_i) m \\
  \end{split}
\end{equation}
with the credibility factors
\begin{align*}
  z_{ij}
  &= \frac{w_{ij\pt}}{w_{ijk\pt} + s^2/a}, &
  w_{ij\pt}
  &= \sum_{t = 1}^{n_{ij}} w_{ijt} \\
  z_{i}
  &= \frac{z_{i\pt}}{z_{i\pt} + a/b}, &
  z_{i\pt}
  &= \sum_{j = 1}^{J_i} z_{ij}
\end{align*}
and the weighted averages
\begin{align*}
  X_{ijw} &= \sum_{t = 1}^{n_{ij}} \frac{w_{ijt}}{w_{ij\pt}}\, X_{ijt} \\
  X_{izw} &= \sum_{j = 1}^{J_i} \frac{z_{ij}}{z_{i\pt}}\, X_{ijw}.
\end{align*}

The estimator of $s^2$ is
\begin{equation}
  \label{eq:s2}
  \hat{s}^2 = \frac{1}{\sum_{i = 1}^I \sum_{j = 1}^{J_i} (n_{ij} - 1)}
  \sum_{i = 1}^I \sum_{j = 1}^{J_i} \sum_{t = 1}^{n_{ij}}
  w_{ijt} (X_{ijt} - X_{ijw})^2.
\end{equation}

The three types of estimators for parameters $a$ and $b$ are the
following. First, let
\begin{align*}
  A_i
  &= \sum_{j = 1}^{J_i} w_{ij\pt} (X_{ijw} - X_{iww})^2 -
  (J_i - 1) s^2 &
  c_i
  &= w_{i\pt\pt} - \sum_{j = 1}^{J_i} \frac{w_{ij\pt}^2}{w_{i\pt\pt}} \\
  B
  &= \sum_{i = 1}^I z_{i\pt} (X_{izw} - \bar{X}_{zzw})^2 - (I - 1) a &
  d
  &= z_{\pt\pt} - \sum_{i = 1}^I \frac{z_{i\pt}^2}{z_{\pt\pt}},
\end{align*}
with
\begin{equation}
  \label{eq:Xbzzw}
  \bar{X}_{zzw} = \sum_{i = 1}^I \frac{z_{i\pt}}{z_{\pt\pt}}\, X_{izw}.
\end{equation}
(Hence, $\E{A_i} = c_i a$ and $\E{B} = d b$.) Then, the
Bühlmann--Gisler estimators are
\begin{align}
  \label{eq:ac-BG}
  \hat{a}
  &= \frac{1}{I} \sum_{i = 1}^I \max \left( \frac{A_i}{c_i}, 0 \right) \\
  \label{eq:bc-BG}
  \hat{b}
  &= \max \left( \frac{B}{d}, 0 \right),
\end{align}
the Ohlsson estimators are
\begin{align}
  \label{eq:ac-Ohl}
  \hat{a}^\prime
  &= \frac{\sum_{i = 1}^I A_i}{\sum_{i = 1}^I c_i} \\
  \label{eq:bc-Ohl}
  \hat{b}^\prime
  &= \frac{B}{d}
\end{align}
and the iterative (pseudo-)estimators are
\begin{align}
  \label{eq:at}
  \tilde{a}
  &= \frac{1}{\sum_{i = 1}^I (J_i - 1)}
  \sum_{i = 1}^I \sum_{j = 1}^{J_i} z_{ij} (X_{ijw} - X_{izw})^2 \\
  \label{eq:bt}
  \tilde{b}
  &= \frac{1}{I - 1}
  \sum_{i = 1}^I z_i (X_{izw} - X_{zzw})^2,
\end{align}
where
\begin{equation}
  \label{eq:Xzzw}
  X_{zzw} = \sum_{i = 1}^I \frac{z_i}{z_\pt}\, X_{izw}.
\end{equation}
Note the difference between the two weighted averages \eqref{eq:Xbzzw}
and \eqref{eq:Xzzw}. See \cite{cm} for further discussion on this
topic.

Finally, the estimator of the collective mean $m$ is $\hat{m} =
X_{zzw}$.

The credibility modeling function \code{cm} assumes that data is available in
the format most practical applications would use, namely a rectangular
array (matrix or data frame) with entity observations in the rows and
with one or more classification index columns (numeric or character).
One will recognize the output format of \code{simul} and its summary
methods.

Then, function \code{cm} works much the same as \code{lm}. It takes in
argument: a formula of the form \code{\~{} terms} describing the
hierarchical interactions in a data set; the data set containing the
variables referenced in the formula; the names of the columns
where the ratios and the weights are to be found in the data set. The
latter should contain at least two nodes in each level and more than
one period of experience for at least one entity. Missing values are
represented by \code{NA}s. There can be entities with no experience
(complete lines of \code{NA}s).

In order to give an easily reproducible example, we group states 1 and
3 of the Hachemeister data set into one cohort and states 2, 4 and 5
into another. This shows that data does not have to be sorted by
level. The fitted model using the iterative estimators is:
<<echo=TRUE>>=
X <- cbind(cohort = c(1, 2, 1, 2, 2), hachemeister)
fit <- cm(~cohort + cohort:state, data = X,
          ratios = ratio.1:ratio.12,
          weights = weight.1:weight.12,
          method = "iterative")
fit
@

The function returns a fitted model object of class \code{"cm"}
containing the estimators of the structure parameters. To compute the
credibility premiums, one calls a method of \code{predict} for this
class:
<<echo=TRUE>>=
predict(fit)
@

One can also obtain a nicely formatted view of the most important
results with a call to \code{summary}:
<<echo=TRUE>>=
summary(fit)
@

The methods of \code{predict} and \code{summary} can both report for a
subset of the levels by means of an argument \code{levels}. For
example:
<<echo=TRUE>>=
summary(fit, levels = "cohort")
predict(fit, levels = "cohort")
@

The results above differ from those of \cite{LivreVert} for the same
example because the formulas for the credibility premiums are
different.


\section{Bühlmann and Bühlmann--Straub models}
\label{sec:buhlmann}

As mentioned above, the Bühlmann and Bühlmann--Straub models are
simply one-level hierarchical models. In this case, the
Bühlmann--Gisler and Ohlsson estimators of the between variance
parameters are both identical to the usual \cite{BS_70} estimator
\begin{equation}
  \label{eq:a-hat}
  \hat{a} = \frac{w_{\pt\pt}}{w_{\pt\pt}^2 - \sum_{i=1}^I
    w_{i\pt}^2}
  \left(
    \sum_{i=1}^I w_{i\pt} (X_{iw} - X_{ww})^2 - (I - 1) \hat{s}^2
  \right),
\end{equation}
and the iterative estimator
\begin{equation}
  \label{eq:a-tilde}
  \tilde{a} = \frac{1}{I - 1} \sum_{i = 1}^I z_i (X_{iw} - X_{zw})^2
\end{equation}
is better known as the Bichsel--Straub estimator.

To fit the Bühlmann model using \code{cm}, one simply does not
specify any weights:
<<echo=TRUE>>=
cm(~state, hachemeister, ratios = ratio.1:ratio.12)
@

In comparison, the results for the Bühlmann--Straub model using the
Bichsel--Straub estimator are:
<<echo=TRUE>>=
cm(~state, hachemeister, ratios = ratio.1:ratio.12,
   weights = weight.1:weight.12)
@


\section{Regression model of Hachemeister}
\label{sec:regression}

The regression model of \cite{Hachemeister_75} is a generalization of
the Bühlmann--Straub model. If data shows a systematic trend, the
latter model will typically under- or over-estimate the true premium
of an entity. The idea of \citeauthor{Hachemeister_75} was to fit to
the data a regression model where the parameters are a credibility
weighted average of an entity's regression parameters and the group's
parameters.

In order to use \code{cm} to fit a credibility regression model to a
data set, one simply has to supply as additional arguments
\code{regformula} and \code{regdata}. The first one is a formula of
the form \code{\~{} terms} describing the regression model and the
second is a data frame of regressors. That is, arguments
\code{regformula} and \code{regdata} are in every respect equivalent
to arguments \code{formula} and \code{data} of \code{lm}, with the
minor difference that \code{regformula} does not need to have a left
hand side (and is ignored if present). For example, fitting the model
\begin{displaymath}
  X_{it} = \beta_0 + \beta_1 t + \varepsilon_t, \quad
  t = 1, \dots, 12
\end{displaymath}
to the original data set of \cite{Hachemeister_75} is done with
<<echo=TRUE>>=
fit <- cm(~state, hachemeister, regformula = ~ time,
          regdata = data.frame(time = 1:12),
          ratios = ratio.1:ratio.12,
          weights = weight.1:weight.12)
fit
@

Computing the credibility premiums requires to give the ``future''
values of the regressors as in \code{predict.lm}:
<<echo=TRUE>>=
predict(fit, newdata = data.frame(time = 13))
@

It is well known that the basic regression model has a major drawback:
there is no guarantee that the credibility regression line will lie
between the collective and individual ones. This may lead to grossly
inadequate premiums, as Figure~\ref{fig:state4} shows.

\begin{figure}[t]
  \centering
<<echo=FALSE, fig=TRUE>>=
plot(NA, xlim = c(1, 13), ylim = c(1000, 2000), xlab = "", ylab = "")
x <- cbind(1, 1:12)
lines(1:12, x %*% fit$means$portfolio,
      col = "blue", lwd = 2)
lines(1:12, x %*% fit$means$state[, 4],
      col = "red", lwd = 2, lty = 2)
lines(1:12, x %*% coefficients(fit$adj.models[[4]]),
      col = "darkgreen", lwd = 2, lty = 3)
points(13, predict(fit, newdata = data.frame(time = 13))[4],
       pch = 8, col = "darkgreen")
legend("bottomright",
       legend = c("collective", "individual", "credibility"),
       col = c("blue", "red", "darkgreen"), lty = 1:3)
@
  \caption{Collective, individual and credibility regression lines for
    State 4 of the Hachemeister data set. The point indicates the
    credibility premium.}
  \label{fig:state4}
\end{figure}

The solution proposed by \cite{Buhlmann:regression:1997} is simply to
position the intercept at the barycenter of time instead of at time
origin \citep[see also][Section~8.4]{Buhlmann_Gisler}. In mathematical
terms, this essentially amounts to using an orthogonal design matrix.
By setting the argument \code{adj.intercept} to \code{TRUE} in the
call, \code{cm} will automatically fit the credibility regression
model with the intercept at the barycenter of time. The resulting
regression coefficients have little meaning, but the predictions are
sensible:
<<echo=TRUE>>=
fit2 <- cm(~state, hachemeister, regformula = ~ time,
           regdata = data.frame(time = 1:12),
           adj.intercept = TRUE,
           ratios = ratio.1:ratio.12,
           weights = weight.1:weight.12)
summary(fit2, newdata = data.frame(time = 13))
@ %
Figure~\ref{fig:state4:2} shows the beneficient effect of the
intercept adjustment on the premium of State 4.

\begin{figure}[t]
  \centering
<<echo=FALSE, fig=TRUE>>=
plot(NA, xlim = c(1, 13), ylim = c(1000, 2000), xlab = "", ylab = "")
x <- cbind(1, 1:12)
R <- fit2$transition
lines(1:12, x %*% solve(R, fit2$means$portfolio),
      col = "blue", lwd = 2)
lines(1:12, x %*% solve(R, fit2$means$state[, 4]),
      col = "red", lwd = 2, lty = 2)
lines(1:12, x %*% solve(R, coefficients(fit2$adj.models[[4]])),
      col = "darkgreen", lwd = 2, lty = 3)
points(13, predict(fit2, newdata = data.frame(time = 13))[4],
       pch = 8, col = "darkgreen")
legend("bottomright",
       legend = c("collective", "individual", "credibility"),
       col = c("blue", "red", "darkgreen"), lty = 1:3)
@
  \caption{Collective, individual and credibility regression lines for
    State 4 of the Hachemeister data set when the intercept is
    positioned at the barycenter of time. The point indicates the
    credibility premium.}
  \label{fig:state4:2}
\end{figure}

\bibliography{actuar}

\end{document}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% coding: utf-8
%%% End:
