0) 安裝與載入套件
1) 參數設定
2) 讀檔與時間欄清洗（支援 2000-M1 這種格式）
3) 建立落後特徵（只用 lag，避免洩漏）
4) 指標函數與儲存容器
5) 建立 rolling 切點並執行滾動預測
6) 彙整、評分表、輸出 CSV
7) 三張樣本外走勢圖
8) Nowcast「最後一期再多推 h 期」
9) 重要變數

0) 安裝與載入套件

說明：一次安裝並載入後面會用到的套件（含 ranger 做 Random Forest）。

# ---- 安裝/載入套件 ----
pkgs <- c("data.table","readxl","lubridate","forecast","glmnet","ipred",
          "ggplot2","scales","ranger")
ins <- pkgs[!pkgs %in% installed.packages()[,"Package"]]
if(length(ins)) install.packages(ins, repos="https://cloud.r-project.org")
invisible(lapply(pkgs, library, character.only=TRUE))

1) 參數設定

說明：train_len 決定 rolling 視窗長度；h 是步數（預設 1）；target_col 換成要預測的欄名。

# ---- 參數 ----
excel_path    <- "~/Documents/文件 - Travis' MBA/OneDrive/Documents/1. R/Research/Oct_2025/scratch_3/existing_excel_files/wide_tab.xlsx"   # Excel 路徑
sheet_name    <- 1                  # 工作表
time_col      <- "ym"               # 時間欄（像 2000-M1）
target_col    <- "A100101010_1"     # <--- 請改成你的目標欄名
h             <- 1                  # 預測步數（1..12）
train_len     <- 60                 # <--- rolling window 長度（例：60 個月）
p_max         <- 6                  # 每個變數要做幾階 lag
elastic_alpha <- 0.5                # Elastic Net 的 alpha
adapt_gamma   <- 1.0                # Adaptive 權重 1/|beta|^gamma
var_exp_cut   <- 0.95               # PCR 主成分保留的累積解釋比例

2) 讀檔與時間欄清洗（支援 2000-M1 這種格式）

說明：把 "YYYY-Mm" 轉成該月第一天的 Date，再依時間排序。

# ---- 讀檔 ----
raw <- as.data.table(readxl::read_excel(excel_path, sheet = sheet_name))
setnames(raw, old=names(raw), new=make.names(names(raw)))
time_col   <- make.names(time_col)
target_col <- make.names(target_col)

stopifnot(time_col %in% names(raw), target_col %in% names(raw))

# ---- 專門解析 "YYYY-Mm" 格式成 Date ----
parse_ym_M <- function(v){
  v  <- as.character(v)
  yr <- as.integer(sub("^(\\d{4})-M.*$", "\\1", v))              # 抓年
  mm <- as.integer(sub("^\\d{4}-M(\\d{1,2})$", "\\1", v))        # 抓月
  if (any(is.na(yr)) || any(is.na(mm))) {
    stop("時間欄位必須長得像 '2025-M1' 或 '2025-M11'")
  }
  as.Date(sprintf("%04d-%02d-01", yr, mm))   # 轉換成該月第一天
}


# 先把 ym 轉成日期欄
raw[, (time_col) := parse_ym_M(get(time_col))]

# 正確的排序寫法（升冪）
data.table::setorderv(raw, cols = time_col, order = 1L, na.last = TRUE)

3) 建立落後特徵（只用 lag，避免洩漏）

說明：對所有數值欄（包含目標 y）都做 L1..Lp 的 lag。後續 ML 只會用 lag 當特徵。

# ---- 產生 lag 特徵 ----
num_vars <- setdiff(names(raw), time_col)
make_lags <- function(dt, vars, pmax){
  dt <- copy(dt)
  for (v in vars){
    for (L in 1:pmax){
      dt[, paste0(v,"_L",L) := shift(get(v), n=L, type="lag")]
    }
  }
  dt
}
dat <- make_lags(raw, vars=num_vars, pmax=p_max)
lag_vars <- grep("(_L[0-9]+)$", names(dat), value=TRUE)

4) 指標函數與儲存容器

說明：評分用 RMSE、MAE、MAD（MAD=median absolute error）；建立三個模型家族的暫存列表。

# ---- 指標 ----
rmse <- function(e) sqrt(mean(e^2, na.rm=TRUE))
mae  <- function(e) mean(abs(e), na.rm=TRUE)
madm <- function(e) median(abs(e), na.rm=TRUE)

# ---- 結果容器 ----
collect_ts    <- list()  # RandomWalk / AR / ARMA
collect_lasso <- list()  # Lasso / Ridge / EN / Ada* / PCR
collect_ml    <- list()  # Bagging Tree / Random Forest

5) 建立 rolling 切點並執行滾動預測

說明：每次只用最近 train_len 期做訓練，預測 h 期後的下一點。

# ---- rolling 索引 ----
n <- nrow(dat)
idx_train_end <- seq(from = train_len, to = n - h, by = 1)
if (!length(idx_train_end)) stop("資料太短，請降低 train_len 或確認資料長度。")

for (te in idx_train_end){
  train_idx <- (te - train_len + 1):te   # 固定長度視窗
  test_idx  <- te + h

  tr     <- dat[train_idx]
  te_row <- dat[test_idx]

  # ===== (A) 時間序列族：RW / AR / ARMA =====
  y_tr <- tr[[target_col]]
  y_te <- te_row[[target_col]]

  yhat_rw <- tail(y_tr, 1)

  fit_ar <- tryCatch(forecast::auto.arima(y_tr, d=0, max.q=0, seasonal=FALSE,
                                          stepwise=TRUE, approximation=TRUE), error=function(e) NULL)
  yhat_ar <- tryCatch(as.numeric(forecast::forecast(fit_ar, h=h)$mean[h]), error=function(e) yhat_rw)

  fit_arma <- tryCatch(forecast::auto.arima(y_tr, d=0, seasonal=FALSE,
                                            stepwise=TRUE, approximation=TRUE), error=function(e) NULL)
  yhat_arma <- tryCatch(as.numeric(forecast::forecast(fit_arma, h=h)$mean[h]), error=function(e) yhat_rw)

  collect_ts[[length(collect_ts)+1]] <- data.table(
    date = te_row[[time_col]],
    model = c("RandomWalk","AR","ARMA"),
    y     = y_te,
    yhat  = c(yhat_rw, yhat_ar, yhat_arma)
  )

  # ===== (B) Lasso 家族 + PCR =====
  X_tr <- as.matrix(tr[, ..lag_vars])
  X_te <- as.matrix(te_row[, ..lag_vars])
  good <- colSums(is.na(X_tr)) == 0
  X_tr <- X_tr[, good, drop=FALSE]
  X_te <- X_te[, good, drop=FALSE]
  y_tr2 <- y_tr
  
  # ---- 新增：過濾近零變異 + glmnet 公共參數 ----
  if(ncol(X_tr) > 0){
    sd_ok <- apply(X_tr, 2, sd, na.rm=TRUE) > 1e-8
    X_tr  <- X_tr[, sd_ok, drop=FALSE]
    X_te  <- X_te[, sd_ok, drop=FALSE]
  }
  
  glmnet_common <- list(
    standardize = TRUE,
    maxit = 1e6,              # 提高收斂上限
    nlambda = 80,             # 多一點 λ
    lambda.min.ratio = 1e-3   # 避免過小 λ
  )

  if (ncol(X_tr) >= 1){
    # Lasso
    fit_lasso <- tryCatch(
      do.call(cv.glmnet, c(list(x=X_tr, y=y_tr2, alpha=1), glmnet_common)),
      error=function(e) NULL
    )
    yhat_lasso <- if(!is.null(fit_lasso)) as.numeric(predict(fit_lasso, X_te, s="lambda.1se")) else NA_real_

    # Ridge
    fit_ridge <- tryCatch(
      do.call(cv.glmnet, c(list(x=X_tr, y=y_tr2, alpha=0), glmnet_common)),
      error=function(e) NULL
    )
    yhat_ridge <- if(!is.null(fit_ridge)) as.numeric(predict(fit_ridge, X_te, s="lambda.1se")) else NA_real_
    
    # Elastic Net
    fit_en <- tryCatch(
      do.call(cv.glmnet, c(list(x=X_tr, y=y_tr2, alpha=elastic_alpha), glmnet_common)),
      error=function(e) NULL
    )
    yhat_en <- if(!is.null(fit_en)) as.numeric(predict(fit_en, X_te, s="lambda.1se")) else NA_real_
    
    # Adaptive 權重（以 Ridge 初估）
    w_al <- tryCatch({
      fit0 <- do.call(cv.glmnet, c(list(x=X_tr, y=y_tr2, alpha=0), glmnet_common))
      b0   <- as.numeric(coef(fit0, s="lambda.1se"))[-1]
      w    <- 1/(abs(b0)^adapt_gamma + 1e-6)
      pmin(w, 1e4)   # 上限裁切，避免極端
    }, error=function(e) rep(1, ncol(X_tr)))
    
    fit_ada_lasso <- tryCatch(
      do.call(cv.glmnet, c(list(x=X_tr, y=y_tr2, alpha=1, penalty.factor=w_al), glmnet_common)),
      error=function(e) NULL
    )
    yhat_ada_lasso <- if(!is.null(fit_ada_lasso)) as.numeric(predict(fit_ada_lasso, X_te, s="lambda.1se")) else NA_real_
    
    fit_ada_en <- tryCatch(
      do.call(cv.glmnet, c(list(x=X_tr, y=y_tr2, alpha=elastic_alpha, penalty.factor=w_al), glmnet_common)),
      error=function(e) NULL
    )
    yhat_ada_en <- if(!is.null(fit_ada_en)) as.numeric(predict(fit_ada_en, X_te, s="lambda.1se")) else NA_real_

    # PCR（保留到 var_exp_cut）
    pca <- tryCatch(prcomp(X_tr, center=TRUE, scale.=TRUE), error=function(e) NULL)
    yhat_pcr <- NA_real_
    if(!is.null(pca)){
      var_exp <- cumsum(pca$sdev^2)/sum(pca$sdev^2)
      k <- which(var_exp >= var_exp_cut)[1]
      Z_tr <- pca$x[,1:k,drop=FALSE]
      Z_te <- scale(X_te, center=pca$center, scale=pca$scale) %*% pca$rotation[,1:k,drop=FALSE]
      fit_pcr <- tryCatch(lm(y_tr2 ~ Z_tr), error=function(e) NULL)
      if(!is.null(fit_pcr)) yhat_pcr <- as.numeric(cbind(1, Z_te) %*% coef(fit_pcr))
    }

    collect_lasso[[length(collect_lasso)+1]] <- data.table(
      date = te_row[[time_col]],
      model = c("Lasso","Ridge","ElasticNet","AdaLasso","AdaElasticNet","PCR"),
      y     = y_te,
      yhat  = c(yhat_lasso,yhat_ridge,yhat_en,yhat_ada_lasso,yhat_ada_en,yhat_pcr)
    )
  }

  # ===== (C) 非線性：Bagging Tree + Random Forest =====
  if (ncol(X_tr) >= 1){
    tr_df <- data.frame(y = y_tr2, X_tr)
    te_df <- data.frame(X_te)

    # Bagging
    fit_bag <- tryCatch(ipred::bagging(y ~ ., data=tr_df, coob=TRUE, nbagg=100),
                        error=function(e) NULL)
    yhat_bag <- if(!is.null(fit_bag)) as.numeric(predict(fit_bag, newdata=te_df)) else NA_real_

    # Random Forest（ranger）
    fit_rf <- tryCatch(
      ranger::ranger(y ~ ., data=tr_df,
                     num.trees=500,
                     mtry=max(1, floor(sqrt(ncol(X_tr)))),
                     min.node.size=5),
      error=function(e) NULL
    )
    yhat_rf <- if(!is.null(fit_rf)) as.numeric(predict(fit_rf, data=te_df)$predictions) else NA_real_

    collect_ml[[length(collect_ml)+1]] <- rbind(
      data.table(date = te_row[[time_col]], model="BaggingTree",  y=y_te, yhat=yhat_bag),
      data.table(date = te_row[[time_col]], model="RandomForest", y=y_te, yhat=yhat_rf)
    )
  }
}
# ---- rolling 索引 ----
n <- nrow(dat)
idx_train_end <- seq(from = train_len, to = n - h, by = 1)
if (!length(idx_train_end)) stop("資料太短，請降低 train_len 或確認資料長度。")

