使用 R 的时间序列分析程序和方法

机器算法验证 r 时间序列 预测 神经网络 有马
2022-02-13 16:27:13

我正在做一个小项目,我们试图预测未来 6 个月的商品价格(石油、铝、锡等)。我有 12 个这样的变量要预测,并且我有 2008 年 4 月至 2013 年 5 月的数据。

我应该如何进行预测?我做了以下事情:

  • 导入的数据作为时间序列数据集
  • 所有变量的季节性往往随趋势变化,所以我将使用乘法模型。
  • 我记录了变量以转换为加法模型
  • 对于每个变量,使用 STL 分解数据

我打算使用 Holt Winters 指数平滑、ARIMA 和神经网络进行预测。我将数据拆分为训练和测试 (80, 20)。计划选择MAE、MPE、MAPE、MASE少的机型。

我做对了吗?

我还有一个问题是,在传递给 ARIMA 或神经网络之前,我应该对数据进行平滑处理吗?如果是,使用什么?数据显示季节性和趋势。

编辑:

附加时间序列图和数据 在此处输入图像描述

Year  <- c(2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2009, 2009, 
           2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2010, 
           2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 
           2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 
           2011, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 
           2012, 2012, 2013, 2013)
Month <- c(4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
           12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 
           8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2) 
Coil  <- c(44000, 44500, 42000, 45000, 42500, 41000, 39000, 35000, 34000, 
           29700, 29700, 29000, 30000, 30000, 31000, 31000, 33500, 33500, 
           33000, 31500, 34000, 35000, 35000, 36000, 38500, 38500, 35500, 
           33500, 34500, 36000, 35500, 34500, 35500, 38500, 44500, 40700, 
           40500, 39100, 39100, 39100, 38600, 39500, 39500, 38500, 39500, 
           40000, 40000, 40500, 41000, 41000, 41000, 40500, 40000, 39300, 
           39300, 39300, 39300, 39300, 39800)
coil <- data.frame(Year = Year, Month = Month, Coil = Coil)

编辑2: 一个问题,你能告诉我我的数据是否有季节性或趋势吗?也请给我一些关于如何识别它们的提示。 在此处输入图像描述 在此处输入图像描述

2个回答

您应该使用forecast 包,它支持所有这些模型(以及更多模型),并且可以快速拟合它们:

library(forecast)
x <- AirPassengers
mod_arima <- auto.arima(x, ic='aicc', stepwise=FALSE)
mod_exponential <- ets(x, ic='aicc', restrict=FALSE)
mod_neural <- nnetar(x, p=12, size=25)
mod_tbats <- tbats(x, ic='aicc', seasonal.periods=12)
par(mfrow=c(4, 1))
plot(forecast(mod_arima, 12), include=36)
plot(forecast(mod_exponential, 12), include=36)
plot(forecast(mod_neural, 12), include=36)
plot(forecast(mod_tbats, 12), include=36)

我建议不要在拟合模型之前对数据进行平滑处理。您的模型本质上会尝试平滑数据,因此预平滑只会使事情复杂化。

在此处输入图像描述

根据新数据编辑:

实际上,看起来 arima 是您可以为此训练和测试集选择的最差模型之一。

我将您的数据保存到文件调用coil.csv中,将其加载到 R 中,并将其拆分为训练和测试集:

library(forecast)
dat <- read.csv('~/coil.csv')
x <- ts(dat$Coil, start=c(dat$Year[1], dat$Month[1]), frequency=12)
test_x <- window(x, start=c(2012, 3))
x <- window(x, end=c(2012, 2))

接下来我拟合了一堆时间序列模型:arima、指数平滑、神经网络、tbats、bats、季节性分解和结构时间序列:

models <- list(
  mod_arima = auto.arima(x, ic='aicc', stepwise=FALSE),
  mod_exp = ets(x, ic='aicc', restrict=FALSE),
  mod_neural = nnetar(x, p=12, size=25),
  mod_tbats = tbats(x, ic='aicc', seasonal.periods=12),
  mod_bats = bats(x, ic='aicc', seasonal.periods=12),
  mod_stl = stlm(x, s.window=12, ic='aicc', robust=TRUE, method='ets'),
  mod_sts = StructTS(x)
  )

