{"word_count": 825, "character_count": 9244}
Après discussion, la date limite pour me renvoyer le second devoir est fixée à vendredi 5 octobre midi, par courriel. Je veux un fichier pdf par équipe, avec une page de garde indiquant le nom des membres du groupe, et le mot clé retenu.Maintenant, je vais en profiter pour mettre en ligne le code tapé de matin, en cours. On avait travaillé sur la modélisation de la série des ventesde voiture, au Québec. La série est en ligne,X=read.table("http://freakonometrics.blog.free.fr/public/data/car-sales-quebec.csv",header=TRUE,sep=";",nrows=108)Xt=ts(X[,2],start=c(1960,1),frequency=12)On peut regarder l'évolution de cette série temporelles (car les dessins, et la visualisation, c'est important),plot(Xt)On note une tendance linéaire (ou qu'on pourrait supposer linéaire), que l'on va estimer,X=as.numeric(Xt)temps=1:length(X)base=data.frame(temps,X)reg=lm(X~temps)et que l'on pourra représenterplot(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 iciY=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 .&#62; fit.ar12=arima(Y,order=c(12,0,0))&#62; fit.ar12Series: YARIMA(12,0,0) with non-zero meanCoefficients:ar1     ar2      ar3      ar4     ar5      ar6     ar7      0.1975  0.0832  -0.1062  -0.1212  0.1437  -0.1051  0.0319s.e.  0.0838  0.0809   0.0826   0.0843  0.0850   0.0833  0.0854ar9     ar10    ar11    ar12  intercept      -0.0332  -0.0616  0.2635  0.4913  -148.3180s.e.   0.0853   0.0840  0.0840  0.0841   384.5095sigma^2 estimated as 2177974:  log likelihood=-946.75AIC=1921.51   AICc=1926.03   BIC=1959.06La 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 paspoints(-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élationsacf(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&#62; Box.test(residuals(fit.ar12),lag=12,+  type='Box-Pierce')Box-Pierce testdata:  residuals(fit.ar12)X-squared = 7.7883, df = 12, p-value = 0.8014Pour 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,&#62; Box.test(residuals(fit.ar12),lag=18,+  type='Box-Pierce')Box-Pierce testdata:  residuals(fit.ar12)X-squared = 20.3861, df = 18, p-value = 0.3115on a la même conclusion,Pour rappels, on peut retrouver à la main la p-value,&#62; 1-pchisq(20.3861,df=18)[1] 0.3115071Faisons 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.valueplot(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,&#62; library(caschrono)&#62; 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.232On peut être tenté d'aller voir ce qui se passerait en rajoutant une composante moyenne mobile à notre modèle précédant,&#62; fit.arma12.1=arima(Y,order=c(12,0,1))&#62; fit.arma12.1Series: YARIMA(12,0,1) with non-zero meanCoefficients:ar1     ar2      ar3      ar4     ar5      ar6     ar7      0.0301  0.1558  -0.0941  -0.1461  0.1063  -0.0688  -0.002s.e.  0.1235  0.0854   0.0757   0.0784  0.0807   0.0774   0.080ar9     ar10    ar11    ar12     ma1  intercept      -0.0646  -0.0798  0.2538  0.5786  0.2231  -131.3495s.e.   0.0802   0.0766  0.0751  0.0861  0.1393   368.8156sigma^2 estimated as 2127759:  log likelihood=-945.65AIC=1921.31   AICc=1926.52   BIC=1961.54Le 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.valueplot(1:24,Vectorize(BP)(1:24),type='b',ylim=c(0,1),col="red")abline(h=.05,lty=2,col="blue")BP=function(h) Box.test(residuals(fit.ar12),lag=h,type='Ljung-Box')$p.valuelines(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 .&#62; fit.ar14=+ arima(Y,order=c(14,0,0),method="CSS")&#62; fit.ar14Series: YARIMA(14,0,0) with non-zero meanCoefficients:ar1     ar2      ar3      ar4     ar5      ar6     ar7      0.2495  0.2105  -0.0584  -0.1569  0.1282  -0.1152  0.0268s.e.  0.0956  0.0972   0.0854   0.0830  0.0838   0.0840  0.0847ar9     ar10    ar11    ar12     ar13     ar14  intercept      -0.0327  -0.1116  0.2649  0.5887  -0.1575  -0.1572    80.5s.e.   0.0855   0.0851  0.0853  0.0886   0.1031   0.0999   338.9sigma^2 estimated as 2218612:  part log likelihood=-942.31C'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=12subY=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)&#62; (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èle397  -4836.2174 -5689.3331 -5885.4486 -6364.247198  -3876.4199 -4274.0391 -4287.2193 -4773.811699   1930.3776  1817.8411  2127.9915  2290.1460100  3435.1751  4089.3598  3736.1110  4039.4150101  7727.9726  6998.9829  7391.6694  7281.4797102  2631.7701  3456.8819  3397.5478  4230.5324103  -509.4324 -2128.6315 -2268.9672 -2258.7216104 -1892.6349 -3877.7609 -3694.9409 -3620.4798105 -4310.8374 -3384.0905 -3430.4090 -2881.4942106  2564.9600  -504.6883  -242.5018   183.2891107 -1678.2425 -1540.9904 -1607.5996  -855.7677108 -4362.4450 -3927.4772 -3928.0626 -3718.3922&#62;  sum((M[,1]-M[,2])^2)[1] 19590931&#62;  sum((M[,1]-M[,3])^2)[1] 17293716&#62;  sum((M[,1]-M[,4])^2)[1] 21242230I.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.1PREDARMA=forecast(fit.arma12.1,12)plot(Y,xlim=c(1,120),type="l")temps=T+1:12polygon(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:120polygon(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:120polygon(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")