Données Personnelles et Assurance

Depuis plusieurs semaines (pour ne pas dire plusieurs mois), je discute avec Delphine Cocteau-Senn et Rodolphe Bigot (dialogue de juristes avec un actuaire) sur la protection des données personnelles en assurance. Tout avait commencé il y 15 mois, avec une table ronde organisée à Amiens, avec un resto la veille, les discussions pendant la table ronde, puis des dizaines et des dizaines d’échanges de messages, de coup de fils, de restos…. Rodolphe Bigot m’avait sollicité pour une autre table ronde sur le droit de l’assurance à Caen, et de mon côté, j’avais demandé à Delphine Cocteau-Senn d’intervenir dans la formation Actuariat – Data Science de l’Institut des Actuaires, depuis l’an dernier. Bref, on a eu de très nombreux échanges, et quand on voit les prises de positions ici et là (pour ne pas parler de brassage d’air), on s’était dit que nos discussions pourraient éclairer les débats qui agitent actuaires et juristes.

C’était d’autant plus intéressant, je pense, qu’on a vraiment appris, mutuellement, de ces échanges. Echanges entre un actuaire et des juristes, tous les trois universitaires. Depuis plusieurs semaines, on tente de mettre par écrit ces échanges, et le résultat est aujourd’hui en ligne. La forme retenue est une présentation juridique de plusieurs points (en lien avec des notions évoquées dans le Réglement Général sur la Protection des Données – n° 2016/679), quelques questions actuarielles, et une ébauche de réponse. On a longuement hésité sur la forme, mais on a gardé la forme de la discussion, probablement plus agréable à lire qu’un texte écrit à trois. Je mets ici l’introduction (le document complet est en ligne)

Continue reading Données Personnelles et Assurance

Actuarial Pricing Game

Within a few days, we will launch the fourth actuarial pricing game.

  1. People have to register for that new game. People have to fill a form online before February 28th to get a training dataset. As in our previous versions, players can be individuals or in a team (identified as one player).
  2. When registration closes, players will receive a training dataset (training.csv). The training dataset will be based on two years of data on household policies. There will be information about the insurance policy from underwriters, as well as, claims data. They will also receive a pricing dataset (pricing.csv) with only underwriters information and no claims. The goal is to propose a premium for all records in the pricing dataset. Note that those registered will also receive additional information and template files after registration closes.
  3. Before  April 9th, 2018, players will have to provide two files: a pricing function (either a formula.py Python function, or a formula.R R function) and a dataset with prices prices.csv). The pricing function should be an understandable explicit function which, once applied to the training dataset, will yield premium prices for each record. There will be limits on the complexity of this function. To ensure this, players will receive instructions on what can be used to construct formulas either in Python or R.
  4. Before May 14th, 2018, players are asked again to submit two files: a dataset with prices (second_prices.csv) and a pricing model (either a model.py Python file, or a model.R R file). For this round there are no restrictions on what can be used in the model, as long as, by using the model file submitted the prices dataset (second_prices.csv) can be reproduced. Players will again receive further instruction on this after registration.

Unicode, UTF-8… la base de l’encodage

