前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >用于时间序列数据的泊松回归模型

用于时间序列数据的泊松回归模型

作者头像
deephub
发布2021-04-30 14:59:59
2K0
发布2021-04-30 14:59:59
举报
文章被收录于专栏:DeepHub IMBADeepHub IMBA

泊松和类泊松回归模型常用于基于计数的数据集,即包含整数计数的数据。例如,每小时走进医院急诊室的人数就是一个这样的数据集。

基于普通最小二乘回归的线性模型或非线性模型(例如基于基于神经网络的回归技术的线性模型)不适用于此类数据集,因为它们可以预测负值。

如果数据集是计数的时间序列,则会产生额外的建模复杂性,因为时间序列数据通常是自相关的。以前的计数会影响将来计数的值。如果回归模型不能充分捕获这些相关性中包含的“信息”,“未解释的”信息将以自相关误差的形式泄漏到模型的残差中。在这种情况下,模型的拟合优度会很差。

解决这个问题的一般补救办法如下:

  1. 在拟合回归模型之前,检查时间序列是否具有季节性,如果具有,则进行季节性调整。这样做,就算有季节性的自相关性,也可以解释得通。
  2. 对所有t进行时间序列的第一次差分,即y_t - y_(t-1),并对差分时间序列进行白噪声测试。如果差分时间序列是白噪声,则原始时间序列是随机游走。在这种情况下,不需要进一步建模。
  3. 在季节性调整后的时间序列上拟合基于Poisson(或相关)计数的回归模型,但包括因变量y的滞后副本作为回归变量。

在本文中,我们将解释如何使用方法(3)在计数的时间序列上拟合泊松或类泊松模型。

MANUFACTURING STRIKES数据集

为了说明模型拟合过程,我们将使用以下在回归建模文献中广泛使用的开源数据集:

该数据是一个月度时间序列,显示了从1968年到1976年,美国制造业活动(偏离趋势线)与美国制造业合同罢工数量之间的关系。

这个数据可以使用statsmodels数据集包来获取。

因变量y是罢工数量。

我们将从导入所有必需的程序包开始:

代码语言:javascript
复制
 import statsmodels.api as sm
 import statsmodels.discrete.discrete_model as dm
 import numpy as np
 from patsy import dmatrices
 import statsmodels.graphics.tsaplots as tsa
 from matplotlib import pyplot as plt

让我们使用statsmodels将数据集加载到内存中:

代码语言:javascript
复制
 strikes_dataset = sm.datasets.get_rdataset(dataname='StrikeNb', package='Ecdat')

打印出数据集:

代码语言:javascript
复制
 print(strikes_dataset.data)

我们看到以下输出:

我们将前92个数据点视为训练集,其余16个数据点视为测试数据集

代码语言:javascript
复制
 strikes_data = strikes_dataset.data.copy()
 strikes_data_train = strikes_data.query('time<=92')
 strikes_data_test = strikes_data.query('time>92').reset_index().drop('index', axis=1)

打印出因变量的统计信息:

代码语言:javascript
复制
 print('Mean='+str(np.mean(strikes_data_train['strikes'])) + ' Variance='+str(np.var(strikes_data_train['strikes'])))

我们得到以下输出:

代码语言:javascript
复制
 Mean=5.5 Variance=14.728260869565217

我们看到y过于分散,因此它违反了Poisson模型的均值=方差假设。为了解决过度分散问题,我们将拟合具有以下NB2方差函数的负二项式回归模型:

这是我们的回归表达式。罢工是因变量,输出是我们的解释变量。假定存在回归截距:

代码语言:javascript
复制
 expr = 'strikes ~ output

我们将使用Patsy来计算X和y矩阵。Patsy会自动向X添加回归截距列:

代码语言:javascript
复制
 y_train, X_train = dmatrices(expr, strikes_data_train, return_type='dataframe')
 print(y_train)
 print(X_train)
 
 y_test, X_test = dmatrices(expr, strikes_data_test, return_type='dataframe')
 print(y_test)
 print(X_test)

使用NB2方差函数构建和训练负二项式回归模型:

代码语言:javascript
复制
 nb2_model = dm.NegativeBinomial(endog=y_train, exog=X_train, loglike_method='nb2')
 nb2_model_results = nb2_model.fit(maxiter=100)
 print(nb2_model_results.summary())

我们得到了拟合模型摘要的以下输出:

如系数的p值(0.034和0.000)所示,输出和色散参数alpha在95%置信度上都具有统计学意义。

贴合度

伪R平方仅为0.9%,说明训练数据集的拟合质量非常差。

对数似然比检验的p值为0.03589,表明该模型在95%的置信水平下比仅截距模型(又称为零模型)做得更好,但在99%或更高的置信水平下表现不佳。

让我们看一下拟合模型的残差的自相关图:

我们可以看到残差误差在时间滞后1、2和3时是自相关的,这表明因变量罢工中存在自相关,因为NB2模型无法完全解释导致泄漏到模型残差中的原因。。

总体而言,此模型的拟合优度非常差。

建立自回归泊松模型

为了解决残差自相关的情况,我们将引入y的滞后副本,具体为y(t-1)、y(t-2)和y(t-3)作为输出变量的回归变量。但我们不直接引入y(t-k)作为回归变量,而是使用ln[y(t-k)]来解决y(t-k)系数为正时的“模型爆炸”问题。

但是使用ln()变换产生了一个问题即如何处理0的y_t对于这个对数是没有定义的。所以我们使用Cameron和Trivedi在他们的书《Regression Analysis of Count Data》(见第7.5节:自回归模型)中概述的以下技巧来解决这个问题:

我们将为每个感兴趣的时延定义一个新的指标变量d_t,如下所示:

  • 当y_t = 0时:d_t=1,我们将y_t设为0。
  • 当y_t > 0: d_t=0时,我们将y_t保持原样。

让我们对dataframe进行这些更改。创建一个新列strikes_adj,如果strikes < 1,则设置为1,否则设置为strikes:

代码语言:javascript
复制
 strikes_data['strikes_adj'] = np.maximum(1, strikes_data['strikes'])

定义一个函数来设置指标变量d_t的值,如上所述。

代码语言:javascript
复制
 def indicator_func(x):
     if x == 0:
         return 1
     else:
         return 0

并使用此函数创建一个新的指标变量列:

代码语言:javascript
复制
 strikes_data['d'] = strikes_data['strikes'].apply(indicator_func)

现在为strikes_adj和d创建滞后变量:

代码语言:javascript
复制
 strikes_data['ln_strikes_adj_lag1'] = strikes_data['strikes_adj'].shift(1)
 strikes_data['ln_strikes_adj_lag2'] = strikes_data['strikes_adj'].shift(2)
 strikes_data['ln_strikes_adj_lag3'] = strikes_data['strikes_adj'].shift(3)
 
 strikes_data['d_lag1'] = strikes_data['d'].shift(1)
 strikes_data['d_lag2'] = strikes_data['d'].shift(2)
 strikes_data['d_lag3'] = strikes_data['d'].shift(3)

删除所有带有空单元格的行:

代码语言:javascript
复制
 strikes_data = strikes_data.dropna()

最后,取一个自然日志ln_strikes_adj_lag1。回想一下,我们想添加滞后变量y(t_1),y(t_2)和y_(t_3)的自然对数。

让我们看看我们的数据现在是什么样子:

代码语言:javascript
复制
 print(strikes_data)

让我们再次将数据分为训练和测试数据集:

代码语言:javascript
复制
 strikes_data_train=strikes_data.query('time<=92')
 strikes_data_test=strikes_data.query('time>92').reset_index().drop('index', axis=1)

我们的回归表达式也需要更新以包含滞后变量:

代码语言:javascript
复制
 expr = 'strikes ~ output + ln_strikes_adj_lag1 + ln_strikes_adj_lag2 + ln_strikes_adj_lag3 + d_lag1 + d_lag2 + d_lag3'

使用Patsy可以得出y和X矩阵:

代码语言:javascript
复制
 y_train, X_train = dmatrices(expr, strikes_data_train, return_type='dataframe')
 print(y_train)
 print(X_train)

最后,我们将在(y_train,X_train)上构建并拟合回归模型。这次,我们将使用直接的Poisson回归模型:

代码语言:javascript
复制
 poisson_model = dm.Poisson(endog=y_train, exog=X_train)
 poisson_model_results = poisson_model.fit(maxiter=100)
 print(poisson_model_results.summary()

我们看到以下结果:

贴合度

首先要注意的是,通过Pseudo-R-squared测量的拟合优度比早期的NB2模型提高了0.9%至15.69%。这是一个很大的进步。这次,LLR测试的p值也很小,为1.295e-15。这意味着我们可以以接近100%的置信度说滞后变量Poisson模型要比截距模型好。回想一下,我们可以说只有95%的置信度。

让我们看一下这个滞后变量泊松模型的残差的自相关图:

代码语言:javascript
复制
 tsa.plot_acf(poisson_model_results.resid, alpha=0.05)
 plt.show()

我们看到以下图:

除了在LAG 13处存在非常轻微的显着相关性之外,残差与所有其他滞后的相关性都很好地位于规定的alpha范围内。

我们在Poisson模型的回归变量中添加滞后罢工副本的策略似乎已经解释了很多罢工变量中的自相关。

变量的意义

我们从训练注意总结滞后变量的泊松模型的系数,输出,ln_strikes_adj_lag1和ln_strikes_adj_lag2是重要的在95%置信水平,第三个滞后ln_strikes_adj_lag3是重要的系数只有大约75%置信水平表示的假定值为0.237。此外,三个滞后指标变量d_lag1、d_lag2和d_lag3在95%置信水平下均不具有统计学意义。

预测

让我们使用拟合的滞后变量Poisson模型来预测我们先前预留的测试数据集的罢工次数。我们不应该寄希望于预测的质量过高。请记住,尽管此模型比以前的NB2模型拟合得好得多,但伪R平方仍然仅为16%。

我们将使用Patsy来查看(y_test, X_test):

代码语言:javascript
复制
 y_test, X_test = dmatrices(expr, strikes_data_test, return_type='dataframe')
 print(y_test)
 print(X_test)

对X_test进行预测:

代码语言:javascript
复制
 poisson_predictions = poisson_model_results.predict(X_test

绘制预测值和实际值:

代码语言:javascript
复制
 predicted_counts=poisson_predictions
 actual_counts = y_test['strikes']
 fig = plt.figure()
 fig.suptitle('Predicted versus actual strike counts')
 predicted, = plt.plot(X_test.index, predicted_counts, 'go-', label='Predicted counts')
 actual, = plt.plot(X_test.index, actual_counts, 'ro-', label='Actual counts')
 plt.legend(handles=[predicted, actual])
 plt.show()

我们得到如下图:

下一步

我们可以尝试通过以下修改来改善滞后变量模型的拟合优度:

除输出外,还将输出变量的前三个时间滞后量作为回归变量。

将输出变量和罢工变量的时滞值作为回归变量。

使用负二项模型(使用NB1或NB2方差函数)代替泊松模型,并将上述类型的滞后变量作为回归变量。

论文和相关连接

Cameron A. Colin, Trivedi Pravin K., Regression Analysis of Count Data, Econometric Society Monograph №30, Cambridge University Press, 1998. ISBN: 0521635675

McCullagh P., Nelder John A., Generalized Linear Models, 2nd Ed., CRC Press, 1989, ISBN 0412317605, 9780412317606

Kennan J., The duration of contract strikes in U.S. manufacturing, Journal of Econometrics, Volume 28, Issue 1, 1985, Pages 5–28, ISSN 0304–4076, https://doi.org/10.1016/0304-4076(85)90064-8. PDF download link

Cameron C. A., Trivedi P. K., Regression-based tests for overdispersion in the Poisson model, Journal of Econometrics, Volume 46, Issue 3, 1990, Pages 347–364, ISSN 0304–4076, https://doi.org/10.1016/0304-4076(90)90014-K.

作者:Sachin Date

原文地址:https://towardsdatascience.com/poisson-regression-models-for-time-series-data-sets-54114e68c46d

deephub翻译组

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2021-04-28,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 DeepHub IMBA 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • MANUFACTURING STRIKES数据集
  • 建立自回归泊松模型
  • 预测
  • 下一步
  • 论文和相关连接
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档