Showing posts with label Actuarial. Show all posts
Showing posts with label Actuarial. Show all posts
Fitting a distribution in Stan from scratch
Last week the French National Institute of Health and Medical Research (Inserm) organised with the Stan Group a training programme on Bayesian Inference with Stan for Pharmacometrics in Paris.
Daniel Lee and Michael Betancourt, who run the course over three days, are not only members of Stan's development team, but also excellent teachers. Both were supported by Eric Novik, who gave an Introduction to Stan at the Paris Dataiku User Group last week as well.
I have been playing around with Stan on and off for some time, but as Eric pointed out to me, Stan is not that kind of girl(boy?). Indeed, having spent three days working with Stan has revitalised my relationship. Getting down to the basics has been really helpful and I shall remember, Stan is not drawing samples from a distribution. Instead, it is calculating the joint distribution function (in log space), and evaluating the probability distribution function (in log space).
Thus, here is a little example of fitting a set of random numbers in R to a Normal distribution with Stan. Yet, instead of using the built-in functions for the Normal distribution, I define the log probability function by hand, which I will use in the model block as well, and even generate a random sample, starting with a uniform distribution. However, I do use pre-defined distributions for the priors.
Why do I want to do this? This will be a template for the day when I have to use a distribution, which is not predefined in Stan, e.g. the actuar package has some interesting candidates.
I then use the Stan script to fit the data, i.e. to find the the parameters \(\mu\) and \(\sigma\), assuming that the data was generated by a Gaussian process.
The posterior parameter distributions include both \(\mu\) and \(\sigma\) in the 95% credible interval. The distribution of posterior predictive check (
Daniel Lee and Michael Betancourt, who run the course over three days, are not only members of Stan's development team, but also excellent teachers. Both were supported by Eric Novik, who gave an Introduction to Stan at the Paris Dataiku User Group last week as well.
 |
| Eric Kramer (Dataiku), Daniel Lee, Eric Novik & Michael Betancourt (Stan Group) |
I have been playing around with Stan on and off for some time, but as Eric pointed out to me, Stan is not that kind of girl(boy?). Indeed, having spent three days working with Stan has revitalised my relationship. Getting down to the basics has been really helpful and I shall remember, Stan is not drawing samples from a distribution. Instead, it is calculating the joint distribution function (in log space), and evaluating the probability distribution function (in log space).
Thus, here is a little example of fitting a set of random numbers in R to a Normal distribution with Stan. Yet, instead of using the built-in functions for the Normal distribution, I define the log probability function by hand, which I will use in the model block as well, and even generate a random sample, starting with a uniform distribution. However, I do use pre-defined distributions for the priors.
Why do I want to do this? This will be a template for the day when I have to use a distribution, which is not predefined in Stan, e.g. the actuar package has some interesting candidates.
Testing
I start off by generating fake data, a sample of 100 random numbers drawn from a Normal distribution with a mean of 4 and a standard deviation of 2. Note, the sample mean of the 100 figures is 4.2 and not 4. |
| Histogram of 100 random numbers drawn from N(4,2). |
 |
| Traceplot of 4 chains, including warm-up phase |
 |
| Histograms of posterior parameter and predictive samples |
 |
| Comparison of the emperical distributions |
y_ppc) is wider, taking into account the uncertainty of the parameters. The interquartile range and mean of my initial fake data and the sample of the posterior predictive distribution look very similar. That's good, my model generates data, which looks like the original data.Bayesian Mixer Meetup
Btw, tonight we have the 4th Bayesian Mixer Meetup in London.Session Info
R version 3.3.1 (2016-06-21)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.12 (Sierra)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] MASS_7.3-45 rstan_2.12.1 StanHeaders_2.12.0 ggplot2_2.1.0
loaded via a namespace (and not attached):
[1] Rcpp_0.12.7 codetools_0.2-14 digest_0.6.10 grid_3.3.1
[5] plyr_1.8.4 gtable_0.2.0 stats4_3.3.1 scales_0.4.0
[9] labeling_0.3 tools_3.3.1 munsell_0.4.3 inline_0.3.14
[13] colorspace_1.2-6 gridExtra_2.2.1
Actuarial
,
Bayesian
,
fit distribution
,
R
,
Stan
Posterior predictive output with Stan
I continue my Stan experiments with another insurance example. Here I am particular interested in the posterior predictive distribution from only three data points. Or, to put it differently I have a customer of three years and I'd like to predict the expected claims cost for the next year to set or adjust the premium.
The example is taken from section 16.17 in Loss Models: From Data to Decisions [1]. Some time ago I used the same example to get my head around a Bayesian credibility model.
Suppose the claims likelihood distribution is believed to follow an exponential distribution for a given parameter \(\Theta\). The prior parameter distribution on \(\Theta\) is assumed to be a gamma distribution with parameters \(\alpha=4, \beta=1000\):
\[\begin{aligned}\Theta & \sim \mbox{Gamma}(\alpha, \beta)\\
\ell_i & \sim \mbox{Exp}(\Theta) , \; \forall i \in N
\end{aligned}\]In this case the predictive distribution is a Pareto II distribution with density \(f(x) = \frac{\alpha \beta^\alpha}{(x+\beta)^{\alpha+1}}\) and a mean of \(\frac{\beta}{\alpha-1}=\,\)333.33.
I have three independent observations, namely losses of $100, $950 and $450. The posterior predictive expected loss is $416.67 and can be derived analytical, as shown in my earlier post. Now let me reproduce the answer with Stan as well.
Implementing the model in Stan is straightforward and I follow the same steps as in my simple example of last week. However, here I am also interested in the posterior predictive distribution, hence I add a
The output shows a simulated predictive mean of $416.86, close to the analytical answer. I can also read out that the 75%ile of the posterior predictive distribution is a loss of $542 vs. $414 from the prior predictive. That means every four years I shouldn't be surprised to observe a loss in excess of $500. Further I note that 90% of losses are expected to be less than $950, or in other words the observation in my data may reflect the outcome of an event with a 1 in 10 return period.
Comparing the sampling output from Stan with the analytical output gives me some confidence that I am doing the 'right thing'.
The example is taken from section 16.17 in Loss Models: From Data to Decisions [1]. Some time ago I used the same example to get my head around a Bayesian credibility model.
Suppose the claims likelihood distribution is believed to follow an exponential distribution for a given parameter \(\Theta\). The prior parameter distribution on \(\Theta\) is assumed to be a gamma distribution with parameters \(\alpha=4, \beta=1000\):
\[\begin{aligned}\Theta & \sim \mbox{Gamma}(\alpha, \beta)\\
\ell_i & \sim \mbox{Exp}(\Theta) , \; \forall i \in N
\end{aligned}\]In this case the predictive distribution is a Pareto II distribution with density \(f(x) = \frac{\alpha \beta^\alpha}{(x+\beta)^{\alpha+1}}\) and a mean of \(\frac{\beta}{\alpha-1}=\,\)333.33.
I have three independent observations, namely losses of $100, $950 and $450. The posterior predictive expected loss is $416.67 and can be derived analytical, as shown in my earlier post. Now let me reproduce the answer with Stan as well.
Implementing the model in Stan is straightforward and I follow the same steps as in my simple example of last week. However, here I am also interested in the posterior predictive distribution, hence I add a
generated quantities code block. The output shows a simulated predictive mean of $416.86, close to the analytical answer. I can also read out that the 75%ile of the posterior predictive distribution is a loss of $542 vs. $414 from the prior predictive. That means every four years I shouldn't be surprised to observe a loss in excess of $500. Further I note that 90% of losses are expected to be less than $950, or in other words the observation in my data may reflect the outcome of an event with a 1 in 10 return period.
Comparing the sampling output from Stan with the analytical output gives me some confidence that I am doing the 'right thing'.
References
[1] Klugman, S. A., Panjer, H. H. & Willmot, G. E. (2004), Loss Models: From Data to Decisions, Wiley Series in Probability and Statistics.Session Info
R version 3.2.0 (2015-04-16)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.3 (Yosemite)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lattice_0.20-31 actuar_1.1-8 rstan_2.6.0 inline_0.3.14
[5] Rcpp_0.11.6
loaded via a namespace (and not attached):
[1] tools_3.2.0 codetools_0.2-11 grid_3.2.0 stats4_3.2.0
First steps with ChainLadder: Import triangle from Excel into R
Taking the first step is often the hardest: getting data from Excel into R.
Suppose you would like to use the
How do you get the triangle into R and execute a reserving function, such as
Well, there are many ways to do this and the ChainLadder package vignette, as well as the R manual on Data Import/Export has all of the details, but here is a quick and dirty solution using a CSV-file.
Open a new Excel workbook and copy your triangle into cell A1, with the first column being the accident or origin period and the first row describing the development period or age. You find an example CSV-file on GitHub.
Make sure that your triangle has no formatting, such a commas to separate thousands, as Excel will save those cell as characters. Now open R and go through the following commands:
To copy the full triangle back into Excel you can use the clipboard:
Go back to Excel and hit <Ctrl> + V on your keyboard to paste the data into R.
For more details see the package vignette and Dan's post on pasting triangles from Excel into R via the clipboard.
If you are after a thorough overview of R in insurance take a look at the book Computational Actuarial Science with R.
Finally, join the Special Interest Group on using R in actuarial science and insurance to share your questions and answers.
Suppose you would like to use the
ChainLadder package to forecast future claims payments for a run-off triangle that you have stored in Excel.
How do you get the triangle into R and execute a reserving function, such as
MackChainLadder?Well, there are many ways to do this and the ChainLadder package vignette, as well as the R manual on Data Import/Export has all of the details, but here is a quick and dirty solution using a CSV-file.
Open a new Excel workbook and copy your triangle into cell A1, with the first column being the accident or origin period and the first row describing the development period or age. You find an example CSV-file on GitHub.
Make sure that your triangle has no formatting, such a commas to separate thousands, as Excel will save those cell as characters. Now open R and go through the following commands:
# The first command will open a window and
# ask you to select your CSV-file
myCSVfile <- file.choose()
# Read file into R
dat <- read.csv(file=myCSVfile)
# use read.csv2 if semicolons are used as a separator
# likely to be the case if you are in continental Europe
dat # to see your data
AY X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
1 1981 5012 8269 10907 11805 13539 16181 18009 18608 18662 18834
2 1982 106 4285 5396 10666 13782 15599 15496 16169 16704 NA
3 1983 3410 8992 13873 16141 18735 22214 22863 23466 NA NA
4 1984 5655 11555 15766 21266 23425 26083 27067 NA NA NA
5 1985 1092 9565 15836 22169 25955 26180 NA NA NA NA
6 1986 1513 6445 11702 12935 15852 NA NA NA NA NA
7 1987 557 4020 10946 12314 NA NA NA NA NA NA
8 1988 1351 6947 13112 NA NA NA NA NA NA NA
9 1989 3133 5395 NA NA NA NA NA NA NA NA
10 1990 2063 NA NA NA NA NA NA NA NA NA
# Load the ChainLadder package
library(ChainLadder)
# Ignore first column which holds accident year information
tri <- dat[,-1]
# Convert to matrix
tri <- as.matrix(tri)
# Add dimension names
dimnames(tri) <- list(origin=dat[,1], dev=1:ncol(tri))
# Convert into a triangle class
tri <- as.triangle(tri)
tri
dev
origin 1 2 3 4 5 6 7 8 9 10
1981 5012 8269 10907 11805 13539 16181 18009 18608 18662 18834
1982 106 4285 5396 10666 13782 15599 15496 16169 16704 NA
1983 3410 8992 13873 16141 18735 22214 22863 23466 NA NA
1984 5655 11555 15766 21266 23425 26083 27067 NA NA NA
1985 1092 9565 15836 22169 25955 26180 NA NA NA NA
1986 1513 6445 11702 12935 15852 NA NA NA NA NA
1987 557 4020 10946 12314 NA NA NA NA NA NA
1988 1351 6947 13112 NA NA NA NA NA NA NA
1989 3133 5395 NA NA NA NA NA NA NA NA
1990 2063 NA NA NA NA NA NA NA NA NA
MackChainLadder function:M <- MackChainLadder(tri, est.sigma = "Mack")
M
Latest Dev.To.Date Ultimate IBNR Mack.S.E CV(IBNR)
1981 18,834 1.000 18,834 0 0 NaN
1982 16,704 0.991 16,858 154 206 1.339
1983 23,466 0.974 24,083 617 623 1.010
1984 27,067 0.943 28,703 1,636 747 0.457
1985 26,180 0.905 28,927 2,747 1,469 0.535
1986 15,852 0.813 19,501 3,649 2,002 0.549
1987 12,314 0.694 17,749 5,435 2,209 0.406
1988 13,112 0.546 24,019 10,907 5,358 0.491
1989 5,395 0.336 16,045 10,650 6,333 0.595
1990 2,063 0.112 18,402 16,339 24,566 1.503
Totals
Latest: 160,987.00
Dev: 0.76
Ultimate: 213,122.23
IBNR: 52,135.23
Mack S.E.: 26,909.01
CV(IBNR): 0.52
To copy the full triangle back into Excel you can use the clipboard:
write.table(M$FullTriangle, file="clipboard", sep="\t"). Go back to Excel and hit <Ctrl> + V on your keyboard to paste the data into R.
For more details see the package vignette and Dan's post on pasting triangles from Excel into R via the clipboard.
If you are after a thorough overview of R in insurance take a look at the book Computational Actuarial Science with R.
Finally, join the Special Interest Group on using R in actuarial science and insurance to share your questions and answers.
Session Info
R version 3.1.2 (2014-10-31)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] ChainLadder_0.1.9 systemfit_1.1-14 lmtest_0.9-33
[4] zoo_1.7-11 car_2.0-21 Matrix_1.1-4
loaded via a namespace (and not attached):
[1] acepack_1.3-3.3 actuar_1.1-6 cluster_1.15.3
[4] foreign_0.8-61 Formula_1.1-2 grid_3.1.2
[7] Hmisc_3.14-5 lattice_0.20-29 latticeExtra_0.6-26
[10] MASS_7.3-35 nnet_7.3-8 plyr_1.8.1
[13] RColorBrewer_1.0-5 Rcpp_0.11.3 reshape2_1.4
[16] rpart_4.1-8 sandwich_2.3-2 splines_3.1.2
[19] statmod_1.4.20 stringr_0.6.2 survival_2.37-7
[22] tools_3.1.2 tweedie_2.2.1
Actuarial
,
ChainLadder
,
Insurance
,
R
Not only verbs but also believes can be conjugated
Following on from last week, where I presented a simple example of a Bayesian network with discrete probabilities to predict the number of claims for a motor insurance customer, I will look at continuous probability distributions today. Here I follow example 16.17 in Loss Models: From Data to Decisions [1].
Suppose there is a class of risks that incurs random losses following an exponential distribution (density \(f(x) = \Theta {e}^{- \Theta x}\)) with mean \(1/\Theta\). Further, I believe that \(\Theta\) varies according to a gamma distribution (density \(f(x)= \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha \,-\, 1} e^{- \beta x } \)) with shape \(\alpha=4\) and rate \(\beta=1000\).
In the same way as I had good and bad driver in my previous post, here I have clients with different characteristics, reflected by the gamma distribution.
The textbook tells me that the unconditional mixed distribution of an exponential distribution with parameter \(\Theta\), whereby \(\Theta\) has a gamma distribution, is a Pareto II distribution (density \(f(x) = \frac{\alpha \beta^\alpha}{(x+\beta)^{\alpha+1}}\)) with parameters \(\alpha,\, \beta\). Its k-th moment is given in the general case by
\[
E[X^k] = \frac{\beta^k\Gamma(k+1)\Gamma(\alpha - k)}{\Gamma(\alpha)},\; -1 < k < \alpha. \] Thus, I can calculate the prior expected loss (\(k=1\)) as \(\frac{\beta}{\alpha-1}=\,\)333.33.
Now suppose I have three independent observations, namely losses of $100, $950 and $450 over the last 3 years. The mean loss is $500, which is higher than the $333.33 of my model.
Question: How should I update my belief about the client's risk profile to predict the expected loss cost for year 4 given those 3 observations?
Visually I can regard this scenario as a graph, with evidence set for years 1 to 3 that I want to propagate through to year 4.
Suppose there is a class of risks that incurs random losses following an exponential distribution (density \(f(x) = \Theta {e}^{- \Theta x}\)) with mean \(1/\Theta\). Further, I believe that \(\Theta\) varies according to a gamma distribution (density \(f(x)= \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha \,-\, 1} e^{- \beta x } \)) with shape \(\alpha=4\) and rate \(\beta=1000\).
In the same way as I had good and bad driver in my previous post, here I have clients with different characteristics, reflected by the gamma distribution.
The textbook tells me that the unconditional mixed distribution of an exponential distribution with parameter \(\Theta\), whereby \(\Theta\) has a gamma distribution, is a Pareto II distribution (density \(f(x) = \frac{\alpha \beta^\alpha}{(x+\beta)^{\alpha+1}}\)) with parameters \(\alpha,\, \beta\). Its k-th moment is given in the general case by
\[
E[X^k] = \frac{\beta^k\Gamma(k+1)\Gamma(\alpha - k)}{\Gamma(\alpha)},\; -1 < k < \alpha. \] Thus, I can calculate the prior expected loss (\(k=1\)) as \(\frac{\beta}{\alpha-1}=\,\)333.33.
Now suppose I have three independent observations, namely losses of $100, $950 and $450 over the last 3 years. The mean loss is $500, which is higher than the $333.33 of my model.
Question: How should I update my belief about the client's risk profile to predict the expected loss cost for year 4 given those 3 observations?
Visually I can regard this scenario as a graph, with evidence set for years 1 to 3 that I want to propagate through to year 4.
Predicting claims with a Bayesian network
Here is a little Bayesian Network to predict the claims for two different types of drivers over the next year, see also example 16.15 in [1].
Let's assume there are good and bad drivers. The probabilities that a good driver will have 0, 1 or 2 claims in any given year are set to 70%, 20% and 10%, while for bad drivers the probabilities are 50%, 30% and 20% respectively.
Further I assume that 75% of all drivers are good drivers and only 25% would be classified as bad drivers. Therefore the average number of claims per policyholder across the whole customer base would be:
To answer the above question I present the data here as a Bayesian Network using the
Next, I set the client's evidence (0 claims in year one and 1 claim in year two) and propagate these back through my network to estimate the probabilities that the customer is either a good (73.68%) or a bad (26.32%) driver. Knowing that a good driver has on overage 0.4 claims a year and a bad driver 0.7 claims I predict the number of claims for my customer with the given claims history as 0.4789.
Alternatively I could have added a third node for year 3 and queried the network for the probabilities of 0, 1 or 2 claims given that the customer had zero claims in year 1 and one claim in year 2. The sum product of the number of claims and probabilities gives me again an expected claims number of 0.4789.
[2] Søren Højsgaard (2012). Graphical Independence Networks with the gRain Package for R. Journal of Statistical Software, 46(10), 1-26. URL http://www.jstatsoft.org/v46/i10/
Let's assume there are good and bad drivers. The probabilities that a good driver will have 0, 1 or 2 claims in any given year are set to 70%, 20% and 10%, while for bad drivers the probabilities are 50%, 30% and 20% respectively.
Further I assume that 75% of all drivers are good drivers and only 25% would be classified as bad drivers. Therefore the average number of claims per policyholder across the whole customer base would be:
0.75*(0*0.7 + 1*0.2 + 2*0.1) + 0.25*(0*0.5 + 1*0.3 + 2*0.2) = 0.475
To answer the above question I present the data here as a Bayesian Network using the
gRain package [2]. I start with the contingency probability tables for the driver type and the conditional probabilities for 0, 1 and 2 claims in year 1 and 2. As I assume independence between the years I set the same probabilities. I can now review my model as a mosaic plot (above) and as a graph (below) as well.
Next, I set the client's evidence (0 claims in year one and 1 claim in year two) and propagate these back through my network to estimate the probabilities that the customer is either a good (73.68%) or a bad (26.32%) driver. Knowing that a good driver has on overage 0.4 claims a year and a bad driver 0.7 claims I predict the number of claims for my customer with the given claims history as 0.4789.
Alternatively I could have added a third node for year 3 and queried the network for the probabilities of 0, 1 or 2 claims given that the customer had zero claims in year 1 and one claim in year 2. The sum product of the number of claims and probabilities gives me again an expected claims number of 0.4789.
References
[1] Klugman, S. A., Panjer, H. H. & Willmot, G. E. (2004), Loss Models: From Data to Decisions, Wiley Series in Proability and Statistics.[2] Søren Højsgaard (2012). Graphical Independence Networks with the gRain Package for R. Journal of Statistical Software, 46(10), 1-26. URL http://www.jstatsoft.org/v46/i10/
Session Info
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] Rgraphviz_2.6.0 gRain_1.2-2 gRbase_1.6-12 graph_1.40.0
loaded via a namespace (and not attached):
[1] BiocGenerics_0.8.0 igraph_0.6.6 lattice_0.20-24 Matrix_1.1-0
[5] parallel_3.0.2 RBGL_1.38.0 stats4_3.0.2 tools_3.0.2Claims Inflation - a known unknown
Over the last year I worked with two colleagues of mine on the subject of inflation and claims inflation in particular. I didn't expect it to be such a challenging topic, but we ended up with more questions than answers. The key question and biggest challenge is to define what inflation, or indeed claims inflation actually is and how to measure it. We published a summary of our thoughts and findings in this month's issue of The Actuary.
Last year's discussion about the differences between the retail price index (RPI) and consumer price index (CPI) in the UK only exemplified the challenge. The economist Tim Harford illustrated the differences between the RPI and CPI with a simple example of price changes for a shirt and blouse in his Radio 4 programme More or Less. The radio podcast is still available from the BBC. Start listening after about 18 minutes into the show.
Last year's discussion about the differences between the retail price index (RPI) and consumer price index (CPI) in the UK only exemplified the challenge. The economist Tim Harford illustrated the differences between the RPI and CPI with a simple example of price changes for a shirt and blouse in his Radio 4 programme More or Less. The radio podcast is still available from the BBC. Start listening after about 18 minutes into the show.
ChainLadder 0.1.5-6 released on CRAN
Last week we released version 0.1.5-6 of the ChainLadder package on CRAN. The ChainLadder package provides statistical models, which are typically used for the estimation of outstanding claims reserves in general insurance. The package vignette gives an overview of the package functionality.
Since the last CRAN release Dan Murphy added new features to the
Please get in touch if you would like to collaborate or find any issues or bugs.
Join me at the first R in Insurance conference at Cass Business School in London, 15 July 2013.
 |
