# Detection of Outliers and Patches in Bilinear Time Series Models

Authors: Ping Chen; Ling Li; Ye Liu; Jin-Guan Lin
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

We propose a Gibbs sampling algorithm to detect additive outliers and patches of outliers in bilinear time series models based on Bayesian view. We first derive the conditional posterior distributions, and then use the results of first Gibbs run to start the second adaptive Gibbs sampling. It is shown that our procedure could reduce possible effects on masking and swamping. At last, some simulations are performed to demonstrate the efficacy of detection and estimation by Monte Carlo methods.


---

## Body

## 1. Introduction

Dynamic systems or engineering time series observations [1] are often perturbed by some interrupting phenomena which generate aberrant data. They may be due to some unusual events, such as sudden disturbed factor, noise, and even recording errors. Such values are usually referred to as outliers, which may be isolated or of patch. Because outliers and patches in a time series can have adverse effects in data analysis, which maybe make the resultant inference unreliable or even invalid, it is important to detect and remove such outlier's effects. Some authors have considered the detection problem in the linear and nonlinear models. Abraham and Box [2] analyzed some outlier problems based on Bayesian methods in time series. Tsay [3] considered time series model specification in the presence of outliers. Ljung [4] presented outlier detection in time series. Justel et al. [5] studied detection of outlier patches in autoregressive time series. Battaglia and Orfei [6] considered outlier detection and estimation in nonlinear time series. Chen et al. [7] studied similar problem in ARMAX time series models. P. Chen and Y. Chen [8] gave the identification of outliers in ARMAX models via genetic algorithm.Based on some different prior distributions, a new Gibbs sampling algorithm is proposed for identifying additive isolated outliers and patches of outliers in bilinear time series. This paper considers detection of outliers and patches in bilinear time series models, which are one of the fractal time series models [9]. Wang and Wei [10] analyzed the probabilistic structure for a rather general class of bilinear models systematically. Wang [11] considered parameter estimation and subset selection for separable lower triangular bilinear models. In particular, Chen [12] proposed a procedure for detecting additive outliers in bilinear time series, which also discusses some major problems encountered in practice, such as how one can distinguish between ARMA model with outliers and a bilinear model without outliers. Some special cases of bilinear models, such as diagonal bilinear model, vector bilinear model and so on, have been studied by many authors on not only probabilistic structure but also statistical inferences; for instance, Subba Rao [13], Subba Rao and Gabr [14], Kim et al. [15], and Sesay and Subba Rao [16], among others.The format of the paper is as follows. Section2 gives the outliers model of bilinear series. Section 3 presents a method to detect outliers in bilinear model via standard Gibbs sampling. In Section 4, we propose an adaptive Gibbs sampling method to detect patches of outliers in bilinear model. Section 5 gives simulated examples and some conclusions.


## 2. Outliers Models of Bilinear Time Series

The time series{yt} is called a bilinear process, if it satisfies the model(2.1)yt=∑i=1pϕiyt-i+ɛt+∑j=1qθjɛt-j+∑i=1r∑j=1kηijyt-iɛt-j,
where  ɛt~IIDN(0,σ2), and ϕi,θj and ηij are unknown parameters. This model is called a bilinear model with order (p,q,r,k). The bilinear time series models were introduced by Granger and Anderson [17]. It is a generalization of the linear ARMA model by incorporating the product term of time series yt-i and innovation  ɛt-j. This bilinear model has attracted much attention in the recent time series literature. In terms of potential applications, bilinear models are known to be able to mode occasional outbursts in time series, which might be useful for modeling seismological data such as records for explosions and earthquakes.When additive outliers (AO) are present,yt is disturbed and unobservable. In this case, it is assumed that the observed series {zt} follows(2.2)zt=yt+βtδt,
where δt's are independent and identically distributed Bernoulli random variables with P(δt=1)=α, and βt's are random variables from normal distribution. That means that the observation zt may be AO with probability α; its magnitude is βt at time t.For simplicity, lets=max{p,r} and assume that y1,…,ys are fixed and zt=yt for t=1,…,s, that is, there exist no outliers in the first s observations. The indicator vector of outliers then becomes  δ=(δs+1,δs+2,…,δn)′ and the size vector is  β=(βs+1,βs+2,…,βn)′. Assume that a vector of residuals ɛ̂=(ɛ̂1,…,ɛ̂n-1)′ of model (2.1) is available. And let   ɛ=(ɛ1,…,ɛn-1)′,Θ=(ϕ1,…,ϕp,θ1,…,θq,η11,…,ηrk)′, z=(z1,z2,…,zn)′, and Φt-1=(yt-1,…,yt-p,ɛ̂t-1,…,ɛ̂t-q,yt-1ɛ̂t-1,…,yt-rɛ̂t-k)′, where ɛ̂j denotes the estimation of ɛj in model (2.1) without AO. Then the model (2.1) and (2.2) may be transformed as the follows(2.3)yt=Θ′Φt-1+ɛt,zt=yt+βtδt.


