R 中的基本时间序列算法和统计假设

2020-04-25 08:00:00 · 飞浪

介绍时间序列算法广泛用于分析和预测基于时间的数据。这些算法建立在底层统计假设之上。在本指南中，您将了解底层统计假设和基本时间序列算法以及如何在R中实现它们。让我们从问题陈述和数据开始

介绍

时间序列算法广泛用于分析和预测基于时间的数据。这些算法建立在底层统计假设之上。在本指南中，您将了解底层统计假设和基本时间序列算法以及如何在 R 中实现它们。

让我们从问题陈述和数据开始。

问题陈述

失业是任何国家面临的重大社会经济和政治问题，管理失业也是任何政府的首要任务。在本指南中，您将预测十二个月内的失业率。

本指南中使用的数据来自https://research.stlouisfed.org/fred2h提供的美国经济时间序列数据。数据包含 574 行和 6 个变量，如下所述：

date：从 1968 年 1 月开始的每个月的第一天
psavert：个人储蓄率
pce：个人消费支出，单位：十亿美元
uempmed：失业持续时间中位数，以周为单位
pop：总人口，以千为单位
失业：失业人数（千人）（因变量）

重点将放在日期和失业变量上，因为感兴趣的领域是单变量时间序列预测。

加载所需的库

下面的代码行加载所需的库和数据。

      library(readr)
library(ggplot2)
library(forecast)
library(fpp2)
library(TTR)
library(dplyr)

dat <- read_csv("timeseriesdata.csv")
glimpse(dat)
    

输出：

      Observations: 564
Variables: 7
$ date     <chr> "01-01-1968", "01-02-1968", "01-03-1968", "01-04-1968", "01-0...
$ pce      <dbl> 531.5, 534.2, 544.9, 544.6, 550.4, 556.8, 563.8, 567.6, 568.8...
$ pop      <int> 199808, 199920, 200056, 200208, 200361, 200536, 200706, 20089...
$ psavert  <dbl> 11.7, 12.2, 11.6, 12.2, 12.0, 11.6, 10.6, 10.4, 10.4, 10.6, 1...
$ uempmed  <dbl> 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4.5, 4.2, 4.6, 4.8, 4.4, 4.4, 4...
$ unemploy <int> 2878, 3001, 2877, 2709, 2740, 2938, 2883, 2768, 2686, 2689, 2...
$ Class    <chr> "Train", "Train", "Train", "Train", "Train", "Train", "Train"...
    

数据分区

为了进行模型验证，创建训练和测试数据集。

      dat_train = subset(dat, Class == 'Train')
dat_test = subset(dat, Class == 'Test')

nrow(dat_train); nrow(dat_test)
    

输出：

      1] 552
[1] 12
    

准备时间序列对象

为了在 R 中运行预测模型，您必须将数据转换为时间序列对象，这在下面的第一行代码中完成。开始和结束参数分别指定第一次和最后一次观察的时间。参数频率指定每单位时间的观察次数。

您还需要一个用于计算平均绝对百分比误差 (MAPE) 的实用函数，该函数将用于评估预测模型的性能。MAPE 值越低，预测模型性能越好。这在第二到第四行代码中完成。

      dat_ts <- ts(dat_train[, 6], start = c(1968, 1), end = c(2013, 12), frequency = 12)

#lines 2 to 4
mape <- function(actual,pred){
  mape <- mean(abs((actual - pred)/actual))*100
  return (mape)
}
    

准备好数据和 MAPE 函数后，您就可以开始学习后续章节中的预测技术了。不过，在开始预测之前，了解时间序列中的白噪声和平稳性这两个重要的统计概念非常重要。

白噪音

白噪声序列是纯随机的时间序列，变量独立且均值为零，呈相同分布。这意味着观测值具有相同的方差，并且没有自相关性。下面的代码行设置了可重复性的种子并生成了序列的图。

      set.seed(100)
wnoise = ts(dat_ts)
autoplot(wnoise)
    

输出：

数据中似乎包含信息，而且它不是纯随机序列。在白噪声序列中，预计自相关为零。要将其可视化，请使用 ggAcf ()函数，如下面的代码所示。

      ggAcf(wnoise) + ggtitle("ACF Plot for white noise")

输出：

