# =============================================================================
# JKP Common Task Framework submission
# Ensemble return forecasts (ridge + LightGBM) mapped into a beta- and
# dollar-neutral mean-variance portfolio using a Barra-style factor risk model
# (industry + algorithmically clustered characteristic factors) estimated from
# daily returns. Volatility-targeted to a constant ex-ante 10% per year.
#
# Look-ahead discipline: every quantity used for the portfolio at month-end t
# is computed from rows with eom <= t, labels realised by t (eom_ret <= t) and
# daily returns dated <= t. Refit dates and the risk-model burn-in are anchored
# to the first test month, so truncating the sample at the end never changes
# earlier weights.
#
# CTF Admin Modifications (2026-10-08):
# --------------------------------------
# 1. Added renv.lock declaring lightgbm 4.7.0 (with R6 and jsonlite).
#    Reason: the script loads lightgbm, which is not pre-installed, and no
#    lockfile was uploaded. The version is the current release, not
#    necessarily the one the authors used.
# =============================================================================

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

CONFIG <- list(
  seed = 20261006L,
  threads = 16L,
  refit_every_years = 1L,      # expanding-window refit every January portfolio
  holdout_months = 60L,        # time-ordered validation block inside training data
  min_train_months = 36L,
  ridge_grid = 10^seq(-3, 2, by = 0.5),
  lgb_max_rounds = 1500L,
  lgb_early_stop = 50L,
  lgb_configs = list(
    list(name = "lgb_small", num_leaves = 31L, min_data_in_leaf = 2000L,
         feature_fraction = 0.5, bagging_fraction = 0.5, lambda_l2 = 10, learning_rate = 0.05),
    list(name = "lgb_large", num_leaves = 127L, min_data_in_leaf = 5000L,
         feature_fraction = 0.3, bagging_fraction = 0.5, lambda_l2 = 50, learning_rate = 0.05)
  ),
  n_clusters = 12L,            # characteristic clusters used as risk factors
  cluster_method = "ward.D2",  # hclust linkage (average linkage chains into one giant cluster)
  cluster_months = 120L,       # months of pre-test data used to form clusters
  risk_burnin_years = 4L,      # daily history used before the first test month
  hl_var = 84,                 # EWMA half-lives in trading days
  hl_cor = 504,
  hl_idio = 84,
  min_idio_obs = 20L,
  vol_target = 0.10,           # annualised ex-ante volatility
  method = "mv_bn",            # "rank", "mv", or "mv_bn" (beta- and dollar-neutral MV)
  ridge_weight = 0,            # weight of the linear model in the ensemble (set on pre-1990 validation)
  lgb_fixed_rounds = 400L,     # fixed number of trees (chosen on pre-1990 validation over early stopping and 1000); NULL = early stopping
  idio_shrink = 0,             # shrink idiosyncratic variances toward the cross-sectional median
  fcov_mult = 1,               # inflation of the factor covariance in the optimiser
  target_arctan = 0,           # 1: target = arctan of within-month z-scored return; 0: within-month rank
  industry_neutral = 1,        # 1: demean the forecast within SIC group before portfolio construction
  book_mix = 0.5,              # weight on Book B (Ledoit-Wolf + Tikhonov covariance); 0 = factor book only
  lw_window = 756L,            # trading days in the Book B covariance window
  lw_min_obs = 250L            # minimum observed days for a stock to enter Book B
)

log_msg <- function(...) {
  cat(format(Sys.time(), "%H:%M:%S"), ..., "\n")
  flush.console()
}

as_flag <- function(x) {
  if (is.logical(x)) return(!is.na(x) & x)
  as.character(x) %in% c("1", "TRUE", "true", "T", "t")
}

month_end <- function(d) {
  d <- as.Date(d)
  first <- as.Date(format(d, "%Y-%m-01"))
  as.Date(sapply(first, function(f) seq(f, by = "month", length.out = 2)[2] - 1), origin = "1970-01-01")
}

