【気象データで時系列解析⑤】MAモデル・ARMAモデル

# 概要
シリーズ「気象データで時系列解析」では、気象データを例に時系列解析の基礎を学びます。今回は、**移動平均（MA）モデル**について扱います。

# MAモデル

$q$次の**移動平均モデル（Moving Average model）** の時刻$t$における出力は、それまでの時刻$t−1,\cdots,t−q$におけるホワイトノイズ$\epsilon_t$の重み付き和でモデリングされます。

```math
x_t = c + \sum_{j=1}^{q} \theta_j\epsilon_{t-j} + \epsilon_t
```

※この重み係数$\epsilon_t$の和が1の場合、これはホワイトノイズの加重移動平均に定数を足した形となることから「移動平均」という名前が付いたそうです。

MAモデルのお気持ちとしては、過去時刻での回帰残差$\epsilon$を説明変数とする線形回帰モデル、です。

なお、$\epsilon$は**ショック**と呼ばれたりします。

## MAモデルの例
下図はMAモデルの例です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/3a98f0e6-3b6e-b522-0639-959b8de6924c.png)

左側は

```math
x_t = \epsilon_{t-1} + 0.1\epsilon_{t-2} + \epsilon_t
```
右側は

```math
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次の自己相関が強く出ます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/e75be7ee-7fc1-9e22-45a7-fe9e7ee4a403.png)

# ARMAモデル

$(p,q)$次の**ARMAモデル**とは、単に$p$次のARモデルと$q$次のMAモデルを組み合わせたものです。

```math
x_t = c + \sum_{i=1}^{p} \phi_ix_{t-i} + \sum_{j=1}^{q} \theta_j\epsilon_{t-j} + \epsilon_t
```

このARMAモデルを使って、前回までと同様、気温データをモデリングしてみます。

使うデータは前回の↓と同じです。以降、データは変数`df`に入っていると考えてください。

https://qiita.com/matsuda_tkm/items/8f0ca38e3b6e03de192d

### モデルの構築
`statsmodels`ライブラリを利用します。ARMAが無いので、ARIMAを使います。

```python
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モデルとなります。

```python
model = ARIMA(df['temp'], order=(1, 0, 1))  # AR次数, I次数, MA次数
```

フィッティングは`.fit()`で簡単にできます。
結果は`.summary()`で見ることができます。

```python
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()`を使用します。

```python
pred = result.predict()
```

モデルの残差を見てみましょう。

```python
# 残差のplot
plt.figure(figsize=(20,5))
plt.plot(result.resid)
plt.axhline(0, color='red', linestyle='--')
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/5c89e874-d127-7cc7-b42a-e0c3f802699b.png)



ARモデル単体のときよりも、外れたときのブレが大きくなってしまっていますが、全体としてはR2スコアは向上しています（0.989→0.992）。

# ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています。

本記事のコードは`01-JMA_data_analytics`の中に入っています。

https://github.com/matsuda-tkm/qiita-calender-2023