# =============================================================================
# MGT 924 - Problem Set 3 - JKP CTF competition
# Team model: "Four-signal ensemble with a Barra-style factor risk model and
#              mean-variance portfolio construction"
# Authors: Sunmeet (ridge + RFF, risk model, portfolio), Antoine (theme
#          clusters), Connor (LightGBM, integration), Anika (data)
#
# Entry point required by the competition:  main(chars, features, daily_ret)
# Returns a data.frame with columns id (integer), eom (Date), w (numeric).
#
# Pipeline (every step uses only information available at formation date t):
#   1. Features: each of the ~400 characteristics is ranked cross-sectionally
#      every month, mapped to [-0.5, 0.5], scaled to unit variance; missing
#      values are set to 0 (= cross-sectional median).
#   2. Target: next-month excess return, winsorised at the 1/99% cross-sectional
#      quantiles and cross-sectionally demeaned.
#   3. Four return signals, all refit every December on an expanding window
#      (only data whose target month has ended):
#        (a) linear ridge regression on the ranked characteristics   [Sunmeet]
#        (b) ridge regression on random Fourier features (nonlinear) [Sunmeet]
#        (c) theme portfolios: one long-short portfolio per characteristic,
#            grouped into 13 themes by hierarchical clustering, tangency
#            weights on the themes, mapped back to stocks              [Antoine]
#        (d) LightGBM gradient-boosted trees on the ranked characteristics,
#            number of trees set by early stopping on the last 5 years [Connor]
#   4. Ensemble: each signal is z-scored within the month, then averaged with
#      equal weights.
#   5. Risk model (Barra style): daily cross-sectional regressions of daily
#      excess returns on 12 industry dummies (FF12 from SIC) + style factors
#      obtained by hierarchical clustering of the characteristics. EWMA factor
#      covariance (half-life 252 days) and EWMA idiosyncratic variances
#      (half-life 63 days).
#   6. Portfolio: w = Sigma^{-1} mu (via Woodbury), scaled to a constant 10%
#      annualised ex-ante volatility.
#
# CTF Admin Modifications (2026-10-07):
# --------------------------------------
# 1. Reworded the trailing comments on lines 73, 349 and 463 to remove
#    an apostrophe.
#    Reason: at the time of the edit the pipeline syntax check rejected an
#    apostrophe in a trailing comment. It now only reports one.
#
# 2. In theme_weights, the number of themes is capped at the number of
#    usable factors.
#    Reason: the validation dataset lists 10 characteristics, fewer than the
#    13 themes requested, which stopped the run. The full dataset lists 402,
#    so the full run is unchanged.
# =============================================================================

suppressPackageStartupMessages({
  library(data.table)
  library(lightgbm)
})