Un petit billet rapide, sans aucune prétention. Juste pour revenir sur des trucs bizarres sur l’encodage… J’avais commencé à m’y intéresser il y a 3 ans, alors que je discutais avec des utilisateurs de R en Inde et au Japon. Mais revenons à la base… Si on regarde la page “unicode” de Wikipedia (ou plus précisément, si on regarde le code source, on peut avoir des informations (dans les métadonnées) : en l’occurrence, on voit que la page est encodée en UTF-8

C’est cohérent car de plus en plus de documents en ligne semblent encodé dans ce format. Si on revient à la source, l’encodage c’est lié au passage par les bytes 0/1. Plein de sites expliquent (en gros) la philosophie. Quand j’étais petit, on nous expliquait les caractères ASCII, qui utilisent 7 bits ce qui permet de disposer  128 caractères uniquement, numérotés de 0 à 127.  Genre “A” était le 65ème caractère, et “a” le 97ème. Il y avait des symboles exotiques (comme “@” qui est le 64ème) ou la ponctuation, mais pas d’accents, comme “à”. Si on veut les lettres avec des accents, il faut plus que 128 caractères. Est-alors arrivée ISO/CEI 8859, qui proposait d’encoder les caractères sur 8 bits (et pas 7). On parle aussi de ISO 8859-1 plus populaire sous le nom “latin-1” très populaires dans les pays latins, avec le “à“, et puis on a continué d’enrichir pour arriver au ISO 8859-15 (latin-9) avec dedans le symbole “€”. Pour les langues asiatiques, il a fallu enrichir encore. Bref, c’est un jeu sans fin. Dans les années 90, est arrivée la norme ISO 10646 correspondant au “jeu universel de caractères” ou UCS pour Universal Character Set (en anglais), qui donnera Unicode. Je crois que les puristes font la différence entre les jeux de caractères (l’ensemble de caractères auxquels on attribue à chacun un point de code unique, la liste des caractères ASCII par exemple) et l’encodage (qui est la façon de représenter ce point de code en mémoire). Mais on va faire simple…

Là où ça se complique c’est que les OS utilisaient (par défaut) des encodages différents. Sous Windows, UTF-16 et Unicode sont majoritairement utilisé. Sous Mac OS,il y avait un truc appelé MacRoman mais il semble que UTF-8 l’emporte aujourd’hui. Enfin, sous linux, historiquement c’était du latin-1 par défaut, mais progressivement l’UTF-8 semble prendre le dessus. Bref, tous ceux qui naviguent entre les OS le savent : les accents c’est l’enfer. Et le principal soucis c’est qu’il est impossible, quand on récupère un fichier d’avoir son encodage a priori. Regardons le fichier suivant BaseANSI.txt

On peut tenter de regarder s’il y a des sortes de méta-données dans le document…

mais ce n’est pas le cas. Le document est un document brut, avec juste les données, et aucune métadonnées. Rien sur l’encodage (en particulier, car c’est ce qui nous intéresse aujourd’hui). Pour l’anecdote, si au lieu d’ouvrir le fichier avec mon bloc-note sous Windows, on tente sur un mac, on a

Et quand je passe par le terminal DOS sous Windows, c’est la même chose (ce qui est un peu rassurant)

Regardons un autre fichier BaseUNICODE.txt. Si vous voyez une différence avec l’autre, dites le moi…

Pourtant il y en a une : la preuve, dans un terminal, on peut lire les accents du fichier,

Un dernier ? Regardons BaseUTF8.txt,Là encore, on a l’impression d’avoir le même fichier. Et pourtant non.

Commençons par récupérer les trois fichiers (via R)

download.file(url="http://freakonometrics.free.fr/BaseUTF8.txt",
destfile="BaseUTF8.txt")
download.file(url="http://freakonometrics.free.fr/BaseANSI.txt",
destfile="BaseANSI.txt")
download.file(url="http://freakonometrics.free.fr/BaseUNICODE.txt",
destfile="BaseUNICODE.txt")

Techniquement, ces documents sont différents, même si dans le bloc note, on ne voyait rien : il n’ont pas la même taille !

file.size("BaseUNICODE.txt")
[1] 322
file.size("BaseANSI.txt")
[1] 163
file.size("BaseUTF8.txt")
[1] 171

On peut essayer de les lire, pour voir…

read.table("BaseANSI.txt",header=TRUE,sep=";")
                       Département No Habitants
1                              Ain 01    631877
2                          Ardèche 07    324209
3                 Bouches-du-Rhône 13   2016622
4                     Corse-du-Sud 2A    152730
5 Côte-dOr;21;533147\nCôtes-dArmor 22    598357

On a un petit soucis avec l’apostrophe, mais rien de grave

read.table("BaseANSI.txt",header=TRUE,sep=";",quote="\"")
       Département No Habitants
1              Ain 01    631877
2          Ardèche 07    324209
3 Bouches-du-Rhône 13   2016622
4     Corse-du-Sud 2A    152730
5        Côte-d'Or 21    533147
6    Côtes-d'Armor 22    598357

R arrive parfaitement a lire (sur mon ordinateur du moins) cette première base de données. Tentons une seconde

read.table("BaseUTF8.txt",header=TRUE,sep=";",quote="\"")
    ï..DÃ.partement No Habitants
1               Ain 01    631877
2          Ardèche 07    324209
3 Bouches-du-Rhône 13   2016622
4      Corse-du-Sud 2A    152730
5        Côte-d'Or 21    533147
6    Côtes-d'Armor 22    598357

On a manifestement un soucis, un soucis d’encodage. Le soucis c’est que, d’ordinaire, on ne sait pas comment le fichier a été encodé. Et comme on l’a vu, il n’y a pas de métadonnées qui nous donne l’encodage. Comment faire ? Heureusement, il existe une fonction sous R qui devine l’encodage. Littéralement.

library(readr)
guess_encoding("BaseUTF8.txt", n_max = 1000)
# A tibble: 3 x 2
  encoding   confidence
 
1 UTF-8           1.00 
2 ISO-8859-1      0.620
3 ISO-8859-2      0.430

Autrement dit, le best-guess est ici un encodage UTF-8. Essayons pour voir,

read.table("BaseUTF8.txt",header=TRUE,sep=";",encoding="UTF-8",quote="\"")
  X.U.FEFF.Département No Habitants
1                  Ain 01    631877
2              Ardèche 07    324209
3     Bouches-du-Rhône 13   2016622
4         Corse-du-Sud 2A    152730
5            Côte-d'Or 21    533147
6        Côtes-d'Armor 22    598357

Ça marche ! ou disons que ça marche presque. J’ai un soucis dans mes libellés de variables. Mais disons que c’est anecdotique ici.

Tentons le dernier fichier

read.table("BaseUNICODE.txt",header=TRUE,sep=";",quote="\"")
  ÿþD
1  NA
2  NA
3  NA
4  NA
5  NA
6  NA
7  NA
8  NA

Comme ça ne marche pas, essayons d’avoir un best-guess sur l’encodage

guess_encoding("BaseUNICODE.txt", n_max = 1000)
# A tibble: 3 x 2
  encoding   confidence
 
1 UTF-16LE        1.00 
2 ISO-8859-1      0.530
3 ISO-8859-2      0.270

On en tient un !

read.table("BaseUNICODE.txt",header=TRUE,encoding="UTF-16LE")
  ÿþD
1  \n
2   B
3  \n
4   C
5  \n
6   C
7  \n
8   C

Damned, ça ne marche pas… Heureusement, Ewen est venu à ma rescousse. La solution est relativement étrange

read.table("BaseUNICODE.txt",header=TRUE,fileEncoding="UTF-16LE", sep=";", quote="")
       Département No Habitants
1              Ain 01    631877
2          Ardèche 07    324209
3 Bouches-du-Rhône 13   2016622
4     Corse-du-Sud 2A    152730
5        Côte-d'Or 21    533147
6    Côtes-d'Armor 22    598357

Oui, on n’utilise pas l’encodage via “encoding” mais “fileEncoding” ! Subtil non… Cela dit, si on voulait être cohérent on devrait utiliser la fonction du package pour lire la base, non ? Car dans le package “readr” il y a une fonction “read_table”

read_table("BaseUNICODE.txt",locale=locale(encoding = "UTF-16LE"))
Error in guess_types_(datasource, tokenizer, locale, n = guess_max) : 
  Incomplete multibyte sequence

Sauf que quand on la lance, on a un message d’erreur… Bref, l’encodage c’est compliqué. C’est ce qui est expliqué dans le forum de discussion du package cela dit. Je ne suis pas sûr que mon billet serve à quoi que ce soit, mais je voulais garder une trace de tout ça !

De la communication scientifique…

Vendredi dernier, j’étais convié à une matinée de discussion autour de la communication scientifique, par des universitaires, à destination des doctorants de la région. Muriel Lefebvre a fait une présentation captivante sur la communication scientifique, et comme promis, je vais mettre par écrit quelques notes prises (je pourrais préciser que j’ai rajouté des réflexions personnelles, mélées à mes notes, car je n’ai pas réussi à déméler). L’idée était de questionner la communication scientifique, et ses pratiques diverses. Le premier point était de noter la variété des formes, la variété des objectifs (écrire permet de réfléchir, de communiquer, de valider, etc), et la variété des publics. En effet, on peut écrire pour soi, pour les pairs (présents, mais aussi passés, en revisitant des études), et pour d’autres sphères. Pour illustrer ce dernier point, Muriel Lefebvre a proposé de revenir sur la rose des vents proposée par Callon, Courtial & Penan (1993), que je reprends sommairement

(bon, le géomètre qui est en moi bout, car une rose des vents est une étoile à 32 directions, pour les 32 directions du vent, car 32=25 or ici on n’a que 5 branches…. mais peu importe). La “communication scientifique” est un objet difficile à identifier parce qu’elle est confrontée à ces 5 axes, qui sont de nature très différente.

En haut à gauche, on retrouve les deux branches classiques du monde académique : la communication liée à l’enseignement, et celle liée à la recherche académique. Dans ces deux sphères, l’enseignant-chercheur est légitime parce que c’est dans son nom : il est enseignant et chercheur. Le terme “média” est à entendre au sens large, il s’agit de toucher un grand public. La “sphère économique” est liée à d’autres termes importants aujourd’hui, comme l’innovation, ou la valorisation. Pour ce dernier, si on pense naturellement aux applications industrielles, on peut toutefois rajouter des activités commerciales au sens large, comme cela peut se voir en droit. Enfin, les “institutions” pourrait être un terme fourre-tour, dès lors que l’enseignant-chercheur est appelé à partager son expertise (par exemple dans des discussions en lien avec la bioéthique).

Il y a aussi une confusion sur les termes, où les termes de publication, de médiation, de vulgarisation, sont utilisés dans tous les sens. Historiquement, la communication scientifique passe par l’écriture, qui légitime le travail. La publication permet de diffuser des travaux, tout en servant d’instrument de validation, d’assurer la paternité d’une découverte, de valoriser un travail, et elle sert aujourd’hui de base pour évaluer la recherche. Cette publication scientifique passe par un processus éditorial et une évaluation par les pairs. La chaîne de l’écriture scientifique est un processus non linéaire, largement décrit par Hans Reichenbach et Karl Popper, mais qui historiquement, restait dans la sphère académique. La vulgarisation (du latin vulgaris) vise à transmettre un contenu (scientifique) d’une manière aisément compréhensible par un large public. Son apogée a été au XIXème siècle, et a été basé sur une relation fondamentalement hiérarchique, du savant vers l’ignorant, du sacré vers le profane. La transmission est alors unidirectionnelle. Ce terme a aujourd’hui une connotation négative, on lui préferera celui d’instruction, ou d’éducation populaire. Cette vulgarisation n’est alors pas une fin en soi, mais une méthode. Puis dans les années 60/70, on voit arriver l’idée de médiation scientifique. Il s’agit alors de mettre du lien, de dialoguer, dans une époque qui remet en question l’autorité de la science, qui n’associe plus science et progrès aussi aisément. Ce dialogue qui doit être rétabli entre science et société se retrouve plus généralement entre la culture et la société. C’est l’époque des maisons de la culture, nées de l’idée d’ouvrir et de rentre accessible l’art (et la culture) à tout le monde. Par la suite, dans les années 80, la notion de culture va s’élargir, en reconnaissant des formes plurielles et diverses. Et la science va bénéficier de ce mouvement. C’est l’idée de “(re)mettre la science en culture“. La culture scientifique est alors une sous partie de la culture commune.

Ces changements vont se traduire par des nouvelles formes de publicisation de la science, pour reprendre la terminologie de Jürgen Habermas. Il s’agit pour le scientifique de rendre public ses travaux, d’être visible. La publicisation se fait par les revues savantes traditionnelles, mais aussi les revues électroniques, les archives ouvertes (hal, arxiv, etc), les réseaux sociaux académiques (research gate, academia, etc), qui restent tous dans la sphère académique traditionnelle. Mais il y a aujourd’hui également de la publicisation vers un publique non-académique, qui peut être atteint par les carnets et blogs de recherche, ou par des chaînes youtube. Depuis 1982, et la redéfinition du rôle de l’enseignant-chercheur, l’importance de la diffusion n’a cessé d’augmenter.

Les expériences de vulgarisation/médiation ne cessent de se multiplier, allant de la fête de la science, à la nuit européenne des chercheurs ou MT180. L’objectif n’est plus très clair, même si les institutions académiques semblent valoriser ces activités, allant d’un souhait de partager des savoirs à une forme d’injonction de visibilité. Pour conclure, Muriel Lefebvre est longuement revenue sur cette double injonction, de qualité et de visibilité. Les fontières semblent se brouiller, Il y a une demande d’humanisation (voire une personnalisation) de la recherche.

On a ensuite eu une table ronde, et ma prise de notes en a pâti (on s’est un peu écouté parler). Il y a eu de nombres interrogations sur la légitimité pour faire de la médiation, et sur l’ouverture des revues savantes. Mais une intervention a aussi souligné l’émergence d’un nouvel acteur depuis quelques années dans la médiation : les associations. Le point de départ semble avoir été les associations en lien avec les personnes malades, atteintes du HIV, dans les années 80, qui ont voulu comprendre. Ces associations peuvent jouer un rôle important dans la médiation. L’histoire amusante est un article de Muriel Lefebvre, intitulé The Circulation of Scientific Articles in the Sphere of Web-Based Media: Citation Practices, Communities of Interests and Local Ties. Dans cet article, elle s’interroge sur la vie d’un article paru en 2012, dans PLoS, intitulé “Freshwater Killer Whales”: Beaching Behavior of an Alien Fish to Hunt Land Birds, et plus particulièrement qui mentionne cet article. Pour information, cet article étudie un gros poisson, la silure (Silurus glanis)

et ce poisson semble avoir tendance à manger des oiseaux, comme ces pigeons (Columbia livia). L’article illustre avec quelques photos (les plus courageux pour aussi trouver en ligne des vidéos).

La question n’était pas de savoir qui dans les revues scientifiques citent cet article (comme peut le faire la bibliométrie traditionnelle) mais plutôt de comprendre au sens large qui a pu mentionner cette étude (publiée en accès ouvert). En fait, cet article a été énormément mentionné sur des sites internet d’associations de pêcheurs (comme seriouslyfish), qui peuvent relayer une information publiée par des chercheurs, pour la diffuser auprès d’un public beaucoup plus large. En fait, le processus de la production scientifique est complexe, car ces associations peuvent aussi faire remonter un problème auprès de scientifiques, qui prendront alors le temps de creuser, et de comprendre. Il s’agit alors véritablement d’échanges. Plus généralement aujourd’hui, on retrouvera beaucoup d’association de malades, qui tentent de comprendre la recherche, mais qui tentent aussi de sensibiliser (et donc de financer). Il n’y a plus de hiérarchie, entre sachants et profanes, mais au contraire une vraie médiation jouée par le chercheur.

Seven Questions Raised by the Bootstrap

Pour le cours de magistère, suite à la première partie du cours, présentant les méthodes de simulation, faite par Magali Fromont, nous allons travailler sur le boostrap, en partant de l’article de Bradly EfronSix Questions Raised by the Bootstrap. L’article date un peu (il a été publié en 1990) mais il pose six questions qui permettront de revenir sur le fonctionnement du bootstrap (et de la statistique mathématique).

  • why do maximum likelihood estimated distributions tend to be short-tailed ?
  • why does the delta method tend to underestimate standard errors ?
  • why are cross-validation estimators so variable ?
  • what is a correct confidence interval ?
  • what is a good nonparametric pivotal quantity ?
  • what are computationally efficients ways to bootstrap ?

Je me permets d’en rajouter une septième

  • why generate samples of size n, why not try smaller sizes ?

Ce dernier point pose la question du subsampling. Je renvoie aux notes de cours de Charles Geyer qui présente ce concept. Je recommande à tous la lecture de What Teachers Should Know about the Bootstrap. Enfin, pour un ouvrage de référence plus moderne que l’article de Bradley Efron, je ne peux que recommander les chapitres 10 et 11 de Computer Age Statistical Inference.

 

Course on “Big Data for Economics”

In July, I will give a lecture 2018 edition of the Summer School at the UB School of Economics, in Barcelona. It will be a four day crash course. Registration should be opened.

Lecture 1: Introduction : Why Big Data brings New Questions
Lecture 2: Simulation Based Techniques & Bootstrap
Lecture 3: Loss Functions : from OLS to Quantile Regression
Lecture 4: Nonlinearities and Discontinuities
Lecture 5: Cross-Validation and Out-of-Sample diagnosis
Lecture 6: Variable and model selection
Lecture 7: New Tools for Classification Problems
Lecture 8: New Tools for Time Series & Forecasting

Prédire l’espérance de vie (à la naissance)

Pour finir les exemples de modèles de prévision, on a parlé de calculs d’espérance de vie (à la naissance). Le point de départ des modèles démographiques, ce sont les tables de mortalité, dont je parlais dans le précédant billet. Mais cette vision est très biaisée, car elle suppose qu’il n’y aura pas d’amélioration des conditions de vie. Pour faire les choses proprement, on utilise des données plus complètes, avec des nombres de décès en fonction de l’âge x, mais aussi de la date t.

url="http://freakonometrics.free.fr/FranceDeaths_1x1.txt"
download.file(url,"FRD.txt")
url="http://freakonometrics.free.fr/FranceExposures_1x1.txt"
download.file(url,"FRE.txt")
DECES=read.table("FRD.txt",skip = 3,header=TRUE)
EXPOS=read.table("FRE.txt",skip = 3,header=TRUE)

On notera D_{x,t} le nombre de décès observé, et E_{x,t} l’exposition. De telle sorte que la probabilité de décéder dans l’année, pour une personne d’âge x à la date t est q_{x,t}=D_{x,t}/E_{x,t}. Classiquement, ces données sont stockées dans une matrice pour la visualisation, et dans une base de données pour les régressions

annees = 1816:2015
ages   = 0:110
MDF=matrix(DECES$Female,length(ages),length(annees))
MDH=matrix(DECES$Male,length(ages),length(annees))
MEF=matrix(EXPOS$Female,length(ages),length(annees))
MEH=matrix(EXPOS$Male,length(ages),length(annees))
MEF[MEF<=0]=1
MEH[MEH<=0]=1
QF=MDF/MEF
QH=MDH/MEH
QF[QF<=0]=NA
QH[QH<=0]=NA

(il faut bricoler un peu pour éviter les problèmes de valeurs nulles car (i) on fait des ratios (ii) on passera ensuite au logarithme). On peut visualiser cette surface, en fonction de x et t.

persp(log(QF))

ou (plus visuel)

library(rgl)
persp3d(ages,annees,log(QH),col="light blue")

Pour modèliser l’évolution de q_{x,t} on peut s’inspirer du modèle de Lee & Carter (1992) qui suppose que \log(q_{x,t})=A_x+B_x\cdot K_t. Les \mathbf{A}=(A_0,A_1,\cdots,A_{110}) c’est en quelques sortes l’évolution moyenne de \log(q_{x,t}). \mathbf{K}=(K_{1816},K_{1817},\cdots,K_{2015}) va permettre de comprendre les gains de conditions de vie, la baisse de la probabilité de décéder dans l’année. Ces gains ne sont pas uniformes, et on utilise \mathbf{B}=(B_0,B_1,\cdots,B_{110}) pour faire dépendre les gains de l’âge.

Pour estimer les paramères, \mathbf{A}, \mathbf{B} et \mathbf{K}, on peut être tentés d’utiliser un modèle binomial, au sens où D_{x,t}\sim \mathcal{B}(E_{x,t},q_{x,t}), qui est le modèle de base en assurance-vie. Autrement dit, ici D_{x,t}\sim \mathcal{B}(E_{x,t},\exp[A_x+B_x\cdot K_t]).

Une autre piste est d’utiliser la loi des petits nombres, qui dit que si la probabilité est faible (et c’est le cas quand on s’intéresse à la probabilité de mourir dans l’année), la loi binomiale peut être approchée par une loi de Poisson. Autrement dit, D_{x,t}\sim \mathcal{P}(E_{x,t}\cdot q_{x,t}), i.e. D_{x,t}\sim \mathcal{P}(\cdot\exp[\log(E_{x,t})+A_x+B_x\cdot K_t]). On reconnait ici une régression de Poisson, avec comme variables explicatives l’âge x, l’année t, et l’exposition apparait comme variable offset. Le seul soucis est que ce n’est pas une régression linéaire. On a ici un modèle non linéaire, puisque \mathbb{E}[D_{x,t}]=(\cdot\exp[\log(E_{x,t})+A_x+B_x\cdot K_t]). On a à la fois une somme et un produit. Bref, on va partir vers un modèle nonlinéaire généralisé.

base=data.frame(DH=DECES$Male,DF=DECES$Female,
EH=EXPOS$Male,EF=EXPOS$Female,
age=rep(ages,length(annees)),
annee=rep(annees,each=length(ages)))
idx=which((base$DH<=0)|(base$DF<=0)|(base$EH<=0)|(base$EF<=0))
base=base[-idx,]
library(gnm)
reg=gnm( DH ~ offset(log(EH)) + as.factor(age) +
   Mult(as.factor(age),as.factor(annee)),
   family = poisson(link="log"),
   data=base)

On a nos estimateurs \hat{\mathbf{A}}, \hat{\mathbf{B}} et \hat{\mathbf{K}}.

Ax=reg$coefficients[2:111]
Bx=reg$coefficients[112:222]
Kt=reg$coefficients[223:length(reg$coefficients)]

On va passer sous silence la non identifiabilité du modèle. On peut représenter nos trois jeux de coefficients. Tout d’abord les \hat{\mathbf{A}} qui nous donnent l’évolution moyenne,

plot(ages[-1],Ax)

On peut aussi tracer \hat{\mathbf{K}} qui donne l’évolution temporelle.

plot(annees,-Kt)

Encore une fois, le modèle n’est pas identifiable. Si (\hat{\mathbf{A}},\hat{\mathbf{B}},\hat{\mathbf{K}}) est solution, il en est de même pour (\hat{\mathbf{A}},-\hat{\mathbf{B}},-\hat{\mathbf{K}}). Bref, le sens d’évolution n’a pas de sens. On peut représenter -\hat{\mathbf{K}} qui présente l’avantage d’être décroissant (et de décrire l’amélioration des conditions de vie). Enfin, traçons aussi -\hat{\mathbf{B}}

plot(ages,-Bx)

La difficulté est que pour prévoir l’espérance de vie, on a besoin de calculer q_{t,x} pour des valeurs importantes de t (non encore observée). Par exemple, on pourrait avoir envie de connaître q_{50,2020} (pour quelqu’un né en 1970). On a envie d’utiliser q_{50,2020}=\exp(\hat{A}_{50}+\hat{B}_{50}\hat{K}_{2020}). Le soucis c’est que \hat{K}_{2020} ne fait pas partie des grandeurs estimées \hat{\mathbf{K}}.

L’idée est de construire un prédicteur \tilde{K}_{2020} à partir de \{\hat{K}_{1816},\cdots,\hat{K}_{2015}\}. Classiquement (et c’est l’idée originale de Lee & Carter (1992)) on peut tenter un modèle ARIMA. Ici, la série est un peu trop régulière. On peut tenter soit un modèle exponentiel, soit un modèle linéaire (sur la série originale \hat{\mathbf{K}}, au delà de 1950)

idx=which(annees>=1950)
regk=lm(log(Kt[idx])~annees[idx])
futur=2016:2125
prevKt=exp(regk$coefficients[1]+regk$coefficients[2]*futur)
plot(annees,Kt,xlim=c(1815,2125),ylim=range(c(Kt,prevKt)))
points(futur,prevKt,col="red")
regk2=lm(Kt[idx]~annees[idx])
prevKt2=(regk2$coefficients[1]+regk2$coefficients[2]*futur)
points(futur,prevKt2,col="blue")

On peut alors construire la série des prédictions sur les données passées, \hat{q}_{x,t}=\exp[\hat{A}_{x}+\hat{B}_{x}\hat{K}_{t}], et sur les données futures \tilde{q}_{x,t}=\exp[\hat{A}_{x}+\hat{B}_{x}\tilde{K}_{t}].

Si on reste sur les données passées, c’est assez simple, avec ici les probabilités de décès pour 1880

BASE=data.frame(
  x=base$age,
  t=base$annee,
  indice_x=1+base$age,
  indice_t=base$annee-1815)
BASE$Ax=Ax[BASE$indice_x]
BASE$Bx=Bx[BASE$indice_x]
BASE$Kt=Kt[BASE$indice_t]
cste=reg$coefficients[1]
BASE$Qpred=exp(cste+BASE$Ax+BASE$Bx*BASE$Kt)
plot(BASE$x[BASE$t==1880],BASE$Qpred[BASE$t==1880],
log="y")

On peut faire la même chose pour le futur (ici pour l’année 2050) avec les deux modèles

BASE2=data.frame(
  x=rep(ages,length(futur)),
  t=rep(futur,each=length(ages)),
  indice_x=1+rep(ages,length(futur)),
  indice_t=rep(futur,each=length(ages))-2015)
BASE2$Ax=Ax[BASE2$indice_x]
BASE2$Bx=Bx[BASE2$indice_x]
BASE2$Kt1=prevKt[BASE2$indice_t]
BASE2$Kt2=prevKt2[BASE2$indice_t]
cste=reg$coefficients[1]
BASE2$Qpred1=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt1)
BASE2$Qpred2=exp(cste+BASE2$Ax+BASE2$Bx*BASE2$Kt2)
head(BASE2)
plot(BASE2$x[BASE2$t==2050],BASE2$Qpred1[BASE2$t==
2050],log="y")

