*

#  Estimating Conditional Value-at-Risk with Nonstationary Quantile Predictive
Regression Models††thanks: Article history: First draft 24 April 2022. Second
draft in revised form 29 August 2023.

Christis Katsouris
University of Helsinki

Working paper Dr. Christis Katsouris is currently a Postdoctoral Researcher at
the Faculty of Social Sciences, University of Helsinki, P.O. Box 17,
Arkadiankatu 7 Helsinki, FI-00014 Finland. E-mail:
<EMAIL_ADDRESS>
(August 31, 2024)

###### Abstract

This paper develops an asymptotic distribution theory for a two-stage
instrumentation estimation approach in quantile predictive regressions when
both generated covariates and persistent predictors are used. The generated
covariates are obtained from an auxiliary quantile regression model and our
main interest is the robust estimation and inference of the primary quantile
predictive regression in which this generated covariate is added to the set of
nonstationary regressors. We find that the proposed doubly IVX corrected
estimator is robust to the abstract degree of persistence regardless of the
presence of generated regressor obtained from the first stage procedure. The
asymptotic properties of the two-stage IVX estimator such as mixed Gaussianity
are established while the asymptotic covariance matrix is adjusted to account
for the first-step estimation error.

Keywords: risk measures; generated regressors; quantile regression;
persistence.

JEL Classification: C12, C22

##  Introduction

Estimating risk measures such as the conditional value at risk entails a
generated regressor problem. Moreover, econometric inference with generated
regressors is a commonly used approach when one needs to construct an
explanatory variable using structural econometric methods. In this article, we
present an approach for estimation and inference in nonstationary time series
regressions with generated regressors, which has not seen much attention in
the literature so far and thus requires to develop and establish our own
asymptotic theory results.

Inferential procedures based on generated regressors are studied in linear and
conditional quantile models as in Chen et al., (2021) (see, Pagan, (1984)111In
particular, Pagan, (1984) examines the statistical properties of regression
estimators involving generated regressors in the parametric context but not
corresponding to a time series environment with nonstationarity. , Doran and
Griffiths, (1983), Hoffman, (1987), Oxley and McAleer, (1993), Dufour and
Jasiak, (2001) and Chen et al., (2023)), although in the case of risk measures
such as $\mathsf{VaR}$ and $\mathsf{CoVaR}$ a conditional quantile
specification is required. The conventional approach for estimating these risk
measures is to employ stationary quantile regressions (as in Adrian and
Brunnermeier, (2016) and Härdle et al., (2016)). Several studies in the
literature consider the estimation of risk measures using stationary time
series environments either with parametric (e.g., see He et al., (2020)) or
semiparametric approaches (e.g., see, Wang and Zhao, (2016)). In this study,
we extend the econometric specifications for the joint estimation of the pair
$(\mathsf{VaR},\mathsf{CoVaR})$ as discussed in the study of Katsouris, 2023c
(see also Katsouris, (2021); Katsouris, 2023b ), to the case of nonstationary
quantile predictive regression models (see, Lee, (2016)), which implies using
a local-to-unity parametrization for modelling the unknown form of
persistence.

Specifically, Katsouris, 2023c  considers the statistical estimation for
covariance structures with tail estimates using nodewise quantile predictive
regression models, extending the framework of Katsouris, (2021) who discusses
identification and estimation issues of conditional quantile risk measures. On
the other hand, in both of these studies the estimation is based on stationary
time series models (see also, Patton et al., (2019)), thereby relying on OLS-
based estimators. Moreover, a unified framework that presents an estimation
approach and the related asymptotic theory under the presence of both a
generated regressor and persistent regressors is still an open problem. In
addition, a vast stream of literature develops estimation and inference
methods in predictive regression models robust to the unknown form of
persistence (see, Phillips, (2023)), with a specific approach being the IVX
filtration of Phillips and Magdalinos, (2009). Motivated by the aforementioned
issues, our research objective is to study the large-sample theory of a doubly
IVX corrected estimator in nonstationary quantile predictive regression models
for the purpose of estimating risk measures. To the best of our knowledge, our
approach that incorporates time series nonstationarity when estimating risk
measures in a low dimensional setting, is novel. Our contributions to the
literature include the study of the properties of doubly corrected estimators
in a two-stage estimation procedure in nonstationary time series regression
models.

As a result, when developing robust estimation and inference methods in
econometric environments with generated regressors, it is crucial to obtain
analytical and computational feasible expressions for the asymptotic variance
of these potentially semiparametric estimators. In particular, Rilstone,
(1996) investigates the aspect of generated covariate in a nonparameteric
stationary time series environment. Moreover, Hahn and Ridder, (2013) study
the asymptotic distribution of three-step estimators of a finite-dimensional
parameter vector where the second step consists of one or more nonparametric
regressions on a regressor that is estimated in the first step. Specifically,
the authors using the path-derivative method of Newey, (1994), derive the
contribution of the first-step estimator to the influence function is used to
investigate the dual role that the first-step estimator plays in the second-
step nonparametric regression. Moreover, Bhattacharya, (2020) consider the
aspects of estimation and inference in stationary quantile regression models
with generated covariates. Recently, Ghosh and Tan, (2022) develop a framework
for semiparametric estimation where a doubly robust estimating function for a
low-dimensional parameter is available, depending on two working models.
Therefore, our framework proposes a computationally tractable two-step
algorithm and an asymptotic theory analysis for the doubly IVX corrected
estimator (which is found to robustify inference to the unknown form of
persistence222The unknown form of persistence in time series data is modelled
via the the local-to-unity parametrization introduced by the seminal study of
Phillips, 1987a ; Phillips, 1987b  who established a unified framework for
unit root and nearly integrated processes in time series models.), when
estimating risk measures under the assumption of statistical elicitability
(see, Fissler and Ziegel, (2016)).

More specifically, our practical two stage estimation procedure permits to
obtain estimates for the forecasts of the pair $(\mathsf{VaR},\mathsf{CoVaR})$
as per the econometric specifications of Adrian and Brunnermeier, (2016) (see
also Härdle et al., (2016)), but with the additional assumption that the
stochastic processes that correspond to the regressors are generated from
nonstationary autoregressive equations via the LUR parametrization as in Lee,
(2016) and Fan and Lee, (2019). During the first-stage procedure using the
nonstationary quantile predictive regression model (QPR) and the local-to-
unity parametrized regressors, we obtain quantile dependent model estimates
that correspond to the one-period ahead forecasted $\mathsf{VaR}$ risk
measure. During the second stage estimation procedure, the generated regressor
that corresponds to the $\mathsf{VaR}$ is employed as an additional covariate
to a different QPR model which is employed to obtain the one-period ahead
forecasted $\mathsf{CoVaR}$ risk measure. Because of the well-known
endogeneity problem (see, Cai and Wang, (2014)) that appears due to the
correlation between the two innovation sequences of the predictive regression
and the autoregressive equation, parameter estimation is achieved using the
IVX filtration (see, Phillips and Magdalinos, (2009) and Kostakis et al.,
(2015)), which is robust to the unknown persistence properties of regressors
and is applied to both stages of the procedure that allows to identify and
jointly estimate the risk measures of $\mathsf{VaR}$ and
$\mathsf{CoVaR}$333Although statistical comparisons between the assumption of
stationary regressors against regressors generated via a LUR parametrization
can be constructed more formally using goodness-of-fit testing (e.g., see Chan
and Zhang, (2013)); we leave these considerations for future research..

We focus on the IVX estimators in each estimation stage, since it is well-
known to be robust to the abstract degree of persistence (e.g., see Lee,
(2016)). As a result, the doubly corrected IVX estimator (which is the IVX
estimator obtained from the second step estimation procedure), is the main
parameter of interest in terms of statistical inference and asymptotic theory.
Specifically, we establish the asymptotic properties of the doubly corrected
IVX estimator which verify the mixed Gaussianity property of the limiting
distribution regardless of the degree of persistence in both estimation
stages. Moreover, we consider a suitable correction to the expression of the
asymptotic variance-covariance matrix in the second estimation step, which is
adjusted in order to account for the first-step estimation error that produces
the generated regressor under nonstationarity.

Regarding the statistical properties of parameters obtained from the first-
stage procedure those are employed for the $\mathsf{VaR}$ risk measure
forecast and in practice parameter consistency is not violated (see, Adrian
and Brunnermeier, (2016)) regardless of the LUR parametrization. In addition,
the particular generated regressor that corresponds to fitted values of a
stationary regressand can be treated as the stationary component of the vector
of regressors used during the second stage procedure ($\mathsf{CoVaR}$).
However, the generated regressors are included in the QPR model of the second
estimation step with some sampling error, which introduces additional noise
into the asymptotic variance-covariance matrix of the model coefficients.
Therefore, the sampling error from the first stage contaminates the second
stage estimation procedure (see, Chen et al., (2021)). Furthermore, during the
second-stage ($\mathsf{CoVaR}$ estimation), model estimates are obtained using
the IVX estimator in which case the generated regressor is exogenously
determined, this implies that it can only be assumed to be stationary
(although is modelled as the same lag as the regressand). As a result, the
usual way of calculating the QR variance-covariance matrix fails to account
for the additional source of error and an adjustment is required. To sum-up,
the generated regressor is a proxy of the estimated $\mathsf{CoVaR}$ when the
nonstationary properties of regressors are included, which implies that this
proxy is treated as a nonstochastic stationary regressor in a quantile
predictive regression model with nonstationary regressors. Thus, accounting
for the presence of generated regressor bias, inference is robustify against
both the abstract degree of persistence that these nonstationary regressors
exhibit as well as to the presence of any generated regressor bias.

All random elements are defined within a probability space, denoted by
$\left(\Omega,\mathcal{F},\mathbb{P}\right)$. All limits are taken as
$n\to\infty$, where $n$ is the sample size. The symbol $"\Rightarrow"$ is used
to denote the weak convergence of the associated probability measures as
$n\to\infty$. The symbol $\overset{d}{\to}$ denotes convergence in
distribution and $\overset{\text{plim}}{\to}$ denotes convergence in
probability, within the probability space (see, Billingsley, (1968), Hall and
Heyde, (1981) and Davidson, (1994)).

The rest of the paper is organized as follows. Section 2, introduce the
econometric model and main assumptions of our framework. Section 3, presents
the estimation procedure for the generated regressor in the case of a linear
predictive regression with a conditional mean functional form while Section 4
discuss the conditional quantile case. Section 5 concludes. All technical
proofs can be found in the Appendix of the paper.

###  Main contributions

  * (i).

We focus on estimating systemic risk forecasts that incorporate information
related to the degree of persistence in macroeconomic regressors using the LUR
parametrization. The proposed dynamic quantile regression models are employed
for forecasting the pair $(\mathsf{VaR},\mathsf{CoVaR})$, based on the
definitions of Adrian and Brunnermeier, (2016) and the statistical
elicitability property of these risk measures (see, Patton et al., (2019),
Frongillo and Kash, (2021) and Fissler and Ziegel, (2016)). Moreover, the loss
function which characterizes M-estimators allows to obtain model estimates and
thus can directly affect the robustness of forecasts (see, Dimitriadis et al.,
(2022)); we leave aspects to forecast evaluation for future research.

  * (ii).

Our first main contribution to the literature is to introduce dynamic co-
quantile predictive regression models with nonstationary regressors
($\mathsf{CoQR-LUR}$). Our generalized class of estimators is employed to
jointly estimate the vector $(\mathsf{VaR},\mathsf{CoVaR})$ by incorporating
the persistence properties of regressors. Our proposed general class of models
nests both the predictive co-quantile regressions in the spirit of Adrian and
Brunnermeier, (2016) as well as the quantile predictive regressions with
persistence covariates in the spirit of Lee, (2016).

  * (iii).

Our second main contribution is to propose an estimator for the model
parameters and derive its large sample properties. The main technical
difficulty to overcome in developing related asymptotic theory, is that the
joint estimation of the risk measures of $(\mathsf{VaR},\mathsf{CoVaR})$ fails
to be elicitable, such that no real-valued scoring function exists that is
uniquely minimized by the true values of these risk measures (see, Dimitriadis
and Hoga, (2022)). This pitfall suggests the implementation of a two-step
M-estimator (see, also Patton et al., (2019)). Our proposed modelling approach
allows to extend to the case of nonstationary data, and thus is a more general
approach than current estimation methods.

  * (iv).

We focus on showing consistency and Asymptotic Mixed Gaussianity of the
proposed doubly IVX robust estimator and propose a valid inference procedure
based on consistent estimation of the asymptotic variance-covariance matrix.
The particular approach allows us to develop a unified asymptotic theory where
we take into account both the persistence properties of regressors as well as
an adjustement for the variance-covariance matrix to account for the
estimation uncertainty of standard errors occured during the first stage
procedure444Measures of information and efficiency are less popular in the
literature when regressions involve nonstationary data. However, formal tests
of the joint statistical significance of these nearly integrated regressors as
predictors for systemic risk are not so well developed in the literature
especially under the presence of a generated covariate..

The use of endogenously generated instruments (IVX filter) is employed for the
purpose of filtering out the unknown form of persistence in the regressors of
each quantile predictive regression used in the two-stage estimation
procedure. In practice, the ideal IV has to meet two requirements: it is
correlated with the endogenous regressor via an explainable and validated
relationship (relevance relationship), yet uncorrelated with the structural
error (exclusion restriction).

Regarding identification issues of the risk measure pair
$(\mathsf{VaR},\mathsf{CoVaR})$ are discussed in several studies such as
Adrian and Brunnermeier, (2016), Härdle et al., (2016), Patton et al., (2019).
Moreover, Fissler and Hoga, (2023) mention that although there is no real-
valued scoring function associated with risk measure pairs and propose to
consider a $\mathbb{R}^{2}-$valued scoring function in the closely related
context of forecast evaluation. Nevertheless, our estimation method does not
rely on the score function estimation approach, thereby the statistical
elicitability property of the risk measure pair facilitates econometric
estimation. We also avoid using copula correction methods (see, Yang et al.,
(2022) among others), since within the setting of nonstationary predictive
regression models, the IVX filtration of Phillips and Magdalinos, (2009) has
been found to be robust to abstract degree of persistence. In addition to the
proposed two-stage algorithm for the joint estimation of the risk measure
pair, as well as asymptotic theory analysis of estimators and their covariance
matrices, our framework allows to study statistical inference issues which we
leave for future research555The related statistical inference aspects we refer
to here is the subvector inference problem which is commonly examined in the
linear instrumental variables regression model. In that case, our objective is
to develop a powerful test such as the asymptotic null rejection probability
is controlled uniformly over a parameter space. The inference problem consists
of making inference on the coefficients of structural models, but in the
context of systemic risk modelling. Such statistical problems usually give
attention to aspects such as developing robust methods in the presence of weak
instruments or powerful subvector testing in IV regressions under conditional
heteroscedasticity; we give attention on the simultaneous modelling of risk
measures which involve the use of a generated covariate. Therefore, drawing
inference on individual coefficients of the parameter vector of the structural
equation for estimating the $\mathsf{CoVaR}$ risk measure, motivates us to
propose a subvector testing approach which is based on marginal conditional
regression estimation and testing with nonstationary regressors (see, Hsiao,
(1997) and De Jong, (2001)). .

We propose a practical two-stage estimation procedure666We study a broad class
of semiparametric problems, where a doubly robust estimating function
$h_{\uptau}(U;\theta,\theta_{1},\theta_{2})$, where $\uptau\in(0,1)$ for the
parameter of interest $\theta$ is available as follows. Consider two working
models denoted with $g_{1}(\boldsymbol{x}_{t-1};\theta_{1})$ and
$g_{2|1}(\boldsymbol{x}_{t-1};\theta_{2})$ for unknown functions
$y_{1t}=g_{1}^{*}(\boldsymbol{x}_{t-1})$ and
$y_{2t}=g_{2}^{*}(\boldsymbol{x}_{t-1})$. Furthermore, the estimating function
(in our case it represents the tail estimates for the $\mathsf{CoVaR}$ risk
measure), $h_{\uptau}(.;.)$ is unbiased such that
$\mathbb{E}\bigg{[}h_{\uptau}\bigg{(}U;\theta,\theta_{1},\theta_{2}\bigg{)}\bigg{]}=0$,
when $\theta$ is set to the true value $\theta^{*}$, and either $\theta_{1}$
or $\theta_{2}$, but not necessarily both, is set to the true value
$\theta_{1}^{*}$ or $\theta_{2}^{*}$ defined respectively such that
$g_{1}^{*}(\boldsymbol{x}_{t-1})\equiv g_{1}(\boldsymbol{x}_{t-1})$ or
$g_{2}^{*}(\boldsymbol{x}_{t-1})\equiv g_{2}(\boldsymbol{x}_{t-1})$ if model
$g_{1}(\boldsymbol{x}_{t-1};\theta_{1})$ or
$g_{2}(\boldsymbol{x}_{t-1};\theta_{2})$ is correctly specified. In general,
doubly robust estimation using the estimating function
$h_{\uptau}(U;\theta,\theta_{1},\theta_{2})$ consists of two stages: some
estimators $\left(\hat{\theta}_{1},\hat{\theta}_{2}\right)$ are first defined,
and then $\hat{\theta}$ is defined by solving the estimating equation
$\tilde{\mathbb{E}}\big{[}h_{\uptau}\big{(}U;\theta,\hat{\theta}_{1},\hat{\theta}_{2}\big{)}\big{]}=0$,
where $\tilde{\mathbb{E}}[\cdot]$ denotes a sample average. for estimating the
model parameters of the quantile predictive regression with nonstationary
regressors. Then, the second step employs the generated regressors in the form
of fitted values, along with a different set of nonstationary regressors to a
quantile predictive regression model. As a result, the doubly IVX corrected
estimator, is the IVX-based estimator corresponding to the second step
estimation procedure and is the main parameter of interest in terms of
statistical inference and limit theory. Since the generated regressor is
assumed to be estimated using a different functional form as well as
nonstationary regressors, then any dependence structure between the generated
regressors is captured by the variance covariance matrix of the second-step
estimation procedure. Within our econometric environment the proposed two-
stage (semiparametric) estimation method, depends on possibly nonstationary
regressors, therefore relevant parameter restrictions on $\vartheta_{1}$ and
$\vartheta_{2}$ (see, also Katsouris, (2021); Katsouris, 2023c ), which
correspond to the "working" models, are necessary to ensure robust statistical
inference.

