# Nonparametric quantile regression for spatio-temporal processes

Soudeep Deb ID
Indian Institute of Management Bangalore
Cláudia Neves ID
King’s College London
Subhrajyoty Roy ID
Indian Statistical Institute Kolkata

###### Abstract

In this paper, we develop a new and effective approach to nonparametric
quantile regression that accommodates ultrahigh-dimensional data arising from
spatio-temporal processes. This approach proves advantageous in staving off
computational challenges that constitute known hindrances to existing
nonparametric quantile regression methods when the number of predictors is
much larger than the available sample size. We investigate conditions under
which estimation is feasible and of good overall quality and obtain sharp
approximations that we employ to devising statistical inference methodology.
These include simultaneous confidence intervals and tests of hypotheses, whose
asymptotics is borne by a non-trivial functional central limit theorem
tailored to martingale differences. Additionally, we provide finite-sample
results through various simulations which, accompanied by an illustrative
application to real-worldesque data (on electricity demand), offer guarantees
on the performance of the proposed methodology.

KEY WORDS AND PHRASES: Asymptotic theory, electricity demand forecasting,
gaussian processes, multivariate quantiles.

MSC2020 SUBJECT CLASSIFICATIONS: Primary 62M20, 62A99; secondary 60G15, 60G70.

##  1 Introduction

Within the context of various economic, environmental and societal research
problems, lies the common aim of wanting to estimate the effect of a covariate
on the behavior of the quantiles of a continuous response variable whilst it
is also well known that a simple regression model explaining the conditional
expectation of the response variable may turn out short in providing a
satisfactory description of the full spread of the targeted conditional
distribution, particularly at the level of tail-related data (cf. Reich,,
2012).

Originally proposed by Koenker and Bassett Jr, (1978), even in its most
simplified setting, quantile regression serves as a valuable alternative to
the familiar regression on conditional expectation, for it is not as heavily
reliant on the typical assumption about normally distributed errors for
assessing significance of the postulated model. Indeed, working off the
inherent nonparametric nature of conditional quantiles, quantile regression
proves as a reliable inference methodology to the extent of readily
accommodating heavy-tailed distributions for the errors, also being robust to
the presence of outliers (see Chaudhuri et al.,, 1997; Ezzahrioui and Ould-
Saïd,, 2008). In more complex situations, especially those concerning time
series modeling or involving space-time changing attributes, quantile
regression recognizably offers a plethora of opportunities to statistically
meaningful analyses of empirical data from real-life applications.

The leading example in this paper has dedicated focus on conditional quantiles
of electricity load profiles obtained from smart-meter records installed in
many UK households. The key motivation behind this illustrative application is
to gain improved understanding of peak-load events so as to help distribution
network providers maintain a stable supply chain at particular times of the
day or week and ward off rare but sometimes inevitable power outages (Jacob et
al.,, 2020). Like with many other fields in applied science, an important
difficulty arises from the dimension of the space-time covariate (thus
infinite-dimensional in theory) becoming increasingly large as the sample size
increases. We tackle this problem through a multivariate time series approach
under minimally restrictive assumptions on the data generating process, hence
making it amenable to translate to other applied settings given the added
flexibility for adjusting to different factors and confounding effects. Before
delving into both rationale and technical aspects of our proposed suite of
statistical methods concerning quantile regression on a spatio-temporal
dimension, a succinct account of extant literature in this context is
pertinent at this point.

There has been primarily two distinct nonparametric approaches towards
modeling conditional quantiles when dealing with a response variable at a
time. One such approach hinges on a pinball type loss function conditional on
the value of the covariate and views the resulting estimation problem from the
optimization modeling paradigm. In order to translate this into a viable and
tractable problem, however, one often needs to make a set of assumptions which
have typically ranged from a continuous functional of the covariate of
interest, to functionals belonging to a reproducible kernel Hilbert space
(Takeuchi et al.,, 2006), and to polynomial functions of the covariate (Zhou
and Wu,, 2009), all of which resulting in smooth quantile curves. A common
caveat to this type of approaches is that they can still suffer from the well-
known quantile crossing effect, which will lead to extra constraints yet to be
factored into the optimization setup. To combat this issue, a nonparametric
kernel estimation of the conditional distribution function was put forward
(cf. Ferraty et al.,, 2005; Ezzahrioui and Ould-Said,, 2008, and references
therein) which, upon a relatively straightforward inversion procedure, has had
proven efficacy in yielding conditional quantile curves that satisfy the non-
crossing criterion. However, because there is a necessary inversion step
involved in the latter, the former challenge now shifts to the smoothness
property of the estimated quantile curve which is not always attainable in
practice, thus posing a critical hurdle to the modeling of real-worldesque
data.

When looking at multivariate time series or spatio-temporal data, a case which
this paper seeks to address, we note that despite the vastly accumulating
literature on conditional quantile function estimation in the independent
setting, the development of flexible statistical methodology for the
underpinning stochastic processes is still relatively in its infancy. Progress
in these domains has been limited, possibly by the challenges faced when
accounting for the core dependency structure in a way that does not break down
at the laws of physics governing the process generating the data and therefore
is capable of producing realistic inference. To the best of our knowledge,
only a few successful attempts in this direction have been made thus far. A
couple of them (e.g. Reich et al.,, 2011; Das and Ghosal,, 2017) are Bayesian
approaches that tend to anchor on basis polynomials, with ensuing
considerations in eliciting prior distributions to the relevant parameters
alongside discussion on which allowances must be brought into practice to make
the method feasible. Because any misspecification involving the latter is
detrimental to the overall quality and feasibility of inference, the
possibility of undermining end-users’ confidence in the statistical methods
being employed cannot be wholly discarded. The approach taken by Reich,
(2012); Duan et al., (2021), among others, considers each quantile curve at
different locations separately, with the levels of the quantile curves not
affecting each other except for the non-crossing condition. However, this
implies that such approaches may fail in the situations where observations at
most locations follow a general trend (thus remaining close to the median)
with a few remainder locations exhibiting disparate behavior from the bulk.

In this paper, we introduce a nonparametric estimator of the conditional
quantile of a continuous response variable. Our statistical procedure has the
advantage of being able to deliver smooth estimates whilst demonstrating to
stave off the above-mentioned quantile crossing effect. The extremity issue is
overcome by modeling the multivariate quantiles of the relevant distribution
through an appropriate random element of the unit euclidean ball. The
asymptotic properties are developed as part of the theoretical study of the
proposed statistical methodology, both aimed at estimation and testing. We
consider the asymptotic framework where both the number of time-points and the
number of locations grow to infinity, which required significant adaptation of
the conventional delta-method used for deriving consistency and asymptotic
normality of estimators, since our approach leverages the sub-Gaussian
properties of the data to maintain higher-order moments under control for
ensuring tightness. Moreover, ensuring consistency has dictated that weak
convergence of the process must hold in a uniform sense in order to account
for the changing vector spaces with varying dimension, hence being the source
of additional complexity. This was handled by employing concentration
inequalities based on a sequence of compact sets for guaranteeing boundedness
in the differences between estimated and true quantiles. Lastly, to prove
asymptotic normality, we make use of Portmanteau’s theorem and show that, for
a rich class of smooth functions with bounded derivatives, the expectation of
the desired functional of the normalized estimators converges to the
expectation of the same functional value applied on the appropriate slice of
Gaussian process, which we also identify.

The structure of this paper is as follows. The statistical methodological
details and the key results developed are described in Sections 2, 3 and 4.
Following the exposition of asymptotic theory, we go on to developing
simultaneous confidence bands for the spatio-temporal quantile curves across
all locations and time-points in Section 5. Associated statistical tests of
hypotheses are also described therein. In order to maintain the flow of the
paper, we defer the proofs of the theoretical results until Section 6. In
Section 7, we present our maiden application of the proposed estimator to
electricity load forecasting based on smart-meter data collected as part of
the Thames Valley Vision project (details in Jacob et al.,, 2020). This
section focuses on inference, encompassing estimation, prediction and
hypothesis testing, with a particular emphasis on two testing procedures for
homogeneity. Numerical experiments that provide empirical evidence supporting
the validity of our statistical methods are conducted as a supplement to this
study. Notably, the fundamental large sample properties elucidated in the
paper are mirrored in the finite sample simulation study detailed in the
supplementary material accompanying this paper.

##  2 Spatio-temporal framework for quantile regression

This section presents the framework for the proposed quantile regression for
spatio-temporal processes, also furthering related works by De Gooijer et al.,
(2006); Zhou, (2010); Deb and Jana, (2021).

Let $\left\\{Y(t,s):\,t\in\mathcal{T},s\in\mathcal{S}\right\\}$ be a spatio-
temporal process indexed by temporal domain consisting of a (dense) temporal
horizon $\mathcal{T}=[0,1]$ and by a countable spatial horizon set
$\mathcal{S}\subset\mathbb{R}^{2}$. For all $d\in\mathbb{N}$, let
$\bm{X}:\mathcal{T}\rightarrow\chi\subseteq\mathbb{R}^{d}$ be a
$d$-dimensional temporal covariate process such that for each
$t\in\mathcal{T}$, $\bm{X}(t)=\left(X_{1}(t),\ldots,X_{d}(t)\right)^{\prime}$
is a predictor variable in $\chi\subseteq\mathbb{R}^{d}$. In the sliding scale
of generality, this covariate process can be allowed to be made fairly generic
so as to encapsulate external factors that might affect the temporal
dependence to the main process $Y(t,s)$, to being set out as specialized into
a given periodic component, for example, aiming to capture inherent
seasonality. Henceforth , our interest lies in the conditional space-time
process $Y(t,s)\mid\\{\bm{X}(t):t\in\mathcal{T}\\}$ with expectation and
covariance structure respectively denoted by

$\mathbb{E}(Y(t,s))=\mu(\bm{X}(t),s)\quad\mbox{ and
}\quad\mathrm{Cov}(Y(t,s_{1}),Y(t,s_{2}))=\Sigma(\bm{X}(t),s_{1},s_{2}),$
(2.1)

for every $s,s_{1},s_{2}\in\mathcal{S}$ and $t\in\mathcal{T}$. Note that the
temporal dependence in the process $Y(t,s)$ is essentially driven by that of
the associated covariate process $\bm{X}(t)$. Of course, it would be beguiling
to assume that $Y(t,s)$ can actually be observed in the way of covering the
whole set $\mathcal{T}\times\mathcal{S}$, so we must settle with $Y(t,s)$
indexed by the discrete subset of time points
$\mathcal{T}_{n}=\\{t_{1},t_{2},\dots t_{n}\\}$ and locations
$\mathcal{S}_{p}=\\{s_{1},s_{2},\dots s_{p}\\}$, where both $n$ and $p$ are
assumed sufficiently large, i.e. both tending to $\infty$ albeit in a
controlled manner for the latter which we impose from the outset as an
intermediate sequence $p=p_{n}$ such that $p/n\rightarrow 0$. This renders a
more realistic framework for inference with a practical resonance: for each
$t\in\mathcal{T}_{n}$, suppose that the data is generated according to the
nonparametric spatial regression model:

$\bm{Y}(t)=\begin{pmatrix}Y(t,s_{1})\\\ Y(t,s_{2})\\\ \vdots\\\
Y(t,s_{p})\end{pmatrix}=\bm{\mu}(\bm{X}(t))+\bm{\Sigma}^{1/2}(\bm{X}(t))\begin{pmatrix}\xi_{t1}\\\
\xi_{t2}\\\ \vdots\\\ \xi_{tp}\end{pmatrix},$ (2.2)

where $\bm{\mu}(\bm{X}(t))$ is a random vector in $\mathbb{R}^{p}$ of the form

$\bm{\mu}(\bm{X}(t)):=\begin{pmatrix}\mu(\bm{X}(t),s_{1})\\\
\mu(\bm{X}(t),s_{2})\\\ \vdots\\\ \mu(\bm{X}(t),s_{p})\end{pmatrix},$

the covariance matrix is given by $\bm{\Sigma}(\bm{X}(t))$, whose $(i,j)$-th
entry is $\Sigma(\bm{X}(t),s_{i},s_{j})$ for $1\leqslant i,j\leqslant p$, and
the (unobserved) errors $\left(\xi_{tj}\right)_{j=1,\dots p}$ are independent
and identically distributed with mean zero, positive density on an interval
containing $0$ (albeit not necessarily that of a normal distribution) and unit
variance. Hereafter, $\bm{\xi}(t)$ denotes the vector of $\xi_{tj}$’s.

The framework (2.2) lays bare that the key problem of estimating a spatial
quantile at any given time $t\in\mathcal{T}$ fits well into the logic of
estimating a $p$-variate quantile for the distribution of $\bm{Y}(t)$ centered
at time-points $t_{1},t_{2},\ldots,t_{n}$. It remains quite general on the
whole, and especially in two ways outlined as follows. Firstly, it
accommodates multivariate time series data even when there is no inherent
spatial structure in the coordinates but rather a rank correlation can be
inferred from domain knowledge, such as daily closing prices of multiple stock
indices. If relevant information about the spatial dependency structure is
available _a priori_ , it can be encapsulated in $\bm{\Sigma}(\cdot)$ through
a set of suitable constraints. These constraints range from being precisely
defined, fundamentally reliant on a certain postulated parametric model, to
being defined by a given class of (unknown) distributions in a manner that
echoes a familiar balancing act: the robustness-efficiency trade-off.
Secondly, such a formulation enables inference using multivariate quantiles
that may describe different degrees of extremity on different coordinates as
it allows for possible different tail behavior at nearby locations. This
justifies the focus on nonparametric estimation of the conditional
multivariate quantiles of $\bm{Y}(t)$, given the temporal covariate process
$\bm{X}(t)$, assumed predictable.

To this date, various definitions for a quantile of a multivariate
distribution have been put forward in the literature, notably those drawing on
convex hulls (Eddy,, 1985), minimum volume sets (Einmahl and Mason,, 1992),
geometric generalization of spatial median (Chaudhuri,, 1996) and center
outward quantile functions (Figalli,, 2018). Among these, the geometric
concept of a quantile in $\mathbb{R}^{d}$, $d\geqslant 2$, originating from
Chaudhuri, (1996) and endorsed in Konen and Paindaveine, (2022), stands out as
the most closely related to the present paper.

Indeed, central to this paper is the loss function
$\Phi(\bm{u},\bm{v})=\left\|\bm{v}\right\|+\langle\bm{u},\bm{v}\rangle$
($\left\|\bm{v}\right\|$ denotes the usual Euclidean norm of $\bm{v}$), for
all $\bm{u}\in B^{d-1}=\\{\bm{x}\in\mathbb{R}^{d}:\left\|\bm{x}\right\|<1\\}$
and the estimated quantile $\bm{v}\in\mathbb{R}^{d}$, in which setting a
recognizable form of the multivariate $\bm{u}$-quantile of a $d$-variate
random vector $\bm{X}$ beckons:

$Q(\bm{u}):=\arg\min_{\bm{q}\in\mathbb{R}^{d}}\,\mathbb{E}\left[\Phi(\bm{u},\bm{X}-\bm{q})\right],$
(2.3)

