import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd
from pandas import read_csv
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.tsa.arima_model import ARMA
from statsmodels.stats.diagnostic import acorr_ljungbox
EXERCICE 1¶
Question 1¶
"CH" est l'acronyme de "Conditional Heteroskedasticity", qui est une façon sophistiquée de dire que la variance conditionnelle $\text{Var}[X_t|X_{t-1},X_{t-2},\ldots]$ dépend du temps $t$. Si $(X_t)$ est un GARCH alors il est stationnaire et en particuler sa variance (tout court) ne dépend pas du temps. En vertu de la Proposition 6.2 du cours, $(X_t)$ est un bruit blanc, et donc son ACF est celle d'un bruit blanc : $\rho(h)=\boldsymbol 1_{h=0}$. Ce n'est pas un bruit blanc fort car la série $(X_t^2)$ n'est pas un bruit blanc (même faible).
Question 2¶
On peut générer une réalisation $x_1,\ldots,x_n$ d'un processus ARCH(1), $X_t=\epsilon_t\sigma_t$, de paramètres $\omega,\alpha$ :
- en simulant $\epsilon_1,\ldots,\epsilon_n$ i.i.d $N(0,1)$ gaussien
- en calculant $\sigma_1^2,\ldots,\sigma_{n}^2$ et $x_1,\ldots,x_{n}$ récursivement : on part de $\sigma_0:=0$ (et donc $x_0=0$) et puis pour tout $t\geq 0$,
- après, comme on l'a fait pour l'AR(1) dans le TP1, il vaut mieux simuler une trajectoire plus longue et supprimer une partie pour faire disparaître la dépendance en la condition initiale (arbitraire) $\sigma_0=0$. Par exemple, on peut simuler une trajectoire $x_1,\ldots,x_{2n}$ de longueur $2n$ et ne garder que les $n$ dernières valeurs $x_{n+1},\ldots,x_{2n}$.
def ARCH1(omega,alpha,n):
eps = npr.randn(2*n)
sigma2 = 0
ARCH = np.zeros(2*n)
for i in range(2*n-1):
sigma2 = omega + alpha * ARCH[i]**2
ARCH[i+1] = eps[i] * np.sqrt(sigma2)
return(ARCH[n:2*n])
Question 3¶
Remarque préliminaire : Comme on l'a vu dans le cours, la condition $\mathbb E[X_t^4]<\infty$ pour un ARCH(1) requiert que $0<\alpha<1/\sqrt{3}$ quand $(\epsilon_t)\sim \text{BB}(0,1)$ gaussien. En fait cette condition aide aussi à ce que la simulation de la question 2 soit représentative d'un ARCH(1) théorique. Ce qui suit dépendra donc de votre choix du paramètre $\alpha$, sachant que :
1/np.sqrt(3)
Xt = ARCH1(1,0.2,1000)
plt.plot(Xt)
plot_acf(Xt)
plt.show()
acorr_ljungbox(Xt,lags=range(60),boxpierce=True)[3] < 0.05
On peut raisonnablement supposer que la série $X_t$ est un bruit blanc. Etudions maintenant $X_t^2$.
plt.plot(Xt**2)
plot_acf(Xt**2)
plot_pacf(Xt**2)
plt.show()
acorr_ljungbox(Xt**2,lags=range(60),boxpierce=True)[3] < 0.05
L'hypothèse que $X_t^2$ soit un bruit blanc est clairement rejetée. De plus, la PACF laisse penser que $X_t^2$ est modélisable par un AR(1) avec un terme constant, ce qu'on vérifie (notez qu'on n'a pas rajouté trend='nc' comme on l'a fait jusqu'à présent) :
model = ARMA(Xt**2,order=(1,0)).fit()
acorr_ljungbox(model.resid,lags=range(60),boxpierce=True)[3] < 0.05
NB : On peut retrouver une approximation des coefficients $\omega,\alpha$ de l'ARCH(1) dans la modélisation AR(1) (pourquoi ?)
model.summary()
Notez qu'on peut isoler ces coefficients en utilisant : model.params
Comme $\sigma_{n+1}^2 = \omega + \alpha \, X_n^2$, on obtient l'estimateur de $\sigma_{1001}$ (imaginons un instant qu'on ne connaisse pas $\alpha,\omega$) :
np.sqrt(model.params[0] + model.params[1]*Xt[-1])
INTC = read_csv('INTC.csv').x
INTC.plot()
plt.show()
plot_acf(INTC)
plot_pacf(INTC)
plt.show()
acorr_ljungbox(INTC,lags=range(60),boxpierce=True)[3] < 0.05
Conclusion : On peut raisonnablement supposer que la série INTC est un bruit blanc.
Question 2¶
On étudie maintenant la série INTC$^2$
INTC2 = INTC**2
INTC2.plot()
plot_acf(INTC2)
plot_pacf(INTC2)
plt.show()
acorr_ljungbox(INTC2,lags=range(60),boxpierce=True)[3] < 0.05
L'hypothèse d'un bruit blanc pour INTC$^2$ est clairement rejetée.
Question 3¶
Essayons de modéliser INTC à l'aide d'un GARCH, c'est-à-dire de modéliser INTC2 avec un ARMA. D'après la PACF de INTC2 un AR(1) ne semble pas approprié pour modéliser cette série (pas de pic isolé en $h=1$), donc ARCH(1) ne semble pas approprié pour INTC (car on fait l'hypothèse que cette série à un moment d'ordre 4). On teste quelques modèles ARMA(p,q) avec petits $p,q$ :
model_20 = ARMA(INTC2,order=(2,0)).fit()
model_30 = ARMA(INTC2,order=(2,0)).fit()
model_40 = ARMA(INTC2,order=(2,0)).fit()
model_11 = ARMA(INTC2,order=(1,1)).fit()
model_21 = ARMA(INTC2,order=(2,1)).fit()
model_20.aic, model_30.aic, model_40.aic, model_11.aic, model_21.aic
model_20.bic, model_30.bic, model_40.bic, model_11.bic, model_21.bic
Le modèle de plus petit critère AIC (et BIC) est l'ARMA(1,1), dont on étudie les résidus maintenant.
residus = model_11.resid
acorr_ljungbox(residus,lags=range(60),boxpierce=True)[3] < 0.05
Conclusion : Les résidus ne semblent pas corrélés et le modèle GARCH(1,1) semble donc approprié pour INTC.
Question 4¶
- Si on note $(X_t)$ la série théorique de INTC, alors sous l'hypothèse du modèle GARCH(1,1) sa volatilité $\sigma_t$ satisfait : $$ \sigma_t^2 = \omega + \alpha X_{t-1}^2 +\beta \sigma_{t-1}^2 $$ pour des paramètres $\omega,\alpha,\beta >0$ et $\alpha+\beta<1$.
- En faisant de plus l'hypothèse que $\mathbb E[X_t^4]<\infty$, on sait que $(X_t^2)$ est un ARMA(1,1) avec une constante additive dont les paramètres sont en lien avec $\omega,\alpha,\beta$ :
$$ X_t^2 = \omega + \varphi X_{t-1}^2+\theta\eta_{t-1} $$ où $(\eta_t)\sim BB(0,\eta^2)$ et $$\varphi=\alpha+\beta,\qquad \theta=-\beta.$$
- On trouve les constantes $\omega,\varphi,\theta$ à l'aide de la modélisation faite en Question 3 : $\omega=0.018229$, $\varphi=0.912309$ et $\theta=-0.797182$
model_11.params
On obtien alors $\alpha = \varphi+\theta$, $\beta = - \theta$ et et on vérifie que ces paramètres satisfont les contraintes.
omega = model_11.params[0]
alpha = model_11.params[1] + model_11.params[2]
beta = - model_11.params[2]
omega>0, alpha>0, beta>0, alpha+beta<1
En prenant la condition iniale $\hat\sigma_0:=0$ (c'est arbitraire), on obtient un estimateur $\hat\sigma_n$ de la volatilité $\sigma_n$ en calculant récursivement pour $n\geq 0$, $$ \hat\sigma_{n+1}^2 := \omega+ \alpha \,x[n]^2 + \beta \,\hat\sigma_{n}^2 $$ et en prenant la racine carrée.
def Est_volatilité(x,n):
est_sigma_2 = np.zeros(n+1)
for i in range(n):
est_sigma_2[i+1] = omega + alpha * x[i]**2 + beta * est_sigma_2[i]
return(np.sqrt(est_sigma_2))
Dessinons un graphe de la volatilité ainsi estimée superposé à INTC.
N = len(INTC)
Vol = Est_volatilité(INTC,N)
plt.plot(Vol,linewidth=2)
plt.plot(INTC,linewidth=0.5)
plt.show()
- On observe que les pics de volatilités se superposent aux endroits où INTC fluctue le plus.
- On observe également des "clusters" de volatilité (les pics raprochés).
- On pourrait changer la condition initiale $\hat\sigma_0=0$ en une autre sans trop changer le modèle à moyen terme (essayez !).
REMARQUE : Concernant l'hypothèse de modélisation $\mathbb E[X_t^4]<\infty$ que nous avons faite ici et là, il n'y a pas de façon propre de tester cette hypothèse. On peut néanmoins se convaincre qu'elle est raisonnable ici car, si $X_t$ est un GARCH(1,1), en particulier elle est stationnaire et on peut prendre $\frac1N\sum_{t=1}^N X_t^4$ comme estimateur de $\mathbb E[X_t^4]$. On vérifie ensuite que ce cette quantité n'est pas "trop grande" ; pour l'ordre de grandeur, on peut se demander si $\frac1N\sum_{t=1}^N X_t^4$ n'est pas beaucoup plus grand que $\frac1N\sum_{t=1}^N X_t^2\simeq \text{Var}[X_t].$ Attention : ce n'est en rien une preuve, c'est juste un critère qualitatif.
np.mean(INTC**4), np.mean(INTC**2)
Complément : Testez ce critère dans l'Exercice 1, que l'on essaiera sur plusieurs réalisations d'un ARCH(1) suivant si $\alpha<1/\sqrt{3}$ ou non.