局部模型未完全优化Irregu

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()

它产生以下收敛参数：

^{pr2}$

但我们将相同的结果与R的^{{cd1>}包进行了比较，得到了水平和β值基本相同，而不规则则与方差^{cd2>}的差异（标准差为0.38）。

然后我们决定计算残差的对数可能性，同时忽略水平分量（因为其方差相对较低），以查看对于不规则的最合适值，如：

^{pr3}$

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

^{pr4}$

结果是，由于这一差异，我们正在观察两个包的冲突结果，因此我想寻求帮助，以了解^{{cd3>}组件产生的值似乎无法完全优化拟合数据的对数可能性。

从表面上看，statsmodels通过最大化观测数据的loglikehood来找到状态空间的最佳值；这样，最终状态不应该收敛到标准偏差0.38，因为它进一步提高了目标成本吗？

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

我们用所有稳定版本的statsmodels进行了测试，但得到了相同的^{cd3>}sd结果。

Answer 1

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

我这样说的原因是，给定来自bsts输出的sigma.obs、sigma.level、和{}，对于它们的MCMC算法的任何迭代，我可以复制它们的一步预测误差以及对数似然通过将Kalman滤波器应用到局部水平+exog模型中，使用该迭代的参数值。在

它们与Statsmodels有一个区别，那就是它们根据对数据集的第一次观测和数据集的方差为未观测状态设置先验值，这可能不是最优的（但不应该引起任何重大问题）。Statsmodels反而使用了漫反射先验，这也不会导致任何重大差异。如上所述，当我使用它们的prior时，我可以匹配它们的过滤输出。在

所以区别肯定在估计方法上，这可能与他们的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\β+\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.unregular，因此适当的方差为sigma2.level+sigma2.unregular=0.00403。如果您这样做：

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

好的，这仍然比你最初的post计算的要小，但是这是因为状态不是一个常数，所以{}不包含对数似然计算的实际相关残差。实际的对数似然计算是使用fm.forecasts_error完成的，其方差来自fm.forecasts_error_cov。在