Journal of Econometrics 203 (2018) 344–358

Contents lists available at ScienceDirect

Journal of Econometrics
journal homepage: www.elsevier.com/locate/jeconom

Consistent estimation of linear regression models using
matched data✩
Masayuki Hirukawa a , Artem Prokhorov b,c,d, *aFaculty of Economics, Setsunan University, 17-8 Ikeda Nakamachi, Neyagawa, Osaka 572-8508, Japan
Discipline of Business Analytics, Business School, University of Sydney, H04-499 Merewether Building, Sydney, NSW 2006, AustraliacSt. Petersburg State University, St. Petersburg, Russian FederationdInnopolis University, Kazan, Russian Federationb

articleinfo

Article history:
Received 6 June 2015
Received in revised form 16 April 2017
Accepted 4 July 2017
Available online 6 January 2018
JEL classification:
C13
C14
C31

a b s t r a c t
Economists often use matched samples, especially when dealing with earnings data where a number ofmissing observations need to be imputed. In this paper, we demonstrate that the ordinary least squares
estimator of the linear regression model using matched samples is inconsistent and has a non-standard
convergence rate to its probability limit. If only a few variables are used to impute the missing data, then it
is possible to correct for the bias. We propose two semiparametric bias-corrected estimators and explore
their asymptotic properties. The estimators have an indirect-inference interpretation, and they attain theparametric convergence rate when the number of matching variables is no greater than four. Monte Carlosimulations confirm that the bias correction works very well in such cases.
© 2017 Elsevier B.V. All rights reserved.

Keywords:
Bias correction
Indirect inferenceLinear regression
Matching estimation
Measurement error bias

1. Introduction
Suppose that we are interested in estimating a linear regressionmodel
Y = β0 + X1′ β1 + X2′ β2 + Z ′ γ + u := W ′ θ + u, E ( u| W ) = 0, (1)using a random sample, where X1 ∈ Rd1 , X2 ∈ Rd2 and Z ∈ Rd3 .
The reason for distinguishing between the regressors X1 , X2 andZ will become clear shortly. In addition, while d1 = 0 is allowed, d2 , d3 > 0 must be the case in our setup. When W =✩ The authors would like to thank the editor Oliver Linton, an associate editor, three anonymous referees, Alberto Abadie, Bruce Hansen, Hidehiko Ichimura,
Seojeong Lee, Arthur Lewbel, Adrian Pagan, and the participants of NZESG2014,
AMES2014, ESAM2014, ESEM2014, CORE@50, CFE2016, IAAE2017, and seminars
at the University of New South Wales, the University of Sydney, the University ofTokyo, and the Development Bank of Japan for their constructive comments andsuggestions. Financial support through grants from Japan Society of the Promotionof Science (M. Hirukawa, grant numbers 23530259 and 15K03405) and the Russian
Science Foundation (A. Prokhorov, Project No. 16-18-10432) for various and nonoverlapping parts of this research is gratefully acknowledged.
Corresponding author at: Discipline of Business Analytics, Business School, University of Sydney, H04-499 Merewether Building, Sydney, NSW 2006, Australia.
E-mail addresses: hirukawa@econ.setsunan.ac.jp (M. Hirukawa),
artem.prokhorov@sydney.edu.au (A. Prokhorov).

*

https://doi.org/10.1016/j.jeconom.2017.07.006
0304-4076/© 2017 Elsevier B.V. All rights reserved.

)′

1, X1′ , X2′ , Z ′ ∈ Rd+1 , where d := d1 + d2 + d3 , is exogenous
and a single random sample of (Y , X1 , X2 , Z ) can( be obtained,)the′ordinary least squares (OLS) estimator of θ = β0 , β1′ , β2′ , γ ′ isconsistent.
In reality, however, we often face the problem that (Y , X1 , X2 , Z )cannot be taken from a single data source. It is not uncommon
that economists who use survey data for empirical analysis mustcollect all necessary variables from more than one source. Examples include Lusardi (1996), Björklund and Jäntti (1997), Currie
and Yelowitz (2000), Dee and Evans (2003), Borjas (2004), Bover(2005), Fujii (2008), Bostic et al. (2009), and Murtazashvili et al.
(2015), to name a few. Ridder and Moffitt (2007) provide an
excellent survey. This is the setting in which we are interested.
Specifically, suppose that instead of observing a complete dataset (Y , X1 , X2 , Z ), we have the following two overlapping subsets
of data, (Y , X1 , Z ) and (X2 , Z ), i.e., some of the regressors are notavailable in the initial data set, where the initial data set is the one
containing observations on the dependent variable along with afew other regressors. In such a setting, it is natural to constructa matched data set via exploiting the proximity of the common
regressor(s) Z across the two samples. This is often called ‘‘probabilistic record linkage’’. Here are two examples of the setting.

(M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

Example 1 (Earnings Data). Matching is currently used for imputing missing records of earnings in important economic data
sets. For example, the U.S. Current Population Survey (CPS) files
use the so-called ‘‘hot deck imputation’’ procedure of the Census(see, e.g., Little and Rubin, 2002; Hirsch and Schumacher, 2004;
Bollinger and Hirsch, 2006), which allocates to nonrespondents
the reported earnings of a matched respondent who has similar
recorded attributes.1 The share of imputed values is as high as
30%. The resulting earnings data have been used to uncover muchof what is known about the labor market dynamics and outcomes.
Example 2 (Returns to Schooling). Let Y denote (the logarithmof) earnings, X1 individual characteristics, X2 ability measured bytest scores, and Z education. Although (Y , X1 , Z ) is available in the
Panel Study of Income Dynamics (PSID), for instance, it is often thecase that (X2 , Z ) can be found only in a different, psychometric data
set. Utilizing the proximity of the common variable Z , we must
construct a matched data set of (Y , X1 , X2 , Z ).
There are many algorithms that can be used to construct
matched data sets (see, e.g., Smith and Todd, 2005; Ridder andMoffitt, 2007). We focus on the nearest neighbor matching (NNM)
because of its simplicity and wide use. Abadie and Imbens (2006,2012) use it in the context of the average treatment effect (ATE) estimation. Chen and Shao (2001) and Shao and Wang (2008) study
the problem of variance estimation after a nearest neighbor-based
imputation. The NNM can be used as a building block in construction of more complicated matching algorithms, most notably thesingle index or propensity score matching, but we do not pursue
these extensions here.
We demonstrate that the OLS estimator from the regression(1) using NNM-based matched samples is inconsistent. The sourceof the inconsistency is a non-vanishing bias term, which can beviewed as a measurement error bias stemming from replacing
unobservable X2 with a proxy in the matched data. In this sense,
the paper is related to the literature on the classical problem
of generated regressors and missing data (see, e.g., Pagan, 1984;
Prokhorov and Schmidt, 2009). Moreover, we show that the rate ofconvergence to the probability limit of OLS depends on the numberof common, matching variables and the divergence patterns of two
sample sizes.
In line with these findings, we propose two semiparametricbias-corrected estimators. The first, one-step estimator is shown toattain the parametric convergence rate for the cases with at most
two matching variables. The second estimator attempts to remedy
the curse of dimensionality with respect to the number of matchingvariables. It is a two-step estimator, and in the second step it
eliminates the second-order bias due to the so-called matching
discrepancy (Abadie and Imbens, 2006) asymptotically in a similarmanner to the one studied by Abadie and Imbens (2011). It isdemonstrated that the estimator can achieve parametric convergence when the number of matching variables is four or less. Bothestimators can be also interpreted as indirect inference estimators(Gouriéroux et al., 1993; Smith, 1993) in the sense that they can beobtained by taking the probability limit of the OLS estimator fromthe regression (1) as the ‘‘binding’’ function.
The paper contributes to three important areas. First, we provide new asymptotic results for regression analysis using matched
data. In particular, we explicitly handle the issue of biases due tomatching errors, which has been often ignored in the literature as if
1 The distinction between hot and cold deck imputation seems to primarily refer
to which sample (of punch cards) to use for matching, a current sample (hot)or an earlier sample (cold). Hence, hot deck imputation often means imputationof missing values of an existing variable, whereas cold deck imputation means
imputation of entire missing variables. In this respect, this paper may be closer tocold rather than hot deck imputation.

345

there were no mismatches; see Ridder and Moffitt (2007, p. 5480)
for a discussion and Bover (2005) and Bostic et al. (2009) for regression analysis using matched data. Available results are limited
to the case of matching in the ATE estimation. For example, Abadie
and Imbens (2006) show that when there is only one matchingcovariate, the bias in NNM-based matching estimators of the ATE
may be asymptotically ignored; they attain the parametric convergence rate in that case. To the best of our knowledge, bias-corrected
estimation using matched data and the convergence properties of
estimators within the framework of regression analysis have notbeen explored in the literature before.
Second, the estimation theory we develop provides guidanceon repeated survey sampling when some covariates are found to
be completely or partially missing after the initial survey. Our
theory suggests (approximately) how many observations should
be collected in a follow-up survey and how to estimate the linear
regression model of interest consistently using the matched datafrom two surveys.
Finally, the paper offers an alternative to some well-known
estimation methods based on two samples. A number of suchmethods have been designed within the framework of instrumental variables (IV) or generalized method of moments (GMM)
estimation, where we can construct required moments from thetwo samples individually so that no matching is required (e.g.,
Angrist and Krueger, 1992, 1995; Arellano and Meghir, 1992; Inoueand Solon, 2010; Murtazashvili et al., 2015). These approaches are
not applicable to the setting of a linear regression where some
regressors are missing and two-sample moment based estimation
is infeasible.
Throughout we assume that the two samples jointly identify
the regression models. There are other two-sample estimators thatcover the cases where the first sample alone identifies the models
and the second sample is used for efficiency gains (see, e.g., Imbens
and Lancaster, 1994; Hellerstein and Imbens, 1999). These are not
the settings we consider.
The remainder of this paper is organized as follows. Section 2shows inconsistency of the OLS estimation from the regression(1) using matched samples. Section 3 proposes two bias-corrected
estimators and explores their convergence properties. We also discuss consistent estimation of their asymptotic covariance matrices.
Section 4 contains selected Monte Carlo simulations, examining
how the bias correction works in finite samples. As an empiricalexample, in Section 5, we apply the bias-corrected two-sample estimation to a version of Mincer’s (1974) wage regression. Section 6
concludes with a few questions for future research. All proofs are
given in Appendix A.2{ The( paper)}1/2 adopts the following notational conventions: ∥A∥ =tr A′ Ais the Euclidean norm of matrix A; 1 {·} denotes anindicator function; 0p×q signifies the p × q zero matrix, where thesubscript may be suppressed if q = 1; and the symbol > applied
to matrices means positive definiteness.
2. Inconsistency of OLS estimation using matched samples
2.1. Setup
In order to explain how a matched sample is constructed, weneed more notations. Denote the two random samples by S1 andS2 . Also let n and m be sample sizes of S1 and S2 , respectively. Then,the two samples can
as S1 = S1n = {(Yi , X1i , Zi )}ni=1{( be expressed)}mand S2 = S2m = X2j , Zj j=1 . A natural way of matching basedon Z is to use the NNM based on some metric. For a vector x
2 An online Supplement contains some lengthy derivations and comprehensive
simulation results. GAUSS codes implementing the estimators are also available
from the authors upon request.

346

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

and some(symmetric)1/2 matrix A > 0, a vector norm is denoted by
x′ Ax. While there may be numerous choices of A,
following Abadie
andImbens (2011), we adopt{} the Mahalanobis∥x∥A =(1/N )

metric AM =

∑N (
i=1

Zi − Z̄)(Zi − Z̄

1
ized Euclidean metric ANE = diag A−
M

(

)−1

)′ −1

and the normal-

∑N

A

j=1

Also let JK (i) = {j1 (i) , . . . , jK (i)} denote the set of indices for the
first K matches for the unit i. The NNM constructs the matched dataset
S = Sn =θ̂OLS := Q̂W−1 R̂W :=

)−1

n

1∑
n

Wi,j(i) Wi′,j(i)

i=1

n

1∑
n

Wi,j(i) Yi

i=1

is referred to as the matched-sample OLS (MSOLS) estimator hereinafter.
2.2. Regularity conditions

