|
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.
MODEL
Array of length 4 containing the order 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)
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)
Generic: CALL TS_OUTLIER_IDENTIFICATION (MODEL, W, X [, ])
Specific: The specific interface names are S_TS_OUTLIER_IDENTIFICATION and D_TS_OUTLIER_IDENTIFICATION.
Consider a univariate time series that can be described by the following multiplicative seasonal
ARIMA model of order
:
Here, ,
.
is the lag operator,
,
is a white
noise process, and
denotes the mean of
the series
.
In general, 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 Lius algorithm is unable
to determine the outliers 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 into account, Chen and
Liu consider the following model:
Here, is the observed
outlier contaminated series, and
and
denote the magnitude and dynamic pattern of outlier
, respectively.
is an indicator function that determines the temporal
course of the outlier effect,
,
otherwise. Note that
operates on
via
.
The last formula shows that the outlier free series can be obtained from the original series
by removing all occurring outlier effects:
.
The different types of outliers are characterized by different values for :
1. for an innovational
outlier,
2. for an additive
outlier,
3. for a level shift
outlier and
4. for a temporary change
outlier.
TS_OUTLIER_IDENTIFICATION
is an implementation of Chen and Lius algorithm. It determines the coefficients
in 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 Lius original paper (1993).
Intermediate and final estimates for the coefficients in
and
are computed by routines NSPE and NSLSE.
If the roots of
or
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.
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,
,
Gaussian White noise, to the given series. TS_OUTLIER_IDENTIFICATION
computes parameters
and
and identifies a LS outlier 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(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
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
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
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
PHONE: 713.784.3131 FAX:713.781.9260 |