*! tsadvroot 1.0.0  03jul2026
*! Advanced time-series unit-root tests:
*!   qadf  - quantile ADF (Koenker & Xiao 2004)
*!   fqadf - Fourier quantile ADF (Li & Zheng 2018) with residual bootstrap
*!   npadf - two-break ADF (Narayan & Popp 2010)
*!   cisur - GLS unit-root tests with multiple structural breaks
*!           (Carrion-i-Silvestre, Kim & Perron 2009)
*! Exact Stata translation of the GAUSS routines by Saban Nazlioglu (tspdlib)
*! and Josep Lluis Carrion-i-Silvestre (based on Ng & Perron 2001 code).
*! Author: Merwan Roudane, merwanroudane920@gmail.com
*! https://github.com/merwanroudane

program define tsadvroot, rclass
    version 14.0
    gettoken sub rest : 0, parse(" ,")
    local sub = lower(`"`sub'"')
    if ("`sub'" == "qadf") {
        tsadvroot_qadf `rest'
    }
    else if inlist("`sub'", "fqadf", "qfadf", "fourierqadf") {
        tsadvroot_fqadf `rest'
    }
    else if inlist("`sub'", "npadf", "np", "adf2b") {
        tsadvroot_npadf `rest'
    }
    else if inlist("`sub'", "cisur", "cis", "gls") {
        tsadvroot_cisur `rest'
    }
    else if ("`sub'" == "") {
        di as err "subcommand required: {bf:qadf} | {bf:fqadf} | {bf:npadf} | {bf:cisur}"
        di as err "see {helpb tsadvroot} for details"
        exit 198
    }
    else {
        di as err `"unknown tsadvroot subcommand "`sub'""'
        di as err "valid subcommands: {bf:qadf} | {bf:fqadf} | {bf:npadf} | {bf:cisur}"
        exit 199
    }
    return add
end

*==============================================================================
* Common checks: tsset, no panel, contiguous sample
* (touse is created by the CALLER and passed by name; see help gotchas)
*==============================================================================
program define _tsav_check, rclass
    version 14.0
    args touse
    capture qui tsset
    if _rc {
        di as err "data must be {helpb tsset} with a time variable"
        exit 459
    }
    if "`r(panelvar)'" != "" {
        di as err "tsadvroot works on a single time series; data are xtset with panel variable {bf:`r(panelvar)'}"
        di as err "use {cmd:if} to select one unit, or {cmd:tsset} the series"
        exit 459
    }
    local tvar "`r(timevar)'"
    qui tsreport if `touse'
    if r(N_gaps) > 0 {
        di as err "the estimation sample contains `r(N_gaps)' gap(s) in `tvar'"
        di as err "unit-root tests require a contiguous series (no gaps, no interior if/in holes)"
        exit 498
    }
    qui count if `touse'
    return scalar T = r(N)
    return local tvar "`tvar'"
end

