# backtest_levetf.R: what leveraged ETF "decay" actually is.
#
# The retail claim: 3x ETFs bleed value over time, so they're day-trading
# vehicles only. The counter-claim: TQQQ crushed QQQ over the last decade, so
# decay is a myth. Both are wrong in the same way. They treat a
# path-dependent compounding effect as if it were a constant fee.
#
# Decomposition, per pair:
#   (a) naive     = L x (underlying total return)      <- what buyers expect
#   (b) mechanism = prod(1 + L * daily underlying ret) <- daily rebalance, no costs
#   (c) actual    = prod(1 + daily leveraged ETF ret)  <- what you get
#
#   (a) -> (b) is the compounding/path effect. Cuts both ways.
#   (b) -> (c) is fees + financing + tracking error. Only ever costs you.
#
#   Rscript backtest_levetf.R
#
# This is the published copy of the script behind
# adjustedclose.com/leveraged-etf-decay. It differs from the research
# original in exactly two ways, so it runs standalone: the repo-internal
# data loader (a thin on-disk-cached wrapper around quantmod::getSymbols)
# is inlined below, and the results/ output directory is created if
# missing. Every number and figure on the page comes from this analysis.
#
# Dependencies: R with the quantmod and xts packages.
# Data: split/dividend-adjusted daily bars from Yahoo Finance via
# quantmod. Free, no API key. Yahoo's endpoint is unofficial; if it
# changes, point yahoo_daily() at any source of adjusted daily OHLC.

suppressPackageStartupMessages({
  library(quantmod)
  library(xts)
})

yahoo_daily <- function(symbol, from = "2000-01-01") {
  getSymbols(toupper(symbol), src = "yahoo", from = from, auto.assign = FALSE)
}

PAIRS <- list(
  list(lev = "TQQQ", base = "QQQ",  L = 3),
  list(lev = "UPRO", base = "SPY",  L = 3),
  list(lev = "SOXL", base = "SOXX", L = 3),
  list(lev = "TNA",  base = "IWM",  L = 3),
  list(lev = "SSO",  base = "SPY",  L = 2),
  list(lev = "QLD",  base = "QQQ",  L = 2)
)

FROM <- "2009-01-01"
fmt  <- function(x, d = 1) formatC(x, format = "f", digits = d, big.mark = ",")

# ---- load ------------------------------------------------------------------

get_ret <- function(sym) {
  x <- yahoo_daily(sym, from = FROM)
  r <- dailyReturn(Ad(x))
  colnames(r) <- sym
  r
}

cache <- new.env()
ret_of <- function(sym) {
  if (is.null(cache[[sym]])) cache[[sym]] <- get_ret(sym)
  cache[[sym]]
}

# ---- per-pair decomposition ------------------------------------------------

rows <- list()
detail <- list()

for (p in PAIRS) {
  rl <- ret_of(p$lev)
  rb <- ret_of(p$base)
  m  <- merge(rl, rb, join = "inner")
  m  <- m[complete.cases(m), ]
  m  <- m[-1, ]   # first bar of dailyReturn is 0

  lr <- as.numeric(m[, 1])
  br <- as.numeric(m[, 2])
  n  <- length(lr)
  yrs <- n / 252

  actual    <- prod(1 + lr) - 1
  mechanism <- prod(1 + p$L * br) - 1
  base_tot  <- prod(1 + br) - 1
  naive     <- p$L * base_tot

  # annualized geometric
  ann <- function(tot) (1 + tot)^(1 / yrs) - 1

  # daily tracking: is the ETF actually delivering Lx per day?
  fit  <- lm(lr ~ br)
  beta <- unname(coef(fit)[2])
  r2   <- summary(fit)$r.squared

  sig  <- sd(br) * sqrt(252)
  # theoretical annual drag of Lx vs L * (underlying geometric return)
  drag_theory <- 0.5 * p$L * (p$L - 1) * sig^2
  drag_actual <- p$L * ann(base_tot) - ann(mechanism)

  rows[[length(rows) + 1]] <- data.frame(
    pair = sprintf("%s/%s", p$lev, p$base), L = p$L, yrs = yrs,
    base_tot = base_tot, naive = naive, mechanism = mechanism, actual = actual,
    ann_base = ann(base_tot), ann_actual = ann(actual),
    beta = beta, r2 = r2, vol = sig,
    drag_theory = drag_theory, drag_actual = drag_actual,
    cost_gap = ann(mechanism) - ann(actual),
    stringsAsFactors = FALSE
  )
  detail[[p$lev]] <- list(lr = lr, br = br, dates = index(m), L = p$L)
}

res <- do.call(rbind, rows)

cat("\n================ DECOMPOSITION (inception -> today) ================\n\n")
cat(sprintf("%-11s %3s %5s %10s %10s %10s %10s\n",
            "pair", "L", "yrs", "base tot", "naive Lx", "mechanism", "ACTUAL"))
