GLS UNDER MONOTONE HETEROSKEDASTICITYYOICHI ARAI, TAISUKE OTSU, AND MENGSHAN XU
Abstract. The generalized least square (GLS) is one of the most basic tools in regression
analyses. A major issue in implementing the GLS is estimation of the conditional variance

arXiv:2210.13843v2 [econ.EM] 23 Jan 2024

function of the error term, which typically requires a restrictive functional form assumption forparametric estimation or smoothing parameters for nonparametric estimation. In this paper, we
propose an alternative approach to estimate the conditional variance function under nonparametric monotonicity constraints by utilizing the isotonic regression method. Our GLS estimatoris shown to be asymptotically equivalent to the infeasible GLS estimator with knowledge of theconditional error variance, and involves only some tuning to trim boundary observations, notonly for point estimation but also for interval estimation or hypothesis testing. Our analysis
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 infor semiparametric objects. Simulation studies illustrate excellent finite sample performancesof the proposed method. As an empirical example, we revisit Acemoglu and Restrepo’s (2017)
study on the relationship between an aging population and economic growth to illustrate how
our GLS estimator effectively reduces estimation errors.

1. 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 beenstudied extensively in econometrics and statistics literature; see e.g., Wooldridge (2010, Chapter7) for a review. A major issue in implementing the GLS is that the optimal weights given by theconditional error variance function (say, σ 2 (·)) are typically unknown to researchers and need
to be estimated. One way to estimate σ 2 (·) is to specify its parametric functional form andestimate 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 formsof σ 2 (·), and the feasible GLS using misspecified σ 2 (·) is no longer asymptotically efficient(Cragg, 1983). To address this issue, Carroll (1982) and Robinson (1987) proposed to estimate
σ 2 (·) 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 theconditional variance and associated density functions, such as their differentiability orders.
In this paper, we propose an alternative approach to estimate the conditional error variancefunction to implement the GLS by exploring a shape constraint of σ 2 (·) instead of its smoothness1

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 σ 2 (·) is known to be monotone in its argument
even though its exact functional form is unspecified, and propose to estimate σ 2 (·) 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 constantfunction estimates and does not involve any tuning parameters. Although the limiting behaviorof the isotonic regression estimator is less tractable (such as the n1/3 -consistency and complicatedlimiting distribution), we show that our feasible GLS estimator using the optimal weights by theisotonic 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 ′ β + U in the presence of heteroskedasticity σ 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
σ 2 (·). See Remark 5 below. On the other hand, Carroll (1982) and Robinson (1987) estimatedσ 2 (·) 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.1 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 theMonte Carlo simulations, we show that our proposed method outperforms the above-mentioned
nonparametric methods at almost every choice of smoothing parameters, while it performs aswell as parametric feasible GLS estimators with correctly specified conditional error variancefunction.
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 squareestimation with monotonicity constraints. Suppose that the conditional expectation E[Y |X] =m(X) is monotone increasing, for an iid random sample {Yi , Xi }ni=1 , the isotonic estimator isP
the minimizer of the sum of squared errors, min ni=1 {Yi − m(Xi )}2 , where M is the class ofm∈Mmonotone 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 ofP
the cumulative sum diagram {(0, 0), (i, ij=1 Yj ), i = 1, . . . , n}, where the corresponding {Xi }ni=11Monotone heteroskedasticity is often observed in economic literature. For example, Mincer (1974) argued that thevariance 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 empiricalevidence 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.
2

are ordered sequence; see Groeneboom and Jongbloed (2014) for a comprehensive discussionof different aspects of isotonic regression. Moreover, recent developments in the monotone
single index model provide convenient and flexible tools for combining monotonicity and multidimensional covariates. In a monotone single index model, the conditional mean of Y is modeledas E[Y |X] = m(X ′ α), and the monotone link function m(·) is solved with isotonic regression.

Balabdaoui, Durot and Jankowski (2019) studied the monotone single index model with themonotone least square method. Groeneboom and Hendrickx (2018), Balabdaoui, Groeneboom
and Hendrickx (2019), and Balabdaoui and Groeneboom (2021) developed a score-type approachfor the monotone single index model. Their approach can estimate the single index parameter α
and the link function m(·) at n−1/2 -rate and n−1/3 -rate respectively. We employ their approachfor 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 Lopuhaä, 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 andKumar (2023) to our context of the conditional variance estimation. We contribute to thisliterature 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 theoreticaldevelopments substantially different from the existing ones.
This paper is organized as follows. In Section 2, we consider the case where σ 2 (·) is monotone
in one covariate, present our GLS estimator, and study its asymptotic properties. Section
3 extends our GLS approach to the case where σ 2 (·) is specified by a monotone single index
function. Section 4 illustrates the proposed method by a simulation study and empirical example.
2. 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 modelY = α + βX + Z ′ γ + U,

E[U |X, Z] = 0,

(2.1)

where X ∈ X = [xL , xU ] is a scalar covariate with compact support and Z is a vector ofother covariates. In this section, we focus on the case where heteroskedasticity is caused by thecovariate X, i.e.,E[U 2 |X, Z] = E[U 2 |X] =: σ 2 (X),

(2.2)

and σ 2 (·) is a monotone increasing function. The case of monotone decreasing σ 2 (·) is analyzed
analogously (by setting U 2 as −U 2 ). In the setup (2.2), we assume that the researcher knows
which covariate should be included in σ 2 (·) based on economic theory or other prior information.

This setup should be considered as a useful benchmark to provide a clear exposition of the mainconcept and the asymptotic properties of the proposed monotone GLS estimator. Without the3covariates Z, the above model covers a bivariate regression model, and our approach is new evenin such a fundamental setup. Furthermore, this setup covers the case where X contained in(2.2) does not enter the regression model (2.1) by setting β = 0 (such a situation is consideredin our empirical illustration in Section 4.2). Extensions to relax the assumption in (2.2) will be
discussed in Remark 1 and Section 3.
Let θ = (α, β, γ ′ )′ be a vector of the slope parameters and W := (1, X, Z ′ )′ so that the modelin (2.1) can be written as Y = W ′ θ + U . Based on an iid sample {Yi , Xi , Zi }ni=1 , the infeasible

GLS estimator for θ is written asθ̂IGLS =nX

σi−2 Wi Wi′

i=1

!−1

n
X
i=1

σi−2 Wi Yi

!

,

(2.3)

where σi2 = σ 2 (Xi ). 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 σ 2 (·) ismonotone in a particular regressor X but its exact functional form is unspecified. In particular, by utilizing knowledge of the monotonicity of σ 2 (·), we propose to estimate σ 2 (·) by
the isotonic regression from the squared OLS residual on the regressor X. More precisely, letPPθ̂OLS = ( ni=1 Wi Wi′ )−1 ( ni=1 Wi Yi ) be the OLS estimator for (2.1), and Ûj = Yj − Wj′ θ̂OLS be

its residual. Then we estimate σ 2 (·) byσ̂ 2 (·) = isotonic regression function from {Ûj2 }nj=1 on {Xj }nj=1 .

(2.4)

Although this estimator is shown to be consistent for σ 2 (·) in the interior of support [xL , xU ]of X, it is generally biased at the lower boundary xL , which may cause inconsistency of theresulting GLS estimator. Therefore, we propose to trim observations whose Xi ’s are too close
to xL , and develop the following feasible GLS estimator
!
!−1 n
n
X
X
−2
−2
′
I{Xi ≥ qn }σ̂i Wi Yi ,I{Xi ≥ qn }σ̂i Wi Wi
θ̂ =

(2.5)

i=1

i=1

where I{·} is the indicator function, and the trimming term qn is set as the (n−1/3 )-th sample
quantile of {Xi }ni=1 .

Let B(a, R) be a ball around a with radius R; for ε = U 2 − σ 2 (X), define σε2 (x) = E[ε2 |X =x]. To study the asymptotic properties of the proposed estimator θ̂, we impose the following
assumptions.
Assumption.
A1: {Yi , Xi , Zi }ni=1 is an iid sample of (Y, X, Z). The support of (X, Z) is convex with

non-empty interiors and is a subset of B(0, R) for some R > 0. The support of X is acompact interval X = [xL , xU ].

A2: σ 2 : X → R is a monotone increasing function defined on X , and 0 < σ 2 (xL ) <σ 2 (xU ) < ∞. There exist positive constants a0 and M such that E[|U |2s |X = x] ≤a0 s!M s−2 for all integers s ≥ 2 and x ∈ X . For some positive constant δ, σ 2 (·) is
continuously differentiable on (xL , xL + δ), and σε2 (·) is continuous on (xL , xL + δ).
4

A3: X has a continuous density function fX (·) on X , and there exists a positive constantb such that b < fX (x) < ∞ for all x ∈ X .
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 L2 -convergence rate of the isotonic estimator will decrease fromOp (n−1/3 log n) to Op (n−1/3 (log n)5/4 ). Another impact of relaxing the distribution of X (andZ) to a sub-Gaussian one is on the concentration rate of maxj |Ûj2 − Uj2 | (see Appendix A for

more details). This rate, used in proving Lemma 1 and explaining the concentration of T1 andT2 in Appendix A.2, will inflate by a factor of log n. However, even with this change, we still
have maxj |Ûj2 − Uj2 | = op (n−1/3 ), which is the key to show that the impact of substituting

infeasible U 2 with estimated Û 2 on isotonic estimators is asymptotically negligible. Consideringthat 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 theerror term. The monotonicity of σ 2 (·) is the main assumption. The assumption on arbitrary
higher moments, which rules out some fat-tailed distributions, is commonly used to obtainsome maximal inequalities (cf. van der Vaart and Wellner, 1996, Lemma 2.2.11, for a similarassumption). Assumption A3 contains additional mild conditions on the density of X.
We first present asymptotic properties of the conditional error variance estimator σ̂ 2 (·) in (2.4).
L [f ](a) be the left derivative of the greatest
Let qn∗ be the (n−1/3 )-th population quantile of X, DA

convex minorant of a function f (·) evaluated at a ∈ A, and {Wt } be the standard Brownian

motion. Also define c∗ = limn→∞ n1/3 (qn∗ − xL ). Assumption A3 guarantees 0 < c∗ < ∞. Then

we obtain the following lemma for the behavior of σ̂ 2 (·) around the boundary xL , which extends
the result by Babii and Kumar (2023, Theorem 2.1(ii)) by allowing the generated variable Ûi2
as a regressand for σ̂ 2 (·).
2

Lemma 1. Under Assumptions A1-A3 and limx↓xL dσdx(x) > 0, it holds
"s

#
 
2 (x)
2 (x )
dσ
1
σ
d
L
ε
L
Wt + limt2 − t (1).
c∗n1/3 {σ̂ 2 (qn ) − σ 2 (qn )} → D[0,∞)x↓xLc∗ fX (xL )
dx
2

(2.6)

Based on this lemma, the asymptotic distribution of our feasible GLS estimator θ̂ is obtained
as follows.
Theorem 1. Under Assumptions A1-A3, it holds√dn(θ̂ − θ) → N (0, E[σ −2 (X)W W ′ ]−1 ),and the asymptotic variance matrix is consistently estimated by−2
1 Pn′ −1 .
i=1 σ̂i Wi Win

This theorem implies that our estimator θ̂ has the same limiting distribution as the infeasible
GLS estimator θ̂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 with5

generated variables, can be employed as first stage estimates to be plugged in for semiparametricobjects. We re-emphasize that θ̂ involves only a trimming term qn , the (n−1/3 )-th sample quantileof {Xi }ni=1 .2
Remark 1. [Extensions of (2.2)] The benchmark setup E[U 2 |X, Z] = σ 2 (X) considered inthis section can be extended in various ways. First, an extension to a single index model(say, E[U 2 |X, Z] = σ 2 (Xηx + Z ′ ηz )) will be discussed in the next section. Second, the modelin (2.1)-(2.2) can be extended to the case where the conditional variance varies with discrete
covariates Z (or its subvector), say E[U 2 |X, Z = z] = σz2 (X) with monotone functions σz2 (·)

for z ∈ {z (1) , . . . , z (D) }. In this case, we can implement the isotonic regression for each groupcategorized by z, and construct the feasible GLS estimator in an analogous way as (2.5). Third,
our approach may be extended to the additive monotone heteroskedasticity, say E[U 2 |X, Z] =σx2 (X) + σz2 (Z) with monotone functions σx2 (·) and σz2 (·). Although formal analysis is beyond thescope of this paper, the results in Mammen and Yu (2007) suggest that the isotonic estimators
for additive functions converge at similar rates as the univariate case, and we conjecture thata similar result as Theorem 1 can be obtained. Finally, when the conditional error variancefunction is multiplicative, say E[U 2 |X, Z] = σx2 (X)σz2 (Z), and the researcher knows the form ofσz2 (·) (e.g., Z is household size and σz2 (Z) = Z 2 ), then our feasible GLS estimator can be applied
to observations reweighted by 1/σz (Z).

Remark 2. [Monotonicity testing] Monotonicity is an assumption that can be tested. For observable random variables (Y, X), several methods have been developed to test whether E[Y |X]

is monotone increasing in X; see, e.g., Ghosal, Sen and van der Vaart (2000), Hall and Heckman (2000), Dümbgen and Spokoiny (2001), Chetverikov (2019), and Hsu, Liu and Shi (2019),

among others. All these tests can be adapted for our case, testing the monotonicity of σ 2 (·) withgenerated {Ûj2 }nj=1 and observed {Xj }nj=1 . Since Assumptions A1-A2 and θ̂OLS − θ = Op (n−1/2 )imply Ûj2 − Uj2 = Op (n−1/2 log n) uniformly over j = 1, . . . , n, the critical values of these tests
can be adjusted accordingly to maintain a proper asymptotic size.

