時系列解析（AR, ARCH, GARCH）-実装編-

４月に入りました。春休みもおしまいですね。
今回は前回の記事（時系列解析（ARモデルからARCHモデル） - ii_da_ba_shi’s blog）を具体的に実装してみたという記事になります。

ぶっちゃけ長めです。山なし意味なしオチなしみたいな感じですが、ニッチな需要があると信じて上げます。

早速いきましょう。
使うデータはS&P500、期間は2006年から2015年です。

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime, date
import statsmodels.api as sm
from scipy.optimize import minimize
from arch import arch_model
import warnings
warnings.filterwarnings('ignore')

#############
# Variable Input
code = 'SP500'
start_y = 2006
start_m = 1
end_y = 2016
end_m = 1
#############

#############
# stock_plot
#############
name = 'csv/'+code+'.csv'
start = str(start_y)+'-'+str(start_m)
start = datetime.strptime(start, '%Y-%m')
end = str(end_y)+'-'+str(end_m)
end = datetime.strptime(end, '%Y-%m')

data = pd.read_csv(name)
data.date = pd.to_datetime(data.date)
data.set_index('date', inplace=True)
data.columns = [code]
data = data[start:end]
display(data.head())
display(data.tail())

plt.title('S&P Index(2006-2015)')
plt.xlabel('Year')
plt.ylabel('Price($)')
plt.plot(data)

f:id:ii_da_ba_shi:20190403193113p:plain — S&P Index

#############
# Log Return
#############
data = np.log(data)
N = len(data)
for i in range(1,N):
    data.iloc[N-i,0] = data.iloc[N-i,0]-data.iloc[N-i-1,0]
tmp = data.index.values[0]
data = data.drop(tmp)
display(data.head())
plt.style.use("Solarize_Light2")
plt.rcParams["axes.labelcolor"] = "#202020"
plt.rcParams["xtick.color"] = "#202020"
plt.rcParams["ytick.color"] = "#202020"
plt.rcParams["figure.dpi"] = 100
plt.title('S&P Log Return')
plt.xlabel('Year')
plt.ylabel('Daily Return')
plt.plot(data)

f:id:ii_da_ba_shi:20190403193212p:plain — Log Return

変動が激しい時期と変動が穏やかな時期が確認される。
これは分散不均一性（Heteroscedasticity）と呼ばれ、のちに紹介するARCHモデルの理論的な拠り所となっている。

AR(p)モデル

$\begin{eqnarray} R_t &=& \omega + b_1R_{t-1}+b_2R_{t-2}+\cdots +b_pR_{t-p} +\varepsilon_t　\\ \varepsilon_t &\sim& \mathrm{W.N.}(\sigma^2) \end{eqnarray}$

最小二乗法によるパラメータ推定

例：AR(1)モデル

$\begin{eqnarray} \mathrm{SSR} = \sum_{t=1}^{T} \varepsilon_t^2 = \sum_{t=1}^{T} \left(R_t-(\omega + b_1R_{t-1})\right)^2 \end{eqnarray}$

残差平方和（SSR）を最小にする $\omega,b_1$ を求めたい。
そのためにSSRをそれぞれのパラメータで偏微分したものを考える。

