R Notebook

setwd(dir = "~/Desktop/R/R2")

require(readxl)
require(tidyverse)
require(ggplot2)
require(urca)
require(fixest)
require(lmtest)
require(forecast)

data_schularick <- read_excel("~/Desktop/R/data_schularick.xlsx")

dk <- data_schularick %>%
  filter(country == "Denmark") %>%
  dplyr::select(year, cpi, xrusd) %>%
  filter(year > 1960)

head(dk)

PARTIE 1: RÉGRESSION LINÉTAIRE

Introduction

---

L’objectif dans cette première partie de projet, est d’analyser le pass-through du taux de change sur l’inflation au Danemark en étudiant la manière dont les variations du taux de change se transmettent aux prix à la consommation.

Pour répondre à cette question, nous étudions d’abord le comportement des séries CPI et XRUSD, puis nous vérifions leur stationnarité à l’aide de tests ADF. Une fois les séries correctement stationnarisées, nous analysons leur relation à travers un nuage de points et une régression linéaire. Enfin, nous vérifions si les hypothèses classiques du modèle MCO sont respectées grâce à l’étude des résidus afin de savoir si notre estimation statistique est pertinente.

I. Analyse des séries CPI et XRUSD

---

On utulisera la base de donnée datascularick sur laquelle dispose les données pour notre pays, le Danemark. Nous allons donc étudier l’indice des prix à la consommation et le taux de change de la monnaie locale vis à vis du dollar sur une période post 1960.

Rappelons qu’avant toute régression linéaire et/ou analyse économétrique, il faut s’assurer que les variables soient stationnaires. Nous allons donc ici, par représentation graphique, émettre une première hypothèse sur leur stationnarité ou non, avant de l’affirmer ou de l’infirmer via des tests ADF.

1. Étude de la variable CPI

A. Visualisation du CPI

Nous commencerons par analyser la premiere variable: le CPI

plot(dk$year, dk$cpi, type = "l",
  main = "CPI Denmark",
  xlab = "Année", ylab = "Indice de prix")

Tout d’abord, l’indice des prix à la consommation montre une tendance haussière très marqué durant la période 1960-2020. Il n’y a aucun fluctuation autour de la moyenne, au contraire la trajectoire est ascendante.

Rappelons également qu’une série DS est définie comme une série qui présente une racine unitaire, une série dont le niveau dépend fortement de ses valeurs passées. Ce type de série ne revient pas vers sa moyenne mais s’en éloigne durablement ce qui crée une tendance de long terme. C’est exactement que l’on peut remarquer pour notre variable.

Ainsi nous pouvons affirmer que visuellement, le CPI n’est pas stationnaire et est probablement un DS.

De plus, un processus DS étant défini par la présence d’une racine unitaire, on s’attend alors à ce que les valeurs passées influencent fortement les valeurs présentes = autocorrélation

On peut le vérifier avec l’autocorrélogramme:

acf(dk$cpi)

Chaque tiret vertical représente l’influence d’un délai (lag) particulier sur l’évolution de la série à la période \(t\). Le premier tiret peut être négligé ici, c’est l’influence de la série sur elle-même à la même période.

En l’absence d’autocorrélation en revanche, le reste des tirets est censé être non-significatif (donc les tirets ne doivent pas dépasser la ligne horizontale bleue). Ici, on voit bien que les délais jusqu’à 15 périodes (donc 15 ans avant) influencent la série à la période \(t\). C’est une indications très forte que la série n’est pas stationnaire.

Visuellement, on conclu que la variable CPI serait probablement non stationnaire, mais qu’en est il réellement ? Pour le verifier formellement, nous devons utiliser les tests ADF qui nous permettra de savoir si la variable est un TS ou un DS.

B. Tests de Stationnarité (ADF) sur CPI

L’analyse graphique du CPI révèle un trend haussier ce qui suggère une non-stationnarité en niveau. Nous devons tester les trois tests ADF (trend, drift et none). Pour déterminer la nature de la non-stationnarité du CPI, on commence par tester la présence d’une tendance déterministe en utilisant le modèle “trend”.

Rappel 1:

H0: la serie n’est pas stationnaire - Si t-stat < valeur critique → on rejette H0 - Si t-stat > valeur critique → on ne rejette pas H0

summary(ur.df(dk$cpi, type = "trend"))

############################################### 
# 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.86447 -0.48159 -0.07759  0.62536  2.82975 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.76188    0.35925   2.121   0.0386 *  
z.lag.1     -0.02557    0.01911  -1.338   0.1866    
tt           0.06771    0.05727   1.182   0.2423    
z.diff.lag   0.79734    0.08565   9.309  8.1e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9837 on 54 degrees of freedom
Multiple R-squared:  0.6258,    Adjusted R-squared:  0.605 
F-statistic: 30.11 on 3 and 54 DF,  p-value: 1.411e-11


Value of test-statistic is: -1.3378 2.5808 1.4574 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

Ici, aucune tendance n’est détectée, “tt” (0.2423) n’étant pas significative nous pouvons affirmer que le CPI n’est pas un trend.

Si aucune tendance n’est détectée, on teste ensuite le modèle avec drift afin de tester si une constante affecte la dynamique de la série.

summary(ur.df(dk$cpi, type = "drift"))

