# Robust Mean Change-Point Detecting through Laplace Linear Regression Using EM Algorithm

Authors: Fengkai Yang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We proposed a robust mean change-point estimation algorithm in linear regression with the assumption that the errors follow the Laplace distribution. By representing the Laplace distribution as an appropriate scale mixture of
normal distribution, we developed the expectation maximization (EM) algorithm to estimate the position of mean change-point. We investigated the performance of the algorithm through different simulations, finding that our methods is robust to the distributions of errors and is effective to estimate the position of mean change-point. Finally, we applied our method to the classical Holbert data and detected a change-point.


---

## Body

## 1. Introduction

Change-point analysis has been an active research area since the early 1950s. During the following period of sixty-some years, numerous articles have been published in various journals and proceedings. Chen and Gupta [1] summed up the main methods and applications on change-point detection and estimation. Ever since the change-point hypothesis was introduced into statistical analyses, the study of switching regression models has taken place in regression analysis. This made some previously poorly fitted regression models better fitted to some datasets after the change-point has been located in the regression models.The Schwarz information criterion (SIC) proposed by Schwarz [2] has been applied to change-point analysis for different underlying models by many authors in the literature. Chen [3] was the first to use the SIC model selection method to study the mean change-point problem in normal linear regression model, and later Chen and Gupta [4] used the same method to detect both mean and variance change-point in normal model. Chen and Wang [5] developed a statistical change-point model approach for the detection of DNA copy number variations in array CGH data using the SIC method and assuming the error follows the normal distribution.However, in practice, we do not know the real distribution of the data, and it is difficult to determine the real distribution especially when some change-point is present in the data. So, the normal assumption is not always suitable, for a lot of real data usually shows heavy tail and skewness. In such cases, some robust change-point detecting model with heavy-tailed distribution might be better than the normal model. Osorio and Galea [6] developed a mean change-point linear regression model with independent errors distributed according to the Studentt-distribution and located the change-point using the SIC method. Lin et al. [7] studied the variance change-points in the Student t regression model under Bayesian framework and analyzed the U.S. stock market data.The symmetric Laplace distribution, also known as the double exponential distribution or the first law of Laplace, is another heavy-tailed error distribution besides the Studentt-distribution. It is less sensitive to the outlier and more robust than the normal distribution. Kotz et al. [8] gave a systematic overview of the Laplace distribution, in their book, and the authors were devoted to presenting properties, generalizations of the Laplace distribution, and their applications in communications, economics, engineering, and finance.In recent years, statistical models based on Laplace distribution have developed rapidly both in theory and application. Purdom and Holmes [9] found that the error distribution for gene expression data from microarray experiments can be better fitted by Laplace distribution than normal distribution. Pop [10] identified a Laplace distribution in the change in daily sunspot number, and later Noble and Wheatland [11] showed the physical origin of Laplace distribution and its use in daily sunspot numbers. van Sanden and Burzykowski [12] considered the analysis of microarray data by using ANOVA models under the assumption of Laplace-distributed error terms. Phillips [13] developed the expectation maximization (EM) algorithm for the least absolute deviation regression, which is also known as the Laplace regression or median regression. Park and Casella [14] interpreted the Lasso estimate for linear regression parameters as a Bayesian posterior mode estimate when the regression parameters have independent Laplace priors. Song et al. [15] proposed a robust estimation procedure for mixture linear regression models assuming that the error terms follow the Laplace distribution.In this paper, we study the single mean change-point problem in linear regression model assuming that the error follows the Laplace distribution via EM algorithm and use the SIC model selection method to estimate the position of the mean change-point. Then, we investigate the robustness of the algorithm through simulations under different error distributions. Finally, we apply our method to some stock market data set.


## 2. Laplace Linear Regression Model with Mean Change-Point

### 2.1. Laplace Distribution as Scale Mixture of the Normal Distribution

The symmetric Laplace distribution is commonly denoted byL(μ,σ), where μ∈(-∞,∞) is the location parameter and σ>0 is the scale parameter. The density function is given by
(1)fx=12σexp⁡⁡-x-μσ,x∈R,
with mean μ and variance 2σ2. The Laplace distribution is more peaked in the center and more heavy-tailed compared to the normal distribution, and, for a given samples X1,X2,…,Xn from L(μ,σ), the maximum likelihood estimate of μ is the median of the samples, which is robust to the outliers.Andrews and Mallows [16] presented necessary and sufficient conditions under which a random variable X may be generated as the ratio Z/V, where Z and V are independent and Z has a standard normal distribution. Random variable X is referred to as a normal variance mixture distribution or a scale mixture of Gaussian distribution. It was established that when V2/2 is exponential, X is double exponential, so Laplace is some kind of scale mixture of Gaussian distribution, where the mixing distribution is exponential distribution. The statement above can be described by the following proposition.Proposition 1 (representation of Laplace distribution).
Suppose random variableϵ~L(μ,σ), and then there exists a random variable ω~Γ(1,1/2), such that ϵ∣ω~N(μ,ωσ2), where Γ(α,β) denotes the Gamma distribution with parameters α and β.


### 2.2. Laplace Regression Model with Mean Change-Point

