[PYTHON] Analyse des séries chronologiques Partie 2 AR / MA / ARMA

1. Vue d'ensemble

• Dans la continuité de la partie 1 (https://qiita.com/asys/items/d28a0ce0a8f51681a243), j'étudie sur la base de «l'analyse des séries chronologiques de mesure des données économiques et financières».
• Cette fois sur le modèle ARMA, qui correspond au chapitre 2.

2. Données

En continuant de la fois précédente, nous utiliserons les données historiques TOPIX. De plus, j'ai senti qu'il y avait des données qui semblaient avoir une autocorrélation, j'ai donc préparé des données sur le nombre de visiteurs étrangers au Japon chaque mois. J'ai utilisé celui sur le site Web de JTB Research Institute lié ci-dessous. https://www.tourism.jp/tourism-database/stats/inbound/ À première vue, les données sont saisonnières et le yen se déprécie dans les Abenomics, donc la tendance générale est à la hausse.

3. Moving Average Process Commençons par le premier processus MA. MA(1):y_t=\mu+\epsilon_t+\theta_1\epsilon_{t-1}, \quad \epsilon_t \sim W.N.(\sigma^2) Décidons quelques paramètres et tracent-les. Pour le moment, définissez $ \ sigma = 1.0 $ et utilisez random.randn de numpy.

wn = np.random.randn(50)
mu = 0.0
theta_1 = 1.0
map = mu + wn[1:] + theta_1 * wn[:-1]

Les données ont été générées comme ça. En considérant la ligne verte $ \ theta_1 = 0.0 $ comme référence, il est facile de comprendre l'influence du signe et de la valeur de $ \ theta_1 . La ligne bleue ( \ theta_1 = 1.0 ) et la ligne rouge ( \ theta_1 = -1.0 ) sont situées sur les côtés opposés de la ligne verte, et la ligne orange ( \ theta_1 = 0,5 $) est bleue. Il se déplace entre la ligne et la ligne verte.

Les principales statistiques sont les suivantes. mean : \quad E(y_t)=\mu variance : \quad \gamma_0=Var(y_t)=(1+\theta_1^2)\sigma^2 first-order auto-covariance : \quad \gamma_1=Cov(y_t,y_{t-1}) = \theta_1\sigma^2 first-order auto-correlation : \quad \rho_1=\frac{\gamma_1}{\gamma_0}=\frac{\theta_1}{1+\theta_1^2}

Vérifiez le cholérogramme. Dans le graphique ci-dessus, le nombre de données est de 50, mais ici il est créé avec 1000 données.

Jusqu'à présent, c'est une image intuitive. Seule l'autocorrélation du premier ordre est significative, et lorsque $ \ theta_1 <0 $, l'autocorrélation est une valeur négative.

Ensuite, regardons une généralisation du processus MA. MA(q):y_t=\mu+\epsilon_t+\theta_1\epsilon_{t-1}+\theta_2\epsilon_{t-2}+\cdots+\theta_q\epsilon_{t-q}, \quad \epsilon_t \sim W.N.(\sigma^2)

Les principales statistiques sont les suivantes. \quad E(y_t)=\mu \quad \gamma_0=Var(y_t)=(1+\theta_1^2+\theta_2^2+\cdots+\theta_q^2)\sigma^2 \quad\gamma_k=\left\\{\begin{array}{ll}(\theta_k+\theta_1\theta_{k+1}+\cdots+\theta_{q-k}\theta_q)\sigma^2,& 1\leqq k\leqq q\\\0, & k\geqq q+1\end{array}\right. Le processus $ \ quad $ MA est toujours stable \quad\rho_k=\left\\{\begin{array}{ll}\frac{\theta_k+\theta_1\theta{k+1}+\cdots+\theta_{q-k}\theta_q}{1+\theta_1^2+\theta_2^2+\cdots+\theta_q^2},&1\leqq k\leqq q\\\0,&k\geqq q+1\end{array}\right.

[Supplément] Qu'est-ce que la stationnarité? Si moins de $ \ quad $ tient, le processus est une stationnarité faible. \qquad \forall t,k\in\mathbb{N},\quad E(y_t)=\mu \land Cov(y_t,t_{t-k})=E[(y_t-\mu)(y_{t-k}-\mu)]=\gamma_k $ \ quad $ Aussi, à ce moment, l'autocorrélation \qquad Corr(y_t,y_{t-k})=\frac{\gamma_kt}{\sqrt{\gamma_{0,t}\gamma_{0,t-k}}}=\frac{\gamma_k}{\gamma_0}=\rho_k $ \ Quad $ tient. $ \ quad $ L'auto-covariance et l'autocorrélation dépendent uniquement du décalage horaire $ k $. $ \ Quad $ De plus, $ \ gamma_k = \ gamma_ {-k} $ et $ \ rho_k = \ rho_ {-k} $ sont établis.

4. Autoregressive Process Vient ensuite le processus AR. Tout d'abord, à partir du processus AR primaire. AR(1):y_t=c+\phi_1y_{t-1}+\epsilon_t,\quad \epsilon_t\sim W.N.(\sigma^2)

Tracons-le. Ici, $ y_0 = 0 $, $ c = 1 $ et $ \ sigma ^ 2 = 1 $.

AR1 = []
y = 0
c = 1
phi_1 = 0.5
for i in range(200):
    AR1.append(y)
    y = c + phi_1*y + np.random.randn()

\phi_1=1.0À ce moment-là, on peut comprendre visuellement que le processus n'est pas régulier. En réalité,|\phi_1|\leq1À ce moment, le processus devient stable. Quant à la valeur de $ \ phi_1 $, on constate que plus la valeur est petite, plus la valeur moyenne est basse, et plus la valeur est petite, plus la valeur moyenne est élevée ou faible.

Les principales statistiques sont les suivantes. cependant,|\phi_1|\leq1alorsyEst stable. \quad E(y_t)=c+\phi_1E(y_{t-1})=\frac{c}{1-\phi_1} \quad Var(y_t)=\phi_1^2Var(y_{t-1})+\sigma^2=\frac{\sigma^2}{1-\phi_1^2} \quad \gamma_k=\phi_1\gamma_{k-1} $\quad \rho_k=\phi_1\rho_{k-1}\qquad\cdots $(Yule-Walker Equation)

Vient ensuite le processus de RA généralisé. AR(p):y_t=c+\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\epsilon_t, \quad \epsilon_t\sim W.N.(\sigma^2)

Les principales statistiques sont les suivantes. \quad \mu=E(y_t)=\frac{c}{1-\phi_1-\phi_2-\cdots-\phi_p} \quad \gamma_0=Var(y_t)=\frac{\sigma^2}{1-\phi_1\rho_1-\phi_2\rho_2-\cdots-\phi_p\rho_{p}} \quad \gamma_k=\phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}+\cdots+\phi_p\gamma_{k-p},\quad k\geqq1 \quad \rho_k=\phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p},\quad k\geqq1,\quad\cdots(Yule-Walker Equation) L'autocorrélation du processus $ \ quad $ AR décroît de manière exponentielle.

Le cholérogramme est le suivant.

5. Autoregressive Moving Average Process Les termes généraux du processus ARMA, qui combine le processus AR et le processus MA, sont donnés ci-dessous. ARMA(p,q):y_t=c+\phi_1y_{t-1}+\cdots+\phi_py_{t-p}+\epsilon_t+\theta_1\epsilon_t-1+\cdots+\theta_q\epsilon_{t-q},\quad \epsilon_t\sim W.N.(\sigma^2)

$ ARMA (p, q) $ a les propriétés suivantes. \quad \mu=E(y_t)=\frac{c}{1-\phi_1-\phi_2-\cdots-\phi_q} \quad \gamma_k=\phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}+\cdots+\phi_p\gamma_{k-p},\quad k\geqq q+1 \quad \rho_k=\phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\cdots+\phi_p\rho_{k-p},\quad k\geqq q+1 L'autocorrélation du processus $ \ quad $ ARMA décroît de façon exponentielle.

