Robust Procedure for Change-Point Estimation Using Quantile Regression Model with Asymmetric Laplace Distribution

Abstract

The usual mean change-point detecting method based on normal linear regression is not robust to heavy-tailed data with potential outlying points. We propose a robust change-point estimation procedure based on a quantile regression model with asymmetric Laplace error distribution and develop a non-iterative sampling algorithm from a Bayesian perspective. The algorithm can generate independently and identically distributed samples approximately from the posterior distribution of the position of the change-point, which can be used for statistical inferences straightforwardly. The procedure combines the robustness of quantile regression and the computational efficiency of the non-iterative sampling algorithm. A simulation study is conducted to illustrate the performance of the procedure with satisfying findings, and finally, real data is analyzed to show the usefulness of the algorithm by comparison with the usual change-point detection method based on normal regression.

1. Introduction

Regression models incorporating change-points have been popular in many application areas, such as finance (Chen [1]), medicine research (Kim et al. [2]), biology (Chen and Wang [3], Acitas and Senoglu [4]), climate (Reeves et al. [5]). In these documents, change-point detection in the normal linear regression model is probably the most popular. Chen and Gupta [6] first used the Schwarz Information Criterion (SIC) model selection method to detect the mean and variance change-points in the normal linear regression model. Chen and Wang [3] applied the normal mean change-point model to detect DNA copy number variations using the SIC method. Chen et al. [7] compared four methods for estimating single or multiple change-points in a normal regression model when both the error variance and regression coefficients change simultaneously at the unknown point. Meanwhile, many robust change-point detecting models with heavy-tailed distributions have been developed to deal with data with outliers. Osorio and Galea [8] discussed the regression coefficient change-point model with the Student-t error assumption and designed the SIC method to detect the position of the change-point. Moreover, in the Student t regression framework, Lin et al. [9] proposed a kind of Bayesian methodology to detect the variance change-point and analyze the U.S. stock market data. Kang et al. [10] extended Lin’ work to linear models with variance change-point by the assumption of a class of symmetric heavy-tailed errors which includes the Student-t distribution. Lu and Chang [11] developed a robust algorithm to detect regression coefficient change-point by incorporating the Student-t distribution with the expectation-maximization (EM) algorithm and the fuzzy classification procedure. Yang [12] designed a robust mean change-point detection procedure by employing the Laplace regression model via the EM algorithm. Alin et al. [13] proposed a robust F-test for change-point detection in linear regression models with a non-normal error structure.

The quantile regression (QR) model, pioneered by Koenker and Bassett [14], has been powerful in comprehensively assessing the effects of explanatory covariates at different quantiles of response outcomes. In contrast to the usual mean regression model, the QR model is less sensitive to outliers in the data; thus, when the data comes from some heavy-tailed distribution, the inferences based on QR are more credible than the ones derived from the usual mean regression. Since a set of quantiles often provides a more detailed description of the response distribution than the usual mean, QR offers a practically important alternative to the classical mean regression. Koenker [15] provides a detailed introduction to QR in theories, applications and computation. In recent years, many authors have devoted themselves to studying the change-point problem in the QR framework; Qu [16] proposed several tests for structural change in regression quantiles, such as fluctuation type statistic and Wald type statistic. Zhang et al. [17] developed a procedure for testing change-point due to covariate threshold in regression quantiles. Gabriela [18] considered the estimation problem in a nonlinear quantile model with change-points. Guo and He [19] discussed the quantile auto-regression change-point issue from a Bayesian perspective and developed the monte carlo markov chain (MCMC) algorithm. Although MCMC algorithms are widely used in Bayesian inference for complex data, we should mention that it is an iterative sampling algorithm; thus, the samples generated from the full conditional distributions are not independent, and the convergence assessment is first needed before inferences using these samples. Despite several well-known convergence criteria, determining when the markov chain achieves convergence is still a challenging task in Bayesian analysis.

Tan et al. [20] discovered an equation called the Inverse Bayes formula (IBF) and developed it into a sampling algorithm named IBF sampler, which is highly useful in doing bayesian sampling. Compared to the classical MCMC sampling, this sampler possesses two main advantages. First, it is non-iterative and thus fast. Second, it can draw independently and identically distributed (i.i.d) samples exactly or approximately from the target posterior distribution, and the samples can be used for statistical inferences immediately without convergence diagnosis. The main idea of the IBF sampler is to use the EM algorithm to obtain an optimal importance sampling density and then use the sampling/importance resampling algorithm to get i.i.d. samples from the target distribution. To get a comprehensive understanding of the IBF sampler, one can refer to Tan et al. [21], which gave a detailed introduction to IBF, including its extensions and wide applications. In recent years, the applications of the IBF algorithm appear prosperous, to name a few, Yang and Yuan [22,23] applied the idea of IBF to deal with quantile regression model and scale mixture normal regression model, respectively. Tsionas [24] used the IBF sampler to research financial areas with stochastic volatility models. As to change-point detection, Tian et al. [25] was the first to apply the non-iterative Bayesian method to deal with change-point problem, and the authors identified the position of change-points in a sequence of hemolytic uremic syndrome using IBF sampler. Recently, Yang and Yuan [26,27] used the idea of an IBF sampler to detect the position of change-point in normal mean and normal regression coefficient, respectively.

Motivated by Tan et al. [20,21], in this article, a non-iterative Bayesian sampling algorithm to detect the position of regression coefficient change-point based on the QR model is discussed. The procedure involves several key steps. Firstly, the representation of asymmetric Laplace distribution (ALD) makes it possible to transfer the observed likelihood into a missing data structure by augmenting a latent exponential distributed variable, which is crucial for implementing the EM algorithm and IBF sampler. Then, the best posterior importance sampling distribution to overlap the target density is obtained by using the posterior mode from the EM algorithm. Finally, the sampling IBF and sampling/importance resampling algorithm are combined to get i.i.d. samples approximately from the observed posterior distributions. Simulation studies are conducted to evaluate the robustness of the procedure, and real data are analyzed to illustrate the usefulness of the algorithm by comparison with the classical method to detect the position of the change-point in regression models. The novelty of this paper lies in that it is the first time to apply the IBF algorithm to focus on a change-point problem based on a quantile regression model, which combines the robustness of quantile regression in detecting change-point and the computational effectiveness of the IBF sampler.

2. Quantile Regression (QR) Model

Let Y denote a response variable and X the related explanatory vector with p components, and suppose the observations are $(y_i,{\mathbf{x}}_i),i=1,\cdots ,n$. We introduce the linear quantile regression model as

