*! _gvar_solve 1.0.1  21aug2026
*! gvar solve -- build the link matrices, stack the country models into the
*! global model, invert to the reduced form, and check stability.
*! Author: Merwan Roudane
*! merwanroudane920@gmail.com
*! https://github.com/merwanroudane
*
* Step -> source map
*   z_it = W_i x_t, link matrices            <- Toolbox create_linkmatrices.m
*   G0 = stack_i A_i W_i, H_l = stack_i B_il W_i,
*   F_l = G0^-1 H_l, eta = G0^-1 zeta        <- Toolbox solve_GVAR.m
*                                               BGVAR gvar_stacking.cpp
*   companion eigenvalues and stability      <- same
*   number of unit eigenvalues vs K - sum r  <- Toolbox gvar.m section 4.3

program define _gvar_solve, rclass
    version 14.0

    syntax [, TOLerance(real 1e-6) ///
              GRaph                ///
              NAME(string)         ///
              NEIGen(integer 20)   ///
              noSUMmary ]

    * "dominant" implies "estimate"; it adds the check that a dominant block
    * DECLARED at setup has actually been fitted, which is what stacking needs.
    _gvar_require dominant

    mata: gvar_solvemodel()

    mata: st_local("N",    strofreal(gvar_getN()))
    mata: st_local("K",    strofreal(gvar_getK()))
    mata: st_local("pmax", strofreal(gvar_getpmax()))

    tempname EIG RK
    mata: st_matrix("`EIG'", gvar_geteig())
    mata: st_matrix("`RK'",  gvar_getrank())

    local ne = rowsof(`EIG')
    local maxmod = `EIG'[1, 1]

    * number of eigenvalues at (numerically) one
    local nunit 0
    forvalues i = 1/`ne' {
        if (abs(`EIG'[`i', 1] - 1) < `tolerance') local ++nunit
    }
    local nexpl 0
    forvalues i = 1/`ne' {
        if (`EIG'[`i', 1] > 1 + `tolerance') local ++nexpl
    }

    * K - sum(r) is the theoretical number of unit roots
    local sumr 0
    forvalues i = 1/`N' {
        local sumr = `sumr' + `RK'[`i', 1]
    }

    * The dominant block is in K but not in RK: RK has one row per COUNTRY
    * model.  Left out, a block that cointegrates internally is counted as gd
    * common trends instead of gd - r_du, and the diagnostic reports a
    * mismatch that is its own arithmetic rather than the model's.
    *
    * Only the VECM branch has a rank to subtract.  For a univariate dominant
    * variable, whether it contributes a unit root depends on where the roots
    * of its own AR polynomial fall -- an estimated property, not an imposed
    * one -- so nothing is subtracted and the report says the count excludes
    * the block.  Assuming 0 or 1 there would trade a visible mismatch for an
    * invisible one.
    local durk 0
    local duk  0
    local dukn 0
    capture mata: st_local("duk", strofreal(gvar_hasdu()))
    if ("`duk'" == "1") {
        mata: st_local("dukn", strofreal(rows(gvar_getduylist())))
        mata: st_local("durk", strofreal(gvar_getdurank()))
        if (`dukn' > 1) local sumr = `sumr' + `durk'
    }
    local expected = `K' - `sumr'

    if ("`summary'" != "nosummary") {
        _gvar_title "Solved GVAR"
        di as text "  Endogenous variables (K)     " as result %8.0f `K'
        di as text "  GVAR lag order               " as result %8.0f `pmax'
        di as text "  Companion eigenvalues        " as result %8.0f `ne'
        di ""
        di as text "  {hline 66}"
        di as text "  Stability"
        di as text "  {hline 66}"
        di as text "  Largest eigenvalue modulus   " as result %12.6f `maxmod'
        di as text "  Eigenvalues > 1 (tol " %6.0e `tolerance' ")     " ///
                   as result %8.0f `nexpl'
        if (`nexpl' == 0) {
            di as text "  " as result "The GVAR is stable: no eigenvalue exceeds unity."
        }
        else {
            * {err:} markup, not "as err".  This is a WARNING inside the summary
            * block, and "as err" prints through -quietly- -- so a loop doing
            * "qui gvar solve" got instability warnings it had asked to suppress.
            * The error paths in this file are different: they exit, so they
            * should print regardless.
            di as text "  {err:The GVAR is UNSTABLE: `nexpl' eigenvalue(s)" ///
                       " exceed unity.}"
            di as text "  {err:Impulse responses will not die out. Re-check" ///
                       " the lag}"
            di as text "  {err:orders, the ranks and the deterministic case.}"
        }
        di ""
        di as text "  {hline 66}"
        di as text "  Unit roots"
        di as text "  {hline 66}"
        di as text "  Eigenvalues at unity         " as result %8.0f `nunit'
        * A Bayesian VARX in levels imposes NO cointegrating rank, so K - sum(r)
        * has nothing to compare against: m.rnk holds whatever gvar coint or an
        * earlier gvar estimate left behind, and printing a "mismatch" against
        * it would be this command's arithmetic rather than the model's.  BGVAR
        * runs no such check either -- the eigenvalue trim in gvar bayes is what
        * takes its place.
        mata: st_local("et", gvar_getesttype())
        if ("`et'" == "bvarx") {
            di as text "  " as result "The model was fit by gvar bayes as a" ///
               " VARX in levels,"
            di as text "  so no cointegrating rank is imposed and there is no"
            di as text "  K - sum(r) to compare this against.  Stability was"
            di as text "  screened draw by draw instead; see the eigenvalue"
            di as text "  trim reported by {bf:gvar bayes}."
        }
        else {
            di as text "  K - sum(r) expected          " ///
               as result %8.0f `expected'
            if ("`duk'" == "1") {
                if (`dukn' > 1) {
                    di as text "  including the dominant block's rank " ///
                       as result %5.0f `durk'
                }
                else {
                    di as text "  {err:The dominant block is univariate, so" ///
                               " no rank is subtracted}"
                    di as text "  {err:for it; the expected count excludes it" ///
                               " and a difference of one}"
                    di as text "  {err:is what a unit root in that block" ///
                               " looks like.}"
                }
            }
            if (`nunit' == `expected') {
                di as text "  " as result ///
                   "Consistent with the imposed cointegrating ranks."
            }
            else {
                di as text "  {err:These differ; the number of common" ///
                           " stochastic}"
                di as text "  {err:trends implied by the solved model is not" ///
                           " the one implied}"
                di as text "  {err:by the country-by-country ranks. This is" ///
                           " common when the}"
                di as text "  {err:ranks are chosen unit by unit} -- see" ///
                           " {help gvar_methods}."
            }
        }
        di ""
        di as text "  {hline 66}"
        di as text "  Largest " as result `neigen' as text " eigenvalue moduli"
        di as text "  {hline 66}"
        local show = min(`neigen', `ne')
        local line "   "
        forvalues i = 1/`show' {
            local piece : display %8.5f `EIG'[`i', 1]
            if (length("`line'`piece'  ") > 76) {
                di as text "`line'"
                local line "   "
            }
            local line "`line'`piece'  "
        }
        if (trim("`line'") != "") di as text "`line'"
        di as text "  {hline 66}"
        di ""
    }

    if ("`graph'" != "") {
        _gvar_solve_graph `EIG' "`name'"
    }

    return matrix eigen = `EIG', copy
    return scalar maxmod   = `maxmod'
    return scalar nunit    = `nunit'
    return scalar nexpl    = `nexpl'
    return scalar expected = `expected'
    return scalar K        = `K'
    return scalar pmax     = `pmax'
    return scalar stable   = (`nexpl' == 0)