(cf. Chowdhury and Chaudhuri,, 2019). As desired, a choice of $\bm{u}\in
B^{d-1}$ in (2.3) towards the boundary, i.e. $\left\|\bm{u}\right\|\to 1$,
leads to a (high or low) extreme quantile, whereas fixing $\bm{u}$ near the
origin yields spatially centered measurable functionals, including expectiles
such as the mean vector ascribed to the $\mathcal{L}_{2}$-objective function
in (2.3) or the spatial median in connection with $\mathcal{L}_{1}$-norm.
Since the vector $\bm{u}$ depends on the dimension $p$ of the covariate, and
$p=p_{n}$ is such that $p\rightarrow\infty$ as $n\rightarrow\infty$ (in a
controlled manner stipulated before), we shall consider a spatial process $u$
(measurable with respect to the Lebesgue measure on $\mathbb{R}^{p}$)
satisfying

$\int_{\mathcal{S}}|u(s)|^{2}ds\leqslant 1,$

which, for every finite restriction $\mathcal{S}_{p}$ of $\mathcal{S}$, gives
rise to $\bm{u}|_{\mathcal{S}_{p}}=(u(s_{1}),\dots u(s_{p}))\in B^{p-1}$. In a
similar vein, we now define for every $t\in\mathcal{T},\bm{x}\in\chi$, the
$\bm{u}$-conditional spatial quantile of $\bm{Y}(t)$ given $\bm{X}(t)=\bm{x}$
as

$Q(\bm{u},\bm{x})|_{\mathbb{R}^{p}}=\arg\min_{\bm{q}\in\mathbb{R}^{p}}\,\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|+\langle\bm{u}|_{\mathcal{S}_{p}},\bm{Y}(t)-\bm{q}\rangle\mid\bm{X}(t)=\bm{x}\right],$
(2.4)

where $Q(\bm{u},\bm{x})|_{\mathbb{R}^{p}}$ highlights that the conditional
quantile process $Q(\bm{u},\bm{x})$ is indexed by a $p$-dimenstional spatial
domain $\mathcal{S}_{p}$. The set $\mathcal{S}_{p}$ has to be asymptotically
dense in $\mathcal{S}$ as $p\rightarrow\infty$, i.e., for any
$s\in\mathcal{S}$ and however small $\epsilon>0$, we can find a point
$s^{\ast}\in\mathcal{S}_{p}$ within the ball of radius $\epsilon$ centered at
$s$ for sufficiently large $p$. It is important to stress that the inclusion
of the temporal domain adds to a body of contemporary works in a substantive
way since, to the extent of our knowledge, research thus far refrained to the
study of nonparametric spatial quantile regression. Our proposal is
distinctive in that the criterion (loss) function presiding the M-estimation
contemplated in (2.4) is no longer a fixed function, but rather it is allowed
to vary with $p$ whilst preserving many of the fundamental properties for
quantile regression investigated in Chowdhury and Chaudhuri, (2019). This is
the key insight to expanding the notion of a spatial (multivariate) quantile
across a direction, furthering the aim of tackling the spatio-temporal
regression estimation problem stated in (2.4). The precise meaning of this
will become more apparent in Section 3.

In this instance, a natural question that comes to the fore is whether such a
definition of conditional quantile process (as in (2.4)) entails
identifiability across the set $\mathcal{S}_{p}$. For example, if
$\mathcal{S}^{\ast}\subseteq\mathcal{S}$ is a superset of $\mathcal{S}_{p}$
for some fixed $p$, then the restriction of
$Q(\bm{u},\bm{x})|_{\mathcal{S}^{\ast}}$ to the set of points
$\mathcal{S}_{p}$ should be same as $Q(\bm{u},\bm{x})|_{\mathbb{R}_{p}}$. To
see this, we consider
$\mathcal{S}^{\ast}=\mathcal{S}_{p}\cup(\mathcal{S}^{\ast}-\mathcal{S}_{p})=S_{1}\cup
S_{2}$ and correspondingly define $\bm{Y}_{1}(t),\bm{Y}_{2}(t)$ and
$\bm{q}_{1},\bm{q}_{2}$ as the restrictions of $\bm{Y}(t)$ and $\bm{q}$
respectively. Finally, we note that

$Q(\bm{u},\bm{x})|_{\mathcal{S}^{\ast}}=\operatorname*{arg\,min}_{\bm{q}_{1},\bm{q}_{2}}\mathbb{E}\left[\sqrt{\|\bm{Y}_{1}(t)-\bm{q}_{1}\|^{2}+\|\bm{Y}_{2}(t)-\bm{q}_{2}\|^{2}}\right.\\\
\left.+\langle\bm{u}|_{S_{1}},\bm{Y}_{1}(t)-\bm{q}_{1}\rangle+\langle\bm{u}|_{S_{2}},\bm{Y}_{2}(t)-\bm{q}_{2}\rangle\mid\bm{X}(t)=\bm{x}\right].$

It now clearly follows that optimizing with respect to
$\bm{q}_{1}\in\mathbb{R}^{p}$ yields a minimizer of (2.4) for any fixed value
of $\bm{q}_{2}$, thus ensuring identifiability.

##  3 Nonparametric spatio-temporal quantile estimation

For ease of introduction, we first address the case where we look at one
location at a time. This corresponds to setting $p=1$ in the basic model
(2.2), whereby the canonical nonparametric quantile estimator gives rise to
the kernel-weighted version of the quantile $q$ of order $\tau\in(0,1)$ of the
corresponding empirical probability distribution function of $Y(t)$ given
$\bm{X}(t)=\bm{x}$ on $\chi=\\{x_{1},\ldots,x_{n}\\}\subset\mathbb{R}^{d}$.
Formally, for $u\equiv u(\tau)=(2\tau-1)$, $0<\tau<1$, the quantile $q$ of
order $\tau$ is the minimizer of

$\mathcal{M}_{u,n}^{(1)}(q)=\sum_{i=1}^{n}k(\bm{x},\bm{X}(t_{i}))\left(|Y(t_{i})-q|+u\cdot(Y(t_{i})-q)\right),$
(3.1)

where $k(\cdot,\cdot)$ is a suitably chosen kernel function. In this respect,
we note that $\mathcal{M}_{u,n}^{(1)}$ is convex on $\mathbb{R}$. With
$\tau=1/2$ standing for the median as the minimizer of
$\mathcal{M}_{u,n}^{(1)}$, $u$ can be viewed as a gauge for central
displacement or tail-relatedness, since the greater the $|u|$ the more tail-
relatable the quantile. Clearly, finding minimizers of (3.1) configures an
$M$-estimation problem, and with the stated sample criterion function
incorporating a kernel function of bandwidth $b_{n}$, this is one of a
nonparametric nature. By this token, a generalization to the $p$-multivariate
case can be readily achieved by taking suitably normalized matrix-valued
kernel functions $\bm{K}(\cdot,\cdot)$ into the inner product that constitutes
(3.1), thus giving way to the following definition of a nonparametric spatio-
temporal $\bm{u}$-quantile that mirrors (2.4) in the associated sequence of
statistical experiments: for $\bm{u}|_{\mathcal{S}_{p}}\in\mathbb{R}^{p}$
comprised of the set $\\{u(s):s\in\mathcal{S}_{p}\\}$, define the estimator

$\hat{\bm{q}}_{n}(\bm{u},\bm{x}):=\arg\min_{\bm{q}\in\mathbb{R}^{p}}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}),$
(3.2)

where

$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}):=\frac{1}{n}\sum_{i=1}^{n}\left\\{\left\|\bm{Y}(t_{i})-\bm{q}\right\|_{\bm{K}(t_{i})}+\left\langle\bm{u}|_{\mathcal{S}_{p}},\bm{Y}(t_{i})-\bm{q}\right\rangle_{\bm{K}(t_{i})}\right\\},$
(3.3)

for which we write

$\displaystyle\bm{K}(t)$ $\displaystyle:=\bm{K}(\bm{x},\bm{X}(t))$
$\displaystyle\left\|\bm{Y}(t)-\bm{q}\right\|_{\bm{K}(t)}^{2}$
$\displaystyle:=\left\|\bm{K}(t)(\bm{Y}(t)-\bm{q})\right\|^{2},$
$\displaystyle\left\langle\bm{u},\bm{Y}(t)-\bm{q}\right\rangle_{\bm{K}(t)}$
$\displaystyle:=\bm{u}^{\intercal}\bm{K}(t)(\bm{Y}(t)-\bm{q}).$

As the objective function $\mathcal{M}_{\bm{u},n}^{(p)}$ in (3.3) has been
highlighted to depend on both the parameter dimension $p$ and sample size $n$,
in view of (2.4) its counterpart population version (for criterion function)
is defined as

$\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}):=\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|+\langle\bm{u}|_{\mathcal{S}_{p}},\bm{Y}(t)-\bm{q}\rangle\mid\bm{X}(t)=\bm{x}\right].$

In Section 4, this setting will be expanded with $p$ taken as a sequence of
positive integers indexed in $n$.

Before going any further, we spell out the key assumptions underpinning the
basic asymptotic results established in this paper; each of which will be
called upon as deemed appropriate and/or necessary throughout.

  1. (A1)

The covariate process $\bm{X}(t)$ is a predictable process, i.e., there exists
a sequence of random variables $\\{\epsilon_{t}:t=t_{1},t_{2},\dots\\}$, with
its natural filtration $\mathcal{F}_{t}$ independent of the random variables
$\xi_{tj}$, such that $\bm{X}(t_{i})$ is $\mathcal{F}_{t_{(i-1)}}$-measurable
for any $t_{i}\in\mathcal{T}_{n}$.

  2. (A2)

There exists a sequence of convex, open sets
$\mathcal{U}_{1}\subseteq\mathbb{R},\mathcal{U}_{2}\subseteq\mathbb{R}^{2},\ldots,\mathcal{U}_{p}\subseteq\mathbb{R}^{p}$
where the sequence of criterion functions $\mathcal{M}_{\bm{u},n}^{(p)}$,
defined as in (3.3) takes on the vector-valued $\bm{q}\in\mathcal{U}_{p}$, for
all $t_{i}\in\mathcal{T}_{n}$, and for any fixed $p\in\mathbb{N}$ not
depending on $n$. Moreover, we assume that
$\inf_{\bm{q}\in\mathcal{U}_{p}}\inf_{t_{i}\in\mathcal{T}_{n}}\|\bm{Y}(t_{i})-\bm{q}\|>\delta$
for sufficiently small $\delta>0$, independent of $n$ and $p$.

  3. (A3)

The kernel matrix $\bm{K}(\bm{x},\bm{X}(t))$ is a positive definite ($p\times
p$) matrix for all $t\in\mathcal{T}$ and $\bm{x}\in\chi$, and the support of
the conditional distribution of each $\bm{Y}(t_{1}),\dots\bm{Y}(t_{n})$,
provided $\bm{X}(t)=\bm{x}$ is not enclosed in a straight line in
$\mathcal{S}$ (cf. Chowdhury and Chaudhuri,, 2019).

Assumption (A1) on the predictable process $\bm{X}(t)$ makes it possible to
write $\bm{X}(t_{i})=h(\dots,\epsilon_{t_{(i-2)}},\epsilon_{t_{(i-1)}})$ for
some measurable function $h:\mathbb{R}^{\infty}\rightarrow\chi$. It is a
common assumption in related literature (see, e.g., Zhao and Wu,, 2008).
Assumption (A2) is an operational assumption, which in conjunction with (A3),
ensures that $\mathcal{M}_{\bm{u},n}^{(p)}$ consists of a sequence of random
convex functions on $\mathcal{U}_{p}$, for every $p\geqslant 1$; thereby
forming a sequence of continuous functions. This, in turn, ascertains the
existence of M-estimators $\widehat{\bm{q}}_{n}$ in the quality of a uniquely
defined sequence of minimizers of $\mathcal{M}_{\bm{u},n}^{(p)}$. We further
note that the objective function $\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ is
strictly convex in its argument $\bm{q}$, and the same conclusion is true for
$\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})$ regardless of the square-integrable
spatial process $u$, and for every fixed $p\in\mathbb{N}$. This strict
convexity ensures that the true quantile $\bm{q}_{0}$ is a well-separated
minimizer of $\mathcal{M}_{\bm{u}}^{(p)}$, that is, there exists a
$\bm{q}_{0}\in\mathcal{U}_{p}$ for which
$\inf_{\bm{q}:\,\|\bm{q}-\bm{q}_{0}\|\geqslant\delta}\,(\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}_{0}))>0$,
for arbitrarily small $\delta>0$. Assumption (A2) also ensures that the
observed points $\bm{Y}(t_{i})$ are excluded from the feasible set of
$\bm{q}$, to ensure well-definedness of the estimating equation introduced
later in Section 3.3.

###  3.1 Estimation algorithm

The idea behind the Iteratively Re-weighted Least Squares (IRLS) type
estimation as an effective means to achieve minimization in (3.2) is to
replace the objective function with an equivalent version whose iterative
minimization leads to the closest answer to the original problem of obtaining
a quantile estimator $\hat{\bm{q}}_{n}(\bm{u},\bm{x})$. Firstly, we expand the
sample criterion function $\mathcal{M}_{\bm{u},n}$ as

$\frac{1}{n}\sum_{i=1}^{n}\frac{\left(\bm{Y}(t_{i})-\bm{q}\right)^{\intercal}\bm{K}^{2}(t_{i})\left(\bm{Y}(t_{i})-\bm{q}\right)}{\|\bm{Y}(t_{i})-\bm{q}\|_{\bm{K}^{2}(t_{i})}}+\langle\bm{u},(\bm{Y}(t_{i})-\bm{q})\rangle_{\bm{K}(t_{i})},$

where we use the convention that $\bm{y}/\|\bm{y}\|=0$ if $\bm{y}=\bm{0}$. In
addition, we highlight that $\bm{q}$ does not change with $i$. Henceforth, we
shall use $w_{i}(\bm{q}):=\|\bm{Y}(t_{i})-\bm{q}\|^{-1}_{\bm{K}^{2}(t_{i})}$
so as to avoid notational burden. Minimization of the above with respect to
$\bm{q}$ yields the sequence of estimators

$\hat{\bm{q}}_{n}=\left[\sum_{i=1}^{n}w_{i}(\bm{q})\bm{K}^{2}(t_{i})\right]^{-1}\left\\{\frac{1}{2}\sum_{i=1}^{n}\bm{K}(t_{i})\bm{u}+\sum_{i=1}^{n}w_{i}(\bm{q})\bm{K}^{2}(t_{i})\bm{Y}(t_{i})\right\\}.$
(3.4)