*==============================================================================
* SUBCOMMAND 1 : qadf  --  Koenker & Xiao (2004) quantile ADF
* Exact translation of qr_adf.src (QRADF proc, fourier=0 path)
*==============================================================================
program define tsadvroot_qadf, rclass
    version 14.0
    syntax varname(ts) [if] [in] [, Tau(numlist >0 <1 sort) Model(string) ///
        PMax(integer 8) IC(string) GRaph NAme(string) noPRint ]

    local yv `varlist'

    * ---- defaults & validation -------------------------------------------
    if "`tau'" == "" local tau "0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9"
    if "`model'" == "" local model "c"
    local model = lower("`model'")
    if inlist("`model'", "c", "1", "constant") local mnum 1
    else if inlist("`model'", "ct", "2", "trend") local mnum 2
    else {
        di as err "model() must be {bf:c} (constant) or {bf:ct} (constant and trend)"
        exit 198
    }
    if "`ic'" == "" local ic "tstat"
    local ic = lower("`ic'")
    if inlist("`ic'", "aic", "1") local icn 1
    else if inlist("`ic'", "sic", "bic", "2") local icn 2
    else if inlist("`ic'", "tstat", "t", "3") local icn 3
    else {
        di as err "ic() must be {bf:aic}, {bf:sic} or {bf:tstat}"
        exit 198
    }
    if `pmax' < 0 {
        di as err "pmax() must be non-negative"
        exit 198
    }

    * ---- sample -----------------------------------------------------------
    marksample touse
    _tsav_check `touse'
    local T = r(T)
    local tvar "`r(tvar)'"
    if `T' < `pmax' + 15 {
        di as err "insufficient observations (T=`T') for pmax(`pmax')"
        exit 2001
    }

    * protect the user's e() results (we run qreg internally)
    tempname ehold
    capture _estimates hold `ehold', restore nullok

    * ---- lag selection: GAUSS { ADFt, p, cv } = ADF(y, 1, pmax, ic) --------
    mata: _tsav_adfsel_st("`yv'", "`touse'", `pmax', `icn')
    local p = r(adflag)

    * ---- build regressors on the estimation sample -------------------------
    * GAUSS trims the first p+1 observations: sample = positions p+2 .. T
    tempvar pos esamp y1
    qui gen long `pos' = sum(`touse')
    qui gen byte `esamp' = `touse' & (`pos' > `p' + 1)
    qui gen double `y1' = L.`yv' if `touse'
    local xvars "`y1'"
    forvalues j = 1/`p' {
        tempvar dy`j'
        qui gen double `dy`j'' = L`j'.D.`yv' if `touse'
        local xvars "`xvars' `dy`j''"
    }
    if `mnum' == 2 {
        tempvar trnd
        qui gen double `trnd' = `pos' if `touse'
        local xvars "`xvars' `trnd'"
    }
    qui count if `esamp'
    local neff = r(N)

    * ---- per-quantile computation ------------------------------------------
    local ntau : word count `tau'
    tempname R b0 b1 b2
    matrix `R' = J(`ntau', 8, .)
    local row 0
    foreach t of local tau {
        local ++row
        * bandwidth h (GAUSS __get_qr_adf_h)
        mata: st_local("h", strofreal(_tsav_h(`t', `neff'), "%18.0g"))
        local t1 = `t' + `h'
        if `t1' >= 1 local t1 = .9999
        local t2 = `t' - `h'
        if `t2' <= 0 local t2 = .0001
        capture qui qreg `yv' `xvars' if `esamp', quantile(`t')
        if _rc {
            di as err "qreg failed at tau = `t' (rc = " _rc ")"
            exit _rc
        }
        matrix `b0' = e(b)
        capture qui qreg `yv' `xvars' if `esamp', quantile(`t1')
        if _rc {
            di as err "qreg failed at tau+h = `t1' (rc = " _rc ")"
            exit _rc
        }
        matrix `b1' = e(b)
        capture qui qreg `yv' `xvars' if `esamp', quantile(`t2')
        if _rc {
            di as err "qreg failed at tau-h = `t2' (rc = " _rc ")"
            exit _rc
        }
        matrix `b2' = e(b)
        mata: _tsav_qadf_one("`R'", `row', "`yv'", "`xvars'", "`esamp'", ///
            `t', `h', "`b0'", "`b1'", "`b2'", `p', 0, `mnum')
    }
    matrix colnames `R' = tau rho_tau rho_ols delta2 tn cv1 cv5 cv10

    * ---- display ------------------------------------------------------------
    if "`print'" == "" {
        local mlab "constant"
        if `mnum' == 2 local mlab "constant and trend"
        local iclab "AIC"
        if `icn' == 2 local iclab "SIC"
        if `icn' == 3 local iclab "t-stat (1.645)"
        di
        di as text "{hline 82}"
        di as text "Quantile ADF unit-root test" _col(50) "Koenker & Xiao (2004, JASA)"
        di as text "{hline 82}"
        di as text "Variable      : " as result "`yv'" ///
            as text _col(44) "Obs (effective) = " as result %8.0f `neff'
        di as text "Model         : " as result "`mlab'" ///
            as text _col(44) "Max lags        = " as result %8.0f `pmax'
        di as text "Lag selection : " as result "`iclab'" ///
            as text _col(44) "Lags selected   = " as result %8.0f `p'
        di as text "{hline 82}"
        di as text "    tau    rho(tau)    rho(OLS)    delta{c 94}2     t_n(tau)" ///
            "        1%       5%      10%"
        di as text "{hline 82}"
        forvalues i = 1/`ntau' {
            local st ""
            if `R'[`i',5] < `R'[`i',6] local st "***"
            else if `R'[`i',5] < `R'[`i',7] local st "**"
            else if `R'[`i',5] < `R'[`i',8] local st "*"
            di as text %7.2f `R'[`i',1] ///
                as result %12.4f `R'[`i',2] %12.4f `R'[`i',3] ///
                %11.4f `R'[`i',4] %13.3f `R'[`i',5] as text %-4s "`st'" ///
                as result %8.3f `R'[`i',6] %9.3f `R'[`i',7] %9.3f `R'[`i',8]
        }
        di as text "{hline 82}"
        di as text "H0: unit root at quantile tau, i.e. rho(tau) = 1."
        di as text "Critical values: Hansen (1995), interpolated on delta{c 94}2" ///
            " as in Koenker & Xiao (2004)."
        di as text "Rejection: * p<0.10, ** p<0.05, *** p<0.01 (t_n below the critical value)."
        di as text "{hline 82}"
    }

    * ---- graph (before return matrix; gotcha: return matrix moves) ----------
    if "`graph'" != "" {
        if `ntau' < 2 {
            di as text "(graph skipped: needs at least two quantiles in tau())"
        }
        else {
            local gname "`name'"
            if "`gname'" == "" local gname "tsavqadf"
            tempname RG
            matrix `RG' = `R'
            preserve
            qui clear
            qui svmat double `RG', name(_tq)
            local gopt graphregion(color(white)) plotregion(color(white)) ///
                ylabel(, angle(horizontal) grid glcolor(gs14) glwidth(vthin)) ///
                xlabel(, grid glcolor(gs14) glwidth(vthin))
            twoway (line _tq2 _tq1, lcolor(navy) lwidth(medthick)) ///
                   (scatter _tq2 _tq1, mcolor(navy) msymbol(O)), ///
                yline(1, lpattern(dash) lcolor(gs8)) ///
                ytitle("{&rho}({&tau})") xtitle("") ///
                title("Persistence  {&rho}({&tau})", size(medsmall) color(black)) ///
                legend(off) `gopt' name(`gname'_rho, replace) nodraw
            twoway (line _tq6 _tq1, lpattern(shortdash) lcolor(gs6)) ///
                   (line _tq7 _tq1, lpattern(dash) lcolor(maroon)) ///
                   (line _tq8 _tq1, lpattern(longdash_dot) lcolor(gs10)) ///
                   (line _tq5 _tq1, lcolor(navy) lwidth(medthick)) ///
                   (scatter _tq5 _tq1, mcolor(navy) msymbol(O)), ///
                ytitle("t{sub:n}({&tau})") xtitle("Quantile {&tau}") ///
                title("Quantile ADF statistic", size(medsmall) color(black)) ///
                legend(order(4 "t{sub:n}({&tau})" 1 "1% cv" 2 "5% cv" 3 "10% cv") ///
                    rows(1) size(small) region(lstyle(none))) ///
                `gopt' name(`gname'_tn, replace) nodraw
            graph combine `gname'_rho `gname'_tn, cols(1) iscale(0.9) ///
                graphregion(color(white)) ///
                title("Quantile ADF unit-root test: `yv'", size(medium) color(black)) ///
                note("Koenker & Xiao (2004); Hansen (1995) critical values.", ///
                    size(vsmall)) name(`gname', replace)
            restore
        }
    }

    * ---- returns -------------------------------------------------------------
    return scalar T = `T'
    return scalar N = `neff'
    return scalar lags = `p'
    return local model = cond(`mnum' == 1, "c", "ct")
    return local ic "`ic'"
    return local tau "`tau'"
    return local varname "`yv'"
    return local cmd "tsadvroot qadf"
    if `ntau' == 1 {
        return scalar tau = `R'[1,1]
        return scalar rho_tau = `R'[1,2]
        return scalar rho_ols = `R'[1,3]
        return scalar delta2 = `R'[1,4]
        return scalar tn = `R'[1,5]
        return scalar cv1 = `R'[1,6]
        return scalar cv5 = `R'[1,7]
        return scalar cv10 = `R'[1,8]
    }
    return matrix results = `R', copy
end

*==============================================================================
* SUBCOMMAND 2 : fqadf  --  Li & Zheng (2018) Fourier quantile ADF
* Exact translation of qr_fourier_adf.src
* (QR_Fourier_ADF and QR_Fourier_ADF_bootstrap)
*==============================================================================
program define tsadvroot_fqadf, rclass
    version 14.0
    syntax varname(ts) [if] [in] [, Tau(numlist >0 <1 sort) Model(string) ///
        Lags(integer 8) Freq(integer 3) NBoot(integer 1000) SEED(integer 0) ///
        noBOOTstrap GRaph NAme(string) noPRint ]

    local yv `varlist'
    local p = `lags'
    local k = `freq'

    if "`tau'" == "" local tau "0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9"
    if "`model'" == "" local model "c"
    local model = lower("`model'")
    if inlist("`model'", "c", "1", "constant") local mnum 1
    else if inlist("`model'", "ct", "2", "trend") local mnum 2
    else {
        di as err "model() must be {bf:c} or {bf:ct}"
        exit 198
    }
    if `p' < 0 {
        di as err "lags() must be non-negative"
        exit 198
    }
    if `k' < 1 {
        di as err "freq() must be a positive integer"
        exit 198
    }
    if `nboot' < 0 {
        di as err "nboot() must be non-negative"
        exit 198
    }
    if "`bootstrap'" != "" local nboot 0

    marksample touse
    _tsav_check `touse'
    local T = r(T)
    local tvar "`r(tvar)'"
    if `T' < `p' + 20 {
        di as err "insufficient observations (T=`T') for lags(`p')"
        exit 2001
    }

    tempname ehold
    capture _estimates hold `ehold', restore nullok

    * ---- observed statistics (exact QR_Fourier_ADF) --------------------------
    local ntau : word count `tau'
    tempname R
    matrix `R' = J(`ntau', 8, .)
    tsadvroot_fqadf_core `yv' if `touse', taulist(`tau') model(`mnum') ///
        p(`p') k(`k') rmat(`R')
    local neff = r(neff)

    * ---- bootstrap (exact QR_Fourier_ADF_bootstrap) ---------------------------
    tempname BS
    matrix `BS' = J(`ntau', 7, .)
    if `nboot' > 0 {
        if "`print'" == "" {
            di as text "(running " as result `nboot' ///
                as text " bootstrap replications, this may take a while...)"
        }
        * residuals of the null model on the full sample -> __tsav_mu (Mata)
        mata: _tsav_fq_prep("`yv'", "`touse'", `mnum', `p', `k')
        local tt = r(tt)
        tempname tnobs BR
        matrix `tnobs' = `R'[1..., 5]
        matrix `BR' = J(`ntau', 8, .)
        preserve
        qui drop _all
        qui set obs `tt'
        tempvar tb yb
        qui gen long `tb' = _n
        qui tsset `tb'
        qui gen double `yb' = .
        local corecmd "quietly tsadvroot_fqadf_core `yb', taulist(`tau') model(`mnum') p(`p') k(`k') rmat(`BR')"
        mata: _tsav_fq_boot("`yb'", `tt', `nboot', st_local("corecmd"), ///
            "`BR'", `ntau', "`BS'", "`tnobs'", `seed')
        restore
        capture mata: mata drop __tsav_mu
    }
    matrix colnames `R' = tau rho_tau rho_ols delta2 tn cv1 cv5 cv10
    matrix colnames `BS' = cvlt1 cvlt5 cvlt10 cvsrc1 cvsrc5 cvsrc10 pboot

    * ---- display ---------------------------------------------------------------
    if "`print'" == "" {
        local mlab "constant"
        if `mnum' == 2 local mlab "constant and trend"
        di
        di as text "{hline 82}"
        di as text "Fourier quantile ADF unit-root test" _col(48) ///
            "Li & Zheng (2018, Fin Res Letters)"
        di as text "{hline 82}"
        di as text "Variable      : " as result "`yv'" ///
            as text _col(44) "Obs (effective) = " as result %8.0f `neff'
        di as text "Model         : " as result "`mlab'" ///
            as text _col(44) "Lags (fixed)    = " as result %8.0f `p'
        di as text "Fourier freq. : " as result "`k'" ///
            as text _col(44) "Bootstrap reps  = " as result %8.0f `nboot'
        di as text "{hline 82}"
        if `nboot' > 0 {
            di as text "    tau    rho(tau)     F-QADF t_n    boot p" ///
                "        1%       5%      10%   (bootstrap cv)"
        }
        else {
            di as text "    tau    rho(tau)     F-QADF t_n"
        }
        di as text "{hline 82}"
        forvalues i = 1/`ntau' {
            if `nboot' > 0 {
                local st ""
                if `R'[`i',5] < `BS'[`i',1] local st "***"
                else if `R'[`i',5] < `BS'[`i',2] local st "**"
                else if `R'[`i',5] < `BS'[`i',3] local st "*"
                di as text %7.2f `R'[`i',1] ///
                    as result %12.4f `R'[`i',2] %14.3f `R'[`i',5] ///
                    as text %-4s "`st'" ///
                    as result %8.3f `BS'[`i',7] ///
                    %10.3f `BS'[`i',1] %9.3f `BS'[`i',2] %9.3f `BS'[`i',3]
            }
            else {
                di as text %7.2f `R'[`i',1] ///
                    as result %12.4f `R'[`i',2] %14.3f `R'[`i',5]
            }
        }
        di as text "{hline 82}"
        di as text "H0: unit root at quantile tau (with smooth structural change)."
        if `nboot' > 0 {
            di as text "Displayed cv: left-tail (1st/5th/10th percentile)" ///
                " of the bootstrap distribution;"
            di as text "boot p = share of bootstrap t_n below the observed t_n."
            di as text "The GAUSS-source order statistics (0.99/0.95/0.90)" ///
                " are stored in r(boot) cols 4-6."
            di as text "Rejection: * p<0.10, ** p<0.05, *** p<0.01."
        }
        else {
            di as text "No critical values displayed: use the bootstrap" ///
                " (option nboot(#)) for inference."
        }
        di as text "{hline 82}"
    }

    * ---- graph -------------------------------------------------------------------
    if "`graph'" != "" {
        if `ntau' < 2 {
            di as text "(graph skipped: needs at least two quantiles in tau())"
        }
        else {
            local gname "`name'"
            if "`gname'" == "" local gname "tsavfqadf"
            tempname RG
            matrix `RG' = `R'[1..., 1], `R'[1..., 5], `BS'[1..., 1..3]
            preserve
            qui clear
            qui svmat double `RG', name(_tf)
            local gopt graphregion(color(white)) plotregion(color(white)) ///
                ylabel(, angle(horizontal) grid glcolor(gs14) glwidth(vthin)) ///
                xlabel(, grid glcolor(gs14) glwidth(vthin))
            if `nboot' > 0 {
                twoway (line _tf3 _tf1, lpattern(shortdash) lcolor(gs6)) ///
                       (line _tf4 _tf1, lpattern(dash) lcolor(maroon)) ///
                       (line _tf5 _tf1, lpattern(longdash_dot) lcolor(gs10)) ///
                       (line _tf2 _tf1, lcolor(navy) lwidth(medthick)) ///
                       (scatter _tf2 _tf1, mcolor(navy) msymbol(O)), ///
                    ytitle("t{sub:n}({&tau})") xtitle("Quantile {&tau}") ///
                    title("Fourier quantile ADF test: `yv'", size(medium) color(black)) ///
                    subtitle("Fourier frequency k = `k', p = `p' lags", ///
                        size(small)) ///
                    legend(order(4 "F-QADF t{sub:n}({&tau})" 1 "1% cv" ///
                        2 "5% cv" 3 "10% cv") rows(1) size(small) ///
                        region(lstyle(none))) ///
                    note("Bootstrap critical values (left tail), B = `nboot'." ///
                        , size(vsmall)) ///
                    `gopt' name(`gname', replace)
            }
            else {
                twoway (line _tf2 _tf1, lcolor(navy) lwidth(medthick)) ///
                       (scatter _tf2 _tf1, mcolor(navy) msymbol(O)), ///
                    ytitle("t{sub:n}({&tau})") xtitle("Quantile {&tau}") ///
                    title("Fourier quantile ADF test: `yv'", size(medium) color(black)) ///
                    legend(off) `gopt' name(`gname', replace)
            }
            restore
        }
    }

    * ---- returns --------------------------------------------------------------------
    return scalar T = `T'
    return scalar N = `neff'
    return scalar lags = `p'
    return scalar k = `k'
    return scalar nboot = `nboot'
    return local model = cond(`mnum' == 1, "c", "ct")
    return local tau "`tau'"
    return local varname "`yv'"
    return local cmd "tsadvroot fqadf"
    if `ntau' == 1 {
        return scalar tn = `R'[1,5]
        return scalar rho_tau = `R'[1,2]
        return scalar delta2 = `R'[1,4]
        if `nboot' > 0 {
            return scalar pboot = `BS'[1,7]
            return scalar cv1 = `BS'[1,1]
            return scalar cv5 = `BS'[1,2]
            return scalar cv10 = `BS'[1,3]
            return scalar cvsrc1 = `BS'[1,4]
            return scalar cvsrc5 = `BS'[1,5]
            return scalar cvsrc10 = `BS'[1,6]
        }
    }
    if `nboot' > 0 {
        return matrix boot = `BS', copy
    }
    return matrix results = `R', copy
end

*------------------------------------------------------------------------------
* fqadf computational core -- shared by the observed test and the bootstrap.
* Assumes a tsset, contiguous series. Fills rmat (ntau x 8):
*   tau  rho_tau  .  delta2  tn  .  .  .
* Exact translation of QR_Fourier_ADF (print branch excluded):
*   trims p+1 obs, THEN builds the Fourier terms on the trimmed length,
*   x = [y1, dyl(1..p), sin, cos, (trend = original index)].
*------------------------------------------------------------------------------
program define tsadvroot_fqadf_core, rclass
    version 14.0
    syntax varname [if] [, TAUlist(numlist >0 <1) MODel(integer 1) ///
        P(integer 8) K(integer 3) RMAT(string) ]

    local yv `varlist'
    marksample touse
    qui count if `touse'
    local T = r(N)
    local neff = `T' - `p' - 1
    if `neff' < `p' + 8 {
        di as err "fqadf core: series too short (T=`T', lags=`p')"
        exit 2001
    }

    tempvar pos esamp y1 sink cosk
    qui gen long `pos' = sum(`touse')
    qui gen byte `esamp' = `touse' & (`pos' > `p' + 1)
    qui gen double `y1' = L.`yv' if `touse'
    local xvars "`y1'"
    forvalues j = 1/`p' {
        tempvar dy`j'
        qui gen double `dy`j'' = L`j'.D.`yv' if `touse'
        local xvars "`xvars' `dy`j''"
    }
    * Fourier terms: sequence restarts at 1 on the trimmed sample and the
    * denominator is the trimmed length (exact source behaviour)
    qui gen double `sink' = sin(2*_pi*`k'*(`pos'-`p'-1)/`neff') if `esamp'
    qui gen double `cosk' = cos(2*_pi*`k'*(`pos'-`p'-1)/`neff') if `esamp'
    local xvars "`xvars' `sink' `cosk'"
    if `model' == 2 {
        tempvar trnd
        qui gen double `trnd' = `pos' if `touse'
        local xvars "`xvars' `trnd'"
    }

    tempname b0 b1 b2
    local row 0
    foreach t of local taulist {
        local ++row
        mata: st_local("h", strofreal(_tsav_h(`t', `neff'), "%18.0g"))
        local t1 = `t' + `h'
        if `t1' >= 1 local t1 = .9999
        local t2 = `t' - `h'
        if `t2' <= 0 local t2 = .0001
        capture qui qreg `yv' `xvars' if `esamp', quantile(`t')
        if _rc {
            di as err "qreg failed at tau = `t' (rc = " _rc ")"
            exit _rc
        }
        matrix `b0' = e(b)
        capture qui qreg `yv' `xvars' if `esamp', quantile(`t1')
        if _rc {
            di as err "qreg failed at tau+h = `t1' (rc = " _rc ")"
            exit _rc
        }
        matrix `b1' = e(b)
        capture qui qreg `yv' `xvars' if `esamp', quantile(`t2')
        if _rc {
            di as err "qreg failed at tau-h = `t2' (rc = " _rc ")"
            exit _rc
        }
        matrix `b2' = e(b)
        mata: _tsav_qadf_one("`rmat'", `row', "`yv'", "`xvars'", "`esamp'", ///
            `t', `h', "`b0'", "`b1'", "`b2'", `p', 1, `model')
    }
    return scalar neff = `neff'
end

*==============================================================================
* SUBCOMMAND 3 : npadf  --  Narayan & Popp (2010) two-break unit-root test
* Exact translation of narayan pop.src (ADF_2breaks)
*==============================================================================
program define tsadvroot_npadf, rclass
    version 14.0
    syntax varname(ts) [if] [in] [, Model(string) PMax(integer 8) IC(string) ///
        TRim(real 0.10) GRaph NAme(string) noPRint ]

    local yv `varlist'
    if "`model'" == "" local model "1"
    local model = lower("`model'")
    if inlist("`model'", "1", "a", "level", "m1") local mnum 1
    else if inlist("`model'", "2", "c", "both", "m2") local mnum 2
    else {
        di as err "model() must be {bf:1} (M1: breaks in level) or {bf:2}" ///
            " (M2: breaks in level and slope)"
        exit 198
    }
    if "`ic'" == "" local ic "tstat"
    local ic = lower("`ic'")
    if inlist("`ic'", "aic", "1") local icn 1
    else if inlist("`ic'", "sic", "bic", "2") local icn 2
    else if inlist("`ic'", "tstat", "t", "3") local icn 3
    else {
        di as err "ic() must be {bf:aic}, {bf:sic} or {bf:tstat}"
        exit 198
    }
    if `trim' <= 0 | `trim' >= 0.5 {
        di as err "trim() must be strictly between 0 and 0.5"
        exit 198
    }

    marksample touse
    _tsav_check `touse'
    local T = r(T)
    local tvar "`r(tvar)'"
    if `T' < 2*`pmax' + 20 {
        di as err "insufficient observations (T=`T') for pmax(`pmax')"
        exit 2001
    }

    if "`print'" == "" {
        di as text "(grid search over two break dates, this may take a moment...)"
    }
    mata: _tsav_np("`yv'", "`touse'", `mnum', `pmax', `icn', `trim')
    local stat = r(npstat)
    local tb1p = r(tb1)
    local tb2p = r(tb2)
    local plag = r(nplag)
    local cv1 = r(cv1)
    local cv5 = r(cv5)
    local cv10 = r(cv10)

    * map break positions to time values
    tempvar pos
    qui gen long `pos' = sum(`touse')
    qui su `tvar' if `pos' == `tb1p' & `touse', meanonly
    local tb1 = r(min)
    qui su `tvar' if `pos' == `tb2p' & `touse', meanonly
    local tb2 = r(min)
    local tf : format `tvar'
    local tb1s = string(`tb1', "`tf'")
    local tb2s = string(`tb2', "`tf'")

    local st ""
    if `stat' < `cv1' local st "***"
    else if `stat' < `cv5' local st "**"
    else if `stat' < `cv10' local st "*"

    if "`print'" == "" {
        local mlab "M1: two breaks in level"
        if `mnum' == 2 local mlab "M2: two breaks in level and slope"
        local iclab "AIC"
        if `icn' == 2 local iclab "SIC"
        if `icn' == 3 local iclab "t-stat (1.645)"
        di
        di as text "{hline 78}"
        di as text "Narayan & Popp (2010) unit-root test with two structural breaks"
        di as text "{hline 78}"
        di as text "Variable      : " as result "`yv'" ///
            as text _col(44) "Obs             = " as result %8.0f `T'
        di as text "Model         : " as result "`mlab'"
        di as text "Lag selection : " as result "`iclab'" ///
            as text _col(44) "Lags selected   = " as result %8.0f `plag'
        di as text "Trimming      : " as result %5.2f `trim'
        di as text "{hline 78}"
        di as text "   ADF-stat        TB1          TB2          1%       5%      10%"
        di as text "{hline 78}"
        di as result %10.3f `stat' as text "`st'" ///
            as result _col(17) %10s "`tb1s'" _col(30) %10s "`tb2s'" ///
            _col(43) %9.3f `cv1' %9.3f `cv5' %9.3f `cv10'
        di as text "{hline 78}"
        di as text "H0: unit root with two breaks in the DGP."
        di as text "Break dates are the last period of each pre-break regime" ///
            " (DU_t = 1 if t > TB)."
        di as text "Critical values: Narayan & Popp (2010), Table 3, T = `T'."
        di as text "Rejection: * p<0.10, ** p<0.05, *** p<0.01."
        di as text "{hline 78}"
    }

    if "`graph'" != "" {
        local gname "`name'"
        if "`gname'" == "" local gname "tsavnpadf"
        local subt = "ADF = " + string(`stat', "%9.3f") + "`st'" + ///
            ",  TB1 = `tb1s',  TB2 = `tb2s'"
        twoway (line `yv' `tvar' if `touse', lcolor(navy) lwidth(medthin)), ///
            xline(`tb1' `tb2', lpattern(dash) lcolor(maroon)) ///
            ytitle("`yv'") xtitle("") ///
            title("Narayan-Popp two-break test: `yv'", size(medium) color(black)) ///
            subtitle("`subt'", size(small)) ///
            note("Dashed lines: estimated break dates.", size(vsmall)) ///
            graphregion(color(white)) plotregion(color(white)) ///
            ylabel(, angle(horizontal) grid glcolor(gs14) glwidth(vthin)) ///
            name(`gname', replace)
    }

    return scalar T = `T'
    return scalar stat = `stat'
    return scalar lags = `plag'
    return scalar tb1 = `tb1'
    return scalar tb2 = `tb2'
    return scalar tb1pos = `tb1p'
    return scalar tb2pos = `tb2p'
    return scalar frac1 = `tb1p'/`T'
    return scalar frac2 = `tb2p'/`T'
    return scalar cv1 = `cv1'
    return scalar cv5 = `cv5'
    return scalar cv10 = `cv10'
    return local model "`mnum'"
    return local ic "`ic'"
    return local breakdates "`tb1s' `tb2s'"
    return local varname "`yv'"
    return local cmd "tsadvroot npadf"
end

*==============================================================================
* SUBCOMMAND 4 : cisur  --  Carrion-i-Silvestre, Kim & Perron (2009)
* GLS-based unit-root tests with multiple structural breaks.
* Exact translation of carrion silvestre2009.src (brute-force path,
* the source default: sburControlCreate() sets estimation = 0)
*==============================================================================
program define tsadvroot_cisur, rclass
    version 14.0
    syntax varname(ts) [if] [in] [, Model(string) Breaks(integer 1) ///
        BREAKDates(numlist) Penalty(string) KMax(integer 4) KMin(integer 0) ///
        GRaph NAme(string) noPRint ]

    local yv `varlist'
    if "`model'" == "" local model "break"
    local model = lower("`model'")
    if inlist("`model'", "0", "const", "constant") local mnum 0
    else if inlist("`model'", "1", "trend") local mnum 1
    else if inlist("`model'", "2", "slope") local mnum 2
    else if inlist("`model'", "3", "break", "both") local mnum 3
    else {
        di as err "model() must be one of: {bf:const} (0, constant, no breaks)," ///
            " {bf:trend} (1, linear trend, no breaks),"
        di as err "{bf:slope} (2, breaks in trend slope) or {bf:break}" ///
            " (3, breaks in level and slope)"
        exit 198
    }
    if "`penalty'" == "" local penalty "maic"
    local penalty = lower("`penalty'")
    if inlist("`penalty'", "maic", "0") local pen 0
    else if inlist("`penalty'", "bic", "1") local pen 1
    else {
        di as err "penalty() must be {bf:maic} or {bf:bic}"
        exit 198
    }
    if `kmin' < 0 | `kmax' < `kmin' {
        di as err "require 0 <= kmin() <= kmax()"
        exit 198
    }

    marksample touse
    _tsav_check `touse'
    local T = r(T)
    local tvar "`r(tvar)'"
    if `T' < `kmax' + 25 {
        di as err "insufficient observations (T=`T')"
        exit 2001
    }

    * ---- break setup -----------------------------------------------------------
    local known 0
    local m = `breaks'
    tempvar pos
    qui gen long `pos' = sum(`touse')
    tempname TB
    if `mnum' >= 2 {
        if "`breakdates'" != "" {
            local known 1
            local m : word count `breakdates'
            if `m' > 5 {
                di as err "at most 5 known break dates are allowed"
                exit 198
            }
            matrix `TB' = J(`m', 1, .)
            local i 0
            local prev = -1e300
            foreach d of numlist `breakdates' {
                local ++i
                if `d' <= `prev' {
                    di as err "breakdates() must be strictly increasing"
                    exit 198
                }
                local prev = `d'
                qui su `pos' if `touse' & float(`tvar') == float(`d'), meanonly
                if r(N) == 0 {
                    di as err "breakdates(): value `d' not found in `tvar'" ///
                        " within the sample"
                    exit 198
                }
                matrix `TB'[`i', 1] = r(min)
            }
        }
        else {
            if `m' < 1 | `m' > 3 {
                di as err "with unknown break dates the brute-force search" ///
                    " (the GAUSS-source default) supports breaks(1) to breaks(3);"
                di as err "for 4 or 5 breaks supply known dates via breakdates()"
                exit 198
            }
            if `m' == 3 & `T' > 150 & "`print'" == "" {
                di as text "(3 unknown breaks with T=`T': the O(T{c 94}3)" ///
                    " grid search may take several minutes)"
            }
        }
    }
    else {
        local m 0
    }
    if !`known' matrix `TB' = J(1, 1, 0)

    if "`print'" == "" & `mnum' >= 2 & !`known' {
        di as text "(searching for `m' break(s) by GLS-SSR minimization...)"
    }

    tempvar fitv
    qui gen double `fitv' = .
    tempname ST CV TBOUT
    mata: _tsav_cis("`yv'", "`touse'", `mnum', `known', "`TB'", `m', ///
        `pen', `kmax', `kmin', "`fitv'", "`ST'", "`CV'", "`TBOUT'")
    local cbar = r(cbar)
    local krule = r(krule)

    * stats vector: pt mpt adf za mza msb mzt
    local pt = `ST'[1,1]
    local mpt = `ST'[1,2]
    local adf = `ST'[1,3]
    local za = `ST'[1,4]
    local mza = `ST'[1,5]
    local msb = `ST'[1,6]
    local mzt = `ST'[1,7]

    * break dates in time units
    local bdates ""
    local bxline ""
    local nbfound = rowsof(`TBOUT')
    if `mnum' >= 2 {
        local tf : format `tvar'
        forvalues i = 1/`nbfound' {
            local bp = `TBOUT'[`i',1]
            qui su `tvar' if `pos' == `bp' & `touse', meanonly
            local bd = r(min)
            local bds = string(`bd', "`tf'")
            local bdates "`bdates' `bds'"
            local bxline "`bxline' `bd'"
        }
    }

    * ---- display -------------------------------------------------------------------
    if "`print'" == "" {
        local mlab "constant, no breaks (Model 0)"
        if `mnum' == 1 local mlab "linear trend, no breaks (Model 1)"
        if `mnum' == 2 local mlab "breaks in the trend slope (Model 2)"
        if `mnum' == 3 local mlab "breaks in level and slope (Model 3)"
        local blab "estimated (GLS-SSR minimization)"
        if `known' local blab "known (user supplied)"
        local plab "MAIC"
        if `pen' == 1 local plab "BIC"
        di
        di as text "{hline 78}"
        di as text "GLS unit-root tests with multiple structural breaks"
        di as text "Carrion-i-Silvestre, Kim & Perron (2009, Econometric Theory)"
        di as text "{hline 78}"
        di as text "Variable   : " as result "`yv'" ///
            as text _col(44) "Obs           = " as result %8.0f `T'
        di as text "Model      : " as result "`mlab'"
        if `mnum' >= 2 {
            di as text "Breaks     : " as result "`nbfound'" ///
                as text "  `blab'"
            di as text "Break dates: " as result "`bdates'"
        }
        di as text "c-bar      : " as result %8.4f `cbar' ///
            as text _col(44) "Lags (`plab')   = " as result %8.0f `krule'
        di as text "{hline 78}"
        di as text "   Test        Statistic        1%        5%       10%    Decision"
        di as text "{hline 78}"
        * cv rows in CV: 1=msb 2=mza 3=mzt 4=pt ; mapping as in the source
        local names `""PT" "MPT" "ADF" "ZA" "MZA" "MSB" "MZT""'
        local cvrow "4 4 3 2 2 1 3"
        forvalues i = 1/7 {
            local nm : word `i' of `names'
            local cr : word `i' of `cvrow'
            local s = `ST'[1, `i']
            local c1 = `CV'[`cr', 1]
            local c5 = `CV'[`cr', 2]
            local c10 = `CV'[`cr', 3]
            local dec "Cannot reject H0"
            local st ""
            if `s' < `c1' {
                local dec "Reject H0 (1%)"
                local st "***"
            }
            else if `s' < `c5' {
                local dec "Reject H0 (5%)"
                local st "**"
            }
            else if `s' < `c10' {
                local dec "Reject H0 (10%)"
                local st "*"
            }
            di as text %7s "`nm'" as result %14.3f `s' as text %-4s "`st'" ///
                as result %9.3f `c1' %10.3f `c5' %10.3f `c10' ///
                as text "   `dec'"
        }
        di as text "{hline 78}"
        di as text "H0: unit root. All tests reject for values BELOW the critical value."
        di as text "Critical values from the response surfaces in the GAUSS code of"
        di as text "Carrion-i-Silvestre et al. (2009), evaluated at the break fractions."
        di as text "* p<0.10, ** p<0.05, *** p<0.01."
        di as text "{hline 78}"
    }

    * ---- graph -----------------------------------------------------------------------
    if "`graph'" != "" {
        local gname "`name'"
        if "`gname'" == "" local gname "tsavcisur"
        local xl ""
        if "`bxline'" != "" local xl xline(`bxline', lpattern(dash) lcolor(maroon))
        local ngtxt = "Dashed lines: break dates. MZt = " + ///
            string(`mzt', "%9.3f") + ", MSB = " + string(`msb', "%9.3f") + "."
        twoway (line `yv' `tvar' if `touse', lcolor(navy) lwidth(medthin)) ///
               (line `fitv' `tvar' if `touse', lcolor(dkorange) ///
                lwidth(medthick) lpattern(solid)), ///
            `xl' ytitle("`yv'") xtitle("") ///
            title("CiS-Kim-Perron GLS test: `yv'", size(medium) color(black)) ///
            subtitle("Broken deterministic trend (GLS estimates)", size(small)) ///
            legend(order(1 "`yv'" 2 "GLS trend") rows(1) size(small) ///
                region(lstyle(none))) ///
            note("`ngtxt'", size(vsmall)) ///
            graphregion(color(white)) plotregion(color(white)) ///
            ylabel(, angle(horizontal) grid glcolor(gs14) glwidth(vthin)) ///
            name(`gname', replace)
    }

    * ---- returns -----------------------------------------------------------------------
    return scalar T = `T'
    return scalar pt = `pt'
    return scalar mpt = `mpt'
    return scalar adf = `adf'
    return scalar za = `za'
    return scalar mza = `mza'
    return scalar msb = `msb'
    return scalar mzt = `mzt'
    return scalar cbar = `cbar'
    return scalar lags = `krule'
    return scalar nbreaks = cond(`mnum' >= 2, `nbfound', 0)
    return local model "`mnum'"
    return local penalty "`penalty'"
    return local breakdates "`bdates'"
    return local varname "`yv'"
    return local cmd "tsadvroot cisur"
    matrix rownames `CV' = MSB MZA MZT PT
    matrix colnames `CV' = cv1 cv5 cv10
    return matrix cv = `CV', copy
    if `mnum' >= 2 {
        return matrix breakpos = `TBOUT', copy
    }
    tempname STAB
    matrix `STAB' = `ST'
    matrix colnames `STAB' = PT MPT ADF ZA MZA MSB MZT
    return matrix stats = `STAB', copy
end

*==============================================================================
* Mata computational engine
*==============================================================================
version 14.0
mata:

// ---------------------------------------------------------------------------
// tspdlib _get_lag: index (1..pmax+1) of the selected lag+1
// ic 1 = AIC (first minimum), 2 = SIC, 3 = general-to-specific t-stat (1.645)
// ---------------------------------------------------------------------------
real scalar _tsav_getlagidx(real scalar ic, real colvector aicp,
    real colvector sicp, real colvector tstatp)
{
    real scalar pidx, j, n
    real colvector v

    n = rows(aicp)
    pidx = 1
    if (ic == 1) {
        v = aicp
        for (j = 2; j <= n; j++) {
            if (v[j] < v[pidx]) {
                pidx = j
            }
        }
    }
    else if (ic == 2) {
        v = sicp
        for (j = 2; j <= n; j++) {
            if (v[j] < v[pidx]) {
                pidx = j
            }
        }
    }
    else {
        pidx = 1
        j = n
        while (j >= 2) {
            if (tstatp[j] > 1.645) {
                pidx = j
                break
            }
            j = j - 1
        }
    }
    return(pidx)
}

// ---------------------------------------------------------------------------
// tspdlib ADF(y, 1, pmax, ic) lag selection (constant model)
// dep sample: dy trimmed twice (once for diff, once for p+1), regressors
// x = [y_{t-1}, const, dy lags]; aic/sic use the (k+2) penalty of the source
// Returns the selected LAG COUNT via r(adflag)
// ---------------------------------------------------------------------------
void _tsav_adfsel_st(string scalar yname, string scalar touse,
    real scalar pmax, real scalar ic)
{
    real colvector y
    real scalar p

    y = st_data(., yname, touse)
    p = _tsav_adfsel(y, pmax, ic)
    st_numscalar("r(adflag)", p)
}

real scalar _tsav_adfsel(real colvector y, real scalar pmax, real scalar ic)
{
    real scalar T, p, lo, n1, kx, pidx, j
    real colvector dy, y1, dep, bb, ee, sevec, taup, aicp, sicp, tstatp
    real matrix lmat, X, XX

    T = rows(y)
    dy = y[2::T] - y[1::(T-1)]
    y1 = y[1::(T-1)]
    lmat = J(T-1, pmax, 0)
    for (j = 1; j <= pmax; j++) {
        if (T-1-j >= 1) {
            lmat[(j+1)::(T-1), j] = dy[1::(T-1-j)]
        }
    }
    taup = J(pmax+1, 1, .)
    aicp = J(pmax+1, 1, .)
    sicp = J(pmax+1, 1, .)
    tstatp = J(pmax+1, 1, .)
    for (p = 0; p <= pmax; p++) {
        lo = p + 2
        dep = dy[lo::(T-1)]
        X = y1[lo::(T-1)], J(rows(dep), 1, 1)
        if (p > 0) {
            X = X, lmat[lo::(T-1), 1::p]
        }
        n1 = rows(dep)
        kx = cols(X)
        XX = invsym(cross(X, X))
        bb = XX*cross(X, dep)
        ee = dep - X*bb
        sevec = sqrt(diagonal(XX)*(cross(ee, ee)/(n1-kx)))
        taup[p+1] = bb[1]/sevec[1]
        aicp[p+1] = ln(cross(ee, ee)/n1) + 2*(kx+2)/n1
        sicp[p+1] = ln(cross(ee, ee)/n1) + (kx+2)*ln(n1)/n1
        tstatp[p+1] = abs(bb[kx]/sevec[kx])
    }
    pidx = _tsav_getlagidx(ic, aicp, sicp, tstatp)
    return(pidx - 1)
}

// ---------------------------------------------------------------------------
// bandwidth (GAUSS: bandwidth + __get_qr_adf_h) -- Hall-Sheather then Bofinger
// ---------------------------------------------------------------------------
real scalar _tsav_bwhs(real scalar tau, real scalar n)
{
    real scalar x0, f0

    x0 = invnormal(tau)
    f0 = normalden(x0)
    return(n^(-1/3) * invnormal(1 - 0.05/2)^(2/3) * ((1.5*f0^2)/(2*x0^2 + 1))^(1/3))
}

real scalar _tsav_bwbof(real scalar tau, real scalar n)
{
    real scalar x0, f0

    x0 = invnormal(tau)
    f0 = normalden(x0)
    return(n^(-0.2) * ((4.5*f0^4)/(2*x0^2 + 1)^2)^0.2)
}

real scalar _tsav_h(real scalar tau, real scalar n)
{
    real scalar h

    h = _tsav_bwhs(tau, n)
    if (tau <= 0.5 & h > tau) {
        h = _tsav_bwbof(tau, n)
        if (h > tau) {
            h = tau/1.5
        }
    }
    if (tau > 0.5 & h > 1 - tau) {
        h = _tsav_bwbof(tau, n)
        if (h > 1 - tau) {
            h = (1 - tau)/1.5
        }
    }
    return(h)
}

// ---------------------------------------------------------------------------
// Hansen (1995) critical values, interpolated on delta2 (GAUSS crit_QRadf)
// model: 0 = no deterministic, 1 = constant, 2 = constant + trend
// ---------------------------------------------------------------------------
real rowvector _tsav_qadfcv(real scalar r2, real scalar model)
{
    real matrix crt
    real rowvector ct
    real scalar r210, r2a, r2b, wa

    if (model == 0) {
        crt = (-2.4611512, -1.7832090, -1.4189957 \
               -2.4943410, -1.8184897, -1.4589747 \
               -2.5152783, -1.8516957, -1.5071775 \
               -2.5509773, -1.8957720, -1.5323511 \
               -2.5520784, -1.8949965, -1.5418830 \
               -2.5490848, -1.8981677, -1.5625462 \
               -2.5547456, -1.9343180, -1.5889045 \
               -2.5761273, -1.9387996, -1.6020210 \
               -2.5511921, -1.9328373, -1.6128210 \
               -2.5658, -1.9393, -1.6156)
    }
    else if (model == 1) {
        crt = (-2.7844267, -2.1158290, -1.7525193 \
               -2.9138762, -2.2790427, -1.9172046 \
               -3.0628184, -2.3994711, -2.0573070 \
               -3.1376157, -2.5070473, -2.1680520 \
               -3.1914660, -2.5841611, -2.2520173 \
               -3.2437157, -2.6399560, -2.3163270 \
               -3.2951006, -2.7180169, -2.4085640 \
               -3.3627161, -2.7536756, -2.4577709 \
               -3.3896556, -2.8074982, -2.5037759 \
               -3.4336, -2.8621, -2.5671)
    }
    else {
        crt = (-2.9657928, -2.3081543, -1.9519926 \
               -3.1929596, -2.5482619, -2.1991651 \
               -3.3727717, -2.7283918, -2.3806008 \
               -3.4904849, -2.8669056, -2.5315918 \
               -3.6003166, -2.9853079, -2.6672416 \
               -3.6819803, -3.0954760, -2.7815263 \
               -3.7551759, -3.1783550, -2.8728146 \
               -3.8348596, -3.2674954, -2.9735550 \
               -3.8800989, -3.3316415, -3.0364171 \
               -3.9638, -3.4126, -3.1279)
    }
    if (r2 < 0.1) {
        ct = crt[1, .]
    }
    else {
        r210 = r2*10
        if (r210 >= 10) {
            ct = crt[10, .]
        }
        else {
            r2a = floor(r210)
            r2b = ceil(r210)
            if (r2a < 1) {
                r2a = 1
            }
            if (r2a == r2b) {
                ct = crt[r2a, .]
            }
            else {
                wa = r2b - r210
                ct = wa*crt[r2a, .] + (1 - wa)*crt[r2b, .]
            }
        }
    }
    return(ct)
}

// ---------------------------------------------------------------------------
// per-quantile QADF computation (GAUSS QRADF tail / __get_qr_adf_stat /
// __get_qr_adf_delta2), given the three qreg coefficient vectors.
// usef = 0 : plain QADF  (w = dy, xx = [1, dyl], Hansen cvs)
// usef = 1 : Fourier QADF (w = residuals, xx = [1, all x but y1], no cvs)
// Fills row `row' of matrix `rmatname':
//   tau rho_tau rho_ols delta2 tn cv1 cv5 cv10
// ---------------------------------------------------------------------------
void _tsav_qadf_one(string scalar rmatname, real scalar row,
    string scalar yname, string scalar xvars, string scalar esamp,
    real scalar tau, real scalar h,
    string scalar b0name, string scalar b1name, string scalar b2name,
    real scalar p, real scalar usef, real scalar model)
{
    real colvector Y, y1, bg0, bg1, bg2, res, ind, phi, w, bols, tvec
    real matrix X, Xc, xx, ixx, Rm
    real rowvector b0, b1, b2, z1m, cv
    real scalar n, k, rho_tau, rho_ols, q1, q2, dq, fz, y1p, stat
    real scalar mw, mphi, covv, sdw, delta2

    Y = st_data(., yname, esamp)
    X = st_data(., tokens(xvars), esamp)
    n = rows(Y)
    Xc = J(n, 1, 1), X

    // reorder e(b): Stata puts _cons last, GAUSS puts it first
    b0 = st_matrix(b0name)
    k = cols(b0)
    bg0 = b0[k] \ b0[1::(k-1)]'
    b1 = st_matrix(b1name)
    bg1 = b1[k] \ b1[1::(k-1)]'
    b2 = st_matrix(b2name)
    bg2 = b2[k] \ b2[1::(k-1)]'

    rho_tau = bg0[2]

    // OLS rho (GAUSS: beta_ols = y/(ones~x))
    bols = invsym(cross(Xc, Xc))*cross(Xc, Y)
    rho_ols = bols[2]

    // density at the quantile: fz = 2h / (q1 - q2)
    z1m = 1, mean(X)
    q1 = z1m*bg1
    q2 = z1m*bg2
    dq = q1 - q2
    if (dq == 0) {
        fz = 0.01
    }
    else {
        fz = 2*h/dq
    }
    if (fz < 0) {
        fz = 0.01
    }

    // projection: xx = [1, dyl] (plain) or [1, all x except y1] (fourier)
    y1 = X[., 1]
    if (usef == 1) {
        if (cols(X) >= 2) {
            xx = J(n, 1, 1), X[., 2::cols(X)]
        }
        else {
            xx = J(n, 1, 1)
        }
    }
    else {
        if (p > 0) {
            xx = J(n, 1, 1), X[., 2::(p+1)]
        }
        else {
            xx = J(n, 1, 1)
        }
    }
    ixx = invsym(cross(xx, xx))
    tvec = cross(xx, y1)
    y1p = cross(y1, y1) - tvec'*ixx*tvec
    if (y1p < 0) {
        y1p = 0
    }
    stat = fz/sqrt(tau*(1 - tau)) * sqrt(y1p) * (rho_tau - 1)

    // delta2 (GAUSS __get_qr_adf_delta2)
    res = Y - Xc*bg0
    ind = res :< 0
    phi = J(n, 1, tau) - ind
    if (usef == 1) {
        w = res
    }
    else {
        w = Y - y1
    }
    mw = mean(w)
    mphi = mean(phi)
    covv = sum((w :- mw) :* (phi :- mphi))/(n - 1)
    sdw = sqrt(sum((w :- mw) :* (w :- mw))/(n - 1))
    delta2 = (covv/(sdw*sqrt(tau*(1 - tau))))^2

    if (usef == 0) {
        cv = _tsav_qadfcv(delta2, model)
    }
    else {
        cv = (., ., .)
    }

    Rm = st_matrix(rmatname)
    Rm[row, .] = (tau, rho_tau, rho_ols, delta2, stat, cv[1], cv[2], cv[3])
    st_matrix(rmatname, Rm)
}

// ---------------------------------------------------------------------------
// Fourier QADF bootstrap preparation (GAUSS QR_Fourier_ADF_bootstrap, part 1)
// Estimates the null model on the FULL sample (Fourier terms built on full T),
// fits AR(1) to its residuals and stores the centred innovations in the
// external __tsav_mu. Returns rows(mu) via r(tt).
// ---------------------------------------------------------------------------
void _tsav_fq_prep(string scalar yname, string scalar touse,
    real scalar model, real scalar p, real scalar k)
{
    external real colvector __tsav_mu
    real colvector y, s, sink, cosk, dyfull, yt, b, yd, yd1, yd0, mu
    real matrix X, Xt, L
    real scalar T, j, lo, n0, fi

    y = st_data(., yname, touse)
    T = rows(y)
    s = (1::T)
    sink = sin(2*pi()*k*s/T)
    cosk = cos(2*pi()*k*s/T)
    dyfull = J(T, 1, .)
    dyfull[2::T] = y[2::T] - y[1::(T-1)]
    X = J(T, 1, 1)
    if (model == 1) {
        X = X, sink, cosk
    }
    else {
        X = X, s, sink, cosk
    }
    if (p > 0) {
        L = J(T, p, .)
        for (j = 1; j <= p; j++) {
            if (T-j >= 2) {
                L[(j+2)::T, j] = dyfull[2::(T-j)]
            }
        }
        X = X, L
    }
    lo = p + 2
    Xt = X[lo::T, .]
    yt = y[lo::T]
    b = invsym(cross(Xt, Xt))*cross(Xt, yt)
    yd = yt - Xt*b
    n0 = rows(yd)
    yd1 = yd[1::(n0-1)]
    yd0 = yd[2::n0]
    fi = cross(yd1, yd0)/cross(yd1, yd1)
    mu = yd0 - yd1*fi
    mu = mu :- mean(mu)
    __tsav_mu = mu
    st_numscalar("r(tt)", rows(mu))
}

// ---------------------------------------------------------------------------
// Fourier QADF bootstrap loop (GAUSS QR_Fourier_ADF_bootstrap, part 2)
// For each replication: iid resample of mu with replacement, cumulate into a
// pure random walk, run the full Fourier-QADF core on the pseudo-series.
// Writes BSname (ntau x 7):
//   cvlt1 cvlt5 cvlt10  cvsrc1 cvsrc5 cvsrc10  pboot
// cvsrc* reproduce the GAUSS source order statistics (0.99/0.95/0.90 of the
// ascending-sorted bootstrap stats); cvlt* are the left-tail 1/5/10 percent
// order statistics used for the displayed decisions.
// ---------------------------------------------------------------------------
void _tsav_fq_boot(string scalar ystar, real scalar tt, real scalar nboot,
    string scalar corecmd, string scalar rmatname, real scalar ntau,
    string scalar BSname, string scalar tnobsname, real scalar seed)
{
    external real colvector __tsav_mu
    real colvector idx, mus, yd, col, tnobs
    real matrix B, R, BS
    real scalar r, s2, j, nb, i1, i5, i10, j1, j5, j10

    if (seed > 0) {
        rseed(seed)
    }
    tnobs = st_matrix(tnobsname)
    B = J(nboot, ntau, .)
    for (r = 1; r <= nboot; r++) {
        idx = ceil(tt :* runiform(tt, 1))
        for (s2 = 1; s2 <= tt; s2++) {
            if (idx[s2] < 1) {
                idx[s2] = 1
            }
        }
        mus = __tsav_mu[idx]
        yd = runningsum(mus)
        st_store(., ystar, yd)
        stata(corecmd)
        R = st_matrix(rmatname)
        B[r, .] = R[., 5]'
    }
    BS = J(ntau, 7, .)
    nb = nboot
    // source convention: index = alpha*Nboot, GAUSS truncation
    j1 = trunc(0.99*nb)
    j5 = trunc(0.95*nb)
    j10 = trunc(0.90*nb)
    if (j1 < 1) {
        j1 = 1
    }
    if (j5 < 1) {
        j5 = 1
    }
    if (j10 < 1) {
        j10 = 1
    }
    i1 = trunc(0.01*nb)
    i5 = trunc(0.05*nb)
    i10 = trunc(0.10*nb)
    if (i1 < 1) {
        i1 = 1
    }
    if (i5 < 1) {
        i5 = 1
    }
    if (i10 < 1) {
        i10 = 1
    }
    for (j = 1; j <= ntau; j++) {
        col = sort(B[., j], 1)
        BS[j, 1] = col[i1]
        BS[j, 2] = col[i5]
        BS[j, 3] = col[i10]
        BS[j, 4] = col[j1]
        BS[j, 5] = col[j5]
        BS[j, 6] = col[j10]
        BS[j, 7] = sum(col :<= tnobs[j])/nb
    }
    st_matrix(BSname, BS)
}

// ---------------------------------------------------------------------------
// Narayan & Popp (2010) two-break test (GAUSS ADF_2breaks, exact)
// model 1 = M1 (breaks in level), 2 = M2 (breaks in level and slope)
// NOTE the source quirks reproduced here:
//   - both models include a linear trend among the deterministics
//   - the deterministic terms enter LAGGED (z_{t-1})
//   - default trimming 0.10 (dynargsGet default, not the 0.15 of the header)
//   - T1/t1 and T2/t2 are the same (case-insensitive) GAUSS symbols:
//     the effective bounds are T1 = max(3+pmax, ceil(trim*T)) (then pmax+3 if
//     < pmax+2) and T2 = min(T-3-pmax, floor((1-trim)*T))
// ---------------------------------------------------------------------------
void _tsav_np(string scalar yname, string scalar touse, real scalar model,
    real scalar pmax, real scalar ic, real scalar trimm)
{
    real colvector y, dy, y1, dep, bb, ee, sevec, taup, aicp, sicp, tstatp
    real colvector du1, du2, dt1v, dt2v, dc, dtv, cv
    real matrix lmat, z, z1, X, XX
    real scalar T, j, T1, T2, tb1, tb2, tb2s, p, lo, n1, kx, pidx, stat
    real scalar ADFmin, tb1min, tb2min, optlag

    y = st_data(., yname, touse)
    T = rows(y)
    dy = y[2::T] - y[1::(T-1)]
    y1 = y[1::(T-1)]
    lmat = J(T-1, pmax, 0)
    for (j = 1; j <= pmax; j++) {
        if (T-1-j >= 1) {
            lmat[(j+1)::(T-1), j] = dy[1::(T-1-j)]
        }
    }
    T1 = max((3 + pmax, ceil(trimm*T)))
    T2 = min((T - 3 - pmax, floor((1 - trimm)*T)))
    if (T1 < pmax + 2) {
        T1 = pmax + 3
    }
    dc = J(T, 1, 1)
    dtv = (1::T)
    ADFmin = 1000
    tb1min = 0
    tb2min = 0
    optlag = 1
    taup = J(pmax+1, 1, .)
    aicp = J(pmax+1, 1, .)
    sicp = J(pmax+1, 1, .)
    tstatp = J(pmax+1, 1, .)
    for (tb1 = T1; tb1 <= T2; tb1++) {
        if (model == 1) {
            tb2s = tb1 + 2
        }
        else {
            tb2s = tb1 + 3
        }
        for (tb2 = tb2s; tb2 <= T2; tb2++) {
            du1 = J(tb1, 1, 0) \ J(T-tb1, 1, 1)
            du2 = J(tb2, 1, 0) \ J(T-tb2, 1, 1)
            if (model == 1) {
                z = dc, dtv, du1, du2
            }
            else {
                dt1v = J(tb1, 1, 0) \ (1::(T-tb1))
                dt2v = J(tb2, 1, 0) \ (1::(T-tb2))
                z = dc, dtv, du1, du2, dt1v, dt2v
            }
            // z1 = lagged deterministics: z_{t-1} aligned with dy_t
            z1 = z[1::(T-1), .]
            for (p = 0; p <= pmax; p++) {
                lo = p + 2
                dep = dy[lo::(T-1)]
                X = y1[lo::(T-1)], z1[lo::(T-1), .]
                if (p > 0) {
                    X = X, lmat[lo::(T-1), 1::p]
                }
                n1 = rows(dep)
                kx = cols(X)
                XX = invsym(cross(X, X))
                bb = XX*cross(X, dep)
                ee = dep - X*bb
                sevec = sqrt(diagonal(XX)*(cross(ee, ee)/(n1-kx)))
                taup[p+1] = bb[1]/sevec[1]
                aicp[p+1] = ln(cross(ee, ee)/n1) + 2*(kx+2)/n1
                sicp[p+1] = ln(cross(ee, ee)/n1) + (kx+2)*ln(n1)/n1
                tstatp[p+1] = abs(bb[kx]/sevec[kx])
            }
            pidx = _tsav_getlagidx(ic, aicp, sicp, tstatp)
            stat = taup[pidx]
            if (stat < ADFmin) {
                tb1min = tb1
                tb2min = tb2
                ADFmin = stat
                optlag = pidx
            }
        }
    }
    cv = _tsav_np_cv(T, model)
    st_numscalar("r(npstat)", ADFmin)
    st_numscalar("r(tb1)", tb1min)
    st_numscalar("r(tb2)", tb2min)
    st_numscalar("r(nplag)", optlag - 1)
    st_numscalar("r(cv1)", cv[1])
    st_numscalar("r(cv5)", cv[2])
    st_numscalar("r(cv10)", cv[3])
}

// Narayan & Popp (2010), Table 3 critical values by sample size
real colvector _tsav_np_cv(real scalar T, real scalar model)
{
    real colvector cv

    if (model == 1) {
        if (T <= 50) {
            cv = (-5.259 \ -4.514 \ -4.143)
        }
        else if (T <= 200) {
            cv = (-4.958 \ -4.316 \ -3.980)
        }
        else if (T <= 400) {
            cv = (-4.731 \ -4.136 \ -3.825)
        }
        else {
            cv = (-4.672 \ -4.081 \ -3.772)
        }
    }
    else {
        if (T <= 50) {
            cv = (-5.949 \ -5.181 \ -4.789)
        }
        else if (T <= 200) {
            cv = (-5.576 \ -4.937 \ -4.596)
        }
        else if (T <= 400) {
            cv = (-5.318 \ -4.741 \ -4.430)
        }
        else {
            cv = (-5.287 \ -4.692 \ -4.396)
        }
    }
    return(cv)
}

// ---------------------------------------------------------------------------
// Carrion-i-Silvestre et al. (2009): response surface for c_bar
// (GAUSS __sbur_c_bar_rs), lam is a 5x1 vector of break fractions
// ---------------------------------------------------------------------------
real scalar _tsav_cbar(real colvector lam)
{
    real rowvector x
    real colvector prm
    real scalar i, j, idx

    x = J(1, 61, .)
    x[1] = 1
    for (i = 1; i <= 5; i++) {
        x[1+i] = lam[i]
        x[6+i] = lam[i]^2
        x[11+i] = lam[i]^3
        x[16+i] = lam[i]^4
    }
    idx = 21
    for (i = 1; i <= 4; i++) {
        for (j = i+1; j <= 5; j++) {
            idx = idx + 1
            x[idx] = abs(lam[i] - lam[j])
            x[idx+10] = abs(lam[i] - lam[j])^2
            x[idx+20] = abs(lam[i] - lam[j])^3
            x[idx+30] = abs(lam[i] - lam[j])^4
        }
    }
    prm = (-13.12832 \ -36.53045 \ 0 \ 20.2423 \ -4.596202 \ -10.31678 \
        115.2092 \ -29.18712 \ -68.36453 \ 5.873121 \ 0 \ -130.337 \
        74.64396 \ 85.48737 \ 0 \ 0 \ 51.98117 \ -53.03452 \ -36.27221 \
        0 \ 11.27727 \ -23.39517 \ -5.360149 \ 23.99683 \ 4.788676 \
        -27.10002 \ -35.78388 \ 51.12371 \ -29.8518 \ -3.069174 \
        -37.45898 \ 64.95842 \ 5.825729 \ -88.78176 \ -11.54197 \
        83.48645 \ 125.2349 \ -173.1259 \ 80.95821 \ 2.863782 \ 118.2829 \
        -80.1287 \ 0 \ 128.872 \ 6.387147 \ -118.1043 \ -199.0615 \
        247.6469 \ -98.05947 \ 0 \ -160.5713 \ 38.52177 \ 0 \ -65.21576 \
        0 \ 62.86494 \ 117.9976 \ -127.5544 \ 46.2304 \ 0 \ 79.1693)
    return(x*prm)
}

// ---------------------------------------------------------------------------
// CiS et al. (2009): response surfaces for the critical values
// (GAUSS pd_msbur_rsf). Returns a 4x3 matrix, rows MSB / MZA / MZT / PT,
// columns 1% / 5% / 10% (the source keeps columns 1, 3, 4 of each 4-block).
// ---------------------------------------------------------------------------
real matrix _tsav_ciscv(real colvector lam, real scalar cbar)
{
    real rowvector x, crit
    real matrix P
    real scalar i, j, idx, d

    x = J(1, 63, .)
    x[1] = 1
    for (i = 1; i <= 5; i++) {
        x[1+i] = lam[i]
        x[6+i] = lam[i]^2
        x[12+i] = lam[i]*cbar
        x[18+i] = lam[i]*cbar^2
    }
    x[12] = cbar
    x[18] = cbar^2
    idx = 23
    for (i = 1; i <= 4; i++) {
        for (j = i+1; j <= 5; j++) {
            idx = idx + 1
            d = abs(lam[i] - lam[j])
            x[idx] = d*cbar
            x[idx+10] = d^2*cbar
            x[idx+20] = d^3*cbar
            x[idx+30] = d^4*cbar
        }
    }
    P = _tsav_cisparam()
    crit = x*P
    return((crit[1], crit[3], crit[4] \
            crit[5], crit[7], crit[8] \
            crit[9], crit[11], crit[12] \
            crit[13], crit[15], crit[16]))
}

real matrix _tsav_cisparam()
{
    real matrix P

    P = J(63, 16, 0)
    P[1,.] = (0.206065483, 0.247173646, 0.279911696, 0.311573002, -26.31391813, -20.61149374, -12.1438623, -6.08490852, -2.52133657, -1.766570893, -1.46435731, -1.277987954, -3.518835863, -3.305558261, -3.454833615, -3.240058047)
    P[2,.] = (-0.131592168, -0.083176707, -0.079273217, -0.136364352, -129.5317914, -84.29286654, -36.47970616, -31.60984523, -6.668037145, -3.349004828, -3.066463141, -3.217982311, -15.69764073, -15.89838295, -10.46560768, -18.14173976)
    P[3,.] = (-0.018230144, 0, 0, 0, -3.503797177, 0, 0, 0, -0.193154126, 0, 0, 0, 2.698367477, 0, 0, 3.094894401)
    P[4,.] = (-0.001829617, 0, 0.036867994, 0, 0, 0, 31.56762014, 13.63038899, 0, 0, 1.32634005, 0.633301965, 6.412055579, 0, 5.223542808, 4.405653332)
    P[5,.] = (-0.071694008, -0.069876819, -0.098386033, -0.063992057, -22.82788603, -13.10518388, -39.87684614, -23.51617143, -2.651936734, -1.893827718, -2.317691146, -1.235047933, -10.32717062, -1.986102341, -8.649758634, -1.987760716)
    P[6,.] = (-0.113224418, -0.123618939, -0.114531349, -0.171308084, -71.88919188, -56.62886058, -22.2258743, -45.74522263, -4.308941432, -3.028641495, -2.666831955, -3.513398467, -5.232680619, -13.76324196, -6.17789082, -12.35876925)
    P[7,.] = (0.045497638, 0.034139777, 0.033789576, 0.055316264, 47.47172982, 32.22798967, 18.27175882, 11.56516258, 2.242409839, 1.205414124, 1.172595231, 1.065619758, 4.664094355, 5.198014133, 4.311973301, 6.869459093)
    P[8,.] = (0.005667139, 0.007183722, 0.014895671, 0.00832245, 0, 5.942060926, 8.371062059, 6.269767228, 0, 0.307559368, 0.499975826, 0.381693852, 0, 0, 1.28702057, 0)
    P[9,.] = (0, 0, 0, 0.007275262, 0, 0, 0, 0, -0.085060207, 0, 0, 0.194838181, 0, 0, 0, 0)
    P[10,.] = (0.011393725, 0, 0.006925649, 0.007101886, 14.70004959, 11.44618918, 8.81485526, 8.184740859, 0.87936997, 0.56718476, 0.473032483, 0.402464154, 1.183189107, 1.643475816, 0, 1.791921929)
    P[11,.] = (0.041416456, 0.037458345, 0.036987117, 0.053884486, 37.54355821, 29.58007215, 20.19100032, 26.06901605, 2.110275336, 1.487278396, 1.44212234, 1.676575305, 4.970060977, 6.400077634, 6.119151722, 6.360822366)
    P[12,.] = (0.006744983, 0.009229194, 0.01117728, 0.012830135, 0, 0, 0.383835339, 0.85365241, 0.114546234, 0.141951363, 0.144171363, 0.139480635, -0.544592947, -0.53171116, -0.666577377, -0.74225554)
    P[13,.] = (-0.008881686, -0.005163807, -0.00419873, -0.007623439, -8.508298137, -5.471425562, -1.873513518, -2.233163794, -0.457380818, -0.223355446, -0.186935755, -0.219966762, -1.116029433, -1.044782065, -0.624407428, -1.106702061)
    P[14,.] = (-0.000477302, 0.001411259, 0.001490823, 0.000992921, 0, 0.821678805, 0.839344685, 0.693315109, 0, 0.040353055, 0.049102085, 0.044657513, 0.3017767, 0.013704832, 0.128412682, 0.273588757)
    P[15,.] = (0, 0, 0.003166102, 0.000383907, 0.315654917, 0.353092298, 2.84583848, 1.121511646, 0, 0.017148054, 0.120078208, 0.072095127, 0.636962573, 0.022144415, 0.523175329, 0.412172316)
    P[16,.] = (-0.004879334, -0.005905712, -0.007875096, -0.005118819, -0.270010842, 0, -2.545569579, -1.407264853, -0.12770179, -0.112106002, -0.155499347, -0.077321047, -0.756479306, 0, -0.693933651, 0)
    P[17,.] = (-0.005109823, -0.006312857, -0.005745887, -0.008799213, -2.346920289, -1.877066153, 0, -1.439270423, -0.15405538, -0.106423683, -0.083556342, -0.135095954, 0, -0.587410296, 0, -0.447876179)
    P[18,.] = (0.000147113, 0.000180001, 0.000216343, 0.000261777, 0.010036211, 0, 0, 0.019804635, 0.003272583, 0.002943926, 0.002822847, 0.003038937, 0, 0.004379253, 0, 0)
    P[19,.] = (-0.00018911, -0.000109989, -7.2476E-05, -0.000130686, -0.180808958, -0.117234577, -0.036217559, -0.050234142, -0.009870463, -0.004857975, -0.003790936, -0.00464363, -0.023994328, -0.021043324, -0.012868843, -0.021946248)
    P[20,.] = (0, 4.98993E-05, 3.87592E-05, 2.10757E-05, 0, 0.023326443, 0.022083426, 0.01757633, 0, 0.001151668, 0.001220296, 0.001083382, 0.007411666, 0, 0.003373215, 0.005094191)
    P[21,.] = (0, 0, 6.56928E-05, 0, 0.012434324, 0.014488305, 0.063997488, 0.02118732, 0, 0.000643625, 0.002670921, 0.001514238, 0.014861192, 0, 0.01257116, 0.009431941)
    P[22,.] = (-0.000101491, -0.000129658, -0.000170948, -0.000117271, 0, 0, -0.053703036, -0.033337638, -0.002353799, -0.002482217, -0.003375858, -0.001867064, -0.016239629, 0, -0.014642206, 0)
    P[23,.] = (-9.12367E-05, -0.000115959, -0.000107388, -0.000164942, -0.042581947, -0.034973534, 0, -0.028709714, -0.002863078, -0.002011657, -0.001581567, -0.002664683, 0, -0.01272247, -0.000903779, -0.009291186)
    P[24,.] = (0.001502737, 0.001632669, 0.001702841, 0.002021723, 1.632727956, 1.605877644, 1.004404184, 0.868853384, 0.083333993, 0.079478462, 0.065367552, 0.063715851, 0.310198878, 0.360029916, 0.330191659, 0.381965256)
    P[25,.] = (0, 0, 0, 0.001328702, 0, 0, 0.202483643, 0.620890565, 0.006479305, 0, 0.006933813, 0.032275901, 0, 0.134517084, 0, 0)
    P[26,.] = (0, 0, 0, 0, -0.192709346, 0, 0, 0.36221429, -0.021488357, 0, -0.001200112, 0.011828416, 0, 0.099982238, 0.103069901, 0.084964524)
    P[27,.] = (0.000680309, 0, 0.000414141, 0, 0.499555624, 0, 0.269613377, 0.236974663, 0.040862054, 0.022834275, 0.016456528, 0.017878489, 0.127583576, 0, 0.129658797, 0.079531443)
    P[28,.] = (0.001076323, 0.000746542, 0.001534534, 0.001854214, 1.222722601, 0.980994239, 0.952309602, 0.981692971, 0.06103088, 0.048322535, 0.062848187, 0.068704135, 0.34819102, 0.351457151, 0.285248474, 0.389323838)
    P[29,.] = (0.001602307, 0.001718432, 0.000513486, 0.001320012, 1.284000309, 0.790608102, 0.29128284, 0.192243305, 0.073187675, 0.037307962, 0.015603486, 0.008831525, 0, 0.161831707, 0.161667434, 0)
    P[30,.] = (-3.13393E-05, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.178056784)
    P[31,.] = (0.000915733, 0.000465616, 0.000828461, 0.0015141, 0.895732659, 0.635112942, 0.357984838, 0.70758854, 0.051259672, 0.033026078, 0.039391705, 0.054005928, 0.114186355, 0.272130694, 0.173459972, 0.236092255)
    P[32,.] = (0.001297134, 0.001989702, 0.001648169, 0.000516108, 1.486377973, 1.448576409, 1.438138712, 0.44934008, 0.06084729, 0.063368696, 0.060177344, 0.023938672, 0.184586257, 0.01941423, 0.270773963, 0.306240705)
    P[33,.] = (0.000779578, 0.000538182, 0.001171215, 0.001602482, 0.553062822, 0.350422194, 0.468555532, 0.425359716, 0.03353265, 0.028507047, 0.042001183, 0.043941342, 0.138176188, 0.167194529, 0.177010468, 0.210812736)
    P[34,.] = (-0.004345646, -0.005825265, -0.006868541, -0.00796319, -4.479282599, -4.8858522, -3.136191913, -2.458114748, -0.231523802, -0.248286094, -0.213133815, -0.194637365, -0.922468402, -0.985366571, -0.923999903, -1.130047)
    P[35,.] = (0, 0, -5.04974E-05, -0.004738169, 0, 0, -0.403826354, -1.909565827, -0.009229909, 0, -0.008242746, -0.098804196, 0, -0.507406066, 0, -0.021572133)
    P[36,.] = (0, 0, 0, -5.05847E-05, 0, 0, -0.417036823, -1.119339897, 0.043558418, 0, 0, -0.02166457, 0, -0.454778552, -0.405307294, -0.332490735)
    P[37,.] = (-0.004679147, -0.002222758, -0.004695542, -0.003025465, -3.000000235, -0.964674073, -1.860955551, -1.416293256, -0.208817124, -0.134046896, -0.136841302, -0.125359252, -0.631341206, 0, -0.587433653, -0.487244954)
    P[38,.] = (-0.002164419, -0.001967051, -0.005817976, -0.006214813, -2.092813644, -1.850008801, -2.941485078, -2.857336306, -0.104657956, -0.097680824, -0.195504402, -0.201534614, -1.128138371, -0.979191193, -0.703716053, -1.028065809)
    P[39,.] = (-0.005484897, -0.006163019, -0.00124581, -0.004023965, -4.057652598, -2.738110318, -0.61059659, -0.332618017, -0.222016788, -0.126749269, -0.030055771, -0.009635946, 0, -0.599177023, -0.584313329, 0)
    P[40,.] = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.089011851, 0.124842286, 0, -0.508722695)
    P[41,.] = (-0.00165207, -0.001163458, -0.003198713, -0.005190305, -1.276113948, -0.889271056, -0.747636053, -1.839058302, -0.069225581, -0.055577102, -0.112927132, -0.159646685, -0.154761007, -0.597212331, -0.288535459, -0.379766195)
    P[42,.] = (-0.003789148, -0.0064054, -0.004817945, -0.000649124, -4.496633229, -4.594051759, -4.597029802, -0.694146329, -0.167498301, -0.194612003, -0.171395671, -0.037729237, -0.67125797, 0, -0.997655245, -1.102797825)
    P[43,.] = (-0.001096697, -0.000571708, -0.00416902, -0.004537825, 0, 0, -1.150001536, -0.397475842, 0, -0.036557664, -0.117430944, -0.100061587, -0.189720447, -0.156458132, -0.171258066, -0.367949029)
    P[44,.] = (0.004725461, 0.008280992, 0.010451449, 0.011883757, 4.700545573, 6.27119741, 4.178380993, 3.094835252, 0.247801364, 0.322844187, 0.291165672, 0.255602766, 1.114026174, 1.079445866, 1.116574652, 1.435226881)
    P[45,.] = (0, 0, 0, 0.006868754, -0.046355934, 0, 0.220249151, 2.5334366, 0, 0, 0, 0.133924549, 0, 0.776781333, 0, 0.023477765)
    P[46,.] = (0, 0, 0, 0, 0.811913273, 0, 0.998239994, 1.47110097, -0.023485383, 0, 0, 0.010896637, 0, 0.75874924, 0.657353781, 0.52514633)
    P[47,.] = (0.009039167, 0.005334347, 0.009880517, 0.007579264, 5.503001956, 2.327899451, 3.493280873, 2.574956682, 0.361357325, 0.240732181, 0.272754543, 0.241334704, 1.095018463, 0, 0.92651935, 0.916080571)
    P[48,.] = (0.001237718, 0.001401401, 0.008484588, 0.008416308, 1.02572597, 1.046754093, 3.967072173, 3.760263385, 0.052317695, 0.060133444, 0.261308977, 0.263752707, 1.565059276, 1.23081264, 0.745421462, 1.252746367)
    P[49,.] = (0.008139359, 0.009159069, 0.000918836, 0.005481258, 5.865893641, 4.234610015, 0.400719375, 0.166727218, 0.311141107, 0.193402656, 0.018125688, 0, 0, 0.979632255, 0.860652087, 0)
    P[50,.] = (0, -5.0721E-05, 0, -4.79068E-05, -0.038957002, 0, 0, -0.024479515, -0.002278556, 0, 0, -0.001680802, -0.216616021, -0.359131022, 0, 0.639272615)
    P[51,.] = (0.000809364, 0.000809675, 0.004538856, 0.007550505, 0, 0, 0.474813133, 2.435657851, 0, 0.0278107, 0.141916393, 0.221121513, 0, 0.592286581, 0.142276904, 0.175199011)
    P[52,.] = (0.004893033, 0.00886491, 0.006136727, 0, 6.088534864, 6.377018844, 6.286955077, 0.283095346, 0.216868391, 0.263466015, 0.214566444, 0.01661417, 1.025111042, 0, 1.492001945, 1.627703373)
    P[53,.] = (0.000352955, 0, 0.005980789, 0.005904448, -1.299540015, -0.727894226, 1.347458203, 0, -0.07801745, 0.011508158, 0.152800184, 0.120073816, 0.059252271, 0, 0, 0.316314136)
    P[54,.] = (-0.00179525, -0.004086714, -0.005417723, -0.006105744, -1.709027953, -2.961023806, -2.016300404, -1.471404534, -0.092721304, -0.151065777, -0.142110387, -0.122864487, -0.497264715, -0.42244023, -0.500157415, -0.683828054)
    P[55,.] = (-0.000109346, -5.35444E-05, 0, -0.003596853, 0, 0, 0, -1.260438948, 0, 0, 0, -0.069784359, -0.009288362, -0.426512724, 0, 0)
    P[56,.] = (0, 0, -5.57416E-05, 0, -0.674459313, 0, -0.621350492, -0.713530106, 0, 0, 0, 0, 0, -0.427537109, -0.374390888, -0.288974182)
    P[57,.] = (-0.005138561, -0.003190206, -0.00580564, -0.004742335, -3.034884695, -1.395732499, -1.970609696, -1.436365327, -0.195891984, -0.132984979, -0.157536611, -0.138135224, -0.602410686, -0.011191895, -0.488796736, -0.530722743)
    P[58,.] = (0, 0, -0.004283777, -0.004053352, 0, 0, -1.984788663, -1.891122229, 0, 0, -0.126796139, -0.128361041, -0.785227631, -0.594785802, -0.300198458, -0.603224179)
    P[59,.] = (-0.004463549, -0.004931673, 0, -0.002862951, -3.272972482, -2.462761762, 0, 0, -0.170881185, -0.111894795, 0, 0, -0.010121395, -0.588066534, -0.443113015, 0)
    P[60,.] = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.132384584, 0.271088476, 0, -0.311195596)
    P[61,.] = (0, 0, -0.002139895, -0.003923755, 0.528575644, 0.366450912, 0, -1.31877996, 0.026385386, 0, -0.064682374, -0.114348818, 0.067265656, -0.247275032, 0, 0)
    P[62,.] = (-0.00241816, -0.004474527, -0.00297588, 0.00023128, -3.169547182, -3.321085197, -3.169497865, 0, -0.1144538, -0.135673769, -0.103644277, 0, -0.570723203, -0.03959037, -0.788709363, -0.849343883)
    P[63,.] = (0, 0, -0.002984229, -0.002989047, 0.798089652, 0.412611617, -0.64630145, 0, 0.048845232, 0, -0.074881653, -0.061526855, 0, 0, 0, -0.147843764)
    return(P)
}