pour la prévision exponentielle, et

pour la prévision linéaire On va ensuite bricoler pour créer le vecteur \mathbf{q} qui contiendra soit \hat{q}_{h,t+h}, soit \tilde{q}_{h,t+h}. En effet, pour faire un calcul d’espérance de vie, on ne regarde pas à un instant donnée, mais on suit une “cohorte” dans le temps. Par exemple, pour une personne née en 1968, on a la probabilité de décès (dans l’année) suivante

annee = 1968
vt=annee+ages
vx=ages
vq=rep(NA,length(vx))
sbase=data.frame(x=vx,t=vt)
for(i in 1:nrow(sbase)){
  if(sbase$t[i] <= 2015)
{vq[i]=BASE[(BASE$x==sbase$x[i])&(BASE$t==sbase$t[i]),"Qpred"]} 
if(sbase$t[i] > 2015) 
{vq[i]=BASE2[(BASE2$x==sbase$x[i])&(BASE2$t==sbase$t[i]),"Qpred2"]}
}
plot(ages,vq,log="y")
abline(v=2015-1968+.5,lty=2)

avec à gauche les valeurs estimées par notre modèle (on aurait pu utiliser les valeurs brutes), et à droite, les valeurs prédites (puisqu’on parle de probabilités de décès en 2030, entre autres).

