Detection of Outliers and Patches in Bilinear Time Series Models

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.

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. Section 2 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 {𝑦_𝑡} is called a bilinear process, if it satisfies the model

𝑦_𝑡=^𝑝_𝑖=1𝜙_𝑖𝑦_𝑡−𝑖+𝜀_𝑡+^𝑞_𝑗=1𝜃_𝑗𝜀_𝑡−𝑗+^𝑟𝑘_𝑖=1𝑗=1𝜂_𝑖𝑗𝑦_𝑡−𝑖𝜀_𝑡−𝑗, (2.1) where 𝜀_𝑡∼IIDN(0,𝜎²), and 𝜙_𝑖,𝜃_𝑗 and 𝜂_𝑖𝑗 are unknown parameters. This model is called a bilinear model with order (𝑝,𝑞,𝑟,𝑘). 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 𝑦_𝑡−𝑖 and innovation 𝜀_𝑡−𝑗. 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, 𝑦_𝑡 is disturbed and unobservable. In this case, it is assumed that the observed series {𝑧_𝑡} follows

𝑧_𝑡=𝑦_𝑡+𝛽_𝑡𝛿_𝑡, (2.2) where 𝛿_𝑡's are independent and identically distributed Bernoulli random variables with 𝑃(𝛿_𝑡=1)=𝛼, and 𝛽_𝑡's are random variables from normal distribution. That means that the observation 𝑧_𝑡 may be AO with probability 𝛼; its magnitude is 𝛽_𝑡 at time 𝑡.

For simplicity, let 𝑠=max{𝑝,𝑟} and assume that 𝑦₁,…,𝑦_𝑠 are fixed and 𝑧_𝑡=𝑦_𝑡 for 𝑡=1,…,𝑠, that is, there exist no outliers in the first 𝑠 observations. The indicator vector of outliers then becomes 𝜹=(𝛿_𝑠+1,𝛿_𝑠+2,…,𝛿_𝑛)^′ and the size vector is 𝜷=(𝛽_𝑠+1,𝛽_𝑠+2,…,𝛽_𝑛)^′. Assume that a vector of residuals 𝜺=(𝜀₁̂,…,𝜀̂_𝑛−1)^′ of model (2.1) is available. And let 𝜺=(𝜀₁,…,𝜀_𝑛−1)^′,Θ=(𝜙₁,…,𝜙_𝑝,𝜃₁,…,𝜃_𝑞,𝜂₁₁,…,𝜂_𝑟𝑘)^′, 𝐳=(𝑧₁,𝑧₂,…,𝑧_𝑛)^′, and Φ_𝑡−1=(𝑦_𝑡−1,…,𝑦_𝑡−𝑝,𝜀̂_𝑡−1,…,𝜀̂_𝑡−𝑞,𝑦_𝑡−1𝜀̂_𝑡−1,…,𝑦_𝑡−𝑟𝜀̂_𝑡−𝑘)^′, where 𝜀_𝑗̂ denotes the estimation of 𝜀_𝑗 in model (2.1) without AO. Then the model (2.1) and (2.2) may be transformed as the follows

𝑦_𝑡=Θ^′Φ_𝑡−1+𝜀_𝑡,𝑧_𝑡=𝑦_𝑡+𝛽_𝑡𝛿_𝑡. (2.3)

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 Θ,𝜷,𝜹,𝛼,𝜎²,𝛿_𝑗=1, and 𝛽_𝑗.

For computational reason, we give the following conditions: (1) using conjugate prior distribution for parameters Θ and 𝜎², which distributed as multidimensional uniform distribution on [0,1] region and inverted-Gamma distribution 𝐼𝐺(𝜈/2,𝜈𝜆/2), respectively. (2) Assume that the outlier indicator 𝛿_𝑡 and the outlier magnitude 𝛽_𝑡 are independent and distributed as 𝐵𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝛼) and 𝑁(0,𝜏²), respectively for all 𝑡. (3) The prior distribution of the contamination parameter 𝛼 is 𝐵𝑒𝑡𝑎(𝛾₁,𝛾₂), and (𝛽_𝑡)^′𝑠 are i.i.d. for all 𝑡. (4) The hyperparameters in our model are 𝜆,𝜈,𝛾₁,𝛾₂, and 𝜏², 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, 𝑃(𝛿_𝑡=1)=𝛼, for 𝑡=𝑠+1,…,𝑛. 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 𝑁(Θ^∗,𝑉_∗⁻¹), where Θ^∗=𝑉_∗⁻¹(∑^𝑛_𝑡=𝑠+1Φ_𝑡−1𝑦_𝑡𝜎²),𝑉_∗=∑^𝑛_𝑡=𝑠+1Φ_𝑡−1Φ^′_𝑡−1𝜎². (3.1)