for (te in idx_train_end){
  train_idx <- (te - train_len + 1):te   # 固定長度視窗
  test_idx  <- te + h

  tr     <- dat[train_idx]
  te_row <- dat[test_idx]

  # ===== (A) 時間序列族：RW / AR / ARMA =====
  y_tr <- tr[[target_col]]
  y_te <- te_row[[target_col]]

  yhat_rw <- tail(y_tr, 1)

  fit_ar <- tryCatch(forecast::auto.arima(y_tr, d=0, max.q=0, seasonal=FALSE,
                                          stepwise=TRUE, approximation=TRUE), error=function(e) NULL)
  yhat_ar <- tryCatch(as.numeric(forecast::forecast(fit_ar, h=h)$mean[h]), error=function(e) yhat_rw)

  fit_arma <- tryCatch(forecast::auto.arima(y_tr, d=0, seasonal=FALSE,
                                            stepwise=TRUE, approximation=TRUE), error=function(e) NULL)
  yhat_arma <- tryCatch(as.numeric(forecast::forecast(fit_arma, h=h)$mean[h]), error=function(e) yhat_rw)

  collect_ts[[length(collect_ts)+1]] <- data.table(
    date = te_row[[time_col]],
    model = c("RandomWalk","AR","ARMA"),
    y     = y_te,
    yhat  = c(yhat_rw, yhat_ar, yhat_arma)
  )

  # ===== (B) Lasso 家族 + PCR =====
  X_tr <- as.matrix(tr[, ..lag_vars])
  X_te <- as.matrix(te_row[, ..lag_vars])
  good <- colSums(is.na(X_tr)) == 0
  X_tr <- X_tr[, good, drop=FALSE]
  X_te <- X_te[, good, drop=FALSE]
  y_tr2 <- y_tr

  if (ncol(X_tr) >= 1){
    # Lasso / Ridge / Elastic Net
    fit_lasso <- tryCatch(cv.glmnet(X_tr, y_tr2, alpha=1, standardize=TRUE), error=function(e) NULL)
    yhat_lasso <- if(!is.null(fit_lasso)) as.numeric(predict(fit_lasso, X_te, s="lambda.min")) else NA_real_

    fit_ridge <- tryCatch(cv.glmnet(X_tr, y_tr2, alpha=0, standardize=TRUE), error=function(e) NULL)
    yhat_ridge <- if(!is.null(fit_ridge)) as.numeric(predict(fit_ridge, X_te, s="lambda.min")) else NA_real_

    fit_en <- tryCatch(cv.glmnet(X_tr, y_tr2, alpha=elastic_alpha, standardize=TRUE), error=function(e) NULL)
    yhat_en <- if(!is.null(fit_en)) as.numeric(predict(fit_en, X_te, s="lambda.min")) else NA_real_

    # Adaptive 權重（以 Ridge 初估）
    w_al <- tryCatch({
      fit0 <- cv.glmnet(X_tr, y_tr2, alpha=0, standardize=TRUE)
      b0   <- as.numeric(coef(fit0, s="lambda.min"))[-1]
      1/(abs(b0)^adapt_gamma + 1e-6)
    }, error=function(e) rep(1, ncol(X_tr)))

    fit_ada_lasso <- tryCatch(cv.glmnet(X_tr, y_tr2, alpha=1, penalty.factor=w_al, standardize=TRUE), error=function(e) NULL)
    yhat_ada_lasso <- if(!is.null(fit_ada_lasso)) as.numeric(predict(fit_ada_lasso, X_te, s="lambda.min")) else NA_real_

    fit_ada_en <- tryCatch(cv.glmnet(X_tr, y_tr2, alpha=elastic_alpha, penalty.factor=w_al, standardize=TRUE), error=function(e) NULL)
    yhat_ada_en <- if(!is.null(fit_ada_en)) as.numeric(predict(fit_ada_en, X_te, s="lambda.min")) else NA_real_

    # PCR（保留到 var_exp_cut）
    pca <- tryCatch(prcomp(X_tr, center=TRUE, scale.=TRUE), error=function(e) NULL)
    yhat_pcr <- NA_real_
    if(!is.null(pca)){
      var_exp <- cumsum(pca$sdev^2)/sum(pca$sdev^2)
      k <- which(var_exp >= var_exp_cut)[1]
      Z_tr <- pca$x[,1:k,drop=FALSE]
      Z_te <- scale(X_te, center=pca$center, scale=pca$scale) %*% pca$rotation[,1:k,drop=FALSE]
      fit_pcr <- tryCatch(lm(y_tr2 ~ Z_tr), error=function(e) NULL)
      if(!is.null(fit_pcr)) yhat_pcr <- as.numeric(cbind(1, Z_te) %*% coef(fit_pcr))
    }

    collect_lasso[[length(collect_lasso)+1]] <- data.table(
      date = te_row[[time_col]],
      model = c("Lasso","Ridge","ElasticNet","AdaLasso","AdaElasticNet","PCR"),
      y     = y_te,
      yhat  = c(yhat_lasso,yhat_ridge,yhat_en,yhat_ada_lasso,yhat_ada_en,yhat_pcr)
    )
  }

  # ===== (C) 非線性：Bagging Tree + Random Forest =====
  if (ncol(X_tr) >= 1){
    tr_df <- data.frame(y = y_tr2, X_tr)
    te_df <- data.frame(X_te)

    # Bagging
    fit_bag <- tryCatch(ipred::bagging(y ~ ., data=tr_df, coob=TRUE, nbagg=100),
                        error=function(e) NULL)
    yhat_bag <- if(!is.null(fit_bag)) as.numeric(predict(fit_bag, newdata=te_df)) else NA_real_

    # Random Forest（ranger）
    fit_rf <- tryCatch(
      ranger::ranger(y ~ ., data=tr_df,
                     num.trees=500,
                     mtry=max(1, floor(sqrt(ncol(X_tr)))),
                     min.node.size=5),
      error=function(e) NULL
    )
    yhat_rf <- if(!is.null(fit_rf)) as.numeric(predict(fit_rf, data=te_df)$predictions) else NA_real_

    collect_ml[[length(collect_ml)+1]] <- rbind(
      data.table(date = te_row[[time_col]], model="BaggingTree",  y=y_te, yhat=yhat_bag),
      data.table(date = te_row[[time_col]], model="RandomForest", y=y_te, yhat=yhat_rf)
    )
  }
}