$$\begin{array}{c}\hfill y_i={\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_q+{\epsilon }_i,(i=1,\cdots ,n),\end{array}$$

(1)

where ${\mathbf{\beta }}_q$ is a unknown parameters vector of interest, ${\epsilon }_i$ is the error term whose distribution is restricted to have the q-th quantile equal to zero, and ${\mathbf{x}}_i^{\prime }$ denotes the transpose of ${\mathbf{x}}_i$. So, ${\mathbf{\beta }}_q$ quantified the q-th conditional quantile regression coefficient of $y_i$ given ${\mathbf{x}}_i$, for $i=1,\cdots ,n$, and hereafter, ${\mathbf{\beta }}_q$ is simplified as $\mathbf{\beta }$. The estimation $\widehat{\mathbf{\beta }}$ of $\mathbf{\beta }$ is the one minimizing

$$S(\mathbf{\beta })=\sum _{i=1}^n{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }\mathbf{\beta }),$$

with ${\rho }_q(·)$ denoting the loss function ${\rho }_q(u)=u\{q-I(u<0)\},$ and here $I(·)$ is the indicator function. Note that when $q=0.5$, the loss function is the symmetric absolute loss, which leads to the Laplace regression or median regression. Realizing that ${\rho }_q(u)$ is not differentiable at zero, $\widehat{\mathbf{\beta }}$ cannot be obtained in a closed form. In documents, some elaborate methods, such as the interior point algorithm and the Majorization–Minimization algorithm are developed to obtain $\widehat{\mathbf{\beta }}$, see Koenker and Park [28], Hunter and Lange [29]. Many authors have realized that minimizing $S(\mathbf{\beta })$ is equivalent to maximizing the likelihood function of a linear regression model with errors following the asymmetric Laplace distribution (ALD$(0,\sigma ,q)$), whose density function is given by

$$\begin{array}{c}\hfill f(u)=\frac{q(1-q)}{\sigma }\exp \{-{\rho }_q(\frac{u}{\sigma })\},u\in \mathbb{R}.\end{array}$$

(2)

In the literature, Koenker and Machado [30] studied the likelihood ratio test statistics quantile regression coefficients with ALD assumption, Yu and Moyeed [31] proposed the statistical inference of QR from a Bayesian perspective based on the ALD error distribution. Using a novel location-scale mixture representation of ALD, Reed and Yu [32], Kozumi and Kobayashi [33] developed an efficient and easily implemented Gibbs sampling method for quantile regression in the Bayesian framework. This representation has been widely used in the latest parameter quantile regression studies, such as longitudinal data and variable selection. To name a few, Zhou et al. [34] and Tian et al. [35] proposed the EM algorithm for computing the MLE or posterior mode of regression coefficients in the linear quantile regression model, Luo et al. [36] studied the quantile regression model for longitudinal data using a Gibbs sampler. Alhamzawi and Yu [37], Ji et al. [38] studied the variable selection problems in quantile regression, tobit and binary quantile regression. Tian et al. [39] discussed the inferences problems of a mixture of linear quantile regression using the EM algorithm. We notice that recently, Gallardo et al. [40] introduced a new parametric quantile regression model for asymmetric response variables with a power skew–normal distribution, and Reyes et al. [41] proposed a kind of quantile regression model with a generalization of the Student-t distribution as response variables, although these two distributions are different from the asymmetric Laplace distribution in this paper, the ideas of re-parametrization of these two distributions are similar to the mixture representation of ALD.

3. Quantile Regression Model with Change-Point

In this paper, we suppose a change-point exists in the q-th quantile regression coefficients at an unknown position k, such that

$$\begin{array}{c}\hfill \begin{array}{cc}& y_i={\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1+{\epsilon }_i,i=1,2,\cdots ,k,\hfill \\ & y_i={\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2+{\epsilon }_i,i=k+1,\cdots ,n,\hfill \\ & {\epsilon }_i\sim \mathrm{ALD}(0,1,q),i=1,2,\cdots ,n.\hfill \end{array}\end{array}$$

(3)

where ${\mathbf{\beta }}_1={({\beta }_1,\cdots ,{\beta }_p)}^{\prime }$ and ${\mathbf{\beta }}_2={({\beta }_1^{*},\cdots ,{\beta }_p^{*})}^{\prime }$ are the q-th conditional quantile regression coefficients for the first k observations and the last $n-k$ observations, respectively. For simplicity, the scale parameter $\sigma$ in the ALD error distribution is set to 1. Our principal interest involves the determination of the unknown change-point position k and other parameters from a Bayesian perspective. Denote $\mathbf{\theta }=({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k)$ and observation vector $\mathbf{y}={(y_1,\cdots ,y_n)}^{\prime }$, then the likelihood function of $\mathbf{\theta }$ based on $\mathbf{y}$ is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L(\mathbf{\theta }|\mathbf{y})& =\prod _{i=1}^{k}q(1-q)\exp \{-{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1)\}\hfill \\ & \times \prod _{i=k+1}^{n}q(1-q)\exp \{-{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2)\}\hfill \\ & \propto \exp \left\{-[\sum _{i=1}^k{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1)+\sum _{i=k+1}^n{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2)]\right\}.\hfill \end{array}\end{array}$$

(4)

By combining this likelihood with some prior distributions, we can get the kernel of the posterior distribution, but the complexity of the ALD makes the posterior inference difficult. So Reed and Yu [32], Kozumi and Kobayashi [33], and Zhou et al. [34] represented the ALD$(0,\sigma ,q)$ into a normal-exponential mixing distribution as Proposition 1.

Proposition 1

(Mixture Representation of the Asymmetric Laplace Distribution). Let z be an exponential variable with mean σ and u a standard normal variable. If a random variable ϵ follows the asymmetric Laplace distribution with density (2), then we can represent ε as a location-scale mixture of normals given by $\epsilon =\theta z+\tau \sqrt{z}u,$ where

$$\theta =\frac{1-2q}{q(1-q)},{\tau }^2=\frac{2}{q(1-q)}.$$

Based on this representation, the quantile regression model with a single change-point (3) becomes

$$\begin{array}{c}\hfill \begin{array}{cc}& y_i|z_i\sim N({\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1+\theta z_i,{\tau }^2{z}_i),i=1,\cdots ,k.\hfill \\ & y_i|z_i\sim N({\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2+\theta z_i,{\tau }^2{z}_i),i=k+1,\cdots ,n.\hfill \\ & z_i\sim \mathrm{Exp}(1),i=1,2,\cdots ,n.\hfill \end{array}\end{array}$$

(5)

where $z_i,i=1,\cdots ,n$ are mutually independent. By introducing latent vector $\mathbf{z}={(z_1,\cdots ,z_n)}^{\prime }$, the likelihood function based on complete data $(\mathbf{y},\mathbf{z})$ appears as

$$\begin{array}{c}\hfill \begin{array}{cc}& L({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y},\mathbf{z})\hfill \\ & =\prod _{i=1}^k\frac{1}{\sqrt{2\pi }\sqrt{{\tau }^2{z}_i}}\exp \{-\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1-\theta z_i)}^2}{2{\tau }^2{z}_i}-z_i\}\hfill \\ & \times \prod _{i=k+1}^n\frac{1}{\sqrt{2\pi }\sqrt{{\tau }^2{z}_i}}\exp \{-\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2-\theta z_i)}^2}{2{\tau }^2{z}_i}-z_i\}.\hfill \end{array}\end{array}$$

