Freakonometrics

To content | To menu | To search

Tag - tendance

Entries feed - Comments feed

Friday, December 7 2012

Modélisation et prévision, cas d'école

Quelques lignes de code que l'on reprendra au prochain cours, avec une transformation en log, et une tendance linéaire. Considérons la recherche du mot clé headphones, au Canada, la base est en ligne sur http://freakonometrics.blog.free.fr/...

> report=read.table(
+ "report-headphones.csv",
+ skip=4,header=TRUE,sep=",",nrows=464)
> source("http://freakonometrics.blog.free.fr/public/code/H2M.R")
> headphones=H2M(report,lang="FR",type="ts")
> plot(headphones)

Mais le modèle linéaire ne devrait pas convenir, car la série explose,

> n=length(headphones)
> X1=seq(12,n,by=12)
> Y1=headphones[X1]
> points(time(headphones)[X1],Y1,pch=19,col="red")
> X2=seq(6,n,by=12)
> Y2=headphones[X2]
> points(time(headphones)[X2],Y2,pch=19,col="blue")

Il est alors naturel de prendre le logarithme de la série,

> plot(headphones,log="y")

C'est cette série que l'on va modéliser (mais c'est bien entendu la première série, au final, qu'il faudra prévoir). On commence par ôter la tendance (ici linéaire)

> X=as.numeric(headphones)
> Y=log(X)
> n=length(Y)
> T=1:n
> B=data.frame(Y,T)
> reg=lm(Y~T,data=B)
> plot(T,Y,type="l")
> lines(T,predict(reg),col="purple",lwd=2)

On travaille alors sur la série résiduelle. 

> Z=Y-predict(reg)
> acf(Z,lag=36,lwd=6)
> pacf(Z,lag=36,lwd=6)

On peut tenter de différencier de manière saisonnière,

> DZ=diff(Z,12)
> acf(DZ,lag=36,lwd=6)
> pacf(DZ,lag=36,lwd=6)

On ajuste alors un processus ARIMA, sur la série différenciée,

> mod=arima(DZ,order=c(1,0,0),
+ seasonal=list(order=c(1,0,0),period=12))
> mod
 
Coefficients:
ar1     sar1  intercept
0.7937  -0.3696     0.0032
s.e.  0.0626   0.1072     0.0245
 
sigma^2 estimated as 0.0046:  log likelihood = 119.47

Mais comme c'est la série de base qui nous intéresse, on utilise une écriture SARIMA,

> mod=arima(Z,order=c(1,0,0),
+ seasonal=list(order=c(1,1,0),period=12))

On fait alors la prévision de cette série.

> modpred=predict(mod,24)
> Zm=modpred$pred
> Zse=modpred$se

On utilise aussi le prolongement de la tendance linéaire,

> tendance=predict(reg,newdata=data.frame(T=n+(1:24)))

Pour revenir enfin à notre série initiale, on utilise les propriétés de la loi lognormales, et plus particulièrement la forme de la moyenne, pour prédire la valeur de la série,

> Ym=exp(Zm+tendance+Zse^2/2)

Graphiquement, on a 

> plot(1:n,X,xlim=c(1,n+24),type="l",ylim=c(10,90))
> lines(n+(1:24),Ym,lwd=2,col="blue")

Pour les intervalles de confiance, on peut utiliser les quantiles de la loi lognormale,

> Ysup975=qlnorm(.975,meanlog=Zm+tendance,sdlog=Zse)
> Yinf025=qlnorm(.025,meanlog=Zm+tendance,sdlog=Zse)
> Ysup9=qlnorm(.9,meanlog=Zm+tendance,sdlog=Zse)
> Yinf1=qlnorm(.1,meanlog=Zm+tendance,sdlog=Zse)
> polygon(c(n+(1:24),rev(n+(1:24))),
+ c(Ysup975,rev(Yinf025)),col="orange",border=NA)
> polygon(c(n+(1:24),rev(n+(1:24))),
+ c(Ysup9,rev(Yinf1)),col="yellow",border=NA)

Tuesday, September 25 2012

Prévision des vente de voitures au Québec

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 ventes de 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 http://latex.codecogs.com/gif.latex?AR(12).

