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, acovf # acovf = auto-covariance
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.stats.diagnostic import acorr_ljungbox
ar_parms = [1,-1.5,0.75]
ma_parms = [1]
AR2 = ArmaProcess(ar=ar_parms,ma=ma_parms)
Xt = AR2.generate_sample(500)
plt.plot(Xt)
plt.show()
AR2.arroots # on voit que les racines du polynôme de la partie AR sont de module > 1.
AR2.isstationary
AR2.isinvertible
# ACF
plot_acf(Xt,lags=40)
plt.show()
# PACF NB : Python prend pour convention PACF(0)=1
plot_pacf(Xt,lags = 40)
plt.show()
Question (4)¶
On s'intéresse maintenant au modèle ARMA(2,1) de paramètres $\phi_1=0.89$, $\phi_2=-0.5$ et $\theta_1=-0.23$.
ar_parms2 = [1,-0.89,0.5]
ma_parms2 = [1,-0.23]
ARMA21 = ArmaProcess(ar=ar_parms2,ma=ma_parms2)
Yt = ARMA21.generate_sample(1000)
plt.plot(Yt)
plt.show()
plot_acf(Yt,lags=40)
plot_pacf(Yt,lags=40)
plt.show()
# Retrouvons les paramètres de ce modèle
model = ARMA(Yt,order=(2,1)).fit(trend='nc')
model.summary()
residus21 = model.resid
acorr_ljungbox(residus21,lags=range(100),boxpierce=True)[3] < 0.05
Les résidus passent le test de Box-Pierce à $\alpha=5\%$ pour des lags $h$ allant de $1$ à $60$ : pas d'évidence contre l'hypothèse que les résidus forment un bruit blanc.
On cherche maintenant à retrouver l'ordre de $(p,q)$ de l'ARMA en regardant quels les ordres qui ont un critère BIC et/ou AIC minimal pour, mettons, $1\leq p,q\leq 4$.
- NB : Si on considère une classe de modèles paramétrés (comme ARMA(p,q)) de paramètres $\theta$ (ici $\phi_1,\ldots,\phi_p,\theta_1,\ldots,\theta_q$) et qu'on note $k$ le nombre de paramètres (ici $k=p+q$), on définit le critère d'information d'Akaike par :
et le critère d'information Bayésien par
$$ \text{BIC} := \log(N)k-2\max_{\theta}\log\text{Vraisemblance}(\theta|\text{données}) $$où $N$ est le nombre d'observations (données).
- L'AIC pénalise donc le nombre de paramètres moins fortement que le BIC dès que $N\geq 8$.
l = 4
AIC = np.zeros((l,l))
BIC = np.zeros((l,l))
for p in range(l):
for q in range(l):
try:
model = ARMA(Yt,order=(p,q)).fit(trend='nc')
AIC[p,q] = model.aic
BIC[p,q] = model.bic
except: continue
# l'utilisation de "try"/"continue" permet de gérer les situations où ARMA.fit
# n'arrive pas à trouver des paramètres adaptés pour un ARMA de paramètres p,q donnés
# (par exemple si le modèle n'est pas stationnaire)
AIC - 2900 # on soustrait 2900 pour lire plus facilement les résultats.
On trouve que le modèle ARMA(2,1) a un AIC minimal, suivit par un AR(3) et un AR(2)
BIC - 2900
Ici le critère BIC est minimal pour AR(2), puis ARMA(2,1) et AR(3)
On retrouve donc le vrai modèle ARMA(2,1) dans le top 3 de ces deux critères. Il est aussi tentant de tester un AR(2) et AR(3) :
model_AR2 = ARMA(Yt,order=(2,0)).fit(trend='nc')
residus2 = model_AR2.resid
acorr_ljungbox(residus2,lags=range(100),boxpierce=True)[3] < 0.05
model_AR3 = ARMA(Yt,order=(3,0)).fit(trend='nc')
residus3 = model_AR3.resid
acorr_ljungbox(residus3,lags=range(100),boxpierce=True)[3] < 0.05
On voit alors que modéliser notre réalisation d'un ARMA(2,1) par un AR(2) ou un AR(3) est tout aussi raisonnable sur cet exemple : en général il n'y a de "bonne réponse" unique en modélisation.
EXERCICE 2¶
Oil_price = read_csv('oil.csv')
type(Oil_price)
Oil_price.head()
Oil_price.plot()
plt.show()
log_rend_Oil = np.log(Oil_price).diff()
log_rend_Oil.head()
log_rend_Oil = log_rend_Oil[1:545] # on supprime le Nan
log_rend_Oil.plot()
plt.show()
plot_acf(log_rend_Oil,lags=10)
plot_pacf(log_rend_Oil,lags=10)
plt.show()
L'ACF et la PACF ne semblent pas indiquer quelconque préférence vers un modèle AR ou MA. On essaie quand même une modélisation avec un AR pour voir ce qu'il se passe.
Modélisation AR¶
l = 5
# MODELISATION AR
AIC_AR = np.zeros(l)
BIC_AR = np.zeros(l)
for p in range(l):
model = ARMA(log_rend_Oil,order=(p+1,0)).fit(trend='nc')
AIC_AR[p] = model.aic
BIC_AR[p] = model.bic
AIC_AR + 1800, BIC_AR + 1800
A la vue des deux criètes, le modèle AR(3) semble le plus adapté.
model_AR3 = ARMA(log_rend_Oil,order=(3,0)).fit(trend='nc')
model_AR3.summary()
On obtient le modèle : $$ X_t = 0.1582 X_{t-1}- 0.1069 X_{t-2} +0.1621 X_{t-3}+\epsilon_t. $$ Celà conduit à la prédiction (en négligeant le bruit) :
X_pred = 0.1582*log_rend_Oil.x[544] - 0.1069*log_rend_Oil.x[543] + 0.1621*log_rend_Oil.x[542]
X_pred
Mais peut-on négliger le bruit ? L'erreur est de l'ordre de grandeur de sa moyenne plus ou moins sa déviation standard. Il faudrait aussi montrer que $\epsilon_t$ forme un bruit blanc (décorélé du temps) : on fera un test de blancheur par la suite.
# résidus de la modélisation
Err_AR3 = model_AR3.resid
# moyenne et déviation standard
np.mean(Err_AR3), np.std(Err_AR3)
On peut aussi essayer de prédire une valeur passée, pour vérifier. Par exemple la dernières valeur des log-rendements est log_rend_Oil.x[544] que l'on va prédire par
X544_pred = 0.1582*log_rend_Oil.x[543] - 0.1069*log_rend_Oil.x[542] + 0.1621*log_rend_Oil.x[541]
X544_pred, log_rend_Oil.x[544] # pas terrible vu les ordres de grandeur...
# Etude graphique des résidus
Err_AR3.plot()
plot_acf(Err_AR3)
plot_pacf(Err_AR3)
plt.show()
# Test de blancheur ?
acorr_ljungbox(Err_AR3,lags=range(60),boxpierce=True)[3] < 0.05
# L'hypothèse que les résidus forment un bruit blanc est refusée (à 5%) !
La modélisation AR permet une prédiction interessante comme indiqué ci-dessus car on sait que le prédicteur ci-dessus est également le prédicteur linéaire optimal. Malheureusement l'étude des résidus de la modélisation ne permet pas de supposer que l'erreur de modélisation forme un bruit blanc : la modélisation AR et la prédiction associée ne sont pas adaptées.
Modélisation ARMA¶
# le cas MA
l = 5
AIC_MA = np.zeros(l)
BIC_MA = np.zeros(l)
for q in range(l):
model = ARMA(log_rend_Oil,order=(0,q+1)).fit(trend='nc')
AIC_MA[q] = model.aic
BIC_MA[q] = model.bic
# le cas ARMA non AR ni MA
AIC_ARMA = np.zeros((l,l))
BIC_ARMA = np.zeros((l,l))
for p in range(l):
for q in range(l):
try:
model = ARMA(log_rend_Oil,order=(p+1,q+1)).fit(trend='nc')
AIC_ARMA[p,q] = model.aic
BIC_ARMA[p,q] = model.bic
except: continue
AIC_AR + 1900, AIC_MA + 1900 , AIC_ARMA + 1900
Les modèles minimisant AIC sont donc (classés par AIC croissant) : ARMA(4,4), ARMA(5,3), MA(3), ARMA(3,3), ...
NB : On pourrait automatiser cette recherche de valeur minimale (par exemple en utilisant la fonction argmin de numpy)
BIC_AR + 1800, BIC_MA + 1800, BIC_ARMA + 1800
Les modèles minimisant BIC sont donc (classés par BIC croissant) : ARMA(1,1), MA(3), AR(3), ...
En testant tous ces quelques modèles, on trouve que le seul qui passe le test de blancheur pour les résidus est ARMA(4,4).
model_ARMA44 = ARMA(log_rend_Oil,order=(4,4)).fit(trend='nc')
model_ARMA44.summary()
# Test de blancheur :
Err_ARMA44 = model_ARMA44.resid
acorr_ljungbox(Err_ARMA44,lags=range(100),boxpierce=True)[3] < 0.05
# L'hypothèse que les résidus forment un bruit blanc est n'est pas réfusée (à 5%)
Le modèle ARMA(4,4) est donc une modélisation raisonnable des log-rendements. Par contre, pour la prédiction, il n'est pas clair comment procéder avec ce qui a été vu en cours pour l'instant...
Exercice 3¶
On implémente l'algorithme de Durbin-Levinson comme vu dans le cours :
$$ \varphi_{11}= \frac{\gamma(1)}{\gamma(0)},\qquad \sigma_1^2=\gamma(0)-\frac{\gamma(1)^2}{\gamma(0)}=\gamma(0)(1-\varphi_{11}^2) $$puis, pour tout $m\geq 1$ et si on pose $$ \kappa_m:= \frac{\gamma(m+1)-\sum_{k=1}^m\varphi_{k,m}\gamma(m-k+1)}{\sigma_m^2}, $$ alors on a la relation de récurrence $$ \begin{pmatrix}\varphi_{1,m+1}\\\vdots\\\varphi_{m,m+1}\end{pmatrix}=\begin{pmatrix}\varphi_{1,m}\\\vdots\\\varphi_{m,m}\end{pmatrix} -\kappa_m \begin{pmatrix}\varphi_{m,m}\\\vdots\\\varphi_{1,m}\end{pmatrix},\qquad \varphi_{m+1,m+1}=\kappa_m $$ ainsi que $$ \sigma_{m+1}^2=\sigma_m^2(1-\kappa_m^2) $$
def Durbin_Levinson(serie,M):
# retourne phi = [[phi_11],[phi_12,phi_22],...,[phi_1M,...,phi_MM]]
# et sigma2 = [sigma_1^2,...,sigma_M^2]
# initialisation de l'algorithme
gamma = acovf(serie,fft=False)
phi = []
sigma2 = []
for m in range(M):
if m==0:
phi11 = gamma[1]/gamma[0]
phi.append( [phi11] )
sigma2.append( gamma[0]*(1-phi11**2) )
else:
phi_m = phi[m-1] # [phi_1m,...,phi_mm]
gamma_m_inv = np.flip(gamma[1:m+1]) # [gamma(m),...,gamma(1)]
kappa = (gamma[m+1] - phi_m @ gamma_m_inv)/sigma2[m-1]**2
phi_mplus1 = list(phi_m - kappa * np.flip(phi_m))
phi_mplus1.append(kappa)
phi.append(phi_mplus1)
sigma2.append(sigma2[m-1]*(1-kappa**2))
return(phi,sigma2)
Vérifions l'algorithme sur un modèle AR(p),
$$ X_t = \varphi_1 X_{t-1}+\dots+\varphi_pX_{t-p}+\epsilon_t, $$où l'on sait que pour tout $m\geq p$ on a : $$ \varphi_{k,m}=\begin{cases} \varphi_k & \text{ si } k\leq p\\ 0 & \text{ sinon } \end{cases} $$
# Modèle AR(2)
ar_parms_testDL = [1,-0.9,0.4]
ma_parms_testDL = [1]
testDL = ArmaProcess(ar=ar_parms_testDL,ma=ma_parms_testDL).generate_sample(2000)
Durbin_Levinson(testDL,4)[0][1]
Durbin_Levinson(testDL,4)[0][2]
Durbin_Levinson(testDL,4)[0][3]