R Notebook : Analyse des résidus

setwd(dir = "C:/Users/fkraus/Desktop/TD/2025-2026/ECONOMETRIE/séance4")

Les résidus d’une régression sont souvent notés \(\varepsilon_t\) par convention. Ils représentent l’écart entre la variable dépendante de la régression et ses valeurs prédites par le modèles : \[\begin{align} &Y_t = a_0 + \beta X_t + \varepsilon_t \\ &\varepsilon_t = Y_t - \hat Y_t \end{align}\]

Et les résidus représentent réellement à quel point notre modèle est précis dans son estimation. Le modèle de régression est censé représenter correctement la dynamique de la variable dépendante \(Y_t\) en fonction de la dynamique d’autres variables (\(X_t\)), du temps (avec une trend, un cycle, une saisonnalité, etc… si on observe leur présence) et de la structure propre de la série (ARIMA). Cependant, on n’est pas censé représenter la variable dépendante avec une précision parfaite : le modèle peut sur-estimer ou sous-estimer de temps en temps la dynamique.

Cependant, quand on fait le modèle de régression, on part du principe que, si il est correct, on va correctement estimer la dynamique dans l’ensemble. On peut sur-estimer ou sous-estimer, mais en moyenne il faut que ces erreurs soient le plus proche possible de 0 : la moyenne des résidus doit donc être proche de 0. Egalement, on veut se tromper autant dans la sur-estimation que dans la sous-estimation (on ne veut pas sous-estimer la variable en permanence) : la somme des résidus doit être proche de 0. Egalement, on veut se tromper du même ordre peu importe la période : la moyenne et la variance des résidus doivent être constants. De cette condition découle une autre : les résidus doivent être indépendants. Autrement dit, le coefficient d’autocorrélation des résidus doit être nul (si on se trompe à une période, ça ne veut pas dire qu’on va se tromper à la période suivante) et ne doivent pas dépendre de \(Y_t\) (le coefficient de régression entre \(\varepsilon_t\) et \(Y_t\) doit être nul). Au final, les résidus doivent donc suivre une loi normale, centrée autour de zéro.

Si on résume ces cas : \[\begin{align} & \mathbb{E}(\varepsilon) = 0 \\ & \sum_{i=1}^N\varepsilon_i = 0 \\ & \mathbb{E}(\varepsilon)=\mathbb{E}(\varepsilon_{t+h}) \\ & \operatorname{Var}(\varepsilon) = \sigma^2_\varepsilon \\ & \operatorname{Cov}(\varepsilon_t, \varepsilon_{t-h})=0 \\ & \operatorname{Cov}(\varepsilon_t, Y_{t})=0\\ & \rightarrow\varepsilon_t \overset{\text{i.i.d.}}{\sim} \mathcal N(0,\sigma_\varepsilon^2) \end{align}\]

Note : \(i.i.d\) signifie identiquement et indépendamment distribués. Ces cas sont nombreux, mais au final on retombe assez vite sur ce qu’on a vu dans la séance sur la stationnarisation : on veut obtenir une suite de bruits blancs.

En effet, il faut imaginer que les séries temporelles ne sont qu’une suite de bruits blancs, et que ces bruits blancs sont impactés par des chocs exogènes (qui proviennent d’autres variables). Avec la régression, on retire l’influence de ces variables exogènes, et on est censé retomber sur la suite de bruits blancs originelle.

data_schularick <- read_excel("C:/users/fkraus/Desktop/data_schularick.xlsx") 

donnees_france <- data_schularick %>%
  filter(country=="France")%>%
  dplyr::select(year, unemp, cpi) %>%
  mutate(inflation = (log(cpi) - lag(log(cpi)))*100)%>%
  mutate(dinflation = inflation - lag(inflation))%>%
  mutate(dunemp = unemp - lag(unemp))%>%
  na.omit()%>%
  filter(year > 1950)

head(donnees_france)

require(fixest)
reg <- feols(inflation ~ dunemp, data=donnees_france)
etable(reg)

                              reg
