状態空間モデル - 尤度関数からスケールを集約する¶
[1]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
dta = sm.datasets.macrodata.load_pandas().data
dta.index = pd.date_range(start='1959Q1', end='2009Q4', freq='Q')
C:\Users\user\AppData\Local\Temp\ipykernel_22348\2378637424.py:6: FutureWarning: 'Q' is deprecated and will be removed in a future version, please use 'QE' instead.
dta.index = pd.date_range(start='1959Q1', end='2009Q4', freq='Q')
はじめに¶
（この内容の多くはHarvey (1989)に基づいています；特にセクション3.4を参照）
状態空間モデルは一般的に次のように書くことができます（ここでは時間不変の状態空間モデルに焦点を当てていますが、同様の結果は時間変動モデルにも適用されます）：
しばしば、行列 \(Z, H, T, R, Q\) の一部またはすべての値は未知であり、推定する必要があります。statsmodelsでは、推定は通常、最尤関数を最大化するパラメータを見つけることで行います。特に、パラメータをベクトル \(\psi\) に集めると、これらの行列のそれぞれはそのパラメータの関数として考えることができ、例えば \(Z = Z(\psi)\) などです。
通常、最尤関数は数値的に最大化され、例えば準ニュートン法による「ヒルクライミング」アルゴリズムを適用することで行われますが、パラメータが多くなるほど、この作業は難しくなります。多くのケースで、モデルを \([\psi_*', \sigma_*^2]'\) のように再パラメータ化することで、\(\sigma_*^2\) がモデルの「スケール」（通常は誤差分散項の一つを置き換える）であることがわかり、最尤推定値 \(\sigma_*^2\) を解析的に求めることができます。これにより、数値的手法はパラメータ \(\psi_*\) の推定にのみ必要となり、\(\psi\) の次元より1つ少ない次元になります。
例: ローカルレベルモデル¶
（例えば、Harvey（1989）のセクション4.2を参照）
具体的な例として、ローカルレベルモデルを考えます。このモデルは次のように書けます：
このモデルでは、\(Z, T, R\) はすべて \(1\) に固定され、2つの未知のパラメータがあるため、\(\psi = [\sigma_\varepsilon^2, \sigma_\eta^2]\) となります。
一般的なアプローチ¶
まず、statsmodelsの状態空間ライブラリを使用して、スケールを集中させずにこのモデルをどのように定義するかを示します。パラメータは、\(\psi = [\sigma_\varepsilon^2, \sigma_\eta^2]\) となります。
[2]:
class LocalLevel(sm.tsa.statespace.MLEModel):
_start_params = [1., 1.]
_param_names = ['var.level', 'var.irregular']
def __init__(self, endog):
super(LocalLevel, self).__init__(endog, k_states=1, initialization='diffuse')
self['design', 0, 0] = 1
self['transition', 0, 0] = 1
self['selection', 0, 0] = 1
def transform_params(self, unconstrained):
return unconstrained**2
def untransform_params(self, unconstrained):
return unconstrained**0.5
def update(self, params, **kwargs):
params = super(LocalLevel, self).update(params, **kwargs)
self['state_cov', 0, 0] = params[0]
self['obs_cov', 0, 0] = params[1]
このモデルには選択しなければならない2つのパラメータがあります：var.level \((\sigma_\eta^2)\) と var.irregular \((\sigma_\varepsilon^2)\) です。これらは、組み込みの fit メソッドを使用して、尤度関数を数値的に最大化することで選択できます。
私たちの例では、消費者物価指数のインフレーションにローカルレベルモデルを適用しています。
[3]:
mod = LocalLevel(dta.infl)
res = mod.fit(disp=False)
print(res.summary())
Statespace Model Results
==============================================================================
Dep. Variable: infl No. Observations: 203
Model: LocalLevel Log Likelihood -457.632
Date: Thu, 28 Nov 2024 AIC 921.263
Time: 23:11:30 BIC 931.203
Sample: 03-31-1959 HQIC 925.285
- 09-30-2009
Covariance Type: opg
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
var.level 0.7447 0.156 4.766 0.000 0.438 1.051
var.irregular 3.3733 0.315 10.715 0.000 2.756 3.990
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 182.26
Prob(Q): 0.99 Prob(JB): 0.00
Heteroskedasticity (H): 1.75 Skew: -1.02
Prob(H) (two-sided): 0.02 Kurtosis: 7.18
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
数値最適化器の結果は、結果属性 mle_retvals で確認できます。
[4]:
print(res.mle_retvals)
{'fopt': np.float64(2.254343511382184), 'gopt': array([-7.10302928e-06, -9.72857350e-06]), 'fcalls': 27, 'warnflag': 0, 'converged': True, 'iterations': 7}
スケールを集約する¶
このモデルを上記のように再パラメータ化する方法は2通りあります：
最初の方法は、\(\sigma_*^2 \equiv \sigma_\varepsilon^2\) と設定して、\(\psi_* = \psi / \sigma_\varepsilon^2 = [1, q_\eta]\) とし、ここで \(q_\eta = \sigma_\eta^2 / \sigma_\varepsilon^2\) です。
2番目の方法は、\(\sigma_*^2 \equiv \sigma_\eta^2\) と設定して、\(\psi_* = \psi / \sigma_\eta^2 = [h, 1]\) とし、ここで \(h = \sigma_\varepsilon^2 / \sigma_\eta^2\) です。
最初のケースでは、\(q_\eta\) に関して尤度を数値的に最大化する必要があり、2番目のケースでは、\(h\) に関して尤度を数値的に最大化する必要があります。
どちらのアプローチもほとんどの場合にうまく機能しますが、以下の例では2番目の方法を使用します。
モデルを集約された尤度関数を活用できるように再定式化するためには、パラメータベクトル \(\psi_* = [g, 1]\) の観点でモデルを記述する必要があります。このパラメータベクトルは \(\sigma_\eta^2 \equiv 1\) を定義するため、self['state_cov', 0, 0] = 1 という新しい行を追加し、唯一の未知のパラメータは \(h\) です。\(h\) はもはや分散ではないため、ここでは ratio.irregular に名前を変更しました。
スケールがカルマンフィルタの再帰から計算できるようにモデルを定式化するために必要な重要な部分は、フラグ self.ssm.filter_concentrated = True を設定することです。
[5]:
class LocalLevelConcentrated(sm.tsa.statespace.MLEModel):
_start_params = [1.]
_param_names = ['ratio.irregular']
def __init__(self, endog):
super(LocalLevelConcentrated, self).__init__(endog, k_states=1, initialization='diffuse')
self['design', 0, 0] = 1
self['transition', 0, 0] = 1
self['selection', 0, 0] = 1
self['state_cov', 0, 0] = 1
self.ssm.filter_concentrated = True
def transform_params(self, unconstrained):
return unconstrained**2
def untransform_params(self, unconstrained):
return unconstrained**0.5
def update(self, params, **kwargs):
params = super(LocalLevelConcentrated, self).update(params, **kwargs)
self['obs_cov', 0, 0] = params[0]
再度、組み込みの fit メソッドを使用して、\(h\) の最尤推定を求めることができます。
[6]:
mod_conc = LocalLevelConcentrated(dta.infl)
res_conc = mod_conc.fit(disp=False)
print(res_conc.summary())
Statespace Model Results
==================================================================================
Dep. Variable: infl No. Observations: 203
Model: LocalLevelConcentrated Log Likelihood -457.632
Date: Thu, 28 Nov 2024 AIC 921.263
Time: 23:11:31 BIC 931.203
Sample: 03-31-1959 HQIC 925.285
- 09-30-2009 Scale 0.745
Covariance Type: opg
===================================================================================
coef std err z P>|z| [0.025 0.975]
-----------------------------------------------------------------------------------
ratio.irregular 4.5297 1.226 3.694 0.000 2.126 6.933
===================================================================================
Ljung-Box (L1) (Q): 0.00 Jarque-Bera (JB): 182.26
Prob(Q): 0.99 Prob(JB): 0.00
Heteroskedasticity (H): 1.75 Skew: -1.02
Prob(H) (two-sided): 0.02 Kurtosis: 7.18
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
\(h\)の推定値はパラメータの中央のテーブル（ratio.irregular）に提供されており、スケールの推定値は上部のテーブルに提供されています。以下では、これらの推定値が前のアプローチからのものと一致していることを示します。
そして、再び数値最適化器の結果を結果属性 mle_retvals から確認できます。この場合、選択するパラメータが1つ少なかったため、2回少ない反復回数で済みました。さらに、数値的な最大化問題が簡単だったため、最適化器はこのパラメータに対する勾配が以前よりもわずかにゼロに近づくような値を見つけることができました。
[7]:
print(res_conc.mle_retvals)
{'fopt': np.float64(2.2543435111703576), 'gopt': array([-6.71906974e-08]), 'fcalls': 12, 'warnflag': 0, 'converged': True, 'iterations': 5}
推定結果の比較¶
\(h = \sigma_\varepsilon^2 / \sigma_\eta^2\) およびスケールは \(\sigma_*^2 = \sigma_\eta^2\) であることを思い出してください。これらの定義を使用すると、両方のモデルがほぼ同じ結果を生み出していることがわかります。
[8]:
print('Original model')
print('var.level = %.5f' % res.params[0])
print('var.irregular = %.5f' % res.params[1])
print('\nConcentrated model')
print('scale = %.5f' % res_conc.scale)
print('h * scale = %.5f' % (res_conc.params[0] * res_conc.scale))
Original model
var.level = 0.74469
var.irregular = 3.37330
Concentrated model
scale = 0.74472
h * scale = 3.37338
C:\Users\user\AppData\Local\Temp\ipykernel_22348\976557282.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
print('var.level = %.5f' % res.params[0])
C:\Users\user\AppData\Local\Temp\ipykernel_22348\976557282.py:3: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
print('var.irregular = %.5f' % res.params[1])
C:\Users\user\AppData\Local\Temp\ipykernel_22348\976557282.py:7: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
print('h * scale = %.5f' % (res_conc.params[0] * res_conc.scale))
例: SARIMAX¶
SARIMAXモデルでは、デフォルトで分散項は尤度関数を数値的に最大化することで選択されますが、スケールを外向きに集中させるオプションが追加されました。
[9]:
# Typical approach
mod_ar = sm.tsa.SARIMAX(dta.cpi, order=(1, 0, 0), trend='ct')
res_ar = mod_ar.fit(disp=False)
# Estimating the model with the scale concentrated out
mod_ar_conc = sm.tsa.SARIMAX(dta.cpi, order=(1, 0, 0), trend='ct', concentrate_scale=True)
res_ar_conc = mod_ar_conc.fit(disp=False)
これらの2つのアプローチは、ほぼ同じ対数尤度とパラメータを生成しますが、集中スケールを使用したモデルはフィットをわずかに改善することができました。
[10]:
print('Loglikelihood')
print('- Original model: %.4f' % res_ar.llf)
print('- Concentrated model: %.4f' % res_ar_conc.llf)
print('\nParameters')
print('- Original model: %.4f, %.4f, %.4f, %.4f' % tuple(res_ar.params))
print('- Concentrated model: %.4f, %.4f, %.4f, %.4f' % (tuple(res_ar_conc.params) + (res_ar_conc.scale,)))
Loglikelihood
- Original model: -245.8275
- Concentrated model: -245.8264
Parameters
- Original model: 0.4921, 0.0243, 0.9808, 0.6490
- Concentrated model: 0.4864, 0.0242, 0.9809, 0.6492
今回は、集中アプローチの下で最適化の反復回数が約1/3少なくて済みます。
[11]:
print('Optimizer iterations')
print('- Original model: %d' % res_ar.mle_retvals['iterations'])
print('- Concentrated model: %d' % res_ar_conc.mle_retvals['iterations'])
Optimizer iterations
- Original model: 36
- Concentrated model: 22