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

Question

在使用UnobservedComponentsfrom 拟合本地级模型之后statsmodels,我们试图找到用结果模拟新时间序列的方法.就像是:

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

在trained_model给定外生变量的情况下,是否可以用来模拟新的时间序列X[70:]？正如我们所拥有的那样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))

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

当y_post曲折在100左右时,模拟在-400.这种方法总是导致p_value为50%.

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

现在似乎模拟遵循预测的平均值和95%的可信区间(如下面的评论,这实际上是一种错误的方法,它取代了训练模型的水平方差).

我尝试使用这种方法只是为了看看我观察到的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开发的因果影响模型.由于它已在R中​​实现,我们一直在尝试使用Python statsmodels作为处理时间序列的核心来复制Python中的实现.

我们已经有了一个非常酷的WIP 实现,但是我们仍然需要知道p值,实际上我们的影响实际上是无法用随机性来解释的(模拟系列的方法和计算总和超过的y_post.sum()方法也是如此)在谷歌的模型中).

在我的例子中,我使用了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,最终会启动卡尔曼滤波器的表示.默认情况下,它接收参数initial_variance 1e6,它设置对象的相同属性.

当我们运行simulate方法,该initial_state_cov值将创建使用相同的值:

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 std err z P> | z | [0.025 0.975]
-------------------------------------------------- ----------------------------------
sigma2.irregular 0.9278 0.194 4.794 0.000 0.548 1.307
sigma2水平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％，根据您指定的型号。

使用时：

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

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

你写：

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

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

编辑：

感谢您进一步调查（在Edit 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]。