将 ARIMA 模型拟合到多个时间序列并将预测存储到矩阵中答案

【问题标题】：Fitting ARIMA model to multiple time series and storing forecast into a matrix将 ARIMA 模型拟合到多个时间序列并将预测存储到矩阵中
【发布时间】：2016-10-22 18:18:49
【问题描述】：

由于它很大，我不能在这里dput 它。但是假设realmatrix 是一个具有非平凡值的“mts”

realmatrix <- matrix(NA, ncol = 100, nrow = 138)

实际上它存储了 100 个时间序列，长度（行）= 138（从 2005 年 1 月到 2016 年 6 月）。

我想通过以下循环将 Arima 预测（提前 12 个月：即从 2016 年 7 月到 2017 年 6 月）存储在另一个矩阵 farimamatrix（应该有 12 行和 100 列）中：

farimamatrix <- matrix(NA, nrow = 12, ncol = 100)

m <- k <- list()

for (i in 1:100) {
  try(m[[i]] <- Arima(realmatrix[,i], order = c(0,1,0), seasonal = c(1,0,1)))
  k[[i]] <- forecast.Arima(m[[i]], h=12)
  farimamatrix[,i] <- fitted(k[[i]])
  }

但我收到以下消息：

farimatrix[, i]

矩阵上的下标数不正确

怎么了？提前致谢。

已编辑 (24/10)：在哲元的回答下更新/更正，之前的问题消失了

原始数据：

tsdata <- 
structure(c(28220L, 27699L, 28445L, 29207L, 28482L, 28326L, 28322L, 
28611L, 29187L, 29145L, 29288L, 29352L, 28881L, 29383L, 29898L, 
29888L, 28925L, 29069L, 29114L, 29886L, 29917L, 30144L, 30531L, 
30494L, 30700L, 30325L, 31313L, 32031L, 31383L, 30767L, 30500L, 
31181L, 31736L, 32136L, 32654L, 32305L, 31856L, 31731L, 32119L, 
31953L, 32300L, 31743L, 32150L, 33014L, 32964L, 33674L, 33410L, 
31559L, 30667L, 30495L, 31978L, 32043L, 30945L, 30715L, 31325L, 
32262L, 32717L, 33420L, 33617L, 34123L, 33362L, 33731L, 35118L, 
35027L, 34298L, 34171L, 33851L, 34715L, 35184L, 35190L, 35079L, 
35958L, 35875L, 35446L, 36352L, 36050L, 35567L, 35161L, 35419L, 
36337L, 36967L, 36745L, 36370L, 36744L, 36303L, 36899L, 38621L, 
37994L, 36809L, 36527L, 35916L, 37178L, 37661L, 37794L, 38642L, 
37763L, 38367L, 38006L, 38442L, 38654L, 38345L, 37628L, 37698L, 
38613L, 38525L, 39389L, 39920L, 39556L, 40280L, 41653L, 40269L, 
39592L, 39100L, 37726L, 37867L, 38551L, 38895L, 40100L, 40950L, 
39838L, 40643L, 40611L, 39611L, 39445L, 38059L, 37131L, 36697L, 
37746L, 37733L, 39188L, 39127L, 38554L, 38219L, 38497L, 39165L, 
40077L, 38370L, 37174L), .Dim = c(138L, 1L), .Dimnames = list(
    NULL, "Data"), .Tsp = c(2005, 2016.41666666667, 12), class = "ts")

代码

library("forecast")

z <- stl(tsdata[, "Data"], s.window="periodic")

t <- z$time.series[,"trend"]
s <- z$time.series[,"seasonal"]
e <- z$time.series[,"remainder"]

# error matrix
ematrix <- matrix(rnorm(138 * 100, sd = 100), nrow = 138)

# generating a ts class error matrix
ematrixts <- ts(ematrix, start=c(2005,1), freq=12)

# combining the trend + season + error matrix into a real matrix
realmatrix <- t + s + ematrixts

# creating a (forecast) arima matrix
farimamatrix <- matrix(NA, ncol = 100, nrow = 12)

m <- k <- vector("list", length = 100)