############################################### 
# 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.8107 -0.5295 -0.1019  0.5483  2.8707 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.852702   0.352201   2.421   0.0188 *  
z.lag.1     -0.003183   0.002593  -1.227   0.2249    
z.diff.lag   0.776554   0.084130   9.230 9.07e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9873 on 55 degrees of freedom
Multiple R-squared:  0.6162,    Adjusted R-squared:  0.6022 
F-statistic: 44.14 on 2 and 55 DF,  p-value: 3.667e-12


Value of test-statistic is: -1.2274 3.1497 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1  6.70  4.71  3.86

Rappel 2:

Pr(>|t|) = p-value du t-test -> indique si le coefficient est statistiquement significatif (donc présence de drift) z.diff.lag = présence d’autocorrélation ?? t stat = à comparer aux valeurs critiques pour conclure sur la stationnarité

Dans ce modèle, la significativité du drift est évaluée par la constante. Celle-ci est significative (p = 0.0188) ce qui indique qu’un drift est présent. Ainsi, lorsque le drift est détecté, on doit ensuite vérifier si l’autocorrélation est significative. Pour cela, on regarde la variable “z.diff.lag” = 9.07e-13 ** qui est siginificative. -> la série dépend de ses propres valeurs passées.

En plus d’un t-stat superieur aux valeurs critiques (on ne rejette pas H0), dans ce cas nous pouvons affirmer que le CPI est un DS.

2. Étude de la variable XRUSD

A. Visualisation de XRUSD

Visualisation de la variable XRUSD:

plot(dk$year, dk$xrusd, type = "l",
     main = "XRUSD",
     xlab = "Année", ylab = "DKK par USD")

Concernant le taux de change XRUSD, il n’y a pas de trend déterminé mais nous pouvons affirmer que la série fluctue autour d’une moyenne. Or, nous observons que sur certaines périodes notamment en 1985 et en 2000, il y a des fluctuations beaucoup plus fortes nous laissant penser que la variance peut changer dans le temps = variance n’est pas contante (même si l’on observe un retour à la “moyenne” assez vite)

Ce comportement est incompatible avec la stationnarité, car rappelons le; Une série stationnaire est une série dont les propriétés statistiques – moyenne, variance et covariance – sont constantes dans le temps. Formellement, on doit avoir : (Y_t)=, (Y_t)=^2, et (Y_t,Y_{t-h})=(h). Toute violation de ces conditions implique une non-stationnarité.

Cependant ici, nous mesurerons nos propos en affirmant qu’il y a effectivement des chocs, qu’il y a effectiviement quelques amplitudes inégales, pour autant la série revient plus ou moins autour de sa moyenne après ces derniers donc nous ne pouvons pas rejeter visuellement une stationnarité.

B. Tests de Stationnarité (ADF) sur XRUSD

Dans un second temps, afin de vérifier formellement la stationnarité de nos séries, nous appliquons le test de Dickey-Fuller augmenté (ADF). L’idée est de savoir si la série est stationnaire ou non.

Nous commençons par estimer le modèle avec tendance déterministe (type = “trend”), puis, en l’absence de tendance significative, nous testons le modèle avec constante uniquement (type = “drift”), et enfin, si la constante elle-même n’est pas significative, nous estimons le modèle sans constante (type = “none”).

acf(dk$xrusd)

Contrairement au CPI, cette autocorrélation décroît progressivement et disparaît rapidement après les premiers lags. On n’observe pas une autocorrélation sur toute la série. Visuellement, cela suggère que XRUSD pourrait être stationnaire ou proche de l’être car la dépendance existe mais reste modérée et s’éteint rapidement.

Pour trancher formellement, il est nécessaire d’avoir recours aux tests ADF qui permettront de déterminer si XRUSD suit un processus TS ou DS.

Même chose pour XRUSD, on a tester les trois tests ADF (trend, drift, none):

summary(ur.df(dk$xrusd, type="trend"))

############################################### 
# 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 
-1.42662 -0.44158 -0.02216  0.37511  1.74513 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.352715   0.618315   3.805 0.000364 ***
z.lag.1     -0.325384   0.082166  -3.960 0.000221 ***
tt          -0.006841   0.005317  -1.286 0.203762    
z.diff.lag   0.441086   0.123031   3.585 0.000724 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.648 on 54 degrees of freedom
Multiple R-squared:  0.2839,    Adjusted R-squared:  0.2441 
F-statistic: 7.137 on 3 and 54 DF,  p-value: 0.000401


Value of test-statistic is: -3.9601 5.2414 7.8494 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

Pour “trend”, la tendance tt n’est pas significative (t-stat = –1.286, p = 0.2038). On en conclut qu’il n’y a pas de trend déterministe et que le processus n’est pas un TS.

summary(ur.df(dk$xrusd, type="drift"))

############################################### 
# 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 
-1.50179 -0.46566  0.05874  0.30895  1.72165 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.93817    0.53085   3.651 0.000583 ***
z.lag.1     -0.29433    0.07901  -3.725 0.000461 ***
z.diff.lag   0.42369    0.12301   3.444 0.001103 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6518 on 55 degrees of freedom
Multiple R-squared:  0.262, Adjusted R-squared:  0.2351 
F-statistic: 9.761 on 2 and 55 DF,  p-value: 0.0002356