Dependent Var.:         inflation
                                 
Constant        4.060*** (0.4607)
dunemp            1.663* (0.7335)
_______________ _________________
S.E. type                     IID
Observations                   70
R2                        0.07031
Adj. R2                   0.05664
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analyse des résidus

On peut comparer la variable dépendante avec les valeurs prédites par le modèles :

fitted_val <- reg$fitted.values
residuals <- reg$residuals

donnees_france_2 <- donnees_france %>%
  cbind(fitted_val)%>%
  cbind(residuals)
ggplot(donnees_france_2, aes(x=year, y=inflation))+
  geom_line(lwd=1)+
  geom_line(aes(y=fitted_val), color="red", lty=2, lwd=1)+ 
  theme_bw()

On peut voir que notre modèle, bien qu’il utilise des variables correctement stationnarisées, ne représente par correctement la dynamique de la série : les fitted values ne suivent pas bien la série d’origine, particulièrement avant les années 1990.

Homoscedasticité

Pour vérifier l’homoscédasticité (=variance constante dans le temps), on peut déjà fait un graphique des résidus dans le temps.

ggplot(donnees_france_2, aes(x=year, y=residuals))+
  geom_line(lwd=1)+
  geom_hline(yintercept = 0, lty=2)+
  theme_bw()

On peut voir que la variance des résidus change drastiquement à partir de la fin des années 1980, ce qui pourrait indiquer que les résidus sont hétéroscédastiques. Cependant, l’analyse visuelle n’est pas un test formel. Pour vérifier l’homoscédasticité ou non des résidus, on doit utiliser le test de Breusch-Pagan (BP). Le test BP a pour hypothèse nulle l’hétéroscédasticité. Si la p-value est significative, on rejette donc l’hypothèse nulle et on conclut de l’homoscédasticité.

require(lmtest)
bptest(inflation ~ dunemp, data=donnees_france) 

    studentized Breusch-Pagan test

data:  inflation ~ dunemp
BP = 4.3177, df = 1, p-value = 0.03772

Ici, la p-value est significative, donc on conclut que les résidus de la régression sont hétéroscédastiques.

Normalité des résidus

ggplot(donnees_france_2, aes(residuals))+geom_density(fill="grey60", bounds=c(-Inf, Inf))+theme_bw()

Les résidus sont bien distribués autour de zéro, et la moyenne est relativement proche de zéro. En revanche, on peut voir que la queue de distribution à droite est relativement large, donc on tend à sous-estimer grandement le modèle à certains moments. La distribution n’est qu’une vérification visuelle, pas un test statistique. Pour un test formel, il faut se tourner vers le test de Shapiro qui teste la normalité d’une série. Le test de Shapiro a pour hypothèse nulle la normalité des résidus. Si la p-value est significative, on rejette donc l’hypothèse nulle, et on conclut de la non-normalité de la série.

shapiro.test(donnees_france_2$residuals)

    Shapiro-Wilk normality test

data:  donnees_france_2$residuals
W = 0.92324, p-value = 0.0003666

Ici, la p-value est significative, donc les résidus ne suivent pas une loi normale selon le test de Shapiro.

Indépendance

D’après l’hypothèse d’indépendance des résidus, les résidus d’une période ne doivent pas être expliqués par les résidus précédents. On peut tester cette hypothèse avec un simple autocorrélogramme :

acf(donnees_france_2$residuals)

Ici, on observe que les résidus ne semblent pas être indépendants, car les résidus passés expliquent de manière systématique les résidus actuels. Après l’analyse visuel, on peut faire le test formel. Le test Ljung-Box effectue le test d’indépendance des résidus, et a pour hypothèse nulle l’indépendance des résidus. Si la p-value est significative, on rejette l’hypothèse nulle et on conclut de l’autocorrélation (donc, l’absence d’indépendance) des résidus.

Box.test(donnees_france_2$residuals, lag = 10, type = "Ljung-Box")

    Box-Ljung test