CFG <- list(
  seed           = 42L,
  n_rff          = 512L,          # number of random Fourier features (sin+cos)
  rff_gamma      = 1,             # RFF bandwidth (sd of the random projection)
  lambda_mult    = 10^seq(-3, 4, by = 0.5),  # ridge grid, x mean eigenvalue
  max_val_years  = 10L,           # validation window for the ridge penalty
  winsor         = c(0.01, 0.99), # target winsorisation quantiles
  n_clusters     = 13L,           # style factors in the risk model
  cluster_months = 40L,           # months sampled to estimate feature correlations
  hl_factor      = 252,           # EWMA half-life (days), factor covariance
  hl_idio        = 63,            # EWMA half-life (days), idiosyncratic variance
  min_idio_obs   = 20L,           # min daily obs before trusting the idio var of a stock
  idio_floor_q   = 0.05,          # floor idio var at this cross-sectional quantile
  ret_clip       = 0.5,           # clip daily returns for risk estimation only
  target_vol     = 0.10,          # annualised ex-ante portfolio volatility
  portfolio      = "mv_vt",       # "mv_vt" | "mv_gross" | "rank"
  # ensemble
  signals        = c("lin", "rff", "lgb"),  # "clu" dropped: leaving it out raised the pre-1990 validation Sharpe
  # theme clusters (Antoine)
  n_themes       = 13L,           # themes of characteristic portfolios
  theme_min_months = 120L,        # months of factor returns before the first estimate
  theme_min_t    = 1,             # v2: drop factors whose past mean has |t| < 1 (Antoine, chosen on 1962-1989)
  # volatility management (Antoine, v2)
  vol_manage     = TRUE,          # scale the final portfolio by its realised volatility (adopted by rule R2: validation 5.24 -> 5.43)
  vol_window     = 6L,            # months of daily returns in the realised volatility (Antoine, chosen on 1962-1989)
  # LightGBM (Connor)
  lgb_sample     = 0.5,           # share of stock-months kept for training (hash of id, month)
  lgb_min_years  = 10L,           # years of training data before the first fit
  lgb_every      = 1L,            # refit every k-th December
  lgb_val_years  = 5L,            # early-stopping window (last years of the training data)
  lgb_max_rounds = 1000L,
  lgb_stop       = 50L,
  lgb_threads    = 8L,
  lgb_params     = list(objective = "regression", learning_rate = 0.05, num_leaves = 31L,
                        min_data_in_leaf = 500L, feature_fraction = 0.5,
                        bagging_fraction = 0.5, bagging_freq = 1L, lambda_l2 = 10,
                        max_bin = 63L),
  verbose        = TRUE
)

# ---------------------------------------------------------------- helpers ----
as_flag <- function(x) {
  if (is.logical(x)) return(!is.na(x) & x)
  if (is.numeric(x)) return(!is.na(x) & x == 1)
  x <- toupper(trimws(as.character(x)))
  x %in% c("1", "TRUE", "T")
}

# Fama-French 12 industries from 4-digit SIC
ff12 <- function(sic) {
  s <- suppressWarnings(as.integer(sic)); s[is.na(s)] <- -1L
  b <- function(lo, hi) s >= lo & s <= hi
  out <- rep(12L, length(s))
  out[b(100,999)|b(2000,2399)|b(2700,2749)|b(2770,2799)|b(3100,3199)|b(3940,3989)] <- 1L
  out[b(2500,2519)|b(2590,2599)|b(3630,3659)|b(3710,3711)|s==3714|s==3716|b(3750,3751)|
      s==3792|b(3900,3939)|b(3990,3999)] <- 2L
  out[b(2520,2589)|b(2600,2699)|b(2750,2769)|b(3000,3099)|b(3200,3569)|b(3580,3629)|
      b(3700,3709)|b(3712,3713)|s==3715|b(3717,3749)|b(3752,3791)|b(3793,3799)|
      b(3830,3839)|b(3860,3899)] <- 3L
  out[b(1200,1399)|b(2900,2999)] <- 4L
  out[b(2800,2829)|b(2840,2899)] <- 5L
  out[b(3570,3579)|b(3660,3692)|b(3694,3699)|b(3810,3829)|b(7370,7379)] <- 6L
  out[b(4800,4899)] <- 7L
  out[b(4900,4949)] <- 8L
  out[b(5000,5999)|b(7200,7299)|b(7600,7699)] <- 9L
  out[b(2830,2839)|s==3693|b(3840,3859)|b(8000,8099)] <- 10L
  out[b(6000,6999)] <- 11L
  out
}

# Cross-sectional rank transform: [-0.5, 0.5] * sqrt(12); NA -> 0 (median)
rank_std <- function(M) {
  out <- matrix(0, nrow(M), ncol(M))
  for (j in seq_len(ncol(M))) {
    x <- M[, j]; ok <- which(!is.na(x)); k <- length(ok)
    if (k > 1L) out[ok, j] <- (frank(x[ok], ties.method = "average") - 1) / (k - 1) - 0.5
  }
  out * sqrt(12)
}

demean_cols <- function(M) sweep(M, 2L, colMeans(M), check.margin = FALSE)

# z-score within the month (a constant signal carries no information -> 0)
zscore <- function(p) { s <- sd(p); if (!is.finite(s) || s == 0) return(p * 0); (p - mean(p)) / s }