然后我做了一些预测并与测试集进行比较。我包括了一个天真的预测,它总是预测一条平坦的水平线:

forecasts <- lapply(models, forecast, 12)
forecasts$naive <- naive(x, 12)
par(mfrow=c(4, 2))
for(f in forecasts){
  plot(f)
  lines(test_x, col='red')
}

在此处输入图像描述

如您所见,华宇模型弄错了趋势,但我有点喜欢“基本结构模型”的外观

最后,我在测试集上测量了每个模型的准确度:

acc <- lapply(forecasts, function(f){
  accuracy(f, test_x)[2,,drop=FALSE]
})
acc <- Reduce(rbind, acc)
row.names(acc) <- names(forecasts)
acc <- acc[order(acc[,'MASE']),]
round(acc, 2)
                ME    RMSE     MAE   MPE MAPE MASE ACF1 Theil's U
mod_sts     283.15  609.04  514.46  0.69 1.27 0.10 0.77      1.65
mod_bats     65.36  706.93  638.31  0.13 1.59 0.12 0.85      1.96
mod_tbats    65.22  706.92  638.32  0.13 1.59 0.12 0.85      1.96
mod_exp      25.00  706.52  641.67  0.03 1.60 0.12 0.85      1.96
naive        25.00  706.52  641.67  0.03 1.60 0.12 0.85      1.96
mod_neural   81.14  853.86  754.61  0.18 1.89 0.14 0.14      2.39
mod_arima   766.51  904.06  766.51  1.90 1.90 0.14 0.73      2.48
mod_stl    -208.74 1166.84 1005.81 -0.52 2.50 0.19 0.32      3.02

Hyndman, RJ 和 Athanasopoulos, G. (2014) “预测:原则和实践”中描述了使用的指标,他们也恰好是预测包的作者。我强烈建议您阅读他们的文字:它可以在线免费获得。结构化时间序列是几个指标的最佳模型,包括 MASE,这是我倾向于选择模型的指标。

最后一个问题是:结构模型在这个测试集上是否幸运?评估这一点的一种方法是查看训练集错误。训练集错误不如测试集错误可靠(因为它们可能过拟合),但在这种情况下,结构模型仍然是最重要的:

acc <- lapply(forecasts, function(f){
  accuracy(f, test_x)[1,,drop=FALSE]
})
acc <- Reduce(rbind, acc)
row.names(acc) <- names(forecasts)
acc <- acc[order(acc[,'MASE']),]
round(acc, 2)
                ME    RMSE     MAE   MPE MAPE MASE  ACF1 Theil's U
mod_sts      -0.03    0.99    0.71  0.00 0.00 0.00  0.08        NA
mod_neural    3.00 1145.91  839.15 -0.09 2.25 0.16  0.00        NA
mod_exp     -82.74 1915.75 1359.87 -0.33 3.68 0.25  0.06        NA
naive       -86.96 1936.38 1386.96 -0.34 3.75 0.26  0.06        NA
mod_arima  -180.32 1889.56 1393.94 -0.74 3.79 0.26  0.09        NA
mod_stl     -38.12 2158.25 1471.63 -0.22 4.00 0.28 -0.09        NA
mod_bats     57.07 2184.16 1525.28  0.00 4.07 0.29 -0.03        NA
mod_tbats    62.30 2203.54 1531.48  0.01 4.08 0.29 -0.03        NA

(注意神经网络过拟合,在训练集上表现出色,在测试集上表现不佳)

最后,交叉验证所有这些模型是一个好主意,也许通过 2008-2009 年培训/2010 年测试、2008-2010 年培训/2011 年测试、2008-2011 年培训/2012 年测试、培训2008-2012 年/2013 年测试,并在所有这些时间段内平均误差。如果你想走这条路,我有一个部分完整的包,用于在 github 上交叉验证时间序列模型,我希望你尝试一下并给我反馈/请求:

devtools::install_github('zachmayer/cv.ts')
library(cv.ts)

编辑2:让我们看看我是否记得如何使用我自己的包!

首先,从 github 安装并加载包(见上文)。然后交叉验证一些模型(使用完整的数据集):

library(cv.ts)
x <- ts(dat$Coil, start=c(dat$Year[1], dat$Month[1]), frequency=12)
ctrl <- tseriesControl(stepSize=1, maxHorizon=12, minObs=36, fixedWindow=TRUE)
models <- list()

