Advances in Systems Science and Applications (2012) Vol.12 No.4 399-405

The Identification of Outliers in ARMAX Models via
Genetic Algorithm
Ping Chen1 and Ying Chen2

1Department of Mathematics, Southeast University, Nanjing 210096, China
2Department of Forensic Science, Jiangsu Police Institute , Nanjing, 210012,

China

Abstract

This paper proposes a procedure to identify additive and innovational outlier-
s by genetic algorithm in autoregressive moving average with exogenous vari-
able(ARMAX) time series models. We use some methods to delete the influence
of input process in ARMAX model and then detect outliers in time series based
on the previous work, which is an improvement and extension of the detection
method on ARMA models. Empirical and simulation studies show that the pro-
posed procedure is effective.
Keywords dynamic systems, innovational outliers, ARMAX model, genetic algo-
rithm

1 Introduction

Outliers in dynamic systems or engineering time series can have adverse effects on
model identification and parameter estimation. Several procedures are available
in literature to handle outliers in a time series. However, the case of multiple
additive outliers and innovational outliers is very difficult to study because of the
great number of alternatives and of the masking and swamping effects. Compared
to other search algorithms, genetic algorithms allow many candidate solutions to
be considered simultaneously at each step. Baragona et al.[1] showed how to
use genetic algorithms for outlier detection and classification in ARMA series.
Peña and Sánchez[2] presented a new procedure for multifold predictive validation
in ARMAX models. Also, Chen et al.[3] and Chen et al.[4] developed some
methods for detecting outliers, change point and outlier patches in bilinear time
series models. On the other hand, Huang et al.[5] discussed the improved genetic
algorithm for vehicle routing problem with time windows. In this paper, a genetic
algorithm is proposed to identify additive and innovational outliers in ARMAX
series. We are using the standard genetic algorithm with complete replacement
of the past population and elitist strategy. The relationship between inverse
correlations and outliers is helpful to simplify the fitness function, which may be
quickly computed by Trench’s algorithm. For the case of large series, it is better
to divide the series into several parts so that the two neighbor subseries share
a length of same data, and then to detect the subseries respectively. Thus, the
problem of large population in genetic algorithms has been avoided, as well as



400 Ping Chen:The Identification of Outliers in ARMAX Models via Genetic Algorithm

the loss of outliers around the cut point. At last, simulation studies are carried
out, which show promising results.

2 Genetic Algorithm and Outliers Models in ARMAX Series

The genetic algorithm(GA) is known to be able to provide us with a powerful
optimization tool when the solution space happens to be both discrete and large,
and the objective function does not fulfill the usual regularity requirements.The
key feature of a GA is the manipulation of a population whose individuals are
characterized by possessing a chromosome. This latter can be coded as a string
of characters of given length. Each string represents a feasible solution to the
optimization problem. The link between the GA and the problem at hand is
provided by the fitness function (FF). The FF establishes a mapping from the
chromosomes to some set of real numbers. The greater the FF is, the better the
adaptation of the individual. The procedure is iterative. It makes use of three
evolutionary operators: reproduction, crossover and mutation.

An ARMA model with input process is called ARMAX model, which is defined
as

Zt =
d∑

i=1

υi(B)Xi,t + nt,

where υi(B) = (δ−1
i (B) · ωi(B))Bki is the transfer function of ith input process,

nt = (θ(B)/ϕ(B))εt is noise process. {Zt} is called response process. And Xi,t

denotes the ith input process or the difference of ith input process at time t, ki
presents the influence’s time delay of ith input process, εt is normal white nose
process.

θ(B) = 1 + θ1B + · · ·+ θqB
q, ϕ(B) = 1− ϕ1B − · · · − ϕpB

p

where B is the backshift operator.
When υi(B) = 0, i = 1, · · · , d, it is ARMA model, when some υi(B) is nonzero

constant, i = 1, · · · , d, it is regression model with ARMA error.
(1) Additive outliers(AO) model:
Suppose that only the jth point zj be AO, whose influence magnitude is wtj ,

