Overview
The sysbreak package provides methods for detecting,
estimating, and conducting inference on multiple structural breaks in a
linear system of equations, such as multivariate linear regressions,
seemingly unrelated regressions, vector autoregressions, and network
vector autoregressions. Breaks can occur in the regression coefficients,
the covariance matrix of the errors, or both, and can be restricted to a
subset of the coefficients.
For observation \(t\) in regime
\(j\),
\[
y_t = Z_t \beta_j + u_t,
\quad t = T_{j-1}+1,\ldots,T_j,
\quad j = 1,\ldots,m+1.
\]
The unknown break dates are \(T_1,\ldots,T_m\). The package searches over
admissible partitions subject to a trimming rule. It determines the
number of breaks, gives the estimated break dates, the regime-specific
parameters, and confidence intervals.
The main reference is Qu and Perron (2007), which develops estimation
and testing for structural changes in multivariate regressions.
Main Functions
The main user-facing functions are:
sysbreak() for a general linear system of equations,
where the number of regressors can differ across equations.varbreak() for structural breaks in VAR models.netbreak() for structural breaks in Network
models.ordbreak() for ordered breaks in two sets of
equations.brtype() for tests of what changes at each break: the
volatilities, the correlations, or the coefficients.dyspill() for generalized Diebold-Yilmaz spillover
measures from fitted VAR-type models.dynobreak() for the same spillover measures from a VAR
estimated on the whole sample without breaks, for comparison.
Basic Example
This example uses simulated data with two breaks.
library(sysbreak)
set.seed(3)
bigt <- 60
x <- rnorm(bigt)
beta <- rbind(
c(1.0, 0.5),
c(2.0, -0.2),
c(0.5, 0.8)
)
y <- numeric(bigt)
y[1:20] <- beta[1, 1] + beta[1, 2] * x[1:20] + rnorm(20, sd = 0.4)
y[21:40] <- beta[2, 1] + beta[2, 2] * x[21:40] + rnorm(20, sd = 0.4)
y[41:60] <- beta[3, 1] + beta[3, 2] * x[41:60] + rnorm(20, sd = 0.4)
df <- data.frame(y = y, x = x)
The model is given as a formula through model, with the
data supplied in data. The intercept is included
automatically. Which coefficients are allowed to change is controlled by
coef_breaks; the default "all" lets every
coefficient change across regimes. The model is estimated by:
fit <- sysbreak(y ~ x, data = df, m = 2)
The results are summarized by:
summary(fit)
#> sysbreak: 1 equation, T = 60, m = 2, trim = 0.15, robust = none
#> Covariance matrix: constant across regimes
#>
#> Tests for the number of breaks (5% level)
#> Statistic Critical value
#> WDmax 77.595 12.594 *
#> SEQ(2|1) 63.080 13.472 *
#> Breaks selected: 2
#>
#> Break dates (rows of the data) with 95% confidence intervals
#> Estimate Lower Upper
#> Break 1 20 19 21
#> Break 2 40 39 41
#>
#> Regime 1 (rows 1 to 20)
#> Estimate Std. Error z value Pr(>|z|)
#> y:(Intercept) 1.10076 0.06468 17.018 < 2e-16 ***
#> y:x 0.45801 0.08286 5.527 3.25e-08 ***
#>
#> Regime 2 (rows 21 to 40)
#> Estimate Std. Error z value Pr(>|z|)
#> y:(Intercept) 2.04581 0.06330 32.32 < 2e-16 ***
#> y:x -0.35734 0.07308 -4.89 1.01e-06 ***
#>
#> Regime 3 (rows 41 to 60)
#> Estimate Std. Error z value Pr(>|z|)
#> y:(Intercept) 0.51638 0.06426 8.036 9.29e-16 ***
#> y:x 0.75479 0.06037 12.504 < 2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The output reports the tests used to select the number of breaks
(WDmax and the sequential test of one against two breaks), the selected
number of breaks, the break dates with their confidence intervals, and
the coefficient estimates in each regime. Both tests reject, and two
breaks are selected. The break dates and their confidence intervals are
reported as rows of df. They are returned by
confint(), and the regime coefficients by
coef(); their covariance matrix is returned by
vcov():
confint(fit)
#> 2.5 % break 97.5 %
#> break1 19 20 21
#> break2 39 40 41
coef(fit)
#> regime1:y:(Intercept) regime1:y:x regime2:y:(Intercept)
#> 1.1007640 0.4580116 2.0458059
#> regime2:y:x regime3:y:(Intercept) regime3:y:x
#> -0.3573438 0.5163846 0.7547884
The fit can be plotted. Each series is plotted with the estimated
break dates, their confidence intervals and the fitted values:
The SupLR statistics for one and two breaks, together with the other
test results, are stored in fit$testing; see the next
section.
Output Structure
The methods shown above extract the main results;
residuals(), fitted() and nobs()
are also available. The components returned by sysbreak(),
varbreak(), and netbreak() are:
spec: the fitted specification, including
m, trim, sample size, equation and coefficient
dimensions, break options, and coefficient and observation names.testing: the tests, reported at
test_level: statistic, testing no break
against \(k\) breaks for \(k = 1,\ldots,m\), with its
critical_values; statistic_type
("SupLR", "SupF", or missing when a requested
robust test is unavailable); wdmax and
wdmax_critical_values; the sequential tests of
\(l\) against \(l+1\) breaks, each containing
l, a statistic and its
critical_values; and the number of
changing_parameters.br_selected: the selected number of breaks at
test_level, with the WDmax and sequential decisions that
led to that result.estimation: the estimates for the selected number of
breaks: breaks, regime-stacked coefficients,
their covariance matrix vcov, stacked error
covariance matrices, break-date confidence intervals
ci_breaks, and the residuals. It is
NULL when no break is selected; coef(),
vcov(), confint() and residuals()
then stop with a message that points to testing.y: the responses used in estimation.
varbreak() and netbreak() also return
spillover when h is supplied and at least one
break is selected.
Restricted Estimation
Instead of letting every coefficient change, the breaks can be
restricted to a subset of the coefficients. The argument
coef_breaks then takes a list of formulas, one per
equation: the left-hand side identifies the equation and the right-hand
side the coefficients that are allowed to change. All other coefficients
are held constant across regimes.
To illustrate, the data below keep the intercept at 1 in all regimes
and let only the slope on x change, at the same dates as
before:
set.seed(13)
y_r <- 1 + rep(c(0.5, -0.2, 0.8), each = 20) * x + rnorm(bigt, sd = 0.4)
df_r <- data.frame(y = y_r, x = x)
The model with only the slope allowed to change is estimated by:
restricted_fit <- sysbreak(y ~ x, data = df_r, m = 2, coef_breaks = list(y ~ x))
summary(restricted_fit)
#> sysbreak: 1 equation, T = 60, m = 2, trim = 0.15, robust = none
#> Covariance matrix: constant across regimes
#>
#> Tests for the number of breaks (5% level)
#> Statistic Critical value
#> WDmax 33.034 9.567 *
#> SEQ(2|1) 24.074 10.187 *
#> Breaks selected: 2
#>
#> Break dates (rows of the data) with 95% confidence intervals
#> Estimate Lower Upper
#> Break 1 20 17 23
#> Break 2 40 38 42
#>
#> Coefficients allowed to change: y:x
#>
#> Regime 1 (rows 1 to 20)
#> Estimate Std. Error z value Pr(>|z|)
#> y:(Intercept) 1.01180 0.04807 21.048 < 2e-16 ***
#> y:x 0.55702 0.10600 5.255 1.48e-07 ***
#>
#> Regime 2 (rows 21 to 40)
#> Estimate Std. Error z value Pr(>|z|)
#> y:(Intercept) 1.01180 0.04807 21.048 < 2e-16 ***
#> y:x -0.30795 0.09485 -3.247 0.00117 **
#>
#> Regime 3 (rows 41 to 60)
#> Estimate Std. Error z value Pr(>|z|)
#> y:(Intercept) 1.01180 0.04807 21.05 <2e-16 ***
#> y:x 0.80879 0.07757 10.43 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Two breaks are selected, at rows 20 and 40. The critical values are
smaller than in the basic example, because fewer coefficients change.
The line Coefficients allowed to change lists the
coefficients that change; the intercept is estimated from the whole
sample and is therefore the same in the three regimes. The returned
components describe this restricted model; an unrestricted fit is not
returned.
A pure covariance-break model (no coefficient changes) is obtained
with coef_breaks = "none" together with
brv = TRUE. With robust = "hc" or
"hac", a covariance-only robust test is not implemented, so
the test statistic and its type are returned as missing.
Further Options
The examples above use the default settings. The following arguments,
common to sysbreak(), varbreak() and
netbreak(), can be used to alter how the model is estimated
and how inference is carried out.
brv: with brv = TRUE, the covariance
matrix of the errors is allowed to change at the break dates. The
default holds it constant.robust: the default "none" gives the
likelihood-ratio tests, valid when the errors are serially uncorrelated
and conditionally homoskedastic. "hc" gives the robust SupF
test of Qu and Perron (2007, eq. 20) with a
heteroskedasticity-consistent covariance estimate, valid when the errors
are a martingale difference sequence. "hac" also allows for
serial correlation, using the quadratic spectral kernel with the
automatic bandwidth of Andrews (1991). The choice also applies to the
coefficient standard errors and the break-date confidence intervals.
With robust = "hac", prewhit = TRUE prewhitens
the scores by a VAR(1) before the HAC estimation.hetq: with hetq = TRUE, the regressor
moments may differ across regimes when the break-date confidence
intervals are computed. This is recommended when the distribution of the
regressors changes across regimes, as in a VAR.restart: the restricted estimator of the previous
section is iterative and can stop at a local optimum.
restart = "brbp" (resampling) or
"brbp_residual" (residuals) restarts it from
restart_B bootstrap samples (10 by default) and keeps the
fit associated with the highest likelihood. It applies only when some
coefficients are held constant.trim: the minimum regime length is
round(trim * T). At least 3(n + k) observations are
recommended, where n is the number of equations and k the number of
regressors per equation.ci_level: the coverage of the break-date confidence
intervals, 0.95 (the default) or 0.90.
For example, the basic model with inference robust to serial
correlation and 90% confidence intervals is estimated by:
fit_hac <- sysbreak(y ~ x, data = df, m = 2, robust = "hac", ci_level = 0.90)
VAR Example and Types of Breaks
varbreak() estimates structural breaks in a VAR. Its
first argument is the matrix of the series, and the lag order is set by
lags. The data below follow a bivariate VAR(1) whose
coefficients are the same throughout, while at row 100 the volatility of
y1 doubles and the correlation of the errors rises from 0.2
to 0.7:
set.seed(14)
bigt_var <- 200
A <- matrix(c(0.5, 0.1, 0.2, 0.4), 2, 2)
chol1 <- t(chol(matrix(c(1, 0.2, 0.2, 1), 2, 2)))
chol2 <- t(chol(matrix(c(4, 1.4, 1.4, 1), 2, 2)))
y_var <- matrix(0, bigt_var, 2, dimnames = list(NULL, c("y1", "y2")))
for (tt in 2:bigt_var) {
chol_t <- if (tt <= 100) chol1 else chol2
y_var[tt, ] <- A %*% y_var[tt - 1, ] + chol_t %*% rnorm(2)
}
The model lets both the coefficients and the covariance matrix change
(brv = TRUE). hetq = TRUE is used because the
distribution of the lagged series changes with the volatility:
var_fit <- varbreak(y_var, lags = 1, m = 2, brv = TRUE, hetq = TRUE)
confint(var_fit)
#> 2.5 % break 97.5 %
#> break1 99 100 115
One break is selected, at row 100. brtype() then tests,
at each break date, whether the volatilities, the correlations and the
coefficients differ between the two adjacent regimes:
round(brtype(var_fit), 4)
#> breaks volatility volatility_p correlation correlation_p coefficients
#> 1 100 55.5739 0 20.7622 0 4.9234
#> coefficients_p
#> 1 0.5537
The volatility and correlation tests reject, with p-values below
0.0001, while the coefficient test does not (p-value 0.55). The break is
in the covariance matrix only, as in the simulated data.
Network VAR Example
netbreak() estimates structural breaks in a Network VAR
(Zhu et al., 2017) whose transition matrices are restricted by a fixed
network matrix:
\[
y_t =
\mu_j +
\sum_{\ell=1}^p
(\alpha_{\ell,j} W + \beta_{\ell,j} I)y_{t-\ell}
+ u_t .
\]
The argument coef_breaks controls which low-dimensional
Network VAR parameters are allowed to change across regimes.
set.seed(5)
bigt_nvar <- 80
W <- matrix(c(0, 1, 1, 0), 2, 2, byrow = TRUE)
y_nvar <- matrix(0, bigt_nvar, 2)
colnames(y_nvar) <- c("node1", "node2")
for (tt in 2:bigt_nvar) {
alpha <- if (tt <= 40) 0.05 else 0.75
beta <- 0.1
y_nvar[tt, ] <- alpha * as.numeric(W %*% y_nvar[tt - 1, ]) +
beta * y_nvar[tt - 1, ] + rnorm(2, sd = 0.2)
}
nvar_fit <- netbreak(
y = y_nvar,
W = W,
lags = 1,
m = 1,
coef_breaks = c("network", "own"),
trim = 0.2,
h = 10
)
nvar_fit$estimation$breaks
#> [1] 45
nvar_fit$spillover$regimes[[1]]$total
#> [1] 2.193206
The same spillover engine can be called after estimation:
spill <- dyspill(nvar_fit, h = 10)
spill$regimes[[1]]$fevd
spill$regimes[[1]]$from
spill$regimes[[1]]$to
spill$regimes[[1]]$net
Ordered Two-Break Example
ordbreak() estimates an ordered two-break model in which
the equations in set1 break at the first date and those in
set2 at the second. The test is that of Qu and Perron
(2007, Section 6.4) for locally ordered breaks. The test for globally
ordered breaks in their Section 6.5, whose limiting distribution
differs, is not implemented.
set.seed(13)
bigt2 <- 60
x1 <- rnorm(bigt2)
x2 <- rnorm(bigt2)
y1 <- 1 + 0.5 * x1 + rnorm(bigt2, sd = 0.3)
y2 <- 1 + 0.5 * x2 + rnorm(bigt2, sd = 0.3)
y1[21:bigt2] <- y1[21:bigt2] + 1.5 * x1[21:bigt2]
y2[36:bigt2] <- y2[36:bigt2] - 1.5 * x2[36:bigt2]
df2 <- data.frame(y1 = y1, y2 = y2, x1 = x1, x2 = x2)
ordered_fit <- ordbreak(y1 ~ x1, y2 ~ x2, data = df2, trim = 0.2)
ordered_fit$estimation$breaks
#> [1] 20 35
Data Sets
The package includes five data sets, each documented on its help page
(for example, ?housing):
housing: quarterly U.S. data from 1959:Q4 to 2009:Q3 on
real GDP, consumption, real home prices and the federal funds rate, for
the house price application of Ahamada and Diaz Sanchez (2013), to be
used with sysbreak().inflation: monthly inflation for the G7 countries and
the euro area from 1973:04 to 2007:12, from Bataa, Osborn, Sensier and
van Dijk (2013), to be used with varbreak().netsim: one simulated draw from Scenario 1 of Han and
Lee (2024), a Network VAR for 20 nodes over 300 periods with breaks
after periods 100 and 200, to be used with netbreak(); it
contains the series y and the network matrix
W.data1 and data2: the example data of the
original GAUSS programs of Qu and Perron (2007), for
sysbreak() and ordbreak(), respectively.
Critical Values
Critical values are computed internally, following Bai and Perron
(2003), and reported at the significance level test_level
(default 0.05). In the basic example above,
fit$testing$critical_values has one entry for each
candidate number of breaks \(1,\ldots,m\); the critical value of the
WDmax test is fit$testing$wdmax_critical_values, and that
of each sequential test is stored in the corresponding
fit$testing$sequential[[l]] entry. Another level is
obtained by refitting with a different test_level (0.10,
0.05, 0.025, or 0.01).
References
Qu, Z., and Perron, P. (2007). Estimating and testing structural
changes in multivariate regressions. Econometrica, 75(2),
459-502.
Bataa, E., Osborn, D. R., Sensier, M., and van Dijk, D. (2013).
Structural breaks in the international dynamics of inflation. Review
of Economics and Statistics, 95(2), 646-659.
Bai, J., and Perron, P. (2003). Critical values for multiple
structural change tests. The Econometrics Journal, 6(1),
72-78.
Zhu, X., Pan, R., Li, G., Liu, Y., and Wang, H. (2017). Network
vector autoregression. The Annals of Statistics, 45(3),
1096-1123.
Diebold, F. X., and Yilmaz, K. (2012). Better to give than to
receive: Predictive directional measurement of volatility spillovers.
International Journal of Forecasting, 28(1), 57-66.
Ahamada, I., and Diaz Sanchez, J. L. (2013). A retrospective analysis
of the house price macro-relationship in the United States.
International Journal of Central Banking, 9(4), 153-174.
Han, Y., and Lee, T. C. M. (2024). Structural break detection in
non-stationary network vector autoregression models. IEEE
Transactions on Network Science and Engineering, 11(5),
4134-4145.
Andrews, D. W. K. (1991). Heteroskedasticity and autocorrelation
consistent covariance matrix estimation. Econometrica, 59(3),
817-858.
Session Info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] sysbreak_1.1.0 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.39 R6_2.6.1 fastmap_1.2.0 xfun_0.62
#> [5] maketools_1.3.2 cachem_1.1.0 knitr_1.52 htmltools_0.5.9
#> [9] buildtools_1.0.0 lifecycle_1.0.5 cli_3.6.6 sass_0.4.10
#> [13] jquerylib_0.1.4 compiler_4.6.1 sys_3.4.3 tools_4.6.1
#> [17] evaluate_1.0.5 bslib_0.12.0 yaml_2.3.12 otel_0.2.0
#> [21] jsonlite_2.0.0 rlang_1.3.0 MASS_7.3-66