for (i in seq_len(nrow(res))) with(res[i, ], cat(sprintf(
  "%-11s %3d %5.1f %9.0f%% %9.0f%% %9.0f%% %9.0f%%\n",
  pair, L, yrs, base_tot * 100, naive * 100, mechanism * 100, actual * 100)))

cat("\n---- annualized, and where the gap goes ----\n\n")
cat(sprintf("%-11s %9s %9s %8s %8s %9s %9s\n",
            "pair", "ann base", "ann act", "vol", "beta", "drag thr", "cost gap"))
for (i in seq_len(nrow(res))) with(res[i, ], cat(sprintf(
  "%-11s %8.1f%% %8.1f%% %7.1f%% %8.3f %8.1f%% %8.2f%%\n",
  pair, ann_base * 100, ann_actual * 100, vol * 100, beta,
  drag_theory * 100, cost_gap * 100)))

# ---- holding-period analysis ----------------------------------------------
# Over rolling windows of length N: how often does the 3x actually beat
# L x the underlying's return over that same window?

cat("\n\n================ ROLLING WINDOWS: does Lx beat L * base? ================\n")

HORIZONS <- c(21, 63, 126, 252, 504, 756)

roll_tbl <- list()
for (nm in names(detail)) {
  d <- detail[[nm]]
  for (H in HORIZONS) {
    if (length(d$lr) < H + 10) next
    idx <- seq(1, length(d$lr) - H, by = 5)   # every 5th start, overlapping
    win_lev  <- numeric(length(idx))
    win_nv   <- numeric(length(idx))
    win_vol  <- numeric(length(idx))
    for (k in seq_along(idx)) {
      s <- idx[k]; e <- s + H - 1
      bl <- d$br[s:e]
      win_lev[k] <- prod(1 + d$lr[s:e]) - 1
      win_nv[k]  <- d$L * (prod(1 + bl) - 1)
      win_vol[k] <- sd(bl) * sqrt(252)
    }
    edge <- win_lev - win_nv
    roll_tbl[[length(roll_tbl) + 1]] <- data.frame(
      sym = nm, H = H, n = length(idx),
      win_rate = mean(edge > 0),
      med_edge = median(edge),
      p10 = quantile(edge, .10), p90 = quantile(edge, .90),
      stringsAsFactors = FALSE)
  }
}
rt <- do.call(rbind, roll_tbl)

for (nm in unique(rt$sym)) {
  cat(sprintf("\n%s\n", nm))
  cat(sprintf("  %6s %6s %9s %10s %10s %10s\n",
              "days", "n", "beat Lx", "med edge", "p10", "p90"))
  s <- rt[rt$sym == nm, ]
  for (i in seq_len(nrow(s))) with(s[i, ], cat(sprintf(
    "  %6d %6d %8.0f%% %9.1f%% %9.1f%% %9.1f%%\n",
    H, n, win_rate * 100, med_edge * 100, p10 * 100, p90 * 100)))
}

# ---- the actual driver: realized vol during the window --------------------

cat("\n\n================ EDGE BY REALIZED VOL (252-day windows) ================\n")
cat("\nThe claim under test: 'decay' is not a constant. It is a function of\n")
cat("realized volatility over the holding period.\n\n")

qtl <- list()
for (nm in names(detail)) {
  d <- detail[[nm]]
  H <- 252
  if (length(d$lr) < H + 10) next
  idx <- seq(1, length(d$lr) - H, by = 5)
  ev <- data.frame(edge = numeric(length(idx)), vol = numeric(length(idx)),
                   base = numeric(length(idx)))
  for (k in seq_along(idx)) {
    s <- idx[k]; e <- s + H - 1
    bl <- d$br[s:e]
    bt <- prod(1 + bl) - 1
    ev$edge[k] <- (prod(1 + d$lr[s:e]) - 1) - d$L * bt
    ev$vol[k]  <- sd(bl) * sqrt(252)
    ev$base[k] <- bt
  }
  ev$q <- cut(ev$vol, breaks = quantile(ev$vol, seq(0, 1, .2)),
              labels = paste0("Q", 1:5), include.lowest = TRUE)
  cat(sprintf("\n%s  (252-day windows, n=%d)\n", nm, nrow(ev)))
  cat(sprintf("  %-4s %12s %12s %12s %10s\n",
              "vol", "vol range", "med edge", "med base ret", "beat Lx"))
  for (q in levels(ev$q)) {
    s <- ev[ev$q == q, ]
    cat(sprintf("  %-4s %5.0f%%-%5.0f%% %11.1f%% %11.1f%% %9.0f%%\n",
                q, min(s$vol) * 100, max(s$vol) * 100,
                median(s$edge) * 100, median(s$base) * 100,
                mean(s$edge > 0) * 100))
    qtl[[length(qtl) + 1]] <- data.frame(
      sym = nm, q = q, vol_lo = min(s$vol), vol_hi = max(s$vol),
      med_edge = median(s$edge), med_base = median(s$base),
      beat = mean(s$edge > 0), stringsAsFactors = FALSE)
  }
}