$\begin{eqnarray} \left\{ \begin{array}{ll} \frac{\partial \mathrm{SSR}}{\partial \omega} = -2\sum_{t=1}^{T} \left(R_t-(\omega + b_1R_{t-1})\right) = 0\\ \frac{\partial \mathrm{SSR}}{\partial b_1} = -2\sum_{t=1}^{T} R_{t-1}\left(R_t-(\omega + b_1R_{t-1})\right) = 0 \end{array} \right. \end{eqnarray}$

これらの連立方程式を解くと未知パラメータの推定量が求まる(が、実装においてはnp.linalg.lstsqを用いた)。

np.linalg.lstsqによる実装

https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html#numpy.linalg.lstsq を参照した。

numpy.linalg.lstsq(a, b, rcond='warn'）

ここでaは説明変数、bは被説明変数。

今回は定数と $R_{t-1}$ で $R_t$ を説明しようとしているので、下のようなイメージになる。

$\begin{eqnarray} \begin{bmatrix} 1 & R_0 \\ 1 & R_1 \\ \vdots & \vdots \\ 1 & R_{T-1} \end{bmatrix} \begin{bmatrix} \omega \\ b_1 \end{bmatrix} \sim \begin{bmatrix} R_1 \\ R_2 \\ \vdots \\ R_T \end{bmatrix} \end{eqnarray}$

つまり、aに左の行列を、bに右のベクトルを入れれば良さそう。
アウトプットとして得られるもの(のうち必要になるもの)はパラメータの推定値xと残差residuals。
それぞれを求めてあげれば終了となる。

#############
#AR(1) Model
#R_t = w + b_1 * R_{t-1} + eps_t
#############
T = len(data)
c = np.ones(T-1).reshape((-1, 1)) #const
r = data.iloc[0:T-1].values #realized return
a = np.hstack((c, r))
b = data.iloc[1:T].values
x,residuals,rank,s = np.linalg.lstsq(a,b,rcond=None)
#coefficient of determination
R2 = np.sum((np.dot(a,x)-np.mean(b))**2)/np.sum((b-np.mean(b))**2)
display('決定係数R2='+str(R2))
#plot
plt.scatter(r,b)
plt.plot(r,np.dot(a,x),color='orange')

f:id:ii_da_ba_shi:20190403200242p:plain — AR(1)

決定係数 $R^2=0.0104$ と酷い有様である。
収益率データ系列について自己回帰分析をするのは筋が悪いのが見てとれる。
念の為コレログラムも以下に掲載しておく。

#############
#Correlogram
#############
fig = plt.figure()
ax1 = fig.add_subplot(111)
sm.graphics.tsa.plot_acf(data, lags=20, ax=ax1)
plt.show()

f:id:ii_da_ba_shi:20190403200343p:plain — Correlogram(AR)

残差の2乗の時系列分析

決定係数・コレログラムを見るかぎり、収益率をそれ自身の自己回帰モデルで表すのは無理がありそう。
その代わり、残差（言い換えるなら、マーケットの期待から外れた大きさ。カタカナ語で言うならショック）は自己回帰モデルによってモデリングできるかもしれない。
残差の2乗に対してAR(q)モデルを考えてみる。

$\begin{eqnarray} \varepsilon^2_t = \omega + a_1\varepsilon^2_{t-1} + a_2\varepsilon^2_{t-2} + \cdots a_q\varepsilon^2_{t-q} \end{eqnarray}$

#AR(4)
p = 4
eps = (b-np.dot(a,x))**2
T = len(eps)
c = np.ones(T-p).reshape((-1, 1)) #const
e1 = eps[p-1:T-1] #eps1
e2 = eps[p-2:T-2] #eps2
e3 = eps[p-3:T-3] #eps2
e4 = eps[p-4:T-4] #eps2
a = np.hstack([c,e1,e2,e3,e4])
b = eps[p:T]
x,residuals,rank,s = np.linalg.lstsq(a,b,rcond=None)
#coefficient of determination
R2 = np.sum((np.dot(a,x)-np.mean(b))**2)/np.sum((b-np.mean(b))**2)
display('決定係数R2='+str(R2))

決定係数は $R^2=0.2178$ と多少マシになった。
コレログラムからもある程度の相関は確認される。

fig = plt.figure()
ax1 = fig.add_subplot(111)
sm.graphics.tsa.plot_acf(eps, lags=20, ax=ax1)
plt.show()

f:id:ii_da_ba_shi:20190403200924p:plain — Correlogram(ARCH)

ARCH(q)モデル

上記の「残差の2乗に対してARモデルを適用して得られた、残差の推定値」をボラティリティの推定値とする、というのが分散自己回帰モデルの概要だ。
式で言うなら

$\begin{eqnarray} \varepsilon^2_t &=& \omega + a_1\varepsilon^2_{t-1} + a_2\varepsilon^2_{t-2} + \cdots a_q\varepsilon^2_{t-q} \\ \sigma^2_t &:=& \varepsilon^2_t \end{eqnarray}$

ということ。
このモデルではボラティリティクラスタリングを表現できているのが強みである。
T期での収益率は「T-1期での収益率の期待値 $\mathbb{E}(R_{t-1})$ 」と「予期出来ない変動 $\varepsilon_t$ 」の和として表される。
後者に関して、その変動のボラティリティはARCHモデルによって判明している。
予期できない変動は確率項を含むので、平均が0、分散が1でi.i.d.な確率変数 $z_t$ を考えることで $\varepsilon_t = \sigma_t + z_t$ と表せる。
以上をまとめると

$\begin{eqnarray} R_t &=& \mathbb{E}(R_{t-1}) + \varepsilon_t \\ \varepsilon_t &=& \sigma_t z_t \\ z_t &\sim& \mathrm{i.i.d.},\ \mathbb{E}(z_t) = 0, \ \mathrm{Var}(z_t) = 1 \end{eqnarray}$

GARCH(p,q)モデル

ボラティリティの説明変数に「過去の予期できない変動の2乗」だけでなく、「過去の観測されたボラティリティ」を加えたものがG(Generalized)ARCHモデルと呼ばれるモデルだ。
数式で表すと以下のようになる。

$\begin{eqnarray} \sigma^2_t = \omega + \sum_{i=1}^{p}\beta_i\sigma^2_{t-i} +\sum_{j=1}^{q}\alpha_j\varepsilon^2_{t-j} \end{eqnarray}$

GARCHモデルのパラメータ推定

今回は簡単のため（また、実証分析においてもフィッティングがよく採用される頻度が高い）GARCH(1,1)モデルを考える。
念の為、数式で表しておくと以下のようになる。

$\begin{eqnarray} \sigma^2_t = \omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1} \end{eqnarray}$

推定したいパラメータは $\omega,\beta,\alpha$ の３つ
T期ではT-1期でのボラティリティと残差が既知になっている

というのが特徴として挙げられる
「T期で観測された誤差の尤度が最大になるようにパラメータを推定する」ことが今回の肝になる。
もう少し具体的に噛み砕いて言うなら、以下のようになる。

T-1期の既知データ $\sigma_{T-1},\varepsilon_{T-1}$ を用いることで、T期でのボラティリティの予測 $\hat{\sigma}^2_T$ が得られる。
T期では誤差 $\varepsilon_T$ が観測されるが、この値が先ほどの予測からしてどの程度あり得そうなのか（＝尤度 $f(\varepsilon_T|\hat{\sigma}^2_T)$ ）が得られる。
誤差関数の中身である $z_T$ が標準正規分布に従うという仮定をひとまず立てておく。
すると、尤度関数 $f$ は平均0,分散 $\sigma^2_T$ の正規分布の確率密度関数となる。
$T=0$ から順に尤度を求めていき、それらの積を最大にするようなパラメータを求める。

最尤推定においては尤度 $L$ そのものではなく、対数を取った「対数尤度 $\ln L$ 」が使われることが多い。
最後にこれらの定式化をして締めくくる。

$\begin{eqnarray} L &=& \prod_{t=1}^{T} f(\varepsilon_t|\hat{\sigma}^2_t) \\ f(\varepsilon_t|\hat{\sigma}^2_t) &=& \frac{1}{\sqrt{2\pi\hat{\sigma}^2_t}}\exp{\left(-\frac{\varepsilon^2_t}{2\hat{\sigma}^2_t} \right)} \\ \ln L &=& -\frac{T}{2}\ln (2\pi) - \frac{1}{2} \sum_{t=1}^{T} \ln (\hat{\sigma}^2_t) - \frac{1}{2} \sum_{t=1}^{T} \frac{\varepsilon^2_t}{\hat{\sigma}^2_t} \\ \hat{\sigma}^2_t &=& \omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1} \end{eqnarray}$

