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.
71,289 characters · 12 sections · 31 citation commands
Functional Spatial Autoregressive Models
Spatial interdependence among units is an essential element in spatial data analysis. To incorporate spatial interactions into econometric analysis, researchers have extensively utilized the Spatial Auto-Regressive (SAR) model:
where $y_i$ denotes a scalar outcome, $w_{i,j}$ denotes a known spatial weight between $i$ and $j$, $x_i$ denotes a vector of explanatory variables, and $\varepsilon_i$ denotes an error term. The spatial lag term $\sum_{j=1}^n w_{i,j} y_j$ captures the spatial trend of the outcome variable in the neighborhood of $i$, and the scalar parameter $\alpha_0$ measures its impact. The usefulness of SAR modelling (ref) has been demonstrated in various empirical topics, including regional economics, local politics, real estate, crimes, etc. In addition, if we define the weight term $w_{i,j}$ based on social distance or friendship connections instead of geographic distance, then the SAR models can be utilized to analyze social network data, and their applicability is vast.
To further broaden the applications of SAR modelling, this study aims to extend (ref) to a functional SAR model where the dependent variable is a function defined on a common closed interval:
where $q_i: [0,1] \to \mathbb{R}$ denotes the outcome function of interest. Restricting the support to $[0,1]$ is a normalization. In particular, for empirical relevance, this study primarily focuses on the case in which $q_i$ is the quantile function for a scalar dependent variable of interest. Regression models involving functional variables have been widely studied in the literature of functional data analysis (FDA) for several decades (e.g., ramsay2005functional). Our model is essentially different from the existing ones in that we explicitly consider the simultaneous spatial interactions of the outcome functions.
As a motivating example, suppose we intend to investigate the impact of a regional childcare subsidy program in a given city on the age distribution of the city. The policy is likely to attract households with young children from other regions to benefit from the subsidy. Additionally, if childcare facilities and schools need to be newly constructed, inflows of other age groups can also be anticipated as workers. To obtain a comprehensive picture of the shift in the age distribution owing to the subsidy program in its entirety, it would be natural to consider a regression model in which the dependent variable represents the age distribution of each city, such as the quantile function. Meanwhile, when the size of the young population in a given city is in an increasing trend (no matter the cause), which serves as a driver of economic growth of the city, this might also lead to an influx of working-age population into the surrounding regions owing to the spatial spillover of economic activities. The proposed functional SAR model (ref) is able to account for such interdependency between the outcome functions of nearby spatial units.
In the literature, we are not the first to consider an SAR-type modelling in the functional regression context. zhu2022network proposed a social network model similar to ours in a time-series setting, where the response variable is a function of time. They assumed that only concurrent interactions exist at each moment such that the past and future outcomes of others do not affect the present outcome. Consequently, when fixed at each time point, their model can be reduced to the standard SAR model in (ref). In this regard, our model may be considered to be a generalization of theirs such that $\alpha_0(t,s) \neq 0$ is allowed for $t \neq s$ in general.
Another related modelling approach to ours is the SAR quantile regression (e.g., su2011instrumental, malikov2019under, ando2023spatial). When $q_i$ represents a quantile function, our model and theirs are conceptually similar in that both approaches can examine the distributional effects of explanatory variables on the outcome and the spatial interaction of outcomes in a unified framework. However, a fundamental distinction lies in that we consider a model in which each unit has its own unique quantile function as the dependent variable. Consequently, we can explicitly allow for each specific quantile value of an outcome to interact with other quantiles of others' outcomes. For instance, our model can investigate the impacts of median outcome of neighborhoods on a specific (say) 10 percentile value of own outcome. In the time-series context, dong2024functional consider the same type of interaction structure as above.
Notice that our model (ref) is characterized as a simultaneous integral equation system, and to the best of our knowledge, this type of modelling has not been investigated in the econometrics literature. To construct a consistent estimator for our model, the model space should be restricted such that the realized $q_i$'s are uniquely (in some sense) associated with the true parameters. We show that to establish this uniqueness property, as in the standard SAR model (cf. kelejian2010specification), the spatial effects $\alpha_0$ must be bounded within a certain range. In particular, we demonstrate that the tightness of the bound required for $\alpha_0$ depends on the smoothness of the outcome function.
To estimate the model parameters, we need to address the endogeneity issue arising from the simultaneous interaction among the outcome functions. Thus, we propose a regularized two-stage least squares (2SLS) estimator that is based on a series approximation of $\alpha_0(\cdot, s)$ at each evaluation point $s$. Under the availability of a sufficient number of instrumental variables (IVs) and regularity conditions, we prove that both the estimator for $\beta_0$ and that for $\alpha_0(t,s)$ are consistent at certain convergence rates and asymptotically normally distributed. Additionally, we develop a Wald-type test for assessing the presence of any spatial effects at each $s$. We show that the proposed test statistic asymptotically distributes as the standard normal after appropriate normalization. Furthermore, we discuss performing the estimation when the outcome functions are not fully observable on the entire interval $[0,1]$, but are only discretely observed, which is typical in most empirical situations. Our proposed estimator relies on a simple interpolation method, and we derive a set of conditions under which the estimator can achieve the same asymptotic properties as the infeasible counterpart.
As an empirical illustration, we investigate the determinants of age distribution in Japanese cities. Since many Japanese cities are currently rapidly aging, which has emerged as one of the central social problems in the country, understanding the mechanisms underlying the age structure of cities is crucial. Using recent government survey data, including the Census, we apply our estimation and testing method to 1883 Japanese cities. Here, the outcome function $q_i$ represents the quantile function of the age distribution in city $i$, and covariates $x_i$ include variables such as annual commercial sales, unemployment rate, number of childcare facilities, and others. Our results suggest that spatial interaction effects are extremely weak at quantiles close to the boundary points 0 or 1. This may not be surprising as all individuals are born at age 0 and have a life expectancy of approximately 100 years at maximum, resulting in little regional heterogeneity. In contrast, strong spatial effects are observed when both $t$ and $s$ are at approximately the ages of young working population, possibly indicating that economic activities and their spillovers are the main factors in shaping the spatial trend of age structure.
The remainder of this paper is organized as follows: In Section (ref), we formally introduce the model proposed in this study and discuss the condition under which it is well defined with a unique solution. In addition, focusing on the cases where the outcome function is a quantile function, we discuss the motivations and interpretation of such a modelling approach. In Section (ref), we describe our 2SLS method for estimating $\beta_0$ and $\alpha_0$. Thereafter, we study the asymptotic properties of the proposed estimator under a set of assumptions. In this section, we also propose a test statistic for testing the null hypothesis that $\alpha_0(t, s) = 0$ for $t \in \mathcal{I}$, and its asymptotic distribution is derived. In Section (ref), we present the results of Monte Carlo experiments to evaluate the finite sample performance of the proposed estimator and test. Section (ref) presents our empirical analysis on the age distribution of Japanese cities, and Section (ref) concludes the paper.
\paragraph{Notation}
For a natural number $n$, $I_n$ denotes an $n \times n$ identity matrix. For a function $h$ defined on $[0,1]$, the $L^p$ norm of $h$ is written as $||h||_{L^p} \coloneqq (\int_0^1 |h(s)|^p \text{d}s)^{1/p}$, and $L^p(0,1)$ denotes the set of $h$'s such that $||h||_{L^p} < \infty$. For a random variable $x$, the $L^p$ norm of $x$ is written as $||x||_p \coloneqq (\mathbb{E}|x|^p)^{1/p}$. For a matrix $A$, $|| A ||$ and $||A||_\infty$ denote the Frobenius norm and the maximum absolute row sum of $A$, respectively. If $A$ is a square matrix, we use $\rho_{\max} (A)$ and $\rho_{\min} (A)$ to denote its largest and smallest eigenvalues, respectively. In addition, $A^{-}$ is a symmetric generalized inverse of $A$. We write $a \lesssim b$ and $a \lesssim_p b$ if $a = O(b)$ and $a = O_P(b)$, respectively. Finally, we write $a \sim b$ when $a \lesssim b$ and $b \lesssim a$.
Suppose that we have data of size $n$: $\{(q_i, x_i, w_{i,1}, \ldots, w_{i,n})\}_{i = 1}^n$, where $q_i$ denotes a random outcome function of interest with the common support $[0,1]$, $x_i = (x_{i,1}, \ldots, x_{i,d_x})^\top \in \mathbb{R}^{d_x}$ denotes a vector of covariates including a constant term, and $w_{i,j} \in \mathbb{R}$ denotes the $(i,j)$-th element of an $n \times n$ pre-specified spatial weight matrix $W_n = (w_{i,j})_{i,j = 1}^n$. The value of each $w_{i,j}$ is determined non-randomly. As is the convention, we set $w_{i,i} = 0$ for all $i$ for normalization. Note that the spatial configurations of the units generally change with the sample size. Thus, the variables generally form triangular arrays, and model parameters depend on $n$ through spatial interactions. However, when there is no confusion, the dependence on $n$ is suppressed for notational convenience.
As shown in (ref), our working model is
where $\overline q_i$ denotes the spatial lag of the outcome function: $\overline q_i \coloneqq \sum_{j = 1}^n w_{i,j} q_j$. The unknown parameters to be estimated are $\alpha_0$ and $\beta_0 = (\beta_{0 1}, \ldots, \beta_{0 d_x})^\top$. For instance, in our empirical analysis, $q_i(s)$ denotes the $s$-th quantile of the age distribution in city $i$, and $\alpha_0(t,s)$ captures the impacts from the $t$-th quantile ages of neighborhood cities to the $s$-th quantile age of own city. For other examples, $q_i(s)$ could be the $s$-th quantile of the income distribution in city $i$, $s$-th quantile of the daily activity energy expenditure of person $i$, number of available bicycles at the bicycle-sharing station $i$ at time $s$, and so forth. Hereinafter, we assume that $q_i \in L^p(0,1)$ for some $2 \le p < \infty$ and that $\alpha_0 \in C[0,1]^2$, where $C[0,1]^2$ denotes the set of continuous functions on $[0,1]^2$.
Before turning to the estimation of $\alpha_0$ and $\beta_0$, we discuss the completeness of our model, that is, whether model (ref) can be characterized by a unique solution $(q_1, \ldots, q_n)$. As our model comprises a system of $n$ functional equations, the existence and uniqueness of the solution are non-trivial problems. If the system does not have or has multiple solutions, consistently estimating the model parameters without some ad hoc assumptions is generally impossible.
Let $Q(s) = (q_1(s), \ldots, q_n(s))^\top$, $X = (x_1, \ldots, x_n)^\top$, and $\mathcal{E}(s) = (\varepsilon_1(s), \ldots, \varepsilon_n(s))^\top$. Then, we can re-write (ref) in matrix form as
This expression suggests that our model is seen as a system of Fredholm integral equations of the second kind with kernel $\alpha_0(t,s)$. Defining $\overline \alpha_0 \coloneqq \max_{(t,s) \in [0,1]^2} |\alpha_0(t,s)|$, whose existence is ensured under the continuity of $\alpha_0$, assume the following:
Let us denote $\mathcal{H}_{n,p} \coloneqq \{H = (h_1, \ldots, h_n) : h_i \in L^p(0,1) \; \text{for all} \; i\}$, and define a linear operator $\mathcal{T}$ as
whose range is $\mathcal{H}_{n,p}$ under Assumption (ref). Thus, we can write $Q = \mathcal{T} Q + X\beta_0 + \mathcal{E}$. Then, denoting $\text{Id}$ to be the identity operator, if the inverse operator $(\text{Id} - \mathcal{T})^{-1}$ exists, the solution $Q$ of the system can be uniquely determined (as an element of $\mathcal{H}_{n,p}$) as $Q = (\text{Id} - \mathcal{T})^{-1}[X \beta_0 + \mathcal{E}]$.
The next proposition states that Assumption (ref) is sufficient for the existence of $(\text{Id} - \mathcal{T})^{-1}$ and uniqueness of $Q$.
The proof is straightforward. Under Assumption (ref), we have
for any $H \in \mathcal{H}_{n,p}$ by Minkowski's and Jensen's inequalities. This implies that $\mathcal{T}H \in \mathcal{H}_{n,p}$. As is well known, if the operator norm of $\mathcal{T}$ is smaller than one, $(\text{Id} - \mathcal{T})^{-1}$ exists, and we have the Neumann series expansion $(\text{Id} - \mathcal{T})^{-1} = \sum_{\ell = 0}^\infty \mathcal{T}^\ell$ converging in the operator norm (e.g., Theorem 2.14, Kress2014linear). It is immediate from (ref) that $\left\| \mathcal{T} H \right\|_{\infty, p} < 1$ follows for any $H$ such that $||H||_{\infty, p} = 1$, which yields the desired result.
When the spatial weight matrix is row-normalized such that $||W_n||_\infty = 1$, as is often the case in empirical applications, Assumption (ref) can be reduced to $\overline \alpha_0 < 1$, which somewhat resembles the solvability condition $|\alpha_0| < 1$ for the standard linear SAR model (ref).
The Neumann series expansion implies that $Q$ can be expressed as $Q = X\beta_0 + \mathcal{T}X\beta_0 + \mathcal{T}^2 X\beta_0 + \cdots + \mathcal{E} + \mathcal{T}\mathcal{E} + \mathcal{T}^2 \mathcal{E} + \cdots$, that is,
Hence, the marginal effect of increasing $x_{i,j}$ on $Q(\cdot)$ is obtained by
where $\bm{e}_i$ denotes the $i$-th column of $I_n$. This clearly shows that a change in $i$'s covariate affects not only the outcome of $i$ but also those of other units through the spatial interaction - the so-called spatial multiplier effect.
One of the situations in which model (ref) can be most nicely applied empirically would be when the outcome function $q_i$ represents the quantile function for the cumulative distribution function (CDF) of a variable of interest. In our empirical analysis, we study the determinants of the population pyramids of Japanese cities by employing the age quantile function of city $i$ as $q_i$.
Suppose that for each $i$ we can observe a random CDF $F_i$ for an outcome variable $y \in \mathcal{Y}_i \subseteq \mathbb{R}$ of interest. The quantile function of $y$ for $i$ is defined as $q_i(s) \coloneqq \inf\{y \in \mathcal{Y}_i : s \le F_i(y) \}$. In the FDA literature, models where the response variable represents a probability distribution have garnered significant attention, for example, petersen2016functional, han2020additive, yang2020quantile, yang2020random, ghodrati2022distribution, petersen2021wasserstein, chen2023wasserstein. For an excellent review on this topic, refer to petersen2022modeling. A common view in these studies is that performing a regression analysis directly in the space of CDFs (or densities) is often problematic. Hence, we should consider imposing a regression model on the quantile function (rather than on the CDF per se), as in yang2020quantile and yang2020random, enabling us to enjoy several analytically and interpretationally preferable properties as mentioned below.
First, quantile functions can be easily computed without considering the range boundaries, unlike CDFs. Second, the domains of CDFs are typically heterogeneous across individuals, whereas that of quantile functions is always the fixed interval $[0,1]$. Third, the least-squares regression of the quantile function can be nicely interpreted as a Wasserstein distance minimization problem.
More specifically for the third point, denoting $F_i^{\alpha,\beta}$ to be the CDF induced from the quantile function $q_i^{\alpha, \beta}(s) \coloneqq \int_0^1 \overline q_i(t) \alpha(t,s) \text{d}t + x_i^\top \beta(s)$, the squared $2$-Wasserstein distance between $F_i$ and $F_i^{\alpha, \beta}$ is obtained as\footnote{ Formally, the Wasserstein distance is a distance between two probability measures. We abuse the notation using CDFs in its arguments for ease of explanation. For more precise discussions on the properties of Wasserstein distance, see panaretos2020invitation, for instance. }
Thus, minimizing the mean squared Wasserstein distance with respect to $(\alpha, \beta)$: $\min_{\alpha, \beta} n^{-1} \sum_{i = 1}^n \mathcal{W}_2^2(F_i, F_i^{\alpha,\beta})$ is equivalent to performing a functional least squares regression based on model (ref). Note, however, that the resulting least squares estimator does not produce a consistent estimate of $(\alpha_0, \beta_0)$ because of the endogeneity of $\overline q_i$. To circumvent the endogeneity issue, we introduce a penalized 2SLS method in the next section.
It is also worth noting that the spatially lagged quantile $\overline q_i$ corresponds to the quantile function of the spatially weighted Fr\'{e}chet mean $\overline F_i$ of $\{F_1, \ldots, F_n\}$ at the location of $i$: $\overline F_i \coloneqq \arg\min_F \sum_{j=1}^n w_{i,j} \mathcal{W}_2^2(F, F_j)$, which is also referred to as the Wasserstein barycenter, with $w_{i,j} \ge 0$ and $\sum_{j=1}^n w_{i,j} = 1$. Rather than using $\overline q_i$, one might consider employing the quantile function of the spatially lagged CDF: $\sum_{j=1}^n w_{i,j} F_i$ as the spatial trend term. However, a linear mixture of CDFs is generally multimodal and does not inherit the shape properties of the original CDFs. In this regard, the weighted Fr\'{e}chet mean should be more representative and faithful as an indicator of the neighborhood trend. This point is also highlighted in gunsilius2023distributional in the context of synthetic control analysis.
We now discuss the estimation of $\alpha_0$ and $\beta_0$. Let $\{\phi_k: k = 1,2, \ldots\}$ be a series of basis functions, such as Fourier series, B-splines, and wavelets, such that we can expand $\alpha_0(t,s) = \sum_{k = 1}^\infty \phi_k(t) \theta_{0k}(s)$ for each $s$. Then, we have $\int_0^1 q_i(t) \alpha_0(t,s) \text{d} t = \sum_{k = 1}^\infty r_{i,k} \theta_{0k}(s)$, where $r_{i,k} \coloneqq \int_0^1 q_i(t) \phi_k(t) \text{d}t$. Hence, our model (ref) can be re-written as
where $\overline r_{i,k} \coloneqq \sum_{j=1}^n w_{i,j} r_{j,k}$, $u_i(s) \coloneqq \int_0^1 \overline q_i(t) \alpha_0(t,s) \text{d}t - \sum_{k = 1}^K \overline r_{i,k} \theta_{0k}(s)$, and $K \equiv K_n$ is a sequence of integers tending to infinity as $n$ grows. Note that this is just a multiple regression model having $K$ endogenous regressors $\overline r_i = (\overline r_{i,1}, \ldots, \overline r_{i,K})^\top$ with a composite error term $ \varepsilon_i(s) + u_i(s)$. Thus, we can resort to the 2SLS approach to estimate $\theta_0(s) = (\theta_{01}(s), \ldots, \theta_{0K}(s))^\top$ and $\beta_0(s)$ under the availability of a sufficient number of valid IVs for $\overline r_i$.
Suppose that we have an $L \times 1$ vector $z_{1,i}$ of IVs that are correlated with $\overline r_i$ but not with $\varepsilon_i$ such that $L \equiv L_n \ge K_n$. The choice of IVs will be discussed later. Further, let $z_i = (z_{1,i}^\top, x_i^\top)^\top$, $Z = (z_1, \ldots, z_n)^\top$, $\overline R = (\overline r_1, \ldots, \overline r_n)^\top$, $M_z = Z (Z^\top Z)^{-} Z^\top$, $M_x = X (X^\top X)^{-1} X^\top$, and $\overline R_x = (I_n - M_x) \overline R$. Then, our 2SLS estimator is defined as follows:
for a given evaluation point $s \in (0,1)$, where $S = M_z \overline R [ \overline R^\top M_z \overline R ]^{-} \overline R^\top M_z$, $\lambda \equiv \lambda_n$ is a non-negative regularization parameter tending to zero as $n$ increases, and $D$ denotes a $K$-dimensional matrix, which is positive semidefinite, symmetric, and satisfies $\rho_\text{max} (D) \lesssim 1$ uniformly in $K$. Once $\widehat \theta_n(s) = (\widehat \theta_{n1}(s), \ldots, \widehat \theta_{nK}(s))^\top$ is obtained, we can estimate $\alpha_0(\cdot,s)$ by
To recover the entire functional form of $\alpha_0(\cdot, \cdot)$ and $\beta_0(\cdot)$, we can repeat the described estimation procedure over a sufficiently fine grid on $[0,1]$.
In practice, the 2SLS estimator in (ref) would be rarely feasible because function $q_i$ can usually only be incompletely observed. For instance, we might only be able to observe the values of $q_i$ at finite points, ${q_i(s_{i,1}), \ldots, q_i(s_{i,m_i})}$. This is the case of our empirical analysis of the age distribution in Japanese cities. In this empirical analysis, we cannot access the complete age distribution for each city, but we only know the distribution up to every five-year age interval. In such a case, for example, we can apply a linear interpolation method to obtain an approximation of the entire functional form of $q_i$. Without loss of generality, suppose the observations are ordered in an increasing way: $s_{i,1} \le s_{i,2} \le \dots \le s_{i,m_i}$. Then, for each given $s \in [s_{i,l}, s_{i,l + 1}]$, we estimate $q_i(s)$ by
where $\omega_i(s) = (s_{i,l + 1} - s)/(s_{i,l + 1} - s_{i,l})$. When $s < s_{i,1}$ (resp. $s > s_{i,m_i}$), we can set $\widehat q_i(s) = q(s_{i,1})$ (resp. $\widehat q_i(s) = q(s_{i,m_i})$).
When $q_i$ is a quantile function, it is also typical that a finite sample $\{y_{i,1}, \ldots, y_{i,m_i}\}$ randomly drawn from $F_i$ is only available. In this case, a straightforward approach to estimate $q_i$ would be to perform a nonparametric kernel CDF estimation and invert the estimate. Alternatively, we can also use a simple interpolation method as described in yang2020random.
Letting $\widehat q_i$ be any estimator of $q_i$, compute $\widehat r_{i,k} \coloneqq \sum_{j = 1}^n w_{i,j} \int_0^1 \widehat q_j(t) \phi_k(t) \text{d}t$ and let $\widehat r_i = (\widehat r_{i,1}, \ldots, \widehat r_{i,K})^\top$. Now, the feasible version of (ref) is defined as
where $\widehat Q(s) = (\widehat q_1(s), \ldots, \widehat q_n(s))^\top$, $\widehat R = (\widehat r_1, \ldots, \widehat r_n)^\top$, $\widehat R_x = (I_n - M_x)\widehat R$, and $\widehat S = M_z \widehat R [ \widehat R^\top M_z \widehat R]^{-} \widehat R^\top M_z$. The estimator for $\alpha_0(\cdot, s)$ can be obtained by $\widetilde \alpha_n(\cdot, s) \coloneqq \sum_{k=1}^K \phi_k(\cdot) \widetilde \theta_{nk}(s)$.
To derive the asymptotic properties of our estimators, we first need to specify the structure of our sampling space. Following jenish2012spatial, let $\mathcal{D} \subset \mathbb{R}^d$, $1 \le d < \infty$ be a possibly uneven lattice, and $\mathcal{D}_n \subset \mathcal{D}$ be the set of observation locations, which may differ across different $n$. For spatial data, $\mathcal{D}$ would be defined by a geographical space with $d = 2$. Notably, $\mathcal{D}$ does not necessarily have to be exactly observable to us. For example, $\mathcal{D}$ is possibly a complex space of general social and economic characteristics. In this case, we can consider it to be an embedding of individuals in a latent space, instead of their physical locations.
Assumptions (ref)(i) and (ii) together imply that the number of interacting neighbors for each unit is bounded. We believe this is not too restrictive in practice.
Assumption (ref)(i) states that the covariates and instruments are constant. The same type of assumption as this has been often utilized in the literatures on spatial econometrics and many-IV estimation (e.g., kelejian2010specification, hausman2012instrumental). Note that this assumption is essentially equivalent to considering all stochastic arguments as being conditional on $\{z_i\}_{i=1}^n$. Assumption (ref) restricts the distribution of the error functions, which accommodates virtually any form of heteroscedasticity. We might be able to relax the independence assumption in (ii) to some weak dependence condition, but we introduce this for technical simplicity. The $s$ in (iii) is a given interior point of $[0,1]$ at which the estimation is performed.
Assumption (ref) imposes a set of conditions on the basis functions. The $L^2$-convergence rate of the approximation errors for various bases is discussed in belloni2015some, where it is shown that $\ell_K(s) \lesssim K^{-\pi}$ typically holds when $\alpha_0(\cdot, s)$ is a $\pi$-smooth function (i.e., H\"older class of smoothness order $\pi$).
To state the next assumption, we define the following matrices: $\bm{V}_n(s) \coloneqq \text{diag}\{\mathbb{E}[\varepsilon_1^2(s)], \ldots, \mathbb{E}[ \varepsilon_n^2(s)]\}$, $\Omega_{n,x}(s) \coloneqq \Psi_{n,x}^\top \bm{V}_n(s) \Psi_{n,x}/n$,
The next theorem gives the convergence rate of our estimator.
The proofs of Theorem (ref) and those presented below are somewhat similar in several parts to those in hoshino2022sieve, but for completeness, they are all presented in Appendix (ref). Theorem (ref)(i) shows that the coefficients of $x_i$ can be estimated at the root-n rate. Meanwhile, result (ii) indicates that the $L^2$-convergence rate of $\widehat \alpha_n(\cdot, s)$ is not standard owing to the potential weak identification and the presence of the penalty term $\lambda D$. We can observe a trade-off that the first term converges to zero quickly by selecting a large $\lambda$, while the second term can vanish if we select $\lambda$ diminishing at a sufficiently fast rate such that $\nu_{KL}/\lambda \to \infty$. It is clear that the order of $||\theta_0(s)||_D$ is bounded by $\sqrt{K}$. When $\theta_0(s)$ is a sparse vector or it is decaying in the order of basis expansion, $||\theta_0(s)||_D \lesssim 1$ might be possible.
Next, define $\sigma_{n, \lambda}(t, s) \coloneqq \sqrt{\bm{\phi}_K(t)^\top \Sigma_{n, r, \lambda}^{-1} \Omega_{n,r}(s) \Sigma_{n, r, \lambda}^{-1}\bm{\phi}_K(t)}$, $\Sigma_{n, r, \lambda} \coloneqq \mathbb{E} \overline R_x^\top M_z \mathbb{E} \overline R_x / n + \lambda D$, and $\Omega_{n,r}(s) \coloneqq \mathbb{E} \overline R_x^\top M_z \bm{V}_n(s) M_z \mathbb{E} \overline R_x /n$. Moreover, let
where $\widehat{\bm{V}}_n(s) \coloneqq \text{diag}\{\widehat \varepsilon_1^2(s), \ldots , \widehat \varepsilon_n^2(s) \}$, and $\widehat \varepsilon_i(s) \coloneqq q_i(s) - \overline r_i^\top \widecheck \theta_n(s) - x_i^\top \widehat \beta_n(s)$, where $\widecheck \theta_n(s)$ denotes the estimator of $\theta_0(s)$ obtained following (ref) with $\lambda$ set to zero. Then, the limiting distribution of our estimator can be characterized as in the following theorem.
In this section, we consider statistically testing the presence of spatial effects. Specifically, for each given $s$, we test the following null hypothesis:
where $\mathcal{I}$ denotes a non-degenerate sub-interval of $[0,1]$. Then, a natural test statistic for testing $\mathbb{H}_0$ would be the Wald-type statistic given as follows:
where the dependence of $T_n$ on $s$ is suppressed. To derive the asymptotic distribution of $T_n$ under $\mathbb{H}_0$, let $\Xi_n \coloneqq \Sigma_{n, r, \lambda}^{-1} \mathbb{E} (\overline R_x^\top Z / n) (Z^\top Z / n)^{-}$ and $\Phi_\mathcal{I} \coloneqq \int_\mathcal{I} \bm{\phi}_K(t) \bm{\phi}_K(t)^\top \text{d}t$. Further, define
which serve as the mean and variance of $T_n$, respectively.
Here, we introduce the following miscellaneous assumptions.
The next theorem characterizes the asymptotic distribution of our test statistic.
When $\mathbb{H}_0$ does not hold, the standardized test statistic $(T_n - \mu_n)/\sqrt{v_n}$ deviates to a positive value. Thus, considering Theorem (ref), we can reject $\mathbb{H}_0$ at the $100\alpha$% significance level if the realized value of $(T_n - \mu_n)/\sqrt{v_n}$ exceeds the upper $\alpha$-quantile of $\mathcal{N}(0, 1)$. To implement the test in practice, we need to consistently estimate $\mu_n$ and $v_n$, which can be easily performed by the sample analogue estimators, the definitions of which should be clear from the context. The consistency of these estimators is straightforward (refer to Lemmas (ref) and (ref) and Theorem (ref)(iii), (iv)).
Finally, it is important to notice that when $\mathbb{H}_0: \alpha_0(t,s) = 0$ is indeed true over the entire $[0,1]$, higher-order spatially-lagged covariates are not valid IVs, that is, for example, $\overline{\overline{x}}_i$ and $\overline q_i$ are not related to each other. Thus, basically, we need to prepare a sufficient number of IVs using only $\overline x_i$ and possibly its transformations in this case.
Finally, in this section, we examine the cases in which the outcome functions are only discretely observed, and they are linearly interpolated following (ref). Letting $s_{i,0} = 0$ and $s_{i,m_i + 1} = 1$ for all $i$, we introduce the following assumption.
Assumption (ref)(i) determines the overall precision of the linear interpolation approximation. For simplicity of discussion, it assumes that the values of the outcome function are (quasi) uniformly observed such that the distance of any two consecutive observations is of order $\kappa$. In addition, note that we treat each observation point as nonstochastic. Assumption (ref)(ii) requires that the outcome function is H\"{o}lder continuous with exponent $\xi$ for all $i$. This assumption may be somewhat restrictive, but similar assumptions are often considered in the FDA literature (e.g., crambes2009smoothing). Obviously, we need some form of continuity in order for the interpolation approximation to work.
The following theorem states that the approximation errors caused by the linear interpolation are asymptotically negligible if $\kappa^\xi$ is sufficiently small.
Under Assumption (ref), the approximation error $|\widehat q_i(s) - q_i(s)|$ is of order $\kappa^\xi$ uniformly in $s$. The condition $\sqrt{n} \kappa^\xi / \sqrt{\nu_{KL}} \to 0$ states that the interpolation error should shrink to zero faster than $n^{-1/2}$, similar to the basis approximation error $\ell_K(s)$. From this result, it is also straightforward to observe the asymptotic equivalence between the feasible Wald test $\widetilde T_n \coloneqq n \int_{\mathcal{I}} \widetilde \alpha_n^2(t,s) \text{d}t$ and $T_n$ presented in the previous section.
\paragraph{Performance of the 2SLS estimator}
In this section, we first examine the finite sample performance of the proposed 2SLS estimator. We consider the following three data-generating processes (DGPs) for the Monte Carlo experiments:
where
$\beta_{0j}(s) = 1 + 1.2\log(s + 1)$ for $j = 1,2,3$, $\beta_{0j}(s) = \exp(s) - 0.4$ for $j = 4,\ldots, 7$, $x_{i,j} \overset{IID}{\sim} \mathcal{N}(0, 1)$ for all $j$, and $\varepsilon_i(s) = \varepsilon_{1,i} + \sum_{j = 1}^4 s^{j/2} \varepsilon_{2,i,j}$ with $\varepsilon_{1,i} \overset{IID}{\sim} \mathcal{N}(0, 0.3^2)$ and $\varepsilon_{2,i,j} \overset{IID}{\sim} \mathcal{N}(0, 0.6^2)$ for all $j$. When estimating the model, an intercept term is also included. We randomly allocate $n$ units on the lattice of $n/20 \times 40$, where we consider two sample sizes: $n \in \{400, 1600\}$. The spatial weight matrix $W_n$ is defined according to the Rook contiguity with row normalization. Since these three DGPs satisfy the requirements in Assumption (ref), we can generate the outcome functions $Q$ using the Neumann series approximation: $Q \approx Q^{(L)} \coloneqq \sum_{\ell = 0}^L \mathcal{T}^\ell [X\beta_0 + \mathcal{E}]$, where $L$ is increased until $\max_{1 \le i \le n} |q_i^{(L)}(s) - q_i^{(L-1)}(s)| < 0.001$ is met for all $s$. For computing the integrals over $[0, 1]$, we approximate them by finite summations over 199 grid points: 0.005, 0.010, \ldots, 0.995.
For the choice of the basis functions $\{\phi_k\}$, we use the cubic B-splines. We examine two values for the number of the inner knots of the B-splines: $\text{\# knots} \in \{2,3\}$, corresponding to $K = 6$ and $7$, respectively, both of which are equally spaced in $[0,1]$. The IVs used are the first- and second-order spatial lags of $\{1, x_{i,1}, \ldots, x_{i,7}\}$. Note that because there may exist some units that have no neighboring units, the spatial lags of $1$ are not necessarily constants. For the penalty term $\lambda D$, we set $D = I_K$ (i.e., the ridge penalty) and attempt using four values for $\lambda = \lambda_c n^{-3/5}$ with $\lambda_c \in \{0.5, 1, 2, 3\}$. The number of Monte Carlo repetitions for each setup is set to 1000. Throughout, the evaluation point $s$ is fixed at $s = 0.5$.
The performance of the coefficient estimator $\widehat \beta_n$ is evaluated using the average bias (BIAS) and the average root mean squared error (RMSE):
where superscript $(r)$ means that the estimate is obtained from the $r$-th replicated dataset. Similarly, for the estimator $\widehat \alpha_n$ of the spatial effect, we evaluate the performance based on the BIAS and RMSE averaged over the 19 evaluation points $\{t_1, t_2, \ldots, t_{19}\}$ equally spaced on $[0,1]$:
Table (ref) summarizes the simulation results. Our main findings are as follows: First, the results suggest that our estimator works satisfactorily well for all scenarios. The RMSE values for estimating $\beta_0$ are approximately halved when the sample size is increased from 400 to 1600, which is consistent with our theorem. Meanwhile, the RMSE values for estimating $\alpha_0$ do not decrease significantly even when the sample size is increased. This result would be owing to the increased variances caused by employing a smaller penalty parameter $\lambda$ (recall that $\lambda \sim n^{-3/5}$). When comparing the results of the estimators with different $\lambda$ values, our results suggest that when the functional form of the spatial effect $\alpha_0$ is simple as in DGPs 1 and 2, using an estimator with a relatively large penalty is advisable in terms of RMSE. In contrast, when the functional form of $\alpha_0$ is complex as in DGP 3, the estimator with the smallest penalty outperforms the others, which should be a reasonable result. It seems that the number of inner knots has only minute impacts on the estimation performance.
\paragraph{Performance of the Wald test}
Next, we assess the finite sample performance of our test for the presence of spatial effects. In this analysis, we use the same DGP as given above to generate the data, with a slight modification on $\alpha_0$ in DGP 2. Specifically,
where $\varrho \in \{0, 0.1, 0.2\}$. The null hypothesis to be tested is $\mathbb{H}_0: \alpha_0(t,0.5) = 0$ for $t \in [0.1, 0.9]$. Thus, $\mathbb{H}_0$ holds true when $\varrho = 0$.
In Table (ref), we present the rejection frequency over 1000 Monte Carlo repetitions at the 10%, 5%, and 1% significance levels. The results for $\varrho = 0$ demonstrate that the size of our test is reasonably well-controlled, with at most 1--2% deviation from the nominal levels for most cases. When the spatial effect is mild in magnitude ($\varrho = 0.1$), the estimator with a smaller penalty ($\lambda_c = 0.5$) is not sufficiently powerful to detect the effect probably owing to its large estimation variance. However, as expected, the power of the test can be significantly improved by increasing the sample size. In the case of a stronger spatial effect ($\varrho = 0.2$), all tests exhibit nearly perfect power property for all sample sizes.
\paragraph{Simulations under discretely observed outcome functions}
Finally, we evaluate the performance of our estimator and test when the entire shapes of the outcome functions are not perfectly observed but their values are discretely observable at finite points. The DGPs investigated here are identical to those used previously. To recover the entire functional form of the outcome function for each unit, we use the linear interpolation method in (ref). For all units, we assume that $m$ pairs of points $\{(s_{i,j}, q_i(s_{i,j}))\}_{j=1}^m$ are observable, where $s_{i,j}$'s are uniformly randomly drawn from $[0,1]$, and $m$ is selected from two values $m \in \{15, 50\}$.
To save space, the simulation results are omitted here and provided in Tables (ref) and (ref) in Appendix (ref). From these tables, we can observe similar overall tendencies as those shown above. An interesting finding is that, although increasing $m$ from 15 to 50 improves the RMSE for most cases, there are some situations in which the estimator with a smaller $m$ achieves an even slightly better RMSE. Similarly, comparing the results when $m = 15$ with those when the outcome function is fully observable (those reported in Table (ref)), the former occasionally exhibits smaller RMSE values. We conjecture that these phenomena occurred because the linear interpolation “smoothed out” the original, potentially noisier, outcome function, leading to a reduction in estimation variance. A similar discussion can be found in imaizumi2018pca in a different but related context. In contrast, regarding the size property of the Wald test, the linear interpolation seems to introduce certain distortions. Unsurprisingly, these distortions can be somewhat mitigated if $m$ is large. Except when $\lambda_c = 0.5$, the test exhibits a satisfactory power for both values of $m$.
In this section, we apply the proposed estimator and test to analyze the determinants of the age distribution of Japanese cities. While this type of data has been regularly studied in the FDA literature (e.g., delicado2011dimensionality, hron2016simplicial, Bigot2017geodesic), there are few papers attempting a regression-based analysis. In recent decades, many rural Japanese cities have been facing a serious aging population, prompting them to plan campaigns to encourage young people from urban areas to settle in their cities. Thus, investigating the relationship between the regional socioeconomic characteristics and the age structure and the impact of neighborhood trend on it would be meaningful.
Our sample comprises all local municipalities (Shi-ku-cho-son) in Japan. The age distribution data for each city are taken from the 2020 Census. For the covariates to explain the age distribution, we use the ratio of agricultural, forestry, and fishery workers, number of hospital beds per capita, number of childcare facilities per capita, unemployment rate, logarithm of annual commercial sales, and logarithm of average residential landprice. All variables are as of the most recent year before 2020, and they are all publicly available.\footnote{ Landprice data: \url{https://www.lic.or.jp/landinfo/research.html}; all others: \url{https://www.e-stat.go.jp/en}. } In addition to these, we include five regional dummies.\footnote{ They correspond to each of the following: Hokkaido-Tohoku, Chubu, Kinki, Chugoku-Shikoku, and Kyushu-Okinawa regions. } After excluding the observations with missing items, the analysis is performed on 1883 municipalities. Table (ref) in Appendix (ref) summarizes the detailed definitions of the variables used and their basic statistics.
Our age distribution data are not complete; we only have information on the population size at five-year intervals (0 -- 4 years old, 5 -- 9 years old, and so forth). Therefore, when computing the quantile function for each city, we performed the linear interpolation as in (ref). In Figure (ref), we depict the obtained quantile functions for 20 randomly selected cities from our dataset. The figure clearly shows the existence of certain regional heterogeneity in age compositions except those close to the boundary points.
For estimation, we follow the same procedure as in the previous section with $K = 7$ (three inner knots) and $\lambda = 3n^{-3/5}$. The integrals are replaced by summations over 399 equally-spaced grid points on $[0,1]$. For the spatial weight, expecting that the impacts from demographic changes in large cities should be larger than those from small cities, we consider the following specification:
When city $i$ has no neighbors (e.g., islands), we set $w_{i,j} = 0$ for all $j$. The estimation is performed on nine quantile values: $s = 0.1$, $0.2$, \ldots, $0.9$.
To save space, the estimated coefficients $\beta_0(s)$ are presented in Figure (ref) in Appendix (ref). Our major findings from the figure are as follows: Interestingly, for all variables, the impacts on age distribution become prominent around the median ($s = 0.5$), suggesting the residential flexibility of this age group in response to the socioeconomic conditions of a city. The variables considered as indicators of urbanness, such as the commercial sales and the landprice, exhibit negative effects, contributing to population rejuvenation. As expected, cities with a higher rate of agricultural workers exhibit a significant aging trend. Both the number of hospital beds and childcare facilities positively affect age distribution, although the underlying mechanisms are unclear. It is important to recall that in this study, the covariates are treated as fixed, and their potential endogeneity is ignored. To interpret the obtained results as a causal relationship, addressing the endogeneity issue more carefully would be necessary.
The estimated spatial effect function is reported in Figure (ref). The figure includes nine panels, each corresponding to different $s$-values. In the figure, we also report the computed test statistic $(T_n - \mu_n)/\sqrt{v_n}$ for $\mathcal{I} = [0, 1]$. From these results, we can observe the following: First, the values of the test statistic suggest that the spatial effects exist significantly at all nine quantiles. However, when quantile $t$ of the neighbor is close to either of the boundary points 0 or 1, almost no or weak spatial effects are present. This seems reasonable considering Figure (ref); only a little regional heterogeneity in age distribution is present at these extreme quantiles. The spatial interaction effects become particularly strong when both $t$ and $s$ are approximately 0.2 -- 0.5, which roughly correspond to the ages of the younger working population. This result might suggest that the growth of economic activities and their spillovers play main roles in forming the spatial trend of age distribution. Notably, the impacts from these lower-to-middle quantile values somewhat persist even for higher quantile ages. This could be reflecting the indirect effects from positive interactions among younger age groups, rather than a direct causal relationship across different quantiles.
In this study, we developed a new SAR model for analyzing spatial interactions among functional outcomes. For estimation, we developed a penalized 2SLS estimator and established its asymptotic properties under certain regularity conditions. Additionally, we developed a method for statistically testing the presence of spatial interactions. To illustrate the effectiveness of our proposed method, we performed an empirical analysis focusing on the age distribution in Japanese cities.
An important potential limitation of our study is that, while we have treated the covariates as fixed variables to simplify the theoretical exposition, this approach essentially obscures the endogeneity issue underlying the covariates. For instance, in our empirical analysis, it might be reasonable to consider the unemployment rate as an endogenous variable correlated with unobserved regional factors affecting the age distribution as well. One way to mitigate the endogeneity issue would be to extend the current model to a panel data model with functional fixed effects, which should be a promising topic for future studies. Another important future work is how to perform the estimation and inference when $\alpha_0$ and $\beta_0$ are not significant, leading to a weak IV problem. We conjecture that the inclusion of additional moment conditions based on the distribution of the error term might be effective in addressing this issue, as in lee2007gmm. Several other issues that need future investigation include: data-driven selection of tuning parameters and developing methods for uniform inference on the functional parameters.