Somewhere else, part 50

A lot of interesting posts and articles, this week, here and there

avec comme toutes les semaines, quelques billets et articles en français

Did I miss something ?

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…

Somewhere else, part 49

Some extremely interesting posts, this week, again on the Reinhart-Rogoffing story (I do mention many posts and articles related to that story, because I think it is extremely informative about how research in Economics actually works…)

but (probably) more interesting, a series of posts on Excel,

And a usual, a lost of interesting posts and articles, on other topics, e.g.

See also

Avec comme souvent plusieurs billets en français

Did I miss something ?

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…

Somewhere else, part 48

Following previous posts on this blog (# 46 and 47), a couple of articles that are worth reading,

and, as usual, several interesting posts, here and there

et comme toujours, une poignée d’articles en français

Did I miss something ?

Somewhere else, part 47

Still a lot of posts, those past days, about Reinhart-Rogoffing

“Reverse-causality” on felixsalmon ‘s blog http://blogs.reuters.com/felix-salmon/ see

and to conclude (?) http://bloomberg.com/news/…

Of course, there are other posts and articles, on other topics,

Quelques articles intéressants en français,

Did I miss something ?

Somewhere else, part 46

Interesting econometric posts, here and there, during the past two days, about Kenneth Rogoff and Carmen Reinhart’s paper. All those posts were very informative, about economics and politics, econometrics and datasets, and (publicly) admitting mistakes,

But there were also interesting posts, here and there, on other topics,

avec comme souvent des articles en français,

Did I miss something ?

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/…

Somewhere else, part 45

Interesting readings from this week-end, starting with the following two posts :

and a lot of other posts and articles, outside this blog

Comme souvent, quelques articles en français, à commencer par un joli poisson d’avril

et une image pour sourire

avec plusieurs lectures intéressantes

Did I miss something ?

Somewhere else, part 44

Some interesting posts and articles this week,

a very nice tweet

  • People found guilty of not using punctuation deserve the longest sentence possible.” by sixthformpoet

and the only picture I want to share this week

et comme à chaque fois, quelques billets et articles en français,

Did I miss something ?

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