使用滚动 Window 在一个图中绘制预测值和实际值

Question

我有一个代码，它将输入作为收益率利差（相关变量）和远期利率（独立变量）并运行 auto.arima 来获取订单。之后，我预测接下来的 25 个日期 (forc.horizon)。我的训练数据是前600（训练）。然后我将时间 window 移动到 25 个日期，这意味着使用从 26 到 625 的数据，估计 auto.arima 然后预测从 626 到 650 的数据等等。我的数据集是 2298 行（日期）和 30 列（成熟度）。

我想存储所有预测值，然后在同一图中绘制预测值和实际值。

这是我的代码，但它没有以稍后绘制的方式存储预测。

forecast.func <- function(NS.spread, ind.v, maturity, training, forc.horizon){
  
  NS.spread <- NS.spread/100
  forc <- c()
  j <- 0
  
  for(i in 1:floor((nrow(NS.spread)-training)/forc.horizon)){
    
    
    # test data
    y <- NS.spread[(1+j):(training+j) , maturity]
    f <- ind.v[(1+j):(training+j) , maturity]
        
    # auto- arima
    c <- auto.arima(y, xreg = f, test= "adf")

    
    # forecast
    e <- ind.v[(training+j+1):(training+j+forc.horizon) , maturity]
    h <- forecast(c, xreg = lagmatrix(e, -1))
    
    forc <- c(forc, list(h))
    
    j <- j + forc.horizon

  }
  
  return(forc)
}


a <- forecast.func(spread.NS.JPM, Forward.rate.JPM, 10, 600, 25)
lapply(a, plot)

这是我的两个数据集的 link： https://drive.google.com/drive/folders/1goCxllYHQo3QJ0IdidKbdmfR-DZgrezN?usp=sharing

Answer 1

查看最后的完整功能示例，了解如何使用 XREG 和 FOURIER SERIES 具有 ROLLING STARTING TIMES 和交叉验证的训练和测试。

没有可重现的示例，没有人可以帮助您，因为他们无法运行您的代码。你需要提供数据。 :-(

即使它不是 Whosebug 的一部分来讨论统计问题，你为什么不用 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。如果您需要在每个运行.[=38 更改其值，请将其放在外面=]

forecast的输出是classforecast的一个对象，实际上是list的一个类型。因此，您不能有效地使用 cbind。

我个人认为，您应该这样创建 forc：forc <- list()

并以这种方式创建最终结果列表：

forc <- c(forc, list(h)) # instead of forc <- cbind(forc, h)

这将创建 class forecast 的对象列表。

然后您可以通过访问每个对象或使用 lapply.

使用 for 循环绘制它们

lapply(output_of_your_function, plot)

这是我在没有可重现示例的情况下所能做到的。

最终编辑

完整的功能示例

在这里，我试图从我们写的一百万条评论中总结出一个结论。

根据您提供的数据，我构建了一个可以处理您需要的一切的代码。

从训练和测试到建模，再到预测，最后绘制 X 轴和您评论中要求的时间。

我删除了 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]]

post

作者的示例

# 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]]

使用滚动 Window 在一个图中绘制预测值和实际值

Plotting Forecast and Real values in one plot using a Rolling Window

r

plot

regression

forecasting

rolling-computation