1#!/usr/bin/env Rscript2# Cross-check: rugarch (version-locked) log-likelihoods and conditional variances on the shared fixture3# tests/replication/fixture_returns.csv, written by tests/replication/test_rugarch.py. Output: JSON with per-model4# fitted parameters, log-likelihood and the last 5 conditional variances (fixed, no estimation differences involved5# beyond the optimizer). Usage: Rscript rugarch_crosscheck.R fixture.csv out.json6suppressPackageStartupMessages({library(rugarch); library(jsonlite)})7a <- commandArgs(trailingOnly = TRUE)8d <- read.csv(a[1]); r <- d$r; x <- d$x9out <- list(rugarch_version = as.character(packageVersion("rugarch")))10fitone <- function(name, spec, rr, realized = NULL) {11 fit <- tryCatch({ if (is.null(realized)) ugarchfit(spec, rr, solver = "hybrid") else ugarchfit(spec, rr, solver = "hybrid", realizedVol = xts::xts(realized, order.by = as.Date("2000-01-01") + seq_along(realized))) }, error = function(e) NULL)12 if (is.null(fit) || fit@fit$convergence != 0) return(list(converged = FALSE))13 list(converged = TRUE, loglik = likelihood(fit), coef = as.list(coef(fit)), sigma2_last5 = as.numeric(tail(sigma(fit)^2, 5)),14 fc = as.numeric(sigma(ugarchforecast(fit, n.ahead = 22))^2))15}16m <- list(mean.model = list(armaOrder = c(0, 0), include.mean = TRUE))17out$garch11_t <- fitone("garch11", ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)18out$gjr11_t <- fitone("gjr11", ugarchspec(variance.model = list(model = "gjrGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)19out$egarch11_t <- fitone("egarch11", ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)20out$aparch11_t <- fitone("aparch11", ugarchspec(variance.model = list(model = "apARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)21out$figarch11_t <- fitone("figarch11", ugarchspec(variance.model = list(model = "fiGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)22out$igarch11_t <- fitone("igarch11", ugarchspec(variance.model = list(model = "iGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)23out$cgarch11_t <- fitone("cgarch11", ugarchspec(variance.model = list(model = "csGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)24out$nagarch11_t <- fitone("nagarch11", ugarchspec(variance.model = list(model = "fGARCH", submodel = "NAGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)25out$tgarch11_t <- fitone("tgarch11", ugarchspec(variance.model = list(model = "fGARCH", submodel = "TGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "std"), r)26out$realgarch11_n <- fitone("realgarch11", ugarchspec(variance.model = list(model = "realGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "norm"), r, realized = sqrt(x))27out$garch11_skewt <- fitone("garch11skt", ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = m$mean.model, distribution.model = "sstd"), r)28write_json(out, a[2], auto_unbox = TRUE, digits = NA)29# ---- extensions (A-03): TGARCH arch counterpart computed on the Python side; GARCH-X with external regressor (return equation);30# GAS-t with inverse-information scaling (R package GAS); MS-GARCH two-regime normal (R package MSGARCH) with 22-step forecasts.31out2 <- list()32xr <- matrix(c(NA, head(x, -1)), ncol = 1) # lagged realized measure (x_{t-1}); first row NA -> drop first obs for the fit33spec_x <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1), external.regressors = xr[-1, , drop = FALSE]), mean.model = m$mean.model, distribution.model = "norm")34fx <- tryCatch(ugarchfit(spec_x, r[-1], solver = "hybrid"), error = function(e) NULL)35if (!is.null(fx) && fx@fit$convergence == 0) out2$garchx_return_eq_n <- list(converged = TRUE, loglik = likelihood(fx), coef = as.list(coef(fx))) else out2$garchx_return_eq_n <- list(converged = FALSE)36if (requireNamespace("GAS", quietly = TRUE)) {37 suppressPackageStartupMessages(library(GAS))38 gs <- UniGASSpec(Dist = "std", ScalingType = "Inv", GASPar = list(location = FALSE, scale = TRUE, shape = FALSE))39 gf <- tryCatch(UniGASFit(gs, r), error = function(e) NULL)40 out2$gas_t <- tryCatch(if (!is.null(gf) && length(gf@Estimates$lLik) > 0) list(converged = TRUE, loglik = as.numeric(gf@Estimates$lLik)[1], coef = as.list(gf@Estimates$lParList), note = "GAS package parameterizes the t scale, not the variance: inverse-information scaling is not invariant, so this is a reference value, not the same model") else list(converged = FALSE, note = "no likelihood returned"), error = function(e) list(converged = FALSE, error = conditionMessage(e)))41}42if (requireNamespace("MSGARCH", quietly = TRUE)) {43 suppressPackageStartupMessages(library(MSGARCH))44 ms <- CreateSpec(variance.spec = list(model = c("sGARCH", "sGARCH")), distribution.spec = list(distribution = c("norm", "norm")), switch.spec = list(do.mix = FALSE))45 mf <- tryCatch(FitML(ms, data = r - mean(r)), error = function(e) NULL)46 out2$msgarch2_n <- tryCatch(if (!is.null(mf)) { pr <- predict(mf, nahead = 22, do.return.draw = FALSE)47 list(converged = TRUE, loglik = as.numeric(mf$loglik), par = as.list(mf$par), vol_forecast = as.numeric(pr$vol), mean_removed = mean(r)) } else list(converged = FALSE),48 error = function(e) list(converged = FALSE, error = conditionMessage(e)))49}50out <- c(out, out2)51write_json(out, a[2], auto_unbox = TRUE, digits = NA)52