如何在使用 hts::combinef() 协调自定义预测时使用加权最后一个平方？

Question

我正在尝试使用 hts 包中的 combinef 函数来协调自定义预测。 我想比较重组的自动化方法和重组的自定义预测 (https://otexts.com/fpp2/reconciliation.html)

combinef (forecast.gts) 的自动计算部分非常用户友好。 如果需要自下而上的预测，那么您可以设置 method = "comb"。 如果需要对账，则可以通过将参数 weights 设置为 weights = c("wls", "ols", "nseries") 在加权最后一个正方形、普通最后一个正方形和结构缩放之间进行选择。

对于 combinef()，默认的 weights 参数是普通的最后一个正方形，因此这种方法可以很容易地实现。 在另一个线程中已经解释了如何应用自下而上的方法 (How to get top down forecasts using `hts::combinef()`?)。

我现在对使用“wls”以及“nseries”也很感兴趣。 有没有简单的方法来实现这个？

编辑 我深入研究了 forecast.gts() 和 combinef() 函数的原始代码。我最终得到了这个自定义代码：

library(hts)
library(rlist)

#forecast grouped time series by custom function
ally_df <- aggts(htseg1) %>% as.data.frame
forecast_list <- apply(ally_df, 2, function(x){x %>% auto.arima %>% forecast(h = 12)})

ally_fitted <- lapply(forecast_list, function(x){x$fitted %>% as.data.frame}) %>% list.cbind
colnames(ally_fitted) <- colnames(ally)

ally_forecast <- lapply(forecast_list, function(x){x$mean %>% as.data.frame}) %>% list.cbind
colnames(ally_forecast) <- colnames(ally)


#create weights for reconciliation
recomb_approaches <- c("wls", "ols", "nseries", "bu")
recomb_approach <- recomb_approaches[1]

if(recomb_approach == "bu"){
  
  weights <- c(rep(0, ncol(ally_df)-ncol(htseg1$bts)), rep(1, ncol(htseg1$bts)))
  
}else if(recomb_approach == "ols"){
  
  weights <- NULL
  
}else if(recomb_approach == "wls"){
  
  tmp.resid <- ally_df - ally_fitted
  weights <- 1/colMeans(tmp.resid^2, na.rm = TRUE)
  
}else if(recomb_approach == "nseries"){
  
  # A function to calculate No. of groups at each level
  Mlevel <- function(xgroup) {
    m <- apply(xgroup, 1, function(x) length(unique(x)))
    return(m)
  }
  
  # A function to get the inverse of row sums of S matrix
  InvS4g <- function(xgroup) {
    mlevel <- Mlevel(xgroup)
    len <- length(mlevel)
    repcount <- mlevel[len]/mlevel
    inv.s <- 1/unlist(mapply(rep, repcount, mlevel, SIMPLIFY = FALSE))
    return(inv.s)
  }
  weights <- InvS4g(htseg1$groups)
}

ally_forecast_recombined_df <- combinef(ally_forecast
                                      , nodes = get_nodes(htseg1)
                                      , weights = weights
                                      , algorithms = "lu"
                                      , keep = "bottom"
                                      , parallel = TRUE
                                      , num.cores = cores
) 

这能解决问题吗？

Answer 1

WLS 预测调节使用等于剩余方差的一步预测方差。这里有一些代码可以在一个小例子中做到这一点。

library(hts)

h <- 12
ally <- aggts(htseg1)
nseries <- NCOL(ally)
allf <- ts(matrix(NA, nrow = h, ncol = nseries))
resid_var <- numeric(nseries)
for(i in seq(nseries)) {
  fc  <- forecast(auto.arima(ally[,i]), h = h, level=95)
  resid_var[i] <- fc$model$sigma2
  allf[,i] <- fc$mean
}
tsp(allf) <- tsp(fc$mean)
wls <- combinef(allf, get_nodes(htseg1), weights = resid_var, 
                keep = "gts", algorithms = "lu")

^{由reprex package (v0.3.0) 于 2020 年 7 月 3 日创建}