状态空间模型 - 从似然函数中集中尺度¶
[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\) 是模型的“尺度”（通常，它取代了一个误差方差项），并且可以通过对似然函数进行微分，以解析方式找到 \(\sigma_*^2\) 的最大似然估计。这意味着仅需要使用数值方法估计参数 \(\psi_*\)，其维度比 \(\psi\) 的维度少一。
示例：局部水平模型¶
(例如，参见 Harvey (1989) 的第 4.2 节)
作为一个具体的例子，考虑局部水平模型，它可以写成
在这个模型中，\(Z, T,\) 和 \(R\) 都固定为等于 \(1\)，有两个未知参数，因此 \(\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)\)。我们可以使用内置的 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).
我们可以查看结果属性 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}
集中尺度¶
现在，有两种方法可以将此模型重新参数化为上述形式。
第一种方法是设置 \(\sigma_*^2 \equiv \sigma_\varepsilon^2\)，这样 \(\psi_* = \psi / \sigma_\varepsilon^2 = [1, q_\eta]\)，其中 \(q_\eta = \sigma_\eta^2 / \sigma_\varepsilon^2\)。
第二种方法是设置 \(\sigma_*^2 \equiv \sigma_\eta^2\)，这样 \(\psi_* = \psi / \sigma_\eta^2 = [h, 1]\)，其中 \(h = \sigma_\varepsilon^2 / \sigma_\eta^2\)。
在第一种情况下，我们只需要针对 \(q_\eta\) 对似然函数进行数值最大化，在第二种情况下，我们只需要针对 \(h\) 对似然函数进行数值最大化。
在大多数情况下，两种方法都可以很好地工作，在下面的示例中，我们将使用第二种方法。
为了重新制定模型以利用集中似然函数，我们需要用参数向量 \(\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)，而尺度的估计值在上面的表中提供。下面，我们将证明这些估计与先前方法的估计一致。
我们还可以查看结果属性 mle_retvals 中的数值优化器结果。事实证明，在这种情况下，只需要进行少两次迭代，因为需要选择的参数少了一个。此外，由于数值最大化问题更易于解决，优化器能够找到一个值，使该参数的梯度比上面更接近于零。
[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)
这两种方法产生了几乎相同的对数似然函数和参数，尽管集中尺度模型能够非常轻微地改善拟合度。
[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