# Cross-sectional rank to [-0.5, 0.5]; missing values map to the median (0)
rank_center <- function(x) {
  n <- sum(!is.na(x))
  if (n < 2) return(rep(0, length(x)))
  r <- (frank(x, na.last = "keep", ties.method = "average") - 1) / (n - 1) - 0.5
  r[is.na(r)] <- 0
  r
}

zscore <- function(x) {
  s <- sd(x)
  if (!is.finite(s) || s == 0) return(rep(0, length(x)))
  (x - mean(x)) / s
}

zscore_na <- function(x) {
  ok <- !is.na(x)
  out <- rep(NA_real_, length(x))
  if (sum(ok) > 1) out[ok] <- zscore(x[ok])
  out
}

# -----------------------------------------------------------------------------
# 1. Data preparation
# -----------------------------------------------------------------------------
prepare_data <- function(chars, features, cfg = CONFIG) {
  feats <- sort(as.character(features$features))
  dt <- as.data.table(chars)
  feats <- intersect(feats, names(dt))
  keep <- unique(c("id", "eom", "eom_ret", "sic", "size_grp", "ret_exc_lead1m", "ctff_test",
                   "market_equity", feats))
  dt <- dt[, intersect(keep, names(dt)), with = FALSE]
  dt[, id := as.numeric(id)]
  dt[, eom := as.Date(eom)]
  if ("eom_ret" %in% names(dt)) dt[, eom_ret := as.Date(eom_ret)] else dt[, eom_ret := month_end(eom + 1)]
  dt[, test := as_flag(ctff_test)]
  dt[, me_raw := as.numeric(market_equity)]
  setorder(dt, eom, id)
  dt[, (feats) := lapply(.SD, rank_center), by = eom, .SDcols = feats]
  if (isTRUE(cfg$target_arctan == 1)) {
    dt[, y := atan(zscore_na(ret_exc_lead1m)), by = eom]
  } else {
    dt[, y := rank_center(ret_exc_lead1m), by = eom]
  }
  dt[is.na(ret_exc_lead1m), y := NA_real_]
  list(dt = dt, feats = feats)
}

# -----------------------------------------------------------------------------
# 2. Return models
# -----------------------------------------------------------------------------
# Per-month cross-products for ridge: training sums are then cheap to assemble
month_crossprods <- function(X, y, month_rows) {
  lapply(month_rows, function(idx) {
    ok <- idx[!is.na(y[idx])]
    Xm <- X[ok, , drop = FALSE]
    list(xtx = crossprod(Xm), xty = crossprod(Xm, y[ok]), yy = sum(y[ok]^2), n = length(ok))
  })
}

sum_cp <- function(cps) {
  list(xtx = Reduce(`+`, lapply(cps, `[[`, "xtx")), xty = Reduce(`+`, lapply(cps, `[[`, "xty")),
       yy = sum(sapply(cps, `[[`, "yy")), n = sum(sapply(cps, `[[`, "n")))
}

fit_ridge <- function(cp_fit, cp_hold, cp_all, grid) {
  p <- nrow(cp_all$xtx)
  ridge <- function(cp, g) solve(cp$xtx + diag(g * cp$n / 12, p), cp$xty)
  g_best <- grid[ceiling(length(grid) / 2)]
  if (!is.null(cp_hold) && cp_hold$n > 0) {
    mse <- sapply(grid, function(g) {
      b <- ridge(cp_fit, g)
      cp_hold$yy - 2 * sum(b * cp_hold$xty) + drop(t(b) %*% cp_hold$xtx %*% b)
    })
    g_best <- grid[which.min(mse)]
  }
  list(beta = ridge(cp_all, g_best), g = g_best)
}

