Category Archives: Actuarial science

Amsterdam

I will be in Amsterdam for the end of this week. I will be in the jury of the PhD defense of Julien Tomas, entitled “Quantifying Biometric Life Insurance Risks With Non-Parametric Smoothing Methods” (the thesis will probably be online soon). But before, I will give a talk at the actuarial seminar at UvA. My visit last time was a real pleasure, and it should be the same this time too. I will give a talk this Thursday on “R for actuarial science“. The slides can be downloaded from here.

R for actuarial science

As mentioned in the Appendix of Modern Actuarial Risk Theory, “R (and S) is the ‘lingua franca’ of data analysis and statistical computing, used in academia, climate research, computer science, bioinformatics, pharmaceutical industry, customer analytics, data mining, finance and by some insurers. Apart from being stable, fast, always up-to-date and very versatile, the chief advantage of R is that it is available to everyone free of charge. It has extensive and powerful graphics abilities, and is developing rapidly, being the statistical tool of choice in many academic environments.

R is based on the S statistical programming language developed by Joe Chambers at Bell labs in the 80’s. To be more specific, R is an open-source implementation of the S language, developed by Robert Gentlemn and Ross Ihaka. It is a vector based language, which makes it extremely interesting for actuarial computations. For instance, consider some Life Tables,

> TD[39:52,]       > TV[39:52,]
     Age    Lx         Age    Lx
  39  38 95237          38 97753
  40  39 94997          39 97648
  41  40 94746          40 97534
  42  41 94476          41 97413
  43  42 94182          42 97282
  44  43 93868          43 97138
  45  44 93515          44 96981
  46  45 93133          45 96810
  47  46 92727          46 96622
  48  47 92295          47 96424
  49  48 91833          48 96218
  50  49 91332          49 95995
  51  50 90778          50 95752
  52  51 90171          51 95488

Those (French) Life Tables can be found here

> TD <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",sep=";",header=TRUE)
> TV <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",sep=";",header=TRUE)

From those vectors, it is possible to construct the matrix of death probabilities, https://latex.codecogs.com/gif.latex?\boldsymbol{P}=[\text{%20}_{k}p_x], using for instance

>  Lx <- TD$Lx
>  m <- length(Lx)
>  p <- matrix(0,m,m); d <- p
>  for(i in 1:(m-1)){
+  p[1:(m-i),i] <- Lx[1+(i+1):m]/Lx[i+1]
+  d[1:(m-i),i] <- (Lx[(1+i):(m)]-Lx[(1+i):(m)+1])/Lx[i+1]}
>  diag(d[(m-1):1,]) <- 0
>  diag(p[(m-1):1,]) <- 0
>  q <- 1-p

One can compute easily, e.g., the (curtate) expectation of life defined as

https://latex.codecogs.com/gif.latex?e_x%20=\mathbb{E}(K_x)=\sum_{k=1}^\infty%20k\cdot%20\text{%20}_{k|1}q_x%20=%20\sum_{k=1}^\infty%20\text{%20}_{k}p_x

and one can compute the vector of life expectancy, at various ages https://latex.codecogs.com/gif.latex?\boldsymbol{e}=[e_x], as

> life.exp = function(x){sum(p[1:nrow(p),x])}
> e = Vectorize(life.exp)(1:m)

An actually, any kind of actuarial quantity can be derived from those matrices. The expected present value (or actuarial value) of a temporary life annuity-due is, for instance,

https://latex.codecogs.com/gif.latex?\ddot{a}_{x:\overline{n}|}=\sum_{k=0}^{n-1}%20\nu^k%20\cdot%20{}_{k}p_x%20=\frac{1-A_{x:\overline{n}|}}{1-\nu}

The code to compute those functions is here

> for(j in 1:(m-1)){ adots[,j]<-cumsum(1/(1+i)^(0:(m-1))*c(1,p[1:(m-1),j])) }

or consider the expected present value of a term insurance

https://latex.codecogs.com/gif.latex?%20A^1_{x:\overline{n}|}%20=\sum_{k=0}^{n-1}%20\nu^{k+1}%20\cdot%20\text{%20}_{k|}q_x

with the following code

> for(j in 1:(m-1)){ A[,j]<-cumsum(1/(1+i)^(1:m)*d[,j]) }

Some more details can be found in the first part of the notes of the crash courses of last summer, in Meielisalp. Vector – or matrices – are extremely convenient to work with, when dealing with life contingencies. It is also possible to model prospective mortality. Here, the mortality is not only function of the age https://latex.codecogs.com/gif.latex?x, but also time https://latex.codecogs.com/gif.latex?t,

> t(DTF)[1:10,1:10]
    1899  1900  1901  1902  1903  1904  1905  1906  1907  1908
0  64039 61635 56421 53321 52573 54947 50720 53734 47255 46997
1  12119 11293 10293 10616 10251 10514  9340 10262 10104  9517
2   6983  6091  5853  5734  5673  5494  5028  5232  4477  4094
3   4329  3953  3748  3654  3382  3283  3294  3262  2912  2721
4   3220  3063  2936  2710  2500  2360  2381  2505  2213  2078
5   2284  2149  2172  2020  1932  1770  1788  1782  1789  1751
6   1834  1836  1761  1651  1664  1433  1448  1517  1428  1328
7   1475  1534  1493  1420  1353  1228  1259  1250  1204  1108
8   1353  1358  1255  1229  1251  1169  1132  1134  1083   961
9   1175  1225  1154  1008  1089   981  1027  1025   957   885