Letx1,y1,…,(xn,yn) be a sequence of observations obtained in a practical situation, where xi=(xi1,…,xip)T,  i=1,2,…,n, is a nonstochastic p vector variable, and the model we are going to discuss is
(2)yi=xiTβ+ϵi,i=1,2,…,n,
where β=(β1,…,βp)T is a p unknown parameter vector and ϵi,  i=1,…,n  are random errors, which are independent and identically distributed as L(0,σ), with σ unknown. We can see that yi, i=1,2,…,n, are independently distributed and, for i=1,2,…,n, yi~L(xiTβ,σ). To develop EM algorithm for the Laplace linear regression model, we utilize a mixture representation of the Laplace regression, which is given by Proposition 2, derived from Proposition 1.Proposition 2 (representation of Laplace linear regression).
Suppose random variablesyi~L(xiTβ,σ),  i=1,…,n, andthen there exist random latent variables ωi,  i=1,…,n, where ωi~Γ(1,1/2), such that yi∣ωi~N(xiTβ,ωiσ2), for i=1,2,…,n.Proposition2 represents the Laplace linear regression model as a normal mixture model with mixture variables following the Gamma distribution. This kind of representation is crucial to perform the EM algorithm in the following mean change-point model.The single mean (regression coefficients) change-point problem in a Laplacian linear regression model can be formulated as to test the following null hypothesis:(3)H0:yi=xiTβ+ϵi,i=1,2,…,n,
versus the alternative
(4)H1:yi=xiTβ1+ϵi,i=1,2,…,k,yi=xiTβ2+ϵi,i=k+1,…,n,
where
(5)β1=β0,β1,…,βp-1T,β2=β0*,β1*,…,βp-1*T,ϵi~L0,σ,i=1,2,…,n.
That is, a change exists (in the regression coefficients) in an unknown position k, denominated mean change-point. In the following section, we perform EM algorithm to estimate the unknown parameters and use the SIC model selection method to detect the position of change-point.


## 2.1. Laplace Distribution as Scale Mixture of the Normal Distribution

The symmetric Laplace distribution is commonly denoted byL(μ,σ), where μ∈(-∞,∞) is the location parameter and σ>0 is the scale parameter. The density function is given by
(1)fx=12σexp⁡⁡-x-μσ,x∈R,
with mean μ and variance 2σ2. The Laplace distribution is more peaked in the center and more heavy-tailed compared to the normal distribution, and, for a given samples X1,X2,…,Xn from L(μ,σ), the maximum likelihood estimate of μ is the median of the samples, which is robust to the outliers.Andrews and Mallows [16] presented necessary and sufficient conditions under which a random variable X may be generated as the ratio Z/V, where Z and V are independent and Z has a standard normal distribution. Random variable X is referred to as a normal variance mixture distribution or a scale mixture of Gaussian distribution. It was established that when V2/2 is exponential, X is double exponential, so Laplace is some kind of scale mixture of Gaussian distribution, where the mixing distribution is exponential distribution. The statement above can be described by the following proposition.Proposition 1 (representation of Laplace distribution).
Suppose random variableϵ~L(μ,σ), and then there exists a random variable ω~Γ(1,1/2), such that ϵ∣ω~N(μ,ωσ2), where Γ(α,β) denotes the Gamma distribution with parameters α and β.


## 2.2. Laplace Regression Model with Mean Change-Point

Letx1,y1,…,(xn,yn) be a sequence of observations obtained in a practical situation, where xi=(xi1,…,xip)T,  i=1,2,…,n, is a nonstochastic p vector variable, and the model we are going to discuss is
(2)yi=xiTβ+ϵi,i=1,2,…,n,
where β=(β1,…,βp)T is a p unknown parameter vector and ϵi,  i=1,…,n  are random errors, which are independent and identically distributed as L(0,σ), with σ unknown. We can see that yi, i=1,2,…,n, are independently distributed and, for i=1,2,…,n, yi~L(xiTβ,σ). To develop EM algorithm for the Laplace linear regression model, we utilize a mixture representation of the Laplace regression, which is given by Proposition 2, derived from Proposition 1.Proposition 2 (representation of Laplace linear regression).
Suppose random variablesyi~L(xiTβ,σ),  i=1,…,n, andthen there exist random latent variables ωi,  i=1,…,n, where ωi~Γ(1,1/2), such that yi∣ωi~N(xiTβ,ωiσ2), for i=1,2,…,n.Proposition2 represents the Laplace linear regression model as a normal mixture model with mixture variables following the Gamma distribution. This kind of representation is crucial to perform the EM algorithm in the following mean change-point model.The single mean (regression coefficients) change-point problem in a Laplacian linear regression model can be formulated as to test the following null hypothesis:(3)H0:yi=xiTβ+ϵi,i=1,2,…,n,
versus the alternative
(4)H1:yi=xiTβ1+ϵi,i=1,2,…,k,yi=xiTβ2+ϵi,i=k+1,…,n,
where
(5)β1=β0,β1,…,βp-1T,β2=β0*,β1*,…,βp-1*T,ϵi~L0,σ,i=1,2,…,n.
That is, a change exists (in the regression coefficients) in an unknown position k, denominated mean change-point. In the following section, we perform EM algorithm to estimate the unknown parameters and use the SIC model selection method to detect the position of change-point.


## 3. EM Algorithm and Schwarz Information Criterion

### 3.1. EM Algorithm underH0

