greeners 1.4.2

//! Numerical accuracy tests against statsmodels / linearmodels reference values.
//! Reference data generated by tests/generate_reference_data.py

use greeners::{
    CovarianceType, DataFrame, FixedEffects, Formula, InferenceType, Logit, Probit, QuantileReg,
    RandomEffects, FGLS, IV, OLS,
};
use ndarray::{Array1, Array2};
use std::collections::HashMap;
use std::fs;

/// Tolerance for coefficient / SE comparisons (4 decimal places)
const TOL: f64 = 1e-4;
/// Tolerance for values with minor algorithmic differences (3 decimal places)
const TOL_MID: f64 = 1e-3;
/// Looser tolerance for values sensitive to algorithmic differences
const TOL_LOOSE: f64 = 1e-2;

fn load_reference() -> HashMap<String, f64> {
    let content = fs::read_to_string("tests/reference_values.csv")
        .expect("tests/reference_values.csv not found — run generate_reference_data.py first");
    let mut map = HashMap::new();
    for line in content.lines().skip(1) {
        let mut parts = line.splitn(2, ',');
        if let (Some(key), Some(val)) = (parts.next(), parts.next()) {
            if let Ok(v) = val.parse::<f64>() {
                map.insert(key.to_string(), v);
            }
        }
    }
    map
}

fn read_vec(ref_data: &HashMap<String, f64>, prefix: &str, n: usize) -> Vec<f64> {
    (0..n)
        .map(|i| *ref_data.get(&format!("{}.{}", prefix, i)).unwrap())
        .collect()
}

fn assert_close(actual: f64, expected: f64, tol: f64, label: &str) {
    assert!(
        (actual - expected).abs() < tol,
        "{}: actual={:.10}, expected={:.10}, diff={:.2e}",
        label,
        actual,
        expected,
        (actual - expected).abs()
    );
}

/// Relative tolerance: |actual - expected| / max(|expected|, 1e-10) < rel_tol
fn assert_close_rel(actual: f64, expected: f64, rel_tol: f64, label: &str) {
    let denom = expected.abs().max(1e-10);
    let rel_err = (actual - expected).abs() / denom;
    assert!(
        rel_err < rel_tol,
        "{}: actual={:.10}, expected={:.10}, rel_err={:.2e}",
        label,
        actual,
        expected,
        rel_err
    );
}

// ============================================================================
// OLS
// ============================================================================

fn build_ols_data(ref_data: &HashMap<String, f64>) -> (DataFrame, Formula) {
    let n = 50;
    let y = read_vec(ref_data, "ols_data.y", n);
    let x1 = read_vec(ref_data, "ols_data.x1", n);
    let x2 = read_vec(ref_data, "ols_data.x2", n);
    let x3 = read_vec(ref_data, "ols_data.x3", n);

    let mut data = HashMap::new();
    data.insert("y".to_string(), Array1::from(y));
    data.insert("x1".to_string(), Array1::from(x1));
    data.insert("x2".to_string(), Array1::from(x2));
    data.insert("x3".to_string(), Array1::from(x3));

    let df = DataFrame::new(data).unwrap();
    let formula = Formula::parse("y ~ x1 + x2 + x3").unwrap();
    (df, formula)
}

#[test]
fn test_ols_nonrobust_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);
    let result = OLS::from_formula(&formula, &df, CovarianceType::NonRobust).unwrap();

    // Variable order: const, x1, x2, x3
    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_nonrobust.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL, &key);

        let key = format!("ols_nonrobust.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL, &key);

        let key = format!("ols_nonrobust.t.{}", name);
        assert_close(result.t_values[i], ref_data[&key], TOL, &key);

        let key = format!("ols_nonrobust.p.{}", name);
        assert_close(result.p_values[i], ref_data[&key], TOL, &key);
    }

    assert_close(
        result.r_squared,
        ref_data["ols_nonrobust.r_squared"],
        TOL,
        "r_squared",
    );
    assert_close(
        result.adj_r_squared,
        ref_data["ols_nonrobust.adj_r_squared"],
        TOL,
        "adj_r_squared",
    );
    assert_close(
        result.f_statistic,
        ref_data["ols_nonrobust.f_statistic"],
        TOL,
        "f_statistic",
    );
    assert_close(
        result.log_likelihood,
        ref_data["ols_nonrobust.log_likelihood"],
        TOL,
        "log_likelihood",
    );
    assert_close(result.aic, ref_data["ols_nonrobust.aic"], TOL, "aic");
    assert_close(result.bic, ref_data["ols_nonrobust.bic"], TOL, "bic");
    assert_close(result.sigma, ref_data["ols_nonrobust.sigma"], TOL, "sigma");
    assert_eq!(result.n_obs, 50);
    assert_eq!(result.df_resid, 46);
}

