spb/metrika Public

1#!/usr/bin/env Rscript2#3#  generate.R — Metrika4#5#  Author:  Simon-Pierre Boucher6#  Contact: contact@spboucher.ai7#  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8#9#  Generates the golden numerical fixtures for MetrikaKit estimator tests10#  (CLAUDE.md §9). Writes a reference dataset (CSV) plus expected values11#  (TSV, key<TAB>value at 17 significant digits) into the Swift test12#  resources at MetrikaKit/Tests/MetrikaKitTests/Fixtures/.13#14#  All expected values are computed from the CSV **after** a write/read15#  round-trip so Swift and R consume bit-identical inputs.16#1718out_dir <- file.path(19  dirname(dirname(normalizePath(sub("--file=", "", grep("--file=", commandArgs(FALSE), value = TRUE))))),20  "..", "MetrikaKit", "Tests", "MetrikaKitTests", "Fixtures"21)22out_dir <- normalizePath(out_dir, mustWork = FALSE)23dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)2425# ---------------------------------------------------------------- dataset26set.seed(42)27n <- 6028region  <- rep(1:3, length.out = n)29firm_id <- rep(1:12, each = 5)30price   <- round(runif(n, 5, 20), 4)31noise   <- round(rnorm(n, 0, 0.15), 6)32log_rev <- 4 - 0.08 * price + 0.10 * (region == 2) + 0.20 * (region == 3) + noise33revenue <- round(exp(log_rev), 6)3435# GLM outcomes — drawn AFTER all draws above so earlier fixture values36# stay bit-identical when regenerating.37purchase <- rbinom(n, 1, plogis(2 - 0.20 * price))38orders   <- rpois(n, exp(1.6 - 0.09 * price))3940# Instruments for the 2SLS fixture (correlated with price by41# construction), drawn after everything above.42z1 <- round(price + rnorm(n, 0, 2), 4)43z2 <- round(0.5 * price + rnorm(n, 0, 3), 4)4445# Inject missing values to exercise listwise deletion.46price_missing <- price47price_missing[c(7, 23, 41)] <- NA4849data <- data.frame(50  revenue = revenue,51  price = price_missing,52  region = region,53  firm_id = firm_id,54  purchase = purchase,55  orders = orders,56  z1 = z1,57  z2 = z258)59csv_path <- file.path(out_dir, "regression.csv")60write.csv(data, csv_path, row.names = FALSE, quote = FALSE, na = "")6162# Round-trip: recompute everything from what was actually written.63data <- read.csv(csv_path)6465sink_path <- file.path(out_dir, "expected.tsv")66lines <- character(0)67emit <- function(key, value) {68  lines <<- c(lines, sprintf("%s\t%.17g", key, value))69}7071# ------------------------------------------------------------- summarize72x <- data$revenue73emit("sum_revenue_n",    length(x))74emit("sum_revenue_mean", mean(x))75emit("sum_revenue_var",  var(x))76emit("sum_revenue_sd",   sd(x))77emit("sum_revenue_min",  min(x))78emit("sum_revenue_max",  max(x))79m2 <- mean((x - mean(x))^2)80m3 <- mean((x - mean(x))^3)81m4 <- mean((x - mean(x))^4)82emit("sum_revenue_skewness", m3 / m2^1.5)   # Stata definition83emit("sum_revenue_kurtosis", m4 / m2^2)8485# ------------------------------------------------------------ regression86# Estimation sample: listwise deletion on price.87complete <- !is.na(data$price)88d <- data[complete, ]89d$log_rev <- log(d$revenue)90n_est <- nrow(d)9192X <- cbind(price = d$price, `_cons` = 1)93y <- d$log_rev94k <- ncol(X)9596fit <- lm(log_rev ~ price, data = d)97beta <- coef(fit)                      # (Intercept), price98resid <- residuals(fit)99XtXinv <- solve(t(X) %*% X)100df_r <- n_est - k101sigma2 <- sum(resid^2) / df_r102103emit("ols_n", n_est)104emit("ols_dropped", sum(!complete))105emit("ols_b_price", beta["price"])106emit("ols_b_cons",  beta["(Intercept)"])107emit("ols_r2",   summary(fit)$r.squared)108emit("ols_r2a",  summary(fit)$adj.r.squared)109emit("ols_rmse", sqrt(sigma2))110emit("ols_F",    summary(fit)$fstatistic[["value"]])111112V_classical <- sigma2 * XtXinv113emit("ols_se_classical_price", sqrt(V_classical["price", "price"]))114emit("ols_se_classical_cons",  sqrt(V_classical["_cons", "_cons"]))115t_price <- beta["price"] / sqrt(V_classical["price", "price"])116emit("ols_t_price", t_price)117emit("ols_p_price", 2 * pt(-abs(t_price), df_r))118crit <- qt(0.975, df_r)119emit("ols_ci_lower_price", beta["price"] - crit * sqrt(V_classical["price", "price"]))120emit("ols_ci_upper_price", beta["price"] + crit * sqrt(V_classical["price", "price"]))121122# HC0–HC3 sandwich estimators (manual, no packages).123h <- hatvalues(fit)124sandwich_se <- function(w) {125  meat <- t(X * w) %*% X126  V <- XtXinv %*% meat %*% XtXinv127  sqrt(diag(V))128}129se_hc0 <- sandwich_se(resid^2)130se_hc1 <- sandwich_se(resid^2 * n_est / df_r)131se_hc2 <- sandwich_se(resid^2 / (1 - h))132se_hc3 <- sandwich_se(resid^2 / (1 - h)^2)133emit("ols_se_hc0_price", se_hc0["price"])134emit("ols_se_hc0_cons",  se_hc0["_cons"])135emit("ols_se_hc1_price", se_hc1["price"])136emit("ols_se_hc1_cons",  se_hc1["_cons"])137emit("ols_se_hc2_price", se_hc2["price"])138emit("ols_se_hc2_cons",  se_hc2["_cons"])139emit("ols_se_hc3_price", se_hc3["price"])140emit("ols_se_hc3_cons",  se_hc3["_cons"])141142# Cluster-robust with Stata regress small-sample factor.143cl <- d$firm_id144u <- rowsum(X * resid, cl)145G <- length(unique(cl))146meat_cl <- (G / (G - 1)) * ((n_est - 1) / df_r) * (t(u) %*% u)147V_cl <- XtXinv %*% meat_cl %*% XtXinv148emit("ols_G", G)149emit("ols_se_cluster_price", sqrt(V_cl["price", "price"]))150emit("ols_se_cluster_cons",  sqrt(V_cl["_cons", "_cons"]))151t_cl <- beta["price"] / sqrt(V_cl["price", "price"])152emit("ols_p_cluster_price", 2 * pt(-abs(t_cl), G - 1))153154# Factor-variable regression: log_rev ~ price + i.region (base = 1).155fit2 <- lm(log_rev ~ price + factor(region), data = d)156b2 <- coef(fit2)157emit("ols2_b_price",   b2[["price"]])158emit("ols2_b_region2", b2[["factor(region)2"]])159emit("ols2_b_region3", b2[["factor(region)3"]])160emit("ols2_b_cons",    b2[["(Intercept)"]])161emit("ols2_r2", summary(fit2)$r.squared)162163# ------------------------------------------------------------------ GLMs164# Fit tightly (epsilon 1e-12) so R and Swift land on the same optimum to165# well past the 1e-10 test tolerance. Estimation sample = complete cases.166ctrl <- glm.control(epsilon = 1e-12, maxit = 100)167168# Expected information (X'WX)^-1 evaluated AT the converged coefficients.169# R's vcov(fit) instead reuses the weights of the last IRLS step (at the170# second-to-last iterate), so it is only ~1e-6-accurate by its own171# stopping rule — not good enough for 1e-10 fixtures.172glm_bread <- function(fit) {173  Xg <- model.matrix(fit)174  eta <- as.vector(Xg %*% coef(fit))175  fam <- fit$family176  mu <- fam$linkinv(eta)177  W <- fam$mu.eta(eta)^2 / fam$variance(mu)178  solve(t(Xg * W) %*% Xg)179}180181glm_bread_sandwich <- function(fit, score_scale, cl = NULL) {182  Xg <- model.matrix(fit)183  u <- Xg * score_scale184  bread <- glm_bread(fit)185  if (is.null(cl)) {186    meat <- t(u) %*% u187  } else {188    ug <- rowsum(u, cl)189    G <- nrow(ug)190    meat <- (G / (G - 1)) * (t(ug) %*% ug)191  }192  sqrt(diag(bread %*% meat %*% bread))193}194195lg <- glm(purchase ~ price, data = d, family = binomial(), control = ctrl)196V_lg <- glm_bread(lg)197emit("logit_b_price", coef(lg)[["price"]])198emit("logit_b_cons",  coef(lg)[["(Intercept)"]])199emit("logit_se_price", sqrt(V_lg["price", "price"]))200emit("logit_se_cons",  sqrt(V_lg["(Intercept)", "(Intercept)"]))201emit("logit_ll",  as.numeric(logLik(lg)))202lg0 <- glm(purchase ~ 1, data = d, family = binomial(), control = ctrl)203emit("logit_ll0", as.numeric(logLik(lg0)))204emit("logit_chi2", 2 * (as.numeric(logLik(lg)) - as.numeric(logLik(lg0))))205emit("logit_chi2p", pchisq(206  2 * (as.numeric(logLik(lg)) - as.numeric(logLik(lg0))), 1, lower.tail = FALSE207))208emit("logit_se_hc0_price", glm_bread_sandwich(lg, d$purchase - fitted(lg))["price"])209emit("logit_se_cluster_price",210     glm_bread_sandwich(lg, d$purchase - fitted(lg), cl = d$firm_id)["price"])211212pr <- glm(purchase ~ price, data = d,213          family = binomial(link = "probit"), control = ctrl)214emit("probit_b_price", coef(pr)[["price"]])215emit("probit_b_cons",  coef(pr)[["(Intercept)"]])216emit("probit_se_price", sqrt(glm_bread(pr)["price", "price"]))217emit("probit_ll", as.numeric(logLik(pr)))218219ps <- glm(orders ~ price, data = d, family = poisson(), control = ctrl)220emit("pois_b_price", coef(ps)[["price"]])221emit("pois_b_cons",  coef(ps)[["(Intercept)"]])222emit("pois_se_price", sqrt(glm_bread(ps)["price", "price"]))223emit("pois_ll", as.numeric(logLik(ps)))224emit("pois_se_hc0_price", glm_bread_sandwich(ps, d$orders - fitted(ps))["price"])225226# ----------------------------------------------------- fixed effects (FE)227# Within estimator, Stata conventions: demean within firm, add back grand228# means (so _cons is reported), df = N − K − G. Cluster VCE on the panel229# with the G/(G−1) factor and t on G−1 df.230fe_within <- function(v, g) {231  gm <- ave(v, g)             # group means232  v - gm + mean(v)233}234g <- d$firm_id235Nw <- nrow(d)236yt <- fe_within(d$log_rev, g)237xt_p <- fe_within(d$price, g)238Xw <- cbind(price = xt_p, `_cons` = 1)239G <- length(unique(g))240Kw <- 1241df_fe <- Nw - Kw - G242243fe_fit <- lm.fit(Xw, yt)244b_fe <- fe_fit$coefficients245res_fe <- fe_fit$residuals246XtXinv_w <- solve(t(Xw) %*% Xw)247sigma2_fe <- sum(res_fe^2) / df_fe248emit("fe_N", Nw)249emit("fe_G", G)250emit("fe_b_price", b_fe[["price"]])251emit("fe_b_cons",  b_fe[["_cons"]])252emit("fe_se_price", sqrt(sigma2_fe * XtXinv_w["price", "price"]))253yd <- d$log_rev - ave(d$log_rev, g)254emit("fe_r2_within", 1 - sum(res_fe^2) / sum(yd^2))255emit("fe_rmse", sqrt(sigma2_fe))256257# Cluster on the panel: scores use DEMEANED x (constant column as-is).258xd <- d$price - ave(d$price, g)259u_fe <- rowsum(cbind(xd * res_fe, res_fe), g)260meat_fe <- (G / (G - 1)) * (t(u_fe) %*% u_fe)261V_fe_cl <- XtXinv_w %*% meat_fe %*% XtXinv_w262emit("fe_se_cluster_price", sqrt(V_fe_cl[1, 1]))263t_fe_cl <- b_fe[["price"]] / sqrt(V_fe_cl[1, 1])264emit("fe_p_cluster_price", 2 * pt(-abs(t_fe_cl), G - 1))265266# ------------------------------------------------------------- 2SLS (IV)267# ivregress 2sls log_rev (price = z1 z2), Stata `small` convention:268# residuals from the ORIGINAL regressors, sigma2 = u'u/(N-K), t stats.269Xiv <- cbind(price = d$price, `_cons` = 1)270Ziv <- cbind(z1 = d$z1, z2 = d$z2, `_cons` = 1)271PZ <- Ziv %*% solve(t(Ziv) %*% Ziv) %*% t(Ziv)272Xhat <- PZ %*% Xiv273XhXhinv <- solve(t(Xhat) %*% Xhat)274b_iv <- XhXhinv %*% t(Xhat) %*% d$log_rev275u_iv <- d$log_rev - Xiv %*% b_iv276df_iv <- Nw - 2277sigma2_iv <- sum(u_iv^2) / df_iv278V_iv <- sigma2_iv * XhXhinv279emit("iv_N", Nw)280emit("iv_b_price", b_iv[1, 1])281emit("iv_b_cons",  b_iv[2, 1])282emit("iv_se_price", sqrt(V_iv[1, 1]))283emit("iv_se_cons",  sqrt(V_iv[2, 2]))284t_iv <- b_iv[1, 1] / sqrt(V_iv[1, 1])285emit("iv_p_price", 2 * pt(-abs(t_iv), df_iv))286emit("iv_r2", 1 - sum(u_iv^2) / sum((d$log_rev - mean(d$log_rev))^2))287emit("iv_rmse", sqrt(sigma2_iv))288289# Robust (HC1 with the small factor N/(N-K)) built on projected X.290meat_iv <- t(Xhat * as.vector(u_iv^2)) %*% Xhat * (Nw / df_iv)291V_iv_r <- XhXhinv %*% meat_iv %*% XhXhinv292emit("iv_se_hc1_price", sqrt(V_iv_r[1, 1]))293294# ---------------------------------------------------- margins (AME, delta)295# Average marginal effects with delta-method SEs at the converged296# coefficients, mirroring the engine's formulas exactly.297ame_glm <- function(fit, var, fprime, fsecond) {298  Xg <- model.matrix(fit)299  eta <- as.vector(Xg %*% coef(fit))300  fp <- fprime(eta)301  fs <- fsecond(eta)302  b <- coef(fit)[[var]]303  ame <- b * mean(fp)304  grad <- colMeans(Xg * fs) * b305  grad[[var]] <- grad[[var]] + mean(fp)306  V <- glm_bread(fit)307  se <- sqrt(as.numeric(t(grad) %*% V %*% grad))308  c(ame = ame, se = se)309}310311m_lg <- ame_glm(312  lg, "price",313  function(e) plogis(e) * (1 - plogis(e)),314  function(e) plogis(e) * (1 - plogis(e)) * (1 - 2 * plogis(e))315)316emit("margins_logit_price", m_lg[["ame"]])317emit("margins_logit_se",    m_lg[["se"]])318319m_ps <- ame_glm(ps, "price", exp, exp)320emit("margins_pois_price", m_ps[["ame"]])321emit("margins_pois_se",    m_ps[["se"]])322323# ------------------------------------------------------------ elastic net324# glmnet at fixed lambda, tight threshold. price/z1/z2 are correlated by325# construction, so the l1 penalty must actually choose among them.326if (requireNamespace("glmnet", quietly = TRUE)) {327  Xen <- as.matrix(d[, c("price", "z1", "z2", "orders")])328  yen <- d$log_rev329330  fit_lasso <- glmnet::glmnet(331    Xen, yen, alpha = 1, lambda = 0.05,332    thresh = 1e-15, maxit = 1e7, standardize = TRUE333  )334  cl <- as.vector(coef(fit_lasso))335  emit("lasso_b_cons",   cl[1])336  emit("lasso_b_price",  cl[2])337  emit("lasso_b_z1",     cl[3])338  emit("lasso_b_z2",     cl[4])339  emit("lasso_b_orders", cl[5])340341  fit_enet <- glmnet::glmnet(342    Xen, yen, alpha = 0.4, lambda = 0.02,343    thresh = 1e-15, maxit = 1e7, standardize = TRUE344  )345  ce <- as.vector(coef(fit_enet))346  emit("enet_b_cons",   ce[1])347  emit("enet_b_price",  ce[2])348  emit("enet_b_z1",     ce[3])349  emit("enet_b_z2",     ce[4])350  emit("enet_b_orders", ce[5])351} else {352  cat("glmnet not installed — skipping elastic-net fixtures\n")353}354355# ------------------------------------------------------- gradient boosting356# xgboost exact-greedy at fixed settings, no subsampling — deterministic.357# Complete-case features only (missing-value default-direction search358# differs between implementations).359if (requireNamespace("xgboost", quietly = TRUE)) {360  # Full sample: these variables have no missing values, so the engine's361  # listwise deletion keeps all 60 rows.362  Xgb <- as.matrix(data[, c("z1", "z2", "orders")])363  ygb <- log(data$revenue)364  base <- mean(ygb)365  booster <- xgboost::xgboost(366    x = Xgb, y = ygb,367    nrounds = 20, max_depth = 3, learning_rate = 0.3,368    reg_lambda = 1, reg_alpha = 0, min_child_weight = 1,369    subsample = 1, colsample_bytree = 1,370    tree_method = "exact", base_score = base,371    nthread = 1, verbosity = 0372  )373  pred <- predict(booster, Xgb)374  emit("gbm_rmse", sqrt(mean((ygb - pred)^2)))375  emit("gbm_pred_1",  pred[1])376  emit("gbm_pred_17", pred[17])377  emit("gbm_pred_42", pred[42])378  emit("gbm_pred_57", pred[57])379} else {380  cat("xgboost not installed — skipping boosting fixtures\n")381}382383# -------------------------------------------------------------- correlate384emit("corr_rev_price", cor(d$revenue, d$price))385emit("corr_rev_logrev", cor(d$revenue, d$log_rev))386387# --------------------------------------------------------- distributions388emit("dist_pchisq_3p8_1", pchisq(3.8, 1))389emit("dist_pchisq_25_4",  pchisq(25, 4))390emit("dist_qnorm_0p975",  qnorm(0.975))391emit("dist_pt_2p5_df10",   pt(2.5, 10))392emit("dist_pt_m1p3_df3",   pt(-1.3, 3))393emit("dist_pt_0p05_df57",  pt(0.05, 57))394emit("dist_pf_3p7_2_30",   pf(3.7, 2, 30))395emit("dist_pf_0p5_5_100",  pf(0.5, 5, 100))396emit("dist_qt_0p975_df12", qt(0.975, 12))397emit("dist_qt_0p995_df4",  qt(0.995, 4))398emit("dist_pnorm_1p64",    pnorm(1.64))399400# ------------------------------------------------------------- .dta files401# haven-written fixtures for the native DTA reader: same dataset plus a402# string column, in formats 118 (version 14) and 117 (version 13).403if (requireNamespace("haven", quietly = TRUE)) {404  dta <- data405  dta$firm_name <- paste0("firm_", data$firm_id)406  dta$firm_name[5] <- ""   # string missing407  haven::write_dta(dta, file.path(out_dir, "regression_v118.dta"), version = 14)408  haven::write_dta(dta, file.path(out_dir, "regression_v117.dta"), version = 13)409  cat("wrote .dta fixtures (117, 118)\n")410} else {411  cat("haven not installed — skipping .dta fixtures\n")412}413414writeLines(lines, sink_path)415cat("wrote", csv_path, "\n")416cat("wrote", sink_path, "with", length(lines), "expected values\n")417