1#!/usr/bin/env Rscript2# Reference values from R: highfrequency (realized measures, HARmodel) and rugarch (GARCH family) on the committed3# synthetic fixtures. Output: datasets/fixtures/reference_r.json. Versions are recorded in the JSON.4suppressPackageStartupMessages({ library(highfrequency); library(rugarch); library(xts); library(jsonlite) })5root <- normalizePath(file.path(dirname(sub("--file=", "", grep("--file=", commandArgs(), value = TRUE)[1])), ".."))6fx <- file.path(root, "datasets", "fixtures")7out <- list(versions = list(R = R.version.string, highfrequency = as.character(packageVersion("highfrequency")), rugarch = as.character(packageVersion("rugarch"))))89# ---- realized measures per session on the synthetic 1-minute bars (regular session, log returns) -------------------10bars <- read.csv(file.path(fx, "synthetic_bars_1m.csv"), stringsAsFactors = FALSE)11bars$ts <- as.POSIXct(bars$timestamp_utc, format = "%Y-%m-%dT%H:%M:%SZ", tz = "UTC")12bars$et <- format(bars$ts, tz = "America/New_York", format = "%H:%M")13bars$sess <- as.Date(format(bars$ts, tz = "America/New_York", format = "%Y-%m-%d"))14early <- as.Date("2024-11-29")15rm <- list()16for (sym in unique(bars$symbol)) {17 b <- bars[bars$symbol == sym, ]18 for (d in as.character(unique(b$sess))) {19 s <- b[b$sess == as.Date(d), ]; close_min <- if (as.Date(d) == early) "13:00" else "16:00"20 s <- s[s$et >= "09:30" & s$et <= close_min, ]21 p <- xts(s$close, order.by = s$ts); r <- diff(log(p))[-1]22 rm[[paste(sym, d)]] <- list(symbol = sym, session = d, n_returns = length(r),23 rv = as.numeric(rCov(r, makeReturns = FALSE)), bpv = as.numeric(rBPCov(r, makeReturns = FALSE)),24 rsv = as.list(as.numeric(rSVar(r, makeReturns = FALSE))), rq = as.numeric(rQuar(r, makeReturns = FALSE)),25 bns_jump_stat = tryCatch(as.numeric(BNSjumpTest(r, makeReturns = FALSE)$ztest), error = function(e) NA))26 }27}28out$realized <- unname(rm)2930# ---- HARmodel on the daily RV fixture -----------------------------------------------------------------------------31h <- read.csv(file.path(fx, "har_rv.csv")); dates <- seq(as.Date("2018-01-01"), by = "day", length.out = nrow(h))32rvx <- xts(h$rv, order.by = dates)33hm <- HARmodel(rvx, periods = c(1, 5, 22), RVest = c("rCov"), type = "HAR", h = 1, transform = NULL, inputType = "RM")34out$har_h1 <- list(coef = as.list(unname(coef(hm))), names = names(coef(hm)), n = length(hm$residuals), ssr = sum(hm$residuals^2))35hm5 <- HARmodel(rvx, periods = c(1, 5, 22), RVest = c("rCov"), type = "HAR", h = 5, transform = NULL, inputType = "RM")36out$har_h5 <- list(coef = as.list(unname(coef(hm5))), n = length(hm5$residuals), ssr = sum(hm5$residuals^2))37hq <- 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")38out$harq_h1 <- list(coef = as.list(unname(coef(hq))), names = names(coef(hq)), n = length(hq$residuals), ssr = sum(hq$residuals^2))3940# ---- rugarch on the simulated GARCH returns ------------------------------------------------------------------------41r <- read.csv(file.path(fx, "garch_returns.csv"))$r42fit_one <- function(model, dist) {43 spec <- ugarchspec(variance.model = list(model = model, garchOrder = c(1, 1)), mean.model = list(armaOrder = c(0, 0), include.mean = TRUE), distribution.model = dist)44 f <- ugarchfit(spec, r, solver = "hybrid")45 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,46 forecast_h5 = as.numeric(sigma(ugarchforecast(f, n.ahead = 5)))^2)47}48out$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"))49write(toJSON(out, auto_unbox = TRUE, digits = 17, pretty = TRUE), file.path(fx, "reference_r.json"))50cat("wrote reference_r.json\n")51