6) 彙整、評分表、輸出 CSV

說明：把所有樣本外預測合併，計算 RMSE/MAE/MAD，輸出表。

# ---- 合併滾動預測 ----
roll_ts    <- rbindlist(collect_ts,    fill=TRUE)
roll_lasso <- rbindlist(collect_lasso, fill=TRUE)
roll_ml    <- rbindlist(collect_ml,    fill=TRUE)

metrics <- rbindlist(list(roll_ts, roll_lasso, roll_ml), fill=TRUE)
metrics <- metrics[!is.na(yhat)]

score <- metrics[, .(
  RMSE = rmse(y - yhat),
  MAE  = mae(y - yhat),
  MAD  = madm(y - yhat)
), by=.(model)]


# 建立輸出資料夾：gen_excel_files/P6_彙整_評分表_輸出/當日日期
out_root <- "generated_excel_files/P6_彙整_評分表_輸出"
out_dir  <- file.path(out_root, format(Sys.Date(), "%Y-%m-%d"))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

# 分別輸出兩個 xlsx 檔
writexl::write_xlsx(score,   path = file.path(out_dir, "model_metrics_rmse_mae_mad.xlsx"))
writexl::write_xlsx(metrics, path = file.path(out_dir, "oos_predictions_by_model.xlsx"))

# 合併兩個 Sheet 
writexl::write_xlsx(
  list(
    model_metrics_rmse_mae_mad = score,
    oos_predictions_by_model   = metrics
  ),
  path = file.path(out_dir, "forecast_outputs.xlsx")
)

score

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

model	RMSE	MAE	MAD
RandomWalk	154157.9	107254.10	73629.11
AR	124158.5	81819.52	52235.10
ARMA	125746.4	84037.73	53110.68
Lasso	112593.3	64145.27	36259.37
Ridge	109013.4	69685.52	44931.30
ElasticNet	104431.6	60678.63	36190.86
AdaLasso	120591.9	71403.15	43325.47
AdaElasticNet	130447.3	71171.08	42471.68
PCR	106319.7	62967.89	40404.51
BaggingTree	103983.6	67013.72	44623.75
RandomForest	114332.3	81662.60	60914.07

7) 三張樣本外走勢圖

說明：一張時間序列家族、一張 Lasso 家族、一張非線性（Bagging+RF）。

# ---- 作圖 ----
plot_series <- function(dt, title_str){
  ggplot(dt, aes(x=date)) +
    geom_line(aes(y=y), linewidth=0.6) +
    geom_line(aes(y=yhat, color=model), linewidth=0.6) +
    scale_x_date(date_labels="%Y-%m") +
    labs(x=NULL, y=target_col, title=title_str) +
    theme_minimal()
}
p1 <- plot_series(roll_ts,    "Time Series Model (RW/AR/ARMA) OOS")
p2 <- plot_series(roll_lasso, "Lasso Family (Ridge/EN/Adaptive/PCR) OOS")
p3 <- plot_series(roll_ml,    "Nonlinear ML models (Bagging & Random Forest) OOS")

out_dir <- file.path("generated_figures/三張樣本外走勢圖", format(Sys.Date(), "%Y-%m-%d"))
if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE)

ggsave(file.path(out_dir, "plot_timeseries.png"), p1, width = 22, height = 7, units = "in", dpi = 600)
ggsave(file.path(out_dir, "plot_lasso.png"),      p2, width = 22, height = 7, units = "in", dpi = 600)
ggsave(file.path(out_dir, "plot_ml.png"),         p3, width = 22, height = 7, units = "in", dpi = 600)

p1; p2; p3

8) Nowcast「最後一期再多推 h 期」

說明：仍用 rolling 的最後一個視窗來擬合，再預測 last + h；並輸出各模型的單點預測值。

8a) 參數：一次預測未來 H_out 個月

# ---- 一次多步預測的步數 ----
H_out <- 12   # 想要幾個月就寫幾個月（例：12 = 一整年）

8b) 準備最後一個 rolling 訓練窗 + 預測用日期索引

說明：仍用 rolling 的最後視窗來擬合；預測日期用月加法產生。

# ---- 取 rolling 的最後視窗 ----
last_time <- max(dat[[time_col]])
end_idx   <- nrow(dat)
start_idx <- max(1, end_idx - train_len + 1)
final_tr  <- dat[start_idx:end_idx]

# ---- 目標向量與特徵矩陣（同前）----
y_all   <- final_tr[[target_col]]
X_all   <- as.matrix(final_tr[, ..lag_vars])
good    <- colSums(is.na(X_all)) == 0
X_all   <- X_all[, good, drop=FALSE]

# 預測當下可用的特徵（來自全資料最後一列的 lag）
last_row <- tail(dat, 1)
x_origin <- as.matrix(last_row[, ..lag_vars])[, good, drop=FALSE]

# 預測的日期（t+1 ... t+H_out）
future_dates <- sapply(1:H_out, function(hh) last_time %m+% months(hh))

8c) 時間序列模型：一次 12 步

說明：RW/AR/ARMA 可以一口氣輸出 h=H_out 的路徑。

# ---- 時間序列：一次多步 ----
next_rw_path <- rep(tail(y_all,1), H_out)

fit_ar_all <- tryCatch(forecast::auto.arima(y_all, d=0, max.q=0, seasonal=FALSE), error=function(e) NULL)
next_ar_path <- tryCatch(as.numeric(forecast::forecast(fit_ar_all, h=H_out)$mean), error=function(e) next_rw_path)

fit_arma_all <- tryCatch(forecast::auto.arima(y_all, d=0, seasonal=FALSE), error=function(e) NULL)
next_arma_path <- tryCatch(as.numeric(forecast::forecast(fit_arma_all, h=H_out)$mean), error=function(e) next_rw_path)

extra_ts <- data.table(
  date  = as.Date(future_dates, origin = "1970-01-01"),
  horizon = 1:H_out,
  model = rep(c("RandomWalk","AR","ARMA"), each = H_out),
  yhat  = c(next_rw_path, next_ar_path, next_arma_path)
)

8d) Lasso 家族 + PCR：direct strategy（每個 h 一個模型）

說明：對每個 h，把目標改成「h 步之後的 y」（lead h），用當下的 lag 特徵訓練，然後用 x_origin 預測 t+h。這樣不需要知道未來 X 的路徑，也避免遞迴誤差累積。

# ---- 公共設定（穩定收斂）----
glmnet_common <- list(standardize=TRUE, maxit=1e6, nlambda=80, lambda.min.ratio=1e-3)

# ---- 建 H_out 組 training pairs： (X_all_trim, y_h) ----
train_list <- lapply(1:H_out, function(hh){
  # y 需能往前對齊：去掉最後 hh 期
  y_h  <- shift(y_all, type="lead", n=hh)
  keep <- !is.na(y_h)
  list(y = y_h[keep], X = X_all[keep,, drop=FALSE], h = hh)
})

# ---- 對每個 h 訓練與預測 ----
res_lasso <- list()
for (item in train_list){
  y_tr2 <- item$y
  X_tr2 <- item$X
  hh    <- item$h

  # 過濾近零變異
  sd_ok <- apply(X_tr2, 2, sd, na.rm=TRUE) > 1e-8
  X_tr2 <- X_tr2[, sd_ok, drop=FALSE]
  x0    <- x_origin[, sd_ok, drop=FALSE]

  # 如果沒有特徵就跳過
  if (ncol(X_tr2) < 1){
    res_lasso[[length(res_lasso)+1]] <- data.table(
      date = future_dates[hh], horizon=hh,
      model = c("Lasso","Ridge","ElasticNet","AdaLasso","AdaElasticNet","PCR"),
      yhat  = NA_real_)
    next
  }

  # Lasso / Ridge / EN
  fit_lasso <- tryCatch(do.call(cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=1), glmnet_common)), error=function(e) NULL)
  yhat_lasso <- if(!is.null(fit_lasso)) as.numeric(predict(fit_lasso, x0, s="lambda.1se")) else NA_real_

  fit_ridge <- tryCatch(do.call(cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=0), glmnet_common)), error=function(e) NULL)
  yhat_ridge <- if(!is.null(fit_ridge)) as.numeric(predict(fit_ridge, x0, s="lambda.1se")) else NA_real_

  fit_en <- tryCatch(do.call(cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=elastic_alpha), glmnet_common)), error=function(e) NULL)
  yhat_en <- if(!is.null(fit_en)) as.numeric(predict(fit_en, x0, s="lambda.1se")) else NA_real_

  # Adaptive 權重（ridge 初估 + 上限裁切）
  w_al <- tryCatch({
    fit0 <- do.call(cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=0), glmnet_common))
    b0   <- as.numeric(coef(fit0, s="lambda.1se"))[-1]
    w    <- 1/(abs(b0)^adapt_gamma + 1e-6)
    pmin(w, 1e4)
  }, error=function(e) rep(1, ncol(X_tr2)))

  fit_ada_lasso <- tryCatch(do.call(cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=1, penalty.factor=w_al), glmnet_common)), error=function(e) NULL)
  yhat_ada_lasso <- if(!is.null(fit_ada_lasso)) as.numeric(predict(fit_ada_lasso, x0, s="lambda.1se")) else NA_real_

  fit_ada_en <- tryCatch(do.call(cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=elastic_alpha, penalty.factor=w_al), glmnet_common)), error=function(e) NULL)
  yhat_ada_en <- if(!is.null(fit_ada_en)) as.numeric(predict(fit_ada_en, x0, s="lambda.1se")) else NA_real_

  # PCR
  yhat_pcr <- NA_real_
  pca <- tryCatch(prcomp(X_tr2, center=TRUE, scale.=TRUE), error=function(e) NULL)
  if(!is.null(pca)){
    var_exp <- cumsum(pca$sdev^2)/sum(pca$sdev^2)
    k <- which(var_exp >= var_exp_cut)[1]
    Z_tr <- pca$x[,1:k,drop=FALSE]
    fit_pcr <- tryCatch(lm(y_tr2 ~ Z_tr), error=function(e) NULL)
    if(!is.null(fit_pcr)){
      z0 <- scale(x0, center=pca$center, scale=pca$scale) %*% pca$rotation[,1:k,drop=FALSE]
      yhat_pcr <- as.numeric(cbind(1, z0) %*% coef(fit_pcr))
    }
  }

  res_lasso[[length(res_lasso)+1]] <- data.table(
    date = future_dates[hh], horizon = hh,
    model = c("Lasso","Ridge","ElasticNet","AdaLasso","AdaElasticNet","PCR"),
    yhat  = c(yhat_lasso,yhat_ridge,yhat_en,yhat_ada_lasso,yhat_ada_en,yhat_pcr)
  )
}
extra_lasso <- rbindlist(res_lasso, fill=TRUE)

