Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
67,809 characters · 14 sections · 66 citation commands
Partitioned Wild Bootstrap for Panel Data Quantile Regression
\if00 {
} \fi
\if10 {
} \fi
\doublespacing
Since the seminal work of KoenkerBassett78, quantile regression (QR) models have provided a valuable tool as a way of capturing heterogeneous effects that covariates may exert on the outcome of interest, exposing a wide variety of forms of conditional heterogeneity under weak distributional assumptions. QR methods have been widely employed to estimate causal effects as well as structural economic models.
Recently, there has been growing empirical and theoretical literatures on QR models for panel data. Koenker04 introduced a general approach for estimation of fixed effects quantile regression (FE-QR) models which treats the individual effects as parameters to be estimated. The FE-QR estimator is designed to control for individual-specific heterogeneity while exploring heterogeneous covariate effects, and therefore, it provides a flexible method for practical analysis of panel data. Consistency and weak convergence of the FE-QR estimator have been established in Kato2012 and GalvaoGuVolgushev20 under conditions that are similar to the ones considered for nonlinear panel data estimators. Many papers have suggested alternative methods for panel data QR and derived the corresponding statistical properties, including, among others, Lamarche10, Canay11, KimYang11, GalvaoLamarcheLima13, ArellanoBonhomme15, GrahamHahnPoirierPowell2015, GalvaoKato16, ChetverikovLarsenPalmer16, GuVolgushev19, MachadoSantosSilva19, and CHERNOZHUKOV2024. A number of these studies rely on large-$T$ panel data, where the number of time periods ($T$) grows faster than the number of cross-sectional units ($N$).
Practical inference procedures for FE-QR have been a pervasive concern in empirical work, especially when the time-series is large and temporal dependence has to be accommodated. Inference methods, which have been developed mostly for independent data, have mainly used asymptotic approximations for the construction of test statistics and the variances of the estimators. Nevertheless, these variances depend on the conditional density of the innovations, which can be difficult to work with in practice and have bandwidth nuisance parameters. Hence, the use of the bootstrap as an alternative to such asymptotic approximations has been considered, but its properties have not received the same amount of attention as in the QR for cross-section or time-series literatures. The sole exceptions are LamarcheParker23 and GalvaoParkerXiao24. LamarcheParker23 proposed a copula-based approach for a penalized estimator under stationary and $\beta$-mixing conditions, but the implementation of the approach relies on parametric copulas and knowledge of dependence parameters. GalvaoParkerXiao24 propose a bootstrap based on multiplying one non-negative weight with all of a unit's observations. While this approach has the advantage that serial dependence within each cross-sectional unit is preserved, it can lead to size distortions when practitioners use panels with a small number of cross-sectional units. The bootstrap literature on the time-series dimension for conditional average linear panel data models has a longer history Kapetanios08,Goncalves11,GoncalvesKaffo15. The analysis of bootstrap algorithms for the linear mean regression model is aided by its linearity and the within transform. Meanwhile, the model for conditional quantiles with individual effects must be treated as a nonlinear model, and this makes bootstrap inference a challenging problem.
This paper develops a new bootstrap procedure for FE-QR models that is easy to implement. We propose a novel wild residual bootstrap for dependent panel data where one bootstrap weight is applied to each cell in a partition of a unit's observations. Intuitively, reweighting residuals is particularly useful to accommodate conditioning on covariates and fixed effects, while the partition-based approach is aimed at estimating the stationary distribution over a large number of time periods. Hence, we incorporate dependence into bootstrap procedures for panel data FE-QR, improving upon existing results in the literature. This is an important addition, since the statistical properties of panel data QR models are established under large panel approximations. Our results contribute to the broader QR literature on cross-section and time-series models as in Hahn95,Hahn97, Fitzenberger98, Shao10ss, FengHeHu11, SongRitovHardle12, Hagemann17, and GregoryLahiriNordman18.
We establish the asymptotic validity of the procedure using developments that are different than those used in LamarcheParker23, since new results are needed in the dependent case. While the conditions for the marginal weight distribution is standard in the literature FengHeHu11, we require a partition cell length condition which is different from the block length conditions considered by Fitzenberger98, Shao10ss and GregoryLahiriNordman18. We show how to empirically choose cell length, which aims to minimize the difference between the bootstrap variance and the asymptotic variance of the FE-QR estimator. Using an extensive simulation study, we find that empirical coverage rates for the proposed bootstrap are close to the nominal counterparts and coverage improves as the sample size increases.
The remainder of the paper is organized as follows. Section (ref) reviews standard inference procedures for FE-QR. In Section (ref), we present the new bootstrap procedure and discuss its uses for inference. Section (ref) evaluates the finite sample performance of the bootstrap procedure. In Section (ref), we estimate a panel data quantile regression model and apply the method to evaluate how consumers respond to time-of-use electricity pricing. Finally, Section (ref) concludes.
Our model of interest is the fixed effects quantile regression (FE-QR) model. We observe $T$ time periods (indexed by $t$) of jointly stationary data $\{ ( y_{it},\bm{x}_{it}') \}_{t=1}^{T}$ for each of $N$ units (indexed by $i$), where $y_{it} \in \mathbb{R}$ denotes the response and $\bm{x}_{it}$ denotes a $p$-dimensional vector of covariates for unit $i$ at time $t$. For quantile $\tau \in (0,1)$, the FE-QR model is
where $u_{it}$ is a disturbance whose $\tau$-th quantile conditional on $\bm{x}_{it}$ is equal to zero, implying that the conditional quantile of the response variable is $Q_{y_{it}}(\tau | \bm{x}_{it}) = \bm{x}_{it}' \bm{\beta}_0(\tau) + \alpha_{i0}(\tau)$. The parameter of interest is $\bm{\beta}_0(\tau) \in \mathcal{B} \subseteq \mathbb{R}^p$, and the scalar individual specific effect, $\alpha_{i0}(\tau)$, is treated as a nuisance parameter. Because we consider just one quantile value, we suppress the dependence of the parameters on $\tau$ in the sequel.
Let $\bm{\theta} = (\bm{\beta}',\bm{\alpha}')' \in \bm{\Theta} \subseteq \mathbb{R}^{p+N}$, where $\bm{\alpha} = (\alpha_{1},...,\alpha_{N})'$, and let $\bm{\theta}_0 = (\bm{\beta}_0',\bm{\alpha}_0')'$. To estimate $\bm{\theta}_0$, we consider the following FE-QR estimator:
where $\rho _{\tau }(u)=$ $u(\tau -I(u<0))$ is the quantile regression loss function. The estimator defined in (ref) was introduced by Koenker04. Although the parameter of interest is $\bm{\beta}_0$, it is important to note that this estimation strategy estimates all of $\bm{\theta}_0$ and then focuses on $\bm{\beta}_0$ for inference purposes. This is due to the fact that no known data transformation exists that allows one to avoid estimation of $\bm{\alpha}$ for this model as in, for example, linear conditional average panel models. It prevents us from naively applying bootstrap algorithms from the time-series literature unit-wise to these observations, as is shown below in Section (ref).
The asymptotic distribution of $\hat{\bm{\beta}}$ is described below in Lemma (ref). It depends on the following regularity conditions.
Similar conditions are used in the literature. For instance, a version of Assumption (ref) has been used in Kato2012 and GalvaoGuVolgushev20. Assumptions (ref)--(ref) guarantee the convexity of the limiting distribution of the objective function, and therefore, the uniqueness of the fixed effects estimator (see Kato2012 for slightly more minimal assumptions for consistency). Assumption (ref) implies appropriate moment conditions on the covariates and it is similar to (B1) in Kato2012 and (B3) in LamarcheParker23. Lastly, (ref) assumes the extra regularity conditions needed, beyond those for consistency, to establish an asymptotic distribution for $\hat{\bm{\beta}}$. The first condition is technical and used to derive bounds on special functional classes used in the proof of asymptotic normality. The other conditions are for convenience and ensure existence of limiting covariance matrices in the case of dependent data, and are similar to (C3) in LamarcheParker23.
The large sample theory of the estimator (ref) is described in the following result.
In Lemma (ref), $T$ diverges to infinity to show the asymptotic results, and it must grow faster than $N$ for asymptotic normality. The condition on the relative sizes of $T$ and $N$ is close to the standard rates for smooth non-linear panel data models, despite the lack of differentiability of the quantile regression score function $\psi_\tau(u) = \tau - I(u < 0)$.
We now offer some heuristics that help with intuition on the limitations of existing inference approaches applied to panel data quantile regression, as well as to establish the validity of the proposed bootstrap approach in Section (ref) below. Our results rely on the Bahadur representation of the FE-QR estimator. Kato2012 established conditions under which the Bahadur representation of the estimator in equation (ref) holds, specifically
where $\tilde{\bm{x}}_{it}$, and $\bm{D}_N$ are defined in Assumption (ref). The other part of the variance that was defined in that assumption was labeled $\bm{V}_{NT}$ and is the variance of the remaining part of the linear term in (ref). With temporal dependence of unit $i$'s observations the variance of this linear term includes a weighted sum of covariances
Note that $\gamma_{i}(0) = \tau (1 - \tau) \textnormal{E} \left[ \tilde{\bm{x}}_{it} \tilde{\bm{x}}_{it}' \right]$. Then we can write $\bm{V}_{NT}$ as
We are interested in conducting inference for the parameter $\bm{\beta}_{0}$. Most straightforwardly, one could estimate the asymptotic variance-covariance matrix $\bm{\Sigma}$ in (ref) and construct confidence intervals directly from it. However, in the presence of serial dependence, a heteroskedasticity autocorrelation consistent (HAC) estimator of $\bm{V}_{NT}$ is essential, and at the moment, theoretical guidance about such an estimator is lacking in the literature. In addition, $\bm{D}_N$ depends on the conditional density of the error term, which is generally difficult to estimate well. Hence, the next section proposes a bootstrap procedure for inference on $\bm{\beta}_0$ in the FE-QR model under dependence.
Before we provide details on the bootstrap procedure, we introduce a required condition on the data generating process to establish the validity of the bootstrap described in Section (ref) below. Let $\bm{V}_{NT}^0$ describe the variance of the sum in the Bahadur representation (ref), that is, the first term in equation (ref). In the case where errors at different time periods for each unit are independent, we can write (ref) as:
Longitudinal data in the social sciences including economics commonly exhibit positive time-series dependence that makes the variance in (ref) larger than $\bm{V}_{NT}^0$. In this paper, we will propose a bootstrap that is tailored to cases where
converges to a positive definite limit. This regularity condition is stated formally as Assumption (ref) below. To make the results more intuitive, we use the convention that for two conformable matrices $\bm{A}$ and $\bm{B}$, $\bm{A} \geq \bm{B}$ means that $\bm{A} - \bm{B}$ is positive semidefinite. Consider the following condition:
Assumption (ref) imposes two conditions on the data for which the proposed bootstrap algorithm works well. The first is on the sum of the autocovariances in (ref), which allows for positive or negative dependence between different time periods. This condition together with the stationarity and mixing condition on the data in Assumption (ref) imply that the covariance between terms $\tilde{\bm{x}}_{it} \psi_\tau(u_{it})$ and $\tilde{\bm{x}}_{it+k} \psi_\tau(u_{it+k})$ can be bounded, due to the boundedness of the function $\psi_\tau(\cdot)$ that appears in the Bahadur representation (ref). As noted above, $\gamma_i(0)$ has a convenient familiar form. The covariance between periods $t$ and $t+k$ can be similarly bounded. Writing
the mixing condition in Assumption (ref) implies the bound $|\textnormal{P} \left\{ u_{it} < 0, u_{it+k} < 0 \, | \, \bm{x}_i \right\} - \tau^2| \rightarrow 0$, which means that it is less than $\tau(1 - \tau)$ for all $k$ large enough. This fact will be exploited by the partitioned wild bootstrap, which is described in detail in the next section. The second condition in Assumption (ref) limits the dependence between the regressors $\tilde{\bm{x}}_{it}$ and the quantile error terms. Without such a condition, oscillating patterns that offset each other can cause problems with partition selection for the bootstrap.
This section proposes a new wild residual bootstrap approach. For each unit $i$, we partition the time-series observations into $b$ non-overlapping parts of length $l$. We denote blocks of time series observations as partitions (or cells) instead of blocks to distinguish our approach from related methods which resample blocks Fitzenberger98,Goncalves11,GregoryLahiriNordman18,Hounyo23. In contrast, we resample partition invariant weights to construct a bootstrap distribution. The length $l$ of the partition is selected by a data-driven approach introduced below.
For simplicity, we assume that $T = b \times l$. Letting $j = 1, 2, \ldots, b$, and denoting within-partition observations by $s = 1, 2, \ldots, l$, we redefine the FE-QR estimator in (ref):
The corresponding residuals are $\hat{u}_{ijs} = y_{ijs} - \bm{x}_{ijs}' \hat{\bm{\beta}} - \hat{\alpha}_i$. New bootstrap observations are obtained from predicted values and reweighted residuals. Specifically, the bootstrap data generating process relies on the following bootstrap residuals
where the weight $w_{ij}$ is drawn from a pre-specified distribution satisfying conditions (ref)--(ref) below. These within-partition invariant weights are independent and identically distributed (i.i.d.). Using the bootstrap residuals defined in (ref), the bootstrapped dependent variable is
Finally, a bootstrap FE-QR estimate is computed by finding
This procedure is labeled partitioned wild bootstrap (PWB). Theorem (ref) below establishes the asymptotic validity of the method and Section (ref) describes how one can obtain valid confidence intervals using the bootstrap estimator proposed in (ref).
The distribution of the random weights is chosen by the researcher. As mentioned above, they are i.i.d., which makes them easy to generate, but must come from a distribution with CDF $G_W$ that satisfies the following conditions.
Assumption (ref) is used in similar versions of the wild bootstrap with the exception that we require here that the weights are independent across partitions LamarcheParker23. Assumptions (ref)--(ref) on $G_W$ were first proposed in FengHeHu11 and they are satisfied by several weight distributions FengHeHu11, lanWang2018. We follow the literature and adopt a two-point mass distribution in the empirical examples. The distribution generates $w_{ij} = 2 (1-\tau)$ with probability $\tau$ and $w_{ij} = - 2 \tau$ with probability $(1-\tau)$.
The practical implementation of the PWB method is simple. The main algorithm for implementing the methods is as follows.
By choosing the number of bootstrap simulations $B$ in the algorithm above large enough, the distribution of $\sqrt{NT}(\bm{\beta}^{\ast} - \hat{\bm{\beta}})$ can be computed with any desired precision. There are several way of using this distribution for inference on the parameters.
The distribution function of $\bm{\beta}^\ast - \hat{\bm{\beta}}$ can be used to estimate the distribution function of $\hat{\bm{\beta}} - \bm{\beta}_0$. Specifically, suppose that $\beta_0$ is one coordinate of $\bm{\beta}_0$ and we would like to compute a confidence interval for $\beta$ with confidence level $1-\lambda$. Given bootstrap realizations $\{\beta_b^\ast\}_{b=1}^B$ of this coordinate, we can find $\beta_{\lambda/2}^\ast$ and $\beta_{1-\lambda/2}^\ast$. Then one confidence interval for $\beta$ is
These percentiles may be used as to estimate the endpoints of a confidence interval for $\beta_0$.
For a fixed quantile level $\tau$, we define the bootstrap estimate of the asymptotic covariance matrix $\bm{\Sigma}$ given bootstrap realizations $\{\bm{\beta}^{\ast }_b\}_{b=1}^B$ as
Under the regularity conditions in Theorem (ref) below, $\bm{\Sigma}^\ast$ is a consistent estimator of the asymptotic covariance matrix $\bm{\Sigma}$ defined in equation (ref). The estimated standard errors of $\hat{\bm{\beta}}$ are the square roots of the diagonal elements of $\bm{\Sigma}^{\ast}$. Given $\bm{\Sigma}^\ast$, testing general hypotheses $R\bm{\beta}_0=r$ for the vector $\bm{\beta}_0$ can be accommodated by Wald-type tests.
In one dimension we can compare this method with the percentile method. Once again, assume that $\beta_0$ is one coordinate of $\bm{\beta}_0$. Let $\hat{\beta}$ be its estimate and let $se^\ast$ be the square root of the corresponding diagonal element of $\bm{\Sigma}^\ast$. Then a standard-error based confidence interval is
where $z_{\lambda}$ denotes the $\lambda$-th quantile of the standard normal distribution.
The selection of the size of the partition is different for the PWB algorithm than for traditional block bootstrap approaches Fitzenberger98, Shao10ss, GregoryLahiriNordman18. The length here is chosen to match the central term $\bm{V}$ in the asymptotic covariance matrix in equation (ref), in particular the average of the covariances between non-contemporaneous score terms in the Bahadur representation (ref), and the bootstrap variance. The next result indicates how to choose the size of these cells, assuming that Assumptions (ref)--(ref) hold.
The partition length $l^{\text{o}}$ in (ref) can be estimated by numerically solving a finite-sample analog using plug-in estimates. Let $\check{\bm{x}}_{it} = \bm{x}_{it} - \bar{\varphi}_i^{-1} \bar{\bm{g}}_i$, where $\bar{\varphi}_i$ and $\bar{\bm{g}}_i$ are respectively the sample analogs of $\varphi_i$ and $\bm{g}_i$. Then we find $\hat{l}$ such that
where $K(\cdot)$ is a kernel function for variance estimation with bandwidth $h$ to estimate each unit's contribution to the variance GalvaoYoon24. Without kernel weighting, the right-hand side of (ref) would be almost numerically zero -- these sums are sample analogs of the sum of $\gamma_i(k)$ terms for $k > 0$ on the right-hand side of (ref), but they also represent half of the off-diagonal terms of the first derivative of the empirical loss function evaluated at the optimizer, squared. To implement the selection rule, it is convenient to solve for the partition size unit by unit over the panel, leading to $\{\hat{l}_i\}_{i=1}^N$, where for each $i$,
Under the conditions of Lemma (ref), the feasible length selection $\hat{l}$ in (ref) converges in probability to $l^\text{o}$, as stated in the following result.
Length selection based on (ref) is investigated in the simulation study reported in Section (ref). However, to help with the intuition on how the procedure works, we now offer an illustrative example. If we ignore the contributions of $\check{\bm{x}}_{it}$, $\hat{l}$ can be found to satisfy
This approximate equality can be solved for $\hat{l}$, and, recalling that it should be a positive integer, we suggest
where for $a \in \mathbb{R}$, $\lceil a \rceil$ refers to the smallest integer greater than $a$ (so we overestimate cell size slightly in small samples) and $a_+ = \max\{a, 0\}$, in case the term on the right is negative (under Assumption (ref), this occurs with probability decreasing to zero). Although the i.i.d. case is not considered in our investigation, under no temporal dependence when $T$ is sufficiently large, $\hat{l}$ in the last expression should be approximately equal to one, as expected.
As a direct implication of Lemma (ref), the number of partitions is required to grow quickly as $T$ diverges to infinity. This is a minimal condition compared with block sizes conditions in the literature (see, e.g., C7 in GregoryLahiriNordman18 and Theorem 3.3 in Fitzenberger98). The case of fixed block size is discussed in Fitzenberger98 but no formal results are given for quantile regression. Our approach relies on $b$ growing quickly, which is consistent with the large sample theory of FE-QR which is established under large $T$ panel approximations (Lemma (ref)).
The next result describes the validity of our bootstrap approach.
The result in Theorem (ref) implies the consistency of the PWB for the FE-QR estimator in the case of dependent data. In the next sections, we estimate the partition length $l^{\text{o}}$ and document the performance of the feasible version of the proposed PWB estimator.
Under moment inequalities for mixing processes and the finiteness of $l^\text{o}$, the argument for consistent i.i.d. variance estimation in LamarcheParker23 implies the finiteness of this variance estimator, with the resulting consistent/conservative description depending on whether the data obey Assumption (ref). It is interesting to note that conservative bootstrap inference due to nonexistent bootstrap moments was recently considered in HahnLiao21. Also, MachadoParente5 propose an $L$-estimator of the variance that might apply under more general data conditions than those that we consider, but we leave investigation of such an estimator to future research.
In this section, we report results of several simulation experiments designed to evaluate the finite sample performance of the proposed method. We consider a data generating process similar to the one considered in GalvaoGuVolgushev20. The dependent variable is $y_{it} = \alpha_i + x_{it} + (1 + \zeta x_{it}) u_{it}$, where $x_{it} = 0.5 \alpha_i + z_i + \epsilon_{it}$, and $z_i$ is an independent and identically distributed (i.i.d.) random variable distributed as $\chi^2$ with 3 degrees of freedom ($\chi_3^2$). The corresponding quantile regression function is $Q_y (\tau | x_{it}) = \alpha_{i0} + \beta_0 x_{it}$, where $\alpha_{i0} = \alpha_i + F^{-1}(\tau)$, $\beta_0 = 1 + \zeta F^{-1}(\tau)$, and $F(\cdot)$ denotes the CDF of $u_{it}$. By varying $\zeta \in \{0,0.25\}$, we are able to generate data from two variations of the basic model. The location-shift model assumes $\zeta=0$ and then $\beta_0 = 1$. In the location-scale shift version of the model, $\zeta = 0.25$, and $\beta_0 = 1 + 0.25 F^{-1}(\tau)$ varies by quantile $\tau$.
We generate the observations by combining different distributions for $\alpha_i$ and $x_{it}$. As in LamarcheParker23, we assume that the individual intercept $\alpha_i=i/N$ for $1 \leq i \leq N$, or alternatively, $\alpha_i$ is assumed to be an i.i.d. Gaussian random variable.
We depart from the simulation study in GalvaoGuVolgushev20 by allowing the errors to be serially dependent:
We consider innovation terms $v_{it,u}, v_{it,\epsilon}$ following $\mathcal{N}(0,1)$ and $t_3$ distributions and set the auto-regressive parameters $\rho_{1,u} = \rho_{1,\epsilon}=0.7$ and $\rho_{2,u} = \rho_{2,\epsilon}=0.1$ as in GregoryLahiriNordman18. The variance $\sigma_u^2 = (1 - \rho_{2,u}) \sigma_v^2 / ((1 + \rho_{2,u}) (1 - \rho_{1,u} - \rho_{2,u}) ( 1 + \rho_{1,u} - \rho_{2,u} ))$, which is approximately $2.56$ for normal innovations and three times that for $t_3$ innovations.
Lastly, we use different combinations of $N \in \{5,50\}$ and $T\in \{200,400\}$ and quantile levels $\tau \in \{0.25,0.50,0.75\}$. The results are obtained by using 1000 random samples and 400 bootstrap repetitions.
Figures (ref) and (ref) present results assuming that $\alpha_i = i/N$ and $u_{it} \sim \mathcal{N}(0,\sigma_u^2)$. The figures show coverage probabilities for a nominal level of 90% for the slope parameter over a range of values for the size of the partition $l$. The intervals are constructed using the empirical bootstrap distribution following (ref), and consequently, we present coverages obtained using bootstrap methods. The figures show the performance of the moving block bootstrap (MBB) proposed by Fitzenberger98, extended tapered block bootstrap (ETBB) of Shao10ss generalized to quantile regression as in GregoryLahiriNordman18, weighted block bootstrap (WEB) of GalvaoParkerXiao24, and the proposed PWB estimator. The first three procedures are briefly described in the online appendix. While the weights vary by block size in the MBB, ETBB and PWB approaches, the weight to each unit $i$ is constant in the case of WEB. The figures present results for the location-shift case in the upper panels, and the location-scale shift model in the lower panels. To understand the effect of temporal dependence in panel quantiles, we start with $N=5$ in Figure (ref), and then we increase $N$ to 50 to generate the results shown in Figure (ref).
Figure (ref) shows that the coverage of the MBB and ETBB improve as we increase the partition size from $l = 2$, although it remains roughly constant when $l \geq 10$. We observe the same pattern for PWB, but the method reaches levels closer to the target 90% when $5 \leq l \leq 7$. The performance of WEB is not surprising, since the method relies on cross-sectional variation. We note that when $N$ increases, as in Figure (ref), the performance of WEB significantly improves and it is superior to MBB and ETBB. On the other hand, PWB offers the best coverage in all cases and the excellent performance of the approach is not restricted to a single cell size.
The evidence reveals that by judiciously selecting the partition size, the empirical coverage of the PWB method can reach the target 90% coverage level. Naturally, the best partition size is unknown, so we now extend our investigation to the feasible version of the estimator. Tables (ref) and (ref) present coverage probabilities for a nominal level of 90% obtained by estimating the variance of the estimator as in (ref). Table (ref) presents results for errors distributed as Normal, and Table (ref) presents results for errors distributed as $t_3$. Relative to Figures (ref) and (ref), we are able to expand the evidence to consider different sample sizes, alternative methods for inference on the parameters (i.e., $\textnormal{CI}_{P}$ vs $\textnormal{CI}_{SE}$), and different distributions. The first column shows coverage of the Powell's estimator (PO) of the asymptotic covariance matrix of the FE-QR estimator. PO has been shown to be consistent under the assumption that the data is i.i.d. (Proposition 3.1 in Kato2012). In the next columns, we present the bootstrap approaches. For the MBB and ETBB, we use the {\tt R} package {\tt QregBB} for determining the block length as discussed in GregoryLahiriNordman18. In the case of PWB, coverage is obtained for the selected partition length $\hat{l}$ (summaries are shown in Figure (ref) below).
In all the variations of the model considered in Table (ref), the PWB estimator performs better than the other estimators. Considering for instance the location-shift model, PO does not perform well, although this is expected because the estimator does not account for the temporal dependence of the errors. The weights used by the MBB and ETBB methods are time-varying because they are primarily intended to be used with time-series data. They are not very well suited for panel data because of the inclusion of time-invariant indicators of unit membership --- there is no transformation of the data that can be made beforehand to remove individual effects. As a result, the MBB and ETBB methods undercover. WEB performs well only when $N$ is large relative to $T$, as expected from the evidence presented in Figure (ref).
The results for the location-scale shift model presented in the second panel in Table (ref) are similar. We continue to see that PWB performs better than alternative methods. This conclusion holds when we consider the evidence in Table (ref). The evidence confirms two results. The PWB offers better coverage than WEB in models with $N$ small relative to $T$ and the procedure proposed in this paper is valid for approximating the distribution of the FE-QR estimator.
Finally, using Figure (ref), we present the frequency of partition size selections for PWB corresponding to $N \in \{5,50\}$, $T=200$ and $\alpha_i = i/N$ in Table (ref). The frequencies of $\hat{l}$ corresponding to the other entries in Tables (ref) and (ref) are similar and we do not report them to save space. We first estimate the residuals $\hat{u}_{it}$ using the FE-QR estimator and construct $\check{x}_{it} \psi_\tau(\hat{u}_{it}) \check{x}_{it+k} \psi_\tau(\hat{u}_{it+k})$ for all $k$ considering $\check{x}_{it} = x_{it} - \bar{x}_i$, where $\bar{x}_i = \sum_{t=1}^T x_{it}/T$. This choice does not involve estimation of nuisance parameters and performed well in the simulations. We also employed a triangular kernel $(1-|k|/h)$ with $h=1$. For each $1 \leq i \leq N$, equation (ref) was evaluated at each length $l$ in the contained set $\{1,2,\hdots,L\}$, where $L$ was set to 25 to minimize computational time. We then selected the partition size $\hat{l}_i$ that best fit $\eqref{lhats}$. Finally, we obtained $\hat{l}$ by averaging over individual partition sizes.
The excellent performance of PWB in Tables (ref) and (ref) is conditional on estimates of the partition size obtained from this selection rule, and thus, the evidence confirms that the variance can be approximated empirically by selecting the length of the partition. As expected based on the evidence reported in Figure (ref), Figure (ref) shows that the highest relative frequencies are obtained for values of $\hat{l}$ between 6 and 9. Moreover, we observe that the precision of the selection rule quickly increases when $N$ increases. For instance, the figure shows that over 70% of the selections are in the range $7 \leq \hat{l} \leq 8$ when $N=50$, consistent with the values of $l$ that produce coverage probabilities close to the nominal level of 90% in Figure (ref).
In this section, we illustrate the use of the proposed approach by employing data from a randomized control trial on electricity consumption. We estimate a panel data quantile regression model and apply the method to evaluate how consumers respond to time-of-use electricity pricing. We use a data set that includes $N = 268$ customers observed over $T = 2,160$ time intervals of 30 minutes. Despite the increasing number of studies using household-level panel data with large $T$ kJessoe2014,mHarding2016, dependence within household has been ignored in the empirical literature on quantile regression. Our findings suggest that households do not respond to a modest change in the price of electricity that occurs during the day. However, when the price of electricity increases by 85 percent during the evening, high-usage households reduce their consumption by about 10 percent. If we compare existing approaches with the proposed approach, we find that practitioners would tend to conclude incorrectly that the effects of modest and large changes in the price of electricity are significant across the conditional distribution if they ignore the positive temporal dependence in each household's observations.
We use the CER Smart Metering Project data from the Irish Social Science Data Archive (ISSDA). The experiment was conducted from 2008 to 2011 and we employ data from the period first two months of 2010. Electricity consumption was recorded over 30 minute intervals from Monday to Friday at the residential level. We consider only two treatment types in this paper. The customers selected for the control group had a time-invariant rate of \EUR{0.141} per kilowatt hour (kwh). Customers selected for the treatment group were charged \EUR{0.135} per kwh (Day), with the exception of \EUR{0.110} from 23:00 to 8:00 (Night) and \EUR{0.260} from 17:00 to 19:00 (Peak). Households in the treatment group ($N_1=68$) received an in-home display (IHD) device and an energy usage statement in their bills. Households in the control group ($N_0=200$) did not receive an IHD.
We estimate the following panel data model:
where $y_{igt}$ is the natural logarithm of electricity usage for household $i$ during the interval $t$. The variable $d_i$ indicates treatment status and its effect is not identified because the model includes household fixed effects, $\alpha_i$. We can identify, however, the effect of being treated at different times during the day. The variable $d_{t,1}$ denotes whether $t \in [8-17] \cup [19-23]$ and $d_{t,2}$ denotes if $t \in [17-19]$. The coefficients of interest $(\beta_4(\tau), \beta_5(\tau))$ are identified by the time variation associated with time-of-use pricing across households.
The vector of independent variables $\bm{x}_{it}$ includes average temperature and average relative humidity in Ireland, which are household-invariant regressors. We also include the logarithm of electricity consumption 30 minutes earlier and the logarithm of electricity consumption a day earlier, which accounts for the correlation of electricity usage across weekdays due to daily routines. We have household size, size of the house, and other characteristics of the household and home in the CER data, but we do not include them in the model because these variables are time-invariant and the model includes household fixed effects.
Table (ref) presents fixed effects results for the coefficients $\beta_j$ for $j \in \{1,2,3,4\}$ and standard errors corresponding to the point-estimates in parentheses. The standard errors are obtained using the proposed partitioned wild bootstrap (PWB) procedure. The selected partition size, $\hat{l}$, is constant and equal to 5 across quantiles, with the exception of $\tau=0.9$. The first five columns show quantile regression (FE-QR) results at different quantile levels and the last column presents mean fixed effects regression results. To save space, we do not present results on the control variables included in the regression but the results are available upon request. To examine in more detail the performance of the proposed approach in practice, Figure (ref) compares confidence intervals obtained by alternative approaches.
Looking at the mean fixed effects results in the last column, we see, as expected, that electricity consumption increases during the day and peak hours relative to night hours. These increases are significant and large, ranging from $\exp(0.393)-1 = 22\%$ percent during the day to 48 percent during peak hours. When we compare the mean effect with quantile results shown in the first columns, we observe that these estimated effects vary significantly across quantiles. For instance, if we consider consumption during the day, the effect varies from 7 percent at the 0.1 quantile to 48 percent at the 0.9 quantile. When we focus our attention on the parameters of interest, we find results that are consistent with expectations. Recall that the households in the treatment group were charged a slightly lower price during the day than households in the control group, so it is not entirely surprising to find insignificant results across the conditional distribution. On the other hand, when we consider the effect of the price increase during peak hours, we find significant decreases in electricity consumption in the upper half of the distribution, ranging from a modest 2 percent at the 0.5 quantile to 10 percent at the 0.9 quantile. Consistent with the literature, we find that households reduced their usage and are responsive to temporary price increases.
Figure (ref) presents fixed effects quantile regression results for $(\beta_3(\tau), \beta_4(\tau))$ for 10 equally spaced $\tau$ in the interval 0.10-0.90. The upper panels present results for the effect of consumption over the day and the lower panels present results for the effect of consumption over the peak hours. In each row, we compare 95 percent pointwise confidence intervals obtained by the proposed PWB procedure and the alternative methods for inference available to practitioners. The left panels show results obtained by estimating the asymptotic covariance matrix under two different assumptions. PO denotes the Kernel estimator considered in Proposition 3.1 in Kato2012 which assumes i.i.d. data. The next panels show the bootstrap approaches MBB and ETBB applied to panel data. Finally, the panels in the right show confidence intervals obtained by the proposed PWB estimator and the weighted block bootstrap estimator (WEB) proposed in GalvaoParkerXiao24. The WEB method uses one weight per unit, an i.i.d. non-negative random weight with mean and variance both equal to one.
The results in Figure (ref) highlight the importance of the proposed approach in practice. By ignoring temporal dependence, the kernel estimators produce overly optimistic results, suggesting that the increase in the price during peak hours on electricity usage is negative and significant across the conditional distribution. The application of the bootstrap methods MBB and ETBB lead to a similar observation. As previously discussed, they most likely underrepresent the size of the confidence intervals as well. On the other hand, the WEB and PWB methods offer more conservative inference, suggesting that the effects are only significant in the upper tail. Finally, we find significant efficiency improvements from adopting PWB relative to WEB. We believe that this is due to the fact that the treatment varies predictably over time and the absolute price changes are important.
Statistical inference for quantile regression models of panel data has been a challenge in practice. Existing resampling approaches do not incorporate temporal dependence, which is usually an important concern when practitioners analyze data with a large number of time-series observations. To address this issue, we propose a novel bootstrap method that is easy to implement. The approach re-weights cells of a partition of the estimated residuals to simultaneously satisfy quantile moment conditions conditional on the data and account for how temporal dependence affects standard errors.
We demonstrate that the partitioned wild bootstrap procedure is asymptotically valid for approximating the distribution of the fixed effects estimator. We also investigate the finite sample performance of the approach and find that the new procedure works well. Although we offer an important contribution for panel data observed over many time periods, there are some directions that remain to be investigated. In this work, we focus on the fixed effects estimator, leaving aside the penalized estimator for panel data. We conjecture that the approach can be easily extended to cover the case of penalized estimation of individual effects. We do expect important changes in the consistency result, but we do not offer a formal treatment. We also note that this bootstrap relies mainly on the existence of separable residuals that can be weighted, and therefore this technique might be applied to other scenarios beyond estimation of a quantile regression model. We leave these investigations to future work.