end

* ---------------------------------------------------------------------------
* Eigenvalue moduli against the unit circle
* ---------------------------------------------------------------------------
program define _gvar_solve_graph
    version 14.0
    args EIG gname
    if ("`gname'" == "") local gname gvar_eigen

    _gvar_palette
    local reg  "`r(region)'"
    local c1   "`r(c1)'"
    local zero "`r(zero)'"

    preserve
    clear
    local ne = rowsof(`EIG')
    qui set obs `ne'
    qui gen double modulus = .
    forvalues i = 1/`ne' {
        qui replace modulus = `EIG'[`i', 1] in `i'
    }
    qui gen int idx = _n

    twoway ///
        (scatter modulus idx, msymbol(circle) msize(small) ///
             mcolor("`c1'%70") mlcolor("`c1'") mlwidth(vthin)) ///
        , `reg' ///
          yline(1, lcolor("`zero'") lpattern(dash) lwidth(medthin)) ///
          ylabel(, angle(0) labsize(small) grid glcolor(gs15)) ///
          xlabel(, labsize(small)) ///
          ytitle("modulus", size(small)) ///
          xtitle("eigenvalue (sorted)", size(small)) ///
          title("Eigenvalues of the GVAR companion matrix", ///
                size(medium) color(black)) ///
          subtitle("moduli must not exceed the unit circle", ///
                   size(vsmall) color(gs7)) ///
          name(`gname', replace)
    restore
end