The Stata Journal (2021)
21, Number 1, pp. 123–140

DOI: 10.1177/1536867X211000008

msreg: A command for consistent estimation of
linear regression models using matched data
Masayuki Hirukawa
Ryukoku University
Kyoto, Japan
hirukawa@econ.ryukoku.ac.jp

Di LiuStataCorp
College Station, TX
dliu@stata.com

Artem Prokhorov
University of Sydney Business School
Sydney, Australia, and
CEBA, St. Petersburg State University
St. Petersburg, Russia, and
CIREQ, University of Montreal
Montreal, Canada
artem.prokhorov@sydney.edu.au
Abstract. Economists often use matched samples, especially when dealing withearning data where some observations are missing in one sample and need to
be imputed from another sample. Hirukawa and Prokhorov (2018, Journal of
Econometrics 203: 344–358) show that the ordinary least-squares estimator using
matched samples is inconsistent and propose two consistent estimators. We describe a new command, msreg, that implements these two consistent estimators
based on two samples. The estimators attain the parametric convergence rate if
the number of continuous matching variables is no greater than four.
Keywords: st0630, msreg, bias correction, linear regression, matching estimation

1

Introduction

Matching-based imputation is common in economic datasets. For example, the U.S.
Census uses a practice known as “hot-deck imputation”, which is implemented when
census reports for nonresponders the values of important variables such as earnings
and income are borrowed from responders with a few similar characteristics. In somesurveys, the share of such imputed responses reaches 30%.
Hirukawa and Prokhorov (2018) were concerned with this widely used but often ignored practice. The concern was that users of such data in applied econometrics areoften unaware that these are imputed, rather than actual, observations and that the resulting matching discrepancy leads to nonnegligible biases in the ordinary least-squares(OLS) estimator. They list many other settings, usually involving more than one dataset,
where matching is unavoidable and needs to be accounted for.
The goal of this article is to facilitate the use of the consistent estimation approaches
proposed by Hirukawa and Prokhorov (2018) in their numerous applications. Hirukawa
c 2021 StataCorp LLC

st0630

124

msreg

and Prokhorov (2018) derive the imputation bias analytically and propose two biascorrected estimators. In this article, we introduce a new command, msreg, that implements both estimators in Stata.
Section 2 documents theoretical backgrounds for the msreg command. Section 3
discusses the msreg syntax and provides a numerical example. Section 4 contains asimulation study. Section 5 provides an empirical application by estimating the returnto schooling as in Hirukawa and Prokhorov (2018).

2

Setup and estimators

2.1

Setup and assumptions

Suppose that we are interested in fitting a linear regression model,
y = β0 + x01 β 1 + x02 β 2 + z0 γ + u = w0 θ + u

E(u|w) = 0

