[PYTHON] Zeitreihenanalyse Teil 2 AR / MA / ARMA

1. Übersicht

• In Fortsetzung von Teil 1 (https://qiita.com/asys/items/d28a0ce0a8f51681a243) studiere ich basierend auf "Zeitreihenanalyse von Wirtschafts- und Finanzdaten".
• Diesmal über das ARMA-Modell, das Kapitel 2 entspricht.

2. Daten

In Anlehnung an die vorherige Zeit werden wir historische TOPIX-Daten verwenden. Außerdem hatte ich das Gefühl, dass es Daten gab, die eine Autokorrelation zu haben schienen, und so bereitete ich jeden Monat Daten über die Anzahl ausländischer Besucher in Japan vor. Ich habe die auf der Website des JTB Research Institute unten verlinkte verwendet. https://www.tourism.jp/tourism-database/stats/inbound/ Auf den ersten Blick sind die Daten saisonabhängig und der Yen wertet in Abenomics ab, sodass der Gesamttrend nach oben tendiert.

3. Moving Average Process Beginnen wir mit dem ersten MA-Prozess. MA(1):y_t=\mu+\epsilon_t+\theta_1\epsilon_{t-1}, \quad \epsilon_t \sim W.N.(\sigma^2) Lassen Sie uns einige Parameter festlegen und sie zeichnen. Setzen Sie vorerst $ \ sigma = 1.0 $ und verwenden Sie random.randn von numpy.

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

Die Daten wurden so generiert. Wenn man die grüne Linie $ \ theta_1 = 0.0 $ als Referenz betrachtet, ist der Einfluss des Vorzeichens und des Werts von $ \ theta_1 $ leicht zu verstehen. Die blaue Linie ($ \ theta_1 = 1,0 ) und die rote Linie ( \ theta_1 = -1,0 ) befinden sich auf gegenüberliegenden Seiten der grünen Linie, und die orange Linie ( \ theta_1 = 0,5 $) ist blau. Es bewegt sich zwischen der Linie und der grünen Linie.

Die Hauptstatistiken lauten wie folgt. 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}

Überprüfen Sie das Cholerogramm. In der obigen Darstellung beträgt die Anzahl der Daten 50, hier wird sie jedoch mit 1.000 Daten erstellt.

Bis zu diesem Punkt ist es ein intuitives Bild. Nur die Autokorrelation erster Ordnung ist signifikant, und wenn $ \ theta_1 <0 $ ist, ist die Autokorrelation ein negativer Wert.

Schauen wir uns als nächstes eine Verallgemeinerung des MA-Prozesses an. 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)

Die Hauptstatistiken lauten wie folgt. \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. Der $ \ quad $ MA-Prozess ist immer stabil \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.

[Ergänzung] Was ist Stationarität? Wenn weniger als $ \ quad $ gilt, ist der Prozess eine schwache Stationarität. \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 $ Zu diesem Zeitpunkt auch Autokorrelation \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 $ gilt. $ \ quad $ Sowohl Selbstkovarianz als auch Autokorrelation hängen nur von der Zeitdifferenz $ k $ ab. $ \ Quad $ Außerdem werden $ \ gamma_k = \ gamma_ {-k} $ und $ \ rho_k = \ rho_ {-k} $ eingerichtet.

4. Autoregressive Process Als nächstes folgt der AR-Prozess. Zunächst aus dem primären AR-Prozess. AR(1):y_t=c+\phi_1y_{t-1}+\epsilon_t,\quad \epsilon_t\sim W.N.(\sigma^2)

Lassen Sie es uns planen. Hier ist $ y_0 = 0 $, $ c = 1 $ und $ \ 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.0Zu diesem Zeitpunkt ist visuell leicht zu erkennen, dass der Prozess nicht stabil ist. Eigentlich,|\phi_1|\leq1Zu diesem Zeitpunkt wird der Prozess stabil. Was den Wert von $ \ phi_1 $ betrifft, so ist ersichtlich, dass der Durchschnittswert umso höher oder niedriger ist, je kleiner der Wert, desto niedriger der Durchschnittswert und je kleiner der Wert ist.

Die Hauptstatistiken lauten wie folgt. Jedoch,|\phi_1|\leq1damityIst stabil. \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)

Als nächstes folgt der verallgemeinerte AR-Prozess. 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)

Die Hauptstatistiken lauten wie folgt. \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) Die Autokorrelation des $ \ quad $ AR-Prozesses nimmt exponentiell ab.

Das Cholerogramm ist wie folgt.

5. Autoregressive Moving Average Process Die allgemeinen Begriffe des ARMA-Prozesses, der den AR-Prozess und den MA-Prozess kombiniert, sind nachstehend aufgeführt. 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) $ hat die folgenden Eigenschaften. \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 Die Autokorrelation des $ \ quad $ ARMA-Prozesses nimmt exponentiell ab.

Zeichnen wir tatsächlich $ ARMA (3,3) $. Die Parameter sind wie folgt. 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)

• Beständigkeit und Reversibilität des ARMA-Prozesses
  Der MA-Prozess ist stationär, aber der AR-Prozess ist nicht immer stationär.
  Da es mehrere MA-Prozesse mit demselben erwarteten Wert und derselben Autokorrelationsstruktur gibt, wird dies zu einem Problem bei der Modellauswahl.
  Wichtig dabei ist Stabilität und Reversibilität.

Konstanz (Dies entspricht dem stationären Zustand des AR-Prozesses.) AR characteristic equation : 1-\phi_1z-\cdots-\phi_pz^p=0 Der AR-Prozess ist stationär, wenn die Absolutwerte aller Lösungen dieser AR-Kennliniengleichung größer als 1 sind.

