The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
108,274 characters
Functional Linear Projection and Impulse Response Analysis
\title{\Large Functional Linear Projection and Impulse Response Analysis\thanks{We are grateful to \`{O}scar Jord\`{a}, Bonsoo Koo, Benjamin Wong, and seminar participants at the 18th International Symposium on Econometric Theory and Applications (SETA2024) and 33rd Australian New Zealand Econometric Study Group Meeting for their invaluable comments.} }
\author{ Won-Ki Seo\textsuperscript{a} \quad\quad Dakyung Seong\textsuperscript{a}\thanks{Corresponding author. Address: School of Economics, University of Sydney, Camperdown, 2006, NSW, Australia. E-mail addresses: \texttt{[email removed] (D.\ Seong), \texttt{[email removed] (W.-K.\ Seo)}}}\\ \large{\textsuperscript{a} School of Economics, University of Sydney}
}
\maketitle
\vspace{-1em}
\begin{abstract}
This paper proposes econometric methods for studying how economic variables respond to function-valued shocks. Our methods are developed based on linear projection estimation of predictive regression models with a function-valued predictor and other control variables. We show that the linear projection coefficient associated with the functional variable allows for the impulse response interpretation in a functional structural vector autoregressive model under a certain identification scheme, similar to well-known \citeauthor{Sims1972}' (\citeyear{Sims1972}) causal chain, but with nontrivial complications in our functional setup. A novel estimator based on an operator Schur complement is proposed and its asymptotic properties are studied. We illustrate its empirical applicability with two examples involving functional variables: economy sentiment distributions and functional monetary policy shocks. \\
\textbf{Keywords: }Local Projection, Functional linear regression, VARs, Identification
\end{abstract}
\section{Introduction}\label{sec: intro}
After the initial development by \cite{Oscar2005}, the \textit{local projection} approach has become one of the foremost applications of time series analysis in the fields of macroeconomics and public policy for studying dynamic causal effects, known as structural impulse responses. The work by \cite{PW2021} provides more rigorous grounding for the approach as an important complement to structural vector autoregressive (SVAR) models, by establishing the theoretical equivalence between estimation results of the two. Practitioners have benefited from this result, as the local projection approach is easier to implement and interpret, computationally less burdensome, and less sensitive to model misspecification (see \citealp{Oscar2005} for a detailed discussion).
This paper contributes to recent developments in this area by providing statistical inference methodologies to study responses of target variables when shocks are characterized by function-valued random variables, such as the functional monetary policy shocks in \cite{IR2021}. We consider a linear projection of a target variable onto the space spanned by the variables in the conditioning set, including function-valued shocks. Similar to \cite{Oscar2005}, this setup simplifies impulse response analysis into the estimation of a regression model. Our first result is the formal theoretical evidence for interpreting the parameters of interest in our regression model as structural impulse response functions of SVAR models involving a functional variable, under a certain identification condition. Considering the scarcity of the literature, our paper will be a valuable first step toward formally understanding functional structural shocks and their impact. Furthermore, while the equivalence for this special case has been formally established by \cite{Oscar2005} for finite dimensional VAR models, its extension to the case with a functional covariate has not yet been fully explored.
Our benchmark model to be studied is closely related to the \textit{functional local projection} approach proposed by \cite{IR2021}, but a crucial distinction exists. Specifically, in contrast to them, we do not require a parametric assumption on the structure of function-valued shocks and their coefficients. In \cite{IR2021}, functional shocks are represented by a few known functional factors. This assumption crucially reduces the dimensionality of the variables to be analyzed, thereby making existing estimation and inference methodologies valid despite the infinite dimensionality of functional data. However, such a parametric structure is not always available in practice and could result in model misspecification; this point will be further demonstrated in Sections~\ref{sec: model} and~\ref{subsec: para} by using an example. Considering practical challenges of exact identification of those factors in infinite dimensional function spaces, our approach could be an appealing alternative for practitioners.
The current paper is closely related to recently growing studies on the VAR model involving a functional covariate, including \cite{IR2021}, \cite{BHCJ2023}, and \cite{Chang_Chen_Schorfheide_2024}. However, our focus on identification and the relationship between SVAR models and linear projection in the presence of a functional covariate distinguishes the current paper from theirs. Moreover, the aforementioned studies typically approximate a functional covariate using a few factors or basis functions and then apply identification strategies developed in the finite dimensional SVAR literature (e.g.,\ \citealp{kilian2017structural}), assuming that the approximation error is either zero or negligible as the sample size increases. Practitioners might prefer these approaches due to their ease of implementation. However, this simplified method does not always ensure the identification of infinite dimensional structural parameters on the entire space in which the parameters take values, thereby restricting our ability to fully comprehend the features of structural parameters. This point will be detailed in Section~\ref{subsec: para} by using an example.
We propose an estimator based on the operator Schur complement and establish its asymptotic properties, including (local) asymptotic normality. As our model belongs to the so-called scalar-on-function models studied in, e.g., \cite{Hall2007}, \cite{SHIN2009}, \cite{Florence2015}, \cite{imaizumi2018}, and \cite{Babii2022}, our estimator can be understood as a complement to the existing ones developed therein. However, existing estimators (i) first reduce the dimension of functional predictors and then use the dimension-reduced ones in the analysis (e.g., \citealp{AP2006,SHIN2009}), (ii) are designed without scalar-valued covariates (e.g., \citealp{Hall2007, Florence2015, imaizumi2018}), or (iii) do not account for time series dependence, which is crucial in our context (except for \citealp{Babii2022}, among the aforementioned articles). The approaches in (i) may not be preferred because they do not take into account covariation across variables in the dimension reduction procedure, while the latter two may limit practical applicability. Our estimator complements them by using a regularization method that allows us to account for covariation between scalar- and function-valued variables. This point is particularly relevant if the function-valued variable is given by a regressor without a structural interpretation, as in Section~\ref{sec: svar}. Later in Appendix~\ref{sec:est2}, we further extend our method to accommodate an endogenous functional predictor.
As discussed in \cite{Hall2007}, \cite{Florence2015}, and references therein, estimating functional linear regression models intrinsically involves an inverse problem similar to those in nonparametric estimation. Although we focus on a model that is linear in functional variables, these function-valued random variables can be represented by an infinite number of basis functions. In this regard, our paper is also related to studies on nonlinear impulse responses in, e.g., \cite{Oscar2005}, \cite{kilian2017structural}, and \cite{KP2024}.
Although our main focus is on estimating impulse response coefficients when shocks are characterized by functions, our methodology could also be an empirically appealing alternative to prediction methods, such as that in \cite{Barbaglia2023}, in that it allows us to utilize richer information, such as distributional information or observations at different frequencies. In particular, if a functional variable is given by a predictor without a structural interpretation, our model reduces to a predictive regression model with a functional predictor. As studied in \cite{Babii2022} and \cite{seong2021functional}, functional predictors could lead to better forecasting outcomes by exploiting information overlooked in conventional estimation approaches.
In our empirical application, we first employ the quantile curve of the economic sentiment measure proposed by \cite{Barbaglia2023} and study the impact of perturbations given to the sentiment distribution on Total Nonfarm Payrolls in the US, without a structural interpretation. The sentiment measure is observed at a daily frequency and describes the presence of negative or positive tones in relation to specific terms, such as \textsl{economy} or \textsl{inflation}. The proposed estimator shows that the response depends not only on the magnitude of the perturbation, but also on its shape, which has not been observed previously. Moreover, we find that location shifts in the sentiment distribution have significant effects on the target variable in the short term. For instance, if the overall sentiment distribution shifts to the left, the average sentiment on the economy becomes more pessimistic, predicting significantly negative economic growth in the near future.
We then revisit the work by \cite{IR2021} on the impact of monetary policy shocks, which are characterized by shifts in yield curves on monetary policy announcement dates. We find that during conventional periods, inflation responses are generally insignificant; however, their impact tends to vary depending on how the shock affects interest rates at short or long maturities. That is, the impact of functional monetary policy shocks depends on their shapes and directions.
The rest of the paper is organized as follows. Section~\ref{sec: model} motivates our benchmark model and provides examples of functional variables. In Section~\ref{sec: svar}, we consider the SVAR model involving a functional variable and discuss when the parameters in our model can be interpreted as functional structural impulse responses. Section~\ref{sec:est} proposes our estimator and establishes its asymptotic properties. We apply our estimators to the empirical data in Section~\ref{sec:emp}. Section~\ref {sec:sim} summarizes simulation results. Section~\ref{sec:con} concludes. The appendix includes proofs and an extension of the main theoretical results to instrumental variable (IV) estimation. In the Supplementary Material, we provide mathematical preliminaries and extension to functional SVAR models.
\section{Model} \label{sec: model}
We let $y_t$ denote a scalar-valued target (dependent) variable and $\mathbf{w}_t= ( {w}_{1,t},\ldots, {w}_{m,t})'$ be a vector of exogenous (scalar-valued) control variables possibly containing lagged $y_t$'s. The variable $X_t = \{X_t (r) : r \in [0,1]\} $ denotes a \textit{function-valued} variable that takes values in some function space $\mathcal H$ to be detailed shortly.\footnote{For the ease of explanation, we suppose that the domain of $X_t$ is $[0,1]$, but its extension to any arbitrary compact interval $[a,b]$ is straightforward.} Examples of $X_t$ include, but are not limited to, sentiment quantile curves in Example~\ref{example1} and functional monetary policy shocks in \cite{IR2021}. In this paper, we are particularly interested in the response of $y_{t+h}$ when an additional perturbation (or shock) $\zeta$ is introduced to $X_t$; that is, for any $x \in \mathcal H$, we are interested in \begin{equation}
\mathbb E\left[y_{t+h} |X_t = x+ \zeta , \mathbf{w}_{t} \right] - \mathbb E\left[y_{t+h} | X_t = x, \mathbf{w}_{t}\right]. \label{eq: irf: general}
\end{equation}
In previous studies on time series analysis, including \cite{Oscar2005} and \cite{PW2021}, $X_t$, the variable exposed to a perturbation, is characterized by a scalar or a vector, and the linear projection of $y_t$ onto the space spanned by $\{X_t, \mathbf w_t\}$ reduces the impulse response analysis to inference of projection coefficients. However, its extension to the case where $X_t$ is given by a function is not straightforward. For example, in our setup, the perturbation $\zeta$ in \hyperref[{eq: irf: general}]{\textup{\tagform@{\ref*{eq: irf: general}}}} is specified as a function on $[0,1]$. Thus, the response of $y_t$ in \hyperref[{eq: irf: general}]{\textup{\tagform@{\ref*{eq: irf: general}}}} depends not only on the magnitude of $\zeta$, but also on its shape. This feature has not been observed in the conventional studies on SVAR models.
We follow \cite{Oscar2005} and consider the following linear projection of the scalar target variable $y_t$ onto the space spanned by the variables in the conditioning set, including the function-valued $X_t$: \begin{equation}
y_{t+h} = \int_{0}^{1} \beta_h(r)X_t(r) dr + \mathbf{w}_{t} ' \alpha_h
+ u_{h,t}. \label{eq: model: benchmark}
\end{equation}
The functional covariate $X_t$ and its coefficient $\beta_h$ are assumed to take values in $\mathcal H=L^2[0,1]$ (the Hilbert space of square-integrable functions), which allows us to conveniently write $\int_{0}^{1} \beta_h(r)X_t(r) dr$ as the inner product $\langle \beta_h, X_t \rangle$ defined on $\mathcal H$.
Under the linear projection in \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}}, we identify \hyperref[{eq: irf: general}]{\textup{\tagform@{\ref*{eq: irf: general}}}} as follows:\begin{equation}
\mathbb E\left[y_{t+h} |X_t = x+ \zeta , \mathbf{w}_{t} \right] - \mathbb E\left[y_t | X_t = x, \mathbf{w}_{t}\right] = \int_{0}^{1} \beta_h(r)\zeta(r)dr = \langle \beta_h,\zeta\rangle. \label{def: local: irf}
\end{equation}
If $X_t$ is given by a functional structural shock (such as the functional monetary policy shock in \citealp{IR2021}) or if the true data generating process (DGP) follows a special structure in Section~\ref{sec: svar}, the parameter $\beta_h$ reduces to the \textit{structural impulse response function} (SIRF) of $y_t$ when $X_t$ experiences a function-valued shock $\zeta$. In this regard, we extend the seminal work of \cite{Oscar2005} to allow for a functional variable.
\begin{figure}[h!]
\centering
\caption{Quantile Curves of Economy Sentiment}
\begin{subfigure}{.6\textwidth}\subcaption{Overall sentiment}
\includegraphics[width = \textwidth, height = 0.5\textwidth, trim= {0 0 0 5cm},clip]{figure/3dsentiment.png}
\label{fig0a}
\end{subfigure}\\
\begin{subfigure}{.32\textwidth}\subcaption{Overall period}
\includegraphics[width = \textwidth]{figure/sentiment.jpg}
\label{fig1a}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\subcaption{Recessive period}
\includegraphics[width = \textwidth]{figure/sentiment_rec.jpg}
\label{fig1b}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\subcaption{Expansive period}
\includegraphics[width = \textwidth]{figure/sentiment_exp.jpg}
\label{fig1c}
\end{subfigure}\vspace{-1.5em}
\flushleft{\scriptsize{Notes: The monthly quantile curves of economic sentiment are reported for overall, recessive, and expansive periods with each other their mean functions (black). Each curve is obtained by smoothing the daily economic sentiment measure proposed by \cite{Barbaglia2023} for each month, see Section~\ref{sec:emp1} for details. }}
\label{fig: 1}
\end{figure}
\begin{figure}[h!]
\caption{Examples of $\zeta$ in \hyperref[{def: local: irf}]{\textup{\tagform@{\ref*{def: local: irf}}}} and distributional changes}
\begin{subfigure}{.5\textwidth} \flushleft{\subcaption{ positive location shift ($\zeta_1$)\label{fig2a}}}
\includegraphics[width = 0.49\linewidth]{figure/Shock_quant2.jpeg}
\includegraphics[width = 0.49\linewidth]{figure/Shock_dist2.jpeg}
\end{subfigure}
\begin{subfigure}{.5\textwidth} \subcaption{ negative location shift ($\zeta_2$)\label{fig2b}}
\includegraphics[width = 0.49\linewidth]{figure/Shock_quant3.jpeg}
\includegraphics[width = 0.49\linewidth]{figure/Shock_dist3.jpeg}
\end{subfigure}
\begin{subfigure}{.5\textwidth} \subcaption{ greater dispersion ($\zeta_3$) \label{fig2c}}
\includegraphics[width = 0.49\linewidth]{figure/Shock_quant4.jpeg}
\includegraphics[width = 0.49\linewidth]{figure/Shock_dist4.jpeg}
\end{subfigure}
\begin{subfigure}{.5\textwidth} \subcaption{ greater dispersion, negative location shift ($\zeta_4$)\label{fig2d}}
\includegraphics[width = 0.49\linewidth]{figure/Shock_quant5.jpeg}
\includegraphics[width = 0.49\linewidth]{figure/Shock_dist5.jpeg}
\end{subfigure}
\label{fig: 2}\vspace{-.5em}
\flushleft{\scriptsize{Notes: Each panel reports the impact of the perturbation $\zeta$ on the average sentiment quantile (left) and its associated distribution (right). The grey (resp.\ black) lines report the average quantile and its associated probability density function before (resp.\ after) the perturbation is given. }}
\end{figure}
\begin{example}[Economic Sentiment Quantile Functions]\label{example1}
In Section~\ref{sec:emp1}, we consider a quantile function as an example of $X_t$, which represents the distribution of \citepos{Barbaglia2023} sentiment measure about the US economy. The distributional information is expected to provide meaningful insights into how economic variables respond when individuals' sentiment about the economy changes.
To help with intuition, we report the sentiment curves in Figure~\ref{fig: 1}. Figure~\ref{fig: 1}.\ref{fig0a} reports the sentiment quantile data over the sample period from January 1984 to December 2021. Figures~\ref{fig: 1}.\ref{fig1a}, \ref{fig: 1}.\ref{fig1b}, and \ref{fig: 1}.\ref{fig1c} respectively report the quantiles during the overall, recessive, and expansive periods with the probability $r\in[0,1]$ on the x-axis. The black lines in the figures represent the mean functions for each period.\footnote{The recessive and expansive periods are defined following the NBER Business Cycle dates: \url{https://www.nber.org/research/business-cycle-dating}.} The first interesting observation is the different scales of the vertical axis; during recessive periods, individuals tend to hold negative sentiments on the economy and its mean function exhibits a larger slope ($\approx$ -1.96), compared to overall ($\approx$ -0.002) and expansive ($\approx$ 0.2) periods. This suggests that the economy, on average, tends to hold a higher level of heterogeneous beliefs, which may be explained by greater uncertainty during that time. Such distributional information will be lost if we aggregate $X_t$ into a single scalar-valued random variable.
When a random shock is given to the sentiment distribution, its impact on economic variables may depend on whether the shock increases or decreases overall economic sentiment and/or changes the level of disagreement about the economic sentiment. This can be studied by, for example, setting $\zeta (\equiv \zeta(\cdot))$ in \hyperref[{def: local: irf}]{\textup{\tagform@{\ref*{def: local: irf}}}} to $\zeta_1(s) = 2 $, $\zeta_2(s) = -2 $, $\zeta_3(s) = 2 s$, and $\zeta_4(s) = -2 + 2s$, for $s \in[0,1]$.
Figures~\ref{fig2a}--\ref{fig2d} present the average sentiment quantiles and distributions before (grey) and after (black) these perturbations are introduced. As shown in Figures~\ref{fig2a} and \ref{fig2b}, the mean function, tends to shift up or down in its level when $\zeta_1$ and $\zeta_2$ are introduced, which corresponds to the location shifts of the sentiment distribution. Meanwhile, $\zeta_3$ increases the slope of the quantile without changing its average level. Thus, the sentiment distribution exposed to the perturbation will have a greater dispersion, meaning a higher level of disagreements about economy status. The last shock $\zeta_4$ increases the level of disagreement and negatively shifts the overall sentiment distribution, so that it gets similar to the average quantile observed in recessive periods.
\end{example}
\begin{figure}[h!]
\centering
\caption{ Functional Monetary Policy Shocks in \cite{IR2021} }
\begin{subfigure}{.32\textwidth}\subcaption{Overall period}
\includegraphics[width = \textwidth]{figure/IR2/shock_overall.jpg} \label{fig1a:ir}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\subcaption{Conventional period}
\includegraphics[width = \textwidth]{figure/IR2/shock_conv.jpg} \label{fig1b:ir}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\subcaption{Unconventional period}
\includegraphics[width = \textwidth]{figure/IR2/shock_unconv.jpg} \label{fig1c:ir}
\end{subfigure}\vspace{-1.5em}
\flushleft{\scriptsize{Notes: The monthly functional monetary policy shocks proposed by \cite{IR2021} are reported for overall, conventional, and unconventional periods for the period from January 1995 to June 2016. The black lines represent the mean functions for each of these periods. The shocks are characterized by the shifts in yield curves on monetary policy announcement dates, see \cite{IR2021} for details. }}
\label{fig: ir: 1}
\end{figure}
\begin{example}[Functional Monetary Policy Shocks]\label{example2} Another important example of $X_t$ is functional monetary policy shocks proposed by \cite{IR2021}.
The functional monetary policy shocks are reported in Figure~\ref{fig: ir: 1} for the overall (Fig.\ \ref{fig1a:ir}), conventional (Fig.\ \ref{fig1b:ir}), and unconventional (Fig.\ \ref{fig1c:ir}) periods, along with their mean functions (black). The functional monetary policy shocks in conventional periods appear to have regular shapes. Specifically, during conventional periods, the shocks tend to have a greater negative value at short maturities, compared to those at long maturities. Such a regular shape is not observed during unconventional periods in Figure~\ref{fig1c:ir}. The figures suggest that monetary policy shocks realized by shifts in yield curves may allow us to utilize richer information related to their shape and direction at different maturities. This example was formerly studied by \cite{IR2021} and will be further explored in Section~\ref{sec:emp2}.
\end{example}
\begin{figure}[h!]
\centering
\caption{Financial Risk Shocks }
\begin{subfigure}{.32\textwidth}\subcaption{Overall period}
\includegraphics[width = \textwidth]{figure/ted_general.jpg}
\label{fig1a:ted}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\subcaption{Recessive period}
\includegraphics[width = \textwidth]{figure/ted_rec.jpg}
\label{fig1b:ted}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\subcaption{Expansive period}
\includegraphics[width = \textwidth]{figure/ted_exp.jpg}
\label{fig1c:ted}
\end{subfigure} \label{fig: ted} \vspace{-1.5em}
\flushleft{\scriptsize{Notes: Each figure reports the monthly quantiles of TED spread for the period from February 1990 to December 2022. The black lines report the mean function in the overall, recessive, and expansive periods. }}
\end{figure}
\begin{example}[Proxy of Financial Risk Shocks]
Another example is the monthly quantiles of the TED spread that captures the likelihood of unpredictable default risk in financial markets. The TED spread is popularly employed in the literature as a proxy or an IV to measure the structural shock associated with financial markets (see \citealt[p.\ 110]{SW2012}). This quantile function could be used in place of $X_t$ in \hyperref[{eq: irf: general}]{\textup{\tagform@{\ref*{eq: irf: general}}}} as a proxy for the financial or liquidity shock, or as a function-valued instrument for the IV estimator developed in Appendix~\ref{sec:est2}.
\end{example}
\section{Relation to SVAR with a functional covariate} \label{sec: svar}
If $X_t$ is identified as a functional structural shock, as in \cite{IR2021}, the parameter $\beta_h$ in \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} can be certainly interpreted as the SIRF of $y_t$ to that shock. However, its general
link to the SIRF similar to those in \cite{Oscar2005} or \cite{PW2021} has not yet formally established in the functional setup, although \hyperref[{eq: irf: general}]{\textup{\tagform@{\ref*{eq: irf: general}}}} and \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} partly provide an intuitive guidance on its interpretation as a direct estimator of the reduced-form impulse response.
In this section, we show that the parameter $\beta_h$ in \hyperref[{def: local: irf}]{\textup{\tagform@{\ref*{def: local: irf}}}} provides a crucial and interesting interpretation as a SIRF when the true DGP follows an autoregressive structure. This suggests the importance of our benchmark model \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}}, particularly as an extension of \cite{Oscar2005}.
\subsection{Formulation of the functional SVAR model}
For the ease of exposition, we exclude $\mathbf w_t$ in this section. The results given in this section can be easily extended to include other covariates as long as they are finite dimensional.
Thus we do not lose any generality by this simplification.
As mathematical preliminaries, Section~\ref{sec_prelim} of the Supplementary Material provides a brief review of definitions and properties of \(\mathcal{H}\)-valued random variable \(X_t\), linear operators on \(\mathcal{H}\), and the product Hilbert space \(\mathbb{R} \times \mathcal{H}\), where the tuple \(\left[\begin{smallmatrix} y_t\\ X_t\end{smallmatrix}\right]\) takes values.\footnote{We treat the tuple \( \left[\begin{smallmatrix} y\\ x\end{smallmatrix} \right] \), consisting of \(\mathbb{R}\)-valued and \(\mathcal{H}\)-valued random variables, as a “column vector” in the multivariate setting. This notation introduces no confusion or mathematical imprecision in the subsequent discussion. Similar notation will be used for any other tuples appearing later.}
As detailed in the appendix, a linear operator \(\mathcal{D}\) on \(\mathbb{R} \times \mathcal{H}\) can be represented as an operator matrix, say $\mathcal{D} = \left[\begin{smallmatrix} d_{11} & d_{12} \\ d_{21} & d_{22} \end{smallmatrix}\right].$
This operator maps $ \left[\begin{smallmatrix} y\\ x\end{smallmatrix} \right] \in \mathbb{R} \times \mathcal{H} $ to
$\left[ \begin{smallmatrix}
d_{11}y + d_{12}x\\ d_{21}y + d_{22}x
\end{smallmatrix}\right] \in \mathbb{R} \times \mathcal{H}$,
analogous to transforming a \((2 \times 1)\) vector using a \((2 \times 2)\) matrix. For this reason, with a slight abuse of notation, we henceforth treat the \((\mathbb{R} \times \mathcal{H})\)-valued random element \( \left[\begin{smallmatrix} y\\ x\end{smallmatrix} \right]\) as if it were a \((2 \times 1)\) random vector, and any linear operator on \(\mathbb{R} \times \mathcal{H}\) as if it were a \((2 \times 2)\) matrix. This simplification is valid in the considered setup without any loss of technical rigor (see Section~\ref{sec_prelim2} of the Supplementary Material).
In this section, we suppose that the tuple $\left[\begin{smallmatrix} y_t\\ X_t\end{smallmatrix}\right] \in \mathbb R \times \mathcal H$ follows the model below:
\begin{equation}
\underbrace{\begin{bmatrix} I_1& \beta_{12} \\ \beta_{21}& I_2\end{bmatrix} }_{:=\mathcal B }\begin{bmatrix} y_t
\\ X_t \end{bmatrix} =\underbrace{ \begin{bmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & a_{22} \end{bmatrix}}_{:=\mathcal A}\begin{bmatrix} y_{t-1}
\\ X_{t-1} \end{bmatrix} + \begin{bmatrix} u_{1t}
\\ U_{2t}\end{bmatrix}.\label{eq: model: svar}
\end{equation}
The above model is identical to the standard bivariate SVAR model other than that $X_t$ is a $\mathcal H$-valued random variable, and thus the coefficients $\mathcal A$ and $\mathcal B$ are in fact operator matrices.
Specifically, $\alpha_{12}$ and $ \beta_{12}$ are given by linear maps from $ \mathcal H$ to $ \mathbb R$, while $\alpha_{21}$ and $\beta_{21}$ are maps from $\mathbb R$ to $\mathcal H$. The $(1,1)$-th (resp.\ $(2,2)$-th) elements, $I_1$ and $a_{11}$ (resp.\ $I_2$ and $a_{22}$), are linear operators acting on $\mathbb{R}$ (resp.\ $\mathcal H)$, where $I_1$ and $I_2$ denote the identity maps in the relevant spaces. The tuple of structural shocks $\left[\begin{smallmatrix} u_{1t}\\U_{2t}\end{smallmatrix}\right]$ is a $(\mathbb R \times \mathcal H)$-valued random element. Its covariance operator (see \hyperref[{eqcovop}]{\textup{\tagform@{\ref*{eqcovop}}}} in the Supplementary Material) is given as follows:
\begin{equation*}
\mathbf \Sigma
= \begin{bmatrix} \mathbb{E}[u_{1t} \otimes u_{1t}] & \mathbb{E}[U_{2t} \otimes u_{1t}] \\ \mathbb{E}[u_{1t} \otimes U_{2t}] & \mathbb{E}[U_{2t}\otimes U_{2t}]\end{bmatrix} = \begin{bmatrix} \sigma_{11} & 0 \\ 0 & \mathbf{\Sigma}_{22} \end{bmatrix},
\end{equation*}
where $\otimes$ denotes the tensor product, which generalizes the outer product in the Euclidean space, and $\sigma_{11}$ (resp.\ $\mathbf{\Sigma}_{22}$) is a linear operator acting on $\mathbb{R}$ (resp.\ $\mathcal H$).\footnote{Since $\mathbf{\Sigma}$ is an operator matrix whose $(i,j)$-th entry is a map from $\mathcal{H}_i$ to $\mathcal{H}_j$, with $\mathcal{H}_1 = \mathbb R$ and $\mathcal{H}_2=\mathcal{H}$, $\sigma_{11}$ is given by a linear map on $\mathbb{R}$. However, any linear map on $\mathbb{R}$ is nothing but a scalar multiplication, given by $c I_1$ for some $c \in \mathbb{R}$. Thus, there is little risk of confusion even if we understand $\sigma_{11}$ as a real-valued constant. \label{foot1}
} See Section~\ref{sec_prelim} of the Supplementary Material for details.
We first establish the fully functional identification of structural parameters and highlight its differences from existing approaches. Then, we link the SIRF implied by the SVAR structure to the coefficients in our benchmark model.
\subsection{Fully Functional Identification of the SVAR with a functional covariate} \label{sec_svar}
It can be shown that the operator $\mathcal B$ in \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} is invertible and $\Gamma:=\mathcal B^{-1}\mathcal A$ can be well defined under the condition in Proposition \ref{prop: svar: identification} that will appear shortly. Thus \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} can be written into the following reduced-form VAR (RFVAR) model:
\begin{equation}
\begin{bmatrix} y_t
\\ X_t \end{bmatrix} =\underbrace{\begin{bmatrix} \gamma_{11} & \gamma_{12} \\ \gamma_{21} & \gamma_{22} \end{bmatrix}}_{:=\Gamma}\begin{bmatrix} y_{t-1}
\\ X_{t-1} \end{bmatrix} + \begin{bmatrix} \varepsilon_{1t}
\\ \mathcal E_{2t}\end{bmatrix}, \quad\text{ where }\quad \Sigma_{\varepsilon} = \begin{bmatrix} \sigma_{\varepsilon,11} & \Sigma_{\varepsilon,12} \\ \Sigma_{\varepsilon,21} & \Sigma_{\varepsilon,22} \end{bmatrix}.\label{eqreduced1}
\end{equation}
$\Sigma_{\varepsilon}$ denotes the covariance of $\left[\begin{smallmatrix}\varepsilon_{1t}\\\mathcal E_{2t}\end{smallmatrix}\right]$.
In this section, we discuss conditions under which the structural parameters in \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} can be identified from the reduced-form parameters in \hyperref[{eqreduced1}]{\textup{\tagform@{\ref*{eqreduced1}}}}.
Our assumption for achieving such identification is based on the standard causal ordering of the (contemporaneous) variables in the system, as in \cite{Sims1972,Sims1980}. Later, we will show that the coefficients of our benchmark model reduces to the SIRFs under this identification scheme.
\begin{proposition}\label{prop: svar: identification}\normalfont
Suppose that $ \left[\begin{smallmatrix} y_t\\ X_t\end{smallmatrix} \right]$ satisfies \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}}. Then the structural parameters in \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} are identified
under either
\begin{enumerate*}[(i)]
\item\label{prop: svar: identification1} \textsl{$\beta_{12} = 0$} or
\item\label{prop: svar: identification2} \textsl{$\beta_{21}=0$ and $\mathbf\Sigma_{22}$ is injective.}
\end{enumerate*}
\end{proposition}
Although the conditions \ref{prop: svar: identification1} and \ref{prop: svar: identification2} in Proposition \ref{prop: svar: identification} are seemingly similar to each other, they suggest different levels of complexity in estimating the outcome equations. Specifically, if there is no contemporaneous impact of function-valued $X_t$ to scalar-valued $y_t$ (i.e., $\beta_{12} = 0$), the outcome equation of $y_t$ reduces to $y_t = \alpha_{11} y_{t-1} + \alpha_{12}X_{t-1} + u_{1t}$, and thus it includes only one infinite dimensional parameter $\alpha_{12}$. On the other hand, under the second identification condition (i.e., $\beta_{21} = 0$), the outcome equation reduces to $y_t = \alpha_{11} y_{t-1} -\beta_{12}X_t + \alpha_{12}X_{t-1} + u_{1t}$, in which the number of infinite dimensional parameters doubles. Estimating these parameters involves solving an inverse problem that necessitates the use of a regularization scheme. Considering this, the resulting estimator from the second identification scheme is likely to suffer from a larger regularization bias.
Another interesting implication of Proposition~\ref{prop: svar: identification} is that the standard rank-based identification strategy cannot be straightforwardly translated into our setup involving a functional covariate. Instead, depending on the direction of the contemporaneous impact, we may additionally need an injectivity condition
to identify structural parameters from the reduced-form parameters.
For example, under the condition \ref{prop: svar: identification2} in Proposition~\ref{prop: svar: identification}, $\Sigma_{\varepsilon,22} = \Sigma_{22}$ and $\beta_{12}$ satisfies
\begin{equation} \label{eqinject}
\sigma_{e,12} =\beta_{12}\Sigma_{22};
\end{equation}
see our proof of Proposition~\ref{prop: svar: identification} in Appendix~\ref{app_proof}. In this case, the injectivity of $\Sigma_{22}$ is required to uniquely identify $\beta_{12}$ from \hyperref[{eqinject}]{\textup{\tagform@{\ref*{eqinject}}}} (see \citealp{carrasco2007linear,Seo2024}).
The two observations mentioned above contrast with identification conditions and their implications in the existing literature, such as \cite{Sims1972}, where a different ordering of variables primarily affects their economic interpretation and the direction of contemporaneous shocks, rather than the level of asymptotic bias or computational burden. Thus, in the SVAR model with functional covariates, identification conditions need to be cautiously chosen with taking into account the above.
\begin{remark}\label{remchol}
In the existing literature on SVAR models concerning $k$-dimensional vector-valued time series, it is commonly assumed that the $k\times k$ matrix $\Sigma_\varepsilon$ allows the Cholesky factorization such that $\Sigma_\varepsilon = LL'$ for a lower triangular matrix $L
$. Then $\mathcal B $ in \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} simply reduces to $ L^{-1}$. This is one of the most popular identification strategies and similar to the causal ordering in our Proposition~\ref{prop: svar: identification}. However, in the considered setup, where $X_t$ is a function-valued random element, the covariance $\Sigma_\varepsilon$ is not invertible on $\mathbb{R}\times\mathcal H$ (see \citealp{Mas2007}, Section 2.2), which makes it infeasible to apply the existing identification strategy.
\end{remark}
\subsection{SIRF and its relationship with the linear projection} \label{sec_IRF}
To understand the relationship between the coefficient in our benchmark model \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} and the SIRF implied by \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}}, it is convenient to consider the MA($\infty$) representation.
Under the invertibility of $\mathcal B$, the SVAR model in \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} allows the reduced-form representation in \hyperref[{eqreduced1}]{\textup{\tagform@{\ref*{eqreduced1}}}} and furthermore, we have
\begin{equation} \label{eqrfvar}
\begin{bmatrix}y_t\\ X_t
\end{bmatrix} = \sum_{j=0}^\infty \Gamma^{j} \mathcal B^{-1} {\mathbf{u}}_{t-j},
\end{equation}
where $\mathbf{u}_t = \left[\begin{smallmatrix}{u}_{1,t}\\{U}_{2,t}\end{smallmatrix}\right] \in\mathbb R \times \mathcal H$ and $\Gamma ^0 = I$ which is the identity operator on $\mathbb{R}\times \mathcal H$. \hyperref[{eq: irf: general}]{\textup{\tagform@{\ref*{eq: irf: general}}}} and \hyperref[{eqrfvar}]{\textup{\tagform@{\ref*{eqrfvar}}}} suggest that the SIRF of $ \left[\begin{smallmatrix} y_t\\ X_t\end{smallmatrix} \right]$ at horizon~$h$, when the perturbation $\widetilde\zeta \equiv \left[\begin{smallmatrix}1 \\ \zeta\end{smallmatrix}\right] \in \mathbb{R}\times \mathcal H$ is introduced to $\mathbf{u}_t$, is given as follows:
\begin{equation}\label{eq: sirf}
\operatorname{SIRF}_h(\widetilde{\zeta}) :=\mathbb{E}\left[\left.\begin{bmatrix} y_{t+h} \\ X_{t+h} \end{bmatrix} \right| {\mathbf u}_t+\widetilde\zeta,{\mathbf u}_{t-1},\ldots\right] - \mathbb{E}\left[ \left.\begin{bmatrix} y_{t+h} \\ X_{t+h} \end{bmatrix} \right| \mathbf u_t,\mathbf u_{t-1},\ldots\right] = \Gamma^h\mathcal B^{-1} \widetilde{\zeta}.
\end{equation}
Both $\Gamma ^h$ and $\mathcal B^{-1}$ are operator matrices mapping from $\mathbb R \times \mathcal H$ to $\mathbb R \times \mathcal H$, and the same holds for $\Gamma ^h\mathcal B^{-1}$. Therefore, $\Gamma ^h \mathcal B^{-1}$ allows the following representation:
\begin{equation} \nonumber
\Gamma^{h} \mathcal B^{-1} = \begin{bmatrix} P_{\mathbb{R}}\Gamma^{h} \mathcal B^{-1}P_{\mathbb{R}}^\ast & P_{\mathbb{R}} \Gamma^{h} \mathcal B^{-1}P_{\mathcal{H}}^\ast \\ P_{\mathcal{H}}\Gamma^{h} \mathcal B^{-1}P_{\mathbb{R}}^\ast & P_{\mathcal{H}}\Gamma^{h} \mathcal B^{-1}P_{\mathcal{H}}^\ast \end{bmatrix} = \begin{bmatrix} \operatorname{SIRF}_{11,h} & \operatorname{SIRF}_{12,h} \\ \operatorname{SIRF}_{21,h} & \operatorname{SIRF}_{22,h} \end{bmatrix},
\end{equation}
where $P_{\mathbb{R}}$ denotes the projection map given by $P_{\mathbb{R}}\left[\begin{smallmatrix} x_1\\ x_2\end{smallmatrix}\right] = x_1$ and $P_{\mathbb{R}}^\ast$ is its adjoint satisfying $P_{\mathbb{R}}^\ast(x_1) = \left[\begin{smallmatrix} x_1\\ 0\end{smallmatrix}\right]$. Similarly, $P_{\mathcal{H}}$ (resp.\ $P_{\mathcal{H}}^\ast$) is a projection map satisfying $P_{\mathcal{H}} \left[\begin{smallmatrix}
x_1\\x_2
\end{smallmatrix}\right] = x_2$ (resp.\ $P_{\mathcal{H}}^\ast (x_2) =\left[\begin{smallmatrix}
0\\x_2
\end{smallmatrix}\right] $) for $x_2 \in \mathcal H$. Given these projection maps, each element of $\Gamma ^h \mathcal B^{-1}$ is characterized by different linear maps\footnote{$\operatorname{SIRF}_{11,h} : \mathbb R \to \mathbb R$ (i.e., a scalar multiplication), $\operatorname{SIRF}_{12,h}: \mathcal H \to \mathbb R$, $\operatorname{SIRF}_{21,h}: \mathbb R \to \mathcal H$ and $\operatorname{SIRF}_{22,h}:\mathcal H \to \mathcal H$.} and measures an effect of an additional structural shock on $y_{t+h}$ or $X_{t+h}$. For example, the response of $y_{t+h}$ to a function-valued shock $\zeta\in \mathcal H$ is given by a scalar such that \begin{equation}\label{eqpartialirf1}
\operatorname{SIRF}_{12,h}(\zeta) = P_{\mathbb{R}} \Gamma^{h} \mathcal B^{-1}P_{\mathcal{H}}^\ast (\zeta) = \delta_{12,h} \in \mathbb R.
\end{equation}
The response of $X_{t+h}$ to the unit shock on $u_{1,t}$ is characterized by a function such that
\begin{equation}\label{eqpartialirf2}
\operatorname{SIRF}_{21,h}(1)=P_{\mathcal{H}} \Gamma^{h} \mathcal B^{-1} P_{\mathbb{R}}^\ast (1) = \delta_{21,h} \in \mathcal H.
\end{equation}
These are often of interest to practitioners.
The SIRFs in \hyperref[{eqpartialirf1}]{\textup{\tagform@{\ref*{eqpartialirf1}}}} and \hyperref[{eqpartialirf2}]{\textup{\tagform@{\ref*{eqpartialirf2}}}} can be interpreted similarly to those of standard SVAR models. However, because $X_t$ and $U_{2,t}$ are given by $\mathcal H$-valued functions of infinite dimension, the implementation of $\zeta$ to the structural error can be represented with an infinite number of basis functions. That is, for $\{ \xi_j\}_{j \geq 1}$, a set of orthonormal basis functions that span $\mathcal H$, we have \begin{equation*}
\langle U_{2,t},\zeta \rangle =\langle \sum_{j=1} ^\infty \langle U_{2,t}, \xi_j \rangle \xi_j,\zeta \rangle = \sum_{j=1} ^\infty \langle U_{2,t}, \xi_j \rangle \langle \xi_j ,\zeta\rangle.
\end{equation*}
In this regard, the SIRF in \hyperref[{eqpartialirf1}]{\textup{\tagform@{\ref*{eqpartialirf1}}}} will be interpreted as the response of $y_t$ when all the basis functions move \textit{jointly} to the direction of $\zeta$. A similar observation was previously made by \cite{IR2021} under a parametric assumption, and it remains valid even in our setting without that assumption.
\begin{proposition}\label{prop: svar: identification: a} \normalfont Suppose that $\left[\begin{smallmatrix} y_t\\ X_t\end{smallmatrix} \right]$ satisfies \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}} and \hyperref[{eqrfvar}]{\textup{\tagform@{\ref*{eqrfvar}}}}. Then the following hold: below, $v_{1,t}$ (resp.\ $V_{2,t}$) is some $\mathbb{R}$-valued (resp.\ $\mathcal H$-valued) MA($h-1$) process, and
$v_{1,t}$ (resp.\ $V_{2,t}$) is not correlated with $\{y_{t-\ell},X_{t-\ell}\}$ for every $\ell \geq 0$.
\begin{enumerate}[(i)]
\item \label{eq: model: svar1}Under the conditions in Proposition \ref{prop: svar: identification}\ref{prop: svar: identification1}, $\operatorname{SIRF}_{12,h}$ is equivalent to the coefficient map $\langle \beta_h^y, \cdot \rangle:\mathcal H \to \mathbb{R}$ in \begin{equation}
y_{t+h}= \alpha_h ^y y_t + \langle \beta_{h} ^y, X_{t}\rangle + v_{1,t}, \label{prop2: eq1}
\end{equation}
where $\alpha_h ^y \in \mathbb{R}$ and $\beta_{h}^y \in \mathcal H$. Moreover, $\operatorname{SIRF}_{21,h}$ is equivalent to the coefficient $a_h ^X : \mathbb{R} \to \mathcal H$ in
\begin{equation}
X_{t+h} = \alpha_{h} ^X y_t + \beta_{h}^X y_{t-1} + \gamma_h^X X_{t-1} + V_{2,t},\label{prop2: eq2}
\end{equation}
where $\beta_h^X : \mathbb{R}\to\mathcal H$ and $\gamma_h^X : \mathcal H\to \mathcal H$.
\item Under the conditions in Proposition \ref{prop: svar: identification}\ref{prop: svar: identification2}, $\operatorname{SIRF}_{12,h}$ is equivalent to the coefficient map $\langle \alpha_{h} ^y,\cdot \rangle:\mathcal H \to \mathbb{R}$ in
\begin{equation}
y_{t+h} = \langle \alpha_{h} ^y, X_t \rangle + \langle {\beta}_{h} ^y, X_{t-1}\rangle + \gamma_h ^y y_{t-1} + v_{1,t},\label{prop2: eq4}
\end{equation}
where $\alpha_h^y \in \mathcal H$, $\beta_h^y \in \mathcal H$ and $\gamma_h^y \in \mathbb{R}$. Moreover, $\operatorname{SIRF}_{21,h}$ is equivalent to the coefficient $ \beta_h ^X: \mathbb{R} \to \mathcal H$ in
\begin{equation}
X_{t+h} = {\alpha}_h ^X X_t + {\beta}_{h} ^X y_{t} + V_{2,t}, \label{prop2: eq3}
\end{equation}
where $\alpha_h^X : \mathcal H\to \mathcal H$.
\end{enumerate}\end{proposition}
Proposition~\ref{prop: svar: identification: a}\ref{eq: model: svar1} provides a theoretical foundation for interpreting the map $\langle \beta_h, \cdot \rangle$ in the benchmark model \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} as $\operatorname{SIRF}_{12,h}$. Specifically, if the identification condition $\beta_{12} = 0$ holds and the true DGP follows \hyperref[{eq: model: svar}]{\textup{\tagform@{\ref*{eq: model: svar}}}}, then the coefficient map $\langle \beta_h, \cdot \rangle$
represents the SIRF at horizon $h$ when the functional structural error associated with $X_t$ experiences a shock $\zeta$. It may be possible to estimate the SIRFs directly from the RFVAR model by applying an appropriate regularization scheme required for an infinite dimensional setup. We provide a brief outline of this approach in Section~\ref{sec: svar: est} of the Supplementary Material. However, as detailed in the appendix, this approach makes statistical inference on the SIRFs much more challenging,
while our benchmark model simplifies the inference procedure by exclusively focusing on the outcome equation of interest.
Proposition~\ref{prop: svar: identification: a}\ref{eq: model: svar1} implies that $\operatorname{SIRF}_{21,h}$ can be characterized by a certain coefficient map in a function-on-function regression model with scalar control variables. We study how statistical inference can be implemented for this quantity in Section~\ref{sec_extension} of the Supplementary Material.
\subsection{Caveats of finite dimensional approximation\label{subsec: para}}
In empirical studies involving functional random variables, it is a common practice to first approximate the functional variable as a finite dimensional vector and then apply existing estimation or inference methods. However, as \cite{Nielsen2023} pointed out in the context of cointegration tests for functional time series, this finite dimensional approximation can lead to misspecification errors, causing misleading estimation and interpretation. To illustrate this, suppose
\begin{equation}
\widetilde \Upsilon_t = A \widetilde\Upsilon_{t-1} + \mathbf{u}_t, \label{eq: model: svar3}
\end{equation} where $ \widetilde \Upsilon_t = \left[\begin{smallmatrix} y_t\\ X_t\end{smallmatrix} \right] \in \mathbb{R}\times \mathcal H$. To study the above model, a popular approach is to transform $X_t$ in $\widetilde\Upsilon_t$ into a finite dimensional vector, say $\mathbf{x}_t$, in advance and then estimate a VAR model similar to \hyperref[{eq: model: svar3}]{\textup{\tagform@{\ref*{eq: model: svar3}}}} with $(y_t,\mathbf{x}_t')'$. The transformation from $X_t$ to $\mathbf{x}_t$ may be expressed as a projection of $X_t$ onto a finite dimensional space such that $P_0{\widetilde\Upsilon}_t = (y_t, \langle X_t, \xi_1 \rangle,\ldots, \langle X_t, \xi_K \rangle)'$ with some finite $K$. $\xi_j$ may be either a pre-specified parametric function (\citealp{IR2021}) or a principal component (\citealp{BHCJ2023}). However, as shown by \cite{Nielsen2023}, this projection operation does not preserve the VAR structure. Specifically, from \hyperref[{eq: model: svar3}]{\textup{\tagform@{\ref*{eq: model: svar3}}}}, we have
\begin{equation}
P_0 \widetilde \Upsilon_t = P_0 A P_0\widetilde\Upsilon_{t-1} + \mathbf{v}_t, \qquad \mathbf{v}_t = P_0 A (I-P_0) \widetilde\Upsilon_{t-1} + P_0 \mathbf{u}_t, \nonumber
\end{equation} If $P_0 \neq I$, $P_0\widetilde\Upsilon_{t}$ does not follow the VAR(1) structure in \hyperref[{eq: model: svar3}]{\textup{\tagform@{\ref*{eq: model: svar3}}}} since it is generally correlated with $\mathbf{v}_t$, which potentially invalidates many existing estimation and inference methodologies developed for VAR models; for instance, \cite{Nielsen2023} show that the cointegration rank test of \cite{Johansen1991,Johansen1995} is generally misleading in this setup. This finite dimensional approximation can be justified only when $P_0$ is close enough to $I$ (unless we consider the special case where $X_t$ can be fully expressed by a finite number of known basis functions, which is essentially equivalent to a finite dimensional setup). This implies that \( K \) must be sufficiently large. In functional linear models, it is well known that an increase in \( K \) can significantly raise the variance of coefficient estimators while reducing the regularization bias. Thus, choosing \( K \) in advance without accordance with the desired inferential methods may introduce an unbalanced bias-variance trade-off.
\section{Estimation}\label{sec:est}
\subsection{Representation of the proposed model}
To facilitate our discussion, it is useful to consider an alternative representation of our benchmark model in \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}}. We let $\widetilde{\mathcal H}$ denote the product Hilbert space, whose inner product is given by the sum of the inner products in $\mathbb{R}^m$ and $\mathcal H$ (see Section~\ref{sec_prelim2} of the Supplementary Material). With a slight abuse of notation, we let $\langle \cdot, \cdot \rangle$ denote the inner product defined on any of $\mathbb{R}^m$, $\mathcal H$, and $\widetilde{\mathcal H}$ (e.g., if $h_1, h_2 \in \mathbb{R}^m$ then $\langle h_1,h_2 \rangle=h_1'h_2$). Because the inner product is inherently defined for elements in the same space, there is little risk of confusion following this simplification. For any elements $h_j \in \mathcal H_j$ and $h_k\in \mathcal H_k$, where $\mathcal H_j$ and $\mathcal H_k$ can be any of $\mathbb{R}^m$, $\mathcal H$, and $\widetilde{\mathcal H}$, we let $\otimes$ denote the tensor product, defined by $h_j\otimes h_k (\cdot) = \langle h_j,\cdot \rangle h_k$ (see Section~\ref{sec_prelim} of the Supplementary Material). In particular, if $h_1 = \left[\begin{smallmatrix}h_{11}\\ h_{12}\end{smallmatrix}\right] \in \widetilde{\mathcal H}$ and $h_{2} = \left[\begin{smallmatrix}h_{21}\\ h_{22}\end{smallmatrix}\right] \in \widetilde{\mathcal H}$, then $h_1\otimes h_2$ can be understood as an operator matrix given by $\left[\begin{smallmatrix} h_{11} \otimes h_{21} & h_{12} \otimes h_{21} \\ h_{11} \otimes h_{22} & h_{12} \otimes h_{22}\end{smallmatrix}\right]$. We define the following $\widetilde{\mathcal H}$-valued random element $\Upsilon_t$ and its coefficient $\theta_h$:
\begin{equation} \label{eqdefY}
\Upsilon_t = \begin{bmatrix}\mathbf{w}_t\\ X_t\end{bmatrix} \quad\text{ and }\quad \theta_h = \begin{bmatrix}
\alpha_h \\\beta_h
\end{bmatrix}.
\end{equation}
Under this representation, $ \langle\theta_h, \Upsilon_t\rangle = \langle \alpha_h,\mathbf{w}_t \rangle + \langle \beta_h, X_t \rangle$, and \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} is simplified into the following regression model:
\begin{equation} \label{eq: model: benchmark: reduced2}
y_{t+h} = \langle\theta_h, \Upsilon_t\rangle + u_{h,t}.
\end{equation}
\subsection{Estimator}\label{sec:est1}
We hereafter assume that the variables $y_{t+h}$, $X_t$ and $\mathbf w_t$ have zero means. Extending to the case where these means are unknown is straightforward by considering their demeaned values, under the assumptions detailed shortly. Furthermore, we assume that $X_t$ and $\mathbf{w}_t$ are exogenous with respect to $u_{h,t}$. Then, by its definition in \hyperref[{eqdefY}]{\textup{\tagform@{\ref*{eqdefY}}}}, $\Upsilon_t$ satisfies the exogeneity condition such that $C_{\Upsilon u} = \mathbb{E}[u_{h,t}\Upsilon_t] = 0$. This in turn implies the following moment condition:
\begin{equation}
C_{\Upsilon y}\equiv \mathbb{E}[y_{t+h}\Upsilon_t] = \mathbb E[ \langle\Upsilon_t ,\theta_h \rangle \Upsilon_t ] = \mathbb E[\Upsilon_t \otimes \Upsilon_t] \theta_h \equiv C_{\Upsilon\Upsilon}\theta_h. \label{eqpopmoment}
\end{equation}
Under the identification condition to be detailed shortly, $\theta_h$ can be estimated from the sample counterpart of \hyperref[{eqpopmoment}]{\textup{\tagform@{\ref*{eqpopmoment}}}} given by
\begin{equation}
\widehat{C}_{\Upsilon y} = \widehat{C}_{\Upsilon\Upsilon}\bar{\theta}_h, \label{eqsammoment}
\end{equation}where $\widehat{C}_{\Upsilon y}$ and $\widehat{C}_{\Upsilon\Upsilon}$ are given by $
\widehat{C}_{\Upsilon y} = T^{-1}\sum_{t=1}^T y_{t+h} \Upsilon_t$ and $ \widehat{C}_{\Upsilon \Upsilon} = T^{-1} \sum_{t=1}^T \Upsilon_t \otimes \Upsilon_t$. However, in our setup, $\bar{\theta}_h$ is not generally obtainable from \hyperref[{eqsammoment}]{\textup{\tagform@{\ref*{eqsammoment}}}}, since $ \widehat{C}_{\Upsilon\Upsilon}$, an operator acting on $\widetilde{\mathcal H}$, is not invertible. We circumvent this issue by considering an estimator constructed using a regularized inverse of $\widehat{C}_{\Upsilon\Upsilon}$ based on its operator Schur complement. Specifically, we note that whenever convenient, ${C}_{\Upsilon\Upsilon}$ and $\widehat{C}_{\Upsilon\Upsilon}$ can be understood as operator matrices on $\widetilde{\mathcal H}$ such that
\begin{equation} \label{covblock}
{C}_{\Upsilon\Upsilon} = \begin{bmatrix}
\Gamma_{11} &\Gamma_{12} \\ \Gamma_{21}& \Gamma_{22}
\end{bmatrix} \quad\text{and}\quad
\widehat{C}_{\Upsilon\Upsilon} = \begin{bmatrix}
\widehat{\Gamma}_{11} & \widehat{\Gamma}_{12} \\ \widehat{\Gamma}_{21} & \widehat{\Gamma}_{22}
\end{bmatrix},
\end{equation}
where $\Gamma_{ij} = \mathbb{E}[g_{jt}\otimes g_{it}]$, $\widehat{\Gamma}_{ij} = T^{-1}\sum_{t=1}^T g_{jt}\otimes g_{it}$, $g_{1t} = \mathbf{w}_t$, and $g_{2t} = X_t$. That is, for any $\mathbf{a} \in \mathbb R^m$, $\Gamma_{11}\mathbf{a} = \mathbb E[ \langle \mathbf w_t, \mathbf a \rangle \mathbf{w}_t]$ and $\Gamma_{21}\mathbf{a} = \mathbb E[ \langle \mathbf w_t, \mathbf a \rangle X_t]$. Similarly, for any $ h\in \mathcal{H}$, $\Gamma_{12} h = \mathbb E[\langle X_t, h\rangle \mathbf w_t]$ and $\Gamma_{22} h = \mathbb E[\langle X_t, h \rangle X_t]$. Note that $\Gamma_{11}$ is the covariance of $ \mathbf{w}_t$ which is assumed to be invertible throughout the paper (see Assumption \ref{assum2}).
We then define the following operator Schur complements $\mathrm{S}$ (of ${C}_{\Upsilon\Upsilon}$) and $\widehat{\mathrm{S}}$ (of $\widehat{C}_{\Upsilon\Upsilon}$) (see \citealp[Section 2.2]{Bart2007}):
\begin{equation} \label{eqopschur}
\mathrm{S} ={\Gamma}_{22}-{\Gamma}_{21}{\Gamma}_{11}^{-1}{\Gamma}_{12}\quad\text{and}\quad \widehat{\mathrm{S}}=\widehat{\Gamma}_{22}-\widehat{\Gamma}_{21}\widehat{\Gamma}_{11}^{-1}\widehat{\Gamma}_{12}.
\end{equation}
{The operator Schur complement $\mathrm{S}$ can be understood as the covariance of the residuals obtained by regressing $X_t$ on $\mathbf{w}_t$, and $\widehat{\mathrm{S}}$ is its sample counterpart (see \citealp{fukumizu2004dimensionality}).} Lastly, to introduce our regularization scheme, we further represent $\mathrm{S}$ (resp.\ $\widehat{\mathrm{S}}$) with respect to its eigenvalues and eigenvectors $\{ \lambda_j,\nu_j \}_{j \geq 1}$ (resp.\ $\{ \widehat\lambda_j, \widehat \nu_j \}_{j \geq 1}$) as follows:
\begin{equation}\label{eqshur0}
{\mathrm{S}}= \sum_{j=1}^\infty {\lambda}_j {v}_j \otimes {v}_j\quad\text{and}\quad \widehat{\mathrm{S}}= \sum_{j=1}^\infty \widehat{\lambda}_j \widehat{v}_j \otimes \widehat{v}_j,
\end{equation}
where $\lambda_1\geq\lambda_2\geq\ldots \geq 0$ and $\widehat{\lambda}_1\geq\widehat{\lambda}_2\geq\ldots\geq 0$.
This representation is possible as $\mathrm{S}$ and $\widehat{\mathrm{S}}$ are self-adjoint, nonnegative, and compact (see \citealp{Bosq2000}, p.\ 34). The empirical eigenelements $\{\widehat{\lambda}_j , \widehat{v}_j\}$ can be obtained using the functional principal component analysis (FPCA). From \hyperref[{eqshur0}]{\textup{\tagform@{\ref*{eqshur0}}}}, noninvertibility of $\widehat{\mathrm{S}}$ is evident as its partial inverse $\widehat{\mathrm{S}}_K^{-1} = \sum_{j=1}^K \widehat{\lambda}_j ^{-1}\widehat{v}_j\otimes \widehat{v}_j$ grows without bound in the operator norm as $K$ increases. This consequently leads to noninvertibility of $\widehat{C}_{\Upsilon\Upsilon}$ (see \citealp{Bart2007}, p.\ 29). Therefore, to construct a regularized inverse of $\widehat{C}_{\Upsilon\Upsilon}$, we first construct regularized inverses of the operator Schur complements $\mathrm{S}$ and $\widehat{\mathrm{S}}$ as follows:
\begin{equation} \label{eqshur1}
\mathrm{S}_{\operatorname{\mathrm{K}_{\tau}}}^{-1} = \sum_{j=1}^{\operatorname{\mathrm{K}_{\tau}}} \lambda_j^{-1}v_j\otimes v_j \quad\text{and}\quad \widehat{\mathrm{S}}_{\operatorname{\mathrm{K}_{\tau}}}^{-1}= \sum_{j=1}^{\operatorname{\mathrm{K}_{\tau}}} \widehat{\lambda}_j^{-1} \widehat{v}_j \otimes \widehat{v}_j,\quad \text{ where}\quad \operatorname{\mathrm{K}_{\tau}} = \max\{j: \widehat{\lambda}_j^2 \geq \tau\} ,
\end{equation}
for the regularization parameter $\tau$, decaying to zero as $T$ increases. Then, we define the following regularized inverse $\widehat{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1}$ of $\widehat{C}_{\Upsilon\Upsilon}$:
\begin{equation} \label{eqreginv}
\widehat{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1} = \begin{bmatrix}
\widehat{\Gamma}_{11}^{-1} + \widehat{\Gamma}_{11}^{-1}\widehat{\Gamma}_{12} \widehat{\mathrm{S}}^{-1}_{\operatorname{\mathrm{K}_{\tau}}} \widehat{\Gamma}_{21}\widehat{\Gamma}_{11}^{-1} & -\widehat{\Gamma}_{11}^{-1} \widehat{\Gamma}_{12} \widehat{\mathrm{S}}^{-1}_{\operatorname{\mathrm{K}_{\tau}}} \\ -\widehat{\mathrm{S}}^{-1}_{\operatorname{\mathrm{K}_{\tau}}} \widehat{\Gamma}_{21}\widehat{\Gamma}_{11}^{-1} &
\widehat{\mathrm{S}}^{-1}_{\operatorname{\mathrm{K}_{\tau}}}
\end{bmatrix}.
\end{equation}
By using the regularized inverse \hyperref[{eqreginv}]{\textup{\tagform@{\ref*{eqreginv}}}} and the moment condition \hyperref[{eqsammoment}]{\textup{\tagform@{\ref*{eqsammoment}}}}, we define our estimator $\widehat\theta_h$ as follows:
\begin{equation}
\widehat{\theta}_h = \widehat{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1} \widehat{C}_{\Upsilon y}. \nonumber
\end{equation}
\begin{remark}
Similar to the literature in Section~\ref{subsec: para}, \cite{AP2006} and \cite{SHIN2009} resolve the ill-posed inverse problem similar to \hyperref[{eqsammoment}]{\textup{\tagform@{\ref*{eqsammoment}}}} by approximating $X_t$ with a finite number of eigenvectors of $\widehat\Gamma_{22}$.
On the other hand, we in this paper utilize the fact that the ill-posedness is directly associated by the Schur complement $\widehat{\mathrm{S}}$, so we regularize it to solve the inverse problem. The operator to be regularized $\widehat{\mathrm{S}}$ is the the covariance of the residuals obtained by regressing $X_t$ on $\mathbf{w}_t$, which is differentiated from $\widehat{\Gamma}_{22}$. Moreover, in contrast to the aforementioned literature whose focus lies on consistency, we provide a formal statistical inference procedure for $\theta_h$ based on the asymptotic properties of our proposed estimators, which is another novel contribution of this paper.
\end{remark}
\subsection{Identification and consistency} \label{sec:est1a}
Let $\ker A$ denote the kernel of $A$. To uniquely identify $\theta_h$ from \hyperref[{eqpopmoment}]{\textup{\tagform@{\ref*{eqpopmoment}}}}, we assume the following:
\begin{assumFLS} \label{assum1} $\langle x,\theta_h \rangle=0$ for all $x\in \ker C_{\Upsilon\Upsilon}$.
\end{assumFLS}
Note that $\Upsilon_t$ involves the $\mathcal H$-valued random variable $X_t$. In the literature on functional data analysis, it is commonly assumed that the covariance of such a functional random variable allows infinitely many nonzero eigenvalues. Therefore, as deduced from \citet[Proposition 2.1]{Mas2007} and the Riesz representation theorem (\citealp{Conway1994}, p.\ 13), the parameter of interest in \hyperref[{eqpopmoment}]{\textup{\tagform@{\ref*{eqpopmoment}}}}, $\theta_h$, is not uniquely identified as an element of $\widetilde{\mathcal H}$ if $\ker C_{\Upsilon \Upsilon} \neq \{0\}$. The failure of identification occurs because, for any $\psi \in \ker C_{\Upsilon\Upsilon}$, $C_{\Upsilon\Upsilon} ( \theta_h+\psi ) = C_{\Upsilon\Upsilon}\theta_h$. Thus, if $\theta_h$ satisfies \hyperref[{eqpopmoment}]{\textup{\tagform@{\ref*{eqpopmoment}}}}, $\theta_h+\psi$ also satisfies it. Assumption~\ref{assum1} prevents this failure of identification.
Our next assumption is related to asymptotic properties of $\widehat\theta_h$. Below, $\widehat{C}_{\Upsilon u} = T^{-1}\sum_{i=1} ^T u_{h,t} \Upsilon_t$ and $\Vert \cdot \Vert _{\operatorname{op}}$ is the operator norm (see Section~\ref{sec_prelim} of the Supplementary Material for its formal definition).
\begin{assumFLS} \label{assum2}
\begin{enumerate*}[(i)]
\item \label{assum2a} \hyperref[{eq: model: benchmark: reduced2}]{\textup{\tagform@{\ref*{eq: model: benchmark: reduced2}}}} holds with $\mathbb{E}[u_{h,t}\Upsilon_t]=0$;
\item\label{assum2b} $\{\Upsilon_t\}$, $\{u_{h,t}\}$ and $\{u_{h,t}\Upsilon_t\}$ are stationary and $L^4$-$m$-approximable;
\item \label{assum2c} $\|\widehat{C}_{\Upsilon u}\|_{\operatorname{op}} = O_p(T^{-1/2})$, $\|\widehat{C}_{\Upsilon\Upsilon} -C_{\Upsilon\Upsilon}\|_{\operatorname{op}} = O_p(T^{-1/2})$, and $\|\widehat{\Gamma}_{11}^{-1} - \Gamma_{11}^{-1}\|_{\operatorname{op}} = O_p(T^{-1/2})$.
\end{enumerate*}
\end{assumFLS}
The $L^4$-$m$-approximability in Assumption \ref{assum2}\ref{assum2b}, whose formal definition is provided in Section~\ref{Section_AFTS} of the Supplementary Material, is employed to use existing limit theorems. The assumption is not only widely adopted in the literature on stationary functional time series but also inclusive of many practical and interesting examples, such as the SVAR model in Section \ref{sec: svar}. Assumption~\ref{assum2}\ref{assum2c} contains high-level conditions on the limiting behavior of $\widehat{C}_{\Upsilon \Upsilon}$ and $\widehat{C}_{\Upsilon u}$, and is not restrictive given the stationarity of $\{\Upsilon_t\}$ and $\{ u_{h,t}\Upsilon_t\}$. Some primitive sufficient conditions can be found in \citet[Chapter 2]{Bosq2000}.
We impose the following conditions to characterize the rate of convergence of $\widehat\theta_h$.
\begin{assumFLS} \label{assum3} For a generic constant $\mathtt{C}>0$, the following holds: \begin{enumerate*}[(i)]
\item\label{assum3a}
$\lambda_j^2 \leq \mathtt{C} j^{-\rho}$ and $\lambda_j^2-\lambda_{j+1}^2 \geq \mathtt{C} j^{-\rho-1}$ for $\rho>2$;
\item\label{assum3b} $|\langle \beta_h, v_j \rangle| \leq \mathtt{C} j^{-\varsigma}$ for $\varsigma > 1/2$.
\end{enumerate*}
\end{assumFLS}
When $\mathrm{S}$ is given by the standard covariance of the functional explanatory variable, Assumption~\ref{assum3}\ref{assum3a} reduces to a standard assumption commonly used in the literature on functional linear models.
Given that $\sum_{j=1}^\infty \lambda_j < \infty$ must hold
(see Lemma \ref{lem1}), the first condition $\lambda_j^2 \leq \mathtt{C} j^{-\rho}$ for some $\rho>2$ is natural (obviously, $\rho\leq 2$ may not result in the summability of $\{\lambda_j\}_{j\geq 1}$). By the latter condition, we require that the eigenvalues of $\mathrm{S}$ are well separated. As documented in a similar context concerning functional linear models (see e.g., \citealp{Hall2007,imaizumi2018,seong2021functional}), this separation is crucial for achieving sufficient accuracy in the estimation of eigenelements. Assumption~\ref{assum3}\ref{assum3b} can be understood as a smoothness condition on $\beta_h$ with respect to the eigenvectors $\{v_j\}_{j \geq 1} $. In contrast to $\beta_h$, a similar restriction is not necessary for $\alpha_h$ associated with the vector-valued random variable $\mathbf{w}_t$. Considering that $\beta_{h}$ is an element in a Hilbert space satisfying $\sum_{j=1}^\infty \langle \beta_h, v_j \rangle^2 <\infty$, Assumption \ref{assum3}\ref{assum3b} does not seem too restrictive.
The following theorem states consistency and the rate of convergence of $\widehat\theta_h$.
\begin{theorem} \label{thm1} Suppose Assumptions \ref{assum1}--\ref{assum3} hold and $T\tau^{1+4/\rho} \to \infty$. Then, \begin{equation} \|\widehat{\theta}_h - \theta_h\| = O_p(T^{-1/2}\tau^{-1/2-2/\rho} + \tau^{(2\varsigma-1)/2\rho}).\nonumber\end{equation}
\end{theorem}
In Theorem \ref{thm1}, the choice of the regularization parameter $\tau$ to ensure the consistency of the proposed estimator depends on $\rho$, allowing $\tau$ to decay at a faster rate as $\rho$ increases. Assuming that the same regularization parameter is employed, an increase in $\rho$ implies a faster convergence of $\widehat{\theta}_h$ to $\theta_h$.
It is worth noting that the consistency is established if $\tau$ satisfies $T\tau^{3} \to \infty$ (such as e.g., $\tau = T^{-1/3 + \epsilon}$ or $T^{-1/3} \log^{\epsilon} T$ for small $\epsilon>0$) as long as $\rho >2$ as assumed in Assumption \ref{assum3}; of course, $\tau$ decaying at a slower rate can also be considered in practice, without affecting consistency. Nevertheless, this naive selection of $\tau$ may be particularly advantageous for practitioners with little knowledge of these eigenvalues.
\subsection{Statistical inference based on local asymptotic normality}\label{sec:est1b}
Given that $\langle \beta_h, \zeta\rangle$ can naturally be interpreted as the response to a perturbation applied to $X_t$, as in \hyperref[{def: local: irf}]{\textup{\tagform@{\ref*{def: local: irf}}}}, and in some special cases, such as in Proposition \ref{prop: svar: identification: a}, it can further be interpreted as $\operatorname{SIRF}_{12,h}(\zeta)$, it is of our interest to conduct statistical inference on this quantity. More generally, we consider inference on \( \langle \theta_h, \zeta \rangle \) for any \( \zeta \in \widetilde{\mathcal H} \). To clarify the perturbation \( \zeta \) in the subsequent discussion and avoid potential confusion regarding the space in which \( \zeta \) takes values, we let
\begin{equation} \label{eqzetadecom}
\zeta =\left[\begin{matrix}\zeta_{1}\\ \zeta_{2}\end{matrix}\right] \in \widetilde{\mathcal H}, \quad\text{where}\quad \zeta_1 \in \mathbb{R}^m \quad\text{and}\quad \zeta_2 \in \mathcal H .
\end{equation}
By setting \( \zeta_1 = 0 \), inference on \( \langle \theta_h, \zeta \rangle \) reduces to inference on \( \langle \beta_h, \zeta_2 \rangle \) for \( \zeta_2 \in \mathcal H \). Meanwhile, if $\zeta_2 = 0$, it reduces to inference on $\langle \alpha_h, \zeta_1\rangle \equiv \alpha_h ' \zeta_1$.
Let $\widehat{u}_{h,t}= y_{t+h}- \langle \Upsilon_t, \widehat{\theta}_h\rangle$ and $\mathrm{k}(\cdot)$ be a standard weight function to be specified shortly. Then, under Assumption \ref{assum2},
we define the following long-run covariance operator and its sample counterpart (see \citealt[Theorem 2]{berkes2013weak}): below, ${\mathcal U}_{h,t} = {u}_{h,t}\Upsilon_{t}$ and $\widehat{\mathcal U}_{h,t} = \widehat{u}_{h,t}\Upsilon_{t}$.
\begin{equation*}
\Lambda_{\mathcal U} =\sum_{s=-\infty}^\infty \mathbb{E}[\mathcal U_{h,t} \otimes \mathcal U_{h,t-s}] \quad \text{and}\quad \widehat{\Lambda}_{\mathcal U} = \frac{1}{T} \sum_{s=-\mathsf{b}}^{\mathsf{b}}\mathrm{k}\left(\frac{s}{\mathsf{b}}\right)\left( \sum_{1 \leq t, t-s \leq T} \widehat{\mathcal U}_{h,t} \otimes \widehat{\mathcal U}_{h,t-s} \right).
\end{equation*}
We also define the quantity $c_{m,j}$ as
\begin{equation} \label{eqcm}
c_{m,j}(\zeta) = \bar{\lambda}_j^{-1}\langle \zeta,\bar{v}_j \rangle\bigg/\sqrt{\sum_{j=1}^{m} \bar{\lambda}_j^{-2}\langle \zeta,\bar{v}_j \rangle^2} ,
\end{equation}
where $\{ \bar{\lambda}_j \}$ and $ \{ \bar{v}_j\}$ denote the eigenvalues and eigenvectors of $C_{\Upsilon\Upsilon}$, which are generally different from the eigenelements of $\mathrm{S}$. The quantity satisfies $\sum_{j=1}^{m} c_{m,j}^2(\zeta) = 1$ for every $m$ unless $\langle \zeta,\bar{v}_j \rangle \neq 0$ for at least one $j$.
To establish local asymptotic normality of $\widehat\theta_h$, we employ the following assumption:
\begin{assumFLS} \label{assum4}
\begin{enumerate*}[(i)]
\item \label{assum4a} $\mathsf{b}=\mathsf{b}(T) \to \infty$, $\mathsf{b}(T)/T \to 0$ and $\mathrm{k}(\cdot)$ is an even function with $\mathrm{k}(0)=1$, $\mathrm{k}(s) = 0$ if $|s| > c$ for some $c>0$ and $\mathrm{k}(\cdot)$ is continuous on $[-c,c]$;
\item \label{assum4aa} $\zeta \notin \ker C_{\Upsilon\Upsilon}$ and $\Lambda_{\mathcal U} \bar{v}_j \neq 0$ for $\bar{v}_j$ corresponding to $\bar{\lambda}_j >0$;
\item \label{assum4bb} $\sum_{j=1}^{m} \sum_{\ell=1}^{m} c_{m,j}(\zeta) c_{m,\ell}(\zeta) \langle \Lambda_{\mathcal U}\bar{v}_j,\bar{v}_{\ell} \rangle \to \mathtt{C}>0$ as $m\to \infty$; \item \label{assum4cc} $\sup_{1\leq t\leq T} \|\Upsilon_t\| = O_p(1)$.
\end{enumerate*}
\end{assumFLS}
Assumption \ref{assum4}\ref{assum4a} is adopted from \cite{horvath2013estimation} and imposed for convenience in our proof of consistency of $\widehat{\Lambda}_{\mathcal U}$. Assumption \ref{assum4}\ref{assum4aa} is imposed to ensure nondegenerate convergence rate in our asymptotic normality result; in the special case where $\Lambda_{\mathcal U} = \sigma_u^2 C_{\Upsilon\Upsilon}$ for some $\sigma_u^2>0$,\footnote{This can happen when $\{u_t\}$ is a homogeneous martingale difference with respect to the filtration $\mathfrak F_t=\sigma(\{u_{s}\}_{s\leq t-1},\{\Upsilon_{s}\}_{s\leq t})$ as in the case considered by \cite{seong2021functional}} the former condition $\zeta \notin \ker C_{\Upsilon\Upsilon}$ implies the latter, and hence Assumption \ref{assum4}\ref{assum4aa} can be simplified. Assumptions~\ref{assum4}\ref{assum4bb} and \ref{assum4}\ref{assum4cc} are technical requirements facilitating our asymptotic analysis, and are not quite restrictive. In particular, as $\Lambda_{\mathcal U}$ is a covariance operator (\citealt[Theorem 1.7]{Bosq2000}), the quantity in Assumption~\ref{assum4}\ref{assum4bb} converges to a nonnegative constant for every $\zeta \in \widetilde{\mathcal H}$. Thus, what we require of $\zeta$ is solely to ensure the positivity of the limit. Assumption \ref{assum4}\ref{assum4cc} is similar to standard assumptions found in the literature (see \citealt[Theorems 2.12-2.14]{Bosq2000}).
We employ the following assumption: below, we let $X_{w,t} = X_t-\Gamma_{21}\Gamma_{11}^{-1}\mathbf{w}_t$, $r_t(j,\ell) = \langle X_t,v_j \rangle \langle X_{w,t},v_{\ell} \rangle - \mathbb E[\langle X_t,v_j \rangle \langle X_{w,t},v_{\ell} \rangle]$ for $j,\ell \geq 1$, and $\mathtt{C}$ be a generic positive constant.
\begin{assumFLS} \label{assum5}
\begin{enumerate*}[(i)]
\item\label{assum5a} $\mathbb{E}[\langle X_t,v_j \rangle^4] \leq \mathtt{C} \lambda_j^2$, $\mathbb{E}[\langle X_{w,t},v_j \rangle^4] \leq \mathtt{C} \lambda_j^2$, and for some $\widetilde{\mathtt{C}} > 1$ and $s\geq 1$, $\mathbb{E}[r_t(j,\ell)r_{t-s}(j,\ell)]\leq \mathtt{C} s^{-\widetilde{\mathtt{C}}}\mathbb{E}[r_t^2(j,\ell)]$;
\item\label{assum5b} for some $\delta>1/2$, $|\langle \zeta_2,v_j \rangle| \leq \mathtt{C} j^{-\delta}$ and $|\langle \zeta_1, {\Gamma}_{11}^{-1}{\Gamma}_{12} v_j \rangle| \leq \mathtt{C} j^{-\delta}$.
\end{enumerate*}
\end{assumFLS}
A similar but slightly different assumption can be found in \cite{seong2021functional} and references therein. In particular, Assumption~\ref{assum5}\ref{assum5b} pertains to the smoothness of $\zeta$; similar to Assumption~\ref{assum3}\ref{assum3b}, it is natural to consider $\delta >1/2$ as $\zeta_2\in\mathcal H$ and $ {\Gamma}_{11}^{-1}{\Gamma}_{12}$ is a bounded linear functional.\footnote{There exists unique element $h \in \mathcal H$ such that ${\Gamma}_{11}^{-1}{\Gamma}_{12}v_j = \langle h, v_j \rangle$ by the Riesz representation theorem (see \citealp{Conway1994}, p.\ 13), and $\sum_{j=1}^\infty\|{\Gamma}_{11}^{-1}{\Gamma}_{12}v_j\|^2 < \infty$ as $\{v_j\}$ is an orthonormal basis in $\mathcal H$.}
Lastly, we let $\widehat{P}_{\operatorname{\mathrm{K}_{\tau}}}=\widehat{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1} \widehat{C}_{\Upsilon\Upsilon} $ and $P_{\operatorname{\mathrm{K}_{\tau}}} = {C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1} {C}_{\Upsilon\Upsilon}$, where ${C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1}$ is defined by replacing $\widehat\Gamma_{ij}$ and $\widehat\mathrm{S}_{\operatorname{\mathrm{K}_{\tau}}}^{-1}$ in \hyperref[{eqreginv}]{\textup{\tagform@{\ref*{eqreginv}}}} with $\Gamma_{ij}$ and $\mathrm{S}_{\operatorname{\mathrm{K}_{\tau}}}^{-1}$ respectively. Then, $\widehat{P}_{\operatorname{\mathrm{K}_{\tau}}} $ and $P_{\operatorname{\mathrm{K}_{\tau}}} $ allow the following representations:
\begin{equation}\label{eqprojections}
\widehat{P}_{\operatorname{\mathrm{K}_{\tau}}}
= \begin{bmatrix}
I_1 & \widehat{\Gamma}_{11}^{-1}\widehat{\Gamma}_{12}(I_2-\widehat{\Pi}_{\operatorname{\mathrm{K}_{\tau}}}) \\ 0 & \widehat{\Pi}_{\operatorname{\mathrm{K}_{\tau}}}
\end{bmatrix}\quad\text{ and }\quad {P}_{\operatorname{\mathrm{K}_{\tau}}}
= \begin{bmatrix}
I_1 & {\Gamma}_{11}^{-1}{\Gamma}_{12}(I_2-{\Pi}_{\operatorname{\mathrm{K}_{\tau}}}) \\ 0 & {\Pi}_{\operatorname{\mathrm{K}_{\tau}}}
\end{bmatrix},
\end{equation}
where $I_1$ (resp.\ $I_2$) is the identity map on $\mathbb{R}^{m}$ (resp.\ $\mathcal H$) and
\begin{equation}
\widehat{\Pi}_{\operatorname{\mathrm{K}_{\tau}}}=\sum_{j=1}^{\operatorname{\mathrm{K}_{\tau}}} \widehat{v}_j \otimes \widehat{v}_j\quad\text{and}\quad \Pi_{\operatorname{\mathrm{K}_{\tau}}}=\sum_{j=1}^{\operatorname{\mathrm{K}_{\tau}}} v_j \otimes v_j. \nonumber
\end{equation}
That is, $\Pi_{\operatorname{\mathrm{K}_{\tau}}}$ is the projection onto the span of $\{v_j\}_{j=1}^{\operatorname{\mathrm{K}_{\tau}}}$ and $\widehat{\Pi}_{\operatorname{\mathrm{K}_{\tau}}}$ is its sample counterpart.
Using $\widehat{P}_{\operatorname{\mathrm{K}_{\tau}}}$ and ${P}_{\operatorname{\mathrm{K}_{\tau}}}$, we decompose $\langle \widehat{\theta}_h-\theta_h,\zeta\rangle$ into $ \widehat{\Theta}_1+\widehat{\Theta}_{2A}+\widehat{\Theta}_{2B}$ such that
\begin{equation} \label{eqdecomtheta}
\widehat{\Theta}_1 = \langle \widehat{\theta}_h-\widehat{P}_{\operatorname{\mathrm{K}_{\tau}}}\theta_h, \zeta \rangle,\quad \widehat{\Theta}_{2A}= \langle\widehat{P}_{\operatorname{\mathrm{K}_{\tau}}}\theta_h - P_{\operatorname{\mathrm{K}_{\tau}}}{\theta}_h,\zeta\rangle, \text{ and } \widehat{\Theta}_{2B}= \langle {P}_{\operatorname{\mathrm{K}_{\tau}}}\theta_h - {\theta}_h,\zeta\rangle.
\end{equation}
\begin{theorem} \label{thm2} Suppose that Assumptions \ref{assum1}-\ref{assum4} hold and $T\tau^{2+4/\rho} \to \infty$. Then the following holds:
\begin{enumerate}[(i)]
\item\label{thm2i0} $\sqrt{{T}/{\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}}\widehat{\Theta}_1 \to_d N(0,1)$ for $ {\psi}_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) = \langle {\Lambda}_{\mathcal U}{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1}\zeta,{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1}\zeta \rangle$.
\end{enumerate}
If Assumption \ref{assum5} is additionally satisfied with $\rho/2 + 2 < \varsigma+ \delta$,
the following hold:
\begin{enumerate}[(i)]\addtocounter{enumi}{1}
\item\label{thm2i1} If $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) \to_p \infty$, $
\sqrt{{T}/{\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}}\widehat{\Theta}_{2A} \to_p 0$.
\item \label{thm2i2} If $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) \to_p \infty$ and $T^{1/2}\tau^{(\delta+\varsigma-1)/\rho} \to 0$,
$\sqrt{{T}/{\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}}\widehat{\Theta}_{2B} \to_p 0.$
\end{enumerate}
The results in \ref{thm2i0}-\ref{thm2i2} hold when $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)$ is replaced by $\widehat{\psi}_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) = \langle \widehat{\Lambda}_{\mathcal U} \widehat{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1}\zeta, \widehat{C}_{\Upsilon\Upsilon,\operatorname{\mathrm{K}_{\tau}}}^{-1}\zeta \rangle$.
\end{theorem}
The quantity $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)$ in Theorem~\ref{thm2} plays an important role as a normalizing factor in our asymptotic normality result.
With an obvious adaptation of the discussion given by \citet[Remarks 4 and 11]{seong2021functional}, we know that it is convergent only on a strict subspace of $\widetilde{\mathcal H}$, and ${\psi}_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)$ is likely to diverge (at a rate slower than $T$; see Remark \ref{remadd2}) for $\zeta$ arbitrarily chosen by practitioners.
Similar to Theorem \ref{thm1}, the decay rate of the regularization parameter in Theorem~\ref{thm2} depends on $\rho$, and a larger value of $\rho$ leads to a faster convergence. Note that $\tau$ satisfying $T\tau^{4} \to \infty$ (e.g., $\tau = T^{-1/4 + \epsilon}$ or $T^{-1/4} \log^{\epsilon} T$ for small $\epsilon>0$) meets the requirement for Theorem~\ref{thm2}\ref{thm2i0}, provided that $\rho >2$ as assumed in Assumption~\ref{assum3}. Such a choice of $\tau$ may be preferred for practitioners. Of course, as in Theorem~\ref{thm1}, we may also consider $\tau$ decaying to zero at a slower rate, but this is not recommended due to the condition $T^{1/2}\tau^{(\delta+\varsigma-1)/\rho} \to 0$ for the result given in Theorem~\ref{thm2}\ref{thm2i2}. This condition can be understood as requiring sufficiently large $\varsigma$ and $\delta$ for a given $\tau$; moreover, as the decay rate of $\tau$ decreases, stricter requirements are imposed on $\varsigma$ and $\delta$.
\begin{remark} \label{remadd2}
The normalizing factor $\sqrt{{T}/{\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}}$ can diverge at a different rate depending on the choice of $\zeta$, and the factor is not stochastically bounded for any choice of $\zeta$. From \hyperref[{addeqrem}]{\textup{\tagform@{\ref*{addeqrem}}}} and Assumption \ref{assum3}, it is straightforward to see that ${\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}$ is bounded above by $O_p(\operatorname{\mathrm{K}_{\tau}}^{\rho})$, which is in turn $O_p(\tau^{-1})$ as shown in our proof of Theorem \ref{thm1}. Since $T\tau \to \infty$, this implies that ${\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}/T$ must decay to zero (i.e., $\sqrt{{T}/{\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta)}}$ is divergent) regardless of the choice of $\zeta$.
\end{remark}
\begin{remark}\normalfont
If $\tau = T^{-1/4} \log^{\epsilon} T$ for $\epsilon>0$, then $T^{1/2}\tau^{(\delta+\varsigma-1)/\rho}
=T^{(2\rho+1-\delta-\varsigma)/(4\rho)}\times$ $(\log T)^{\epsilon(\delta+\varsigma-1)/\rho}$. Since $\varsigma>1/2$ and $\delta>1/2$ under the employed assumptions, either (i) $\delta+\varsigma > 2\rho+1$ or (ii) $\varsigma \geq 2\rho + 1/2$ ensures the condition that $T^{1/2}\tau^{(\delta+\varsigma-1)/\rho} =o(1)$. In case (ii), Theorems \ref{thm2}\ref{thm2i1}-\ref{thm2i2} hold for any $\delta > 1/2$. More generally, if $\tau = T^{-\rho/(2\rho+4)} \log^{\epsilon}T$, it can be shown that the aforementioned condition is satisfied under either of (i)$'$ $\delta+\varsigma>\rho+3$ or (ii)$'$ $\varsigma \geq \rho + 5/2$.
Given that $\rho>2$, (i)$'$ (resp.\ (ii)$'$) is weaker than (i) (resp.\ (ii)), but they are nearly identical if $\rho$ is close to $2$.
\end{remark}
If $\zeta$ is restricted to have a nonzero element only in $\mathbb{R}^m$ (i.e., $\zeta_2 = 0 \in \mathcal{H}$) and $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) \to_p C_{\zeta} <\infty$, it can be shown that $\langle \widehat{\theta}_h - \theta_h, \zeta \rangle$ converges in distribution to a normal random variable at the rate of $\sqrt{T}$. However, we here emphasize that the former condition, $\zeta_2 = 0$, does not guarantee the latter condition, $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) \to_p C_{\zeta}$. Thus, it is common to observe a nonparametric rate of convergence even for the parameters associated with the finite dimensional random variables in this setup (see Remark~\ref{remmas} for more details).
\begin{corollary} \label{coradd}
Suppose that Assumptions \ref{assum1}-\ref{assum5} hold, $\psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) \to_p C_{\zeta} < \infty$, and $\zeta = \left[\begin{smallmatrix}\zeta_1\\0\end{smallmatrix}\right]$ where $\zeta_1 \in \mathbb{R}^m$ and $0\in \mathcal H$. Then,
\begin{equation}
\sqrt{T}\langle \widehat{\theta}_h-\theta_h,\zeta\rangle = \sqrt{T} \langle \widehat \alpha_h - \alpha_h, \zeta_1 \rangle \to_d N(0, C_{\zeta}).
\end{equation}
\end{corollary}
\begin{remark} \label{remmas}
As deduced from our results in Theorem \ref{thm2} and the discussions provided by \cite{Mas2007}, the convergence rate of the proposed estimator depends on the choice of $\zeta$ and is generally slower than $\sqrt{T}$. Moreover, even when we concern the coefficients associated with finite dimensional elements, the parametric $\sqrt{T}$-convergence rate of $\widehat\alpha_h$ requires $\lim_{\operatorname{\mathrm{K}_{\tau}} \to \infty} \psi_{\operatorname{\mathrm{K}_{\tau}}}(\zeta) < \infty$, which, in turn, necessitates sufficient smoothness of $\zeta$ with respect to the eigenvectors ${\bar{v}_j}$ of ${C}_{\Upsilon\Upsilon}$.
A similar condition can be found in the literature on partially linear functional regression models, including \cite{AP2006} and \citet[Theorem 3.1 and eqn.\ (20)]{SHIN2009}. Given its nontrivial nature, this condition reflects the cost of implementing inference in a model involving functional variables.
\end{remark}
\section{Empirical Application \label{sec:emp}}
\subsection{US Economic Sentiment Quantile Curves\label{sec:emp1}}
In this section, we study how economic variables respond to random perturbations introduced to economic sentiment distribution in the US. To measure sentiment, we adopt \citepos{Barbaglia2023} data, which quantify sentiment on specific economic subjects from sentences published in major US news articles. The sentiment measure is defined by each token (i.e., word) and its neighboring sentences, and includes both intensity and tone for each subject. We refer interested readers to \cite{Barbaglia2023} for a detailed definition of the sentiment measure.
Among others, we particularly focus on the daily sentiment measure on the \textsl{economy}, and use it to construct monthly quantile sentiment curves. The quantile curves are represented by 31 Fourier basis functions, which become our functional predictor, and are reported in Figure~\ref{fig: 1}. As mentioned earlier in Section~\ref{sec: model} and also by \cite{Barbaglia2023}, the sentiment quantile and its associated distribution seem to be closely related to business cycle fluctuations: during recessive periods, the overall quantiles tend to shift downward and exhibit a steeper slope coefficient, indicating a larger dispersion in the sentiment distribution. This could be interpreted as a higher level of disagreement about economic sentiment during recessive periods.
The dependent variable $y_t$ is specified to monthly Total Nonfarm Payroll (PAYEMS) in the US. The data are provided by the Federal Reserve Bank. We follow \cite{Barbaglia2023} and let $\mathbf{w}_t$ in \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} be the vector consisting of the Chicago Fed National Activity Index (CFNAI), the National Financial Conditions Index (NFCI), and \citepos*{ADS2009} measure of economic activity (ADS). Because NFCI and ADS are observed at a higher frequency, we average them over months to match the frequency of the target variable. The sample runs from January 1984 to December 2021, with a total size of 456 observations.
The regularization parameter of our estimator is chosen by roughly considering $\rho$ satisfying Assumption~\ref{assum3}. Specifically, if $\mathtt{C}=1$, Assumption~\ref{assum3}\ref{assum3b} tells us that $\rho \geq \rho^\ast \equiv -(\log (\lambda_j ^2 - \lambda_{j+1} ^2) /\log j)-1$. Thus, we set $\widetilde \rho$ to $\lceil 100 \rho ^\ast \rceil/100$, where $\lceil \cdot \rceil$ is the ceiling function. Then the regularization parameter is set to $0.01\Vert \widehat C_{\Upsilon\Upsilon}\Vert_{\text{HS}} T^{-\widetilde \rho/(\widetilde \rho+2)}$, where $\Vert \cdot \Vert_{\text{HS}}$ denotes the Hilbert-Schmidt norm. This approach may provide practical insights into the magnitude of $\rho$, thereby allowing us to choose $\tau$ without violating the assumption. We investigate the finite-sample performance of our estimator computed with this regularization parameter in Section~\ref{sec:sim}.
\begin{figure}[h!]
\caption{ Responses of Growth in Total Nonfarm Payroll to Sentiment Shocks} \label{fig: emp-1}
\includegraphics[width = \textwidth, height= 0.2\textheight]{figure/IRF_PAYEMS_A4.jpeg}
\flushleft{\scriptsize{Notes: The figure reports the impulse response estimates when the shocks are specified to $\zeta_1, \zeta_2$, $\zeta_3$, and $\zeta_4$ in Figure~\ref{fig: 2} respectively. The shocks are normalized to a unit norm before analysis. The solid lines are reporting the pointwise estimates and grey areas report their 90\% confidence interval.}}
\end{figure}
We first estimate the impulse responses of PAYEMS to shocks given to the sentiment distribution using our benchmark model in \hyperref[{eq: model: benchmark}]{\textup{\tagform@{\ref*{eq: model: benchmark}}}} for $ h \in\{1,2,\ldots,12\}$, without a structural interpretation. We let $\zeta$ be the shocks considered in Figure \ref{fig: 2}, with their norms normalized to one. The normalization is taken to ensure reasonable comparison between their effects. The estimation results are reported in Figure~\ref{fig: emp-1}, where solid lines indicate $\{ \langle \beta_h, \zeta \rangle \}_{h=1}^{12}$ in \hyperref[{eq: sirf}]{\textup{\tagform@{\ref*{eq: sirf}}}} and the shaded area represents their 90\% confidence intervals computed using the local asymptotic normality result in Theorem~\ref{thm2}.
As mentioned earlier, the first two shocks are associated with positive and negative shifts in the sentiment distribution. Meanwhile, in the last two columns, we consider cases where sentiment about the economy exhibits greater dispersion, with and without a negative shift, respectively. The responses in Figure~\ref{fig: emp-1} have different vertical scales, which may require caution in their interpretation.
Figure~\ref{fig: emp-1} suggests that the effects of the distributional shocks tend to reach their peak after 6 to 8 months, and their magnitude tends to decrease as the forecasting horizon increases.
In particular, when the average sentiment about the economy becomes more optimistic (resp.\ pessimistic), payroll growth appears to significantly increase (resp.\ decrease), especially during the first six months. The change in disagreement levels in sentiment appears to have a smaller effect on PAYEMS growth compared to locational shifts in the distribution.
Lastly, we study if our estimation approach produces any improvement compared to the existing estimators, by using the empirical median absolute prediction error (MAPE). The MAPEs are computed based on a rolling window with two different test sets, selected to include approximately 60\% and 70\% of the observations from the total sample. In addition to ours, denoted by SCInv, we consider two alternative estimators for comparison: PCA-FR and PCA-SVAR. The PCA-FR is obtained by applying \citepos{SHIN2009} estimation approach to our benchmark model. The PCA-SVAR is computed using the standard estimation methods developed for the recursive SVAR model in a finite dimensional setting, after pre-applying the FPCA-based dimension reduction to $X_t$.
\begin{table}[tbp]
\caption{Estimated Median Absolute Prediction Errors \label{tab: emp-1}}
\vskip -8pt
\small \centering
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcccccccccc}
\toprule
& \multicolumn{5}{c}{2006 Q4$\sim$} & \multicolumn{5}{c}{2010 Q3$\sim$} \\ \cmidrule{2-6}\cmidrule{7-11}
h & 1 & 2 & 3 & 4 & 5 & 1 & 2 & 3 & 4 & 5 \\ \midrule
SCInv & 0.055 & 0.060 & 0.067 & 0.067 & 0.063 & 0.054 & 0.058 & 0.054 & 0.059 & 0.053 \\
PCA-FR & 0.055 & 0.060 & 0.067 & 0.066 & 0.062 & 0.054 & 0.057 & 0.054 & 0.059 & 0.053 \\
PCA-SVAR & 0.060 & 0.065 & 0.061 & 0.069 & 0.072 & 0.056 & 0.064 & 0.059 & 0.058 & 0.056 \\
\bottomrule
\end{tabular*}
\flushleft{\scriptsize{Notes: The table reports the empirical median absolute prediction errors computed using our estimator (SCInv) and two FPCA-based estimators: PCA-FR and PCA-SVAR. PCA-FR is obtained by applying \citepos{SHIN2009} estimator to our benchmark model. PCA-SVAR is derived by applying the SVAR model to the vector consisting of $\mathbf w_t$ and $X_t$, with its dimension reduced using the FPCA approach. The regularization parameters are set to choose the first two score functions for all the estimators. }}
\end{table}
Table~\ref{tab: emp-1} summarizes MAPEs. To mitigate a potential effect of different choices of regularization parameters on estimation results, we fix $\tau$ so that all the estimators select the first two score functions, regardless of the forecasting horizon. This choice is based on empirical data, which suggests that approximately 99\% of the variations are explained by the first two scores. Overall, it seems that our estimator outperforms the PCA-SVAR in terms of smaller MAPEs, and the superior performance is more significantly observed as the forecasting horizon increases. An interesting observation is that the PCA-FR approach also outperforms the PCA-SVAR, although both estimators are computed with the same empirical eigenelements, and thus the PCA-FR can be understood as a single-equation estimation of SVAR models. This may be attributed to the fact that, in the PCA-SVAR estimation, the sentiment quantile curve serves as both the target (i.e., dependent) and prediction (i.e., explanatory) variables, whereas it is only used as a predictor in the PCA-FR approach. Therefore, the PCA-SVAR estimator is likely to be subject to a larger bias associated with \textsl{double truncation}. Meanwhile, in this particular example, the two functional approaches, the SCInv and the PCA-FR, produce similar MAPEs regardless of forecasting horizons.
\subsection{Impact of Functional Monetary Policy Shocks \label{sec:emp2}}
In this section, we employ \citepos{IR2021} functional monetary policy shocks and study the monetary policy impact on the inflation growth by using the proposed estimator. The data span the period from January 1995 to June 2016 and the functional shock is reported in Figure~\ref{fig: ir: 1}. Under Assumption~I in \cite{IR2021}, the slope coefficient of the shock ($\beta_h$ in our notation) can be understood as the impulse response to monetary policy shocks, with the effect represented by a function. In this section, we overall follow \cite{IR2021}. The control vector $\mathbf w_t$ is given by the set of the first two lagged dependent variables. $y_t$ is given by the inflation growth rate. The impulse response $\{ \langle \beta_h,\zeta\rangle \}_{h=1}^{20}$ is estimated with $\zeta$ corresponding to shocks observed on three specific dates of interest identified by the authors: 9/1998, 2/1999, and 1/2007. Note that the shocks and the changes in yield curves induced by them in the current study may differ slightly from those in \cite{IR2021} because we do not represent them into level and curvature components.
\begin{figure}[h!]
\centering
\caption{Inflation Response to Functional Shocks in \cite{IR2021}}\label{fig: emp: ir}
\begin{subfigure}{.32\textwidth}\caption{Change in Sep 1998}
\includegraphics[width = \textwidth ]{figure/IR2/shock_yield_daySep1998.jpeg}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\caption{Change in Feb 1999}
\includegraphics[width = \textwidth ]{figure/IR2/shock_yield_dayFeb1999.jpeg}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\caption{Change in Jan 2007}
\includegraphics[width = \textwidth ]{figure/IR2/shock_yield_dayJan2007.jpeg}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\caption{Inflation IRF (Sep 1998)}\label{figIRd}
\includegraphics[width = \textwidth ]{figure/IR2/IRF_SChInv_ShockSep1998.jpeg}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\caption{Inflation IRF (Feb 1999)}\label{figIRe}
\includegraphics[width = \textwidth ]{figure/IR2/IRF_SChInv_ShockFeb1999.jpeg}
\end{subfigure}
\begin{subfigure}{.32\textwidth}\caption{Inflation IRF (Jan 2007)}\label{figIRf}
\includegraphics[width = \textwidth ]{figure/IR2/IRF_SChInv_ShockJan2007.jpeg}
\end{subfigure}
\flushleft{\scriptsize{Notes: The figure reports the change in the yield curves induced by shocks in 9/1998, 2/1999, and 1/2007, specified in \cite{IR2021}, and inflation response estimates. The top panel shows the yield curves before (grey) and after (black) the shock is introduced. In the bottom panel, the solid lines represent the pointwise estimates, while the dotted gray lines indicate their 90\% confidence interval.}}
\label{fig:enter-label}
\end{figure}
The estimated impulse responses and shifts in the yield curve caused by the shocks are reported in Figure~\ref{fig: emp: ir}.
The first observation is that the response seems to depend on the shape of the shock. For example, in September 1998 and January 2007, the yield curves experience decreases in overall maturities. Their effects reported in Figures~\ref{figIRd} and \ref{figIRf} align with economic theory, which suggests that decreases in interest rates for all maturities lead to higher inflation growth. On the other hand, if the shock raises interest rates in general (e.g., February 1999), inflation tends to respond negatively, with the effect becoming more pronounced over a longer horizon. Although such a similar observation can be found in \cite{IR2021}, it is worth noting that the results in Figure~\ref{fig: emp: ir} are not directly comparable with those reported in \cite{IR2021}. This is because of the nonparametric nature of our estimation approach, which does not require a parametric restriction on either the parameter of interest or the functional variables. Therefore, our estimator should be understood as a complement to the estimator considered by the authors.
\section{Simulation\label{sec:sim}}
In this section, we study finite sample performance of estimators studied in Section~\ref{sec:emp} using Monte Carlo experiments. Throughout the section, we consider the case with $T=250$ and $T=500$, and
the total number of replications is set to 1,000.
\subsection{Experiment A: Linear projection DGP\label{sec:simA}}
In our simulation, the variable $X_t $ is designed to mimic the economic sentiment quantile functions in Section~\ref{sec:emp1}. To this end, let $X_t = \sum_{j=1} ^{31} x_{j,t}\xi_j$ where $\xi_j$ denotes the $j$-th eigenvector of $X_t$'s variance and the $j$-th coordinate process $x_{j,t} $ is given by $ \langle X_t, \xi_j\rangle$. Assume that $x_{j,t}$ follows AR(1) such that \begin{equation*}
x_{j,t} = {\alpha}_{j} ^x x_{j,t-1} + c_{e,j} { \sigma}_{j} e_{j,t},
\end{equation*} for $j=1, \ldots, 31$, where $e_{j,t}\sim_{iid} N(0,1)$. The parameter values $\alpha_j ^x$ and $\sigma_j$ are set to the estimates obtained from the economic sentiment quantile curves, with $c_{e,j}$ set to one for all $j$. Later in the section, we keep or exaggerate the variance of each coordinate process by setting $c_{e,j} = c_{e,-1} 1\{ j\neq 1\} + 1\{j=1\}$ and considering two different values of $c_{e,-1}$: (i) $c_{e,-1} = 1$ and (ii) $c_{e,-1} = 4$. In the latter case, the ordered eigenvalues (from largest to smallest) are scaled up, except for the first one, without altering their order. Hence, $\lambda_j$ exhibits a significantly slower diminishing tendency for $j$ that is not large, compared to the former case.
We assume $\{ y_t, X_t \}_{t=1} ^T$ follows the DGP such that
\begin{equation*}
y_{t+h} = \alpha_{h} y_{t} + \int_0 ^1 X_t(s) \beta_h (s) ds + c_u \sigma_{u,h} u_t, \quad \text{for }h=1,\ldots, H,
\end{equation*}
where $u_{t}\sim_{iid} N(0,1)$ and $\beta_h(\cdot) = \sum_{j=1} ^{J} \beta_{j,h} \xi_j(\cdot)$. In our empirical data, the first two coordinate processes ($x_{1,t}$ and $x_{2,t}$) explain more than 99\% of the variations. Thus, for each $h$, the parameters $\{\alpha_h, \beta_{1,h}, \beta_{2,h}, \sigma_{u,h}\}$ are replaced by the estimates from the empirical data when $c_u=1$ and $J=2$. In the simulation, we increase $J$ to 31 and let the remaining coefficients $\{\beta_{j+2,h}\}_{j=1} ^{29}$ be given by $ \min \{ |\beta_{1,h}|, |\beta_{2,h}| \} \times 0.7^j$. This simulation design allows us to avoid the inverse problem in obtaining the coefficient estimates while keeping most variations in $X_t$, so that realizations from the DGP can mimic the empirical data. In the simulation, the constant $c_u$ is set to 0.5. We note that the coefficients in this section should not be interpreted as the coefficient estimates in Section~\ref{sec:emp1}.
We consider three estimators: SCInv, SCInv$_1$ and PCA-FR.
SCInv and SCInv$_1$ are differentiated in the choice of the regularization parameter. SCInv is computed as in Section~\ref{sec:emp}. The second estimator SCInv$_1$ and also \citepos{SHIN2009} PCA-FR are computed by using the AIC criterion detailed in \cite{SHIN2009}.
\begin{table}[tbp]
\caption{ Relative Bias and Variance Estimates (Experiment 1) } \label{tab: simA-1}
\vskip -8pt
\small
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lllccccccccc}
\toprule & & & \multicolumn{3}{c}{$h= 1$} & \multicolumn{3}{c}{ $h=3$} & \multicolumn{3}{c}{$h= 5$} \\\cmidrule{4-6}\cmidrule{7-9}\cmidrule{10-12}
T & $c_{e,-1}$ & & SCInv &SCInv$_1$ & PCA-FR & SCInv & SCInv$_1$ & PCA-FR & SCInv & SCInv$_1$ & PCA-FR \\ \midrule
\multirow{4}{*}{250} & \multirow{2}{*}{4} & Bias & 1.00 & 1.07 & 1.07 & 1.00 & 2.03 & 2.03 & 1.00 & 1.53 & 1.53 \\
& & Var & 1.00 & 1.80 & 1.79 & 1.00 & 2.15 & 2.15 & 1.00 & 0.30 & 0.30 \\ \cmidrule{2-12}
& \multirow{2}{*}{1} & Bias & 1.00 & 3.13 & 3.13 & 1.00 & 4.44 & 4.44 & 1.00 & 1.48 & 1.48 \\
& & Var & 1.00 & 1.83 & 1.83 & 1.00 & 1.24 & 1.24 & 1.00 & 0.23 & 0.23 \\ \midrule
\multirow{4}{*}{500} & \multirow{2}{*}{4} & Bias & 1.00 & 1.11 & 1.11 & 1.00 & 1.20 & 1.20 & 1.00 & 1.64 & 1.64 \\
& & Var & 1.00 & 0.90 & 0.90 & 1.00 & 1.68 & 1.68 & 1.00 & 0.27 & 0.27 \\ \cmidrule{2-12}
& \multirow{2}{*}{1} & Bias & 1.00 & 1.82 & 1.82 & 1.00 & 4.27 & 4.27 & 1.00 & 1.49 & 1.49 \\
& & Var & 1.00 & 2.57 & 2.57 & 1.00 & 2.11 & 2.11 & 1.00 & 0.15 & 0.15 \\
\bottomrule
\end{tabular*}
\flushleft{\scriptsize{Notes: Based on 1,000 replications. The table reports the $L_2$ norm of the bias (Bias) and the variance (Var) relative to those of SCInv. As $h$ increases, the estimators exhibit smaller bias and variance, which is expected, as the function values of $\beta_h$, computed from empirical data, tend to diminish in scale as $t$ increases.}}
\end{table}
\begin{table}[h!]
\caption{Coverage Probability (Experiment 1) } \label{tab: simA-2}
\vskip -8pt
\small
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcccccc}
\toprule & \multicolumn{3}{c}{$T= 250$} & \multicolumn{3}{c}{ $T=500$} \\\cmidrule{2-4}\cmidrule{5-7}
$c_{e,-1} \backslash h $ & 1&3&5&1&3&5 \\ \midrule
4 & 0.95 & 0.93 & 0.93 & 0.95 & 0.94 & 0.94 \\
1 & 0.94 & 0.93 & 0.93 & 0.95 & 0.94 & 0.94 \\
\bottomrule
\end{tabular*}
\flushleft{\scriptsize{Notes: Based on 1,000 replications. The 95\% coverage probabilities are computed using SCInv and Theorem~\ref{thm2}. }}
\end{table}
The estimation results are summarized in Table~\ref{tab: simA-1} for three forecasting horizons: $h=1,3,5$. We report the bias (Bias) and variance (Var) of each estimator relative to those computed using SCInv. Overall, the proposed estimator based on the naive choice of $\rho$ (SCInv) tends to produce the smallest bias and variance, particularly when $c_{e,-1}$ is large, which is related to the relative decay rate of $\lambda_j$.
Meanwhile, the two estimators based on the same regularization parameter, SCInv$_{1}$ and PCA-FR, perform similarly.
We then study Theorem~\ref{thm2} by the means of local asymptotic normality to study the impulse responses $\{\langle \beta_h, \zeta \rangle; h=1,3,5 \}$, where $\zeta$ is set to a constant function for simplicity. Table~\ref{tab: simA-2} reports 95\% coverage probabilities computed with SCInv for each forecasting horizon. Overall, the coverage probability is very close to the nominal level in both sample sizes. As discussed in Section~\ref{sec:est}, to the best of the authors' knowledge, the assumptions on $\rho$, $\varsigma$ and $\delta$ in Theorem~\ref{thm2} are not directly testable in practice, and thus the asymptotic bias terms $\widehat \Theta_{2A}$ and $\widehat \Theta_{2B}$ may not be asymptotically negligible. Nevertheless, the simulation results reported in Table~\ref{tab: simA-2} suggest that those asymptotic bias would be small and thus do not distort testing results based on the asymptotic approximation in Theorem~\ref{thm2}.
Figure~\ref{fig: simA-1} reports the pointwise estimates of the functional coefficient at $h=1$ and the impulse response estimates for $h=1,\ldots, 6$ when $c_{e,-1}=1$. The simulation results for $c_{e,-1}=4$ are similar; thus, we omit the figures to save space. The reported coefficient estimates are obtained by averaging the estimates across 1,000 simulations. Consistent with the observations in the tables, our estimator produces both functional coefficient estimates and the pointwise impulse response estimates close to the true values, demonstrating the practical applicability of our approach.
\begin{figure}
\centering\caption{ $\beta_1$ and $\{\langle \zeta, \beta_h\rangle \}_{h=1} ^{6}$ Estimates (Experiment 1; $T=250$) \label{fig: simA-1}}
\begin{subfigure}{.45\linewidth}\subcaption{Functional coefficient estimates at $h=1$\label{fig: simA-1a}}
\includegraphics[width = \textwidth]{figure/Estcurve\_Lsigfac100lnfac2Ctmp60Cexp70ModelEx.jpeg}
\end{subfigure}
\begin{subfigure}{.45\linewidth} \subcaption{Impulse response estimates for $h=1,\ldots , 6$.\label{fig: simA-1b}}
\includegraphics[width = \textwidth]{figure/IRF\_Lsigfac100lnfac2Ctmp60Cexp70ModelSC.jpeg}
\end{subfigure}
\flushleft{\scriptsize{Notes: Averaged across 1,000 replications. Figure~\ref{fig: simA-1a} reports the functional coefficient estimate for $h=1$ computed with SCInv, SCInv1 and PCA-FR. Figure~\ref{fig: simA-1b} reports $\langle \zeta, \beta_h \rangle $ (blue), its pointwise estimates based on the SCInv (dashed) and its 90\% confidence interval (shaded area) for $h=1,\ldots, 6$, when $\zeta$ is specified to the constant function. $c_{e,-1} = 1$. }}
\end{figure}
\subsection{Experiment B: SVAR-based DGP\label{sec:simB}}
Lastly, we examine the performance of our estimator as an estimator of the SIRF suggested in Propositions~\ref{prop: svar: identification} and \ref{prop: svar: identification: a}, using a small-scale Monte Carlo simulation based on the following SVAR(1) process:\begin{align*}
y_t &= \alpha_{11} y_{t-1} + \alpha_{12} x_{1,t-1} + \alpha_{13} x_{2,t-1} + u_{1,t},\\
x_{1,t} & = \beta_{1} y_{t} + \alpha_{21} y_{t-1} + \alpha_{22} x_{1,t-1} + \alpha_{23} x_{2,t-1} + u_{2,t},\\
x_{2,t} & = \beta_{2} y_{t} + \alpha_{31} y_{t-1} + \alpha_{32} x_{1,t-1} + \alpha_{33} x_{2,t-1} + u_{3,t},
\end{align*} where $(u_{1,t}, u_{2,t}, u_{3,t} )' \sim_{iid} \mathcal N(0, \text{diag}(\sigma_1^2, \sigma_2^2, \sigma_3^2))$ and $x_{j,t}$ denotes the $j$-th coordinate process of economic sentiment quantiles. The functional predictor $X_t$ is assumed to be generated as follows. \begin{equation}
X_t = \sum_{j=1} ^{2} \mathbf{x}_{t, 2 }\xi_j + \sum_{j= 3} ^{31} u_{j+1, t} \xi_j.
\end{equation}
where $ u_{j+1, t} \sim_{iid} N( {0}, \sigma_j ^2 ) $ across $j$ and $t$, and $\sigma_j ^2 = \sigma_3 0.8^j$ for $j = 3, \ldots, 31$, and $\xi_j$ are the eigenvectors of $X_t$'s covariance. The parameters and $\{\xi_j\}_{j=1} ^{31}$ are estimated from empirical data as in Section \ref{sec:simA}. As before, this setup is designed to keep the structure of $X_t$ whose largest variations are determined by the first two coordinate processes. Then, in the simulation, we replace the variance structure of $\mathbf {u}_t= (u_{1,t},u_{2,t}, \ldots , u_{32,t})'$ with $c_{\ast}\text{diag}(c_1 \sigma_1^2, \sigma_2^2,\sigma_3^2, \ldots, \sigma_{32}^2 )$, where $c_{\ast}$ is the normalizing constant that allows us to keep the norm of the variance of $\mathbf {u}_t$ to be equal to one. The other constant $c_1$ determines the relative magnitude of the structural error associated with $y_t$, and thus it is inversely related to the signal-to-noise ratio. We consider three values of $c_1$: 1, 0.5, and 0.2. Along with our estimators, the SCInv and the SCInv$_{1}$, we consider the PCA-SVAR for comparison. As the PCA-SVAR produces the same estimation results with the PCA-FR in this setup, we omit the latter.
\begin{table}[tbp]
\caption{Relative Bias and Variance Estimates (Experiment 2) \label{tab: simB-1}}
\vskip -8pt
\small
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}llccccccccc}
\toprule
& $c_{1}$ & \multicolumn{3}{c}{ 1} & \multicolumn{3}{c}{ 0.5} & \multicolumn{3}{c}{ 0.2} \\\cmidrule{3-5}\cmidrule{6-8}\cmidrule{9-11}
T & & SCInv & SCInv$_{1}$ & SVAR & SCInv & SCInv$_{1}$ & SVAR & SCInv & SCInv$_{1}$ & SVAR \\\midrule
\multirow{2}{*}{250} & Bias & 1.00 & 8.98 & 8.98 & 1.00 & 10.27 & 10.25 & 1.00 & 9.35 & 9.32 \\
& Var & 1.00 & 0.46 & 0.46 & 1.00 & 0.74 & 0.74 & 1.00 & 1.24 & 1.24 \\ \midrule
\multirow{2}{*}{500}& Bias & 1.00 & 29.57 & 29.53 & 1.00 & 25.11 & 25.06 & 1.00 & 13.61 & 13.56 \\
& Var & 1.00 & 0.65 & 0.65 & 1.00 & 1.00 & 1.00 & 1.00 & 1.17 & 1.17 \\ \bottomrule
\end{tabular*}
\flushleft{\scriptsize{Notes: Based on 1,000 replications. The table reports the $L_2$ norm of the bias (Bias) and the variance (Var) of each estimator relative to those of SCInv. }}
\end{table}
\begin{figure}[h!]
\caption{Estimated Functional Coefficient (Experiment 2; $T=250$) \label{fig: simB-1}}
\includegraphics[width = \textwidth, height= .2\textheight]{figure/Estcurve_Lsigfac45lnfac2E180SVar.jpeg}
\flushleft{\scriptsize{Notes: Averaged across 1,000 replications. The figure reports the true functional coefficient (black) and its estimates based on SCInv (dotted) and SVAR (dashed) when $c_1$ is 1 (left), 0.5 (middle), and 0.2 (right). $h=1$. }}
\end{figure}
Table~\ref{tab: simB-1} summarizes simulation results. Overall, the SCInv produces the smallest bias at the cost of a relatively large variance. As observed previously, the other two estimators, the SCInv$_1$ and the PCA-SVAR, report similar estimation results with each other, regardless of the value of $c_1$ and the sample size, as the difference between the two regularization schemes is not significant in this simulation design.
Figure~\ref{fig: simB-1} reports the true functional coefficient and its estimates (averaged across 1,000 replications) for each value of $c_1$ when $h=1$. In the figure, it is noticeable that the PCA-SVAR estimator tends to get closer to the true functional coefficient as $c_1$ decreases. On the other hand, our estimator, the SCInv (dotted), produces very robust estimation results close to the true regardless of the value of $c_1$.
\section{Conclusion \label{sec:con}}
In this paper, we study impulse response analysis with functional predictors and other scalar-valued covariates and propose new estimation and inference methodologies. We show that the proposed estimator allows for an interesting interpretation as a structural impulse response in some special cases. In our empirical application, we study how economic variables respond to a certain distributional shock on sentiment. The results are consistent with existing observations and show negative (resp.\ positive) responses of economic growth variables when sentiment distributions shift to the left (resp.\ right). Monte Carlo simulation results also confirm our theoretical findings.
\makeatletter
\def\@seccntformat#1{
\csname the#1\endcsname.\quad
}
\makeatother
\newpage