带非负约束的MinT层级预测调和：fable与hts工具方案问询

最优方案思路 #

直接在MinT调和的优化过程中加入非负约束（通过二次规划实现），既保留fable框架的自定义模型灵活性，又保证调和结果非负，同时支持分布型预测输出。

这种方案比你提出的事后截断修正更严谨——事后截断会破坏MinT的最优迹最小化特性，而约束优化是从根源上满足非负要求，同时严格保留层级一致性。

具体实现步骤 #

1. 用fable生成自定义基准预测（带分布） #

通过fable构建任意你需要的基准模型（比如指定ARIMA(d=1)、ETS、甚至机器学习模型），确保基准预测的点估计非负（可通过对数变换、Box-Cox变换，或在模型中直接约束输出），同时保留预测分布信息（残差抽样、参数化分布均可）。

示例代码：

library(tidyverse)
library(fable)
library(fabletools)
library(hts)
library(quadprog)

# 构建层级时间序列示例数据
hier_data <- aus_retail %>%
  filter(Industry == "Department stores") %>%
  mutate(State = factor(State)) %>%
  as_hierarchical(State ~ .)

# 自定义基准模型：指定ARIMA(0,1,1)，通过对数变换保证基准预测非负
bench_forecasts <- hier_data %>%
  model(
    arima = ARIMA(log(Turnover) ~ pdq(0,1,1) + PDQ(0,0,0))
  ) %>%
  forecast(h = 12) %>%
  mutate(Turnover = exp(Turnover))  # 转换回原始尺度，确保基准预测非负

2. 实现带非负约束的MinT调和 #

MinT的核心是求解二次规划问题：$\hat{y}h^* = \arg\min{y_h \in \mathcal{H}} (y_h - \hat{y}_h)^T W^{-1} (y_h - \hat{y}_h)$，其中$\mathcal{H}$是层级一致约束集合。我们只需在这个问题中加入$y_h \geq 0$的约束即可。

以下是封装好的调和函数，支持非负约束和分布生成：

# 带非负约束的MinT调和函数
min_trace_nonneg <- function(forecasts, weights = "mint_shrink", bootstrap = 1000) {
  # 提取层级结构与基准预测
  gts_obj <- as.gts(forecasts)
  y_hat <- gts_obj$bts
  S <- gts_obj$nodes$S  # 层级求和矩阵
  n_bottom <- nrow(S)   # 底层序列数量
  
  # 构建MinT权重矩阵W
  W <- create_weights(gts_obj, method = weights)
  W_inv <- solve(W)
  
  # 二次规划参数：将问题转化为底层变量的优化，加入非负约束
  Dmat <- 2 * t(S) %*% W_inv %*% S
  dvec <- 2 * t(S) %*% W_inv %*% y_hat
  Amat <- cbind(t(S), diag(n_bottom))  # 层级等式约束 + 非负不等式约束
  bvec <- c(gts_obj$agg, rep(0, n_bottom))  # 顶层预测值 + 非负下限
  meq <- nrow(gts_obj$agg)  # 层级约束为等式约束
  
  # 求解调和点估计
  qp_sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = meq)
  reconciled_point <- S %*% qp_sol$solution
  
  # 生成分布预测（基于残差自举）
  if(bootstrap > 0){
    # 提取基准模型残差
    resids <- residuals(model(bench_forecasts)) %>%
      pull(.resid) %>%
      na.omit()
    
    # 自举抽样生成调和分布样本
    reconciled_dist <- replicate(bootstrap, {
      y_hat_boot <- y_hat + sample(resids, length(y_hat), replace = TRUE)
      dvec_boot <- 2 * t(S) %*% W_inv %*% y_hat_boot
      qp_sol_boot <- solve.QP(Dmat, dvec_boot, Amat, bvec, meq = meq)
      S %*% qp_sol_boot$solution
    }) %>% t()
    
    # 整理分布结果为tidy格式
    dist_tbl <- tibble(
      .rep = 1:bootstrap,
      .value = as.vector(reconciled_dist),
      .key = rep(gts_obj$labels$key, bootstrap)
    ) %>%
      pivot_wider(names_from = .key, values_from = .value)
  }
  
  # 返回结果：点估计 + 分布样本
  list(
    point = tibble(
      key = gts_obj$labels$key,
      value = as.vector(reconciled_point)
    ),
    distribution = if(bootstrap > 0) dist_tbl else NULL
  )
}

# 调用调和函数
reconciled_result <- min_trace_nonneg(bench_forecasts, weights = "mint_shrink", bootstrap = 1000)

3. 结果验证 #

• 调和后的所有层级预测值均非负
• 严格满足层级一致性（底层求和等于顶层）
• 分布样本保留了原始基准预测的不确定性，同时符合调和约束

方案优势对比 #

对比你的临时方案：

1. 最优性保留：直接在MinT的优化过程中加入非负约束，不会破坏迹最小化的最优性准则
2. 层级一致性严格满足：避免事后截断导致的层级不一致问题（比如底层设0后顶层手动调整引入的误差）
3. 分布合理性：自举过程中每一次抽样都经过约束调和，生成的分布更符合实际业务逻辑

内容的提问来源于stack exchange，提问作者lowndrul