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

1. Vue d'ensemble

2. Données

En continuant de la fois précédente, nous utiliserons les données historiques TOPIX. save.png 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/ save.png À 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.

  1. 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. save.png 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. save.png

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.

  1. 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()

save.png

\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. save.png

  1. 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)

save.png

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

À 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.

modèle Autocorrélation 偏Autocorrélation
AR(p) Pourriture p+1Après le prochain0
MA(q) q+1Après le prochain0 Pourriture
ARMA(p,q) Pourriture Pourriture

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)
month visits
0 1996-01-01 276086
1 1996-02-01 283667
2 1996-03-01 310702
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()

save.png 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()

save.png

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.

0 1 2
0 3368.221521 3363.291271 3350.327397
1 3365.779849 3348.257023 3331.293926
2 3365.724955 3349.083663 3328.831252
3 3361.660853 3347.156390 3329.447773
4 3332.100417 3290.984260 3292.800604

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()

save.png

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

Recommended Posts

Analyse des séries chronologiques Partie 2 AR / MA / ARMA
Analyse des séries chronologiques partie 4 VAR
Analyse de séries chronologiques Partie 3 Prévisions
Analyse de séries chronologiques Partie 1 Autocorrélation
Analyse des séries chronologiques 2 Stabilité, modèle ARMA / ARIMA
Python: analyse des séries chronologiques
Analyse des séries chronologiques RNN_LSTM1
Analyse des séries chronologiques 1 Principes de base
Python: Analyse des séries temporelles: Constantity, modèle ARMA / ARIMA
Python: analyse des séries chronologiques: prétraitement des données des séries chronologiques
Analyse des séries chronologiques
Analyse des séries chronologiques 3 Prétraitement des données des séries chronologiques
[Statistiques] [Analyse des séries chronologiques] Tracez le modèle ARMA et saisissez la tendance.
Analyse des séries chronologiques 4 Construction du modèle SARIMA
Analyse des séries chronologiques n ° 6 Faux retour et partie républicaine
Python: analyse des séries temporelles: création d'un modèle SARIMA
Décomposition des séries temporelles
série pandas partie 1
Python 3.4 Créer un environnement Windows7-64bit (pour l'analyse des séries chronologiques financières)
Kaggle ~ Analyse du logement ③ ~ Part1
Question sur la série chronologique Python
Afficher les séries chronologiques TOPIX
Diagramme de séries chronologiques / Matplotlib
Défi des prévisions de ventes futures: ② Analyse des séries chronologiques à l'aide de PyFlux
Une méthode d'étude pour les débutants pour apprendre l'analyse des séries chronologiques
Défi des prévisions de ventes futures: ⑤ Analyse des séries chronologiques par Prophet
Défi pour les prévisions de ventes futures: ① Qu'est-ce que l'analyse des séries chronologiques?