## 3. Outliers Detection via Standard Gibbs Sampling

In order to detect outliers in the bilinear model, it is essential to derive the conditional posterior distributions of parametersΘ,β,δ,α,σ2,δj=1, and βj.For computational reason, we give the following conditions: (1) using conjugate prior distribution for parametersΘ and σ2, which distributed as multidimensional uniform distribution on [0,1] region and inverted-Gamma distribution IG(ν/2,νλ/2), respectively. (2) Assume that the outlier indicator δt and the outlier magnitude βt are independent and distributed as Bernoulli(α) and N(0,τ2), respectively for all t. (3) The prior distribution of the contamination parameter α is Beta(γ1,γ2), and (βt)′s are i.i.d. for all t. (4) The hyperparameters in our model are λ,ν,γ1,γ2, and τ2, all of which are assumed to be known.Note that the prior probability of being contaminated by an outlier is the same for all observations, namely,P(δt=1)=α, for t=s+1,…,n. Then, under these conditions and by using the standard Bayesian method, we can obtain the following results.Theorem 3.1.
The conditional posterior distribution ofΘ is N(Θ*,V*-1),  where
(3.1)Θ*=V*-1(∑t=s+1nΦt-1ytσ2),V*=∑t=s+1nΦt-1Φt-1′σ2.Proof.
Under the conditions above, we have that(3.2)p(Θ∣z,δ,β,α,σ2)∝L(Θ,σ2,δ,β,α∣z)p(Θ)∝exp{-12σ2∑t=s+1n(yt-Θ′Φt-1)2}∝exp{-12σ2[Θ′∑t=s+1nΦt-1Φt-1′Θ-2Θ′∑t=s+1nΦt-1yt+∑t=s+1nyt2]}∝exp{-12[Θ′∑t=s+1nΦt-1Φt-1′σ2Θ-2Θ′∑t=s+1nΦt-1ytσ2+Θ*′Φt-1Φt-1′σ2Θ*]}⋮∝exp{-12(Θ-Θ*)′V*(Θ-Θ*)}~N(Θ*,V*-1).
This completes the proof of the theorem.Theorem 3.2.
(1) The conditional posterior distribution ofσ2 is IG((ν+(n-s))/2,(νλ+𝕊2)/2), where 𝕊2=∑t=s+1n(yt-Θ′Φt-1)2.
(2) The conditional distribution ofα given the sample and the other parameters is Beta(γ1+k,γ2+n-s-k), where k denotes the number of δj=1,j=1,…,n.
(3) Whenδj=0, there is no new information about the posterior distribution of βj, namely, βj distributed as N(0,τ2). When δj=1, since zt contains information of βj, we have that
(3.3)p(βj∣z,δ,β(-j),Θ,σ2,α)~N(βj*,τj*),
where β(-j) is obtained β by eliminating the element βj, and βj*=𝔸τ2/(𝔹τ2+σ2),τj*=σ2τ2/(𝔹τ2+σ2),  where 𝔸=-∑t=jTjat(1)φt-j and 𝔹=∑t=jTjφt-j2. The proof of this Theorem is very tedious and similar to that in MuCulloch and Tsay [18]; we omit it here.Theorem 3.3.
For the conditional posterior distribution ofδj=1, we have
(3.4)p(δj=1∣z,δ(-j),β,Θ,σ2,α)=[1+1-ααB10(j)]-1,
where δ(-j) is obtained δ by eliminating the element δj, Tj=min{n,j+s}, and
(3.5)B10(j)=exp{12σ2[∑t=jTjat2(1)-∑t=jTjat2(0)]};
here at(1)=(yt-Θ′Φt-1)δj=1, at(0)=(yt-Θ′Φt-1)δj=0, and at(0)=at(1)+φt-jβj. When s=0, then φi=0 for all i; when s>0, we have that
(3.6)φi={-1,i=0,ϕi,i=1,…,s,0,i≥s+1.Proof.
The conditional posterior distribution ofδj(j=s+1,…,n) is Bernoulli with probability
(3.7)p(δj=1∣z,δ(-j),β,Θ,σ2,α)=αL(Θ,σ2,δj=1,β,α∣z)αL(Θ,σ2,δj=1,β,α∣z)+(1-α)L(Θ,σ2,δj=0,β,α∣z)=αexp{-(1/2σ2)∑t=jTj(yt-Θ′Φt-1)2}δj=1αexp{-(1/2σ2)∑t=jTj(yt-Θ′Φt-1)2}δj=1+(1-α)exp{-(1/2σ2)∑t=jTj(yt-Θ'Φt-1)2}δj=0=αα+(1-α)exp{(1/2σ2)[∑t=jTj[at(1)]2-∑t=jTj[at(0)]2]}=[1+1-ααB10(j)]-1.
So the conclusion is obtained.


## 4. Detection of Outlier Patches via Adaptive Gibbs Sampling

Similar to Justel et al. [5], our procedure also consists of two Gibbs runs. In the first run, the standard Gibbs sampling based on the results of Section 3 is carried out. The results of this Gibbs run are then used to implement a second Gibbs sampling that is adaptive in treating identified outliers and in using block interpolation to reduce possible masking and swamping effects. Let Θ̂(m),σ̂(m),β̂(m), and α̂(m) be the posterior means of Θ,σ2,β, and α respectively based on the m iterations of the first Gibbs run. First, we select an appropriate critical value c1 to identify potential outliers. An observation zj is identified as an outlier if the posterior probability p̂j(m)>c1. Let {t1,…,tg} be the collection of time indexes of outliers identified by the first Gibbs run. Second, let c2(c2≤c1) be another appropriate critical value to specify the beginning and end points of a potential outlier patch. We select a window of length 2h around the identified outlier to search for the boundary points of a possible outlier patch by a forward-backward method. Exampling an identified outlier zti, for the h observations before zti if their posterior probabilities p̂j(m)>c2, then these points are regarded as possible outlier patch associated with zti. We then select the farthest point from zti as the beginning point of the outlier patch. Denote the point by zti-ki. Then we do the same for the h observations after zti and select the farthest point from zti with p̂j(m)>c2 as the end point of the outlier patch. Denote the end point by zti+vi. Combine the two blocks to form a possible outlier patch associated with zti, which is denoted by (zti-ki,…,zti+vi). Consecutive or overlapping patches should be merged to form a larger patch. Lastly, draw Gibbs samples jointly within a patch. Suppose that a patch of d outliers starting at time index j is specified. Denote the vectors of outlier indicators and magnitudes by δj,d=(δj,…,δj+d-1)′ and βj,d=(βj,…,βj+d-1)′, respectively, associated with the patch. We give the conditional posterior distributions of δj,d and βj,d in the next proposition; its derivation is similar to Theorem 1 of Justel et al. [5].Proposition 4.1.
Letz=(z1,z2,…,zn) be a vector of observations according to (2.3), with no outliers in the first s data points. Assume δj~Bernoulli(α),j=s+1,s+2,…,n, and
(4.1)p(Θ,σ2,α,β)∝(1σ2)ν/2+1αγ1-1(1-α)γ2-1exp{-12τ2∑t=s+1n(βt)2-υλ21σ2},
where the parameters γ1,γ2,υ,λ, and τ2 are known. Let et(0)=yt-∑i=1pϕiyt-i-∑j=1qθjɛ̂t-j-∑i=1r∑j=1kηijyt-iɛ̂t-j be the residual at time t when the series is adjusted for all identified outliers not in the interval [j,j+d-1] with the outliers identified in δj,d, Tj,d=min{n,j+d+s-1}, and
(4.2)et(δj,d)={et(0)+∑i=0t-jπiδt-iβt-i,t=j,…,j+d-1,et(0)+∑i=t-j-d+1t-jπiδt-iβt-i,t>j+d-1,
where π0=-1,πi=ϕi if i=1,…,s, and πi=0 if i≥s+1. Then we have the following.(a)
The conditional posterior probability of a block configurationδj,d, given the sample and the other parameters, is
(4.3)p(δj,d∣z,θδj,d)=Cαsj,d(1-α)d-sj,dexp{-12σ2∑t=jTj,det(δj,d)2},
where sj,d=∑t=jj+d-1δt,θδj,d denotes all the other parameters except δj,d, and C is a normalization constant so that the total probability of the 2d possible configurations of δj,d is one.(b)
The conditional posterior distribution ofβj,d given the sample and other parameters is Nd(βj,d*,Ωj,d-1), where(4.4)βj,d*=-1σ2Ωj,d-1∑t=jTj,det(0)Dj,dΠt-j,Ωj,d=1σ2Dj,d(∑t=jTj,dΠt-jΠt-j′)Dj,d+1τ2I,Dj,d=diag{δj,…,δj+d-1}, and Πt=(πt,…,πt-d+1)′ with π0=-1,πi=ϕi for i=1,…,s, and πi=0 for i<0 or i>s.For the second adaptive Gibbs sampling, we use the results of the first Gibbs run to start the second Gibbs sampling and to specify prior distributions of the parameters. For each outlier patch, we use the results of Proposition4.1 to draw δj,d and βj,d in the second Gibbs sampling, which is also run for m iterations.The starting values ofδt are as follows: δt(0)=1 if p̂t(m)>c2, otherwise, δt(0)=0. Then the prior distributions of βt are as follows.(a)
Ifzt is identified as an isolated outlier, then the prior distribution of βt is N(β̂t(m),τ2), where β̂t(m) is the Gibbs estimate of βt from the first Gibbs run.(b)
Ifzt belongs to an outlier patch, then the prior distribution of βt is N(β̃t(m),τ2), where β̃t(m) is the conditional posterior mean as follows:
(4.5)β̃j,d=(Dj,d∑t=jTj,dΠt-jΠt-j′Dj,d)-1(-∑t=jTj,det(0)Dj,dΠt-j).(c)
Ifzt does not belong to any outlier patch and is not an isolated outlier, then the prior distribution of βt is N(0,τ2).


## 5. Simulated Example and Conclusions

Example 5.1.
In the simulations, the designed bilinear model is BL(2,2,1,1) model:(5.1)yt=0.5yt-1+0.3yt-2+ɛt-0.34ɛt-1-0.3ɛt-2+0.4yt-1ɛt-1,zt=yt+5δt,31-4δt,32+5δt,33-4δt,34+3δt,35-6δt,50,
where δt,k is Kronecker symbol: if t=k, then δt,k=1, else δt,k=0, and {ɛt} is standard normal white noise.We produce a set of observations of the bilinear model (5.1) by simulation. In order to obtain stationary data, we produce 1000 observations of yt firstly; then we use the last 100 observations as our simulated sample, in which a patch of five consecutive additive outliers has been introduced from t=31 to t=35, a single AO has been added at t=50, and the outlier magnitudes are β31=5,β32=-4,β33=5,β34=-4,β35=3, and β50=-6, respectively. Figure 1 shows the trend of simulated series {zt}. It is obvious that the curve of observations had large volatility; it would be very difficult to distinguish between “outliers” and normal points of nonlinear model.Figure 1
The curve of the simulated observations{zt}Letγ1=5,γ2=95,ν=3,λ=σ2/2,α=0.5,τ=25σ, and c1=0.5. Here γ1=5 means that we believe the prior probability of each point is an outlier approximate to 0.05. First, we detect the outliers in {zt} using the standard Gibbs sampling. We can obtain the posterior probabilities that each observation is an outlier by using the given methods in Section 3, which are shown in Figure 2. It displayed that the posterior probabilities of being outliers only at t=31 and t=50 were larger than 0.5. Meanwhile, the outlying posterior probabilities of other observations were lower; the posterior probability of being an outlier at t=35 was even smaller than 0.2. We see that the standard Gibbs sampling failed to detect the inner and border points at t=32,33,35, that resulted in the masking effect.Figure 2
The posterior probabilities that each point is an outlier via standard Gibbs sampling.Second, in order to avoid the presence of masking and swamping problem, we use the method given in Section4, and take 900 iterations by the adaptive Gibbs algorithm. Figure 3 gives the plot of the estimated posterior probability that each point is seen as an outlier via adaptive Gibbs sampling in simulation, where the window width of search was 6. It is obvious that the posterior probabilities of being an outlier obtained by adaptive Gibbs algorithm for data points at t=31,…,35 and t=50 are larger than 0.5, which meant that these points were possible outliers. Meanwhile, the outlying posterior probabilities of other observations are very small. Actually, if we select the critical values c1=0.5 and c2=0.15, then we could identify the patch. On the other hand, many normal points like outliers in Figure 1 were not misspecified as outliers because the outlying posterior probabilities of these points were smaller than 0.2, which shows that the adaptive Gibbs sampling was more effective than the standard Gibbs sampling in mining the additive outlier patch for bilinear time series. Third, the posterior means of the sizes of these outliers are β̂31=4.6714, β̂32=-3.4534, β̂33=4.3176, β̂34=-6.9497, β̂35=5.1142, and β̂50=-2.2523, respectively.Figure 3
The posterior probability that each point is an outlier via adaptive Gibbs algorithm.Furthermore, we did 200 replications for the above simulation. Table1 gives the average values of the posterior means of β31,…,β35 and β50 and their standard deviations. In these replications, the isolated outliers were identified by standard and adaptive Gibbs sampling procedures. Through computing, we find that the standard Gibbs sampling failed to detect all the outlier patches in 200 replications, but in 179 replications, the adaptive Gibbs sampling successfully could specify all outliers and patches. On the other hand, we note that the last three parameters are actually not estimated well by the adaptive Gibbs algorithm in many replications. Investigating its reason, we think it may be the influence of nonlinearity or the frontal outliers of patch. Such situation deserves careful investigation in the future.Table 1
The simulated results of the model BL(2,2,1,1) via adaptive Gibbs algorithm.
β31β32β33β34β35β50True value5-45-43-6Posterior mean4.7139-3.65294.1193-6.94874.7747-2.6025Standard deviation0.62110.65550.60440.63740.67110.5813Since the bilinear series will by itself produce what appears to be outliers, the detection of the locations of outliers and patches is more difficult, and the estimation of the magnitude of outliers and patches is more complex than autoregressive series. By a number of simulations of detecting the outlier patches in bilinear model by adaptive Gibbs sampling, we discovered that the critical valuec2 should be selected smaller than that in ARMA model. It may be that the variation of bilinear series is larger than that of ARMA series. On the other hand, in the process of running, the Gibbs sampler was repeated several times with different hyper-parameters and different numbers of iterations to reanalyze the data. The results show that the locations of possible outliers and patch are stable, even though the estimated outlying probabilities may vary slightly between the Gibbs samples. Some other case studies also show that our procedure could reduce possible masking and swamping effects, which is an improvement and extension on bilinear time series model over the existing outliers detection methods.



---

*Source: 580583-2010-04-20.xml*