(6)

4. Bayesian Inferences and IBF Algorithm

4.1. The Prior and Conditional Distributions

To emphasize the influence of data on the position of change-point and other parameters, we adopt a set of independent non-informative priors as follows

$$\pi ({\mathbf{\beta }}_1)\propto 1,\pi ({\mathbf{\beta }}_2)\propto 1,\pi (k)\propto \frac{1}{n-2p+1},k=p,\cdots ,n-p,$$

here the prior of k is a discrete uniform distribution on $\{p,p+1,\cdots ,n-p\}$. The implementation of the IBF algorithm involves several conditional distributions, which are $\pi ({\mathbf{\beta }}_1|\mathbf{y},\mathbf{z},k)$, $\pi ({\mathbf{\beta }}_2|\mathbf{y},\mathbf{z},k)$, $\pi (k|\mathbf{y},\mathbf{z})$, and $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k)$. These distributions are all proportional to the joint conditional distribution of the parameters, that is

$$\begin{array}{c}\hfill \begin{array}{cc}& \pi ({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y},\mathbf{z})=\pi ({\mathbf{\beta }}_1|\mathbf{y},\mathbf{z},k)\pi ({\mathbf{\beta }}_2|\mathbf{y},\mathbf{z},k)\pi (k|\mathbf{y},\mathbf{z})\hfill \\ & \propto L({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y},\mathbf{z})\pi ({\mathbf{\beta }}_1)\pi ({\mathbf{\beta }}_2)\pi (k)\hfill \\ & \propto \exp \{-\frac{1}{2}\sum _{i=1}^k\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1-\theta z_i)}^2}{{\tau }^2{z}_i}\}\hfill \\ & \times \exp \{-\frac{1}{2}\sum _{i=k+1}^n\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2-\theta z_i)}^2}{{\tau }^2{z}_i}\}\hfill \\ & \times \exp \{-\sum _{i=1}^n{z}_i\}\mathbf{I}(p\le k\le n-p),\hfill \end{array}\end{array}$$

(7)

with $\mathbf{I}(·)$ representing the indicator function. Denote

$${\mathbf{X}}_1={({\mathbf{x}}_1,\cdots ,{\mathbf{x}}_k)}^{\prime },{\mathbf{W}}_1=\mathrm{diag}({\tau }^2{z}_1,\cdots ,{\tau }^2{z}_k),$$

$${\mathbf{z}}_1={(z_1,\cdots ,z_k)}^{\prime },{\mathbf{y}}_1={(y_1,\cdots ,y_k)}^{\prime }$$

and

$${\mathbf{b}}_1={({\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{-1}{\mathbf{X}}_1)}^{-1}{\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{-1}({\mathbf{y}}_1-\theta {\mathbf{z}}_1),$$

we have

$$\begin{array}{c}\hfill \pi ({\mathbf{\beta }}_1|\mathbf{y},\mathbf{z},k)\propto \exp \left\{-\frac{{({\mathbf{\beta }}_1-{\mathbf{b}}_1)}^{\prime }{\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{-1}{\mathbf{X}}_1({\mathbf{\beta }}_1-{\mathbf{b}}_1)}{2}\right\},\end{array}$$

(8)

consequently,

$$\begin{array}{c}\hfill {\mathbf{\beta }}_1|\mathbf{y},\mathbf{z},k\sim \mathbb{N}({\mathbf{b}}_1,{({\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{-1}{\mathbf{X}}_1)}^{-1}).\end{array}$$

(9)

Likewise, let

$${\mathbf{X}}_2={({\mathbf{x}}_{k+1},\cdots ,{\mathbf{x}}_n)}^{\prime },{\mathbf{W}}_2=\mathrm{diag}({\tau }^2{z}_{k+1},\cdots ,{\tau }^2{z}_n),$$

$${\mathbf{z}}_2={(z_{k+1},\cdots ,z_n)}^{\prime },{\mathbf{y}}_2={(y_{k+1},\cdots ,y_n)}^{\prime },$$

and

$${\mathbf{b}}_2={({\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{-1}{\mathbf{X}}_2)}^{-1}{\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{-1}({\mathbf{y}}_2-\theta {\mathbf{z}}_2),$$

we have

$$\begin{array}{c}\hfill \pi ({\mathbf{\beta }}_2|\mathbf{y},\mathbf{z},k)\propto \exp \left\{-\frac{{({\mathbf{\beta }}_2-{\mathbf{b}}_2)}^{\prime }{\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{-1}{\mathbf{X}}_2({\mathbf{\beta }}_2-{\mathbf{b}}_2)}{2}\right\},\end{array}$$

(10)

consequently,

$$\begin{array}{c}\hfill {\mathbf{\beta }}_2|\mathbf{y},\mathbf{z},k\sim \mathbb{N}({\mathbf{b}}_2,{({\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{-1}{\mathbf{X}}_2)}^{-1}).\end{array}$$

(11)

Next, by integrating ${\mathbf{\beta }}_1,{\mathbf{\beta }}_2$ from the joint distribution $\pi ({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y},\mathbf{z})$ in Equation (4), one can get the conditional distribution of the position of change-point k based on complete data $(\mathbf{y},\mathbf{z})$, that is for $k=p,p+1,\cdots ,n-p$,

$$\begin{array}{c}\hfill \pi (k|\mathbf{y},\mathbf{z})\propto \frac{|{\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{-1}{\mathbf{X}}_1{|}^{-1/2}{|{\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{-1}{\mathbf{X}}_2|}^{-1/2}}{exp(\frac{S_1+S_2+2{\sum }_{i=1}^n{z}_i}{2})},\end{array}$$

(12)

with

$$S_i={({\mathbf{y}}_i-\theta {\mathbf{z}}_i-{\mathbf{X}}_i{\mathbf{b}}_i)}^{\prime }{\mathbf{W}}_i^{-1}({\mathbf{y}}_i-\theta {\mathbf{z}}_i-{\mathbf{X}}_i{\mathbf{b}}_i),i=1,2.$$

Finally, we present the conditional distribution of latent vector $\mathbf{z}$, that is

$$\begin{array}{c}\hfill \begin{array}{cc}& \pi (\mathbf{z}|\mathbf{y},{\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k)\propto L({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y},\mathbf{z})\hfill \\ & \propto \prod _{i=1}^k{z}_i^{-\frac{1}{2}}\exp \{-\frac{1}{2}[\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1)}^2}{{\tau }^2}\frac{1}{z_i}+\frac{2{\tau }^2+{\theta }^2}{{\tau }^2}{z}_i]\}\hfill \\ & \times \prod _{i=k+1}^n{z}_i^{-\frac{1}{2}}\exp \{-\frac{1}{2}[\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2)}^2}{{\tau }^2}\frac{1}{z_i}+\frac{2{\tau }^2+{\theta }^2}{{\tau }^2}{z}_i]\},\hfill \end{array}\end{array}$$

(13)

explicitly, $z_1,\cdots ,z_n$ are mutually conditional independent, and for $i=1,2,\cdots ,k$,

$$\begin{array}{c}\hfill z_i|y_i,{\mathbf{\beta }}_1\sim \mathrm{GIG}(\frac{1}{2},\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1)}^2}{{\tau }^2},\frac{2{\tau }^2+{\theta }^2}{{\tau }^2}),\end{array}$$

(14)

for $i=k+1,\cdots ,n$,

$$\begin{array}{c}\hfill z_i|y_i,{\mathbf{\beta }}_2\sim \mathrm{GIG}(\frac{1}{2},\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2)}^2}{{\tau }^2},\frac{2{\tau }^2+{\theta }^2}{{\tau }^2}),\end{array}$$