Value of test-statistic is: -3.7254 6.9518 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1  6.70  4.71  3.86

Pour ce test de regression drift, la constante est significative (p < 0.01), ce qui révèle la présence d’un drift. En plus de ça, il y a de l’autocorrélation significative 5% (0.001103 **). Cependant, la précense d’un drift + autocorrélation n’implique pas directement la présence d’un DS. Si le drift est significatif avec une autocorrélation significative (comme ici) mais que t-stat = –3.7254 est inférieure à la valeur critique à 1 % (–3.51), cela conduit à rejeter très fortement l’H0.

Ainsi, bien que le graphique du taux de change XRUSD laisse entrevoir des variations d’amplitude inégales, XRUSD est stationnaire.

Pour résumer, nous avons donc:

-CPI: DS -XRUSD: Stationnaire

Afin d’effectuer notre regression linéaire, il convient ici de stationnariser la variable “CPI”.

II. Stationnarisation des séries CPI et XRUDS

---

1. Mise en différence

Dans le cas des processus stationnaires en différence (DS), il suffit de passer la série en différence \(\Delta Y_t = Y_t - Y_{t-1}\) pour la rendre stationnaire.

dk <- dk %>%
  mutate(dcpi = cpi - lag(cpi)) %>%
  na.omit()

2. Vérification de la Stationnarité

Après avoir identifié que le CPI est un processus DS et l’avoir “stationnarisé” en appliquant une différence première, il faut vérifier si cette transformation a bien rendu la série stationnaire. Comme la série différenciée ne contient ni tendance déterministe ni constante, le test ADF approprié est le modèle “none” qui teste uniquement la présence d’une racine unitaire (stationnaire ou non).

dk <- dk %>%
  mutate(dcpi = c(NA, diff(cpi))) %>%
  na.omit()

summary(ur.df(dk$dcpi, type = "none"))

############################################### 
# 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.6849 -0.3856  0.1642  0.6789  2.8243 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)  
z.lag.1    -0.04604    0.04579  -1.006   0.3191  
z.diff.lag -0.22395    0.13294  -1.685   0.0978 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.024 on 54 degrees of freedom
Multiple R-squared:  0.07814,   Adjusted R-squared:  0.04399 
F-statistic: 2.288 on 2 and 54 DF,  p-value: 0.1112


Value of test-statistic is: -1.0055 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

On remarque encore une fois que la t-stat est superieur aux valeurs critiques, donc la difference première n’a pas suffit pour rendre CPI stationnaire = il faut passer en difference seconde.

dk <- dk %>%
  mutate(
    dcpi  = c(NA, diff(cpi)),    
    ddcpi = c(NA, diff(dcpi))    
  )

Il faut supprimer les NA (valeurs manquantes) générées par les différences successives puisqu’une régression ou un test ADF ne peut pas être estimé avec des observations manquantes. Nous utilisons donc la série sans NA avec “ddcpi_clean <- na.omit(dk$ddcpi)” afin de tester correctement la stationnarité.

ddcpi_clean <- na.omit(dk$ddcpi)
summary(ur.df(ddcpi_clean, type = "none"))

############################################### 
# 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.80753 -0.57796  0.05711  0.58136  2.73969 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1     -1.3745     0.2180  -6.304 6.35e-08 ***
z.diff.lag   0.1017     0.1381   0.736    0.465    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.046 on 52 degrees of freedom
Multiple R-squared:  0.6266,    Adjusted R-squared:  0.6122 
F-statistic: 43.63 on 2 and 52 DF,  p-value: 7.534e-12


Value of test-statistic is: -6.304 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Après avoir établi que le CPI était non-stationnaire en niveau et en première différence, nous avons appliqué une seconde différence. Le test ADF “none” sur la série ddCPI montre un t-stat de –6.20 largement inférieur aux valeurs critiques (–2.6 au seuil de 1 %). Nous rejetons donc fortement l’hypothèse nulle de racine unitaire. Ainsi ddCPI est stationnaire. Enfin, nous allons représenter le CPI avant et après la difference seconde:

# Graphique en niveau: CPI non stationnaire
ggplot(dk, aes(x = year, y = cpi)) +
  geom_line(lwd = 1) +
  labs(title = "CPI Denmark - Niveau",
       x = "Année", y = "Indice des prix") +
  theme_bw()

# Graphique de CPI après la difference seconde
ggplot(dk, aes(x = year, y = ddcpi)) +
  geom_line() +
  labs(title = "CPI Denmark - Différence seconde (Δ²CPI)",
       x = "Année", y = "Δ²CPI") +
  theme_bw()

->Le CPI est donc un processus intégré d’ordre 2, I(2).

III. Relation entre les deux variables CPI et XRUSD

---

1. Nuage de points

Pour estimer graphiquement le lien entre les variables, il est nécessaire de les visualier via un nuage de points.

