TS_OUTLIER_IDENTIFICATION
Detects and determines outliers and simultaneously estimates the model parameters in a time series whose underlying outlier free series follows a general seasonal or nonseasonal ARMA model.
Required Arguments
MODEL — Array of length 4 containing the order (p, 0, q) x (0, d, 0)s of the ARIMA model the outlier free series is following. Specifically, MODEL(1) = p, MODEL(2) = q, MODEL(3) = s, MODEL(4) = d. (Input)
W — Array of length NOBS containing the original time series. (Input)
X — Array of length NOBS containing the outlier free series. (Output)
Optional Arguments
NOBS — Number of observations in time series W. (Input)
Default: NOBS = size (W).
DELTA — The dynamic dampening effect parameter used in the detection of a Temporary Change Outlier (TC), 0.0 < DELTA < 1.0. (Input)
Default: DELTA = 0.7 .
CRITICAL — Critical value used as a threshold for outlier detection, CRITICAL > 0. (Input)
Default: CRITICAL = 3.0.
EPSILON — Positive tolerance value controlling the accuracy of parameter estimates during outlier detection. (Input)
Default: EPSILON = 0.001.
RELERR — Stopping criterion for use in the nonlinear equation solver used by NSPE, see routine NSPE for more details. (Input)
Default: RELERR = 1.0e‑10.
TOLSS — Tolerance level used to determine convergence of the nonlinear least‑squares algorithm used by NSLSE, see routine NSLSE for more details. (Input)
TOLSS must be greater than zero and less than one.
Default: TOLSS = 0.9 × AMACH(4).
RESIDUAL — Array of length NOBS containing the residuals for the outlier free series. (Output)
RESSIGMA — Residual standard error of the outlier free series. (Output)
NOUTLIERS — The number of outliers detected. (Output)
IOUTLIERSTATS — Pointer to an array of size NOUTLIERS by 2 containing outlier statistics. The first column contains the time at which the outlier was observed (
) and the second column contains an identifier indicating the type of outlier observed. Outlier types fall into one of five categories:
|
IOUTLIERSTATS |
Category |
|
0 |
Innovational Outliers (IO) |
|
1 |
Additive Outliers (AO) |
|
2 |
Level Shift Outliers(LS) |
|
3 |
Temporary Change Outliers(TC) |
|
4 |
Unable to Identify(UI) |
If no outliers are detected, an array of size 0 is returned. (Output)
TAUSTAT — Pointer to an array of length NOUTLIERS containing the t value for each detected outlier. (Output)
OMEGA — Pointer to an array of length NOUTLIERS containing the computed 
weights for the detected outliers. (Output)
PARAMS — Array of length 1+p+q containing the estimated constant, AR and MA parameters, respectively. (Output)
AIC — Akaike's Information Criterion (AIC) for the fitted model. (Output)
AICC— Akaike's Corrected Information Criterion (AICC) for the fitted model. (Output)
BIC — Bayesian Information Criterion (BIC) for the fitted model. (Output)
FORTRAN 90 Interface
Generic: CALL TS_OUTLIER_IDENTIFICATION (MODEL, W, X [, …])
Specific: The specific interface names are S_TS_OUTLIER_IDENTIFICATION and D_TS_OUTLIER_IDENTIFICATION.
Description
Consider a univariate time series {Yt} that can be described by the following multiplicative seasonal ARIMA model of order (p, 0, q) x (0, d, 0)s:
Here, 
, 
. B is the lag operator, 
, 
is a white noise process, and μ denotes the mean of the series {Yt}.
In general, {Yt} is not directly observable due to the influence of outliers. Chen and Liu (1993) distinguish between four types of outliers: innovational outliers (IO), additive outliers (AO), temporary changes (TC) and level shifts (LS). If an outlier occurs as the last observation of the series, then Chen and Liu’s algorithm is unable to determine the outlier’s classification. In TS_OUTLIER_IDENTIFICATION, such an outlier is called a UI (unable to identify) and is treated as an innovational outlier.
In order to take the effects of multiple outliers occurring at time points t1, t2, …, tm into account, Chen and Liu consider the following model:
Here, {Yt*} is the observed outlier contaminated series, and ωj and Lj(B) denote the magnitude and dynamic pattern of outlier j, respectively. It(tj) is an indicator function that determines the temporal course of the outlier effect, Itj(tj) = 1, It(tj) = 0 otherwise. Note that Lj(B) operates on It via BkIt = It-k, k = 0, 1, … .
The last formula shows that the outlier free series {Yt} can be obtained from the original series {Yt*} by removing all occurring outlier effects:
The different types of outliers are characterized by different values for 
:
-
for an innovational outlier, -
for an additive outlier, -
for a level shift outlier and -
for a temporary change outlier.
TS_OUTLIER_IDENTIFICATION is an implementation of Chen and Liu’s algorithm. It determines the coefficients inφ (B), θ (B), and the outlier effects in the model for the observed series jointly in three stages. The magnitude of the outlier effects is determined by least squares estimates. Outlier detection itself is realized by examination of the maximum value of the standardized statistics of the outlier effects. For a detailed description, see Chen and Liu’s original paper (1993).
Intermediate and final estimates for the coefficients in φ (B) and θ (B) are computed by routines NSPE and NSLSE. If the roots of φ (B) or θ (B) lie on or within the unit circle, then the algorithm stops with an appropriate error message. In this case, different values for p and q should be tried.
Examples
Example 1
This example is based on estimates of the Canadian lynx population. TS_OUTLIER_IDENTIFICATION is used to fit an ARIMA(2,2,0) model of the form 
, t = 1, 2, ..., 144, {at} Gaussian White noise, to the given series. TS_OUTLIER_IDENTIFICATION computes parameters φ1 = 0.123609 and φ2 = ‑0.178963, and identifies a LS outlier at time point t = 16.
use umach_int
use ts_outlier_identification_int
implicit none
integer :: i, nout
integer :: noutliers
integer, dimension(4) :: model
integer, dimension(:,:), pointer :: outlierstat
real(kind(1e0)) :: ressigma, aic
real(kind(1e0)), dimension(3) :: parameters
real(kind(1e0)), dimension(114) :: w, x
w = (/ 0.24300E01,0.25060E01,0.27670E01,0.29400E01,0.31690E01,&
0.34500E01,0.35940E01,0.37740E01,0.36950E01,0.34110E01,&
0.27180E01,0.19910E01,0.22650E01,0.24460E01,0.26120E01,&
0.33590E01,0.34290E01,0.35330E01,0.32610E01,0.26120E01,&
0.21790E01,0.16530E01,0.18320E01,0.23280E01,0.27370E01,&
0.30140E01,0.33280E01,0.34040E01,0.29810E01,0.25570E01,&
0.25760E01,0.23520E01,0.25560E01,0.28640E01,0.32140E01,&
0.34350E01,0.34580E01,0.33260E01,0.28350E01,0.24760E01,&
0.23730E01,0.23890E01,0.27420E01,0.32100E01,0.35200E01,&
0.38280E01,0.36280E01,0.28370E01,0.24060E01,0.26750E01,&
0.25540E01,0.28940E01,0.32020E01,0.32240E01,0.33520E01,&
0.31540E01,0.28780E01,0.24760E01,0.23030E01,0.23600E01,&
0.26710E01,0.28670E01,0.33100E01,0.34490E01,0.36460E01,&
0.34000E01,0.25900E01,0.18630E01,0.15810E01,0.16900E01,&
0.17710E01,0.22740E01,0.25760E01,0.31110E01,0.36050E01,&
0.35430E01,0.27690E01,0.20210E01,0.21850E01,0.25880E01,&
0.28800E01,0.31150E01,0.35400E01,0.38450E01,0.38000E01,&
0.35790E01,0.32640E01,0.25380E01,0.25820E01,0.29070E01,&
0.31420E01,0.34330E01,0.35800E01,0.34900E01,0.34750E01,&
0.35790E01,0.28290E01,0.19090E01,0.19030E01,0.20330E01,&
0.23600E01,0.26010E01,0.30540E01,0.33860E01,0.35530E01,&
0.34680E01,0.31870E01,0.27230E01,0.26860E01,0.28210E01,&
0.30000E01,0.32010E01,0.34240E01,0.35310E01 /)
model = (/ 2,0,1,2 /)
call ts_outlier_identification(model, w, x, CRITICAL=3.5, &
RESSIGMA=ressigma, NOUTLIERS=noutliers,&
IOUTLIERSTATS=outlierstat, PARAMS=parameters,&
AIC=aic)
call umach(2, nout)
write(nout,FMT="(T2,A,I2)") 'Number of outliers:', noutliers
write(nout,FMT="(T2,A)") 'Outlier statistics:'
write(nout,FMT="(T4,A,TR10,A)") 'Time point','Outlier type'
write(nout,FMT="(I8,TR20,I2)") (outlierstat(i,:),i=1,noutliers)
write(nout,FMT="(/,T2,A)") 'ARMA parameters:'
write(nout,FMT="(T3,I2,TR5,f10.6)") (i,parameters(i),i=1,3)
write(nout,FMT="(/,T2,A,f10.6)") 'RSE: ', ressigma
write(nout,FMT="(T2,A,f10.6)") 'AIC: ', aic
write(nout,FMT="(/,T2,A)") 'Extract from the series:'
write(nout,FMT="(T2,A,TR6,A,TR6,A)") 'time point',&
'original series', 'outlier free series'
do i=1,36
write(nout, FMT="(T5,I2,T15,F15.6,T35,F19.6)") i, w(i), x(i)
end do
end
Output
Number of outliers: 1
Outlier statistics:
Time point Outlier type
16 2
ARMA parameters:
1 0.000000
2 0.124390
3 -0.179959
RSE: 0.319650
AIC: 282.995575
Extract from the series:
time point original series outlier free series
1 2.430000 2.430000
2 2.506000 2.506000
3 2.767000 2.767000
4 2.940000 2.940000
5 3.169000 3.169000
6 3.450000 3.450000
7 3.594000 3.594000
8 3.774000 3.774000
9 3.695000 3.695000
10 3.411000 3.411000
11 2.718000 2.718000
12 1.991000 1.991000
13 2.265000 2.265000
14 2.446000 2.446000
15 2.612000 2.612000
16 3.359000 2.701984
17 3.429000 2.771984
18 3.533000 2.875984
19 3.261000 2.603984
20 2.612000 1.954984
21 2.179000 1.521984
22 1.653000 0.995984
23 1.832000 1.174984
24 2.328000 1.670984
25 2.737000 2.079984
26 3.014000 2.356984
27 3.328000 2.670984
28 3.404000 2.746984
29 2.981000 2.323984
30 2.557000 1.899984
31 2.576000 1.918984
32 2.352000 1.694984
33 2.556000 1.898984
34 2.864000 2.206984
35 3.214000 2.556984
36 3.435000 2.777984
Example 2
This example is an artificial realization of an ARMA(1,1) process via formula 
Gaussian white noise, 
.
An additive outlier with 
was added at time point 
, a temporary change outlier with 
was added at time point 
.
use umach_int
use ts_outlier_identification_int
implicit none
integer :: i, nout
integer :: noutliers
integer, dimension(4) :: model
integer, dimension(:,:), pointer :: outlierstat
real(kind(1e0)) :: ressigma, aic
real(kind(1e0)), dimension(3) :: parameters
real(kind(1e0)), dimension(300) :: w, x
real(kind(1e0)), dimension(:), pointer :: omega
w = (/ 50.0000000,50.2728081,50.6242599,51.0373917,51.9317627,&
50.3494759,51.6597252,52.7004929,53.5499802,53.1673279,&
50.2373505,49.3373871,49.5516472,48.6692696,47.6606636,&
46.8774185,45.7315445,45.6469727,45.9882355,45.5216560,&
46.0479660,48.1958656,48.6387749,49.9055367,49.8077278,&
47.7858467,47.9386749,49.7691956,48.5425873,49.1239853,&
49.8518791,50.3320694,50.9146347,51.8772049,51.8745689,&
52.3394470,52.7273712,51.4310036,50.6727448,50.8370399,&
51.2843437,51.8162918,51.6933670,49.7038231,49.0189247,&
49.455703,50.2718010,49.9605980,51.3775749,50.2285385,&
48.2692299,47.6495590,49.2938499,49.1924858,49.6449242,&
50.0446815,51.9972496,54.2576981,52.9835434,50.4193535,&
50.3617897,51.8276901,53.1239929,54.0682144,54.9238319,&
55.6877632,54.8896332,54.0701065,52.2754097,52.2522354,&
53.1248703,51.1287193,50.5003815,49.6504173,47.2453079,&
45.4555626,45.8449707,45.9765129,45.7682228,45.2343674,&
46.6496811,47.0894432,49.3368340,50.8058052,49.9132500,&
49.5893288,48.2470627,46.9779968,45.6760864,45.7070389,&
46.6158409,47.5303612,47.5630417,47.0389214,46.0352287,&
45.8161545,45.7974396,46.0015373,45.3796463,45.3461685,&
47.6444016,49.3327446,49.3810692,50.2027817,51.4567032,&
52.3986320,52.5819206,52.7721825,52.6919098,53.3274345,&
55.1345940,56.8962631,55.7791634,55.0616989,52.3551178,&
51.3264084,51.0968323,51.1980476,52.8001442,52.0545082,&
50.8742943,51.5150337,51.2242050,50.5033989,48.7760124,&
47.4179192,49.7319527,51.3320541,52.3918304,52.4140434,&
51.0845947,49.6485748,50.6893463,52.9840813,53.3246994,&
52.4568024,51.9196091,53.6683121,53.4555359,51.7755814,&
49.2915611,49.8755112,49.4546776,48.6171913,49.9643021,&
49.3766441,49.2551308,50.1021881,51.0769119,55.8328133,&
52.0212708,53.4930801,53.2147255,52.2356453,51.9648819,&
52.1816330,51.9898071,52.5623627,51.0717278,52.2431946,&
53.6943054,54.3752098,54.1492615,53.8523254,52.1093712,&
52.3982697,51.2405128,50.3018112,51.3819618,49.5479546,&
47.5024452,47.4447708,47.8939056,48.4070015,48.2440681,&
48.7389755,49.7309227,49.1998024,49.5798340,51.1196213,&
50.6288414,50.3971405,51.6084099,52.4564743,51.6443901,&
52.4080658,52.4643364,52.6257210,53.1604691,51.9309731,&
51.4137230,52.1233368,52.9867249,53.3180733,51.9647636,&
50.7947655,52.3815842,50.8353729,49.4136009,52.8355217,&
52.2234840,51.1392517,48.5245132,46.8700218,46.1607285,&
45.2324257,47.4157829,48.9989090,49.6230736,50.4352913,&
51.1652985,50.2588654,50.7820129,51.0448799,51.2880516,&
49.6898804,49.0288200,49.9338837,48.2214432,46.2103348,&
46.9550171,47.5595894,47.7176018,48.4502945,50.9816895,&
51.6950073,51.6973495,52.1941261,51.8988075,52.5617599,&
52.0218391,49.5236053,47.9684906,48.2445183,48.8275146,&
49.7176971,51.5649338,52.5627213,52.0182419,50.9688835,&
51.5846901,50.9486771,48.8685837,48.5600624,48.4760094,&
48.5348396,50.4187813,51.2542381,50.1872864,50.4407692,&
50.6222687,50.4972000,51.0036087,51.3367500,51.7368202,&
53.0463791,53.6261253,52.0728683,48.9740753,49.3280830,&
49.2733917,49.8519020,50.8562126,49.5594254,49.6109200,&
48.3785629,48.0026474,49.4874268,50.1596375,51.8059540,&
53.0288620,51.3321075,49.3114815,48.7999306,47.7201881,&
46.3433914,46.5303612,47.6294632,48.6012459,47.8567657,&
48.0604057,47.1352806,49.5724792,50.5566483,49.4182968,&
50.5578079,50.6883736,50.6333389,51.9766159,51.0595245,&
49.3751640,46.9667702,47.1658173,47.4411278,47.5360374,&
48.9914742,50.4747620,50.2728043,51.9117165,53.7627792 /)
model = (/ 1,1,1,0 /)
call ts_outlier_identification(model, w, x, CRITICAL=3.5,&
RESSIGMA=ressigma, NOUTLIERS=noutliers,&
IOUTLIERSTATS=outlierstat,OMEGA=omega,&
PARAMS=parameters, AIC=aic, RELERR=1.0e-05)
call umach(2, nout)
write(nout,FMT="(/,T2,A)") 'ARMA parameters:'
write(nout,FMT="(T3,I2,TR5,f10.6)") (i,parameters(i),i=1,3)
write(nout,FMT="(/,T2,A,I2)") 'Number of outliers:', noutliers
write(nout,FMT="(/,T2,A)") 'Outlier statistics:'
write(nout,FMT="(T4,A,TR10,A)") 'Time point','Outlier type'
write(nout,FMT="(I8,TR20,I2)") (outlierstat(i,:),i=1,noutliers)
write(nout,FMT="(/,T2,A)") 'Omega statistics:'
write(nout,FMT="(T4,A,TR10,A)") 'Time point','Omega'
write(nout,FMT="(I9,TR10,f10.6)") (outlierstat(i,1),omega(i),&
i=1,noutliers)
write(nout,FMT="(/,T2,A,f10.6)") 'RSE: ', ressigma
write(nout,FMT="(T2,A,f12.6)") 'AIC: ', aic
end
Output
ARMA parameters:
1 10.700689
2 0.787765
3 -0.498039
Number of outliers: 2
Outlier statistics:
Time point Outlier type
150 1
200 3
Omega statistics:
Time point Omega
150 4.478198
200 3.381984
RSE: 1.007209
AIC: 1417.036377
