Change-Point Detection in the Volatility of Conditional Heteroscedastic Autoregressive Nonlinear Models

Abstract

This paper studies single change-point detection in the volatility of a class of parametric conditional heteroscedastic autoregressive nonlinear (CHARN) models. The conditional least-squares (CLS) estimators of the parameters are defined and are proved to be consistent. A Kolmogorov–Smirnov type-test for change-point detection is constructed and its null distribution is provided. An estimator of the change-point location is defined. Its consistency and its limiting distribution are studied in detail. A simulation experiment is carried out to assess the performance of the results, which are compared to recent results and applied to two sets of real data.

1. Introduction

Detecting jumps in a series of real numbers and determining their number and locations is known in statistics as a change-point problem. This is usually solved by testing for the stationarity of the series and estimating the change locations when the null hypothesis of stationarity is rejected.

Change-point problems can be encountered in a wide range of disciplines, such as quality control, genetic data analysis and bioinformatics (see, e.g., [1,2]) and financial analysis (see, e.g., [3,4,5]). The study of conditional variations in financial and economic data receives particular attention as a result of its interest in hedging strategies and risk management.

The literature on change-points is vast. Parametric or non-parametric approaches are used for independent and identically distributed (iid) data as well as for dependent data. The pioneering works of [6,7] proposed tests for identifying a deviation in the mean of iid Gaussian variables from industrial quality control. Ref. [8] proposed a test for detecting changes in the mean. This criterion was later generalized by [9], which envisioned a much more general model allowing incremental changes. Refs. [10,11,12] presented a documentary analysis of several non-parametric procedures.

A popular alternative to using the likelihood ratio test was employed in [13,14]. Ref. [15] reviewed the asymptotic behavior of likelihood ratio statistics for testing a change in the mean in a series of iid Gaussian random variables. Ref. [16] came up with statistics based on linear rank statistical processes with quantum scores. Ref. [17] looked at detection tests and change-point estimation methods for models based on the normal distribution. The contribution of [18] is related to the change in mean and variance. Ref. [19] proposed permutation tests for the location and scale parameters of a law. Ref. [20] developed a change-point test using the empirical characteristic functions. Ref. [21] proposed several CUSUM approaches. Ref. [22] proposed tests for change detection in the mean, variance, and autoregressive parameters of a p-order autoregressive model. Ref. [23] used a weighted CUSUM procedure to identify a potential change in the mean and covariance structure of linear processes.

There are several types of changes depending on the temporal behavior of the series studied. The usual ones are abrupt change, gradual change, and intermittent change. These are studied in the frameworks of on-line or off-line data. Some papers dealing with this problem are [24,25]. In this paper, we focus on abrupt change in the conditional variance of off-line data issue from a class of CHARN models (see [26,27]). These models are of the most famous and significant ones in finance, which include many financial time series models. We suggest a hybrid estimation procedure, which combines CLS and non-parametric methods to estimate the change location. Indeed, conditional least-squares estimators own a computational advantage and require no knowledge of the innovation process.

The rest of the paper is organized as follows. Section 2 presents the class of models studied, the notation, the main assumptions and the main result on the CLS estimators of the parameters. Section 3 presents the change-point test and the change location LS estimation. The asymptotic distribution of the test statistic under the null hypothesis is investigated. The consistency rates are obtained for the change location estimator and its limit distribution is derived. Section 4 presents the simulation results from a few simple time series models. Here, our results are compared to some recent methods. They are also applied to two real data sets. Our work ends with a conclusion in Section 5. The proofs and auxiliary results are given in Appendix A.

2. Model and Assumptions

2.1. Notation and Assumptions

Let l and r be positive integers. For given real functions $\mathcal{A}(\alpha ;z)$ defined on a non-empty subset of ${\mathbb{R}}^{\ell }\times {\mathbb{R}}^p$, $\ell =l,r$ and $\mathcal{K}(\psi ;z)$ defined on a non-empty subset of ${\mathbb{R}}^r\times {\mathbb{R}}^l\times {\mathbb{R}}^p$, $\alpha =\left({\alpha }_1,\dots ,{\alpha }_{\ell }\right)$, $\psi =\left({\rho }^{⊺},{\theta }^{⊺}\right)$, ${\rho }^{⊺}=\left({\rho }_1,\dots ,{\rho }_r\right)$, ${\theta }^{⊺}=\left({\theta }_1,\dots ,{\theta }_l\right)$. We denote:

$${\partial }_{\alpha }\mathcal{A}(\alpha ;z)={\left({\displaystyle \frac{\partial \mathcal{A}}{\partial {\alpha }_1}}(\alpha ;z),\dots ,{\displaystyle \frac{\partial \mathcal{A}}{\partial {\alpha }_{\ell }}}(\alpha ;z)\right)}^{⊺}$$

$${\partial }_{{\alpha }^2}^2\mathcal{A}(\alpha ;z)=\left({\displaystyle \frac{{\partial }^2\mathcal{A}}{\partial {\alpha }_i\partial {\alpha }_j}}(\alpha ;z);1⩽i,j⩽\ell \right)$$

$${\partial }_{\rho }\mathcal{K}(\psi ;z)={\left({\displaystyle \frac{\partial \mathcal{K}}{\partial {\rho }_1}}(\psi ;z),\dots ,{\displaystyle \frac{\partial \mathcal{K}}{\partial {\rho }_r}}(\psi ;z)\right)}^{⊺}$$

$${\partial }_{\theta }\mathcal{K}(\psi ;z)={\left({\displaystyle \frac{\partial \mathcal{K}}{\partial {\theta }_1}}(\psi ;z),\dots ,{\displaystyle \frac{\partial \mathcal{K}}{\partial {\theta }_l}}(\psi ;z)\right)}^{⊺}$$

$${\partial }_{\rho \theta }^2\mathcal{K}(\psi ;z)=\left({\displaystyle \frac{{\partial }^2\mathcal{K}}{\partial {\rho }_i\partial {\theta }_j}}(\psi ;z);1⩽i⩽r,1⩽j⩽l\right)$$

$${\partial }_{\theta \rho }^2\mathcal{K}(\psi ;z)=\left({\displaystyle \frac{{\partial }^2\mathcal{K}}{\partial {\theta }_i\partial {\rho }_j}}(\psi ;z);1⩽i⩽l,1⩽j⩽r\right)$$

$${\partial }_{{\rho }^2}^2\mathcal{K}(\psi ;z)=\left({\displaystyle \frac{{\partial }^2\mathcal{K}}{\partial {\rho }_i\partial {\rho }_j}}(\psi ;z);1⩽i,j⩽r\right)$$

$${\partial }_{{\theta }^2}^2\mathcal{K}(\psi ;z)=\left({\displaystyle \frac{{\partial }^2\mathcal{K}}{\partial {\theta }_i\partial {\theta }_j}}(\psi ;z);1⩽i,j⩽l\right).$$

For a vector or matrix function $\zeta \left(x\right)$, we denote by ${\partial }^{⊺}\zeta \left(x\right)$, the transpose of $\partial \zeta \left(x\right)$. We define:

$$\partial \mathcal{K}(\psi ;z)=\left(\begin{array}{c}{\partial }_{\rho }\mathcal{K}(\psi ;z)\\ {\partial }_{\theta }\mathcal{K}(\psi ;z)\end{array}\right)\mathrm{and}{\partial }^2\mathcal{K}(\psi ;z)=\left(\begin{array}{cc}{\partial }_{{\rho }^2}^2\mathcal{K}(\psi ;z)& {\partial }_{\rho \theta }^2\mathcal{K}(\psi ;z)\\ {\partial }_{\theta \rho }^2\mathcal{K}(\psi ;z)& {\partial }_{{\theta }^2}^2\mathcal{K}(\psi ;z)\end{array}\right).$$

All along the text, the notations $\stackrel{d}{⟶}$ and $\stackrel{\mathcal{D}}{⟶}$ denote, respectively, the weak convergence in functional spaces and the convergence in distribution.

We place ourselves in the framework where the observations at hand are assumed to be issued from the following CHARN ($p,p$) model:

$$X_t=m\left(\rho ;Z_{t-1}\right)+\sigma \left(\theta ;Z_{t-1}\right){\epsilon }_t,t\in \mathbb{Z},$$

(1)