models$arima = cv.ts(
  x, auto.arimaForecast, tsControl=ctrl,
  ic='aicc', stepwise=FALSE)

models$exp = cv.ts(
  x, etsForecast, tsControl=ctrl,
  ic='aicc', restrict=FALSE)

models$neural = cv.ts(
  x, nnetarForecast, tsControl=ctrl,
  nn_p=6, size=5)

models$tbats = cv.ts(
  x, tbatsForecast, tsControl=ctrl,
  seasonal.periods=12)

models$bats = cv.ts(
  x, batsForecast, tsControl=ctrl,
  seasonal.periods=12)

models$stl = cv.ts(
  x, stl.Forecast, tsControl=ctrl,
  s.window=12, ic='aicc', robust=TRUE, method='ets')

models$sts = cv.ts(x, stsForecast, tsControl=ctrl)

models$naive = cv.ts(x, naiveForecast, tsControl=ctrl)

models$theta = cv.ts(x, thetaForecast, tsControl=ctrl)

(请注意,我降低了神经网络模型的灵活性,以防止它过度拟合)

一旦我们拟合了模型,我们就可以通过 MAPE 来比较它们(cv.ts 还不支持 MASE):

res_overall <- lapply(models, function(x) x$results[13,-1])
res_overall <- Reduce(rbind, res_overall)
row.names(res_overall) <- names(models)
res_overall <- res_overall[order(res_overall[,'MAPE']),]
round(res_overall, 2)
                 ME    RMSE     MAE   MPE MAPE
naive     91.40 1126.83  961.18  0.19 2.40
ets       91.56 1127.09  961.35  0.19 2.40
stl     -114.59 1661.73 1332.73 -0.29 3.36
neural     5.26 1979.83 1521.83  0.00 3.83
bats     294.01 2087.99 1725.14  0.70 4.32
sts     -698.90 3680.71 1901.78 -1.81 4.77
arima  -1687.27 2750.49 2199.53 -4.23 5.53
tbats   -476.67 2761.44 2428.34 -1.23 6.10

哎哟。看来我们的结构性预测很幸运。从长远来看,幼稚的预测是最好的预测,在 12 个月的范围内进行平均(arima 模型仍然是最差的模型之一)。让我们比较 12 个预测范围中的每个模型,看看它们中是否有任何一个击败了幼稚模型:

library(reshape2)
library(ggplot2)
res <- lapply(models, function(x) x$results$MAPE[1:12])
res <- data.frame(do.call(cbind, res))
res$horizon <- 1:nrow(res)
res <- melt(res, id.var='horizon', variable.name='model', value.name='MAPE')
res$model <- factor(res$model, levels=row.names(res_overall))
ggplot(res, aes(x=horizon, y=MAPE, col=model)) +
  geom_line(size=2) + theme_bw() +
  theme(legend.position="top") +
  scale_color_manual(values=c(
    "#1f78b4", "#ff7f00", "#33a02c", "#6a3d9a",
    "#e31a1c", "#b15928", "#a6cee3", "#fdbf6f",
    "#b2df8a")
    )

模型比较

很明显,指数平滑模型总是选择朴素模型(橙色线和蓝色线重叠 100%)。换句话说,“下个月的卷材价格将与本月的卷材价格相同”的幼稚预测(几乎在每个预测范围内)比 7 个极其复杂的时间序列模型更准确。除非你有一些卷材市场还不知道的秘密信息,否则要打破幼稚的卷材价格预测将非常困难

这绝不是任何人想听到的答案,但如果预测准确性是您的目标,您应该使用最准确的模型。使用朴素模型。

你采取的方法是合理的。如果您不熟悉预测,那么我会推荐以下书籍:

  1. Makridakis、Wheelright 和 Hyndman 的预测方法和应用
  2. 预测: Hyndman 和 Athanasopoulos 的原则和实践。

第一本书是我强烈推荐的经典。第二本书是一本开源书籍,您可以参考它来了解预测方法以及如何使用 R开源软件包预测来应用它。这两本书都为我使用的方法提供了很好的背景。如果您对预测很认真,那么我会推荐Armstrong的《预测原理》,它收集了大量的预测研究,从业者可能会发现它非常有帮助。