DenoteX=(x1,…,xn)T, y=(y1,…,yn)T and random latent vector ω=(ω1,…,ωn)T. The likelihood function of y under H0 is
(6)Lβ,σ∣y=∏i=1n12σexp⁡⁡-yi-xiTβσ.
From the stochastic representation in Section 2.2, the regression model under H0 becomes
(7)yi∣ωi~NxiTβ,ωiσ2;ωi~Γ1,12,i=1,2,…,n.
So, the likelihood function of the complete data (y,ω) has the following form:
(8)Lβ,σ∣y,ω=∏i=1n12πωiσexp⁡⁡-yi-xiTβ22σ2ωi×12exp⁡⁡-12ωi.
Consequently, the complete log-likelihood function is
(9)lβ,σ∣y,ω=-n2log⁡⁡π-nlog⁡⁡2-n2log⁡⁡σ2-12σ2∑i=1nyi-xiTβ2ωi-12∑i=1nωi-12∑i=1nlog⁡⁡ωi.Given the initial valuesβ(0), σ(0), we can obtain the maximum likelihood estimate of β, σ based on EM algorithm via the following two steps.(i) E Step. Given the (t-1)th iteration values of (β(t-1),σ(t-1)) and omitting the terms having no relation to (β,σ), the Q function of the (t)th iteration is
(10)Qβ,σ∣y,βt-1,σt-1=Elβ,σ∣y,ω∣y,βt-1,σt-1=-n2ln⁡σ2-12σ2∑i=1nδiyi-xiTβ2,
where δi=E(1/ωi∣y,β(t-1),σ(t-1)).In order to obtainδi, we firstly compute the conditional probability density function (pdf) of f(ωi∣yi,β,σ). Note that the joint pdf of y and ω is
(11)fy,ω∣β,σ=∏i=1n12πωiσexp⁡⁡-yi-xiTβ22σ2ωi×12exp⁡⁡-12ωi.
Then, we have
(12)fω∣y,β,σ∝f(y,ω∣β,σ)∝∏i=1n‍ωi-1/2exp⁡-121ωiyi-xiTβ2σ2+ωi.
Due to the conditionally independent of ω1,ω2,…,ωn, given y, we obtain the marginal conditional pdf of ωi given yi as
(13)fωi∣yi,β,σ∝ωi1/2-1exp⁡⁡-121ωiyi-xiTβ2σ2+ωi~GIG12,yi-xiTβ2σ2,1,
where GIG(λ,a,b) is the generalized inverse Gaussian distribution and its pdf is given by
(14)f(x∣λ,a,b)=(b/a)λ/22Kλabxλ-1exp⁡⁡-12ax-1+bx,
where Kλ(·) is the modified Bessel function of the third kind in Barndorff-Nielsen and Shephard [17].Therefore, we can obtain(15)δi=E1ωi∣y,βt-1,σt-1=σ(t-1)yi-xiTβt-1.(ii)
M
Step. we maximize the Q function in E step with respect to β, σ.DenoteΣ=diag⁡(δ1,…,δn), and we get
(16)β(t)=XTΣX-1XTΣy,(17)σt=1n∑i=1nδiyi-xiTβt2-1/2.
It can be seen from (15) and (16) that the updating formulae of β is independent of the estimated value of σ; in other words, we need not to update σ at each EM iteration, and, after the final estimate β~ of β is obtained, the estimate of σ can then be gotten by maximizing the original log-likelihood function which produces
(18)σ~=1n∑i=1nyi-xiTβ~.
Finally, the maximum of the log-likelihood function under H0, is obtained as
(19)l0β~,σ~=-nlog⁡⁡2-nlog⁡⁡σ~-n=-nlog⁡⁡2n-nlog⁡⁡Sβ~-n
with
(20)Sβ~=∑i=1nyi-xiTβ~.


### 3.2. EM Algorithm underH1