# ---- drawdown reality ------------------------------------------------------

cat("\n\n================ DRAWDOWN AND RECOVERY ================\n\n")
cat(sprintf("%-6s %10s %10s %14s %14s\n",
            "sym", "maxDD", "base maxDD", "days underwtr", "worst day"))
for (nm in names(detail)) {
  d  <- detail[[nm]]
  eq <- cumprod(1 + d$lr); pk <- cummax(eq); dd <- eq / pk - 1
  eb <- cumprod(1 + d$br); pb <- cummax(eb); db <- eb / pb - 1
  under <- sum(dd < -0.001)
  cat(sprintf("%-6s %9.1f%% %9.1f%% %13d %13.1f%%\n",
              nm, min(dd) * 100, min(db) * 100, under, min(d$lr) * 100))
}

# ---- common-window control -------------------------------------------------
# The pairs have different inception dates, so the headline table compares
# funds over different windows. TNA (Nov 2008) is 15 months OLDER than TQQQ
# (Feb 2010) and its extra window contains the March 2009 bottom, so the
# unequal samples flatter TNA rather than penalize it. Re-run everything from
# the latest inception in the set (SOXL, 2010-03-11) so all six share one
# identical window and the comparison is apples to apples.

COMMON <- as.Date("2010-03-11")

cat("\n\n================ COMMON WINDOW (all pairs, %s -> today) ================\n\n")
cat(sprintf("Every pair over the identical window starting %s.\n\n", COMMON))

crows <- list()
for (p in PAIRS) {
  m <- merge(ret_of(p$lev), ret_of(p$base), join = "inner")
  m <- m[complete.cases(m), ]
  m <- m[index(m) >= COMMON, ]

  lr <- as.numeric(m[, 1]); br <- as.numeric(m[, 2])
  yrs <- length(lr) / 252
  ann <- function(tot) (1 + tot)^(1 / yrs) - 1

  base_tot <- prod(1 + br) - 1
  actual   <- prod(1 + lr) - 1
  sig      <- sd(br) * sqrt(252)

  crows[[length(crows) + 1]] <- data.frame(
    pair = sprintf("%s/%s", p$lev, p$base), L = p$L,
    base_tot = base_tot, naive = p$L * base_tot, actual = actual,
    ratio = actual / (p$L * base_tot),
    ann_base = ann(base_tot), ann_actual = ann(actual),
    vol = sig, sharpe_ish = ann(base_tot) / sig,
    stringsAsFactors = FALSE)
}
cw <- do.call(rbind, crows)

cat(sprintf("%-11s %3s %10s %10s %10s %8s %9s %8s %8s\n",
            "pair", "L", "base tot", "naive Lx", "ACTUAL", "act/nv",
            "ann base", "vol", "ret/vol"))
for (i in seq_len(nrow(cw))) with(cw[i, ], cat(sprintf(
  "%-11s %3d %9.0f%% %9.0f%% %9.0f%% %8.2f %8.1f%% %7.1f%% %8.2f\n",
  pair, L, base_tot * 100, naive * 100, actual * 100, ratio,
  ann_base * 100, vol * 100, sharpe_ish)))

cat("\nact/nv > 1 = compounding helped; < 1 = compounding hurt.\n")
cat("Sorted by underlying return/vol, the ordering should be monotonic:\n\n")
cwo <- cw[order(-cw$sharpe_ish), ]
for (i in seq_len(nrow(cwo))) with(cwo[i, ], cat(sprintf(
  "  %-11s ret/vol %.2f  ->  act/nv %.2f\n", pair, sharpe_ish, ratio)))

cat("\n")

# ---- exports for the site figures -----------------------------------------
# Growth-of-$1 curves for the hero/anti-hero contrast (TQQQ vs TNA), sampled
# monthly so the rendered SVG stays small. naive as a curve is 1 + L*(cum-1).

crv <- list()
for (nm in c("TQQQ", "TNA")) {
  d  <- detail[[nm]]
  cb <- cumprod(1 + d$br)
  i  <- unique(c(seq(1, length(cb), by = 21), length(cb)))
  crv[[nm]] <- data.frame(
    fund = nm, date = format(d$dates[i]),
    naive = 1 + d$L * (cb[i] - 1),
    mechanism = cumprod(1 + d$L * d$br)[i],
    actual = cumprod(1 + d$lr)[i], stringsAsFactors = FALSE)
}

dir.create("results", showWarnings = FALSE)
write.csv(res, "results/levetf_decomposition.csv", row.names = FALSE)
write.csv(cw,  "results/levetf_common_window.csv", row.names = FALSE)
write.csv(rt,  "results/levetf_rolling.csv", row.names = FALSE)
write.csv(do.call(rbind, crv), "results/levetf_curves.csv", row.names = FALSE)
write.csv(do.call(rbind, qtl), "results/levetf_quintiles.csv", row.names = FALSE)