m
∑


 }
{
1 Zj − Zi  ≤ Zj (i) − Zi  = k.
k(

, where N := n + m andZ̄ = (1/N ) i=1 Zi .
Furthermore, let jk (i) be the index of the kth match in S2 to theunit i in S1 , i.e., for each i ∈ {1, . . . , n}, jk (i) satisfies
A

We start our analysis from running OLS for the regression of Yi
on Wi,j(i) . The OLS estimator

Yi , X1i , X2j1 (i) , . . . , X2jK (i) , Zi , Zj1 (i) , . . . , ZjK (i)

{(

)}n
i=1

.

In what follows, we develop the asymptotic theory of estimation of θ in the regression (1) as n and m diverge while K is fixed. Allof the estimation theory, including the bias-corrected estimation
methods and their convergence properties, is new to the literature.
All of the results in the paper, including the inconsistency proofof the MSOLS estimator and the consistency proof of the new biascorrected estimators of θ , require the following conditions.

We also define
X2j(i) :=

( ) ∑
1
K

j∈JK (i)X2j and Zj(i) :=

( ) ∑
1
K

Assumption 1. Two random samples (S1 , S2 ) = (S1n , S2m ) are
drawn independently from the joint distribution of (Y , X1 , X2 , Z )with finite fourth-order moments.

Zj .

j∈JK (i)It is worth noting that X2 is missing entirely but only from thefirst sample. When the problem is considered in the context of
the imputed sample, only the values corresponding to the firstsample are missing. Thus formally, the problem can be viewed
as both value imputation and variable imputation. However, inwhat follows we view the problem as missing variable (rather than
missing values) imputation.
In our NNM, the number of matches K remains fixed, as in
Abadie and Imbens (2006). While it is possible to achieve consistency as in the K -nearest neighbor (K -NN) method by letting K
diverge at a slower rate than n and m, there are two reasons why
we keep K fixed. First, this is what is done in practice. In many
applications, the NNM is implemented with small values of K , andK = 1 (i.e., NNM with a single match) is often chosen even forlarge n and m. Second, it can be argued that the analysis with fixed
K provides a better approximation to the finite sample behavior of
the estimator than under K → ∞.
A few additional remarks on NNM are in order. First, matching is
made with replacement, and each element of the matching vector
Z is assumed to be continuous. Hence, our setting can be viewed
as a foundation for more complicated methods of kernel-basedmatching (see, e.g., Busso et al., 2014). Second, matching withreplacement, allowing each unit to be used as a match more thanonce, seems to be standard in the econometric literature. Third,
inclusion of discrete matching variables with a finite number of
support points does not affect the subsequent asymptotic results
but raises the question of how to treat ties. For simplicity, we ignore
ties in the NNM, which happen with probability zero as long as Z
is continuous.
Throughout( it is assumed )that we estimate θ by regressing Yi′′ ′on Wi,j(i) := 1, X1i′ , X2j(i) , Zi . It is possible to use Zj(i) in place
†

)′

′
′
of Zi and run the regression of Yi on Wi,j(i) := 1, X1i′ , X2j(i) , Zj(i) .
However, we focus exclusively on the former scenario becauseof the following two reasons. First, the two scenarios yield firstorder asymptoticallyequivalent results. To[( that)()′ ]′ see this, observe†Wi,j(i) = Wi,j(i) + 01×(d1 +d2 +1) Zj(i) − Zi
= Wi,j(i) + Op m−1/d3
by Lemma A.1, i.e., the second term serves merely as an extra
second-order bias term. It is noteworthy that the identification
condition is derived from the latter scenario. Second, as illustrated
in Section 4, bias-corrected estimators based on Wi,j(i) exhibits
better finite-sample properties.

(

Assumption 2. The matching variable Z is continuously distributed with a convex and compact support Z, and its density isbounded and bounded away from zero on Z.
Assumption 3.
(i) Theerror u satisfies E ( u| W ) = 0 and σu2 (W ) :=( ⏐regression)
E u2 ⏐ W ∈ (0, ∞).

[
]′
[
]′
g1 (Z )′ g2 (Z )′ := E ( X1 | Z )′ E ( X2 | Z )′ and
[
]′
[
]′
:= ( X1′ −)g1 (Z )′ X(2′ − g)2 (Z )′ . Then, Σ1 :=let η := η1′ η2′( ′)E η1 η1 > 0, Σ2 := E η2 η2′ > 0, E η1 η2′ = 0d1 ×d2 , and g2 (·)(ii) Let g (Z ) :=is a first-order Lipschitz continuous, strictly nonlinear function
on Z.
These regularity conditions are largely inspired by those in
the literature on semiparametric, partial linear regression models(e.g., Robinson, 1988; Yatchew, 1997), matching estimators for the
ATE (e.g., Abadie and Imbens, 2006), and regression estimation
based on two samples (e.g., Angrist and Krueger, 1992; Inoue andSolon, 2010). In particular, equivalents to Assumption 1 (the common distribution assumption) are often imposed in the literature(e.g., Assumption 3 of Abadie and Imbens, 2006; Assumption a ofInoue and Solon, 2010). This is a strong assumption which simplifies the subsequent derivations considerably. It implies that thematched sample S behaves as a pseudo-population, from which
the two samples are drawn. Assumption 2 plays a key role incontrolling the order of magnitude in the matching discrepancy.
Nonlinearity of g2 (·) in Assumption 3(ii) will be discussed in
Remark 1 in relation to identification. Zero correlation between
η1 and η2 in Assumption 3(ii) is also a key assumption. Because
we observe X1 and X2 in (two independent
samples, there seems to)
be no way to estimate E η1 η2′ consistently. Therefore we assume
uncorrelatedness between η1 and η2 .
2.3. Inconsistency of MSOLS
Our asymptotic analysis is built on reformulating the regression(1) in a ‘partial linear’-like format. A straightforward calculation
yieldsYi := Wi′,j(i) θ + λi,j(i) + ϵi,j(i) , i = 1, . . . , n,(2)

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

restrictions somewhat and/or improve efficiency. For example, inthe presence of an outside matching variable, g2 (·) can be allowed
to be linear. But we do not pursue this point here.

where

(
)
λi,j(i) = λ Zi , Zj(i) =

⎧
⎨

1 ∑

⎩

K

g2 (Zi ) −

⎫′
( )⎬

g2 Zj

j∈JK (i)β2 , and⎭⎞′⎛ϵi,j(i) = ui + ⎝η2i −

1 ∑
K

)′

η2j ⎠ β2 := ui + η2i − η2j(i) β2 .

(j∈JK (i)The reason why this is not exactly a partial linear model is that
there is a common regressor Zj(i) included in Wi,j(i) and λi,j(i) . In thisformulation, Wi,j(i) is employed as the regressor of the fully parametric part Wi′,j(i) θ . On the other hand, the semiparametric partλi,j(i) generates the second-order bias that will be discussed shortly,and thus it could be viewed as an analogue to the conditional bias
discussed in Abadie and Imbens (2006). A key difference from thepartial linear regression models studied in Robinson (1988) andYatchew (1997) is that the matched regressor X2j(i) is endogenous,i.e., X2j(i) and the composite error ϵi,j(i) are correlated. The theorem
below is established for the model in (2); it provides the probability
limit of θ̂OLS and its associated rate of convergence.
Theorem 1. If Assumptions 1–3 hold, then()()θ̂OLS = QW−1 PW θ + O m−1/d3 + Op n−1/2)(
as n, m → ∞, where QW := E Wi,j(i) Wi′,j(i) , PW := QW − (1/K ) Σ ,and Σ is a{ (d + 1) × (d + 1) block-diagonal
matrix of the form}
Σ := diag 0(d1 +1)×(d1 +1) , Σ2 , 0d3 ×d3 .
Remark 1. Basic identification assumptions for MSOLS follow fromthose of the standard OLS. Fundamentally, they require that η1and η2 are not in the linear span of each other and that X1 and X2
are not in the linear span of Z . As in the standard OLS, we needE(WW ′ ) to be of full rank. In our setting, the additional issue iswhether Q̂W and QW are invertible. While we implicitly assume
non-singularity of the former, the invertibility of the latter can be
examined explicitly.
The invertibilityof QW can be()′ analyzed from the regression of†
Yi on Wi,j(i) = 1, X1i , X2j(i) , Zj(i) . The result below turns out to be)′valid for the regression of Yi on Wi,j(i) = 1, X1i , X2j(i) , Zi , because(††′
i=1 Wi,j(i) Wi,j(i) and Q̂W are first-order asymptotically equivalent in that Q̂W
Q̂W †Op m−1/d3 by Lemma A.1.
†′†
Let QW †E Wi,j(i) Wi,j(i) . Then, a straightforward calculation
yields its determinant as QW †K −(d2 +d3 ) ΩX1 ,X1 ΩX2 ,X2

Q̂W † := (1/n)

:=

∑n
(

)=

+

347

(

Remark 2. Theorem 1 implies that MSOLS is inconsistent in general. The term (1/K ) Σ in PW , which is the source of inconsistency,
is generated by misspecifying the regression of Yi on Wi as the oneof Yi on Wi,j(i) , or equivalently, employing X2j(i) as a proxy of thelatent variable X2i . Therefore, the non-vanishing bias in MSOLS canbe thought of as a measurement error bias. The measurement error
interpretation is revisited in Section 2.4.
A straightforward calculation shows that the MSOLS estimatorof β2 is biased toward zero (an attenuation bias) in the limit
and that all other parameter estimates are also biased. Perhaps
surprisingly, the entire set of MSOLS estimates is inconsistent. Thisis contrary to the textbook regression with proxy variables where
only a part of the OLS estimates is biased (see, e.g., Wooldridge,
2013, Ch. 9). It is also easy to demonstrate that the bias in θ̂OLS
is a function of β2 , Σ2 and K . Therefore, the estimator would beconsistent if either (i) β2 = 0, i.e., X2 were irrelevant in the correctly specified model; or (ii) Σ2 = 0, i.e., X2 were a deterministic
function of Z .
Remark 3. The convergence rate of θ̂OLS is affected by the()
Op m−1/d3 term, which corresponds to the second-order bias
term λi,j(i) due to the matching discrepancy. The rate can be determined by three different divergence patterns of (n, m), namely,
n/m → κ ∈ (0, ∞), n/m → 0, and n/m → ∞ as n, m → ∞, andthere exists a curse of dimensionality with respect to the matchingvariable Z for each divergence pattern. ()−1
When n/m → κ , θ̂OLS = QW
PW θ + Op n− min{(1/2,1/d3 } . For d3 =)1, a central limit theorem (CLT) implies that

√

⏐⏐
⏐⏐

⏐
⏐

|⏐ | =⏐( )⏐⏐−1
× ⏐ΩZ ,Z − ΩZ ,X2 ΩX2 ,X2 ΩX2 ,Z ⏐, where Ωa,b := E ab′ − E (a) E (b)′between two random
vectors a and b. Since⏐⏐is the covariance⏐⏐⏐ΩX ,X ⏐ ⏐⏐ΩZ ,Z − ΩZ ,X Ω −1 ΩX ,Z ⏐⏐ is the determinant of the coX2 ,X2
2 2
2
2
(
)′
variance matrix of X2′ , Z ′ , |QW † | > 0 holds with no additional

restrictions. Hence, QW † is invertible. Furthermore, the identification of bias-corrected estimators that will be proposed in the next
section requires us to ensure non-singularity of PW † := QW † −(1/K ) Σ , where |PW † |, the determinant of PW † , can be obtainedby replacing( X2 in |)Q′ W † | with g2 (Z ). Considering the covariancematrix of g2 (Z )′ , Z ′ , we find that |PW † | > 0 if and only if g2 (·) is
strictly nonlinear, as assumed in Assumption 3(ii).
So far we have maintained the assumption that the vector ofcommon variables Z is employed for both matching and estimation. It is possible that at least one common variable is used exclusively for matching (and thus not included in the regression (1)).3
In this case the variable can be used to form yet another identification condition, which would allow us to relax our identification
3 We thank an anonymous referee for pointing out this possibility to us.

−1n θ̂OLS − QW
PW θ

has a normal limit. For d3 = 2, θ̂OLS is still n-convergent, but we
can only demonstrate asymptotic normality of θ̂OLS after subtracting the second-order bias(term, i.e., the best we can) do in this caseis to apply the CLT to

)

⏐
⏐

√

√

−1
n θ̂OLS − QW
PW θ − BOLS2 , wheren−1−1
BOLS2 := Q̂W
BRW 2 := Q̂W

1∑n

Wi,j(i) λi,j(i) .

i=1

These limiting distributions would reduce to the usual one of OLS
if a complete data set of (Y , X1 , X2 , Z ) were available. For d3 ≥ 3,
the convergence rate of θ̂OLS is slower than the parametric one, andit becomes slower as d3 increases.
()≤ 2. Hence, θ̂OLS =When n/m → 0, m−1/d3 =(o n−1/2 for d3 )−1
QW
PW θ + Op n−1/2 , and

(

)

√

−1
PW θ has a normal limitn θ̂OLS − QW

in this case. However, for d3 ≥ 3, the convergence rate of θ̂OLS canbe determined only if an extra divergence patternof (n, m) is√
imposed. For instance, when d3 = 3, θ̂OLS is n-convergent if()
n3 = O m2 and its convergence rate is a nonparametric one if
n3 /m2 → ∞.
√When n/m → ∞, a n-convergent of θ̂OLS( can be attained)()

only if d3 = 1 and n = O m2 . Moreover,√−1n θ̂OLS − QW
PW θ

has a normal limit when d3 = 1 and n/m2 → 0. On the other
hand, if d3 = 1 and n/m2 → ∞ or if d3 ≥ 2, then θ̂OLS =
(
)
−1
QW
PW θ + Op m−1/d3 , and the convergence bound m1/d3 is slower√
than n, regardless of whether or not the bound is sharp.

348

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

2.4. A measurement error interpretation
Before moving to our proposal for bias-corrected estimation,
it is helpful to consider the problem of imputation as a measurement error problem arising from using a proxy.4 Write the modelin (1) as
Y = β0 + X1 β1 + g2 (Z ) β2 + Z γ + e,′′′

where e := {X2 − g2 (Z )}′ β2 + u. Then, g2 (Z ) can be viewed as aproxy for X2 and if we could observe g2 (Z ) then the model could be
estimated by OLS as long as X1 is uncorrelated with {X2 − g2 (Z )}and g2 (Z ) is not in the linear span of Z .
However, g2 (Z ) is not observed and needs to be estimated.
There are two complications here. One is that we need to use
an estimator ĝ2 (Z ) based on another sample. The other is thatthe estimator uses matched values of X2 obtained using nearestneighbors of Z from the other sample, not the Z itself. Supposethat ĝ2 (Z ) is the estimate via the K -NN method for the moment.5
Rewriting the model as
Y = β0 + X1′ β1 + ĝ2 (Z )′ β2 + Z ′ γ + v,}′

where v := X2 − ĝ2 (Z ) β2 + u, we attempt to estimate thisregression by OLS. If ĝ2 (Z ) were estimated{from the same
} sample,then the correlation between ĝ2 (Z ) and X2 − ĝ2 (Z ) would benear zero because of orthogonality of g2 (Z ) and {X2 − g2 (Z )}. We
actually employ a different sample to estimate (or impute) g2 (Z ),and thus the correlation does not equal zero, which causes a nonnegligible bias in the OLS estimator. This can be interpreted as aclassical measurement error problem.
As is well known in the literature on measurement error problems, the bias of OLS can be corrected if the variance of themeasurement error can be obtained analytically, given that thematching discrepancy from K -NN is bounded. Our bias correction
methods in the next section basically follow this idea, although thenearest-neighbor algorithm that we use is intended only to find K
closest matches to Z and not to estimate g2 (Z ).

{

p

−1
of θ using strategy (1a) be θ̂(1a) = P̂W
R̂W , where P̂W → PW . Givenθ̂(1a) , we obtain the second-step estimator asθ̂(1b) := Q̂W

−1

(
R̂W +

1
K

Σ̂ θ̂(1a)

)

−1

(

= Q̂W

Id+1 +

1KΣ̂ P̂W−1)

R̂W , (3)where Σ̂ is a consistent estimate of Σ . Post-multiplying both sides−1−1of P̂W + (1/K ) Σ̂ = Q̂W by P̂W
yields Id+1 + (1/K ) Σ̂ P̂W=−1Q̂W P̂W . Substituting this into the right-hand side of (3) immediately establishes that θ̂(1b) = θ̂(1a) . Therefore, there is no pointin pursuing strategy (1b) separately; strategy (1b) is interesting
only if an alternative consistent estimator of θ (other than θ̂(1a) )is chosen.
Now we turn to the bias correction based on strategy (1a).
The idea behind the strategy comes from indirect inference (II)estimation by Gouriéroux et al. (1993) and Smith (1993). Take theprobability limit of θ̂OLS as the binding function b (θ), i.e., b (θ) =−1−1
QW
PW θ .6 Because PW
exists as discussed in Remark 1, the II
estimator can be built on the inverse mapping of θ̂OLS = b (θ),−1i.e., θ = PW
QW θ̂OLS . The interpretation then follows from replacing√
PW with its n-consistent estimator P̂W and regarding R̂W as a‘sample analog’ of QW θ̂OLS . Accordingly, we call this estimation
method the matched-sample indirect inference (MSII) estimation.
We formally define the MSII estimator as−1θ̂II := P̂WR̂W ,which has been called θ̂(1a) before.7
Our remaining task is to deliver a consistent estimator of PW .
Obviously, Q̂W is a natural estimator of}{ QW . Furthermore, it turns
out that when estimating Σ = diag 0(d1 +1)×(d1 +1) , Σ2 , 0d3 ×d3 ,we can do without a nonparametric estimation of g2 (·). To do so,
we first reorder S2 with respect to Z by the following recursion:
1. Define Z(1) as the observation that has the smallest firstelement, i.e., (1) = arg min1≤j≤m Zj1 .

2. For j = 2, . . . , m, choose (j) = arg minj̸=(1),...,(j−1) Zj − Z(j−1) .8

3. Bias-corrected estimation

Given the reordered sample S2 =consistently estimated by

3.1. One-step bias-corrected estimator

Σ̂2 =

1

m
∑

2 (m − 1)j=2}m

X2(j) , Z(j) j=1 , Σ2 can be{∆X2(j) ∆X2′ (j) ,(4)This section develops bias-corrected estimation of θ . As suggested by the proof of Theorem 1, inconsistency of MSOLS comesppfrom the fact that Q̂W → QW whereas R̂W → PW θ ={QW − (1/K ) Σ } θ . Therefore, the non-vanishing bias in MSOLS canbe eliminated if either

where ∆X2(j) := X2(j) − X2(j−1) . This is known as the differencebased variance estimator; see von Neumann (1941) and Rice
(1984) for univariate and Yatchew (1997) and Horowitz andSpokoiny (2001) for multivariate cases. In the end, the estimatorof PW is given by(1a) the denominator Q̂W is replaced by a consistent estimator of
PW with the numerator R̂W left unchanged; or(1b) an extra term consistent for (1/K ) Σ θ is added to R̂W withQ̂W held as it is.

P̂W := Q̂W −Bias correction in each strategy is semiparametric in that a
consistent estimate of Σ2 (covariance matrix of the nonparametric regression error η2 ) is required. Moreover, implementing (1b)inevitably leads to a two-step estimation with an initial consistentestimate of θ plugged in. However, if the plug-in estimator is theone using strategy (1a), then the two-step estimation will produce
a numerically identical result. To see why, let an initial estimator
4 We thank an anonymous referee for suggesting this interpretation to us.
5 We adopt a power-series approximation to estimate g (Z ); see Section 3.3 for
2details.

1KΣ̂ = Q̂W −

1
K

{

}

diag 0(d1 +1)×(d1 +1) , Σ̂2 , 0d3 ×d3 .

6 Typically the binding function is unknown, and it must be approximated via
simulations. However, when the function has a closed form, there is no need for
simulations; see Carrasco and Florens (2002) for another example. As suggested bya referee, this also gives the estimator the interpretation of a classical minimum
distance estimation based on the distance b(θ ) − θ̂OLS but we do not pursue thispoint further.
7 The estimator θ̂ also has a method-of-moment interpretation, where theIImoment is
E Wi,j(i) ϵi,j(i) = −

(

)

1
K

Σ θ.

From the viewpoint of likelihood-based methods MSII may leave some information (or moment restrictions) unused, and thus there may be room for efficiency
improvement. But pursuing this point is beyond the scope of this paper.
{ }8 If Z is a scalar, then the recursion reduces to rearranging Z m in an ascending
j j=1

order Z(1) ≤ · · · ≤ Z(m) .

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

3.2. Convergence properties of the MSII estimator+κ

We first provide a consistency result for the MSII estimator θ̂II .
The result holds regardless of the number of matching variables d3
and of the divergence patterns of (n, m).
p

Theorem 2. If Assumptions 1–3 hold, then θ̂II → θ as n, m → ∞.
Next we establish asymptotic normality of θ̂II after correcting
for the two bias terms, namely, BMD (second-order bias term dueto the matching discrepancy) and BΣ (third-order bias term due tothe nonparametric estimation of Σ ). Before proceeding, we makean additional assumption. Like Assumption c of Inoue and Solon(2010), Assumption 4 makes derivations of asymptotic variancesin the limiting distributions of θ̂II easier.
Assumption 4. In the nonparametric regression X2 = g2 (Z ) + η2 ,g2 (Z ) and η2 are independent, and third-order moments of η2 arezeros.
The next theorem establishes the limiting distributions of θ̂II
under a variety of divergence patterns of (n, m). Again, these results
hold regardless of the number of matching variables d3 .
Theorem 3. If Assumptions 1–4 hold, then, as n, m → ∞,

) d
⎧√ (
⎪
n
θ̂
−
θ
−
B
−
B
→ N (0, VI )
II
MD
Σ
⎪
⎪
(
)
⎪
⎪
−
1
−
1
⎪
:=( N 0, PW Ω PW⎪⎪
) ifd n/m → κ ∈ (0, ∞)⎪⎨√n θ̂II − θ − BMD − BΣ → N (0, VII ),⎪ := N (0, P −1 Ω11A P −1 ) if n/m → 0
⎪
⎪
W
W
) d
⎪
⎪√ (
⎪
⎪
m θ̂II − θ − BMD − BΣ → N (0, VIII )
⎪
⎪
⎩
(
)
−1
−1
:= N 0, PWΩ22 PW
if n/m → ∞∑n−1where BMD = 1/nP̂W
i=1 Wi,j(i) λi,j(i) and⎧⎨−1 1
BΣ = P̂W
diag 0(d1 +1)×(d1 +1) ,
⎩
K
⎫
m
⎬
∑
( )′
( )
1
∆g2 Z(j) ∆g2 Z(j) , 0d3 ×d3 θ⎭
2 (m − 1)
j=2

349

[

( ′
)
β2 Σ2 β2 E (W ) E (W )′
{
}]
)
(
1
1
+ 2 diag 0(d1 +1)×(d1 +1) , β2′ Σ2 β2 Vg2 + Ψ , 0d3 ×d3 .
K

2This suggests that VI is greater than VII in the positive definite
sense, i.e., amassing more information in the form of a quicker
growing m provides asymptotic efficiency gains.
Another relative efficiency implication is that when we include
an extra matching variable, i.e., d3 increases, the additional variability captured by PW leads to a smaller VIII (in the positive definite
sense) since the non-zero elements of Ω22 remain unchanged. Thisis similar to the effect of adding a regressor to a regression. The
effect of a larger d3 on VI and VII is unclear as Ω11A and Ω will bealso inflated by the inclusion.
More importantly, the limiting distributions of θ̂II vary acrossdivergence patterns of (n, m). The convergence rate of θ̂II is deter√
mined by the sample size of the smaller sample. The rate is nwhen n/m → κ , where κ ∈ [0, ∞), i.e., S2 is no smaller than S1
in order of magnitude. In contrast, when n/m → ∞√or S1 is (much√ )n .
larger than S2 , the best possible rate slows down to m = o
The slow rate can be thought of as the price paid for estimating θ
using a considerably smaller sample S2 via the NNM. As a result of
the bias correction, the order of magnitude in the estimation errorof Σ2 dominates.
While the MSII estimator is consistent, its convergence is affected by the asymptotic bias terms. For a large d3 the second-orderbias term BMD dominates and the convergence rate becomes inferior. This is the curse of dimensionality of NNM that is commonly
observed in other applications. With regards to the ATE estimation,
Abadie and Imbens (2006, Corollary 1), for instance, show thatthe matching discrepancy bias can be safely ignored only whenmatching is made on a single variable.
A similar result applies in our setting for small values of d3 . The
next corollary illustrates such cases.
Corollary 1. If Assumptions 1–4 hold, then, as n, m → ∞,

⎧√ (
) d
⎪
n θ̂II − θ → N (0, VI ) if n/m → κ ∈ (0, ∞) and d3 = 1
⎪
⎪
⎨√ (
) d
,
n θ̂II − θ → N (0, VII ) if n/m → 0 and d3 = 1, 2
⎪
(
) d
⎪
⎪
⎩√m θ̂II − θ →N (0, VIII ) if n/m → ∞ and d3 = 1

are the two bias terms,)√ (′κ Ω12 + Ω12+ κ Ω22)√ (′:= (Ω11A + Ω11B ) + κ Ω12 + Ω12+ κ Ω22 ,
{(
)(
)′ }
1
1
Ω11A := E
Wi,j(i) ϵi,j(i) + Σ θ
Wi,j(i) ϵi,j(i) + Σ θ
,
K
K
[
)
(
Ω11B := κ β2′ Σ2 β2 E (W ) E (W )′Ω := Ω11 +]()}+ 2 diag 0(d1 +1)×(d1 +1) , β2′ Σ2 β2 Vg2 + Ξ , 0d3 ×d3 ,
K
√
{
}
κ
Ω12 := − 2 diag 0(d1 +1)×(d1 +1) , Ξ , 0d3 ×d3 ,
K
{
}
1
1
Ω22 := 2 diag 0(d1 +1)×(d1 +1) , Ξ + Ψ , 0d3 ×d3 ,
1

{

K

2

Vg2 := Var {g2 (Z )} ,

{(
)
(
)}
Ξ := E η2 η2′ − Σ2 β2 β2′ η2 η2′ − Σ2 , and()Ψ := β2′ Σ2 β2 Σ2 + Σ2 β2 β2′ Σ2 .
Observe that Ω in VI can be simplified as

Ω=E

{(

Wi,j(i) ϵi,j(i) +

1
K

Σθ

)(

Wi,j(i) ϵi,j(i) +

1
K

Σθ

)′ }

where VI , VII and VIII are defined in Theorem 3.
In the cases covered by Corollary 1, the bias terms BMD as BΣ arealready of smaller order than n−1/2 or m−1/2 , and thus we obtainparametric convergence rates without bias correction.
3.3. Two-step bias-corrected estimator
If a parametric rate of convergence cannot be obtained withoutbias correction as in Corollary 1, then we need to find a way of
eliminating the second-order bias BMD , or equivalently, the effectof λi,j(i) asymptotically from (2). There are two possible strategies,
namely:Strategy 1: taking the first-order difference of (2); andStrategy 2: subtracting a consistent estimate of λi,j(i) from thedependent variable Yi .
Yatchew (1997) advocates Strategy 1 in semiparametric regression estimation, whereas Robinson (1988) and Abadie and Imbens(2011) adopt a similar strategy to Strategy 2 in semiparametric regression and ATE estimations, respectively. In our settings, we havefound that Strategy 1 has a few disadvantages. First, differencing(2) leaves β0 and γ unidentified. Second, our preliminary Monte

350

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

Carlo study suggests that MSII estimates from the differenced
regression are numerically quite unstable. For these reasons wefocus on Strategy 2.
Estimating λi,j(i) requires consistent estimates of θ and g2 (·).
For the former, it suffices to employ the MSII estimate θ̂II . For thelatter, as in Abadie and Imbens (2011),( we adopt) a nonparametricpower-series estimation. Let υ = υ1 , . . . , υd3 be a multi-index
of dimension d3 , whichof nonnegative∏∑isd3a d3 -dimensional vector
d3
υl
υ
=
integers with |υ| =l=1 zl , wherel=1 υl . Also denote z∞
zl is the lth element of z. Consider a series {υ (K)}K=1 containing
distinct vectors such that |υ (K)| is non-decreasing. Let pK (z ) =z υ(K) and pK (z ) = (p1 (z ) , . . . , pK (z ))′ . Then, a nonparametricseries estimator of the regression function g2r (z ) , r = 1, . . . , d2 ,is given byĝ2r (z ) := pK(m) (z )′

⎧
m
⎨∑
⎩

⎫−
m
( )⎬ ∑

pK(m) Zj pK(m) Zj

( )

′

⎭

j=1

pK(m) Zj X2r ,j ,

( )

j= 1

where X2r ,j is the rth element of X2j in S2 , (·)− denotes the generalized inverse, and K = K (m) signifies the dependence of K on thesample size of S2 .
The entire estimation procedure based on Strategy 2 can besummarized in the following two steps:
1. Run MSII using the original( matched sample S) to obtain the(1)initial estimate θ̂II′

= β̂II(1,0) , β̂II(1,1)′ , β̂II(1,2)′ , γ̂II(1)′ .

Assumption 7. There is a constant C such that for each multi-index
υ , the υ th partial derivative of g2 (z ) exists and its norm is bounded
by C |υ| .
The theorem below refers to the limiting distributions of θ̂II −FM
under a variety of divergence patterns of (n, m). It is worth emphasizing that asymptotic variances of the MSII-FM estimator take thesame forms as in Theorem 3 and Corollary 1, i.e., the FM procedure
asymptotically removes the second- and third-order bias terms
without inflating the variance.
Theorem 4. If Assumptions 1–7 hold, then, as n, m → ∞,

⎧√ (
) d
⎪
n
θ̂
−
θ
→ N (0, VI ) if n/m → κ ∈ (0, ∞) and d3 = 2, 3
⎪
II
−
FM
⎪ (
⎨
) d
√
,
n θ̂II −FM − θ → N (0, VII ) if n/m → 0 and d3 = 3, 4
⎪
(
) d
⎪
√
⎪
⎩ m θ̂II −FM − θ → N (0, VIII ) if n/m → ∞ and d3 = 2, 3
where VI , VII and VIII are defined in Theorem 3.
An important practical question when implementing MSII-FM
is how to choose K (m), the number of terms in the series approximation. We will address this issue in Section 4 and in Section C of
the Supplement.

2. Construct adjusted dependent variables{ + }nYi

i=1

3.4. Covariance estimation}n{,:= Yi − λ̂i,j(i)
i=1

where

λ̂i,j(i) =⎧⎨ĝ2 (Zi ) −

⎩

1 ∑
K

j∈JK (i)⎫′( )⎬ (1)ĝ2 Zjβ̂⎭ II ,2)′and ĝ2 (z ) = ĝ21 (z ) , . . . , ĝ2d2 (z ) , and rerun MSII{( +usingthe modified matched sample S + = Sn+ :=Yi , X1i ,)}n
X2j1 (i) , . . . , X2jK (i) , Zi , Zj1 (i) , . . . , ZjK (i) i=1 to obtain the final
estimator

(

n

−1 +
−1
θ̂II −FM := P̂W
R̂W := P̂W

1∑nWi,j(i) Yi+ .

i=1

We conclude this section by discussing covariance estimation,
which is essential for inference. Our focus is on the MSII estimatorof Corollary 1 and the MSII-FM estimator of Theorem 4, which are
first-order asymptotically equivalent. Because P̂W is consistent for
PW , the problem of estimating VI , VII and VIII consistently is boiled
down to proposing consistent estimators of Ω , Ω11A and Ω22 .
The consistent estimators are presented in the proposition below. Notice that the proposition is built on the assumption that θ̂IIis employed as a consistent estimator for θ ; it is easy to see that theresult equally holds after it is replaced by θ̂II −FM .
Proposition 1. Let the estimators of Ω11A , Ω22 and Ω be

The idea behind the above procedure is as follows. The initial(1)
MSII estimate θ̂II is consistent but inefficient, because the slow
convergence rate m1/d3 of the second-order bias term BMD dominates. Then, in the second step, we (asymptotically) eliminate thesource of the inferior rate by subtracting λ̂i,j(i) from the dependent
variable andθ by MSII using the bias-adjusted data
√ reestimate√
to ensure n- or m-consistency. The entire procedure is reminiscent of the fully-modified least squares estimation for cointegrating regressions by Phillips and Hansen (1990). In this sense,
we call the estimator the fully-modified MSII (MSII-FM) estimatorhereinafter.
In order to deliver convergence results for θ̂II −FM , we must
additionally impose the following regularity conditions. These areanalogous to conditions (i)-(iii) in Theorem 2 of Abadie and Imbens(2011).

nΩ̂11A =
Ω̂22 =

1∑n(

Wi,j(i) ϵ̂i,j(i) +

i=1

1
K2

Assumption 6. K (m) ≍ m for some constant(ν ∈ 0, min

{

2

,

2

4d3 + 3 4d23 − d3

})

.

K

Σ̂ θ̂II)(

Wi,j(i) ϵ̂i,j(i) +{1KΣ̂ θ̂II)′,}

diag 0(d1 +1)×(d1 +1) , Γ̂ (−1) + Γ̂ (0) + Γ̂ (1) , 0d3 ×d3 , andΩ̂ = Ω̂11A +
1

n [(
m{

)
β̂2′ ,II Σ̂2 β̂2,II W̄ W̄ ′+ 2 diag 0(d1 +1)×(d1 +1) ,
K
(
)
{
}
}]
β̂2′ ,II Σ̂2 β̂2,II V̂g2 + Γ̂ (0) − Γ̂ (−1) + Γ̂ (1) , 0d3 ×d3 ,where ϵ̂i,j(i)= Yi − Wi′,j(i) θ̂II is the MSII residual, β̂2,II is the
MSII{( estimator) of β2}, Γ̂ (ℓ) is the ℓth sample autocovariance of∆X2j ∆X2j′ /2 − Σ̂2 β̂2,II , i.e.,

Assumption 5. Z is a Cartesian product of compact intervals.
ν1Γ̂ (ℓ) =1
m−1

min{m,m+ℓ}∑(∆X2j ∆X2j′
2

j=max{2,2+ℓ}× β̂2,II β̂2′ ,II(∆X2j−ℓ ∆X2j−ℓ
′

2

− Σ̂2

)
− Σ̂2 ,)

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

⎡

⎤

1

Nonetheless we investigate them to see how the violation affects
finite-sample properties of MSII-FM. In addition, Model C is (a mirror image of) the Doppler function, which is a rapidly oscillating,spatially inhomogeneous function, as illustrated in Fig. 1 of Donoho
and Johnstone (1994). Therefore, the model may be thought of as
the most difficult case among the three. This is the model for which
we report the results here.9
Finally, Y is generated by setting all coefficients in (5) equal to

n
⎢1∑
⎥
⎡ ⎤ ⎢
X1i ⎥
⎢
⎥
1
⎢ n i=1⎥⎢⎥⎢X̄1 ⎥ ⎢ ∑
m
⎥
⎥
⎢
1
W̄ = ⎣ ⎦ = ⎢
, and
X2j ⎥
⎢
⎥
X̄2
⎢ m j=1 ⎥
⎢
⎥
Z̄
⎢ ∑
⎥
N
⎣ 1
⎦

N
V̂g2 =

1
m−1

Zi

iid

i=1

m
∑
(

X2j − X̄2)(

X2j − X̄2)′

− Σ̂2 .

j=1ppThen, under Assumptions 1–4, Ω̂11A → Ω11A , Ω̂22 → Ω22 andpΩ̂ → Ω as n, m → ∞.
4. Finite-sample performance
4.1. Monte Carlo setup
In this section we conduct Monte Carlo simulations to examine
finite-sample properties of proposed bias-corrected estimators.
The simulation study takes a unified approach in the sense thatthe same regression model is employed regardless of the numberof matching variables d3 . The model considered throughout is
Y = β0 + X1′ β1 + X2′ β2 + Z ′ γ + u,

(5)

where X1 = (X11 , X12 ) , β1 = (β11 , β12 () ∈ R , )X2 =′(X21 , X22 )′ , β2 = (β21 , β22 )′ ∈ R2 , and Z = Z1 , . . . , Zd3 , γ =()′
that two
γ1 , . . . , γd3 ∈ Rd3 for d3 = 1, 2, 3. It is assumed)}{(msamples, namely, S1 = {(Yi , X1i , Zi )}ni=1 and S2 = X2j , Zj j=1 , areonly observable. The complete sample S ∗ = {(Yi , X1i , X2i , Zi )}ni=1 is′′

2

the sample that would not be observed in practice.
∗( ∗The∗ data
) are generated in the following manner. First, Z =∗ ′Z1 , Z2 , Z3 is generated by

√

√ ⎤⎞

2/ 3

1

1/ 2 √1/ √3iidZ ∗ ∼ N ⎝⎣0⎦ , ⎣1/ 2
2/ 3⎦⎠ .
√ √ 1√⎛⎡ ⎤ ⎡
0

1
√

0

1/ 3

351

Each Zp (p = 1, 2, 3) is transformed to Zp = 4Φ Zp∗ − 2, whereΦ (·) is the cdf of N (0, 1). Observe that the Zp are mutuallycorrelated U [−2, 2] random variables. Then, for a given d3 , theZp (p ≤ d3 ) are used as matching variables.
∑d3
Second, X1 = (X11 , X12 )′ is generated by X1q =p=1 Zp +

( )

∗

iid

η1q (q = 1, 2), where η1 = (η11 , η12 )′ ∼ N ((02×)1 , I2 ). Third, X2 =∑d3
(X21 , X22 )′ is generated by X2r =p=1 g2r Zp + η2r (r = 1, 2)
iid

for some nonlinear function g2r (·), where η2 = (η21 , η22 ) ∼N (02×1 , I2 ). While g21 (z ) = z + (5/τ ) φ (z /τ ) , τ = 0.25 isemployed throughout, where φ (·) is the pdf of N (0, 1), one of thefollowing three functional forms is chosen as g22 (z ):′⎧z + (5/τ ) φ (z /τ ) , τ = 0.75 [Model A]⎪⎪⎪⎨
2 |z |[Model B]√
g22 (z ) =⎪4 |z /2| (1 − |z /2|) sin {2π (1 + ϵ) /(|z /2| + ϵ) },⎪⎪⎩ϵ = 0.05[Model C].
Both g21 (·) and Model A, which are inspired by the Monte Carlodesign of Horowitz and Spokoiny (2001), can be viewed as a linear
function with a bump. Model A is a smooth function, whereas
Models B and C have a kink at the origin. Strictly speaking, thesemodels violate the smoothness condition given in Assumption 7.

1 with u ∼ N (0, 1).
The above procedure provides us with two{(observable)}m samX2j , Zj j=1 , andples S1 = {(Yi , X1i , Zi )}ni=1 and S2 =one complete sample S ∗ . Finally, the matchedsample S ={()}nYi , X1i , X2j1 (i) , . . . , X2jK (i) , Zi , Zj1 (i) , . . . , ZjK (i) i=1 is constructed
via the NNM with respect to Z , where the NNM is based on the
Mahalanobis metric. We focus only on small numbers of matchesand examine K ∈ {1, 2, 4, 8}.10
With regards to sample sizes, for each of n ∈ {1000, 2000},
m is chosen as one of m ∈ {n/2, n, 2n} so that the values ofκ are κ = 2, 1 and 1/2, respectively. For each combination ofsample sizes (n, m) and the functional form of g22 (z ), we generate 1000 Monte Carlo samples. The following five estimators
are examined: (i) the infeasible OLS estimator using the complete
sample S ∗ [OLS*]; (ii) the MSOLS estimator using the matchedsample S and Wi,j(i) [MSOLS-A]; (iii) the MSOLS estimator using†
the matched sample S and Wi,j(i) [MSOLS-B]; (iv) the MSII(-FM)
estimator using the matched sample S and Wi,j(i) [MSII(-FM)-A];and (v) the MSII(-FM) estimator using the matched sample S and†Wi,j(i) [MSII(-FM)-B]. Second-, third- and fourth-order polynomials
are investigated in the power-series approximation for MSII-FM,
and these specifications are denoted as ‘‘2nd’’, ‘‘3rd’’ and ‘‘4th’’ inthe row ‘‘Poly.’’, respectively. Results on the initial MSII are alsoavailable as ‘‘initial ’’ for reference. Moreover, the consistent es†
timator of the second-order bias term for MSII-FM-B is λ̂i,j(i) ={}′∑λ̂i,j(i) + Zi − (1/K ) j∈JK (i) Zj γ̂II(1)′ .

We focus on finite-sample properties of estimators of β22 andγ1 . For each estimator, the following performance measures arecomputed: (i) Mean (simulation average of the parameter estimate); (ii) SD (simulation standard deviation of the parameter
estimate); (iii) RMSE (root mean-squared error of the parameter
estimate); (iv) SE (simulation average of the standard error); and(v) CR (coverage rate for the nominal 95% confidence interval).
Since MSOLS is inconsistent and limiting distributions of the initial
MSII for d3 = 2, 3 are not available, their standard errors are notwell defined. Accordingly, SE and CR are not computed for theseestimators.
4.2. Results
Simulation results are summarized in Table 1. To save space, we
present only the results from the most difficult case (Model C) for(n, m) = (1000, 1000) and (2000, 2000).
(a) For d3 = 1:. Panel (a) reports the results for a single matchingvariable. Because of conditional homoskedasticity of the error termu, OLS* is the best linear unbiased estimator. The results indicatethat it is unbiased and yields small standard deviations. However,
OLS* is an infeasible, oracle estimator. Instead, we should makea realistic comparison between MSOLS and MSII and use OLS* as
the benchmark to measure the efficiency loss when all variables
cannot be taken from a single data source.
9 Comprehensive simulation results, available in Section B of the Supplement, areeven more favorable for Models A and B.
10 In our preliminary Monte Carlo study larger values of matches (e.g., K =16, 32, 64, 128) have been also investigated. However, the results are quite poor.

352

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

Table 1Monte Carlo results for Model C.
Panel (a): d3 = 1(n, m)

Estimator(1000, 1000)OLS*

MSOLS-AMSOLS-B

MSII-A

MSII-B(2000, 2000)

OLS*

MSOLS-AMSOLS-B

MSII-A

MSII-B

β22

γ1

Mean
SD
RMSE

1.00030.02020.02020.99700.05290.0529SECR

0.020796%

0.0527
95%

K
Mean
SD
RMSE
Mean
SD
RMSE

1
0.5556
0.05120.44740.5556
0.05120.44742
0.71480.05460.29030.71480.05460.290340.83550.05820.17450.83550.05820.174580.92030.06110.10040.92030.06110.10051
1.05130.11340.12451.05130.11350.12452
1.02720.10520.10871.02710.10520.108741.01450.10190.10291.01450.10200.10308
1.00910.10080.1012
1.00920.10080.1012

K
Mean
SD
RMSE

1
1.02510.11410.11682
1.01420.0906
0.091741.01260.07740.07848
1.02210.07110.07441
0.99700.12310.1231

2
0.99800.10980.10984
0.99930.10400.104081.00130.10190.1019

SE
CR
Mean
SD
RMSE

0.104094%
1.02510.11410.11680.074089%
1.01420.0906
0.09170.063388%
1.01260.07740.07840.060990%
1.02210.07110.07440.119995%0.99700.12310.12310.106494%
0.99790.10980.10990.099493%
0.99930.10400.10400.096193%
1.00130.10190.1019SECR

0.104094%

0.074089%

0.063388%

0.060989%

0.119995%

0.106494%

0.099493%

0.0962
93%

Mean
SD
RMSE
SE
CR

0.99950.01450.01450.0147
96%

K
Mean
SD
RMSE
Mean
SD
RMSE

1
0.56020.03590.44130.56020.03590.441320.72040.03800.28210.72040.03800.282140.84060.03990.16430.84060.03990.164380.9191
0.04160.09090.9191
0.04160.090911.05020.08140.09561.05020.08140.09562
1.02630.07580.0803
1.02630.07580.080341.01370.07290.0742
1.01370.07290.07428
1.00630.0716
0.07191.00630.0716
0.0719

K
Mean
SD
RMSE
SE
CR
Mean
SD
RMSE
SE
CR

1
1.01440.07450.07580.071294%
1.01440.07450.07580.071294%

21.01000.06140.06220.051290%
1.01000.06140.06220.051290%

41.00990.05190.05280.043690%
1.00990.05190.05280.043690%

81.01350.04780.0496
0.041391%
1.01350.04780.0496
0.041391%

10.99610.08790.08800.084394%
0.99610.08790.08800.084395%

20.99750.07900.07900.075094%
0.99750.07900.07900.075094%

40.99860.07440.07450.070094%
0.99870.07440.07450.070094%

80.99850.07240.07240.067693%
0.99850.07240.07240.067693%

0.99880.03720.03720.037494%(continued on next page)†

For MSOLS, whether Wi,j(i) or Wi,j(i) is used as the regressor
has almost no difference; this reflects the fact that the( −extra)second-orderbias
inducedbyreplacing
Z
with
Z
is
O
n 1 =
i
j
i
p
(
)
( −1/2 )op n. As predicted in Theorem 1, the bias of the MSOLS estimate decreases with the number of matches K . However, it isinconsistent in that its bias does not vanish with the sample size n.
Also observe that the standard deviation of each MSOLS estimate
shrinks with n, as Theorem 1 suggests.
Now we turn to MSII. At a glance, we can find that the proposed bias-correction method works remarkably well, and that the
choice of the regressor again does not change the results. However,unlike MSOLS, increasing K has little effect at best, which suggeststhat MSII works well across small values of K . The results also
confirm consistency of MSII; as n increases, the simulation averageof each MSII estimate gets closer to the truth and its standard
deviation shrinks. In addition, SE is reasonably close to SD, whichindicates that the (properly-scaled) covariance estimator Ω̂ yields
good estimates of standard deviations of MSII. Coverage rates arealso close to the nominal level of confidence, and the single matchcase appears to have advantage from the viewpoint of coverage
accuracy.
Comparing MSII with OLS*, we have the following two findings.
First, unlike OLS*, MSII is not unbiased. However, it is nearly unbiased for large sample sizes. Second, standard deviations of thelatter are always greater than those of the former. The relative
efficiency loss can be thought of as the price to pay for identifying and estimating the regression using two samples jointly.
It is worth noting that while standard deviations of MSOLS are
greater than those of OLS*, they are smaller than those of MSII.
This(can be explained by the) fact that the asymptotic variance of√−1n θ̂OLS − QW
PW θ − BOLS2

−1
is QW
Ω11 QW−1 , which tends to be−1−1smaller (in the matrix sense) than PWΩ PW.

(b) For d3 = 2:. Next, we look into Panel (b), which presents theresults from two matching variables. Only results of MSII-FM for
K = 1 are provided, because those for K ≥ 2 are quite poor.

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

353

Table 1 (continued)
Panel (b): d3 = 2(n, m)

Estimator(1000, 1000)OLS*

MSOLS-AMSOLS-B

MSII-FM-A

MSII-FM-B(2000, 2000)

OLS*

MSOLS-AMSOLS-B

MSII-FM-A

MSII-FM-B

β22

γ1

Mean
SD
RMSE

0.99860.01650.0166

0.99770.05710.0572

SECR

0.016495%

0.0588
95%

K
Mean
SD
RMSE
Mean
SD
RMSE

1
0.47330.05280.52940.47350.05290.52922
0.63370.05710.37070.63400.05730.3705

40.7856
0.06620.22430.78580.06640.224280.94590.08470.10050.94610.08500.10061
1.05970.17670.18651.01230.17820.17862
1.02910.17250.17490.99310.17310.173241.01000.17660.17690.97950.17860.179880.97800.1967
0.19790.94570.19910.2064

Poly.
Mean
SD
RMSE

(initial)
1.17850.17680.25122nd
1.18030.17720.25283rd
1.18050.17730.25304th
1.15880.17500.2363(initial)
0.97400.21000.2116

2nd
0.97230.21230.2141

3rd
0.97250.21330.21504th
0.96670.21650.2191

SE
CR
Mean
SD
RMSE

−
−

0.168987%
1.18050.17730.25300.167990%
1.15870.17500.2363−−

1.17910.17700.25180.168887%
1.18030.17720.25280.92720.21140.22360.186992%
0.97100.21530.21730.187192%
0.97140.21600.21790.189192%
0.96540.21850.2212

SECR−−

0.168887%

0.168987%0.167990%−−0.186690%0.186890%0.1887
91%

Mean
SD
RMSE
SE
CR

0.99970.01160.01160.0116
95%

K
Mean
SD
RMSE
Mean
SD
RMSE

1
0.53650.03500.46480.53650.03510.464820.69530.03720.30700.69530.03720.306940.83740.04210.16800.83740.04210.167980.96980.05020.0586
0.96990.05030.05871
1.04290.10950.11761.01920.11050.11212
1.02000.10490.10681.00200.10550.105541.00000.10710.1071
0.98440.10770.108980.98110.11600.11750.9656
0.11710.1220

Poly.
Mean
SD
RMSE
SE
CR
Mean
SD
RMSE
SE
CR

(initial)
1.12290.09320.15432nd
1.1242
0.09330.15530.089463%
1.1242
0.09330.15530.089463%

3rd
1.1243
0.09330.15540.089463%
1.1243
0.09330.15540.089463%

4th
1.1132
0.09240.14610.089269%
1.1131
0.09240.14610.089269%(initial)
0.97870.12500.12682nd
0.97780.12540.12740.118387%
0.97740.12670.12870.118286%

3rd
0.97760.12560.12760.118387%
0.97720.12690.12890.118286%

4th
0.97520.12740.12980.119888%
0.97530.12880.13120.119687%

1.00090.04050.0406
0.041595%−−
1.12300.09330.1544−−−−
0.95480.12580.1337−−(continued on next page)†

As in the case for d3 = 1, employing Wi,j(i) or Wi,j(i) has little
effect on MSOLS or MSII-FM; although(the extra) second-order bias
generated by switching Zi to Zj(i) is Op n−1/2 , its effect appears tobe minor at best.
Even after the number of matching variables increases, the general tendency remains unchanged. Performance of MSOLS varies
with K . MSII-FM successfully corrects the bias generated by MSOLS,
at the expense of precision in estimation. Standard deviations of
MSII-FM are close to that of the initial MSII, which reflects that
the FM procedure corrects the second-order bias of MSII without
inflating the variance. However, FM works only for K = 1. The
rationale could be that FM requires both the initial MSII andsecond-order bias estimates to be of good quality. This requirement is unlikely to be satisfied with many matches, which include
poor ones and thus inevitably affect the performance of MSIIFM. In terms of the power-series approximation, results from thesecond- and third-order polynomials look similar, and those fromthe fourth-order polynomial differ slightly. Coverage accuracy inestimates of β22 may be a concern. However, based on the results

for larger samples (reported in the Supplement), it seems that theunder-coverage is due to finite-sample bias of MSII-FM.
(c) For d3 = 3:. In Panel (c), only results of MSII-FM for K = 1
are provided again in view of quality. An apparent difference isthat once the number of matching variables increases to three,†results from using Wi,j(i) or Wi,j(i) differ substantially for each of
MSOLS and MSII-FM. Observe that MSII-FM using Wi,j(i) exhibits
much better finite-sample properties. In contrast, MSII-FM based†on Wi,j(i) generates considerable biases in estimates of γ . The extra
second-order( bias) when Zj(i) is used in place of Zi becomes as
slow as Op n−1/3 , and its adverse effect is no longer negligible
in finite samples. Coverage rates of MSII-FM are improved fromthose for d3 = 2. In terms of the series approximation, resultsfrom the second-and third-order polynomials are again similar.
However, those from the fourth-order polynomial look inferior
in the presence of non-smoothness in g22 (·), in particular, forModel B.

354

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

Table 1 (continued)
Panel (c): d3 = 3(n, m)

Estimator(1000, 1000)OLS*

MSOLS-AMSOLS-B

MSII-FM-A

MSII-FM-B(2000, 2000)

OLS*

MSOLS-AMSOLS-B

MSII-FM-A

MSII-FM-B

β22

γ1

Mean
SD
RMSE

0.99940.01350.01350.99970.05800.0580SECR

0.013996%

0.0585
96%

K
Mean
SD
RMSE
Mean
SD
RMSE

1
0.21930.07480.78430.22050.07550.78322
0.36870.08870.63750.37030.08950.63604
0.57580.11630.43980.57880.11680.437180.89420.15280.18590.89940.15420.1842

1
1.14980.31030.34450.64390.31760.47722
1.06580.30500.31210.68350.3146
0.4463

4
0.99780.32460.3246
0.67750.3406
0.469180.93330.36010.36630.62990.38370.5332

Poly.
Mean
SD
RMSE

(initial)
1.11510.40640.42242nd
1.08890.40090.4106

3rd
1.09010.40050.41054th
1.06510.39530.4007(initial)
0.97630.36980.3705

2nd0.95500.37510.37783rd
0.95340.37700.37994th
0.94040.37120.3760

SE
CR
Mean
SD
RMSE

−
−

0.372692%
1.09030.40090.41100.366991%
1.06490.39540.4007−−1.12170.40990.42750.371892%
1.08900.40120.41100.47090.37770.6501

0.328885%0.83580.42100.45190.330485%
0.83180.42490.45700.324586%0.82000.41360.4511SECR−−

0.372292%

0.373092%

0.367391%−−

0.327377%

0.329077%0.3200
76%

Mean
SD
RMSE
SE
CR

1.00020.00960.00960.0099
96%

K
Mean
SD
RMSE
Mean
SD
RMSE

1
0.29940.04540.70210.30030.04590.70122
0.46530.05410.5374
0.46640.0546
0.5364

40.66320.06570.3432
0.66480.06640.341780.91490.08770.12220.91750.08860.12111
1.10370.20070.22590.75340.20330.31952
1.04920.19040.19670.79110.19310.2845

4
0.99100.1946
0.19480.78040.19770.295580.93470.21690.22650.73660.22350.3454

Poly.
Mean
SD
RMSE
SE
CR
Mean
SD
RMSE
SE
CR

(initial)
1.05760.18260.19152nd
1.04770.18160.18770.180097%
1.04770.18170.18790.180196%

3rd
1.04810.18180.18810.180297%
1.04810.18190.18820.180296%

4th
1.01910.17870.17970.177196%
1.01890.17880.17980.177296%(initial)
0.98000.22770.2286

2nd
0.96670.23050.23280.198590%0.90940.2526
0.26830.199684%

3rd
0.96630.23070.2332
0.198990%0.90940.25300.26870.200085%

4th
0.96170.22390.22710.196191%
0.90320.24380.26230.195585%−−
1.06080.18400.1938−−0.9991
0.04190.04190.041594%−−
0.63000.23280.4371−−

Note: Mean = simulation average of the parameter estimate; SD = simulation standard deviation of the parameter estimate; RMSE = root mean-squared error of theparameter estimate; SE = simulation average of the standard error; and CR = coverage rate for the nominal 95% confidence interval.

(d) Summary:. Simulation results confirm that the bias-corrected
estimation proposed in this paper works reasonably well. Simulation averages of MSII(-FM) for d3 = 1 (d3 = 2, 3) tend to becloser to the truths as n increases, even in the most difficult case.
Judging from the Monte Carlo evidence, we recommend setting
K = 1, employing Wi,j(i) as the regressor, and applying the secondor third-order polynomials for the series approximation in MSIIFM. It follows that making MSOLS consistent by use of K -NN
method (i.e., by letting K diverge at a slower rate than n and m)
does not appear to be a solution in the setting of matched sample
estimation. Rather, it looks promising to pursue the strategy ofconstructing a matched sample based on a single match and thencorrecting the non-negligible bias of the estimate analytically.
5. An empirical application: Returns to schooling
We now apply our proposed estimation methods to a versionof Mincer’s (1974) wage regression. As argued in Card (1995),

the estimation result may suffer from the ‘‘ability bias’’ unless it
includes a variable representing ability as a regressor. Therefore,
we consider the following wage regressionlog (w age) = β0 + β1 educ + β2 exper + β3 exper 2 + β4 abil

+ β5 feduc + β6 meduc + β7 black + β8 smsa + β9 south + u,(6)

where educ is years of education, exper is work experience, abil
is an ability measure, feduc and meduc are years of father’s andmother’s education, and black, smsa and south are indicator variables that take one if the individual is black, lives in the urban areaand south, respectively.
We estimate regression (6) using three data sets, namely, thoseused in Card (1995), Blackburn and Neumark (1992), and Heckman
et al. (2000). The data sets are available under the names ‘‘card’’,‘‘wage2’’ and ‘‘htv ’’, respectively, as supplemental materials for
Wooldridge (2013). Each of the three data sets is drawn fromthe National Longitudinal Survey (NLS) and contains some ability

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

355

Table 2Estimation results of wage regressions with ability measures.
Dependent variable: log (w age)
Regressors
educ
exper
exper 2
abil
feduc
meduc
black
smsa
south
intercept
abil?
Matching?
(n, m)

(1)

( 2)

(3)

(4)

OLS*

MSOLS

MSII-FM

MSOLS

0.0612(0.0054)0.0787(0.0084)−0.0022(0.0004)
0.0056(0.0013)−0.0018(0.0031)0.0071(0.0037)−0.1321(0.0258)
0.1517(0.0179)−0.1111(0.0178)
4.6861(0.0841)

0.0736(0.0050)0.0875(0.0082)−0.0023(0.0004)−0.0007(0.0014)−0.0006(0.0031)0.0080(0.0037)−0.1664(0.0259)
0.1602(0.0183)−0.1126(0.0179)
4.6491(0.0861)0.08470.0876(0.0083)−0.0022(0.0004)(0.0081)−0.0023(0.0004)

0.0016(0.0046)−0.0007(0.0031)0.0073(0.0039)−0.1559(0.0331)0.1576(0.0186)−0.1125(0.0178)4.6540(0.1107)0.0006(0.0049)−0.0007(0.0032)
0.0079(0.0037)−0.1630(0.0249)
0.1595(0.0181)−0.1125(0.0180)
4.6425(0.0849)

0.0693(0.0165)0.0876(0.0082)−0.0023(0.0004)
0.0070(0.0356)−0.0010(0.0038)
0.0072(0.0041)−0.1607(0.0283)
0.1612(0.0198)−0.1104(0.0216)
4.6818(0.1945)kww
No(2191, −)kww
Yes(2191, 466)

kww
Yes
(2191, 466)g
Yes(2191, 589)g
Yes
(2191, 589)

0.06900.0724(0.0074)(0.0050)(5)MSII-FM

Note: Numbers in parentheses are standard errors. White’s (1980) heteroskedasticity-robust standard errors are calcu−1
lated for OLS*, whereas ‘standard errors’ for MSOLS are square-roots of diagonal elements of Q̂WΩ̂11 Q̂W−1 /n.

measure; to be precise, while both card and wage2 include scoresof IQ and Knowledge of the World of Work (kww ) tests, htv has the‘‘g’’ measure constructed from 10 component tests of the Armed
Services Vocational Aptitude Battery.
We conduct two exercises that address the following questions:(Q1) How would the estimation result change if kww in card
were missing and instead taken from wage2?
(Q2) What would happen if kww in card were replaced by gfrom htv?
For these exercises, the OLS result using 2191 male observations
in card with kww chosen as abil can be viewed as the benchmark
result from the infeasible OLS*. Because each of Q1 and Q2 requires
a matched sample, we regard card as S1 and wage2 or htv as S2 .
The NNM is made in the following manner. When wage2 is employed as S2 , (educ , feduc , meduc , black, smsa, south) are chosen
as matching variables, where the first three variables are treatedas continuous. On the other hand, htv contains only white-maleobservations. Accordingly, when using it as S2 , we choose five
matching variables excluding black. Not surprisingly, there are several ties of the matching variables in S2 . Then, we take an averageof kww or g within ties and assign the average as the unique valueof the ability measure to each combination of matching variables.
As a consequence, 466 and 589 distinct combinations of matchingvariables remain in male samples of wage2 and htv, respectively.
In both cases, the NNM is based on the Mahalanobis metric, andwe set the number of matches K = 1 (single match) based on thesimulation results.
Given the matched sample, we estimate (6) by MSOLS andMSII-FM. Specifically, MSOLS-A and MSII-FM-A (i.e., estimatorswith Wi,j(i) used as the regressor) are chosen, and the third-orderpolynomial is applied for the power-series approximation of MSIIFM, again based on the simulation results; estimation results fromsecond- and fourth-order polynomials are qualitatively similar.
Table 2 presents estimation results and standard errors (inparentheses). White’s (1980) heteroskedasticity-robust standard

errors are computed for OLS*, whereas ‘standard errors’ for MSOLS−1−1
are square-roots()of diagonal elements of Q̂W Ω̂11 Q̂W /n :=−1
Q̂WΩ̂11A + Ω̂11B Q̂W−1 /n, wheren [()β̂2′ ,II Σ̂2 β̂2,II W̄ W̄ ′m{
1
+ 2 diag 0(d1 +1)×(d1 +1) ,
}
}]
{
(K
)
β̂2′ ,II Σ̂2 β̂2,II V̂g2 + 2 Γ̂ (−1) + Γ̂ (1) , 0d3 ×d3Ω̂11B =for (n, m) given in the corresponding column of Table 2. The lattershould be interpreted with caution; because θ̂OLS is inconsistent(and even its convergence rate is slower than the parametric one),the numbers merely indicate measures of dispersion at the same
scale as other estimates and are not intended for inference.
The benchmark OLS* result using card is provided in the firstcolumn. Signs of the coefficient estimates on educ, exper, exper 2 ,and abil (= kww) are as expected, and they are significant at the
5% level. To answer Q1, we run MSOLS and MSII-FM using thematched sample with wage2. The results are reported in columns
2 and 3. Signs of the coefficient estimates by MSII-FM are thesame as those by OLS*. On the other hand, MSOLS overestimates
returns to schooling due to failure to correct for matching results.
It also yields a negative estimate of the ability effect, whereas theone from MSII-FM is positive (but insignificant due to the large
standard error).
Furthermore, to answer Q2, we replace the ability measure withg by constructing the matched sample with htv. Results fromMSOLS and MSII-FM using this sample are presented in columns 4
and 5. There is still the tendency that MSII-FM estimates are closer
to those of OLS*. MSOLS again tends to inflate returns to schooling.
The estimated ability effect turns positive, but its magnitude ismuch smaller than the one from MSII-FM.
6. Conclusion
Regression estimation using samples constructed via the NNM
from two sources is not uncommon in applied economics. This

356

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

paper has demonstrated that such OLS estimators are generally
inconsistent, and thus an appropriate bias correction is required.
It has also been shown that the convergence rate to the probability
limit of the OLS depends on the number of matching variables andthe divergence pattern of two sample sizes.
Two versions of bias-corrected estimators have been proposed,and each can be interpreted as a variant of indirect inference
estimators. The MSII estimator attains the parametric convergence
rate for the cases with at most two matching variables, whereas the
MSII-FM estimator achieves the parametric convergence rate when
the number of matching variables does not exceed four. Monte
Carlo results suggest that a small number of matches work wellin practice, and in particular, we should consider the single match
when the number of matching variables is two or three.
The paper aims at providing corrections for an establishedpractice, which is to run (parametric) OLS ignoring imputation. In
particular, our proposal for MSII-FM is based on a nonparametricseries estimation of g2 (Z ) = E(X2 |Z ). The nonparametric estimator
is employed only when the curse of dimensionality in matchingvariables prevents MSII from attaining parametric convergence.
Alternatively, it is possible to use a nonparametric estimate ofg2 (Z ) as a (generated) regressor in place of X2 in regression (1)from the beginning. As illustrated in Section 2.4, this would resultin a partially linear semiparametric model with a measurement
error problem. There could be several different (nonparametric)
estimators available for the model of this class. However, such
estimators are not as widely used in practice as the OLS, and hence
we leave the development of the estimators for future work.
Several other extensions would be fruitful. First, we may adopt
propensity score matching as a means of dimension reduction
using multiple matching variables. This would involve using theobservable variables to estimate a selection model for observations that are imputed, and obtaining the (imputation) propensity
score. In a closely related paper, Abadie and Imbens (2016) deliver
asymptotic properties of the matching estimators of ATE using an
estimated propensity score as a plug-in. It may be worth pursuinga similar idea for matched-sample regression estimation.
Second, combining our matched-sample estimation theory
with IV/GMM estimation would be also of interest in the presenceof endogeneity in regressors. This is particularly relevant to empirical studies using earnings data, which are thought to include
measurement errors and imputation biases.
Finally, the estimation theory may be extended to kernel estimation of varying coefficient models using matched samples. It is
not difficult to see that kernel estimators of the varying coefficientsare also inconsistent, and appropriate bias-correction methods
similar to those proposed in this paper are worth investigating.

Lemma A.1 (Abadie and Imbens, 2006, Lemma 2). If Assumptions
1–2 hold, then all the moments of m1/d3 ∥Uk ∥ are uniformly boundedin m and z ∈ Z.
A.2. Proof of Theorem 1
The proof requires the following lemma.
Lemma
If Assumptions 1–3 hold, then BRW 2 = O m−1/d3 +( −1/2A.2.
)
op n
.

(

)

A.2.1. Proof of Lemma A.2
(
) {
(
)}
Consider the identity BRW 2 ≡ E BRW 2 + BRW 2 − E BRW 2 . It
follows from Assumption 1, Lipschitz continuity
Cauchy– ( 1of/dg2 , the)
E m 3 BR 2  < ∞,Schwarz inequality,
and
Lemma
A.1
that
W(
(
)
( −1/d )3 . Similarly, we have Var n1/2 m1/d3
and thus) E BRW 2 = O m
BR 2  < ∞ so thatW
BRW 2 = O m−1/d3 + Op n−1/2 m−1/d3 = O m−1/d3 + op n−1/2 .

(

)

(

)

(

)

(

)

■

A.2.2. Proof of Theorem 1
It is easy to see from (2) that R̂W := Q̂W θ + BRW 1 + BRW 2 + ERW ,where
BRW 1 = E Wi,j(i) ϵi,j(i) ,

(

)

n

BRW 2 =E RW =

1∑n

Wi,j(i) λi,j(i) , and

i=1
n

(
)}
1 ∑{
Wi,j(i) ϵi,j(i) − E Wi,j(i) ϵi,j(i) .
n
i=1

It follows that θ̂OLS := θ + BOLS1 + BOLS2 + EOLS , where BOLS1 =−1−1−1
BRW 1 , BOLS2 = Q̂W
BRW 2 and EOLS = Q̂W
ERW correspond toQ̂W
the first-order (or leading) bias, the second-order bias due to thematching discrepancy and the weighted average of errors, respectively.
)
(
′
=
with
evaluating
BOLS1 . (First note) that E X1i η2i)}({ We begin
′
=
1
/
K
Σ
,
and
E g1 (Z ) η2′ + E η1 η2′ = 0d1 ×d2 , E X2j(i) η2j
(
)
2
(i)
that the ith and jk (i)th observations are independent. Then,⎡0(d1 +1)×1⎤

BRW 1 = ⎣− (1/K ) Σ2 β2 ⎦0d3 ×1=−1K1diag 0(d1 +1)×(d1 +1) , Σ2 , 0d3 ×d3 θ{}:= − Σ θ.
K

Appendix A. Technical proofs
A.1. A useful lemma
Before proceeding, we present a lemma about the error bounds
from NNM, which is repeatedly applied in the technical proofs
below. To do so, we provide the formal definition of the matching
discrepancy from Abadie and Imbens (2006).
Let z ∈ Z be a fixed value of the matching variable Z , where,
in practice, z is one of {Zi }ni=1 in S1 . Then, the kth closest matching
discrepancy Uk = Uk (z ) , k = 1, . . . , K is defined{as}Uk := Zjk (z ) − zm
if Zjk (z ) is the kth closest match to z among all Zj j=1 in S2 . The
following lemma states uniform moment bounds of the matching
discrepancy.

Because Q̂W = QW + Op n−1/2 , we obtain()−1
BOLS1 = − (1/K() QWΣ)θ + Op( n−1/2) . Next, Lemma A.2 implies()
that BOLS2 = O m−1/d3 + op n−1/2 . Finally, ERW = Op n−1/2( −1/2 )−1by( CLT, and=) Op n. Therefore, θ̂OLS = QW
PW θ +) thus(EOLS
−1/d3
−1/2
by denoting PW := QW − (1/K ) Σ . ■O m+ Op n()

A.3. Proof of Theorem 2
The proof requires the following lemma.
Lemma A.3. If Assumptions 1–3(hold, then) Σ̂2 and BΣ(2 admit
)the
−1/2
−2/d3
expansions
Σ̂
−
B
=
Σ
+
O
m
and
B
=
O
m
+
2
Σ
2
p
Σ
2
2
(
)
op m−1/2 .

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

A.3.1. Proof of Lemma A.3
Observe that

Σ̂2 =

=

1

m
∑

2 (m − 1)

j=2

1

m
∑

2 (m − 1)j=2+∆X2(j) ∆X2′ (j)2 (m − 1)n/m→κ∆η2(j) ∆η2′ (j) +

m
∑

1

1

m
∑

2 (m − 1)j=2∆g2(j) ∆η2′ (j) +j=2∆η2(j) ∆g2′ (j)1m∑

2 (m − 1)j=2∆g2(j) ∆g2′ (j)

2
∑m 

 /{2 (m − 1)}. Bein ∥E (S4 )∥ is determined by E
j=2 ∆Z(j){ }mcause the re-ordered sample Z(j) j=1 is constructed by the nearestneighbor algorithm in accordance with Yatchew (1997), we may
apply the moment {approximation in} Theorem 5.4 of Evans et al.

2
(
)
∑m 

 /{2 (m − 1)} = O m−2/d3(2002) to obtain E
j=2 ∆Z(j)()( 1/2+2/d)
3 ∥S ∥
so that E (S4 ) = O m−2/d3 . Alsom< ∞) by
4
(
( Var
)
the
and( thus )S4 = O m−2/d3 + Op m−(1/2+2/d3 ) =( theorem,)
O m−2/d3 + op m−1/2 . Third, for S2 and S3 , we have E (S2 ) =E (S3() = 0d2 ×d2 .) Because,( by a similar) argument to above,
Var (m1/2+1/d3 ∥S3)∥ , Var( m1/2+)1/d3 ∥S3 ∥ < ∞, each of S2 and S3
is Op m−(1/2+1/d3 ) = op m−1/2 . Finally, rewriting S4 as BΣ2 yields

R̂W =

Q̂W −

K

−1
VII := PW
lim VarFinally, if n/m

K

m→∞

K

Σ̂ − Σ θ + ERW .

K

diag 0(d1 +1)×(d1 +1) , BΣ2 , 0d3 ×d3 θ{√ (m θ̂II − θ − BMD − BΣ−1mEΣ2 PW
.

)

)

d

→

■

A.6. Proof of Theorem 4
The proof requires the following lemma.
Lemma A.4. If Assumptions 1–7 hold, then

⏐
⏐

⏐
⏐

max ⏐λ̂i,j(i) − λi,j(i) ⏐

1≤i≤n

=

⎧ ( −1/2 )⎨ op nif n/m → κ and d3 = 2, 3 or

−1/2

if n/m → ∞ and d3 = 2, 3)Σ̂ − Σ θ{()}
1= diag 0(d1 +1)×(d1 +1) , Σ̂2 − BΣ2 − Σ2 , 0d3 ×d3 θ
1∞, thenBecause the corollary is a straightforward outcome of Theorem 3, its proof is omitted. ■

R1i =

)

