\documentclass{article}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
  \usepackage{amsmath,amsthm,color,hyperref}
  \usepackage[round]{natbib}
  \usepackage[english]{babel}
  \usepackage[scaled=0.9]{helvet}
  \usepackage[sc]{mathpazo}
  \usepackage{paralist}
  \usepackage[noae]{Sweave}

  %\VignetteIndexEntry{Simulation of compound hierarchical models}
  %\VignettePackage{actuar}

  \title{Simulation of compound hierarchical models}
  \author{Christophe Dutang \\ ISFA, Université Claude Bernard Lyon 1 \\[3ex]
    Vincent Goulet \\ École d'actuariat, Université Laval \\[3ex]
    Mathieu Pigeon \\ École d'actuariat, Université Laval \\[3ex]
    Louis-Philippe Pouliot \\ École d'actuariat, Université Laval}
  \date{}

  %% Some new commands
  \newcommand{\proglang}[1]{\textsf{#1}}
  \newcommand{\pkg}[1]{\textbf{#1}}
  \newcommand{\code}[1]{\texttt{#1}}

  %% Redefine Sweave environments with smaller font
  \RecustomVerbatimEnvironment{Sinput}{Verbatim}{%
    fontshape=sl,fontsize=\small,xleftmargin=0pt}
  \RecustomVerbatimEnvironment{Soutput}{Verbatim}{%
    fontsize=\small,xleftmargin=0pt}

  \bibliographystyle{plainnat}

  \definecolor{Red}{rgb}{0.7,0,0}
  \definecolor{Blue}{rgb}{0,0,0.8}
  \hypersetup{%
    hyperindex = {true},
    colorlinks = {true},
    linktocpage = {true},
    plainpages = {false},
    linkcolor = {Blue},
    citecolor = {Blue},
    urlcolor = {Red},
    pdfstartview = {Fit},
    pdfpagemode = {UseOutlines},
    pdfview = {XYZ null null null}
  }

  \theoremstyle{definition}
  \newtheorem{example}{Example}

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

\begin{document}

\maketitle

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

Hierarchical probability models are widely used for data classified in
a tree-like structure and in Bayesian inference. The main
characteristic of such models is to have the probability law at some
level in the classification structure be conditional on the outcome in
previous levels. For example, adopting a bottom to top description of
the model, a simple hierarchical model could be written as
\begin{equation}
  \label{eq:basic_model}
  \begin{split}
    X_t|\Lambda, \Theta
    &\sim \text{Poisson}(\Lambda) \\
    \Lambda|\Theta
    &\sim \text{Gamma}(3, \Theta) \\
    \Theta &\sim \text{Gamma}(2, 2),
  \end{split}
\end{equation}
where $X_t$ represents actual data. The random variables $\Theta$ and
$\Lambda$ are generally seen as uncertainty, or risk, parameters in
the actuarial literature; in the sequel, we refer to them as
mixing parameters.

The example above is merely a multi-level mixture of models, something
that is simple to simulate ``by hand''. The following \proglang{R}
expression will yield $n$ variates of the random variable $X_t$:
<<echo=TRUE, eval=FALSE>>=
rpois(n, rgamma(n, 3, rgamma(n, 2, 2)))
@

However, for categorical data common in actuarial applications there
will usually be many categories --- or \emph{nodes} --- at each level.
Simulation is then complicated by the need to always use the correct
parameters for each variate. Furthermore, actuaries often need to
simulate both the frequency and the severity of claims for compound
models of the form
\begin{equation}
  \label{eq:definition-S}
  S = C_1 + \dots + C_N,
\end{equation}
where $C_1, C_2, \dots$ are mutually independent and identically
distributed random variables each independent of $N$.

The package provides function \code{simul} to simulate data from
compound models like \eqref{eq:definition-S} where both the frequency
and the severity components can have a hierarchical structure. The
function also supports weights (or volumes) in the model. We give here
a brief description of the function and its usage; see
\cite{Goulet:simpf:2008} for details about the supported models and
more thorough examples.


\section{Description of hierarchical models}
\label{sec:description}

We need a method to describe hierarchical models in \proglang{R} that
will meet the following criteria:
\begin{enumerate}
\item simple and intuitive to go from the mathematical formulation of
  the model to the \proglang{R} formulation and back;
\item allows for any number of levels and nodes;
\item at any level, allows for any use of parameters higher in the
  hierarchical structure.
\end{enumerate}

A hierarchical model is completely specified by the number of nodes at
each level ($I$, $J_1, \dots, J_I$ and $n_{11}, \dots, n_{IJ}$, above)
and by the probability laws at each level. The number of nodes is
passed to \code{simul} by means of a named list where each element is
a vector of the number of nodes at a given level. Vectors are recycled
when the number of nodes is the same throughout a level. Probability
models are expressed in a semi-symbolic fashion using an object of
mode \code{"expression"}. Each element of the object must be named ---
with names matching those of the number of nodes list --- and should
be a complete call to an existing random number generation function,
but with the number of variates omitted. Hierarchical models are
achieved by replacing one or more parameters of a distribution at a
given level by any combination of the names of the levels above. If no
mixing is to take place at a level, the model for this level can be
\code{NULL}.

\begin{example}
  Consider the following expanded version of model \eqref{eq:basic_model}:
  \begin{align*}
    X_{ijt}|\Lambda_{ij}, \Theta_i
    &\sim \text{Poisson}(\Lambda_{ij}), & t &= 1, \dots, n_{ij} \\
    \Lambda_{ij}|\Theta_i
    &\sim \text{Gamma}(3, \Theta_i), & j &= 1, \dots, J_i \\
    \Theta_i
    &\sim \text{Gamma}(2, 2), & i &= 1, \dots, I,
  \end{align*}
  with
  $I = 3$,
  $J_1 = 4$,
  $J_2 = 5$,
  $J_3 = 6$ and
  $n_{ij} \equiv n = 10$.
  Then the number of nodes is specified by
\begin{verbatim}
list(Theta = 3, Lambda = c(4, 5, 6), Data = 10)
\end{verbatim}
  and the probability model is expressed as
\begin{verbatim}
expression(Theta = rgamma(2, 2),
           Lambda = rgamma(3, Theta),
           Data = rpois(Lambda))
\end{verbatim}
  \qed
\end{example}

Storing the probability model requires an expression object in order
to avoid evaluation of the incomplete calls to the random number
generation functions. Function \code{simul} builds and executes the
calls to the random generation functions from the top of the
hierarchical model to the bottom. At each level, the function
\begin{inparaenum}[1)]
\item infers the number of variates to generate from the number of
  nodes list, and
\item appropriately recycles the mixing parameters simulated
  previously.
\end{inparaenum}

The actual names in the list and the expression object can be
anything; they merely serve to identify the mixing parameters.
Furthermore, any random generation function can be used. The only
constraint is that the name of the number of variates argument is
\code{n}.

Function \code{simul} also supports usage of weights in models. These
usually modify the frequency parameters to take into account the
``size'' of an entity. Weights are used in simulation wherever the
name \code{weights} appears in a model.


\section{Usage}
\label{sec:usage}

Function \code{simul} can simulate data for structures where both the
frequency model and the severity model are hierarchical. It has four
main arguments:
\begin{inparaenum}[1)]
\item \code{nodes} for the number of nodes list;
\item \code{model.freq} for the frequency model;
\item \code{model.sev} for the severity model;
\item \code{weights} for the vector of weights in lexicographic
  order, that is all weights of entity 1, then all weights of
  entity 2, and so on.
\end{inparaenum}

The function returns the variates in a list of class
\code{"portfolio"} with a \code{dim} attribute of length two. The list
contains all the individual claim amounts for each entity. Since every
element can be a vector, the object can be seen as a three-dimension
array with a third dimension of potentially varying length. The
function also returns a matrix of integers giving the classification
indexes of each entity in the portfolio.

The package also defines methods for four generic functions to easily
access key quantities for each entity of the simulated portfolio:
\begin{enumerate}
\item a method of \code{aggregate} to compute the aggregate claim
  amounts $S$;
\item a method of \code{frequency} to compute the number of claims
  $N$;
\item a method of \code{severity} (a generic function introduced by
  the package) to return the individual claim amounts $C_j$;
\item a method of \code{weights} to extract the weights matrix.
\end{enumerate}

In addition, all methods have a \code{classification} and a
\code{prefix} argument. When the first is \code{FALSE}, the
classification index columns are omitted from the result. The second
argument overrides the default column name prefix; see the
\code{simul.summaries} help page for details.

The following example illustrates these concepts in detail.

\begin{example}
  Consider the following compound hierarchical model:
  \begin{equation*}
    S_{ijt} = C_{ijt1} + \dots + C_{ijt N_{ijt}},
  \end{equation*}
  for
  $i = 1, \dots, I$,
  $j = 1, \dots, J_i$,
  $t = 1, \dots, n_{ij}$
  and with
  \begin{align*}
    N_{ijt}|\Lambda_{ij}, \Phi_i &\sim \text{Poisson}(w_{ijt} \Lambda_{ij}) &
    C_{ijtu}|\Theta_{ij}, \Psi_i &\sim \text{Lognormal}(\Theta_{ij}, 1)
    \notag \\
    \Lambda_{ij}|\Phi_i &\sim \text{Gamma}(\Phi_i, 1) &
    \Theta_{ij}|\Psi_i     &\sim N(\Psi_i, 1) \\
    \Phi_i &\sim \text{Exponential}(2) &
    \Psi_i &\sim N(2, 0.1). \notag
  \end{align*}
  Using as convention to number the data level 0, the above is a
  two-level compound hierarchical model.

  Assuming that $I = 2$, $J_1 = 4$, $J_2 = 3$, $n_{11} = \dots =
  n_{14} = 4$ and $n_{21} = n_{22} = n_{23} = 5$ and that weights are
  simply simulated from a uniform distribution on $(0.5, 2.5)$, then
  simulation of a data set with \code{simul} is achieved with:
<<echo=FALSE>>=
set.seed(3)
@
<<echo=TRUE>>=
nodes <- list(cohort = 2,
              contract = c(4, 3),
              year = c(4, 4, 4, 4, 5, 5, 5))
mf <- expression(cohort = rexp(2),
                 contract = rgamma(cohort, 1),
                 year = rpois(weights * contract))
ms <- expression(cohort = rnorm(2, sqrt(0.1)),
                 contract = rnorm(cohort, 1),
                 year = rlnorm(contract, 1))
wijt <- runif(31, 0.5, 2.5)
pf <- simul(nodes = nodes, model.freq = mf, model.sev = ms, weights = wijt)
@

  Object \code{pf} is a list of class \code{"portfolio"} containing,
  among other things, the aforementioned two-dimension list as element
  \code{data} and the classification matrix (subscripts $i$ and $j$)
  as element \code{classification}:
<<echo=TRUE>>=
class(pf)
pf$data
pf$classification
@


  The output of \code{pf\$data} is not much readable. Printing the
  results of \code{simul} like this would bring many users to wonder
  what \code{Numeric,\emph{n}} means. It is actually \proglang{R}'s
  way to specify that a given element in the list is a numeric vector
  of length $n$ --- the third dimension mentioned above. To ease
  reading, the \code{print} method for objects of class
  \code{"portfolio"} only prints the simulation model and the number
  of claims in each node:
<<echo=TRUE>>=
pf
@

  By default, the method of \code{aggregate} returns the values of
  $S_{ijt}$ in a regular matrix (subscripts $i$ and $j$ in the rows,
  subscript $t$ in the columns). The method has a \code{by} argument
  to get statistics for other groupings and a \code{FUN}
  argument to get statistics other than the sum:
<<echo=TRUE>>=
aggregate(pf)
aggregate(pf, by = c("cohort", "year"), FUN = mean)
@

  The method of \code{frequency} returns the values of $N_{ijt}$. It
  is mostly a wrapper for the \code{aggregate} method with the default
  \code{sum} statistic replaced by \code{length}. Hence, arguments
  \code{by} and \code{FUN} remain available:
<<echo=TRUE>>=
frequency(pf)
frequency(pf, by = "cohort")
@

  The method of \code{severity} returns the individual
  variates $C_{ijtu}$ in a matrix similar to those above, but with a
  number of columns equal to the maximum number of observations per
  entity,
  \begin{displaymath}
    \max_{i, j} \sum_{t = 1}^{n_{ij}} N_{ijt}.
  \end{displaymath}
  Thus, the original period of observation (subscript $t$) and the
  identifier of the severity within the period (subscript $u$) are
  lost and each variate now constitute a ``period'' of observation. For
  this reason, the method provides an argument \code{splitcol} in case
  one would like to extract separately the individual severities of
  one or more periods:
<<echo=TRUE>>=
severity(pf)
severity(pf, splitcol = 1)
@

  Finally, the weights matrix corresponding to the data in object
  \code{pf} is
<<echo=TRUE>>=
weights(pf)
@

  Combined with the argument \code{classification = FALSE}, the above
  methods can be used to easily compute loss ratios:
<<echo=TRUE>>=
aggregate(pf, classif = FALSE) / weights(pf, classif = FALSE)
@
  \qed
\end{example}

\begin{example}
  \cite{Scollnik:2001:MCMC} considers the following model for the
  simulation of claims frequency data in a Markov Chain Monte Carlo
  (MCMC) context:
  \begin{align*}
    S_{it}|\Lambda_i, \alpha, \beta &\sim \text{Poisson}(w_{ij} \Lambda_i) \\
    \Lambda_i|\alpha, \beta         &\sim \text{Gamma}(\alpha, \beta) \\
    \alpha &\sim \text{Gamma}(5, 5) \\
    \beta  &\sim \text{Gamma}(25, 1)
  \end{align*}
  for $i = 1, 2, 3$, $j = 1, \dots, 5$ and with weights $w_{it}$
  simulated from
  \begin{align*}
    w_{it}|a_i, b_i &\sim \text{Gamma}(a_i, b_i) \\
    a_i &\sim U(0, 100) \\
    b_i &\sim U(0, 100).
  \end{align*}
  Strictly speaking, this is not a hierarchical model since the random
  variables $\alpha$ and $\beta$ are parallel rather than nested.
  Nevertheless, with some minor manual intervention, function
  \code{simul} can simulate data from this model.

  First, one simulates the weights (in lexicographic order) with
<<echo=FALSE>>=
set.seed(123)
options(width = 55)
@
<<echo=TRUE>>=
wit <- rgamma(15, rep(runif(3, 0, 100), each = 5),
              rep(runif(3, 0, 100), each = 5))
@

  Second, one calls \code{simul} to simulate the frequency data. The key
  here consists in manually inserting the simulation of the shape and rate
  parameters of the gamma distribution in the model for $\Lambda_i$.
  Finally, wrapping the call to \code{simul} in \code{frequency} will
  immediately yield the matrix of observations:
<<echo=TRUE>>=
frequency(simul(list(entity = 3, year = 5),
      expression(entity = rgamma(rgamma(1, 5, 5), rgamma(1, 25, 1)),
          year = rpois(weights * entity)),
      weights = wit))
@
  \qed
\end{example}

One will find more examples of \code{simul} usage in the
\code{simulation} demo file. Function \code{simul} was used to
simulate the data in \cite{Goulet_cfs}.


\bibliography{actuar}

\end{document}