The above expression suggests that, given an initial estimate
$\hat{\bm{q}}^{(0)}$, it is possible to use (3.4) iteratively in order to
improve estimation at the next iteration. Specifically,
$\widehat{\bm{q}}^{(k+1)}$ is obtained by calculating
$w_{i}(\widehat{\bm{q}}^{(k)})$ for all $i=1,2,\dots n$, using the current
estimate and then using the above equation to go to the next step of the
iteration. Despite the matrix $\sum_{i=1}^{n}w_{i}(\bm{q})\bm{K}(t_{i})$
guaranteed to be positive definite with probability tending to one, in
practice it may become numerically singular at some stage of the iterative
algorithm, particularly when the number of locations $p$ is large. By adding a
small positive number $\rho$ to the diagonal elements of this matrix, such a
numerical complication can be easily averted. The hyperparameter $\rho$
subsequently becomes the Lagrangian coefficient associated with a new
constraint for the objective function (3.2): $\|\bm{q}\|\leqslant C$, for some
fixed positive number $C$. This iterative procedure can also be viewed as a
Newton-Raphson type procedure, ensuring its convergence. In particular, it is
of straightforward verification that this iteration scheme reduces the value
of the objective function $\mathcal{M}_{\bm{u},n}^{(p)}(\widehat{\bm{q}})$ at
the next iteration, i.e.,
$\mathcal{M}_{\bm{u},n}^{(p)}(\widehat{\bm{q}}^{(k+1)})\leqslant\mathcal{M}_{\bm{u},n}^{(p)}(\widehat{\bm{q}}^{(k)})$
for all $k\in\mathbb{N}$. Combining this fact with the strict convexity of the
objective function, it follows that the solution of the algorithm converges to
the global minimizer of $\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ since it is
inherently assumed that the set of global minimizers is non-empty. This result
is formulated in the next Lemma.

######  Lemma 3.1.

Let $\bm{K}(t)=\bm{K}(\bm{x},\bm{X}(t))$ be a matrix-valued kernel function
such that the maximum eigenvalue of $\bm{K}(t)$ is uniformly bounded for all
$t\in[0,1]$. Then, with an initial estimate $\widehat{\bm{q}}^{(0)}$, the
sequence of improved estimators obtained by the iteration equation

$\widehat{\bm{q}}^{(k+1)}=\left[\sum_{i=1}^{n}w_{i}(\widehat{\bm{q}}^{(k)})\bm{K}^{2}(t_{i})\right]^{-1}\left\\{\frac{1}{2}\sum_{i=1}^{n}\bm{K}(t_{i})\bm{u}+\sum_{i=1}^{n}w_{i}(\widehat{\bm{q}}^{(k)})\bm{K}^{2}(t_{i})\bm{Y}(t_{i})\right\\},\
k\geqslant 0,$

converges as $k\rightarrow\infty$ and the ultimate estimator
$\widehat{\bm{q}}^{\ast}$ is the solution
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ of (3.2).

###  3.2 Choice of the kernel and associated bandwidth

The matrix-valued kernel function $\bm{K}(\bm{x},\bm{X}(t))$ involving the
bandwidth $b_{n}>0$ requires suitable normalizing constants $c_{n}>0$.
Specifically, we consider the matrix

$\bm{K}(\bm{x},\bm{X}(t))=c_{n}^{-1}\mathcal{K}\left(\dfrac{\bm{x}-\bm{X}(t)}{b_{n}}\right),$

given in terms of a symmetric and continuous kernel
$\mathcal{K}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p\times p}$, with compact
support on $B^{d-1}=\\{\bm{x}\in\mathbb{R}^{d}:\|\bm{x}\|\leqslant 1\\}$, and
in keeping with the formulation we typically encounter for such nonparametric
adjacency matrices of weights (e.g. Deb et al.,, 2017). Specific choices for
constants $b_{n},c_{n}>0$ must be based not only on theoretical
considerations, but also informed by empirical evidence on the behavior of
competing estimators in situations of repeated sampling, all considered
against the backdrop of assumptions threaded in the relevant subject
literature. For example, univariate nonparametric kernel estimators for a
trend in time series, which have been proposed in Zhao and Wu, (2008), appear
suitably normalized via $c_{n}=\sqrt{n\,b_{n}}$. We also assume that
$\mathcal{K}(0)=c_{n}\bm{I}_{p}$, similar to the univariate case where the
kernel function at origin is taken to be equal to $1$. Given the choice of a
particular form of kernel $\mathcal{K}$, one may choose $c_{n}$ to be a
function of the bandwidth $b_{n}$ such that the eigenvalues of the normalized
kernel matrix $\bm{K}(\bm{x},\bm{X}(t))$ remain bounded away from $0$ and
infinity.

Since the problem of estimating the conditional kernel tackled in here is
closely related to that of conditional density estimation, existing methods
with proven efficacy can be invoked (see e.g. Einbeck and Tutz,, 2006),
alongside the use of graphical approaches such as an empirical scatterplot of
the maximum eigenvalue of kernel matrix over different choices of $b_{n}$ to
obtain the rate of growth of the eigenvalues. It is also possible to couple
our results with the theoretical findings on the spectrum of kernel matrices
for popular kernels (e.g. Gaussian kernel, Matérn kernel) presented in
Barthelmé and Usevich, (2021), albeit the choice of the bandwidth $b_{n}$
should be moderated by a cross-validation approach that verifies the
theoretical conditions set out, which in this particular research are outlined
in Section 4. A complete treatment of these and analogous conditions requires
technical apparatus of which a detailed study can be found in Fan et al.,
(1996).

###  3.3 Key properties of the nonparametric quantile estimator

The results encompassing this section, summarized in Lemma 3.2 and ensuing
considerations, are viewed as fundamental steps to establishing performance
guarantees of the proposed nonparametric methodology for quantile estimation.
Indeed, a desirable property of a quantile estimator is affine equivariance:
if the samples $\bm{Y}(t)$ are changed to $\bm{A}\bm{Y}(t)+\bm{b}$ for some
non-singular matrix $\bm{A}$, the estimator should also exhibit the same
affine transformation. We note that the proposed quantile estimator
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ resulting from (3.2) is equivariant
under location transformation, i.e., employing a shift to observations through
$\bm{Y}(t_{i})+\bm{a}$ changes the sequence of estimators accordingly to
$\hat{\bm{q}}_{n}(\bm{u},\bm{x})+\bm{a}$. Although the equivariance also holds
in general with respect to a scale transformation $c\bm{Y}(t_{i})$, for some
scalar $c\in(0,\infty)$, verification of the orthogonal equivariance property
(or more general affine equivariance property) is dependent on the
incorporating kernel matrix $\bm{K}(t_{i})$ (see Chaudhuri,, 1996; Chen,,
1996; Konen and Paindaveine,, 2022).

As highlighted in Section 1, quantile estimators should generally exhibit two
key properties: non-crossing of the quantiles and smoothness of the estimator
as a function of the data. In the following discussion, we demonstrate how
upon a very minimal set of assumptions, both these properties follow in
relation to $\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$.

For each $n\geqslant 1$, the quantile estimator
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ is continuous with respect to
observables $\bm{Y}=\left(\bm{Y}(t_{1}),\dots\bm{Y}(t_{n})\right)$. In view of
assumption (A2), we can without loss of generality restrict the optimization
to a compact set $\mathcal{Q}\subseteq\mathbb{R}^{p}$ containing the true
targeted quantile. Consequently,
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})=\bm{q}(\bm{Y})=\operatorname*{arg\,min}_{\bm{q}\in\mathcal{Q}}\mathcal{M}^{(p)}_{\bm{u},n}(\bm{q},\bm{Y})$.
Now, let $\bm{Y}_{r}\rightarrow\bm{Y}_{0}$ be a converging sequence with
$r=1,2,\dots$, and let $\bm{q}_{r}=\bm{q}(\bm{Y}_{r})$ for every
$r\in\mathbb{N}$. We proceed with a subsequencing argument: for any
$\\{\bm{q}_{r_{k}}\\}_{k\in\mathbb{N}}$, the compactness of $\mathcal{Q}$
ensures that $\\{\bm{q}_{r_{k}}\\}$ has an accumulation point
$\bm{q}^{\prime}$, say. Also, by definition of $\bm{q}_{r_{k}}$, we have that
$\mathcal{M}^{(p)}_{\bm{u},n}(\bm{q}_{r_{k}},\bm{Y}_{r_{k}})\leqslant\mathcal{M}^{(p)}_{\bm{u},n}(\bm{q}_{0},\bm{Y}_{r_{k}}),\
r_{k}\in\mathbb{N}$. Since $\mathcal{M}^{(p)}_{\bm{u},n}$ is a continuous
function of $\bm{q}$ and $r_{k}\rightarrow\infty$, then
$\mathcal{M}_{\bm{u},n}(\bm{q}^{\prime},\bm{Y}_{0})\leqslant\mathcal{M}_{\bm{u},n}(\bm{q}_{0},\bm{Y}_{0})$,
thus entailing $\bm{q}^{\prime}=\bm{q}_{0}$ for any subsequence, from which
the continuity of the estimator ensues.

The next lemma outlines the desirable non-crossing property in any quantile
estimation methodology (e.g. in line with Theorem 2.1.2 of Chaudhuri,, 1996).
It is mainly a technical lemma aiming to justify that the M-estimator
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$, viewed as the approximate (asymptotic)
solution of the convex optimization problem (3.2), gains the configuration of
a $Z$-estimation problem for every fixed $(\bm{u},\bm{x})$. This is envisaged
to simplify substantially the proofs for asymptotic normality. Nevertheless,
because we will need to take into account the spatial dependence structure on
the basis of an infinitely growing number of locations, the rate of
convergence will differ from the familiar $\sqrt{n}$ for tightness mustered by
predominant asymptotic theory for $Z$-estimators (see van der Vaart,, 1998;
Kosorok,, 2008).

######  Lemma 3.2.

Under assumptions (A1) to (A3), for every $n,p\in\mathbb{N}$,
$\bm{x}\in\mathbb{R}^{d}$ and $\bm{u}\in B^{p-1}$, the sequence of minimizers
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ of the objective function
$\mathcal{M}_{\bm{u},n}^{(p)}$ (as in (3.3)) satisfies

$\frac{1}{n}\sum_{i=1}^{n}\bm{K}(t_{i})\left[\dfrac{\bm{K}(t_{i})(\bm{Y}(t_{i})-\hat{\bm{q}}_{n}(\bm{u},\bm{x}))}{\|\bm{K}(t_{i})(\bm{Y}(t_{i})-\hat{\bm{q}}_{n}(\bm{u},\bm{x}))\|}+\bm{u}\right]=\bm{0},$
(3.5)

for
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})\notin\\{\bm{Y}(t_{1}),\dots,\bm{Y}(t_{n})\\}$.
The true quantile $Q(\bm{u},\bm{x})$ in the quality of minimizer of
$\mathcal{M}_{\bm{u}}^{(p)}$ verifies the population version of (3.5),
specifically

$\mathbb{E}_{\bm{Y}(t)\mid\bm{X}(t)=\bm{x}}\left[\dfrac{\bm{Y}(t)-Q(\bm{u},\bm{x})}{\|\bm{Y}(t)-Q(\bm{u},\bm{x})\|}+\bm{u}\right]=\bm{0}.$
(3.6)

The actionable idea stemming from Lemma 3.2 – that we may regard (3.5) as an
estimating equation for obtaining $\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ from a
finite sample data with $n$ time-points and $p$ spatial points – gives rise to
an alternative Newton-Raphson type algorithm for obtaining the estimator.
Notice that, for $\bm{u}=0$, the estimating equation simply corresponds to a
weighted spatial median estimate of $\bm{Y}(t_{i})$ values, weighted by the
covariate process $\\{\bm{X}(t):t\in\mathcal{T}\\}$ through the kernel matrix
$\bm{K}(t)=\bm{K}(\bm{x},\bm{X}(t))$. However, for $\bm{u}\neq\bm{0}$, there
is an additional term creeping up in the way of an approximation bias, thereby
increasing or decreasing the coordinates of
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ as appropriate in a data-driven manner.
This is found pivotal to ensuring that the estimated quantile curves returned
by our methodology enjoy the non-crossing property. In view of (3.6), it is
clear that if a single coordinate $u_{i}$ of $\bm{u}$ is positive and
increases in value, then the corresponding coordinate of
$\bm{Y}(t)-Q(\bm{u},\bm{x})$ must be negative and decrease further in
expectation, resulting in $Q(\bm{u},\bm{x})$ turning away from the center
$\mathbb{E}(\bm{Y}(t))$. This justifies the non-crossing property.

##  4 Asymptotic results

The nonparametric quantile estimator introduced in Section 3 enjoys a number
of desirable large-sample properties. This section is aimed at addressing
these by allowing the number of spatial locations $p$ tend to infinity, albeit
in a controlled manner relatively to $n\rightarrow\infty$, as mentioned
earlier. The obtained asymptotic properties provide the theoretical
justification for the distinctive features of the quantile estimation and
ensuing testing methodology proposed in this paper.

We now fix some more notations in order to progress. A bounded sequence of
numbers $\\{x_{n}\\}$ is throughout written as $x_{n}=\mathcal{O}(1)$,
springing out to the notation $x_{n}=\mathcal{O}(y_{n})$ if
$x_{n}/y_{n}=\mathcal{O}(1)$. In general, if the random sequence $X_{n}$ is
tight, we write $\|X_{n}\|=\mathcal{O}_{p}(1)$, and if
$\|X_{n}\|/\|Y_{n}\|\,\displaystyle{\mathop{\longrightarrow}^{P}}\,0$ as
$n\rightarrow\infty$, we denote this fact by $X_{n}=o_{p}(Y_{n})$. In case of
deterministic sequence, $x_{n}=o(y_{n})$ denotes the convergence
$\limsup_{n\rightarrow\infty}|x_{n}/y_{n}|=0$.

Let $\mathcal{T}_{n}=\\{t_{1},\dots,t_{n}\\}$ denote the set of observable
time-points that satisfy the smoothing condition $\max_{2\leqslant i\leqslant
n}|t_{i}-t_{i-1}|\rightarrow 0$ as $n\rightarrow\infty$. We also require that
$\mathcal{T}_{n}$ is asymptotically dense in the temporal horizon
$\mathcal{T}$. In other words, for any $\tilde{t}\in\mathcal{T}$, there exists
a sequence of points $\tilde{t}_{n}\in\mathcal{T}_{n}$ such that
$|\tilde{t}-\tilde{t}_{n}|\rightarrow 0$.

######  Theorem 4.1.

Assume conditions (A1)-(A3) hold. Suppose that the number of spatial locations
$p$ is such that $p=p_{n}=o(n^{1/2})$ and $p_{n}^{-1}\log n=\mathcal{O}(1)$,
as $n\rightarrow\infty$. Then, for all $\bm{x}$ lying in the compact set
$\chi$, the sequence of estimators $\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ of
the population (theoretical) quantile $Q(\bm{u},\bm{x})$ satisfies

$\|\widehat{\bm{q}}_{n}(\bm{u},\bm{x})-Q(\bm{u},\bm{x})\|=\mathcal{O}_{p}(p/\sqrt{n}),$
(4.1)

provided a square-integrable spatial process $\bm{u}$, i.e.
$\int_{\mathcal{S}}|u(s)|^{2}ds<1$.

