spb/metrika Public

1//2//  StatsTests.swift3//  Metrika4//5//  Author:  Simon-Pierre Boucher6//  Contact: contact@spboucher.ai7//  Copyright © 2026 Simon-Pierre Boucher. All rights reserved.8//910import Foundation11import Testing12import ZQData13import ZQStats1415/// Numerical tests against R-generated fixtures (Tests/Fixtures/generate.R)16/// at 1e-10 relative tolerance (CLAUDE.md §0).17@Suite("ZQStats vs R fixtures")18struct StatsTests {19    let fixtures: Fixtures20    let y: [Double]                 // log(revenue), estimation sample21    let price: [Double]22    let region: [Double]23    let firmID: [Int]2425    init() async throws {26        self.fixtures = try Fixtures()2728        // Load the fixture CSV through the real data path (DuckDB).29        let store = try ZQDataStore()30        let frame = try await store.load(contentsOf: fixtures.datasetURL)31        let (revenue, revenueMissing) = try frame.requireNumeric("revenue")32        let (priceAll, priceMissing) = try frame.requireNumeric("price")33        let (regionAll, _) = try frame.requireNumeric("region")34        let (firmAll, _) = try frame.requireNumeric("firm_id")3536        // Listwise deletion, mirroring generate.R.37        var y: [Double] = [], price: [Double] = []38        var region: [Double] = [], firm: [Int] = []39        for i in 0..<frame.rowCount where !priceMissing[i] && !revenueMissing[i] {40            y.append(Foundation.log(revenue[i]))41            price.append(priceAll[i])42            region.append(regionAll[i])43            firm.append(Int(firmAll[i]))44        }45        self.y = y46        self.price = price47        self.region = region48        self.firmID = firm49    }5051    @Test("summarize matches R moments")52    func summarize() async throws {53        let store = try ZQDataStore()54        let frame = try await store.load(contentsOf: fixtures.datasetURL)55        let (values, missing) = try frame.requireNumeric("revenue")56        let summary = try #require(ZQSummarize.summary(57            name: "revenue", values: values, missing: missing, detail: true58        ))59        #expect(summary.observationCount == Int(fixtures["sum_revenue_n"]))60        expectClose(summary.mean, fixtures["sum_revenue_mean"], "mean")61        expectClose(summary.variance, fixtures["sum_revenue_var"], "variance")62        expectClose(summary.standardDeviation, fixtures["sum_revenue_sd"], "sd")63        expectClose(summary.minimum, fixtures["sum_revenue_min"], "min")64        expectClose(summary.maximum, fixtures["sum_revenue_max"], "max")65        let detail = try #require(summary.detail)66        expectClose(detail.skewness, fixtures["sum_revenue_skewness"], "skewness")67        expectClose(detail.kurtosis, fixtures["sum_revenue_kurtosis"], "kurtosis")68    }6970    @Test("OLS point estimates and classical inference")71    func classicalOLS() throws {72        let result = try ZQOLS.fit(73            y: y, predictors: [("price", price)], variance: .classical74        )75        #expect(result.observationCount == Int(fixtures["ols_n"]))76        let priceCoef = result.coefficients[0]77        let constant = result.coefficients[1]78        expectClose(priceCoef.estimate, fixtures["ols_b_price"], "b[price]")79        expectClose(constant.estimate, fixtures["ols_b_cons"], "b[_cons]")80        expectClose(priceCoef.standardError, fixtures["ols_se_classical_price"], "se[price]")81        expectClose(constant.standardError, fixtures["ols_se_classical_cons"], "se[_cons]")82        expectClose(priceCoef.tStatistic, fixtures["ols_t_price"], "t[price]")83        expectClose(priceCoef.pValue, fixtures["ols_p_price"], "p[price]")84        expectClose(priceCoef.confidenceLower, fixtures["ols_ci_lower_price"], "ci lower")85        expectClose(priceCoef.confidenceUpper, fixtures["ols_ci_upper_price"], "ci upper")86        expectClose(result.rSquared, fixtures["ols_r2"], "R²")87        expectClose(result.adjustedRSquared, fixtures["ols_r2a"], "adj. R²")88        expectClose(result.rootMSE, fixtures["ols_rmse"], "root MSE")89        expectClose(try #require(result.fStatistic), fixtures["ols_F"], "F")90    }9192    @Test("HC0–HC3 sandwich standard errors", arguments: [93        ("hc0", ZQVarianceEstimator.hc0),94        ("hc1", ZQVarianceEstimator.hc1),95        ("hc2", ZQVarianceEstimator.hc2),96        ("hc3", ZQVarianceEstimator.hc3),97    ])98    func robustSE(_ label: String, _ variance: ZQVarianceEstimator) throws {99        let result = try ZQOLS.fit(100            y: y, predictors: [("price", price)], variance: variance101        )102        expectClose(103            result.coefficients[0].standardError,104            fixtures["ols_se_\(label)_price"], "\(label) se[price]"105        )106        expectClose(107            result.coefficients[1].standardError,108            fixtures["ols_se_\(label)_cons"], "\(label) se[_cons]"109        )110    }111112    @Test("cluster-robust standard errors (Stata small-sample factor)")113    func clusterSE() throws {114        let result = try ZQOLS.fit(115            y: y, predictors: [("price", price)], variance: .cluster(firmID)116        )117        #expect(result.clusterCount == Int(fixtures["ols_G"]))118        expectClose(119            result.coefficients[0].standardError,120            fixtures["ols_se_cluster_price"], "cluster se[price]"121        )122        expectClose(123            result.coefficients[1].standardError,124            fixtures["ols_se_cluster_cons"], "cluster se[_cons]"125        )126        expectClose(127            result.coefficients[0].pValue,128            fixtures["ols_p_cluster_price"], "cluster p[price] (t on G−1 df)"129        )130    }131132    @Test("OLS invariance: permuting observations leaves estimates unchanged")133    func permutationInvariance() throws {134        let base = try ZQOLS.fit(y: y, predictors: [("price", price)])135        var order = Array(0..<y.count)136        // Deterministic shuffle (LCG) — no seeding dependency in tests.137        var state: UInt64 = 88_172_645_463_325_252138        for i in stride(from: order.count - 1, through: 1, by: -1) {139            state = state &* 6_364_136_223_846_793_005 &+ 1_442_695_040_888_963_407140            order.swapAt(i, Int(state % UInt64(i + 1)))141        }142        let permuted = try ZQOLS.fit(143            y: order.map { y[$0] },144            predictors: [("price", order.map { price[$0] })]145        )146        for (a, b) in zip(base.coefficients, permuted.coefficients) {147            expectClose(a.estimate, b.estimate, "permuted b[\(a.name)]")148            expectClose(a.standardError, b.standardError, "permuted se[\(a.name)]")149        }150    }151152    @Test("OLS invariance: scaling a regressor rescales its coefficient")153    func scalingInvariance() throws {154        let base = try ZQOLS.fit(y: y, predictors: [("price", price)])155        let scaled = try ZQOLS.fit(156            y: y, predictors: [("price", price.map { $0 * 100 })]157        )158        expectClose(159            scaled.coefficients[0].estimate * 100,160            base.coefficients[0].estimate,161            "scaled coefficient"162        )163        expectClose(scaled.rSquared, base.rSquared, "scaled R²")164    }165}166167@Suite("Distribution functions vs R")168struct DistributionTests {169    let fixtures = try! Fixtures()170171    @Test("Student t CDF")172    func studentT() {173        expectClose(174            ZQDistributions.studentTCDF(2.5, df: 10),175            fixtures["dist_pt_2p5_df10"], rtol: 1e-12, "pt(2.5, 10)"176        )177        expectClose(178            ZQDistributions.studentTCDF(-1.3, df: 3),179            fixtures["dist_pt_m1p3_df3"], rtol: 1e-12, "pt(-1.3, 3)"180        )181        expectClose(182            ZQDistributions.studentTCDF(0.05, df: 57),183            fixtures["dist_pt_0p05_df57"], rtol: 1e-12, "pt(0.05, 57)"184        )185    }186187    @Test("F CDF")188    func fisher() {189        expectClose(190            ZQDistributions.fCDF(3.7, df1: 2, df2: 30),191            fixtures["dist_pf_3p7_2_30"], rtol: 1e-12, "pf(3.7, 2, 30)"192        )193        expectClose(194            ZQDistributions.fCDF(0.5, df1: 5, df2: 100),195            fixtures["dist_pf_0p5_5_100"], rtol: 1e-12, "pf(0.5, 5, 100)"196        )197    }198199    @Test("t quantile")200    func tQuantile() {201        expectClose(202            ZQDistributions.studentTQuantile(0.975, df: 12),203            fixtures["dist_qt_0p975_df12"], rtol: 1e-10, "qt(0.975, 12)"204        )205        expectClose(206            ZQDistributions.studentTQuantile(0.995, df: 4),207            fixtures["dist_qt_0p995_df4"], rtol: 1e-10, "qt(0.995, 4)"208        )209    }210211    @Test("normal CDF")212    func normal() {213        expectClose(214            ZQDistributions.normalCDF(1.64),215            fixtures["dist_pnorm_1p64"], rtol: 1e-14, "pnorm(1.64)"216        )217    }218}219