Freakonometrics

To content | To menu | To search

Teaching › ACT6420-H2012

Entries feed - Comments feed

Friday, September 28 2012

Fin de session !

La session d'hiver est enfin terminée (pour les cours en tous cas). L'examen du cours ACT6420 méthodes de prévisions est maintenant en ligne, avec l'énoncé, les annexes (qui correspondent à celles mises en ligne en début de semaine, avec quelques nombres entourés, qu'il fallait commenter expliquer) et quelques éléments (succincts) de correction (qui seront en ligne à partir de 15 heures).

Il reste encore le deuxième devoir à rendre, pour vendredi prochain dernier délai (comme indiqué mardi, après le cours). Merci aux étudiants qui ont joué le jeu, et qui ont mis une prévision de leur propre note.

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

Monday, September 24 2012

Séries temporelles, dernier cours, cas d'école

Quelques séries temporelles (mensuelles pour la plupart) pour le dernier cours,

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/boston-robberies.csv",
header=TRUE,sep=",",nrows=118)
Xt=ts(X[,2],start=c(1966,1),frequency=12)

Monthly Boston armed robberies Jan.1966-Oct.1975, via http://datamarket.com/

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/ozone-arosa.csv",
header=TRUE,sep=";",nrows=480)
Xt=ts(X[,2],start=1932,frequency=12)

Ozone, Arosa, 1932-1972, via http://datamarket.com

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/uk-deaths.csv",
header=TRUE,sep=";",nrows=204)
Xt=ts(X[,2],start=c(1974,1),frequency=12)

U.K. deaths from bronchitis, emphysema and asthma, via http://datamarket.com/

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)

Monthly car sales in Quebec 1960-1968, via http://datamarket.com/

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/auto-registration.csv",
header=TRUE,sep=",",nrows=118)
Xt=as.ts(X[,2],start=1966,frequency=12)

Monthly U.S. auto registration (thousands) 1947 – 1968, via http://datamarket.com/

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/traffic-fatalities-ontario.csv",
header=TRUE,sep=";",nrows=180)
Xt=ts(X[,2],start=c(1960,1),frequency=12)

Monthly traffic fatalities in Ontario 1960-1974, via http://datamarket.com/

X=read.table(
"http://freakonometrics.blog.free.fr/public/
data/airline.csv",
header=TRUE,sep=",",nrows=118)
Xt=ts(X[,2],start=c(1966,1),frequency=12)

International airline passengers: monthly totals in thousands. Jan 49 – Dec 60, via http://datamarket.com/

Examen final ACT6420, énoncé

L'examen ACT6420 aura lieu vendredi midi, et durera trois heures. Les Annexes (contenant 23 pages de sorties informatiques) sont en ligne ici. La moitié des questions porteront sur l'analyse de ces sorties.

Friday, September 21 2012

Transformation logarithmique de séries temporelles

Pour poursuivre une discussion amorcée en fin de cours, dans certains cas, on peut avoir l'impression que modéliser une série pourrait être compliqué,

plot(X,xlim=c(1,length(X)+20))

Mais on peut avoir l'intuition que modéliser le logarithme de la série pourrait être plus simple,

> X=log(Y)
> plot(X,xlim=c(1,length(X)+20))

On va alors tenter une modélisation par un processus ARMA de cette dernière série,

> md=arima(X,c(12,0,1))
> P=predict(md,24)
> E=P$pred
> V=P$se^2

On peut alors faire une prévision sur cette série plus simple à modéliser, et visualiser cette prévision.

> temps=length(X)+1:24
> ciu=(E+2*sqrt(V))
> cil=(E-2*sqrt(V))
> polygon(c(temps,rev(temps)),c(ciu,
+ rev(cil)),col="yellow",border=NA)
> lines(temps,E,col="red",lwd=2)
> lines(temps,ciu,col="red",lty=2)
> lines(temps,cil,col="red",lty=2)

Maintenant, on va devoir remonter. On va utiliser un résultat que l'on a vu sur la transformation logarithmique dans une régression: si après avoir pris le logarithme, on a un modèle simple, Gaussien, c'est que le modèle initial était log-normal. On peut alors utiliser les propriétés de la loi lognormale, dont on connait les moments à partir de ceux de la loi Gaussienne sous-jacente. Pour la prévision, on n'a pas trop le choix,

> mu=exp(E+.5*V)

Par contre, pour construire un intervalle de confiance, soit on utilise la variance de notre loi lognormale pour avoir la variance de notre processus, et on oublie cette histoire de loi lognormale pour construire un intervalle Gaussien,

> sig2=(exp(V)-1)*exp(2*E+V)
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)

ou alors on utilise le fait que comme la transformation est monotone, l'intervalle de confiance peut etre vu comme une transformation du précédant intervalle de confiance,

> ci2u=exp(E+2*sqrt(V))
> ci2l=exp(E-2*sqrt(V))