#[test]
fn test_ols_hc1_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);
    let result = OLS::from_formula(&formula, &df, CovarianceType::HC1).unwrap();

    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_hc1.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL, &key);

        let key = format!("ols_hc1.t.{}", name);
        assert_close(result.t_values[i], ref_data[&key], TOL, &key);
    }
}

#[test]
fn test_ols_hc3_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);
    let result = OLS::from_formula(&formula, &df, CovarianceType::HC3).unwrap();

    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_hc3.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL, &key);

        let key = format!("ols_hc3.t.{}", name);
        assert_close(result.t_values[i], ref_data[&key], TOL, &key);
    }
}

#[test]
fn test_ols_normal_inference_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);
    let result = OLS::from_formula(&formula, &df, CovarianceType::NonRobust)
        .unwrap()
        .with_inference(InferenceType::Normal)
        .unwrap();

    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_normal.p.{}", name);
        assert_close(result.p_values[i], ref_data[&key], TOL, &key);
    }
}

// ============================================================================
// IV / 2SLS
// ============================================================================

#[test]
fn test_iv_2sls_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "iv_data.y", n));
    let x_endog = read_vec(&ref_data, "iv_data.x_endog", n);
    let z1 = read_vec(&ref_data, "iv_data.z1", n);
    let z2 = read_vec(&ref_data, "iv_data.z2", n);

    // X = [const, x_endog]
    let mut x_flat = Vec::with_capacity(n * 2);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x_endog[i]);
    }
    let x = Array2::from_shape_vec((n, 2), x_flat).unwrap();

    // Z = [const, z1, z2]
    let mut z_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        z_flat.push(1.0);
        z_flat.push(z1[i]);
        z_flat.push(z2[i]);
    }
    let z = Array2::from_shape_vec((n, 3), z_flat).unwrap();

    let result = IV::fit(&y, &x, &z, CovarianceType::NonRobust).unwrap();

    let names = ["const", "x_endog"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("iv_2sls.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL, &key);

        let key = format!("iv_2sls.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL, &key);
    }
}

// ============================================================================
// Logit
// ============================================================================

fn build_binary_data(ref_data: &HashMap<String, f64>) -> (DataFrame, Formula) {
    let n = 200;
    let y = read_vec(ref_data, "binary_data.y", n);
    let x1 = read_vec(ref_data, "binary_data.x1", n);
    let x2 = read_vec(ref_data, "binary_data.x2", n);

    let mut data = HashMap::new();
    data.insert("y".to_string(), Array1::from(y));
    data.insert("x1".to_string(), Array1::from(x1));
    data.insert("x2".to_string(), Array1::from(x2));

    let df = DataFrame::new(data).unwrap();
    let formula = Formula::parse("y ~ x1 + x2").unwrap();
    (df, formula)
}

#[test]
fn test_logit_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_binary_data(&ref_data);
    let result = Logit::from_formula(&formula, &df).unwrap();

    let names = ["const", "x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("logit.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL, &key);

        let key = format!("logit.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL, &key);

        let key = format!("logit.z.{}", name);
        assert_close(result.z_values[i], ref_data[&key], TOL, &key);

        let key = format!("logit.p.{}", name);
        assert_close(result.p_values[i], ref_data[&key], TOL, &key);
    }

    assert_close(
        result.log_likelihood,
        ref_data["logit.log_likelihood"],
        TOL,
        "logit.log_likelihood",
    );
    assert_close(
        result.pseudo_r2,
        ref_data["logit.pseudo_r2"],
        TOL,
        "logit.pseudo_r2",
    );
}

