GLS under Monotone Heteroskedasticity

GLS under monotone heteroskedasticity

\address{School of Social Sciences, Waseda University, 1-6-1 Nishiwaseda, Shinjuku-ku, Tokyo 169-8050, Japan.} \email{[email removed]} \address{Department of Economics, London School of Economics, Houghton Street, London, WC2A 2AE, UK.} \email{[email removed]} \address{Department of Economics, University of Mannheim, L7 3-5, 68161, Mannheim, Germany. } \email{[email removed]}

Introduction

The generalized least square (GLS) is one of the most basic tools in regression analyses. It yields the best linear unbiased estimator in the classical linear regression model, and has been studied extensively in econometrics and statistics literature; see e.g., Wooldridge (2010, Chapter 7) for a review. A major issue in implementing the GLS is that the optimal weights given by the conditional error variance function (say, $\sigma^{2}(\cdot)$) are typically unknown to researchers and need to be estimated. One way to estimate $\sigma^{2}(\cdot)$ is to specify its parametric functional form and estimate it by a parametric regression for the squared OLS residuals of the original regression on the specified covariates. However, economic theory rarely provides exact functional forms of $\sigma^{2}(\cdot)$, and the feasible GLS using misspecified $\sigma^{2}(\cdot)$ is no longer asymptotically efficient (Cragg, 1983). To address this issue, Carroll (1982) and Robinson (1987) proposed to estimate $\sigma^{2}(\cdot)$ nonparametrically and established the asymptotic equivalence of the resulting feasible GLS estimator with the infeasible one under certain regularity conditions. This is a remarkable result, but it requires theoretically and practically judicious choices of smoothing parameters, such as bandwidths, series lengths, or numbers of neighbors. It should be noted that such smoothing parameters appear in not only the point estimator but also its standard error for inference, and their choices typically require some assumption or knowledge of the smoothness of the conditional variance and associated density functions, such as their differentiability orders.

In this paper, we propose an alternative approach to estimate the conditional error variance function to implement the GLS by exploring a shape constraint of $\sigma^{2}(\cdot)$ instead of its smoothness as in Robinson (1987). As argued by Matzkin (1994), economic theory often provides shape constraints for functions of economic variables, such as monotonicity, concavity, or symmetry. In particular, we focus on situations where $\sigma^{2}(\cdot)$ is known to be monotone in its argument even though its exact functional form is unspecified, and propose to estimate $\sigma^{2}(\cdot)$ by utilizing the method of isotonic regression (see a review by Groeneboom and Jongbloed, 2014). It is known that the conventional isotonic regression estimator typically yields piecewise constant function estimates and does not involve any tuning parameters. Although the limiting behavior of the isotonic regression estimator is less tractable (such as the $n^{1/3}$-consistency and complicated limiting distribution), we show that our feasible GLS estimator using the optimal weights by the isotonic estimator with some trimming for boundary observations is asymptotically equivalent to the infeasible GLS estimator. Furthermore, we can plug in this isotonic estimator to estimate the asymptotic variance of the GLS estimator for statistical inference.

For the linear model $Y=X^{\prime}\beta+U$ in the presence of heteroskedasticity $\sigma^{2}(X)=E[U^{2}|X]$, using feasible GLS to improve the estimation efficiency has a long history. On the one hand, several parametric models have been proposed to estimate conditional error variance function $\sigma^{2}(\cdot)$. See Remark (ref) below. On the other hand, Carroll (1982) and Robinson (1987) estimated $\sigma^{2}(\cdot)$ with kernel and nearest neighbor estimator, respectively, and they showed their semiparametric GLS estimators are asymptotically equivalent to the infeasible GLS estimator and thus efficient. Compared to existing parametric methods, our proposed method imposes monotonicity, a feature implied by many parametric models, but it is nonparametric and does not rely on any specific parametric function form.\footnote{Monotone heteroskedasticity is often observed in economic literature. For example, Mincer (1974) argued that the variance of wages, when conditioned on education, should increase with the level of education because individuals with higher education have a broader array of job choices. Ruud (2000) cited this argument and provided empirical evidence in his Figure 18.1 based on the CPS data from March 1995. Another example can be found in Example 8.6 of Wooldridge (2013, pp. 283-284), where he employed a univariate conditional variance function of log income to explain the heteroskedasticity observed in net total financial wealth of people in the United States.} Compared to existing nonparametric methods, our proposed method involves only some tuning to trim boundary observations which does not require knowledge of the smoothness of the conditional variance and associated density functions. In the Monte Carlo simulations, we show that our proposed method outperforms the above-mentioned nonparametric methods at almost every choice of smoothing parameters, while it performs as well as parametric feasible GLS estimators with correctly specified conditional error variance function.

