Freakonometrics

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 .

> 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='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 .

> 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")

Suite aux (nombreuses) questions de ce matin, quelques éléments de réponse...

• l'estimation ne "converge pas"

quelques exemples d'erreurs qui peuvent apparaître lors de l'estimation...  Car quand on travaille sur une série non-stationnaire (ou presque) et que l'on ajuste un ARIMA, on peut avoir des soucis...

> arima(cumsum(rnorm(100)),c(1,0,0))
 
Call:
arima(x = cumsum(rnorm(100)), order = c(1, 0, 0))
 
Coefficients:
ar1  intercept
        1     0.5096
s.e.  NaN  1302.3084
 
sigma^2 estimated as 0.8942:
log likelihood = -136.81,  aic = 279.61
Message d'avis :
In sqrt(diag(x$var.coef)) : production de NaN

cette erreur n'est pas trop gênante: il n'arrive pas à calculer l'écart-type du coefficient d'autorégression, qui vaut 1. Bref, on a une racine unité, et le mieux (après avoir fait un test de racine unité pour confirmer) est alors de différencier. Maintenant, on peut avoir une erreur un peu plus vicieuse,

> arima(cumsum(rnorm(100)),c(5,0,0))
Erreur dans arima(cumsum(rnorm(100)), c(5, 0, 0)) :
partie AR non stationnaire dans CSS

