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

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 的 cmets 后，我尝试实现以下内容：

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

但这导致模拟似乎无法将 predicted_mean 的预测跟踪 X_post 作为 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)

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

我们已经有了一个非常酷的 WIP implementation，但我们仍然需要 p 值来知道实际上我们何时产生了不能由单纯随机性解释的影响（模拟系列并计算那些summation超过y_post.sum()也在Google的模型中实现）。

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

只有当我给y_post加上+5时，模型返回的p_value为0.02并且低于0.05，我们认为y_post有影响。

我正在使用 python3，statsmodels 版本 0.9.0

[EDIT2]

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

当我们创建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模拟。

但由于某种原因，模拟结果总是低于 y_post。

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

coef 标准错误 z P>|z| [0.025 0.975] -------------------------------------------------- ---------------------------------- sigma2.不规则 0.9278 0.194 4.794 0.000 0.548 1.307 sigma2.level 0.0021 0.008 0.270 0.787 -0.013 0.018 Beta.x1 1.1882 0.058 20.347 0.000 1.074 1.303

水平的方差非常小，这意味着根据您指定的模型，水平在单个时期内向上移动近 10% 的可能性极。 p>

您使用时：

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

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

你写：

p 值计算更接近，但仍然不正确。

我不确定这意味着什么 - 为什么您希望模型认为这种跳跃很可能发生？您期望达到什么 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] 给出。