(15)

where $\mathrm{GIG}(\lambda ,a,b)$ denotes the generalized inverse Gaussian distribution with density

$$f(x|\lambda ,a,b)=\frac{{(b/a)}^{\lambda /2}}{2K_{\lambda }(\sqrt{ab})}{x}^{\lambda -1}\exp \{-\frac{1}{2}(a{x}^{-1}+bx)\},$$

and $K_{\lambda }(·)$ presents the modified Bessel function of the third kind in Barndorff-Nielsen and Shephard [42].

4.2. EM Algorithm to Compute the Posterior Modes

To implement the IBF sampler, we first use the EM algorithm to find the observed posterior mode ${\mathbf{\theta }}^0=({\mathbf{\beta }}_1^0,{\mathbf{\beta }}_2^0,k^0)$ of $\pi ({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y})$ and get the optimal importance sampling density $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^0)$. The log-posterior based on complete data $(\mathbf{y},\mathbf{z})$ is that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill l({\mathbf{\beta }}_1,{\mathbf{\beta }}_2,k|\mathbf{y},\mathbf{z})& =-\frac{1}{2}\sum _{i=1}^k\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1-\theta z_i)}^2}{{\tau }^2{z}_i}\hfill \\ & -\frac{1}{2}\sum _{i=k+1}^n\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2-\theta z_i)}^2}{{\tau }^2{z}_i}-\sum _{i=1}^n{z}_i.\hfill \end{array}\end{array}$$

Based on Zhou et al. [34] and Tian et al. [35], we get

$$\begin{array}{ccc}& & {\gamma }_{ji}=E(z_i|y_i,{\mathbf{\beta }}_j)=\frac{{\tau }^2}{{\theta }^2+2{\tau }^2}+\frac{|y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_j|}{\sqrt{{\theta }^2+2{\tau }^2}},\hfill \\ & & {\delta }_{ji}=E(\frac{1}{z_i}|y_i,{\mathbf{\beta }}_j)=\frac{\sqrt{{\theta }^2+2{\tau }^2}}{|y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_j|},j=1,2.\hfill \end{array}$$

Suppose the initial value of $\mathbf{\theta }$ is $({\mathbf{\beta }}_1^{(0)},{\mathbf{\beta }}_2^{(0)},k^{(0)})$, and the $(t-1)$-th iteration values is ${\mathbf{\theta }}^{(t-1)}=({\mathbf{\beta }}_1^{(t-1)},{\mathbf{\beta }}_2^{(t-1)},k^{(t-1)})$, the EM algorithm in the t-th iteration involves two steps:

E step: Omitting the terms without relation to $\mathbf{\theta }$, the Q function in the E step is the conditional expectation of $l(\mathbf{\theta }|\mathbf{y},\mathbf{z})$ with respect to the conditional distribution $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^{(t-1)})$, which is presented as

$$\begin{array}{ccc}& & Q(\mathbf{\theta }|\mathbf{y},{\mathbf{\theta }}^{(t-1)})=E[l(\mathbf{\theta }|\mathbf{y},\mathbf{z})|\mathbf{y},{\mathbf{\theta }}^{(t-1)}]\hfill \\ & & =-\frac{1}{{\tau }^2}[\sum _{i=1}^k{\delta }_{1i}^{(t-1)}{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1)}^2-2\theta \sum _{i=1}^k(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1)\hfill \\ & & +{\theta }^2\sum _{i=1}^k{\gamma }_{1i}^{(t-1)}]-\frac{1}{{\tau }^2}[\sum _{i=k+1}^n{\delta }_{2i}^{(t-1)}{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2)}^2\hfill \\ & & -2\theta \sum _{i=k+1}^n(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2)+{\theta }^2\sum _{i=k+1}^n{\gamma }_{2i}^{(t-1)}]\hfill \\ & & -[\sum _{i=1}^k{\gamma }_{1i}^{(t-1)}+\sum _{i=k+1}^n{\gamma }_{2i}^{(t-1)}],\hfill \end{array}$$

with

$$\begin{array}{ccc}& & {\gamma }_{ji}^{(t-1)}=\frac{{\tau }^2}{{\theta }^2+2{\tau }^2}+\frac{|y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_j^{(t-1)}|}{\sqrt{{\theta }^2+2{\tau }^2}},\hfill \\ & & {\delta }_{ji}^{(t-1)}=\frac{\sqrt{{\theta }^2+2{\tau }^2}}{|y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_j^{(t-1)}|},j=1,2.\hfill \end{array}$$

