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

在使用来自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值仍然不合适.也许这是进一步调整输入以提高其准确性的问题.

最佳答案

@Josef是正确的,你做对了：

mod1 = UnobservedComponents(np.zeros(y_post),len(X_post))

模拟方法根据所讨论的模型模拟数据,这就是为什么你不能直接使用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.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％.

当你使用：

mod1.simulate(f_model.params,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),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,initial_state=initial_state)

现在,解释：

在Statsmodels 0.9中,我们还没有对弥散初始化的状态进行精确处理(从那时起它就被合并了,但这是我无法复制结果的一个原因,直到我测试你的例子为止0.9代码库).这些“初始扩散”状态不具有我们可以解决的长期均值(例如随机游走过程),并且本地级情况中的状态是这样的状态.

“近似”漫反射初始化包括将初始状态均值设置为零,将初始状态方差设置为大数(如您所发现的).

对于模拟,默认情况下,初始状态是从给定的初始状态分布中采样的.由于此模型使用近似漫反射初始化进行初始化,因此可以解释为什么您的进程是围绕某个随机数进行初始化的.

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

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

相关推荐