*! version 1.0.0 29mar2026
help xtbreakmodel version 1.0.0

Title

xtbreakmodel -- Heterogeneous Structural Breaks in Panel Data Models

Version

Version 1.0.0, 29 March 2026

Author: Dr Merwan Roudane (merwanroudane920@gmail.com)

Implements four state-of-the-art panel structural break estimators: GAGFL
(Okui and Wang, 2021), AGFL (Qian and Su, 2016), BFK (Baltagi, Feng and Kao,
2016), and SaRa (Li, Xiao and Chen, 2025).

Syntax

xtbreakmodel depvar indepvars [if] [in]
, method(string) [options]

Model specification
method(string) estimation method: gagfl, pls, bfk, or sara (required)

groups(#) number of latent groups; required for method(gagfl), minimum 2

GAGFL/PLS tuning
maxlambda(#) upper bound of lambda grid; default is maxlambda(100)

minlambda(#) lower bound of lambda grid; default is minlambda(0.01)

ngrid(#) number of log-spaced grid points; default is ngrid(40)

nsim(#) number of random starts for GFE initialization; default is nsim(50)

maxiter(#) maximum GAGFL iterations; default is maxiter(20)

tolerance(#) convergence tolerance; default is tolerance(1e-4)

SaRa tuning
bandwidths(numlist) bandwidth values for local estimation; default is bandwidths(3 5 8)

c1(#) IC penalty constant for static SaRa; default is c1(0.1)

c2(#) IC penalty constant for dynamic SaRa; default is c2(0.025)

Reporting
nograph suppress coefficient path graphs

level(#) confidence level for CI bands; default is level(95)

You must xtset your data before using xtbreakmodel; see xtset.

The panel must be strongly balanced (no gaps or missing values).

All four methods support one or more regressors (K >= 1).

Description

xtbreakmodel estimates structural breaks in panel data regression models
of the form:

y_it = x'_it * beta(t) + alpha_i + epsilon_it, i = 1,...,N, t = 1,...,T

where alpha_i are individual fixed effects and beta(t) is a K-dimensional
coefficient vector that is piecewise constant over time with unknown break
dates. The command offers four different estimation methods, each with
different assumptions about the heterogeneity structure:

Method 1: GAGFL — Grouped Adaptive Group Fused Lasso

Okui and Wang (2021, Journal of Econometrics). The most general method.
Allows units to belong to G latent groups, where each group has its own
break dates and regime-specific coefficients. Simultaneously estimates group
membership and group-specific structural breaks via an iterative algorithm
that alternates between Adaptive Group Fused Lasso (AGFL) estimation within
groups and unit reassignment across groups.

Method 2: PLS — Adaptive Group Fused Lasso (Common Breaks)

Qian and Su (2016, Journal of Econometrics). Assumes all units share the
same break dates but coefficients may change at those dates. All K coefficients
are assumed to break simultaneously (vector break). The method uses Block
Coordinate Descent (BCD) with adaptive fused L1 penalty to determine breaks
and select the optimal tuning parameter via BIC-type information criterion.

Method 3: BFK — Sequential Least Squares

Baltagi, Feng and Kao (2016). Detects common breaks via sequential least
squares minimization. For each candidate break date, the pooled SSR is
computed; the date minimizing total SSR is selected. Breaks are detected
sequentially: after finding the first break, the algorithm searches the
larger remaining segment for the next break.

Method 4: SaRa — Screening and Ranking Algorithm

Li, Xiao and Chen (2025). Uses local kernel estimation of left and right
regression coefficients at each time point. The local statistic measuring
the difference between left and right estimates is computed at multiple
bandwidths. Break candidates are identified as local maximizers of this
statistic, filtered by a median threshold, and the final number is selected
via BIC-type information criterion.

Assumptions

Common assumptions (all methods):

(A1) Balanced panel: The panel must be strongly balanced with no gaps
or missing observations within units.

(A2) Fixed effects: Individual-specific intercepts alpha_i are
eliminated via within-group demeaning or first-period deviation.

(A3) Exogeneity: E(epsilon_it | x_it, alpha_i) = 0 for static
panels. For dynamic specifications, appropriate instruments are required.

(A4) Regularity: The regressors x_it must have sufficient within-period
variation: E(x_it * x'_it) is positive definite for each t.

(A5) Minimum segment length: Each regime must contain at least
one time period. In practice, T >= 3 is required.

Method-specific assumptions:

GAGFL: (i) The number of latent groups G is finite and known or
estimated separately (e.g., via BIC over a range of G). (ii) Group
membership is time-invariant: unit i belongs to group g_i for all t.
(iii) N -> infinity asymptotics with T fixed.

PLS: (i) All N units share the same break dates
(common break assumption). (ii) Under Qian and Su (2016), the initial
estimator is sqrt(N)-consistent for each t, which requires N to be
large relative to T. (iii) The adaptive weight exponent kappa = 2 provides
oracle efficiency.

BFK: (i) Common breaks across units. (ii) The sequential
procedure assumes the true number of breaks is bounded. (iii) The method
is robust to heteroskedasticity in the cross-section via pooled SSR
minimization.

SaRa: (i) Common breaks. (ii) The bandwidths must be chosen
such that h << T^(1/2). (iii) The local kernel estimator requires sufficient
observations in each window (Nh > K). (iv) The method is non-parametric and
robust to distributional assumptions on the errors.

Requirements

Software: Stata 14.0 or higher. No external packages required.

Data:

(i) Minimum dimensions: N >= 10, T >= 3.

(ii) Panel structure: Must be xtset with both panel and time variables.

(iii) No gaps: Strongly balanced panel; use tsfill if needed.

(iv) Non-collinearity: Regressors must not be perfectly collinear.

(v) GAGFL: For G groups, N must be sufficiently large that each group
has at least 2 units. Rule of thumb: N >= 5*G.

(vi) SaRa: T must be large enough to support local estimation: T >= 2*max(bandwidth) + 1.

Options

Model specification

method(string) specifies the estimation method. One of:

gagfl — Grouped AGFL: heterogeneous breaks + latent groups (Okui and Wang, 2021)

pls — AGFL: common breaks across all units (Qian and Su, 2016)

bfk — Sequential least squares (Baltagi, Feng and Kao, 2016)

sara — Screening and ranking algorithm (Li, Xiao and Chen, 2025)

groups(#) specifies the number of latent groups G. Required when
method(gagfl) is specified. Must be >= 2. The choice of G can be guided
by the BIC criterion: estimate for G = 2, 3, 4, ... and choose the G that
minimizes BIC = SSR/(NT) + sigma^2 * (np_G + N) * ln(NT) / (NT).

GAGFL/PLS tuning

maxlambda(#) upper bound of the penalty parameter grid. Default is 100.
Increase if the algorithm selects the boundary lambda.

minlambda(#) lower bound of the penalty parameter grid. Default is 0.01.
Decrease for finer break detection (risks over-segmentation).

ngrid(#) number of log-spaced grid points between minlambda and maxlambda.
Default is 40. Higher values improve lambda selection precision but increase
computation time proportionally.

nsim(#) number of random initializations for the GFE step. Default
is 50. Higher values reduce sensitivity to initial conditions but increase
computation. The MATLAB reference uses 100.

maxiter(#) maximum number of GAGFL outer iterations (alternating between
group-specific AGFL and group reassignment). Default is 20. The algorithm
terminates earlier if the group assignment stabilizes.

tolerance(#) convergence threshold for the BCD solver. Default is 1e-4.
Smaller values give more precise estimates but take longer.

SaRa tuning

bandwidths(numlist) specifies the bandwidth(s) for local estimation.
Default is bandwidths(3 5 8). Larger bandwidths smooth more and detect
only large breaks; smaller bandwidths detect finer breaks but are noisier.
Multiple bandwidths are recommended for robustness: candidates from all
bandwidths are pooled and filtered.

c1(#) IC penalty constant for break number selection in static SaRa.
Default is 0.1. Larger values are more conservative (fewer breaks).

Reporting

nograph suppresses the coefficient path graphs. By default,
xtbreakmodel produces regime-specific coefficient plots with confidence
bands and vertical dashed lines at estimated break dates.

level(#) specifies the confidence level for confidence intervals in
tables and graphs. Default is level(95).

Methodology

1. Data Generating Process

Consider the panel model:

y_it = alpha_i + x'_it * beta_(t) + epsilon_it

where g_i in {1,...,G} denotes unit i's group. Within group g,
the coefficient beta_g(t) is piecewise constant:

beta_g(t) = a_ for tau_ < t <= tau_, j = 1,...,m_g+1

with m_g breaks at dates tau_ < ... < tau_.

When G = 1, this reduces to the common-break model of Qian and Su (2016).

2. GAGFL Algorithm (Okui and Wang, 2021)

The GAGFL estimator consists of three stages:

Stage 1 — Grouped Fixed Effects (GFE) initialization: Following
Bonhomme and Manresa (2015), units are assigned to G groups via iterative
k-means clustering on time-varying coefficients. For each time t and unit i,
an initial beta_t is estimated by cross-sectional OLS. Multiple random starts
(nsim) are used to avoid local minima.

Stage 2 — Group-specific AGFL: For each group g, the Adaptive
Group Fused Lasso is applied to the group's pooled data. This uses Block
Coordinate Descent (BCD) to solve:

minimize (1/Ng) * SUM_i SUM_t [y_it - x'_it * beta_t]^2
+ lambda * SUM_^T w_t * ||beta_t - beta_||

where w_t = ||beta_dot_t - beta_dot_||^(-2) are adaptive weights
computed from initial estimates, and ||.|| is the L2 norm (group lasso).
The penalty lambda is selected by minimizing a BIC-type IC over a grid.

Stage 3 — Group reassignment: Units are reassigned to the group
whose estimated coefficient path minimizes the individual's residual sum of
squares. Stages 2-3 alternate until convergence.

3. BCD Inner Solver (mirrors plsbcd.m / cplsbcd.c)

The BCD algorithm iteratively updates beta_t for t = 1,...,T. For t >= 2, the
update uses soft-thresholding combined with scalar optimization:

d_t = argmin (1/2) * d' * A_t * d + g'_t * d + lambda * w_ * ||d||

where A_t = SUM_ X'_s X_s / N and g_t combines the gradient. If
||g_t|| <= lambda * w_, the solution is d_t = 0 (fusing with previous
period). Otherwise, a line search with Brent's method determines the optimal
step size gamma*, giving d_t = -gamma* * (gamma* A_t + (lambda w_)^2/2 * I)^(-1) * g_t.

4. BFK Sequential Detection (Baltagi, Feng and Kao, 2016)

For each candidate break date k in {1,...,T-1}, the pooled SSR is:

SSR(k) = SUM_^N [SSR_i(1:k) + SSR_i(k+1:T)]

where SSR_i(t1:t2) is the residual sum of squares from OLS of y_i on
(x_i, z_i) where z_i captures the regime indicator. The first break is
khat_1 = argmin SSR(k). Subsequent breaks are found by searching within
the largest remaining segment, sequentially subdividing.

5. SaRa Algorithm (Li, Xiao and Chen, 2025)

For each bandwidth h and time t, the local statistic is:

Ds(t,h) = max_k |sqrt(N) * (beta_R_(t,h) - beta_L_(t,h))|

where beta_L and beta_R are local OLS estimates from the left window
(t-h, t] and right window (t, t+h]. Break candidates are h-local maximizers
of Ds(t,h) across all bandwidths. These are filtered by a median threshold
and the final number is selected by BIC.

6. Post-Estimation (all methods)

After break detection, regime-specific coefficients are re-estimated by
restricted OLS. For each regime j (from tau_+1 to tau_j):

alpha_j = [SUM_i SUM_ x_it x'_it]^(-1) * [SUM_i SUM_ x_it y_it]

Standard errors use the Eicker-Huber-White sandwich estimator:

V(alpha_j) = Phi^(-1) * Omega_j * Phi^(-1) / (N * tau_j)

where Phi = E(x x') and Omega_j = E(x x' * epsilon^2) in regime j.

Output Tables

xtbreakmodel produces publication-quality output tables:

1. Header: Model specifications, panel dimensions, method-specific
settings.

2. Break detection: Number of breaks detected, break dates. For
GAGFL, this is reported per group.

3. Group membership (GAGFL only): Number of units and breaks
per group.

4. Coefficient table: Regime-specific estimates with standard
errors, z-statistics, p-values, and significance stars (* p<0.10,
** p<0.05, *** p<0.01). Each regressor is shown with its regime
interval.

5. Graphs: Coefficient paths with confidence bands and vertical
dashed lines at break dates. GAGFL produces one panel per group,
combined into a single figure.

Stored Results

xtbreakmodel stores the following in e():

Scalars

e(N) number of cross-sectional units

e(T) number of time periods

e(p) number of regressors

e(G) number of groups (1 for non-GAGFL methods)

e(ssr) total residual sum of squares

e(rmse) root mean squared error

e(niter) number of iterations (GAGFL only)

Matrices

e(nbreaks) 1 x G vector of break counts per group

e(regime) (max_regimes) x G matrix of regime start points

e(alpha) (max_regimes) x (K*G) matrix of regime-specific coefficients

e(se) (max_regimes) x (K*G) matrix of standard errors

e(group) N x G group membership indicator matrix (GAGFL only)

Macros

e(cmd) xtbreakmodel

e(cmdline) full command as typed

e(depvar) dependent variable name

e(indepvars) independent variable names

e(method) estimation method used

e(title) estimation title

Examples

Example 1: Basic PLS — Common breaks with single regressor

. webuse grunfeld, clear

. xtbreakmodel invest mvalue, method(pls)

Detects common structural breaks in the investment-market value relationship
across all firms. Uses default lambda grid [0.01, 100] with 40 grid points.

Example 2: PLS with multiple regressors

. webuse grunfeld, clear

. xtbreakmodel invest mvalue kstock, method(pls) nograph

Estimates common breaks where both market value and capital stock coefficients
shift simultaneously (vector break assumption).

Example 3: GAGFL — Heterogeneous breaks with latent groups

. webuse grunfeld, clear

. xtbreakmodel invest mvalue, method(gagfl) groups(2)

Classifies firms into 2 latent groups and estimates group-specific break
dates. This captures heterogeneity: e.g., some firms may have experienced
structural shifts in their investment behavior at different times.

Example 4: BFK — Sequential break detection

. webuse grunfeld, clear

. xtbreakmodel invest mvalue kstock, method(bfk) nograph

Finds up to 3 breaks sequentially via SSR minimization. Simple and fast.

Example 5: SaRa — Nonparametric screening

. webuse grunfeld, clear

. xtbreakmodel invest mvalue, method(sara) bandwidths(2 3 4) c1(0.1)

Screening and ranking with bandwidths 2, 3, 4. Robust to distributional
assumptions. The c1 parameter controls the BIC penalty for number of breaks.

Example 6: Simulated DGP with known breaks

. clear all

. set seed 12345

. local N = 100

. local T = 20

. set obs `=`N'*`T''

. gen id = ceil(_n/`T')

. bysort id: gen time = _n

. xtset id time

. gen x = rnormal()

. gen group0 = cond(id <= 33, 1, cond(id <= 66, 2, 3))

. gen beta_true = .

. * Group 1: breaks at t=10, t=17; beta = (1, 2, 3)

. replace beta_true = 1 if group0 == 1 & time < 10

. replace beta_true = 2 if group0 == 1 & time >= 10 & time < 17

. replace beta_true = 3 if group0 == 1 & time >= 17

. * Group 2: breaks at t=7, t=17; beta = (3, 4, 5)

. replace beta_true = 3 if group0 == 2 & time < 7

. replace beta_true = 4 if group0 == 2 & time >= 7 & time < 17

. replace beta_true = 5 if group0 == 2 & time >= 17

. * Group 3: no breaks; beta = 1.5

. replace beta_true = 1.5 if group0 == 3

. gen y = x * beta_true + 0.5 * rnormal()

. xtbreakmodel y x, method(gagfl) groups(3)

This DGP from Okui and Wang (2021, Section 4) has 3 groups with heterogeneous
breaks. GAGFL correctly recovers group assignments, break dates (7, 10, 17),
and regime coefficients.

Example 7: Multiple regressors with GAGFL

. clear all

. set seed 54321

. local N = 80

. local T = 15

. set obs `=`N'*`T''

. gen id = ceil(_n/`T')

. bysort id: gen time = _n

. xtset id time

. gen x1 = rnormal()

. gen x2 = rnormal()

. gen b1 = cond(time <= 8, 1, 2)

. gen b2 = cond(time <= 10, 3, -1)

. gen y = x1 * b1 + x2 * b2 + 0.5 * rnormal()

. xtbreakmodel y x1 x2, method(pls)

With two regressors having breaks at different times, PLS (common breaks)
will detect the union of break dates. For coefficient-specific analysis,
consider xtcbc.

Example 8: Comparing methods

. * Run all four methods on the same data

. xtbreakmodel y x, method(pls) nograph

. est store pls

. xtbreakmodel y x, method(bfk) nograph

. est store bfk

. xtbreakmodel y x, method(sara) nograph

. est store sara

. xtbreakmodel y x, method(gagfl) groups(3)

. est store gagfl

Example 9: Accessing stored results

. xtbreakmodel y x1 x2, method(pls) nograph

. display "Number of breaks: " e(nbreaks)[1,1]

. matrix list e(regime)

. matrix list e(alpha)

. display "RMSE = " e(rmse)

. display "Method: " e(method)

Example 10: Custom tuning for large panels

. xtbreakmodel y x, method(gagfl) groups(4) maxlambda(200) minlambda(0.001) ngrid(80) nsim(100) maxiter(30)

For large panels (N > 500), increase ngrid for finer lambda selection and nsim
for better GFE initialization. Use wider lambda range if the algorithm selects
boudary values.

Notes

Choosing the number of groups (G).

For GAGFL, the user must specify G. If unknown, estimate for G = 2, 3, ...,
G_max and compare the BIC = SSR/(NT) + sigma^2 * (np_G + N) * ln(NT) / (NT),
where np_G is the total number of group-specific parameters. Choose G
minimizing BIC. Typically, G <= 5 is sufficient for most applications.

When to use which method.

GAGFL: Use when heterogeneity across units is expected. This is the
most general and powerful method but requires specifying G and is the most
computationally intensive.

PLS: Use when all units are expected to share the same break dates.
Good default choice for homogeneous panels. Fast and well-understood
asymptotics.

BFK: Use as a quick diagnostic. Simple, fast, but limited to sequential
detection (may miss nearby breaks). Good for up to 3 breaks.

SaRa: Use when robustness to distributional assumptions is important.
Non-parametric approach. Good for detecting breaks in the presence of heavy-tailed
errors or outliers.

Sensitivity to tuning parameters.

For GAGFL/PLS, the key parameter is the IC constant (currently fixed at
c = 0.05 in the information criterion). If too many breaks are detected,
the result may be over-segmented; if too few, under-segmented. The IC is
consistent as N -> infinity with T fixed, so for small N, results should be
interpreted carefully.

Computational considerations.

GAGFL complexity: O(nGrid * maxIter * G * T * N * K). For N=100, T=20, G=3,
K=1, nGrid=40: approximately 30 seconds. For larger panels, reduce nGrid or
use PLS first for exploratory analysis.

Comparison with xtcbc.

xtbreakmodel uses vector breaks: all K coefficients break simultaneously
at each break date. xtcbc (Kaddoura, 2025) allows coefficient-by-coefficient
breaks where each regressor has its own break schedule. If you suspect different
regressors break at different times, use xtcbc instead.

Balanced panels required.

The current implementation requires strongly balanced panels. For unbalanced
panels, balance first using tsfill or imputation methods. See
xtmispanel for panel imputation tools.

Time-series operators.

xtbreakmodel supports Stata time-series operators in the variable list,
including L. (lag), D. (difference), and F. (forward). Ensure
that applying operators does not create gaps in the panel.

References

Okui, R. and Wang, Y. (2021). Heterogeneous structural breaks in panel data
models. Journal of Econometrics, 220(2), 447-473.

Qian, J. and Su, L. (2016). Shrinkage estimation of common breaks in panel
data models via adaptive group fused lasso. Journal of Econometrics,
191(1), 86-109.

Baltagi, B. H., Feng, Q. and Kao, C. (2016). Estimation of heterogeneous
panels with structural breaks. Journal of Econometrics, 191(1), 176-195.

Li, J., Xiao, Z. and Chen, J. (2025). Screening and ranking algorithm for
change-point detection in panel data. Working Paper.

Bonhomme, S. and Manresa, E. (2015). Grouped patterns of heterogeneity in
panel data. Econometrica, 83(3), 1147-1184.

Bai, J. and Perron, P. (1998). Estimating and testing linear models with
multiple structural changes. Econometrica, 66(1), 47-78.

Author

Dr Merwan Roudane

merwanroudane920@gmail.com

Please cite as:

Roudane, M. (2026). xtbreakmodel: Stata module for heterogeneous structural
breaks in panel data models.

Also see

xtcbc, xtpmg, xtlmbreak, xtset, regress