JEL Classification: C12, C13, C14.
Introduction
Two-way clustering is widely used in empirical economics to account for dependence among observations that share a cluster in either of two dimensions. In a firm panel, for example, outcomes may be correlated both within a firm over time and across firms in the same period. The conventional approach estimates a parameter from the full array and adjusts its standard error for dependence within rows and columns [3 ] . Even for the sample mean, however, the dependence structure can affect the form of the limiting distribution. When row or column effects dominate, the normalized full-array mean typically has a Gaussian limit. As these effects weaken, interactions between rows and columns can produce a non-Gaussian limit. Menzel [24 ] shows that the limiting law of the normalized full-array mean cannot be estimated uniformly consistently over such sequences. We approach the inference problem by constructing a different estimator of the same population mean, for which Gaussian critical values are uniformly valid.
Our estimator averages the means obtained from repeated random matchings. Taking the dimension with fewer clusters as the rows, each matching assigns a distinct column to every row, uniformly at random and independently of the data. No two selected observations share a row or a column. Under separately exchangeable, dissociated sampling, they therefore form an i.i.d. sample from the marginal distribution. A single matching permits conventional inference for the population mean but uses only one observation per row. Additional matchings recover precision: conditional on the observed array, their means are independent and centered at the full-array mean. The choice of how many matchings to use matters because averaging too many can cause the estimator to inherit the non-Gaussian limit of the full-array mean.
We make four contributions. First, we characterize the limiting distribution of the matching estimator as the number of matchings varies. An exact variance decomposition separates the sampling variance of the full-array mean from the additional variance due to matching. Increasing the number of matchings improves precision, but it can also give the full-array error enough weight to produce a non-Gaussian limit. This trade-off motivates letting the number of matchings grow more slowly than the smaller cluster count to recover precision while retaining a Gaussian limit across dependence regimes.
Second, we provide a closed-form variance estimator that makes the confidence interval straightforward to compute. It combines estimates of row, column, and intersection variation with nonnegative weights. The estimate is positive semidefinite in every sample by construction. Inference then requires only the estimated standard error and a standard normal critical value, without a bootstrap or a preliminary test of the dependence regime.
Third, we establish uniformly asymptotically exact inference for the population mean. Under our moment conditions, the studentized matching estimator admits a uniform standard normal approximation whenever the number of matchings is fixed or grows at any rate slower than the smaller cluster count. If the number of matchings grows more slowly than the square root of that count, a uniform bound on the standardized fourth moment suffices for the same conclusion. The approximation holds uniformly over sequences of data-generating processes that may vary with sample size. The strength of clustering in either dimension need not converge, and the two cluster counts may grow at different rates.
Fourth, we extend the procedure to smooth functions of generalized method of moments (GMM) parameters, including regression coefficients. Under two-way clustering, different moments can converge at different rates. Estimation error in a slowly converging coefficient can then enter the limiting distribution of a faster one through a quadratic term (Davezies et al. [13 ] , Section 3). We show that the same restriction on the number of matchings makes this term negligible relative to the standard error of the target. Applying the variance estimator to the estimated influence function therefore gives uniformly asymptotically exact inference even when the moment components follow different dependence regimes.
Our analysis complements the literature on inference based on the full array. Davezies et al. [12 ] develop empirical process theory for exchangeable arrays, and Chiang and Sasaki [11 ] give conditions under which conventional two-way standard errors yield valid Gaussian inference. MacKinnon et al. [23 ] study wild bootstrap inference for regression, while Chiang, Matsushita, and Otsu (2025) construct an empirical likelihood statistic using estimates that leave out rows or columns. Davezies et al. [13 ] develop analytic tests that control size uniformly over their class. Their tests are asymptotically exact in Gaussian regimes and can be conservative in non-Gaussian regimes. By changing the estimator, we obtain exact asymptotic coverage even for laws under which the full-array mean has a non-Gaussian limit, at a cost in precision that we quantify.
Menzel [24 ] shows that, when the observation pattern is independent of the data, observing only a vanishing fraction of intersections can yield a Gaussian limit for the sample mean. Chiang, Li, and Sasaki [9 ] establish asymptotic normality under Bernoulli subsampling for each fixed data-generating process. We select intersections by repeated matching and establish asymptotically exact coverage uniformly over laws in which row and column effects may weaken or vanish as the sample size grows. The role of the number of sampled terms in determining the limiting distribution also arises in the theory of incomplete U U U -statistics [19 , 6 ] .
Section 2 introduces the two-way clustering framework and describes how the dependence structure affects the limiting distribution of the full-array mean. Section 3 presents the matching estimator, shows how the number of matchings affects its precision and limiting distribution, and establishes uniformly asymptotically exact coverage for the resulting interval. Section 4 extends the results to smooth functions of GMM parameters. Section 5 reports Monte Carlo evidence on confidence interval coverage and length. Section 6 concludes. The appendix contains the proofs.
Two-way clustering
Sampling framework
We observe a complete array Y N T = { Y N T , i t : 1 ≤ i ≤ N , 1 ≤ t ≤ T } Y_{NT} = \{Y_{NT,it} : 1 \le i \le N, 1 \le t \le T\} Y N T = { Y N T , i t : 1 ≤ i ≤ N , 1 ≤ t ≤ T } . Label the smaller dimension as the rows, so that 2 ≤ N ≤ T 2 \le N \le T 2 ≤ N ≤ T , and let N , T → ∞ N, T \to\infty N , T → ∞ jointly.
Assumption 1 (Two-way clustered sampling). At each ( N , T ) (N,T) ( N , T ) , the observations admit the representation
Y N T , i t = f N T ( ξ i , η t , ζ i t ) . (1) Y_{NT,it}=f_{NT}(\xi_i,\eta_t,\zeta_{it}). \tag*{(1)} Y N T , i t = f N T ( ξ i , η t , ζ i t ) . (1) The row shocks { ξ i } \{\xi_i\} { ξ i } , column shocks { η t } \{\eta_t\} { η t } , and intersection shocks { ζ i t } \{\zeta_{it}\} { ζ i t } form three mutually independent families, with i.i.d. shocks within each family. The marginal mean is θ N T = E [ Y N T , i t ] \theta_{NT}=E[Y_{NT,it}] θ N T = E [ Y N T , i t ] , and the marginal variance satisfies σ N T 2 = Var ( Y N T , i t ) ∈ ( 0 , ∞ ) \sigma^2_{NT}=\operatorname{Var}(Y_{NT,it})\in(0,\infty) σ N T 2 = Var ( Y N T , i t ) ∈ ( 0 , ∞ ) .
Equation (1) is the Aldous-Hoover-Kallenberg representation for a separately exchangeable, dissociated array [1 , 16 , 22 ] . Observations may be dependent within a row or column, while collections of observations with disjoint row and column indices are independent. Allowing f N T f_{NT} f N T to vary with ( N , T ) (N,T) ( N , T ) permits both the marginal distribution and the strength of clustering to change with sample size. For related work on multiway clustering, see, e.g., (Davezies et al. 2021, 2025; Menzel 2021; Chiang et al. 2023, 2024; Graham et al. 2024; Chen and Vogelsang 2024; Hounyo and Lin 2026a,b; Juodis 2026; Jochmans 2026).
To describe how each shared shock affects the conditional mean, define the row projection and column projection by
a N T , i = E [ Y N T , i t ∣ ξ i ] − θ N T , b N T , t = E [ Y N T , i t ∣ η t ] − θ N T . a_{NT,i}=E[Y_{NT,it}\mid\xi_i]-\theta_{NT},\qquad b_{NT,t}=E[Y_{NT,it}\mid\eta_t]-\theta_{NT}. a N T , i = E [ Y N T , i t ∣ ξ i ] − θ N T , b N T , t = E [ Y N T , i t ∣ η t ] − θ N T .
The remaining variation is the intersection residual r N T , i t r_{NT,it} r N T , i t , giving the following decomposition
Y N T , i t − θ N T = a N T , i + b N T , t + r N T , i t . (2) Y_{NT,it}-\theta_{NT}=a_{NT,i}+b_{NT,t}+r_{NT,it}. \tag*{(2)} Y N T , i t − θ N T = a N T , i + b N T , t + r N T , i t . (2)
All three components are centered and mutually uncorrelated. The intersection residual satisfies E [ r N T , i t ∣ ξ i ] = E [ r N T , i t ∣ η t ] = 0 E[r_{NT,it}\mid\xi_i]=E[r_{NT,it}\mid\eta_t]=0 E [ r N T , i t ∣ ξ i ] = E [ r N T , i t ∣ η t ] = 0 , although it can still depend on the interaction between ξ i \xi_i ξ i and η t \eta_t η t . Distinct intersection residuals are uncorrelated.
Averaging (2) gives the corresponding decomposition of the full-array mean Y ˉ N T = \bar{Y}_{NT}= Y ˉ N T =
( N T ) − 1 ∑ i , t Y N T , i t , (NT)^{-1}\sum_{i,t}Y_{NT,it}, ( N T ) − 1 i , t ∑ Y N T , i t ,
Y ˉ N T − θ N T = 1 N ∑ i = 1 N a N T , i + 1 T ∑ t = 1 T b N T , t + 1 N T ∑ i = 1 N ∑ t = 1 T r N T , i t . \bar{Y}_{NT}-\theta_{NT}=\frac{1}{N}\sum_{i=1}^{N}a_{NT,i}+\frac{1}{T}\sum_{t=1}^{T}b_{NT,t}+\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}r_{NT,it}. Y ˉ N T − θ N T = N 1 i = 1 ∑ N a N T , i + T 1 t = 1 ∑ T b N T , t + N T 1 i = 1 ∑ N t = 1 ∑ T r N T , i t .
Let σ N T , R 2 = E [ a N T , i 2 ] \sigma^{2}_{NT,R}=E[a^{2}_{NT,i}] σ N T , R 2 = E [ a N T , i 2 ] , σ N T , C 2 = E [ b N T , t 2 ] \sigma^{2}_{NT,C}=E[b^{2}_{NT,t}] σ N T , C 2 = E [ b N T , t 2 ] , and σ N T , I 2 = E [ r N T , i t 2 ] \sigma^{2}_{NT,I}=E[r^{2}_{NT,it}] σ N T , I 2 = E [ r N T , i t 2 ] , where the subscripts R, C, and I denote the row, column, and intersection components, respectively. By orthogonality,
σ N T 2 = σ N T , R 2 + σ N T , C 2 + σ N T , I 2 , V N T , F = Var ( Y ˉ N T ) = σ N T , R 2 N + σ N T , C 2 T + σ N T , I 2 N T . (3) \sigma^{2}_{NT}=\sigma^{2}_{NT,R}+\sigma^{2}_{NT,C}+\sigma^{2}_{NT,I},\qquad V_{NT,F}=\operatorname{Var}(\bar{Y}_{NT})=\frac{\sigma^{2}_{NT,R}}{N}+\frac{\sigma^{2}_{NT,C}}{T}+\frac{\sigma^{2}_{NT,I}}{NT}. \tag*{(3)} σ N T 2 = σ N T , R 2 + σ N T , C 2 + σ N T , I 2 , V N T , F = Var ( Y ˉ N T ) = N σ N T , R 2 + T σ N T , C 2 + N T σ N T , I 2 . (3)
Limiting distributions across dependence regimes
By (3), the row and column projections account for the following share of the variance of the full-array mean:
ν N T = σ N T , R 2 / N + σ N T , C 2 / T V N T , F , 1 − ν N T = σ N T , I 2 / ( N T ) V N T , F . (4) \nu_{NT}=\frac{\sigma^{2}_{NT,R}/N+\sigma^{2}_{NT,C}/T}{V_{NT,F}},\qquad1-\nu_{NT}=\frac{\sigma^{2}_{NT,I}/(NT)}{V_{NT,F}}. \tag*{(4)} ν N T = V N T , F σ N T , R 2 / N + σ N T , C 2 / T , 1 − ν N T = V N T , F σ N T , I 2 / ( N T ) . (4)
We call a sequence nondegenerate if ν N T → 1 \nu_{NT}\to1 ν N T → 1 : the projection averages dominate and have a Gaussian limit after standardization under the Lindeberg condition for their sum. In the degenerate case, ν N T → 0 \nu_{NT}\to0 ν N T → 0 and the intersection term dominates. Exact degeneracy , where σ N T , R 2 = σ N T , C 2 = 0 \sigma^{2}_{NT,R}=\sigma^{2}_{NT,C}=0 σ N T , R 2 = σ N T , C 2 = 0 at every sample size, is a special case. In the intermediate case, ν N T → c ∈ ( 0 , 1 ) \nu_{NT}\to c\in(0,1) ν N T → c ∈ ( 0 , 1 ) : both terms contribute to the variance. These cases describe the full-array mean; our uniform inference result does not require ν N T \nu_{NT} ν N T to converge.
Even under exact degeneracy, the limiting law can be Gaussian or non-Gaussian. Consider two designs with θ N T = 0 \theta_{NT}=0 θ N T = 0 and σ N T , I 2 = σ N T 2 = 1 \sigma^{2}_{NT,I}=\sigma^{2}_{NT}=1 σ N T , I 2 = σ N T 2 = 1 , so that V N T , F = 1 / ( N T ) V_{NT,F}=1/(NT) V N T , F = 1/ ( N T ) in both. For Y N T , i t = ζ i t Y_{NT,it}=\zeta_{it} Y N T , i t = ζ i t , with ζ i t \zeta_{it} ζ i t standard normal distributed, N T Y ˉ N T \sqrt{NT}\bar{Y}_{NT} N T Y ˉ N T is standard normal. For Y N T , i t = ξ i η t Y_{NT,it}=\xi_i\eta_t Y N T , i t = ξ i η t with mutually independent standard normal row and column shocks, Y ˉ N T \bar{Y}_{NT} Y ˉ N T is the product of the two shock averages. Hence
N T Y ˉ N T = ( N ξ ˉ N ) ( T η ˉ T ) = d Z 1 Z 2 , \sqrt{NT}\bar{Y}_{NT}=(\sqrt{N}\bar{\xi}_{N})(\sqrt{T}\bar{\eta}_{T})\mathrel{\overset{d}{=}} Z_{1}Z_{2}, N T Y ˉ N T = ( N ξ ˉ N ) ( T η ˉ T ) = d Z 1 Z 2 ,
where Z 1 Z_1 Z 1 and Z 2 Z_2 Z 2 are independent standard normals.
A local shift moves the product design into the intermediate regime. Let
Y N T , i t = ξ i ( η t + δ / T ) . (5) Y_{NT,it} = \xi_i(\eta_t + \delta/\sqrt{T}). \tag*{(5)} Y N T , i t = ξ i ( η t + δ / T ) . (5)
For fixed δ ≠ 0 \delta\ne0 δ = 0 , a N T , i = δ ξ i / T a_{NT,i} = \delta\xi_i/\sqrt{T} a N T , i = δ ξ i / T , b N T , t = 0 b_{NT,t} = 0 b N T , t = 0 , and r N T , i t = ξ i η t r_{NT,it} = \xi_i\eta_t r N T , i t = ξ i η t . Although σ N T , R 2 = δ 2 / T → 0 \sigma^2_{NT,R} = \delta^2/T \to0 σ N T , R 2 = δ 2 / T → 0 , the row projection contributes δ 2 / ( N T ) \delta^2/(NT) δ 2 / ( N T ) to V N T , F V_{NT,F} V N T , F , compared with 1 / ( N T ) 1/(NT) 1/ ( N T ) from the intersection residual. Thus ν N T = δ 2 / ( 1 + δ 2 ) \nu_{NT} = \delta^2/(1+\delta^2) ν N T = δ 2 / ( 1 + δ 2 ) , and
Y ˉ N T − θ N T V N T , F = d Z 1 ( Z 2 + δ ) 1 + δ 2 . \frac{\bar{Y}_{NT}-\theta_{NT}}{\sqrt{V_{NT,F}}} \overset{d}{=} \frac{Z_1(Z_2+\delta)}{\sqrt{1+\delta^2}}. V N T , F Y ˉ N T − θ N T = d 1 + δ 2 Z 1 ( Z 2 + δ ) .
Thus a column shift of order T − 1 / 2 T^{-1/2} T − 1/2 changes the law of the normalized full-array mean. For example, δ = 0 \delta= 0 δ = 0 and δ = 1 \delta= 1 δ = 1 yield different limiting distributions, even though the data cannot consistently distinguish between the two sequences. Consequently, no estimator of this distribution is uniformly consistent over a class containing both sequences [24 Proposition 4.1] . The next section uses random matching to estimate the same θ N T \theta_{NT} θ N T with a studentized statistic that admits a uniform Gaussian approximation.
Matching and uniform inference
One matching
A matching selects one observation from each row by assigning a different column to every row. To draw the first matching, start with all T T T columns available. For row 1, choose uniformly from the T T T columns. For each subsequent row i = 2 , … , N i = 2,\ldots,N i = 2 , … , N , choose uniformly from the T − i + 1 T-i+1 T − i + 1 remaining columns. These draws use randomness independent of the array. Write π 1 ( i ) \pi_1(i) π 1 ( i ) for the column assigned to row i i i , leaving the array dimensions ( N , T ) (N,T) ( N , T ) implicit in the matching notation. By construction,
i ≠ j ⟹ π 1 ( i ) ≠ π 1 ( j ) . (6) i \ne j \;\Longrightarrow\; \pi_1(i) \ne\pi_1(j). \tag*{(6)} i = j ⟹ π 1 ( i ) = π 1 ( j ) . (6)
Conditional on the draw, the N N N selected observations come from distinct rows and distinct columns, so they use disjoint triples of shocks in (1) . They are therefore conditionally i.i.d. with the same marginal law as Y N T , 11 Y_{NT,11} Y N T , 11 .
The matching mean and its sample variance are
M N T , 1 = 1 N ∑ i = 1 N Y N T , i π 1 ( i ) , S N T , 1 2 = 1 N − 1 ∑ i = 1 N ( Y N T , i π 1 ( i ) − M N T , 1 ) 2 . (7) M_{NT,1} = \frac{1}{N}\sum_{i=1}^{N}Y_{NT,i\pi_1(i)}, \qquad S^2_{NT,1} = \frac{1}{N-1}\sum_{i=1}^{N}\left(Y_{NT,i\pi_1(i)}-M_{NT,1}\right)^2. \tag*{(7)} M N T , 1 = N 1 i = 1 ∑ N Y N T , i π 1 ( i ) , S N T , 1 2 = N − 1 1 i = 1 ∑ N ( Y N T , i π 1 ( i ) − M N T , 1 ) 2 . (7)
For q ≥ 2 q \ge2 q ≥ 2 and fixed K < ∞ K < \infty K < ∞ , let Q N T , q ( K ) Q_{NT,q}(K) Q N T , q ( K ) denote the array laws satisfying Assumption 1 and the marginal moment bound E ∣ Y N T , 11 − θ N T ∣ q ≤ K σ N T q E|Y_{NT,11}-\theta_{NT}|^q \le K\sigma_{NT}^q E ∣ Y N T , 11 − θ N T ∣ q ≤ K σ N T q . Let Φ \Phi Φ be the standard normal distribution function and z p = Φ − 1 ( p ) z_p = \Phi^{-1}(p) z p = Φ − 1 ( p ) . We set a studentized statistic to zero whenever its denominator is zero.
Proposition 1 (One matching). As N , T → ∞ N,T \to\infty N , T → ∞ with N ≤ T N \le T N ≤ T ,
sup P ∈ Q N T , 3 ( K ) sup x ∈ R ∣ Pr P ( N ( M N T , 1 − θ N T ) S N T , 1 ≤ x ) − Φ ( x ) ∣ → 0. \sup_{P\in Q_{NT,3}(K)}\sup_{x\in\mathbb{R}}\left|\Pr_P\left(\frac{\sqrt{N}(M_{NT,1}-\theta_{NT})}{S_{NT,1}}\le x\right)-\Phi(x)\right|\to0. P ∈ Q N T , 3 ( K ) sup x ∈ R sup P Pr ( S N T , 1 N ( M N T , 1 − θ N T ) ≤ x ) − Φ ( x ) → 0.
By Proposition 1, the coverage probability of the one-matching interval M N T , 1 ± z 1 − α / 2 S N T , 1 / N M_{NT,1}\pm z_{1-\alpha/2}S_{NT,1}/\sqrt{N} M N T , 1 ± z 1 − α /2 S N T , 1 / N converges to 1 − α 1-\alpha 1 − α uniformly over the stated class of laws, whether or not the projection share ν N T \nu_{NT} ν N T converges. To assess its precision, write V N T , 1 = Var ( M N T , 1 ) = σ N T 2 / N V_{NT,1}=\operatorname{Var}(M_{NT,1})=\sigma_{NT}^2/N V N T , 1 = Var ( M N T , 1 ) = σ N T 2 / N . Comparing this variance with (3) gives
V N T , 1 − V N T , F = T − N N T σ N T , C 2 + T − 1 N T σ N T , I 2 ≥ 0. (8) V_{NT,1}-V_{NT,F}=\frac{T-N}{NT}\sigma_{NT,C}^2+\frac{T-1}{NT}\sigma_{NT,I}^2\ge0. \tag*{(8)} V N T , 1 − V N T , F = N T T − N σ N T , C 2 + N T T − 1 σ N T , I 2 ≥ 0. (8)
The row variance cancels because one matching uses every row. On a square array (N = T N=T N = T ) it also uses every column, so V N T , 1 = V N T , F V_{NT,1}=V_{NT,F} V N T , 1 = V N T , F exactly when σ N T , I 2 = 0 \sigma_{NT,I}^2=0 σ N T , I 2 = 0 . If N < T N<T N < T , equality requires both σ N T , I 2 = 0 \sigma^{2}_{NT,I}=0 σ N T , I 2 = 0 and σ N T , C 2 = 0 \sigma^{2}_{NT,C}=0 σ N T , C 2 = 0 .
The precision loss is particularly clear when the intersections are independent. One matching then averages N N N observations, while the full-array mean averages all N T NT N T , giving V N T , 1 = T V N T , F V_{NT,1}=T V_{NT,F} V N T , 1 = T V N T , F . A second matching can use observations omitted by the first. We next examine how much precision averaging several matching means can recover.
Repeated matching and trade-off
Fix the number of matchings B = B N T ≥ 1 B=B_{NT}\ge1 B = B N T ≥ 1 in advance. Draw B B B independent matchings by the rule in Section 3.1. Write π b ( i ) \pi_b(i) π b ( i ) for the column assigned to row i i i in matching b b b . The matching means and their average are
M N T , b = 1 N ∑ i = 1 N Y N T , i π b ( i ) , θ ^ N T , B = 1 B ∑ b = 1 B M N T , b . M_{NT,b}=\frac{1}{N}\sum_{i=1}^{N}Y_{NT,i\pi_b(i)},\qquad\hat{\theta}_{NT,B}=\frac{1}{B}\sum_{b=1}^{B}M_{NT,b}. M N T , b = N 1 i = 1 ∑ N Y N T , i π b ( i ) , θ ^ N T , B = B 1 b = 1 ∑ B M N T , b .
Conditional on the observed array Y N T Y_{NT} Y N T , the matching means are i.i.d. For each row, every column has probability 1 / T 1/T 1/ T of being selected, so E [ M N T , b ∣ Y N T ] = Y ˉ N T E[M_{NT,b}\mid Y_{NT}]=\bar{Y}_{NT} E [ M N T , b ∣ Y N T ] = Y ˉ N T . The error of the matching estimator therefore separates into two terms,
θ ^ N T , B − θ N T = ( Y ˉ N T − θ N T ) + 1 B ∑ b = 1 B ( M N T , b − Y ˉ N T ) . \hat{\theta}_{NT,B}-\theta_{NT}=(\bar{Y}_{NT}-\theta_{NT})+\frac{1}{B}\sum_{b=1}^{B}(M_{NT,b}-\bar{Y}_{NT}). θ ^ N T , B − θ N T = ( Y ˉ N T − θ N T ) + B 1 b = 1 ∑ B ( M N T , b − Y ˉ N T ) .
The first term is the sampling error of the full-array mean. It has variance V N T , F V_{NT,F} V N T , F from (3) and can have a non-Gaussian limit. The second term is the error from drawing the matchings. Conditional on Y N T Y_{NT} Y N T , its summands are independent and centered, so averaging B B B draws divides the variance of one draw by B B B . Under the conditions below, this error contributes a Gaussian term whenever it enters the limit.
Proposition 2 (Variance of repeated matching). For every N ≤ T N\le T N ≤ T and B ≥ 1 B\ge1 B ≥ 1 , E [ θ ^ N T , B ] = θ N T E[\hat{\theta}_{NT,B}]=\theta_{NT} E [ θ ^ N T , B ] = θ N T and
V N T , B : = Var ( θ ^ N T , B ) = V N T , F + V N T , 1 − V N T , F B . (9) V_{NT,B}:=\operatorname{Var}(\hat{\theta}_{NT,B})=V_{NT,F}+\frac{V_{NT,1}-V_{NT,F}}{B}. \tag*{(9)} V N T , B := Var ( θ ^ N T , B ) = V N T , F + B V N T , 1 − V N T , F . (9)
Proposition 2 holds for any N ≤ T N \le T N ≤ T . For the rest of this subsection, we use square array (N = T N = T N = T ) to show the trade-off in its simplest form. Every matching then uses every row and every column exactly once. Only the average of the selected intersection residuals varies across matchings. By (3) , (8) , and (9) ,
V N T , B = σ N T , R 2 + σ N T , C 2 N + σ N T , I 2 N 2 ⏟ full-array variance + σ N T , I 2 N B ( 1 − 1 N ) ⏟ matching variance . (10) V_{NT,B} =
\underbrace{\frac{\sigma^{2}_{NT,R}+\sigma^{2}_{NT,C}}{N}+\frac{\sigma^{2}_{NT,I}}{N^{2}}}_{\text{full-array variance}}
+
\underbrace{\frac{\sigma^{2}_{NT,I}}{NB}\left(1-\frac{1}{N}\right)}_{\text{matching variance}}. \tag*{(10)} V N T , B = full-array variance N σ N T , R 2 + σ N T , C 2 + N 2 σ N T , I 2 + matching variance N B σ N T , I 2 ( 1 − N 1 ) . (10)
The matching variance is proportional to 1 / B 1/B 1/ B , while the full-array variance is the same for every B B B . By (10) and (4) , the full-array share of the total variance is
w N T , B = V N T , F V N T , B = B B + ( N − 1 ) ( 1 − ν N T ) . (11) w_{NT,B}=\frac{V_{NT,F}}{V_{NT,B}}=\frac{B}{B+(N-1)(1-\nu_{NT})}. \tag*{(11)} w N T , B = V N T , B V N T , F = B + ( N − 1 ) ( 1 − ν N T ) B . (11)
Suppose the marginal fourth-moment bound of Proposition 3 holds. In the degenerate and intermediate cases (ν N T → c < 1 \nu_{NT} \to c < 1 ν N T → c < 1 ), choosing B = o ( N ) B = o(N) B = o ( N ) makes the full-array variance share w N T , B w_{NT,B} w N T , B vanish, so the standardized matching estimator has a standard normal limit. If B / N → λ ∈ ( 0 , ∞ ] B/N \to\lambda\in(0,\infty] B / N → λ ∈ ( 0 , ∞ ] , the full-array error retains a positive share of the variance and can contribute a non-Gaussian component to the limit. In the nondegenerate case (ν N T → 1 \nu_{NT} \to1 ν N T → 1 ), the standardized full-array mean is itself asymptotically Gaussian under the projection Lindeberg condition, so the standardized matching estimator has a standard normal limit for any deterministic choice of B B B .
The choice of B B B also determines precision. In the degenerate and intermediate cases, (10) gives
V N T , B ≍ σ N T 2 ( 1 N B + 1 N 2 ) . V_{NT,B}\asymp\sigma^{2}_{NT}\left(\frac{1}{NB}+\frac{1}{N^{2}}\right). V N T , B ≍ σ N T 2 ( N B 1 + N 2 1 ) .
When σ N T → σ ∈ ( 0 , ∞ ) \sigma_{NT} \to\sigma\in(0,\infty) σ N T → σ ∈ ( 0 , ∞ ) , the convergence rate is N B \sqrt{NB} N B for B = o ( N ) B = o(N) B = o ( N ) . Letting B → ∞ B \to\infty B → ∞ within this range therefore improves precision over one matching while preserving the Gaussian approximation. At B = N B = N B = N , the matching estimator attains the full-array rate N N N , although V N T , B / V N T , F → 2 − c V_{NT,B}/V_{NT,F} \to2-c V N T , B / V N T , F → 2 − c . Increasing B B B further can still reduce the variance: if B / N → ∞ B/N \to\infty B / N → ∞ , then V N T , B / V N T , F → 1 V_{NT,B}/V_{NT,F} \to1 V N T , B / V N T , F → 1 , with no further gain in the convergence rate.
Example 1: Pure product model. Consider the pure product model Y N T , i t = ξ i η t Y_{NT,it} = \xi_i\eta_t Y N T , i t = ξ i η t on a square array. Here ν N T = 0 \nu_{NT} = 0 ν N T = 0 , so w N T , B = B / ( B + N − 1 ) w_{NT,B} = B/(B + N - 1) w N T , B = B / ( B + N − 1 ) , and the standardized full-array mean has the limiting law of Z 1 Z 2 Z_1Z_2 Z 1 Z 2 . Holding the model fixed and varying only B B B gives
θ ^ N T , B − θ N T V N T , B ⇒ { Z 3 , B / N → 0 , ( Z 1 Z 2 + Z 3 ) / 2 , B = N , Z 1 Z 2 , B / N → ∞ . \frac{\widehat{\theta}_{NT,B} - \theta_{NT}}{\sqrt{V_{NT,B}}}
\Rightarrow
\begin{cases}
Z_3, & B/N \to0, \\
(Z_1Z_2 + Z_3)/\sqrt{2}, & B = N, \\
Z_1Z_2, & B/N \to\infty.
\end{cases} V N T , B θ N T , B − θ N T ⇒ ⎩ ⎨ ⎧ Z 3 , ( Z 1 Z 2 + Z 3 ) / 2 , Z 1 Z 2 , B / N → 0 , B = N , B / N → ∞.
where Z 1 , Z 2 , Z 3 Z_1, Z_2, Z_3 Z 1 , Z 2 , Z 3 are independent standard normals. The full-array error Y ‾ N T − θ N T \overline{Y}_{NT} - \theta_{NT} Y N T − θ N T accounts for none, one half, or all of the limiting variance, respectively.
The limiting distribution as B B B varies
In the square-array example, the share w N T , B w_{NT,B} w N T , B determines how much of the non-Gaussian full-array error survives. We now return to the general setting N ≤ T N \le T N ≤ T . When N < T N < T N < T , a matching uses every row but only N N N of the T T T columns, so the choice of columns can also contribute to the matching variance, as (8) shows. For any array dimensions, (9) gives the full-array variance share
w N T , B = V N T , F V N T , B = B B + V N T , 1 / V N T , F − 1 . (12) w_{NT,B} = \frac{V_{NT,F}}{V_{NT,B}} = \frac{B}{B + V_{NT,1}/V_{NT,F} - 1}. \tag*{(12)} w N T , B = V N T , B V N T , F = B + V N T , 1 / V N T , F − 1 B . (12)
Proposition 3 shows how this share weights the two errors in the limiting distribution.
Proposition 3 (Limiting law for any number of matchings). Let N , T → ∞ N, T \to\infty N , T → ∞ with N ≤ T N \le T N ≤ T . Suppose the array laws belong to Q N T , 4 ( K ) Q_{NT,4}(K) Q N T , 4 ( K ) for a common K K K , and let B = B N T ≥ 1 B = B_{NT} \ge1 B = B N T ≥ 1 be deterministic. Suppose w N T , B → w w_{NT,B} \to w w N T , B → w .
If w = 0 w = 0 w = 0 , then ( θ ^ N T , B − θ N T ) / V N T , B ⇒ N ( 0 , 1 ) (\widehat{\theta}_{NT,B} - \theta_{NT})/\sqrt{V_{NT,B}} \Rightarrow\mathcal{N}(0,1) ( θ N T , B − θ N T ) / V N T , B ⇒ N ( 0 , 1 ) .
If w > 0 w > 0 w > 0 and ( Y ‾ N T − θ N T ) / V N T , F ⇒ L F (\overline{Y}_{NT} - \theta_{NT})/\sqrt{V_{NT,F}} \Rightarrow L_F ( Y N T − θ N T ) / V N T , F ⇒ L F , then
θ ^ N T , B − θ N T V N T , B ⇒ w L F + 1 − w Z , \frac{\widehat{\theta}_{NT,B} - \theta_{NT}}{\sqrt{V_{NT,B}}} \Rightarrow\sqrt{w}L_F + \sqrt{1-w}Z, V N T , B θ N T , B − θ N T ⇒ w L F + 1 − w Z ,
where Z Z Z is standard normal and independent of L F L_F L F .
When w = 0 w = 0 w = 0 , the matching error alone determines the limit. For 0 < w < 1 0 < w < 1 0 < w < 1 , an independent Gaussian matching error enters alongside the full-array error. At w = 1 w = 1 w = 1 , the matching estimator inherits the full-array law. Thus, if the standardized full-array mean is asymptotically standard normal, so is the standardized matching estimator for any deterministic sequence B B B , by applying the proposition along subsequences. In the nondegenerate case (ν N T → 1 \nu_{NT} \to1 ν N T → 1 ), a Lindeberg condition on the row and column projection sums suffices for this conclusion.
The role of the longer dimension. In the degenerate and intermediate cases, it is T T T , the number of columns, that determines the transition. To see why, observe that
V N T , 1 T V N T , F = σ N T 2 N T V N T , F = 1 − ν N T + σ N T , R 2 + σ N T , C 2 N T V N T , F , 0 ≤ σ N T , R 2 + σ N T , C 2 N T V N T , F ≤ ν N T N . \frac{V_{NT,1}}{T V_{NT,F}} = \frac{\sigma^2_{NT}}{NTV_{NT,F}} = 1 - \nu_{NT} + \frac{\sigma^2_{NT,R} + \sigma^2_{NT,C}}{NTV_{NT,F}}, \qquad0 \le\frac{\sigma^2_{NT,R} + \sigma^2_{NT,C}}{NTV_{NT,F}} \le\frac{\nu_{NT}}{N}. T V N T , F V N T , 1 = N T V N T , F σ N T 2 = 1 − ν N T + N T V N T , F σ N T , R 2 + σ N T , C 2 , 0 ≤ N T V N T , F σ N T , R 2 + σ N T , C 2 ≤ N ν N T .
If ν N T → c < 1 \nu_{NT} \to c < 1 ν N T → c < 1 , it follows that ( V N T , 1 / V N T , F − 1 ) / T → 1 − c (V_{NT,1}/V_{NT,F} - 1)/T \to1 - c ( V N T , 1 / V N T , F − 1 ) / T → 1 − c . Hence B = o ( T ) B = o(T) B = o ( T ) gives w N T , B → 0 w_{NT,B} \to0 w N T , B → 0 , whereas B / T → ∞ B/T \to\infty B / T → ∞ gives w N T , B → 1 w_{NT,B} \to1 w N T , B → 1 and the matching estimator attains the variance of the full-array mean to first order. When B / T → λ ∈ ( 0 , ∞ ) B/T \to\lambda\in(0,\infty) B / T → λ ∈ ( 0 , ∞ ) , both errors contribute, with w N T , B → λ / ( λ + 1 − c ) w_{NT,B} \to\lambda/(\lambda+ 1 - c) w N T , B → λ / ( λ + 1 − c ) . Table 1 gives the corresponding rates and limits. For example, with B = N B = N B = N the pure product model has a non-Gaussian limit on a square array, but a Gaussian limit if T / N → ∞ T/N \to\infty T / N → ∞ : in the latter case N N N matchings still use only a vanishing fraction of the N T NT N T intersections.
Table 1. Rates and limiting laws of the matching estimator for N ≤ T N \le T N ≤ T .
Table 1 gives two routes to a Gaussian limit. When ν N T → c < 1 \nu_{NT} \to c < 1 ν N T → c < 1 , choosing B = o ( T ) B = o(T) B = o ( T ) makes the full-array error negligible. When ν N T → 1 \nu_{NT} \to1 ν N T → 1 , Table 1 gives a Gaussian limit for any B B B provided the projection sums satisfy the projection Lindeberg condition. The $T = N^2
Notes. The scale multiplies θ ^ N T , B − θ N T \hat{\theta}_{NT,B}-\theta_{NT} θ ^ N T , B − θ N T . The last two rows assume convergence of the standardized full-array mean to L F L_F L F as in Proposition 3; Z Z Z is standard normal and independent of L F L_F L F .
example in Appendix C illustrates the role of the Lindeberg condition in this second case.
There, ν N T → 1 \nu_{NT} \to1 ν N T → 1 , yet the standardized row projection sum converges to a centered Poisson law.
With B = N B=N B = N , the share of the matching estimator’s variance attributable to the row projection tends to one half. The standardized matching estimator therefore has a non-Gaussian limit even though B = o ( T ) B=o(T) B = o ( T ) .
Variance estimation and uniform Gaussian inference
We take B = o ( N ) B=o(N) B = o ( N ) . Under the moment bounds below, this choice yields asymptotically exact coverage uniformly over the admissible laws, without requiring ν N T \nu_{NT} ν N T to converge or imposing a separate Lindeberg condition on the full-array projection sums.
To accommodate sequences along which the full-array error remains important, we estimate the total variance V N T , B = ( 1 − 1 / B ) V N T , F + σ N T 2 / ( N B ) V_{NT,B}=(1-1/B)V_{NT,F}+\sigma^2_{NT}/(NB) V N T , B = ( 1 − 1/ B ) V N T , F + σ N T 2 / ( N B ) from Proposition 2. Let Y ˉ N T , i ⋅ \bar{Y}_{NT,i\cdot} Y ˉ N T , i ⋅ and Y ˉ N T , ⋅ t \bar{Y}_{NT,\cdot t} Y ˉ N T , ⋅ t denote the row and column means. The sample counterparts of the variance components in (3) are
σ ^ N T , R 2 = 1 N ∑ i ( Y ˉ N T , i ⋅ − Y ˉ N T ) 2 , σ ^ N T , C 2 = 1 T ∑ t ( Y ˉ N T , ⋅ t − Y ˉ N T ) 2 , r ^ N T , i t = Y N T , i t − Y ˉ N T , i ⋅ − Y ˉ N T , ⋅ t + Y ˉ N T , σ ^ N T , I 2 = 1 N T ∑ i , t r ^ N T , i t 2 . (13) \begin{aligned}
\hat{\sigma}^{2}_{NT,R} &= \frac{1}{N}\sum_i\left(\bar{Y}_{NT,i\cdot}-\bar{Y}_{NT}\right)^2, &
\hat{\sigma}^{2}_{NT,C} &= \frac{1}{T}\sum_t\left(\bar{Y}_{NT,\cdot t}-\bar{Y}_{NT}\right)^2,\\
\hat{r}_{NT,it} &= Y_{NT,it}-\bar{Y}_{NT,i\cdot}-\bar{Y}_{NT,\cdot t}+\bar{Y}_{NT}, &
\hat{\sigma}^{2}_{NT,I} &= \frac{1}{NT}\sum_{i,t}\hat{r}^{2}_{NT,it}. \tag*{(13)}
\end{aligned} σ ^ N T , R 2 r ^ N T , i t = N 1 i ∑ ( Y ˉ N T , i ⋅ − Y ˉ N T ) 2 , = Y N T , i t − Y ˉ N T , i ⋅ − Y ˉ N T , ⋅ t + Y ˉ N T , σ ^ N T , C 2 σ ^ N T , I 2 = T 1 t ∑ ( Y ˉ N T , ⋅ t − Y ˉ N T ) 2 , = N T 1 i , t ∑ r ^ N T , i t 2 . (13)
The marginal variance estimate σ ^ N T 2 = ( N T ) − 1 ∑ i , t ( Y N T , i t − Y ˉ N T ) 2 \hat{\sigma}^{2}_{NT}=(NT)^{-1}\sum_{i,t}(Y_{NT,it}-\bar{Y}_{NT})^2 σ ^ N T 2 = ( N T ) − 1 ∑ i , t ( Y N T , i t − Y ˉ N T ) 2 equals σ ^ N T , R 2 + σ ^ N T , C 2 + σ ^ N T , I 2 \hat{\sigma}^{2}_{NT,R}+\hat{\sigma}^{2}_{NT,C}+\hat{\sigma}^{2}_{NT,I} σ ^ N T , R 2 + σ ^ N T , C 2 + σ ^ N T , I 2 .
Substitution in (9) gives variance estimator
V ^ N T , B = ( 1 − 1 B ) ( σ ^ N T , R 2 N + σ ^ N T , C 2 T + σ ^ N T , I 2 N T ) + σ ^ N T 2 N B . (14) \widehat{V}_{NT,B}=\left(1-\frac{1}{B}\right)\left(\frac{\widehat{\sigma}^{2}_{NT,R}}{N}+\frac{\widehat{\sigma}^{2}_{NT,C}}{T}+\frac{\widehat{\sigma}^{2}_{NT,I}}{NT}\right)+\frac{\widehat{\sigma}^{2}_{NT}}{NB}. \tag*{(14)} V N T , B = ( 1 − B 1 ) ( N σ N T , R 2 + T σ N T , C 2 + N T σ N T , I 2 ) + N B σ N T 2 . (14)
The estimator V ^ N T , B \widehat{V}_{NT,B} V N T , B is positive semidefinite by construction in every sample. The standard two-way cluster-robust variance estimator for full-array mean subtracts the intersection component and can therefore yield a negative variance estimate [3 , 7 ] .
The matching estimator uses the selected observations, while its standard error V ^ N T , B \sqrt{\widehat{V}_{NT,B}} V N T , B uses the full array. The resulting confidence interval is
I N T , B ( 1 − α ) = [ θ ^ N T , B − z 1 − α / 2 V ^ N T , B , θ ^ N T , B + z 1 − α / 2 V ^ N T , B ] . (15) I_{NT,B}(1-\alpha)=\left[\widehat{\theta}_{NT,B}-z_{1-\alpha/2}\sqrt{\widehat{V}_{NT,B}},\ \widehat{\theta}_{NT,B}+z_{1-\alpha/2}\sqrt{\widehat{V}_{NT,B}}\right]. \tag*{(15)} I N T , B ( 1 − α ) = [ θ N T , B − z 1 − α /2 V N T , B , θ N T , B + z 1 − α /2 V N T , B ] . (15)
For any deterministic B B B such that w N T , B → 0 w_{NT,B}\to0 w N T , B → 0 , the proof of Proposition 3 also gives V ^ N T , B / V N T , B → 1 \widehat{V}_{NT,B}/V_{NT,B}\to1 V N T , B / V N T , B → 1 in probability under the marginal fourth-moment bound alone. Studentization therefore preserves the Gaussian limit along those sequences. For uniform coverage we must also allow the full-array error to retain a positive variance share. Every matching includes every row, so an unusually large row projection can affect every matching mean. We control such projections through the following moment bounds.
Assumption 2 (Moment bounds). For a constant K K K common to all laws and array sizes under consideration , E [ a N T , 1 4 ] ≤ K σ N T 2 σ N T , R 2 E[a_{NT,1}^{4}]\le K\sigma_{NT}^{2}\sigma_{NT,R}^{2} E [ a N T , 1 4 ] ≤ K σ N T 2 σ N T , R 2 , and E [ b N T , 1 4 ] ≤ K σ N T 2 σ N T , C 2 E[b_{NT,1}^{4}]\le K\sigma_{NT}^{2}\sigma_{NT,C}^{2} E [ b N T , 1 4 ] ≤ K σ N T 2 σ N T , C 2 .
Let P N T ( K ) P_{NT}(K) P N T ( K ) denote the subclass of Q N T , 4 ( K ) Q_{NT,4}(K) Q N T , 4 ( K ) that also satisfies the row and column projection bounds in Assumption 2. The moment bounds allow either projection variance to approach zero.
Theorem 1 (Uniform Gaussian inference). As N , T → ∞ N,T\to\infty N , T → ∞ with N ≤ T N\le T N ≤ T , let B = B N T ≥ 1 B=B_{NT}\ge1 B = B N T ≥ 1 be
deterministic with B = o ( N ) B=o(N) B = o ( N ) . For every ϵ > 0 \epsilon>0 ϵ > 0 ,
sup P ∈ P N T ( K ) Pr P ( ∣ V ^ N T , B V N T , B − 1 ∣ > ϵ ) → 0 , (16) \sup_{P\in\mathcal{P}_{NT}(K)}\Pr_{P}\left(\left|\frac{\widehat{V}_{NT,B}}{V_{NT,B}}-1\right|>\epsilon\right)\to0, \tag*{(16)} P ∈ P N T ( K ) sup P Pr ( V N T , B V N T , B − 1 > ϵ ) → 0 , (16)
sup P ∈ P N T ( K ) sup x ∈ R ∣ Pr P ( θ ^ N T , B − θ N T V ^ N T , B ≤ x ) − Φ ( x ) ∣ → 0. (17) \sup_{P\in\mathcal{P}_{NT}(K)}\sup_{x\in\mathbb{R}}\left|\Pr_{P}\left(\frac{\widehat{\theta}_{NT,B}-\theta_{NT}}{\sqrt{\widehat{V}_{NT,B}}}\leq x\right)-\Phi(x)\right|\to0. \tag*{(17)} P ∈ P N T ( K ) sup x ∈ R sup P Pr V N T , B θ N T , B − θ N T ≤ x − Φ ( x ) → 0. (17)
If B = o ( N ) B=o(\sqrt{N}) B = o ( N ) , (16) and (17) also hold uniformly over Q N T , 4 ( K ) \mathcal{Q}_{NT,4}(K) Q N T , 4 ( K ) without the projection bounds in Assumption 2.
For B = o ( N ) B=o(N) B = o ( N ) , the interval in (15) has uniformly asymptotically exact coverage over P N T ( K ) \mathcal{P}_{NT}(K) P N T ( K ) . The projection share ν N T \nu_{NT} ν N T may vary without converging. The theorem also requires no separate Lindeberg condition of the kind used in Table 1 to obtain a Gaussian limit for unrestricted B B B when ν N T → 1 \nu_{NT}\to1 ν N T → 1 .1
Remark 1. Theorem 1 permits any deterministic sequence B = o ( N ) B=o(N) B = o ( N ) . On a square array, one can therefore choose B B B so that the matching estimator attains the convergence rate of the full-array mean up to an arbitrarily slowly diverging factor. For example, under exact degeneracy, B = ⌊ N / log N ⌋ B=\lfloor N/\log N\rfloor B = ⌊ N / log N ⌋ gives a standard error asymptotic to σ N T log N / N \sigma_{NT}\sqrt{\log N}/N σ N T log N / N , compared with σ N T / N \sigma_{NT}/N σ N T / N for the full-array mean. The main cost of uniform inference in this case is therefore an additional factor of ( log N ) − 1 / 2 (\log N)^{-1/2} ( log N ) − 1/2 in the convergence rate.
The weaker moment condition in Theorem 1 places a separate restriction on B B B . With only a bounded standardized marginal fourth moment, an influential row can be rare enough that its sample count has a Poisson limit. Appendix C constructs a counterexample for every sequence with B = o ( N ) B=o(N) B = o ( N ) for which B / N B/\sqrt{N} B / N fails to converge to zero. Thus B = o ( N ) B=o(\sqrt{N}) B = o ( N ) is sharp for the uniform Gaussian approximation under a marginal fourth-moment bound alone.
Extension to generalized method of moments
The population mean solves the moment condition E ( Y N T , i t − θ N T ) = 0 E(Y_{NT,it}-\theta_{NT})=0 E ( Y N T , i t − θ N T ) = 0 . GMM replaces this one equation by a vector of equations. At the true parameter, each row-column intersection contributes a vector of moments. The components of this vector may exhibit different dependence regimes. We show that the Gaussian inference developed in § 3 extends to smooth functions of the GMM parameter even when the sample moments converge at different rates.
Let X N T , i t = f N T ( ξ i , η t , ζ i t ) \mathbf{X}_{NT,it}=f_{NT}(\xi_i,\eta_t,\zeta_{it}) X N T , i t = f N T ( ξ i , η t , ζ i t ) follow the sampling model in (1). Let m N T ( X N T , i t , β ) \mathbf{m}_{NT}(\mathbf{X}_{NT,it},\beta) m N T ( X N T , i t , β ) be a vector of r r r moment functions. Suppose a d d d -dimensional parameter β 0 , N T \beta_{0,NT} β 0 , N T satisfies E [ m N T ( X N T , i t , β 0 , N T ) ] = 0 E[\mathbf{m}_{NT}(\mathbf{X}_{NT,it},\beta_{0,NT})]=0 E [ m N T ( X N T , i t , β 0 , N T )] = 0 , where r ≥ d r\ge d r ≥ d and both dimensions are fixed. We consider inference on a scalar quantity ℓ ( β 0 , N T ) \ell(\beta_{0,NT}) ℓ ( β 0 , N T ) . In the GMM framework of Hansen [15 ] , the population moments identify β 0 , N T \beta_{0,NT} β 0 , N T . We replace their usual sample average by the matching average. Using the same B B B matchings at every candidate value β \beta β , compute
m ^ N T , B ( β ) = 1 N B ∑ b = 1 B ∑ i = 1 N m N T ( X N T , i π b ( i ) , β ) , \widehat{\mathbf{m}}_{NT,B}(\beta)=\frac{1}{NB}\sum_{b=1}^{B}\sum_{i=1}^{N}\mathbf{m}_{NT}(\mathbf{X}_{NT,i\pi_b(i)},\beta), m N T , B ( β ) = N B 1 b = 1 ∑ B i = 1 ∑ N m N T ( X N T , i π b ( i ) , β ) ,
β ^ N T , B ∈ arg min β ∈ Θ m ^ N T , B ( β ) T W ^ N T m ^ N T , B ( β ) . (18) \widehat{\beta}_{NT,B}\in\mathop{\arg\min}_{\beta\in\Theta}\widehat{\mathbf{m}}_{NT,B}(\beta)^{\mathsf{T}}\widehat{W}_{NT}\widehat{\mathbf{m}}_{NT,B}(\beta). \tag*{(18)} β N T , B ∈ arg min β ∈ Θ m N T , B ( β ) T W N T m N T , B ( β ) . (18)
The weight W ^ N T \widehat{W}_{NT} W N T is held fixed in this minimization. Exactly identified and overidentified models use the same matching average.
Let m ˙ N T \dot{\mathbf{m}}_{NT} m ˙ N T be the r × d r\times d r × d derivative of the moment vector and define
D N T = E m ˙ N T ( X N T , 11 , β 0 , N T ) , A N T = ( D N T T W N T D N T ) − 1 D N T T W N T , D_{NT}=E\dot{\mathbf{m}}_{NT}(\mathbf{X}_{NT,11},\beta_{0,NT}),\qquad A_{NT}=(D_{NT}^{\mathsf{T}}W_{NT}D_{NT})^{-1}D_{NT}^{\mathsf{T}}W_{NT}, D N T = E m ˙ N T ( X N T , 11 , β 0 , N T ) , A N T = ( D N T T W N T D N T ) − 1 D N T T W N T ,
where W N T W_{NT} W N T is the population weight. Linearizing the moment conditions gives the contribution of intersection ( i , t ) (i,t) ( i , t ) to the estimation error in ℓ ( β ^ N T , B ) \ell(\widehat{\beta}_{NT,B}) ℓ ( β N T , B ) ,
u N T , i t = − ∇ ℓ ( β 0 , N T ) T A N T m N T ( X N T , i t , β 0 , N T ) . (19) u_{NT,it}=-\nabla\ell(\beta_{0,NT})^{\mathsf{T}}A_{NT}\mathbf{m}_{NT}(\mathbf{X}_{NT,it},\beta_{0,NT}). \tag*{(19)} u N T , i t = − ∇ ℓ ( β 0 , N T ) T A N T m N T ( X N T , i t , β 0 , N T ) . (19)
For the population mean, this contribution equals Y N T , i t − θ N T Y_{NT,it}-\theta_{NT} Y N T , i t − θ N T . For a general GMM parameter, we apply the variance estimator in (14) to the scalar contributions in (19) . To obtain their fitted values u ^ N T , i t \widehat{u}_{NT,it} u N T , i t , use β ^ N T , B \widehat{\beta}_{NT,B} β N T , B , the matching average of m ˙ N T \dot{m}_{NT} m ˙ N T evaluated at β ^ N T , B \widehat{\beta}_{NT,B} β N T , B , and W ^ N T \widehat{W}_{NT} W N T in place of β 0 , N T \beta_{0,NT} β 0 , N T , D N T D_{NT} D N T , and W N T W_{NT} W N T , respectively. Denote the fitted derivative by D ^ N T , B \widehat{D}_{NT,B} D N T , B . Writing V ^ N T , B ( x ) \widehat{V}_{NT,B}(x) V N T , B ( x ) for (14) applied to a scalar array x x x , the resulting standard error is V ^ N T , B ( u ^ N T ) \sqrt{\widehat{V}_{NT,B}(\widehat{u}_{NT})} V N T , B ( u N T ) .
Assumption 3 (GMM regularity). The following conditions hold uniformly over array sizes and laws.
The parameter space Θ \Theta Θ is fixed, compact, and convex. The distance from β 0 , N T \beta_{0,NT} β 0 , N T to the boundary of Θ \Theta Θ and the smallest singular value of D N T D_{NT} D N T are both bounded away from zero. For every ϵ > 0 \epsilon> 0 ϵ > 0 , there is a c ϵ > 0 c_{\epsilon} > 0 c ϵ > 0 such that ∥ E m N T ( X N T , 11 , β ) ∥ ≥ c ϵ \lVert E m_{NT}(X_{NT,11},\beta)\rVert\ge c_{\epsilon} ∥ E m N T ( X N T , 11 , β )∥ ≥ c ϵ for all β ∈ Θ \beta\in\Theta β ∈ Θ with ∥ β − β 0 , N T ∥ ≥ ϵ \lVert\beta-\beta_{0,NT}\rVert\ge\epsilon ∥ β − β 0 , N T ∥ ≥ ϵ .
For every x x x , each m N T , j ( x , ⋅ ) m_{NT,j}(x,\cdot) m N T , j ( x , ⋅ ) , j = 1 , … , r j=1,\ldots,r j = 1 , … , r , is twice continuously differentiable on a neighborhood of Θ \Theta Θ . On Θ \Theta Θ , the moment functions and their first two derivatives are bounded in norm by F N T ( x ) F_{NT}(x) F N T ( x ) , with E [ F N T ( X N T , 11 ) 2 ] ≤ K E[F_{NT}(X_{NT,11})^2] \le K E [ F N T ( X N T , 11 ) 2 ] ≤ K . The gradient and Hessian of ℓ \ell ℓ are continuous and bounded on Θ \Theta Θ .
The weight W ^ N T \widehat{W}_{NT} W N T is symmetric positive definite and satisfies ∥ W ^ N T − W N T ∥ = O P ( N − 1 / 2 ) \lVert\widehat{W}_{NT}-W_{NT}\rVert=O_P(N^{-1/2}) ∥ W N T − W N T ∥ = O P ( N − 1/2 ) , where W N T W_{NT} W N T is deterministic and symmetric with eigenvalues bounded above and away from zero.
Assumption 3 ensures that the population moments uniformly identify β 0 , N T \beta_{0,NT} β 0 , N T and that their Jacobian has full column rank at β 0 , N T \beta_{0,NT} β 0 , N T . The smoothness and moment bounds justify linearizing the GMM estimator and its scalar target, while the weighting condition permits either a fixed weight or a weight estimated at the N \sqrt{N} N rate.
For q ≥ 2 q \ge2 q ≥ 2 and fixed constants K K K and v > 0 v > 0 v > 0 , let Q N T , G M M , q ( K ) Q_{NT,\mathrm{GMM},q}(K) Q N T , GMM , q ( K ) denote the class of GMM laws P P P satisfying Assumption 3 such that Var P ( u N T , 11 ) ≥ v \operatorname{Var}_P(u_{NT,11}) \ge v Var P ( u N T , 11 ) ≥ v and the law of the influence array u N T u_{NT} u N T belongs to Q N T , q ( K ) Q_{NT,q}(K) Q N T , q ( K ) . Let P N T , G M M ( K ) \mathcal{P}_{NT,\mathrm{GMM}}(K) P N T , GMM ( K ) be the subclass of Q N T , G M M , 4 ( K ) Q_{NT,\mathrm{GMM},4}(K) Q N T , GMM , 4 ( K ) for which the row and column projections of u N T u_{NT} u N T also satisfy Assumption 2.
Theorem 2 (Matching GMM). Let N , T → ∞ N,T \to\infty N , T → ∞ with N ≤ T N \le T N ≤ T , and let B = o ( N ) B=o(N) B = o ( N ) be deterministic. Then
sup P ∈ P N T , G M M ( K ) sup x ∈ R ∣ Pr P ( ℓ ( β ^ N T , B ) − ℓ ( β 0 , N T ) V ^ N T , B ( u ^ N T ) ≤ x ) − Φ ( x ) ∣ → 0. (20) \sup_{P\in\mathcal{P}_{NT,\mathrm{GMM}}(K)}\sup_{x\in\mathbb{R}}\left|\Pr_{P}\left(\frac{\ell(\widehat{\beta}_{NT,B})-\ell(\beta_{0,NT})}{\sqrt{\widehat{V}_{NT,B}(\widehat{u}_{NT})}}\le x\right)-\Phi(x)\right|\to0. \tag*{(20)} P ∈ P N T , GMM ( K ) sup x ∈ R sup P Pr V N T , B ( u N T ) ℓ ( β N T , B ) − ℓ ( β 0 , N T ) ≤ x − Φ ( x ) → 0. (20) If B = o ( N ) B=o(\sqrt{N}) B = o ( N ) , the convergence in (20) holds uniformly over the larger class Q N T , G M M , 4 ( K ) \mathcal{Q}_{NT,\mathrm{GMM},4}(K) Q N T , GMM , 4 ( K ) .
The lower variance bound v v v concerns variation in the contribution from a single intersection u N T , 11 u_{NT,11} u N T , 11 . Both its projection variances may equal zero. Appendix D also shows that, when the stated bounds hold for all unit linear contrasts, the joint Wald statistic for β 0 , N T \beta_{0,NT} β 0 , N T has a uniform χ d 2 \chi^2_d χ d 2 limit.
Remark 2. For nonlinear GMM, a quadratic term in the estimation error can enter the limiting distribution under two-way clustering [13 ] . In the present construction the estimation remainder is O P ( N − 1 ) O_P(N^{-1}) O P ( N − 1 ) , whereas the standard deviation of the matching average of u N T u_{NT} u N T is at least Var ( u N T , 11 ) / ( N B ) \sqrt{\operatorname{Var}(u_{NT,11})/(NB)} Var ( u N T , 11 ) / ( N B ) . Thus B = o ( N ) B=o(N) B = o ( N ) makes the remainder negligible even when the quantity of interest converges faster than other components of β ^ N T , B \widehat{\beta}_{NT,B} β N T , B .
Example 2: GMM moments with different dependence regimes. On a square array, observe X N T , i t = ( x 1 , i , x 2 , i t ) ⊤ X_{NT,it}=(x_{1,i},x_{2,it})^\top X N T , i t = ( x 1 , i , x 2 , i t ) ⊤ , where { x 1 , i } i = 1 N \{x_{1,i}\}_{i=1}^{N} { x 1 , i } i = 1 N and { x 2 , i t } i , t = 1 N \{x_{2,it}\}_{i,t=1}^{N} { x 2 , i t } i , t = 1 N are mutually independent standard normal draws. Define the moment function by m N T ( X N T , i t ; β ) = ( x 1 , i − β 1 x 2 , i t + β 1 2 − β 2 ) \boldsymbol{m}_{NT}(X_{NT,it};\beta)=\begin{pmatrix} x_{1,i}-\beta_1 \\ x_{2,it}+\beta_1^2-\beta_2 \end{pmatrix} m N T ( X N T , i t ; β ) = ( x 1 , i − β 1 x 2 , i t + β 1 2 − β 2 ) .
The first population moment gives β 1 = E [ x 1 , i ] = 0 \beta_1=E[x_{1,i}]=0 β 1 = E [ x 1 , i ] = 0 . Substituting into the second gives β 2 = E [ x 2 , i t ] + β 1 2 = 0 \beta_2=E[x_{2,it}]+\beta_1^2=0 β 2 = E [ x 2 , i t ] + β 1 2 = 0 . Thus β 0 , N T = 0 \beta_{0,NT}=0 β 0 , N T = 0 , and we take ℓ ( β ) = β 2 \ell(\beta)=\beta_2 ℓ ( β ) = β 2 .
Let x ˉ 1 = N − 1 ∑ i x 1 , i \bar{x}_1=N^{-1}\sum_i x_{1,i} x ˉ 1 = N − 1 ∑ i x 1 , i and x ˉ 2 = N − 2 ∑ i , t x 2 , i t \bar{x}_2=N^{-2}\sum_{i,t}x_{2,it} x ˉ 2 = N − 2 ∑ i , t x 2 , i t . The full-array moment equations give x ˉ 1 \bar{x}_1 x ˉ 1 for the first coefficient and x ˉ 2 + x ˉ 1 2 \bar{x}_2+\bar{x}_1^2 x ˉ 2 + x ˉ 1 2 for the second. Since N x ˉ 2 N\bar{x}_2 N x ˉ 2 and N x ˉ 1 \sqrt{N}\bar{x}_1 N x ˉ 1 are independent standard normals,
N β ^ N T , F , 2 = N x ˉ 2 + ( N x ˉ 1 ) 2 + o P ( 1 ) ⇒ Z 0 + Z 1 2 , N\widehat{\beta}_{NT,F,2}=N\bar{x}_2+(\sqrt{N}\bar{x}_1)^2+o_P(1)\Rightarrow Z_0+Z_1^2, N β N T , F , 2 = N x ˉ 2 + ( N x ˉ 1 ) 2 + o P ( 1 ) ⇒ Z 0 + Z 1 2 ,
where Z 0 Z_0 Z 0 and Z 1 Z_1 Z 1 are independent standard normals. This illustrates a key difficulty of the full-array estimator identified by Davezies et al. [13 ] Section 3): the estimation error in the slower first coefficient enters the limiting law of the faster second coefficient.
Every matching uses all N N N rows, so the first moment equation still gives x ˉ 1 \bar{x}_1 x ˉ 1 . In the second equation, x ˉ 2 \bar{x}_2 x ˉ 2 is replaced by x ˉ 2 , B = 1 N B ∑ b = 1 B ∑ i = 1 N x 2 , i π b ( i ) \bar{x}_{2,B} = \frac{1}{NB} \sum_{b=1}^{B} \sum_{i=1}^{N} x_{2,i}\pi_b(i) x ˉ 2 , B = N B 1 ∑ b = 1 B ∑ i = 1 N x 2 , i π b ( i ) . For B = o ( N ) B = o(N) B = o ( N ) , Table 1 gives N B x ˉ 2 , B ⇒ N ( 0 , 1 ) \sqrt{NB}\,\bar{x}_{2,B} \Rightarrow\mathcal{N}(0,1) N B x ˉ 2 , B ⇒ N ( 0 , 1 ) . The quadratic term is now negligible because B / N → 0 \sqrt{B/N} \to0 B / N → 0 . Hence
N B β ^ N T , B , 2 = N B x ˉ 2 , B + B N ( N x ˉ 1 ) 2 + o P ( 1 ) ⇒ N ( 0 , 1 ) . \sqrt{NB}\,\hat{\beta}_{NT,B,2} = \sqrt{NB}\,\bar{x}_{2,B} + \sqrt{\frac{B}{N}}\left(\sqrt{N}\,\bar{x}_1\right)^2 + o_P(1) \Rightarrow\mathcal{N}(0,1). N B β ^ N T , B , 2 = N B x ˉ 2 , B + N B ( N x ˉ 1 ) 2 + o P ( 1 ) ⇒ N ( 0 , 1 ) .
The Gaussian approximation remains valid even when different components of the GMM moment vector exhibit different dependence regimes.
Monte Carlo evidence for the mean
We return to the sample mean to measure how much repeated matching shortens an interval and how closely each method attains nominal coverage. Each design has N = T N = T N = T and θ N T = 0 \theta_{NT} = 0 θ N T = 0 . Let ξ i \xi_i ξ i , η t \eta_t η t , and ζ i t \zeta_{it} ζ i t be mutually independent standard normal shocks, independent across indices within each family. The four designs specify the outcome at each intersection as follows.
Independent intersections \quad Y N T , i t = ζ i t Y_{NT,it} = \zeta_{it} Y N T , i t = ζ i t , \qquad Additive effects \quad Y N T , i t = ξ i + η t Y_{NT,it} = \xi_i + \eta_t Y N T , i t = ξ i + η t ,
Pure product \quad Y N T , i t = ξ i η t Y_{NT,it} = \xi_i\eta_t Y N T , i t = ξ i η t , \qquad Local product \quad Y N T , i t = ξ i ( η t + N − 1 / 2 ) Y_{NT,it} = \xi_i(\eta_t + N^{-1/2}) Y N T , i t = ξ i ( η t + N − 1/2 ) .
Independent intersections and additive effects give two Gaussian benchmarks with quite different costs of matching. The two product designs have the same intersection variance, σ N T , I 2 = 1 \sigma^2_{NT,I} = 1 σ N T , I 2 = 1 , and zero column projection variance. Their row projection variances are σ N T , R 2 = 0 \sigma^2_{NT,R} = 0 σ N T , R 2 = 0 for the pure product and σ N T , R 2 = 1 / N \sigma^2_{NT,R} = 1/N σ N T , R 2 = 1/ N for the local product. The pure product is therefore exactly degenerate, whereas the local product is intermediate, with ν N T = 1 / 2 \nu_{NT} = 1/2 ν N T = 1/2 . They are the δ = 0 \delta= 0 δ = 0 and δ = 1 \delta= 1 δ = 1 cases of (5) , the local comparison in Menzel [24 ] .
For each array we construct five nominal 95% confidence intervals. For one matching we use the ordinary sample standard error of the N N N selected observations. For repeated matching we use (15) with B = ⌊ N 1 / 3 ⌋ B = \lfloor N^{1/3} \rfloor B = ⌊ N 1/3 ⌋ and B = ⌊ N / log N ⌋ B = \lfloor N/\log N \rfloor B = ⌊ N / log N ⌋ . Both choices satisfy B = o ( N ) B = o(N) B = o ( N ) , and the cube-root choice also satisfies B = o ( N ) B = o(\sqrt{N}) B = o ( N ) . All three matching intervals use the standard normal critical value and the same sequence of matching draws, taking the first B B B draws for each choice.
The other two intervals use the full-array mean. For the two-way cluster-robust interval of Cameron et al. [3 ] (CGM), the degrees-of-freedom-adjusted variance estimate on the square array is
V ^ C G M = σ ^ N T , R 2 + σ ^ N T , C 2 N − 1 − σ ^ N T 2 N 2 − 1 . \widehat{V}_{\mathrm{CGM}} = \frac{\widehat{\sigma}^{2}_{NT,R} + \widehat{\sigma}^{2}_{NT,C}}{N-1} - \frac{\widehat{\sigma}^{2}_{NT}}{N^{2}-1}. V CGM = N − 1 σ N T , R 2 + σ N T , C 2 − N 2 − 1 σ N T 2 .
The interval is Y ‾ N T ± z 0.975 max { 0 , V ^ C G M } \overline{Y}_{NT} \pm z_{0.975}\sqrt{\max\{0,\widehat{V}_{\mathrm{CGM}}\}} Y N T ± z 0.975 max { 0 , V CGM } . A nonpositive variance estimate yields a zero-length interval centered at Y ‾ N T \overline{Y}_{NT} Y N T . We include these replications in both coverage and mean length.
The bootstrap interval uses the model-selection procedure (BS-S) of Menzel [24 ] and its equal-tailed, nonstudentized interval (REG). Table 2 uses the selection threshold log ( N ) / N \log(N)/N log ( N ) / N for both corrected projection variances in every design.2
Table 2. Coverage and length of nominal 95% confidence intervals. Notes. One uses a single matching. The repeated matching uses the N 1 / 3 N^{1/3} N 1/3 rule and the N / log N N/\log N N / log N rule, respectively. CGM and BS-S use the full-array mean. Each design and size uses 10,000 Monte Carlo replications, with 999 bootstrap draws per replication for BS-S.
Coverage with one matching approaches the nominal 95% level across the sample sizes in Table 2 , in line with Proposition 1. Repeated matching recovers much of the precision lost by using only one matching. Under independence and in both product designs, the cube-root rule gives intervals about half as long as those from one matching at N = 100 N = 100 N = 100 , with coverage close to nominal at all three sample sizes. Increasing B B B to ⌊ N / log N ⌋ \lfloor N/\log N \rfloor ⌊ N / log N ⌋ shortens the intervals further but leads to greater overcoverage. Under additive effects, the two rules give identical intervals because every matching averages the same row and column effects.
Adding the local row effect changes coverage little for either repeated-matching rule, as the full-array error is asymptotically negligible in both product designs. BS-S gives shorter intervals in both designs and achieves close to nominal coverage under the pure product. Under the local product, its coverage still appears reasonably close to nominal in these samples, although it falls to about 92% at N = 100 N = 100 N = 100 . The row projection accounts for half the variance of the full-array mean, yet its variance 1 / N 1/N 1/ N lies below the BS-S selection threshold log ( N ) / N \log(N)/N log ( N ) / N . Consequently, BS-S omits the row effect with probability tending to one and continues to approximate the pure-product law, so its limiting coverage under the local product is below 95%. CGM undercovers more severely in both designs and frequently produces zero-length intervals.
Based on these results, we recommend B = ⌊ N 1 / 3 ⌋ B = \lfloor N^{1/3} \rfloor B = ⌊ N 1/3 ⌋ as a default. It delivers a substantial gain in precision over one matching while keeping finite-sample coverage close to nominal. When shorter intervals are the priority, B = ⌊ N / log N ⌋ B = \lfloor N/\log N \rfloor B = ⌊ N / log N ⌋ is a reasonable choice, with the trade-off of more conservative coverage in the independent and product designs.
Conclusion
A single matching produces an i.i.d. sample because every selected observation has a different row index and a different column index. We average matching means and use a closed-form standard error with the usual Gaussian critical value. Coverage is uniformly asymptotically exact even when the limiting law of the normalized full-array mean cannot be estimated uniformly. Repeated matching recovers much of the precision of the full array. The same variance formula yields inference for smooth functions of GMM parameters. With B = o ( N ) B = o(N) B = o ( N ) , a nonlinear remainder is negligible on the matching scale even when the moment components exhibit different dependence regimes. In such settings, the remainder can contribute to the limiting distribution of the full-array GMM estimator [13 ] .
One matching and the exact variance
Proof of Proposition 1. Condition on the random matching π 1 \pi_1 π 1 . By (6) , the selected observations use disjoint triples of shocks. Since the draw is independent of the array, for any z 1 , … , z N z_1,\ldots,z_N z 1 , … , z N ,
Pr ( ⋂ i = 1 N { Y N T , i , π 1 ( i ) ≤ z i } | π 1 ) = ∏ i = 1 N Pr ( Y N T , 11 ≤ z i ) . \Pr\left(\bigcap_{i=1}^{N}\{Y_{NT,i,\pi_1(i)} \le z_i\}\middle|\pi_1\right)=\prod_{i=1}^{N}\Pr(Y_{NT,11}\le z_i). Pr ( i = 1 ⋂ N { Y N T , i , π 1 ( i ) ≤ z i } π 1 ) = i = 1 ∏ N Pr ( Y N T , 11 ≤ z i ) .
The right-hand side is the same for every draw, so the product law also holds unconditionally. Write X N T , i = ( Y N T , i , π 1 ( i ) − θ N T ) / σ N T X_{NT,i}=(Y_{NT,i,\pi_1(i)}-\theta_{NT})/\sigma_{NT} X N T , i = ( Y N T , i , π 1 ( i ) − θ N T ) / σ N T . Along any sequence of laws in Q N T , 3 ( K ) \mathcal{Q}_{NT,3}(K) Q N T , 3 ( K ) , the X N T , i X_{NT,i} X N T , i within each sample are i.i.d. with mean zero, variance one, and E ∣ X N T , i ∣ 3 ≤ K E|X_{NT,i}|^3\le K E ∣ X N T , i ∣ 3 ≤ K . The Berry-Esseen inequality gives a standard normal approximation for N − 1 / 2 ∑ i X N T , i N^{-1/2}\sum_i X_{NT,i} N − 1/2 ∑ i X N T , i with error at most C K / N CK/\sqrt{N} C K / N . For sample variance consistency, truncate X N T , i 2 X_{NT,i}^2 X N T , i 2 at a fixed L 2 L^2 L 2 . The truncated sample average concentrates around its expectation by Chebyshev’s inequality, and E [ X N T , i 2 1 { ∣ X N T , i ∣ > L } ] ≤ K / L E[X_{NT,i}^2\mathbf{1}\{|X_{NT,i}|>L\}]\le K/L E [ X N T , i 2 1 { ∣ X N T , i ∣ > L }] ≤ K / L . Letting L → ∞ L\to\infty L → ∞ gives N − 1 ∑ i X N T , i 2 → 1 N^{-1}\sum_i X_{NT,i}^2\to1 N − 1 ∑ i X N T , i 2 → 1 in probability.
Since E [ X ˉ N T 2 ] = 1 / N E[\bar{X}_{NT}^{2}]=1/N E [ X ˉ N T 2 ] = 1/ N ,
S N T , 1 2 σ N T 2 = N N − 1 ( 1 N ∑ i X N T , i 2 − X ˉ N T 2 ) → P 1. \frac{S_{NT,1}^{2}}{\sigma_{NT}^{2}}=\frac{N}{N-1}\left(\frac{1}{N}\sum_i X_{NT,i}^{2}-\bar{X}_{NT}^{2}\right)\xrightarrow{\mathrm{P}}1. σ N T 2 S N T , 1 2 = N − 1 N ( N 1 i ∑ X N T , i 2 − X ˉ N T 2 ) P 1.
In particular, the probability of a zero denominator vanishes. Slutsky’s theorem and continuity of Φ \Phi Φ give convergence in Kolmogorov distance along every admissible sequence, proving uniformity. For a direct bound on the studentized statistic, see Bentkus and Götze [2 ] .
Proof of Proposition 2. Each column is equally likely to be assigned to row i i i , so E [ M N T , b ∣ Y N T ] = N − 1 ∑ i T − 1 ∑ t Y N T , i t = Y ˉ N T E[M_{NT,b}\mid Y_{NT}]=N^{-1}\sum_i T^{-1}\sum_t Y_{NT,it}=\bar{Y}_{NT} E [ M N T , b ∣ Y N T ] = N − 1 ∑ i T − 1 ∑ t Y N T , i t = Y ˉ N T . Thus E θ ^ N T , B = E Y ˉ N T = θ N T E\hat{\theta}_{NT,B}=E\bar{Y}_{NT}=\theta_{NT} E θ ^ N T , B = E Y ˉ N T = θ N T . Conditional independence of the B B B draws and the law of total variance give
Var ( θ ^ N T , B ) = Var ( Y ˉ N T ) + 1 B E Var ( M N T , 1 ∣ Y N T ) , \operatorname{Var}(\hat{\theta}_{NT,B})=\operatorname{Var}(\bar{Y}_{NT})+\frac{1}{B}E\operatorname{Var}(M_{NT,1}\mid Y_{NT}), Var ( θ ^ N T , B ) = Var ( Y ˉ N T ) + B 1 E Var ( M N T , 1 ∣ Y N T ) ,
E Var ( M N T , 1 ∣ Y N T ) = Var ( M N T , 1 ) − Var ( Y ˉ N T ) . E\operatorname{Var}(M_{NT,1}\mid Y_{NT})=\operatorname{Var}(M_{NT,1})-\operatorname{Var}(\bar{Y}_{NT}). E Var ( M N T , 1 ∣ Y N T ) = Var ( M N T , 1 ) − Var ( Y ˉ N T ) .
Substitution gives (9).
Proof of Theorem 1
Fix any sequence of laws in P N T ( K ) \mathcal{P}_{NT}(K) P N T ( K ) as N , T → ∞ N,T\to\infty N , T → ∞ with N ≤ T N\leq T N ≤ T . Both conclusions are invariant to centering each array at θ N T \theta_{NT} θ N T and dividing its outcomes by σ N T \sigma_{NT} σ N T . We work with the standardized array throughout this proof, so θ N T = 0 \theta_{NT}=0 θ N T = 0 and σ N T 2 = 1 \sigma_{NT}^{2}=1 σ N T 2 = 1 at each ( N , T ) (N,T) ( N , T ) . Constants in the bounds below depend only on K K K . From (9),
V N T , B ≥ V N T , F , V N T , B ≥ 1 N B , V N T , F = σ N T , R 2 N + σ N T , C 2 T + σ N T , I 2 N T . (21) V_{NT,B}\geq V_{NT,F},\qquad V_{NT,B}\geq\frac{1}{NB},\qquad V_{NT,F}=\frac{\sigma_{NT,R}^{2}}{N}+\frac{\sigma_{NT,C}^{2}}{T}+\frac{\sigma_{NT,I}^{2}}{NT}. \tag*{(21)} V N T , B ≥ V N T , F , V N T , B ≥ N B 1 , V N T , F = N σ N T , R 2 + T σ N T , C 2 + N T σ N T , I 2 . (21)
Write a ˉ N T \bar{a}_{NT} a ˉ N T , b ˉ N T \bar{b}_{NT} b ˉ N T and r ˉ N T \bar{r}_{NT} r ˉ N T for the sample averages of the row projections, column projections and residuals, respectively.
Variance estimation. Distinct residuals are uncorrelated, so their row, column and full-array averages satisfy
E r ˉ N T , i ⋅ 2 = σ N T , I 2 T , E r ˉ N T , ⋅ t 2 = σ N T , I 2 N , E r ˉ N T 2 = σ N T , I 2 N T . E\bar{r}_{NT,i\cdot}^{2} = \frac{\sigma_{NT,I}^{2}}{T}, \qquad E\bar{r}_{NT,\cdot t}^{2} = \frac{\sigma_{NT,I}^{2}}{N}, \qquad E\bar{r}_{NT}^{2} = \frac{\sigma_{NT,I}^{2}}{NT}. E r ˉ N T , i ⋅ 2 = T σ N T , I 2 , E r ˉ N T , ⋅ t 2 = N σ N T , I 2 , E r ˉ N T 2 = N T σ N T , I 2 .
To bound the row variance estimate, write its centered row means in terms of the true projections and residuals:
σ ^ N T , R 2 = 1 N ∑ i [ ( a N T , i − a ˉ N T ) + ( r ˉ N T , i ⋅ − r ˉ N T ) ] 2 . \widehat{\sigma}_{NT,R}^{2} = \frac{1}{N}\sum_i \left[(a_{NT,i}-\bar{a}_{NT})+(\bar{r}_{NT,i\cdot}-\bar{r}_{NT})\right]^2. σ N T , R 2 = N 1 i ∑ [ ( a N T , i − a ˉ N T ) + ( r ˉ N T , i ⋅ − r ˉ N T ) ] 2 .
The mean square of the residual contrasts is bounded by N − 1 ∑ i r ˉ N T , i ⋅ 2 = O P ( σ N T , I 2 / T ) N^{-1}\sum_i\bar{r}_{NT,i\cdot}^{2}=O_{\mathrm{P}}(\sigma_{NT,I}^{2}/T) N − 1 ∑ i r ˉ N T , i ⋅ 2 = O P ( σ N T , I 2 / T ) . The mean square of the projection contrasts is O P ( σ N T , R 2 ) O_{\mathrm{P}}(\sigma_{NT,R}^{2}) O P ( σ N T , R 2 ) , so Cauchy-Schwarz bounds the cross term by O P ( σ N T , R 2 σ N T , I 2 / T ) O_{\mathrm{P}}(\sqrt{\sigma_{NT,R}^{2}\sigma_{NT,I}^{2}/T}) O P ( σ N T , R 2 σ N T , I 2 / T ) . For the projection term, independence and Assumption 2 give
Var ( N − 1 ∑ i a N T , i 2 ) ≤ K σ N T , R 2 N , E a ˉ N T 2 = σ N T , R 2 N . \operatorname{Var}\left(N^{-1}\sum_i a_{NT,i}^{2}\right)\leq\frac{K\sigma_{NT,R}^{2}}{N}, \qquad E\bar{a}_{NT}^{2}=\frac{\sigma_{NT,R}^{2}}{N}. Var ( N − 1 i ∑ a N T , i 2 ) ≤ N K σ N T , R 2 , E a ˉ N T 2 = N σ N T , R 2 .
Applying the same argument to the columns gives the bounds
∣ σ ^ N T , R 2 − σ N T , R 2 ∣ = O P ( σ N T , R 2 N + σ N T , R 2 σ N T , I 2 T + σ N T , I 2 T ) , \left|\widehat{\sigma}_{NT,R}^{2}-\sigma_{NT,R}^{2}\right|=O_{\mathrm{P}}\left(\sqrt{\frac{\sigma_{NT,R}^{2}}{N}}+\sqrt{\frac{\sigma_{NT,R}^{2}\sigma_{NT,I}^{2}}{T}}+\frac{\sigma_{NT,I}^{2}}{T}\right), σ N T , R 2 − σ N T , R 2 = O P N σ N T , R 2 + T σ N T , R 2 σ N T , I 2 + T σ N T , I 2 ,
∣ σ ^ N T , C 2 − σ N T , C 2 ∣ = O P ( σ N T , C 2 T + σ N T , C 2 σ N T , I 2 N + σ N T , I 2 N ) . (22) \left|\widehat{\sigma}_{NT,C}^{2}-\sigma_{NT,C}^{2}\right|=O_{\mathrm{P}}\left(\sqrt{\frac{\sigma_{NT,C}^{2}}{T}}+\sqrt{\frac{\sigma_{NT,C}^{2}\sigma_{NT,I}^{2}}{N}}+\frac{\sigma_{NT,I}^{2}}{N}\right). \tag*{(22)} σ N T , C 2 − σ N T , C 2 = O P T σ N T , C 2 + N σ N T , C 2 σ N T , I 2 + N σ N T , I 2 . (22)
By Jensen’s inequality and the marginal fourth-moment bound, E r N T , 11 4 Er_{NT,11}^{4} E r N T , 11 4 is bounded uniformly. The squares r N T , i t 2 r_{NT,it}^{2} r N T , i t 2 and r N T , i ′ t ′ 2 r_{NT,i't'}^{2} r N T , i ′ t ′ 2 are independent when i ≠ i ′ i\ne i' i = i ′ and t ≠ t ′ t\ne t' t = t ′ . Counting pairs with a common index gives Var ( ( N T ) − 1 ∑ i , t r N T , i t 2 ) = O ( 1 / N + 1 / T ) \operatorname{Var}((NT)^{-1}\sum_{i,t}r_{NT,it}^{2})=O(1/N+1/T) Var (( N T ) − 1 ∑ i , t r N T , i t 2 ) = O ( 1/ N + 1/ T ) . Double demeaning the outcome has the same effect as double demeaning the residual. Thus
0 ≤ 1 N T ∑ i , t r N T , i t 2 − σ ^ N T , I 2 = 1 N ∑ i r ˉ N T , i ⋅ 2 + 1 T ∑ t r ˉ N T , ⋅ t 2 − r ˉ N T 2 = O P ( 1 / N ) . 0\leq\frac{1}{NT}\sum_{i,t}r_{NT,it}^{2}-\widehat{\sigma}_{NT,I}^{2}=\frac{1}{N}\sum_i\bar{r}_{NT,i\cdot}^{2}+\frac{1}{T}\sum_t\bar{r}_{NT,\cdot t}^{2}-\bar{r}_{NT}^{2}=O_{\mathrm{P}}(1/N). 0 ≤ N T 1 i , t ∑ r N T , i t 2 − σ N T , I 2 = N 1 i ∑ r ˉ N T , i ⋅ 2 + T 1 t ∑ r ˉ N T , ⋅ t 2 − r ˉ N T 2 = O P ( 1/ N ) .
Consequently, σ ^ N T , I 2 − σ N T , I 2 = O P ( N − 1 / 2 ) \widehat{\sigma}^{2}_{NT,I}-\sigma^{2}_{NT,I}=O_P(N^{-1/2}) σ N T , I 2 − σ N T , I 2 = O P ( N − 1/2 ) . The same pair count applies to Y N T , i t 2 Y_{NT,it}^{2} Y N T , i t 2 , and E [ Y ‾ N T 2 ] = V N T , F ≤ 1 / N E[\overline{Y}_{NT}^{2}]=V_{NT,F}\leq1/N E [ Y N T 2 ] = V N T , F ≤ 1/ N . Thus σ ^ N T 2 − 1 = O P ( N − 1 / 2 ) \widehat{\sigma}^{2}_{NT}-1=O_P(N^{-1/2}) σ N T 2 − 1 = O P ( N − 1/2 ) as well.
By (21) , N V N T , B ≥ max { σ N T , R 2 , 1 / B } NV_{NT,B}\geq\max\{\sigma^{2}_{NT,R},1/B\} N V N T , B ≥ max { σ N T , R 2 , 1/ B } and T V N T , B ≥ max { σ N T , C 2 , 1 / B } TV_{NT,B}\geq\max\{\sigma^{2}_{NT,C},1/B\} T V N T , B ≥ max { σ N T , C 2 , 1/ B } . In particular x / max { x , 1 / B } ≤ B \sqrt{x}/\max\{x,1/B\}\leq\sqrt{B} x / max { x , 1/ B } ≤ B for all x ≥ 0 x\geq0 x ≥ 0 . Dividing the first and second bounds in (22) by N V N T , B NV_{NT,B} N V N T , B and T V N T , B TV_{NT,B} T V N T , B , respectively, therefore gives
∣ σ ^ N T , R 2 − σ N T , R 2 ∣ N V N T , B + ∣ σ ^ N T , C 2 − σ N T , C 2 ∣ T V N T , B = O P ( B N + B N ) . \frac{\left|\widehat{\sigma}^{2}_{NT,R}-\sigma^{2}_{NT,R}\right|}{NV_{NT,B}}+\frac{\left|\widehat{\sigma}^{2}_{NT,C}-\sigma^{2}_{NT,C}\right|}{TV_{NT,B}}=O_P\left(\sqrt{\frac{B}{N}}+\frac{B}{N}\right). N V N T , B σ N T , R 2 − σ N T , R 2 + T V N T , B σ N T , C 2 − σ N T , C 2 = O P ( N B + N B ) .
The estimation errors for σ N T , I 2 \sigma^{2}_{NT,I} σ N T , I 2 and σ ^ N T 2 \widehat{\sigma}^{2}_{NT} σ N T 2 in (14), divided by V N T , B V_{NT,B} V N T , B , are o P ( 1 ) o_P(1) o P ( 1 ) since N T V N T , B ≥ T / B NTV_{NT,B}\geq T/B N T V N T , B ≥ T / B and N B V N T , B ≥ 1 NBV_{NT,B}\geq1 N B V N T , B ≥ 1 . This proves (16).
Variation in the full-array mean. The identity Y ‾ N T = a ‾ N T + b ‾ N T + r ‾ N T \overline{Y}_{NT}=\overline{a}_{NT}+\overline{b}_{NT}+\overline{r}_{NT} Y N T = a N T + b N T + r N T separates its independent projection sums from the residual average. By (21), E r ‾ N T 2 / V N T , B ≤ B / T → 0 E\overline{r}_{NT}^{2}/V_{NT,B}\leq B/T\to0 E r N T 2 / V N T , B ≤ B / T → 0 . Write V N T , P = σ N T , R 2 / N + σ N T , C 2 / T = ν N T V N T , F V_{NT,P}=\sigma^{2}_{NT,R}/N+\sigma^{2}_{NT,C}/T=\nu_{NT}V_{NT,F} V N T , P = σ N T , R 2 / N + σ N T , C 2 / T = ν N T V N T , F . The sum of fourth moments of the summands of ( a ‾ N T + b ‾ N T ) / V N T , B (\overline{a}_{NT}+\overline{b}_{NT})/\sqrt{V_{NT,B}} ( a N T + b N T ) / V N T , B satisfies
1 V N T , B 2 ( E a ‾ N T , 1 4 N 3 + E b ‾ N T , 1 4 T 3 ) ≤ K V N T , P N 2 V N T , B 2 ≤ K B N → 0. \frac{1}{V_{NT,B}^{2}}\left(\frac{E\overline{a}_{NT,1}^{4}}{N^{3}}+\frac{E\overline{b}_{NT,1}^{4}}{T^{3}}\right)\leq\frac{KV_{NT,P}}{N^{2}V_{NT,B}^{2}}\leq\frac{KB}{N}\to0. V N T , B 2 1 ( N 3 E a N T , 1 4 + T 3 E b N T , 1 4 ) ≤ N 2 V N T , B 2 K V N T , P ≤ N K B → 0.
From any subsequence take a further subsequence with V N T , P / V N T , B → ℓ ∈ [ 0 , 1 ] V_{NT,P}/V_{NT,B}\to\ell\in[0,1] V N T , P / V N T , B → ℓ ∈ [ 0 , 1 ] . If ℓ > 0 \ell>0 ℓ > 0 , Lyapunov’s theorem gives Y ‾ N T / V N T , B ⇒ N ( 0 , ℓ ) \overline{Y}_{NT}/\sqrt{V_{NT,B}}\Rightarrow N(0,\ell) Y N T / V N T , B ⇒ N ( 0 , ℓ ) . If ℓ = 0 \ell=0 ℓ = 0 , the normalized full-array mean converges to zero in mean square.
Variation from drawing the matchings. Put D N T , b = M N T , b − Y ‾ N T D_{NT,b}=M_{NT,b}-\overline{Y}_{NT} D N T , b = M N T , b − Y N T . Given the full array Y N T Y_{NT} Y N T , the D N T , b D_{NT,b} D N T , b are i.i.d. with conditional mean zero. The conditional variance has the observable form
ω N T 2 ( Y N T ) : = N Var ( D N T , 1 ∣ Y N T ) = T T − 1 ( σ ^ N T , I 2 + ( 1 − N / T ) σ ^ N T , C 2 ) . (23) \omega^{2}_{NT}(Y_{NT}):=N\operatorname{Var}(D_{NT,1}\mid Y_{NT})=\frac{T}{T-1}\left(\widehat{\sigma}^{2}_{NT,I}+(1-N/T)\widehat{\sigma}^{2}_{NT,C}\right). \tag*{(23)} ω N T 2 ( Y N T ) := N Var ( D N T , 1 ∣ Y N T ) = T − 1 T ( σ N T , I 2 + ( 1 − N / T ) σ N T , C 2 ) . (23)
To verify (23) , pad the N N N row-centered observations in each column with T − N T-N T − N zeros to obtain a T × T T \times T T × T matrix A A A . Write g t = Y ˉ N T , t − Y ˉ N T g_t=\bar{Y}_{NT,t}-\bar{Y}_{NT} g t = Y ˉ N T , t − Y ˉ N T . The mean of column t t t of A A A is ( N / T ) g t (N/T)g_t ( N / T ) g t .
Let C C C be A A A after subtracting these column means, so all its row and column sums are zero. Since
( N T ) − 1 ∑ i , t ≤ T A i t 2 = σ ^ N T , I 2 + σ ^ N T , C 2 (NT)^{-1}\sum_{i,t\leq T} A_{it}^{2}=\hat{\sigma}_{NT,I}^{2}+\hat{\sigma}_{NT,C}^{2} ( N T ) − 1 i , t ≤ T ∑ A i t 2 = σ ^ N T , I 2 + σ ^ N T , C 2
and
T − 1 ∑ t g t 2 = σ ^ N T , C 2 , T^{-1}\sum_{t}g_t^{2}=\hat{\sigma}_{NT,C}^{2}, T − 1 t ∑ g t 2 = σ ^ N T , C 2 ,
1 N T ∑ i , t ≤ T C i t 2 = σ ^ N T , I 2 + σ ^ N T , C 2 − 1 N ∑ t = 1 T ( N T g t ) 2 = σ ^ N T , I 2 + ( 1 − N / T ) σ ^ N T , C 2 . \frac{1}{NT}\sum_{i,t\leq T}C_{it}^{2}=\hat{\sigma}_{NT,I}^{2}+\hat{\sigma}_{NT,C}^{2}-\frac{1}{N}\sum_{t=1}^{T}\left(\frac{N}{T}g_t\right)^{2}=\hat{\sigma}_{NT,I}^{2}+\left(1-N/T\right)\hat{\sigma}_{NT,C}^{2}. N T 1 i , t ≤ T ∑ C i t 2 = σ ^ N T , I 2 + σ ^ N T , C 2 − N 1 t = 1 ∑ T ( T N g t ) 2 = σ ^ N T , I 2 + ( 1 − N / T ) σ ^ N T , C 2 .
Complete the matching to a uniform permutation of the T T T columns. Its sum over C C C equals N D N T , 1 ND_{NT,1} N D N T , 1 , since the subtracted column means sum to zero. Put S C = ∑ i , t ≤ T C i t 2 S_C=\sum_{i,t\leq T}C_{it}^{2} S C = ∑ i , t ≤ T C i t 2 . The zero row and column sums give ∑ j ≠ i , s ≠ t C j s = C i t \sum_{j\ne i,\ s\ne t}C_{js}=C_{it} ∑ j = i , s = t C j s = C i t . A single permutation assignment has probability 1 / T 1/T 1/ T , and a specified pair in distinct rows and columns has probability 1 / [ T ( T − 1 ) ] 1/[T(T-1)] 1/ [ T ( T − 1 )] . Hence
Var ( N D N T , 1 ∣ Y N T ) = S C T + ∑ i , t ≤ T C i t ∑ j ≠ i , s ≠ t C j s T ( T − 1 ) = S C T − 1 . \operatorname{Var}(ND_{NT,1}\mid Y_{NT})=\frac{S_C}{T}+\frac{\sum_{i,t\leq T}C_{it}\sum_{j\ne i,\ s\ne t}C_{js}}{T(T-1)}=\frac{S_C}{T-1}. Var ( N D N T , 1 ∣ Y N T ) = T S C + T ( T − 1 ) ∑ i , t ≤ T C i t ∑ j = i , s = t C j s = T − 1 S C .
which proves (23) . Total variance gives E ω N T 2 ( Y N T ) = 1 − N V N T , F = ( 1 − N / T ) σ N T , C 2 + ( 1 − 1 / T ) σ N T , I 2 E\omega_{NT}^{2}(Y_{NT})=1-NV_{NT,F}=(1-N/T)\sigma_{NT,C}^{2}+(1-1/T)\sigma_{NT,I}^{2} E ω N T 2 ( Y N T ) = 1 − N V N T , F = ( 1 − N / T ) σ N T , C 2 + ( 1 − 1/ T ) σ N T , I 2 . By (23) and variance component consistency, ω N T 2 ( Y N T ) − E ω N T 2 ( Y N T ) = o P ( 1 ) \omega_{NT}^{2}(Y_{NT})-E\omega_{NT}^{2}(Y_{NT})=o_P(1) ω N T 2 ( Y N T ) − E ω N T 2 ( Y N T ) = o P ( 1 ) . Consequently, along the chosen subsequence,
ω N T 2 ( Y N T ) N B V N T , B = 1 − V N T , F V N T , B + o P ( 1 ) ⟶ P ρ : = 1 − ℓ . \frac{\omega_{NT}^{2}(Y_{NT})}{NBV_{NT,B}}=1-\frac{V_{NT,F}}{V_{NT,B}}+o_P(1)\mathrel{\mathop{\longrightarrow}^{\mathrm{P}}}\rho:=1-\ell. N B V N T , B ω N T 2 ( Y N T ) = 1 − V N T , B V N T , F + o P ( 1 ) ⟶ P ρ := 1 − ℓ .
If ρ = 0 \rho=0 ρ = 0 , the normalized matching error ( θ ^ N T , B − Y ˉ N T ) / V N T , B (\hat{\theta}_{NT,B}-\bar{Y}_{NT})/\sqrt{V_{NT,B}} ( θ ^ N T , B − Y ˉ N T ) / V N T , B vanishes in mean square. If ρ > 0 \rho>0 ρ > 0 , (21) gives Pr ( ω N T 2 ( Y N T ) > ρ / 2 ) → 1 \Pr(\omega_{NT}^{2}(Y_{NT})>\rho/2)\to1 Pr ( ω N T 2 ( Y N T ) > ρ /2 ) → 1 .
Consider first B → ∞ B\to\infty B → ∞ with ρ > 0 \rho>0 ρ > 0 . Because the N N N observations in one matching are i.i.d.,
E [ ( N M N T , 1 ) 4 ] ≤ K / N + 3. E[(\sqrt{N}M_{NT,1})^{4}]\leq K/N+3. E [( N M N T , 1 ) 4 ] ≤ K / N + 3.
As Y ˉ N T = E [ M N T , 1 ∣ Y N T ] \bar{Y}_{NT}=E[M_{NT,1}\mid Y_{NT}] Y ˉ N T = E [ M N T , 1 ∣ Y N T ] , Jensen’s and Markov’s inequalities yield E [ N 2 D N T , 1 4 ∣ Y N T ] = O P ( 1 ) E[N^{2}D_{NT,1}^{4}\mid Y_{NT}]=O_P(1) E [ N 2 D N T , 1 4 ∣ Y N T ] = O P ( 1 ) . The conditional Lyapunov ratio is therefore
∑ b = 1 B E [ D N T , b 4 ∣ Y N T ] { ∑ b = 1 B Var ( D N T , b ∣ Y N T ) } 2 = E [ N 2 D N T , 1 4 ∣ Y N T ] B ( ω N T 2 ( Y N T ) ) 2 = o P ( 1 ) . \frac{\sum_{b=1}^{B}E[D_{NT,b}^{4}\mid Y_{NT}]}{\left\{\sum_{b=1}^{B}\operatorname{Var}(D_{NT,b}\mid Y_{NT})\right\}^{2}}
=\frac{E[N^{2}D_{NT,1}^{4}\mid Y_{NT}]}{B(\omega_{NT}^{2}(Y_{NT}))^{2}}=o_P(1). { ∑ b = 1 B Var ( D N T , b ∣ Y N T ) } 2 ∑ b = 1 B E [ D N T , b 4 ∣ Y N T ] = B ( ω N T 2 ( Y N T ) ) 2 E [ N 2 D N T , 1 4 ∣ Y N T ] = o P ( 1 ) .
For a fixed number of matchings, apply Theorem 1.1 of Chen and Fang [5 ] to a single permutation of C \mathcal{C} C . With s 2 = N ω N T 2 ( Y N T ) s^2=N\omega_{NT}^2(Y_{NT}) s 2 = N ω N T 2 ( Y N T ) , its conditional normal approximation satisfies
sup x ∣ Pr ( N D N T , 1 / s ≤ x ∣ Y N T ) − Φ ( x ) ∣ ≤ C 0 T s 3 ∑ i , t ≤ T ∣ C i t ∣ 3 = o P ( 1 ) . \sup_x \left|\Pr(ND_{NT,1}/s\le x\mid Y_{NT})-\Phi(x)\right|\le\frac{C_0}{Ts^3}\sum_{i,t\le T}|C_{it}|^3=o_P(1). x sup ∣ Pr ( N D N T , 1 / s ≤ x ∣ Y N T ) − Φ ( x ) ∣ ≤ T s 3 C 0 i , t ≤ T ∑ ∣ C i t ∣ 3 = o P ( 1 ) .
Here C 0 C_0 C 0 is a universal constant. The entries of the first N N N rows before padding have uniformly bounded third moments. Adding the zero rows leaves the sum of absolute cubes unchanged, and subtracting the column means increases it by at most a constant factor by Jensen’s inequality. Hence T − 1 ∑ ∣ C i t ∣ 3 = O P ( N ) T^{-1}\sum|C_{it}|^3=O_P(N) T − 1 ∑ ∣ C i t ∣ 3 = O P ( N ) , whereas s 3 / N 3 / 2 s^3/N^{3/2} s 3 / N 3/2 is bounded away from zero in probability. Conditional independence gives the Gaussian approximation for the average of any fixed number of draws.
Every integer sequence B B B has a further subsequence along which it is constant or diverges. On either subsequence, for every real u u u , the preceding limits give
E [ exp { i u ( θ ^ N T , B − Y ‾ N T ) / V N T , B } ∣ Y N T ] ⟶ P exp ( − ρ u 2 / 2 ) . E\left[\exp\left\{iu(\widehat{\theta}_{NT,B}-\overline{Y}_{NT})/\sqrt{V_{NT,B}}\right\}\mid Y_{NT}\right]\overset{\mathrm{P}}{\longrightarrow}\exp(-\rho u^2/2). E [ exp { i u ( θ N T , B − Y N T ) / V N T , B } ∣ Y N T ] ⟶ P exp ( − ρ u 2 /2 ) .
The conditional characteristic functions are bounded, so convergence also holds in L 1 L^1 L 1 . Multiply by exp ( i u Y ‾ N T / V N T , B ) \exp(iu\overline{Y}_{NT}/\sqrt{V_{NT,B}}) exp ( i u Y N T / V N T , B ) and take expectations. The limiting characteristic function is exp [ − ( ℓ + ρ ) u 2 / 2 ] = exp ( − u 2 / 2 ) \exp[-(\ell+\rho)u^2/2]=\exp(-u^2/2) exp [ − ( ℓ + ρ ) u 2 /2 ] = exp ( − u 2 /2 ) . Equation (16) and Slutsky’s theorem give (17) , since a continuous limiting distribution gives convergence in Kolmogorov distance. The argument applies to every sequence of laws and thus establishes uniformity.
Uniformity under a marginal fourth-moment bound. For P ∈ Q N T , 4 ( K ) P\in Q_{NT,4}(K) P ∈ Q N T , 4 ( K ) , Jensen’s inequality gives E a N T , 1 4 , E b N T , 1 4 ≤ K Ea_{NT,1}^4,Eb_{NT,1}^4\le K E a N T , 1 4 , E b N T , 1 4 ≤ K . In the variance calculation this gives absolute errors O P ( N − 1 / 2 ) O_P(N^{-1/2}) O P ( N − 1/2 ) for σ ^ N T , R 2 \widehat{\sigma}^2_{NT,R} σ N T , R 2 and σ ^ N T , C 2 \widehat{\sigma}^2_{NT,C} σ N T , C 2 , hence relative contributions O P ( B / N ) O_P(B/\sqrt{N}) O P ( B / N ) . For the projection sums the fourth-moment bound becomes
1 V N T , B 2 ( E a N T , 1 4 N 3 + E b N T , 1 4 T 3 ) ≤ 2 K N 3 V N T , B 2 ≤ 2 K B 2 N . \frac{1}{V_{NT,B}^2}\left(\frac{Ea_{NT,1}^4}{N^3}+\frac{Eb_{NT,1}^4}{T^3}\right)\le\frac{2K}{N^3V_{NT,B}^2}\le\frac{2KB^2}{N}. V N T , B 2 1 ( N 3 E a N T , 1 4 + T 3 E b N T , 1 4 ) ≤ N 3 V N T , B 2 2 K ≤ N 2 K B 2 .
Both bounds vanish when B = o ( N ) B=o(\sqrt{N}) B = o ( N ) . The conditional approximation requires only the marginal fourth moment and consistency of the variance components established above. This proves the claim under the marginal fourth-moment bound. □ \square □
Limiting distributions and a boundary example
Proof of Proposition 3. Center the observations and set σ N T 2 = 1 \sigma_{NT}^{2}=1 σ N T 2 = 1 by rescaling. Write D N T , b = M N T , b − Y ˉ N T D_{NT,b}=M_{NT,b}-\bar{Y}_{NT} D N T , b = M N T , b − Y ˉ N T and ω N T 2 = N Var ( D N T , 1 ∣ Y N T ) \omega_{NT}^{2}=N\operatorname{Var}(D_{NT,1}\mid Y_{NT}) ω N T 2 = N Var ( D N T , 1 ∣ Y N T ) as in (23) . Total variance and (9) give
E ω N T 2 = N ( V N T , 1 − V N T , F ) = N B V N T , B ( 1 − w N T , B ) , N B V N T , B ≥ 1. E\omega_{NT}^{2}=N(V_{NT,1}-V_{NT,F})=NBV_{NT,B}(1-w_{NT,B}), \qquad NBV_{NT,B}\geq1. E ω N T 2 = N ( V N T , 1 − V N T , F ) = N B V N T , B ( 1 − w N T , B ) , N B V N T , B ≥ 1.
The marginal fourth-moment bound implies σ ^ N T , C 2 − σ N T , C 2 = o P ( 1 ) \hat{\sigma}_{NT,C}^{2}-\sigma_{NT,C}^{2}=o_{P}(1) σ ^ N T , C 2 − σ N T , C 2 = o P ( 1 ) and σ ^ N T , I 2 − σ N T , I 2 = o P ( 1 ) \hat{\sigma}_{NT,I}^{2}-\sigma_{NT,I}^{2}=o_{P}(1) σ ^ N T , I 2 − σ N T , I 2 = o P ( 1 ) by the variance-component bounds in Appendix B. These bounds do not involve B B B . Consequently (23) gives ω N T 2 − E ω N T 2 = o P ( 1 ) \omega_{NT}^{2}-E\omega_{NT}^{2}=o_{P}(1) ω N T 2 − E ω N T 2 = o P ( 1 ) on rectangles as well.
If w < 1 w<1 w < 1 , then lim inf E ω N T 2 ≥ 1 − w > 0 \liminf E\omega_{NT}^{2}\geq1-w>0 lim inf E ω N T 2 ≥ 1 − w > 0 . Conditional on the array, the matching errors are i.i.d. For B → ∞ B\to\infty B → ∞ , the bounded conditional fourth moment E [ N 2 D N T , 1 4 ∣ Y N T ] = O P ( 1 ) E[N^{2}D_{NT,1}^{4}\mid Y_{NT}]=O_{P}(1) E [ N 2 D N T , 1 4 ∣ Y N T ] = O P ( 1 ) gives a conditional Lyapunov ratio of order O P ( 1 / B ) O_{P}(1/B) O P ( 1/ B ) . For fixed B B B , the single-permutation bound following (23) gives a conditional normal approximation to each draw, since ω N T 2 \omega_{NT}^{2} ω N T 2 is bounded away from zero in probability. This argument also covers bounded or oscillating B B B by taking subsequences. Since ω N T 2 / E ω N T 2 → 1 \omega_{NT}^{2}/E\omega_{NT}^{2}\to1 ω N T 2 / E ω N T 2 → 1 in probability, for every fixed real t t t we obtain
E [ exp { i t θ ^ N T , B − Y ˉ N T V N T , B } ∣ Y N T ] ⟶ exp { − t 2 ( 1 − w ) / 2 } in L 1 . E\left[\left.\exp\left\{it\frac{\hat{\theta}_{NT,B}-\bar{Y}_{NT}}{\sqrt{V_{NT,B}}}\right\}\right|Y_{NT}\right]\longrightarrow\exp\{-t^{2}(1-w)/2\}\quad\text{in }L^{1}. E [ exp { i t V N T , B θ ^ N T , B − Y ˉ N T } Y N T ] ⟶ exp { − t 2 ( 1 − w ) /2 } in L 1 .
Here convergence in probability implies L 1 L^{1} L 1 convergence because the characteristic functions are bounded. If w = 1 w=1 w = 1 , the same display holds because the matching error on this scale converges to zero in mean square. Multiplying the conditional characteristic function by exp { i t ( Y ˉ N T − θ N T ) / V N T , B } \exp\{it(\bar{Y}_{NT}-\theta_{NT})/\sqrt{V_{NT,B}}\} exp { i t ( Y ˉ N T − θ N T ) / V N T , B } and taking expectations gives
E e i t ( θ ^ N T , B − θ N T ) / V N T , B = e − t 2 ( 1 − w ) / 2 E e i t w N T , B ( Y ˉ N T − θ N T ) / V N T , F + o ( 1 ) . E e^{it(\widehat{\theta}_{NT,B}-\theta_{NT})/\sqrt{V_{NT,B}}}=e^{-t^{2}(1-w)/2}E e^{it\sqrt{w_{NT,B}}(\bar{Y}_{NT}-\theta_{NT})/\sqrt{V_{NT,F}}}+o(1). E e i t ( θ N T , B − θ N T ) / V N T , B = e − t 2 ( 1 − w ) /2 E e i t w N T , B ( Y ˉ N T − θ N T ) / V N T , F + o ( 1 ) .
For w > 0 w>0 w > 0 this gives the stated independent sum. For w = 0 w=0 w = 0 the full-array error on the total standard deviation scale vanishes in mean square, so no full-array limit is required. If the variance share itself does not converge, the same argument applies along each subsequence on which it does.
To establish studentization when w = 0 w=0 w = 0 , let V ^ N T , F = σ ^ N T , R 2 / N + σ ^ N T , C 2 / T + σ ^ N T , I 2 / ( N T ) \widehat{V}_{NT,F}=\widehat{\sigma}^{2}_{NT,R}/N+\widehat{\sigma}^{2}_{NT,C}/T+\widehat{\sigma}^{2}_{NT,I}/(NT) V N T , F = σ N T , R 2 / N + σ N T , C 2 / T + σ N T , I 2 / ( N T ) . Centering the row means gives E σ ^ N T , R 2 = ( 1 − 1 / N ) ( σ N T , R 2 + σ N T , I 2 / T ) E\widehat{\sigma}^{2}_{NT,R}=(1-1/N)(\sigma^{2}_{NT,R}+\sigma^{2}_{NT,I}/T) E σ N T , R 2 = ( 1 − 1/ N ) ( σ N T , R 2 + σ N T , I 2 / T ) ; centering the column means gives E σ ^ N T , C 2 = ( 1 − 1 / T ) ( σ N T , C 2 + σ N T , I 2 / N ) E\widehat{\sigma}^{2}_{NT,C}=(1-1/T)(\sigma^{2}_{NT,C}+\sigma^{2}_{NT,I}/N) E σ N T , C 2 = ( 1 − 1/ T ) ( σ N T , C 2 + σ N T , I 2 / N ) . Double demeaning the residuals gives E σ ^ N T , I 2 ≤ σ N T , I 2 E\widehat{\sigma}^{2}_{NT,I}\leq\sigma^{2}_{NT,I} E σ N T , I 2 ≤ σ N T , I 2 . Adding the three contributions gives E V ^ N T , F ≤ σ N T , R 2 / N + σ N T , C 2 / T + 3 σ N T , I 2 / ( N T ) ≤ 3 V N T , F E\widehat{V}_{NT,F}\leq\sigma^{2}_{NT,R}/N+\sigma^{2}_{NT,C}/T+3\sigma^{2}_{NT,I}/(NT)\leq3V_{NT,F} E V N T , F ≤ σ N T , R 2 / N + σ N T , C 2 / T + 3 σ N T , I 2 / ( N T ) ≤ 3 V N T , F . Consequently, for every ϵ > 0 \epsilon>0 ϵ > 0 ,
Pr ( V ^ N T , F / V N T , B > ϵ ) ≤ 3 V N T , F ϵ V N T , B = 3 w N T , B ϵ ⟶ 0. \Pr\left(\widehat{V}_{NT,F}/V_{NT,B}>\epsilon\right)\leq\frac{3V_{NT,F}}{\epsilon V_{NT,B}}=\frac{3w_{NT,B}}{\epsilon}\longrightarrow0. Pr ( V N T , F / V N T , B > ϵ ) ≤ ϵ V N T , B 3 V N T , F = ϵ 3 w N T , B ⟶ 0.
By Appendix § 8 , σ ^ N T 2 / σ N T 2 → 1 \widehat{\sigma}^{2}_{NT}/\sigma^{2}_{NT}\to1 σ N T 2 / σ N T 2 → 1 in probability. Since
σ N T 2 N B V N T , B = 1 − ( 1 − 1 / B ) w N T , B → 1 , \frac{\sigma^{2}_{NT}}{NBV_{NT,B}}=1-(1-1/B)w_{NT,B}\to1, N B V N T , B σ N T 2 = 1 − ( 1 − 1/ B ) w N T , B → 1 ,
(14) gives V ^ N T , B / V N T , B → 1 \widehat{V}_{NT,B}/V_{NT,B}\to1 V N T , B / V N T , B → 1 in probability. □
Scales in Table 1 . If ν N T → c < 1 \nu_{NT}\to c<1 ν N T → c < 1 , § 3.3 gives V N T , 1 / ( T V N T , F ) → 1 − c V_{NT,1}/(TV_{NT,F})\to1-c V N T , 1 / ( T V N T , F ) → 1 − c . Using V N T , 1 = σ N T 2 / N V_{NT,1}=\sigma^{2}_{NT}/N V N T , 1 = σ N T 2 / N in (9) , we obtain
N T V N T , F σ N T 2 ⟶ 1 1 − c , N T ( V N T , B − V N T , F ) σ N T 2 = T B ( 1 − N V N T , F σ N T 2 ) . \frac{NTV_{NT,F}}{\sigma^{2}_{NT}}\longrightarrow\frac{1}{1-c},\qquad\frac{NT(V_{NT,B}-V_{NT,F})}{\sigma^{2}_{NT}}=\frac{T}{B}\left(1-\frac{NV_{NT,F}}{\sigma^{2}_{NT}}\right). σ N T 2 N T V N T , F ⟶ 1 − c 1 , σ N T 2 N T ( V N T , B − V N T , F ) = B T ( 1 − σ N T 2 N V N T , F ) .
Thus for B / T → λ ∈ ( 0 , ∞ ) B/T\to\lambda\in(0,\infty) B / T → λ ∈ ( 0 , ∞ ) the two standard deviations on the N T / σ N T \sqrt{NT}/\sigma_{NT} N T / σ N T scale tend to 1 / 1 − c 1/\sqrt{1-c} 1/ 1 − c and 1 / λ 1/\sqrt{\lambda} 1/ λ , respectively. For B / T → ∞ B/T\to\infty B / T → ∞ the second tends to zero. For B = o ( T ) B=o(T) B = o ( T ) use the N B / σ N T \sqrt{NB}/\sigma_{NT} N B / σ N T scale instead: N B V N T , F / σ N T 2 → 0 NBV_{NT,F}/\sigma_{NT}^{2} \to0 N B V N T , F / σ N T 2 → 0 and N B V N T , B / σ N T 2 = 1 + ( B − 1 ) N V N T , F / σ N T 2 → 1 NBV_{NT,B}/\sigma_{NT}^{2} = 1 + (B - 1)NV_{NT,F}/\sigma_{NT}^{2} \to1 N B V N T , B / σ N T 2 = 1 + ( B − 1 ) N V N T , F / σ N T 2 → 1 .
When the projections dominate. Let V N T , P = σ N T , R 2 / N + σ N T , C 2 / T = ν N T V N T , F V_{NT,P} = \sigma_{NT,R}^{2}/N + \sigma_{NT,C}^{2}/T = \nu_{NT}V_{NT,F} V N T , P = σ N T , R 2 / N + σ N T , C 2 / T = ν N T V N T , F and suppose ν N T → 1 \nu_{NT} \to1 ν N T → 1 . The intersection average divided by V N T , F \sqrt{V_{NT,F}} V N T , F has variance 1 − ν N T → 0 1-\nu_{NT} \to0 1 − ν N T → 0 . The independent projection summands enter the mean with weights 1 / N 1/N 1/ N and 1 / T 1/T 1/ T . Their Lindeberg condition requires, for every ϵ > 0 \epsilon> 0 ϵ > 0 ,
E [ a N T , 1 2 1 { ∣ a N T , 1 ∣ > ϵ N V N T , P } ] N V N T , P + E [ b N T , 1 2 1 { ∣ b N T , 1 ∣ > ϵ T V N T , P } ] T V N T , P → 0. \frac{E[a_{NT,1}^{2}1\{|a_{NT,1}| > \epsilon N\sqrt{V_{NT,P}}\}]}{NV_{NT,P}} + \frac{E[b_{NT,1}^{2}1\{|b_{NT,1}| > \epsilon T\sqrt{V_{NT,P}}\}]}{TV_{NT,P}} \to0. N V N T , P E [ a N T , 1 2 1 { ∣ a N T , 1 ∣ > ϵ N V N T , P }] + T V N T , P E [ b N T , 1 2 1 { ∣ b N T , 1 ∣ > ϵ T V N T , P }] → 0.
Under this condition the normalized projection sum is asymptotically standard normal. Since V N T , P / V N T , F → 1 V_{NT,P}/V_{NT,F} \to1 V N T , P / V N T , F → 1 and the residual average is negligible on this scale, the same limit holds for the standardized full-array mean. In particular, for laws in P N T ( K ) \mathcal{P}_{NT}(K) P N T ( K ) , the Lyapunov ratio satisfies
E [ a N T , 1 4 ] / N 3 + E [ b N T , 1 4 ] / T 3 V N T , P 2 ≤ K σ N T 2 N 2 V N T , P ≤ K T N ( 1 N + 1 − ν N T ν N T ) . \frac{E[a_{NT,1}^{4}]/N^{3} + E[b_{NT,1}^{4}]/T^{3}}{V_{NT,P}^{2}} \le\frac{K\sigma_{NT}^{2}}{N^{2}V_{NT,P}} \le K\frac{T}{N}\left(\frac{1}{N} + \frac{1-\nu_{NT}}{\nu_{NT}}\right). V N T , P 2 E [ a N T , 1 4 ] / N 3 + E [ b N T , 1 4 ] / T 3 ≤ N 2 V N T , P K σ N T 2 ≤ K N T ( N 1 + ν N T 1 − ν N T ) .
This tends to zero if T / N T/N T / N is bounded. Without such a restriction the conclusion need not hold. For example, take T = N 2 T=N^{2} T = N 2 and Y N T , i t = I N T , i − 1 / N + ζ i t Y_{NT,it}=I_{NT,i}-1/N+\zeta_{it} Y N T , i t = I N T , i − 1/ N + ζ i t , with independent Bernoulli( 1 / N ) (1/N) ( 1/ N ) row shocks and independent standard normal intersection shocks, independent of the row shocks. Here σ N T , R 2 = ( 1 − 1 / N ) / N \sigma_{NT,R}^{2}=(1-1/N)/N σ N T , R 2 = ( 1 − 1/ N ) / N , σ N T , C 2 = 0 \sigma_{NT,C}^{2}=0 σ N T , C 2 = 0 , and σ N T , I 2 = 1 \sigma_{NT,I}^{2}=1 σ N T , I 2 = 1 . The moment bounds defining P N T ( K ) \mathcal{P}_{NT}(K) P N T ( K ) hold for a common K K K , since the row projection is bounded by one and σ N T 2 ≥ 1 \sigma_{NT}^{2}\ge1 σ N T 2 ≥ 1 . Yet N 2 V N T , F = 1 N^{2}V_{NT,F}=1 N 2 V N T , F = 1 and ν N T = 1 − 1 / N → 1 \nu_{NT}=1-1/N\to1 ν N T = 1 − 1/ N → 1 , while N Y ‾ N T ⇒ Poisson ( 1 ) − 1 N\overline{Y}_{NT}\Rightarrow\operatorname{Poisson}(1)-1 N Y N T ⇒ Poisson ( 1 ) − 1 : the row count converges to a Poisson law and the intersection average has variance 1 / N 1/N 1/ N on this scale.
Sharpness under a marginal fourth-moment bound. Suppose B = o ( N ) B=o(N) B = o ( N ) and B / N B/\sqrt{N} B / N does not tend to zero. Work along a subsequence on which B ≥ c N B\ge c\sqrt{N} B ≥ c N for some c > 0 c>0 c > 0 . For these array dimensions N ≤ T N\le T N ≤ T , let
Y N T , i t = a N T , i + e i t , a N T , i = N / B ( I N T , i − 1 / N ) , I N T , i ∼ i . i . d . Bernoulli ( 1 / N ) , Y_{NT,it}=a_{NT,i}+e_{it}, \qquad a_{NT,i}=\sqrt{N/B}\left(I_{NT,i}-1/N\right), \qquad I_{NT,i}\mathrel{\overset{\mathrm{i.i.d.}}{\sim}}\operatorname{Bernoulli}(1/N), Y N T , i t = a N T , i + e i t , a N T , i = N / B ( I N T , i − 1/ N ) , I N T , i ∼ i.i.d. Bernoulli ( 1/ N ) ,
where the e i t e_{it} e i t are independent standard normals independent of the I N T , i I_{NT,i} I N T , i . Here θ N T = 0 \theta_{NT}=0 θ N T = 0 , σ N T , R 2 = ( 1 − 1 / N ) / B \sigma^2_{NT,R}=(1-1/N)/B σ N T , R 2 = ( 1 − 1/ N ) / B , σ N T , C 2 = 0 \sigma^2_{NT,C}=0 σ N T , C 2 = 0 , σ N T , I 2 = 1 \sigma^2_{NT,I}=1 σ N T , I 2 = 1 , and E [ a N T , i 4 ] ∼ N / B 2 = O ( 1 ) E[a^4_{NT,i}]\sim N/B^2=O(1) E [ a N T , i 4 ] ∼ N / B 2 = O ( 1 ) . Hence the standardized marginal fourth moment is bounded, while E [ a N T , i 4 ] / ( σ N T 2 σ N T , R 2 ) ∼ N / B → ∞ E[a^4_{NT,i}]/(\sigma^2_{NT}\sigma^2_{NT,R})\sim N/B\to\infty E [ a N T , i 4 ] / ( σ N T 2 σ N T , R 2 ) ∼ N / B → ∞ .
Each matching uses every row once, so the row contribution is
N B N ∑ i a N T , i = ∑ i = 1 N I N T , i − 1 ⇒ Q − 1 , Q ∼ Poisson ( 1 ) . \frac{\sqrt{NB}}{N}\sum_i a_{NT,i}=\sum_{i=1}^{N}I_{NT,i}-1\Rightarrow Q-1,\qquad Q\sim\operatorname{Poisson}(1). N N B i ∑ a N T , i = i = 1 ∑ N I N T , i − 1 ⇒ Q − 1 , Q ∼ Poisson ( 1 ) .
Let K N T , i t K_{NT,it} K N T , i t count the number of times intersection ( i , t ) (i,t) ( i , t ) is selected. For any row i i i and distinct matchings b , b ′ b,b' b , b ′ , the assigned columns coincide with probability 1 / T 1/T 1/ T . Summing over N N N rows and B ( B − 1 ) B(B-1) B ( B − 1 ) ordered pairs of matchings gives
E ∑ i , t K N T , i t ( K N T , i t − 1 ) = N B ( B − 1 ) T . E\sum_{i,t}K_{NT,it}(K_{NT,it}-1)=\frac{NB(B-1)}{T}. E i , t ∑ K N T , i t ( K N T , i t − 1 ) = T N B ( B − 1 ) .
Since ∑ i , t K N T , i t = N B \sum_{i,t}K_{NT,it}=NB ∑ i , t K N T , i t = N B and B = o ( N ) B=o(N) B = o ( N ) with N ≤ T N\le T N ≤ T ,
E [ ∑ i , t K N T , i t 2 N B − 1 ] = B − 1 T → 0. E\left[\frac{\sum_{i,t}K^2_{NT,it}}{NB}-1\right]=\frac{B-1}{T}\to0. E [ N B ∑ i , t K N T , i t 2 − 1 ] = T B − 1 → 0.
The quantity in brackets is nonnegative, so ∑ i , t K N T , i t 2 / ( N B ) → 1 \sum_{i,t}K^2_{NT,it}/(NB)\to1 ∑ i , t K N T , i t 2 / ( N B ) → 1 in probability by Markov’s inequality. Conditional on the matchings, the intersection error in N B θ ^ N T , B \sqrt{NB}\hat{\theta}_{NT,B} N B θ ^ N T , B is Gaussian with this variance and is independent of the row projections. Thus N B θ ^ N T , B ⇒ Q − 1 + Z \sqrt{NB}\hat{\theta}_{NT,B}\Rightarrow Q-1+Z N B θ ^ N T , B ⇒ Q − 1 + Z , where Z Z Z is standard normal and independent of Q Q Q .
The variance estimate V ^ N T , B \hat{V}_{NT,B} V ^ N T , B is positive almost surely, so the studentized statistic has the same sign as θ ^ N T , B \hat{\theta}_{NT,B} θ ^ N T , B . Its distribution function at zero converges to
E Φ ( 1 − Q ) ≥ e − 1 { Φ ( 1 ) + 1 2 + 1 2 Φ ( − 1 ) } > 0.52. E\Phi(1-Q)\ge e^{-1}\left\{\Phi(1)+\frac{1}{2}+\frac{1}{2}\Phi(-1)\right\}>0.52. E Φ ( 1 − Q ) ≥ e − 1 { Φ ( 1 ) + 2 1 + 2 1 Φ ( − 1 ) } > 0.52.
Since Φ ( 0 ) = 1 / 2 \Phi(0)=1/2 Φ ( 0 ) = 1/2 , a marginal fourth-moment bound alone cannot guarantee uniform Gaussian approximation for these choices of B B B . This includes the boundary choice B = ⌊ N ⌋ B=\lfloor\sqrt{N}\rfloor B = ⌊ N ⌋ .
Proof of Theorem 2
Work along any admissible sequence of laws. Applying (9) to a scalar function x x x of an intersection observation bounds the variance of its matching average by Var ( x N T , 11 ) / N \operatorname{Var}(x_{NT,11})/N Var ( x N T , 11 ) / N . At β 0 , N T \beta_{0,NT} β 0 , N T , this bound applied componentwise gives
m ^ N T , B ( β 0 , N T ) = O P ( N − 1 / 2 ) , ∂ m ^ N T , B ( β 0 , N T ) ∂ β ⊤ − D N T = O P ( N − 1 / 2 ) . \widehat{m}_{NT,B}(\beta_{0,NT}) = O_{\mathrm{P}}(N^{-1/2}), \qquad\frac{\partial\widehat{m}_{NT,B}(\beta_{0,NT})}{\partial\beta^{\top}} - D_{NT} = O_{\mathrm{P}}(N^{-1/2}). m N T , B ( β 0 , N T ) = O P ( N − 1/2 ) , ∂ β ⊤ ∂ m N T , B ( β 0 , N T ) − D N T = O P ( N − 1/2 ) .
Applied on a finite net of Θ \Theta Θ , the same bound gives convergence to the population moments on that net. The sampled and population averages of the derivative envelope are O P ( 1 ) O_{\mathrm{P}}(1) O P ( 1 ) and O ( 1 ) O(1) O ( 1 ) , respectively. Taking finer nets therefore gives uniform convergence over Θ \Theta Θ . The weight bounds and uniform identification imply β ^ N T , B − β 0 , N T = o P ( 1 ) \widehat{\beta}_{NT,B} - \beta_{0,NT} = o_{\mathrm{P}}(1) β N T , B − β 0 , N T = o P ( 1 ) , so the minimizer is interior with probability tending to one.
Write Δ N T , B = β ^ N T , B − β 0 , N T \Delta_{NT,B} = \widehat{\beta}_{NT,B} - \beta_{0,NT} Δ N T , B = β N T , B − β 0 , N T , and let D ‾ N T , B \overline{D}_{NT,B} D N T , B be the average sampled derivative along the segment joining these parameter values. The second-derivative bound and consistency give D ^ N T , B − D N T = o P ( 1 ) \widehat{D}_{NT,B} - D_{NT} = o_{\mathrm{P}}(1) D N T , B − D N T = o P ( 1 ) and D ‾ N T , B − D N T = o P ( 1 ) \overline{D}_{NT,B} - D_{NT} = o_{\mathrm{P}}(1) D N T , B − D N T = o P ( 1 ) . With probability tending to one, the solution is interior and the derivative product below is invertible. Substitution in the first-order condition gives
0 = D ^ N T , B ⊤ W ^ N T { m ^ N T , B ( β 0 , N T ) + D ‾ N T , B Δ N T , B } , 0 = \widehat{D}_{NT,B}^{\top}\widehat{W}_{NT}\left\{\widehat{m}_{NT,B}(\beta_{0,NT}) + \overline{D}_{NT,B}\Delta_{NT,B}\right\}, 0 = D N T , B ⊤ W N T { m N T , B ( β 0 , N T ) + D N T , B Δ N T , B } ,
Δ N T , B = − A ~ N T , B m ^ N T , B ( β 0 , N T ) = O P ( N − 1 / 2 ) . \Delta_{NT,B} = -\widetilde{A}_{NT,B}\widehat{m}_{NT,B}(\beta_{0,NT}) = O_{\mathrm{P}}(N^{-1/2}). Δ N T , B = − A N T , B m N T , B ( β 0 , N T ) = O P ( N − 1/2 ) .
where A ~ N T , B = ( D ^ N T , B ⊤ W ^ N T D ‾ N T , B ) − 1 D ^ N T , B ⊤ W ^ N T \widetilde{A}_{NT,B} = (\widehat{D}_{NT,B}^{\top}\widehat{W}_{NT}\overline{D}_{NT,B})^{-1}\widehat{D}_{NT,B}^{\top}\widehat{W}_{NT} A N T , B = ( D N T , B ⊤ W N T D N T , B ) − 1 D N T , B ⊤ W N T . The derivative errors are now O P ( N − 1 / 2 ) O_{\mathrm{P}}(N^{-1/2}) O P ( N − 1/2 ) , so A ~ N T , B − A N T = O P ( N − 1 / 2 ) \widetilde{A}_{NT,B} - A_{NT} = O_{\mathrm{P}}(N^{-1/2}) A N T , B − A N T = O P ( N − 1/2 ) . A Taylor expansion of ℓ \ell ℓ and the definition of u N T , i u_{NT,i} u N T , i in (19) give
ℓ ( β ^ N T , B ) − ℓ ( β 0 , N T ) = 1 N B ∑ b , i u N T , i π b ( i ) + R N T , B , \ell(\widehat{\beta}_{NT,B}) - \ell(\beta_{0,NT}) = \frac{1}{NB}\sum_{b,i}u_{NT,i}\pi_{b(i)} + R_{NT,B}, ℓ ( β N T , B ) − ℓ ( β 0 , N T ) = N B 1 b , i ∑ u N T , i π b ( i ) + R N T , B ,
∣ R N T , B ∣ ≤ C ( ∥ Δ N T , B ∥ 2 + ∥ A ~ N T , B − A N T ∥ ∥ m ^ N T , B ( β 0 , N T ) ∥ ) = O P ( N − 1 ) . \lvert R_{NT,B}\rvert\leq C\left(\lVert\Delta_{NT,B}\rVert^{2} + \lVert\widetilde{A}_{NT,B} - A_{NT}\rVert\lVert\widehat{m}_{NT,B}(\beta_{0,NT})\rVert\right) = O_{\mathrm{P}}(N^{-1}). ∣ R N T , B ∣ ≤ C ( ∥ Δ N T , B ∥ 2 + ∥ A N T , B − A N T ∥ ∥ m N T , B ( β 0 , N T )∥ ) = O P ( N − 1 ) .
By (9) , V N T , B ( u N T ) ≥ Var ( u N T , 11 ) / ( N B ) ≥ v / ( N B ) V_{NT,B}(u_{NT}) \ge\operatorname{Var}(u_{NT,11})/(NB) \ge v/(NB) V N T , B ( u N T ) ≥ Var ( u N T , 11 ) / ( N B ) ≥ v / ( N B ) . Thus R N T , B / V N T , B ( u N T ) = O P ( B / N ) = o P ( 1 ) R_{NT,B}/\sqrt{V_{NT,B}(u_{NT})} = O_P(\sqrt{B/N}) = o_P(1) R N T , B / V N T , B ( u N T ) = O P ( B / N ) = o P ( 1 ) .
For feasible studentization, the full-array average of ( u ^ N T , i t − u N T , i t ) 2 (\widehat{u}_{NT,it}-u_{NT,it})^2 ( u N T , i t − u N T , i t ) 2 is O P ( N − 1 ) O_P(N^{-1}) O P ( N − 1 ) . Indeed, the coefficient vector in (19) has estimation error O P ( N − 1 / 2 ) O_P(N^{-1/2}) O P ( N − 1/2 ) , and the full-array average of the squared envelope is O P ( 1 ) O_P(1) O P ( 1 ) . Formula (14) is a sum of squared linear functions with nonnegative coefficients, so its square root is a seminorm. The row, column, and residual sums of squares add to the centered marginal sum of squares. With N ≤ T N \le T N ≤ T , this gives V ^ N T , B ( x ) ≤ N − 1 ( N T ) − 1 ∑ i , t x N T , i t 2 \widehat{V}_{NT,B}(x) \le N^{-1}(NT)^{-1}\sum_{i,t}x_{NT,it}^2 V N T , B ( x ) ≤ N − 1 ( N T ) − 1 ∑ i , t x N T , i t 2 . Hence
∣ V ^ N T , B ( u ^ N T ) − V ^ N T , B ( u N T ) ∣ ≤ V ^ N T , B ( u N T − u N T ) ≤ ( 1 N 2 T ∑ i , t ( u ^ N T , i t − u N T , i t ) 2 ) 1 / 2 = O P ( N − 1 ) . \begin{aligned}
\left\lvert\sqrt{\widehat{V}_{NT,B}(\widehat{u}_{NT})}-\sqrt{\widehat{V}_{NT,B}(u_{NT})}\right\rvert\le\sqrt{\widehat{V}_{NT,B}(u_{NT}-u_{NT})} \\
&\le\left(\frac{1}{N^2T}\sum_{i,t}(\widehat{u}_{NT,it}-u_{NT,it})^2\right)^{1/2}=O_P(N^{-1}).
\end{aligned} V N T , B ( u N T ) − V N T , B ( u N T ) ≤ V N T , B ( u N T − u N T ) ≤ ( N 2 T 1 i , t ∑ ( u N T , i t − u N T , i t ) 2 ) 1/2 = O P ( N − 1 ) .
Applying Theorem 1 to u N T u_{NT} u N T proves (20) , including the assertion under B = o ( N ) B=o(\sqrt{N}) B = o ( N ) .
Joint inference. For a unit vector c c c , take ℓ ( β ) = c ⊤ β \ell(\beta)=c^\top\beta ℓ ( β ) = c ⊤ β and denote its influence contribution by u N T ( c ) u_{NT}(c) u N T ( c ) . Suppose the variance and moment bounds of Theorem 2 hold with common constants for all such c c c . Define covariance matrices S N T , B S_{NT,B} S N T , B and S ^ N T , B \widehat{S}_{NT,B} S N T , B by the quadratic forms c ⊤ S N T , B c = V N T , B ( u N T ( c ) ) c^\top S_{NT,B}c=V_{NT,B}(u_{NT}(c)) c ⊤ S N T , B c = V N T , B ( u N T ( c )) and c ⊤ S ^ N T , B c = V ^ N T , B ( u ^ N T ( c ) ) c^\top\widehat{S}_{NT,B}c=\widehat{V}_{NT,B}(\widehat{u}_{NT}(c)) c ⊤ S N T , B c = V N T , B ( u N T ( c )) . Since the preceding argument applies to every deterministic sequence of unit directions, Cramér–Wold applied in directions proportional to S N T , B − 1 / 2 c S_{NT,B}^{-1/2}c S N T , B − 1/2 c gives S N T , B − 1 / 2 ( β ^ N T , B − β 0 , N T ) ⇒ N ( 0 , I d ) S_{NT,B}^{-1/2}(\widehat{\beta}_{NT,B}-\beta_{0,NT})\Rightarrow\mathcal{N}(0,I_d) S N T , B − 1/2 ( β N T , B − β 0 , N T ) ⇒ N ( 0 , I d ) . Directional variance consistency and polarization on the same scale give S N T , B − 1 / 2 S ^ N T , B S N T , B − 1 / 2 → I d S_{NT,B}^{-1/2}\widehat{S}_{NT,B}S_{NT,B}^{-1/2}\to I_d S N T , B − 1/2 S N T , B S N T , B − 1/2 → I d in probability. Thus the inverse exists with probability tending to one and ( β ^ N T , B − β 0 , N T ) ⊤ S ^ N T , B − 1 ( β ^ N T , B − β 0 , N T ) ⇒ χ d 2 (\widehat{\beta}_{NT,B}-\beta_{0,NT})^\top\widehat{S}_{NT,B}^{-1}(\widehat{\beta}_{NT,B}-\beta_{0,NT})\Rightarrow\chi_d^2 ( β N T , B − β 0 , N T ) ⊤ S N T , B − 1 ( β N T , B − β 0 , N T ) ⇒ χ d 2 . The arbitrary-sequence argument establishes uniformity. □