Après discussion, la date limite pour me renvoyer le second devoir est fixée à vendredi 5 octobre midi, par courriel. Je veux un fichier pdf par équipe, avec une page de garde indiquant le nom des membres du groupe, et le mot clé retenu.
Maintenant, je vais en profiter pour mettre en ligne le code tapé de matin, en cours. On avait travaillé sur la modélisation de la série des ventesde voiture, au Québec. La série est en ligne,
X=read.table( "http://freakonometrics.blog.free.fr/public/ data/car-sales-quebec.csv", header=TRUE,sep=";",nrows=108) Xt=ts(X[,2],start=c(1960,1),frequency=12)
On peut regarder l’évolution de cette série temporelles (car les dessins, et la visualisation, c’est important),
plot(Xt)

On note une tendance linéaire (ou qu’on pourrait supposer linéaire), que l’on va estimer,
X=as.numeric(Xt) temps=1:length(X) base=data.frame(temps,X) reg=lm(X~temps)
et que l’on pourra représenter
plot(temps,X,type="l",xlim=c(0,120)) abline(reg,col="red")
Si on veut une série stationnaire (et c’est effectivement ce que l’on cherche), on va retrancher cette tendance à notre série, et travailler sur la série résiduelle.

La série résiduelle est ici
Y=X-predict(reg) plot(temps,Y,type="l")

que l’on devrait pouvoir supposer stationnaire. Classiquement, on regarde les autocorrélations,
acf(Y,36,lwd=4)

où on repère un beau cycle annuel, mais on va surtout utiliser les autocorrélations partielles pour identifier l’ordre de la composante autorégressive.
pacf(Y,36,lwd=4)

La 12ème est non-nulle, et tous les autres ensuite sont significativement nulles (ou presque). On va donc tenter un .
> fit.ar12=arima(Y,order=c(12,0,0)) > fit.ar12 Series: Y ARIMA(12,0,0) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.1975 0.0832 -0.1062 -0.1212 0.1437 -0.1051 0.0319 s.e. 0.0838 0.0809 0.0826 0.0843 0.0850 0.0833 0.0854 ar9 ar10 ar11 ar12 intercept -0.0332 -0.0616 0.2635 0.4913 -148.3180 s.e. 0.0853 0.0840 0.0840 0.0841 384.5095 sigma^2 estimated as 2177974: log likelihood=-946.75 AIC=1921.51 AICc=1926.03 BIC=1959.06
La douzième composante semble significative, si on repense au test de Student,
u=seq(-6,6,by=.1) plot(u,dnorm(u),type="l") points(0.4913/0.0841,0,pch=19,col="red")
En revanche, la huitième ne l’est pas
points(-0.1018/0.0847,0,pch=19,col="blue")

(voilà pour le petit retour sur la lecture des tests). Vérifions maintenant que le bruit résiduel est bien blanc… Si on regarde les autocorrélations
acf(residuals(fit.ar12),36,lwd=4)

Mais ce n’est pas un test de bruit blanc. En particulier, si une autocorrélation semblerait significative, on pourrait accepter que – globalement – les autocorrélations soient très proches de 0. On va alors faire des tests de Box-Pierce
> Box.test(residuals(fit.ar12),lag=12, + type='Box-Pierce') Box-Pierce test data: residuals(fit.ar12) X-squared = 7.7883, df = 12, p-value = 0.8014
Pour l’interpétation du test, c’est toujours pareil: on a une statistique de l’ordre de 7, et la loi sous-jacente (la somme des carrés des 12 premières autocorrélations du bruit) est une loi du chi-deux, à 12 degrés de liberté,
u=seq(0,30,by=.1) plot(u,dchisq(u,df=12),type="l") points(7.7883,0,pch=19, col="red")

On accepte l’hypothèse que les 12 premières autocorrélations soient nulles. Si on va un peu plus loin,
> Box.test(residuals(fit.ar12),lag=18, + type='Box-Pierce') Box-Pierce test data: residuals(fit.ar12) X-squared = 20.3861, df = 18, p-value = 0.3115
on a la même conclusion,

Pour rappels, on peut retrouver à la main la p-value,
> 1-pchisq(20.3861,df=18) [1] 0.3115071
Faisons notre petite dessin de toutes les p-value, pour s’assurer que le bruit est bien un bruit blanc,
BP=function(h) Box.test(residuals(fit.ar12),lag=h, type='Box-Pierce')$p.value plot(1:24,Vectorize(BP)(1:24),type='b', ylim=c(0,1),col="red") abline(h=.05,lty=2,col="blue")

