概要

シリーズ「気象データで時系列解析」では，気象データを例に時系列解析の基礎を学びます．今回は，自己回帰（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$はホワイトノイズ¹です．

つまり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$ ↩