# Ridge on sufficient statistics: A = sum X'X, cc = sum X'y, n = #obs
ridge_eig <- function(A, cc, n) {
  e <- eigen(A / n, symmetric = TRUE)
  list(V = e$vectors, d = pmax(e$values, 0), q = drop(crossprod(e$vectors, cc / n)))
}
ridge_beta <- function(eg, lam) drop(eg$V %*% (eg$q / (eg$d + lam)))

choose_lambda <- function(tr, va, mult) {
  eg  <- ridge_eig(tr$A, tr$c, tr$n)
  lam <- mult * mean(eg$d)
  mse <- vapply(lam, function(l) {
    b <- ridge_beta(eg, l)
    (va$yy - 2 * sum(b * va$c) + drop(crossprod(b, va$A %*% b))) / va$n
  }, numeric(1))
  list(mult = mult[which.min(mse)], mse = mse)
}

sum_stats <- function(store, yrs) {
  A <- Reduce(`+`, lapply(yrs, function(y) store[[y]]$A))
  cc <- Reduce(`+`, lapply(yrs, function(y) store[[y]]$c))
  list(A = A, c = cc,
       yy = sum(vapply(yrs, function(y) store[[y]]$yy, 0)),
       n  = sum(vapply(yrs, function(y) store[[y]]$n, 0)))
}

fit_ridge_cv <- function(store, yrs, cfg) {
  L <- length(yrs)
  nv <- min(cfg$max_val_years, floor(L / 3))
  if (nv >= 1L) {
    sel <- choose_lambda(sum_stats(store, yrs[seq_len(L - nv)]),
                         sum_stats(store, yrs[(L - nv + 1L):L]), cfg$lambda_mult)
    mult <- sel$mult
  } else mult <- 1
  full <- sum_stats(store, yrs)
  eg <- ridge_eig(full$A, full$c, full$n)
  list(beta = ridge_beta(eg, mult * mean(eg$d)), mult = mult)
}

# Style clusters from characteristic correlations (no returns used), estimated
# on months before the first test month only.
make_clusters <- function(chars, feats, idx_by_month, months, test_start, cfg) {
  pre <- months[months < test_start]
  if (length(pre) < 12L) pre <- months[seq_len(max(1L, floor(length(months) / 2)))]
  pick <- unique(pre[round(seq(1, length(pre), length.out = min(cfg$cluster_months, length(pre))))])
  S <- matrix(0, length(feats), length(feats)); n <- 0
  for (i in seq_along(pick)) {
    X <- rank_std(as.matrix(chars[idx_by_month[[as.character(pick[i])]], ..feats]))
    S <- S + crossprod(X); n <- n + nrow(X)
  }
  sdv <- sqrt(diag(S) / n)
  live <- sdv > 0.3                       # features observed for a reasonable share of stocks
  C <- (S / n)[live, live] / tcrossprod(sdv[live])
  hc <- hclust(as.dist(sqrt(pmax(2 * (1 - C), 0))), method = "ward.D2")
  cl <- rep(NA_integer_, length(feats))
  cl[live] <- cutree(hc, k = min(cfg$n_clusters, sum(live)))
  cl
}

style_exposures <- function(X, cl) {
  K <- max(cl, na.rm = TRUE)
  S <- vapply(seq_len(K), function(k) rowMeans(X[, which(cl == k), drop = FALSE]), numeric(nrow(X)))
  S <- matrix(S, nrow = nrow(X))
  sdv <- apply(S, 2, sd); sdv[!is.finite(sdv) | sdv == 0] <- 1
  sweep(demean_cols(S), 2, sdv, "/")
}

# w = Sigma^{-1} mu with Sigma = E F E' + diag(v)   (Woodbury, no F inverse)
mv_weights <- function(mu, E, F, v) {
  Dm <- mu / v
  DE <- E / v
  K  <- ncol(E)
  M  <- diag(K) + F %*% crossprod(E, DE)
  Dm - drop(DE %*% solve(M, F %*% crossprod(E, Dm)))
}
port_var <- function(w, E, F, v) { b <- crossprod(E, w); drop(crossprod(b, F %*% b)) + sum(v * w^2) }