Tracons en fait $ ARMA (3,3) $. Les paramètres sont les suivants. c=1.5, \phi_1=-0.5, \phi_2=0.5,\phi_3=0.8, \theta_1=0.8, \theta_2=0.5, \theta_3=-0.5, \sigma=1

arma = [0,0,0]
wn = [0,0,0]
c = 1.5
phi_1 = -0.5
phi_2 = 0.5
phi_3 = 0.8
theta_1 = 0.8
theta_2 = 0.5
theta_3 = -0.5
for i in range(1000):
    eps = np.random.randn()
    y = c + phi_1*arma[-1] + phi_2*arma[-2] + phi_3*arma[-3] + eps + theta_1*wn[-1] + theta_2*wn[-2] + theta_3*wn[-3]
    wn.append(eps)
    arma.append(y)

• Stabilité et réversibilité du processus ARMA
  Le processus MA est stationnaire, mais le processus AR n'est pas toujours stationnaire.
  Puisqu'il existe plusieurs processus MA avec la même valeur attendue et la même structure d'autocorrélation, cela devient un problème dans la sélection du modèle.
  Ce qui est important en considérant cela, c'est la stabilité et la réversibilité.

La constance (C'est la même chose que l'état stable du processus AR.) AR characteristic equation : 1-\phi_1z-\cdots-\phi_pz^p=0 Le processus AR est stationnaire lorsque les valeurs absolues de toutes les solutions de cette équation caractéristique AR sont supérieures à 1.

