局部模型未完全优化不规则状态

Question

我们使用局部模型来拟合此数据：

import pandas as pd
from scipy import stats
from statsmodels.tsa.statespace.structural import UnobservedComponents


data = {
    'Basque': [3.85318463, 3.9456583, 4.03356173, 4.0234219, 4.01378197, 4.2859184, 4.5743361, 4.89895735, 5.19701498, 5.33890298, 5.46515301, 5.54591563, 5.61489573, 5.85218493, 6.08140542, 6.17009424, 6.2836334, 6.5555554, 6.81076856, 7.1051843, 7.37789168],
    'another': [3.54662963, 3.69044557, 3.826835, 3.87567838, 3.92173673, 4.2417882, 4.57533548, 4.83804641, 5.0813341, 5.15809788, 5.22365053, 5.33247651, 5.42944892, 5.67437885, 5.91552394, 6.06683787, 6.22764921, 6.53906013, 6.83797506, 6.98736082, 7.12489303]
}

data = pd.DataFrame(data)

model = UnobservedComponents(data['Basque'], level='llevel', exog=data['another'])
fm = model.fit()

产生以下收敛参数：

sigma2.irregular    2.365834e-11
sigma2.level        4.032277e-03
beta.another        9.704585e-01
dtype: float64

但随后我们将相同的结果与来自 R 的 bsts 包进行了比较，得到了基本相同的水平和 beta 值，而不规则的则完全不同：0.16 方差（0.38 标准差）。

然后我们决定计算残差的对数似然，同时忽略水平分量（因为它的方差相对较低），看看什么是不规则的最合适的值，如下所示：

residuals = data['Basque'] - fm.params[-1] * data['another']

然后我们计算了两个标准差的对数似然值（statsmodels 为 5e-6，R 为 0.38），如下所示：

stats.norm.logpdf(residuals, 0, scale=5e-6).sum() # statsmodels sd
-> -39342125822.669464

stats.norm.logpdf(residuals, 0, scale=0.38).sum() # bsts sd
-> -5.789754559185148

事实证明，正是由于这种差异，我们从两个包中观察到相互矛盾的结果，所以我想寻求帮助以了解为什么 irregular 组件产生的值看起来不像完全优化拟合数据的对数似然。

看起来，statsmodels 通过最大化 loglikehood 观察数据来找到状态空间的最佳值；在这种情况下，最终状态不应该收敛到 0.38 的标准偏差，因为它进一步提高了 objective 成本吗？

作为参考，this issue 在我们基于 statsmodels 开发的库中打开；这是一个与因果影响推断相关的非常有趣的问题，不规则组件让 Python 包确定存在因果影响，而 R 的结论恰恰相反。

我们测试了所有稳定版本的 statsmodels，但得到了相同的 irregular sd 结果。

Answer 1

这是一个有趣的问题和两个包之间的比较。差异显然来自不同的估计方法：Statsmodels 中的最大似然和 bsts 中的 Bayesian MCMC。在这种情况下出现差异并不奇怪，因为时间序列非常短。

我这样说的原因是，给定 bsts 输出的 sigma.obs、sigma.level 和 coefficients，对于他们的 MCMC 算法的任何迭代，我可以复制他们的 one.step.prediction_errors 和 log.likelihood 通过使用该迭代的参数值将卡尔曼滤波器应用于本地级别 + exog 模型来进行该迭代。

它们确实与 Statsmodels 有一个区别，那就是它们根据数据集的第一次观察和数据集的方差为未观察到的状态设置先验，这可能不是最优的（但不应该是导致任何重大问题）。 Statsmodels 改为使用弥散先验，这同样不会造成任何重大差异。正如我上面提到的，当我使用他们的先验时，我可以匹配他们的过滤输出。

所以区别肯定是在估计方法上，这可能与他们的MCMC算法的细节和他们设置的先验有关。您可以跟进他们，看看他们是否对他们的设置可能如何影响结果有任何直觉。

讨论残差和对数似然计算

We decided then to compute the loglikelihood of the residuals while ignoring the level component (as its variance is relatively low) to see what would be the most appropriate value for irregular, like so:
residuals = data['Basque'] - fm.params[-1] * data['another']
And then we computed loglikelihood for both standard deviations (5e-6 for statsmodels, 0.38 for R), like so:
stats.norm.logpdf(residuals, 0, scale=5e-6).sum() # statsmodels sd
-> -39342125822.669464

stats.norm.logpdf(residuals, 0, scale=0.38).sum() # bsts sd
-> -5.789754559185148

问题是不规则分量的方差估计与状态不是常数的模型中残差的方差不同。相反，模型是：

y_t = x_t \beta + \alpha_t + \varepsilon_t

和 Var(y_t - x_t \beta | t-1) = Var(\alpha_t | t-1) + Var(\varepsilon_t | t-1)

因为 Var(\alpha_t | t-1) = sigma2.level 和 Var(\varepsilon_t | t-1) = sigma2.irregular，合适的方差是sigma2.level + sigma2.irregular = 0.00403。如果你这样做：

stats.norm.logpdf(residuals, 0, scale=0.00403).sum() # statsmodels sd
-> -205.4612464569425

好吧，这仍然比 bsts 计算的原始 post 小，但这是因为状态不是常量，所以 residuals 不是包含对数似然计算的实际相关残差。实际的对数似然计算是使用 fm.forecasts_error 和 fm.forecasts_error_cov.

的方差完成的

局部模型未完全优化不规则状态

Local Level Model Not Fully Optimizing Irregular State

python

statsmodels