DenoteX1=(x1,…,xk)T, y1=(y1,…,yk)T,  X2=(xk+1,…,xn)T, y2=(yk+1,…,yn)T, and random latent vectors ω1=(ω1,…,ωk)T, ω2=(ωk+1,…,ωn)T.The likelihood function ofy under H1 is
(21)Lβ1,β2,σ∣y=∏i=1k‍12σexp⁡⁡-yi-xiTβ1σ×∏i=k+1n‍12σexp⁡⁡-yi-xiTβ2σ.From the stochastic representation in Section2.2, the regression model under H1 becomes
(22)yi∣ωi~NxiTβ1,ωiσ2,i=1,…,k,yi∣ωi~NxiTβ2,ωiσ2,i=k+1,…,n,ωi~Γ1,12,i=1,2,…,n,
where ωi,i=1,…,n are independent.So, the likelihood function of the complete data(y,ω) under H1 has the following form:
(23)Lβ1,β2,σ∣y,ω=∏i=1k12πωiσexp⁡⁡-yi-xiTβ122σ2ωi×12exp⁡⁡-12ωi×∏i=k+1n12πωiσexp⁡⁡-yi-xiTβ122σ2ωi×12exp⁡⁡-12ωi.Consequently, the complete log-likelihood function is(24)lβ1,β2,σ∣y,ω=-n2log⁡⁡2π-nlog⁡⁡2-n2log⁡⁡σ2-12∑i=1n‍ωi-12∑i=1n‍log⁡⁡ωi-12σ2∑i=1k‍yi-xiTβ12ωi+∑i=k+1n‍yi-xiTβ22ωi.Given initial valuesβ1(0), β2(0), and σ(0), we can obtain the maximum likelihood estimate of β1, β2, and σ based on EM algorithm via the following two steps.(i) E Step. Given the (t-1)th iteration values of (β(t-1),σ(t-1)) and omitting the terms having no relation to (β,σ), the Q function of the (t)th iteration is
(25)Q(β1,β2,σ∣y,β1(t-1),β2(t-1),σ(t-1))=Elβ1,β2,σ∣y,ω∣y,β1t-1,β2t-1,σt-1=-n2log⁡σ2-12σ2∑i=1k‍γiyi-xiTβ12+∑i=k+1n‍γiyi-xiTβ22,
where γi=E(1/ωi∣y,β1(t-1),β2(t-1),σ(t-1)).The marginal conditional pdf ofωi given yi is that, for i=1,…,k,
(26)fωi∣yi,β1,σ∝ωi1/2-1exp⁡⁡-121ωiyi-xiTβ12σ2+ωi~GIG12,yi-xiTβ12σ2,1,
and, for i=k+1,…,n,
(27)fωi∣yi,β2,σ∝ωi1/2-1exp⁡⁡-121ωiyi-xiTβ22σ2+ωi~GIG12,yi-xiTβ22σ2,1.
Therefore, for i=1,…,k,
(28)γi=E1ωi∣y,β1t-1,σt-1=σ(t-1)yi-xiTβ1t-1,
and, for i=k+1,…,n,
(29)γi=E1ωi∣y,β2t-1,σt-1=σ(t-1)yi-xiTβ2t-1.(ii) M Step. We maximize the Q function in E step with respect to β1, β2, and σ.DenoteΣ1=diag⁡(γ1,…,γk), Σ2=diag⁡(γk+1,…,γn), and we get
(30)β1t=X1TΣ1X1-1X1TΣ1y1,β2(t)=X2TΣ2X2-1X2TΣ2y2,σ(t)=1n∑i=1k‍γiyi-xiTβ1t2+∑i=k+1n‍γiyi-xiTβ2t2-1/2.It can be seen from (28)–(30) that the updating formulae of β1, β2 is independent of the estimated value of σ, so we need not to update σ at each EM iteration. After the final estimates β^1 of β1 and β^2 of β2 are obtained, the estimate of σ can then be gotten by maximizing the original log-likelihood function, which produces
(31)σ^=1n∑i=1kyi-xiTβ^1+∑i=k+1nyi-xiTβ^2‍.
Finally, the maximum of the log-likelihood function under H1 is obtained as
(32)lkβ^1,β^2,σ^=-nlog⁡⁡2-nlog⁡⁡σ^-n=-nlog⁡⁡2n-nlog⁡⁡Sβ^1+Sβ^2-n
with
(33)Sβ^1=∑i=1kyi-xiTβ^1,Sβ^2=∑i=k+1nyi-xiTβ^2.


### 3.3. SIC Algorithm for Laplacian Regression Model

The work by Chen [3] proposes transforming the process of hypothesis testing in a procedure of model selection using the Schwarz Information Criterion (SIC) defined by
(34)SIC=-2lθ^+slog⁡⁡n,
where l(θ^) corresponds to the log-likelihood function evaluated on the maximum likelihood estimate of the parameters and s is the number of model parameters and n is the sample size. Note that maximizing the log-likelihood function is equivalent to minimizing the Schwarz information criterion.In the Laplacian linear regression model, the Schwarz information criterion underH0, denoted by SIC(n), is given by
(35)SICn=-2l0(β~,σ~)+p+1log⁡⁡n=2nlog⁡⁡Sβ~+2nlog⁡⁡2n+2n+p+1log⁡⁡n,
where l0(β~,σ~) and S(β~) are given by (19) and (20). The Schwarz information criterion under H1, denoted by SIC(k), is given by
(36)SICk=-2lk(β^1,β^2,σ^)+2p+1log⁡⁡n=2nlog⁡⁡Sβ^1+Sβ^2+2nlog⁡⁡2n+2n+2p+1log⁡⁡n
with lk(β^1,β^2,σ^), S(β^1), and S(β^2) given by (32) and (33).The selection criteria are to choose a model with a change-point in thek position, if, for some k,
(37)SIC(n)>SIC(k).
When the null hypothesis is rejected, the maximum likelihood estimate of the change-point in the regression coefficients, denoted by k^, must satisfy
(38)SICk^=min⁡⁡SICk:p≤k≤n-p=max⁡⁡lkβ^1,β^2,σ^:p≤k≤n-p=min⁡⁡Sβ^1+Sβ^2.


## 3.1. EM Algorithm underH0