Pour calculer l’espérance de vie à la naissance, on utilise le code suivant

sum(cumprod(exp(-vq[1:110])))
[1] 77.62047

On peut alors faire une petite fonction, afin de visualiser l’évolution de cette espérance de vie

Ex=function(annee = 1968,past="obs",fut="2"){
vt=annee+ages
vx=ages
sbase=data.frame(x=vx,t=vt,vq=NA)
for(i in 1:nrow(sbase)){
  if(sbase$t[i] <= 2015) {sbase$vq[i]=BASE[(BASE$x==sbase$x[i])&(BASE$t==sbase$t[i]),paste0("Q",past)]} if(sbase$t[i] > 2015) 
  {sbase$vq[i]=BASE2[(BASE2$x==sbase$x[i])&(BASE2$t==sbase$t[i]),paste0("Qpred",fut)]}
}
vP = cumprod(exp(-(sbase$vq[1:110])))
sum(vP)}

Aussi

ANNEES=1930:2010
E2=Vectorize(Ex)(ANNEES)
plot(ANNEES,E2)

Si on regarde en variation, on voit qu’il y a (en gros) une variation de 0.25 d’année en année (d’où le fameux “on gagne un trimestre d’espérance de vie tous les ans”)

plot(ANNEES[-1],diff(E2),ylim=c(0,1))

