概要
シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、移動平均（MA）モデルについて扱います。
MAモデル
$q$次の移動平均モデル（Moving Average model） の時刻$t$における出力は、それまでの時刻$t−1,\cdots,t−q$におけるホワイトノイズ$\epsilon_t$の重み付き和でモデリングされます。
x_t = c + \sum_{j=1}^{q} \theta_j\epsilon_{t-j} + \epsilon_t
※この重み係数$\epsilon_t$の和が1の場合、これはホワイトノイズの加重移動平均に定数を足した形となることから「移動平均」という名前が付いたそうです。
MAモデルのお気持ちとしては、過去時刻での回帰残差$\epsilon$を説明変数とする線形回帰モデル、です。
なお、$\epsilon$はショックと呼ばれたりします。
MAモデルの例
下図はMAモデルの例です。
左側は
x_t = \epsilon_{t-1} + 0.1\epsilon_{t-2} + \epsilon_t
右側は
x_t = 0.1\epsilon_{t-1} + \epsilon_{t-2} + \epsilon_t
というモデルで作成しました。
左のMAモデルは1期前のショック$\epsilon_{t-1}$を強く受けているため、1次の自己相関が強く出ます（つまり$x_t$と$x_{t-1}$がよく似ている）。一方、右のMAモデルは2期前のショック$\epsilon_{t-2}$を強く受けているため、2次の自己相関が強く出ます。
ARMAモデル
$(p,q)$次のARMAモデルとは、単に$p$次のARモデルと$q$次のMAモデルを組み合わせたものです。
x_t = c + \sum_{i=1}^{p} \phi_ix_{t-i} + \sum_{j=1}^{q} \theta_j\epsilon_{t-j} + \epsilon_t
このARMAモデルを使って、前回までと同様、気温データをモデリングしてみます。
使うデータは前回の↓と同じです。以降、データは変数dfに入っていると考えてください。
モデルの構築
statsmodelsライブラリを利用します。ARMAが無いので、ARIMAを使います。
from statsmodels.tsa.arima.model import ARIMA
# statsmodels: 0.14.0
次数を指定する引数orderには、ARIMAモデルの次数$(p,d,q)$を入れます。今回はARMAなので、「I」に相当する次数$d$はゼロにしておき、$(p,d,q)=(1,0,1)$とします。
※ARIMAモデルについては次回詳しく説明します。
つまりこれで$(1,1)$次のARMAモデルとなります。
model = ARIMA(df['temp'], order=(1, 0, 1)) # AR次数, I次数, MA次数
フィッティングは.fit()で簡単にできます。
結果は.summary()で見ることができます。
result = model.fit()
print(result.summary())
result.summary()の結果がこちら。
SARIMAX Results
==============================================================================
Dep. Variable: temp No. Observations: 8784
Model: ARIMA(1, 0, 1) Log Likelihood -10826.150
Date: Thu, 21 Dec 2023 AIC 21660.300
Time: 10:54:10 BIC 21688.623
Sample: 0 HQIC 21669.950
- 8784
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const 17.3195 1.632 10.610 0.000 14.120 20.519
ar.L1 0.9919 0.001 716.608 0.000 0.989 0.995
ma.L1 0.4545 0.006 72.151 0.000 0.442 0.467
sigma2 0.6883 0.007 102.929 0.000 0.675 0.701
===================================================================================
Ljung-Box (L1) (Q): 281.40 Jarque-Bera (JB): 4450.60
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 1.23 Skew: 0.20
Prob(H) (two-sided): 0.00 Kurtosis: 6.46
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
いろいろ情報がありますが、
-
No. Observations：データ数 -
coef：モデルの係数。-
const：定数項$c$ -
ar.L1：AR部分の係数$\phi_1$ -
ma.L1：MA部分の係数$\theta_1$ -
sigma2：$\epsilon_t$の分散（多分）
-
が主な情報でしょうか。精度などを見たい場合は、AICとかBICとかP値などを見ます。
モデルの出力
構築したモデルの出力値を見るには、.predict()を使用します。
pred = result.predict()
モデルの残差を見てみましょう。
# 残差のplot
plt.figure(figsize=(20,5))
plt.plot(result.resid)
plt.axhline(0, color='red', linestyle='--')
plt.show()
ARモデル単体のときよりも、外れたときのブレが大きくなってしまっていますが、全体としてはR2スコアは向上しています（0.989→0.992）。
ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています。
本記事のコードは01-JMA_data_analyticsの中に入っています。