|
Automatically identifies time series outliers, determines
parameters of a multiplicative seasonal ARIMA model and produces forecasts that incorporate the effects
of outliers whose effects persist beyond the end of the series.
ITIME_POINTS
Array of length NOBS containing the
time points at which the time
series was observed. Time points must be integer and in strictly ascending
order. It is assumed that the time points of the time series after estimation of
the missing values are equidistant with distance 1 between two consecutive time
points. (Input)
W Array of
length NOBS
containing the observed time series values . This series can contain outliers and missing observations.
Outliers are identified by this routine and missing values are identified by the
time values in array ITIME_POINTS. If the
time interval between two consecutive time points is greater than one,
i.e.
, then
missing values are assumed to exist between
and
at times
. Therefore, the gap free series is assumed to be defined for
equidistant time points . Missing values are automatically estimated prior to
identifying outliers and producing forecasts. Forecasts are generated for
both missing and observed values. (Input)
PARAMS
Allocatable array of length 1+p+q containing the estimated
constant, AR and MA parameters of the adjusted optimum seasonal ARIMA model. If
, then an ARMA(p, q)
model is fitted to the outlier-free version of the observed series
. If
, these parameters are
computed for an ARMA(p,q) representation of the seasonally
adjusted series
, where
and
. (Output)
NOBS Number of
observations in the time series. Assuming that the series is defined at time
points , the actual length of
the series, including missing values, is
. (Input)
Default: NOBS = size (W).
IMETH
Method to be used in model selection. (Input)
1 Automatic selection
2 Grid search (requires arguments IAR and IMA)
3 Specified
model (Requires
argument MODEL)
Default:
IMETH = 1.
MAXLAG Maximum
number of autoregressive parameters allowed in fitting an AR model to the
series. (Input)
Default: MAXLAG = 10.
INFOCRIT The information criterion used for optimum model selection. (Input)
INFOCRIT |
selected information criterion |
0 |
Akaikes Information Criterion (AIC) |
1 |
Akaikes Corrected Information Criterion (AICC) |
2 |
Bayesian Information Criterion (BIC) |
Default: INFOCRIT = 0.
DELTA Dynamic
dampening effect parameter used in the detection of a Temporary Change Outlier
(TC). (Input)
It is required that 0.0 < DELTA <
1.0.
Default: DELTA = 0.7.
CRITICAL
Critical value used as a threshold for outlier detection.
(Input)
CRITICAL must be
greater than zero.
Default: CRITICAL = 3.0.
EPSILON
Positive tolerance value controlling the accuracy of parameter estimates during
outlier detection. (Input)
Default: EPSILON =
0.001.
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).
IAR Array containing the candidate values for the AR order p from which the optimum is being selected. All candidate values in IAR must be non-negative and IAR must contain at least one value. If IMETH = 2 then IAR must be defined. Otherwise, IAR is ignored. (Input)
IMA Array containing the candidate values for the MA order q from which the optimum is being selected. All candidate values in IMA must be non-negative and IMA must contain at least one value. If IMETH = 2 then IMA must be defined. Otherwise, IMA is ignored. (Input)
IPER Array
containing the candidate values for s from which the optimum is being
selected. All candidate values in IPER must be positive
and IPER must
contain at least one value. (Input)
Default: IPER(:) = 1
IORD Array
containing the candidate values for d from which the optimum is being
selected. All candidate values in IORD must be
non-negative and IORD must contain at
least one value. (Input)
Default: IORD (:) = 0
ALPHA Value in
the exclusive interval (0,1) used to specify the 100(1-ALPHA)% probability
limits of the forecast. (Input)
Default: ALPHA = 0.05.
MXLEAD Maximum
lead time for forecasts. (Input)
Default: MXLEAD = 0.
MODEL Array
of length 4 containing values for p, q, s, and
d. (Input/Output)
For IMETH = 1 or IMETH = 2 ,
MODEL is ignored
on input. If IMETH = 3 then
p and q must be defined on input. If IPER and IORD are not defined
then s and d must also be defined on input. On output, MODEL contains optimum
values for p, q, s, and d in MODEL(1), MODEL(2), MODEL(3) and MODEL(4),
respectively.
RESIDUAL Array
of length containing estimates
for the white noise in the gap-free and outlier-free original series.
(Output)
RSE Residual standard error (RSE) of the outlier-free and gap-free original series. (Output)
NOUTLIERS Number of outliers detected. (Output)
IOUTLIERSTATS
Pointer to an array of size NOUTLIERS by 2
containing the outlier statistics. The first column contains the time point at
which the outlier was observed (time points ranging from to
) and the second
column contains an identifier indicating the type of outlier observed. Outlier
types fall into one of five categories:
IOUTLIERSTATS |
Catigory |
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, then an array of size 0 is returned. (Output)
AIC Akaike's Information Criterion (AIC) for the fitted optimum model. Uses an approximation of the maximum log-likelihood based on an estimate of the innovation variance of the series. (Output)
AICC Akaike's Corrected Information Criterion (AICC) for the fitted optimum model. Uses an approximation of the maximum log-likelihood based on an estimate of the innovation variance of the series. (Output)
BIC Bayesian Information Criterion (BIC) for the fitted optimum model. Uses an approximation of the maximum log-likelihood based on an estimate of the innovation variance of the series. (Output)
OUTFREESERIES
Array of size N
by 2 containing the adjusted time series. The first column contains the NOBS observations from
the original series, , plus estimated values for
any time gaps. The second column contains the same values as the first column
adjusted by removing any outlier effects. In effect, the second column contains
estimates of the underlying outlier-free series,
. If no outliers are detected then both columns will
contain identical values. (Output)
OUTFREEFCST
Array of size MXLEAD by 3 containing
the forecasted values for the original outlier and gap free series at origin
for lead times
in the first column. The second column contains standard
errors for these forecasts, and the third column contains the
weights of the infinite order moving average form of the
model. (Output)
OUTLIERFCST
Array of size MXLEAD by 3 containing
the forecasted values for the original and gap free series for in the first column. The second column contains standard errors
for these forecasts, and the third column contains the
weights of the infinite order moving average form of the
model. (Output)
Generic: CALL AUTO_ARIMA (ITIME_POINTS, W, PARAMS [, ])
Specific: The specific interface names are S_AUTO_ARIMA and D_AUTO_ARIMA.
Routine AUTO_ARIMA
determines the parameters of a multiplicative seasonal ARIMA model, and then uses the fitted model to identify
outliers and prepare forecasts. The order of this model can be specified or
automatically determined.
The ARIMA model handled by AUTO_ARIMA
has the following form:
where
and
It is assumed that all roots of and
lie outside the unit
circle. Clearly, if
this reduces to the
traditional ARIMA(p, d, q) model.
is the unobserved, gap-free and outlier-free time series with
mean
, and white noise
. This model is referred to as the underlying, outlier-free
model. Routine AUTO_ARIMA
does not assume that this series is observable. It assumes that the observed
values might be contaminated by one or more outliers, whose effects are added to
the underlying outlier-free series:
Outlier identification uses the algorithm developed by Chen and Liu (1993). Outliers are classified into 1 of 5 types:
1. innovational
2. additive
3. level shift
4. temporary change and
5. unable to identify
Once outliers
are identified, AUTO_ARIMA
estimates , the outlier-free series
representation of the data, by removing the estimated outlier effects.
Using the
information about the adjusted ARIMA model and the removed outliers, forecasts are then
prepared for the outlier-free series. Outlier effects are added to these
forecasts to produce a forecast for the observed series,
. If there are no outliers, then the forecasts for the
outlier-free series and the observed series will be identical.
Users have an option of either specifying specific values for p, q, s and d or have AUTO_ARIMA automatically select best fit values. Model selection can be conducted in one of three methods listed below depending upon the value of optional argument IMETH.
This method initially searches for the AR(p) representation with minimum AIC for the noisy data, where p = 0,..., MAXLAG.
If IORD
is defined then the values in IPER
and IORD
are included in the search to find an optimum ARIMA representation of the series. Here, every possible
combination of values for p, s in IPER
and d in IORD
is examined. The best found ARIMA
representation is then
used as input for the outlier detection routine.
The optimum values for p, q, s and d are returned in MODEL(1), MODEL(2), MODEL(3)and MODEL(4), respectively.
The second automatic method conducts a grid search for p and q using all possible combinations of candidate values in IAR and IMA. Therefore, for this method the definition of IAR and IMA is required.
If IORD is defined, the grid search is extended to include the candidate values for s and d given in IPER and IORD, respectively.
If IORD is not defined, no seasonal adjustment is attempted, and the grid search is restricted to searching for optimum values of p and q only.
The optimum values of p, q, s and d are returned in MODEL(1), MODEL(2), MODEL(3) and MODEL(4), respectively.
In the third
method, specific values for p, q, s and d are given. The values for p and q must be defined in MODEL(1) and MODEL(2),
respectively. If IPER
and IORD
are not defined, then values and
must be specified in
MODEL(3)
and MODEL(4). If IPER
and IORD are defined, then a grid search for the
optimum values of s and
d is conducted using all
possible combinations of input values in IPER
and IORD.
The optimum values of s
and d can be found in
MODEL(3)
and MODEL(4),
respectively.
The algorithm of Chen and Liu (1993) is used to identify outliers. The number of outliers identified is returned in NOUTLIERS. Both the time and classification for these outliers are returned in IOUTLIERSTATS. Outliers are classified into one of five categories based upon the standardized statistic for each outlier type. The time at which the outlier occurred is given in the first column of IOUTLIERSTATS. The outlier identifier returned in the second column is according to the descriptions in the following table:
Outlier Identifier |
Name |
General Description |
0 |
(IO) Innovational Outlier |
Innovational
outliers persist. That is, there is an initial impact at the time
the outlier occurs. This effect continues in a lagged fashion with
all future observations. The lag coefficients are determined by the
coefficient of the underlying ARIMA |
1 |
(AO) Additive Outlier |
Additive outliers do not persist. As the name implies, an additive outlier affects only the observation at the time the outlier occurs. Hence additive outliers have no effect on future forecasts. |
2 |
(LS) |
Level shift outliers persist. They have the effect of either raising or lowering the mean of the series starting at the time the outlier occurs. This shift in the mean is abrupt and permanent. |
3 |
(TC) |
Temporary change outliers persist and are similar to level shift outliers with one major exception. Like level shift outliers, there is an abrupt change in the mean of the series at the time this outlier occurs. However, unlike level shift outliers, this shift is not permanent. The TC outlier gradually decays, eventually bringing the mean of the series back to its original value. The rate of this decay is modeled using the parameter DELTA. The default of DELTA= 0.7 is the value recommended for general use by Chen and Liu (1993). |
4 |
(UI) Unable to Identify |
If an outlier is identified as the last observation, then the algorithm is unable to determine the outliers classification. For forecasting, a UI outlier is treated as an IO outlier. That is, its effect is lagged into the forecasts. |
Except for additive outliers (AO), the effect of an outlier persists to observations following that outlier. Forecasts produced by AUTO_ARIMA take this into account.
This example uses time series LNU03327709 from the US Department of Labor, Bureau of Labor Statistics. It contains the unadjusted special unemployment rate, taken monthly from January 1994 through September 2005. The values 01/2004 03/2005 are used by AUTO_ARIMA for outlier detection and parameter estimation. In this example , Method 1 without seasonal adjustment is chosen to find an appropriate AR(p) model. A forecast is done for the following six months and compared with the actual values 04/2005 09/2005.
use auto_arima_int
use wrrrl_int
use umach_int
implicit none
! Specifications for parameters
integer :: nobs, mxlead, i, noutliers, nout
integer, dimension(4) :: model
integer, dimension(:,:), pointer :: outlierstat
integer, dimension(141) :: times
real(kind(1e0)) :: aic, rse
real(kind(1e0)), dimension(141) :: x
real(kind(1e0)), dimension(6,3) :: outlierfcst
real(kind(1e0)), dimension(6,4) :: forecast_table
real(kind(1e0)), dimension(:), allocatable :: parameters
character (len=10), dimension(1) :: rlabel
character (len = 14), dimension(5) :: clabel
character (len = 10) :: fmt
! Time series data
x = (/ 12.8,12.2,11.9,10.9,10.6,11.3,11.1,10.4,10.0,9.7,9.7,&
9.7,11.1,10.5,10.3,9.8,9.8,10.4,10.4,10.0,9.7,9.3,&
9.6,9.7,10.8,10.7,10.3,9.7,9.5,10.0,10.0,9.3,9.0,&
8.8,8.9,9.2,10.4,10.0,9.6,9.0,8.5,9.2,9.0,8.6,&
8.3,7.9,8.0,8.2,9.3,8.9,8.9,7.7,7.6,8.4,8.5,&
7.8,7.6,7.3,7.2,7.3,8.5,8.2,7.9,7.4,7.1,7.9,&
7.7,7.2,7.0,6.7,6.8,6.9,7.8,7.6,7.4,6.6,6.8,&
7.2,7.2,7.0,6.6,6.3,6.8,6.7,8.1,7.9,7.6,7.1,&
7.2,8.2,8.1,8.1,8.2,8.7,9.0,9.3,10.5,10.1,9.9,&
9.4,9.2,9.8,9.9,9.5,9.0,9.0,9.4,9.6,11.0,10.8,&
10.4,9.8,9.7,10.6,10.5,10.0,9.8,9.5,9.7,9.6,10.9,&
10.3,10.4,9.3,9.3,9.8,9.8,9.3,8.9,9.1,9.1,9.1,&
10.2,9.9,9.4,8.7,8.6,9.3,9.1,8.8,8.5 /)
times = (/ 1,2,3,4,5,6,7,8,9,10,11,12,&
13,14,15,16,17,18,19,20,21,22,23,24,&
25,26,27,28,29,30,31,32,33,34,35,36,&
37,38,39,40,41,42,43,44,45,46,47,48,&
49,50,51,52,53,54,55,56,57,58,59,60,&
61,62,63,64,65,66,67,68,69,70,71,72,&
73,74,75,76,77,78,79,80,81,82,83,84,&
85,86,87,88,89,90,91,92,93,94,95,96,&
97,98,99,100,101,102,103,104,105,106,107,108,&
109,110,111,112,113,114,115,116,117,118,119,120,&
121,122,123,124,125,126,127,128,129,130,131,132,&
133,134,135,136,137,138,139,140,141 /)
rlabel(1) = 'NUMBER'
clabel(1) = ' '
clabel(2) = 'Series'
clabel(3) = 'Forecast'
clabel(4) = 'Prob. Limits'
clabel(5) = 'Psi Weights'
mxlead = 6
nobs = 135
call auto_arima (times, x, parameters, NOBS=nobs, MAXLAG=5,&
MODEL=model, AIC=aic, CRITICAL=4.0,&
NOUTLIERS=noutliers, IOUTLIERSTATS=outlierstat,&
RSE=rse, MXLEAD=mxlead, OUTLIERFCST=outlierfcst)
call umach (2,nout)
write (nout,*) 'Method 1: Automatic ARIMA model selection,'//&
' no differencing'
write (nout,FMT="(T2,4(A,I2))") 'Model chosen: p =', model(1),&
', q =', model(2), ', s =', model(3), ', d =', model(4)
write (nout,*)
write (nout,FMT="(T2,A,I2)") 'Number of outliers: ', noutliers
write (nout,*) 'Outlier statistics:'
write (nout,*) 'Time point Outlier type'
do i=1,noutliers
write (nout,FMT="(I11,T15,I13)") outlierstat(i,1),&
outlierstat(i,2)
end do
write (nout,*)
write (nout,*) 'AIC = ', aic
write (nout,*) 'RSE = ', rse
write (nout,*)
write(nout,FMT="(/,T2,A)") 'ARMA parameters:'
do i=1, 1+model(1)+model(2)
write (nout,FMT="(T3,I2,TR5,f10.6)") i, parameters(i)
end do
forecast_table(1:mxlead,1) = x(nobs+1:nobs+mxlead)
forecast_table(1:mxlead, 2:4) = outlierfcst(1:mxlead, 1:3)
write (nout,*)
fmt = '(F11.4)'
call wrrrl('* * * Forecast Table * * *', forecast_table,&
rlabel, clabel, FMT = fmt)
end
Method 1: Automatic ARIMA model selection, no differencing
Model chosen: p = 5, q = 0, s = 1, d = 0
Number of outliers: 4
Outlier statistics:
Time point Outlier type
8 2
13 0
97 0
109 0
AIC = 397.5339
RSE = 0.3966153
ARMA parameters:
1 0.481449
2 0.813321
3 -0.043181
4 -0.220261
5 0.199172
6 0.199179
* * * Forecast Table * * *
Series Forecast Prob. Limits Psi Weights
1 8.7000 9.0273 0.7774 0.8133
2 8.6000 9.0309 1.0020 0.6183
3 9.3000 9.3195 1.1113 0.2475
4 9.1000 9.4767 1.1278 0.1946
5 8.8000 9.4176 1.1379 0.3726
6 8.5000 9.2256 1.1742 0.5253
This is the
same as Example 1, except now AUTO_ARIMA uses Method 2 with a possible seasonal adjustment. As a result, the
unadjusted model with is chosen as
optimum.
use auto_arima_int
use wrrrl_int
use umach_int
implicit none
! Specifications for parameters
integer :: nobs, mxlead, i, noutliers, nout
integer, dimension(4) :: model
integer, dimension(2) :: iper = (/ 1,2 /)
integer, dimension(3) :: iord = (/ 0,1,2 /)
integer, dimension(4) :: ima = (/ 0,1,2,3 /)
integer, dimension(4) :: iar = (/ 0,1,2,3 /)
integer, dimension(141) :: times
integer, dimension(:,:), pointer :: outlierstat
real(kind(1e0)) :: aic, rse
real(kind(1e0)), dimension(:), allocatable :: parameters
real(kind(1e0)), dimension(6,3) :: outlierfcst
real(kind(1e0)), dimension(6,4) :: forecast_table
real(kind(1e0)), dimension(141) :: x
character (len = 10), dimension(1) :: rlabel
character (len = 14), dimension(5) :: clabel
character (len = 10) :: fmt
! Time series data
x = (/ 12.8,12.2,11.9,10.9,10.6,11.3,11.1,10.4,10.0,9.7,9.7,&
9.7,11.1,10.5,10.3,9.8,9.8,10.4,10.4,10.0,9.7,9.3,&
9.6,9.7,10.8,10.7,10.3,9.7,9.5,10.0,10.0,9.3,9.0,&
8.8,8.9,9.2,10.4,10.0,9.6,9.0,8.5,9.2,9.0,8.6,&
8.3,7.9,8.0,8.2,9.3,8.9,8.9,7.7,7.6,8.4,8.5,&
7.8,7.6,7.3,7.2,7.3,8.5,8.2,7.9,7.4,7.1,7.9,&
7.7,7.2,7.0,6.7,6.8,6.9,7.8,7.6,7.4,6.6,6.8,&
7.2,7.2,7.0,6.6,6.3,6.8,6.7,8.1,7.9,7.6,7.1,&
7.2,8.2,8.1,8.1,8.2,8.7,9.0,9.3,10.5,10.1,9.9,&
9.4,9.2,9.8,9.9,9.5,9.0,9.0,9.4,9.6,11.0,10.8,&
10.4,9.8,9.7,10.6,10.5,10.0,9.8,9.5,9.7,9.6,10.9,&
10.3,10.4,9.3,9.3,9.8,9.8,9.3,8.9,9.1,9.1,9.1,&
10.2,9.9,9.4,8.7,8.6,9.3,9.1,8.8,8.5 /)
times = (/ 1,2,3,4,5,6,7,8,9,10,11,12,&
13,14,15,16,17,18,19,20,21,22,23,24,&
25,26,27,28,29,30,31,32,33,34,35,36,&
37,38,39,40,41,42,43,44,45,46,47,48,&
49,50,51,52,53,54,55,56,57,58,59,60,&
61,62,63,64,65,66,67,68,69,70,71,72,&
73,74,75,76,77,78,79,80,81,82,83,84,&
85,86,87,88,89,90,91,92,93,94,95,96,&
97,98,99,100,101,102,103,104,105,106,107,108,&
109,110,111,112,113,114,115,116,117,118,119,120,&
121,122,123,124,125,126,127,128,129,130,131,132,&
133,134,135,136,137,138,139,140,141 /)
rlabel(1) = 'NUMBER'
clabel(1) = ' '
clabel(2) = 'Series'
clabel(3) = 'Forecast'
clabel(4) = 'Prob. Limits'
clabel(5) = 'Psi Weights'
mxlead = 6
nobs = 135
call auto_arima (times, x, parameters, NOBS=nobs, MAXLAG=5,&
MODEL=model, AIC=aic,CRITICAL=4.0,&
NOUTLIERS=noutliers, IMETH=2,&
IAR=iar, IMA=ima, IPER=iper, IORD=iord,&
IOUTLIERSTATS=outlierstat, RSE=rse,&
MXLEAD=mxlead, OUTLIERFCST=outlierfcst)
call umach (2, nout)
write (nout,*) 'Method 2: Grid search, differencing allowed'
write (nout,FMT="(T2,4(A,I2))") 'Model chosen: p =', model(1),&
', q =', model(2), ', s =', model(3), ', d =', model(4)
write (nout,*)
write (nout,FMT="(T2,A,I2)") 'Number of outliers: ', noutliers
write (nout,*) 'Outlier statistics:'
write (nout,*) 'Time point Outlier type'
do i=1,noutliers
write (nout, FMT="(I11, T15, I13)") outlierstat(i,1),&
outlierstat(i,2)
end do
write (nout,*)
write (nout,*) 'AIC = ', aic
write (nout,*) 'RSE = ', rse
write (nout,*)
write (nout,*) 'ARMA parameters:'
do i=1, 1+model(1)+model(2)
write (nout,FMT="(T3,I2,TR5,f10.6)") i, parameters(i)
end do
write (nout,*)
forecast_table(1:mxlead,1) = x(nobs+1:nobs+mxlead)
forecast_table(1:mxlead,2:4) = outlierfcst(1:mxlead, 1:3)
fmt = '(F11.4)'
call wrrrl('* * * Forecast Table * * *', forecast_table,&
rlabel, clabel, FMT = fmt)
end
Method 2: Grid search, differencing allowed
Model chosen: p = 3, q = 2, s = 1, d = 0
Number of outliers: 1
Outlier statistics:
Time point Outlier type
109 0
AIC = 408.0768
RSE = 0.4124085
ARMA parameters:
1 0.509427
2 1.944686
3 -1.901132
4 0.901670
5 1.113016
6 -0.915008
* * * Forecast Table * * *
Series Forecast Prob. Limits Psi Weights
1 8.7000 9.1109 0.8083 0.8317
2 8.6000 9.1811 1.0513 0.6312
3 9.3000 9.5185 1.1686 0.5481
4 9.1000 9.7804 1.2497 0.6157
5 8.8000 9.7117 1.3452 0.7245
6 8.5000 9.3842 1.4671 0.7326
This example is the same as Example 2 but
now Method 3 with the optimum model parameters from Example 2 is chosen for outlier detection and
forecasting.
use auto_arima_int
use wrrrl_int
use umach_int
implicit none
! Parameter specifications
integer :: nobs, mxlead, i, noutliers, nout
integer, dimension(4) :: model
integer, dimension(141) :: times
integer, dimension(:,:), pointer :: outlierstat
real(kind(1e0)) :: aic, rse
real(kind(1e0)), dimension(141) :: x
real(kind(1e0)), dimension(6,3) :: outlierfcst
real(kind(1e0)), dimension(6,4) :: forecast_table
real(kind(1e0)), dimension(:), allocatable :: parameters
character (len = 10), dimension(1) :: rlabel
character (len = 14), dimension(5) :: clabel
character (len = 10) :: fmt
! Time series data
x = (/ 12.8,12.2,11.9,10.9,10.6,11.3,11.1,10.4,10.0,9.7,9.7,&
9.7,11.1,10.5,10.3,9.8,9.8,10.4,10.4,10.0,9.7,9.3,&
9.6,9.7,10.8,10.7,10.3,9.7,9.5,10.0,10.0,9.3,9.0,&
8.8,8.9,9.2,10.4,10.0,9.6,9.0,8.5,9.2,9.0,8.6,&
8.3,7.9,8.0,8.2,9.3,8.9,8.9,7.7,7.6,8.4,8.5,&
7.8,7.6,7.3,7.2,7.3,8.5,8.2,7.9,7.4,7.1,7.9,&
7.7,7.2,7.0,6.7,6.8,6.9,7.8,7.6,7.4,6.6,6.8,&
7.2,7.2,7.0,6.6,6.3,6.8,6.7,8.1,7.9,7.6,7.1,&
7.2,8.2,8.1,8.1,8.2,8.7,9.0,9.3,10.5,10.1,9.9,&
9.4,9.2,9.8,9.9,9.5,9.0,9.0,9.4,9.6,11.0,10.8,&
10.4,9.8,9.7,10.6,10.5,10.0,9.8,9.5,9.7,9.6,10.9,&
10.3,10.4,9.3,9.3,9.8,9.8,9.3,8.9,9.1,9.1,9.1,&
10.2,9.9,9.4,8.7,8.6,9.3,9.1,8.8,8.5/)
times = (/ 1,2,3,4,5,6,7,8,9,10,11,12,&
13,14,15,16,17,18,19,20,21,22,23,24,&
25,26,27,28,29,30,31,32,33,34,35,36,&
37,38,39,40,41,42,43,44,45,46,47,48,&
49,50,51,52,53,54,55,56,57,58,59,60,&
61,62,63,64,65,66,67,68,69,70,71,72,&
73,74,75,76,77,78,79,80,81,82,83,84,&
85,86,87,88,89,90,91,92,93,94,95,96,&
97,98,99,100,101,102,103,104,105,106,107,108,&
109,110,111,112,113,114,115,116,117,118,119,120,&
121,122,123,124,125,126,127,128,129,130,131,132,&
133,134,135,136,137,138,139,140,141/)
mxlead = 6
nobs = 135
model = (/ 3, 2, 1, 0 /)
rlabel(1) = 'NUMBER'
clabel(1) = ' '
clabel(2) = 'Series'
clabel(3) = 'Forecast'
clabel(4) = 'Prob. Limits'
clabel(5) = 'Psi Weights'
call auto_arima (times, x, parameters, NOBS=nobs, MODEL=model,&
AIC=aic, CRITICAL=4.0, NOUTLIERS=noutliers, IMETH=3,&
IOUTLIERSTATS=outlierstat, RSE=rse, MXLEAD=mxlead,&
OUTLIERFCST=outlierfcst)
call umach (2, nout)
write (nout,*) 'Method 3: Specified ARIMA model'
write (nout,FMT="(T2,4(A,I2))") 'Model: p =',model(1),', q =',&
model(2), ', s =', model(3), ', d =', model(4)
write (nout,*)
write (nout,FMT="(T2,A,I2)") 'Number of outliers: ', noutliers
write (nout,*) 'Outlier statistics:'
write (nout,*) 'Time point Outlier type'
do i=1,noutliers
write (nout,FMT="(I11,T15,I12)") outlierstat(i,1),&
outlierstat(i,2)
end do
write (nout,*)
write (nout,*) 'AIC=', aic
write (nout,*) 'RSE=', rse
write (nout,*)
write (nout,*) 'ARMA parameters:'
do i=1, 1+model(1)+model(2)
write (nout,FMT="(T3,I2,TR5,f10.6)") i, parameters(i)
end do
forecast_table(1:mxlead,1) = x(nobs+1:mxlead)
forecast_table(1:mxlead,2:4) = outlierfcst(1:mxlead,1:3)
write (nout,*)
fmt = '(F11.4)'
call wrrrl ('* * * Forecast Table * * *', forecast_table,&
rlabel, clabel, FMT = fmt)
end
Method 3: Specified ARIMA model
Model: p = 3, q = 2, s = 1, d = 0
Number of outliers: 1
Outlier statistics:
Time point Outlier type
109 0
AIC= 408.0768
RSE= 0.4124085
ARMA parameters:
1 0.509427
2 1.944686
3 -1.901132
4 0.901670
5 1.113016
6 -0.915008
* * * Forecast Table * * *
Series Forecast Prob. Limits Psi Weights
1 8.7000 9.1109 0.8083 0.8317
2 8.6000 9.1811 1.0513 0.6312
3 9.3000 9.5185 1.1686 0.5481
4 9.1000 9.7804 1.2497 0.6157
5 8.8000 9.7117 1.3452 0.7245
6 8.5000 9.3842 1.4671 0.7326
PHONE: 713.784.3131 FAX:713.781.9260 |