Remark 3. [Misspecification of E[U 2 |X, Z]] We want to note that even if the assumption in(2.2) is violated (e.g., E[U 2 |X, Z] varies with Z or E[U 2 |X, Z] = σ 2 (X) with non-monotone

σ 2 (·)), our feasible GLS estimator θ̂ in (2.5) is still consistent for θ due to E[U |X, Z] = 0, and√asymptotically normal at the n-rate with the limiting distribution√dn(θ̂ − θ) → N (0, E[ρ(X)−1 W W ′ ]−1 E[ρ(X)−2 E[U 2 |X, Z]W W ′ ]E[ρ(X)−1 W W ′ ]−1 ),where ρ(·) = arg minm∈M E[{U 2 − m(X)}2 ] for the class of monotone increasing functions M.
Since σ̂ 2 (·) can estimate ρ(·), then the asymptotic variance matrix can be consistently estimated

2Although our estimator θ̂ in (2.5) does not involve any tuning constant, the trimming term qn should be understood as the c · (n−1/3 )-th sample quantile of {Xi }n
i=1 , where the tuning constant is set as c = 1. Indeed Theorem
1 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 = c1 n−c2 , researchers need to choose c1 and c2 . The constant c1 , which is analogous to c
above, can be any positive number. However, they also need to choose a positive constant c2 whose upper boundtypically depends on (unknown) smoothness of underlying functions.

6

by
n

1 X −2σ̂ (Xi )Wi Wi′
n
i=1

!−1

n

1 X −4σ̂ (Xi )Ûi2 Wi Wi′
n
i=1

!

n

1 X −2σ̂ (Xi )Wi Wi′
n
i=1

!−1

.

(2.7)

This misspecification robust variance estimator is analogous to the one proposed by Cragg(1992) for the feasible GLS estimator with parametrically specified models for the conditional
error variance E[U 2 |X, Z].3
Remark 4. [Endogenous regressor] The result of Theorem 1 can also be extended to some linear
instrumental variable (IV) regression model. For notational simplicity, consider the following
univariate IV regression:
Y = α + βX + U,

E[U |Z] = 0,

where X is a scalar endogenous regressor and Z is a scalar IV, and we further assume E[X|Z] =η+γZ for some parameters (η, γ). This linearity assumption on E[X|Z] is not essential, and may
be relaxed by some nonparametric estimator of E[X|Z]. In this setup, the optimal instrument
for estimating (α, β)′ is given by (see, e.g., Newey, 1993)
!


1
0
∂(Y − α − βX)
E
Z E[U 2 |Z]−1 = −∂(α, β)′
η γ

1
Z

!

v −2 (Z),where v 2 (·) = E[U 2 |Z = ·]. Under the assumption of γ 6= 0 (i.e., the IV is relevant), theoptimal IV estimator is obtained by the method of moments estimator of the following moment
condition:
E

"

1
Z

!

v

−2

#

(Z)(Y − α − βX) = 0.

(2.8)

Under the monotonicity assumption of v 2 (·), we can obtain the isotonic estimator v̂ 2 (·) for v 2 (·)
by regressing the squared residuals ê2 = (Y − α̃ − β̃X)2 for an initial estimator (α̃, β̃) (e.g., thetwo-stage least squares estimator) on Z. The resulting estimator, v̂ 2 (·), should have the same
properties as those of σ̂ 2 (·) presented in Lemma 1, where qn is replaced with the (n−1/3 )-thsample quantile of {Zi }ni=1 . Based on this isotonic estimator, a feasible optimal IV estimatorθ̂IV = (α̂IV , β̂IV )′ is given byθ̂IV =nX
i=1

I{Zi ≥ qn }v̂

−2

(Zi )

1
Zi

!

!−1

(1, Xi )nX
i=1

I{Zi ≥ qn }v̂

−2

(Zi )

1
Zi

!

Yi

!

.

By applying the same arguments for Theorem 1, we can show that θ̂IV is asymptotically equivalent to the infeasible optimal IV estimator based on (2.8) with known v 2 (·).

3Based on simulation studies, Cragg (1992) recommended to use his misspecification robust variance estimatoreven when the parametric form of heteroskedasticity is correctly specified. Although a similar analysis is beyond
the scope of this paper, we also recommend to employ the variance estimator (2.7) in practice due to its consistency
regardless of the assumption in (2.2).
7

3. Heteroskedasticity by multivariate covariates
We now consider the modelY = α + X ′ β + Z ′ γ + U,

E[U |X, Z] = 0,

(3.1)

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 η0 , i.e.,E[U 2 |X, Z] = E[U 2 |X] = σ 2 (X ′ η0 ) for a monotone increasing function σ 2 (·). Single index

models are known to be more flexible than parametric models and achieve dimension reduction
relative to nonparametric models.
Remark 5. First, the monotone index model σ 2 (X ′ η0 ) covers several existing parametricmodels. Popular examples include σ 2 (X) = C(X ′ η0 )2−2λ (Box and Hill, 1974), σ 2 (X) =C exp(λ(X ′ η0 )) (Bickel, 1978), σ 2 (X) = C{1 + λ(X ′ η0 )2 } (Fuller, 1980) for some constants

C > 0 and λ; interestingly, all these parametric functions are monotone increasing (or de-creasing) in the index of X. Second, although the setup E[U 2 |X, Z] = σ 2 (X ′ η0 ) assumes thatthe researcher knows which (sub-)vector of covariates should be included in σ 2 (·), researchers
do not have to select those covariates in the case where such prior information is unavailable.
They can simply re-define the model in (3.1) without covariates Z (or equivalently specify as
E[U 2 |X, Z] = σ 2 (X ′ η0 + Z ′ ηz0 )). Our asymptotic theory below applies even if some covariates
are irrelevant for E[U 2 |X, Z].

For identification, η0 is normalized as ||η0 || = 1. Defineση2 (a) = E[σ 2 (X ′ η0 )|X ′ η = a].

(3.2)

We show in Lemma 4 that σ 2 (·) and η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 η, define theisotonic regression of {Ûi2 }ni=1 on {Xi′ η}ni=1 as1σ̂η2 = arg min
m∈M n

n
X
{Ûi2 − m(Xi′ η)}2 ,

(3.3)

i=1

where M is the set of monotone increasing functions defined on R. Based on this, η̂ 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 by2n

1X
Xi {Ûi2 − σ̂η2 (Xi′ η)} ,η̂ = arg min
η
n
i=1
qP
k
′
k
2
where k·k is the Euclidean norm: kak =j=1 aj for a = (a1 , . . . , ak ) ∈ R .

Letting σ̂i2 = σ̂η̂2 (Xi′ η̂) and W = (1, X ′ , Z ′ )′ , we propose the following GLS estimator
!
!−1 n
n
X
X
I{Xi′ η̂ ≥ qn }σ̂i−2 Wi Yi ,I{Xi′ η̂ ≥ qn }σ̂i−2 Wi Wi′
θ̂ =
i=1

i=1

8

(3.4)

(3.5)

where qn is the (n−1/3 )-th sample quantile of {Xi′ η̂}ni=1 .

To avoid unnecessarily heavy notations, in the multivariate case, we redefine some notations,which have similar meanings to those used in Section 2. Define ε = U 2 − σ 2 (X ′ η0 ), σε2 (·) =E[ε2 |X ′ η0 = ·], xL = inf x∈X (x′ η0 ), and xU = supx∈X (x′ η0 ). Let fX (·) be the density functionof the random variable X ′ η0 . Let qn∗ be the (n−1/3 )-th population quantile of X ′ η0 , qn be theL [f ](a) be the left(n−1/3 )-th sample quantile of {Xi′ η̂}ni=1 , c∗ = limn→∞ n1/3 (qn∗ − xL ), and DA

derivative of the greatest convex minorant of function f (·) evaluated at a ∈ A. Let dim(w) be
the dimension of a vector w.

Assumption.
M1: {Yi , Xi , Zi }ni=1 is an iid sample of (Y, X, Z). The support of (X, Z), X ×Z , is convex
with non-empty interiors and is a subset of B(0, R) for some R > 0.

M2: (i) There exists δ0 > 0 such that the function a 7→ ση2 (a) defined in (3.2) is monotone

increasing on Iη = {x′ η, x ∈ X } for each η ∈ B(η0 , δ0 ). (ii) 0 < inf a∈Iη ση2 (a) <supa∈Iη ση2 (a) < ∞ for each η ∈ B(η0 , δ0 ). (iii) There exist positive constants a0 andM such that E[|U |2s |X = x] ≤ a0 s!M s−2 for all integers s ≥ 2 and x ∈ X . (iv) ση2 (·)

is continuously differentiable on Iη for each η ∈ B(η0 , δ0 ). (v) σε2 (·) is continuous on(xL , xL + δ1 ) for some δ1 > 0.

M3: The random variable X ′ η0 has a density function fX (·) that is continuous on Iη0 .
There exists some real positive numbers b and b, such that 0 < b < fX (a) < b < ∞ holds

for all a ∈Iη0 .

M4: For each η ∈ B(η0 , δ0 ), the mapping a 7→ E[X|X ′ η = a] defined on Iη is bounded andhas a finite total variation.

M5: Cov[X ′ (η0 − η), σ 2 (X ′ η0 )|X ′ η] 6= 0 almost surely for each η 6= η0 .
R
2
M6: B := (x − E[X|x′ η0 ])(x − E[X|x′ η0 ])′ dσda(a)dP (x) has rank dim(η0 )−1.
′
a=x η0Assumptions M1-M3 are analogs of Assumptions A1-A3, respectively. The main assumption
is the monotonicity of ση2 (·). Assumptions M4-M6 are additional regularity conditions for themonotone index model. By Assumption M1, we have −∞ < xL < xU < ∞. Then similar toLemma 1, we obtain the following lemma for the behavior of σ̂η̂2 (·) around xL .
2

Lemma 2. Under Assumptions M1-M6 and lima↓xL dσda(a) > 0, it holds
"s

#
 
σε2 (xL )dσ 2 (a) ∗ 1 2
d
1/3
2
2
L
n {σ̂η̂ (qn ) − σ (qn )} → D[0,∞)Wt + limt − t (1).
ca↓xLc∗ fX (xL )da2Based on this lemma, the asymptotic distribution of the GLS estimator θ̂ in (3.5) is obtained
as follows. Let σi2 = σ 2 (Xi′ η0 ).
Theorem 2. Under Assumptions M1-M6, it holds√dn(θ̂ − θ) → N (0, E[σ −2 (X ′ η0 )W W ′ ]−1 ),−1P.
and the asymptotic variance matrix is consistently estimated by n1 ni=1 σ̂i−2 Wi Wi′
9

Similar comments to Theorem 1 apply here. Our estimator θ̂ is asymptotically equivalent
to the infeasible GLS estimator θ̂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 3, even when themonotonicity assumption of ση2 (·) is violated, θ̂ is still consistent for θ and asymptotically normal√
at the n-rate with certain robust asymptotic variance. Furthermore, endogenous regressors
can be accommodated as in Remark 4.
Remark 6. We can suggest two informal robustness checks for the monotone index assumptionin (3.2). One is to compute the standard errors robust to possible misspecification obtainedin the same manner as Remark 3 and compare them to those in Theorem 2. This can serve
as a robustness check for the monotone specification given variables of the conditional error
variance functions. Another is to report the results for the specification where all exogenous
variables are included to σ 2 (·) in addition to those for the chosen specifications. A large difference
between these results can be a sign of the misspecification of the chosen ones. See Section 4.2
for illustration.
Remark 7. In this section, we employ the monotone single index structure to model the multivariate conditional variance function. This strategy allows us to strike a balance between robustness and mitigating the curse of dimensionality. Indeed, the current specification can be extended
P
′
′ ′
′
′
to the multiple index model E[U 2 |X = x] = x′0 η0 + M
i=1 Gi (xi ηi ), for X = (X0 , X1 . . . , XM ) ,where {Gi (·)}M
i=1 are unknown monotone increasing functions. For the case of M = 1, this model

simplifies to a monotone partially linear single index model whose properties have been studied
by Xu and Otsu (2020). We are optimistic that, under certain regularity conditions, similar results as in this section can be obtained. To the best of our knowledge, we have not come across
any works that discuss the multiple monotone index model with M > 1 even for the conven-tional regression setup for E[Y |X = x]. A possible solution could be derived by combining theexisting literature on the monotone single index model (as cited in Section 1) with the literatureon the monotone additive model (for instance, Mammen and Yu, 2007). Another potential extension involves employing the nonparametric framework of Fang, Guntuboyina and Sen (2021)
to model the multivariate conditional variance function. This framework is free of parametricstructure, and it requires the true conditional variance to be entirely monotone increasing in its
arguments, i.e., σ 2 (x1 , z1 ) ≤ σ 2 (x2 , z2 ) if only if x1 ≤ x2 and z1 ≤ z2 . Explorations of these
extensions exceed the scope of this paper, and we leave them for future research.

4. Numerical illustrations
4.1. 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
10

with a univariate covariate and normally distributed disturbance:4
Yi = β0 + β1 Xi + ui ,β0 = β1 = 1,

ui = σi εi ,log(Xi ) ∼ N (0, 1),εi ∼ N (0, 1),
Xi and εi are independent,

σi2 = .1 + .2Xi + .3Xi2 .

(4.1)

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 (4.1), 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 ofthe disturbance, censoring the result below 0.04.
Table 4.1 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 thecorresponding results with n = 100 and the last two columns with n = 500. Two rows for eachestimator show the results for β0 and β1 , respectively. The inefficiency and inaccuracy of OLS
are apparent. FGLS performs quite well, and this is natural when the conditional error variancefunctions are correctly specified. The performance of k-NN varies with the choice of k and isin between OLS and FGLS. We observe that the performance of MGLS is better than k-NN inevery 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 notethat MGLS performs well even for n = 50.
The last four columns of Table 4.1 present the results for DGP2 with a homoskedastic error:
Yi = β0 + β1 Xi + ui ,β0 = β1 = 1,

