Le livre Méthodes de prévision en finance est sorti, en France, cette semaine. L’application que je propose pour montrer la modélisation par copules (par opposition aux modèles GARCH multivariés en particulier) est basée sur un vieux billet sur le blog (ou les codes permettent de reproduire l’analyse que j’en fait), remis en forme dans un court article, Prévision avec des copules en finance
La troisième (et dernière) partie de mon intervention sur la science des données, dans le cadre du cours INF7100 portera sur une introduction aux modèles prédictifs. Le plan sera le suivant
La seconde partie de transparents pour le cours ACT6420 est en ligne. Il s’agit d’un document qui permettra à tout le monde garder des traces, mais je n’entends pas forcément les suivre… Par exemple, ce soir, je propose de jouer un peu avec les données d’incendies à Chicago
> chicago=read.table("http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)
ou une base sur les états américains (cf un ancien billet)
> US=read.table("http://freakonometrics.free.fr/US.txt", + header=TRUE,sep=";")
Sinon, la semaine passée, on avait évoqué la matrice chapeau (dont on reparlera davantage lorsque l’on évoquera la phase de diagnostique). L’idée était de noter que
avec
On peut visualiser cette matrice avec le code suivant
> X=cbind(1,cars$speed)
> M=X %*% solve( t(X) %*% X) %*% t(X)
> Q=matrix(as.numeric(cut(M,c(-10,quantile(M,(1:5)/6),10),labels=1:6)),50,50)
> plot(1:50,1:50,col="white",axes=FALSE,xlab="",ylab="",ylim=c(50,0))
> library(RColorBrewer)
> couleur=rev(brewer.pal(6,"RdBu"))
> for(i in 1:50){for(j in 1:50){ points(i,j,col=couleur[Q[i,j]],pch=15,cex=1)}}
> axis(1)
> axis(2)
avec des cases rouges pour les valeurs importantes (et positive) et des case bleues pour des valeurs importantes (mais négatives)
Cet session d’été, je redonne le cours ACT6420, modèles de prévision. Le plan de cours est maintenant en ligne (comme souvent, les salles de cours ne sont pas encore attribuées définitivement mais elles seront disponibles sur le site de l’UQAM). Les transparents pour la première partie de cours sont en ligne (il s’agit de rappels de statistique, et d’une introduction aux modèles de régression).
Des exercices et des codes informatiques seront mis en ligne sur /courses/modeles-de-prevision.
Durant le cours, j’essayerais d’illustrer avec le maximum d’exemples sur des données réelles, à partir de codes R. Quelques séances de démonstrations se feront en salle informatique, afin de découvrir les rudiments de R. Pour aller plus loin, je recommande
ainsi que quelques tutoriels en ligne, par exemple, dans un ancien billet, mais aussi sur https://youtube.com/.
Hier, sur Twitter, @JF_Godbout partageait un joli graphique relatif aux élections québécoises, avec les nombres de votes obtenus (ici en pourcentage des votes totaux) et le pourcentage de sièges que cela donne,
Il faut dire qu’hier, c’était jour d’élection. Ce sont des élections à un tour, avec plusieurs partis (disons 4 grands partis si on se limite aux élections les plus récentes). Jusqu’à lundi, le PQ (Parti Québécois) avait une majorité (relative), et depuis lundi, le PLQ (Parti Libéral) est revenu au pouvoir. Et il y a à côté de ces deux grands partis d’autres partis, que l’on va regrouper. @JF_Godbout m’ayant envoyé ses données, j’ai pu reproduire les graphiques.
La grande différence ici, entre nos approches est – je pense – la méthode de lissage: alors que @JF_Godbout faisait des régressions locales, j’ai opté pour des splines.
reg=glm(Seat~bs(Election),data=d,family=binomial) yp=predict(reg,newdata=data.frame(Election=1970:2015),type="response")
J’utilise des splines parce que ça se généralise facilement à d’autres lois, et qu’on peut aussi les utiliser quand on régresse sur plusieurs variables continues. Si on rajoute les résultats de l’élection de lundi, on obtient le graphique suivant
Ce qui me gêne dans cette modélisation, c’est qu’on considère nos trois dernières régressions comme indépendantes…. alors que ce n’est pas le cas: on a trois modalités possibles, et la somme doit forcément faire 1. D’ailleurs, quand on somme les prédictions, on dépasse 1.
Avec deux partis, notre régression logistique nous garantie d’avoir une somme à un, mais plus ici. A moins d’utiliser un modèle multinomial. Si on oublie le lissage, l’idée est d’écrire
qui est tout simplement l’extension du modèle logistique. La difficulté technique est que, sous R, quand on fait une régression logistique, la variable d’intérêt doit être une variable prenant trois facteurs, et non pas un vecteur de probabilités. Au lieu d’utiliser des taux, ou des pourcentages, je devrais utiliser des nombres d’observations de chaque modalités. On va donc crééer une fausse base
> vote=data.frame(
+ plq=plq$Vote,
+ pq=pq$Vote,
+ autre=autre$Vote,
+ Election=pq$Election)
> vote
plq pq autre Election
1 0.5465 0.3022 0.1513 1972
2 0.3378 0.4137 0.2485 1976
3 0.4608 0.4926 0.0466 1981
4 0.5599 0.3869 0.0532 1985
5 0.4995 0.4016 0.0989 1989
6 0.4440 0.4475 0.1085 1994
7 0.4355 0.4287 0.1358 1998
8 0.4599 0.3324 0.2077 2003
9 0.3308 0.2835 0.3857 2007
10 0.4208 0.3517 0.2275 2008
11 0.3120 0.3195 0.3685 2012
12 0.4140 0.2538 0.3312 2014
> M=t(round(vote[,1:3]*1000))
> M
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
plq 546 338 461 560 500 444 436 460 331 421 312 414
pq 302 414 493 387 402 448 429 332 284 352 320 254
autre 151 249 47 53 99 108 136 208 386 227 368 331
On réplique alors autant de fois le nom du parti qu’on a de (pseudo) votes
> Y=X=NULL
> for(i in 1:3){
+ for(j in 1:ncol(M)){
+ Y=c(Y,rep(c("PLQ","PQ","Autre")[i],M[i,j]))
+ X=c(X,rep(pq$Election[j],M[i,j]))
+ }}
> basev=data.frame(Y=as.factor(Y),X=X)
C’est un peu artificiel, j’en conviens, mais ça marche… On peut alors faire la régression
> library(nnet) > regv=multinom(Y~bs(X),data=basev) # weights: 15 (8 variable) initial value 13186.643301 iter 10 value 12241.931805 final value 12221.591172 converged
et des prévisions
> ypv=predict(regv,newdata=
+ data.frame(X=1973:2014),type="probs")
> ypv
Autre PLQ PQ
1 0.18426281 0.4705415 0.3451957
2 0.16051060 0.4767408 0.3627486
3 0.14149153 0.4805514 0.3779571
4 0.12631826 0.4826100 0.3910718
5 0.11426236 0.4834165 0.4023212
6 0.10474046 0.4833543 0.4119053
7 0.09729304 0.4827117 0.4199953
8 0.09156246 0.4817016 0.4267360
9 0.08727323 0.4804776 0.4322491
10 0.08421535 0.4791477 0.4366370
11 0.08223095 0.4777835 0.4399855
12 0.08120364 0.4764292 0.4423672
13 0.08105030 0.4751066 0.4438431
14 0.08171457 0.4738203 0.4444651
Et on notera que cette fois, nos prévisions somment à un.
Cela dit, visuellement, on obtient des choses proches de ce qu’on avait auparavant,
mais intellectuellement, je trouve plus logique d’utiliser autre chose que des régressions logistiques quand on est dans une élection où plusieurs partis sont en compétition.
Pour la fin de cours de modèles de prévision , quelques transparents sur l’identification et l’estimation des modèles SARIMA, quelques compléments sur les tests (racines unités et non-stationnarité, ainsi que saisonnalité), et enfin, quelques pistes pour construire des prédictions (avec une quantification de l’incertitude), avec des codes R. Les transparents sont en ligne ici (même si la page de garde est identique aux autres, il s’agit de nouveau matériel). Maintenant que les transparents sont finis (et en ligne) les prochains billets seront orientés autour de la modélisation et des aspects computationnels.
Tous les profs qui surveillent leur examen se posent la même question: que signifie un départ au bout d’une heure et demi, si l’examen est supposé durer trois heures? Que l’examen était trop facile? trop dur?
Loin de moi l’idée de blâmer les étudiants qui partent tôt ! Tout d’abord parce que j’ai toujours souvent été dans les premiers à quitter la salle lors d’examens, ou de concours. Et parce que je ne pense pas qu’il y ait un lien évident entre le fait de partir tôt, ou tard, et la note obtenue.
Mais on peut regarder (et répondre au commentaire de @kl suite à l’examen de la session passée). Car l’examen du cours ACT6420 avait lieu la semaine passée, et il était court. Autrement dit, 95% des étudiants étaient sortis avant que j’annonce la fin de l’épreuve. Et j’ai noté sur toutes les copies l’heure de sortie (pour ceux qui veulent plus d’information sur le protocole de l’expérience, les étudiants ne savaient pas que je notais l’heure).
Il s’agissait d’un examen à choix multiples, ce qui garantit une forme d’impartialité dans la correction (sinon on pourrait objecter que si je note en prenant le paquet de copie tel qu’il est arrivé, les notes pour les élèves partis tôt – corrigés en dernier – et ceux partis à la fin – corrigés en premier – pourraient ne pas être vraiment comparables). Et comme à la dernière session, je demandais aux étudiants de prédire leur nombre de réponses. La question était facultative, mais les élèves pouvaient gagner un point en répondant (à condition de prédire exactement le bon nombre).
Les données sont (partiellement) en ligne ici
> act6420=read.table( + "http://freakonometrics.blog.free.fr/public/ data/baseact6420a1.txt", + sep=";",header=FALSE,skip=1) > names(act6420)=c("n","heure","obs","pred")
On peut alors comparer les notes en fonction de l’heure de sortie, ainsi que la note prédites (les notes sont sur 30, ou plutôt les nombres de bonnes réponses), en faisant des régressions linéaires simples,
> summary(lm(obs~heure,data=act6420)) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.1519 6.0546 0.521 0.6039 heure 1.3209 0.5391 2.450 0.0162 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.819 on 93 degrees of freedom Multiple R-squared: 0.06063, Adjusted R-squared: 0.05053 F-statistic: 6.002 on 1 and 93 DF, p-value: 0.01616 > summary(lm(pred~heure,data=act6420)) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.1178 8.8373 0.126 0.8996 heure 1.7234 0.7889 2.185 0.0317 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.999 on 85 degrees of freedom (8 observations deleted due to missingness) Multiple R-squared: 0.05316, Adjusted R-squared: 0.04202 F-statistic: 4.773 on 1 and 85 DF, p-value: 0.03167
On observe qu’à la fois les bons élèves (ceux avec une bonne note) et les élèves qui se pensent bons (ceux qui ont prédit une note élevée) sortent plutôt tard. Avec dans les deux cas, une significativité de la variable heure de sortie sur la prédiction de la note. Mais peut-être est-ce plus complexe que ça, et qu’il se cache un effet non-linéaire. On essayer de lisser un peu, avec des splines
> D=data.frame(Y=note,X=heure) > library(splines) > rg=lm(Y~bs(X,df=5),data=D) > newheure=seq(10.2,12.1,by=.01) > notep=predict(rg,newdata=data.frame(X= + newheure),interval="confidence")
On obtient les régressions (non-linéaires) suivantes
ou, si on souhaite rajouter des intervalles de confiance autour de la prédiction faite par nos modèles,
(avec en abscisse l’heure de sortie, l’examen finissant un peu après midi). L’analyse est un peu plus complexe qu’il n’y paraissait, avec quelques bons élèves qui finissent très tôt. Si on souhaite d’ailleurs mieux comprendre le lien entre la note obtenue (ou plutôt, encore une fois le nombre de bonnes réponses données, car aucun bonus n’est pris en compte ici), on peut faire un graphique proche de celui fait à la session passée
Le trait bleu est ici la régression obtenue à la session passée. On retrouve exactement la même chose sur les deux sessions, étonnant, non?
Mais j’arrête là pour l’analyse de ces notes. Pour ceux qui veulent en savoir plus, je mettrais bientôt en ligne une correction, avec des statistiques par réponse. On en profitera pour faire une discussion sur l’utilisation de règles démocratiques, à savoir, la majorité a-t-elle toujours raison? Et pour ceux qui n’auraient pas encore compris, on a fait ici une analyse de corrélations entre les heures de sortie, et les notes, et non pas une analyse de causalité. Donc inutile de se forcer à attendre la fin de l’examen pour sortir la prochaine fois, ça ne fera a pas (a priori) monter la note.
Après discussion, la date limite pour me renvoyer le second devoir est fixée à vendredi 5 octobre midi, par courriel. Je veux un fichier pdf par équipe, avec une page de garde indiquant le nom des membres du groupe, et le mot clé retenu.
Maintenant, je vais en profiter pour mettre en ligne le code tapé de matin, en cours. On avait travaillé sur la modélisation de la série des ventesde voiture, au Québec. La série est en ligne,
X=read.table( "http://freakonometrics.blog.free.fr/public/ data/car-sales-quebec.csv", header=TRUE,sep=";",nrows=108) Xt=ts(X[,2],start=c(1960,1),frequency=12)
On peut regarder l’évolution de cette série temporelles (car les dessins, et la visualisation, c’est important),
plot(Xt)
On note une tendance linéaire (ou qu’on pourrait supposer linéaire), que l’on va estimer,
X=as.numeric(Xt) temps=1:length(X) base=data.frame(temps,X) reg=lm(X~temps)
et que l’on pourra représenter
plot(temps,X,type="l",xlim=c(0,120)) abline(reg,col="red")
Si on veut une série stationnaire (et c’est effectivement ce que l’on cherche), on va retrancher cette tendance à notre série, et travailler sur la série résiduelle.
La série résiduelle est ici
Y=X-predict(reg) plot(temps,Y,type="l")
que l’on devrait pouvoir supposer stationnaire. Classiquement, on regarde les autocorrélations,
acf(Y,36,lwd=4)
où on repère un beau cycle annuel, mais on va surtout utiliser les autocorrélations partielles pour identifier l’ordre de la composante autorégressive.
pacf(Y,36,lwd=4)
La 12ème est non-nulle, et tous les autres ensuite sont significativement nulles (ou presque). On va donc tenter un .
> fit.ar12=arima(Y,order=c(12,0,0)) > fit.ar12 Series: Y ARIMA(12,0,0) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.1975 0.0832 -0.1062 -0.1212 0.1437 -0.1051 0.0319 s.e. 0.0838 0.0809 0.0826 0.0843 0.0850 0.0833 0.0854 ar9 ar10 ar11 ar12 intercept -0.0332 -0.0616 0.2635 0.4913 -148.3180 s.e. 0.0853 0.0840 0.0840 0.0841 384.5095 sigma^2 estimated as 2177974: log likelihood=-946.75 AIC=1921.51 AICc=1926.03 BIC=1959.06
La douzième composante semble significative, si on repense au test de Student,
u=seq(-6,6,by=.1) plot(u,dnorm(u),type="l") points(0.4913/0.0841,0,pch=19,col="red")
En revanche, la huitième ne l’est pas
points(-0.1018/0.0847,0,pch=19,col="blue")
(voilà pour le petit retour sur la lecture des tests). Vérifions maintenant que le bruit résiduel est bien blanc… Si on regarde les autocorrélations
acf(residuals(fit.ar12),36,lwd=4)
Mais ce n’est pas un test de bruit blanc. En particulier, si une autocorrélation semblerait significative, on pourrait accepter que – globalement – les autocorrélations soient très proches de 0. On va alors faire des tests de Box-Pierce
> Box.test(residuals(fit.ar12),lag=12, + type='Box-Pierce') Box-Pierce test data: residuals(fit.ar12) X-squared = 7.7883, df = 12, p-value = 0.8014
Pour l’interpétation du test, c’est toujours pareil: on a une statistique de l’ordre de 7, et la loi sous-jacente (la somme des carrés des 12 premières autocorrélations du bruit) est une loi du chi-deux, à 12 degrés de liberté,
u=seq(0,30,by=.1) plot(u,dchisq(u,df=12),type="l") points(7.7883,0,pch=19, col="red")
On accepte l’hypothèse que les 12 premières autocorrélations soient nulles. Si on va un peu plus loin,
> Box.test(residuals(fit.ar12),lag=18, + type='Box-Pierce') Box-Pierce test data: residuals(fit.ar12) X-squared = 20.3861, df = 18, p-value = 0.3115
on a la même conclusion,
Pour rappels, on peut retrouver à la main la p-value,
> 1-pchisq(20.3861,df=18) [1] 0.3115071
Faisons notre petite dessin de toutes les p-value, pour s’assurer que le bruit est bien un bruit blanc,
BP=function(h) Box.test(residuals(fit.ar12),lag=h, type='Box-Pierce')$p.value plot(1:24,Vectorize(BP)(1:24),type='b', ylim=c(0,1),col="red") abline(h=.05,lty=2,col="blue")
Cette fois c’est bon, on tient notre premier modèle, un . Si on veut aller plus loin, on peut regarder ce que donnerait de la sélection automatique,
> library(caschrono) > armaselect(Y,nbmod=5) p q sbc [1,] 14 1 1635.214 [2,] 12 1 1635.645 [3,] 15 1 1638.178 [4,] 12 3 1638.297 [5,] 12 4 1639.232
On peut être tenté d’aller voir ce qui se passerait en rajoutant une composante moyenne mobile à notre modèle précédant,
> fit.arma12.1=arima(Y,order=c(12,0,1)) > fit.arma12.1 Series: Y ARIMA(12,0,1) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.0301 0.1558 -0.0941 -0.1461 0.1063 -0.0688 -0.002 s.e. 0.1235 0.0854 0.0757 0.0784 0.0807 0.0774 0.080 ar9 ar10 ar11 ar12 ma1 intercept -0.0646 -0.0798 0.2538 0.5786 0.2231 -131.3495 s.e. 0.0802 0.0766 0.0751 0.0861 0.1393 368.8156 sigma^2 estimated as 2127759: log likelihood=-945.65 AIC=1921.31 AICc=1926.52 BIC=1961.54
Le nouveau coefficient, est à peine significatif. Éventuellement avec un seuil à 10%… Pourquoi pas? Si on regarde les résidus, sans grande surprise, on a toujours un bruit blanc, encore plus blanc qu’auparavant (en violet sur le dessin ci-dessous)
BP=function(h) Box.test( residuals(fit.arma12.1),lag=h, type='Box-Pierce')$p.value plot(1:24,Vectorize(BP)(1:24),type='b', ylim=c(0,1),col="red") abline(h=.05,lty=2,col="blue") BP=function(h) Box.test(residuals(fit.ar12),lag=h, type='Ljung-Box')$p.value lines(1:24,Vectorize(BP)(1:24),col="purple", type="b")
On peut aussi aller voir parmi les modèles proposés, en particulier le modèle .
> fit.ar14= + arima(Y,order=c(14,0,0),method="CSS") > fit.ar14 Series: Y ARIMA(14,0,0) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.2495 0.2105 -0.0584 -0.1569 0.1282 -0.1152 0.0268 s.e. 0.0956 0.0972 0.0854 0.0830 0.0838 0.0840 0.0847 ar9 ar10 ar11 ar12 ar13 ar14 intercept -0.0327 -0.1116 0.2649 0.5887 -0.1575 -0.1572 80.5 s.e. 0.0855 0.0851 0.0853 0.0886 0.1031 0.0999 338.9 sigma^2 estimated as 2218612: part log likelihood=-942.31
C’est un peu tiré par les cheveux, mais on pourrait accepter l’hypothèse que soit significativement non-nuls. Mais on est encore limite… Allez, on l’accepte aussi dans notre gang de modèle.
On a finalement trois modèles. Si on fait un peu de backtesting, sur les 12 derniers mois,
T=length(Y) backtest=12 subY=Y[-((T-backtest+1):T)] subtemps=1:(T-backtest) plot(temps,Y,type="l") lines(subtemps,subY,lwd=4) fit.ar12.s=arima(subY, order=c(p=12,d=0,q=0),method="CSS") fit.arma12.1.s=arima(subY, order=c(p=12,d=0,q=1),method="CSS") fit.ar14.s=arima(subY, order=c(p=14,d=0,q=0),method="CSS") p.ar12=predict(fit.ar12.s,12) pred.ar12=as.numeric(p.ar12$pred) p.arma12.1=predict(fit.arma12.1.s,12) pred.arma12.1=as.numeric(p.arma12.1$pred) p.ar14=predict(fit.ar14.s,12) pred.ar14=as.numeric(p.ar14$pred)
on obtient les prévisions suivantes (avec les valeurs observées dans la première colonne)
> (M=cbind(observé=Y[((T-backtest+1):T)], + modèle1=pred.ar12, + modèle2=pred.arma12.1, + modèle3=pred.ar14)) observé modèle1 modèle2 modèle3 97 -4836.2174 -5689.3331 -5885.4486 -6364.2471 98 -3876.4199 -4274.0391 -4287.2193 -4773.8116 99 1930.3776 1817.8411 2127.9915 2290.1460 100 3435.1751 4089.3598 3736.1110 4039.4150 101 7727.9726 6998.9829 7391.6694 7281.4797 102 2631.7701 3456.8819 3397.5478 4230.5324 103 -509.4324 -2128.6315 -2268.9672 -2258.7216 104 -1892.6349 -3877.7609 -3694.9409 -3620.4798 105 -4310.8374 -3384.0905 -3430.4090 -2881.4942 106 2564.9600 -504.6883 -242.5018 183.2891 107 -1678.2425 -1540.9904 -1607.5996 -855.7677 108 -4362.4450 -3927.4772 -3928.0626 -3718.3922 > sum((M[,1]-M[,2])^2) [1] 19590931 > sum((M[,1]-M[,3])^2) [1] 17293716 > sum((M[,1]-M[,4])^2) [1] 21242230
I.e. on aurait envie de retenir le second modèle, . On va maintenant l’utiliser pour faire un peu de prévision,
library(forecast) fit.arma12.1= arima(Y,order=c(12,0,1)) fit.arma12.1 PREDARMA=forecast(fit.arma12.1,12) plot(Y,xlim=c(1,120),type="l") temps=T+1:12 polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2], rev(PREDARMA$upper[,2])),col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1], rev(PREDARMA$upper[,1])),col="orange",border=NA) lines(temps,PREDARMA$mean,col="red")
Cela dit, on peut aussi aller beaucoup plus loin dans la prévision,
PREDARMA=forecast(fit.arma12.1,120) plot(Y,xlim=c(1,210),type="l") temps=T+1:120 polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2], rev(PREDARMA$upper[,2])),col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1], rev(PREDARMA$upper[,1])),col="orange",border=NA) lines(temps,PREDARMA$mean,col="red")
Bon, on y est presque, car on a modélisé , la série obtenue en enlevant la tendance linéaire. Pour remonter sur
, on va rajouter la tendance à la prévision faite auparavant,
X=as.numeric(Xt) temps=1:length(X) plot(temps,X,type="l",xlim=c(0,210), ylim=c(5000,30000)) base=data.frame(temps,X) reg=lm(X~temps) abline(reg,col="red") PREDTENDANCE=predict(reg,newdata= data.frame(temps=T+1:120)) temps=T+1:120 polygon(c(temps,rev(temps)),c(PREDTENDANCE+ PREDARMA$lower[,2],rev(PREDTENDANCE+PREDARMA$upper[,2])), col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDTENDANCE+ PREDARMA$lower[,1],rev(PREDTENDANCE+PREDARMA$upper[,1])), col="orange",border=NA) lines(temps,PREDTENDANCE+PREDARMA$mean,col="red")
