TS_OUTLIER_IDENTIFICATION Function
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.
Usage
result = TS_OUTLIER_IDENTIFICATION (model, w)
Input Parameters
model—Four element array containing the numbers p, q, s, d of the ARIMA 
model the outlier free series is following.
model the outlier free series is following.w—Array containing all the observations (n_obs) in the time series.
Returned Value
result—Array of length n_obs containing the outlier free time series.
Input Keywords
Double—If present and nonzero, double precision is used.
Delta—The dampening effect parameter used in the detection of a Temporary Change Outlier (TC), 0 < Delta < 1. Default: Delta = 0.7
Critical—Critical value used as a threshold for outlier detection, Critical > 0. Default: Critical = 3.0
Epsilon—Positive tolerance value controlling the accuracy of parameter estimates during outlier detection. Default: Epsilon = 0.001
Relative_error—Stopping criterion for the nonlinear equation solver used in function ARMA. Default: Relative_error = 10–10.
Output Keywords
Residual—Array of length n_obs containing the residuals for the outlier free series.
Res_sigma—Residual standard error of the outlier free series.
Num_outliers—Scalar value indicating the number of outliers detected.
Outlier_statistics—2D Array Num_outliers by 2 containing outlier statistics. The first column contains the time at which the outlier was observed (t = 1, 2, ..., n_obs) and the second column contains an identifier indicating the type of outlier observed. Outlier types fall into one of five categories:
0—Innovational Outliers (IO)1—Additive outliers (AO)2—Level Shift Outliers (LS)3—Temporary Change Outliers (TC)4—Unable to Identify (UI).If Num_outliers = 0, a scalar –1 is returned.
Tau_statistics—Array of length Num_outliers containing the t value for each detected outlier. If Num_outliers = 0, a scalar –1 is returned.
Omega_weights—Array of length Num_outliers containing the computed ω weights for the detected outliers. If Num_outliers = 0, a scalar –1 is returned.
Arma_param—Array of length 1 + p + q containing the estimated constant, AR and MA parameters.
Aic—Scalar value containing Akaike’s information criterion (AIC).
Discussion
Consider a univariate time series {Yt} that can be described by the following multiplicative seasonal ARIMA model of order 
:
:
Here, 
, 
, and 
. B is the lag operator, 
, {at} is a white noise process, and μ denotes the mean of the series {Yt}.
, , and . B is the lag operator, , {at} 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, 
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, 
, It(tj) = 0 otherwise.
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, , It(tj) = 0 otherwise. Note that Lj(B) operates on It via:
The last formula shows that the outlier free series {Yt} can be obtained from the original series 
by removing all occurring outlier effects:
by removing all occurring outlier effects:
The different types of outliers are charaterized by different values for Lj(B):
1. 
for an innovational outlier
for an innovational outlier2. Lj(B) = 1 for an additive outlier
3. Lj(B) = (1 – B)–1 for a level shift outlier
4. 
for a temporary change outlier
for a temporary change outlier Function 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 functions ARMA and MAX_ARMA. 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.
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 (1 – B)2(1 – φ1B – φ2B2)Yt = at, 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 an LS outlier at time point t = 16.
series = [ 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]
n_obs = N_ELEMENTS(series) ; 114
model = [2, 0, 1, 2]
result = TS_OUTLIER_IDENTIFICATION( $
model, series, Critical=3.5, $
Num_outliers=num_outliers, $
Outlier_statistics=outlier_stat, $
Arma_param=parameters, $
Res_sigma=res_sigma, $
Aic=aic)
PRINT, "ARMA parameters:"
p = model(0) + model(1)
PRINT, parameters(0:p), Format='(F11.6)'
PRINT, ''
PRINT, num_outliers, Format="('Number of outliers: ', I1)"PRINT, ''
PRINT, "Outlier statistics:"
PRINT, "Time point Outlier type"
FOR i=0L, num_outliers-1 DO $
PRINT, outlier_stat(2*i), outlier_stat(2*i+1), $
Format="(I6, 10X, I3)"
PRINT, ''
PRINT, res_sigma, Format="('RSE: ', F10.6)"PRINT, aic, Format="('AIC: ', F10.6)"PRINT, ''
; Print out the first 36 values of result and series
PRINT, "Extract from the series:"
PRINT, ''
PRINT, "Time point Original Series Outlier free series"
FOR i=0L, 35 DO $
PRINT, (i+1), series(i), result(i), $
Format='(I6, F17.6, F18.6)'
Output
ARMA parameters:
0.000000
0.106532
-0.195856
Number of outliers: 1
Outlier statistics:
Time point Outlier type
16 2
RSE: 0.319542
AIC: 282.918152
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.699728
17 3.429000 2.769728
18 3.533000 2.873728
19 3.261000 2.601728
20 2.612000 1.952728
21 2.179000 1.519728
22 1.653000 0.993728
23 1.832000 1.172728
24 2.328000 1.668728
25 2.737000 2.077728
26 3.014000 2.354728
27 3.328000 2.668728
28 3.404000 2.744728
29 2.981000 2.321728
30 2.557000 1.897728
31 2.576000 1.916728
32 2.352000 1.692728
33 2.556000 1.896728
34 2.864000 2.204728
35 3.214000 2.554728
36 3.435000 2.775728
Example 2
This example is an artificial realization of an ARMA(1,1) process via formula Yt – 0.8Yt–1 = 10.0 + at + 0.5at–1, t = 1, ..., 300, {at} Gaussian white noise, E[Yt] = 50.0. An additive outlier with ω 1 = 4.5 was added at time point t = 150, a temporary change outlier with ω 2 = 3.0 was added at time point t = 200.
n_obs = 300
series = [ 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]
result = TS_OUTLIER_IDENTIFICATION( $
model, series, $
Num_outliers=num_outliers, $
Outlier_statistics=outlier_stat, $
Omega_weights=omega_user, $
Arma_param=parameters_user, $
Res_sigma=res_sigma, $
Aic=aic, $
Relative_Error=1.0e-5)
PRINT, "ARMA parameters:"
p = N_ELEMENTS(parameters_user)
PRINT, parameters_user(0:(p-1)), Format='(F11.6)'
PRINT, ''
PRINT, num_outliers, Format="('Number of outliers: ', I1)"PRINT, ''
PRINT, "Outlier statistics:"
PRINT, "Time point Outlier type"
FOR i=0L, num_outliers-1 DO $
PRINT, outlier_stat(i,0), outlier_stat(i,1), $
Format="(I6, 10X, I3)"
PRINT, ''
PRINT, "Omega statistics:"
PRINT, "Time point Omega"
FOR i=0L, num_outliers-1 DO $
PRINT, outlier_stat(i,0), omega_user(i), $
Format="(I6, 6X, F11.6)"
PRINT, ''
PRINT, res_sigma, Format="('RSE: ', F11.6)"PRINT, aic, Format="('AIC: ', F11.6)"Output
ARMA parameters:
10.828934
0.785221
-0.496502
Number of outliers: 2
Outlier statistics:
Time point Outlier type
150 1
200 3
Omega statistics:
Time point Omega
150 4.477869
200 3.381565
RSE: 1.007222
AIC: 1417.044067