如果这是白噪声时间序列，95% 或更多的滞后将位于 ACF 图表的边界之间。如上图蓝色虚线所示。在这种情况下，您可以看到大多数线都在蓝色虚线上方，这表明这不是白噪声时间序列。

Ljung-Box 检验

除了上面的可视化，您还可以使用Ljung-Box 检验来确定该序列是否为白噪声序列。显着性检验（由较小的 p 值表示）确认该序列可能不是白噪声。下面的代码执行此测试。

      Box.test(dat_ts, lag = 50, fitdf = 0, type = "Lj")

输出：

      Box-Ljung test

data:  dat_ts
X-squared = 10051, df = 50, p-value < 2.2e-16
    

p 值小于 0.05，表明该序列不是白噪声序列。这意味着数据中存在可用于预测未来值的信息。

文具系列

流行的时间序列算法之一是自回归综合移动平均法 (ARIMA)，它被定义为平稳序列。平稳序列是指其属性不会随时间而变化的序列。简单来说，序列的水平和方差随时间大致保持不变。您可以使用下面的代码可视化该序列。

      plot.ts(dat_ts)

输出：

上图显示该序列不是平稳的，这是构建 ARIMA 模型所必需的。要使序列平稳，请使用R 中的diff()函数执行统计操作差分。如下所示。

      dat_tsdiff1 <- diff(dat_ts, differences=1)
plot.ts(dat_tsdiff1)
    

输出：

上图显示，一阶差分的时间序列在均值和方差上确实看起来大致平稳。因此，我们似乎有一个 ARIMA(p,1,q) 模型。了解了时间序列的基本统计概念后，您现在将构建一些时间序列预测模型。

简单指数平滑法

在指数平滑法中，预测是使用过去观测值的加权平均值得出的，随着观测值的变旧，权重呈指数衰减。水平的平滑参数值由参数alpha决定。以下代码创建了简单的指数平滑模型并打印摘要。

      se_model <- ses(dat_ts, h = 12)
summary(se_model)
    

输出：

      Forecast method: Simple exponential smoothing

Model Information:
Simple exponential smoothing 

Call:
 ses(y = dat_ts, h = 12) 

  Smoothing parameters:
    alpha = 0.9999 

  Initial states:
    l = 2849.6943 

  sigma:  215.1798

     AIC     AICc      BIC 
9419.182 9419.226 9432.123 

Error measures:
                   ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
Training set 13.63606 214.7896 160.8971 0.1946177 2.146727 0.1687594 0.2017391

Forecasts:
         Point Forecast     Lo 80    Hi 80    Lo 95    Hi 95
Jan 2014       10376.04 10100.280 10651.81 9954.299 10797.79
Feb 2014       10376.04  9986.074 10766.01 9779.637 10972.45
Mar 2014       10376.04  9898.438 10853.65 9645.609 11106.48
Apr 2014       10376.04  9824.557 10927.53 9532.618 11219.47
May 2014       10376.04  9759.466 10992.62 9433.070 11319.02
Jun 2014       10376.04  9700.619 11051.47 9343.071 11409.02
Jul 2014       10376.04  9646.504 11105.58 9260.308 11491.78
Aug 2014       10376.04  9596.134 11155.95 9183.274 11568.81
Sep 2014       10376.04  9548.826 11203.26 9110.923 11641.17
Oct 2014       10376.04  9504.080 11248.01 9042.490 11709.60
Nov 2014       10376.04  9461.521 11290.57 8977.402 11774.69
Dec 2014       10376.04  9420.857 11331.23 8915.212 11836.88
    

上面的输出显示，简单指数平滑对所有预测都具有相同的值。由于 alpha 值接近 1，因此预测更接近最近的观测值。现在，评估测试数据上的模型性能。

下面的代码将模型的输出存储在数据框中，并在测试数据中添加一个新变量simplexp ，其中包含简单指数模型的预测值。最后，使用mape函数生成测试数据的 MAPE 误差，结果为 8.5%。

      df_fc = as.data.frame(se_model)
dat_test$simplexp = df_fc$`Point Forecast`
mape(dat_test$unemploy, dat_test$simplexp)
    

输出：

      1] 8.486869

霍尔特趋势法

这是简单指数平滑法的扩展，在生成预测时考虑了趋势成分。此方法涉及两个平滑方程，一个用于水平，一个用于趋势成分。