Proof. Under the conditions above, we have that 𝑝(Θ∣𝐳,𝜹,𝜷,𝛼,𝜎²)∝𝐿(Θ,𝜎²1,𝜹,𝜷,𝛼∣𝐳)𝑝(Θ)∝exp{−2𝜎^2𝑛_𝑡=𝑠+1(𝑦_𝑡−Θ^′Φ_𝑡−1)²1}∝exp{−2𝜎²[Θ^′𝑛_𝑡=𝑠+1Φ_𝑡−1Φ^′_𝑡−1Θ−2Θ^′𝑛_𝑡=𝑠+1Φ_𝑡−1𝑦_𝑡+^𝑛_𝑡=𝑠+1𝑦²_𝑡1]}∝exp{−2[Θ^′𝑛_𝑡=𝑠+1Φ_𝑡−1Φ^′_𝑡−1𝜎²Θ−2Θ^′∑^𝑛_𝑡=𝑠+1Φ_𝑡−1𝑦_𝑡𝜎²+Θ^∗′Φ_𝑡−1Φ^′_𝑡−1𝜎²Θ^∗1]}⋮∝exp{−2(Θ−Θ^∗)^′𝑉_∗(Θ−Θ^∗)}∼𝑁(Θ^∗,𝑉_∗⁻¹). (3.2)(1)(2)(3)(4)(5)(6)(7) This completes the proof of the theorem.

Theorem 3.2. (1) The conditional posterior distribution of 𝜎² is 𝐼𝐺((𝜈+(𝑛−𝑠))/2,(𝜈𝜆+𝕊²)/2), where 𝕊²=∑^𝑛_𝑡=𝑠+1(𝑦_𝑡−Θ^′Φ_𝑡−1)².
(2) The conditional distribution of 𝛼 given the sample and the other parameters is 𝐵𝑒𝑡𝑎(𝛾₁+𝑘,𝛾₂+𝑛−𝑠−𝑘), where 𝑘 denotes the number of 𝛿_𝑗=1,𝑗=1,…,𝑛.
(3) When 𝛿_𝑗=0, there is no new information about the posterior distribution of 𝛽_𝑗, namely, 𝛽_𝑗 distributed as 𝑁(0,𝜏²). When 𝛿_𝑗=1, since 𝑧_𝑡 contains information of 𝛽_𝑗, we have that
𝑝(𝛽_𝑗∣𝐳,𝜹,𝛽_(−𝑗),Θ,𝜎²,𝛼)∼𝑁(𝛽^∗_𝑗,𝜏^∗_𝑗), (3.3) where 𝛽_(−𝑗) is obtained 𝜷 by eliminating the element 𝛽_𝑗, and 𝛽^∗_𝑗=𝔸𝜏²/(𝔹𝜏²+𝜎²),𝜏^∗_𝑗=𝜎²𝜏²/(𝔹𝜏²+𝜎²), where 𝔸=−∑^𝑇𝑗_𝑡=𝑗𝑎_𝑡(1)𝜑_𝑡−𝑗 and 𝔹=∑^𝑇𝑗_𝑡=𝑗𝜑²_𝑡−𝑗. 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 𝛿_𝑗=1, we have 𝑝(𝛿_𝑗=1∣𝐳,𝛿_(−𝑗),𝜷,Θ,𝜎²,𝛼)=[1+1−𝛼𝛼𝐵₁₀(𝑗)]⁻¹, (3.4) where 𝛿_(−𝑗) is obtained 𝜹 by eliminating the element 𝛿_𝑗, 𝑇_𝑗=min{𝑛,𝑗+𝑠}, and 𝐵₁₀1(𝑗)=exp{2𝜎²[^𝑇𝑗_𝑡=𝑗𝑎²_𝑡(1)−^𝑇𝑗_𝑡=𝑗𝑎²_𝑡(0)]}; (3.5) here 𝑎_𝑡(1)=(𝑦_𝑡−Θ^′Φ_𝑡−1)_𝛿𝑗=1, 𝑎_𝑡(0)=(𝑦_𝑡−Θ^′Φ_𝑡−1)_𝛿𝑗=0, and 𝑎_𝑡(0)=𝑎_𝑡(1)+𝜑_𝑡−𝑗𝛽_𝑗. When 𝑠=0, then 𝜑_𝑖=0 for all 𝑖; when 𝑠>0, we have that 𝜑_𝑖=𝜙−1,𝑖=0,_𝑖,𝑖=1,…,𝑠,0,𝑖≥𝑠+1. (3.6)