In Theorem 4.1, the number of spatial locations $p$ is allowed to tend to
infinity but at a vanishing rate relatively to $n^{1/2}$. It turns out that
such a mild restriction on the growth rate of $p_{n}$ is just what it takes
for ensuring the consistency of the quantile estimator paving the way for
attaining a non-degenarate weak limit.

For ascertaining the asymptotic normality of the proposed quantile estimator,
however, the number of spatial locations must align with the initial
requirement that $\mathcal{S}_{p}$ becomes asymptotically dense in
$\mathcal{S}$. Analogously to the temporal domain, we fix the setting where
the observable spatial points $\mathcal{S}_{p}=\\{s_{1},s_{2},\dots,s_{p}\\}$
satisfy the condition $\max_{1\leqslant j\neq k\leqslant
n}d(s_{j},s_{k})\rightarrow 0$ as $p\rightarrow\infty$. Here, $d(s_{j},s_{k})$
represents the geodesic distance between the spatial coordinates $s_{j}$ and
$s_{k}$. This ensures that $\mathcal{S}_{p}$ becomes dense in $\mathcal{S}$
with a growing number of spatial locations $p_{n}$ based on the geodesic
metric. Also, the bandwidth $b_{n}$ of the kernel matrix, architected and
discussed in Section 3.2 as $b_{n}\rightarrow 0$, must be such that
$nb_{n}\rightarrow\infty$ as $n\rightarrow\infty$ so that the maximum
eigenvalue of $\bm{K}(\bm{x},\bm{X}(t))$ is uniformly tight.

###  4.1 Assumptions

While assumptions (A1) to (A3) in Section 3 lay the groundwork for this paper,
the next set of assumptions (B1) to (B4) are deemed sufficient for attaining a
Gaussian limiting distribution for the proposed quantile estimator.

Let $f_{t}$ denote the joint probability density function of $\bm{X}(t)$ and
$\mathcal{G}_{k}$ and $\mathcal{G}_{k-}$ be the $\sigma$-algebra generated by
$\\{\bm{Y}(t,s),\bm{X}(t):t\leqslant t_{k},s\in\mathcal{S}\\}$ and
$\\{\bm{Y}(t,s):t<t_{k},s\in\mathcal{S}\\}\cup\\{\bm{X}(t):t\leqslant
t_{k},s\in\mathcal{S}\\}$ respectively. Additionally, let $\mathcal{P}_{k}$ be
the projection operator defined as
$\mathcal{P}_{k}(Z):=\mathbb{E}(Z\mid\mathcal{G}_{k})-\mathbb{E}(Z\mid\mathcal{G}_{k-1})$,
for all $k\in\mathbb{Z}$.

  1. (B1)

Suppose that the joint probability density function $f_{t}$ is differentiable
with respect to each coordinate of $\bm{X}(t)$, for each $t\in\mathcal{T}$.
This takes care of finite-dimensional distributions of the observable process.

  2. (B2)

The bandwidth $b_{n}>0$ satisfies the following boundary condition as
$n\rightarrow\infty$:

$\max\left(b_{n}^{4/3}\log n,\,\frac{(\log n)^{3}}{nb_{n}^{3}},\frac{(\log
n)^{2}}{b_{n}^{4/3}}\frac{\Xi_{n}}{n^{2}}\right)=o(1),$

with the cumulative effect of the first error $\bm{\xi}(0)$ in predicting the
covariate process $\bm{X}(n)$ defined as
$\Xi_{n}:=n\Theta_{2n}^{2}+\alpha_{n}$ (see Zhao and Wu,, 2008) where

$\displaystyle\Theta_{n}$
$\displaystyle=\sum_{i=1}^{n}\sup_{x\in\chi}\|\mathcal{P}_{0}(f_{t}(\bm{x}\mid\mathcal{G}_{i}))\|+\sup_{x\in\chi}\|\mathcal{P}_{0}(f_{t}^{\prime}(\bm{x}\mid\mathcal{G}_{i}))\|,$
$\displaystyle\alpha_{n}$
$\displaystyle=\lim_{m\rightarrow\infty}\sum_{k=n}^{m}\left(\Theta_{n+k}-\Theta_{k}\right)^{2}.$

  3. (B3)

There exists a constant $c>1$ such that

$\dfrac{1}{c}<\sup_{t\in\mathcal{T}}\dfrac{\lambda_{1}(\bm{K}(t))}{\lambda_{p}(\bm{K}(t))}<c,$

where $\bm{K}(t)=\bm{K}(\bm{x},\bm{X}(t))$ is the $p\times p$ kernel matrix as
before and $\lambda_{1}(\bm{A})$ and $\lambda_{p}(\bm{A})$ respectively denote
the largest and the smallest eigenvalues of the matrix $\bm{A}$.

  4. (B4)

The number of spatial points $p_{n}$ is such that $p_{n}b_{n}\rightarrow 0$,
$p_{n}^{9/8}b_{n}\rightarrow\infty$ and $p_{n}^{17/8}/n\rightarrow 0$, as
$n\rightarrow\infty$.

Assumption (B1) is natural and indeed a customary assumption that ensures
smoothness of the probability density functions, allowing to employ the
standard techniques in asymptotic analysis. Assumption (B2) is a further
restriction to the growth of the bandwidth $b_{n}$ of the kernel matrix, which
exerts influence on the strength of dependence of the process. From (B2), it
results that $b_{n}^{-1}(n^{-1/3}\log n)=o(1)$, thus imposing a lower bound on
$b_{n}$, and $b_{n}(\log n)^{3/4}=o(1)$, which provides an upper bound for
$b_{n}$. A popular choice for $b_{n}>0$ is that of the type $c_{0}n^{-\beta}$
for some (fixed) constants $c_{0}>0$ and $\beta>0$. If it is used in tandem
with the kernel proposed in Section 3.2, we find that $b_{n}=n^{-1/5}$, a
sequence that has been singled out in Zhao and Wu, (2008) for producing
satisfactory results.

Next, we note that the third term embedded in the triplet in (B2) implies
$\Xi_{n}/n^{2}=o(1)$ due to the earlier condition $b_{n}^{-1}(n^{-1/3}\log
n)=o(1)$ in conjunction with $nb_{n}\rightarrow\infty$. Because each summand
of the partial $\Theta_{n}$ gives the contribution of $\epsilon_{0}$ (as in
(A1)) to the associated predictor $\bm{X}(t_{i+1})$, this term plays a
relevant role in moderating the dependence of $\bm{X}(t)$. In particular, the
random vectors $\bm{X}(t_{i})$ must satisfy a summability property designed to
maintain the integrity of the predictor process $\bm{X}(t)$ by way of ensuring
that the covariance structure of the spatio-temporal process is well-defined
and can always be estimated. On the whole, this assumption has a number of
parallels with published works on inference for spatio-temporal processes (see
e.g. Zhao and Wu,, 2008; Zhou and Wu,, 2009; Deb and Jana,, 2021)). According
to Zhao and Wu, (2008), a condition such as this corrals virtually the full
spectrum for the strength of dependence, from the boundary case of the
predictor process $\bm{X}(t)$ exhibiting only short-range dependence (SRD,
determined by $\lim_{n\rightarrow\infty}\Xi_{n}<\infty$) to a weakening long-
range dependence (LRD).

Assumption (B3) on the other hand is required in order to control random
oscillations of the kernel matrix. In keeping with the random matrix theory,
Wigner’s law asserts that the eigenvalues of the $p\times p$ covariance matrix
$\bm{\Sigma}(\bm{X}(t))$ are of order $\sqrt{p}$ (Vershynin,, 2010) for
sufficiently large $p$, hence monotonically increasing with $n$ and $p$. This
means that for any fixed $t\in\mathcal{T}$, the ratio of the largest and the
smallest eigenvalue of $\bm{K}(t)$ remains bounded for sufficiently large $p$.
For the first finitely many terms, boundedness of this ratio is ensured by
tightness and positive definiteness of the kernel matrix $\bm{K}(t)$.
Assumption (B3) is a stronger version of this statement that requires uniform
boundedness over all choices of $t\in\mathcal{T}$. This assumption guarantees
that there exists a suitable kernel matrix such that the eigenvalues of the
covariance matrices of the kernel-weighted observables $\bm{Y}(t)$ remain
bounded. Finally, assumption (B4) in conjunction with the bounds on bandwidth
$b_{n}$ that ensue from (B2) gives that $(\log n)^{-2/3}\ll
p_{n}\ll\left(n/\log n\right)^{1/3}$, for large enough $n$. This speaks to the
ultra-dimension element in our methods.

###  4.2 Asymptotic distribution

In this section, we establish the Gaussian limiting distribution attained by
the proposed quantile estimator. The consistency of the proposed quantile
estimator expounded in Theorem 4.1 and accompanying text herald the rate at
which the number of locations $p$ can grow in relation to the number of time-
points $n$ so as to make this nonparametric quantile estimation feasible.

Continuing along the lines of (3.3) and adding to the notation introduced as
part of Lemma 3.2, we define

$\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}):=-\dfrac{1}{n}\sum_{i=1}^{n}\bm{K}(t_{i})\left[\dfrac{\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q})}{\|\bm{Y}(t_{i})-\bm{q}\|_{\bm{K}(t_{i})}}+\bm{u}\right].$

Let $\nabla^{2}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})$ be the matrix whose
$i$-th row is comprised of the gradient of
$\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ with respect to the $i$-th
coordinate of $\bm{q}$. Also, let us denote by
$\nabla^{3}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ the third order tensor whose
$(i,j,k)$-th entry is the third order derivative of
$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ with respect to the $i,j$ and $k$-th
coordinate of $\bm{q}$. We have the following proposition.

######  Proposition 1.

Under assumptions (A1)-(A3) and (B1)-(B4), the following hold.

  1. (P1)

There exists a $p$-variate Gaussian random variable $\bm{Z}(\bm{q}_{0})$ with
mean parameter $\eta_{t}(\bm{q}_{0})$ and dispersion matrix
$\Omega_{t}(\bm{q}_{0})$, such that

$\sup_{\bm{u}\in\mathcal{B}_{p}}\sup_{h\in\mathcal{H}}\left\|\mathbb{E}\left[h\left((pb_{n})^{-1}\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})\right)\right]-\mathbb{E}\left[h(\bm{Z}(\bm{q}_{0})\right]\right\|=\mathcal{O}_{p}(p^{2}/n),$

where $\mathcal{B}_{p}$ is a sequence of compact subsets of $B^{p-1}$,
$\mathcal{H}$ is the set of bounded thrice differentiable functions from
$\mathbb{R}^{p}$ to $\mathbb{R}$ with a bounded derivative up to the third
order.

  2. (P2)

There exists a collection of positive definite matrices
$\\{\Psi_{t}(\bm{q}_{0})\\}$ for all possible $\bm{q}_{0}\in\mathcal{Q}_{p}$
as a minimizer of (2.4) with respect to the spatial points $\mathcal{S}_{p}$,
such that

$\sup_{\bm{q}_{0}\in\mathcal{Q}_{p}}\left\|\nabla^{2}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})-\Psi_{t}(\bm{q}_{0})\right\|\xrightarrow{P}0.$

  3. (P3)

There exists $M>0$ (not depending on $n$) such that

$\sup_{\bm{q}\in\mathcal{Q}_{p}}\left\|\nabla^{3}(\mathcal{M}_{\bm{u},n}^{(p)})(\bm{q})\right\|\leqslant
M.$

The three parts (P1)-(P3) embodying Proposition 1, whose proof is deferred to
the supplementary material, mirror each of the usual necessary steps to
demonstrating the asymptotic normality of $M$-estimators (Huber,, 2011). The
asymptotic development for the proposed estimator of the conditional quantile
of $\bm{Y}(t)$, given $\bm{X}(t)=\bm{x}$ for any fixed $t\in\mathcal{T}$,
stays within this realm. Although general in their principle, the innovation
in Proposition (P1) is that it amplifies complementary works on quantile
regression by enabling the increasing dimension $p=p_{n}\rightarrow\infty$, as
$n\rightarrow\infty$. The simpler i.i.d. version analogous to this proposition
can be traced back to Portnoy, (1986), and contemporary related results
accommodating a changing $p_{n}$ can be found in Das and Lahiri, (2020). Our
infinite-dimensional framework shall build on a combination of (pre-
asymptotic) moment-type results approached by a functional central limit
theorem for martingale differences (cf. Belloni and Oliveira,, 2018).

In order to lessen the notational burden throughout the remainder of this
paper, we collect the relevant building blocks in the following:

$\displaystyle\bm{\eta}(t_{i})$
$\displaystyle=\mathbb{E}\left[\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})\mid\mathcal{G}_{t_{i}-}\right]=\bm{K}(t_{i})\left(\bm{\mu}(\bm{X}(t_{i}))-\bm{q}_{0}\right),$
(4.2) $\displaystyle\bm{V}_{1}(t_{i})$
$\displaystyle=\mathrm{Var}\left[\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})\mid\mathcal{G}_{t_{i}-}\right]=\bm{K}(t_{i})\bm{\Sigma}(\bm{X}(t_{i}))\bm{K}(t_{i}),$
(4.3) $\displaystyle\bm{V}_{2}(t_{i})$
$\displaystyle=\left(\dfrac{1}{\|\bm{\eta}(t_{i})\|}\bm{I}_{p}-\dfrac{\bm{\eta}(t_{i})\bm{\eta}(t_{i})^{\intercal}}{\|\bm{\eta}(t_{i})\|^{3}}\right),$
(4.4)

where $\bm{I}_{p}$ denotes the identity matrix of order $p$. Additionally, the
expectation and covariance pertaining to the finite-dimensional distributions
of the quantile process are encapsulated in:

$\displaystyle\bm{\eta}_{t}(\bm{q}_{0}):=$
$\displaystyle\dfrac{1}{npb_{n}}\sum_{i=1}^{n}\bm{K}(t_{i})\left[\dfrac{\bm{\mu}(\bm{X}(t))-\bm{q}_{0}}{\|\bm{\mu}(\bm{X}(t))-\bm{q}_{0}\|}-\dfrac{\bm{\eta}(t_{i})}{\|\bm{\eta}(t_{i})\|}\right],$
(4.5) $\displaystyle\bm{\Omega}_{t}(\bm{q}_{0}):=$
$\displaystyle\dfrac{1}{n^{2}p^{2}b_{n}^{2}}\sum_{i=1}^{n}\bm{K}(t_{i})\bm{V}_{2}(t_{i})\bm{V}_{1}(t_{i})\bm{V}_{2}(t_{i})^{\intercal}\bm{K}(t_{i}),$
(4.6) $\displaystyle\bm{\Psi}_{t}(\bm{q}_{0}):=$
$\displaystyle\dfrac{1}{n}\sum_{i=1}^{n}\bm{K}(t_{i})\bm{V}_{2}(t_{i})\bm{K}(t_{i}).$
(4.7)

We arrive at the main result in this paper, stated in the following theorem.

######  Theorem 4.2.