下面的代码创建了 holt 的模型并打印摘要。

      holt_model <- holt(dat_ts, h = 12)
summary(holt_model)
    

输出：

      Forecast method: Holt's method

Model Information:
Holt's method 

Call:
 holt(y = dat_ts, h = 12) 

  Smoothing parameters:
    alpha = 0.8743 
    beta  = 0.2251 

  Initial states:
    l = 2941.2241 
    b = -0.9146 

  sigma:  200.7121


     AIC     AICc      BIC 
9344.330 9344.440 9365.898 

Error measures:
                    ME     RMSE      MAE        MPE     MAPE    MASE         ACF1
Traini

_{免责声明：本内容来源于第三方作者授权、网友推荐或互联网整理，旨在为广大用户提供学习与参考之用。所有文本和图片版权归原创网站或作者本人所有，其观点并不代表本站立场。如有任何版权侵犯或转载不当之情况，请与我们取得联系，我们将尽快进行相关处理与修改。感谢您的理解与支持！}

_查看原文

技术指南

R语言

阅读全文

R 中的基本时间序列算法和统计假设

杭州电子商务研究院

5年前 · 面向社会、服务行业、政产学研结合、整合资源、和谐发展

Observations: 564 Variables: 7 $ date <chr> "01-01-1968", "01-02-1968", "01-03-1968", "01-04-1968", "01-0... $ pce <dbl> 531.5, 534.2, 544.9, 544.6, 550.4, 556.8, 563.8, 567.6, 568.8... $ pop <int> 199808, 199920, 200056, 200208, 200361, 200536, 200706, 20089... $ psavert <dbl> 11.7, 12.2, 11.6, 12.2, 12.0, 11.6, 10.6, 10.4, 10.4, 10.6, 1... $ uempmed <dbl> 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4.5, 4.2, 4.6, 4.8, 4.4, 4.4, 4... $ unemploy <int> 2878, 3001, 2877, 2709, 2740, 2938, 2883, 2768, 2686, 2689, 2... $ Class <chr> "Train", "Train", "Train", "Train", "Train", "Train", "Train"...

dat_ts <- ts(dat_train[, 6], start = c(1968, 1), end = c(2013, 12), frequency = 12) #lines 2 to 4 mape <- function(actual,pred){ mape <- mean(abs((actual - pred)/actual))*100 return (mape) }

Forecast method: Simple exponential smoothing Model Information: Simple exponential smoothing Call: ses(y = dat_ts, h = 12) Smoothing parameters: alpha = 0.9999 Initial states: l = 2849.6943 sigma: 215.1798 AIC AICc BIC 9419.182 9419.226 9432.123 Error measures: ME RMSE MAE MPE MAPE MASE ACF1 Training set 13.63606 214.7896 160.8971 0.1946177 2.146727 0.1687594 0.2017391 Forecasts: Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 Jan 2014 10376.04 10100.280 10651.81 9954.299 10797.79 Feb 2014 10376.04 9986.074 10766.01 9779.637 10972.45 Mar 2014 10376.04 9898.438 10853.65 9645.609 11106.48 Apr 2014 10376.04 9824.557 10927.53 9532.618 11219.47 May 2014 10376.04 9759.466 10992.62 9433.070 11319.02 Jun 2014 10376.04 9700.619 11051.47 9343.071 11409.02 Jul 2014 10376.04 9646.504 11105.58 9260.308 11491.78 Aug 2014 10376.04 9596.134 11155.95 9183.274 11568.81 Sep 2014 10376.04 9548.826 11203.26 9110.923 11641.17 Oct 2014 10376.04 9504.080 11248.01 9042.490 11709.60 Nov 2014 10376.04 9461.521 11290.57 8977.402 11774.69 Dec 2014 10376.04 9420.857 11331.23 8915.212 11836.88

Forecast method: Holt's method Model Information: Holt's method Call: holt(y = dat_ts, h = 12) Smoothing parameters: alpha = 0.8743 beta = 0.2251 Initial states: l = 2941.2241 b = -0.9146 sigma: 200.7121 AIC AICc BIC 9344.330 9344.440 9365.898 Error measures: ME RMSE MAE MPE MAPE MASE ACF1 Traini