> 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 http://latex.codecogs.com/gif.latex?\phi_{12} 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 http://latex.codecogs.com/gif.latex?\phi_{8} 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 http://latex.codecogs.com/gif.latex?AR(12). 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, http://latex.codecogs.com/gif.latex?\theta_1 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='Box-Pierce')$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 http://latex.codecogs.com/gif.latex?AR(14).

> 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 http://latex.codecogs.com/gif.latex?\phi_{14} 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, http://latex.codecogs.com/gif.latex?ARMA(12,1). 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é http://latex.codecogs.com/gif.latex?(Y_t), la série obtenue en enlevant la tendance linéaire. Pour remonter sur http://latex.codecogs.com/gif.latex?(X_t), 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")

Friday, September 21 2012

modélisation des données du Nil, et prévision

Un billet pour remettre le code de ce soir, sur la modélisation du Nil (sans le détail de tous les graphiques). Tout d'abord, dans un premier temps, on notait qu'il devait y avoir une rupture dans la série, et on optait pour un modèle basé sur les dernières observations (on enlevait les 30 premières),

seuil=30
library(datasets)
X=as.numeric(Nile)
Xc=X[-(1:seuil)]
plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
lines(Xc,lwd=3)

On va alors commencer à la date 0 à partir de la trentième observation. On va translater la série pour avoir une série centrée sur 0 (dans un premier temps).

PP.test(Xc)
Y=Xc-mean(Xc)
plot(acf(Y))	

On regardant un peu la série, on rejette l'hypothèse de bruit blanc, et de cycle. Donc on se lance dans un modélisation par un processus http://latex.codecogs.com/gif.latex?ARMA(p,q).

pacf(Y,30)
model1=arima(Y,order=c(4,0,0))
model1
acf(residuals(model1))
Box.test(residuals(model1),lag=12,type='Box-Pierce')

On peut commencer par valider une modélisation http://latex.codecogs.com/gif.latex?AR(4). Ensuite,

acf(Y,30)
model2=arima(Y,order=c(0,0,4))
model2
acf(residuals(model2))
Box.test(residuals(model2),lag=12,type='Box-Pierce')

on va valider une modélisation http://latex.codecogs.com/gif.latex?MA(4). Enfin, on peut voir parmi la sélection automatique de modèle,

library(caschrono)
armaselect(Y,nbmod=5)

Le première modèle que l'on va tester, un processus http://latex.codecogs.com/gif.latex?MA(1) est acceptée.

model3=arima(Y,order=c(0,0,1))
model3
acf(residuals(model3))
Box.test(residuals(model3),lag=6,type='Box-Pierce')
 
