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.
57,536 characters
\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}
\begin{center}
{\Large \bf Time-Varying Multivariate Causal Processes}
\bigskip
$^{\ast}${\sc Jiti Gao} and $^{\ast}${\sc Bin Peng} and $^{\dag}${\sc Wei Biao Wu} and $^{\ast}${\sc Yayi Yan}
\medskip
$^{\ast}$Department of Econometrics and Business Statistics, Monash University,\\
and $^{\dag}$Department of Statistics, University of Chicago
\end{center}
\medskip
\begin{abstract}
In this paper, we consider a wide class of time-varying multivariate causal processes which nests many classic and new examples as special cases. We first prove the existence of a weakly dependent stationary approximation for our model which is the foundation to initiate the theoretical development. Afterwards, we consider the QMLE estimation approach, and provide both point-wise and simultaneous inferences on the coefficient functions. In addition, we demonstrate the theoretical findings through both simulated and real data examples. In particular, we show the empirical relevance of our study using an application to evaluate the conditional correlations between the stock markets of China and U.S. We find that the interdependence between the two stock markets is increasing over time.
\medskip
\noindent
{\it Keywords:} Local Linear Quasi-Maximum Likelihood Estimation; Multivariate Causal Process; Simultaneous Confidence Interval.
\smallskip
\noindent{\it JEL Classification:} C14, C32, G15.
\end{abstract}
\vfill
\newpage
\spacingset{1.9}
\section{Introduction}\label{Sec1}
The family of vector autoregressive (VAR) models and the family of multivariate (G)ARCH models are among some of the most popular frameworks for modelling dynamic interactions of multiple variables. The VAR family usually captures the dynamic by imposing structures on the time series itself, while the (G)ARCH family imposes restrictions on the conditional second moments. We acknowledge the vast literature of both families, and have no intention to exhaust all relevant studies in this paper for the sake of space. We refer interested readers to \cite{stock2001vector} and \cite{bauwens2006multivariate} for excellent review on both families.
Although both families have rich literature on their own, to the best of the authors' knowledge not many works have been done to bridge them. Among limited attempts (e.g., \citealp{ling2003asymptotic,bardet2009asymptotic}), most (if not all) of these studies rely on the stationarity assumption. While the stationarity assumption comes in handy when deriving asymptotic properties, it may not be very realistic in practice (\citealp{preuss2015detection,chen2021inference}). For example, economic and financial data always include different macro shocks, as a consequence the behaviour can be quite volatile; the climate data may contain certain time trend which recently has attracted lots of attention due to greenhouse emission; etc. Anyway, certain nonstationarity may always occur.
To account for nonstationarity, locally stationary processes have received considerable attention since the seminal work of \cite{dahlhaus1996kullback}, \cite{dette2011measure}, \cite{zhang2012inference}, \cite{truquet2017parameter}, \cite{dahlhaus2019towards}, among others. In contrast to the unit root process, the locally stationary process nicely balances stationarity and nonstationarity by allowing for the simultaneous presence of both types of behaviours in one time series process. In a very recent paper, \cite{karmakar2021simultaneous} consider simultaneous inference for a general class of univariate $p$-Markov processes with time-varying coefficients, which covers several time-varying versions of the classical univariate models (e.g., AR, ARCH, AR-ARCH) as special cases. Despite its generality, their study still rules out the time-varying versions of some widely used models (e.g., ARMA, GARCH, ARMA-GARCH). Also, it is worth mentioning this line of research heavily focuses on univariate time series, which somewhat limits the popularity of locally stationary processes.
That said, it is reasonable to call for a framework which can marry the VAR family and the (G)ARCH family while allowing for nonstationarity. To provide a concrete example, consider a time-varying multivariate GARCH model, which can model the co-movements of financial returns. Detailed investigation on such a model can help answer research questions like (i). Is the volatility of a market leading the volatility of other markets? (ii) Whether the correlations between asset returns change over time? (iii). Are they increasing in the long run, perhaps because of the globalization of financial markets? These are of great practical importance for both investors and policymakers (\citealp{bauwens2006multivariate,diebold2009measuring}).
To allow for flexibility as much as possible from the modelling perspective, we consider a class of multivariate causal processes as follows:
\begin{eqnarray}\label{Eq2.1}
\mathbf{x}_t =\left\{\begin{array}{ll}
\bm{\mu}\left(\mathbf{x}_{t-1},\mathbf{x}_{t-2},\ldots;\bm{\theta}(\tau_t)\right) + \mathbf{H}\left(\mathbf{x}_{t-1},\mathbf{x}_{t-2},\ldots;\bm{\theta}(\tau_t)\right)\bm{\varepsilon}_t, & \text{for}\quad t=1,\ldots, T\\
\bm{\mu}\left(\mathbf{x}_{t-1},\mathbf{x}_{t-2},\ldots;\bm{\theta}(0)\right) + \mathbf{H}\left(\mathbf{x}_{t-1},\mathbf{x}_{t-2},\ldots;\bm{\theta}(0)\right)\bm{\varepsilon}_t & \text{for}\quad t\le 0
\end{array}\right. ,
\end{eqnarray}
where $\tau_t= t/T$, $\bm{\mu}\left(\cdot\right)$ is an $m$-dimensional random vector, $\mathbf{H}\left(\cdot\right)$ is an $m\times m$-dimensional random matrix, $\bm{\theta}(\tau)$ is a $d\times 1$ time-varying parameter of interest with each element belonging to $C^3[0,1]$, and $\{\bm{\varepsilon}_t\}$ is a sequence of independent and identically distributed (i.i.d.) random vectors. Note that the value of $d$ usually depends on the value of $m$, and the connection becomes clear once a specific model is considered. As far as we are concerned, both of $m$ and $d$ are fixed throughout the paper. Notably, both $\bm{\mu}(\cdot)$ and $ \mathbf{H}(\cdot)$ are known, and share the same unknown parameter $\bm{\theta}(\cdot)$. The setting for $t\le 0$ regulates the time series for the periods that we do not observe, which is commonly adopted when certain nonstationarity gets involved (e.g., \citealp{vogt2012nonparametric}). Essentially, it requires the initial time period does not have a diverging behaviour.
Before proceeding further, we provide two examples to briefly illustrate the rationality behind \eqref{Eq2.1}, and leave the detailed investigation on these examples to Section \ref{Sec2.4}. We refer interested readers to \cite{ling2003adaptive}, \cite{ling2003asymptotic} and \cite{bardet2009asymptotic} for extensive investigation on the parametric counterparts of these examples.
\noindent \textbf{Example 1}: Consider the time-varying VARMA($p,q$) model
\begin{equation}\label{Eq4.1}
\mathbf{x}_t = \mathbf{a}(\tau_t) + \sum_{j=1}^{p}\mathbf{A}_j(\tau_t)\mathbf{x}_{t-j} + \bm{\eta}_t + \sum_{j=1}^{q}\mathbf{B}_j(\tau_t)\bm{\eta}_{t-j} \quad \text{with}\quad \bm{\eta}_t = \bm{\omega}(\tau_t)\bm{\varepsilon}_t.
\end{equation}
It is not hard to show that \eqref{Eq4.1} admits a presentation in the form of \eqref{Eq2.1}, and
\begin{eqnarray}\label{Eq4.1.1}
\bm{\theta}(\tau)=\mathrm{vec}(\mathbf{a}(\tau),\mathbf{A}_1(\tau),\ldots,\mathbf{A}_p(\tau),\mathbf{B}_1(\tau),\ldots,\mathbf{B}_q(\tau),\bm{\Omega}(\tau)),
\end{eqnarray}
where $\bm{\Omega}(\cdot):=\bm{\omega}(\cdot)\bm{\omega}^\top(\cdot)$.
\smallskip
\noindent \textbf{Example 2}: Consider the time-varying multivariate GARCH($p,q$) model
\begin{eqnarray}\label{Eq4.4}
\mathbf{x}_t &=& \mathrm{diag} (h_{1,t}^{1/2},\ldots,h_{m,t}^{1/2} ) \bm{\eta}_t,\nonumber \\
\mathbf{h}_t &=& \mathbf{c}_0(\tau_t) + \sum_{j=1}^{p} \mathbf{C}_j(\tau_t) \left(\mathbf{x}_{t-j}\odot\mathbf{x}_{t-j}\right) + \sum_{j=1}^{q} \mathbf{D}_j(\tau_t)\mathbf{h}_{t-j},
\end{eqnarray}
where $h_{j,t}$ stands for the $j^{th}$ element of $\mathbf{h}_t$, and $\bm{\eta}_t = \bm{\Omega}^{1/2}(\tau_t) \bm{\varepsilon}_t$. The model \eqref{Eq4.4} generalizes the models of \cite{bollerslev1990modelling} and \cite{jeantheau1998strong}. Similar to Example 1, we show that \eqref{Eq4.4} admits a representation in the form of \eqref{Eq2.1}, and
\begin{eqnarray}\label{Eq4.4.1}
\bm{\theta}(\tau)=\mathrm{vec}(\mathbf{c}_0(\tau),\mathbf{C}_1(\tau),\ldots,\mathbf{C}_p(\tau),\mathbf{D}_1(\tau),\ldots,\mathbf{D}_q(\tau),\bm{\Omega}(\tau)).
\end{eqnarray}
\smallskip
In view of the development of Example 1 and Example 2 in Section \ref{Sec2.4}, one may further show the time-varying counterparts of the parametric models mentioned in \cite{bardet2009asymptotic} are also covered by \eqref{Eq2.1}. To this end, we argue that \eqref{Eq2.1} does not only allows for nonstationarity and conditional heteroskedasticity, but also provides sufficient flexibility to cover many well adopted models in the literature.
In this paper, our contributions are in the following four-fold: (1). we consider a wide class of time-varying multivariate causal processes which nests many classic and new examples as special cases; (2). we prove the existence of a weakly dependent stationary approximation for the model \eqref{Eq2.1} at any given time of interest (i.e., $\forall\tau\in[0,1]$), which is the foundation in order to establish asymptotic properties associated with the model; (3). we establish the estimation theory, and provide both point-wise and simultaneous inferences on the coefficient functions of which both are important for practical works (\citealp{zhou2010simultaneous}); (4). we demonstrate the theoretical findings through both simulated and real data examples.
The paper is organized as follows. Section \ref{Sec2} presents the theoretical findings associated with the stationary approximation, estimation and inferences. In Section \ref{Sec3}, we conduct extensive simulation studies to examine the theoretical findings, and further investigate the time-varying conditional correlations between the Chinese and U.S. Stock market. Section \ref{Sec4} concludes. Due to space limit, we give the proofs of the main results to the online appendices of the paper.
Before proceeding further, it is convenient to introduce some notation: the symbol $|\cdot|$ denotes the Euclidean norm of a vector or the spectral norm for a matrix; $\|\mathbf{v}\|_q:=\left(E|\mathbf{v}|^q\right)^{1/q}$ and $\|\cdot\|:=\|\cdot\|_2$ for short; $\otimes$ denotes the Kronecker product; $\odot$ denotes the Hadamard product; $\mathbf{I}_a$ stands for an $a\times a$ identity matrix; $\mathbf{0}_{a\times b}$ stands for an $a\times b$ matrix of zeros, and we write $\mathbf{0}_a$ for short when $a=b$; for a function $g(w)$, let $g^{(j)}(w)$ be the $j^{th}$ derivative of $g(w)$, where $j\ge 0$ and $g^{(0)}(w) \equiv g(w)$; $K_h(\cdot) =K(\cdot/h)/h$, where $K(\cdot)$ and $h$ stand for a nonparametric kernel function and a bandwidth respectively; let $\tilde{c}_k =\int_{-1}^{1} u^k K(u) \mathrm{d}u$ and $\tilde{v}_k= \int_{-1}^{1} u^k K^2(u) \mathrm{d}u$ for integer $k\ge 0$; $\mathrm{diag}(\mathbf{a})$ is a diagonal matrix with the vector $\mathbf{a}$ on its main diagonal, while $\mathrm{diag}(\mathbf{A})$ creates a vector from the diagonal of matrix $\mathbf{A}$; finally, let $\to_P$ and $\to_D$ denote convergence in probability and convergence in distribution, respectively.
\section{Estimation and Asymptotics}\label{Sec2}
In this section, we first prove the existence of a weakly dependent stationary approximation for the model \eqref{Eq2.1} in Section \ref{Sec2.1}; we then provide the estimation approach using the local linear quasi-maximum-likelihood estimation and establish the asymptotic properties of the proposed estimator in Section \ref{Sec2.2}; Section \ref{Sec2.3} provides results on both point-wise and simultaneous inferences; Section \ref{Sec2.4} gives some detailed examples to justify the usefulness of our study.
\subsection{Stationary Approximation}\label{Sec2.1}
To study \eqref{Eq2.1}, the first challenge lies in the fact that the model may not be stationary. Therefore, for $\forall\tau\in [0,1]$, we initial our analysis by finding a stationary approximation for each $\mathbf{x}_t$ with $t\ge 1$. By doing so, we are able to measure the weak dependence of $\{\mathbf{x}_t \}$ using the nonlinear system theory in \cite{wu2005nonlinear}, which then provides us a framework to derive the asymptotic properties accordingly.
To be clear on the dependence measure, consider an example in which $\mathbf{e}_t$ is a stationary process, and admits a causal representation $\mathbf{e}_t = \mathbf{J}(\bm{\varepsilon}_t,\bm{\varepsilon}_{t-1},\ldots)$ with $\mathbf{J}(\cdot)$ being a measurable function. See \cite{tong1990non} for discussion on nonlinear time series of this kind. For $k\geq 0$, we define the following dependence measure:
\begin{eqnarray}\label{DefDM}
\delta_{r}^{\mathbf{e}}(k)=\left\|\mathbf{J}(\bm{\varepsilon}_k,\bm{\varepsilon}_{k-1},\ldots\bm{\varepsilon}_1,\bm{\varepsilon}_{0},\bm{\varepsilon}_{-1},\ldots)-\mathbf{J}(\bm{\varepsilon}_k,\ldots, \bm{\varepsilon}_1,\bm{\varepsilon}_{0}^*,\bm{\varepsilon}_{-1},\ldots)\right\|_r,
\end{eqnarray}
where $\bm{\varepsilon}_0^*$ is an independent copy of $\{\bm{\varepsilon}_j\}$. Being able to measure the time series dependence such as \eqref{DefDM} is the starting point for time series analyses.
We now introduce some basic assumptions.
\begin{assumption} \label{Ass1}
\item
\begin{enumerate}
\item $\{\bm{\varepsilon}_t\}$ is a sequence of i.i.d. random vectors with $E(\bm{\varepsilon}_1)=\mathbf{0}$, $E (\bm{\varepsilon}_1\bm{\varepsilon}_1^\top )=\mathbf{I}_m$, and $ \|\bm{\varepsilon}_1 \|_r<\infty$ for some $r\ge 2$.
\item For $\forall\mathbf{z},\mathbf{z}^\prime \in (\mathbb{R}^m)^\infty$ and $\forall\bm\vartheta\in \mathbb{R}^d$, there exist nonnegative sequences $\{\alpha_j(\bm{\bm\vartheta})\}_{j=1}^{\infty}$ and $\{\beta_j(\bm{\bm\vartheta})\}_{j=1}^{\infty}$ such that
\begin{eqnarray*}
&&|\bm{\mu}(\mathbf{z};\bm{\bm\vartheta})-\bm{\mu}(\mathbf{z}^\prime;\bm{\bm\vartheta})| \leq \sum_{j=1}^{\infty}\alpha_j(\bm{\bm\vartheta}) |\mathbf{z}_j-\mathbf{z}_j^\prime|,\\
&&|\mathbf{H}(\mathbf{z};\bm{\bm\vartheta})-\mathbf{H}(\mathbf{z}^\prime;\bm{\bm\vartheta})| \leq \sum_{j=1}^{\infty}\beta_j(\bm{\bm\vartheta})|\mathbf{z}_j-\mathbf{z}_j^\prime|,
\end{eqnarray*}
where $\mathbf{z}_j$ and $\mathbf{z}_j^\prime$ are the $j^{th}$ columns of $\mathbf{z}$ and $\mathbf{z}^\prime$ respectively.
\item For $\forall \tau\in [0,1]$, $\bm{\theta}(\tau)$ lies in the interior of $\bm{\Theta}_r$, where
\begin{eqnarray*}
\bm{\Theta}_r:=\left\{\bm{\vartheta}\in \bm{\Theta} \mid \sum_{j=1}^{\infty}\alpha_j(\bm{\vartheta}) + \left\|\bm{\varepsilon}_1\right\|_r \sum_{j=1}^{\infty}\beta_j(\bm{\vartheta}) < 1 \right\}
\end{eqnarray*}
and $\bm{\Theta}$ is a compact set of $\mathbb{R}^d$.
\end{enumerate}
\end{assumption}
Assumption \ref{Ass1}.1 is standard when studying dynamic time series model (\citealp{lutkepohl2005new}). In Assumption \ref{Ass1}.2, $\bm\vartheta$ is a generic $d\times 1$ vector, and has the same length as $\bm \theta(\cdot)$. This assumption imposes Lipschitz-type conditions on $\bm{\mu}(\cdot)$ and $\mathbf{H}(\cdot)$, which are rather minor, and can be easily fulfilled by a variety of models such as those mentioned in Section \ref{Sec1}. See Propositions \ref{proposition4.1}-\ref{proposition4.2} below for details. Assumption \ref{Ass1}.3 does not only guarantee a stationary approximation for each $\mathbf{x}_t$, but also ensures the approximated process has some proper moments. Similar conditions have also been adopted in \cite{bardet2009asymptotic}.
With these conditions in hand, we present the following proposition which facilitates the development in what follows.
\begin{proposition}\label{proposition2.1}
Let Assumption \ref{Ass1} hold. For any $\tau \in [0,1]$, there exists a stationary process
\begin{eqnarray*}
\widetilde{\mathbf{x}}_t(\tau) = \bm{\mu}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right) + \mathbf{H}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right)\bm{\varepsilon}_t
\end{eqnarray*}
such that
\begin{enumerate}
\item $\sup_{\tau\in[0,1]}\left\|\widetilde{\mathbf{x}}_t(\tau)\right\|_r < \infty$,
\item $\delta_{r}^{\widetilde{\mathbf{x}}(\tau)}(k) \leq O(1)\inf_{1\leq p \leq k}\{\rho(\tau)^{k/p} + \sum_{j=p+1}^{\infty}\left[\alpha_j(\bm{\theta}(\tau)) + \beta_j(\bm{\theta}(\tau))\right]\} \to 0$ as $k\to \infty$,
\end{enumerate}
where $\rho(\tau) := \sum_{j=1}^{\infty}\alpha_j(\bm{\theta}(\tau)) + \left\|\bm{\varepsilon}_1\right\|_r \sum_{j=1}^{\infty}\beta_j(\bm{\theta}(\tau))$.
\end{proposition}
It is worth mentioning that for a univariate $p$-Markov process
\begin{eqnarray*}
\widetilde{x}_{p,t}(\tau) = \mu\left(\widetilde{x}_{t-1}(\tau),\ldots,\widetilde{x}_{t-p}(\tau);\bm{\theta}(\tau)\right) + H\left(\widetilde{x}_{t-1}(\tau),\ldots,\widetilde{x}_{t-p}(\tau);\bm{\theta}(\tau)\right)\varepsilon_t,
\end{eqnarray*}
\cite{karmakar2021simultaneous} show that there exists $0 < \rho < 1$ such that $\sup_{\tau \in[0,1]}\delta_{r}^{\widetilde{x}_p(\tau)}(k) =O(\rho^k)$ based on the development of \cite{wu2004limit}. From a methodological viewpoint, we give a set of new proofs which allow us to measure the dependence of multivariate causal processes with infinity memory. The term $\sum_{j=p+1}^{\infty}\left[\alpha_j(\bm{\theta}(\tau)) + \beta_j(\bm{\theta}(\tau))\right]$ in the second result of Proposition \ref{proposition2.1} arises due to the infinity memory structure of $\widetilde{\mathbf{x}}_t(\tau)$. Thus, the dependence $\delta_{r}^{\widetilde{\mathbf{x}}(\tau)}(k)$ relies on the choice of $p$ and the decay rates of the coefficients $\alpha_j(\bm{\theta}(\tau))$ and $\beta_j(\bm{\theta}(\tau))$.
To ensure $\widetilde{\mathbf{x}}_t(\tau) $ can approximate $\mathbf{x}_t$ reasonably well, we impose more structure below.
\begin{assumption}\label{Ass2}
\item
\begin{enumerate}
\item There exists a nonnegative sequence $\{\chi_j\}$ with $\sum_{j=1}^\infty\chi_j < \infty$ such that for $\forall\mathbf{z}\in (\mathbb{R}^m)^\infty$ and $\forall\bm{\vartheta},\bm{\vartheta}^\prime\in\bm{\Theta}_r$
\begin{eqnarray*}
&&|\bm{\mu}(\mathbf{z};\bm{\vartheta})-\bm{\mu}(\mathbf{z};\bm{\vartheta}^\prime)|+ |\mathbf{H}(\mathbf{z};\bm{\vartheta})-\mathbf{H}(\mathbf{z};\bm{\vartheta}^\prime)| \leq |\bm{\vartheta}-\bm{\vartheta}^\prime| \sum_{j=1}^{\infty}\chi_j|\mathbf{z}_j| .
\end{eqnarray*}
\item Let $\sup_{\tau \in [0,1]} \alpha_j(\bm{\theta}(\tau)) = O(j^{-(2+s)})$ and $\sup_{\tau \in [0,1]} \beta_j(\bm{\theta}(\tau)) = O(j^{-(2+s)})$ for some $s > 0$.
\end{enumerate}
\end{assumption}
Assumption \ref{Ass2}.1 imposes another Lipschitz-type condition with respect to the parameter space. Assumption \ref{Ass2}.2 further restricts the decay rates of $\alpha_j(\bm{\theta}(\tau))$ and $\beta_j(\bm{\theta}(\tau))$.
Using Assumptions \ref{Ass1}--\ref{Ass2}, we can measure the distance between $\widetilde{\mathbf{x}}_t(\tau) $ and $\mathbf{x}_t$ as follows.
\begin{proposition}\label{proposition2.2}
Suppose Assumptions \ref{Ass1}--\ref{Ass2} hold. Then
\begin{enumerate}
\item $\left\|\widetilde{\mathbf{x}}_1(\tau)-\widetilde{\mathbf{x}}_1(\tau^\prime)\right\|_r = O(|\tau - \tau^\prime|)$ for $\forall\tau,\tau^\prime \in [0,1]$,
\item $\max_{t\ge 1}\left\|\mathbf{x}_t-\widetilde{\mathbf{x}}_t(\tau_t)\right\|_r = O(T^{-1})$.
\end{enumerate}
\end{proposition}
We can consider Proposition \ref{proposition2.2} as the stochastic version of the H\"older continuity. Having established the stationary approximation in Proposition \ref{proposition2.2}, we move on to investigate the estimation theory in the next subsection.
\subsection{Estimation}\label{Sec2.2}
We point out a few facts to facilitate the setup of the likelihood function. First, let $\mathbf{z}_{t} =(\mathbf{x}_t,\mathbf{x}_{t-1},\ldots)$ include all the information of $\mathbf{x}_t$ up to the time period $t$. However, in practice, our observation on $\mathbf{x}_t$ only starting from $t=1$, so we have to work with the truncated version of $\mathbf{z}_{t}$ for each $t\ge 1$:
\begin{eqnarray}
\mathbf{z}_{t}^c =(\mathbf{x}_t,\ldots, \mathbf{x}_1, \mathbf{0},\ldots).
\end{eqnarray}
Second, we note that when $\tau_t$ is sufficiently close to $\tau$,
\begin{eqnarray}
\bm{\theta}(\tau_t) \approx \bm{\theta}(\tau ) +h\bm{\theta}^{(1)}(\tau) \cdot \frac{\tau_t-\tau}{h}.
\end{eqnarray}
Therefore, we are able to parametrize $\bm{\theta}(\cdot)$, and consider the maximum-likelihood estimation for each given $\tau$. Finally, since $\bm{\varepsilon}_t$ may not be normally distributed, we consider the local linear quasi-maximum-likelihood estimation (QMLE) method.
Thus, our likelihood function is specified as follows:
\begin{equation}\label{Eq2.3}
\mathcal{L}_\tau(\bm{\eta}_1,\bm{\eta}_2)=\frac{1}{T}\sum_{t=1}^{T}\mathcal{l}(\mathbf{x}_t,\mathbf{z}_{t-1}^c;\bm{\eta}_1+\bm{\eta}_2\cdot(\tau_t-\tau)/h)K_h(\tau_t-\tau),
\end{equation}
where
\begin{eqnarray*}
\mathcal{l}(\mathbf{x}_t,\mathbf{z}_{t-1}^c;\bm\vartheta) &=&- \frac{1}{2}(\mathbf{x}_t - \bm{\mu}(\mathbf{z}_{t-1}^c;\bm\vartheta))^\top\left(\mathbf{H}(\mathbf{z}_{t-1}^c;\bm\vartheta)\mathbf{H}(\mathbf{z}_{t-1}^c;\bm\vartheta)^\top\right)^{-1}(\mathbf{x}_t - \bm{\mu}(\mathbf{z}_{t-1}^c;\bm\vartheta))\\
&&-\frac{1}{2}\log\det\left(\mathbf{H}(\mathbf{z}_{t-1}^c;\bm\vartheta)\mathbf{H}(\mathbf{z}_{t-1}^c;\bm\vartheta)^\top\right).
\end{eqnarray*}
Accordingly, for $\forall \tau$, $(\bm{\theta}(\tau),h\bm{\theta}^{(1)}(\tau))$ is estimated by
\begin{equation}\label{Eq2.4}
(\widehat{\bm{\theta}} (\tau),\widehat{\bm{\theta}}^\star(\tau)) = \operatorname*{\arg\!\max}_{(\bm{\eta}_1,\bm{\eta}_2)\in\mathbf{E}_T(r)}\mathcal{L}_\tau (\bm{\eta}_1,\bm{\eta}_2),
\end{equation}
where $\mathbf{E}_T(r) = \bm{\Theta}_r\times(h\cdot\bm{\Theta}^{(1)})$ and $\bm{\Theta}^{(1)}$ is a compact set.
\medskip
We impose more structures in order to derive the asymptotic distribution.
\begin{assumption} \label{Ass3}
\item
\begin{enumerate}
\item $\inf_{\bm{\vartheta}\in\bm{\Theta}_r, \mathbf{z}\in (\mathbb{R}^m)^\infty}\lambda_{\min}\left(\mathbf{H}(\mathbf{z};\bm{\vartheta})\mathbf{H}(\mathbf{z};\bm{\vartheta})^\top\right) \geq \underline{c}$ for some $\underline{c}>0$.
\item For any $\bm{\vartheta}\in\bm{\Theta}_r$, $\bm{\mu}(\widetilde{\mathbf{z}}_t(\tau);\bm{\theta}(\tau)) =\bm{\mu}(\widetilde{\mathbf{z}}_t(\tau);\bm{\vartheta}) $ and $\mathbf{H}(\widetilde{\mathbf{z}}_t(\tau);\bm{\theta}(\tau)) =\mathbf{H}(\widetilde{\mathbf{z}}_t(\tau);\bm{\vartheta})$ a.s. imply $\bm{\vartheta} = \bm{\theta}(\tau)$ for some $t$, where $\widetilde{\mathbf{z}}_t(\tau)= \left(\widetilde{\mathbf{x}}_{t}(\tau),\widetilde{\mathbf{x}}_{t-1}(\tau),\ldots\right)$.
\end{enumerate}
\end{assumption}
\begin{assumption} \label{Ass4}
\item
\begin{enumerate}
\item $\bm{\mu}(\cdot;\bm\vartheta)$ and $\mathbf{H}(\cdot;\bm\vartheta)$ are twice continuously differentiable with respect to $\bm{\vartheta}$.
\item There exists a nonnegative sequence $\{\chi_j\}_{j=1}^{\infty}$ with $\chi_j = O(j^{-(2+s)})$ and some $s > 0$ such that for any $\mathbf{z},\mathbf{z}^\prime\in (\mathbb{R}^m)^\infty$ and any $\bm{\vartheta},\bm{\vartheta}^\prime\in\bm{\Theta}_r$:
\begin{eqnarray*}
&&|\gradient_{\bm{\vartheta}}^k\bm{\mu}(\mathbf{z};\bm{\vartheta})-\gradient_{\bm{\vartheta}}^k\bm{\mu}(\mathbf{z};\bm{\vartheta}^\prime)|+|\gradient_{\bm{\vartheta}}^k\mathbf{H}(\mathbf{z};\bm{\vartheta})-\gradient_{\bm{\vartheta}}^k\mathbf{H}(\mathbf{z};\bm{\vartheta}^\prime)| \leq |\bm{\vartheta}-\bm{\vartheta}^\prime| \sum_{j=1}^{\infty}\chi_j|\mathbf{z}_j|,\\
&&|\gradient_{\bm{\vartheta}}^k\bm{\mu}(\mathbf{z};\bm{\vartheta})-\gradient_{\bm{\vartheta}}^k\bm{\mu}(\mathbf{z}^\prime;\bm{\vartheta})|+|\gradient_{\bm{\vartheta}}^k\mathbf{H}(\mathbf{z};\bm{\vartheta})-\gradient_{\bm{\vartheta}}^k\mathbf{H}(\mathbf{z}^\prime;\bm{\vartheta})| \leq \sum_{j=1}^{\infty}\chi_j|\mathbf{z}_j-\mathbf{z}_j^\prime|,
\end{eqnarray*}
where $\gradient_{\bm{\vartheta}}= \left(\frac{\partial}{\partial \vartheta_1},\ldots,\frac{\partial}{\partial \vartheta_d}\right)^\top$, and $k=1,2$.
\end{enumerate}
\end{assumption}
\begin{assumption}\label{Ass5}
Let $K(\cdot)$ be a symmetric and positive kernel function defined on $[-1,1]$ with $\int_{-1}^{1}K(u)\mathrm{d}u = 1$. Moreover, $K(\cdot)$ is Lipschitz continuous on $[-1,1]$. As $(T,h) \to (\infty, 0)$, $Th\to \infty$.
\end{assumption}
Assumption \ref{Ass3}.1 ensures the positive definiteness of the covariance matrix of the likelihood function, and is widely adopted when studying the multivariate time series (e.g., page 2736 of \citealp{bardet2009asymptotic}). In fact, the validity of this assumption is easy to justify in view of \eqref{Eq4.32} and \eqref{Eq4.52} for Example 1 and Example 2 below. Assumption \ref{Ass3}.2 imposes an standard identification condition in the literature of M-estimation (e.g., Proposition 3.4 of \citealp{jeantheau1998strong}). It is noteworthy that the current form of Assumption \ref{Ass3} accommodates the flexibility of the model \eqref{Eq2.1}, which is in fact unnecessary if we have a detailed model in practice. See Section \ref{Sec2.4} for example.
Assumption \ref{Ass4} imposes the Lipschitz-type conditions on the first and second order derivatives of $\bm{\mu}(\cdot)$ and $\mathbf{H}(\cdot)$ to ensure the smoothness of their functional components.
Assumption \ref{Ass5} is a set of regular conditions on the kernel function and the bandwidth.
\medskip
With these conditions in hand, we summarize the first theorem of this paper below.
\begin{theorem}\label{Thm3.1}
Suppose Assumptions \ref{Ass1}--\ref{Ass5} hold with $r\geq 6$.
(1). If $Th^7 \to 0$, then for any $\tau \in (0,1)$
\begin{eqnarray*}
\sqrt{Th}\left(\widehat{\bm{\theta}}(\tau) - \bm{\theta}(\tau) - \frac{1}{2}h^2\widetilde{c}_2\bm{\theta}^{(2)}(\tau)\right) \to_D N\left(\mathbf{0}, \widetilde{v}_0 \bm{\Sigma}_{\bm{\theta}}(\tau) \right),
\end{eqnarray*}
where $\bm{\Sigma}_{\bm{\theta}}(\tau) = \bm{\Sigma}^{-1}(\tau)\bm{\Omega}(\tau)\bm{\Sigma}^{-1}(\tau)$, $\bm{\Sigma}(\tau)=E\left(\gradient_{\bm{\vartheta}}^2\mathcal{l}(\widetilde{\mathbf{x}}_1(\tau),\widetilde{\mathbf{z}}_{0}(\tau);\bm{\theta}(\tau)) \right)$ and
\begin{eqnarray*}
\bm{\Omega}(\tau)=E\left(\gradient_{\bm{\vartheta}}\mathcal{l}(\widetilde{\mathbf{x}}_1(\tau),\widetilde{\mathbf{z}}_0(\tau);\bm{\theta}(\tau))\cdot \gradient_{\bm{\vartheta}}\mathcal{l}(\widetilde{\mathbf{x}}_1(\tau),\widetilde{\mathbf{z}}_0(\tau);\bm{\theta}(\tau))^\top\right).
\end{eqnarray*}
(2). In addition, if $\bm{\varepsilon}_t$ is normally distributed, we have $\bm{\Omega}(\tau)=-\bm{\Sigma}(\tau)$ and thus $\bm{\Sigma}_{\bm{\theta}}(\tau) = \bm{\Omega}^{-1}(\tau)$.
\end{theorem}
After deriving the asymptotic distribution, we will establish both the point-wise inference and the simultaneous inference in the following.
\subsection{Inference}\label{Sec2.3}
In this section, we first discuss how to conduct point-wise inference, and then move on to derive the asymptotic results associated with the simultaneous inference. Specifically, for some preassigned significance level $\alpha \in (0,1)$, we shall construct a $100(1-\alpha)\%$ asymptotic simultaneous confidence band (SCB) $\{ \Upsilon(\tau), 0\leq \tau\leq 1 \}$ for $\bm{\theta}(\cdot)$ in the sense that
$$
\lim_{T\to\infty} \Pr\left(\bm{\theta}(\tau) \in \Upsilon(\tau), 0\leq \tau\leq 1\right) =1-\alpha.
$$
Notably, the simultaneous inference nests the traditional constancy test as a special case. It does not only allow one to examine whether a time-varying model should be preferred to its parametric counterpart, but also allows one to test any particular functional form of interest. For example, if a horizontal line can be embedded in the SCB $\{ \Upsilon(\tau)\}$, then we accept the hypothesis that some elements of $\bm{\theta}(\tau)$ are constant.
\medskip
\noindent \textbf{Point-wise Inference:} First, we construct a bias-corrected estimator in order to remove the asymptotic bias of Theorem \ref{Thm3.1}. Specifically, we let
\begin{equation}\label{Eq3.1}
\widetilde{\bm{\theta}} (\tau) = 2 \widehat{\bm{\theta}}_{h/\sqrt{2}}(\tau)-\widehat{\bm{\theta}} (\tau),
\end{equation}
where $\widehat{\bm{\theta}}_{h/\sqrt{2}}(\tau)$ is defined in the same way as $\widehat{\bm{\theta}} (\tau)$ but using the bandwidth $h/\sqrt{2}$.
After tedious development (Lemma B.7 of Appendix B), we have uniformly over $\tau \in[h,1-h]$
\begin{eqnarray*}
\widetilde{\bm{\theta}}(\tau) - \bm{\theta}(\tau) &=& -\bm{\Sigma}^{-1}(\tau)\frac{1}{Th}\sum_{t=1}^{T}\widetilde{K}((\tau_t-\tau)/h) \gradient_{\bm{\vartheta}}\mathcal{l}(\widetilde{\bm{x}}_t(\tau_t),\widetilde{\bm{z}}_{t-1}(\tau_t);\bm{\theta}(\tau_t))\\
&& +O_P((Th)^{-1/2}h^{3/2}(\log T)^{1/2}) + o(h^3),
\end{eqnarray*}
where $\widetilde{K}(x)=2\sqrt{2}K(\sqrt{2}x)-K(x)$ that is essentially a fourth-order kernel. It then infers that under the conditions of Theorem \ref{Thm3.1},
\begin{eqnarray}
\sqrt{Th}(\widetilde{\bm{\theta}} (\tau) - \bm{\theta}(\tau) ) \to_D N\left(\mathbf{0}, {v}_0 \bm{\Sigma}_{\bm{\theta}}(\tau) \right),
\end{eqnarray}
where $v_0 = \int_{-1}^{1}\widetilde{K}^2(u)\mathrm{d}u$.
It is noteworthy that the construction of \eqref{Eq3.1} is different from directly using the fourth-order kernel in the regression. In terms of bandwidth selection, the traditional methods (e.g., cross-validation) still remain valid for \eqref{Eq3.1} (\citealp{richter2019cross}). However, if one directly employs the fourth-order kernel in the regression, it remains unclear how to select the optimal bandwidth in practice.
\medskip
Now we discuss how to estimate $\bm{\Sigma}_{\bm{\theta}}(\tau) $ which is constructed by $\bm{\Sigma}(\tau)$ and $\bm{\Omega}(\tau)$. Intuitively, we consider the following estimator
\begin{equation}\label{Eq3.4}
\widehat{\bm{\Sigma}}_{\bm{\theta}}(\tau) = \widehat{\bm{\Sigma}}^{-1}(\tau) \widehat{\bm{\Omega}}(\tau) \widehat{\bm{\Sigma}}^{-1}(\tau),
\end{equation}
where
\begin{eqnarray*}
\widehat{\bm{\Sigma}}(\tau) &=&A_T(\tau)^{-1} \sum_{t=1}^{T}\gradient_{\bm{\vartheta}}^2\mathcal{l}(\mathbf{x}_t,\mathbf{z}_{t-1}^c;\widehat{\bm{\theta}}(\tau))K_h(\tau_t-\tau),\\
\widehat{\bm{\Omega}}(\tau) &=& A_T(\tau)^{-1} \sum_{t=1}^{T}\gradient_{\bm{\vartheta}}\mathcal{l}(\mathbf{x}_t,\mathbf{z}_{t-1}^c;\widehat{\bm{\theta}}(\tau)) \cdot \gradient_{\bm{\vartheta}}\mathcal{l}(\mathbf{x}_t,\mathbf{z}_{t-1}^c;\widehat{\bm{\theta}}(\tau))^\top K_h(\tau_t-\tau),\\
A_T(\tau)&=&\sum_{t=1}^{T}K_h(\tau_t-\tau).
\end{eqnarray*}
Note that we consider a local constant estimator in \eqref{Eq3.4} rather than a local linear one, that is to avoid an implementation issue for finite sample studies (i.e., nonpositive definite covariance may occur when the local linear approach is employed). Such a numerical problem has been well explained and investigated in the literature. See \cite{chen2015local} for example.
The following corollary summarizes the asymptotic property of \eqref{Eq3.4}.
\begin{corollary}\label{proposition3.1}
Under the conditions of Theorem \ref{Thm3.1}.1, suppose further that
\begin{eqnarray*}
\sup_{\tau \in [0,1]} [\alpha_j(\bm{\theta}(\tau))+\beta_j(\bm{\theta}(\tau))] = O(j^{-(5/2+s)})
\end{eqnarray*}
for some $s > 0$. In addition, let $h (\log T)^2\to 0$ and $T^{1-6/r}h\to \infty$. Then
\begin{eqnarray*}
\sup_{\tau \in [0,1]}|\widehat{\bm{\Sigma}}_{\bm{\theta}}(\tau) -\bm{\Sigma}_{\bm{\theta}}(\tau)|=o_P(1).
\end{eqnarray*}
\end{corollary}
\medskip
\noindent \textbf{Simultaneous Inference:} We now consider the simultaneous inference. To allow for flexibility, we first introduce a selection matrix $\mathbf{C}$ with full row rank, which selects the parameters of interest as follows:
\begin{eqnarray}
\bm{\theta}_{\mathbf{C}}(\tau):=\mathbf{C}\bm{\theta}(\tau).
\end{eqnarray}
Accordingly, the estimator and the corresponding asymptotic covariance matrix become
\begin{eqnarray}
\widehat{\bm{\theta}}_{\mathbf{C}}(\tau):=\mathbf{C}\widehat{\bm{\theta}} (\tau) \quad \text{and}\quad \bm{\Sigma}_{\mathbf{C}}(\tau) = \mathbf{C}\bm{\Sigma}_{\bm{\theta}}(\tau)\mathbf{C}^\top.
\end{eqnarray}
\begin{theorem}\label{Thm3.2}
Under the conditions of Theorem \ref{Thm3.1}.1, suppose further that
\begin{eqnarray*}
\sup_{\tau \in [0,1]} [\alpha_j(\bm{\theta}(\tau)) + \beta_j(\bm{\theta}(\tau))]= O(j^{-(3+s)})
\end{eqnarray*}
for some $s > 0$. In addition, let $(\log T)^4/(T^{\nu}h)\to 0$ with $\nu = \frac{1}{2}-\frac{r-6}{4rs/3+2r-4}$ and $Th^7\log T \to 0$. Then
\begin{eqnarray*}
\lim_{T\to\infty}\mathrm{Pr}&&\left(\sqrt{\frac{Th}{\widetilde{v}_0}}\sup_{\tau\in[h,1-h]}\left|\bm{\Sigma}_{\mathbf{C}}^{-1/2}(\tau)\left\{\widehat{\bm{\theta}}_{\mathbf{C}}(\tau) - \bm{\theta}_{\mathbf{C}}(\tau) - \frac{1}{2}h^2\widetilde{c}_2\bm{\theta}_{\mathbf{C}}^{(2)}(\tau) \right\} \right| \right.\\
&& \left.- B(1/h)\leq\frac{u}{\sqrt{2\log(1/h)}} \right) = \exp(-2\exp(-u)),
\end{eqnarray*}
where
\begin{eqnarray*}
B(1/h) &=& \sqrt{2 \log(1/h)} + \frac{\log (C_K) + (k/2-1/2)\log(\log(1/h))-\log(2)}{\sqrt{2 \log(1/h)}},\\
C_K &=& \frac{\{\int_{-1}^{1}|K^{(1)}(u)|^2\mathrm{d}u/\widetilde{v}_0\pi\}^{1/2}}{\Gamma(k/2)},
\end{eqnarray*}
and $\Gamma(\cdot)$ is the Gamma function.
\end{theorem}
In Theorem \ref{Thm3.2}, $\nu$ is slightly smaller than $1/2$ as we only require $r$ to be slightly larger than $6$. Hence, the usual optimal bandwidth $h_{opt} = O(T^{-1/5})$ satisfies the conditions $(\log T)^4/(T^{\nu}h)\to 0$ and $Th^7\log T \to 0$.
As shown in Theorem \ref{Thm3.2}, the convergence rate of the simultaneous confidence intervals for $\bm{\theta}_{\mathbf{C}}(\cdot)$ is of logarithmic rate and is therefore slow. In order to improve the rate, we consider a bootstrap method which shows a much better finite sample performance. We summarize the result in the following corollary.
\begin{corollary}\label{proposition3.4}
Under the conditions of Theorem \ref{Thm3.2}. Suppose that $h = O(T^{-\kappa})$ with $1/7<\kappa<\nu$. Then, on a richer probability space, there exists i.i.d. $k$-dimensional standard normal variables $\mathbf{v}_1,\ldots,\mathbf{v}_T$ such that
\begin{equation*}
\sup_{\tau\in[0,1]}|\widehat{\bm{\theta}}_{\mathbf{C}}(\tau)-\bm{\theta}_{\mathbf{C}}(\tau) -\frac{1}{2}h^2b_h(\tau)\bm{\theta}_{\mathbf{C}}^{(2)}(\tau) -\bm{\Sigma}_{\mathbf{C}}^{1/2}(\tau)\mathbf{V}_h^*(\tau)| =O_P\left(\frac{T^{-\alpha}}{\sqrt{Th\log T}} \right),
\end{equation*}
where $\alpha = \min\{(\nu-\kappa)/2,(7\kappa-1)/2,\kappa/2\}$, $\widetilde{c}_{k,h}(\tau) = \int_{-\tau/h}^{(1-\tau)/h} u^k K(u) \mathrm{d}u $, $\mathbf{V}_h^*(\tau) = T^{-1}\sum_{t=1}^{T}\mathbf{v}_t\omega_{t,h}(\tau)$,
$$
b_h(\tau) = \frac{\widetilde{c}_{2,h}^2(\tau) - \widetilde{c}_{1,h}(\tau)\widetilde{c}_{3,h}(\tau)}{ \widetilde{c}_{0,h}(\tau)\widetilde{c}_{2,h}(\tau)-\widetilde{c}_{1,h}^2(\tau)}\quad \text{and}\quad \omega_{t,h}(\tau) = K_h(\tau_t-\tau)\frac{ \widetilde{c}_{2,h}(\tau)- \frac{\tau_t-\tau}{h}\widetilde{c}_{1,h}(\tau)}{\widetilde{c}_{0,h}(\tau)\widetilde{c}_{2,h}(\tau)-\widetilde{c}_{1,h}^2(\tau)}.
$$
\end{corollary}
By Corollary \ref{proposition3.4}, we propose the following numerical procedure to construct the SCB of $\bm{\theta}_{\mathbf{C}}(\tau)$:
\begin{itemize}
\item[Step 1] Use the sample $\{\mathbf{x}_t\}_{t=1}^T$ to estimate $\widehat{\bm{\theta}}_{\mathbf{C}}(\tau)$ by \eqref{Eq2.4}, and compute $\widetilde{\bm{\theta}}_{\mathbf{C}}(\tau)$ based on \eqref{Eq3.1}.
\item[Step 2] Generate i.i.d. $k$-dimensional standard normal variables $\{\mathbf{v}_t^*\}$ and calculate the quantity $\sup_{\tau \in [0,1]}|\mathbf{V}_{h}^*(\tau)|$, in which $\mathbf{V}_{h}^*(\tau)= T^{-1}\sum_{t=1}^{T}\mathbf{v}_t^*(2\omega_{t,h/\sqrt{2}}(\tau)-\omega_{t,h}(\tau))$.
\item[Step 3] Repeat Step 2 $R$ times to obtain the empirical $(1-\alpha)^{th}$ quantile $\widehat{q}_{1-\alpha}$ of $\sup_{\tau \in [0,1]}|\mathbf{V}_{h}^*(\tau)|$.
\item[Step 4] Calculate $\widehat{\bm{\Sigma}}_{\mathbf{C}}(\tau)$ using \eqref{Eq3.4}, and construct the SCB of $\bm{\theta}_{\mathbf{C}}(\tau)$ by $\widetilde{\bm{\theta}}_{\mathbf{C}}(\tau) + \widehat{\bm{\Sigma}}_{\mathbf{C}}^{1/2}(\tau) \widehat{q}_{1-\alpha} \mathbb{B}_k$, where $\mathbb{B}_k = \{\mathbf{u}\in \mathbb{R}^k:|\mathbf{u}|\leq 1\}$ is the unit ball, and $k$ is the rank of $\mathbf{C}$.
\end{itemize}
\subsection{Examples} \label{Sec2.4}
Below, we demonstrate the usefulness of the aforementioned results by considering Example 1 and Example 2 of Section \ref{Sec1}.
\smallskip
\noindent \textbf{Example 1 (Cont.)} --- For $\forall\tau\in[0,1]$, simple algebra shows that the approximated stationary process is defined by
\begin{eqnarray}\label{Eq4.3}
\widetilde{\bm{x}}_t(\tau) = \bm{\mu}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right) +\mathbf{H}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right)\bm{\varepsilon}_t,
\end{eqnarray}
where $\bm{\theta}(\tau)$ has been defined in \eqref{Eq4.1.1}, and
\begin{eqnarray}\label{Eq4.32}
\bm{\mu}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right) &=& \mathbf{B}_{\tau}^{-1}(1)\mathbf{a}(\tau)+\sum_{j=1}^{\infty}\bm{\Gamma}_j(\tau)\widetilde{\bm{x}}_{t-j}(\tau),\nonumber \\
\mathbf{H}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right) &=& \bm{\omega}(\tau).
\end{eqnarray}
Additionally, in \eqref{Eq4.32}, $\bm{\Gamma}_j(\tau)$ is yielded as follows:
\begin{eqnarray}
\mathbf{I}_m - \sum_{j=1}^{\infty}\bm{\Gamma}_j(\tau)L^j = \mathbf{B}_{\tau}^{-1}(L)\mathbf{A}_{\tau}(L),
\end{eqnarray}
where $\mathbf{A}_{\tau}(L):=\mathbf{I}_m - \mathbf{A}_1(\tau)L-\cdots-\mathbf{A}_p(\tau)L^p$
and $\mathbf{B}_{\tau}(L):=\mathbf{I}_m + \mathbf{B}_1(\tau)L+\cdots+\mathbf{B}_q(\tau)L^q$.
Then we are able to present the following proposition.
\begin{proposition}\label{proposition4.1}
Let $\|\bm{\varepsilon}_t\|_r < \infty$ for some $r>4$. Suppose that there is a compact set
\begin{eqnarray*}
\bm{\Theta} =\{\bm{\vartheta} = \mathrm{vec}(\mathbf{a} ,\mathbf{A}_1 ,\ldots,\mathbf{A}_p ,\mathbf{B}_1 ,\ldots,\mathbf{B}_q ,\bm{\Omega} ) \mid \bm{\vartheta} \in \mathbb{R}^d\}
\end{eqnarray*}
such that (1). for $\forall \tau\in [0,1]$, $\bm{\theta}(\tau)$ lies in the interior of $\bm{\Theta}$, (2). $\mathrm{det}(\mathbf{A}(L)\mathbf{B}(L)) \neq 0$ for all $|L|\leq 1$, (3). $\bm{\Omega}>0$, where $\mathbf{A}(L):=\mathbf{I}_m - \mathbf{A}_1L-\cdots-\mathbf{A}_pL^p$ and $\mathbf{B}(L):=\mathbf{I}_m + \mathbf{B}_1L+\cdots+\mathbf{B}_qL^q$ are coprime and satisfy some necessary identification conditions. Then, the results of Theorems \ref{Thm3.1} and \ref{Thm3.2} hold for model \eqref{Eq4.1}.
\end{proposition}
We note that the detailed identification conditions required for VARMA processes (e.g., the final equations form or echelon form) can be found in \cite{lutkepohl2005new}. We no longer discuss them here in order not to derivative from our main goal.
\medskip
\noindent \textbf{Example 2 (Cont.)} --- We further let
\begin{eqnarray}\label{GARCH_rho}
\bm{\Omega}(\tau) = \left[\begin{matrix}
1 & \rho_{1,2}(\tau) & \cdots &\rho_{1,m}(\tau) \\
\rho_{1,2}(\tau) & 1 & \ddots &\vdots \\
\vdots & \ddots & \ddots &\rho_{m-1,m}(\tau) \\
\rho_{1,m}(\tau) & \rho_{m-1,m}(\tau) & \ddots &1 \\
\end{matrix} \right].
\end{eqnarray}
For $\forall\tau\in[0,1]$, the corresponding approximated stationary process is defined as
\begin{eqnarray}\label{Eq4.5}
\widetilde{\bm{x}}_t(\tau) = \mathbf{H}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right)\bm{\varepsilon}_t,
\end{eqnarray}
where
\begin{eqnarray}\label{Eq4.52}
&&\mathbf{H}\left(\widetilde{\mathbf{x}}_{t-1}(\tau),\widetilde{\mathbf{x}}_{t-2}(\tau),\ldots;\bm{\theta}(\tau)\right)\nonumber \\
& =& \mathrm{diag}^{1/2}\left( \mathbf{D}_{\tau}^{-1}(1)\mathbf{c}_0(\tau)+ \sum_{j=1}^{\infty} \mathbf{\Psi}_j(\tau) \left(\widetilde{\mathbf{x}}_{t-j}(\tau)\odot\widetilde{\mathbf{x}}_{t-j}(\tau)\right)\right).
\end{eqnarray}
Note that $ \mathbf{\Psi}_j(\tau)$ is generated as follows:
\begin{eqnarray}
\bm{\Psi}_{\tau}(L):=\mathbf{I}_m - \sum_{j=1}^{\infty}\bm{\Psi}_j(\tau)L^j=\mathbf{D}_{\tau}^{-1}(L)\mathbf{C}_{\tau}(L),
\end{eqnarray}
where $\mathbf{C}_{\tau}(L):=\mathbf{C}_1(\tau)L+\cdots+\mathbf{C}_p(\tau)L^p$ and $\mathbf{D}_{\tau}(L):=\mathbf{I}_m - \mathbf{D}_1(\tau)L-\cdots-\mathbf{D}_q(\tau)L^q$.
Consequently, we can present the following proposition.
\begin{proposition}\label{proposition4.2}
Suppose that there is a compact set
\begin{eqnarray*}
\bm{\Theta} =\{\bm{\vartheta} = \mathrm{vec}(\mathbf{c}_0,\mathbf{C}_1,\ldots,\mathbf{C}_p ,\mathbf{D}_1,\ldots,\mathbf{D}_q,\bm{\Omega}) \mid \bm{\vartheta} \in \mathbb{R}^d\}
\end{eqnarray*}
such that (1). for $\forall \tau\in [0,1]$, $\bm{\theta}(\tau)$ lies in the interior of $\bm{\Theta}$, (2). $\|\bm{\Omega}^{1/2}\bm{\varepsilon}_t\|_r^2 \sum_{j=1}^\infty|\mathbf{\Psi}_j|< 1$ for some $r > 6$, (3). all the roots of $|\mathbf{I}_m - \sum_{j=1}^{p} \mathbf{C}_j - \sum_{j=1}^{q} \mathbf{D}_j|$ are outside the unit circle with $\mathbf{C}_j$'s and $\mathbf{D}_j$'s being squared matrices of nonnegative elements, (4). $\mathbf{c}_0$ is a vector of positive elements, (5). $\mathbf{C}(L)$ and $\mathbf{D}(L)$ are coprime and the formulation of the GARCH part is minimal, where $\mathbf{C}(L):=\mathbf{C}_1 L+\cdots+\mathbf{C}_p L^p$ and $\mathbf{D} (L):=\mathbf{I}_m - \mathbf{D}_1 L-\cdots-\mathbf{D}_q L^q$. Then the results Theorems \ref{Thm3.1} and \ref{Thm3.2} hold for model \eqref{Eq4.4}.
\end{proposition}
For the identification conditions of the GARCH process, we refer readers to Proposition 3.4 of \cite{jeantheau1998strong}, who proves that assuming the minimal representation is enough for ensuring Assumption \ref{Ass3} holds.
In the following section, we conduct numerical studies using both simulated and real data to evaluate the finite-sample performance of the proposed estimation and inferential methods.
\section{Numerical Studies}\label{Sec3}
In this section, we first present the details of the numerical implementations in Section \ref{Sec3.1}, and then conduct extensive simulations in Section \ref{Sec3.2}. Section \ref{Sec3.3} presents a real data example on the conditional correlations between the Chinese and U.S. stock markets.
\subsection{Numerical Implementation}\label{Sec3.1}
Throughout the numerical studies, the Epanechnikov kernel $K(u) = 0.75(1-u^2)I(|u|\leq1)$ is adopted. Following \cite{zhou2010simultaneous}, we use $\widetilde{h} = 2\widehat{h}$ for the biased corrected estimator, where $\widehat{h}$ is the bandwidth selected by the cross-validation method of \cite{richter2019cross}.
Specifically, define the leave-one-out local linear QMLE
\begin{equation}\label{Eq5.1}
(\widehat{\bm{\theta}}_{h,-t}(\tau),h\widehat{\bm{\theta}}_{h,-t}^{(1)}(\tau)) = \operatorname*{\arg\!\max}_{(\bm{\eta}_1,\bm{\eta}_2)\in\mathbf{E}_T(r)}\mathcal{L}_{T,-t}^c(\tau,\bm{\eta}_1,\bm{\eta}_2),
\end{equation}
where
$$
\mathcal{L}_{T,-t}^c(\tau,\bm{\eta}_1,\bm{\eta}_2)=\frac{1}{T}\sum_{s=1,\neq t}^{T}\mathcal{l}(\mathbf{x}_s,\mathbf{z}_{s-1}^c;\bm{\eta}_1+\bm{\eta}_2\cdot(\tau_s-\tau)/h)K_h(\tau_s-\tau).
$$
Then, the bandwidth is chosen by
\begin{equation}\label{Eq5.2}
\widehat{h} = \operatorname*{\arg\!\max}_{h}T^{-1}\sum_{t=1}^{T}\mathcal{l}(\mathbf{x}_t,\mathbf{z}_{t-1}^c;\widehat{\bm{\theta}}_{h,-t}(\tau_t)).
\end{equation}
As shown in \cite{richter2019cross}, this cross validation method works well as long as $\gradient\mathcal{l}$ is uncorrelated, which implies that this desirable property should hold in our case.
Notably, when considering some specific models, the implementation may be further simplified. We provide more discussions along this line in Appendix B.4.
\subsection{Simulation Results}\label{Sec3.2}
In the simulation studies, we examine the empirical coverage probabilities of simultaneous confidence intervals for nominal levels $\alpha =90\%,\ 95\%$. We consider the time-varying VARMA($2,1$) and multivariate GARCH($1,1$) model as follows:
\begin{enumerate}
\item $\text{DGP 1}: \mathbf{x}_t = a_1(\tau_t)\mathbf{x}_{t-1}+a_2(\tau_t)\mathbf{x}_{t-2} + \bm{\eta}_t + \mathbf{B}_1(\tau_t)\bm{\eta}_{t-1},\quad \bm{\eta}_t = \bm{\omega}(\tau_t) \bm{\varepsilon}_t$, where $\{\bm{\varepsilon}_t\}$ are i.i.d. draws from $N(\mathbf{0}_{2\times 1},\mathbf{I}_2)$, $a_1(\tau) = 0.6\exp(\tau-1)$, $a_2(\tau) = -0.3\exp(\tau-1)$,
\begin{eqnarray*}
\mathbf{B}_1(\tau)&=&\left[\begin{matrix}
0.5\exp{\tau-0.5} & -0.8(\tau-0.5)^2 \\
-0.8(\tau-0.5)^2 & 0.5+0.3\sin(\pi \tau)
\end{matrix} \right], \nonumber \\
\bm{\omega}(\tau)&=&\left[\begin{matrix}
1.5+0.2\exp{0.5-\tau}& 0 \\
0.2\exp{0.5-\tau} & 1.5+0.5(\tau-0.5)^2
\end{matrix}\right].
\end{eqnarray*}
Here we use final equations form to ensure the uniqueness of the VARMA representation.
\item $\text{DGP 2}: \mathbf{x}_t = \mathrm{diag}(h_{1,t}^{1/2},\ldots,h_{m,t}^{1/2}) \bm{\eta}_t$, where $\bm{\eta}_t = \bm{\Omega}^{1/2}(\tau_t) \bm{\varepsilon}_t$, $\mathbf{h}_{t} = \mathbf{c}_0(\tau_t) + \mathbf{C}_1(\tau_t) \left(\mathbf{x}_{t-1}\odot\mathbf{x}_{t-1}\right) + \mathbf{D}_1(\tau_t)\mathbf{h}_{t-1}$, $\{\bm{\varepsilon}_t\}$ are i.i.d. draws from $N(\mathbf{0}_{2\times 1},\mathbf{I}_2)$, $\mathbf{c}_0(\tau) = [2\exp{0.5\tau-0.5}, 3+0.2\cos(\tau)]^\top$,
\begin{eqnarray*}
\mathbf{C}_1(\tau)&=&\left[\begin{matrix}
0.4 + 0.05\cos(\tau) & 0.05(\tau-0.5)^2 \\
0.05(\tau-0.5)^2 & 0.4+0.05\sin(\tau)
\end{matrix} \right], \nonumber \\
\mathbf{D}_1(\tau)&=&\left[\begin{matrix}
0.4 - 0.1\cos(\tau) & 0 \\
0 & 0.3 - 0.1\sin(\tau)
\end{matrix} \right], \nonumber \\
\bm{\Omega}(\tau)&=&\left[\begin{matrix}
1 & 0.3\sin(\tau) \\
0.3\sin(\tau) & 1
\end{matrix}\right].
\end{eqnarray*}
\end{enumerate}
Let the sample size be $T \in\{500,1000\}$ ($T \in\{1000,2000,4000\}$) for the VARMA model (the GARCH model). We conduct $1000$ replications for each choice of $T$. Several different bandwidths close to $\widetilde{h}$ are reported to check the sensitivity of bandwidth selection.
We present the empirical coverage probabilities associated with the SCB
in Tables \ref{table_sim1}--\ref{table_sim2}. For the vector- or matrix-valued unknown coefficients, we take an average across the elements. A few facts emerge from the tables. First, the finite sample coverage probabilities are smaller than their nominal level when $T = 500$ ($T = 1000,2000$) for the VARMA model (the GARCH model), but are fairly close to their nominal level as $T = 1000$ ($T=4000$) for the VARMA model (the GARCH model). Second, the behaviour of the estimated simultaneous confidence intervals is not sensitive to the choices of bandwidths. Third, the GARCH model requires more data to reach a reasonable finite sample performance.
\begin{table}[h]
\caption{Empirical Coverage Probabilities of the SCB for DGP 1}\label{table_sim1}
\begin{center}
\begin{tabular}{c c c cccc c cccc}
\hline
& & &\multicolumn{4}{c}{$90\%$}& &\multicolumn{4}{c}{$95\%$}\\
\cline{4-7} \cline{9-12}
&\text{$\widetilde{h}$}& &$\alpha_1(\cdot)$ &$\alpha_2(\cdot)$ &$\mathbf{B}_1(\cdot)$ & $\bm{\Omega}(\cdot)$ & &$\alpha_1(\cdot)$ &$\alpha_2(\cdot)$ &$\mathbf{B}_1(\cdot)$ & $\bm{\Omega}(\cdot)$\\
\hline
\multirow{4}{*}{\shortstack{$T=500$}}
&$0.35$ & & 0.845 & 0.877 & 0.821&0.847 & & 0.905 & 0.915 & 0.889 &0.905\\
&$0.4$ & & 0.865& 0.875 & 0.847 &0.876 & & 0.912 & 0.930 & 0.897& 0.909\\
&$0.45$ & & 0.862 & 0.895 & 0.847 &0.878 & &0.915 &0.930 & 0.898 &0.919\\
&$0.5$ & & 0.875 &0.895 & 0.847 &0.876 & &0.905 &0.945 & 0.901 &0.920\\
\hline
\multirow{4}{*}{\shortstack{$T=1000$}}
&$0.3$ & & 0.895 & 0.925 & 0.887& 0.884& & 0.960 & 0.960 & 0.947& 0.947\\
&$0.35$ & & 0.910 & 0.927 & 0.886&0.890& & 0.940 & 0.967 & 0.940&0.930 \\
&$0.4$ & & 0.917 & 0.939 & 0.901&0.899 & & 0.947 & 0.959 & 0.948 &0.939\\
&$0.45$ & &0.937 & 0.932 &0.908 &0.895 & &0.957 &0.957 &0.947 &0.937 \\
\hline
\end{tabular}
\end{center}
\end{table}
\begin{table}[h]
\caption{Empirical Coverage Probabilities of the SCB for DGP 2}\label{table_sim2}
\begin{center}
\begin{tabular}{c c c cccc c cccc}
\hline
& & &\multicolumn{4}{c}{$90\%$}& &\multicolumn{4}{c}{$95\%$}\\
\cline{4-7} \cline{9-12}
&\text{$\widetilde{h}$}& &$\mathbf{c}_0(\cdot)$ &$\mathbf{C}_1(\cdot)$ &$\mathbf{D}_1(\cdot)$ & $\bm{\Omega}(\cdot)$ & &$\mathbf{c}_0(\cdot)$ &$\mathbf{C}_1(\cdot)$ &$\mathbf{D}_1(\cdot)$ & $\bm{\Omega}(\cdot)$\\
\hline
\multirow{4}{*}{\shortstack{$T=1000$}}
&$0.55$ & & 0.802& 0.810 & 0.784 &0.889 & & 0.869 & 0.876 & 0.838 & 0.945\\
&$0.60$ & & 0.824 & 0.820 & 0.791 &0.879 & &0.882 &0.866 & 0.843 &0.945\\
&$0.65$ & & 0.832 &0.820 & 0.796 &0.874 & &0.889 &0.872 & 0.859 &0.945\\
&$0.70$ & & 0.820 &0.823 & 0.792 &0.879 & &0.892 &0.881 & 0.871 &0.950\\
\hline
\multirow{4}{*}{\shortstack{$T=2000$}}
&$0.50$ & & 0.827 & 0.835 & 0.841 &0.889 & & 0.897 & 0.881 & 0.901&0.950\\
&$0.55$ & & 0.829 & 0.825 & 0.843 &0.884 & & 0.892 & 0.881 & 0.903 & 0.940 \\
&$0.60$ & & 0.849 & 0.833 & 0.871 &0.900 & &0.900 &0.888 & 0.910 &0.950\\
&$0.65$ & & 0.852 &0.835 & 0.873 &0.910 & &0.907 &0.889 & 0.910 &0.950\\
\hline
\multirow{4}{*}{\shortstack{$T=4000$}}
&$0.35$ & & 0.879 & 0.879 & 0.882 &0.869& & 0.929 & 0.932 & 0.943 &0.920 \\
&$0.4$ & & 0.899 & 0.879 & 0.882 &0.859 & & 0.950 & 0.944 & 0.943 &0.919\\
&$0.45$ & &0.904 & 0.899 &0.879 &0.838 & &0.950 &0.947 &0.946 &0.950 \\
&$0.50$ & &0.867 & 0.857 &0.884 &0.898 & &0.929 &0.944 &0.946 &0.960 \\
\hline
\end{tabular}
\end{center}
\end{table}
\subsection{A Real Data Example}\label{Sec3.3}
In this subsection, we investigate the time-varying conditional correlations between the Chinese and U.S. stock markets using the time-varying multivariate GARCH model. Recently, there is a growing literature to study the relationship of the two stock markets (e.g., \citealp{zhang2014has,pan2022modeling}), as the Chinese stock market has become the world's second largest stock market after 2009. Understanding the interactions among different financial markets is important for investors and policymakers \cite[]{diebold2009measuring,bensaida2019good}. For example, high equity market interdependence implies poor diversification benefits from portfolios, but highlights the possibility of better hedging benefits.
Previous research documents a strong positive link between the degree of globalization and equity market interdependence \cite[]{baele2005volatility}. Along this line of research, one important question is that whether the interdependence between the Chinese and U.S. stock markets has increased over time due to globalization so that estimates from historical data are unreliable for modern policy analysis, asset pricing and risk management. The existing results present many discrepancies, which may be due to the fact that the relationship evolves with time. Apparently, the results also indicate that one should use time-varying GARCH model to accommodate potential nonstationarity inherited in these financial variables. In addition, as pointed out by \cite{caporin2013ten}, dynamic conditional correlation (DCC) GARCH model represents the dynamic conditional covariances of the standardized residuals, and hence does not yield dynamic
conditional correlations; DCC yields inconsistent two step estimators; DCC has no asymptotic properties. In what follows, we address these issues using the newly proposed approach. The estimation is conducted in exactly the same way as in Section \ref{Sec3.1}, so we no longer repeat the details.
We calculate the Chinese and U.S. stock returns based on weekly Shanghai Stock Exchange (SSE) Composite Index and S\&P 500 Index as they are the most comprehensive and diversified stock indices. The sample employed in this study spanning from January 2000 to February 2022 provides $1119$ observations\footnote{The data are collected from Yahoo Finance at \url{https://finance.yahoo.com/}.}. Figure \ref{Fg1} plots the two weekly returns as well as sample autocorrelation functions of squared data, which shows the typical ``volatility clustering'' phenomenon.
\begin{figure}[h]
\centering
{\includegraphics[width=16cm]{timeseriesplot.pdf}}
\caption{\small S\&P 500 and SSE Index returns as well as sample autocorrelation functions of squared data}\label{Fg1}
\end{figure}
We next fit the data to a time-varying multivariate GARCH(1,1) model and are particularly interested in the estimates of time-varying conditional correlations, i.e.,
$$
E\left(x_{1,t}x_{2,t}\mid \mathcal{F}_{t-1}\right)/\sqrt{E\left(x_{1,t}^2\mid \mathcal{F}_{t-1}\right)E\left(x_{2,t}^2\mid \mathcal{F}_{t-1}\right)} = \rho_{1,2}(\tau_t),
$$
where $\rho_{1,2}(\cdot)$ is defined in \eqref{GARCH_rho}. Figure \ref{Fg2} plots the estimates (black solid line) of time-varying conditional correlations between the two stock markets as well as 95\% simultaneous confidence intervals (red dashed line) and 95\% pointwise confidence intervals (black dashed line). Based on the simultaneous confidence intervals, apparently, the conditional correlations vary with respect to time. Moreover, as clearly presented in Figure \ref{Fg2}, the interdependence between the two stock markets is increasing over time. By examining the pointwise confidence intervals, we can conclude that the two stock markets are not significantly correlated before 2005, but the relationship has been greatly enhanced in recent years. These results have important implications for investment and risk management. For example, it implies that the Chinese and U.S. investors who use cross-country portfolio strategies to eliminate country specific risks may be benefit from hedging. However, all types of investors should be cautious since the relations between the Chinese and U.S. stock markets are time-varying.
\begin{figure}[h]
\centering
{\includegraphics[width=14cm]{TVCC.pdf}}
\caption{\small Time-varying conditional correlations between the Chinese and U.S. stock markets}\label{Fg2}
\end{figure}
\section{Conclusions}\label{Sec4}
In this paper, we consider a wide class of time-varying multivariate causal processes which nests many classic and new examples as special cases. We first prove the existence of a weakly dependent stationary approximation for the model \eqref{Eq2.1} which is the foundation to establish the corresponding asymptotic properties. Afterwards, we consider the QMLE estimation approach, and provide both point-wise and simultaneous inferences on the coefficient functions. In addition, we demonstrate the theoretical findings through both simulated and real data examples. In particular, we show the empirical relevance of our study using an application to evaluate the conditional correlations between the stock markets of China and U.S. We find that the interdependence between the two stock markets is increasing over time.
There are several directions for possible extensions. The first one is to consider quantile regression methods for such locally stationary multivariate causal processes. The second one is to propose a more powerful $L_2$ test based on the weighted integrated squared errors for testing whether some coefficients are time-invariant. We wish to leave such issues for future study.
\section{Acknowledgements}
The authors of this paper would like to thank George Athanasopoulos, David Frazier and Gael Martin for their constructive comments on earlier versions of this paper. Thanks also go to seminar participants for their insightful suggestions. Gao and Peng would also like to acknowledge the Australian Research Council Discovery Projects Program for its financial support under Grant Numbers: DP170104421 \& DP210100476.
{\footnotesize
\bibliography{Bibliography-MM-MC}
}
\bigskip
{\small
\begin{center}
{\large \textbf{Online Supplementary Appendices to \\``Time-Varying Multivariate Causal Processes"}}
\bigskip
{\sc Jiti Gao$^\ast$ and Bin Peng$^{\ast}$ and Wei Biao Wu$^{\dag}$ and Yayi Yan$^{\ast}$}
\medskip
$^\ast$Monash University and $^{\dag}$University of Chicago
\medskip
\today
\end{center}
The file includes Appendix A and Appendix B. We first present some technical tools in Appendix \ref{App.A1}, which will be repeatedly used in the development. We then provide the proofs of main results in Appendix \ref{App.A3}. We provide several preliminary lemmas in Appendix \ref{AppB.1} as well as some secondary lemmas in Appendix \ref{AppB.2}, and then present the proofs of preliminary lemmas in Appendix \ref{AppB.3}. Appendix \ref{AppB.4} discusses several computational issues of the local linear ML estimation.
In what follows, $M$ and $O(1)$ always stand for some bounded constants, and may be different at each appearance.
{\small