Si on compare visuellement les deux, on a dans le premier cas,

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> polygon(c(temps,rev(temps)),c(ci1u,
> rev(ci1l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci1u,col="red",lty=2)
> lines(temps,ci1l,col="red",lty=2)

(qui est symétrique et centré sur notre prévision) et dans le second

> plot(Y,xlim=c(1,length(X)+20))
> temps=length(X)+1:24
> ci1u=mu+2*sqrt(sig2)
> ci1l=mu-2*sqrt(sig2)
> polygon(c(temps,rev(temps)),c(ci2u,
> rev(ci2l)),col="yellow",border=NA)
> lines(temps,mu,col="red",lwd=2)
> lines(temps,ci2u,col="red",lty=2)
> lines(temps,ci2l,col="red",lty=2)

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

Thursday, September 20 2012

FAQ prévision et ARIMA avec R

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
Il semble que la fonction que j'avais mise en ligne dans un précédant billet, pour transformer les données hebdomadaire ne fonctionne pas... sous windows. En fait, la fonction marche, il semble que l'on récupère une série temporelle, mais on récupère un message d'erreur
> plot(X)
Erreur dn[[2]]
N'ayant pas windows, je vais avancer à tâtons... La fonction suivante devrait faire la job
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)}
Dans le pire des cas, il est possible de récupérer une série sur laquelle travailler en utilisant
> X=as.numeric(H2M(base))
> temps=seq(2004,length=length(X),by=1/12)
> plot(temps,X)
Le vecteur contient toutes les informations pour faire la modélisation, c'est juste que quand on fera un graphique, en abscisse on aura l'indice du vecteur, et pas l'année (je renvoie à la discussion en cours sur ce sujet).

  • le fichier de Google est en anglais...
Là c'est un peu plus vicieux... mais on va y arriver, il suffit de changer un peu la fonction que l'on va utiliser... Mais tout d'abord, je vais me contredire un peu: on commencer par ouvrir la base sous excel, et on enregistre en csv (avec un séparateur de colonne en point virgule),
> 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%
Et cette fois, on modifie un peu la fonction, car la date a un format différent,
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)}
Normalement, ce code marche...

  • un peu de lecture pour finir...
Suite au cours de mardi dernier, un petit complément sur les suites définies par une récurrence (linéaire) à l'ordre 2, ou je peux renvoyer vers quelques notes de Djamal Rebaïne.

Wednesday, September 19 2012

Prévision, choix de modèle(s) et backtesting

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 http://latex.codecogs.com/gif.latex?X_t-\overline%20X

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

On avait vu qu'un modèle http://latex.codecogs.com/gif.latex?AR(12) 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 http://latex.codecogs.com/gif.latex?{}_T\widehat{X}_{T+h}, i.e. http://latex.codecogs.com/gif.latex?EL(X_{T+h}|\underline{X}_T), où http://latex.codecogs.com/gif.latex?EL(\cdot) 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 http://latex.codecogs.com/gif.latex?{}_T\text{mse}_{T+h}=\mathbb{E}([{}_T\widehat{X}_{T+h}-X_{T+h}]^2|\underline{X}_T),

> 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,

Tuesday, September 18 2012

Examen final ACT6420

Chose promise, chose due: les deux bases de données qui seront utilisées pour l'examen sont désormais en ligne. 10 questions porteront sur l'analyse de la série temporelle suivante

serie=read.table("http://freakonometrics.free.fr/exam-ts.txt")
X=serie$x
XTS=serie=ts(X,start=c(1999,9),frequency=12)

(qui peut être modélisée soit comme un objet série temporel, soit comme un vecteur numérique).

Enfin, 10 questions porteront sur la base de données suivante

examen=read.table(
"http://freakonometrics.blog.free.fr/public/
data/basket-exam-v2.csv"
, header=TRUE,sep=";")

(si l'importation ne marche pas, veillez à ce que l'adresse du fichier tienne sur une ligne). Cette base devait servir initialement pour l'examen intra. Elle contient des matchs de baskets universitaires (aux États-Unis). Les variables sont les suivantes

  • EQUIPE le nom de l'équipe qui sert de référence
  • HALFDIFF le nombre de points de différence à la mi-temps (négatif signifie que l'EQUIPE perdait)
  • WINNER vaut 1 si l'équipe a gagné le match, 0 sinon,
  • TOTALPOINTSDIFF le nombre de points de différence à la fin du match
  • SECONDHALFPOINTS est le nombre de points marqués par l'EQUIPE pendant la seconde mi-temps
  • NODIFF_PTS_1000-LEFT la différence de points entre les équipes 10 minutes avant la fin du match (au milieu de la seconde période)
  • HOMETEAM vaut 0 si l'équipe joue à l'extérieur, 1 si elle joue à domicile
  • HALFWINNER vaut 1 si l'EQUIPE gagnait (supérieur ou nul) à la mi-temps
  • HALFLOOSER vaut 1 si l'EQUIPE perdait (inférieur ou nul) à la mi-temps
  • SECONDHALFPOINTSOTHER est le nombre de points marqués par l'adversaire au cours de la seconde mi-temps

La variable d'intérêt sera le nombre de points d'écart à la fin du match.

Enfin, 20 questions supplémentaires porteront sur la compréhension générale du cours. Les calculatrices seront autorisées. Maintenant, pour tous ceux qui ont encore du mal avec R, je recommande

  • “R pour les débutants”d'Emmanuel Paradis (PDF)
  • “Brise Glace-R” d'Andrew Robinson er Arnaud Schloesing (PDF)
  • “Introduction à la programmation en R” de Vincent Goulet (PDF)

(les questions d'analyse de sortie ne porteront pas sur la programmation en R, mais davantage sur l'analyse critique des sorties, et leur interprétation).

Monday, September 17 2012

Les tests de bruit-blanc

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 http://freakonometrics.blog.free.fr/public/perso6/ar12.gif

> 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 http://freakonometrics.blog.free.fr/public/perso6/epsilonhat.gif, obtenus par partir de l'estimation précédente,

http://freakonometrics.blog.free.fr/public/perso6/hatepsilon.gif
> plot(acf(residuals(fit)))

Les lignes horizontales sont les valeurs critiques d'un test http://freakonometrics.blog.free.fr/public/perso6/rhotest1.gif contre http://freakonometrics.blog.free.fr/public/perso6/rhotest2.gif, à http://freakonometrics.blog.free.fr/public/perso6/rhotest3.gif donné. Ce n'est pas un test de bruit blanc, qui cherche à tester si toutes les autocorrélations sont nulles,

http://freakonometrics.blog.free.fr/public/perso6/rhotest4.gif

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

http://freakonometrics.blog.free.fr/public/perso6/rhotest6.gif

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

http://freakonometrics.blog.free.fr/public/perso6/box.gif
> 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),

