*! rbfmvar_simulate — Monte Carlo simulation for RBFM-VAR
*! Reproduces the DGP from Section 5 of Chang (2000)
*! Version 1.0.0, February 2026
*! Author: Dr. Merwan Roudane
*! Email: merwanroudane920@gmail.com

capture program drop rbfmvar_simulate
program define rbfmvar_simulate, rclass
    version 14.0

    syntax , [                                    ///
        CASE(string)                              /// a b c
        NOBS(integer 150)                         /// sample size T
        REPS(integer 1000)                        /// Monte Carlo replications
        KERnel(string)                            /// kernel type
        SAVing(string)                            /// save dataset
        NOPRint                                   /// suppress table
        SEED(integer 12345)                       /// random seed
        ]

    *--------------------------------------------------------------------------
    * Defaults and validation
    *--------------------------------------------------------------------------
    if "`case'" == "" local case "a"
    local case = lower("`case'")
    if !inlist("`case'", "a", "b", "c") {
        di as error "case() must be {bf:a}, {bf:b}, or {bf:c}"
        exit 198
    }

    if "`kernel'" == "" local kernel "bartlett"
    local kernel = lower("`kernel'")

    if `nobs' < 30 {
        di as error "nobs() must be at least 30"
        exit 198
    }

    * DGP parameters from Chang (2000) Section 5
    if "`case'" == "a" {
        local rho1 = 1
        local rho2 = 0
        local case_label "Case A: (ρ₁,ρ₂) = (1,0) — Both I(2), no causality"
        local case_desc "Both Y₁ and Y₂ are I(2) with no cointegration"
    }
    else if "`case'" == "b" {
        local rho1 = 0.5
        local rho2 = 0
        local case_label "Case B: (ρ₁,ρ₂) = (0.5,0) — Y₁ I(1), Y₂ I(2), no causality"
        local case_desc "Y₁ is I(1), Y₂ is I(2), no Granger causality"
    }
    else {
        local rho1 = -0.3
        local rho2 = -0.15
        local case_label "Case C: (ρ₁,ρ₂) = (-0.3,-0.15) — Y₂ Granger-causes Y₁"
        local case_desc "Y₁ is I(1), Y₂ is I(2), Y₂ Granger-causes Y₁"
    }

    * Error covariance: Σ = [[1, 0.5], [0.5, 1.5]]
    local sig11 = 1
    local sig12 = 0.5
    local sig22 = 1.5

    *--------------------------------------------------------------------------
    * Display header
    *--------------------------------------------------------------------------
    if "`noprint'" == "" {
        di
        di as txt "{hline 78}"
        di as txt "{col 5}{bf:RBFM-VAR Monte Carlo Simulation}"
        di as txt "{col 5}{it:DGP from Chang (2000), Section 5, Eq. 24}"
        di as txt "{hline 78}"
        di
        di as txt "{col 5}`case_label'"
        di as txt "{col 5}`case_desc'"
        di
        di as txt "{col 5}Sample size (T):{col 30}" as res `nobs'
        di as txt "{col 5}Replications:{col 30}" as res `reps'
        di as txt "{col 5}Kernel:{col 30}" as res "`kernel'"
        di as txt "{col 5}Seed:{col 30}" as res `seed'
        di
        di as txt "{col 5}DGP: Δ²Y₁ₜ = ρ₁·ΔY₁,ₜ₋₁ + ρ₂·(Y₁,ₜ₋₁ - ΔY₂,ₜ₋₁) + ε₁ₜ"
        di as txt "{col 5}     Δ²Y₂ₜ = ε₂ₜ"
        di as txt "{col 5}     Σ = [[`sig11', `sig12'], [`sig12', `sig22']]"
        di
        di as txt "{col 5}Testing H₀: Y₂ does NOT Granger-cause Y₁"
        di
        di as txt "{hline 78}"
        di as txt "{col 5}Running simulations..."
    }

    *--------------------------------------------------------------------------
    * Load Mata engine
    *--------------------------------------------------------------------------
    capture findfile _rbfmvar_mata.ado
    if _rc {
        di as error "Required file _rbfmvar_mata.ado not found."
        exit 601
    }
    qui run `"`r(fn)'"'

    *--------------------------------------------------------------------------
    * Monte Carlo simulation via Mata
    *--------------------------------------------------------------------------
    set seed `seed'

    preserve
    clear

    * Run simulation
    mata: _rbfmvar_mc_run(`nobs', `reps', `rho1', `rho2', ///
        `sig11', `sig12', `sig22', "`kernel'")

    * Retrieve results
    local bias_ols_11   = r(bias_ols_11)
    local bias_ols_12   = r(bias_ols_12)
    local bias_ols_21   = r(bias_ols_21)
    local bias_ols_22   = r(bias_ols_22)
    local sd_ols_11     = r(sd_ols_11)
    local sd_ols_12     = r(sd_ols_12)
    local sd_ols_21     = r(sd_ols_21)
    local sd_ols_22     = r(sd_ols_22)

    local bias_rbfm_11  = r(bias_rbfm_11)
    local bias_rbfm_12  = r(bias_rbfm_12)
    local bias_rbfm_21  = r(bias_rbfm_21)
    local bias_rbfm_22  = r(bias_rbfm_22)
    local sd_rbfm_11    = r(sd_rbfm_11)
    local sd_rbfm_12    = r(sd_rbfm_12)
    local sd_rbfm_21    = r(sd_rbfm_21)
    local sd_rbfm_22    = r(sd_rbfm_22)

    local rej_ols_01    = r(rej_ols_01)
    local rej_ols_05    = r(rej_ols_05)
    local rej_ols_10    = r(rej_ols_10)

    local rej_rbfm_01   = r(rej_rbfm_01)
    local rej_rbfm_05   = r(rej_rbfm_05)
    local rej_rbfm_10   = r(rej_rbfm_10)

    restore

    *--------------------------------------------------------------------------
    * Display Table 1: Bias and Standard Deviations
    *--------------------------------------------------------------------------
    if "`noprint'" == "" {
        di
        di as txt "{hline 78}"
        di as txt "{col 5}{bf:Table 1: Bias and Standard Deviation of Π₁ Estimates}"
        di as txt "{col 5}{it:`case_label'}, T = `nobs'"
        di as txt "{hline 78}"
        di
        di as txt "{col 5}" _col(22) "OLS-VAR" _col(48) "RBFM-VAR"
        di as txt "{col 5}" _col(18) "Bias" _col(30) "S.D." _col(44) "Bias" _col(56) "S.D."
        di as txt "{col 5}{hline 60}"

        di as txt "{col 5}π₁₁" ///
           _col(15) as res %10.4f `bias_ols_11' ///
           _col(27) as res %10.4f `sd_ols_11' ///
           _col(41) as res %10.4f `bias_rbfm_11' ///
           _col(53) as res %10.4f `sd_rbfm_11'

        di as txt "{col 5}π₁₂" ///
           _col(15) as res %10.4f `bias_ols_12' ///
           _col(27) as res %10.4f `sd_ols_12' ///
           _col(41) as res %10.4f `bias_rbfm_12' ///
           _col(53) as res %10.4f `sd_rbfm_12'

        di as txt "{col 5}π₂₁" ///
           _col(15) as res %10.4f `bias_ols_21' ///
           _col(27) as res %10.4f `sd_ols_21' ///
           _col(41) as res %10.4f `bias_rbfm_21' ///
           _col(53) as res %10.4f `sd_rbfm_21'

        di as txt "{col 5}π₂₂" ///
           _col(15) as res %10.4f `bias_ols_22' ///
           _col(27) as res %10.4f `sd_ols_22' ///
           _col(41) as res %10.4f `bias_rbfm_22' ///
           _col(53) as res %10.4f `sd_rbfm_22'

        di as txt "{col 5}{hline 60}"

        *----------------------------------------------------------------------
        * Display Table 2: Size / Rejection Probabilities
        *----------------------------------------------------------------------
        di
        di as txt "{hline 78}"
        if "`case'" == "c" {
            di as txt "{col 5}{bf:Table 2: Rejection Probabilities (Power)}"
        }
        else {
            di as txt "{col 5}{bf:Table 2: Finite Sample Sizes (Rejection under H₀)}"
        }
        di as txt "{col 5}{it:`case_label'}, T = `nobs'"
        di as txt "{hline 78}"
        di
        di as txt "{col 5}" _col(25) "1%" _col(40) "5%" _col(55) "10%"
        di as txt "{col 5}{hline 60}"

        di as txt "{col 5}OLS-VAR (W_F)" ///
           _col(22) as res %10.4f `rej_ols_01' ///
           _col(37) as res %10.4f `rej_ols_05' ///
           _col(52) as res %10.4f `rej_ols_10'

        di as txt "{col 5}RBFM-VAR (W_F⁺)" ///
           _col(22) as res %10.4f `rej_rbfm_01' ///
           _col(37) as res %10.4f `rej_rbfm_05' ///
           _col(52) as res %10.4f `rej_rbfm_10'

        di as txt "{col 5}{hline 60}"
        di as txt "{col 5}Nominal size:" _col(22) "0.0100" _col(37) "0.0500" _col(52) "0.1000"
        di
        if "`case'" != "c" {
            di as txt "{col 5}Note: Closer to nominal = better size control."
            di as txt "{col 5}OLS-VAR typically shows severe size distortions."
        }
        else {
            di as txt "{col 5}Note: Higher = better power."
        }
        di as txt "{hline 78}"
    }

    *--------------------------------------------------------------------------
    * Return results
    *--------------------------------------------------------------------------
    return scalar reps       = `reps'
    return scalar nobs       = `nobs'

    return scalar bias_ols_11  = `bias_ols_11'
    return scalar bias_ols_12  = `bias_ols_12'
    return scalar bias_ols_21  = `bias_ols_21'
    return scalar bias_ols_22  = `bias_ols_22'
    return scalar sd_ols_11    = `sd_ols_11'
    return scalar sd_ols_12    = `sd_ols_12'
    return scalar sd_ols_21    = `sd_ols_21'
    return scalar sd_ols_22    = `sd_ols_22'

    return scalar bias_rbfm_11 = `bias_rbfm_11'
    return scalar bias_rbfm_12 = `bias_rbfm_12'
    return scalar bias_rbfm_21 = `bias_rbfm_21'
    return scalar bias_rbfm_22 = `bias_rbfm_22'
    return scalar sd_rbfm_11   = `sd_rbfm_11'
    return scalar sd_rbfm_12   = `sd_rbfm_12'
    return scalar sd_rbfm_21   = `sd_rbfm_21'
    return scalar sd_rbfm_22   = `sd_rbfm_22'

    return scalar rej_ols_01   = `rej_ols_01'
    return scalar rej_ols_05   = `rej_ols_05'
    return scalar rej_ols_10   = `rej_ols_10'
    return scalar rej_rbfm_01  = `rej_rbfm_01'
    return scalar rej_rbfm_05  = `rej_rbfm_05'
    return scalar rej_rbfm_10  = `rej_rbfm_10'

    return local  case         "`case'"
    return local  kernel       "`kernel'"
    return local  cmd          "rbfmvar_simulate"