ui ∼ N (0, 1),log(Xi ) ∼ N (0, 1),

Xi and ui are independent.

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

4Normal 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.
11

Table 4.1. Simulation: Estimation with univariate covariate

12

DGP1
DGP2
n = 50n = 100
n = 500
n = 50n = 100
n = 500
Estimator
RMSE MAE RMSE MAE RMSE MAE RMSE MAE RMSE MAE RMSE MAE
GLS (infeasible)
0.132
0.085 0.093 0.059 0.041 0.028 0.194 0.122 0.133 0.088 0.057 0.039
0.157
0.100 0.108 0.073 0.048 0.032 0.083 0.046 0.055 0.034 0.021 0.014
OLS
3.103
2.856 3.831 3.479 5.574 4.495 1.000 1.000 1.000 1.000 1.000 1.000
2.072
2.098 2.543 2.370 3.377 2.971 1.000 1.000 1.000 1.000 1.000 1.000
FGLS
1.279
1.210 1.245 1.233 1.598 1.152 1.032 1.041 1.026 1.033 1.024 1.067
1.427
1.268 1.406 1.280 1.271 1.242 1.090 1.092 1.075 1.036 1.088 1.090
k-NN (Automatic) 1.630
1.373 1.633 1.511 1.355 1.167 1.123 1.081 1.130 1.081 1.181 1.138
1.535
1.427 1.606 1.431 1.424 1.267 1.092 1.074 1.065 1.006 1.197 1.097
k-NN (k = 6)
1.554
1.361 1.525 1.498 1.474 1.417 1.274 1.243 1.253 1.155 1.359 1.276
1.466
1.421 1.472 1.462 1.454 1.459 1.178 1.143 1.177 1.114 1.350 1.344
k-NN (k = 15)1.600
1.386 1.566 1.365 1.251 1.108 1.037 1.076 1.079 1.059 1.081 1.140
1.520
1.398 1.546 1.408 1.247 1.197 1.003 1.046 1.037 1.012 1.066 1.053
k-NN (k = 24)
1.781
1.568 1.685 1.457 1.291 1.160 1.011 1.039 1.039 0.980 1.044 1.098
1.630
1.560 1.673 1.471 1.312 1.246 1.002 1.026 1.015 0.994 1.038 1.025
MGLS
1.379
1.285 1.326 1.279 1.113 1.129 1.039 1.091 1.049 1.075 1.027 1.075
1.327
1.214 1.332 1.249 1.113 1.144 1.043 1.051 1.051 1.058 1.055 1.066
Note: “RMSE” and “MAE” stand for the root mean squared error and mean absolute error, respectively. The results for GLS report
the levels of the RMSE and MAE, and those for others are their ratios relative to GLS.

Next, we consider the heteroskedastic linear models with multivariate covariates, denoted by
DGP3:
Yi = β0 + β1 X1i + β2 X2i + ui ,β0 = β1 = β2 = 1,

ui = σi εi ,εi ∼ N (0, 1),log(X1i ), log(X2i ) ∼ N (0, 1),

σi2 = .2(X1i + X2i )2 .

X1i , X2i and εi are independent,(4.2)

The conditional error variance function of DPG3 is of a monotone single index structure. Using
the notation in (3.2), DGP3 corresponds to the structure with σ 2 (a) = a2 , X ′ = (X1 , X2 ), and√ √η0 = ( .2, .2)′ . The left panel of Table 4.2 shows the results of DGP3 in the same manner as
Table 4.1. For each estimation method, two rows show the results for β0 and β1 , and those forβ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 fora few cases.
To see the potential applicability of MGLS to a non-single index structure, we consider anotherheteroskedastic linear model denoted by DGP4:
Yi = β0 + β1 X1i + β2 X2i + ui ,β0 = β1 = β2 = 1,

ui = σi εi ,εi ∼ N (0, 1),log(X1i ), log(X2i ) ∼ N (0, 1),

σi2 = .1 + .2X̃i + .3X̃i2 ,log(X̃i ) =X1i , X2i and εi are independent,log(X1i ) + log(X2i )
√
.
2

(4.3)

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

13

Table 4.2. Simulation: Estimation with multivariate covariates

14

DGP 3
DGP 4
n = 50n = 100
n = 500
n = 50n = 100
n = 500
Estimator
RMSE MAE RMSE MAE RMSE MAE RMSE MAE RMSE MAE RMSE MAE
GLS (infeasible)
0.162
0.103 0.110 0.071 0.045 0.028 0.165 0.108 0.115 0.076 0.049 0.033
0.163
0.107 0.109 0.072 0.048 0.033 0.108 0.067 0.071 0.046 0.029 0.020
OLS
3.401
3.589 4.255 4.168 6.650 6.695 3.069 2.653 3.792 2.980 4.897 4.186
1.942
1.950 2.317 2.260 3.051 2.610 2.318 2.170 2.914 2.338 3.809 3.198
FGLS
2.531
2.239 2.516 2.141 2.731 2.037 1.381 1.189 1.427 1.233 1.699 1.219
1.606
1.441 1.709 1.486 1.732 1.358 1.359 1.227 1.344 1.281 1.326 1.271
k-NN (Automatic) 1.952
1.925 2.108 1.709 1.786 1.537 1.868 1.638 1.778 1.488 1.709 1.390
1.546
1.429 1.680 1.489 1.516 1.318 1.766 1.771 1.865 1.763 2.009 1.626
k-NN (k = 6)
1.827
1.766 1.787 1.587 1.670 1.666 1.719 1.541 1.674 1.521 1.594 1.537
1.458
1.397 1.514 1.486 1.497 1.362 1.717 1.704 1.769 1.764 1.813 1.654
k-NN (k = 15)
1.914
1.957 1.850 1.669 1.490 1.385 1.769 1.611 1.669 1.491 1.373 1.246
1.468
1.401 1.511 1.428 1.313 1.248 1.712 1.727 1.769 1.639 1.588 1.517
k-NN (k = 24)
2.182
2.203 2.008 1.816 1.570 1.510 1.952 1.888 1.799 1.581 1.392 1.254
1.562
1.455 1.611 1.571 1.371 1.267 1.825 1.807 1.890 1.729 1.626 1.511
MGLS
2.144
1.977 1.993 1.659 1.667 1.481 1.839 1.549 1.647 1.422 1.320 1.251
1.486
1.467 1.477 1.401 1.238 1.186 1.670 1.533 1.604 1.451 1.448 1.360
Note: “RMSE” and “MAE” stand for the root mean squared error and mean absolute error, respectively. The results for GLS report
the levels of the RMSE and MAE, and those for others are their ratios relative to GLS.

Next, we turn to the simulation results on inference. Tables 4.3 and 4.4 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 heteroskedasticityrobust 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 naturalthat OLS-U performs poorly since it is invalid except for DGP2. The performance of k-NN isworse 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 lengthof 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 aresmaller 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.

15

Table 4.3. Simulation: Inference with univariate covariate

16

n = 50EstimatorEC
AL
GLS (infeasible)
0.956 0.528
0.946 0.627
OLS-U
0.798 1.008
0.492 0.409
OLS-R
0.766 0.962
0.730 0.761
OLS-Boot
0.740 0.885
0.690 0.681
FGLS
0.800 0.451
0.737 0.504
k-NN (Automatic) 0.708 0.410
0.576 0.351
k-NN (k = 6)
0.732 0.410
0.576 0.351
k-NN (k = 15)
0.735 0.418
0.582 0.353
k-NN (k = 24)
0.711 0.440
0.512 0.324
MGLS
0.779 0.499
0.725 0.523
MGLS (Robust)
0.762 0.483
0.725 0.465

DGP 1
n = 100
ECAL
0.955 0.370
0.962 0.441
0.742 0.740
0.421 0.290
0.805 0.862
0.772 0.689
0.845 0.517
0.856 0.527
0.847 0.328
0.812 0.395
0.659 0.258
0.574 0.251
0.666 0.258
0.574 0.251
0.704 0.266
0.592 0.258
0.688 0.269
0.537 0.244
0.812 0.363
0.744 0.392
0.791 0.354
0.744 0.359

n = 500
ECAL
0.956 0.164
0.960 0.196
0.636 0.349
0.348 0.131
0.884 0.648
0.880 0.488
0.907 0.356
0.894 0.345
0.872 0.162
0.885 0.195
0.701 0.102
0.650 0.115
0.621 0.102
0.650 0.115
0.717 0.105
0.677 0.118
0.721 0.107
0.668 0.118
0.905 0.165
0.888 0.188
0.903 0.163
0.888 0.181

n = 50ECAL
0.939 0.749
0.948 0.316
0.939 0.749
0.948 0.316
0.933 0.733
0.874 0.272
0.917 0.718
0.846 0.270
0.916 0.709
0.761 0.231
0.902 0.711
0.927 0.306
0.845 0.711
0.927 0.306
0.929 0.725
0.944 0.310
0.945 0.744
0.953 0.316
0.885 0.640
0.951 0.333
0.879 0.635
0.951 0.258

DGP 2
n = 100
ECAL
0.947 0.516
0.952 0.207
0.947 0.516
0.952 0.207
0.941 0.507
0.881 0.185
0.947 0.516
0.952 0.207
0.925 0.493
0.758 0.159
0.884 0.483
0.917 0.197
0.845 0.483
0.917 0.197
0.907 0.492
0.931 0.200
0.921 0.504
0.942 0.204
0.907 0.468
0.968 0.222
0.902 0.463
0.968 0.177

n = 500
ECAL
0.955 0.224
0.970 0.085
0.955 0.224
0.970 0.085
0.949 0.222
0.935 0.081
0.955 0.224
0.970 0.085
0.935 0.216
0.844 0.072
0.883 0.205
0.881 0.079
0.819 0.205
0.881 0.079
0.914 0.210
0.919 0.081
0.935 0.216
0.939 0.083
0.937 0.219
0.972 0.092
0.933 0.216
0.972 0.079

Note: “EC” and “AL” stand for the empirical coverage probability and average length, respectively. “OLS-U”, “OLS-R”, and“OLS-Boot” use the normal approximation with the usual OLS standard error, the heteroskedasticity robust standard error,
and the percentile bootstrap interval, respectively. “MGLS (Robust)” is based on the variance formula presented in Remark 2.

Table 4.4. Simulation: Inference with multivariate covariates17

n = 50EstimatorEC
AL
GLS (infeasible)
0.944 0.611
0.951 0.632
OLS-U
0.824 1.535
0.589 0.526
OLS-R
0.787 1.411
0.729 0.767
OLS-Boot
0.756 1.319
0.688 0.708
FGLS
0.831 1.069
0.658 0.517
k-NN (Automatic) 0.571 0.481
0.471 0.296
k-NN (k = 6)
0.597 0.481
0.471 0.296
k-NN (k = 15)
0.592 0.504
0.491 0.301
k-NN (k = 24)
0.582 0.557
0.459 0.287
MGLS
0.803 0.956
0.687 0.524
MGLS (Robust)
0.755 0.855
0.687 0.548

DGP 3
n = 100
ECAL
0.951 0.413
0.961 0.439
0.786 1.108
0.549 0.369
0.815 1.232
0.782 0.660
0.781 1.133
0.749 0.593
0.845 0.759
0.722 0.395
0.557 0.289
0.472 0.205
0.574 0.289
0.472 0.205
0.590 0.299
0.484 0.212
0.561 0.313
0.450 0.204
0.863 0.623
0.756 0.401
0.833 0.573
0.756 0.415

n = 500
ECAL
0.946 0.175
0.960 0.194
0.632 0.511
0.491 0.164
0.869 0.873
0.891 0.441
0.839 0.797
0.826 0.400
0.897 0.336
0.826 0.198
0.587 0.105
0.534 0.091
0.549 0.105
0.534 0.091
0.639 0.108
0.570 0.096
0.618 0.111
0.562 0.096
0.938 0.248
0.897 0.198
0.920 0.234
0.897 0.199

n = 50ECAL
0.944 0.636
0.949 0.414
0.819 1.222
0.639 0.420
0.797 1.197
0.762 0.596
0.752 1.115
0.719 0.569
0.801 0.596
0.797 0.382
0.672 0.526
0.596 0.298
0.670 0.526
0.596 0.298
0.690 0.551
0.607 0.309
0.681 0.590
0.598 0.297
0.801 0.805
0.707 0.404
0.762 0.750
0.707 0.422

DGP 4
n = 100
ECAL
0.943 0.440
0.950 0.282
0.780 0.893
0.611 0.298
0.843 1.068
0.810 0.517
0.785 0.951
0.780 0.460
0.823 0.424
0.862 0.262
0.646 0.323
0.599 0.200
0.639 0.323
0.599 0.200
0.675 0.338
0.629 0.210
0.664 0.347
0.608 0.207
0.844 0.526
0.776 0.305
0.833 0.511
0.776 0.308

n = 500
ECAL
0.956 0.192
0.968 0.123
0.675 0.412
0.521 0.133
0.906 0.718
0.914 0.334
0.860 0.654
0.866 0.301
0.834 0.191
0.855 0.112
0.609 0.122
0.605 0.082
0.571 0.122
0.605 0.082
0.689 0.129
0.668 0.088
0.702 0.132
0.662 0.091
0.902 0.216
0.908 0.154
0.902 0.219
0.908 0.148

Note: “EC” and “AL” stand for the empirical coverage probability and average length, respectively. “OLS-U”, “OLS-R”, and“OLS-Boot” use the normal approximation with the usual OLS standard error, the heteroskedasticity robust standard error,
and the percentile bootstrap interval, respectively. “MGLS (Robust)” is based on the variance formula presented in Remark 2.

