QR.break: An R Package for Structural Breaks in Quantile Regression

2.1 Econometric models and issues of interest

2.2 Estimating break locations based on a single quantile

Suppose the \(\tau\)th quantile is affected by \(m\) structural changes. For a set of candidate break dates \(T^b = (T_1,\dots,T_m)\), the break dates and the coefficients are estimated jointly by minimizing the check function over admissible partitions:

\[ (\hat{\beta}(\tau ),\hat{T}^{b})=\arg \min_{\beta (\tau ),\,T^{b}\in \Lambda _{\varepsilon }}\ \sum_{j=0}^{m}\ \sum_{t=T_{j}+1}^{T_{j+1}}\ \sum_{i=1}^{N}\rho _{\tau }(y_{it}-x_{it}^{\prime }\beta _{j+1}(\tau )), \]

where \(\rho_\tau(u) = u(\tau - 1(u < 0))\) is the check function (Koenker, 2005), \(\beta(\tau) = (\beta_1(\tau)',\dots,\beta_{m+1}(\tau)')'\), \(T_0 = 0\), \(T_{m+1} = T\), and the summation over \(i\) is absent for a single time series. The set \(\Lambda_\varepsilon\) collects the admissible partitions, ensuring each estimated regime is a positive fraction of the sample: every regime must be at least \(\varepsilon T\) periods long, and no break may fall within \(\varepsilon T\) of either end of the sample. The user specifies \(\varepsilon\) through the argument trim.e.

The computation is carried out using a dynamic programming algorithm as in Bai and Perron (2003), so that the computation is of order \(O(T^2)\) irrespective of the number of breaks allowed in the model.

2.3 Estimating break locations based on multiple quantiles

Suppose the quantiles in \(\mathcal{T}_\omega = [\omega_1, \omega_2]\) with \(0 < \omega_1 < \omega_2 < 1\) are affected by structural changes. A natural approach is to consider a partition of this interval, examine a set of quantiles \(\tau_h\), \(h = 1,\dots,q\), and minimize the check function summed over them as well:

\[ (\hat{\beta}(\mathcal{T}_{\omega }),\hat{T}^{b})=\arg \min_{\beta (\mathcal{T}% _{\omega }),\,T^{b}\in \Lambda _{\varepsilon }}\ \sum_{h=1}^{q}\ \sum_{j=0}^{m}\ \sum_{t=T_{j}+1}^{T_{j+1}}\ \sum_{i=1}^{N}\rho _{\tau _{h}}(y_{it}-x_{it}^{\prime }\beta _{j+1}(\tau _{h})), \]

with \(\Lambda_\varepsilon\) as before. The estimated break dates are now common across quantiles, while the coefficients remain quantile-specific. Regarding the partition, a coarse partition, such as quantiles spaced by 0.1, is sufficient to deliver informative results. The computation is once more \(O(T^2)\) irrespective of the number of breaks allowed.

For all cases above, confidence intervals for the break dates are computed from the limiting distribution of the break point estimates derived in Oka and Qu (2011). That distribution has an analytical density function, so no simulation is needed to obtain the critical values, which keeps the computational cost down; see Bai (1995) and Oka and Qu (2011) for details.

2.4 Determining the number of breaks

The package uses the two test statistics \(SQ_\tau\) (single quantile) and \(DQ\) (multiple quantiles) proposed in Qu (2008).

2.5 Inference on quantile regression coefficients after estimating the breaks

After determining the break dates, the coefficients are estimated conditional on those dates. This is equivalent to partitioning the sample using the break dates and running standard quantile regression on each subsample. These results are reported automatically by the main function.

Oka and Qu (2011) showed that if the break size is of higher order than \(T^{-1/2}\), so that breaks are not confounded with estimation uncertainty of order \(T^{-1/2}\), then the break dates are estimated fast enough that the asymptotic distribution of the estimated coefficients is the same as if the break dates were known — a familiar result in the structural break literature. This permits computing confidence intervals treating the estimated break dates as known. The asymptotic variance involves the within-regime second moment matrix of the regressors and the same matrix weighted by the conditional density of the dependent variable at the quantile of interest; the density is estimated by a difference quotient, as described in Qu (2008, pp. 176–177). Instead of using the package, users can obtain the same confidence intervals directly from the quantreg package with the standard error option set to "nid".

2.6 Normalization of the subgradient process (version 1.0.3)

The \(SQ\) and \(DQ\) statistics normalize the subgradient process by an inverse square root of \(W = \sum_t \sum_i x_{it} x_{it}'\) (the sums over \(i\) being absent for a single time series). Such a matrix is not unique, and the argument norm.method selects which one is used.

"cholesky" (the default) uses \((R')^{-1}\), where \(R\) is the Cholesky factor returned by chol(), with \(R'R = W\). This reproduces the results of versions 1.0.2 and earlier exactly. Because the factor is constructed sequentially, it implicitly gives more weight to regressors appearing earlier in the regression. It is suitable when there is a natural ordering — for instance a quantile autoregression in which low lags are listed first, as in the GDP application. It is invariant to the units in which the regressors are measured and to their signs.

"spectral" uses \(C^{-1/2}D^{-1}\), where \(D = \mathrm{diag}(\sqrt{\mathrm{diag}(W)})\) and \(C = D^{-1}WD^{-1}\) is the correlation matrix of the regressors, with \(C^{-1/2}\) its symmetric square root from the spectral decomposition (Kessy, Lewin and Strimmer, 2018). This treats the regressors symmetrically, yielding invariance to their order as well as to their units and signs. It is the appropriate choice when the regressors have no natural ordering, as with the age, gender and winter regressors of the second application.

Both choices give the same limiting null distribution, so the same critical values apply and both tests are valid. The second option might yield lower power. Estimated break dates are unaffected by this argument, since they are obtained by minimizing the check function, which does not depend on the normalization. A numerical comparison is given in Section 4.7.

These invariance properties refer to the normalization itself. Because the quantile regression is re-estimated, the computed statistics can still differ slightly when the regressors are reordered or rescaled, by an amount that decreases with the sample size.

3.1 Arguments

Whether or not verbose is set, all estimation and testing results are saved in the returned object.

3.2 How the data must be arranged

y is a single vector of length \(NT\), stacked by period: the first \(N\) elements are the units observed in the first period, the next \(N\) the second period, and so forth. x is stacked identically. For a time series \(N = 1\) and y is simply the series in time order.

vec.time is different: it has length \(T\), not \(NT\), being one label per period. For repeated cross-sections use unique() to reduce the period column accordingly.

A break date is the final date of the existing regime, not the starting date of a new regime. An estimate of 1984 Q1 means the first regime ends in 1984 Q1 and the second begins in 1984 Q2.

3.3 The return object

The result is a list with two components: $s.out, based on individual quantiles, and $m.out, based on all quantiles jointly. Entries in $s.out are ordered according to vec.tau, and for each quantile the break testing results come first; if at least one significant break is detected, the break locations and parameter estimates follow. Element names carry the quantile they refer to:

4.1 The data

After loading the package, the U.S. GDP data are loaded with data(gdp). The yq column contains the dates, gdp the dependent variable, and the remaining two columns the first and second lagged values of gdp, used as regressors. The estimation sample runs from 1947 Q4 to 2009 Q2 — the raw series begins in 1947 Q2, with the first two observations consumed by the lags.

data(gdp)
str(gdp)
#> 'data.frame':    247 obs. of  4 variables:
#>  $ yq  : chr  "1947 Q4" "1948 Q1" "1948 Q2" "1948 Q3" ...
#>  $ gdp : num  6.066 6.357 7.363 2.15 0.614 ...
#>  $ lag1: num  -0.329 6.066 6.357 7.363 2.15 ...
#>  $ lag2: num  -0.613 -0.329 6.066 6.357 7.363 ...
head(gdp, 3)
#>        yq     gdp     lag1     lag2
#> 1 1947 Q4 6.06552 -0.32868 -0.61348
#> 2 1948 Q1 6.35664  6.06552 -0.32868
#> 3 1948 Q2 7.36332  6.35664  6.06552

4.2 Setting the inputs

The following commands define y and x:

y <- gdp[, "gdp"]
x <- gdp[, c("lag1", "lag2")]

A column of ones is always added to the regressors, so the model has three parameters allowed to be affected by structural breaks: the intercept, the coefficient on the first lag, and that on the second lag.

The next command specifies the quantiles of interest:

vec.tau <- seq(0.20, 0.80, by = 0.15)
vec.tau
#> [1] 0.20 0.35 0.50 0.65 0.80

Given these inputs the function performs two sets of calculations. First it analyzes the quantiles in vec.tau independently, so the number of breaks and their locations may differ across quantiles. Then it analyzes all quantiles simultaneously, assuming breaks are common across quantiles and pooling information to estimate the break dates.

Since this is a time series regression, N = 1. The minimum regime length matters because if a regime is too short the model fit may be non-unique and the estimation may pick up spurious breaks; trim.e = 0.15 implies a regime is at least 15% of the sample, or 37 quarters. Related to this is the maximum number of breaks, m.max = 3, allowing at most three breaks and hence four regimes. Finally v.a = 2 selects the 5% level for determining the number of breaks and v.b = 2 selects 95% confidence intervals for the break dates.

Setting vec.time = gdp[, "yq"] makes the program report break dates using those labels; setting it to NULL reports integers instead.

4.3 Running the function

res <- rq.break(y, x,
                vec.tau     = vec.tau,
                N           = 1,
                trim.e      = 0.15,
                vec.time    = gdp[, "yq"],
                m.max       = 3,
                v.a         = 2,
                v.b         = 2,
                verbose     = TRUE,
                norm.method = "cholesky")

4.4 Output based on separate quantiles

For \(\tau = 0.2\) the results are:

$s.out$test_0.2                     # testing results at the chosen level
                1 Breaks 2 Breaks 3 Breaks
SQ test         1.423269 1.373012        0
Critical values 1.529859 1.637547        0

$s.out$nbreak_0.2                   # number of breaks detected
[1] 0

For this quantile the test of no break against a single break equals 1.423269, while the critical value is 1.529859. The value is insignificant at the chosen level, so no break is detected. The final column is zero because the test is not computed when the previous tests are insignificant at the 10% level. No break estimation results are produced.

The results for the other quantiles also show insignificance until \(\tau = 0.65\). At that point one significant break is detected, followed by its confidence interval reported first as index values and then as dates:

$s.out$test_0.65
                1 Breaks 2 Breaks 3 Breaks
SQ test         1.817933 1.023126        0
Critical values 1.529859 1.637547        0

$s.out$nbreak_0.65
[1] 1

$s.out$br_est_0.65
        Estimate CI_Lower_Bound CI_Upper_Bound
Break 1      147             83            161

$s.out$br_est_time_0.65
        Estimate  CI_Lower_Bound CI_Upper_Bound
Break 1 "1984 Q2" "1968 Q2"      "1987 Q4"

One break is detected with a point estimate of \(t = 147\), corresponding to the second quarter of 1984, with a confidence interval of \([83, 161]\), or equivalently [1968 Q2, 1987 Q4]. The output also includes the estimated coefficients for the two regimes:

$s.out$coef_0.65$Regime_1
               Value Std. Error   t value     Pr(>|t|)
Intercept  4.5169837 0.63003316  7.169438 3.625833e-11
x1         0.4199469 0.09643037  4.354924 2.512308e-05
x2        -0.1051057 0.09458095 -1.111277 2.683010e-01

$s.out$coef_0.65$Regime_2
              Value Std. Error  t value     Pr(>|t|)
Intercept 2.2855100  0.4866286 4.696620 8.711282e-06
x1        0.1724133  0.1218581 1.414869 1.603091e-01
x2        0.2386392  0.1320517 1.807164 7.383669e-02

The sum of the autoregressive coefficients changes little, but the intercept decreases significantly, indicating a notably lower 0.65 quantile after the break:

$s.out$bsize_0.65_Regime_2_minus_Regime_1
               Value Std. Error   t value    Pr(>|t|)
Intercept -2.2314737  0.7957371 -2.804285 0.005453863
x1        -0.2475336  0.1375265 -1.799898 0.073127343
x2         0.3437448  0.1418386  2.423492 0.016109478

The method similarly detects a break in the 0.80 quantile, where the pattern is sharper still:

$s.out$test_0.8
                1 Breaks 2 Breaks 3 Breaks
SQ test         2.170968 1.064567        0
Critical values 1.529859 1.637547        0

$s.out$br_est_time_0.8
        Estimate  CI_Lower_Bound CI_Upper_Bound
Break 1 "1984 Q1" "1978 Q4"      "1990 Q1"

$s.out$bsize_0.8_Regime_2_minus_Regime_1
               Value Std. Error   t value     Pr(>|t|)
Intercept -3.0893051  0.7719650 -4.001872 8.367631e-05
x1        -0.2108298  0.1159821 -1.817778 7.033963e-02
x2         0.4047554  0.1420980  2.848425 4.773147e-03

Collecting the number of breaks across quantiles:

At the median the statistic is 1.001 against a critical value of 1.530 — no evidence of a break at all — while the upper quantiles reject clearly.

4.5 Output based on multiple quantiles

The results are structured similarly, now based on all chosen quantiles:

$m.out$test_joint
                 1 Breaks  2 Breaks 3 Breaks
DQ test         1.0275870 0.5892746        0
Critical values 0.9098714 0.9584567        0

$m.out$nbreak_joint
[1] 1

$m.out$br_est_joint
        Estimate CI_Lower_Bound CI_Upper_Bound
Break 1      146            120            147

$m.out$br_est_joint_time
        Estimate  CI_Lower_Bound CI_Upper_Bound
Break 1 "1984 Q1" "1977 Q3"      "1984 Q2"

A single break is detected, as in the analysis based on the individual upper quantiles. Note that the joint confidence interval, [1977 Q3, 1984 Q2], is considerably tighter than the single-quantile interval at \(\tau = 0.65\), [1968 Q2, 1987 Q4]: pooling information across quantiles buys precision in the date. The rest of the output contains the coefficient estimates and their confidence intervals as in the single-quantile case.

In summary, the findings shed light on the “Great Moderation” debate on U.S. GDP growth. The decline in volatility mainly affected the upper tail, with the median and lower quantiles remaining stable. This implies that expansions became less rapid, while recessions remained as severe. A test on the conditional mean, which averages over the distribution, would understate the change and misdescribe its nature.

4.6 A picture

tt  <- seq_len(nrow(gdp))
brk <- 146; lo <- 120; hi <- 147     # DQ estimate and its 95% interval

op <- par(mar = c(3.5, 4, 2.5, 1))
plot(tt, gdp$gdp, type = "n", xaxt = "n", bty = "n", ylim = c(-12, 18),
     xlab = "", ylab = "Real GDP growth (%, annualized)")
rect(lo, -12, hi, 18, col = "#EDE7F6", border = NA)
abline(h = 0, col = "#CFCFD4")
lines(tt, gdp$gdp, col = "#5A5A66", lwd = 1.4)
segments(brk, -12, brk, 15, col = "#6C4FB8", lwd = 2)
text(brk, 16.5, "  break: 1984 Q1", adj = c(0, 0.5), col = "#6C4FB8", cex = 0.85)
text(lo,  16.5, "95% CI  ",         adj = c(1, 0.5), col = "#8E7BC6", cex = 0.8)
at <- seq(2, nrow(gdp), by = 40)
axis(1, at = at, labels = sub(" Q[1-4]$", "", gdp$yq[at]), col = "#CFCFD4")

par(op)

The shaded band is the 95% confidence interval for the break date and the vertical rule the point estimate. The visible narrowing of the series after it is the compression that the upper-quantile results quantify.

4.7 Comparing the two normalizations

Re-running the analysis with norm.method = "spectral" changes the statistics somewhat but not the conclusions:

Every rejection decision is unchanged, and the estimated break date remains 1984 Q1 with the identical confidence interval [1977 Q3, 1984 Q2] — as it must, since break dates do not depend on the normalization.

4.8 Computational time

The dynamic programming algorithm determines the globally optimal break partitions at a cost growing with the square of the sample size, regardless of the number of breaks allowed. In this example the program finishes within a few minutes on a typical desktop computer with a single processor. Because the quantile trimming here is symmetric (\(\omega_1 = 0.2\), \(\omega_2 = 0.8\)) and there are only three coefficients, the \(DQ\) critical values come from the response surface and no simulation is required.

4.9 Potential error messages

The function displays error messages if the inputs are not configured properly. Suppose we set trim.e = 0.2 and m.max = 6. The product exceeds 1, because it is not possible to allow six breaks when each regime is at least 20% of the sample:

rq.break(y, x, vec.tau, N = 1, trim.e = 0.2, vec.time = gdp[, "yq"],
         m.max = 6, v.a = 2, v.b = 2)
#> Error in rq.break(y, x, vec.tau, N = 1, trim.e = 0.2, vec.time = gdp[, : Error: m.max*trim.e exceeds 1. This occurs because too many regimes are allowed or the minimum length of a regime is too large. Consider decreasing m.max, trim.e, or both.

As another example, if trim.e is too small the message suggests increasing it:

rq.break(y, x, vec.tau, N = 1, trim.e = 0.01, vec.time = gdp[, "yq"],
         m.max = 3, v.a = 2, v.b = 2)
#> Error in rq.break(y, x, vec.tau, N = 1, trim.e = 0.01, vec.time = gdp[, : Error: trim.e * nrow(x) must be at least the number of regressors; otherwise, the estimation results are not unique. Consider increasing trim.e.

When an error message is produced the program exits with no saved results, and the user can modify the inputs and restart.

5.1 The data

The main steps and the output structure are similar to the time series case. The data record the blood alcohol concentration of drivers under 21 involved in motor vehicle accidents in California from 1983 Q1 to 2007 Q4, together with age, a gender dummy, and a dummy for the fourth quarter. There are \(N = 108\) drivers per quarter over \(T = 100\) quarters, giving \(NT = 10{,}800\) observations. The data are organized so that the first \(N\) rows are for the first time period, the next \(N\) rows for the next period, and so forth.

data(driver)
str(driver)
#> 'data.frame':    10800 obs. of  5 variables:
#>  $ yq    : chr  "1983 Q1" "1983 Q1" "1983 Q1" "1983 Q1" ...
#>  $ bac   : num  0.23 0 0 0 0 0 0.04 0.22 0 0 ...
#>  $ age   : int  20 17 20 20 18 16 20 19 20 20 ...
#>  $ gender: int  1 1 1 1 1 1 1 1 0 1 ...
#>  $ winter: int  0 0 0 0 0 0 0 0 0 0 ...

y        <- driver[, "bac"]
x        <- driver[, c("age", "gender", "winter")]
vec.time <- unique(driver[, "yq"])    # length T = 100, one label per quarter
length(vec.time)
#> [1] 100

The choice of quantiles is dictated by the data. BAC is exactly zero for most of these drivers:

mean(driver$bac == 0)
#> [1] 0.6187037
quantile(driver$bac, c(0.50, 0.60, 0.65, 0.70, 0.80, 0.85))
#>  50%  60%  65%  70%  80%  85% 
#> 0.00 0.00 0.04 0.07 0.13 0.16

About 62% of the observations sit at zero, so every quantile below roughly \(\tau = 0.62\) falls inside that point mass and carries no information about drinking behaviour. Attention is therefore restricted to \(\tau \in \{0.70, 0.75, 0.80, 0.85\}\) — high enough to clear the mass point, low enough to be estimated reliably.

5.2 Running the function

res.d <- rq.break(y, x,
                  vec.tau     = seq(0.70, 0.85, by = 0.05),
                  N           = 108,
                  trim.e      = 0.05,
                  vec.time    = vec.time,
                  m.max       = 3,
                  v.a         = 2,
                  v.b         = 2,
                  verbose     = TRUE,
                  norm.method = "cholesky")

Two features differ from the time series case. First, N = 108 changes the meaning of the sample: regime lengths, trimming and break dates are counted in periods. Second, the quantile range is asymmetric (\(\omega_1 = 0.70\), \(\omega_2 = 0.85\), so \(\omega_2 \neq 1 - \omega_1\)), so the response surface does not apply and the function prints

The critical values of the DQ test is obtained via simulations

before generating them by simulation, with 50,000 replications, for this specific design matrix. This is accurate but noticeably slower, and is the dominant cost of this example.

5.3 Output based on separate quantiles

The first quantile is \(\tau = 0.7\), for which two breaks are detected along with their confidence intervals:

$s.out$test_0.7
                1 Breaks 2 Breaks 3 Breaks
SQ test         5.179373 2.205150 1.339123
Critical values 1.574681 1.679331 1.737955

$s.out$nbreak_0.7
[1] 2

$s.out$br_est_0.7
        Estimate CI_Lower_Bound CI_Upper_Bound
Break 1       10              6             15
Break 2       38             35             39

$s.out$br_est_time_0.7
        Estimate  CI_Lower_Bound CI_Upper_Bound
Break 1 "1985 Q2" "1984 Q2"      "1986 Q3"
Break 2 "1992 Q2" "1991 Q3"      "1992 Q3"

The tests reject far more decisively than in the time series case — the statistics are several times their critical values — which is what a cross-section of 108 per period buys in power. The output continues with the parameter estimates for the three regimes and then moves to the next quantile.

The number of breaks varies across the range:

At \(\tau = 0.80\) and \(0.85\) only the 1992 break is found; the mid-1980s break appears only at the two lower quantiles. This is the pattern behind the substantive conclusion below.

5.4 Output based on multiple quantiles

$m.out$test_joint
                 1 Breaks  2 Breaks  3 Breaks
DQ test         2.3734870 1.0105268 0.5995825
Critical values 0.7717347 0.8162768 0.8428049

$m.out$nbreak_joint
[1] 2

$m.out$br_est_joint
        Estimate CI_Lower_Bound CI_Upper_Bound
Break 1        9              5             13
Break 2       38             34             39

$m.out$br_est_joint_time
        Estimate  CI_Lower_Bound CI_Upper_Bound
Break 1 "1985 Q1" "1984 Q1"      "1986 Q1"
Break 2 "1992 Q2" "1991 Q2"      "1992 Q3"

The joint analysis recovers both breaks. The rest of the output contains the coefficient estimates for each quantile followed by the break size estimates at each quantile. At \(\tau = 0.7\):

$m.out$bsize_0.7_Regime_2_minus_Regime_1
                  Value  Std. Error       t value   Pr(>|t|)
Intercept -2.200000e-01 0.086392487 -2.546518e+00 0.01089414
x1         1.000000e-02 0.004616669  2.166064e+00 0.03032812
x2         1.000000e-02 0.015634579  6.396079e-01 0.52244114
x3        -2.674365e-17 0.013993342 -1.911170e-15 1.00000000

$m.out$bsize_0.7_Regime_3_minus_Regime_2
                 Value  Std. Error    t value   Pr(>|t|)
Intercept  0.103333333 0.054260043  1.9044093 0.05688347
x1        -0.008333333 0.003023686 -2.7560181 0.00586080
x2         0.003333333 0.010360097  0.3217473 0.74765036
x3         0.020000000 0.010266885  1.9480105 0.05143960

Note that in $m.out the coefficients and break sizes are reported for every quantile in vec.tau at the \(DQ\) break dates, including \(\tau = 0.80\) and \(0.85\) where the single-quantile analysis found only one break. This is by construction: the joint procedure decides the number and location of the breaks once, using all quantiles together, and then estimates each quantile’s coefficients subject to that common partition.

The two dates are consistent with the 1984 National Minimum Drinking Age Act and a 1991 beer tax increase. The changes are negative and meaningful in magnitude, but smaller for higher quantiles — the policies were more effective for “light drinkers” than for “heavy drinkers”. This is encouraging but falls short of expectations, since heavy drinkers are more likely to cause accidents, suggesting that additional policies are needed to deter heavy drinking.

6.1 A single test at a single quantile

sq.test.0vs1() implements \(SQ_\tau\) of Section 2.4.1 with no search over the number of breaks:

y <- gdp[, "gdp"]
x <- gdp[, c("lag1", "lag2")]

sq.test.0vs1(y, x, v.tau = 0.8, n.size = 1)
#> [1] 2.170968

At the 5% level with three coefficients the critical value is 1.530, so this rejects. dq.test.0vs1() is the corresponding \(DQ\) test over a quantile range:

dq.test.0vs1(y, x, q.L = 0.2, q.R = 0.8, n.size = 1)
#> [1] 1.027587

For a symmetric range its critical values come from the response surface:

res.surface(p = 3, l = 0, q.L = 0.2, q.R = 0.8, d.Sym = TRUE)   # 10%, 5%, 1%
#>           [,1]
#> [1,] 0.8580936
#> [2,] 0.9098714
#> [3,] 1.0223330

The sequential versions of Section 2.4.2 take the break dates estimated under the null and test for one more:

sq.test.lvsl_1(y, x, v.tau = 0.8, n.size = 1, vec.date = 146)
#> [1] 1.064567
dq.test.lvsl_1(y, x, q.L = 0.2, q.R = 0.8, n.size = 1, vec.date = 146)
#> [1] 0.5892746

Both fall below their critical values (1.638 and 0.958 at the 5% level), confirming a single break.

6.2 Estimation given break dates

rq.est.regime() fits each regime separately:

rq.est.regime(y, x, v.tau = 0.8, vec.date = 146, n.size = 1)
#> $Regime_1
#>                Value Std. Error   t value     Pr(>|t|)
#> Intercept  6.1285339 0.58621285 10.454452 0.000000e+00
#> x1         0.3744357 0.07818645  4.789010 4.142365e-06
#> x2        -0.0906033 0.07957956 -1.138525 2.568052e-01
#> 
#> $Regime_2
#>               Value Std. Error  t value     Pr(>|t|)
#> Intercept 3.0392288  0.4220644 7.200865 1.231277e-10
#> x1        0.1636059  0.1187198 1.378084 1.713156e-01
#> x2        0.3141521  0.1122296 2.799192 6.170306e-03

rq.est.full() fits all regimes in a single quantile regression, which is what makes the differences between adjacent regimes and their standard errors available:

fit <- rq.est.full(y, x, v.tau = 0.8, vec.date = 146, n.size = 1)
summary(fit, se = "nid", covariance = TRUE)

ci.date.m() builds the confidence intervals of Section 2.3 for dates supplied by the user; the three columns are the estimate and the lower and upper bounds:

ci.date.m(y, x, vec.tau = 0.8, vec.date = 146, n.size = 1, v.b = 2)
#>      [,1] [,2] [,3]
#> [1,]  146  125  170

Passing a vector to vec.tau gives the multi-quantile interval instead.

6.3 The search itself

The expensive step is gen.long(), which evaluates the objective function over every admissible segment. brdate() then runs the dynamic program of Bai and Perron (2003) over those values, and sq() / dq() apply the sequential tests. Splitting them apart allows the gen.long() output to be reused across several settings of m.max, v.a or norm.method:

T.size    <- length(y)
trim.size <- round(T.size * 0.15)

out.long   <- gen.long(y, x, vec.tau, n.size = 1, trim.size)   # the expensive step
mat.long.s <- out.long$mat.long    # one column per quantile
vec.long.m <- out.long$vec.long    # quantiles combined

mat.date <- brdate(y, x, n.size = 1, m = 3, trim.size, vec.long.m)

dq(y, x, vec.tau, q.L = 0.2, q.R = 0.8, n.size = 1, m.max = 3,
   trim.size, mat.date, d.Sym = TRUE, table.cv = NULL)

d.Sym = TRUE asserts a symmetric quantile range so that the response surface applies, in which case table.cv = NULL is correct. For an asymmetric range set d.Sym = FALSE and supply simulated critical values.