コードで理解するVARモデル - えんじにあのじゆうちょう

はじめに

VARモデルのとりあえず推定までをやってみることで、どういうモデルなのか、ということの理解を試みます。

今回も経済・ファイナンスデータの軽量時系列分析を勉強しながら書いたものなので、数式などの詳細はそちらを読むこと。

やってみる

VARモデルの数式

VARモデルはAR系のモデルでは単変量の時系列データしか扱えませんでしたが、複数のものを関係のあるものとして同時に扱いたい場合に使います。
例えば今回、時刻tの売上を $y^{(1)}_t$ , 同時刻の客数を $y^{(2)}_t$ とすると、 $y^{(1)}_t$ を予測するときに $y^{(1)}_{t-1}, \cdots y^{(1)}_{t-p}$ $y^{(2)}_{t-1}, \cdots y^{(2)}_{t-p}$ を使い、同様に $y^{(2)}_t$ を $y^{(1)}_{t-1}, \cdots y^{(1)}_{t-p}$ $y^{(2)}_{t-1}, \cdots y^{(2)}_{t-p}$ を使い表現します。
今回簡単のためにp=1とすると、数式で真面目に書くと以下のようになります。

$y^{(1)}_{t} = c_1 + \phi_{11} y^{(1)}_{t-1} + \phi_{12} y^{(2)}_{t-1}$
$y^{(2)}_{t} = c_2 + \phi_{21} y^{(1)}_{t-1} + \phi_{22} y^{(2)}_{t-1}$

$\bf{\phi}$ を求めるとき、実際の測定値との差分の二乗和誤差(Sum of Squared Residual, SSR)を計算してそれを最小にするようにします。
SSRは以下の式で求めることができます。

$SSR = \sum (y_t - \hat{y_{t}})^2$

これってE資格のときに散々やった2次関数なので微分して傾き分動かしていけばパラメータを勾配降下法で求められそうなので微分しておきます。

$\frac{\partial SSR}{\partial c_1} = -2\sum(y_t - \hat{y_{t}}) = -2\sum(y_t^{(1)} - c_1 - \phi_{11} y^{(1)}_{t-1} - \phi_{12} y^{(2)}_{t-1})$
$\frac{\partial SSR}{\partial c_2} = -2\sum(y_t - \hat{y_{t}}) = -2\sum(y_t^{(2)} - c_2 - \phi_{21} y^{(1)}_{t-1} - \phi_{22} y^{(2)}_{t-1})$
$\frac{\partial SSR}{\partial \phi_{11}} = -2\sum y^{(1)}_{t-1} (y_t^{(1)} - c_1 - \phi_{11} y^{(1)}_{t-1} - \phi_{12} y^{(2)}_{t-1})$
$\frac{\partial SSR}{\partial \phi_{12}} = -2\sum y^{(2)}_{t-1} (y_t^{(1)} - c_1 - \phi_{11} y^{(1)}_{t-1} - \phi_{12} y^{(2)}_{t-1})$
$\frac{\partial SSR}{\partial \phi_{21}} = -2\sum y^{(1)}_{t-1} (y_t^{(2)} - c_2 - \phi_{21} y^{(1)}_{t-1} - \phi_{22} y^{(2)}_{t-1})$
$\frac{\partial SSR}{\partial \phi_{22}} = -2\sum y^{(2)}_{t-1} (y_t^{(2)} - c_2 - \phi_{21} y^{(1)}_{t-1} - \phi_{22} y^{(2)}_{t-1})$

学習率を $\eta$ として、

$c_1 = c_1 - \eta \frac{\partial SSR}{\partial c_1}$
$c_2 = c_2 - \eta \frac{\partial SSR}{\partial c_2}$
$\phi_{11} = \phi_{11}- \eta \frac{\partial SSR}{\partial \phi_{11}}$
$\phi_{12} = \phi_{12}- \eta \frac{\partial SSR}{\partial \phi_{12}}$
$\phi_{21} = \phi_{21}- \eta \frac{\partial SSR}{\partial \phi_{21}}$
$\phi_{22} = \phi_{22}- \eta \frac{\partial SSR}{\partial \phi_{22}}$

このようにして、逐次更新していけば求められそうです。

実装の前に使うデータ

使用するデータは以下からダウンロードしたmcsi_day.xlsをcsv化して使うことにしました。

www.asakura.co.jp

実装

データの読み込みと整形をまず実施します。
今回はjp, uk, usのみ使って、対数差分系列を作ります。

df = pd.read_csv("./msci.csv",skiprows=1, names=['Date','jp','uk','us']).dropna()
data = pd.DataFrame(index=df.index)
data['jp'] = np.log(df['jp']) - np.log(df['jp'].shift(1))
data['uk'] = np.log(df['uk']) - np.log(df['uk'].shift(1))
data['us'] = np.log(df['us']) - np.log(df['us'].shift(1))

先に、現在のパラメータを受け取って予測をするpred関数と、SSRを計算するerror関数を定義します。

def pred(y, phi, c):
    return c + np.dot(phi, y)
def error(y_t, y_t1, phi, c):
    return np.sum((y -pred(y_t1, phi, c)) ** 2)

そして学習のコードです。前処理が長ったらしい感じで、学習自体は数行です。
なお、今回みたいな幅の大きいデータの場合先に正規化しないとうまく学習できなかったので正規化もいれています。

# パラメータ設定
N = 3 # 説明変数の数
p = 3 #VAR(p)
eta = 0.0001
iter_num = 100000

y = df[["jp", "uk", "us"]].to_numpy()
print(y.shape)

# 正規化
mean = y.mean(axis=0)
std = y.std(axis=0)
y = (y - mean ) / std

# 初期化
phi = np.random.normal(size=(N,N*p))
c = np.random.normal(size=(N, 1))

# 1, 2, ..., p時点前のデータ列を生成
y_diff = []
for j in range(N):
    for i in range(1, p+1):
        y_diff.append(y[p-i:-i, j])
y_diff = np.array(y_diff)
y = y[p:].T
print(y_diff.shape)
print(y.shape)
    
best = [None, None]
best_error = np.inf

# 推定はここだけ
for i in range(iter_num):
    dc = (- np.sum((y - (c + np.dot(phi, y_diff))), axis=1))[:, np.newaxis]
    dphi = - np.dot(y_diff, (y - (c + np.dot(phi, y_diff))).T).T
    c = c - eta * dc
    phi = phi - eta * dphi
    e = error(y, y_diff, phi, c)
    if  e < best_error:
        best = [phi, c]
        best_error = e
phi, c = best
preds = pred(y_diff, phi, c) * std[:, np.newaxis] + mean[:, np.newaxis]

とりあえずできたということで可視化していきます。

plt.plot(np.arange(1391), df["jp"], color="blue", label="real")
plt.plot(np.arange(p,1391), preds[0, :], color="red", label="pred")
plt.legend()

plt.plot(np.arange(1391), df["uk"], color="blue", label="real")
plt.plot(np.arange(p,1391), preds[1, :], color="red", label="pred")
plt.legend()

plt.plot(np.arange(1391), df["us"], color="blue", label="real")
plt.plot(np.arange(p,1391), preds[2, :], color="red", label="pred")
plt.legend()