BP=function(h) Box.test(residuals(model3),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")

En revanche les autres modèles testés sont rejetés.

model4=arima(Y,order=c(1,0,1))
model4
model5=arima(Y,order=c(2,0,1))
model5

On a trois modèles, et c'est déjà pas mal. Faisons un peu de backtesting pour comparer.

AIC(model1)
AIC(model2)
AIC(model3)
b=10
T=length(Y)
sY=Y[-((T-b):T)]
smodel1=arima(sY,order=c(4,0,0))
smodel2=arima(sY,order=c(0,0,4))
smodel3=arima(sY,order=c(0,0,1))
Yp1=forecast(model1,b+1)
Yp2=forecast(model2,b+1)
Yp3=forecast(model3,b+1)
M=cbind(Y[((T-b):T)],Yp1$mean,Yp2$mean,Yp3$mean)
sum((M[,1]-M[,2])^2)
sum((M[,1]-M[,3])^2)
sum((M[,1]-M[,4])^2)

Le premier modèle semble plus intéressant. On va donc le retenir pour faire de la prévision.

plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
abline(h=mean(Xc),col="green")
lines(Xc,lwd=3)
polygon(c(temps,rev(temps)),
c(Yhat$lower[,2],
rev(Yhat$upper[,2]))+mean(Xc),col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(Yhat$lower[,1],
rev(Yhat$upper[,1]))+mean(Xc),col="orange",border=NA)
lines(temps,Yhat$mean+mean(Xc),col="red")

Maintenant, on peut se demander s'il n'y avait pas une tendance linéaire (on va faire comme si, pour l'exercice)

plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
lines(Xc,lwd=3)
temps=1:length(Xc)-seuil
abline(lm(Xc~temps),lwd=2,col="red")

Dans ce cas, on modèlise la série à laquelle on a retranché cette tendance linéaire. On reprend tous les étapes précédentes et cette fois (pour changer un peu) on va faire comme si on retenait un modèle moyenne-mobile.

Y=Xc-predict(lm(Xc~temps))
plot(Y,type="l")
pacf(Y,30)
acf(Y,30)
model2=arima(Y,order=c(0,0,4))
model2
acf(residuals(model2))
Box.test(residuals(model2),lag=12,type='Box-Pierce')

Pour la prévision, on va projeter la série stationnaire et tout remonter, en tenant compte de la tendance,

Xc=X[-(1:seuil)]
tps=1:length(Xc)
base=data.frame(Xc,tps)
tendance=lm(Xc~tps,data=base)
T=length(Y)
temps=T+1:20
pred=predict(tendance,newdata=
data.frame(tps=c(temps,rev(temps))))
plot(1:length(X)-seuil,X,xlim=c(-seuil,90),type="l")
lines(Xc,lwd=3)
abline(tendance,col="green")
polygon(c(temps,rev(temps)),c(Yhat$lower[,2],
rev(Yhat$upper[,2]))+pred,col="yellow",border=NA)
polygon(c(temps,rev(temps)),c(Yhat$lower[,1],
rev(Yhat$upper[,1]))+pred,col="orange",border=NA)
lines(temps,Yhat$mean+pred[1:length(temps)],col="red")

L'étape suivante sera du faire du backtesting avec ces deux types de modèles, avec constante versus avec tendance... Mais je ne vais pas tout faire non plus...

Wednesday, March 21 2012

La maudite constante des modèles ARIMA

Dans les modèles ARIMA, autant le dire tout de suite, les constantes c'est pénible ! Pour comprendre un peu mieux ce qui se passe, considérons ici un modèle ARMA avec constante. Ou pour commencer, juste un AR(1)
Supposons que la série  soit stationnaire, de moyenne . Alors en prenant l’espérance de part et d'autre dans l'expression précédente, , i.e. , ou . Plus généralement, avec un ARMA, , en prenant la encore l’espérance des deux cotés, on en déduit .

Simulons un processus AR(1),

Pour simuler des processus ARIMA, on pourrait utiliser la commande

> X=arima.sim(list(order=c(1,0,0),ar=1/3),n=1000)+2
> mean(X)
[1] 1.931767
mais comme le but est de comprendre ce qui se passe, autant faire les choses calmement,
> X=rep(NA,1010)
> X[1]=0
> for(t in 2:1010){X[t]=4/3+X[t-1]/3+rnorm(1)}
> X=X[-(1:10)]
> mean(X)
[1] 2.03397
I.e. ici le processus est de moyenne qui vaut 2.


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

> arima(X, order = c(1, 0, 0))
 
Call:
arima(x = X, order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
0.3738     2.0334
s.e.  0.0294     0.0487
 
sigma^2 estimated as 0.9318:  log likelihood = -1383.68
De manière un peu surprenante, le coefficient appelé intercept n'est pas la constante dans le modèle AR(1), mais la moyenne du processus. Autrement dit, R n'ajuste pas un processus
comme le laisserait penser l'intuition, mais un processus
Ces deux formes sont bien entendu équivalentes. Mais les coefficients estimés ne sont pas tout à fait ce que l'on attendait...
Plaçons nous maintenant dans le cas d'un processus non-stationnaire. L'extension naturelle serait de considérer un processus ARIMA(1,1,0), ou bien un processus tel que la différence soit un processus AR(1). Un processus ARIMA(1,1,0) avec une constante s'écrirait
en utilisant l'opérateur retard. Ceci fait penser à un modèle avec une tendance linéaire. Posons  (afin de se débarrasser de cette tendance). Alors
i.e.
ou encore
On note que  ou encore  suit alors un processus ARIMA(1,1,0) sans constante cette fois.
Afin de visualiser ce que donnerait l'inférence, simulons le processus suivant, qui est un processus AR(1) avec constante que l'on intègre:  avec
Peut-on retrouver les différents paramètres du processus avec R ?
Commençons (là encore) par simuler un tel processus,
> U=rep(NA,1010)
> U[1]=0
> for(t in 2:1010){U[t]=4/3+U[t-1]/3+rnorm(1)}
> U=U[-(1:10)]
> X=cumsum(U)
Sous R, on obtient l'estimation suivant si l'on tente de calibrer un modèle ARIMA(1,1,0)
> arima(X, order = c(1, 1, 0))
 
Call:
arima(x = X, order = c(1, 1, 0))
 
Coefficients:
ar1
0.8616
s.e.  0.0160
  
sigma^2 estimated as1.343:log likelihood = -1565.63

Ça ne convient pas du tout.... On peut tenter un processus AR(1) (avec constante) sur la série différenciée...

> arima(diff(X), order = c(1, 0, 0))
 
Call:
arima(x = diff(X), order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
0.3564     2.0200
s.e.  0.0295     0.0486
 
sigma^2 estimated as 0.9782:  log likelihood = -1406.6
On progresse, sauf que comme auparavant, le terme qui est donne n'est pas la constante dans le modèle ARIMA, mais la moyenne du processus différencié. Mais cette fois, on a un interprétation, c'est que la constante est la pente de la tendance ! Si on estime la pente associée a , on récupère la même valeur:
> arima(X, order = c(1, 1, 0), xreg=1:length(X))
 
Call:
arima(x = X, order = c(1, 1, 0), xreg = 1:length(X))
 
Coefficients:
ar1  1:length(X)
0.3566       2.0519
s.e.  0.0296       0.0487
 
sigma^2 estimated as 0.9787:  log likelihood = -1406.82
Sur la figure ci-dessous, on retrouve le fait qu'en enlevant la tendance linéaire à  donneune série intégrée, sans constante (qui fait penser à une marche aléatoire)

Autrement dit

  • dans le cas d'une série stationnaire, la constante estimée n'est pas du tout la constante, mais la moyenne de la série temporelle
  • dans le cas d'une série non-stationnaire, la constante estimée dans la série différenciée a du sens, au sens ou il s'agit de la pente de la tendance (linéaire) du processus
Avant de conclure, une petite remarque. Quid de la prévision ? Si on commencer par une révision à l'aide du premier processus ARIMA, on obtient une prédiction pour avec un énorme intervalle de confiance (sans aucun bon sens)
> ARIMA1=arima(X, order = c(1, 1, 0))
> ARIMA2=arima(X, order = c(1, 1, 0), xreg=1:length(X))
> Xp1=predict(ARIMA1,20)
> Xp2=predict(ARIMA2,20,newxreg=
+ (length(X)+1):(length(X)+20))
> plot(960:1000,X[960:1000],xlim=c(960,1020),type="l")
> polygon(c(1001:1020,rev(1001:1020)),
+ c(Xp1$pred+2*Xp1$se,rev(Xp1$pred-2*Xp1$se)),
+ col=CL[3],border=NA)
> lines(1001:1020,Xp1$pred,col="red",lwd=2)
intervalle qui se visualise sur le graphique ci-dessous,

Si on regarde pour l'autre modèle,
> lines(1001:1020,Xp2$pred,col="blue",lwd=2)

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

Friday, August 6 2010

Quand surviennent les accidents de la route ? (partie 2)

Je poursuis un peu mon étude des accidents de la route (ou pour être plus précis, des accidents corporels, ayant causé des blessés, et ayant fait l'objet d'un rapport de police).

J'avais parlé ici de la tendance de long terme, sur 6 ou 7 ans, montrant qu'il y avait une tendance à la baisse aussi bien du nombre de tués que du nombre de blessés grave. On avait observé que les cycles annuels avaient tendance à diminuer également. Mais qu'en est-il des cycles à court terme ?
En particulier, on peut se demander s'il n'existe pas des cycles dans la journée. En semaine, l'hiver (entre novembre et début février), on observe les tendances suivantes,

(les points sont les moyennes brutes, par tranche d'une demi heure, et la courbe est un lissage de ces points). On distinguera les années 2002-2004 (en rouge) et 2004-2007 (en bleu) Si l'on regarde les semaines l'été (entre mai et début août), on a

Notons que si l'on compte non plus les accidents, mais les blessés graves ou les tués (il y a généralement plusieurs victimes dans ces accidents), on observe les cycles suivants,

pour l'hiver, alors que l'été,

Autrement dit, les tendances sont les mêmes: les accidents du matin ne sont pas moins meutriers que ceux du soir... Les week-ends, les tendances sont assez différentes entre l'hiver et l'été, avec respectivement

et

Si l'on regarde les blessés graves et les tués, on a les cycles suivants,

et