Consider $\bm{q}_{0}:=Q(\bm{u},\bm{x})$. Under the assumptions (B1)-(B4) and
in the setup of Theorem 4.1, suppose $p^{17/8}/n\rightarrow 0$, as
$n\rightarrow\infty$. Then, with $\mathcal{H}$ from (P1), it holds that

$\sup_{\bm{u}\in\mathcal{B}_{p}}\sup_{h\in\mathcal{H}}\left\|\mathbb{E}\Bigl{[}h\left((pb_{n})^{-1}(\widehat{\bm{q}}_{n}(\bm{u},\bm{x})-\bm{q}_{0}\right)\Bigr{]}-\mathbb{E}\Bigl{[}h\left(-\bm{\Psi}_{t}(\bm{q}_{0})^{-1}\bm{Z}(\bm{q}_{0})\right)\Bigr{]}\right\|\xrightarrow{P}0,$

where $\mathcal{B}_{p}$ is any sequence of compact subsets of $B^{p-1}$,
$\bm{Z}(\bm{q}_{0})$ is a $p$-variate Gaussian random vector with mean vector
$\bm{\eta}_{t}(\bm{q}_{0})$ and dispersion matrix
$\bm{\Omega}_{t}(\bm{q}_{0})$ furnished in (4.5) and (4.6).

At this point, we tactically turn focus to sketching a proof for Proposition 1
since this forms the basis for the proof of Theorem 4.2. For the first part in
reference, we start by noting that the challenge in obtaining the asymptotic
distribution of $(pb_{n})^{-1}\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})$
stems from the non-linearity of the normalized term
$\bigl{(}\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})\bigr{)}/\|\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})\|$
with respect to the random variables $\bm{Y}(t_{i})$. A standard trick to
thwart this issue is to employ a suitable Taylor’s expansion for obtaining a
sufficiently sharp approximation to a linear functional of $\bm{Y}(t_{i})$ or
$\bm{\xi}(t_{i})$. Indeed, the expressions (4.2)-(4.4) have already been
carved out in such a way to prepare for such an expansion: writing
$\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})=\bm{\eta}(t_{i})+\bm{V}_{1}^{1/2}(t_{i})\bm{\xi}(t_{i})$.
Application of Taylor’s expansion on the errors $\xi_{tj}$ thus precipitates
the development

$\dfrac{\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})}{\|\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q}_{0})\|}=\dfrac{\bm{\eta}(t_{i})}{\|\bm{\eta}(t_{i})\|}+\bm{V}_{2}(t_{i})\bm{V}_{1}^{1/2}(t_{i})\bm{\xi}(t_{i})+\mathcal{O}\left(\dfrac{\|\bm{V}_{1}^{1/2}(t_{i})\|}{\|\bm{\eta}(t_{i})\|^{3}}\|\bm{\xi}(t_{i})\|^{2}\right),$

i.e.,

$\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})=\dfrac{1}{n}\sum_{i=1}^{n}\bm{K}(t_{i})\left(\dfrac{\bm{\eta}(t_{i})}{\|\bm{\eta}(t_{i})\|}+\bm{u}\right)+\dfrac{1}{n}\sum_{i=1}^{n}\bm{K}(t_{i})\bm{V}_{2}(t_{i})\bm{V}_{1}^{1/2}(t_{i})\bm{\xi}(t_{i})\\\
+\dfrac{1}{n}\sum_{i=1}^{n}\mathcal{O}(\|\bm{K}(t_{i})\|\|\bm{\eta}(t_{i})\|^{-3}\|\bm{V}_{1}^{1/2}(t_{i})\|\|\bm{\xi}(t_{i})\|^{2}).$

The proof now unfolds in three steps. Owing to Lemma 3.2, the first term can
be shown to converge to a constant which constitutes the bias. Next,
application of a bespoke central limit theorem for the multidimensional
martingale provided by Belloni and Oliveira, (2018), ensures that the
asymptotic weak limit to the second term is a Gaussian random variable with
mean zero. Finally, by assumption (B4), the third term vanishes to zero at the
rate $\mathcal{O}(\max\\{p^{-9/4}b_{n}^{-2},p^{-3/2}b_{n}^{-1}\\})$. Further
details concerning this proof for Proposition (P1) are expounded in Appendix
A1 as part of the supplementary material accompanying this paper.

As for the proof of Proposition (P2), it is here sketched as follows (the
proof in its entirety is deferred to Appendix A2). The second order derivative
of the objective function $\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ is calculated
in the form:

$\nabla^{2}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})=\dfrac{1}{n}\sum_{i=1}^{n}\bm{K}(t_{i})\biggl{[}\dfrac{1}{\left\|\bm{Y}(t_{i})-\bm{q}\right\|_{\bm{K}(t_{i})}}\bm{I}_{p}-\dfrac{\bm{K}(t_{i})(\bm{Y}(t_{i})-\bm{q})(\bm{Y}(t_{i})-\bm{q})^{\intercal}\bm{K}(t_{i})}{\left\|\bm{Y}(t_{i})-\bm{q}\right\|^{3}_{\bm{K}(t_{i})}}\biggr{]}\bm{K}(t_{i}).$

After splitting the above expression into a sum of martingale difference terms
and a remainder term, the former sum, suitably normalized, can be shown to
converge to zero via application of a martingale central limit theorem (see
Belloni and Oliveira,, 2018). The dominant component in the remainder term
boils down to $\bm{\Psi}_{t}(\bm{q}_{0})$ with everything else being
asymptotically negligible. This is achieved through a dedicated Taylor’s
expansion.

Finally, regarding Proposition (P3), we find it convenient and useful to prove
the stronger statement that $\nabla^{2}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$
is Lipschitz with a Lipschitz constant $M$ independent of $n$ and $p$.
Starting with assumption (A2), we note that since $\mathcal{Q}_{p}$ is a
compact subset of $\mathcal{U}_{p}$,
$\|\bm{Y}(t_{i})-\bm{q}\|_{\bm{K}(t_{i})}\geqslant\lambda_{p}(\bm{K}(t_{i}))\delta$
for all $\bm{q}\in\mathcal{Q}_{p}$. Since the kernel matrix is positive
definite and has bounded entries, it readily follows that
$\nabla^{2}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ is a Lipschitz function of
$\bm{q}$ with a Lipschitz constant
$\bigl{(}\delta^{-1}\sup_{t_{i}\in\mathcal{T}_{p}}\lambda_{1}(\bm{K}(t_{i}))/\lambda_{p}(\bm{K}(t_{i}))\bigr{)}^{-3}$.
By assumption (B3), this quantity remains bounded.

##  5 Statistical inference: interval estimation and hypotheses testing

While the estimation of quantile curves yielding from (3.2) possesses several
desirable qualities, including the consistency and asymptotic normality, the
underpinning convergence is pointwise in the temporal domain, that is, it is
obtained for each value of the covariate process $\\{\bm{X}(t)=\bm{x}\\}$.
From the practical stance, an interesting approach resides in obtaining
simultaneous confidence band. Inevitably, this will require extending the
asymptotic properties of the estimator to the entire temporal domain,
uniformly on the index range of $\bm{X}(\cdot)$. The means to achieve this has
a number of parallels with the techniques exploited in Zhao and Wu, (2008).

To this end, we shall consider a set of values for the covariates,
sufficiently apart in order that at every time-point $t$, the kernel estimator
is affected by at most one of the values in that set. This is to ensure that
no correlation between estimates are at play for sufficiently distant time-
points. Upon the already obtained asymptotic normality results, this
consideration opens way to making uniform statistical inference possible:
asymptotic behavior of the proposed nonparametric estimator’s sequence across
multiple time-points, uniformly, is now within reach. The ultimate aim is to
devise meaningful asymptotic results pertaining to partial maxima of this
sequence of estimators by drawing on extreme value theory.

We define $\chi_{n}=\\{\bm{X}(t):t\in\mathcal{T}_{n}\\}\subseteq\chi$ as the
range of the covariate process $\\{\bm{X}(t)\\}$. We assume that, for any
$\bm{x},\bm{x}^{\prime}\in\chi_{n}$,
$\|\bm{x}-\bm{x}^{\prime}\|>2b_{n}\rightarrow 0$ thus implying that $\chi_{n}$
becomes dense within $\chi$ as $n$ tends to infinity. The cardinality of
$\chi_{n}$ is encapsulated in $m_{n}$.

###  5.1 Simultaneous confidence bounds

For the sake of simplicity, the results comprising in this section refrain to
$\bm{u}$-quantiles with $\bm{u}=(2\tau-1)\bm{1}_{p}/\sqrt{p}$ for some
$\tau\in(0,1)$. Henceforth, we denote this as $\bm{u}_{\tau}$. This framework
essentially renders quantiles that spread evenly among the considered spatial
locations. We shall denote the corresponding quantile estimator by
$\widehat{\bm{q}}_{n}(\tau,\bm{x})$ (or sometimes simply by
$\widehat{\bm{q}}(\tau,\bm{x})$) whereas the true targeted quantile is denoted
by $Q(\tau,\bm{x})$. We emphasize that this slight restriction to the
framework is imposed for convenience, but all results derived in this section
can be adapted to incorporate the more general multivariate case as before, in
a relatively straightforward manner.

By construction, it follows from the definition of $\chi_{n}$ that, for any
two points $\bm{x}_{i_{1}},\bm{x}_{i_{2}}\in\chi_{n}$, $i_{1}\neq i_{2}$, and
every $\delta>0$, the exists $n_{0}\in\mathbb{N}$ such that
$\left\|\bm{x}_{i_{1}}-\bm{x}_{i_{2}}\right\|<2\delta$, $i_{1},\,i_{2}>n_{0}$.
Hence, for every $t\in\mathcal{T}_{n}$, $\|\bm{X}(t)-\bm{x}_{j}\|<\delta$,
$j=i_{1},i_{2}$, which results from the triangle inequality. Since
$\mathcal{K}$ has bounded support in the unit hypersphere $B^{d-1}$, it thus
follows that the cross products of the weights intervening in the objective
function resolve to zero, i.e.

$\bm{K}(\bm{x}_{i_{1}},\bm{X}(t))\bm{K}(\bm{x}_{i_{2}},\bm{X}(t))=c_{n}^{-2}\mathcal{K}\left(\dfrac{\bm{X}(t)-\bm{x}_{i_{1}}}{b_{n}}\right)\mathcal{K}\left(\dfrac{\bm{X}(t)-\bm{x}_{i_{2}}}{b_{n}}\right)=0$
(5.1)

(the reader may wish to recall the discussion on $c_{n}$ as part of Section
3.2). Equation (5.1) lays bare that fixing two quantile levels $\tau$ and
$\tau^{\prime}$ ($0<\tau,\tau^{\prime}<1$), so that objective function if
optimized for $\widehat{\bm{q}}(\tau,\bm{x}_{i_{1}})$ contains the terms with
weights $\bm{K}(\bm{x}_{i_{1}},\bm{X}(t))$ for some $t\in\mathcal{T}_{n}$,
whereas for the sequence of estimators
$\widehat{\bm{q}}(\tau^{\prime},\bm{x}_{i_{2}})$ these weights amount to
$\bm{K}(\bm{x}_{i_{2}},\bm{X}(t))$, results in the two objective functions
(both for $\widehat{\bm{q}}(\tau,\bm{x}_{i_{1}})$ and
$\widehat{\bm{q}}(\tau^{\prime},\bm{x}_{i_{2}})$) being uncorrelated, thus
implying that the estimators $\widehat{\bm{q}}(\tau,\bm{x}_{i_{1}})$ and
$\widehat{\bm{q}}(\tau^{\prime},\bm{x}_{i_{2}})$ are also uncorrelated.

Therefore, by taking a multitude of points on the set
$\chi_{n}=\\{\bm{x}_{i_{1}},\bm{x}_{i_{2}},\dots\bm{x}_{i_{m_{n}}}\\}$, we
obtain from Theorem 4.2 that the estimated process
$\\{\hat{q}_{n}(\tau,\bm{x})(s):s\in\mathcal{S}\\}$ becomes asymptotically a
Gaussian process on both temporal and spatial domains, for every
$\bm{x}\in\chi_{n}$. More concretely, on the finite set of spatial points
$S_{d}=\\{s_{1},s_{2},\dots s_{d}\\}$ and with unit vector
$\bm{l}\in\mathbb{R}^{d}$, we have for every $\bm{x}\in\chi_{n}$, as
$n\rightarrow\infty$, that

$R_{n}(\bm{x},\bm{l},S_{d})=(pb_{n})^{-1}\dfrac{\bm{l}^{\intercal}(\hat{q}_{n}(\tau,\bm{x})|_{S_{d}}-Q(\tau,\bm{x})|_{S_{d}}-\bm{\mu}_{q}(\bm{x})|_{S_{d}})}{\sqrt{\bm{l}^{\intercal}\bm{\Sigma}_{q}(\bm{x})|_{S_{d}}\bm{l}}}\xrightarrow{d}\mathcal{N}(0,1).$
(5.2)

We highlight that $\hat{q}_{n}(\tau,\bm{x})|_{S_{d}}$ and
$Q(\tau,\bm{x})|_{S_{d}}$ respectively denote spatial process
$\hat{q}_{n}(\tau,\bm{x})$ and $Q(\tau,\bm{x})$ restricted to the spatial
locations in $S_{d}$. Letting $q_{0}=Q(\tau,\bm{x})$, then
$\bm{\mu}_{q}(\bm{x})|_{S_{d}}$ and $\bm{\Sigma}_{q}(\bm{x})|_{S_{d}}$
respectively denote mean and covariance structure, evaluated at the spatial
locations in the set $S_{d}$, of the Gaussian process
$\\{G(\bm{x},s):\bm{x}\in\chi,s\in\mathcal{S}\\}$ characterized by

$\begin{split}\mathbb{E}\left[G(\bm{x},s)\right]&=\bm{\Psi}_{\bm{X}(t)=\bm{x}}^{-1}({q}_{0})\bm{\eta}_{\bm{X}(t)=\bm{x}}({q}_{0})\\\
\text{Cov}\left(G(\bm{x},s),G(\bm{x}^{\prime},s)\right)&=\begin{cases}\bm{\Psi}_{\bm{X}(t)=\bm{x}}^{-1}({q}_{0})\bm{\Omega}_{\bm{X}(t)=\bm{x}}({q}_{0})\bm{\Psi}_{\bm{X}(t)=\bm{x}}^{-1}({q}_{0}),&\bm{x}=\bm{x}^{\prime};\\\
0,&\bm{x}\neq\bm{x}^{\prime}.\end{cases}\end{split}$ (5.3)

where the terms $\bm{\eta}({q}_{0}),\bm{\Omega}({q}_{0})$ and
$\bm{\Psi}({q}_{0})$ are evaluated at $\bm{X}(t)=\bm{x}$ as defined in
(4.5)-(4.7).