end


// ========================================================================
// Mata: Monte Carlo simulation engine
// ========================================================================

capture mata: mata drop _rbfmvar_mc_run()
capture mata: mata drop _rbfmvar_mc_dgp()
capture mata: mata drop _rbfmvar_mc_wald_ols()

version 14.0
mata:
mata set matastrict on

// ========================================================================
// DGP: Generate one replication of the bivariate process (Eq. 24)
// ========================================================================

void _rbfmvar_mc_dgp(
    real scalar T,
    real scalar rho1,
    real scalar rho2,
    real matrix Sigma,
    real matrix Y_out)
{
    real matrix eps, L_chol
    real scalar t
    real colvector y1, y2, dy1, dy2, d2y1, d2y2

    // Cholesky for correlated errors
    L_chol = cholesky(Sigma)

    // Generate i.i.d. N(0, Sigma) errors
    eps = rnormal(T + 50, 2, 0, 1) * L_chol'

    // Initialize
    y1  = J(T + 50, 1, 0)
    y2  = J(T + 50, 1, 0)
    dy1 = J(T + 50, 1, 0)
    dy2 = J(T + 50, 1, 0)

    // Generate process
    for (t = 3; t <= T + 50; t++) {
        // Δ²Y₂ₜ = ε₂ₜ  =>  ΔY₂ₜ = ΔY₂,ₜ₋₁ + ε₂ₜ
        dy2[t] = dy2[t-1] + eps[t, 2]
        y2[t]  = y2[t-1] + dy2[t]

        // Δ²Y₁ₜ = ρ₁·ΔY₁,ₜ₋₁ + ρ₂·(Y₁,ₜ₋₁ - ΔY₂,ₜ₋₁) + ε₁ₜ
        d2y1 = rho1 * dy1[t-1] + rho2 * (y1[t-1] - dy2[t-1]) + eps[t, 1]
        dy1[t] = dy1[t-1] + d2y1
        y1[t]  = y1[t-1] + dy1[t]
    }

    // Trim burn-in
    Y_out = (y1[51::(T+50)], y2[51::(T+50)])
}