where $p\in {\mathbb{N}}^{\ast }\bigcup \left\{\infty \right\}$; $m(·)$ and $\sigma (·)$ are two real-valued functions of known forms depending on unknown parameters $\rho$ and $\theta$, respectively; for all $t\in \mathbb{Z}$, $Z_{t-1}={\left(X_{t-1},X_{t-2},\dots ,X_{t-p}\right)}^{⊺}$; ${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$ is a sequence of stationary random variables with $\mathbb{E}\left({\epsilon }_t\mid Z_{t-1}\right)=0$ and $\mathbb{V}\mathrm{ar}\left({\epsilon }_t\mid Z_{t-1}\right)=1$ such that ${\epsilon }_t$ is independent of the $\sigma –$algebra ${\mathcal{F}}_{t-1}=\sigma \left(Z_k,k<t\right)$. The case $p=\infty$ is treated in [28,29,30] where the stationarity and the ergodicity of the process ${\left(X_t\right)}_{t\in \mathbb{Z}}$ is studied. Although we restrict to $p<\infty$, all the results stated here also hold for $p=\infty$.

Let $\psi ={\left({\rho }^{⊺},{\theta }^{⊺}\right)}^{⊺}\in \mathsf{\Psi }=\mathrm{int}(\Theta )\times \mathrm{int}\left(\tilde{\Theta }\right)\subset {\mathbb{R}}^r\times {\mathbb{R}}^l$, the vector of the parameters of the model (1) and ${\psi }_0={\left({\rho }_0^{⊺},{\theta }_0^{⊺}\right)}^{⊺}$ the true parameter vector. Denote by $∥M∥$ an appropriate norm of a vector or a matrix M. We assume that all the random variables in the whole text are defined on the same probability space $\left(\Omega ,\mathcal{F},\mathbb{P}\right)$. We make the following assumptions:

(A₁)

The common fourth order moment of the ${\epsilon }_t$ is finite.

(A₂)

• The function $m(·)$ is twice continuously differentiable, a.e., with respect to $\rho$ in some neighborhood ${\mathcal{B}}_1$ of ${\rho }_0$.
• The function $\sigma (.)$ is twice continuously differentiable, a.e., with respect to $\theta$ in some neighborhood ${\mathcal{B}}_2$ of ${\theta }_0$.
• There exists a positive function $\omega$ such that $\mathbb{E}\left({\omega }^4\left(Z_0\right)\right)<\infty$, and

  $$max\{\underset{\rho \in \mathrm{int}(\Theta )}{sup}\left|m(\rho ;z)\right|,\underset{\rho \in \mathrm{int}(\Theta )}{sup}∥{\partial }_{\rho }m(\rho ;z)∥,\underset{\rho \in \mathrm{int}(\Theta )}{sup}∥{\partial }_{{\rho }^2}^{2}m(\rho ;z)∥\}⩽\omega \left(z\right)$$

  $$max\left\{\underset{\theta \in \mathrm{int}\left(\tilde{\Theta }\right)}{sup}\left|\sigma (\theta ;z)\right|,\underset{\theta \in \mathrm{int}\left(\tilde{\Theta }\right)}{sup}∥{\partial }_{\theta }\sigma (\theta ;z)∥,\underset{\theta \in \mathrm{int}\left(\tilde{\Theta }\right)}{sup}∥{\partial }_{{\theta }^2}^2\sigma (\theta ;z)∥\right\}⩽\omega \left(z\right).$$

(A₃)

There exists a positive function $\beta$ such that $\mathbb{E}\left({\beta }^4\left(Z_0\right)\right)<\infty$, and for all ${\rho }_1,{\rho }_2\in \mathrm{int}(\Theta )$, and ${\theta }_1,{\theta }_2\in \mathrm{int}(\tilde{\Theta })$,

$$\begin{array}{cc}\hfill max& \{\left|m\left({\rho }_1;z\right)-m\left({\rho }_2;z\right)\right|,∥{\partial }_{\rho }m\left({\rho }_1;z\right)-{\partial }_{\rho }m\left({\rho }_2;z\right)∥,\hfill \\ \hfill & ∥{\partial }_{{\rho }^2}^{2}m\left({\rho }_1;z\right)-{\partial }_{{\rho }^2}^{2}m\left({\rho }_2;z\right)∥,\left|\sigma \left({\theta }_1;z\right)-\sigma \left({\theta }_2;z\right)\right|,\hfill \\ \hfill & ∥{\partial }_{\theta }\sigma \left({\theta }_1;z\right)-{\partial }_{\theta }\sigma \left({\theta }_2;z\right)∥,∥{\partial }_{{\theta }^2}^2\sigma \left({\theta }_1;z\right)-{\partial }_{{\theta }^2}^2\sigma \left({\theta }_2;z\right)∥\}\hfill \\ \hfill & ⩽\beta \left(z\right)min\left\{{∥{\rho }_1-{\rho }_2∥}_2,{∥{\theta }_1-{\theta }_2∥}_2\right\}.\hfill \end{array}$$

(A₄)

The sequence ${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$ is stationary and satisfies either of the following two conditions:

• $\alpha$-mixing with mixing coefficient satisfying ${\sum }_{n\ge 1}{\left[\alpha \left(n\right)\right]}^{\delta /(2+\delta )}<\infty$ and $\mathbb{E}|{\epsilon }_0{|}^{2+\delta }<\infty$ for some $\delta >0$;
• $\varphi$-mixing with mixing coefficient satisfying ${\sum }_{n⩾1}^{}{\left[\varphi \left(n\right)\right]}^{1/2}<\infty$ and $\mathbb{E}{\left|{\epsilon }_0\right|}^{4+\delta }<\infty$ for some $\delta >0$.

2.2. Parameter Estimation

2.2.1. Conditional Least-Squares Estimation

The conditional mean and the conditional variance of $X_t$ are given, respectively, by $\mathbb{E}\left(X_t\mid {\mathcal{F}}_{t-1}\right)$=$m\left(\rho ;Z_{t-1}\right)$ and $\mathbb{V}\mathrm{ar}\left(X_t\mid {\mathcal{F}}_{t-1}\right)={\sigma }^2\left(\theta ;Z_{t-1}\right)$. From these, one has that for all $z\in {\mathbb{R}}^p$,

$$\mathbb{E}\left(X_1\mid Z_0=z\right)=m\left(\rho ;z\right)\mathrm{and}\mathbb{E}\left({\left(X_1-m\left(\rho ;Z_0\right)\right)}^2\mid Z_0=z\right)={\sigma }^2\left(\theta ;z\right).$$

Therefore, for any bounded measurable functions $g(·)$ and $k(·)$, we have

$$\mathbb{E}\left\{\left[X_1-m\left(\rho ;Z_0\right)\right]g\left(Z_0\right)\right\}=0\mathrm{and}\mathbb{E}\left(\left\{{\left[X_1-m\left(\rho ;Z_0\right)\right]}^2-{\sigma }^2\left(\theta ;Z_0\right)\right\}k\left(Z_0\right)\right)=0.$$

Without a loss of generality, in the following we take, for all $z\in {\mathbb{R}}^p,g\left(z\right)=k\left(z\right)=1$. Now, given $X_{-p+1}, \dots , X_{-1}, X_0, X_1, \dots , X_n$ with $n\gg p$, we let ${\mathbb{X}}_n=(X_{-p+1},$ $\dots ,$ $X_{-1},$ $X_0,$ $X_1,$ $\dots ,$ $X_n)$ and consider the sequences of random functions

$$\begin{array}{ccc}\hfill Q_n\left(\rho \right)=Q_n(\rho ;{\mathbb{X}}_n)& =& \sum _{t=1}^n{\left(X_t-\mathbb{E}\left(X_t\mid {\mathcal{F}}_{t-1}\right)\right)}^2=\sum _{t=1}^n{\left(X_t-m\left(\rho ;Z_{t-1}\right)\right)}^2\hfill \\ \hfill S_n(\rho ,\theta )=S_n(\rho ,\theta ;{\mathbb{X}}_n)& =& \sum _{t=1}^n{\left({\left(X_t-m\left(\rho ;Z_{t-1}\right)\right)}^2-{\sigma }^2\left(\theta ;Z_{t-1}\right)\right)}^2.\hfill \end{array}$$

We have the following theorem:

Theorem 1.

Under assumptions $\left(A_1\right)$–$\left(A_3\right)$, there exists a sequence of estimators ${\widehat{\psi }}_n={\left({\widehat{\rho }}_n^{⊺},{\widehat{\theta }}_n^{⊺}\right)}^{⊺}$ such that ${\widehat{\psi }}_n⟶{\psi }_0$ almost surely, and for any $ϵ>0$, there exists an event E with $\mathbb{P}\left(E\right)>1-ϵ$, and a non-negative integer $n_0$ such that on E, for $n>n_0$,

• ${\displaystyle \frac{\partial Q_n}{\partial \rho }}\left({\widehat{\rho }}_n;{\mathbb{X}}_n\right)=0$ and $Q_n(\rho ;{\mathbb{X}}_n)$ attains a relative minimum at $\rho ={\widehat{\rho }}_n$;
• assuming ${\widehat{\rho }}_n$ fixed, ${\displaystyle \frac{\partial S_n}{\partial \theta }}\left({\widehat{\psi }}_n;{\mathbb{X}}_n\right)=0$ and $S_n\left(\left({\widehat{\rho }}_n,\theta \right);{\mathbb{X}}_n\right)$ attains a relative minimum at $\theta ={\widehat{\theta }}_n$.

Proof. 

This result is an extension of [31] to the case ${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$ is a mixing martingale difference. The proof can be handled in the same lines and is left to the reader.    □

3. Change-Point Study

3.1. Change-Point Test and Change Location Estimation

We essentially use the techniques of [32], who studied the estimation of the shift in the mean of a linear process by a LS method. We first consider the model (1) for known $\rho$, and $\sigma \left(\theta ;Z_{t-1}\right)=\underline{\theta }{\delta }_0\left(Z_{t-1}\right)$, for some known positive real-valued function ${\delta }_0(·)$ defined on ${\mathbb{R}}^p$ and for an unknown positive real number $\underline{\theta }$. We wish to test

$$H_0:\underline{\theta }={\vartheta }_1={\vartheta }_2\mathrm{over}t⩽n$$

against

$$H_1:\underline{\theta }=\left\{\begin{array}{cc}{\vartheta }_1,\hfill & t=1,\dots ,t^{\ast }\hfill \\ {\vartheta }_2,\hfill & t=t^{\ast }+1,\dots ,n.\hfill \end{array}\right.\left({\vartheta }_1\ne {\vartheta }_2\right)$$

where ${\vartheta }_1,{\vartheta }_2$ and $t^{\ast }$ are unknown parameters.

We are also interested in estimating ${\vartheta }_1,{\vartheta }_2$ and the change location $t^{\ast }$, when $H_0$ is rejected. It is assumed that $t^{\ast }=\left[n\tau \right]$ for some $\tau \in \left(0,1\right)$, with $\left[x\right]$ standing for the integer part of any real number x. From (1), one can easily check that

$${\left(X_t-m\left(\rho ;Z_{t-1}\right)\right)}^2={\delta }^2\left(Z_{t-1}\right)+{\delta }^2\left(Z_{t-1}\right)\left({\epsilon }_t^2-1\right),t\in \mathbb{Z}$$

(2)

from which we define the LS estimator ${\widehat{t}}^{\ast }$ of $t^{\ast }$ as follows:

$${\widehat{t}}^{\ast }:=\underset{1⩽k<n}{arg\; min}\left[\underset{{\vartheta }_1,{\vartheta }_2}{min}\left\{\sum _{t=1}^k{\left(W_t^2-{\vartheta }_1^2\right)}^2+\sum _{t=k+1}^n{\left(W_t^2-{\vartheta }_2^2\right)}^2\right\}\right],$$

(3)

where $W_t=\left(X_t-m(\rho ;Z_{t-1})\right)/{\delta }_0\left(Z_{t-1}\right)$. Thus, the change location is estimated by minimizing the sum of squares of residuals among all possible sample slits.

Letting

$${\overline{W}}_k={\displaystyle \frac{1}{k}}\sum _{t=1}^k{W}_t^2,{\overline{W}}_{n-k}={\displaystyle \frac{1}{n-k}}\sum _{t=k+1}^n{W}_t^2\mathrm{and}\overline{W}={\displaystyle \frac{1}{n}}\sum _{t=1}^n{W}_t^2,$$

it is easily seen that for some k, the LS estimator of ${\vartheta }_1^2\left(t⩽k\right)$ and ${\vartheta }_2^2\left(t>k\right)$ are ${\overline{W}}_k$ and ${\overline{W}}_{n-k}$, respectively, and that (3) can be written as

$$\begin{array}{ccc}{\widehat{t}}^{\ast }\hfill & =\hfill & {\displaystyle \underset{1⩽k<n}{arg\; min}\left\{\sum _{t=1}^k{\left(W_t^2-{\overline{W}}_k\right)}^2+\sum _{t=k+1}^n{\left(W_t^2-{\overline{W}}_{n-k}\right)}^2\right\}}\hfill \\ \hfill & =\hfill & {\displaystyle \underset{1⩽k<n}{arg\; min}{S}_k^2.}\hfill \end{array}$$

(4)

Let ${\displaystyle S^2=\sum _{t=1}^n{\left(W_t^2-\overline{W}\right)}^2}$. A simple algebra gives

$$S^2=S_k^2+U_k,$$

(5)

where

$$U_k=k{\left({\overline{W}}_k-\overline{W}\right)}^2+\left(n-k\right){\left({\overline{W}}_{n-k}-\overline{W}\right)}^2.$$

(6)

From (4) and (5), we have

$$\begin{array}{ccc}{\widehat{t}}^{\ast }\hfill & =\hfill & {\displaystyle \underset{1⩽k<n}{arg\; min}\left(S^2-U_k\right)}\hfill \\ \hfill & =\hfill & {\displaystyle \underset{1⩽k<n}{arg\; max}{U}_k.}\hfill \end{array}$$

(7)

From (6), a simple algebraic computation gives the following alternative expression for $U_k$:

$$\begin{array}{ccc}{U}_k\hfill & =\hfill & {\displaystyle {\displaystyle \frac{n}{k(n-k)}}{\left(\sum _{t=1}^k\left(W_t^2-\overline{W}\right)\right)}^2}\hfill \\ \hfill & =\hfill & {\displaystyle {\left(\sqrt{{\displaystyle \frac{n}{k(n-k)}}}\sum _{t=1}^k\left(W_t^2-\overline{W}\right)\right)}^2}\hfill \\ \hfill & =\hfill & T_k^2.\hfill \end{array}$$

(8)

It results from (7) and (8) that

$$\begin{array}{ccc}\hfill {\widehat{t}}^{\ast }& =& \underset{1⩽k<n}{arg\; max}{T}_k^2\hfill \\ & =& \underset{1⩽k<n}{arg\; max}\left|T_k\right|.\hfill \end{array}$$

(9)

Writing $T_k^2=n{\Delta }_k^2$, it is immediate that

$$\begin{array}{ccc}\hfill {\Delta }_k^2& =& {\displaystyle \frac{1}{k(n-k)}}{\left(\sum _{t=1}^k\left(W_t^2-\overline{W}\right)\right)}^2\hfill \\ & =& {\displaystyle \frac{k}{(n-k)}}{\left({\overline{W}}_k-\overline{W}\right)}^2.\hfill \end{array}$$

Simple computations give

$${\Delta }_k^2={\displaystyle \frac{k(n-k)}{n^2}}{\left({\overline{W}}_{n-k}-{\overline{W}}_k\right)}^2,$$

from which we have

$${\widehat{t}}^{\ast }=\underset{1⩽k<n}{arg\; max}{\Delta }_k^2=\underset{1⩽k<n}{arg\; max}\left|{\Delta }_k\right|.$$

(10)

The test statistic we use for testing $H_0$ against $H_1$ is a scale version of ${\displaystyle \underset{1⩽k⩽n-1}{max}}\left|T_k\right|$.

One can observe that under some conditions (e.g., ${\epsilon }_t$ i.i.d. with ${\epsilon }_t\sim \mathcal{N}(0,1)$), this statistic is the equivalent likelihood based test statistic for testing $H_0$ against $H_1$ (see, e.g., [33]).

Let

$$C_k=\sum _{t=1}^k{W}_t^2,C_{n-k}=\sum _{t=k+1}^n{W}_t^2\mathrm{and}{C}_n=\sum _{t=1}^n{W}_t^2.$$

(11)

By simple calculations, we obtain

$$\begin{array}{ccc}{T}_k\hfill & =\hfill & {\displaystyle \sqrt{{\displaystyle \frac{n}{k(n-k)}}}\sum _{t=1}^k\left(W_t^2-\overline{W}\right)}\hfill \\ \hfill & =\hfill & {\displaystyle {\left(q\left(\frac{k}{n}\right)\right)}^{-1}\left({\displaystyle \frac{1}{\sqrt{n}}}\left(C_k-{\displaystyle \frac{k}{n}}{C}_n\right)\right),}\hfill \end{array}$$

(12)

where $q(·)$ is a positive weight function defined for any $x\in \left(0,1\right)$ by $q\left(x\right)=\sqrt{x\left(1-x\right)}$.

3.2. Asymptotics

3.2.1. Asymptotic Distribution of the Test Statistic

The study of the asymptotic distribution of the test statistic under $H_0$, is based on that of the the process ${\xi }_n(·)$ defined for any $s\in [0,1]$ by

$${\xi }_n\left(s\right)=C_n\left(s\right)-s{C}_n\left(1\right),$$

(13)

where

$$C_n\left(s\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}0⩽s<{\displaystyle \frac{1}{n}}\mathrm{and}1-{\displaystyle \frac{1}{n}}<s<1\hfill \\ {\displaystyle \sum _{t=1}^{\left[ns\right]}{W}_t^2}\hfill & \mathrm{if}{\displaystyle \frac{1}{n}}⩽s⩽1-{\displaystyle \frac{1}{n}}\hfill \\ {\displaystyle \sum _{t=1}^n{W}_t^2}\hfill & \mathrm{if}s=1,\hfill \end{array}\right.$$

(14)

where we recall that $\left[ns\right]$ is the integer part of $ns$. For some $\delta \in (1/n,1/2)$ and for any s in $\left[\delta ,1-\delta \right]$, we define

$$\begin{array}{c}\hfill T_n\left(s\right)={\displaystyle \frac{{\xi }_n\left(s\right)}{\sqrt{n}q\left(s\right)}}\mathrm{and}{\Lambda }_n=\underset{\delta ⩽s⩽1-\delta }{max}\frac{|T_n\left(s\right)|}{{\widehat{\sigma }}_w},\end{array}$$

(15)

where $q\left(s\right)=\sqrt{s(1-s)}$ and ${\widehat{\sigma }}_w$ is any consistent estimator of

$${\displaystyle {\sigma }_w^2=\mathbb{E}{\left(W_1^2-\mathbb{E}\left(W_1^2\right)\right)}^2+2\sum _{t⩾2}^{}\mathbb{E}\left(\left(W_1^2-\mathbb{E}\left(W_1^2\right)\right)\left(W_t^2-\mathbb{E}\left(W_t^2\right)\right)\right).}$$

For $\delta \in (0,1/2)$, we denote by $D_{\delta }\equiv D\left(\left[\delta ,1-\delta \right]\right)$ the space of all right continuous functions with left limits on $\left[\delta ,1-\delta \right]$ endowed with the Skorohod metric. It is clear that $C_n(·),{\xi }_n(·)\in D_0$ and $T_n(·)\in D_{\delta }$.

Theorem 2.

Assume that the assumptions $\left(A_1\right)$–$\left(A_4\right)$ hold. Then, under $H_0$, we have

1. 

${\displaystyle \frac{{\xi }_n\left(s\right)}{{\sigma }_w\sqrt{n}}}\stackrel{d}{⟶}\tilde{B}\left(s\right)\mathrm{in}{D}_0\mathrm{as}n⟶\infty$;

2. 

${\Lambda }_n\stackrel{\mathcal{D}}{⟶}{\displaystyle \underset{\delta ⩽s⩽1-\delta }{sup}}{\displaystyle \frac{|\tilde{B}\left(s\right)|}{q\left(s\right)}}\mathrm{as}n⟶\infty$,

where $\left\{\tilde{B}\left(s\right),0⩽s⩽1\right\}$ is a Brownian Bridge on $[0,1]$.

Proof. 

See Appendix A.    □

It is worth noting that if the change occurs at the very beginning or at the very end of the data, we may not have sufficient observations to obtain consistent LSE estimators of the parameters or these may not be unique. This is why we stress on the truncated version of the test statistic given in [21] that we recall:

$${\Lambda }_n=\underset{\frac{\nu }{n}⩽s⩽1-\frac{\nu }{n}}{max}\frac{|T_n\left(s\right)|}{{\widehat{\sigma }}_w},\mathrm{for}\mathrm{any}1⩽\nu <{\displaystyle \frac{n}{2}}.$$

By Theorem 2, it is easy to see that for any $1⩽\nu <n/2$,

$$\underset{\frac{\nu }{n}⩽s⩽1-\frac{\nu }{n}}{sup}{\displaystyle \frac{|T_n\left(s\right)|}{{\widehat{\sigma }}_w}}-\underset{\frac{\nu }{n}⩽s⩽1-\frac{\nu }{n}}{sup}{\displaystyle \frac{|\tilde{B}\left(s\right)|}{q\left(s\right)}}\stackrel{\mathcal{D}}{⟶}0\mathrm{as}n⟶\infty ,$$

which yields the asymptotic null distribution of the test statistic. With this, at level of significance $\alpha \in (0,1)$, $H_0$ is rejected if ${\Lambda }_n>C_{\alpha ,n}$, where $C_{\alpha ,n}$ is the $(1-\alpha )$-quantile of the distribution of the above limit. This quantile can be computed by observing that under $H_0$, for larger values of n, one has

$$\begin{array}{ccc}\hfill \alpha & =& \mathbb{P}\left(\underset{\frac{\nu }{n}⩽s⩽1-\frac{\nu }{n}}{sup}\frac{|T_n\left(s\right)|}{{\widehat{\sigma }}_w}>C_{\alpha ,n}\right)\hfill \\ & \approx & \mathbb{P}\left(\underset{h_{\nu }\left(n\right)⩽s⩽1-h_{\nu }\left(n\right)}{sup}{\displaystyle \frac{|\tilde{B}\left(s\right)|}{q\left(s\right)}}>C_{\alpha ,n}\right),\mathrm{where}{h}_{\nu }\left(n\right)=\frac{\nu }{n}.\hfill \end{array}$$

From the following relation (1.3.26) of [34], for each $h_{\nu }\left(n\right)>0$, and for larger real number x, we have

$$\begin{array}{ccc}{\displaystyle \mathbb{P}\left\{\underset{h_{\nu }\left(n\right)⩽s⩽1-h_{\nu }\left(n\right)}{sup}{\displaystyle \frac{|\tilde{B}\left(s\right)|}{q\left(s\right)}}⩾x\right\}}\hfill & =\hfill & {\displaystyle {\displaystyle \frac{1}{\sqrt{2\pi }}}xexp\left({\displaystyle \frac{-x^2}{2}}\right)[ln\left({\displaystyle \frac{{\left(1-h_{\nu }\left(n\right)\right)}^2}{h_{\nu }^2\left(n\right)}}\right)}\hfill \\ \hfill & \hfill & {\displaystyle -{\displaystyle \frac{1}{x^2}}ln\left({\displaystyle \frac{{\left(1-h_{\nu }\left(n\right)\right)}^2}{h_{\nu }^2\left(n\right)}}\right)+{\displaystyle \frac{4}{x^2}}+\mathcal{O}\left({\displaystyle \frac{1}{x^4}}\right)],}\hfill \end{array}$$

(16)

which gives an approximation of the tail distribution of ${\displaystyle \underset{h_{\nu }\left(n\right)⩽s⩽1-h_{\nu }\left(n\right)}{sup}}\left|\tilde{B}\left(s\right)\right|/q\left(s\right)$. Thus, using ${\widehat{\sigma }}_w$, an estimation of $C_{\alpha ,n}$ can be obtained from this approximation. Monte Carlo simulations are often carried out to obtain accurate approximations of $C_{\alpha ,n}$. In this purpose, it is necessary to make a good choice of $\nu$. We selected $\nu =0.9\times n^{4/5}$ as our option, which we found to be a suitable choice for all the cases we examined. But, to avoid the difficulties associated with the computation of $C_{\alpha ,n}$, a decision can also be taken by using the p-value method as in [35]. That is using the approximation (16), reject $H_0$ if

$$\mathbb{P}\left({\displaystyle \frac{|\tilde{B}\left(s\right)|}{q\left(s\right)}}>{\Lambda }_n\right)\le \alpha .$$

This idea is used in the simulation section.

3.2.2. Rate of Convergence of the Change Location Estimator

For the study of the estimator ${\widehat{t}}^{\ast }$, we let $\kappa ={\kappa }_n={\vartheta }_2^2-{\vartheta }_1^2$ and assume without loss of generality that ${\kappa }_n>0\left({\vartheta }_2>{\vartheta }_1\right),{\kappa }_n⟶0\mathrm{as}n⟶\infty$ (see, e.g., [36]) and that the unknown change point $t^{\ast }$ depends on the sample size n. We have the following result:

Theorem 3.

Assume that $\left(A_4\right)$ is satisfied, $t^{\ast }/n\in \left(a,1-a\right)$ for some $0<a<1/2$, $t^{\ast }=\left[n\tau \right]$ for some $\tau \in \left(0,1\right)$ and as $n⟶\infty$, ${\kappa }_n⟶0$ and ${\displaystyle \frac{{\kappa }_n\sqrt{n}}{\sqrt{lnn}}}⟶\infty .$ Then, we have

$${\widehat{t}}^{\ast }-t^{\ast }=O_{\mathbb{P}}\left({\displaystyle \frac{1}{{\kappa }_n^2}}\right),$$

where $O_{\mathbb{P}}$ denotes a “big-O” of Landau in probability.

Proof. 

See Appendix A.    □

3.2.3. Limit Distribution of the Location Estimator

In this section, we study the asymptotic behavior of the location estimator. We make the additional assumptions that ${\kappa }_n>>n^{-\frac{1}{2}}$ and that as $n⟶\infty$,

$${\displaystyle \frac{{\kappa }_n\sqrt{n}}{\sqrt{lnn}}}⟶\infty \mathrm{and}{n}^{\frac{1}{2}-\zeta }{\kappa }_n⟶\infty \mathrm{for}\mathrm{some}\zeta \in \left(0,\frac{1}{2}\right).$$

By (10), we have

$$\begin{array}{ccc}\hfill {\widehat{t}}^{\ast }& =& \underset{1⩽k<n}{arg\; max}n\left({\Delta }_k^2-{\Delta }_{t^{\ast }}^2\right).\hfill \end{array}$$

(17)

To derive the limiting distribution of ${\widehat{t}}^{\ast }$, we study the behavior of $n\left({\Delta }_k^2-{\Delta }_{t^{\ast }}^2\right)$ for those ks in the neighborhood of $t^{\ast }$ such that $k=\left[t^{\ast }+r{\kappa }_n^{-2}\right]$, where r varies in an arbitrary bounded interval $\left[-N,N\right]$. For this purpose, we define

$$\begin{array}{c}\hfill P_n\left(r\right):=n\left\{{\Delta }_n^2\left(\left[t^{\ast }+r{\kappa }_n^{-2}\right]\right)-{\Delta }_n^2\left(t^{\ast }\right)\right\},\end{array}$$

where ${\Delta }_n\left(r\right)={\Delta }_{\left[r\right]}$. In addition, we define the two-sided standard Wiener process $\left\{B^{\ast }\left(r\right),r\in \mathbb{R}\right\}$ as follows:

$$B^{\ast }\left(r\right):=\left\{\begin{array}{ccc}{B}_1\left(-r\right)\hfill & \mathrm{if}\hfill & r<0\hfill \\ B_2\left(r\right)\hfill & \mathrm{if}\hfill & r⩾0,\hfill \end{array}\right.$$

where $B_i\left(r\right),i=1,2$ are two independent standard Wiener processes defined on $\left[0,\infty \right)$ with $B_i\left(0\right)=0,i=1,2$.

First, we identify the limit of the process $P_n\left(r\right)$ on $\left|r\right|⩽N$ for every given $N>0$. We denote by $C\left(\left[-N,N\right]\right)$ the space of all continuous functions on $\left[-N,N\right]$ endowed with the uniform metric.

Proposition 1.

Assume that $\left(A_4\right)$ holds, that $t^{\ast }=\left[n\tau \right]$ for some $\tau \in \left(0,1\right)$ and that as $n⟶\infty$, ${\kappa }_n⟶0\mathrm{and}{\displaystyle \frac{{\kappa }_n\sqrt{n}}{\sqrt{lnn}}}⟶\infty$. Then, for every $0<N<\infty$, the process $P_n\left(r\right)$ converges weakly in $C\left(\left[-N,N\right]\right)$ to the process $P\left(r\right)=2\left\{{\sigma }_w{B}^{\ast }\left(r\right)-{\displaystyle \frac{1}{2}}\left|r\right|\right\}$, where $B^{\ast }(·)$ is the two-sided standard Wiener process defined above.

Proof. 

See Appendix A.    □

The above results make it possible to achieve a weak convergence result for $n\left({\Delta }_k^2-{\Delta }_{t^{\ast }}^2\right)$ and then apply the Argmax-Continuous Mapping Theorem (Argmax-CMT). We have:

Theorem 4.

Assume that $\left(A_4\right)$ is satisfied, that $t^{\ast }=\left[n\tau \right]$ for some $\tau \in \left(0,1\right)$ and as $n⟶\infty$, ${\kappa }_n⟶0\mathrm{and}{\displaystyle \frac{{\kappa }_n\sqrt{n}}{\sqrt{lnn}}}⟶\infty$. Then we have

$${\displaystyle \frac{{\kappa }_n^2\left({\widehat{t}}^{\ast }-t^{\ast }\right)}{{\sigma }_w^2}}\stackrel{\mathcal{D}}{⟶}\mathcal{S},$$

where $\mathcal{S}:=arg\; max\left\{B^{\ast }\left(u\right)-{\displaystyle \frac{1}{2}}\left|u\right|,u\in \mathbb{R}\right\}$.

Proof. 

See Appendix A.    □

This result yields the asymptotic distribution of the change location estimator. Refs. [37,38,39] investigated the density function of the random variable $\mathcal{S}$ (see Lemma 1.6.3 of [34] for more details). They also showed that $\mathcal{S}$ has a symmetric (with respect to 0) probability density function $\gamma (·)$ defined for any $x\in \mathbb{R}$ by

$$\gamma \left(x\right)={\displaystyle \frac{3}{2}}exp\left(\right|x\left|\right)\Phi \left({\displaystyle \frac{-3}{2}}\sqrt{\left|x\right|}\right)-{\displaystyle \frac{1}{2}}\Phi \left({\displaystyle \frac{-1}{2}}\sqrt{\left|x\right|}\right),$$

where $\Phi (·)$ is the cumulative distribution function of the standard normal variable. From this result, a confidence interval for the change-point location can be obtained, if one has consistent estimates of ${\kappa }_n^2$ and ${\sigma }_w^2$. With ${\widehat{t}}^{\ast }$, consistent estimates of ${\vartheta }_1^2$ and ${\vartheta }_2^2$ are given, respectively, by

$${\widehat{\vartheta }}_1^2={\overline{W}}_{{\widehat{t}}^{\ast }}={\displaystyle \frac{1}{{\widehat{t}}^{\ast }}}\sum _{t=1}^{{\widehat{t}}^{\ast }}{W}_t^2,\mathrm{and}{\widehat{\vartheta }}_2^2={\overline{W}}_{n-{\widehat{t}}^{\ast }}={\displaystyle \frac{1}{n-{\widehat{t}}^{\ast }}}\sum _{t={\widehat{t}}^{\ast }+1}^n{W}_t^2.$$

Thus, a consistent estimate of ${\kappa }_n^2$ is given by

$${\widehat{\kappa }}_n^2={\displaystyle \frac{1}{n-{\widehat{t}}^{\ast }}}\sum _{t={\widehat{t}}^{\ast }+1}^n{W}_t^2-{\displaystyle \frac{1}{{\widehat{t}}^{\ast }}}\sum _{t=1}^{{\widehat{t}}^{\ast }}{W}_t^2.$$

A consistent estimator of ${\sigma }_w^2$ that we denote by ${\widehat{\sigma }}_w^2$ can be easily obtained by taking its empirical counterpart. So, at risk $\alpha \in (0,1)$, letting $q_{1-\frac{\alpha }{2}}$ be the quantile of order $1-\frac{\alpha }{2}$ of the distribution of the random variable $\mathcal{S}$, an asymptotic confidence interval for $t^{\ast }$ is given by

$$\mathrm{CI}={\widehat{t}}^{\ast }\pm \left(\left[q_{1-\frac{\alpha }{2}}{\displaystyle \frac{{\widehat{\sigma }}_w^2}{{\widehat{\kappa }}_n^2}}\right]+1\right).$$

Remark 1.

In the case that the parameter ρ is unknown, it can be estimated by the CLS method (see Section 2.2.1), and be substituted for its estimator in $W_t$. Indeed, one can easily show that

$${\displaystyle \frac{1}{k}}\sum _{t=1}^k{W}_t^2={\displaystyle \frac{1}{k}}\sum _{t=1}^k{\widehat{W}}_t^2+o_{\mathbb{P}}\left(1\right)\mathrm{and}{\displaystyle \frac{1}{n-k}}\sum _{t=k+1}^n{W}_t^2={\displaystyle \frac{1}{n-k}}\sum _{t=k+1}^n{\widehat{W}}_t^2+o_{\mathbb{P}}\left(1\right),$$

where for any $t=1,\dots ,n$, ${\widehat{W}}_t=\left(X_t-m\left({\widehat{\rho }}_n;Z_{t-1}\right)\right)/{\delta }_0\left(Z_{t-1}\right)$ and ${\widehat{\rho }}_n$ is the conditional least squares estimators of ρ obtained from Theorem 1. Hence, the same techniques as in the case where ρ is known can be used.

4. Practical Consideration

In this section we perform numerical simulations to evaluate the performances of our methods and these are applied to two sets of real data. We start with the presentation of the results of numerical simulations found with the software R 4.1.1. The trials are based on 1000 replications of observations of lengths $n=500$, 1000, 5000 and 10,000 generated from the model (1) for $\rho ={\left({\rho }_0,{\rho }_1,{\rho }_2\right)}^{⊺};$ $\theta ={\left({\theta }_0,{\theta }_1,{\theta }_2\right)}^{⊺};$ $m\left(\rho ;x\right)$ = $\left({\rho }_0+{\rho }_{1}exp\left(-{\rho }_2{x}^2\right)\right)x;$ $\sigma \left(\theta ;x\right)=\underline{\theta }{\delta }_0\left(x\right)$ with ${\delta }_0\left(x\right)$ = $\sqrt{{\theta }_0^2+{\theta }_1^2{x}^{2}exp\left(-{\theta }_2{x}^2\right)};$ ${\rho }_2>0,{\rho }_0{\rho }_1⩾0,{\theta }_2⩾0$ and $0<{\underline{\theta }}^2{\theta }_1^2<1$; ${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$ is a white noise with density function f. We also assume the sufficient condition $\left|{\rho }_0\right|+\left|{\rho }_1\right|+\left|\underline{\theta }{\theta }_1\right|+2{\rho }_0{\rho }_1<1$, to ensure the strict stationarity and ergodicity of the process ${\left(X_t\right)}_{t\in \mathbb{Z}}$ (see, e.g., Theorem 3.2.11 of [40], p. 86 and [41], p. 5). The noise densities f that we employed were Gaussian.

Note that in the application to real data, only the points 2–7 of the above Algorithm 1 are considered, so that the change location estimation is given by point 7.

The change-point location is estimated using the following algorithm:

Algorithm 1 Change-point location estimation

1: for  $i=1,\dots ,1000$ do 2:       for $t=1,\dots ,n$ do $W_t=\left(X_t-m(\rho ;X_{t-1})\right)/{\delta }_0\left(X_{t-1}\right)$ 3:       end for 4:       ${\displaystyle \overline{W}={\displaystyle \frac{1}{n}}\sum _{t=1}^n{W}_t^2}$ 5:       for $k=1,\dots ,n-1$ do ${\displaystyle T_k=\sqrt{{\displaystyle \frac{n}{k(n-k)}}}\sum _{t=1}^k\left(W_t^2-\overline{W}\right)}$ 6:       end for 7:       Compute ${\widehat{t}}_i^{\ast }=\underset{1⩽k<n}{arg\; max}\left|T_k\right|$ (that is the value of k for which $|T_k|$ is the largest) 8: end for 9: Compute $L=\left({\widehat{t}}_1^{\ast }+{\widehat{t}}_2^{\ast }+\dots +{\widehat{t}}_{1000}^{\ast }\right)/1000$ 10: Change-point location estimation is given by ${\widehat{t}}^{\ast }=\left[L\right]$, the integer part of L

4.1. Example 1

We consider the model (1) for ${\rho }_0={\rho }_1=0$, ${\theta }_2=0$, ${\delta }_0\left(X_{t-1}\right)$=$\sqrt{0.04+0.36X_{t-1}^2}$, ${\vartheta }_1=1$, ${\vartheta }_2=1+\varphi$ and ${\epsilon }_t\sim \mathcal{N}(0,1)$. The resulting model is an ARCH(1). The change location estimators are calculated for $\varphi =0.3,0.8\mathrm{and}1.5$ at the locations $t^{\ast }=\tau \times n$ for $\tau =0.25,0.5\mathrm{and}0.75$. In each case, we compute the bias and the standard error SE (SE = SD/$\sqrt{n}$, where SD denotes the standard deviation) of the change location estimator.

Table 1. Change location estimation, its bias and SE for several values of $\varphi$, n and $\tau$ for iid ${\epsilon }_t\sim \mathcal{N}(0,1)$.

ϕ	n	${\mathit{t}}^{\ast }=0.25\mathit{n}(\mathit{\tau }=0.25)$			${\mathit{t}}^{\ast }=0.5\mathit{n}(\mathit{\tau }=0.5)$			${\mathit{t}}^{\ast }=0.75\mathit{n}(\mathit{\tau }=0.75)$
		${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias
0.3	500	181	4.9667	0.1120	277	3.4284	0.0540	384	3.3993	0.0180
	1000	287	3.8961	0.0370	522	2.4946	0.0220	767	2.8024	0.0170
	5000	1264	0.6270	0.0028	2516	0.6260	0.0032	3765	0.8172	0.0030
	10,000	2517	0.4398	0.0017	5015	0.4131	0.0015	7515	0.4701	0.0015
0.80	500	137	1.8286	0.0240	258	0.9659	0.0160	383	1.0514	0.0160
	1000	257	0.5079	0.0070	507	0.6687	0.0070	757	0.6378	0.0070
	5000	1256	0.1874	0.0012	2506	0.1750	0.0012	3755	0.1602	0.0010
	10,000	2506	0.1230	0.0006	5006	0.1388	0.0006	7505	0.1169	0.0005
1.5	500	130	0.8538	0.0100	254	0.4724	0.0080	379	0.4611	0.0080
	1000	253	0.2662	0.0030	504	0.2884	0.0040	753	0.2562	0.0030
	5000	1254	0.1344	0.0008	2503	0.1053	0.0006	3754	0.1174	0.0008
	10,000	2504	0.0880	0.0004	5004	0.0842	0.0004	7504	0.0753	0.0004

shows that the bias declines rapidly as $\varphi$ increases. Also, as the sample size n increases, the bias and the SE decrease. This tends to show the consistency of ${\widehat{t}}^{\ast }$, as expected from the asymptotic results.

We also consider the case ${\epsilon }_t=\beta {\epsilon }_{t-1}+{\gamma }_t$, where $\left|\beta \right|<1\mathrm{and}{\gamma }_t\sim \mathcal{N}(0,\sqrt{1-{\beta }^2})$. It is easy to check that with this ${\left({\epsilon }_t\right)}_{t\in \mathbb{Z}}$ is stationary and strongly mixing, and that $\mathbb{E}\left({\epsilon }_t\right)=0$ and $\mathbb{V}\mathrm{ar}\left({\epsilon }_t\right)=1$. In this case, we only study the SE for $n=5000,$ 10,000 and the results are compared to those obtained for ${\epsilon }_t\sim \mathcal{N}(0,1)$, for the same values of $\varphi$ as above but for $\tau =0.25$ and $0.75$. These results listed in

Table 2. Change location estimation, its bias and SE for several values of $\varphi$, n and $\tau$ for iid ${\epsilon }_t\sim \mathcal{N}(0,1)$ and for ${\epsilon }_t\sim$ AR(1).

ϕ	n	${\mathit{t}}^{\ast }=0.25\mathit{n}(\mathit{\tau }=0.25)$						${\mathit{t}}^{\ast }=0.75\mathit{n}(\mathit{\tau }=0.75)$
		${\mathit{\epsilon }}_{\mathit{t}}\sim \mathcal{N}(\mathbf{0},\mathbf{1})$			${\mathit{\epsilon }}_{\mathit{t}}\sim \mathbf{AR}\left(\mathbf{1}\right)$			${\mathit{\epsilon }}_{\mathit{t}}\sim \mathcal{N}(\mathbf{0},\mathbf{1})$			${\mathit{\epsilon }}_{\mathit{t}}\sim \mathbf{AR}\left(\mathbf{1}\right)$
		${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias
0.3	5000	1265	0.6091	0.0030	1286	2.6225	0.0072	3766	0.6596	0.0032	3776	1.2243	0.0052
	10,000	2516	0.4471	0.0016	2525	0.7136	0.0025	7515	0.4182	0.0015	7523	0.6563	0.0023
0.8	5000	1256	0.1701	0.0012	1260	0.2760	0.0020	3756	0.1818	0.0012	3760	0.2870	0.0020
	10,000	2506	0.1482	0.0006	2510	0.1907	0.0010	7506	0.1338	0.0006	7509	0.1835	0.0009
1.5	5000	1254	0.1165	0.0008	1256	0.1835	0.0012	3754	0.1154	0.0008	3756	0.1784	0.0012
	10,000	2503	0.0807	0.0003	2506	0.1284	0.0006	7504	0.0776	0.0004	7506	0.1263	0.0006

show that for ${\epsilon }_t\sim \mathcal{N}(0,1)$, the location estimator is more accurate and the SE decreases slightly compared to the case ${\epsilon }_t\sim \mathrm{AR}\left(1\right)$. It seems from these results that the nature of the white noise ${\epsilon }_t$ does not much affect the location estimator for larger values of n and $\varphi$.

We present two graphs showing a change in volatility at a time ${\widehat{t}}^{\ast }$. This is indicated by a vertical red line on both graphics where one can easily see the evolution of the time series before and after the change location estimator ${\widehat{t}}^{\ast }$. The series in both figures are obtained for $m(\rho ;x)=0$, ${\delta }_0\left(x\right)=\sqrt{1+0.036x^2}$, $n=500,\tau =0.65,$ and $\varphi =0.8$. That in Figure 1a is obtained for standard iid Gaussian ${\epsilon }_t$s. In this case, using our method, the change location $t^{\ast }=0.65\times 500=325$ is estimated by ${\widehat{t}}^{\ast }=326$. The time series in Figure 1b is obtained with ${\epsilon }_t\sim$ AR(1). In this case, $t^{\ast }$ is estimated by ${\widehat{t}}^{\ast }=325$.

4.2. Example 2

We generate n observations from the model (1) for $m(\rho ;x)=\left(0.5exp\left(-0.03x^2\right)\right)x$ and $\sigma (\theta ;x)=\underline{\theta }{\delta }_0\left(x\right)$, ${\delta }_0\left(x\right)=\sqrt{1+0.02x^2}$, ${\vartheta }_1=1$, ${\vartheta }_2=1+\varphi$ and ${\epsilon }_t\sim \mathcal{N}(0,1)$. We assume $\rho ={({\rho }_0,{\rho }_1,{\rho }_2)}^{⊺}$ is unknown and estimated by CLS method and $\sigma \left(\theta ;.\right)$ is an unknown function that depends on the unknown parameter $\underline{\theta }$. We made 1000 replications for the lengths $n=500,1000,5000\mathrm{and}$ 10,000 from this model. The change location estimator, its bias and SE are calculated for the same values of $\varphi$ and locations $t^{\ast }$ as in the preceding example. The results given in

Table 3. Change location estimation, its bias and SE for several values of $\varphi$, n and $\tau$ for iid ${\epsilon }_t\sim \mathcal{N}(0,1)$.

ϕ	n	${\mathit{t}}^{\ast }=0.25\mathit{n}(\mathit{\tau }=0.25)$			${\mathit{t}}^{\ast }=0.5\mathit{n}(\mathit{\tau }=0.5)$			${\mathit{t}}^{\ast }=0.75\mathit{n}(\mathit{\tau }=0.75)$
		${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias	${\widehat{\mathit{t}}}^{\ast }$	SE	Bias
0.3	500	182	4.9959	0.1140	280	3.2396	0.0600	387	3.0412	0.0180
	1000	298	4.3560	0.0480	525	2.6278	0.0250	770	2.3343	0.0200
	5000	1267	1.7941	0.0034	2517	0.7852	0.0034	3767	0.7948	0.0034
	10,000	2517	0.4716	0.0017	5016	0.4454	0.0016	7513	0.4245	0.0013
0.8	500	139	2.0945	0.0280	259	1.0386	0.0180	384	0.9061	0.0180
	1000	259	1.1755	0.00900	506	0.4205	0.0060	757	0.5427	0.0070
	5000	1256	0.1780	0.0012	2506	0.1713	0.0012	3757	0.2107	0.0014
	10,000	2506	0.1304	0.0006	5007	0.1375	0.0007	7506	0.1236	0.0006
1.5	500	135	1.8053	0.0200	256	0.9248	0.0120	382	0.7279	0.0140
	1000	255	0.3217	0.0050	505	0.3138	0.0050	755	0.4666	0.0050
	5000	1254	0.1469	0.0008	2505	0.1378	0.0010	3754	0.1325	0.0008
	10,000	2505	0.1010	0.0005	5004	0.0912	0.0004	7504	0.0915	0.0004

are very similar to those displayed in Table 1.

As in the previous example, we present two graphs illustrating our method’s ability to detect the change-point in the time series considered. On both graphics, one can easily see the evolution of the time series before and after the change location estimator. The series in both figures are obtained for $m(\rho ;x)=0.5exp(-0.03x^2)x$, ${\delta }_0\left(x\right)=\sqrt{1+0.02x^2}$, $n=500,\tau =0.65,$ and $\varphi =1.5$. That in Figure 1c is obtained for standard i.i.d. Gaussian ${\epsilon }_t$s. In this case, using our method, the change location $t^{\ast }=0.65\times 500=325$ is estimated by ${\widehat{t}}^{\ast }=326$. The time series in Figure 1d is obtained with AR(1) ${\epsilon }_t$’s. In this case, $t^{\ast }$ is estimated by ${\widehat{t}}^{\ast }=326$.

When performing our test, we used the p-value method. In other words, for the nominal level $\alpha =5\%$, we simulated 1000 samples each of length $n=100,200,500$ and 1000 from the model (1) for $m(\rho ;x)=0$ and $\sigma (\theta ;x)=\underline{\theta }{\delta }_0\left(x\right)$, ${\delta }_0\left(x\right)=\sqrt{0.99+0.2x^2}$, ${\vartheta }_1=1$, ${\vartheta }_2=1+\varphi$ and ${\epsilon }_t\sim \mathcal{N}(0,1)$. We then calculated ${\Lambda }_n$ and counted the number of samples for which

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\sqrt{2\pi }}}{\Lambda }_{n}exp\left({\displaystyle \frac{-{\Lambda }_n^2}{2}}\right)\left[ln\left({\displaystyle \frac{{\left(1-h_{\nu }\left(n\right)\right)}^2}{h_{\nu }^2\left(n\right)}}\right)-{\displaystyle \frac{1}{{\Lambda }_n^2}}ln\left({\displaystyle \frac{{\left(1-h_{\nu }\left(n\right)\right)}^2}{h_{\nu }^2\left(n\right)}}\right)+{\displaystyle \frac{4}{{\Lambda }_n^2}}\right]\le \alpha ,\end{array}$$

and we divided this number by 1000. This ratio corresponds to the empirical power of our statistical test for change in volatility. The results obtained are listed in

Table 4. Statistical test powers for different $\varphi$ values at different locations $t^{\ast }=\left[n\tau \right]$.

$\mathit{\tau }$		$\mathit{n}=100$			$\mathit{n}=200$			$\mathit{n}=500$			$\mathit{n}=1000$
		0.25	0.50	0.75	0.25	0.50	0.75	0.25	0.50	0.75	0.25	0.50	0.75
	0	0.051			0.048			0.05			0.05
$\varphi$	0.03	0.145	0.147	0.131	0.121	0.115	0.109	0.093	0.089	0.083	0.071	0.077	0.056
	0.05	0.150	0.151	0.150	0.124	0.120	0.115	0.098	0.090	0.085	0.085	0.084	0.060
	0.08	0.157	0.177	0.158	0.147	0.153	0.131	0.118	0.100	0.096	0.090	0.090	0.070
	0.1	0.191	0.198	0.177	0.170	0.174	0.136	0.145	0.167	0.101	0.100	0.110	0.098
	0.3	0.249	0.296	0.214	0.271	0.358	0.248	0.458	0.530	0.371	0.685	0.750	0.610
	0.5	0.344	0.465	0.315	0.421	0.609	0.433	0.765	0.891	0.780	0.974	0.998	0.992
	0.7	0.413	0.561	0.422	0.591	0.803	0.616	0.932	0.978	0.971	0.995	0.998	0.998
	0.9	0.477	0.710	0.532	0.708	0.887	0.787	0.971	0.996	0.993	0.998	0.999	0.999
	1.1	0.577	0.806	0.654	0.808	0.946	0.897	0.985	0.998	0.999	0.999	1.000	1.000
	1.3	0.634	0.838	0.721	0.863	0.964	0.952	0.990	0.999	0.999	1.000	1.000	1.000
	1.5	0.640	0.860	0.800	0.907	0.967	0.967	0.997	1.000	1.000	1.000	1.000	1.000

. We can clearly see that, when $\varphi =0$, the empirical power of the test is almost the same as the nominal level $\alpha$ for all n sizes (see Table 4).

4.3. Comparison with Some Recent Algorithms

We compare our method, referred to as LS, with the Wild Binary Segmentation (WBS) method studied in [42], and one of its variants, called Narrowest-Over-Threshold (NOT), proposed by [43], as well as, the Iterative Cumulative Sum of Squares (ICSS) algorithm suggested by [24]. All these methods are implemented under R software, and can, respectively, be found in the packages wbs, not and ICSS.

Our comparison is based on n observations simulated from (1) for ${\rho }_0={\rho }_1=0$, ${\theta }_2=0$, ${\delta }_0\left(X_{t-1}\right)=\sqrt{0.04+0.36X_{t-1}^2}$, ${\vartheta }_1=1$, ${\vartheta }_2=1+\varphi$ and ${\epsilon }_t\sim \mathcal{N}(0,1)$. The change location estimators are calculated for $\varphi =0.3,0.8\mathrm{and}1.5$ at the locations $t^{\ast }=\tau \times n$ for $\tau =0.25\mathrm{and}0.75$.

From the results obtained (see

Table 5. Estimates of change location derived from LS, WBS, NOT and ICSS for a sample with a single break.

		$\mathit{\varphi }=0.3$		$\mathit{\varphi }=0.8$		$\mathit{\varphi }=1.5$
		${\mathit{t}}^{\ast }=\mathbf{0}.\mathbf{25}\mathit{n}$	${\mathit{t}}^{\ast }=\mathbf{0}.\mathbf{75}\mathit{n}$	${\mathit{t}}^{\ast }=\mathbf{0}.\mathbf{25}\mathit{n}$	${\mathit{t}}^{\ast }=\mathbf{0}.\mathbf{75}\mathit{n}$	${\mathit{t}}^{\ast }=\mathbf{0}.\mathbf{25}\mathit{n}$	${\mathit{t}}^{\ast }=\mathbf{0}.\mathbf{75}\mathit{n}$
Methods	$\mathit{n}$	${\widehat{\mathit{t}}}^{\ast }$	${\widehat{\mathit{t}}}^{\ast }$	${\widehat{\mathit{t}}}^{\ast }$	${\widehat{\mathit{t}}}^{\ast }$	${\widehat{\mathit{t}}}^{\ast }$	${\widehat{\mathit{t}}}^{\ast }$
LS	100	39	70	31	79	28	78
WBS		46∣43	87∣90	38∣39	82∣80	29∣30	78∣76
NOT		51∣55	96	38∣42	79∣82∣91∣94	38∣35	78∣82
ICSS		42∣48	68	31∣43	78∣81	33∣41	78∣97
LS	200	57	145	54	157	52	155
WBS		70∣66	175∣177	62∣59	171∣167	58∣57	159∣158
NOT		66∣70∣77	175∣190∣194	67∣70	158∣161	64∣68	155∣162
ICSS		66∣87	80∣159	61∣76∣87	159	59∣66∣170∣187	156∣173

), WBS, NOT and ICSS generate sequences of change-point location estimates, some of which have values close to the true locations. For $n=100$ and $n=200$, LS generally provides more accurate estimates ${\widehat{t}}^{\ast }$ of the true change-point location than WBS, NOT and ICSS, for different $\varphi$ and locations $t^{\ast }$. Among all, it is generally the best, especially for larger values of n and $\varphi$.

4.4. Application to Real Data

In this section, we apply our procedure to two sets of genuine time series, namely, the USA stock market prices and the Brent crude oil prices. As these have large values, we take their logarithms and differentiate the resulting series to remove their trends.

4.4.1. USA Stock Market

These data of length 2022 from the American stock market were recorded daily from 2 January 1992 to 31 December 1999. They represent the daily stock prices of the S&P 500 stock market (SPX). They are among the most closely followed stock market indices in the world and are considered as an indicator of the USA economy. They have also been recently examined by [44] and can be found at www.investing.com (accessed on 1 August 2023).

In Figure 2, we observe that the trend of the SPX daily stock price series is not constant over time. We also observe that stock prices have fallen sharply, especially in the time interval between the two vertical dashed blue lines (the period of the 1997–1998 Asian financial crisis).

Denote by $D_t$ the value of the stock price for the SPX index at day t, and the first difference of the logarithm of stock price, $X_t$ as

$$X_t=log\left(D_t\right)-log\left(D_{t-1}\right)=log\left({\displaystyle \frac{D_t}{D_{t-1}}}\right).$$

$X_t$ is the logarithmic return of stock price for the SPX index at day t.

The series $\left(X_t\right)$ is approximately piece-wise stationary on two segments and symmetric around zero (see Figure 3). This brought us to consider a CHARN model with $m(\rho ;x)=0$, ${\delta }_0\left(x\right)=\sqrt{{\theta }_0+{\theta }_1{x}^2}$ for ${\theta }_0=1,{\theta }_1=0$, ${\vartheta }_1$ and ${\vartheta }_2$ estimated by CLS described in Section 2.

Using our procedure, we found an important change point in stock price volatility on 26 March 1997, which is consistent with the date found by [44] (see Figure 3).

The vertical dashed blue line of Figure 3 represents the date at which the change occurred. It should be noted that the change in volatility coincides with the Asian crisis in 1997 when Thailand devalued its currency, the baht, against the US dollar. This decision led to a fall in the currencies and financial markets of several countries in its surroundings. The crisis then spread to other emerging countries with important social and political consequences and repercussions on the world economy.

4.4.2. Brent Crude Oil

These data of length 585 are based on Brent oil futures. They represent the prices of Brent oil (USD/barrel) on a daily basis between 4 January 2021 and 6 April 2023. They are available at www.investing.com (accessed on 1 August 2023).

Figure 4 shows that the evolution of the daily series of Brent oil prices is non-stationary. It also shows that stock prices fell sharply, especially in early March 2022 (the date of the conflict between OPEC and Russia, when the latter refused to reduce its oil production in the face of declining global demand caused by the COVID-19 pandemic).

We follow the same procedure as in the previous example and obtain the logarithmic transformation of the daily rate of return series for Brent oil (see Figure 5 below). Proceeding as for the first data, the application of our procedure allows to find a change at 25 February 2022. The break date is marked by a dashed blue vertical line (see Figure 5).

It also performs well on these data. Indeed, oil volatility was very high in March 2022 due to the COVID-19 pandemic and the conflict between OPEC and Russia. The health crisis led to a significant drop in global oil demand, while Russia refused to cut oil production as proposed by OPEC, which caused oil prices to fall. Brent crude oil fell from over 60 dollars a barrel to less than 20 dollars in a month.

4.4.3. Comparison with WBS and ICSS

This time, we compare our LS method with both the WBS and ICSS algorithms for real S&P 500 stock market and Brent crude oil data, just as we carried out for simulated data. The results shown in

Table 6. Estimates of change location derived from LS, WBS, and ICSS for real data S&P 500 stock prices and Brent crude oil prices.

	S&P 500 Stock Prices	Brent Crude Oil Prices
Methods	${\widehat{\mathit{t}}}^{\ast }$	${\widehat{\mathit{t}}}^{\ast }$
LS	26 March 1997	25 February 2022
WBS	12 May 1997	26 January 2022
ICSS	22 December 1998 ∣ 26 January 1998	14 January 2022
	20 April 1995 ∣ 14 November 1996

, indicate that WBS and ICSS produce sequences of change-point location estimates, some of which are close to the actual locations. In contrast, LS provides more accurate ${\widehat{t}}^{\ast }$ estimates of the actual change-point location than WBS and ICSS. This clearly demonstrates that our method remains robust for both real and simulated data. Of all these methods, LS generally stands out as the best for the examples we have considered.

5. Conclusions

In this article, we have presented a CUSUM test based on a least-squares statistics for detecting an abrupt change in the volatility of CHARN models. We have shown that the test statistic converges to a supremum of a weighted standard Brownian Bridge. We have constructed a change location estimator and obtained its rate of convergence and its limiting distribution whose density function is given. A simulation experiment shows that on the example tested, our method generally performs better than the competitors considered. Applied to real equity price data and real Brent crude oil data, our procedure finds the change locations found in the literature. This is not the case with WBS and ICSS. The next step to this work is its extension to the detection of multiple changes in volatility of multivariate CHARN models.