For notational convenience consider the following two functional forms:

$\displaystyle g^{*}_{1}(\boldsymbol{x}_{t-1})$
$\displaystyle=g_{1}\big{(}\boldsymbol{x}_{t-1};\boldsymbol{\beta}_{1}\big{)}=\psi_{(1)\uptau}\left(\boldsymbol{\beta}_{1}^{\top}\xi(\boldsymbol{x}_{t})\right),$
(1.1) $\displaystyle g^{*}_{2}(\boldsymbol{x}_{t-1})$
$\displaystyle=g_{2}\big{(}\boldsymbol{x}_{t-1};\boldsymbol{\beta}_{2}\big{)}=\psi_{(2)\uptau}\left(\boldsymbol{\beta}_{2}^{\top}\xi(\boldsymbol{x}_{t})\right),$
(1.2)

where $\xi(\boldsymbol{x}_{t})$ is $p\times 1$ vector of known functions on
$\mathcal{X}$ such as
$\xi(\boldsymbol{x}_{t})=\big{(}1,\boldsymbol{x}_{t-1}^{\top}\big{)}$, and
$\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}$ are vectors of unknown
coefficients, which correspond to the two econometric specifications presented
in the framework proposed by Adrian and Brunnermeier, (2016). Within the
framework of LUR parametrization,
$(\hat{\boldsymbol{\beta}}_{1},\hat{\boldsymbol{\beta}}_{2})$ we focus on the
least-squares and the endogenous instrumentation (IVX) method of Phillips and
Magdalinos, (2009). Our implementation is, however, a plug-in estimation
approach and a novel methodology to existing approaches discussed in the
literature. Therefore, we consider the asymptotic theory analysis for
$\big{(}\hat{\boldsymbol{\beta}}_{1},\hat{\boldsymbol{\beta}}_{2}\big{)}$ and
$\widehat{\vartheta}\big{(}\hat{\boldsymbol{\beta}}_{1},\hat{\boldsymbol{\beta}}_{2}\big{)}$,
where $\boldsymbol{\beta}_{1}$ and $\boldsymbol{\beta}_{2}$ are fixed-
dimensional as the sample size $n$ grows. In particular, we assume that
$\hat{\boldsymbol{\beta}}_{1}$ is consistent for $\boldsymbol{\beta}^{*}_{1}$
when the first working model is correctly specified, and
$\hat{\boldsymbol{\beta}}_{2}$ is consistent for $\boldsymbol{\beta}^{*}_{2}$
if the second working model is correctly specified (absence of functional form
misspecification).

Assume the existence of any two pair of random variables
$\left\\{y_{1t},y_{2t}\right\\}$ and a set of predictors $\boldsymbol{x}_{t}$.
Another way of viewing our econometric environment is to consider the example
of two-countries (or firms across two countries) where we have a balanced set
of $p$ nonstationary macroeconomic variables for both economies (or firms),
respectively. In particular, let $\boldsymbol{x}_{1t}$ be the $(p\times 1)$
vector of domestic variables and $\boldsymbol{x}_{2t}$ be the $(p\times 1)$
vector of the same variables for the foreign country where each of these
possible nonstationary processes have the following LUR parametrizations

$\displaystyle\boldsymbol{x}_{1t}$
$\displaystyle=\boldsymbol{R}_{1n}\boldsymbol{x}_{1t-1}+\boldsymbol{u}_{1t}$
(1.3) $\displaystyle\boldsymbol{x}_{2t}$
$\displaystyle=\boldsymbol{R}_{2n}\boldsymbol{x}_{2t-1}+\boldsymbol{u}_{2t}$
(1.4)

In addition, we have two output variables such that $(y_{1t},y_{2t})$ which
can be predicted based on these set of nonstationary predictors. However, in
addition to the regressors that correspond to each economy, we assume that we
are also interested in another regressor which is proxy for the $\mathsf{VaR}$
which represents a quantile risk event on the macroeconomic state of the
foreign economy when predicting the current state of the home economy.
Furthermore, the IVX instrumentation is applied only to the nonstationary
regressors $\boldsymbol{x}_{1t}$ and $\boldsymbol{x}_{2t}$ for each
autoregressive equation separately. Therefore, when developing the asymptotic
theory of the (long-run) quantile-dependent parameter we can use a self-
instrumentation approach for the generated regressors that correspond to the
systemic risk proxy of each economy (as it is usually done via the
specifications of Adrian and Brunnermeier, (2016) and Härdle et al., (2016)).

The steps below should not be confused with the first or second stage
estimators in instrumental variable regressions. These procedures refer to the
two different predictive regression models. The main challenge in this step,
is that to obtain a forecast for the $\mathsf{CoVaR}$, using a quantile
predictive regression model, we need to use both the set of nearly integrated
regressors as well as the generated regressor which corresponds to the
estimated $\mathsf{VaR}$ under nonstationarity. However, the presence of this
additional covariate, implies that we have more variables than the available
instruments (endogenously generated instruments via the IVX filtration). In
particular, to overcome this issue, Demetrescu, (2014) propose a framework for
inference in an augmented predictive regression. In that case, the author
propose to instrument that additional covariate with the same process itself.
We follow a similar approach in our setting although with major necessary
modifications.

###  Two-Stage Estimation Algorithm

Our proposed two-stage estimation procedure is summarized as below:

  * Step 1.

Let $\boldsymbol{\theta}_{j}$ be the parameter vector that corresponds to both
the nonstationar regressors and the generated covariate. Estimate
$\boldsymbol{\theta}_{j}$ and compute the fitted values such that
$\widehat{x}_{tj}=g_{j}\left(\boldsymbol{w}_{tj},\widehat{\boldsymbol{\theta}}_{j}\right)$
for $j\in\left\\{1,...,N\right\\}$ and then obtain the generated regressors,
where $N$ denotes the number of nodes in the network or the number of firms in
a cross-section777Note that the main idea of modelling systemic risk based in
this approach is proposed by Katsouris, (2021) and Katsouris, 2023c  which is
also related to various studies in the literature of financial connectedness
(see, Billio et al., (2012), Diebold and Yilmaz, (2012) and Diebold and
Yılmaz, (2014) among others. . We denote with $g(\cdot)$ the optimization
function which corresponds to the conditional quantile functional form.

  * Step 2.

Compute the unknown parameter vector $\boldsymbol{\beta}$ from the following
QR model

$\displaystyle\widehat{\boldsymbol{\beta}}(\uptau)=\underset{\beta\in\mathbb{R}^{p}}{\mathsf{arg\
min}}\sum_{i=1}^{n}\rho_{\uptau}\bigg{(}y_{t}-\boldsymbol{x}_{t-1}^{\prime}\boldsymbol{\beta}\bigg{)}$
(1.5)

where the check function is defined as
$\rho_{\uptau}(\mathsf{u}):=\mathsf{u}\big{[}\uptau-\mathbf{1}\left\\{\mathsf{u}\geq
0\right\\}\big{]}$ and $\boldsymbol{x}_{t}$ are LUR processes. One uses the
estimates of $\boldsymbol{\theta}_{j}$, denoted by
$\widehat{\boldsymbol{\theta}}_{j}$ for $j\in\left\\{1,...,N\right\\}$, during
the first-stage, to obtain the generated regressor
$\widehat{x}^{\mathsf{VaR}}_{t}$.

  * Step 3.

Compute the unknown parameter vector $\boldsymbol{\beta}^{\star}$ from the
following QR model

$\displaystyle\widehat{\boldsymbol{\beta}}^{\star}(\uptau)=\underset{\boldsymbol{\beta}^{\star}\in\mathbb{R}^{p+1}}{\mathsf{arg\
min}}\sum_{i=1}^{n}\rho_{\uptau}\bigg{(}y_{t}-\tilde{\boldsymbol{x}}_{t-1}^{\prime}\boldsymbol{\beta}^{\star}\bigg{)},\
\ \ \tilde{\boldsymbol{x}}_{t-1}=\big{[}\boldsymbol{x}_{t-1}\ \ \
\widehat{x}^{\mathsf{VaR}}_{t}\big{]}$ (1.6)

where $\boldsymbol{\beta}^{\star}$ corresponds to the parameter vector of the
nonstationary regressor of the second-stage as well as the unknown parameter
that correspond to the fitted-values of the $\mathsf{VaR}$888In the framework
of Härdle et al., (2016), this coefficient corresponds to the unknown model
parameter $\hat{\beta}_{j|i}$ in their expression (6). See also expressions
(7) and (8) of Härdle et al., (2016) that correspond to the econometric
specifications used to obtain the estimates of the risk measure pair under the
assumption of stationary regressors. that have been estimated during the
first-stage estimation procedure.

Recall that we can consistently estimate the conditional quantile based on the
quantile-dependent model parameter
$\widehat{\boldsymbol{\beta}}^{\star}(\uptau)$ such that
$\widehat{Q}_{\uptau}(\uptau|\boldsymbol{x}_{t-1},\widehat{x}^{\mathsf{VaR}}_{t})=\tilde{\boldsymbol{x}}_{t-1}^{\prime}\widehat{\boldsymbol{\beta}}^{\star}(\uptau)$
(see, Koenker and Bassett, (1978)).

Therefore, one uses the estimates of $\boldsymbol{\theta}_{j}$, denoted by
$\widehat{\boldsymbol{\theta}}_{j}$, in the first step, to obtain the
generated regressor $\widehat{\boldsymbol{x}}_{t}$. Our aim is to demonstrate
that our proposed methodology achieve rate and model double robustness
simultaneously, provided that the parameter $\boldsymbol{\theta}$ satisfies
certain regularity conditions. In the current study the doubly IVX corrected
estimation approach corresponds to conditional quantile specification forms
for nonstationary data.

The remained of this paper is organized as follows. In Section 2, we introduce
the econometric model and the main assumptions of our framework. In Section 3
we present our proposed doubly corrected estimation method in the case of
conditional mean functional forms while in Section 4 discusses the conditional
quantile case which corresponds to the estimation of the risk measures under
time series nonstationarity. We conclude in Section 5.

Throughout the paper, for any vector $\boldsymbol{a}$, we denote with
$\left\lVert\boldsymbol{a}\right\rVert_{1}=\sum_{j}|a_{j}|$ and with
$\left\lVert\boldsymbol{a}\right\rVert_{2}=\sqrt{\sum_{j}|a_{j}|^{2}}$. For
any real arbitrary matrix $\boldsymbol{A}$, its Frobenius norm is defined by
$\left\lVert\boldsymbol{A}\right\rVert=\sqrt{\mathsf{trace}(\boldsymbol{A}^{\prime}\boldsymbol{A})}$.
The spectral norm of $\boldsymbol{A}$ is denoted by
$\left\lVert\boldsymbol{A}\right\rVert_{2}$, such that,
$\left\lVert\boldsymbol{A}\right\rVert_{2}=\sqrt{\lambda_{1}(\boldsymbol{A}^{\prime}\boldsymbol{A})}$,
is maximum column sum norm is
$\left\lVert\boldsymbol{A}\right\rVert_{1}=\left\lVert\boldsymbol{A}\right\rVert_{1}=\mathsf{max}_{1\leq
j\leq n}\sum_{i=1}^{m}|A_{ij}|$ and its maximum row sum norm is
$\left\lVert\boldsymbol{A}\right\rVert_{\infty}=\left\lVert\boldsymbol{A}\right\rVert_{\infty}=\mathsf{max}_{1\leq
i\leq n}\sum_{i=1}^{m}|A_{ij}|$. The operator $\overset{P}{\to}$ denotes
convergence in probability, and $\overset{D}{\to}$ denotes convergence in
distribution.

##  Econometric Model and Assumptions

###  Main Assumptions

Consider the linear predictive regression model formulated as below

$\displaystyle y_{t}$
$\displaystyle=\mu+\boldsymbol{\beta}^{\prime}\boldsymbol{x}_{t-1}+u_{t},\ \ \
1\leq t\leq n$ (2.1) $\displaystyle\boldsymbol{x}_{t}$
$\displaystyle=\boldsymbol{R}_{n}\boldsymbol{x}_{t-1}+\boldsymbol{v}_{t}$
(2.2)

where $\boldsymbol{x}_{t}\in\mathbb{R}^{p}$ is a $p-$dimensional vector and
$\boldsymbol{R}_{n}=\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}_{p}}{n^{\gamma}}\right)$
with $\gamma=1$.

We discuss the predictive regression model with a conditional mean functional
form, although for estimating the risk measures pair
$(\mathsf{VaR},\mathsf{CoVaR})$ a conditional quantile specification form is
necessary. The asymptotic behaviour of estimators under moderate deviations
from unity is beyond the scope of our study (e.g., see, Phillips and
Magdalinos, (2007), Buchmann and Chan, (2007)). In addition, several studies
in the literature consider methodologies for uniform inference in predictive
regression models such as Mikusheva, (2007), Phillips and Magdalinos, (2007),
Magdalinos, (2020), Magdalinos and Petrova, (2022) and references therein. On
the other hand, the use of the quantile predictive regression is necessary for
estimating risk measures.

######  Assumption 1.

Let
$\boldsymbol{\epsilon}_{t}=\left(u_{t},\boldsymbol{e}^{\prime}_{t}\right)^{\prime}$
be an $\mathbb{R}^{p+1}-$valued martingale difference sequence with respect to
the filtration
$\mathcal{F}_{t}=\sigma\big{(}\boldsymbol{\epsilon}_{t},\boldsymbol{\epsilon}_{t-1},...\big{)}$
with
$\mathbb{E}\big{[}\boldsymbol{\epsilon}_{t}\boldsymbol{\epsilon}_{t}^{\prime}\mid\mathcal{F}_{t-1}\big{]}=\boldsymbol{\Sigma}_{\epsilon}>0$.
Let $\boldsymbol{v}_{t}$ be an $\mathbb{R}^{p}-$valued stationary linear
process such that

$\displaystyle\boldsymbol{v}_{t}=\sum_{j=0}^{\infty}\boldsymbol{F}_{j}\boldsymbol{e}_{t-j},\
\ \ \text{with}\ \ \
\sum_{j=0}^{\infty}\left\lVert\boldsymbol{F}_{j}\right\rVert<\infty,$

where $\boldsymbol{F}_{j}$ is a sequence of constant matrices such that
$\boldsymbol{I}_{0}=\boldsymbol{I}_{p}$,
$\displaystyle\sum_{j=0}^{\infty}\boldsymbol{F}_{j}$ has full rank.

Borrowing from the literature in the case when $\sigma_{uv}\neq 0$ (which is
the case we consider), the exact OLS bias of $\hat{\beta}$ is computed from
the predictive regression model,
$\mathbb{E}\left[\hat{\beta}-\beta\right]=\delta\mathbb{E}\left[\hat{\rho}-\rho\right]$
where $\hat{\rho}$ is the OLS estimate of $\rho$ and
$\delta:=\sigma_{uv}/\sigma_{v}^{2}$ is the slope coefficient of in a
regression of $u_{t}$ on $v_{t}$. Furthermore, since $\hat{\rho}$ is known to
be downward biased in small-samples, and $(u_{t},v_{t})^{\prime}$ are
typically strongly negatively contemporaneously correlated, the autoregressive
OLS bias feeds into the small-sample distribution of $\hat{\beta}$ causing
over-rejections of the null hypothesis of no predictability, $H_{0}:\beta=0$.
Moreover, notice even with the use of a possible finite-sample-bias correction
on the OLS estimate, this reduces the noncentrality of the limiting
distribution of the OLS t-statistic, but the distribution remains nonstandard
in the near-integrated case.

######  Remark 1.

Assumption 1 imposes a condition of a martingale difference sequence
(regressor predeterminedness) for the innovation term $u_{t}$ and allows to
impose further assumptions regarding the modelling of conditional
heteroscedasticity via volatility processes such as ARCH/GARCH (see, Phillips
and Magdalinos, (2009)). A Garch process is a typical example of a strong
mixing process with exponential tails (see, Meitz and Saikkonen, (2023)).

Under Assumption 1 a functional central limit theorem holds (Phillips and
Solo, (1992))

$\displaystyle\frac{\boldsymbol{x}_{\lfloor nr\rfloor}}{\sqrt{n}}\Rightarrow
J_{c}(r):=\int_{0}^{r}e^{(r-s)\boldsymbol{C}}d\boldsymbol{W}(s),$ (2.3)

which satisfies the following stochastic differential equation

$\displaystyle
d\boldsymbol{J}_{c}(r)=c\boldsymbol{J}_{c}(r)dr+d\boldsymbol{W}(r),\ \ \
r\in(0,1).$

Relevant notation can be found in several studies in the literature such as
Dickey and Fuller, (1981), Jansson and Moreira, (2006), Mikusheva, (2012) and
Demetrescu and Rodrigues, (2020) among others. Furthermore, to avoid confusion
regarding the properties of estimators based on a conditional mean functional
form vis-a-vis a conditional quantile functional form, we include the relevant
assumptions in this section. A similar approach of comparing estimation
approaches across these two different modelling methods is presented by Liu
and Phillips, (2023), although in our framework this is might be practical
given the gap in the literature on the relevant econometric issues we aim to
tackle.