DenoteX=(x1,…,xn)T, y=(y1,…,yn)T and random latent vector ω=(ω1,…,ωn)T. The likelihood function of y under H0 is
(6)Lβ,σ∣y=∏i=1n12σexp⁡⁡-yi-xiTβσ.
From the stochastic representation in Section 2.2, the regression model under H0 becomes
(7)yi∣ωi~NxiTβ,ωiσ2;ωi~Γ1,12,i=1,2,…,n.
So, the likelihood function of the complete data (y,ω) has the following form:
(8)Lβ,σ∣y,ω=∏i=1n12πωiσexp⁡⁡-yi-xiTβ22σ2ωi×12exp⁡⁡-12ωi.
Consequently, the complete log-likelihood function is
(9)lβ,σ∣y,ω=-n2log⁡⁡π-nlog⁡⁡2-n2log⁡⁡σ2-12σ2∑i=1nyi-xiTβ2ωi-12∑i=1nωi-12∑i=1nlog⁡⁡ωi.Given the initial valuesβ(0), σ(0), we can obtain the maximum likelihood estimate of β, σ based on EM algorithm via the following two steps.(i) E Step. Given the (t-1)th iteration values of (β(t-1),σ(t-1)) and omitting the terms having no relation to (β,σ), the Q function of the (t)th iteration is
(10)Qβ,σ∣y,βt-1,σt-1=Elβ,σ∣y,ω∣y,βt-1,σt-1=-n2ln⁡σ2-12σ2∑i=1nδiyi-xiTβ2,
where δi=E(1/ωi∣y,β(t-1),σ(t-1)).In order to obtainδi, we firstly compute the conditional probability density function (pdf) of f(ωi∣yi,β,σ). Note that the joint pdf of y and ω is
(11)fy,ω∣β,σ=∏i=1n12πωiσexp⁡⁡-yi-xiTβ22σ2ωi×12exp⁡⁡-12ωi.
Then, we have
(12)fω∣y,β,σ∝f(y,ω∣β,σ)∝∏i=1n‍ωi-1/2exp⁡-121ωiyi-xiTβ2σ2+ωi.
Due to the conditionally independent of ω1,ω2,…,ωn, given y, we obtain the marginal conditional pdf of ωi given yi as
(13)fωi∣yi,β,σ∝ωi1/2-1exp⁡⁡-121ωiyi-xiTβ2σ2+ωi~GIG12,yi-xiTβ2σ2,1,
where GIG(λ,a,b) is the generalized inverse Gaussian distribution and its pdf is given by
(14)f(x∣λ,a,b)=(b/a)λ/22Kλabxλ-1exp⁡⁡-12ax-1+bx,
where Kλ(·) is the modified Bessel function of the third kind in Barndorff-Nielsen and Shephard [17].Therefore, we can obtain(15)δi=E1ωi∣y,βt-1,σt-1=σ(t-1)yi-xiTβt-1.(ii)
M
Step. we maximize the Q function in E step with respect to β, σ.DenoteΣ=diag⁡(δ1,…,δn), and we get
(16)β(t)=XTΣX-1XTΣy,(17)σt=1n∑i=1nδiyi-xiTβt2-1/2.
It can be seen from (15) and (16) that the updating formulae of β is independent of the estimated value of σ; in other words, we need not to update σ at each EM iteration, and, after the final estimate β~ of β is obtained, the estimate of σ can then be gotten by maximizing the original log-likelihood function which produces
(18)σ~=1n∑i=1nyi-xiTβ~.
Finally, the maximum of the log-likelihood function under H0, is obtained as
(19)l0β~,σ~=-nlog⁡⁡2-nlog⁡⁡σ~-n=-nlog⁡⁡2n-nlog⁡⁡Sβ~-n
with
(20)Sβ~=∑i=1nyi-xiTβ~.


## 3.2. EM Algorithm underH1

DenoteX1=(x1,…,xk)T, y1=(y1,…,yk)T,  X2=(xk+1,…,xn)T, y2=(yk+1,…,yn)T, and random latent vectors ω1=(ω1,…,ωk)T, ω2=(ωk+1,…,ωn)T.The likelihood function ofy under H1 is
(21)Lβ1,β2,σ∣y=∏i=1k‍12σexp⁡⁡-yi-xiTβ1σ×∏i=k+1n‍12σexp⁡⁡-yi-xiTβ2σ.From the stochastic representation in Section2.2, the regression model under H1 becomes
(22)yi∣ωi~NxiTβ1,ωiσ2,i=1,…,k,yi∣ωi~NxiTβ2,ωiσ2,i=k+1,…,n,ωi~Γ1,12,i=1,2,…,n,
where ωi,i=1,…,n are independent.So, the likelihood function of the complete data(y,ω) under H1 has the following form:
(23)Lβ1,β2,σ∣y,ω=∏i=1k12πωiσexp⁡⁡-yi-xiTβ122σ2ωi×12exp⁡⁡-12ωi×∏i=k+1n12πωiσexp⁡⁡-yi-xiTβ122σ2ωi×12exp⁡⁡-12ωi.Consequently, the complete log-likelihood function is(24)lβ1,β2,σ∣y,ω=-n2log⁡⁡2π-nlog⁡⁡2-n2log⁡⁡σ2-12∑i=1n‍ωi-12∑i=1n‍log⁡⁡ωi-12σ2∑i=1k‍yi-xiTβ12ωi+∑i=k+1n‍yi-xiTβ22ωi.Given initial valuesβ1(0), β2(0), and σ(0), we can obtain the maximum likelihood estimate of β1, β2, and σ based on EM algorithm via the following two steps.(i) E Step. Given the (t-1)th iteration values of (β(t-1),σ(t-1)) and omitting the terms having no relation to (β,σ), the Q function of the (t)th iteration is
(25)Q(β1,β2,σ∣y,β1(t-1),β2(t-1),σ(t-1))=Elβ1,β2,σ∣y,ω∣y,β1t-1,β2t-1,σt-1=-n2log⁡σ2-12σ2∑i=1k‍γiyi-xiTβ12+∑i=k+1n‍γiyi-xiTβ22,
where γi=E(1/ωi∣y,β1(t-1),β2(t-1),σ(t-1)).The marginal conditional pdf ofωi given yi is that, for i=1,…,k,
(26)fωi∣yi,β1,σ∝ωi1/2-1exp⁡⁡-121ωiyi-xiTβ12σ2+ωi~GIG12,yi-xiTβ12σ2,1,
and, for i=k+1,…,n,
(27)fωi∣yi,β2,σ∝ωi1/2-1exp⁡⁡-121ωiyi-xiTβ22σ2+ωi~GIG12,yi-xiTβ22σ2,1.
Therefore, for i=1,…,k,
(28)γi=E1ωi∣y,β1t-1,σt-1=σ(t-1)yi-xiTβ1t-1,
and, for i=k+1,…,n,
(29)γi=E1ωi∣y,β2t-1,σt-1=σ(t-1)yi-xiTβ2t-1.(ii) M Step. We maximize the Q function in E step with respect to β1, β2, and σ.DenoteΣ1=diag⁡(γ1,…,γk), Σ2=diag⁡(γk+1,…,γn), and we get
(30)β1t=X1TΣ1X1-1X1TΣ1y1,β2(t)=X2TΣ2X2-1X2TΣ2y2,σ(t)=1n∑i=1k‍γiyi-xiTβ1t2+∑i=k+1n‍γiyi-xiTβ2t2-1/2.It can be seen from (28)–(30) that the updating formulae of β1, β2 is independent of the estimated value of σ, so we need not to update σ at each EM iteration. After the final estimates β^1 of β1 and β^2 of β2 are obtained, the estimate of σ can then be gotten by maximizing the original log-likelihood function, which produces
(31)σ^=1n∑i=1kyi-xiTβ^1+∑i=k+1nyi-xiTβ^2‍.
Finally, the maximum of the log-likelihood function under H1 is obtained as
(32)lkβ^1,β^2,σ^=-nlog⁡⁡2-nlog⁡⁡σ^-n=-nlog⁡⁡2n-nlog⁡⁡Sβ^1+Sβ^2-n
with
(33)Sβ^1=∑i=1kyi-xiTβ^1,Sβ^2=∑i=k+1nyi-xiTβ^2.