Proof. The conditional posterior distribution of 𝛿_𝑗(𝑗=𝑠+1,…,𝑛) is 𝐵𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖 with probability 𝑝(𝛿_𝑗=1∣𝐳,𝛿_(−𝑗),𝜷,Θ,𝜎²=,𝛼)𝛼𝐿(Θ,𝜎²,𝛿_𝑗=1,𝜷,𝛼∣𝐳)𝛼𝐿(Θ,𝜎²,𝛿_𝑗=1,𝜷,𝛼∣𝐳)+(1−𝛼)𝐿(Θ,𝜎²,𝛿_𝑗==0,𝜷,𝛼∣𝐳)𝛼exp{−(1/2𝜎²)∑^𝑇𝑗_𝑡=𝑗(𝑦_𝑡−Θ^′Φ_𝑡−1)²}_𝛿𝑗=1𝛼exp{−(1/2𝜎²)∑^𝑇𝑗_𝑡=𝑗(𝑦_𝑡−Θ^′Φ_𝑡−1)²}_𝛿𝑗=1+(1−𝛼)exp{−(1/2𝜎²)∑^𝑇𝑗_𝑡=𝑗(𝑦_𝑡−Θ^′Φ_𝑡−1)²}_𝛿𝑗=0=𝛼𝛼+(1−𝛼)exp{(1/2𝜎²)[∑^𝑇𝑗_𝑡=𝑗[𝑎_𝑡(1)]²−∑^𝑇𝑗_𝑡=𝑗[𝑎_𝑡(0)]²]}=[1+1−𝛼𝛼𝐵₁₀(𝑗)]⁻¹. (3.7) 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 Θ^(𝑚),𝜎^(𝑚),𝜷^(𝑚), and 𝛼^(𝑚) be the posterior means of Θ,𝜎²,𝜷, and 𝛼 respectively based on the 𝑚 iterations of the first Gibbs run. First, we select an appropriate critical value 𝑐₁ to identify potential outliers. An observation 𝑧_𝑗 is identified as an outlier if the posterior probability 𝑝_𝑗^(𝑚)>𝑐₁. Let {𝑡₁,…,𝑡_𝑔} be the collection of time indexes of outliers identified by the first Gibbs run. Second, let 𝑐₂(𝑐₂≤𝑐₁) be another appropriate critical value to specify the beginning and end points of a potential outlier patch. We select a window of length 2ℎ around the identified outlier to search for the boundary points of a possible outlier patch by a forward-backward method. Exampling an identified outlier 𝑧_𝑡𝑖, for the ℎ observations before 𝑧_𝑡𝑖 if their posterior probabilities 𝑝_𝑗^(𝑚)>𝑐₂, then these points are regarded as possible outlier patch associated with 𝑧_𝑡𝑖. We then select the farthest point from 𝑧_𝑡𝑖 as the beginning point of the outlier patch. Denote the point by 𝑧_{𝑡𝑖−𝑘𝑖}. Then we do the same for the ℎ observations after 𝑧_𝑡𝑖 and select the farthest point from 𝑧_𝑡𝑖 with 𝑝_𝑗^(𝑚)>𝑐₂ as the end point of the outlier patch. Denote the end point by 𝑧_{𝑡𝑖+𝑣𝑖}. Combine the two blocks to form a possible outlier patch associated with 𝑧_𝑡𝑖, which is denoted by (𝑧_{𝑡𝑖−𝑘𝑖},…,𝑧_{𝑡𝑖+𝑣𝑖}). 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 𝑑 outliers starting at time index 𝑗 is specified. Denote the vectors of outlier indicators and magnitudes by 𝜹_𝑗,𝑑=(𝛿_𝑗,…,𝛿_{𝑗+𝑑−1})^′ and 𝜷_𝑗,𝑑=(𝛽_𝑗,…,𝛽_{𝑗+𝑑−1})^′, respectively, associated with the patch. We give the conditional posterior distributions of 𝜹_𝑗,𝑑 and 𝜷_𝑗,𝑑 in the next proposition; its derivation is similar to Theorem 1 of Justel et al. [5].

