【気象データで時系列解析④】ARモデル

# 概要
シリーズ「気象データで時系列解析」では，気象データを例に時系列解析の基礎を学びます．今回は，**自己回帰（AR）モデル**について扱います．

# ARモデル
$p$次の**自己回帰モデル（AutoRegressive model）** では，時刻$t$における値を，その前の時刻$t-1,\cdots,t-p$までの値の線形和でモデリングします．すなわち，

```math
x_t = c + \sum_{i=1}^p \phi_i x_{t-i} + \epsilon_t \tag{1}
```

というモデルで，$\phi_i\ (i=1,\cdots,p)$はモデルパラメータ，$c$は定数項，$\epsilon_t$はホワイトノイズ[^1]です．

つまりARモデルは，説明変数を前時刻までのデータとする重回帰モデルということになります．

なお，$p$次のARモデルのことをよく「$AR(p)$」と表記するので覚えておきましょう．

# ARモデルの性質
ARモデルが弱定常となる条件を調べてみます．

AR(1)モデル，つまり$p=1$としたときの式$(1)$を繰り返し適用して展開してみましょう．

```math
\begin{align*}
x_t &= c + \phi_1 x_{t-1} + \epsilon_t\\
&= c + \phi_1 (c + \phi_1 x_{t-2} + \epsilon_{t-1}) + \epsilon_t\\
&= c(1+\phi_1) + \phi_1^2 x_{t-2} + \epsilon_t + \phi_1\epsilon_{t-1}\\
&= c(1+\phi_1) + \phi_1^2(c + \phi_1 x_{t-3} + \epsilon_{t-2}) + \epsilon_t + \phi_1\epsilon_{t-1}\\
&= c(1+\phi_1+\phi_1^2) + \phi_1^3x_{t-3} + \epsilon_t + \phi_1\epsilon_{t-1} + \phi_1^2\epsilon_{t-2}\\
&= \cdots\\
&= c(1+\phi_1+\phi_1^2+\cdots+\phi_1^{t-1}) + \phi_1^tx_0 + (\epsilon_t + \phi_1\epsilon_{t-1} + \phi_1^2\epsilon_{t-2} + \cdots + \phi_1^{t-1}\epsilon_1)\\
&= c\cdot\frac{1-\phi_1^t}{1-\phi_1} + \phi_1^tx_0 + (\epsilon_t + \phi_1\epsilon_{t-1} + \phi_1^2\epsilon_{t-2} + \cdots + \phi_1^{t-1}\epsilon_1)
\end{align*}
```
$\phi_1= 1$なら，

```math
x_t = ct + x_0 + (\epsilon_t + \epsilon_{t-1} + \epsilon_{t-2} + \cdots + \epsilon_1)
```
となりますが，$c\ne 0$のとき$t\to\infty$の極限で$x_t\to\pm\infty$と発散してしまいます．$c=0$のときは**ランダムウォーク**となります（ランダムウォークは非定常過程）．

一方$\phi_1\ne 1$なら，

```math
x_t = c\cdot\frac{1-\phi_1^t}{1-\phi_1} + \phi_1^tx_0 + (\epsilon_t + \phi_1\epsilon_{t-1} + \phi_1^2\epsilon_{t-2} + \cdots + \phi_1^{t-1}\epsilon_1)
```
となります．$t\to\infty$の極限では$|\phi_1|>1$なら発散しますが，$|\phi_1|<1$なら収束します．

よって，AR(1)モデルが定常であるための必要条件は

```math
|\phi_1|<1
```
です．一般には，AR(p)が定常であるための必要条件は，$z$に関する方程式

```math
z^p - \sum_{i=1}^p \phi_i z^{p-i} = 0
```
のすべての解$z_i$が$|z_i|<1$を満たすことである，ということを示せます．

[^1]:$\epsilon_t$が**ホワイトノイズ**$\iff$ 任意の$t,j$に対し$\mathbb{E}[\epsilon_t]=0, \mathbb{V}[\epsilon_t]=\sigma^2, \text{Cov}[\epsilon_t, \epsilon_{t-j}]=0$


# 気温データをARモデルでモデリングしてみる

気象庁の「[過去の気象データダウンロード](https://www.data.jma.go.jp/gmd/risk/obsdl/index.php)」から取得した気温データを，ARモデルでモデリングしてみましょう．

>地点：京都
>期間：2022年11月26日～2023年11月26日
>データの種類：時別値

前処理方法は割愛します．
変数`df`に下表のようなデータが入っていると思ってください．
`df`はPandasのDataFrame型です．

| datetime | temp |
|:-:|:-:|
| 2022-11-26 01:00:00 | 10.9 |
| 2022-11-26 02:00:00 | 10.9 |
| 2022-11-26 03:00:00 | 10.7 |
| 2022-11-26 04:00:00 | 10.7 |
| 2022-11-26 05:00:00 | 10.6 |
| ... | ... |

### 時間をずらしたデータを作る
`df['temp']`を1～pだけずらしたデータ`df['lag1']`～`df['lagp']`を作成します．

これを行うには，Pandasの`.shift()`メソッドが便利です．

```python
P = 1  # ARモデルの次数

df_lag = df[['temp']].copy()  # dfからカラム「temp」を取り出し，df_lagへ格納
for lag in range(1, P+1):
    df_lag[f'lag{lag}'] = df_lag.temp.shift(lag)  # カラム「temp」をlagだけ下へずらす
df_lag = df_lag.dropna()  # ずらすと先頭にnanが出るのでdropする
```

これで，1～P個ずらしたデータがまとめて`df_lag`へDataFrameとして格納されました．


### モデル構築
ARモデルは結局重回帰モデルと同等なので，`sklearn`の`LinearRegression`で構築します．

```python
from sklearn.linear_model import LinearRegression

# 線形回帰モデルの構築
ar = LinearRegression()
ar.fit(df_lag.drop('temp', axis=1), df_lag['temp'])
```

### モデリング結果

モデリング結果を見てみましょう．ここでは残差をプロットします．

```python
# 予測値の計算
pred = ar.predict(df_lag.drop('temp', axis=1))

# 評価
print(f'R^2 score: {ar.score(df_lag.drop("temp", axis=1), df_lag["temp"])}')

# 残差のプロット
plt.figure(figsize=(20, 5))
plt.plot(pred-df_lag['temp'], label='residual')
plt.legend()
plt.show()
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1900599/96e14b1b-49b6-e68c-2bfb-305b930c8cc1.png)

どの時刻でも残差はおおむね0まわりを同じような振幅で振動しており，ホワイトノイズに近いように見えます（ざっくりですが）．

### statsmodelsでARモデルを楽に構築
実は，`statsmodels`というライブラリを使えば一発でARモデルを構築できます．

```python
from statsmodels.tsa.ar_model import AutoReg

# ARモデルの作成
P = 1
ar2 = AutoReg(df['temp'], lags=P, old_names=False)
ar2_fit = ar2.fit()

# 予測
pred2 = ar2_fit.predict(start=P, end=len(df)-1)
```

念のため，自分の実装とライブラリの結果を比較しましたが，きちんと合致しているようです．
```python
(pred - pred2).min(), (pred - pred2).max()
# (-6.217248937900877e-15, 1.4210854715202004e-14)
```

# ソースコード

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

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

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