8e) 非線性：Bagging Tree + Random Forest（direct strategy）

說明：跟上面相同概念，每個 h 一個模型；特徵是 t 時點的 lag。

res_ml <- list()
for (item in train_list){
  y_tr2 <- item$y
  X_tr2 <- item$X
  hh    <- item$h

  # 過濾近零變異
  sd_ok <- apply(X_tr2, 2, sd, na.rm=TRUE) > 1e-8
  X_tr2 <- X_tr2[, sd_ok, drop=FALSE]
  x0    <- x_origin[, sd_ok, drop=FALSE]

  if (ncol(X_tr2) < 1){
    res_ml[[length(res_ml)+1]] <- data.table(
      date=future_dates[hh], horizon=hh,
      model=c("BaggingTree","RandomForest"),
      yhat=c(NA_real_, NA_real_))
    next
  }

  tr_df <- data.frame(y=y_tr2, X_tr2)
  te_df <- data.frame(x0)

  # Bagging
  fit_bag <- tryCatch(ipred::bagging(y ~ ., data=tr_df, coob=FALSE, nbagg=100), error=function(e) NULL)
  yhat_bag <- if(!is.null(fit_bag)) as.numeric(predict(fit_bag, newdata=te_df)) else NA_real_

  # Random Forest
  fit_rf <- tryCatch(
    ranger::ranger(y ~ ., data=tr_df, num.trees=500,
                   mtry=max(1, floor(sqrt(ncol(X_tr2)))), min.node.size=5),
    error=function(e) NULL
  )
  yhat_rf <- if(!is.null(fit_rf)) as.numeric(predict(fit_rf, data=te_df)$predictions) else NA_real_

  res_ml[[length(res_ml)+1]] <- data.table(
    date=future_dates[hh], horizon=hh,
    model=c("BaggingTree","RandomForest"),
    yhat=c(yhat_bag,yhat_rf)
  )
}
extra_ml <- rbindlist(res_ml, fill=TRUE)

8f) 整併輸出（一次包含 1..H_out 步）

# 排序不變
extra_all <- rbindlist(list(extra_ts, extra_lasso, extra_ml), fill=TRUE)
setorderv(extra_all, c("horizon","model"))

# 輸出到 generated_excel_files/P8_Nowcast/YYYY-MM-DD/extra_multi_step_forecasts_H{H_out}.xlsx
if (!requireNamespace("writexl", quietly = TRUE)) install.packages("writexl")

out_root <- file.path("generated_excel_files/P8_Nowcast")
out_dir  <- file.path(out_root, format(Sys.Date(), "%Y-%m-%d"))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

writexl::write_xlsx(
  extra_all,
  path = file.path(out_dir, sprintf("extra_multi_step_forecasts_H%d.xlsx", H_out))
)

extra_all

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

date	horizon	model	yhat
2025-09-01	1	AR	250525.95
2025-09-01	1	ARMA	216171.96
2025-09-01	1	AdaElasticNet	311322.72
2025-09-01	1	AdaLasso	311394.35
2025-09-01	1	BaggingTree	311695.53
2025-09-01	1	ElasticNet	311391.84
2025-09-01	1	Lasso	276140.08
2025-09-01	1	PCR	359919.03
2025-09-01	1	RandomForest	401311.40
2025-09-01	1	RandomWalk	153710.02
2025-09-01	1	Ridge	295751.00
2025-10-01	2	AR	271017.53
2025-10-01	2	ARMA	276167.50
2025-10-01	2	AdaElasticNet	280886.50
2025-10-01	2	AdaLasso	283950.99
2025-10-01	2	BaggingTree	327087.31
2025-10-01	2	ElasticNet	289208.32
2025-10-01	2	Lasso	291638.56
2025-10-01	2	PCR	261090.93
2025-10-01	2	RandomForest	412557.74
2025-10-01	2	RandomWalk	153710.02
2025-10-01	2	Ridge	300581.36
2025-11-01	3	AR	275354.67
2025-11-01	3	ARMA	276167.50
2025-11-01	3	AdaElasticNet	287740.47
2025-11-01	3	AdaLasso	285390.30
2025-11-01	3	BaggingTree	295354.59
2025-11-01	3	ElasticNet	288774.03
2025-11-01	3	Lasso	295737.05
2025-11-01	3	PCR	346900.20
2025-11-01	3	RandomForest	431130.33
2025-11-01	3	RandomWalk	153710.02
2025-11-01	3	Ridge	320669.30
2025-12-01	4	AR	276272.65
2025-12-01	4	ARMA	276167.50
2025-12-01	4	AdaElasticNet	181891.08
2025-12-01	4	AdaLasso	173956.73
2025-12-01	4	BaggingTree	237168.54
2025-12-01	4	ElasticNet	247146.59
2025-12-01	4	Lasso	241773.20
2025-12-01	4	PCR	220852.37
2025-12-01	4	RandomForest	395123.46
2025-12-01	4	RandomWalk	153710.02
2025-12-01	4	Ridge	304274.49
2026-01-01	5	AR	276466.95
2026-01-01	5	ARMA	276167.50
2026-01-01	5	AdaElasticNet	306531.56
2026-01-01	5	AdaLasso	300510.85
2026-01-01	5	BaggingTree	298958.86
2026-01-01	5	ElasticNet	282541.63
2026-01-01	5	Lasso	302399.81
2026-01-01	5	PCR	353849.04
2026-01-01	5	RandomForest	435389.10
2026-01-01	5	RandomWalk	153710.02
2026-01-01	5	Ridge	354078.69
2026-02-01	6	AR	276508.07
2026-02-01	6	ARMA	276167.50
2026-02-01	6	AdaElasticNet	274153.57
2026-02-01	6	AdaLasso	279686.79
2026-02-01	6	BaggingTree	175807.27
2026-02-01	6	ElasticNet	269860.30
2026-02-01	6	Lasso	260875.68
2026-02-01	6	PCR	153710.34
2026-02-01	6	RandomForest	465582.25
2026-02-01	6	RandomWalk	153710.02
2026-02-01	6	Ridge	315384.61
2026-03-01	7	AR	276516.77
2026-03-01	7	ARMA	276167.50
2026-03-01	7	AdaElasticNet	277929.06
2026-03-01	7	AdaLasso	265252.30
2026-03-01	7	BaggingTree	327321.06
2026-03-01	7	ElasticNet	282349.89
2026-03-01	7	Lasso	281967.82
2026-03-01	7	PCR	301926.20
2026-03-01	7	RandomForest	467540.40
2026-03-01	7	RandomWalk	153710.02
2026-03-01	7	Ridge	345250.04
2026-04-01	8	AR	276518.62
2026-04-01	8	ARMA	276167.50
2026-04-01	8	AdaElasticNet	283794.68
2026-04-01	8	AdaLasso	268351.07
2026-04-01	8	BaggingTree	216947.55
2026-04-01	8	ElasticNet	271526.00
2026-04-01	8	Lasso	268351.07
2026-04-01	8	PCR	252561.11
2026-04-01	8	RandomForest	485590.68
2026-04-01	8	RandomWalk	153710.02
2026-04-01	8	Ridge	292230.46
2026-05-01	9	AR	276519.01
2026-05-01	9	ARMA	276167.50
2026-05-01	9	AdaElasticNet	292211.07
2026-05-01	9	AdaLasso	297202.97
2026-05-01	9	BaggingTree	522624.00
2026-05-01	9	ElasticNet	293842.67
2026-05-01	9	Lasso	306885.18
2026-05-01	9	PCR	448977.64
2026-05-01	9	RandomForest	541002.13
2026-05-01	9	RandomWalk	153710.02
2026-05-01	9	Ridge	292211.07
2026-06-01	10	AR	276519.09
2026-06-01	10	ARMA	276167.50
2026-06-01	10	AdaElasticNet	291191.82
2026-06-01	10	AdaLasso	306402.28
2026-06-01	10	BaggingTree	565933.06
2026-06-01	10	ElasticNet	291191.82
2026-06-01	10	Lasso	306402.28
2026-06-01	10	PCR	413978.72
2026-06-01	10	RandomForest	554060.76
2026-06-01	10	RandomWalk	153710.02
2026-06-01	10	Ridge	291191.82
2026-07-01	11	AR	276519.11
2026-07-01	11	ARMA	276167.50
2026-07-01	11	AdaElasticNet	266201.01
2026-07-01	11	AdaLasso	274136.15
2026-07-01	11	BaggingTree	461237.15
2026-07-01	11	ElasticNet	326400.58
2026-07-01	11	Lasso	311730.00
2026-07-01	11	PCR	409758.29
2026-07-01	11	RandomForest	477473.68
2026-07-01	11	RandomWalk	153710.02
2026-07-01	11	Ridge	372981.39
2026-08-01	12	AR	276519.11
2026-08-01	12	ARMA	276167.50
2026-08-01	12	AdaElasticNet	185341.00
2026-08-01	12	AdaLasso	288748.68
2026-08-01	12	BaggingTree	209415.06
2026-08-01	12	ElasticNet	258509.77
2026-08-01	12	Lasso	238416.35
2026-08-01	12	PCR	97408.65
2026-08-01	12	RandomForest	433694.84
2026-08-01	12	RandomWalk	153710.02
2026-08-01	12	Ridge	284697.64