// ========================================================================
// OLS Wald statistic for Granger non-causality (for comparison)
// ========================================================================

void _rbfmvar_mc_wald_ols(
    real matrix F_ols,
    real matrix Sigma_e,
    real matrix X,
    real scalar T,
    real scalar n,
    real scalar ncols_z,
    real scalar wald_stat)
{
    real matrix R, XpXinv, V_mat, V_inv
    real colvector r_vec, vecF, Rf_r
    real scalar np_total, j

    np_total = cols(X)
    R = J(2, n * np_total, 0)
    r_vec = J(2, 1, 0)

    // Restriction: F[1, ncols_z + 2] = 0  (Π₁: effect of y2 on y1 equation)
    // vec(F') position for F[eq,reg] = (eq-1)*np_total + reg
    j = ncols_z + 2
    R[1, (1-1)*np_total + j] = 1

    // Restriction: F[1, ncols_z + n + 2] = 0  (Π₂: effect of y2 on y1 equation)
    j = ncols_z + n + 2
    R[2, (1-1)*np_total + j] = 1

    vecF = vec(F_ols')
    Rf_r = R * vecF - r_vec

    XpXinv = invsym(cross(X, X))
    V_mat = R * (Sigma_e # XpXinv) * R'
    V_inv = invsym(V_mat)

    wald_stat = Rf_r' * V_inv * Rf_r
}


// ========================================================================
// Main MC runner
// ========================================================================

void _rbfmvar_mc_run(
    real scalar T_obs,
    real scalar n_reps,
    real scalar rho1,
    real scalar rho2,
    real scalar sig11,
    real scalar sig12,
    real scalar sig22,
    string scalar kernel)
{
    real matrix Sigma, Y_sim, Y, Z, W, X
    real matrix F_ols, F_plus, E_hat, V_hat
    real matrix Omega_ev, Omega_vv, Delta_vdw, dW
    real matrix Gamma_plus, A_plus, Sigma_e
    real scalar bw, r, T_eff, n, ncols_z
    real scalar wald_ols, wald_rbfm

    // True Π₁ and Π₂ values from DGP (Eq. 24):
    //   Δ²y₁ₜ = ρ₁·Δy₁,ₜ₋₁ + ρ₂·(y₁,ₜ₋₁ - Δy₂,ₜ₋₁) + ε₁ₜ
    //         = ρ₁·Δy₁,ₜ₋₁ + (-ρ₂)·Δy₂,ₜ₋₁ + ρ₂·y₁,ₜ₋₁ + ε₁ₜ
    // So in ECM form y_t = Π₁·Δy_{t-1} + Π₂·y_{t-1} + ε_t:
    //   Π₁ = [[ρ₁, -ρ₂], [0, 0]]  and  Π₂ = [[ρ₂, 0], [0, 0]]
    real matrix Pi1_true
    Pi1_true = (rho1, -rho2 \ 0, 0)

    Sigma = (sig11, sig12 \ sig12, sig22)
    n = 2

    // Storage for estimates
    real matrix Pi1_ols_store, Pi1_rbfm_store
    real colvector wald_ols_store, wald_rbfm_store
    Pi1_ols_store  = J(n_reps, 4, 0)
    Pi1_rbfm_store = J(n_reps, 4, 0)
    wald_ols_store  = J(n_reps, 1, 0)
    wald_rbfm_store = J(n_reps, 1, 0)

    real scalar p_lags, valid_reps
    p_lags = 1
    valid_reps = 0

    for (r = 1; r <= n_reps; r++) {
        // Generate DGP
        _rbfmvar_mc_dgp(T_obs, rho1, rho2, Sigma, Y_sim)

        // ECM reparameterize
        _rbfm_ecm_reparameterize(Y_sim, p_lags, Y, Z, W)
        T_eff = rows(Y)
        X = (Z, W)
        ncols_z = cols(Z)

        // OLS estimation
        _rbfm_ols_var(Y, X, Z, W, F_ols, E_hat)

        // Construct v_hat
        _rbfm_construct_vhat(Y_sim, p_lags, F_ols, X, V_hat)

        // Bandwidth
        bw = _rbfm_andrews_bandwidth(V_hat, kernel)

        // LRV
        _rbfm_kernel_lrv(E_hat, V_hat, kernel, bw, Omega_ev, Omega_vv)

        // One-sided LRV
        dW = W[2::T_eff, .] - W[1::(T_eff-1), .]
        dW = (J(1, cols(W), 0) \ dW)
        _rbfm_kernel_olrv(V_hat, dW, kernel, bw, Delta_vdw)

        // RBFM correction
        _rbfm_correct(Y, Z, W, X, E_hat, V_hat, Omega_ev, Omega_vv, Delta_vdw, T_eff, F_plus, Gamma_plus, A_plus)

        // Residual covariance
        Sigma_e = (E_hat' * E_hat) / T_eff

        // Extract Π₁ estimates (cols ncols_z+1 to ncols_z+n)
        real matrix Pi1_ols_r, Pi1_rbfm_r
        Pi1_ols_r  = F_ols[., (ncols_z+1)::(ncols_z+n)]
        Pi1_rbfm_r = F_plus[., (ncols_z+1)::(ncols_z+n)]

        // Store: vec(Π₁) = (π₁₁, π₂₁, π₁₂, π₂₂)
        Pi1_ols_store[r, .]  = (Pi1_ols_r[1,1], Pi1_ols_r[1,2], Pi1_ols_r[2,1], Pi1_ols_r[2,2])
        Pi1_rbfm_store[r, .] = (Pi1_rbfm_r[1,1], Pi1_rbfm_r[1,2], Pi1_rbfm_r[2,1], Pi1_rbfm_r[2,2])

        // Wald tests
        _rbfmvar_mc_wald_ols(F_ols, Sigma_e, X, T_eff, n, ncols_z, wald_ols)
        wald_ols_store[r] = wald_ols

        // RBFM Wald — reuse the same R matrix structure
        real scalar wald_tmp, wald_p_tmp, wald_df_tmp
        real matrix R_gc
        real colvector r_gc
        real scalar j_col, np_total

        np_total = cols(X)
        R_gc = J(2, n * np_total, 0)
        r_gc = J(2, 1, 0)

        j_col = ncols_z + 2
        R_gc[1, (1-1)*np_total + j_col] = 1
        j_col = ncols_z + n + 2
        R_gc[2, (1-1)*np_total + j_col] = 1

        _rbfm_wald_test(R_gc, r_gc, F_plus, Sigma_e, X, T_eff, wald_tmp, wald_p_tmp, wald_df_tmp)
        wald_rbfm_store[r] = wald_tmp

        valid_reps = valid_reps + 1
    }

    // Compute bias and s.d.
    real rowvector true_vals, mean_ols, mean_rbfm, sd_ols, sd_rbfm, bias_ols, bias_rbfm
    true_vals = (Pi1_true[1,1], Pi1_true[1,2], Pi1_true[2,1], Pi1_true[2,2])

    mean_ols  = mean(Pi1_ols_store)
    mean_rbfm = mean(Pi1_rbfm_store)

    bias_ols  = mean_ols - true_vals
    bias_rbfm = mean_rbfm - true_vals

    sd_ols  = J(1, 4, 0)
    sd_rbfm = J(1, 4, 0)
    for (r = 1; r <= 4; r++) {
        sd_ols[r]  = sqrt(variance(Pi1_ols_store[., r]))
        sd_rbfm[r] = sqrt(variance(Pi1_rbfm_store[., r]))
    }

    // Rejection rates (χ²_2 critical values)
    real scalar cv01, cv05, cv10
    cv01 = invchi2(2, 0.99)  // 9.210
    cv05 = invchi2(2, 0.95)  // 5.991
    cv10 = invchi2(2, 0.90)  // 4.605

    real scalar rej_ols_01, rej_ols_05, rej_ols_10
    real scalar rej_rbfm_01, rej_rbfm_05, rej_rbfm_10

    rej_ols_01  = mean(wald_ols_store :> cv01)
    rej_ols_05  = mean(wald_ols_store :> cv05)
    rej_ols_10  = mean(wald_ols_store :> cv10)

    rej_rbfm_01 = mean(wald_rbfm_store :> cv01)
    rej_rbfm_05 = mean(wald_rbfm_store :> cv05)
    rej_rbfm_10 = mean(wald_rbfm_store :> cv10)

    // Post to Stata
    st_numscalar("r(bias_ols_11)", bias_ols[1])
    st_numscalar("r(bias_ols_12)", bias_ols[2])
    st_numscalar("r(bias_ols_21)", bias_ols[3])
    st_numscalar("r(bias_ols_22)", bias_ols[4])
    st_numscalar("r(sd_ols_11)", sd_ols[1])
    st_numscalar("r(sd_ols_12)", sd_ols[2])
    st_numscalar("r(sd_ols_21)", sd_ols[3])
    st_numscalar("r(sd_ols_22)", sd_ols[4])

    st_numscalar("r(bias_rbfm_11)", bias_rbfm[1])
    st_numscalar("r(bias_rbfm_12)", bias_rbfm[2])
    st_numscalar("r(bias_rbfm_21)", bias_rbfm[3])
    st_numscalar("r(bias_rbfm_22)", bias_rbfm[4])
    st_numscalar("r(sd_rbfm_11)", sd_rbfm[1])
    st_numscalar("r(sd_rbfm_12)", sd_rbfm[2])
    st_numscalar("r(sd_rbfm_21)", sd_rbfm[3])
    st_numscalar("r(sd_rbfm_22)", sd_rbfm[4])

    st_numscalar("r(rej_ols_01)", rej_ols_01)
    st_numscalar("r(rej_ols_05)", rej_ols_05)
    st_numscalar("r(rej_ols_10)", rej_ols_10)
    st_numscalar("r(rej_rbfm_01)", rej_rbfm_01)
    st_numscalar("r(rej_rbfm_05)", rej_rbfm_05)
    st_numscalar("r(rej_rbfm_10)", rej_rbfm_10)
}


end