The isotonic estimator can date back to the middle of the last century. Earlier work includes Ayer et al. (1955), Grenander (1956), Rao (1969, 1970), and Barlow and Brunk (1972), among others. The isotonic estimator of a regression function can be formulated as a least square estimation with monotonicity constraints. Suppose that the conditional expectation $E[Y|X]=m(X)$ is monotone increasing, for an iid random sample $\{Y_{i},X_{i}\}_{i=1}^{n}$, the isotonic estimator is the minimizer of the sum of squared errors, $\underset{m\in\mathcal{M}}{\min}\sum_{i=1}^{n}\{Y_{i}-m(X_{i})\}^{2},$ where $\mathcal{M}$ is the class of monotone increasing functions. The minimizer can be calculated with the pool adjacent violators algorithm (Barlow and Brunk, 1972), or equivalently by solving the greatest convex minorant of the cumulative sum diagram $\{(0,0),(i,\sum_{j=1}^{i}Y_{j}),i=1,\ldots,n\}$, where the corresponding $\{X_{i}\}_{i=1}^{n}$ are ordered sequence; see Groeneboom and Jongbloed (2014) for a comprehensive discussion of different aspects of isotonic regression. Moreover, recent developments in the monotone single index model provide convenient and flexible tools for combining monotonicity and multi-dimensional covariates. In a monotone single index model, the conditional mean of $Y$ is modeled as $E[Y|X]=m(X^{\prime}\alpha)$, and the monotone link function $m(\cdot)$ is solved with isotonic regression. Balabdaoui, Durot and Jankowski (2019) studied the monotone single index model with the monotone least square method. Groeneboom and Hendrickx (2018), Balabdaoui, Groeneboom and Hendrickx (2019), and Balabdaoui and Groeneboom (2021) developed a score-type approach for the monotone single index model. Their approach can estimate the single index parameter $\alpha$ and the link function $m(\cdot)$ at $n^{-1/2}$-rate and $n^{-1/3}$-rate respectively. We employ their approach for the estimation of the conditional variance function in the multivariate case. Recently, Babii and Kumar (2023) applied the isotonic regression to their analysis of regression discontinuity designs. To this end, Babii and Kumar (2023) extended existing results concerning the boundary properties of Grenander's estimator (e.g., those from Woodroofe and Sun, 1993, and Kulikov and Lopuha�, 2006) to derive the asymptotic distribution of their trimmed isotonic regression discontinuity estimator. To regularize the isotonic estimator in the weights of our proposed GLS estimator, we employ a similar trimming strategy while adapting the theory of Babii and Kumar (2023) to our context of the conditional variance estimation. We contribute to this literature on isotonic regression by showing that the isotonic estimates can be employed as first stage estimates to be plugged in for semiparametric objects. Furthermore, we note that our isotonic estimator involves generated variables (i.e., OLS residuals), which make theoretical developments substantially different from the existing ones.

This paper is organized as follows. In Section (ref), we consider the case where $\sigma^{2}(\cdot)$ is monotone in one covariate, present our GLS estimator, and study its asymptotic properties. Section (ref) extends our GLS approach to the case where $\sigma^{2}(\cdot)$ is specified by a monotone single index function. Section (ref) illustrates the proposed method by a simulation study and empirical example.

Heteroskedasticity by univariate covariate

We first consider the case where monotone heteroskedasticity is caused by a single covariate. In particular, consider the following multiple linear regression model

where $X\in\mathscr{X}=[x_{L},x_{U}]$ is a scalar covariate with compact support and $Z$ is a vector of other covariates. In this section, we focus on the case where heteroskedasticity is caused by the covariate $X$, i.e.,

