[PYTHON] Analyse de séries chronologiques Partie 1 Autocorrélation

Objectif

Données de séries chronologiques, RNN pour le moment! , LSTM! !! C'était comme ça, mais j'ai découvert ARIMA et SARIMA et je voulais les ajouter à mes tiroirs.
Il serait facile d'utiliser la bibliothèque et de s'y plonger, mais j'ai acheté un livre parce que je voulais connaître la théorie.
Un livre intitulé "Measurement Time Series Analysis of Economic and Finance Data" par le professeur Tatsuyoshi Okimoto. Ceci est un livre de 2010.

Les données

J'ai décidé d'utiliser les données historiques TOPIX pour le moment car je voulais bouger mes mains en lisant. https://quotes.wsj.com/index/JP/TOKYO%20EXCHANGE%20(TOPIX)/180460/historical-prices J'ai pu le télécharger à partir du site Web du WSJ.

Il semble qu'il n'y ait aucune valeur manquante pour le moment. La période est du 30 décembre 2008 au 29 novembre 2019. Données historiques de 4 valeurs sans ajustement ex-dividende. Lorsque vous tracez le cours de clôture, il ressemble à la figure ci-dessous. C'est une période entre le moment où les 1000 points ont été cassés immédiatement après le choc de Lehman et le moment où Avenomics a été adopté. Les rendements quotidiens sont les suivants. Comme c'est quotidien, cela ne fera pas une grande différence, mais j'ai utilisé un retour logarithmique ($ \ Delta \ log {y_t} = \ log {(y_t)} - \ log {(y_ {t-1})} $). La plus forte baisse a été enregistrée le 15 mars 2011, à $ \ Delta \ log {y_t} = -0,0995 $. La cause était que le tremblement de terre au large de l'océan Pacifique dans la région de Tohoku s'est produit environ 15 minutes avant la fermeture de la veille.

Diverses statistiques

moyenne \bar{y} = \frac{1}{T}\displaystyle{\sum_{t=1}^{T}y_t}

tpx_return = np.log(tpx['close'].values)[1:] - np.log(tpx['close'].values)[:-1]
tpx_return.mean()

0.00025531962222667643

** auto-covariance ** \hat{\gamma}_k = \frac{1}{T} \displaystyle{\sum\_{t=k+1}^{T}}(y\_t-\bar{y})(y\_{t-k}-\bar{y}),\quad k = 0,1,2,...

import statsmodels.api as sm
sm.tsa.stattools.acovf(tpx_return, fft=True, nlag=5)

array([ 1.57113176e-04,  1.16917913e-06,  3.48846296e-06, -4.67502133e-06, -5.31500282e-06, -2.82855375e-06])

Je ne suis pas habitué à la bibliothèque statsmodels, alors vérifiez-la à la main.

# k=0
((tpx_return-tpx_return.mean())**2).sum() / len(tpx_return)

0.00015711317609153836

# k=4
((tpx_return-tpx_return.mean())[4:]*(tpx_return-tpx_return.mean())[:-4]).sum() / len(tpx_return)

-5.315002816332674e-06

** Auto-corrélation ** \hat{\rho}_k = \frac{\hat{\gamma}\_k}{\hat{\gamma}\_0},\quad k=1,2,3,...

sm.tsa.stattools.acf(tpx_return)[:5]

array([ 1.        ,  0.00744164,  0.0222035 , -0.02975576, -0.03382913])

Si vous vérifiez en utilisant les résultats de $ k = 0 $ et $ k = 4 $ plus tôt,

-5.315002816332674e-06 / 0.00015711317609153836

-0.03382913482212345

Il semble donc que la bibliothèque fasse le calcul comme je l'imaginais.

Je vais également dessiner un cholérogramme. Un cholérogramme est un graphique du coefficient d'autocorrélation.

autocorr = sm.tsa.stattools.acf(tpx_return)
ci = 1.96 / np.sqrt(len(tpx_return))
plt.bar(np.arange(len(autocorr)), autocorr)
plt.hlines([ci,-ci],0,len(autocorr), linestyle='dashed')
plt.title('Correlogram')
plt.ylim(-0.05,0.05)
plt.show()