Thus, we now have a force of mortality matrix https://latex.codecogs.com/gif.latex?\boldsymbol{\mu}=[\mu_{x,t}], or surface

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/Capture-d%E2%80%99e%CC%81cran-2013-01-10-a%CC%80-14.29.04.png

It is also possible to use R packages to estimate a Lee-Carter model of the mortality rate,

https://latex.codecogs.com/gif.latex?\log%20\mu%20_{x,t}%20=\alpha%20_{x}%20+\beta%20_{x}%20\cdot%20\kappa_{t}%20+\varepsilon%20_{x,t}

> library(demography)
> MUH =matrix(DEATH$Male/EXPOSURE$Male,nL,nC)
> POPH=matrix(EXPOSURE$Male,nL,nC)
> BASEH <- demogdata(data=MUH, pop=POPH, ages=AGE, years=YEAR, type="mortality",
+ label="France", name="Hommes", lambda=1)
> RES=residuals(LCH,"pearson")

One can easily study residuals, for instance as a function of the age,

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/Capture-d%E2%80%99e%CC%81cran-2013-01-10-a%CC%80-14.29.15.png

or a function of the year,

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/Capture-d%E2%80%99e%CC%81cran-2013-01-10-a%CC%80-14.29.22.png

Some more details can be found in the second part of the notes of the crash courses of last summer, in Meielisalp.

R is also interesting because of its huge number of libraries, that can be used for predictive modeling. One can easily use smoothing functions in regression, or regression trees,

> TREE = tree((nbr>0)~ageconducteur,data=sinistres,split="gini",mincut = 1)
> age = data.frame(ageconducteur=18:90)
> y1 = predict(TREE,age)
> reg = glm((nbr>0)~bs(ageconducteur),data=sinistres,family="binomial")
> y = predict(reg,age,type="response")

http://freakonometrics.hypotheses.org/files/2013/01/predictive-gam-tree.png

Some practitioners might be scared because the legend claims that R is not as good as SAS to handle large databases. Actually, a lot of functions can be used to import datasets. The most convenient one is probably

> baseCOUT = read.table("http://freakonometrics.free.fr/baseCOUT.csv",
+  sep=";",header=TRUE,encoding="latin1")
>  tail(baseCOUT,4)
     numeropol  debut_pol    fin_pol freq_paiement langue  type_prof alimentation type_territoire
6512     87291 2002-10-16 2003-01-22       mensuel      A Professeur   Vegetarien          Urbain
6513     87301 2002-10-01 2003-09-30       mensuel      A Technicien   Vegetarien          Urbain
6514     87417 2002-10-24 2003-10-21       mensuel      F Technicien   Vegetalien     Semi-urbain
6515     88128 2003-01-17 2004-01-16       mensuel      F     Avocat   Vegetarien     Semi-urbain
             utilisation presence_alarme marque_voiture sexe exposition age duree_permis age_vehicule i   coutsin
6512 Travail-occasionnel             oui           FORD    M  0.2684932  47           29           28 1 1274.5901
6513              Loisir             oui          HONDA    M  0.9972603  44           24           25 1  278.0745
6514 Travail-occasionnel             non     VOLKSWAGEN    F  0.9917808  23            3           11 1  403.1242
6515              Loisir             non           FIAT    F  0.9972603  23            4           11 1  230.9565

But if the dataset is too large, it is also possible to specify which variables might be interesting, using

> mycols = rep("NULL", 18)
> mycols[c(1,4,5,12,13,14,18)] <- NA
> baseCOUTsubC = read.table("http://freakonometrics.free.fr/baseCOUT.csv",
+  colClasses = mycols,sep=";",header=TRUE,encoding="latin1")
> head(baseCOUTsubC,4)
  numeropol freq_paiement langue sexe exposition age    coutsin
1         6        annuel      A    M  0.9945205  42   279.5839
2        27       mensuel      F    M  0.2438356  51   814.1677
3        27       mensuel      F    M  1.0000000  53   136.8634
4        76       mensuel      F    F  1.0000000  42   608.7267

It is also possible (before running a code on the entire dataset) to import only the first lines of the dataset.

> baseCOUTsubCR = read.table("http://freakonometrics.free.fr/baseCOUT.csv",
+  colClasses = mycols,sep=";",header=TRUE,encoding="latin1",nrows=100)
> tail(baseCOUTsubCR,4)
    numeropol freq_paiement langue sexe exposition age   coutsin
97       1193       mensuel      F    F  0.9972603  55  265.0621
98       1204       mensuel      F    F  0.9972603  38 9547.7267
99       1231       mensuel      F    M  1.0000000  40  442.7267
100      1245        annuel      F    F  0.6767123  48  179.1925

It is also possible to import a zipped file. The file itself has a smaller size, and it can usually be imported faster.

> import.zip = function(file){
+ temp = tempfile()
+ download.file(file,temp);
+ read.table(unz(temp, "baseFREQ.csv"),sep=";",header=TRUE,encoding="latin1")}
> system.time(import.zip("http://freakonometrics.free.fr/baseFREQ.csv.zip"))
trying URL 'http://freakonometrics.free.fr/baseFREQ.csv.zip'
Content type 'application/zip' length 692655 bytes (676 Kb)
opened URL
==================================================
downloaded 676 Kb
   user  system elapsed 
      0.762       0.029       4.578 