and $\sigma^{2}(\cdot)$ is a monotone increasing function. The case of monotone decreasing $\sigma^{2}(\cdot)$ is analyzed analogously (by setting $U^{2}$ as $-U^{2}$). In the setup ((ref)), we assume that the researcher knows which covariate should be included in $\sigma^{2}(\cdot)$ based on economic theory or other prior information. This setup should be considered as a useful benchmark to provide a clear exposition of the main concept and the asymptotic properties of the proposed monotone GLS estimator. Without the covariates $Z$, the above model covers a bivariate regression model, and our approach is new even in such a fundamental setup. Furthermore, this setup covers the case where $X$ contained in ((ref)) does not enter the regression model ((ref)) by setting $\beta=0$ (such a situation is considered in our empirical illustration in Section (ref)). Extensions to relax the assumption in ((ref)) will be discussed in Remark (ref) and Section (ref).

Let $\theta=(\alpha,\beta,\gamma^{\prime})^{\prime}$ be a vector of the slope parameters and $W:=(1,X,Z^{\prime})^{\prime}$ so that the model in ((ref)) can be written as $Y=W^{\prime}\theta+U$. Based on an iid sample $\{Y_{i},X_{i},Z_{i}\}_{i=1}^{n}$, the infeasible GLS estimator for $\theta$ is written as

where $\sigma_{i}^{2}=\sigma^{2}(X_{i})$. In order to make this estimator feasible, various approaches have been proposed in the literature.

In this paper, we are concerned with the situation where the researcher knows $\sigma^{2}(\cdot)$ is monotone in a particular regressor $X$ but its exact functional form is unspecified. In particular, by utilizing knowledge of the monotonicity of $\sigma^{2}(\cdot)$, we propose to estimate $\sigma^{2}(\cdot)$ by the isotonic regression from the squared OLS residual on the regressor $X$. More precisely, let $\hat{\theta}_{\mathrm{OLS}}=\left(\sum_{i=1}^{n}W_{i}W_{i}^{\prime}\right)^{-1}\left(\sum_{i=1}^{n}W_{i}Y_{i}\right)$ be the OLS estimator for ((ref)), and $\hat{U}_{j}=Y_{j}-W_{j}^{\prime}\hat{\theta}_{\mathrm{OLS}}$ be its residual. Then we estimate $\sigma^{2}(\cdot)$ by

Although this estimator is shown to be consistent for $\sigma^{2}(\cdot)$ in the interior of support $[x_{L},x_{U}]$ of $X$, it is generally biased at the lower boundary $x_{L}$, which may cause inconsistency of the resulting GLS estimator. Therefore, we propose to trim observations whose $X_{i}$'s are too close to $x_{L}$, and develop the following feasible GLS estimator

where $\mathbb{I}\{\cdot\}$ is the indicator function, and the trimming term $q_{n}$ is set as the $(n^{-1/3})$-th sample quantile of $\{X_{i}\}_{i=1}^{n}$.

Let $\mathcal{B}(a,R)$ be a ball around $a$ with radius $R$; for $\varepsilon=U^{2}-\sigma^{2}(X)$, define $\sigma_{\varepsilon}^{2}(x)=E[\varepsilon^{2}|X=x]$. To study the asymptotic properties of the proposed estimator $\hat{\theta}$, we impose the following assumptions.

Assumption A1 is standard. As pointed out in Balabdaoui, Groeneboom and Hendrickx (2019, p.13), the compact support assumption can be relaxed when $X$ follows a sub-Gaussian distribution. In this case, the $L^{2}$-convergence rate of the isotonic estimator will decrease from $O_{p}(n^{-1/3}\log n)$ to $O_{p}(n^{-1/3}(\log n)^{5/4})$. Another impact of relaxing the distribution of $X$ (and $Z$) to a sub-Gaussian one is on the concentration rate of $\max_{j}|\hat{U}_{j}^{2}-U_{j}^{2}|$ (see Appendix (ref) for more details). This rate, used in proving Lemma 1 and explaining the concentration of $T_{1}$ and $T_{2}$ in Appendix (ref), will inflate by a factor of $\log n$. However, even with this change, we still have $\max_{j}|\hat{U}_{j}^{2}-U_{j}^{2}|=o_{p}(n^{-1/3})$, which is the key to show that the impact of substituting infeasible $U^{2}$ with estimated $\hat{U}^{2}$ on isotonic estimators is asymptotically negligible. Considering that the convergence rates of these aforementioned terms are slowed down by a factor of $\log n$ at most, the validity of the main results in this paper is preserved with sub-Gaussian covariates, but the analytical derivation would become more cumbersome. For a clearer and more concise exposition, we maintain the compact support assumption on $X$. Assumption A2 is on the error term. The monotonicity of $\sigma^{2}(\cdot)$ is the main assumption. The assumption on arbitrary higher moments, which rules out some fat-tailed distributions, is commonly used to obtain some maximal inequalities (cf. van der Vaart and Wellner, 1996, Lemma 2.2.11, for a similar assumption). Assumption A3 contains additional mild conditions on the density of $X$.