for (i in 1:100) {
  try(m[[i]] <- Arima(realmatrix[,i], order = c(0,1,0), seasonal = c(1,0,1)))
  print(i)
  k[[i]] <- forecast.Arima(m[[i]], h = 12)
  farimamatrix[,i] <- k[[i]]$mean
  }

# ts.plot(farimamatrix[,1:100],col = c(rep("gray",100),rep("red",1)))

循环似乎可以工作，但由于Arima 失败，在几次迭代后就崩溃了：

统计错误::arima(x = x, order = order,seasonal =season, include.mean = include.mean, : "来自 CSS 的非平稳季节性 AR 部分

【问题讨论】：

(1) 您能否粘贴来自dput(head(realmatrix, 20)) 的输出（或给我们一个假矩阵来处理仍会产生错误消息），以便我们重现您的问题？见How to make a great R reproducible example。 (2) fitted() 会给你一个长度为nrow(realmatrix) 的向量，但farimamatrix 有12 行。 farimamatrix 的期望内容是什么？是预测吗？
@WeihuangWong 确定。 (1) 我已将 dput 包含在已编辑的帖子中。 (2) 是的，farimatrix 的目的是存储每个拟合的 arima 的预测（12 行/月）。
代码不完整，没有库语句。
realmatrix 似乎是一个向量，基于dput——对吗？这让我感到困惑，因为脚本中的对象名称和 realmatrix[,i] 都向我表明它应该是一个矩阵。
确实是矩阵。它存储 100 个时间序列，长度（行）= 138（从 2005 年 1 月到 2016 年 6 月）。我想将 arima 预测（提前 12 个月：即从 2016 年 7 月到 2017 年 6 月）存储在另一个矩阵中（自然应该有 12 行和 100 列）

标签： r list matrix time-series

【解决方案1】：

是的，之前的问题已经解决了，现在你有一个新问题，关于Arima 的失败。严格来说，你应该就此提出一个新问题。不过我还是在这里回答吧。

错误信息非常具有说明性。当您拟合模型ARIMA(0,1,0)(1,0,1) 时，有时季节性部分是非平稳的，因此需要进一步的季节性差异。

通过查看ts.plot(realmatrix)，我看到realmatrix 的所有100 列都非常相似。因此，我将取出第一列进行一些分析。

x <- realmatrix[,1]

显然非季节性差异是必须的，但我们是否也需要季节性差异？与 ACF 核对一下

acf(diff(x))

我们实际上发现了季节性模式的有力证据。所以是的，需要季节性差异。

现在让我们检查两个差分后的 ACF：

acf(diff(diff(x, lag = 12)))  ## first do seasonal diff, then non-seasonal diff

季节之间似乎出现负峰值，表明存在季节性 MA 过程。所以ARIMA(0,1,0)(0,1,1)[12] 是个不错的选择。

fit <- arima(x, order = c(0,1,0), seasonal = c(0,1,1))

检查残差：

acf(fit$residuals)

我实际上会对这个结果感到非常高兴，因为根本没有滞后 1 甚至滞后 2 自相关，而且也没有季节性自相关。您实际上可以尝试进一步添加季节性和/或非季节性 AR(1)，但不会有任何改进。所以这是我们的最终模型。

所以使用下面的循环：

farimamatrix <- matrix(NA, ncol = 100, nrow = 12)

m <- k <- vector("list", length = 100)

for (i in 1:100) {
  m[[i]] <- Arima(realmatrix[,i], order = c(0,1,0), seasonal = c(0,1,1))
  print(i)
  k[[i]] <- forecast.Arima(m[[i]], h = 12)
  farimamatrix[,i] <- k[[i]]$mean
  }

现在100个模型全部拟合成功。

---------

回顾性反思

也许我应该在最初的答案中解释为什么ARIMA(0,1,0)(1,0,1)[12] 模型适用于我的模拟数据。因为请注意我如何模拟我的数据：

seasonal <- rep_len(sin((1:12) * pi / 6), 138)

是的，潜在的季节性模式是真实的复制，当然是固定的。

【讨论】：