Comme on plafonne ici les âges (à 110 ans) on comprend que ce gain ne peut que dimunuer, ce qu’on observe ici.

Par contre, si on prend une prévision qui conserve la variation exponentielle de K_t, on obtient

E1=Vectorize(function(a) Ex(a,"obs","1"))(ANNEES)
plot(ANNEES,E1)

ce qui n’est pas réaliste. Encore moins en considérant la variation de cette courbe (beaucoup trop optimiste)

plot(ANNEES[-1],diff(E1),ylim=c(0,1))

Modéliser la consommation électrique

A partir de nos prévisions pour la température, on peut tenter de prévoir la consommation électrique. Rappelons que la série de consommation électrique ressemble à ca

plot(electricite[passe,"Load"],type="l")

On peut tenter un modèle assez simple, où la consommation à la date Y_t est fonction d’une tendance linéaire a+bt, de la position dans l’année (sous une forme non-paramétrique), de la température T_t (sous une forme non-linéaire, disons un polynôme de degré 3), de la températeure en t-1, T_{t-1} et de l’indice de production industriel IPI_t (sous une forme linéaire pour faire simple).

model1 = lm(Load~1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

L’idée de la fonction polynomial pour l’impact de la température venait du graphique suivant (où on représente la série de la consommation, une fois enlevée la tendance linéaire)

On peut aussi supposer une forme autorégressive, où Y_{t} sera fonction de Y_{t-1}

model2 = lm(Load~1+Load1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

On peut alors tenter de faire de la prévision. Le soucis avec ce second modèle, c’est la composante autorégressive. Pour prévoir Y_{t+h}, il faut utiliser la prévision faite en t+h-1, \hat Y_{t+h-1}.

base_prevision=data.frame(
Load_prev = NA,  Load1=c(electricite[length(passe),"Load"],rep(NA,110),
Time = electricite[futur,"Time"],
NumWeek = electricite[futur,"NumWeek"],
Temp = MODELES[,"y4"], 
Temp1 = c(electricite[length(passe),"Temp"],
MODELES[1:110,"y4"]),
IPI = electricite[futur,"IPI"])
for(t in 1:110){
 p = predict(model2,newdata=base_prevision)[t]
 base_prevision[t,"Load_prev"] = p
 base_prevision[t+1,"Load1"] = p}

On peut alors comparer la prévision \hat Y_{t} avec l’observation Y_{t}.

plot(1:111,electricite[futur,"Load"],type="l")
lines(1:111,base_prevision[,"Load_prev"],col="red")

On est bon les mois d’été (on capture bien le pic de la première quinzaine du mois d’août) par contre, on sous-estime fortement la consommation l’hiver. Et c’est gênant…

Pour finir, on peut oublier le passage par les variables explicatives, et tenter directement un modèle de série temporelles.

plot(electricite[passe,"Load"],type="l")

Le soucis, c’est le côté hétéroscédastique de la série, avec une pente sur le minimum plus faible que sur le maximum.

n=length(passe)
plot(1:n,electricite[passe,"Load"],type="l")
m=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),min)
M=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),max)
vm=which(electricite[passe,"Load"]%in%m$x)
vM=which(electricite[passe,"Load"]%in%M$x)
points(sort((1:n)[vm]),m$x,pch=19,col="blue")
points(sort((1:n)[vM]),M$x,pch=19,col="blue")
xm=((1:n)[vm])
xM=((1:n)[vM])
regm=lm(m$x~xm,col="blue")
regM=lm(M$x~xM,col="blue")
abline(regm,lty=2)
abline(regM,lty=2)

