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é

> PP.test(diff(X))

Phillips-Perron Unit Root Test

data:  diff(X)
Dickey-Fuller = -15.9522, Trunc lag parameter = 4, p-value = 0.01