Réversibilité Le processus MA est réversible lorsque le processus MA peut être réécrit en un processus AR (∞). MA characteristic equation : 1+\theta_1z+\theta_2z^2+\cdots+\theta_pz^p=0 Le processus MA est inversible lorsque toutes les solutions de cette équation caractéristique MA sont supérieures à 1. Des processus MA inversibles sont souhaitables pour la sélection de modèles.

6. Estimation du modèle ARMA

• Carrés minimum (OLS: moindres carrés ordinaires)
  OLS sélectionne les paramètres du modèle ARMA pour minimiser la somme des carrés des résidus (SSR) entre les valeurs réellement observées et les valeurs estimées par le modèle ARMA.
  \quad\displaystyle{\min_{\Theta}}\quad SSR=(y_t-ARMA(\Theta))^2
  Vous pouvez résoudre l'équation normale qui est partiellement différenciée par le paramètre de chaque paramètre.
  \quad \frac{\partial SSR}{\partial \Theta}=0

• Estimation du maximum de vraisemblance (MLE)
  La méthode la plus probable est l'approche consistant à sélectionner le paramètre qui maximise la vraisemblance logarithmique.
  \quad\displaystyle{\max_{\Theta}}\quad \mathcal{L}(\theta) = \displaystyle\sum_{t=1}^{T}\log f_{Y_t|Y_{t-1}}(y_t|y_{t-1};\Theta)

À partir de là, estimons les paramètres pour $ ARMA (3,3) $ réellement créés en 5. Comme il est compliqué de calculer à la main, utilisez la bibliothèque statsmodels. Pour plus de détails, accédez à l'URL suivante. https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima_model.ARMA.html Il semble que si vous donnez l'échantillon échantillon et l'ordre du modèle avec order = (p, q) à l'argument d'ARMA, le paramètre sera calculé par l'estimation la plus probable.

arma_model = sm.tsa.ARMA(arma, order=(3,3))
result = arma_model.fit()
print(result.summary())

                              ARMA Model Results                              
==============================================================================
Dep. Variable:                      y   No. Observations:                 1003
Model:                     ARMA(3, 3)   Log Likelihood               -1533.061
Method:                       css-mle   S.D. of innovations              1.113
Date:                Sun, 08 Dec 2019   AIC                           3082.123
Time:                        14:51:34   BIC                           3121.409
Sample:                             0   HQIC                          3097.052
                                                                              
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          6.8732      0.410     16.773      0.000       6.070       7.676
ar.L1.y       -0.4986      0.024    -21.149      0.000      -0.545      -0.452
ar.L2.y        0.5301      0.019     27.733      0.000       0.493       0.568
ar.L3.y        0.8256      0.020     41.115      0.000       0.786       0.865
ma.L1.y        0.7096      0.036     19.967      0.000       0.640       0.779
ma.L2.y        0.3955      0.039     10.176      0.000       0.319       0.472
ma.L3.y       -0.4078      0.033    -12.267      0.000      -0.473      -0.343
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.0450           -0.0000j            1.0450           -0.0000
AR.2           -0.8436           -0.6689j            1.0766           -0.3933
AR.3           -0.8436           +0.6689j            1.0766            0.3933
MA.1           -0.6338           -0.8332j            1.0469           -0.3535
MA.2           -0.6338           +0.8332j            1.0469            0.3535
MA.3            2.2373           -0.0000j            2.2373           -0.0000
-----------------------------------------------------------------------------

Etc. Il est possible d'estimer une valeur raisonnablement proche du paramètre d'origine. Quelques notes, S.D. des innovations représente la dispersion du bruit blanc, et pour trouver la valeur de c à partir de const,

print(result.params[0] * (1-result.arparams.sum()))

0.9824998883509097

Besoin de faire comme.

7. Sélection du modèle ARMA

Pour sélectionner $ ARMA (p, q) $, commencez par réduire de manière conservatrice les valeurs de $ p, q $ à partir des valeurs d'autocorrélation et d'autocorrélation partielle, puis déterminez le modèle optimal en fonction de la quantité d'informations. prendre.