We first present asymptotic properties of the conditional error variance estimator $\hat{\sigma}^{2}(\cdot)$ in ((ref)). Let $q_{n}^{*}$ be the $(n^{-1/3})$-th population quantile of $X$, $D_{A}^{L}[f](a)$ be the left derivative of the greatest convex minorant of a function $f(\cdot)$ evaluated at $a\in A$, and $\{\mathcal{W}_{t}\}$ be the standard Brownian motion. Also define $c^{*}=\lim_{n\to\infty}n^{1/3}(q_{n}^{*}-x_{L})$. Assumption A3 guarantees $0<c^{*}<\infty$. Then we obtain the following lemma for the behavior of $\hat{\sigma}^{2}(\cdot)$ around the boundary $x_{L}$, which extends the result by Babii and Kumar (2023, Theorem 2.1(ii)) by allowing the generated variable $\hat{U}_{i}^{2}$ as a regressand for $\hat{\sigma}^{2}(\cdot)$.

Based on this lemma, the asymptotic distribution of our feasible GLS estimator $\hat{\theta}$ is obtained as follows.

This theorem implies that our estimator $\hat{\theta}$ has the same limiting distribution as the infeasible GLS estimator $\hat{\theta}_{\mathrm{IGLS}}$ and thus achieves the semiparametric efficiency bound. This result extends the scope of the isotonic regression method by showing that the isotonic estimates, possibly with generated variables, can be employed as first stage estimates to be plugged in for semiparametric objects. We re-emphasize that $\hat{\theta}$ involves only a trimming term $q_{n}$, the $(n^{-1/3})$-th sample quantile of $\{X_{i}\}_{i=1}^{n}$.\footnote{Although our estimator $\hat{\theta}$ in ((ref)) does not involve any tuning constant, the trimming term $q_{n}$ should be understood as the $c\cdot(n^{-1/3})$-th sample quantile of $\{X_{i}\}_{i=1}^{n}$, where the tuning constant is set as $c=1$. Indeed Theorem (ref) holds true with any $c>0$. If we compare with other nonparametric methods, smoothing parameters, such as bandwidths, series lengths, and neighbors, typically require two constants to implement. For example, for the bandwidth parameter $b=c_{1}n^{-c_{2}}$, researchers need to choose $c_{1}$ and $c_{2}$. The constant $c_{1}$, which is analogous to $c$ above, can be any positive number. However, they also need to choose a positive constant $c_{2}$ whose upper bound typically depends on (unknown) smoothness of underlying functions.}

Heteroskedasticity by multivariate covariates

We now consider the model

where $X$ is a vector of covariates. This section focuses on the case where heteroskedasticity takes the form of a monotone single index function of $X$ with unknown parameters $\eta_{0}$, i.e., $E[U^{2}|X,Z]=E[U^{2}|X]=\sigma^{2}(X^{\prime}\eta_{0})$ for a monotone increasing function $\sigma^{2}(\cdot)$. Single index models are known to be more flexible than parametric models and achieve dimension reduction relative to nonparametric models.

For identification, $\eta_{0}$ is normalized as $||\eta_{0}||=1$. Define

We show in Lemma (ref) that $\sigma{}^{2}(\cdot)$ and $\eta_{0}$ can be consistently estimated by extending the method proposed in Balabdaoui, Groeneboom and Hendrickx (2019) (BGH hereafter) and Balabdaoui and Groeneboom (2021) to allow generated variables. In particular, for a given $\eta$, define the isotonic regression of $\{\hat{U}_{i}^{2}\}_{i=1}^{n}$ on $\{X_{i}^{\prime}\eta\}_{i=1}^{n}$ as