Proposition 4.1. Let z =(𝑧₁,𝑧₂,…,𝑧_𝑛) be a vector of observations according to (2.3), with no outliers in the first 𝑠 data points. Assume 𝛿_𝑗∼𝐵𝑒𝑟𝑛𝑜𝑢𝑙𝑙𝑖(𝛼),𝑗=𝑠+1,𝑠+2,…,𝑛, and 𝑝(Θ,𝜎²1,𝛼,𝜷)∝(𝜎²)^𝜈/2+1𝛼^𝛾1−1(1−𝛼)^𝛾2−11exp{−2𝜏^2𝑛_𝑡=𝑠+1(𝛽_𝑡)²−𝜐𝜆21𝜎²}, (4.1) where the parameters 𝛾₁,𝛾₂,𝜐,𝜆, and 𝜏² are known. Let 𝑒_𝑡(0)=𝑦_𝑡−∑^𝑝_𝑖=1𝜙_𝑖𝑦_𝑡−𝑖−∑^𝑞_𝑗=1𝜃_𝑗𝜀̂_𝑡−𝑗−∑^𝑟_𝑖=1∑^𝑘_𝑗=1𝜂_𝑖𝑗𝑦_𝑡−𝑖𝜀̂_𝑡−𝑗 be the residual at time 𝑡 when the series is adjusted for all identified outliers not in the interval [𝑗,𝑗+𝑑−1] with the outliers identified in 𝜹_𝑗,𝑑, 𝑇_𝑗,𝑑=min{𝑛,𝑗+𝑑+𝑠−1}, and 𝑒_𝑡(𝜹_𝑗,𝑑𝑒)=^𝑡(0)+^𝑡−𝑗_𝑖=0𝜋_𝑖𝛿_𝑡−𝑖𝛽_𝑡−𝑖𝑒,𝑡=𝑗,…,𝑗+𝑑−1,_𝑡(0)+^𝑡−𝑗_{𝑖=𝑡−𝑗−𝑑+1}𝜋_𝑖𝛿_𝑡−𝑖𝛽_𝑡−𝑖,𝑡>𝑗+𝑑−1, (4.2) where 𝜋₀=−1,𝜋_𝑖=𝜙_𝑖 if 𝑖=1,…,𝑠, and 𝜋_𝑖=0 if 𝑖≥𝑠+1. Then we have the following.(a)The conditional posterior probability of a block configuration 𝜹_𝑗,𝑑, given the sample and the other parameters, is 𝑝(𝜹_𝑗,𝑑∣𝑧,𝜃_{𝛿𝑗,𝑑})=𝐶𝛼^{𝑠𝑗,𝑑}(1−𝛼)^{𝑑−𝑠𝑗,𝑑}1exp{−2𝜎^{2𝑇𝑗,𝑑}_𝑡=𝑗𝑒_𝑡(𝜹_𝑗,𝑑)²}, (4.3) where 𝑠_𝑗,𝑑=∑^{𝑗+𝑑−1}_𝑡=𝑗𝛿_𝑡,𝜃_{𝛿𝑗,𝑑} denotes all the other parameters except 𝜹_𝑗,𝑑, and 𝐶 is a normalization constant so that the total probability of the 2^𝑑 possible configurations of 𝜹_𝑗,𝑑 is one.(b)The conditional posterior distribution of 𝜷_𝑗,𝑑 given the sample and other parameters is 𝑁_𝑑(𝜷^∗_𝑗,𝑑,Ω⁻¹_𝑗,𝑑), where 𝜷^∗_𝑗,𝑑1=−𝜎²Ω^𝑇−1_𝑗,𝑑^𝑗,𝑑_𝑡=𝑗𝑒_𝑡(0)𝐷_𝑗,𝑑Π_𝑡−𝑗,Ω_𝑗,𝑑=1𝜎²𝐷_𝑗,𝑑(^{𝑇𝑗,𝑑}_𝑡=𝑗Π_𝑡−𝑗Π^′_𝑡−𝑗)𝐷_𝑗,𝑑+1𝜏²𝐼, (4.4)𝐷_𝑗,𝑑=diag{𝛿_𝑗,…,𝛿_{𝑗+𝑑−1}}, and Π_𝑡=(𝜋_𝑡,…,𝜋_{𝑡−𝑑+1})^′ with 𝜋₀=−1,𝜋_𝑖=𝜙_𝑖 for 𝑖=1,…,𝑠, and 𝜋_𝑖=0 for 𝑖<0 or 𝑖>𝑠.

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 Proposition 4.1 to draw 𝜹_𝑗,𝑑 and 𝜷_𝑗,𝑑 in the second Gibbs sampling, which is also run for 𝑚 iterations.

