Functional instrumental variable regression with an application to estimating the impact of immigration on native wages

Functional instrumental variable regression with an application to estimating the impact of immigration on native wages

Introduction

The recent developments in data collection and storage technologies ignite studies on how to use more complicated observations such as curves, probability density functions, or images. This area of study, commonly called functional data analysis, has become popular in statistics, and researchers in various fields, including economics, have been benefited from advances in this area. In particular, for practitioners who are interested in studying the relationship between two or more such variables, functional linear models are of central importance, and crucial contributions on this topic include Bosq2000, Yao2005, Mas2007, Hall2007, Park2012, crambes2013asymptotics, Benatia2017, imaizumi2018, Sun2018, and Chen_et_al_2020 to name only a few.

The existing statistical approaches for estimating the functional linear model, including those proposed in the aforementioned literature, are mostly established under the assumption that the explanatory variable of interest is exogenous, meaning that it is uncorrelated with the regression error. However, this assumption is not likely to hold in practice; that is, explanatory variables are often {\it{endogenous}}. The issue of endogeneity is particularly relevant in the context of functional linear models because functional observations used in the analysis are typically constructed by smoothing their discrete, and often sparsely observed, realizations (see, e.g., Yao2005). If this being the case, the functional observations may inevitably involve small or large measurement errors, which leads to the violation of the exogeneity condition at least to some degree (see, e.g., Sections (ref) and (ref)). This issue may hinder practitioners from applying the functional linear model.

If the volume of the literature on functional linear regression under the exogeneity condition provides any indication, the developments made so far to deal with endogeneity do not seem to sufficiently meet practitioners' needs. Although a few papers, such as Benatia2017 and Chen_et_al_2020, study the issue of endogeneity in the functional linear model, not much is known about asymptotic distributions of their estimators and how to implement statistical inference on the parameter of interest; this may limit the practical applicability of the functional endogenous linear model. We will fill this gap to some extent by providing new estimators and inferential methods based on their asymptotic properties. This is a crucial point where the present paper is differentiated from the existing ones concerning the issue of endogeneity in the functional linear model.

Specifically, this paper provides new estimation results for the functional endogenous linear model based on (i) the functional principal component analysis (FPCA) and (ii) the instrumental variable (IV) approach. The former has been widely adopted by researchers dealing with functional data (see e.g., Ramsay2005; shang2014survey), and the latter has also been widely adopted in order to address endogeneity not only in the conventional Euclidean space setting (see e.g., Bekker1994; Chao2005; Newey2009a), but also in the setting involving functional observations (see, e.g., Carrasco2012; Florence2015; Benatia2017; Chen_et_al_2020; Babii2021). However, the application of the FPCA to the functional endogenous linear model has not been fully explored.

We consider the case where the response variable $y_t$, explanatory variable $x_t$, and instrumental variable $z_t$ are all function-valued; of course, with slight modifications, our results to be subsequently given can be adjusted for the case where $y_t$ is scalar- or vector-valued. Unlike in the most of the papers mentioned above, we do not require the variables of interest to be independent and identically distributed (iid), but allow those to exhibit some weak dependence so that our methodology can be applied to various empirical examples. Given that many functional observations considered in the literature for applications in fields of energy, environmental and financial economics tend to involve time dependence, this extension may be attractive to practitioners.

Among the aforementioned papers, the study by Chen_et_al_2020 is most closely related to the present paper in the sense that they consider FPCA-based consistent estimation of function-on-function regression models with endogeneity introduced by measurement errors. Benatia2017 earlier considered a similar model and proposed a consistent estimation method, but their theoretical results are obtained from a quite different theoretical methodology (ridge-type regularization). We complement these studies by providing new FPCA-based estimators and in-depth discussion on their asymptotic properties.

Technically, we view the function-valued variables of interest as random variables taking values in a Hilbert space of square-integrable functions, and then propose our FPCA-based functional IV estimator (FIVE). As is well known in the literature, estimation of a model involving function-valued random variables is not straightforward because some important sample operators, such as the covariance of such a random variable, are not invertible over the entire Hilbert space(s). We circumvent this issue by employing a rank-regularized inverse of such an operator, and this is the point where we make use of the FPCA. The reason why we focus on this regularzation scheme comes not from its theoretical superiority, but merely from its popularity in the literature. Other schemes such as ridge-type regularization (e.g., Florence2015; Benatia2017) may be alternatively adopted, and are expected to have their own merits (see Remark (ref)). \phantomsection{\leavevmode\color{black}It is worth summarizing some crucial differences between our estimators and the alternative estimator proposed by Benatia2017 based on ridge regularization. To the best of our knowledge, there has not been exploration of asymptotic inference on specific characteristics of the regression (coefficient) operator for their estimator; this seems to be because of a nontrivial challenge associated with a particular asymptotic bias (see Section 4 of Benatia2017). In contrast, in this paper, we tackle the issue by employing the FPCA augmented with a proper extension of the asymptotic approach introduced by Hall2007. As a result, this paper provides mathematical conditions that support a valid asymptotic inference on the regression operator. Moreover, \citepos{Benatia2017} methodology does not take into account for the more general presence of weak dependence, although we believe their results can be extended to such a setting.}

This paper studies in depth the asymptotic properties of the proposed estimators. It is first shown that, under some mild conditions on the data generating process of $\{y_t,x_t,z_t\}_{t\geq1}$, the FIVE achieves the weak (convergence in probability) and strong (almost sure convergence) consistencies as long as the regularization parameter, which is introduced for a rank-regularized inverse of a certain sample operator used to construct the estimator, decays to zero at an appropriate rate. We then establish more detailed asymptotic properties of the FIVE under some nonrestrictive assumptions on the eigenstructure of the cross-covariance operator of the explanatory variable $x_t$ and the IV $z_t$. By doing so, we can see how the cross-covariance structure of $x_t$ and $z_t$ and the choice of the regularization parameter affect the convergence of the FIVE toward its true counterpart. In addition to these results, we show that the FIVE is asymptotically normal in a pointwise sense if it is centered at a certain operator that is slightly biased from the true parameter of interest; moreover, if certain additional conditions are satisfied, such a bias becomes asymptotically negligible and thus, in this case, the FIVE centered at the true parameter becomes asymptotically normal. The asymptotic normality results given in this paper are quite different from similar results given in a finite dimensional setting in the sense that the convergence rate is (i) possibly random and (ii) not uniformly given over the entire Hilbert space on which our estimator is defined. This result implies that the proposed estimator does not weakly converge to any elements in the usual operator topology, which {\leavevmode\color{black}generalizes} what Mas2007 earlier found in the context of functional autoregressive (AR) models of order 1. Based on our study of the FIVE, we also propose a different but closely related estimator, called the functional two-stage least square estimator (F2SLSE) and obtain its asymptotic properties in a similar manner. We discuss how our estimators and their asymptotic properties can be used to implement usual statistical inference on the parameter of interest.

To see how the asymptotic properties of our estimators are revealed in finite samples, we implement Monte Carlo experiments under various simulation designs. The simulation results are quite satisfactory. Overall, it seems that our estimators can be good alternatives or sometimes complements to some existing estimators that are closely related to ours.

As an empirical illustration, we study the impact of immigration on native wages. Specifically, we employ a model that is similar to those considered by, e.g., Dustmann2012 and Sharpe2020101902. The previous literature in this area, including Ottaviano2011, Card2009, and the aforementioned articles, show that an inflow of immigrants differently affects native wages depending on skill levels (captured by, e.g., years of education and experience) of both natives and immigrants. We, in this paper, investigate such heterogeneous effects using our functional linear model, which is initiated by viewing both the labor supply and the native wage as functions of a certain measure of workers' skill (will be detailed in Section (ref)). This approach has a couple of advantages compared to that taken in the earlier literature. For example, in the previous literature, workers of various levels of skill are often classified into a few skill groups before analysis, which is necessitated to reduce dimensionality of the considered model (see Example (ref) and Section (ref)). However, such a pre-classification, which may affect estimation results and their interpretation, is not required in our approach. Moreover, our methodology allows for studying if an inflow of immigrants in a particular skill group heterogenously affects workers equipped with different skill levels. Using the methodology developed in this paper, we find evidence supporting the presence of heterogeneous effects of immigration.

This paper is organized as follows. Section (ref) introduces a functional endogenous linear model and provides motivating examples. In Section (ref), we define the FIVE and discuss its asymptotic properties. Section (ref) introduces the F2SLSE and discusses its asymptotic properties. Section (ref) reports simulation results and details our empirical example. Section (ref) concludes. The mathematical proofs of the theoretical results can be found in the Supplementary Material.

Functional endogenous linear model

Endogeneity and motivating examples

We suppose that a stationary sequence of random functions $\{y_t, x_t,u_t\}_{t\geq 1}$ satisfies the following:

where $c_y$ is the intercept function and $\mathcal A$ is a linear operator satisfying certain conditions to be clarified. In \hyperref[{eqmodel0}]{\tagform@{\ref*{eqmodel0}}}, $y_t$, $x_t$ and $u_t$ will be technically understood as random variables taking values in separable Hilbert spaces. Appendix (ref) briefly introduces the definitions of a Hilbert-valued random variable $X$, its expectation (denoted $\mathbb{E}[X]$), covariance operator (denoted $\mathcal C_{XX}\coloneqq\mathbb{E}[(X-\mathbb{E}[X]) \otimes (X-\mathbb{E}[X])]$), and cross-covariance operator with another Hilbert-valued random variable $Y$ (denoted $\mathcal C_{XY} \coloneqq \mathbb{E}[(X-\mathbb{E}[X]) \otimes (Y-\mathbb{E}[Y])]$), where $\otimes$ signifies the tensor product defined {\leavevmode\color{black}by $X\otimes Y(\cdot) = \langle X, \cdot \rangle Y $ for any random or nonrandom $X$ and $Y$ taking values in $\mathcal H$ (see \hyperref[{eqtensor}]{\tagform@{\ref*{eqtensor}}})}.

We say that the explanatory variable $x_t$ is endogenous if the cross-covariance of $x_t$ and $u_t$, given by the operator $\mathbb{E}[(x_t-\mathbb{E}[x_t])\otimes (u_t-\mathbb{E}[u_t])]$, is nonzero. The present paper focuses on estimation and inference of the functional linear model in the presence of endogeneity. Below we provide specific and practical examples that motivate this model of interest.

It is expected from Example (ref) that endogeneity can arise in many practical applications of the functional linear model, where $x_t$ is incompletely observed. In such a case, the exogeneity condition is likely violated. \phantomsection{\leavevmode\color{black}In particular, due to advancements in data collection techniques, it is now possible to construct density- or curve-valued economic variables from large datasets. As a result, the analysis of these variables has gained popularity, as evidenced by various empirical examples in the literature (e.g., Benatia2017; Babii2021; Chen_et_al_2020; Nielsen2019; seo2020functional). In economic functional data, observations are often incomplete, and the discrete and finite realizations used to construct functional observations may not be enough to fully capture the entire function. Therefore, practitioners may remain cautious about potential endogeneity, even if the considered regressor $x_t$ is presumed to be exogenous. This limitation could hinder the practical use of the functional linear model.} As well expected from the literature on the standard linear simultaneous equation model, endogeneity should be properly addressed for consistent estimation of the regression operator. A widely used strategy to do this is the IV approach, which will be pursued in our Hilbert space setting.

