使用 R 进行时间序列预测

2019-07-12 08:00:00 · 飞浪

介绍在本指南中，您将学习如何使用统计编程语言“R”实现以下时间序列预测技术：朴素方法简单指数平滑法霍尔特趋势法 ARIMA 丁基对苯二甲酸我们将从探索数据开始。问

介绍

在本指南中，您将学习如何使用统计编程语言“R”实现以下时间序列预测技术：

朴素方法
简单指数平滑法
霍尔特趋势法
ARIMA
丁基对苯二甲酸

我们将从探索数据开始。

问题陈述

失业是任何国家的主要社会经济和政治问题，因此，管理失业是任何政府的主要任务。在本指南中，我们将尝试预测十二个月内的失业率。本指南中使用的数据来自美联储经济数据提供的美国经济时间序列数据。数据包含 574 行和 6 个变量，如下所述：

日期 - 表示从 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...
$ pce  	<dbl> 531.5, 534.2, 544.9, 544.6, 550.4, 556.8, 563.8, 567....
$ pop  	<int> 199808, 199920, 200056, 200208, 200361, 200536, 20070...
$ psavert  <dbl> 11.7, 12.2, 11.6, 12.2, 12.0, 11.6, 10.6, 10.4, 10.4,...
$ uempmed  <dbl> 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4.5, 4.2, 4.6, 4.8, 4.4...
$ unemploy <int> 2878, 3001, 2877, 2709, 2740, 2938, 2883, 2768, 2686,...
$ Class	<chr> "Train", "Train", "Train", "Train", "Train", "Train",...
    

数据分区

我们将创建用于模型验证的训练和测试数据集。下面的前几行代码使用已经具有 Train 和 Test 标签的“Class”变量执行此任务。

第三行代码打印训练和测试数据的行数。请注意，测试数据集有 12 行，因为我们将构建模型来预测未来 12 个月的情况。

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

输出：

      1] 552
[1] 12
    

准备时间序列对象

要在“R”中运行预测模型，我们需要将数据转换为时间序列对象，这在下面的第一行代码中完成。“start”和“end”参数分别指定第一次和最后一次观察的时间。参数“frequency”指定每单位时间的观察次数。

我们还将创建一个用于计算平均绝对百分比误差（或 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 函数后，我们将转到后续章节中的预测技术。

朴素预测方法

最简单的预测方法是使用最近的观察结果作为下一次观察结果的预测。这称为朴素预测，可以使用“naive()”函数实现。这种方法可能不是最好的预测技术，但它通常为其他更先进的预测方法提供了有用的基准。

下面的第一行代码读取时间序列对象“dat_ts”并创建简单预测模型。第二个参数“h”指定要预测的值的数量，在我们的例子中设置为 12。第二行打印模型摘要以及未来 12 个月的预测值。

      naive_mod <- naive(dat_ts, h = 12)
summary(naive_mod)
    

输出：

      Forecast method: Naive method
 
Model Information:
Call: naive(y = dat_ts, h = 12)
 
Residual sd: 214.9766
 
Error measures:
                	ME 	RMSE      MAE   	MPE 	MAPE      MASE
 Training set 13.60799 214.9766 161.1361 0.1931667 2.148826 0.1690101
              	ACF1
 Training set 0.201558
 
Forecasts:
     	Point Forecast 	Lo 80    Hi 80	Lo 95	Hi 95
Jan 2014      	10376 10100.496 10651.50 9954.654 10797.35
Feb 2014      	10376  9986.379 10765.62 9780.126 10971.87
Mar 2014      	10376  9898.814 10853.19 9646.206 11105.79
Apr 2014      	10376  9824.993 10927.01 9533.307 11218.69
May 2014      	10376  9759.955 10992.04 9433.841 11318.16
Jun 2014      	10376  9701.157 11050.84 9343.916 11408.08
Jul 2014      	10376  9647.086 11104.91 9261.222 11490.78
Aug 2014      	10376  9596.758 11155.24 9184.252 11567.75
Sep 2014      	10376  9549.489 11202.51 9111.961 11640.04
Oct 2014      	10376  9504.781 11247.22 9043.585 11708.41
Nov 2014      	10376  9462.258 11289.74 8978.552 11773.45
Dec 2014      	10376  9421.627 11330.37 8916.413 11835.59
    

上面的输出表明，朴素方法在整个预测范围内预测相同的值。现在让我们使用预测值并评估测试数据上的模型性能。

下面的第一行代码在测试数据中添加了一个新变量 naive，其中包含从 naive 方法获得的预测值。第二行使用 mape 函数在测试数据上生成 MAPE 误差，结果为 8.5%。

      dat_test$naive = 10376
mape(dat_test$unemploy, dat_test$naive)  ## 8.5%
    

输出：

      1] 8.486551

简单指数平滑法

指数平滑法是简单法的延伸，其中预测是使用过去观测值的加权平均值得出的，权重随着观测值的变旧而呈指数衰减。简而言之，最近的观测值赋予更高的权重，反之亦然。水平的平滑参数值由参数“alpha”决定。

下面的第一行代码读取时间序列对象“dat_ts”并创建简单指数平滑模型。第二行打印模型摘要以及未来 12 个月的预测值。

      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
 Training set 13.63606 214.7896 160.8971 0.1946177 2.146727 0.1687594
               	ACF1
 Training set 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_model”中。第二行打印摘要和未来 12 个月的预测。

      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
Training set -1.775468 199.9836 152.9176 0.04058417 2.056684 0.16039
                	ACF1
