*! xtkpybreak version 1.0.0 09jul2026 *! CCE estimation with non-stationary factors (Kapetanios, Pesaran & Yamagata 2011) *! and multiple structural changes in non-stationary heterogeneous panels *! (Baltagi, Feng & Wang 2025). *! Author: Dr Merwan Roudane (merwanroudane920@gmail.com) *! GitHub: https://github.com/merwanroudane *! *! Implements, faithfully to the two papers: *! subcommand cce : CCEMG and CCEP estimators, valid under I(1) common *! factors. Proxy f_t by cross-section averages *! (Pesaran 2006; KPY 2011, eqs 14-20, 38-44). *! subcommand break : least-squares estimation of multiple structural *! breaks (Bai & Perron 1998, 2003 dynamic programming) *! in slopes (model 4) or in slopes AND error-factor *! loadings (model 5) of non-stationary heterogeneous *! panels with common factors (BFW 2025, eqs 4-19). *! *! Data: Stata long format, xtset panelvar timevar, balanced panel required. program define xtkpybreak, eclass version 14.0 if replay() { if ("`e(cmd)'" != "xtkpybreak") error 301 _xtkpb_display exit } gettoken sub 0 : 0, parse(" ,") local l = length("`sub'") if ("`sub'" == substr("cce",1,max(2,`l')) & `l'>=2 & "`sub'"!="") { _xtkpb_cce `0' } else if ("`sub'" == substr("break",1,max(2,`l')) & `l'>=2 & "`sub'"!="") { _xtkpb_break `0' } else { di as error "unknown subcommand '`sub''" di as error "syntax: xtkpybreak {bf:cce} ... or xtkpybreak {bf:break} ..." exit 198 } end *----------------------------------------------------------------------------- * CCE (Kapetanios, Pesaran & Yamagata 2011) *----------------------------------------------------------------------------- program define _xtkpb_cce, eclass syntax varlist(min=2 numeric ts) [if] [in] [ , /// PROXY(string) /// ESTimator(string) /// NOCONstant /// Level(cilevel) /// COEFplot FACTORplot /// name(string) ] gettoken dv xvars : varlist local k : word count `xvars' local nocons = ("`noconstant'"!="") * ---- proxy choice ---- if ("`proxy'"=="") local proxy "yx" if !inlist("`proxy'","yx","x") { di as error "proxy() must be {bf:yx} (default) or {bf:x}" exit 198 } local pyx = ("`proxy'"=="yx") * ---- which estimator is posted in e(b) ---- if ("`estimator'"=="") local estimator "mg" if !inlist("`estimator'","mg","pooled","pool") { di as error "estimator() must be {bf:mg} or {bf:pooled}" exit 198 } if ("`estimator'"=="pool") local estimator "pooled" * ---- panel setup (inline, per skill note on tempvar scope) ---- qui xtset local ivar "`r(panelvar)'" local tvar "`r(timevar)'" if ("`ivar'"=="" | "`tvar'"=="") { di as error "data must be {bf:xtset} panelvar timevar first" exit 459 } marksample touse markout `touse' `dv' `xvars' qui count if `touse' if (r(N)==0) error 2000 * ---- Mata engine ---- tempname bMG seMG bP seP Bi Bise VMG VP mata: xtkpb_cce_engine("`dv'","`xvars'","`ivar'","`tvar'","`touse'", /// `pyx', `nocons', "`bMG'","`seMG'","`bP'","`seP'","`Bi'","`Bise'", /// "`VMG'","`VP'") if (`r(balanced)'==0) { di as error "xtkpybreak requires a balanced panel " /// "(N*T = `r(nobs)' does not match `r(N_g)' x `r(Tbar)')" exit 451 } local N = r(N_g) local T = r(Tbar) * ---- display ---- _xtkpb_cce_tab, dv(`dv') xvars(`xvars') n(`N') t(`T') proxy(`proxy') /// level(`level') bmg(`bMG') semg(`seMG') bp(`bP') sep(`seP') * ---- post e() (chosen estimator) ---- tempname b V if ("`estimator'"=="mg") { matrix `b' = `bMG' matrix `V' = `VMG' } else { matrix `b' = `bP' matrix `V' = `VP' } matrix colnames `b' = `xvars' matrix colnames `V' = `xvars' matrix rownames `V' = `xvars' ereturn post `b' `V', esample(`touse') depname(`dv') obs(`=`N'*`T'') ereturn local cmd "xtkpybreak" ereturn local subcmd "cce" ereturn local estimator "`estimator'" ereturn local proxy "`proxy'" ereturn local depvar "`dv'" ereturn local indepvars "`xvars'" ereturn local ivar "`ivar'" ereturn local tvar "`tvar'" ereturn local title "CCE estimation with non-stationary factors (KPY 2011)" ereturn scalar N_g = `N' ereturn scalar Tbar = `T' ereturn scalar N = `N'*`T' ereturn scalar k = `k' matrix colnames `bMG' = `xvars' matrix colnames `bP' = `xvars' ereturn matrix b_mg = `bMG' ereturn matrix se_mg = `seMG' ereturn matrix b_pooled = `bP' ereturn matrix se_pooled= `seP' ereturn matrix b_i = `Bi' ereturn matrix se_i = `Bise' ereturn matrix V_mg = `VMG' ereturn matrix V_pooled = `VP' * ---- plots (isolated frames; require Stata 16; never abort estimation) ---- if ("`name'"=="") local name "xtkpb" if ("`coefplot'`factorplot'"!="" & c(stata_version)<16) { di as text "(plots require Stata 16 or later; estimation results are unaffected)" } if ("`coefplot'"!="" & c(stata_version)>=16) { capture noisily _xtkpb_cce_coefplot, xvars(`xvars') level(`level') name(`name') if (_rc) di as error "(coefplot skipped: graph error rc=`=_rc')" } if ("`factorplot'"!="" & c(stata_version)>=16) { capture noisily _xtkpb_factorplot, dv(`dv') xvars(`xvars') ivar(`ivar') /// tvar(`tvar') touse(`touse') pyx(`pyx') name(`name') if (_rc) di as error "(factorplot skipped: graph error rc=`=_rc')" } end program define _xtkpb_cce_tab syntax , dv(string) xvars(string) n(string) t(string) proxy(string) /// level(string) bmg(name) semg(name) bp(name) sep(name) local k : word count `xvars' local zc = invnormal(1-(100-`level')/200) di "" di as text "Common Correlated Effects estimation, non-stationary factors" /// _col(58) as text "(KPY 2011)" di as text "{hline 78}" di as text "Dependent variable : " as result "`dv'" di as text "Panels (N) : " as result %-9.0g `n' /// as text " Time periods (T): " as result %-9.0g `t' if ("`proxy'"=="yx") di as text "Factor proxy : " as result /// "cross-section averages of (`dv', regressors)" else di as text "Factor proxy : " as result /// "cross-section averages of regressors" di as text "{hline 78}" di as text "{col 1}Variable{col 16}{col 22}Coef.{col 34}Std.Err." /// "{col 47}z{col 56}P>|z|{col 65}[`level'% C.I.]" foreach blk in MG POOLED { if ("`blk'"=="MG") { di as text "{hline 78}" di as text "CCEMG (mean-group)" local bb `bmg' local ss `semg' } else { di as text "{hline 78}" di as text "CCEP (pooled)" local bb `bp' local ss `sep' } di as text "{hline 78}" forvalues j = 1/`k' { local vn : word `j' of `xvars' local co = `bb'[1,`j'] local se = `ss'[1,`j'] local zv = `co'/`se' local pv = 2*normal(-abs(`zv')) local lo = `co' - `zc'*`se' local hi = `co' + `zc'*`se' local st "" if (`pv'<0.10) local st "*" if (`pv'<0.05) local st "**" if (`pv'<0.01) local st "***" di as text %-15s abbrev("`vn'",14) /// as result "{col 18}" %9.4f `co' "`st'" /// as result "{col 33}" %9.4f `se' /// as result "{col 44}" %7.2f `zv' /// as result "{col 54}" %6.3f `pv' /// as result "{col 62}" %9.4f `lo' /// as result "{col 72}" %9.4f `hi' } } di as text "{hline 78}" di as text "Significance: * 10%, ** 5%, *** 1%. " /// "MG s.e. from cross-panel dispersion (KPY eq.38);" di as text "pooled s.e. from KPY eq.42-44. Individual b_i in e(b_i)." end *----------------------------------------------------------------------------- * BREAK (Baltagi, Feng & Wang 2025) *----------------------------------------------------------------------------- program define _xtkpb_break, eclass syntax varlist(min=2 numeric ts) [if] [in] [ , /// NBReaks(integer 1) /// LOADings /// PROXY(string) /// TRIM(real 0.10) /// NONStationary /// NOCONstant /// HAC(integer -1) /// Level(cilevel) /// BREAKplot COEFEvolution /// name(string) ] gettoken dv xvars : varlist local k : word count `xvars' if (`nbreaks'<1) { di as error "nbreaks() must be >= 1" exit 198 } if ("`proxy'"=="") local proxy "x" if !inlist("`proxy'","yx","x") { di as error "proxy() must be {bf:x} (default, BFW) or {bf:yx}" exit 198 } local pyx = ("`proxy'"=="yx") local ldg = ("`loadings'"!="") local nocons = ("`noconstant'"!="") if (`trim'<=0 | `trim'>=0.5) { di as error "trim() must be in (0,0.5); BFW use 0.10" exit 198 } qui xtset local ivar "`r(panelvar)'" local tvar "`r(timevar)'" if ("`ivar'"=="" | "`tvar'"=="") { di as error "data must be {bf:xtset} panelvar timevar first" exit 459 } marksample touse markout `touse' `dv' `xvars' qui count if `touse' if (r(N)==0) error 2000 tempname MGco MGse Bdate Bidx MGco_p MGse_p ssr Bi Bise mata: xtkpb_break_engine("`dv'","`xvars'","`ivar'","`tvar'","`touse'", /// `nbreaks', `ldg', `pyx', `trim', `hac', `nocons', /// "`MGco'","`MGse'","`Bdate'","`Bidx'","`MGco_p'","`MGse_p'", /// "`ssr'","`Bi'","`Bise'") if (`r(balanced)'==0) { di as error "xtkpybreak requires a balanced panel" exit 451 } if (`r(feasible)'==0) { di as error "trim()/nbreaks() leave too few observations per regime " /// "for T=`r(Tbar)' and " `k' " regressor(s); reduce nbreaks() or trim()" exit 198 } local N = r(N_g) local T = r(Tbar) local hacw = r(hacw) local m = `nbreaks' local nreg = `m'+1 local nsflag = ("`nonstationary'"!="") _xtkpb_break_tab, dv(`dv') xvars(`xvars') n(`N') t(`T') m(`m') /// proxy(`proxy') loadings(`ldg') nonstat(`nsflag') hac(`hacw') /// level(`level') mgco(`MGco') mgse(`MGse') bdate(`Bdate') ssr(`ssr') * ---- build stacked e(b): eq r1..r(m+1) each with xvars ---- tempname b V local names "" forvalues r = 1/`nreg' { foreach v of local xvars { local names "`names' r`r':`v'" } } mata: xtkpb_break_post("`MGco'","`MGse'","`b'","`V'", `nreg', `k') matrix colnames `b' = `names' matrix colnames `V' = `names' matrix rownames `V' = `names' ereturn post `b' `V', esample(`touse') depname(`dv') obs(`=`N'*`T'') ereturn local cmd "xtkpybreak" ereturn local subcmd "break" ereturn local proxy "`proxy'" ereturn local depvar "`dv'" ereturn local indepvars "`xvars'" ereturn local ivar "`ivar'" ereturn local tvar "`tvar'" if (`ldg') ereturn local model "breaks in slopes and factor loadings (BFW model 5)" else ereturn local model "breaks in slopes only (BFW model 4)" if ("`nonstationary'"!="") ereturn local regressors "I(1) - Case 2 (T-consistent)" else ereturn local regressors "I(0) - Case 1 (root-T-consistent)" ereturn local title "Multiple structural breaks, non-stationary panels (BFW 2025)" ereturn scalar N_g = `N' ereturn scalar Tbar = `T' ereturn scalar N = `N'*`T' ereturn scalar k = `k' ereturn scalar nbreaks = `m' ereturn scalar ssr = `ssr'[1,1] ereturn scalar trim = `trim' ereturn scalar hac = `hacw' matrix colnames `Bdate' = break ereturn matrix breakdates = `Bdate' ereturn matrix breakobs = `Bidx' ereturn matrix b_regime = `MGco' ereturn matrix se_regime = `MGse' ereturn matrix bp_regime = `MGco_p' ereturn matrix sep_regime = `MGse_p' ereturn matrix b_i = `Bi' ereturn matrix se_i = `Bise' * ---- plots (isolated frames; require Stata 16; never abort estimation) ---- if ("`name'"=="") local name "xtkpb" if ("`breakplot'`coefevolution'"!="" & c(stata_version)<16) { di as text "(plots require Stata 16 or later; estimation results are unaffected)" } if ("`breakplot'"!="" & c(stata_version)>=16) { capture noisily _xtkpb_breakplot, dv(`dv') ivar(`ivar') tvar(`tvar') /// touse(`touse') name(`name') if (_rc) di as error "(breakplot skipped: graph error rc=`=_rc')" } if ("`coefevolution'"!="" & c(stata_version)>=16) { capture noisily _xtkpb_coefevo, xvars(`xvars') tvar(`tvar') ivar(`ivar') /// touse(`touse') m(`m') level(`level') name(`name') if (_rc) di as error "(coefevolution skipped: graph error rc=`=_rc')" } end program define _xtkpb_break_tab syntax , dv(string) xvars(string) n(string) t(string) m(string) /// proxy(string) loadings(string) nonstat(string) hac(string) /// level(string) mgco(name) mgse(name) bdate(name) ssr(name) local k : word count `xvars' local nreg = `m'+1 local zc = invnormal(1-(100-`level')/200) di "" di as text "Multiple structural breaks in non-stationary heterogeneous panels" /// _col(66) as text "(BFW 2025)" di as text "{hline 78}" di as text "Dependent variable : " as result "`dv'" di as text "Panels (N) : " as result %-9.0g `n' /// as text " Time periods (T): " as result %-9.0g `t' if ("`loadings'"=="1") di as text "Break specification: " as result /// "slopes and error-factor loadings (model 5)" else di as text "Break specification: " as result "slopes only (model 4)" if ("`nonstat'"=="1") di as text "Regressors : " as result /// "I(1) after CCE transform - Case 2 (T-consistent)" else di as text "Regressors : " as result /// "I(0) after CCE transform - Case 1 (root-T-consistent)" if ("`proxy'"=="yx") di as text "Factor proxy : " as result /// "cross-section averages of (`dv', regressors)" else di as text "Factor proxy : " as result /// "cross-section averages of regressors" * ---- estimated break dates ---- di as text "{hline 78}" di as text "Estimated break points (Bai-Perron dynamic programming):" forvalues j = 1/`m' { local bd = `bdate'[`j',1] di as text " break " %2.0f `j' as text " at " as result /// "`tvar'" as text " = " as result %9.0g `bd' } local ss = `ssr'[1,1] di as text " total SSR (all panels) = " as result %12.4f `ss' * ---- regime mean-group coefficients ---- di as text "{hline 78}" di as text "Mean-group slope estimates by regime (95%% level = `level')" di as text "{hline 78}" di as text "{col 1}Regime{col 12}Variable{col 26}Coef.{col 38}Std.Err." /// "{col 51}z{col 60}P>|z|" di as text "{hline 78}" forvalues r = 1/`nreg' { di as text "regime `r'" forvalues j = 1/`k' { local vn : word `j' of `xvars' local co = `mgco'[`r',`j'] local se = `mgse'[`r',`j'] local zv = `co'/`se' local pv = 2*normal(-abs(`zv')) local st "" if (`pv'<0.10) local st "*" if (`pv'<0.05) local st "**" if (`pv'<0.01) local st "***" di as text "{col 12}" %-14s abbrev("`vn'",13) /// as result "{col 24}" %9.4f `co' "`st'" /// as result "{col 37}" %9.4f `se' /// as result "{col 49}" %7.2f `zv' /// as result "{col 58}" %6.3f `pv' } } di as text "{hline 78}" di as text "Significance: * 10%, ** 5%, *** 1%. " /// "MG s.e. from cross-panel dispersion (BFW Prop.2)." di as text "Break points estimated jointly with slopes (BFW eq.13)." di as text "Individual b_i in e(b_i); Newey-West s.e. (window = `hac') " /// "in e(se_i) (BFW Prop.1)." end *----------------------------------------------------------------------------- * replay *----------------------------------------------------------------------------- program define _xtkpb_display if ("`e(subcmd)'"=="cce") { tempname bmg semg bp sep matrix `bmg' = e(b_mg) matrix `semg' = e(se_mg) matrix `bp' = e(b_pooled) matrix `sep' = e(se_pooled) _xtkpb_cce_tab, dv(`e(depvar)') xvars(`e(indepvars)') /// n(`e(N_g)') t(`e(Tbar)') proxy(`e(proxy)') level(95) /// bmg(`bmg') semg(`semg') bp(`bp') sep(`sep') } else { di as text "(replay of xtkpybreak break; full table via stored e() results)" ereturn display } end *----------------------------------------------------------------------------- * Graphs *----------------------------------------------------------------------------- program define _xtkpb_cce_coefplot syntax , xvars(string) level(string) name(string) local vn : word 1 of `xvars' local zc = invnormal(1-(100-`level')/200) tempname M bmg smg matrix `M' = e(b_i) matrix `bmg' = e(b_mg) matrix `smg' = e(se_mg) local mg = `bmg'[1,1] local se = `smg'[1,1] local lo = `mg' - `zc'*`se' local hi = `mg' + `zc'*`se' local NN = colsof(`M') forvalues i = 1/`NN' { local v`i' = `M'[1,`i'] } * build the graph in an isolated frame (never touches the caller's data) tempname fr frame create `fr' frame `fr' { qui set obs `NN' qui gen double bcoef = . forvalues i = 1/`NN' { qui replace bcoef = `v`i'' in `i' } qui gen long panel = _n twoway (scatter bcoef panel, msymbol(oh) mcolor(navy)) /// , yline(`mg', lcolor(cranberry) lwidth(medthick)) /// yline(`lo', lpattern(dash) lcolor(gs8)) /// yline(`hi', lpattern(dash) lcolor(gs8)) /// title("Heterogeneous CCE slopes: `vn'") /// subtitle("per-panel b{sub:i} with CCEMG (line) and `level'% band") /// ytitle("b{sub:i}(`vn')") xtitle("panel index") /// graphregion(color(white)) legend(off) /// name(`name'_coef, replace) } frame drop `fr' end program define _xtkpb_factorplot syntax , dv(string) xvars(string) ivar(string) tvar(string) /// touse(string) pyx(string) name(string) tempname fr frame put `dv' `xvars' `tvar' if e(sample), into(`fr') frame `fr' { tempvar ybar qui bysort `tvar': egen double `ybar' = mean(`dv') local avlist "`ybar'" local anames `"1 "CS avg `dv'""' local pos 2 foreach v of local xvars { tempvar xb`pos' qui bysort `tvar': egen double `xb`pos'' = mean(`v') local avlist "`avlist' `xb`pos''" local anames `"`anames' `pos' "CS avg `v'""' local ++pos } qui bysort `tvar': keep if _n==1 local plots "" foreach a of local avlist { local plots "`plots' (line `a' `tvar', lwidth(medthick))" } twoway `plots', title("Cross-section-average factor proxies") /// subtitle("proxies for the unobserved I(1) common factors") /// ytitle("cross-section average") xtitle("`tvar'") /// graphregion(color(white)) /// legend(order(`anames') size(small) rows(1)) /// name(`name'_factor, replace) } frame drop `fr' end program define _xtkpb_breakplot syntax , dv(string) ivar(string) tvar(string) touse(string) name(string) tempname BD matrix `BD' = e(breakdates) local m = rowsof(`BD') local xl "" forvalues j = 1/`m' { local bd = `BD'[`j',1] local xl "`xl' xline(`bd', lcolor(cranberry) lpattern(dash) lwidth(medthick))" } tempname fr frame put `dv' `tvar' if e(sample), into(`fr') frame `fr' { tempvar ybar qui bysort `tvar': egen double `ybar' = mean(`dv') qui bysort `tvar': keep if _n==1 twoway (line `ybar' `tvar', lwidth(medthick) lcolor(navy)), /// `xl' /// title("Cross-section average of `dv' and estimated breaks") /// subtitle("dashed lines = Bai-Perron break dates") /// ytitle("CS average `dv'") xtitle("`tvar'") /// graphregion(color(white)) legend(off) /// name(`name'_break, replace) } frame drop `fr' end program define _xtkpb_coefevo syntax , xvars(string) tvar(string) ivar(string) touse(string) /// m(string) level(string) name(string) tempname MGco MGse Bdate matrix `MGco' = e(b_regime) matrix `MGse' = e(se_regime) matrix `Bdate' = e(breakdates) local vn : word 1 of `xvars' local nreg = `m'+1 local zc = invnormal(1-(100-`level')/200) qui su `tvar' if e(sample), meanonly local tmin = r(min) local tmax = r(max) * pull all matrix values into locals forvalues r = 1/`nreg' { local co`r' = `MGco'[`r',1] local se`r' = `MGse'[`r',1] } local xl "" forvalues j = 1/`m' { local bk`j' = `Bdate'[`j',1] local xl "`xl' xline(`bk`j'', lcolor(gs9) lpattern(dash))" } * build the step series in an isolated frame tempname fr frame create `fr' frame `fr' { qui set obs `nreg' qui gen int regime = _n qui gen double b = . qui gen double lo = . qui gen double hi = . qui gen double tstart = . qui gen double tstop = . forvalues r = 1/`nreg' { qui replace b = `co`r'' in `r' qui replace lo = `co`r''-`zc'*`se`r'' in `r' qui replace hi = `co`r''+`zc'*`se`r'' in `r' } local prev = `tmin' forvalues r = 1/`nreg' { qui replace tstart = `prev' in `r' if (`r'<=`m') { qui replace tstop = `bk`r'' in `r' local prev = `bk`r'' } else qui replace tstop = `tmax' in `r' } qui expand 2 sort regime by regime: gen double tt = cond(_n==1, tstart, tstop) twoway (rarea lo hi tt, sort color(navy%20)) /// (line b tt, sort lwidth(medthick) lcolor(navy)), /// `xl' /// title("Regime evolution of MG slope: `vn'") /// subtitle("step = CCEMG slope by regime, band = `level'% CI") /// ytitle("MG slope(`vn')") xtitle("`tvar'") /// graphregion(color(white)) /// legend(order(2 "MG slope" 1 "`level'% CI") rows(1) size(small)) /// name(`name'_evo, replace) } frame drop `fr' end *============================================================================= * Mata engine (compiles at load; use // comments only, one stmt per line) *============================================================================= version 14.0 mata: mata set matastrict off // ---- reshape long balanced panel into T x N and per-panel X ---------------- // returns success in ok (1/0); fills Ymat (T x N) and pointer array Xp real scalar xtkpb_reshape(string scalar yv, string scalar xv, /// string scalar idv, string scalar tv, string scalar tousev, /// real matrix Ymat, pointer(real matrix) rowvector Xp, /// real scalar N, real scalar T, real scalar k) { real colvector y, id, tt, ids, times, perm real matrix X real scalar i, nobs, a, b y = st_data(., yv, tousev) X = st_data(., xv, tousev) id = st_data(., idv, tousev) tt = st_data(., tv, tousev) k = cols(X) nobs = rows(y) // robust: sort by (panel, time) so panel blocks are contiguous/ordered perm = order((id, tt), (1, 2)) y = y[perm] X = X[perm, .] id = id[perm] tt = tt[perm] ids = uniqrows(id) times = uniqrows(tt) N = rows(ids) T = rows(times) if (nobs != N*T) { return(0) } Ymat = J(T, N, .) Xp = J(1, N, NULL) for (i=1; i<=N; i++) { a = (i-1)*T + 1 b = i*T Ymat[., i] = y[a::b] Xp[i] = &(X[a::b, .]) } return(1) } // ---- CCE engine (KPY 2011) ------------------------------------------------- void xtkpb_cce_engine(string scalar yv, string scalar xv, string scalar idv, /// string scalar tv, string scalar tousev, real scalar pyx, /// real scalar nocons, /// string scalar bMGn, string scalar seMGn, string scalar bPn, /// string scalar sePn, string scalar Bin, string scalar Bisen, /// string scalar VMGn, string scalar VPn) { real matrix Ymat, Hbar, Mbar, Xi, XM, A, B, Ap_i, Sig, Psi, R, Q real matrix SigMG, SigP, VMG, VP, Bse real colvector ybar, yi, bi, c, bMG, bP, d, seMG, seP, ehat, resid real scalar N, T, k, i, ok, dof, s2 pointer(real matrix) rowvector Xp pointer(real matrix) rowvector Ap real matrix xbar, cons, sumXMX, sumXMy N = 0 T = 0 k = 0 ok = xtkpb_reshape(yv, xv, idv, tv, tousev, Ymat, Xp, N, T, k) st_numscalar("r(N_g)", N) st_numscalar("r(Tbar)", T) st_numscalar("r(nobs)", N*T) if (ok==0) { st_numscalar("r(balanced)", 0) return } st_numscalar("r(balanced)", 1) // cross-section averages ybar = J(T, 1, 0) xbar = J(T, k, 0) for (i=1; i<=N; i++) { ybar = ybar + Ymat[., i] xbar = xbar + *Xp[i] } ybar = ybar :/ N xbar = xbar :/ N cons = J(T, 1, 1) if (nocons==1) { if (pyx==1) { Hbar = ybar, xbar } else { Hbar = xbar } } else { if (pyx==1) { Hbar = cons, ybar, xbar } else { Hbar = cons, xbar } } Mbar = I(T) - Hbar * invsym(Hbar'Hbar) * Hbar' // per-panel CCE B = J(k, N, .) Bse = J(k, N, .) Ap = J(1, N, NULL) sumXMX = J(k, k, 0) sumXMy = J(k, 1, 0) dof = T - cols(Hbar) - k if (dof < 1) { dof = 1 } for (i=1; i<=N; i++) { Xi = *Xp[i] yi = Ymat[., i] XM = Xi' * Mbar A = XM * Xi c = XM * yi bi = invsym(A) * c B[., i] = bi // store a FRESH copy of A (do not alias the reused loop var) Ap[i] = &(XM*Xi) sumXMX = sumXMX + A sumXMy = sumXMy + c // individual s.e. (KPY Thm 3, eq.49-50) resid = yi - Xi*bi s2 = (resid' * Mbar * resid) / dof Bse[., i] = sqrt(diagonal(s2 * invsym(A))) } // CCEMG and its variance (eq.14, 38) bMG = J(k, 1, 0) for (i=1; i<=N; i++) { bMG = bMG + B[., i] } bMG = bMG :/ N Sig = J(k, k, 0) for (i=1; i<=N; i++) { d = B[., i] - bMG Sig = Sig + d*d' } SigMG = Sig :/ (N-1) VMG = SigMG :/ N seMG = sqrt(diagonal(VMG)) // CCEP and its variance (eq.20, 42-44) bP = invsym(sumXMX) * sumXMy Psi = J(k, k, 0) for (i=1; i<=N; i++) { Psi = Psi + (*Ap[i]) :/ T } Psi = Psi :/ N R = J(k, k, 0) for (i=1; i<=N; i++) { Q = (*Ap[i]) :/ T d = B[., i] - bMG R = R + Q*d*d'*Q } R = R :/ (N-1) SigP = invsym(Psi) * R * invsym(Psi) VP = SigP :/ N seP = sqrt(diagonal(VP)) st_matrix(bMGn, bMG') st_matrix(seMGn, seMG') st_matrix(bPn, bP') st_matrix(sePn, seP') st_matrix(Bin, B) st_matrix(Bisen, Bse) st_matrix(VMGn, VMG) st_matrix(VPn, VP) } // ---- segment SSR summed across panels (direct accumulate) ------------------ void xtkpb_costmat(pointer(real matrix) rowvector Wp, real matrix YY, /// real scalar N, real scalar T, real scalar q, real scalar h, /// real matrix C) { real matrix Wi, XtX real colvector Xty, yi real scalar i, a, b, cnt, yty, ssr real rowvector z C = J(T, T, 0) for (i=1; i<=N; i++) { Wi = *Wp[i] yi = YY[., i] for (a=1; a<=T; a++) { XtX = J(q, q, 0) Xty = J(q, 1, 0) yty = 0 cnt = 0 for (b=a; b<=T; b++) { z = Wi[b, .] XtX = XtX + z'z Xty = Xty + z'*yi[b] yty = yty + yi[b]*yi[b] cnt = cnt + 1 if (cnt>=h) { ssr = yty - Xty' * invsym(XtX) * Xty C[a, b] = C[a, b] + ssr } } } } } // ---- Bai-Perron dynamic programming --------------------------------------- real rowvector xtkpb_dp(real matrix C, real scalar T, real scalar m, real scalar h) { real matrix OPT, POS real scalar bigM, s, j, lo, hi, l, best, bestl, cand, jj real rowvector brk bigM = 1e300 OPT = J(T, m+1, bigM) POS = J(T, m+1, 0) for (j=1; j<=T; j++) { if (j>=h) { OPT[j, 1] = C[1, j] } } for (s=2; s<=m+1; s++) { for (j=1; j<=T; j++) { lo = (s-1)*h hi = j - h if (hi>=lo) { best = bigM bestl = 0 for (l=lo; l<=hi; l++) { if (OPT[l, s-1] < bigM) { cand = OPT[l, s-1] + C[l+1, j] if (cand < best) { best = cand bestl = l } } } OPT[j, s] = best POS[j, s] = bestl } } } brk = J(1, m, 0) jj = T for (s=m+1; s>=2; s=s-1) { l = POS[jj, s] brk[s-1] = l jj = l } return(brk) } // ---- Newey-West HAC sandwich for one panel (BFW Prop.1, eq.17-18) ---------- // b_i is OLS of (M Y_i) on Xt = M X_i(K0); Var(b_i) = A^-1 S A^-1 with A=Xt'Xt // and S the Bartlett-weighted long-run variance of the scores Xt_t * e_t. real matrix xtkpb_nwhac(real matrix Xt, real colvector e, real matrix Ainv, real scalar w) { real scalar T, p, jj, t real matrix S, Lam real rowvector zt T = rows(Xt) p = cols(Xt) S = J(p, p, 0) for (t=1; t<=T; t++) { zt = Xt[t, .] S = S + (e[t]*e[t]) * (zt' * zt) } for (jj=1; jj<=w; jj++) { Lam = J(p, p, 0) for (t=jj+1; t<=T; t++) { Lam = Lam + (e[t]*e[t-jj]) * (Xt[t,.]' * Xt[t-jj,.]) } S = S + (1 - jj/(w+1)) * (Lam + Lam') } return(Ainv * S * Ainv) } // ---- BREAK engine (BFW 2025) ---------------------------------------------- void xtkpb_break_engine(string scalar yv, string scalar xv, string scalar idv, /// string scalar tv, string scalar tousev, real scalar m, /// real scalar ldg, real scalar pyx, real scalar trim, real scalar hacw, /// real scalar nocons, /// string scalar MGcon, string scalar MGsen, string scalar Bdaten, /// string scalar Bidxn, string scalar MGcopn, string scalar MGsepn, /// string scalar ssrn, string scalar Bin, string scalar Bisen) { real matrix Ymat, xbar, ybar, cons, Hprox, Xi, Wi, Wseg real matrix C, MGco, MGse, MGcop, MGsep, Bslp real matrix Xblk, K1blk, MK1, Abig, Breg, Psi, R, SigP, VP, Q, sumA real matrix Bse_i, Xtil, Vi real colvector times, yi, coef, bounds, bd, bidx, dd, bmgr, cbig, bhat real colvector bMGfull, bP, sumc, dfull, seP, ei real scalar N, T, k, i, ok, q, h, r, nreg, a, b, j, feas, p2, w pointer(real matrix) rowvector Xp, Wp, Ap real rowvector brk real scalar totssr N = 0 T = 0 k = 0 ok = xtkpb_reshape(yv, xv, idv, tv, tousev, Ymat, Xp, N, T, k) st_numscalar("r(N_g)", N) st_numscalar("r(Tbar)", T) st_numscalar("r(nobs)", N*T) if (ok==0) { st_numscalar("r(balanced)", 0) st_numscalar("r(feasible)", 0) return } st_numscalar("r(balanced)", 1) times = uniqrows(st_data(., tv, tousev)) // cross-section averages ybar = J(T, 1, 0) xbar = J(T, k, 0) for (i=1; i<=N; i++) { ybar = ybar + Ymat[., i] xbar = xbar + *Xp[i] } ybar = ybar :/ N xbar = xbar :/ N cons = J(T, 1, 1) // proxy block (BFW eq.9 uses x-bar only; the constant is the KPY D term // and is included unless noconstant is specified) if (nocons==1) { if (pyx==1) { Hprox = ybar, xbar } else { Hprox = xbar } } else { if (pyx==1) { Hprox = cons, ybar, xbar } else { Hprox = cons, xbar } } // ---- break-DATE estimation: fully-broken z=(x,proxy), BFW eq.(11)-(13) ---- // Unrestricted segment-additive SSR (per-regime OLS of y on [x,proxy]); // the DP objective is identical whether or not loadings break. p2 = cols(Hprox) Wp = J(1, N, NULL) for (i=1; i<=N; i++) { Xi = *Xp[i] Wp[i] = &(Xi, Hprox) } q = k + p2 // minimum segment length: each regime needs > (k+p2) free parameters h = floor(trim*T) if (h < q+1) { h = q + 1 } feas = 1 if ((m+1)*h > T) { feas = 0 } st_numscalar("r(feasible)", feas) if (feas==0) { return } // Newey-West window (auto = floor(4 (T/100)^(2/9)) if hacw < 0) if (hacw < 0) { w = floor(4 * (T/100)^(2/9)) } else { w = hacw } if (w < 0) { w = 0 } if (w > T-1) { w = T - 1 } st_numscalar("r(hacw)", w) xtkpb_costmat(Wp, Ymat, N, T, q, h, C) brk = xtkpb_dp(C, T, m, h) nreg = m + 1 bounds = J(nreg+1, 1, 0) bounds[1] = 0 for (j=1; j<=m; j++) { bounds[j+1] = brk[j] } bounds[nreg+1] = T // total SSR of the fully-broken z model at the optimum (BFW eq.13 objective) totssr = 0 for (i=1; i<=N; i++) { Wi = *Wp[i] for (r=1; r<=nreg; r++) { a = bounds[r] + 1 b = bounds[r+1] Wseg = Wi[a::b, .] yi = Ymat[a::b, i] coef = invsym(Wseg'Wseg) * Wseg'yi totssr = totssr + (yi - Wseg*coef)' * (yi - Wseg*coef) } } // ---- loading (proxy) projection M_X(K1), common across panels ---- // model 4 (loadings constant): X(K1) = global proxy -> BFW eq.(14) // model 5 (loadings break) : X(K1) = block-diagonal -> BFW eq.(15) if (ldg==0) { K1blk = Hprox } else { K1blk = J(T, p2*nreg, 0) for (r=1; r<=nreg; r++) { a = bounds[r] + 1 b = bounds[r+1] K1blk[a::b, ((r-1)*p2+1)::(r*p2)] = Hprox[a::b, .] } } MK1 = I(T) - K1blk * invsym(K1blk'K1blk) * K1blk' // ---- SLOPE estimation: partitioned CCE regression, BFW eq.(14)-(15) ---- // bhat_i = [X_i(K0)' M_X(K1) X_i(K0)]^-1 X_i(K0)' M_X(K1) Y_i // X_i(K0) is X_i made block-diagonal across the m+1 regimes. Breg = J(nreg*k, N, .) Bse_i = J(nreg*k, N, .) Ap = J(1, N, NULL) for (i=1; i<=N; i++) { Xi = *Xp[i] yi = Ymat[., i] Xblk = J(T, k*nreg, 0) for (r=1; r<=nreg; r++) { a = bounds[r] + 1 b = bounds[r+1] Xblk[a::b, ((r-1)*k+1)::(r*k)] = Xi[a::b, .] } Abig = Xblk' * MK1 * Xblk cbig = Xblk' * MK1 * yi bhat = invsym(Abig) * cbig Breg[., i] = bhat Ap[i] = &(Xblk' * MK1 * Xblk) // individual Newey-West s.e. (BFW Prop.1, eq.17-18) Xtil = MK1 * Xblk ei = MK1*yi - Xtil*bhat Vi = xtkpb_nwhac(Xtil, ei, invsym(Abig), w) Bse_i[., i] = sqrt(diagonal(Vi)) } // full-length mean group (stack of all regime slopes) bMGfull = J(nreg*k, 1, 0) for (i=1; i<=N; i++) { bMGfull = bMGfull + Breg[., i] } bMGfull = bMGfull :/ N // pooled point estimate (BFW footnote 13; eq.19) sumA = J(nreg*k, nreg*k, 0) sumc = J(nreg*k, 1, 0) for (i=1; i<=N; i++) { sumA = sumA + *Ap[i] sumc = sumc + (*Ap[i]) * Breg[., i] } bP = invsym(sumA) * sumc // pooled variance (sandwich, analogous to CCEP eq.42-44) Psi = J(nreg*k, nreg*k, 0) for (i=1; i<=N; i++) { Psi = Psi + (*Ap[i]) :/ T } Psi = Psi :/ N R = J(nreg*k, nreg*k, 0) for (i=1; i<=N; i++) { Q = (*Ap[i]) :/ T dfull = Breg[., i] - bMGfull R = R + Q*dfull*dfull'*Q } R = R :/ (N-1) SigP = invsym(Psi) * R * invsym(Psi) VP = SigP :/ N seP = sqrt(diagonal(VP)) // ---- regime tables: mean group (Prop.2) and pooled (footnote 13) ---- MGco = J(nreg, k, .) MGse = J(nreg, k, .) MGcop = J(nreg, k, .) MGsep = J(nreg, k, .) for (r=1; r<=nreg; r++) { Bslp = Breg[((r-1)*k+1)::(r*k), .] bmgr = J(k, 1, 0) for (i=1; i<=N; i++) { bmgr = bmgr + Bslp[., i] } bmgr = bmgr :/ N dd = J(k, 1, 0) for (i=1; i<=N; i++) { dd = dd + (Bslp[., i]-bmgr):^2 } MGco[r, .] = bmgr' MGse[r, .] = sqrt(dd :/ (N*(N-1)))' MGcop[r, .] = (bP[((r-1)*k+1)::(r*k)])' MGsep[r, .] = (seP[((r-1)*k+1)::(r*k)])' } bidx = J(m, 1, .) bd = J(m, 1, .) for (j=1; j<=m; j++) { bidx[j] = brk[j] bd[j] = times[brk[j]] } st_matrix(MGcon, MGco) st_matrix(MGsen, MGse) st_matrix(MGcopn, MGcop) st_matrix(MGsepn, MGsep) st_matrix(Bdaten, bd) st_matrix(Bidxn, bidx) st_matrix(ssrn, totssr) st_matrix(Bin, Breg) st_matrix(Bisen, Bse_i) } // ---- stack regime coefficients into b (1 x p) and V (p x p block) --------- void xtkpb_break_post(string scalar MGcon, string scalar MGsen, /// string scalar bn, string scalar Vn, real scalar nreg, real scalar k) { real matrix MGco, MGse, V real rowvector b real scalar r, j, p, pos MGco = st_matrix(MGcon) MGse = st_matrix(MGsen) p = nreg*k b = J(1, p, 0) V = J(p, p, 0) pos = 0 for (r=1; r<=nreg; r++) { for (j=1; j<=k; j++) { pos = pos + 1 b[pos] = MGco[r, j] V[pos, pos] = MGse[r, j]^2 } } st_matrix(bn, b) st_matrix(Vn, V) } end