Model, assumptions, and notation

To facilitate the subsequent discussions, it may be helpful to introduce some additional notation. We first let $\mathcal H$ denote the Hilbert space of square-integrable functions defined on the unit interval $[0,1]$, where the inner product $\langle \cdot, \cdot \rangle$ is defined by $\langle \zeta_1,\zeta_2 \rangle = \int_{0}^1 \zeta_1(s) \zeta_2(s) ds$ for $\zeta_1,\zeta_2 \in \mathcal H$ and $\|\cdot\|= \langle \cdot,\cdot\rangle^{1/2}$ defines the norm of $\mathcal H$. $\mathcal{L}_{\mathcal H}$ denotes the space of bounded linear operators acting on $\mathcal H$, equipped with the operator norm $\|\mathcal T\|_{\mathrm{op}} = \sup_{ \|\zeta\|\leq 1} \|\mathcal T\zeta\|$. For any $\mathcal T \in \mathcal L_{\mathcal H}$, we let $\mathcal T^\ast$, $\operatorname{ran} \mathcal T$, $\ker \mathcal T$, and $\|\mathcal T\|_{\operatorname{HS}}$ denote the adjoint, range, kernel, and the Hilbert-Schmidt norm of $\mathcal T$, which are briefly reviewed in Appendix (ref); in that section, various properties of $\mathcal T \in \mathcal L_{\mathcal H}$ such as nonnegativity, positivity, self-adjointness, compactness and Hilbert-Schmidtness are also reviewed. For any nonnegative, self-adjoint and compact $\mathcal T$, we may write $\mathcal T = \sum_{j=1}^\infty a_j\zeta_{j} \otimes \zeta_{j}$ for some nonnegative sequence $\{a_j\}_{j \geq 1}$ and some orthonormal basis $\{\zeta_{j}\}_{j \geq 1}$. Then $\mathcal T^{1/2}$ can be well defined by replacing $a_j$ with $\sqrt{a_j}$.

This paper concerns the case where the response variable $y_t$ and the endogenous explanatory variable $x_t$ are infinite-dimensional random variables taking values in separable Hilbert spaces. We hereafter conveniently assume that all of such variables take values in $\mathcal H$. This setup in fact encompasses a more general scenario where $y_t$ and $x_t$ take values in different separable Hilbert spaces of infinite dimension, say $\mathcal H_y$ and $\mathcal H_x$. This is because that these spaces are all isomorphic to $\mathcal H$ (see e.g., Corollary 5.5 in Conway1994, Conway1994), and thus there is no loss of generality by assuming that $\mathcal H_y=\mathcal H_x =\mathcal H$. We further assume for convenience that $y_t$ and $x_t$ have zero means, i.e., $\mathbb{E}[y_t]=\mathbb{E}[x_t]=0$; this assumption naturally makes $c_y$ in \hyperref[{eqmodel0}]{\tagform@{\ref*{eqmodel0}}} be suppressed to zero. The extension to the case where the means are unknown and needed to be estimated is straightforward. After adopting all such simplifying assumptions, the functional endogenous linear model, which will be subsequently considered, is given as follows: for a linear operator $\mathcal A:\mathcal H\mapsto \mathcal H$,

\phantomsection{\leavevmode\color{black}We then let $z_t$ (to be called the IV) be another zero-mean $\mathcal H$-valued random variable satisfying $\mathbb{E}[z_t\otimes u_t]=0$.} For notational convenience, we use $\mathcal C_{zz}$, $\mathcal C_{xz}$, $\mathcal C_{yz}$ and $\mathcal C_{uz}$ to denote the following operators,

Similarly, let $\widehat{\mathcal C}_{zz}$, $\widehat{\mathcal C}_{xz}$, $\widehat{\mathcal C}_{yz}$ and $\widehat{\mathcal C}_{uz}$ denote their sample counterparts that are computed as follows:

We will employ the following assumptions throughout the paper: below, $\mathfrak F_t$ denotes the filtration given by $\mathfrak F_t = \sigma\left( \{z_s\}_{s\leq t+1}, \{u_s\}_{s\leq t} \right)$, and $\widehat{\mathcal C}_{uu} = T^{-1} \sum_{t=1}^T u_t \otimes u_t$.

By Assumption (ref).(ref), we allow $\{x_t,z_t\}_{t\geq 1}$ to be a weakly dependent sequence; this is because (i) we want to accommodate various empirical examples such as those given in HK2012 by not restricting our attention to the iid case and (ii) the variables to be considered in our empirical application (Section (ref)) naturally exhibit time series dependence. In Assumptions (ref).(ref) and (ref).(ref), the error term $u_t$ is assumed to be a homoscedastic martingale difference sequence. Assumption (ref).(ref) states the conditions required for the IV $z_t$ in this setting (see Example (ref) for a possible IV presented for the model in Example (ref)), and this condition implies that $u_t$ is uncorrelated with $z_t$. In Assumption (ref).(ref), we impose some requirements on the moments of $u_t$. We here note that, if $\{z_t,u_t\}_{t\geq 1}$ is an iid sequence, as often assumed in the literature, Assumptions (ref).(ref) and (ref).(ref) reduce to the following: $$\mathbb{E}[u_t|z_t] = 0, \quad \mathbb{E}[u_t\otimes u_t|z_{t}] = \mathcal C_{uu},\quad \text{and}\quad \mathbb{E}[\|u_t\|^{2+\delta}|z_t]<\infty \text{ for some } \delta>0.$$ The Hilbert-Schmidt condition of $\mathcal A$ given in Assumption (ref).(ref) would become redundant if we considered a finite dimensional Hilbert space, but in our setting it imposes a nontrivial mathematical condition on $\mathcal A$. In Assumptions (ref).(ref) and (ref).(ref), high-level conditions on limiting behaviors of some sample operators are given, and these are for mathematical convenience. We first note that $\{x_t\otimes z_t-\mathcal C_{xz}\}_{t\geq 1}$, $\{ z_t \otimes z_t-\mathcal C_{zz} \}_{t\geq1}$, and $\{u_t \otimes z_t -\mathcal C_{uz}\}_{t\geq1}$ are sequences in the Hilbert space of Hilbert-Schmidt operators, denoted by $\mathcal S_{\mathcal H}$ (see Section (ref) of the Supplementary Material). If those sequences are iid (resp.\ geometrically strongly mixing), then Assumption (ref).(ref) holds once $\mathbb E [(\Vert x_t \Vert \Vert z_t\Vert)^{\upsilon}]$, $\mathbb E [\Vert z_t \Vert ^{2\upsilon}]$, and $\mathbb E [(\Vert u_t \Vert \Vert z_t \Vert)^{\upsilon}]$ are finite for some $\upsilon \geq 2$ (resp.\ $\upsilon \geq 2 + \delta$ for some $\delta>0$) Bosq2000; such primitive sufficient conditions can also be found for martingale differences Bosq2000 and weakly stationary sequences Bosq2000. We also observe that $\{u_t \otimes u_t-\mathcal C_{uu}\}_{t\geq1}$ is a martingale difference sequence in $\mathcal S_{\mathcal H}$, and some primitive sufficient conditions for Assumption (ref).(ref) can be found in e.g.,\ Theorems 2.11 and 2.14 of Bosq2000. {\leavevmode\color{black}Lastly, Assumption (ref).(ref) enables us to identify the unique bounded linear operator $\mathcal A$ satisfying \hyperref[{eqmodel}]{\tagform@{\ref*{eqmodel}}} using the IV $z_t$, which will be discussed in more detail in Section (ref)}.

Functional IV estimator

This section discusses estimation of the model \hyperref[{eqmodel}]{\tagform@{\ref*{eqmodel}}} given observations $\{y_t,x_t,z_t\}_{t=1}^T$. We first propose the FIVE in detail and study its asymptotic properties.

The proposed estimator

We find from \hyperref[{eqmodel}]{\tagform@{\ref*{eqmodel}}} that $\mathcal C_{yz}^\ast = \mathbb E[z_t \otimes y_t] = \mathcal A \mathbb E[z_t \otimes x_t] = \mathcal A \mathcal C_{xz}^\ast$ and hence,

As discussed in Mas2007 and Benatia2017, $\mathcal A$ is a uniquely identified bounded linear operator if and only if $ \ker \mathcal C_{xz} = \{0\}$ (see Assumption (ref).(ref)), and we note that all the eigenvalues of $\mathcal C_{xz}^\ast \mathcal C_{xz}$ are positive under the condition Mas2007. In the sequel, we thus let $\{\lambda_j^2\}_{j\geq 1}$ denote the collection of the eigenvalues of $\mathcal C_{xz}^\ast \mathcal C_{xz}$ ordered from the largest to the smallest, and represent $\mathcal C_{xz}^\ast \mathcal C_{xz}$ as its spectral decomposition given by

where $f_j$ is the eigenfunction corresponding to $\lambda_j^2$. Given \hyperref[{popmoment1}]{\tagform@{\ref*{popmoment1}}}, it may be natural to consider an estimator $\bar{\mathcal A}$ that satisfies the equation $\widehat{\mathcal C}_{yz}^\ast\widehat{\mathcal C}_{xz}=\bar{\mathcal A} \widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz}$, obtained by replacing $\mathcal C_{yz}$ and $\mathcal C_{xz}$ with their sample counterparts. However, it is generally impossible to directly compute the estimator $\bar{\mathcal A}$ from this equation since $\widehat {\mathcal C}_{xz} ^\ast \widehat {\mathcal C}_{xz}$ is not invertible over the entire Hilbert space $\mathcal H$. We circumvent this issue by employing a regularized inverse of $\widehat {\mathcal C}_{xz} ^\ast \widehat {\mathcal C}_{xz}$ which may be understood as the well defined inverse on a strict subspace of $\mathcal H$.

To this end, we first note that $\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz}$ is nonnegative, self-adjoint and compact and hence allows the following representation:

where $\{ \widehat{\lambda}_j^2 , \widehat{f}_j \}_{j\geq 1}$ are the pairs of eigenvalues and eigenfunctions, and $\widehat{\lambda}_1^2 \geq \ldots \geq \widehat{\lambda}_T^2 \geq 0 = \widehat{\lambda}_{T+1}^2=\ldots$. We then define $\operatorname{{\mathrm{K}}}$ as the random integer determined by the threshold parameter $\upalpha >0$ such that

Using the first $\operatorname{{\mathrm{K}}}$ eigenfunctions of $\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz}$, its rank-regularized inverse, denoted $(\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz})_{\operatorname{{\mathrm{K}}}} ^{-1}$, and the FIVE, denoted $\widehat { \mathcal A}$, are defined as follows:

The largest eigenvalue of the regularized inverse $(\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz})_{\operatorname{{\mathrm{K}}}}^{-1}$ is bounded above by $\upalpha$ and thus the regularized inverse is a well defined bounded linear operator for every $\upalpha>0$. It is worth mentioning that the FIVE becomes equivalent to the estimator proposed by Park2012 in the case where $z_t = x_t$ and $\operatorname{{\mathrm{K}}}$ is deterministically chosen by researchers (see Remark (ref) below), so our estimator may be understood as an extension of their estimator. We also note that the FIVE $\widehat{\mathcal A}$ may be viewed as a sample-analogue of $\mathcal A$ satisfying \hyperref[{popmoment1}]{\tagform@{\ref*{popmoment1}}} in the sense that $\widehat{\mathcal A}$ is the solution to $\widehat{\mathcal C}_{yz}^\ast\widehat{\mathcal C}_{xz}=\widehat{\mathcal A} \widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz}$ on the restricted domain given by $\widehat{\mathcal H}_{\operatorname{{\mathrm{K}}}} = \operatorname{span}\{\widehat{f_j}\}_{j=1}^{\operatorname{{\mathrm{K}}}}$. {Section (ref) of the Supplementary Material discusses on how $\widehat{\mathcal A}$ can be computed from the data using the FPCA.}

General asymptotic properties

As may be deduced from our construction of the FIVE, $\widehat{\mathcal A} = 0$ is imposed outside a subspace whose dimension increases as $\upalpha$ gets larger. Thus, for $\widehat{\mathcal A}$ to be a consistent estimator of $\mathcal A$ defined on the entire space $\mathcal H$, the {regularization} parameter $\upalpha$ given in \hyperref[{eqdef2}]{\tagform@{\ref*{eqdef2}}} needs to diverge to infinity. Taking this into consideration, we investigate the asymptotic properties of the FIVE when $T \to \infty$ and $\upalpha\to \infty$ jointly. We will employ the following assumption:

That is, the eigenvalues of $\mathcal C_{xz}^\ast \mathcal C_{xz}$ are required to be distinct. This is employed to see asymptotic properties of the FIVE in detail and does not seem to be restrictive in practice; in fact, similar assumptions have been employed in the literature on functional linear models, see e.g., Bosq2000, Mas2007, Hall2007 and Park2012 to name only a few.

We now provide the asymptotic properties of the estimator $\widehat{\mathcal A}$ when $\upalpha$ and $T$ grow jointly without bound. To this end, we consider the following decomposition of $\widehat {\mathcal A} - \mathcal A$:

where $\widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$ denotes the orthogonal projection defined by $ \widehat{\Pi}_{\operatorname{{\mathrm{K}}}} = \sum_{j=1}^{\operatorname{{\mathrm{K}}}} \widehat{f}_j \otimes \widehat{f}_j$ and $\mathcal I$ is the identity operator acting on $\mathcal H$. Given that the FIVE is computed on the restricted domain $\operatorname{ran} \widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$ (note that $\widehat{\mathcal A}=0$ on $\operatorname{ran} (\mathcal I-\widehat{\Pi}_{\operatorname{{\mathrm{K}}}})$ by construction) the first term of \hyperref[{eqdecom}]{\tagform@{\ref*{eqdecom}}} may be understood as the deviation of $\widehat{\mathcal A}$ from $\mathcal A$ on $\operatorname{ran} \widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$. Thus, this term is hereafter called the deviation component on the restricted domain (the DR component). On the other hand, the second term $\mathcal A (\mathcal I-\widehat{\Pi}_{\operatorname{{\mathrm{K}}}})$ may be understood as the bias induced by the fact that $\widehat{\mathcal A}$ is enforced to zero on $\operatorname{ran} (\mathcal I-\widehat{\Pi}_{\operatorname{{\mathrm{K}}}})$. We thus call this term the regularization bias component (the RB component). Our first result below shows that both the DR and RB components are asymptotically negligible and thus $\widehat{\mathcal A}$ becomes weakly consistent once the regularization parameter $\upalpha$ diverges to infinity at an appropriate rate; in the next theorem, we let $\tau(\upalpha)$ be a random function that increases without bound as $\upalpha\to \infty$, which is defined by $\tau(\upalpha) = \sum_{j=1}^{\operatorname{{\mathrm{K}}}} \tau_j$, where $\tau_j = 2\sqrt{2} \max\{(\lambda_{j-1}^2-\lambda_{j}^2)^{-1}, (\lambda_{j}^2-\lambda_{j+1}^2)^{-1}\}$.

The following is an immediate consequence of Theorem (ref) and Assumption (ref).(ref).

The condition imposed on $\tau(\upalpha)$ in Theorem (ref) does not place any essential restrictions on the eigenvalues of $\mathcal C_{xz}^\ast \mathcal C_{xz}$. Given that $\tau(\upalpha)$ increases as $\upalpha$ (and thus $\operatorname{{\mathrm{K}}}$) gets larger, the condition, together with the requirement that $T^{-1}\upalpha \to 0$, merely tells us that $\upalpha$ needs to grow with a sufficiently slower rate than $T$ for the weak consistency of the FIVE. In fact, under some additional conditions, the strong (almost sure) consistency of the estimator can also be derived; we need more mathematical preliminaries to present this result, and thus leave the discussion to Section (ref) of the Supplementary Material.

Under stronger assumptions than what we require for the weak consistency of $\widehat {\mathcal A}$, we can further find that (i) the decaying rate of $ \widehat{\mathcal A}-\mathcal A$ is not uniform over the entire Hilbert space $\mathcal H$ and (ii) the choice of $\upalpha$ can affect the decaying rates of the DR and RB components in different directions. These results are given as consequences of the following asymptotic normality result of the DR component; in the theorem below, $({\mathcal C}_{xz}^\ast{\mathcal C}_{xz})_{\operatorname{{\mathrm{K}}}}^{-1}$ denotes the operator {given by $\sum_{j=1}^{\operatorname{{\mathrm{K}}}}\lambda_j^{-2}f_j\otimes f_j$ and $N(0,\mathcal G)$ denotes Gaussian random element in $\mathcal H$ with covariance operator $\mathcal G$.}

Depending on the choice of $\zeta$, $\theta_{\operatorname{{\mathrm{K}}}}(\zeta)$ may be convergent or divergent in probability (see Remark (ref) below), and thus the convergence rate of the DR component $(\widehat{\mathcal A}- \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}})\zeta$ depends on $\zeta$. This finding is not completely new; similar results were formerly observed by Mas2007 and Hu2016 in the context of functional AR(1) models. If $\theta_{\operatorname{{\mathrm{K}}}} (\zeta)$ is convergent {in probability}, then $(\widehat{\mathcal A}- \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}})\zeta$ converges at $\sqrt{T}$-rate, otherwise it converges at a slower rate given by $\sqrt{{{T}}/{{\theta_{\operatorname{{\mathrm{K}}}}(\zeta)}}}$ which is random (because of the randomness of $\operatorname{{\mathrm{K}}}$). As noted by Mas2007, this discrepancy in convergence rates implies that (i) there exists no sequence of normalizing constants $c_T$ such that $c_T (\widehat{\mathcal A}- \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}})\zeta$ weakly converges to a well defined limiting distribution uniformly in $\zeta \in \mathcal H$, and therefore it is impossible that $\widehat{\mathcal A}- \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$ weakly converges to a well defined bounded linear operator in the topology of $\mathcal L_{\mathcal H}$, and (ii) this statement is also true if $\mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$ is replaced by $\mathcal A$ (see Theorem 3.1 of Mas2007). \phantomsection{\leavevmode\color{black}Moreover, it can be deduced from Theorem (ref) that as $\operatorname{{\mathrm{K}}}$ increases, the regularization parameter $\upalpha$ induces a trade-off between the decaying rates of the DR and RB components when $\theta_{\operatorname{{\mathrm{K}}}}(\zeta)$ is not convergent. If $\operatorname{{\mathrm{K}}}$ relative to $T$ increases by a larger choice of $\upalpha$, then the operator norm of the RB component tends to shrink to zero at a faster rate. On the other hand, this change results in a faster divergence of ${\theta}_{\operatorname{{\mathrm{K}}}}(\zeta)$, and Theorem (ref) shows that the DR component will decay at a slower rate in such a case.}

In applications involving economic or statistical time series, practitioners are often interested in the marginal effect of some additive and hypothetical perturbation, say $\zeta$, in $x_t$ on $y_t$. In the considered linear model, this marginal effect is simply given by $\mathcal A\zeta$, which can be consistently estimated by $\widehat{\mathcal A}\zeta$. Let $\{\zeta_{\operatorname{{\mathrm{K}}}}\}$ be a sequence of random elements given by $\zeta_{\operatorname{{\mathrm{K}}}} = \widehat{\Pi}_{\operatorname{{\mathrm{K}}}}\zeta$. Then $\zeta_{\operatorname{{\mathrm{K}}}}$ is given by the orthogonal projection of the new perturbation $\zeta$ onto the subspace on which the sample cross-covariance of $x_t$ and $z_t$ is the most explained, in a certain sense (see Remark (ref)), among all the subspaces of dimension $\operatorname{{\mathrm{K}}}$; that is, $\zeta_{\operatorname{{\mathrm{K}}}}$ is the best linear approximation of $\zeta$ based on the covariation of the explanatory and instrumental variables. Therefore, $\zeta_{\operatorname{{\mathrm{K}}}}$ may be interpreted as a nice approximation showing how a hypothetical perturbation $\zeta$ can be revealed given the dataset, and we thus call $\zeta_{\operatorname{{\mathrm{K}}}}$ a data-supporting approximation of $\zeta$. The following is an immediate consequence of Theorem (ref): under the assumptions employed in Theorem (ref),

Based on this result, we may implement standard statistical inference on various characteristics of $\mathcal A \zeta_{\operatorname{{\mathrm{K}}}}$, which may provide a practical and interpretable insight for practitioners. We illustrate this by constructing a confidence interval for the random variable given by $\langle {\mathcal A}\zeta_{\operatorname{{\mathrm{K}}}}, \psi\rangle$ for some $\psi\in \mathcal H$. In fact, various characteristics of ${\mathcal A}\zeta_{\operatorname{{\mathrm{K}}}}$ may be written in this form; for example, if $\psi(s)=1\{s_1\leq s\leq s_2\}$ then $\langle {\mathcal A}\zeta_{\operatorname{{\mathrm{K}}}}, \psi\rangle = \int_{s_1}^{s_2} {\mathcal A}\zeta_{\operatorname{{\mathrm{K}}}}(s)ds$ means the locally (if $s_1 \neq 0$ or $s_2 \neq 1$) or globally (if $s_1=0$ and $s_2 = 1$) aggregated marginal effect on $y_t$. We then consider the interval whose endpoints are given as follows:

where $\Phi^{-1}(\cdot)$ is the quantile function of the standard normal distribution and $\widehat{\mathcal C}_{\widehat{u}\widehat{u}}=T^{-1}\sum_{t=1}^T \widehat{u}_t \otimes \widehat{u}_t$. Based on Theorem (ref) and Corollary (ref), the intervals that are repeatedly constructed as in \hyperref[{eq:conf}]{\tagform@{\ref*{eq:conf}}} are expected to include $\langle \mathcal A \zeta_{\operatorname{{\mathrm{K}}}} , \psi \rangle$ with $100(1-{\varpi}$)$\%$ of probability for a large $T$. Of course, \hyperref[{eq:conf}]{\tagform@{\ref*{eq:conf}}} may not be quite satisfactory for practitioners who want to consider a purely hypothetical perturbation $\zeta$ without data-supporting approximation; however, given that (i) the discrepancy between $\zeta$ and $\zeta_{\operatorname{{\mathrm{K}}}}$ caused by the noninvertibility of $\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz}$, which inevitably arises in our functional setting, and (ii) the magnitude of discrepancy is expected to be small since it is, anyhow, asymptotically negligible, a small bias caused by the data-supporting approximation may be understood as a cost to implement standard inference based on asymptotic normality in our setting. Furthermore and more importantly, it will be shown in Section (ref) that, if certain conditions, which are not that restrictive, are satisfied, then the convergence result given in Theorem (ref).(ref) holds even if $\mathcal A\widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$ is replaced by $\mathcal A$ (see Remark (ref)); this, of course, implies that \hyperref[{eq:conf}]{\tagform@{\ref*{eq:conf}}} can be understood as a confidence interval for $\langle \mathcal A\zeta,\psi \rangle$ with no data-supporting approximation.

Refinements of the general asymptotic results

In Section (ref), we established some general asymptotic properties of the FIVE, which do not require any specific assumptions on the eigenstructure of the cross-covariance of $x_t$ and $z_t$ other than the assumption of distinct eigenvalues. The results given in the previous section tell us that the FIVE is a reasonable estimator in this functional setting. However, what can be learned from Theorems (ref) and (ref) is not rich enough; we only know that the FIVE is consistent (Theorem (ref)) and its DR component is asymptotically normal in a pointwise sense (Theorem (ref)) if $\upalpha$ diverges to infinity at a sufficiently slow rate. We in this section investigate the asymptotic behavior of the FIVE in more detail under a set of assumptions {\leavevmode\color{black}which is stronger than Assumption (ref) but not restrictive in practice}. By doing so, we will obtain useful refinements of Theorems (ref) and (ref). The specific assumptions that we need are given as follows: in the assumption below, {we note that $\mathcal C_{xz}\mathcal C_{xz}^\ast$ allows the spectral decomposition $\mathcal C_{xz} \mathcal C_{xz}^\ast=\sum_{j=1}^\infty \lambda_j^2 \xi_{j} \otimes \xi_j$ for some orthonormal basis $\{\xi_j\}_{j\geq 1}$, and let} $\upsilon_{t}(j,\ell) = \langle x_t,f_j \rangle\langle z_t,\xi_{\ell}\rangle -\mathbb{E}[\langle x_t,f_j \rangle\langle z_t,\xi_{\ell}\rangle]$ for $j,\ell \geq 1$.

Assumptions (ref).(ref) and (ref).(ref) restrict the eigenstructure of $\mathcal C_{xz}$ (or equivalently $\mathcal C_{xz}^\ast\mathcal C_{xz}$), and these are adapted from similar conditions in Hall2007 and imaizumi2018.\footnote{In fact, Assumptions (ref).(ref) and (ref).(ref) can be replaced by the following conditions: $|\lambda_j|\leq c_{\circ} j^{-\rho/2}$ and $|\lambda_j|-|\lambda_{j+1}|\geq c_{\circ} j^{-\rho/2-1}$ for $\rho>2$. These conditions are directly comparable with similar conditions given for the eigenvalues of $\mathbb{E}[x_t\otimes x_t]$ in Hall2007 and imaizumi2018.} Assumption (ref).(ref) is a very natural condition given that $\langle \mathcal Af_j, \xi_{\ell} \rangle$ must be square-summable with respect to both $j$ and $\ell$; in this assumption, it is worth mentioning that $\varsigma$ is the parameter determining the smoothness of $\mathcal A$ on $\operatorname{ran} \mathcal C_{xz}^\ast \mathcal C_{xz} $. As may be deduced from the definition of $\upsilon_{t}(j,\ell)$ and Assumption (ref).(ref), $\{\upsilon_{t}(j,\ell)\}_{t\geq 1}$ is a stationary sequence in $\mathbb{R}$ for each $j$ and $\ell$, and the former condition of Assumption (ref).(ref) states that its lag-$s$ autocovariance function decays at a sufficiently fast rate; this condition is satisfied for a wide class of stationary processes. Note that both $\mathbb{E}[\|\langle x_t,f_j\rangle z_t\|^2]$ and $\mathbb{E}[\|\langle z_t,\xi_j\rangle x_t\|^2]$ naturally decrease as $j$ gets larger and its decaying rate is restricted by Assumption (ref).(ref). {Specifically, we require the second moments of $\|\langle x_t,f_j\rangle z_t\|$ and $\|\langle z_t,\xi_j\rangle x_t\|$ as functions of $j$ have a constant multiple of $\lambda_j^2$ as their upper envelope}; a similar condition for the iid case can be found in e.g., Hall2007.

The following theorem refines the result given in Theorem (ref) under Assumption (ref).

Some comments on the requirement $\upalpha= o(T^{\rho/(2\rho+2)})$ are first in order. This condition is needed in our proof of Theorem (ref) to deal with estimation errors associated with $\widehat{\lambda}_j$ (see Remark (ref)). This may be replaced by a sufficient and more convenient condition given by $\upalpha = o(T^{1/3})$, which does not depend on the value of $\rho$ under Assumption (ref) requiring $\rho>2$.

Theorem (ref) not only gives us a more detailed consistency result than that given in Theorem (ref), but also better clarifies how certain parameters appearing in Assumption (ref) can affect the convergence rate of the FIVE. Specifically, in the above theorem, the convergence rate is characterized by the regularization parameter $\upalpha$, the smoothness $\varsigma$ of $\mathcal A$ on $\operatorname{ran} \mathcal C_{xz}^\ast \mathcal C_{xz}$, and the decaying rate $\rho$ of $\lambda_j^2$ (as a function of $j$). From \hyperref[{eqthmconvrate}]{\tagform@{\ref*{eqthmconvrate}}}, it is evident that if $\upalpha$ grows to infinity at a sufficiently slow rate, then the convergence rate of the RB term will be dominantly determined by the second term (appearing in \hyperref[{eqthmconvrate}]{\tagform@{\ref*{eqthmconvrate}}}), whose convergence rate is positively related to $\upalpha$. Therefore, in this case, we expect a slower convergence rate of the RB component; this is a quite natural property that can also be deduced from our earlier discussion following Theorem (ref) in Section (ref). Moreover, it can be shown that the convergence rate of the FIVE is generally positively (resp.\ negative) related to $\varsigma$ (resp.\ $\rho$); the former is immediately seen from \hyperref[{eqthmconvrate}]{\tagform@{\ref*{eqthmconvrate}}}, and the latter is discussed in detail in Remark (ref).

We next refine our pointwise asymptotic normality result under Assumption (ref). To this end, it is convenient to decompose the RB component again as follows:

where $\Pi_{\operatorname{{\mathrm{K}}}} = \sum_{j=1}^{\operatorname{{\mathrm{K}}}} f_j \otimes f_j$, and this may be understood as the population counterpart of $\widehat{\Pi}_{\operatorname{{\mathrm{K}}}}$. The next theorem refines the results in Theorem (ref), {\leavevmode\color{black}but in order to simplify the subsequent discussion, we for now only consider the case where $\rho/2+2<\varsigma+\delta_{\zeta}$; the result without this condition is given in Appendix (ref).}

Theorem (ref) refines Theorem (ref) by providing a detailed asymptotic order of the RB component. Some remarks on the theorem are given in Remarks (ref) and (ref) below; particularly, in the latter remark, an improvement of the asymptotic normality result in Theorem (ref) is discussed.

Functional two-stage least square estimator

The proposed estimator

\phantomsection{\leavevmode\color{black}In the context of Euclidean space setting, our FIVE reduces to a certain IV estimator as in Remark (ref). As may be expected from the existing results, this estimator generally exhibits a larger asymptotic variance compared to the two-stage least square estimator (2SLSE), which achieves asymptotic efficiency within some context. In this section, we explore an extension of the conventional 2SLSE tailored to our functional setting.} To the best of the authors' knowledge, such an estimator has not been considered in the context of function-on-function regression, although a similar estimator was studied by Florence2015 in the case where the dependent variable is scalar-valued. For this estimator, it is assumed that $x_t$ and $z_t$ satisfy the so-called first-stage relationship:

If we consider the case $\mathcal H = \mathbb{R}^n$, then the 2SLSE is defined as follows:

and it is widely known that $\widetilde{\mathcal A}^\circ$ has many desirable properties as an estimator of $\mathcal A$. Coming back to our functional setting, it is not difficult to see that the use of the standard 2SLSE is problematic since it involves $\widehat{\mathcal C}_{zz}^{-1}$ and $(\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{zz}^{-1} \widehat{\mathcal C}_{xz})^{-1}$ which are not well defined as bounded linear operators.

To have a well-behaved analogue of the 2SLSE in our setting, we regularize those inverses as we did in Section (ref) and propose an alternative estimator. To this end, we hereafter let $\mathcal T_{\operatorname{{\mathrm{K}}}}^{-1}$ denote the regularized inverse of a compact operator $\mathcal T$ based on its first $\operatorname{{\mathrm{K}}}$ eigenelements (this is defined in the same way as $(\widehat{\mathcal C}_{xz}^\ast \widehat{\mathcal C}_{xz})_{\operatorname{{\mathrm{K}}}}^{-1}$ given in \hyperref[{eqdef3}]{\tagform@{\ref*{eqdef3}}}). Our proposed estimator is defined as follows:

and, if we let $\{\widehat{\mu}_j\}_{j\geq 1}$ (resp.\ $\{\widehat{\nu}_j\}_{j\geq 1}$) be the ordered (from the largest to the smallest) eigenvalues of $\widehat{\mathcal C}_{zz}$ (resp.\ $\widehat{\mathcal Q}$),\footnote{The eigenvalues of $\widehat{\mathcal C}_{zz}$ and $\widehat{\mathcal Q}$ are almost surely positive since they are nonnegative self-adjoint by construction.} then $\operatorname{{\mathrm{K}}}_1$ and $\operatorname{{\mathrm{K}}}_2$ are defined as

Note that by definition, $\operatorname{{\mathrm{K}}}_2 \leq \operatorname{{\mathrm{K}}}_1$ holds almost surely because $(\widehat{\mathcal C}_{zz})_{\operatorname{{\mathrm{K}}}_1}^{-1}$ is of finite rank $\operatorname{{\mathrm{K}}}_1$. We conveniently call $\widetilde{\mathcal A}$ the F2SLSE.