fit_lgb <- function(X, y, idx_fit, idx_hold, idx_all, cfg, lc, seed) {
  params <- list(objective = "regression", learning_rate = lc$learning_rate, num_leaves = lc$num_leaves,
                 min_data_in_leaf = lc$min_data_in_leaf, feature_fraction = lc$feature_fraction,
                 bagging_fraction = lc$bagging_fraction, bagging_freq = 1L, lambda_l2 = lc$lambda_l2,
                 max_bin = 63L, num_threads = cfg$threads, deterministic = TRUE, force_row_wise = TRUE,
                 seed = seed, verbose = -1L)
  n_rounds <- if (is.null(cfg$lgb_fixed_rounds)) 300L else as.integer(cfg$lgb_fixed_rounds)
  if (is.null(cfg$lgb_fixed_rounds) && length(idx_hold) > 0) {
    d_fit <- lgb.Dataset(X[idx_fit, , drop = FALSE], label = y[idx_fit], params = list(max_bin = 63L))
    d_hold <- lgb.Dataset.create.valid(d_fit, X[idx_hold, , drop = FALSE], label = y[idx_hold])
    m <- lgb.train(params, d_fit, nrounds = cfg$lgb_max_rounds, valids = list(hold = d_hold),
                   early_stopping_rounds = cfg$lgb_early_stop, verbose = -1L)
    n_rounds <- max(50L, m$best_iter)
    rm(m, d_fit, d_hold)
  }
  d_all <- lgb.Dataset(X[idx_all, , drop = FALSE], label = y[idx_all], params = list(max_bin = 63L))
  m <- lgb.train(params, d_all, nrounds = n_rounds, verbose = -1L)
  list(model = m, rounds = n_rounds)
}

predict_returns <- function(dt, feats, test_start, cfg) {
  X <- as.matrix(dt[, ..feats])
  y <- dt$y
  months <- sort(unique(dt$eom))
  month_rows <- split(seq_len(nrow(dt)), factor(dt$eom, levels = months))
  log_msg("ridge cross-products for", length(months), "months")
  cps <- month_crossprods(X, y, month_rows)
  names(cps) <- as.character(months)
  month_ret <- dt[, .(eom_ret = max(eom_ret)), by = eom]
  eom_ret_of <- setNames(month_ret$eom_ret, as.character(month_ret$eom))

  pred_months <- months[months >= test_start]
  start_year <- as.integer(format(test_start, "%Y"))
  is_refit <- function(m) {
    format(m, "%m") == format(test_start, "%m") &&
      (as.integer(format(m, "%Y")) - start_year) %% cfg$refit_every_years == 0
  }
  refits <- pred_months[vapply(pred_months, is_refit, logical(1))]
  if (length(refits) == 0 || refits[1] != pred_months[1]) refits <- c(pred_months[1], refits)
  refits <- sort(unique(refits))

  model_names <- c("ridge", vapply(cfg$lgb_configs, `[[`, "", "name"))
  out <- vector("list", length(refits))
  diag_rows <- vector("list", length(refits))
  for (k in seq_along(refits)) {
    r <- refits[k]
    next_r <- if (k < length(refits)) refits[k + 1] else as.Date("9999-12-31")
    use_months <- pred_months[pred_months >= r & pred_months < next_r]
    train_months <- months[eom_ret_of[as.character(months)] <= r]
    if (length(train_months) < cfg$min_train_months) {
      log_msg("refit", format(r), ": only", length(train_months), "training months, skipped")
      next
    }
    n_hold <- min(cfg$holdout_months, floor(length(train_months) / 3))
    hold_months <- tail(train_months, n_hold)
    fit_months <- head(train_months, length(train_months) - n_hold)
    rows_of <- function(ms) {
      idx <- unlist(month_rows[as.character(ms)], use.names = FALSE)
      idx[!is.na(y[idx])]
    }
    idx_fit <- rows_of(fit_months); idx_hold <- rows_of(hold_months); idx_all <- c(idx_fit, idx_hold)
    idx_pred <- unlist(month_rows[as.character(use_months)], use.names = FALSE)

    rg <- fit_ridge(sum_cp(cps[as.character(fit_months)]), sum_cp(cps[as.character(hold_months)]),
                    sum_cp(cps[as.character(train_months)]), cfg$ridge_grid)
    preds <- data.table(row = idx_pred, ridge = drop(X[idx_pred, , drop = FALSE] %*% rg$beta))
    rounds <- integer(0)
    for (j in seq_along(cfg$lgb_configs)) {
      lc <- cfg$lgb_configs[[j]]
      fit <- fit_lgb(X, y, idx_fit, idx_hold, idx_all, cfg, lc, seed = cfg$seed + j)
      set(preds, j = lc$name, value = predict(fit$model, X[idx_pred, , drop = FALSE]))
      rounds <- c(rounds, fit$rounds)
      rm(fit)
    }
    out[[k]] <- preds
    diag_rows[[k]] <- data.table(refit = r, n_train = length(idx_all), ridge_g = rg$g,
                                 t(setNames(rounds, vapply(cfg$lgb_configs, `[[`, "", "name"))))
    log_msg("refit", format(r), ": train rows", length(idx_all), "ridge g", rg$g,
            "lgb rounds", paste(rounds, collapse = "/"))
    invisible(gc())
  }
  preds <- rbindlist(out)
  if (nrow(preds) == 0) return(list(preds = NULL, diag = NULL, models = model_names))
  preds[, `:=`(id = dt$id[row], eom = dt$eom[row])]
  preds[, row := NULL]
  # Ensemble of the linear and tree families, z-scored within month
  lgb_names <- vapply(cfg$lgb_configs, `[[`, "", "name")
  preds[, (model_names) := lapply(.SD, zscore), by = eom, .SDcols = model_names]
  preds[, lgb := rowMeans(.SD), .SDcols = lgb_names]
  preds[, ens := zscore(cfg$ridge_weight * ridge + (1 - cfg$ridge_weight) * lgb), by = eom]
  list(preds = preds, diag = rbindlist(diag_rows, fill = TRUE), models = c(model_names, "lgb", "ens"))
}

