*! ckptest 1.0.0  07jul2026
*! GLS-based unit root tests with multiple structural breaks
*! Carrion-i-Silvestre, Kim, and Perron (2009, Econometric Theory 25, 1754-1792)
*! Line-by-line Stata/Mata port of the original GAUSS code msbur.src by
*! Josep Lluis Carrion-i-Silvestre (based on Ng and Perron 2001 code),
*! last GAUSS modification July 1st 2006, distributed with CKP (2009).
*! Both estimation paths of the source are implemented:
*!   method(dp)    = sbur_multiple_gls_algorithm (estimation = 1, default)
*!   method(brute) = sbur_multiple_gls_brute     (estimation = 0)
*! Single documented deviation from the source: in the unknown-break brute
*! search min_ssra is initialized with maxdouble() instead of y'y; this is
*! outcome-identical whenever the GAUSS code runs and avoids its singular
*! matrix crash on strongly stationary series (see the Mata comment).
*! Author of the Stata port: Ozan Eruygur, AHBV University
*! GAUSS-to-Mata primitive mapping used throughout:
*!   inv()    -> luinv()      (GAUSS general LU inverse)
*!   invpd()  -> cholinv()    (GAUSS positive definite inverse)
*!   olsqr2() -> qrsolve()    (GAUSS QR least squares)
*!   y/x      -> qrsolve(x,y) (GAUSS matrix division, least squares)
*!   minc/minindc -> min with first-occurrence index helper
*! Option defaults follow the author-distributed usage and tspdlib control
*! structure: model(break), penalty(maic), kmax(4), kmin(0), method(dp),
*! maxiter(100).