M step: We maximize the Q function with respect to $\mathbf{\theta }$. For fixed k, let ${\mathbf{W}}_1^{(t-1)}=\mathrm{diag}({\delta }_{11}^{(t-1)},\cdots ,{\delta }_{1k}^{(t-1)}),$ and

$${\mathbf{W}}_2^{(t-1)}=\mathrm{diag}({\delta }_{2,k+1}^{(t-1)},\cdots ,{\delta }_{2,n}^{(t-1)}),$$

then we get

$${\mathbf{\beta }}_1^{(t)}(k)={({\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{(t-1)}{\mathbf{X}}_1)}^{-1}({\mathbf{X}}_1^{\prime }{\mathbf{W}}_1^{(t-1)}{\mathbf{y}}_1-\theta {\mathbf{1}}_k),$$

$${\mathbf{\beta }}_2^{(t)}(k)={({\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{(t-1)}{\mathbf{X}}_2)}^{-1}({\mathbf{X}}_2^{\prime }{\mathbf{W}}_2^{(t-1)}{\mathbf{y}}_2-\theta {\mathbf{1}}_{n-k}),$$

with ${\mathbf{1}}_k$ indicating the k dimension vector with each element equal 1. By denoting the convergence values as ${\widehat{\mathbf{\beta }}}_1(k)$ and ${\widehat{\mathbf{\beta }}}_2(k)$, we can obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \widehat{k}& =\arg \underset{p\le k\le n-p}{max}l({\widehat{\mathbf{\beta }}}_1(k),{\widehat{\mathbf{\beta }}}_2(k),k|\mathbf{y})\hfill \\ & =\arg \underset{p\le k\le n-p}{min}\{S_1(k)+S_2(k)\},\hfill \end{array}\end{array}$$

with

$$S_1(k)=\sum _{i=1}^k{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathbf{\beta }}}_1(k)),S_2(k)=\sum _{i=k+1}^n{\rho }_q(y_i-{\mathbf{x}}_i^{\prime }{\widehat{\mathbf{\beta }}}_2(k)).$$

and then, we use samples $(y_i,{\mathbf{x}}_i),i=1,\cdots ,\widehat{k}$ and $(y_i,{\mathbf{x}}_i),i=\widehat{k}+1,\cdots ,n$ to obtain the posterior mode ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$ via another EM algorithm, respectively. After convergence, we get ${\mathbf{\theta }}^0=({\mathbf{\beta }}_1^0,{\mathbf{\beta }}_2^0,k^0)$, the mode of posterior distribution $\pi (\mathbf{\theta }|\mathbf{y})$.

4.3. IBF Sampler

Based on the following expression

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \pi (\mathbf{\theta }|\mathbf{y})& =\int \pi (\mathbf{\theta }|\mathbf{y},\mathbf{z})\pi (\mathbf{z}|\mathbf{y})\mathrm{d}\mathbf{z}\hfill \\ & =\int \pi ({\mathbf{\beta }}_1|\mathbf{y},\mathbf{z},k)\pi ({\mathbf{\beta }}_2|\mathbf{y},\mathbf{z},k)\pi (k|\mathbf{y},\mathbf{z})\pi (\mathbf{z}|\mathbf{y})\mathrm{d}\mathbf{z},\hfill \end{array}\end{array}$$

in order to get samples from $\pi (\mathbf{\theta }|\mathbf{y})$, we first need to draw samples from $\pi (\mathbf{z}|\mathbf{y})$. The IBF sampler proposed by Tan et al. [20,21] shows that

$$\pi (\mathbf{z}|\mathbf{y})\propto \frac{\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^0)}{\pi ({\mathbf{\theta }}^0|\mathbf{y},\mathbf{z})},$$

where $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^0)$ is an importance sampling density or proposal density in the sampling/importance resampling algorithm, with ${\mathbf{\theta }}^0$ denoting the posterior mode of $\pi (\mathbf{\theta }|\mathbf{y})$ obtained by EM algorithm in Section 4.2. The IBF sampling algorithm follows by

Step 1. Based on (14), for $i=1,\cdots ,k^0$, draw i.i.d. samples

$${\{z_i^{(j)}\}}_{j=1}^J\sim \mathrm{GIG}(\frac{1}{2},\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_1^0)}^2}{{\tau }^2},\frac{2{\tau }^2+{\theta }^2}{{\tau }^2}),$$

and based on (15), for $i=k^0+1,\cdots ,n$, draw i.i.d. samples

$${\{z_i^{(j)}\}}_{j=1}^J\sim \mathrm{GIG}(\frac{1}{2},\frac{{(y_i-{\mathbf{x}}_i^{\prime }{\mathbf{\beta }}_2^0)}^2}{{\tau }^2},\frac{2{\tau }^2+{\theta }^2}{{\tau }^2}).$$

Denote ${\mathbf{z}}^{(j)}={(z_1^{(j)},z_2^{(j)},\cdots ,z_n^{(j)})}^{\prime }$, then ${\{{\mathbf{z}}^{(j)}\}}_{j=1}^J$ are the samples from $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^0).$

Step 2. Calculate the weights

$${\lambda }_j=\frac{\frac{1}{\pi ({\mathbf{\theta }}^0|\mathbf{y},{\mathbf{z}}^{(j)})}}{{\sum }_{l=1}^J\frac{1}{\pi ({\mathbf{\theta }}^0|\mathbf{y},{\mathbf{z}}^{(l)})}},j=1,2\cdots ,J,$$

where

$$\begin{array}{c}\hfill \pi ({\mathbf{\theta }}^0|\mathbf{y},{\mathbf{z}}^{(j)})=\pi ({\mathbf{\beta }}_1^0|\mathbf{y},{\mathbf{z}}^{(j)},k^0)\pi ({\mathbf{\beta }}_2^0|\mathbf{y},{\mathbf{z}}^{(j)},k^0)\pi (k^0|\mathbf{y},{\mathbf{z}}^{(j)}),\end{array}$$

is obtained from (8), (10) and (12) with $\mathbf{z}$ replaced by ${\mathbf{z}}^{(j)}$.

