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.
88,603 characters · 19 sections · 38 citation commands
Network Autoregressive Models for Functional Panel Data
\textsuperscript{1}Melbourne Business School, University of Melbourne, 200 Leicester Street, Carlton, Victoria 3053, Australia. Email: [email removed] \\
\textsuperscript{2}Corresponding author. School of Political Science and Economics, Waseda University, 1-6-1 Nishi-Waseda, Shinjuku-ku, Tokyo 169-8050, Japan. Email: [email removed]
Keywords: functional data analysis, panel data, endogeneity, network autoregressive models, bike-sharing systems
The availability of functional data has been rapidly expanding across all fields of research, leading to a growing need for statistical tools that appropriately account for the unique characteristics of each type of functional data. In the analysis of socioeconomic data, there are at least two key aspects that should be addressed. The first is that an individual's decision or behavioral pattern may influence that of others through social networks -- interactions between individuals. The second is that individuals are intrinsically heterogeneous, even after controlling for observable characteristics -- unobserved heterogeneity of individuals. Therefore, analyzing socioeconomic functional data requires functional models that jointly capture both of these aspects, which is the aim of this study.
More specifically, to account for the interactions among units, we extend the network (or spatial) autoregressive (NAR) modeling approach to a functional response model. To address the unobserved individual heterogeneity, we introduce the functional fixed effects approach, given the availability of panel data. When the response variable is a scalar rather than a function, there already exists a vast body of studies investigating fixed-effect NAR models for panel data, such as
and its variants (e.g., yu2008quasi, lee2010estimation, lee2014efficient, kuersteiner2020dynamic, among others). Here, $Y_{it}$ is a scalar outcome, $w_{i,j}$ denotes a known weight term measuring the social or geographical proximity between units $i$ and $j$, $X_{it}$ is a vector of covariates, $f_{0i}$ represents a fixed effect specific to each $i$, and $\varepsilon_{it}$ denotes an error term. The term $\sum_{j=1}^n w_{i,j} Y_{jt}$ captures the local trend of the outcome variable in the neighborhood of $i$. Model (ref) is typically applied in fields such as health, real estate, transportation, education, and municipal data. However, with the increasing availability of functional data in these fields, such as real-time activity recognition, real-time population mobility and congestion patterns, and regional wealth distributions, scalar models like (ref) may fail to appropriately capture the complex nature of these interactions. For an illustrative example, we refer the reader to the empirical analysis in Subsection (ref).
The above discussion motivates us to extend (ref) to the following model: for $s \in [0,1]$,
where $Y_{it}$ represents the outcome function of interest, which may or may not be a smooth function of $s$, $\alpha_0$ is the interaction effect function, $\beta_0$ is a vector of functional coefficients, $ f_{0i}$ is the fixed effect function, and $\varepsilon_{it}$ is the functional error term with mean zero at each $s$. Here, $A(\cdot , s)$ denotes a known functional, whose functional form may differ according to the research interest.
Since the response variable is a function, we can consider various types of interaction patterns. The most typical form of interaction would be the "concurrent" interaction, where only the responses of others at the same evaluation point $s$ are influential. In this case, $A(\cdot, s)$ is a point-evaluation functional at $s$: $A(Y_{jt}, s) = Y_{jt}(s)$. When $s$ represents time, then the past outcome should affect the future outcome (but the converse should not), which motivates us to employ $A(Y_{jt}, s) = \int_0^1 Y_{jt}(u)\nu(u,s)\text{d}u$, where $\nu(u,s)$ is a user-chosen function that is non-negative, increasing in $u$ up to $s$, and $\nu(u,s) = 0$ for $u > s$. For other examples, if others' responses at all evaluation points are equivalently influential, we may use $A(Y_{jt}, s) = \int_0^1 Y_{jt}(u) \text{d}u$. \Copy{const_A}{In this case, $A(Y_{jt}, s)$ is independent of $s$.} These examples can be represented as an integral operator $A(Y_{jt}, s) = \int_0^1 Y_{jt}(u) \nu(u, s)\text{d}u$ with some kernel weight function $\nu(u, s)$. For example, in the case of the point-evaluation functional, we can set $\nu(u,s) = \delta(u - s)$, where $\delta$ denotes the Dirac delta function.
Here, we provide three typical empirical topics to which model (ref) would be nicely applied.
In model ((ref)), the parameters of our primary interest to be estimated are $\alpha_0$ and $\beta_0$. With the total time periods $T$ possibly large or small, we apply a first-differencing transformation to eliminate the individual fixed effects from the model. For the transformed model, rather than estimating $\alpha_0(s)$ and $\beta_0(s)$ pointwise at many different $s$-values separately, we approximate them using orthonormal basis expansions and estimate their entire functional forms jointly in a single estimation. Our proposed estimator is based on the generalized method of moments (GMM). Specifically, we first derive a set of moment conditions at each $s$, in a similar manner to lin2010gmm, and then integrate these conditions numerically over $s \in [0,1]$. These integrated moment functions define our GMM objective function, and the resulting estimator is referred to as the integrated-GMM estimator. Once $\alpha_0(s)$ and $\beta_0(s)$ are estimated, if necessary, we can estimate the fixed effects $f_{01}(s), \ldots, f_{0n}(s)$ simply by taking the individual-level mean of the residuals.
Note that in model (ref), the outcome functions appear on both the left- and right-hand sides, implying that it is formulated as a system of simultaneous functional equations. Depending on the true values of the functional parameters, the model may exhibit an explosive network interaction process, leading to inconsistency of the proposed estimator. Thus, we first derive a condition for the completeness of our model -- a condition ensuring that the model is uniquely characterized by a stable solution. We consequently show that the magnitude of network interactions must reside within a certain range.
Then, under the uniqueness condition on network interactions, along with some regularity conditions, we derive the convergence rates of the integrated-GMM estimators for $\alpha_0$ and $\beta_0$. In addition, we prove that the estimators are asymptotically normal at each evaluation point $s$. These theoretical results are numerically corroborated through a series of Monte Carlo experiments.
As an empirical illustration, we apply our method to demand analysis for a bike-sharing system in the San Francisco Bay Area, U.S. Using publicly available data from Bay Area Bike Share, we study spatial interactions in bike availability across 70 stations from May 2014 to August 2015. Our results reveal significantly positive spatial interactions in bike availability during the morning hours, while negative interactions emerge in the early evening. \Copy{nonFDA_intro}{ Importantly, when we estimate a non-functional model (ref) by averaging bike availability over time, the estimated $\alpha_0$ is negative and close to zero. This highlights that ignoring the functional nature of the data may lead to misleading policy implications. } Furthermore, we conduct an impulse response analysis to demonstrate how a reduction in bike availability at a given station propagates over time to nearby stations. These findings underscore the importance of spatial interactions in shared mobility services and demonstrate the practical applicability of our method.
Our paper relates to a broad range of theoretical and empirical literature. From a theoretical perspective, our study contributes to both the FDA literature and the network/spatial interactions literature by proposing a new model that connects the two. In this sense, one of the most closely related studies to ours is zhu2022network. They proposed a functional NAR model similar to (ref) but not in panel data settings. In contrast to zhu2022network, our GMM estimator requires neither parametric assumptions nor IID conditions for the disturbance term. This weak requirement arises from the fact that we treat the individual effects as parameters, whereas zhu2022network perform functional principal component analysis to control them based on some homogeneity condition. Moreover, they considered only a concurrent interaction case (i.e., $A(Y_{jt}, s) = Y_{jt}(s)$). As described earlier, the variable $s$ typically represents a time on some scale.\footnote{ In this case, both $s$ and $T$ represent time, but they capture different time scales. Typically, $s$ is used to describe short-term cycles, such as intraday stock price movements, or daily variation in temperature. In contrast, $T$ captures more macro-level trends and thus corresponds to a longer time scale. } The concurrent interaction rules out interactions even with immediate past outcomes and allows only strictly simultaneous interactions, which should limit the interpretability of the model. Computationally, our estimator can recover the full functional forms of the functional parameters in a single step, while the estimator in zhu2022network must be repeatedly applied at each evaluation point $s$.
On the empirical side, demand forecasting for bike-sharing systems has been an active topic in the data science literature (e.g., faghih2016incorporating, lin2018predicting, eren2020review, torti2021modelling, among others). Among these studies, faghih2016incorporating is most closely related to our study in that they employed a spatial panel model similar to (ref) to analyze the spatial and temporal interaction structure for the bike-sharing system in New York City, CitiBike. In their approach, however, the data are not treated as functional, and thus the model parameters are not allowed to vary over time. By contrast, torti2021modelling analyzed the flow of bikes in the bike-sharing system in Milan, BikeMi, through a functional linear model with functional coefficients; however, they did not account for the spatial interactions of mobility. Thus, our empirical study can be viewed as combining the strengths of these two papers.
\paragraph{Paper organization}
The rest of the paper is organized as follows. In Section (ref), we present our model and discuss its uniqueness condition. Section (ref) introduces our integrated-GMM estimator and investigates its asymptotic properties. In Section (ref), we discuss additional topics related to our model, including the estimation of marginal effects and network impulse response analysis. Section (ref) conducts a set of Monte Carlo simulations to numerically demonstrate the properties of our estimator. Section (ref) discusses several extensions that should be useful in empirical applications. Section (ref) presents our empirical analysis on the U.S. bike-sharing data, and Section (ref) concludes. \Copy{appendix}{Proofs of technical results, additional empirical material, and the programming codes used in the numerical studies are provided in the online supplementary material.}
\paragraph{Notation}
For a function $h$ defined on $[0,1]$ and $p \in [1,\infty)$, 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 $M$, $|| M ||$, $||M||_1$, and $||M||_\infty$ denote the Frobenius norm, the maximum absolute column sum, and the maximum absolute row sum of $M$, respectively. If $M$ is a symmetric matrix, we use $\lambda_{\max} (M)$ and $\lambda_{\min} (M)$ to denote its largest and smallest eigenvalues, respectively. For a positive integer $Z$, we denote $[Z] \coloneqq \{1, \ldots, Z\}$. We use $I_Z$ to denote an identity matrix of dimension $Z$. Finally, $X \lesssim Y$ if $X = O(Y)$ almost surely, and $X \lesssim_P Y$ if $X = O_P(Y)$.
Suppose that we have balanced panel data of size $(n, T)$: $\{(Y_{it}, X_{it}, w_{i,1}, \ldots, w_{i,n}): i \in [n], \: t \in [T]\}$. The number of time periods $T \ge 2$ can be either fixed or tending to infinity jointly with the sample size $n$. Here, $Y_{it}: [0,1] \to \mathbb{R}$ denotes a random outcome function of interest with the common support $[0,1]$, $X_{it} = (X_{it}^1, \ldots, X_{it}^{d_x})^\top$ denotes a vector of covariates, and $w_{i,j} \in \mathbb{R}$ is the $(i,j)$-th element of an $n \times n$ time-invariant interaction matrix $W_n = (w_{i,j})$. \Copy{Y}{Throughout, we assume that the functional outcome $Y_{it}(s)$ is fully observed over $s \in [0,1]$. In our empirical application, $Y_{it}(s)$ represents the number of available bikes at station $i$, where $s$ denotes the time of day. The data record the timing of all pickups and returns at each station, which implies the complete observability of $Y_{it}$ as a function of $s$. As in the case of bike availability, we do not rule out situations in which $Y_{it}(s)$ is discontinuous in $s$. } \Copy{X}{The covariates $X_{it}$ do not include a constant term and vary with both $i$ and $t$.} The value of each $w_{i,j}$ is pre-determined non-randomly. In social network analysis, it is common to set $w_{i,j} = c_{i,j} \bm{1}\{\text{$i$ and $j$ are peers}\}$, where $c_{i,j}$ is some normalizing constant. Similarly, if each $i$ represents a spatial unit, one may use $w_{i,j} = c_{i,j} \bm{1}\{\Delta(i,j) \le \overline \Delta\}$, where $\Delta(i,j)$ is the distance between $i$ and $j$, and $\overline \Delta$ is a given threshold. As is the convention, we set $w_{i,i} = 0$ for all $i$ for normalization.
As shown in (ref), our working model is given as follows: for $s \in [0,1]$,
where $\overline Y_{it} = \sum_{j=1}^n w_{i,j} Y_{jt}$. Considering analytical transparency and interpretational simplicity in empirical applications, we assume $A(\cdot, \cdot)$ is linear in its first argument so that we have $\sum_{j=1}^n w_{i,j} A(Y_{jt}, \cdot) = A(\overline Y_{it}, \cdot)$. While this linear specification provides a reasonable starting point, more complicated forms are studied in Subsection (ref).
The parameters of primary interest are the interaction effect function $\alpha_0(s)$ and the coefficient functions $\beta_0 (s)= (\beta_{01}(s), \ldots, \beta_{0d_x}(s))^\top$. The functional individual effects $f_{01}(s), \ldots, f_{0n}(s)$ are treated as nuisance parameters. Restricting the support of $s$ to be a unit interval is a normalization, which is harmless as long as the response functions have the identical interval support. For simplicity, we do not explicitly assume that $X_{it}$ is a function of $s$, which can be relaxed easily at the expense of more complicated notation and proofs. In the rest of this section, we assume that $Y_{it} \in L^2(0,1)$, $\varepsilon_{it} \in L^2(0,1)$, and that $\alpha_0$ and $\beta_0$ are continuous functions. To derive the asymptotic properties of our estimator, we impose stronger regularity conditions later.
We discuss the completeness of our model (in the sense of tamer2003incomplete, lewbel2007coherency). That is, since the model is a system of simultaneous functional equations, it may not have a unique interior solution in general, depending on the parameter values. Thus, we need to impose some conditions on the model to ensure that the outcome functions follow a unique data-generating process and prevent explosive behavior.
Let $Y_t(s) = (Y_{1t}(s), \ldots, Y_{nt}(s))^\top$, $\bm{A}(Y_t, s) = (A(Y_{1t}, s), \ldots, A(Y_{nt}, s))^\top$, $X_t = (X_{1t}, \ldots, X_{nt})^\top$, $F_0(s) = (f_{01}(s), \ldots, f_{0n}(s))^\top$, and $\mathcal{E}_t(s) = (\varepsilon_{1t}(s), \ldots, \varepsilon_{nt}(s))^\top$. Then, we can re-write (ref) in matrix form as
This expression clearly indicates that our model is characterized as $T$ independent systems of functional equations of size $n$. We introduce the following assumption to ensure that the model has a unique stable solution.
Assumption (ref)(i) requires that the magnitude of the network interaction is not too strong. With Assumption (ref)(ii), we have for any $h,h' \in L^2(0,1)$ that $||A(h - h',\cdot)||_{L^2} \le ||h - h'||_{L^2}$, which indicates the non-expansive property of the operator $A$. This assumption still accommodates many empirically interesting interaction patterns. For example, in the case of the point-evaluation functional $A(h,s) = h(s)$, it trivially satisfies $||A(h,\cdot)||_{L^2} = ||h||_{L^2}$. For another example, suppose $A(h,s) = \int_0^1 h(u) \nu(u,s) \text{d}u$ for some continuous $\nu$. Since
where $\overline \nu \coloneqq \max_{u,s \in [0,1]^2} |\nu(u,s)|$, Assumption (ref)(ii) is implied if $\overline \nu \le 1$ holds.
Now, denote $\mathcal{H}_{n,p} \coloneqq \{H = (h_1, \ldots, h_n) : h_i \in L^p(0,1), \; i \in [n]\}$, and define a linear operator $\mathcal{A}$ as
Then, we can write our model symbolically as follows:
Further, denoting $\text{Id}$ to be the identity operator, if the inverse operator $(\text{Id} - \mathcal{A})^{-1}$ exists, the solution $Y_t$ of the system can be uniquely determined up to an equivalence class in $\mathcal{H}_{n,2}$ as
The next proposition states that Assumption (ref) is sufficient for the existence of $(\text{Id} - \mathcal{A})^{-1}$.
The proof of Proposition (ref) is given in Appendix (ref). Note that the explicit form of the inverse operator $(\text{Id} - \mathcal{A})^{-1}$ cannot be derived in general, except for some simple cases such as $(\mathcal{A} Y_t)(s) = \alpha_0(s) W_n Y_t(s)$. In this case, $(\text{Id} - \mathcal{A})^{-1}$ is obtained as $(I_n - \alpha_0(\cdot) W_n)^{-1}$. However, in practice, we can approximate it with arbitrary precision by truncating the Neumann series expansion $(\text{Id} - \mathcal{A})^{-1} = \sum_{\ell = 0}^\infty \mathcal{A}^\ell$ at a sufficiently large order (see, e.g., Kress2014linear). See Remark 1 in zhu2022network for a related discussion.
\Copy{explosive_alpha}{ Finally, we note that Assumption (ref) is only one type of sufficient condition for the existence of $(\text{Id} - \mathcal{A})^{-1}$. Condition (i) can be relaxed by allowing $\alpha_0(s)$ to exhibit locally explosive network dependence for some $s \in [0,1]$, provided that the operator $A$ satisfies a condition more restrictive than condition (ii); see Subsection (ref) in Appendix (ref) for details. Accordingly, unlike standard non-functional NAR models, the completeness of our model can still be satisfied even when the network effect $\alpha_0(s)$ takes locally extreme values. That said, we adopt the current form of Assumption (ref) for technical convenience in the subsequent analysis. }
To estimate the unknown functional parameters $\alpha_0$ and $\beta_0$, there are broadly two approaches. The first is a "local" approach that estimates the values of these functions at specific $s$-values, repeating the estimation across different points to recover the full functional forms. The second is a "global" approach that estimates the entire functional forms in a single step using a series approximation method. Although both approaches are theoretically valid, the local approach typically requires more computation time and often leads to larger variance (but smaller bias) because it does not exploit information from nearby evaluation points. This study adopts the global approach.
Let $\{\phi_k: k = 1,2, \ldots\}$ be a series of orthonormal basis functions. We throughout assume that $\phi_k$'s are continuous on $[0,1]$. Then, if the functions $\alpha_0$ and $\beta_0$ are sufficiently smooth, we can approximate
uniformly in $s \in [0,1]$, for some coefficient vectors $\theta_{0\alpha} = (\theta_{0\alpha, 1}, \ldots, \theta_{0\alpha, K})^\top$ and $\theta_{0j} = (\theta_{0j, 1}, \ldots, \theta_{0j, K})^\top$, $j \in [d_x]$. Here, $K \equiv K_{nT}$ is a sequence of positive integers tending to infinity as $nT$ increases. For simplicity of presentation, we use the same basis function $\phi_k$ and the same basis order $K$ to approximate both $\alpha_0$ and $\beta_0$. Define $\theta_0 = (\theta_{0\alpha}^\top, \theta_{01}^\top, \ldots, \theta_{0d_x}^\top)^\top$, $\phi^K(s) = (\phi_1(s), \ldots, \phi_K(s))^\top$,
Then, we can further re-write the model in (ref) as
Here, $V_t(s) = (v_{1t}(s), \ldots, v_{nt}(s))^\top$ is an $n \times 1$ vector of series approximation errors:
Under the assumptions we will introduce, this approximation error diminishes to zero at a certain rate as $K$ goes to infinity. How to choose an appropriate $K$ will be discussed in Remark (ref).
Further, let
and $\underset{N \times nT}{\bm{D}} = (d_{ij})$ be the one-period lag operator, whose $(i,j)$-th element is defined as
Then, we can remove the unknown fixed effects from the model in the following manner:
We estimate $\theta_0$ based on this expression. In order to consistently estimate $\theta_0$, we need to address the endogeneity issue caused by the simultaneous interactions of the response functions; that is, since $A(\overline Y_{it}, s)$ is correlated with the error term $\varepsilon_{it}(s)$ in general, simply regressing $\bm{D} \bm{Y}(s)$ on $\bm{D} \bm{H}(s)$ does not yield a consistent estimate of $\theta_0$. To tackle this issue, we employ an instrumental variable (IV) approach.
Suppose we have a $d_q \times 1$ vector of IVs $Q_{it} = (Q_{it}^1, \ldots, Q_{it}^{d_q})^\top$ for $A(\overline Y_{it}, s)$. Define
Then, we have the linear moment conditions
\Copy{IVchoice}{ Noting that $W_n Y_t = W_n \mathcal{A} Y_t + W_n X_t \beta_0 + W_n F_0 + W_n \mathcal{E}_t$, we can find that the network lagged covariates $\overline X_{it} \coloneqq \sum_{j = 1}^n w_{i,j} X_{jt}$ (and also their lags) are valid candidates for $Q_{it}$. However, when $X_{it}$ is locally homogeneous, $\overline X_{it}$ and $X_{it}$ become strongly correlated, which leads to a weak IV problem. Moreover, even when they are not strongly correlated, if the impact of $X_{it}$ on $Y_{it}$ is small (i.e., $\beta_0 \approx 0$), a weak IV problem may still occur. }
Although one can estimate $\theta_0$ based on the linear moment conditions only, which results in a two-stage least squares (2SLS) type estimator, we can utilize additionally the quadratic moment conditions to improve the efficiency of estimation (see, e.g., lin2010gmm). That is, under the independence assumption on the error terms $\{\varepsilon_{it}(s)\}_{i \in [n], t \in [T]}$ (Assumption (ref)(i) below), for any $N \times N$ matrices $P_m \coloneqq I_{T-1} \otimes P_{m,1}$, where $P_{m,1}$ ($m = 1, \ldots , M$) is an $n \times n$ matrix whose diagonal elements are all zero, we have\footnote{ \Copy{M}{ Throughout the paper, we assume that $M$ is fixed. However, as pointed out by a referee, it would be interesting to consider an increasing number of quadratic moment conditions by setting, for example, $P_{m,1} = W_n^m - \text{diag}(W_n^m)$. The quadratic moments constructed in this way are strongly correlated with each other, and the contribution of higher-order quadratic moments may be marginal. This setting corresponds to GMM with many potentially weak moments (e.g., han2006gmm), for which a separate and careful discussion is required. Thus, we leave this issue for future research. } }
Some examples of $P_{m,1}$ include $P_{m,1} = W_n$ and $P_{m,1} = W_n^\top W_n - \text{diag}(W_n^\top W_n)$. \Copy{quadratic}{ The use of quadratic moment conditions is particularly important when the weakness of the IV $\bm Z(s)$ is a concern. As discussed in yang2025estimation and also demonstrated in our numerical results, quadratic moment conditions can help mitigate weak identification when $\bm Z(s)$ are not strong IVs. }
Combining the linear and quadratic moment conditions, we can construct our estimator based on the following $d_g \coloneqq (d_q + d_x)K + M$ moment conditions: for $s \in [0,1]$,
As the empirical counterpart of these moment conditions, given a candidate value $\theta$ for $\theta_0$, we define
where $\bm{E}(s; \theta) \coloneqq \bm{Y}(s) - \bm{H}(s)\theta$. \Copy{continuum}{ Since we have a set of moment conditions that holds continuously over $s \in [0,1]$, to construct our GMM estimator, we numerically integrate these conditions. Specifically, for a set of pre-specified $L \equiv L_{nT} \ge K$ grid points in $[0,1]$, denoted by $0 \le s_1 \le \cdots \le s_L \le 1$, we define
} Remark (ref) below provides practical guidance on how to choose the grid $\{s_l\}$. Situations where the response functions are not fully observable are discussed in Subsection (ref).
Now, we are ready to introduce our estimator:
$\Omega_{nT}$ is a $d_g \times d_g$ positive definite symmetric weight matrix. \Copy{compactness}{For each $K \ge 1$, $\Theta_K \subset \mathbb{R}^{(d_x + 1)K}$ is a compact parameter space containing $\theta_0$ in its interior.} Moreover, we assume that $\phi^K(s)^\top \theta_\alpha$ and $\phi^K(s)^\top \theta_j$ are bounded for any $\theta = (\theta_\alpha^\top, \theta_1^\top, \ldots, \theta_{d_x}^\top)^\top \in \Theta_K$, uniformly in $s \in [0,1]$. For one example of the weight matrix, we can use $\Omega_{nT} = I_{d_g}$. Another example is the 2SLS-type weight matrix:
Once $\widehat \theta_{nT}$ is obtained, the estimators of $\alpha_0(s)$ and $\beta_0(s)$ are given as
which we refer to as the integrated-GMM estimators. Additionally, if one is interested in the estimation of individual fixed effect functions, the following estimator can be used:
where $\widehat \beta_{nT}(s) = (\widehat \beta_{nT,1}(s), \ldots, \widehat \beta_{nT,d_x}(s))^\top$. Consistent estimation of $f_{0i}(s)$ by $\widehat f_{ni}(s)$ requires $T$ to increase to infinity, while $\alpha_0(s)$ and $\beta_0(s)$ can be consistently estimated even when $T$ is fixed. Unlike $\widehat \alpha_{nT}(s)$ and $\widehat \beta_{nT}(s)$, the estimator $\widehat f_{ni}(s)$ is not necessarily continuous, as we do not preclude cases where $Y_{it}(s)$ and $A(\overline Y_{it}, s)$ are discontinuous in $s$.
To derive the asymptotic properties of our estimator, we first need to specify the structure of our sampling space. Following jenish2012spatial, let $\mathcal{D} \subset \mathbb{R}^d$ be a possibly uneven lattice, and $\mathcal{D}_n \subset \mathcal{D}$ be the set of observation locations. Once the observation locations are determined for a given sample of size $n$, we assume that they do not vary over time $t$. For spatial data, $\mathcal{D}$ would be defined by a geographical space with $d = 2$.\footnote{ Note that $\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, rather than their physical locations. }
We first derive the rates of convergence of our estimator under the following set of assumptions.
Assumptions (ref)(i) and (ii) together imply that the number of interacting partners for each unit is bounded (i.e., the network must be sparse). These assumptions play a crucial role in characterizing the stochastic process of the outcome functions. In Assumption (ref), part (i) assumes that the covariates are non-stochastic and bounded. This type of assumption is frequently utilized in the spatial and network literature and can be interpreted as viewing the analysis conditional on the realized values of the covariates. Meanwhile, part (ii) is introduced to ensure some convergence results for the quadratic moments.
Assumption (ref)(i) allows the error terms to be fully heteroskedastic. Part (ii) should be standard. Part (iii) is a high-level condition, which plays an important role to obtain the parametric convergence rate for the GMM estimator. This assumption essentially requires the covariance operator $\Gamma_{it}$ to be trace class (see, e.g., hsing2015theoretical). More specifically, if $\Gamma_{it}$ belongs to $L^2([0,1]^2)$, it admits the following series expansion:
for some sequence of constants $\{ \kappa_{it,k_1,k_2}\}$. By the orthonormality of $\phi_k$,
Since $L^{-2} \sum_{l=1}^L \sum_{l'=1}^L \Gamma_{it}(s_l, s_{l'}) \phi_k(s_l) \phi_k(s_{l'})$ can be seen as a numerical approximation of the left-hand side of the above expression, Assumption (ref)(iii) essentially requires that $\sum_{k=1}^K \kappa_{it,k,k} \lesssim 1$ uniformly in $K$. This type of assumption is commonly made in the FDA literature.
Assumption (ref) is not restrictive in most empirically relevant situations. For example, in the case where $A(h,s) = h(s)$, we can set $\omega_p(u,s) = \delta(u - s)$ for any $p$. For another example, when $A(h,s) = \int_0^1 h(u) \nu(u, s) \text{d}u$ for some kernel $\nu(u, s)$, since $|A(h,s)|^p \le \int_0^1 |h(u)|^p |\nu(u, s)|^p \text{d}u$, we can set $|\nu(u, s)|^p = \omega_p(u,s)$ in this case.
In Assumption (ref), we assume that the weight matrices in the quadratic moments are symmetric. Note that this assumption does not lose any generality because $A^\top P_{m,1} A = A^\top P_{m,1}^\top A$ for any $n \times 1$ vector $A$. If $P_{m,1}$ is not symmetric in practice, we can always symmetrize it as $(P_{m,1} + P_{m,1}^\top)/2$. The assumption for the existence of a threshold distance $\overline \Delta_m$ may be non-standard, but it simplifies the proof. Since $P_{m,1}$'s are usually created from the interaction matrix $W_n$ and its powers, this assumption is consistent with Assumption (ref)(ii).
Assumption (ref) is a regularity condition to ensure the identifiability of $\theta_0$. \Copy{pi}{ Assumption (ref) is standard in the series estimation literature. As shown in Theorem (ref), we require $\pi > 3/2$ in order to establish the asymptotic normality of the estimator. This condition is satisfied, for example, when spline basis functions are used and the functions $\alpha_0$ and $\beta_{0j}$ are twice continuously differentiable (see, e.g., CHEN20075549, belloni2015some).
}
The proof of Theorem (ref) is given in Appendix (ref). Result (i) shows that, if the functional parameters are sufficiently smooth so that $K^{-\pi} \lesssim (nT)^{-1/2}$, the series coefficient estimator is consistent and converges at the parametric rate. The same convergence rate applies to the $L^2$-convergence rate of the functional estimators, as shown in (ii) and (iii). \Copy{parametric_rate}{ Although these results may appear somewhat unexpected since the dimension of $\theta_0$ diverges as $K \to \infty$, parametric convergence rates are commonly observed in dense functional data settings (e.g., huang2004polynomial, li2010uniform, hu2021robust), our setup corresponding to an extreme case in which the entire functions are observed. The key condition underlying this rate is Assumption (ref)(iii), which ensures that, once $K$ is sufficiently large, adding further basis terms does not inflate the variance of the series coefficient estimator.} The uniform convergence rate for these estimators is $K^{1/2}$ slower than the $L^2$-convergence rate.\footnote{ Note that the uniform convergence results obtained here are not necessarily the sharpest, and the theoretically optimal convergence rates under our setup are also unknown. These points are left for future research. }
We next present the limiting distribution of our estimators. To this end, we introduce the following notations and additional assumption:
More explicit forms of the matrices $J_{nT}(s; \theta)$ and $\mathcal{V}_{nT}$ can be found in (ref) and in (ref) in Appendix (ref), respectively. Further, let $\mathbb{S}_\alpha$ and $\mathbb{S}_j$ be the $K \times (d_x + 1)K$ selection matrices such that $\theta_{0\alpha} = \mathbb{S}_\alpha \theta_0$ and $\theta_{0j} = \mathbb{S}_j \theta_0$ hold.
Theorem (ref) and Corollary (ref) establish the pointwise asymptotic normality of the integrated-GMM estimators. We impose additional undersmoothing conditions to ensure that the bias terms vanish sufficiently quickly. See also Proposition (ref) in Appendix (ref), where additional asymptotic results for consistent variance estimation are established.
For the interpretation of these results, there are several points that should be noted. \Copy{grid_choice}{ First, the standard deviations $\sigma_{nT,\alpha}(s)$ and $\sigma_{nT,j}(s)$ depend on the choice of the grid $\{s_l\}$ and its size $L$, as shown in their definitions. In this sense, one could consider a data-driven procedure for selecting the grid so as to improve the efficiency. We leave the investigation of such procedures for future research. } \Copy{root_n_normal}{ Second, a straightforward calculation shows that the order of the standard deviations and their estimators is governed by $\|\phi^{K}(s)\|$, which is typically of order $\sqrt{K}$ for $s \in (0,1)$. Therefore, the theorem does not necessarily assert the $\sqrt{N}$-asymptotic normality for our functional estimator (similar to Theorem 2 of newey1997convergence, for example). } \Copy{boundary}{ Third, the asymptotic normality results hold for any evaluation point $s \in [0,1]$, including the boundary points $s = 0$ and $s = 1$. Nevertheless, both the required conditions on $K$ and the structure of the standard deviations depend on $s$ through the behavior of $\phi^{K}(s)$. For example, consider $\phi^{K}(0) = (1,0,\ldots,0)^\top$, as is typical in practice. In this case, $\|\phi^{K}(0)\|^{2} = 1$. Since Theorem (ref) requires $K^{2}/(\sqrt{nT}\|\phi^{K}(s)\|^{2}) \to 0$, this implies that achieving the normality at the boundary needs a more restrictive condition on $K$ than interior points. }
Once the model is estimated, as a next step, one might be interested in computing the marginal effects of covariates on the outcome. In a standard linear regression model without network interaction, the estimated coefficients directly represent the marginal effects of their corresponding covariates. However, in the presence of intricate functional interaction, this is no longer the case.
As shown in Section (ref), under Assumption (ref), we have the following moving-average type representation:
This expression indicates that the marginal effect of increasing $X_{it}^j$ by one unit on $Y_t$ is given by $\partial Y_t/(\partial X_{it}^j) = \lim_{c \to 0} \sum_{\ell = 0}^\infty [\mathcal{A}^\ell (X_t^j + \bm{e}_i c) \beta_{0j} - \mathcal{A}^\ell X_t^j \beta_{0j}]/c = \sum_{\ell = 0}^\infty \mathcal{A}^\ell \bm{e}_i \beta_{0j}$ by the linearity of $\mathcal{A}^\ell$, where $X_t^j$ is the $j$-th column of $X_t$, and $\bm{e}_i$ denotes the $i$-th column of $I_n$. Alternatively, a little more informative expression can be obtained as follows: letting $\gamma(h,s) \coloneqq \alpha_0(s) A(h, s)$,
where $\gamma^0(\beta_{0j}, s) = \beta_{0j}(s)$, and $\gamma^\ell(\beta_{0j}, s) = \gamma(\gamma^{\ell - 1}(\beta_{0j}, \cdot), s)$ for $\ell \geq 1$. From this, we can clearly see that the marginal effects $M(i,j,s)$ of increasing $X_{it}^j$ consist of the direct effect on unit $i$, the indirect effect on $i$'s immediate neighbors, the second-order indirect effect on $i$'s neighbors' neighbors, and so forth, highlighting the presence of the network multiplier effect. More specifically, recall that when $W_n$ represents a (weighted) adjacency matrix, the $(i,j)$-th element of $W_n^\ell$ corresponds to the number of (weighted) walks between $i$ and $j$ of length $\ell$. Thus, the $k$-th element of $M(i,j,s)$ is interpreted as the weighted sum of the number of walks from $i$ to $k$, where the contribution of each length-$\ell$ walk to the sum decays exponentially at $\gamma^\ell(\beta_{0j}, s)$.
To estimate the marginal effects, not just replacing the unknown parameters with their estimators, the infinite sum generally needs to be approximated by a truncated sum: for some positive integer $S$,
where $\widehat \gamma_{nT}(h,s) \coloneqq \widehat \alpha_{nT}(s) A(h,s)$. Meanwhile, in the special case of concurrent interaction such that $\gamma(h,s) = \alpha_0(s)h(s)$, $\gamma^2(h, s) = (\alpha_0(s))^2 h(s)$, \ldots, it is easy to see that $M(i,j,s) = \sum_{\ell = 0}^\infty (\alpha_0(s) W_n)^\ell \bm{e}_i \beta_{0j}(s) = (I_n - \alpha_0(s) W_n)^{-1}\bm{e}_i \beta_{0j}(s)$ holds. Thus, in this case, we can estimate $M(i,j,s)$ directly as $(I_n - \widehat \alpha_{nT}(s) W_n)^{-1}\bm{e}_i \widehat \beta_{nT,j}(s)$, without computing the infinite sum.
In the above discussion, we have demonstrated how the impacts of shifting one's covariate propagate to others. Similarly, just like the impulse response analysis in time-series vector autoregression, we can consider network impulse responses when an external shock occurs at a given unit. In particular, in a similar spirit to koop1996impulse, we define
where $\eta$ is a given "function" representing the external shock. By a similar calculation as above, we obtain
When plotting each element of $W_n^\ell \bm{e}_i \gamma^\ell(\eta, s)$ against $\ell = 0, 1, 2, \ldots$, it can be interpreted as a network version of the impulse response function (as a function of $\ell$), similarly to denbee2021network.
When assuming a concurrent interaction model, the impulse responses at $s$ take the following form: $I(i,\eta,s) = (I_n - \alpha_0(s) W_n)^{-1}\bm{e}_i \eta(s)$. Thus, if there is no exogenous shock at $s$, i.e., if $\eta(s) = 0$, the expected outcome at $s$ remains unaffected. This implies, for instance, that a travel demand shock that occurred five minutes ago has no impact on current mobility availability, which is unrealistic. On the other hand, if the interaction structure is given by $A(h,s) = \int_0^1 h(u) \nu(u,s) \text{d}u$ with $\nu(s', s) \neq 0$ for $s' < s$, then a shock occurring at $s'$ can transmit to the outcome at $s$, leading to nonzero impulse responses at $s$ even when $\eta(s) = 0$.
The estimation of $I(i,\eta,s)$ can be performed in the same manner as above. For some positive integer $S$, we estimate $I(i,\eta,s)$ by $\widehat I^S_{nT}(i,\eta,s) \coloneqq \sum_{\ell = 0}^S W_n^\ell \bm{e}_i \widehat \gamma_{nT}^\ell(\eta, s)$. The next proposition provides the convergence rate of $\widehat M^S_{nT}(i,j,s)$ and that of $\widehat I^S_{nT}(i,\eta,s)$.
This proposition indicates that the uniform convergence rates for the marginal effect and the impulse response estimators depend on the uniform convergence rate of the integrated-GMM estimator and the summation order $S$. \Copy{S}{Since the approximation error from truncating the infinite sum decreases geometrically as $S$ increases, if $\overline{\alpha}_0$ is close to zero, setting $S = 4$ or $5$ would be sufficient. Although $\overline{\alpha}_0$ is unknown in practice, it can be naturally estimated by $\overline{\widehat{\alpha}}_{nT} \coloneqq \max_{s \in [0,1]} |\widehat{\alpha}_{nT}(s)|$. Then, if the value of $\overline{\widehat{\alpha}}_{nT}$ is relatively large, it would be preferable to use a larger truncation order, say $S = 9$ or $10$.}
In this section, we conduct a series of Monte Carlo experiments to evaluate the finite-sample performance of the integrated-GMM estimator. To save space, we only summarize the main findings from the Monte Carlo simulations here; a complete description of the simulation design and full numerical results are provided in Appendix (ref). We compare the performance of three estimators: GMM 1: the integrated-GMM estimator using the weight matrix given in (ref), GMM 2: the integrated-GMM estimator using the identity weight matrix, and (integrated) 2SLS: GMM 1 estimator without utilizing the quadratic moment conditions.
Overall, all estimators perform well in terms of bias; however, GMM 1 consistently achieves the lowest RMSE (root-mean-squared error) across almost all settings, highlighting the efficiency gains from incorporating quadratic moment conditions. Interestingly, GMM 2 underperforms compared to 2SLS in many cases, indicating that the choice of the GMM weight matrix is equally (or potentially more) crucial to the inclusion of additional moment conditions. Meanwhile, as expected, 2SLS is more sensitive to the strength of IV than the others, with its RMSE for estimating $\alpha_0$ decreasing more drastically as the IV strength increases. For the estimation of $\beta_0$, the strength of IV has almost no impact, which is intuitive. Overall, the choice of $L$ seems to have only minor impacts on performance. Lastly, when we increase the sample size from $nT = 240$ to $nT = 960$, the RMSE values are roughly halved, demonstrating $\sqrt{nT}$-consistency of the estimators, which numerically corroborates our theoretical result in Theorem (ref)(ii).
In this section, we present several extensions of the current framework that are useful in empirical applications and can be implemented with relatively minor modifications to our method.
So far, we have assumed that the error terms are independent over time, which may be restrictive in practice. A practical approach to accommodating time effects is to extend the present model (ref) to the following two-way fixed effects model: for $s \in [0,1]$,
To facilitate the analysis of this model, we introduce the following additional conditions.
Assumption (ref)(i) rules out isolated units. In condition (ii), we assume that $P_{m,1}$ satisfies Assumption (ref)(i), while the restriction $\text{diag}(P_{m,1}) = \bm 0_n$ is now no longer necessary.
Under Assumption (ref)(i), to eliminate the time effects $c_{0t}(s)$ from the model, we can consider de-meaning by the neighborhood average rather than the conventional de-meaning by the global average. Specifically, letting $\mathcal{R}_n \coloneqq I_n - W_n$ and $\underset{N \times N}{\bm{\mathcal{R}}} \coloneqq I_{T-1} \otimes \mathcal{R}_n$, we obtain
It is straightforward to verify that the linear moment conditions
hold. In addition, by Assumption (ref)(ii), the quadratic moment condition
also holds. Based on these moment conditions, we can construct a GMM estimator for $\alpha_0$ and $\beta_0$ as in Subsection (ref). With an abuse of notation, let $\widehat \alpha_{nT}(s)$ and $\widehat \beta_{nT, j}(s)$ be the resulting estimators of $\alpha_0(s)$ and $\beta_{0j}(s)$, respectively.
The proof of Theorem (ref) is outlined in Appendix (ref). It is worth noting that constructing a weight matrix $P_{m,1}$ satisfying $\text{diag}(\mathcal{R}_n^\top P_{m,1} \mathcal{R}_n)=\bm 0_n$ is relatively straightforward. For an arbitrarily chosen $\widetilde P_{m,1} = (\widetilde p_{m,i,j})$, let $\mathcal{L}=\text{diag}(\ell_{1,1},\ldots,\ell_{n,n})$ be a diagonal matrix and set $P_{m,1} = \widetilde P_{m,1} + \mathcal{L}$. Then, for each $i \in [n]$,
where $r_{j,i}$ denotes the $(j,i)$-th element of $\mathcal{R}_n$. Hence, the constraint $\text{diag}(\mathcal{R}_n^\top P_{m,1} \mathcal{R}_n)=\bm 0_n$ is equivalent to the linear system
where $\bm y = - \text{diag}(\mathcal{R}_n^\top \widetilde P_{m,1} \mathcal{R}_n)$, $\bm C = (c_{i,j})$, $c_{i,j}= r_{j,i}^2$, and $\bm \ell = (\ell_{1,1},\ldots,\ell_{n,n})^\top$. Therefore, one may choose $\bm \ell = \bm C^{-} \bm y$, where $\bm C^{-}$ is a generalized inverse of $\bm C$.
As a straightforward extension of our current NAR model, one may consider a nonlinear NAR model:
where $\mu_{it}(s) \coloneqq X_{it}^\top \beta_0(s) + f_{0i}(s) + \varepsilon_{it}(s)$, and $m_s$ is a known transformation (up to parameters). For example, consider $m_s(Y_{jt}) = b_1(s)Y_{jt} + b_2(s) Y_{jt}^2 + \cdots + b_P(s) Y_{jt}^P$. Then, setting $\alpha_{0p}(s) \coloneqq \alpha_0(s) b_p(s)$, we can write the model as
where $\overline Y_{it}^{(p)} \coloneqq \sum_{j = 1}^n w_{i,j} Y_{jt}^p$.
Another potentially useful extension is a multiple-network model:
where $W_n^{(p)} = (w^{(p)}_{i,j})$, $p \in [P]$, are different interaction matrices. In this case, if we define $\overline Y_{it}^{(p)} \coloneqq \sum_{j = 1}^n w^{(p)}_{i,j} Y_{jt}$, the model can again be expressed as (ref).
For both cases, as long as valid IVs $Q_{it}$ for the multiple endogenous components $A(\overline Y_{it}^{(p)}, s)$, $p \in [P]$, are available, essentially the same estimation and inference procedure as in Section (ref) can be applied, while we need to impose more restrictive conditions to ensure completeness. Additionally, models with multiple kernels $\nu^{(p)}$, $p \in [P]$, or heterogeneous interaction effects of the form $\alpha_0(s) = \sum_{p = 1}^P \text{x}_{it}^{(p)} \alpha_{0p}(s)$ can also be analyzed in the same way.
The integrated-GMM estimator is often infeasible because the response functions are typically observed only at a finite set of points in $[0,1]$. Even in such cases, we can approximate the entire functional form of $Y_{it}$ using a linear interpolation method.
Suppose that for each $(i,t)$, $Y_{it}$ is observed at $L_{it}$ distinct points $0 \le s_{it,1} \le s_{it,2} \le \dots \le s_{it,L_{it}} \le 1$. Then, for each given $s \in [s_{it,l}, s_{it,l + 1}]$, define
When $s < s_{it,1}$ (resp. $s > s_{it,L_{it}}$), we set $Y^{\text{int}}_{it}(s) \coloneqq Y_{it}(s_{it,1})$ (resp. $Y^{\text{int}}_{it}(s) \coloneqq Y_{it}(s_{it,L_{it}})$).
Other than linear interpolation, one may also use a kernel method, as in zhu2022network, to obtain $Y^{\text{int}}_{it}(s)$. Then, using $Y^{\text{int}}_{it}(s)$ in place of $Y_{it}(s)$, we can write
where $u_{it}(s)$ is the interpolation error: $u_{it}(s) \coloneqq Y^{\text{int}}_{it}(s) - Y_{it}(s) + \alpha_0(s) A(\overline Y_{it} - \overline Y^{\text{int}}_{it}, s)$. Thus, if $u_{it}(s)$ converges to zero sufficiently quickly for all $s \in [0,1]$, $i \in [n]$, and $t \in [T]$, we can apply the same estimation and inference strategy as above.
As an empirical application, we analyze spatial interactions in the demand for a bike-sharing system in the U.S. Demand analysis of shared mobility has been a highly active research topic in recent years across various areas, including transportation research, marketing, economics, and environmental studies. In particular, bike-sharing systems have attracted increasing attention. For a comprehensive review of this literature, see eren2020review.
The dataset used in this analysis comes from the Bay Area Bike Share in San Francisco, which was established in August 2013 and is now known as Bay Wheels. The dataset is publicly available on the Kaggle website.\footnote{\url{https://www.kaggle.com/datasets/benhamner/sf-bay-area-bike-share}} It contains detailed information about the system from August 2013 to August 2015, including station locations, the number of available bicycles at each station over time, and all trip-level data during this time period. The trip data include details such as start and end times and stations, as well as the user type (subscriber or casual user). In this dataset, there are 70 bike stations in total; for a map of all 70 station locations, see Figure (ref) in Appendix (ref).
Since the initial installation of stations in August 2013, the 70th station (Ryland Park station) was added in April 2014. Accordingly, we use data from May 2014 to August 2015 for this analysis, which represents the largest balanced panel dataset that can be extracted from the raw data.
\Copy{rebalance}{ One concern in the analysis is that shared mobility services often relocate bikes across stations in order to rebalance bike availability. To mitigate the effects of rebalancing operations in demand analysis, several approaches have been considered in the literature. For example, using the number of trips, pickups, or returns at each station as the outcome variable is a common strategy (e.g., el2017effects, sun2018promoting). However, this approach is not straightforwardly applicable in the present functional data framework. Another approach is to detect potential rebalancing operations using heuristic thresholding rules (e.g., faghih2014land, faghih2017empirical). Although it is generally impossible to access formal records of relocation operations, as noted by gammelli2022predictive, the majority of bike-sharing systems are rebalanced during nighttime (this is called static rebalancing). Then, if we observe sudden jumps or drops in bike counts during nighttime hours, they are likely due to relocation operations conducted by the service provider. We adopt this approach in our analysis. Specifically, we first identify instances where the number of available bicycles jumps up/down by more than or equal to 10 all at once. We then examine the distribution of such events across hours and days, as shown in Figure (ref) in Appendix (ref). This figure indicates that sudden drops or increases in bike availability tend to occur between midnight and early morning, particularly on Sundays. Based on this, we exclude these time periods from our analysis. }
Another concern is the enormous size of the dataset. Because the original data are recorded in minutes every day, using the raw data directly can lead to a memory problem. Moreover, daily data tend to fluctuate and to be noisy due to random events.
To address the aforementioned issues, we first rounded the trip data to 15-minute intervals and then averaged over Monday through Friday at each interval, discarding data from Saturdays and Sundays. Furthermore, to avoid potential bike relocation events in weekdays, we restrict the analysis to the time period from 6 AM to 9 PM. Consequently, our final dataset is a weekly-level panel with $n = 70$ stations and $T = 69$ weeks. The outcome of interest is the number of bicycles at each station for $s \in [0,1]$, where $s=0$ corresponds to 6 AM and $s=1$ corresponds to 9 PM.
Figure (ref) presents the trajectories of average bike availability for all 70 stations during the first week in our panel. It clearly shows that most of the variation in bike availability occurs between 6 AM and 9 PM.
Based on the dataset constructed as described above, we estimate model (ref), where
Here, $\text{dist}(i,j)$ denotes the Euclidean distance between stations $i$ and $j$. \Copy{KL}{ The estimation procedure is the same as the GMM 1 estimator in Section (ref). The number of basis terms is selected via the cross-validation procedure described in Appendix (ref), using the last 11 weeks as the validation sample. As a result, we set $K = 8$ and $L = 3K$. } The rainy day dummy and month dummy variables are not used as IVs. All integrals are approximated by finite summations over grid points at 15-minute intervals.
The estimation result for the interaction effect function $\alpha_0$ is presented in Figure (ref). In the figure, the shaded area depicts the (pointwise) 95% confidence interval. From the figure, we observe that positive spatial interaction in bicycle availability exists during the morning hours. It is plausible that as bike-sharing becomes more popular particularly among commuters, it encourages further use of the service, thereby reinforcing demand during the morning. Meanwhile, interestingly, negative interaction appears around 5--7 PM. In the evening, main users may include not only returning commuters but also individuals going out for dining, shopping, concerts, etc. As a result, bicycles might accumulate at certain popular stations while nearby less-popular stations experience lower availability, leading to the negative interaction. \Copy{negative}{ Alternatively, as kindly suggested by a referee, the negative correlation may also be driven by rebalancing operations conducted by the service provider, which artificially accumulate bikes at popular origin stations. However, it should be noted that autoregressive models of the current type are generally able to identify the presence of autocorrelation but are not designed to distinguish its underlying mechanisms. Investigating the detailed mechanisms generating such autocorrelation would require a different empirical framework. }
To save space, the estimation results for $\beta_0(s)$ are presented in Figure (ref) in Appendix (ref), excluding the coefficients for the month dummies. Among the key covariates, we observe that only the ratio of arriving subscribers has a statistically significant positive impact on bike availability. This result is intuitive, as stations with a higher number of regular users arriving are expected to hold a richer stock of bikes. For other variables, for instance, the rainy day dummy has a positive effect on bike availability, which is consistent with previous studies, though the effect is not statistically significant. One possible explanation is that the rainy "day" dummy does not capture detailed temporal variations (i.e., it is not a function of $s$), and since our dataset is averaged over weeks, these may have diluted its impact.
Lastly, we conduct an impulse response analysis. The figures summarizing the results are presented in Figure (ref). For illustration, we arbitrarily select the Embarcadero at Folsom station as the target station receiving an external shock. Specifically, we consider a hypothetical scenario in which the bike stock at this station is reduced by 2 at the peak of 9 AM (panel (a)). Panels (b) and (c) illustrate how the shock propagates to its two nearest stations, Spear at Folsom and Temporary Transbay Terminal. These figures indicate that the external shock spills over to these stations with a slight time delay, peaking just before 10 AM. Since the magnitude of both the external shock and spatial interaction is moderate in this analysis, the impulse responses for both stations are relatively mild.
\Copy{rbcheck}{ For robustness checks, we also estimated the model using $K = 7$ and $K = 9$. The corresponding estimates of $\alpha_0$ are presented in Figure (ref) in Appendix (ref). In these figures, we can observe that the estimated interaction function exhibits the same pattern, with positive interaction in the morning hours and negative interaction at night. }
\Copy{nonfunctional}{ Finally, for comparison, we also estimate a non-functional panel spatial autoregressive model with individual fixed effects in lee2010estimation. The dependent variable is the average weekly bike availability at each station, obtained by averaging $Y_{it}(s)$ over $s \in [\text{6 AM}, \text{9 PM}]$. To save space, the result is reported in Table (ref) in Appendix (ref). Notably, the estimated spatial autoregressive parameter is slightly negative and statistically significant, possibly capturing the negative interaction during nighttime, as observed in Figure (ref). However, as the figure also shows, positive interaction occurs during the morning hours. This illustrates that ignoring the functional nature of time-varying bike availability may result in a misinterpretation of the magnitude and direction of spatial interactions that change dynamically over time. }
In this paper, we proposed a novel functional regression framework to analyze spatial and network interactions in functional panel data settings. By extending the standard NAR model to accommodate functional outcomes and individual fixed effects, we developed an integrated-GMM estimator that can estimate the functional parameters potentially more efficiently than 2SLS-based estimators. Under certain conditions, we established the theoretical properties of our estimator, including consistency, convergence rates, and asymptotic normality, and confirmed its finite-sample performance through Monte Carlo simulations. As an empirical application, we analyzed the demand for a bike-sharing system in the San Francisco Bay Area and discovered significant spatial interactions in bike availability that vary over the time of day. This complex interaction pattern cannot be observed if the functional nature of the response variable is ignored, and thus our finding highlights the importance of accounting for functional spatial dependencies in modeling the demand for shared mobility services.
\Copy{futurework}{ There are various avenues to extend our methodological and theoretical contributions. For example, following beyaztas2024spatial and hoshino2024functional, one may instead consider the specification $\int_0^1 Y_{jt}(u)\alpha_0(u,s)\text{d}u$ as an alternative to $\alpha_0(s)\int_0^1 Y_{jt}(u)\nu(u,s)\text{d}u$, which corresponds to the standard functional linear form in the FDA literature. \Copy{illposed}{ When $m_s$ in (ref) is unknown and estimated nonparametrically, this formulation leads to an ill-posed estimation problem (e.g., newey2013nonparametric,hoshino2022sieve,yang2025estimation), which substantially complicates the asymptotic analysis.} Also, one might want to consider a dynamic interaction structure in which the outcomes in the previous time period affect the current outcome through $A(Y_{j,t-1}, s)$. Such a model can be viewed as a functional extension of zhu2017network. The choice between our simultaneous model and the dynamic alternative should be decided by the data type. For instance, in the analysis of local interactions of mobility usage, where demand propagation occurs almost instantaneously, the simultaneous model may be more appropriate. In summary, many different types of functional interaction structures, as well as their combinations, can be considered. We leave these extensions to future research. }
We would like to thank co-editor, an associate editor, and four anonymous referees for their insightful comments that greatly improved our paper. The authors thank the participants at the seminars of National Cheng Kung University and National Taiwan University for their valuable feedback. Hoshino's work was supported by JSPS KAKENHI Grant Number 23KK0226. Most parts of this paper were written during Hoshino's research visit at the Melbourne Business School (MBS), University of Melbourne. He is deeply grateful to MBS for their hospitality.