program define ckptest, rclass
    version 14.0
    syntax varname(ts) [if] [in] [, Model(string) Breaks(integer 1) BREAKDates(string) METHod(string) MAXiter(integer 100) Penalty(string) KMax(integer 4) KMin(integer 0) noPRint noPREtest PRETESTAll BG BGLags(integer 0) BGDf(integer 20) ]

    local dispname `varlist'

    * ---- model option: 0 const / 1 trend / 2 slope / 3 level+slope ----
    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", "mixed") 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), {bf:slope} (2, breaks in the trend slope) or {bf:break} (3, breaks in level and slope)"
        exit 198
    }

    * ---- penalty option: 0 maic / 1 bic (GAUSS input penalty) ----
    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
    }

    * ---- method option: 1 dp (default) / 0 brute (GAUSS input estimation) ----
    if "`method'" == "" local method "dp"
    local method = lower("`method'")
    if inlist("`method'", "brute", "bruteforce", "0") local est 0
    else if inlist("`method'", "dp", "algorithm", "algo", "dynamic", "1") local est 1
    else {
        di as err "method() must be {bf:dp} (dynamic programming, GAUSS estimation=1, the default) or {bf:brute} (grid search, GAUSS estimation=0)"
        exit 198
    }
    if `maxiter' < 1 {
        di as err "maxiter() must be a positive integer"
        exit 198
    }
    if `kmin' < 0 | `kmax' < `kmin' {
        di as err "require 0 <= kmin() <= kmax()"
        exit 198
    }

    * ---- sample checks (wrapper-level validation, not part of the GAUSS algorithm) ----
    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 "ckptest works on a single time series; data are xtset with panel variable {bf:`r(panelvar)'}"
        exit 459
    }
    local tvar "`r(timevar)'"
    marksample touse
    * resolve time-series operators (D. L. F. S.) into a temporary plain variable
    tsrevar `varlist'
    local yv `r(varlist)'
    markout `touse' `yv'
    qui tsreport if `touse'
    if r(N_gaps) > 0 {
        di as err "the estimation sample contains `r(N_gaps)' gap(s) in `tvar'; unit-root tests require a contiguous series"
        exit 498
    }
    qui count if `touse'
    local T = r(N)
    if `T' < `kmax' + 25 {
        di as err "insufficient observations (T=`T')"
        exit 2001
    }

    * ---- break setup: known dates (model[2]=0) or unknown count (model[2]=1) ----
    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 "the procedure is designed for up to m <= 5 structural breaks"
                exit 198
            }
            matrix `TB' = J(`m', 1, .)
            local tfmt : format `tvar'
            local i 0
            local prev = -1e300
            foreach tok of local breakdates {
                local ++i
                capture local d = `tok'
                if _rc | missing(`"`d'"') {
                    if strpos("`tfmt'", "tq") local d = quarterly("`tok'", "YQ")
                    else if strpos("`tfmt'", "tm") local d = monthly("`tok'", "YM")
                    else if strpos("`tfmt'", "th") local d = halfyearly("`tok'", "YH")
                    else if strpos("`tfmt'", "tw") local d = weekly("`tok'", "YW")
                    else if strpos("`tfmt'", "td") local d = daily("`tok'", "DMY")
                    else local d = .
                }
                if missing(`d') {
                    di as err "breakdates(): could not interpret `tok' as a value of `tvar' or an observation number"
                    exit 198
                }
                * interpretation order: (1) a value of the time variable within the
                * sample; (2) otherwise an integer observation number 1..T (position
                * within the estimation sample). For t = 1..T the two coincide.
                local posval = .
                qui su `pos' if `touse' & float(`tvar') == float(`d'), meanonly
                if r(N) > 0 local posval = r(min)
                else if `d' == int(`d') & `d' >= 1 & `d' <= `T' local posval = `d'
                if missing(`posval') {
                    di as err "breakdates(): `tok' is neither a value of `tvar' in the sample nor an observation number between 1 and `T'"
                    exit 198
                }
                if `posval' <= `prev' {
                    di as err "breakdates() must be strictly increasing"
                    exit 198
                }
                local prev = `posval'
                matrix `TB'[`i', 1] = `posval'
            }
        }
        else {
            if `m' < 1 | `m' > 5 {
                di as err "the procedure is designed for up to m <= 5 structural breaks"
                exit 198
            }
            if `est' == 0 & `m' > 3 {
                di as err "with unknown break dates the brute-force search of the GAUSS source supports breaks(1) to breaks(3) only; use {bf:method(dp)} (the default) for 4 or 5 unknown breaks, or supply known dates via breakdates()"
                exit 198
            }
            if `est' == 1 {
                local h = int(0.10*`T')
                if `h' < 5 | `T' < (`m'+2)*`h' {
                    di as err "sample too short for the dynamic programming (dp) with breaks(`m') and 10 percent trimming"
                    exit 2001
                }
            }
        }
    }
    else {
        local m 0
    }
    if !`known' matrix `TB' = J(1, 1, 0)

    * ---- Perron-Yabu (2009, JBES) pretest for a break in trend ----
    * ---- (run by default; suppressed by the nopretest option) ----
    if "`pretest'" == "" {
        * defaults of qfgls.prg: BIC (criteria=2), eps=0.15, kmax rule;
        * PY model 2 (slope shift) for model(slope), PY model 3 otherwise
        local pym = cond(`mnum' == 2, 2, 3)
        local pycrit 2
        local pyeps = 0.15
        local pyk = int(12*(`T'/100)^0.25)
        if max(int(`pyeps'*`T'), `pyk'+2) > int((1-`pyeps')*`T') {
            di as text "(pretest skipped: insufficient observations)"
            local pretest "nopretest"
        }
    }
    if "`pretest'" == "" {
        tempname PYOUT
        mata: _ckp_py("`yv'", "`touse'", `pym', `pycrit', `pyeps', `pyk', "`PYOUT'")
        local pyexpw = `PYOUT'[1,1]
        local pybp = `PYOUT'[1,2]
        local pycv10 = `PYOUT'[1,3]
        local pycv5 = `PYOUT'[1,4]
        local pycv1 = `PYOUT'[1,5]
        qui su `tvar' if `pos' == `pybp' & `touse', meanonly
        local pybd = r(min)
        local tf2 : format `tvar'
        local pybds = string(`pybd', "`tf2'")
        if "`print'" == "" {
            local pymlab "breaks in the trend slope"
            if `pym' == 3 local pymlab "breaks in level and slope"
            local pystars ""
            if `pyexpw' > `pycv10' local pystars "*"
            if `pyexpw' > `pycv5' local pystars "**"
            if `pyexpw' > `pycv1' local pystars "***"
            di
            di as text "{hline 78}"
            di as text "Pretest for a break in trend (Perron and Yabu 2009, JBES 27, 369-396)"
            di as text "{hline 78}"
            di as text "Variable   : " as result "`dispname'" as text _col(44) "Obs           = " as result %8.0f `T'
            di as text "Model      : " as result "`pymlab'"
            if "`pretestall'" != "" {
                di as text "Lags       : " as result "BIC" as text _col(44) "kmax          = " as result %8.0f `pyk'
                di as text "Trimming   : " as result %4.2f `pyeps' as text _col(44) "Break date    = " as result "`pybds'" as text " (obs `pybp')"
            }
            di as text "{hline 78}"
            di as text "  Exp-W(FS)" as result %14.6f `pyexpw' as text %-4s "`pystars'" as text _col(36) "cv:  10% " as result %6.2f `pycv10' as text "  5% " as result %6.2f `pycv5' as text "  1% " as result %6.2f `pycv1'
            di as text "{hline 78}"
            di as text "H0: no break in the trend function."
            if `pyexpw' > `pycv1' {
                di as text "Conclusion: Reject at the 1% level: evidence of a break in the trend function."
            }
            else if `pyexpw' > `pycv5' {
                di as text "Conclusion: Reject at the 5% level: evidence of a break in the trend function."
            }
            else if `pyexpw' > `pycv10' {
                di as text "Conclusion: Reject at the 10% level: evidence of a break in the trend function."
            }
            else {
                di as text "Conclusion: Do not reject at the 10% level: no evidence of a trend break. With no"
                di as text "break, the break unit root tests can over-reject (CKP 2009, Section 6)."
            }
            di as text "{hline 78}"
        }
    }

    * ---- Breusch-Godfrey setup (mkbg mechanics: zero-filled lagged residuals,
    * ---- no constant, uncentered R2, LM = n*R2 ~ chi2(p)) ----
    if `bgdf' < 0 {
        di as err "bgdf() must be a nonnegative integer"
        exit 198
    }
    local bgsel = ("`bg'" != "")
    if `bglags' < 0 {
        di as err "bglags() must be a positive integer"
        exit 198
    }
    if `bglags' > 0 {
        * explicit horizon: used as given, never capped
        local bgh = `bglags'
    }
    else {
        * automatic horizon: about two years of orders by frequency,
        * capped at 5*floor(4*(T/100)^0.25) as in narayanp
        qui tsset
        local tv `r(timevar)'
        local fmt : format `tv'
        local bgh 2
        if strpos("`fmt'", "%ty") local bgh 2
        if strpos("`fmt'", "%tq") local bgh 8
        if strpos("`fmt'", "%tm") local bgh 24
        if strpos("`fmt'", "%tw") local bgh 52
        if strpos("`fmt'", "%td") local bgh 100
        local bgcap = 5*floor(4*(`T'/100)^0.25)
        if `bgh' > `bgcap' local bgh = `bgcap'
    }
    if `bgh' < 1 local bgh 1
    scalar __ckpbg_sel = `bgsel'
    scalar __ckpbg_h = `bgh'
    scalar __ckpbg_df = `bgdf'

    tempname ST CV TBOUT
    mata: _ckp_main("`yv'", "`touse'", `mnum', `known', "`TB'", `m', `pen', `kmax', `kmin', `est', `maxiter', "`ST'", "`CV'", "`TBOUT'")
    local cbar = r(cbar)
    local krule = r(krule)
    local bgminp = scalar(__ckpbg_minp)
    local bgmino = scalar(__ckpbg_mino)
    local bgwarn = scalar(__ckpbg_warn)
    scalar drop __ckpbg_sel __ckpbg_h __ckpbg_df __ckpbg_minp __ckpbg_mino __ckpbg_warn

    * stats vector order as in the GAUSS retp: 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 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' (`bp')"
        }
    }

    * ---- display ----
    if "`print'" == "" {
        local mlab "constant, no breaks"
        if `mnum' == 1 local mlab "linear trend, no breaks"
        if `mnum' == 2 local mlab "breaks in the trend slope (Model I)"
        if `mnum' == 3 local mlab "breaks in level and slope (Model II)"
        local blab "estimated, brute-force GLS-SSR grid search"
        if `est' == 1 local blab "estimated, dynamic programming (dp)"
        if `known' local blab "known (user supplied)"
        local plab "MAIC"
        if `pen' == 1 local plab "BIC"
        if `bgsel' local plab "BG"
        di
        di as text "{hline 78}"
        di as text "GLS-based unit root tests with multiple structural breaks"
        di as text "Carrion-i-Silvestre, Kim, and Perron (2009, Econometric Theory 25, 1754-1792)"
        di as text "{hline 78}"
        di as text "Variable   : " as result "`dispname'" 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'
        if `bgminp' < . {
            di as text "BG min p   : " as result %6.4f `bgminp' as text "   (spectral AR residuals, orders 1 to `bgh', df floor `bgdf')"
        }
        else {
            di as text "BG min p   : no order with df >= `bgdf' available (orders 1 to `bgh')"
        }
        di as text "{hline 78}"
        di as text "   Test        Statistic        1%      2.5%        5%       10%"
        di as text "{hline 78}"
        * cv row mapping follows the author's sample.gss: PT,MPT -> critical_pt ; ADF,MZT -> critical_mzt ; ZA,MZA -> critical_mza ; MSB -> critical_msb
        * rows of CV: 1=msb 2=mza 3=mzt 4=pt ; columns: 1% 2.5% 5% 10%
        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 c25 = `CV'[`cr', 2]
            local c5 = `CV'[`cr', 3]
            local c10 = `CV'[`cr', 4]
            local st ""
            if `s' < `c1' local st "***"
            else if `s' < `c5' local st "**"
            else if `s' < `c10' local st "*"
            di as text %7s "`nm'" as result %14.3f `s' as text %-4s "`st'" as result %9.3f `c1' %10.3f `c25' %10.3f `c5' %10.3f `c10'
        }
        di as text "{hline 78}"
        di as text "H0: unit root.  * p<0.10, ** p<0.05, *** p<0.01."
        if `bgwarn' {
            di as text "Warning: autocorrelation could not be eliminated with kmax = `kmax' lags;"
            di as text "kmax was used. Consider increasing kmax()."
        }
        if !`bgsel' & `bgminp' < 0.05 {
            di as text "Note: the minimum BG p-value is below 0.05: the spectral AR residuals show"
            di as text "remaining autocorrelation. Consider the {bf:bg} option or a larger kmax()."
        }
        di as text "{hline 78}"
    }

    * ---- saved results ----
    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 bgminp = `bgminp'
    return scalar bglags = `bgh'
    return scalar bgdf = `bgdf'
    return scalar nbreaks = cond(`mnum' >= 2, `nbfound', 0)
    return local model "`mnum'"
    return local method = cond(`est' == 1, "dp", "brute")
    return local penalty "`penalty'"
    return local breakdates "`bdates'"
    return local varname "`dispname'"
    if "`pretest'" == "" {
        return scalar pyexpw = `pyexpw'
        return scalar pybreakpos = `pybp'
        return scalar pybreakdate = `pybd'
        return scalar pycv10 = `pycv10'
        return scalar pycv5 = `pycv5'
        return scalar pycv1 = `pycv1'
    }
    return local cmd "ckptest"
    matrix rownames `CV' = MSB MZA MZT PT
    matrix colnames `CV' = cv1 cv2_5 cv5 cv10
    return matrix cv = `CV', copy
    if `mnum' >= 2 {
        return matrix breakpos = `TBOUT', copy
    }
end

mata:
// ===========================================================================
// ckptest Mata library: line-by-line port of msbur.src (CKP 2009 GAUSS code)
// Every function below carries the name of the GAUSS procedure it replicates.
// ===========================================================================

// ---------------------------------------------------------------------------
// GAUSS primitives
// ---------------------------------------------------------------------------

// GAUSS trimr(x,a,b): drop a rows from the top and b rows from the bottom
real matrix _ckp_trimr(real matrix x, real scalar a, real scalar b)
{
    return(x[(a+1)::(rows(x)-b), .])
}

// GAUSS seqa(a,b,n): additive sequence a, a+b, ..., a+(n-1)b (column)
real colvector _ckp_seqa(real scalar a, real scalar b, real scalar n)
{
    return(a :+ b:*(0::(n-1)))
}

// GAUSS lagn(x,n)
real matrix _ckp_lagn(real matrix x, real scalar n)
{
    if (n > 0) {
        return(J(n, cols(x), 0) \ _ckp_trimr(x, 0, n))
    }
    else {
        return(_ckp_trimr(x, abs(n), 0) \ J(abs(n), cols(x), 0))
    }
}

// GAUSS diff(x,k)
real matrix _ckp_diff(real matrix x, real scalar k)
{
    if (k == 0) {
        return(x)
    }
    return(J(k, cols(x), 0) \ (_ckp_trimr(x, k, 0) - _ckp_trimr(_ckp_lagn(x, k), k, 0)))
}

// GAUSS minindc: index of the first occurrence of the column minimum
real scalar _ckp_minindc(real colvector v)
{
    real scalar i, n, w
    w = min(v)
    n = rows(v)
    for (i = 1; i <= n; i++) {
        if (v[i] == w) {
            return(i)
        }
    }
    return(.)
}

// GAUSS olsqr2(y,X): QR least squares; returns coefficients and residuals
void _ckp_olsqr2(real colvector y, real matrix X, real colvector b, real colvector e)
{
    b = qrsolve(X, y)
    e = y - X*b
}

// ---------------------------------------------------------------------------
// GAUSS proc(2)=glsd(y,z,cbar): quasi-GLS detrending
// ---------------------------------------------------------------------------
void _ckp_glsd(real colvector y, real matrix z, real scalar cbar, real colvector yt, real scalar ssr)
{
    real scalar nt, abar
    real colvector ya, bhat, eg
    real matrix za

    nt = rows(y)
    abar = 1 + cbar/nt
    ya = J(nt, 1, 0)
    za = J(nt, cols(z), 0)
    ya[1] = y[1]
    za[1, .] = z[1, .]
    ya[2::nt] = y[2::nt] - abar*y[1::(nt-1)]
    za[2::nt, .] = z[2::nt, .] - abar*z[1::(nt-1), .]
    // GAUSS: bhat=inv(za'za)*za'ya
    bhat = luinv(cross(za, za))*cross(za, ya)
    yt = y - z*bhat
    eg = ya - za*bhat
    ssr = cross(eg, eg)
}

// ---------------------------------------------------------------------------
// GAUSS proc(2)=olsd(y,z): OLS detrending, only the residual vector is used
// downstream (krule=s2ar(yt_ols,...))
// ---------------------------------------------------------------------------
real colvector _ckp_olsd(real colvector y, real matrix z)
{
    real colvector ahat, yd
    _ckp_olsqr2(y, z, ahat, yd)
    return(yd)
}

// ---------------------------------------------------------------------------
// GAUSS proc(1)=s2ar(yts,penalty,kmax,kmin): lag choice for the long-run
// variance (MAIC penalty=0, BIC penalty=1)
// ---------------------------------------------------------------------------
real scalar _ckp_s2ar(real colvector yts, real scalar penalty, real scalar kmax, real scalar kmin)
{
    real scalar nt, nef, k, i, kopt, sumy
    real colvector dyts, dyts0, b, e, s2e, tau, mic, kk
    real matrix reg, reg0

    nt = rows(yts)
    tau = J(kmax+1, 1, 0)
    s2e = 999*J(kmax+1, 1, 1)

    dyts = _ckp_diff(yts, 1)
    reg = _ckp_lagn(yts, 1)

    i = 1
    while (i <= kmax) {
        reg = reg, _ckp_lagn(dyts, i)
        i = i + 1
    }

    dyts0 = _ckp_trimr(dyts, kmax+1, 0)
    reg0 = _ckp_trimr(reg, kmax+1, 0)
    sumy = sum(reg0[., 1]:*reg0[., 1])
    nef = nt - kmax - 1

    k = kmin
    while (k <= kmax) {
        // GAUSS: b=dyts0/reg0[.,1:k+1] (matrix division = least squares)
        b = qrsolve(reg0[., 1::(k+1)], dyts0)
        e = dyts0 - reg0[., 1::(k+1)]*b
        s2e[k+1] = cross(e, e)/nef
        tau[k+1] = (b[1]*b[1])*sumy/s2e[k+1]
        k = k + 1
    }

    kk = _ckp_seqa(0, 1, kmax+1)

    if (penalty == 0) {
        mic = ln(s2e) + 2.0:*(kk + tau):/nef
    }
    else {
        mic = ln(s2e) + ln(nef):*kk:/nef
    }

    kopt = _ckp_minindc(mic) - 1
    return(kopt)
}

// ---------------------------------------------------------------------------
// GAUSS proc(3)=adfp(yt,kstar): ADF regression on the detrended series
// returns (adf, rho[1]+1, s2vec)
// ---------------------------------------------------------------------------
void _ckp_adfp(real colvector yt, real scalar kstar, real scalar adf, real scalar a1, real scalar s2vec)
{
    real scalar i, nef, s2e, sre, sumb
    real colvector dyt, rho, ee
    real matrix reg, xx

    reg = _ckp_lagn(yt, 1)
    dyt = _ckp_diff(yt, 1)

    i = 1
    while (i <= kstar) {
        reg = reg, _ckp_lagn(dyt, i)
        i = i + 1
    }

    dyt = _ckp_trimr(dyt, kstar+1, 0)
    reg = _ckp_trimr(reg, kstar+1, 0)

    _ckp_olsqr2(dyt, reg, rho, ee)

    nef = rows(dyt)
    s2e = cross(ee, ee)/nef
    // GAUSS: xx=inv(reg'reg)
    xx = luinv(cross(reg, reg))
    sre = xx[1, 1]*s2e
    adf = rho[1]/sqrt(sre)

    if (kstar > 0) {
        sumb = sum(rho[2::(kstar+1)])
    }
    else {
        sumb = 0
    }

    s2vec = s2e/(1 - sumb)^2
    a1 = rho[1] + 1
}

// ---------------------------------------------------------------------------
// GAUSS proc(1)=c_bar_rs(lam): response surface for the c_bar parameter,
// (5x1) vector of break fractions; regressor order exactly as in the source:
// 1, lam1..lam5, lam^2, lam^3, lam^4, |lami-lamj| (pairs (1,2)(1,3)(1,4)(1,5)
// (2,3)(2,4)(2,5)(3,4)(3,5)(4,5)), then the same pairs squared, cubed, ^4
// ---------------------------------------------------------------------------
real scalar _ckp_cbar_rs(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)
}

// ---------------------------------------------------------------------------
// GAUSS proc(4)=msbur_rs(lam,c_bar): response surfaces for the critical
// values; returns a 4x4 matrix, rows MSB / MZA / MZT / PT, columns
// 1% / 2.5% / 5% / 10% (exactly the four levels of the source).
// Regressor order exactly as in the source: 1, lam, lam^2, c_bar, lam*c_bar,
// c_bar^2, lam*c_bar^2, |di-dj|*c_bar, |di-dj|^2*c_bar, ^3*c_bar, ^4*c_bar
// ---------------------------------------------------------------------------
real matrix _ckp_msbur_rs(real colvector lam, real scalar cbar)
{
    real rowvector x, critical
    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 = _ckp_cvparam()
    critical = x*P
    // GAUSS: critical_msb=critical[1:4]; _mza=[5:8]; _mzt=[9:12]; _pt=[13:16]
    return((critical[1], critical[2], critical[3], critical[4] \ critical[5], critical[6], critical[7], critical[8] \ critical[9], critical[10], critical[11], critical[12] \ critical[13], critical[14], critical[15], critical[16]))
}

real matrix _ckp_cvparam()
{
    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)
}

// ---------------------------------------------------------------------------
// Deterministic regressor matrix for models 2 and 3, column order exactly as
// in the GAUSS source: cns ~ tend ~ [du1 ~] dt1 ~ [du2 ~] dt2 ~ ...
// (model 3 interleaves du_j and dt_j per break; model 2 has only dt_j)
// ---------------------------------------------------------------------------
real matrix _ckp_zbuild(real scalar nt, real scalar model1, real colvector tb)
{
    real matrix z
    real colvector du, dt
    real scalar j, m

    z = J(nt, 1, 1), _ckp_seqa(1, 1, nt)
    m = rows(tb)
    for (j = 1; j <= m; j++) {
        dt = J(tb[j], 1, 0) \ _ckp_seqa(1, 1, nt-tb[j])
        if (model1 == 3) {
            du = J(tb[j], 1, 0) \ J(nt-tb[j], 1, 1)
            z = z, du, dt
        }
        else {
            z = z, dt
        }
    }
    return(z)
}

// ---------------------------------------------------------------------------
// Common statistics tail, identical in sbur_multiple_gls_brute (lines
// following the break search) and sbur_multiple_gls_algorithm; transcribed
// once and used by both paths.
// ---------------------------------------------------------------------------
void _ckp_tail(real colvector y, real matrix z, real scalar cbar, real scalar model1, real scalar penalty, real scalar kmax, real scalar kmin, real rowvector stats, real scalar krule)
{
    real scalar nt, ssra, s2u, sumyt2, adf, a1, sar, bt, za, mza, msb, mzt, ssr1, pt, mpt, bgwarn, bgminp, bgmino
    real colvector yt, r, yt_ols, ytf
    real scalar ahat_s
    real colvector ahat

    nt = rows(y)

    _ckp_glsd(y, z, cbar, yt, ssra)
    // GAUSS: {ahat,r,fit}=olsqr2(yt[2:nt,1],yt[1:nt-1,1])
    _ckp_olsqr2(yt[2::nt], yt[1::(nt-1)], ahat, r)
    ahat_s = ahat[1]
    s2u = cross(r, r)/(nt-1)
    sumyt2 = sum(yt[1::(nt-1)]:^2)/(nt-1)^2

    yt_ols = _ckp_olsd(y, z)
    if (st_numscalar("__ckpbg_sel") == 1) {
        bgwarn = 0
        krule = _ckp_bgsel(yt, kmax, kmin, st_numscalar("__ckpbg_h"), st_numscalar("__ckpbg_df"), bgwarn)
        st_numscalar("__ckpbg_warn", bgwarn)
    }
    else {
        krule = _ckp_s2ar(yt_ols, penalty, kmax, kmin)
        st_numscalar("__ckpbg_warn", 0)
    }
    _ckp_bgmin(yt, krule, st_numscalar("__ckpbg_h"), st_numscalar("__ckpbg_df"), bgminp, bgmino)
    st_numscalar("__ckpbg_minp", bgminp)
    st_numscalar("__ckpbg_mino", bgmino)
    _ckp_adfp(yt, krule, adf, a1, sar)

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

    _ckp_glsd(y, z, 0, ytf, ssr1)
    pt = (ssra - (1 + cbar/nt)*ssr1)/sar

    if (model1 == 0) {
        mpt = (cbar*cbar*sumyt2 - cbar*(yt[nt]^2)/nt)/sar
    }
    else {
        mpt = (cbar*cbar*sumyt2 + (1 - cbar)*(yt[nt]^2)/nt)/sar
    }

    stats = (pt, mpt, adf, za, mza, msb, mzt)
}

// ---------------------------------------------------------------------------
// GAUSS proc(9)=sbur_multiple_gls_brute(y,model,penalty,kmax,kmin)
// ---------------------------------------------------------------------------
void _ckp_brute(real colvector y, real scalar model1, real scalar known, real colvector tbin, real scalar nbrk, real scalar penalty, real scalar kmax, real scalar kmin, real rowvector stats, real colvector min_tb, real scalar cbar, real scalar krule)
{
    real scalar nt, j, jj, jjj, min_ssra, ssra
    real colvector yt, tbc
    real matrix z

    nt = rows(y)
    min_tb = J(1, 1, 0)

    if (model1 == 0) {
        z = J(nt, 1, 1)
        cbar = -7.0
    }
    else if (model1 == 1) {
        z = J(nt, 1, 1), _ckp_seqa(1, 1, nt)
        cbar = -13.5
    }
    else {
        if (known == 1) {
            // GAUSS known-break branch: min_tb = model[3:rows(model)]
            min_tb = tbin
            z = _ckp_zbuild(nt, model1, min_tb)
            if (rows(min_tb) < 5) {
                cbar = _ckp_cbar_rs((min_tb:/nt) \ J(5-rows(min_tb), 1, 0))
            }
            else {
                cbar = _ckp_cbar_rs(min_tb:/nt)
            }
        }
        else {
            // GAUSS unknown-break branch. SINGLE DOCUMENTED DEVIATION from the
            // source: GAUSS initializes min_ssra = y'y, which is not an upper
            // bound of the quasi-differenced SSR for strongly stationary series;
            // in that case no candidate is ever accepted, min_tb stays 0, and the
            // original code crashes on a singular regressor matrix (dt duplicates
            // tend). Initializing with maxdouble() returns the identical argmin
            // whenever the GAUSS code works (the minimizing candidate is then
            // below y'y anyway) and a valid result where the GAUSS code fails.
            min_ssra = maxdouble()
            if (nbrk == 1) {
                j = 3
                while (j <= nt-3) {
                    tbc = (j)
                    z = _ckp_zbuild(nt, model1, tbc)
                    cbar = _ckp_cbar_rs((j/nt) \ 0 \ 0 \ 0 \ 0)
                    _ckp_glsd(y, z, cbar, yt, ssra)
                    if (ssra < min_ssra) {
                        min_tb = tbc
                        min_ssra = ssra
                    }
                    j = j + 1
                }
                z = _ckp_zbuild(nt, model1, min_tb)
                cbar = _ckp_cbar_rs((min_tb[1]/nt) \ 0 \ 0 \ 0 \ 0)
            }
            else if (nbrk == 2) {
                j = 3
                while (j <= nt-3-2) {
                    jj = j + 2
                    while (jj <= nt-3) {
                        tbc = (j \ jj)
                        z = _ckp_zbuild(nt, model1, tbc)
                        cbar = _ckp_cbar_rs((j/nt) \ (jj/nt) \ 0 \ 0 \ 0)
                        _ckp_glsd(y, z, cbar, yt, ssra)
                        if (ssra < min_ssra) {
                            min_tb = tbc
                            min_ssra = ssra
                        }
                        jj = jj + 1
                    }
                    j = j + 1
                }
                z = _ckp_zbuild(nt, model1, min_tb)
                cbar = _ckp_cbar_rs((min_tb[1]/nt) \ (min_tb[2]/nt) \ 0 \ 0 \ 0)
            }
            else {
                j = 3
                while (j <= nt-3-4) {
                    jj = j + 2
                    while (jj <= nt-3-2) {
                        jjj = jj + 2
                        while (jjj <= nt-3) {
                            tbc = (j \ jj \ jjj)
                            z = _ckp_zbuild(nt, model1, tbc)
                            cbar = _ckp_cbar_rs((j/nt) \ (jj/nt) \ (jjj/nt) \ 0 \ 0)
                            _ckp_glsd(y, z, cbar, yt, ssra)
                            if (ssra < min_ssra) {
                                min_tb = tbc
                                min_ssra = ssra
                            }
                            jjj = jjj + 1
                        }
                        jj = jj + 1
                    }
                    j = j + 1
                }
                z = _ckp_zbuild(nt, model1, min_tb)
                cbar = _ckp_cbar_rs((min_tb[1]/nt) \ (min_tb[2]/nt) \ (min_tb[3]/nt) \ 0 \ 0)
            }
        }
    }

    _ckp_tail(y, z, cbar, model1, penalty, kmax, kmin, stats, krule)
}

// ---------------------------------------------------------------------------
// GAUSS proc(1)=ssr(start,y,z,h,last): recursive SSR, OLS detrending stage
// ---------------------------------------------------------------------------
real colvector _ckp_ssr(real scalar start, real colvector y, real matrix z, real scalar h, real scalar last)
{
    real colvector vecssr, delta1, delta2, invz, res
    real matrix inv1, inv2
    real scalar v, f, r

    vecssr = J(last, 1, 0)

    inv1 = luinv(cross(z[start::(start+h-1), .], z[start::(start+h-1), .]))
    delta1 = inv1*cross(z[start::(start+h-1), .], y[start::(start+h-1)])
    res = y[start::(start+h-1)] - z[start::(start+h-1), .]*delta1
    vecssr[start+h-1] = cross(res, res)

    r = start + h
    while (r <= last) {
        v = y[r] - z[r, .]*delta1
        invz = inv1*z[r, .]'
        f = 1 + z[r, .]*invz
        delta2 = delta1 + (invz*v)/f
        inv2 = inv1 - (invz*invz')/f
        inv1 = inv2
        delta1 = delta2
        vecssr[r] = vecssr[r-1] + v*v/f
        r = r + 1
    }

    return(vecssr)
}

// ---------------------------------------------------------------------------
// GAUSS proc(1)=ssr_gls(start,y,z,h,last): recursive SSR, GLS stage, with an
// initial impulse dummy that equals 1 at the first observation of the segment
// ---------------------------------------------------------------------------
real colvector _ckp_ssr_gls(real scalar start, real colvector y, real matrix z, real scalar h, real scalar last)
{
    real colvector vecssr, delta1, delta2, invz, res, initial_impulse
    real matrix inv1, inv2, zb
    real rowvector zr
    real scalar v, f, r

    vecssr = J(last, 1, 0)

    initial_impulse = 1 \ J(rows(z)-1, 1, 0)

    zb = z[start::(start+h-1), .], initial_impulse[1::h]
    inv1 = luinv(cross(zb, zb))
    delta1 = inv1*cross(zb, y[start::(start+h-1)])
    res = y[start::(start+h-1)] - zb*delta1
    vecssr[start+h-1] = cross(res, res)

    r = start + h
    while (r <= last) {
        zr = z[r, .], initial_impulse[r-start+1]
        v = y[r] - zr*delta1
        invz = inv1*zr'
        f = 1 + zr*invz
        delta2 = delta1 + (invz*v)/f
        inv2 = inv1 - (invz*invz')/f
        inv1 = inv2
        delta1 = delta2
        vecssr[r] = vecssr[r-1] + v*v/f
        r = r + 1
    }

    return(vecssr)
}

// ---------------------------------------------------------------------------
// GAUSS proc(2)=parti(start,b1,b2,last,bigvec,bigt)
// ---------------------------------------------------------------------------
void _ckp_parti(real scalar start, real scalar b1, real scalar b2, real scalar last, real colvector bigvec, real scalar bigt, real scalar ssrmin, real scalar dx)
{
    real colvector dvec
    real scalar j, jj, k, ini

    dvec = J(bigt, 1, 0)
    ini = (start-1)*bigt - (start-2)*(start-1)/2 + 1

    j = b1
    while (j <= b2) {
        jj = j - start
        k = j*bigt - (j-1)*j/2 + last - j
        dvec[j] = bigvec[ini+jj] + bigvec[k]
        j = j + 1
    }

    ssrmin = min(dvec[b1::b2])
    dx = (b1-1) + _ckp_minindc(dvec[b1::b2])
}

// ---------------------------------------------------------------------------
// GAUSS proc(3)=dating(y,z,h,m,q,bigt): Bai-Perron dynamic programming, OLS
// stage (the q argument is kept to mirror the GAUSS signature; it is unused
// in the GAUSS body as well)
// ---------------------------------------------------------------------------
void _ckp_dating(real colvector y, real matrix z, real scalar h, real scalar m, real scalar q, real scalar bigt, real colvector glob, real matrix datevec, real colvector bigvec, real scalar glsflag)
{
    real matrix optdat, optssr
    real colvector dvec, vecssr
    real scalar i, ssrmin, datx, j1, ib, jlast, jb, xx

    datevec = J(m, m, 0)
    optdat = J(bigt, m, 0)
    optssr = J(bigt, m, 0)
    dvec = J(bigt, 1, 0)
    glob = J(m, 1, 0)
    bigvec = J(bigt*(bigt+1)/2, 1, 0)

    i = 1
    while (i <= bigt-h+1) {
        if (glsflag == 0) {
            vecssr = _ckp_ssr(i, y, z, h, bigt)
        }
        else {
            vecssr = _ckp_ssr_gls(i, y, z, h, bigt)
        }
        bigvec[((i-1)*bigt+i-(i-1)*i/2)::(i*bigt-(i-1)*i/2)] = vecssr[i::bigt]
        i = i + 1
    }

    if (m == 1) {
        _ckp_parti(1, h, bigt-h, bigt, bigvec, bigt, ssrmin, datx)
        datevec[1, 1] = datx
        glob[1] = ssrmin
    }
    else {
        j1 = 2*h
        while (j1 <= bigt) {
            _ckp_parti(1, h, j1-h, j1, bigvec, bigt, ssrmin, datx)
            optssr[j1, 1] = ssrmin
            optdat[j1, 1] = datx
            j1 = j1 + 1
        }
        glob[1] = optssr[bigt, 1]
        datevec[1, 1] = optdat[bigt, 1]

        ib = 2
        while (ib <= m) {
            if (ib == m) {
                jlast = bigt
                jb = ib*h
                while (jb <= jlast-h) {
                    dvec[jb] = optssr[jb, ib-1] + bigvec[(jb+1)*bigt - jb*(jb+1)/2]
                    jb = jb + 1
                }
                optssr[jlast, ib] = min(dvec[(ib*h)::(jlast-h)])
                optdat[jlast, ib] = (ib*h-1) + _ckp_minindc(dvec[(ib*h)::(jlast-h)])
            }
            else {
                jlast = (ib+1)*h
                while (jlast <= bigt) {
                    jb = ib*h
                    while (jb <= jlast-h) {
                        dvec[jb] = optssr[jb, ib-1] + bigvec[jb*bigt - jb*(jb-1)/2 + jlast - jb]
                        jb = jb + 1
                    }
                    optssr[jlast, ib] = min(dvec[(ib*h)::(jlast-h)])
                    optdat[jlast, ib] = (ib*h-1) + _ckp_minindc(dvec[(ib*h)::(jlast-h)])
                    jlast = jlast + 1
                }
            }

            datevec[ib, ib] = optdat[bigt, ib]

            i = 1
            while (i <= ib-1) {
                xx = ib - i
                datevec[xx, ib] = optdat[datevec[xx+1, ib], xx]
                i = i + 1
            }
            glob[ib] = optssr[bigt, ib]

            ib = ib + 1
        }
    }
}

// ---------------------------------------------------------------------------
// GAUSS proc(1)=ssr2_gls(y,z,b,q,br): SSR contributions for all segments
// under estimated (restricted) coefficients b
// ---------------------------------------------------------------------------
real colvector _ckp_ssr2_gls(real colvector y, real matrix z, real colvector b, real scalar q, real colvector br)
{
    real scalar m, i, bigt
    real colvector bigvec2
    real matrix initial_impulse

    m = rows(br)
    bigt = rows(y)
    bigvec2 = J(bigt*(m+1), 1, 0)

    initial_impulse = 1 \ J(bigt-1, 1, 0)
    i = 1
    while (i <= m) {
        initial_impulse = initial_impulse, (J(br[i], 1, 0) \ 1 \ J(bigt-br[i]-1, 1, 0))
        i = i + 1
    }

    i = 1
    while (i <= m+1) {
        bigvec2[((i-1)*bigt+1)::(i*bigt)] = (y - (z, initial_impulse[., i])*b[((i-1)*(q+1)+1)::(i*(q+1))]):^2
        i = i + 1
    }

    return(bigvec2)
}

// ---------------------------------------------------------------------------
// GAUSS proc(2)=parti2(start,b1,b2,last,bigvec2,bigt)
// ---------------------------------------------------------------------------
void _ckp_parti2(real scalar start, real scalar b1, real scalar b2, real scalar last, real colvector bigvec2, real scalar bigt, real scalar ssrmin, real scalar dx)
{
    real colvector dvec
    real scalar j, ssr1, ssr2

    dvec = J(bigt, 1, 0)
    j = b1
    while (j <= b2) {
        ssr1 = sum(bigvec2[start::j])
        ssr2 = sum(bigvec2[(1*bigt+j+1)::(1*bigt+last)])
        dvec[j] = ssr1 + ssr2
        j = j + 1
    }

    ssrmin = min(dvec[b1::b2])
    dx = (b1-1) + _ckp_minindc(dvec[b1::b2])
}

// ---------------------------------------------------------------------------
// GAUSS proc(2)=dating2_gls(bigvec2,h,m,bigt): dynamic programming under the
// restricted coefficient estimates
// ---------------------------------------------------------------------------
void _ckp_dating2_gls(real colvector bigvec2, real scalar h, real scalar m, real scalar bigt, real colvector glob, real matrix datevec)
{
    real matrix optdat, optssr
    real colvector dvec
    real scalar i, ssrmin, datx, j1, ib, jlast, jb, xx

    datevec = J(m, m, 0)
    optdat = J(bigt, m, 0)
    optssr = J(bigt, m, 0)
    dvec = J(bigt, 1, 0)
    glob = J(m, 1, 0)

    if (m == 1) {
        _ckp_parti2(1, h, bigt-h, bigt, bigvec2, bigt, ssrmin, datx)
        datevec[1, 1] = datx
        glob[1] = ssrmin
    }
    else {
        j1 = 2*h
        while (j1 <= bigt) {
            _ckp_parti2(1, h, j1-h, j1, bigvec2[1::(2*bigt)], bigt, ssrmin, datx)
            optssr[j1, 1] = ssrmin
            optdat[j1, 1] = datx
            j1 = j1 + 1
        }

        glob[1] = optssr[bigt, 1]
        datevec[1, 1] = optdat[bigt, 1]

        ib = 2
        while (ib <= m) {
            if (ib == m) {
                jlast = bigt
                jb = ib*h
                while (jb <= jlast-h) {
                    dvec[jb] = optssr[jb, ib-1] + sum(bigvec2[(bigt*m+jb+1)::(bigt*(m+1))])
                    jb = jb + 1
                }
                optssr[jlast, ib] = min(dvec[(ib*h)::(jlast-h)])
                optdat[jlast, ib] = (ib*h-1) + _ckp_minindc(dvec[(ib*h)::(jlast-h)])
            }
            else {
                jlast = (ib+1)*h
                while (jlast <= bigt) {
                    jb = ib*h
                    while (jb <= jlast-h) {
                        dvec[jb] = optssr[jb, ib-1] + sum(bigvec2[(bigt*ib+jb+1)::(bigt*ib+jlast)])
                        jb = jb + 1
                    }
                    optssr[jlast, ib] = min(dvec[(ib*h)::(jlast-h)])
                    optdat[jlast, ib] = (ib*h-1) + _ckp_minindc(dvec[(ib*h)::(jlast-h)])
                    jlast = jlast + 1
                }
            }

            datevec[ib, ib] = optdat[bigt, ib]

            i = 1
            while (i <= ib-1) {
                xx = ib - i
                datevec[xx, ib] = optdat[datevec[xx+1, ib], xx]
                i = i + 1
            }
            glob[ib] = optssr[bigt, ib]

            ib = ib + 1
        }
    }
}

// ---------------------------------------------------------------------------
// Restriction matrix R for est2_gls, models 2 and 3, exactly as in the source
// ---------------------------------------------------------------------------
real matrix _ckp_est2_R(real scalar model1, real scalar alpha, real colvector brdates, real scalar m)
{
    real matrix R, temp
    real scalar i

    R = (-alpha, 0, 1), J(1, m*3, 0)
    temp = J(m, (m+1)*3, 0)
    i = 1
    while (i <= m) {
        if (model1 == 2) {
            temp[i, (3*i-2)::((i+1)*3)] = (-1, -brdates[i], 0, 1, brdates[i], 0)
        }
        else {
            temp[i, (3*i-2)::((i+1)*3)] = (-1, -brdates[i], 0, 1, brdates[i], -1/alpha)
        }
        i = i + 1
    }
    R = R \ temp
    if (model1 == 2) {
        temp = J(m, (m+1)*3, 0)
        i = 1
        while (i <= m) {
            temp[i, 3*(i+1)] = 1
            i = i + 1
        }
        R = R \ temp
    }
    return(R)
}

// ---------------------------------------------------------------------------
// Block diagonal regressor matrix zbar for est2_gls, exactly as in the source
// ---------------------------------------------------------------------------
real matrix _ckp_est2_zbar(real matrix z, real colvector initial_impulse, real colvector dvm, real scalar q, real scalar m, real scalar bigt)
{
    real matrix zbar, zi
    real scalar i, t1, t2, tbar

    tbar = dvm[1]
    zbar = _ckp_trimr((z, initial_impulse), 0, bigt-tbar) \ J(bigt-tbar, q+1, 0)

    i = 2
    while (i <= m) {
        t1 = dvm[i-1]
        t2 = dvm[i]
        zi = J(t1, q+1, 0) \ (_ckp_trimr(z, t1, bigt-t2), initial_impulse[1::(t2-t1)]) \ J(bigt-t2, q+1, 0)
        zbar = zbar, zi
        i = i + 1
    }

    t2 = dvm[m]
    zi = J(t2, q+1, 0) \ (_ckp_trimr(z, t2, 0), initial_impulse[1::(bigt-t2)])
    zbar = zbar, zi

    return(zbar)
}

// ---------------------------------------------------------------------------
// GAUSS proc(3)=est2_gls(y,model,q,m,bigt,trm,datevec): restricted estimation
// of breaks and coefficients (Perron and Qu 2004 restricted SSR code)
// ---------------------------------------------------------------------------
void _ckp_est2_gls(real colvector y, real scalar model1, real scalar q, real scalar m, real scalar bigt, real scalar trm, real colvector datevec_in, real colvector dx, real colvector delta, real scalar ssr_out)
{
    real scalar h, tstar, c_bar, alpha, count, ssr_iter_prev, rr, repeatflag
    real colvector br, y_gls, cns_gls, tend_gls, initial_impulse, b, e, glob2, bigvec2
    real matrix z, zbar, zz, R, datevec2

    tstar = bigt
    h = round(trm*bigt)

    br = datevec_in

    if (m < 5) {
        c_bar = _ckp_cbar_rs(sort(br:/bigt, 1) \ J(5-m, 1, 0))
    }
    else {
        c_bar = _ckp_cbar_rs(sort(br:/bigt, 1))
    }

    alpha = 1 + c_bar/bigt

    cns_gls = (1-alpha)*J(bigt, 1, 1)
    tend_gls = 1 \ ((-c_bar/bigt)*_ckp_seqa(1, 1, bigt-1) :+ 1)
    z = cns_gls, tend_gls
    y_gls = y[1] \ (y[2::bigt] - alpha*y[1::(bigt-1)])

    initial_impulse = 1 \ J(rows(z)-1, 1, 0)

    zbar = _ckp_est2_zbar(z, initial_impulse, br, q, m, bigt)

    _ckp_olsqr2(y_gls, zbar, b, e)
    // GAUSS: zz=invpd(zbar'zbar)
    zz = cholinv(cross(zbar, zbar))

    R = _ckp_est2_R(model1, alpha, br, m)

    rr = 0

    // GAUSS: delta=b+zz*R'*invpd(R*zz*R')*(rr-R*b)
    delta = b + zz*R'*cholinv(R*zz*R')*(rr :- R*b)

    // iterative estimation procedure (GAUSS goto start_iter loop, max 10)
    count = 0
    ssr_iter_prev = cross(y_gls - zbar*delta, y_gls - zbar*delta)

    repeatflag = 1
    while (repeatflag) {
        bigvec2 = _ckp_ssr2_gls(y_gls, z, delta, q, br)
        _ckp_dating2_gls(bigvec2, h, m, tstar, glob2, datevec2)

        br = datevec2[., m]

        if (m < 5) {
            c_bar = _ckp_cbar_rs(sort(br:/bigt, 1) \ J(5-m, 1, 0))
        }
        else {
            c_bar = _ckp_cbar_rs(sort(br:/bigt, 1))
        }

        alpha = 1 + c_bar/bigt

        cns_gls = (1-alpha)*J(bigt, 1, 1)
        tend_gls = 1 \ ((-c_bar/bigt)*_ckp_seqa(1, 1, bigt-1) :+ 1)
        z = cns_gls, tend_gls
        y_gls = y[1] \ (y[2::bigt] - alpha*y[1::(bigt-1)])

        zbar = _ckp_est2_zbar(z, initial_impulse, datevec2[., m], q, m, tstar)
        _ckp_olsqr2(y_gls, zbar, b, e)
        zz = cholinv(cross(zbar, zbar))

        R = _ckp_est2_R(model1, alpha, br, m)

        delta = b + zz*R'*cholinv(R*zz*R')*(rr :- R*b)
        e = y_gls - zbar*delta

        // GAUSS: if (count < 10) and (abs(e'e-ssr_iter_prev) > 1e-3)
        if (count < 10 & abs(cross(e, e) - ssr_iter_prev) > 1e-3) {
            count = count + 1
            ssr_iter_prev = cross(e, e)
        }
        else {
            repeatflag = 0
        }
    }

    dx = br
    ssr_out = cross(e, e)
}

// ---------------------------------------------------------------------------
// GAUSS proc(9)=sbur_multiple_gls_algorithm(y,model,penalty,kmax,kmin,maxiter)
// ---------------------------------------------------------------------------
void _ckp_algorithm(real colvector y, real scalar model1, real scalar nbrk, real scalar penalty, real scalar kmax, real scalar kmin, real scalar maxiter, real rowvector stats, real colvector min_tb, real scalar cbar, real scalar krule)
{
    real scalar t, trm, bigt, h, m, q, count, differe, alpha, ssr_prev_iter, ssr_est_prev, ssr_est, contflag
    real colvector datevec, bigvec, y_gls, cns_gls, tend_gls, dx, delta, glob
    real matrix z, z_gls, datevec_mat, datevec_iteration

    real scalar verbose, ndots

    t = rows(y)
    ndots = 0

    // progress feedback: dp iterations can take a while for several breaks
    verbose = (st_local("print") == "" & nbrk >= 2)
    if (verbose) {
        printf("{txt}(dynamic programming for %g breaks ", nbrk)
        displayflush()
    }

    // Computing the break points using OLS
    z = J(t, 1, 1), _ckp_seqa(1, 1, t)

    trm = 0.10
    bigt = t
    h = trunc(trm*bigt)
    m = nbrk
    q = cols(z)

    _ckp_dating(y, z, h, m, q, bigt, glob, datevec_mat, bigvec, 0)

    // Computing the break points using GLS
    datevec = select(datevec_mat[., m], datevec_mat[., m] :> 0)
    datevec = sort(datevec, 1)

    count = 0
    ssr_prev_iter = glob[m]

    differe = 0
    while (differe != -1) {
        if (m < 5) {
            cbar = _ckp_cbar_rs((datevec:/t) \ J(5-m, 1, 0))
        }
        else {
            cbar = _ckp_cbar_rs(datevec:/t)
        }

        alpha = 1 + cbar/t

        cns_gls = 1 \ ((1-alpha)*J(t-1, 1, 1))
        tend_gls = 1 \ ((-cbar/t)*_ckp_seqa(2, 1, t-1) :+ 1)
        z_gls = cns_gls, tend_gls
        y_gls = y[1] \ (y[2::t] - alpha*y[1::(t-1)])

        _ckp_dating(y_gls, z_gls, h, m, q, bigt, glob, datevec_iteration, bigvec, 1)

        if (count < maxiter & abs(glob[m] - ssr_prev_iter) > 1e-3) {
            ssr_prev_iter = glob[m]
            datevec = datevec_iteration[., m]
            count = count + 1
            if (verbose) {
                printf(".")
                ndots = ndots + 1
                if (mod(ndots, 50) == 0) {
                    printf("\n")
                }
                displayflush()
            }
        }
        else {
            differe = -1
        }
    }

    // Computing using restrictions on the parameters
    count = 0
    ssr_est_prev = cross(y, y)

    contflag = 1
    while (contflag) {
        _ckp_est2_gls(y, model1, q, m, bigt, trm, datevec, dx, delta, ssr_est)

        if (count < maxiter & abs(ssr_est - ssr_est_prev) > 1e-3) {
            count = count + 1
            ssr_est_prev = ssr_est
            datevec = sort(dx, 1)
            if (verbose) {
                printf(".")
                ndots = ndots + 1
                if (mod(ndots, 50) == 0) {
                    printf("\n")
                }
                displayflush()
            }
        }
        else {
            contflag = 0
        }
    }
    if (verbose) {
        printf(" done)\n")
        displayflush()
    }

    // Construct the matrix of regressors from the final datevec
    z = _ckp_zbuild(t, model1, datevec)

    // common statistics tail; note that, exactly as in the GAUSS source, cbar
    // is the value from the GLS dating loop and is not recomputed after the
    // restricted estimation stage
    _ckp_tail(y, z, cbar, model1, penalty, kmax, kmin, stats, krule)

    min_tb = datevec
}

// ---------------------------------------------------------------------------
// Dispatcher, GAUSS proc(9)=sbur_multiple_gls, plus Stata interface
// ---------------------------------------------------------------------------
void _ckp_main(string scalar yname, string scalar touse, real scalar model1, real scalar known, string scalar tbmat, real scalar nbrk, real scalar penalty, real scalar kmax, real scalar kmin, real scalar method, real scalar maxiter, string scalar STname, string scalar CVname, string scalar TBname)
{
    real colvector y, min_tb, tbin, lam
    real rowvector stats
    real matrix CV
    real scalar T, cbar, krule

    y = st_data(., yname, touse)
    T = rows(y)

    tbin = st_matrix(tbmat)
    if (cols(tbin) > 1) {
        tbin = tbin'
    }

    // dispatcher exactly as in sbur_multiple_gls: models 0/1 -> brute;
    // models 2/3 -> brute (estimation=0) or algorithm (estimation=1);
    // known break dates are handled by the brute proc as in the source
    if (model1 == 0 | model1 == 1) {
        _ckp_brute(y, model1, 0, tbin, 0, penalty, kmax, kmin, stats, min_tb, cbar, krule)
    }
    else if (known == 1) {
        _ckp_brute(y, model1, 1, tbin, rows(tbin), penalty, kmax, kmin, stats, min_tb, cbar, krule)
    }
    else if (method == 0) {
        _ckp_brute(y, model1, 0, tbin, nbrk, penalty, kmax, kmin, stats, min_tb, cbar, krule)
    }
    else {
        _ckp_algorithm(y, model1, nbrk, penalty, kmax, kmin, maxiter, stats, min_tb, cbar, krule)
    }

    // critical values, exactly the sample.gss mechanics:
    // lam=(min_tb/T)|zeros(5-m,1); {cv_msb,cv_mza,cv_mzt,cv_pt}=msbur_rs(lam,cbar)
    if (model1 <= 1) {
        lam = J(5, 1, 0)
    }
    else {
        lam = min_tb:/T
        if (rows(lam) < 5) {
            lam = lam \ J(5-rows(lam), 1, 0)
        }
    }
    CV = _ckp_msbur_rs(lam, cbar)

    st_matrix(STname, stats)
    st_matrix(CVname, CV)
    st_matrix(TBname, min_tb)
    st_numscalar("r(cbar)", cbar)
    st_numscalar("r(krule)", krule)
}


// ---------------------------------------------------------------------
// Perron-Yabu (2009, JBES 27, 369-396) pretest for a shift in trend:
// Exp-W_FS. Line-by-line port of qfgls.prg (Perron and Yabu, version 2,
// March 2009), reusing the GAUSS-verbatim helpers above
// (_ckp_trimr, _ckp_lagn, _ckp_diff, _ckp_seqa, _ckp_minindc).
// ---------------------------------------------------------------------

// qfgls.prg proc IC: lag length by AIC (criteria=1) or BIC (criteria=2);
// the deterministic regressors are augmented with lagged LEVELS of x
real scalar _ckp_py_ic(real colvector x, real matrix reg_in, real scalar kmax, real scalar criteria)
{
    real matrix reg, rego
    real colvector ICV, dep, e, b
    real scalar i, khat

    reg = reg_in
    ICV = J(kmax+1, 1, 0)
    dep = _ckp_trimr(x, kmax, 0)
    rego = _ckp_trimr(reg, kmax, 0)
    e = dep - rego*luinv(rego'rego)*(rego'dep)
    ICV[1] = ln((e'e)/rows(e))
    for (i=1; i<=kmax; i++) {
        reg = reg, _ckp_lagn(x, i)
        rego = _ckp_trimr(reg, kmax, 0)
        b = qrsolve(rego, dep)
        e = dep - rego*b
        if (criteria == 1) {
            ICV[i+1] = ln((e'e)/rows(e)) + (2*i)/rows(e)
        }
        else {
            ICV[i+1] = ln((e'e)/rows(e)) + (ln(rows(e))*i)/rows(e)
        }
    }
    khat = _ckp_minindc(ICV) - 1
    khat = khat*(khat>=1) + 1*(khat==0)
    return(khat)
}

// qfgls.prg proc ACV: demeaned autocovariances R(0),...,R(T-1), divisor T
real colvector _ckp_py_acv(real colvector x)
{
    real scalar T, i, j, xbar
    real colvector R

    T = rows(x)
    R = J(T, 1, 0)
    xbar = mean(x)
    for (i=1; i<=T; i++) {
        j = i - 1
        R[i] = ((x[|1,1 \ T-j,1|] :- xbar)'(x[|1+j,1 \ T,1|] :- xbar))/T
    }
    return(R)
}

// qfgls.prg proc h0W: Andrews (1991) long-run variance, Quadratic
// Spectral window with AR(1) plug-in bandwidth (no intercept in the
// plug-in regression, as in the GAUSS code)
real scalar _ckp_py_h0w(real colvector x)
{
    real scalar T, b, a, m, h0
    real colvector R, s, delta, lamda

    T = rows(x)
    b = qrsolve(x[|1,1 \ T-1,1|], x[|2,1 \ T,1|])
    a = (4*(b^2))/((1-b)^4)
    R = _ckp_py_acv(x)
    s = _ckp_seqa(1, 1, T-1)
    m = 1.3221*((a*T)^(1/5))
    delta = (6*pi()*s)/(5*m)
    lamda = 3*((sin(delta):/delta) - cos(delta)) :/ (delta:^2)
    h0 = R[1] + 2*(lamda'R[|2,1 \ T,1|])
    return(h0)
}

// qfgls.prg proc breakdate: OLS minimum-SSR single break date over the
// same trimmed range as the main loop
real scalar _ckp_py_breakdate(real colvector x, real scalar model, real scalar kmax, real scalar eps)
{
    real scalar T, TB0, TBN, TBi
    real colvector vsum, constant, trend, DUi, DTi, e
    real matrix reg

    T = rows(x)
    constant = J(T, 1, 1)
    trend = _ckp_seqa(1, 1, T)
    TB0 = max((trunc(eps*T), kmax+2))
    TBN = trunc((1-eps)*T)
    vsum = J(TBN-TB0+1, 1, 1000)
    for (TBi=TB0; TBi<=TBN; TBi++) {
        DUi = (trend :> TBi)
        DTi = (trend :> TBi) :* (trend :- TBi)
        if (model == 1) {
            reg = constant, trend, DUi
        }
        else if (model == 2) {
            reg = constant, trend, DTi
        }
        else {
            reg = constant, DUi, trend, DTi
        }
        e = (I(T) - reg*luinv(reg'reg)*reg')*x
        vsum[TBi-TB0+1] = e'e
    }
    return(TB0 + _ckp_minindc(vsum) - 1)
}

// main body of qfgls.prg: Exp-W_FS over the trimmed candidate range;
// writes (wald, TB, cv10, cv5, cv1) to OUTname
void _ckp_py(string scalar yname, string scalar touse, real scalar model, real scalar criteria, real scalar eps, real scalar kmax, string scalar OUTname)
{
    real colvector y, v_t, constant, trend, vect1, DUi, DTi, u, du, depu, gdep, v, depv, beta, e, b, ehat, DUki, DTki, gdepki
    real matrix cv, VR, reg, regu, greg, VCV, regv, vbeta, regki, gregki
    real scalar T, TBi, lam1, khat, j, ahat, vahat, tau1, t05, IP, r, k, c1, c2, ctau, rhomd1, amu, CR, amus, h0, ki, sige, wald, TB, ie

    y = st_data(., yname, touse)
    T = rows(y)

    // v_t and critical values embedded in qfgls.prg, by model
    if (model == 1) {
        VR = (0, 0, 1)
        v_t = (-4.30 \ -4.39 \ -4.39 \ -4.34 \ -4.32 \ -4.45 \ -4.42 \ -4.33 \ -4.27 \ -4.27)
        cv = (1.60, 2.07, 3.33 \ 1.52, 1.97, 3.24 \ 1.41, 1.88, 3.05 \ 1.26, 1.74, 3.12 \ 0.91, 1.33, 2.83)
    }
    else if (model == 2) {
        VR = (0, 0, 1)
        v_t = (-4.27 \ -4.41 \ -4.51 \ -4.55 \ -4.56 \ -4.57 \ -4.51 \ -4.38 \ -4.26 \ -4.26)
        cv = (1.52, 2.02, 3.37 \ 1.40, 1.93, 3.27 \ 1.28, 1.86, 3.20 \ 1.13, 1.67, 3.06 \ 0.74, 1.28, 2.61)
    }
    else {
        VR = (0, 1, 0, 0 \ 0, 0, 0, 1)
        v_t = (-4.38 \ -4.65 \ -4.78 \ -4.81 \ -4.90 \ -4.88 \ -4.75 \ -4.70 \ -4.41 \ -4.41)
        cv = (2.96, 3.55, 5.02 \ 2.82, 3.36, 4.78 \ 2.65, 3.16, 4.59 \ 2.48, 3.12, 4.47 \ 2.15, 2.79, 4.57)
    }

    constant = J(T, 1, 1)
    trend = _ckp_seqa(1, 1, T)
    vect1 = J(trunc((1-2*eps)*T)+2, 1, 0)

    for (TBi = max((trunc(eps*T), kmax+2)); TBi <= trunc((1-eps)*T); TBi++) {

        lam1 = TBi/T
        DUi = (trend :> TBi)
        DTi = (trend :> TBi) :* (trend :- TBi)
        if (model == 1) {
            reg = constant, trend, DUi
        }
        else if (model == 2) {
            reg = constant, trend, DTi
        }
        else {
            reg = constant, DUi, trend, DTi
        }

        // ---- estimation of alpha ----
        khat = _ckp_py_ic(y, reg, kmax, criteria)
        u = (I(T) - reg*luinv(reg'reg)*reg')*y
        du = _ckp_diff(u, 1)
        depu = u
        regu = _ckp_lagn(u, 1)
        for (j=1; j<=khat-1; j++) {
            regu = regu, _ckp_lagn(du, j)
        }
        depu = _ckp_trimr(depu, khat, 0)
        regu = _ckp_trimr(regu, khat, 0)
        b = luinv(regu'regu)*(regu'depu)
        ehat = depu - regu*b
        VCV = ((ehat'ehat)/rows(ehat))*luinv(regu'regu)
        ahat = b[1]
        vahat = VCV[1,1]
        tau1 = (ahat-1)/sqrt(vahat)

        // ---- upper-biased estimator ----
        if (lam1 <= 0.1) {
            t05 = v_t[1]
        }
        else if (lam1 <= 0.2) {
            t05 = v_t[2]
        }
        else if (lam1 <= 0.3) {
            t05 = v_t[3]
        }
        else if (lam1 <= 0.4) {
            t05 = v_t[4]
        }
        else if (lam1 <= 0.5) {
            t05 = v_t[5]
        }
        else if (lam1 <= 0.6) {
            t05 = v_t[6]
        }
        else if (lam1 <= 0.7) {
            t05 = v_t[7]
        }
        else if (lam1 <= 0.8) {
            t05 = v_t[8]
        }
        else if (lam1 <= 0.9) {
            t05 = v_t[9]
        }
        else {
            t05 = v_t[10]
        }
        IP = trunc((khat+1)/2)
        r = cols(reg)
        k = 10
        c1 = sqrt((1+r)*T)
        c2 = ((r+1)*T - t05*t05*(IP+T))/(t05*(t05+k)*(IP+T))
        if (tau1 > t05) {
            ctau = -tau1
        }
        else if (tau1 > -k) {
            ctau = IP*(tau1/T) - (r+1)/(tau1 + c2*(tau1+k))
        }
        else if (tau1 > -c1) {
            ctau = IP*(tau1/T) - (r+1)/tau1
        }
        else {
            ctau = 0
        }
        rhomd1 = ahat + ctau*sqrt(vahat)
        rhomd1 = 1*(rhomd1>=1) + rhomd1*(abs(rhomd1)<1) - 0.99*(rhomd1<=-1)
        amu = rhomd1

        // ---- super-efficient truncation and quasi-GLS ----
        CR = (T^0.5)*abs(amu-1)
        amus = amu*(CR>1) + 1*(CR<=1)
        gdep = y[1] \ (_ckp_trimr(y, 1, 0) - amus*_ckp_trimr(y, 0, 1))
        greg = reg[1,.] \ (_ckp_trimr(reg, 1, 0) - amus*_ckp_trimr(reg, 0, 1))
        b = luinv(greg'greg)*(greg'gdep)
        v = gdep - greg*b

        // ---- h(0) ----
        if (khat == 1) {
            h0 = (v'v)/rows(v)
        }
        else {
            if (amus == 1) {
                for (ki=1; ki<=khat-1; ki++) {
                    if (ki == 1) {
                        regv = _ckp_lagn(v, ki)
                    }
                    else {
                        regv = regv, _ckp_lagn(v, ki)
                    }
                }
                depv = _ckp_trimr(v, khat-1, 0)
                regv = _ckp_trimr(regv, khat-1, 0)
                beta = luinv(regv'regv)*(regv'depv)
                e = depv - regv*beta
                if (model == 1) {
                    vbeta = J(3, khat-1, 0)
                    for (ki=1; ki<=khat-1; ki++) {
                        DUki = (trend :> TBi-ki)
                        regki = constant, trend, DUki
                        gdepki = y[1] \ (_ckp_trimr(y, 1, 0) - amus*_ckp_trimr(y, 0, 1))
                        gregki = reg[1,.] \ (_ckp_trimr(regki, 1, 0) - amus*_ckp_trimr(regki, 0, 1))
                        vbeta[.,ki] = luinv(gregki'gregki)*(gregki'gdepki)
                    }
                    b[3] = b[3] - vbeta[3,.]*beta
                    h0 = (e'e)/(T-khat)
                }
                else if (model == 2) {
                    h0 = ((e'e)/(T-khat))/((1 - sum(beta))^2)
                }
                else {
                    vbeta = J(4, khat-1, 0)
                    for (ki=1; ki<=khat-1; ki++) {
                        DUki = (trend :> TBi-ki)
                        DTki = (trend :> TBi-ki) :* (trend :- TBi)
                        regki = constant, DUki, trend, DTki
                        gdepki = y[1] \ (_ckp_trimr(y, 1, 0) - amus*_ckp_trimr(y, 0, 1))
                        gregki = reg[1,.] \ (_ckp_trimr(regki, 1, 0) - amus*_ckp_trimr(regki, 0, 1))
                        vbeta[.,ki] = luinv(gregki'gregki)*(gregki'gdepki)
                    }
                    sige = (e'e)/(T-khat)
                    h0 = sige/((1 - sum(beta))^2)
                    b[2] = sqrt(h0)*(b[2] - vbeta[2,.]*beta)/sqrt(sige)
                }
            }
            else if (abs(amus) < 1) {
                h0 = _ckp_py_h0w(v)
            }
        }

        VCV = h0*luinv(greg'greg)
        vect1[TBi - trunc(eps*T) + 1] = (VR*b)'*luinv(VR*VCV*VR')*(VR*b)
    }

    // Exp functional exactly as in qfgls.prg (including the division by T)
    wald = ln(sum(exp(vect1:/2))/T)
    TB = _ckp_py_breakdate(y, model, kmax, eps)

    if (eps == 0.01) {
        ie = 1
    }
    else if (eps == 0.05) {
        ie = 2
    }
    else if (eps == 0.10) {
        ie = 3
    }
    else if (eps == 0.15) {
        ie = 4
    }
    else {
        ie = 5
    }
    st_matrix(OUTname, (wald, TB, cv[ie,1], cv[ie,2], cv[ie,3]))
}


// ---------------------------------------------------------------------
// Breusch-Godfrey LM machinery for the spectral AR regression, exactly
// as in narayanp (cross-validated there against estat bgodfrey, nomiss0):
// the auxiliary regression for order m drops the first m observations
// (no zero-filling), has no constant (the spectral regression has none,
// so the uncentered R-squared is used), LM = (n-m)*R2 ~ chi2(m), and its
// residual degrees of freedom are df = n - k - 2m.
// ---------------------------------------------------------------------

// LM p-value for a single BG order p on residuals ehat of the regression
// of dep on X (no constant), nomiss0 variant
real scalar _ckp_bgp(real colvector ehat, real matrix X, real scalar p)
{
    real matrix XA
    real colvector dep, ba, ea
    real scalar j, n, na, r2, lm

    n = rows(ehat)
    na = n - p
    dep = ehat[|p+1,1 \ n,1|]
    XA = X[|p+1,1 \ n,.|]
    for (j=1; j<=p; j++) {
        XA = XA, ehat[|p+1-j,1 \ n-j,1|]
    }
    ba = qrsolve(XA, dep)
    ea = dep - XA*ba
    r2 = 1 - cross(ea, ea)/cross(dep, dep)
    lm = na*r2
    return(chi2tail(p, lm))
}

// minimum BG p-value over orders 1..H on the adfp regression at lag kstar
// (orders with df = n - cols(X) - p < dffloor are skipped); returns
// (minp, minorder), missing if no admissible order
void _ckp_bgmin(real colvector yt, real scalar kstar, real scalar H, real scalar dffloor, real scalar minp, real scalar mino)
{
    real colvector dyt, ehat, b
    real matrix reg
    real scalar i, p, n, pv

    reg = _ckp_lagn(yt, 1)
    dyt = _ckp_diff(yt, 1)
    i = 1
    while (i <= kstar) {
        reg = reg, _ckp_lagn(dyt, i)
        i = i + 1
    }
    dyt = _ckp_trimr(dyt, kstar+1, 0)
    reg = _ckp_trimr(reg, kstar+1, 0)
    _ckp_olsqr2(dyt, reg, b, ehat)

    n = rows(dyt)
    minp = .
    mino = .
    for (p=1; p<=H; p++) {
        if (n - cols(reg) - 2*p < dffloor) {
            continue
        }
        pv = _ckp_bgp(ehat, reg, p)
        if (minp >= . | pv < minp) {
            minp = pv
            mino = p
        }
    }
}

// general-to-specific BG lag selection: start at kmax, walk down while
// the regression stays clean (min BG p >= 0.05 over orders 1..H); at the
// first dirty lag return to the previous clean one. If kmax itself is
// dirty, kmax is used and the warning flag is set.
real scalar _ckp_bgsel(real colvector yt, real scalar kmax, real scalar kmin, real scalar H, real scalar dffloor, real scalar warn)
{
    real scalar k, minp, mino

    warn = 0
    _ckp_bgmin(yt, kmax, H, dffloor, minp, mino)
    if (minp < 0.05) {
        warn = 1
        return(kmax)
    }
    k = kmax
    while (k > kmin) {
        _ckp_bgmin(yt, k-1, H, dffloor, minp, mino)
        if (minp < 0.05) {
            return(k)
        }
        k = k - 1
    }
    return(kmin)
}

end
*==============================================================================
* End of ckptest.ado
*==============================================================================