Il semble dire que l'on a une racine unité, mais on ne voit rien... Une stratégie peut être de changer la méthode d'estimation (qui, encore une fois, est une méthode numérique d'optimisation), par exemple, on peut se contenter des moindres carrés conditionnels (par défaut, on commence par les moindres carrés, puis une routine de maximisation de la vraisemblance est lancée),

> arima(cumsum(rnorm(100)),c(5,0,0),method="CSS")
 
Call:
arima(x=cumsum(rnorm(100)), order=c(5,0,0), method="CSS")
 
Coefficients:
ar1     ar2     ar3     ar4      ar5  intercept
      0.7468  0.1532  0.0244  0.1836  -0.1255    -8.4584
s.e.  0.0997  0.1231  0.1247  0.1219   0.1003     8.3181
 
sigma^2 estimated as 1.153:  part log likelihood = -149

Normalement, cette méthode permet d'avoir des estimateurs même si on nous encourage à aller faire un test de racine unité...

• la fonction H2M ne marche pas

> plot(X)
Erreur dn[[2]]

H2M=function(BASE){
X=BASE[,2]
Y=BASE[,1]
date1=substr(as.character(Y),1,10)
date2=substr(as.character(Y),14,23)
D1=as.Date(date1,"%Y-%m-%d")
D2=as.Date(date2,"%Y-%m-%d")
vm=vy=N=NA
for(t in 1:length(D1)){
mois=seq(D1[t],D2[t],length=7)
vm=c(vm,as.POSIXlt(mois)$mon+1)
vy=c(vy,as.POSIXlt(mois)$year+1900)
N=c(N,rep(X[t],7))}
N=N[-1]; vm=vm[-1]; vy=vy[-1]
YM=vy*100+vm
Z=tapply(N,as.factor(YM),mean)
Zts=as.ts(as.numeric(Z),start=c(2004,1),frequency=12)
return(Zts)}

> X=as.numeric(H2M(base))
> temps=seq(2004,length=length(X),by=1/12)
> plot(temps,X)

• le fichier de Google est en anglais...

> base=read.table("trends2.csv",skip=5,
+ header=FALSE,nrows=458,sep=";")
> head(base)
V1   V2   V3
1  Jan 4 2004 0.00 >10%
2 Jan 11 2004 1.20 >10%
3 Jan 18 2004 1.40 >10%
4 Jan 25 2004 1.64 >10%
5  Feb 1 2004 1.40 >10%
6  Feb 8 2004 1.27 >10%

H2M=function(BASE){
Mois = c("Jan","Fév", "Mar", "Avr",
"Mai", "Jui", "Jul", "Aoû", "Sep", "Oct",
"Nov", "Déc")
Months = c("Jan", "Feb", "Mar", "Apr",
"May", "Jun", "Jul", "Aug",
"Sep", "Oct", "Nov", "Dec")
X=BASE[,2]
jour=BASE[,1]
D1=as.Date(as.character(jour),"%b %d %Y")
if(sum(is.na(D1)>0)){
moisGB=substr(as.character(jour),1,3)
V=1:12; names(V)=Months
NoMois=as.numeric(V[moisGB])
moisFR=Mois[NoMois]
jourFR=paste(moisFR,substr(as.character(
jour),4,nchar(as.character(jour))),sep="")
D1=as.Date(jourFR,"%b %d %Y")
}
D2=D1+6
vm=vy=N=NA
for(t in 1:length(D1)){
mois=seq(D1[t],D2[t],length=7)
vm=c(vm,as.POSIXlt(mois)$mon+1)
vy=c(vy,as.POSIXlt(mois)$year+1900)
N=c(N,rep(X[t],7))}
N=N[-1]; vm=vm[-1]; vy=vy[-1]
YM=vy*100+vm
Z=tapply(N,as.factor(YM),mean)
Zts=ts(as.numeric(Z),start=c(2004,1),frequency=12)
return(Zts)}

• un peu de lecture pour finir...

Un billet pour reprendre (et compléter) le code vu hier, en cours de modèles de prévision. On a toujours pour but de prévoir l'évolution du trafic autoroutier, en France, sur l'autoroute A7.

> autoroute=read.table(
+ "http://freakonometrics.blog.free.fr/public/data/autoroute.csv",
+ header=TRUE,sep=";")
> a7=autoroute$a007
> A7=ts(a7,start = c(1989, 9), frequency = 12)

Comme je l'ai dit en cours, les constantes dans les modèles ARIMA, a complexifie inutilement les notations, donc pour faire des choses proches de ce qu'ai pu faire dans le cours, on va centrer la série, et considérer

> A7d=A7-mean(A7)
> plot(A7d,xlim=c(1989,2000))

On avait vu qu'un modèle pouvait convenir, donc on va l'utiliser pour faire de la prévision.

>  fit=arima(A7d,order=c(p=12,d=0,q=0),
+      include.mean =FALSE)
>  summary(fit)
Series: A7d
ARIMA(12,0,0) with zero mean
 
Coefficients:
ar1      ar2      ar3     ar4      ar5     ar6
      0.0981  -0.0353  -0.0375  0.0530  -0.0934  0.0598
s.e.  0.0630   0.0623   0.0600  0.0606   0.0618  0.0568
ar7     ar8     ar9     ar10    ar11    ar12
      -0.0839  0.0278  0.0072  -0.1080  0.1578  0.7887
s.e.   0.0619  0.0629  0.0627   0.0629  0.0637  0.0623
 
sigma^2 estimated as 12823951:  log likelihood=-827.24
AIC=1678.48   AICc=1683.6   BIC=1710.23
 
In-sample error measures:
ME         RMSE          MAE          MPE
257.8363135 3581.0544202 2335.0686684    1.7258390
MAPE         MASE
34.7942409    0.2603139 

La fonction permettant de faire de la prévision est alors

> horizon=48
> library(forecast)
> A7hat=forecast(fit,h=horizon)

On a alors un ensemble de séries temporelles, la première étant la prévision proprement dite , i.e. , où désigne une espérance linéaire (ou en terme technique, une projection sur l'adhérence de l'ensemble des des combinaisons linéaires du passé). On a ensuite les intervalles de confiance estimés, à partir d'une hypothèse de normalité, en utilise un estimateur de ,

> A7hat
Point    Lo 80   Hi 80    Lo 95    Hi 95
Oct 1996    -7276.44 -11865.7 -2687.1 -14295.1  -257.70
Nov 1996   -10453.79 -15065.1 -5842.4 -17506.1 -3401.39
Dec 1996    -6970.67 -11583.5 -2357.8 -14025.3    84.02
Jan 1997   -13404.09 -18021.2 -8786.9 -20465.4 -6342.78

On peut alors représenter ces valeurs graphiquement, comme on l'a fait en cours

>  (TPS=tsp(A7hat$mean))
[1] 1996.750 2000.667   12.000
>  lines(A7hat$mean,col="red",lwd=2)
>  temps=seq(from=TPS[1],to=TPS[2],by=1/TPS[3])
>  lines(temps,A7hat$lower[,1],col="blue",lwd=2)
>  lines(temps,A7hat$lower[,2],col="purple",lwd=2)
>  lines(temps,A7hat$upper[,1],col="blue",lwd=2)
>  lines(temps,A7hat$upper[,2],col="purple",lwd=2)
>  abline(v=TPS[1]-.5/TPS[3],col="green")

En allant un peu plus loin, on peut faire plus joli, en utilisant des aires colorées afin d'illustrer une région de confiance

> plot(A7d,xlim=c(1989,2000))
> polygon(c(temps,rev(temps)),c(A7hat$lower[,2],
+ rev(A7hat$upper[,2])),col="yellow",border=NA)
> polygon(c(temps,rev(temps)),c(A7hat$lower[,1],
+ rev(A7hat$upper[,1])),col="orange",border=NA)
> lines(temps,A7hat$mean,col="red")
> abline(v=TPS[1]-.5/TPS[3],col="green")

Maintenant, si on revient sur notre série initiale, l'autocorrélogramme pourrait laisser penser qu'il y a une racine unité saisonnière,

> plot(acf(A7d,lag=36))

Donc on peut tenter de différencier à l'ordre 12,

> plot(acf(diff(A7d,12),lag=36))

On pourrait avoir l'impression que la série ainsi différenciée est un bruit blanc. On peut faire un test de Box-Pierce pour tester cette hypothèse,

> BP=function(h) Box.test(diff(A7d,12),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")

On pourrait valider ici l'hypothèse de bruit blanc. On va donc ajuster un (pur) modèle ARIMA saisonnier,

> fit2=arima(A7d, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

Si on compare les deux modèles,

> AIC(fit1)
[1] 1511.815
> AIC(fit2)
[1] 1245.374

Le second modèle ayant un critère d'Akaike plus faible, on aurait tendance  à préférer ce modèle au premier. Mais cette méthode de comparaison n'est pas forcément la plus pertinent, si l'objet est de faire une prévision à court terme (disons quelques mois). Pour cela, on va plutôt enlever les dernières observations, ajuster un modèle sur la série restreinte, et comparer les prévisions, avec ce qui a été réellement observé. Disons qu'on va faire du backtesting sur les 9 derniers mois,

> T=length(A7d)
> backtest=9
> subA7d=ts(A7d[-((T-backtest+1):T)],
+ start = c(1989, 9), frequency = 12)
> fit1s=arima(subA7d,order=c(p=12,d=0,q=0),
+ include.mean =FALSE)

pour le premier modèle, et

> fit2s=arima(subA7d, order = c(p=0, d=0, q=0),
+ seasonal = list(order = c(0, 1, 0), period = 12))

pour le second. Ensuite, on peut faire de la prévision, et comparer les prévisions avec les données écartées pour l'estimation du modèle.,

> XA7p1=forecast(fit1s,h=backtest)
> XA7p2=forecast(fit2s,h=backtest)
> A7dobs=A7d[((T-backtest+1):T)]
> (M=cbind(observé=A7dobs,modèle1=as.numeric(XA7p1$mean),
+ modèle2=as.numeric(XA7p2$mean)))
observé     modèle1    modèle2
[1,] -13610.071 -12264.7138 -13558.071
[2,] -10108.071 -12236.6227 -10410.071
[3,]  -8034.071  -9333.8365  -9504.071
[4,]   3899.929   4071.9313   4108.929
[5,]   2009.929    901.3745   2661.929
[6,]   5253.929   7610.9954   5866.929
[7,]  25357.929  29232.2317  30317.929
[8,]  35229.929  29641.5172  31737.929
[9,]   4341.929   5603.8008   4940.929

On peut alors retenir un critère classique, comme l'erreur quadratique, et donc comparer la somme des carrées des erreurs,

> sum((M[,1]-M[,2])^2)
[1] 62677238
> sum((M[,1]-M[,3])^2)
[1] 40253827

Le second modèle semble ici avoir été meilleur, sur les 9 derniers mois. On peut aussi visualiser les erreurs de prévision,

> (TPS=tsp(subA7d))
[1] 1989.667 1995.917   12.000
> temps=seq(TPS[2]+1/TPS[3],by=1/TPS[3],length=backtest)
> rangex=c(TPS[1],TPS[2]+backtest/TPS[3])
> plot(subA7d,type="l",xlim=rangex)
> abline(v=T-11.5,lty=2)
> polygon(c(temps,rev(temps)),c(A7dobs,rev(XA7p1$mean)),
+ col="yellow",border=NA)
> lines(temps,A7dobs)
> lines(temps,XA7p1$mean,col="red")

pour le premier modèle,

et

> polygon(c(temps,rev(temps)),c(A7dobs,rev(XA7p2$mean)),
+ col="yellow",border=NA)
> lines(temps,XA7p2$mean,col="red")

pour le second,

Histoire de presque finir la modélisation des processus ARMA, un mot sur les tests de bruit blanc, basé sur ce que nous avons pu voire en cours la semaine passée. Poursuivons à partir du billet précédant, sur la série du trafic autoroutier, ou plus précisément la série

> a7=as.numeric(A7)

qui est un simple vecteur, et non plus une série temporel, au sens défini par R. On peut tenter une modélisation , 

> fit=arima(a7,order=c(12,0,0))
> summary(fit)
Series: a7
ARIMA(12,0,0) with non-zero mean
 
Coefficients:
ar1      ar2      ar3     ar4      ar5
      0.0979  -0.0355  -0.0377  0.0528  -0.0935
s.e.  0.0630   0.0623   0.0600  0.0606   0.0618
ar6      ar7     ar8     ar9     ar10
      0.0596  -0.0839  0.0276  0.0072  -0.1081
s.e.  0.0568   0.0618  0.0629  0.0627   0.0629
ar11    ar12  intercept
      0.1577  0.7885  38653.234
s.e.  0.0637  0.0623   1507.591
 
sigma^2 estimated as 12820052:  log likelihood=-827.22
AIC=1680.44   AICc=1686.44   BIC=1714.64
 
In-sample error measures:
ME         RMSE          MAE
3.152304e+02 3.580510e+03 2.334799e+03
MPE         MAPE
3.321677e-02 6.232293e+00 

et on peut s'intéresser un peu à la série des autocorrélations de la série des résidus empiriques , obtenus par partir de l'estimation précédente,

> plot(acf(residuals(fit)))

Les lignes horizontales sont les valeurs critiques d'un test contre , à donné. Ce n'est pas un test de bruit blanc, qui cherche à tester si toutes les autocorrélations sont nulles,

Naturellement, on a envie d'utiliser la distance euclidienne, pour voir si la distance est petit, i.e. quelque chose du genre

C'est la statistique proposée par Box-Pierce,

> rho=as.vector(acf(residuals(fit))$acf)
> (Q=length(residuals(fit))*sum(rho[1:6]^2))
[1] 5.27944

On peut l'avoir automatiquement à l'aide de la fonction

> Box.test(residuals(fit),lag=6,
+ type='Box-Pierce')
data:  residuals(fit)
X-squared = 5.2794, df = 6, p-value = 0.5085

Sous l'hypothèse nulle, la statistique doit suivre une loi du chi-deux, d'où la p-value renvoyée ci-dessus,

> 1-pchisq(Q,df=6)
[1] 0.5085045

Une autre statistique possible est celle de Ljung-Box, qui tient compte du fait que les autocorrélations n'ont pas été calculées sur le même nombre d'observations (d'où la petite correction,  qui devient importante si on a peu de données),

> Box.test(residuals(fit),lag=6,
+ type='Ljung-Box')
data:  residuals(fit)
X-squared = 5.5278, df = 6, p-value = 0.4781

Mais le test est fait ici pour un nombre de retard fixé. Il convient de le faire pour plusieurs retards possible. Ou mieux, de faire un petit dessin, permettant de visualiser la p-value pour différents retards,

> BP=function(h) Box.test(residuals(fit),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")

On valide ici l'hypothèse de bruit blanc pour notre série. Si on avait tenté sur un modèle ,

> fit=arima(a7,order=c(7,0,0))

cette fois, en revanche, on ne peut accepter l'hypothèse de bruit blanc sur la série des résidus induits,

Un billet rapide pour reprendre le code tapé en cours, la semaine passée. Considérons  un processus autorégressif d'ordre 1,  où  est un bruit blanc, stationnaire, i.e.  appartient à l'intervalle . Le code pour simuler un tel processus est

n=1000
bruit=rnorm(n)
phi1= .85
X=rep(NA,n)
X[1]=0
for(t in 2:n){X[t]=phi1*X[t-1]+bruit[t]}
plot(acf(X),lwd=5,col='blue')
plot(pacf(X),lwd=5,col='blue')

ou avec un autocorrélation au premier ordre négative,

phi1= -0.7

On peut aussi regarder un processus autorégressif au second ordre, 

sur la figure ci-dessous (avec en haut à gauche le triangle de stationnarité du couple de paramètres).

phi1=  0.3
phi2=  0.5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=phi1*X[t-1]+phi2*X[t-2]+bruit[t]}

Histoire de changer un peu, on peut regarder un processus moyenne mobile au premier ordre,  où  est un paramètre dans .

theta1=  .8
X=rep(NA,n)
X[1]=0
for(t in 2:n){
X[t]=bruit[t]+theta1*bruit[t-1]}

ou une moyenne mobile du second ordre,

theta1= -.6
theta2=  .5
X=rep(NA,n)
X[1:2]=0
for(t in 3:n){
X[t]=bruit[t]+theta1*bruit[t-1]+
theta2*bruit[t-2]}