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.
81,045 characters · 13 sections · 32 citation commands
Weighted asymmetric least squares regression with fixed-effects
{\bf Keywords:} Expectile regression, quantile regression, fixed-effects, within-transformation, endogenous model, panel data.
The fixed-effects $(\operatorname{\textit{FE}})$ model is commonly used in econometric to analyze panel data. The $\operatorname{\textit{FE}}$ model has the ability to account for the correlation between the regressors and the omitted (unmeasured) factors which is common in many applications. For example, in econometrics, the education level is known to be correlated with the individual unobserved ability cardEstimatingReturnSchooling2001. In perinatal studies, the birth weight is influenced by maternal genetic warringtonMaternalFetalGenetic2019, which is usually a missing information. Therefore, in such context{\textemdash}where the unmeasured factors are correlated with the regressors, the $\operatorname{\textit{FE}}$ estimator (within-estimator) is unbiased, consistent and computationally efficient CornwellRupert1988.
Several quantile regression $(\operatorname{\textit{QR}})$-based methods koenker_quantile_2004, galvao_penalized_2010, lamarcheRobustPenalizedQuantile2010a have been proposed to overcome the heteroscedasticity problem in the $\operatorname{\textit{FE}}$ framework. However, they fail to extend the favorable properties of the within-estimator and suffer from two significant limitations. First, the fixed-effects $\operatorname{\textit{QR}}$-based methods do not extend the within-transformation strategy to solve the incidental parameter problem. Thus, the fixed-effects $\operatorname{\textit{QR}}$-based methods simultaneously estimate the parameter of interest and the incidental parameter which results in a computationally demanding algorithm. Additionally, the covariance of the $\operatorname{\textit{QR}}$-based methods is based on the random error density function which further adds a computational burden and certain numerical issues Chen2004, YinCai2005a, kocherginskyPracticalConfidenceIntervals2005. Second, the fixed-effects $\operatorname{\textit{QR}}$-based method do not control for the correlation between the individual effects and the regressors. Thus, in the presence of such correlations, the fixed-effects $\operatorname{\textit{QR}}$-based method yields biased and inconsistent estimates.
In this paper, we rely on expectiles to successfully generalize the within-estimator and take into account the heteroscedasticity present in the panel data under the $\operatorname{\textit{FE}}$ framework. To the best of our knowledge, this is the first approach that estimates the marginal effect of the regressors on the response distribution, and generalizes the within transformation strategy in the $\operatorname{\textit{FE}}$ framework.
The expectiles are statistics that characterize the distribution function of a random variable girardFunctionalEstimationExtreme2021. The expectiles and the expectile regression $(\operatorname{\textit{ER}})$ were introduced in the seminal paper by newey_asymmetric_1987. The expectiles and quantiles play similar statistical roles, except that expectiles are weighted averages while quantiles are order statistics. This interpretation difference offers significant computational advantages. In other words, quantiles focus on the ordering of the observations while the expectiles target their values. For instance, the mean is a particular expectile as the median is a particular quantile. The research on expectiles is very active and for further details we refer to GEEE_Barry2020.
Typically quantiles are more robust than expectiles, but as mentioned earlier, the proposed $\operatorname{\textit{QR}}$-based fixed-effects models can not extend the within transformation strategy to solve the computational challenges raised by the incidental parameter problem efficiently. Further, the $\operatorname{\textit{QR}}$-based fixed-effects models fail to control for the correlation between the individual effects and the regressors. Therefore, the expectile-based approach could be an effective alternative for inference in the $\operatorname{\textit{FE}}$ framework.
In this paper, we combine the weighted asymmetric least squares regression $(\operatorname{\textit{ER}})$ and the $\operatorname{\textit{FE}}$ model to propose a new panel model that we call: expectile regression for fixed-effects model $(\operatorname{\textit{ERFE}}).$ The $\operatorname{\textit{ERFE}}$ model retains the attractive properties of the $\operatorname{\textit{FE}}$ model, while accounting for the heteroskedasticity present in the panel data. We derive its asymptotic properties and propose a heterogeneous, consistent, and robust estimator of its variance-covariance matrix. We share our code as a free R package available on GitHub (\url{github.com/AmBarry/erfe}) to simplify its implementation.
Our main contributions are: i. Extension of the within-transformation strategy in the $\operatorname{\textit{ER}}$ framework to solve the incidental parameter problem, offering a significant computational advantage particularly with the advent of high dimensional data, where the sample size can be very large; ii. Elimination of any bias that might result from the correlation between the individual effects and the regressors; iii. Derivation of the asymptotic properties of the $\operatorname{\textit{ERFE}}$ estimators; iv. Proposition of an estimator of its variance-covariance matrix for inference.
Our $\operatorname{\textit{ERFE}}$ model accounts for the omitted time-invariant factors and their correlation with the regressors present in the model. It also captures the heteroskedasticity present in the panel data by estimating the effects of the regressors at the conditional expectiles of the response distribution. Indeed, in the presence of heteroskedasticity, the parameters of the model are function of the asymmetric points, and in this case the $\operatorname{\textit{ERFE}}$ model captures the heteroskedasticity by estimating several regression coefficient vectors at different locations of the conditional response distribution. The $\operatorname{\textit{ERFE}}$ model provides a detailed overview of the regressor effects on the response distribution without making any assumption about the random error distribution. Our $\operatorname{\textit{ERFE}}$ model is computationally efficient and easy to implement, with its available R package. We believe that it will be a useful tool for addressing the heteroskedasticity present in the panel data.
In Section (ref), we introduce the expectile function and the expectile regression model, and then present the expectile regression with fixed-effects (ERFE) model. In Section (ref), we derive the asymptotic properties of the ERFE estimator, and propose an estimator of its variance-covariance (VC) matrix. We present the sample performance of the ERFE estimator in Section (ref) and its application to a real dataset in Section (ref). The conclusions is in Section (ref) and detail of the proofs are in the Supplementary material.
The expectile of level \(\tau\in[0,1]\) of a random variable \(Y\) is defined as the unique solution of
where \(\rho_\tau(t)=\lvert \tau-\mathds{1}(t\leq 0)\rvert \cdot t^2\) is the asymmetric square loss function that assigns weights \(\tau\) and \(1-\tau\) to positive and negative deviations, respectively.
The expectiles summarize the cumulative distribution function of a random variable. In this regard, the expectiles play a similar statistical role to the quantiles, except that quantiles are order statistics while expectiles are weighted averages, and this interpretation difference is accompanied by significant computational advantages. The expectiles generalize the mean which corresponds to the expectile of level $\tau=0.5$ and assigns the same weight to positive and negative deviations. The expectiles are location and scale equivariant, that is for \(s>0 \mbox{ and } t\in\mathbb{R}, \ \mu_{\tau}(sY+t)=s\mu_{\tau}(Y)+t.\) A detailed study of the expectiles can be found in newey_asymmetric_1987.
Once the optimization problem of equation ((ref)) is solved, for a fixed $\tau,$ the expectile of the random variable $Y$ can be defined as a weighted average:
where \(\psi_{\tau}(t)=\lvert \tau-\mathds{1}(t\leq 0)\rvert\) is the check function. The only subtlety is that the weights are random. Given a random sample, \(\lbrace(y_i)\rbrace_{i=1}^{n},\) the corresponding \(\tau\)-th sample expectile
is the weighted mean, where the weights depend on the sample data. For a fixed \(\theta,\) equation ((ref)) is derived as the solution which minimizes the following empirical risk function:
In addition to the expectiles, newey_asymmetric_1987 introduced the expectile linear regression $(\operatorname{\textit{ER}})$ as a tool to study the regressor effects on the response distribution and capture the heteroscedasticity present in the data. Consider the following linear regression model
where \(\boldsymbol{x}_i\) is a \(p\times 1\) vector of regressors, \(y_i\) and \(\varepsilon_i\) are respectively the response variable and the random error with unspecified distribution function. The parameter of interest \(\boldsymbol{\beta} \in \mathbb{R}^p\) is unknown and needs to be estimated. The assumption, \(\mu_{\tau}(\varepsilon_{i})=0,\) ensures that the random error is centered on the \(\tau\)-th expectile. The corresponding $\operatorname{\textit{ER}}$ model, for a fixed \(\tau \in (0,\ 1),\) is given as:
Therefore, the $\operatorname{\textit{ER}}$ estimator \(\widehat{\boldsymbol{\beta}}_{\tau},\) for a fixed \(\tau \in (0,1),\) can be derived by minimizing the following objective function:
over \(\boldsymbol{\beta}_{\tau}\;\in\;\mathbb{R}^p.\) Since the loss function \(\rho_\tau(t)\) is continuously differentiable, we have through the first order condition:
where \(\psi_\tau(t)=\lvert \tau-\mathds{1}(t\leq 0)\rvert\) is the check function. The $\operatorname{\textit{ER}}$ estimator can be computed as an iterated weighted least squares estimators. For the special case of \(\tau=0.5, \ \widehat{\boldsymbol{\beta}}_{0.5}\) is the classical ordinary least squares (OLS) estimator and this makes the $\operatorname{\textit{ER}}$ a natural complement of the OLS regression.
Consider the standard linear fixed-effects model for panel data
where \(y_{ij}\) is the scalar response variable, the vector \(\boldsymbol{x}_{ij}=(x_{ij}^1,x_{ij}^2, \ldots, x_{ij}^p){}^{\text{\sffamily T}} \in \mathbb{R}^p\) is the vector of regressors measured on subject \(i\) at time \(j,\) the parameter \(\alpha_{i}\) is the subject-specific effects parameter, and the variable \(\varepsilon_{ij}\) is a random error with unspecified distribution function. The equation model ((ref)) is conveniently represented in individual notation as:
where \(\boldsymbol{y}_i \mbox{ and } \boldsymbol{\varepsilon}_i\) are \(m \times 1\) vectors, \(\boldsymbol{X}_i\) is an \(m \times p\) design matrix and \(\boldsymbol{Z}_i\) is an \(m \times n\) incidence matrix and \(\boldsymbol{\alpha}\) is an \(n \times 1\) subject-specific effects vector. We can also stack all the data and represent the equation model ((ref)) as:
where \(\boldsymbol{y}\) and \(\boldsymbol{\varepsilon}\) are \(N\times 1\) vectors, \(\boldsymbol{X}\) and \(\boldsymbol{Z}\) are respectively \(N\times p\) and \(N\times n\) matrices with \(N=mn.\) The incidence matrix \(\boldsymbol{Z}\) identifies the \(n\) distinct subjects of the sample.
The fixed-effects $\boldsymbol{Z}\boldsymbol{\alpha}$ of model equation ((ref)) is infinite in nature and is potentially correlated with the regressors of the model. The traditional estimation method used to overcome this issue is the within-transformation strategy. This technique consists of pre-multiplying both sides of equation model ((ref)) by the idempotent matrix \(\boldsymbol{M}_{\boldsymbol{Z}}=\mathbb{I}_N-\boldsymbol{Z}(\boldsymbol{Z}{}^{\text{\sffamily T}}\boldsymbol{Z})^{-1}\boldsymbol{Z}{}^{\text{\sffamily T}}\) to eliminate the infinite-dimensional parameter, and then applies the ordinary least squares (OLS) regression to the transformed data. The model that results from this transformation is represented as:
where \(\boldsymbol{y}^{*}=\boldsymbol{M}_{\boldsymbol{Z}}\boldsymbol{y}\) and \(\boldsymbol{X}^{*} \mbox{ and } \boldsymbol{\varepsilon}^{*}\) are defined similarly. The OLS estimator of the fixed-effects model, known as the within-transformation estimator, is given as:
The within-transformation estimator is consistent and asymptotically normally distributed Greene2011. The within-transformation estimator is computationally efficient and is not affected by any bias resulting from the correlation between the individual effects and the regressors. The within-transformation technique does not allow estimation of the time-invariant regressors which could be seen as a limitation. However, this can be a strength when the number of time-invariant confounders is large, and when the collection of some of these variables (genetic factor) is complex and costly bruderlFixedEffectsPanelRegression2014. In the following subsection, we present the expectile regression for fixed-effects $(\operatorname{\textit{ERFE}})$ model and derive the iterative-within-transformation $\operatorname{\textit{ERFE}}$ estimator.
The $\operatorname{\textit{ERFE}}$ model of the linear fixed-effects model is defined, for fixed \(\tau \in (0,1),\) as:
In this setting the parameter \(\boldsymbol{\beta}_{\tau} \in \mathbb{R}^p\) captures the influence of the regressors \(\boldsymbol{x}_{ij}\) on the location, scale, and shape of the conditional distribution of the response variable \(y_{ij}.\) The subject-specific effects \(\boldsymbol{\alpha}\) is assumed to be independent of \(\tau\) across the percentiles and to have a pure location-shift effect on the conditional expectile of the response. Assuming a \(\tau\)-dependency of the subject-specific effects implies estimating its distribution with \(m\) number of within-subject observations, which is relatively small in most applications. Take note that no assumption is made about the shape of the response distribution.
The corresponding $\operatorname{\textit{ERFE}}$ estimator of model equation ((ref)) is defined as the vector minimizing the following objective function:
Since the loss function \(\rho_{\tau}(\cdot)\) is continuously differentiable, we can apply the first-order condition and derive the resulting $\operatorname{\textit{ERFE}}$ estimator \(\widehat{\boldsymbol{\beta}}_{\tau}\), which is defined as:
where the diagonal check function matrix is:
The projection matrix \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\tau)\) and its complement \(\widehat{\boldsymbol{P}}_{\boldsymbol{Z}}(\tau)\) are idempotent matrices and are defined as:
The function \(\boldsymbol{\Psi}_{\tau}(\cdot)\) defined in equation ((ref)) depends on the subject-specific parameter estimator \(\widehat{\boldsymbol{\alpha}}\ \) which, by the first-order condition of equation ((ref)), verifies the relationship:
Now, using equation ((ref)), the argument of the check function matrix can be written as
Therefore, the incidental parameter estimate is eliminated from the expression of equation ((ref)) of the $\operatorname{\textit{ERFE}}$ estimator. Now, using the following relationship:
and the idempotent property of the projection matrix \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\tau),\) we can rewrite the $\operatorname{\textit{ERFE}}$ estimator as:
In summary, the within-estimator is extended to the $\operatorname{\textit{ER}}$ framework. The strategy is derived by applying the projection matrix \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\tau)\) to the initial data \([\boldsymbol{y},\boldsymbol{X}],\) to eliminate the subject-specific effects parameter. Additionally, like the $\operatorname{\textit{ER}}$ estimator in equation ((ref)), the within $\operatorname{\textit{ERFE}}$ estimator can be computed iteratively using the iterative weighted least squares algorithm. The detailed algorithm for computing the iterative-within-transformation $\operatorname{\textit{ERFE}}$ estimator is summarized in the following stepwise procedures.
The parameter $\zeta$ is the convergence tolerance and the default value in our code implementation is set to $10^{-7}. \ $ In practice, Algorithm (ref) is computationally efficient and usually the number of iterations required to achieve convergence is between 3 and 5. Note that when \(\tau=0.5\) we have \(\boldsymbol{\Psi}_{\tau}=0.5\mathbb{I}_{N}\) and the iterative within-transformation $\operatorname{\textit{ERFE}}$ estimator is nothing else than the OLS within-transformation estimator.
From the above development, the multiplication of a vector (say \(\boldsymbol{y}\)) by the matrix \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\tau)\) deviates that vector from its expectile as shown by the following expression:
We can see, from this expression, how the projection matrix \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\tau)\) eliminates the subject-specific effects parameter and any other time-invariant regressors from the initial model. For a matrix, like the design matrix \(\boldsymbol{X},\) the transformation is applied column-wise.
ERFE model for a sequence of expectiles
The preceding development shows that the classical within-transformation strategy can be generalized in the $\operatorname{\textit{ERFE}}$ framework. Now, we present the $\operatorname{\textit{ERFE}}$ estimator for a sequence of expectiles using the transformed data. The sequence of expectiles, for example the mean and a few expectiles above and below it, is necessary in the description of the conditional distribution of the response variable and for capturing the data heteroscedasticity. In addition, the simultaneous estimation allows the multiple expectiles to share strength among each other and to gain better estimation accuracy than individually estimated expectile functions LiuWu2011.
The iterative within-transformation $ERFE$ estimator \(\widehat{\boldsymbol{\beta}}_{\boldsymbol{\tau}}=[\widehat{\boldsymbol{\beta}}_{\tau_1}{}^{\text{\sffamily T}}, \ldots, \widehat{\boldsymbol{\beta}}_{\tau_q}{}^{\text{\sffamily T}}]{}^{\text{\sffamily T}}\) for a sequence of asymmetric points \(\boldsymbol{\tau}=(\tau_1,\ldots,\tau_q){}^{\text{\sffamily T}}\) is defined as the minimum of the following objective function:
The vector \(\boldsymbol{v}=(v_1, \ldots, v_q){}^{\text{\sffamily T}}\) is the vector of weights controlling the relative influence of the q asymmetric points \(\lbrace \tau_1, \ldots, \tau_q\rbrace\) and it choice depends on the research question. For a sequence of expectiles, the iterative within-transformation $\operatorname{\textit{ERFE}}$ estimator is defined as:
where \(\boldsymbol{\Psi}_{\boldsymbol{\tau}}(\widehat{\boldsymbol{\varepsilon}_{\boldsymbol{\tau}}^{*}})=\Big[\operatorname{diag}\big(\boldsymbol{\Psi}_{\tau_k}(\widehat{\boldsymbol{y}^{*}}-\widehat{\boldsymbol{X}^{*}} \widehat{\boldsymbol{\beta}}_{\tau_{k}})\big)\Big]_{k=1}^{q}, \; \boldsymbol{V}=[\operatorname{diag}(v_k)]_{k=1}^{q}\) and the transformed data \([(\mathds{1}_q\otimes\widehat{\boldsymbol{y}^{*}}),(\mathbb{I}_q\otimes\widehat{\boldsymbol{X}^{*}})]\) is obtained by pre-multiplying the matrix \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\boldsymbol{\tau})\) to the initial data \([\boldsymbol{y},\boldsymbol{X}].\) The projection matrix is defined as \(\widehat{\boldsymbol{M}}_{\boldsymbol{Z}}(\boldsymbol{\tau})=\mathbb{I}_{nmq}-\widehat{\boldsymbol{P}}_{\boldsymbol{Z}} (\boldsymbol{\tau})\) and
In this section, the asymptotic properties of the $\operatorname{\textit{ERFE}}$ estimator are presented. As stated by koenker_quantile_2004, the presence of the incidence parameter, which has an infinite dimension, can raise some challenges. For this reason, we present first the asymptotic results of the $\operatorname{\textit{ERFE}}$ estimator in the simplest case, namely for a single \(\tau.\) We then present the asymptotic properties of the $\operatorname{\textit{ERFE}}$ estimator for a simultaneous sequence of asymmetric points \(\boldsymbol{\tau}=(\tau_1,\ldots,\tau_q).\) The section ends with the suggestion of an estimator of the covariance matrix for the $\operatorname{\textit{ERFE}}$ estimator. All the proofs are available in the Supplementary file.
Asymptotics for a single expectile
In the following, the asymmetric square-loss function of the $\operatorname{\textit{ERFE}}$ model,
is replaced in term of optimization by the equivalent loss function
where \(\mu_{ij\tau}=\boldsymbol{x}_{ij}{}^{\text{\sffamily T}} \boldsymbol{\beta}_{\tau} + \boldsymbol{z}_{ij}{}^{\text{\sffamily T}}\boldsymbol{\alpha}.\) Now, observe that the following estimator
minimizes the new objective function
The asymptotic theory of the $\operatorname{\textit{ERFE}}$ estimator is derived using this new objective function ((ref)) and under the following assumptions.
A1. The data \(\lbrace (\boldsymbol{y}_i,\boldsymbol{X}_i)\rbrace_{i=1}^{n}\) are independent across \(i,\) and,
where $\boldsymbol{\varepsilon}_{i\tau}=(\varepsilon_{i1\tau},\ldots,\varepsilon_{im\tau}){}^{\text{\sffamily T}}, \ \varepsilon_{ij\tau}=y_{ij}-\boldsymbol{x}_{ij}{}^{\text{\sffamily T}}\boldsymbol{\beta}_{\tau}\ $ and $ \ \boldsymbol{\Psi}_{\tau }(\boldsymbol{\varepsilon}_{i\tau})=\Big[\operatorname{diag}(\psi_{\tau}(\varepsilon_{ij\tau}))\Big]_{j=1}^{m}.$
A2. The limiting forms of the following matrices are positive definite
where \(\boldsymbol{\Sigma}_{\tau}=\operatorname{\mathbb{V}\text{ar}}[\boldsymbol{\Psi}_{\tau}(\boldsymbol{\varepsilon}_{\tau})\boldsymbol{\varepsilon}_{\tau}]=\operatorname{diag}[\boldsymbol{\Sigma}_{i\tau }]_{i=1}^{n}.\)
A3. The norm of the regressors is bounded by a positive constant $M, \ $ \(\max_{i,j} \left\lVert\boldsymbol{x}_{ij}\right\rVert< M.\)
The stated assumptions A1-A3 are standard for panel data models koenker_quantile_2004. Condition A1 ensures independence across individuals, but allows a within-subject dependency and heterogeneity across individuals. Condition A2 is a full rank condition and is used to invoke the Lindeberg-Feller Central Limit Theorem. We observe that, when \(\ \tau =1/2\ \) then \(\boldsymbol{D}_1 (\tau)\) simplifies and Condition A2 reduces to a condition on the matrices \( \ \boldsymbol{X}{}^{\text{\sffamily T}}\boldsymbol{X}/nm \ \) and \( \ \boldsymbol{Z}{}^{\text{\sffamily T}}\boldsymbol{Z}/m.\ \) Condition A3 is useful both for the application of the Lindeberg-Feller Central Limit Theorem and for ensuring the finite dimensional convergence of the objective function.
To show the closed form of the above matrix \(\Big[\boldsymbol{D}_1^{-1}(\tau)\boldsymbol{D}_0(\tau)\boldsymbol{D}_1^{-1}(\tau)\Big]_{22}\) assume that the limiting forms of the following matrices are positive definite
where \(\boldsymbol{M}_{\boldsymbol{Z}}(\tau)=\mathbb{I}-\boldsymbol{P}_{\boldsymbol{Z}}(\tau) \mbox{ and } \boldsymbol{P}_{\boldsymbol{Z}}(\tau)=\boldsymbol{Z}\Big[\boldsymbol{Z}{}^{\text{\sffamily T}}\operatorname{\mathbb{E}}[\boldsymbol{\Psi}_{\tau}(\boldsymbol{\varepsilon})] \boldsymbol{Z}\Big]^{-1} \boldsymbol{Z}{}^{\text{\sffamily T}}\operatorname{\mathbb{E}}[\boldsymbol{\Psi}_{\tau}(\boldsymbol{\varepsilon})].\)
Under the above conditions and the conditions of Theorem (ref) it follows that:
Asymptotics for several expectiles
The asymptotic properties of the $\operatorname{\textit{ERFE}}$ estimator for a sequence of asymmetric points \(\boldsymbol{\tau}=(\tau_1,\cdots,\tau_q)\) are derived using the transformed data, \(\ [\boldsymbol{y}^{*};\boldsymbol{X}^{*}],\ \) where \(\boldsymbol{X}^{*}=\boldsymbol{M}_{\boldsymbol{Z}}(\tau)\boldsymbol{X} \ \mbox{ and } \ \boldsymbol{y}^{*} = \boldsymbol{M}_{\boldsymbol{Z}}(\tau)\boldsymbol{y}.\) Both projection matrices \(\boldsymbol{M}_{\boldsymbol{Z}}(\tau) \mbox{ and } \boldsymbol{P}_{\boldsymbol{Z}}(\tau)\) are idempotent and are defined as:
where \(\boldsymbol{\varepsilon}^{*}_{\tau}=\boldsymbol{y}^{*}- \boldsymbol{X}^{*}\boldsymbol{\beta}_{\tau}.\)
A robust estimator of the covariance matrix is also proposed. Assume the following conditions.
B1. The data \(\lbrace (\boldsymbol{y}_i,\boldsymbol{X}_i) \rbrace_{i=1}^{n}\) are independent across \(i\) and,
where $\boldsymbol{\varepsilon}_{i\boldsymbol{\tau}}^{*}=\Big(\boldsymbol{\varepsilon}_{i\tau_1}^{*}{}^{\text{\sffamily T}},\ldots,\boldsymbol{\varepsilon}_{i\tau_q}^{*}{}^{\text{\sffamily T}}\Big){}^{\text{\sffamily T}}, \; \boldsymbol{\varepsilon}_{i\tau_k}^{*}=(\varepsilon_{i1\tau_k}^{*},\ldots,\varepsilon_{im\tau_k}^{*}){}^{\text{\sffamily T}}, \; \varepsilon_{ij\tau_k}^{*}=y_{ij}^{*}-\boldsymbol{x}_{ij}^{*}{}^{\text{\sffamily T}}\boldsymbol{\beta}_{\tau_k} \ $ and $ \ \boldsymbol{\Psi}_{\boldsymbol{\tau}}(\boldsymbol{\varepsilon}_{i\boldsymbol{\tau}}^{*})= [\operatorname{diag}(\boldsymbol{\Psi}_{\tau_k}(\boldsymbol{\varepsilon}_{i{\tau_k}}^{*}))]_{k=1}^{q}.$
B2. The limiting forms of the following matrices are positive definite
B3. The norm of the regressors is bounded by a positive constant $M, \ $ \(\max_{\substack{ 1\leq i \leq n \\ 1 \leq j \leq m} }\left\lVertx_{ij}^{*}\right\rVert < M.\)
In order to use the $\operatorname{\textit{ERFE}}$ estimator to make inference, an estimator of its covariance matrix is presented in Theorem (ref). This will make it possible to construct large sample confidence intervals or hypothesis tests. The proposed covariance matrix estimator is robust and consistent, and is a generalization of the commonly advocated covariance matrix estimator proposed by White1980.
We end this section with the result for a single \(\tau.\)
In this section we conducted a simulation study to evaluate the performance of the $\operatorname{\textit{ERFE}}$ estimator. We started by presenting the simulation design, then the metrics to evaluate the estimators and the results.
The random samples were generated from the following linear model:
We considered two versions of model equation ((ref)) according to the heteroscedastic parameter \(\gamma\in\lbrace 0, \ 3/10 \rbrace.\) The value of \(\gamma=0\) corresponds to a location shift model \((M_0)\) where the regressors are uncorrelated to the random error. The model \((M_0)\) is used to assess the performance of the estimators for a homoscedastic scenario. In contrast, when the value of \(\gamma = 3/10,\) then there is a correlation between the predictor $x_{2}$ and the random error. In that case, the model is a location-scale shift model \((M_{3/10})\) and is set to assess the performance of the estimators in the presence of heteroscedasticity.
In the location shift scenario, the $\operatorname{\textit{ERFE}}$ model corresponds to \(\mu_{\tau}(y_{ij})=x_{ij1}\beta_1 + x_{ij2}\beta_2 + \alpha_i + \mu_{\tau}(\varepsilon_{ij})\) where only the intercept term, \(\beta_{0\tau}=\alpha_i +\mu_{\tau}(\varepsilon_{ij}),\) varies with \(\tau\) and the expectile functions are parallel lines. In the location-scale shift scenario, the related $\operatorname{\textit{ERFE}}$ model is defined as: \(\mu_{\tau}(y_{ij})= x_{ij1}\beta_1 + x_{ij2}\beta_{2\tau} + \alpha_i + \mu_{\tau}(\varepsilon_{ij})\) where the intercept $\beta_{0\tau}= \alpha_i + \mu_{\tau}(\varepsilon_{ij})$ and $\beta_{2\tau}= \beta_{2} + \gamma\mu_{\tau}(\varepsilon_{ij}).$ Therefore, in the presence of heteroscedasticity both the intercept and the slope of the predictor $x_2$ vary with \(\tau.\)
The parameters are set to $\beta_1=0.6$ and $\beta_2=1,$ and the corresponding regressors are generated from a non-central student distribution with 3 degree of freedom $(\mathcal{T}_2(1.3))$ and a normal distribution \((\mathcal{N}(2,\ 1.5)),\) respectively. The individual-specific effects parameter \(\alpha\) is generated from a normal distribution \((\mathcal{N}(1,\ 1)),\) and is correlated $(\rho=0.5)$ to the predictor $x_2.$ Indeed, in real data applications it is more likely that omitted factors are correlated with regressors in the model. The random error \(\varepsilon\) of the model equation ((ref)) is generated from three different distributions: normal distribution \((\mathcal{N}(0,1)),\) Student distribution $(\mathcal{T}_3)$ with 3 degrees of freedom, and chi-squared distribution \((\chi^2_3)\) with \(3\) degrees of freedom. We have set the sample size and the repeated measurements to \(n\times m\in \ \lbrace 100,\ 250,\ 500 \rbrace \times \lbrace 5,\ 15,\ 30 \rbrace.\) The extensive simulation was carried out with 400 replications. In each case the focus is on the regressor effects at the asymmetric points $\tau\in\lbrace 0.1, \ 0.3, 0.5, \ 0.8, \ 0.9 \rbrace.$
All simulations were conducted using high performance computing clusters provided by Calcul Quebec and Compute Canada. All computations were performed with the R (v3.6.0) statistical programming language rcran. The implemented R package erfe that comes with this manuscript is publicly available on GitHub at \url{https://github.com/AmBarry/erfe}.
We compared our $\operatorname{\textit{ERFE}}$ model to the quantile regression with fixed-effects $(\operatorname{\textit{QRFE}})$ model proposed by koenker_quantile_2004. The $\operatorname{\textit{QRFE}}$ model estimated the parameter of interest and the nuisance parameter of the model which could be computationally demanding as the sample size increased. We also considered the expectile regression model $(\operatorname{\textit{ER}})$ and the quantile regression model $(\operatorname{\textit{QR}}),$ which ignored the individual fixed-effects parameter. Given that expectile and quantile of the same level $\tau$ were generally different, we carried out the appropriate conversions between the asymmetric points and the percentiles to ensure that the expectile-based regressions and the quantile-based regressions estimated the same statistics (that is quantiles and expectiles are identical). For example, the Gaussian quantiles of level \(\tau=(0.33, \ 0.5, \ 0.67)\) are identical to the Gaussian expectiles of level \(\tau=(0.25,\ 0.5,\ 0.75).\) In other words, the $\operatorname{\textit{ER}}$ based-model and the $\operatorname{\textit{QR}}$ based-model estimate the same locations of the response distribution.
We evaluated the quality of the estimators by reporting the distribution of their coefficient estimate as box-plots. We also evaluated the performance of the asymptotic standard error $(\operatorname{\textit{SE}})$ presented in Theorem (ref) by reporting the distribution of the ratio between the asymptotic standard error $(\operatorname{\textit{SE}})$ and the Monte Carlo standard deviation $(\operatorname{\textit{SD}})$ defined as:
where $\overline{\widehat{\beta}}_{k\tau} = \frac{1}{400} \sum_{j=1}^{400} \widehat{\beta}^{(j)}_{k\tau}.$
We estimated the $\operatorname{\textit{ERFE}}$ model with the erfe package and the $\operatorname{\textit{QRFE}}$ model with the rqpd package rqpd. The $\operatorname{\textit{ER}}$ model and the $\operatorname{\textit{QR}}$ model was obtained from the well-known packages: expectreg expectreg and quantreg quantreg, respectively.
We present here the results related to the Gaussian random error and we brought the results for the Student and Chi-square random errors in the Supplementary file. Figure (ref) and Figure (ref) report the distribution of the coefficient estimates in the location-shift and location-scale-shift scenarios, respectively.
In the location-shift scenario, we observe that the coefficient estimates of our $\operatorname{\textit{ERFE}}$ model are centered around the true value of the parameters with a small interquartile range. We also observe that the coefficient estimate of the $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ models are centered around the true value for the parameter $\beta_1$ only. We notice that the coefficient estimates of the $\operatorname{\textit{QRFE}}$ model are not close to the true value of the parameters except when $\tau=0.5$ for the parameter $\beta_1$ only. In other words our $\operatorname{\textit{ERFE}}$ model performs well in estimating the parameter coefficients of the model in the location-shift scenario. The $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ models perform similarly in the location-shift scenario, with an unbiased estimator for the parameter $\beta_1$ and a biased estimator for the parameter $\beta_2.$ The $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ models do not account for the individual fixed-effects which are correlated to the regressor $x_2,$ which could explain the bias for the parameter $\beta_2.$ In contrast, the $\operatorname{\textit{QRFE}}$ model performed poorly in estimating the parameter coefficients of the model in the location-shift scenario. The $\operatorname{\textit{QRFE}}$ model includes the individual fixed-effects in its specification but, similarly to the random-effect model, it did not account for the dependence between the individual fixed-effects and the regressors of the model. This could explained the poor performance of the $\operatorname{\textit{QRFE}}.$
Indeed, similar to the within-estimator, the $\operatorname{\textit{ERFE}}$ model transforms the data by subtracting the person-specific expectile of level $\tau$ from the observed values of each variable and then applied the $\operatorname{\textit{ER}}$ method to the de-expectilized model given by:
where $y_{ij}^{*} = y_{ij} - \widehat{\mu}_{\tau}(y_{ij}), \ x_{ij}^{*}$ and $\varepsilon_{ij}^{*}$ are defined similarly. This transformation concentrated out the individual fixed-effects and any bias that could result from its association with the regressors.
The $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ models do not take into account the individual fixed-effects parameter, which is included in the random error component. Since, the individual fixed-effects parameter is correlated to the predictor $x_2,$ then the random error of the model equation ((ref)) is also correlated to the predictor $x_2$ of the model. Hence, the coefficient estimate of the $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ methods for the parameter $\beta_2$ is biased.
Consider, the reformulation of model equation ((ref)) in the location-shift scenario:
where $Q_{\tau}(\alpha_i| x_{ij1}, x_{ij2})$ is the quantile of the individual fixed-effects $\alpha_i$ of level $\tau$ and $\eta_{ij}$ the new random variable. The corresponding $\operatorname{\textit{QRFE}}$ model, for a fixed $\tau,$ can be specified as: $ Q_{\tau}(y_{ij}| x_{ij1}, x_{ij2})=x_{ij1}\beta_1 + x_{ij2}\beta_2 + f(x_{ij1}, x_{ij2}), \ $ where $ \ f(x_{ij1}, x_{ij2})=Q_{\tau}(\alpha_i| x_{ij1}, x_{ij2})$ (since the individual fixed-effects are correlated to the regressors). Thus, in this context, the coefficient estimates of the $\operatorname{\textit{QRFE}}$ model would be biased.
Figure (ref) report the distribution of the coefficient estimates in the location-scale-shift scenario. Again, we observe that the $\operatorname{\textit{ERFE}}$ model performs well in estimating the parameter coefficients of the model and outperformed its competitors. The apparent bias of the $\operatorname{\textit{ERFE}}$ estimator for the parameter $\beta_2$ is due to the effect of the heteroscedasticity in this formulation of the model and is not surprising. Indeed, in the location-scale-shift scenario, because of the correlation between the predictor $x_2$ and the error term, the parameter of the predictor $x_2$ is function of the asymmetric point and then different to $\beta_2$ except when $\tau=0.5$ for the symmetric distributions (Normal and Student), where the expectile of level $\tau=0.5$ is zero. The same remark could be applied to the other methods which in addition did not account for the individual fixed-effects $(\operatorname{\textit{ER}} \ \mbox{ and } \ \operatorname{\textit{QR}})$ or its correlation with the regressors $(\operatorname{\textit{QRFE}}).$ We observed similar results for the Student and Chi-square random errors (results are available in the Supplementary file).
Overall, the $\operatorname{\textit{ERFE}}$ model outperforms its competitor and extends the favorable properties of the fixed-effects model. The $\operatorname{\textit{ERFE}}$ model accounts for the time-invariant omitted variables and for the heteroscedasticity present in the data.
To evaluate the asymptotic standard error $(\operatorname{\textit{SE}})$ of the $\operatorname{\textit{ERFE}}$ parameter estimates, we use the Monte Carlo standard deviation $(\operatorname{\textit{SD}})$ as a benchmark and present the distribution of the ratio $\frac{\operatorname{\textit{SE}}}{\operatorname{\textit{SD}}}$ as an error plot centered at the mean, Figure (ref) and Figure (ref). In general the error plots of the $\operatorname{\textit{ERFE}}$ model and the $\operatorname{\textit{QRFE}}$ model are centered around 1, which means that on average the asymptotic standard error $\operatorname{\textit{SE}}$ and the Monte Carlo standard deviation $\operatorname{\textit{SD}}$ are identical. However, we observe that the error plots of the $\operatorname{\textit{ER}}$ model and the $\operatorname{\textit{QR}}$ model are not centered around the mean for the $\beta_2$ parameter and the range of their error plot is generally larger. Similar performances were observed for the Student and Chi-Squared random error which results can be found in the Supplementary file.
We end this section by comparing the run-times of the $\operatorname{\textit{ER}}$-based algorithms and the $\operatorname{\textit{QR}}$-based algorithms. We fitted the methods to a dataset $(n=300, \ m=10)$ generated by a location-shift model with a Gaussian random error. We used the microbenchmark package microbenchmark2019 with 100 replications to evaluate the computation time of the different algorithm. The results in Figure (ref) show that the cross-sectional algorithms $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ are the fastest algorithms, and our $\operatorname{\textit{ERFE}}$ algorithm is faster than the $\operatorname{\textit{QRFE}}$ algorithm. We also performed the comparison for a larger sample size $(n>2500),$ but the algorithm stopped due to a shortage of memory for the $\operatorname{\textit{QRFE}}$ algorithm. This problem has also been reported by canaySimpleApproachQuantile2011a.
Returns to schooling also known as returns to education is a topic widely studied in empirical economics. It is often presented in standard econometric textbooks baltagi2008, Greene2011, CameronTrivedi2005 as an example of an endogeneity model. Indeed, there is a potential correlation between individual's ability and the other regressors such as education. In the presence of endogeneity, the $\operatorname{\textit{FE}}$ model is often preferred than other mean regression models for panel data. Despite the fact that it does not estimate the effect of the time-invariant regressors, the $\operatorname{\textit{FE}}$ estimator is consistent even if the individual effects are correlated with the regressors of the model baltagi2008.
In this section, we replicated BaltagiKhantiAkom1990's study using the Panel Study of Income Dynamics (PSID) dataset. The dataset is a cohort of 595 individuals observed over the period 1976--1982. The respondents, aged between 18 and 65 in 1976, are those who reported a positive wage in private non-farm employment for all 7 years, CornwellRupert1988.
The log wage is the dependent variable and is regressed on weeks worked (WKS), years of full-time work experience (EXP), occupation (OCC=1, if the individual is in a blue-collar occupation), residence (SOUTH = 1, SMSA = 1, if the individual resides in the South, or in a standard metropolitan statistical area), marital status (MS = 1, if the individual is married), industry (IND = 1, if the individual works in a manufacturing industry), and union coverage (UNION = 1, if the individual's wage is set by a union contract).
We fitted the $\operatorname{\textit{ERFE}}$ model to the PSID dataset. In addition to the regressor effects on the average salary baltagi2008, BaltagiKhantiAkom1990, CornwellRupert1988, the $\operatorname{\textit{ERFE}}$ model captures the regressor effects on the entire wage distribution. Consequently, the $\operatorname{\textit{ERFE}}$ model controls for the endogeneity resulting from unmeasured factors and captures the heterogeneity present in the data. The corresponding Mincer equation of the $\operatorname{\textit{ERFE}}$ model, for a fixed $\tau\in (0,1),$ is specified as:
where the initial model is transformed to eliminate the individual effects.
We estimated the conditional expectiles of the log wage distribution using 91 asymmetric points \( (\tau \in (0.05,\ 0.06,\ 0.07,\ \ldots,\ 0.95)).\) We generated the confidence intervals using the asymptotic standard error of the $\operatorname{\textit{ERFE}}$ model. For comparison, we also fitted the $\operatorname{\textit{ER}}$ model, the $\operatorname{\textit{QR}}$ model and the $\operatorname{\textit{QRFE}}$ model. Notice that the covariance matrix of the $\operatorname{\textit{QR}}$-based method depends on the random error density function which add a computational burden and some numerical issues Chen2004, YinCai2005a, kocherginskyPracticalConfidenceIntervals2005. We used a kernel estimate of the sandwich as proposed by barnettNonparametricSemiparametricMethods1991 to compute the standard error of the $\operatorname{\textit{QR}}$ estimates and the generalized bootstrap of boseGeneralizedBootstrapEstimators2003 to compute the standard error of the $\operatorname{\textit{QRFE}}$ estimates. Moreover, since an expectile of level \(\tau\) is not necessarily equal to a quantile of the same level, the comparison between the $\operatorname{\textit{ER}}$-based results and the $\operatorname{\textit{QR}}$-based results must be done globally.
Figure (ref) and Figure (ref) display the coefficient estimates of the regressors obtained by fitting the $\operatorname{\textit{ER}}$-based methods while Figure (ref) and Figure (ref) display the coefficient estimates of the regressors obtained by fitting the $\operatorname{\textit{QR}}$-based methods. The overall results show the potential of both $\operatorname{\textit{ER}}$-based and $\operatorname{\textit{QR}}$-based methods to reveal the heterogeneous regressor effects on the response distribution and therefore to capture the heteroscedasticity present in the data. We observe that the parameter estimates of some regressors (UNION, IND and SOUTH, for example) vary with respect to the asymmetric points or percentiles suggesting the presence of heteroscedasticity in the data. For example, we observe that the parameter estimates of the UNION variable decrease with respect to the asymmetric points or percentiles suggesting that individuals with low salary have more advantage of being unionized than individuals with high salary. We also observe that the parameter estimates of some regressors may vary a little or not at all with respect to the asymmetric points, suggesting that the mean effect of these regressors would be enough to summarize their relationship with the response variable.
We also observe that the curves of the $\operatorname{\textit{ER}}$-based results are smoother than those from the $\operatorname{\textit{QR}}$-based method which seem to be more wiggly and unstable. Indeed, the $\operatorname{\textit{QR}}$-based results is more volatile and it is more difficult to identify an overall trend of the heterogeneity of the regressor effects. For example, the $\operatorname{\textit{QRFE}}$ parameter estimates of the IND variable is decreasing between the percentiles 0.1 and 0.25, and then increasing between the percentiles 0.25 and 0.9.
Despite the similar trend, the parameter estimates of the different methods have different statistical properties. The coefficient estimates of the $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ have similar range and are biased upward. For example, the $\operatorname{\textit{ER}}$ and $\operatorname{\textit{QR}}$ coefficient estimates of the WKS variable fluctuate between 0.0025 and 0.005. While the $\operatorname{\textit{QRFE}}$ coefficient estimates of the WKS variable vary between 0.06 and 0.07, 10 times higher than that of the $\operatorname{\textit{ERFE}}$ parameter estimates. Therefore, the $\operatorname{\textit{QRFE}}$ coefficient estimates is severely biased because of its inability to account for the correlation between the individual fixed-effects and some regressors in the model.
This results are in line with the simulation results, where we observed that the $\operatorname{\textit{ER}}$ and the $\operatorname{\textit{QR}}$ estimates have similar and lower bias than the $\operatorname{\textit{QRFE}}$ estimates which have higher bias particularly when the individual fixed-effects is correlated to the regressors in the model.
In summary, the data analysis shows that some parameter estimates vary according to the asymmetric points or the percentiles. Therefore, we need to consider beyond the mean or median regression in order to capture the heterogeneity present in the data. The $\operatorname{\textit{FE}}$ model, like other methods that estimate the mean effect, is not sufficient to analyze the returns to schooling because the impact of most of the regressors vary across the wage distribution.
We introduced the $\operatorname{\textit{ERFE}}$ model which inherits the attractive properties of the weighted asymmetric least squares regression $(\operatorname{\textit{ER}})$ and the $\operatorname{\textit{FE}}$ model. As with the $\operatorname{\textit{FE}}$ model, the ERFE model is an endogenous model that takes into account the possible correlation between the omitted time-invariant variables and the regressors included in the model. In addition, the $\operatorname{\textit{ERFE}}$ model estimates the regressor effects on the conditional expectiles of the response distribution allowing to study the influence of the regressors on the location, scale, and shape of the conditional response distribution.
We derived the asymptotic properties of the $\operatorname{\textit{ERFE}}$ estimator and suggest an estimator of its variance covariance matrix. We showed that the $\operatorname{\textit{ERFE}}$ estimator is an iterative-within-transformation estimator. That is, the $\operatorname{\textit{ERFE}}$ estimator can be derived by using iteratively the within-transformation strategy to concentrate out the incidental parameter from the model. The $\operatorname{\textit{ERFE}}$ model is computationally efficient and easy to implement. See our GitHub for a free R package that simplifies the implementation (\url{github.com/AmBarry/erfe}).
The exhaustive simulations showed that the $\operatorname{\textit{ERFE}}$ estimator outperformed its competitors, including the $\operatorname{\textit{QRFE}}$ estimator in the location-shift and location-scale-shift scenarios. These results are not surprising because our $\operatorname{\textit{ERFE}}$ estimator inherits the properties of the within-estimator which is simply an $\operatorname{\textit{ERFE}}$ estimator of level $\tau=0.5.$ The real data application showed that some parameter estimates vary according to the asymmetric points signaling the presence of heteroscedasticity in the data. Therefore we need to go beyond the mean regression to capture unobserved heterogeneity of the data and provide an overview of the relationship between the regressors and the dependent variables for a better decision making.
Our $\operatorname{\textit{ERFE}}$ model suffers from the same limitations as the $\operatorname{\textit{FE}}$ model which corresponds to the $\operatorname{\textit{ERFE}}$ model of level $\tau = 0.5.$ The $\operatorname{\textit{ERFE}}$ model estimates only the effects of the time-variant regressors. The $\operatorname{\textit{ERFE}}$ model ignores also the between-subject variations which can affect the efficiency of its standard error. The $\operatorname{\textit{ERFE}}$ model is a weighted mean regression and as such it is sensitive to aberrant values. Fortunately, there is a large number of regression diagnostic tools available to mitigate their influence.
There are alternatives in the literature that have been proposed to circumvent the lack of inference for the time-invariant regressors CornwellRupert1988, BaltagiKhantiAkom1990, while keeping the favorable properties of the $\operatorname{\textit{FE}}$ model. Future research should investigate the possibility of adapting these methods to the $\operatorname{\textit{ERFE}}$ framework.
In addition to the research avenues mentioned above we are currently exploring different alternatives such as penalizing the individual fixed-effects parameter to solve the incidental parameter problem while allowing inference on the time-invariant regressors.