# -----------------------------------------------------------------------------
# 3. Risk model
# -----------------------------------------------------------------------------
sic_group <- function(sic) {
  s <- suppressWarnings(as.integer(substr(sprintf("%04d", as.integer(sic)), 1, 2)))
  g <- cut(s, breaks = c(-Inf, 9, 14, 17, 29, 39, 48, 49, 51, 59, 67, 79, 89, Inf),
           labels = c("agri", "mining", "constr", "manuf_nd", "manuf_d", "transp", "util",
                      "wholesale", "retail", "finance", "services_a", "services_b", "public"))
  g <- as.character(g)
  g[is.na(g)] <- "other"
  g
}

# Clusters formed once from pre-test data only: correlation of ranked
# characteristics, Ward hclust on 1 - |rho|, oriented by PC1 sign.
form_clusters <- function(dt, feats, test_start, cfg) {
  months <- sort(unique(dt$eom[dt$eom < test_start]))
  months <- tail(months, cfg$cluster_months)
  X <- as.matrix(dt[eom %in% months, ..feats])
  C <- suppressWarnings(cor(X))
  C[!is.finite(C)] <- 0
  diag(C) <- 1
  hc <- hclust(as.dist(1 - abs(C)), method = cfg$cluster_method)
  cl <- cutree(hc, k = min(cfg$n_clusters, length(feats)))
  sgn <- numeric(length(feats))
  for (k in sort(unique(cl))) {
    j <- which(cl == k)
    if (length(j) == 1) { sgn[j] <- 1; next }
    v <- eigen(C[j, j], symmetric = TRUE)$vectors[, 1]
    if (sum(v) < 0) v <- -v
    sgn[j] <- ifelse(v >= 0, 1, -1)
  }
  list(cluster = setNames(cl, feats), sign = setNames(sgn, feats), hc = hc)
}

build_exposures <- function(dt, feats, clusters) {
  X <- as.matrix(dt[, ..feats])
  K <- sort(unique(clusters$cluster))
  E <- sapply(K, function(k) {
    j <- which(clusters$cluster == k)
    drop(X[, j, drop = FALSE] %*% clusters$sign[j]) / length(j)
  })
  colnames(E) <- paste0("clu", K)
  ex <- data.table(id = dt$id, eom = dt$eom, me = dt$me_raw, ind = sic_group(dt$sic), size_grp = dt$size_grp)
  ex <- cbind(ex, E)
  cl_cols <- colnames(E)
  ex[, (cl_cols) := lapply(.SD, zscore), by = eom, .SDcols = cl_cols]
  ind_levels <- sort(unique(ex$ind))
  for (g in ind_levels) set(ex, j = paste0("ind_", g), value = as.numeric(ex$ind == g))
  list(ex = ex, cols = c(paste0("ind_", ind_levels), cl_cols))
}