Output of plot(MackChainLadder(GenIns)) |
Since the last CRAN release Dan Murphy added new features to the
MackChainLadder function and we fixed a bug in BootChainLadder. Here are he details:New Features
- The list output of the MackChainLadder function now includes the parameter risk and process risk breakdowns of the total risk estimate for the sum of projected losses across all origin years by development age.
- The Mack Method's recursive parameter risk calculation now enables Mack's original two-term formula (the default) and optionally the three-term formula found in Murphy's 1994 paper and in the 2006 paper by Buchwalder, Bühlmann, Merz, and Wüthrich.
- A few more Mack Method examples.
Bug Fixes
- The phi-scaling factor in BootChainLadder was incorrect. Instead of calculating the number of data items in the upper left triangle as n*(n+1)/2, n*(n-1)/2 was used. Thanks to Thomas Girodot for reporting this bug.
Please get in touch if you would like to collaborate or find any issues or bugs.
Join me at the first R in Insurance conference at Cass Business School in London, 15 July 2013.
Submit a talk for the first R in Insurance conference
The registration for the first R in Insurance is open and there is still time to submit a talk / lightning talk.
The conference will take place at Cass Business School in London on Monday, 15 July 2013. This is the Monday following the useR! 2013 conference in Spain. Thus, if you come from overseas to Spain, why not stop in London on your way back?
All further information and registration details are available on the Cass Business School conference site.
 |
