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

在使用来自statsmodel的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的预测的预测值作为exog.这是一个例子：

虽然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开发的Causal Impact模型.由于它已在R中??实现,我们一直在尝试使用statsmodels作为处理时间序列的核心来复制Python中的实现.

我们已经有了一个非常酷的WIP implementation,但是我们仍然需要知道p值,实际上我们的影响实际上并不仅仅是随机性的解释(模拟系列的方法和计算总和超过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的评论后,我决定完全调试代码,看看发生了什么.这是我发现的：

当我们创建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值仍然不合适.也许这是进一步调整输入以提高其准确性的问题.