*! tnardll version 1.0.0  03jun2026
*! Threshold (Nonlinear) Autoregressive Distributed Lag model -- TARDL
*!
*! Implements the Threshold ARDL of Cho, Greenwood-Nimmo & Shin (2020c,d),
*! surveyed in Cho, Greenwood-Nimmo & Shin (2021, "Recent Developments of the
*! Autoregressive Distributed Lag Modelling Framework"). The threshold/momentum
*! decomposition is applied to the FIRST differences of one regressor around one
*! or more UNKNOWN thresholds (tau), estimated by profile (concentrated) least
*! squares.  Number of regimes S in {1,2,3} chosen by information criteria
*! (six criteria: AIC/SIC/HQIC and the Pitarakis 2006 modified pAIC/pSIC/pHQIC)
*! or fixed by the user.  Threshold existence is tested by a quasi-likelihood
*! ratio (QLR) statistic with a fixed-regressor wild bootstrap p-value
*! (the asymptotic null law is non-pivotal: Davies 1977/87 problem).
*!
*! Author: Dr Merwan Roudane <merwanroudane920@gmail.com>
*! https://github.com/merwanroudane

program define tnardll, eclass sortpreserve
    version 17.0

    if replay() {
        if "`e(cmd)'" != "tnardll" {
            error 301
        }
        Display `0'
        exit
    }

    syntax varlist(min=2 numeric ts)  [if] [in] [ ,    ///
            THReshold(varname numeric)                 ///
            Lags(numlist integer >=1 max=2)            ///
            REGimes(integer 0)                         ///
            MAXReg(integer 3)                          ///
            IC(string)                                 ///
            TRim(real 0.15)                            ///
            GRid(integer 50)                           ///
            QLR                                        ///
            Breps(integer 499)                         ///
            SEED(string)                               ///
            noCONStant                                 ///
            TRend                                      ///
            Robust                                     ///
            Level(cilevel)                             ///
            noDOTS ]

    * ---------------------------------------------------------------
    * Parse variable lists
    * ---------------------------------------------------------------
    gettoken depvar indeps : varlist
    local indeps : list retokenize indeps
    if "`indeps'" == "" {
        di as err "at least one regressor is required"
        exit 102
    }

    * threshold (decomposed) variable
    if "`threshold'" == "" {
        gettoken thrvar : indeps
    }
    else {
        local thrvar "`threshold'"
        local ck : list thrvar in indeps
        if !`ck' {
            di as err "threshold() variable {bf:`thrvar'} must be among the regressors"
            exit 198
        }
    }
    * other (linear) regressors
    local othervars : list indeps - thrvar

    * ---------------------------------------------------------------
    * Lags
    * ---------------------------------------------------------------
    if "`lags'" == "" {
        local p 1
        local q 1
    }
    else {
        tokenize `lags'
        local p `1'
        local q = cond("`2'"=="", "`1'", "`2'")
    }

    * ---------------------------------------------------------------
    * Information criterion
    * ---------------------------------------------------------------
    if "`ic'" == "" local ic "sic"
    local ic = lower("`ic'")
    local okic "aic sic hqic paic psic phqic"
    if !`: list ic in okic' {
        di as err "ic() must be one of: `okic'"
        exit 198
    }

    * regimes / maxreg bounds (S in {1,2,3})
    if `maxreg' > 3 {
        di as txt "note: maximum supported regimes is 3 (two thresholds); maxreg() set to 3"
        local maxreg 3
    }
    if `maxreg' < 1 local maxreg 1
    if `regimes' != 0 & (`regimes' < 1 | `regimes' > 3) {
        di as err "regimes() must be 1, 2 or 3 (0 = choose by information criterion)"
        exit 198
    }
    if `trim' <= 0 | `trim' >= 0.5 {
        di as err "trim() must lie strictly between 0 and 0.5"
        exit 198
    }
    if `grid' < 5 local grid 5

    * vce
    local vce = cond("`robust'"!="", "robust", "ols")

    * constant / trend flags
    local consflag = cond("`constant'"=="", 1, 0)
    local trflag   = cond("`trend'"!="",   1, 0)

    * seed
    if "`seed'" != "" {
        set seed `seed'
    }

    * ---------------------------------------------------------------
    * tsset & sample
    * ---------------------------------------------------------------
    qui tsset
    if "`r(panelvar)'" != "" {
        di as err "tnardll requires a pure time series (use tsset {it:timevar}); panel data not supported"
        exit 198
    }
    local timevar "`r(timevar)'"
    local tdelta = cond("`r(tdelta)'"=="", 1, `r(tdelta)')

    marksample touse
    markout `touse' `depvar' `indeps'
    qui count if `touse'
    if r(N) == 0 {
        di as err "no observations"
        exit 2000
    }

    tempvar used
    qui gen byte `used' = 0

    * ---------------------------------------------------------------
    * Hand off to Mata engine
    * ---------------------------------------------------------------
    mata: _tnardll_run()

    * scalars / locals returned via st_* in Mata:
    *   __err (string, nonempty => abort)
    if "`__err'" != "" {
        di as err "`__err'"
        exit 459
    }

    * ---------------------------------------------------------------
    * Build coefficient names (order MUST match Mata design column order)
    * ---------------------------------------------------------------
    local S = __S
    local m : word count `othervars'
    local names ""
    * 1. error-correction term (rho on y_{t-1})
    local names "`names' ec"
    * 2. long-run levels: threshold regimes then other vars
    forvalues s = 1/`S' {
        local names "`names' lr_`thrvar'_r`s'"
    }
    foreach w of local othervars {
        local names "`names' lr_`w'"
    }
    * 3. short-run Dy lags 1..p-1
    forvalues j = 1/`=`p'-1' {
        local names "`names' D`j'.`depvar'"
    }
    * 4. short-run Dx regime lags 0..q-1
    forvalues s = 1/`S' {
        forvalues j = 0/`=`q'-1' {
            local names "`names' D`j'.`thrvar'_r`s'"
        }
    }
    * 5. short-run other-var diffs 0..q-1
    foreach w of local othervars {
        forvalues j = 0/`=`q'-1' {
            local names "`names' D`j'.`w'"
        }
    }
    * 6. trend
    if `trflag' local names "`names' trend"
    * 7. constant
    if `consflag' local names "`names' _cons"

    * sanitise names to valid Stata coefficient names
    local cnames ""
    foreach nm of local names {
        local cnames "`cnames' `=strtoname("`nm'")'"
    }

    matrix colnames __b = `cnames'
    matrix colnames __V = `cnames'
    matrix rownames __V = `cnames'

    * long-run coefficient names (returned by Mata in __lrnames)
    local lrc ""
    foreach nm of local __lrnames {
        local lrc "`lrc' `=strtoname("`nm'")'"
    }
    matrix colnames __lrb = `lrc'
    matrix colnames __lrV = `lrc'
    matrix rownames __lrV = `lrc'

    * ---------------------------------------------------------------
    * Post results
    * ---------------------------------------------------------------
    ereturn post __b __V , esample(`used') depname(`depvar') obs(`=__N')

    ereturn local cmd       "tnardll"
    ereturn local cmdline   "tnardll `0'"
    ereturn local title     "Threshold (Nonlinear) ARDL -- TARDL"
    ereturn local depvar    "`depvar'"
    ereturn local indepvars "`indeps'"
    ereturn local thrvar    "`thrvar'"
    ereturn local othervars "`othervars'"
    ereturn local timevar   "`timevar'"
    ereturn local ic        "`ic'"
    ereturn local vce       "`vce'"
    ereturn local properties "b V"
    ereturn local predict   "tnardll_p"
    ereturn local marginsnotok "_ALL"

    ereturn scalar N        = __N
    ereturn scalar k        = __k
    ereturn scalar df_r     = __df_r
    ereturn scalar S        = __S
    ereturn scalar nthr     = __S - 1
    ereturn scalar p        = `p'
    ereturn scalar q        = `q'
    ereturn scalar rho      = __rho
    ereturn scalar ssr      = __ssr
    ereturn scalar rmse     = __rmse
    ereturn scalar r2       = __r2
    ereturn scalar r2_a     = __r2_a
    ereturn scalar trim     = `trim'
    ereturn scalar grid     = `grid'
    ereturn scalar level    = `level'
    ereturn scalar cons     = `consflag'
    ereturn scalar trend    = `trflag'
    ereturn scalar pss_k    = __pssk
    ereturn scalar Fpss     = __Fpss
    ereturn scalar tBDM     = __tBDM
    ereturn scalar lr_asym_chi2 = __lachi2
    ereturn scalar lr_asym_df   = __ladf
    ereturn scalar lr_asym_p    = __lap
    ereturn scalar sr_asym_chi2 = __sachi2
    ereturn scalar sr_asym_df   = __sadf
    ereturn scalar sr_asym_p    = __sap

    if "`qlr'" != "" {
        ereturn scalar qlr   = __qlr
        ereturn scalar qlr_p = __qlrp
        ereturn scalar breps = `breps'
    }

    if `S' >= 2 {
        matrix colnames __thr = tau
        ereturn matrix thresholds = __thr
    }
    ereturn matrix lr_b   = __lrb
    ereturn matrix lr_V   = __lrV
    ereturn matrix ictable = __ictab
    capture confirm matrix __profile
    if !_rc {
        ereturn matrix profile = __profile
    }
    * structural pieces for multipliers
    ereturn matrix theta = __theta
    capture confirm matrix __phi
    if !_rc  ereturn matrix phi = __phi
    ereturn matrix pimat = __pimat

    * ---------------------------------------------------------------
    * PSS bounds-test critical values via Kripfganz & Schneider (2020)
    * surface regression (ships with -ardl-).  Mapped to a determ. case:
    *   no constant -> 1 ; unrestricted constant -> 3 ; +trend -> 5.
    * k = number of long-run forcing terms (S regimes + other levels);
    * bounds are an approximate reference for the threshold decomposition.
    * ---------------------------------------------------------------
    local pcase = cond(`consflag'==0, 1, cond(`trflag'==1, 5, 3))
    ereturn scalar pss_case = `pcase'
    local pk = e(pss_k)
    capture qui ardlbounds , case(`pcase') stat(f) k(`pk') siglevels(10 5 1)
    if !_rc {
        tempname Fcv tcv
        matrix `Fcv' = r(cvmat)
        ereturn matrix F_critval = `Fcv'
        capture qui ardlbounds , case(`pcase') stat(t) k(`pk') siglevels(10 5 1)
        if !_rc {
            matrix `tcv' = r(cvmat)
            ereturn matrix t_critval = `tcv'
        }
    }

    Display , level(`level')
end

* =====================================================================
* Display
* =====================================================================
program define Display
    syntax [, Level(cilevel) ]
    if "`level'" == "" local level = e(level)

    local depvar  "`e(depvar)'"
    local thrvar  "`e(thrvar)'"
    local S       = e(S)
    local p       = e(p)
    local q       = e(q)

    di ""
    di as txt "{hline 78}"
    di as txt "Threshold (Nonlinear) ARDL                       " ///
              _col(50) "Number of obs " _col(66) "= " as res %9.0f e(N)
    di as txt "Dependent variable: " as res "`depvar'" ///
              _col(50) as txt "Regimes (S)   " _col(66) "= " as res %9.0f e(S)
    di as txt "Threshold variable: " as res "D.`thrvar'" ///
              _col(50) as txt "ARDL(p,q)     " _col(66) "= " as res "(`p',`q')"
    di as txt "Adjustment (rho)  : " as res %9.4f e(rho) ///
              _col(50) as txt "Adj. R-squared" _col(66) "= " as res %9.4f e(r2_a)

    * thresholds
    local nthr = e(nthr)
    if `nthr' > 0 {
        matrix b = e(thresholds)
        local tlist ""
        forvalues i = 1/`nthr' {
            local tv = b[`i',1]
            local tlist "`tlist' `:di %7.4f `tv''"
        }
        di as txt "Estimated threshold(s) in D.`thrvar':" as res "`tlist'"
    }
    else {
        di as txt "Estimated threshold(s): " as res "(none -- linear ARDL)"
    }
    di as txt "{hline 78}"

    * ---------------- Long-run table -----------------
    di ""
    di as txt "Long-run coefficients (beta = -theta/rho), delta-method SE"
    _coltab_header , level(`level')

    matrix lb = e(lr_b)
    matrix lV = e(lr_V)
    local lrnames : colnames lb
    local zc = invnormal(1 - (100-`level')/200)
    local i 0
    foreach nm of local lrnames {
        local ++i
        local est = lb[1,`i']
        local se  = sqrt(lV[`i',`i'])
        _coltab_row , name("`nm'") b(`est') se(`se') z(`zc')
    }
    di as txt "{hline 78}"

    * ---------------- Short-run / ECM table -----------------
    di ""
    di as txt "Short-run dynamics and error-correction (unrestricted ECM)"
    _coltab_header , level(`level')
    matrix bb = e(b)
    matrix VV = e(V)
    local bn : colnames bb
    local i 0
    foreach nm of local bn {
        local ++i
        local est = bb[1,`i']
        local se  = sqrt(VV[`i',`i'])
        _coltab_row , name("`nm'") b(`est') se(`se') z(`zc')
    }
    di as txt "{hline 78}"

    * ---------------- Tests -----------------
    di ""
    di as txt "Asymmetry / threshold tests"
    di as txt "{hline 78}"
    di as txt "  Long-run regime equality (H0: beta^(s) all equal)"
    di as txt "      chi2(" as res e(lr_asym_df) as txt ") = " ///
        as res %8.3f e(lr_asym_chi2) as txt "   Prob > chi2 = " as res %6.4f e(lr_asym_p)
    di as txt "  Short-run regime equality (H0: sum pi^(s) all equal)"
    di as txt "      chi2(" as res e(sr_asym_df) as txt ") = " ///
        as res %8.3f e(sr_asym_chi2) as txt "   Prob > chi2 = " as res %6.4f e(sr_asym_p)
    if !missing(e(qlr)) {
        di as txt "  QLR threshold-existence test (H0: linear ARDL vs S=2 TARDL)"
        di as txt "      QLR = " as res %8.3f e(qlr) ///
            as txt "   bootstrap p = " as res %6.4f e(qlr_p) ///
            as txt "   (B = " as res e(breps) as txt ")"
    }
    local pk = e(pss_k)
    local pcase = e(pss_case)
    di as txt "  PSS bounds test of no levels relationship (H0: no cointegration)"
    di as txt "      F = " as res %8.3f e(Fpss) ///
        as txt "    t(rho) = " as res %8.3f e(tBDM) ///
        as txt "    [k = " as res `pk' as txt ", case " as res `pcase' as txt "]"
    capture confirm matrix e(F_critval)
    if !_rc {
        matrix Fcv = e(F_critval)
        di as txt "      Kripfganz-Schneider (2020) asymptotic critical bounds:"
        di as txt "                  10%               5%                1%"
        di as txt "                I(0)    I(1)     I(0)    I(1)     I(0)    I(1)"
        di as txt "        F   " as res %8.3f Fcv[1,1] %8.3f Fcv[1,2] ///
            "  " %8.3f Fcv[1,3] %8.3f Fcv[1,4] ///
            "  " %8.3f Fcv[1,5] %8.3f Fcv[1,6]
        capture confirm matrix e(t_critval)
        if !_rc {
            matrix tcv = e(t_critval)
            di as txt "        t   " as res %8.3f tcv[1,1] %8.3f tcv[1,2] ///
                "  " %8.3f tcv[1,3] %8.3f tcv[1,4] ///
                "  " %8.3f tcv[1,5] %8.3f tcv[1,6]
        }
        local fi0 = Fcv[1,3]
        local fi1 = Fcv[1,4]
        local Fv  = e(Fpss)
        if `Fv' > `fi1'      local fdec "F > I(1) => reject H0: a levels relationship exists"
        else if `Fv' < `fi0' local fdec "F < I(0) => do not reject H0: no levels relationship"
        else                 local fdec "I(0) <= F <= I(1) => inconclusive"
        di as txt "      Decision (5%): " as res "`fdec'"
    }
    else {
        di as txt "      compare F,t to Pesaran/Shin/Smith (2001) I(0)/I(1) bounds"
        di as txt "      (install {bf:ardl} (Kripfganz/Schneider) for automatic bounds)"
    }
    di as txt "{hline 78}"

    * ---------------- IC table -----------------
    matrix ict = e(ictable)
    di ""
    di as txt "Regime selection -- information criteria (chosen S minimises {bf:`e(ic)'})"
    di as txt "{hline 78}"
    di as txt "    S  " _col(12) "AIC" _col(24) "SIC" _col(36) "HQIC" ///
        _col(48) "pAIC" _col(60) "pSIC" _col(72) "pHQIC"
    local nr = rowsof(ict)
    forvalues r = 1/`nr' {
        local sval = ict[`r',1]
        di as res %5.0f `sval' as txt "  " ///
            _col(8)  as res %10.4f ict[`r',2] ///
            _col(20) as res %10.4f ict[`r',3] ///
            _col(32) as res %10.4f ict[`r',4] ///
            _col(44) as res %10.4f ict[`r',5] ///
            _col(56) as res %10.4f ict[`r',6] ///
            _col(68) as res %10.4f ict[`r',7]
    }
    di as txt "{hline 78}"
    di as txt "Note: threshold p-value bootstrap is fixed-regressor wild (Rademacher);"
    di as txt "      see {help tnardll##refs:references}. Long-run SEs are FM-type asymptotic."
end

* small helpers for coefficient tables
program define _coltab_header
    syntax , Level(cilevel)
    di as txt "{hline 78}"
    di as txt %-22s "" _col(23) "{c |}" ///
        _col(27) "Coef." _col(39) "Std.err." _col(51) "z" _col(58) "P>|z|" ///
        _col(66) "[`level'% c.i.]"
    di as txt "{hline 23}{c +}{hline 54}"
end

program define _coltab_row
    syntax , name(string) b(real) se(real) z(real)
    if `se' <= 0 | missing(`se') {
        di as txt %-22s abbrev("`name'",22) _col(23) "{c |}" ///
            _col(25) as res %10.4f `b' _col(37) as txt "  (dropped/exact)"
        exit
    }
    local zst = `b'/`se'
    local pv  = 2*normal(-abs(`zst'))
    local lo  = `b' - `z'*`se'
    local hi  = `b' + `z'*`se'
    di as txt %-22s abbrev("`name'",22) _col(23) "{c |}" ///
        _col(25) as res %10.4f `b' ///
        _col(37) as res %10.4f `se' ///
        _col(49) as res %7.2f `zst' ///
        _col(57) as res %6.3f `pv' ///
        _col(64) as res %9.4f `lo' " " %9.4f `hi'
end

* =====================================================================
* Mata engine
* =====================================================================
version 17.0
mata:
mata set matastrict off

// ---- helpers ---------------------------------------------------------
real colvector _tn_lag(real colvector v, real scalar k)
{
    real scalar n
    n = rows(v)
    if (k <= 0) return(v)
    if (k >= n) return(J(n,1,.))
    return( J(k,1,.) \ v[|1 \ n-k|] )
}

real colvector _tn_cumsum(real colvector v)
{
    real scalar i, n
    real colvector c
    n = rows(v)
    c = J(n,1,0)
    if (n == 0) return(c)
    c[1] = (v[1]==. ? 0 : v[1])
    for (i=2; i<=n; i++) c[i] = c[i-1] + (v[i]==. ? 0 : v[i])
    return(c)
}

// segment the difference series dx into S regime-columns given thresholds tau
real matrix _tn_seg(real colvector dx, real rowvector tau)
{
    real scalar n, S, t, s
    real matrix seg
    n = rows(dx)
    S = cols(tau) + 1
    seg = J(n, S, 0)
    for (t=1; t<=n; t++) {
        if (dx[t]==.) continue
        if (S==1) {
            seg[t,1] = dx[t]
            continue
        }
        if (dx[t] <= tau[1]) {
            seg[t,1] = dx[t]
        }
        else if (dx[t] > tau[S-1]) {
            seg[t,S] = dx[t]
        }
        else {
            for (s=2; s<=S-1; s++) {
                if (dx[t] > tau[s-1] && dx[t] <= tau[s]) seg[t,s] = dx[t]
            }
        }
    }
    return(seg)
}

// build design; returns dependent (dyv) and regressors (X) over usable rows,
// plus column-index bookkeeping and the surviving row indices (keep)
struct tndes {
    real colvector dyv
    real matrix    X
    real colvector keep
    real scalar    S, m, p, q, irho
    real rowvector ilr_xT, ilr_w, iphi
    real matrix    ipi
}

struct tndes scalar _tn_design(real colvector y, real colvector xT,
        real matrix W, real scalar p, real scalar q,
        real rowvector tau, real scalar consf, real scalar trf)
{
    struct tndes scalar d
    real scalar n, S, m, s, j, jj, col
    real colvector dy, dx, dw
    real matrix seg, xseg, X, M
    real rowvector pirow

    n = rows(y)
    S = cols(tau) + 1
    m = cols(W)
    dy = y - _tn_lag(y,1)
    dx = xT - _tn_lag(xT,1)
    seg  = _tn_seg(dx, tau)
    xseg = J(n,S,0)
    for (s=1; s<=S; s++) xseg[.,s] = _tn_cumsum(seg[.,s])

    X = J(n,0,.)
    col = 0
    // 1. rho : y_{t-1}
    X = _tn_lag(y,1)
    col = 1
    d.irho = 1
    // 2. long-run levels: threshold regimes
    d.ilr_xT = J(1,S,0)
    for (s=1; s<=S; s++) {
        X = X, _tn_lag(xseg[.,s],1)
        col++
        d.ilr_xT[s] = col
    }
    // long-run levels: other vars
    d.ilr_w = (m>0 ? J(1,m,0) : J(1,0,0))
    for (j=1; j<=m; j++) {
        X = X, _tn_lag(W[.,j],1)
        col++
        d.ilr_w[j] = col
    }
    // 3. short-run Dy lags 1..p-1
    d.iphi = (p>1 ? J(1,p-1,0) : J(1,0,0))
    for (j=1; j<=p-1; j++) {
        X = X, _tn_lag(dy,j)
        col++
        d.iphi[j] = col
    }
    // 4. short-run Dx regime lags 0..q-1
    d.ipi = J(S,q,0)
    for (s=1; s<=S; s++) {
        for (j=0; j<=q-1; j++) {
            X = X, _tn_lag(seg[.,s],j)
            col++
            d.ipi[s,j+1] = col
        }
    }
    // 5. short-run other-var diffs 0..q-1
    for (jj=1; jj<=m; jj++) {
        dw = W[.,jj] - _tn_lag(W[.,jj],1)
        for (j=0; j<=q-1; j++) {
            X = X, _tn_lag(dw,j)
            col++
        }
    }
    // 6. trend
    if (trf) {
        X = X, (1::n)
        col++
    }
    // 7. constant
    if (consf) {
        X = X, J(n,1,1)
        col++
    }

    M = dy, X
    d.keep = selectindex(rowmissing(M):==0)
    d.dyv = dy[d.keep]
    d.X   = X[d.keep, .]
    d.S = S; d.m = m; d.p = p; d.q = q
    return(d)
}

// OLS fit: returns ssr; optionally fills b, V
real scalar _tn_ols(real colvector yv, real matrix X, real scalar robust,
                    real colvector b, real matrix V)
{
    real scalar n, k, ssr, s2
    real colvector e
    real matrix XtXi, meat
    n = rows(X); k = cols(X)
    XtXi = invsym(cross(X,X))
    b = XtXi * cross(X, yv)
    e = yv - X*b
    ssr = cross(e,e)
    if (n > k) s2 = ssr/(n-k)
    else s2 = ssr/n
    if (robust) {
        meat = cross(X, e:^2, X)
        V = (n/(n-k)) * XtXi * meat * XtXi
    }
    else {
        V = s2 * XtXi
    }
    return(ssr)
}

// candidate thresholds from the difference series
real colvector _tn_cand(real colvector dx, real scalar trim, real scalar G)
{
    real colvector v, sv, sub, idx
    real scalar nd, lo, hi
    v = dx[selectindex(dx:!=.)]
    sv = sort(v,1)
    nd = rows(sv)
    lo = floor(trim*nd)+1
    hi = ceil((1-trim)*nd)
    if (lo < 1) lo = 1
    if (hi > nd) hi = nd
    if (hi < lo) {
        lo = 1
        hi = nd
    }
    sub = sv[|lo \ hi|]
    sub = uniqrows(sub)
    if (rows(sub) > G) {
        idx = round(rangen(1, rows(sub), G))
        sub = sub[idx]
        sub = uniqrows(sub)
    }
    return(sub)
}

// profile least squares for a given number of regimes S
// returns min ssr; fills besttau (rowvector length S-1) and (optionally) profile
real scalar _tn_profile(real colvector y, real colvector xT, real matrix W,
        real scalar p, real scalar q, real scalar S,
        real colvector cand, real scalar minobs,
        real scalar consf, real scalar trf,
        rowvector besttau, real matrix profile)
{
    real scalar best, ng, i, j, ssr, nb, na, nm, ok
    real rowvector tau
    real colvector dx, b
    struct tndes scalar d
    real matrix V

    best = .
    besttau = J(1, S-1, .)
    ng = rows(cand)
    dx = xT - _tn_lag(xT,1)

    if (S == 1) {
        d = _tn_design(y, xT, W, p, q, J(1,0,.), consf, trf)
        best = _tn_ols(d.dyv, d.X, 0, b, V)
        profile = J(0,0,.)
        return(best)
    }
    else if (S == 2) {
        profile = J(ng, 2, .)
        for (i=1; i<=ng; i++) {
            tau = cand[i]
            nb = sum(dx :<= tau)
            na = sum(dx :> tau)
            profile[i,1] = cand[i]
            if (nb < minobs | na < minobs) continue
            d = _tn_design(y, xT, W, p, q, tau, consf, trf)
            ssr = _tn_ols(d.dyv, d.X, 0, b, V)
            profile[i,2] = ssr
            if (ssr < best) {
                best = ssr
                besttau = tau
            }
        }
        return(best)
    }
    else {
        // S == 3 : two thresholds tau1 < tau2
        profile = J(0,0,.)
        for (i=1; i<=ng-1; i++) {
            for (j=i+1; j<=ng; j++) {
                tau = (cand[i], cand[j])
                nb = sum(dx :<= tau[1])
                nm = sum((dx :> tau[1]) :& (dx :<= tau[2]))
                na = sum(dx :> tau[2])
                if (nb < minobs | nm < minobs | na < minobs) continue
                d = _tn_design(y, xT, W, p, q, tau, consf, trf)
                ssr = _tn_ols(d.dyv, d.X, 0, b, V)
                if (ssr < best) {
                    best = ssr
                    besttau = tau
                }
            }
        }
        return(best)
    }
}

// QLR statistic + fixed-regressor wild bootstrap p-value (S=1 null vs S=2 alt)
void _tn_qlr(real colvector y, real colvector xT, real matrix W,
        real scalar p, real scalar q, real colvector cand, real scalar minobs,
        real scalar consf, real scalar trf, real scalar B,
        real scalar nN, real scalar qlr_obs, real scalar pval, real scalar dots)
{
    struct tndes scalar d0, d
    real scalar ng, i, bcount, bi, ssr0, ssrT, ssr, nb, na
    real scalar ssr0b, ssrTb, qlrb
    real colvector b0, fit0, e0, dyb, bb, dx
    real matrix V, X0
    pointer(real matrix) rowvector Xs
    pointer(real matrix) rowvector Hs
    real rowvector v

    dx = xT - _tn_lag(xT,1)
    // null design (linear ARDL)
    d0 = _tn_design(y, xT, W, p, q, J(1,0,.), consf, trf)
    X0 = d0.X
    ssr0 = _tn_ols(d0.dyv, X0, 0, b0, V)
    fit0 = X0 * b0
    e0   = d0.dyv - fit0
    nN   = rows(d0.dyv)

    // alternative: profile over candidates, store admissible designs
    ng = rows(cand)
    ssrT = .
    Xs = J(1,0, NULL)
    Hs = J(1,0, NULL)
    for (i=1; i<=ng; i++) {
        nb = sum(dx :<= cand[i])
        na = sum(dx :>  cand[i])
        if (nb < minobs | na < minobs) continue
        d = _tn_design(y, xT, W, p, q, cand[i], consf, trf)
        if (rows(d.dyv) != nN) continue   // alignment guard
        ssr = _tn_ols(d.dyv, d.X, 0, bb, V)
        if (ssr < ssrT) ssrT = ssr
        Xs = Xs, &(d.X)
        Hs = Hs, &(invsym(cross(d.X,d.X)) * d.X')
    }
    qlr_obs = nN * (1 - ssrT/ssr0)

    // bootstrap
    real scalar H0pre
    real matrix H0
    H0 = invsym(cross(X0,X0)) * X0'
    bcount = 0
    real scalar nstore
    nstore = cols(Xs)
    for (bi=1; bi<=B; bi++) {
        v = (runiform(1,nN) :> 0.5) :* 2 :- 1     // Rademacher +/-1
        dyb = fit0 + e0 :* v'
        // null
        ssr0b = quadcross(dyb,dyb) - quadcross(dyb, X0*(H0*dyb))
        // alternative: min over stored designs
        ssrTb = .
        for (i=1; i<=nstore; i++) {
            ssr = quadcross(dyb,dyb) - quadcross(dyb, (*Xs[i])*((*Hs[i])*dyb))
            if (ssr < ssrTb) ssrTb = ssr
        }
        qlrb = nN * (1 - ssrTb/ssr0b)
        if (qlrb >= qlr_obs) bcount++
        if (dots) {
            if (mod(bi,50)==0) printf(".")
            displayflush()
        }
    }
    if (dots) printf("\n")
    pval = (bcount+1)/(B+1)
}

// =====================================================================
// main driver -- reads options & data from Stata, writes results back
// =====================================================================
void _tnardll_run()
{
    real scalar p, q, S, fixedS, maxreg, trim, G, consf, trf, robust
    real scalar doqlr, B, dots, minobs, level
    real scalar i, s, j, nN, k, ssr, sigma2, rmse, tss, r2, r2a, nthr
    real scalar ssr1, ssrS, qlr_obs, qlr_p
    string scalar depvar, thrvar, othervars, touse, ic
    real colvector y, xT, tv, b, cand, rowsused, keepflag, lrb_v, dx
    real matrix W, X, V, Wdummy, profile, ictab, lrV, thetaM, piM, phiM
    struct tndes scalar d
    real rowvector besttau, tauS, lrnames_idx
    real scalar tdelta, ndiff

    // ---- read options ----
    depvar    = st_local("depvar")
    thrvar    = st_local("thrvar")
    othervars = st_local("othervars")
    touse     = st_local("used")               // we will fill esample here
    p   = strtoreal(st_local("p"))
    q   = strtoreal(st_local("q"))
    fixedS  = strtoreal(st_local("regimes"))
    maxreg  = strtoreal(st_local("maxreg"))
    trim    = strtoreal(st_local("trim"))
    G       = strtoreal(st_local("grid"))
    consf   = strtoreal(st_local("consflag"))
    trf     = strtoreal(st_local("trflag"))
    robust  = (st_local("vce")=="robust")
    doqlr   = (st_local("qlr")!="")
    B       = strtoreal(st_local("breps"))
    dots    = (st_local("dots")!="nodots")
    level   = strtoreal(st_local("level"))
    ic      = st_local("ic")

    // ---- pull data over the marked sample (already time-ordered) ----
    real colvector tvar0
    tvar0 = st_data(., st_local("touse"))           // selection vector full length
    rowsused = selectindex(tvar0)                   // dataset rows in the sample

    y  = st_data(., depvar, st_local("touse"))
    xT = st_data(., thrvar, st_local("touse"))
    if (othervars != "") W = st_data(., othervars, st_local("touse"))
    else W = J(rows(y), 0, .)
    tv = st_data(., st_local("timevar"), st_local("touse"))

    // contiguity check
    tdelta = strtoreal(st_local("tdelta"))
    if (rows(tv) > 1) {
        ndiff = sum( (tv[|2 \ rows(tv)|] - tv[|1 \ rows(tv)-1|]) :!= tdelta )
        if (ndiff > 0) {
            st_local("__err", "gaps in the time series are not supported; sample must be contiguous after tsset")
            return
        }
    }

    if (rows(y) < p+q+10) {
        st_local("__err", "too few observations for the requested lag structure")
        return
    }

    // ---- candidate thresholds & minimum regime size ----
    dx = xT - _tn_lag(xT,1)
    cand = _tn_cand(dx, trim, G)
    minobs = max((ceil(0.05*rows(y)), p+q+2))

    // ---- information criteria over S = 1..maxreg ----
    real scalar Smax
    Smax = maxreg
    ictab = J(Smax, 7, .)
    real colvector ssrS_vec
    real matrix tau_store
    ssrS_vec = J(Smax,1,.)
    tau_store = J(Smax, 2, .)

    for (s=1; s<=Smax; s++) {
        ssr = _tn_profile(y, xT, W, p, q, s, cand, minobs, consf, trf, besttau, profile)
        if (ssr==.) continue
        // effective N for this S (rows used)
        d = _tn_design(y, xT, W, p, q, (s==1 ? J(1,0,.) : besttau), consf, trf)
        nN = rows(d.dyv)
        sigma2 = ssr/nN
        // parameter counts
        real scalar npar, nthrp, npar0
        npar  = cols(d.X)            // regression coefficients
        nthrp = s - 1                // thresholds
        // IC = log(sigma2) + (c_T/T)*npars
        real scalar cAIC, cSIC, cHQ
        cAIC = 2
        cSIC = log(nN)
        cHQ  = 2*log(log(nN))
        // standard ICs count thresholds; Pitarakis (p*) exclude them
        ictab[s,1] = s
        ictab[s,2] = log(sigma2) + (cAIC/nN)*(npar+nthrp)   // AIC
        ictab[s,3] = log(sigma2) + (cSIC/nN)*(npar+nthrp)   // SIC
        ictab[s,4] = log(sigma2) + (cHQ /nN)*(npar+nthrp)   // HQIC
        ictab[s,5] = log(sigma2) + (cAIC/nN)*(npar)         // pAIC
        ictab[s,6] = log(sigma2) + (cSIC/nN)*(npar)         // pSIC
        ictab[s,7] = log(sigma2) + (cHQ /nN)*(npar)         // pHQIC
        ssrS_vec[s] = ssr
        if (s >= 2) {
            tau_store[s,1] = besttau[1]
            if (s == 3) tau_store[s,2] = besttau[2]
        }
    }

    // ---- choose S ----
    real scalar iccol
    if (ic=="aic") iccol = 2
    else if (ic=="sic") iccol = 3
    else if (ic=="hqic") iccol = 4
    else if (ic=="paic") iccol = 5
    else if (ic=="psic") iccol = 6
    else iccol = 7

    if (fixedS != 0) {
        S = fixedS
    }
    else {
        real scalar bestic, bs
        bestic = .
        bs = 1
        for (s=1; s<=Smax; s++) {
            if (ictab[s,iccol] != . && ictab[s,iccol] < bestic) {
                bestic = ictab[s,iccol]
                bs = s
            }
        }
        S = bs
    }

    // ---- final fit at chosen S ----
    nthr = S - 1
    if (S == 1) tauS = J(1,0,.)
    else if (S == 2) tauS = tau_store[2,1]
    else tauS = (tau_store[3,1], tau_store[3,2])

    // if fixedS was given we may not have profiled it -- ensure tau available
    // (Mata does not short-circuit &&, so guard the subscript with a nested if:
    //  for S==1 tauS is the 1x0 empty vector and tauS[1] would be invalid)
    if (S >= 2) {
        if (tauS[1]==.) {
            ssr = _tn_profile(y, xT, W, p, q, S, cand, minobs, consf, trf, besttau, profile)
            tauS = besttau
        }
    }

    d = _tn_design(y, xT, W, p, q, tauS, consf, trf)
    ssr = _tn_ols(d.dyv, d.X, robust, b, V)
    nN = rows(d.dyv)
    k  = cols(d.X)
    sigma2 = ssr/(nN-k)
    rmse = sqrt(sigma2)
    tss = quadcross(d.dyv :- mean(d.dyv), d.dyv :- mean(d.dyv))
    r2  = 1 - ssr/tss
    r2a = 1 - (1-r2)*(nN-1)/(nN-k)

    real scalar rho
    rho = b[d.irho]

    // ---- long-run coefficients beta = -theta/rho (delta method) ----
    real scalar nlr, idx
    nlr = S + d.m
    lrb_v = J(nlr,1,.)
    lrV   = J(nlr,nlr,0)
    real rowvector lridx
    lridx = (d.ilr_xT, d.ilr_w)         // indices of long-run level coeffs
    // build via delta method using submatrix over (rho, theta_j)
    real scalar t1, t2
    real colvector grad
    real matrix Vsub
    for (i=1; i<=nlr; i++) {
        idx = lridx[i]
        lrb_v[i] = -b[idx]/rho
    }
    // covariance via stacked Jacobian: beta_i = -b[idx_i]/rho
    // d beta_i/d b[idx_i] = -1/rho ; d beta_i/d rho = b[idx_i]/rho^2
    real matrix Jac
    real rowvector allidx
    allidx = (d.irho, lridx)            // rho first, then thetas
    Vsub = V[allidx, allidx]
    Jac = J(nlr, 1+nlr, 0)
    for (i=1; i<=nlr; i++) {
        Jac[i,1]   = b[lridx[i]]/rho^2     // wrt rho
        Jac[i,1+i] = -1/rho                // wrt theta_i
    }
    lrV = Jac * Vsub * Jac'

    // ---- asymmetry tests ----
    // long-run: H0 theta^(s) all equal across regimes (s=1..S) -> on level coeffs of xT
    real scalar lachi2, ladf, lap, sachi2, sadf, sap
    if (S >= 2) {
        real matrix R, Vr, RbR
        real colvector Rb
        // long-run: differences of consecutive theta (xT regimes)
        R = J(S-1, k, 0)
        for (i=1; i<=S-1; i++) {
            R[i, d.ilr_xT[i]]   =  1
            R[i, d.ilr_xT[i+1]] = -1
        }
        Rb  = R*b
        RbR = R*V*R'
        lachi2 = (Rb' * invsym(RbR) * Rb)
        ladf = S-1
        lap = chi2tail(ladf, lachi2)

        // short-run: sum_j pi^(s) equal across regimes
        R = J(S-1, k, 0)
        for (i=1; i<=S-1; i++) {
            for (j=1; j<=q; j++) {
                R[i, d.ipi[i,j]]   =  1
                R[i, d.ipi[i+1,j]] = -1
            }
        }
        Rb  = R*b
        RbR = R*V*R'
        sachi2 = (Rb' * invsym(RbR) * Rb)
        sadf = S-1
        sap = chi2tail(sadf, sachi2)
    }
    else {
        lachi2=.; ladf=0; lap=.
        sachi2=.; sadf=0; sap=.
    }

    // ---- PSS bounds test : F for rho=0 and all long-run levels =0 ----
    real scalar pssk, Fpss, tBDM
    real rowvector bidx
    real matrix Rp, Vp
    real colvector Rbp
    pssk = nlr
    bidx = (d.irho, lridx)
    Rp = J(cols(bidx), k, 0)
    for (i=1; i<=cols(bidx); i++) Rp[i, bidx[i]] = 1
    Rbp = Rp*b
    Vp  = Rp*V*Rp'
    Fpss = (Rbp' * invsym(Vp) * Rbp) / cols(bidx)
    tBDM = b[d.irho]/sqrt(V[d.irho,d.irho])

    // ---- QLR (optional) ----
    if (doqlr) {
        real scalar nNq
        if (dots) printf("{txt}bootstrapping QLR (B=%g) ", B)
        _tn_qlr(y, xT, W, p, q, cand, minobs, consf, trf, B, nNq, qlr_obs, qlr_p, dots)
        st_numscalar("__qlr", qlr_obs)
        st_numscalar("__qlrp", qlr_p)
    }

    // ---- structural pieces for multipliers ----
    // theta : S-vector (level coeffs on xT regimes)
    thetaM = J(S,1,.)
    for (s=1; s<=S; s++) thetaM[s] = b[d.ilr_xT[s]]
    // pi : S x q
    piM = J(S,q,.)
    for (s=1; s<=S; s++) for (j=1; j<=q; j++) piM[s,j] = b[d.ipi[s,j]]
    // phi : (p-1) vector
    if (p > 1) {
        phiM = J(p-1,1,.)
        for (j=1; j<=p-1; j++) phiM[j] = b[d.iphi[j]]
        st_matrix("__phi", phiM)
    }

    // ---- write esample flag ----
    keepflag = J(st_nobs(),1,0)
    keepflag[rowsused[d.keep]] = J(rows(d.keep),1,1)
    st_store(., st_local("used"), keepflag)

    // ---- write results back ----
    st_matrix("__b", b')
    st_matrix("__V", V)
    st_matrix("__lrb", lrb_v')
    st_matrix("__lrV", lrV)
    st_matrix("__ictab", ictab[selectindex(ictab[.,1]:!=.), .])
    st_matrix("__theta", thetaM)
    st_matrix("__pimat", piM)
    if (S >= 2) st_matrix("__thr", tauS')
    if (S == 2 && rows(profile) > 0) st_matrix("__profile", profile)

    // long-run coefficient names (for colnames in ado)
    string scalar lrnm
    lrnm = ""
    for (s=1; s<=S; s++) lrnm = lrnm + " lr_" + thrvar + "_r" + strofreal(s)
    if (othervars != "") {
        string rowvector ow
        ow = tokens(othervars)
        for (j=1; j<=cols(ow); j++) lrnm = lrnm + " lr_" + ow[j]
    }
    st_local("__lrnames", lrnm)

    st_numscalar("__S", S)
    st_numscalar("__N", nN)
    st_numscalar("__k", k)
    st_numscalar("__df_r", nN-k)
    st_numscalar("__rho", rho)
    st_numscalar("__ssr", ssr)
    st_numscalar("__rmse", rmse)
    st_numscalar("__r2", r2)
    st_numscalar("__r2_a", r2a)
    st_numscalar("__lachi2", lachi2)
    st_numscalar("__ladf", ladf)
    st_numscalar("__lap", lap)
    st_numscalar("__sachi2", sachi2)
    st_numscalar("__sadf", sadf)
    st_numscalar("__sap", sap)
    st_numscalar("__pssk", pssk)
    st_numscalar("__Fpss", Fpss)
    st_numscalar("__tBDM", tBDM)
    st_local("__err", "")
}
end