Cette fois c’est bon, on tient notre premier modèle, un . Si on veut aller plus loin, on peut regarder ce que donnerait de la sélection automatique,
> library(caschrono) > armaselect(Y,nbmod=5) p q sbc [1,] 14 1 1635.214 [2,] 12 1 1635.645 [3,] 15 1 1638.178 [4,] 12 3 1638.297 [5,] 12 4 1639.232
On peut être tenté d’aller voir ce qui se passerait en rajoutant une composante moyenne mobile à notre modèle précédant,
> fit.arma12.1=arima(Y,order=c(12,0,1)) > fit.arma12.1 Series: Y ARIMA(12,0,1) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.0301 0.1558 -0.0941 -0.1461 0.1063 -0.0688 -0.002 s.e. 0.1235 0.0854 0.0757 0.0784 0.0807 0.0774 0.080 ar9 ar10 ar11 ar12 ma1 intercept -0.0646 -0.0798 0.2538 0.5786 0.2231 -131.3495 s.e. 0.0802 0.0766 0.0751 0.0861 0.1393 368.8156 sigma^2 estimated as 2127759: log likelihood=-945.65 AIC=1921.31 AICc=1926.52 BIC=1961.54
Le nouveau coefficient, est à peine significatif. Éventuellement avec un seuil à 10%… Pourquoi pas? Si on regarde les résidus, sans grande surprise, on a toujours un bruit blanc, encore plus blanc qu’auparavant (en violet sur le dessin ci-dessous)
BP=function(h) Box.test( residuals(fit.arma12.1),lag=h, type='Box-Pierce')$p.value plot(1:24,Vectorize(BP)(1:24),type='b', ylim=c(0,1),col="red") abline(h=.05,lty=2,col="blue") BP=function(h) Box.test(residuals(fit.ar12),lag=h, type='Ljung-Box')$p.value lines(1:24,Vectorize(BP)(1:24),col="purple", type="b")

On peut aussi aller voir parmi les modèles proposés, en particulier le modèle .
> fit.ar14= + arima(Y,order=c(14,0,0),method="CSS") > fit.ar14 Series: Y ARIMA(14,0,0) with non-zero mean Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 0.2495 0.2105 -0.0584 -0.1569 0.1282 -0.1152 0.0268 s.e. 0.0956 0.0972 0.0854 0.0830 0.0838 0.0840 0.0847 ar9 ar10 ar11 ar12 ar13 ar14 intercept -0.0327 -0.1116 0.2649 0.5887 -0.1575 -0.1572 80.5 s.e. 0.0855 0.0851 0.0853 0.0886 0.1031 0.0999 338.9 sigma^2 estimated as 2218612: part log likelihood=-942.31
C’est un peu tiré par les cheveux, mais on pourrait accepter l’hypothèse que soit significativement non-nuls. Mais on est encore limite… Allez, on l’accepte aussi dans notre gang de modèle.
On a finalement trois modèles. Si on fait un peu de backtesting, sur les 12 derniers mois,
T=length(Y) backtest=12 subY=Y[-((T-backtest+1):T)] subtemps=1:(T-backtest) plot(temps,Y,type="l") lines(subtemps,subY,lwd=4) fit.ar12.s=arima(subY, order=c(p=12,d=0,q=0),method="CSS") fit.arma12.1.s=arima(subY, order=c(p=12,d=0,q=1),method="CSS") fit.ar14.s=arima(subY, order=c(p=14,d=0,q=0),method="CSS") p.ar12=predict(fit.ar12.s,12) pred.ar12=as.numeric(p.ar12$pred) p.arma12.1=predict(fit.arma12.1.s,12) pred.arma12.1=as.numeric(p.arma12.1$pred) p.ar14=predict(fit.ar14.s,12) pred.ar14=as.numeric(p.ar14$pred)
on obtient les prévisions suivantes (avec les valeurs observées dans la première colonne)
> (M=cbind(observé=Y[((T-backtest+1):T)], + modèle1=pred.ar12, + modèle2=pred.arma12.1, + modèle3=pred.ar14)) observé modèle1 modèle2 modèle3 97 -4836.2174 -5689.3331 -5885.4486 -6364.2471 98 -3876.4199 -4274.0391 -4287.2193 -4773.8116 99 1930.3776 1817.8411 2127.9915 2290.1460 100 3435.1751 4089.3598 3736.1110 4039.4150 101 7727.9726 6998.9829 7391.6694 7281.4797 102 2631.7701 3456.8819 3397.5478 4230.5324 103 -509.4324 -2128.6315 -2268.9672 -2258.7216 104 -1892.6349 -3877.7609 -3694.9409 -3620.4798 105 -4310.8374 -3384.0905 -3430.4090 -2881.4942 106 2564.9600 -504.6883 -242.5018 183.2891 107 -1678.2425 -1540.9904 -1607.5996 -855.7677 108 -4362.4450 -3927.4772 -3928.0626 -3718.3922 > sum((M[,1]-M[,2])^2) [1] 19590931 > sum((M[,1]-M[,3])^2) [1] 17293716 > sum((M[,1]-M[,4])^2) [1] 21242230
I.e. on aurait envie de retenir le second modèle, . On va maintenant l’utiliser pour faire un peu de prévision,
library(forecast) fit.arma12.1= arima(Y,order=c(12,0,1)) fit.arma12.1 PREDARMA=forecast(fit.arma12.1,12) plot(Y,xlim=c(1,120),type="l") temps=T+1:12 polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2], rev(PREDARMA$upper[,2])),col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1], rev(PREDARMA$upper[,1])),col="orange",border=NA) lines(temps,PREDARMA$mean,col="red")