######  Assumption 2.

The distribution function of innovations $u_{t}$, $F(u)$ admits a continuous
density function $f(u)$ away from zero on the domain
$\mathcal{U}:=\left\\{u:0<F(u)<1\right\\}$.

######  Assumption 3.

Let
$\big{\\{}\mathbb{P}\left(Y_{t}<y|\boldsymbol{X}_{t}=x\right),y\in\mathbb{R},x\in\mathbb{R}^{p}\big{\\}}$
be the conditional distribution function, denoted by $F(\cdot)$ which is
absolutely continuous with respect to Lebesgue measure on $\mathbb{R}$ with
continuous density functions $f_{i}(\cdot)\equiv f(\cdot|\boldsymbol{x}_{i})$,
where for each $i=1,2,...$ $f_{i}(\cdot)$ is uniformly bounded away from zero
and infinity.

Then, the conditional quantile function of $y_{t}$ denoted with
$\mathsf{Q}_{y_{t}}\left(\uptau|\mathcal{F}_{t-1}\right)$, replaces the
conditional mean function of the predictive regression which implies the
following model specification

$\displaystyle\mathsf{Q}_{y_{t}}\left(\uptau|\mathcal{F}_{t-1}\right):=F^{-1}_{y_{t}|\boldsymbol{x}_{t-1}}(\uptau)\equiv\alpha(\uptau)+\boldsymbol{\beta}(\uptau)^{\prime}\boldsymbol{x}_{t-1}.$
(2.4)

such that
$F_{y_{t}|\boldsymbol{x}_{t-1}}(\uptau):=\mathbb{P}\big{(}y_{t}\leq\mathsf{Q}_{y_{t}}\left(\uptau|\mathcal{F}_{t-1}\right)\big{|}\mathcal{F}_{t-1}\big{)}\equiv\uptau$,
where $\uptau\in(0,1)$ is some quantile level in the compact set $(0,1)$.
Therefore, in order to define the innovation structure that corresponds to the
quantile predictive regression, we employ the piecewise derivative of the loss
function such that
$\psi_{\uptau}(\mathsf{u})=\big{[}\uptau-\boldsymbol{1}\left\\{\mathsf{u}<0\right\\}\big{]}$.
Consequently, this implies that $u_{t}(\uptau):=u_{t}-F^{-1}_{u}(\uptau)$
where $F^{-1}_{u}(\tau)$ denotes the unconditional $\uptau-$quantile of the
error term $u_{t}$. Then, the corresponding invariance principle for the
nonstationary quantile predictive regression model is formulated as below

$\displaystyle\frac{1}{\sqrt{n}}\sum_{t=1}^{\lfloor
nr\rfloor}\begin{bmatrix}\psi_{\uptau}\big{(}u_{t}(\uptau)\big{)}\\\
\boldsymbol{v}_{t}\end{bmatrix}\Rightarrow\begin{pmatrix}B_{\psi_{\uptau}}(r)_{(1\times
n)}\\\ \boldsymbol{B}_{v}(r)_{(p\times
n)}\end{pmatrix}\equiv\mathcal{BM}\begin{bmatrix}\uptau(1-\uptau)&\boldsymbol{\sigma}_{\psi_{\uptau}v}^{\prime}\\\
\boldsymbol{\sigma}_{v\psi_{\uptau}}&\boldsymbol{\Omega}_{vv}\end{bmatrix}$
(2.5)

######  Assumption 4.

Under Assumption 2 and 3, then the following conditions for the innovation
sequence hold (see, Katsouris, 2023d ):

  * (i)

The sequence of stationary conditional probability distribution functions
(pdf) denoted with $\big{\\{}f_{u_{t}(\uptau),t-1}(.)\big{\\}}$ evaluated at
zero with a non-degenerate mean function such that
$f_{u_{t}(\uptau)}(0):=\mathbb{E}\left[f_{u_{t}(\uptau),t-1}(0)\right]>0$
satisfies a FCLT given as below

$\displaystyle\frac{1}{\sqrt{n}}\sum_{t=1}^{\lfloor
nr\rfloor}\big{(}f_{u_{t}(\uptau),t-1}(0)-\mathbb{E}\left[f_{u_{t}(\uptau),t-1}(0)\right]\big{)}\Rightarrow
B_{f_{u_{t}(\uptau)}}(r).$ (2.6)

  * (ii)

For each $t$ and $\uptau\in(0,1)$, $f_{u_{t}(\uptau),t-1}(.)$ is uniformly
bounded away from zero with a corresponding conditional distribution function
$F_{t}(.)$ which is absolutely continuous with respect to Lebesgue measure on
$\mathbb{R}$ (see, Goh and Knight, (2009) and Lee, (2016)).

where

$\displaystyle\underset{0\leq r\leq
1}{\mathsf{sup}}\left|\frac{1}{n^{1-\delta}}\sum_{t=1}^{\lfloor
nr\rfloor}\left[f_{u_{t\uptau},t-1}(0)-f_{u_{\uptau}}(0)\right]\right|=o_{p}(1).$

###  IVX filtration

A solution to the endogeneity problem discussed above, which provides robust
statistical inference, is the IVX filtration proposed by Phillips and
Magdalinos, (2009). The IVX instrumentation ensures instrument relevance while
controlling for the abstract degree of persistence (see, Kostakis et al.,
(2015)) and is constructed as below

$\displaystyle z_{t}:=\sum_{j=0}^{t-1}\rho_{z}^{j}\Delta x_{t-j},\ \
\rho_{z}=\left(1-\frac{c_{z}}{n^{\gamma_{z}}}\right),c_{z}>0\ \text{and}\
\gamma_{z}\in(0,1).$ (2.7)

The IVX filtration for $x_{t}$, is chosen so that $z_{t}$ corresponds to a
mildly integrated process when the predictor $x_{t}$ is nearly integrated.
Therefore, the IVX estimator of $\beta$ is found to have a slower convergence
rate than the conventional OLS estimator under near integration, such that
$n^{\frac{1+\gamma_{z}}{2}}$ (see, also Katsouris, 2023b  and Duffy and
Kasparis, (2021)). The asymptotic theory of the IVX estimator can be found to
be mixed Gaussian irrespectively of the degree of endogeneity implied by
$\delta$, leading to standard inference in $t$ and Wald tests. Moreover, it
can be proved that under low persistence, the IVX estimator is asymptotically
equivalent to the OLS procedure. Our main focus in this paper is the
estimation of the risk measures of $\mathsf{VaR}$ and $\mathsf{CoVaR}$ using
nearly integrated regressors, which is a novel aspect in the literature as
conventional estimation approaches mainly consider that macroeconomic
variables are stationary. Specifically, we start by considering the estimation
of the $\mathsf{VaR}$ with a nearly integrated regressor. A vast literature
has previously proposed robust methodologies for estimating the parameter of
interest using the OLS approach and found evidence of significant size
distortions when the model includes high peristence regressors.

To obtain the IVX estimator we write the corresponding model with the demeaned
variates denoted with
$y^{\mu}_{t}=\boldsymbol{\beta}^{\prime}\boldsymbol{x}^{\mu}_{t-1}+v^{\mu}_{t}$,
where the sample moments are estimated as
$\bar{y}_{t}=\left(y_{t}-\frac{1}{n-1}\sum_{t=1}^{n}y_{t}\right)$. Then, the
IVX estimator is constructed

$\displaystyle\widehat{\boldsymbol{\beta}}^{ivx}=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}y_{t}\boldsymbol{z}_{t-1}^{\prime}\right).$
(2.8)

Denote with $\widehat{\boldsymbol{\beta}}^{ols}$ the OLS estimator and define
the residuals
$\widehat{u}_{t}=\bar{y}_{t}-\widehat{\boldsymbol{\beta}}^{ols}\bar{\boldsymbol{x}}_{t-1}$
and the residual variance estimator
$\widehat{\sigma}_{u}^{2}=\frac{1}{n-1}\sum_{t=1}^{n}\widehat{u}_{t}^{2}$.
Then, the IVX-Wald statistic under the null hypothesis is
$\mathcal{W}=\widehat{\boldsymbol{\beta}}^{ivx\prime}\widehat{\boldsymbol{V}}^{-1}\widehat{\boldsymbol{\beta}}^{ivx}$,
where a feasible estimator for the covariance matrix is given by

$\displaystyle\widehat{\boldsymbol{V}}=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\right)^{-1}\left(\widehat{\sigma}_{u}^{2}\sum_{t=1}^{n}\boldsymbol{z}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\right)\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\right)^{-1}$
(2.9)

if our main interest is testing the null hypothesis of no predictability.
Notice also, that an alternative to the IVX filtration which also ensures
unified inference in the case of quantile-dependent parameters is proposed by
Liu et al., (2023). Although the extension to the case of an augmented
quantile predictive regression is currently work in progress by Katsouris,
2023e .

###  Two-Stage Estimation Procedure

######  Assumption 5.

Under the prespecified Assumptions 1 to 3 and the null hypothesis
$\mathcal{H}_{0}$, the coefficients except for the intercept of the
conditional quantile are constant across different quantiles in the compact
subset $(0,1)$.

####  First-Stage Procedure

During the first stage of our procedure we obtain a consistent estimator of
the parameter vector. We focus on the IVX estimator which is found to be
robust to the abstract degree of persistence (nearly integrated regressors).
The econometric model is the linear predictive regression as below

$\displaystyle y_{1,t}$
$\displaystyle=\beta_{01}(\uptau)+\boldsymbol{\beta}_{11}^{\prime}(\uptau)\boldsymbol{x}_{1,t-1}+u_{1t}(\uptau),\
\ \ \text{for}\ \ \ t=1,...,n$ (2.10) $\displaystyle\boldsymbol{x}_{1,t}$
$\displaystyle=\boldsymbol{R}_{1,n}\boldsymbol{x}_{1,t-1}+\boldsymbol{v}_{1,t}$
(2.11)

Denote with
$\boldsymbol{\beta}(\uptau)=\left(\beta_{01}(\uptau),\boldsymbol{\beta}_{11}^{\prime}(\uptau)\right)^{\prime}$
and the corresponding IVX estimator with $\boldsymbol{\beta}^{ivx}(\uptau)$.
The asymptotic behaviour of the IVX estimator that corresponds to the linear
predictive regression model under abstract degree of regressors persistence is
studied by Kostakis et al., (2015) while the asymptotic properties of the IVX
estimator for the quantile predictive regression model is studied by Lee,
(2016). Since when estimating the risk measure pair
$(\mathsf{VaR},\mathsf{CoVaR})$, we use the conditional quantile distribution,
i.e., to capture the effect of $\boldsymbol{x}_{t-1}$ on the conditional
quantile of $\boldsymbol{y}_{t}$; however with the additional assumption of
possibly nonstationary regressors (see, also Koenker and Xiao, (2004) and
Xiao, (2009)), then estimation relies on the nonstationary quantile predictive
regression.

####  Second-Stage Procedure

During the second stage of our procedure we consider a consistent estimator of
the parameter vector for the set of regressors which includes the
nonstationary regressors as well as the generated regressor from the first
stage procedure. More precisely, the generated regressor in our study
corresponds to the fitted values of the predictive regression model based on
the IVX estimator obtained in the first stage procedure.

$\displaystyle y_{2,t}$
$\displaystyle=\beta_{02}(\uptau)+\boldsymbol{\beta}_{12}^{\prime}(\uptau)\boldsymbol{x}_{2,t-1}+\delta(\uptau)\widehat{y}_{1,t}\big{(}\widehat{\boldsymbol{\beta}}^{ivx}(\uptau)\big{)}+u_{2t}(\uptau),\
\ \ \text{for}\ \ \ t=1,...,n$ (2.12) $\displaystyle\boldsymbol{x}_{2,t}$
$\displaystyle=\boldsymbol{R}_{2,n}\boldsymbol{x}_{2,t-1}+\boldsymbol{v}_{2,t},\
\
\boldsymbol{R}_{jn}=\left(\boldsymbol{I}_{p}-\frac{\boldsymbol{C}_{jp}}{n^{\gamma}}\right),\
\gamma=1.$ (2.13)

where $\boldsymbol{C}_{jp}=\mathsf{diag}\left\\{c_{j1},...,c_{jp}\right\\}$
and $j\in\left\\{1,2\right\\}$. Define with
$\widetilde{\boldsymbol{\beta}}(\uptau)=\big{(}\beta_{02}(\uptau),\boldsymbol{\beta}_{12}^{\prime}(\uptau),\delta(\uptau)\big{)}^{\prime}$
the parameter vector of the second stage procedure and the corresponding IVX
estimator with $\widetilde{\boldsymbol{\beta}}^{ivx}(\uptau)$.

Specifically, we observe that this extended parameter vector includes both the
nearly integrated regressors as well as the generated regressor from the first
stage estimation step.

In this case, we need to develop the asymptotic distribution theory for both
the OLS and IVX estimators when the generated regressor is included in the set
of regressors in the second stage predictive regression model that corresponds
to the econometric specification of the CoVaR risk measure. A modified IVX
estimator is necessary to be developed in order to account for the presence of
the particular effect in the setting of the quantile predictive regression
models. The nuisance coefficient of persistence is defined such that $c_{i}>0$
and is either $\gamma=1$, which corresponds to near unit root regressors or
$\gamma\in(0,1)$ that corresponds to mildly integrated regressors.

Therefore, the focus of this paper is the implementation of the above
econometric environment and the development of the corresponding asymptotic
theory in the case of the quantile predictive regression model. Consider the
piecewise derivative of the loss function which is defined as below
$\psi_{\uptau}\left(u\right)=\uptau-\mathbf{1}\left(u<0\right)$. Then the
innovation sequence of the quantile predictive regression model,
$u_{t}(\uptau)=u_{t}-F_{u}^{-1}(\uptau)$, and $F_{u}^{-1}(\uptau)$ is the
unconditional $\uptau$-quantile of $u_{t}$, where $\uptau\in(0,1)$ is a fixed
quantile (see, Lee, (2016), Fan and Lee, (2019) and Katsouris, (2022)). Thus,
the $\boldsymbol{\beta}^{ivx}$ estimator is obtained from the first stage
procedure using the IVX estimation method999The reason that we apply the IVX
instrumentation procedure only to the estimate from the first stage quantile
predictive regression model rather to the corresponding generated covariate is
to ensure that those fitted values can preserve their definition as the
estimated Value-at-Risk but in our case adjusted based on the presence of
persistent predictors and corrected accordingly using the IVX methodology. for
the model coefficients based on the nonstationary regressors which implies
that the generated regressor is defined

$\displaystyle\widehat{y}_{1,t}\left(\widehat{\boldsymbol{\beta}}^{ivx}(\uptau)\right)=\boldsymbol{x}_{1,t-1}^{\prime}\widehat{\boldsymbol{\beta}}^{ivx}(\uptau)$
(2.14)

Thus to conduct inference we need to obtain a consistent estimator for the
covariance matrix of $\tilde{\boldsymbol{\beta}}_{2}^{ivx}$ using the usual
"sandwich" formula given below (see, Demetrescu and Rodrigues, (2020)).

##  Doubly Corrected Estimation Method: Conditional Mean Case

###  Residual Augmented IVX Estimation Example

We begin our asymptotic theory analysis by considering the framework for the
residual augmented IVX estimation methodology previously examined in the
literature (see, Demetrescu and Rodrigues, (2020)). The IVX instrumentation
provides a methodology such that the asymptotic distribution of the IVX-Wald
based test statistic converges to a nuisance-parameter free distribution.
Consider the feasible predictive regression defined with the following form

$\displaystyle y^{\mu}_{t}=\beta
x_{t-1}^{\mu}+\delta\hat{v}_{t}+\varepsilon_{t},\ \ \
\varepsilon_{t}\sim_{\textit{i.i.d}}\ \text{and}\ \ t=1,...,n.$ (3.1)

where $y^{\mu}_{t}$ and $x_{t-1}^{\mu}$ denote the demeaned variates,
following the notation employed in several studies in the literature.
Moreover, $\hat{v}_{t}$ represents a generated covariate which in our study
corresponds to the estimated $\mathsf{VaR}$ under nonstationarity.

The parameter of interest for inference purposes is the slope coefficient of
the quantile predictive regression of the second-stage procedure. Based on the
specification above, the IVX estimator of $\beta$, has a slower convergence
rate than the conventional OLS estimator under near integration. In addition,
the IVX estimator is mixed Gaussian in the limit regardless of the degree of
endogeneity which implies standard inference based on the Wald test. Under low
persistence it is asymptotically equivalent to the OLS procedure which implies
that the IVX estimator is easy to implement. The null hypothesis of no
predictability implies the presence of no serial correlation. Furthermore, we
conjecture that the generated regressor since it is stationary it will not
affect the convergence of the IVX estimator to its asymptotic limit, given the
lower convergence rate of IVX compared to the OLS estimator. Although our
procedure has some similarities with residual-based augmentation procedures
proposed in the literature, our objective is to demonstrate that in finite-
samples and asymptotically has similar desirable properties. All directly
observable variables and demeaned. Let
$\bar{y}_{t}=y_{t}-\bar{y}_{n}=\left(y_{t}-\frac{1}{n}\sum_{t=1}^{n}y_{t}\right)$,
$\bar{x}_{t}=x_{t}-\bar{x}_{n}=\left(x_{t}-\frac{1}{n}\sum_{t=1}^{n}x_{t}\right)$
and
$\bar{x}_{t-1}=x_{t-1}-\bar{x}_{n-1}=\left(x_{t-1}-\frac{1}{n}\sum_{t=1}^{n}x_{t-1}\right)$.
Then the IVX estimator is defined as below

