Computational Actuarial Science

Last week, we’ve been through the book, completely, one last time, before sending it back to the publisher, with some comments and remarks, before publication ! So, this is it, the book will finally appear soon ! It was schedule for this week actually, but… you know. It should appear sometime by the end of May, or beginning of June. I will keep you posted on this blog.

A few months ago, we published with Christophe Dutang an ebook on the same topic, in French, online on cran.r-project.org/doc/contrib/. This contribution was based on lecture notes we had. When I’ve been asked to publish an English version, by John Kimmel, I was honored, but I thought it would be some kind of fraud if I write a book on that topic. I do know a bit of actuarial science, and a bit of R, but most of the advanced computations rely on others packages. Because I am extremely lazy, I did not try (so far) to edit my own package. I frequently publish some lines of codes on my blog, but nothing too serious.

So, for this book, I decided to ask those who actually did publish a package used in actuarial computations (or who did work previously on packages comparison for instance) to write a chapter, in this book. I am usually not a big fan of books with twenty contributors, because there is no coherence. So here, my task was to link all those chapters together, to make sure that notations are coherent, etc. Over 700 pages, that was difficult. And I asked all of them not only to illustrate actuarial concepts with some R code, but also give – if possible – some self written function, to understand the algorithm, but also some built-in functions. The goal was to explain the core of the algorithm. Some codes might not be efficient, but they help to understand how it could be possible to compute some actuarial quantities.

The first chapter is an Introduction (to the R language) I wrote with Rob Kaas (everyone in the actuarial community knows Rob ! not only as the Editor of Insurance: Mathematics & Economics. but also as a prolific author, including the popular textbook Modern Actuarial Risk Theory – Using R). The aim is to help those who might use another language for actuarial computation to understand the basics of the grammar, to read and write in R. I will probably publish a longer post, to explain the structure of that chapter. And show some codes.

• Methodology

The first section is a very general methodology section. It starts with Standard Statistical Inference by Christophe Dutang (Christophe is extremely active in the R community, as the maintainer of the Distributions task view page, for instance).Then, Ben Escoto (Ben works in the insurance industry, and launched the actuarial vignettes in R a few years ago) and myself, wrote a chapter which can be seen as an introduction to the Bayesian Philosophy for actuaries (I will give a talk on that topic at the R in Insurance conference this summer, so additional material will be online soon). With Stéphane Tufféry, we Statistical Learning (Stéphane published a Data Mining and Statistics for Decision Making a few years ago). Then, I wanted a chapter dedicated to Spatial Analysis. I did ask Renato Assunção, Marcelo Azevedo Costa, Marcos Oliveira Prates, and Luís Gustavo Silva e Silva to write that chapter (I met them a few years ago while I was visiting Renato in Belo Horizonte, while that started to work on spatial aspects of actuarial science). And finally, Eric Gilleland and Mathieu Ribatet wrote the chapter on Reinsurance and Extremal Events (both of them work on climate and extreme values, and they did publish a very interesting software review for extreme value analysis a few years ago).

• Life Insurance

For the section on life insurance, I asked Giorgio Spedicato to write the chapter on Life Contingencies (Giorgio is the author of the lifecontingencies package). Then Heather Booth, Rob J. Hyndman, and Leonie Tickle agreed to write the chapter on Prospective Life Tables (here we have a great match, with a demographer, an actuary, and… Rob. Every one who studied time series knows Rob. He is the author of the amazing forecast package, as well as the demography package, among many others. And he has a great blog too). To go further, Julien Tomas and Frédéric Planchet wrote the chapter on Prospective Mortality Tables and Portfolio Experience (both of them published the ELT – experience life tables – package a few months ago). And finally, there is a chapter on Survival Analysis by Frédéric Planchet  and Pierre-E. Thérond (they did publish a book – in French – on survival analysis for actuarial science, with examples in R).

• Finance

For the section on financial computations, Yohan Chalabi and Diethelm Würtz wrote two chapters, one on Stock Prices and Time Series and one on Portfolio Allocation (both of them have worked on the Rmetrics project, with the timeSeriesfArmafGarchfPortfolio packages). And Sergio S. Guirreri wrote a chapter on Yield Curves and Interest Rates Models (Sergio is the author of the YieldCurve package).

• Non-Life Insurance

Last, but not least, there is a section on non-life insuranceJean-Philippe Boucher (who published several articles on counts models) and myself, wrote the chapter on General Insurance Pricing. Then, I asked Katrien Antonio, Peng Shi, and Frank van Berkum to go further, with a chapter on Longitudinal Data and Experience Rating (I know Katrien from my PhD, and she was already working on that topic by that time… she did publish great surveys on that topic). And finally, Claims Reserving and IBNR is a chapter I wanted to write, because it’s a topic I love, but I asked Markus Gesmann to write it (Markus is known not only for his googleVis package, but also for the ChainLadder package – not to mention his awesome blog).

I will try to post some additional material on this blog, with R code (of course), graphs, and slides. And probably some pdfs with answers for the exercises. And all the datasets will be available in a CASdatasets package (online soon).

Stationarity of ARCH processes

In the context of AR(1) processes, we spent some time to explain what happens when $\phi$ is close to 1.

• if $\phi<1$ the process is stationary,
• if $\phi=1$ the process is a random walk
• if $\phi>1$ the process will explode

Again, random walks are extremely interesting processes, with puzzling properties. For instance,

$\text{Var}(X_{t+h}\vert X_t)\sim \sigma^2 h\rightarrow\infty$

as $h\rightarrow\infty$, and the process will cross the x-axis an infinite number of times…