http://freakonometrics.blog.free.fr/public/perso6/box2.gif
> 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 http://latex.codecogs.com/gif.latex?AR(7),

> 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,

Ordres d'un processus ARMA

Dans la méthodologie de Box & Jenkins, une étape qui arrive très rapidement est le choix des ordres du processus http://freakonometrics.blog.free.fr/public/perso6/ARMA__latex.gif, une fois que l'on validé l'hypothèse de stationnarité de la série, comme on l'a vu en cours la semaine passée. Considérons la série du trafic autoroutier,

> 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)
Un outils pratique de sélections des ordres http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif dans un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMAlatex.gif est la fonction d'autocorrélation étendue. La définition est donnée dans les notes de cours (Def. 223) à partir des statistiques proposées par Tsay & Tiao (1984),
>  EACF=eacf(A7)
AR/MA
  0 1 2 3 4 5 6 7 8 9 10 11 12 13
0 x o o x x x x x o o x  x  x  o
1 x x o o x x x o o o x  x  x  x
2 o x o o o o x o o o o  x  x  o
3 o x o x o o o o o o o  x  x  x
4 x x x x o o o o o o o  x  o  o
5 x x o o o x o o o o o  x  o  x
6 x x o o o x o o o o o  x  o  o
7 o x o o o x o o o o o  x  o  o
>  EACF
$eacf
           [,1]       [,2]        [,3]        [,4]
[1,]  0.6476234  0.2124105 -0.02413173 -0.24234535
[2,]  0.4889076  0.2797152 -0.02494135 -0.06094037
[3,]  0.1006541 -0.2285771  0.03514148  0.08000588
[4,]  0.1390240 -0.2788742  0.04386746  0.28307260
[5,] -0.5091680  0.3144899 -0.34572269  0.31450865
[6,] -0.4224571 -0.4877505  0.16054232 -0.09130728
[7,] -0.4731353 -0.4324857 -0.04847184 -0.10500350
[8,] -0.2129591 -0.4072901  0.09487899 -0.06493243
             [,5]        [,6]        [,7]        [,8]
[1,] -0.514330187 -0.61634046 -0.52314403 -0.28008661
[2,] -0.254912957 -0.28966664 -0.33963243 -0.21863077
[3,] -0.156624357 -0.01199786 -0.25116738  0.13079231
[4,] -0.183283544  0.03651508 -0.08711829  0.11626377
[5,] -0.190885091 -0.09786463 -0.09182557  0.08818875
[6,] -0.072904071  0.29271777 -0.09334712  0.01972648
[7,] -0.009873289  0.36909726  0.01698660 -0.03317456
[8,]  0.020485930  0.38342158  0.16981715 -0.02592442
             [,9]       [,10]       [,11]     [,12]
[1,] -0.082607058  0.15178887  0.56583179 0.8368975
[2,] -0.085026400  0.02731460  0.26158357 0.7844748
[3,]  0.001866994 -0.11658312  0.03038621 0.7026361
[4,]  0.025183793 -0.21608692  0.05660781 0.6674301
[5,] -0.006831894 -0.02514440 -0.07390257 0.4762563
[6,]  0.010058718 -0.03888613 -0.04382043 0.5338091
[7,]  0.032124905 -0.07022090 -0.04427400 0.4674165
[8,]  0.024189179 -0.20818201  0.01459933 0.4830369
          [,13]       [,14]
[1,]  0.5637439  0.17862571
[2,]  0.4530716  0.24413569
[3,]  0.2534178 -0.20160890
[4,]  0.2409861 -0.29462510
[5,] -0.1517324 -0.14763294
[6,] -0.1701182 -0.28771495
[7,] -0.1651515  0.05466457
[8,] -0.1403198 -0.04095030
 
$ar.max
[1] 8
 
$ma.ma
[1] 14
Les ronds dans la matrice désignent des valeurs non-significatives. Par défaut, le nombre de retards, pris en compte pour la composante autorégressive est faible, mais on peut l'augmenter.
> EACF=eacf(A7,13,13)
AR/MA
   0 1 2 3 4 5 6 7 8 9 10 11 12 13