9) 重要變數

9a) 參數與小工具（可調整）

# ── 重要性參數（可調） ─────────────────────────────────────────────
h               <- 1          # 這一輪先做 h=1（預測下一期）
k_top           <- 8          # 取前幾個重要變數
var_exp_cut     <- 0.95       # PCR 累積解釋比例
elastic_alpha   <- 0.5        # Elastic Net 的 alpha
adapt_gamma     <- 1.0        # Adaptive 權重 1/|beta|^gamma
agg_method      <- "sum"      # 彙總到原始變數時的方式："sum" 或 "max"

# 名稱對照（meat_tab.xlsx）設定
map_path        <- "/Users/traviscua/Library/CloudStorage/OneDrive-個人/Documents/1. R/Research/Oct_2025/scratch_3/existing_excel_files/meta_tab.xlsx"
map_sheet       <- 1
short_max_chars <- 28
short_use_dataset_prefix <- TRUE
dataset_prefix_nchar     <- 6

# ── glmnet 共用設定（提升收斂穩定） ───────────────────────────────
glmnet_common <- list(
  standardize = TRUE,
  maxit = 1e6,
  nlambda = 80,
  lambda.min.ratio = 1e-3
)

# ── 輔助：取前 k 名（依重要性排序的名稱） ─────────────────────────
topk_names <- function(imp_vec, k = k_top){
  imp_vec <- imp_vec[is.finite(imp_vec)]
  if(length(imp_vec)==0) return(rep(NA_character_, k))
  ord <- order(imp_vec, decreasing=TRUE)
  picks <- names(imp_vec)[ord][seq_len(min(k, length(ord)))]
  out <- rep(NA_character_, k)
  out[seq_along(picks)] <- picks
  out
}

9b) 建最後一個 rolling 視窗資料與 direct strategy 的訓練組

# 取得最後一個 rolling 視窗
end_idx   <- nrow(dat)
start_idx <- max(1, end_idx - train_len + 1)
final_tr  <- dat[start_idx:end_idx]

# 目標與特徵（lag 特徵）
y_all <- final_tr[[target_col]]
X_all <- as.matrix(final_tr[, ..lag_vars])

# 去掉含 NA 的欄、近零變異欄
good_na <- colSums(is.na(X_all)) == 0
X_all   <- X_all[, good_na, drop=FALSE]
if(ncol(X_all) == 0) stop("此視窗內沒有可用的 lag 特徵。")
sd_ok   <- apply(X_all, 2, sd, na.rm=TRUE) > 1e-8
X_all   <- X_all[, sd_ok, drop=FALSE]

# direct strategy：用 (X_t, y_{t+h})
y_lead <- data.table::shift(y_all, type="lead", n=h)
keep   <- !is.na(y_lead)
X_tr2  <- X_all[keep,, drop=FALSE]
y_tr2  <- y_lead[keep]

9c) Lasso / Ridge / EN / Adaptive / Adaptive EN 的重要性（|β|）

# Lasso
fit_lasso <- tryCatch(do.call(glmnet::cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=1), glmnet_common)), error=function(e) NULL)
imp_lasso <- if(!is.null(fit_lasso)) abs(as.numeric(coef(fit_lasso, s="lambda.1se"))[-1]) else setNames(numeric(0), character(0))
names(imp_lasso) <- colnames(X_tr2)

# Ridge
fit_ridge <- tryCatch(do.call(glmnet::cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=0), glmnet_common)), error=function(e) NULL)
imp_ridge <- if(!is.null(fit_ridge)) abs(as.numeric(coef(fit_ridge, s="lambda.1se"))[-1]) else setNames(numeric(0), character(0))
names(imp_ridge) <- colnames(X_tr2)

# Elastic Net
fit_en <- tryCatch(do.call(glmnet::cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=elastic_alpha), glmnet_common)), error=function(e) NULL)
imp_en <- if(!is.null(fit_en)) abs(as.numeric(coef(fit_en, s="lambda.1se"))[-1]) else setNames(numeric(0), character(0))
names(imp_en) <- colnames(X_tr2)

# Adaptive 權重（Ridge 初估 + 上限裁切）
w_al <- tryCatch({
  fit0 <- do.call(glmnet::cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=0), glmnet_common))
  b0   <- as.numeric(coef(fit0, s="lambda.1se"))[-1]
  pmin(1/(abs(b0)^adapt_gamma + 1e-6), 1e4)
}, error=function(e) rep(1, ncol(X_tr2)))

fit_ada_lasso <- tryCatch(do.call(glmnet::cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=1, penalty.factor=w_al), glmnet_common)), error=function(e) NULL)
imp_ada_lasso <- if(!is.null(fit_ada_lasso)) abs(as.numeric(coef(fit_ada_lasso, s="lambda.1se"))[-1]) else setNames(numeric(0), character(0))
names(imp_ada_lasso) <- colnames(X_tr2)

fit_ada_en <- tryCatch(do.call(glmnet::cv.glmnet, c(list(x=X_tr2, y=y_tr2, alpha=elastic_alpha, penalty.factor=w_al), glmnet_common)), error=function(e) NULL)
imp_ada_en <- if(!is.null(fit_ada_en)) abs(as.numeric(coef(fit_ada_en, s="lambda.1se"))[-1]) else setNames(numeric(0), character(0))
names(imp_ada_en) <- colnames(X_tr2)

9d) PCR 的重要性（loading × 迴歸係數的加權和）

imp_pcr <- setNames(numeric(0), character(0))
pca <- tryCatch(prcomp(X_tr2, center=TRUE, scale.=TRUE), error=function(e) NULL)
if(!is.null(pca)){
  var_exp <- cumsum(pca$sdev^2)/sum(pca$sdev^2)
  k_pc <- which(var_exp >= var_exp_cut)[1]
  Z_tr <- pca$x[,1:k_pc, drop=FALSE]
  fit_pcr <- tryCatch(lm(y_tr2 ~ Z_tr), error=function(e) NULL)
  if(!is.null(fit_pcr)){
    beta <- coef(fit_pcr)[-1]
    rot  <- pca$rotation[,1:k_pc, drop=FALSE]   # loading: var × PC
    sc   <- as.numeric(abs(rot %*% abs(beta)))  # 變數層級加權
    imp_pcr <- setNames(sc, rownames(rot))
  }
}

9e) 樹模型重要性：Bagging（proxy）與 Random Forest（Permutation importance）

library(ranger)
p <- ncol(X_tr2)

# Bagging 的 proxy：ranger + mtry = p (不抽樣特徵) + permutation importance
fit_bag_imp <- tryCatch(
  ranger(y = y_tr2, x = data.frame(X_tr2),
         num.trees = 500, mtry = p, min.node.size = 5,
         importance = "permutation", seed = 123),
  error=function(e) NULL
)
imp_bag <- if(!is.null(fit_bag_imp)) fit_bag_imp$variable.importance else setNames(numeric(0), character(0))

# Random Forest：mtry = sqrt(p)
fit_rf_imp <- tryCatch(
  ranger(y = y_tr2, x = data.frame(X_tr2),
         num.trees = 500, mtry = max(1, floor(sqrt(p))), min.node.size = 5,
         importance = "permutation", seed = 123),
  error=function(e) NULL
)
imp_rf <- if(!is.null(fit_rf_imp)) fit_rf_imp$variable.importance else setNames(numeric(0), character(0))

9e)【Lag 層級】組成前 k 名表（代碼名）

library(data.table)