$\displaystyle\widehat{\beta}^{ivx}:=\frac{\displaystyle\sum_{t=1}^{n}\bar{y}_{t}z_{t-1}}{\displaystyle\sum_{t=1}^{n}\bar{x}_{t-1}z_{t-1}},\
\ \text{with}\ \
se\left(\widehat{\beta}^{ivx}\right):=\left(\sum_{t=1}^{n}\bar{x}_{t-1}z_{t-1}\right)^{-1}\widehat{\sigma}_{u}\left(\sum_{t=1}^{n}z^{2}_{t-1}\right)^{1/2}.$

where $se\left(\widehat{\beta}^{ivx}\right)$ is the standard error of the IVX
estimator. To improve the finite-sample behaviour Kostakis et al., (2015)
suggest the use of OLS residuals $\widehat{u}_{t}$ for the computation of
$\widehat{\sigma}_{u}$, do not demean the instrument $z_{t}$ and correct the
standard errors by subtracting from
$\widehat{\sigma}_{u}\left(\sum_{t=1}^{n}z^{2}_{t-1}\right)^{1/2}$. In
addition, we define with
$\tilde{\varepsilon}_{t}:=\bar{y}_{t}-\big{(}\widehat{\beta}^{ols}\bar{x}_{t-1}+\widehat{\delta}\widehat{v}_{t}\big{)}$,
and use those instead of the IVX residuals due to the superconsistency
properties of the OLS-based residuals in the near-integrated context. Recall
that the above expressions correspond to the model estimates from a linear
predictive regression model with conditional mean functional form. Based on
the above illustrative example, we replace the regressor $\hat{v}_{t}$ with
the generated dependent variable, $\widehat{\bar{y}}_{1,t}$, which is used as
an additional covariate. We do this first in the case of a linear predictive
regression measure, so the estimation procedures does not correspond to the
risk measure pair yet.

###  Generated Regressor in Nonstationary Linear Predictive Regression

In compact form we have that

$\displaystyle\bar{y}_{2,t}=\boldsymbol{X}_{t-1}\boldsymbol{\Gamma}^{\prime}+\varepsilon_{t},\
\ \ \varepsilon_{t}\sim_{\textit{i.i.d}}\ \text{and}\ \ t=1,...,n.$ (3.2)

where $\boldsymbol{X}_{t-1}=\big{(}\bar{\boldsymbol{x}}_{t-1}\
\widehat{\bar{y}}_{1,t}\big{)}$ and
$\boldsymbol{\Gamma}=\big{(}\boldsymbol{\beta}\ \delta\big{)}$. Moreover,
define with $\boldsymbol{Z}_{t-1}=\big{(}\boldsymbol{z}_{t-1}\
\widehat{\bar{y}}_{1,t}\big{)}$ where $\boldsymbol{z}_{t-1}$ is the IVX
instrument for $\boldsymbol{x}_{t-1}$ and the generated regressor is self-
instrumented.

Notice that the generated covariate $\widehat{\bar{y}}_{1,t}$ aims to resemble
the estimated VaR under nonstationarity. In other words, the generated
regressor, denoted with $\widehat{\bar{y}}_{1,t}$ is obtained from a
predictive regression model which includes only the nonstationary predictors.
In practise, these nonstationary predictors can be different than the
nonstationary predictors corresponding to the regressand $\widehat{y}_{2,t}$
due to the proposed dependence structure. However, for simplicity we can also
consider the case in which we use the same nonstationary regressors for both
predictive regression models. In this section, we obtain derivations for the
linear predictive regression model. However, when estimating the risk measures
of $\mathsf{VaR}$ and $\mathsf{CoVaR}$, the nonstationary quantile predictive
regression models are used which implies that all parameters are quantile-
dependent for a fixed quantile $\uptau\in(0,1)$. Therefore, we are interested
for the consistent estimation of the $\beta_{2}$ estimator (e.g., via IVX) as
well as constructing hypothesis testing for: (i) only $\beta$ parameter (e.g.,
no predictability), (ii) only $\delta$ parameter (e.g., no presence of
generated regressor) and (iii) linear restrictions on all the parameters of
the predictive regression model that includes both nonstationary regressors
and the generated regressor (e.g., univariate generated regressor). The
covariance estimator for $\tilde{\boldsymbol{\beta}}_{2}^{ivx}$ expressed as

$\displaystyle\widehat{\mathsf{Cov}\left(\tilde{\boldsymbol{\beta}}_{2}^{ivx}\right)}$
$\displaystyle:=\left(\boldsymbol{B}_{n}^{-1}\right)\mathbb{V}_{n}\left(\boldsymbol{B}_{n}^{-1}\right)^{\prime},\
\ \
\boldsymbol{B}_{n}:=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\right)$
$\displaystyle\mathbb{M}_{n}$
$\displaystyle:=\left(\sum_{t=1}^{n}\boldsymbol{z}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\tilde{\varepsilon}^{2}_{t}\right)+\left[\boldsymbol{\gamma}^{\top}\otimes\left(\frac{1}{n}\sum_{t=1}^{n}\boldsymbol{z}_{t-1}\bar{\boldsymbol{x}}_{t-1}^{\top}\right)\left(\frac{1}{n}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\bar{\boldsymbol{x}}_{t-1}^{\top}\right)^{-1}\right]\left(\sum_{t=1}^{n}\boldsymbol{\nu}_{t}\boldsymbol{\nu}_{t}^{\top}\otimes\bar{\boldsymbol{x}}_{t-1}\bar{\boldsymbol{x}}_{t-1}^{\top}\right)$
$\displaystyle\times\left[\boldsymbol{\gamma}\otimes\left(\frac{1}{n}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\bar{\boldsymbol{x}}_{t-1}^{\top}\right)^{-1}\left(\frac{1}{n}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\boldsymbol{z}_{t-1}^{\top}\right)\right].$

Therefore, the just-identified model for the IV estimator of
$\boldsymbol{\Gamma}$ is given by

$\displaystyle\hat{\boldsymbol{\Gamma}}=\left(\sum_{t=1}^{n}\boldsymbol{Z}_{t-1}^{\prime}\boldsymbol{X}_{t-1}\right)^{-1}\left(\sum_{t=1}^{n}\boldsymbol{Z}_{t-1}^{\prime}\bar{y}_{2,t}\right).$
(3.3)

Furthermore, notice that the nonstationary regressors $x_{t-1}$ and the
generated regressor are orthogonal since the generated regressor can be
considered as exogenous variable to the particular model, then we can compute
the generated regressor augmented estimator and test statistic in two steps.
Notice the constructed test statistic refers to the IVX estimator. Therefore,
we focus on deriving the asymptotic distribution for the IVX estimator under
the presence of the generated regressor in the quantile predictive regression
model. The estimation procedure is implemented into two-stages, and our
asymptotic theory analysis demonstrates that although the limit converges into
a nonstandard distribution the particular limit is nuisance-parameter free.

######  Remark 2.

The key point in this estimation step, is that the variance-covariance needs
to be adjusted to account for the first-step estimation error. In other words,
when estimating the VaR using the quantile predictive regression of the first
stage, the forecast at the current time period $t$, is estimated with some
additional source of error. For this reason, an adjustement is necessary.
Consider the estimated coefficient of the quantile predictive regression model
of the first stage as

$\displaystyle Q_{\uptau}\big{(}y_{2t}|x_{t-1},y_{1t}\big{)}=\theta X_{t-1},\
\ \ \uptau\in(0,1).$ (3.4)

In other words, the estimation step needs to account for two sources of
errors, that is, the usaul estimation error in obtaining a consistent
estimator that corresponds to the IVX instrumentation of the nearly integrated
regressors, and the second source of error is the sampling error in generating
the forecast for the VaR from the first stage of the process. More
specifically, this implies that using a Bahadur representation of the QR-IVX
estimator we need to determine the precise stochastic order of the remainder
term when the generated regressor is included in the conditional quantile
specification of the model. Thus the main focus of the estimation step should
be how to correct for the effect to the overall variance due to the presence
of this generated regressor. In particular, the parameter vector from the
first stage has the usual convergence rate that the IVX estimator has
(regardless of the degree of persistence). However, we need to check whether
this convergence rate remains the same for the parameter vector of the second
stage estimation under the presence of the generated regressor. Furthermore,
the property of stochastic equicontinuity still holds in this case, since we
have a plug-in estimator, regardless if the underline stochastic processes
correspond to nonstationary time series data.

####  Estimation Procedure in Nonstationary Linear Predictive Regression Case

  * Step 1.

OLS-regress $\bar{y}_{2,t}$ on the generated regressor
$\widehat{\bar{y}}_{1,t}$ and obtain the OLS estimator below

$\displaystyle\widehat{\delta}^{ols}=\left(\sum_{t=1}^{n}\widehat{\bar{y}}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}\widehat{\bar{y}}_{1,t}\bar{y}_{2,t}\right).$
(3.5)

and compute the residual such that
$\widehat{\eta}_{t}=\bar{y}_{2,t}-\widehat{\delta}^{ols}\widehat{\bar{y}}_{1,t}$.
Here we assume that the elements of the dependent vector are not serially
correlated, which implies that $y_{1t}$ and $y_{2t}$ are uncorrelated. This
simplifies the development of the asymptotic theory, despite the presence of
endogeneity and highly persistent regressors.

  * Step 2.

IVX-regress $\bar{y}_{2,t}$ on $\bar{\boldsymbol{x}}_{t-1}$, leading to

$\displaystyle\widehat{\boldsymbol{\beta}}_{2}^{ivx}=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{t-1}\boldsymbol{z}_{t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\tilde{y}_{t}\boldsymbol{z}_{t-1}^{\prime}\right).$
(3.6)

Recall the linear predictive regression model (for node 1) is given by

$\displaystyle y_{1,t}$
$\displaystyle=\boldsymbol{\beta}_{1}^{\prime}\boldsymbol{x}_{t-1}+u_{1,t}$
(3.7) $\displaystyle\boldsymbol{x}_{t}$
$\displaystyle=\boldsymbol{R}\boldsymbol{x}_{t-1}+\boldsymbol{v}_{1,t}$ (3.8)

Then we estimate the following model

$\displaystyle\bar{y}_{2,t}$
$\displaystyle=\boldsymbol{\beta}_{2}^{\prime}\bar{\boldsymbol{x}}_{t-1}+\delta\widehat{\bar{y}}_{1,t}+u_{2,t}$
(3.9)

Moreover, by OLS-regressing $\bar{y}_{2,t}$ on $\widehat{\bar{y}}_{1,t}$ such
that $\bar{y}_{2,t}=\theta\widehat{\bar{y}}_{1,t}+\eta_{t}$, where $\eta_{t}$
is a disturbance term, we construct the transformed dependent variable
$\tilde{y}_{2,t}$ given by the following expression

$\displaystyle\tilde{y}_{2,t}:=\big{(}\bar{y}_{2,t}-\widehat{\theta}^{ols}\widehat{\bar{y}}_{1,t}\big{)}$
(3.10)

  * •

In the second stage of the estimation procedure, we need to re-estimate the
quantile kernel function when the generated regressor is included in the
design matrix, especially when constructing the Wald test under the null
hypothesis of no systemic risk in the network. To obtain a measure of the
performance of the testing procedure, we need to compare the standard
estimation procedure without the variance-covariance matrix correction against
the estimation procedure that includes the adjustment, in both cases under the
null hypothesis. Then, under the alternative hypothesis of deviations, that
is, existence of non-zero systemic risk effect in the network (see, Katsouris,
(2021); Katsouris, 2023c ), then the larger this distance is, under the
adjustement then the larger the power function of the test would be when all
the individual quantile predictive regression models of the system correspond
to a non-zero coefficient of node-specific systemic risk.

  * •

When the nonstationary regressors for the two nodes are the same then second-
order effects can contribute to the covariance of the estimator. However, due
to the proposed dependence structure (see, Katsouris, (2021); Katsouris, 2023c
), in practice we ensure that these nonstationary (or near unit roots) are not
equivalent across the different predictive regression models. Notice that the
quantity $\widehat{y}_{1,t}$ is the forecaster VaR which is obtained using
time series observations from $t$ to $(t-1)$. Therefore, it might have an
effect on the convergence of the parameter estimate $\delta$. In summary, we
need to control for two effects, that is, a bias effect on the estimate of
$\delta$ and a variance effect. In particular, the bias effect will disappear
when $\boldsymbol{x}_{t-1}$ and the regressors used to estimate
$\widehat{y}_{1,t}$ are uncorrelated. The variance effect appears due to the
choice of the generated regressors. Standard OLS estimation methods do not
correct for the particular variance effect. On the other hand, the variance
effect does not influence the inference methodology for the estimator
$\hat{\beta}_{2}$ in which case an important property to hold is to obtain a
consistent estimator regardless of the presence of nonstationary regressors.
In the derivations below, we use $\widehat{y}_{1,t}$ e.g., the estimated
$\mathsf{VaR}$ under nonstationarity (when we employ the quantile predictive
regression model).

Therefore, the IVX estimator for the second econometric specification
($\mathsf{CoVaR}$ quantile predictive regression) is obtained by the following
expression:

$\displaystyle\tilde{\boldsymbol{\beta}}_{2}^{ivx}$
$\displaystyle=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\tilde{y}_{2,t}\boldsymbol{z}_{2,t-1}^{\prime}\right)$
$\displaystyle=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bigg{[}\boldsymbol{\beta}_{2}^{\prime}\bar{\boldsymbol{x}}_{2,t-1}+\delta\widehat{y}_{1,t}+u_{2,t}-\widehat{\delta}^{ols}\widehat{y}_{1,t}\bigg{]}\boldsymbol{z}_{2,t-1}^{\prime}\right)$
$\displaystyle=\boldsymbol{\beta}_{2}-\left(\widehat{\delta}^{ols}-\delta\right)\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\widehat{y}_{1,t}\boldsymbol{z}_{2,t-1}^{\prime}\right)+\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{z}_{2,t-1}^{\prime}\right)$

which gives

$\displaystyle\bigg{(}\tilde{\boldsymbol{\beta}}_{2}^{ivx}-\boldsymbol{\beta}_{2}\bigg{)}=-\left(\widehat{\delta^{ols}}-\delta\right)\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\widehat{y}_{1,t}\boldsymbol{z}_{2,t-1}^{\prime}\right)+\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{z}_{2,t-1}^{\prime}\right)$

Notice that for the first term we have that

$\displaystyle\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\widehat{\boldsymbol{\beta}}_{1}^{ivx}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)\equiv\widehat{\boldsymbol{\beta}}_{1}^{ivx}$
(3.11)

Notice that the above expression only holds when we have identical
nonstationary regressors such that
$\bar{\boldsymbol{x}}_{1,t-1}\equiv\bar{\boldsymbol{x}}_{2,t-1}$. Therefore,
we obtain the following expression

$\displaystyle\bigg{(}\tilde{\boldsymbol{\beta}}_{2}^{ivx}-\boldsymbol{\beta}_{2}\bigg{)}$
$\displaystyle=-\left(\widehat{\delta}^{ols}-\delta\right)\times\widehat{\boldsymbol{\beta}}_{1}^{ivx}+\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{z}_{2,t-1}^{\prime}\right)$
(3.12)

######  Remark 3.

Notice that when $\boldsymbol{x}_{t-1}$ and $\widehat{y}_{1,t}$ are not
correlated, then $\widehat{\theta}^{ols}$ is a consistent estimator of
$\delta$. However, in the case we employ the same nonstationary predictors to
estimate the value of $\widehat{y}_{1,t}$ then these two quantities are
correlated and then it is not so clear to determine the correlation structure
due to the presence of different effects. On the other hand, the dependence
structure proposed in Katsouris, (2021); Katsouris, 2023c  ensures that this
case doesn’t occur, so that we can consistently estimate these two parameters
which are important for inference purposes.

Note that $\widehat{\boldsymbol{\beta}}_{1}^{ivx}$ corresponds to the IVX
estimator of the quantile predictive regression model for the $\mathsf{VaR}$
risk measure. Additionally, note that from Lemma B4 of Kostakis et al., (2015)
it holds,

$\displaystyle\frac{1}{n^{\frac{1+\gamma_{z}}{2}}}\sum_{t=1}^{n}\boldsymbol{z}_{2,t-1}u_{2,t}\Rightarrow\mathcal{N}\left(0,\sigma_{u_{2}}^{2}\times\boldsymbol{V}_{c_{z}}\right).$
(3.13)

Therefore, by incorporating the convergence rate we have that

$\displaystyle
n^{\frac{1+\gamma_{z}}{2}}\bigg{(}\tilde{\boldsymbol{\beta}}_{2}^{ivx}-\boldsymbol{\beta}_{2}\bigg{)}=n^{\frac{1+\gamma_{z}}{2}}\left(\delta-\widehat{\delta}^{ols}\right)\times\widehat{\boldsymbol{\beta}}_{1}^{ivx}+\left(\frac{1}{n^{1+\gamma_{z}}}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\frac{1}{n^{\frac{1+\gamma_{z}}{2}}}\sum_{t=1}^{n}\boldsymbol{z}_{2,t-1}u_{2,t}\right).$

Thus, for the LUR regressor case we have that the second term of expression
(3.2.1) converges to

$\displaystyle\left(\frac{1}{n^{1+\gamma_{z}}}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\frac{1}{n^{\frac{1+\gamma_{z}}{2}}}\sum_{t=1}^{n}\boldsymbol{z}_{2,t-1}u_{2,t}\right)\Rightarrow-\left(\boldsymbol{\Omega}_{uu}+\int_{0}^{1}\boldsymbol{J}^{\mu}_{c}dJ^{\prime}_{c}\right)\boldsymbol{C}_{z}^{-1}\times\mathcal{N}\left(0,\sigma_{u_{2}}^{2}\times\boldsymbol{V}_{C_{z}}\right).$