Step 3. Choose a subset from ${\{{\mathbf{z}}^{(j)}\}}_{j=1}^J$ via resampling without replacement from the Step discrete distribution on ${\{{\mathbf{z}}^{(j)}\}}_{j=1}^J$ with probabilities ${\{{\lambda }_j\}}_{j=1}^J$ to obtain an i.i.d. sample of size $L(<J)$ approximately from $\pi (\mathbf{z}|\mathbf{y})$, denoted by ${\{{\mathbf{z}}^{(j_l)}\}}_{l=1}^L.$

Step 4. Denote

$${\mathbf{z}}_1^{(j_l)}={(z_1^{(j_l)},\cdots ,z_k^{(j_l)})}^{\prime },{\mathbf{z}}_2^{(j_l)}={(z_{k+1}^{(j_l)},\cdots ,z_n^{(j_l)})}^{\prime },$$

$${\mathbf{W}}_1^{(j_l)}=\mathrm{diag}({\tau }^2{\mathbf{z}}_1^{(j_l)}),{\mathbf{W}}_2^{(j_l)}=\mathrm{diag}({\tau }^2{\mathbf{z}}_2^{(j_l)}),$$

$${\mathbf{b}}_1^{(j_l)}={({{\mathbf{X}}_1}^{\prime }{{\mathbf{W}}_1^{(j_l)}}^{-1}{\mathbf{X}}_1)}^{-1}{{\mathbf{X}}_1}^{\prime }{{\mathbf{W}}_1^{(j_l)}}^{-1}({\mathbf{y}}_1-\theta {\mathbf{z}}_1^{(j_l)}),$$

$${\mathbf{b}}_2^{(j_l)}={({{\mathbf{X}}_2}^{\prime }{{\mathbf{W}}_2^{(j_l)}}^{-1}{\mathbf{X}}_2)}^{-1}{{\mathbf{X}}_2}^{\prime }{{\mathbf{W}}_2^{(j_l)}}^{-1}({\mathbf{y}}_2-\theta {\mathbf{z}}_2^{(j_l)}),$$

and

$$S_1^{(j_l)}={({\mathbf{y}}_1-\theta {\mathbf{z}}_1^{(j_l)}-{\mathbf{X}}_1{\mathbf{b}}_1^{(j_l)})}^{\prime }{{\mathbf{W}}_1^{(j_l)}}^{-1}({\mathbf{y}}_1-\theta {\mathbf{z}}_1^{(j_l)}-{\mathbf{X}}_1{\mathbf{b}}_1^{(j_l)}),$$

$$S_2^{(j_l)}={({\mathbf{y}}_2-\theta {\mathbf{z}}_2^{(j_l)}-{\mathbf{X}}_2{\mathbf{b}}_2^{(j_l)})}^{\prime }{{\mathbf{W}}_2^{(j_l)}}^{-1}({\mathbf{y}}_2-\theta {\mathbf{z}}_2^{(j_l)}-{\mathbf{X}}_2{\mathbf{b}}_2^{(j_l)}).$$

For $k=p,p+1,\cdots ,n-p$, let

$$w_k^l=\frac{|{\mathbf{X}}_1^{\prime }{{\mathbf{W}}_1^{j_l}}^{-1}{\mathbf{X}}_1{|}^{-1/2}{|{\mathbf{X}}_2^{\prime }{{\mathbf{W}}_2^{j_l}}^{-1}{\mathbf{X}}_2|}^{-1/2}}{exp(\frac{S_1^{(j_l)}+S_2^{(j_l)}+2{\sum }_{i=1}^n{z}_i^{j_l}}{2})},$$

and

$$p_k^l=\frac{w_k^l}{{\sum }_{k=p}^{n-p}{w}_k^l},k=p,\cdots ,n-p.$$

For $l=1,\cdots ,L$, sample $k^{(l)}$ from $\{p,\cdots ,n-p\}$ with probabilities $\{p_k^l,k=p,\cdots ,n-p\}$ based on (12).

Step 5. Denote

$${\mathbf{X}}_1^l={({\mathbf{x}}_1,\cdots ,{\mathbf{x}}_{k^{(l)}})}^{\prime },{\mathbf{X}}_2^l={({\mathbf{x}}_{k^{(l)}+1},\cdots ,{\mathbf{x}}_n)}^{\prime },$$

$${\mathbf{y}}_1^l={(y_1,\cdots ,y_{k^{(l)}})}^{\prime },{\mathbf{y}}_2^l={(y_{k^{(l)}+1},\cdots ,y_n)}^{\prime },$$

$${\mathbf{z}}_1^{(l)}={(z_1^{(j_l)},\cdots ,z_{k^{(l)}}^{(j_l)})}^{\prime },{\mathbf{z}}_2^{(l)}={(z_{k^{(l)}+1}^{(j_l)},\cdots ,z_n^{(j_l)})}^{\prime },$$

$${\mathbf{W}}_1^{(l)}=\mathrm{diag}({\mathbf{z}}_1^{(l)}),{\mathbf{W}}_2^{(l)}=\mathrm{diag}({\mathbf{z}}_2^{(l)}),$$

$${\mathbf{b}}_1^{(l)}={({{\mathbf{X}}_1^l}^{\prime }{{\mathbf{W}}_1^{(l)}}^{-1}{\mathbf{X}}_1^l)}^{-1}{{\mathbf{X}}_1^l}^{\prime }{{\mathbf{W}}_1^{(l)}}^{-1}({\mathbf{y}}_1^l-\theta {\mathbf{z}}_1^{(l)}),$$

$${\mathbf{b}}_2^{(l)}={({{\mathbf{X}}_2^l}^{\prime }{{\mathbf{W}}_2^{(l)}}^{-1}{\mathbf{X}}_2^l)}^{-1}{{\mathbf{X}}_2^l}^{\prime }{{\mathbf{W}}_2^{(l)}}^{-1}({\mathbf{y}}_2^l-\theta {\mathbf{z}}_2^{(l)}),$$

for $l=1,\cdots ,L$, based on (9), generate

$${\mathbf{\beta }}_1^{(l)}\sim \mathbb{N}\left({\mathbf{b}}_1^{(l)},{({{\mathbf{X}}_1^l}^{\prime }{{\mathbf{W}}_1^{(l)}}^{-1}{\mathbf{X}}_1^l)}^{-1}\right),$$

and based on (11), generate

$${\mathbf{\beta }}_2^{(l)}\sim \mathbb{N}\left({\mathbf{b}}_2^{(l)},{({{\mathbf{X}}_2^l}^{\prime }{{\mathbf{W}}_2^{(l)}}^{-1}{\mathbf{X}}_2^l)}^{-1}\right).$$

