#!/usr/bin/env Rscript
# Reference values from R: highfrequency (realized measures, HARmodel) and rugarch (GARCH family) on the committed
# synthetic fixtures. Output: datasets/fixtures/reference_r.json. Versions are recorded in the JSON.
suppressPackageStartupMessages({ library(highfrequency); library(rugarch); library(xts); library(jsonlite) })
root <- normalizePath(file.path(dirname(sub("--file=", "", grep("--file=", commandArgs(), value = TRUE)[1])), ".."))
fx <- file.path(root, "datasets", "fixtures")
out <- list(versions = list(R = R.version.string, highfrequency = as.character(packageVersion("highfrequency")), rugarch = as.character(packageVersion("rugarch"))))

# ---- realized measures per session on the synthetic 1-minute bars (regular session, log returns) -------------------
bars <- read.csv(file.path(fx, "synthetic_bars_1m.csv"), stringsAsFactors = FALSE)
bars$ts <- as.POSIXct(bars$timestamp_utc, format = "%Y-%m-%dT%H:%M:%SZ", tz = "UTC")
bars$et <- format(bars$ts, tz = "America/New_York", format = "%H:%M")
bars$sess <- as.Date(format(bars$ts, tz = "America/New_York", format = "%Y-%m-%d"))
early <- as.Date("2024-11-29")
rm <- list()
for (sym in unique(bars$symbol)) {
  b <- bars[bars$symbol == sym, ]
  for (d in as.character(unique(b$sess))) {
    s <- b[b$sess == as.Date(d), ]; close_min <- if (as.Date(d) == early) "13:00" else "16:00"
    s <- s[s$et >= "09:30" & s$et <= close_min, ]
    p <- xts(s$close, order.by = s$ts); r <- diff(log(p))[-1]
    rm[[paste(sym, d)]] <- list(symbol = sym, session = d, n_returns = length(r),
      rv = as.numeric(rCov(r, makeReturns = FALSE)), bpv = as.numeric(rBPCov(r, makeReturns = FALSE)),
      rsv = as.list(as.numeric(rSVar(r, makeReturns = FALSE))), rq = as.numeric(rQuar(r, makeReturns = FALSE)),
      bns_jump_stat = tryCatch(as.numeric(BNSjumpTest(r, makeReturns = FALSE)$ztest), error = function(e) NA))
  }
}
out$realized <- unname(rm)

# ---- HARmodel on the daily RV fixture -----------------------------------------------------------------------------
h <- read.csv(file.path(fx, "har_rv.csv")); dates <- seq(as.Date("2018-01-01"), by = "day", length.out = nrow(h))
rvx <- xts(h$rv, order.by = dates)
hm <- HARmodel(rvx, periods = c(1, 5, 22), RVest = c("rCov"), type = "HAR", h = 1, transform = NULL, inputType = "RM")
out$har_h1 <- list(coef = as.list(unname(coef(hm))), names = names(coef(hm)), n = length(hm$residuals), ssr = sum(hm$residuals^2))
hm5 <- HARmodel(rvx, periods = c(1, 5, 22), RVest = c("rCov"), type = "HAR", h = 5, transform = NULL, inputType = "RM")
out$har_h5 <- list(coef = as.list(unname(coef(hm5))), n = length(hm5$residuals), ssr = sum(hm5$residuals^2))
hq <- HARmodel(xts(cbind(h$rv, h$rq), order.by = dates), periods = c(1, 5, 22), periodsQ = c(1), RVest = c("rCov", "rQuar"), type = "HARQ", h = 1, transform = NULL, inputType = "RM")
out$harq_h1 <- list(coef = as.list(unname(coef(hq))), names = names(coef(hq)), n = length(hq$residuals), ssr = sum(hq$residuals^2))

# ---- rugarch on the simulated GARCH returns ------------------------------------------------------------------------
r <- read.csv(file.path(fx, "garch_returns.csv"))$r
fit_one <- function(model, dist) {
  spec <- ugarchspec(variance.model = list(model = model, garchOrder = c(1, 1)), mean.model = list(armaOrder = c(0, 0), include.mean = TRUE), distribution.model = dist)
  f <- ugarchfit(spec, r, solver = "hybrid")
  list(coef = as.list(coef(f)), se = as.list(f@fit$matcoef[, 2]), robust_se = as.list(f@fit$robust.matcoef[, 2]), loglik = as.numeric(likelihood(f)), converged = f@fit$convergence == 0,
       forecast_h5 = as.numeric(sigma(ugarchforecast(f, n.ahead = 5)))^2)
}
out$garch <- list(sGARCH.norm = fit_one("sGARCH", "norm"), gjrGARCH.norm = fit_one("gjrGARCH", "norm"), eGARCH.norm = fit_one("eGARCH", "norm"), sGARCH.std = fit_one("sGARCH", "std"))
write(toJSON(out, auto_unbox = TRUE, digits = 17, pretty = TRUE), file.path(fx, "reference_r.json"))
cat("wrote reference_r.json\n")