To investigate the asymptotic properties of the F2SLSE, it is necessary to establish some preliminary results and fix notation. First, we note that the operator $\mathcal Q$ defined by $\mathcal Q = \mathcal C_{xz}^\ast \mathcal C_{zz}^{-1} \mathcal C_{xz}$ can be understood as a well defined compact operator. To be more specific, Lemma (ref) (and the following discussion given in Section (ref) of the Supplementary Material) shows that $\mathcal C_{zz}^{-1/2} \mathcal C_{xz} = \mathcal R_{xz}\mathcal C_{xx}^{1/2}$ for some unique bounded linear operator $\mathcal R_{xz}$, which may be understood as the correlation operator of $x_t$ and $z_t$, and thus $\mathcal Q = \mathcal C_{xx}^{1/2}\mathcal R_{xz}^\ast\mathcal R_{xz} \mathcal C_{xx}^{1/2}$. From similar arguments, it can be easily shown that the operator $\mathcal P = \mathcal C_{yz}^\ast \mathcal C_{zz}^{-1} \mathcal C_{xz}$ is also well defined. We then let $\{\mu_j,g_j\}_{j\geq 1}$ (resp.\ $\{\nu_j,h_j\}_{j\geq 1}$) be the eigenelements of $\mathcal C_{zz}$ (resp.\ $\mathcal Q$), i.e.,

Since $\mathcal C_{zz}$ and $\mathcal Q$ are self-adjoint and nonnegative, $\mu_j$ and $\nu_j$ are all nonnegative. We know from \hyperref[{eqmodel}]{\tagform@{\ref*{eqmodel}}} that the population relationship $\mathcal P=\mathcal A \mathcal Q$ holds.

General asymptotic properties

This section discusses the asymptotic properties of the F2SLSE under the following assumption:

It should be noted that the operator $\mathcal A$ satisfying $\mathcal P=\mathcal A\mathcal Q$ is uniquely identified if Assumptions (ref) and (ref) are satisfied. To see this in detail, note that $\mathcal A$ satisfying $\mathcal C_{yz}^\ast \mathcal C_{xz} = \mathcal A \mathcal C_{xz}^\ast \mathcal C_{xz}$ (resp.\ $\mathcal P=\mathcal A \mathcal Q$) is identified if and only if $\mathcal C_{xz} $ (resp.\ $\mathcal C_{zz} ^{-1/2} \mathcal C_{xz}$) is injective. If Assumption (ref) is satisfied and hence $\mathcal C_{zz} $ is injective, then $\ker \mathcal C_{xz} (=\ker \mathcal C_{xz}^\ast \mathcal C_{xz})$ becomes identical to $ \ker \mathcal C_{zz} ^{-1/2} \mathcal C_{xz} (= \ker \mathcal Q)$.

As in Section (ref), we also consider the following decomposition:

where $\widetilde{\Pi}_{\operatorname{{\mathrm{K}}}_2}$ denotes the orthogonal projection defined by $ \widetilde{\Pi}_{\operatorname{{\mathrm{K}}}_2} = \sum_{j=1}^{\operatorname{{\mathrm{K}}}_2} \widehat{h}_j \otimes \widehat{h}_j$ and $\{\widehat{h}_j\}_{j=1}^{\operatorname{{\mathrm{K}}}_2}$ is the collection of the eigenvectors of $\widehat{\mathcal Q}$ corresponding to the first $\operatorname{{\mathrm{K}}}_2$ leading eigenvalues. The two terms in \hyperref[{eqdecom2}]{\tagform@{\ref*{eqdecom2}}} are similarly interpreted as in the case of the FIVE (see \hyperref[{eqdecom1}]{\tagform@{\ref*{eqdecom1}}}), and we thus call the first (resp.\ the second) term the DR (resp.\ RB) component.

We first show that both the DR and RB components are asymptotically negligible (and thus $\widehat{\mathcal A}$ is weakly consistent) if the regularization parameters $\upalpha_1$ and $\upalpha_2$ diverge to infinity at appropriate rates: in the theorem below, we let $\tau_{1,j} = 2\sqrt{2} \max\{(\mu_{j-1}-\mu_{j})^{-1},(\mu_{j}-\mu_{j+1})^{-1}\}$ and $\tau_{2,j} = 2\sqrt{2} \max\{(\nu_{j-1}-\nu_{j})^{-1}, (\nu_{j}-\nu_{j+1})^{-1}\}$.

An immediate consequence of Theorem (ref) is given as follows:

As in the case of the FIVE, the conditions imposed on the quantities $(\sum_{j=1}^{\operatorname{{\mathrm{K}}}_1} \mu_j \tau_{1,j})(\sum_{j=1}^{\operatorname{{\mathrm{K}}}_2}\tau_{2,j})$ and $(\sum_{j=\operatorname{{\mathrm{K}}}_1+1}^{\infty}\mu_j)(\sum_{j=1}^{\operatorname{{\mathrm{K}}}_2}\tau_{2,j} )$ are understood not as special restrictions on the eigenvalues, but as requirements on the growing rates of $\upalpha_1$ and $\upalpha_2$ in our asymptotic theory. Specifically, the condition on the former quantity merely requires $\upalpha_1$ and $\upalpha_2$ to increase slowly so that $\operatorname{{\mathrm{K}}}_1$ and $\operatorname{{\mathrm{K}}}_2$ tend to grow with sufficiently slower rates than $T$. Moreover, given that, for fixed $\upalpha_2$, $(\sum_{j=\operatorname{{\mathrm{K}}}_1+1}^{\infty}\mu_j)(\sum_{j=1}^{\operatorname{{\mathrm{K}}}_2}\tau_{2,j} )$ can be arbitrarily small by choosing $\upalpha_1$ large enough, the condition on the latter quantity merely tells us that the growing rate of $\upalpha_1$ needs to be sufficiently higher than that of $\upalpha_2$. In addition to the weak consistency given by Theorem (ref), the strong consistency of the F2SLSE can be derived under additional conditions; this is discussed in Section (ref) of the Supplementary Material.

We also obtain an asymptotic normality result similar to that given by Theorem (ref) for the FIVE:

As shown, we need more stringent requirements on the growing rates of $\upalpha_1$ and $\upalpha_2$. This is mainly due to that the F2SLSE involves the doubly regularized inverse $\widehat{\mathcal Q}_{\operatorname{{\mathrm{K}}}_2}^{-1}$, which may be ill-behaved if $\upalpha_1$ and $\upalpha_2$ do not grow at sufficiently slow rates. This implies that there is no reason why the F2SLSE is generally preferred to the FIVE in this functional setting, unlike what we can expect from the general preference for the 2SLSE by practitioners in the Euclidean space setting.\footnote{Moreover, the two estimators have different RB terms whose magnitudes depend on various parameters (e.g., the eigenvalues of $\mathcal C_{xz}^\ast\mathcal C_{xz}$ and $\mathcal Q$), and thus the RB component of the FIVE can have a smaller asymptotic order. }

Refinements of the general asymptotic results

We provide refinements of Theorems (ref) and (ref) {\leavevmode\color{black}under the following set of assumptions which is stronger than Assumption (ref)}: below, we let $\tilde{\upsilon}_{t}(j,\ell) = \langle z_t,g_j \rangle\langle z_t,g_{\ell}\rangle -\mathbb{E}[\langle z_t,g_j \rangle\langle z_t,g_{\ell}\rangle]$ for $j,\ell \geq 1$.

The conditions are somewhat similar to those in Assumption (ref), and thus we omit detailed comments except the following two points: (i) from a technical point of view, Assumption (ref).(ref) is similar to Assumption (ref).(ref) employed for our study of the FIVE and helps us obtain convergence rates of the eigenelements of $\widehat{\mathcal Q}$ (which are crucial inputs to our main results of the F2SLSE), and {\leavevmode\color{black}(ii) we require a smoothness condition on \(\mathcal B\) which characterizes the linear relationship between \(x_t\) and \(z_t\) whereas such a condition is not necessary in the case of the FIVE. This reveals that the data generating process (DGP) is more restricted for our asymptotic analysis of the F2SLSE.}

Our next result refines Theorem (ref) by providing a more detailed result on the RB component.

The convergence rate of the RB component, described in the above theorem, depends not only on the regularization parameters, but also on smoothness of $\mathcal A$ as in the case of the FIVE. However, the convergence rate described in \hyperref[{eqthmconvrate2}]{\tagform@{\ref*{eqthmconvrate2}}} is generally slower than that of the FIVE, and this is somewhat expected from the fact that the F2SLSE involves a doubly regularized (and thus less stable) inverse. Despite this disadvantage of the F2SLSE over the FIVE, our simulation results support that the F2SLSE performs comparably well among a set of the competing estimators (including the FIVE), and thus this estimator can also be used in practice.

Using Assumption (ref), the next theorem refines Theorem (ref), but as in Section (ref), we for now only focus on the case where $\rho_{\nu}/2+2<\varsigma_{\nu}+\delta_{\zeta}$. The result without this condition is provided in the Supplementary Material (see Section (ref)). In the theorem below, we, as in \hyperref[{eqdecom1}]{\tagform@{\ref*{eqdecom1}}}, consider the decomposition of the RB component given by

where $\Pi_{\operatorname{{\mathrm{K}}}_2} = \sum_{j=1}^{\operatorname{{\mathrm{K}}}_2} h_j \otimes h_j$ is understood as the population counterpart of $\widetilde{\Pi}_{\operatorname{{\mathrm{K}}}_2}$.

Obtaining a pointwise asymptotic normality result that is not dependent on the RB component as in Remark (ref) requires more stringent conditions, which will be detailed in Remark (ref) below.

Numerical Studies

We first investigate the finite sample performance of our estimators via Monte Carlo studies. In Sections (ref)--(ref), the number of replications is set to 1,000 and all the considered random variables are demeaned before computing the estimators of $\mathcal A$. Section (ref) provides an empirical application.

Experiment 1: Functional linear simultaneous equation model

We consider the following functional linear simultaneous model: for $t\geq 1$,

where $u_t =0.8 v_t + 0.6 \varepsilon_t$, $ \{v_t\}_{t \geq 1}$ and $\{\varepsilon_t\}_{t \geq 1} $ are mutually independent iid sequences of standard Brownian bridges satisfying $\mathbb{E}[v_t\otimes \varepsilon_{\ell}] =0$ for all $t,\ell \geq 1$. The constant $\vartheta$ is chosen in such a way that the first-stage functional coefficient of determination (see, Yao2005), defined by $\mathbb E[\Vert \vartheta \mathcal B z_t \Vert ^2]/\mathbb E [\Vert x_t\Vert^2]$, has a specific value of $\mathtt{r}^2$. In this section, we will focus on empirical MSEs of a few estimators at various levels of $\mathtt{r}^2$, and in particular we consider $\mathtt{r}^2 \in \{ 0.1,0.2,\ldots,0.5 \}$.

The DGP here is specially designed to examine the performance of our estimators when all the employed assumptions (Assumptions (ref), (ref), (ref), (ref), (ref) and (ref)) are satisfied (see Section (ref) of the Supplementary Material). Specifically, we let $\{z_t\}_{t \geq 1}$ be an iid sequence of standard Brownian bridges satisfying $\mathbb{E}[z_t\otimes v_t]=\mathbb{E}[z_t\otimes u_t]=0$. Then, we have $\mu_j = (j\pi)^{-2} $ and $g_j(s) = \sqrt{2}\sin(j\pi s)$ for $s \in [0,1]$, see, e.g., JaimezBonnet. The operators $\mathcal A$ and $\mathcal B$ are defined as follows:

In this setup, $f_j = g_j$. In view of the fact that function-valued random variables are only partially observed in practice, we assume that the discrete realizations of $y_t$, $x_t$ and $z_t$ at 50 equally-spaced points of $[0,1]$ are available. Then, following the literature, e.g., Ramsay2005, we represent functional variables $y_t$, $x_t$ and $z_t$ by using 31 Fourier basis functions.

We will compare the performance of our estimators with the ridge regularized IV estimator (RIVE) of Benatia2017 with denoting their regularization parameter to $\upalpha^{-1}$ to keep notational consistency. To compute the FIVE and RIVE, we consider $ \delta_{\upalpha}{T^{-0.4}}\Vert \widehat{\mathcal C}_{xz}\Vert_{\operatorname{HS}}^2 $ as candidates for the inverse of $\upalpha$. This candidate value is calculated at 20 equidistant points of $\delta_{\upalpha}$ ranging from $0.1$ to $T^{0.2}$. Among such candidates, we choose the value that minimizes the empirical MSE of each estimator. The F2SLSE needs two regularization parameters: $\upalpha_{1}$ and $\upalpha_{2}$. The parameter $\upalpha_1$ is chosen as the FIVE and RIVE with $\Vert \widehat{\mathcal C}_{xz}\Vert_{\operatorname{HS}} ^2$ being replaced by $\Vert \widehat{\mathcal C}_{zz} \Vert_{\operatorname{HS}}^2$. Once $\upalpha_1$ is chosen, we similarly choose the inverse of $\upalpha_2$ from $\delta_{\upalpha_2} (\upalpha_1^{-1}\Vert \mathcal C_{zz}\Vert_{\operatorname{HS}}^2)^{1/2}\Vert \widehat{\mathcal Q}_{\operatorname{{\mathrm{K}}}_1}\Vert_{\operatorname{HS}} ^2 $ with $\delta_{\upalpha_2}$ being 20 equidistant points between $T^{0.05}$ and $ T^{0.2}$. This setup enforces $\upalpha_1$ to grow at a faster rate than that of $\upalpha_2$.

To save space, we report estimation results only for the case with $T=500$; the results with a smaller sample size are qualitatively similar and are reported in Section (ref) of the Supplementary Material. Figure (ref) reports boxplots (without outliers) of the empirical MSE estimated with the FIVE (red), the F2SLSE (blue) and the RIVE (green). The first interesting observation in the figure is that our estimators tend to produce smaller MSEs when the signal from $x_t$ to $y_t$ is more concentrated on the first few components, i.e., when $n_a =5$. This observation is consistent regardless of the values of $n_b$ and $\mathtt{r}^2$. This may not be surprising because when $n_a$ is large, the first few $f_j$'s (=$g_j$'s) summarize the most significant information of $\mathcal A$.

In the figure, as the value of $\mathtt{r}^2$ decreases, the considered estimators tend to exhibit larger MSEs. Similar observations can be found in the standard IV literature, in which the so-called concentration parameter is used to measure the strength of IVs. Given that the coefficient of determination $\mathtt{r}^2$ is closely related to the concentration parameter in the IV literature, such a larger MSE may be understood as the distortion related to weak instruments.

In the subsequent sections, we will consider a more general setting to investigate the robustness of our estimators when the assumptions are unlikely to hold. Unlike the DGP considered in this section, verifying if all the required conditions are satisfied for the DGPs under consideration becomes nearly impossible. Given the practical challenges of confirming these conditions, practitioners may find it valuable to observe the performance of our estimators in the presence of potential violations of the required conditions.

Experiment 2: a modification of \citepos{Benatia2017} simulation DGP

In this section, we consider a simulation DGP similar to that in \citepos{Benatia2017}. Specifically, we let $\mathcal B$ in \hyperref[{eqsimdgp}]{\tagform@{\ref*{eqsimdgp}}} be the identity operator $\mathcal I$ and let $\mathcal A$ be the integral operator with kernel $\kappa_{\mathcal A} (s_1, s_2) =1-|s_1-s_2|^2$ for $s_1,s_2 \in [0,1]$. In this setup, the first-stage signal is solely determined by the constant $\vartheta$. The IV $z_t$ is given as follows:

where $a_t$ and $b_t$ are randomly drawn from the uniform distribution $\text{U}[2,5]$ for each $t$. That is, $z_t$ is obtained by adding an additive noise $\eta_t$ to the beta density function with parameters $a_t$ and $b_t$. The IV in \hyperref[{eqsimu2}]{\tagform@{\ref*{eqsimu2}}} is analogous to that used for the simulation experiments in Benatia2017, in which the additive noise $\eta_t(s)$ is given by $q_t$ for all $s\in [0,1]$ with $q_t$ being randomly drawn from $\text{N}(0,1)$. In this section, we allow a more general form of $\eta_t$ by letting $\eta_t = \sum_{j=1} ^{n_J} \sigma_{j} q_{t,j} \xi_j $ where $n_J=31$, $\{\xi_j\}_{j\geq1}$ is the Fourier basis functions with the constant basis function $\xi_1$, and $q_{t,j}\sim_{\text{iid}}\text{N}(0,1)$ across $t$ and $j$. \citepos{Benatia2017} IV can be understood to the case $n_J = 1$ under our notation. Then, we consider three different designs of $\{\sigma_j\}_{j\geq1}$. Firstly, we consider the case $\sigma_{j} = c_1 \sigma_\eta$ for $j \leq 2$ and $\sigma_{j} = c_1 \sigma_\eta (0.1)^{j-2}$ for $j > 2$; this is called the sparse design. Secondly, we set $\sigma_{j} = \sigma_\eta (0.9)^{j-1}$ and call this setting the exponential design. In the last design, which we call the geometric design, we let $\sigma_j =c_2 \sigma_\eta j^{-1}$. The parameter $\sigma_\eta$ is set to $0.5$ and $0.9$, and the constants $c_1$ and $c_2$ are chosen in such a way as to have the same Hilbert-Schmidt norm of $\mathbb E[\eta_t \otimes \eta_t]$ in all three designs. Lastly, the parameter $\vartheta$ is chosen as in Section (ref) with $\mathtt{r}^2$ being set to 0.5; this can be done by using that $\mathbb{E}[\|v_t\|^2] = 1/6$, $\mathbb{E}[\|\eta_t\|^2] = \sum_{j=1}^{31}\sigma_j^2$, and the value of $\mathbb E [\Vert \tilde z_t(s_i; a_t, b_t)\Vert ^2]$ can be approximated from a large number of simulations.

Table (ref) summarizes simulation results. Overall, the MSEs of our estimators and those of the RIVE are similar to each other. However, in the exponential design, our estimators tend to have smaller MSEs compared to the RIVE.

We note that our estimators and related asymptotic results can be used to discuss the coverage probability of the interval \hyperref[{eq:conf}]{\tagform@{\ref*{eq:conf}}} that is computed from the FIVE. The interval is expected to contain the random quantity $\langle \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}} \zeta, \psi \rangle$ with ($100-\varpi$)% of probability; moreover, if certain conditions are satisfied (see Remark (ref)) the interval \hyperref[{eq:conf}]{\tagform@{\ref*{eq:conf}}} can be understood as the ($100-\varpi$)% confidence interval for $\langle \mathcal A \zeta, \psi \rangle$ which is nonrandom. Based on Theorem (ref) (and also Remark (ref)), we may construct a similar interval with the F2SLSE and the interval is expected to include $\langle \mathcal A \widetilde{\Pi}_{\operatorname{{\mathrm{K}}}_2} \zeta, \psi \rangle$ (and also $\langle \mathcal A \zeta, \psi \rangle$ under certain conditions) with ($100-\varpi$)% of probability; the coverage of this confidence interval will also be examined in this experiment. In order to compute the coverage probabilities, we let $\psi=\ell_1$ and let $\zeta$ be randomly generated by $\zeta = \sum_{j=1} ^{11} \ddot q_{1,j} \ell_j$ for each realization of the DGP, where $\{ \ell_j \}_{j\geq1}$ is the polynomial basis with the constant basis function $\ell_1$ and $\ddot q_{j}\sim_{\text{iid}}\text{N}(0,j^{-4})$ across $j$. The simulation results are reported at the bottom of Table (ref). We note that, in all the considered cases, the coverage probabilities for $\langle \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}} \zeta, \psi \rangle$ or $\langle \mathcal A \widetilde{\Pi}_{\operatorname{{\mathrm{K}}}_2} \zeta, \psi \rangle$ are close to the nominal level, which supports our findings in Theorems (ref) and (ref). Moreover, even if the reported coverage probabilities for $\langle \mathcal A \zeta, \psi \rangle$ tend to be worse than those for $\langle \mathcal A \widehat{\Pi}_{\operatorname{{\mathrm{K}}}} \zeta, \psi \rangle$ or $\langle \mathcal A \widetilde{\Pi}_{\operatorname{{\mathrm{K}}}_2} \zeta, \psi \rangle$, they are still reasonably close to the nominal level of 95%. This is what can be expected from Remarks (ref) and (ref). In unreported simulations, we further experimented with different choices of $\zeta$ and $\psi$, but found no significant difference.

Experiment 3: AR(1) model of probability density functions

In this section, we examine the performance of the proposed estimators in the AR(1) model of probability density functions. What is mainly different from the earlier experiments given in Sections (ref) and (ref) is that endogeneity is not explicitly imposed, but implicitly introduced by estimation errors.

We let $\{p_t ^\circ\}_{t \geq 1}$ be a sequence of probability densities supported on $[0,1]$, and consider the linear prediction model of $p_t ^\circ$ given $p_{t-1}^\circ$. Each density may be treated as a random variable taking values in $\mathcal H$, but the collection of probability densities in $\mathcal H$ is not a linear subspace. As a result, a direct application of the statistical methods developed in a Hilbert space setting may not be recommended; see e.g., delicado2011dimensionality, petersen2016, Hron2016330, kokoszka2019forecasting, and zhang2020wasserstein. As a way to circumvent such issues, we consider the centered-log-ratio (clr) transformation $ y_t^\circ(s) = \log p_t^\circ(s) - \int \log p_t^\circ(s) ds$, $s \in [0,1]$ (see e.g., Egozcue2006). Then, $\{y_t^\circ\}_{t\geq 1}$ turns out to be a sequence in $\mathcal H_c$, the collection of all $\zeta \in \mathcal H$ satisfying $\int_{0}^1 \zeta(s)ds = 0$, and $\mathcal H_c$ is obviously a Hilbert space. Any element in $\mathcal H_c$ may be understood as a probability density via the inverse transformation $y_t^{\circ}(s) \mapsto \exp(y_t^{\circ}(s))/\int_{0}^{1} \exp(y_t^{\circ}(s))ds$. Thus, the linear prediction model of $p_t^\circ$ given $p_{t-1}^\circ$ may be recast into that of $y_t ^\circ$ given $y_{t-1} ^\circ$ in $\mathcal H_c$. We thus consider the following prediction model:

where $y_{t-1}^\circ$ and $\varepsilon_t$ are uncorrelated. To mimic situations commonly encountered in practice, we assume that $p_t^\circ$ (and thus $y_t^\circ$) is not observed, but only random samples $\{s_{i,t}\}_{i=1}^{n_t}$ drawn from $p_t^\circ$ are available. If so, by replacing the density $p_t^\circ$ or the log-density $\log p_t^\circ$ with its proper nonparametric estimate, we may obtain an estimate $y_t$ of $y_t^\circ$, and then, as shown in Example (ref), $\{y_t\}_{t\geq 1}$ satisfies $y_t = c_y + \mathcal A (y_{t-1}-c_y) + u_t$, but now $y_{t-1}$ and $u_t$ are generally correlated due to errors arising from the nonparametric estimation. We will compute the FIVE and F2SLSE by assuming that $y_{t-2}$ is a proper IV, as in Example (ref). Of course, this assumption may not be true depending on how estimation errors are generated. Even with this possibility, it may be of interest to practitioners, who are very often have no choice but to replace $p_t^\circ$ or $\log p_t^\circ$ with a standard nonparametric estimate, to see if a naive use of our estimators can make any actual improvements in estimating $\mathcal A$. This is the purpose of simulation experiments in this section.

Specifically, we first estimate $y_t^\circ$ from $n$ random samples that are generated from $p_t^\circ$ by (i) the local likelihood density estimation method proposed by loader1996 (see Appendix (ref) for more details) and (ii) the standard kernel density estimation method with the Gaussian kernel and Silverman's rule-of-thumb bandwidth silverman2018density. Even if the former is more suitable for estimation of $y_t^\circ$ SEO2019, the latter is considered as well because of its popularity in empirical studies. Once $y_t$ is computed, it is represented by the first 30 nonconstant Fourier basis functions for implementation of the FPCA in $\mathcal H_c$. We let $c_y$ be the clr transformation of the normal density function with mean $0.5$ and variance $0.25^2$ that is truncated on $[0,1]$. In addition, $\varepsilon_t = \sum_{j=1} ^{\infty} \sigma_{j} q_{t,j} \xi_j^c$, where $\{\xi^c_j\}_{j\geq1}$ is the Fourier basis functions except for the constant basis function, and $q_{t,j}\sim_{\text{iid}}\text{N}(0,1)$ across $t$ and $j$.\footnote{In actual computation, $\varepsilon_t$ can be approximated by $ \sum_{j=1}^{L} \sigma_{j} q_{t,j} \xi_j^c$ for some large $L$. We set $L$ to 50 in this example and found no significant difference even from big changes in $L$ as long as $L \geq 50$.} Below we consider two different specifications of $\sigma_j$, which are respectively called the exponential design and the sparse design; in the exponential design, $\sigma_{j} = 0.1 (0.9)^{j-1}$, and in the sparse design, $\sigma_{j} = c_{\sigma}$ for $j \leq 2$ and $\sigma_{j} =c_{\sigma}(0.1)^{j-2}$ for $j > 2$, where $c_{\sigma}$ is chosen so that the Hilbert-Schmidt norms of $\mathbb{E}[\varepsilon_t \otimes \varepsilon_t]$ in both designs are equal. These two designs are respectively obtained by setting $\sigma_{\eta}$ to $0.1$ in the sparse and exponential designs considered for $\eta_t$ in Section (ref), and the reason why we choose a relatively smaller scale of $\sigma_j$ in this experiment is only to avoid as much as possible that the simulated densities have shapes that are rarely observed in practice (e.g., densities that are U-shaped or highly multimodal). We let $\mathcal A$ be defined by $\sum_{j=1}^\infty a_j \xi_j ^c \otimes \xi_j ^c$,\footnote{$\mathcal A$ is approximated by $\sum_{j=1}^{50} a_j \xi_j ^c \otimes \xi_j ^c$ in actual computation as in the case of $\varepsilon_t$} and, for each realization of the DGP, the coefficients $\{a_j\}_{j\geq 1}$ are independently determined across $j$ as follows,