The starting values of 𝛿_𝑡 are as follows: 𝛿_𝑡⁽⁰⁾=1 if 𝑝_𝑡^(𝑚)>𝑐₂, otherwise, 𝛿_𝑡⁽⁰⁾=0. Then the prior distributions of 𝛽_𝑡 are as follows.

(a)If 𝑧_𝑡 is identified as an isolated outlier, then the prior distribution of 𝛽_𝑡 is 𝑁(𝛽_𝑡^(𝑚),𝜏²), where 𝛽_𝑡^(𝑚) is the Gibbs estimate of 𝛽_𝑡 from the first Gibbs run.(b)If 𝑧_𝑡 belongs to an outlier patch, then the prior distribution of 𝛽_𝑡 is 𝑁(𝛽_𝑡^(𝑚),𝜏²), where 𝛽_𝑡^(𝑚) is the conditional posterior mean as follows: 𝜷_𝑗,𝑑=(𝐷^𝑇_𝑗,𝑑^𝑗,𝑑_𝑡=𝑗Π_𝑡−𝑗Π^′_𝑡−𝑗𝐷_𝑗,𝑑)⁻¹(−^{𝑇𝑗,𝑑}_𝑡=𝑗𝑒_𝑡(0)𝐷_𝑗,𝑑Π_𝑡−𝑗). (4.5)(c)If 𝑧_𝑡 does not belong to any outlier patch and is not an isolated outlier, then the prior distribution of 𝛽_𝑡 is 𝑁(0,𝜏²).

5. Simulated Example and Conclusions

Example 5.1. In the simulations, the designed bilinear model is BL(2,2,1,1) model: 𝑦_𝑡=0.5𝑦_𝑡−1+0.3𝑦_𝑡−2+𝜀_𝑡−0.34𝜀_𝑡−1−0.3𝜀_𝑡−2+0.4𝑦_𝑡−1𝜀_𝑡−1,𝑧_𝑡=𝑦_𝑡+5𝛿_𝑡,31−4𝛿_𝑡,32+5𝛿_𝑡,33−4𝛿_𝑡,34+3𝛿_𝑡,35−6𝛿_𝑡,50, (5.1) where 𝛿_𝑡,𝑘 is Kronecker symbol: if 𝑡=𝑘, then 𝛿_𝑡,𝑘=1, else 𝛿_𝑡,𝑘=0, and {𝜀_𝑡} 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 𝑦_𝑡 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 𝑡=31 to 𝑡=35, a single AO has been added at 𝑡=50, and the outlier magnitudes are 𝛽₃₁=5,𝛽₃₂=−4,𝛽₃₃=5,𝛽₃₄=−4,𝛽₃₅=3, and 𝛽₅₀=−6, respectively. Figure 1 shows the trend of simulated series {𝑧_𝑡}. 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 $\{z_t\}$

Let 𝛾₁=5,𝛾₂=95,𝜈=3,𝜆=𝜎²/2,𝛼=0.5,𝜏=2√5𝜎, and 𝑐₁=0.5. Here 𝛾₁=5 means that we believe the prior probability of each point is an outlier approximate to 0.05. First, we detect the outliers in {𝑧_𝑡} 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 𝑡=31 and 𝑡=50 were larger than 0.5. Meanwhile, the outlying posterior probabilities of other observations were lower; the posterior probability of being an outlier at 𝑡=35 was even smaller than 0.2. We see that the standard Gibbs sampling failed to detect the inner and border points at 𝑡=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 Section 4, 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 𝑡=31,…,35 and 𝑡=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 𝑐₁=0.5 and 𝑐₂=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 𝛽₃₁=4.6714, 𝛽₃₂=−3.4534, 𝛽₃₃=4.3176, 𝛽₃₄=−6.9497, 𝛽₃₅=5.1142, and 𝛽₅₀=−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. Table 1 gives the average values of the posterior means of 𝛽₃₁,…,𝛽₃₅ and 𝛽₅₀ 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.

Since 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 value 𝑐₂ 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.

Acknowledgments

The authors thank the editor and three referees for their careful and constructive comments that led to substantial improvements of the paper. The project is supported by NSFCJS BK2008284.