概要

シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、移動平均（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の中に入っています。