--- title: "QR.break: An R Package for Structural Breaks in Quantile Regression" author: "Zhongjun Qu, Tatsushi Oka, and Samuel Messer" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{QR.break: An R Package for Structural Breaks in Quantile Regression} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 3.6, fig.align = "center" ) options(warn = -1) # the quantile regression fit is occasionally non-unique ``` > This vignette follows Qu, Oka and Messer, *"QR.break: An R Package for > Structural Breaks in Quantile Regression"*, which should be consulted for > further detail. It outlines the methods implemented in the package and then > works through the two datasets: a time series on U.S. real GDP > growth and a repeated cross-section on youth drinking and driving. All numerical > output shown here was produced with version `r utils::packageVersion("QR.break")` > of the package. ## 1. Introduction The issue of structural breaks has been extensively studied. Much of the literature has focused on the conditional mean, but in many cases structural change in the conditional quantile function is more informative. For example, in studying income inequality, it is important to assess whether the wage gap between racial groups, conditional on covariates, has changed over time. Inequality may increase conditional dispersion without affecting the mean. Similarly, when evaluating a policy aimed at improving outcomes for low-performing students, attention should focus on lower quantiles. In both cases it is desirable to estimate break dates from the data: in the former, the cause of change may be unclear a priori; in the latter, the policy effect may exhibit an unknown time lag. To address these issues, Qu (2008) and Su and Xiao (2008) developed tests for detecting structural change in conditional quantile functions with unknown break dates, but did not consider estimation and inference for the number of breaks or their locations. Oka and Qu (2011) studied the estimation of multiple structural breaks at unknown dates in conditional quantile functions for two models: a time series model, useful for analyzing macroeconomic data, and a repeated cross-section model, relevant for evaluating social programs and policy effects. Their framework allows for structural change in single or multiple quantiles. Analyzing multiple quantiles requires stronger assumptions but can improve estimation efficiency. Key aspects of Oka and Qu's (2011) procedure are as follows. Assuming a known number of breaks, the methods construct estimates of break dates and coefficients as global minimizers of the check function over admissible break points. For multiple quantiles, the check function is integrated over the set of quantiles chosen by the user. The assumptions permit dynamic models and impose restrictions only in neighborhoods around the quantiles of interest, leaving other quantiles unspecified. This flexibility allows researchers to examine slices of the conditional distribution without imposing global distributional assumptions. The distributions of the break estimators were derived following Picard (1985) and Yao (1987); they involve consistently estimable parameters, which enables confidence interval construction without simulation. Oka and Qu (2011) also propose a test for the number of breaks based on the subgradient methods of Qu (2008). These tests do not require variance estimation and have monotonic power even with multiple breaks. The package's main function, `rq.break()`, returns all testing and estimation results based on user-specified quantiles of interest, the maximum number of allowed breaks, and the minimum length of a single regime: ```r result <- rq.break(y, x, vec.tau, N, trim.e, vec.time, m.max, v.a, v.b, verbose, norm.method) ``` We illustrate the function through the two empirical examples of Oka and Qu (2011). The first revisits the "Great Moderation" in U.S. GDP growth using quarterly data; the results suggest the decline in volatility mainly affected the upper tail, with the median and lower quantiles remaining stable, implying that expansions became less rapid while recessions remained as severe. The second analyzes blood alcohol levels of young drivers in California from 1983 to 2007; two breaks are detected, consistent with the 1984 National Minimum Drinking Age Act and a 1991 beer tax increase, with effects more pronounced at lower quantiles, indicating a greater impact on lighter drinkers than on heavier ones. ```{r load} library(QR.break) ``` ## 2. Methods for estimation and inference This section explains (1) the model and the econometric issues of interest; (2) estimating break locations and regression coefficients when the number of breaks is known, based on a single quantile; (3) the same based on multiple quantiles; (4) a procedure for determining the number of breaks; and (5) inference on quantile regression coefficients after estimating the breaks. Section 2.6 covers the normalization option introduced in package version 1.0.3. ### 2.1 Econometric models and issues of interest #### 2.1.1 The model Let $y_{it}$ be the outcome for unit $i$ in period $t$, $x_{it}$ a $p$-dimensional vector of covariates, and $Q_{y_{it}}(\tau \mid x_{it})$ the conditional $\tau$-quantile of $y_{it}$ given $x_{it}$. Let $N$ be the number of cross-sectional units, assumed constant over time, and $T$ the number of periods. For a single time series $N = 1$ and the index $i$ is dropped throughout. Assume the conditional quantile function is linear and potentially affected by $m$ structural breaks: $$ Q_{y_{it}}(\tau \mid x_{it})= \begin{cases} x_{it}^{\prime }\beta _{1}^{0}(\tau ), & t=1,\dots ,T_{1}^{0}, \\ x_{it}^{\prime }\beta _{2}^{0}(\tau ), & t=T_{1}^{0}+1,\dots ,T_{2}^{0}, \\ \qquad \vdots & \qquad \vdots \\ x_{it}^{\prime }\beta _{m+1}^{0}(\tau ), & t=T_{m}^{0}+1,\dots ,T, \end{cases} $$ where $\tau \in (0,1)$, the $\beta_j^0(\tau)$ $(j = 1,\dots,m+1)$ are unknown parameters, and the $T_j^0$ $(j = 1,\dots,m)$ are unknown break dates. The regressors may include discrete as well as continuous variables. **A column of ones is automatically added to the regression when applying the methods.** #### 2.1.2 Examples The following examples, taken from Oka and Qu (2011), illustrate the model. **Example 1 (time series).** Oka and Qu (2011) studied a quantile autoregression for U.S. quarterly real GDP growth rates over 1947:Q2 to 2009:Q2, in which the regressors are lagged values of the growth rate and the intercept and autoregressive coefficients are allowed to shift across regimes. They detect a structural break in 1984 that affects only the upper quantiles of the distribution. The coefficient estimates suggest that growth was slower during expansions, while recessions remained just as severe when they occurred. This dataset is included in the package and is used in Section 4. **Example 2 (repeated cross-sections).** Motor vehicle crashes are the leading cause of death among youth aged 15–20, a high proportion of which involve drunk driving. Blood alcohol concentration (BAC) is a key measure of alcohol impairment, and changes in BAC among young drivers provide useful information on how their drinking behavior has changed over time. Oka and Qu (2011) studied structural change in BAC among young drivers involved in traffic accidents, using age, gender, and a dummy for the fourth quarter as regressors. They detect breaks in 1985 and 1992. The changes are negative and meaningful in magnitude; however, the change is smaller for higher quantiles, suggesting the policies are more effective for "light drinkers" than for "heavy drinkers" in the sample. This is encouraging but falls short of expectations, as heavy drinkers are more likely to cause accidents, suggesting that additional policies are needed to deter heavy drinking. This dataset is included in the package and is used in Section 5. #### 2.1.3 Econometric issues of interest The methods address the following issues. 1. **Estimation based on a single quantile when the number of breaks is known.** The method estimates both the break locations and the regression coefficients. If the user specifies more than one quantile level, the analysis is performed independently for each quantile, allowing break locations to differ across quantiles. The program returns the estimated break locations, their confidence intervals, and the corresponding estimates and intervals for the regression coefficients. 2. **Estimation based on multiple quantiles when the number of breaks is known.** Here the break locations are assumed common across quantiles and are estimated using information from all specified quantiles. The program returns the same set of objects. 3. **Selection of the number of breaks.** The user specifies the maximum number of breaks, and the program determines the number of breaks using a dynamic programming algorithm. A single main function performs all of these tasks. ### 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.4.1 Testing for a single structural break The $SQ_\tau$ test detects the presence of a structural break in a given quantile $\tau$: $$ SQ_{\tau }=\sup_{\lambda \in \lbrack 0,1]}\left\Vert (\tau (1-\tau ))^{-1/2}\left[ H_{\lambda ,T}(\hat{\beta}(\tau ))-\lambda H_{1,T}(\hat{\beta}% (\tau ))\right] \right\Vert _{\infty }, $$ where $$ H_{\lambda ,T}(\hat{\beta}(\tau ))=\left( \sum_{t=1}^{T}\sum_{i=1}^{N}x_{it}x_{it}^{\prime }\right) ^{-1/2}\sum_{t=1}^{% [\lambda T]}\sum_{i=1}^{N}x_{it}\psi _{\tau }(y_{it}-x_{it}^{\prime }\hat{\beta}% (\tau )), $$ with $\psi_\tau(u) = \tau - 1(u < 0)$ and, as before, the summation over $i$ absent for a single time series. Here $\hat\beta(\tau)$ is the full-sample estimate assuming no structural change, and $\|\cdot\|_\infty$ is the sup norm, revealing the strongest evidence against the null. The $DQ$ test detects structural changes in quantiles in an interval $\mathcal{T}_\omega$, taking a further supremum over $\tau$: $$ DQ=\sup_{\tau \in \mathcal{T}_{\omega }}\sup_{\lambda \in \lbrack 0,1]}\left\Vert H_{\lambda ,T}(\hat{\beta}(\tau ))-\lambda H_{1,T}(\hat{\beta}% (\tau ))\right\Vert _{\infty }. $$ These correspond to the exported functions `sq.test.0vs1()` and `dq.test.0vs1()`. #### 2.4.2 Testing $l$ against $l+1$ breaks Suppose a model with $l$ breaks has been estimated, with break estimates $\hat{T}_1,\dots,\hat{T}_l$. These partition the sample into $l+1$ segments, the $j$th being $[\hat{T}_{j-1}+1, \hat{T}_j]$. The strategy tests each of the $l+1$ segments for the presence of an additional break. Letting $SQ_{\tau,j}$ and $DQ_j$ denote the tests applied to the $j$th segment, $$ SQ_{\tau }(l+1\mid l)=\max_{1\leq j\leq l+1}SQ_{\tau ,j},\qquad DQ(l+1\mid l)=\max_{1\leq j\leq l+1}DQ_{j}, $$ where the segment statistics replace the full-sample sums by sums running from $T_{j-1}+1$ to $T_j$, and $\hat\beta_j(\tau)$ is estimated using the $j$th regime only. We reject the null in favor of a model with $l+1$ breaks if the resulting value exceeds the corresponding critical value. These correspond to `sq.test.lvsl_1()` and `dq.test.lvsl_1()`. #### 2.4.3 Critical values These tests are asymptotically nuisance parameter free, and tables of critical values are provided in Qu (2008). They do not require the estimation of any variance parameter, hence have monotonic power even when multiple breaks are present. The package has built-in critical values for common configurations: for the $SQ$ test when the number of coefficients is below 100, and for the $DQ$ test when the number of coefficients is below 20 and the quantile trimming is symmetric ($\omega_2 = 1 - \omega_1$), via a response surface (`res.surface()`). Outside those cases — asymmetric trimming, more than 20 coefficients, or `m.max` above 5 — the $DQ$ critical values are computed by simulation, which can add a few minutes or more to the running time. #### 2.4.4 The recommended procedure Consider the interval $\mathcal{T}_\omega$ and the quantile grid $\tau_1,\dots,\tau_q \in \mathcal{T}_\omega$. * **Step 1.** Apply the $DQ$ test. If it does not reject, conclude there is no break and terminate. If it rejects, estimate the model allowing one break, save the estimated break date, and proceed. * **Step 2.** Apply the $DQ(l+1 \mid l)$ tests starting with $l = 1$. Increase $l$ if the test rejects. At each stage the model is re-estimated and the break dates are global minimizers of the objective function allowing $l$ breaks. Continue until the test fails to reject. * **Step 3.** Let $\hat{l}$ denote the first value for which the test fails to reject. Estimate the model allowing $\hat{l}$ breaks. Save the estimated break dates and confidence intervals. * **Step 4.** Treat the $q$ quantiles separately: for every $\tau_h$, apply the $SQ_\tau$ and $SQ_\tau(l+1 \mid l)$ tests, carrying out Steps 1–3. Examine whether the estimated breaks agree with those from Step 3. `rq.break()` prints the single-quantile analysis first and the joint analysis second, which is the order of computation rather than the order of reading. When interpreting output, consider follow the steps above. ### 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. Using the main function ### 3.1 Arguments | Argument | Meaning | |---|---| | `y` | The dependent variable, a vector of length $NT$ | | `x` | The regressors, an $NT \times p$ matrix or data frame. A column of ones should **not** be included; it is added automatically | | `vec.tau` | The quantiles of interest, e.g. `seq(0.1, 0.9, by = 0.1)`. A coarse grid is generally sufficient | | `N` | The size of the cross-section; set to 1 for time series data | | `trim.e` | The minimum length of any regime as a fraction of the total time span. **Values between 0.10 and 0.20 are recommended** | | `vec.time` | The time indices, of length $T$, used for reporting the estimated break dates; or `NULL` to report indices | | `m.max` | The maximum number of breaks allowed (at most 10) | | `v.a` | Significance level for determining the number of breaks: 1, 2 or 3 for 10%, 5% or 1% | | `v.b` | Coverage level for the confidence intervals of break dates: 1 or 2 for 90% and 95% | | `verbose` | `TRUE` displays results in the console; the default is `FALSE` | | `norm.method` | `"cholesky"` (default) or `"spectral"`; see Section 2.6 | 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: | Element | Contents | |---|---| | `test_`, `test_joint` | Test statistics and critical values | | `nbreak_`, `nbreak_joint` | Number of breaks detected | | `br_est_`, `br_est_joint` | Break dates and confidence intervals, as indices | | `br_est_time_`, `br_est_joint_time` | The same in calendar format | | `coef_` | Regime-by-regime coefficient estimates | | `bsize__Regime_j+1_minus_Regime_j` | Break sizes between adjacent regimes | ## 4. Application 1: a time series — U.S. real GDP growth ### 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. ```{r gdp-data} data(gdp) str(gdp) head(gdp, 3) ``` ### 4.2 Setting the inputs The following commands define `y` and `x`: ```{r gdp-setup} 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: ```{r gdp-tau} vec.tau <- seq(0.20, 0.80, by = 0.15) vec.tau ``` 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 `r round(0.15 * 247)` 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 ```{r gdp-run, eval = FALSE} 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: | $\tau$ | 0.20 | 0.35 | 0.50 | 0.65 | 0.80 | |---|---|---|---|---|---| | $SQ_\tau$ (1 break) | 1.423 | 1.479 | 1.001 | 1.818 | 2.171 | | 5% critical value | 1.530 | 1.530 | 1.530 | 1.530 | 1.530 | | Breaks detected | 0 | 0 | 0 | 1 | 1 | 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 ```{r gdp-plot} 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: | $\tau$ | $SQ_\tau$, `"cholesky"` | $SQ_\tau$, `"spectral"` | 5% critical value | |---:|---:|---:|---:| | 0.20 | 1.423 | 1.158 | 1.530 | | 0.35 | 1.479 | 1.461 | 1.530 | | 0.50 | 1.001 | 0.962 | 1.530 | | 0.65 | 1.818 | 1.845 | 1.530 | | 0.80 | 2.171 | 2.332 | 1.530 | | $DQ$ | 1.028 | 1.043 | 0.910 | 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: ```{r err1, error = TRUE} 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) ``` As another example, if `trim.e` is too small the message suggests increasing it: ```{r err2, error = TRUE} 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) ``` When an error message is produced the program exits with no saved results, and the user can modify the inputs and restart. ## 5. Application 2: repeated cross-sections — youth drinking and driving ### 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. ```{r driver-data} data(driver) str(driver) ``` ```{r driver-setup} y <- driver[, "bac"] x <- driver[, c("age", "gender", "winter")] vec.time <- unique(driver[, "yq"]) # length T = 100, one label per quarter length(vec.time) ``` The choice of quantiles is dictated by the data. BAC is exactly zero for most of these drivers: ```{r driver-zeros} mean(driver$bac == 0) quantile(driver$bac, c(0.50, 0.60, 0.65, 0.70, 0.80, 0.85)) ``` 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 ```{r driver-run, eval = FALSE} 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: | $\tau$ | 0.70 | 0.75 | 0.80 | 0.85 | |---|---|---|---|---| | $SQ_\tau$ (1 break) | 5.179 | 4.613 | 3.503 | 3.258 | | $SQ_\tau(2 \mid 1)$ | 2.205 | 1.931 | 1.392 | 1.056 | | 5% critical values | 1.575 / 1.679 | 1.575 / 1.679 | 1.575 / 1.679 | 1.575 / 1.679 | | Breaks detected | 2 | 2 | 1 | 1 | 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. Working with the individual functions `rq.break()` is a wrapper around exported components, each of which can be used on its own — to test at a date decided in advance, to reuse an expensive computation, or to inspect intermediate quantities. ### 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: ```{r sq-test} y <- gdp[, "gdp"] x <- gdp[, c("lag1", "lag2")] sq.test.0vs1(y, x, v.tau = 0.8, n.size = 1) ``` 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: ```{r dq-test} dq.test.0vs1(y, x, q.L = 0.2, q.R = 0.8, n.size = 1) ``` For a symmetric range its critical values come from the response surface: ```{r res-surface} res.surface(p = 3, l = 0, q.L = 0.2, q.R = 0.8, d.Sym = TRUE) # 10%, 5%, 1% ``` The sequential versions of Section 2.4.2 take the break dates estimated under the null and test for one more: ```{r seq-tests} sq.test.lvsl_1(y, x, v.tau = 0.8, n.size = 1, vec.date = 146) dq.test.lvsl_1(y, x, q.L = 0.2, q.R = 0.8, n.size = 1, vec.date = 146) ``` 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: ```{r est-regime} rq.est.regime(y, x, v.tau = 0.8, vec.date = 146, n.size = 1) ``` `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: ```r 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: ```{r ci-date} ci.date.m(y, x, vec.tau = 0.8, vec.date = 146, n.size = 1, v.b = 2) ``` 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`: ```r 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. ## 7. Practical notes **Runtime.** The cost is dominated by `gen.long()`, which is $O(T^2)$ per quantile and does not depend on `m.max`. Raising `m.max` adds only the cheap search over precomputed values and a few more sequential tests, so there is little to gain by keeping it artificially small. What drives the cost is $T$ (quadratically) and the number of quantiles (linearly). **Warnings about non-unique fits.** `rq()` warns when a quantile regression solution is not unique, which happens routinely on short segments and on data with many tied values — the `driver` BAC data especially. These warnings are benign here; `options(warn = -1)` suppresses them. **"confidence interval is out of the range".** This warning means an interval extends past the start or end of the sample. The date estimate is still reported, but is not converted to calendar form, so `br_est_time_*` will be missing for that quantile. It usually signals a break estimated close to a sample boundary, or one estimated imprecisely. **Reproducibility.** When the $DQ$ critical values come from tables or the response surface they are fixed and repeated runs give identical output. When they are simulated instead, they carry Monte Carlo noise, and the package does not set a seed internally. In the `driver` example the simulated 5% critical values move by roughly $\pm 0.003$ across runs — far too small to change the number of breaks detected, but enough that the printed critical values will not match to the last digit from one session to the next. Call `set.seed()` before `rq.break()` if byte-identical output is required. Everything else — break dates, coefficients, and the $SQ$ and $DQ$ statistics themselves — is deterministic. ## References Bai, J. (1995). Least Absolute Deviation Estimation of a Shift. *Econometric Theory*, 11, 403–436. Bai, J. and P. Perron (2003). Computation and Analysis of Multiple Structural Change Models. *Journal of Applied Econometrics*, 18, 1–22. Kessy, A., A. Lewin and K. Strimmer (2018). Optimal Whitening and Decorrelation. *The American Statistician*, 72(4), 309–314. Koenker, R. (2005). *Quantile Regression*. Cambridge University Press. Oka, T. and Z. Qu (2011). Estimating Structural Changes in Regression Quantiles. *Journal of Econometrics*, 162, 248–267. [doi:10.1016/j.jeconom.2011.01.005](https://doi.org/10.1016/j.jeconom.2011.01.005) Picard, D. (1985). Testing and Estimating Change-Points in Time Series. *Advances in Applied Probability*, 17, 841–867. Qu, Z. (2008). Testing for Structural Change in Regression Quantiles. *Journal of Econometrics*, 146, 170–184. [doi:10.1016/j.jeconom.2008.08.006](https://doi.org/10.1016/j.jeconom.2008.08.006) Su, L. and Z. Xiao (2008). Testing for Parameter Stability in Quantile Regression Models. *Statistics and Probability Letters*, 78, 2768–2775. Yao, Y. C. (1987). Approximating the Distribution of the Maximum Likelihood Estimate of the Change-Point in a Sequence of Independent Random Variables. *The Annals of Statistics*, 15, 1321–1328.