> system.time(read.table("http://freakonometrics.free.fr/baseFREQ.csv", 
+ sep=";",header=TRUE,encoding="latin1"))
   user  system elapsed 
      0.591       0.072       9.277

Finally, note that it is possible to import any kind of dataset, not only a text file. Even a Microsoft Excel folder. On a Windows computer, one can use SQL queries

> sheet = "c:\\Documents and Settings\\user\\excelsheet.xls"
> connection = odbcConnectExcel(sheet)
> spreadsheet = sqlTables(connection)
> query = paste("SELECT * FROM",spreadsheet$TABLE_NAME[1],sep=" ")
> result = sqlQuery(connection,query)

Then, once the dataset is imported, several functions can be used,

> cost = aggregate(coutsin~ AgeSex,mean, data=baseCOUT)
> frequency = merge(aggregate(nbsin~ AgeSex,sum, data=baseFREQ),
+ aggregate(exposition~ AgeSex,sum, data=baseFREQ))
> frequency$freq = frequency$nbsin/frequency$exposition
> base.freq.cost = merge(frequency, cost)

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/cost-freq-qc.png

Finally, R is interesting for its graphical interface. “If you can picture it in your head, chances are good that you can make it work in R. R makes it easy to read data, generate lines and points, and place them where you want them. Its very flexible and super quick. When youve only got two or three hours until deadline, R can be brilliant” as said Amanda Cox, a graphics editor at the New York Times. “R is particularly valuable in deadline situations when data is scant and time is precious.”.
Several cases were considered on the blog http ://chartsnthings.tumblr.com/…. First, we start with a simple graph, here State Government control in the US

http://freakonometrics.hypotheses.org/files/2013/01/nyt-chartsnthings-1.png

Then try to find a nice visual representation, e.g.

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/nyt-chartsnehings-2.png

And finally, you can just print it in your favorite newspaper,

http://freakonometrics.hypotheses.org/files/2013/01/nyt-chartsnthings-3.jpg

And you can get any kind of graphs,

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/nyt-6.png

And not only about politics,

http://freakonometrics.hypotheses.org/files/2013/01/nyt-7-b.jpg Graphs are important. “Its not just about producing graphics for publication. Its about playing around and making a bunch of graphics that help you explore your data. This kind of graphical analysis is a really useful way to help you understand what you’re dealing with, because if you cant see it, you cant really understand it. But when you start graphing it out, you can really see what you’ve got” as said Peter Aldhous, San Francisco bureau chief of New Scientist magazine. Even for actuaries. “The commercial insurance underwriting process was rigorous but also quite subjective and based on intuition. R enables us to communicate our analytic results in appealing and innovative ways to non-technical audiences through rapid development lifecycles. R helps us show our clients how they can improve their processes and effectiveness by enabling our consultants to conduct analyses efficiently”, as explained by John Lucker, team of advanced analytics professionals at Deloitte Consulting Principal, in http://blog.revolutionanalytics.com/r-is-hot/. See also Andrew Gelman’s view, on graphs, http://www.stat.columbia.edu/…

So yes, actuaries might be interested to use R for actuarial communication, as mentioned in http ://www.londonr.org/…

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/mango-R-4.png

