Andrew Mirii Kinyita¹, Anthony Waititu², Antony Ngunyi^1
1Department of Statistics and Actuarial Science, Dedan Kimathi University of Technology, Nyeri, Kenya.
²Department of Statistics and Actuarial Science, Jomo Kenyatta University of Agriculture and Technology, Nairobi, Kenya.
DOI: 10.4236/ojs.2026.165021   PDF    HTML   XML   1 Downloads   33 Views   Citations

Abstract

This article presents an analytical framework that utilizes deep neural network-based quantile regression to detect structural changes in high-dimensional, time dependent data. By leveraging on deep neural networks for quantile regression, the model combines the robustness of quantile regression with the computational efficiency of deep learning to estimate conditional quantiles. These estimators can identify abrupt shifts in data distribution without imposing strict structural assumptions, making them more resilient to heavy-tailed data compared to mean-based methods. The study derives the asymptotic properties of the deep quantile regression estimator for change point detection, particularly within high-dimensional, non-linear, and complex datasets. Under specific regularity conditions, the asymptotic properties of the proposed estimator are established. The simulation results demonstrate the convergence of the deep quantile regression estimator to the true conditional quantile function as the sample size. The deep quantile regression estimator achieves the minimax optimal convergence rate, effectively alleviating the curse of dimensionality when the target function adheres to a compositional structure. The empirical process theory and the method of sieve estimation are employed to demonstrate the consistency of the deep quantile regression estimator. The results show that the normalized estimator is asymptotically normal and exhibits inherent robustness to outliers and heavy-tailed noise compared to least-squares-based deep learning models. This theoretical framework provides a foundation for the construction of confidence intervals and hypothesis testing for complex, non-parametric quantile models in fields such as finance, risk management, and econometrics.

Keywords

Change Point, Single Change Point, Multiple Change Points, Deep Quantile Regression Estimator, Deep Neural Network

Share and Cite:

1. Introduction

Change happens all the time. So it is natural to become interested in detecting and estimating changes. The statistical literature on change-point analysis dates back to the 1950s [1], when it was originally used to detect changes in the mean values of independently and identically distributed Gaussian variables for quality control in the industry. As changes are not a problem of a specific scientific field, the theory of change point methods has been extensively studied across various fields, including industrial manufacturing, healthcare [2], financial data analysis [3] and environmental monitoring [4].

Over the decades, many statistical model-based methods for change point detection (CPD) have been developed to address the needs of various applications. Many publications analysing different models or adopting different points of view have since followed. Besides testing for a change, the estimation of the time of the change is also of interest. [5] considered the estimation of the change-point for the mean change model based on independent and identically distributed (i.i.d.) random variables. Most of the existing statistical models developed require some form of distributive pattern in the data and estimate the efficiency of a detected change through the theoretical asymptotic presentation [6]. Among the various statistical methods is the Bayesian change-point model, which integrates prior knowledge of changes within a probabilistic framework. This approach facilitates more informed decision-making by allowing for the systematic analysis of data where change points are anticipated [7].

Quantile regression, introduced by [8], is a vital econometric tool for modelling the conditional distribution of a response variable. Unlike traditional least squares methods that focus on the conditional mean, quantile regression estimates various points (quantiles) in the distribution, providing a deeper insight into the relationship between variables. Furthermore, quantile regression serves as a robust alternative to traditional mean-based forecasting methods, as it directly estimates conditional quantiles, thereby providing a more nuanced understanding of potential outcomes. Despite these advancements, a significant challenge remains: ensuring the consistency and reliability of change point estimators within these sophisticated high-dimensional models. This challenge underscores the need for continued research and development in this area to enhance predictive accuracy and risk management strategies. Recent developments have advanced quantile regression into more complex settings. High-dimensional quantile regression utilises regularisation techniques to manage large sets of predictors [9] [10]. Additionally, semi-parametric and non-parametric quantile models provide flexibility in functional forms [11]. However, these approaches often depend on assumptions of sparsity or smoothness and may still have difficulty capturing highly non-linear interactions in financial data.

Deep Neural Networks (DNNs) have emerged as powerful non-parametric function approximators capable of modelling highly non-linear relationships. Universal approximation theorems establish that feed-forward neural networks with sufficient width can approximate any continuous function on compact domains [12]. Later work by [13] and others refined these results by characterising approximation error rates in terms of network depth, width, and smoothness of the target function. Deep Neural Networks (DNNs) are popular for their ability to approximate complex nonlinear functions. According to the universal approximation theorem, feed-forward neural networks can approximate any continuous function under mild conditions. This makes DNNs well-suited for modelling the non-linear and high-dimensional dependencies often found in financial time series.

In recent years, advancements in deep learning methodologies have shown remarkably improved results in several time series problems, such as modelling and forecasting of time series data, allowing for the identification of complex non-linear dependencies and interactions among diverse risk factors. This capability enables practitioners to gain deeper insights into market dynamics. Deep learning models designed especially for the processing of sequential data, including but not limited to Recurrent Neural Networks (RNNs) and Auto-encoder models, are being increasingly applied to overcome certain difficulties presented in high-dimensional data sets [14]. Although these deep learning models do not have such a strong theoretical basis compared with conventional statistical models, they demonstrate great capability in better handling complex change point detection tasks and in deeply excavating data features and transmitting temporal information within time series data. By doing so, this also provides them with an extended capability of finding out patterns and dependencies that may not be considered in traditional statistical methods. Thus, this model further develops the detecting ability of significant shifts that are difficult for statistical models to detect with ease, thus giving higher accuracy in practical applications without explicit model assumptions. Moreover, rapid development in hardware technology, such as parallelisation in GPUs and CPUs, has also ensured that the computational requirements for deep learning methods remain within feasible limits for practical applications.

Recent advances in deep learning have led to the development of various methods for change-point detection (CPD). These methods have enhanced the ability of CPD techniques to handle complex, high-dimensional data, making them more adaptable and less dependent on strict assumptions about the distributions of data. With increased complexity and volume of data, however, the development of robust and scalable methods that could effectively detect change points has increased. Despite their theoretical robustness, applications to high-dimensional time series data have mostly demonstrated limitations of classical methods. These usually manifest as model complexity and generalisation challenges across diverse, data-intensive environments [15].

The concept of deep quantile estimator that utilizes neural networks to predict Value-at-Risk (VaR) by modelling non-linear relationships between asset returns and predictors was introduced by [16]. This methodology capitalizes on the universal approximation property of neural networks to effectively capture general, non-linear associations for conditional quantiles, thereby advancing beyond traditional linear quantile regression techniques. The framework is adaptable, accommodating various penalty functions and high-dimensional covariates. Monte Carlo simulations reveal that the deep quantile estimator demonstrates commendable finite-sample performance. Moreover, empirical analyses indicate that it yields substantial improvements over linear quantile regression and other forecasting models for VaR, as substantiated by a range of statistical tests.

While deep learning methods have shown impressive empirical performance, their theoretical properties are still under active development. Recent contributions include [17], who provides convergence rates for neural network estimators; [18] established the asymptotic theory for deep learning in econometrics; and [19] analysed the neural networks as sieve estimators.

The consistency of deep quantile regression estimators under assumptions on network complexity, smoothness of the conditional quantile function, and uniform convergence of the empirical quantile loss was established by [20], [21] derived oracle inequalities for neural network estimators under general convex but non-smooth loss functions, including the quantile loss. These studies lay important foundations for understanding the statistical properties of deep quantile regression models. Theoretical work in this area remains relatively sparse.

The rest of the paper is organised as follows. Section 2 introduces the methodology for estimating multiple change points within deep quantile regression settings. Section 2.3.3 presents the statistical theory for quantile regression with ReLU networks and theoretical properties of the proposed deep quantile regression estimator. Section 3 presents the simulation study. Finally, Section 4 concludes the study.

2. Methodology

In this section, the underlying theory of a quantile regression is introduced for both linear quantile and non-linear quantile regression model. This section presents some fundamental theoretical results on the universal approximation property of neural networks that can be used to approximate a non-linear relationship inherent in financial time series data.

2.1. Quantile Regression

In statistical modelling, the concept of quantile regression is used in estimating uncertainty of the predicted variable expressed using conditional quantiles. Quantile regression estimates a quantile of a random variable of the $\tau$ -th $(0<\tau <1)$ quantile of random variable $X$ is defined as follows:

Definition 2.1 Suppose $X$ is a random variable on a probability space $\left(\Omega ,ℱ,P\right)$ is characterised by a (right-continuous) distribution function denoted by $F(\cdot )$ , For any $0<\tau <1$ , the $\tau$ -th quantile of $X$ is given by

$F^{-1}\left(\tau \right)=\text{inf}\left\{x:F\left(x\right)\ge \tau \right\}.$ (1)

Linear Quantile Regression Model

[8] developed the theory for the estimation of the quantiles of a variable $Y_t$ which is assumed to be a linear function of other variables. The linear quantile regression model can be defined as

$Q_y\left(\tau |x_t\right)={x^{\prime }}_t\beta \left(\tau \right),\tau \in \left(0,1\right),$ (2)

such that $y_t$ satisfies the quantile constraint $P\left(y_t\le {x^{\prime }}_t\beta \left(\tau \right)|X_t\right)=\tau$ , where $\beta \left(\tau \right)$ are regression coefficients that depend on $\tau$ . Quantile regression estimation determines an estimate for the ${\tau }^{\text{th}}$ conditional quantile, ${\widehat{Q}}_y\left(\tau ,X_t\right)=\widehat{\beta }\left(\tau \right)$ , by minimizing the following function

$\widehat{\beta }\left(\tau \right)=\text{arg}\underset{\beta }{\text{min}}\frac{1}{n}{\displaystyle \sum }_{t=1}^n{\rho }_{\tau }\left(y_t-{x^{\prime }}_t\beta \left(\tau \right)\right),$ (3)

where ${\rho }_{\tau }(\cdot )$ is the quantile loss function defined as