L'IC (intervalle de confiance) est un intervalle de confiance. Lorsque les données sont indépendantes les unes des autres et suivent la même distribution, $ \ hat {\ rho} _k $ fait progressivement la moyenne de $ 0 $ et les variances $ \ frac {1} {T. Nous calculons 95% des deux côtés en utilisant la propriété de suivre la distribution normale de} $.

Cela ressemble à ceci lors de l'utilisation de la bibliothèque. À la mode.

sm.graphics.tsa.plot_acf(tpx_return)
plt.ylim(-0.05,0.05)
plt.show()

Pensez un peu aux chiffres que vous avez. \hat{\rho}\_{k=11}=-0.0421, \quad \hat{\rho}_{k=16}=0.0415 Il est légèrement supérieur à $ CI = 0,0379 $ deux fois. En particulier, $ \ hat {\ rho} \ _ {k = 12} = -0.0323 $ est également une valeur négative, ce qui suggère que même si le prix du marché commence à augmenter ou baisser une fois, l'élan prendra fin dans environ 2 semaines. Vous pourrez peut-être penser que vous le faites.

Test du sac (test du portemanteau)

Une méthode pour tester l'hypothèse nulle selon laquelle plusieurs coefficients d'autocorrélation sont tous égaux à 0. H_0 : \rho_1 = \rho_2 =\quad ... \quad= \rho_m = 0 Ici, nous allons effectuer le test en utilisant les statistiques de Ljung et Box (1978) présentées dans le livre. L'approche consiste à comparer le $ Q (m) $ ci-dessous avec le point à 95% de la distribution du chi carré. Q(m) = T(T+2)\displaystyle{\sum\_{k=1}^{m}}\frac{\hat{\rho}^2_k}{T-k} \sim \chi^2(m) De plus, une statistique appelée valeur P est également définie, qui indique la probabilité qu'une variable de probabilité suivant la distribution du chi carré prenne une valeur supérieure à $ Q (m) $. Autrement dit, en supposant un niveau de signification de 5%, $ H_0 $ est rejeté lorsque la valeur P est inférieure à 0,05. Quant à la valeur de $ m $, $ m \ approx \ log {(T)} $ semble être un guide, mais il semble qu'il soit courant de tester plusieurs $ m $ et de porter un jugement complet. Est.

Tout d'abord, à partir du modèle qui utilise la bibliothèque de modèles de statistiques. $ m \ approx \ log {(T)} = 7.89 $ Donc, pour le moment, considérez le décalage dans la plage allant jusqu'à 16.

lvalue, pvalue = sm.stats.diagnostic.acorr_ljungbox(tpx_return)

C'est tout. Très facile.

La valeur P n'est pas tombée en dessous de 0,05 quelle que soit la valeur de $ m $, et le résultat est que le taux de changement quotidien de TOPIX ne peut être considéré comme ayant une autocorrélation. Eh bien, le prix du marché n'est pas si simple.

Enfin, essayez sans bibliothèque pour approfondir votre compréhension.

from scipy.stats import chi2

def Q_func(data, max_m):
    T = len(data)
    auto_corr = sm.tsa.stattools.acf(data)[1:]
    lvalue = T*(T+2)*((auto_corr**2/(T-np.arange(1,len(auto_corr)+1)))[:max_m]).cumsum()
    pvalue = 1 - chi2.cdf(lvalue, np.arange(1,len(lvalue)+1))
    return lvalue, pvalue

Confirmez que le même résultat est obtenu.

l_Q_func, p_Q_func = Q_func(tpx_return,max_m=16)
l_sm, p_sm = sm.stats.diagnostic.acorr_ljungbox(tpx_return, lags=16)
((l_Q_func-l_sm)**2).mean(), ((p_Q_func-p_sm)**2).mean()

(0.0, 7.824090399073146e-34)

Si vous regardez de plus près le cas de $ m = 8 $,

T = len(tpx_return)
auto_corr = sm.tsa.stattools.acf(tpx_return)[1:]
lvalue = T*(T+2)*((auto_corr**2/(T-np.arange(1,len(auto_corr)+1)))[:8]).sum()
print(lvalue)

8.604732778577853

1-chi2.cdf(lvalue,8)

0.37672860496603844

Il s'avère donc qu'il est assez difficile de rejeter l'hypothèse nulle.