where x1 ∈ Rd1 , x2 ∈ Rd2 , z ∈ Rdz , w = (1, x01 , x02 , z0 ), and θ = (β0 , β 1 0 , β 2 0 , γ 0 )0 .
If we can observe all the variables in one sample, OLS is a consistent estimator for θ.
However, in reality, we often encounter a situation where the variables are taken fromtwo different samples. To be precise, we need more notation to distinguish between thetwo samples. The first sample is denoted by S1 = {(yi , x1i , zi )}ni=1 . The second sample
is denoted by S2 = {(x2j , zj )}m
j=1 . For inference, we also denote d3 as the number ofcontinuous common variables in z hereafter, which is not always equal to dz .
Estimation theories in Hirukawa and Prokhorov (2018) are built on a set of assumptions that are required for identification, consistency, and asymptotic normality of theirestimators. Some of them are quite common. For example, assumption 2 imposes compactness of the support of continuous common variables. In our empirical analysis in
section 5, educ, feduc, and meduc are such variables, and it is natural to think of theirsupport as compact. On the other hand, there are more subtle assumptions in Hirukawa
and Prokhorov (2018) that may or may not hold in a given application. Examples include a common joint distribution for S1 and S2 (assumption 1), E (η 1 η 02 ) = 0, and strict
nonlinearity in g2 (·) [assumption 3(ii)], where η ` := x` − g` (z) and g` (z) := E ( x` | z)for ` = 1, 2. It is difficult to test the validity of these assumptions because x1 and x2
belong to two distinct samples.

2.2

Nearest-neighbor matching

The matched sample can be constructed via the nearest-neighbor matching (NNM) using
a vector of common variables z across two samples. Note that z must contain at least
one continuous variable for valid inference; inclusion of discrete common variables witha finite number of support points (for example, binary variables) in z does not affectthe asymptotic results that will be stated shortly.

M. Hirukawa, D. Liu, and A. Prokhorov

125

To specify the NNM, we need to first define a matrix norm to measure the distance
between two vectors. For a vector x and some symmetric and positive definite matrix
A, the vector norm is defined as ||x||A = (x0 Ax)1/2 . Following Abadie and Imbensno−1
PN(2011), we use the Mahalanobis metric AM = (1/N ) i=1 (zi − z)(zi − z)0and
−1
the normalized Euclidean metric AN E = diag A−1, where N = n + m and z =M
PN(1/N ) i=1 zi .
Let jk (i) be the index of the kth match in S2 to the unit i in S1 ; that is, for eachi ∈ {1, . . . , n}, jk (i) satisfiesmX
1 ||zj − zi ||A ≤ ||zjk (i) − zi ||A = k
j=1

In other words, zjk (i) is the kth nearest neighbor in S2 to the unit i in S1 .
For each unit i, let JK (i) = {ji (1), . . . , ji (K)} denote the K matches from S2 . The
NNM-based matched sample isn
S = (yi , x1i , x2j1 (i) , . . . , x2jK (i) , zi ) i=1
P
We also write x2j(i) = (1/K) j∈JK (i) x2j .
For estimation, we use a transformation of the matched sample S.
nS∗ = (yi , x1i , x2j(i) , zi ) i=1
In contrast with the original matched sample S, we replace the individual matchedvariable x2j1 (i) , . . . , x2jK (i) by its mean x2j(i) .
Throughout, it is assumed that we estimate θ by regressing yi on
wi,j(i) = (1, x01i , x2j(i) 0 , zi )0 .

2.3

Inconsistency of matched-sample OLS

We start by using an OLS estimator on the matched sample S∗ . The OLS estimator is−1
bMSOLS = Q
d
d
θ
R
W
W

Pn
Pn
0
d
d
where Q
W = (1/n)i=1 wi,j(i) wi,j(i) and RW = (1/n)
i=1 wi,j(i) yi . It is referred toas the matched-sample OLS (MSOLS) estimator.
Theorem 1 (Hirukawa and Prokhorov [2018], theorem 1). Under some regularity conditions,bMSOLS = Q−1 PW θ + Op m−1/d3 + Op n−1/2θW0where QW = E(wi,j(i) wi,j(i)), PW = QW − (1/K)Σ, Σ is a (d + 1) × (d + 1) blockdiagonal matrix of the form Σ = diag{0(d1 +1)×(d1 +1) , Σ2 , 0dz ×dz }, and Σ2 = E(η 2 η 02 ).

126

msreg

Theorem 1 implies that MSOLS is inconsistent in general. The inconsistency
is atP
tributed to correlation between the imputed regressor x2j(i) and (1/K) j∈JK(i) η2j in
the composite error term i,j(i) . All asymptotic analyses in Hirukawa and Prokhorov(2018) are based on letting n and m diverge while keeping K fixed. It is in principle possible to restore consistency by letting K diverge at a rate slower than n andm. However, a fixed K is what researchers are likely to do in practice; Abadie andImbens (2006) also adopt this setup. Moreover, the Op m−1/d3 term corresponds tothe second-order bias term λi,j(i) because of the matching discrepancy of Abadie andbMSOLS to its
Imbens (2006). Observe that this term affects the convergence rate of θ
probability limit; see remark 3 of Hirukawa and Prokhorov (2018) for more details.

2.4

One-step bias-corrected estimatorpd
The source of the inconsistency of the MSOLS estimator is that Q
W −→ QW , whereaspdRW −→ PW θ = {QW − (1/K)Σ} θ. To eliminate the nonvanishing bias, we replacedd
the denominator Q
W by a consistent estimator of PW and leave the numerator RW unchanged. Because this bias correction has an indirect inference interpretation, Hirukawa
and Prokhorov (2018) call this estimator the matched-sample indirect inference (MSII)destimator. Let P
W be some consistent estimator of PW . Then, the MSII estimator isdefined as−1bMSII = P
d
d
θ
R
W
W

To consistently estimate PW , we need consistent estimators for QW and Σ. Apd
parently, Q
W is a natural estimator for QW . Furthermore, it turns out that we canconsistently estimate Σ without nonparametric estimation of E(x2 |z). To do so, we
first reorder S2 with respect to z in ascending order.
1. Define z(1) as the observation that has the smallest first element; that is, (1) =arg min1≤j≤m zj1 .
2. For j = 2, . . . , m, choose (j) = arg minj6=(1),...,(j−1) ||zj − z(j−1) ||, where the normof a matrix ||A|| is defined as ||A|| = {tr(A0 A)}1/2 .
Given the reordered sample S2 = {(x2(j) , z(j) )}m
j=1 , Σ2 can be consistently estimated
by
c2 =
Σ

m
X
1
∆x2(j) ∆x02(j)2(m − 1) j=2

where ∆x2(j) = x2(j) − x2(j−1) . This is known as the difference-based variance estimator; see von Neumann (1941), Yatchew (1997), and Horowitz and Spokoiny (2001) for
references.

M. Hirukawa, D. Liu, and A. Prokhorov

127

The estimator of PW is given by
n
o
1 b
1
d
d
d
c2 , 0d ×dP
Σ=Qdiag 0(d1 +1)×(d1 +1) , Σ
W = QW −
W −
z
z
K
K
bMSII . The theorem applies
Theorem 2 below documents asymptotic normality of θ
only when the number of continuously distributed matching variables is so small thatthe second-order, matching discrepancy bias can be safely ignored. Observe that bothbMSII and its asymptotic variance depend on the divergencethe convergence rate of θ
pattern of (n, m).
Theorem 2 (Hirukawa and Prokhorov [2018], Corollary 1). Under some regularity
conditions, as n, m → ∞,
√ 


d
−1
bMSII − θ −→n θ
N (0, VI ) = N 0, P−1
if n/m → κ ∈ (0, ∞)

W ΩPW





and d3 = 1











d
√n θ
bMSII − θ −→N (0, VII ) = N 0, P−1 Ω11A P−1
if n/m → 0
W

W












√ b

d
−1


m θ MSII − θ −→ N (0, VIII ) = N 0, P−1
W Ω22 PW
and d3 = 1, 2
if n/m → ∞and d3 = 1

where the definitions of Ω, Ω11A , and Ω22 can be found in the appendix, along with
their consistent estimates.
As demonstrated in this theorem and theorem 3 below, the bias-corrected estimatorsof Hirukawa and Prokhorov (2018) attain the parametric convergence rate only when
the number of continuous common variables is four or less. It may be tempting toinclude as many continuous common variables as possible in the NNM. However, thisresults in slowing down the convergence rate, and we do not recommend it.

2.5

Two-step bias-corrected estimator

The one-step bias-corrected estimator can attain the parametric rate of convergencewith at most two matching variables. To overcome this curse of dimensionality, we
should eliminate the second-order bias λi,j(i) . The entire procedure is reminiscent of
the fully modified least-squares estimation for cointegrating regressions by Phillips andHansen (1990). In this sense, Hirukawa and Prokhorov (2018) call the estimator the
fully modified MSII (MSII-FM) estimator.
Estimating λi,j(i) requires consistent estimates of θ and g2 (·). For θ, we can use
bMSII . For g2 (·), we use a nonparametric power-series estimation
the MSII estimate θ
as in Abadie and Imbens (2011). Let v = (v1 , . . . , vdz ) be a multi-index of dimension

128

msreg

Pdz
dz , which is a dz -dimensional vector of nonnegative integers with |v| = l=1
vl . Also,
vl
dz
v
∞
denote z = Πl=1 zl . Consider a series {vQ }Q=1 containing distinct vectors such that|v(Q)| is nondecreasing. Let pQ (z) = zv(Q) and pQ (z) = {p1 (z), . . . , pQ (z)}0 . Then, anonparametric series estimator of the regression function g2r (z), r = 1, . . . , d2 is

−
m
m
X
 X
Q(m)
gc
(z)0
pQ(m) (zj )pQ(m) (zj )0
pQ(m) (zj )x2r,j2r = p


j=1

j=1

where x2r,j denotes the rth element of x2j in S2 , (·)− denotes the generalized inverse,
and Q = Q(m) implies the dependence of Q on the sample size of S2 .
The entire estimation procedure can be summarized in the following three steps:
1. Run MSII using the matched sample S∗ to obtain

0
0
0
0
0
bMSII = β[[ [ dθII,0 , β II,1 , β II,2 , γIIbi,j(i) }n , where
2. Construct adjusted dependent variables {yi+ }ni=1 = {yi − λ
i=1

0
X
bi,j(i) = gb2 (zi ) − 1[
λgb2 (zj ) β
II,2K
j∈JK (i)n
3. Rerun MSII using the modified matched sample S+ = (yi+ , x1i , x2j(i) , zi ) i=1 toobtain the final estimatorn−1+−1 1 X 0bdddθ MSII−FM = PW RW = PWwy+n i=1 i,j(i) i
Theorem 3 (Hirukawa and Prokhorov [2018], theorem 4). Under some regularity conditions, as n, m → ∞,√ dbMSII−FM − θ −→n θ
N (0, VI )if n/m → κ ∈ (0, ∞) and d3 = 2, 3






√ 

d
bMSII−FM − θ −→n θ
N (0, VII )if n/m → 0 and d3 = 3, 4







√ 
d

bMSII−FM − θ −→m θ
N (0, VIII )if n/m → ∞ and d3 = 2, 3
where the definitions of VI , VII , and VIII are the same as in theorem 2.
In practice, the standard errors (SEs) resulting from each of the three cases may
be quite different. The relative magnitudes of n and m determine which case applies.
Hirukawa and Prokhorov (2018) did not provide any generic comparisons for the variancematrices. Besides the scaling factor, the differences can be attributed to the specificfeatures of the datasets and model specification. In borderline cases, it is advisable touse larger SEs for conservative inference.

M. Hirukawa, D. Liu, and A. Prokhorov

3

The msreg command

3.1

Syntax129

msreg has the following syntax.
   msreg depvar varlist X1 (varlist X2 = varlist Z ) using filename ifin, vce(vce spec) estimator(est spec) nneighbor(#) metric(metric spec)order(#) noconstant level(#) display options coeflegend

3.2

Optionsvce(vce spec) specifies the type of variance–covariance matrix used in computation.
vce spec can be one of vi, vii, or viii. The default is vce(vi). The definition ofvi, vii, and viii can be found in theorem 2.
estimator(est spec) specifies the type of estimator. est spec can be either onestep or
twostep. onestep specifies to use the one-step bias-corrected estimator. twostep
specifies to use the two-step bias-corrected estimator. The default isestimator(twostep).1
nneighbor(#) specifies the number of matches per observation. The default isnneighbor(1). The maximum allowed number of matches is 10. Each observation
is matched with the mean of the specified number of observations from the otherdataset.
metric(metric spec) specifies the distance matrix that is used as the weight matrix ina quadratic form that transforms the multiple distances into a single distance measure. metric spec can be either mahalanobis or euclidean. metric(mahalanobis)
specifies to use the inverse of the sample covariance matrix of matching variables,
which is the default. metric(euclidean) specifies to use the inverse of only diagonal
elements of the sample covariance matrix of matching variables.
order(#) specifies the order of polynomials in the power-series approximation for MSIIFM. The default is order(2). The maximum allowed number of order is 5.
noconstant suppresses the constant term.
level(#) specifies the level of significance for the output table.
display options: noci, nopvalues, noomitted, vsquish, noemptycells, baselevels,
allbaselevels, nofvlabel, fvwrap(#), fvwrapon(style), cformat(%fmt),
pformat(%fmt), sformat(%fmt), and nolstretch; see [R] Estimation options.
coeflegend specifies that the legend of the coefficients and how to specify them in an
expression be displayed rather than displaying the statistics for the coefficients.
1. We intentionally leave the option estimator(msols) undocumented. MSOLS is inconsistent and
we include it in the article for simulations only. We do not recommend using the option estimator(msols).

130

msreg

3.3

Stored results

msreg stores the following in e():Scalarse(N 1)e(N 2)e(nneighbor)e(chi2)e(p)e(df m)e(order)
Macrose(cmd)e(title)e(chi2type)e(vce)e(vcetype)e(properties)e(footnote)e(common)e(matched)e(metric)e(estimator)

msreg
title in estimation output
Wald
vce spec specified in vce()
title used to label Std. Err.
b V
footnote displayed under the output tablecommon variables that exist in both samples
matched variables
type of distance matrix
type of estimator

Matricese(b)e(V)

coefficient vector
variance–covariance matrix of the estimators

Functionse(sample)marks estimation sample

3.4

number of observations in the first sample
number of observations in the second sample
number of nearest neighbors matchedchi-squared
p-value for test of variables
degree of freedom in the modelorder of polynomial in the power-series estimation

A numerical example

We illustrate the use of msreg with a numerical example.
For illustration, we simulate two datasets: s1.dta and s2.dta. The first sample,
s1.dta, contains the dependent variable y and some independent variables, x11, x12,
z1, and z2. The second sample, s2.dta, contains some other dependent variables x21,
x22, z1, and z2. Notice that variables z1 and z2 exist in both samples. In contrast, thevariables x21 and x22 exist only in the second sample, s2.dta. The data-generatingprocess is described in section 4.
We want to fit the following regression model:
y = β0 + β11 x11 + β12 x12 + β21 x21 + β22 x22 + γ1 z1 + γ2 z2 + u(1)

The true values of all the coefficients are set to 1.
Apparently, we cannot estimate (1) using just the first sample, s1.dta, because thevariables x21 and x22 are missing in this dataset. Instead, we want to use the common
variables that exist in both samples, that is, z1 and z2, to construct matched variables
x22 and x21 from the second sample, s2.dta.

M. Hirukawa, D. Liu, and A. Prokhorov131

We now use msreg to estimate the coefficients in (1). We can use the default two-stepbias-corrected estimator and the default vi type variance if we assume n/m converges
to a nonzero constant.
. use s1
. msreg y x11 x12 (x21 x22 = z1 z2) using s2, nneighbor(2) order(3)Matched sample regression
Number of obs in sample 1 =Number of obs in sample 2 =Number of matches=
Order of power-series=
Estimator: twostep
Wald chi2(6)=Distance metric : mahalanobis
Prob > chi2
=

x11
x12
x21
x22
z1
z2
_cons

250
750
2
3
5060.150.0000

Coef.

VI
Std. Err.

z

P>|z|

[95% Conf. Interval].8993908.9877756
1.103668
1.063374
1.105623.7304006.7213349.1602669.1594351.0568935.161514.2745877.2866126.2237863

5.61
6.20
19.40
6.58
4.03
2.55
3.22

0.000
0.000
0.000
0.000
0.000
0.011
0.001

.5852736.6752885.9921588.7468123.5674414.1686502.2827219

1.213508
1.300263
1.215177
1.379936
1.643805
1.292151
1.159948

Common variables : z1 z2
Matched variables: x21 x22
Note: VI Std. Err. assumes that sample size ratio between sample 1 and sample
2 converges to a nonzero constant and there are at most three continuous
common variables.

Here are some comments on the syntax.
• The (x21 x22 = z1 z2) specifies that the variables x21 and x22 are the variables
to be matched, and the variables z1 and z2 are the common variables that exist
in both samples.
• using s2 specifies that the variables x21 and x22 come from s2.dta.
• We use the default two-step bias-corrected estimator.
• The default option vce(vi) specifies to use the vi-type variance matrix as specified in theorem 2 because we assume the sample-size ratio between the two samples
converges to a nonzero constant and there are only two continuous common variables used for matching. Actually, the note states a similar explanation aboutvce(vi).
• Option nneighbor(2) specifies to pick out 2 matches via the NNM.
• Option order(3) specifies to fit a third-order polynomial in the power-series approximation for MSII-FM.
• The output shows the point estimates of coefficients and their SEs, and they canbe interpreted as in a regular linear regression framework.

132

msreg

4

A simulation study

4.1

Simulation design

We conduct a Monte Carlo simulation study for two purposes: first, we want to see
the finite-sample performance of MSII-FM in contrast to MSOLS; second, the simulation
results can serve as a verification of the numerical implementation of our commandmsreg. The simulation study replicates that of Hirukawa and Prokhorov (2018).
The model considered throughout isy = β0 + x01 β 1 + x02 β 2 + z0 γ + u(2)where x1 = (x11 , x12 )0 , β 1 = (β11 , β12 )0 ∈ R2 , x2 = (x21 , x22 )0 , β 2 = (β21 , β22 )0 ∈ R2 ,and z = (z1 , . . . , zd3 )0 , γ = (γ1 , . . . , γd3 )0 ∈ Rd3 for d3 = 1, 2, 3. We assume that
two samples, namely, S1 = {(yi , x1i , zi )}ni=1 and S2 = {(x2j , zj )}m
j=1 are observable inpractice.
Here is how to generate the data. First, z∗ = (z1∗ , z2∗ , z3∗ ) is generated by√√   
1√
1/ 2 √1/ √30i.i.d.
z∗ ∼ N 0 , 1/√2 √ 1√
2/ 3
0
1/ 3
2/ 3
1
Each zp∗ (p = 1, 2, 3) is transformed to zp = 4Φ(zp∗ ) − 2, where Φ(·) is the cdf of N (0, 1).
Notice that zp are mutually correlated U [−2, 2] random variables. For a given d3 , thezp (p ≤ d3 ) are used as matching variables.
Pd3
Second, x1 = (x11 , x12 )0 is generated by x1q =p=1 zp + ηq (q = 1, 2), wherePd30ηq ∼ N (0, 1). Third, x2 = (x21 , x22 ) is generated by x2r = p=1
g2r (zp )+η2r (r = 1, 2)
for some nonlinear function g2r (·), where η2r ∼ N (0, 1). Specifically,
p g21 (z) = z +(5/τ ) φ(z/τ ), τ = 0.25, where φ(·) is the pdf of N (0, 1), and g22 (z) = 4 |z/2|(1 − |z/2|)
sin{2π(1 + )/(|z/2| + )}, with  = 0.05.
i.i.d.

Finally, y is generated by setting all coefficients equal to 1 with u ∼ N (0, 1). The
sample sizes are set to (n, m) = (1000, 1000). The number of replications is 1,000.
We focus on the finite-sample properties of estimators of β22 and γ1 . For eachestimator, the following performance measures are computed: i) Bias (1 − Mean), where
Mean is the simulation average of the parameter estimate; ii) standard deviation (SD)(simulation SD of the parameter estimate); iii) SE (simulation average of the SE); andiv) Rej. rate (rejection rate for the test of parameter estimate equal to its true value 1
against the nominal 5% level of significance).
For d3 = 1, 2, 3, we estimate the coefficients in (2) using both MSOLS and MSII-FM.
For MSOLS, the number of matches is K = 1, 2, 4, 8. For MSII-FM, the number of matches
K is fixed at 1, and orders of polynomials in the power-series approximation are 2, 3,
and 4. For a more complete simulation study, see section 4 of Hirukawa and Prokhorov(2018).

M. Hirukawa, D. Liu, and A. Prokhorov

4.2

133

Results

The simulation results are summarized in tables 1 and 2 for MSOLS and MSII-FM, respectively.
Table 1 shows that regardless of the number of matches, there is a big bias of β22 ,and there is a large rejection rate, which indicate inconsistency of MSOLS as implied intheorem 1.
Table 1. Monte Carlo results for MSOLS
K

β22
1

2

4

8

γ1
1

2

4

8

0.44860.04260.05071.00000.28660.04550.05231.00000.16340.04740.05460.91100.07730.04920.06030.3800−0.16800.1122
0.11740.3210−0.09190.10520.11410.1550−0.04420.10050.11190.1070−0.0216
0.09800.11090.09100.52800.04080.04621.00000.37240.04600.05291.00000.22390.05230.06270.96600.08280.06190.07540.3110−0.1289
0.14680.15480.1640−0.0723
0.13980.15020.1030−0.0372
0.13690.15340.08800.00930.13790.16670.11000.75320.04720.05141.00000.61490.05750.06481.00000.44720.07070.07941.00000.23050.08930.10820.6860−0.2206
0.19930.2015
0.1980−0.1306
0.1906
0.19900.1210−0.0475
0.18810.20780.08600.02810.19220.22700.0920d3 = 1Bias
SD
SE
Rej. rate
d3 = 2Bias
SD
SE
Rej. rate
d3 = 3Bias
SD
SE
Rej. rate

Table 2 shows that 1) the bias is small; that is, the mean of the point estimates isvery close to its true value; 2) the SD of the point estimate is very close to the meanof the SEs; and 3) the overall rejection rate is close to the nominal 5% level. Noticethat for the case d3 = 2, the rejection rate is a little bit off for β22 . However, based on
results for larger samples, it seems that the over-rejection rate is due to the finite-sample
bias of MSII-FM (reported in the Supplement of Hirukawa and Prokhorov [2018]). The
simulation result shows that MSII-FM performs well in a finite sample as predicted bytheorem 3 and that it also numerically verifies the implementation of msreg.

134

msreg
Table 2. Simulation results for MSII-FM

Order

β22
2

3

4

γ1
2

3

4

−0.0305
0.10490.11300.0620−0.0305
0.10490.11480.0610−0.03230.10470.11490.06400.01100.12570.13530.07300.01100.12590.13610.07300.01190.12580.13590.0720−0.17400.15390.16370.1400−0.1735
0.15410.16360.1380−0.16410.15400.16430.12700.03180.17500.19120.08200.03820.17540.19410.08400.04010.17650.19300.0760−0.09480.28840.30410.0370−0.0904
0.29250.31520.0350−0.0866
0.29040.31320.03700.03720.24810.26800.06100.04080.24990.27490.06900.0526
0.24950.2736
0.0650d3 = 1Bias
SD
SE

Rej. rate
d3 = 2Bias
SD
SE

Rej. rate
d3 = 3Bias
SD
SE

Rej. rate

5

An empirical application: Return to schooling

We now apply msreg to a version of Mincer’s (1974) wage equation. We consider thefollowing wage regression,log(wage) = β0 + β1 expr + β2 expr2 + β3 kww + β4 educ + β5 feduc + β6 meduc+ β7 smsa + β8 south + β9 black + u(3)
where expr is years of experience, educ is years of education, kww is Knowledge of Worldof Work test score, feduc and meduc are years of father’s and mother’s education, andblack, smsa, and south are dummy variables to indicate whether an individual is black,lives in an urban area, and lives in the south, respectively.

M. Hirukawa, D. Liu, and A. Prokhorov

135

We can estimate (3) using only card.dta from Card (1995), as in the benchmark
OLS result below.
. use card
. regress lnwage exper c.exper#c.exper> kww educ feduc meduc i.smsa i.south i.black
SourceSS
df
MS

Number of obs
F(9, 2181)
Prob > F
R-squared
Adj R-squared
Root MSE

=
=
=
=
=
=

2,191
89.96
0.00000.27070.2677.37603

Model
Residual

114.487495
308.394118

9
2,181

12.7208327.141400329

Total

422.881613

2,190.193096627

lnwageCoef.

Std. Err.

t

P>|t|

[95% Conf. Interval]

exper.0786594.0082641

9.52

0.000

.062453.0948658c.exper#c.exper-.0021546.0004049

-5.32

0.000

-.0029486-.0013607

kww
educ
feduc
meduc
1.smsa1.south
1.black_cons.0055563.061174-.0018021.0070962.1516693-.1111491-.1320785
4.68611.0013232.0051257.0030073.0035015.0186575.0177506.0254406.0815838

4.20
11.93
-0.60
2.03
8.13
-6.26
-5.19
57.44

0.000
0.000
0.549
0.043
0.000
0.000
0.000
0.000

.0029615.0511222-.0076996.0002295.1150809-.1459589-.1819688
4.52612.0081511.0712258.0040955.013963.1882577-.0763393-.0821882
4.8461. estimate store ols

The estimation result is stored as ols.
Nonetheless, we pretend that the variable kww is missing in this dataset. In accordance with this scenario, we use yet another dataset, wage2.dta, from Blackburn and
Neumark (1992). The dataset contains six variables—educ, feduc, meduc, smsa, south,and black—other than kww. All six variables are used as matching variables to imputethe missing kww, where educ, feduc, and meduc are assumed to be continuous. Our aim
is to see how the estimation result of (3) changes if kww is imputed from wage2.dta.

136

msreg
. msreg lnwage exper c.exper#c.exper>(kww = educ feduc meduc i.smsa i.south i.black)>using wage2, order(3)Matched sample regression
Number of obs in sample 1
Number of obs in sample 2
Number of matches
Order of power-series
Estimator: twostep
Wald chi2(9)Distance metric : mahalanobis
Prob > chi2

=
=
=
=
=
=

2,191
470
1
3
785.110.0000

Coef.

VI
Std. Err.

z

P>|z|

[95% Conf. Interval]

exper.0840744.0082734

10.16

0.000

.0678589.1002898

c.exper#c.exper-.0021476.0004202

-5.11

0.000

-.0029713-.001324

kww
educ
feduc
meduc
1.smsa1.south
1.black_cons.0020755.0680349-.0005824.0070538.1557546-.1114039-.1517787
4.654032.0052253.0079751.003143.0039251.0186934.0178276.0354906.1187133

0.40
8.53
-0.19
1.80
8.33
-6.25
-4.28
39.20

0.691
0.000
0.853
0.072
0.000
0.000
0.000
0.000

-.0081658.0524041-.0067427-.0006393.1191161-.1463453-.221339
4.421359.0123168.0836658.0055778.0147469.192393-.0764626-.0822183
4.886706

Common variables : educ feduc meduc 0b.smsa 1.smsa 0b.south 1.south 0b.black
1.black
Matched variables: kww
Note: VI Std. Err. assumes that sample size ratio between sample 1 and sample
2 converges to a nonzero constant and there are at most three continuous
common variables.
. estimate store twostep_vi

We use the default vi type covariance estimation assuming the sample-size ratio
between the two datasets converge to a nonzero constant. We use a third-order polynomial in power-series estimation to remove the second-order bias. The estimation result
is stored as twostep vi.

M. Hirukawa, D. Liu, and A. Prokhorov

137

We can now compare these two estimation results.
. estimate table ols twostep_vi, b star
Variableols

twostep_vi

exper.0786594***.08407435***

c.exper#c.exper-.00215464***-.00214765***

kww
educ
feduc
meduc

.00555626***.06117396***-.00180207.00709624*.00207551.06803494***-.00058241.00705379

smsa
1.15166928***.15575459***

south
1

-.11114906***-.11140395***

black
1

-.13207849***-.15177868***_cons

4.6861104***

4.6540324***legend: * p<0.05; ** p<0.01; *** p<0.001

The first column shows the benchmark OLS results. The signs of expr, expr2 , kww,
and educ are as expected, and the estimates are significant at the 5% level.
The second column shows the results from the two-step bias-corrected estimatorwith the default vi-type covariance. All the point estimates have the same sign as in
the OLS benchmark. However, the coefficient on kww is insignificant because of a large
SE.

6Conclusion

In this article, we described a new command, msreg, that implements two estimators
proposed in Hirukawa and Prokhorov (2018). The command allows users to obtainconsistent estimators of linear regression models after imputing missing regressors via
the NNM. We illustrated the use of msreg through a numerical example and an empiricalapplication.

138

7

msreg

Programs and supplemental materials

To install a snapshot of the corresponding software files as they existed at the time of
publication of this article, type. net sj 21-1
. net install st0630. net get st0630

8(to install program files, if available)(to install ancillary files, if available)

References

Abadie, A., and G. W. Imbens. 2006. Large sample properties of matching estimators foraverage treatment effects. Econometrica 74: 235–267. https://doi.org/10.1111/j.14680262.2006.00655.x.
. 2011. Bias-corrected matching estimators for average treatment effects. Journalof Business & Economic Statistics 29: 1–11. https://doi.org/10.1198/jbes.2009.07333.
Blackburn, M., and D. Neumark. 1992. Unobserved ability, efficiency wages, andinterindustry wage differentials. Quarterly Journal of Economics 107: 1421–1436.
https: // doi.org / 10.2307 / 2118394.
Card, D. E. 1995. Using geographic variation in college proximity to estimate thereturn to schooling. In Aspects of Labour Market Behaviour: Essays in Honour ofJohn Vanderkamp, ed. L. N. Christofides, E. K. Grant, and R. Swidinsky, 201–222.
Toronto, Canada: University of Toronto Press.
Hirukawa, M., and A. Prokhorov. 2018. Consistent estimation of linear regression models
using matched data. Journal of Econometrics 203: 344–358. https://doi.org/10.1016/
j.jeconom.2017.07.006.
Horowitz, J. L., and V. G. Spokoiny. 2001. An adaptive, rate-optimal test of a parametric mean-regression model against a nonparametric alternative. Econometrica 69:
599–631. https: // doi.org / 10.1111 / 1468-0262.00207.
Mincer, J. 1974. Schooling, Experience, and Earnings. New York: National Bureau of
Economic Research.
von Neumann, J. 1941. Distribution of the ratio of the mean square successive difference
to the variance. Annals of Mathematical Statistics 12: 367–395. https://doi.org /10.
1214 / aoms / 1177731677.
Phillips, P. C. B., and B. E. Hansen. 1990. Statistical inference in instrumental variables
regression with I(1) processes. Review of Economic Studies 57: 99–125. https: // doi.
org / 10.2307 / 2297545.
Yatchew, A. 1997. An elementary estimator of the partial linear model. Economic
Letters 57: 135–143. https: // doi.org / 10.1016 / S0165-1765(97)00218-8.

M. Hirukawa, D. Liu, and A. Prokhorov

139

About the authors
Masayuki Hirukawa is a professor in the Faculty of Economics at Ryukoku University in Japan.
Research for this article was supported by a grant from the Japan Society for the Promotionof Science (Project No. 19K01595).
Di Liu is a senior econometrician at StataCorp.
Artem Prokhorov is a professor in the Discipline of Business Analytics at the University ofSydney, Sydney, Australia, and a research fellow at CEBA, St. Petersburg State University,
St. Petersburg, Russia, and CIREQ, University of Montreal, Montreal, Canada. Research
for this article was supported by a grant from the Russian Science Foundation (Project No.
20-18-00365).

A

AppendixbMSII and θ
bMSII−FM , respectively.
Theorems 2 and 3 give the asymptotic distributions of θ
The covariance matrices VI , VII , and VIII depend on the definitions of Ω, Ω11A , andΩ22 . We define Ω, Ω11A , and Ω22 as follows:Ω = Ω11A + κ (β 02 Σ2 β 2 )E(w)E(w)0
1
1+ 2 diag 0(d1 +1)×(d1 +1) , (β 02 Σ2 β 2 )Vg2 + Ψ, 0dz ×dz
K
2
(

0 )
1
1
Ω11A = E
wi,j(i) i,j(i) + Σθwi,j(i) i,j(i) + Σθ
K
K


1
1
Ω22 = 2 diag 0(d1 +1)×(d1 +1) , Ξ + Ψ, 0dz ×dz
K2Vg2 = Var {g2 (z)}
Ξ = E (η 2 η 02 − Σ2 ) β 2 β 02 (η 2 η 02 − Σ2 )
Ψ = (β 02 Σ2 β 2 )Σ2 + Σ2 β 2 β 02 Σ2
We present consistent estimators of Ω11A , Ω22 , and Ω for MSII below. Because
MSII-FM is first-order asymptotically equivalent to MSII as documented in theorem 3,bMSII in Ω
bMSII−FM yields the estimators forb 11A , Ωb 22 , and Ωb with θ
simply replacing θ
MSII-FM,

140

msreg


0
n 
X
1 bb
1 bb
b 11A = 1Ω
wi,j(i) [Σθ MSIIwi,j(i) [Σθ MSIIi,j(i) +i,j(i) +
n i=1
K
K
n
o
b 22 = 1 diag 0(d +1)×(d +1) , Γ(−1)
b
b
b
Ω
+
Γ(0)
+
Γ(1),
0
d
×d
z
z
1
1
K2

n  \ 0b \ bbΩ = Ω11A +β 2,MSII Σ2 β 2,MSII ww0
m

0
1
b \ b
+ 2 diag 0(d1 +1)×(d1 +1) , β\
2,MSII Σ2 β 2,MSII Vg2
K

n
o
b
b
b
+ Γ(0) − Γ(−1) + Γ(1) , 0dz ×dz
0bb[
where [i,j(i) = yi − wi,j(i) θ MSII , β 2,II is the MSII estimate of β 2 , and Γ(l) is the lthnob 2 ; that is,
sample autocovariance of ∆x2j ∆x0 /2 − Σ
2j




min{m,m+l} X0∆x2j−l ∆x02j−l∆x2j ∆x02j
1
b 2 β\
b
\
−Σ
β
−
Σ
2
2,MSII 2,MSII
m−1
22
j=max{2,2+l}
  

1
1
  

   1 Pn
x1   n i=1 x1i 
  

 

w=
  =  1 Pm
x2   m j=1 x2j 
  

  

PN
1
z
z
i=1 i
N

b =
Γ(l)

bg =
V
2

m
1 X
b2
(x2j − x2j )(x2j − x2j )0 − Σ
m − 1 j=1