${\rho }_{\tau }\left(u_t\right)=\{\begin{array}{ll}\tau u_t\left(\tau \right),\hfill & \text{if}{u}_t\left(\tau \right)\ge 0\hfill \\ \left(1-\tau \right)u_t\left(\tau \right),\hfill & \text{if}{u}_t\left(\tau \right)<0\hfill \end{array}$ (4)

and $y_t\left(\tau \right)=y_t-{x^{\prime }}_t\left(\tau \right)$ . The quantile estimator in Equation (2) provides more information on the whole conditional distribution of $y_t$ as a function of the $X_t$ , and more robust estimates under the presence of outliers and non-linearities, when compared to the ordinary least squares estimator.

The linear association assumption for the conditional quantile of the dependent variable given the covariates, $Q_y\left(\tau |x_t\right)={x^{\prime }}_t\beta \left(t\right)$ , can be generally restrictive when non-linearities and outliers exist [22]. There is extensive literature on non-parametric quantile regression that relaxes the strict linearity assumption in the covariates. As non-parametric quantile regression models become more complex, local fitting and back-fitting to accommodate new components become more computationally demanding, and the literature has evolved towards non-linear quantile regression models.

In this paper, the focus is on the case of the following non-linear is on the case of the following non-linear association, $Q_y\left(\tau |x_t\right)=h_{\tau }\left(x_t\right)$ where $h_{\tau }(\cdot )$ is some unknown non-linear function. The estimation strategy proposed is to approximate $h_{\tau }\left(x_t\right)$ with neural networks using the universal approximation property. Specifically, we assume that there exists a neural network with a function $G_{\tau }\left(x_t,w\right)$ , that can approximate $h_{\tau }\left(x_t\right)$ accurately. Before we illustrate how the methodology is implemented, the next section presents a discussion on how neural networks can approximate $h_{\tau }\left(x_t\right)$ .

2.2. Deep Neural Networks

In this paper, the feed-forward neural network is utilized to approximate $h_{\tau }\left(x_t\right)$ . This architecture consists of an input layer of covariates, the hidden layer(s) where non-linear transformation of covariates occur, and the output layer that gives the final prediction. Each hidden layer has several interconnected neurons relating it to both the previous and the next ones. Specifically, information flows from one layer to the other, via neurons only in one direction, and the connections correspond to weights. Optimizing a loss function with respect to these weights makes the neural network capable of learning.

Given that $L$ denote the total number of hidden layers, a measure for the depth of a neural network, and $J^{\left(l\right)}$ denotes the total number of neurons at layer $l$ , a measure for its width. First, a general definition of a deep (multi-layer) feed-forward neural network is presented. Let ${\sigma }_l(\cdot ),l=0,\cdots ,L$ be the activation function used at the $l$ -th layer, that is applied element-wise and induces non-linearity. The ReLU activation function, ${\sigma }_l(\cdot )=\text{max}\left(\cdot ,0\right)$ , for $l=1,\cdots ,L-1$ is utilized and a linear one for the output layer, $l=L$ . The ReLU activation function, despite its utter simplicity, seems to be the most popular choice in practical applications [23]. The output of the $l$ -th layer denoted by $g^{\left(l\right)}$ , is a vector of length equal to the number of $J^{\left(L\right)}$ neurons in that layer, such that $g^{\left(0\right)}=x_t$ . The overall structure of a deep (multi-layer) feed-forward neural network is given by:

$G_{\tau }\left(x_t,w\right)=g^{\left(L\right)}\left(g^{\left(L-1\right)}\left(\cdots \left(g^{\left(l\right)}(\cdot )\right)\right)\right),$ (5)

where

$g^{\left(l\right)}\left(x_t\right)={\sigma }_l\left(W^{\left(l-1\right)}{g}^{\left(l-1\right)}+b^{\left(l\right)}\right),l=1,\cdots ,L,$ (6)

$W^{\left(l\right)}$ is a $J^{\left(l\right)}\times J^{\left(l-1\right)}$ matrix of weights, $b^{\left(l\right)}$ is a $J^{\left(l\right)}\times 1$ vector of biases giving an overall vector $w=\left(\text{vec}{\left(W^{\left(0\right)}\right)}^{\prime },\cdots ,\text{vec}{\left(W^{\left(L\right)}\right)}^{\prime },b^{\left(1\right)}{}^{\prime },\cdots ,b^{\left(L\right)}{}^{\prime }\right)$ of trainable parameters of dimensions $j^{\left(l\right)}\left(1+J^{\left(l-1\right)}\right)$ total number of parameters in each hidden layer $l$ , $J^{\left(0\right)}=p$ and $j^{\left(L\right)}=1$ .

The increase in complexity can occur, either by letting $L\to \infty$ , which stands for deep learning, or by letting $J^{\left(l\right)}\to \infty$ . While both ways asymptotically deliver the same results [18], the approximation error has been shown to decline exponentially with $L$ , [24] but only polynomially with $J^{\left(l\right)}$ , providing some evidence for the prevalent use of deep learning. The work of [17] also proposes an alternative approximation theory for sparse deep learning. According to various approximation theorems (see the theoretical results in [25] and [26]), $G_{\tau }\left(x_t,w\right)$ can approximate arbitrarily well $h_{\tau }\left(x_t\right)$ , such that, for any $\epsilon >0$ ,

$\underset{t}{\text{sup}}\left|G_{\tau }\left(x_t,w\right)-h_{\tau }\left(x_t\right)\right|<\epsilon$ (7)

The $\epsilon$ -approximation can be seen as a sieve type non-parametric estimation bond, where $\epsilon$ can become arbitrarily small by increasing the complexity of $G_{\tau }\left(x_t,w\right)$ .

2.3. Deep Neural Network Quantile Regression Model

Assume that the conditional quantile follows a non-linear relationship $Q_y\left(\tau |x_t\right)=h_{\tau }\left(x_t\right)$ and there exists a function $G_{\tau }\left(x_t,w\right)$ , that can $\epsilon$ -approximate $h_{\tau }\left(x_t\right)$ , in Equation (7). Specifically, it is assumed that there exists a neural network with a function $G_{\tau }\left(x_t,w\right)$ , that can approximate $h_{\tau }\left(x_t\right)$ . Under this assumption, the conditional quantile function can be defined as the following approximation

$Q_y\left(\tau |x_t\right)=G_{\tau }\left(x_t,w\right)+O\left(\epsilon \right),$ (8)

where $G_{\tau }\left(x_t,w\right)$ is the unknown non-linear function to be estimated in order to approximate $h_{\tau }\left(x_t\right)$ . The deep neural network conditional quantile estimate is obtained from the solution of the following minimization problem:

$Q_y\left(\tau |x_t\right)=\arg \underset{w}{\min }\frac{1}{n}{\displaystyle \sum }_{t=1}^n{\rho }_{\tau }\left(y_t-G_{\tau }\left(x_t,w\right)\right),$ (9)

where $w=\left(\text{vec}{\left(W^{\left(0\right)}\right)}^{\prime },\cdots ,\text{vec}{\left(W^{\left(L\right)}\right)}^{\prime }\right)$ , $b^{\left(1\right)}{}^{\prime },\cdots ,b^{\left(L\right)}{}^{\prime }$ contains all model parameters, and $G_{\tau }\left(x_t,w\right)$ will govern whether the model is parametric or non-parametric. If the number of neurons and layers is small, then the model is parametric, if the above number becomes large, then the model becomes non-parametric, since the number of estimated parameters increases with the sample size, similar to sieve non-parametric approximations.

The network parameters $w$ are estimated by minimizing the quantile regression loss function (also called the pinball loss) for a given quantile level $\tau \in \left(0,1\right)$ .

$S_n\left(w\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\rho }_{\tau }\left(Y_i-g_{\tau }\left(x_t,w\right)\right)$ (10)

where ${\rho }_{\tau }\left(u\right)=u\left(\tau -\mathbb{I}\left\{u<0\right\}\right)$ , and $\mathbb{I}\{\cdot \}$ is the indicator function.

At the population level, the non-parametric quantile estimation is to find a measurable function $f^{\text{*}}:R^d\to R$ satisfying

$f^{\text{*}}:=\text{arg}\underset{f}{\text{min}}{R}^r\left(f\right)=\text{arg}\underset{f}{\text{min}}{E}_Z\left\{{\rho }_{\tau }\left(\left(x_t,w\right)\right)\right\}$ (11)

where $E_Z$ means that the expectation is taken with respect to the distribution of $Z$ . If the conditional $\tau$ -th quantile of $\nu$ given $X$ is 0 and $E\left(\left|\nu \right||X=x\right)<\infty$ for all $x\in X$ , then the true regression function $f_0$ is the optimal solution of $f^{\text{*}}$ on $X$ .

For a given random function $g_{\tau }\left(x_t,w\right):R^d\to R$ , where $Z\equiv \left(X,Y\right)$ is a random vector independent of $g_{\tau }\left(x_t,w\right)$ , the risk of $g_{\tau }\left(x_t,w\right)$ under the loss function ${\rho }_{\tau }(\cdot )$ is defined by

$R^{\tau }\left(g_{\tau }\left(x_t,w\right)\right)=E_Z\left\{{\rho }_{\tau }\left(Y_i-g_{\tau }\left(x_t,w\right)\right)\right\}$ (12)

where $\theta \in \Theta \subset {ℝ}^d$ .

Let ${\mathcal{D}}_n=\left\{Z_1,\cdots ,Z_n\right\}$ denote the observed data, where $Z_i=\left(X_i,Y_i\right)\in {\left[0,1\right]}^p\times ℝ$ . In applications, when a random sample $S\equiv {\left\{\left(X_i,Y_i\right)\right\}}_{i=1}^n$ is available, the Deep Quantile Regression estimator ${\widehat{q}}_{\tau }\left(x\right)$ is defined as the minimizer of the empirical risk:

$R_n^{\tau }\left(f\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\rho }_{\tau }\left(x_t,w\right)$ (13)

The goal is to construct an estimator of $f_0$ within a certain class of functions ${ℱ}_n$ by minimising the empirical risk, that is

${\widehat{f}}_n\in \text{arg}\underset{f\in {ℱ}_n}{\text{min}}{R}_n^{\tau }\left(f\right),$ (14)

where ${\widehat{f}}_n$ is the deep quantile regression (DQR) estimator, which serves as the empirical risk minimizer (ERM). The function class ${ℱ}_n$ consists of deep neural networks (DNNs). The architecture of a DQR generally comprises multiple hidden layers with non-linear activation functions, similar to a standard feed-forward neural network.

2.3.1. The Deep Quantile Regression Model

Let ${\left\{\left(X_t,Y_t\right)\right\}}_{t=1}^T$ be a strict stationary and mixing stochastic process defined on a probability space $\left(\Omega ,ℱ,\mathbb{P}\right)$ , where $Y_t\in ℝ$ denotes the target financial variable such as asset returns, $X_t\in X\subset {ℝ}^d$ denotes the lagging conditioning information like past returns, realized variance and macroeconomic factors. The target conditional $\tau$ -th quantile of $Y_t$ given $X_t=x$ is denoted as $g_0\left(x\right)=Q_{\tau }\left(Y_t|X_t=x\right)$ . Let ${ℱ}_{DNN}\left(L,W,S,B\right)$ define a class of deep neural networks parametrized by $L$ -denote number of hidden layers (depth), $W$ -denote maximum width per later, $S$ -denote sparsity constraint (total number of non-zero parameters) and $B$ -denote bounded parameter size $\left({\Vert \theta \Vert }_{\infty }\le B\right)$ . This conditional quantile is parametrized via a deep neural network (DNN) architecture $f\left(x;\theta \right)$ , where $\theta \in \Theta \subset {ℝ}^{P_T}$ represents the neural network parameter vector (weights and biases).

The DQR estimator ${\widehat{g}}_T(\cdot )=f\left(\cdot ;{\widehat{\theta }}_T\right)={\widehat{\theta }}_T$ is obtained by solving the empirical risk minimization (ERM):

${\widehat{\theta }}_T=\arg \underset{\widehat{\theta }\in \Theta }{\min }\frac{1}{T}{\displaystyle \sum }_{t=1}^T{\rho }_{\tau }\left(Y_t-f\left(X_t;\theta \right)\right)$ (15)

where ${\rho }_{\tau }\left(u\right)=u\left(\tau -I(u<0)\right)$ is the classical pinball loss function. The population pseudo-true parameter vector is defined as:

${\theta }^{*}=\text{arg}\underset{\theta \in \text{Θ}}{\text{min}}E\left[{\rho }_{\tau }\left(Y_t-f\left(X_t;\theta \right)\right)\right]$ (16)

2.3.2. Assumptions

The following regularity assumptions are assumed because financial time series data possess heavy tails and conditional heteroscedasticity:

Assumption 1. The true conditional quantile function $g_0$ belongs to a Holder continuous function space ${\mathscr{H}}^{\beta }\left(\mathcal{X}\right)$ with smoothness index $\beta >0$ . That is, all partial derivatives of $g_0$ up to the order $\lfloor \beta \rfloor$ exist and are Lipschitz continuous of order $\beta -\lfloor \beta \rfloor$ .

Assumption 2. The stochastic process $\left\{\left(X_t,Y_t\right)\right\}$ is strictly stationary and absolutely regular ($\beta$ -mixing) or strongly mixing ($\alpha$ -mixing). The mixing coefficients $\beta \left(m\right)$ decaying exponentially: $\beta \left(m\right)\le c_1\exp \left(-c_{2}m\right)$ for some constants $c_1,c_2>0$ and $\alpha$ -mixing coefficients satisfying $\alpha \left(m\right)=O\left(m^{-\gamma }\right)$ for $\gamma >\frac{2+\delta }{\delta }$ where $\delta >0$ , ensuring the memory of the shocks decays at a geometric or high polynomial rate.

The financial time series must satisfy a form of weak dependence (such as $\alpha$ -mixing or $\beta$ -mixing) to ensure the consistency and asymptotic normality of the deep quantile regression estimators. This implies that data points far apart in time become approximately independent. This bounds the long-term memory property of asset volatility.

Assumption 3. The conditional distribution $F_{Y|X}\left(y|x\right)$ is continuous and has a continuous, differentiable density $f_{Y|X}\left(y|x\right)$ that is strictly bounded away from zero and infinity in a neighbourhood around the true conditional quantile

$0<\underset{\_}{f}\le f_{Y|X}\left(f\left(x;{\theta }^{*}\right)|x\right)\le \overline{f}<\infty .$

Assumption 4. The parameter space $\Theta$ is a compact subset of ${ℝ}^P$ . The activation functions (like ReLU) ensure that $\theta \to f\left(x;\theta \right)$ is twice continuously differentiable with uniformly bounded gradients:

$E\left[{\Vert {\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right)\Vert }^{2+\delta }\right]<\infty .$

Assumption 5. The positive definite matrices $J\left({\theta }^{*}\right)$ (Hessian counterpart) and $\Omega \left({\theta }^{*}\right)$ (information matrix) are defined as follows:

$J\left({\theta }^{*}\right)=E\left[f_{Y|X}\left(f\left(X_t;{\theta }^{*}\right)|X_t\right)\cdot {\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right){\nabla }_{\theta }f{\left(X_t;{\theta }^{*}\right)}^{\text{T}}\right]\ge 0$ (17)

$\Omega \left({\theta }^{*}\right)=\underset{T\to \infty }{\text{lim}}\text{Var}\left(\frac{1}{\sqrt{T}}{\displaystyle \sum }_{t=1}^T{\varphi }_{\tau }\left(Y_t-f\left(X_t;{\theta }^{\text{*}}\right)\right){\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right)\right)>0$ (18)

where ${\psi }_{\tau }\left(u\right)=\tau -I\left(u<0\right)$ is the sub-gradient score function.

Assumption 6. The neural network utilizes a continuous, piece-wise linear activation function (such as ReLU). The parameters scale with the sample size $T$ such that

$L\asymp \log T,W\asymp T^{\frac{d}{2\beta +d}}\log T,S\asymp T^{\frac{d}{2\beta +d}}\log T$

The assumption guarantees that the temporal memory of shocks decays exponentially over time.

2.3.3. Asymptotic Properties of the DQR Estimator

This section outlines the asymptotic properties of the Deep Quantile Regression (DQR) estimator, which merges the flexibility of deep neural networks with the robustness of quantile regression. The assumption of independent and identically distributed (i.i.d.) observations does not hold for financial market data due to the presence of volatility clustering, regime shifts, and strong temporal dependence. To address these dependencies, we replace i.i.d. assumptions with regularized dynamic constraints. Consequently, the framework for deriving the asymptotic properties of the DQR estimator relies on strict stationarity and either $\beta$ -mixing (absolutely regular) or $\alpha$ -mixing (strongly mixing) conditions to effectively capture the data dynamics.

Consistency

The consistency of the DQR estimator under temporal dependence is given by the following theorem.

Theorem 1 (Consistency). Under assumptions A1-A3, if the VC-dimension or covering number of the deep neural network function class $ℱ=\left\{f\left(\cdot ;\theta \right):\theta \in \Theta \right\}$ grows at rate slower than sample size $\left(o\left(T\right)\right)$ , then:

$\widehat{\theta }\stackrel{p}{\to }{\theta }^{*}$ as $T\to \infty$ (19)

Proof. The convexity Lemma ([27] [28]): Point wise convergence in probability of convex functions defined on an open set to a convex limit function implies uniform convergence on compact subsets. Therefore, for any compact set $\Theta \subset ℝ$ containing $q_{\tau }\left(x\right)$ :

$\underset{\theta \in \Theta }{\sup }\left|M_{n,B}\left(\theta \right)-M\left(\theta \right)\right|\stackrel{p}{\to }0$ (20)

Standard M-estimation arguments Argmin continuous mapping Theorem [29] imply consistency:

${\widehat{q}}_{\tau }\left(x\right)=\arg \underset{\theta }{\min }{M}_{n,B}\left(\theta \right)\stackrel{p}{\to }{q}_{\tau }\left(x\right)=\arg \underset{\theta }{\min }M\left(\theta \right)$ (21)

The following M-estimation lemma is utilized in deriving the consistency of the deep quantile regression estimator under temporal dependence.

Lemma 1. Let ${\mathbb{M}}_T\left(\theta \right)=\frac{1}{T}{\displaystyle \sum }_{t=1}^T{\rho }_{\tau }\left(Y_t-f\left(X_t;\theta \right)\right)$ and $\mathbb{M}\left(\theta \right)=\mathbb{E}\left[{\rho }_{\tau }\left(Y_t-f\left(X_t;\theta \right)\right)\right]$ .

If ${\sup }_{\theta \in \Theta }\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|\stackrel{p}{\to }0$ as $T\to \infty$ (under the Uniform Law of Large Numbers), ${\theta }^{\text{*}}$ is the unique minimizer of $\mathbb{M}\left(\theta \right)$ , and for any $\epsilon >0$ , ${\inf }_{\Vert \theta -{\theta }^{*}\Vert \ge \epsilon }\mathbb{M}\left(\theta \right)>\mathbb{M}\left({\theta }^{*}\right)$ . Then any sequence of empirical minimizers ${\widehat{\theta }}_T$ satisfying the inequality ${\mathbb{M}}_T\left({\widehat{\theta }}_T\right)\le {\mathbb{M}}_T\left({\theta }^{\text{*}}\right)+o_p\left(1\right)$ converges in probability to ${\theta }^{\text{*}}$ , i.e., ${\widehat{\theta }}_T\stackrel{p}{\to }{\theta }^{\text{*}}$.

By the definition of the estimator ${\widehat{\theta }}_T$ as the empirical minimizer, we have the initial condition:

${\mathbb{M}}_T\left({\widehat{\theta }}_T\right)\le {\mathbb{M}}_T\left({\theta }^{\text{*}}\right)$ (22)

This statement can be manipulated by adding and subtracting the population limits $\mathbb{M}\left({\widehat{\theta }}_T\right)$ and $\mathbb{M}\left({\theta }^{\text{*}}\right)$ :

$\begin{array}{c}\mathbb{M}\left({\widehat{\theta }}_T\right)-\mathbb{M}\left({\theta }^{\text{*}}\right)=\left(\mathbb{M}\left({\widehat{\theta }}_T\right)-{\mathbb{M}}_T\left({\widehat{\theta }}_T\right)\right)+\left({\mathbb{M}}_T\left({\widehat{\theta }}_T\right)-{\mathbb{M}}_T\left({\theta }^{\text{*}}\right)\right)\\ +\left({\mathbb{M}}_T\left({\theta }^{\text{*}}\right)-\mathbb{M}\left({\theta }^{\text{*}}\right)\right)\end{array}$ (23)

Using the empirical minimizer inequality (${\mathbb{M}}_T\left({\widehat{\theta }}_T\right)-{\mathbb{M}}_T\left({\theta }^{\text{*}}\right)\le 0$), we bound the equation from Equation (22), to obtain

$\mathbb{M}\left({\widehat{\theta }}_T\right)-\mathbb{M}\left({\theta }^{\text{*}}\right)\le \left(\mathbb{M}\left({\widehat{\theta }}_T\right)-{\mathbb{M}}_T\left({\widehat{\theta }}_T\right)\right)+\left({\mathbb{M}}_T\left({\theta }^{\text{*}}\right)-\mathbb{M}\left({\theta }^{\text{*}}\right)\right)$ (24)

To make this statement uniform across the parameter space Θ, we take the supremum over all possible parameters $\theta \in \Theta$ :

$0\le \mathbb{M}\left({\widehat{\theta }}_T\right)-\mathbb{M}\left({\theta }^{*}\right)\le 2\underset{\theta \in \Theta }{\sup }\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|$ (25)

Therefore, to prove that $\mathbb{M}\left({\widehat{\theta }}_T\right)\stackrel{p}{\to }\mathbb{M}\left({\theta }^{\text{*}}\right)$, it is necessary and sufficient to prove that the uniform empirical process term converges to zero:

$\underset{\theta \in \Theta }{\sup }\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|\stackrel{p}{\to }0$ (26)

Given that the sample data ${\left\{\left(X_t,Y_t\right)\right\}}_{t=1}^T$ is strictly stationary and $\alpha$ -mixing (Assumption A1), the classical Glivenko-Cantelli theorem designed for i.i.d. random variables cannot be directly applied. A a sequence of independent blocks is constructed by dividing the time horizons $1,\cdots ,T$ into $2k_T$ blocks of length $m_T$ , such that $T=2k_T{m}_T$ . Let the odd blocks be denoted by $H_j=\left\{\left(j-1\right)2m_T+1,\cdots ,\left(j-1\right)2m_T+m_T\right\}$ for $j=1,\cdots ,k_T$ . Let the even blocks be denoted by $I_j=\left\{\left(2j-1\right)m_T+1,\cdots ,2j{m}_T\right\}$ for $j=1,\cdots ,k_T$ . The empirical process can be rewritten as the sum of its odd and even structural components given by

${\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)=\frac{1}{2}\left({\mathbb{M}}_{T,\text{odd}}\left(\theta \right)-\mathbb{M}\left(\theta \right)\right)+\frac{1}{2}\left({\mathbb{M}}_{T,\text{even}}\left(\theta \right)-\mathbb{M}\left(\theta \right)\right)$ (27)

By Berbee’s Lemma, there exists a new sequence of coupled variables ${\left\{\left({\tilde{X}}_t,{\tilde{Y}}_t\right)\right\}}_{t=1}^T$ such that: Within each block $H_j$ , the joint distribution of $\left({\tilde{X}}_t,{\tilde{Y}}_t\right)$ is identical to the original time series $\left(X_t,Y_t\right)$ . The blocks $H_1,H_2,\cdots ,H_{k_T}$ are mutually independent. The total variation distance between the true dependent joint distribution and the decoupled independent block sequence satisfies:

$\mathbb{P}\left({\displaystyle \underset{t=1}{\overset{T}{\cup }}\left\{\left(X_t,Y_t\right)\ne \left({\tilde{X}}_t,{\tilde{Y}}_t\right)\right\}}\right)\le 2k_T\cdot \alpha \left(m_T\right)$ (28)

Under Assumption A1, the mixing coefficients decay geometrically ($\alpha \left(m\right)\le c_1{\text{e}}^{-c_{2}m}$ ). By setting the tracking block size parameters to scale logarithmically, $m_T=\frac{2}{c_2}\log T$ , the coupling failure bound simplifies:

$2k_T\alpha \left(m_T\right)\le \frac{T}{m_T}{c}_1{\text{e}}^{-2\log T}=\frac{c_1}{T\cdot m_T}\to 0\text{as}T\to \infty$ (29)

Thus, with probability approaching 1 ($1-o\left(1\right)$ ), we can execute the rest of the uniform convergence proof as if the blocks were perfectly independent.

Let $\mathcal{G}=\left\{\left(x,y\right)\mapsto {\rho }_{\tau }\left(y-f\left(x;\theta \right)\right):\theta \in \Theta \right\}$ . Given that the pinball loss ${\rho }_{\tau }(\cdot )$ is Lipschitz continuous with parameter $\text{max}\left(\tau ,1-\tau \right)\le 1$ , the covering number of the transformed loss class is bounded by the covering number of the network function class $ℱ$ :

$N\left(\epsilon ,\mathcal{G},L_1\left({\mathbb{P}}_T\right)\right)\le N\left(\epsilon ,ℱ,L_1\left({\mathbb{P}}_T\right)\right)$ (30)

Using the standard symmetrization theorem via Rademacher variables ${\sigma }_j\in \left\{-1,1\right\}$ over the $k_T$ independent blocks, we bound the expected supremum of the empirical process:

$\begin{array}{l}\mathbb{E}\left[\underset{\theta \in \Theta }{\text{sup}}\left|{\mathbb{M}}_{T,\text{odd}}\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|\right]\\ \le 2{\mathbb{E}}_{\tilde{X},\tilde{Y}}\left[{\mathbb{E}}_{\sigma }\left[\underset{f\in ℱ}{\text{sup}}\left|\frac{1}{k_T}{\displaystyle \sum }_{j=1}^{k_T}{\sigma }_j\cdot \left(\frac{1}{m_T}{\displaystyle \sum }_{t\in H_j}{\rho }_{\tau }\left({\tilde{Y}}_t-f\left({\tilde{X}}_t\right)\right)\right)\right|\right]\right]\end{array}$ (31)

Applying Dudley’s Entropy Integral Theorem to control the sub-Gaussian process via the log covering numbers yields:

${\mathbb{E}}_{\sigma }\left[\underset{f\in ℱ}{\sup }\cdots \right]\le \frac{C}{\sqrt{k_T}}{\displaystyle {\int }_0^{2B}\sqrt{\log N\left(\epsilon ,ℱ,L_1\left({\mathbb{P}}_T\right)\right)}\text{d}\epsilon }$ (32)

By the explicit theorem condition, the VC-dimension or covering capacity satisfies $\log N\left(\epsilon ,ℱ,L_1\left({\mathbb{P}}_T\right)\right)=o\left(T/m_T\right)=o\left(k_T\right)$ . As a result, the integral scales bounded under the square root, and the scaling coefficient ensures:

$\frac{1}{\sqrt{k_T}}\cdot \sqrt{o\left(k_T\right)}=o\left(1\right)\to 0\text{as}T\to \infty$ (33)

This establishes the Uniform Law of Large Numbers for our time series process:

$\underset{\theta \in \Theta }{\sup }\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|\stackrel{p}{\to }0$ (34)

Finally, the convergence of the objective functions is translated into the convergence of the parameter vector ${\widehat{\theta }}_T$ . Suppose ${\widehat{\theta }}_T$ does not converge to ${\theta }^{\text{*}}$ . Then there exists a subsequence and a fixed radius ${\epsilon }_0>0$ such that $\Vert {\widehat{\theta }}_T-{\theta }^{\text{*}}\Vert \ge {\epsilon }_0$ infinitely often. Define the compact subset ${\Theta }_{{\epsilon }_0}=\left\{\theta \in \Theta :\Vert \theta -{\theta }^{*}\Vert \ge {\epsilon }_0\right\}$ . By Assumption A1 (${\theta }^{\text{*}}$ is a unique, interior minimum) and Assumption A2 ($f_{Y|X}$ is continuous and bounded away from zero, implying strict convexity of the population risk), we have a positive gap value $\delta >0$ :

$\underset{\theta \in {\Theta }_{{\epsilon }_0}}{\inf }\mathbb{M}\left(\theta \right)=\mathbb{M}\left({\theta }^{*}\right)+\delta$ (35)

Now, trace the following sequence of events:

$\begin{array}{l}\mathbb{M}\left({\widehat{\theta }}_T\right)\ge \mathbb{M}\left({\theta }^{*}\right)+\delta \Rightarrow \mathbb{M}\left({\widehat{\theta }}_T\right)-{\mathbb{M}}_T\left({\widehat{\theta }}_T\right)+{\mathbb{M}}_T\left({\widehat{\theta }}_T\right)-\mathbb{M}\left({\theta }^{*}\right)\ge \delta \\ \Rightarrow \underset{\theta \in \Theta }{\sup }\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|+{\mathbb{M}}_T\left({\theta }^{*}\right)-\mathbb{M}\left({\theta }^{*}\right)\ge \delta \\ \Rightarrow 2\underset{\theta \in \Theta }{\sup }\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|\ge \delta \end{array}$ (36)

Taking the probability limits of both sides:

$\mathbb{P}\left(\Vert {\widehat{\theta }}_T-{\theta }^{\text{*}}\Vert \ge {\epsilon }_0\right)\le \mathbb{P}\left(\underset{\theta \in \text{Θ}}{\text{sup}}\left|{\mathbb{M}}_T\left(\theta \right)-\mathbb{M}\left(\theta \right)\right|\ge \frac{\delta }{2}\right)$ (37)

Since we proved that the right-hand side of Equation (37) converges to zero as $T\to \infty$ , the probability of the parameter falling outside the epsilon ball must also vanish:

$\underset{T\to \infty }{\text{lim}}\mathbb{P}\left(\Vert {\widehat{\theta }}_T-{\theta }^{\text{*}}\Vert \ge {\epsilon }_0\right)=0$ (38)

Hence,

${\widehat{\theta }}_T\stackrel{p}{\to }{\theta }^{\text{*}}$

Given the $\alpha$ -mixing property in assumption (A1) and compactness in assumption (A3), the classical i.i.d. Uniform Law of Large numbers (ULLN) for mixing processes is replaced by an empirical process central limit theorem for dependent sequences. The empirical risk uniformly converges to the population risk given by:

${\displaystyle \sum }_{\theta \in \text{Θ}}\left|\frac{1}{T}{\displaystyle \sum }_{t=1}^T{\rho }_{\tau }\left(Y_t-f\left(X_t;\theta \right)\right)-E\left[{\rho }_{\tau }\left(Y_t-f\left(X_t;\theta \right)\right)\right]\right|\stackrel{p}{\to }0$ (39)

Since ${\theta }^{*}$ uniquely minimizes the continuous, strictly convex expected pinball loss profile over Θ, any sequence minimizing the empirical objective must merge in probability with ${\theta }^{*}$ .

Asymptotic Normality

The following theorem describes the asymptotic normality property of the DQR estimator.

Theorem 2 (Asymptotic Normality). Under assumptions A1-A4, as $T\to \infty$ , the normalized DQR estimator satisfies:

$\sqrt{T}\left({\widehat{\theta }}_T-{\theta }^{*}\right)\stackrel{d}{\to }\mathcal{N}\left(0,J{\left({\theta }^{*}\right)}^{-1}\Omega \left({\theta }^{*}\right)J{\left({\theta }^{*}\right)}^{-1}\right)$ (40)

Proof. The pinball loss function ${\rho }_{\tau }\left(u\right)$ is non-differentiable at zero and as a result the standard Taylor series expansion fails. The Knight’s Identify is utilized to isolate the parameter deviation vector $u=\sqrt{T}\left(\theta -{\theta }^{*}\right)$ :

${\rho }_{\tau }\left(x-\nu \right)-{\rho }_{\tau }\left(x\right)=-\nu {\psi }_{\tau }\left(x\right)+{\displaystyle {\int }_0^{\nu }\left(I\left(x\le s\right)-I\left(x\le 0\right)\right)\text{d}s}$ (41)

Let ${\nu }_t\left(\theta \right)=f\left(X_t;\theta \right)-f\left(X_t;{\theta }^{*}\right)$ , applying the identity to our empirical objective gives a localized objective function $Z_T\left(u\right)$ :

$\begin{array}{c}{Z}_T\left(u\right)={\displaystyle \sum }_{t=1}^T\left[{\rho }_{\tau }\left(Y_t-f\left(X_t;{\theta }^{*}+T^{-1/2}u\right)\right)-{\rho }_{\tau }\left(Y_t-f\left(X_t;{\theta }^{*}\right)\right)\right]\\ =-\frac{1}{\sqrt{T}}{\displaystyle \sum }_{t=1}^T{u}^{\text{T}}{\nabla }_{\theta }f\left(X_t;{\theta }^{\text{*}}\right){\psi }_{\tau }\left({\epsilon }_t\right)\\ +{\displaystyle \sum }_{t=1}^T{\displaystyle {\int }_0^{T^{-1/2}{u}^T{\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right)}\left(I\left({\epsilon }_t\le s\right)-I\left({\epsilon }_t\le 0\right)\right)\text{d}s}+o_p\left(1\right)\end{array}$ (42)

where ${\epsilon }_t=Y_t-f\left(X_t;{\theta }^{*}\right)$ represents the true conditional quantile residual.

Taking the conditional expectation of the integral term in Equation (42) with respect to the filtration ${ℱ}_{t-1}$ . and applying a localized Taylor series expansion to the conditional distribution function $F_{Y|X}$ :

$\begin{array}{l}E\left[{\displaystyle {\int }_0^{T^{-1/2}{u}^{\text{T}}{\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right)}\left(I\left({\epsilon }_t\le s\right)-I\left({\epsilon }_t\le 0\right)\right)\text{d}s}|X_t\right]\\ \approx \frac{1}{2T}{u}^{\text{T}}{\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right)f_{Y|X}\left(f\left(X_t;{\theta }^{*}\right)|X_t\right){\nabla }_{\theta }f{\left(X_t;{\theta }^{*}\right)}^{\text{T}}u\end{array}$ (43)

Summing the right-hand side of Equation (43) over $t=1$ to $T$ , by the Ergodic Theorem for stationarity processes, this converges in probability to:

$\frac{1}{2}{u}^{\text{T}}J\left({\theta }^{*}\right)u$ (44)

Define the score vector $W_t={\psi }_{\tau }\left({\epsilon }_t\right){\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right)$ . By definition of a conditional quantile, $E\left[{\psi }_{\tau }\left({\epsilon }_t\right)|X_t\right]=0$ . meaning $W_t$ forms a martingale difference sequence (MDS) if the conditioning set matches the filtration perfectly. However, under $\alpha$ -mixing, $W_t$ is a weakly dependent stationary process. Applying the Gordin-Root Functional Central Limit Theorem for $\alpha$ -mixing sequences:

$\frac{1}{\sqrt{T}}{\displaystyle \sum }_{t=1}^T{W}_t\stackrel{d}{\to }\mathcal{N}\left(0,\Omega \left({\theta }^{*}\right)\right)$ (45)

where $\Omega \left({\theta }^{*}\right)$ is the Long-Run Variance Covariance matrix (LRV) given by

$\omega \left({\theta }^{*}\right)={\Sigma }_0+{\displaystyle \sum }_{k=1}^{\infty }\left({\Sigma }_k+{\Sigma }_k^{\text{T}}\right)\text{with}{\Sigma }_k=E\left[W_t{W}_{t-k}^{\text{T}}\right]$ (46)

when the neural network successfully captures all temporal dynamics such that residuals are conditionally independent (no remaining autocorrelation in the signs of the residuals), $\Omega \left({\theta }^{*}\right)$ simplifies to

$\Omega \left({\theta }^{*}\right)=\tau \left(1-\tau \right)E\left[{\nabla }_{\theta }f\left(X_t;{\theta }^{*}\right){\nabla }_{\theta }f{\left(X_t;{\theta }^{*}\right)}^{\text{T}}\right]$

However, if quantile correlation remains, the full long-run variance structure is required.

Substituting the components in Equation (42), the objective function behaves asymptotically as:

$Z_T\left(u\right)=-u^{\text{T}}\left(\frac{1}{\sqrt{T}}{\displaystyle \sum }_{t=1}^T{W}_t\right)+\frac{1}{2}{u}^{\text{T}}J\left({\theta }^{*}\right)u+o_p\left(1\right)$ (47)

By the Convexity Lemma, the minimizer $\widehat{u}=\sqrt{T}\left({\widehat{\theta }}_T-{\theta }^{*}\right)$ converges to the minimizer of the limit quadratic form given by:

$\widehat{u}=J{\left({\theta }^{*}\right)}^{-1}\left(\frac{1}{\sqrt{T}}{\displaystyle \sum }_{t=1}^T{W}_t\right)+o_p\left(1\right)$ (48)

Substituting the normal limit of the score vector gives the result:

$\sqrt{T}\left({\widehat{\theta }}_T-{\theta }^{*}\right)\stackrel{d}{\to }\mathcal{N}\left(0,J{\left({\theta }^{*}\right)}^{-1}\Omega \left({\theta }^{*}\right)J{\left({\theta }^{*}\right)}^{-1}\right)$ (49)

To obtain the convergence rate for deep quantile regression (DQR) estimator, Theorem 3.4.1 in [30] (use the oracle inequality framework) is applied. The goal is to find the optimal balance between the network’s ability to approximate a smooth function and the statistical penalty of using a high-capacity model.

Rate of Convergence

The following theorem describes the rate of convergence property of the DQR estimator.

Theorem 3 (Rate of Convergence). Under Assumptions A1-A4, if the network hyperparameters are optimally balanced relative to the dataset length $T$ and dimension $d$ , the excess risk and the $L_2\left({\mathbb{P}}_X\right)$ convergence rate of the deep quantile estimator satisfies

${\Vert {\widehat{g}}_T-g_0\Vert }_{L_2\left(P_X\right)}^2=\mathcal{O}\left(T^{-\frac{2\beta }{2\beta +d}}{\log }^{\xi }T\right)$ (50)

where $\xi >0$ is a structural logarithmic multiplier. This rate matches the information-theoretic minimax-optimal rate for non-parametric estimation, up to a log factor.

Proof. The proof of the rate of convergence for the DQR estimator relies on the convexity of the objective function and the uniform convergence of the empirical risk. The Convexity Lemma (Pollard, 1991) is utilized to establish uniform convergence from point wise convergence, avoiding the need for smoothness of the score function.

The metric error is decomposed into approximation error (neural expressivity) and stochastic error (variance)

$\begin{array}{l}E\left[{\rho }_{\tau }\left(Y-{\widehat{g}}_T\right)\right]-E\left[{\rho }_{\tau }\left(Y-g_0\right)\right]\\ =\underset{\text{Approximation Error}\left(\text{Bias}\right)}{\underbrace{\underset{f\in {ℱ}_{\text{DNN}}}{\text{inf}}\left(E\left[{\rho }_{\tau }\left(Y-f\right)\right]-E\left[{\rho }_{\tau }\left(Y-g_0\right)\right]\right)}}+\underset{\text{Stochastic Error}\left(\text{Variance}\right)}{\underbrace{{\epsilon }_{\text{stoch}}}}\end{array}$ (51)

By performing a second-order Taylor series expansion of the population objective function around $g_0\left(X_t\right)$ and applying Assumption A3, the excess quantile risk is linked to the $L_2$ function norm:

$\begin{array}{l}E\left[{\rho }_{\tau }\left(Y_t-f\left(X_t\right)\right)|X_t\right]-E\left[{\rho }_{\tau }\left(Y_t-g_0\left(X_t\right)\right)|X_t\right]\\ ={\displaystyle {\int }_{g_0\left(X_t\right)}^{f\left(X_t\right)}\left(f\left(X_t\right)-y\right)f_{Y|X}\left(y|X_t\right)\text{d}y}\ge \frac{\underset{\_}{f}}{2}{\left(f\left(X_t\right)-g_0\left(X_t\right)\right)}^2\end{array}$ (52)

Integrating over the distribution of $X$ yields the key geometric margin inequality given by

$E\left[{\rho }_{\tau }\left(Y-f\right)\right]-E\left[{\rho }_{\tau }\left(Y-g_0\right)\right]\ge C_0{\Vert f-g_0\Vert }_{L_2\left(P_X\right)}^2$

By deep network approximation results for Holder spaces ${\mathscr{H}}^{\beta }$ , for a network satisfying Assumption (A4), there exists an oracle network function $f_T^{\text{*}}\in {ℱ}_{\text{DNN}}$ satisfying the following inequality:

${\Vert f_T^{*}-g_0\Vert }_{\infty }\le c\cdot {\left(W\cdot L\right)}^{-\frac{2\beta }{d}}\asymp T^{-\frac{\beta }{2\beta +d}}$

Thus, the approximation error contribution scales bounded as $\mathcal{O}\left(T^{-\frac{\beta }{2\beta +d}}\right)$ .

To bound ${\epsilon }_{stoch}$ given dependent time series observations, the blocking technique is utilized. The sequence $\left\{1,\cdots ,T\right\}$ is divided into $2k_t$ blocks of size $m_T$ , where $T=2k_T{m}_T$ . Let $H_j$ capture the odd blocks and $I_j$ capture the even blocks. By Berbee’s Lemma, we can construct an independent copy sequence $\left\{\left({\tilde{X}}_t,{\tilde{Y}}_t\right)\right\}$ such that the joint distribution of each block matches the original, but blocks are independent of each other. The probability of the coupling failing is bounded by the mixing coefficient sum: $2k_T\beta \left(m_T\right)$ .

Using Assumption (A2), we set $m_T\asymp \text{log}T$ to drive the coupling failure rate to zero exponentially fast:

$2k_T\beta \left(m_T\right)\to 0$ as $T\to \infty$

Applying the Rademecher complexity/VC-dimension bound $\left(\text{VCdim}\left({ℱ}_{DNN}\right)\le c\cdot SL\log \left(S\right)\right)$ over the coupled independent blocks, the stochastic empirical process variation scales with the capacity metric:

${\epsilon }_{\text{stoch}}=\mathcal{O}\left(\frac{\text{VCdim}\left({ℱ}_{\text{DNN}}\right)\cdot m_T}{T}\right)$ (53)

Substituting the neural network size profiles from Assumption (A4) in Equation (53)

${\mathcal{O}}_P\left(\frac{T^{\frac{d}{2\beta +d}}\cdot {\log }^{2}T\cdot \log T}{T}\right)=\mathcal{O}\left(T^{-\frac{2\beta }{2\beta +d}}{\log }^{3}T\right)$ (54)

Combining the symmetric bounds derived from Equation (53) and Equation (54), the overall error is minimized when the bias and variance convergence speeds are matched as given by

${\Vert {\widehat{g}}_T-g_0\Vert }_{L_2{P}_X}^2=\mathcal{O}\left(T^{-\frac{2\beta }{2\beta +d}}\right)+{\mathcal{O}}_P\left(T^{-\frac{2\beta }{2\beta +d}}{\log }^{3}T\right)$ (55)

2.4. Change-Point Detection

In this section, the Deep Quantile Regression (DQR) framework for detecting change points in financial time series is discussed. The focus shifts from estimating a static conditional quantile to identifying structural breaks in the network parameters or in the underlying quantile function. The mathematical formulation of the change-point detection problem can be expressed as follows.

The Change-Point Model

Let ${\left\{\left(X_t,Y_t\right)\right\}}_{t=1}^T$ be a observed piecewise stationary stochastic process defined on a probability space $\left(\Omega ,ℱ,\mathbb{P}\right)$ , observed over a total time horizon $T$ , where $Y_t\in ℝ$ represents asset returns, and $X_t\in \mathcal{X}\subset {ℝ}^d$ represents lagged market features. Assume the existence of a single true structural break at time index $T^{\text{*}}$ corresponding to the true change-point fraction ${\lambda }^{\text{*}}=T^{\text{*}}/T\in \left(0,1\right)$ . The conditional $\tau$ -quantile function transitions between two distinct operational states governed by a deep neural network function class ${ℱ}_{\text{DNN}}$ :

$Q_{\tau }\left(Y_t|X_t=x\right)=\{\begin{array}{ll}f\left(x;{\theta }_1^{\text{*}}\right)\hfill & \text{for}1\le t\le T^{\text{*}}\hfill \\ f\left(x;{\theta }_2^{\text{*}}\right)\hfill & \text{for}{T}^{\text{*}}+1\le t\le T\hfill \end{array}$ (56)

where ${\theta }_1^{\text{*}}\ne {\theta }_2^{\text{*}}$ are interior points of a compact parameter space $\Theta \subset {ℝ}^p$ .

Assuming the time series can be partitioned into $K_0+1$ distinct local segments, separated by $K_0$ true change-points denoted by the time indices $T_1^{\text{*}},T_2^{\text{*}},\cdots ,T_{K_0}^{\text{*}}$ . For consistency in notations, let $T_0^{\text{*}}=0$ and $T_{K_0+1}^{\text{*}}=T$ . The fraction of the timeline at which the structural breaks occur as change-point fractions is defined as:

${\lambda }_k^{\text{*}}=\frac{T_k^{\text{*}}}{T},\text{where}0={\lambda }_0^{\text{*}}<{\lambda }_1^{\text{*}}<\cdots <{\lambda }_{K_0}^{\text{*}}<{\lambda }_{K_0+1}^{\text{*}}=1$ (57)

Within each segment $t\in \left(T_{k-1}^{\text{*}},T_k^{\text{*}}\right]$ , the true conditional $\tau$ -th quantile function of $Y_t$ given $X_t=x$ is parameterized by a segment-specific deep neural network parameter vector ${\theta }_k^{\text{*}}$ :

$Q_{\tau }\left(Y_t|X_t=x\right)=f\left(x;{\theta }_k^{\text{*}}\right)\text{for}{T}_{k-1}^{\text{*}}<t\le T_k^{\text{*}}$ (58)

where $f\left(\cdot ;\theta \right)$ is a deep neural network mapping $\mathcal{X}\to ℝ$ , parameterized by $\theta \in \Theta \subset {ℝ}^p$ . The structural break is active at $\tau$ if ${\theta }_k^{\text{*}}\ne {\theta }_{k+1}^{\text{*}}$

To locate a single change-point ($K_0=1$ ) at a location $t$ , we optimize the split empirical check profile. The DQR-CPD estimator minimizes the segmented empirical pinball loss profile ${\rho }_{\tau }\left(u\right)=u\left(\tau -\mathbb{I}\left(u<0\right)\right)$ :

$\left(\widehat{T},{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)=\text{arg}\underset{\begin{array}{c}t\in \left[\lfloor ϵT\rfloor ,\lfloor \left(1-ϵ\right)T\rfloor \right]\\ {\theta }_1,{\theta }_2\in \Theta \end{array}}{\text{min}}\left\{{\displaystyle \sum }_{i=1}^t{\rho }_{\tau }\left(Y_i-f\left(X_i;{\theta }_1\right)\right)+{\displaystyle \sum }_{i=t+1}^T{\rho }_{\tau }\left(Y_i-f\left(X_i;{\theta }_2\right)\right)\right\}$

The estimated change-point fraction is defined as $\widehat{\lambda }=\widehat{T}/T$ . Let the segmented empirical risk for any candidate split point $t$ be defined using the pinball loss ${\rho }_{\tau }\left(u\right)=u\left(\tau -\mathbb{I}\left(u<0\right)\right)$ :

${\mathbb{M}}_T\left(t;{\theta }_1,{\theta }_2\right)=\frac{1}{T}{\displaystyle \sum }_{i=1}^t{\rho }_{\tau }\left(Y_i-f\left(X_i;{\theta }_1\right)\right)+\frac{1}{T}{\displaystyle \sum }_{i=t+1}^T{\rho }_{\tau }\left(Y_i-f\left(X_i;{\theta }_2\right)\right)$ (59)

The estimated change-point $\widehat{T}$ and its fraction $\widehat{\lambda }=\widehat{T}/T$ are obtained by minimizing this joint profile:

$\left(\widehat{T},{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)=\arg \underset{t\in \left[\lfloor ϵT\rfloor ,\lfloor \left(1-ϵ\right)T\rfloor \right],{\theta }_1,{\theta }_2\in \Theta }{\min }{\mathbb{M}}_T\left(t;{\theta }_1,{\theta }_2\right)$ (60)

where $ϵ>0$ is a trimming parameter preventing the candidate break from getting too close to the sample edges.

2.5. Asymptotic Properties of the DQR-CPD Estimator

This section presents the asymptotic properties for the Deep Quantile Regression Change-Point Detection (DQR-CPD) estimator under financial time-series dependence. The derivation of the asymptotic properties of the DQR-CPD estimator for a change point under serial dependence, relies on a strong-mixing concentration inequality that bounds the non-differentiable empirical risk process theory for dependent data.

2.5.1. Assumptions

Establishing consistency and rates of convergence for change-point locations under deep learning quantile models requires imposing strict decay properties on the memory of the system and standard constraints on density profiles.

Assumption 7. Within each structural segment $\left(T_{k-1}^{*},T_k^{*}\right]$ , the sub-process $\left\{\left(X_t,Y_t\right)\right\}$ is strictly stationary and strongly mixing ($\alpha$ -mixing) with exponentially decaying coefficients: $\alpha \left(m\right)\le c_1\exp \left(-c_{2}m\right)$ for constants $c_1,c_2>0$ . The mixing coefficients satisfy $\alpha \left(m\right)=O\left(m^{-\gamma }\right)$ for $\gamma >\frac{2+\delta }{\delta }$ , ensuring that past market shocks decay quickly enough to avoid confounding localized breaks.

Definition 2.2 (mixing (strong mixing)). The sequence is $\alpha$ -mixing if the mixing coefficients $\alpha \left(m\right)$ decay to zero as $m\to \infty$ :

$\alpha \left(m\right)=\underset{A\in {ℳ}_{-\infty }^t,B\in {ℳ}_{t+m}^{\infty }}{\text{sup}}\left|P\left(A\cap B\right)-P\left(A\right)P\left(B\right)\right|\to 0$

where ${ℳ}_a^b$ is the $\sigma$ -algebra generated by ${\left\{X_t,Y_t\right\}}_{t=a}^b$ .

Assumption 8. The conditional distribution function $F_{Y|X}\left(y|x\right)$ possesses a continuous conditional density $f_{Y|X}\left(y|x\right)$ that is uniformly bounded away from zero and infinity in a neighbourhood around all segment profiles:

$0<\underset{\_}{f}\le f_{Y|X}\left(f\left(x;{\theta }_k^{\text{*}}\right)+v|x\right)\le \overline{f}<\infty$ for $\left|v\right|\le \delta ,\forall k$

Assumption 9. (Quantile Multi-Segment Identifiability): For every segment $k$ , the parameter ${\theta }_k^{*}$ is a unique interior point of the compact space $\Theta$ . Additionally, the magnitudes of the parameter jumps at each true break point are non-vanishing with sample size:

$\underset{1\le k\le K_0}{\text{min}}{\Vert {\theta }_{k+1}^{\text{*}}-{\theta }_k^{\text{*}}\Vert }_2\ge \underset{\_}{v}>0$

Assumption 10. The network function $\theta \mapsto f\left(x;\theta \right)$ is twice continuously differentiable with uniformly bounded gradients: ${\sup }_{x,\theta }{\Vert {\nabla }_{\theta }f\left(x;\theta \right)\Vert }_2\le C_g<\infty$ . Furthermore, the true change-point fractions are well-separated:

$\underset{1\le k\le K_0+1}{\text{min}}\left({\theta }_k^{\text{*}}-{\theta }_{k-1}^{\text{*}}\right)\ge {ϵ}_0>0$

This guarantees that each regime contains enough data points proportional to $T$ to allow the neural networks to learn the conditional distributions.

Assumption 11. (Network Complexity Growth): The VC-dimension or covering capacity of ${ℱ}_{\text{DNN}}$ scales sub-linearly relative to the sample size:

$\log N\left(\delta ,{ℱ}_{\text{DNN}},L_1\left(\mathbb{P}\right)\right)=o\left(T\right).$

To bound the capacity of the network space under dependence, its metric entropy is computed. Let $\mathcal{N}\left(ϵ,{ℱ}_n,{\Vert \cdot \Vert }_{\infty }\right)$ be the $ϵ$ -covering number of ${ℱ}_n$ under the supremum norm.

To derive the asymptotic properties theorems, first we establish three structural lemmas handling uniform convergence over mixing series and loss metric localization.

Lemma 2. (Uniform Concentration via Berbee’s Blocking) Let ${\mathbb{M}}_T\left(\lambda ,{\theta }_1,{\theta }_2\right)$ be the normalized empirical objective function. Under Assumptions A7 and A11, the empirical objective converges uniformly to its population counterpart:

$\underset{\lambda ,{\theta }_1,{\theta }_2}{\text{sup}}\left|{\mathbb{M}}_T\left(\lambda ,{\theta }_1,{\theta }_2\right)-\mathbb{M}\left(\lambda ,{\theta }_1,{\theta }_2\right)\right|\stackrel{p}{\to }0$ as $T\to \infty$ (61)

Proof. Divide the time horizon into $2k_T$ alternating blocks of size $m_T$ . Let $m_T\asymp \text{log}T$ , driving the total number of blocks $k_T=T/\left(2m_T\right)$ . By Berbee’s Lemma, we couple the dependent sequence to an independent block sequence $\left\{\left({\tilde{X}}_t,{\tilde{Y}}_t\right)\right\}$ . The total variation distance between their joint distributions is bounded by $2k_T\alpha \left(m_T\right)$ . Given the exponential decay in Assumption A7:

$2k_T\alpha \left(m_T\right)\le \frac{T}{\log T}{c}_1{\text{e}}^{-c_2\log T}\to 0$

Given that the blocks are mutually independent with high probability, we apply Dudley’s entropy integral over the network capacity defined in A4. Since $\log N\left(\delta ,{ℱ}_{\text{DNN}},L_1\right)=o\left(T/m_T\right)$ , the uniform law of large numbers holds over the parameter space.

Lemma 3. (Knight’s Loss Identity) For any real numbers $x$ and $v$ , the pinball loss function satisfies:

${\rho }_{\tau }\left(x-v\right)-{\rho }_{\tau }\left(x\right)=-v{\psi }_{\tau }\left(x\right)+{\displaystyle {\int }_0^v\left(\mathbb{I}\left(x\le s\right)-\mathbb{I}\left(x\le 0\right)\right)\text{d}s}$ (62)

where ${\psi }_{\tau }\left(x\right)=\tau -\mathbb{I}\left(x<0\right)$ represents the sub-gradient score.

Lemma 4. (DNN Covering Bounds): For a neural network class ${ℱ}_n$ with piecewise linear activations (e.g., ReLU), depth $L_n$ , and at most $W_n$ parameters, the covering number satisfies:

$\log \mathcal{N}\left(ϵ,{ℱ}_n,{\Vert \cdot \Vert }_{\infty }\right)\le c\cdot W_n{L}_n\log \left(\frac{M_n}{ϵ}\right)$

2.5.2. Consistency

To establish the consistency of the Deep Quantile Regression-based Change-Point Detection (DQR-CPD) estimator under temporal dependence, we cannot rely on classical independent and identically distributed (i.i.d.) empirical process theory. Financial data requires handling non-smooth loss functions under mixing dependencies, where the network parameters and the structural break fraction must be resolved simultaneously.

The following are the necessary foundational lemmas and theorems, to prove both global parameter consistency and the super-consistent oracle rate ${\mathcal{O}}_P\left(T^{-1}\right)$ for the change-point fraction.

Theorem 4. Under Assumptions A7-A11, the joint DQR-CPD estimator is globally consistent:

$\left(\widehat{\lambda },{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)\stackrel{p}{\to }\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)$ as $T\to \infty$ (63)

Proof. By definition of the argmin configuration, the empirical minimizers satisfy ${\mathbb{M}}_T\left(\widehat{\lambda },{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)\le {\mathbb{M}}_T\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)$. Standard M-estimation framing yields:

$\begin{array}{c}0\le \mathbb{M}\left(\widehat{\lambda },{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)-\mathbb{M}\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)\\ \le \mathbb{M}\left(\widehat{\lambda },{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)-{\mathbb{M}}_T\left(\widehat{\lambda },{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)+{\mathbb{M}}_T\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)-\mathbb{M}\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)\\ \le 2\underset{\lambda ,{\theta }_1,{\theta }_2}{\text{sup}}\left|{\mathbb{M}}_T\left(\lambda ,{\theta }_1,{\theta }_2\right)-\mathbb{M}\left(\lambda ,{\theta }_1,{\theta }_2\right)\right|\end{array}$ (64)

By Lemma 2, the right-hand side converges to zero in probability. Thus:

$\mathbb{M}\left(\widehat{\lambda },{\widehat{\theta }}_1,{\widehat{\theta }}_2\right)\stackrel{p}{\to }\mathbb{M}\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)$ (65)

Under Assumption A8 and A9, the expected loss function $\mathbb{M}$ is uniquely and strictly minimized at the true parameter structure $\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)$ . The well-separation condition dictates that for any small neighbourhood radius $\eta >0$ :

$\underset{\Vert \left(\lambda ,{\theta }_1,{\theta }_2\right)-\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)\Vert \ge \eta }{\text{inf}}\mathbb{M}\left(\lambda ,{\theta }_1,{\theta }_2\right)>\mathbb{M}\left({\lambda }^{\text{*}},{\theta }_1^{\text{*}},{\theta }_2^{\text{*}}\right)$ (66)

Therefore, the parameter vector configuration must converge in probability to the true target values:

$\widehat{\lambda }\stackrel{p}{\to }{\lambda }^{\text{*}},{\widehat{\theta }}_1\stackrel{p}{\to }{\theta }_1^{\text{*}},{\widehat{\theta }}_2\stackrel{p}{\to }{\theta }_2^{\text{*}}.$ (67)

2.5.3. Rate of Convergence

Theorem 5. Fast Rate of Convergence for the Break Fraction (Oracle Rate): Under assumptions A7-A11, the estimated change-point fraction converges to the true break fraction at the information-theoretic super-consistent oracle rate:

$\left|\widehat{\lambda }-{\lambda }^{\text{*}}\right|={\mathcal{O}}_P\left(T^{-1}\right)$ (68)

Proof. To derive the oracle convergence rate $\left|\widehat{\lambda }-{\lambda }^{\text{*}}\right|={\mathcal{O}}_P\left(T^{-1}\right)$ for the change-point fraction in a deep quantile regression framework, the arguments are formalized using localized empirical process theory and random walk dynamics with negative drift.

The estimated change-point $\widehat{T}$ is the global minimizer of the profile empirical risk:

${\mathbb{V}}_T\left(t\right)=\underset{{\theta }_1,{\theta }_2\in \text{Θ}}{\text{min}}\left\{\frac{1}{T}{\displaystyle \sum }_{t=1}^t{\rho }_{\tau }\left(Y_t-f\left(X_t;{\theta }_1\right)\right)+\frac{1}{T}{\displaystyle \sum }_{t=t+1}^T{\rho }_{\tau }\left(Y_t-f\left(X_t;{\theta }_2\right)\right)\right\}$ (69)

Let $\widehat{\lambda }=\widehat{T}/T$ , and let ${\lambda }^{\text{*}}=T^{\text{*}}/T$ be the true break fraction. To establish that $\left|\widehat{\lambda }-{\lambda }^{\text{*}}\right|={\mathcal{O}}_P\left(T^{-1}\right)$, it is mathematically equivalent to proving that the absolute time-step error is stochastically bounded:

$\left|\widehat{T}-T^{\text{*}}\right|={\mathcal{O}}_P\left(1\right)$ (70)

Given that Theorem 4 establishes global consistency, The search window can be restricted to a shrinking neighbourhood $I_T\left(\eta \right)=\left\{t:\left|t-T^{\text{*}}\right|\le \eta T\right\}$ for an arbitrarily small $\eta >0$ . To show that $\left|\widehat{\lambda }-{\lambda }^{\text{*}}\right|={\mathcal{O}}_P\left(T^{-1}\right)$, it is mathematically equivalent to proving that the discrete time error is bounded in probability: $\left|\widehat{T}-T^{\text{*}}\right|={\mathcal{O}}_P\left(1\right)$.

Without loss of generality, assume the candidate split point falls to the left of the true transition point: $t\in \left[T^{\text{*}}-\eta T,T^{\text{*}}\right)$ . We isolate the profile empirical risk difference ${\mathbb{V}}_T\left(t\right)-{\mathbb{V}}_T\left(T^{\text{*}}\right)$ across the misclassified observation domain $i\in \left\{t+1,\cdots ,T^{\text{*}}\right\}$ :

${\mathbb{V}}_T\left(t\right)-{\mathbb{V}}_T\left(T^{\text{*}}\right)=\frac{1}{T}{\displaystyle \sum }_{i=t+1}^{T^{\text{*}}}\left[{\rho }_{\tau }\left(Y_i-f\left(X_i;{\widehat{\theta }}_2\left(t\right)\right)\right)-{\rho }_{\tau }\left(Y_i-f\left(X_i;{\theta }_1^{\text{*}}\right)\right)\right]+{\mathcal{O}}_P\left(T^{-1}\right)$ (71)

Applying Lemma 3 (Knight’s Identity) to the summation terms where ${\epsilon }_i=Y_i-f\left(X_i;{\theta }_1^{\text{*}}\right)$ and ${\Delta }_i=f\left(X_i;{\widehat{\theta }}_2\left(t\right)\right)-f\left(X_i;{\theta }_1^{*}\right)$ , the excess objective function is decomposed into its deterministic drift and stochastic noise components:

${\mathbb{V}}_T\left(t\right)-{\mathbb{V}}_T\left(T^{\text{*}}\right)={\mathbb{D}}_T\left(t\right)+{\mathbb{W}}_T\left(t\right)+o_P\left(T^{-1}\right)$ (72)

The drift term collects the conditional expectations of the integrated indicator terms:

${\mathbb{D}}_T\left(t\right)=\frac{1}{T}{\displaystyle \sum }_{i=t+1}^{T^{\text{*}}}\mathbb{E}\left[{\displaystyle {\int }_0^{{\Delta }_i}\left(\mathbb{I}\left({\epsilon }_i\le s\right)-\mathbb{I}\left({\epsilon }_i\le 0\right)\right)\text{d}s}|X_i\right]$ (73)

Applying a Taylor expansion to the conditional distribution function $F_{Y|X}$ combined with the density threshold $\underset{\_}{f}$ from Assumption A8:

${\mathbb{D}}_T\left(t\right)\ge \frac{\underset{\_}{f}}{2T}{\displaystyle \sum }_{i=t+1}^{T^{\text{*}}}{\left(f\left(X_i;{\widehat{\theta }}_2\left(t\right)\right)-f\left(X_i;{\theta }_1^{\text{*}}\right)\right)}^2$ (74)

Since $t\in I_T\left(\eta \right)$ , ${\widehat{\theta }}_2\left(t\right)\stackrel{p}{\to }{\theta }_2^{\text{*}}$. Assumption A9 guarantees a non-vanishing parameter jump ($\Vert {\theta }_1^{\text{*}}-{\theta }_2^{\text{*}}\Vert \ge \underset{\_}{v}$ ). Therefore, the expected squared functional difference behaves like a strictly positive constant $C_1>0$ . Summing across the $m=T^{\text{*}}-t$ misclassified observations shows that the drift grows strictly linearly with the step count:

${\mathbb{D}}_T\left(t\right)\ge C_0\frac{T^{\text{*}}-t}{T}\text{for}{C}_0=\frac{\underset{\_}{f}{C}_1}{2}>0$ (75)

The noise process is defined by the centred sub-gradient scores:

${\mathbb{W}}_T\left(t\right)=-\frac{1}{T}{\displaystyle \sum }_{i=t+1}^{T^{\text{*}}}{\psi }_{\tau }\left({\epsilon }_i\right)\left(f\left(X_i;{\widehat{\theta }}_2\left(t\right)\right)-f\left(X_i;{\theta }_1^{\text{*}}\right)\right)$ (76)

Given that $\mathbb{E}\left[{\psi }_{\tau }\left({\epsilon }_i\right)|X_i\right]=0$ , this forms a zero-mean weakly dependent process under $\alpha$ -mixing (Assumption A7). By applying a Hájek-Rényi type maximal inequality for mixing sequences, its variance scales linearly with the window width $m=T^{\text{*}}-t$ , meaning its supremum across the sequence behaves sub-linearly

$\underset{t:T^{*}-t=m}{\max }\left|{\mathbb{W}}_T\left(t\right)\right|={\mathcal{O}}_P\left(\frac{\sqrt{m\log \log m}}{T}\right)$ (77)

For the estimated change-point $\widehat{T}$ to be selected over the true break point $T^{\text{*}}$ , the empirical objective condition requires

${\mathbb{V}}_T\left(\widehat{T}\right)-{\mathbb{V}}_T\left(T^{\text{*}}\right)\le 0,$

which requires the stochastic noise to completely overwhelm the deterministic drift (${\mathbb{D}}_T\left(\widehat{T}\right)\le \left|{\mathbb{W}}_T\left(\widehat{T}\right)\right|$ ). A sequence of Dyadic shells is constructed to evaluate this probability outside a bounded distance $M>0$ :

$\mathbb{P}\left(\left|\widehat{T}-T^{*}\right|>M\right)\le {\displaystyle \sum }_{k=\lfloor {\log }_{2}M\rfloor }^{\infty }\mathbb{P}\left(\underset{2^k\le T^{*}-t<2^{k+1}}{\max }\frac{\left|{\mathbb{W}}_T\left(t\right)\right|}{{\mathbb{D}}_T\left(t\right)}\ge 1\right)$ (78)

Substituting the linear drift boundary and applying Doob’s maximal inequality to the mixing sequence sums across each shell yields:

${\displaystyle \sum }_{k=\lfloor {\log }_{2}M\rfloor }^{\infty }\frac{\mathbb{E}\left[{\left({{\displaystyle \sum }}_{i=T^{*}-2^{k+1}}^{T^{*}}{\psi }_{\tau }\left({\epsilon }_i\right){\Delta }_i\right)}^2\right]}{{\left(C_0{2}^k\right)}^2}\le {\displaystyle \sum }_{k=\lfloor {\log }_{2}M\rfloor }^{\infty }\frac{C_2{2}^{k+1}}{C_0^2{2}^{2k}}=\frac{4C_2}{C_0^2}\cdot \frac{1}{M}$ (79)

By choosing the tracking boundary variable $M$ to be sufficiently large, this probability is driven below any arbitrary error threshold $ϵ>0$ :

$\frac{4C_2}{C_0^2\cdot M}<ϵ\Rightarrow \left|\widehat{T}-T^{\text{*}}\right|={\mathcal{O}}_P\left(1\right)$ (80)

Dividing the discrete temporal parameter boundary by the global time horizon $T$ :

$\left|\widehat{\lambda }-{\lambda }^{\text{*}}\right|=\left|\frac{\widehat{T}}{T}-\frac{T^{\text{*}}}{T}\right|=\frac{1}{T}\left|\widehat{T}-T^{\text{*}}\right|=\frac{1}{T}\cdot {\mathcal{O}}_P\left(1\right)={\mathcal{O}}_P\left(T^{-1}\right)$ (81)

2.5.4. Asymptotic Normality

To derive the asymptotic normality of the deep neural network parameters inside each stable segment, first we linearize the non-differentiable pinball loss function. This is achieved via a generalized Bahadur representation extended to non-linear sieve spaces (like deep neural networks) under serial dependence.

Let ${\theta }_n$ represent the vectorized weights and biases of the deep neural network $f\left(X_t,\theta \right)\in {ℱ}_n$ . Let ${\theta }_n^0$ be the sieve parameter vector that optimizes the population conditional quantile within a stable segment. The parameter deviation is scaled by defining a local parameter $u$ :

$\theta ={\theta }_n^0+\frac{u}{\sqrt{n}}$

The localized empirical objective function using the pinball loss ${\rho }_q\left(z\right)=z\left(q-I\left\{z<0\right\}\right)$ is written as:

${\mathbb{G}}_n\left(u\right)={\displaystyle \sum }_{t=1}^n\left[{\rho }_q\left(Y_t-f\left(X_t,{\theta }_n^0+\frac{u}{\sqrt{n}}\right)\right)-{\rho }_q\left(Y_t-f\left(X_t,{\theta }_n^0\right)\right)\right]$ (82)

The goal is to approximate ${\mathbb{G}}_n\left(u\right)$ by a quadratic form in $u$ .

While the global pinball loss is non-differentiable, the neural network activation function $f\left(X_t,\theta \right)$ itself is assumed to be almost everywhere differentiable with respect to its weights (smoothed ReLU or standard GeLU/Sigmoid activations). Let $J_t\left(\theta \right)={\nabla }_{\theta }f\left(X_t,\theta \right)$ be the $d_n\times 1$ gradient vector of the network with respect to its parameters, where $d_n=\text{dim}\left({\theta }_n\right)$ grows with $n$ . By Taylor expansion:

$f\left(X_t,{\theta }_n^0+\frac{u}{\sqrt{n}}\right)-f\left(X_t,{\theta }_n^0\right)=\frac{1}{\sqrt{n}}{J}_t^{\text{T}}u+R_{nt}\left(u\right)$ (83)

where the remainder satisfies ${\text{max}}_t\left|R_{nt}\left(u\right)\right|=o_p\left(n^{-1/2}\right)$ due to weight bounds $M_n$ .

The Knight’s Identity is utilized in decomposing the objective function for the pinball loss function:

${\rho }_q\left(x-y\right)-{\rho }_q\left(x\right)=-y{\psi }_q\left(x\right)+{\displaystyle {\int }_0^y\left(I\left\{x\le s\right\}-I\left\{x\le 0\right\}\right)\text{d}s}$ (84)

where ${\psi }_q\left(x\right)=q-I\left\{x<0\right\}$ yields the structural quantile residuals. Substituting $x=Y_t-f\left(X_t,{\theta }_n^0\right)$ and $y=\frac{1}{\sqrt{n}}{J}_t^{\text{T}}u+R_{nt}\left(u\right)$ into the identity transforms ${\mathbb{G}}_n\left(u\right)$ into two main components:

${\mathbb{G}}_n\left(u\right)=-\frac{1}{\sqrt{n}}{\displaystyle \sum }_{t=1}^n{\psi }_q\left(Y_t-f\left(X_t,{\theta }_n^0\right)\right)J_t^{\text{T}}u+{\displaystyle \sum }_{t=1}^n{\mathbb{I}}_{nt}\left(u\right)+o_p\left(1\right)$ (85)

where the stochastic integral component is defined as:

${\mathbb{I}}_{nt}\left(u\right)={\displaystyle {\int }_0^{\frac{1}{\sqrt{n}}{J}_t^{\text{T}}u}\left(I\left\{Y_t-f\left(X_t,{\theta }_n^0\right)\le s\right\}-I\left\{Y_t-f\left(X_t,{\theta }_n^0\right)\le 0\right\}\right)\text{d}s}$ (86)

To linearize the integral term ${\displaystyle {\sum }_{t=1}^n{\mathbb{I}}_{nt}\left(u\right)}$ , we separate it into its conditional expectation and an empirical discrepancy process. Taking the expectation conditional on the historical $\sigma$ -algebra ${ℳ}_{-\infty }^{t-1}$ :

$E\left[{\mathbb{I}}_{nt}\left(u\right)|X_t\right]={\displaystyle {\int }_0^{\frac{1}{\sqrt{n}}{J}_t^{\text{T}}u}\left(F_{Y|X}\left(f\left(X_t,{\theta }_n^0\right)+s|X_t\right)-F_{Y|X}\left(f\left(X_t,{\theta }_n^0\right)|X_t\right)\right)\text{d}s}$ (87)

Applying a first-order Taylor expansion to the conditional distribution function $F_{Y|X}$ around the true quantile, where $f_{Y|X}$ is the conditional continuous density function:

$\begin{array}{c}E\left[{\mathbb{I}}_{nt}\left(u\right)|X_t\right]\approx {\displaystyle {\int }_0^{\frac{1}{\sqrt{n}}{J}_t^{\text{T}}u}{f}_{Y|X}\left(f\left(X_t,{\theta }_n^0\right)|X_t\right)\cdot s\cdot \text{d}s}\\ =\frac{1}{2n}{u}^{\text{T}}\left[f_{Y|X}\left(f\left(X_t,{\theta }_n^0\right)|X_t\right)J_t{J}_t^{\text{T}}\right]u\end{array}$ (88)

Summing over $t$ , the deterministic portion of the representation becomes:

${\displaystyle \sum }_{t=1}^{n}E\left[{\mathbb{I}}_{nt}\left(u\right)|X_t\right]=\frac{1}{2}{u}^{\text{T}}\left(\frac{1}{n}{\displaystyle \sum }_{t=1}^n{f}_{Y|X}\left(f\left(X_t,{\theta }_n^0\right)|X_t\right)J_t{J}_t^{\text{T}}\right)u=\frac{1}{2}{u}^{\text{T}}{D}_{n}u$ (89)

where $D_n$ is the generalized design matrix weighted by the conditional density of the time series error process. To complete the Bahadur linearization, we prove that the centered remainder process converges uniformly to zero:

$\underset{\Vert u\Vert \le C}{\text{sup}}\left|{\displaystyle \sum }_{t=1}^n\left({\mathbb{I}}_{nt}\left(u\right)-E\left[{\mathbb{I}}_{nt}\left(u\right)|X_t\right]\right)\right|=o_p\left(1\right)$ (90)

This is where the geometric $\alpha$ -mixing condition is critically deployed. We construct a sequence of growing parameter subsets and split the dependent sum into blocks via a coupling argument (Bosq’s inequality). The variance of the indicator variations over a local neighbourhood scales at a rate of $O\left(n^{-1/2}\Vert u\Vert \right)$ .Under mixing conditions, the covariance terms decay exponentially: ${\displaystyle {\sum }_{m=1}^{\infty }\alpha {\left(m\right)}^{\delta }}<\infty$ . Applying a localized version of the Bernstein inequality for mixing processes, the chaining bound yields a convergence rate for the supremum of $O_p\left(n^{-1/4}{d}_n^{3/4}\right)$ , which vanishes as long as the neural network dimension growth satisfies $d_n^3/n\to 0$ .

The results yield the uniform Bahadur representation for the deep quantile network parameters:

${\mathbb{G}}_n\left(u\right)=-u^{\text{T}}{W}_n+\frac{1}{2}{u}^{\text{T}}{D}_{n}u+o_p\left(1\right)$ (91)

where $W_n=\frac{1}{\sqrt{n}}{\displaystyle {\sum }_{t=1}^n{\psi }_q\left(Y_t-f\left(X_t,{\theta }_n^0\right)\right)J_t}$ . By minimizing this quadratic limit with respect to $u$ (setting the derivative to zero), we obtain the explicit linear mapping of the estimation error:

$\sqrt{n}\left({\theta }_n-{\theta }_n^0\right)=D_n^{-1}{W}_n+o_p\left(1\right)$ (92)

${\theta }_n-{\theta }_n^0=\frac{1}{n}{D}_n^{-1}{\displaystyle \sum }_{t=1}^n{\psi }_q\left(Y_t-f\left(X_t,{\theta }_n^0\right)\right){\nabla }_{\theta }f\left(X_t,{\theta }_n^0\right)+o_p\left(n^{-1/2}\right)$ (93)

By the Central Limit Theorem (CLT) for $\alpha$ -mixing triangular arrays, $W_n\stackrel{d}{\to }\mathcal{N}\left(0,{\Omega }_n\right)$ , where ${\Omega }_n$ is the long-run covariance matrix that accounts for serial correlation in the quantile residuals:

${\Omega }_n=\underset{n\to \infty }{\lim }\frac{1}{n}{\displaystyle \sum }_{t=1}^n{\displaystyle \sum }_{s=1}^{n}E\left[{\psi }_q\left(Y_t-f\left(X_t\right)\right){\psi }_q\left(Y_s-f\left(X_s\right)\right)J_t{J}_s^{\text{T}}\right]$ (94)

Thus, the localized network parameters follow the limiting normal distribution:

$\sqrt{n}\left({\theta }_n-{\theta }_n^0\right)\stackrel{d}{\to }\mathcal{N}\left(0,D_n^{-1}{\Omega }_n{D}_n^{-1}\right)$ (95)

This direct linearization allows the construction of valid Wald-type tests or sandwich covariance estimators to verify if a structural break has truly occurred in the asset’s quantile dynamics.

3. Simulation Results

This section presents the results of a simulation study of the deep quantile regression estimator for change-point estimation. Subsection 4.1 presents the simulation study on consistency of the deep quantile estimator for single change point detection. Subsection 4.2 presents a simulation study on consistency of the deep quantile estimator for multiple change points detection. Subsection 4.3 presents the convergence rate of the deep quantile regression change point estimator.

3.1. Simulation Settings

In this section, we outline the simulation settings used to assess the asymptotic properties of the Deep Quantile Regression Change Point Detection (DQR-CPD) estimator. The asymptotic properties of this estimator were investigated using simulated data generated from commonly used models in financial time series analysis: the ARMA(1,1)-GARCH(1,1) model and the ARMA(1,1)-GJR-GARCH(1,1) model, both with errors following a Student’s t-distribution.

For the simulation study on change-point detection in volatile ARMA-GARCH(1,1) financial data, the DQR network was designed as a fully connected feed-forward architecture. To detect change points in the simulated time series, the generated data was input into a Bidirectional Long Short-Term Memory (BiLSTM) Deep Neural Network Quantile Regression Architecture. This network effectively models non-linear structural shifts and asymmetric financial shocks by utilizing the ReLU (Rectified Linear Unit) activation function across the hidden layers. The output layer employs a linear activation function to compute the targeted pinball loss. Network optimization was primarily driven by the Adam optimizer, set with an initial learning rate between 0.001 and 0.01, which may optionally decay using a learning rate scheduler to stabilize convergence.

The performance of the deep quantile estimator was evaluated by varying the sample size and then determining the location and magnitude of change points. The simulations were conducted for the following quantile values: $\tau \in \left(1\text{\%},2.5\text{\%},5\text{\%},50\text{\%},75\text{\%},90\text{\%}\right)$ . Sample sizes varied between $n=1000,2000,3000$ and 5000, with a total of 1000 Monte Carlo replications. The evaluation metrics included Mean Absolute Error, Median Absolute Deviation, Mean Squared Error, and Root Mean Squared Error for each targeted quantile threshold ($\tau$ ) and for different sample sizes.

3.2. Consistency of the Deep Quantile Estimator for Single Change Point Detection

The following two data generating models are considered to assess the finite sample properties of the deep quantile estimator:

3.2.1. ARMA(1,1)-GARCH(1,1) Model

The simulated data $\left\{X_t\right\}$ is generated using the following ARMA(1,1)-GARCH(1,1) model

$\{\begin{array}{l}{X}_t=0.01+0.5X_{t-1}+{\epsilon }_t-0.2{\epsilon }_{t-1}\\ {\sigma }_t^2=0.05+0.4{\epsilon }_{t-1}+0.25{\sigma }_{t-1}^2\text{where}{\epsilon }_t={\sigma }_t{z}_t,z_t~N\left(0,1\right)\end{array}$ (96)

Figure 1 displays the simulated ARMA(1,1)-GARCH(1,1) model for a sample size of 1000 observations, along with the corresponding plot of the conditional volatility trajectory (illustrating volatility clustering). Parameter for the fitted models were chosen to satisfy the constraints of the GARCH models. The simulated data were split, a change point was introduced, and the two time series were merged to create a combined time series that includes a change point. This was subsequently extended to include multiple change points, with several change points introduced in the datasets. The conditional volatility plot visually illustrates the patterns of volatility clustering, where significant shocks in returns are followed by prolonged periods of elevated baseline variance. This phenomenon is a well-known characteristic of financial time series modelling.

Figure 1. Simulated ARMA(1,1)-GARCH(1,1) process for a sample size of $n=1000$ .

In the simulated data, a single change point is introduced, dividing the model into two distinct regimes. At $t=500$ , the model transitions smoothly from a stable baseline to a high-volatility regime. The model used to implement this structural change is given as follows:

Pre-change point region

$\{\begin{array}{l}{X}_t=0.01+0.4X_{t-1}+{\epsilon }_t-0.1{\epsilon }_{t-1}\\ {\sigma }_t^2=0.02+0.10{\epsilon }_{t-1}^2\end{array}$ (97)

Post-change point regime

$\{\begin{array}{l}{X}_t=0.08+0.6X_{t-1}+{\epsilon }_t-0.3{\epsilon }_{t-1}\\ {\sigma }_t^2=0.15+0.45{\epsilon }_{t-1}^2\end{array}$ (98)

Figure 2 illustrates the time plot of the ARMA(1,1)-GARCH(1,1) model with a structural break at $t=500$ . The regime shift is highlighted in the conditional volatility plot following the introduction a single change point introduced at $t=500$ . The change point in the simulated time series was detected using the BiLSTM Deep Neural Network with a Quantile Regression Architecture. The time series data generated was fed into this BiLSTM DNN for analysis. To evaluate the performance of the Deep Quantile Change Point Estimator, a Monte Carlo simulation was conducted with 1000 replications, focusing on the ARMA(1,1)-GARCH(1,1) dataset, which contains a true structural break at $t=500$ . Table 1 presents the empirical results, including the estimated location of the change point and various performance metrics.

Figure 2. Simulated ARMA(1,1)-GARCH(1,1) process for a sample size of $n=1000$ with a change point at $t=500$ .

Table 1. Simulation results on the consistency of the deep quantile regression change point estimator using the ARMA(1,1)-GARCH(1,1) model.

$n=1000$
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	500	502.214	2.214	4.858	4.000	37.276	6.105
0.25	500	501.793	1.793	3.625	3.000	20.817	4.563
0.50	500	501.209	1.209	2.897	2.000	13.431	3.665
0.75	500	501.417	1.417	3.659	2.500	20.689	4.549
0.95	500	501.810	1.810	4.870	4.000	36.652	6.054
$n=2000$
0.05	1000	1001.084	1.084	2.592	2.000	11.241	3.353
0.25	1000	1000.672	0.672	1.905	1.000	6.274	2.505
0.50	1000	1000.381	0.381	1.465	1.000	3.803	1.950
0.75	1000	1000.695	0.695	1.918	1.000	6.311	2.512
0.95	1000	1001.053	1.053	2.571	2.000	11.082	3.329
$n=5000$
0.05	2500	2500.412	0.412	1.348	1.000	3.514	1.875
0.25	2500	2500.226	0.226	0.941	1.000	1.802	1.342
0.50	2500	2500.114	0.114	0.684	0.000	0.916	0.957
0.75	2500	2500.231	0.231	0.953	1.000	1.825	1.351
0.95	2500	2500.398	0.398	1.332	1.000	3.468	1.862

In each simulation, the estimates of the change point depend on the location of the change point, the sample size, and the magnitude of the change. The estimates are most accurate when the true change point, is situated near the median of the sample. This positioning ensures that the pre-change and post-change samples are of approximately equal size. As the sample size increases, the estimation error significantly decreases. Additionally, the bias in the estimates diminishes with larger sample sizes and greater shifts in change, further confirming the consistency of the deep quantile regression estimator for change point detection.

The estimated change point is a random variable that follows the sampling distribution of the estimator. Figure 3 shows the histogram of the estimated change points for the simulated ARMA(1,1)-GARCH(1,1) process with a sample size of $n=1000$ , where the change point occurs at $t=500$ .

Figure 3. Histogram of the estimated change point for the simulated ARMA(1,1)-GARCH(1,1) process for a sample size of $n=1000$ with a change point at $t=500$ .

3.2.2. ARMA(1,1)-GJR-GARCH(1,1) Model with Student t Error Distribution

The simulated data $\left\{X_t\right\}$ is generated using the following ARMA(1,1)-GJR-GARCH(1,1) model with Student t error distribution given as follows:

$\{\begin{array}{l}{X}_t=0.01+0.5X_{t-1}+{\epsilon }_t-0.2{\epsilon }_{t-1}\\ {\sigma }_t^2=0.05+0.4{\epsilon }_{t-1}+0.25{\sigma }_{t-1}^2\text{where}{\epsilon }_t={\sigma }_t{z}_t,z_t~N\left(0,1\right)\end{array}$ (99)

Figure 4 displays the simulated ARMA(1,1)-GJR-GARCH(1,1) model for a sample size of 1000 observations, along with the corresponding plot of the conditional volatility trajectory (illustrating volatility clustering). In the simulated data, a single change point is introduced, dividing the model into two distinct regimes. At $t=500$ , the model transitions smoothly from a stable baseline to a high-volatility regime. The mathematical model used to implement this structural change is given as follows:

Pre-change point region ($t<500$ ):

$\{\begin{array}{l}{X}_t=0.01+0.4X_{t-1}+{\epsilon }_t-0.1{\epsilon }_{t-1}\\ {\sigma }_t^2=0.02+0.10{\epsilon }_{t-1}^2\end{array}$ (100)

Post-change point regime ($t\ge 500$ )

$\{\begin{array}{l}{X}_t=0.08+0.6X_{t-1}+{\epsilon }_t-0.3{\epsilon }_{t-1}\\ {\sigma }_t^2=0.15+0.45{\epsilon }_{t-1}^2\end{array}$ (101)

Figure 4. Simulated ARMA(1,1)-GJR-GARCH(1,1) process for a sample size of $n=1000$ .

Figure 5 shows a time series plot of the ARMA(1,1)-GJR-GARCH(1,1) model, which includes a structural break at $t=500$ . This regime shift is highlighted in the conditional volatility plot, indicating a significant change at that specific point. By incorporating asymmetry (leverage effects) and heavy Student-$t$ tails, the estimation behaviour is altered. In this model, large negative shocks prompt an immediate and pronounced reaction in volatility, as indicated by the sum of the parameters ($\alpha +\gamma =0.45$ ). Additionally, lower quantiles, such as $\tau =0.05$ , are more responsive and can detect the structural shift faster and more accurately than higher quantiles, such as $\tau =0.95$ . The change point was identified using a BiLSTM Deep Neural Network with a Quantile Regression Architecture. Time series data generated from this analysis was input into the BiLSTM DNN for further evaluation. To assess the effectiveness of the Deep Quantile Change Point Estimator, we conducted a Monte Carlo simulation with 1000 replications, focusing on the ARMA(1,1)-GJR-GARCH(1,1) dataset, which has a known structural break at $t=500$ .

Figure 5. Simulated ARMA(1,1)-GJR-GARCH(1,1) process for a sample size of $n=1000$ with a change point at $t=500$ .

Table 2 presents the empirical results of the Deep Quantile Regression Change Point Estimator’s performance across these 1000 replications of the ARMA(1,1)-GJR-GARCH(1,1) model, which includes Student-$t$ innovations and the confirmed structural break at $t=500$ . The table summarizes the estimated location of the change point along with various performance metrics. These metrics include Mean Absolute Error, Median Absolute Deviation, Mean Squared Error, and Root Mean Squared Error for each targeted quantile threshold and for different sample sizes.

Table 2. Simulation results on the consistency of the deep quantile regression change point estimator using the ARMA(1,1)-GJR-GARCH(1,1) model with Student-t error distribution.

$n=1000$
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	500	501.255	1.255	4.065	3.000	26.497	5.148
0.25	500	501.223	1.223	3.417	3.000	18.699	4.324
0.50	500	501.582	1.582	3.370	3.000	18.028	4.246
0.75	500	501.515	1.515	3.927	3.000	23.733	4.872
0.95	500	502.287	2.287	5.511	4.000	46.903	6.849
$n=2000$
0.05	1000	1000.592	0.592	2.115	1.000	8.142	2.853
0.25	1000	1000.614	0.614	1.954	1.000	6.819	2.611
0.50	1000	1000.418	0.418	1.512	1.000	4.103	2.026
0.75	1000	1000.743	0.743	2.052	1.500	7.426	2.725
0.95	1000	1001.216	1.216	3.148	2.000	17.514	4.185
$n=5000$
0.05	2500	2500.183	0.183	0.814	1.000	1.942	1.394
0.25	2500	2500.191	0.191	0.732	0.000	1.411	1.188
0.50	2500	2500.089	0.089	0.541	0.000	0.724	0.851
0.75	2500	2500.215	0.215	0.784	1.000	1.513	1.230
0.95	2500	2500.312	0.312	1.154	1.000	3.102	1.761

Figure 6 displays the histogram of the estimated change points for the simulated ARMA(1,1)-GJR-GARCH(1,1) process with sample sizes of $n=1000$ and $n=5000$ . In this simulation, the change point occurs at $t=500$ for the $n=1000$ and at $t=2500$ for the $n=5000$ . As the sample size increases, the estimated change point moves closer to the actual change point.

The Deep Quantile Regression estimator for single change point estimation consistently demonstrates that the estimator is statistically consistent. As the sample size $n\to \infty$ , the estimated change point ${\tau }_{cp}$ converges in probability to the true change point ${\tau }_0$ , and the regression parameters converge to their true values.

Figure 6. Histogram of change point estimator estimator for a sample size of $n=1000$ with a change point at $t=500$ .

3.3. Consistency of the Deep Quantile Estimator for Multiple Change Point Detection

The consistency of the deep quantile estimator for multiple point detection was investigated similarly by simulating data from two conditional heteroscedastic models commonly used in the analysis of financial time series data: the ARMA(1,1)-GARCH(1,1) model and the ARMA(1,1)-GJR-GARCH(1,1) model, with errors following a Student-t distribution. The performance of the deep quantile estimator was evaluated based on sample size, change point location and the magnitude of change. All the simulations are based on the following values: $\tau \in \left(1\text{\%},2.5\text{\%},5\text{\%},50\text{\%},75\text{\%},90\text{\%}\right)$ , sample size varies as follows: 1000, 2000, 3000, 5000 and the number of Monte Carlo replications is 1000.

Figure 7 illustrates the simulated ARMA (1,1)-GARCH (1,1) model based on a sample of 2000 observations, along with the corresponding trajectory of conditional volatility, which demonstrates the phenomenon of volatility clustering. The time series is divided into four distinct regimes at the true structural breaks occurring at t = 500, t = 1000, and t = 1500. To effectively manage multiple change points using Deep Quantile Regression, the network employs an iterative multi-headed sequence tracking architecture. The model’s estimation performance is evaluated independently at each break point across 1000 replications. Figure 8 presents the histogram of the estimated change points for the simulated ARMA (1,1)-GARCH (1,1) process with a sample size of n = 2000. In this simulation, the change points are identified at t = 500, t = 1000, and t = 1500. As the sample size increases, the estimated change points converge closer to the actual change points.

Table 3 presents performance metrics for three different change point scenarios: volatility burst, negative mean shift, and mixed high aggression. These metrics are provided for various quantiles and include average estimated change points, bias, and error measures. The performance data includes specific metrics such as the average estimated change point ($\widehat{\tau }$ ) and empirical bias, detailing true change points at 1250, 2500, and 3750. The results in the rows highlight the best quantile performance for each event.

Figure 7. Simulated ARMA (1,1)-GARCH (1,1) process for a sample size of 2000, with change points at t = 500, t = 1000 and t = 1500.

Figure 8. Histograms of simulated ARMA (1,1)-GARCH (1,1) process for a sample size of n = 2000 with a change point at t = 500, n = 1000, and t = 1500.

Table 3. Simulation results on the consistency of the deep quantile regression multiple change point estimator using the ARMA (1,1)-GARCH (1,1) model for a sample size of n = 2000 and n = 5000.

Sample size ($n=2000$ )
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	500	501.185	1.185	2.942	2.000	14.112	3.757
	1000	1000.912	0.912	2.511	2.000	10.145	3.185
	1500	1500.495	0.495	1.614	1.000	4.819	2.195
0.25	500	500.841	0.841	2.226	1.5000	8.134	2.852
	1000	100.584	0.584	1.849	1.000	5.812	2.411
	1500	1500.312	0.312	1.152	1.000	2.401	1.550
0.5	500	500.524	0.524	1.711	1.000	5.053	2.248
	1000	1000.295	0.295	1.328	1.000	3.114	1.765
	1500	1500.147	0.147	0.825	0.000	1.09	1.053
0.75	500	500.683	0.668	2.259	1.500	8.241	2.871
	1000	1000.511	0.511	1.815	1.000	5.608	2.368
	1500	1500.284	0.284	1.191	1.000	2.513	1.585
0.95	500	501.214	1.214	3.012	2.000	14.502	3.808
	1000	1001.054	1.054	2.603	2.000	11.217	3.349
	1500	1500.518	0.518	1.683	1.000	5.106	2.260
Sample size ($n=5000$ )
0.05	1250	1250.72	0.72	2.11	1.0	11.42	3.38
	2500	2500.44	0.44	1.85	1.0	7.41	2.72
	3750	3750.31	0.31	1.40	1.0	4.85	2.20
0.25	1250	1259.40	9.40	12.30	9.0	385.12	19.62
	2500	2502.11	2.11	6.40	5.0	112.40	10.60
	3750	3755.80	5.80	10.25	8.0	211.20	14.53
0.50	1250	1584.12	334.12	482.11	412.00	384112.50	619.76
	2500	2500.82	0.82	1.42	1.0	4.92	2.22
	3750	4110.45	360.45	520.14	448.0	481512.90	693.91
0.75	1250	1258.10	8.10	10.45	8.0	290.41	17.04
	2500	2505.44	5.44	8.12	6.0	185.10	13.61
	3750	3754.91	4.91	8.44	6.0	154.12	12.41
0.95	1250	1250.88	0.88	2.25	2.0	12.96	3.60
	2500	2508.33	8.33	11.20	9.0	312.44	17.68
	3750	3750.54	0.54	1.62	1.0	5.91	2.43

ARMA(1,1)-GJR-GARCH(1,1) model with Student-t error distribution

In Tables 3-5, the simulation results obtained using a bidirectional BiLSTM deep neural network architecture combined with Deep Quantile Regression (DQR) for change point detection show that the estimated change points ($k$ ) are significantly different from the observed true change points ($k_0$ ). This discrepancy is often caused by specific structural and statistical phenomena. Financial return series tend to be highly chaotic and exhibit volatility clustering, which refers to periods of sustained high or low variance. A bidirectional network processes information from both ends of a time window, which can lead to misinterpretation of temporary, localized clusters of extreme conditional variance or sequences of outliers as permanent structural change points. If an ARMA-GARCH shock or outlier occurs near the boundaries of the time series, the resulting asymmetry can cause spikes in the network’s error metrics at the edges. Consequently, when the global pinball loss optimization seeks to minimize total error, it may mistakenly drop the estimated change point $k$ near the boundaries of the sequence. This results in an estimate that is markedly distorted from an observed change point that is actually situated in the middle of the timeline. To address these challenges, we propose incorporating additional deep neural network architectures in our next article to explore how different models can mitigate the weaknesses of the bidirectional BiLSTM architecture.

Figure 9. Simulated ARMA (1,1)-ARCH (1,1) process for a sample size of 5000 and corresponding histograms at change points at t = 1250 and t =2500. and t =3750.

Figure 9 illustrates the simulated ARMA(1,1)-GARCH(1,1) model based on a sample of 5000 observations, segmented into four regimes at the actual structural breaks occurring at t = 1250, t = 2500, and t = 3750. High volatility clustering is evident in the second and fourth regimes of the simulated data, characterised by sharp peaks. In contrast, the first and third regimes display low volatility movements, remaining relatively calm and stable in terms of fluctuations. The second part of the plot features a histogram of the estimated change points for the simulated ARMA (1,1)-GARCH (1,1) process.

Table 4. Simulation results on consistency of the deep quantile regression multiple change point estimator using ARMA(1,1)-GARCH(1,1) model.

Sample size ($n=2000$ )
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	500	501.24	1.24	3.82	2.0	31.14	5.58
	1000	1000.62	0.62	2.14	1.0	10.12	3.18
	1500	1500.51	0.51	2.05	1.0	9.42	3.07
0.25	500	509.41	9.41	14.52	11.0	412.30	20.31
	1000	100.94	0.94	1.82	1.0	6.18	2.49
	1500	1506.18	6.18	11.45	9.0	284.10	16.86
0.50	500	712.18	212.18	392.40	314.0	215410.20	464.12
	1000	1000.94	0.94	1.82	1.0	6.18	2.49
	1500	1341.12	−158.88	412.33	342.0	291402.50	539.82
0.75	500	508.62	8.62	12.14	9.0	321.44	17.93
	1000	1005.12	5.12	9.40	7.0	210.45	14.51
	1500	1505.20	5.20	9.18	7.0	192.40	13.87
0.95	500	501.55	1.55	4.10	3.0	38.25	6.18
	1000	1009.44	9.44	14.12	11.0	402.18	20.05
	1500	1500.82	0.82	2.40	1.0	12.18	3.49
Sample size ($n=5000$ )
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	1250	1248.51	−1.49	6.12	5.00	62.45	7.90
	2500	2505.34	5.34	19.12	15.00	572.40	23.92
	3750	3748.21	−1.79	9.14	7.00	128.50	11.34
0.25	1250	1249.12	−0.88	5.04	4.00	41.80	6.47
	2500	2511.40	11.40	32.18	27.00	1422.30	37.71
	3750	3747.65	−2.35	12.44	10.00	225.80	15.03
0.50	1250	1249.65	−0.35	3.12	2.00	18.22	4.27
	2500	2506.11	6.11	24.89	20.00	884.22	29.74
	3750	3746.88	−3.12	14.10	11.00	288.42	16.98
0.75	1250	1249.34	−0.66	4.88	4.00	39.10	6.25
	2500	2506.11	6.11	24.89	20.00	884.22	29.74
	3750	3747.92	−2.08	11.85	9.00	204.60	14.30
0.95	1250	1249.58	−0.42	5.92	5.00	58.71	7.66
	2500	2502.43	2.43	17.52	13.00	496.11	22.27
	3750	3749.03	−0.97	8.81	7.00	118.25	10.87

Table 5. Simulation results on consistency of the deep quantile regression multiple change point estimator using ARMA(1,1)-GJR-GARCH(1,1) model with Student-t error distribution.

Sample size ($n=2000$ )
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	500	498.42	−1.58	5.82	4.00	56.12	7.49
	1000	1001.95	1.95	10.42	8.00	178.50	13.36
	1500	1498.12	−1.88	7.95	6.00	98.40	9.92
0.25	500	499.18	−0.82	4.12	3.00	31.45	5.61
	1000	1004.12	4.12	16.85	13.00	412.33	20.31
	1500	1497.35	−2.65	10.42	8.00	162.15	12.73
0.50	500	499.74	−0.26	2.54	2.00	11.20	3.35
	1000	1007.84	7.84	22.10	18.00	722.14	26.82
	1500	1492.42	−3.58	12.18	10.00	224.60	14.99
0.75	500	499.31	−0.69	4.35	3.00	33.18	5.76
	1000	1005.11	5.11	18.24	14.00	498.60	22.33
	1500	1497.64	−2.36	10.15	8.00	154.22	12.42
0.95	500	498.92	−1.08	5.44	4.00	48.91	6.99
	1000	1003.18	3.18	13.56	10.00	284.15	16.86
	1500	1498.54	−1.46	7.42	6.00	87.14	9.33
Sample size ($n=5000$ )
Quantile (τ)	True CP	Est. CP	Bias	MAE	MAD	MSE	RMSE
0.05	1250	1247.62	−2.38	7.94	6.00	98.42	9.92
	2500	2502.84	2.84	14.12	11.00	324.50
	3750	3747.88	−2.12	10.65	8.00	182.14	13.50
0.25	1250	1248.81	−1.19	5.85	5.00	54.10	7.36
	2500	2506.15	6.15	22.40	18.00	711.20	26.67
	3750	3746.51	−3.49	15.12	12.00	344.20	18.55
0.50	1250	1249.68	−0.32	3.65	3.00	22.84	4.78
	2500	2514.92	14.92	38.65	331.00	2104.18	45.87
	3750	3744.18	−5.82	21.40	17.00	682.15	26.12
0.75	1250	1249.11	−0.89	5.12	4.00	43.15	6.57
	2500	2508.31	8.31	25.18	21.00	915.60	30.26
	3750	3747.12	−2.88	13.84	11.00	291.60	17.08
0.95	1250	1248.34	−1.66	7.21	6.00	85.34	9.24
	2500	2504.12	4.12	18.25	15.00	488.94	22.11
	3750	3748.42	−1.58	9.94	8.00	154.20	12.42

3.4. Convergence Rate

The asymptotic convergence rate of a deep quantile change-point estimator depends on the magnitude of the structural break, the complexity of the network architecture, and the intrinsic dimension of the data. To implement a Monte Carlo study of the asymptotic convergence rate of the Deep Quantile Regression (DQR) estimator for change-point detection in an ARIMA(1,1)-GARCH(1,1) process, multiple time-series trajectories for different sample sizes $n$ were generated. The deep quantile regression neural network was fitted to estimate the conditional quantile, and the change-point location $\tau$ was then determined for each dataset to empirically calculate the convergence rate. To assess the estimator’s asymptotic convergence rate, the estimator error must satisfy the typical parametric rate ${\mathcal{O}}_p\left(n^{-1}\right)$ or the localized structural rate ${\mathcal{O}}_p\left(1\right)$ . This was tested by varying sample sizes $n\in \left\{250,500,1000,2000\right\}$ across iterations. To determine the empirical convergence speed, the results were log-transformed and the regression slope was extracted from

$\text{log}\left(\text{MSE}\right)={\beta }_0+{\beta }_1\log \left(n\right)$

Table 6 presents the relative localization error matrix across different quantiles and sample sizes.

Figure 10 presents the change-point estimator convergence profile simulation results that confirm the deep quantile regression change-point estimator achieves a parametric convergence rate of ${\mathcal{O}}_p\left(n^{-1}\right)$ under stable, regular conditions. This means the estimation error scales inversely with the sample size, making it highly accurate for identifying structural breaks in the conditional distribution of complex time-series data.

Table 6. Relative localization error matrix ($\Delta t/n$ ) across quantiles and sample sizes.

Sample Size (n)	$\tau =0.05$	$\tau =0.25$	$\tau =0.50$	$\tau =0.75$	$\tau =0.90$	Combined Multi-Loss
1000	13.42%	11.15%	9.80%	11.08%	13.91%	11.35%
2000	8.85%	6.94%	5.92%	7.12%	8.71%	7.20%
4000	4.61%	3.25%	2.68%	3.40%	4.88%	3.54%
Empirical Rate (α)	0.77	0.89	0.93	0.85	0.76	0.84

Figure 10. Change-point estimator convergence profile (Max $N=4000$ ).

4. Conclusion

This article derived the statistical asymptotic properties of the deep quantile regression estimator for change point detection. Specifically, the study proves both the consistency and the rate of convergence of the estimator, providing a solid theoretical foundation for a field that has largely been driven by empirical research. These properties ensure that as the sample size increases and the depth and width of the network are scaled appropriately, the estimated quantile curves converge to the underlying conditional quantiles. The established asymptotic properties bridge the gap between empirical deep learning applications and non-parametric high-dimensional statistical inference. Consequently, the deep quantile regression estimator can be applied with greater confidence, using the derived asymptotic properties to construct confidence intervals and perform hypothesis testing on complex, non-linear data structures. Future research should extend these asymptotic properties to non-independent and identically distributed mixing processes to accommodate heavy-tailed, dependent time-series data. Further, there is a need to explore non-asymptotic boundary cases through finite-sample error bounds.