Financial model complexity

Today, Olivier Scaillet gave a great talk on fast recursive projections. The idea was great, and the talk was amazing. A great plenary session talk actually. And after lunch, while we were having a coffee, we started to discuss about financial model complexity, mentioning that sometimes, traders and quants are lost, and it might be good to spend more time on basics than on very advanced stuff (which was, in fact, the general idea of my own talk, yesterday). And Olivier recalled that story, on how a rookie excel error led JPMorgan to misreport its VaR for years, published on the blog http://zerohedge.com/…. The short story is that the JPM’s reported VaR did rise by some 93% year over a year, from 2011 to 2012 (while it was decreasing for all competitors). The reason is explained in the very last page of its JPM task force report

… a decision was made to stop using the Basel II.5 model and not to rely on it for purposes of reporting CIO VaR in the Firm’s first-quarter Form 10-Q. Following that decision, further errors were discovered in the Basel II.5 model, including, most significantly, an operational error in the calculation of the relative changes in hazard rates and correlation estimates. Specifically, after subtracting the old rate from the new rate, the spreadsheet divided by their sum instead of their average, as the modeler had intended. This error likely had the effect of muting volatility by a  factor of two and of lowering the VaR…. It also remains unclear when this error was introduced in the calculation.

(Tyler Durden did highlight some parts, and I keep it like that). Let’s admit it: we did have fun about practitioners (actually, there was also a quant sitting with the two of us). But on the other hand, it is a bit scary, to see that we spend so much time to implement complex algorithms, to faster computations, and finally, we end up with a mistake in a spreadsheet….

Advanced methods in trees

I will give a talk tomorrow morning at the Mathematical  Finance Days, organized in HEC Montréal Monday and Tuesday, on Advanced methods in trees with (as mentioned in the subtitle of the first slide) a some thoughts on teaching mathematical finance. It is mainly a survey on advanced tools, based on the idea expressed in Price (1996),

The paper that showed that European option pricing could be put on a rational mathematical basis was Black and Scholes published in 1973. It was so revolutionary that the authors had to submit it to a number of journals before it was accepted. Although there are now numerous approaches to the result, they mostly require specialized methods, including Ito calculus and partial differential  equations, and perhaps Girsanov theory and Feynman-Kac methods. But it is the binomial method due initially to Sharpe and substantially extended by Cox, Ross, and Rubinstein that made the theory of option pricing accessible to everyone with limited mathematical background.  Even though it requires only routine algebraic manipulations, the method is still able to elucidate many of the ideas behind the full theory. Furthermore, all the surprising results mentioned in the opening can be located in this approach. For these reasons it is usually the first method presented in text books and finance courses; we shall follow this trend and step through it. The binomial method is, however, much more than a pedagogical breakthrough, since it allows for the development of numerical approximation methods for a wide range of options for which there are no  known analytic solutions.

Some recent results, obtained in work in progress with colleagues in combinatorial analysis will also be mentioned at the end of the talk (slides can be downloaded in a pdf format, with animations)

I will also be chairing the Numerical Methods session.

Loi de Stigler, langue et publication

Lundi, je dois faire un exposé aux Journées de Finance Mathématiques de Montréal, qui sont organisées comme tous les ans à HEC. Avec les collègues, on a – à nouveau – relancé, entre nous, autour d’un café, le débat du “pourquoi on ne présente pas en français ?” cher aux québécois (le site est d’ailleurs à la fois en français, et en anglais). Mais comme toujours, on râle, on râle, puis on finit par faire des exposés en anglais (transparents et présentation orale). Pourtant, même si j’ai l’impression de pouvoir faire passer les grandes lignes de mon message, je vois bien que je n’arriverais pas à rentrer dans des subtilités, comme je pourrais le faire en français… Mais le monde de la recherche est ainsi fait, et les personnes dont l’anglais est la langue maternelle ont un net avantage. Car si on souhaite diffuser nos travaux, on n’a pas vraiment le choix, il faut passer par l’anglais.

Bien entendu, il y a des contre-exemples. Le plus connu (dans les disciplines que je connais le mieux) est Maurice Allais qui publiait en 1953 en français, dans Econometrica, le comportement de l’homme rationnel devant le risque: critique des postulats et axiomes de l’école Américaine. Et le papier a été abondamment cité. Mais le plus souvent, ne pas publier en anglais, c’est condamner une publication à ne pas avoir l’écho qu’elle pourrait mériter. Par exemple, j’ai découvert (très) récemment qu’en 1940, Bruno de Finetti publiait il problema dei « pieni », dans le Giornale dell’Istituto Italiano degli Attuari, avec des choses qui ressemblent à la section suivante

Tous ceux qui ont vu passer des choses sur les mesures de risques reconnaîtront une axiomatique connue. Mais le papier est encore plus fabuleux, et j’encourage toutes les personnes qui s’intéressent aux risques à aller le lire. Car il a malheureusement été peu lu, et (très) peu cité. En fait, il a fallu attendre que le Journal of Investment Management le traduise, sous le titre the Problem of Full-Risk Insurances, en 2006, avec un préambule intitulé de Finetti scoops Markowitz écrit justement par Harry Markowitz, pour que le monde académique (re)découvre ce papier. Car effectivement, avec plus de 10 ans d’avance, Bruno de Finetti invente la théorie du portefeuille, que l’on attribue toujours à Harry  Markowitz. C’est toutefois Mark Rubinstein qui a (re)découvert la portée du papier de de Finetti, comme il le raconte dans Bruno de Finetti and Mean-Variance Portfolio Selection.

it has recently came to the attention of economists in the English speaking world that among de Finetti’s papers is a treasurer trove of results in economics and finance written well before the works of the scholars that are traditionally credited with these ideas. de Finetti’s 1940 paper anticipating much of mean variance portfolio’s theory later developed by Harry Markowitz

(et qui a obtenu un prix Nobel pour ce travail). Cela dit, le papier de de Finetti est probablement moins limpide que celui de Markowitz, car, si on retrouve une équation de la forme