// ---------------------------------------------------------------------------
// GLS detrending (GAUSS __sbur_glsd): returns only the quasi-differenced SSR
// (fast path used inside the break-search loops)
// ---------------------------------------------------------------------------
real scalar _tsav_gls_ssr(real colvector y, real matrix z, real scalar cbar)
{
    real scalar T, abar
    real colvector ya, e, bhat
    real matrix za

    T = rows(y)
    abar = 1 + cbar/T
    ya = y[1] \ (y[2::T] - abar*y[1::(T-1)])
    za = z[1, .] \ (z[2::T, .] - abar*z[1::(T-1), .])
    bhat = invsym(cross(za, za))*cross(za, ya)
    e = ya - za*bhat
    return(cross(e, e))
}

// GLS detrending returning the detrended series (yt = y - z*bhat)
real colvector _tsav_glsd_yt(real colvector y, real matrix z, real scalar cbar)
{
    real scalar T, abar
    real colvector ya, bhat
    real matrix za

    T = rows(y)
    abar = 1 + cbar/T
    ya = y[1] \ (y[2::T] - abar*y[1::(T-1)])
    za = z[1, .] \ (z[2::T, .] - abar*z[1::(T-1), .])
    bhat = invsym(cross(za, za))*cross(za, ya)
    return(y - z*bhat)
}