Recently, in the MAT8181 course, we studied carefully properties of the ARCH(1) process, especially when $\alpha\sim 1$. And again, what we get might be puzzling.

Consider some ARCH(1) process $(\varepsilon_t)$, with a Gaussian noise, i.e.

$\varepsilon_t=\sigma_t\cdot \eta_t$

where

$\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2$

and $(\eta_t)$ is a sequence of i.i.d. $\mathcal{N}(0,1)$ variables. Here both $\omega$ and $\alpha$ have to be positive.

Recall that $\mathbb{E}(\varepsilon_t)=0$ since $\mathbb{E}(\eta_t)=0$. Further

$\text{var}(\varepsilon_t)=\mathbb{E}(\varepsilon_t^2)=\omega+\alpha \mathbb{E}(\varepsilon_{t-1}^2)$

since $\mathbb{E}(\eta_t^2)=1$, so the variance exists, and is constant only if $\alpha\in(0,1)$, and in that case

$\sigma^2=\text{var}(\varepsilon_t)=\frac{\omega}{1-\alpha}\in[0,+\infty)$

Further, if $3\alpha^2<1$, then the fourth moment can be obtained,

$\mathbb{E}(\varepsilon_t^4)=\frac{3\omega^2}{1-\alpha^2}\frac{1-\alpha^2}{1-3\alpha^2}$

since$\mathbb{E}(\eta_t^4)=3$. Now, if we get back on the property obtained while studying the variance, what does that mean if $\alpha=1$, or $\alpha>1$ ?

If we look at simulations, we can generate an ARCH(1) process with $\alpha=2$ for instance.

> n=600
> a=2
> w=0.2
> set.seed(1)
> eta=rnorm(n)
> epsilon=rnorm(n)
> sigma2=rep(w,n)
> for(t in 2:n){
+ sigma2[t]=w+a*epsilon[t-1]^2
+ epsilon[t]=eta[t]*sqrt(sigma2[t])
+ }
> plot(epsilon,type="l")

In order to understand what’s going on, we should keep in mind that, what we good is that $\alpha$ has to lie in $(0,1)$ to be able to compute the second moment of $(\varepsilon_t)$. But it is possible to have a stationary process with infinite variance. And actually, this is what we have here.

Write

$\sigma^2_t=\omega+\alpha \varepsilon_{t-1}^2 = \omega+[\alpha \eta_{t-1}^2] \sigma_{t-1}^2$

and them, iterate

$\sigma^2_t= \omega+[\alpha \eta_{t-1}^2] \left( \omega+[\alpha \eta_{t-2}^2] \sigma_{t-2}^2\right)$

and iterate again, and again, and again…

$\sigma^2_t= \underbrace{\omega\left[ 1+\sum_{i=1}^h[\alpha \eta_{t-1}^2]\cdots [\alpha \eta_{t-i}^2] \right]}_{\Sigma_t(h)} +[\alpha \eta_{t-1}^2]\cdots[\alpha \eta_{t-h-1}^2]\sigma_{t-h-1}^2$

where

$\Sigma_t(h)=\sum_{i=1}^h\underbrace{[\alpha \eta_{t-1}^2]\cdots [\alpha \eta_{t-i}^2]}_{u_i}$

Here, we have a sum of positive terms, and we can use the so-called Cauchy rule: define

$\lambda=\text{limsup}\{ u_n^{1/n}\}$

then, if $\lambda<1$, the series $\sum u_n$ converges. Here,

$u_n^{1/n}=\left[[\alpha \eta_{t-1}^2]\cdots [\alpha \eta_{t-n}^2]\right]^{1/n}$

which can also be written

$u_n^{1/n}=\exp\left[\frac{1}{n}\sum_{i=1}^n\log[\alpha \eta_{t-i}^2]\right]$

and from the law of large numbers, since we have here a sum of i.i.d. terms,

$u_n^{1/n}\rightarrow\exp\left[\mathbb{E}(\log[\alpha \eta^2])\right]$

So, if $\exp\left[\mathbb{E}(\log[\alpha \eta^2])\right]<1$, then $\Sigma_t(h)$ will have a limit when $h$ goes to infinity.

The condition above can be written

$\gamma=\mathbb{E}(\log[\alpha \eta^2]<0$

which is called Lyapunov coefficient.

The equation

$\exp\left[\mathbb{E}(\log[\alpha \eta^2])\right]=\alpha\exp\left[-\mathbb{E}(\log[\eta^2])\right]<1$

is a condition on $\alpha$.

In the case where $\eta\sim\mathcal{N}(0,1)$, the numerical value of this upper bound is 3.56.

> 1/exp(mean(log(rnorm(1e7)^2)))
[1] 3.562517

In that case ($\gamma<0$), the variance may be infinite, but the series is stationary. On the other hand, if $\gamma>0$, then $\varepsilon_t^2$ will go to infinity almost surely, as $t$ goes to infinity.

But in order to observe this difference, we need a lot of observations. For instance, with $\alpha=0.8$,

and $\alpha=1.2$,

we can easily see a difference. I do not say that it’s easy to see that the distribution above has an infinite variance, but still. Actually, if we consider Hill’s plot on the series above, on the tails of positive $\varepsilon_t$‘s

> library(evir)
> hill(epsilon)

or on the tails of negative $\varepsilon_t$‘s

> hill(-epsilon)

we can see that the tail index is (strictly) smaller than 2 (meaning that the moment of order 2 does not exist).

Why is it puzzling? Maybe because here, $(\varepsilon_t)$ is not weakly stationary (in the $L^2$ sense), but it is strongly stationary. Which is not the usual way weak and strong are related. This might be why we will not call this strong stationarity, but strict.