Note here that we let the first two coefficients $a_{1}$ and $a_{2}$ be bounded below by $0.4$, which is to ensure that the operator norm of the cross-covariance operator of $y_{t-1}^{\circ}$ and $y_{t-2}^\circ$ is bounded away from zero. If this quantity is close to zero, then the employed IV {may }become `weak' and this case is not considered in the present paper.

Table (ref) reports the empirical MSEs of the proposed estimators and \citepos{Park2012} FLSE when $n=100$ and $150$ (recall that $n$ is the number of random samples drawn from the distribution $p_t ^\circ$ to estimate $\log p_t^\circ$ or $p_t^\circ$). The IV estimators tend to exhibit smaller empirical MSEs than the FLSE. The superior comparative performance of the IV estimators is more noticeable when $n$ is small and $T$ is large. This is what can be conjectured from our earlier discussion; as $n$ gets smaller, $y_t$ becomes a less accurate estimate of $y_t^\circ$, and hence the estimators that address the possible endogeneity caused by estimation errors will work better. The FIVE or the F2SLSE exhibits the smallest MSE in most of the cases (see also Table (ref) in Appendix reporting the simulation results for the case where the lower bound of $a_{1}$ and $a_{2}$ increases to $0.6$). However, it is hard to conclude the relative performance between the IV estimators; this may depend on various factors such as the DGP and the method of density estimation. Thus, it would be advisable to use those IV estimators complementarily in practice.

Empirical application: effect of immigration on native wages

In this section, we use our estimation methods to investigate the effect of immigrant inflows on the labor market outcomes of workers with heterogeneous skills, which has received due attention from both researchers and policymakers, see, e.g., Card2009, borjas2011, Ottaviano2011, and Glitz2012. To begin with, we use national level data and generalize a widely used empirical model by viewing the variables of interest as functions depending on a measure of relative communication skill provision. Our measure of relative communication skill provision is similar to \citepos{Peri2009} measure of occupation-specific relative provision of communication versus manual skills. The number of distinct skill levels, denoted $s_j$, is 223, and, by construction, each occupation is uniquely identified by the skill score $s_j\in[0,1]$. Its formal definition is provided in the Supplementary Material.

We merge the percentile scores of relative communication skill provision to individuals in the monthly CPS data running from January 1996 to December 2019. The CPS data, which can be downloaded from the Integrated Public Use Microdata Series (IPUMS)\nocite{Flood2020}, provide information on various characteristics of individuals: hourly wage, citizenship status, age, employment status, and occupation. We focus on individuals who (i) are aged between 18 and 64 years, (ii) are not self-employed, and (iii) have positive income. Immigrants are defined by those who are not a citizen or are a naturalized citizen. The skill-dependent labor supply of immigrants ($\ell^\circ_{it}(s_j)$) and that of natives ($\ell^\circ_{nt}(s_j)$) are computed by the total hours of work per week (weighted by the variable WTFINL) provided by the foreign- and native-born workers for each $s_j$. The skill-dependent native wage is computed by weighted averaging weekly wages of native workers\footnote{The weekly wage of a native worker is computed as $(\text{hourly wage})\times (\text{usual hours of work})$, and the variables required to compute this quantity are also available in the CPS. We use the variable EARNWT as a weight.} in the occupation corresponding to $s_j$, and its logged value ($w_t^\circ(s_j)$) is used for the analysis.

The empirical models used in the labor economics literature (e.g., Dustmann2012; Sharpe2020101902) can be written as follows: $ \Delta w_{t} ^\circ (s_j) = \beta_j^\circ \Delta h_t ^\circ(s_j) + u_t^\circ (s_j), $ where $\Delta w_t ^\circ (s_j) = w_t^\circ(s_j)-w_{t-1}^\circ(s_j)$, $u_t^\circ (s_j)$ denotes the disturbance term, $\beta_j^\circ$ is the parameter of interest, the explanatory variable $\Delta h_t ^\circ(s_j)$ is the first difference of $h_t ^\circ (s_j)$, and $h_t^\circ (s_j)= \ell_{it}^\circ (s_j) /( \ell_{nt}^\circ (s_j)+\ell_{it}^\circ (s_j))$. In the above model, an inflow of immigrants in the occupation with $s_j$ is assumed to affect only the wages of natives in the occupation requiring the same skill level, which seems to be restrictive. To resolve this issue, one may instead allow spillover effects across occupations, but this requires researchers to estimate too many parameters; for example, if we allow a spillover effect from the occupation corresponding to $s_i$ to another occupation corresponding to $s_j$ for any arbitrary $i,j\in \{1,\ldots,223\}$, then there are $223^2$ elements to be estimated. As an alternative, we view observations $w_t ^\circ(s_j)$ and $h_t^\circ (s_j)$ for each $t$ as imperfect realizations of curves $w_t$ and $h_t$, and use our methodology developed in the previous sections. To this end, we first estimate each of those curves with the standard Nadaraya-Watson estimator employing the second-order Gaussian kernel and the bandwidth minimizing the least square cross validation criterion. The smoothed curves are represented by 15 cubic B-Spline functions and are denoted by $ w_t $ and $ h_t $, respectively. Then, we estimate the following model:

where $\Delta w_t=w_t-w_{t-1}$, $\Delta h_t=h_t-h_{t-1}$, and $\Delta h_t$ is likely to be correlated with $u_t$ due to, e.g., the self-selection bias pointed out by borjas1987 and llull2018. Thus, we use the changes in the imputed share of immigrants as an IV, which has been employed in various contexts, including Card2009, Peri2009, and david2013growth. Specifically, the imputed share of immigrants in the occupation corresponding to $s_j$, denoted $z_t ^\circ (s_j)$, is defined as follows:

where $b$ denotes the country of birth of immigrants, $\ell ^\circ _{it1994,b} (s_j)$ is the labor supply of immigrants in the occupation corresponding to $s_j$ from the country $b$ in the month $t$ of the year 1994, and $\ell ^\circ_{it1994,b}$ is its aggregation over $s_j$. The curve of imputed shares of immigrants, denoted $ z_t$, is obtained by smoothing $ z_t ^\circ (s_j)$, and the instrument, denoted $\Delta z_t$, is the first difference of $ z_t$.

The smoothed curves are reported in Figure (ref). The solid lines in Figure (ref) indicate the mean functions of $ w_t$, $ h_t$, and $ z_t$. Figure (ref) shows that native workers tend to be better paid if they are in occupations needing relatively higher communication skills. On the other hand, the share of immigrants tends to decrease in such occupations, and so does the imputed share of immigrants; this may be because natives have a comparative advantage in communication intensive tasks.

Then, we apply our estimation method to study if an inflow of immigrants has a heterogenous impact to native workers depending on the value of $s$. For the ease of interpretation, we focus on the case where immigrants are fully concentrated in a group of occupations with low, medium or high communication skill intensity, but they are evenly distributed within the group; that is, we set $\zeta$ in Theorem (ref) to $1\{ 0 \leq s < 1/3 \}$, $1\{ 1/3\leq s < 2/3 \}$, and $1\{ 2/3 \leq s < 1 \}$. The regulariation parameter is chosen as in Section (ref). The results computed from the FIVE are summarized in Figure (ref). The estimation results from the F2SLSE are similar and thus omitted. Overall, our findings in Figure (ref) reconfirm the existing evidence that an inflow of immigrants heterogeneously affects the labor market outcomes of native workers according to workers' skills. For example, in Figure (ref), if the share of immigrants increases in occupations with a low value of $s_j$, then native wages are overall positively affected, although the size of this effect depends on the value of native workers' $s_j$ as well. In particular, Figure (ref) suggests that native workers in the occupations with $s \in[0.1,0.4]$ will experience the most significant positive wage effects. This is somewhat consistent with \citepos{Peri2009} finding that native workers in occupations intensive in manual skills take advantage of having better-paid jobs when similarly skilled immigrants enter into the market. On the other hand, in Figure (ref), it seems that the natives in occupations with $s \in [0.1,0.4]$ are negatively affected if the share of immigrants increases in occupations of which $s\in[1/3,2/3]$.

Before concluding this section, recall that the earlier literature mostly relies on the strategy of reducing dimensionality of the model by classifying workers into a few groups according to a measure of their skill. We will make a comparison of our estimation results with those based on this strategy. To this end, we let $\zeta_{j(J)} = J \times 1\{ {j -1}< J s \leq {j } \}$. Then, the inner product $\langle \Delta w_t, \zeta_{j(J)}\rangle$ computes the average of changes in the log wages of natives in the occupations of which $s$ is between $({j}-1)/J$ and ${j}/J$, at time $t$. We then estimate the following using the standard 2SLSE,

where ${\Delta w_{ t(J) }} = ( \langle \Delta w_t, \zeta_{1(J)}\rangle , \ldots, \langle \Delta w_t , \zeta_{J(J)}\rangle )' $ and ${\Delta h_{ t(J) }} = ( \langle \Delta h_t, \zeta_{1(J)}\rangle, \ldots, \langle \Delta h_t , \zeta_{J(J)}\rangle )'$. The IV that is used to compute the 2SLSE is $ {\Delta z_{ t(J) } } = ( \langle \Delta z_t, \zeta_{1(J)}\rangle , \ldots, \langle \Delta z_t , \zeta_{J(J)}\rangle )'$. For example, if $J=3$, we have $\langle \Delta w_t , \zeta_{1(3)}\rangle $, $\langle \Delta w_t , \zeta_{2(3)}\rangle $, and $\langle \Delta w_t, \zeta_{3(3)} \rangle$, which are plotted in Figure (ref). In the following, we consider three different values for $J$: $3,$ $7$, and $11$.

We compare the performance of estimators by using the root mean squared prediction error (RMSPE), which is computed with a rolling window for three test sets, with setting their starting points respectively as 2013/01, 2015/01, and 2017/01. Specifically, if we let \(\bar{u}_h\) be the forecasting error computed from the FIVE, F2SLSE or RIVE, then the RMSPE is given by \((H^{-1}\sum_{h=1}^H \int \bar{u}_h(s)^2 ds)^{1/2}\) and the regularization parameters are chosen in such a way as to minimize the RMSPE. Let $\widehat{u}_{h,j(J)}$ denote the $j$th element of the forecasting error $\widehat{u}_{h(J)}$ computed from the 2SLSE for each $J$. The standard RMSPE of the 2SLSE is given by $(H^{-1}\sum_{h=1} ^H \sum_{j=1} ^J\widehat{u}_{h,j(J)}^2)^{1/2}$. {Because this measure is nondecreasing in $J$ and thus does not provide a fair comparison between RMSPEs,} we instead consider the normalized RMSPE, given by $((JH)^{-1} \sum_{h=1} ^H \sum_{j=1} ^J\widehat{u}_{h,j(J)}^2)^{1/2}$, which can be reasonably compared to the RMSPEs of our estimators since, for each $h\geq 1$, (i) both $\widehat{u}_{h,j(J)}$ and $\langle \bar{u}_h, \zeta_{j(J)} \rangle$ are estimates of the local average of $u_h$ over the interval $[(j-1)/J,j/J]$ and (ii) $\int \bar{u}_h(s)^2 ds$ may be approximated by $J^{-1} \sum_{j=1}^J \langle \bar{u}_h, \zeta_{j(J)} \rangle^2$.

Estimation results are reported in Table (ref). We first note that the results from our estimators and those from the RIVE are very similar to each other. These estimators report smaller RMSPEs than those of the 2SLSE except for the case $J=3$. Even if the 2SLSE reports the smallest RMSPEs when $J=3$, we should note that in this case, 223 different levels of skill are aggregated into only three groups, resulting in a lot of information loss. Moreover, the normalized RMSPE of the 2SLSE rapidly increases as we consider more finely defined skill groups. This may be because, as $J$ gets larger, the number of parameters to be estimated rapidly increases. In addition, for a large $J$, the sample (cross-)covariance matrices for computing the 2SLSE tend to be singular, and thus the 2SLSE is expected to perform poorly. {This result also suggests that the pre-classification strategy can have a significant effect on the estimation results and our interpretation of them.} Therefore, the results given by Table (ref) imply that our functional IV methodology can be an appealing alternative to practitioners.

{Conclusion}

\phantomsection{\leavevmode\color{black}This paper extends the existing results on the endogenous functional linear model to a more general setup, allowing for weakly dependent errors and without assuming a specific type of endogeneity. Additionally, it suitably extends the asymptotic approach of Hall2007 to encompass cases with endogeneity, a common occurrence in practical applications. Consequently, this paper proposes two novel estimators and provides their detailed asymptotic properties under this broader setting. Notably, even in the case where there is no endogeneity and hence the FIVE reduces to the FLSE estimator of Park2012 (see Section (ref)), most of the asymptotic results that we obtain in the present paper have not been explored in the literature, to the best of the authors' knowledge. Given the potential prevalence of endogeneity in the functional linear model (see Section (ref)), we believe that the theoretical results presented in this paper hold value beyond the existing findings.} \phantomsection{\leavevmode\color{black}From a practical perspective, our methodology can be applied to study relationships between economic functional variables of which availability has been being increasing; potential examples include density-on-density regression model (see e.g., Park2012). }

\makeatletter \setcounter{page}{1} \setcounter{section}{0} \setcounter{rmk}{0} \makeatother

{Supplementary Material}

Preliminaries

Let $(\mathbb S,\mathbb F, \mathbb P)$ denote the underlying probability space and let $\mathcal H$ be a separable Hilbert space equipped with the inner product $\langle \cdot,\cdot\rangle$ and the usual Borel $\sigma$-field.

Random elements of Hilbert spaces

A $\mathcal H$-valued random variable $X$ is defined by a measurable map from $\mathbb S$ to $\mathcal H$. We say that such a random variable $X$ is integrable (resp.\ square-integrable) if $\mathbb{E} [\|X\|] < \infty$ (resp.\ $\mathbb{E}[ \|X\|^2] < \infty$), where $\|\cdot\|$ is the norm induced by the inner product. If $X$ is integrable, there exists a unique element $\mathbb{E}[X] \in \mathcal H$ satisfying $\mathbb{E}[\langle X,\zeta\rangle] = \langle \mathbb{E}[X],\zeta\rangle$ for every $\zeta \in \mathcal H$. The element $\mathbb{E}[X]$ is called the expectation of $X$.

Let $Y$ be another $\mathcal H$-valued random variable. We let $\otimes$ denote the tensor product defined as follows: for all $\zeta_1, \zeta_2\in \mathcal H$,

Note that $\zeta_1 \otimes \zeta_2$ is a linear map from $\mathcal H$ to $\mathcal H$. If $\mathbb{E}[\|X\|\|Y\|] < \infty$, we may well define a linear map $\mathcal C_{XY}$ from $\mathcal H$ to $\mathcal H$ as follows: $\mathcal C_{XY} = \mathbb{E}[(X-\mathbb{E}[X]) \otimes (Y-\mathbb{E}[Y])]$. $\mathcal C_{XY}$ is called the cross-covariance operator of $X$ and $Y$. If $X=Y$ and $X$ is square-integrable, we then may define $\mathcal C_{XX}$ similarly, and this is called the covariance operator of $X$. If the cross-covariance operator of two random variables $X$ and $Y$ are a nonzero operator, $X$ is said to be correlated with $Y$.

Bounded linear operators on $\mathcal H$

Let $\mathcal{L}_{\mathcal H}$ denote the space of bounded linear operators acting on $\mathcal H$, equipped with the operator norm $\|\mathcal T\|_{\mathrm{op}} = \sup_{ \|\zeta\|\leq 1} \|\mathcal T\zeta\|$. For any $\mathcal T \in \mathcal L_{\mathcal H}$, the adjoint of $\mathcal T$ (denoted $\mathcal T^\ast$) is the unique element of $\mathcal L_{\mathcal H}$ satisfying that $\langle \mathcal T \zeta_1 ,\zeta_2 \rangle=\langle \zeta_1 , \mathcal T^\ast \zeta_2 \rangle$ for all $\zeta_1,\zeta_2 \in \mathcal H$. {The range (denoted $\operatorname{ran} \mathcal T$) and kernel (denoted $\ker \mathcal T$) of $\mathcal T\in \mathcal{L}_{\mathcal H}$ are respectively given by $\operatorname{ran} \mathcal T = \{\mathcal T\zeta : \zeta \in \mathcal H\}$ and $\ker \mathcal T = \{\zeta \in \mathcal H : \mathcal T\zeta = 0\}$.} $\mathcal T$ is said to be nonnegative if $\langle \mathcal T \zeta,\zeta \rangle \geq 0$ for all $\zeta\in \mathcal H$, and positive if the inequality is strict. An element $\mathcal T \in \mathcal L_{\mathcal H}$ is called compact if there exist two orthonormal bases $\{\zeta_{1j}\}_{j \geq 1}$ and $\{\zeta_{2j}\}_{j \geq 1}$ of $\mathcal H$, and a sequence of real numbers $\{a_j\}_{j \geq 1}$ tending to zero, such that $\mathcal T = \sum_{j=1}^\infty a_j\zeta_{1j} \otimes \zeta_{2j} $. In this expression, it may be assumed that $\zeta_{1j} = \zeta_{2j}$ and $a_1 \geq a_2 \geq \ldots \geq 0$ if $\mathcal T$ is self-adjoint (i.e., $\mathcal T = \mathcal T^\ast$) and nonnegative Bosq2000. In this case, $a_j$ becomes an eigenvalue of $\mathcal T$ and hence $\zeta_{1j}$ is the corresponding eigenfunction, and moreover, we may define $\mathcal T^{1/2}$ by replacing $a_j$ with $\sqrt{a_j}$. It is well known that the covariance of a $\mathcal H$-valued random variable is self-adjoint, nonnegative and compact if exists. A linear operator $\mathcal T$ is called a Hilbert-Schmidt operator if its Hilbert-Schmidt norm $\|\mathcal T\|_{\operatorname{HS}} = (\sum_{j=1} ^\infty \Vert \mathcal T \zeta_j \Vert ^2)^{1/2}$ is finite, where $\{ \zeta_j \}_{j\geq 1}$ is an arbitrary orthonormal basis of $\mathcal H$. It is well known that $\|\mathcal T\|_{\mathrm{op}}\leq \|\mathcal T\|_{\operatorname{HS}}$ holds and the collection of Hilbert-Schmidt operators consists of a strict subclass of $\mathcal L_{\mathcal H}$; see Section 1.5 of Bosq2000.