// GLS coefficient vector (for the fitted broken trend in the graph)
real colvector _tsav_glsd_b(real colvector y, real matrix z, real scalar cbar)
{
    real scalar T, abar
    real colvector ya
    real matrix za

    T = rows(y)
    abar = 1 + cbar/T
    ya = y[1] \ (y[2::T] - abar*y[1::(T-1)])
    za = z[1, .] \ (z[2::T, .] - abar*z[1::(T-1), .])
    return(invsym(cross(za, za))*cross(za, ya))
}

// ---------------------------------------------------------------------------
// GAUSS __sbur_s2ar: lag choice for the long-run variance on OLS-detrended
// data, MAIC (penalty=0) or BIC (penalty=1)
// ---------------------------------------------------------------------------
real scalar _tsav_s2ar(real colvector yts, real scalar penalty,
    real scalar kmax, real scalar kmin)
{
    real scalar T, nef, k, j, kopt, sumy
    real colvector dyf, dep, b, e, s2e, tauv, mic, kk
    real matrix reg, regk

    T = rows(yts)
    dyf = yts[2::T] - yts[1::(T-1)]
    // dep and regressors on the common sample t = kmax+2 .. T
    dep = dyf[(kmax+1)::(T-1)]
    reg = yts[(kmax+1)::(T-1)]
    for (j = 1; j <= kmax; j++) {
        reg = reg, dyf[(kmax+1-j)::(T-1-j)]
    }
    nef = T - kmax - 1
    sumy = cross(reg[., 1], reg[., 1])
    s2e = J(kmax+1, 1, 999)
    tauv = J(kmax+1, 1, 0)
    for (k = kmin; k <= kmax; k++) {
        regk = reg[., 1::(k+1)]
        b = invsym(cross(regk, regk))*cross(regk, dep)
        e = dep - regk*b
        s2e[k+1] = cross(e, e)/nef
        tauv[k+1] = (b[1]*b[1])*sumy/s2e[k+1]
    }
    kk = (0::kmax)
    if (penalty == 0) {
        mic = ln(s2e) + 2:*(kk + tauv):/nef
    }
    else {
        mic = ln(s2e) + ln(nef):*kk:/nef
    }
    kopt = 1
    for (j = 2; j <= kmax+1; j++) {
        if (mic[j] < mic[kopt]) {
            kopt = j
        }
    }
    return(kopt - 1)
}