data:  donnees_france_2$residuals
X-squared = 87.576, df = 10, p-value = 1.621e-14

Ici, il semble que nos résidus ne sont pas indépendants car la p-value est significative.

Note : certains de ces tests sont résumés dans la fonction checkresiduals() du package forecast:

require(forecast)

Le chargement a nécessité le package : forecast
Avis : le package ‘forecast’ a été compilé avec la version R 4.5.1Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

checkresiduals(reg)

    Ljung-Box test

data:  Residuals from feols
Q* = 87.576, df = 10, p-value = 1.621e-14

Model df: 0.   Total lags used: 10

Ainsi, dans l’ensemble, notre modèle est peu performant dans l’explication de la dynamique de l’inflation. En revanche, il explique correctement ce qu’il explique ! C’est un bon début, car cela veut dire qu’on peut maintenant ajouter d’autres variables explicatives dans le modèles (en testant, à chaque ajout, les effets sur les résidus pour vérifier qu’on ne crée par de problèmes qui seraient détectés dans les résidus) pour améliorer le pouvoir explicatif, tout en respectant le principe de parsimonie : on ne cherche pas à expliquer la totalité des variations de \(Y_t\), simplement à l’expliquer suffisamment pour que l’estimation de l’impact de \(\beta\) soit correct.

Overfitting

Parfois, on peut ajouter des variables qui vont créer un biais de collinéarité : les variables explicatives peuvent s’influencer entre elles, ou avoir un effet circulaire avec la variable dépendante, ce qui va créer un problème dans la régression.

reg_colin <- feols(inflation ~dunemp + lag(inflation), data= donnees_france)
etable(reg_colin)

                             reg_colin
Dependent Var.:              inflation
                                      
Constant              0e-16 (5.11e-17)
dunemp          4.44e-16*** (5.77e-17)
lag(inflation)         1*** (9.19e-18)
_______________ ______________________
S.E. type                          IID
Observations                        70
R2                                   1
Adj. R2                              1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Et quand on regard l’écart entre les fitted valued et la variable dépendante, on obtient :

fitted_val_colin <- reg_colin$fitted.values
residuals_colin <- reg_colin$residuals

donnees_france_3 <- donnees_france %>%
  cbind(fitted_val_colin)%>%
  cbind(residuals_colin)
ggplot(donnees_france_3, aes(x=year, y=inflation))+
  geom_line(lwd=1)+
  geom_line(aes(y=fitted_val_colin), color="red", lty=2, lwd=1)+ 
  theme_bw()

Ici, on peut être tenté de dire que notre modèle est parfait parce qu’il explique bien notre variable dépendante. Mais ici, on a un problème de collinéarité : lag(dinflation) explique quasi-parfaitement dinflation, et l’ajouter dans le modèle créé de la collinéarité parfaite.

On peut tester la collinéarité avec un test VIF

require(car)

Le chargement a nécessité le package : car
Le chargement a nécessité le package : carData

Attachement du package : ‘car’

L'objet suivant est masqué depuis ‘package:dplyr’:

    recode

L'objet suivant est masqué depuis ‘package:purrr’:

    some

vif(reg_colin)

                  GVIF Df GVIF^(1/(2*Df))
dunemp         1.07563  0             Inf
lag(inflation) 1.07563  0             Inf

Ici, le GVIF corrigé (colonne à droite) tend vers l’infini. On essaie en général de garder un VIF au maximum a 6 ou 7, au delà il y a un problème.

Donc ici, ajouter le lag de l’inflation va créer plus de problème qu’il n’en résout. Cependant, cette régression a pour mérite de montrer que l’inflation, même si elle est considérée comme étant stationnaire, continue d’avoir une structure auto-régressive importante, au point où si on la néglige, on perd un pouvoir explicatif très important. C’est la raison pour laquelle il faut s’intéresser à la structure Auto-Régressive (AR) des séries temporelles, à leur degré d’intégration (I) et à leur moyenne mobile (Moving Average - MA), dans le cadre des modèles ARIMA que l’on verra dans la séance 4.