查看末尾的完整功能示例，了解如何使用 XREG 和 处理带有 DAILY DATA 的 AUTO.ARIMA MODEL具有滚动开始时间的傅里叶级数以及经过交叉验证的训练和测试。
- 如果没有可重现的示例，没有人可以帮助您，因为他们无法运行您的代码。您需要提供数据。 :-(
- 即使讨论统计问题不是 StackOverflow 的一部分，为什么不使用
xreg 而不是 lm + auto.arima 对残差进行 auto.arima 处理？尤其是考虑到你最后的预测方式，这种训练方法看起来确实是错误的。考虑使用：
fit <- auto.arima(y, xreg = lagmatrix(f, -1))
h <- forecast(fit, xreg = lagmatrix(e, -1))
auto.arima 会自动按照最大似然计算最佳参数。
- 关于您的编码问题..
forc <- c() 应该在for 循环之外，否则在每次运行时都会删除之前的结果。
j <- 0 也一样：每次运行时都将其设置回 0。如果您需要在每次运行时更改其值，请将其放在外面。
forecast 的输出是forecast 类的对象，实际上是list 的一个类型。因此，您不能有效地使用cbind。
我的意见是，你应该这样创建forc：forc <- list()
并以这种方式创建最终结果列表：
forc <- c(forc, list(h)) # instead of forc <- cbind(forc, h)
这将创建forecast 类的对象列表。
然后，您可以使用for 循环通过访问每个对象或使用lapply 来绘制它们。
lapply(output_of_your_function, plot)
如果没有可重现的例子，我可以做到这一点。
最终编辑
完整的功能示例
在这里，我试图从我们写的百万 cmets 中总结出一个结论。
使用您提供的数据，我构建了一个可以处理您需要的所有内容的代码。
从训练和测试到模型，直到预测并最终绘制 X 轴和您的一个 cmets 中所需的时间。
我删除了for 循环。 lapply 更适合您的情况。
如果你愿意，你可以离开傅立叶级数。 That's how Professor Hyndman 建议处理每日时间序列。
需要的函数和库：
# libraries ---------------------------
library(forecast)
library(lubridate)
# run model -------------------------------------
.daily_arima_forecast <- function(init, training, horizon, tt, ..., K = 10){
# create training and test
tt_trn <- window(tt, start = time(tt)[init] , end = time(tt)[init + training - 1])
tt_tst <- window(tt, start = time(tt)[init + training], end = time(tt)[init + training + horizon - 1])
# add fourier series [if you want to. Otherwise, cancel this part]
fr <- fourier(tt_trn[,1], K = K)
frf <- fourier(tt_trn[,1], K = K, h = horizon)
tsp(fr) <- tsp(tt_trn)
tsp(frf) <- tsp(tt_tst)
tt_trn <- ts.intersect(tt_trn, fr)
tt_tst <- ts.intersect(tt_tst, frf)
colnames(tt_tst) <- colnames(tt_trn) <- c("y", "s", paste0("k", seq_len(ncol(fr))))
# run model and forecast
aa <- auto.arima(tt_trn[,1], xreg = tt_trn[,-1])
fcst <- forecast(aa, xreg = tt_tst[,-1])
# add actual values to plot them later!
fcst$test.values <- tt_tst[,1]
# NOTE: since I modified the structure of the class forecast I should create a new class,
# but I didnt want to complicate your code
fcst
}
daily_arima_forecast <- function(y, x, training, horizon, ...){
# set up x and y together
tt <- ts.intersect(y, x)
# set up all starting point of the training set [give it a name to recognize them later]
inits <- setNames(nm = seq(1, length(y) - training, by = horizon))
# remove last one because you wouldnt have enough data in front of it
inits <- inits[-length(inits)]
# run model and return a list of all your models
lapply(inits, .daily_arima_forecast, training = training, horizon = horizon, tt = tt, ...)
}
# plot ------------------------------------------
plot_daily_forecast <- function(x){
autoplot(x) + autolayer(x$test.values)
}
关于如何使用先前功能的可重现示例
# create a sample data
tsp(EuStockMarkets) <- c(1991, 1991 + (1860-1)/365.25, 365.25)
# model
models <- daily_arima_forecast(y = EuStockMarkets[,1],
x = EuStockMarkets[,2],
training = 600,
horizon = 25,
K = 5)
# plot
plots <- lapply(models, plot_daily_forecast)
plots[[1]]
文章作者示例
# your data
load("BVIS0157_Forward.rda")
load("BVIS0157_NS.spread.rda")
spread.NS.JPM <- spread.NS.JPM / 100
# pre-work [out of function!!!]
set_up_ts <- function(m){
start <- min(row.names(m))
end <- max(row.names(m))
# daily sequence
inds <- seq(as.Date(start), as.Date(end), by = "day")
ts(m, start = c(year(start), as.numeric(format(inds[1], "%j"))), frequency = 365.25)
}
mts_spread.NS.JPM <- set_up_ts(spread.NS.JPM)
mts_Forward.rate.JPM <- set_up_ts(Forward.rate.JPM)
# model
col <- 10
models <- daily_arima_forecast(y = mts_spread.NS.JPM[, col],
x = stats::lag(mts_Forward.rate.JPM[, col], -1),
training = 600,
horizon = 25,
K = 5) # notice that K falls between ... that goes directly to the inner function
# plot
plots <- lapply(models, plot_daily_forecast)
plots[[5]]