then we have

Zt =
d∑

i=1

υi(B)Xi,t + wtjδt,tj +
θ(B)

ϕ(B)
εt

where δt,tj is Kronecker symbol: If t = tj, then δt,tj = 1, else δt,tj = 0.
(2) Innovational outliers(IO) model
Suppose that only the jth point zj be IO, whose influence magnitude is wtj ,



Advances in Systems Science and Applications (2012) Vol.12 No.4 401

then we have

Zt =
d∑

i=1

υi(B)Xi,t +
θ(B)

ϕ(B)
(εt + wtjδt,tj)

=

d∑
i=1

υi(B)Xi,t + wtj
θ(B)

ϕ(B)
δt,tj +

θ(B)

ϕ(B)
εt.

3 Identification of ARMAX Models

Suppose that the ARMAX model of only one input process is as follows:

Zt = δ−1(B)ω(B)Xt−b + nt = υ(B)Xt + nt, (1)

Where

δ(B) = 1− δ1(B)− · · · − δr1B
r1, ω(B) = ω0 − ω1(B)− · · · − ωr2B

r2

and
υ(B) = δ−1(B)ω(B)Bb

Suppose the input process Xt is stationary and is able to be represented by
some member of the general linear class of autoregressive-moving average mod-
els. Given a set of data, similar to Box et al.[6], then we can carry out our
usual identification and estimation methods to obtain a model for the Xt process
ϕ(B)θ−1(B)Xt = αt which, to a close approximation, transforms the correlated
input series Xt to the uncorrelated white nose series αt. At the same time, we can
obtain an estimate s2α of σ2α from the sum of squares of the α̂′s. If we now apply
this same transformation to Zt to obtain βt = ϕ(B)θ−1(B)Zt, then the model(1)
may be written βt = υ(B)αt + εt, multiplying αt−k on both sides and taking
expectations, we obtain γαβ(k) = υkσ

2
α, where γαβ(k) = E[αt−kβt]is the cross

covariance at lag k between α and β. Thus υk = [ραβ(k)σβ]/[σα], k = 0, 1, 2 · · · .
Hence, after “prewhitening” the input, the cross correlation function between

the prewhitened input and correspondingly transformed output is directly pro-
portional to the response function. In practice, we do not know the theoretical
function ραβ(k), so we must substitute estimates in υk to give

υ̂k = [rαβ(k)sβ]/[sα], k = 0, 1, 2...

where
rαβ(k) = cαβ(k)/[sα/sβ],

cαβ(k) =
1

n

n−k∑
n=1

(αt − ᾱ)(βt+k − β̄),

sα =
√
cαα(0),

sβ =
√
cββ(0), k = 0, 1, 2, ...



402 Ping Chen:The Identification of Outliers in ARMAX Models via Genetic Algorithm

The preliminary estimates υ̂k can provide a rough basis for selecting suitable
transfer function model. First, we may use the estimates υ̂k so obtained to make
guesses of the order r1 and r2 of δ(B) and ω(B), and of the delay parameter
b. Second, we do not consider the noise nt now, substituting Zt = υ̂(B)Xt in
the equation δ(B)Zt = ω(B)BbXt, based on equating coefficients of B, to obtain
initial estimates of the parameters δ(B) and ω(B).
4 The Identification of Outliers Via Genetic Algorithm

We let {Yt, t = 0, 1, 2, · · · } be a zero mean and stationary time series: Yt =
∞∑
j=0

ψjαt−j , where {αt} is Gaussian zero mean white noise and V ar(αt) = σ2,

{ψj , j = 0, 1, 2, · · ·} form a absolutely summable sequence. Let γih denote the
inverse autocovariance function of the process, for integer h. Also let ρih =
γih/γi0 denote the inverse autocorrelations. We have

γik + ψ1γik−1 + ψ2γik−2 + · · · = 0, k > 0 (2)