ewma_bias_corrected <- function(x, hl) {
  lam <- 0.5^(1 / hl)
  num <- stats::filter((1 - lam) * x, lam, method = "recursive")
  den <- 1 - lam^seq_along(x)
  as.numeric(num) / den
}

estimate_risk <- function(expo, daily_ret, risk_months, test_start, cfg) {
  ex <- expo$ex
  cols <- expo$cols
  start <- as.Date(paste0(as.integer(format(test_start, "%Y")) - cfg$risk_burnin_years,
                          format(test_start, "-%m-01")))
  dr <- as.data.table(daily_ret)
  setnames(dr, names(dr)[1:3], c("id", "date", "ret"))
  dr <- dr[as.Date(date) >= start & as.Date(date) <= max(risk_months),
           .(id = as.numeric(id), date = as.Date(date), ret = pmin(pmax(ret, -0.5), 0.5))]
  dr[, eom_prev := as.Date(format(date, "%Y-%m-01")) - 1]   # exposures from previous month-end
  setorder(dr, date, id)
  ex <- ex[eom >= min(dr$eom_prev) & is.finite(me) & me > 0]
  setkey(ex, eom, id)
  p <- length(cols)
  days <- sort(unique(dr$date))
  log_msg("risk model: daily cross-sectional regressions,", length(days), "days")

  # Month by month: exposures known at the previous month-end, one WLS per day
  fret <- matrix(0, length(days), p, dimnames = list(NULL, cols))
  res <- rep(NA_real_, nrow(dr))
  day_index <- match(dr$date, days)
  month_rows <- split(seq_len(nrow(dr)), dr$eom_prev)
  for (mp in names(month_rows)) {
    i_m <- month_rows[[mp]]
    E <- ex[J(as.Date(mp))]
    if (nrow(E) == 0 || is.na(E$id[1])) next
    pos <- match(dr$id[i_m], E$id)
    i_m <- i_m[!is.na(pos)]; pos <- pos[!is.na(pos)]
    Xm <- as.matrix(E[, ..cols]); wm <- sqrt(E$me)
    for (dd in split(seq_along(i_m), day_index[i_m])) {
      i <- i_m[dd]; xi <- Xm[pos[dd], , drop = FALSE]
      fit <- lm.wfit(xi, dr$ret[i], wm[pos[dd]])
      b <- fit$coefficients
      b[is.na(b)] <- 0
      fret[day_index[i[1]], ] <- b
      res[i] <- dr$ret[i] - drop(xi %*% b)
    }
  }
  dr <- dr[!is.na(res), .(id, date, res2 = res[!is.na(res)]^2)]
  rm(res); invisible(gc())

  # Idiosyncratic variance: per-stock bias-corrected EWMA of squared residuals,
  # sampled on each stock's last trading day of every month
  setorder(dr, id, date)
  dr[, `:=`(ivar = ewma_bias_corrected(res2, cfg$hl_idio), nobs = seq_len(.N)), by = id]
  dr[, ym := format(date, "%Y-%m")]
  n <- nrow(dr)
  last <- c(dr$id[-1] != dr$id[-n] | dr$ym[-1] != dr$ym[-n], TRUE)
  idio_m <- dr[last & nobs >= cfg$min_idio_obs, .(id, ym, ivar = 21 * ivar)]
  rm(dr); invisible(gc())

  # Factor covariance: EWMA variances (short half-life) x EWMA correlations (long)
  lam_v <- 0.5^(1 / cfg$hl_var); lam_c <- 0.5^(1 / cfg$hl_cor)
  Sv <- matrix(0, p, p); Sc <- matrix(0, p, p); wv <- 0; wc <- 0
  ym_days <- format(days, "%Y-%m")
  is_snap <- c(ym_days[-1] != ym_days[-length(days)], TRUE)
  fcov <- list()
  for (d in seq_along(days)) {
    ff <- tcrossprod(fret[d, ])
    Sv <- lam_v * Sv + (1 - lam_v) * ff; wv <- lam_v * wv + (1 - lam_v)
    Sc <- lam_c * Sc + (1 - lam_c) * ff; wc <- lam_c * wc + (1 - lam_c)
    if (is_snap[d]) {
      V <- Sv / wv; Cc <- Sc / wc
      sd_v <- sqrt(pmax(diag(V), 1e-12)); sd_c <- sqrt(pmax(diag(Cc), 1e-12))
      R <- Cc / tcrossprod(sd_c)
      fcov[[ym_days[d]]] <- 21 * (R * tcrossprod(sd_v))
    }
  }
  list(fcov = fcov, idio = idio_m, fret = data.table(date = days, fret), cols = cols)
}