where $\mathcal{M}$ is the set of monotone increasing functions defined on $\mathbb{R}$. Based on this, $\hat{\eta}$ can be estimated by minimizing the square sum of a score function. For example, the simple score estimator in the spirit of BGH and Balabdaoui and Groeneboom (2021) is given by

where $\left\Vert \cdot\right\Vert $ is the Euclidean norm: $\left\Vert a\right\Vert =\sqrt{\sum_{j=1}^{k}a_{j}^{2}}$ for $a=(a_{1},\dots,a_{k})^{\prime}\in\mathbb{R}^{k}$.

Letting $\hat{\sigma}_{i}^{2}=\hat{\sigma}_{\hat{\eta}}^{2}(X_{i}^{\prime}\hat{\eta})$ and $W=(1,X^{\prime},Z^{\prime})^{\prime}$, we propose the following GLS estimator

where $q_{n}$ is the $(n^{-1/3})$-th sample quantile of $\{X_{i}^{\prime}\hat{\eta}\}_{i=1}^{n}$.

To avoid unnecessarily heavy notations, in the multivariate case, we redefine some notations, which have similar meanings to those used in Section (ref). Define $\varepsilon=U^{2}-\sigma^{2}(X^{\prime}\eta_{0})$, $\sigma_{\varepsilon}^{2}(\cdot)=E[\varepsilon^{2}|X^{\prime}\eta_{0}=\cdot]$, $x_{L}=\inf_{x\in\mathscr{X}}(x^{\prime}\eta_{0})$, and $x_{U}=\sup_{x\in\mathscr{X}}(x^{\prime}\eta_{0})$. Let $f_{X}(\cdot)$ be the density function of the random variable $X^{\prime}\eta_{0}$. Let $q_{n}^{*}$ be the $(n^{-1/3})$-th population quantile of $X^{\prime}\eta_{0}$, $q_{n}$ be the $(n^{-1/3})$-th sample quantile of $\{X_{i}^{\prime}\hat{\eta}\}_{i=1}^{n}$, $c^{*}=\lim_{n\to\infty}n^{1/3}(q_{n}^{*}-x_{L})$, and $D_{A}^{L}[f](a)$ be the left derivative of the greatest convex minorant of function $f(\cdot)$ evaluated at $a\in A$. Let $\text{dim}(w)$ be the dimension of a vector $w$.

Assumptions M1-M3 are analogs of Assumptions A1-A3, respectively. The main assumption is the monotonicity of $\sigma_{\eta}^{2}(\cdot)$. Assumptions M4-M6 are additional regularity conditions for the monotone index model. By Assumption M1, we have $-\infty<x_{L}<x_{U}<\infty$. Then similar to Lemma (ref), we obtain the following lemma for the behavior of $\hat{\sigma}_{\hat{\eta}}^{2}(\cdot)$ around $x_{L}$.

Based on this lemma, the asymptotic distribution of the GLS estimator $\hat{\theta}$ in ((ref)) is obtained as follows. Let $\sigma_{i}^{2}=\sigma^{2}(X_{i}^{\prime}\eta_{0})$.

Similar comments to Theorem (ref) apply here. Our estimator $\hat{\theta}$ is asymptotically equivalent to the infeasible GLS estimator $\hat{\theta}_{\mathrm{IGLS}}$. In terms of technical contribution, our theoretical analysis generalizes existing ones in, e.g., Babii and Kumar (2023), BGH, and Balabdaoui and Groeneboom (2021) to accommodate generated variables. Similar to Remark (ref), even when the monotonicity assumption of $\sigma_{\eta}^{2}(\cdot)$ is violated, $\hat{\theta}$ is still consistent for $\theta$ and asymptotically normal at the $\sqrt{n}$-rate with certain robust asymptotic variance. Furthermore, endogenous regressors can be accommodated as in Remark (ref).

Numerical illustrations

Simulation

We now investigate the finite sample properties of the proposed GLS estimator by a Monte Carlo experiment. We follow the simulation design by Cragg (1983) and Newey (1993). The first data generating process, denoted by DGP1, is the heteroskedastic linear model with a univariate covariate and normally distributed disturbance:\footnote{Normal random variables are not compactly supported, and hence it violates Assumption A1. However, as discussed in the remark on Assumption A1, this assumption can be relaxed.}

