*! caustests v1.1.0  22feb2026
*! Granger Causality Tests: TY, Fourier (Single/Cumulative), and Quantile Variants
*! Author: Dr Merwan Roudane (merwanroudane920@gmail.com)
*!
*! References:
*!  Toda & Yamamoto (1995) J.Econometrics 66:225-250
*!  Enders & Jones (2016) Stud.Nonlinear Dyn.Econometrics
*!  Nazlioglu, Gormus, Soytas (2016) Energy Economics 60:168-175
*!  Nazlioglu, Soytas, Gormus (2019) Emerg.Mkts Finance Trade 55:105-117
*!  Cai, Chang, Xiang, Chang (2023) Finance Research Letters 58:104327
*!  Cheng, Hsueh, Ranjbar, Wang, Chang (2021) Letters Spatial Resource Sci 14:31-49

program define caustests, eclass
    version 14
    syntax varlist(min=2 ts) [if] [in] , ///
        TEST(integer)            ///
        [ PMAX(integer 8)        ///
          IC(integer 1)          ///
          NBOOT(integer 1000)    ///
          KMAX(integer 3)        ///
          DMAX(integer -1)       ///
          Quantiles(numlist >0 <1 ascending) ///
          GRAPH                  ///
          noHEADer               ///
          SCHEME(string)         ]

    //????????? Validate ???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
    if !inrange(`test',1,7) {
        di as err "test() must be between 1 and 7"
        exit 198
    }
    if !inlist(`ic',1,2) {
        di as err "ic() must be 1 (AIC) or 2 (SBC)"
        exit 198
    }
    if `nboot' < 99 {
        di as err "nboot() must be >= 99"
        exit 198
    }

    // Auto dmax: 0 for Granger in differences, 1 for TY in levels
    if `dmax' == -1 local dmax = cond(inlist(`test',2,4), 0, 1)

    // Default quantiles
    if "`quantiles'" == "" local quantiles "0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90"

    // Default graph scheme
    if "`scheme'" == "" local scheme "s1color"

    //????????? Get time setting ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
    qui tsset
    local timevar `r(timevar)'

    //????????? Prepare data ???????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
    marksample touse
    preserve
    qui keep if `touse'
    qui keep `varlist' `timevar'
    qui tsset `timevar'
    // Remove listwise missing
    foreach v of local varlist {
        qui drop if missing(`v')
    }
    local T = _N
    if `T' < `=2*`pmax'+4' {
        di as err "Insufficient observations (T=`T') for pmax=`pmax'"
        restore
        exit 2001
    }

    //????????? Note kmax only relevant for tests 2-5, 7 ????????????????????????????????????????????????????????????????????????????????????
    if !inlist(`test',2,3,4,5,7) local kmax 0

    //????????? Run tests ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
    if inlist(`test',1,2,3,4,5) {
        // OLS path via Mata
        mata: _ct_data   = st_data(., "`varlist'")
        mata: _ct_vnames = tokens("`varlist'")'
        mata: _caustests_ols(_ct_data, _ct_vnames, `pmax', `ic', `test', `nboot', `kmax', `dmax')
        mata: mata drop _ct_data _ct_vnames
    }
    else {
        // Quantile path via Stata
        _ct_qrun `varlist' , test(`test') pmax(`pmax') ic(`ic') ///
            nboot(`nboot') kmax(`kmax') dmax(`dmax') ///
            quantiles(`quantiles')
    }

    restore

    //????????? Display ??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
    if "`header'" == "" _ct_header `test'

    if inlist(`test',1,2,3,4,5) {
        _ct_table_ols `test'
    }
    else {
        _ct_table_q `test' "`quantiles'"
        if "`graph'" != "" _ct_graph_q `test' "`quantiles'" "`scheme'"
    }

end // caustests


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// HEADER
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_header
    args test
    local sep "{hline 70}"
    if      `test'==1 local ttl "Toda & Yamamoto (1995) Granger Causality Test"
    else if `test'==2 local ttl "Single Fourier-Frequency Granger Causality Test"
    else if `test'==3 local ttl "Single Fourier-Frequency Toda-Yamamoto Causality Test"
    else if `test'==4 local ttl "Cumulative Fourier-Frequency Granger Causality Test"
    else if `test'==5 local ttl "Cumulative Fourier-Frequency Toda-Yamamoto Causality Test"
    else if `test'==6 local ttl "Toda-Yamamoto Granger Causality Test in Quantiles"
    else              local ttl "Bootstrap Fourier Granger Causality in Quantiles (BFGC-Q)"

    if      `test'==1 local ref "(Toda & Yamamoto, 1995)"
    else if `test'==2 local ref "(Enders & Jones, 2016)"
    else if `test'==3 local ref "(Nazlioglu, Gormus & Soytas, 2016)"
    else if `test'==4 local ref "(Enders & Jones, 2019)"
    else if `test'==5 local ref "(Nazlioglu, Soytas & Gormus, 2019)"
    else if `test'==6 local ref "(Cai, Chang, Xiang & Chang, 2023)"
    else              local ref "(Cheng, Hsueh, Ranjbar, Wang & Chang, 2021)"

    di ""
    di as txt "`sep'"
    di as res _col(4) "`ttl'"
    di as txt _col(4) "`ref'"
    di as txt "`sep'"