# -----------------------------------------------------------------------------
# 4. Portfolio construction
# -----------------------------------------------------------------------------
# Sigma^{-1} V via Woodbury: Sigma = X F X' + D
sigma_inv <- function(V, X, F, dvec) {
  DiV <- V / dvec
  DiX <- X / dvec
  M <- diag(ncol(X)) + F %*% crossprod(X, DiX)
  DiV - DiX %*% solve(M, F %*% crossprod(X, DiV))
}

portfolio_month <- function(sub, ym, risk, cols, method, vol_target, idio_shrink = 0, fcov_mult = 1) {
  n <- nrow(sub)
  F <- risk$fcov[[ym]]
  if (!is.null(F)) F <- fcov_mult * F
  if (method == "rank" || is.null(F)) {
    r <- frank(sub$signal, ties.method = "average")
    w <- r - mean(r)
    return(w / sum(abs(w)) * 2)
  }
  X <- as.matrix(sub[, ..cols])
  dvec <- sub$ivar
  med <- median(dvec, na.rm = TRUE)
  dvec[!is.finite(dvec)] <- med
  dvec <- pmax(dvec, quantile(dvec, 0.05, na.rm = TRUE))   # floor to limit leverage on tiny-variance names
  dvec <- (1 - idio_shrink) * dvec + idio_shrink * median(dvec)
  mu <- sub$signal
  w <- drop(sigma_inv(matrix(mu), X, F, dvec))
  if (method == "mv_bn") {
    me <- sub$me; me[!is.finite(me)] <- 0
    wm <- me / sum(me)
    Sw <- drop(X %*% (F %*% crossprod(X, wm))) + dvec * wm
    beta <- Sw / sum(wm * Sw)
    A <- cbind(1, beta)
    SiA <- sigma_inv(A, X, F, dvec)
    lam <- solve(crossprod(A, SiA), crossprod(A, w))
    w <- w - drop(SiA %*% lam)
  }
  vol <- sqrt(12 * (sum(drop(crossprod(X, w)) * drop(F %*% crossprod(X, w))) + sum(dvec * w^2)))
  if (!is.finite(vol) || vol <= 0) return(rep(0, n))
  w * vol_target / vol
}

build_portfolios <- function(preds, signal_col, expo, risk, method, vol_target, idio_shrink = 0, fcov_mult = 1) {
  sig <- preds[, .(id, eom, signal = get(signal_col))]
  ex <- expo$ex[, c("id", "eom", "me", expo$cols), with = FALSE]
  sig <- merge(sig, ex, by = c("id", "eom"), all.x = TRUE)
  for (cname in expo$cols) set(sig, which(is.na(sig[[cname]])), cname, 0)
  sig[, ym := format(eom, "%Y-%m")]
  sig <- merge(sig, risk$idio, by = c("id", "ym"), all.x = TRUE)
  setorder(sig, eom, id)
  sig[, w := portfolio_month(.SD, ym[1], risk, expo$cols, method, vol_target, idio_shrink, fcov_mult), by = eom,
      .SDcols = c("signal", "me", "ivar", expo$cols)]
  sig[, .(id, eom, w)]
}

