*! _gvar_gc 1.0.1  21aug2026
*! gvar gc -- Granger causality within a country model, with and without the
*! foreign block.
*! Author: Merwan Roudane
*! merwanroudane920@gmail.com
*! https://github.com/merwanroudane
*
* Step -> source map
*   two tests per unit, one on a plain VAR of the unit's own variables and
*   one on a VARX carrying the foreign block as exogenous regressors
*                                             <- GVARX .grangerGVAR
*   the test itself                           <- vars::causality
*       STAT = (R b)' [R (Sigma kron (Z'Z)^-1) R']^-1 (R b) / N
*       df1  = p * #cause * #(non-cause)
*       df2  = K * obs - length(PI)
*   foreign block entered at lags 0..FLag-1   <- embed(Ft, FLag)
*   instantaneous causality                   <- vars::causality, result2
*       lambda = obs * v' C' [2 C D+ (S kron S) D+' C']^-1 C v,  v = vech(S)
*
* causality() has NO effect argument: R2[g,-j] restricts the cause lags in
* every equation EXCEPT the cause's own, so the effect set is all non-cause
* variables and df1 follows from that.  {bf:effect()} narrows it, and is a
* generalisation of the source rather than part of it.
*
* {bf:Note on fidelity.}  vars::causality was read from the package source;
* sandwich::vcovHC, which .grangerGVAR passes by default, is not part of the
* supplied sources.  {bf:vce(oim)} is causality()'s own default,
* Sigma kron (Z'Z)^-1, and is exact.  {bf:vce(robust)} is the standard HC0
* sandwich rather than a port of sandwich::vcovHC, and is labelled as such.
*
* Reporting both the VAR and the VARX test is the point of the exercise: it
* shows whether conditioning on the rest of the world overturns a
* within-country causal reading.

program define _gvar_gc, rclass
    version 14.0

    syntax [,                                   ///
        UNITs(string)                           ///
        CAUSE(string)                           ///
        EFFect(string)                          ///
        LAGS(integer 0)                         ///
        FLAGS(integer 3)                        ///
        VCE(string)                             ///
        noSUMmary                               ///
        SAVing(name)                            ///
        GRaph                                   ///
        NAME(string)                            ///
    ]

    _gvar_require foreign

    if ("`cause'" == "") {
        di as err "cause() is required: the variable(s) whose lags are"
        di as err "excluded under the null, e.g. {bf:cause(r)}"
        exit 198
    }
    if (`flags' < 1) {
        di as err "flags() must be at least 1"
        exit 198
    }

    if ("`vce'" == "") local vce oim
    local vce = lower("`vce'")
    if      ("`vce'" == "oim")    local rb 0
    else if ("`vce'" == "robust") local rb 1
    else {
        di as err "vce() must be {bf:oim} or {bf:robust}"
        exit 198
    }

    mata: st_local("N",  strofreal(gvar_getN()))
    mata: st_local("cn", invtokens(gvar_getcname()'))
    tempname LG
    mata: st_matrix("`LG'", gvar_getlags())

    * ---- which units --------------------------------------------------------
    local usel ""
    if ("`units'" == "") {
        forvalues i = 1/`N' {
            local usel "`usel' `i'"
        }
    }
    else {
        foreach u of local units {
            local p : list posof "`u'" in cn
            if (`p' == 0) {
                di as err "unit {bf:`u'} is not in the model"
                exit 111
            }
            local usel "`usel' `p'"
        }
    }

    * ---- run ----------------------------------------------------------------
    * Stata matrices cannot have zero rows, so J(0, 9, .) is an error rather
    * than an empty accumulator.  Count the eligible units first.
    local elig ""
    local nskip 0
    local skipped ""
    foreach i of local usel {
        local u : word `i' of `cn'
        mata: st_local("yl", invtokens(gvar_getylist(`i')'))
        local ok 1
        local hitc 0
        foreach v of local cause {
            local q : list posof "`v'" in yl
            if (`q' > 0) local hitc 1
        }
        if (`hitc' == 0) local ok 0
        if ("`effect'" != "") {
            local hite 0
            foreach v of local effect {
                local q : list posof "`v'" in yl
                if (`q' > 0) local hite 1
            }
            if (`hite' == 0) local ok 0
        }
        if (`ok') local elig "`elig' `i'"
        else {
            local ++nskip
            local skipped "`skipped' `u'"
        }
    }
    local nelig : word count `elig'
    if (`nelig' == 0) {
        di as err "no unit has {bf:`cause'} in its own endogenous block"
        if ("`effect'" != "") {
            di as err "together with {bf:`effect'}"
        }
        exit 459
    }

    tempname R ROW CS ES
    matrix `R' = J(`nelig', 12, .)
    local rowq 0
    foreach i of local elig {
        local ++rowq
        local u : word `i' of `cn'
        mata: st_local("yl", invtokens(gvar_getylist(`i')'))
        local ki : word count `yl'

        local cpos ""
        foreach v of local cause {
            local q : list posof "`v'" in yl
            if (`q' > 0) local cpos "`cpos' `q'"
        }

        * causality() takes no effect argument: R2[g,-j] restricts the cause
        * lags in EVERY equation except the cause's own, so the effect set is
        * all non-cause variables.  effect() narrows that, and is a
        * generalisation of the source rather than part of it.
        local epos ""
        if ("`effect'" == "") {
            forvalues z = 1/`ki' {
                local in : list posof "`z'" in cpos
                if (`in' == 0) local epos "`epos' `z'"
            }
        }
        else {
            foreach v of local effect {
                local q : list posof "`v'" in yl
                if (`q' > 0) local epos "`epos' `q'"
            }
        }
        local nc : word count `cpos'
        local ne : word count `epos'

        matrix `CS' = J(`nc', 1, 0)
        local q 0
        foreach z of local cpos {
            local ++q
            matrix `CS'[`q', 1] = `z'
        }
        matrix `ES' = J(`ne', 1, 0)
        local q 0
        foreach z of local epos {
            local ++q
            matrix `ES'[`q', 1] = `z'
        }

        local pu = `LG'[`i', 1]
        if (`lags' > 0) local pu = `lags'
        mata: st_matrix("`ROW'", gvar_gcrun(`i', st_matrix("`CS'"), ///
                                            st_matrix("`ES'"), `pu', ///
                                            `flags', `rb'))
        forvalues z = 1/12 {
            matrix `R'[`rowq', `z'] = `ROW'[1, `z']
        }
    }

    local nr = rowsof(`R')

    * ---- counters live outside the summary guard ----------------------------
    local nv 0
    local nx 0
    local nflip 0
    forvalues q = 1/`nr' {
        local pa = `R'[`q', 5]
        local pb = `R'[`q', 9]
        if (`pa' < . & `pa' < 0.05) local ++nv
        if (`pb' < . & `pb' < 0.05) local ++nx
        if (`pa' < . & `pb' < .) {
            if ((`pa' < 0.05) != (`pb' < 0.05)) local ++nflip
        }
    }

    * the label is needed by the graph too, so it cannot live inside the
    * display guard
    local elab "`effect'"
    if ("`effect'" == "") local elab "every other variable"

    if ("`summary'" != "nosummary") {
        _gvar_title "Granger causality: `cause' -> `elab'"
        di as text "  H0: the lags of " as result "`cause'" as text ///
                   " do not enter the " as result "`elab'" as text ///
                   " equation(s)."
        di as text "  Left block: a plain VAR of the unit's own variables."
        di as text "  Right block: the same VAR with the foreign block added"
        di as text "  as exogenous regressors at lags 0 to " ///
                   as result `=`flags'-1' as text "."
        di as text "  Covariance: " as result "`vce'" _continue
        if (`rb' == 1) {
            di as text "  (HC0 sandwich, not a port of sandwich::vcovHC)"
        }
        else {
            di as text "  (Sigma kron (Z'Z)^-1, causality()'s own default)"
        }
        di ""
        di as text "{hline 100}"
        di as text %-10s "  Unit" _col(14) "{bf:VAR}" _col(22) "F" ///
                   _col(34) "df1" _col(42) "df2" _col(52) "p" ///
                   _col(60) "{bf:VARX}" _col(68) "F" _col(78) "p" ///
                   _col(88) "instant p"
        di as text "{hline 100}"
        forvalues q = 1/`nr' {
            local i = `R'[`q', 1]
            local u : word `i' of `cn'
            local pa = `R'[`q', 5]
            local pb = `R'[`q', 9]
            local ma " "
            if (`pa' < . & `pa' < 0.05) local ma "*"
            local mb " "
            if (`pb' < . & `pb' < 0.05) local mb "*"
            local fl " "
            if (`pa' < . & `pb' < .) {
                if ((`pa' < 0.05) != (`pb' < 0.05)) local fl "<"
            }
            di as text "  " %-8s abbrev("`u'", 8)          ///
               _col(16) as result %10.3f `=`R'[`q',2]'     ///
               _col(30) as result %6.0f  `=`R'[`q',3]'     ///
               _col(38) as result %7.0f  `=`R'[`q',4]'     ///
               _col(46) as result %8.3f  `pa' as text "`ma'" ///
               _col(62) as result %10.3f `=`R'[`q',6]'     ///
               _col(72) as result %8.3f  `pb' as text "`mb'" ///
               _col(84) as result %9.3f  `=`R'[`q',12]'      ///
               _col(95) as text "`fl'"
        }
        di as text "{hline 100}"
        di as text "  instant p is the instantaneous-causality test that"
        di as text "  causality() returns alongside the Granger one: a Wald"
        di as text "  chi2 on the contemporaneous covariance block."
        di as text "  * rejects at 5%.  VAR " as result `nv' as text "/" ///
           as result `nr' as text ",  VARX " as result `nx' as text "/" ///
           as result `nr' as text "."
        if (`nflip' > 0) {
            di as text "  < marks the " as result `nflip' as text ///
               " unit(s) where the two disagree: conditioning on"
            di as text "  the foreign block changes the conclusion, which is"
            di as text "  the whole reason for running both."
        }
        else {
            di as text "  The two blocks agree everywhere."
        }
        if (`nskip' > 0) {
            di as text "  Skipped " as result `nskip' as text ///
               " unit(s) lacking one of the variables:"
            di as text "  " as result "`skipped'"
        }
        di ""
    }

    if ("`graph'" != "") {
        _gvar_gc_graph `R' `nr' "`cn'" "`cause'" "`elab'" "`name'"
    }

    if ("`saving'" != "") {
        matrix `saving' = `R'
    }
    return matrix gc = `R', copy
    return scalar nvar  = `nv'
    return scalar nvarx = `nx'
    return scalar nflip = `nflip'
    return scalar nunits = `nr'
    return scalar nskip  = `nskip'
end

* ---------------------------------------------------------------------------
* Granger causality with and without the foreign block, one point per unit.
*
* The interesting thing about a GVAR Granger test is not the level of either
* p-value but whether conditioning on x* changes the answer, so the two are
* plotted against each other.  A unit off the diagonal in the top-left or
* bottom-right quadrant is one whose verdict FLIPS when the foreign variables
* are added -- which for the demo model happens for a quarter of the units, and
* is invisible in either column read on its own.
*
* Both axes are on a log scale because the interesting p-values are small and a
* linear axis piles them all into one corner.
* ---------------------------------------------------------------------------
program define _gvar_gc_graph
    version 14.0
    args R nr cn clab elab name

    _gvar_palette
    local c1   "`r(c1)'"
    local c2   "`r(c2)'"
    local grid "`r(grid)'"
    local reg  "`r(region)'"
    local neg  "`r(neg)'"
    local band "`r(band)'"

    local nm "gvar_gc"
    if ("`name'" != "") local nm "`name'"

    preserve
    quietly {
        clear
        svmat double `R', names(col)
        rename c1 unit
        rename c5 pvar
        rename c9 pvarx

        * floor the p-values so a log axis has something to show for an exact
        * zero, and say so in the note rather than dropping the point
        gen double lp_var  = max(pvar,  1e-6)
        gen double lp_varx = max(pvarx, 1e-6)

        gen str12 ulab = ""
        forvalues q = 1/`nr' {
            local u = unit[`q']
            local l : word `u' of `cn'
            _gvar_ablab "`l'" 10
            replace ulab = "`_ablab'" in `q'
        }

        * a flip is a unit that crosses 5% in one specification but not the
        * other, in either direction
        gen byte flip = 0
        replace flip = 1 if (pvar < 0.05) != (pvarx < 0.05)
        count if flip == 1
        local nflip = r(N)

        gen double f_var  = lp_var  if flip == 1
        gen double f_varx = lp_varx if flip == 1
        gen double k_var  = lp_var  if flip == 0
        gen double k_varx = lp_varx if flip == 0

        twoway (function y = x, range(1e-6 1) lcolor("`grid'")              ///
                    lpattern(solid) lwidth(thin))                          ///
               (scatter k_varx k_var, msymbol(o) msize(small)              ///
                    mcolor("`c1'%70"))                                     ///
               (scatter f_varx f_var, msymbol(O) msize(small)              ///
                    mcolor("`neg'") mlcolor("`c2'") mlwidth(vthin)         ///
                    mlabel(ulab) mlabsize(vsmall) mlabcolor(black)         ///
                    mlabposition(3))                                       ///
               , `reg'                                                     ///
                 xscale(log) yscale(log)                                   ///
                 xlabel(0.000001 "1e-6" 0.0001 "1e-4" 0.01 "0.01"          ///
                        0.05 "0.05" 0.5 "0.5" 1 "1", labsize(vsmall))      ///
                 ylabel(0.000001 "1e-6" 0.0001 "1e-4" 0.01 "0.01"          ///
                        0.05 "0.05" 0.5 "0.5" 1 "1", labsize(vsmall)       ///
                        angle(0))                                          ///
                 xline(0.05, lcolor("`c2'") lpattern(dash) lwidth(thin))   ///
                 yline(0.05, lcolor("`c2'") lpattern(dash) lwidth(thin))   ///
                 xtitle("p-value, VAR (no foreign block)", size(small))    ///
                 ytitle("p-value, VARX* (conditioning on x*)", size(small)) ///
                 title("Does `clab' Granger-cause `elab'?",                ///
                       size(medium) color(black))                          ///
                 subtitle("labelled points change verdict at 5% when x* is added", ///
                          size(small) color(black))                        ///
                 note("dashed lines: 5%; solid line: no change."            ///
                      "  p-values are floored at 1e-6 for the log axis.",   ///
                      size(vsmall) color(black))                            ///
                 legend(off) name(`nm', replace)
    }
    restore

    di as text "  graph saved as {bf:`nm'}  (" as result `nflip' ///
       as text " unit(s) change verdict)"
end