谈到你关于线圈的具体例子,它让我想起了大多数教科书经常忽略的可预测性概念。某些系列(例如您的系列)根本无法预测,因为它的模式较少,因为它没有表现出趋势或季节性模式或任何系统变化。在这种情况下,我会将一个系列归类为不可预测的系列。在尝试外推方法之前,我会查看数据并提出以下问题:我的序列是否可预测?在这个具体示例中,使用预测的最后一个值的随机游走预测等简单外推被发现是最准确的.

还有关于神经网络的另一条评论:众所周知,神经网络在经验竞赛中会失败。在尝试使用神经网络进行时间序列预测任务之前,我会尝试传统的时间序列统计方法。

我试图在 中对您的数据进行建模R's forecast package,希望这些评论是不言自明的。

coil <- c(44000, 44500, 42000, 45000, 42500, 41000, 39000, 35000, 34000, 
          29700, 29700, 29000, 30000, 30000, 31000, 31000, 33500, 33500, 
          33000, 31500, 34000, 35000, 35000, 36000, 38500, 38500, 35500, 
          33500, 34500, 36000, 35500, 34500, 35500, 38500, 44500, 40700, 
          40500, 39100, 39100, 39100, 38600, 39500, 39500, 38500, 39500, 
          40000, 40000, 40500, 41000, 41000, 41000, 40500, 40000, 39300, 
          39300, 39300, 39300, 39300, 39800)


coilts <- ts(coil,start=c(2008,4),frequency=12)

library("forecast")

# Data for modeling
coilts.mod <- window(coilts,end = c(2012,3))

#Data for testing
coil.test <- window(coilts,start=c(2012,4))

# Model using multiple methods - arima, expo smooth, theta, random walk, structural time series

#arima
coil.arima <- forecast(auto.arima(coilts.mod),h=11)

#exponential smoothing
coil.ets <- forecast(ets(coilts.mod),h=11)

#theta
coil.tht <- thetaf(coilts.mod, h=11)

#random walk
coil.rwf <- rwf(coilts.mod, h=11)

#structts
coil.struc <- forecast(StructTS(coilts.mod),h=11)


##accuracy 

arm.acc <- accuracy(coil.arima,coil.test)
ets.acc <- accuracy(coil.ets,coil.test)
tht.acc <- accuracy(coil.tht,coil.test)
rwf.acc <- accuracy(coil.rwf,coil.test)
str.acc <- accuracy(coil.struc,coil.test)

在保留数据上使用 MAE,我会选择 ARIMA 进行短期预测(1 - 12 个月)。从长远来看,我会依赖随机游走预测。请注意,ARIMA 选择了具有漂移(0,1,0)+漂移的随机游走模型,在这类问题中,特别是短期内,它往往比纯随机游走模型准确得多。见下图。这是基于上面代码中所示的精度函数。

在此处输入图像描述

您的具体问题的具体答案:我还有一个问题是,在传递给 ARIMA 或神经网络之前,我应该对数据进行平滑处理吗?如果是,使用什么?

  • 不,预测方法自然会平滑您的数据以适合模型。

数据显示季节性和趋势。

  • 以上数据不显示趋势或季节性。如果您确定数据表现出季节性和趋势,则选择适当的方法。

提高准确性的实用技巧:

结合多种预测方法: - 您可以尝试使用非外推方法,例如类比预测、判断预测或其他技术,并将它们与您的统计方法相结合,以提供准确的预测。有关组合的好处,请参阅本文我尝试将上述 5 种方法结合起来,但作为个体方法,预测并不准确,一个可能的原因是个体预测相似。当您结合统计和判断预测等多种方法时,您将获得结合预测的好处。

检测和理解异常值: - 现实世界的数据充满了异常值。识别并适当处理时间序列中的异常值。推荐阅读这篇文章在查看您的线圈数据时,2009 年之前的下降是异常值吗?

编辑

这些数据似乎遵循某种类型的宏观经济趋势。我的猜测是,在 2009 年之前出现的下降趋势是在 2008 年至 2009 年之间出现的经济衰退之后,并在 2009 年之后开始回升。如果是这种情况,那么我将一起避免使用任何外推方法,而是依靠合理的理论来了解如何这些经济趋势的表现与@GraemeWalsh所引用的一样。

希望这可以帮助