# Industry-neutral signal: demean the forecast within SIC group each month
neutralize_signal <- function(preds, signal_col, expo) {
  p <- merge(preds[, c("id", "eom", signal_col), with = FALSE], expo$ex[, .(id, eom, ind)],
             by = c("id", "eom"), all.x = TRUE)
  p[is.na(ind), ind := "other"]
  p[, (signal_col) := get(signal_col) - mean(get(signal_col)), by = .(eom, ind)]
  p[, (signal_col) := zscore(get(signal_col)), by = eom]
  p[, ind := NULL]
  p
}

# -----------------------------------------------------------------------------
# 4b. Book B: Ledoit-Wolf shrunk daily covariance with a Tikhonov ridge
# -----------------------------------------------------------------------------
# For month-end t: X = trailing lw_window trading days (<= t) of daily excess
# returns for that month's stocks with >= lw_min_obs observations (missing days
# zero, columns demeaned). S = X'X/T is shrunk toward nu*I (Ledoit-Wolf 2004),
# and A = Sigma_LW + nu*I is the ridge matrix used for the mean-variance solve.
# All solves use Woodbury on the T x T matrix XX', never an N x N matrix.
lw_month <- function(X, mu, me, vol_target) {
  T <- nrow(X); N <- ncol(X)
  G <- tcrossprod(X)
  g <- diag(G)
  nu <- sum(g) / (T * N)
  s2 <- sum(G^2) / T^2                       # ||S||_F^2
  d2 <- s2 - nu^2 * N                        # ||S - nu I||_F^2
  b2 <- min((sum(g^2) - T * s2) / T^2, d2)
  delta <- if (d2 > 0) min(max(b2 / d2, 0), 1) else 1
  k <- (1 - delta) / T
  cc <- nu * (1 + delta)                     # shrinkage target nu plus ridge lambda = nu
  if (k > 0) R <- chol(G + diag(cc / k, T))
  a_inv <- function(V) {
    V <- as.matrix(V)
    if (k <= 0) return(V / cc)
    (V - crossprod(X, backsolve(R, backsolve(R, X %*% V, transpose = TRUE)))) / cc
  }
  sig_lw <- function(v) (1 - delta) * drop(crossprod(X, X %*% v)) / T + delta * nu * v
  wm <- pmax(me, 0); wm[!is.finite(wm)] <- 0
  beta <- if (sum(wm) > 0) { wm <- wm / sum(wm); sw <- sig_lw(wm); sw / sum(wm * sw) } else rep(1, N)
  C <- cbind(1, beta)
  a <- a_inv(mu)
  AC <- a_inv(C)
  w <- drop(a - AC %*% solve(crossprod(C, AC), crossprod(C, a)))
  vol <- sqrt(252 * sum(w * sig_lw(w)))
  if (!is.finite(vol) || vol <= 0) return(rep(0, N))
  w * vol_target / vol
}

lw_book <- function(sig, daily_ret, cfg) {
  dr <- as.data.table(daily_ret)
  setnames(dr, names(dr)[1:3], c("id", "date", "ret"))
  months <- sort(unique(sig$eom))
  start <- min(months) - ceiling(cfg$lw_window * 1.6)
  dr <- dr[as.Date(date) >= start & as.Date(date) <= max(months),
           .(id = as.numeric(id), date = as.Date(date), ret = pmin(pmax(ret, -0.5), 0.5))]
  setorder(dr, date, id)
  days <- sort(unique(dr$date))
  dr[, d := match(date, days)]
  rng <- dr[, .(r1 = .I[1], r2 = .I[.N]), by = d]
  r1 <- rng$r1; r2 <- rng$r2
  log_msg("Book B: Ledoit-Wolf + Tikhonov covariance for", length(months), "months")
  out <- vector("list", length(months))
  for (m in seq_along(months)) {
    t <- months[m]
    d_end <- findInterval(as.numeric(t), as.numeric(days))
    if (d_end < cfg$lw_min_obs) next
    d_start <- max(1L, d_end - cfg$lw_window + 1L)
    s <- sig[eom == t]
    rows <- r1[d_start]:r2[d_end]
    sub <- dr[rows][id %in% s$id]
    col <- match(sub$id, s$id)
    elig <- which(tabulate(col, nrow(s)) >= cfg$lw_min_obs)
    if (length(elig) < 50) next
    keep <- col %in% elig
    X <- matrix(0, d_end - d_start + 1L, length(elig))
    X[cbind(sub$d[keep] - d_start + 1L, match(col[keep], elig))] <- sub$ret[keep]
    X <- sweep(X, 2, colMeans(X))
    w <- lw_month(X, s$signal[elig], s$me[elig], cfg$vol_target)
    out[[m]] <- data.table(id = s$id[elig], eom = t, w = w)
  }
  rbindlist(out)
}