#[test]
fn test_logit_ame_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 200;
    let y = Array1::from(read_vec(&ref_data, "binary_data.y", n));
    let x1 = read_vec(&ref_data, "binary_data.x1", n);
    let x2 = read_vec(&ref_data, "binary_data.x2", n);

    // Build X with intercept column
    let mut x_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 3), x_flat).unwrap();

    let result = Logit::fit(&y, &x).unwrap();
    let ame = result.average_marginal_effects(&x).unwrap();

    // AME includes all columns (const, x1, x2); statsmodels only reports non-const
    // In Greeners, AME[0] is for const, AME[1] for x1, AME[2] for x2
    let names = ["x1", "x2"];
    for (j, name) in names.iter().enumerate() {
        let key = format!("logit.ame.{}", name);
        assert_close(ame[j + 1], ref_data[&key], TOL_LOOSE, &key);
    }
}

// ============================================================================
// Probit
// ============================================================================

#[test]
fn test_probit_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_binary_data(&ref_data);
    let result = Probit::from_formula(&formula, &df).unwrap();

    let names = ["const", "x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("probit.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL_MID, &key);

        let key = format!("probit.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL_MID, &key);

        let key = format!("probit.z.{}", name);
        assert_close_rel(result.z_values[i], ref_data[&key], 0.01, &key);

        // p-values near zero amplify tiny z-value differences
        let key = format!("probit.p.{}", name);
        assert_close(result.p_values[i], ref_data[&key], TOL_MID, &key);
    }

    assert_close(
        result.log_likelihood,
        ref_data["probit.log_likelihood"],
        TOL_MID,
        "probit.log_likelihood",
    );
    assert_close(
        result.pseudo_r2,
        ref_data["probit.pseudo_r2"],
        TOL_MID,
        "probit.pseudo_r2",
    );
}

// ============================================================================
// Panel: Fixed Effects
// ============================================================================

#[test]
fn test_panel_fe_vs_linearmodels() {
    let ref_data = load_reference();
    let n = 100; // 20 entities * 5 periods
    let y = Array1::from(read_vec(&ref_data, "panel_data.y", n));
    let x1 = read_vec(&ref_data, "panel_data.x1", n);
    let x2 = read_vec(&ref_data, "panel_data.x2", n);
    let entities: Vec<i64> = read_vec(&ref_data, "panel_data.entity", n)
        .iter()
        .map(|&v| v as i64)
        .collect();

    let mut x_flat = Vec::with_capacity(n * 2);
    for i in 0..n {
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 2), x_flat).unwrap();

    let result = FixedEffects::fit(&y, &x, &entities).unwrap();

    let names = ["x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("panel_fe.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL, &key);

        // SE/t differ due to df correction: Greeners uses N*T-N-k, linearmodels may differ
        let key = format!("panel_fe.se.{}", name);
        assert_close_rel(result.std_errors[i], ref_data[&key], 0.01, &key);

        let key = format!("panel_fe.t.{}", name);
        assert_close_rel(result.t_values[i], ref_data[&key], 0.01, &key);
    }

    assert_close(
        result.r_squared,
        ref_data["panel_fe.r_squared"],
        TOL_MID,
        "panel_fe.r_squared",
    );
    assert_eq!(result.n_obs, 100);
    assert_eq!(result.n_entities, 20);
}

// ============================================================================
// Panel: Random Effects
// ============================================================================

#[test]
fn test_panel_re_vs_linearmodels() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "panel_data.y", n));
    let x1 = read_vec(&ref_data, "panel_data.x1", n);
    let x2 = read_vec(&ref_data, "panel_data.x2", n);
    let entities = Array1::from(
        read_vec(&ref_data, "panel_data.entity", n)
            .iter()
            .map(|&v| v as i64)
            .collect::<Vec<_>>(),
    );

    // RE includes intercept
    let mut x_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 3), x_flat).unwrap();

    let result = RandomEffects::fit(&y, &x, &entities).unwrap();

    let names = ["const", "x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("panel_re.params.{}", name);
        // RE implementations can vary slightly due to Swamy-Arora vs other variance estimators
        assert_close(result.params[i], ref_data[&key], TOL_LOOSE, &key);

        let key = format!("panel_re.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL_LOOSE, &key);
    }
}

// ============================================================================
// Quantile Regression
// ============================================================================