0  x o o x x x x x o o x  x  x  o
1  x x o o x x x o o o x  x  x  x
2  o x o o o o x o o o o  x  x  o
3  o x o x o o o o o o o  x  x  x
4  x x x x o o o o o o o  x  o  o
5  x x o o o x o o o o o  x  o  x
6  x x o o o x o o o o o  x  o  o
7  o x o o o x o o o o o  x  o  o
8  o x o o o x o o o o o  x  o  o
9  x x o o x o x o o o o  x  o  o
10 x x o o x o x o o o o  x  o  o
11 x x o o x x x o o o x  x  o  o
12 x x o o o o o o o o o  o  o  o
13 x x o o o o o o o o o  o  o  o
> EACF$eacf
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
[1,]  0.64  0.21 -0.02 -0.24 -0.51 -0.61 -0.52 -0.28 -0.08  0.15
[2,]  0.48  0.27 -0.02 -0.06 -0.25 -0.28 -0.33 -0.21 -0.08  0.02
[3,]  0.10 -0.22  0.03  0.08 -0.15 -0.01 -0.25  0.13  0.00 -0.11
[4,]  0.13 -0.27  0.04  0.28 -0.18  0.03 -0.08  0.11  0.02 -0.21
[5,] -0.50  0.31 -0.34  0.31 -0.19 -0.09 -0.09  0.08  0.00 -0.02
[6,] -0.42 -0.48  0.16 -0.09 -0.07  0.29 -0.09  0.01  0.01 -0.03
[7,] -0.47 -0.43 -0.04 -0.10  0.00  0.36  0.01 -0.03  0.03 -0.07
[8,] -0.21 -0.40  0.09 -0.06  0.02  0.38  0.16 -0.02  0.02 -0.20
[9,] -0.14 -0.50 -0.10 -0.06  0.11  0.38 -0.01 -0.02 -0.08 -0.20
[10,] -0.29  0.48  0.00  0.18  0.27 -0.12  0.41  0.00  0.16  0.15
[11,]  0.24  0.48 -0.04  0.02  0.46  0.10  0.37 -0.10  0.05  0.16
[12,] -0.59  0.49 -0.18  0.05  0.31 -0.32  0.34 -0.16  0.07  0.16
[13,] -0.31  0.31 -0.12  0.16 -0.11  0.04 -0.04  0.12 -0.09  0.00
[14,]  0.47  0.26  0.11  0.13 -0.01 -0.01  0.00  0.07 -0.03  0.02
On peut aussi visualiser graphiquement les différentes valeurs, avec en ordonnée les ordres autoréregressifs (http://freakonometrics.blog.free.fr/public/perso6/platex.gif) et en abscisse les ordres moyenne mobile (http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif). 
> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> ceacf=matrix(as.numeric(cut(EACF$eacf,
+ ((-3):3)/3,labels=1:6)),nrow(EACF$eacf),
+ ncol(EACF$eacf))
> for(i in 1:ncol(EACF$eacf)){
+ for(j in 1:nrow(EACF$eacf)){
+ polygon(c(i-1,i-1,i,i)-.5,c(j-1,j,j,j-1)-.5,
+ col=CL[ceacf[j,i]])
+ }}

Un bruit blanc est en bas à gauche (http://freakonometrics.blog.free.fr/public/perso6/platex.gif et http://freakonometrics.blog.free.fr/public/perso6/qlatex.gif nuls). Dans cette méthode, dite méthode des coins, on cherche un coin telle que dans le quadrant supérieur, les valeurs soient non-significatives (ternes sur le dessin ci-dessus). La figure ci-desssous correspond à l'analyse d'un modèle http://freakonometrics.blog.free.fr/public/perso6/ARMA126latex.gif

Les fortes valeurs positives sont en bleu foncé, les fortes valeurs négatives sont en rouge foncé. On peut aussi regarder la fonction suivante, qui utilise une identification de modèle par MINIC (Minimum Information Criterion) à l'aide du critère de Shcwarz (BIC, ou SBC Schwarz's Bayesian Criterion)
> armaselect(A7,nbmod=5)
p q      sbc
[1,] 12 1 1441.798
[2,] 13 0 1442.628
[3,] 12 0 1443.188
[4,] 12 2 1443.362
[5,] 14 0 1445.069
Enfin, une dernière fonction possible est évoquée dans la section 6.5. du livre de Cryer & Chan (2008),
> ARMA.SELECTION=
+ armasubsets(A7,nar=14,nma=14,ar.method='ols')
> plot(ARMA.SELECTION)

basé sur le critère de Schwarz.
Je vais me répéter, mais ces méthodes ne sont que des outils, histoire d'avoir des pistes si on ne sait trop dans quelle direction partir. Et compte tenu de la saisonnalité de la série, je serais pour ma part parti sur un modèle avec http://freakonometrics.blog.free.fr/public/perso6/platex.gif=12, histoire de voir si on ne pourrait pas avoir un modèle simple, et facilement interprétable en plus...

Monday, September 10 2012

Unit root, or not ? is it a big deal ?

Consider a time series, generated using

set.seed(1)
E=rnorm(240)
X=rep(NA,240)
rho=0.8
X[1]=0
for(t in 2:240){X[t]=rho*X[t-1]+E[t]}

The idea is to assume that an autoregressive model can be considered, but we don't know the value of the parameter. More precisely, we can't choose if the parameter is either one (and the series is integrated), or some value strictly smaller than 1 (and the series is stationary). Based on past observations, the higher the autocorrelation, the lower the variance of the noise.

rhoest=0.9; H=260
u=241:(240+H)
P=X[240]*rhoest^(1:H)
s=sqrt(1/(sum((rhoest^(2*(1:300))))))*sd(X)
Now that we have a model, consider the following forecast, including a confidence interval,
 
plot(1:240,X,xlab="",xlim=c(0,240+H),
ylim=c(-9.25,9),ylab="",type="l")
V=cumsum(rhoest^(2*(1:H)))*s
polygon(c(u,rev(u)),c(P+1.96*sqrt(V),
rev(P-1.96*sqrt(V))),col="yellow",border=NA)
polygon(c(u,rev(u)),c(P+1.64*sqrt(V),
rev(P-1.64*sqrt(V))),col="orange",border=NA)
lines(u,P,col="red")
Here, forecasts can be derived, with any kind of possible autoregressive coefficient, from 0.7 to 1. I.e. we can chose to model the time series either with a stationary, or an integrated series,
.
As we can see above, assuming either that the series is stationary (parameter lower - strictly - than 1) or integrated (parameter equal to 1), the shape of the prediction can be quite different. So yes, assuming an integrated model is a big deal, since it has a strong impact on predictions.

Les tests de non-stationnarité (racine unité)

Pour commencer la modélisation d'une série temporelle, la première étape est de savoir si on peut la considérer comme stationnaire. L'alternative étant qu'elle soit intégrée. Le test le plus classique est le test de Dickey-Fuller.De manière générale, considérons un modèle autorégressive à l'ordre 1,

http://freakonometrics.blog.free.fr/public/perso6/adf-1.gif

On peut éventuellement réécrire ce processus sous la forme

http://freakonometrics.blog.free.fr/public/perso6/acf-18.gif

http://freakonometrics.blog.free.fr/public/perso6/adf-21.gif. On aura un processus non-stationnaire dès lors que http://freakonometrics.blog.free.fr/public/perso6/adf-20.gif ou encore http://freakonometrics.blog.free.fr/public/perso6/adf-16.gif. Si on n'a pas de tendance linéaire, tester http://freakonometrics.blog.free.fr/public/perso6/acf-22.gif ou http://freakonometrics.blog.free.fr/public/perso6/adf-26.gif peut se faire avec un test de Student. L'idée de Dickey-Fuller est de généraliser le test de Student, en rajoutant cette tendance linéaire (pour commencer).

Afin d'illustrer l'utilisation de ce test, commençons par un cas "simple": on va intégrer un bruit blanc simulé.

>  set.seed(1)
>  E=rnorm(240)
>  X=cumsum(E)
>  plot(X,type="l")

On peut regarder s'il y a une constante non nulle (on parlera de modèle avec drift) voire une tendance linéaire (on parlera de trend),

La première étape pourrait être de tenter un modèle sans rien.
> df=ur.df(X,type="none",lags=1)
> summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression none
 
 
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
 
Residuals:
Min       1Q   Median       3Q      Max
-2.87492 -0.53977 -0.00688  0.64481  2.47556
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1    -0.005394   0.007361  -0.733    0.464
z.diff.lag -0.028972   0.065113  -0.445    0.657
 
Residual standard error: 0.9666 on 236 degrees of freedom
Multiple R-squared: 0.003292,	Adjusted R-squared: -0.005155
F-statistic: 0.3898 on 2 and 236 DF,  p-value: 0.6777
 
 
Value of test-statistic is: -0.7328
 
Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62
La région critique du test est ici (pour un seuil à 5%) l'ensemble des valeurs inférieures à -1.95. Or ici la statistique de test est ici -0.73, on est alors dans la région d'acceptation du test, et on va retenir l'hypothèse de racine unité, i.e. la série n'est pas stationnaire. Mais peut-être faudrait-il juste prendre en compte une constante ?
> df=ur.df(X,type="drift",lags=1)
> summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression drift
 
 
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
 
Residuals:
Min       1Q   Median       3Q      Max
-2.91930 -0.56731 -0.00548  0.62932  2.45178
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.29175    0.13153   2.218   0.0275 *
z.lag.1     -0.03559    0.01545  -2.304   0.0221 *
z.diff.lag  -0.01976    0.06471  -0.305   0.7603
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.9586 on 235 degrees of freedom
Multiple R-squared: 0.02313,	Adjusted R-squared: 0.01482
F-statistic: 2.782 on 2 and 235 DF,  p-value: 0.06393
 
 
Value of test-statistic is: -2.3039 2.7329
 
Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81
Deux statistiques de test sont calculées, ici: la première relative à la racine unité, la seconde relative à la constante.  On observe ici que la statistique de test pour le test de racine unité (-2.3) est ici supérieure à toutes les valeurs critiques associées (données sur la première ligne). On va encore accepter l'hypothèse de racine unité. Mais le modèle était peut-être faux, et peut-être avait-on en fait en tendance linéaire ?
> df=ur.df(X,type="trend",lags=1)
> summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression trend
 
 
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
 
Residuals:
Min       1Q   Median       3Q      Max
-2.87727 -0.58802 -0.00175  0.60359  2.47789
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3227245  0.1502083   2.149   0.0327 *
z.lag.1     -0.0329780  0.0166319  -1.983   0.0486 *
tt          -0.0004194  0.0009767  -0.429   0.6680
z.diff.lag  -0.0230547  0.0652767  -0.353   0.7243
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.9603 on 234 degrees of freedom
Multiple R-squared: 0.0239,	Adjusted R-squared: 0.01139
F-statistic:  1.91 on 3 and 234 DF,  p-value: 0.1287
 
 
Value of test-statistic is: -1.9828 1.8771 2.7371
 
Critical values for test statistics:
1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47
On obtient cette fois trois statistique, la première relative au test de racine unité, et les deux suivantes aux tests sur la constante, et sur la tendance (la pente de l'ajustement linéaire). Là encore, la valeur de test (-1.98) excède les valeurs critiques: la p-value serait ici supérieure à 10%. On va là encore accepter l'hypothèse de racine unité.
Mais la modélisation autorégressive à l'ordre 1 était peut-être elle aussi fausse. Aussi, il existe un test de Dickey-Fuller augmenté. L'idée est de considérer, de manière beaucoup plus générale, un processus autorégressif à un ordre supérieur,
http://freakonometrics.blog.free.fr/public/perso6/adf-2.gif
Comme auparavant, on peut réécrire ce modèle en fonction des variations (on parlera de modèle à correction d'erreur)
http://freakonometrics.blog.free.fr/public/perso6/acf-29.gif
où on construit les nouveaux coefficients par une relation de récurrence, du genre http://freakonometrics.blog.free.fr/public/perso6/adf-6.gif. Bref, à nouveau, on est ramené à tester http://freakonometrics.blog.free.fr/public/perso6/acf-22.gif ou http://freakonometrics.blog.free.fr/public/perso6/adf-26.gif . Et là encore, on a le choix sur la tendance.
Si on commence par supposer qu'il n'y a pas de tendance,
> df=ur.df(X,type="none",lags=2)
> summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression none
 
 
Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
 
Residuals:
Min       1Q   Median       3Q      Max
-2.86738 -0.53887 -0.02009  0.67058  2.45970
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
z.lag.1     -0.005168   0.007389  -0.699    0.485
z.diff.lag1 -0.029619   0.065289  -0.454    0.650
z.diff.lag2 -0.037856   0.065436  -0.579    0.563
 
Residual standard error: 0.9685 on 234 degrees of freedom
Multiple R-squared: 0.004704,	Adjusted R-squared: -0.008056
F-statistic: 0.3687 on 3 and 234 DF,  p-value: 0.7757
 
 
Value of test-statistic is: -0.6994
 
Critical values for test statistics:
1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62
la valeur de la statistique de test excède là encore les valeurs critiques, i.e. la p-value dépasse 10%. Si on rajoute une constante,
> df=ur.df(X,type="drift",lags=2)
> summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression drift
 
 
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
 
Residuals:
Min       1Q   Median       3Q      Max
-2.91609 -0.56702  0.01025  0.62109  2.43970
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.31105    0.13332   2.333   0.0205 *
z.lag.1     -0.03744    0.01565  -2.392   0.0175 *
z.diff.lag1 -0.01917    0.06483  -0.296   0.7677
z.diff.lag2 -0.02794    0.06496  -0.430   0.6676
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.9594 on 233 degrees of freedom
Multiple R-squared: 0.02664,	Adjusted R-squared: 0.0141
F-statistic: 2.125 on 3 and 233 DF,  p-value: 0.09774
 
 
Value of test-statistic is: -2.3924 2.9709
 
Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81
on valide encore l’hypothèse de racine unité et de même avec une tendance linéaire,
> df=ur.df(X,type="trend",lags=2)
> summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression trend
 
 
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
 
Residuals:
Min       1Q   Median       3Q      Max
-2.85835 -0.58826 -0.00867  0.61407  2.47280
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.3530591  0.1522373   2.319   0.0213 *
z.lag.1     -0.0338831  0.0168523  -2.011   0.0455 *
tt          -0.0005674  0.0009879  -0.574   0.5663
z.diff.lag1 -0.0237595  0.0654158  -0.363   0.7168
z.diff.lag2 -0.0328215  0.0656102  -0.500   0.6174
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.9608 on 232 degrees of freedom
Multiple R-squared: 0.02802,	Adjusted R-squared: 0.01126
F-statistic: 1.672 on 4 and 232 DF,  p-value: 0.1573
 
 
Value of test-statistic is: -2.0106 2.0849 3.0185
 
Critical values for test statistics:
1pct  5pct 10pct
tau3 -3.99 -3.43 -3.13
phi2  6.22  4.75  4.07
phi3  8.43  6.49  5.47
On notera toutefois que ces trois modèles nous suggèrent de ne pas retenir autant de retard, qui ne semblent pas significatifs.
On peut bien entendu faire ces tests sur de vraies données par exemple le niveau du Nil
> library(datasets)
> NILE=Nile

ou encore les taux d'intérêt américains,
> base=read.table(
+ "http://freakonometrics.free.fr/basedata.txt",
+ header=TRUE)
> Y=base[,"R"]
> Y=Y[(base$yr>=1960)&(base$yr<=1996.25)]
> TAUX=ts(Y,frequency = 4, start = c(1960, 1))

Par exemple, sur cette dernière, on rejette l'hypothèse de stationnarité

>  df=ur.df(y=TAUX,lags=3,type="drift")
>  summary(df)
 
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 
 
Test regression drift
 
 
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
 
Residuals:
Min      1Q  Median      3Q     Max
-3.1982 -0.2947 -0.0629  0.3524  3.1899
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.40609    0.16706   2.431 0.016362 *
z.lag.1     -0.06339    0.02500  -2.535 0.012354 *
z.diff.lag1  0.32613    0.08145   4.004 0.000102 ***
z.diff.lag2 -0.31102    0.08027  -3.875 0.000165 ***
z.diff.lag3  0.28712    0.08103   3.543 0.000541 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.7805 on 137 degrees of freedom
Multiple R-squared: 0.2051,	Adjusted R-squared: 0.1819
F-statistic: 8.839 on 4 and 137 DF,  p-value: 2.227e-06
 
 
Value of test-statistic is: -2.5354 3.2457
 
Critical values for test statistics:
1pct  5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1  6.52  4.63  3.81

Mais toutes sortes d'autres tests (plus robustes) peuvent être faits. Les plus connus sont le test de Philips-Perron et le test dit KPSS. Pour ce dernier, il faut spécifier s'il l'on suppose que la série est de moyenne constante, ou si une tendance doit être prise en compte. Si l'on suppose une constante non-nulle, mais pas de tendance, on utilise

> summary(ur.kpss(X,type="mu"))
 
####################### 
# KPSS Unit Root Test # 
####################### 
 
Test is of type: mu with 4 lags.
 
Value of test-statistic is: 0.972
 
Critical value for a significance level of:
10pct  5pct 2.5pct  1pct
critical values 0.347 0.463  0.574 0.739

alors que pour prendre en compte une tendance

> summary(ur.kpss(X,type="tau"))
 
####################### 
# KPSS Unit Root Test # 
####################### 
 
Test is of type: tau with 4 lags.
 
Value of test-statistic is: 0.5057
 
Critical value for a significance level of:
10pct  5pct 2.5pct  1pct
critical values 0.119 0.146  0.176 0.216

Derrière se cache un teste du multiplicateur de Lagrange. L'hypothèse nulle correspond à l'absence de racine unité: plus la statistique de test est grande, plus on s'éloigne de la stationnarité (hypothèse nulle). Avec ces deux tests, on rejette là encore l'hypothèse de stationnarité de notre marche aléatoire simulée.

Pour le test de Philips-Perron, on a un test de type Dickey-Fuller,

> PP.test(X)
 
Phillips-Perron Unit Root Test
 
data:  X
Dickey-Fuller = -2.0116, Trunc lag parameter = 4, p-value = 0.571
Là encore, la p-value nous recommande de valider l'hypothèse de non-stationnarité. A titre de comparaison, si on avait travaillé sur la série différenciée, on aurait accepté l'hypothèse de stationnarité
> PP.test(diff(X))
 
Phillips-Perron Unit Root Test
 
data:  diff(X)
Dickey-Fuller = -15.9522, Trunc lag parameter = 4, p-value = 0.01

Friday, September 7 2012

De l'hebdomadaire au mensuel

Pour le devoir de série temporelles, les données fournies par Google Insight sont hebdomadaires, ce qui peut rendre la modélisation compliquée. Comme on l'a évoqué en fin de cours, il peut être plus simple de travailler sur des données mensuelles. La petite fonction suivante permet de transformer les données pour avoir des données mensuelles (avec des moyennes par mois, ce qui fait qu'un mois de 28 jours et un mois de 31 jours sont comparables). Histoire d'illustrer, prenons un mot au hasard, disons.... épilation. On commence par sauver le fichier csv, et on le lit sous R,

EPILATION=read.table(
"/home/charpentier/report-epilation.csv",
skip=4,header=TRUE,sep=",",nrows=426)

La petite fonction suivant va nous aider à convertir la base en une série temporelle mensuelle,

H2M=function(BASE){
X=BASE[,2]
date1=substr(as.character(BASE$Semaine),1,10)
date2=substr(as.character(BASE$Semaine),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=ts(Z,start=c(2004,1),frequency=12)
return(Zts)}

Si on travaille sur la série hebdomadaire, on a la série suivante

hebdo=ts(EPILATION$épilation,start=2004,
frequency=52)

La fonction d'autocorrélation est alors

acf(hebdo,160)

Maintenant, on peut travailler sur les données mensuelles. On utilise alors

mensuel=H2M(EPILATION)

Le graphique est alors le suivant

La fonction d'autocorrélation est cette fois la suivante

acf(mensuel,40)

On retrouve le comportement cyclique, avec une saisonnalité annuelle. Mais avec 12 retards, on devrait avoir des modèles plus simples qu'avec 52 retards. Bref, il peut être plus simple de travailler sur des données mensuelles qu'hebdomadaires. Mais chacun est libre de la série qu'il ou elle analysera...

That damn R-squared !

Another post about the R-squared coefficient, and about why, after some years teaching econometrics, I still hate when students ask questions about it. Usually, it starts with "I have a _____ R-squared... isn't it too low ?" Please, feel free to fill in the blanks with your favorite (low) number. Say 0.2. To make it simple, there are different answers to that question:

  1. if you don't want to waste time understanding econometrics, I would say something like "Forget about the R-squared, it is useless" (perhaps also "please, think twice about taking that econometrics course")
  2. if you're ready to spend some time to get a better understanding on subtle concepts, I would say "I don't like the R-squared. I might be interesting in some rare cases (you can probably count them on the fingers of one finger), like comparing two models on the same dataset (even so, I would recommend the adjusted one). But usually, its values has no meaning. You can compare 0.2 and 0.3 (and prefer the 0.3 R-squared model, rather than the 0.2 R-squared one), but 0.2 means nothing". Well, not exactly, since it means something, but it is not a measure tjat tells you if you deal with a good or a bad model. Well, again, not exactly, but it is rather difficult to say where bad ends, and where good starts. Actually, it is exactly like the correlation coefficient (well, there is nothing mysterious here since the R-squared can be related to some correlation coefficient, as mentioned in class)
  3. if you want some more advanced advice, I would say "It's complicated..." (and perhaps also "Look in a textbook write by someone more clever than me, you can find hundreds of them in the library !")
  4. if you want me to act like people we've seen recently on TV (during electoral debate), "It's extremely interesting, but before answering your question, let me tell you a story..."
Perhaps that last strategy is the best one, and I should focus on the story. I mean, this is exactly why I have my blog: to tell (nice) stories. With graphs, and math formulas inside. First of all, consider a regression model

so that the R-squared is defined as

Let us generate datasets, and then run regressions, to see what's going on...
For instance, consider 20 observations, with one variable of interest, one explanatory variable, and some low variance noise (to start with)
> set.seed(1)
> n=20
> X=runif(n)
> E=rnorm(n)
> Y=2+5*X+E*.5
> base=data.frame(X,Y)
> reg=lm(Y~X,data=base)
> summary(reg)
 
Call:
lm(formula = Y ~ X, data = base)
 
Residuals:
Min       1Q   Median       3Q      Max
-1.15961 -0.17470  0.08719  0.29409  0.52719
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)   2.4706     0.2297   10.76 2.87e-09 ***
X             4.2042     0.3697   11.37 1.19e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.461 on 18 degrees of freedom
Multiple R-squared: 0.8778,	Adjusted R-squared: 0.871
F-statistic: 129.3 on 1 and 18 DF,  p-value: 1.192e-09 
The R-squared is high (almost 0.9). What if the underlying model is exactly the same, but now, the noise has a much higher variance ?
> Y=2+5*X+E*4
> base=data.frame(X,Y)
> reg=lm(Y~X,data=base)
> summary(reg)
 
Call:
lm(formula = Y ~ X, data = base)
 
Residuals:
Min      1Q  Median      3Q     Max
-9.2769 -1.3976  0.6976  2.3527  4.2175
 
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)    5.765      1.837   3.138  0.00569 **
X             -1.367      2.957  -0.462  0.64953
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 3.688 on 18 degrees of freedom
Multiple R-squared: 0.01173,	Adjusted R-squared: -0.04318
F-statistic: 0.2136 on 1 and 18 DF,  p-value: 0.6495 
Now, the R-squared is rather low (around 0.01). Thus, the quality of the regression depends clearly on the variance of the noise. The higher the variance, the lower the R-squared. And usually, there is not much you can do about it ! On the graph below, the noise is changing, from no-noise, to extremely noisy, with the least square regression in blue (and a confidence interval on the prediction)
If we compare with the graph below, one can observe that the quality of the fit depends on the sample size, with now 100 observations (instead of 20),
So far, nothing new (if you remember what was said in class). In those two cases, here is the evolution of the R-squared, as a function of the variance of the noise (more precisely, here, the standard deviation of the noise) 
> S=seq(0,4,by=.2)
> R2=rep(NA,length(S))
> for(s in 1:length(S)){
+ Y=2+5*X+E*S[s]
+ base=data.frame(X,Y)
+ reg=lm(Y~X,data=base)
+ R2[s]=summary(reg)$r.squared}
with 20 obervations in blue, 100 in black. One important point is that in econometrics, we rarely choose the number of observations. If we have only 100 observations, we have to deal with it. Similarly, if observations are quite noisy there is (usually) not much we can do about it. All the more if you don't have any explanatory variable left. Perhaps you migh play (or try to play) with nonlinear effect...

Nevertheless, it looks like some econometricians really care about the R-squared, and cannot imagine looking at a model if the R-squared is lower than - say - 0.4. It is always possible to reach that level ! you just have to add more covariates ! If you have some... And if you don't, it is always possible to use polynomials of a continuous variate. For instance, on the previous example,

> S=seq(1,25,by=1)
> R2=rep(NA,length(S))
> for(s in 1:length(S)){
+ reg=lm(Y~poly(X,degree=s),data=base)
+ R2[s]=summary(reg)$r.squared}
If we plot the R-squared as a function of the degree of the polynomial regression, we have the following graph. Once again, the higher the degree, the more covariates, and the more covariates, the higher the R-squared,
I.e. with 22 degrees, it is possible to reach a 0.4 R-squared. But it might be interesting to the prediction we have with that model,
So, was it worth adding so much polynomial parts ? I mean, 22 is quite a large power... Here, the linear regression was significant, but not great. So what ? The R-squared was small ? so what ? sometimes, there's not much you can do about it... When dealing with individual observations (so called micro-econometrics), the variable of interest might be extremely noisy, and there is not much you can do. So your R-squared can be small, but your regression is perhaps still better than doing nothing... The good thing with a low R-squared is perhaps that it will recall us that we have to remain modest when we build a model. And always be cautious about conclusion. At least, provide a confidence interval...

- page 1 of 3