Quelques transparents en plus, qui devraient correspondre aux deux prochains cours de séries temporelles, sur les processus autorégressifs (AR) et moyennes mobiles (MA), les ARMA, les ARIMA (intégrés) et les SARIMA (saisonniers). J’ai mis des notes sur les tests de racine unités, je rajouterais quelques transparents la semaine prochaine sur les tests de saisonnalité, et quelques exemples pratiques de prévision. Les transparents sont en ligne ici,
Tag Archives: AR
Inference and autoregressive processes
Consider a (stationary) autoregressive process, say of order 2,
for some white noise with variance . Here is a code to generate such a process,
> phi1=.5 > phi2=-.4 > sigma=1.5 > set.seed(1) > n=240 > WN=rnorm(n,sd=sigma) > Z=rep(NA,n) > Z[1:2]=rnorm(2,0,1) > for(t in 3:n){Z[t]=phi1*Z[t-1]+phi2*Z[t-2]+WN[t]}
Here, we have to estimate two sets of parameters: the autoregressive coefficients, and the variance of the innovation process . There are (at least) three techniques to estimate those parameters.
- using least square regression
A natural idea is to see here a regression model, and thus, if we consider a matrix formulation,
Here we can run (conditional) ordinary least squares estimation,
> base=data.frame(Y=Z[3:n],X1=Z[2:(n-1)],X2=Z[1:(n-2)]) > regression=lm(Y~0+X1+X2,data=base) > summary(regression) Call: lm(formula = Y ~ 0 + X1 + X2, data = base) Residuals: Min 1Q Median 3Q Max -4.3491 -0.8890 -0.0762 0.9601 3.6105 Coefficients: Estimate Std. Error t value Pr(>|t|) X1 0.45107 0.05924 7.615 6.34e-13 *** X2 -0.41454 0.05924 -6.998 2.67e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.449 on 236 degrees of freedom Multiple R-squared: 0.2561, Adjusted R-squared: 0.2497 F-statistic: 40.61 on 2 and 236 DF, p-value: 6.949e-16 > regression$coefficients X1 X2 0.4510703 -0.4145365 > summary(regression)$sigma [1] 1.449276
- using Yule-Walker equations
As we’ve seen in class, we can easily get the following equations for the autocovariance functions,
which can also be written
So we just have to solve a simple linear system of equations. Note that if we divide by the variance, those equations can be written in terms of the autocorrelation functions
The code is the following
> rho1=cor(Z[1:(n-1)],Z[2:n]) > rho2=cor(Z[1:(n-2)],Z[3:n]) > A=matrix(c(1,rho1,rho1,1),2,2) > b=matrix(c(rho1,rho2),2,1) > (PHI=solve(A,b)) [,1] [1,] 0.4517579 [2,] -0.4155920
Now, we need to extract the estimated innovation process, from this set of parameters (note that it could be possible to include the variance term in Yule-Walker equations, to get a three dimensional linear equation)
> estWN=base$Y-(PHI[1]*base$X1+PHI[2]*base$X2) > sd(estWN) [1] 1.445706
This estimator is probably not the best one (we can take into account that we’ve lost two degrees of freedom), but as a starting point, let us consider this one.
- using (conditional) likelihood estimators
Finally, we can assume some distribution for the innovation process. Thestandard model is a Gaussian model, i.e.
In that case, the conditional log likelihood (conditional since we set the first two observations here) is
> CondLogLik=function(A,TS){ + phi1=A[1]; phi2=A[2] + sigma=A[3] ; L=0 + for(t in 3:length(TS)){ + L=L+dnorm(TS[t],mean=phi1*TS[t-1]+ + phi2*TS[t-2],sd=sigma,log=TRUE)} + return(-L)}
Now, we can run standard optimization procedures,
> LogL=function(A) CondLogLik(A,TS=Z) > optim(c(0,0,1),LogL) $par [1] 0.4509685 -0.4144938 1.4430930 $value [1] 425.0164 $counts function gradient 88 NA $convergence [1] 0 $message NULL
Here, our three estimators are rather close. Actually, if we generate 1,000 time series (of size 240), those are the Box-plots of our three estimators, for the first order autoregressive coefficient
for the second one,
and finally for the standard deviation of the innovation process
All those estimators behave nicely, and are rather close. Note that they all might be biased, but they are consistent (see Davidson and MacKinnon for instance, in their book, for more details).
La maudite constante des modèles ARIMA
Dans les modèles ARIMA, autant le dire tout de suite, les constantes c’est pénible ! Pour comprendre un peu mieux ce qui se passe, considérons ici un modèle ARMA avec constante. Ou pour commencer, juste un AR(1)
Supposons que la série soit stationnaire, de moyenne . Alors en prenant l’espérance de part et d’autre dans l’expression précédente, , i.e. , ou . Plus généralement, avec un ARMA, , en prenant la encore l’espérance des deux cotés, on en déduit .
Simulons un processus AR(1),
Pour simuler des processus ARIMA, on pourrait utiliser la commande
> X=arima.sim(list(order=c(1,0,0),ar=1/3),n=1000)+2 > mean(X) [1] 1.931767
mais comme le but est de comprendre ce qui se passe, autant faire les choses calmement,
> X=rep(NA,1010) > X[1]=0 > for(t in 2:1010){X[t]=4/3+X[t-1]/3+rnorm(1)} > X=X[-(1:10)] > mean(X) [1] 2.03397
I.e. ici le processus est de moyenne qui vaut 2.
Regardons maintenant ce que donnerait l’estimation du processus AR(1),
> arima(X, order = c(1, 0, 0)) Call: arima(x = X, order = c(1, 0, 0)) Coefficients: ar1 intercept 0.3738 2.0334 s.e. 0.0294 0.0487 sigma^2 estimated as 0.9318: log likelihood = -1383.68
De manière un peu surprenante, le coefficient appelé intercept n’est pas la constante dans le modèle AR(1), mais la moyenne du processus. Autrement dit, R n’ajuste pas un processus
comme le laisserait penser l’intuition, mais un processus
Ces deux formes sont bien entendu équivalentes. Mais les coefficients estimés ne sont pas tout à fait ce que l’on attendait…
Plaçons nous maintenant dans le cas d’un processus non-stationnaire. L’extension naturelle serait de considérer un processus ARIMA(1,1,0), ou bien un processus tel que la différence soit un processus AR(1). Un processus ARIMA(1,1,0) avec une constante s’écrirait
en utilisant l’opérateur retard. Ceci fait penser à un modèle avec une tendance linéaire. Posons (afin de se débarrasser de cette tendance). Alors
i.e.
ou encore
On note que ou encore suit alors un processus ARIMA(1,1,0) sans constante cette fois.
Afin de visualiser ce que donnerait l’inférence, simulons le processus suivant, qui est un processus AR(1) avec constante que l’on intègre: avec
Commençons (là encore) par simuler un tel processus,
> U=rep(NA,1010) > U[1]=0 > for(t in 2:1010){U[t]=4/3+U[t-1]/3+rnorm(1)} > U=U[-(1:10)] > X=cumsum(U)
Sous R, on obtient l’estimation suivant si l’on tente de calibrer un modèle ARIMA(1,1,0)
> arima(X, order = c(1, 1, 0)) Call: arima(x = X, order = c(1, 1, 0)) Coefficients: ar1 0.8616 s.e. 0.0160
sigma^2 estimated as1.343:loglikelihood = -1565.63Ça ne convient pas du tout…. On peut tenter un processus AR(1) (avec constante) sur la série différenciée…
> arima(diff(X), order = c(1, 0, 0)) Call: arima(x = diff(X), order = c(1, 0, 0)) Coefficients: ar1 intercept 0.3564 2.0200 s.e. 0.0295 0.0486 sigma^2 estimated as 0.9782: log likelihood = -1406.6
On progresse, sauf que comme auparavant, le terme qui est donne n’est pas la constante dans le modèle ARIMA, mais la moyenne du processus différencié. Mais cette fois, on a un interprétation, c’est que la constante est la pente de la tendance ! Si on estime la pente associée a , on récupère la même valeur:
> arima(X, order = c(1, 1, 0), xreg=1:length(X)) Call: arima(x = X, order = c(1, 1, 0), xreg = 1:length(X)) Coefficients: ar1 1:length(X) 0.3566 2.0519 s.e. 0.0296 0.0487 sigma^2 estimated as 0.9787: log likelihood = -1406.82
Sur la figure ci-dessous, on retrouve le fait qu’en enlevant la tendance linéaire à donneune série intégrée, sans constante (qui fait penser à une marche aléatoire).
Autrement dit
- dans le cas d’une série stationnaire, la constante estimée n’est pas du tout la constante, mais la moyenne de la série temporelle
- dans le cas d’une série non-stationnaire, la constante estimée dans la série différenciée a du sens, au sens ou il s’agit de la pente de la tendance (linéaire) du processus
Avant de conclure, une petite remarque. Quid de la prévision ? Si on commencer par une révision à l’aide du premier processus ARIMA, on obtient une prédiction pour avec un énorme intervalle de confiance (sans aucun bon sens)
> ARIMA1=arima(X, order = c(1, 1, 0)) > ARIMA2=arima(X, order = c(1, 1, 0), xreg=1:length(X)) > Xp1=predict(ARIMA1,20) > Xp2=predict(ARIMA2,20,newxreg= + (length(X)+1):(length(X)+20)) > plot(960:1000,X[960:1000],xlim=c(960,1020),type="l") > polygon(c(1001:1020,rev(1001:1020)), + c(Xp1$pred+2*Xp1$se,rev(Xp1$pred-2*Xp1$se)), + col=CL[3],border=NA) > lines(1001:1020,Xp1$pred,col="red",lwd=2)
intervalle qui se visualise sur le graphique ci-dessous,
Tous ceux qui auront reconnu ici des problèmes qui se posent dans la modélisation de la série dans le modèle de Lee & Carter (1992) auront compris qu’il existe une vraie application a tout ce que je viens de raconter (je pense au particulier au commentaire poste par @ClaudeT en début de semaine) car la série des ressemble très fortement à ce genre de série, non-stationnaire avec une tendance linéaire. Et qu’un modèle mal spécifié laisse à penser que l’incertitude est beaucoup beaucoup plus grande que ce qu’elle est vraiment (en plus de donner une prédiction surprenante à long terme !). Donc avant de faire tourner rapidement les codes, il vaut mieux regarder attentivement ce qu’ils font réellement…