The Actuarial Toolkit (see http ://www.actuaries.org.uk/…) stresses the interest of R, “The power of the language R lies with its functions for statistical modelling, data analysis and graphics ; its ability to read and write data from various data sources; as well as the opportunity to embed R in excel or other languages like VBA. In the way SAS is good for data manipulations, R is superior for modelling and graphical output“.

From 2011, Asia Capital Reinsurance Group (ACR) uses R to Solve Big Data Challenges (see http ://www.reuters.com/…). And Lloyd’s uses motion charts created with R to provide analysis to investors (as discussed on http ://blog.revolutionanalytics.com/…)

A lot of information can be found on http ://jeffreybreen.wordpress.com/…

http://freakonometrics.hypotheses.org/files/2013/01/6a010534b1db25970b01538fea1796970b-800wi.png

Markus Gesmann mentioned on his blog a lot of interesting graphs used for actuarial reporting, http ://lamages.blogspot.ca/…

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/Capture-d%E2%80%99e%CC%81cran-2013-01-10-a%CC%80-15.37.33.png

Further, R is free. Which can be compared with SAS, $6,000 per PC, or $28,000 per processor on a server (as mentioned on http ://en.wikipedia.org/…)

It is also becoming more and more popular, as a programming language. As mentioned on this month Transparent Language Popularity (see http ://lang-index.sourceforge.net/), R is ranked 12. Far away after C or Java, but before Matlab (22) or SAS (27). On StackOverFlow (see http ://stackoverflow.com/) is also far being C++ (399,232 occurrences) or Java (348,418), but with 21,818 occurrences, it appears before Matlab (14,580) and SAS (899). As mentioned on http ://r4stats.com/articles/popularity/ R is becoming more and more popular, on listserv discussion traffic

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/fig_1_listserv.png

It is clearly the most popular software in data analysis, as mentioned by the Rexer Analytics survey, in 2009

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/fig_3_rexersurvey.png

What about actuaries ? In a survey (see http ://palisade.com/…), R was not extremely popular.

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/mango-R-1.png

If we consider only statistical softwares, SAS is still far ahead, among UK and CAS actuaries

http://freakonometrics.hypotheses.org/wp-content/blogs.dir/253/files/2013/01/mango-R-2.png

But, as mentioned by Mike King, Quantitative Analyst, Bank of America, “I cant think of any programming language that has such an incredible community of users. If you have a question, you can get it answered quickly by leaders in the field. That means very little downtime.” This was also mentioned by Glenn Meyers, in the Actuarial Review “The most powerful reason for using R is the community” (in http ://nytimes.com/…). For instance, http ://r-bloggers.com/ has contributions from more than 425 R users.

As said by Bo Cowgill, from Google “The best thing about R is that it was developed by statisticians. The worst thing about R is that it was developed by statisticians.

Segmentation en tarification, compléments

Dans le premier cours d’actuariat IARD, nous avons vu l’importance de la ségmentation, et son implication sur le calcul des primes (passer d’une espérance mathématique à une espérance conditionnelle). Pour aller un peu plus loin, quelques compléments,

pour une lecture plus économique de la problématique de la segmentation en assurance

ou pour une lecture plus légale

Sinon, plusieurs articles de vulgarisation peuvent être lu sur internet,

La première démo aura lieu lundi, en salle informatique. Karim sera une introduction au langage R, à la manipulation des variables (qualitatives et quantitatives). Je mettrais en ligne les transparents en fin de semaine, et les codes seront mis en ligne dans le courant de la semaine prochaine.

Comme annoncé hier, il n’y aura pas cours mercredi prochain. Le mercredi suivant, nous verrons la modélisation des variables indicatrices, i.e. la régression logistique, et les arbres de régression. On supposera que le modèle linéaire aura été vu, je mets un lien vers les transparents du cours ACT6420 de la session passée, notes de cours transparents1 et transparents2. Il est aussi possible de relire Frees (2010), chapitres 3, 4, 5 et 6.

Pour commencer à pratiquer la régression logistique, on utilisera la petite fonction suivante

logit = function(formula, lien="logit", data=NULL) {
glm(formula,family=binomial(link=lien),data)
}

Sinon, la Casualty Actuarial Society a mis en ligne plusieurs documents en ligne sur les arbres de régression (qui sont peu abordés dans les livres mentionnés auparavant),

pour une comparaison de toutes les méthodes

Les transparents seront mis en ligne en fin de semaine prochaine. A suivre donc…

Generating a non-homogeneous Poisson process

Consider a Poisson process gif.latex (54×20), with non-homogeneous intensity . Here, we consider a deterministic function, not a stochastic intensity. Define the cumulated intensity

in the sense that the number of events that occurred between time gif.latex (8×13) and gif.latex (6×12) is a random variable that is Poisson distributed with parameter  .

For example, consider here a cyclical Poisson process, with intensity

   lambda=function(x) 100*(sin(x*pi)+1)

To compute the cumulated intensity, consider a very general function

   Lambda=function(t) integrate(f=lambda,lower=0,upper=t)$value

The idea is to generate a Poisson process on a finite interval .

The first code is based on a proposition from Çinlar (1975),

  1. start with https://latex.codecogs.com/gif.latex?s=0
  2. generate gif.latex (96×19)
  3. set gif.latex (112×19)
  4. set gif.latex (6×12) denote gif.latex (124×19)
  5. deliver
  6. go to step 2.

In order to get the infinimum of gif.latex (12×13), consider a code as

   v=seq(0,Tmax,length=1000)
   t=min(v[which(Vectorize(Lambda)(v)>=s)])

(it might not be very efficient…. but it should work). Here, the code to generate that Poisson process is

   s=0; v=seq(0,Tmax,length=1000)
   X=numeric(0)
   while(X[length(X)]<=Tmax){
     u=runif(1)
     s=s-log(u)
     t=min(v[which(Vectorize(Lambda)(v)>=s)])
     X=c(X,t)
   }

Here, we get the following histogram,

   hist(X,breaks=seq(0,max(X)+1,by=.1),col="yellow")
   u=seq(0,max(X),by=.02)
   lines(u,lambda(u)/10,lwd=2,col="red")

Consider now another strategy. The idea is to use the conditional distribution before the next event, given that one occurred at time ,

  1. start with
  2. generate gif.latex (51×16)
  3. set gif.latex (74×14)
  4. deliver
  5. go to step 2.

Here the algorithm is simple. For the computational side, at each step, we have to compute and then http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=F_t%5E%7B-1%7D. To do so, since is increasing with values in , we can use a dichotomic algorithm,

   Ft=function(x) 1-exp(-Lambda(t+x)+Lambda(t))
   Ftinv=function(u){
     a=0
     b=Tmax
     for(j in 1:20){
       if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}
       if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
       a=binf
       b=bsup
     }
   return((a+b)/2)
   }

Here the code is the following

   t=0; X=t
   while(X[length(X)]<=Tmax){
     Ft=function(x) 1-exp(-Lambda(t+x)+Lambda(t))
     Ftinv=function(u){
      a=0
      b=Tmax
      for(j in 1:20){
        if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}
        if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
        a=binf
        b=bsup
      }
      return((a+b)/2)
     }
     x=Ftinv(runif(1))
     t=t+x
     X=c(X,t)
   }

The third code is based on a classical algorithm to generate an homogeneous Poisson process on a finite interval: first, we generate the number of events, then, we draw uniform variates, and we sort them. Here, the strategy is closed, except that is won’t be uniform any longer.

  1. generate the number of events on the time interval gif.latex (101×19)
  2. generate independently gif.latex (114×17) where 
  3. set gif.latex (60×15) i.e. the ordered values  gif.latex (136×16)
  4. deliver http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=t_i‘s

This algorithm is extremely simple, and also very fast. This is one function to inverse, and it is not in the loop,

   n=rpois(1,Lambda(Tmax))
   Ft=function(x) Lambda(x)/Lambda(Tmax)
   Ftinv=function(u){
     a=0
     b=Tmax
     for(j in 1:20){
       if(Ft((a+b)/2)<=u){binf=(a+b)/2;bsup=b}
       if(Ft((a+b)/2)>=u){bsup=(a+b)/2;binf=a}
       a=binf
       b=bsup
     }
     return((a+b)/2)
     }
   X0=rep(NA,n)
   for(i in 1:n){
     X0[i]=Ftinv(runif(1))
    }
   X=sort(X0)

Here is the associated histogram,

An alternative is based on a rejection technique. Actually, it was the algorithm mentioned a few years ago on this blog (well, the previous one). Here, we need an upper bound for the intensity, so that computations might be much faster. Here, consider

  1. start with
  2. generate gif.latex (96×19)
  3. set gif.latex (137×19)
  4. generate gif.latex (95×19) (independent of http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=u)
  5. if gif.latex (90×19) then deliver http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=t
  6. go to step 2.

Here, consider a constant upper bound,

   lambdau=function(t) 200
   Lambdau=function(t) lambdau(t)*t

The code to generate a Poisson process is

   t=0
   X=numeric(0)
   while(X[length(X)]<=Tmax){
     u=runif(1)
     t=t-log(u)/lambdau
     if(runif(1)<=lambda(t)/lambdau) X=c(X,t)
  }

The histogram is here

Finally, the last one is also based on a rejection technique, mixed with the second one. I.e. define

gif.latex (433×20)

The good thing is that this function can easily be inverted

gif.latex (215×21)

  1. start (as usual) with
  2. generate gif.latex (63×19)
  3. set gif.latex (74×14)
  4. generate gif.latex (96×19)
  5. if gif.latex (124×19) then deliver http://www.forkosh.com/cgi-bin/mathtex.cgi?formdata=t
  6. goto step 2.

Here, the algorithm is simply

   t=0
   while(X[length(X)]<=Tmax){
     Ftinvu=function(u) -log(1-x)/lambdau
     x=Ftinvu(runif(1))
     t=t+x
     if(runif(1)<=lambda(t+x)/lambdau(t+x)) X=c(X,t)
   }

Obviously those five codes work, the first one being much slower than the other three. But it might be because my strategy to seek the infimum is not great. And the latter worked well since there were not much rejection, I guess it can be worst…

All those algorithms were mentioned in a nice survey written by Raghu Pasupathy and can be downloaded from http://web.ics.purdue.edu/~pasupath/…. In the paper, non-homogeneous spatial Poisson processes are also mentioned…

 

Actuariat IARD

Cet hiver (même si la nouvelle ne sera officielle qu’à la rentrée), je devrais donner le cours ACT2040, actuariat IARD. Le plan de cours sera bientôt en ligne, mais je peux déjà dire que le cours sera basé sur le Tome 2 du livre écrit avec Michel Denuit il y a quelques années, mathématiques de l’assurance non-vie. Le cours est une suite du cours ACT6420 méthodes de prévisions, donné cet automne (qui est un prérequis indiqué sur le site du registrariat http://websysinfo.uqam.ca/…): je partirais donc du fait que le modèle linéaire de régression est connu (et compris) et que tout le monde sait utiliser R, et lire des sorties de régression. Mais les premières démonstrations reviendront sur l’utilisation de R, et sur l’analyse de la variance, que l’on n’a pas vraiment eu le temps d’aborder dans le cours de régression. Pour des références sur R, je conseille

  • “R pour les débutants” d’Emmanuel Paradis, (PDF)
  • “Introduction à la programmation en S” par Vincent Goulet, (PDF)

pour les documents en français, ou pour des documents plus complets, mais en anglais

  • “R for Beginners” d’Emmanuel Paradis (PDF),
  • “An Introduction to R” par Longhow Lam (PDF)
  • “The R language — a short companion” par Marc Vandemeulebroecke (PDF),
  • “The R Guide” par Jason Owen (PDF),
  • “Econometrics in R” par Grant Farnsworth (PDF) pour aller plus loin sur les régressions,
  • “Practical Regression and Anova using R” by Julian Faraway (PDF) sur le meme sujet
  • “Statistics with R and S-Plus” d’Hugo Quené (PDF)
  • “Statistical Computing and Graphics Course Notes” par Frank Harrell, (PDF).
  • “Using R for Data Analysis and Graphics – Introduction, Examples and Commentary” par John Maindonald (PDF).

Sinon, les transparents du premier cours sont en ligne ici, et le plan de cours est

et je mettrais bientôt en ligne des liens vers des bases de données que l’on utilisera tout au long du cours, ou en démonstration.

Pour les références, je citerais deux livres sur lesquels je m’appuierai beaucoup car je les connais presque par cœur. Ils sont disponibles à la Coopuqam

Jeux et assurance

Félicitations à Christophe Dutang, qui a obtenu il y a quelques heures le prix Scor de la meilleure thèse de doctorat en actuariat. Christophe avait soutenu sa thèse à Lyon sur “Étude des marchés d’assurance non-vie à l’aide d’équilibres de Nash et de modèles de risques avec dépendance” (les transparents sont en ligne sur sa page, et la thèse est en ligne sur http://tel.archives-ouvertes.fr/…)

La cérémonie avait lieu mercredi soir, et je ne pouvais pas y être car je fêtais mon anniversaire en famille (et accessoirement, je suis à Montréal pour finir la session). En fait, je n’ai jamais eu l’occasion de me rendre à cette cérémonie annuelle (même lorsque le prix m’avait été attribué, voilà quelques années, car j’étais alors à Valparaiso). Je pourrais aussi souligner que si je trouve que le prix est une très bonne idée, le lieu n’est pas idéal.

Histoire de jouer un peu les vieux cons, Le Cercle de l’Union Interallié est un club très select, très sexiste (« les femmes n’ont pas accès aux instances dirigeantes du Cercle » comme le précise http://fr.wikipedia.org/…), et en plus la cravate est de rigueur pour les hommes (sinon, on ne peut pas entrer). Si je sors mon costume à l’occasion (certains collègues ont eu l’occasion de se moquer), je dois avouer que mes cravates sont depuis fort longtemps dans la malle à déguisement des enfants. Et elles ne sont pas prêt d’en sortir !

En attendant, félicitations à Christophe qui méritait le prix… Et félicitations à Aymric Kamega qui a obtenu le même soir une mention spéciale !

Actuariat en Afrique subsaharienne francophone

Depuis que le blog (ou les précédents) existe, je sais (par les messages que je reçois par courriel) qu’il est beaucoup lu en Afrique (francophone). J’ai reçu avant hier un livre publié par Aymric Kamega et Frédéric Planchet, tiré de la thèse de doctorat d’Aymric (en ligne sur http://halshs.archives-ouvertes.fr/) défendue en décembre dernier. J’étais alors rapporteur extérieur, et j’avais souligné l’intérêt des travaux d’Aymric dans le cadre de la volonté de la CIMA (autorité de contrôle régionale des marchés d’assurance pour l’Afrique subsaharienne francophone) de fournir aux assureurs de la région des outils adaptés. En particulier des tables de mortalité d’expérience, propres à la région. En effet, (comme le rappelle Aymric), , suite aux états généraux (de l’assurance vie)  en 2007, le principe de la construction de tables de mortalité d’expérience a été adopté en remplacement des tables de mortalité de la population générale française entre 1960 et 1964 (tables dites PM 60-64 et PF 60-64) jusqu’alors imposées. Aymric avait travaillé sur ce sujet, et a soutenu une thèse de doctorat qui avait fait l’unanimité. J’avais pris beaucoup de plaisir à lire la thèse, j’en pense que j’en aurais à lire le livre (disponible sur le site de l’éditeur http://seddita.com/).

Compound Poisson and vectorized computations

Yesterday, I was asked how to write a code to generate a compound Poisson variables, i.e. a series of random variables  where  is a counting random variable (here Poisson disributed) and where the ‘s are i.i.d (and independent of ), with the convention  when . I came up with the following algorithm, but I was wondering if it was possible to get a better one…

>  rcpd=function(n,rN,rX){
+  N=rN(n)
+  X=rX(sum(N))
+  I=as.factor(rep(1:n,N))
+  S=tapply(X,I,sum)
+  V=as.numeric(S[as.character(1:n)])
+  V[is.na(V)]=0
+  return(V)}

Here, consider – to illustrate – the case where  and ,

>  rN.P=function(n) rpois(n,5)
>  rX.E=function(n) rexp(n,2)

We can generate a sample

>  S=rcpd(1000,rN=rN.P,rX=rX.E)

and check (using simulation) than 

> mean(S)
[1] 2.547033
> mean(rN.P(1000))*mean(rX.E(1000))
[1] 2.548309

and that 

> var(S)
[1] 2.60393
> mean(rN.P(1000))*var(rX.E(1000))+
+ mean(rX.E(1000))^2*var(rN.P(1000))
[1] 2.621376

If anyone might think of a faster algorithm, I’d be glad to hear about it…

Longevity and mortality dynamics with R

Following the previous post on life contingencies and actuarial models in life insurance, I upload additional material for the short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The second part of the talk (on Actuarial models with R) will be dedicated to longevity and mortality. A complete set of slides can be downloaded from the blog, but again, only some part will be presented.

As mentioned earlier, the codes are from a book on actuarial science in R, written with Christophe Dutang (so far in French) that should appear, some day… The code used in the slides above can be downloaded from here, and datasets are the following,

> DEATH <- read.table(
+ "http://freakonometrics.free.fr/Deces-France.txt",
+ header=TRUE)
> EXPO  <- read.table(
+ "http://freakonometrics.free.fr/Exposures-France.txt",
+ header=TRUE,skip=2)

For additional resources, I will use Rob Hyndman‘s package on demography, Heather Turner and David Firth’s package on generalized nonlinear models (e.g. the slides of the short course Heather gave in Rennes at the UseR! conference in 2009), as well as functions developed by JPMorgan’s LifeMetrics (functions are  fully documented in the LifeMetrics Technical Document). All those functions can be obtained using

> library(demography)
> library(gnm)
> source("http://freakonometrics.free.fr/fitModels.R")

Life contingencies with R

I will be giving in less than four weeks a short course at the 6th R/Rmetrics Meielisalp Workshop & Summer School on Computational Finance and Financial Engineering organized by ETH Zürich, https://www.rmetrics.org/. The talk will be on Actuarial models with R, and first part will be dedicated to life insurance. A complete set of slides can be downloaded from the blog, but in the talk, only some part will be presented.

The codes are from a book on actuarial science in R, written with Christophe Dutang (so far in French) that should appear, some day… The code used in the slides can be downloaded from here, and datasets are the following,

> TD <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TD8890.csv",sep=";",header=TRUE)
> TV <- read.table(
+ "https://perso.univ-rennes1.fr/arthur.charpentier/TV8890.csv",sep=";",header=TRUE)

For additional resources, I recommend Emiliano’s website, http://www.math.uconn.edu/, with great lectures on life insurance mathematics, and the (new) lifecontinfencies vignette on http://cran.r-project.org/,

> library(lifecontingencies)

Ruin probability and infinite time

A couple of weeks ago, I had a discussion with a practitioner, working in some financial company, about ruin, and infinite time. And it reminded me a weird result. Well, not a weird result, but a result I found disturbing, at first, when I was a student (that I rediscovered with the eyes of someone dealing with computational issues, seeing here a difficult theoretical question). Consider a simple ruin problem. A player has wealth . Then he flips a coin: tails he has a gain of 1, heads he experiences a loss of 1. At time , his wealth is where  is associated to the th coin:  is equal to 1 with probability (tails), and -1 with probability  (heads). It is also possible to write

where  can be interpreted as the net gain of the player. In order to get a good understanding of results that can be obtained. Assume  to be given. Let denote the number of heads and  the number of tails. Then , while . Let  denote the number of paths to go from point A (wealth  at time ) to point B (wealth  at time ). Note that this is a Markovian problem, that can be modeled using Markov chains

But here, we will focus on combinatorial results. Hence,

In order to derive probabilities to reach , let  denote the number of paths going from  to . And let denote the number of paths going from  to  that do reach  at some point between  and . Using a simple reflexion property, then if  and  are positive,

Based on those reflexions, two results can be derived (focusing on probability, instead of counting paths). First, we can obtain that

(given that n and x have the same parity). The second result we can obtain is that

Based on those two expressions, if  denotes the first time  become null, given ,

then

This can be computed easily,

> x=10
> p=.55
> ProbN=function(n){
+ pb=0
+ if(abs(n-x) %% 2 == 0)
+ pb=x/n*choose(n,(n+x)/2)*(1-p)^((n+x)/2)*(p)^((n-x)/2)
+ return(pb)}
> plot(Vectorize(ProbN)(1:1000),type="s")

That looks nice… But if we look closer, we can wonder what

would be ? Since we have the distribution of a probabilty measure, we might expect one. But here

> sum(Vectorize(ProbN)(1:1000))
[1] 0.134385

And this is not due to calculation mistakes that we do not get 1 here. Actually, we should write

which might be interpreted as the probability of ruin, starting from , that we denote  from now on. The term on the left can be approximated using monte-carlo simulations

> p=.55
> x=10
> m=1000
> simul=10000
> S=sample(c(-1,1),size=m*simul,replace=TRUE,prob=c(1-p,p))
> MS=matrix(S,simul,m)
> for(k in 2:m) MS[,k]=MS[,k]+MS[,k-1]
> T0=function(vm) which(vm<=(-x))[1]
> MTmin=apply(MS,1,T0)
> mean(is.na(MTmin)==FALSE)
[1] 0.1328

To check the validity of the relationship above, a simple (theoretical) recursive formula can be derived for the term on the right (ruin probability), namely

with a boundary conditions , and . Then is comes that

Note that it might be tricky to check using monte carlo simulation… since we cannot have an infinite number of runs. And we’re dealing precisely with things that do occur when time is infinite. Actually, we can still check convergence, considering an upper limit  for the number of runs, and then letting  go to infinity. Note that an explicit formula can then be derived (using additional border condition )

Using the following code, it is possible to calculate ruin probability, in order to estimate .

> MSmin=apply(MS,1,min)
> mean(MSmin<=(-x))
[1] 0.1328
> (((1-p)/p)^x-((1-p)/p)^m)/(1-((1-p)/p)^m)
[1] 0.1344306

The following graph shows the evolution of ruin probability as a function of initial wealth (with monte carlo simulation, with a fixed horizon – including a confidence interval – versus the analytical expression)

Hence, with stopping times, one should remember that

and that those two terms can be approximated simply using simulations or standard approximations.

La maudite constante des modèles ARIMA

Dans les modèles ARIMA, autant le dire tout de suite, les constantes c’est pénible ! Pour comprendre un peu mieux ce qui se passe, considérons ici un modèle ARMA avec constante. Ou pour commencer, juste un AR(1)

Supposons que la série  soit stationnaire, de moyenne . Alors en prenant l’espérance de part et d’autre dans l’expression précédente, , i.e. , ou . Plus généralement, avec un ARMA, , en prenant la encore l’espérance des deux cotés, on en déduit .

Simulons un processus AR(1),

Pour simuler des processus ARIMA, on pourrait utiliser la commande

> X=arima.sim(list(order=c(1,0,0),ar=1/3),n=1000)+2
> mean(X)
[1] 1.931767

mais comme le but est de comprendre ce qui se passe, autant faire les choses calmement,

> X=rep(NA,1010)
> X[1]=0
> for(t in 2:1010){X[t]=4/3+X[t-1]/3+rnorm(1)}
> X=X[-(1:10)]
> mean(X)
[1] 2.03397

I.e. ici le processus est de moyenne qui vaut 2.


Regardons maintenant ce que donnerait l’estimation du processus AR(1),

> arima(X, order = c(1, 0, 0))

Call:
arima(x = X, order = c(1, 0, 0))

Coefficients:
ar1  intercept
0.3738     2.0334
s.e.  0.0294     0.0487

sigma^2 estimated as 0.9318:  log likelihood = -1383.68

De manière un peu surprenante, le coefficient appelé intercept n’est pas la constante dans le modèle AR(1), mais la moyenne du processus. Autrement dit, R n’ajuste pas un processus

comme le laisserait penser l’intuition, mais un processus

Ces deux formes sont bien entendu équivalentes. Mais les coefficients estimés ne sont pas tout à fait ce que l’on attendait…
Plaçons nous maintenant dans le cas d’un processus non-stationnaire. L’extension naturelle serait de considérer un processus ARIMA(1,1,0), ou bien un processus tel que la différence soit un processus AR(1). Un processus ARIMA(1,1,0) avec une constante s’écrirait

en utilisant l’opérateur retard. Ceci fait penser à un modèle avec une tendance linéaire. Posons  (afin de se débarrasser de cette tendance). Alors

i.e.

ou encore

On note que  ou encore  suit alors un processus ARIMA(1,1,0) sans constante cette fois.
Afin de visualiser ce que donnerait l’inférence, simulons le processus suivant, qui est un processus AR(1) avec constante que l’on intègre:  avec

Peut-on retrouver les différents paramètres du processus avec R ?
Commençons (là encore) par simuler un tel processus,

> U=rep(NA,1010)
> U[1]=0
> for(t in 2:1010){U[t]=4/3+U[t-1]/3+rnorm(1)}
> U=U[-(1:10)]
> X=cumsum(U)

Sous R, on obtient l’estimation suivant si l’on tente de calibrer un modèle ARIMA(1,1,0)

> arima(X, order = c(1, 1, 0))

Call:
arima(x = X, order = c(1, 1, 0))

Coefficients:
ar1
0.8616
s.e.  0.0160

sigma^2 estimated as1.343:loglikelihood = -1565.63Ça ne convient pas du tout…. On peut tenter un processus AR(1) (avec constante) sur la série différenciée…

> arima(diff(X), order = c(1, 0, 0))

Call:
arima(x = diff(X), order = c(1, 0, 0))

Coefficients:
ar1  intercept
0.3564     2.0200
s.e.  0.0295     0.0486

sigma^2 estimated as 0.9782:  log likelihood = -1406.6

On progresse, sauf que comme auparavant, le terme qui est donne n’est pas la constante dans le modèle ARIMA, mais la moyenne du processus différencié. Mais cette fois, on a un interprétation, c’est que la constante est la pente de la tendance ! Si on estime la pente associée a , on récupère la même valeur:

> arima(X, order = c(1, 1, 0), xreg=1:length(X))

Call:
arima(x = X, order = c(1, 1, 0), xreg = 1:length(X))

Coefficients:
ar1  1:length(X)
0.3566       2.0519
s.e.  0.0296       0.0487

sigma^2 estimated as 0.9787:  log likelihood = -1406.82

Sur la figure ci-dessous, on retrouve le fait qu’en enlevant la tendance linéaire à  donneune série intégrée, sans constante (qui fait penser à une marche aléatoire).

Autrement dit

  • dans le cas d’une série stationnaire, la constante estimée n’est pas du tout la constante, mais la moyenne de la série temporelle
  • dans le cas d’une série non-stationnaire, la constante estimée dans la série différenciée a du sens, au sens ou il s’agit de la pente de la tendance (linéaire) du processus

Avant de conclure, une petite remarque. Quid de la prévision ? Si on commencer par une révision à l’aide du premier processus ARIMA, on obtient une prédiction pour avec un énorme intervalle de confiance (sans aucun bon sens)

> ARIMA1=arima(X, order = c(1, 1, 0))
> ARIMA2=arima(X, order = c(1, 1, 0), xreg=1:length(X))
> Xp1=predict(ARIMA1,20)
> Xp2=predict(ARIMA2,20,newxreg=
+ (length(X)+1):(length(X)+20))
> plot(960:1000,X[960:1000],xlim=c(960,1020),type="l")
> polygon(c(1001:1020,rev(1001:1020)),
+ c(Xp1$pred+2*Xp1$se,rev(Xp1$pred-2*Xp1$se)),
+ col=CL[3],border=NA)
> lines(1001:1020,Xp1$pred,col="red",lwd=2)

intervalle qui se visualise sur le graphique ci-dessous,

Si on regarde pour l’autre modèle,

> lines(1001:1020,Xp2$pred,col="blue",lwd=2)

Tous ceux qui auront reconnu ici des problèmes qui se posent dans la modélisation de la série http://freakonometrics.blog.free.fr/public/maths/viekt.png dans le modèle de Lee & Carter (1992) auront compris qu’il existe une vraie application a tout ce que je viens de raconter (je pense au particulier au commentaire poste par @ClaudeT en début de semaine) car la série des http://freakonometrics.blog.free.fr/public/maths/viekt.png ressemble très fortement à ce genre de série, non-stationnaire avec une tendance linéaire. Et qu’un modèle mal spécifié laisse à penser que l’incertitude est beaucoup beaucoup plus grande que ce qu’elle est vraiment (en plus de donner une prédiction surprenante à long terme !). Donc avant de faire tourner rapidement les codes, il vaut mieux regarder attentivement ce qu’ils font réellement…