We now proceed to obtaining the asymptotic behavior of the supremum of
$G(\bm{x},s)$ over the temporal index through
$\sup_{\bm{x}\in\chi_{n}}\hat{q}_{n}(\tau,\bm{x})(s)$ for any
$s\in\mathcal{S}$. From the case in point that
$\hat{q}_{n}(\tau,\bm{x}_{i_{1}})$ and
$\hat{q}_{n}(\tau^{\prime},\bm{x}_{i_{2}})$ are asymptotically uncorrelated,
the next step is to apply the extreme value theorem upon the univariate
martingale difference sequence $R_{n}(\bm{x}_{i_{k}},\bm{l},S_{d})$,
$k=0,1,\dots,m_{n}-1$. Owing to the rich theory of extreme values, there
exists normalizing functions $C_{t}$ specific to the folded normal
distribution,

$C_{t}(z):=\frac{1}{\sqrt{2\log t}}\left(z-\frac{\log\log
t+\log(4\pi)}{2}\right)$ (5.4)

(see Resnick,, 1987, p.72), such that the weak convergence to the Gumbel
distribution function is attained in the following sense: for all
$z\in\mathbb{R}$,

$\
\lim_{n\rightarrow\infty}\mathbb{P}\left(\sup_{\bm{x}\in\chi_{m_{n}}}\sup_{\left\|\bm{l}\right\|=1}|R_{n}(\bm{x},\bm{l},S_{d})|\leqslant
C_{m_{n}}(z+2\log m_{n})\right)=\exp\\{-e^{-z}\\}.$ (5.5)

Next, we consider the identification

$\sup_{\left\|\bm{l}\right\|=1}\dfrac{\bm{l}^{\intercal}\bm{A}\bm{l}}{\bm{l}^{\intercal}\bm{B}\bm{l}}=\lambda_{1}\left(\bm{B}^{-1/2}\bm{A}\bm{B}^{-1/2}\right),$

where $\bm{A}$ and $\bm{B}$ are positive definite matrices. In conjunction
with the limiting distribution statement in (5.5), we subsequently find that

$\lim_{n\rightarrow\infty}\mathbb{P}\biggl{(}\sup_{\bm{x}\in\chi_{m_{n}}}\dfrac{1}{pb_{n}}\lambda_{1}\left[\bm{\Sigma}_{q}(\bm{x})^{-1/2}|_{S_{d}}\bm{e}(\bm{x})|_{S_{d}}\,\bm{e}(\bm{x})|_{S_{d}}^{\intercal}\bm{\Sigma}_{q}(\bm{x})^{-1/2}|_{S_{d}}\right]\\\
<C_{m_{n}}\left(z+2\log m_{n}\right)\biggr{)}=\exp\\{-e^{-z}\\},$ (5.6)

by taking supremum over all possible unit vectors $\bm{l}$ and whereby
$\bm{e}(\bm{x})|_{S_{d}}=(\hat{q}_{n}(\tau,\bm{x})|_{S_{d}}-Q(\tau,\bm{x})|_{S_{d}}-\bm{\mu}_{q}(\bm{x})|_{S_{d}})$.
Hence, under the assumptions (A1)-(A3) and (B1)-(B4), with a judicious choice
of $\chi_{m_{n}}$, i.e., so long as $\log m_{n}=o\left((pb_{n})^{-2}\right)$,
the sequence of estimators $\hat{q}_{n}(\tau,\bm{x})$ can be guaranteed to
fulfill (5.5) and (5.6) on a finite set of spatial locations
$S_{d}=\\{s_{1},s_{2},\dots s_{d}\\}$ from the “larger” set $\mathcal{S}$.
This result, which resonates with that in Theorem 2 of Zhao and Wu, (2008),
serves as basis to introducing simultaneous confidence bands for the
$\tau$-quantile curves $Q(\tau,\bm{x})$ ascribed to multiple values
$\bm{x}\in\chi_{n}$. Because $\chi_{n}$ becomes increasingly dense in $\chi$
as $n$ tends to $\infty$, an asymptotic confidence band for the quantile curve
$\\{Q(\tau,\bm{x}):\bm{x}\in\chi\\}$ arises from (5.6).

A contrast-like linear combination of the quantile curves across different
locations stems from (5.5) as follows: with a proper specification of the set
$S_{d}$ and vector $\bm{l}\in\mathbb{R}^{d}$ so that
$\bm{l}^{\intercal}Q(\tau,\bm{x})$ becomes the desired contrast-like quantile
curve, we find that the asymptotic simultaneous
$100(1-\alpha^{*})\%$-confidence bounds for
$\bm{l}^{\intercal}Q(\tau,\bm{x})|_{S_{d}}$ consist of

$\bm{l}^{\intercal}\hat{q}_{n}(\tau,\bm{x})|_{S_{d}}-\bm{l}^{\intercal}\widehat{\bm{\mu}}_{q}(\bm{x})|_{S_{d}}\pm\left[C_{m_{n}}\left(-\log\left(-\log(1-\alpha^{*}\right))+2\log
m_{n}\right)\right](pb_{n})\left(\bm{l}^{\intercal}\widehat{\bm{\Sigma}}_{q}(\bm{x})|_{S_{d}}\bm{l}\right)^{1/2},$

with embedded consistent estimators $\widehat{\bm{\mu}}_{q}(\bm{x})$ and
$\widehat{\bm{\Sigma}}_{q}(\bm{x})$ for the bias term $\bm{\mu}_{q}(\bm{x})$
and covariance $\bm{\Sigma}_{q}(\bm{x})$, respectively. We note that the
symmetry in the confidence interval is justified by the same limiting extreme
value distribution of $\inf_{\bm{x}}R_{n}(\bm{x},\bm{l},S_{d})$ and
$-\sup_{\bm{x}}\left\lvert R_{n}(\bm{x},\bm{l},S_{d})\right\rvert$, a remark
that will prove useful to constructing one-sided hypothesis testing
procedures.

###  5.2 Two tests

The statement in (5.6) can be adapted to the purpose of hypothesis testing
aimed at statistically validating structural assumptions on the quantile
curves across all locations simultaneously. Under the null hypothesis about
the structural restriction of the quantile curve $Q(\tau,\bm{x})$ at some
specific quantile $\tau\in[0,1]$, we can calculate the observed value of the
test statistic

$\bm{e}_{0}(\bm{x})|_{\mathcal{S}_{p}}=(\hat{q}_{n}(\tau,\bm{x})|_{\mathcal{S}_{p}}-\hat{q}_{0}(\tau,\bm{x})|_{\mathcal{S}_{p}}-\widehat{\bm{\mu}}_{q}(\bm{x})|_{\mathcal{S}_{p}})$
(5.7)

where $\hat{q}_{0}(\tau,\bm{x})$ stands as the estimate of $Q(\tau,\bm{x})$
under the postulated null hypothesis $H_{0}$ regarding a stable regime, and
$\widehat{\bm{\mu}}_{q}(\bm{x})$ is the estimate of the mean of the Gaussian
process introduced in (5.3) under the null hypothesis. This can be done by
replacing $\bm{\mu}(\bm{X}(t))$ and $\bm{q}_{0}$ by the estimated median
$\widehat{\bm{q}}_{n}(\bm{0},\bm{X}(t))$ and the estimated quantile
$\widehat{\bm{q}}_{n}(\bm{u},\bm{X}(t))$ under $H_{0}$ in the corresponding
terms. Then, the null hypothesis is rejected if (5.7) exhibits large
oscillation, i.e., at the nominal level of the test $\alpha^{*}$, $H_{0}$ is
rejected if

$\dfrac{1}{pb_{n}}\sup_{\bm{x}\in\chi_{n}}\lambda_{1}\left(\widehat{\bm{\Sigma}}_{q}(\bm{x})^{-1/2}|_{\mathcal{S}_{p}}\bm{e}_{0}(\bm{x})|_{\mathcal{S}_{p}}\bm{e}_{0}(\bm{x})|_{\mathcal{S}_{p}}^{\intercal}\widehat{\bm{\Sigma}}_{q}(\bm{x})^{-1/2}|_{\mathcal{S}_{p}}\right)$

exceeds the critical value
$c_{\alpha^{*}}=C_{m_{n}}\left(-\log\left(-\log(1-\alpha^{*})\right)+2\log
m_{n}\right)$.

A structural assumption that is often put to testing is that of homogeneity of
the covariate effect across the spatial domain, which, if true, will
automatically narrow down the spatio-temporal structure of the data to a
temporal one. This is statistically parlayed in the null hypothesis
$H_{0}^{(1)}:\,Q(\tau,\bm{x})=\beta(\bm{x})\bm{u}_{\tau}$, for all
$\tau\in[0,1]$. Under $H_{0}^{(1)}$, the estimator $\widehat{\beta}(\bm{x})$
arises from (3.5) as a solution of:

$\sum_{i=1}^{n}\bm{K}(\bm{x},\bm{X}(t_{i}))\left[\dfrac{\bm{K}(\bm{x},\bm{X}(t_{i}))(\bm{Y}(t_{i})-\widehat{\beta}(\bm{x})\bm{u}_{\tau})}{\|\bm{K}(\bm{x},\bm{X}(t_{i}))(\bm{Y}(t_{i})-\widehat{\beta}(x)\bm{u}_{\tau})\|}+\bm{u}_{\tau}\right]=0,$

and the null hypothesis is rejected for large values of the test statistic,
i.e., if

$\dfrac{1}{pb_{n}}\sup_{\bm{x}\in\chi_{n}}\lambda_{1}\left(\widehat{\bm{\Sigma}}_{q}(\bm{x})^{-1/2}|_{\mathcal{S}_{p}}\bm{e}_{1}(\bm{x})|_{\mathcal{S}_{p}}\bm{e}_{1}(\bm{x})|_{\mathcal{S}_{p}}^{\intercal}\widehat{\bm{\Sigma}}_{q}(\bm{x})^{-1/2}|_{\mathcal{S}_{p}}\right)>c_{\alpha^{*}}$

on the relevant set $\mathcal{S}_{p}$ of spatial locations, where
$\bm{e}_{1}(\bm{x})|_{\mathcal{S}_{p}}=(\hat{q}_{n}(\tau,\bm{x})|_{\mathcal{S}_{p}}-\widehat{\beta}(x)\bm{u}_{\tau}-\widehat{\bm{\mu}}_{q}(x)|_{\mathcal{S}_{p}})$
can be defined.

The second testing procedure with practical importance under the realm of
(2.2) ascribes time to a covariate through $\bm{X}(t)=t$. The hypothesis of
homogeneity over time then translates into
$H_{0}^{(2)}:\,Q(\bm{u}_{\tau},t)/Q(\bm{u}_{\tau},1)=t$, for all $t\in(0,1]$
and $\tau\in[0,1]$. At the significance level $\alpha^{*}$, $H_{0}^{(2)}$ is
rejected if

$\dfrac{1}{pb_{n}}\sup_{t\in(0,1]}\lambda_{1}\left(\widehat{\bm{\Sigma}}_{q}(t)^{-1/2}|_{\mathcal{S}_{p}}\bm{e}_{2}(t)|_{\mathcal{S}_{p}}\bm{e}_{2}(t)|_{\mathcal{S}_{p}}^{\intercal}\widehat{\bm{\Sigma}}_{q}(t)^{-1/2}|_{\mathcal{S}_{p}}\right)>c_{\alpha^{*}}\,,$

with associated
$\bm{e}_{2}(t)|_{\mathcal{S}_{p}}=(\hat{q}_{n}(\tau,t)|_{\mathcal{S}_{p}}-t\hat{q}_{n}(\bm{u}_{\tau},1)|_{\mathcal{S}_{p}}-\widehat{\bm{\mu}}_{q}(x)|_{\mathcal{S}_{p}})$.

##  6 Proofs

This section is encompassed of the proofs of the two main theorems in this
paper as well as of their ancillary results. We begin by addressing Lemma 3.2,
which consolidates the estimating equation for the proposed nonparametric
quantile estimator.

Proof of Lemma 3.2: Under assumption (A3),
$\bm{K}(t)=\bm{K}(\bm{x},\bm{X}(t))$ is a positive definite matrix for any
$t\in\mathcal{T}$. Hence,
$\Phi_{K}(\bm{u},\bm{Y}(t)-\bm{q})=\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}$
is a strictly convex function of $\bm{q}$. Since
$\widehat{\bm{q}}=\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ is the minimizer of
$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})=\frac{1}{n}\sum_{i=1}^{n}\Phi_{K}(\bm{u},\bm{Y}(t_{i})-\bm{q})$,
then for every $\bm{v}\in\mathbb{R}^{d}$,

$\displaystyle\lim_{\epsilon\downarrow
0}\dfrac{1}{n\epsilon}\sum_{i=1}^{n}\left[\Phi_{K}(\bm{u},\bm{Y}(t_{i})-\widehat{\bm{q}}+\epsilon\bm{v})-\Phi_{K}(\bm{u},\bm{Y}(t_{i})-\widehat{\bm{q}})\right]\geqslant
0$ $\displaystyle\Rightarrow\,$
$\displaystyle\dfrac{1}{n}\sum_{i=1}^{n}\lim_{\epsilon\downarrow
0}\dfrac{\|\bm{Y}(t_{i})-\widehat{\bm{q}}+\epsilon\bm{v}\|_{\bm{K}(t_{i})}-\|\bm{Y}(t_{i})-\widehat{\bm{q}}\|_{\bm{K}(t_{i})}}{\epsilon}+\lim_{\epsilon\downarrow
0}\dfrac{\langle\bm{u},\epsilon\bm{v}\rangle_{\bm{K}(t_{i})}}{\epsilon}\geqslant
0$ $\displaystyle\Rightarrow\,$
$\displaystyle\dfrac{1}{n}\sum_{i=1}^{n}\left[\dfrac{\langle\bm{Y}(t_{i})-\widehat{\bm{q}},\bm{v}\rangle_{\bm{K}(t_{i})}}{\|\bm{Y}(t_{i})-\widehat{\bm{q}}\|_{\bm{K}(t_{i})}}+\langle\bm{u},\bm{v}\rangle_{\bm{K}(t_{i})}\right]\geqslant
0.$

Using the same steps with $(-\bm{v})$ in place of $\bm{v}$ we arrive (3.5).
Similarly, (3.6) stems from (2.4). $\square$

Proof of Theorem 4.1: Defining
$\widehat{\bm{q}}=\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$ and
$\bm{q}_{0}=Q(\bm{u},\bm{x})$, we write the chain of inequalities:

$\displaystyle\mathcal{M}_{\bm{u}}^{(p)}(\widehat{\bm{q}})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}_{0})$
$\displaystyle\leqslant\mathcal{M}_{\bm{u}}^{(p)}(\widehat{\bm{q}})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}_{0})-\mathcal{M}_{\bm{u},n}^{(p)}(\widehat{\bm{q}})+\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})$
$\displaystyle\leqslant
2\sup_{\bm{q}\in\mathcal{Q}_{p}}|\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})|\,,$

where $\mathcal{Q}_{p}\subseteq\mathbb{R}^{p}$ is a set containing both
$\bm{q}_{0}$ and the sequence of estimators
$\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$. Hence,