尤度最大化

続いて尤度関数を最大化するパラメータ $\omega,\beta,\alpha$ を求めることになる。

まずは対数尤度をそれぞれのパラメータで偏微分する。

$\begin{eqnarray} \frac{\partial \ln L}{\partial \omega} &=& -\frac{1}{2}\sum_{t=1}^{T}\frac{1}{\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1}} + \frac{1}{2}\sum_{t=1}^{T} \frac{\varepsilon^2_t}{(\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1})^2} = \frac{1}{2}\sum_{t=1}^{T}\frac{\varepsilon^2_t - \omega - \beta\sigma^2_{t-1} - \alpha\varepsilon^2_{t-1}}{(\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1})^2} \\ \frac{\partial \ln L}{\partial \beta} &=& -\frac{1}{2}\sum_{t=1}^{T}\frac{\sigma^2_{t-1}}{\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1}} + \frac{1}{2}\sum_{t=1}^{T} \frac{\varepsilon^2_t\sigma^2_{t-1}}{(\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1})^2} = \frac{1}{2}\sum_{t=1}^{T}\frac{\sigma^2_{t-1}(\varepsilon^2_t - \omega - \beta\sigma^2_{t-1} - \alpha\varepsilon^2_{t-1})}{(\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1})^2}\\ \\ \frac{\partial \ln L}{\partial \alpha} &=& -\frac{1}{2}\sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}}{\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1}} + \frac{1}{2}\sum_{t=1}^{T} \frac{\varepsilon^2_t\varepsilon^2_{t-1}}{(\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1})^2}=\frac{1}{2}\sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}(\varepsilon^2_t - \omega - \beta\sigma^2_{t-1} - \alpha\varepsilon^2_{t-1})}{(\omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1})^2} \end{eqnarray}$

