*! version 1.5.9  06aug2026  Ozan Eruygur
*! lumpapell: Lumsdaine-Papell (1997) unit root test with structural breaks
*! Port of the RATS procedure lpunit.src (Tom Doan, Estima, revision 05/2017)
*! plus an article mode replicating Lumsdaine-Papell (1997): model CA, the
*! article grid, per-combination general-to-specific lag selection, and the
*! article finite-sample critical values
*! two structural breaks at unknown dates

program define lumpapell, rclass
    version 11.0
    syntax varname(numeric ts) [if] [in] [, Break(string) NBreaks(integer 2) Method(string) Lags(integer -1) MAXLags(integer -1) SIGnif(real 0.10) PI(real 0.15) CV(string) TCrit(real -1) PAPER RATSOut Graph noPRint Title(string asis)]

    * ------------------------------------------------------------------
    * Parse break() : intercept (default) / trend / both / ca -> 1/2/3/4
    * ------------------------------------------------------------------
    local brk = lower(trim(`"`break'"'))
    if `"`brk'"'=="" local brk "intercept"
    local breakn 0
    if `"`brk'"'=="intercept" local breakn 1
    if `"`brk'"'=="trend" local breakn 2
    if `"`brk'"'=="both" local breakn 3
    if `"`brk'"'=="ca" local breakn 4
    if `breakn'==0 {
        di as error "break() must be intercept, trend, both, or ca"
        exit 198
    }

    * ------------------------------------------------------------------
    * Parse method() : input (default) / aic / bic / hq / ttest / gtos
    * gtos and ttest are identical (as in the RATS code, lmethod 6 -> 5)
    * ------------------------------------------------------------------
    local mth = lower(trim(`"`method'"'))
    if `"`mth'"'=="" local mth "input"
    local lmethod 0
    if `"`mth'"'=="input" local lmethod 1
    if `"`mth'"'=="aic" local lmethod 2
    if `"`mth'"'=="bic" local lmethod 3
    if `"`mth'"'=="hq" local lmethod 4
    if `"`mth'"'=="ttest" local lmethod 5
    if `"`mth'"'=="gtos" local lmethod 5
    if `lmethod'==0 {
        di as error "method() must be input, aic, bic, hq, ttest, or gtos"
        exit 198
    }

    if `nbreaks'<0 {
        di as error "nbreaks() must be a nonnegative integer"
        exit 198
    }
    if `pi'<=0 | `pi'>=1 {
        di as error "pi() must be strictly between 0 and 1"
        exit 198
    }
    if `signif'<=0 | `signif'>=1 {
        di as error "signif() must be strictly between 0 and 1"
        exit 198
    }
    if `lags'<-1 | `maxlags'<-1 {
        di as error "lags() and maxlags() must be nonnegative integers"
        exit 198
    }

    * ------------------------------------------------------------------
    * Article mode switches: paper, cv(), tcrit(), break(ca)
    * ------------------------------------------------------------------
    local papermode = cond("`paper'"!="", 1, 0)
    if `breakn'==4 & `papermode'==0 {
        di as error "break(ca) requires the paper option"
        exit 198
    }
    if `breakn'==4 & `nbreaks'!=2 {
        di as error "break(ca) requires nbreaks(2)"
        exit 198
    }
    local cvsrc = lower(trim(`"`cv'"'))
    if `papermode'==1 & `"`cvsrc'"'=="" local cvsrc "paper"
    if `"`cvsrc'"'=="" local cvsrc "rats"
    if `"`cvsrc'"'!="rats" & `"`cvsrc'"'!="paper" {
        di as error "cv() must be rats or paper"
        exit 198
    }
    if `"`cvsrc'"'=="paper" {
        if `nbreaks'!=2 | `breakn'==2 {
            di as error "cv(paper) is available only for nbreaks(2) with break(intercept), break(ca), or break(both)"
            exit 198
        }
    }
    if `breakn'==4 & `"`cvsrc'"'=="rats" {
        di as error "no RATS critical values exist for model CA; use cv(paper)"
        exit 198
    }
    if `tcrit'!=-1 & `tcrit'<=0 {
        di as error "tcrit() must be a positive number"
        exit 198
    }
    if `tcrit'>0 & `tcrit'<. & `lmethod'!=5 {
        di as error "tcrit() may be combined only with method(ttest) or method(gtos)"
        exit 198
    }
    if `papermode'==1 & `lmethod'==5 & `tcrit'==-1 local tcrit = 1.60

    * ------------------------------------------------------------------
    * Time-series settings
    * ------------------------------------------------------------------
    capture tsset
    if _rc {
        di as error "the data must be tsset before using lumpapell"
        exit 111
    }
    if "`r(panelvar)'"!="" {
        di as error "lumpapell works on a single time series; panel data are not allowed"
        exit 198
    }
    local tvar `r(timevar)'
    local dlt = r(tdelta)
    local tsfmt "`r(tsfmt)'"
    if `dlt'>=. local dlt 1

    tsrevar `varlist'
    local tv `r(varlist)'

    * entry numbering: entry 1 = first observation in the dataset
    * (matches a RATS workspace whose CALENDAR begins at the start of the data)
    quietly summarize `tvar', meanonly
    local tmin = r(min)

    * ------------------------------------------------------------------
    * [in] maps to the RATS start/end parameters (range restriction)
    * [if] maps to the RATS SMPL option (rows skipped inside the range)
    * ------------------------------------------------------------------
    tempvar inrng smplv okv wsel
    quietly gen byte `inrng' = 0
    quietly replace `inrng' = 1 `in'
    quietly gen byte `smplv' = 0
    quietly replace `smplv' = 1 `if'

    * range available for the series and one lag (RATS: inquire(reglist) # series{0 1})
    quietly gen byte `okv' = `inrng' & !missing(`tv') & !missing(L.`tv')
    quietly count if `okv'
    if r(N)==0 {
        di as error "no observations with both the series and its first lag available"
        exit 2000
    }
    quietly summarize `tvar' if `okv', meanonly
    local startl = round((r(min)-`tmin')/`dlt') + 1
    local endl = round((r(max)-`tmin')/`dlt') + 1

    * the time grid must be gap-free on [startl-1, endl]
    local t_lo = `tmin' + (`startl'-2)*`dlt'
    local t_hi = `tmin' + (`endl'-1)*`dlt'
    quietly count if `tvar'>=`t_lo' & `tvar'<=`t_hi'
    if r(N) != `endl'-`startl'+2 {
        di as error "gaps found in the time variable within the estimation window; run tsfill first"
        exit 498
    }
    quietly gen byte `wsel' = `tvar'>=`t_lo' & `tvar'<=`t_hi'

    * ------------------------------------------------------------------
    * maxlag resolution (RATS lines 135-141) and the raw MAXLAGS value
    * used for pinobs (RATS line 124 uses the option value, 0 if not given)
    * ------------------------------------------------------------------
    local maxlag0 = -1
    if `lags'>=0 local maxlag0 = `lags'
    else if `maxlags'>=0 local maxlag0 = `maxlags'
    local mlraw = cond(`maxlags'>=0, `maxlags', 0)

    * ------------------------------------------------------------------
    * Core computation in Mata
    * ------------------------------------------------------------------
    mata: lumpapell_main("`tv'", "`tvar'", "`smplv'", "`wsel'", `tmin', `dlt', `startl', `endl', `breakn', `nbreaks', `lmethod', `maxlag0', `mlraw', `signif', `pi', `papermode', `tcrit')

    local mrc = __lp_rc
    if `mrc'==1 {
        capture scalar drop __lp_rc
        di as error "too few usable observations for the requested specification"
        exit 2001
    }
    if `mrc'==2 {
        capture scalar drop __lp_rc
        di as error "pi() and the sample length imply a trimming window of zero observations; increase pi() or the sample"
        exit 498
    }
    if `mrc'==3 {
        capture scalar drop __lp_rc
        di as error "too few observations for nbreaks() breaks with the requested grid"
        exit 2001
    }

    tempname smint
    scalar `smint' = __lp_mint
    local bestlag = __lp_bestlag
    local maxlag = __lp_maxlag
    local pinobs = __lp_pinobs
    local N = __lp_N
    local K = __lp_K
    tempname bmat tmat bpsmat
    matrix `bmat' = __lp_b
    matrix `tmat' = __lp_t
    if `nbreaks'>0 matrix `bpsmat' = __lp_bps
    capture scalar drop __lp_rc __lp_mint __lp_bestlag __lp_maxlag __lp_pinobs __lp_N __lp_K
    capture matrix drop __lp_b __lp_t __lp_bps

    * ------------------------------------------------------------------
    * Break dates as time values and display labels
    * ------------------------------------------------------------------
    local dlist ""
    local blabels ""
    if `nbreaks'>0 {
        forvalues i = 1/`nbreaks' {
            local e_i = `bpsmat'[1,`i']
            local tv_`i' = `tmin' + (`e_i'-1)*`dlt'
            local dl_`i' = trim(`"`: display `tsfmt' `tv_`i'''"')
            local dlist "`dlist' `dl_`i''"
            local blabels "`blabels'`dl_`i'' "
        }
    }

    * ------------------------------------------------------------------
    * Break dummies and trend written to the dataset (lumpapell_ prefix)
    * and the manual regress command that reproduces the break regression
    * ------------------------------------------------------------------
    capture drop lumpapell_*
    tempvar entv
    quietly gen double `entv' = round((`tvar'-`tmin')/`dlt') + 1
    quietly gen int lumpapell_trend = `entv'
    label variable lumpapell_trend "lumpapell: trend (entry number, first obs = 1)"
    local dvars ""
    if `nbreaks'>0 {
        forvalues i = 1/`nbreaks' {
            local e_i = `bpsmat'[1,`i']
            local mkdu = (`breakn'==1 | `breakn'==3 | `breakn'==4)
            local mkdt = (`breakn'==2 | `breakn'==3 | (`breakn'==4 & `i'==1))
            if `mkdu' {
                quietly gen byte lumpapell_du`i' = `entv' > `e_i'
                label variable lumpapell_du`i' "lumpapell: D(`dl_`i'')"
                local dvars "`dvars' lumpapell_du`i'"
            }
            if `mkdt' {
                quietly gen int lumpapell_dt`i' = (`entv' - `e_i') * (`entv' > `e_i')
                label variable lumpapell_dt`i' "lumpapell: DT(`dl_`i'')"
                local dvars "`dvars' lumpapell_dt`i'"
            }
        }
    }
    local fromt2 = `tmin' + (`startl'+`bestlag'-1)*`dlt'
    local endt2 = `tmin' + (`endl'-1)*`dlt'
    local lagpart ""
    if `bestlag'>0 local lagpart " L(1/`bestlag').D.`varlist'"
    * sample condition needed only when the user restricted with [in] or [if];
    * on a full series the lag operators already reproduce the sample
    local conds ""
    if `"`in'"'!="" {
        local fam = substr("`tsfmt'",1,3)
        local tfun ""
        if "`fam'"=="%tq" local tfun "tq"
        if "`fam'"=="%tm" local tfun "tm"
        if "`fam'"=="%tw" local tfun "tw"
        if "`fam'"=="%th" local tfun "th"
        if "`fam'"=="%td" local tfun "td"
        if "`tfun'"!="" {
            local fl = trim(`"`: display `fam' `fromt2''"')
            local el = trim(`"`: display `fam' `endt2''"')
            local conds "`tvar'>=`tfun'(`fl') & `tvar'<=`tfun'(`el')"
        }
        else local conds "`tvar'>=`fromt2' & `tvar'<=`endt2'"
    }
    if `"`if'"'!="" {
        local uexp0 = subinstr(`"`if'"', "if ", "", 1)
        if `"`conds'"'=="" local conds "(`uexp0')"
        else local conds "`conds' & (`uexp0')"
    }
    if `"`conds'"'=="" local regcmd `"regress D.`varlist' L.`varlist'`dvars' lumpapell_trend`lagpart'"'
    else local regcmd `"regress D.`varlist' L.`varlist'`dvars' lumpapell_trend`lagpart' if `conds'"'

    * ------------------------------------------------------------------
    * Critical values: RATS table (default) or the article's finite-sample
    * table (cv(paper); T=125, includes the 2.5 percent level)
    * ------------------------------------------------------------------
    tempname cv1 cv25 cv5 cv10
    scalar `cv1' = .
    scalar `cv25' = .
    scalar `cv5' = .
    scalar `cv10' = .
    if `"`cvsrc'"'=="paper" {
        if `breakn'==1 {
            scalar `cv1' = -6.94
            scalar `cv25' = -6.53
            scalar `cv5' = -6.24
            scalar `cv10' = -5.96
        }
        if `breakn'==4 {
            scalar `cv1' = -7.24
            scalar `cv25' = -7.02
            scalar `cv5' = -6.65
            scalar `cv10' = -6.33
        }
        if `breakn'==3 {
            scalar `cv1' = -7.34
            scalar `cv25' = -7.02
            scalar `cv5' = -6.82
            scalar `cv10' = -6.49
        }
    }
    else if `nbreaks'==0 {
        local n0 = `N' + 1
        scalar `cv1' = -3.9638 - 8.353/`n0' - 47.44/(`n0'^2)
        scalar `cv5' = -3.4126 - 4.039/`n0' - 17.83/(`n0'^2)
        scalar `cv10' = -3.1279 - 2.418/`n0' - 7.58/(`n0'^2)
    }
    else if `nbreaks'==1 {
        if `breakn'==1 {
            scalar `cv1' = -5.34
            scalar `cv5' = -4.80
            scalar `cv10' = -4.58
        }
        if `breakn'==2 {
            scalar `cv1' = -4.93
            scalar `cv5' = -4.42
            scalar `cv10' = -4.11
        }
        if `breakn'==3 {
            scalar `cv1' = -5.57
            scalar `cv5' = -5.08
            scalar `cv10' = -4.82
        }
    }
    else if `nbreaks'==2 {
        if `breakn'==1 {
            scalar `cv1' = -6.74
            scalar `cv5' = -6.16
            scalar `cv10' = -5.89
        }
        if `breakn'==2 {
            scalar `cv1' = -7.19
            scalar `cv5' = -6.62
            scalar `cv10' = -6.37
        }
        if `breakn'==3 {
            scalar `cv1' = -7.19
            scalar `cv5' = -6.75
            scalar `cv10' = -6.48
        }
    }

    * significance tag: ** below the 1% value, * below the 5% value
    local stars ""
    if `cv5'<. & `smint'<`cv5' local stars "*"
    if `cv1'<. & `smint'<`cv1' local stars "**"

    * ------------------------------------------------------------------
    * Graph (RATS lines 254-257; drawn regardless of print, as in RATS)
    * ------------------------------------------------------------------
    if "`graph'"!="" {
        twoway tsline `tv', tline(`dlist') note("Series with Breaks Highlighted") name(lumpapell, replace)
    }

    * ------------------------------------------------------------------
    * Coefficient labels: Y{1}, then break dummies, Constant, Trend, lags
    * Model CA order: D(TB1) DT(TB1) D(TB2), as in Table 4 of the article
    * ------------------------------------------------------------------
    local vlabs `""Y{c -(}1{c )-}""'
    local cnames "Y1"
    if `breakn'==4 {
        local vlabs `"`vlabs' "D(`dl_1')" "DT(`dl_1')" "D(`dl_2')""'
        local cnames "`cnames' D1 DT1 D2"
    }
    else if `nbreaks'>0 {
        forvalues i = 1/`nbreaks' {
            if `breakn'!=2 {
                local vlabs `"`vlabs' "D(`dl_`i'')""'
                local cnames "`cnames' D`i'"
            }
            if `breakn'!=1 {
                local vlabs `"`vlabs' "DT(`dl_`i'')""'
                local cnames "`cnames' DT`i'"
            }
        }
    }
    local vlabs `"`vlabs' "Constant" "Trend""'
    local cnames "`cnames' Constant Trend"
    if `bestlag'>0 {
        forvalues j = 1/`bestlag' {
            local vlabs `"`vlabs' "dy{c -(}`j'{c )-}""'
            local cnames "`cnames' dyL`j'"
        }
    }
    matrix colnames `bmat' = `cnames'
    matrix colnames `tmat' = `cnames'

    * readable labels for the default report: lags as y(-1), dy(-1), dy(-2)
    local vlabs2 `""y(-1)""'
    if `breakn'==4 {
        local vlabs2 `"`vlabs2' "D(`dl_1')" "DT(`dl_1')" "D(`dl_2')""'
    }
    else if `nbreaks'>0 {
        forvalues i = 1/`nbreaks' {
            if `breakn'!=2 local vlabs2 `"`vlabs2' "D(`dl_`i'')""'
            if `breakn'!=1 local vlabs2 `"`vlabs2' "DT(`dl_`i'')""'
        }
    }
    local vlabs2 `"`vlabs2' "Constant" "Trend""'
    if `bestlag'>0 {
        forvalues j = 1/`bestlag' {
            local vlabs2 `"`vlabs2' "dy(-`j')""'
        }
    }

    * ------------------------------------------------------------------
    * Report: default format below; ratsout reproduces the RATS layout
    * ------------------------------------------------------------------
    if "`ratsout'"!="" {
        if "`print'"!="noprint" {
            if `"`title'"'=="" local title "Lumsdaine-Papell Unit Root Test, Series `varlist'"
            local frome = `startl' + `bestlag' + 1
            local fromt = `tmin' + (`frome'-1)*`dlt'
            local endt = `tmin' + (`endl'-1)*`dlt'
            local froml = trim(`"`: display `tsfmt' `fromt''"')
            local endll = trim(`"`: display `tsfmt' `endt''"')
            di as text ""
            di as text `"`title'"'
            di as text "Regression Run From `froml' to `endll'"
            di as text "Observations " as result `N'
            if `nbreaks'==0 local descript "None -- Dickey-Fuller Test"
            if `nbreaks'>0 & `breakn'==1 local descript "Intercept Only"
            if `nbreaks'>0 & `breakn'==2 local descript "Trend Only"
            if `nbreaks'>0 & `breakn'==3 local descript "Intercept and Trend"
            if `breakn'==4 local descript "Model CA (Intercept and Trend; Intercept Only)"
            di as text "Breaks in `descript'"
            if `nbreaks'>0 di as text "Breaks at `blabels'"
            if `lmethod'==1 di as text "Using fixed lags " as result `bestlag'
            else di as text "With " as result `bestlag' as text " lags chosen from " as result `maxlag'
            if `lmethod'==1 local llagmethod "User"
            if `lmethod'==2 local llagmethod "AIC"
            if `lmethod'==3 local llagmethod "BIC"
            if `lmethod'==4 local llagmethod "HQ"
            if `lmethod'==5 {
                if `tcrit'>0 & `tcrit'<. {
                    local tcs = trim(`"`: display %5.2f `tcrit''"')
                    local llagmethod "GTOS/t-tests(|t|>=`tcs')"
                }
                else {
                    local sg = trim(`"`: display %5.2f `signif''"')
                    local llagmethod "GTOS/t-tests(`sg')"
                }
            }
            di as text "Selected by `llagmethod'"
            di as text "Sig Level" _col(13) %14s "Crit Value"
            if `cv1'<. di as text "1%(**)" _col(13) as result %14.4f `cv1'
            else di as text "1%(**)" _col(13) %14s "NA"
            if `cv25'<. di as text "2.5%" _col(13) as result %14.4f `cv25'
            if `cv5'<. di as text "5%(*)" _col(13) as result %14.4f `cv5'
            else di as text "5%(*)" _col(13) %14s "NA"
            if `cv10'<. di as text "10%" _col(13) as result %14.4f `cv10'
            else di as text "10%" _col(13) %14s "NA"
            di as text ""
            di as text "Variable" _col(15) %14s "Coefficient" %14s "T-Stat"
            local pos 1
            foreach lab of local vlabs {
                local bb = `bmat'[1,`pos']
                local tt = `tmat'[1,`pos']
                if `pos'==1 di as text `"`lab'"' _col(15) as result %14.4f `bb' %14.4f `tt' " `stars'"
                else di as text `"`lab'"' _col(15) as result %14.4f `bb' %14.4f `tt'
                local pos = `pos' + 1
            }
            di as text ""
        }
    }
    else if "`print'"!="noprint" {
        local endt = `tmin' + (`endl'-1)*`dlt'
        local endll = trim(`"`: display `tsfmt' `endt''"')
        local truee = `startl' + `bestlag'
        local truet = `tmin' + (`truee'-1)*`dlt'
        local truel = trim(`"`: display `tsfmt' `truet''"')
        if `nbreaks'==0 local mdesc "Augmented Dickey-Fuller, no breaks"
        if `nbreaks'>0 & `breakn'==1 local mdesc "Intercept breaks (DU)"
        if `nbreaks'==2 & `breakn'==1 local mdesc "AA - both breaks in the intercept (DU)"
        if `nbreaks'>0 & `breakn'==2 local mdesc "Trend breaks (DT)"
        if `nbreaks'>0 & `breakn'==3 local mdesc "Intercept and trend breaks (DU and DT)"
        if `nbreaks'==2 & `breakn'==3 local mdesc "CC - both breaks in intercept and trend (DU and DT)"
        if `breakn'==4 local mdesc "CA - TB1: intercept and trend; TB2: intercept only"
        local modetxt = cond(`papermode'==1, "article (Lumsdaine-Papell 1997)", "RATS (lpunit.src)")
        local pis = trim(`"`: display %5.2f `pi''"')
        if `papermode'==1 local gdesc "article grid (2nd to next-to-last obs, min gap 2)"
        else local gdesc "RATS grid (pi = `pis', min gap `pinobs' obs)"
        if `lmethod'==1 local ldesc "Fixed (input)"
        if `lmethod'==2 local ldesc "AIC"
        if `lmethod'==3 local ldesc "BIC"
        if `lmethod'==4 local ldesc "Hannan-Quinn"
        if `lmethod'==5 {
            local tcs = trim(`"`: display %5.2f `tcrit''"')
            local sg = trim(`"`: display %5.2f `signif''"')
            if `tcrit'>0 & `tcrit'<. local ldesc "Sequential t-test (|t| >= `tcs')"
            else local ldesc "Sequential t-test (p <= `sg')"
        }
        di as text ""
        if `"`title'"'!="" di as text `"`title'"'
        di as text "Lumsdaine-Papell test for unit root" _col(46) "Number of obs    = " as result %6.0f `N'
        di as text "Variable: " as result "`varlist'" as text _col(46) "Number of breaks = " as result %6.0f `nbreaks'
        di as text ""
        di as text "  Model" _col(18) ": " as result "`mdesc'"
        di as text "  Mode" _col(18) ": " as result "`modetxt'"
        di as text "  Break grid" _col(18) ": " as result "`gdesc'"
        di as text "  Lag selection" _col(18) ": " as result "`ldesc'"
        di as text "  Selected lag" _col(18) ": " as result `bestlag' as text "  (maximum " as result `maxlag' as text ")"
        di as text "  Sample" _col(18) ": " as result "`truel' - `endll'"
        di as text ""
        di as text "H0: " as result "`varlist'" as text " has a unit root"
        di as text ""
        if `cv25'<. {
            di as text _col(30) "Lumsdaine-Papell"
            di as text _col(12) "Test" _col(22) "----------- critical value -----------"
            di as text _col(9) "statistic" _col(27) "1%" _col(35) "2.5%" _col(44) "5%" _col(52) "10%"
            di as text "{hline 58}"
            di as text " t(alpha)" _col(12) as result %8.3f `smint' _col(21) %8.3f `cv1' _col(30) %8.3f `cv25' _col(39) %8.3f `cv5' _col(47) %8.3f `cv10'
            di as text "{hline 58}"
        }
        else if `cv1'<. {
            di as text _col(28) "Lumsdaine-Papell"
            di as text _col(12) "Test" _col(22) "------- critical value -------"
            di as text _col(9) "statistic" _col(27) "1%" _col(36) "5%" _col(44) "10%"
            di as text "{hline 50}"
            di as text " t(alpha)" _col(12) as result %8.3f `smint' _col(21) %8.3f `cv1' _col(30) %8.3f `cv5' _col(38) %8.3f `cv10'
            di as text "{hline 50}"
        }
        else {
            di as text "  Test statistic t(alpha) = " as result %9.3f `smint' as text "  (no tabulated critical values)"
        }
        if `cv1'<. {
            local rj1 = cond(`smint'<`cv1', "Reject", "Do not reject")
            di as text " - `rj1' H0 at the 1% significance level"
        }
        if `cv25'<. {
            local rj25 = cond(`smint'<`cv25', "Reject", "Do not reject")
            di as text " - `rj25' H0 at the 2.5% significance level"
        }
        if `cv5'<. {
            local rj5 = cond(`smint'<`cv5', "Reject", "Do not reject")
            di as text " - `rj5' H0 at the 5% significance level"
        }
        if `cv10'<. {
            local rj10 = cond(`smint'<`cv10', "Reject", "Do not reject")
            di as text " - `rj10' H0 at the 10% significance level"
        }
        di as text ""
        if `nbreaks'>0 {
            forvalues i = 1/`nbreaks' {
                local e_i = `bpsmat'[1,`i']
                if `breakn'==1 local comp "intercept"
                if `breakn'==2 local comp "trend"
                if `breakn'==3 local comp "intercept and trend"
                if `breakn'==4 local comp = cond(`i'==1, "intercept and trend", "intercept")
                di as text "  Break `i' (TB`i')" _col(18) ": " as result "`dl_`i''" as text " (obs: " as result `e_i' as text ") - `comp'"
            }
            di as text ""
        }
        di as text "Break regression at the selected break dates:"
        di as text "Dependent variable: " as result "dy = D.`varlist'"
        di as text "Variable" _col(15) %14s "Coefficient" %14s "T-Stat"
        local pos 1
        foreach lab of local vlabs2 {
            local bb = `bmat'[1,`pos']
            local tt = `tmat'[1,`pos']
            di as text `"`lab'"' _col(15) as result %14.4f `bb' %14.4f `tt'
            local pos = `pos' + 1
        }
        di as text ""
        di as text "Replicate the break regression manually with:"
        di as text "  . " as result `"`regcmd'"'
        di as text ""
    }

    * ------------------------------------------------------------------
    * Stored results
    * ------------------------------------------------------------------
    return scalar cdstat = `smint'
    return scalar N = `N'
    return scalar autop = `bestlag'
    return scalar lags = `bestlag'
    return scalar maxlag = `maxlag'
    return scalar pinobs = `pinobs'
    return scalar pi = `pi'
    return scalar signif = `signif'
    return scalar nbreaks = `nbreaks'
    return scalar cv1 = `cv1'
    return scalar cv5 = `cv5'
    return scalar cv10 = `cv10'
    if `cv25'<. return scalar cv25 = `cv25'
    if `tcrit'>0 & `tcrit'<. return scalar tcrit = `tcrit'
    if `nbreaks'>=1 {
        return scalar minent = `tv_1'
        tempname brtimes
        matrix `brtimes' = J(1, `nbreaks', .)
        forvalues i = 1/`nbreaks' {
            matrix `brtimes'[1,`i'] = `tv_`i''
            return scalar break`i' = `tv_`i''
        }
        return matrix breaks = `brtimes'
    }
    if `nbreaks'>=2 return scalar maxent = `tv_2'
    return local mode = cond(`papermode'==1, "paper", "rats")
    return local cvsource "`cvsrc'"
    return local breakdates "`blabels'"
    return local break "`brk'"
    return local method "`mth'"
    return local tsfmt "`tsfmt'"
    return local varname "`varlist'"
    return local regcmd `"`regcmd'"'
    return local cmd "lumpapell"
    return matrix b = `bmat'
    return matrix t = `tmat'
end

version 11.0
mata:

// ---------------------------------------------------------------------
// OLS via quad-precision cross moments (mirrors the moment-matrix OLS
// used by RATS CMOM/LINREG). Returns (K+1) x 2:
// rows 1..K = (b, t), row K+1 = (rss, ndf)
// ---------------------------------------------------------------------
real matrix lumpapell_ols(real matrix X, real colvector y)
{
    real matrix XX, Vi, R
    real colvector Xy, b, se, t
    real scalar yy, rss, n, k, s2

    n = rows(X)
    k = cols(X)
    XX = quadcross(X, X)
    Xy = quadcross(X, y)
    yy = quadcross(y, y)
    Vi = invsym(XX)
    b = Vi * Xy
    rss = yy - quadcross(Xy, b)
    s2 = rss / (n - k)
    se = sqrt(s2 :* diagonal(Vi))
    t = b :/ se
    R = (b, t) \ (rss, n - k)
    return(R)
}

// ---------------------------------------------------------------------
// Break dummy columns for a candidate break vector
// brk: 1 = D only, 2 = DT only, 3 = D and DT per break, 4 = model CA
// (D and DT at the first break, D only at the second)
// ---------------------------------------------------------------------
real matrix lumpapell_dummies(real colvector T, real colvector tbv, real scalar brk, real scalar nbreaks)
{
    real matrix B
    real scalar i

    if (brk == 4) {
        B = ((T :> tbv[1]), (T :- tbv[1]) :* (T :> tbv[1]), (T :> tbv[2]))
        return(B)
    }
    B = J(rows(T), 0, .)
    for (i=1; i<=nbreaks; i++) {
        if (brk != 2) B = (B, T :> tbv[i])
        if (brk != 1) B = (B, (T :- tbv[i]) :* (T :> tbv[i]))
    }
    return(B)
}

// ---------------------------------------------------------------------
// Design bundle for a given lag count k: returns (yv, y1, ones, T, Lg)
// with rows startl+k..endl restricted to usable observations
// ---------------------------------------------------------------------
real matrix lumpapell_design(real scalar k, real colvector dyv, real colvector w, real colvector s, real scalar startl, real scalar endl)
{
    real colvector T, posv, keep
    real matrix Lg
    real scalar j, n

    T = (startl + k)::endl
    posv = T :- startl :+ 2
    keep = s[posv] :& (dyv[posv] :< .) :& (w[posv :- 1] :< .)
    for (j=1; j<=k; j++) keep = keep :& (dyv[posv :- j] :< .)
    T = select(T, keep)
    posv = T :- startl :+ 2
    n = rows(T)
    Lg = J(n, 0, .)
    for (j=1; j<=k; j++) Lg = (Lg, dyv[posv :- j])
    return((dyv[posv], w[posv :- 1], J(n, 1, 1), T, Lg))
}

// ---------------------------------------------------------------------
// Evaluate one break combination in article mode: select the lag for
// this combination (per-combination selection, as implied by the k
// columns of Tables 2-4 of the article), run the selected regression,
// and return a header (t1, ksel) \ (n, K) stacked over (b, t)
// ---------------------------------------------------------------------
real matrix lumpapell_eval(real colvector tbv, real scalar brk, real scalar nbreaks, real scalar lmethod, real scalar kmin, real scalar kmax, real scalar signif, real scalar tcrit, real colvector dyv, real colvector w, real colvector s, real scalar startl, real scalar endl)
{
    real matrix B, X, R, Rsel, Z, Zf
    real scalar k, ksel, tl, ndf, sig, ic, icmin, n0, icmult, nsel

    ksel = kmin
    Rsel = J(0, 2, .)
    nsel = .
    if (lmethod == 5) {
        ksel = 0
        for (k=kmax; k>=1; k--) {
            Z = lumpapell_design(k, dyv, w, s, startl, endl)
            B = lumpapell_dummies(Z[., 4], tbv, brk, nbreaks)
            X = (Z[., 2], B, Z[., 3], Z[., 4], Z[., 5..4+k])
            if (rows(X) <= cols(X)) continue
            R = lumpapell_ols(X, Z[., 1])
            if (hasmissing(R[|1,2 \ cols(X),2|])) continue
            tl = R[cols(X), 2]
            ndf = R[rows(R), 2]
            if (tcrit > 0 & tcrit < .) sig = (abs(tl) >= tcrit)
            else sig = (2 * ttail(ndf, abs(tl)) <= signif)
            if (sig) {
                ksel = k
                Rsel = R
                nsel = rows(X)
                break
            }
        }
        if (ksel == 0) {
            Z = lumpapell_design(0, dyv, w, s, startl, endl)
            B = lumpapell_dummies(Z[., 4], tbv, brk, nbreaks)
            X = (Z[., 2], B, Z[., 3], Z[., 4])
            if (rows(X) <= cols(X)) return(J(2, 2, .))
            Rsel = lumpapell_ols(X, Z[., 1])
            if (hasmissing(Rsel[|1,2 \ cols(X),2|])) return(J(2, 2, .))
            nsel = rows(X)
        }
    }
    else if (lmethod >= 2 & lmethod <= 4) {
        // criterion selection on the fixed sample starting at kmax
        Zf = lumpapell_design(kmax, dyv, w, s, startl, endl)
        n0 = rows(Zf)
        B = lumpapell_dummies(Zf[., 4], tbv, brk, nbreaks)
        if (lmethod == 2) icmult = 2 / n0
        else if (lmethod == 3) icmult = ln(n0) / n0
        else icmult = 2 * ln(ln(n0)) / n0
        X = (Zf[., 2], B, Zf[., 3], Zf[., 4])
        if (rows(X) <= cols(X)) return(J(2, 2, .))
        R = lumpapell_ols(X, Zf[., 1])
        if (hasmissing(R[|1,2 \ cols(X),2|])) return(J(2, 2, .))
        icmin = ln(R[rows(R), 1] / n0) + icmult * cols(X)
        ksel = 0
        for (k=kmax; k>=1; k--) {
            X = (Zf[., 2], B, Zf[., 3], Zf[., 4], Zf[., 5..4+k])
            if (rows(X) <= cols(X)) continue
            R = lumpapell_ols(X, Zf[., 1])
            if (hasmissing(R[|1,2 \ cols(X),2|])) continue
            ic = ln(R[rows(R), 1] / n0) + icmult * cols(X)
            if (ic < icmin) {
                icmin = ic
                ksel = k
            }
        }
        // final regression at the selected lag on its own sample
        Z = lumpapell_design(ksel, dyv, w, s, startl, endl)
        B = lumpapell_dummies(Z[., 4], tbv, brk, nbreaks)
        if (ksel > 0) X = (Z[., 2], B, Z[., 3], Z[., 4], Z[., 5..4+ksel])
        else X = (Z[., 2], B, Z[., 3], Z[., 4])
        if (rows(X) <= cols(X)) return(J(2, 2, .))
        Rsel = lumpapell_ols(X, Z[., 1])
        if (hasmissing(Rsel[|1,2 \ cols(X),2|])) return(J(2, 2, .))
        nsel = rows(X)
    }
    else {
        // fixed lag count
        Z = lumpapell_design(kmin, dyv, w, s, startl, endl)
        B = lumpapell_dummies(Z[., 4], tbv, brk, nbreaks)
        if (kmin > 0) X = (Z[., 2], B, Z[., 3], Z[., 4], Z[., 5..4+kmin])
        else X = (Z[., 2], B, Z[., 3], Z[., 4])
        if (rows(X) <= cols(X)) return(J(2, 2, .))
        Rsel = lumpapell_ols(X, Z[., 1])
        if (hasmissing(Rsel[|1,2 \ cols(X),2|])) return(J(2, 2, .))
        nsel = rows(X)
    }
    if (rows(Rsel) < 2) return(J(2, 2, .))
    return((Rsel[1,2], ksel) \ (nsel, rows(Rsel)-1) \ Rsel[1..rows(Rsel)-1, .])
}

// ---------------------------------------------------------------------
// Main routine: RATS path replicates lpunit.src line by line; the
// article path implements the grid and lag selection of the article
// ---------------------------------------------------------------------
void lumpapell_main(string scalar yname, string scalar tname, string scalar sname, string scalar wname, real scalar tmin, real scalar dlt, real scalar startl, real scalar endl, real scalar brk, real scalar nbreaks, real scalar lmethod, real scalar maxlag0, real scalar mlraw, real scalar signif, real scalar pi, real scalar papermode, real scalar tcrit)
{
    real matrix D, R, X, B, Lg, Lg0, C, H
    real colvector w, s, dyv, y1, T, posv, keep, bps, ub, bestb, bestt, bestbps, yv, tbv
    real scalar M, i, j, p, znobs, maxlag, bestlag, lower, pinobs, n0, icmult, ic, icmin, lag, tst, pv, ndf, mint, done, first, K, kmax, kmin, lowc, upc, gap, tb1, tb2, bestN, bestK, bestk

    // load the window [startl-1, endl], sort by time, map to positions
    D = st_data(., (tname, yname, sname), wname)
    D = sort(D, 1)
    M = endl - startl + 2
    w = J(M, 1, .)
    s = J(M, 1, 0)
    for (i=1; i<=rows(D); i++) {
        p = round((D[i,1] - tmin)/dlt) + 1 - (startl - 2)
        w[p] = D[i,2]
        s[p] = D[i,3]
    }

    // dy over positions (position 1 = entry startl-1 has no dy)
    dyv = J(M, 1, .)
    for (p=2; p<=M; p++) dyv[p] = w[p] - w[p-1]

    // pinobs from the raw MAXLAGS option value (RATS lines 124-129)
    lower = startl + mlraw + 1
    pinobs = trunc(pi * (endl - lower + 1))

    // default maxlag (RATS lines 133-141)
    znobs = endl - startl + 1
    if (maxlag0 >= 0) maxlag = maxlag0
    else maxlag = trunc(znobs^0.25)

    // -----------------------------------------------------------------
    // RATS path lag selection on the base model, fixed sample
    // startl+maxlag..endl (RATS lines 147-177); the article path
    // selects the lag per break combination inside the grid instead
    // -----------------------------------------------------------------
    if (lmethod == 1) bestlag = maxlag
    else if (!papermode) {
        T = (startl + maxlag)::endl
        posv = T :- startl :+ 2
        keep = s[posv] :& (dyv[posv] :< .) :& (w[posv :- 1] :< .)
        for (j=1; j<=maxlag; j++) keep = keep :& (dyv[posv :- j] :< .)
        T = select(T, keep)
        posv = T :- startl :+ 2
        n0 = rows(T)
        if (n0 <= 3 + maxlag) {
            st_numscalar("__lp_rc", 1)
            return
        }
        yv = dyv[posv]
        y1 = w[posv :- 1]
        C = (J(n0, 1, 1), T)
        Lg0 = J(n0, 0, .)
        for (j=1; j<=maxlag; j++) Lg0 = (Lg0, dyv[posv :- j])
        if (lmethod == 2) icmult = 2 / n0
        else if (lmethod == 3) icmult = ln(n0) / n0
        else if (lmethod == 4) icmult = 2 * ln(ln(n0)) / n0
        else icmult = 0
        // lag 0 baseline (RATS lines 159-162)
        R = lumpapell_ols((C, y1), yv)
        ic = ln(R[4,1] / n0) + icmult * 3
        icmin = ic
        bestlag = 0
        // loop from maxlag down to 1 (RATS lines 163-177)
        for (lag=maxlag; lag>=1; lag--) {
            X = (C, y1, Lg0[., 1..lag])
            R = lumpapell_ols(X, yv)
            if (lmethod == 5) {
                tst = R[3 + lag, 2]
                ndf = R[rows(R), 2]
                pv = 2 * ttail(ndf, abs(tst))
                if (pv <= signif) {
                    bestlag = lag
                    break
                }
            }
            else {
                ic = ln(R[rows(R), 1] / n0) + icmult * (3 + lag)
                if (ic < icmin) {
                    icmin = ic
                    bestlag = lag
                }
            }
        }
    }
    else bestlag = -1

    if (!papermode) {
        // -------------------------------------------------------------
        // RATS grid over break combinations (RATS lines 193-246)
        // -------------------------------------------------------------
        lower = startl + bestlag + 1
        if (nbreaks > 0 & pinobs < 1) {
            st_numscalar("__lp_rc", 2)
            return
        }
        if (nbreaks > 0) {
            if ((lower - 1) + pinobs * (nbreaks + 1) >= endl + 1) {
                st_numscalar("__lp_rc", 3)
                return
            }
        }
        bps = J(nbreaks, 1, 0)
        ub = J(nbreaks, 1, 0)
        for (i=1; i<=nbreaks; i++) {
            bps[i] = (lower - 1) + pinobs * i
            ub[i] = endl + 1 - pinobs * (nbreaks + 1 - i)
        }

        // estimation rows for the final regressions
        T = (startl + bestlag)::endl
        posv = T :- startl :+ 2
        keep = s[posv] :& (dyv[posv] :< .) :& (w[posv :- 1] :< .)
        for (j=1; j<=bestlag; j++) keep = keep :& (dyv[posv :- j] :< .)
        T = select(T, keep)
        posv = T :- startl :+ 2
        n0 = rows(T)
        yv = dyv[posv]
        y1 = w[posv :- 1]
        C = (J(n0, 1, 1), T)
        Lg = J(n0, 0, .)
        for (j=1; j<=bestlag; j++) Lg = (Lg, dyv[posv :- j])
        K = 1 + nbreaks * (1 + (brk == 3)) + 2 + bestlag
        if (n0 <= K) {
            st_numscalar("__lp_rc", 1)
            return
        }

        mint = .
        first = 1
        bestb = J(K, 1, .)
        bestt = J(K, 1, .)
        bestbps = J(nbreaks, 1, .)
        done = 0
        while (!done) {
            // dummies grouped by break point: D then DT for each break
            B = lumpapell_dummies(T, bps, brk, nbreaks)
            X = (y1, B, C, Lg)
            R = lumpapell_ols(X, yv)
            tst = R[1, 2]
            // keep the first value, afterwards strict minimum (RATS line 230)
            if (first | tst < mint) {
                mint = tst
                bestb = R[1..K, 1]
                bestt = R[1..K, 2]
                bestbps = bps
                first = 0
            }
            // odometer (RATS lines 235-245)
            done = 1
            for (i=nbreaks; i>=1; i--) {
                bps[i] = bps[i] + 1
                if (bps[i] >= ub[i]) continue
                for (j=i+1; j<=nbreaks; j++) bps[j] = bps[j-1] + pinobs
                done = 0
                break
            }
        }
        bestN = n0
        bestK = K
        bestk = bestlag
    }
    else {
        // -------------------------------------------------------------
        // Article grid: candidates from the second observation of the
        // series to the next-to-last observation, break dates at least
        // two periods apart, model CA unordered; the lag is selected
        // for each break combination
        // -------------------------------------------------------------
        lowc = startl
        upc = endl - 1
        gap = 2
        kmax = (lmethod == 1 ? bestlag : maxlag)
        kmin = (lmethod == 1 ? bestlag : 0)
        if (nbreaks > 0 & lowc + gap * (nbreaks - 1) > upc) {
            st_numscalar("__lp_rc", 3)
            return
        }

        mint = .
        first = 1
        bestb = J(0, 1, .)
        bestt = J(0, 1, .)
        bestbps = J(max((nbreaks,1)), 1, .)
        bestk = .
        bestN = .
        bestK = .
        if (nbreaks == 0) {
            H = lumpapell_eval(J(0,1,.), brk, 0, lmethod, kmin, kmax, signif, tcrit, dyv, w, s, startl, endl)
            mint = H[1,1]
            if (rows(H) > 2) {
                bestk = H[1,2]
                bestN = H[2,1]
                bestK = H[2,2]
                bestb = H[3..rows(H), 1]
                bestt = H[3..rows(H), 2]
            }
        }
        else if (brk == 4) {
            for (tb1=lowc; tb1<=upc; tb1++) {
                for (tb2=lowc; tb2<=upc; tb2++) {
                    if (abs(tb1 - tb2) < gap) continue
                    tbv = (tb1 \ tb2)
                    H = lumpapell_eval(tbv, brk, nbreaks, lmethod, kmin, kmax, signif, tcrit, dyv, w, s, startl, endl)
                    tst = H[1,1]
                    if ((first & tst < .) | tst < mint) {
                        mint = tst
                        bestk = H[1,2]
                        bestN = H[2,1]
                        bestK = H[2,2]
                        bestb = H[3..rows(H), 1]
                        bestt = H[3..rows(H), 2]
                        bestbps = tbv
                        first = 0
                    }
                }
            }
        }
        else {
            bps = J(nbreaks, 1, 0)
            ub = J(nbreaks, 1, 0)
            for (i=1; i<=nbreaks; i++) {
                bps[i] = lowc + gap * (i - 1)
                ub[i] = upc - gap * (nbreaks - i)
            }
            done = 0
            while (!done) {
                H = lumpapell_eval(bps, brk, nbreaks, lmethod, kmin, kmax, signif, tcrit, dyv, w, s, startl, endl)
                tst = H[1,1]
                if ((first & tst < .) | tst < mint) {
                    mint = tst
                    bestk = H[1,2]
                    bestN = H[2,1]
                    bestK = H[2,2]
                    bestb = H[3..rows(H), 1]
                    bestt = H[3..rows(H), 2]
                    bestbps = bps
                    first = 0
                }
                done = 1
                for (i=nbreaks; i>=1; i--) {
                    bps[i] = bps[i] + 1
                    if (bps[i] > ub[i]) continue
                    for (j=i+1; j<=nbreaks; j++) bps[j] = bps[j-1] + gap
                    done = 0
                    break
                }
            }
        }
        if (mint >= .) {
            st_numscalar("__lp_rc", 1)
            return
        }
        bestlag = bestk
    }

    st_numscalar("__lp_rc", 0)
    st_numscalar("__lp_mint", mint)
    st_numscalar("__lp_bestlag", bestlag)
    st_numscalar("__lp_maxlag", maxlag)
    st_numscalar("__lp_pinobs", pinobs)
    st_numscalar("__lp_N", bestN)
    st_numscalar("__lp_K", bestK)
    st_matrix("__lp_b", bestb')
    st_matrix("__lp_t", bestt')
    if (nbreaks > 0) st_matrix("__lp_bps", bestbps')
}

end