$\displaystyle\mathbb{P}\left(\|\widehat{\bm{q}}-\bm{q}_{0}\|>Ap^{2}/n\right)$
$\displaystyle\leqslant$
$\displaystyle\mathbb{P}\left(\|\widehat{\bm{q}}-\bm{q}_{0}\|>Ap^{2}/n,\,\mathcal{M}_{\bm{u}}^{(p)}(\widehat{\bm{q}})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}_{0})<Bp^{2}/n\right)$
(6.2) $\displaystyle\mbox{\hskip
128.0374pt}+\mathbb{P}\left(\mathcal{M}_{\bm{u}}^{(p)}(\widehat{\bm{q}})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}_{0})\geqslant
Bp^{2}/n\right)$ $\displaystyle\leqslant$
$\displaystyle\mathbb{P}\Bigl{(}\|\widehat{\bm{q}}-\bm{q}_{0}\|>Ap^{2}/n,\,\mathcal{M}_{\bm{u}}^{(p)}(\widehat{\bm{q}})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q}_{0})<Bp^{2}/n\Bigr{)}$
$\displaystyle\mbox{\hskip
128.0374pt}+\mathbb{P}\biggl{(}\sup_{\bm{q}\in\mathcal{Q}_{p}}|\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})|\geqslant
Bp^{2}/2n\biggr{)}.$

Owing to the convexity lemma by Pollard, (1991), the objective function
$\mathcal{M}_{\bm{u}}$ as part of (2.4) is convex on $\mathcal{Q}_{p}$.
Consequently, the sequence of conditional regression-quantile estimators
$\hat{q}_{n}$ along $\bm{u}$ is well-defined and the true quantile $q_{0}$ is
a well-separated minimizer of $\mathcal{M}_{\bm{u}}$, that is, there exists a
$q_{0}\in\mathcal{U}_{p}$ for which
$\inf_{\bm{q}:\,\|\bm{q}-\bm{q}_{0}\|\geqslant\delta}\,\left(\mathcal{M}_{\bm{u}}(\bm{q})-\mathcal{M}_{\bm{u}}(\bm{q}_{0})\right)>0$,
for arbitrarily small $\delta>0$. This ensures the boundedness of the first
term in the above sum. In particular, since $p=o(n^{1/2})$, taking
$\delta=p^{2}/n$ ensures that the difference
$\left(\mathcal{M}_{\bm{u}}(\bm{q})-\mathcal{M}_{\bm{u}}(\bm{q}_{0})\right)$
is strictly positive. It is then possible to fix $B=A\epsilon$ with
arbitrarily small $\epsilon>0$, which makes the first term equal to $0$. This
takes care of (6.2). Consequently, the gist of the proof is in verifying the
uniform convergence of (6.2). From the sub-Gaussian property of $\bm{Y}(t)$
together with the law of large numbers, it results that

$\left|\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})-\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\right]\right|=\mathcal{O}_{p}(\sqrt{p/n}).$
(6.3)

In order to connect the weighted expectation affecting the kernel $\bm{K}(t)$
with $\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})$, we invoke the law of iterated
expectations upon the conditioning random variable
$R=\left\|\bm{X}(t)-\bm{x}\right\|$:

$\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\right]=\mathbb{E}\left[\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\mid
R\right]\right].$

There is a one-to-one mapping between $R$ and the continuous measurable
functional that is the kernel $\bm{K}(t)$. In particular, the event $R=0$,
amounts to $\bm{K}(t)=\bm{I}_{p}$ due to $\mathcal{K}(0)=c_{n}\bm{I}_{p}$.
Since there exists $\delta_{R}>0$ such that for every $R<\delta_{R}$, all
eigenvalues of the matrix $(\bm{K}(t)-\bm{I}_{p})$ are bounded by
$\sqrt{p/n}$, thus rendering valid the chain of inequalities:

$\displaystyle\left|\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\right]-\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle\mid\bm{X}(t)=\bm{x}\right]\right|$
$\displaystyle=$
$\displaystyle\left|\mathbb{E}\left[\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\mid
R\right]\right]-\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle\mid
R=0\right]\right|$ $\displaystyle=$
$\displaystyle\left|\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\mid
R\leqslant\delta_{R}\right]\mathbb{P}(R\leqslant\delta_{R})\right.$
$\displaystyle\qquad\left.+\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\mid
R>\delta_{R}\right]\mathbb{P}(R>\delta_{R})\right.$
$\displaystyle\qquad\left.-\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle\mid
R=0\right]\right|$ $\displaystyle\leqslant$
$\displaystyle\left|\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}-\|\bm{Y}(t)-\bm{q}\|\mid
R\leqslant\delta_{R}\right]\right|$
$\displaystyle\qquad+\left|\mathbb{E}\left[\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}-\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle\mid
R\leqslant\delta_{R}\right]\right|$
$\displaystyle\qquad+\left|\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\mid
R>\delta_{R}\right]\right|$

The first term in the latter is bounded by the product of the absolute value
of the maximum eigenvalue of $(\bm{K}(t)-\bm{I}_{p})$ and the norm of
$\|\bm{Y}(t)-\bm{q}\|$. When $R\leqslant\delta_{R}$, combining with the sub-
Gaussian behavior of $\bm{Y}(t)$, the first term can be bounded by
$\mathcal{O}(p/\sqrt{n})$. For the second term, an application of the Cauchy-
Schwartz inequality combined with the fact that $\|\bm{u}\|\leqslant 1$ yields
that it is also bounded by $\mathcal{O}(p/\sqrt{n})$. For the third term, note
that as the bandwidth $b_{n}\rightarrow 0$, for sufficiently large $n$,
$\delta_{R}>b_{n}$. Since the kernel $\mathcal{K}$ has bounded support in the
unit hypersphere $B^{d-1}$, it follows that $\bm{K}(t)$ can be made equal to
the null matrix $\bm{0}$ in the third term. Therefore,

$\left|\mathbb{E}\left[\|\bm{Y}(t)-\bm{q}\|_{\bm{K}(t)}+\langle\bm{u},\bm{Y}(t)-\bm{q}\rangle_{\bm{K}(t)}\right]-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})\right|=\mathcal{O}(p/\sqrt{n})$

Combining this with (6.3), we obtain that for each fixed
$\bm{q}\in\mathbb{R}^{p}$,

$\|\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})-\mathcal{M}_{\bm{u}}^{(p)}(\bm{q})\|=\mathcal{O}_{p}(\max\\{p^{2}/n,p/\sqrt{n}\\})=\mathcal{O}_{p}(p/\sqrt{n}),$

where the last equality follows from the initial assumption that
$p=o(n^{1/2})$. The uniform convergence holds if we can show that
$\mathcal{Q}_{p}$ is a compact set. Thus, it is enough to show that the
sequence of estimators $\bm{q}_{n}$ minimizing the sequence of objective
functions $\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ are uniformly bounded.
Without the loss of generality, consider that $0$ belongs to the convex hull
of the datapoints $\\{\bm{Y}(t_{i}):i=1,2,\dots n\\}$, otherwise we can
translate all datapoints suitably. By applying the sub-Gaussian property of
$\bm{Y}(t_{i})$ and union bound, it follows that there exists some $c>0$ such
that

$\mathbb{P}\left(\sup_{1\leqslant i\leqslant
n}\|\bm{Y}(t_{i})\|>M\right)\leqslant ne^{-cp}.$

As $p^{-1}\log(n)=\mathcal{O}(1)$, by choosing sufficiently large $p\geqslant
C\log(n)$ for some absolute constant $C$, we can now make
$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{0})<M$ with probability arbitrarily close to
$1$, by bounding each $\bm{Y}(t_{i})$. Also, since $\|\bm{u}\|<1$, there
exists a $\delta\in(0,1)$ such that $\|\bm{u}\|<(1-\delta)$, and a vector
$\bm{v}$ such that $\|\bm{v}\|>2M\delta^{-1}$. Consider the objective function
at a point $\bm{Y}(t_{j})+\bm{v}$ for any $j=1,2,\dots n$;

$\displaystyle\mathcal{M}_{\bm{u},n}^{(p)}(\bm{Y}(t_{j})+\bm{v})$
$\displaystyle=\dfrac{1}{n}\sum_{i=1}^{n}\left[\|\bm{Y}(t_{i})-\bm{Y}(t_{j})-\bm{v}\|_{\mathcal{K}(t_{i})}+\langle\bm{u},\bm{Y}(t_{i})-\bm{Y}(t_{j})-\bm{v}\rangle_{\mathcal{K}(t_{i})}\right]$
$\displaystyle\geqslant\|\bm{v}\|/2-(1-\delta)\|\bm{v}\|/2,\mbox{ (by Cauchy-
Schwartz inequality)}$ $\displaystyle\geqslant
M\geqslant\mathcal{M}_{\bm{u},n}^{(p)}(\bm{0})$

By general positioning of these vectors,
$\\{\bm{Y}(t_{j})+\bm{v},\bm{Y}(t_{j})-\bm{v}:j=1,\dots n\\}$ do not lie in a
subspace with a dimension strictly less than $p$. Combining this with the
convex nature of the objective function
$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$, it follows that along the boundary of
the convex hull of these points also, the value of the objective function is
more than $M$. Using continuity of the objective function, we find that there
is a $(p-1)$-dimensional hypersphere, independent of the choice of $n$,
containing the origin $\bm{0}$ such that the
$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})>\mathcal{M}_{\bm{u},n}^{(p)}(\bm{0})$
for any $\bm{q}$ lying on the boundary of that hypersphere for any $n$. Now it
follows from the strict convexity of each
$\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ that they have a local minimum within
that hypersphere. Choosing the $\widehat{\bm{q}}_{n}$ as those local minima
establishes the boundedness of them. $\square$

Proof of Theorem 4.2: We note that under the conditions of the theorem,
Propositions (P1)-(P3) are valid. We start by employing Taylor’s expansion of
$\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})$ about the true parameter
$\bm{q}_{0}:=Q(\bm{u},\bm{x})$ while evaluated at
$\widehat{\bm{q}}:=\widehat{\bm{q}}_{n}(\bm{u},\bm{x})$, to obtain:

$\displaystyle\bm{0}$
$\displaystyle=\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\widehat{\bm{q}})$
$\displaystyle=\nabla\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})+\nabla^{2}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})(\widehat{\bm{q}}-\bm{q}_{0})+\mathcal{O}\Bigl{(}\|\widehat{\bm{q}}-\bm{q}_{0}\|^{2}\left\|\nabla^{3}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q})\right\|\Bigr{)}$

Denoting $\nabla^{i}\mathcal{M}_{\bm{u},n}^{(p)}(\bm{q}_{0})=\Delta_{i}$, we
have, with $\delta_{n}=p/\sqrt{n}$, that

$a_{n}\left(\widehat{\bm{q}}-\bm{q}_{0}\right)=-a_{n}\Delta_{2}^{-1}\Delta_{1}+\mathcal{O}(\|\Delta_{2}^{-1}\|\delta_{n}^{2}/pb_{n}).$

Here, we use the consistency result of $\widehat{\bm{q}}$ coupled with
proposition (P3) to bound $\|\Delta_{3}\|$. By proposition (P2), it follows
that for sufficiently large $n$,

$\left\|\Delta_{2}^{-1}-\bm{\Psi}(\bm{q}_{0})^{-1}\right\|\leqslant\dfrac{2}{\|\bm{\Psi}(\bm{q}_{0})\|^{2}}\left\|\bm{\Psi}(\bm{q}_{0})-\Delta_{2}\right\|\leqslant\dfrac{Cp^{2}}{n\|\bm{\Psi}(\bm{q}_{0})\|^{2}}.$

Therefore,

$\displaystyle\quad\left|\mathbb{E}\left(h\left(a_{n}\left(\widehat{\bm{q}}-\bm{q}_{0}\right)\right)\right)-\mathbb{E}\left(h\left(\bm{\Psi}^{-1}(\bm{q}_{0})Z(\bm{q}_{0})\right)\right)\right|$
$\displaystyle\leqslant$
$\displaystyle\quad\left|\mathbb{E}\left(h\left(-a_{n}\Delta_{2}^{-1}\Delta_{1}+\mathcal{O}(\|\Delta_{2}^{-1}\|p/nb_{n})\right)\right)-\mathbb{E}\left(h\left(\bm{\Psi}^{-1}(\bm{q}_{0})Z(\bm{q}_{0})\right)\right)\right|$
$\displaystyle\leqslant$
$\displaystyle\quad\left|\mathbb{E}\left(h\left(-a_{n}\Delta_{2}^{-1}\Delta_{1}+\mathcal{O}(\|\Delta_{2}^{-1}\|p/nb_{n})\right)\right)-\mathbb{E}\left(h\left(-a_{n}\bm{\Psi}^{-1}(\bm{q}_{0})\Delta_{1}\right)\right)\right|+$
$\displaystyle\qquad\qquad\left|\mathbb{E}\left(h\left(-a_{n}\bm{\Psi}^{-1}(\bm{q}_{0})\Delta_{1}\right)\right)-\mathbb{E}\left(h\left(\bm{\Psi}^{-1}(\bm{q}_{0})Z(\bm{q}_{0})\right)\right)\right|$
$\displaystyle\leqslant$
$\displaystyle\quad\mathbb{E}\left[\left\|h^{\prime}(\widetilde{\bm{q}})\right\|\left(\mathcal{O}\left(\|\Delta_{2}^{-1}\|\frac{p}{nb_{n}}\right)+|a_{n}|\|\Delta_{1}\|\left\|\Delta_{2}^{-1}-\bm{\Psi}(\bm{q}_{0})^{-1}\right\|\right)\right]+\mathcal{O}(p^{2}/n)\,,$

where $\widetilde{\bm{q}}$ is some vector lying between
$-a_{n}\bm{\Psi}^{-1}(\bm{q}_{0})\Delta_{1}$ and
$-a_{n}\Delta_{2}^{-1}\Delta_{1}$, ensured by the mean value theorem. Here, we
use the fact that the function $g(x):=\bm{\Psi}^{-1}(\bm{q}_{0})f(x)$ is also
an element of $\mathcal{H}$, hence
$\left|\mathbb{E}\left(g\left(-a_{n}\Delta_{1}\right)\right)-\mathbb{E}\left(g\left(Z(\bm{q}_{0})\right)\right)\right|$
is $\mathcal{O}(p^{2}/n)$ by virtue of (P1). On the other hand, due to the
choice of $\mathcal{H}$ and $h\in\mathcal{H}$, the quantity
$\left\|h^{\prime}(\widetilde{\bm{q}})\right\|$ remains bounded. Noting that
$\left\|\bm{\Psi}(\bm{q}_{0})\right\|$ is positive definite and has bounded
entries, its eigenvalues are of order $\mathcal{O}(\sqrt{p})$. This implies
that the first term is at most $\mathcal{O}(p^{17/8}/n)$, whence

$\left|\mathbb{E}\left(h\left(a_{n}\left(\widehat{\bm{q}}-\bm{q}_{0}\right)\right)\right)-\mathbb{E}\left(h\left(\bm{\Psi}^{-1}(\bm{q}_{0})Z(\bm{q}_{0})\right)\right)\right|=\mathcal{O}(p^{17/8}/n).$