Training set -0.001011626
 
Forecasts:
    	 Point Forecast	Lo 80     Hi 80	Lo 95	Hi 95
Jan 2014  	10194.656 9937.433 10451.879 9801.267 10588.04
Feb 2014   	9973.102 9590.814 10355.390 9388.443 10557.76
Mar 2014   	9751.549 9239.463 10263.634 8968.381 10534.72
Apr 2014   	9529.995 8881.048 10178.943 8537.515 10522.48
May 2014   	9308.442 8514.997 10101.887 8094.972 10521.91
Jun 2014   	9086.888 8141.265 10032.512 7640.682 10533.10
Jul 2014   	8865.335 7759.990  9970.680 7174.856 10555.81
Aug 2014   	8643.782 7371.378  9916.185 6697.808 10589.76
Sep 2014   	8422.228 6975.652  9868.804 6209.881 10634.58
Oct 2014   	8200.675 6573.036  9828.313 5711.416 10689.93
Nov 2014   	7979.121 6163.744  9794.498 5202.741 10755.50
Dec 2014   	7757.568 5747.979  9767.157 4684.167 10830.97
    

上面的输出显示训练数据的 MAPE 为 2.1%。现在让我们评估测试数据的模型性能，这在下面的代码行中完成。测试数据的 MAPE 误差为 6.6%，这比以前的模型有所改进。

      df_holt = as.data.frame(holt_model)
dat_test$holt = df_holt$`Point Forecast`
mape(dat_test$unemploy, dat_test$holt)
    

输出：

      1] 6.597904

ARIMA

ARIMA 建模是时间序列预测最流行的方法之一。指数平滑模型基于对数据趋势和季节性的描述，而 ARIMA 模型旨在描述数据中的自相关性。

“R” 中的“auto.arima()”函数用于使用 Hyndman-Khandakar 算法的变体构建 ARIMA 模型，该算法结合单位根检验、AICc 最小化和 MLE 以获得 ARIMA 模型。

下面的第一行代码创建了 ARIMA 模型并将其存储在对象“arima_model”中。第二行打印摘要和未来 12 个月的预测。

      <code class="

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

_查看原文

技术指南

R语言

阅读全文

使用 R 进行时间序列预测

杭州电子商务研究院

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

Observations: 564 Variables: 7 $ date <chr> "01-01-1968", "01-02-1968", "01-03-1968", "01-04-1968... $ pce <dbl> 531.5, 534.2, 544.9, 544.6, 550.4, 556.8, 563.8, 567.... $ pop <int> 199808, 199920, 200056, 200208, 200361, 200536, 20070... $ psavert <dbl> 11.7, 12.2, 11.6, 12.2, 12.0, 11.6, 10.6, 10.4, 10.4,... $ uempmed <dbl> 5.1, 4.5, 4.1, 4.6, 4.4, 4.4, 4.5, 4.2, 4.6, 4.8, 4.4... $ unemploy <int> 2878, 3001, 2877, 2709, 2740, 2938, 2883, 2768, 2686,... $ Class <chr> "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: Naive method Model Information: Call: naive(y = dat_ts, h = 12) Residual sd: 214.9766 Error measures: ME RMSE MAE MPE MAPE MASE Training set 13.60799 214.9766 161.1361 0.1931667 2.148826 0.1690101 ACF1 Training set 0.201558 Forecasts: Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 Jan 2014 10376 10100.496 10651.50 9954.654 10797.35 Feb 2014 10376 9986.379 10765.62 9780.126 10971.87 Mar 2014 10376 9898.814 10853.19 9646.206 11105.79 Apr 2014 10376 9824.993 10927.01 9533.307 11218.69 May 2014 10376 9759.955 10992.04 9433.841 11318.16 Jun 2014 10376 9701.157 11050.84 9343.916 11408.08 Jul 2014 10376 9647.086 11104.91 9261.222 11490.78 Aug 2014 10376 9596.758 11155.24 9184.252 11567.75 Sep 2014 10376 9549.489 11202.51 9111.961 11640.04 Oct 2014 10376 9504.781 11247.22 9043.585 11708.41 Nov 2014 10376 9462.258 11289.74 8978.552 11773.45 Dec 2014 10376 9421.627 11330.37 8916.413 11835.59

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 Training set 13.63606 214.7896 160.8971 0.1946177 2.146727 0.1687594 ACF1 Training set 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 Training set -1.775468 199.9836 152.9176 0.04058417 2.056684 0.16039 ACF1 Training set -0.001011626 Forecasts: Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 Jan 2014 10194.656 9937.433 10451.879 9801.267 10588.04 Feb 2014 9973.102 9590.814 10355.390 9388.443 10557.76 Mar 2014 9751.549 9239.463 10263.634 8968.381 10534.72 Apr 2014 9529.995 8881.048 10178.943 8537.515 10522.48 May 2014 9308.442 8514.997 10101.887 8094.972 10521.91 Jun 2014 9086.888 8141.265 10032.512 7640.682 10533.10 Jul 2014 8865.335 7759.990 9970.680 7174.856 10555.81 Aug 2014 8643.782 7371.378 9916.185 6697.808 10589.76 Sep 2014 8422.228 6975.652 9868.804 6209.881 10634.58 Oct 2014 8200.675 6573.036 9828.313 5711.416 10689.93 Nov 2014 7979.121 6163.744 9794.498 5202.741 10755.50 Dec 2014 7757.568 5747.979 9767.157 4684.167 10830.97