それぞれの偏微分の値を0にするような連立方程式を考え、それを解く。
分母に現れている未知パラメータ項については、これが真のボラティリティに等しいとみなして計算しやすくする。
つまり、 $\sigma^2_t := \hat{\sigma}^2_t = \omega + \beta\sigma^2_{t-1} + \alpha\varepsilon^2_{t-1}$ ということ。

$\begin{eqnarray} \sum_{t=1}^{T}\frac{\varepsilon^2_t}{\sigma^4_t}-\sum_{t=1}^{T}\frac{1}{\sigma^4_t}\omega-\sum_{t=1}^{T}\frac{\sigma^2_{t-1}}{\sigma^4_t}\beta-\sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}}{\sigma^4_t}\alpha = 0 \\ \sum_{t=1}^{T}\frac{\varepsilon^2_t\sigma^2_{t-1}}{\sigma^4_t}-\sum_{t=1}^{T}\frac{\sigma^2_{t-1}}{\sigma^4_t}\omega-\sum_{t=1}^{T}\frac{\sigma^4_{t-1}}{\sigma^4_t}\beta-\sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}\sigma^2_{t-1}}{\sigma^4_t}\alpha = 0 \\ \sum_{t=1}^{T}\frac{\varepsilon^2_t\varepsilon^2_{t-1}}{\sigma^4_t}-\sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}}{\sigma^4_t}\omega-\sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}\sigma^2_{t-1}}{\sigma^4_t}\beta-\sum_{t=1}^{T}\frac{\varepsilon^4_{t-1}}{\sigma^4_t}\alpha = 0 \end{eqnarray}$

行列の形式に書き直して終わる。

$\begin{eqnarray} \displaystyle \begin{pmatrix} \displaystyle \sum_{t=1}^{T}\frac{1}{\sigma^4_t}& \displaystyle \sum_{t=1}^{T}\frac{\sigma^2_{t-1}}{\sigma^4_t} & \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}}{\sigma^4_t} \\ \displaystyle \sum_{t=1}^{T}\frac{\sigma^2_{t-1}}{\sigma^4_t}& \displaystyle \sum_{t=1}^{T}\frac{\sigma^4_{t-1}}{\sigma^4_t} & \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}\sigma^2_{t-1}}{\sigma^4_t} \\ \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}}{\sigma^4_t}& \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_{t-1}\sigma^2_{t-1}}{\sigma^4_t} & \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^4_{t-1}}{\sigma^4_t} \\ \end{pmatrix} \begin{pmatrix} \omega \\ \beta \\ \alpha \end{pmatrix} = \begin{pmatrix} \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_t}{\sigma^4_t} \\ \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_t\sigma^2_{t-1}}{\sigma^4_t} \\ \displaystyle \sum_{t=1}^{T}\frac{\varepsilon^2_t\varepsilon^2_{t-1}}{\sigma^4_t} \end{pmatrix} \end{eqnarray}$