imp_list_lag <- list(
  "LASSO"            = imp_lasso,
  "EN"               = imp_en,
  "Adaptive LASSO"   = imp_ada_lasso,
  "Adaptive EN"      = imp_ada_en,
  "Ridge"            = imp_ridge,
  "PCR"              = imp_pcr,
  "Bagging (proxy)"  = imp_bag,
  "Random Forest"    = imp_rf
)

rank_vec <- as.character(seq_len(k_top))
tbl_lag_code <- data.table(排序 = rank_vec)
for (nm in names(imp_list_lag)){
  tbl_lag_code[[nm]] <- topk_names(imp_list_lag[[nm]], k_top)
}

out_root <- file.path("generated_excel_files/P9_重要變數/9e")
out_dir  <- file.path(out_root, format(Sys.Date(), "%Y-%m-%d"))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

writexl::write_xlsx(
  tbl_lag_code,
  path = file.path(out_dir, sprintf("variable_importance_LAG_top%d_h%d_CODE.xlsx", k_top, h))
)

9f) 載入 `meta_tab.xlsx` 並把 Lag 層級表映射成短/長名

# 讀對照表
map_raw <- as.data.table(readxl::read_excel(map_path, sheet = map_sheet))
setnames(map_raw, old=names(map_raw), new=make.names(names(map_raw)))
stopifnot(all(c("var_label","dataset_title","series_name") %in% names(map_raw)))
setorder(map_raw, var_label)
if (any(duplicated(map_raw$var_label))){
  warning("var_label 在對照表中有重複，將取第一筆。")
  map_raw <- map_raw[!duplicated(var_label)]
}

# 剝 lag 後綴
split_lag <- function(x){
  lag <- ifelse(grepl("_L\\d+$", x), sub(".*_L(\\d+)$","L\\1", x), NA_character_)
  base <- ifelse(is.na(lag), x, sub("_L\\d+$","", x))
  list(base=base, lag=lag)
}

# 產生長/短名
make_display_names <- function(base, lag, mapDT,
                               short_max=short_max_chars,
                               use_prefix=short_use_dataset_prefix,
                               prefix_n=dataset_prefix_nchar){
  m <- mapDT[.(base), on=.(var_label)]
  # 長名
  long <- ifelse(!is.na(m$dataset_title) & !is.na(m$series_name),
                 paste0(m$dataset_title, "｜", m$series_name),
                 base)
  long <- ifelse(!is.na(lag), paste0(long, " (", lag, ")"), long)

  # 短名
  short_core <- ifelse(!is.na(m$series_name), m$series_name, base)
  if (use_prefix){
    pref <- ifelse(!is.na(m$dataset_title), substr(m$dataset_title, 1L, prefix_n), "")
    short_core <- ifelse(pref!="", paste0(pref, "｜", short_core), short_core)
  }
  short_core <- ifelse(!is.na(lag), paste0(short_core, " (", lag, ")"), short_core)
  short <- ifelse(nchar(short_core) > short_max,
                  paste0(substr(short_core, 1L, short_max-1L), "…"),
                  short_core)
  list(long=long, short=short)
}

map_var_labels <- function(v, mapDT){
  sp <- split_lag(v)
  nm <- make_display_names(sp$base, sp$lag, mapDT)
  data.table(var_label=v, base=sp$base, lag=sp$lag,
             long_name=nm$long, short_name=nm$short)
}

# 套用到 lag 層級表
model_cols <- setdiff(names(tbl_lag_code), "排序")
all_vars <- unique(unlist(lapply(tbl_lag_code[, ..model_cols], unlist), use.names=FALSE))
all_vars <- all_vars[!is.na(all_vars) & all_vars!=""]
dict_lag <- if(length(all_vars)) map_var_labels(all_vars, map_raw) else data.table()

tbl_lag_short <- copy(tbl_lag_code)
tbl_lag_long  <- copy(tbl_lag_code)
for (mc in model_cols){
  m <- merge(data.table(var_label = tbl_lag_code[[mc]]), dict_lag, by="var_label", all.x=TRUE)
  tbl_lag_short[[mc]] <- ifelse(is.na(m$short_name), paste0(m$var_label, " (未對應)"), m$short_name)
  tbl_lag_long[[mc]]  <- ifelse(is.na(m$long_name),  paste0(m$var_label, " (未對應)"), m$long_name)
}


out_root <- file.path("generated_excel_files", "P9_重要變數", "9f")
out_dir  <- file.path(out_root, format(Sys.Date(), "%Y-%m-%d"))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

# 三個檔各自輸出
writexl::write_xlsx(
  dict_lag[, .(var_label, base, lag, short_name, long_name)],
  path = file.path(out_dir, sprintf("var_name_dictionary_LAG_h%d.xlsx", h))
)

writexl::write_xlsx(
  tbl_lag_short,
  path = file.path(out_dir, sprintf("variable_importance_LAG_top%d_h%d_SHORT.xlsx", k_top, h))
)

writexl::write_xlsx(
  tbl_lag_long,
  path = file.path(out_dir, sprintf("variable_importance_LAG_top%d_h%d_LONG.xlsx", k_top, h))
)

writexl::write_xlsx(
  list(
    var_name_dictionary_LAG = dict_lag[, .(var_label, base, lag, short_name, long_name)],
    LAG_top_SHORT           = tbl_lag_short,
    LAG_top_LONG            = tbl_lag_long
  ),
  path = file.path(out_dir, sprintf("P9_重要變數_9f_top%d_h%d.xlsx", k_top, h))
)

tbl_lag_short

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

排序	LASSO	EN	Adaptive LASSO	Adaptive EN	Ridge	PCR	Bagging (proxy)	Random Forest
1	消費者物價基｜總指數 (L1)	受僱員工每人｜工業及服務業 \| 合計 (L3)	消費者物價基｜總指數 (L1)	消費者物價基｜總指數 (L1)	人力資源主要｜勞動力參與率(%) \| 合計 (L1)	每人每月總薪｜工業及服務業 \| 合計 (L4)	每人每月總薪｜工業及服務業 \| 合計 (L4)	每人每月總薪｜工業及服務業 \| 合計 (L4)
2	消費者物價基｜總指數 (L2)	受僱員工每人｜工業及服務業 \| 合計 (L4)	消費者物價基｜總指數 (L2)	消費者物價基｜總指數 (L2)	人力資源主要｜勞動力參與率(%) \| 合計 (L6)	每人每月總薪｜工業及服務業 \| 合計 (L5)	受僱員工每人｜工業及服務業 \| 合計 (L3)	受僱員工每人｜工業及服務業 \| 合計 (L4)
3	消費者物價基｜總指數 (L3)	受僱員工進退｜進入率 \| 工業及服務業 (L4)	消費者物價基｜總指數 (L3)	消費者物價基｜總指數 (L3)	人力資源主要｜失業率(%) \| 合計 (L4)	受僱員工每人｜工業及服務業 \| 合計 (L4)	受僱員工每人｜工業及服務業 \| 合計 (L4)	受僱員工進退｜進入率 \| 工業及服務業 (L4)
4	消費者物價基｜總指數 (L4)	受僱員工進退｜進入率 \| 工業及服務業 (L6)	消費者物價基｜總指數 (L4)	消費者物價基｜總指數 (L4)	人力資源主要｜失業率(%) \| 合計 (L5)	受僱員工進退｜進入率 \| 工業及服務業 (L6)	受僱員工進退｜進入率 \| 工業及服務業 (L1)	工業生產指數｜總指數 (L3)
5	每人每月總薪｜工業及服務業 \| 合計 (L4)	工業生產指數｜總指數 (L4)	消費者物價基｜總指數 (L5)	消費者物價基｜總指數 (L5)	人力資源主要｜失業率(%) \| 合計 (L6)	工業生產指數｜總指數 (L4)	受僱員工進退｜進入率 \| 工業及服務業 (L4)	批發零售及餐｜零售業 (L5)
6	受僱員工進退｜進入率 \| 工業及服務業 (L4)	全體銀行放款｜貼現 (L2)	消費者物價基｜總指數 (L6)	消費者物價基｜總指數 (L6)	受僱員工進退｜進入率 \| 工業及服務業 (L4)	工業生產指數｜製造業 (L4)	批發零售及餐｜零售業 (L5)	利率統計-五｜金額(新臺幣百萬元) \| 購屋貸款 (L…
7	受僱員工進退｜進入率 \| 工業及服務業 (L6)	全體銀行放款｜貼現 (L3)	受僱員工進退｜進入率 \| 工業及服務業 (L4)	受僱員工進退｜進入率 \| 工業及服務業 (L4)	受僱員工進退｜進入率 \| 工業及服務業 (L5)	全體銀行放款｜貼現 (L4)	利率統計-五｜金額(新臺幣百萬元) \| 購屋貸款 (L…	全國賦稅實徵｜總計 (L3)
8	批發零售及餐｜零售業 (L5)	上市股票平均｜本益比(倍) (L3)	受僱員工進退｜進入率 \| 工業及服務業 (L6)	受僱員工進退｜進入率 \| 工業及服務業 (L6)	受僱員工進退｜進入率 \| 工業及服務業 (L6)	全國賦稅實徵｜總計 (L4)	全國賦稅實徵｜總計 (L4)	全國賦稅實徵｜總計 (L4)

tbl_lag_long

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

排序	LASSO	EN	Adaptive LASSO	Adaptive EN	Ridge	PCR	Bagging (proxy)	Random Forest
1	消費者物價基本分類指數｜總指數 (L1)	受僱員工每人每月工時｜工業及服務業 \| 合計 (L3)	消費者物價基本分類指數｜總指數 (L1)	消費者物價基本分類指數｜總指數 (L1)	人力資源主要指標｜勞動力參與率(%) \| 合計 (L1)	每人每月總薪資｜工業及服務業 \| 合計 (L4)	每人每月總薪資｜工業及服務業 \| 合計 (L4)	每人每月總薪資｜工業及服務業 \| 合計 (L4)
2	消費者物價基本分類指數｜總指數 (L2)	受僱員工每人每月工時｜工業及服務業 \| 合計 (L4)	消費者物價基本分類指數｜總指數 (L2)	消費者物價基本分類指數｜總指數 (L2)	人力資源主要指標｜勞動力參與率(%) \| 合計 (L6)	每人每月總薪資｜工業及服務業 \| 合計 (L5)	受僱員工每人每月工時｜工業及服務業 \| 合計 (L3)	受僱員工每人每月工時｜工業及服務業 \| 合計 (L4)
3	消費者物價基本分類指數｜總指數 (L3)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)	消費者物價基本分類指數｜總指數 (L3)	消費者物價基本分類指數｜總指數 (L3)	人力資源主要指標｜失業率(%) \| 合計 (L4)	受僱員工每人每月工時｜工業及服務業 \| 合計 (L4)	受僱員工每人每月工時｜工業及服務業 \| 合計 (L4)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)
4	消費者物價基本分類指數｜總指數 (L4)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L6)	消費者物價基本分類指數｜總指數 (L4)	消費者物價基本分類指數｜總指數 (L4)	人力資源主要指標｜失業率(%) \| 合計 (L5)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L6)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L1)	工業生產指數｜總指數 (L3)
5	每人每月總薪資｜工業及服務業 \| 合計 (L4)	工業生產指數｜總指數 (L4)	消費者物價基本分類指數｜總指數 (L5)	消費者物價基本分類指數｜總指數 (L5)	人力資源主要指標｜失業率(%) \| 合計 (L6)	工業生產指數｜總指數 (L4)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)	批發零售及餐飲業營業額｜零售業 (L5)
6	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)	全體銀行放款餘額統計-科目別｜貼現 (L2)	消費者物價基本分類指數｜總指數 (L6)	消費者物價基本分類指數｜總指數 (L6)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)	工業生產指數｜製造業 (L4)	批發零售及餐飲業營業額｜零售業 (L5)	利率統計-五大銀行新承做放款金額與利率-月｜金額(新臺幣百萬元) \| 購屋貸款 (L4)
7	受僱員工進退狀況｜進入率 \| 工業及服務業 (L6)	全體銀行放款餘額統計-科目別｜貼現 (L3)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L4)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L5)	全體銀行放款餘額統計-科目別｜貼現 (L4)	利率統計-五大銀行新承做放款金額與利率-月｜金額(新臺幣百萬元) \| 購屋貸款 (L4)	全國賦稅實徵淨額-按稅目別分｜總計 (L3)
8	批發零售及餐飲業營業額｜零售業 (L5)	上市股票平均每股市值及獲利概況｜本益比(倍) (L3)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L6)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L6)	受僱員工進退狀況｜進入率 \| 工業及服務業 (L6)	全國賦稅實徵淨額-按稅目別分｜總計 (L4)	全國賦稅實徵淨額-按稅目別分｜總計 (L4)	全國賦稅實徵淨額-按稅目別分｜總計 (L4)

