はじめに

時系列解析については以前にMCMCを用いた状態空間モデルの推定を行なったのですが、状態空間モデルの推定方法としてカルマンフィルタも知っておいた方が良さそうだったので、今回はカルマンフィルタの実装を行なっていきます。

状態空間モデルの推定方法について

状態空間モデルの推定方法はいくつかあるのですが、その内訳は下記に詳細に記載されていました。

状態空間モデルの推定方法の分類 | Logics of Blue

MCMCとカルマンフィルタについての大枠の特徴はこのような形になります。

カルマンフィルタ
- 逐次処理
- ノイズがガウス分布に従う仮定
- 推定が早い
MCMC
- 一括処理
- ノイズがガウス分布に従ってなくても適用可
- 推定が遅い

推定の速さが肝になる自動車やロボット制御ではカルマンフィルタが使われ、特に時間の制約や逐次処理の必要がない時系列解析ではMCMCが使用されることが多いのはこの特性のためですね。

実験

では実際にカルマンフィルタを用いた状態空間モデルの推定を実施していきます。

今回はイチローの年別安打数を予測するモデルを作成したいと思います。データソースは以下になります。

図録▽イチローの年度別安打（ヒット）数の推移

データは以下のような形になります。0に収束してしまうと面白くないので最後の五年は削ってモデルに入力します。

f:id:rmizutaa:20190409084445p:plain

カルマンフィルタのコードも下記サイトを参考にさせていただきました。

最も簡単なローカルレベルモデルに対してカルマンフィルタを適用します。ローカルトレンドモデルで推定が必要なパラメータは状態方程式のノイズの分散(sigmaW)と観測方程式のノイズの分散(sigmaV)の２つで、これらを最尤推定を用いて推定します。

def localLevelModel(y, xPre, pPre, sigmaW, sigmaV):
    # 状態の予測(ローカルレベルモデルなので、予測値は、前期の値と同じ)
    xForecast = xPre

    # 状態の予測誤差の分散
    pForecast = pPre + sigmaW

    # カルマンゲイン
    kGain = pForecast / (pForecast + sigmaV)

    # カルマンゲインを使って補正された状態
    xFiltered = xForecast + kGain * (y - xForecast)

    # 補正された状態の予測誤差の分散
    pFiltered = (1 - kGain) * pForecast

    # 観測値の予測誤差
    v = y - xForecast

    # 観測値の予測誤差の分散
    F = pForecast + sigmaV
    
    return xFiltered,pFiltered,v,F

#尤度の計算
def norm_dens(val, s):
    #val:値
    #s:標準偏差
    return np.log((1/np.sqrt(2*np.pi*s**2))*np.exp(-0.5*(val)**2/s**2))


def kalmanfilter(df,sigmaW=1000,sigmaV=10000,optimize=False):
    x=np.zeros(df.shape[0]+1)
    P=np.zeros(df.shape[0]+1)
    v=np.zeros(df.shape[0])
    F=np.zeros(df.shape[0])
    #状態の予測誤差の分散」の初期値
    P[0]=10000000
    for i in range(df.shape[0]):
        x[i + 1],P[i + 1],v[i],F[i] = localLevelModel(df[1].values[i], x[i], P[i], sigmaW = sigmaW, sigmaV = sigmaV)

    #尤度の計算
    if optimize==True:
        return -norm_dens(v, np.sqrt(F)).sum()
    else:
        return x,P,v,F

#最尤推定の実施
from scipy.optimize import fmin
import math
yu = lambda X: kalmanfilter(df,X[0],X[1],optimize=True)
result=fmin(yu, [1000, 1000])

この結果、sigmaW,sigmaV=[1635.79921897, 594.56116963]と推定できました。

このパラメータを用いて予測した値と実際の値を比較します。

f:id:rmizutaa:20190410081104p:plain

ローカルトレンドモデルなので一期遅れて値が反映されるような形ですね。

一応ライブラリを用いた場合との比較をしておきます。 pythonだとstatsmodelのUnobservedComponentsでできるようです。

import statsmodels.api as sm

# Fit a local level model
mod_ll = sm.tsa.UnobservedComponents(df[1].values, 'local level')
res_ll = mod_ll.fit()

print(res_ll.summary())
# Show a plot of the estimated level and trend component series
fig_ll = res_ll.plot_components(legend_loc="upper left", figsize=(12,8))
plt.tight_layout()

この結果のパラメータは以下のように最尤推定で算出した値とほぼ同じ値となったので、正しく推定できているようです。 sigmaW,sigmaV=[1636.5036, 594.1718]