# Blend of the factor book (A) and the Ledoit-Wolf book (B), each at the same ex-ante volatility
blend_books <- function(wA, wB, mix) {
  b <- merge(wA, wB, by = c("id", "eom"), all = TRUE, suffixes = c("_a", "_b"))
  b[is.na(w_a), w_a := 0]; b[is.na(w_b), w_b := 0]
  b[, .(id, eom, w = (1 - mix) * w_a + mix * w_b)]
}

# -----------------------------------------------------------------------------
# 5. Driver
# -----------------------------------------------------------------------------
run_model <- function(chars, features, daily_ret, cfg = CONFIG, test_start = NULL, keep = FALSE) {
  set.seed(cfg$seed)
  t0 <- Sys.time()
  prep <- prepare_data(chars, features, cfg)
  dt <- prep$dt; feats <- prep$feats
  if (is.null(test_start)) test_start <- min(dt$eom[dt$test])
  test_start <- as.Date(test_start)
  log_msg("prepared", nrow(dt), "rows,", length(feats), "features; first prediction month", format(test_start))

  pr <- predict_returns(dt, feats, test_start, cfg)
  if (is.null(pr$preds)) {
    log_msg("no training data available: falling back to equal weights")
    out <- dt[eom >= test_start, .(id = as.integer(id), eom, w = 1 / .N), by = eom]
    return(list(weights = out[, .(id, eom, w)]))
  }
  clusters <- form_clusters(dt, feats, test_start, cfg)
  expo <- build_exposures(dt, feats, clusters)
  rm(prep); dt[, (feats) := NULL]; invisible(gc())
  risk <- estimate_risk(expo, daily_ret, sort(unique(pr$preds$eom)), test_start, cfg)
  log_msg("building portfolios:", cfg$method)
  sig_preds <- if (isTRUE(cfg$industry_neutral == 1)) neutralize_signal(pr$preds, "ens", expo) else pr$preds
  w <- build_portfolios(sig_preds, "ens", expo, risk, cfg$method, cfg$vol_target, cfg$idio_shrink, cfg$fcov_mult)
  if (cfg$book_mix > 0) {
    sigB <- merge(sig_preds[, .(id, eom, signal = ens)], expo$ex[, .(id, eom, me)], by = c("id", "eom"), all.x = TRUE)
    wB <- lw_book(sigB, daily_ret, cfg)
    w <- blend_books(w, wB, cfg$book_mix)
  }
  w[!is.finite(w), w := 0]
  w[, id := as.integer(id)]
  log_msg("done in", format(round(difftime(Sys.time(), t0, units = "mins"), 1)))
  res <- list(weights = as.data.frame(w[, .(id, eom, w)]))
  if (keep) res <- c(res, list(preds = pr$preds, diag = pr$diag, clusters = clusters,
                               expo = expo, risk = risk, outcomes = dt[, .(id, eom, ret_exc_lead1m, me_raw, size_grp)]))
  res
}

main <- function(chars, features, daily_ret) {
  set.seed(CONFIG$seed)
  run_model(chars, features, daily_ret, CONFIG)$weights
}

# Local testing only (not executed by the CTF pipeline)
if (interactive()) {
  library(arrow)
  chars <- read_parquet("ctff_chars.parquet")
  features <- read_parquet("ctff_features.parquet")
  daily_ret <- read_parquet("ctff_daily_ret.parquet")
  pf <- main(chars, features, daily_ret)
  write.csv(pf, "output.csv", row.names = FALSE)
}