To see this, from the Appendix of KMS expression (27) gives that

$\displaystyle\left(\frac{1}{n^{1+\gamma_{z}}}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)=\bigg{(}\boldsymbol{\Omega}_{uu}+\boldsymbol{V}_{C}\boldsymbol{C}\bigg{)}\boldsymbol{C}_{z}^{-1}+o_{p}(1).$
(3.14)

where

$\displaystyle\boldsymbol{V}_{C}=\int_{0}^{\infty}e^{r\boldsymbol{C}}\boldsymbol{\Omega}_{uu}e^{r\boldsymbol{C}}dr\
\ \ \text{and}\ \ \
\boldsymbol{V}_{C_{z}}=\int_{0}^{\infty}e^{r\boldsymbol{C}_{z}}\boldsymbol{\Omega}_{uu}e^{r\boldsymbol{C}_{z}}dr.$
(3.15)

Notice that by Kostakis et al., (2015) we have that
$\boldsymbol{\Psi}_{uu}=\left(\boldsymbol{\Omega}_{uu}+\displaystyle\int_{0}^{1}\boldsymbol{J}^{\mu}_{c}dJ^{\prime}_{c}\right)$
since we consider local unit root regressors. Therefore, similar to Theorem A
(i) of KMS it holds that

$\displaystyle\left(\frac{1}{n^{1+\gamma_{z}}}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\frac{1}{n^{\frac{1+\gamma_{z}}{2}}}\sum_{t=1}^{n}\boldsymbol{z}_{2,t-1}u_{2,t}\right)\Rightarrow\mathcal{MN}\bigg{(}0,\sigma_{u_{2}}^{2}\times\left(\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\bigg{)}$

which is a mixed Gaussian random variate since the covariance matrix is a
function of the OU process. Next we investigate the asymptotic behaviour of
the first term of expression (3.2.1):

$\displaystyle\mathcal{A}:=n^{\frac{1+\gamma_{z}}{2}}\left(\delta-\widehat{\delta}^{ols}\right)\times\widehat{\boldsymbol{\beta}}_{1}^{ivx}$
(3.16)

Specifically, assuming that node 1 and 2 have a set of nonstationary
regressors that are generated from non-identical stochastic processes (i.e.,
not identical regressors), then it follows that

$\displaystyle\mathcal{A}$
$\displaystyle=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)\times\widehat{\boldsymbol{\beta}_{1}^{ivx}}$
$\displaystyle=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{2,t-1}^{\prime}\right)\times\left[\boldsymbol{\beta}_{1}+\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{1,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bar{u}_{1,t}\boldsymbol{z}_{1,t-1}^{\prime}\right)\right].$

Regarding the asymptotic behaviour of the first term of the expression we have
that the IVX estimator from the quantile predictive regression that
corresponds to the VaR risk measure is

$\displaystyle\widehat{\boldsymbol{\beta}}_{1}^{ivx}$
$\displaystyle=\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{1,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}y_{1,t}\boldsymbol{z}_{1,t-1}^{\prime}\right),\
\ \
\bar{y}_{1,t}=\boldsymbol{\beta}_{1}^{\prime}\bar{\boldsymbol{x}}_{1,t-1}+\bar{u}_{1,t}$
$\displaystyle=\boldsymbol{\beta}_{1}+\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{1,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bar{u}_{1,t}\boldsymbol{z}_{1,t-1}^{\prime}\right)$

Moreover, the number of nonstationary regressors for both econometric
specifications is the same although these two sets of regressors possibly have
different nuisance parameters of persistence. The limiting distributions for
their IVX estimators have the same dimensions such that

$\displaystyle
n^{\frac{1+\gamma_{z}}{2}}\left(\widehat{\boldsymbol{\beta}}_{1}^{ivx}-\boldsymbol{\beta}_{1}\right)$
$\displaystyle=\left(\frac{1}{n^{1+\gamma_{z}}}\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{1,t-1}^{\prime}\right)^{-1}\left(\frac{1}{n^{\frac{1+\gamma_{z}}{2}}}\sum_{t=1}^{n}\bar{u}_{1,t}\boldsymbol{z}_{1,t-1}^{\prime}\right)$
$\displaystyle\Rightarrow\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}.$

Rearranging the above expressions, we obtain the following limit result

$\displaystyle
n^{\frac{1+\gamma_{z}}{2}}\bigg{[}\bigg{(}\tilde{\boldsymbol{\beta}}_{2}^{ivx}-\boldsymbol{\beta}_{2}\bigg{)}-\boldsymbol{\beta}_{1}\bigg{]}$
$\displaystyle=\left(\delta-\widehat{\delta}^{ols}\right)\times\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}$
$\displaystyle\ \ \ \
+\mathcal{MN}\bigg{(}0,\sigma_{u_{2}}^{2}\times\left(\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\bigg{)}.$
(3.17)

######  Remark 4.

In this paper we are not specifically attempting to solve the endogeneity
issue in quantile predictive regressions models. We employ a well-investigated
method in the literature that tackles the endogeneity problem, meaning it
produces weak convergence results into standard asymptotic distributions (see,
Lee, (2016), Fan and Lee, (2019)). Alternative approaches that tackle the
endogeneity problem and produce uniform valid inference regardless of the
unknown persistence properties are proposed in the frameworks of Cai et al.,
(2023) and Liu et al., (2023). In contrast, our study provides a general
framework and establish the asymptotic properties for quantile predictive
regression models with a generated regressor, which is particularly useful
when jointly estimating the risk measure pair of
$(\mathsf{VaR},\mathsf{CoVaR})$ under the presence of time series
nonstationarity.

We suspect that regardless of the estimation method employed to robustify the
quantile-based model parameters to the unknown persistence, when estimating
the $\mathsf{CoVaR}$ based on the specifications we follow ((Adrian and
Brunnermeier,, 2016), Härdle et al., (2016)) and using the LUR
parametrization, similar challenges are needed to be tackled to ensure
statistical properties hold.

####  Asymptotic Properties

Next, consider for a moment the limiting distribution of
$\left(\widehat{\delta}^{ols}-\delta\right)$ where
$\bar{y}_{2,t}=\delta\widehat{y}_{1,t}+\eta_{t}$. Therefore, the OLS estimator
of $\delta$ is given by
$\widehat{\delta}^{ols}=\displaystyle\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}\bar{y}_{2,t}\widehat{y}_{1,t}\right)$,
where $\widehat{y}_{1,t}:=\widehat{\boldsymbol{\beta}}_{1}^{\prime
ivx}\boldsymbol{x}_{1,t-1}$ which implies that it can be expressed as below

$\displaystyle\widehat{\delta}^{ols}=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}\bar{y}_{2,t}\widehat{\boldsymbol{\beta}_{1}^{\prime
ivx}\boldsymbol{x}_{1,t-1}}\right)=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}\bar{y}_{2,t}\boldsymbol{x}_{1,t-1}\right)\widehat{\boldsymbol{\beta}_{1}^{\prime
ivx}}.$ (3.18)

However,
$\bar{y}_{2,t}=\boldsymbol{\beta}_{2}^{\prime}\bar{\boldsymbol{x}}_{2,t-1}+\delta\widehat{y}_{1,t}+u_{2,t}$
which implies that

$\displaystyle\widehat{\delta}^{ols}$
$\displaystyle=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}\bigg{[}\boldsymbol{\beta}_{2}^{\prime}\bar{\boldsymbol{x}}_{2,t-1}+\delta\widehat{y}_{1,t}+u_{2,t}\bigg{]}\boldsymbol{x}_{1,t-1}^{\prime}\right)\widehat{\boldsymbol{\beta}_{1}^{ivx}}$
$\displaystyle=\left(\sum_{t=1}^{n}\widehat{y}_{1,t}\widehat{y}_{1,t}^{\prime}\right)^{-1}\times\left\\{\boldsymbol{\beta}_{2}^{\prime}\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\right)+\delta\left(\sum_{t=1}^{n}\widehat{y}_{1,t}\boldsymbol{x}_{1,t-1}^{\prime}\right)\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{x}_{1,t-1}^{\prime}\right)\right\\}\widehat{\boldsymbol{\beta}_{1}^{ivx}}$

where

$\displaystyle\left(\sum_{t=1}^{n}\widehat{y}_{1,t}\widehat{y}_{1,t}^{\prime}\right)^{-1}=\left(\sum_{t=1}^{n}\widehat{\boldsymbol{\beta}_{1}^{\prime
ivx}\boldsymbol{x}_{1,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\widehat{\boldsymbol{\beta}}_{1}^{ivx}}\right)^{-1}=\left(\widehat{\boldsymbol{\beta}}_{1}^{ivx}\right)^{-1}\left(\sum_{t=1}^{n}\boldsymbol{x}_{1,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\right)^{-1}\left(\widehat{\boldsymbol{\beta}}_{1}^{\prime
ivx}\right)^{-1}$

Therefore, we obtain that

$\displaystyle\widehat{\delta}^{ols}=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\boldsymbol{\beta}_{2}^{\prime}\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\right)\widehat{\boldsymbol{\beta}_{1}^{ivx}}+\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\delta\left(\sum_{t=1}^{n}\widehat{y}_{1,t}\boldsymbol{x}_{1,t-1}^{\prime}\right)\widehat{\boldsymbol{\beta}_{1}^{ivx}}+\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{x}_{1,t-1}^{\prime}\right)\widehat{\boldsymbol{\beta}_{1}^{ivx}}$

Notice that the second term above becomes:

$\displaystyle\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\delta\left(\sum_{t=1}^{n}\widehat{\boldsymbol{\beta}_{1}^{\prime
ivx}\boldsymbol{x}_{1,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\widehat{\boldsymbol{\beta}}_{1}^{ivx}}\right)=\delta\left(\sum_{t=1}^{n}\widehat{\boldsymbol{\beta}_{1}^{\prime
ivx}\boldsymbol{x}_{1,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\widehat{\boldsymbol{\beta}}_{1}^{ivx}}\right)^{-1}\left(\sum_{t=1}^{n}\widehat{\boldsymbol{\beta}_{1}^{\prime
ivx}\boldsymbol{x}_{1,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\widehat{\boldsymbol{\beta}}_{1}^{ivx}}\right)=\delta$

Thus, we have that

$\displaystyle\left(\widehat{\delta}^{ols}-\delta\right)=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\boldsymbol{\beta}_{2}^{\prime}\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\right)\widehat{\boldsymbol{\beta}}_{1}^{ivx}+\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{x}_{1,t-1}^{\prime}\right)\widehat{\boldsymbol{\beta}_{1}^{ivx}}.$

Furthermore, since the nonstationary regressors of node 1 are not correlated
with the nonstationary regressors of node 2, then the term
$\underset{n\to\infty}{\mathsf{plim}}\left(\sum_{t=1}^{n}u_{2,t}\boldsymbol{x}_{1,t-1}^{\prime}\right)=0$,
which implies that

$\displaystyle\left(\widehat{\delta}^{ols}-\delta\right)$
$\displaystyle=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\boldsymbol{\beta}_{2}^{\prime}\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\right)\times\widehat{\boldsymbol{\beta}_{1}^{ivx}}$
$\displaystyle=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)\times\left[\boldsymbol{\beta}_{1}+\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{1,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bar{u}_{1,t}\boldsymbol{z}_{1,t-1}^{\prime}\right)\right]$
$\displaystyle=\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)\boldsymbol{\beta}_{1}+\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)$
$\displaystyle\ \ \ \ \ \ \ \
\times\left[\left(\sum_{t=1}^{n}\bar{\boldsymbol{x}}_{1,t-1}\boldsymbol{z}_{1,t-1}^{\prime}\right)^{-1}\left(\sum_{t=1}^{n}\bar{u}_{1,t}\boldsymbol{z}_{1,t-1}^{\prime}\right)\right]+o_{p}(1).$

Since the term $\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)$ is always
positive then we assume that
$\underset{n\to\infty}{\mathsf{plim}}\left(\sum_{t=1}^{n}\widehat{y}^{2}_{1,t}\right)^{-1}=\mathcal{K}$
where $\mathcal{K}$ is some positive constant.

Therefore, putting all related expressions together we obtain

$\displaystyle
n^{\frac{1+\gamma_{z}}{2}}\bigg{[}\bigg{(}\tilde{\boldsymbol{\beta}}_{2}^{ivx}-\boldsymbol{\beta}_{2}\bigg{)}-\boldsymbol{\beta}_{1}\bigg{]}$
$\displaystyle=\bigg{\\{}\mathcal{K}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)\boldsymbol{\beta}_{1}+\mathcal{K}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)$
$\displaystyle\
\times\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}\bigg{\\}}$
$\displaystyle\
\times\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}$
$\displaystyle+\mathcal{MN}\bigg{(}0,\sigma_{u_{2}}^{2}\times\left(\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\bigg{)}$
$\displaystyle=\mathcal{K}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)\boldsymbol{\beta}_{1}\times\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}$
$\displaystyle+\mathcal{K}\boldsymbol{\beta}_{2}^{\prime}\left(\mathbb{E}\big{[}\boldsymbol{x}_{2,t-1}\boldsymbol{x}_{1,t-1}^{\prime}\big{]}\right)\times\bigg{\\{}\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}\bigg{\\}}^{2}$
$\displaystyle+\mathcal{MN}\bigg{(}0,\sigma_{u_{2}}^{2}\times\left(\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{2}u_{2}}^{-1}\bigg{)}.$

The second term above
$\bigg{\\{}\mathcal{MN}\bigg{(}0,\sigma_{u_{1}}^{2}\times\left(\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\right)^{\prime}\boldsymbol{C}_{z}\boldsymbol{V}_{C_{z}}\boldsymbol{C}_{z}\boldsymbol{\Psi}_{u_{1}u_{1}}^{-1}\bigg{)}\bigg{\\}}^{2}$,
corresponds to a squared mixed Gaussian distribution which gives a form of a
$\chi^{2}$ distribution. Thus, the limiting distribution of (3.2.2) is a
linear combination of a $\chi^{2}$ distribution and a mixed Gaussian
distribution, which is a Generalized $\chi^{2}$ distribution. Thus the $y_{t}$
inherits the properties of $x_{t}$, through a cointegrating relation,
especially our goal is to have a stationary innovation sequence such that
$u_{t}\sim I(0)$. Under the assumption of a consistent estimator for
$\hat{\delta}^{ols}$ then we have convergence in probability to zero such that
$\left(\hat{\delta}^{ols}-\delta\right)\overset{p}{\to}0$. Thus, the first
term of the expression asymptotically tends to zero (negligible) which implies
that the second term converges into a mixed Gaussian distribution.

##  Doubly Corrected Estimation Method: Conditional Quantile Case

In this section we consider the double corrected conditional quantile
estimation methodology which is the main focus of our research study. The
derivations presented in the previous sections were useful to shed light on
some challenges we have to overcome to develop a robust framework in the
proposed setting of quantile predictive regression models. In terms of
estimation approach we focus on the QR estimators under nonstationarity which
are adjusted using the IVX filtration. During the first stage of the
estimation procedure in the same spirit as several studies in the literature
related to the joint estimation of the risk measures of
$(\mathsf{VaR},\mathsf{CoVaR})$ (see, Adrian and Brunnermeier, (2016), Härdle
et al., (2016) and Patton et al., (2019) among others), the $\mathsf{VaR}$ is
estimated given a fixed quantile $\tau\in(0,1)$ level but using information on
the nonstationary properties of regressors as given by the quantile predictive
regression system below

  * Stage 1:

$\displaystyle y_{t}^{(1)}$
$\displaystyle=\alpha_{1}(\tau)+\boldsymbol{\beta}_{1}^{\prime}(\tau)\boldsymbol{x}_{t-1}^{(1)}+u_{t}^{(1)}(\tau),\
\ \ \text{for}\ \ \ t=1,...,n$ (4.1) $\displaystyle\boldsymbol{x}_{t}^{(1)}$
$\displaystyle=\boldsymbol{R}_{n}^{(1)}\boldsymbol{x}^{(1)}_{t-1}+\boldsymbol{v}^{(1)}_{t}$
(4.2)

such that
$\alpha_{1}(\tau)+\boldsymbol{\beta}^{\prime}(\tau)\boldsymbol{x}_{t-1}^{(1)}$
is the $\tau-$conditional quantile of $y_{t}$ given $\boldsymbol{x}_{t-1}$,
where the unknown parameters $(\alpha_{1}(\tau),\boldsymbol{\beta}_{1}(\tau))$
are estimated using the QR objective function.

  * Stage 2:

$\displaystyle y_{t}^{(2)}$
$\displaystyle=\alpha_{2}(\tau)+\boldsymbol{\beta}_{2}^{\prime}(\tau)\boldsymbol{x}_{t-1}^{(2)}+\delta(\tau)\widehat{y}_{1,t}\left(\hat{\boldsymbol{\theta}}^{ivx}_{1}(\uptau)\right)+u_{t}^{(2)}(\tau),\
\ \ \text{for}\ \ \ t=1,...,n$ (4.3) $\displaystyle\boldsymbol{x}_{t}^{(2)}$
$\displaystyle=\boldsymbol{R}_{n}^{(2)}\boldsymbol{x}^{(2)}_{t-1}+\boldsymbol{v}^{(2)}_{t}$
(4.4)

where the parameter of interest corresponds to the estimator
$\boldsymbol{\vartheta}^{ivx}_{n}=\big{(}\alpha_{2}(\tau),\boldsymbol{\beta}_{2}^{\prime}(\tau),\delta(\tau)\big{)}$.