All further information and registration details are available on the Cass Business School conference site.
Actuarial
,
Conference
,
Insurance
,
R
,
R in Insurance
Reserving based on log-incremental payments in R, part III
This is the third post about Christofides' paper on Regression models based on log-incremental payments [1]. The first post covered the fundamentals of Christofides' reserving model in sections A - F, the second focused on a more realistic example and model reduction of sections G - K. Today's post will wrap up the paper with sections L - M and discuss data normalisation and claims inflation.
I will use the same triangle of incremental claims data as introduced in my previous post. The final model had three parameters for origin periods and two parameters for development periods. It is possible to reduce the model further as Christofides illustrates in section L onwards by using an inflation index to bring all claims payments to current value and a claims volume adjustment or weight for each origin period to normalise the triangle.
In his example Christofides uses claims volume adjustments for the origin years and an earning or inflation index for the different payment calendar years. The claims volume adjustments aims to normalise the triangle for similar exposures across origin periods, while the earnings index, which measures largely wages and other forms of compensations, is used as a first proxy for claims inflation. Note that the earnings index shows significant year on year changes from 5% to 9%. Barnett and Zehnwirth [2] would probably recommend to add further parameters for the calendar year effects to the model.
I will use the same triangle of incremental claims data as introduced in my previous post. The final model had three parameters for origin periods and two parameters for development periods. It is possible to reduce the model further as Christofides illustrates in section L onwards by using an inflation index to bring all claims payments to current value and a claims volume adjustment or weight for each origin period to normalise the triangle.
In his example Christofides uses claims volume adjustments for the origin years and an earning or inflation index for the different payment calendar years. The claims volume adjustments aims to normalise the triangle for similar exposures across origin periods, while the earnings index, which measures largely wages and other forms of compensations, is used as a first proxy for claims inflation. Note that the earnings index shows significant year on year changes from 5% to 9%. Barnett and Zehnwirth [2] would probably recommend to add further parameters for the calendar year effects to the model.
# Page D5.36
ClaimsVolume <- data.frame(origin=0:6,
volume.index=c(1.43, 1.45, 1.52, 1.35, 1.29, 1.47, 1.91))
# Page D5.36
EarningIndex <- data.frame(cal=0:6,
earning.index=c(1.55, 1.41, 1.3, 1.23, 1.13, 1.05, 1))
# Year on year changes
round((1-EarningIndex$earning.index[-1]/EarningIndex$earning.index[-7]),2)
# [1] 0.09 0.08 0.05 0.08 0.07 0.05
dat <- merge(merge(dat, ClaimsVolume), EarningIndex)
# Normalise data for volume and earnings
dat$logvalue.ind.inf <- with(dat, log(value/volume.index*earning.index))
with(dat, interaction.plot(dev, origin, logvalue.ind.inf))
points(1+dat$dev, dat$logvalue.ind.inf, pch=16, cex=0.8)
Actuarial
,
Barnett Zehnwirth
,
ChainLadder
,
IBNR
,
Insurance
,
linear model
,
log-incremental
,
R
,
reserving
,
Risk
,
Tutorials
Reserving based on log-incremental payments in R, part II
Following on from last week's post I will continue to go through the paper Regression models based on log-incremental payments by Stavros Christofides [1]. In the previous post I introduced the model from the first 15 pages up to section F. Today I will progress with sections G to K which illustrate the model with a more realistic incremental claims payments triangle from a UK Motor Non-Comprehensive account:
Before I plot the data, I will transform the triangle into a data frame and add extra columns:
Indeed the origin years 1 to 4 (rows 2 to 5) look quite similar and the decay of claims in development year direction appears to be linear on a log-scale from development year 1 onwards.
Based on those observations Christofides suggests two models; the first one will have a unique level for each origin year and a unique level for the zero development period. The parameters for development periods 1 to 6 are assumed to follow a linear relationship with the same slope \(s\):
\begin{align}
\ln(P_{ij}) & = Y_{ij} = a_i + d_j + \epsilon_{ij}
&\mbox{for } i,\,j \mbox{ from } 0 \mbox{ to } 6\\
\mbox{where } d_0 &= d,\quad d_j = s \cdot j
&\mbox{for } j > 0
\end{align}and \(\epsilon_{ij} \sim N(0, \sigma^2)\). The second model will be a reduced version of the above with only two levels for the origin years 5 and 6. Hence, I add four more columns to my data frame:
# Page D5.17
tri <- t(matrix(
c(3511, 3215, 2266, 1712, 1059, 587, 340,
4001, 3702, 2278, 1180, 956, 629, NA,
4355, 3932, 1946, 1522, 1238, NA, NA,
4295, 3455, 2023, 1320, NA, NA, NA,
4150, 3747, 2320, NA, NA, NA, NA,
5102, 4548, NA, NA, NA, NA, NA,
6283, NA, NA, NA, NA, NA, NA), nc=7))
Before I plot the data, I will transform the triangle into a data frame and add extra columns:
m <- dim(tri)[1]; n <- dim(tri)[2]
dat <- data.frame(
origin=rep(0:(m-1), n),
dev=rep(0:(n-1), each=m),
value=as.vector(tri))
## Add dimensions as factors
dat <- with(dat, data.frame(origin, dev, cal=origin+dev,
value, logvalue=log(value),
originf=factor(origin),
devf=as.factor(dev),
calf=as.factor(origin+dev)))
interaction.plot of the stats package does an excellent job for this:op <- par(mfrow=c(2,1), mar=c(4,4,2,2))
with(dat, interaction.plot(x.factor=dev, trace.factor=origin,
response=value))
points(dat$devf, dat$value, pch=16, cex=0.5)
with(dat, interaction.plot(x.factor=dev, trace.factor=origin,
response=logvalue))
points(dat$devf, dat$logvalue, pch=16, cex=0.5)
par(op)
Based on those observations Christofides suggests two models; the first one will have a unique level for each origin year and a unique level for the zero development period. The parameters for development periods 1 to 6 are assumed to follow a linear relationship with the same slope \(s\):
\begin{align}
\ln(P_{ij}) & = Y_{ij} = a_i + d_j + \epsilon_{ij}
&\mbox{for } i,\,j \mbox{ from } 0 \mbox{ to } 6\\
\mbox{where } d_0 &= d,\quad d_j = s \cdot j
&\mbox{for } j > 0
\end{align}and \(\epsilon_{ij} \sim N(0, \sigma^2)\). The second model will be a reduced version of the above with only two levels for the origin years 5 and 6. Hence, I add four more columns to my data frame:
Actuarial
,
Barnett Zehnwirth
,
ChainLadder
,
IBNR
,
Insurance
,
linear model
,
log-incremental
,
R
,
reserving
,
Risk
,
Tutorials
Reserving based on log-incremental payments in R, part I
A recent post on the PirateGrunt blog on claims reserving inspired me to look into the paper Regression models based on log-incremental payments by Stavros Christofides [1], published as part of the Claims Reserving Manual (Version 2) of the Institute of Actuaries.
The paper is available together with a spread sheet model, illustrating the calculations. It is very much based on ideas by Barnett and Zehnwirth, see [2] for a reference. However, doing statistical analysis in a spread sheet programme is often cumbersome. I will go through the first 15 pages of Christofides' paper today and illustrate how the model can be implemented in R.
Let's start with the example data of an incremental claims triangle:
Christofides model assumes the following structure for the incremental paid claims \(P_{ij}\):
\begin{align}
\ln(P_{ij}) & = Y_{ij} = a_i + b_j + \epsilon_{ij}
\end{align}where i and j go from 0 to 3, \(b_0=0\) and \(\epsilon_{ij} \sim N(0, \sigma^2)\). Unlike the basic chain-ladder method, this is a stochastic model that allows me to test my assumptions and calculate various statistics, e.g. standards errors of my predictions.
The paper is available together with a spread sheet model, illustrating the calculations. It is very much based on ideas by Barnett and Zehnwirth, see [2] for a reference. However, doing statistical analysis in a spread sheet programme is often cumbersome. I will go through the first 15 pages of Christofides' paper today and illustrate how the model can be implemented in R.
Let's start with the example data of an incremental claims triangle:
## Page D5.4
tri <- t(matrix(
c(11073, 6427, 1839, 766,
14799, 9357, 2344, NA,
15636, 10523, NA, NA,
16913, NA, NA, NA),
nc=4, dimnames=list(origin=0:3, dev=0:3)))
Christofides model assumes the following structure for the incremental paid claims \(P_{ij}\):
\begin{align}
\ln(P_{ij}) & = Y_{ij} = a_i + b_j + \epsilon_{ij}
\end{align}where i and j go from 0 to 3, \(b_0=0\) and \(\epsilon_{ij} \sim N(0, \sigma^2)\). Unlike the basic chain-ladder method, this is a stochastic model that allows me to test my assumptions and calculate various statistics, e.g. standards errors of my predictions.
Actuarial
,
Barnett Zehnwirth
,
ChainLadder
,
IBNR
,
Insurance
,
linear model
,
log-incremental
,
R
,
reserving
,
Risk
,
Tutorials
Claims reserving in R: ChainLadder 0.1.5-4 released
Last week we released version 0.1.5-4 of the
The chain-ladder method which is a popular method in the insurance industry to forecast future claims payments gave the package its name. However, the
Since we published version 0.1.5-2 in March 2012 additional functionality has been added to the package, see the change log, but in particular the vignette has come a long way.
Many thanks to my co-authors Dan Murphy and Wayne Zhang.
ChainLadder package on CRAN. The R package provides methods which are typically used in insurance claims reserving. If you are new to R or insurance check out my recent talk on Using R in Insurance.The chain-ladder method which is a popular method in the insurance industry to forecast future claims payments gave the package its name. However, the
ChainLadder package has many other reserving methods and models implemented as well, such as the bootstrap model demonstrated below. It is a great starting point to learn more about stochastic reserving.Since we published version 0.1.5-2 in March 2012 additional functionality has been added to the package, see the change log, but in particular the vignette has come a long way.
Many thanks to my co-authors Dan Murphy and Wayne Zhang.
Actuarial
,
ChainLadder
,
Insurance
,
News
,
R
,
reserving
,
stochastic reserving
,
vignette
Using R in Insurance, Presentation at GIRO 2012
Every year the UK’s general insurance actuarial community organises a big conference, which they call GIRO, short for General Insurance Research Organising committee.
This year's conference is in Brussels from 18 - 21 September 2012. Despite the fact that Brussels is actually in Belgium the UK actuaries will travel all the way to enjoy good beer and great talks.
On Wednesday morning I will run a session on Using R in insurance. It would be great to see some of you there.
I prepared the slides with R, RStudio, knitr, pandoc and slidy again. My title page shows a word cloud about the GIRO conference. It uses the wordcloud package and was inspired by Ian Fellows' post on FellStats.
The last slide shows the output of sessionInfo(). I am sure it will become helpful one day, when I have to remind myself how I actually created the slides and which packages and versions I used.
This year's conference is in Brussels from 18 - 21 September 2012. Despite the fact that Brussels is actually in Belgium the UK actuaries will travel all the way to enjoy good beer and great talks.
On Wednesday morning I will run a session on Using R in insurance. It would be great to see some of you there.
I prepared the slides with R, RStudio, knitr, pandoc and slidy again. My title page shows a word cloud about the GIRO conference. It uses the wordcloud package and was inspired by Ian Fellows' post on FellStats.
The last slide shows the output of sessionInfo(). I am sure it will become helpful one day, when I have to remind myself how I actually created the slides and which packages and versions I used.
Actuarial
,
ChainLadder
,
GIRO
,
Insurance
,
knitr
,
pandoc
,
Presentations
,
R
,
slidy
,
Tutorials