il faudra attendre le théorème de Kuhn-Tucker datant de 1951 (http://projecteuclid.org/… même si on pourrait rappeler que William Karush avait obtenu les mêmes conditions dans son mémoire de maîtrise, en 1939) pour avoir les outils mathématiques pour résoudre ce problème de manière élégante. Cela dit, de Finetti propose une résolution intéressante, compte tenu des outils mathématiques à sa disposition avant guerre (comme note dans Pressacco & Serafini (2007) à l’aide d’excursions dans le simplexe), et introduit une notion d’optimalité (au sens de Pareto) qui n’est pas sans rappeler la notion de portefeuille efficient. Presque plus intéressant: de Finetti ayant une formation d’actuaire, il évoque l’idée de mesurer le risque non pas en calculant une variance, mais en calculant une probabilité de ruine, idée que l’on retrouvera plusieurs années plus tard chez Arthur Roy, dans sa notion de safety first. On y retrouve aussi une discussion axiomatique sur les mesures de risques, on y voit même passer des utilités exponentielles, car elles présentent une forme d’aversion pour le risque constante. Mais c’est surtout dans les travaux des années 50, en particulier sulla preferibilità que l’on voit analysé les utilités CARA et HARA (souvent attribuées a Arrow (1971) et Pratt (1970), une dizaine d’années après les travaux de Bruno de Finetti). Là encore, Rubinstein le notait lorsqu’il écrivait

in 1952, anticipating Kenneth Arrow and John Pratt by over a decade, he formulated the notion of absolute risk aversion, used it in connection with the risk premia for small bets and discussed the special case of constant risk aversion“.

Dois-je mentionner sul significato soggettivo della probabilità publié en 1931, sur des calculs de majorants de probabilités de ruine (où il obtient les mêmes résultats que ceux d’Harald Cramér, publié en 1930 dans on the mathematical theory of risk) ? ou la teoria del rischio e il problema della “rovina dei giocatori” publié en 1939 qui discute des notions de martingales et de temps d’arrêt ?

de Finetti est incontestablement un des plus grands mathématiciens et économistes du XXème siècle. Mais beaucoup de ses travaux ont été ignorés, et sont cités ceux publiés par des américains. Ou des chercheurs qui ont fait l’effort de publier en anglais. Au delà de la langue (proprement dite), on pourrait arguer qu’il s’agit d’un problème de discipline, les économistes (et les financiers) connaissant souvent peu les travaux de sciences actuarielles. De plus, Bruno de Finetti avait passé du temps dans les années 50 aux Etats-Unis, pour travailler avec Leonard Savage. Donc dire que les américains ne le connaissaient pas serait une fausse réponse. En fait, il me semble que tous les travaux de de Finetti semblent corroborer la loi de Stigler, qui proclame que « une découverte scientifique ne porte jamais le nom de son auteur » (Steven Stigler se demandant d’ailleurs qui avait trouvé la formule de Bayes, dans http://www.jstor.org/…, même si les mauvaises langues disent que Robert Merton avait enoncé cette loi avant lui).

Dans le contexte actuel où beaucoup de monde s’interroge sur le rôle des éditeurs, ne pourrait-on pas y voir là une solution, qui tirerait le monde de l’édition vers le haut ? Pourquoi les auteurs ne seraient-il pas autorisés à soumettre l’article dans la langue qu’ils souhaitent, et le rôle d’un editeur serait alors de faire traduire puis de diffuser ces travaux. Sur le site, on pourrait lire les articles dans la langue originale, mais aussi en anglais, ou, pourquoi pas, en français ? Tout le monde y gagnerait me semble-t-il ? En tous les cas, nous serions nombreux à décomplexer, et à écrire des articles plus subtiles… Mais si je pense à voix haute, il n’en reste pas moins que régulièrement, des débats sur la langue apparaissent dans le monde académique, et que les réponses sont toujours un peu les mêmes. Par exemple, avoir une ouverture internationale dans une formation se traduit toujours par “dispenser des cours en anglais“. Et répondre “ben fais donc un effort pour parler anglais” me semble une réponse qui pourrait être repensée…

In three months, I’ll be in Vegas (trying to win against the house)

In fact, I’m going there with my family and some friends, including two probabilists (I mean professionals, I am merely an amateur), with this incredible challenge: will I be able to convince  probabilists to go to play at the Casino?

Actually, I also want to study them carefully, to understand how we should play optimally. For example, I hope I can make them play the roulette. Roulette is simple. With a French (or European) roulette, it is probably the simplest: if I bet on black, I win if one of 18 black numbers is out, and I lose if one of the 18 red numbers – or zero (which is green) – is out. This gives a winning probability of 18/37 i.e. a 48.64% chance. But in Vegas, I think it’s mostly American Roulette that can be found in casinos, in which there is a zero and a double zero (both favorable to the bank). Here, the  probability of winning is 18/38, i.e. 47.36% chance. The two roulettes are

Now, let us discuss a little bit about optimal strategy. For instance, suppose I go to Las Vegas with an initial wealth  (say $100). The goal is to find the strategy which maximizes the probability to leave Las Vegas with https://latex.codecogs.com/gif.latex%20?2s (here $200). Should I play big, or small ?

Assume that I can bet https://latex.codecogs.com/gif.latex%20?x (that will be, here, for convenience, a fraction of ). With probability https://latex.codecogs.com/gif.latex%20?p, I will get https://latex.codecogs.com/gif.latex%20?2x, and with probability https://latex.codecogs.com/gif.latex%20?1-p, I will get https://latex.codecogs.com/gif.latex%20?0 (and lose my https://latex.codecogs.com/gif.latex%20?x). As mentioned above, https://latex.codecogs.com/gif.latex%20?p is (a little) smaller than 50%. The casino must win (actually, we will see that this assumption has a very strong impact on the strategy).

Suppose my goal is to double my initial sum, as mentioned in the introduction of this post. Maybe there is an optimum value for https://latex.codecogs.com/gif.latex%20?x, to maximize the probability of doubling my bet. To make it simple, the game ends either because I did not, or because I did, manage to double my initial wealth… Assume further that https://latex.codecogs.com/gif.latex%20?x is fixed, and that I do not revise my bets. One can use monte carlo simulations, to get an intuitive idea…

> bet=function(s=1,t=2*s,x=s/4,p=.4736,nsim=100000){
+     vp=rep(0,nsim); #vw=s
+     for(i in 1:nsim){
+       w=s;
+       while((w>0)&(w<t)){
+          ux=sample(c(min(x,t-w),-x),size=1,prob=c(p,1-p))
+          w=w+ux
+       }
+       vp[i]=(w>=t)}
+     return(mean(vp))
+ }

If we plot this probability as a function of , we have the following

> BET=function(x) bet(x=x)
> vx=1/(1:20)
> px= Vectorize(BET)(vx)
> plot(vx,px,log="x")

Let us see if we can do the maths, and actually compute those probabilities.

For example, if https://latex.codecogs.com/gif.latex%20?x%20=%20s, I play everything I have, and I double with probability https://latex.codecogs.com/gif.latex%20?p. That one was simple.  And indeed, on the graph above, the point on the right is probability  https://latex.codecogs.com/gif.latex%20?p (the red horizontal line).

Assume now that I can bet https://latex.codecogs.com/gif.latex%20?x%20=%20s%20/%202, and I will play, at least, two rounds

  • with probability  I will lose both rounds (and the game is over)
  • with probability , I will win both rounds, and I double my bet (and the game is also over)
  • with probability , I will lose once, and double once. Anyway, I will find myself again with my (initial) wealth . So the game will start again….

To make the story short the probability of doubling my earnings is

which is

https://latex.codecogs.com/gif.latex%20?p%20^%202%20\left%20(1%20+2%20p%20(1-p)%20+%20[2p%20(1-p)]%20^%202%20+%20\cdots%20\right)%20=%20\frac%20{p%20^%202}%20{1-2p%20(1-p)}

Let’s try something more general: I have initial wealth , I can bet https://latex.codecogs.com/gif.latex%20?x and the goal is to reach  (or, more generally, say, ). Now, the probability to reach  from  betting (always) https://latex.codecogs.com/gif.latex%20?x is exactly the same as the probability to reach  from  betting only 1. Let  denote the probability to go from  to  betting 1 (let us use generic parameters). We can easily get the following equation

https://latex.codecogs.com/gif.latex%20?P_b(a)%20=%20p\cdot%20P_b(a+1)%20+%20(1-p)%20\cdot%20P_b(a-1)

Thus, we can write

https://latex.codecogs.com/gif.latex%20?p\cdot%20(P_b(a+1)-P_b(a))%20=%20(1-p)\cdot%20(P_b(a)-P_b(a-1))

or equivalently

https://latex.codecogs.com/gif.latex%20?(P_b(a+1)-P_b(a))%20=\frac{1-p}{p}\cdot%20(P_b(a)-P_b(a-1))

https://latex.codecogs.com/gif.latex%20?\left(\frac{1-p}{p}\right)^a\cdot%20(P_b(1)-P_b(0))

Now, observe that https://latex.codecogs.com/gif.latex%20?P_b(0)=0 (since I cannot have a gain without any money).

Let us write https://latex.codecogs.com/gif.latex%20?P_b(a+1)-P_b(0) using a domino technique :

https://latex.codecogs.com/gif.latex%20?[P_b(a+1)-P_b(a)]+[P_b(a)-P_b(a-1)]+\cdots+[P_b(1)-P_b(0)]

i.e.

https://latex.codecogs.com/gif.latex%20?\left(\frac{1-p}{p}\right)^a%20P_b(1)+\left(\frac{1-p}{p}\right)^{a-1}%20P_b(1)+\cdots+%20\left(\frac{1-p}{p}\right)^0%20P_b(1)

so this geometric sum can also be written

https://latex.codecogs.com/gif.latex%20?\left(1%20-\left[\frac{1-p}{p}\right]^{a+1}%20\right)%20\left(1%20-\left[\frac{1-p}{p}\right]%20\right)^{-1}

Finally, we can write

https://latex.codecogs.com/gif.latex%20?P_b(a)=\left(1%20-\left[\frac{1-p}{p}\right]^{a}%20\right)\left(1%20-\left[\frac{1-p}{p}\right]%20\right)^{-1}\cdot%20P_b(1)

Here, there is still  that I have to explicit. The idea is to observe that , thus

https://latex.codecogs.com/gif.latex%20?P_b(a)=\left(1%20-\left[\frac{1-p}{p}\right]^{a}%20\right)\left(1%20-\left[\frac{1-p}{p}\right]^{b}%20\right)^{-1}

So finally,

https://latex.codecogs.com/gif.latex%20?\mathbb{P}(gain)=\left(1%20-\left[\frac{1-p}{p}\right]^{s/x}%20\right)\left(1%20-\left[\frac{1-p}{p}\right]^{2s/x}%20\right)^{-1}

Nice isn’t it? But to be honest, there is nothing new here. This is actually an old theorem discovered by Christiaan Huygens in 1657, then extended by Jacob Bernoulli in 1680 and finally properly established by Abraham de Moivre in 1711. It is possible to plot this graph, as a function of ,

> bet2=function(s=1,t=2*s,x=s/4,p=.4736){
+     vp=(1-((1-p)/p)^(s/x))/(1-((1-p)/p)^(t/x))
+     return(vp)
+ }

The graph is the same as the one with monte carlo simulation (hopefully). Observe, looking carefully at the function above, that the probability is decreasing with . Which makes sense… Further, the probability is decreasing with : the more hungry, the less chance of winning I have.

Now, the interesting part is what is plotted on the graphs above: the smaller  (the size of the bets at each round), the less chances to win: if I want to win, it is important not to play being little player ! I must bet everything I have ! Actually, the funny thing is that if the probability of winning was (slightly) larger than 1/2, on the contrary, I should bet as small as possible

So far, there is nothing new. Everything mentioned in this post can be related to a fundamental result of Lester Dubins and Leonard Savage, in “How to Gamble if You Must : Inequalities for Stochastic Processes” (published in 1965), see also Sudderth (1972). Of course, I can try another strategy, a little less reasonable, I think, which is sometimes called Martingale of D’Alembert. I believe more in luck than coincidence, so, when I win, I drop my bet (do not tempt fate) but when I lose, I increase my bet (I must win someday). But let’s keep it for another post, someday…

Again, that’s a theory. I guess we should try, and see how it works. I’ll try to upload pictures on the blog during the road trip, so if by the beginning in August nothing has been posted on the blog, please send a rescue team to save me at the Bellagio…

Pavage et magie, rotations et symétries

(oui, mon fils voulait une rime dans le titre… promis on va tenter de faire une comédie musicale à partir de mes billets). La semaine passée, un article intéressant était publié (http://www.bbc.co.uk/news/…) sur sciences et magie. Cet article m’a fait penser à une réflexion que je m’étais faite lors de notre visite au MoMaths, cet hiver. Pour revenir sur l’histoire que j’avais raconté il y a quelques mois, sur la fin de la visite, les enfants avaient joué avec les pavages, à la Escher (les pièces étaient magnétiques, et s’assemblaient sur un grand tableau blanc)

Le soucis c’est qu’à aucun moment je n’ai pu voir d’explication du  pourquoi ça marche ? ou pour être plus explicite comment a été construit le singe pour que tout s’emboîte aussi bien. Car il y a un truc… La magie ça n’existe pas ! ce sont juste des maths (c’est un peu ce que racontent Persi Diaconis et Ron Graham dans Magical Mathematics, même si j’avoue ne pas l’avoir fini, faute de temps, et d’avoir laissé des chapitres en attente). Mon fils travaillant les rotations et les symétries, à l’école, j’ai voulu qu’on prenne le temps de revenir sur les pavages. Mais avant de me lancer avec lui, j’ai essayé de m’entraîner un peu tout seul (dimanche après midi, pendant le temps calme).

La difficulté (ou la magie) est qu’il y a des symétries dans tous les sens. Regardons un exemple simple (ou disons pas trop compliqué) avec des chinois (identifiés par les chapeaux dits chinois). En regardant rapidement, on note que le triangle ci-dessous est reproduit à l’infini, par un jeu de rotations. Il suffit donc d’avoir un tel triangle. Le soucis est que le découpage (parfaitement symétrique) n’est pas trivial (il est délicat de faire un personnage en faisant un découpage dans un triangle, même si c’est la méthode retenue dans l’activité mentionnée sur https://irem.univ-lille1.fr/…)

Le secret part de l’observation suivante: sur le graphique de gauche, en dessous, les trois points (qui sont des intersections des trois personnages, ou disons des trois couleurs) sont les sommets d’un triangle équilatéral. C’est là que la magie va opérer, car c’est cette figure (avec une autre, juxtaposée) que l’on va tenter de reproduire. En fait, un personnage (disons un petit chinois bleu) c’est deux triangles collés (on doit pouvoir parler de parallélogramme isocèle).

On ne le voit pas très bien, mais pourtant, au départ, on avait un parallélogramme isocèle (dont les contours sont tracés en noir) qui était  tout bleu. Bien entendu, on a découpé des parties et en les a envoyé ailleurs (par rotation autour d’un sommet). Le plus simple est de construire la figure, à la main. On part donc du parallélogramme suivant, tout bleu :

Oui, j’ai repris les fonctions de base utilisées l’autre jour pour expliquer R aux enfants,

source("http://freakonometrics.free.fr/RforKIDZ.R")

En fait, j’ai souhaité définir une fonction créant un polynôme à partir d’une liste de points (ou de sommets).

polyg=function(L,couleur="blanc",tour=couleur,ep=1,st=1,...){
  b=unlist(L)
  n=length(b)
  polygon(b[seq(1,n,by=2)],b[seq(2,n,by=2)],col=couleurc(couleur),border=couleurc(tour),lwd=ep,lty=st,...)}

La fonction pour faire une rotation d’une liste de points est ici

rotation=function(L,angle=90,centre=L[[1]]){
  b=unlist(L)
  n=length(b)
  X=cbind(b[seq(1,n,by=2)]-centre[1],b[seq(2,n,by=2)]-centre[2])
  Y=matrix(c(cos(angle/180*pi),sin(angle/180*pi),-sin(angle/180*pi),cos(angle/180*pi)),2,2)%*%t(X)
  Z=Y+centre
  L2=L
  for(i in 1:(n/2)){L2[[i]]=Z[,i]}
  return(L2)}

On a alors ici les quatre sommets

A=c(4,0)
B=c(4,10)
C=c(4+sqrt(3)*5,5)
D=c(4+sqrt(3)*5,15)

et on peut tracer le parallélogramme avec

L1=list(A,B,D,C)
polyg(L1,"bleu",tour="noir",ep=3)

On peut alors jouer à découper des parties, puis à les pivoter (on va faire des rotations, quoi). Par exemple, le prend le triangle en haut à droite, et je pivote de 240 degrés autour d’un des sommets du parallélogramme. C’est la figure ci-dessus à droite.  Le triangle, il est défini par le code suivant

E=c(10,10)
F=c(7,10+sqrt(3)) 
L1=list(B,E,F)
polyg(L1,tour="rouge",ep=3)

ce qui enlève la portion – en traçant un triangle blanc – et pour faire la rotation, j’utilise le code suivant

L2=rotation(L1,240)
polyg(L2,"bleu","rouge",ep=3)

(la seconde ligne permet de tracer le triangle bleu). Il va de soi qu’on peut découper des objets plus complexes que des triangles, voire non connexes. Par exemple, je peux découper la forme verte ci-dessous, puis la pivoter par rapport au sommet en haut à gauche, de 240 degrés.

G=c(10,2*sqrt(3)+10)
H=c(4+sqrt(3)*5,2*sqrt(3)+10)
F1=c(8,10+2/sqrt(3))
F2=c(8,10+4/sqrt(3))
L1=list(B,E,F1,F2,B,G,H,D,G)
polyg(L1,"vert",tour=NA)
L2=rotation(L1,240)
polygL(L2,"bleu",tour=NA)

(il n’est pas nécessaire de revenir au point B, mais peu importe). Ou la figure jaune, qui va pivoter autour du sommet en bas à droite, de 120 degrés.

F3=c(8,10-2/sqrt(3))
F4=c(8,10-4/sqrt(3))
G2=c(10,-2*sqrt(3)+10)
H2=c(4+sqrt(3)*5,-2*sqrt(3)+10)
L1=list(H,G,F2,F1,E,F3,F4,G2,H2,H)
polyg(L1,"jaune",tour=NA)
L2=rotation(L1,120,centre=C)
polygL(L2,"bleu",tour=NA)

Avec des ciseaux, ça sera plus simple, mais le but est aussi de voir que programmer, c’est facile. En gros, on a deux fonctions de base: une qui trace un polynôme et l’autre qui fait des rotations.

Si on fait les deux en même temps (à peu près), on enlève toute la partie en haut à droite, et on la répartie: une partie à gauche, et une partie en bas à droite.  On a alors une forme qui ressemble à un chinois. Non ? On peut aller creuser un peu sous les bras, pour améliorer les mains, mais on y est presque, me semble-t-il…

On peut continuer à améliorer, par exemple en enlevant des parties, puis en faisant des rotations par rapport à un des sommets initiaux. Comme le petit triangle jaune, qui va pivoter de 120 degrés par rapport au sommet en bas à droite. Puis le triangle vert par rapport au sommet en bas à gauche (on note qu’on pivote par rapport aux sommets du triangle équilatéral qu’on avait repéré, initialement)

S=c(8-.4,-2*sqrt(3)+10-.4)
L1=list(F4,G2,S,F4)
polyg(L1,"jaune")
L2=rotation(L1,120,centre=C)
polyg(L2,"bleu")

Bon, là on fait une pause, parce qu’on a, par morceaux, notre chinois ! On un paquet de lettre (voire un paquet de sommets qu’il va falloir nommer, suite à nos rotations). On les nomme avec un code qui ressemble à (oui, on aurait pu pu les nommer dès le départ, lors des rotations)

rE=rotation(list(E),240,centre=B)[[1]]
rF1=rotation(list(F1),240,centre=B)[[1]]
rH=rotation(list(H),240,centre=B)[[1]]

Une fois qu’on a identifié nos sommets, on a notre chinois,

chinois=list(B,rE,rF1,LG,S2,LF2,rH,A,r2H,r2G,r2S2,
        r2F2,r2F1,r2E,rF3,rF4,r2S,rG2,rH2,C,H2,G2,S,F4,F3,E,B)

Bon, ensuite, reste à paver, ce qui est est le plus simple. Par exemple, on peut faire des rotations,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")

voire définir une translation,

translation=function(L,h){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]+h[1],b[seq(2,n,by=2)]++h[2])
L2=L
for(i in 1:(n/2)){L2[[i]]=X[i,]}
return(L2)}

puis translater notre chinois,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")
polyg(translation(r1chinois,15/sqrt(3),15),couleur="jaune")

Il manque l’autre, le vert. On peut vérifier qu’il se superpose parfaitement :

Bon, ben ça y est, on a constitué notre pavage. Reste à travailler un peu le code, essentiellement sur la forme, et on pourra jouer avec les enfants, à faire les lézards, des singes, et toutes sortes de pavages possibles, à la Escher. En notant qu’on peut utiliser d’autres formes de bases plus complexes, comme expliqué sur http://xavier.hubaut.info/…

Overdispersed Poisson et bootstrap

Pour le dernier cours sur les méthodes de provisionnement, on s’est arrête aux méthodes par simulation. Reprenons là où on en était resté au dernier billet où on avait vu qu’en faisant une régression de Poisson sur les incréments, on obtenait exactement le même montant que la méthode Chain Ladder,

> Y
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   53   22   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

> y=as.vector(as.matrix(Y))
> base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
> reg2=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson) 
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson, 
    data = base)

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        8.05697    0.01551 519.426  < 2e-16 ***
as.factor(ai)2001  0.06440    0.02090   3.081  0.00206 ** 
as.factor(ai)2002  0.20242    0.02025   9.995  < 2e-16 ***
as.factor(ai)2003  0.31175    0.01980  15.744  < 2e-16 ***
as.factor(ai)2004  0.44407    0.01933  22.971  < 2e-16 ***
as.factor(ai)2005  0.50271    0.02079  24.179  < 2e-16 ***
as.factor(bj)1    -0.96513    0.01359 -70.994  < 2e-16 ***
as.factor(bj)2    -4.14853    0.06613 -62.729  < 2e-16 ***
as.factor(bj)3    -5.10499    0.12632 -40.413  < 2e-16 ***
as.factor(bj)4    -5.94962    0.24279 -24.505  < 2e-16 ***
as.factor(bj)5    -5.01244    0.21877 -22.912  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
  (15 observations deleted due to missingness)