When outliers are present, {Yt, t = 0, 1, 2, · · ·} is unobservable. Instead the
time series {Zt, t = 0, 1, 2, · · ·} is observed which follows the model: Zt = Yt + dt,
where dt is a deterministic perturbation. Let ψj = 0 if j < 0, then we have
dt = ψt−t0w0, if time is IO; or dt = w0δt,t0 , if time t0 is AO , where w0 denotes the
outlier’s magnitude at t = t0. In the present setting, a chromosome ξ is a string of
characters of assigned length n that can be evaluated in terms of the FF, where n
is the number of observations of the time series, as each locus ξj is corresponding
to an observation zj where an outlier may occur. So, ξ = (ξ1, ξ2, · · · , ξn). Then a
gene ξj = 0, if the locus is an outlier-free time point, ξj = 1 if the observation at
this time point is an additive outlier, and ξj = 2 if it is an innovational outlier.

If k outliers are located at t1, t2, · · · , tk and denoting by Z = (z1, z2, · · · , zn)′
the observed time series and by Y = (y1, y2, · · · , yn)′ the unobserved realization
of Yt, then Z = ΨW + Y , where W = (w1, w2, · · · , wk)

′ is the vector of the
outliers’ magnitudes at t1, t2, · · · , tk and Ψ is the matrix of the n × k elements
Ψjh defined as follows: If th is the timing of an IO, then Ψjh = ψj−th if j > th and
0, otherwise. If an AO is occurring in th, then Ψjh = 1 if j = th and 0, otherwise.
The idea is to seek for the matrix Ψ which maximizes the likelihood function.
When n is large, Γ−1 may be replaced by the matrix of inverse autocovariances
Γi . Then we have the likelihood function of the time series Z

L = p(Z | ξ,W ) = 2π−n/2(detΓi)1/2 exp{−1

2
(Z −ΨW )′Γi(Z −ΨW )} (3)

The joint maximum likelihood estimate W ∗ of W , given the pattern of the
outlying observations and Γi, is

W ∗ = (Ψ′ΓiΨ)−1Ψ′ΓiZ and Y ∗ = Z −ΨW ∗ (4)



Advances in Systems Science and Applications (2012) Vol.12 No.4 403

In practice, however, Γi is seldom known, so we have to estimate it from the
data. Similar to Baragona et al.[1], using the interpolator estimate of Yt. we may
obtain inverse autocorrelations interpolation estimates ρ̂ij , j = 1, · · · ,m and the
estimate γ̂i0 of inverse variance and the estimates of inverse autocovariance γ̂ij =
γ̂i0 · ρ̂ij , j = 1, · · · ,m. Once the inverse autocovariances have been estimated,
the ψj , j = 1, · · · , q may be computed using Equation (2). Nevertheless, these
estimates are biased because of the presence of the outliers. So we have to resort to
an iterative scheme for (4), which is repeated until convergences. If convergence
is attained, and since the likelihood is maximized with respect to the inverse
autocorrelations, the argument of the exponential in likelihood’s Formula (3) is
simply −n/2 , which is shown by the following proposition(the proof is omited).

Proposition 1 Using the former representation and estimations, then we have

−1

2
(Z −ΨW )′Γ(Z −ΨW ) = −1

2
n

Therefore, we have

logL = −n
2
log(2π) +

1

2
log(detΓi)− 1

2
n

Thus, the likelihood function basically depends on detΓi. Because that the
fitness function(FF) has to be positive, we let F (ξ) = a · blog(detΓi)−ck, where b
is a real constant such that b > 1, and the constant terms in the exponent were
dropped. We use a = 10, b = 1.0001 and c = 14 in this paper. We adopted (1)
Population size s = n + 1; (2) Probability of crossover pc = 0.8; (3) Probability
of mutation pm = 0.05; (4) The maximum number of outliers within a chromo-
some g = 10; (5) The number of iterations of the series’ adjustment/parameters’
estimation procedure was 3;(6) we perform as many iterations of the genetic al-
gorithm as possible within a reasonable time period.

5 Simulation Studies

