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,
On peut éventuellement réécrire ce processus sous la forme
où . 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).
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,
Comme auparavant, on peut réécrire ce modèle en fonction des variations (on parlera de modèle à correction d’erreur)
où on construit les nouveaux coefficients par une relation de récurrence, du genre . Bref, à nouveau, on est ramené à tester ou . 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é