Premièrement, l'autocorrélation partielle d'ordre k supprime l'influence de $ y_ {t-1}, \ cdots, y_ {t-k + 1} $ de $ y_t $ et $ y_ {tk} $. Il est défini comme la corrélation des choses.

\left\\{\begin{array}{lll}MA(q):&y_t=\mu+\epsilon_t+\theta_1\epsilon_{t-1}+\theta_2\epsilon_{t-2}+\cdots+\theta_q\epsilon_{t-q}\\\AR(p):&y_t=c+\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\epsilon_t\\\ARIMA(p,q):&y_t=c+\phi_1y_{t-1}+\cdots+\phi_py_{t-p}+\epsilon_t+\theta_1\epsilon_t-1+\cdots+\theta_q\epsilon_{t-q} \end{array}\right.

Compte tenu de la définition ci-dessus, les propriétés indiquées dans le tableau ci-dessous sont claires.

Cela permet d'affiner dans une certaine mesure le modèle en examinant la structure de l'autocorrélation de l'échantillon et de l'autocorrélation partielle de l'échantillon.

Vient ensuite le critère d'information (IC). \quad IC=-2\mathcal{L}(\hat{\Theta})+p(T)k \quad where T : number of samples, k : number of parameters Il existe les deux normes de quantité d'information couramment utilisées suivantes. $ \ left \ {\ begin {array} {ll} Akaike Information Amount Standard (AIC): \ quad & AIC = -2 \ mathcal (L) (\ hat {\ Theta}) + 2k \\ Schwarz Information Amount Standard (AIC) SIC): & SIC = -2 \ mathcal (L) (\ hat {\ Theta}) + \ log (T) k \ end {array} \ right. $ SIC a tendance à choisir un modèle plus petit, mais si le nombre d'échantillons est suffisant, l'estimation la plus probable du paramètre de la partie excessive de l'AIC peut converger vers 0, c'est donc une règle générale qui est la meilleure norme. Ce n'est pas possible.

8. Modélisation des données réelles

Nous avons préparé des données sur TOPIX et le nombre de visiteurs étrangers au Japon, mais les deux sont en train de ne pas être stables. Par conséquent, je n'avais pas d'autre choix que de créer des données qui excluent les tendances en soustrayant la moyenne mobile pour la partie 1996-2007 du nombre de visiteurs étrangers au Japon.

visit.head(3)

visit['trend'] = visit['visits'].rolling(window=13).mean().shift(-6)
visit['residual'] = visit['visits'] - visit['trend']
v = visit[visit['month']<'2008-01-01']
plt.plot(v['month'].values, v['visits'].values, label='original')
plt.plot(v['month'].values, v['trend'].values, label='trend')
plt.plot(v['month'].values, v['residual'].values, label='residual')
plt.legend()
plt.grid()
plt.title('Foreign Visitors')
plt.show()

Comme tendance, j'ai décidé de prendre une moyenne mobile de 6 mois avant et après. Ceci est traité en roulant (window = 13) .mean () et shift (-6).

Tout d'abord, examinons l'autocorrélation et l'autocorrélation partielle. Très facile avec la bibliothèque.

fig = plt.figure(figsize=(8,5))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(v['residual'].dropna().values, lags=40, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(v['residual'].dropna().values, lags=40, ax=ax2)
plt.tight_layout()
fig.show()

Il existe une forte saisonnalité qui ne peut être cachée. .. .. La corrélation tous les 12 mois est très forte. Vous pouvez utiliser arma_order_select_ic de statsmodels pour estimer l'ordre optimal de votre modèle ARMA. Pour plus de détails, accédez à ce lien (https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.arma_order_select_ic.html)

sm.tsa.arma_order_select_ic(v['residual'].dropna().values, max_ar=4, max_ma=2, ic='aic')

Ici, pour $ ARMA (p, q) $, la norme AIC est utilisée avec la valeur maximale de $ p $ étant de 4 et la valeur maximale de $ q $ étant de 2.

Je pense qu'aucun d'entre eux ne change beaucoup, mais le résultat est que $ p = 4, q = 1 $ est bon.

arma_model = sm.tsa.ARMA(v['residual'].dropna().values, order=(4,1))
result = arma_model.fit()
plt.figure(figsize=(10,4))
plt.plot(v['residual'].dropna().values, label='residual')
plt.plot(result.fittedvalues, label='ARMA(4,1)')
plt.legend()
plt.grid()
plt.title('ARMA(4,1)')
plt.show()

C'est comme ça. Personnellement, je sentais qu'il faisait de son mieux malgré le fait qu'il était un modèle relativement simple.