// ---------------------------------------------------------------------------
// GAUSS __sbur_adfp: ADF regression on the GLS-detrended series, no constant
// returns (adf t-stat, alpha-hat, autoregressive long-run variance)
// ---------------------------------------------------------------------------
real rowvector _tsav_adfp(real colvector yt, real scalar kstar)
{
    real scalar T, nef, s2e, sre, adf, sumb, j
    real colvector dyf, dep, rho, e
    real matrix reg, XX

    T = rows(yt)
    dyf = yt[2::T] - yt[1::(T-1)]
    dep = dyf[(kstar+1)::(T-1)]
    reg = yt[(kstar+1)::(T-1)]
    for (j = 1; j <= kstar; j++) {
        reg = reg, dyf[(kstar+1-j)::(T-1-j)]
    }
    XX = invsym(cross(reg, reg))
    rho = XX*cross(reg, dep)
    e = dep - reg*rho
    nef = rows(dep)
    s2e = cross(e, e)/nef
    sre = XX[1, 1]*s2e
    adf = rho[1]/sqrt(sre)
    if (kstar > 0) {
        sumb = sum(rho[2::(kstar+1)])
    }
    else {
        sumb = 0
    }
    return((adf, rho[1] + 1, s2e/(1 - sumb)^2))
}

// build the deterministic matrix z for the CiS models with breaks
real matrix _tsav_cis_z(real scalar T, real scalar model, real colvector tb)
{
    real matrix z
    real colvector du, dtj
    real scalar j, m

    z = J(T, 1, 1), (1::T)
    m = rows(tb)
    for (j = 1; j <= m; j++) {
        dtj = J(tb[j], 1, 0) \ (1::(T-tb[j]))
        if (model == 3) {
            du = J(tb[j], 1, 0) \ J(T-tb[j], 1, 1)
            z = z, du, dtj
        }
        else {
            z = z, dtj
        }
    }
    return(z)
}