We consider three sample sizes, $n=50$, 100, and 500. The number of replications is set to 1,000.

In addition to the feasible GLS estimator with monotone heteroskedasticity (MGLS), we consider the ordinary least squares (OLS), infeasible generalized least squares (GLS), feasible GLS (FGLS), and nearest neighbor estimators (k-NN). GLS requires knowledge of the conditional error variance function ((ref)), including the values of the coefficients. In contrast, FGLS proceeds with the known functional form, but the coefficients are estimated. The “k-NN automatic” chooses the number of neighbors by a cross-validation procedure suggested by Newey (1993). All the estimators except OLS are the weighted least squares estimators, and their differences come from how the weights are calculated. Following Newey (1993), we calculate the weights for each method by taking a ratio of the predicted squared residual to the estimated variance of the disturbance, censoring the result below $0.04$.

Table (ref) presents the simulation results for estimation. The first column shows the estimation methods, and the following two columns show the root mean-squared error (RMSE) and mean absolute error (MAE) for DGP1 with $n=50$. The results for GLS report the levels of the RMSE and MAE, and those for others are their ratios relative to GLS. The next two columns give the corresponding results with $n=100$ and the last two columns with $n=500$. Two rows for each estimator show the results for $\beta_{0}$ and $\beta_{1}$, respectively. The inefficiency and inaccuracy of OLS are apparent. FGLS performs quite well, and this is natural when the conditional error variance functions are correctly specified. The performance of k-NN varies with the choice of $k$ and is in between OLS and FGLS. We observe that the performance of MGLS is better than k-NN in every choice of smoothing parameters. The result of MGLS is comparable to that of FGLS if not better. MGLS's independence of a smoothing parameter is clearly desirable. We also note that MGLS performs well even for $n=50$.

The last four columns of Table (ref) present the results for DGP2 with a homoskedastic error:

For DGP2, all estimators work reasonably well although the performance of k-NN with $k=6$ is worse than others.

Next, we consider the heteroskedastic linear models with multivariate covariates, denoted by DGP3:

The conditional error variance function of DPG3 is of a monotone single index structure. Using the notation in ((ref)), DGP3 corresponds to the structure with $\sigma^{2}(a)=a^{2}$, $X^{\prime}=(X_{1},X_{2})$, and $\eta_{0}=(\sqrt{.2},\sqrt{.2})^{\prime}$. The left panel of Table (ref) shows the results of DGP3 in the same manner as Table (ref). For each estimation method, two rows show the results for $\beta_{0}$ and $\beta_{1}$, and those for $\beta_{2}$ are omitted to avoid redundancy. k-NNs and MGLS perform better than FGLS, and this is in contrast to the performance of DGP1. In general, MGLS works better than k-NNs except for a few cases.

To see the potential applicability of MGLS to a non-single index structure, we consider another heteroskedastic linear model denoted by DGP4:

The right panel of Table (ref) shows the results. The results for DGP 4 are overall similar to those of DGP3. An exception is FGLS, which performs poorly for DGP3. MGLS works remarkably well for the heteroskedasticity of a non-single index structure.

Next, we turn to the simulation results on inference. Tables (ref) and (ref) show empirical coverages (EC) and average lengths (AL) for the 95% confidence intervals under DGPs 1-4. Again we consider GLS, OLS, FGLS, k-NN, and MGLS. For OLS, three types of confidence intervals are considered. They are based on the usual OLS standard error (OLS-U), the heteroskedasticity-robust standard error (OLS-R), and the wild bootstrap standard error (OLS-boot). For MGLS, we also present the results for its robust version. We observe that the empirical coverages are smaller than the nominal coverage 0.95 for all DGPs and all methods except GLS. It is natural that OLS-U performs poorly since it is invalid except for DGP2. The performance of k-NN is worse than others for all DGPSs in terms of empirical coverage. OLS-R, OLS-boot, FGLS, and MGLS work similarly in terms of empirical coverage, however, we note that the average length of OLS-R is much larger than those of FGLS and MGLS except for DGP2. While the empirical coverages of OLS-Boot are similar to those of OLS-R, the average lengths of OLS-Boot are smaller than those of OLS-R but still larger than those of MGLS. MGLS works quite well for all DGPs, especially for $n=500$. The results of MGLS (Robust) are similar to those of MGLS especially when $n=100$ and 500. Finally, we note that the empirical coverages tend to be lower when $n=50$ than when $n=100$ and $500$. Careful interpretation of results is recommended when the sample size is small.