######  Remark 5.

Similar to the remark found in Wang and Zhao, (2016), we point out that it is
up to the practitioner to determine the specific model and parameter
estimation method, which are the starting point to carry out any subsequent
$\mathsf{CoVaR}$ estimation and inference. In our study, we use the
nonstationary quantile predictive regression when obtaining an estimate for
$\mathsf{CoVaR}$. Thus, our approach gives a $\sqrt{n^{1+\delta}}-$consistent
estimator for the unknown coefficient of systemic risk based on the IVX
filter. Recall that
$\widehat{\boldsymbol{\beta}}^{\star}(\uptau):=\underset{\boldsymbol{\beta}^{\star}\in\mathbb{R}^{p+1}}{\mathsf{arg\
min}}\sum_{i=1}^{n}\rho_{\uptau}\bigg{(}y_{t}-\tilde{\boldsymbol{x}}_{t-1}^{\prime}\boldsymbol{\beta}^{\star}\bigg{)}$
where the set of regressors
$\tilde{\boldsymbol{x}}_{t-1}=\big{[}\boldsymbol{x}_{t-1}\ \ \
\widehat{x}^{\mathsf{VaR}}_{t}\big{]}$. The particular optimization function
corresponds to the second-stage procedure which gives the model estimates to
construct the one-period ahead forecasted $\mathsf{CoVaR}$ risk measure.
However, our approach differs from the estimation method of Adrian and
Brunnermeier, (2016) and Härdle et al., (2016) since nonstationary time series
data are modelled via the LUR parametrization and thus we employ the QR-IVX
estimator proposed by Lee, (2016) but with the additional generated regressor
from the first-stage procedure.

Next, we focus on the required preliminary theory for establishing the
consistency and asymptotic mixed Gaussianity of the doubly IVX corrected
estimator $\boldsymbol{\beta}^{ivx}(\uptau)$. Notice that in contrast to the
approach proposed by Bhattacharya, (2020) our generated variable corresponds
to the estimated $\mathsf{VaR}$ under time series nonstationarity which is
used as an additional regressor in a nonstationary quantile predictive
regression model rather a generated dependent variable. In practice, we do
obtain a generated dependent random variable which is used as an additional
covariate in a quantile predictive regression model with persistence
regressors. On the other hand, in a sense our framework is closer to the
approach of Chen et al., (2021) although we have a different model. In
addition, to ensure robustness to the abstract degree of persistence we use
the corresponding IVX estimator rather than the OLS approach thereby avoiding
the necessity to use the control function approach as a method to correct for
the presence of the nuisance parameter of persistence.

###  On Generated (quantile) regressor in stationarity time series models

Denote the conditional CDF of $Y$ given $X=x$ by $F_{Y|X}(.|x)$ and its
conditional quantile at $\uptau\in(0,1)$ by $Q(\uptau|x)$, such that
$Q(\uptau|x)=F_{Y|X}^{-1}(\uptau|x)=\mathsf{inf}\big{\\{}s:F_{Y|X}(s|x)\geq\uptau\big{\\}}$.
Suppose that the quantile $Q(\uptau|x)$ is modelled as a general nonlinear
function of $x$ and $\uptau$. We fix $x$ and treat $Q(\uptau|x)$ as a
stochastic process in $\uptau$, where
$\uptau\in\mathcal{T}=[\lambda_{1},\lambda_{2}]$ with
$0<\lambda_{1}\leq\lambda_{2}<1$.

A key quantity for establishing asymptotic theory and inferential methods is
the regression-quantile process such that
$\boldsymbol{M}_{n}(\uptau)\equiv\left\\{\sqrt{n}\left(\widehat{\boldsymbol{\beta}}_{n}(\uptau)-\boldsymbol{\beta}(\uptau)\right)\right\\}$
for some $\uptau\in(0,1)$ which is regarded as a stochastic process in the
space $\big{(}D[0,1]\big{)}^{d}$ of $\mathbb{R}^{d}-$valued right-continuous
functions with left-hand limits on $[0,1]$ (see, Goh and Knight, (2009), Kato,
(2009) and Portnoy, (2012) among others101010To establish the limiting
behaviour of such stochastic processes several studies in the literature
employ the Bahadur representation which holds uniformly in $\theta$ and
$\uptau\in(0,1)$. Moreover, estimating the parameter $\theta$ induces some
important changes compared to a known $\theta_{0}$ and required finding an
expansion for
$\sqrt{n}\left(\widehat{\boldsymbol{\beta}}(\uptau;\widehat{\theta})-\boldsymbol{\beta}(\uptau;\widehat{\theta})\right)$.).

######  Definition 1.

(Weak Uniform Convergence on $\Theta$) Let $(\Omega,\mathcal{F},\mathcal{P})$
be a suitable probability space. Then, $Q_{n}(.,\theta)-\bar{Q}_{n}(\theta)\to
0$ as $n\to\infty$ in probability law uniformly on $\Theta$ if for any
$\epsilon>0$ and any $\delta>0$ there exists an integer
$N(\epsilon,\delta)<\infty$ such that for all $n>N(\epsilon,\delta)$ then

$\displaystyle\mathbb{P}\left[\underset{\Theta}{\mathsf{sup}}\left|Q_{n}(.,\theta)-\bar{Q}_{n}(\theta)\right|<\epsilon\right]>1-\delta.$
(4.5)

######  Remark 6.

If $Q_{n}$ and $\bar{Q}_{n}$ are vector or matrix valued functions then we use
identical terminology to mean that uniform convergence holds element by
element. Moreover, the assumption of compactness of $\Theta$ is crucial since
it ensures that the supremum of $\Theta$ exists and that the functions
$\underset{\Theta}{\mathsf{sup}}\left|Q_{n}(.,\theta)-\bar{Q}_{n}(\theta)\right|$
are measurable (see, Saikkonen, (1995)). One then would need to consider the
asymptotic behaviour of the following functional (see, Bhattacharya, (2020))

$\displaystyle\widehat{\boldsymbol{\mathcal{E}}}\left(\uptau;\theta\right)=\sqrt{n}\left(\widehat{\boldsymbol{\beta}}\left(\uptau;\theta\right)-\boldsymbol{\beta}\left(\uptau;\theta\right)\right)-\bigg{(}\boldsymbol{H}^{-1}\left(\uptau;\theta\right)\widehat{\boldsymbol{S}}\left(\uptau;\theta\right)\bigg{)}.$
(4.6)

A uniform order for
$\widehat{\boldsymbol{\mathcal{E}}}\left(\uptau;\theta\right)$ relies on a
uniform order study for the remainder term which crucially depends on the
properties of the quantile regression of interest. To do this, one can employ
concepts of maximal inequality under bracketing conditions as given in the
related literature as well as other linearization techniques (e.g., see Goh
and Knight, (2009), Kato, (2009) and Brian, (2021)).

######  Example 1.

As an illustrative example, consider the econometric model studied by
Bhattacharya, (2020). In particular, it holds that for any $\uptau\in(0,1)$

$\displaystyle\sqrt{n}\left(\widehat{\beta}(\uptau)-\beta(\uptau)\right)\overset{d}{\to}\mathcal{N}\bigg{(}0,V(\uptau)\bigg{)}.$
(4.7)

Following the derivations in Proposition 1 of Bhattacharya, (2020), it yields
that

$\displaystyle\sqrt{n}\left(\widehat{\beta}\big{(}\uptau;\widehat{\theta}\big{)}-\beta(\uptau;\theta_{0})\right)$
$\displaystyle=\sqrt{n}\bigg{\\{}\widehat{\beta}(\uptau;\widehat{\theta})-\beta(\uptau;\widehat{\theta})\bigg{\\}}+\sqrt{n}\bigg{\\{}\beta(\uptau;\widehat{\theta})-\beta(\uptau;\theta_{0})\bigg{\\}}$
$\displaystyle=\sqrt{n}\bigg{\\{}\widehat{\beta}(\uptau;\widehat{\theta})-\beta(\uptau;\widehat{\theta})\bigg{\\}}+\left(\frac{\partial\beta(\uptau;\theta_{0})}{\partial\theta}+\mathcal{O}_{p}(1)\right)\sqrt{n}\left(\widehat{\theta}-\theta_{0}\right)$
$\displaystyle=\sqrt{n}\bigg{\\{}\widehat{\beta}(\uptau;\widehat{\theta})-\beta(\uptau;\widehat{\theta})\bigg{\\}}\left(\frac{\partial\beta(\uptau;\theta_{0})}{\partial\theta}\right)^{\prime}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\Psi(Z_{i})+\mathcal{O}_{p}(1)$

Therefore, it can be proved that

$\displaystyle\sqrt{n}\left(\widehat{\beta}\big{(}\uptau;\widehat{\theta}\big{)}-\beta(\uptau;\theta_{0})\right)=H^{-1}\left(\uptau;\widehat{\theta}\right)\widehat{S}\left(\uptau;\widehat{\theta}\right)+\left(\frac{\partial\beta(\uptau;\theta_{0})}{\partial\theta}\right)^{\prime}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\Psi(Z_{i})+\mathcal{O}_{p}(1).$
(4.8)

######  Remark 7.

Suitable Bahadur-based representations for establishing asymptotic theory
results are considered in several studies (see, Portnoy, (2012) and Wang and
Zhao, (2016)). However, within our setting we have the additional complexity
that the quantile-based estimators in each stage has to be robustify against
the unknown persistence properties that appear due to the LUR parametrization
of regressors (e.g., see Lee, (2016)). The LUR parametrization includes the
nuisance parameter of persistence, so we expect that conventional estimators,
such as the ordinary least squares estimator, to asymptotically depend on the
coefficient of persistence. Thus, we employ the IVX filtration of Phillips and
Magdalinos, (2009) which is found to robustify estimation and inference in
predictive regression models (see, also Kostakis et al., (2015)).

######  Example 2.

In this example, we discuss the econometric framework of Chen et al., (2021),
who consider quantile regression with generated regressors. Denote with
$F_{t}\left(u|\boldsymbol{x}_{t}\right)$, the conditional distribution
function of the error term with continuous densities
$f_{t}\left(u|\boldsymbol{x}_{t}\right)$ and a unique conditional $\uptau-$th
quantile which are uniformly bounded away from 0 and $\infty$. It holds that

$\displaystyle\sqrt{n}\left(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}\right)=n^{-1/2}\sum_{t=1}^{n}\boldsymbol{r}_{t}(\boldsymbol{\theta})+o_{p}(1)$
(4.9)

where $\boldsymbol{r}_{t}(.)$ is a continuous function which satisfies
$\mathbb{E}\left[\boldsymbol{r}_{t}(\boldsymbol{\theta})\right]=0$ and
Var$\left[\boldsymbol{r}_{t}(\boldsymbol{\theta})\right]=\boldsymbol{V}$ and
thus the generated regressor has the form
$\frac{1}{n}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}\hat{\boldsymbol{x}}_{t}^{\prime}-\frac{1}{n}\sum_{t=1}^{n}\boldsymbol{x}_{t}\boldsymbol{x}_{t}^{\prime}=o_{p}(1)$.

Moreover, for any fixed $\uptau\in(0,1)$, we have the following results

$\displaystyle\sqrt{n}\left(\hat{\boldsymbol{\beta}}(\uptau)-\boldsymbol{\beta}(\uptau)\right)\overset{d}{\to}\mathcal{N}\big{(}\boldsymbol{0},\boldsymbol{\Omega}(\uptau)\big{)}.$
(4.10)

Under the null hypothesis, it holds that

$\displaystyle\sqrt{n}\left(\boldsymbol{R}\hat{\boldsymbol{\beta}}(\uptau)-\boldsymbol{r}\right)\overset{d}{\to}\mathcal{N}\big{(}\boldsymbol{0},\boldsymbol{R}\boldsymbol{\Omega}(\uptau)\boldsymbol{R}^{\prime}\big{)}.$
(4.11)

Then, since we proved that $\hat{\boldsymbol{\Omega}}(\uptau)$ is a consistent
estimator of $\boldsymbol{\Omega}(\uptau)$, by Slutsky’s theorem,

$\displaystyle\mathcal{W}_{T}=T\left(\boldsymbol{R}\widehat{\boldsymbol{\beta}}(\uptau)-\boldsymbol{r}\right)^{\prime}\big{[}\boldsymbol{R}\widehat{\boldsymbol{\Omega}}(\uptau)\boldsymbol{R}^{\prime}\big{]}^{-1}\left(\boldsymbol{R}\widehat{\boldsymbol{\beta}}(\uptau)-\boldsymbol{r}\right)\sim\chi^{2}_{q}.$
(4.12)

Using the mean value expansion theorem we have that

$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}u_{t}=\frac{1}{n}\sum_{t=1}^{n}\boldsymbol{x}_{t}u_{t}+\left\\{\frac{1}{n}\sum_{t=1}^{n}\nabla_{\delta}h\left(\boldsymbol{w}_{t},\boldsymbol{\delta}\right)^{\prime}u_{t}\right\\}\big{(}\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}\big{)}+o_{p}(1).$
(4.13)

Since we have that $\mathbb{E}\left(u|\boldsymbol{w}\right)=0$ and
$\mathbb{E}\left[\nabla_{\delta}h\left(\boldsymbol{w}_{t},\boldsymbol{\delta}\right)^{\prime}u_{t}\right]=o_{p}(1)$.
Furthermore, since
$\big{(}\hat{\boldsymbol{\delta}}-\boldsymbol{\delta}\big{)}=o_{p}(1)$ and
$\mathbb{E}\left(\boldsymbol{x}_{t}u_{t}\right)=0$, it follows that

$\displaystyle\frac{1}{n}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}u_{t}=\frac{1}{n}\sum_{t=1}^{n}\boldsymbol{x}_{t}u_{t}+o_{p}(1)=o_{p}(1).$
(4.14)

Therefore, it can be proved that
$\big{(}\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0}\big{)}=o_{p}(1)$ and
the asymptotic normality of the OLS-GR estimator can be also established. In
particular, we have that

$\displaystyle\sqrt{n}\big{(}\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}_{0}\big{)}$
$\displaystyle=\left(\frac{1}{n}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}\hat{\boldsymbol{x}}_{t}^{\prime}\right)^{-1}\left\\{\frac{1}{\sqrt{n}}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}\left[\left(\boldsymbol{x}_{t}-\hat{\boldsymbol{x}}_{t}\right)^{\prime}\boldsymbol{\beta}_{0}+u_{t}\right]\right\\},$
$\displaystyle=\widehat{\boldsymbol{H}}^{-1}\left\\{\frac{1}{\sqrt{n}}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}\left[\left(\boldsymbol{x}_{t}-\hat{\boldsymbol{x}}_{t}\right)^{\prime}\boldsymbol{\beta}_{0}+u_{t}\right]\right\\}$

where
$\widehat{\boldsymbol{H}}=\left(\frac{1}{n}\sum_{t=1}^{n}\hat{\boldsymbol{x}}_{t}\hat{\boldsymbol{x}}_{t}^{\prime}\right)$.
Moreover, to establish the consistency of the covariance matrix
$\boldsymbol{\Omega}(\uptau)$ we assume that
$f_{t}\left(\lambda|\boldsymbol{x}_{t}\right)$ satisfies the Lipschitz
condition that
$\left|f_{t}\left(\lambda_{1}|\boldsymbol{x}_{t}\right)-f_{t}\left(\lambda_{2}|\boldsymbol{x}_{t}\right)\right|\leq
L_{t}\left|\lambda_{1}-\lambda_{2}\right|$ for some constant $L_{t}<\infty$
for all $t$.

######  Remark 8.

Notice that each of the GRs can be related to different observable variables
in different functional forms. Therefore, for simplicity we assume that the GR
equations can be estimated separately, and thus any dependence between the GRs
can be captured by the variance covariance matrix. However notice that within
our setting it has been proved by Adrian and Brunnermeier, (2016), in the
stationary case, that the quantile regression with generated regressors that
the standard quantile estimation based on the OLS method is still consistent.
Therefore, our focus is on the implementation of the IVX estimator in the
nonstationary case, that is the case in which we assume the presence of the
LUR specification for the lagged regressors of the model.

###  On Generated (quantile) regressor under nonstationarity

In relation to technical conditions, the notion of stochastic equicontinuity
can be employed in nonstationary time series environments. These conditions
are not uncommon in the cointegrating regression literature such as the
restricted cointegrated regression model of Saikkonen, (1995) and the
constrained least squares estimator proposed by Nagaraj and Fuller, (1991)
(see, Moon and Schorfheide, (2002)). Furthermore, as pointed out by Saikkonen,
(1995), the classical stochastic equicontinuity conditions can be too strong
in cointegrated regression models. In particular, due to the stationarity of
the generated regressor we shall assume that the parameter space in the
nonstationary quantile predictive regression model in the second-stage
procedure corresponds to a short-run parameter while the coefficient of the
nonstationary regressors corresponds to the long-run parameters of the model.
This can be done by partitioning $\theta=(\theta_{1},\theta_{2})$ and
considering a standardization matrix such that

$\displaystyle
N_{n}(\theta_{1},\delta)=\big{\\{}\boldsymbol{\theta}_{1}^{*}\in\Theta_{1}:\left\lVert\boldsymbol{D}_{1n}\left(\boldsymbol{\theta}_{1}^{*}-\boldsymbol{\theta}_{1}\right)\right\rVert<\delta\big{\\}}$
(4.15)