## 3.3. SIC Algorithm for Laplacian Regression Model

The work by Chen [3] proposes transforming the process of hypothesis testing in a procedure of model selection using the Schwarz Information Criterion (SIC) defined by
(34)SIC=-2lθ^+slog⁡⁡n,
where l(θ^) corresponds to the log-likelihood function evaluated on the maximum likelihood estimate of the parameters and s is the number of model parameters and n is the sample size. Note that maximizing the log-likelihood function is equivalent to minimizing the Schwarz information criterion.In the Laplacian linear regression model, the Schwarz information criterion underH0, denoted by SIC(n), is given by
(35)SICn=-2l0(β~,σ~)+p+1log⁡⁡n=2nlog⁡⁡Sβ~+2nlog⁡⁡2n+2n+p+1log⁡⁡n,
where l0(β~,σ~) and S(β~) are given by (19) and (20). The Schwarz information criterion under H1, denoted by SIC(k), is given by
(36)SICk=-2lk(β^1,β^2,σ^)+2p+1log⁡⁡n=2nlog⁡⁡Sβ^1+Sβ^2+2nlog⁡⁡2n+2n+2p+1log⁡⁡n
with lk(β^1,β^2,σ^), S(β^1), and S(β^2) given by (32) and (33).The selection criteria are to choose a model with a change-point in thek position, if, for some k,
(37)SIC(n)>SIC(k).
When the null hypothesis is rejected, the maximum likelihood estimate of the change-point in the regression coefficients, denoted by k^, must satisfy
(38)SICk^=min⁡⁡SICk:p≤k≤n-p=max⁡⁡lkβ^1,β^2,σ^:p≤k≤n-p=min⁡⁡Sβ^1+Sβ^2.


## 4. Simulation Studies