Une solution classique (et simple) est de passe au logarithme

plot(electricite[passe,"Load"],type="l",log="y")
y=electricite[passe,"Load"]
z = log(y)
plot(z,type="l")
B = data.frame(z=z,t=1:length(z))
abline(lm(z~t,data=B),lty=2,col="red")

Il faut ensuite enlever la tendance linéaire, pour stationnariser la série

z = residuals(lm(z~t,data=B))
Z = ts(z,start=1996,frequency=52)
modelz1 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
modelz1 = arima(Z,order = c(4,0,0), 
seasonal = list(order = c(1,0,1)))

Ce premier modèle est stationnaire, sans racine unité. On peut tenter d’introduire une racine unité saisonnière

modelz2 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(0,1,0)))
modelz2 = arima(Z,order = c(1,0,3), 
seasonal = list(order = c(1,1,0)))

Enfin, le dernier sera un peu plus simple

modelz3 = arima(Z,order = c(1,0,0), 
seasonal = list(order = c(2,0,0)))

On stocke alors toutes les prévisions dans une base

DONNEES_u = data.frame(
  u1 = predict(modelz1,n.ahead = 111)$pred,
  u2 = predict(modelz2,n.ahead = 111)$pred,
  u3 = predict(modelz3,n.ahead = 111)$pred)

On rajoute alors la tendance linéaire à la prévision du résidu