AIC: 209.52

Number of Fisher Scoring iterations: 4

> base$py2=predict(reg2,newdata=base,type="response")
> round(matrix(base$py2,n,n),1)
       [,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.7 1202.1 49.8 19.1  8.2 21.0
[2,] 3365.6 1282.1 53.1 20.4  8.8 22.4
[3,] 3863.7 1471.8 61.0 23.4 10.1 25.7
[4,] 4310.1 1641.9 68.0 26.1 11.2 28.7
[5,] 4919.9 1874.1 77.7 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.4 31.6 13.6 34.7

> sum(base$py2[is.na(base$y)])
[1] 2426.985

Le plus intéressant est peut être de noter que la loi de Poisson présente probablement trop peu de variance…

> reg2b=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = quasipoisson, 
    data = base)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        8.05697    0.02769 290.995  < 2e-16 ***
as.factor(ai)2001  0.06440    0.03731   1.726 0.115054    
as.factor(ai)2002  0.20242    0.03615   5.599 0.000228 ***
as.factor(ai)2003  0.31175    0.03535   8.820 4.96e-06 ***
as.factor(ai)2004  0.44407    0.03451  12.869 1.51e-07 ***
as.factor(ai)2005  0.50271    0.03711  13.546 9.28e-08 ***
as.factor(bj)1    -0.96513    0.02427 -39.772 2.41e-12 ***
as.factor(bj)2    -4.14853    0.11805 -35.142 8.26e-12 ***
as.factor(bj)3    -5.10499    0.22548 -22.641 6.36e-10 ***
as.factor(bj)4    -5.94962    0.43338 -13.728 8.17e-08 ***
as.factor(bj)5    -5.01244    0.39050 -12.836 1.55e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasipoisson family taken to be 3.18623)

    Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
  (15 observations deleted due to missingness)
AIC: NA

Number of Fisher Scoring iterations: 4

Mais on en reparlera dans un instant. Ensuite, on avait commencé à regarder erreurs commises, sur la partie supérieure du triangle. Classiquement, par construction, les résidus de Pearson sont de la forme

https://latex.codecogs.com/gif.latex?\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\text{Var}(Y_i)}}

On avait vu dans le cours de tarification que la variance au dénominateur pouvait être remplacé par le prévision, puisque dans un modèle de Poisson, l’espérance et la variance sont identiques. Donc on considérait

https://latex.codecogs.com/gif.latex?\varepsilon_i=\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}

> base$erreur=(base$y-base$py2)/sqrt(base$py2)
> round(matrix(base$erreur,n,n),1)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]  0.9 -1.1 -1.5 -0.5 -0.4    0
[2,]  0.0  0.3 -2.2  0.8  0.4   NA
[3,]  0.1  0.1 -1.0 -0.3   NA   NA
[4,] -1.1  0.9  4.2   NA   NA   NA
[5,]  0.1 -0.2   NA   NA   NA   NA
[6,]  0.0   NA   NA   NA   NA   NA

Le soucis est que si https://latex.codecogs.com/gif.latex?\widehat{Y}_i est – asymptotiquement – un bon estimateur pour https://latex.codecogs.com/gif.latex?\text{Var}(Y_i), ce n’est pas le cas à distance finie, car on a alors un estimateur biaisé pour la variance, et donc la variance des résidus n’a que peu de chances d’être de variance unitaire. Aussi, il convient de corriger l’estimateur de la variance, et on pose alors

https://latex.codecogs.com/gif.latex?\varepsilon_i=\sqrt{\frac{n}{n-k}}\cdot\frac{Y_i-\widehat{Y}_i}{\sqrt{\widehat{Y}_i}}

qui sont alors les résidus de Pearson tel qu’on doit les utiliser.

> E=base$erreur[is.na(base$y)==FALSE]*sqrt(21/(21-11))
> E
 [1]  1.374976e+00  3.485024e-02  1.693203e-01 -1.569329e+00  1.887862e-01
 [6] -1.459787e-13 -1.634646e+00  4.018940e-01  8.216186e-02  1.292578e+00
[11] -3.058764e-01 -2.221573e+00 -3.207593e+00 -1.484151e+00  6.140566e+00
[16] -7.100321e-01  1.149049e+00 -4.307387e-01 -6.196386e-01  6.000048e-01
[21] -8.987734e-15
> boxplot(E,horizontal=TRUE)

En rééchantillonnant dans ces résidus, on peut alors générer un pseudo triangle. Pour des raisons de simplicités, on va générer un peu rectangle, et se restreindre à la partie supérieure,

> Eb=sample(E,size=36,replace=TRUE)
> Yb=base$py2+Eb*sqrt(base$py2)
> Ybna=Yb
> Ybna[is.na(base$y)]=NA
> Tb=matrix(Ybna,n,n)
> round(matrix(Tb,n,n),1)
       [,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3115.8 1145.4 58.9 46.0  6.4 26.9
[2,] 3179.5 1323.2 54.5 21.3 12.2   NA
[3,] 4245.4 1448.1 61.0  7.9   NA   NA
[4,] 4312.4 1581.7 68.7   NA   NA   NA
[5,] 4948.1 1923.9   NA   NA   NA   NA
[6,] 4985.3     NA   NA   NA   NA   NA

Cette fois, on a un nouveau triangle ! on va alors pouvoir faire plusieurs choses,

  1. compléter le triangle para la méthode Chain Ladder, c’est à dire calculer les montants moyens que l’on pense payer dans les années futures
  2. générer des scénarios de paiements pour les années futurs, en générant des paiements suivant des lois de Poisson (centrées sur les montants moyens que l’on vient de calculer)
  3. générer des scénarios de paiements avec des lois présentant plus de variance que la loi de Poisson. Idéalement, on voudrait simuler des lois qusi-Poisson, mais ce ne sont pas de vraies lois. Par contre, on peut se rappeler que dans ce cas, la loi Gamma devrait donner une bonne approximation.

Pour ce dernier point, on va utiliser le code suivant, pour générer des quasi lois,

> rqpois = function(n, lambda, phi, roundvalue = TRUE) {
+ b = phi
+ a = lambda/phi
+ r = rgamma(n, shape = a, scale = b)
+ if(roundvalue){r=round(r)}
+ return(r)
+ }

Je renvois aux diverses notes de cours pour plus de détails sur la justification, ou à un vieux billet. On va alors faire une petite fonction, qui va soit somme les paiements moyens futurs, soit sommer des générations de scénarios de paiements, à partir d’un triangle,

> CL=function(Tri){
+ y=as.vector(as.matrix(Tri))
+ base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
+ reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=quasipoisson)
+ py2=predict(reg,newdata=base,type="response")
+ pys=rpois(36,py2)
+ pysq=rqpois(36,py2,phi=3.18623)
+ return(list(
+ cl=sum(py2[is.na(base$y)]),
+ sc=sum(pys[is.na(base$y)]),
+ scq=sum(pysq[is.na(base$y)])))
+ }