In this section, we investigate the performance of the proposed approach to detect mean change-point through simulations, and we compare our procedure with the change-point detecting procedure proposed by Chen [3], who assumed that the errors follow the normal distribution. A data set of n=200 observations is generated from the model
(39)yi=β1x1i+β2x2i+β3x3i+ϵi,i=1,…,k;yi=β1*x1i+β2*x2i+β3*x3i+ϵi,i=k+1,…,n,
where β1=(β1,β2,β3)T=(1,1,1)T, β2=(β1*,β2*,β3*)T=(2,3,4)T, and xji~U(-1,1) for j=1,…,3 and i=1,…,n. For demonstrating the performances of the algorithm in different cases, we choose k=40,60,80,100,120,140,160, respectively, and consider the following six error distributions:(i)
Simulation 1: normal distribution:ϵi~N(0,1);(ii)
Simulation 2: Laplace distribution:ϵi~L(0,1);(iii)
Simulation 3:t distribution with three degree of freedom: ϵi~t(3);(iv)
Simulation 4:χ2 distribution with three degree of freedom: ϵi~χ2(3);(v)
Simulation 5: log-normal distribution:ϵi~LN(0,1);(vi)
Simulation 6: Cauchy distribution.In order to evaluate the finite sample performance of the proposed method, 500 replications are conducted for different error distribution, respectively. In each replication, the initial values ofβ and β1, β2 are set to be their ordinary least squares estimates, and the EM procedure is stopped when the absolute differences between two successive values of maximum log-likelihood functions are less than a preassigned small number (e.g., 10-6). The final values of β and β1, β2 are taken as the median regression coefficient and the estimate of σ is given by (31). The position k of change-point is estimated by the SIC method denoted as k^(i) in the ith replication.Finally, the mean and standard difference of the estimated change-point positionk are given by
(40)k^=1500∑i=1500k^i,sd⁡=1499∑i=1500k^i-k^2.
In order to compare the results of our proposed method with those of the the normal method proposed in Chen [3], the same simulations are conducted using normal method. The results are presented in Tables 1 and 2, where k is the true value of change-point position in simulation, estimate.L and diff.L are the estimate of k and the absolute difference of the estimate and the true value of k, and finally, sd.L is the standard difference of the estimate. The result by normal method is denoted as estimate.N, diff.N, and sd.N, respectively. We compare diff.L with diff.N and sd.L with sd.N, and the smaller the better. From Tables 1 and 2, we have the following findings.(1)
In theN(0,1) and L(0,1) cases, both the Laplace method and normal method behave quite well with small diff. Besides that, the normal method is better than the Laplace method with smaller sd in the normal error case, and the result is inverse in the Laplace error case.(2)
In the skewχ2(3) and heavy-tailed t(3), LN(0,1), and Cauchy cases, the Laplace method is better than the normal method with smaller diff. and smaller sd.(3)
In theχ2(3) and LN(0,1) cases, when the true position k of change-point is in the middle of the data (80,100,120), both of the two methods have good performance, but when k is in the lower part of the data (40,60) or upper part (140,160), the normal method can not detect the change-point effectively with large diff and large sd.(4)
An interesting phenomenon appears in the case of Cauchy distribution, where the normal method simply uses integers around 100 to estimate the true values ofk whenever k is small or large; on the contrary, the performance of the Laplace method is better relatively.Table 1
Comparison of results for Laplace model and normal model based on 500 replications in Simulations 1–3.
Error distribution
k
Estimate.L
Estimate.N
diff.L
diff.N
sd.L
sd.N
N
(
0,1
)
40
40.07
40.11
0.07
0.11
1.92
1.77
60
59.85
59.81
0.14
0.18
1.83
1.93
80
79.92
79.91
0.07
0.08
1.88
1.85
100
99.92
99.90
0.07
0.10
1.78
1.74
120
120.00
119.97
0.00
0.02
1.78
1.74
140
139.87
139.90
0.12
0.09
1.64
1.53
160
159.92
159.85
0.08
0.14
1.86
1.77
L
(
0,1
)
40
39.98
38.89
0.02
0.11
2.50
2.96
60
59.94
60.08
0.06
0.08
2.47
2.83
80
80.06
80.14
0.06
0.14
2.95
3.19
100
100.21
99.99
0.21
0.01
2.83
3.18
120
120.03
120.07
0.03
0.07
2.23
2.91
140
140.13
140.23
0.13
0.23
2.29
2.52
160
160.03
159.81
0.03
0.19
2.55
5.29
t
(
3
)
40
40.14
40.79
0.14
0.79
2.54
11.26
60
59.96
61.10
0.04
1.10
2.51
13.22
80
80.22
80.80
0.22
0.80
2.56
12.67
100
100.05
99.32
0.05
0.68
2.99
12.68
120
120.03
119.53
0.03
0.47
2.99
10.06
140
140.11
139.51
0.11
0.49
2.82
12.15
160
159.67
159.27
0.33
0.73
7.54
11.82Table 2
Comparison of results for Laplace model and normal model based on 500 replications in Simulations 4–6.
Error distribution
k
Estimate.L
Estimate.N
diff.L
diff.N
sd.L
sd.N
χ
2
(
3
)
40
42.00
46.42
2.00
6.42
16.74
28.94
60
60.67
61.58
0.67
1.58
9.62
19.85
80
80.37
81.44
0.37
1.44
7.76
17.26
100
99.18
99.03
0.82
0.97
8.76
14.99
120
119.82
120.36
0.18
0.36
11.13
16.98
140
139.26
137.25
0.74
2.75
10.98
22.95
160
158.21
152.70
1.79
7.30
15.74
30.94
L
N
(
0,1
)
40
41.23
49.78
1.23
9.78
13.49
38.92
60
60.40
63.87
0.40
3.87
4.83
27.14
80
80.58
83.08
0.58
3.08
4.66
25.60
100
100.00
101.35
0.00
1.35
6.50
21.20
120
120.44
117.74
0.44
2.26
6.03
24.65
140
139.92
136.16
0.08
3.84
5.48
27.81
160
158.80
152.27
1.20
7.73
13.71
35.13
Cauchy
40
51.75
98.10
11.75
58.10
46.99
67.04
60
64.81
98.77
4.81
38.77
34.23
66.56
80
81.63
81.44
1.63
1.44
7.76
17.26
100
98.61
99.34
1.39
0.66
29.76
64.17
120
119.83
105.36
0.17
14.64
28.86
65.80
140
134.43
107.23
5.57
32.78
36.47
63.55
160
149.95
105.12
10.05
54.90
46.53
66.94Another simulations design withβ1=(1,1,1)T and β2=(1,2,3)T are conducted and the results are not presented in our paper. The results also justify our former findings; that is, our Laplace method is more robust than the normal method to estimate the position of mean change-point when the errors follow skew and heavy-tailed distributions, especially when the true position of change-point is in the head part or in the tail part of the data.


## 5. Application to Stock Market Data

