首页 > 解决方案 > 模拟置信区间不等于 conf_int 结果

问题描述

鉴于此模拟数据:

import numpy as np
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.statespace.structural import UnobservedComponents


np.random.seed(12345)
ar = np.r_[1, 0.9]
ma = np.array([1])
arma_process = ArmaProcess(ar, ma)

X = 100 + arma_process.generate_sample(nsample=100)
y = 1.2 * X + np.random.normal(size=100)

我们用前 70 个点构建一个UnobservedComponents模型,对后 30 个点进行推理,如下所示:

model = UnobservedComponents(y[:70], level='llevel', exog=X[:70])
f_model = model.fit()

forecaster = f_model.get_forecast(
    steps=30,
    exog=X[70:].reshape(-1, 1)
)
conf_int = forecaster.conf_int()

如果我们观察 95% 置信区间的平均值,我们会得到以下结果:

conf_int.mean(axis=0)
array([118.19789195, 122.14101161])

但是当试图通过模型模拟获得相同的值时,我们并不能完全得到相同的结果。这是我们为模拟边界运行的脚本:

sim_model = UnobservedComponents(np.zeros(30), level='llevel', exog=X[70:])
res = []
predicted_state = f_model.predicted_state[..., -1]
predicted_state_cov = f_model.predicted_state_cov[..., -1]   
for i in range(1000):
     init_state = np.random.multivariate_normal(
         predicted_state,
         predicted_state_cov
     )
     sim = sim_model.simulate(
         f_model.params,
         30,
         initial_state=init_state)
     res.append(sim.mean())

打印我们得到的下 2.5 和上 97.5 百分位数:

np.percentile(res, [2.5, 97.5])
array([119.06735028, 121.26810407])

当我们使用模型模拟来区分数据中的信号和噪声时,这种差异最终大到足以导致相互矛盾的结论。如果我们做例如:

y[70:] += 1

然后根据第一种技术,我们得出结论,新y的没有信号,因为它的平均值低于122.14。但如果我们使用第二种技术,情况就不一样了:由于上边界是121.2,我们得出结论,有信号。

我们现在试图了解的是这是否是预期的。两种技术的上下 95% 置信区间不应该相等吗?

标签: pythonscipystatsmodels

解决方案


推荐阅读