reg = lm(z~t,data=B)
dfutur = data.frame(t=futur)
DONNEES_z = data.frame(
  z1 = DONNEES_u$u1+predict(reg,newdata=dfutur),
  z2 = DONNEES_u$u2+predict(reg,newdata=dfutur),
  z3 = DONNEES_u$u3+predict(reg,newdata=dfutur)

On a besoin de la variance pour la prévision. En effet, classiquement, la prévision sera l’espérance de notre grandeur d’intérêt. Ici, on a fait un modèle linéaire sur \log Y, i.e. \log Y \sim\mathcal{N}(\mu,\sigma^2) et donc \mathbb{E}[Y]=\exp\left(\mu+\sigma^2/2\right)

sigma=summary(reg)$sigma
DONNEES_seu = data.frame(
  seu1=sqrt(predict(modelz1,n.ahead = 111)$se^2+sigma^2),
  seu2=sqrt(predict(modelz2,n.ahead = 111)$se^2+sigma^2),
  seu3=sqrt(predict(modelz3,n.ahead = 111)$se^2+sigma^2))

On fait ici l’hypothèse que les estimateurs des prévisions des deux modèles (la tendance linéaire et celle du modèle autorégressif) sont indépendants et donc on peut sommer les termes de variance. Aussi, les prévisions pour Y sont

DONNEES_y = data.frame(
  y1 = exp(DONNEES_z$z1+1/2*DONNEES_seu$seu1^2),
  y2 = exp(DONNEES_z$z2+1/2*DONNEES_seu$seu2^2),
  y3 = exp(DONNEES_z$z3+1/2*DONNEES_seu$seu3^2))

On peut représenter les prévisions des trois modèles (avec les observations)

plot(futur,electricite[futur,"Load"],type="l")
lines(futur,DONNEES_y$y1,col="red")
lines(futur,DONNEES_y$y2,col="blue")
lines(futur,DONNEES_y$y3,col="purple")

On peut comparer avec la prévision précédente,

lines(futur,base_prevision[,"Load_prev"],
col="orange")

On y perd un peu l’été, par contre, on s’est un peu amélioré l’hiver.

Forecasting Techniques for Demographics

Tomorrow, we will discuss in our cours forecasting tools for demographics. But first, we will see basic static tools. Before playing with longitudinal dataset, let us use “standard” life tables. Some (French) datasets are available on the INED website, but let us use the most popular ones, the TV8890 and TD8890. Those two tables can be downloaded from wikipedia,

url="https://fr.wikipedia.org/wiki/Table_de_mortalit%C3%A9"
download.file(url,"mortalite.html")
library(XML)
tables=readHTMLTable("mortalite.html")

The code to get the dataset is the following (here, we should be careful since there is a space in the number, as a separator for thousands)

TV8890_0=tables[[2]]
a1=as.numeric(as.character(TV8890_0[,1]))
a2=as.numeric(as.character(TV8890_0[,3]))
espace=function(x) gsub("[[:space:]]", "", x)
b1=espace(as.character(TV8890_0[,2]))
b2=espace(as.character(TV8890_0[,4]))
TV8890=data.frame(x=c(a1,a2),lx=as.numeric(c(b1,b2))

One can also read the second life table (note that there is a typo in wikipedia since here, the two tables are exactly the same)

TD8890_0=tables[[1]]
a1=as.numeric(as.character(TD8890_0[,1]))
a2=as.numeric(as.character(TD8890_0[,3]))
b1=espace(as.character(TD8890_0[,2]))
b2=espace(as.character(TD8890_0[,4]))
TD8890=data.frame(x=c(a1,a2),lx=as.numeric(c(b1,b2)))

It is possible to use that survival function to compute some sort to life expectancy at birth

sum(TV8890$lx)/100000-1
[1] 72.01518

One can visualize the survival probability (up to a 100,000 scaling constant)

plot(TV8890,type="l")

or the death probability, i.e. the probability to die at some specific age x, given that you did reach age x, also called the force of mortality

n=nrow(TV8890)
px=(TV8890$lx[1:(n-1)]-TV8890$lx[2:n])/
TV8890$lx[1:(n-1)]
x=TV8890$x[1:(n-1)]
plot(x,px,type="l",xlab="age")

A more popular visualization is obtained with a log scale for the probability

plot(x,px,type="l",log="y")

Finally, we can compute the density of the age at death

pbx=TV8890$lx[1:(n-1)]*px/100000
plot(x,pbx,type="l")

that can also be used to compute life expectancy

sum(x*pbx)
[1] 72.01518

That is for the static case. For longitudinal tables, we can use those from the Human Mortality Database. For instance, for France, we can get information of the number of death at age x during year t, as well as the exposure (number of people alive). For France, there are datasets available here

url="http://freakonometrics.free.fr/FranceDeaths_1x1.txt"
download.file(url,"FRD.txt")
url="http://freakonometrics.free.fr/FranceExposures_1x1.txt"
download.file(url,"FRE.txt")

and for Canada

url="http://freakonometrics.free.fr/CanadaDeaths_1x1.txt"
download.file(url,"CAD.txt")
url="http://freakonometrics.free.fr/CanadaExposures_1x1.txt"
download.file(url,"CAE.txt")

The following code can be use to read those files.

FRD=read.table("FRD.txt",skip = 3,header=TRUE)
tail(FRD)
      Year  Age Female  Male  Total
22195 2015  105 297.96 35.86 333.82
22196 2015  106 182.95 20.39 203.34
22197 2015  107 104.87 10.50 115.37
22198 2015  108  57.27  5.07  62.34
22199 2015  109  31.93  2.59  34.52
22200 2015 110+  33.03  1.61  34.64

Il est trop nul ce prof…

Que celui qui n’a jamais pensé ca une fois dans sa scolarité me jette la première pierre…

Cela dit, il y a ceux qui le pensent, et ceux qui l’affirment haut et fort, genre sur des sites comme rate my professors. J’ai toujours eu du mal avec ce genre de sites (et plus généralement extraire une information utile des évaluations des élèves) car on ne connaît pas les motivations des élèves. Quand on lit “l’examen final est impossible” puis le commentaire d’après “l’examen, c’est juste des exercices vus en classe, c’est vraiment trop facile“, on finit par être perplexe. Oui, assez souvent une classe est composée d’éléments très hétérogènes, et il est impossible de contenter tout le monde. Le plus étonnant quand on enseigne, c’est l’équilibre qui peut en résulter. Le même cours peut être une réussite, et le semestre suivant un fiasco, alors que le contenu et la forme peuvent être identiques. Je me souviens quand j’étais étudiant de ce prof, annoncé par la promo précédante comme le meilleur prof au monde, que je trouvais finalement so-so.

Et puis un jour, certains d’entre nous décident d’aller un cran plus loin, de sauter le pas, de passer d’élève à prof. Sans pour autant que ce soit une vocation, juste un concours de circonstances. Pour ma part, j’avais travaillé en salle de marché, puis comme actuaire à Hong Kong, avant de revenir à la Fédération Francaise des Assurances sur Paris. Et juste avant les vacances de Noël (il y a maintenant une quinzaine d’années) on me propose de donner un cours trois semaines plus tard, au retour des vacances (prof qui lâche son cours au dernier moment, impossible de trouver un remplacant). Je pensais faire de la recherche en entreprise, et j’ai découvert l’enseignement.

Et je n’ai pas arrêté depuis. Avec toutes sortes de formats, toutes sortes de publics… En particulier, j’ai fait beaucoup de formation continue. Et j’en profite pour la minute culturelle du blog : en avril 1792, l’Assemblée législative nationale discute de l’organisation générale de l’instruction publique en France, et Jean Antoine Nicolas Caritat, plus connu sous le nom de marquis de Condorcet, fait un discours très avant-gardiste

Nous avons observé que l’instruction ne devait pas abandonner les individus au moment où il sortent des écoles ; qu’elle devait embrasser tous les âges ; qu’il n’y en avait aucun où il ne fût utile et possible d’apprendre, et que cette seconde instruction est d’autant plus nécessaire, que celle de l’enfance a été resserrée dans des bornes plus étroites. (…) En continuant l’instruction pendant toute la durée de la vie, on empêchera les connaissances acquises dans les écoles de s’effacer promptement de la mémoire, on entretiendra dans les esprits l’activité utile (…), on pourra montrer enfin l’art de s’instruire par soi-même

Or justement, il y a plusieurs mois, le centre de formation Caritat (dédié au monde de l’actuariat de l’assurance et de la protection sociale) lancait un prix pédagogique. L’idée est intéressante : vous pensez qu’il manque de formation sur tel ou tel thème ? et vous pensez pouvoir animer une telle formation ? eh bien, lancez-vous ! on vous monte ce cours sur mesure ! Quand le prix a été lancé, on m’a demandé de faire partie du jury, et j’ai tout de suite accepté. Il s’agissait, pour des personnes sans aucune expérience dans le domaine de la formation (continue) de proposer un plan pour une journée de 7 heures, sur un sujet, avec un fil rouge, de préparer quelques slides. On en a sélectionné certains, et on les a convié pour un oral, en fin d’année… et la remise des prix aura lieu cette semaine. Les vainqueurs pourront se jeter à l’eau dans quelques mois, et monter leur formation, leur cours, pour de vrai. Je trouve l’idée merveilleuse, et j’espère qu’on aura d’autres candidatures intéressantes pour la prochaine édition ! A suivre donc..

Devoir, Méthodes de Prévisions (consommation électrique)

Pour le projet de l’UE de prévision, je demande une prévision de la consommation d’électricité. Pour construire les modèles, je donne accès à un gros fichier qui donne, pour toute la France, la consommation d’électricité, entre 2012 et 2014, heure par heure.

  • comme en cours, la première étape sera de proposer une série de 52 prévisions pour chaque semaine de 2015,
  • la seconde étape sera de proposer deux séries de 31 prévisions, pour chaque jour d’août 2015 et de décembre 2015
  • enfin, la dernière étape sera une prévision de 168 prévisions, pour chaque heure de la semaine du 19 au 25 octobre 2015

Pour aider (?), je mets un lien vers la température en France  dans quelques villes, paris, toulouse, grenoble, nantes, dunkerque (via ftp.ncdc.noaa.gov). Il sera possible d’utiliser la vraie température de 2015 pour faire les prévisions de consommation électrique.  Idéalement, je veux un fichier html (markdown) et 4 fichiers csv avec les 4 séries de prévision.

Moyennes de Prédicteurs

La dernière fois, nous avions vu 4 modèles de séries temporelles, pour modéliser la série hebdomadaire de températures. Le premier était obtenu avec auto.arima, les deux étaient ensuite ses modèles SARIMA (le second était une forte autocorélation à l’ordre 52, et le troisième avait une racine unité saisonnière), et le dernier était obtenu par la méthode de Buys-Ballot (cyclique, avec un bruit autocorrélé).

Au lieu de se demander lequel des quatre est le “meilleur”, on peut tenter de trouver une moyenne pondérée, \hat{y}_t = \sum_{i}\omega_i \hat{y}_t^{(i)}\omega=(\omega_i) est un point du simplexe, au sens où \omega_i\geq 0 et \omega_1+...+\omega_k=1. On veut ensuite trouver les “meilleurs” pondérations, par exemple au sens de la norme L^2. On va calibrer nos quatre modèles sur les m premières valeurs, et on compare ensuite la combinaison des prévions sur les 111 prochaines valeurs (et les vraies valeurs), \omega^\star=\text{argmin}\{\sum_{h=1}^{111}\left(y_{t+h}- \sum_{i}\omega_i \hat{y}_{t+h}^{(i)}\right)^2\}

Par exemple, on peut utiliser les 200 premières valeurs.

m=200
Ym=ts(y[1:m],start=1996,frequency=52)

On va ensuite estimer les 4 modèles sur ces 200 valeurs

model1=arima(Ym,order = c(4,0,0), 
seasonal = list(order = c(1,0,0))) 
model2=arima(Ym,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
model3=arima(Ym,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
regL=lm(Temp~0+as.factor(NumWeek),
data=electricite[1:m,])
E=residuals(regL)
model4=arima(ts(E,start=1996,
frequency=52),order = c(1,0,0))

On peut ensuite collecter les prédictions

DONNEES = data.frame(
  y1=predict(model1,n.ahead = 111)$pred,
  y2=predict(model2,n.ahead = 111)$pred,
  y3=predict(model3,n.ahead = 111)$pred,
  y4=predict(model4,n.ahead = 111)$pred+
     predict(regL,newdata=electricite[m+(1:111),]),
  vraie=electricite[m+(1:111),"Temp"])

Pour créer une combinaison linéaire des prévisions, on utilise

a=rep(1/4,4)
y_prev = as.matrix(DONNEES[,1:4]) %*% a

On peut ainsi visualiser les 4 prévisions, leur combinaison linéaire (ici avec des équi-pondérations), et leurs observations

plot(1:111,DONNEES[,5],type="l",lwd=2)
for(i in 1:4) lines(1:111,DONNEES[,i],col="blue")
lines(1:111,yp,col="red")

Pour trouver la valeur “optimale” des pondérations, qui minimisera la somme des carrés des erreurs, on utilise le code suivant (on fera une transformation des données, ici sur le simplexe, pour s’assurer que l’optimisation donnera un élément du simplexe)

erreur=function(a) sum((DONNEES[,1:4] 
%*% a-DONNEES[,5])^2)
erreur2=function(b){
  a=rep(NA,length(b)+1)
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  erreur(a)
}

On obtient l’optimum, et on peut obtenir une prévision

b=optim(par=c(0,0,0),erreur2)$par
a=rep(NA,length(b)+1)
a[1:length(b)]=exp(b)/(1+sum(exp(b)))
a[length(b)+1]=1/(1+sum(exp(b)))
yp=DONNEES[,1:4] %*% a
lines(1:111,yp,col="red",lwd=4)

On peut alors faire changer le nombre d’observations utilisées

Techniquement, ca pose un petit soucis, car il faut s’assurer la convergence de deux algorithmes : celui d’estimation des paramètres des SARIMA, celui de recherche du paramètre de poids.

poids_optimaux=function(m){
  Ym=ts(y[1:m],start=1996,frequency=52)
  TRY=try(model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0))),FALSE)
  if(inherits(TRY, "try-error")) model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0)),method="CSS")
  TRY=try(model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0))),FALSE)
  if(inherits(TRY, "try-error")) model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0)),method="CSS")
  TRY=try(model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1))),FALSE)
  if(inherits(TRY, "try-error")) model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1)),method="CSS")
