Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
78,327 characters · 13 sections · 35 citation commands
0.5mm Model Checks in a Kernel Ridge Regression Framework 0.25mm
Keywords: Conditional moment restriction, Kernel ridge regression, Model checking, Reproducing kernel Hilbert space
\pagenumbering{arabic} \setcounter{page}{1} \setcounter{footnote}{0}
A wide range of models in the social sciences can be described by conditional moment restrictions of the form
where the random vector $S = (Y, X)$ is assumed to be strictly stationary with probability measure $\mathbb{P}_S$, $Y \in \mathcal{Y} \subset \mathbb{R}^q$, and $X \in \mathcal{X} \subset \mathbb{R}^d$. The residual function $\varepsilon: \mathcal{S} \times \Theta \to \mathbb{R}^q$ is known up to the finite-dimensional parameter $\theta_0 \in \Theta$, with $\Theta \subset \mathbb{R}^d$. For clarity, we focus on the univariate case $q=1$; the extension to $q > 1$ is straightforward, as the relevant test statistics are simply summed over the components of $\varepsilon_0$ (see Remark (ref) in Section (ref)). Throughout, we denote the residual for observation $i$ as $\varepsilon(s_i; \theta)$, or simply $\varepsilon_{\theta, i}$ when the context is clear.
This framework encompasses many important models, including linear and nonlinear conditional mean regression, quantile regression, treatment effect models, and instrumental variables regression, among others.
There has been longstanding interest in developing goodness-of-fit tests and model checks for such models in statistics and econometrics. The objective is to test \[ H_0: \mathbb{E}(\varepsilon_0 \mid X) = 0 \quad \text{vs.} \quad H_1: \mathbb{P}\left(\mathbb{E}(\varepsilon_0 \mid X) = 0\right) < 1, \quad \forall \theta \in \Theta. \] A straightforward approach to testing this hypothesis is to construct a chi-square statistic based on a finite set of moment conditions. For instance, one could regress the estimated residuals on the covariate $X$: \[ \varepsilon_{\hat \theta} = X^\top \omega + e \] and test the joint null hypothesis $\omega = \boldsymbol{0}$ versus the alternative $\omega \neq \boldsymbol{0}$. Here, $\hat \theta$ is a consistent estimator.
However, this approach has two major drawbacks: (i) the test is inconsistent, as it can only detect linear deviations represented by a finite set of moments; and (ii) the use of the estimated parameter $\hat \theta$ complicates the null distribution, making the resulting test statistic non-pivotal. As a result, modern goodness-of-fit testing techniques,which will be reviewed later, have moved beyond this simple idea, often introducing additional complexity that can make practical implementation and interpretation more challenging for practitioners.
In this paper, we build on the simple regression idea, but instead of relying on a finite set of moment conditions, we map the covariates into an infinite-dimensional reproducing kernel Hilbert space (RKHS) $\mathcal{H}_k$ using a reproducing kernel $k(\cdot,\cdot)$, i.e., $k: \mathcal{X} \to \mathcal{H}_k$. The corresponding regression model becomes: \[ \varepsilon_{\hat \theta} = \langle k(X,\cdot), w^*(\hat \theta) \rangle_{\mathcal{H}_k} + e, \] where $\mathbb{E}(e|X)=0$ by construction, $k(X,\cdot)$ is the kernel function evaluated at the observed data point, $w^*(\hat \theta) \in \mathcal{H}_k$ is the population coefficient function indexed by $\hat \theta$, and $\langle \cdot, \cdot \rangle_{\mathcal{H}_k}$ denotes the RKHS inner product. As shown in the proof of Lemma 1 (see Appendix (ref)), \[ w^*(\hat \theta) = C^{-1} \mathbb{E}\left[\varepsilon_{\hat \theta} k(X,\cdot)\right], \] where $C = \mathbb{E}[k(X,\cdot) \otimes k(X,\cdot)]$ is the second moment operator, with $\otimes$ being the tensor product. However, due to the infinite-dimensionality, $C$ is generally not invertible, so we regularize it in the sense of Tikhonov regularization: \[ C_\lambda = C + \lambda \mathcal{I}, \] where $\mathcal{I}$ is the identity operator on $\mathcal{H}_k$ and $\lambda>0$ is fixed.
The regularized population coefficient function $w^*_\lambda(\hat \theta)$ is defined as \[ w^*_\lambda(\hat \theta) = C_\lambda^{-1} \mathbb{E}\left[\varepsilon_{\hat \theta} k(X,\cdot)\right], \] This function coincides with the coefficient in kernel ridge regression (KRR) with fixed $\lambda$, a widely used method in machine learning for modeling complex relationships. Given data $\{(x_i, \varepsilon_{\hat \theta,i})\}_{i=1}^n$, the KRR estimator is \[ \hat w = \boldsymbol{\varepsilon}_{\hat \theta}^\top \left(\boldsymbol{K} + n \lambda \boldsymbol{I}\right)^{-1} \boldsymbol{\Phi}, \] where $\boldsymbol{I}$ is the $n \times n$ identity matrix, $\boldsymbol{\varepsilon}_{\hat \theta} = (\varepsilon_{\hat \theta,1}, \ldots, \varepsilon_{\hat \theta,n})^\top$ is the vector of estimated residuals, $\boldsymbol{\Phi} = (k(x_1, \cdot), \ldots, k(x_n, \cdot))^\top$ is the vector of kernel functions evaluated at the observed data points, and $\boldsymbol{K} = (k(x_i, x_j))_{i,j=1}^n$ is the kernel matrix.
We will show in Section (ref) that the null hypothesis holds if and only if $w^*_\lambda(\theta_0) = w^*_\lambda = \boldsymbol{0}$ in $\mathcal{H}_k$. The intuition is as follows: suppose the true relationship between the residual and the covariates is \[ \varepsilon_{0} = f^*(X) + e^*, \] where $f^* \in C(\mathcal{X})$ (the space of continuous functions on $\mathcal{X}$) and $e^*$ is a random error term. If the kernel is universal, the associated RKHS $\mathcal{H}_k$ is dense in $C(\mathcal{X})$ micchelli2006universal, so $w^*_\lambda$ can approximate $f^*$ arbitrarily well. Under the null, $f^* = \boldsymbol{0}$, which implies $w^*_\lambda = \boldsymbol{0}$ in the RKHS. Conversely, if $w^*_\lambda \neq \boldsymbol{0}$, this signals model misspecification and a nonzero expected residual.
Focusing on $w^*_\lambda$ offers two main advantages: (i) as a natural infinite-dimensional analogue of the regression coefficient, it provides a powerful means to capture and test for deviations from the model; and (ii) it yields an interpretable and visualizable summary of model inadequacy, as the structure of $w^*_\lambda$ in the RKHS reflects the nature of the misspecification.
To illustrate the first point, we use Lemma 1 in Section (ref), which shows that $w^*_\lambda$ can be represented as \[ w^*_\lambda = \sum_{i=1}^\infty \frac{\mu_i}{\mu_i + \lambda} \mathbb{E}(\varepsilon_{0} \phi_i(X)) \phi_i, \] where $\{\phi_i\}_{i=1}^\infty$ and $\{\mu_i\}_{i=1}^\infty$ are the eigenfunctions and eigenvalues of the integral operator associated with the kernel. The eigenvalues are positive and decay to zero at a certain rate. The term $\mathbb{E}(\varepsilon_0 \phi_i(X))$ quantifies the deviation in the direction of the $i$-th eigenfunction. When deviations occur in high-frequency directions (i.e., those with small eigenvalues), the corresponding weights $\mu_i / (\mu_i + \lambda)$ can still be substantial if the regularization parameter $\lambda$ is chosen appropriately. This contrasts with the test statistics proposed by muandet2020kernel,escanciano2024gaussian, where the weights for each directional deviation are simply the eigenvalues themselves (see Equation (18) in escanciano2024gaussian). However, this does not imply that tests based on $w^*_\lambda$ are always superior. For example, if deviations occur in low-frequency directions (with large eigenvalues), the weights $\mu_i / (\mu_i + \lambda)$ will be less than one, whereas the tests of muandet2020kernel,escanciano2024gaussian could assign larger weights to these directions. Thus, the relative performance depends on the nature of the deviation.
The second point, that $w^*_\lambda$ is interpretable and can be visualized, is particularly valuable for practitioners. While most existing goodness-of-fit tests yield only a binary decision (reject or not reject the null hypothesis), researchers are often interested in diagnosing the nature of model misspecification: specifically, identifying regions in the covariate space where the model deviates most or least from the data. Lemmas 2 and 3 in Section (ref) show that if the kernel is analytic, one can construct test statistics by evaluating $w^*_\lambda$ at randomly chosen locations drawn from any distribution with a Lebesgue density. By sampling a sufficient number of such locations, one can generate visualizations that reveal the “bumps” or patterns of deviation in the covariate space, as illustrated in Figure (ref). This approach is analogous to using residual plots to assess model fit, providing an intuitive and informative diagnostic tool.
Building on $w^*_\lambda$, we introduce two classes of model checking procedures: (i) projection-based tests, which assess the projection of $w^*_\lambda$ onto functions in $\mathcal{H}_k$ to capture global deviations from the null; and (ii) random location tests, which exploit the property that if $w^*_\lambda$ is analytic and equals zero in $\mathcal{H}_k$, then $w^*_\lambda(v) = 0$ almost surely for any $v \in \mathcal{X}$ drawn from a distribution with a Lebesgue density. Here, $w^*_\lambda(v) $ serves as a metric that captures local deviations from the null.
For the projection-based tests, we consider two specific projections: (i) onto $w^*_\lambda$ itself: $\langle w^*_\lambda, w^*_\lambda\rangle_{\mathcal{H}_k} $ and (ii) onto the mean embedding of the residual $m^* = \mathbb{E}(\varepsilon_0 k(X,\cdot))$: $\langle w^*_\lambda, m^* \rangle_{\mathcal{H}_k}$. The resulting two test statistics share identical structure expect for the weights applied to the directional deviations. For the random location tests with $J$ location points, we propose two test statistics that are based on $\sum_{j=1}^{J}(w^*_\lambda(v_j))^2$ and $(\sum_{j=1}^{J}w^*_\lambda(v_j))^2$.
When the residuals are computed using the estimated parameter $\hat \theta$, the resulting test statistics would have additional complexity in their null distributions due to the estimation effect. To address this, we follow the Neyman orthogonality approach, as advocated by escanciano2014specification,escanciano2024gaussian,sant2019specification. Specifically, we introduce a projection operator $\boldsymbol{\Pi}$ that acts on the residual function $\varepsilon_{\theta}$, and redefine the key functions $w^*_\lambda$ and $m^*$ as $w^*_{\lambda,\perp}$ and $m^*_{\perp}$. This modification ensures that the resulting functions are locally insensitive to small perturbations in the nuisance parameter $\theta$ around $\theta_0$, thereby eliminating the first-order impact of parameter estimation on the test statistics.
We conclude this section with a final remark regarding kernel choice. For reproducing kernel-based tests, the selection of the kernel is critical for finite-sample performance. In the machine learning literature, particularly for nonparametric two-sample testing, considerable effort has been devoted to identifying optimal kernels liu2021learningdeepkernelsnonparametric,sutherland2021generativemodelsmodelcriticism,gretton2012optimal. These studies typically focus on maximizing the signal-to-noise ratio of the test statistic, leveraging the asymptotic properties of non-degenerate V- or U-statistics. While similar principles likely apply to model checking tests based on RKHS methods, the literature on optimal kernel selection in this context remains limited. Our KRR-based approach offers a practical alternative: kernel selection can be guided by minimizing the regression error using standard cross-validation procedures for KRR.
The remainder of the paper is organized as follows. Section (ref) reviews the relevant literature and situates our approach within the context of existing methods. In Section (ref), we present the main theoretical results, establishing the equivalence between the null hypothesis and the vanishing of $w^*_\lambda$, introducing the proposed test statistics, and deriving their asymptotic properties under the assumption that the nuisance parameter $\theta_0$ is known. Section (ref) extends these results to the more realistic setting where $\theta_0$ is unknown and must be estimated, and formally introduces the projection operator to eliminate the first-order effects of parameter estimation. Section (ref) provides simulation evidence on the finite-sample performance of the proposed tests, demonstrating that they are competitive with existing methods in moderate dimensions and outperform them as the covariate dimension increases. This section also includes an empirical application to the well-known National Supported Work (NSW) dataset. Section (ref) concludes. All proofs are provided in the Appendix.
Omnibus Tests for Model Checks: Among practitioners, the most widely used specification tests are M-tests, as introduced by newey1985generalized,newey1985maximum,tauchen1985diagnostic,wooldridge1990unified. These tests assess a finite set of unconditional moment restrictions implied by the conditional moment model. However, M-tests are inherently “directional”: they may fail to detect certain types of misspecification, as they only probe a limited set of directions in the space of alternatives.
Omnibus tests, by contrast, are designed to be consistent against any form of misspecification and are particularly valuable when the nature of potential model failure is unknown. There are two main approaches to constructing omnibus tests. The first compares the fitted parametric model to a nonparametric regression estimate, typically using smoothing techniques. Representative works in this category include eubank1990testing,hardle1993comparing,hong1995consistent,zheng1996consistent,ellison2000simple, among others. The second approach is based on integral transforms of the residuals, rather than the residuals themselves. Tests adopt this principle are often called the Integrated Conditional Moment (ICM) tests. Notable contributions in this direction include bierens1982consistent,bierens1997asymptotic,stute1997nonparametric,delgado2006consistent. For a comprehensive review of these two strands of the literature, see gonzalez2013updated.
Our work is closely related to the orthogonal series regression approach of eubank1990testing, who proposed a test statistic based on regressing the estimated residuals $\varepsilon_{\hat \theta}$ on an orthogonal basis expansion of the covariates. This yields a nonparametric estimator of $\mathbb{E}(\varepsilon_0|X)$. Their test can be interpreted as a joint F-test on the coefficients of the first $p_n$ terms in the series, with $p_n$ growing with the sample size. A limitation of this method is that it can only detect local alternatives of order $O(p_n^{1/4} / n^{1/2})$, which is slower than the $n^{-1/2}$ rate detectable by our proposed tests. However, our approach comes at the cost of non-pivotal null distributions, necessitating the use of multiplier bootstrap procedures for inference.
Our methodology is also related to the integral transformation approach. In fact, the RKHS is isometrically isomorphic to the space of square-integrable functions, and the RKHS inner product corresponds to the $L_2$ integral carrasco2007linear. For shift-invariant kernels, this connection is made explicit via Bochner's theorem (see Theorem 2.1 in muandet2020kernel for details).
Tests based on RKHS Tools: Classical ICM test statistics are formulated as V- or U-statistics using a Gaussian kernel with a fixed bandwidth (typically $0.5$). muandet2020kernel extended the ICM framework by allowing for general kernel choices, but their approach does not account for the estimation effect of nuisance parameters. escanciano2024gaussian addressed this limitation by employing a Gaussian process approach, deriving similar statistics but with a Neyman orthogonal kernel.
Our projection-based test statistics can be viewed as a modification of those proposed in these works. Specifically, existing kernel-based statistics assign weights to directional deviations according to the eigenvalues $\mu_i$ of the associated integral operator. In contrast, our approach assigns weights based on the ratios $\mu_i / (\mu_i + \lambda)$ or $\mu_i / (\mu_i + \lambda)^2$. This distinction is important, as it enables our tests to be more sensitive to deviations in high-frequency directions, potentially improving detection power in complex alternatives.
Tests Based on Random Locations: Our random location test statistics are closely related to approaches developed for nonparametric two-sample testing chwialkowski2015fast,jitkrittum2016interpretable and goodness-of-fit testing for distributional models jitkrittum2017linear. A key advantage of these methods, which carries over to our framework, is their interpretability: random location tests not only provide a powerful means of detecting model misspecification, but also offer intuitive, visual diagnostics that highlight where and how the model fails to capture the underlying data structure.
This section is organized into three parts. First, we establish the formal relationship between the population coefficient function $w^*_\lambda$ and the null hypothesis. Second, we analyze test statistics derived from projecting $w^*_\lambda$ onto functions in $\mathcal{H}_k$, examining two specific projections: (i) onto $w^*_\lambda$ itself and (ii) onto the mean embedding of the residual $m^* = \mathbb{E}(\varepsilon_0 k(X,\cdot))$. Third, we investigate random location test statistics, which leverage the idea that when $w^*_\lambda=0$, then under mild regularity conditions, $w^*_\lambda(y)=0$ almost surely for any $y$ sampled from a distribution with a Lebesgue density function.
Throughout this section, we assume the true parameters $\theta_0$ are known, and denote the true residual vector as $\boldsymbol{\varepsilon_0} = \boldsymbol{Y} - \mathcal{M}_{\theta_0}(X) \in \mathbb{R}^n$. In the next section, we will introduce the Neyman orthogonality approach that eliminates the estimation effects.
A central insight for model checking in the KRR framework is the equivalence between the null hypothesis $\mathbb{E}(\varepsilon_0 | X) = 0$ and the condition that the population coefficient function $w^*_\lambda = \boldsymbol{0}$ in $\mathcal{H}_k$. This equivalence underpins the use of its KRR estimator $\hat{w}$ as a diagnostic tool: if $\hat{w}$ is found to be significantly different from zero in the RKHS, it provides evidence against the adequacy of the specified model.
To formalize this, we introduce the following notation. Let $\boldsymbol{\Phi}$ and $\boldsymbol{K}$ be as defined in the Introduction. Denote by $\{\boldsymbol{u}_i\}_{i=1}^n$ and $\{\sigma^2_i\}_{i=1}^n$ the eigenvectors and eigenvalues of $\boldsymbol{K}$, respectively. We also define the integral operator $\boldsymbol{L}_k$ as \[ (\boldsymbol{L}_k f)(y) = \int_{\mathcal{X}} k(x,y) f(x) \, \mathrm{d} \mathbb{P}_X(x) \] where $\mathbb{P}_X$ is a probability measure on $\mathcal{X}$. Let $\{\phi_i\}_{i=1}^n$ and $\{\mu_i\}_{i=1}^n$ be the eigenfunctions and eigenvalues of $\boldsymbol{L}_k$, respectively. The eigenfunctions are orthonormal in $L_2(\mathbb{P}_X)$, the square-integral measurable function space with respect to the probability measure $\mathbb{P}_X$, and the eigenvalues are non-negative.
The vector $\boldsymbol{\Phi}$ admits the following singular value decomposition (SVD): \[ \boldsymbol{\Phi} = \sum_{i=1}^n \boldsymbol{u}_i \sigma_i (\sqrt{\mu_i} \phi_i) \] Note that $\{\sqrt{\mu_i} \phi_i\}_{i \geq 1}$ are orthonormal bases in $\mathcal{H}_k$. This SVD result can be verified by a sanity check:
here we use the RKHS inner product definition for $\phi_i$: \[ \langle \phi_i, \phi_i \rangle_{\mathcal{H}_k} = \frac{1}{\mu_i} \]
The KRR estimator $\hat{w}$ can be written as
To analyze $w^*_\lambda$, we assume the following conditions hold:
Assumption A. (i) The reproducing kernel is universal and bounded in the sense that $\sup_{x \in \mathcal{X}} k(x,x) < \infty$; (ii) the eigenvalues of $L_k$ satisfies $\sum_{i \geq 1}\mu_i <\infty$; (iii) $|\mathbb{E}(\varepsilon_0 \phi_i(X))| <M, \forall i \geq 1$; (iv) the kernel ridge regularization parameter satisfies $\lambda>0$ and is fixed.
Condition (i) is a standard assumption in the kernel literature. Universality implies that the kernel can approximate any continuous function on $\mathcal{X}$ (i.e., the RKHS is dense in the space of continuous functions). Boundedness ensures that each eigenfunction of $L_k$ is also bounded. Condition (ii) is needed to ensure that the RKHS $\mathcal{H}_k$ is well-defined and that the eigenvalues of $L_k$ decay sufficiently fast. This condition is satisfied by most commonly used kernels, including the Gaussian kernel $k(x, y) = \exp(-\gamma\|x - y\|^2 )$ for some $\gamma > 0$. Condition (iii) is a mild regularity condition on the residuals, ensuring that the projections onto the eigenfunctions are bounded. Condition (iv) would simplify our theoretical analysis. For the KRR estimator $\hat w$ to be consistent, the regularization parameter $\lambda$ should be shrinking to zero as the sample size $n$ increases at certain speed. However, we are not interested in the consistency of $\hat w$ in this paper, but rather in the “flatness” of this population coefficient function in the RKHS $\mathcal{H}_k$. Under the null, we expect $w^*_\lambda$ to be “flat” in the RKHS, i.e., $w^*_\lambda = \boldsymbol{0} \in \mathcal{H}_k$ (see Theorem 1). Fixing $\lambda$ simplifies the analysis and allows us to focus on the shape of the population coefficient function without the added complexity of a shrinking regularization parameter.
We now examine the structure of $w^*_\lambda$, which is provided in the following Lemma.
The parameters defining $w^*_\lambda$ are the expectations of the residuals projected onto the eigenfunctions $\phi_i$ of the operator $\boldsymbol{L}_k$. Each projection $\mathbb{E}(\varepsilon_0 \phi_i(X))$ is weighted by the factor $\mu_i / (\mu_i + \lambda)$, illustrating how the regularization parameter $\lambda$ controls the influence of each eigenfunction in the expansion of $w^*_\lambda$.
The null hypothesis $\mathbb{E}(\varepsilon_0 | X) = 0$ implies that the expected projection values $\mathbb{E}(\varepsilon_0 \phi_i(X))$ vanish for all $i \geq 1$. This leads to the conclusion that $w^*_\lambda = \boldsymbol{0} \in \mathcal{H}_k$ almost surely. Conversely, if $w^*_\lambda = \boldsymbol{0} \in \mathcal{H}_k$, then the expected projection values must also be zero, confirming the null hypothesis. This intuition is formalized in the following theorem, which establishes the equivalence between the null hypothesis and the condition that $w^*_\lambda = \boldsymbol{0} \in \mathcal{H}_k$.
Under the null, we expect \[ \langle w^*_\lambda, f \rangle_{\mathcal{H}_k} =0, \quad \forall f \in \mathcal{H}_k \] while under the alternative, \[ \langle w^*_\lambda, f \rangle_{\mathcal{H}_k} \neq 0, \quad \forall f \in \mathcal{H}_k, \text{ and } f \neq \boldsymbol{0} \] In this paper, we focus on two of these functions: (i) $f = w^*_\lambda$ and (ii) $f = m^* = \mathbb{E}(\varepsilon_0 k(X,\cdot))$. Two test statistics are then defined as
where \[ \hat m = \frac{1}{n}\sum_{i=1}^n \varepsilon_{0,i} k(X_i,\cdot) =\frac{1}{n} \boldsymbol{\varepsilon_0}^\top \boldsymbol{\Phi} \]
Using the SVD representation, we can rewrite the test statistics as
It is worthy noting that ICM test statistic of bierens1982consistent and its kernel conditional moment (KCM) extension muandet2020kernel take the form of
The derivations of the SVD representations above are provided in the proof of Theorem 2.
While these three test statistics share a similar structure—differing primarily in the weights applied to the terms $\{(1 / \sqrt{n}) \boldsymbol{\varepsilon_0}^\top \boldsymbol{u}_i\}_{i=1}^n$, which estimate the deviation signals $\{\mathbb{E}(\varepsilon_0 \phi_i(X))\}_{i=1}^n$ (see the proof of Lemma 1 for details)—these seemingly minor differences can have a substantial impact on testing power. To illustrate, consider the case of a Gaussian kernel, and suppose the deviation occurs in a specific direction: \[ \tilde \varepsilon_0 = \varepsilon_0 + \phi_l(x) \] then these two projection functions $w^*_\lambda$ and $m^*$ become
To compare the magnitude of deviations under different projections, consider the following ratios: \[ \frac{|| m^* ||_{\mathcal{H}_k}^2}{|| w^*_\lambda ||_{\mathcal{H}_k}^2} = (\mu_l+\lambda)^2, \qquad \frac{|| m^* ||_{\mathcal{H}_k}^2}{\langle w^*_\lambda, m^* \rangle_{\mathcal{H}_k}} = \frac{\langle w^*_\lambda, m^* \rangle_{\mathcal{H}_k}}{|| w^*_\lambda ||_{\mathcal{H}_k}^2} = \mu_l + \lambda. \] For a Gaussian kernel, it holds that $\sum_{i \geq 1} \mu_i = 1$\footnote{This follows from $1=\mathbb{E}(k(X,X)) = \sum_{i\geq 1}\mu_i \mathbb{E}(\phi_i(X)^2) = \sum_{i\geq 1}\mu_i$, since $\mathbb{E}(\phi_i(X)^2) = 1$ by orthonormality.}, so each $\mu_i \in (0,1)$ for all $i \geq 1$. Therefore, in this example, $|| w^*_\lambda ||_{\mathcal{H}_k}$ will exceed $|| m^* ||_{\mathcal{H}_k}$ whenever $\lambda < 1$. For the other two ratios, the comparison depends on the specific values of the eigenvalues $\mu_l$ and the regularization parameter $\lambda$.
Another important observation is that the deviation signals corresponding to $|| m^* ||_{\mathcal{H}_k}$ and $\langle w^*_\lambda, m^* \rangle_{\mathcal{H}_k}$, which are proportional to $\mu_i$ and $\mu_i / (\mu_i+\lambda)$, respectively, tend to decrease as the frequency index $i$ increases. In contrast, the deviation signal associated with $|| w^*_\lambda ||_{\mathcal{H}_k}$, given by $\mu_i / (\mu_i+\lambda)^2$, does not necessarily decrease monotonically with $i$; its behavior depends on the interplay between $\lambda$ and $\mu_i$.
These observations suggest that the choice of projection can significantly influence the power of the test, particularly in cases where the deviations are concentrated in specific directions. The next theorem, under the Assumption B, establishes the asymptotic distributions of these test statistics under both the null and fixed alternative hypotheses.
Assumption B. (i) The random variables $S=(Y,X)$ forms a strictly stationary process with probability measure $\mathbb{P}_{S}$; (ii) Reularity Conditions. (1) the residual function $\varepsilon: \mathcal{S} \times \Theta \longrightarrow \mathbb{R}$ is continuous on $\Theta$ for each $s \in \mathcal{S}$; (2) $\mathbb{E}(\varepsilon(S;\theta)|X=x)$ exists and is finite for every $\theta \in \Theta$ and $x \in \mathcal{X}$ for which $\mathbb{P}_X(x) >0$; (3) $\mathbb{E}(\varepsilon(S;\theta)|X=x)$ is continuous on $\Theta$ for all $x \in \mathcal{X}$ for which $\mathbb{P}_X(x) >0$.
These conditions are standard in the literature for conditional moment models and model checks, see muandet2020kernel,hall2003generalized for example.
We now discuss the asymptotic power property of the proposed test statistics. We consider the local (Pitman) alternatives: \[ H_{1,n}: Y - \mathcal{M}_{\theta_0}(X) = \boldsymbol{\varepsilon_0} + \frac{R(X)}{\sqrt{n}} a.s., \] where $\mathbb{E}(\varepsilon_0 |X) =0$ and $R(X)$ is a nonzero square integrable measurable function of $X$ with respect to $\mathbb{P}_X$. The asymptotic result under $H_{1,n}$ is given by
The population coefficient function $w^*_\lambda$ acts as a “witness” to deviations from the null hypothesis, such that $|w^*_\lambda(v)|$ becomes large when there is a discrepancy in the region around some point $v \in \mathcal{X}$. The projection-based test statistics discussed earlier quantify the overall “flatness” of $w^*_\lambda$ using the RKHS inner product, providing a global measure of model adequacy. In contrast, the random location test statistics soon to be introduced are designed to assess the “flatness” of $w^*_\lambda$ in a more localized and interpretable way, enabling the detection and visualization of deviations at specific regions in the input space.
A central idea behind random location test statistics is to leverage the analytic properties of the kernel function $k$, such as those of the Gaussian kernel sun2008reproducing, which ensure that all functions in the corresponding RKHS is analytic. Hence, $w^*_\lambda$ is a real analytic function. For real analytic functions, if $w^*_\lambda \neq \boldsymbol{0}$, the set of points where $w^*_\lambda(v) = 0$ has Lebesgue measure zero. This statement is formalized as the following lemma:
To formally ensure that $w^*_\lambda$ is a real analytic function in $\mathcal{H}_k$, we require the following assumption:
Assumption C. (i) The kernel $k$ is shift-invariant and Mercer; (ii) $k(x, y)$ is real analytic in both $x$ and $y$.
A Mercer kernel is continuous, symmetric, and yields a positive semidefinite kernel matrix for any finite set of points. A kernel is shift-invariant if $k(x, y) = k(x - y)$ for all $x, y \in \mathcal{X}$. The Gaussian kernel is a prototypical example that satisfies both conditions in Assumption C.
The following Lemma implies that $w^*_\lambda$ is real analytic in $\mathcal{H}_k$ under Assumption C.
Lemmas (ref) and (ref) imply that, by evaluating $w^*_\lambda$ at a finite set of location points $\boldsymbol{V} = \{v_1, \ldots, v_J\}$ sampled i.i.d. from a distribution with a Lebesgue density $\eta$, we can almost surely detect any nonzero deviation: if $w^*_\lambda \neq \boldsymbol{0}$, then $\{w^*_\lambda(v_j)\}_{j=1}^J$ will be nonzero with probability one.
The following two test statistics are developped based on this simple idea:
where $$\boldsymbol{k}(v_j) = \left(k(x_1,v_j), \ldots, k(x_n,v_j)\right)^\top$$ is the vector of kernel functions evaluated at $\{(x_i,v_j)\}_{i=1}^n$.
We now establish the asymptotic distributions of the random location test statistics under both the null and fixed alternative hypotheses.
The asymptotic properties of the random location test statistics under the local alternative $H_{1,n}$ mirror those of the projection-based statistics. The following theorem summarizes these results.
To demonstrate the interpretability and visualization capabilities of the coefficient function (and the random location test statistics), consider the following three data-generating processes (DGPs):
where $\beta$ is a $2 \times 1$ vector of ones, $X$ is a $2 \times 1$ vector with entries independently drawn from the standard normal distribution, and $e$ is an independent standard normal random variable. DGP0 represents the null model, while DGP1 and DGP2 correspond to alternative models.
For these examples, we use the true residuals: $\varepsilon_{0} = e$ for DGP0, $\varepsilon_{0} = 4.5 (X^\top \beta)^2 + e$ for DGP1, and $\varepsilon_{0} = 4.5 \exp\left(- (X^\top \beta)^2 \right) + e$ for DGP2. The Gaussian kernel is employed, with both the kernel parameter $\gamma$ and the regularization parameter $\lambda$ selected via cross-validation. Figure (ref) displays the estimated coefficient function for these three DGPs.
Under the null model (DGP0), values of the coefficient function are close to zero, indicating no evidence against the model. In contrast, under the alternative models (DGP1 and DGP2), the coefficient function values are substantially larger, reflecting the presence of model misspecification. Furthermore, the spatial patterns of the coefficient function align with the structure of the underlying alternatives: for DGP1, larger values are observed in regions where both components of $X$ have large magnitudes and the same sign, while for DGP2, the coefficient function peaks near regions where $X^\top \beta$ is close to zero. This demonstrates that the coefficient function not only detect deviations from the null but also provide interpretable insights into where the model fails.
So far, we have assumed that the value of $\theta_0$ is known. In practice, $\theta_0$ is estimated by a consistent estimator $\hat{\theta}$. In this section, we discuss how to deal with the estimation effect when $\theta_0$ is estimated.
We adopt the approach advocated by escanciano2014specification,escanciano2024gaussian,sant2019specification, which employs Neyman orthogonality projections to remove the impact of parameter estimation. This idea has been further extended by sancetta2022testing, who consider subspace restrictions with high-dimensional nuisance parameters estimated via penalized methods in an RKHS. In contrast to the conventional wild bootstrap methods used in, for example, delgado2006consistent, this projection-based approach does not depend on a linear Bahadur representation for the estimator of $\theta_0$, nor does it require re-estimating parameters within each bootstrap sample.
The following assumptions are required for this purpose.
Assumption D. (i) The parameter space $\Theta$ is a compact subset of $\mathbb{R}^d$; (ii) the true parameter $\theta_0$ is an interior point of $\Theta$; and (iii) the consistent estimator $\hat{\theta}$ satisfies $\lVert \hat{\theta} - \theta_0 \rVert = O_p (n^{-\alpha})$, with $\alpha > 1/4$.
Assumption D is weaker than related conditions in the literature. We only impose $\hat{\theta}$ converges in probability at a slower rate than usual. In addition, we also do not require it to admit an asymptotically linear representation. This could be useful in the context of non-standard estimation procedures, such as the LASSO.
Additional regular conditions on the smoothness of the residual function is also required.
Assumption E. (i) The residual $\varepsilon(s;\theta)$ is twice continuously differentiable with respect to $\theta$, with its first derivative $g_{\theta}(x) = \mathbb{E}(\nabla_{\theta} \varepsilon(s;\theta)|X=x)$ satisfying $\mathbb{E}\left(\sup_{\theta \in \Theta} \lVert g_{\theta}(X) \rVert\right) < \infty$ and its second derivative satisfying $\mathbb{E}\left(\sup_{\theta \in \Theta} \lVert \nabla g_{\theta}(X) \rVert\right) < \infty$; (ii) the matrix $\Gamma_{\theta} = \mathbb{E}\left[g_{\theta}(X) g_{\theta}(X) ^\top\right]$ is nonsigular in a neighborhood of $\theta_0$.
Under Assumptions D and E, we now introduce a projection operator $\boldsymbol{\Pi}$, acting on a random variable (or random vector) $\omega(S)$ and its realization $\omega(s)$, defined as follows: \[ \boldsymbol{\Pi} \omega(s) = \omega(s) - \mathbb{E}(\omega(S) g_{\theta_0}^\top(X)) \Gamma_{\theta_0}^{-1} g_{\theta_0}(x) \] Applying this projection operator to the residual function leads to modified versions of $w^*_\lambda$ and $m^*$:
To analyze the local behavior of these functions in a neighborhood of $\theta_0$, consider their derivatives with respect to $\theta$ evaluated at $\theta_0$:
and
This establishes that the model checks are locally robust to small perturbations in $\theta$ around $\theta_0$. The vanishing derivatives follow from the dominated convergence theorem and the fact that
The matrix estimator of this projection operator is given by \[ \hat{\boldsymbol{\Pi}} = \boldsymbol{I}_n - \mathbb{G}\left(\mathbb{G}^\top \mathbb{G}\right)^{-1} \mathbb{G}^\top \] where $\mathbb{G}$ is a $n \times d$ matrix of scores whose $i$th row is given by $\hat g_{i}^\top = (\nabla_{\theta} \varepsilon(s_i;\theta)|_{\theta = \hat \theta})^\top$, and $\boldsymbol{I}_n$ is the $n \times n$ identity matrix.
The following lemma states how this projection operator eliminates the estimation effect in finite sample.
As long as the convergence speed satisfies $\alpha>1/4$, we have
Here $\boldsymbol{o_p(1)}$ is a function in $\mathcal{H}_k$ whose RKHS norm converges to zero in probability.
Similarly,
Thus, we can replace $\boldsymbol{\varepsilon}_{\hat \theta}$ with $\boldsymbol{\varepsilon_0}$ in the test statistics, leading to
The same projection can be applied to the random location test statistics:
and
As a result, all asymptotic results established in the previous sections continue to hold after applying the projection, with the following modifications: in the null and local alternative distributions, $Z_i$ is replaced by $Z_{i,\perp}$, and $\mathbb{E}(R(X)\phi_i(X))$ is replaced by $\mathbb{E}(\boldsymbol{\Pi} R(X)\phi_i(X))$ in the local alternative case. Here, $Z_{i,\perp} \sim \mathcal{N}(0, S_{i,\perp}^2)$ are independent normal random variables with variance $S_{i,\perp}^2 = \mathrm{Var}((\boldsymbol{\Pi} \varepsilon_0) \phi_i(X))$.
We employ the multiplier bootstrap to approximate the null distributions of the test statistics. For theoretical justification, we use the notion of almost surely (a.s.) consistency, denoted by $\overset{d^*}{\longrightarrow}$; see Chapter 2.9 in vaart1997weak. The general idea is to generate a sequence of i.i.d. random variables $\{v_i\}_{i=1}^n$ with zero mean, unit variance, bounded support, and independent of the data $\{s_i\}_{i=1}^n$. These are used to construct the bootstrap sample $\{\varepsilon(s_i;\hat \theta) v_i\}_{i=1}^n$. A classical choice for such random variables is the Mammen two-point distribution: \[ \mathbb{P}(V_i = 0.5(1-\sqrt{5}))=b, \quad \mathbb{P}(V_i = 0.5(1+\sqrt{5}))=1-b, \] where $b = (1+\sqrt{5})/(2\sqrt{5})$; see mammen1993bootstrap for details.
Let $\mathbf{a} \odot \mathbf{b}$ denote the element-wise (Hadamard) product. The bootstrap analogues of our test statistics are then given by
where $\boldsymbol{V} = (v_1, \ldots, v_n)^\top$.
We consider the following simulation settings. The null model is given by \[ \mathrm{DGP}_0: \quad Y =\alpha+ X^\top \beta + e \] The fixed alternative models under homoskedasticity condition are given by:
All of the above DGPs (from $0$ to $4$) share the same covariate and error structure: $X$ is a $d \times 1$ vector with entries independently drawn from the standard normal distribution, and $e$ is an independent standard normal random variable. The true parameter vector $\beta$ is of dimension $d \times 1$ with each entry equal to $0.5$, and the intercept term is set to $\alpha = 1$. In the simulation studies, we consider $d = 10$ and $d=20$.
For the fixed alternative models under heteroskedastic condition, we consider the following settings:
For $d=10$, the covariates $\{X_{l}\}_{l=1}^{5}$ are independently drawn from uniform distributions on $[0,\, 1+0.1(l-1)]$, while $\{X_{l}\}_{l=6}^{10}$ are drawn from normal distributions with mean $0$ and standard deviation $1+0.1(l-5)$. For $d=20$, $\{X_{l}\}_{l=1}^{10}$ are drawn from uniform distributions on $[0,\, 1+0.1(l-1)]$, and $\{X_{l}\}_{l=11}^{20}$ are drawn from normal distributions with mean $0$ and standard deviation $1+0.1(l-10)$. The heteroskedasticity functions are defined as $c_1(X) = |X^\top \boldsymbol{1}|$ and $c_2(X) = |X^\top \boldsymbol{1}|$, where $\boldsymbol{1}$ is a $d \times 1$ vector of ones. In both cases, $\beta = \boldsymbol{1}$, $\alpha = 1$, and the error term $e$ is standard normal.
For local alternatives, we consider the following models:
For $d=10$, the covariates $\{X_{l}\}_{l=1}^{5}$ are independently drawn from uniform distributions on $[0,\, l]$, while $\{X_{l}\}_{l=6}^{10}$ are drawn from normal distributions with mean $0$ and standard deviation $1+0.1(l-5)$. For $d=20$, $\{X_{l}\}_{l=1}^{10}$ are drawn from uniform distributions on $[0,\, l]$, and $\{X_{l}\}_{l=11}^{20}$ are drawn from normal distributions with mean $0$ and standard deviation $1+0.1(l-10)$. The heteroskedasticity functions are defined as $d_1(X) = \sqrt{0.1 + \sum_{l=1}^{5} X_{l} + \sum_{l=6}^{10} X_{l}^2}$ and $d_2(X) = \sqrt{0.1 + \sum_{l=1}^{5} X_{l} + \sum_{l=6}^{20} X_{l}^2}$. In both cases, $\beta$ is a vector of ones, $\alpha = 1$, and the error term $e$ is standard normal.
For all models, we use the ordinary least squares (OLS) estimator $\hat \beta$ to estimate the parameter $\beta$.
In addition to the proposed KRR-based test statistics ($\hat T_{\mathrm{proj},\perp}^{(1)}$, $\hat T_{\mathrm{proj},\perp}^{(2)}$, $\hat T_{\mathrm{rand},\perp}^{(1)}$, and $\hat T_{\mathrm{rand},\perp}^{(2)}$), we also include two benchmark methods for comparison: the integrated conditional moment (ICM) test statistic ($\hat T_{\mathrm{ICM}}$) of bierens1982consistent, and the Gaussian process (GP) test statistic ($\hat T_{\mathrm{GP}}$) of escanciano2014specification. Both alternatives have the quadratic form: \[ n \hat T = \frac{1}{n} \boldsymbol{\varepsilon}_{\hat \theta}^\top \boldsymbol{K} \boldsymbol{\varepsilon}_{\hat \theta} \] The ICM test statistic employs the wild bootstrap for inference, following delgado2006consistent, while the GP test statistic uses the multiplier bootstrap as in escanciano2024gaussian, with a modified orthogonal kernel $\boldsymbol{K}_{\perp} = \hat{\boldsymbol{\Pi}}^\top \boldsymbol{K} \hat{\boldsymbol{\Pi}}$. We set the number of bootstrap replications to $B = 500$, and each simulation scenario is repeated $R = 1000$ times. Empirical sizes and powers are calculated as the proportion of rejections over the $R$ replications. The significance level is fixed at $5\%$.
All test statistics are implemented using the Gaussian kernel $k(x, y) = \exp(-\gamma\|x - y\|_2^2)$. For the ICM test, the kernel parameter is fixed at $\gamma = 0.5$. For the GP test, following escanciano2024gaussian, we set $\gamma = 1 / \mathrm{median}(\{\|x_i-x_j\|_2\}_{i\neq j})$. For the KRR-based test statistics, both the kernel parameter and the regularization parameter $\lambda$ are selected via 5-fold cross-validation.
Ideally, one would split the data into training, validation, and testing sets, using the first two to tune parameters and the last for testing. In practice, we find that splitting the data into training and validation sets for parameter selection, and then applying the chosen parameters to the full dataset, yields reliable results.
A distinctive feature of the random location test statistics is the choice of the number of random locations, denoted by $J$. In the main simulation studies, we set $J = 3$. To assess the robustness of the results, we also examine the impact of varying $J$ by considering values from $J = 1$ up to $J = 15$. Theoretically, the random locations can be sampled from any distribution with a Lebesgue density supported on the domain of $X$. In practice, we find that sampling from a multivariate normal distribution fitted to the observed data provides reliable performance.
Tables (ref) and (ref) summarize the empirical sizes and powers of the proposed KRR-based test statistics and benchmark methods for dimensions $d=10$ and $d=20$, respectively. Each table reports results for two sample sizes ($n=200$ and $n=400$). The "SIZE" rows correspond to the null DGP ($DGP_0$), showing the empirical type I error rates. The "POWER" rows correspond to alternative DGPs ($DGP_1$ to $DGP_6$), showing the empirical power of each method. The best-performing method for each scenario is highlighted in bold.
For both dimensions, the KRR-based projection tests ($\hat T_{\mathrm{proj},\perp}^{(1)}$, $\hat T_{\mathrm{proj},\perp}^{(2)}$) and random location tests ($\hat T_{\mathrm{rand},\perp}^{(1)}$, $\hat T_{\mathrm{rand},\perp}^{(2)}$) maintain empirical sizes close to the nominal 5% level, while the benchmark GP and ICM tests tend to be conservative, especially as $d$ increases.
In terms of power, the projection-based KRR tests are most effective for alternatives with strong global nonlinear deviations (e.g., $DGP_1$, $DGP_3$, $DGP_4$), with $\hat T_{\mathrm{proj},\perp}^{(2)}$ often achieving the highest power. The random location tests excel in detecting alternatives with localized or heteroskedastic effects (e.g., $DGP_2$, $DGP_5$, $DGP_6$), and their power increases with sample size. Notably, for $d=20$, the benchmark methods show little to no power, while the KRR-based tests retain substantial power across all alternatives. These results highlight the flexibility and robustness of the proposed KRR-based tests in a variety of challenging scenarios.
To further assess the robustness of the random location tests, we conduct additional simulations varying the number of random locations ($J$). Figure (ref) summarizes these results. The left panel displays the power of the random location tests as a function of $J$ for $\mathrm{DGP}_2$ with $d=10$, while the right panel presents the corresponding results for $\mathrm{DGP}_6$ with $d=20$. The findings demonstrate that the power of the random location tests remains stable across a wide range of $J$ values, indicating that the performance of these tests is not overly sensitive to the specific choice of the number of random locations.
The impact of training programs on labor market outcomes has been a longstanding focus in economics. To rigorously assess such effects, the National Supported Work (NSW) Demonstration was conducted from 1975 to 1979, funded by both private and federal sources. The program randomly assigned participants at various sites across the United States to either a treatment group (receiving supported work interventions) or a control group.
The NSW dataset has become a standard benchmark in the causal inference literature. The original analysis by lalonde1986evaluating and subsequent studies using propensity score methods, such as dehejia1999causal, have provided valuable insights into the effectiveness of the program. In this paper, we use the version of the dataset compiled by Dehejia and Wahba, which consists of 445 observations: 185 in the treatment group and 260 in the control group. The dataset includes covariates such as age, education, and prior earnings, with the primary outcome being real earnings in 1978. The data are publicly available at \url{https://users.nber.org/ rdehejia/nswdata2.html}.
To mitigate numerical instability, we rescale the “age” and “education” variables by dividing them by 10, and apply a logarithmic transformation to all earnings variables. The final dataset comprises 8 covariates: 4 continuous variables and 4 binary indicator (dummy) variables.
We aim to test both the specification of the propensity score model and the hypothesis of a zero Conditional Average Treatment Effect (CATE), both individually and jointly. Let $Y(1)$ and $Y(0)$ denote the potential outcomes under treatment and control, respectively, and let $X$ be a vector of observed covariates. The treatment indicator is $T$, so the observed outcome is $Y = Y(1)T + Y(0)(1-T)$.
The propensity score is modeled as $P(T=1|X) = \Phi(\beta^\top X)$, where $\Phi$ denotes the standard normal cumulative distribution function. The CATE function is defined as $\tau(X) = \mathbb{E}[Y(1) - Y(0) \mid X]$. Under standard unconfoundedness assumptions (i.e., $Y(1)$ and $Y(0)$ are independent of $T$ given $X$), and if the propensity score model is correctly specified, it follows that \[ \mathbb{E}\left( \frac{Y(T-\Phi(\beta^\top X))}{\Phi(\beta^\top X)(1-\Phi(\beta^\top X))} \,\middle|\, X \right) = \tau(X) \quad \text{a.s.} \] The residual for testing the propensity score model alone is \[ \varepsilon_0 = T - \Phi(\beta^\top X). \] For testing the joint hypothesis of correct propensity score specification and zero CATE (i.e., $\tau(X) = 0$), the residual vector is \[ \varepsilon_0 =
. \] For joint testing, the test statistic is constructed by first computing the individual test statistics for each component of the residual vector $\varepsilon_0$, and then summing these statistics to obtain the overall test statistic.
We report results for the test statistics discussed in the simulation studies, replacing the ICM test statistic with the GP test using a kernel parameter of $\gamma = 0.5$. This substitution is made because implementing the ICM test with a wild bootstrap is not straightforward in the context of a probit model. The null distributions of the test statistics are approximated using the multiplier bootstrap method, with $B=500$ bootstrap replications. The results are presented in Table (ref).
The results show that none of the test statistics reject the null hypothesis of correct propensity score specification, as all $p$-values are well above conventional significance levels. In contrast, all test statistics strongly reject the joint null hypothesis of both correct propensity score specification and zero CATE, with $p$-values equal to zero. This suggests that the probit model provides a good fit, and there is strong evidence against the hypothesis of no treatment effect.
In this paper, we introduce two new classes of reproducing kernel-based tests for model checking in conditional moment restriction models. The central idea is to perform kernel ridge regression (KRR) of the estimated residuals on kernel functions indexed by the observed data. The population analogue of the KRR estimator is a function in the RKHS, which serves as a diagnostic metric for model adequacy.
We introduce two classes of tests. The first class, projection-based tests, involves projecting the KRR estimator onto (i) itself and (ii) the mean embedding estimator of the residuals. The second class, random location tests, evaluates the KRR estimator at $J$ randomly chosen points drawn from any distribution with a Lebesgue density, assuming the kernel is analytic. In practice, sampling these locations from a multivariate normal distribution fitted to the observed data yields good results.
We establish the asymptotic properties of the proposed tests under the null, fixed alternatives, and local alternatives. The tests are consistent against fixed alternatives and can detect local alternatives at the optimal $n^{-1/2}$ rate.
Extensive simulation studies demonstrate that the proposed tests maintain nominal size and exhibit strong power across a range of alternatives. In particular, they outperform existing methods in both size control and power, especially as the covariate dimension increases. An empirical application to the National Supported Work dataset further illustrates the practical utility and effectiveness of the proposed approach.