Denote ${\mathbf{\theta }}^{(l)}=({\mathbf{\beta }}_1^{(l)},{\mathbf{\beta }}_2^{(l)},k^{(l)})$, then, ${\{{\mathbf{\theta }}^{(l)}\}}_{l=1}^L$ are the i.i.d. samples approximately from $\pi (\mathbf{\theta }|\mathbf{y})$, which can be used for statistical inferences immediately without convergence diagnosis.

5. Simulation Studies

In this section, we conduct simulations to investigate the performance of the proposed approach to detect the position of the change-point by comparing it with the normal mean change-point detecting method proposed by Chen [1]. We generate $n=200$ observations with a change-point at position k from the following model:

$$\begin{array}{c}\hfill \begin{array}{cc}& y_i={\beta }_1{x}_{1i}+{\beta }_2{x}_{2i}+{\epsilon }_i,i=1,\cdots ,k,\hfill \\ & y_i={\beta }_1^{*}{x}_{1i}+{\beta }_2^{*}{x}_{2i}+{\epsilon }_i,i=k+1,\cdots ,n.\hfill \end{array}\end{array}$$

where ${\mathbf{\beta }}_1={({\beta }_1,{\beta }_2)}^{\prime }={(1,1)}^{\prime }$, ${\mathbf{\beta }}_2={({\beta }_1^{*},{\beta }_2^{*})}^{\prime }={(2,3)}^{\prime }$, and $x_{ji}\sim U(-2,2)$ for $j=1,2$ and $i=1,\cdots ,n$. Here k is set as 40, 60, 80, 100, 120, 140, and 160, respectively, and the following four error distributions are considered:

• Normal distribution: ${\epsilon }_i\sim N(0,1)$,
• Laplace distribution: ${\epsilon }_i\sim L(0,1)$,
• t distribution with three degree of freedom: ${\epsilon }_i\sim t(3)$,
• Cauchy distribution.

In these distributions, cases 2–4 are heavy-tailed and often used to describe the outlier situations, and especially case 4, the Cauchy distribution is extremely heavy-tailed. The quantile levels to estimate the position of change-point are set as 0.25, 0.5 and 0.75, and the experiment is replicated 500 times. In each replication, firstly, a position $k^{(0)}=100$ and the related least square estimators ${\mathbf{\beta }}_1^{(0)}$ and ${\mathbf{\beta }}_2^{(0)}$ based on the first 100 observations and later 100 observations respectively are used as initial values to implement the EM algorithm with aim to obtain the posterior mode ${\mathbf{\theta }}^0=({\mathbf{\beta }}_1^0,{\mathbf{\beta }}_2^0,k^0)$, which is then used in the IBF sampling to get an importance sampling density $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^0)$. Secondly, ${\{{\mathbf{z}}^{(j)}\}}_{j=1}^J$ are generated from $\pi (\mathbf{z}|\mathbf{y},{\mathbf{\theta }}^0)$ with $J=6000$ by SIR algorithm, and finally, $L=3000$ i.i.d. posterior samples ${\{{\mathbf{\theta }}^{(l)}\}}_{l=1}^L$ are generated from $\pi (\mathbf{\theta }|\mathbf{y})$. The estimation of change-point position k is denoted as ${\widehat{k}}^{(i)}$ in the i-th replication, which is the mean of L posterior samples of k, and the final estimate of k is calculated as $\widehat{k}=\frac{1}{500}{\sum }_{i=1}^{500}{\widehat{k}}^{(i)}$. For comparison purposes, the absolute deviance (AD), the root of the mean square error (RMSE), and mean absolute deviance (MAD) of the estimated change-point position k are given as

$$\mathrm{AD}=|\widehat{k}-k|,\mathrm{MAD}=\frac{1}{500}\sum _{i=1}^{500}|{\widehat{k}}^{(i)}-k|,$$

and

$$\mathrm{RMSE}=\sqrt{\frac{1}{500}\sum _{i=1}^{500}{({\widehat{k}}^{(i)}-k)}^2}.$$

The results based on the normal method proposed by Chen [1] are also recorded.

Table 1. Simulation results.

Error	k	Normal Method			$\mathit{q}=0.25$			$\mathit{q}=0.5$			$\mathit{q}=0.75$
		AD	MAD	RMSE	AD	MAD	RMSE	AD	MAD	RMSE	AD	MAD	RMSE
N(0,1)	40	0.27	2.42	4.36	0.37	3.91	7.51	0.06	3.62	6.32	0.83	4.92	10.66
	60	0.45	2.28	3.69	1.71	4.52	9.47	1.51	4.64	9.48	1.50	4.23	8.75
	80	0.50	2.42	4.06	0.12	4.00	7.16	0.17	4.23	7.68	0.95	5.24	12.98
	100	0.24	1.87	3.12	0.10	3.88	8.10	0.04	3.86	8.08	0.08	3.98	7.84
	120	0.09	2.44	3.92	0.15	3.48	5.35	0.18	3.24	4.98	0.08	3.70	5.95
	140	0.13	2.68	5.35	0.56	4.32	9.56	1.34	4.85	10.52	1.54	4.97	10.65
	160	0.35	1.99	3.25	0.32	4.66	8.66	0.52	4.68	7.98	0.35	4.46	7.98
t(3)	40	8.92	14.36	35.26	6.19	11.23	29.04	5.86	10.46	27.69	5.1	9.89	25.68
	60	0.61	10.27	20.41	0.68	7.88	14.08	0.50	7.74	14.54	0.24	8.08	14.85
	80	0.54	9.88	21.95	0.68	6.84	16.31	0.81	5.01	9.11	0.9	4.96	9.06
	100	1.64	7.5	18.83	1	5.5	11.08	0.35	6.27	13.83	0.77	5.25	11.47
	120	3.39	9.99	23.12	1.01	7.07	16.39	1	6.18	15.06	1.08	5.98	15.01
	140	4.61	9.71	23.68	0.25	4.99	8.38	0.25	5.29	10.22	0.59	5.21	10.11
	160	4.95	11.29	25.36	1.37	6.27	11.63	1.42	6.35	12.80	1.58	6.65	13.67