Holbert [18] studied the switching simple linear regression models with changes in the coefficients from a Bayesian point of view. Later on, Chen [3] analyzed this data using the SIC method to detect changes in the mean for normal linear regression models, and Osorio and Galea [6] analyzed the same data also using the SIC method to detect mean change-point but assuming the errors that follow the Student t-distribution. In this section, we reanalyze the Holbert data using the proposed Laplace linear regression model to detect the mean change-point based on the SIC method. The monthly dollar volume of sales (in millions) on the Boston Stock Exchange (BSE) is considered as the response variable y, and the combined New York American Stock Exchange (NYAMSE) is considered as the regressor x. There are n=35 observations corresponding to the period between January 1967 and November 1969. The model we consider is that
(41)yi=β0+β1xi+ϵi,i=1,…,k,yi=β0*+β1*xi+ϵi,i=k+1,…,n,
where
(42)ϵi~L0,σ,i=1,2,…,n.
The computed SIC values are listed in Table 3 along with the original BSE and NYAMSE values given in Holbert [18]. The SIC values calculated using the normal method by Chen [3] are also listed in the last column.Table 3
SIC calculated by Laplace method and normal method for Holbert data.
Time point
Calendar month
NYAMSE
BSE
SIC Laplace
SIC normal
1
Jan. 1967
10581.6
78.8
—
—
2
Feb. 1967
10234.3
69.1
364.2829
368.5739
3
Mar. 1967
13299.5
87.6
363.3368
367.8817
4
Apr. 1967
10746.5
72.8
363.3110
367.7757
5
May 1967
13310.7
79.4
361.7661
366.4980
6
Jun. 1967
12835.5
85.6
359.5049
365.7947
7
Jul. 1967
12194.2
75.0
357.4092
364.8795
8
Aug. 1967
12860.4
85.3
354.8215
363.9410
9
Sep. 1967
11955.6
86.9
353.8327
363.5574
10
Oct. 1967
13351.5
107.8
354.3397
363.5818
11
Nov. 1967
13285.9
128.7
357.2243
364.6607
12
Dec. 1967
13784.4
134.5
360.1891
365.4162
13
Jan. 1968
16336.7
148.7
361.0186
365.3077
14
Feb. 1968
11040.5
94.2
362.3622
365.5670
15
Mar. 1968
11525.3
128.1
363.4817
366.6527
16
Apr. 1968
16056.4
154.1
364.0982
366.8008
17
May 1968
18464.3
191.3
362.8076
366.9825
18
Jun. 1968
17092.2
191.9
360.4358
367.2177
19
Jul. 1968
15178.8
159.6
359.1133
367.3715
20
Aug. 1968
12774.8
185.5
359.5155
368.4097
21
Sep. 1968
12377.8
178.0
359.5954
368.3030
22
Oct. 1968
16856.3
271.8
356.8656
363.5156
23
Nov. 1968
14635.3
212.3
355.3476
358.1847
24
Dec. 1968
17436.9
139.4
357.0433
361.1139
25
Jan. 1969
16482.2
106.0
360.8022
364.8916
26
Feb. 1969
13905.4
112.1
360.8808
365.1567
27
Mar. 1969
11973.7
103.5
360.4217
365.0086
28
Apr. 1969
12573.6
92.5
360.9285
365.3012
29
May. 1969
16566.8
116.9
363.5041
367.3072
30
Jun. 1969
13558.7
78.9
364.3949
368.2468
31
Jul. 1969
11530.9
57.4
362.7110
368.2235
32
Aug. 1969
11278.0
75.9
360.3198
367.7685
33
Sep. 1969
11263.7
109.8
362.6290
368.1350
34
Oct. 1969
15649.5
129.2
—
—
35
Nov. 1969
12197.1
115.1
358.0474
361.4956The bold SIC value in Table3 is the minimum SIC value, which corresponds to time point 9, using the Laplace method and time point 23 using the normal method in Chen [3]. So, Chen [3] believed the regression model change starts at the time point 24, which is December 1968. Osorio and Galea [6] obtained the minimum values of SIC using the Student-t model, which are at time point 9 for degree of freedom 1 and 23 for degree of freedom 4. Osorio and Galea [6] finally drew the same conclusion as Chen [3] that the regression model change starts at the time point 24. We list the minimum values of SIC in different regression models in Table 4. Table 4 tells us that the minimum SIC value is 353.8327 for time point 9 in our Laplace model; hence, we draw the conclusion that the regression model change starts at the time point 10, which is October 1967. The median estimates for the parameters are β^0=10.500, β^1=0.0058 for observations 1–9 and β^0*=-37.646, β^1*=0.0119 for observations 10–35. Finally, the predictive regression lines are plotted in Figure 1. The scatter plot shows that there are outliers in the response observations, so the median regression line estimated form Laplace model is more robust than the ordinary least square regression line estimated from the normal model.Table 4
Minimum SIC values in Holbert data based on different models.
Method
Laplace
Normal
t(4)
t(1)
SIC
(
n
)
358.0474
361.4956
358.082
363.358
minSIC(k)
353.8327
358.1847
355.035
357.416
k
^
9
23
23
9Figure 1
Scatter plot and regression lines of BSE versus NYAMSE. The red points and line correspond to the first 9 observations; the green points and line correspond to the 10–35 observations.


## 6. Summary

In this paper, we proposed the Laplace linear regression model with a mean change-point and developed the EM algorithm with SIC model selection criterion to estimate the position of mean change-point. We investigated the performance of the algorithm for different simulations, finding that the algorithm behaved quite well when the errors follow the Laplace distribution. Besides that, our Laplace method is more robust than the normal method to estimate the position of mean change-point when the errors follow skew and heavy-tailed distributions, especially when the true position of 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 mean change-point. Considering the difficulty in estimating the unknown degree of freedom in the Studentt-distribution, we did not compare our Laplace model with the Studen t model in Osorio and Galea [6], where the degree of freedom was predetermined by the authors. As to the multiple mean change-points problem in Laplace linear regression, the stepwise and dichotomy methods Vostrikova [19] may work and further research is needed in the future.



---

*Source: 856350-2014-11-04.xml*