where $\boldsymbol{D}_{1n}$ is $(k_{1}\times k_{1})$ nonstochastic diagonal
matrix whose diagonal elements are positive and increasing functions of $n$.
Moreover, these SE conditions are considered to be sufficient conditions for
establishing the weak convergence of such estimators while necessary
conditions imply the inverse information matrix has eigenvalues bounded away
from zero in probability. Therefore, in the case of integrated and
cointegrated variables, Saikkonen, (1995) show that the standardized Hessian
matrix
$\displaystyle\boldsymbol{D}_{\boldsymbol{\theta}n}^{-1}\frac{\partial^{2}\bar{Q}_{n}(\boldsymbol{\theta}_{0})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\boldsymbol{D}_{\boldsymbol{\theta}n}^{-1}$
converges weakly and due to the strict exogeneity assumption, the weak limit
is block diagonal between the parameters $\boldsymbol{\theta}_{1}$, $\psi$ and
$\alpha$. The next step would be to show that the standardized Hessian
satisfies an appropriate stochastic equicontinuity condition, that along with
a suitable consistency result, makes it possible to obtain the limiting
distribution.

######  Example 3 (see, Cho et al., (2015)).

Consider the ARDL process denoted by

$\displaystyle
Y_{t}=\alpha+\sum_{j=1}^{p}\phi_{j}Y_{t-j}+\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}\boldsymbol{X}_{t-j}+U_{t},$

where $\boldsymbol{X}_{t}\in\mathbb{R}^{k}$ is a $k-$dimensional integrated
process of stationary and ergodic processes. Moreover, let the $\tau-$quantile
of $Y_{t}$ conditional on the $\sigma-$algebra $\mathcal{F}_{t-1}$ denoted by
$Q_{y_{t}}(\tau|\mathcal{F}_{t-1})$ be

$\displaystyle
Q_{y_{t}}(\tau|\mathcal{F}_{t-1})=\alpha(\tau)+\sum_{j=1}^{p}\phi_{j}(\tau)Y_{t-j}+\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}(\tau)\boldsymbol{X}_{t-j}$

We consider that $Y_{t}$ is represented as a quantile autoregressive
distributed lag (QARDL) process

$\displaystyle
Q_{y_{t}}(\tau|\mathcal{F}_{t-1})=\alpha(\tau)+\sum_{j=1}^{p}\phi_{j}(\tau)Y_{t-j}+\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}(\tau)\boldsymbol{X}_{t-j}+U_{t}(\tau)$
(4.16)

where $U_{t}(\tau)=Y_{t}-Q_{y_{t}}(\tau|\mathcal{F}_{t-1})$. We reformulate
the above expression so that it can capture the long-run relationship between
$Y_{t}$ and $\boldsymbol{X}_{t}$ using the following long-run quantile process

$\displaystyle
Y_{t}=\mu(\tau)+\boldsymbol{X}_{t}^{\prime}\boldsymbol{\beta}(\tau)+\mathcal{R}_{t}(\tau)$
(4.17)

We explain below the steps and the coefficients to obtain the above
reformulation.

  * Step 1.

We decompose the term
${\color[rgb]{0,0,1}\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}(\tau)\boldsymbol{X}_{t-j}}\equiv{\color[rgb]{0,0,1}\boldsymbol{X}_{t}^{\prime}\boldsymbol{\gamma}(\tau)+\sum_{j=0}^{q-1}\boldsymbol{Z}_{t-j}^{\prime}\boldsymbol{\delta}(\uptau)}$
which implies that

$\displaystyle
Y_{t}=\alpha(\tau)+\sum_{j=1}^{p}\phi_{j}(\tau)Y_{t-j}+\underbrace{{\color[rgb]{0,0,1}\boldsymbol{X}_{t}^{\prime}\boldsymbol{\gamma}(\tau)+\sum_{j=0}^{q-1}\boldsymbol{Z}_{t-j}^{\prime}\boldsymbol{\delta}(\uptau)}}_{\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}(\tau)\boldsymbol{X}_{t-j}}+U_{t}(\tau)$

where
$\boldsymbol{\gamma}(\tau):=\sum_{j=0}^{q}\boldsymbol{\theta}_{j}(\tau)$,
$\boldsymbol{\delta}(\tau):=-\sum_{i=j+1}^{q}\boldsymbol{\theta}_{i}(\tau)$
and $\boldsymbol{Z}_{t}=\Delta\boldsymbol{X}_{t}$.

  * Step 2.

Suppose that

$\displaystyle
Y_{t}=\alpha(\tau)+\sum_{j=1}^{p}\phi_{j}(\tau)Y_{t-j}+\underbrace{{\color[rgb]{0,0,1}\boldsymbol{X}_{t}^{\prime}\boldsymbol{\gamma}(\tau)+\sum_{j=0}^{q-1}\boldsymbol{Z}_{t-j}^{\prime}\boldsymbol{\delta}(\uptau)}}_{\sum_{j=0}^{q}\boldsymbol{\theta}_{j}^{\prime}(\tau)\boldsymbol{X}_{t-j}}+U_{t}(\tau)\equiv\mu(\tau)+\boldsymbol{X}_{t}^{\prime}\boldsymbol{\beta}(\tau)+\mathcal{R}_{t}(\tau)$

where it holds that

$\displaystyle\boldsymbol{\beta}(\tau):=\frac{\boldsymbol{\gamma}(\tau)}{\displaystyle\left(1-\sum_{j=1}^{p}\phi_{j}(\tau)\right)},\
\
\mu(\tau):=\frac{\alpha(\tau)}{\displaystyle\left(1-\sum_{j=1}^{p}\phi_{j}(\tau)\right)}$
(4.18)

and the remainder term can be written such that

$\displaystyle
R_{t}(\tau):=\sum_{j=0}^{\infty}\boldsymbol{Z}^{\prime}_{t-j}\boldsymbol{\xi}_{j}(\tau)+\sum_{j=0}^{\infty}\eta_{j}(\tau)U_{t-j}(\tau)$

where each $\boldsymbol{\xi}_{j}(\tau)$ corresponds to the "upper" unbounded
sum of some positive sequence of coefficients $\boldsymbol{\pi}_{\ell}(\tau)$
such that
$\boldsymbol{\xi}_{j}(\tau)\equiv\sum_{\ell=j+1}^{\infty}\boldsymbol{\pi}_{\ell}(\tau)$
where the set of coefficients
$\big{\\{}\boldsymbol{\pi}_{0}(\tau),\boldsymbol{\pi}_{1}(\tau),...\big{\\}}$
and $\big{\\{}\eta_{0}(\tau),\eta_{1}(\tau),...\big{\\}}$ which implies that

$\displaystyle\sum_{j=0}^{\infty}\eta_{j}(\tau)L^{j}$
$\displaystyle\equiv\left(1-\sum_{j=1}^{p}\phi_{j}(\tau)L^{j}\right)^{-1},$
$\displaystyle\sum_{j=0}^{\infty}\boldsymbol{\xi}_{j}(\tau)L^{j}$
$\displaystyle\equiv(1-L)^{-1}\left(\frac{\displaystyle\sum_{j=0}^{q}\boldsymbol{\theta}_{j}(\tau)L^{j}}{\displaystyle
1-\sum_{j=1}^{p}\phi_{j}(\tau)L^{j}}-\frac{\displaystyle\sum_{j=0}^{q}\boldsymbol{\theta}_{j}(\tau)}{\displaystyle
1-\sum_{j=1}^{p}\phi_{j}(\tau)}\right).$

Our main interest lie in developing the estimation theory for the long-run
parameter $\boldsymbol{\beta}(\tau)$ which is a function of
$\boldsymbol{\gamma}(\tau)$ and
$\boldsymbol{\phi}_{j}:=\big{(}\phi_{1}(\tau),...,\phi_{p}(\tau)\big{)}^{\prime}$.

######  Remark 9.

Recall that a key condition for estimating the quantile-dependent
cointegrating vector is that there exists a zero correlation between
regressors and errors (see, Xiao, (2009)). Moreover, cases where estimators
are derived under the assumption large $N,T$ asymptotics include the
frameworks of Zhu et al., (2019), Chen, (2022) and Feng et al., (2024) among
others. Nevertheless, the above derivations provide a good example on the
possible system representation we can employ in the case we consider the
generated regressors that correspond to the estimated $\mathsf{VaR}$ (under
time series nonstationarity) across the cross-section. The above example which
corresponds to the ARDL process as in the study of Cho et al., (2015) is
presented also by Katsouris, 2023b .

######  Proposition 1.

Suppose that the parameter vector of interest is given by

$\displaystyle\boldsymbol{\vartheta}(\uptau):=\big{(}\hat{\alpha}(\uptau),\hat{\boldsymbol{\beta}}^{\prime}_{x}(\uptau),\hat{\boldsymbol{\beta}}^{\prime}_{\hat{y}}(\uptau)\big{)}^{\prime}\equiv\underset{\boldsymbol{\vartheta}\in\Theta}{\mathsf{arg\
min}}\
\sum_{t=1}^{n}\rho_{\uptau}\big{(}y_{t}-\boldsymbol{x}_{t}^{\prime}\boldsymbol{\vartheta}\big{)}$
(4.19)

Moreover, if $\boldsymbol{D}_{0}$ exists and is positive definite, then as
$N,T\to\infty$ for any quantile level $\uptau\in(0,1)$, the following limit
results holds

$\displaystyle\sqrt{N}\left(\hat{\boldsymbol{\phi}}(\cdot)-\boldsymbol{\phi}_{0}(\cdot)\right)\overset{d}{\to}\boldsymbol{D}_{1NT}(\cdot)^{-1}\mathbb{G}_{N}\big{(}\psi(\theta_{i}-\boldsymbol{x}_{i}^{\prime}\boldsymbol{\phi}_{0}(\cdot))\big{)}+o_{p}(1)\Rightarrow\mathbb{G}_{N}(\cdot)\in\mathcal{\ell}^{\infty}(\mathcal{F}),$
(4.20)

where $\mathbb{G}(\cdot)$ is a zero mean Gaussian process with covariance
kernel
$\mathbb{E}\big{[}\mathbb{G}(\uptau)\mathbb{G}(\uptau)^{\prime}\big{]}$.

Moreover, we shall consider the weak convergence of the quantile regression
process $\hat{\boldsymbol{\phi}}(\cdot)$, where for any fixed quantile level
$\uptau\in(0,1)$ it holds that as $(N,T)\to\infty$,

$\displaystyle\sqrt{N}\left(\hat{\boldsymbol{\phi}}(\uptau)-\boldsymbol{\phi}_{0}(\uptau)\right)\overset{d}{\to}\mathcal{N}\left(0,\uptau(1-\uptau)\boldsymbol{\Xi}^{-2}\boldsymbol{D}_{0}^{-1}\right).$

Define with $\mathbb{G}(\uptau)$ is a zero mean Gaussian process with
covariance kernel

$\displaystyle\mathbb{E}\big{[}\mathbb{G}(\uptau)\mathbb{G}^{\prime}(\uptau^{\prime})\big{]}=\boldsymbol{\Xi}(\uptau)^{-1}(\uptau)S(\uptau,\uptau^{\prime})\boldsymbol{\Xi}(\uptau)^{-1}(\uptau)\boldsymbol{D}_{0}^{-1}.$
(4.21)

####  Misspecification Analysis of VaR-CoVaR risk matrix

In order to obtain some further insights in the implications of the
theoretical aspects discussed in this paper, we briefly discuss the extension
of the framework proposed by Katsouris, 2023c  into the econometric
environment of this study. Suppose that $\boldsymbol{y}_{it}$ is generated by
the true adjacency matrix $A$ (for simplicity we think of the novel
$\mathsf{VaR}-\Delta\mathsf{CoVaR}$ risk matrix proposed by Katsouris, (2021)
as an adjacency matrix to simplify the notation). Related studies from the
literature of network regression estimation and inference include Zhu et al.,
(2017), Zhu et al., (2019) and Zhu and Pan, (2020) among others).

Roughly speaking, the consistency result for the population counterpart of the
risk matrix might not be held when the adjacency matrix $A$ is misspecified to
be $A^{\star}=\left(a_{ij}^{\star}\right)$. Denote with $W^{\star}$ to be the
row-normalized $A^{\star}$ and $X_{it}^{\star}$ to be the matrix of model
regressors. Then, the estimator is

$\displaystyle\hat{\boldsymbol{\beta}}^{\star}(\uptau)=\underset{\boldsymbol{\beta}^{\star}}{\mathsf{arg\
min}}\sum_{i=1}^{N}\sum_{t=1}^{T}\uprho_{\uptau}\left(y_{it}-\boldsymbol{X}_{(i)t-1}^{\star\top}\boldsymbol{\beta}^{\star}(\uptau)\right).$
(4.22)

Define with

$\displaystyle\widehat{\boldsymbol{D}}_{0}$
$\displaystyle=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\boldsymbol{X}_{(i)t-1}^{\star}\boldsymbol{X}_{(i)t-1}^{\star\top},$
(4.23) $\displaystyle\widehat{\boldsymbol{D}}_{1}(\uptau)$
$\displaystyle=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}f_{it}\left(\boldsymbol{X}_{(i)t-1}^{\star\top}\boldsymbol{\beta}^{\star}(\uptau)\right)\boldsymbol{X}_{(i)t-1}^{\star}\boldsymbol{X}_{(i)t-1}^{\star\top}.$
(4.24)

######  Remark 10.

Notice that for the VaR-CoVaR risk matrix, misspecification has a different
interpretation than for the usual binary or weighted adjacency matrix for a
graph. The reason is that our risk matrix is constructed based on tail
forecasts using a regression-based methodology that incorporates the
persistence properties of the regressors (see, Katsouris, (2021); Katsouris,
2023c ).

The conditions associated with the misspecified coefficients are listed below:

######  Assumption 6.

(Eigenvalue-Bound)

Let
$\widehat{\boldsymbol{D}}_{1}(\uptau)\overset{p}{\to}\boldsymbol{D}_{1}(\uptau)$
as $\mathsf{min}\left\\{N,T\right\\}\to\infty$ for any $\uptau\in(0,1)$, where
$\boldsymbol{D}_{1}(\uptau)\in\mathbb{R}^{N\times N}$ is a positive definite
matrix. Then, there exists positive constants, $0<c_{1}<c_{2}<\infty$ such
that
$c_{1}\leq\lambda_{\mathsf{min}}\left(\boldsymbol{D}_{1}(\uptau)\right)\leq\lambda_{\mathsf{max}}\left(\boldsymbol{D}_{1}(\uptau)\right)\leq
c_{2}$ for any $\uptau\in(0,1)$.

######  Assumption 7.

(Monotonicity)

It is assumed that
$\boldsymbol{X}_{(i)t-1}^{\star\top}\boldsymbol{\beta}^{\star}(\uptau)$ such
that $1\leq i\leq N$ and $1\leq t\leq T$ is a monotone increasing function
with respect to $\uptau\in(0,1)$.

######  Corollary 1.

Define with
$d\left(W^{\star},W\right)=\displaystyle\sum_{i,t}\left|w_{ij}^{\star}-w_{ij}\right|$
to be the total amount of misspecification of $W$. Then, we have that

$\displaystyle\left(\hat{\boldsymbol{\beta}}^{\star}(\uptau)-\boldsymbol{\beta}^{\star}(\uptau)\right)=-\frac{1}{NT}\left\\{\boldsymbol{D}_{1}(\uptau)\right\\}^{-1}\sum_{i=1}^{N}\sum_{t=1}^{T}\boldsymbol{X}_{(i)t-1}^{\star}\psi_{\uptau}\left(\boldsymbol{\mathcal{E}}_{it\uptau}\right)+\mathcal{R}_{NT}\left(\uptau\right),$
(4.25)

with $\mathsf{sup}_{\uptau}\left\lVert
R_{NT}(\uptau)\right\rVert=o_{p}\left(\left(NT\right)^{-1/2}\right)$.

Moreover, assuming that the amount of misspecification has a bounded
convergence rate such that
$d\left(W^{\star},W\right)=o_{p}\left(\sqrt{N/T}\right)$, then we have that
$\left(\hat{\boldsymbol{\beta}}^{\star}(\uptau)-\boldsymbol{\beta}^{\star}(\uptau)\right)=o_{p}\left(\sqrt{N/T}\right)$.
Consequently, it can be shown that when the misspecification amount is under
control, that is,
$d\left(W^{\star},W\right)=\mathcal{O}_{p}\left(\sqrt{N/T}\right)$, the
resulting estimator is still $\sqrt{NT}-$consistent.

###### Proof.

We divide the proof into two parts. In the first part of the proof, we show
the representation given by expression (4.25) while in the second part we show
the consistency result for $\hat{\boldsymbol{\beta}}^{\star}(\uptau)$. Denote
with
$\Delta^{\star}\left(\uptau\right)=\hat{\boldsymbol{\beta}}^{\star}(\uptau)-\boldsymbol{\beta}^{\star}(\uptau)$.
Then, we have that

$\displaystyle\uprho_{\uptau}\left(y_{it}-\boldsymbol{X}_{(i)t-1}^{\star\top}\hat{\boldsymbol{\beta}}^{\star}(\uptau)\right)=\uprho_{\uptau}\left(\boldsymbol{\mathcal{E}}_{it\uptau}-\boldsymbol{X}_{(i)t-1}^{\star}\Delta^{\star}\left(\uptau\right)\right).$
(4.26)

where
$\boldsymbol{\mathcal{E}}_{it\uptau}=y_{it}-\boldsymbol{X}_{(i)t-1}^{\star\top}\boldsymbol{\beta}^{\star}(\uptau)$.
Then, the minimization of (4.26) is equivalent to minimizing for a fixed
$\uptau\in(0,1)$ such that