9g)【彙總到原始變數】把多個 lag 壓回 base 變數並重新排名

# 彙總函數：把 "A_B_C_L3" → "A_B_C" 為 key，對 importance 做 sum 或 max
aggregate_importance <- function(imp_vec, method = c("sum","max")){
  method <- match.arg(method)
  if(length(imp_vec)==0) return(imp_vec)
  nm <- names(imp_vec)
  base <- sub("_L\\d+$","", nm)
  DT <- data.table(base=base, imp=as.numeric(imp_vec))
  agg <- DT[, .(imp = if (method=="sum") sum(imp, na.rm=TRUE) else max(imp, na.rm=TRUE)), by=base]
  setNames(agg$imp, agg$base)
}

imp_list_base <- lapply(imp_list_lag, aggregate_importance, method = agg_method)

# 組表（代碼名）
tbl_base_code <- data.table(排序 = rank_vec)
for (nm in names(imp_list_base)){
  tbl_base_code[[nm]] <- topk_names(imp_list_base[[nm]], k_top)
}

# 目錄：generated_excel_files/P9_重要變數/9g/當日/
out_root <- file.path("generated_excel_files", "P9_重要變數", "9g")
out_dir  <- file.path(out_root, format(Sys.Date(), "%Y-%m-%d"))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

# 輸出 xlsx（沿用你的檔名規則）
writexl::write_xlsx(
  tbl_base_code,
  path = file.path(out_dir, sprintf("variable_importance_BASE_top%d_h%d_CODE_%s.xlsx",
                                    k_top, h, toupper(agg_method)))
)

9h) 把【彙總版】映射成短/長名（不顯示 lag）

# 建 base 對照字典（不帶 lag）
all_bases <- unique(unlist(lapply(tbl_base_code[, ..model_cols], unlist), use.names=FALSE))
all_bases <- all_bases[!is.na(all_bases) & all_bases!=""]
dict_base <- if(length(all_bases)){
  # 這裡直接把 base 當 var_label 去 map（沒有 lag）
  m <- map_raw[.(all_bases), on=.(var_label)]
  data.table(var_label = all_bases,
             short_name = {
               core <- ifelse(!is.na(m$series_name), m$series_name, all_bases)
               if (short_use_dataset_prefix){
                 pref <- ifelse(!is.na(m$dataset_title), substr(m$dataset_title, 1L, dataset_prefix_nchar), "")
                 core <- ifelse(pref!="", paste0(pref, "｜", core), core)
               }
               ifelse(nchar(core) > short_max_chars,
                      paste0(substr(core, 1L, short_max_chars-1L), "…"),
                      core)
             },
             long_name  = ifelse(!is.na(m$dataset_title) & !is.na(m$series_name),
                                 paste0(m$dataset_title, "｜", m$series_name),
                                 all_bases))
} else data.table()

tbl_base_short <- copy(tbl_base_code)
tbl_base_long  <- copy(tbl_base_code)
for (mc in model_cols){
  m <- merge(data.table(var_label = tbl_base_code[[mc]]), dict_base, by="var_label", all.x=TRUE)
  tbl_base_short[[mc]] <- ifelse(is.na(m$short_name), paste0(m$var_label, " (未對應)"), m$short_name)
  tbl_base_long[[mc]]  <- ifelse(is.na(m$long_name),  paste0(m$var_label, " (未對應)"), m$long_name)
}


# 目錄：generated_excel_files/P9_重要變數/9h/當日/
out_root <- file.path("generated_excel_files", "P9_重要變數", "9h")
out_dir  <- file.path(out_root, format(Sys.Date(), "%Y-%m-%d"))
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

# 三個表分別輸出為 xlsx
writexl::write_xlsx(
  dict_base[, .(var_label, short_name, long_name)],
  path = file.path(out_dir, sprintf("var_name_dictionary_BASE_h%d.xlsx", h))
)

writexl::write_xlsx(
  tbl_base_short,
  path = file.path(out_dir, sprintf("variable_importance_BASE_top%d_h%d_SHORT_%s.xlsx",
                                    k_top, h, toupper(agg_method)))
)

writexl::write_xlsx(
  tbl_base_long,
  path = file.path(out_dir, sprintf("variable_importance_BASE_top%d_h%d_LONG_%s.xlsx",
                                    k_top, h, toupper(agg_method)))
)

writexl::write_xlsx(
  list(
    dictionary_BASE = dict_base[, .(var_label, short_name, long_name)],
    BASE_top_SHORT  = tbl_base_short,
    BASE_top_LONG   = tbl_base_long
  ),
  path = file.path(out_dir, sprintf("P9_重要變數_9h_top%d_h%d_%s.xlsx",
                                    k_top, h, toupper(agg_method)))
)

tbl_base_short

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