#[test]
fn test_quantile_median_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "quantile_data.y", n));
    let x1 = read_vec(&ref_data, "quantile_data.x1", n);
    let x2 = read_vec(&ref_data, "quantile_data.x2", n);

    let mut x_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 3), x_flat).unwrap();

    // Median regression (tau=0.50) — point estimates should match closely
    let result = QuantileReg::fit(&y, &x, 0.50, 200).unwrap();

    let names = ["const", "x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("quantile_0.50.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL_LOOSE, &key);
    }
}

#[test]
fn test_quantile_025_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "quantile_data.y", n));
    let x1 = read_vec(&ref_data, "quantile_data.x1", n);
    let x2 = read_vec(&ref_data, "quantile_data.x2", n);

    let mut x_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 3), x_flat).unwrap();

    let result = QuantileReg::fit(&y, &x, 0.25, 200).unwrap();

    let names = ["const", "x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("quantile_0.25.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL_LOOSE, &key);
    }
}

#[test]
fn test_quantile_075_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "quantile_data.y", n));
    let x1 = read_vec(&ref_data, "quantile_data.x1", n);
    let x2 = read_vec(&ref_data, "quantile_data.x2", n);

    let mut x_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 3), x_flat).unwrap();

    let result = QuantileReg::fit(&y, &x, 0.75, 200).unwrap();

    let names = ["const", "x1", "x2"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("quantile_0.75.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL_LOOSE, &key);
    }
}

// ============================================================================
// WLS
// ============================================================================

#[test]
fn test_wls_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 50;
    let y = Array1::from(read_vec(&ref_data, "wls_data.y", n));
    let x1 = read_vec(&ref_data, "wls_data.x1", n);
    let weights = Array1::from(read_vec(&ref_data, "wls_data.weights", n));

    let mut x_flat = Vec::with_capacity(n * 2);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
    }
    let x = Array2::from_shape_vec((n, 2), x_flat).unwrap();

    let result = FGLS::wls(&y, &x, &weights).unwrap();

    let names = ["const", "x1"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("wls.params.{}", name);
        assert_close(result.params[i], ref_data[&key], TOL, &key);

        let key = format!("wls.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL_MID, &key);

        let key = format!("wls.t.{}", name);
        assert_close(result.t_values[i], ref_data[&key], TOL_MID, &key);
    }

    assert_close(
        result.r_squared,
        ref_data["wls.r_squared"],
        TOL_MID,
        "wls.r_squared",
    );
}

// ============================================================================
// OLS HC2
// ============================================================================

#[test]
fn test_ols_hc2_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);
    let result = OLS::from_formula(&formula, &df, CovarianceType::HC2).unwrap();

    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_hc2.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL, &key);

        let key = format!("ols_hc2.t.{}", name);
        assert_close(result.t_values[i], ref_data[&key], TOL, &key);
    }
}

// ============================================================================
// OLS Newey-West
// ============================================================================

#[test]
fn test_ols_neweywest_vs_statsmodels() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);
    let result = OLS::from_formula(&formula, &df, CovarianceType::NeweyWest(3)).unwrap();

    // Newey-West can differ slightly due to kernel weighting implementation details
    let tol_nw = 0.05;
    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_nw.se.{}", name);
        assert_close_rel(result.std_errors[i], ref_data[&key], tol_nw, &key);

        let key = format!("ols_nw.t.{}", name);
        assert_close_rel(result.t_values[i], ref_data[&key], tol_nw, &key);
    }
}

// ============================================================================
// OLS Clustered
// ============================================================================

#[test]
fn test_ols_clustered_vs_statsmodels() {
    let ref_data = load_reference();
    let n = 50;
    let (df, formula) = build_ols_data(&ref_data);

    // Read cluster IDs from reference data
    let cluster_ids: Vec<usize> = (0..n)
        .map(|i| *ref_data.get(&format!("ols_data.cluster.{}", i)).unwrap() as usize)
        .collect();

    let result = OLS::from_formula(&formula, &df, CovarianceType::Clustered(cluster_ids)).unwrap();

    let names = ["const", "x1", "x2", "x3"];
    for (i, name) in names.iter().enumerate() {
        let key = format!("ols_clustered.se.{}", name);
        assert_close(result.std_errors[i], ref_data[&key], TOL_LOOSE, &key);

        let key = format!("ols_clustered.t.{}", name);
        assert_close(result.t_values[i], ref_data[&key], TOL_LOOSE, &key);
    }
}

