AUTO_ARIMA
Automatically identifies time series outliers, determines parameters of a multiplicative seasonal ARIMA (p,0,q) × (0,d,0)s model and produces forecasts that incorporate the effects of outliers whose effects persist beyond the end of the series.
Required Arguments
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 (p,0,q) × (0,d,0)s model. If d = 0, then an ARMA(p, q) model is fitted to the outlier-free version of the observed series . If d > 0, these parameters are computed for an ARMA(p,q) representation of the seasonally adjusted series , where and s ≥ 1…. (Output)
Optional Arguments
NOBS — Number of observations in the time series. Assuming that the series is defined at time points t1, …, tNOBS, the actual length of the series, including missing values, is . (Input)
Default: NOBS = size (W).
Default: NOBS = size (W).
IMETH — Method to be used in model selection. (Input)
1 –Automatic selection
1 –Automatic selection
2 – Grid search (requires arguments IAR and IMA)
3 – Specified model (Requires argument MODEL)
Default: IMETH = 1.
Default: IMETH = 1.
MAXLAG — Maximum number of autoregressive parameters allowed in fitting an AR model to the series. (Input)
Default: MAXLAG = 10.
Default: MAXLAG = 10.
INFOCRIT — The information criterion used for optimum model selection. (Input)
INFOCRIT | selected information criterion |
|---|---|
0 | Akaike’s Information Criterion (AIC) |
1 | Akaike’s 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.
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.
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.
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).
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
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
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.
Default: ALPHA = 0.05.
MXLEAD — Maximum lead time for forecasts. (Input)
Default: MXLEAD = 0.
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.
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 | 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, 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)
FORTRAN 90 Interface
Generic: CALL AUTO_ARIMA (ITIME_POINTS, W, PARAMS [, …])
Specific: The specific interface names are S_AUTO_ARIMA and D_AUTO_ARIMA.
Description
Routine AUTO_ARIMA determines the parameters of a multiplicative seasonal ARIMA (p,0,q) × (0,d,0)s 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 (p,0,q) × (0,d,0)s model handled by AUTO_ARIMA has the following form:
where
and
It is assumed that all roots of φ(B) and θ(B) lie outside the unit circle. Clearly, if s = 1 this reduces to the traditional ARIMA(p, d, q) model.
Yt is the unobserved, gap-free and outlier-free time series with mean μ, and white noise at. 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 Yt, the outlier-free series representation of the data, by removing the estimated outlier effects.
Using the information about the adjusted ARIMA (p,0,q) × (0,d,0)s 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.
Model Selection
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.
Method 1: Automatic ARIMA (p,0,0) x (0,d,0)s Selection
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 (p,0,0) × (0,d,0)s 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 (p,0,0) × (0,d,0)s 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.
Method 2: Grid Search
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.
Method 3: Specified ARIMA (p,0,q) x (0,d,0)s Model
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 s > 0 and d ≥ 0 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.
Outliers
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 (p,0,q) × (0,d,0)s model. |
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 | 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 | 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 outlier’s 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.
Examples
Example 1
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
Output
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
Example 2
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 p = 3, q = 2, s = 1, d = 0 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, dime nsion(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
Output
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
Example 3
This example is the same as Example 1 but now Method 3 with the optimum model parameters p = 3, q = 2, s = 1, d = 0 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
Output
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