Cela dit, on peut aussi aller beaucoup plus loin dans la prévision,
PREDARMA=forecast(fit.arma12.1,120) plot(Y,xlim=c(1,210),type="l") temps=T+1:120 polygon(c(temps,rev(temps)),c(PREDARMA$lower[,2], rev(PREDARMA$upper[,2])),col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDARMA$lower[,1], rev(PREDARMA$upper[,1])),col="orange",border=NA) lines(temps,PREDARMA$mean,col="red")

Bon, on y est presque, car on a modélisé , la série obtenue en enlevant la tendance linéaire. Pour remonter sur
, on va rajouter la tendance à la prévision faite auparavant,
X=as.numeric(Xt) temps=1:length(X) plot(temps,X,type="l",xlim=c(0,210), ylim=c(5000,30000)) base=data.frame(temps,X) reg=lm(X~temps) abline(reg,col="red") PREDTENDANCE=predict(reg,newdata= data.frame(temps=T+1:120)) temps=T+1:120 polygon(c(temps,rev(temps)),c(PREDTENDANCE+ PREDARMA$lower[,2],rev(PREDTENDANCE+PREDARMA$upper[,2])), col="yellow",border=NA) polygon(c(temps,rev(temps)),c(PREDTENDANCE+ PREDARMA$lower[,1],rev(PREDTENDANCE+PREDARMA$upper[,1])), col="orange",border=NA) lines(temps,PREDTENDANCE+PREDARMA$mean,col="red")






, 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,
et
dans un modèle
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 






. On aura un processus non-stationnaire dès lors que
ou encore
. Si on n’a pas de tendance linéaire, tester
ou
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).



. Bref, à nouveau, on est ramené à tester 



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



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


with variance
. Here is a code to generate such a process,
, and the variance of the innovation process 


, those equations can be written in terms of the autocorrelation functions






soit stationnaire, de moyenne
. Alors en prenant l’espérance de part et d’autre dans l’expression précédente,
, i.e.
, ou
. Plus généralement, avec un ARMA,
, en prenant la encore l’espérance des deux cotés, on en déduit
.



(afin de se débarrasser de cette tendance). Alors


ou encore
suit alors un processus ARIMA(1,1,0) sans constante cette fois.
avec



dans le modèle de
désigne un poids attribué à la nouvelle observation dans la fonction de lissage
. Le code pour lisser une série est alors assez simple,
sur le lissage. La prévision que l’on fait à la date
, pour un horizon
est alors
. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l’erreur de prédiction commise à un horizon de 1




. Le code pour faire un lissage double est là encore assez simple,
(
, où on suppose que le bruit
est ici un bruit indépendant (imprévisible). On va alors supposer que

correspondant à la variation du processus
qui correspond à un cycle (éventuel)
, à la date
(avec les notations du 
Via 

sept mots, le journaliste a réussi a utiliser trois clichés qui m’agacent.

















is the survival copula associated to
, i.e.





Yesterday, I wrote a post (
But as @
Now both stocks are expressed in the same currency. To compare returns volatility, a first idea can be to use GARCH models,