$\displaystyle\mathcal{Z}_{NT}^{o}\left(\Delta,\uptau\right)=\left(NT\right)^{-1}\sum_{i=1}^{N}\sum_{t=1}^{T}\left\\{\uprho_{\uptau}\left(V_{it\uptau}-\left(NT\right)^{-1/2}X^{\top}_{i(t-1)}\Delta\right)-\uprho_{\uptau}\left(V_{it\uptau}\right)\right\\}.$
(4.27)

∎

Suppose that regularity conditions hold for all $\uptau\in(0,1)$. Then for any
$\pi\in(0,1/2)$ we have that

$\displaystyle\underset{\uptau\in[\pi,1-\pi]}{\mathsf{sup}}\overset{d}{\to}\underset{\uptau\in[\pi,1-\pi]}{\mathsf{sup}}\left\\{\boldsymbol{M}^{\prime}(\uptau)\big{[}\boldsymbol{D}_{1,m}^{-1}(\uptau)\boldsymbol{D}_{0,m}\boldsymbol{D}_{1,m}^{-1}(\uptau)\big{]}\boldsymbol{M}(\uptau)\right\\}$
(4.28)

Monte Carlo experiments can be employed to evaluate the performance of the
nonstandard limiting distribution of the Wald statistic
$\mathcal{W}_{n,m}(\uptau)$ against the conventional $\chi^{2}-$asymptotic
approximation.

##  Conclusion

In this article, we study estimation and inference for nonstationary quantile
predictive regression models with generated regressors. We propose a two-stage
estimation procedure, where the generated regressor that corresponds to the
one-step ahead prediction of $\mathsf{VaR}$ based on persistent data is
generated in the first stage of the estimation procedure. In the second-stage
of the estimation procedure the two-stage estimator is employed to construct
the one-step ahead prediction of $\mathsf{CoVaR}$ based on the estimated
$\mathsf{VaR}$ and another set of persistent data not necessarily identical to
the first-stage estimation procedure. We investigate the asymptotic properties
of the two-stage estimator based on quantile predictive regression models. We
show that the variance-covariance matrix needs to be adjusted to account for
the estimation error in the first-stage estimation procedure.

In general, accounting for the presence of generated regressors when obtaining
forecasts of the pair $(\mathsf{VaR},\mathsf{CoVaR})$, in the proposed two-
stage estimation methodology, results in improved inference especially in the
presence of systemic risk in financial markets. This is done by adjusting the
variance-covariance matrix to to ensure that tests obtained from this two-
stage estimation procedure deliver accurate inference. In this article, we
present the main estimation steps of our proposed approach, and discuss main
challenges when deriving the asymptotic theory of estimators under regressors
nonstationarity (via the LUR parametrization).

In a further study we examine the econometric identification and estimation of
a novel SUR representation that corresponds to the node-specific quantile
predictive regression models proposed by Katsouris, (2021); Katsouris, 2023b ;
Katsouris, 2023c . Considering these node-specific interactions especially
under the assumption of an underline graph dependence structure (see,
Katsouris, 2023a ) is exactly that we can use SUR techniques to model tail
interdependencies in a network which has some interesting properties
especially when regressors are modelled using the local-to-unity
parametrization.

##### Conflicts of interest

The author declares that there are no known conflicts of interest.

##### Data availability

No data was used for the research described in this article.

##### Acknowledgements

Financial support from the Research Council of Finland (grant 347986) is
gratefully acknowledged. The author sincerely thank Dr. Ji Hyung Lee for
generously sharing the replication code for the IVX-QR procedure which has
been adapted in this article. The author also wish to thank Jose Olmo, Tassos
Magdalinos and Jayeeta Bhattacharya from the Department of Economics,
University of Southampton for helpful conversations and comments that
significantly improve this paper as well as Markku Lanne and Mika Meitz from
the Faculty of Social Sciences, University of Helsinki. All remaining errors
are my own responsibility.

##  Appendix

###  Appendix A: Auxiliary Results

######  Theorem 1.

(see Theorem 2 in Kato, (2009)) Suppose $f_{n}(x,\uptau)$ with $n\geq 1$ are
convex in $x$ for each $\uptau$ and bounded in $\uptau$ for each $x$.
Moreover, let
$g_{n}(x,\uptau)=-x^{\top}W_{n}(\uptau)+\frac{1}{2}x^{\top}Q(\uptau)x$, where
$\left\\{W_{n}(.)\right\\}$ is a sequence of bounded stochastic processes and
$Q(\uptau)$ is a $d\times d$ non-stochastic symmetric positive definite matrix
for each $\uptau$. Suppose that the maximum eigenvalue of $Q(\uptau)$ is
bounded from above and the minimum eigenvalue of $Q(\uptau)$ is bounded away
from 0 over $\uptau\in\mathcal{B}$. If it holds that

$\displaystyle\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\left|f_{n}(x,\uptau)-g_{n}(x,\uptau)\right|\overset{p}{\to}0,$
(6.1)

for each $x$ and of for every $\eta>0$, there exists a constant $M>0$ such
that

$\displaystyle\underset{n\to\infty}{\mathsf{lim\ sup}}\
\mathbb{P}\left(\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\left\lVert
W_{n}(\uptau)\right\rVert>M\right)\leq\eta,$ (6.2)

then
$x_{n}(\uptau)=\left\\{Q(\uptau)\right\\}^{-1}W_{n}(\uptau)+r_{n}(\uptau)$,
where $\mathsf{sup}_{\uptau\in\mathcal{B}}\left\lVert
r_{n}(\uptau)\right\rVert=o_{p}(1)$.

###### Proof.

Let $y_{n}(\uptau)=\left\\{Q(\uptau)\right\\}^{-1}W_{n}(\uptau)$, which is the
unique minimum point of $g_{n}(.,\uptau)$ for each $\uptau$. Then, it can be
shown that

$\displaystyle
g_{n}(x,\uptau)-g_{n}\left(y_{n}(\uptau),\uptau\right)=\frac{1}{2}\left(x-y_{n}(\uptau)\right)^{\top}Q_{n}(\uptau)\left(x-y_{n}(\uptau)\right)\geq
c\left\lVert x-y_{n}(\uptau)\right\rVert,$ (6.3)

for some constant $c>0$. Denote with

$\displaystyle\Delta_{n}=\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\
\underset{x:\left\lVert
x-y_{n}(\uptau)\right\rVert\leq\delta}{\mathsf{sup}}\left|f_{n}(x,\uptau)-g_{n}(x,\uptau)\right|.$
(6.4)

Therefore, letting $r_{n}(\uptau)=x_{n}(\uptau)-y_{n}(\uptau)$, we have that
$\mathbb{P}^{*}\left(\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\left\lVert
r_{n}(\uptau)\right\rVert>\delta\right)\leq\mathbb{P}^{*}\left(\Delta_{n}\geq(c\delta^{2})/2\right)$.
So it suffices to show that $\Delta_{n}\overset{p}{\to}0$. Let $\eta>0$ be an
arbitrary positive constant. Take $M>0$ such that the following bound holds

$\displaystyle\underset{n\to\infty}{\mathsf{lim\ sup}}\
\mathbb{P}^{*}\left(\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\ \left\lVert
y_{n}(\uptau)\right\rVert>M\right)\leq\eta.$ (6.5)

Define $K:=\left\\{x:\left\lVert x-y\right\rVert\leq\delta,\left\lVert
y\right\rVert\leq M\right\\}=\left\\{x:\left\lVert
x\right\rVert\leq\delta+M\right\\}$. Then, for every $\epsilon>0$,

$\displaystyle\mathbb{P}^{*}\left(\Delta_{n}>\epsilon\right)\leq\mathbb{P}^{*}\left(\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\
\underset{x\in K}{\mathsf{sup}}\
\big{|}f_{n}(x,\uptau)-g_{n}(x,\uptau)\big{|}>\epsilon\right)+\mathbb{P}^{*}\left(\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\left\lVert
y_{n}(\uptau)\right\rVert>M\right).$ (6.6)

Since by Lemma 1, we have that

$\displaystyle\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\ \underset{x\in
K}{\mathsf{sup}}\
\big{|}f_{n}(x,\uptau)-g_{n}(x,\uptau)\big{|}=\underset{\uptau\in\mathcal{B}}{\mathsf{sup}}\
\underset{x\in K}{\mathsf{sup}}\
\big{|}f_{n}(x,\uptau)+x^{\top}W_{n}(\uptau)-\frac{1}{2}x^{\top}Q(\uptau)x\big{|}\overset{p}{\to}0,$
(6.7)

We can conclude that $\underset{n\to\infty}{\mathsf{lim\
sup}}\mathbb{P}^{*}\left(\Delta_{n}>\epsilon\right)\leq\eta$. Since $\eta>0$
is arbitrary, then the proof ends. ∎

######  Example 4.

Notice that it can be shown that
$\hat{\beta}(\uptau)-\beta(\uptau)=O_{p}(n^{-1/2})$ uniformly in
$\uptau\in[\alpha,1-\alpha]$ for any $\alpha\in(0,1/2)$. Then, the
computational property of regression quantile processes is employed to derive
the uniform asymptotic representation of
$n^{1/2}\left(\hat{\beta}(\uptau)-\beta(\uptau)\right)$. Notice also that the
functional class is P-Donsker for any $B\subset\mathbb{R}^{p}$

$\displaystyle\big{\\{}g(y,x)=\big{(}\uptau-I(y\leq
x^{\top}\beta)\big{)}x_{j},\beta\in B,\uptau\in[\alpha,1-\alpha],1\leq j\leq
p\big{\\}}$ (6.8)

In the remaining of this section, we borrow some relevant results from the
framework proposed by Bhattacharya, (2020) (see Appendix of the paper).

######  Lemma 1 (Bhattacharya, (2020)).

Suppose that regularity conditions hold. Then, the following inequalities and
unique representations can be derived.

  * (i)

$Q^{(2)}(\beta;\uptau,\theta)$ is continuous with respect to its three
arguments, with

$\displaystyle\left\lVert
Q^{(2)}\big{(}\beta_{1};\uptau,\theta\big{)}-Q^{(2)}\big{(}\beta_{0};\uptau,\theta\big{)}\right\rVert\leq
C\left\lVert\beta_{1}-\beta_{2}\right\rVert$ (6.9)

for all $\beta_{0}$ and $\beta_{1},\theta\in\Theta$ and $\uptau\in[0,1]$.

  * (ii)

$Q^{(2)}(\beta;\uptau,\theta)$ is strictly positive for all $\beta\in
B(\theta),\theta\in B(\theta)$ and $\uptau\in[0,1]$.

  * (iii)

For $\theta\in\Theta$ and $\uptau\in[0,1]$ has a unique minimizer
$\beta(\uptau;\theta)$ which is continuously differentiable in $\theta$ and
$\uptau$ with

$\displaystyle\frac{\partial\beta(\uptau;\theta)}{\partial\theta^{\prime}}$
$\displaystyle=H(\uptau;\theta)^{-1}D(\uptau;\theta),$ (6.10)
$\displaystyle\frac{\partial\beta(\uptau;\theta)}{\partial\theta^{\prime}}$
$\displaystyle=H(\uptau;\theta)^{-1}\mathbb{E}\left[X(\theta)\right]$ (6.11)

where $H(\uptau;\theta)$ and $D(\uptau;\theta)$ are defined above.

###### Proof.

Consider a uniform order study of the Bahadur error term such that

$\displaystyle\widehat{\boldsymbol{\mathcal{E}}}\left(\uptau;\theta\right)$
$\displaystyle=\sqrt{n}\left(\widehat{\boldsymbol{\beta}}\left(\uptau;\theta\right)-\boldsymbol{\beta}\left(\uptau;\theta\right)\right)-\bigg{(}\boldsymbol{H}^{-1}\left(\uptau;\theta\right)\widehat{\boldsymbol{S}}\left(\uptau;\theta\right)\bigg{)}$
(6.12)
$\displaystyle\underset{(\uptau,\theta)\in[\uptau_{1},\uptau_{2}]\times\Theta}{\mathsf{sup}}\left\lVert\widehat{\boldsymbol{\mathcal{E}}}\left(\uptau;\theta\right)\right\rVert$
$\displaystyle=\mathcal{O}_{p}\left(\frac{\mathsf{log}^{3/4}n}{n^{1/4}}\right).$
(6.13)

Define the following optimization function

$\displaystyle\mathcal{L}_{n}\left(\gamma,\uptau;\theta\right)=\sum_{t=1}^{n}\left\\{\rho_{\uptau}\left(Y_{t}(\theta)-X_{t}(\theta)^{\top}\left(\frac{\gamma}{\sqrt{n}}+\boldsymbol{\beta}\left(\uptau;\theta\right)\right)\right)-\rho_{\uptau}\bigg{(}Y_{t}(\theta)-X_{t}(\theta)^{\top}\boldsymbol{\beta}\left(\uptau;\theta\right)\bigg{)}\right\\},$
(6.14)

This implies that a unique minimizer can be established given by the
expression below

$\displaystyle\sqrt{n}\left(\widehat{\boldsymbol{\beta}}\left(\uptau;\theta\right)-\boldsymbol{\beta}\left(\uptau;\theta\right)\right)=\underset{\gamma}{\mathsf{arg\
min}}\ \mathcal{L}_{n}\left(\gamma,\uptau;\theta\right).$ (6.15)

Relevant studies with the topological convergence as well as existence and
uniqueness proofs are mentioned in the Appendix of Katsouris, (2022) (see,
also Kato, (2009)). Thus, following the argumentation in Bhattacharya, (2020)
a uniform order for
$\widehat{\boldsymbol{\mathcal{E}}}\left(\uptau;\theta\right)$ relies on a
uniform order study for the remainder term
$\mathbb{R}\big{(}\gamma,\epsilon,\uptau;\theta\big{)}$. To do this, one can
employ concepts of maximal inequality under bracketing conditions as given in
the related literature as well as other linearisation techniques. There are
two techniques commonly used in the literature; the bracketing and set
covering approach as presented by Bhattacharya, (2020) and the partitioning
approach presented in several studies (see, Brian, (2021) and Goh and Knight,
(2009)). The remainder term can be expressed as

$\displaystyle\mathbb{R}_{n}\left(\gamma,\epsilon,\uptau;\theta\right)=\mathbb{L}_{n}\left(\gamma,\epsilon,\uptau;\theta\right)-\mathbb{L}_{n}^{0}\left(\gamma,\epsilon,\uptau;\theta\right)=\sum_{t=1}^{n}\mathbb{R}_{t}\left(\gamma,\epsilon,\uptau;\theta\right)$
(6.16)

∎

######  Lemma 2 (Bhattacharya, (2020)).

For real numbers $t_{\gamma}$ and $t_{\epsilon}>0$ with $t_{\epsilon}\geq 1$,
$t_{\epsilon}=\big{(}t\ \mathsf{log}^{3/4}n\big{)}/n^{1/4}$ for some $t>0$,
such that
$\left(t_{\gamma}+t_{\epsilon}\right)^{1/2}/t_{\epsilon}\leq\mathcal{O}\left(n^{1/4}/\mathsf{log}^{1/2}n\right)$,
for large $n$,

$\displaystyle\mathbb{E}\left[\underset{\left(\gamma,\epsilon,\uptau;\theta\right)\in\mathcal{B}\left(0,t_{\gamma}\right)\times\mathcal{B}\left(0,t_{\epsilon}\right)\times[\uptau_{1},\uptau_{2}]\times\Theta}{\mathsf{sup}}\
\bigg{|}\mathbb{R}_{n}^{1}\big{(}\gamma,\epsilon,\uptau;\theta\big{)}\bigg{|}\right]\leq
C\frac{\mathsf{log}^{1/2}n}{n^{1/4}}t_{\epsilon}(t_{\gamma}+t_{\epsilon})^{1/2}.$
(6.17)

######  Remark 11.

Within our framework the auxiliary results briefly discussed above can be
useful to derive error bounds to establish the order of convergence of
remainder terms when implementing a Bahadur representation to quantile-based
functionals. In addition, these result will require suitable modifications to
account for the nonstationar properties of regressors. For instance, OLS-based
residuals from the first-stage procedure depend on the nuisance parameter of
persistence thus stochastic equicontinuity arguments will need to be extended
to a suitable probability space that accounts for dynamic misspecification
(see, Kasparis and Phillips, (2012)).

A process $W_{NT}\left(\uptau\right)$ is said to be tight if and only if for
any $\delta>0$ there exists a compact set $E$ with

$\displaystyle\underset{\uptau\in E}{\mathsf{sup}}\
\mathbb{P}\big{(}W_{NT}\left(\uptau\right)\in E\big{)}>1-\delta.$ (6.18)

###### Proof.

We define with

$\displaystyle\psi_{1}(D)=-\frac{1}{\sqrt{NT}}\sum_{i,t}\big{[}\psi_{\uptau_{2}}\left(V_{it\uptau_{2}}\right)-\psi_{\uptau_{1}}\left(V_{it\uptau_{1}}\right)\big{]}$
(6.19)

for any interval $D=(\uptau_{1},\uptau_{2}]$. Therefore, to show the tightness
property, we adopt Theorem 15.6 from Billingsley, (1968) and prove a
sufficient Chentsov-Billingsley type of inequality (see, Bickel and Wichura,
(1971)) as follows. ∎

######  Lemma 3.

For any two intervals $D_{1}=(\uptau_{1},\uptau_{2}]$ and
$D_{2}=(\uptau_{2},\uptau_{3}]$ we have that

$\displaystyle\mathbb{E}\left[\big{\\{}\eta^{\top}\xi_{1}(D_{1})\big{\\}}^{2}\big{\\{}\eta^{\top}\xi_{1}(D_{2})\big{\\}}^{2}\right]\leq
C(\uptau_{3}-\uptau_{1}),$ (6.20)