EconBase
← Back to paper

Wald inference on varying coefficients

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.

108,263 characters · 21 sections · 82 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Wald inference on varying coefficients

\affil[1]{Department of Economics, Queen's University, Dunning Hall, 94 University Avenue, Kingston, Ontario K7L 3N6, Canada, and Department of Economics, University of Essex, Wivenhoe Park, Colchester, CO4 3SQ, UK.} \affil[2]{Department of Economics, Antai College of Economics and Management, Shanghai Jiao Tong University, 1954 Huashan Road, Shanghai, 200030, China PRC. } \affil[3]{Department of Economics, National University of Singapore, 1 Arts Link, 117570, Singapore.} \affil[4]{International Business School, Shanghai University of International Business and Economics, Shanghai, 201620, China PRC.}

abstractWe present simple to implement Wald-type statistics that deliver a general nonparametric inference theory for linear restrictions on varying coefficients in a range of regression models allowing for cross-sectional or spatial dependence. We provide a general central limit theorem that covers a broad range of error spatial dependence structures, allows for a degree of misspecification robustness via nonparametric spatial weights and permits inference on both varying regression and spatial dependence parameters. Using our method, we first uncover evidence of constant returns to scale in the Chinese nonmetal mineral industry’s production function, and then show that Boston house prices respond nonlinearly to proximity to employment centers. A simulation study confirms that our tests perform very well in finite samples.

{\it Keywords: Spatial autoregression, varying coefficients, inference}\\ {\it JEL classification: C21, C31}

\spacing{1.8}

Introduction

This paper develops a general framework for inference on varying coefficient regression models that is as easy to implement in practice as familiar Wald tests for finite-dimensional parameters. We also allow for spatial autoregressive (SAR) structure to permit interaction across the cross-section of economic agents. SAR models are a popular and parsimonious tool for researchers working in settings where economic agents interact via economic or social links. Such links can be geographic proximity, or proximity in some more general sense such as a social network. The SAR model can accommodate these types of general links easily, and this goes some way to explain its appeal, see e.g. Helmers2014 and hsieh2018 for examples of applications to peer effects. No surprise then that the baseline SAR model, due to cliff1981spatial, has become much studied by econometricians. This model has long been popular with regional scientists, see anselin1988spatial for an early book length treatment; important econometric contributions that established the foundation for rigorous econometric analysis of such models include kelejian1998generalized, lee2004asymptotic and Kuersteiner2020.

In keeping with the growing importance of SAR models and their ability to parsimoniously control for cross-sectional dependence in multiple regression designs, various strands of the literature have explored how to widen their scope. One such is the focus on varying coefficient SAR models, wherein model parameters are allowed to be unknown functions of some observed economic variable(s). This permits a flexible effect of both the regressors and spatial dependence on the outcome variable, and indeed varying coefficient models are abundant in statistics and econometrics (see e.g. Fan1999,Fan2008). In the SAR context, varying coefficient models have been studied by Sun2014 and a sequence of important papers by Sun2016, malikov2017semiparametric, Sun2018 and Sun2024semiparametric, for example. More generally, semiparametric SAR models have attracted much interest, see e.g. Su2010, Su2012, Robinson2012c and Zhang2013, but their focus tends to be on inference on the parametric component.

In this paper, we show how inference on varying coefficients in a range of SAR models can be conducted just like familiar parametric Wald tests for linear restrictions. Standard varying coefficient regression without cross-sectional dependence is naturally covered as a special case. The key idea is that if series approximations are used for the nonparametric varying coefficients then the series coefficients can be employed for very easy practical inference. This is motivated by the scientific aim of simplification of nonparametric/semiparametric inference. Indeed, a growing literature stresses such ideas, see e.g Ackerberg2012, Gupta2018c and Korolev2019, amongst others. Naturally our methods are also applicable to the usual varying coefficient regression model with no cross-sectional dependence, see e.g. Ahmad2005, which is a simply a particular case of our theory.

To explain the intuition of our approach, suppose a sample of $n$ observations is available. With our method, the researcher can reduce the econometric problem of inference on an infinite-dimensional object to inference on a growing number of series coefficients, say $p$. In practice, this simply means writing down the usual Wald test statistic $\mathscr{W}$ for $p$ linear restrictions, observing that $p\rightarrow\infty$ as $n\rightarrow\infty$ if it is the length of a series approximation, and appealing to a theorem that establishes $(\mathscr{W}-p)/\sqrt{2p}\overset{d}{\rightarrow} N(0,1)$ as $n\rightarrow\infty$, under the null hypothesis. The idea stems from the fact that as $p\rightarrow\infty$, a $(\chi^2_p-p)/\sqrt{2p}$ random variable approaches a standard normal variate, see e.g. DeJong1994, Hong1995, Gupta2018c and Gupta2023.

Using this series expansion approach, we theoretically justify Wald-type test statistics for inference on varying coefficients. We first cover a baseline high-order SAR model with varying regression coefficients but constant spatial lag parameters and extend the model to allow for spatial error dependence. The baseline model is similar to the one estimated by Sun2014, but also features a growing number of spatial lags and covariates. We then show how to allow for a degree of misspecification robustness by incorporating nonparametric spatial weights \`a la pinkse2002. Our approach is also shown to work for models with varying coefficients on the spatial lags, in the spirit of malikov2017semiparametric.

We then apply our methods in two different settings where varying coefficients have been shown to play an important role. In the first application, we test for constant returns to scale (CRS) in the production function of the Chinese nonmetal mineral industry. This re-visits the empirical study in li2002semiparametric, who showed that output elasticities of capital and labor vary with managerial expense. They also noted that the returns to scale, which is the sum of output elasticities of inputs, are close to one at many managerial expense levels suggesting CRS, but they did not have a test for this hypothesis. Using a newer dataset, our methods can test for CRS while allowing for spatial dependence. We do not reject the CRS hypothesis in all settings, which suggests that CRS technology is a salient feature of this industry.

Our second study re-examines the Boston house price data of harrison1978hedonic. We estimate a model similar to Sun2014, who used a model selection procedure to recommend that location should have a varying effect on house prices with respect to some other features of the property. We formally test for the varying location effects of these features using our statistics and indeed find them to be statistically significant.

In Monte Carlo simulation studies, we experiment with a range of SAR models and spatial dependence structures to demonstrate the applicability and generality of our approach. Results show that our tests have excellent finite-sample performances.

The rest of the paper is as follows: Section (ref) introduces our method in a baseline regression model allowing for higher-order SAR structure. Section (ref) extends the model to allow for spatial error dependence of the Kelejian2007 type. Section (ref) allows for nonparametric spatial weights, as in pinkse2002, and thereby permits some degree of weight matrix robustness. Section (ref) shows how our approach can be adapted to the setting where the spatial lag coefficient varies, like in malikov2017semiparametric. Section (ref) applies our methods to the Chinese nonmetal mineral industry and Boston house prices. Section (ref) presents the simulation study, and Section (ref) concludes. Proofs and additional simulation results are in the online appendix.

Baseline higher-order spatial autoregression

\setlength\abovedisplayskip{6pt} \setlength\belowdisplayskip{6pt} Consider a vector of unknown functions $\delta\left(\cdot\right)=\left(\delta_{1}\left(\cdot\right),\ldots,\delta_{d_{\delta}}\left(\cdot\right)\right)'$, and the model

equation[equation omitted — 159 chars of source]

where $y_n$ is the $n\times 1$ vector with typical element $y_{in}$, $w_{in,j}'$ is the $i$-th row of the spatial weight matrix $W_{jn}$, $j=1,\ldots,d_\lambda$, $d_\lambda\rightarrow\infty$ as $n\rightarrow\infty$, and $\epsilon_{in}$ is $i.i.d.$ with mean 0 and unit variance. Also, $\left(x_{in},p_{in},z_{in}\right)\in\mathcal{X}\times \mathcal{P}\times\mathcal{Z} \subseteq\mathbb{R}^{d_{\beta} \times d_{\delta} \times d_z}$, and $d_{\beta}\rightarrow\infty$ as $n\rightarrow\infty$ but $d_z$ and $d_\delta$ are fixed. For all $i=1,\ldots,n$, $z_{in}$ is throughout the paper uncorrelated with $\epsilon_{in}$ while $x_{in}$ is allowed to be correlated with $\epsilon_{in}$.

It is worth noting that our results also hold when $d_\lambda$ and $d_\beta$ are held fixed but we state our theorems for the more complex case where they diverge. Allowing $d_\lambda$ and $d_\beta$ to diverge with sample size allows a flexible modeling approach, which has been studied in statistics since at least the work of huber1973robust, but in spirit the sequence of experiments considered by LeCam1960 provides an even earlier reference. Econometricians have also studied such models in least squares and GMM settings, see e.g. Andrews1985, Koenker1999 and Cattaneo2018, to name a few. In a SAR setting, Gupta2013 and Gupta2018c emphasize the advantages of this approach, especially when clustered data imply asymptotic regimes where $d_\lambda$ diverges with $n$.

We now drop $n$ subscripting, but occasionally remind the reader of the dependence on sample size of various quantities. We take spatial weight matrices $W_j$ to be non-stochastic everywhere except in Section (ref). These weight matrices are also uniformly bounded in row and column sums, a commonly used restriction to control spatial dependence kelejian1998generalized, lee2004asymptotic. This property is assumed to hold almost surely in Section (ref), where the spatial weights are stochastic.

Our aim is to do inference on $\delta(\cdot)$, specifically we wish to test a fixed number $m$ restrictions of the form $S\delta(z)=s, z\in\mathcal{Z}$, where $S$ is a known, constant $m\times d_\delta$ matrix and $s$ is a known, constant $m\times 1$ vector. Because we can convert any linear restriction to an exclusion restriction, we will focus on tests of the null hypothesis

equation[equation omitted — 89 chars of source]

Our test statistics will approximate $\delta(\cdot)$ by a series expansion and test if the coefficients in the expansion are jointly zero, thus converting the restriction in $H_0^{true}$ to an increasing-dimensional one, up to some suitably negligible approximation error.

Test statistic

We approximate $\delta_{k}(z)$ by $\psi_{k}^{h_k}(z)'\alpha_{k}^{h_k}$, where $\psi_{k}^{h_{k}}(z)=\left(\psi_{k1}(z),\ldots,\psi_{kh_k}(z)\right)'$ and $\alpha_{k}^{h_{k}}=\left(\alpha_{k1},\ldots,\alpha_{kh_k}\right)'$, for some basis functions $\psi_{k\ell}(\cdot)$, $\ell=1,\ldots,h_k$, $k=1,\ldots,d_{\delta}$. Thus each $\delta_{k}(z)$ is approximated by a linear combination of $h_k$ basis functions $\psi_{k}^{h_{k}}$ with coefficients $\alpha_{k}^{h_{k}}$. Define the $d_\alpha\times 1$ vector $\psi_{i}\equiv\psi_{i}\left(p_{i},z_{i}\right)=\left(p_{i1}\psi_{1}^{h_1}\left(z_{i}\right)',\ldots,p_{i d_{\delta}}\psi_{d_{\delta}}^{h_{d_{\delta}}}\left(z_{i}\right)'\right)'$ and $\alpha=\left(\alpha_{1}^{h_1}{'},\ldots,\alpha_{d_{\delta}}^{h_{d_{\delta}}}{'}\right)'$, with $d_\alpha=\sum_{k=1}^{d_{\delta}}h_k$. Then we can write ((ref)) as

equation[equation omitted — 136 chars of source]

with $u_i=r_{i}+\epsilon_i$, where $r_{i}=p_{i}'\delta\left(z_{i}\right)-\psi_{i}'\alpha$ is the approximation error. Let $X$ and $\Psi$ be matrices with typical rows $x_i'$ and $\psi_i'$, and $u$ be with elements $u_i$. Writing $Q$ for the $n\times d_\lambda$ matrix with typical columns $W_jy$, (ref) can be written in matrix notation as

equation[equation omitted — 86 chars of source]

where $L=[Q, X,\Psi]$ and $\xi=\left(\lambda',\beta',\alpha'\right)'$. The null hypothesis to approximate $H_0^{true}$ is

equation[equation omitted — 46 chars of source]

and our test will be implemented using the 2-stage least squares (2SLS) estimator

equation[equation omitted — 144 chars of source]

where $P_{K_{1}}=K_{1}\left(K_{1}'K_{1}\right)^{-1}K_{1}'$ and $K_{1}$ is an $n\times J_{1}$ instrument matrix with $J_{1}\geq d_{\xi}=d_{\lambda}+d_{\beta}+d_{\alpha}$ but the same asymptotic order, i.e. ${J_1}/{d_{\xi}}$ tends to a constant at least unity, as $n\rightarrow \infty$. Let $T_{1}=\frac{1}{n} L'{P_{K_{1}}}L$ and $\mathscr{D}_{1}=\frac{1}{n} RT_{1}^{-1}R'$, where $R=\left[0_{d_\alpha\times\left(d_\lambda+d_\beta\right)}, I_{d_\alpha}\right]$, and $I_{d_\alpha}$ is the $d_\alpha\times d_\alpha$ identity matrix. Then, the test statistic is

equation[equation omitted — 123 chars of source]

Asymptotic properties

We begin with restrictions on moments of various objects, noting that throughout the paper $C$ denotes a generic positive constant, arbitrarily large but independent of $n$. Let $k_{ri}$ be elements of the instrumental variable matrix, corresponding to $K_{1}, K_{2}$ and $K_{3}$ as defined in the following sections, and $l_{ri}$ denote the elements of $L$.

assumption$\sup_{i\geq 1}\mathrm{E}\left(r_i^2\right)=o(n^{-1})$.
assumption$\mathrm{E}\left\vert\epsilon_i\right\vert^q<C$, for some $q>4$.
assumption$\mathrm{E}\left(l_{ri}^2\right)<C$ and $\mathrm{E}\left(k_{ri}^2\right)<C$.

Assumption (ref) controls the approximation error. Sufficient conditions for various cases can be found in Chen2007. We avoid using specific rates because many types of approximation errors occur in the paper and explicitly allowing different rates of decay introduces extra notation without adding much insight beyond our general conditions. Assumption (ref) is fairly standard when establishing a central limit theorem for quadratic forms and is used to check a Lyapunov condition. Assumption (ref) imposes moment conditions on regressors and instruments.

To ease notation, we introduce the following definition that captures boundedness and non-multicollinearity of random or fixed matrices of growing dimension.

definitionFor any square, symmetric and positive semi-definite matrix $A$, let $\overline\alpha(A)$ and $\underline\alpha(A)$ denote its largest and smallest eigenvalues, respectively. If $A$ is random, we say that $A$ has Property G if \[ \overline{\alpha}(A)=O_p(1)\text{ and }\left\{\underline{\alpha}(A)\right\}^{-1}=O_p(1). \] If $A$ is non-random, we say that $A$ has Property G \protected@edef\@currentlabel{Property G} if \[ \limsup_{n\rightarrow\infty}\overline{\alpha}(A)<\infty\text{ and }\liminf_{n\rightarrow\infty}\underline{\alpha}(A)>0. \]
assumption$n^{-1}K_{1}'K_{1}$, $n^{-1}L'K_{1}K_{1}'L$ and $n^{-1}\Psi'\Psi$ have (ref).

We first show that the statistic $\mathbb{W}_1$ can be approximated by a quadratic form in $\epsilon$. Define $\mathcal{M}_1=\frac{1}{n}K_{1}\mathcal{V}_{1}K_{1}'$, with $\mathcal{V}_{1}=(K_{1}'K_{1})^{-1}K_{1}'LT_{1}^{-1}R'\mathscr{D}_{1}^{-1}RT_{1}^{-1}L'K_{1}(K_{1}'K_{1})^{-1}$. Then, we have the following theorem.

theoremUnder $H_0$, Assumptions (ref)-(ref) and \begin{flalign} \frac{1}{d_{\lambda}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+\frac{d_{\xi}}{n}\rightarrow 0, \end{flalign} as $n\rightarrow\infty$, we have \[ \mathbb{W}_1-\frac{\epsilon'\mathcal{M}_1\epsilon-d_\alpha}{\sqrt{2d_\alpha}}=o_p(1). \]

We remind the reader that results stated in this paper also hold when $d_\lambda$ and $d_\beta$ are held fixed and we do not state separate theorems for that case.

assumption$\left(I_n-\sum_{j=1}^{d_\lambda} \lambda_j W_j\right)^{-1}$ exists and is uniformly bounded in row and column sums for all sufficiently large $n$.

Define $H_{\ell 1}\equiv H_{\ell 1,n}: \alpha=\alpha^{*}\equiv\nu_{1n}d_{\alpha}^{\frac{1}{4}}/(n\nu_{1n}'\Gamma_{1n}\nu_{1n})^{\frac{1}{2}}$, with $\nu_{1n}$ a $d_{\alpha}\times 1$ non-zero vector and $\Gamma_{1n}$ a $d_{\alpha}\timesd_{\alpha}$ matrix defined in detail in the proofs. This sequence of local alternatives features a $d_\alpha^{1/4}$ damping factor that accounts for the cost of our nonparametric approach, and has been found in similar problems by Hong1995 and Gupta2018c, amongst others. We can now state the main theorem of the section.

theoremUnder Assumptions (ref)-(ref) and \begin{flalign} \frac{1}{d_{\lambda}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+\frac{d_{\xi}^3}{n}\rightarrow 0, \end{flalign} as $n\rightarrow\infty$, the following hold: (i) Under $H_0$, $\mathbb{W}_1\overset{d}{\rightarrow}N(0,1)$. (ii) $\mathbb{W}_1$ provides a consistent test. (iii) Under the sequence of local alternatives $H_{\ell 1}$, $\mathbb{W}_1\overset{d}{\rightarrow}N\left(2^{-1/2}, 1\right)$.

SAR with error spatial dependence

Having demonstrated our methods with a baseline model, we now consider ((ref)) but relax the $i.i.d$ assumption on $\epsilon_i$ to

equation[equation omitted — 102 chars of source]

where the $v_j$ are $i.i.d$ with mean 0 and variance 1, and $b_{ij}\equiv b_{ijn}$ can depend on $n$. This is the Kelejian2007 approach to error spatial dependence, and indeed we will employ the spatial heteroskedasticity and autocorrelation consistent (SHAC) method pioneered in that paper. This type of linear process assumption is widely employed in the literature, see e.g. Robinson2011, Robinson2012c, Hidalgo2017 and Conley2023.

Test statistic

Let $\mathscr{D}_2=\frac{1}{n}R\mathcal{U}_1R'$ with $\mathcal{U}_1= T_{1}^{-1}L'K_{1}(K_{1}'K_{1})^{-1}\Xi_1(K_{1}'K_{1})^{-1}K_{1}'LT_{1}^{-1}$ and $\widehat\mathscr{D}_2=\frac{1}{n}R\widehat\mathcal{U}_1R'$ with $\widehat\mathcal{U}_1=T_{1}^{-1}L'K_{1}(K_{1}'K_{1})^{-1}\hat\Xi_1(K_{1}'K_{1})^{-1}K_{1}'LT_{1}^{-1}$, where $\hat\Xi_1$ is the Kelejian2007 SHAC estimate of $\Xi_1=\frac{1}{n}K_{1}'\Sigma K_{1}$, with $\Sigma=E(\epsilon\epsilon')$ being the error variance matrix. We still use the 2SLS estimator $\hat\xi$ and the test statistic is now

equation[equation omitted — 131 chars of source]

We introduce the following assumptions:

assumption$\Sigma$ has (ref).
assumption$\sup_{i\geq 1}\sum_{j=1}^n\left\vert b_{ij}\right\vert+\sup_{j\geq 1}\sum_{i=1}^n\left\vert b_{ij}\right\vert<\infty$.
assumption$\mathscr{K}(\cdot): \mathbb{R}\rightarrow [-1,1]$ is a kernel function that satisfies $\mathscr{K}(0)=1$, $\mathscr{K}(x)=\mathscr{K}(-x)$, $\mathscr{K}(x)=0$ for $\vert x \vert>1$ and $\vert \mathscr{K}(x)-1 \vert\leq C\vert x \vert^\varrho, \vert x \vert\leq 1$, for some $\varrho\geq 1$.

Assumption (ref) restricts the spatial dependence in the errors to a manageable degree, see e.g. Kelejian2007 and Delgado2015 for similar assumptions. Assumption (ref) is a standard assumption on kernels in the HAC setting, see e.g. Kelejian2007.

Now we introduce distance measures $d_{ij,m}=d_{ji,m},m=1,\ldots,M$. As in Kelejian2007, we allow for measurement errors and so, in the following, let $d^*_{ij,m}=d^*_{ji,m}\geq 0$ be the actual distance measures used in practice. Corresponding to each measure, assume that the researcher can select a distance $d_m>0$ satisfying $d_m\uparrow \infty$ as $n\rightarrow\infty$. For each unit $i=1,\ldots,n$, let $\ell_{i}=\sum_{j=1}^n \left(1-\prod_{m=1}^M\mathbf{1}(d^*_{ij,m}>d_m)\right)$ and set $\ell=\max_{i}\ell_i$. Observe that $\ell_i$ is the number of units $j$ for which $d^*_{ij,m} \leq d_m$ for at least one $m=1,\ldots,M$. Also let $\sigma_{ij}$ be a typical element of $\Sigma$.

assumption(a) $\mathrm{E}(\ell^2)=o\left(n^{2\eta}\right)$ where $\eta<\frac{1}{2}(q-2) /(q-1)$ with $q>4$ in Assumption (ref). (b) $\sum_{j=1}^n \left\vert \sigma_{ij}\right\vert d^{\chi}_{ij,1}<C$ for some $\chi\geq 1$. (c) $d^*_{ij,m}=d_{ij,m}+\nu_{ij,m}\geq 0$, with $\left\vert \nu_{ij,m}\right\vert<C$ and $\nu_{ij,m}$ independent of $v_i$ for all $m=1,\ldots,M$.

Then, the $(r,s)$-th element of $\hat\Xi_1$ is

equation[equation omitted — 200 chars of source]

where $\hat{u}_i$ are elements of the estimation residual vector $\hat u=y-L\hat\xi$.

Asymptotic properties

In all subsequent lemmas/theorems, $q$ is defined in Assumption (ref) and $\eta$ is defined in Assumption (ref). In Lemma (ref) and Theorems (ref)-(ref), $d_{\xi}$ is defined as $d_{\xi}=d_{\lambda}+d_{\beta}+d_{\alpha}$. Furthermore, for any matrix $A$, let $\Vert A \Vert=\left\{\overline\alpha\left(A'A\right)\right\}^{\frac{1}{2}}$ i.e. the spectral norm of $A$.

lemmaLet Assumptions (ref), (ref), (ref), (ref)-(ref) hold, and \begin{flalign} \frac{1}{d_{\lambda}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+n^{\eta-1}d_{\xi}^2+n^{\eta(q-1)/q-(q-2)/2q}d_{\xi}\rightarrow 0, as n\rightarrow\infty. \end{flalign} Then, $\big\Vert\hat{\Xi}_1-\Xi_1\big\Vert=o_p(1)$.

Define $\mathcal{M}_2=\frac{1}{n}B'K_{1}\mathcal{V}_2K_{1}'B$, with $\mathcal{V}_2=(K_{1}'K_{1})^{-1}K_{1}'LT_{1}^{-1}R'\mathscr{D}_{2}^{-1}RT_{1}^{-1}L'K_{1}(K_{1}'K_{1})^{-1}$. We have the following theorem.

theoremLet ((ref)) hold. Then, under $H_0$, Assumptions (ref), (ref), (ref), (ref)-(ref), as $n\rightarrow\infty$, \[ \mathbb{W}_2-\frac{v'\mathcal{M}_2v-d_\alpha}{\sqrt{2d_\alpha}}=o_p(1). \]

To derive the asymptotic distribution of $\mathbb{W}_2$, we need the following assumption, as also in Delgado2015.

assumption$\mathrm{E}\left\vert v_i\right\vert^s<C$, for some $s\geq 8$.

Define $H_{\ell 2}\equiv H_{\ell 2,n}: \alpha=\alpha^{*}\equiv\nu_{2n}d_{\alpha}^{\frac{1}{4}}/(n\nu_{2n}'\Gamma_{2n}\nu_{2n})^{\frac{1}{2}}$, with $\nu_{2n}$ a $d_{\alpha}\times 1$ non-zero vector and $\Gamma_{2n}$ a $d_{\alpha}\timesd_{\alpha}$ matrix defined in the proof of the next theorem.

theoremUnder Assumptions (ref), (ref)-(ref), with \begin{flalign} \frac{1}{d_{\lambda}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+\frac{d_{\xi}^3}{n}+n^{\eta-1}d_{\xi}^2+ \ n^{\eta(q-1)/q-(q-2)/2q}d_{\xi}\rightarrow 0, \end{flalign} as $n\rightarrow\infty$, the following hold: (i) Under $H_0$, $\mathbb{W}_2\overset{d}{\rightarrow}N(0,1)$. (ii) $\mathbb{W}_2$ provides a consistent test. (iii) Under the sequence of local alternatives $H_{\ell 2}$, $\mathbb{W}_2\overset{d}{\rightarrow}N(2^{-1/2}, 1)$.

Observe that the condition $n^{\eta-1}d_{\xi}^2\rightarrow 0$ implies ${d_{\xi}^3}/{n}\rightarrow 0$ if $\eta\geq 1/3$, so depending on the value of $q$, the rate condition ((ref)) can be simplified.

Misspecification robustness

Consider a degree of misspecification robustness in ((ref)), where $\epsilon_i$ has the structure in Section (ref). Focusing on the case $d_\lambda=1$ for notational simplicity, and writing $W_1=W$, ((ref)) imposed the parametric form $\lambda Wy$ and took $W$ as known in this “spatial lag” term. We now allow this to take the nonparametric form $Gy$, where $G$ has elements $g\left(d_{rs}\right)$, $r,s=1,\ldots,n$, for some unknown function $g(\cdot)$, and a vector of exogenous (independent of $v_j$, $j=1,\ldots,n$) economic distance measures $d_{rs}$. This idea was introduced by pinkse2002, and has been much used since (see e.g. Sun2016 and Gupta2024).

Test statistic

Following pinkse2002, we approximate $g\left(d_{ij}\right)$ of $G$ with a series of basis functions

equation[equation omitted — 80 chars of source]

where $\tau_l$ are unknown coefficients, and $e_{l}$ form a basis of the function space to which $g(\cdot)$ belongs. Let $c_i'$ be the $i$-th row of the $n \times d_{\tau}$ matrix $\mathfrak{C}$ with typical $(i,l)$-th element $\sum_{j\neq i}e_l\left(d_{ij}\right)y_j$ and $\tau=\left(\tau_1,\ldots,\tau_{d_{\tau}}\right)'$, $d_{\tau}\rightarrow\infty$ as $n\rightarrow\infty$. Then ((ref)) is extended to

equation[equation omitted — 127 chars of source]

and ((ref)) to

equation[equation omitted — 105 chars of source]

with $u_i=r_{iG}+r_{i\delta}+\epsilon_i$, where $r_{iG}=\sum_{l=d_{\tau}+1}^{\infty}\tau_l\sum_{j\neq i}e_l\left(d_{ij}\right)y_j$ and $r_{i\delta}=p_{i}'\delta\left(z_{i}\right)-\psi_{i}'\alpha$. In the matrix notation, we now have

equation[equation omitted — 101 chars of source]

where $F=[\mathfrak{C}, X,\Psi]$ and $\theta=\left(\tau',\beta',\alpha'\right)'$. Our test is based on the 2SLS estimator

equation[equation omitted — 95 chars of source]

where $P_{K_{2}}=K_{2}\left(K_{2}'K_{2}\right)^{-1}K_{2}'$, $K_{2}$ being an $n\times J_2$ instrument matrix with $J_2\geq d_\theta=d_{\tau}+d_{\beta}+d_{\alpha}$, but with $J_2$ and $d_\theta$ having the same asymptotic order.

Let $\mathscr{D}_3=\frac{1}{n}R\mathcal{U}_2R'$ with $\mathcal{U}_2=T_{2}^{-1}F'K_{2}(K_{2}'K_{2})^{-1}\Xi_2(K_{2}'K_{2})^{-1}K_{2}'FT_{2}^{-1}$ and $T_2=\frac{1}{n} F'{P_{K_{2}}}F$, and let $\widehat\mathscr{D}_3=\frac{1}{n}R\widehat\mathcal{U}_2R'$ with $\widehat\mathcal{U}_2= T_{2}^{-1}F'K_{2}(K_{2}'K_{2})^{-1}\hat\Xi_2(K_{2}'K_{2})^{-1}K_{2}'FT_{2}^{-1}$, where $R=\left[0_{d_\alpha\times\left(d_\tau+d_\beta\right)}, I_{d_\alpha}\right]$ and $\hat\Xi_2$ is the Kelejian2007 estimate of $\Xi_2=\frac{1}{n}K_{2}'\Sigma K_{2}$, with $\Sigma=E(\epsilon\epsilon')$. Then, the test statistic is

equation[equation omitted — 135 chars of source]

Denote the elements of $F$ as $f_{ri}$. We have following assumptions.

assumption$\mathrm{E}\left(f^2_{ri}\right)<C$ and $\mathrm{E}\left(k_{ri}^2\right)<C$.
assumption$\sup_{i\geq 1}\mathrm{E}\left(r_{iG}^2\right)=o(n^{-1})$ and $\sup_{i\geq 1}\mathrm{E}\left(r_{i\delta}^2\right)=o(n^{-1})$.
assumption$n^{-1}K_{2}'K_{2}$ and $ n^{-1}F'K_{2}K_{2}'F$ have (ref).
lemmaLet Assumptions (ref)-(ref), (ref)-(ref) hold, and \begin{flalign} \frac{1}{d_{\tau}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+n^{\eta-1}d_\theta^2+\ n^{\eta(q-1)/q-(q-2)/2q}d_\theta\rightarrow 0, as n\rightarrow\infty. \end{flalign} Then, $\big\Vert\hat{\Xi}_2-\Xi_2\big\Vert=o_p(1)$.

Asymptotic properties

Define $\mathcal{M}_3=\frac{1}{n}B'K_{2}\mathcal{V}_3K_{2}'B$, with $\mathcal{V}_3=(K_{2}'K_{2})^{-1}K_{2}'FT_{2}^{-1}R'\mathscr{D}_{3}^{-1}RT_{2}^{-1}F'K_{2}(K_{2}'K_{2})^{-1}$. Then, we have the following theorem.

theoremLet ((ref)) hold. Then, under $H_0$, Assumptions (ref)-(ref), (ref)-(ref), as $n\rightarrow\infty$, \[ \mathbb{W}_3-\frac{v'\mathcal{M}_3v-d_\alpha}{\sqrt{2d_\alpha}}=o_p(1). \]

Define $H_{\ell 3}\equiv H_{\ell 3,n}: \alpha=\alpha^{*}\equiv\nu_{3n}d_{\alpha}^{\frac{1}{4}}/(n\nu_{3n}'\Gamma_{3n}\nu_{3n})^{\frac{1}{2}}$, with $\nu_{3n}$ a $d_{\alpha}\times 1$ non-zero vector and $\Gamma_{3n}$ a $d_{\alpha}\timesd_{\alpha}$ matrix defined in the proofs. The following assumption generalizes Assumption (ref).

assumption$\left(I_n-G\right)^{-1}$ exists and is uniformly bounded in row and column sums for all sufficiently large $n$, almost surely and uniformly over the support of the distance measure vector.
theoremUnder Assumptions (ref)-(ref), with \begin{flalign} \frac{1}{d_{\tau}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+\frac{d_\theta^3}{n}+n^{\eta-1}d_\theta^2+n^{\eta(q-1)/q-(q-2)/2q}d_\theta\rightarrow 0, \end{flalign} as $n\rightarrow\infty$, the following hold: (1) Under $H_0$, $\mathbb{W}_3\overset{d}{\rightarrow}N(0,1)$. (2) $\mathbb{W}_3$ provides a consistent test. (3) Under the sequence of local alternatives $H_{\ell 3}$, $\mathbb{W}_3\overset{d}{\rightarrow}N\left(2^{-1/2}, 1\right)$.

Varying spatial coefficients

Now we revert to the case of known $W_j, j=1,\ldots,d_\lambda$, but allow spatial lag parameters $\lambda_j$ to vary. The model is closely related to malikov2017semiparametric and Sun2018, with the difference that we do not study panel data settings. Instead, we allow $d_\beta\rightarrow\infty$ and general spatial dependence in the errors. Our model also relates to Sun2024semiparametric, wherein the model includes a nonparametric nuisance regression function but does not feature dependent errors. Thus, our model complements the extant literature.

Test statistic

Specifically, we have

equation[equation omitted — 169 chars of source]

with $\lambda_j(\cdot)$ unknown scalar-valued functions and $\epsilon_i=\sum_{j=1}^n b_{ij}v_j$, i.e. error with spatial dependence. Introduce the approximation $\lambda_j\left(z_i\right)=\sum_{k=1}^{l_j}\mu_{kj}\phi_{kj}\left(z_i\right)+r_{ij\lambda},j=1,\ldots,d_\lambda$. Set $d_{\mu}=\sum_{j=1}^{d_\lambda} l_j$, $\phi_j=\left(\phi_{1j},\ldots,\phi_{l_jj}\right)'$, and let $H(z)$ be the $n\times d_{\mu}$ matrix with $i$-th row $\left(y'w_{i,1}\phi_1\left(z_i\right),\ldots,y'w_{i,d_\lambda}\phi_{d_\lambda}\left(z_i\right)\right)$. Then,

equation[equation omitted — 65 chars of source]

where $G(z)=\left[H(z),X,\Psi\right]$, $\gamma=\left(\mu',\beta',\alpha'\right)'$, $\mu=\left(\mu_1',\ldots,\mu_{d_{\mu}}'\right)$, and $u$ has $i$-th element $u_i=r_{i\lambda}+r_{i\delta}+\epsilon_i$ with $r_{i\lambda}=\sum_{j=1}^{d_\lambda}r_{ij\lambda}w_{i,j}'y$. Defining $\Lambda_j(z)=diag\left(\lambda_j\left(z_1\right),\ldots,\lambda_j\left(z_n\right)\right)$ and $S(z)=I-\sum_{j=1}^{d_\lambda}\Lambda_j(z)W_j$, under the usual invertibility conditions,

equation[equation omitted — 78 chars of source]

with $P$ having $i$-th row $p_i'$ and $t(z)$ serving to stack $p_i'\delta\left(z_i\right)$ into an $n\times 1$ vector with conformable subscript to $y_i$ and $\epsilon_i$.

The null hypothesis in this setting is

equation[equation omitted — 48 chars of source]

where $\vartheta$ is either $\mu$ or $\alpha$ and our test will again be implemented using the 2SLS estimator

equation[equation omitted — 90 chars of source]

where $P_{K_{3}}=K_{3}\left(K_{3}'K_{3}\right)^{-1}K_{3}'$ with $K_{3}$ being an $n\times J_3$ instrument matrix with $J_3\geqd_{\gamma}=d_{\mu}+d_{\beta}+d_{\alpha}$, but with $J_3$ and $d_{\gamma}$ having the same asymptotic order. Note that we have omitted the $z$ argument for brevity. Let $\mathscr{D}_4=\frac{1}{n}R\mathcal{U}_3R'$ with \sloppy $\mathcal{U}_3=T_{3}^{-1}G'K_{3}(K_{3}'K_{3})^{-1}\Xi_3(K_{3}'K_{3})^{-1}K_{3}'LT_{3}^{-1}$ and $T_3=\frac{1}{n} G'{P_{K_{3}}}G$, and let $\widehat\mathscr{D}_4=\frac{1}{n}R\widehat\mathcal{U}_3R'$ with $\widehat\mathcal{U}_3=T_{3}^{-1}G'K_{3}(K_{3}'K_{3})^{-1}\hat\Xi_3(K_{3}'K_{3})^{-1}K_{3}'LT_{3}^{-1}$, where $R=\left[I_{d_{\mu}},0_{d_{\mu}\times\left(d_\beta+d_\alpha\right)} \right]$ or $R=\left[0_{d_{\alpha}\times\left(d_{\mu}+d_{\beta}\right)}, I_{d_{\alpha}}\right]$, and $\hat\Xi_3$ is the Kelejian2007 estimate of $\Xi_3=\frac{1}{n}K_{3}'\Sigma K_{3}$, with $\Sigma=E(\epsilon\epsilon')$. Then, the test statistic is

equation[equation omitted — 143 chars of source]

where $d_{\vartheta}$ is either $d_{\mu}$ or $d_{\alpha}$. \fussy Denote the elements of $G$ as $g_{ri}(z)$. We have the following assumptions:

assumption$\sup_{z\in\mathcal Z} \mathrm{E}\left(g^2_{ri}(z)\right)<C$ and $\mathrm{E}\left(k_{ri}^2\right)<C$.
assumption$\sup_{i\geq 1}\mathrm{E}\left(r_{i\lambda}^2\right)=o(n^{-1})$ and $\sup_{i\geq 1}\mathrm{E}\left(r_{i\delta}^2\right)=o(n^{-1})$.
assumption$n^{-1}K_{3}'K_{3}$ and $ n^{-1}G(z)'K_{3}K_{3}'G(z)$ have (ref), the latter uniformly in $z\in \mathcal{Z}$.
lemmaLet Assumptions (ref)-(ref), and (ref) -(ref) hold, and \begin{flalign} \frac{1}{d_{\mu}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+ n^{\eta-1}d_{\gamma}^2+ n^{\eta(q-1)/q-(q-2)/2q}d_{\gamma}\rightarrow 0, as n \rightarrow \infty. \end{flalign} Then, we have $\big\Vert\hat{\Xi}_3-\Xi_3\big\Vert=o_p(1)$.

Asymptotic properties

Define $\mathcal{M}_{4}=\frac{1}{n}B'K_{3}\mathcal{V}_{4}K_{3}'B$, with $\mathcal{V}_{4}=(K_{3}'K_{3})^{-1}K_{3}'GT_{3}^{-1}R'\mathscr{D}_{4}^{-1}RT_{3}^{-1}G'K_{3}(K_{3}'K_{3})^{-1}$. Then, we have the following theorem.

theoremLet ((ref)) hold. Then, under $H_0$, Assumptions (ref)-(ref), (ref)-(ref), as $n\rightarrow\infty$, \[ \mathbb{W}_4-\frac{v'\mathcal{M}_4v-d_{\vartheta}}{\sqrt{2d_{\vartheta}}}=o_p(1). \]

The following assumption is the appropriate version of Assumptions (ref) and (ref) for the model considered in this section.

assumption$\left(I_n-S(z)\right)^{-1}$ exists and is uniformly bounded in row and column sums for all sufficiently large $n$, almost surely, uniformly in $z$.

Define $H_{\ell 4}\equiv H_{\ell 4,n}: \vartheta=\vartheta^{*}\equiv\nu_{4n}d_{\vartheta}^{\frac{1}{4}}/(n\nu_{4n}'\Gamma_{4n}\nu_{4n})^{\frac{1}{2}}$, with $\nu_{4n}$ a $d_{\vartheta}\times 1$ non-zero vector and $\Gamma_{4n}$ a $d_{\vartheta}\timesd_{\vartheta}$ matrix defined in the proof of the next theorem.

theorem\sloppy Under Assumptions (ref)-(ref) and (ref)-(ref), suppose that \begin{flalign} \frac{1}{d_{\mu}}+\frac{1}{d_{\beta}}+\frac{1}{d_{\alpha}}+\frac{d_{\gamma}^3}{n}+ n^{\eta-1}d_{\gamma}^2+ n^{\eta(q-1)/q-(q-2)/2q}d_{\gamma}\rightarrow 0, \end{flalign} as $n\rightarrow\infty$, the following hold: (i) Under $H_0$, $\mathbb{W}_4\overset{d}{\rightarrow}N(0,1)$. (ii) $\mathbb{W}_4$ provides a consistent test. (iii) Under the sequence of local alternatives $H_{\ell 4}$, $\mathbb{W}_4\overset{d}{\rightarrow}N\left(2^{-1/2},1\right)$.

\allowdisplaybreaks

Empirical applications

Testing for the CRS in a production function

In this subsection, we are interested in testing whether the production function in China's nonmetal mineral manufacturing industry has a constant returns to scale (CRS) technology. CRS refers to a scenario where a simultaneous proportional increase in all inputs results in an identical proportional increase in output. Returns to scale is a classical concept in economics to assess the efficiency of a production function, and it has accordingly received much attention (see ackerberg2015identification,basu2017uncertainty,attanasio2020estimating,combes2021production ). Existing works, however, mostly do not account for cross-sectional dependence between firms. The role of cross-firm correlation has recently been highlighted empirically by iyoha2023estimating, who shows that this can capture spillover effects and would lead to model misspecification if left unaccounted for.

Our models generalize li2002semiparametric's production function by incorporating spatial dependence. Their analysis centers on a Cobb-Douglas production function that allows elasticities of capital and labor to be nonparametrically varying with management expenses. Management expenses are costs that are indirectly associated with output production, such as research and development (R&D), equipment upgrades and employee training. CRS corresponds to the elasticities of capital and labor summing to one.

Using data from the Third Industrial Census of China (conducted by the National Statistical Bureau in 1995), li2002semiparametric find that ignoring the varying nature of these coefficients leads to an underestimation of the returns to scale. They suggest their estimate of returns to scale of the technology is close to being constant for most values of managerial expense, but they do not have a test for the CRS hypothesis. Our data comes from the Chinese Industrial Enterprise Database in 2014 as it is the most recent information available. We use all firms in the nonmetal mineral manufacturing industry with the survey code 31, resulting in a sample of 7,355 observations. As advocated by li2002semiparametric, production functions of firms in this industry are expected to be homogeneous since it has a very small proportion of foreign entity ownership. This is important because evidence suggests that production performance may vary across ownership types (murakami1994technical).

Variables are in logs: output ($y_{i}$), capital ($p_{i1}$), labor ($p_{i2}$), and management expenses ($z_{i}$). As in li2002semiparametric, we include a parametric benchmark model:

flaligny_{i}=\delta_{0}+\delta_{1}p_{i1}+\delta_{2}p_{i2}+\beta_{1}z_{i}+\beta_{2}z_{i}^2+\epsilon_{i}.

Elasticities of capital and labor are $\delta _{1}$ and $\delta _{2}$, respectively. The semiparametric model where the elasticities of of capital and labor vary with management expenses is: \setlength\abovedisplayskip{3pt} \setlength\belowdisplayskip{3pt}

flalign*y_{i}=\delta_{0}(z_i)+\delta_{1}(z_i)p_{i1}+\delta_{2}(z_i)p_{i2}+\epsilon_{i}.

Based on the CRS hypothesis, $H_{0}^{true}:\delta _{1}(z_{i})+\delta _{2}(z_{i})=1$, we reparameterize the model as

flalign*y_{i}^{*}=\delta_{0}(z_i)+\delta_{2}(z_i)(p_{i2}-p_{i1})+(\delta_{1}(z_i)+\delta_{2}(z_i)-1)p_{i1}+\epsilon_{i},

where $y_{i}^{\ast }=y_{i}-p_{i1}$. Then, we estimate the following model with no SAR structure:

flalign*\begin{split} & y_{i}^{*}=\beta_0+\sum\limits_{k=1}^{h}\alpha_{1k}\psi_{ik}+\sum\limits_{k=1}^{h}\alpha_{2k}(p_{i2}-p_{i1})\psi_{ik}+\sum\limits_{k=1}^{h}\alpha_{3k}p_{i1}\psi_{ik}+\epsilon_{i}, \end{split}

where $\psi _{ik}$ is a basis function, with argument $z_i$, for $k=1,...,h/2$. We approximate the unknown functions by polynomial and trigonometric functions of different orders. We define the series approximation based Wald statistic derived from this model as $\mathbb{W}_{0}$. Next, we estimate a number of models and perform tests based on the various configurations and statistics discussed in previous sections: \setlength\abovedisplayskip{5pt} \setlength\belowdisplayskip{5pt}

flaligny_{i}^{*} = \beta_{0} + \lambda\sum\limits_{j=1}^{n} w_{ij} y_j + \sum\limits_{k=1}^{h}\alpha_{1k}\psi_{ik} + \sum\limits_{k=1}^{h}\alpha_{2k}(p_{i2}-p_{i1})\psi_{ik} + \sum\limits_{k=1}^{h}\alpha_{3k}p_{i1}\psi_{ik} + \epsilon_{i}.

We compute the Wald statistics as follows: $\mathbb{W}_{1}$ from (ref) with $i.i.d$ errors $\epsilon_{i}$; $\mathbb{W}_{2}$ from (ref) with $\epsilon_{i}=\sum_{j=1}^{n}b_{ij}v_{j};$ $\mathbb{W}_{3}$ from (ref) replacing $\lambda\sum_{j=1}^{n}w_{ij}y_j$ by $\sum_{l=1}^{d_\tau}\tau_{l} c_{il}$ and setting $\epsilon_{i}=\sum_{j=1}^{n} b_{ij} v_{j};$ $\mathbb{W}_{4}$ from (ref) with $\lambda\sum_{j=1}^{n}w_{ij}y_j$ replaced by $\sum_{m=1}^{d_\mu}\mu_m \mathfrak{H}_{im}$ and $\varepsilon_i=\sum_{j=1}^{n} b_{ij}v_j$, where $\mathfrak{H}_{im}$ is the $(i,m)$-th entry of the $n\times d_{\mu}$ matrix $H(z)$ defined in Section (ref). The coefficient on $p_{i1}$ remains the measure of returns to scale technology even with a SAR feature, although this notion differs from the overall returns to scale that combines technological and other spillovers effects. Let $d_{ij,m}^{\ast }$ be the geographical distance between firms $i$ and $j$, and $d_{m}$ represent the 10th percentile of all geographical distances. We use two different row-normalized spatial weight matrices, denoted

flalign*\begin{split} w_{ij}^{p}=[W_{n}^{p}]_{ij}=\begin{cases} 1 \quad if \ i \ and \ j \ in same city,\\ 0 \quad otherwise. \end{cases} w_{ij}^{d}=[W_{n}^{d}]_{ij}=\begin{cases} \frac{1}{d_{ij,m}^{*}} \quad if \ d_{ij,m}^{*}<d_m,\\ 0 \quad otherwise. \end{cases} \end{split}

For the test statistics, $\mathbb{W}_{i}$, $i=0,...,3$, $\alpha _{3}=(\alpha _{31},...\alpha _{3h})^{\prime }$, and the null is $H_{0}:\alpha _{3}=0$. For $\mathbb{W} _{4} $, the null is $H_{0}:\alpha _{3}=0$ or $H_{0}:\mu =0$, where $\mu =(\mu_{1},...,\mu_{d_\mu})^{\prime}$. Throughout our empirical analysis, we use $\psi _{ik}(z)=z_{i}^{k}$ as the polynomial function for $k=1,...,h$, and $\psi _{ik}(z)=[\sin (kz_{i}),\cos (kz_{i})]$ as the trigonometric function $k=1,...,h/2$ for $h=2,4$. To compute $\mathbb{W}_{2}$, $\mathbb{W}_{3}$ and $\mathbb{W}_{4}$, we use the Epanechnikov kernel function. For $\mathbb{W}_{3}$, let $\mathcal{E}_{l}$ be the matrix whose typical entry is, for $i\neq j$ and $l=1,...,d_\tau$, $e_{l}(d_{ij})=[d_{ij,m}^{\ast }]^{l}\mathbf{1}(d_{ij,m}^{\ast }<d_{m})$. Then, $c_{il}=\sum_{j=1}^{n}e_{l}(d_{ij})y_{j}$. For $\mathbb{W}_{4}$, we use $\phi_{im}(z)=\frac{1}{d_\mu}[\frac{2}{\pi}\tanh(z_{i})]^{m}$ as the polynomial function and $ \phi_{im}(z)=\frac{1}{d_{\mu}}\sin(\frac{z_{i}}{2m})$ as the trigonometric function for $ m=1,...,d_\mu$. Then, $\mathfrak{H}_{im}=(\sum_{j=1}^{n}w_{ij}y_{j})\phi_{im}(z)$, where $w_{ij}$ denotes the $(i,j)$-th entry of either $W_n^{p}$ or $W_n^{d}$. For the IVs, set $K_{1}=K_{3}=[X,\Psi,W_n^{p}X]$ when the specification uses $W_n^{p}$, and $K_{1}=K_{3}=[X,\Psi,W_n^{d}X]$ when it uses $W_n^{d}$, for $\mathbb{W}_{1}$, $\mathbb{W}_{2}$, and $\mathbb{W}_{4}$, and $K_{2}=[X, \Psi, \mathcal{E}_{1}X,..., \mathcal{E}_{d_\tau}X]$ with respect to $\mathbb{W}_{3}$. We set $d_\mu=d_\tau=h$, and the results with different basis functions are in Table (ref). In all models, the hypothesis of CRS cannot be rejected. These findings corroborate those in li2002semiparametric. These non-rejection results are economically meaningful. One can conclude that the CRS technology is a salient feature of the production function in this industry, as it is present in the parametric model as well as in semiparametric models with varying coefficients and cross-firm dependence.

We end this application with some remarks: (1) In addition to inference on CRS, our analysis also reveals other similarities with prior works. Table (ref) shows that the parametric method underestimates the returns to scale compared to the semiparametric method.\footnote{ The returns to scale in a semiparametric model are computed by summing $\frac{1}{n}\sum_{i=1}^{n} \hat{\delta}_{1}(z_{i})$ and $\frac{1}{n}\sum_{i=1}^{n}\hat{\delta}_{2}(z_{i})$.} The returns to scale with respect to different input factors also differ, which highlights the crucial role of management expenses. Figures (ref) (a)-(b) in the online appendix illustrate this with the semiparametric estimates with spatial weights $W_{n}^{d}$. They show that the output elasticity of capital is increasing in management expenses, while it is decreasing for labor. The decreasing elasticity of labor is a known empirical fact in this industry due to \textquotedblleft concealed unemployment\textquotedblright\, where the government did not allow state-owned firms to lay off extra employees as a strategy to avoid social unrest. (2) The returns to scale are relatively flat in the semiparametric model and are within the 95% bounds for most $ z_{i}$, as illustrated in Figure (ref) (c) in the online appendix. Further discussions on these findings can be found in li2002semiparametric. (3) Our results are not unique to the estimates from $W_{n}^{d}$. Other choices of spatial weights give the same conclusion.

table[table omitted — 2,221 chars of source]

Testing the impact of the distance on house prices

We use the Boston house price dataset \footnote{ The dataset consists the median value of owner-occupied homes in 506 census tracts in the Boston Standard Metropolitan Statistical Area in 1970 (harrison1978hedonic), which is available at \url{https://github.com/simonbrewer/geog6000/blob/main/datafiles/boston.tr.zip} .} to explore how some factors affect the median value of owner-occupied homes in $1000$s (MEDV), and whether the effects of these factors vary over location. Factors include the per capita crime rate by town (denoted by CRIM), average number of rooms per dwelling (RM), full-value property-tax rate per \$10,000 dollar (TAX), percentage of lower socioeconomic status population (LSTAT) and index of accessibility to radial highways (RAD). The location variable DIS, which represents the weighted distances between a property and the five Boston employment centers, will drive the varying coefficient. Let $y_{i},x_{i1},x_{i2},p_{i1},p_{i2},p_{i3},z_{i}$ denote the logarithm of MEDV, RAD, LSTAT, CRIM, RM, TAX and DIS respectively.

As house prices in a specific location or region are commonly influenced by the prices of nearby houses due to factors such as environmental and geographic considerations, it is natural to consider the spatial autocorrelation.\footnote{ There are also some studies that use the SAR semiparametric model with the Boston house price data, for example, malikov2017semiparametric and li2019tests. However, the former considers the covariate $z_{i}$ as NO$_{2}$, which is more about the impact of environmental pollution on housing prices. The latter concentrates on the geographical weighted regression (GWR) method.} Sun2014 proposed a semiparametric coefficient varying spatial model and derived the distribution theory of their estimator. While they do not have inference results on functions of the random coefficients, they have a model selection procedure to determine the parametric and nonparametric components. We thus take the model they selected as our starting point.

Guided by model selection results, Sun2014 argue that DIS doesn't affect LSTAT and RAD. Intuitively, for LSTAT, the geographical distribution of employment opportunities does not necessarily correlate with a community's socio-economic status; a community might be far from employment centers but could have good public transportation links or a high number of remote work opportunities, mitigating the impact on the community's lower socioeconomic status. For RAD, the accessibility to highways might depend more on urban planning and infrastructure investment rather than the proximity to employment centers. The baseline model, based on Sun2014, is:

flalign*y_i=\lambda\sum\limits_{j=1}^{n} w_{ij}y_j+\sum\limits_{k=0}^{2}\beta_{k}x_{ik}+\sum\limits_{m=1}^{3}\delta_{m}(z_{i})p_{im}+\epsilon_{i}.

We consider a more general version where the spatial lag coefficient varies with DIS:

flalign*y_i=\lambda(z_{i})\sum\limits_{j=1}^{n} w_{ij}y_j+\sum\limits_{k=0}^{2}\beta_{k}x_{ik}+\sum\limits_{m=1}^{3}\delta_{m}(z_{i})p_{im}+\epsilon_{i}.

Let $d_\tau=d_\mu=h$, we estimate:

flalign*\begin{split} &y_i=\lambda\sum\limits_{j=1}^{n} w_{ij}y_j+\sum\limits_{k=0}^{2}\beta_{k}x_{ik}+\sum\limits_{m=1}^{3}\sum\limits_{k=1}^{h}\alpha_{mk}p_{im}\psi_{ik}+\epsilon_{i}, \ and using \ \mathbb{W}_{1}; \\ &y_i=\lambda\sum\limits_{j=1}^{n} w_{ij}y_j+\sum\limits_{k=0}^{2}\beta_{k}x_{ik}+\sum\limits_{m=1}^{3}\sum\limits_{k=1}^{h}\alpha_{mk}p_{im}\psi_{ik}+\epsilon_{i}, \ \ \epsilon_{i}=\sum\limits_{j=1}^{n}b_{ij}v_{j} , \ using \ \mathbb{W}_{2}; \\ &y_i=\sum\limits_{l=1}^{h}\tau_{l}c_{il}+\sum\limits_{k=0}^{2}\beta_{k}x_{ik}+\sum\limits_{m=1}^{3}\sum\limits_{k=1}^{h}\alpha_{mk}p_{im}\psi_{ik}+\epsilon_{i}, \ \ \epsilon_{i}=\sum\limits_{j=1}^{n}b_{ij}v_{j} , \ using \ \mathbb{W}_{3}; \\ &y_i=\sum\limits_{l=1}^{h}\mu_{l}\mathfrak{H}_{il}+\sum\limits_{k=0}^{2}\beta_{k}x_{ik}+\sum\limits_{m=1}^{3}\sum\limits_{k=1}^{h}\alpha_{mk}p_{im}\psi_{ik}+\epsilon_{i}, \ \ \epsilon_{i}=\sum\limits_{j=1}^{n}b_{ij}v_{j} , \ using \ \mathbb{W}_{4}. \end{split}

We consider two row-normalized spatial weight matrices. The first is the first-order queen matrix, $W_{n}^{q}$, discussed in malikov2017semiparametric. The second, $W_{n}^{t}$, indicates whether two locations are in the same tract. For $\mathbb{W}_{i}$, $i=0,...,3$, $\alpha=(\alpha_{11},...\alpha_{1h},...,\alpha_{3h})^{\prime }$, and the null is $H_{0}:\alpha =0$. For $\mathbb{W}_{4}$, the null is $ H_{0}:\alpha =0$ or $H_{0}:\mu=0$, where $\mu=(\mu_{1},...,\mu_{h})^{\prime}$. Results are listed in Table (ref). All tests reject the null $H_{0}:\alpha =0$. In addition, testing with $ \mathbb{W}_{4}$ leads to a rejection of $H_{0}:\mu =0$, indicating a significant nonlinear relationship between $log(DIS)$ and neighboring house prices. We further estimate the varying $\lambda (z_{i})$ and evaluate its empirical mean $\hat{\lambda}=\frac{1}{n}\sum_{i=1}^{n}\lambda (z_{i})$. As shown in Panel E of Table (ref), the parametric method underestimates the spatial coefficient. Thus, our results based on more general models complement those in Sun2014.

Monte Carlo Simulation

We present a basic set of Monte Carlo results in this section and a wider range of simulations in the online appendix. All experiments use 1000 replications.

Basic setting

Taking $n = 200, 500, 900,$ we choose two specifications to generate $y$ with $d_{\delta}=1$:

flalign*y_{i}=\sum_{k=1}^{d_{\lambda}}\lambda_{k}\sum\limits_{j=1}^{n}w_{kij} y_j+x_{i}^{\prime}\beta+p_{i}\delta\left(z_{i}\right)+\epsilon_{i}, i=1, \ldots, n,

and

flalign*y_{i}=\sum_{k=1}^{d_{\lambda}}\lambda_{k}\sum\limits_{j=1}^{n}w_{kij} y_j+x_{i}^{\prime}\beta+p_{i}\delta\left(z_{i}\right)+v_{i}, \ v_i=\sum\limits_{j=1}^{n}b_{ij}v_j+\epsilon_{i}, i=1, \ldots, n,

where $\epsilon_i$ are $i.i.d$ with three zero mean and unit variance distributions: (V1) \protected@edef\@currentlabel{(V1)} $N(0,1)$, (V2) \protected@edef\@currentlabel{(V2)} $\sqrt{\frac{5}{4}}t(10)$ and (V3) \protected@edef\@currentlabel{(V3)} $\frac{1}{4}(\chi_{8}^{2}-8)$. We generate $z_i\stackrel{i.i.d}{\sim}U[0,1]$, $p_i\stackrel{i.i.d}{\sim}U[-2,2]$, $x_{i1}=1$, $x_{i2}\stackrel{i.i.d}{\sim}N(1, 2)$ and set $\beta=(-1,1)'$. For $\lambda$, we fix the total weight at 0.9 and then define a decreasing vector $\lambda^{*}=(d_{\lambda},d_{\lambda}-1,...,1)'$. We then set $\lambda=0.9\left(\sum_{k=1}^{d_{\lambda}}\lambda_{k}^{*}\right)^{-1}\lambda^{*}$. For the spatial weights matrices, we generate $W_{k}$ as circulants. Specifically, let $W_{k}^*$ be the symmetric circulant matrix whose first-row entries are $ w_{k1j}^*= 0$ if $j=1$ or $j=k+2, \ldots, n-k$ and $w_{k1j}^*= 1$ if $j=2, \ldots, k+1$ or $j=n-k+1, \ldots, n$. Thus, the weight matrix $W_{k}^*$ encapsulates a binary neighbourhood criterion for $k$ neighbours on either `side' of a unit. Now define the normalized matrix $W_{k}=\{\overline{\alpha}(W_{k}^*)\}^{-1} W_{k}^*,$ recalling that $\overline{\alpha}(W_{k}^*)=2k$ for a circulant matrix. Furthermore, $b_{ij}$ is the typical element of $(I_n-\sum_{k=1}^{d_{\lambda}}\lambda_{k}W_{k})^{-1}$. To construct the SHAC statistics $\mathbb{W}_{2}, \mathbb{W}_{3}$ and $\mathbb{W}_{4}$, we proceed as follows. Generate $M=d_{\lambda}$ distinct distance measures by setting $\ell=[n^{\eta}]+1$ with $\eta=3/7$ and taking $q=8$, so that each unit has at least $\ell$ neighbors, where $\eta$ is defined in Assumption (ref) (a).\footnote{From the symmetric circulant matrix $W_{d_{\lambda}}$, we build a graph and compute its pairwise shortest‐path distances in MATLAB, then apply cmdscale to generate $n\times 2$ coordinate matrix. We next compute the full $n\times n$ Euclidean distance matrix $D$ with the typical elements $d_{ij}$.} For each $m=1, ..., M$, we add the actual distance matrix, $D_{m}^{*}$, used in practice with unobserved symmetric measurement noise. Its typical entries are $ d_{ij,m}^{*}=d_{ij}+\nu_{ij,m},$ where $\nu_{ij,m}=1/2(\mu_{ij,m}+\mu_{ji,m})$ with $\mu_{ij,m}\stackrel{i.i.d}{\sim} U[0,1]$. Define for each unit $i$ the $\ell$-th nearest‐neighbor distance by $ d_{(i,\ell),m} =$ the $\ell$-th smallest entry of row $i$ of $D^{*}_{m}$, excluding the diagonal. We then set $ d_{m}=\max_{1\leq i\leq n} d_{(i,\ell),m}, $ so that every unit has at least $\ell$ neighbors within distance $d_{m}$. Collecting these over $m=1,\dots,M$ yields the $M\times1$ vector $d_{m}$. We set $\delta=0$ for the null hypothesis and $\delta(\mathrm{x})=1-\mathrm{x}^2$ for any $\mathrm{x}\in R$ for the alternatives. We choose polynomials $\psi_j(\mathrm{x})=\mathrm{x}^j$ for $j=1,\cdots, h$ and the design of trigonometric functions is in the online appendix. The IV matrix for $\mathbb{W}_{1}$ and $\mathbb{W}_{2}$ is $K_{1}=[X, W_{1}X,..., W_{d_{\lambda}}X,\Psi]$. Furthermore, asy-p represents the standard-normal p-values, whereas chi-p represents the chi-square calibration $(\chi^2_{d_{\alpha}}-d_{\alpha})/\sqrt{2d_{\alpha}}$, which is to improve testing performance in small samples. This may also help in empirical applications, although in both applications in the previous section critical values from either distribution yield the same results. We discuss size performance in this section and present power performance in the online appendix, since under global alternatives the power of the asymptotic tests tends to one as $n\to\infty$. Table (ref) reports the empirical sizes of $\mathbb{W}_{1}$, and Table (ref) reports those of the SHAC-corrected statistic $\mathbb{W}_{2}$, at nominal 1%, 5% and 10% levels. When $n=200$, $\mathbb{W}_{2}$ exhibits substantially larger size distortion than $\mathbb{W}_{1}$. The reason is: Theorem (ref) requires $n^{-1}\,d_{\xi}\to 0,$ whereas Theorem (ref) requires $n^{\eta-1}\,d_{\xi}\to 0$ with $\eta=3/7$. In our design $d_{\xi}=d_{\lambda}+d_{\beta}+d_{h}$ with $d_{\lambda}\in\{2,4\}$, $d_{\beta}=2$, $d_{h}\in\{2,4,8\}$, so $d_{\xi}/n$ is already small at $n=200$, but $n^{\eta-1}d_{\xi}^2$ is relatively large. As $n$ increases to 500 and 900, $n^{\eta-1}d_{\xi}^2$ shrinks and the size of $\mathbb{W}_{2}$ converges rapidly to its nominal levels. Table (ref)-(ref) also show a clear pattern across the polynomial order with $h\in\{2,4,8\}$. For both $\mathbb W_{1}$ and $\mathbb W_{2}$, the smallest sieve yields the most accurate sizes, the largest sieve the greatest distortion, particularly for $\mathbb W_{2}$ at $n=200$. This accords with the fact that larger $h$ increases $d_{\alpha}$ and hence $d_{\xi}=d_{\lambda}+d_{\beta}+d_{\alpha}$, which in turn slows the convergence rates $n^{-1}d_{\xi}\to0$ for $\mathbb W_{1}$ and $n^{\eta-1}d_{\xi}^2\to0$ for $\mathbb W_{2}\,$. A similar scenario obtains when increasing $d_{\lambda}$, since it also raises $d_{\xi}$ and worsens finite-sample size accuracy.

Nonparametric spatial matrix

For the nonparametric spatial matrix, we generate $g^{*}_{ij}=\Phi\left(-b_{i j}\right) \mathbf{1}\left(c_{ij}<0.1\right)$ if $i \neq j$, and $w_{1ii}=0$, where $\Phi(\cdot)$ is the standard normal cdf, $b_{ij}\stackrel{i.i.d}{\sim} U[-3,3]$, and $c_{i j} \stackrel{i.i.d}{\sim} U[0,1]$. From this construction, we ensure that $G^*$ is sparse with no more than $10\%$ elements being nonzero. Then, define $G=G^{*} / 1.2 \bar{\alpha}(G^{*})$, ensuring the existence of $(I-G)^{-1}$. For $\mathbb{W}_{3}$, we approximate elements in $G$ by $\widehat{g}_{ij}=\sum_{l=0}^{d_{\tau}}\tau_{l}e_{l}(d_{ij})=\sum_{l=0}^{d_{\tau}}\tau_{l}[d_{ij,m}^{*}]^{l}\mathbf{1}(d_{ij,m}^{*}<d_{m})$ if $i \neq j$ with $d_{\tau}=2,4$, where $d_{ij,m}^{*}$ and $d_{m}$ are discussed in Section (ref), and the IV matrix for $\mathbb{W}_{3}$ is $K_{2}=[X,\mathcal{E}_{1}X,...,\mathcal{E}_{d_{\tau}}X,\Psi]$, where $e_{l}(d_{ij})$ are the typical elements of $\mathcal{E}_{l}$. Table (ref) reports the empirical sizes of $\mathbb{W}_{3}$ at nominal 1%, 5% and 10% levels. When $n=200$, $\mathbb{W}_{3}$ with $d_{\tau}=2$ yields rejection rates very close to the nominal levels, whereas the larger basis $d_{\tau}=8$ produces mild over-rejection. As $n$ increases to 500 and 900, both choices of $d_{\tau}$ rapidly converge to their nominal sizes. This scenario is similar to that of $\mathbb{W}_{2}$ and originates from the rate condition $n^{\eta-1}d_{\tau}^2\to0$ in Theorem (ref). In other words, it also requires an accurate approximation of the SHAC estimator.

Varying spatial coefficient

We now consider varying $\lambda$ such that $\lambda=\lambda(z_i)$ and $d_{\lambda}=1,2$. When $d_{\lambda}=1$, we set $\lambda(z_i)=0.9\sin(\pi z_i)$, when $d_{\lambda}=2$, we set $\lambda_1(z_i)=0.6\sin(\pi z_i)$ and $\lambda_2(z_i)=0.3\sin(\pi z_i)$. Others are the same as those in Section (ref). For estimation, we use $\phi_{il}(z)=\frac{1}{h}[\frac{2}{\pi}\tanh(z_{i})]^{l}$ as the polynomial function and $\phi_{il}(z) =\frac{1}{h}\sin(\frac{z_i}{2l})$ as the trigonometric function for $l=1,...,h$. The IV matrix $K_{3}=K_{1}$ with $d_{\lambda}=1,2$. Table (ref) reports the rejection probabilities of $\mathbb{W}_{4}$ at nominal 1%, 5% and 10% levels when testing $\delta(z)$. At $n=200$, the smallest sieve $(h=2,\,d_{\lambda}=1)$ achieves sizes close to nominal, whereas larger sieves, especially $h=8,\,d_{\lambda}=2$, over‐reject noticeably. As $n$ increases to 900, size distortion vanishes rapidly. Table (ref) shows the sizes of $\mathbb W_{4}$ when testing $\lambda(z)$. In particular, at $n=200$ with $d_{\lambda}=2$ and $h=8$, the over‐rejection is especially severe. By $n=900$, size accuracy is largely restored. This pattern reflects the rate condition in Theorem (ref), $n^{\eta-1}d_{\gamma}^2 \to 0$, since both $h$ and $d_{\lambda}$ enter the total sieve dimension $d_{\gamma}$, slowing the approximation when they are large. Across all designs, the chi-p delivers empirical sizes closer to nominal levels than asy-p in small samples, while the difference dissipates as $n$ grows. Results for trigonometric functions are similar, and left for the online appendix.

Conclusion

We provide a machinery for conducting tests on nonparametrically varying coefficients using ideas from parametric testing, and with similar ease of implementation. We allow for spatial dependence and provide for testing nonparametrically varying spatial dependence coefficients. It seems reasonable to conjecture that the methodology is applicable in other settings where one may be interested in testing linear restrictions on nonparametric functions. These provide directions for future research.

{0.3ex}

spacing{1} \normalem
table[table omitted — 8,500 chars of source]
table[table omitted — 9,477 chars of source]
table[table omitted — 5,801 chars of source]
table[table omitted — 5,800 chars of source]
table[table omitted — 5,794 chars of source]
table[table omitted — 5,775 chars of source]