概要
シリーズ「気象データで時系列解析」では，気象データを例に時系列解析の基礎を学びます．今回は，自己回帰（AR）モデルについて扱います．
ARモデル
$p$次の自己回帰モデル（AutoRegressive model） では，時刻$t$における値を，その前の時刻$t-1,\cdots,t-p$までの値の線形和でモデリングします．すなわち，
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)$を繰り返し適用して展開してみましょう．
\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$なら，
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$なら，
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)モデルが定常であるための必要条件は
|\phi_1|<1
です．一般には，AR(p)が定常であるための必要条件は，$z$に関する方程式
z^p - \sum_{i=1}^p \phi_i z^{p-i} = 0
のすべての解$z_i$が$|z_i|<1$を満たすことである，ということを示せます．
気温データをARモデルでモデリングしてみる
気象庁の「過去の気象データダウンロード」から取得した気温データを，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()メソッドが便利です．
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で構築します．
from sklearn.linear_model import LinearRegression
# 線形回帰モデルの構築
ar = LinearRegression()
ar.fit(df_lag.drop('temp', axis=1), df_lag['temp'])
モデリング結果
モデリング結果を見てみましょう．ここでは残差をプロットします．
# 予測値の計算
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()
どの時刻でも残差はおおむね0まわりを同じような振幅で振動しており，ホワイトノイズに近いように見えます（ざっくりですが）．
statsmodelsでARモデルを楽に構築
実は，statsmodelsというライブラリを使えば一発でARモデルを構築できます．
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)
念のため，自分の実装とライブラリの結果を比較しましたが，きちんと合致しているようです．
(pred - pred2).min(), (pred - pred2).max()
# (-6.217248937900877e-15, 1.4210854715202004e-14)
ソースコード
Qiitaアドベントカレンダーで使ったコードはこちらのレポジトリにまとめています．
本記事のコードは01-JMA_data_analyticsの中に入っています．
-
$\epsilon_t$がホワイトノイズ$\iff$ 任意の$t,j$に対し$\mathbb{E}[\epsilon_t]=0, \mathbb{V}[\epsilon_t]=\sigma^2, \text{Cov}[\epsilon_t, \epsilon_{t-j}]=0$ ↩