状態空間モデル - 尤度関数からスケールを集中化¶
[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')
/tmp/ipykernel_3280/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\)はモデルの「スケール」（通常、誤差分散項の1つを置き換えます）であり、尤度関数を微分することによって\(\sigma_*^2\)の最尤推定値を解析的に見つけることができます。これは、数値的方法は\(\psi\)の次元よりも1次元小さい\(\psi_*\)のパラメータを推定するためにのみ必要であることを意味します。
例：局所レベルモデル¶
(たとえば、Harvey (1989)の4.2節を参照してください)
具体的な例として、次のように記述できる局所レベルモデルを考えてみましょう。
このモデルでは、\(Z, T,\)および\(R\)はすべて\(1\)に固定されており、2つの未知のパラメータがあるため、\(\psi = [\sigma_\varepsilon^2, \sigma_\eta^2]\)となります。
一般的なアプローチ¶
まず、statsmodelsの状態空間ライブラリを使用して、スケールを集中化せずにこのモデルを定義する方法を示します。
[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]
このモデルには、var.level \((\sigma_\eta^2)\)とvar.irregular \((\sigma_\varepsilon^2)\)という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, 03 Oct 2024 AIC 921.263
Time: 15:44:49 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).
数値オプティマイザの結果は、results属性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, 03 Oct 2024 AIC 921.263
Time: 15:44:49 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）に表示され、スケールの推定値は上の表に表示されます。以下では、これらの推定値が以前のアプローチからの推定値と一致していることを示します。
そして、数値オプティマイザの結果は、再びresults属性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
/tmp/ipykernel_3280/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])
/tmp/ipykernel_3280/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])
/tmp/ipykernel_3280/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