regL=lm(Temp~0+as.factor(NumWeek),data=electricite[1:m,])
  E=residuals(regL)
model4=arima(ts(E,start=1996,frequency=52),order = c(1,0,0))
  DONNEES = data.frame(
    y1=predict(model1,n.ahead = 111)$pred,
    y2=predict(model2,n.ahead = 111)$pred,
    y3=predict(model3,n.ahead = 111)$pred,
    y4=predict(model4,n.ahead = 111)$pred+
      predict(regL,newdata=electricite[m+(1:111),]),
    vraie=electricite[m+(1:111),"Temp"])
  erreur=function(a) sum((as.matrix(DONNEES[,1:4]) %*% a-DONNEES[,5])^2)
  erreur2=function(b){
    a=rep(NA,length(b)+1)
    a[1:length(b)]=exp(b)/(1+sum(exp(b)))
    a[length(b)+1]=1/(1+sum(exp(b)))
    erreur(a)
  }
  TRY=try(b <- optim(par=c(0,0,0),erreur2)$par, FALSE)
  a=rep(NA,length(b)+1)
  if(!inherits(TRY, "try-error")) {
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  }
  return(a)
}

On peut alors regarder l’évolution des pondérations au cours du temps. On peut vectoriser, mais si l’algorithme aboutit sur une erreur, c’est plus simple d’identifier

VM=seq(200,480,by=10)
VA=matrix(NA,length(VM),4)
for(i in 1:length(VM)) VA[i,]=poids_optimaux(VM[i])
colnames(VA)=c("auto","SARMA","SARIMA","Buys Ballot")

(pour débugger, c’est plus simple les boucles). On peut alors obtenir le graphique suivant, avec en rose Buys-Ballot, et en vert le modèle SARIMA, avec une racine unité saisonnière,

tva=t(VA)
barplot(tva,legend = rownames(counts),
 col=c("yellow","orange","green","pink"))

On se rend compte qu’assez vite, le modèle le plus intéressant est le modèle de Buys Ballot.

Il est possible de changer de fonction de perte, par exemple si on s’intéresse à un quantile à 90%,

tau=.9
loss=function(e) e*(tau-(e<=0)*1)

Dans la fonction, on utilise alors

erreur=function(a) sum(loss(as.matrix(
DONNEES[,1:4]) %*% a-DONNEES[,5]))

C’est, cette fois, beaucoup plus incertain, même si les deux modèles les plus intéressants sont toujours le SARIMA et Buys-Ballot.

An Open Lab-Notebook Experiment