end


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// OLS RESULT TABLE  (tests 1-5)
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_table_ols
    args test
    tempname W PA PB LG FR
    mat `W'  = e(wald)
    mat `PA' = e(pval_a)
    mat `PB' = e(pval_b)
    mat `LG' = e(lags)
    mat `FR' = e(freq)
    local nr = rowsof(`W')
    local rn : rownames `W'

    di as txt _col(2) "{hline 72}"
    di as txt _col(2) %-30s "Null Hypothesis"          ///
               _col(34) %8s "Wald"                     ///
               _col(44) %9s "Asym.p"                   ///
               _col(55) %9s "Boot.p"                   ///
               _col(66) %3s "Lag"                      ///
               _col(71) %4s "Freq"
    di as txt _col(2) "{hline 72}"

    forval r = 1/`nr' {
        // direction labels may contain "=>" ??? use tokenize instead of :word
        local h  : word `r' of `rn'
        local w  = `W'[`r',1]
        local pa = `PA'[`r',1]
        local pb = `PB'[`r',1]
        local lg = `LG'[`r',1]
        local fr = `FR'[`r',1]
        local st = cond(`pb'<0.01,"***",cond(`pb'<0.05,"**",cond(`pb'<0.10,"*","")))
        // truncate label to 30 chars if needed
        if length("`h'") > 30 local h = substr("`h'",1,30)

        di as txt _col(2) %-30s "`h'"                  ///
           as res _col(32) %8.3f `w'                   ///
           as txt _col(42) %9.3f `pa'                  ///
           as res _col(53) %9.3f `pb'                  ///
                  _col(64) "`st'"                      ///
           as txt _col(68) %3.0f `lg'                  ///
                  _col(73) %4.0f `fr'
    }
    di as txt _col(2) "{hline 72}"
    di as txt _col(2) "Note: *p<0.10, **p<0.05, ***p<0.01 based on bootstrap p-values"
    if inlist(`test',2,3,4,5) di as txt _col(2) "      Freq = optimal Fourier frequency (k*)"
    if `test'==1 di as txt _col(2) "      Freq = 0 (no Fourier terms)"
    di ""
end


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// QUANTILE RESULT TABLE  (tests 6-7)
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_table_q
    args test quantiles
    local nq : word count `quantiles'

    tempname Wq PBq LGq FRq
    mat `Wq'  = e(wald_q)
    mat `PBq' = e(pval_bq)
    mat `LGq' = e(lags_q)
    if `test'==7 {
        capture mat `FRq' = e(freq_q)
        if _rc mat `FRq' = J(`ndirs',1,0)
    }

    local ndirs = rowsof(`Wq')
    local dirs  : rownames `Wq'

    forval d = 1/`ndirs' {
        local hn : word `d' of `dirs'
        local lg = `LGq'[`d',1]
        if `test'==7 {
            local fr = `FRq'[`d',1]
            di as txt _col(2) "H0: `hn'" as txt "  (Lag=`lg', Freq=`fr')"
        }
        else {
            di as txt _col(2) "H0: `hn'" as txt "  (Lag=`lg')"
        }

        di as txt _col(2) "{hline 46}"
        di as txt _col(2) %10s "Quantile " _col(14) %9s "Wald" _col(25) %11s "Boot p-val" _col(40) "Sig."
        di as txt _col(2) "{hline 46}"

        forval q = 1/`nq' {
            local tau : word `q' of `quantiles'
            local w   = `Wq'[`d',`q']
            local pb  = `PBq'[`d',`q']
            local st  = cond(`pb'<0.01,"***",cond(`pb'<0.05,"**",cond(`pb'<0.10,"*","")))
            di as txt _col(2) %10.2f `tau' ///
               as res _col(14) %9.3f  `w'  ///
               as txt _col(25) %11.3f `pb' ///
               as res _col(40) "`st'"
        }
        di as txt _col(2) "{hline 46}"
        di ""
    }
    di as txt "Note: *p<0.10, **p<0.05, ***p<0.01 based on bootstrap p-values"
    if `test'==7 di as txt "      Freq = optimal Fourier frequency (k*)"
    di ""
end


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// QUANTILE GRAPH  (tests 6-7)
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_graph_q
    args test quantiles scheme

    local nq    : word count `quantiles'
    local ndirs = rowsof(e(wald_q))
    local dirs  : rownames e(wald_q)

    tempname Wq PBq
    mat `Wq'  = e(wald_q)
    mat `PBq' = e(pval_bq)

    preserve
    qui drop _all
    qui set obs `nq'

    gen double tau = .
    forval q = 1/`nq' {
        local t : word `q' of `quantiles'
        qui replace tau = `t' in `q'
    }

    local colors navy maroon dkgreen orange purple teal
    local glist_wald ""
    local glist_pval ""

    forval d = 1/`ndirs' {
        local dir : word `d' of `dirs'
        local col : word `d' of `colors'
        if "`col'" == "" local col navy
        local safedn = subinstr("`dir'","=>"," => ",.)

        gen double w_d`d' = .
        gen double p_d`d' = .
        forval q = 1/`nq' {
            qui replace w_d`d' = `Wq'[`d',`q']  in `q'
            qui replace p_d`d' = `PBq'[`d',`q'] in `q'
        }
        // Wald graph
        twoway ///
            (bar w_d`d' tau, barwidth(0.07) fcolor(`col'%50) lcolor(`col') lwidth(vthin)) ///
            (connected w_d`d' tau, mcolor(`col') lcolor(`col') msymbol(circle) msize(small) lwidth(medthick)), ///
            xlabel(`quantiles', format(%4.2f) angle(45)) ///
            xtitle("Quantile ({it: {&tau}})", size(small)) ///
            ytitle("Wald statistic", size(small)) ///
            title("`safedn'", color(`col') size(small)) ///
            yline(3.841, lcolor(red%70) lwidth(thin) lpattern(dash)) ///
            note("Horizontal dashed line: chi2(p) 5% critical value (3.84)", size(vsmall)) ///
            legend(off) scheme(`scheme') name(ct_w`d', replace) nodraw

        // P-value graph
        twoway ///
            (area p_d`d' tau, fcolor(`col'%15) lwidth(none)) ///
            (connected p_d`d' tau, mcolor(`col') lcolor(`col') msymbol(circle) msize(small) lwidth(medthick)), ///
            xlabel(`quantiles', format(%4.2f) angle(45)) ///
            xtitle("Quantile ({it: {&tau}})", size(small)) ///
            ytitle("Bootstrap p-value", size(small)) ///
            title("`safedn'", color(`col') size(small)) ///
            yline(0.05, lcolor(red%70)    lwidth(thin) lpattern(dash)) ///
            yline(0.10, lcolor(orange%70) lwidth(thin) lpattern(shortdash)) ///
            yscale(range(0 1)) ylabel(0(0.1)1, format(%3.1f)) ///
            note("Red: 5% level | Orange: 10% level", size(vsmall)) ///
            legend(off) scheme(`scheme') name(ct_p`d', replace) nodraw

        local glist_wald "`glist_wald' ct_w`d'"
        local glist_pval "`glist_pval' ct_p`d'"
    }

    // Combine Wald graphs
    local wald_title = cond(`test'==6,"Quantile TY - Wald Statistics","BFGC-Q - Wald Statistics")
    graph combine `glist_wald', ///
        title("`wald_title'", size(medsmall) margin(b=2)) ///
        scheme(`scheme') ycommon name(ct_wald_combined, replace)

    // Combine p-value graphs
    local pval_title = cond(`test'==6,"Quantile TY - Bootstrap p-values","BFGC-Q - Bootstrap p-values")
    graph combine `glist_pval', ///
        title("`pval_title'", size(medsmall) margin(b=2)) ///
        scheme(`scheme') name(ct_pval_combined, replace)

    di as txt ""
    di as txt "  Graphs available:"
    di as txt "    {cmd:graph display ct_wald_combined}   (Wald statistics)"
    di as txt "    {cmd:graph display ct_pval_combined}   (Bootstrap p-values)"

    restore
end


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// QUANTILE TEST RUNNER  (dispatch all directions & quantiles)
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_qrun, eclass
    syntax varlist(min=2 ts) , TEST(integer) PMAX(integer) IC(integer) ///
        NBOOT(integer) KMAX(integer) DMAX(integer) Quantiles(numlist)

    local nvar : word count `varlist'
    local nq   : word count `quantiles'

    // Build direction list in same order: for each depvar i, for each cause j???i
    local ndirs = `nvar' * (`nvar'-1)
    local dirlist ""
    forval i = 1/`nvar' {
        forval j = 1/`nvar' {
            if `i'==`j' continue
            local xn : word `j' of `varlist'
            local yn : word `i' of `varlist'
            local dirlist "`dirlist' `xn'=>`yn'"
        }
    }

    // Result matrices
    tempname Wq PBq LGq FRq
    mat `Wq'  = J(`ndirs',`nq',.)
    mat `PBq' = J(`ndirs',`nq',.)
    mat `LGq' = J(`ndirs',1,.)
    mat `FRq' = J(`ndirs',1,0)
    mat rownames `Wq'  = `dirlist'
    mat rownames `PBq' = `dirlist'
    mat rownames `LGq' = `dirlist'
    mat rownames `FRq' = `dirlist'

    // Build col names from quantiles
    local qcols ""
    foreach tau of local quantiles {
        local qs = subinstr("`tau'",".","p",.)
        local qcols "`qcols' `qs'"
    }
    mat colnames `Wq'  = `qcols'
    mat colnames `PBq' = `qcols'

    local dd = 0
    forval i = 1/`nvar' {
        local yn : word `i' of `varlist'
        forval j = 1/`nvar' {
            if `i'==`j' continue
            local dd = `dd'+1
            local xn : word `j' of `varlist'

            // Controls
            local ctrls ""
            forval c = 1/`nvar' {
                if `c'==`i' | `c'==`j' continue
                local cn : word `c' of `varlist'
                local ctrls "`ctrls' `cn'"
            }

            // Progress message
            di as txt "  Testing: `xn' => `yn' ..."

            // Step 1: Select optimal p (& k for test 7)
            _ct_qselect `yn' `xn' `ctrls', pmax(`pmax') ic(`ic') ///
                test(`test') kmax(`kmax') dmax(`dmax')
            local popt = r(popt)
            local kopt = r(kopt)
            mat `LGq'[`dd',1] = `popt'
            if `test'==7 mat `FRq'[`dd',1] = `kopt'

            // Step 2: Loop over quantiles
            local qq = 0
            foreach tau of local quantiles {
                local qq = `qq'+1
                _ct_qtest `yn' `xn' `ctrls', ///
                    popt(`popt') kopt(`kopt') dmax(`dmax') ///
                    tau(`tau') nboot(`nboot') test(`test')
                mat `Wq'[`dd',`qq']  = r(wald)
                mat `PBq'[`dd',`qq'] = r(pval_b)
            }
        }
    }

    ereturn clear
    ereturn matrix wald_q  = `Wq'
    ereturn matrix pval_bq = `PBq'
    ereturn matrix lags_q  = `LGq'
    if `test'==7 ereturn matrix freq_q = `FRq'
    ereturn scalar test = `test'
    ereturn local  cmd  "caustests"
end


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// OPTIMAL LAG (& FREQUENCY) FOR QUANTILE TESTS
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_qselect, rclass
    syntax varlist(min=1), PMAX(integer) IC(integer) TEST(integer) ///
        KMAX(integer) DMAX(integer)

    local depv : word 1 of `varlist'
    local nv   : word count `varlist'
    local best_ic  1e15
    local popt 1
    local kopt 0

    if `test'==6 {
        // No Fourier: select p only
        forval p = 1/`pmax' {
            local plg = `p'+`dmax'
            local reg ""
            foreach v of local varlist {
                forval l = 1/`plg' {
                    local reg "`reg' L`l'.`v'"
                }
            }
            qui capture regress `depv' `reg'
            if _rc continue
            local n = e(N)
            local k = e(df_model)+1
            local rss = e(rss)
            local aic = log(`rss'/`n') + 2*`k'/`n'
            local sbc = log(`rss'/`n') + `k'*log(`n')/`n'
            local icv = cond(`ic'==1,`aic',`sbc')
            if `icv' < `best_ic' {
                local best_ic `icv'
                local popt `p'
            }
        }
    }
    else {
        // Test 7: select (k,p) jointly with cumulative Fourier terms
        forval k = 1/`kmax' {
            forval p = 1/`pmax' {
                local plg = `p'+`dmax'
                // Generate cumulative Fourier terms: sin/cos for j=1,...,k
                local fourier_sel ""
                forval j = 1/`k' {
                    tempvar sf_`j' ck_`j'
                    qui gen double `sf_`j'' = sin(2*_pi*`j'*_n/_N)
                    qui gen double `ck_`j'' = cos(2*_pi*`j'*_n/_N)
                    local fourier_sel "`fourier_sel' `sf_`j'' `ck_`j''"
                }
                local reg ""
                foreach v of local varlist {
                    forval l = 1/`plg' {
                        local reg "`reg' L`l'.`v'"
                    }
                }
                qui capture regress `depv' `reg' `fourier_sel'
                forval j = 1/`k' {
                    cap drop `sf_`j'' `ck_`j''
                }
                if _rc continue
                local n = e(N)
                local k2 = e(df_model)+1
                local rss = e(rss)
                local aic = log(`rss'/`n') + 2*`k2'/`n'
                local sbc = log(`rss'/`n') + `k2'*log(`n')/`n'
                local icv = cond(`ic'==1,`aic',`sbc')
                if `icv' < `best_ic' {
                    local best_ic `icv'
                    local popt `p'
                    local kopt `k'
                }
            }
        }
    }
    return scalar popt = `popt'
    return scalar kopt = `kopt'
end


// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
// QUANTILE WALD + BOOTSTRAP FOR ONE (DIRECTION, TAU)
// ????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????
program define _ct_qtest, rclass
    syntax varlist(min=2), POPT(integer) KOPT(integer) DMAX(integer) ///
        TAU(real) NBOOT(integer) TEST(integer)

    local depv  : word 1 of `varlist'
    local causv : word 2 of `varlist'
    local nv    : word count `varlist'
    local ctrls ""
    forval c = 3/`nv' {
        local ctrls "`ctrls' `: word `c' of `varlist''"
    }

    local plg = `popt'+`dmax'

    // Regressor lists
    local dep_lags ""
    forval l = 1/`plg' {
        local dep_lags "`dep_lags' L`l'.`depv'"
    }
    local cause_test ""
    forval l = 1/`popt' {
        local cause_test "`cause_test' L`l'.`causv'"
    }
    local cause_extra ""
    if `dmax' > 0 {
        forval l = `=`popt'+1'/`plg' {
            local cause_extra "`cause_extra' L`l'.`causv'"
        }
    }
    local ctrl_lags ""
    foreach v of local ctrls {
        forval l = 1/`plg' {
            local ctrl_lags "`ctrl_lags' L`l'.`v'"
        }
    }

    // Fourier terms (test 7 only)
    local fourvars ""
    if `test'==7 & `kopt'>0 {
        forval j = 1/`kopt' {
            tempvar sf_`j' ck_`j'
            qui gen double `sf_`j'' = sin(2*_pi*`j'*_n/_N)
            qui gen double `ck_`j'' = cos(2*_pi*`j'*_n/_N)
            local fourvars "`fourvars' `sf_`j'' `ck_`j''"
        }
    }

    // Convert time-series operators to tempvars for qreg
    local all_ts "`dep_lags' `cause_test' `cause_extra' `ctrl_lags'"
    if "`all_ts'" != "" {
        qui tsrevar `all_ts'
        local all_tmp `r(varlist)'
    }
    else local all_tmp ""

    // Extract cause_test tempvars
    local c_start = `plg' + 1
    local c_end   = `plg' + `popt'
    local cause_tmp ""
    forval idx = `c_start'/`c_end' {
        local v : word `idx' of `all_tmp'
        local cause_tmp "`cause_tmp' `v'"
    }

    // Extract rest variables (everything except cause_tmp)
    local rest_tmp ""
    local n_all : word count `all_tmp'
    forval idx = 1/`n_all' {
        if `idx' >= `c_start' & `idx' <= `c_end' continue
        local v : word `idx' of `all_tmp'
        local rest_tmp "`rest_tmp' `v'"
    }

    local fullreg "`all_tmp' `fourvars'"
    local restreg "`rest_tmp' `fourvars'"

    // ── Observed Wald from qreg at tau ────────────────────────
    qui capture qreg `depv' `fullreg', quantile(`tau')
    if _rc {
        return scalar wald   = .
        return scalar pval_b = .
        exit
    }
    tempname bvec Vsub W_mat
    matrix `bvec' = J(`popt',1,.)
    matrix `Vsub' = J(`popt',`popt',.)
    forval ii = 1/`popt' {
        local cv : word `ii' of `cause_tmp'
        matrix `bvec'[`ii',1] = _b[`cv']
        forval jj = 1/`popt' {
            local cv2 : word `jj' of `cause_tmp'
            matrix `Vsub'[`ii',`jj'] = e(V)["`cv'","`cv2'"]
        }
    }
    matrix `W_mat' = `bvec'' * invsym(`Vsub') * `bvec'
    local W_obs = `W_mat'[1,1]
    if missing(`W_obs') | `W_obs' < 0 local W_obs = 0

    // ── Restricted qreg for bootstrap residuals ──────────────
    qui capture qreg `depv' `restreg', quantile(`tau')
    if _rc {
        return scalar wald   = `W_obs'
        return scalar pval_b = .
        exit
    }
    tempvar yhat0 resid0 ystar
    qui predict double `yhat0', xb
    qui gen double `resid0' = `depv' - `yhat0'
    qui summ `resid0', meanonly
    qui replace `resid0' = `resid0' - r(mean)
    qui gen double `ystar' = `depv'

    // Count valid observations
    qui count if !missing(`yhat0') & !missing(`resid0')
    local Te = r(N)
    if `Te' < 10 {
        return scalar wald   = `W_obs'
        return scalar pval_b = .
        exit
    }

    // Load valid data into Mata (using yhat0 as touse selector)
    mata: _ct_er = st_data(., "`resid0'", "`yhat0'")
    mata: _ct_yh = st_data(., "`yhat0'", "`yhat0'")

    // ── Bootstrap: qreg-based Wald at tau ────────────────────
    local Wboot_ge = 0
    forval b = 1/`nboot' {
        // Resample residuals and create bootstrap y*
        mata: _ct_Te = rows(_ct_er)
        mata: _ct_idx = ceil(_ct_Te :* runiform(_ct_Te, 1))
        mata: _ct_es = _ct_er[_ct_idx, .]
        mata: _ct_es = _ct_es :- mean(_ct_es)
        mata: _ct_ys = _ct_yh + _ct_es
        mata: st_store(., "`ystar'", "`yhat0'", _ct_ys)

        // Qreg unrestricted on bootstrap y* at same tau
        qui capture qreg `ystar' `fullreg', quantile(`tau')
        if _rc continue

        // Compute bootstrap Wald manually from b and V
        local Wb = 0
        capture {
            matrix `bvec' = J(`popt',1,.)
            matrix `Vsub' = J(`popt',`popt',.)
            forval ii = 1/`popt' {
                local cv : word `ii' of `cause_tmp'
                matrix `bvec'[`ii',1] = _b[`cv']
                forval jj = 1/`popt' {
                    local cv2 : word `jj' of `cause_tmp'
                    matrix `Vsub'[`ii',`jj'] = e(V)["`cv'","`cv2'"]
                }
            }
            matrix `W_mat' = `bvec'' * invsym(`Vsub') * `bvec'
            local Wb = `W_mat'[1,1]
        }
        if missing(`Wb') | `Wb' < 0 local Wb = 0
        if `Wb' >= `W_obs' local Wboot_ge = `Wboot_ge' + 1
    }

    capture mata: mata drop _ct_er _ct_yh _ct_Te _ct_idx _ct_es _ct_ys

    return scalar wald   = `W_obs'
    return scalar pval_b = `Wboot_ge'/`nboot'