排序	LASSO	EN	Adaptive LASSO	Adaptive EN	Ridge	PCR	Bagging (proxy)	Random Forest
1	消費者物價基｜總指數	進口物價基本｜總指數	消費者物價基｜總指數	消費者物價基｜總指數	人力資源主要｜勞動力參與率(%) \| 合計	人力資源主要｜失業率(%) \| 合計	每人每月總薪｜工業及服務業 \| 合計	每人每月總薪｜工業及服務業 \| 合計
2	消費者物價商｜核心商品	每人每月總薪｜工業及服務業 \| 合計	消費者物價商｜核心商品	消費者物價商｜核心商品	人力資源主要｜失業率(%) \| 合計	每人每月總薪｜工業及服務業 \| 合計	受僱員工每人｜工業及服務業 \| 合計	受僱員工每人｜工業及服務業 \| 合計
3	進口物價基本｜總指數	受僱員工每人｜工業及服務業 \| 合計	進口物價基本｜總指數	進口物價基本｜總指數	受僱員工進退｜進入率 \| 工業及服務業	受僱員工每人｜工業及服務業 \| 合計	受僱員工進退｜進入率 \| 工業及服務業	受僱員工進退｜進入率 \| 工業及服務業
4	進口物價基本｜總指數	受僱員工進退｜進入率 \| 工業及服務業	進口物價基本｜總指數	進口物價基本｜總指數	利率統計-五｜利率(年息百分比率) \| 購屋貸款	受僱員工進退｜進入率 \| 工業及服務業	批發零售及餐｜零售業	工業生產指數｜總指數
5	人力資源主要｜勞動力參與率(%) \| 合計	工業生產指數｜總指數	人力資源主要｜勞動力參與率(%) \| 合計	人力資源主要｜勞動力參與率(%) \| 合計	利率統計-貨｜金融業拆款	工業生產指數｜總指數	全體銀行放款｜貼現	批發零售及餐｜批發業
6	每人每月總薪｜工業及服務業 \| 合計	批發零售及餐｜零售業	人力資源主要｜失業率(%) \| 合計	人力資源主要｜失業率(%) \| 合計	利率統計-貨｜商業本票-次級市場-31-90天	工業生產指數｜製造業	利率統計-五｜金額(新臺幣百萬元) \| 購屋貸款	批發零售及餐｜零售業
7	受僱員工進退｜進入率 \| 工業及服務業	全體銀行放款｜貼現	每人每月總薪｜工業及服務業 \| 合計	每人每月總薪｜工業及服務業 \| 合計	外匯統計-我｜新台幣 (NTD/USD)	全體銀行放款｜貼現	外匯統計-我｜日圓 (JPY/USD)	利率統計-五｜金額(新臺幣百萬元) \| 購屋貸款
8	批發零售及餐｜零售業	上市股票平均｜本益比(倍)	受僱員工進退｜進入率 \| 工業及服務業	受僱員工進退｜進入率 \| 工業及服務業	外匯統計-我｜人民幣 (CNY/USD)	全國賦稅實徵｜總計	全國賦稅實徵｜總計	全國賦稅實徵｜總計

tbl_base_long

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

排序	LASSO	EN	Adaptive LASSO	Adaptive EN	Ridge	PCR	Bagging (proxy)	Random Forest
1	消費者物價基本分類指數｜總指數	進口物價基本分類指數(按HS)(美元計價)｜總指數	消費者物價基本分類指數｜總指數	消費者物價基本分類指數｜總指數	人力資源主要指標｜勞動力參與率(%) \| 合計	人力資源主要指標｜失業率(%) \| 合計	每人每月總薪資｜工業及服務業 \| 合計	每人每月總薪資｜工業及服務業 \| 合計
2	消費者物價商品性質分類指數｜核心商品	每人每月總薪資｜工業及服務業 \| 合計	消費者物價商品性質分類指數｜核心商品	消費者物價商品性質分類指數｜核心商品	人力資源主要指標｜失業率(%) \| 合計	每人每月總薪資｜工業及服務業 \| 合計	受僱員工每人每月工時｜工業及服務業 \| 合計	受僱員工每人每月工時｜工業及服務業 \| 合計
3	進口物價基本分類指數(按HS)(美元計價)｜總指數	受僱員工每人每月工時｜工業及服務業 \| 合計	進口物價基本分類指數(按HS)(美元計價)｜總指數	進口物價基本分類指數(按HS)(美元計價)｜總指數	受僱員工進退狀況｜進入率 \| 工業及服務業	受僱員工每人每月工時｜工業及服務業 \| 合計	受僱員工進退狀況｜進入率 \| 工業及服務業	受僱員工進退狀況｜進入率 \| 工業及服務業
4	進口物價基本分類指數(按HS)(新臺幣計價)｜總指數	受僱員工進退狀況｜進入率 \| 工業及服務業	進口物價基本分類指數(按HS)(新臺幣計價)｜總指數	進口物價基本分類指數(按HS)(新臺幣計價)｜總指數	利率統計-五大銀行新承做放款金額與利率-月｜利率(年息百分比率) \| 購屋貸款	受僱員工進退狀況｜進入率 \| 工業及服務業	批發零售及餐飲業營業額｜零售業	工業生產指數｜總指數
5	人力資源主要指標｜勞動力參與率(%) \| 合計	工業生產指數｜總指數	人力資源主要指標｜勞動力參與率(%) \| 合計	人力資源主要指標｜勞動力參與率(%) \| 合計	利率統計-貨幣市場利率｜金融業拆款	工業生產指數｜總指數	全體銀行放款餘額統計-科目別｜貼現	批發零售及餐飲業營業額｜批發業
6	每人每月總薪資｜工業及服務業 \| 合計	批發零售及餐飲業營業額｜零售業	人力資源主要指標｜失業率(%) \| 合計	人力資源主要指標｜失業率(%) \| 合計	利率統計-貨幣市場利率｜商業本票-次級市場-31-90天	工業生產指數｜製造業	利率統計-五大銀行新承做放款金額與利率-月｜金額(新臺幣百萬元) \| 購屋貸款	批發零售及餐飲業營業額｜零售業
7	受僱員工進退狀況｜進入率 \| 工業及服務業	全體銀行放款餘額統計-科目別｜貼現	每人每月總薪資｜工業及服務業 \| 合計	每人每月總薪資｜工業及服務業 \| 合計	外匯統計-我國與主要貿易對手通貨對美元之匯率｜新台幣 (NTD/USD)	全體銀行放款餘額統計-科目別｜貼現	外匯統計-我國與主要貿易對手通貨對美元之匯率｜日圓 (JPY/USD)	利率統計-五大銀行新承做放款金額與利率-月｜金額(新臺幣百萬元) \| 購屋貸款
8	批發零售及餐飲業營業額｜零售業	上市股票平均每股市值及獲利概況｜本益比(倍)	受僱員工進退狀況｜進入率 \| 工業及服務業	受僱員工進退狀況｜進入率 \| 工業及服務業	外匯統計-我國與主要貿易對手通貨對美元之匯率｜人民幣 (CNY/USD)	全國賦稅實徵淨額-按稅目別分｜總計	全國賦稅實徵淨額-按稅目別分｜總計	全國賦稅實徵淨額-按稅目別分｜總計

稅收預測

柯龍

0) 安裝與載入套件

1) 參數設定

2) 讀檔與時間欄清洗（支援 2000-M1 這種格式）

3) 建立落後特徵（只用 lag，避免洩漏）

4) 指標函數與儲存容器

5) 建立 rolling 切點並執行滾動預測

6) 彙整、評分表、輸出 CSV

7) 三張樣本外走勢圖

8) Nowcast「最後一期再多推 h 期」

8a) 參數：一次預測未來 H_out 個月

8b) 準備最後一個 rolling 訓練窗 + 預測用日期索引

8c) 時間序列模型：一次 12 步

8d) Lasso 家族 + PCR：direct strategy（每個 h 一個模型）

8e) 非線性：Bagging Tree + Random Forest（direct strategy）

8f) 整併輸出（一次包含 1..H_out 步）

9) 重要變數

9a) 參數與小工具（可調整）

9b) 建最後一個 rolling 視窗資料與 direct strategy 的訓練組

9c) Lasso / Ridge / EN / Adaptive / Adaptive EN 的重要性（|β|）

9d) PCR 的重要性（loading × 迴歸係數的加權和）

9e) 樹模型重要性：Bagging（proxy）與 Random Forest（Permutation importance）

9e)【Lag 層級】組成前 k 名表（代碼名）

9f) 載入 `meta_tab.xlsx` 並把 Lag 層級表映射成短/長名

9g)【彙總到原始變數】把多個 lag 壓回 base 變數並重新排名

9h) 把【彙總版】映射成短/長名（不顯示 lag）

稅收預測

柯龍

0) 安裝與載入套件

1) 參數設定

2) 讀檔與時間欄清洗（支援 2000-M1 這種格式）

3) 建立落後特徵（只用 lag，避免洩漏）

4) 指標函數與儲存容器

5) 建立 rolling 切點並執行滾動預測

6) 彙整、評分表、輸出 CSV

7) 三張樣本外走勢圖

8) Nowcast「最後一期再多推 h 期」

8a) 參數：一次預測未來 H_out 個月

8b) 準備最後一個 rolling 訓練窗 + 預測用日期索引

8c) 時間序列模型：一次 12 步

8d) Lasso 家族 + PCR：direct strategy（每個 h 一個模型）

8e) 非線性：Bagging Tree + Random Forest（direct strategy）

8f) 整併輸出（一次包含 1..H_out 步）

9) 重要變數

9a) 參數與小工具（可調整）

9b) 建最後一個 rolling 視窗資料與 direct strategy 的訓練組

9c) Lasso / Ridge / EN / Adaptive / Adaptive EN 的重要性（|β|）

9d) PCR 的重要性（loading × 迴歸係數 的加權和）

9e) 樹模型重要性：Bagging（proxy）與 Random Forest（Permutation importance）

9e)【Lag 層級】組成前 k 名表（代碼名）

9f) 載入 meta_tab.xlsx 並把 Lag 層級表映射成短/長名

9g)【彙總到原始變數】把多個 lag 壓回 base 變數並重新排名

9h) 把【彙總版】映射成短/長名（不顯示 lag）

9d) PCR 的重要性（loading × 迴歸係數的加權和）

9f) 載入 `meta_tab.xlsx` 並把 Lag 層級表映射成短/長名