# ------------------------------------------------- theme clusters (Antoine) --
# Characteristic portfolios: demeaned ranks scaled to $1 long and $1 short.
# Z holds ranks in [-0.5, 0.5] with NA -> 0.
factor_portfolios <- function(Z) {
  D <- demean_cols(Z)
  g <- colSums(abs(D)); g[g == 0] <- 1
  2 * sweep(D, 2L, g, "/")
}

# Orient, cluster into themes, tangency weights on the themes, back to factors.
# Fx: months x factors, past factor returns only.
theme_weights <- function(Fx, k, min_t = 0) {
  nf <- ncol(Fx)
  direction <- numeric(nf); theme <- integer(nf)
  sdv <- sqrt(colMeans(sweep(Fx, 2L, colMeans(Fx))^2))
  usable <- sdv > 0
  direction[usable] <- ifelse(colMeans(Fx[, usable, drop = FALSE]) >= 0, 1, -1)
  Fo <- sweep(Fx[, usable, drop = FALSE], 2L, direction[usable], "*")
  C <- cor(Fo)
  hc <- hclust(as.dist(sqrt(pmax(2 * (1 - C), 0))), method = "ward.D2")
  # CTF-FIX: never ask for more themes than there are usable factors (validation dataset)
  theme[usable] <- cutree(hc, k = min(k, sum(usable)))
  # v2: after clustering, drop factors whose past mean has |t| < min_t
  tst <- numeric(nf); tst[usable] <- colMeans(Fx[, usable, drop = FALSE]) / sdv[usable] * sqrt(nrow(Fx))
  direction[abs(tst) < min_t] <- 0
  used <- direction != 0
  themes <- sort(unique(theme[used]))
  Mb <- matrix(0, nf, length(themes))
  for (j in seq_along(themes)) {
    inj <- theme == themes[j] & used
    Mb[inj, j] <- direction[inj] / sum(inj)
  }
  G <- Fx %*% Mb
  tw <- solve(cov(G), colMeans(G))
  fw <- drop(Mb %*% tw)
  fw * (0.10 / sqrt(12)) / sqrt(drop(crossprod(fw, cov(Fx) %*% fw)))
}

# --------------------------------------------------------- LightGBM (Connor) --
# Deterministic subsample keyed on (id, month): the same stock-month is kept
# whatever the row order or the length of the sample.
lgb_keep <- function(id, t, share) {
  ym <- year(t) * 100L + month(t)
  ((id * 0.6180339887498949 + ym * 0.7548776662466927) %% 1) < share
}

fit_lgb <- function(store, yrs, cfg) {
  L  <- length(yrs)
  nv <- min(cfg$lgb_val_years, floor(L / 3))
  stack <- function(ys) list(X = do.call(rbind, lapply(ys, function(y) store[[y]]$X)) * 1.0,
                             y = unlist(lapply(ys, function(y) store[[y]]$y), use.names = FALSE))
  par <- c(cfg$lgb_params, list(num_threads = cfg$lgb_threads, seed = cfg$seed,
                                deterministic = TRUE, force_row_wise = TRUE, verbose = -1L))
  tr <- stack(yrs[seq_len(L - nv)]); va <- stack(yrs[(L - nv + 1L):L])
  dtr <- lgb.Dataset(tr$X, label = tr$y, params = list(max_bin = par$max_bin, verbose = -1L))
  dva <- lgb.Dataset.create.valid(dtr, va$X, label = va$y)
  m1 <- lgb.train(par, dtr, nrounds = cfg$lgb_max_rounds, valids = list(va = dva),
                  early_stopping_rounds = cfg$lgb_stop, verbose = -1L)
  n_trees <- max(m1$best_iter, 10L)
  rm(m1, dtr, dva)
  full <- stack(yrs)
  dfu <- lgb.Dataset(full$X, label = full$y, params = list(max_bin = par$max_bin, verbose = -1L))
  list(model = lgb.train(par, dfu, nrounds = n_trees, verbose = -1L), n_trees = n_trees)
}