ggplot(data = dk, aes(x = xrusd, y = ddcpi)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(
    title = "d²CPI en fonction de XRUSD",
    x = "XRUSD",
    y = "Δ²CPI"
  ) +
  theme_bw()

Le nuage de points représentant la relation entre le taux de change XRUSD et la difference seconde du CPI (d²CPI) ne laisse apparaître aucune relation linéaire nette. Les points sont largement dispersés sans structure claire ce qui suggère qu’il n’existe pas de corrélation visuelle entre les deux variables. La droite de régression présente une pente légèrement négative (si XRUSD s’apprécie alors CPI diminue légèrement)

Une régression linéaire est nécessaire pour confirmer ou infirmer l’existence d’un pass-through significatif entre XRUSD et l’inflation.

2. Regréssion linéaire

reg_lin <- feols(ddcpi ~ xrusd, data = dk)
etable(reg_lin)

                         reg_lin
Dependent Var.:            ddcpi
                                
Constant         0.5015 (0.8181)
xrusd           -0.0757 (0.1220)
_______________ ________________
S.E. type                    IID
Observations                  56
R2                       0.00708
Adj. R2                 -0.01131
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On voit ici que le coefficient associé au taux de change XRUSD n’est pas significatif à 5%. Il n’existe pas de relation statistiquement détectable entre les variations du taux de change et les variations du niveau des prix stationnarisé. Dans ce cadre, le pass-through du taux de change vers l’inflation n’apparaît pas, les fluctuations du taux de change ne se transmettent donc pas directement aux prix à la consommation.

IV. Vérification des hypothèses des MCO

---

1. Analyse des valeurs ajustées et des résidus

Les résidus représentent l’erreur de prédiction du modèle : ils mesurent à quel point la régression parvient — ou échoue — à reproduire la dynamique réelle de la variable expliquée. Dans l’idéal, un bon modèle MCO produit des résidus « neutres » : en moyenne proches de zéro, sans tendance particulière, avec une variance constante et aucune dépendance dans le temps. Autrement dit, si le modèle est correctement spécifié, les résidus doivent se comporter comme un bruit blanc. L’analyse des résidus permet donc de vérifier si les hypothèses classiques du MCO sont respectées, et donc si la régression peut être considérée comme fiable et interprétable.

On peut imaginer qu’une série temporelle est, à l’origine, une simple suite de bruits blancs perturbée au fil du temps par divers chocs exogènes. Le rôle de la régression est précisément de retirer l’influence de ces chocs : si le modèle est correctement spécifié, les résidus doivent donc retrouver le comportement d’un bruit blanc =série sans structure, centrée, indépendante et de variance constante.

Pour évaluer la qualité de notre modèle, nous comparons les valeurs observées de Δ²CPI aux valeurs ajustées par la régression afin de visualiser la capacité explicative réelle du modèle

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

dk_clean <- dk %>%
  dplyr::select(year, xrusd, ddcpi) %>%
  na.omit()

dk2 <- dk_clean %>%
  cbind(fitted_val) %>%
  cbind(residuals)

ggplot(dk2, aes(x = year, y = ddcpi)) +
  geom_line(lwd = 1) +
  geom_line(aes(y = fitted_val), color = "red", lty = 2, lwd = 1) +
  labs(
    title = "Valeurs observées vs valeurs ajustées",
    x = "Année",
    y = "Δ²CPI"
  ) +
  theme_bw()

Le graphique montre que les valeurs ajustées par le modèle restent presque constantes tandis que les valeurs observées de Δ²CPI fluctuent fortement. ->le modèle explique très peu la dynamique réelle de la série.

Nous devons maintenant vérifier si les hypothèses classiques du modèle des MCO sont respectées afin de garantir la validité statistique de notre estimation.

2. Tests de diagnostic des résidus

dk_clean <- dk %>%
  dplyr::select(year, xrusd, ddcpi) %>%
  na.omit()

dk2 <- dk_clean %>%
  mutate(
    residuals     = reg_lin$residuals,
    fitted_values = reg_lin$fitted.values
  )

head(dk2)

Juste avant de réaliser un test formel, nous commençons par examiner graphiquement les résidus au cours du temps afin d’évaluer visuellement si leur variance est constante. -Une variance stable indiquerait une homoscédasticité -des variations marquées de l’amplitude des résidus suggéreraient une hétéroscédasticité potentielle.

ggplot(dk2, aes(x = year, y = residuals)) +
  geom_line(lwd = 1) +
  geom_hline(yintercept = 0, lty = 2) +
  labs(title = "Résidus de la régression dans le temps",
       x = "Année", y = "Résidus") +
  theme_bw()

On peut voir que la variance des résidues modérée dans les années 1970–1980, mais à partir du début des années 2008-2010, l’amplitude des résidus augmente fortements. Cette évolution amène à penser que les résidus pourraient être 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 utilise donc le test de Breusch-Pagan.

1. Hétéroscédasticité: Breusch-Pagan

bptest(ddcpi ~ xrusd, data = dk)

    studentized Breusch-Pagan test

data:  ddcpi ~ xrusd
BP = 4.9343, df = 1, p-value = 0.02633

Dans notre cas, le test BP donne une statistique de 4.9343 pour 1 degré de liberté avec une p-value = 0.02633. Comme la p-value est significative à 5%, on rejette l’hypothèse nulle d’homoscédasticité et on conclut que les résidus sont hétéroscédastiques.

2. Normalité: Shapiro-Wilk

ggplot(dk2, 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. 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(dk2$residuals)

    Shapiro-Wilk normality test

data:  dk2$residuals
W = 0.94218, p-value = 0.009632

La p-value est significative, donc les résidus ne suivent pas une loi normale selon le test de Shapiro. Cela signifie que les erreurs du modèle ne sont pas distribuées de façon symétrique autour de zéro.

3. Indépendance: Ljung-Box

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(dk2$residuals)

Ici, on observe un pic significatif au lag 1. Les autres lags ne sont pas significatifs. Après l’analyse visuelle, il convient d’utiliser 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(dk2$residuals, lag = 10, type = "Ljung-Box")

    Box-Ljung test

data:  dk2$residuals
X-squared = 8.401, df = 10, p-value = 0.5897

La p-value du test de Ljung-Box n’est pas significative. On ne rejette donc pas l’hypothèse nulle d’indépendance des résidus. Ce test confirme ce que nous avons vu; les résidus apparaissent indépendants et il n’y a pas d’autocorrélation entre eux.

Concernant les MCO, plusieurs hypothèses sont violées: -les résidus sont hétéroscédastiques -ils ne suivent pas une loi normale -mais ils apparaissent indépendants

V. Conclusion

---

Dans l’ensemble, notre modèle est peu performant dans l’explication de la dynamique de l’inflation au Danemark. Le coefficient du taux de change (β) n’est pas significatif donc les variations du taux de change n’expliquent pas les variations de Δ²CPI. Pour répondre à notre question, on ne détecte aucun pass-through dans notre modèle: les variations du taux de change n’ont aucun impact significatif sur le CPI ###= le pass-through est nul.

Enfin, pour clôre cette première partie, nous comparerons notre étude à celle de la BIS. Dans notre étude nous utilisons une régression linéaire simple en séries temporelles (plusieurs “t” pour un seul individu, ici le Danemark), où l’inflation stationnarisée (Δ²CPI) est expliquée par le taux de change XRUSD. Il s’agit donc d’un modèle univarié qui tient compte d’une seule relation directe. Dans l’étude de la BIS, les auteurs utilisent des modèles multivariés estimés en panel. Leur approche permet donc de mieux expliquer le “pass-through”, avec une richesse de données beaucoup plus importante.

Pour rapprocher notre étude de la leur, il faudrait avant tout travailler avec des données plus fréquentes ce qui permettrait fortement d’enrichir notre analyse du pass-through.

---

---

PARTIE 2: PRÉVISIONS DU PIB

---

Pour continuer maintenant, nous devons prévoir la croissance future de notre pays d’intérêt sur les 5 prochaines années. Pour cela, vous vous servirez des modèles ARMA afin de faire des prévisions. Vous choisirez le modèle optimal et effectuerez une prévision unique de la croissance du PIB. Vous commenterez ensuite les résultats obtenus, et conseillerez le responsable politique sur la politique appropriée. Commentez également les limites liées à votre prévision.

donnees_dk <- data_schularick %>%
  filter(country=="Denmark")%>%
  select(year, gdp)%>%
  mutate(gdp_growth = (log(gdp) - lag(log(gdp)))*100   )%>%
  na.omit()%>%
  filter(year >= 1950 & year < 2020)
head(donnees_dk)

On commence à partir de 1950 car avant c’était la Seconde Guerre Mondiale, et on a retiré 2020 car c’est une année exceptionnelle, avec un choc purement exogène (donc impossible de le relier à des variations passées du PIB).

donnees_dk %>%
  ggplot(aes(x=year, y=gdp_growth))+
  geom_line()

Nous remarquons qu’à partir des années 60 jusqu’en 1985 environ, le PIB du Danemak est a son apogé avec une croissance moyenne par an d’envirion 10%. Cependant après 1985, il y a une chute drastique du PIB chutant à environ 2.5%. Nous observons un pic du PIB qui est baissier juste avant les années 2010 qui est dû à la crise des subprimes.

df <- donnees_dk %>% arrange(year)

y_ts <- ts(
  df$gdp_growth,
  start     = min(df$year),
  frequency = 1        
)

ici, on met les données en ordre chronologique puis on convertit notre variable du PIB en une serie temporelle annuelle afin de pouvoir appliquer le modèle temporel ARIMA.

require(forecast)
fit_arima <- auto.arima(y_ts)
fit_arima

Series: y_ts 
ARIMA(0,1,1) 

Coefficients:
          ma1
      -0.6236
s.e.   0.0880

sigma^2 = 6.778:  log likelihood = -163.67
AIC=331.35   AICc=331.53   BIC=335.82

On estime ensuite le modèle ARIMA. La sélection automatique du modèle nous donne le modèle ARIMA(0,1,1) -> MA(1) comme celui minimisant le critère AIC. Ce modèle est donc considéré comme le plus adapté pour capturer la dynamique du PIB danois.

On peut donc regarder la dynamique des résidus du modèle MA(1) qui correspondent à la composante du PIB non expliquée par le modèle et les comparer aux valeurs observées pour évaluer la qualité de l’ajustement.

Juste avant, il faut extraire les valeurs ajustées (fitted values) du modèle ARIMA correspondant aux valeurs que le modèle ARIMA aurait prédites pour chaque année de la période d’estimation. Cela permet de comparer ce que le modèle prévoit contre ce qui s’est réellement passé.

fitted_arima <- fitted(fit_arima)

Nous construisons ensuite un tableau comprenant les valeurs observées du PIB et les valeurs ajustées par ARIMA afin de comparer la qualité d’ajustement du modèle sur l’échantillon:

df_in_sample <- data.frame(
  year         = as.numeric(time(y_ts)),
  gdp_growth   = as.numeric(y_ts),
  arima_fitted = as.numeric(fitted_arima)
)
head(df_in_sample)

Ainsi nous obtenons:

ggplot(df_in_sample, aes(x = year)) +
  geom_line(aes(y = gdp_growth), linetype = "solid") +
  geom_line(aes(y = arima_fitted), linetype = "dashed", color="red") +
  labs(
    title = "In-sample: actual vs ARIMA fitted",
    y = "Value"
  )

Le graphique compare les valeurs observées du PIB (trait noir) aux valeurs ajustées par le modèle ARIMA (trait rouge). On voit que le modèle capture la tendance centrale de l’évolution générale du PIB, mais reste tout de même légèrement éloigné des fluctuations réelles. En effet, les valeurs ajustées apparaissent plus “lisses” et moins volatiles que les observations, exemple avec le pic baissier causé par la GCF.

On va, afin de réaliser nos prévisions, établir la taille de notre échantillon sur lequel on appliquera le modèle de prévision ARIMA. Ici, la taille de l’échantillon est de 15 ans.

n <- length(y_ts)
h <- 15  # tu fixes ici

Ici, on va venir diviser la partie en deux: “train” sert à estimer le modèle ARIMA, et “test” sert à vérifier la qualité des prévisions hors-échantillon.

y_train <- window(y_ts, end = time(y_ts)[n - h])
y_test  <- window(y_ts, start = time(y_ts)[n - h + 1])

On estime ensuite un modèle ARIMA seulement sur la partie train, puis on génère des prévisions sur les 15 années du test accompagnées des intervalles de confiance à 80 % et 95 %. Enfin, on construit un tableau comparant les valeurs réelles du PIB avec les prévisions ARIMA et leurs intervalles ce qui permet de visualiser et d’évaluer la performance prédictive du modèle.

fit_arima_oos <- auto.arima(y_train)
fc_arima_oos  <- forecast(fit_arima_oos, h = h, level = c(80, 95))

years_test <- as.numeric(time(y_test))

arima_insample <- data.frame(
  year           = years_test,
  gdp_growth     = as.numeric(y_test),
  arima_forecast = as.numeric(fc_arima_oos$mean),
  arima_lower_95 = fc_arima_oos$lower[, "95%"],
  arima_upper_95 = fc_arima_oos$upper[, "95%"]
)
arima_insample

NA

Représentation graphique:

ggplot(arima_insample, aes(x = year)) +
  geom_line(aes(y = gdp_growth), linetype = "solid") +
  geom_line(aes(y = arima_forecast), linetype = "dashed") +
  geom_ribbon(aes(ymin = arima_lower_95, ymax = arima_upper_95),
              alpha = 0.2) +
  labs(
    title = "Out-of-sample: ARIMA forecast vs actual",
    y = "Value"
  )

Le graphique compare les valeurs réelles du PIB aux prévisions du modèle ARIMA (courbe pointillée), accompagnées de leur intervalle de confiance à 95 % correspondant à la zone grisée.

Nous remarquons que le modèle prévoit une trajectoire relativement stable alors que les valeurs réelles présentent des fluctuations importantes notamment le pic négatif très marqué au moment de la crise financière de 2008. Ce choc écononomique n’est pas anticipé par le modèle ARIMA ce qui montre sa principale limite. Pour autant, ce modèle capture relativement bien la tendance générale mais ne peut prévoir à l’avance les crises (ce qui va de soi).

Enfin, nous allons procéder aux prévisions concernant la croissance du PIB danois sur les 5 prochaines années.

h <- 5
fc_arima  <- forecast(fit_arima, h = h)
plot(fc_arima)

Avec le modèle ARIMA, on peut anticiper que la croissance du PIB va rester entre à peu près -1/2% et 7-8% jusqu’en 2024, ce qui reste un intervalle très large. Dans la réalité, la crise du Covid a créé une baisse importante du PIB danois en 2020 avec une décroissance de -1.8%, contre un rebond assez important en 2021 s’élevant à 7.4%. Rien qu’avec ce rebond, on pourrait être tenté d’invalider nos prévisions. Finalement, une fois de plus les chocs imprévisibles démontrent les limites de ce modèle.

Concernant la mise en place d’une politique appropriée, au vu de nos estimations “plates”, il est difficile de conseiller sur une politique en particulier. Nous pouvons, par ailleurs, affirmer qu’en se basant uniquement sur nos prévisions, il n’y aurait a priori pas de crise à prévoir. (Rappelons que le modèle ARIMA ne prévoit pas les crises, ici le but est simplement de répondre à la question sur la politique appropriée en se basant simplement sur nos prévisions.)

---

---

PARTIE 3: MODÈLE À CORRECTION D’ERREURS

---

votre_pays <- "Denmark"

donnees_df <- data_schularick %>%
  filter(country == votre_pays) %>%
  filter(year > 1960) %>%
  select(year, revenue, expenditure) %>%
  mutate(
    revenue_growth     = c(NA, diff(log(revenue)) ),
    expenditure_growth = c(NA, diff(log(expenditure)) )
  ) %>%
  na.omit()

head(donnees_df)

Tout d’abord, nous tester si nos deux variables sont stationnaires.

summary(ur.df(donnees_df$revenue_growth, type="none"))

############################################### 
# 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 
-0.138896 -0.023615 -0.002589  0.046756  0.167424 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -0.11866    0.07170  -1.655    0.104    
z.diff.lag -0.68040    0.09532  -7.138 2.24e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05348 on 55 degrees of freedom
Multiple R-squared:  0.5599,    Adjusted R-squared:  0.5439 
F-statistic: 34.98 on 2 and 55 DF,  p-value: 1.577e-10


Value of test-statistic is: -1.655 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Le test ADF donne un tstat de –1.655. Comme –1.655 est supérieur à la valeur critique de 5% (–1.95) mais à peine plus faible que celle de 10% (–1.61), revenue_growth apparaît non stationnaire. Nous devons donc réaliser un second test ADF, cette fois-ci en passant en difference première:

summary(ur.df(diff(donnees_df$revenue_growth), type="none"))

############################################### 
# 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 
-0.147019 -0.032890 -0.005967  0.030594  0.140596 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1     -2.2542     0.2403  -9.382 6.23e-13 ***
z.diff.lag   0.2935     0.1292   2.272   0.0271 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05239 on 54 degrees of freedom
Multiple R-squared:  0.8794,    Adjusted R-squared:  0.8749 
F-statistic: 196.8 on 2 and 54 DF,  p-value: < 2.2e-16


Value of test-statistic is: -9.3815 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Le test ADF nous dévoile une tstat très inférieure aux valeurs critiques (t = –9.38 < –2.6 à 1 %). Ainsi la série différenciée est stationnaire. Nous obtenons donc: revenue_growth est un processus intégré d’ordre 1, I(1).

Ensuite, il convient ici de faire la même chose pour la variable “expenditure_growth”

summary(ur.df(donnees_df$expenditure_growth, type="none"))

############################################### 
# 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 
-0.070330 -0.021074  0.006501  0.018314  0.161241 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -0.08708    0.06252  -1.393 0.169267    
z.diff.lag -0.44478    0.12227  -3.638 0.000608 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04714 on 55 degrees of freedom
Multiple R-squared:  0.2481,    Adjusted R-squared:  0.2207 
F-statistic: 9.072 on 2 and 55 DF,  p-value: 0.0003935


Value of test-statistic is: -1.3929 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Le test ADF montre que la tstat est de –1.39 qui est supérieure aux seuils critiques, ce qui conduit à affirmer que la série n’est pas stationnaire et doit être différenciée pour le devenir.

summary(ur.df(diff(donnees_df$expenditure_growth), type="none"))

############################################### 
# 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 
-0.099124 -0.026948 -0.000348  0.015206  0.143290 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1     -1.7549     0.2362  -7.429 8.29e-10 ***
z.diff.lag   0.1774     0.1369   1.296      0.2    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04751 on 54 degrees of freedom
Multiple R-squared:  0.737, Adjusted R-squared:  0.7273 
F-statistic: 75.67 on 2 and 54 DF,  p-value: < 2.2e-16


Value of test-statistic is: -7.429 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Après passage en différence première, la tstat (–7.43) devient largement inférieure a ux valeurs critiques ce qui confirme que la difference premiere de notre série expenditure_growth est stationnaire.

Ainsi, les deux séries étudiées; revenue_growth et expenditure_growth, apparaissent non stationnaires en niveau car leurs tstat sont supérieures aux valeurs critiques.

En revanche, après passage en différence première les deux séries deviennent clairement stationnaires : -ADF(Δrevenue_growth) = –9.38 -ADF(Δexpenditure_growth) = –7.4 Ainsi, les deux variables sont intégrées d’ordre 1, I(1).

Les deux modèles optimaux pour nos variables sont des modèles ARMA avec un ordre d’intégration égal à 1. Il y a donc un risque de régression fallacieuse si on fait la régression classique. Etant donné la présence de cointégration dans notre modèle, on teste la validité du modèle à correction d’erreurs:

fs <- feols(revenue_growth ~ expenditure_growth, data = donnees_df)
donnees_df$resid <- residuals(fs)
summary(ur.df(donnees_df$resid, type = "none"))

############################################### 
# 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 
-0.140695 -0.031834 -0.004368  0.028508  0.155290 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
z.lag.1     -0.8870     0.1928  -4.601 2.52e-05 ***
z.diff.lag  -0.1349     0.1385  -0.974    0.334    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05276 on 55 degrees of freedom
Multiple R-squared:  0.4997,    Adjusted R-squared:  0.4815 
F-statistic: 27.46 on 2 and 55 DF,  p-value: 5.371e-09


Value of test-statistic is: -4.6005 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Comme la tstat est –4.60 < –2.6, nous pouvons affirmer que les résidus sont stationnaire et que donc il existe bien une relation de long terme entre revenue_growth et expenditure_growth.

Ainsi, comme les deux séries sont intégrées du même ordre et leurs résidus sont stationnaires on peut alors utiliser un modèle à correction d’erreur (ECM) à l’étape 2.

On crée ensuite la variable lresiduals qui représente l’écart de la relation long terme à la période précédente, car c’est cette déviation qui permet au modèle de mesurer la vitesse de retour vers l’équilibre. -> pour mesurer la vitesse à laquelle le système corrige les déséquilibres de long terme.

donnees_df_2 <- donnees_df %>%
  mutate(lresiduals = dplyr::lag(resid, 1))

Puis, en sachant que nos deux variables sont stationnaires une fois la difference premiere réalisée, il faut donc creer les “variables” différences premières des deux séries car le modèle à correction d’erreur nécessite de travailler avec des variables stationnaires en variation (ΔYₜ et ΔXₜ), afin de capturer correctement la dynamique de court terme.

Par ailleurs, puisque nos deux variables deviennent stationnaires après passage en différence première, il est nécessaire de créer les series en difference premiere “drevenue_growth” et “dexpenditure_growth” car le modèle à correction d’erreur doit utiliser nos variables stationnaires (ΔYₜ et ΔXₜ) afin de capturer la dynamique de court terme.

donnees_df_2 <- donnees_df_2 %>%
  mutate(
    d_revenue_growth     = c(NA, diff(revenue_growth)),
    d_expenditure_growth = c(NA, diff(expenditure_growth))
  ) %>%
  na.omit()

Ainsi on réalise donc le modèle à correction d’erreur permettant d’étudier simultanément la relation court et long terme (corréction des écarts) entre les variations des dépenses publiques et celles des revenus publics.

ecm <- feols(
  d_revenue_growth ~ d_expenditure_growth + lresiduals,
  data = donnees_df_2
)

etable(ecm)

                                    ecm
Dependent Var.:        d_revenue_growth
                                       
Constant               -0.0011 (0.0070)
d_expenditure_growth 0.6273*** (0.1302)
lresiduals           -1.006*** (0.1380)
____________________ __________________
S.E. type                           IID
Observations                         58
R2                              0.57356
Adj. R2                         0.55805
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Le coefficient de d_expenditure_growth 0.6273*** est positif et très significatif, cela signifie qu’à court terme une hausse des dépenses publiques s’accompagne d’une hausse des revenus publics. Le coefficient de lresiduals -1.006*** est négatif, très significatif et légèrement inférieur à –1 ce qui indique une forte force de rappel vers l’équilibre de long terme. Cela signifie que lorsque les revenus s’éloignent de leur trajectoire “théorique”, l’ajustement se fait rapidement.

Par ailleurs, les résidus de la première étape sont stationnaires, donc l’ECM est valide. De plus, puisque le coefficient de lagresiduals est négatif très significatif et compris entre –1 et 0 (ici nous dirons pour des raions de simplicité que -1.006*** est environ égal à -1), il est préférable d’utiliser le modèle à correction d’erreur plutôt qu’un OLS classique pour représenter la relation entre revenue_growth et expenditure_growth car l’ECM intègre à la fois l’ajustement court terme et la relation d’équilibre de long terme alors que l’OLS non.

Avant de terminer, il faut garantir la fiabilité de notre modèle à correction d’erreur en regardant si les résidus verifient trois hypothèses. Nous testerons donc: -l’homoscédasticité à l’aide du test de Breusch-Pagan qui permet de vérifier que la variance des résidus reste constante dans le temps.

-ensuite la normalité des résidus via le test Shapiro-Wilk qui teste la normalité d’une série

-Et enfin, l’indépendance des résidus grâce au test Ljung-Box afin de s’assurer qu’il n’y a pas d’autocorrélation “les résidus d’une période ne doivent pas être expliqués par les résidus précédents.”

res_ecm <- residuals(ecm)
bptest(d_revenue_growth ~ d_expenditure_growth + lresiduals, data = donnees_df_2)

    studentized Breusch-Pagan test

data:  d_revenue_growth ~ d_expenditure_growth + lresiduals
BP = 0.64861, df = 2, p-value = 0.723

Dans notre cas, le test BP donne une statistique de 0.64861 pour 1 degré de liberté avec une p-value = 0.723. Comme la p-value n’est pas significative, on ne rejette pas l’hypothèse nulle d’homoscédasticité et on conclut que les résidus sont homoscédastiques.

shapiro.test(res_ecm)

    Shapiro-Wilk normality test

data:  res_ecm
W = 0.97386, p-value = 0.2427

La p-value est significative, donc les résidus ne suivent pas une loi normale selon le test de Shapiro. Cela signifie que les erreurs du modèle ne sont pas distribuées de façon symétrique autour de zéro.

Box.test(res_ecm, lag = 10, type = "Ljung-Box")

    Box-Ljung test

data:  res_ecm
X-squared = 6.9619, df = 10, p-value = 0.729

La p-value du test de Ljung-Box n’est pas significative. On ne rejette donc pas l’hypothèse nulle d’indépendance des résidus. Les résidus apparaissent indépendants et il n’y a pas d’autocorrélation entre eux.

Pour terminer, le test de Breusch-Pagan indique une homoscédasticité des résidus (variance des erreurs reste stable dans le temps). En revanche, le test de Shapiro-Wilk montre que les résidus ne suivent pas une loi normale. Enfin, le test de Ljung-Box révèle que les résidus sont indépendants donc sans autocorrélation. Ainsi, malgré une non-normalité des erreurs l’ECM reste globalement correcte, car les propriétés essentielles de variance constante et d’absence de dépendance temporelle sont respectées.