// pad break fractions to a 5x1 vector
real colvector _tsav_lam5(real colvector tb, real scalar T)
{
    real colvector lam
    real scalar m

    m = rows(tb)
    lam = tb :/ T
    if (m < 5) {
        lam = lam \ J(5-m, 1, 0)
    }
    return(lam)
}

// ---------------------------------------------------------------------------
// CiS main routine (GAUSS sbur_gls + __sbur_multiple_gls_brute, exact)
// model: 0 const / 1 trend / 2 slope breaks / 3 level+slope breaks
// known = 1 with tbmat holding known break positions; else brute search
// Fills: STname (1x7: pt mpt adf za mza msb mzt), CVname (4x3), TBname (mx1)
//        r(cbar), r(krule); stores the fitted GLS trend in fitname
// ---------------------------------------------------------------------------
void _tsav_cis(string scalar yname, string scalar touse, real scalar model,
    real scalar known, string scalar tbmat, real scalar nbrk,
    real scalar penalty, real scalar kmax, real scalar kmin,
    string scalar fitname, string scalar STname, string scalar CVname,
    string scalar TBname)
{
    real colvector y, mintb, lam, yt, ydols, bo, bgls, tbc, adfrow
    real matrix z, zc, CV
    real scalar T, cbar, cb, minssra, ssra, j, jj, jjj
    real scalar ahat, s2u, sumyt2, krule, adf, sar, bt, za, mza, msb, mzt
    real scalar ssr1, pt, mpt, ssrafin
    real colvector x1, r

    y = st_data(., yname, touse)
    T = rows(y)
    mintb = J(1, 1, 0)

    if (model == 0) {
        z = J(T, 1, 1)
        cbar = -7
    }
    else if (model == 1) {
        z = J(T, 1, 1), (1::T)
        cbar = -13.5
    }
    else {
        if (known == 1) {
            mintb = st_matrix(tbmat)
            if (cols(mintb) > 1) {
                mintb = mintb'
            }
            z = _tsav_cis_z(T, model, mintb)
            cbar = _tsav_cbar(_tsav_lam5(mintb, T))
        }
        else {
            minssra = cross(y, y)
            if (nbrk == 1) {
                for (j = 3; j <= T-3; j++) {
                    tbc = (j)
                    zc = _tsav_cis_z(T, model, tbc)
                    cb = _tsav_cbar(_tsav_lam5(tbc, T))
                    ssra = _tsav_gls_ssr(y, zc, cb)
                    if (ssra < minssra) {
                        mintb = tbc
                        minssra = ssra
                    }
                }
            }
            else if (nbrk == 2) {
                for (j = 3; j <= T-3-2; j++) {
                    for (jj = j+2; jj <= T-3; jj++) {
                        tbc = (j \ jj)
                        zc = _tsav_cis_z(T, model, tbc)
                        cb = _tsav_cbar(_tsav_lam5(tbc, T))
                        ssra = _tsav_gls_ssr(y, zc, cb)
                        if (ssra < minssra) {
                            mintb = tbc
                            minssra = ssra
                        }
                    }
                }
            }
            else {
                for (j = 3; j <= T-3-4; j++) {
                    for (jj = j+2; jj <= T-3-2; jj++) {
                        for (jjj = jj+2; jjj <= T-3; jjj++) {
                            tbc = (j \ jj \ jjj)
                            zc = _tsav_cis_z(T, model, tbc)
                            cb = _tsav_cbar(_tsav_lam5(tbc, T))
                            ssra = _tsav_gls_ssr(y, zc, cb)
                            if (ssra < minssra) {
                                mintb = tbc
                                minssra = ssra
                            }
                        }
                    }
                }
            }
            z = _tsav_cis_z(T, model, mintb)
            cbar = _tsav_cbar(_tsav_lam5(mintb, T))
        }
    }

    // final GLS detrending and unit-root statistics (source common tail)
    ssrafin = _tsav_gls_ssr(y, z, cbar)
    yt = _tsav_glsd_yt(y, z, cbar)
    x1 = yt[1::(T-1)]
    ahat = cross(x1, yt[2::T])/cross(x1, x1)
    r = yt[2::T] - ahat*x1
    s2u = cross(r, r)/(T - 1)
    sumyt2 = cross(x1, x1)/((T - 1)^2)

    // long-run variance: lags chosen on OLS-detrended data (Perron's note)
    bo = invsym(cross(z, z))*cross(z, y)
    ydols = y - z*bo
    krule = _tsav_s2ar(ydols, penalty, kmax, kmin)
    adfrow = _tsav_adfp(yt, krule)'
    adf = adfrow[1]
    sar = adfrow[3]

    bt = T - 1
    za = bt*(ahat - 1) - (sar - s2u)/(2*sumyt2)
    mza = ((yt[T]^2)/bt - sar)/(2*sumyt2)
    msb = sqrt(sumyt2/sar)
    mzt = mza*msb

    // PT and MPT
    ssr1 = _tsav_gls_ssr(y, z, 0)
    pt = (ssrafin - (1 + cbar/T)*ssr1)/sar
    if (model == 0) {
        mpt = (cbar*cbar*sumyt2 - cbar*((yt[T]^2)/T))/sar
    }
    else {
        mpt = (cbar*cbar*sumyt2 + (1 - cbar)*((yt[T]^2)/T))/sar
    }

    // critical values from the response surfaces
    if (model <= 1) {
        lam = J(5, 1, 0)
    }
    else {
        lam = _tsav_lam5(mintb, T)
    }
    CV = _tsav_ciscv(lam, cbar)

    // fitted broken trend for the graph
    bgls = _tsav_glsd_b(y, z, cbar)
    st_store(., fitname, touse, z*bgls)

    st_matrix(STname, (pt, mpt, adf, za, mza, msb, mzt))
    st_matrix(CVname, CV)
    st_matrix(TBname, mintb)
    st_numscalar("r(cbar)", cbar)
    st_numscalar("r(krule)", krule)
}

end
*==============================================================================
* End of tsadvroot.ado
*==============================================================================