# Ensemble portfolio variants reported in the local backtest (keep_all = TRUE)
VARIANTS <- list(team = NULL,                         # = cfg$signals
                 all4 = c("lin", "rff", "lgb", "clu"), sunmeet = c("lin", "rff"),
                 no_lin = c("rff", "clu", "lgb"), no_rff = c("lin", "clu", "lgb"),
                 no_clu = c("lin", "rff", "lgb"), no_lgb = c("lin", "rff", "clu"),
                 only_lin = "lin", only_rff = "rff", only_clu = "clu", only_lgb = "lgb")

# Volatility management (Antoine, v2; Moreira and Muir 2017): scale a portfolio
# by target / realised volatility of its base weights over the last k months.
vm_scale <- function(rv, cfg) {
  if (!length(rv)) return(1)
  (cfg$target_vol / sqrt(12)) / sqrt(mean(tail(rv, cfg$vol_window)))
}

# ---------------------------------------------------------------- engine ----
run_model <- function(chars, features, daily_ret, cfg = CFG, keep_all = FALSE) {
  set.seed(cfg$seed)
  t0 <- Sys.time()
  say <- function(...) if (isTRUE(cfg$verbose)) cat(format(Sys.time(), "%H:%M:%S"), ..., "\n")

  setDT(chars); setDT(daily_ret)
  feats <- intersect(as.character(features[[1]]), names(chars))
  chars[, eom := as.Date(eom)]
  if ("eom_ret" %in% names(chars)) chars[, eom_ret := as.Date(eom_ret)]
  else chars[, eom_ret := seq(eom + 1, by = "month", length.out = 2)[2] - 1, by = eom]
  chars[, test_flag := as_flag(ctff_test)]
  chars[, ind := ff12(sic)]
  chars[, id := as.numeric(id)]

  months <- sort(unique(chars$eom))
  idx_by_month <- split(seq_len(nrow(chars)), chars$eom)
  test_months <- sort(unique(chars$eom[chars$test_flag]))
  test_start <- if (length(test_months)) min(test_months) else max(months)
  say("rows", nrow(chars), "| features", length(feats), "| months", length(months),
      "| test months", length(test_months))

  # daily returns, keyed by year-month
  daily_ret[, date := as.Date(date)]
  daily_ret[, ym := year(date) * 100L + month(date)]
  setkey(daily_ret, ym)
  daily_ret[, id := as.numeric(id)]

  # random Fourier projection (fixed for the whole run -> deterministic)
  P <- length(feats)
  W <- matrix(rnorm(P * cfg$n_rff / 2, sd = cfg$rff_gamma / sqrt(P)), P, cfg$n_rff / 2)

  cl <- make_clusters(chars, feats, idx_by_month, months, test_start, cfg)
  say("style clusters:", paste(table(cl), collapse = " "))

  # state
  st_lin <- list(); st_rff <- list(); st_lgb <- list()
  fit <- NULL; fit_year <- NA_integer_
  lgb_fit <- NULL; lgb_year <- NA_integer_
  Fx <- matrix(0, length(months), P); theta <- NULL
  dF <- 0.5^(1 / cfg$hl_factor); dI <- 0.5^(1 / cfg$hl_idio)
  Kf <- 12L + max(cl, na.rm = TRUE)
  SF <- matrix(0, Kf, Kf); wF <- 0
  idio <- data.table(id = numeric(0), s2 = numeric(0), ws = numeric(0), nobs = integer(0))
  setkey(idio, id)
  E_prev <- NULL; id_prev <- NULL
  variants <- VARIANTS; variants$team <- cfg$signals
  Wprev <- list(); RVh <- list()                     # base weights of the previous month and realised variances

  out_w <- vector("list", length(months)); out_p <- vector("list", length(months))
  out_r <- vector("list", length(months)); lam_log <- list(); lgb_log <- list()

  for (m in seq_along(months)) {
    t <- months[m]; ix <- idx_by_month[[as.character(t)]]
    d <- chars[ix, .(id, eom_ret, ret_exc_lead1m, size_grp, ind, test_flag)]

    # ---- 1. risk model: daily returns in month t, exposures from t-1 -------
    if (!is.null(E_prev)) {
      ymt <- year(t) * 100L + month(t)
      dm <- daily_ret[J(ymt), nomatch = NULL]
      if (nrow(dm)) {
        dm[, pos := match(id, id_prev)]
        dm <- dm[!is.na(pos) & !is.na(ret_exc)]
        setorder(dm, date)
        dm[, e := NA_real_]
        for (nm in names(Wprev)) {
          Rd <- rowsum(Wprev[[nm]][dm$pos] * dm$ret_exc, dm$date, reorder = FALSE)
          RVh[[nm]] <- c(RVh[[nm]], sum(Rd^2))
        }
        days <- split(seq_len(nrow(dm)), dm$date)
        for (dd in days) {
          if (length(dd) < 2L * Kf) next
          Ed <- E_prev[dm$pos[dd], , drop = FALSE]
          r  <- pmin(pmax(dm$ret_exc[dd], -cfg$ret_clip), cfg$ret_clip)
          G  <- crossprod(Ed); G <- G + diag(1e-6 * mean(diag(G)) + 1e-10, Kf)
          f  <- solve(G, crossprod(Ed, r))
          set(dm, dd, "e", drop(r - Ed %*% f))
          SF <- dF * SF + (1 - dF) * tcrossprod(drop(f)); wF <- dF * wF + (1 - dF)
        }
        rs <- dm[!is.na(e), .(id, date, e)]
        if (nrow(rs)) {
          setorder(rs, id, date)
          rs[, `:=`(k = .N, j = seq_len(.N)), by = id]
          rs[, wt := (1 - dI) * dI^(k - j)]
          upd <- rs[, .(k = .N, a = sum(wt * e^2), b = sum(wt)), by = id]
          idio <- merge(idio, upd, by = "id", all = TRUE)
          idio[is.na(s2), `:=`(s2 = 0, ws = 0, nobs = 0L)]
          idio[is.na(k), `:=`(k = 0L, a = 0, b = 0)]
          idio[, `:=`(s2 = dI^k * s2 + a, ws = dI^k * ws + b, nobs = nobs + k)]
          idio[, c("k", "a", "b") := NULL]
          setkey(idio, id)
        }
      }
    }

    # ---- 2. features --------------------------------------------------------
    X0 <- rank_std(as.matrix(chars[ix, ..feats]))      # undemeaned ranks (LightGBM)
    X  <- demean_cols(X0)
    Z <- X %*% W; Z <- demean_cols(cbind(sin(Z), cos(Z)))
    Xf <- factor_portfolios(X0 / sqrt(12))             # characteristic portfolios
    Xq <- round(X0 * 1000); storage.mode(Xq) <- "integer"  # compact copy for LightGBM

    # ---- 3. refit every December on blocks whose targets are realised -----
    avail <- if (month(t) == 12L) year(t) else year(t) - 1L
    yrs <- sort(as.integer(names(st_lin))); yrs <- yrs[yrs <= avail]
    if (length(yrs) && (is.na(fit_year) || max(yrs) > fit_year)) {
      ky <- as.character(yrs)
      fl <- fit_ridge_cv(st_lin, ky, cfg); fr <- fit_ridge_cv(st_rff, ky, cfg)
      fit <- list(lin = fl$beta, rff = fr$beta); fit_year <- max(yrs)
      lam_log[[length(lam_log) + 1L]] <- data.table(eom = t, last_train_year = fit_year,
                                                    lam_lin = fl$mult, lam_rff = fr$mult)
      say("refit", format(t), "| train yrs", min(yrs), "-", fit_year,
          "| lambda mult lin", fl$mult, "rff", fr$mult)
    }
    ly <- sort(as.integer(names(st_lgb))); ly <- ly[ly <= avail]
    if (length(ly) >= cfg$lgb_min_years && (is.na(lgb_year) || max(ly) > lgb_year) &&
        (is.null(lgb_fit) || max(ly) %% cfg$lgb_every == 0L)) {
      t1 <- Sys.time()
      lgb_fit <- fit_lgb(st_lgb, as.character(ly), cfg); lgb_year <- max(ly)
      lgb_log[[length(lgb_log) + 1L]] <- data.table(eom = t, last_train_year = lgb_year, trees = lgb_fit$n_trees)
      say("lgb refit", format(t), "| train yrs", min(ly), "-", lgb_year, "| trees", lgb_fit$n_trees,
          "|", round(as.numeric(difftime(Sys.time(), t1, units = "secs"))), "s")
    }
    # theme weights: every December, from factor returns of earlier formations
    if (m - 1L >= cfg$theme_min_months && (is.null(theta) || month(t) == 12L))
      theta <- theme_weights(Fx[seq_len(m - 1L), , drop = FALSE], cfg$n_themes, cfg$theme_min_t)

    # ---- 4. predict + 5. portfolio -----------------------------------------
    Wnew <- list()
    if (!is.null(fit)) {
      sig <- list(lin = drop(X %*% fit$lin), rff = drop(Z %*% fit$rff))
      if (!is.null(theta)) sig$clu <- drop(Xf %*% theta)
      if (!is.null(lgb_fit)) sig$lgb <- predict(lgb_fit$model, Xq * 1.0, params = list(num_threads = cfg$lgb_threads))
      zs <- lapply(sig, zscore)
      blend <- function(s) { s <- intersect(s, names(zs)); if (!length(s)) return(NULL); Reduce(`+`, zs[s]) / length(s) }
      mu <- blend(cfg$signals)
      E <- cbind(outer(d$ind, 1:12, `==`) * 1, style_exposures(X, cl))
      w_rank <- frank(mu) - (length(mu) + 1) / 2; w_rank <- w_rank / sum(abs(w_rank))
      ws <- list(rank = w_rank)
      if (wF > 0 && nrow(idio)) {
        F  <- SF / wF * 21
        ii <- idio[J(d$id)]
        v  <- ifelse(!is.na(ii$nobs) & ii$nobs >= cfg$min_idio_obs, ii$s2 / ii$ws * 21, NA_real_)
        vmed <- ave(v, d$size_grp, FUN = function(z) if (all(is.na(z))) NA_real_ else median(z, na.rm = TRUE))
        v[is.na(v)] <- vmed[is.na(v)]; v[is.na(v)] <- median(v, na.rm = TRUE)
        if (all(is.na(v))) v[] <- 0.01
        v <- pmax(v, quantile(v, cfg$idio_floor_q, names = FALSE))
        mv_vt <- function(mu) { w <- mv_weights(mu, E, F, v); w * (cfg$target_vol / sqrt(12)) / sqrt(port_var(w, E, F, v)) }
        w_mv <- mv_weights(mu, E, F, v)
        ws$mv_vt <- mv_vt(mu)
        ws$mv_gross <- w_mv / sum(abs(w_mv))
        Wnew$team <- ws$mv_vt
        ws$mv_vm <- ws$mv_vt * vm_scale(RVh$team, cfg)
        if (keep_all) {
          for (nm in names(variants)) {
            s <- variants[[nm]]
            if (all(s %in% names(zs))) {
              vn <- paste0("v_", nm); ws[[vn]] <- mv_vt(blend(s)); Wnew[[vn]] <- ws[[vn]]
              ws[[paste0("vm_", nm)]] <- ws[[vn]] * vm_scale(RVh[[vn]], cfg)
            }
          }
          if (!is.null(sig$clu)) {                       # stand-alone portfolio by Antoine
            ws$antoine <- sig$clu; Wnew$antoine <- sig$clu
            ws$antoine_vm <- sig$clu * vm_scale(RVh$antoine, cfg)
          }
        }
      }
      wsel <- if (isTRUE(cfg$vol_manage) && !is.null(ws$mv_vm)) ws$mv_vm else ws[[cfg$portfolio]]
      if (is.null(wsel)) wsel <- ws$rank
      if (any(d$test_flag)) out_w[[m]] <- data.table(id = d$id, eom = t, w = wsel)[d$test_flag]
      if (keep_all) {
        pr <- data.table(id = d$id, eom = t, pred = mu, ret_exc_lead1m = d$ret_exc_lead1m)
        for (nm in names(sig)) set(pr, j = paste0("pred_", nm), value = sig[[nm]])
        out_p[[m]] <- pr
        out_r[[m]] <- data.table(eom = t, variant = names(ws),
                                 ret = vapply(ws, function(w) sum(w * d$ret_exc_lead1m, na.rm = TRUE), 0),
                                 gross = vapply(ws, function(w) sum(abs(w)), 0),
                                 n = length(mu))
      }
    } else if (any(d$test_flag)) {
      out_w[[m]] <- data.table(id = d$id, eom = t, w = 0)[d$test_flag]
    }

    # ---- 6. accumulate training data for block year(eom_ret) ---------------
    r0 <- d$ret_exc_lead1m; r0[is.na(r0)] <- 0
    Fx[m, ] <- drop(crossprod(Xf, r0))                 # factor returns of formation t
    ok <- !is.na(d$ret_exc_lead1m)
    if (sum(ok) > 10L) {
      y <- d$ret_exc_lead1m[ok]
      q <- quantile(y, cfg$winsor, names = FALSE); y <- pmin(pmax(y, q[1]), q[2]); y <- y - mean(y)
      yb <- as.character(year(d$eom_ret[ok][1]))
      Xo <- X[ok, , drop = FALSE]; Zo <- Z[ok, , drop = FALSE]
      add <- function(s, M) {
        if (is.null(s)) s <- list(A = 0, c = 0, yy = 0, n = 0)
        list(A = s$A + crossprod(M), c = s$c + drop(crossprod(M, y)), yy = s$yy + sum(y^2), n = s$n + length(y))
      }
      st_lin[[yb]] <- add(st_lin[[yb]], Xo); st_rff[[yb]] <- add(st_rff[[yb]], Zo)
      kp <- lgb_keep(d$id[ok], t, cfg$lgb_sample)
      sl <- st_lgb[[yb]]
      st_lgb[[yb]] <- list(X = rbind(sl$X, Xq[ok, , drop = FALSE][kp, , drop = FALSE]), y = c(sl$y, y[kp]))
    }

    # ---- exposures at t for next month's risk regressions ------------------
    E_prev <- cbind(outer(d$ind, 1:12, `==`) * 1, style_exposures(X, cl)); id_prev <- d$id
    Wprev <- Wnew
    if (m %% 60L == 0L) say("month", format(t), "| elapsed", round(as.numeric(difftime(Sys.time(), t0, units = "mins")), 1), "min")
  }

  weights <- rbindlist(out_w)
  if (!nrow(weights)) weights <- data.table(id = numeric(0), eom = as.Date(character(0)), w = numeric(0))
  weights[!is.finite(w), w := 0]
  res <- list(weights = data.frame(id = as.integer(weights$id), eom = as.Date(weights$eom), w = as.numeric(weights$w)),
              lambdas = rbindlist(lam_log), lgb_trees = rbindlist(lgb_log), clusters = data.table(feature = feats, cluster = cl))
  if (keep_all) { res$preds <- rbindlist(out_p, fill = TRUE); res$returns <- rbindlist(out_r) }
  say("done in", round(as.numeric(difftime(Sys.time(), t0, units = "mins")), 1), "min")
  res
}

main <- function(chars, features, daily_ret) {
  set.seed(CFG$seed)
  run_model(chars, features, daily_ret, cfg = CFG, keep_all = FALSE)$weights
}