+

)

A.5. Proof of Corollary 1

R2i =It also follows from Lemma A.3 that
1 (

(√

−1
nERW PW
.

(

op m

)

if n/m → 0 and d3 = 3, 4.

A.6.1. Proof of Lemma A.4
It is easy to see that λ̂i,j(i) := R1i + R2i + R3i + λi,j(i) , whereΣ θ + BRW 2 + ERW= P̂W θ + BRW 2 +→−1
VIII := PW
lim Var

)

1 (

(√

n,m→∞
n/m→0

⎩

■

A.3.2. Proof of Theorem 2
By the proof of Theorem 1,
1

N (0, VII ), where

N (0, VIII ), where

As in the proof of Lemma A.2, consider the identity Si ≡ E (Si ) +{Si − E (Si()} for i )= 1, 2, 3, 4. First, S1 = E (S1 ) + {S1 − E (S1 )} =
Σ2 + Op m−1/2 . Second, for S4 , it follows from Lipschitz conof magnitude
tinuity of g2 (·) in Assumption{3(ii) that the order
}

(

) d
√ (
n θ̂II − θ − BMD − BΣ→

Alternatively, if n/m → 0, then

= S1 + S2 + S3 + S4 (say).

the stated result.

357

) d
√ (
n θ̂II − θ − BMD − BΣ → N (0, VI ), where
{√ (
)} −1
−1
VI := PW
lim Varn ERW + EΣ2 PW
.
n,m→∞

}

K
(
) { (
)
(
)}
= Op m−1/2 + O m−2/d3 + op m−1/2 ,where the first term on the right-hand side is denoted as EΣ2hereinafter. Therefore,−1θ̂II = P̂W
R̂W
−1
−1
= θ + BMD + BΣ + P̂W
ERW + P̂WEΣ{ ( −1/d )( −1/2 )} {2 ( −2/d )( −1/2 )}
3 + o
3 + o
=θ+ O m
n
+
O m
p
p m
(
)
(
) p
+ Op n−1/2 + Op m−1/2 → θ. ■(A.1)

A.4. Proof of Theorem 3
To save space, we only demonstrate asymptotic normality results depending on the divergence pattern of (n, m); see Section Aof the Supplement for detailed derivation of asymptotic variances
VI , VII and VIII . It follows from (A.1) that if n/m → κ , thenR3i =1 ∑ [{
K

}

{

( )

j∈JK (i)1 ∑ [{
K)( )}]′ ( (1)β̂II ,2 − β2 ,ĝ2 (Zi ) − g2 (Zi ) − ĝ2 Zj − g2 Zj( )}]′ĝ2 (Zi ) − g2 (Zi ) − ĝ2 Zj − g2 Zj

}

{

( )

β2 , andj∈JK (i))( )}′ ( (1)1 ∑ {g2 (Zi ) − g2 Zj
β̂II ,2 − β2 .
K
j∈JK (i)Hence, the proof is boiled down to demonstrating that( −each
) of
1/2
max
or
( 1≤−i1≤/n2 )|Rℓi | , ℓ = 1, 2, 3 is bounded by either op n
op m
, depending on the divergence pattern of (n, m) and d3 .
We first work on R3i . To derive the bounds for R1i and R3i , we
may apply (A.1) to obtain

⎧ (
)
(
)
⎪
O n−1/d3 + Op n−1/2⎪⎪⎪⎪
if n/m → κ and d3 = 2, 3⎪⎪()
⎨ ( −1/d3 )O m+ Op n−1/2θ̂II(1) = θ +.
⎪ if n/m → 0 and d3 = 3, 4
⎪
(
)
(
)
⎪
⎪
⎪
O m−1/d3 + Op m−1/2⎪⎪⎩
if n/m → ∞ and d3 = 2, 3
It follows from Lemma A.1 and Lipschitz continuity of g2 (·) thatmax1≤i≤n |R3i | is bounded by⎧ ( −1/d ) { ( −1/d )( −1/2 )}
3
3 + O
Op m
O n
⎪
p n
⎪
(
)
⎪
⎪ = o n−1/2 if n/m → κ and d = 2, 3
⎪
p
3
⎪
⎪
(
)}
⎨ ( −1/d3 ) { ( −1/d3 )Op m
O m+ Op n−1/2( −1/2 ).
⎪= op nif n/m → 0 and d3 = 3, 4
⎪
⎪
(
)
{
(
)
(
)}
⎪
⎪
Op m−1/d3 O m−1/d3 + Op m−1/2
⎪
⎪
(
)
⎩
= op m−1/2
if n/m → ∞ and d3 = 2, 3

358

M. Hirukawa, A. Prokhorov / Journal of Econometrics 203 (2018) 344–358

The remaining task is to demonstrate that for k = 1, . . . , K ,

{
} { (
)
(
)}
max  ĝ2 (Zi ) − g2 (Zi ) − ĝ2 Zj (i) − g2 Zj (i) 
k

1≤i≤n

=

k

⎧ ( −1/2 )⎨ op nif n/m → κ and d3 = 2, 3 orop m−1/2

if n/m → ∞ and d3 = 2, 3⎩(

if n/m → 0 and d3 = 3, 4

)

.

(A.2)

However, Lemma A.2 of Abadie and Imbens (2011) holds under
Assumptions 1–7. Therefore,max ⏐ ĝ2r (Zi ) − ĝ2r Zjk (i)

⏐{

(

1≤i≤n

)}

(
)}⏐
{
− g2r (Zi ) − g2r Zjk (i) ⏐)(= op m−1/2 , r = 1, . . . , d2 ,and thus (A.2) immediately follows. Then, each of max( 1−≤1i≤/2n)|R1i |
|
|
is
also
boundedbyeither
o
n
or
R
and
max
p
2i
1
≤
i
≤
n
(
)
op m−1/2 . This completes the proof. ■A.6.2. Proof of Theorem 4
To save space, we focus on the case in which n/m → κ ∈ (0, ∞)(1)(1)
and d3 = 2, 3. Observe that θ̂II −FM :=( θ̂II − B̂MD := θ̂II) −√d(1)n θ̂II − θ − BMD − BΣ →


(
)


N (0, VI ). Lemma A.4 implies that BMD − B̂MD  = op n−1/2 . Then,)}√ {(n BMD − B̂MD + BΣ = op (1), and thus the result immedi-

(1/n)

∑n

λ̂

i=1 Wi,j(i) i,j(i) . By Theorem 3,

ately follows. ■A.7. Proof of Proposition 1
The proof is obvious in light of Section A of the Supplement andthus omitted. ■Appendix B. Supplementary data
Supplementary material related to this article can be foundonline at https://doi.org/10.1016/j.jeconom.2017.07.006.
References
Abadie, A., Imbens, G.W., 2006. Large sample properties of matching estimators foraverage treatment effects. Econometrica 74, 235–267.
Abadie, A., Imbens, G.W., 2011. Bias-corrected matching estimators for averagetreatment effects. J. Bus. Econom. Statist. 29, 1–11.
Abadie, A., Imbens, G.W., 2012. A martingale representation for matching estimators. J. Amer. Statist. Assoc. 107, 833–843.
Abadie, A., Imbens, G.W., 2016. Matching on the estimated propensity score.
Econometrica 84, 781–807.
Angrist, J.D., Krueger, A.B., 1992. The effect of age at school entry on educational
attainment: An application of instrumental variables with moments from two
samples. J. Amer. Statist. Assoc. 87, 328–336.
Angrist, J.D., Krueger, A.B., 1995. Split-sample instrumental variables estimates of
the return to schooling. J. Bus. Econom. Statist. 13, 225–235.
Arellano, M., Meghir, C., 1992. Female labour supply and on the job search: An
Empirical model estimated using complementary data sets. Rev. Econom. Stud.
59, 537–559.
Björklund, J., Jäntti, M., 1997. Intergenerational income mobility in Sweden compared to the United States. Am. Econ. Rev. 87, 1009–1018.
Blackburn, M., Neumark, D., 1992. Unobserved ability, efficiency wages, and interindustry wage differentials. Quart. J. Econ. 1077, 1421–1436.
Bollinger, C.R., Hirsch, B.T., 2006. Match bias from earnings imputation in thecurrent population survey: The case of imperfect matching. J. Labor Econ. 24,
483–519.
Borjas, G.J., 2004. Food insecurity and public assistance. J. Publ. Econ. 88,
1421–1443.
Bostic, R., Gabriel, S., Painter, G., 2009. Housing wealth, financial wealth, andconsumption: New evidence from micro data. Reg. Sci. Urban Econ. 39, 79–89.
Bover, O., 2005. Wealth Effects on Consumption: Microeconometric Estimates fromthe Spanish Survey of Household Finances. Documentos de Trabajo No.0522,
Banco de España.

Busso, M., DiNardo, J., McCrary, J., 2014. New evidence on the finite sample properties of propensity score reweighting and matching estimators. Rev. Econ. Stat.
96, 885–897.
Card, D., 1995. Using geographic variation in college proximity to estimate thereturn to schooling. In: Christophides, L.N., Grant, E.K., Swidinsky, R. (Eds.),
Aspects of Labour Market Behavior: Essays in Honour of John Vanderkamp.
University of Toronto Press, Toronto, pp. 201–222.
Carrasco, M., Florens, J.-P., 2002. Simulation-based method of moments and efficiency. J. Bus. Econom. Statist. 20, 482–492.
Chen, J., Shao, H., 2001. Jackknife variance estimation for nearest-neighbor imputation. J. Amer. Statist. Assoc. 96, 260–269.
Currie, J., Yelowitz, A., 2000. Are public housing projects good for kids?. J. Publ. Econ.
75, 99–124.
Dee, T.S., Evans, W.N., 2003. Teen drinking and educational attainment: Evidence
from two-sample instrumental variables estimates. J. Labor Econ. 21, 178–209.
Donoho, D.L., Johnstone, I.M., 1994. Ideal spacial adaptation by wavelet shrinkage.
Biometrika 81, 425–455.
Evans, D., Jones, A.J., Schmidt, W.M., 2002. Asymptotic moments of near-neighbour
distance distributions. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 458,
2839–2849.
Fujii, T., 2008. Two-Sample Estimation of Poverty Rates for Disabled People: An Application to Tanzania. Singapore Management University Economics & Statistics
Working Paper No.02-2008.
Gouriéroux, C., Monfort, A., Renault, E., 1993. Indirect inference. J. Appl. Econometrics 8, S85–S118.
Heckman, J., Tobias, J.L., Vytlacil, E., 2000. Simple Estimators for Treatment Parameters in a Latent Variable Framework with an Application to Estimating the
Return to Schooling. NBER Working Paper No.7950.
Hellerstein, J.K., Imbens, G.W., 1999. Imposing moment restrictions from auxiliary
data by weighting. Rev. Econ. Stat. 81, 1–14.
Hirsch, B.T., Schumacher, E.J., 2004. Match bias in wage gap estimates due toearnings imputation. J. Labor Econ. 22, 689–722.
Horowitz, J.L., Spokoiny, V.G., 2001. An adaptive, rate-optimal test of a parametricmean-regression model against a nonparametric alternative. Econometrica 69,
599–631.
Imbens, G.W., Lancaster, T., 1994. Combining micro and macro data in microeconometric models. Rev. Econom. Stud. 61, 655–680.
Inoue, A., Solon, G., 2010. Two-sample instrumental variables estimators. Rev. Econ.
Stat. 92, 557–561.
Little, R.J.A., Rubin, D.B., 2002. Statistical Analysis with Missing Data, second ed..
John Wiley & Sons, New York.
Lusardi, A., 1996. Permanent income, current income, and consumption: Evidence
from two panel data sets. J. Bus. Econom. Statist. 14, 81–90.
Mincer, J.A., 1974. Schooling, Experience and Earnings. National Bureau of Economic
Research, New York.
Murtazashvili, I., Liu, D., Prokhorov, A., 2015. Two-sample nonparametric estimation of intergenerational income mobility in the United States and Sweden. Can.
J. Econ. 48, 1733–1761.
Pagan, A., 1984. Econometric issues in the analysis of regressions with generated
regressors. Internat. Econom. Rev. 25, 221–247.
Phillips, P.C.B., Hansen, B.E., 1990. Statistical inference in instrumental variables
regression with I(1) processes. Rev. Econom. Stud. 57, 99–125.
Prokhorov, A., Schmidt, P., 2009. GMM redundancy results for general missing data
problems. J. Econometrics 151, 47–55.
Rice, J., 1984. Bandwidth choice for nonparametric regression. Ann. Statist. 12,
1215–1230.
Ridder, G., Moffitt, R., 2007. The econometrics of data combination. In: Heckman,
J.J., Leamer, E.E. (Eds.), Handbook of Econometrics, Vol. 6, Part B. Elsevier,
Amsterdam, pp. 5469–5547 (Chapter 75).
Robinson, P.M., 1988. Root-N-consistent semiparametric regression. Econometrica
56, 931–954.
Shao, J., Wang, H., 2008. Confidence intervals based on survey data with nearest
neighbor imputation. Statist. Sinica 18, 281–297.
Smith, J.A., Todd, P.E., 2005. Does matching overcome LaLonde’s critique of nonexperimental estimators. J. Econometrics 125, 305–353.
Smith Jr., A.A., 1993. Estimating nonlinear time-series models using simulated
vector autoregressions. J. Appl. Econometrics 8, S63– S84.
von Neumann, J., 1941. Distribution of the ratio of the mean square successive
difference to the variance. Ann. Math. Stat. 12, 367–395.
White, H., 1980. A heteroskedasticity-consistent covariance matrix estimator and a
direct test for heteroskedasticity. Econometrica 48, 817–838.
Wooldridge, J.M., 2013. Introductory Econometrics: A Modern Approach, fifth ed..
South-Western Cengage Learning, Mason, OH.
Yatchew, A., 1997. An elementary estimator of the partial linear model. Econom.
Lett. 57, 135–143.