4.2. 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 theperspective, 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 changein (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 4.5. 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 tothose 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 the1960, 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 estimatesindicate 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 canbe negative without sufficient capital. Hence it would be reasonable to consider Aging as a sourceof heteroskedasticity. The upper panel of Figure 4.2 shows the relationship between the residualfrom 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 panelof Figure 4.2 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 4.5 show the results of MGLS. Panels B and C present theresults 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
18Theorems 1-2 and their analogous versions for IV estimators are reported in parentheses, whilerobust 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, thetwo 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 theconditional 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 resultsof 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 nolarger than those of OLS, which demonstrates the increased precision of MGLS.

19

2

Residual

1

0

−1

−0.2

0.0

0.2

0.4

Change in the ratio of old to young workers from 1990 to 2015

2

Residual

1

0

−1

0.2

0.4

0.6

0.8

The ratio of old to young workers in 1990

Figure 4.1. Plots for residual and aging (upper) and residual and ratio of old
to young workers in 1990 (lower)
Note: For both panels, the residuals are obtained from the regression of the change in GDP per capita from1990 to 2015 (GDP) on the population aging measured by the ratio of the population above 50 to those between
the ages of 20 and 49 (Aging). For the upper panel, the variable on the X-axis represents the change in the ratioof old to young workers from 1990 to 2015. For the lower panel, it represents the ratio of old to young workers in
1990.

20

Table 4.5. Effects of the Aging on GDP by OLS and MGLS
Specification
Panel A: OLS
Aging
Initial GDP
Panel B: MGLS (covariate of σ 2 (·) =Aging)
Aging

Initial GDP

21

Panel C: MGLS (covariate of σ 2 (·) =Initial Ratio)
Aging

Initial GDP
Panel D: MGLS (covariates of σ 2 (·) =All)
Aging

Initial GDP

Sample of all countries (n = 169)
(1)
(2)
(3)
(4)
(5)

OECD countries (n = 35)
(6)
(7)
(8)

0.335
1.036
1.162
0.773 1.703-0.2620.042(0.210) (0.257) (0.276) (0.322) (0.411) (0.352) (0.346)-0.153 -0.138 -0.156 -0.190-0.205(0.039) (0.042) (0.046) (0.045)(0.072)1.186(0.458)-0.260(0.092)

0.387
1.098
1.191
0.751 0.414-0.391 -0.029(0.189) (0.187) (0.205) (0.267) (0.101) (0.247) (0.284)[0.150] [0.179] [0.198] [0.310] [0.472] [0.190] [0.340]-0.164 -0.155 -0.168 -0.079-0.190(0.027) (0.029) (0.030) (0.009)(0.069)[0.031] [0.032] [0.029] [0.046][0.069]-0.458(0.092)[0.451]-0.297(0.025)[0.141]

0.065
0.771
0.894
0.574 0.483-0.501 -0.344(0.196) (0.223) (0.231) (0.235) (0.142) (0.270) (0.219)[0.196] [0.249] [0.262] [0.272] [0.603] [0.231] [0.213]-0.164 -0.141 -0.159 -0.080-0.148(0.031) (0.035) (0.041) (0.012)(0.056)[0.035] [0.037] [0.046] [0.055][0.065]-0.585(0.226)[0.758]-0.379(0.096)[0.288]

0.285
1.064
1.188
0.810 0.494-0.391
0.062(0.221) (0.265) (0.281) (0.289) (0.100) (0.247) (0.274)[0.206] [0.249] [0.271] [0.323] [0.442] [0.190] [0.340]-0.152 -0.136 -0.146 -0.079-0.203(0.030) (0.033) (0.041) (0.009)(0.072)[0.033] [0.036] [0.044] [0.051][0.072]-0.268(0.186)[0.972]-0.250(0.049)[0.197]

Note: For all specifications from (1) to (8), GDP is the dependent variable. Column 1 shows the result of the simple regression of Aging on GDP. Column 2 shows theresult 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. Columns (6), (7) and (8)report the result for OECD countries using specifications (1), (3) and (5), respectively. Panel A reproduces the results by Acemoglu and Restrepo (2017). For Panel A,
heteroskedasticity robust standard errors are presented in parentheses. Panels B, C, and D present the results by MGLS. Panels B and C show the results where theconditional 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 variables except the regional dummies. For Columns (1)-(4) and (6)-(7) of Panels B and C, standard errors based on the formula in Theorem 1 arepresented in parentheses, while those based on the formula in Remark 3 are presented in square brackets. For Columns (1)-(4) and (6)-(7) of Panel D, standard errors
based on the formula in Theorem 2 are presented in parentheses, while those based on the formula analogous to Remark 3 are presented in square brackets. For Columns(5) and (8) of Panels B, C, and D, standard errors are based on the formulae analogous to Remark 3.

Appendix A. Proof of lemma and theorem in Section 2
Notation. In this section, we use the followingqR notation. For a function f (·), we let kf k∞ =|f (x)|2 dP be the L2 (P ) norm; given there issupx∈X |f (x)| be the sup-norm and kf k2,P =no confusion in the context, we use the same set of notations
for a vector a = (a1 , . . . , ak )′ : weqPk2
let kak∞ = maxj∈{1,...,k} |ak | be the sup-norm and kak =j=1 aj be the Euclidean norm.
L [f ](a) be the left derivative of the greatest convex minorant of a function f evaluated atDAa ∈ A, Pn be the empirical measure of {Yi , Xi , Zi }ni=1 , Gn be the empirical process, i.e., Gn f =Pn√1i=1 {f (Xi ) − E[f (Xi )]}, kGn kF = supf ∈F |Gn f |, and IA (x) = I{x ∈ A}. Let τ0 (x) =nσ 2 (x), τ0′ (xL ) be the right derivative of τ0 at xL , τ̂ (x) = σ̂ 2 (x), W be the support of W :=P
(1, X, Z ′ )′ , F (x) be the distribution function of X, Fn (x) = n1 ni=1 I{Xi ≤ x}, and Mn (x) =1 Pn
2i=1 Ûi I{Xi ≤ x}. For a, b ∈ R, let a ∧ b denote min{a, b}, and a . b denote that there exists
na positive constant C such that a ≤ C · b. Let dim(w) be the dimension of a vector w.

A.1. Proof of Lemma 1. Since Ûj = Yj −Wj′ θ̂OLS is the OLS residual, Assumptions A1-A2 andθ̂OLS − θ = Op (n−1/2 ) imply Ûj2 − Uj2 = Op (n−1/2 log n) = op (n−1/3 ) uniformly over j = 1, . . . , n.

To see this, decomposeÛj2 − Uj2=(Yj − Wj′ θ̂OLS )2 − (Yj − Wj′ θ)2=

Wj′ (θ̂OLS + θ) · Wj′ (θ̂OLS − θ) − 2Wj′ θ · Wj′ (θ̂OLS − θ) − 2Uj Wj′ (θ̂OLS − θ)=: Ij + IIj + IIIj .
For Ij , note that
Ij

= [Wj′ (θ̂OLS − θ)]2 + 2Wj′ θ · Wj′ (θ̂OLS − θ)
≤ kWj k2 θ̂OLS − θ≤ R2 θ̂OLS − θ22+ 2 kWj k · kθk · kWj k · θ̂OLS − θ+ 2R2 kθk · θ̂OLS − θ = Op (n−1/2 ),(A.1)

where R is the constant defined in Assumption A1. The first inequality follows from the CauchySchwarz inequality, the second inequality follows from kWj k ≤ R (by Assumption A1), andthe last equality follows from θ̂OLS − θ = Op (n−1/2 ). Note that in the second inequality, theupper bound no longer depends on the index j, so we have maxj |Ij | = Op (n−1/2 ). For IIj ,using the same reasoning as for the first inequality in (A.1), we have maxj |IIj | = Op (n−1/2 ).
Note that here we only consider the second term following the first inequality of (A.1). For

IIIj , the same argument as above yields maxj |Wj′ (θ̂OLS − θ)| = Op (n−1/2 ). Furthermore, by

Assumption A2 and a similar argument after equation (7.11) on p.3297 of Balabdaoui, Durot
and Jankowski (2019) (BDJ hereafter), we have max1≤j≤n |Uj2 | = Op (log n). By the fact thatmax1≤j≤n |Uj | ≤ max1≤j≤n |Uj2 — if max1≤j≤n |Uj | ≥ 1, we havemax |Uj | = Op (log n).

1≤j≤n

(A.2)

In the case of max1≤j≤n |Uj | < 1, max1≤j≤n |Uj | = Op (log n) holds trivially. Combining (A.2)

and maxj |Wj′ (θ̂OLS − θ)| = Op (n−1/2 ), we have maxj |IIIj | = Op (n−1/2 log n). Consequently, we
22

have
max |Ûj2 − Uj2 | ≤ max |Ij | + max |IIj | + max |IIIj |
j

j

j

j

= Op (n−1/2 log n) = op (n−1/3 ).

(A.3)

Furthermore, Assumption A3 guarantees qn∗ − xL = O(n−1/3 ) (by an expansion of qn∗ =F −1 (n−1/3 ) for the quantile function F −1 (·) of X), and we can define c∗ = limn→∞ n1/3 (qn∗ −−1xL ) = dF dq (q)q↓0∈ (0, ∞).

Now, we analyze n1/3 {τ̂ (qn∗ ) − τ0 (xL )}. The term n1/3 {τ̂ (qn ) − τ0 (qn )} will be addressed in

the final step of this subsection. Pick any m > 0. Let

Zn1 (t) = n2/3 [{n−1/3 m + τ0 (xL )}Fn (xL + t(qn∗ − xL )) − Mn (xL + t(qn∗ − xL ))].
Observe that




1/3
∗
−1/3
∗
P n {τ̂ (qn ) − τ0 (xL )} ≤ m = P arg max [{nm + τ0 (xL )}Fn (s) − Mn (s)] ≥ qns∈[xL ,xU ]

= P



arg

max

∗ −x )]
t∈[0,(xU −xL )/(qn
L

n

−2/3


Zn1 (t) ≥ 1 ,(A.4)

where the first equality follows from the switch relation (see a review by Groeneboom andJongbloed, 2014), and the second equality follows from a change of variables s = xL + t(qn∗ − xL )and its implication, s ≥ qn∗ ⇔ t ≥ 1. Let Û (y, w) = y − w′ θ̂OLS and

gn,t (y, w) = n1/6 {τ0 (xL ) − Û (y, w)2 }I[xL ,xL +t(qn∗ −xL )] (x).
We decompose√
Zn1 (t) =n(Pn − P )gn,t + n2/3 E[{τ0 (xL ) − Û (Y, W )2 }I[xL ,xL +t(qn∗ −xL )] (X)]+n1/3 m{Fn (xL + t(qn∗ − xL )) − F (xL + t(qn∗ − xL ))} + n1/3 mF (xL + t(qn∗ − xL ))abcd=: Zn1(t) + Zn1(t) + Zn1(t) + Zn1(t).
a (t). We verify the conditions of van der Vaart (2000, Theorem 19.28). Define
Analysis of Zn1

the class of random functions (depending on θ̂OLS ):Gn1 = {gn,t (y, w) = n1/6 (τ0 (xL ) − Û (y, w)2 )I[xL ,xL +t(qn∗ −xL )] (x) : t ∈ [0, K]},
for K ∈ (0, ∞), where n in the subscript indicates the dependence on both the scaling parametern1/6 and θ̂OLS . By van der Vaart (2000, Example 19.6) we know that for a bracket size ǫ, Gn1

has the entropy with bracketing of order log(1/ǫ). Thus, Gn1 satisfies the entropy condition forvan der Vaart (2000, Theorem 19.28).

23

For each t, s ∈ [0, K], note that
Cov(gn,t , gn,s ) = n1/3 E[{Û (Y, W )2 − τ0 (xL )}2 I[xL ,xL +(t∧s)(qn∗ −xL )] (X)] + op (1)= n1/3 E[{U 2 − τ0 (xL )}2 I[xL ,xL +(t∧s)(qn∗ −xL )] (X)] + op (1)= n1/3 E[[ε2 + {τ0 (X) − τ0 (xL )}2 ]I[xL ,xL +(t∧s)(qn∗ −xL )] (X)] + op (1)Z xL +(t∧s)(qn∗ −xL )1/3= n[σε2 (x) + {τ0 (x) − τ0 (xL )}2 ]fX (x)dx + op (1)=xL
2[σε (ξn ) + {τ0 (ξn ) − τ0 (xL )}2 ]fX (ξn )c∗ (t ∧ s) + op (1)= σε2 (xL )fX (xL )c∗ (t ∧ s) + op (1),(A.5)

for ξn ∈ (xL , xL + (t ∧ s)qn∗ ). The first equality follows from qn∗ − xL = O(n−1/3 ). In the second

equality, we replace the estimated Û 2 with the unobservable U 2 . By (A.3), the discrepancy

between Û 2 and U 2 converges more rapidly than n−1/3 , and the factor I[xL ,xL +(t∧s)(qn∗ −xL )] (X)further refines this rate. Consequently, under Assumptions A1 and A2, the impact of substitutingÛ 2 with U 2 in the second line is op (1). The third equality follows from the definition ε =U 2 − τ0 (X) and E[ε|X] = 0, the fourth equality follows from the law of iterated expectations,the fifth equality follows from a Taylor expansion, and the last equality follows from c∗ =limn→∞ n1/3 (qn∗ − xL ) and the continuity of σε2 (·) and τ0 (·) at xL from right. Similarly, we have
Var(gn,t ) = σε2 (xL )fX (xL )c∗ t + op (1).

We next consider the envelop function of the class Gn1 , that is
Gn1 (y, w) = n1/6 |τ0 (xL ) − Û (y, w)2 | · I[xL ,xL +K(qn∗ −xL )] (x).
We can see that Gn1 is square integrable since similar arguments to (A.5) yield
E[G2n1 (Y, W )] = n1/3 E[|τ0 (xL ) − Û (Y, W )2 | · I[xL ,xL +K(qn∗ −xL )] (X)]= n1/3 E[|τ0 (xL ) − U 2 | · I[xL ,xL +K(qn∗ −xL )] (X)] + op (1)= n1/3 E[[ε2 + {τ0 (X) − τ0 (xL )}2 ] · I[xL ,xL +K(qn∗ −xL )] (X)] + op (1)Z xL +K(qn∗ −xL )1/3[σε2 (x) + {τ0 (x) − τ0 (xL )}2 ]fX (x)dx + op (1)= nxL= Op (1),(A.6)

and thus the Lindeberg condition can be verified by Assumption A2: for any ζ > 0 and someδ > 0,√
E[G2n1 I{Gn1 > ζ n}] ≤n(2+δ)1/6
E[|τ0 (xL ) − Û (Y, W )2 |2+δ · I[xL ,xL +K(qn∗ −xL )] (X)]ζ δ nδ/2=n(2+δ)1/6
E[|τ0 (xL ) − U 2 |2+δ · I[xL ,xL +K(qn∗ −xL )] (X)] + op (1)ζ δ nδ/2= O(n−δ/3 ) + op (1) = op (1),(A.7)

where the inequality follows from the same arguments that are used in the proof of Markov’s
inequality, the first equality follows from θ̂OLS − θ = Op (n−1/2 ) and Assumptions A1-A2, andthe second equality follows from a similar argument to (A.6).
24

Furthermore, as δn → 0, we obtainsup E|gn,t − gn,s |2 = n1/3 sup E[{Û (Y, W )2 − τ0 (xL )}2 I[xL ,xL +|t−s|qn∗ ] (X)]|t−s|≤δn|t−s|≤δn= n

1/3

sup E[[ε2 + {τ0 (X) − τ0 (xL )}2 ] · I[xL ,xL +|t−s|qn∗ ] (X)] + op (δn )|t−s|≤δn= Op (δn ) = op (1).

(A.8)

By (A.5)-(A.8), we can apply van der Vaart (2000, Theorem 19.28), which implies for eachK ∈ (0, ∞),

d

a
Zn1
(t) →

p

σε2 (xL )fX (xL )c∗ Wt in l∞ [0, K].

(A.9)

b (t). Observe that
Analysis of Zn1b
Zn1(t) = n2/3 E[{τ0 (xL ) − U 2 }I[xL ,xL +t(qn∗ −xL )] (X)] + n2/3 E[(U 2 − Û (Y, W )2 )I[xL ,xL +t(qn∗ −xL )] (X)]Z xL +t(qn∗ −xL )2/3{τ0 (xL ) − τ0 (F −1 (F (x)))}dF (x) + op (1)= nxL= n2/3Z F (xL +t(qn∗ −xL ))F (xL )

= −n

2/3

= −n

2/3

{τ0 (xL ) − τ0 (F −1 (v))}dv + op (1)Z F (xL +t(qn∗ −xL ))F (xL )Z F (xL +t(qn∗ −xL ))F (xL )

= −n2/3 τ0′ (xL )= −τ0′ (xL )τ0′ (xL ){F −1 (v) − F −1 (F (xL ))}dv + op (1)τ0′ (xL )v − F (xL )dv + op (1)fX (xL ){F (xL + t(qn∗ − xL )) − F (xL )}2+ op (1)2fX (xL )t2 (c∗ )2
fX (xL ) + op (1)2(A.10)

holds uniformly over t ∈ [0, K], where the second equality follows from E[{U 2 − Û (Y, W )2 } ·I[xL ,xL +t(qn∗ −xL )] (X)] = op (n−2/3 ), the third equality follows from a change of variables v = F (x),the fourth equality follows from a Taylor expansion, the fifth equality follows from F −1 (v)−xL =1fX (xL ) (v − F (xL )) + o(v − F (xL )), the sixth equality follows from evaluating the integral, andthe last equality follows from a Taylor expansion and c∗ = limn→∞ n1/3 (qn∗ − xL ).
c (t). By Kim and Pollard (1990, Maximal inequality 3.1),
Analysis of Zn1

E

"

#

sup |Fn (xL + t(qn∗ − xL )) − F (xL + t(qn∗ − xL ))|t∈[0,K]≤ n−1/2 J

p

P G2n

holds for some constant J ∈ (0, ∞). Here Gn is the envelope of the set of indicator functions,thus P G2n ≤ 1. As a result,cZn1(t) ≤ n1/3 n−1/2 mJ

uniformly over t ∈ [0, K].

p

P G2n = o(1),(A.11)

d (t). A Taylor expansion yields
Analysis of Zn1dZn1(t) = n1/3 mF (xL + t(qn∗ − xL )) = m · t · fX (xL )c∗ + o(1),25(A.12)

uniformly over t ∈ [0, K], for every K < ∞.

Combining (A.9)-(A.12), it holds that for each 0 < K < ∞,dZn1 (t) → Z1 (t) :=p

σε2 (xL )fX (xL )c∗ Wt − τ0′ (xL )t2 (c∗ )2
fX (xL ) + m · t · fX (xL )c∗ in l∞ [0, K].
2(A.13)

We now verify the conditions of the argmax continuous mapping theorem (Kim and Pollard,1990). Note that for each t 6= s,
Var(Z1 (s) − Z1 (t)) = σε2 (xL )fX (xL )c∗ |t − s| =6 0.
By Kim and Pollard (1990), the process t → Z1 (t) achieves its maximum a.s. at a unique point.

Consider extended versions of Zn1 and Z1 to the real line:Zn1 (t), t ≥ 0
Z1 (t),Z̃n1 (t) =,Z̃1 (t) =
t
t
t<0

t≥0

.

t<0

d

It holds Z̃n1 (t) → Z̃1 (t), and the similar argument to Lemma SM.2.1 (ii) in Babii and Kumar(2023) yields that the maximum of Z̃n1 (t) is uniformly tight. Therefore, by Kim and Pollard(1990, Theorem 2.7),


P n
= P

= P

1/3

{τ̂ (qn∗ ) − τ0 (xL )} ≤ m
s



→P





arg maxZ1 (t) ≥ 1
t≥0

#
!
2 c∗
t
σε2 (xL )Wt − τ0′ (xL )+ mt ≥ 1
arg maxt≥0c∗ fX (xL )
2
s
! #
!
"
2 c∗
2 (x )tσLεL(1) ≤ m ,Wt + τ0′ (xL )D[0,∞)
c∗ fX (xL )2"where the second equality follows from the switch relation and symmetry of the process Wt .
Thus, we have

n

1/3

d
L
{τ̂ (qn∗ ) − τ0 (xL )} → D[0,∞)

s

t2 c∗
σε2 (xL )
′
W
+
τ
(x
)
t
0 L
c∗ fX (xL )
2

!

(1),

(A.14)

which also implies
n1/3 {τ̂ (qn∗ ) − τ0 (qn∗ )}s!
2 c∗
2 (x )tσdLεL(1) − lim n1/3 {τ0 (qn∗ ) − τ0 (xL )}→ D[0,∞)
Wt + τ0′ (xL )n→∞c∗ fX (xL )
2
s
!
2 c∗
2 (x )σtdLεL∼ D[0,∞)
Wt + τ0′ (xL )− τ0′ (xL )c∗ t (1),c∗ fX (xL )2(A.15)

L
where the distribution relation follows from the fact that the D[0,∞)is a linear operator for alinear function of t.
Finally, we analyze n1/3 {τ̂ (qn )−τ0 (qn )}. Recall qn is the (n−1/3 )-th sample quantile of X. As-

sumption A3 guarantees qn − qn∗ = Op (n−1/2 ) = op (n−1/3 ), which also implies plimn→∞ n1/3 (qn −xL ) = limn→∞ n1/3 (qn∗ − xL ) = c∗ . Thus, the same argument for (A.14) can be applied to show26

that the result in (A.14) holds true even if we replace qn∗ with qn . Therefore, the conclusion
follows.

A.2. Proof of Theorem 1. By the definitions of the estimators, it holds that
−1 


X
X
√
1
1
σ̂i−2 Wi Wi′   √σ̂i−2 Wi Ui  ,n(θ̂ − θ) = nni:xi >qn√i:xi >qnnn(θ̂IGLS − θ) =1 X −2σi Wi Wi′
n
i=1

!−1

n

1 X −2√σi Wi Ui
n
i=1

!

.

Thus it is sufficient for the conclusion to shown
1 X −2
1 X −2p√√
σ̂i Wi Ui −σi Wi Ui → 0,
n
n

T1 :=

i:xi >qni=1

1
1 X −2σ̂i Wi Wi′ −
n
n

T2 :=

i:xi >qn

n
X
i=1

p

σi−2 Wi Wi′ → 0.

A.2.1. The concentration of T1 . Decompose
1 X
1 X −2T1 = √σi Wi Ui =: T11 − T12 .
(σ̂i−2 − σi−2 )Wi Ui − √nni:xi >qni:xi ≤qn

h be the h-th element of W
We first consider T12 . For any h ∈ {1 : dim(W )}, let Wih and T12ip

h |q ] = 0 by E[U |W ] = 0. Also we have Var(T h |q ) → 0.
and T12 , respectively. Note that E[T12n
12 n

To see this, decompose
where Ih = n1 E

Ih , note that
Ih =

h P

h
Var(T12|qn ) = Ih − n · (IIh )2 ,in−2
hU 2 qand IIh = E[I{Xi ≤ qn }σi−2 Wih Ui |qn ]. For
I{X
≤
q
}σ
W
n
i
n
i
i
i=1
i





1  
E E
n

n
X
i=1

I{Xi ≤ qn }σi−2 Wih Ui

!2





W  qn = E[E[(I{Xi ≤ qn }σi−2 Wih Ui )2 |W]|qn ] = E[I{X ≤ qn }σ −2 (X)(W h )2 |qn ]p≤ R2 σ −2 (xL )E[I{X ≤ qn }|qn ] → 0,

where W = (W1 , . . . , Wn )′ . The first equality follows from the law of iterated expectation andthe fact that qn is a function of W, the second equality follows from E[U |W ] = 0 and {Ui }ni=1being iid, the third equality follows because conditional on W, I{Xi ≤ qn }(σi−2 Wih )2 is treatedas fixed, the inequality follows from Assumptions A1 and A2, and the convergence follows frompqn → xL . For IIh , note thatIIh = E[I{Xi ≤ qn }σi−2 Wih E[Ui |W]|qn ] = 0,
where the first equality follows from the law of iterated expectation and the fact that qn is ah |q ] =function of W, and the second equality follows from E[Ui |W] = E[Ui |Wi ] = 0. Since E[T12npph |q ) → 0 hold for every h, we can conclude that T
0 and Var(T12n
12 → 0.
27

To proceed, we will utilize Lemma 3 below. Its proof can be found at the end of Appendix
A.2. Recall that earlier in this appendix, we relabel σ 2 (·) as τ0 (·), and τ̂ is used to denote theisotonic estimator of σ 2 (·). Additionally, with some abuse of notation, we use wh to denote theh-th element of vector w.

Lemma 3. Under Assumptions A1-A3,(i): kτ̂ k∞ = Op (log n),(ii): kτ̂ − τ0 k22,P = Op ((log n)2 n−2/3 ),(iii): E[kGn kFn ] ≤ Aν
2 holds for any constants A > 0 and ν > 0, and all sufficiently largen, where Fn is the function class defined as





1
1
−
Fn = fn (w, u) = I{x > qn }wh u :
τ (x) τ0 (x)τ ≥ 0 is monotone increasing on X ,kτ k∞ ≤ C log n, kτ − τ0 k22,P ≤ Crn ,I{x > qn }/τ (x) ≤ 1/K0 , h ∈ {1 : dim(w)}(A.16)

with C and K0 being some positive constants, and rn = (log n)2 n−2/3 .







Now we focus on T11 . Since the proof is similar, we only present the proof for the h-th elementof T11 , i.e., for any constant A > 0,
P {|Gn fˆ| ≥ A} → 0,
(A.17)

1where fˆ(w, u) = I{x > qn } τ̂ (x)− τ01(x) wh u. To this end, we set τ0 (xL ) = C0 = 2K0 > 0. It

holds that for any A > 0 and ν > 0, there exists a positive constant C such that1
νI{x > qn }
2
ˆ
ˆ
≤
+
P {|Gn f | ≥ A} ≤ P |Gn f | ≥ A, kτ̂ k∞ ≤ C log n, kτ̂ − τ0 k2,P ≤ Crn ,τ̂ (x)
K02

E kGn kFnν≤+ ≤ ν,(A.18)
A2

for all sufficiently large n, where the first inequality follows from Lemma 1 and Lemma 3 (i)(ii), and the fact that τ̂ is monotone increasing (so that the lower bound at the truncation
point is the uniform lower bound). Specifically, for any ν > 0, we can find C > 0 and apositive integer n0 such that for any integer n >that (a) P {kτ̂ k∞ > C log n} <n n0 , it holds o2ν
6 , (b) P {kτ̂ − τ0 k2,P

> Crn } < ν6 , and (c) P

I{x>qn }> K10τ̂ (x)< ν6 . Parts (a) and (b) areensured by Lemma 3 (i) and (ii), respectively; part (c)n is guaranteedo by Lemma 1. As aresult, P {kτ̂ k∞ > C log n} or {kτ̂ − τ0 k22,P > Crn } orI{x>qn }> K10τ̂ (x)< ν2 . In the case of2limx↓xL dσdx(x) = 0, part (c) remains valid since τ̂ (x) will converge to τ0 at a faster rate (the√n-rate), then the first inequality of (A.18) holds without invoking Lemma 1.

The second inequality of (A.18) follows from Markov’s inequality and the definition of Fn ,which is given by (A.16). The last inequality follows from Lemma 3 (iii). Since ν can bearbitrarily small, we obtain (A.17) and the conclusion follows.
28,

A.2.2. Proof of T2 = op (1). Note thatT2=n
1 X −2
1 X −2′σi Wi Wi′σ̂i Wi Wi −nni:xi >qn

=

i=1

1 X
1 X −2σi Wi Wi′(σ̂i−2 − σi−2 )Wi Wi′ −nni:xi >qni:xi ≤qn=: T21 − T22 .
pp

First, we have T22 → 0 since qn → xL . For T21 , let sn be the (1 − n−1/3 )-th sample quan-tile of {Xi }ni=1 . By employing arguments similar to those in the proof of Lemma 1, we haveσ̂ 2 (sn ) − σ 2 (sn ) = Op (n−1/3 ). Using reasoning akin to, yet simpler than, those in the proof of
Lemma 1, we can establish that for any x ∈ (qn , sn ), it holds that σ̂ 2 (x) − σ 2 (x) = Op (n−1/3 ).

Combining the aforementioned results with the monotonicity of both σ̂ 2 (·) and σ 2 (·), we canconclude that supx∈[qn,sn ] |σ̂ 2 (x) − σ 2 (x)| = Op (n−1/3 ), i.e., σ̂ 2 (x) is uniformly consistent within

trimmed domain [qn , sn ] (the proof here resembles the one given for the Glivenko-Cantelli Theorem regarding the uniform consistency of the empirical distribution function; see, for example,
the proof of Theorem 19.1 in van der Vaart, 2000). Therefore, we have
X
1
1 X
T21 =
(σ̂i−2 − σi−2 )Wi Wi′ +(σ̂i−2 − σi−2 )Wi Wi′ = op (1),nni:qn <xi <sn(A.19)

i:xi ≥snwhere the second equality follows from the preceding argument, |sn − xU | = Op (n−1/3 ), andppLemma 3 (i). Combining T22 → 0 and (A.19), we have T2 → 0.

A.2.3. Proof of Lemma 3 (i). The min-max formula of the isotonic regression says
Pn
Pk
2
2
j=k Ûj
j=1 Ûj≤ τ̂ (x) ≤ max,
min1≤k≤n n − k + 1
1≤k≤nkfor each x ∈ X , which implies min1≤j≤n Ûj2 ≤ τ̂ (x) ≤ max1≤j≤n Ûj2 for each x ∈ X . Thus, it issufficient for the conclusion to show thatmax Ûj2 = Op (log n).

1≤j≤n(A.20)

Observe thatmax Ûj2 ≤ max Uj2 + 2Rk||θ̂OLS − θ||∞ max |Uj | + R2 k2 ||θ̂OLS − θ||2∞ .

1≤j≤n

1≤j≤n

1≤j≤n

From BDJ (2019, eq. (7.11) on p.3297), Assumption A2 guarantees max1≤j≤n Uj2 = Op (log n).
Since θ̂OLS is the OLS estimator, it holds that ||θ̂OLS − θ||∞ = Op (n−1/2 ). By (A.2), we also
have max1≤j≤n |Uj | = Op (log n). Combining these results with Assumption A1, we have (A.20).

A.2.4. Proof of Lemma 3 (ii). The proof is based on that of Proposition 4 of BGH (p.8 ofP
BGH-supp). Recall that τ̂ (·) is the solution of minτ ∈{all monotone funcitons} nj=1 {Ûj2 − τ (Xj )}2 ,or equivalentlymaxτ ∈{all monotone funcitons}n
X{2Ûj2 τ (Xj ) − τ (Xj )}.
j=1
29(A.21)

On the other hand, τ0 (·) is the solution of minτ ∈{all monotone funcitons} E[{U 2 − τ (X)}2 ], or equivalentlymaxτ ∈{all monotone funcitons}

E[2U 2 τ (X) − τ (X)2 ].

(A.22)

By (A.21), it holdsnnX
X{2Ûj2 τ0 (Xj ) − τ0 (Xj )2 },{2Ûj2 τ̂ (Xj ) − τ̂ (Xj )2 } ≥
j=1

j=1

or equivalently (by plugging in Ûj = Uj − Wj′ (θ̂OLS − θ)),
n
X
j=1

≥

n
X
j=1

{2Uj2 τ̂ (Xj ) − τ̂ (Xj )2 } + 2

n 
X
j=1


−2Uj Wj′ (θ̂OLS − θ) + {Wj′ (θ̂OLS − θ)}2 {τ̂ (Xj ) − τ0 (Xj )}{2Uj2 τ0 (Xj ) − τ0 (Xj )2 }.

(A.23)

Define d22 (τ1 , τ2 ) = −E[2τ1 τ2 − τ12 − τ22 ]. Note that for any monotone function τ ,E[2U 2 τ (X) − τ (X)2 ] − E[2U 2 τ0 (X) − τ0 (X)2 ]= E[2E[U 2 |X]τ (X) − τ (X)2 − 2E[U 2 |X]τ0 (X) + τ0 (X)2 ]= E[2τ0 (X)τ (X) − τ (X)2 − τ0 (X)2 ] = −d22 (τ, τ0 ),(A.24)

where the first equality follows from the law of iterated expectation, the second equality follows
from the definition τ0 (x) = E[U 2 |X = x], and the last equality follows from the definition of
d22 (·, ·).

Definegτ (u, x) = {2u2 τ (x) − τ (x)2 } − {2u2 τ0 (x) − τ0 (x)2 },
n

2 X
Rn =
−2Uj Wj′ (θ̂OLS − θ) + {Wj′ (θ̂OLS − θ)}2 {τ̂ (Xj ) − τ0 (Xj )}.
n
j=1

From (A.23) and (A.24), it holds
Zgτ̂ (u, x)d(Pn − P )(u, x) + Rn ≥ d22 (τ̂ , τ0 ).

(A.25)

Note that Rn is bounded as|Rn | ≤ −(θ̂OLS − θ)′

n

n

j=1

j=1

2X
4X
Wj Uj {τ̂ (Xj ) − τ0 (Xj )} +{Wj′ (θ̂OLS − θ)}2 {τ̂ (Xj ) − τ0 (Xj )} .
nnThe second term is of order Op (n−1 log n) (because θ̂OLS − θ = Op (n−1/2 ) and Lemma 3 (i)).

By similar arguments in p.22 of BGH-supp and in the proof of Lemma 3 (i), the first term is oforder Op n−1 (log n)2 .
ThenRn = Op n−1 (log n)2 .

30(A.26)

Thus, for some constants C, K > 0 and a shrinking sequence ǫn , set inclusion relationships yield


Z
2
2
2
P (d2 (τ̂ , τ0 ) ≥ ǫn ) = P d2 (τ̂ , τ0 ) ≥ ǫn ,gτ̂ (u, x)d(Pn − P )(u, x) + Rn ≥ d2 (τ̂ , τ0 )≤d2 (τ̂ , τ0 ) ≥ ǫn , |Rn | ≤ Cn−1 (log n)2 , kτ̂ k∞ ≤ K log nRgτ̂ (u, x)d(Pn − P )(u, x) + Rn − d22 (τ̂ , τ0 ) ≥ 0

P!

+P (|Rn | > Cn−1 (log n)2 ) + P (kτ̂ k∞ > K log n)=: P1 + P2 + P3 ,where the first equality follows from (A.25). For P2 and P3 , (A.26) and Lemma 3 (i) imply that
we can choose C and K to make these terms arbitrarily small. Thus, we focus on the first term
P1 .
Now letT= {τ : τ is positive and monotone increasing on X , kτ k∞ ≤ K log n},

G = {gτ (u, x) = {2u2 τ (x) − τ (x)2 } − {2u2 τ0 (x) − τ0 (x)2 } : τ ∈ T },Gv = {g ∈ G : d2 (τ, τ0 ) ≤ v}.

Set inclusion relationships and Markov’s inequality yield
!
Z
2
−1
2
{ gτ (u, x)d(Pn − P )(u, x) − d2 (τ, τ0 )} ≥ −Cn (log n)
P1 ≤ Psupτ ∈T ,d2 (τ,τ0 )≥ǫn

≤
≤
≤

∞
X

s=0
∞
X

s=0
∞
X
s=0

P

√

sup
τ ∈T ,2s ǫn ≤d2 (τ,τ0 )≤2s+1 ǫn
P kGn gkG s+1
2

ǫn

≥

√

!
Z

√
n{ gτ (u, x)d(Pn − P )(u, x)} ≥ n 22s ǫ2n − Cn−1 (log n)2n 22s ǫ2n − Cn−1 (log n)2
√
E[kGn gkG s+1 ]/{ n 22s ǫ2n − Cn−1 (log n)2 }.
2

ǫn

2

For a sufficiently large constant C̃ > 0, the sequence ǫ2n := C̃(log n)2 n− 3 dominates Cn−1 (log n)2 , √√so it holds n 22s ǫ2n − Cn−1 (log n)2 = n22s ǫ2n (1 + o(1)). Therefore, the standard result for

the L2 -convergence of the isotonic estimator under Assumption A2 (e.g., pp. 8-11 in BGH-supp)

implies that the last term can be made arbitrarily small by appropriately selecting C̃. Thus,
the proof is concluded.
A.2.5. Proof of Lemma 3 (iii). We show E[kGn kFn ] ≤ Aν
2 by using van der Vaart and Wellner(1996, Lemma 3.4.3). First we introduce some notation for this part. Let N[] (ε, F, ||·||) be the ε-bracketing number of the function class F under the norm ||·||, HB (ε, F, ||·||) = logN[] (ε, F, ||·||)
Rδp
be the entropy, Jn (δ, F, || · ||) = 0 1 + HB (ε, F, || · ||)dε, and kf kB,P = (2E[e|f | − |f | − 1])1/2
be the Bernstein norm.

Lemma 3.4.3 in van der Vaart and Wellner (1996): Let F be a class of measurable

functions such that kf k2B,P ≤ δ for every f in F. Then√
E[kGn kF ] . Jn (δ, F, || · ||B,P ){1 + Jn (δ, F, || · ||B,P )/( nδ2 )}.
31

To apply this lemma, we need to compute HB (ǫ, F˜n , k · kB,P ) and kf˜k2B,P , where F˜n = {f˜ =D −1 f : f ∈ Fn }, the function class Fn is defined below in (A.27), and the constant D > 0 will
be chosen later to guarantee that the Bernstein norm of f˜ is finite. Moreover, let us define thefollowing function class:
TI,K1 = {τ monotone non-decreasing on the interval I and 0 < τ < K1 }.
Assumption A2 implies that there exist positive constants, C and C, such that 0 < C < τ0 <C < ∞. Also let(

Fn =fn (w, u) = I{x > qn }



1
1
−
τ (x) τ0 (x)wh u :

τ ∈ TX ,K1 , kτ − τ0 k22,P ≤ v 2 ,I{x > qn }/τ (x) ≤ 1/K0 , h ∈ {1 : dim(w)}(A.27))

where wh is the h-th component of vector w. We set 2K0 = C, K1 = K2 log n, and v =K3 (log n)n−1/3 for some constants K2 , K3 > 0.
Consider ǫ-brackets (τ L , τ U ) under the L2 (P )-norm for the functions in TI,K1 . According tovan der Vaart and Wellner (1996, Theorem 2.7.5), there exists some constant C ′ > 0 such that
HB (ǫ, TX ,K1 , k · k2,P ) ≤C ′ K1,ǫfor each ǫ ∈ (0, K1 ).

(A.28)

Without loss of generality, we can choose those bracket functions that satisfy I{x > qn }/τ L (x) ≤1/K0 .5 DefineI{x > qn } U1 − 1wh u τ (x) τ0 (x) f L (w, u) =1
I{x > q }− τ01(x) wh un
τ L (x)
I{x > qn } L1 − 1wh u τ (x) τ0 (x) f U (w, u) =I{x > q }1− τ01(x) wh un
τ U (x)

if wh u ≥ 0,
if wh u < 0,
if wh u ≥ 0,
if wh u < 0.

Note that (f L , f U ) is a bracket of f ∈ Fn for every qn ∈ [xL , xU ].
Now we compute the bracket size of (f˜L , f˜U ) := (D −1 f L , D −1 f U ) with respect to the Bernstein norm. Note thatkf˜U − f˜L k2B,P = kD −1 f U − D −1 f L k2B,PZ∞
kXτ U (x) − τ L (x)
1
w
u
≤ 2
dP (w, u)hk!D k W ×R τ L (x)τ U (x)
k=2
)
(

2 X

∞
∞ 
X
1
R
2RM0 K1 k 2
Rk k!M0k−2 a0 (2K1 )k−2 UL 2kτ − τ k2,P ≤ 2a0
≤ 2
ǫ ,
k!D k
DK02
DK02
K02k
k=2
k=0
where the first inequality follows from the definition of k · k2B,P and I{x > qn } ≤ 1, the secondn}≤ K10 .
inequality follows from Assumption A2 (where we can choose a0 , M0 > 1) and I{x>qτ L (x)5By definition (A.27), the τ (·) associated to F

n must satisfy I{x > qn }/τ (x) ≤ 1/K0 . Since TX ,K1 is a classof monotone increasing function, any ǫ-brackets of TX ,K1 can be modified to be a ǫ-bracket of the “Fn -subset”
of TX ,K1 , satisfying I{x > qn }/τ (x) ≤ 1/K0 by leveling-up certain part of lower bounds functions τ L , withoutchanging the bracket numbers, and the size of each modified bracket can only be smaller.

32,Thus, by setting D = 4M0 RK1 /K02 , we obtainkf˜U − f˜L k2B,P ≤ 4Ma20K 2 ǫ2 , which implies
01

kf˜U − f˜L kB,P ≤ K̃ǫ,
1/2

afor K̃ = 2M00 K1 . Note that (f˜L , f˜U ) is: (a) a set of brackets in F̃n , (b) one-to-one induced by(τ L , τ U ), an ǫ-bracket in TX ,K with the entropy HB (ǫ, TX ,K , k · k2,P ), and (c) kf˜U − f˜L kB,P ≤1

1K̃ǫ. Based on these facts, (A.28) yields

HB (K̃ǫ, F̃n , k · kB,P ) ≤ HB (ǫ, TX ,K1 , k · k2,P ) ≤C ′ K1,ǫwhich implies (by a change-of-variable argument)
HB (ǫ, F̃n , k · kB,P ) ≤K̃C ′ K1
B̃
= ,
ǫ
ǫ

1/2

for B̃ =C ′ a0.
2M0(A.29)

We now characterize the Bernstein norm of f˜,Z∞kXτ (x) − τ0 (x)1dP (w, u)wu
kf˜k2B,P ≤ 2
hk!D k W ×R τ (x)τ0 (x)
k=2
)
(
∞
X
Rk k!M0k−2 a0 (2K1 )k−2
1
kτ − τ0 k22,P
≤ 2
2k
k!D k
K
0
k=2

2 X

∞ 
R
a0 1 2
2RM0 K1 k 2
≤ 2a0
v ≤
v ,
DK02
DK02
4M02 K12
k=0n}≤ K10 , and the third inequality follows from
where the second inequality follows from I{x>qτ (x)(A.27) and some rearrangements. Then, we havekf˜kB,P ≤

1/2

Bv
,
K1

for B =a0.
2M0

(A.30)

Combining (A.29) and (A.30), van der Vaart and Wellner (1996, Lemma 3.4.3) implies
Jn (BK −1 v)
E[kGn kF̃n ] . Jn (BK1−1 v) 1 + √ 2 21 2 ,nB v /K1
where Jn (·) is the abbreviation of Jn (·, F̃n , k · kB,P ). By the arguments used in the proof of

Proposition 7.9 of BDJ, it holds−1/2 1/2

Jn (BK1−1 v) ≤ BK1−1 v + 2B̃ 1/2 B 1/2 K1v−1/2 1/2. B1 K1v,

for some B1 > 0 and sufficiently small v. This implies−1/2 1/2v
E[kGn kF˜n ] . B1 K1

−1/2

B1 K
v 1/2
1 + K12 √ 1 2 2
nB v

!

−1/2
. B1 K1 v 1/2

3/2

B2 K1
1 + √ 3/2
nv!

,

for some B2 > 0. By the definition of the class F̃n = {f˜ = D −1 f : f ∈ Fn }, it follows that
!
!
3/2
3/2
B2 K1
B2 K1
−1/2 1/2−2 1/2 1/2
1 + √ 3/2 . B3 K0 K1 v
1 + √ 3/2 ,v
E[kGn kFn ] = D·E[kGn kF̃n ] . DB1 K1
nv
nv

33

for some B3 > 0. The conclusion follows by observing that with v = K3 (log n)n−1/3 , K1 =K2 log n, and all sufficiently large n, we haveE[kGn kFn ] . C3 (log n)n−1/6 (1 + C4 ) .
1/2

1/2

where C3 = B3 K0−2 K2 K3Aν,
2and C4 = B2 (K2 /K3 )3/2 .

Appendix B. Proof of lemma and theorem in Section 3
Notation. To avoid heavy notations, some of them are used in Appendix A but redefined
here. Define τη (a) = E[σ 2 (X ′ η0 )|X ′ η = a] and τη0 (a) = τ0 (a) (note that τ0 (x′ η0 ) = σ 2 (x′ η0 )).
Let τ̂η = τ̂η (x′ η) be the isotonic estimator obtained by (3.3) for a given η, W be the support ofP
P
W := (1, X ′ , Z ′ )′ , Fn (t) = n1 ni=1 I{Xi′ η̂ ≤ t}, and Mn (t) = n1 ni=1 Ûi2 I{Xi′ η̂ ≤ t}.
B.1. Proof of Lemma 2. The main part of the proof is similar to that of Lemma 1. Recall

that qn∗ is the (n−1/3 )-th population quantile of (X ′ η0 ) and qn is the (n−1/3 )-th sample quantileof {Xi′ η̂}ni=1 with η̂ estimated by (3.4). To proceed, we use the following lemma:
Lemma 4. Under Assumptions M1-M6, it holds(i): η̂ − η0 = Op (n−1/2 ),(ii): τη̂ (a) − τ0 (a) = Op (n−1/2 ) for each a, and kτη̂ − τ0 k2,P = Op (n−1/2 ).
The proof of this lemma is in Appendix B.3. Based on Lemma 4 (i), Assumptions M2-M3,
and properties of the sample quantile, we obtain qn − qn∗ = Op (n−1/2 ) = op (n−1/3 ), which impliesc∗ = limn→∞ n1/3 (qn∗ − xL ) = plimn→∞ n1/3 (qn − xL ) < ∞. By Assumption M2, Lemma 4 (ii),and similar arguments in Appendix A.1, we haven1/3 {τ̂η̂ (qn ) − τ0 (qn )} = n1/3 {τ̂η̂ (qn ) − τη̂ (qn )} + op (1)n1/3 [{τ̂η̂ (qn ) − τη̂ (xL )} − {τη̂ (qn ) − τ0 (xL )}] + op (1)s!
2 c∗
2 (x )tσdε LL→ D[0,∞)(1) − plim n1/3 {τ0 (qn ) − τ0 (xL )}Wt + τ0′ (xL )c∗ fX (xL )
2
n→∞
s
!
t2 c∗
σε2 (xL )d′L(1) − lim n1/3 {τ0 (qn∗ ) − τ0 (xL )}
W
+
τ
(x
)
∼ D[0,∞)tL0n→∞c∗ fX (xL )
2
s
!
2 c∗
2 (x )σtdLεL∼ D[0,∞)
Wt + τ0′ (xL )− τ0′ (xL )c∗ t (1),c∗ fX (xL )2=where the first and second equalities follow from Lemma 4 (ii), the convergence follows from asimilar argument to (A.15), the first distribution relation follows from Lemma 4 (ii), Assumption
M2(iv), and qn∗ − qn = op (n−1/3 ), and the second distribution relation follows from the fact thatL
the D[0,∞)is a linear operator for a linear function of t.

B.2. Proof of Theorem 2. Similar to Theorem 1, it is sufficient for the conclusion to prove
the following lemma.
Lemma 5. Under Assumptions M1-M6, it holds
34

(i): kτ̂η k∞ = Op (log n) uniformly over η ∈ B(η0 , δ0 ),(ii): kτ̂η̂ − τ0 k22,P = Op ((log n)2 n−2/3 ),(iii): E kGn kFn ≤ Aν
2 holds for any constants A > 0 and ν > 0, and all sufficiently largen, where Fn is the function class defined asτ ≥ 0 is monotone increasing on Iη ,






kτ k∞ ≤ C log n, kτ − τη k22,P ≤ Crn ,1
1′−
Fn = fn (w, u) = I{x η > qn }wh u :
τ (x′ η) τη (x′ η)I(x′ η > qn )/τ (x′ η) ≤ 1/K0 ,h ∈ {1 : dim(w)}with C and K0 being some positive constants, and rn = (log n)2 n−2/3 .













B.2.1. Proof of Lemma 5 (i). The proof is adapted from BDJ (2019, eq. (7.11) on p.3297).
2 }n2 n
For fixed η, let {Ûη,i
i=1 be a permutation of {Ûj }j=1 , which is arranged according to themonotonically ordered series {Xi′ η}ni=1 . The min-max formula of the isotonic regression says
Pn
Pk
2
2
i=k Ûη,i
i=1 Ûη,i′≤ τ̂η (x η) ≤ max,
min1≤k≤n n − k + 1
1≤k≤nkfor each x ∈ X and η ∈ B(η0 , δ0 ), which implies min1≤j≤n Ûj2 ≤ τ̂η (x′ η) ≤ max1≤j≤n Ûj2 foreach x ∈ X . Thus, it is sufficient for the conclusion to show thatmax Ûj2 = Op (log n).

1≤j≤n

(B.1)

Observe thatmax Ûj2 ≤ max Uj2 + 2Rk||θ̂OLS − θ||∞ max |Uj | + R2 k2 ||θ̂OLS − θ||2∞ ,

1≤j≤n

1≤j≤n

1≤j≤n

where k is the dimension of θ. From Lemma 7.1 of BDJ, Assumption M2 guarantees max1≤j≤n Uj2 =Op (log n). By the same reasoning for the proof of Lemma 3, we have max1≤j≤n |Uj | = Op (log n)

and ||θ̂OLS − θ||∞ = Op (n−1/2 ). Thus, we have kτ̂η k∞ = Op (log n). Since different η only2 }n2changes the permutation {Ûη,i
i=1 but not max1≤j≤n Ûj , we have kτ̂η k∞ = Op (log n) uniformly

over η ∈ B(η0 , δ0 ).

B.2.2. Proof of Lemma 5 (ii). The main part of the proof is similar to those of Lemma 3 (ii)and Proposition 4 of BGH-supp. Definegη,τ (u, x) = {2u2 τ (x′ η) − τ (x′ η)2 } − {2u2 τη (x′ η) − τη (x′ η)2 },n
2 X−2Uj Wj (θ̂OLS − θ) + {Wj (θ̂OLS − θ)}2 {τ̂η (Xj′ η) − τη (Xj′ η)},Rn,η =nj=1

d22 (τ1 , τ2 )= −E[2τ1 τ2 − τ12 − τ22 ],35,

Following reasoning similar to that presented for (A.21)-(A.26), we have for some C and K,P(

supη∈B(η0 ,δ0 )≤d22 (τ̂η , τη ) ≥ ǫ2n )
supη∈B(η0 ,δ0 ) d2 (τ̂η , τη ) ≥ ǫn , supη∈B(η0 ,δ0 ) kτ̂η k∞ ≤ K log n,RP  supη∈B(η0 ,δ0 ) gη,τ̂ (u, x)d(Pn − P )(u, x) + Rn,η − d22 (τ̂η , τη ) ≥ 0, supη∈B(η0 ,δ0 ) |Rn,η | ≤ Cn−1 (log n)2+P (|Rn,η | > Cn−1 (log n)2 ) + Psupη∈B(η0 ,δ0 )kτ̂η k∞ > K log n!

=: P1 + P2 + P3 .
Lemma 5 (i) implies P3 → 0, and P2 → 0 follows from similar arguments for (A.26). For P1 ,we defineT= {τ : τ is positive and monotone increasing function on Iη , kτ k∞ ≤ K log n},

G = {g(x, u) = {2u2 τ (x′ η) − τ (x′ η)2 } − {2u2 τη (x′ η) − τη (x′ η)2 } : τ ∈ T },Gv = {g ∈ G : d2 (τ, τη ) ≤ v},

for each η ∈ B(η0 , δ0 ). By similar arguments for Lemma 3 (ii) and Proposition 4 of BGH-supp,
we can obtain

P1 ≤

∞
X

i √
h
/{ n22s ǫ2n − Cn−1/2 (log n)2 },E kGn gkG s+1

supZ{τ̂η (x′ η) − τη (x′ η)}2 dF (x) = Op ((log n)2 n−2/3 ).

s=0

and
η∈B(η0 ,δ0 )

2

ǫn

(B.2)

By combining (B.2), Lemma 4, and the triangle inequality, we obtain kτ̂η̂ − τ0 k22,P = Op ((log n)2 n−2/3 ).

Proof of Lemma 5 (iii). To avoid heavy notation, we use the same notation as in the proof of
Lemma 3 (iii), but some notation is redefined here. Let
TI,K1 = {τ monotone non-decreasing on some interval I and 0 < τ < K1 }.
Assumption M2 guarantees 0 < C < τ0 < C < ∞. Similar to the proof of Lemma 3 (iii), we
calculate HB (ǫ, F̃ , k · kB,P ) and kf˜k2 , with F̃ = {f˜ = D −1 f : f ∈ F}, where the constantB,P

D > 0 is determined later. Define Iη∗ = (aL , aU ) with aL = inf x∈X ,η∈B(η0 ,δ0 ) x′ η and aU =supx∈X ,η∈B(η0 ,δ0 ) x′ η. Define





1
1
′
−
wh u :
Fn = fn (w, u) = I{x η > qn }τ (x′ η) τη (x′ η)τ ∈ TIη∗ ,K1 , η ∈ B(η0 , δ0 ),kτ − τ0 k22,P ≤ v 2 , h ∈ {1 : dim(w)},I(x′ η > qn )/τ (x) ≤ 1/K0

where wh is the h-th component of w. We set 2K0 = C, K1 = K2 log n, and v = K3 (log n)n−1/3
for some positive constants K2 and K3 .
36







,

By van der Vaart and Wellner (1996, Theorem 2.7.5), it holds for each ǫ ∈ (0, K1 ),HB (ǫ, TIη∗ ,K1 , k · kP ) ≤C ′ K1.
ǫ

Similarly to the univariate case, we can choose those bracket functions (τ L , τ U ), which satisfy
I{x′ η > qn }/τ L (x′ η) ≤ 1/K0 . Then, we defineI{x′ η > qn } U 1 ′ − 1 ′wh u τ (x η) τη (x η) f L (w, u) =1
I{x′ η > q }− τη (x1 ′ η) wh un
τ L (x′ η)I{x′ η > qn } L 1 ′ − 1 ′wh u τ (x η) τη (x η) f U (w, u) =
I{x′ η > q }1− τη (x1 ′ η) wh un
τ U (x′ η)

if wh u ≥ 0,
if wh u < 0,
if wh u ≥ 0,
if wh u < 0.

Note that (f L , f U ) is a bracket for f ∈ Fn . The bracket size is

kf˜U − f˜L k2B,P = kD −1 f U − D −1 f L k2B,P


Z
∞
k
X
1
1
1
′
I{x η > qn }−wh u dP (w, u)= 2
k!D k W ×R
τ L (x′ η) τ U (x′ η)
k=2
(
)
∞
X
1
Rk k!M0k−2 a0 (2K1 )k−2 U
L 2≤ 2kτ − τ kP
2k
k!D k
K
0
k=2
2 X


∞ 
2RM0 K1 k 2
R
ǫ ,
≤ 2a0
DK02
DK02
k=0

where the first inequality follows from Assumption M2 (where we can choose a0 , M0 > 1) and1/2
aI{x′ η>qn }≤ 1 . Setting D = 4M0 RK1 /K 2 yields kf˜U − f˜L kB,P ≤ K̃ǫ for K̃ = 0 , and thusL ′
τ (x η)0

K0

2M0 K1B̃
HB (ǫ, F̃ , k · kB,P ) ≤ ,ǫ

1/2

C 2 a0for B̃ =
.
2M0

(B.3)

Now we compute the Bernstein norm of f˜:


Z
∞
k
X
1
1
1
′
2
˜
I{x
η
>
q
}
−
kf kB,P = 2wudP (w, u)nhk!D k W ×R
τ (x′ η) τη (x′ η)
k=2
(
)
∞
X
1
Rk k!M0k−2 a0 (2K1 )k−22≤ 2kτ − τ0 kP
k!D k
K02k
k=2
2 X


∞ 
a0 1 2
2RM0 K1 k 2
R
v ≤
v ,
≤ 2a0
2
2
DK0
DK0
4M02 K12
k=0′

η>qn }
where the first inequality follows from I{xτ (x≤ K10 . This implies′ η)

kf˜kB,P ≤ B

v
,
K1

1/2

for B =

a0
.
2M0

(B.4)

Combining (B.3) and (B.4), the remaining steps are the same as those in the proof of Lemma
3 (iii).

37

P
B.3. Proof of Lemma 4. Recall for fixed η, we first obtain τ̂η = arg minτ ∈M n1 ni=1 {Ûi2 −P
τ (Xi′ η)}2 and then obtain η̂ by η̂ = argmin|| n1 ni=1 Xi′ {Ûi2 − τ̂η (Xi′ η)}||2 . We denote E[X|X ′ η =ηx′ η] by E[X|x′ η]. The proof is similar to the ones in BGH and Balabdaoui and Groeneboom(2021) except that we need to handle the influence of the estimated dependent variables Ûi2 .
The proof of consistency of η̂ is similar to pp.16-17 of BGH-supp. By a similar argument in
Balabdaoui and Groeneboom (2021, Lemma 3.2), under Assumptions M1-M3, we have
n

n

i=1

i=1

1X
1X ′ 2
Xi {Ûi − τ̂η (Xi′ η)} =(Xi − E[X|Xi′ η]){Ûi2 − τη (Xi′ η)} + op (n−1/2 ),nnfor each η, where we also use (B.2). Thus, it holdsn

1X
Xi {Ûi2 − τ̂η̂ (Xi′ η̂)}
n
i=1

n

1X
= min
Xi {Ûi2 − τ̂η (Xi′ η)}
η
n
≤ min
η

1
n

i=1
n
X
i=1

(Xi − E[X|Xi′ η]){Ûi2 − τη (Xi′ η)} + op (n−1/2 ) .

The leading term inside the norm k·k of the last expression does not depend on the potentiallynon-smooth τ̂η ; it is a smooth function of η. Thus, under standard conditions for the method ofP
moments, we have minη n1 ni=1 (Xi − E[Xi |Xi′ η]){Ûi2 − τη (Xi′ η)} = 0, and
n

op (n

−1/2

)

=
=
=

1X
Xi {Ûi2 − τ̂η̂ (Xi′ η̂)}
n
i=1
n
X

1
(Xi − E[X|Xi′ η̂]){Ûi2 − τ̂η̂ (Xi′ η̂)} + op (n−1/2 + (η̂ − η))ni=1Z(x − E[X|x′ η̂]){û2 − τ0 (x′ η0 )}d(Pn − P )(x, û)Z+ (x − E[X|x′ η̂]){û2 − τη̂ (x′ η̂)}dP (x, û) + op (n−1/2 + (η̂ − η))=: I + II + op (n−1/2 + (η̂ − η)),(B.5)

where the second equality follows from similar arguments to pp.18-20 of BGH-supp and (B.2),
and the third equality follows from a similar argument in pp.21-23 of BGH-supp.
Let Û (w, u) = u − w′ (θ̂OLS − θ) andê(w, u) := Û (w, u)2 − u2 = −2w′ (θ̂OLS − θ)u + {w′ (θ̂OLS − θ)}2 .

(B.6)

For I, we haveI =Z(x − E[X|x′ η̂]){u2 + ê(w, u) − τ0 (x′ η0 )}d(Pn − P )(w, u)

Z(x − E[X|x′ η0 ]){u2 − τ0 (x′ η0 )}d(Pn − P )(x, u)Z+ (x − E[X|x′ η̂])ê(w, u)d(Pn − P )(w, u) + op (n−1/2 )Z=(x − E[X|x′ η0 ]){u2 − τ0 (x′ η0 )}d(Pn − P )(x, u) + op (n−1/2 ),
=

38

(B.7)

where the second equality follows from p.21 of BGH-supp, and the third equality follows fromthe facts that (a) θ̂OLS − θ = Op (n−1/2 ), (b) ê(w, u) is a parametric function of w and u in achanging class indexed by θ̂OLS (see (B.6)), so its ǫ-entropy is of order log(1/ǫ) ≤ 1/ǫ (see, e.g.,

Example 19.7 of van der Vaart and Wellner, 2000), and (c) similar arguments in pp.22-23 of
BGH-supp. By Lemma 17 of BGH-supp we haveτη (x′ η) = τ0 (x′ η0 ) + (η − η0 )(x − E[X|X ′ η0 = x′ η0 ])τ0′ (x′ η0 ) + op (η − η0 ).

(B.8)

For II, observe that
Z
Z
′
2
′
II =
(x − E[X|x η̂]){u − τη̂ (x η̂)}dP (x, u) + (x − E[X|x′ η̂])ê(w, u)dP (w, u)Z=(x − E[X|x′ η0 ])(x − E[X|X ′ η0 = x′ η0 ])τ0′ (x′ η0 )dP (x) (η̂ − η0 )Z
+ (x − E[X|x′ η̂])ê(w, u)dP (w, u) + op (η̂ − η0 )

Z
′
′
′ ′
(x − E[X|x η0 ])(x − E[X|x η0 ])τ0 (x η0 )dP (x) (η̂ − η0 ) + Op (n−1/2 ) + op (η̂ − η0 )== B(η̂ − η0 ) + Op (n−1/2 ) + op (η̂ − η0 ),(B.9)

where the third equality follows from (B.8) and (E[X|x′ η̂] − E[X|x′ η0 ])(η̂ − η0 ) = op (η̂ − η0 ), thefourth equality follows from θ̂OLS − θ = Op (n−1/2 ) and the definition of B in Assumption M6.

Combining (B.5), (B.7), and (B.9), we haveZ−η̂ − η0 = B(x − E[X|x′ η0 ]){u2 − τ0 (x′ η0 )}d(Pn − P )(x, u) + Op (n−1/2 ) + op (n−1/2 + (η̂ − η)),where B − is the Moore-Penrose inverse of B (see p.17 of BGH for more details). Therefore,

we have η̂ − η0 = Op (n−1/2 ). This result, combined with (B.8) and Assumptions M1 and M2,

implies τη̂ (a) − τ0 (a) = Op (n−1/2 ) and kτη̂ − τ0 k2,P = Op (n−1/2 ).
References[1] Acemoglu, D. and P. Restrepo (2017) Secular stagnation? The effect of aging on economic growth in the ageof automation, American Economic Review, 107, 174-179.
[2] Ayer, M., Brunk, H. D., Ewing, G. M., Reid, W. T. and E. Silverman (1955) An empirical distributionfunction for sampling with incomplete information, Annals of Mathematical Statistics, 26, 641-647.
[3] Babii, A. and R. Kumar (2023) Isotonic regression discontinuity designs. Journal of Econometrics, 234(2),371-393.
[4] Balabdaoui, F., Durot, C. and H. Jankowski (2019) Least squares estimation in the monotone single index
model, Bernoulli, 25, 3276-3310.
[5] Balabdaoui, F., Groeneboom, P. and K. Hendrickx (2019) Score estimation in the monotone single-index
model, Scandinavian Journal of Statistics, 46, 517-544.
[6] Balabdaoui, F. and P. Groeneboom (2021) Profile least squares estimators in the monotone single index
model, in Advances in Contemporary Statistics and Econometrics, pp. 3-22, Springer.
[7] Barlow, R. and H. Brunk (1972) The isotonic regression problem and its dual, Journal of the American
Statistical Association, 67, 140-147.
[8] Bickel, P. J. (1978) Using residuals robustly I: Tests for heteroscedasticity, nonlinearity, Annals of Statistics,
6, 266-291.
39

[9] Box, G. E. and W. J. Hill (1974) Correcting inhomogeneity of variance with power transformation weighting,Technometrics, 16, 385-389.
[10] Carroll, R. J. (1982) Adapting for heteroscedasticity in linear models, Annals of Statistics, 10, 1224-1233.
[11] Chetverikov, D. (2019) Testing regression monotonicity in econometric models, Econometric Theory, 35,
729-776.
[12] Cragg, J. G. (1983) More efficient estimation in the presence of heteroscedasticity of unknown form, Econometrica, 51, 751-764.
[13] Cragg, J. G. (1992) Quasi-Aitken estimation for heteroskedasticity of unknown form, Journal of Econometrics, 54, 179-201.
[14] Dümbgen, L. and V. G. Spokoiny (2001) Multiscale testing of qualitative hypotheses, Annals of Statistics,
29, 124-152.
[15] Fang, B., Guntuboyina, A., and Sen, B. (2021) Multivariate extensions of isotonic regression and total
variation denoising via entire monotonicity and Hardy–Krause variation. Ann. Statist. 49(2): 769-792.
[16] Ghosal, S., Sen, A. and A. van der Vaart (2000) Testing monotonicity of regression, Annals of Statistics, 28,1054-1082.
[17] Grenander, U. (1956) On the theory of mortality measurement. II., Skand. Aktuarietidskr, 39, 125-153.
[18] Groeneboom, P. and G. Jongbloed (2014) Nonparametric Estimation under Shape Constraints, Cambridge
University Press.
[19] Hall, P. and N. E. Heckman (2000). Testing for monotonicity of a regression mean by calibrating for linear
functions, Annals of Statistics, 28, 20-39.
[20] Hansen, A. H. (1939) Economic progress and declining population growth, American Economic Review, 29,
1-15.
[21] Hsu, Y. C., Liu, C. A. and X. Shi (2019) Testing generalized regression monotonicity, Econometric Theory,35, 1146-1200.
[22] Jobson, J. D. and W. A. Fuller (1980) Least squares estimation when the covariance matrix and parameter
vector are functionally related, Journal of the American Statistical Association, 75, 176-181.
[23] Kim, J. and D. Pollard (1990) Cube root asymptotics, Annals of Statistics, 18, 191-219.
[24] Kulikov, V. N. and Lopuhaä, H. P. (2006) The behavior of the NPMLE of a decreasing density near theboundaries of the support. Annals of Statistics, 34(2), 742-768.
[25] Matzkin, R. L. (1994) Restrictions of economic theory in nonparametric methods, in Engle R. F. and D. L.
McFadden (eds.), Handbook of Econometrics, vol. IV, pp. 2523-2558, Elsevier.
[26] Mincer, J. (1974) Schooling, Experience, and Earnings. Human Behavior & Social Institutions, No. 2.
[27] Newey, W. K. (1993) Efficient estimation of models with conditional moment restrictions, in Maddala, G. S.,
Rao, C. R. and H. D. Vinod (eds.), Handbook of Statistics, vol. 11, pp. 419-454, Elsevier.
[28] Rao, B. P. (1969) Estimation of a unimodal density, Sankhyā, A 31, 23-36.
[29] Rao, B. P. (1970) Estimation for distributions with monotone failure rate, Annals of Mathematical Statistics,
41, 507-519.
[30] Robinson, P. M. (1987) Asymptotically efficient estimation in the presence of heteroskedasticity of unknown
form, Econometrica, 55, 875-891.
[31] Ruud, P. A. (2000) An introduction to classical econometric theory. OUP Catalogue.
[32] van de Geer, S. (2000) Empirical Processes in M−estimation. Cambridge University Press.
[33] van der Vaart, A. W. (2000) Asymptotic Statistics, Cambridge University press.
[34] van der Vaart, A. W. and J. A. Wellner (1996) Weak Convergence and Empirical Processes, Springer.
[35] Woodroofe, M. and Sun, J. (1993) A penalized maximum likelihood estimate of f (0+) when f is nonincreasing. Statistica Sinica, 501-515.
[36] Wooldridge, J. M. (2010) Econometric Analysis of Cross Section and Panel Data, 2nd ed., MIT Press.
[37] Wooldridge, J. M. (2013) Introductory econometrics: A modern approach, 5th Ed. Cengage learning.

40

School of Social Sciences, Waseda University, 1-6-1 Nishiwaseda, Shinjuku-ku, Tokyo 169-8050,
Japan.
Email address: yarai@waseda.jp
Department of Economics, London School of Economics, Houghton Street, London, WC2A 2AE,
UK.
Email address: t.otsu@lse.ac.uk
Department of Economics, University of Mannheim, L7 3-5, 68161, Mannheim, Germany.
Email address: mengshan.xu@uni-mannheim.de41