Since this bound does not depend on any given $h\in\mathcal{H}$, the above
conclusion holds with the supremum over $h\in\mathcal{H}$. Finally, the choice
of $\bm{u}$ is restricted within the compact set $\mathcal{B}_{p}$, giving way
to a further supremum on this set. $\square$

##  7 Forecasting energy demand

In order to experimentally verify the good overall performance of the proposed
quantile estimator, we conduct a simulation study that encompasses two
distinct setups: one without any explicit spatial information and another
setup with proper spatial information included. In both cases, we find that
the proposed estimator achieves reasonable prediction accuracy. The simulation
results are compiled in Appendix B featuring in the supplementary material.

In this section, the proposed quantile regression methodology is applied to
produce inferences on the Thames Valley Vision (TVV) data, collected by
Scottish and Southern Electricity Network (SSEN), the primary industrial
partner of the TVV project funded by UK gas and electricity regulator Ofgem
through the Low Carbon Networks Fund and Network Innovation Competition. The
dataset handled here has been partly studied in Jacob et al., (2020). The
project’s overall aim was to evaluate a typical low voltage network using
monitoring in households and substations in order to simulate future demand
scenarios.

The dataset under study consists of the energy profiles for $226$ households
with half-hourly resolution, between 19th July 2015 and 20th September 2015,
thus spanning $9$ weeks of data. In our analysis, the first $8$ weeks of data
are used as a training set whereas our proposed estimator is employed to
predict different quantiles for last week’s electricity demand. Given the
temporal resolution of the data, the number of observations for every location
in the training set is $n=2688$. In line with the theory underpinning the
quantile estimation introduced in the present paper (see Theorem 4.1), the
number of spatial points $p=226$ far exceeds the benchmark of $\sqrt{n}\approx
52$. This provides opportunity to split the households into clusters of size
not larger than $52$ based on their average daily demand (ADD) profiles. To
determine these clusters, we calculated the pairwise dynamic time warping
(DTW) distance between every two household’s ADD profiles using the FastDTW
algorithm from Gold and Sharir, (2018) with Hierarchical clustering on the
households. Figure 1 displays the heatmap for the pairwise DTW distances (log-
transformed) between households. Through this procedure, the households were
assigned to clusters $j=1,2,\ldots,12$ of sizes 21, 52, 42, 17, 10, 42, 28, 4,
2, 4, 2 and 2, respectively.

Figure 1: Heatmap of DTW distances (transformed to log-scale) between average
daily household profiles, arranged according to their assigned clusters.

We find that the pairwise within-cluster distances tend to be smaller than the
between-cluster distances, which is reflected in the lighter shaded blocks
containing the reverse diagonal of the heatmap. Hence, it is reasonable to
assume that electricity demand profiles of two households residing in two
different clusters are independent or weakly dependent (the latter is not a
deterrent to our methods). Consequently, we applied our proposed quantile
estimation to each cluster separately. In a nutshell, our statistical analysis
revolves broadly around two goals:

  * (i)

the estimation of relevant quantiles for each cluster of households and
prediction of the quantiles for profiles in the $9$-th week across all $226$
households. As covariates, we take day of the week and hour of the day, as
well as a trend over time (not necessarily monotonic).

  * (ii)

the statistical evaluation of homogeneity of the effects of the day and hour
covariates for different households. In addition, discerning significance of a
potential trend over time in the quantiles by relying on the methodology
described in Section 5.

Throughout this analysis, for ease of implementation and interpretation, we
fix the levels of 0.5, 0.75, 0.8, 0.9, 0.95, 0.99 across all the cases
examined.

###  7.1 One-week-ahead quantile estimation

For tackling the estimation problem described in (i), we estimated quantile
functions from the training data relating the electricity demand profiles per
cluster. The one-week ahead predicted quantiles are displayed in Figure 2. The
plots seem to indicate that the estimated quantiles enjoy a smoothness
property due to the continuous and smooth nature of the kernel function
adopted in the estimator. There is no apparent quantile crossing issue as
well. Unsurprisingly, two quantile curves draw closer when the general trend
in demand is approximately constant, but on certain time of the day quantile
curves associated with the larger probabilities shoot away from the median
level, in what can be identified as peak hours. The systematic cyclical
pattern inherent to electricity demand due to the daily day-night cycle is
also captured in the prediction. In particular, clusters 2, 4, 8 and 12
indicate high electricity loads during the weekends. While for most clusters,
the quantile curves provide a satisfactory description of how a proper
probability distribution would spread, the quantile curves for cluster 9 crop
up very tightly together. This is mainly because there are only two households
(with numbers 106 207) in this cluster, making the number of spatial points
clearly insufficient for inference on the true range of the quantiles and not
the least of their temporal evolution. On the same note, the pattern of the
cyclical nature of the electricity demand for Cluster 12 is in stark contrast
with other clusters, where peak values tend to occur at the beginning of the
day. The households in this cluster are 202 and 210.

Figure 2: Predicted quantile curves for electricity demand during the 9th week
out of 9 weeks of daily profile data, for cluster $j=1,2,\ldots,12$ of size
21, 52, 42, 17, 10, 42, 28, 4, 2, 4, 2 and 2.

###  7.2 Hypotheses testing

Before getting underway, it is important to caution the reader that the two
testing procedures illustrated in this section are computationally
challenging. Since computational efficiency is not the main focus of this
paper, we flag this up as a pertinent aspect for further improvement and limit
the present study to two clusters. The selected clusters 2 and 12 to this
undertaking are the two extreme cases in terms of size, the second cluster
enclosing the maximum allowed of 52 households while the 12th cluster is made
up of only two. Table 1 contains the $p$-values of the two tests of
homogeneity introduced in Section 5 conducted for several quantile levels,
$\tau=0.5,0.75,0.8,0.9,0.95,0.99$.

Table 1: Computed $p$-values for the results of the two tests for assessing homogeneity. Test | Quantile | Cluster 2 | Cluster 12
---|---|---|---
Homogeneity of | 0.5 | $8.07\times 10^{-6}$ | $1.05\times 10^{-6}$
effect of covariates | 0.75 | $6.45\times 10^{-6}$ | $1.04\times 10^{-6}$
| 0.9 | $8.10\times 10^{-6}$ | $1.05\times 10^{-6}$
| 0.95 | $7.96\times 10^{-6}$ | $1.05\times 10^{-6}$
| 0.99 | $8.07\times 10^{-6}$ | $1.05\times 10^{-6}$
Homogeneity of | 0.5 | $2.79\times 10^{-7}$ | $1.88\times 10^{-8}$
effect of time | 0.75 | $2.03\times 10^{-7}$ | $1.88\times 10^{-8}$
| 0.9 | $2.78\times 10^{-7}$ | $1.88\times 10^{-8}$
| 0.95 | $2.79\times 10^{-7}$ | $1.88\times 10^{-8}$
| 0.99 | $2.79\times 10^{-7}$ | $1.88\times 10^{-8}$

It is clear that all $p$-values yielded are highly significant, leading to the
rejection of the null hypothesis of homogeneity in all cases. Interestingly,
this finding is unchallenged by the size the clusters and is also consistent
across quantile levels. Therefore, the is sufficient evidence to conclude that
the covariates have significantly different effects for different household’s
demands. There is also significant evidence of a temporal trend in the
quantiles of all levels considered. Given the magnitude of the $p$-values,
these conclusions hold for any nominal level of test, i.e., both null
hypotheses are rejected notably at the minimum of the usual significance
levels ($\alpha=1\%$).

## Acknowledgements

Cláudia Neves gratefully acknowledges support from EPSRC-UKRI Innovation
Fellowship grant EP/S001263/2. Her work is also partly supported by CEAUL,
Faculty of Science, University of Lisbon, DOI: 10.54499/UIDB/00006/2020,
https://doi.org/10.54499/UIDB/00006/2020.

##  8 Supplementary documentation

The supplementary material is comprised of two sections, namely: Appendix A
with the proof for Proposition 1 in full; Appendix B containing simulation
results based on a comprehensive set of numerical experiments.

## References

  * Barthelmé and Usevich, (2021) Barthelmé, S. and Usevich, K. (2021).  Spectral properties of kernel matrices in the flat limit.  SIAM Journal on Matrix Analysis and Applications, 42(1):17–57.
  * Belloni and Oliveira, (2018) Belloni, A. and Oliveira, R. I. (2018).  A high dimensional central limit theorem for martingales, with applications to context tree models.  arXiv preprint arXiv:1809.02741.
  * Chaudhuri, (1996) Chaudhuri, P. (1996).  On a geometric notion of quantiles for multivariate data.  Journal of the American Statistical Association, 91(434):862–872.
  * Chaudhuri et al., (1997) Chaudhuri, P., Doksum, K., and Samarov, A. (1997).  On average derivative quantile regression.  The Annals of Statistics, 25(2):715 – 744.
  * Chen, (1996) Chen, Z. (1996).  Conditional lp-quantiles and their application to the testing of symmetry in non-parametric regression.  Statistics & Probability Letters, 29(2):107–115.
  * Chowdhury and Chaudhuri, (2019) Chowdhury, J. and Chaudhuri, P. (2019).  Nonparametric depth and quantile regression for functional data.  Bernoulli, 25(1):395 – 423.
  * Das and Lahiri, (2020) Das, D. and Lahiri, S. (2020).  Central limit theorem in high dimensions : The optimal bound on dimension growth rate.
  * Das and Ghosal, (2017) Das, P. and Ghosal, S. (2017).  Analyzing ozone concentration by bayesian spatio-temporal quantile regression.  Environmetrics, 28(4):e2443.
  * De Gooijer et al., (2006) De Gooijer, J. G., Gannoun, A., and Zerom, D. (2006).  A multivariate quantile predictor.  Communications in Statistics—Theory and Methods, 35(1):133–147.
  * Deb and Jana, (2021) Deb, S. and Jana, K. (2021).  Nonparametric quantile regression for time series with replicated observations and its application to climate data.  arXiv preprint arXiv:2107.02091.
  * Deb et al., (2017) Deb, S., Pourahmadi, M., and Wu, W. B. (2017).  An asymptotic theory for spectral analysis of random fields.  Electronic Journal of Statistics, 11:4297–4322.
  * Duan et al., (2021) Duan, Q., McGrory, C. A., Brown, G., Mengersen, K., and Wang, Y.-G. (2021).  Spatio-temporal quantile regression analysis revealing more nuanced patterns of climate change: a study of long-term daily temperature in australia.  arXiv preprint arXiv:2103.05791.
  * Eddy, (1985) Eddy, W. F. (1985).  Ordering of multivariate data.  Computer Science and Statistics: The Interface. L. Billard, Amesterdam: North-Holland, pages 25–30.
  * Einbeck and Tutz, (2006) Einbeck, J. and Tutz, G. (2006).  Modelling beyond regression functions: an application of multimodal regression to speed–flow data.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 55(4):461–475.
  * Einmahl and Mason, (1992) Einmahl, J. H. and Mason, D. M. (1992).  Generalized quantile processes.  The Annals of Statistics, pages 1062–1078.
  * Ezzahrioui and Ould-Said, (2008) Ezzahrioui, M. and Ould-Said, E. (2008).  Asymptotic normality of the kernel estimators of the conditional quantile in the normed space.  Far East J. Theoretical Statist, 25:15–38.
  * Ezzahrioui and Ould-Saïd, (2008) Ezzahrioui, M. and Ould-Saïd, E. (2008).  Asymptotic results of a nonparametric conditional quantile estimator for functional time series.  Communications in Statistics—Theory and Methods, 37(17):2735–2759.
  * Fan et al., (1996) Fan, J., Gijbels, I., Hu, T.-C., and Huang, L.-S. (1996).  A study of variable bandwidth selection for local polynomial regression.  Statistica Sinica, 6(1):113–127.
  * Ferraty et al., (2005) Ferraty, F., Rabhi, A., and Vieu, P. (2005).  Conditional quantiles for dependent functional data with application to the climatic" el niño" phenomenon.  Sankhyā: The Indian Journal of Statistics, pages 378–398.
  * Figalli, (2018) Figalli, A. (2018).  On the continuity of center-outward distribution and quantile functions.  Nonlinear Analysis, 177:413–421.
  * Gold and Sharir, (2018) Gold, O. and Sharir, M. (2018).  Dynamic time warping and geometric edit distance: Breaking the quadratic barrier.  ACM Transactions on Algorithms (TALG), 14(4):1–17.
  * Huber, (2011) Huber, P. J. (2011).  Robust statistics.  In International encyclopedia of statistical science, pages 1248–1251. Springer.
  * Jacob et al., (2020) Jacob, M., Neves, C., and Vukadinović Greetham, D. (2020).  Forecasting and assessing risk of individual electricity peaks.  Springer Nature.
  * Koenker and Bassett Jr, (1978) Koenker, R. and Bassett Jr, G. (1978).  Regression quantiles.  Econometrica: journal of the Econometric Society, pages 33–50.
  * Konen and Paindaveine, (2022) Konen, D. and Paindaveine, D. (2022).  Multivariate $\rho$-quantiles: A spatial approach.  Bernoulli, 28(3):1912 – 1934.
  * Kosorok, (2008) Kosorok, M. R. (2008).  Introduction to Empirical Processes and Semiparametric Inference.  Springer.
  * van der Vaart, (1998) van der Vaart, A. W. (1998).  Asymptotic Statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Pollard, (1991) Pollard, D. (1991).  Asymptotics for least absolute deviation regression estimators.  Econometric Theory, 7(2):186–199.
  * Portnoy, (1986) Portnoy, S. (1986).  On the central limit theorem in rp when p→ $\infty$.  Probability Theory and Related Fields, 73(4):571–583.
  * Reich, (2012) Reich, B. J. (2012).  Spatiotemporal quantile regression for detecting distributional changes in environmental processes.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 61(4):535–553.
  * Reich et al., (2011) Reich, B. J., Fuentes, M., and Dunson, D. B. (2011).  Bayesian spatial quantile regression.  Journal of the American Statistical Association, 106(493):6–20.  PMID: 23459794.
  * Resnick, (1987) Resnick, S. I. (1987).  Extreme Values, Regular Variation and Point Processes.  Springer.
  * Takeuchi et al., (2006) Takeuchi, I., Le, Q., Sears, T., Smola, A., et al. (2006).  Nonparametric quantile estimation.
  * Vershynin, (2010) Vershynin, R. (2010).  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027.
  * Zhao and Wu, (2008) Zhao, Z. and Wu, W. B. (2008).  Confidence bands in nonparametric time series regression.  The Annals of Statistics, 36(4):1854–1878.
  * Zhou, (2010) Zhou, Z. (2010).  Nonparametric inference of quantile curves for nonstationary time series.  The Annals of Statistics, 38(4):2187.
  * Zhou and Wu, (2009) Zhou, Z. and Wu, W. B. (2009).  Local linear quantile estimation for nonstationary time series.  The Annals of Statistics, 37(5B):2696 – 2729.