CODEWITHSUNDEEP
rtime-serieshierarchyforecast
user9898927
user9898927

Reputation: 83

How to use weighted last squares while reconciling custom forecasts using hts::combinef()?

I am trying to reconcile custom forecasts using the combinef function from the hts package. I want to compare recombined automated approaches and recombined custom forecasts (https://otexts.com/fpp2/reconciliation.html)

The automated coutertpart for combinef (forecast.gts) is very user friendly. If a bottom up forecast is required then you can set method = "comb". If reconciliation is required then you can for example choose between weighted last squares, ordinary last squares and structural scaling by setting the parameter weights to weights = c("wls", "ols", "nseries").

For combinef() the default of the weights parameter are the ordinary last squares, therefore this approach can be easily implemented. In another thread it was already explained how to apply the bottom up approach (How to get top down forecasts using `hts::combinef()`?).

I am now interested in also using "wls" as well als "nseries". Is there an easy way to implement this?

Edit I had a deeper look into the raw code of the forecast.gts() and the combinef() function. I ended up with this custom code:

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
)

Will this do the trick?

Upvotes: 1

Views: 185

Answers (1)

Rob Hyndman
Rob Hyndman

Reputation: 31800

WLS forecast reconciliation uses the one-step forecast variances which are equal to the residual variances. Here is some code to do it in a small example.

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")

Created on 2020-07-03 by the reprex package (v0.3.0)

Upvotes: 1

Related Questions