Reversibilität Der MA-Prozess ist reversibel, wenn der MA-Prozess in einen AR (∞) -Prozess umgeschrieben werden kann. MA characteristic equation : 1+\theta_1z+\theta_2z^2+\cdots+\theta_pz^p=0 Der MA-Prozess ist invertierbar, wenn alle Lösungen dieser MA-Kennliniengleichung größer als 1 sind. Invertierbare MA-Prozesse sind für die Modellauswahl wünschenswert.

6. Schätzung des ARMA-Modells

• Minimale Quadrate (OLS: gewöhnliche kleinste Quadrate)
  OLS wählt die Parameter des ARMA-Modells aus, um die Summe der quadratischen Residuen (SSR) zwischen den tatsächlich beobachteten Werten und den vom ARMA-Modell geschätzten Werten zu minimieren.
  \quad\displaystyle{\min_{\Theta}}\quad SSR=(y_t-ARMA(\Theta))^2
  Sie können die normale Gleichung lösen, die teilweise durch den Parameter für jeden Parameter differenziert ist.
  \quad \frac{\partial SSR}{\partial \Theta}=0

• Maximum Likelihood Estimation (MLE)
  Die wahrscheinlichste Methode ist der Ansatz zur Auswahl des Parameters, der die logarithmische Wahrscheinlichkeit maximiert.
  \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)

Von hier aus schätzen wir die Parameter für $ ARMA (3,3) $, die tatsächlich in 5 erstellt wurden. Verwenden Sie die Statistikmodellbibliothek, da die manuelle Berechnung kompliziert ist. Weitere Informationen finden Sie unter der folgenden URL. https://www.statsmodels.org/stable/generated/statsmodels.tsa.arima_model.ARMA.html Es scheint, dass, wenn Sie dem Argument von ARMA die Stichprobenstichprobe und die Reihenfolge des Modells mit der Reihenfolge = (p, q) geben, der Parameter durch die wahrscheinlichste Schätzung berechnet wird.

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

Und so weiter. Es ist möglich, einen Wert zu schätzen, der dem ursprünglichen Parameter angemessen nahe kommt. Ein paar Anmerkungen, S. D. von Innovationen, stellen die Streuung von weißem Rauschen dar, und um den Wert von c aus const zu finden,

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

0.9824998883509097

Müssen als tun.

7. ARMA-Modellauswahl

Um $ ARMA (p, q) $ auszuwählen, müssen Sie zunächst die Werte von $ p, q $ konservativ aus den Werten der Autokorrelation und der partiellen Autokorrelation eingrenzen und dann das optimale Modell basierend auf der Informationsmenge bestimmen. nehmen.

Erstens ist die partielle Autokorrelation k-ter Ordnung $ y_t $ und $ y_ {tk} $ abzüglich der Auswirkungen von $ y_ {t-1}, \ cdots, y_ {t-k + 1} $. Es ist definiert als die Korrelation von Dingen.

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

In Anbetracht der obigen Definition sind die in der folgenden Tabelle gezeigten Eigenschaften klar.

Dies ermöglicht es, das Modell in gewissem Maße einzugrenzen, indem die Struktur der Autokorrelation der Probe und der partiellen Autokorrelation der Probe betrachtet wird.

Als nächstes folgt das Informationskriterium (IC). \quad IC=-2\mathcal{L}(\hat{\Theta})+p(T)k \quad where T : number of samples, k : number of parameters Es gibt die folgenden zwei häufig verwendeten Standards für die Informationsmenge. $ \ left \ {\ begin {array} {ll} Akaike-Informationsmengenstandard (AIC): \ quad & AIC = -2 \ mathcal (L) (\ hat {\ Theta}) + 2k \\ Schwarz-Informationsmengenstandard (AIC) SIC): & SIC = -2 \ mathcal (L) (\ hat {\ Theta}) + \ log (T) k \ end {array} \ right. $ SIC tendiert dazu, ein kleineres Modell zu wählen, aber wenn die Anzahl der Stichproben ausreichend ist, kann die wahrscheinlichste Schätzung des Parameters des übermäßigen Teils von AIC gegen 0 konvergieren, so dass es eine allgemeine Regel ist, die der bessere Standard ist. Es ist nicht möglich.

8. Modellierung der tatsächlichen Daten

Wir haben Daten zu TOPIX und der Anzahl ausländischer Besucher in Japan vorbereitet, aber beide sind nicht stabil. Daher hatte ich keine andere Wahl, als Daten zu erstellen, die Trends ausschließen, indem ich den gleitenden Durchschnitt für den Teil 1996-2007 von der Anzahl ausländischer Besucher in Japan abziehe.

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

Als Trend habe ich mich für einen gleitenden Durchschnitt von 6 Monaten vorher und nachher entschieden. Dies wird durch Rollen (Fenster = 13), Mittelwert () und Verschieben (-6) verarbeitet.

Betrachten wir zunächst die Autokorrelation und die teilweise Autokorrelation. Sehr einfach mit der Bibliothek.

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

Es gibt eine starke Saisonalität, die nicht verborgen werden kann. .. .. Die Korrelation alle 12 Monate ist sehr stark. Sie können arma_order_select_ic von statsmodels verwenden, um die optimale Reihenfolge für Ihr ARMA-Modell zu schätzen. Weitere Informationen finden Sie unter diesem Link (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')

Hier wird für $ ARMA (p, q) $ der AIC-Standard verwendet, wobei der Maximalwert von $ p $ 4 und der Maximalwert von $ q $ 2 beträgt.

Ich habe das Gefühl, dass sich keiner von ihnen viel ändert, aber das Ergebnis ist, dass $ p = 4, q = 1 $ gut ist.

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

Es ist so. Persönlich hatte ich das Gefühl, dass er sein Bestes gab, obwohl er ein relativ einfaches Modell war.