Reste alors à générer des paquets de triangles. Toutefois, il est possible de générer des triangles avec des incréments négatifs. Pour faire simple, on mettra des valeurs nulles quand on a un paiement négatif. L’impact sur les quantiles sera alors (a priori) négligeable.

> for(s in 1:1000){
+ Eb=sample(E,size=36,replace=TRUE)*sqrt(21/(21-11))
+ Yb=base$py2+Eb*sqrt(base$py2)
+ Yb=pmax(Yb,0)
+ scY=rpois(36,Yb)
+ Ybna=Yb
+ Ybna[is.na(base$y)]=NA
+ Tb=matrix(Ybna,6,6)
+ C=CL(Tb)
+ VCL[s]=C$cl
+ VR[s]=C$sc
+ VRq[s]=C$scq
+ }

Si on regarde la distribution du best estimate, on obtient

> hist(VCL,proba=TRUE,col="light blue",border="white",ylim=c(0,0.003))
> boxplot(VCL,horizontal=TRUE,at=.0023,boxwex = 0.0006,add=TRUE,col="light green")
> D=density(VCL)
> lines(D)
> I=which(D$x<=quantile(VCL,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
> I=which(D$x>=quantile(VCL,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

Mais on peut aussi visualiser des scénarios basés sur des lois de Poisson (équidispersé) ou des scénarios de lois quasiPoisson (surdispersées), ci-dessous

Dans ce dernier cas, on peut en déduire le quantile à 99% des paiements à venir.

> quantile(VRq,.99)
    99% 
2855.01

Il faut donc augmenter le montant de provisions de l’ordre 15% pour s’assurer que la compagnie pourra satisfaire ses engagements dans 99% des cas,

> quantile(VRq,.99)-2426.985
    99% 
428.025

Réassurrance

Mercredi aura lieu le dernier cours d’actuariat IARD.

Parmi les compléments, Introduction à la réassurance, publié par Swiss Re, ou ainsi que quelques documents plus techniques, comme The Pareto model in property reinsurance , Exposure rating, ou Designing property reinsurance programmes encore Introduction to reinsurance accounting. Plusieurs réassureurs (et courtiers en réassurance) publient des études techniques sur leurs sites, http://swissre.com/http://munichre.com/, http://aon.com/, http://scor.com/ ou encore http://guycarp.com/. Sinon je renvois aux notes de cours de Peter Antal, quantitative methods in reinsurance.

Les transparents sont en ligne ici,

Sur les sinistres les plus chers, pour les compagnies d’assurance et de réassurance, http://businessinsider.com/… donne le classement suivant en dollars de 2010 (on pourra aussi consulter http://media.swissre.com/…)
  1. Hurricane Katrina (US, Bahamas, Cuba, Aug. 2005), $ 72.3 billion
  2. Tōhoku earthquake and tsunami (Japan, March 2011), $ 35 billion
  3. Hurricane Andrew (US, Bahamas, August 1992), $ 25 billion
  4. September 11 attacks (US) $ 23.1 billion
  5. Northridge earthquake (US) $ 20.6 billion
  6. Hurricane Ike (US, Haiti, Dominican Republic, Sept. 2005) $ 20.5 billion
  7. Hurricane Ivan (US, Barbados, Sept. 2004) $ 14.9 billion
  8. Hurrican Wilman (US, Mexico, Jamaica, Oct. 2005), $ 14 billion
  9. Hurricane Rita (US, Cuba, Sept. 2005) $ 11.3 billion
  10. Hurricane Charley (US, Cuba, Jamaica) $ 9.3. billion

A titre de comparaison, les chiffres d’affaires des plus gros réassureurs (prime émise en 2010) étaient, selon http://www.insurancenetworking.com/…

  1. Munich Reinsurance Company $ 31.3 billion
  2. Swiss Reinsurance Company Limited $ 24.7 billion
  3. Hannover Rueckversicherung AG $ 15.1 billion
  4. Berkshire Hathaway Inc. $ 14.4 billion
  5. Lloyd’s $ 13 billion
  6. SCOR S.E. $  8.8 billion
  7. Reinsurance Group of America Inc. $ 7.2 billion
  8. Allianz S.E. $ 5.7 billion
  9. PartnerRe Ltd. $ 4.9 billion
  10. Everest Re Group Ltd. $ 4.2 billion

Reserving with negative increments in triangles

A few months ago, I did published a post on negative values in triangles, and how to deal with them, when using a Poisson regression (the post was published in French). The idea was to use a translation technique:

  1. Fit a model not on https://latex.codecogs.com/gif.latex?Y_i‘s but on https://latex.codecogs.com/gif.latex?Y_i^{(k)}=Y_i+k, for some https://latex.codecogs.com/gif.latex?k\geq%200,
  2. Use that model to make predictions, and then translate those predictions, https://latex.codecogs.com/gif.latex?\widehat{Y}_i^{(k)}-k

This is what was done to get the following graph, where a Poisson regression was fitted. Black points are https://latex.codecogs.com/gif.latex?Y_i‘s while blue points are https://latex.codecogs.com/gif.latex?\widehat{Y}_i^{(k)}‘s, for some https://latex.codecogs.com/gif.latex?k\geq%200. We fit a model to get the blue prediction, and then translate it to get the red prediction (on the https://latex.codecogs.com/gif.latex?Y_i‘s).
http://freakonometrics.blog.free.fr/public/perso4/glm-translation.gif

In this example, there were no negative values, but it is possible to use it get a better understanding on the impact of this technique. The prediction, here, is the red line. And clearly, the value of https://latex.codecogs.com/gif.latex?k has an impact on the prediction (since we do not consider, here, a linear model: with a linear model, translating has not impact at all, except on the intercept).

The alternative mentioned in the previous post was to use this technique on several https://latex.codecogs.com/gif.latex?k‘s, and them interpolate

  1. For a given https://latex.codecogs.com/gif.latex?k, fit a model not on https://latex.codecogs.com/gif.latex?Y_i‘s but on https://latex.codecogs.com/gif.latex?Y_i^{(k)}=Y_i+k, use that model to make predictions, and then translate those predictions, https://latex.codecogs.com/gif.latex?\widehat{Y}_i^{(k)}-k.
  2. Do it for several https://latex.codecogs.com/gif.latex?k‘s.
  3. Use it to extrapolate when https://latex.codecogs.com/gif.latex?k is https://latex.codecogs.com/gif.latex?0 (which is the case we are interested in).

In the context of loss reserving, the idea is extremely simple. Consider a triangle with incremental payments

> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")
> Y=T=PAID
> n=ncol(T)
> Y[,2:n]=T[,2:n]-T[,1:(n-1)]   
> Y
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   53   22   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

Now, we do not have negative values, here, but we can still see is translation techniques can be used. The benchmark is the Poisson regression, since we can run it :

> y=as.vector(as.matrix(Y))
> base=data.frame(y,ai=rep(2000:2005,n),bj=rep(0:(n-1),each=n))
> reg=glm(y~as.factor(ai)+as.factor(bj),data=base,family=poisson)

Here, the amount is reserve is the sum of predicted values in the lower part of the triangle,

> py=predict(reg,newdata=base,type="response")
> sum(py[is.na(base$y)])
[1] 2426.985

which is exactly Chain Ladder’s estimate.

Now, let us use a translation technique to compute the amount of reserves. The code will be

> decal=function(k){
+ reg=glm(y+k~as.factor(ai)+as.factor(bj),data=base,family=poisson)
+ py=predict(reg,newdata=base,type="response")
+ return(sum(py[is.na(base$y)]-k))

For instance, if we translate of +5, we would get

> decal(5)
[1] 2454.713

while a translation of +10 would return

> decal(10)
[1] 2482.29

Clearly, translations do have an impact on the estimation. Here, just to check, if we do not translate, we do have Chain Ladder’s estimate,

> decal(0)
[1] 2426.985

The idea mentioned in the previous post was to try several translations, and then extrapolate, to get the value in 0. Here, translations will give the following estimates

> K=10:20
> (V=Vectorize(decal)(K))
 [1] 2482.290 2487.788 2493.279 2498.765 2504.245 2509.719 2515.187 2520.649
 [9] 2526.106 2531.557 2537.001

We can plot those values, and run a regression

> plot(K,V,xlim=c(0,20),ylim=c(2425,2540))
> abline(h=decal(0),col="red",lty=2)

the dotted horizontal line is Chain Ladder. Now, let us extrapolate

> b=data.frame(K=K,D=V)
> rk=lm(D~K,data=b)
> predict(rk,newdata=data.frame(K=0))
       1 
2427.623

On has to admit that it is not that bad. But yesterday evening, Karim asked me why I did use a linear regression, for my extrapolation. And to be honest, I do not know. I mean, the only answer might be that points are almost on a straight line. So the first time I saw it, I was exited, and I ran a linear regression.

Now, let us see if we can do better. Because here, we do use a translation of +10 or +20 (which might be rather small). What if we use much larger values ? (because we might have large negative incremental values). With the following code, we try, each time 11 consecutive values, the smallest one going from 0 to 50,

> hausse=1:50; res=rep(NA,50)
> for(k in hausse){
+ VK=k:(10+k)
+ b=data.frame(K=VK,D=Vectorize(decal)(VK))
+ rk=lm(D~K,data=b)
+ res[k]=predict(rk,newdata=data.frame(K=0))
+ }     
> plot(hausse,res,type="l",col="red",ylim=c(2422,2440))
> abline(rk,col="blue")

Here, we compute reserves when extrapolations were done after 11 translations, from https://latex.codecogs.com/gif.latex?k to https://latex.codecogs.com/gif.latex?k+10.  With different values of https://latex.codecogs.com/gif.latex?k. The case where https://latex.codecogs.com/gif.latex?k is ten was the one mentioned above,

> res[hausse==10]
[1] 2427.623

Actually, it might also be possible to consider not 11 translations, but 26, from https://latex.codecogs.com/gif.latex?k to https://latex.codecogs.com/gif.latex?k+25. Here, we get

> hausse=1:50; res=rep(NA,50)
> for(k in hausse){
+ VK=k:(25+k)
+ b=data.frame(K=VK,D=Vectorize(decal)(VK))
+ rk=lm(D~K,data=b)
+ res[k]=predict(rk,newdata=data.frame(K=0))
+ }   
> lines(hausse,res,type="l",col="blue",lty=2)

We now have the dotted line

Here, it is getting worst. So let us keep here 11 translations. Perhaps, we can try something different. For instance a Poisson regression, with a log like (i.e. we consider an exponential extrapolation),

> hausse=1:50; res=rep(NA,50)
> for(k in hausse){
+ VK=k:(10+k)
+ b=data.frame(K=VK,D=Vectorize(decal)(VK))
+ rk=glm(D~K,data=b,family=poisson)
+ res[k]=predict(rk,newdata=data.frame(K=0),type="response")
+ }         
> lines(hausse,res,type="l",col="purple")

The purple line will be a Poisson model, with a log link. Perhaps we can try another link function, like a quadratic one

> hausse=1:50; res=rep(NA,50)
> for(k in hausse){
+ VK=k:(10+k)
+ b=data.frame(K=VK,D=Vectorize(decal)(VK))
+ rk=glm(D~K,data=b,family=poisson(link=
+ power(lambda = 2)))
+ res[k]=predict(rk,newdata=data.frame(K=0),type="response")
+ }     
> lines(hausse,res,type="l",col="orange")

That would be the orange line,

Here, we need a link function between identity (the linear model, the blue line) and the quadratic one (the orange one), for instance a power function 3/2,

> hausse=1:50; res=rep(NA,50)
> for(k in hausse){
+ VK=k:(10+k)
+ b=data.frame(K=VK,D=Vectorize(decal)(VK))
+ rk=glm(D~K,data=b,family=poisson(link=
+ power(lambda = 1.5)))
+ res[k]=predict(rk,newdata=data.frame(K=0),type="response")
+ }         
> lines(hausse,res,type="l",col="green")

Here, it looks like we can use that model for any kind of translation, from +10 till +50, even +100 ! But I do not have any intuition about the use of this power function…

Les experts

Il y a un an, presque jour pour jour, je finissais un rapport que je remettais pour un litige qui devait passer devant un tribunal. J’avais été sollicité pour faire des simulations de scénarios, et étudier des lois de temps d’arrêt et des problèmes de probabilités de ruine. Bref, j’étais un peu devenu un expert. Peut-être pas exactement un expert judiciaire au sens strict, car j’étais sollicité par une des parties, et non pas pas un juge (est-ce toujours le cas ?). Et je n’ai pas non plus prêté (solennellement) le serment “Je jure d’apporter mon concours à la Justice, d’accomplir ma mission, de faire mon rapport, et de donner mon avis en mon honneur et en ma conscience”. Même si c’est ainsi que j’ai mené ma mission. Mais peu importe (pour des raisons de confidentialité, je n’en dirais pas plus sur cette mission). C’est toujours valorisant de se dire que les choses théoriques que l’on enseigne peuvent servir dans la vraie vie !

Hier, au détour d’une conversation avec mart1oeil, goulu ‏et tomroud, j’ai découvert zythom et son blog, http://zythom.blogspot.ca/, qui est le “blog d’un informaticien expert judiciaire“, comme l’indique l’entête. Et hier, zythom posait la question un expert judiciaire peut-il être mauvais dans votre domaine ? Je reviens sur son billet car le billet commence par une recherche de définition de ce que peut être un expert. Et comme toujours la page de Wikipédia consacrée à l’expert propose une lumière intéressante. “L’expert n’est pas simplement celui qui sait, sur un champ délimité de savoir. Son expérience reconnue lui permet d’apporter une réponse argumentée à une demande d’expertise. Il faut le différencier du savant et aussi du spécialiste.” Et zythom prend alors l’image du médecin de famille ou du médecine généraliste, en opposition au spécialiste. Bref, tout d’un coup, hier un monde s’est effondré ! Je suis probablement pas un expert !

Cela dit, j’aurais pu m’en douter plus tôt. En fait, je crois que je le savais depuis des années, mais que je me voilais la face.

Il y a quelques années maintenant, j’enseignais un cours de méthodes numériques en finance (j’avais même tapé des notes de cours). A un moment, on montre que le prix https://latex.codecogs.com/gif.latex?V vérifie une belle équation, qui permettra d’arriver à la formule de Black-Scholes, qui permet d’avoir le prix d’une option dite européenne. L’équation est de la forme

https://latex.codecogs.com/gif.latex?\displaystyle%20\frac{\partial%20V}{\partial%20t}%20+%20\frac{1}{2}\sigma^2%20S^2\frac{\partial^2%20V}{\partial%20S^2}%20+%20rS\frac{\partial%20V}{\partial%20S}%20-%20rV%20=%200

En bricolant un peu (avec quelques changements de variables), on peut se ramener une équation de la forme

https://latex.codecogs.com/gif.latex?\frac{\partial%20u(t,x)}{\partial%20t}-\kappa%20\frac{\partial%20^2%20u(t,x)}{\partial%20x^2}=f(t,x)

Tous ceux qui ont fait un peu de maths reconnaîtrons une équation classique (datant de plus de 200 ans) appelée équation de la chaleur (les plus curieux pourront aller voir les notes de cours de Nicolas Vandewalle sur le sujet, entre autres). Et dans mon cours, je montrais comment résoudre numériquement cette équation, avec différents schémas, comme un schéma explicite

https://latex.codecogs.com/gif.latex?\frac{u_{i+1,j}-u_{i,j}}{\Delta%20t}-\kappa%20\frac{u_{i,j+1}-2%20u_{i,j}+u_{i,j-1}}{h^2}=f_{i,j},

en insistant sur les conditions numériques pour que la résolution n’explose pas (comme sur le dessin ci-dessous)

Il fallait pour cela que le ratio

https://latex.codecogs.com/gif.latex?\lambda=\kappa\frac{\Delta%20t}{h^2}

soit inférieur à 1/2. Etc.

Bref, je suis capable de résoudre numériquement l’équation de la chaleur depuis des années, voire d’enseigner comment le faire. N’empêche que je lorsque je réchauffe mon lunch au micro-onde, tous les midis, je me brule, en ayant en même temps un morceau qui est encore congelé dans le plat ! C’est pénible ! Vivement que je devienne un expert !

Martingale et journalisme scientifique

Être journaliste scientifique ne doit pas être facile. J’imagine qu’il faut être à l’écoute des nouvelles scientifiques, et d’informer, aussi justement que possible. C’est ce que fait avec brio Pierre Barthélémy (aka @PasseurSciences sur Twitter) dans sa chronique hebdomadaire sélection scientifique de la semaine. Il essaye ainsi de parler de sciences dans un journal qui a (trop) souvent confondu science et technologie (technologie étant aussi souvent un mot savant utilisé pour masquer de la publicité pour des appareils technologiques avancés). Et comme l’espace qui lui est imparti est restreint (c’est le moins qu’on puisse dire), la plupart des informations passent par son blog http://passeurdesciences.blog.lemonde.fr/, qui est riche ! Sylvestre Huet fait aussi ce travail pour Libération, avec son blog http://sciences.blogs.liberation.fr/ (là encore, on peut lire bien plus de choses en ligne que que dans la version papier, ce qui est d’autant plus intéressant pour les lecteurs de l’autre coté de l’océan atlantique). Sylvestre Huet arrive à parler géologie, puis démographie, avant d’évoquer dans un article truffé de références les réformes des instances qui régentent la recherche en France. Oui, les deux blogs sont malgré tout très français. Pour une vision plus internationale, on pourra citer l’admirable travail de Pascal Lapointe, par exemple (aka @paslap) sur http://sciencepresse.qc.ca/.

Mais si je commence à mettre les pieds dans le plat, je poserais la question de la légitimité d’un journaliste scientifique, à l’heure où la blogosphère scientifique explose. Par exemple, si je veux apprendre des choses en physique, ou en biologie, je sais que le blog de Tom Roud (que l’on peut suivre sur Twitter @tomroud) http://tomroud.cafe-sciences.org est une source infinie d’information (ok, peut-être pas “infinie” car comme nous tous, il a une vie en dehors de son blog, mais disons qu’il y a de quoi lire). En sciences humaines, il existe des centaines de blogs hébergés sur http://hypotheses.org/ tenus par des universitaires, ou http://www.cafe-sciences.org/ pour des blogs de sciences, en français. Les bloggers scientifiques s’expriment – le plus souvent – en restant dans leur champ d’expertise. Et c’est tant mieux. Le journaliste scientifique, lui, se doit de parler de tous les sujets scientifiques. Et c’est là que l’exercice devient délicat, car le journaliste se doit d’être critique. Compte tenu de la concurrence féroce qui existe dans le monde académique, on nous pousse à faire croire qu’on vient d’écrire l’article qui va révolutionner la science, sinon le monde. Que notre approche est novatrice, et qu’en plus, on vient de montrer ce qui pourrait ce qui pourrait être le sain Graal dans notre communauté. C’est le travail du journaliste de se demander si c’est vrai….

J’en arrive au point de mon billet. Le dernier article  le hasard, martingale boursière ? (en ligne sur http://lemonde.fr/sciences/…) de Pierre Barthélémy m’a un peu agacé. Pas sur le thème abordé, mais l’impression générale qui m’est restée après avoir lu l’article (samedi matin). Pour l’histoire complète, c’est ma femme qui est tombée la première sur l’article, et qui a été étonnée de lire un article pareil dans Le Monde, en tant que mathématicienne (et probabiliste). Et quand je l’ai lu, c’est en tant qu’économiste que je suis resté sans voix. En résumé (trop succinct, j’en convient), on nous résume l’article de physiciens italiens qui “expliquent qu’au grand dam des traders et autres analystes financiers, en quête perpétuelle de justifications rationnelles aux fluctuations des cours boursiers, les marchés demeurent obstinément imprévisibles“. Pour parcourir l’article mentionné, mais non cité – probablement car il s’agit de la version papier ? – l’article est en ligne sur http://arxiv.org/1303.4351  (trouvé via http://improbable.com/…). On retrouve cette idée dans le texte sous la forme suivante “our main result, which is independent of the market considered, is that standard trading strategies and their algorithms, based on the past history of the time series, although have occasionally the chance to be successful inside small temporal windows, on a large temporal scale, perform on average not better than the purely random strategy, which, on the other hand, is also much less volatile. In this respect, for the individual trader, a purely random strategy represents a costless alternative to expensive professional financial consulting, being at the same time also much less risky, if compared to the other trading strategies.” On apprend dans l’article paru dans Le Monde que “pour Pluchino et compagnie, le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre“. Damned, rien que ça ? Pour les amateurs de lyrisme (il doit y en avoir dans les lecteurs du monde, moins parmi les lecteurs plus scientifiques qui aiment les phrases moins alambiquées), je passe les moments où Pierre Barthélémy est clairement plus journaliste que scientifique, comme lorsqu’il nous explique que ces quatre chercheurs “viennent d ‘attaquer à l’acide de l’aléatoire” les modèles financiers. Cela dit, cela a le mérite de mettre les points sur les i: ce sont des physiciens. Pas des économistes.

Je pense qu’il y a des milliers de raisons d’attaquer les traders (si ce sont bien les “spécialistes” visés). Mais qu’on ne les prenne pas pour plus bêtes qu’ils ne le sont (je le dis d’autant plus volontiers que j’en ai formé un paquet, à l’ENSAE ou à Polytechnique, au cours des 10 dernières années).

Avant de revenir sur la petite histoire de martingale financière (tel que le résume le titre de l’article) faisons un court détour par celle brillamment racontée dans Mansuy (en ligne sur http://math.harvard.edu/~ctm/… et traduit en anglais dans http://jehps.net/Mansuy.pdf) sur l’origine de la notion de martingale. D’un point de vue mathématique (disons comme propriété des processus stochastiques), il faut remonter aux travaux de Joseph Bernstein, Paul Lévy, Émile Borel, et surtout Joseph Leo Doob, au milieu au XXième siècle. Cela dit, le mot est plus ancien: il entre dans le dictionnaire de l’Académie Française en 1762. Jouer à la martingale, c’est jouer toujours tout ce que l’on a perdu.  Mais quelques années plus tot, on pouvait déjà trouver le mot sous la plume de l’Abbé Prévost (le jeu qu’il décrit en 1750 comme variante du jeu du pharaon est aussi appelé martingale d’Alembert) Cela dit, l’origine que je préfère est de relier le mot martingale à l’expression provençale “jouga a la martegalo” qui signifierait jouer de manière incompréhensible, absurde, comme l’évoque Frédéric Mistral dans Lou Tresor dòu Felibrige ou dictionnaire provençal-français. On retrouve une origine proche dans le dictionarie of the French and English tongues de Randle Cotgrave, datant de 1611, qui mentionne l’expression “à la martingale” avec le sens absurdly, foolishly, untowardly, grossely, rudely, in the homeliest manner. Il cite même l’usage de l’expression philosopher à la martingale (sans citer Bernard Henri Levy, mais il semble que ce soit l’idée). Cela dit, même dans Lapinot on parle de martingale,

Bref, à partir de ce comportement incompréhensible, voire absurde, des économistes vont définir une notion importante en finance, que l’on appellera efficience (qui est un mot dangereux car il évoque aux oreilles de tous les économiste une notion d’optimalité au sens de Pareto, on pourra relire Malkiel (2003) The Efficient Market Hypothesis and Its Critics sur ce point). L’idée d’utiliser une marche aléatoire pour modéliser le cours d’un actif est ancienne. Par exemple (je ne remontrais pas au XVIème siècle, mais http://e-m-h.org/history.html le suggère) on pourra relire les travaux de Bachelier datant du début du XXième siècle, mais comme le dit la légende, Bachelier a été peu lu, à l’époque (en tous les cas par des économistes). Publié quelques années plus tard, on pourra relire Cowles & Jones (1937), Some A Posteriori Probabilities in Stock Market ou surtout Samuelson (1965) Proof That Properly Anticipated Prices Fluctuate Randomly. Cet article est passionnant, si on prend le temps de le lire “There is no way of making an expected profit by extrapolating past changes in the future price, by chart or any other esoteric devices of magic or mathematics.” Paul Samuelson est d’une modestie remarquable pour un chercheur, “I have not here discussed what the basic probability distributions are supposed to come from. In whose mind are they ex ante? Is there any ex post validation of them? Are they supposed to belong to the market as a whole? And what does that mean? Are they supposed to belong to the “representative individual”, and who is he? Are there some defensible or necessitous compromises of divergent anticipations patterns? Do price quotations somehow produce a Pareto-optimal configuration of ex ante subjective probabilities? This paper has not attempted to pronounce on these interesting questions“. Il peut y avoir des bulles, des processus non Gaussien, un peu tout ce qu’on peut imaginer, le point important étant qu’il est impossible d’utiliser le passé pour prédire le futur, en finance.

Oui, depuis 50 ans, les économistes savent que la marche aléatoire peut être un bon modèle pour décrire le prix des actifs. Ou pour être plus précis (et rendre à Paul Samuelson ce que Paul Samuelson a dit le premier) les martingales. Mais il faudra surtout attendre les travaux d’Eugène Fama, en particulier Fama (1965) The Behavior of Stock Market Prices et Fama (1970) Efficient Capital Markets: A Review of Theory and Empirical Work, qui vont de manière définitive marquer le début de l’utilisation des martingales en finance. Eugène Fama y retient trois notions d’efficience (pour aller plus loin que le fameux “A market in which prices at any time “fully reflect” available information is called “efficient”“)

  1. Expected Returns or “Fair Game” Models
  2. Submartingale Models
  3. Random Walk Model

(en 1965, l’efficience devait être reliée, pour Eugène Fama, à la notion de marche aléatoire, alors que Paul Samuelson utilisait déjà la notion de martingale). Ces modèles disent tous (on est en 1970) qu’il est inutile d’utiliser l’information passée ou présente “to predict the future in a way which makes expected profits greater than they would be under a naive buy-and-hold model“. N’en déplaise à Pierre Barthélémy , depuis presque 50 ans ce principe est énoncé dans les cours d’asset pricing. Bachelier disait la meme chose en 1900, “l’espérance mathématique du spéculateur est nulle” dans Théorie de la spéculation. Et cela sera confirmé quelques années plus tard par les résultats empiriques de Cowles et Jones. Un an plus tard sera publié Black (1971) Implications of the random walk hypothesis for portfolio management, et trois ans plus tard, Hagerman & Richmond (1973) Random Walks, Martingales and the OTC. La littérature économique va se multiplier pendant les années 70, à tel point qu’en 1978, Michael Jensen, alors professeur à Harvard écrivait “I believe there is no other proposition in economics which has more solid empirical evidence supporting it than the Efficient Market Hypothesis” (cité dans http://economist.com/14030296). Plus récemment, à la fin des années 80 était publié le remarquable LeRoy (1989) Efficient Capital Markets and Martingales, qui étoffait la critique initiée dans LeRoy (1976). Efficient capital markets: Comment. LeRoy et Samuelson ont été les premiers à parler de martingale pour modéliser les prix des actifs. On continue (avec les articles qui font référence, et que tous les étudiants qui ont suivi un cours d’asset pricing ont lu) ? L’année suivante était publié Lehman (1990) Fads, Martingales, and Market Efficiency. Les martingales sont des outils incroyablement riches. Ils ne sont pas équivalent à une marche aléatoire: on peut avoir une martingale, sans avoir de marche aléatoire (par exemple avec un processus ARCH, on pourra relire sur le sujet les notes de cours Predictability of Asset Returns, en particulier le premier paragraphe, ou en français, avec des cours de licence, par Francis Diener,  http://math.unice.fr/~diener/…). Pour ceux qui veulent moins de formalisme, une martingale (en simplifant outrageusement) c’est, comme l’explique Nicolas Poupon,

Mais ces articles posent essentiellement les bases de la théorie financière. Qu’en est-il des aspects empiriques ? Andrew Lo et Craig MacKinlay ont publié un livre, en 2001, qui recense plusieurs articles sur le sujet, A Non-Random Walk Down Wall Street. Qui remet en cause l’idée de marche aléatoire, moins celle de marginale. On pourra aussi penser à Beechey Gruen & Vickrey (2000) et leur étude The Efficient Markets Hypothesis: A Survey, en ligne sur la Reserve Bank of Australia. Enfin, Pour une méthodologie plus proche de celle évoquée dans l’article en ligne sur arxiv, on pourra lire Martingales, the Efficient Market Hypothesis, and Spurious Stylized Facts de Joseph McCauley, Kevin Bassler et Gemunu Gunaratne (si on souhaite aussi utiliser des propriétés de mémoire longue). Moralité ? Non, ce n’est pas nouveau que l’on teste l’hypothèse d’efficience et de martingale. Et si je cite des articles académiques, il faut  ajouter que les journalistes économiques connaissent également cette littérature. Il y a 10 ans, Justin Fox posait la même question que nos chercheurs dans un article dans Fortune intitulé Is The Market Rational ? No, say the experts. But neither are you–so don’t go thinking you can outsmart it. On pourra aussi relire Efficiency and beyond (mentionné auparavant) qui revient sur la notion d’efficience des marchés. On y retrouve que les spécialistes savent tout ca: “on such ideas, and on the complex mathematics that described them [i.e martingales], was founded the Wall Street profession of financial engineering“. Enfin, en 2009, Richard Thaler faisait quelques rappels sur la difficulté d’interpréter correctement l’hypothèse d’efficience des marchés dans Markets can be wrong and the price is not always right. Pour conclure sur l’histoire des martingales en finance, en 2010, Fama répondait d’ailleurs de manière très claire à la question dans une interview – interview with Eugene Fama – lorsque John Cassidy lui demandait “the fundamental insight of the efficient market hypothesis [is] that you can’t beat the market ?“, et qu’il répondait “Right—that’s the practical insight. No matter what research gets done, that one always looks good“.

Pour revenir sur l’article de Pierre Barthélémy, je trouve qu’écrire “le CAC 40 est une loterie, ce que les spécialistes ne veulent pas admettre” est incroyablement méprisant envers des centaines, pour ne pas dire des milliers de chercheurs en mathématiques financières, et en économétrie de la finance. Mais rassure toi Pierre, je vais continuer à lire tes chroniques. C’est juste que ma grand mère est encore abonnée au Monde, et j’aimerais bien qu’elle continue de croire que je peux sauver l’humanité (et non pas que je refuse d’admettre ce que des chercheurs italiens viennent d’établir, même si cela confirmerait – à ses yeux – que je peux etre têtu comme une mule quand je veux).

Maintenant, pour conclure sur l’article mentionné, et sur l’analyse faite, on ne peut pas dire que les conclusions soient révolutionnaires. Disons, pour quiconque ayant lu un peu de littérature économique dans les 50 dernières années. En fait, pour être honnête, je pense qu’on peut dire que Alessandro Pluchino et Andrea Rapisarda se sont amusé (on inclura Biondo et Helbing). Sans plus de prétention. Pour ceux qui ne se souviennent pas, ils avaient eu un IGNobel en 2011, pour leur étude sur organizations would become more efficient if they promoted people at random. J’avais aussi beaucoup aimé leur article sur les aspects computationnels du principe de Peter, en ligne sur http://arxiv.org/0907.0455. Quand je lis leurs études, j’ai l’impression qu’ils s’amusent. Et le terme n’est absolument pas péjoratif dans ma bouche, loin de là. C’est un peu ce que disait Lucile Quillet la semaine passée dans http://etudiant.lefigaro.fr/…  lorsqu’en mentionnant une étude de Baptiste Coulmont, elle écrivait “c’est le jeu auquel c’est amusé Baptiste Coulmont sur son blog“. Je ne sais pas si le terme était voulu, si c’était basé sur le fait qu’un blog, c’est fait pour jouer… mais j’ai toujours revendiqué que je faisais de la recherche pour mon plaisir, pour m’amuser. Et comme me le rappelait Baptiste, le mot schola signifie école, en latin, et son étymologie grecque est le mot σχολή, signifiant loisir.

Diary of an addict

After four days offline (at least off my blog, see the previous post for more details), I have to face the truth: I am a computer addict. For sure. Here is the diary of the last four days, that were supposed to be without touching my computer, at work and at home. I tried to keep tracks of everytime I had to go on my computer. At home, that was fine (I have decided a few weeks ago that I should not check my email at home, in the evening, and in the morning, so unless I want to read the news, check on Twitter what’s going on, or write a post on my blog, I do not usually spend much time on our computer). But at the office, that’s another story…

  • Tuesday, April 2nd

6:17 Wake up, first day of the challenge.

8:17 Time to check if the code used for datascraping (on some websites) did run… It’s not like using my computer. It was just fixing problems for a future research. One minute, just cheking. Well, there was a problem in my code, I have to fix it (it takes much longer than I thought) and then, I ran it again (in order to exact some figures out of almost 200,000 internet pages). One code has been scaping a website for almost a week (and it looks like I have only one third of the data), and the other one, I have to run it everyday, to backup some daily figures from several websites.

8:32 While going to get a coffee, JF shows me the web site of the Antartica Journal of Mathematics (where he was kindly invited to submit an article, and also to apply if he’s willing to join the board) and we go through Rob Hyndman’s warning on his blog, about junk journals. No offense to the webdesigner hired by this journal, but we had a lot of fun on that website… so… nineties.

8:47 I have to go online, on my email account to download a paper I have to review. I had a reminder by an editor this weekend. Print the article. While seeking for the email, I notice that I did receive during the night an invitation to give my opinion about another article, for another journal. Just go briefly through the paper. Even if my (personal) quota for 2013 is already exceeded, I decided to accept to write a referee report. Also check something with a co-author. Quickly. Also read four comments on the blog submited during the night, and approve all of them. Damned! one of them mentioned a preprint related so something discussed in a post. Try to avoid to read the preprint. Go offline. Definitively. Turn on the music. Sacred music, namely In Seculum Longum.

9:39 After alsmost 30 minutes reading articles, I give up. I open LaTeX to type corrections for a couple of chapters I have to write down (and send them as soon as possible, the deadline was last week). So far, no internet. I take a gum, I feel nervous…

9:44 Go online to check how to use enquote{} in latex (it looks like I have some damned babel problems !) Go on several forums. Spent the morning working on my LaTeX files.

12:54 Back from lunch, I turn off the sacred music, swith to Sandinista, by The Clash,

13:14 JF want to check with me the schedule for my graduate courses for Winter 2014 session, one on extreme values and copulas, and one on time series. Have to go on the website of the university, to find what has been planed. Spend also some time seeking for old emails, since the information is only partially online.

14:02 Have to go online, one more time, to buy a ticket for the French railway, for a colleague of mine.

14:17 Mathieu asks me go check the emails received this morning, and to send an email in order to confirm that I will join a meeting about the organisation of a workshop.

14:40 Discussion with Mathieu about resubmiting an article. We need to go on my email account to check several (interesting) comments from the referees.

15:03 This time, I have to launch the internet browser. I have to find out how to avoid fist names in a bibliography, using a bib file, and get only initials. Some more time on LaTeX forums.

15:08 Remember that I have to book a room for a colleague, visiting me this August. Need to check the price on the webiste, and send emails.

15:19 Wrote a recommendation for a former student of mine. Send it.

15:26 Discussion by email with two co-authors, since we did plan to work on a joint paper this month.

15:48 Time to work on slides for a conference at the end of this month, and write codes to produce some graphs. So, still on my computer. But not (really) online.

Finally, I went back home early, to cook, and spend some time with the kids.

  • Wednesday, April 3rd

Second day.

8:17 Check again my R code. Run again one of them. Looks like there was a problem.

8:50 Go to give my course. Until noon, I spent the morning producing codes, to show how to compute chain ladder estimates, and explain roots of bootstrap in regression models.

13:03 Buzy red light on my phone when I get back to my office: it looks like some people at the faculty tried to reach me. Usually I do not answer the phone. I hate my phone: you want to reach me, send me an email. OK, in my mailbox, there are a couple of emails from the faculty “please, call us back, regarding the conference you plan to organize”. Damned, how can I tell them I am on a mission, that I try to avoid using my email (and that I have to deal with my phone-phobia at the same time).

13:16 Go online to check code for an R package I want to use to produce nice graphs

13:30 Check quickly my emails, delete 90% of them, answer to one, postpone others. Decide to go to work in a coffee shop, the whole afternoon. While taking the elevator, I started to discuss with a colleague. Looks like I missed an email about a meeting that will take place in the afternoon. I want to work on my chapters. I go to the coffee shop. Nothing serious in the meeting, as far as I understood. I spent the afternoon reading, checking typos in chapters that should be sent soon, and articles on advanced methodes in finance, based on trees… No computer for a few hours ! I did it !

17:49 I have to go online at home, time to pay my HydroQuébec bill. Damned !

17:56 While I am about to log off, I received an email from a former student of mine, with a link to a nice article (entitled “l’informatique, ça s’apprend“) I really want to share it. But I can’t. Less than 36 hours after my mission, that would be a defeat ! Do not go on Twitter !

18:02 Cooking for the kids. Will miss the Montreal Hackathon organized by R users. Wednesdy evening is not a great day to join those social events (could geek meetings be called social ?). My wife has a late course in the evening, and I still try to see how deeply addicted I can be. I finally decided that I will check my emails twice a day. But just to remove spams (or messages like “I am a student in a engineering program from India and I would like to start a PhD with you“, or “the back door in one of the building will be locked during the week-end from 22:15 till 23:42 for security reasons“) and to see if there might be important ones. Probably have to answer some of them, but I’ll try to postpone for most of them.

18:41 started a game of kid audio with the girls. Almost became a pugilat when they started to argue about kettledrum and snare-drum, I wanted to show them on youtube, but finally I gave up (there is an old saying: never interfere in a girl-and-girl fight). Decided to ask the elder to read a story to the youngest, while I was washing the dishes (which is usually the perfect timing for a dvd). Meanwhile, it looks like my son went online, for his music assignment: his teacher is using online videos to help them practicing. Argggg.

  • Thursday, April 4th

Third day.

8:29 As usual, checking the code for datascraping… reload the one that did crash (again) during the night

Error in substr(html, 8, 12) :
  invalid multibyte input string at '<e9>lair,'

Damned. Moving around 200,000 pages without being caught is difficult. Have to play some music. Gonzales, piano solo.

8:47 Have to check quickly my emails. The problem is that, on average, per week day, I have a bit more than 100 emails (excluding official spam). If I do now scan them, I end up with on thousand emails very quickly… Need to moderate a comment on the blog.

8:50 Still online, checking my emails, bad news about fundings for a student of mine, have to send a couple of emails to find a backup solution.

8:54 Quick discussion by email about copyrights for a chapter in a book

8:54 Have to send also emails to book a room for a colleague who will visit me in August

8:55 Postpone a Skype discussion with a co-author, still trying to avoid unecessary use of my computer.

8:57 Answer an email to schudele a meeting because a student asked to get a grade revision, and an adhoc committee is necessary. Looks like I am part of it.

9:17 Start to write recommendation letters for a Christophe who’s applying for positions in several universities, in France.

9:34 Back on the slides and the R code. On my computer, but offline.

10:11 Email, brief answer to a former student of mine, who might be interesting to share some datasets, but it looks like there might some confidentiality issues. I wanted to work on those data with a student, in September. Have to find an alternative.

10:45 Discussion with student in master program (face to face this time). Have to go on Dropbox to download a pdf file he wrote, and to download a couple of paper to check the proofs.

13:05 Work with Amadou, my phd student, need to find a pdf version of a book, since the property is clearly given there, but the book is out of print (no way to get a paper copy). Also go online to find a reference on a complicated model.

14:40 Email from Fred, about a reserving technique that seems new, from a paper he just discovered.

15:03 Upload on slideshare some (old) slides that do the same thing that this new paper (aren’t there anyone checking before publishing papers that results are really new ?)

15:14 Play Rodrigo y Gabriela, need something punchy to finish my day

16:56 Received an email from the immigration department, and I have to go to their website to find a doctor for a medical examination of the whole family.

19:26 Request by email from financial services to get the exact amount (in Euros) taken from my credit card. Have to go on my bank account, online.

19:43 Disccussion by email with Frédéric about one year uncertainty and bootstrap with overdispersed Poisson models

  • Friday, April 5th

Last day of the test. Fourth day.

06:18 My daughter wakes up and tells me it is unfair to have snow for her birthday. Have to go on http://meteomedia.com/weather/… to check for weather forecasts. Hopefully, we should have nice weather this afternoon…

08:16 Once again, checking R codes, still running this time ! Great ! Time to play some music.  Air, Premiers Symptômes.

8:38 Email regarding next week jury, checking legal aspects

9:04 I have to print bank informations, that I did download yesterday, check orders placed on Amazon in the last 3 months, scan documents, send them to financial services.

9:13 Have to check for an account number with financial services

9:15 Brief email to some contributors of a book that I should edit this year.

9:21 Lauch Skype, have to discuss with coauthors in France.

11:54 Received an email about courses in September, looks like a quick answer is needed.

13:10 Update the syllabus for the course I will give in September. Decided to write that cellphone will not be allowed during the class.

14:12 Work with Ben, a master student, on a paper. Need to scan notes I have written down to send him by email (this time, we did work together without using my computer).

14:38 Finished my recommendation letters for Christophe (for positions in France, more than a dozen recommendations). Have to send them individualy by email. Have also to find some email adresses that I do not have.

14:49 Received an email claiming that I cannot give my (graduate) course on extreme value in 2014, not in the official programm. Have to spend some time checking why. It seems that it has been removed from the list, and that the code has changed (hopefully JF was online to check that information much more efficiently than I would).

14:59 Received an answer from one of the colleagues I just sent recommendation to. We used to be students at the same time, a few years ago. Write a short email back to give (personal) news.

16:19 Started to type sketches of what can be the final exam for my course on GLM for actuarial science. So far 5 questions. Need to find about 40…. Will take a while.

Mission aborted.

I finally left the office later on, to pick up my son and bring him to his fencing course (to prepare for the Jeux de Montréal that will take place tomorrow). Went back home, then, for a cake with candles for my daughter’s birthday. Later on, had to spend some time online, on the blog, for my students. And started to type this post. So, here was the story of the past four days. I have to admit that looking back at those four days is quite informative :

  1. I cannot work without a computer, I can hardly work offline. No only for my research, to get help from forums on R and LaTeX, or to seek articles (the time I spent in Paris at the library making photocopies of old articles is clearly over).
  2. I do not only need a computer to write R codes, and produce LaTeX slides and artices. I need a computer… for everything. To plan meetings, for social interactions with colleagues, to find some help, to find a theorem, to book an hotel, etc.
  3. I understand more clearly why I am so unproductive in terms of research ! it looks like I spend (I was about to say waste) a lot of time on administrative tasks. Small tasks. But adding a lot of small tasks, finally, it is difficult to have 4 or 5 consecutive hours to work exclusively on some research project.

Bootstrap et régression

Lors du dernier cours, on a évoqué l’utilisation du bootstrap pour obtenir des intervalles de confiance sur des prévisions. Je mets en ligne les codes tapés en cours (très sommairement commentés, je peux renvoyer vers des vieux billets du cours ACT6420 pour des compléments). On va travailler sur ma base préférée pour évoquer la régression linéaire (avant de parler triangles de provisionnement, revenons cinq minutes sur des choses simples).

> plot(cars)
> reg=lm(dist~speed,data=cars)
> abline(reg,col="red")
> n=nrow(cars)
> x=21
> points(x,predict(reg,newdata= data.frame(speed=x)),pch=19,col="red")

On cherche ici à faire une prédiction en un point. Comme rappelé en cours (mais aussi dans le cours de modèles de prévision), quand on veut donner un intervalle de confiance pour la prévision, il convient de distinger l’intervalle de confiance pour le prédicteur (qui va dépendre de l’erreur d’estimation des paramétres) et l’intervalle de confiance pour une potentielle valeur (on peut parler de génération de scénarios, qui va dépendre en plus de l’erreur de modèle, c’est à dire de la dispersion des résidus). Commençons par l’intervalle de confiance sur la prédiction, sur le best estimate comme on dit en provisionnement

> Yx=rep(NA,500)
> B=matrix(NA,500,2)
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ base=cars[indice,]
+ #points(base,pch=3)
+ reg=lm(dist~speed,data=base)
+ abline(reg,col="light blue")
+ points(x,predict(reg,newdata=data.frame(speed=x)),pch=19,col="blue")
+ Yx[s]=predict(reg,newdata=data.frame(speed=x))
+ B[s,]=coefficients(reg)
+ }

Les valeurs bleues sont ici des prévisions possibles, obtenues en rééchantillonnant dans notre base d’observations. Pour rappel, l’intervalle de confiance (à 90%), sous hypothèse de normalité des résidus (et donc des estimateurs de la pente et de la constante de la droite de régression) s’obtient de la manière suivante

> reg=lm(dist~speed,data=cars)
> U=predict(reg,interval ="confidence",
+ level=.9,newdata=
+ data.frame(speed=0:30))
> lines(0:30,U[,2],col="red",lwd=2)
> lines(0:30,U[,3],col="red",lwd=2)

On peut comparer ici la distribution des valeurs obtenues sur nos 500 jeux de données générées, et comparer les quantiles empiriques, avec les quantiles sous hypothèse de normalité,

> hist(Yx,proba=TRUE,col="light blue",border="white")
> boxplot(Yx,horizontal=TRUE,at=.07,boxwex = 0.02,add=TRUE,col="light green")
> abline(v=U[x+1,2:3],col="red",lwd=2)
> D=density(Yx)
> lines(D)
> I=which(D$x<=quantile(Yx,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
> I=which(D$x>=quantile(Yx,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

On peut noter que les ordres de grandeur sont comparables.

> reg=lm(dist~speed,data=cars)
> quantile(Yx,c(.05,.95))
      5%      95% 
58.63689 70.31281 
> predict(reg,interval ="confidence",
+ level=.9,newdata=data.frame(speed=x)) 
       fit      lwr      upr
1 65.00149 59.65934 70.34364

Regardons maintenant l’autre type d’intervalle de confiance, sur la valeur possible de la variable d’intérêt. Cette fois, en plus de tirer des nouveaux échantillons et calculer des prédictions, on va en plus rajouter un bruit à chaque tirage, qui permettra d’obtenir une valeur possible.

> Yx=rep(NA,500)
> for(s in 1:500){
+ indice=sample(1:n,size=n,
+ replace=TRUE)
+ base=cars[indice,]
+ #points(base,pch=3)
+ reg=lm(dist~speed,data=base)
+ erreur=residuals(reg)
+ #abline(reg,lty=2)
+ E=sample(erreur,size=1)
+ Yx[s]=predict(reg,newdata=data.frame(speed=x))+E
+ points(x,Yx[s],pch=19,col="red")
+ }

Là encore, on peut comparer (graphiquement pour commencer) les valeurs obtenues par rééchantillonnage, et celle obtenues sous hypothèse de normalité,

> hist(Yx,proba=TRUE,col="light blue",border="white")
> boxplot(Yx,horizontal=TRUE,at=.025,boxwex = 0.005,add=TRUE,col="light green")
> abline(v=U[2:3],col="red",lwd=2)
> D=density(Yx)
> lines(D)
> I=which(D$x<=quantile(Yx,.05))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)
> I=which(D$x>=quantile(Yx,.95))
> polygon(c(D$x[I],rev(D$x[I])),c(D$y[I],rep(0,length(I))),col="blue",border=NA)

Ce qui donne, numériquement, les comparaisons suivantes

> quantile(Yx,c(.05,.95))
      5%      95% 
44.43468 96.01357 
> (U=predict(reg,interval ="prediction",level=.9,newdata=data.frame(speed=x)))
       fit      lwr      upr
1 67.63136 45.16967 90.09305

On observe cette fois une légère asymméytrie vers la droite. Manifestement, on ne peut pas supposer les résidus Gaussien, car il y a plus de grandes valeurs positives, que négatives. Ce qui fait du sens compte tenu de la nature des données (une distance ne peut être négative).

On avait ensuite commencé à discuter de l’utilisation des modèles de régression en provisionnement. Afin d’avoir des données présentant de l’indépendance, on avait rappelé qu’il fallait travailler avec les incréments de paiments, et non pas les paiements cumulés.

> T
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730   NA
[3,] 3871 5345 5398 5420   NA   NA
[4,] 4239 5917 6020   NA   NA   NA
[5,] 4929 6794   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA
> n=ncol(T)
> Y=T
> Y[,2:n]=T[,2:n]-
+         T[,1:(n-1)]
> Y
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   53   22   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

On peut alors constituer une base de données, avec comme variables explicatives la ligne et la colonne.

> y=as.vector(as.matrix(Y))
> base=data.frame(
+ y,
+ ai=rep(2000:2005,n),
+ bj=rep(0:(n-1),each=n))
> 
> head(base,12)
      y   ai bj
1  3209 2000  0
2  3367 2001  0
3  3871 2002  0
4  4239 2003  0
5  4929 2004  0
6  5217 2005  0
7  1163 2000  1
8  1292 2001  1
9  1474 2002  1
10 1678 2003  1
11 1865 2004  1
12   NA 2005  1
> tail(base,12)
    y   ai bj
25  7 2000  4
26 10 2001  4
27 NA 2002  4
28 NA 2003  4
29 NA 2004  4
30 NA 2005  4
31 21 2000  5
32 NA 2001  5
33 NA 2002  5
34 NA 2003  5
35 NA 2004  5
36 NA 2005  5

On peut alors commencer par utiliser le modèle Regression models based on log-incremental payments de Stavros Christofides, basé sur une modélisation lognormale, introduite initialement par Etienne de Vylder en 1978 (Markus en parle, en trois parties, sur son blog http://lamages.blogspot.ca/Barnett%20Zehnwirth)

> reg1=lm(log(y)~
+ as.factor(ai)+
+ as.factor(bj),data=base)
> summary(reg1)

Call:
lm(formula = log(y) ~ as.factor(ai) + as.factor(bj), data = base)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26374 -0.05681  0.00000  0.04419  0.33014 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         7.9471     0.1101  72.188 6.35e-15 ***
as.factor(ai)2001   0.1604     0.1109   1.447  0.17849    
as.factor(ai)2002   0.2718     0.1208   2.250  0.04819 *  
as.factor(ai)2003   0.5904     0.1342   4.399  0.00134 ** 
as.factor(ai)2004   0.5535     0.1562   3.543  0.00533 ** 
as.factor(ai)2005   0.6126     0.2070   2.959  0.01431 *  
as.factor(bj)1     -0.9674     0.1109  -8.726 5.46e-06 ***
as.factor(bj)2     -4.2329     0.1208 -35.038 8.50e-12 ***
as.factor(bj)3     -5.0571     0.1342 -37.684 4.13e-12 ***
as.factor(bj)4     -5.9031     0.1562 -37.783 4.02e-12 ***
as.factor(bj)5     -4.9026     0.2070 -23.685 4.08e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1753 on 10 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared: 0.9975,	Adjusted R-squared: 0.9949 
F-statistic: 391.7 on 10 and 10 DF,  p-value: 1.338e-11 

> base$py=exp(predict(reg1,
+ newdata=base)+summary(reg1)$sigma^2/2)
> round(matrix(base$py,n,n),1)
       [,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 2871.2 1091.3 41.7 18.3  7.8 21.3
[2,] 3370.8 1281.2 48.9 21.5  9.2 25.0
[3,] 3768.0 1432.1 54.7 24.0 10.3 28.0
[4,] 5181.5 1969.4 75.2 33.0 14.2 38.5
[5,] 4994.1 1898.1 72.5 31.8 13.6 37.1
[6,] 5297.8 2013.6 76.9 33.7 14.5 39.3
> sum(base$py[is.na(base$y)])
[1] 2481.857

On obtient un montant un peu différent de celui obtenu par la méthode Chain Ladder, mais néanmoins comparable. On peut aussi tenter une régression de Poisson (avec un lien logarithmique), comme suggéré en 1975 par Hachemeister et Stanard,

> reg2=glm(y~
+ as.factor(ai)+
+ as.factor(bj),data=base,
+ family=poisson)
> summary(reg2)

Call:
glm(formula = y ~ as.factor(ai) + as.factor(bj), family = poisson, 
    data = base)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3426  -0.4996   0.0000   0.2770   3.9355  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)        8.05697    0.01551 519.426  < 2e-16 ***
as.factor(ai)2001  0.06440    0.02090   3.081  0.00206 ** 
as.factor(ai)2002  0.20242    0.02025   9.995  < 2e-16 ***
as.factor(ai)2003  0.31175    0.01980  15.744  < 2e-16 ***
as.factor(ai)2004  0.44407    0.01933  22.971  < 2e-16 ***
as.factor(ai)2005  0.50271    0.02079  24.179  < 2e-16 ***
as.factor(bj)1    -0.96513    0.01359 -70.994  < 2e-16 ***
as.factor(bj)2    -4.14853    0.06613 -62.729  < 2e-16 ***
as.factor(bj)3    -5.10499    0.12632 -40.413  < 2e-16 ***
as.factor(bj)4    -5.94962    0.24279 -24.505  < 2e-16 ***
as.factor(bj)5    -5.01244    0.21877 -22.912  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 46695.269  on 20  degrees of freedom
Residual deviance:    30.214  on 10  degrees of freedom
  (15 observations deleted due to missingness)
AIC: 209.52

Number of Fisher Scoring iterations: 4

> base$py2=predict(reg2,
+ newdata=base,type="response")
> 
> round(matrix(base$py2,n,n),1)
       [,1]   [,2] [,3] [,4] [,5] [,6]
[1,] 3155.7 1202.1 49.8 19.1  8.2 21.0
[2,] 3365.6 1282.1 53.1 20.4  8.8 22.4
[3,] 3863.7 1471.8 61.0 23.4 10.1 25.7
[4,] 4310.1 1641.9 68.0 26.1 11.2 28.7
[5,] 4919.9 1874.1 77.7 29.8 12.8 32.7
[6,] 5217.0 1987.3 82.4 31.6 13.6 34.7
> 
> sum(base$py2[is.na(base$y)])
[1] 2426.985

La prédiction coïncide avec l’estimateur obtenu par la méthode Chain Ladder. Le lien avec les méthodes de biais minimal a été établi par  Klaus Schmidt et Angela Wünsche en 1998, dans Chain ladder, marginal sum and maximum likelihood estimation. La semaine prochaine, on parlera des méthodes de bootstrap pour obtenir des intervalles de confiance, ou des quantiles, sur les montants de réserve. Je ne sais pas si j’aurais le temps de taper des transparents, je préfère, sur cette partie du cours taper au fur et à mesure, et écrire au tableau. Je renvoie au Chapitre 3 du livre avec Christophe Dutang – en ligne sur http://cran.r-project.org/doc/contrib/ – pour le détail. C’est le code que je tape en cours, tout en essayant de répondre aussi aux questions.

Off for a couple of days

This morning, the elder kids start at school a program called ten days with screens off. Since we do not have a TV, they were wondering what would it change for them. Because they start to spend some time on our computer at home, they’re not supposed to use it for the next ten days. And to be consistent, I have decided to follow their tracks. The problem is that I do not think I can survive at work working offline for ten days. But I want to try. At least, no more computer for the next two days (except Wednesday morning since I have my class). If possible three. No emails (can we survive in the academic world without an internet connection for two days?), and no more posts on the blog (and activities related to the comments) for the next days…  I’ll be back in … a few days (if I do not – finally – find my blogging activity useless and meaningless).

“I am sorry to inform you…”

I just received another decision letter for a paper submitted.But this time, it was fast. In less than 24 hours, I received a notification by email (the name of the journal is hidden with xxxxx’s for obvious reasons):

Dear Professor Charpentier,

Thank you for submitting your paper to xxxxx.

In order to provide a quick reply to your submission, an automatic process, has been used according. The process is based on proximity measures between the article you submitted, and recent articles published in our journal.

1) The correlation between your keys words and standard keywords in xxxxx is rather high (0.926).

2) 66% of your keywords do match with the top 10 most cited papers in 2012 in xxxxx, so your paper has 90% chance of getting at least 4 citations in less than a year and at least 17 citations in less than 5 years.

3) The rate of references cited in your article that were published in xxxxx is high enough (37%) to consider a publication in xxxxx

The academic and citation based score is 0.857 which makes the paper suitable for publication. Nevertheless, before publishing the paper, major revisions are necessary,

1) Your style does not fit well with our journal: your sentenses have (on average) 13.3 words, while the standard in this journal is 16.2. We recommend that you use longer sentences or delete some periods in your text.

2) The words you use are not close enough to the scope of the journal: for example, estimator appears only twice  (the average is 4.34) and methodology appears only once (average 3.67).

3) The affiliation is not appropriate: this year, 12% of the articles received were from authors in Canada. In order to reinforce positive discrimination, it would be better if your affiliation were from the Eastern Europe (currently only 4% of the articles we receive are from Eastern Europe). We can provide you with a list of potential universities if necessary.

4) Your name is not valid for xxxxx:

  • For a single authored paper, the most common first letter in the last name of the first author is H not C and your name which contains 11 letters is too long (the average is 7.91 letters). Lastly, more than 62% authors have a middle name, we suggest that you add one,
  • For a paper with two authors, the most common first letter of the first author is D, hence C is fine, but your name is a bit too long since the average number of letters in the last name for the first author is 9.57.

Thus, we recommend either that you change your name (we can provide you with a list if necessary) or that you find another author to add on your paper and delete the last two letters of your last name.

For the reasons listed above, we have to ask for major revisions. When the paper is resubmitted, it will be send to a (human) referee for proofread,

Please note that we receive almost 1,000 submissions a year; our target acceptance rate is 10% and many good papers get rejected.

Thank you for giving us the opportunity to consider your work.

Sincerely,

(etc)

Damned ! That was close ! Let me find a co-author to add, and I resubmit it !