t(1)	40	54.8	70.6	88.4	12.80	27.77	53.88	13.34	28.80	55.08	12.6	28.58	54.76
	60	38.83	62.03	76.41	11.00	29.49	54.00	11.43	29.54	54.02	11.05	29.61	53.60
	80	14.08	55.33	64.97	0.43	22.15	42.45	0.45	21.52	41.43	1.38	21.44	41.26
	100	8.86	56.28	63.54	1.16	17.78	35.79	1.40	17.60	35.87	1.33	17.13	34.79
	120	21.5	59.33	69.33	4.40	17.10	36.49	4.37	18.81	38.29	5.09	18.82	38.75
	140	38.76	63.06	76.94	12.68	24.03	48.07	11.55	23.68	47.98	11.55	23.64	47.36
	160	70.27	80.9	98.16	12.64	30.20	57.57	12.54	29.72	57.35	11.51	29.05	56.14
L(0,1)	40	2.42	6.02	16.11	0.81	4.83	12.17	0.26	3.6	6.59	0.6	4.14	6.95
	60	1.54	5.3	11.10	0.34	4.3	7.50	0.34	4.37	7.56	0.06	4.28	7.14
	80	0.79	4.45	9.80	0.45	3.94	7.14	0.11	3.89	6.40	0.03	3.98	6.54
	100	0.88	5.04	10.57	0.26	3.89	7.18	0.16	3.97	7.40	0.09	4.00	7.06
	120	0.01	4.21	7.42	0.12	4.12	7.02	0.41	4.22	7.28	0.07	4.02	6.82
	140	1.69	5.71	13.42	1.03	5.40	9.85	0.61	5.14	9.05	0.97	5.14	8.96
	160	2.48	6.56	17.97	2.32	6.94	18.73	2.13	6.57	16.80	1.43	6.44	16.29

presents the results. It is not surprising that when the errors come from $N(0,1)$ distribution, the estimators based on the normal model to detect the position of the change-point are more accurate than the ones via quantile regression method, with smaller AD, MAD and RMSE. However, in the heavy-tailed $t(3)$, $L(0,1)$ and Cauchy cases, overwhelmingly, the estimators based on the quantile regression approach with IBF sampler are considerately superior to the ones according to the normal method. Specifically, in the extremely heavy-tailed Cauchy situation, the normal method simply uses values around 100 to estimate the position of the change-point whenever the real k is small or large; on the contrary, the performance of the quantile regression method is relatively much better. Clearly, the quantile regression method based on the IBF sampler is more robust than the usual normal method to estimate the position of regression coefficient change-point when the data are from heavy-tailed distributions.

6. American Stock Market Data Analysis

Holbert [43] conducted a Bayesian analysis of a switching linear model and analyzed the Stock Exchange data in America, which includes $n=35$ monthly observations from January 1967 to November 1969. Response variable Y is the monthly dollar volume of sales (in millions) on the Boston Stock Exchange (BSE), and regressor X is the combined New York American Stock Exchange (NYAMSE). The scatter plot shows that a suspicious change-point may exist at some time position, so a regression model with a change-point may be better to fit this data than a simple linear regression model. Chen [1] analyzed this data using the normal change-point model and detected a possible position of change-point, which is $k=23$ with minimal SIC=358.18. Therefore, Chen [1] believed the regression model change starts at the time point 24, December 1968. Osorio and Galea [8] studied this data based on the Student-t model with the SIC technique, finding that $k=9$ for $t(1)$ and $k=23$ when the degree of freedom is bigger than 4. Osorio and Galea [8] finally drew the same conclusion as Chen [1] that the regression model change starts at time point 24. Yang [12] detected change-point position 9 based on the Laplace regression model with a minimal SIC=353.83. Carefully examining the scatter plot, it can be seen that point 22 at the top-right is far away from other points, and $y_{22}$ is much bigger than others; thus, point 22 appears as an outlier. We consider the following robust QR change-point model

$$\begin{array}{c}\hfill \begin{array}{cc}& y_i={\beta }_{10}+{\beta }_{11}{x}_i+{ϵ}_i,i=1,\cdots ,k,\hfill \\ & y_i={\beta }_{20}+{\beta }_{21}{x}_i+{ϵ}_i,i=k+1,\cdots ,n.\hfill \\ & {\epsilon }_i\sim \mathrm{ALD}(0,1,q),i=1,2,\cdots ,n.\hfill \end{array}\end{array}$$

Using the above QR change-point model with 3000 generated IBF samples, we find out that the change-point possibly occurs at $k=9$ for quantile level q inside interval [0.45, 0.78], and $k=23$ for q outside this interval.

Figure 1 presents the scatter plots and estimated regression lines based on the normal model and QR change-point model with five quantile levels. Figure 1a presents the result by the normal model with detected change-point position 23, clearly point 22 has a high effect on the estimation, for that point 22 increases the slope of the black line by dragging this line towards itself. Figure 1b,c,f show the estimated quantile regression lines with relatively low level $q=0.2$, $q=0.4$ and high level $q=0.8$ respectively, and the estimated change-points are all 23, but the recognized regression lines can not effectively reflect the whole dependent relationship between BSE and NYAMSE. Figure 1d,e gives the estimated median and level $q=0.6$ regression lines with change-point 9, which are robust to outlier 22 and reasonable in describing the relationship between the response and regressor. Therefore, we consider time point 9 as the position of the change-point and believe that the regression model change starts at time point 10, which is October of 1967. The median and $q=0.6$ quantile regression models with change-point 9 are both reasonable. Take $q=0.6$ for example, in Figure 1e, the estimated regression coefficients (black solid line) are ${\widehat{\beta }}_{10}=17.93,{\widehat{\beta }}_{11}=0.52$ for observations 1–9, and ${\widehat{\beta }}_{20}=-32.801,{\widehat{\beta }}_{21}=1.213$ for observations 10–35 (red dashed line). This example illustrates that the quantile regression model with the IBF sampler effectively detects the position of possible change-point. As a by-product, the density histograms of the posterior samples generated by the IBF sampler from the observed posterior distribution of parameters are presented in Figure 2, which shows us how the parameters are distributed given the observed data graphically.

7. Discussion

In this paper, we proposed a change-point detection procedure based on a quantile regression model and developed a non-iterative sampling algorithm to estimate the position of the change-point. We investigated the performance of the procedure through simulations under different error distributions, finding that estimators are more accurate than ones based on the normal method when the errors follow heavy-tailed distributions, especially when the true position of the change-point is in the head part or in the tail part of the data. Finally, we applied our method to the Holbert data and detected a change-point different from the result based on the normal method. The procedure combines the robustness of quantile regression and the computational effectiveness of the non-iterative IBF sampler, so the estimations are more credible and effective than the ones based on the normal and likelihood methods. All work is done using the free R package version 4.2.2 [44]. The current work can be extended to multiple change-point problems, and further research is needed in the future.