python - 使用未观察到的组件模型模拟时间序列

Question

在使用 statsmodels 中的 UnobservedComponents 拟合本地级别模型后，我们正在尝试寻找使用结果模拟新时间序列的方法。像这样的东西:

import numpy as np
import statsmodels as sm
from statsmodels.tsa.statespace.structural import UnobservedComponents

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

X = 100 + arma_process.generate_sample(nsample=100)
y = 1.2 * x + np.random.normal(size=100)
y[70:] += 10

plt.plot(X, label='X')
plt.plot(y, label='y')
plt.axvline(69, linestyle='--', color='k')
plt.legend();

ss = {}
ss["endog"] = y[:70]
ss["level"] = "llevel"
ss["exog"] = X[:70]

model = UnobservedComponents(**ss)
trained_model = model.fit()

在给定外生变量 X[70:] 的情况下，是否可以使用 trained_model 来模拟新的时间序列？就像我们有 arma_process.generate_sample(nsample=100) 一样，我们想知道我们是否可以做类似的事情:

trained_model.generate_random_series(nsample=100, exog=X[70:])

其背后的动机是，我们可以计算时间序列与观察到的 y[70:] 一样极端的概率(用于识别响应的 p 值大于预测值一个)。

[编辑]

阅读 Josef 和 cfulton 的评论后，我尝试执行以下操作:

mod1 = UnobservedComponents(np.zeros(y_post), 'llevel', exog=X_post)
mod1.simulate(f_model.params, len(X_post))

但这导致模拟似乎无法跟踪 X_post 预测的 predicted_mean 作为 exog。这是一个例子:

虽然 y_post 在 100 左右徘徊，但模拟结果为 -400。这种方法总是导致 p_value 为 50%。

所以当我尝试使用 initial_sate=0 和随机冲击时，结果如下:

现在看来，模拟遵循了预测平均值及其 95% 可信区间(正如 cfulton 在下面评论的那样，这实际上是一种错误的方法，而且它正在取代训练模型的水平方差)。

我尝试使用这种方法只是为了看看我会观察到什么 p 值。以下是我计算 p 值的方法:

samples = 1000
r = 0
y_post_sum = y_post.sum()
for _ in range(samples):
    sim = mod1.simulate(f_model.params, len(X_post), initial_state=0, state_shocks=np.random.normal(size=len(X_post)))
    r += sim.sum() >= y_post_sum
print(r / samples)

对于上下文，这是 Causal Impact谷歌开发的模型。由于它已在 R 中实现，我们一直在尝试使用 statsmodels 作为处理时间序列的核心来复制 Python 中的实现。

我们已经有了一个非常酷的 WIP implementation但是我们仍然需要 p 值来知道实际上我们什么时候产生了不能仅仅用随机性解释的影响(模拟系列并计算总和超过 y_post.sum() 的方法) > 也在 Google 的模型中实现)。

在我的示例中，我使用了 y[70:] += 10。如果我只添加一个而不是十个，Google 的 p 值计算将返回 0.001(有一个y 中的影响)，而在 Python 的方法中它返回 0.247(无影响)。

只有当我将 +5 添加到 y_post 时，模型返回的 p_value 为 0.02，因为它低于 0.05，我们认为对 y_post 有影响。

我正在使用 python3，statsmodels 版本 0.9.0

[EDIT2]

阅读 cfulton 的评论后，我决定全面调试代码以查看发生了什么。这是我发现的:

当我们创建类型为 UnobservedComponents 的对象时，卡尔曼滤波器的表示最终会启动。默认情况下，它 receives the parameter initial_variance作为 1e6 设置 same property对象的。

当我们运行simulate 方法时，initial_state_cov 值is created使用相同的值:

initial_state_cov = (
        np.eye(self.k_states, dtype=self.ssm.transition.dtype) *
        self.ssm.initial_variance
    )

最后，这个相同的值用于查找 initial_state :

initial_state = np.random.multivariate_normal(
    self._initial_state, self._initial_state_cov)

这导致标准差为 1e6 的正态分布。

然后我尝试运行以下命令:

mod1 = UnobservedComponents(np.zeros(len(X_post)), level='llevel', exog=X_post, initial_variance=1)
sim = mod1.simulate(f_model.params, len(X_post))
plt.plot(sim, label='simul')
plt.plot(y_post, label='y')
plt.legend();
print(sim.sum() > y_post.sum())

结果是:

然后我测试了 p 值，最后在 y_post 中测试了 +1 的变化，模型现在可以正确识别添加的信号。

不过，当我使用与 R 的 Google 包中相同的数据进行测试时，p 值仍然存在偏差。也许需要进一步调整输入以提高其准确性。

Answer 1

@Josef 是正确的，你做了正确的事情:

mod1 = UnobservedComponents(np.zeros(y_post), 'llevel', exog=X_post)
mod1.simulate(f_model.params, len(X_post))

simulate 方法根据问题模型模拟数据，这就是为什么当你有外生变量时不能直接使用 trained_model 来模拟。

But for some reason the simulations always ended up being lower than y_post.

我认为这应该是预料之中的 - 运行您的示例并查看估计的系数，我们得到:

                       coef    std err          z      P>|z|      [0.025      0.975]------------------------------------------------------------------------------------sigma2.irregular     0.9278      0.194      4.794      0.000       0.548       1.307sigma2.level         0.0021      0.008      0.270      0.787      -0.013       0.018beta.x1              1.1882      0.058     20.347      0.000       1.074       1.303

The variance of the level is very small, which means that it is extremely unlikely that the level would shift upwards by nearly 10 percent in a single period, based on the model you specified.

When you used:

mod1.simulate(f_model.params, len(X_post), initial_state=0, state_shocks=np.random.normal(size=len(X_post))

发生的事情是水平项是这里唯一未观察到的状态，并且通过提供您自己的方差等于 1 的冲击，您基本上覆盖了模型实际估计的水平方差。 我不认为将初始状态设置为 0 在这里有很大的影响。(参见编辑)。

你写:

the p-value computation was closer, but still is not correct.

我不确定这意味着什么 - 为什么您会期望模型认为这种跳跃很可能发生？您期望达到什么 p 值？

编辑:

感谢您进一步调查(在编辑 2 中)。首先，我认为你应该做的是:

mod1 = UnobservedComponents(np.zeros(y_post), 'llevel', exog=X_post)
initial_state = np.random.multivariate_normal(
    f_model.predicted_state[..., -1], f_model.predicted_state_cov[..., -1])
mod1.simulate(f_model.params, len(X_post), initial_state=initial_state)

现在，解释:

在 Statsmodels 0.9 中，我们还没有对具有扩散初始化的状态进行精确处理(不过从那时起它就被合并了，这是我在测试之前无法复制你的结果的原因之一您使用 0.9 代码库的示例)。这些“初始扩散”状态没有我们可以解决的长期均值(例如随机游走过程)，而局部级别情况下的状态就是这样的状态。

“近似”漫射初始化涉及将初始状态均值设置为零，将初始状态方差设置为一个大数(如您所见)。

对于模拟，默认情况下，初始状态是从给定的初始状态分布中采样的。由于此模型是使用近似扩散初始化进行初始化的，这就解释了为什么您的过程是围绕某个随机数进行初始化的。

您的解决方案是一个很好的补丁，但它不是最优的，因为它没有将模拟期间的初始状态基于估计模型/数据的最后状态。这些值由 f_model.predicted_state[..., -1] 和 f_model.predicted_state_cov[..., -1] 给出。