// ============================================================================
// with_inference tests
// ============================================================================

#[test]
fn test_ols_with_inference_normal() {
    let ref_data = load_reference();
    let (df, formula) = build_ols_data(&ref_data);

    let result_t = OLS::from_formula(&formula, &df, CovarianceType::NonRobust).unwrap();
    let result_z = result_t.with_inference(InferenceType::Normal).unwrap();

    // p-values should differ between StudentT and Normal inference
    for i in 0..4 {
        // With Normal inference, p-values are typically slightly smaller for significant coefficients
        // They should at least be different
        assert!(
            (result_z.p_values[i]
                - ref_data[&format!("ols_nonrobust.p.{}", ["const", "x1", "x2", "x3"][i])])
                .abs()
                > 1e-10
                || result_z.p_values[i] < 1e-10,
            "Normal inference p-values should differ from StudentT"
        );
    }
}

#[test]
fn test_iv_with_inference_normal() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "iv_data.y", n));
    let x_endog = read_vec(&ref_data, "iv_data.x_endog", n);
    let z1 = read_vec(&ref_data, "iv_data.z1", n);
    let z2 = read_vec(&ref_data, "iv_data.z2", n);

    let mut x_flat = Vec::with_capacity(n * 2);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x_endog[i]);
    }
    let x = Array2::from_shape_vec((n, 2), x_flat).unwrap();

    let mut z_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        z_flat.push(1.0);
        z_flat.push(z1[i]);
        z_flat.push(z2[i]);
    }
    let z = Array2::from_shape_vec((n, 3), z_flat).unwrap();

    let result_t = IV::fit(&y, &x, &z, CovarianceType::NonRobust).unwrap();
    let result_z = result_t.with_inference(InferenceType::Normal).unwrap();

    // Params should be identical
    for i in 0..2 {
        assert_close(
            result_z.params[i],
            result_z.params[i],
            TOL,
            "iv params unchanged",
        );
    }
    // p-values should differ
    assert!(
        result_z.inference_type == InferenceType::Normal,
        "Inference type should be Normal"
    );
}

#[test]
fn test_panel_fe_with_inference_normal() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "panel_data.y", n));
    let x1 = read_vec(&ref_data, "panel_data.x1", n);
    let x2 = read_vec(&ref_data, "panel_data.x2", n);
    let entities: Vec<i64> = read_vec(&ref_data, "panel_data.entity", n)
        .iter()
        .map(|&v| v as i64)
        .collect();

    let mut x_flat = Vec::with_capacity(n * 2);
    for i in 0..n {
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 2), x_flat).unwrap();

    let result_t = FixedEffects::fit(&y, &x, &entities).unwrap();
    let p_values_t = result_t.p_values.clone();
    let result_z = result_t.with_inference(InferenceType::Normal).unwrap();

    // p-values should change
    assert!(
        (result_z.p_values[0] - p_values_t[0]).abs() > 1e-12 || p_values_t[0] < 1e-12,
        "FE Normal inference p-values should differ from StudentT"
    );
}

#[test]
fn test_panel_re_with_inference_normal() {
    let ref_data = load_reference();
    let n = 100;
    let y = Array1::from(read_vec(&ref_data, "panel_data.y", n));
    let x1 = read_vec(&ref_data, "panel_data.x1", n);
    let x2 = read_vec(&ref_data, "panel_data.x2", n);
    let entities = Array1::from(
        read_vec(&ref_data, "panel_data.entity", n)
            .iter()
            .map(|&v| v as i64)
            .collect::<Vec<_>>(),
    );

    let mut x_flat = Vec::with_capacity(n * 3);
    for i in 0..n {
        x_flat.push(1.0);
        x_flat.push(x1[i]);
        x_flat.push(x2[i]);
    }
    let x = Array2::from_shape_vec((n, 3), x_flat).unwrap();

    let result_t = RandomEffects::fit(&y, &x, &entities).unwrap();
    let p_values_t = result_t.p_values.clone();
    let result_z = result_t.with_inference(InferenceType::Normal).unwrap();

    // p-values should change
    assert!(
        (result_z.p_values[0] - p_values_t[0]).abs() > 1e-12 || p_values_t[0] < 1e-12,
        "RE Normal inference p-values should differ from StudentT"
    );
}