Empirical example

We illustrate how the proposed method in this paper can improve the precision of the traditional OLS approach. In doing so, we revisit Acemoglu and Restrepo (2017) that investigate the relationship between an aging population and economic growth. After Hansen (1939), a popular perspective is that countries undergoing faster aging suffer more economically partly because of excessive savings by an aging population. In contrast to the perspective, Acemoglu and Restrepo (2017) find no evidence of a negative relationship between aging and GDP per capita after controlling for initial GDP per capita, initial demographic composition, and trends by region.

Acemoglu and Restrepo (2017) estimated eight specifications for the regression of the change in (log) GDP per capita from 1990 to 2015 (denoted by GDP) on the population aging measured by the change in the ratio of the population above 50 to those between the ages of 20 and 49 (denoted by Aging). The results are reproduced in Panel A of Table (ref). Those in columns 1-5 are based on the sample including 169 countries. Column 1 shows the result of the simple regression. Standard errors robust to heteroskedasticity are reported in square brackets. Column 2 shows the result with an additional regressor, the initial log GDP per worker in 1990. Column 3 in addition includes the initial demographic information, the ratio of the population above 50 to those between 20 and 49 in 1990 (denoted by Initial Ratio), and the population in 1990. Column 4 additionally uses dummies for seven regions, Latin America, East Asia, South Asia, Africa, North Africa and Middle East, Eastern Europe and Central Asia, and Developed Countries. Column 5 estimates the same specification as Column 4 with instruments of birthrates for the 1960, 1965, 1970, 1975, and 1980 cohorts. Columns 6 to 8 report the result for OECD countries using specifications of Columns 1, 3, and 5, respectively. The number of observations for the first five columns is 169, and that for the last three columns is 35. Seven out of eight OLS estimates indicate positive relationships and five of them are statistically significant at the 5 percent level. Acemoglu and Restrepo (2017) discuss that these findings can be explained by the adoption of automation technologies based on a theoretical model.

We estimate the same specifications by MGLS proposed in this paper. Acemoglu and Restrepo (2017) show that the negative effect of aging can be mitigated or reversed by adopting new automation technologies given abundant capital. This also implies that the effect of aging can be negative without sufficient capital. Hence it would be reasonable to consider Aging as a source of heteroskedasticity. The upper panel of Figure (ref) shows the relationship between the residual from the simple regression of column 1 in Panel A and Aging. Heteroskedasticity due to Aging is not easily confirmed visually. We consider Initial Ratio as another source of heteroskedasticity since the low ratio of old to young in 1990 is likely correlated with more aging in 2015, leading to larger variability in GDP per capita by the same reasoning discussed above. The lower panel of Figure (ref) presents the relationship between the residual from the simple regression of column 1 in Panel A and Initial Ratio, and we see that the variability decreases with the growing ratio.

Panels B, C, and D of Table (ref) show the results of MGLS. Panels B and C present the results for cases where the conditional error variance functions depend on Aging and Initial Ratio, respectively. Panel D reports the results where the conditional error variance functions depend on all exogenous regressors except the regional dummies. Standard errors based on Theorems 1-2 and their analogous versions for IV estimators are reported in parentheses, while robust standard errors are reported in square brackets. First, we observe reductions in standard errors for all MGLSs relative to OLS. The differences stand out when $n=169$. Second, the two standard errors are similar for the MGLS estimates under exogeneity while they differ for the IV estimates. These are the supporting evidence for the monotone specification of the conditional error variance function for the MGLS method with exogenous regressors but not for the IV method. Third, the results given in Columns 2, 3, and 4 are stable, while the results of IV estimates and OECD countries contain a lot of variations. Those variations can be due to non-monotone conditional error variance functions and/or small sample sizes, and further investigations will be required. Overall, the standard errors of MGLS tend to be smaller or no larger than those of OLS, which demonstrates the increased precision of MGLS.