end


mata:

void _caustests_ols(real matrix data, string colvector vnames,
    real scalar pmax, real scalar ic, real scalar test,
    real scalar Nboot, real scalar kmax, real scalar dmax)
{
    real scalar nvar, ntests, di_, dj_, idx, dm, km, W, pa, pval_b, r
    real scalar p_opt, k_opt
    real matrix GC, VARy, ylag, VARx, indv, depvar, cause, ctrl, y3col
    real colvector Wboot
    string colvector dirs, indnames
    string scalar rn

    nvar   = cols(data)
    ntests = nvar*(nvar-1)
    GC     = J(ntests, 5, .)
    dirs   = J(ntests, 1, "")
    idx = 0
    for (di_ = 1; di_ <= nvar; di_++) {
        depvar   = data[., di_]
        indv     = _ct_delcol(data, di_)
        indnames = _ct_delrowstr(vnames, di_)
        for (dj_ = 1; dj_ <= nvar-1; dj_++) {
            idx++
            cause = indv[., dj_]
            ctrl  = _ct_delcol(indv, dj_)
            y3col = depvar, cause, ctrl
            if (test==1) {
                dm = 1
                km = 0
            }
            else if (test==2) {
                dm = 0
                km = kmax
            }
            else if (test==3) {
                dm = 1
                km = kmax
            }
            else if (test==4) {
                dm = 0
                km = kmax
            }
            else {
                dm = 1
                km = kmax
            }
            if (dmax >= 0) dm = dmax
            p_opt = 1
            k_opt = 0
            _ct_select_kp(y3col, pmax, km, ic, test, rows(data), p_opt, k_opt)
            _ct_make_lags(y3col, p_opt+dm, VARy, ylag)
            VARx   = (ylag, _ct_det(rows(VARy), test, k_opt, p_opt, dm, rows(data)))
            W      = _ct_wald(VARy[.,1], VARx, p_opt, dm)
            pa     = chi2tail(p_opt, max((W,0)))
            Wboot  = _ct_bootstrap(VARy[.,1], VARx, p_opt, dm, Nboot)
            pval_b = mean(Wboot :>= W)
            GC[idx,.] = W, pa, pval_b, p_opt, k_opt
            dirs[idx] = strtrim(indnames[dj_]) + "=>" + strtrim(vnames[di_])
        }
    }
    st_matrix("_ct_W",  GC[.,1])
    st_matrix("_ct_PA", GC[.,2])
    st_matrix("_ct_PB", GC[.,3])
    st_matrix("_ct_LG", GC[.,4])
    st_matrix("_ct_FR", GC[.,5])
    rn = ""
    for (r=1; r<=ntests; r++) rn = rn + (r>1 ? " " : "") + dirs[r]
    stata(`"matrix rownames _ct_W  = "' + rn)
    stata(`"matrix rownames _ct_PA = "' + rn)
    stata(`"matrix rownames _ct_PB = "' + rn)
    stata(`"matrix rownames _ct_LG = "' + rn)
    stata(`"matrix rownames _ct_FR = "' + rn)
    stata("ereturn clear")
    stata("ereturn matrix wald   = _ct_W")
    stata("ereturn matrix pval_a = _ct_PA")
    stata("ereturn matrix pval_b = _ct_PB")
    stata("ereturn matrix lags   = _ct_LG")
    stata("ereturn matrix freq   = _ct_FR")
    stata(`"ereturn scalar test  = "' + strofreal(test))
    stata(`"ereturn local  cmd   caustests"')
}

void _ct_select_kp(real matrix y, real scalar pmax, real scalar kmax,
    real scalar ic, real scalar test, real scalar T_orig,
    real scalar p_opt, real scalar k_opt)
{
    real scalar p, k, icv, best_ic, T_eff, ki
    real matrix dep, xl, d, z, t_seq, cf
    best_ic = 1e15
    p_opt   = 1
    k_opt   = 0
    if (kmax == 0) {
        for (p=1; p<=pmax; p++) {
            _ct_make_lags(y, p, dep, xl)
            d   = J(rows(dep), 1, 1)
            z   = xl, d
            icv = _ct_ic(dep, z, p, ic)
            if (icv < best_ic) {
                best_ic = icv
                p_opt = p
            }
        }
    }
    else {
        for (k=1; k<=kmax; k++) {
            for (p=1; p<=pmax; p++) {
                _ct_make_lags(y, p, dep, xl)
                T_eff = rows(dep)
                t_seq = (1::T_eff)
                if (test==2 || test==3) {
                    cf = sin(2*pi()*k*t_seq/T_eff), cos(2*pi()*k*t_seq/T_eff)
                    d  = J(T_eff, 1, 1), cf
                }
                else {
                    cf = J(T_eff, 0, .)
                    for (ki=1; ki<=k; ki++) {
                        cf = cf, sin(2*pi()*ki*t_seq/T_eff), cos(2*pi()*ki*t_seq/T_eff)
                    }
                    d = J(T_eff, 1, 1), cf
                }
                z   = xl, d
                icv = _ct_ic(dep, z, p, ic)
                if (icv < best_ic) {
                    best_ic = icv
                    p_opt = p
                    k_opt = k
                }
            }
        }
    }
}

real scalar _ct_ic(real matrix dep, real matrix z, real scalar p, real scalar ic_type)
{
    real scalar T, nk, ldv, aic, sbc
    real matrix b, e, VCV
    T   = rows(dep)
    nk  = cols(dep)
    b   = invsym(z'z)*(z'dep)
    e   = dep - z*b
    VCV = (e'e)/T
    ldv = log(det(VCV))
    aic = ldv + (2  /T)*(nk*nk*p + nk) + nk*(1+log(2*pi()))
    sbc = ldv + (log(T)/T)*(nk*nk*p + nk) + nk*(1+log(2*pi()))
    return(ic_type==1 ? aic : sbc)
}

void _ct_make_lags(real matrix y, real scalar plags, real matrix ya, real matrix yl)
{
    real scalar k, T, i, j, col_
    real matrix ylag
    k    = cols(y)
    T    = rows(y)
    ylag = J(T, k*plags, 0)
    for (i=1; i<=k; i++) {
        for (j=1; j<=plags; j++) {
            col_ = j + plags*(i-1)
            ylag[(j+1)::T, col_] = y[1::(T-j), i]
        }
    }
    ya = y[(plags+1)::T, .]
    yl = ylag[(plags+1)::T, .]
}

real matrix _ct_det(real scalar T_eff, real scalar test,
    real scalar k_freq, real scalar p, real scalar dm, real scalar T_orig)
{
    real scalar p_full, ki
    real matrix d, t_full, cf
    p_full = p + dm
    if (test==1) {
        d = J(T_eff, 1, 1)
    }
    else if (test==2 || test==3) {
        t_full = (1::T_orig)
        cf = sin(2*pi()*k_freq*t_full/T_orig), cos(2*pi()*k_freq*t_full/T_orig)
        d = J(T_eff,1,1), cf[(p_full+1)::T_orig, .]
    }
    else if (test==4 || test==5) {
        t_full = (1::T_orig)
        cf     = J(T_orig, 0, .)
        for (ki=1; ki<=k_freq; ki++) {
            cf = cf, sin(2*pi()*ki*t_full/T_orig), cos(2*pi()*ki*t_full/T_orig)
        }
        d = J(T_eff,1,1), cf[(p_full+1)::T_orig, .]
    }
    else {
        d = J(T_eff, 1, 1)
    }
    return(d)
}

real scalar _ct_wald(real matrix y, real matrix z, real scalar p, real scalar dm)
{
    real scalar pf, nc, df, RSS_ur, RSS_r, F, W
    real matrix b, u, zr, br, ur
    pf = p + dm
    nc = cols(z)
    b      = invsym(z'z)*(z'y)
    u      = y - z*b
    RSS_ur = (u'u)[1,1]
    if (dm == 0) {
        if (2*pf+1 <= nc) {
            zr = z[., 1::pf], z[., (2*pf+1)::nc]
        }
        else {
            zr = z[., 1::pf]
        }
    }
    else {
        if (2*pf <= nc) {
            zr = z[., 1::pf], z[., (2*pf)::nc]
        }
        else {
            zr = z[., 1::pf]
        }
    }
    br    = invsym(zr'zr)*(zr'y)
    ur    = y - zr*br
    RSS_r = (ur'ur)[1,1]
    df = rows(z) - nc
    if (df <= 0 || missing(df)) return(.)
    F = ((RSS_r - RSS_ur)/p) / (RSS_ur/df)
    W = F * p
    return(W)
}

real colvector _ct_bootstrap(real matrix y, real matrix z, real scalar p, real scalar dm, real scalar Nboot)
{
    real scalar pf, nc, Te, s
    real matrix b, u, zr, br, zhat, adje, Wboot, e_star, idx_v, ystar
    pf = p + dm
    nc = cols(z)
    Te = rows(z)
    b  = invsym(z'z)*(z'y)
    u  = y - z*b
    if (dm == 0) {
        if (2*pf+1 <= nc) {
            zr = z[.,1::pf], z[.,(2*pf+1)::nc]
        }
        else {
            zr = z[.,1::pf]
        }
    }
    else {
        if (2*pf <= nc) {
            zr = z[.,1::pf], z[.,(2*pf)::nc]
        }
        else {
            zr = z[.,1::pf]
        }
    }
    br    = invsym(zr'zr)*(zr'y)
    zhat  = zr*br[.,1]
    adje  = u[.,1] :- mean(u[.,1])
    Wboot = J(Nboot, 1, .)
    for (s=1; s<=Nboot; s++) {
        idx_v  = ceil(Te :* runiform(Te,1))
        e_star = adje[idx_v]
        e_star = e_star :- mean(e_star)
        ystar  = zhat + e_star
        Wboot[s] = _ct_wald(ystar, z, p, dm)
    }
    _sort(Wboot, 1)
    return(Wboot)
}

real matrix _ct_delcol(real matrix A, real scalar j)
{
    if (cols(A)==1) return(J(rows(A),0,.))
    if (j==1) return(A[.,2::cols(A)])
    else if (j==cols(A)) return(A[.,1::(cols(A)-1)])
    else return(A[.,1::(j-1)], A[.,(j+1)::cols(A)])
}

string colvector _ct_delrowstr(string colvector v, real scalar j)
{
    if (j==1) return(v[2::rows(v)])
    else if (j==rows(v)) return(v[1::(rows(v)-1)])
    else return(v[1::(j-1)] \ v[(j+1)::rows(v)])
}

end