*! _rardl_simulate — Monte Carlo simulation for ARDL bounds test critical values
*! Version 1.1.0
*! Implements Pesaran, Shin & Smith (2001) DGP for critical value generation
*! Used for both Rolling-Window ARDL and Recursive ARDL bounds tests

capture program drop _rardl_simulate
program define _rardl_simulate, rclass
    version 17

    syntax , tobs(numlist min=1 integer >10) ///
             nregs(integer) ///
           [ nsim(integer 50000) ///
             maxlag(integer 4) ///
             ic(string) ///
             seed(integer -1) ///
             NOIsily ///
             notable ]

    // =========================================================================
    // 1. VALIDATE INPUTS
    // =========================================================================
    if `nregs' < 1 {
        di as err "nregs() must be at least 1"
        exit 198
    }
    if `nsim' < 100 {
        di as err "nsim() must be at least 100"
        exit 198
    }
    if "`ic'" == "" local ic "bic"
    local ic = lower("`ic'")
    if !inlist("`ic'", "aic", "bic", "hqic") {
        di as err "ic() must be aic, bic, or hqic"
        exit 198
    }

    if `seed' > 0 {
        set seed `seed'
    }

    local ntobs : word count `tobs'

    // =========================================================================
    // 2. HEADER
    // =========================================================================
    if "`notable'" == "" {
        di as txt ""
        di as txt "{hline 78}"
        di as txt "{bf:Monte Carlo Simulation for ARDL Bounds Test Critical Values}"
        di as txt "{hline 78}"
        di as txt _col(5) "Replications     : " as res "`nsim'"
        di as txt _col(5) "No. of regressors: " as res "`nregs'"
        di as txt _col(5) "Max lag length   : " as res "`maxlag'"
        di as txt _col(5) "IC for lag select: " as res upper("`ic'")
        di as txt _col(5) "Sample sizes     : " as res "`tobs'"
        di as txt "{hline 78}"
        di as txt ""
    }

    // =========================================================================
    // 3. SIMULATION LOOP — For each sample size
    // =========================================================================
    local burnin = 300

    // Result matrices
    forvalues m = 1/5 {
        tempname cv_m`m'
        mat `cv_m`m'' = J(`ntobs', 7, .)
    }

    local tidx = 0
    foreach T of numlist `tobs' {
        local ++tidx
        local Tgen = `T' + `burnin'

        if "`notable'" == "" {
            di as txt _col(3) "Simulating for T = `T' ..." _continue
        }

        // Temp storage for F-stats
        forvalues m = 1/5 {
            tempname fstats_i0_m`m' fstats_i1_m`m'
            mat `fstats_i0_m`m'' = J(`nsim', 1, .)
            mat `fstats_i1_m`m'' = J(`nsim', 1, .)
        }

        // Create dataset once, reuse via replace
        preserve
        qui clear
        local nobs = max(`Tgen', `nsim')
        qui set obs `nobs'
        qui gen t = _n
        qui tsset t

        // Pre-create all variables
        qui gen double _sim_y = .
        qui gen double _sim_e = .
        forvalues j = 1/`nregs' {
            qui gen double _sim_x`j'_i0 = .
            qui gen double _sim_x`j'_i1 = .
            qui gen double _sim_ex`j' = .
        }

        forvalues sim = 1/`nsim' {

            // Generate y as random walk
            qui replace _sim_e = rnormal()
            qui replace _sim_y = 0 in 1
            qui replace _sim_y = _sim_y[_n-1] + _sim_e in 2/`Tgen'

            // Generate X variables
            forvalues j = 1/`nregs' {
                qui replace _sim_ex`j' = rnormal()
                // I(0): stationary AR(1)
                qui replace _sim_x`j'_i0 = 0 in 1
                qui replace _sim_x`j'_i0 = 0.5 * _sim_x`j'_i0[_n-1] + _sim_ex`j' in 2/`Tgen'
                // I(1): random walk
                qui replace _sim_x`j'_i1 = 0 in 1
                qui replace _sim_x`j'_i1 = _sim_x`j'_i1[_n-1] + _sim_ex`j' in 2/`Tgen'
            }

            // Use fixed lag = 1 for simulation speed
            local plag = 1

            // Compute F-stat for each of the 5 models, both I(0) and I(1)
            foreach iorder in i0 i1 {

                // Build lag terms
                local lagdy "L1.D._sim_y"

                local lxterms ""
                local dxterms ""
                forvalues j = 1/`nregs' {
                    local lxterms "`lxterms' L._sim_x`j'_`iorder'"
                    local dxterms "`dxterms' D._sim_x`j'_`iorder'"
                    local dxterms "`dxterms' L1.D._sim_x`j'_`iorder'"
                }

                // --- Model 1: No intercept, No trend ---
                capture qui reg D._sim_y L._sim_y `lxterms' `lagdy' `dxterms' ///
                    in `=`burnin'+1'/`Tgen', noconstant
                if _rc == 0 {
                    capture qui testparm L._sim_y `lxterms'
                    if _rc == 0 {
                        mat `fstats_`iorder'_m1'[`sim', 1] = r(F)
                    }
                }

                // --- Model 2: Restricted intercept, No trend ---
                capture qui reg D._sim_y L._sim_y `lxterms' `lagdy' `dxterms' ///
                    in `=`burnin'+1'/`Tgen'
                if _rc == 0 {
                    capture qui testparm L._sim_y `lxterms'
                    if _rc == 0 {
                        mat `fstats_`iorder'_m2'[`sim', 1] = r(F)
                    }
                }

                // --- Model 3: Unrestricted intercept, No trend ---
                if _rc == 0 {
                    capture qui testparm L._sim_y `lxterms'
                    if _rc == 0 {
                        mat `fstats_`iorder'_m3'[`sim', 1] = r(F)
                    }
                }

                // --- Model 4: Unrestricted intercept, Restricted trend ---
                capture qui reg D._sim_y L._sim_y `lxterms' t `lagdy' `dxterms' ///
                    in `=`burnin'+1'/`Tgen'
                if _rc == 0 {
                    capture qui testparm L._sim_y `lxterms' t
                    if _rc == 0 {
                        mat `fstats_`iorder'_m4'[`sim', 1] = r(F)
                    }
                }

                // --- Model 5: Unrestricted intercept, Unrestricted trend ---
                capture qui reg D._sim_y t L._sim_y `lxterms' `lagdy' `dxterms' ///
                    in `=`burnin'+1'/`Tgen'
                if _rc == 0 {
                    capture qui testparm L._sim_y `lxterms'
                    if _rc == 0 {
                        mat `fstats_`iorder'_m5'[`sim', 1] = r(F)
                    }
                }
            }

            // Progress
            if mod(`sim', 1000) == 0 & "`notable'" == "" {
                di as txt "." _continue
            }
        }

        if "`notable'" == "" {
            di as txt " done"
        }

        // =====================================================================
        // Extract percentiles for critical values
        // =====================================================================
        forvalues m = 1/5 {
            mat `cv_m`m''[`tidx', 1] = `T'

            foreach iorder in i0 i1 {
                capture drop _fsort
                qui gen double _fsort = .
                forvalues s = 1/`nsim' {
                    qui replace _fsort = el(`fstats_`iorder'_m`m'', `s', 1) in `s'
                }

                qui _pctile _fsort, percentiles(90 95 99)
                local p10 = r(r1)
                local p05 = r(r2)
                local p01 = r(r3)

                if "`iorder'" == "i0" {
                    mat `cv_m`m''[`tidx', 2] = `p01'
                    mat `cv_m`m''[`tidx', 4] = `p05'
                    mat `cv_m`m''[`tidx', 6] = `p10'
                }
                else {
                    mat `cv_m`m''[`tidx', 3] = `p01'
                    mat `cv_m`m''[`tidx', 5] = `p05'
                    mat `cv_m`m''[`tidx', 7] = `p10'
                }
            }
        }

        restore

    }

    // =========================================================================
    // 4. DISPLAY RESULTS TABLE
    // =========================================================================
    if "`notable'" == "" {
        di as txt ""
        di as txt "{hline 78}"
        di as txt "{bf:Simulated Critical Bounds for ARDL F-Test}"
        di as txt "{hline 78}"

        forvalues m = 1/5 {
            if `m' == 1 local mlbl "No intercept, No trend"
            if `m' == 2 local mlbl "Restricted intercept, No trend"
            if `m' == 3 local mlbl "Unrestricted intercept, No trend"
            if `m' == 4 local mlbl "Unrestricted intercept, Restricted trend"
            if `m' == 5 local mlbl "Unrestricted intercept, Unrestricted trend"

            di as txt ""
            di as txt "{bf:Model `m': `mlbl'}"
            di as txt "{hline 78}"
            di as txt _col(3) "{bf:T}" ///
               _col(10) "{bf:{&alpha}=0.01}" ///
               _col(32) "{bf:{&alpha}=0.05}" ///
               _col(54) "{bf:{&alpha}=0.10}"
            di as txt _col(8) "{bf:I(0)}" _col(18) "{bf:I(1)}" ///
               _col(30) "{bf:I(0)}" _col(40) "{bf:I(1)}" ///
               _col(52) "{bf:I(0)}" _col(62) "{bf:I(1)}"
            di as txt "{hline 78}"

            forvalues i = 1/`ntobs' {
                local tt = el(`cv_m`m'', `i', 1)
                local lb1 = el(`cv_m`m'', `i', 2)
                local ub1 = el(`cv_m`m'', `i', 3)
                local lb5 = el(`cv_m`m'', `i', 4)
                local ub5 = el(`cv_m`m'', `i', 5)
                local lb10 = el(`cv_m`m'', `i', 6)
                local ub10 = el(`cv_m`m'', `i', 7)

                di as res _col(3) %4.0f `tt' ///
                   _col(9) %8.4f `lb1' _col(19) %8.4f `ub1' ///
                   _col(31) %8.4f `lb5' _col(41) %8.4f `ub5' ///
                   _col(53) %8.4f `lb10' _col(63) %8.4f `ub10'
            }
            di as txt "{hline 78}"
        }

        // ============ PSS vs Simulated Comparison Table ============
        di as txt ""
        di as txt "{hline 78}"
        di as txt "{bf:Comparison: PSS (2001) Asymptotic vs Simulated Critical Values}"
        di as txt "{hline 78}"

        forvalues m = 1/5 {
            if `m' == 1 local mlbl "No intercept, No trend"
            if `m' == 2 local mlbl "Restricted intercept, No trend"
            if `m' == 3 local mlbl "Unrestricted intercept, No trend"
            if `m' == 4 local mlbl "Unrestricted intercept, Restricted trend"
            if `m' == 5 local mlbl "Unrestricted intercept, Unrestricted trend"

            // PSS (2001) asymptotic CVs for k=`nregs' (hardcoded for k=1)
            if `nregs' == 1 {
                if `m' == 1 {
                    local pss_lb1 = 4.94
                    local pss_ub1 = 6.84
                    local pss_lb5 = 3.62
                    local pss_ub5 = 4.94
                    local pss_lb10 = 3.02
                    local pss_ub10 = 4.04
                }
                if `m' == 2 {
                    local pss_lb1 = 6.84
                    local pss_ub1 = 8.74
                    local pss_lb5 = 4.94
                    local pss_ub5 = 6.56
                    local pss_lb10 = 4.04
                    local pss_ub10 = 5.62
                }
                if `m' == 3 {
                    local pss_lb1 = 6.84
                    local pss_ub1 = 7.84
                    local pss_lb5 = 4.94
                    local pss_ub5 = 5.73
                    local pss_lb10 = 4.04
                    local pss_ub10 = 4.85
                }
                if `m' == 4 {
                    local pss_lb1 = 7.52
                    local pss_ub1 = 8.74
                    local pss_lb5 = 5.59
                    local pss_ub5 = 6.56
                    local pss_lb10 = 4.68
                    local pss_ub10 = 5.62
                }
                if `m' == 5 {
                    local pss_lb1 = 7.52
                    local pss_ub1 = 8.74
                    local pss_lb5 = 5.59
                    local pss_ub5 = 6.56
                    local pss_lb10 = 4.68
                    local pss_ub10 = 5.62
                }
            }
            else {
                // For k>1, use generic approximations
                local pss_lb1 = .
                local pss_ub1 = .
                local pss_lb5 = .
                local pss_ub5 = .
                local pss_lb10 = .
                local pss_ub10 = .
            }

            di as txt ""
            di as txt "{bf:Model `m': `mlbl'}"
            di as txt "{hline 78}"
            di as txt _col(3) "{bf:T}" ///
                _col(10) "{bf:Sig.}" ///
                _col(18) "{bf:PSS I(0)}" ///
                _col(30) "{bf:Sim I(0)}" ///
                _col(41) "{bf:Diff%}" ///
                _col(50) "{bf:PSS I(1)}" ///
                _col(62) "{bf:Sim I(1)}" ///
                _col(73) "{bf:Diff%}"
            di as txt "{hline 78}"

            forvalues i = 1/`ntobs' {
                local tt = el(`cv_m`m'', `i', 1)

                // 1% level
                local slb1 = el(`cv_m`m'', `i', 2)
                local sub1 = el(`cv_m`m'', `i', 3)
                local dlb1 = .
                local dub1 = .
                if `pss_lb1' != . & `pss_lb1' > 0 {
                    local dlb1 = (`slb1' - `pss_lb1') / `pss_lb1' * 100
                    local dub1 = (`sub1' - `pss_ub1') / `pss_ub1' * 100
                }

                // 5% level
                local slb5 = el(`cv_m`m'', `i', 4)
                local sub5 = el(`cv_m`m'', `i', 5)
                local dlb5 = .
                local dub5 = .
                if `pss_lb5' != . & `pss_lb5' > 0 {
                    local dlb5 = (`slb5' - `pss_lb5') / `pss_lb5' * 100
                    local dub5 = (`sub5' - `pss_ub5') / `pss_ub5' * 100
                }

                // 10% level
                local slb10 = el(`cv_m`m'', `i', 6)
                local sub10 = el(`cv_m`m'', `i', 7)
                local dlb10 = .
                local dub10 = .
                if `pss_lb10' != . & `pss_lb10' > 0 {
                    local dlb10 = (`slb10' - `pss_lb10') / `pss_lb10' * 100
                    local dub10 = (`sub10' - `pss_ub10') / `pss_ub10' * 100
                }

                // Display 3 rows per T (one per significance level)
                if `pss_lb1' != . {
                    di as res _col(3) %4.0f `tt' ///
                        _col(11) "1%" ///
                        _col(17) %8.3f `pss_lb1' ///
                        _col(29) %8.3f `slb1' ///
                        _col(40) %6.1f `dlb1' ///
                        _col(49) %8.3f `pss_ub1' ///
                        _col(61) %8.3f `sub1' ///
                        _col(72) %6.1f `dub1'
                    di as res _col(3) "    " ///
                        _col(11) "5%" ///
                        _col(17) %8.3f `pss_lb5' ///
                        _col(29) %8.3f `slb5' ///
                        _col(40) %6.1f `dlb5' ///
                        _col(49) %8.3f `pss_ub5' ///
                        _col(61) %8.3f `sub5' ///
                        _col(72) %6.1f `dub5'
                    di as res _col(3) "    " ///
                        _col(11) "10%" ///
                        _col(17) %8.3f `pss_lb10' ///
                        _col(29) %8.3f `slb10' ///
                        _col(40) %6.1f `dlb10' ///
                        _col(49) %8.3f `pss_ub10' ///
                        _col(61) %8.3f `sub10' ///
                        _col(72) %6.1f `dub10'
                    di as txt "    {hline 74}"
                }
                else {
                    di as txt _col(3) %4.0f `tt' ///
                        _col(10) "  PSS asymptotic CVs not available for k=`nregs'"
                }
            }
            di as txt "{hline 78}"
        }

        di as txt ""
        di as txt _col(2) "PSS = Pesaran, Shin & Smith (2001) asymptotic (T -> infinity)"
        di as txt _col(2) "Sim = Monte Carlo (`nsim' replications)"
        di as txt _col(2) "Diff% = (Sim - PSS)/PSS x 100"
        di as txt _col(2) "Positive Diff% = simulated CV is higher (more conservative)"
        di as txt _col(2) "Negative Diff% = simulated CV is lower (asymptotic too conservative)"
    }

    // =========================================================================
    // 5. STORE RESULTS
    // =========================================================================
    forvalues m = 1/5 {
        mat colnames `cv_m`m'' = T LB_1pct UB_1pct LB_5pct UB_5pct LB_10pct UB_10pct
        return matrix cv_model`m' = `cv_m`m''
    }

    return scalar nsim = `nsim'
    return scalar nregs = `nregs'
    return scalar maxlag = `maxlag'
    return local ic "`ic'"
    return local tobs "`tobs'"

end