いざ実装である。
特筆すべき点は何もないが、分散を求めるにあたって、1次と2次のモーメントの移動平均を用いる方法を採った。（スッキリ書けそうだったので）

##########
#data set
##########
#realized return
display('収益率')
display(data.head())
#Moving Average(1st order moment)
exp = data.rolling(100).mean()
# 2nd order moment
sq = data ** 2
exp_sq = sq.rolling(100).mean()
sig = exp_sq - exp ** 2

#AR(2)
T = len(data)
c = np.ones(T-2).reshape((-1, 1)) #const
r1 = data.iloc[0:T-2].values #realized return
r2 = data.iloc[1:T-1].values #realized return
a = np.hstack((c, r1, r2))
b = data.iloc[2:T].values
x,residuals,rank,s = np.linalg.lstsq(a,b,rcond=None)
eps = data.iloc[2:T] - np.dot(a,x)
eps = eps ** 2

eps = eps.iloc[100-2:,]
sig = sig.iloc[100:,]

display('誤差')
display(eps.head())
display('分散')
display(sig.head())

# Ax = b
L = len(eps)
A = np.zeros((3,3))
b = np.zeros(3)
for i in range(1,L):
    s2 = sig.iloc[i,0]**2
    A[0,0] += 1/s2
    A[0,1] += sig.iloc[i-1,0]/s2
    A[0,2] += eps.iloc[i-1,0]/s2
    A[1,0] += sig.iloc[i-1,0]/s2
    A[1,1] += sig.iloc[i-1,0]**2/s2
    A[1,2] += sig.iloc[i-1,0]*eps.iloc[i-1,0]/s2
    A[2,0] += eps.iloc[i-1,0]/s2
    A[2,1] += sig.iloc[i-1,0]*eps.iloc[i-1,0]/s2
    A[2,2] += eps.iloc[i-1,0]**2/s2
    b[0] += eps.iloc[i,0]/s2
    b[1] += eps.iloc[i,0]*sig.iloc[i-1,0]/s2
    b[2] += eps.iloc[i,0]*eps.iloc[i-1,0]/s2
x = np.linalg.solve(A,b)
display('omega,beta,alpha')
display(x)


#x = np.array([3.4833e-06,0.88,0.1])
lnL = 0
for i in range(1,len(eps)):
    s0 = np.log(2*np.pi)
    s1 = np.log1p(x[0] + x[1]*sig.iloc[i-1,0] + x[2]*eps.iloc[i-1,0])
    s2 = eps.iloc[i,0]/(x[0] + x[1]*sig.iloc[i-1,0] + x[2]*eps.iloc[i-1,0])
    lnL += (s0+s1+s2)/2
display(lnL)

という具合に推定値が求まる。
各パラメータの推定値は以下のようになった。

$(\omega,\beta,\alpha )= (1.195\times 10^{-5}, 0.7514,0.1305)$

また、PythonにはARCHモデルのライブラリもあって、今までやってたようなことを一行で実現してくれる上に、標準誤差とかAICも出してくれる。（以下を参照）
ここまでの努力はなんだったんだ。ぷんぷん。
強いていうなら、ライブラリではイテレーションを用いているけど、今回実装したモデルは最尤推定量を解析的に求めている点で強みかなあとは思う。

garch=arch_model(data.iloc[100:,0],mean='AR',lags=2,vol='GARCH',p=1,o=0,q=1,dist='normal').fit()
display(garch.summary)

同じデータを用いて推定した結果が以下。 $\beta$ の推定を大外ししているが、その他は妥当な範囲だろうか。

f:id:ii_da_ba_shi:20190403230915p:plain — Summary

どうでもしか勝たん

格が違うので卍卍卍