Example A In the following example, we consider the model
(1− 0.8B + 0.3B2)xt = εt

zt = (1− 0.5B)xt − 4δt,30 + 5δt,31 − 5δt,80 − 4δt,90+

6× 1− 0.36B + 0.85B2

1− 0.6B
δt,40 +

1− 0.36B + 0.85B2

1− 0.6B
et

where {εt} and {et} are all normal white noise, their means are zero and variance
σ2 = 1.

We create 101 observations x0, x1, · · · , x100 of xt and 100 observations z1,...,z100
of zt by simulation. Obviously, it is AO at t = 30, 31, 80, 90 singly and IO at



404 Ping Chen:The Identification of Outliers in ARMAX Models via Genetic Algorithm

t = 40 , and outlier magnitudes are w30 = −4, w31 = 5, w80 = −5, w90 = −4 and
w40 = 6, respectively.

Applying our method to the above data and prewhitening the input series.
Making {xt} follows an ARMA model:

(1− 0.88088B + 0.32738B2)xt = εt.

Then we take the same manipulation to prewhiten {zt}. By analyzing filtered
cross correlation coefficient of {zt} and {xt}, we obtain the transfer function
1.03416−0.73967B for {xt}. Delete the influence of input process {xt} in response
process {zt}, and let z∗t = zt − (1.03416 − 0.73967B)xt. We have that {z∗t } is
an ARMA series include outliers. We detect the outliers in {z∗t } by applying the
above method. Let m = 2, q = 9, g = 10, s = 101. Because the 30th point is
very close to 31th point and they influence each other, it is difficult to identify
the outliers. In this case, one needs larger number of iterations. We take 1000
iterations by standard genetic algorithm. and obtain the best individual at 612th
iterations:

t = 30(AO), w∗
30 = −4.7172;

t = 31(AO), w∗
31 = 4.1010;

t = 40(IO), w∗
40 = 4.0487;

t = 80(AO), w∗
80 = −5.6440;

t = 90(AO), w∗
90 = −4.6228

The outcome is consistent with our prearrangement. The outliers in {zt} process
are detected successfully, and there is no misjudgement.

6 Conclusions

There are AO and IO in our model, and also the outliers(AO) present consecu-
tively. It is quite difficult to detect outliers, however, all of the outliers in the
above model have been detected successfully, which shows our method is also
efficient for the AO and IO problem in ARMAX model. Some other case studies
also show that our method is effective in detecting outliers’ location and type and
in estimating their size for ARMAX model.

Acknowledgements
The research is supported by National Natural Science Foundation of China
(No.11171065)

References

[1] Baragona R, Battaglia F, Calzini C. (2001), “Genetic algorithms for the iden-
tification of additive and innovation outliers in time series”, Computational
Statistics & Data Analysis, Vol.37, pp.1-12.



Advances in Systems Science and Applications (2012) Vol.12 No.4 405

[2] Peña D, Sánchez I. (2005), “Multifold predictive validation in armax time
series models”, Journal of the American Statistical Association, Vol.100,
pp.135-146.

[3] Chen P, Li L, Liu Y and Lin J.G. (2010), “Detection of outliers and patch-
es in bilinear time series models”, Mathematical Problems in Engineering,
Vol.2010, pp.1-10.

[4] Chen P, Yang J, Li L.Y. (2013), “Synthetic detection of change point and
outliers in bilinear time series models”, International Journal of Systems
Science, http://dx.doi.org/10.1080/00207721.2013.777983. (In press)

[5] Huang L, Pang W, Wang K.P, Zhou C.G, Xiao Y. (2004), “Improved genet-
ic algorithm for vehicle routing problem with time windows”, Advances in
Systems Science and Applications, Vol.4, pp.118-124.

[6] Box G.E.P, Jenkins G.M, Reinsel G.C. (1994), Time Series Analysis: Fore-
casting and Control, third edition, Prentice-Hall, Englewood Cliffs, NJ.

Corresponding author
Author can be contacted at cp18@263.net.cn