G13AFF (PDF version)
G13 Chapter Contents
G13 Chapter Introduction
NAG Library Manual

NAG Library Routine Document


Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.


    1  Purpose
    7  Accuracy

1  Purpose

G13AFF is an easy-to-use version of G13AEF. It fits a seasonal autoregressive integrated moving average (ARIMA) model to an observed time series, using a nonlinear least squares procedure incorporating backforecasting. Parameter estimates are obtained, together with appropriate standard errors. The residual series is returned, and information for use in forecasting the time series is produced for use in G13AGF and G13AHF.
The estimation procedure is iterative, starting with initial parameter values such as may be obtained using G13ADF. It continues until a specified convergence criterion is satisfied or until a specified number of iterations have been carried out. The progress of the iteration can be monitored by means of an optional printing facility.

2  Specification


3  Description

The time series x1,x2,,xn supplied to the routine is assumed to follow a seasonal autoregressive integrated moving average (ARIMA) model defined as follows:
where dsDxt is the result of applying non-seasonal differencing of order d and seasonal differencing of seasonality s and order D to the series xt, as outlined in the description of G13AAF. The differenced series is then of length N=n-d, where d=d+D×s is the generalized order of differencing. The scalar c is the expected value of the differenced series, and the series w1,w2,,wN follows a zero-mean stationary autoregressive moving average (ARMA) model defined by a pair of recurrence equations. These express wt in terms of an uncorrelated series at, via an intermediate series et. The first equation describes the seasonal structure:
The second equation describes the non-seasonal structure. If the model is purely non-seasonal the first equation is redundant and et above is equated with wt:
Estimates of the model parameters defined by
ϕ1,ϕ2,,ϕp,θ1,θ2,,θq, Φ1,Φ2,,ΦP,Θ1,Θ2,,ΘQ  
and (optionally) c are obtained by minimizing a quadratic form in the vector w=w1,w2,,wN.
The minimization process is iterative, iterations being performed until convergence is achieved (see Section 3 in G13AEF for full details), or until the user-specified maximum number of iterations are completed.
The final values of the residual sum of squares and the parameter estimates are used to obtain asymptotic approximations to the standard deviations of the parameters, and the correlation matrix for the parameters. The ‘state set’ array of information required by forecasting is also returned.
Note:  if the maximum number of iterations are performed without convergence, these quantities may not be reliable. In this case, the sequence of iterates should be checked, using the optional monitoring routine, to verify that convergence is adequate for practical purposes.

4  References

Box G E P and Jenkins G M (1976) Time Series Analysis: Forecasting and Control (Revised Edition) Holden–Day
Marquardt D W (1963) An algorithm for least squares estimation of nonlinear parameters J. Soc. Indust. Appl. Math. 11 431

5  Parameters

1:     MR7 – INTEGER arrayInput
On entry: the orders vector p,d,q,P,D,Q,s of the ARIMA model whose parameters are to be estimated. p, q, P and Q refer respectively to the number of autoregressive ϕ, moving average θ, seasonal autoregressive Φ and seasonal moving average Θ parameters. d, D and s refer respectively to the order of non-seasonal differencing, the order of seasonal differencing and the seasonal period.
  • p, d, q, P, D, Q, s0;
  • p+q+P+Q>0;
  • s1;
  • if s=0, P+D+Q=0;
  • if s>1, P+D+Q>0;
  • d+s×P+Dn;
  • p+d-q+s×P+D-Qn.
2:     PARNPAR – REAL (KIND=nag_wp) arrayInput/Output
On entry: the initial estimates of the p values of the ϕ parameters, the q values of the θ parameters, the P values of the Φ parameters and the Q values of the Θ parameters, in that order.
On exit: contains the latest values of the estimates of these parameters.
3:     NPAR – INTEGERInput
On entry: the total number of ϕ, θ, Φ and Θ parameters to be estimated.
Constraint: NPAR=p+q+P+Q.
4:     C – REAL (KIND=nag_wp)Input/Output
On entry: if KFC=0, C must contain the expected value, c, of the differenced series.
If KFC=1, C must contain an initial estimate of c.
Therefore, if C and KFC are both zero on entry, there is no constant correction.
On exit: if KFC=0, C is unchanged.
If KFC=1, C contains the latest estimate of c.
5:     KFC – INTEGERInput
On entry: must be set to 1 if the constant, c, is to be estimated and 0 if it is to be held fixed at its initial value.
Constraint: KFC=0 or 1.
6:     XNX – REAL (KIND=nag_wp) arrayInput
On entry: the n values of the original undifferenced time series.
7:     NX – INTEGERInput
On entry: n, the length of the original undifferenced time series.
8:     S – REAL (KIND=nag_wp)Output
On exit: the residual sum of squares after the latest series of parameter estimates has been incorporated into the model. If the routine exits with a faulty input parameter, S contains zero.
9:     NDF – INTEGEROutput
On exit: the number of degrees of freedom associated with S, NDF=n-d-D×s-p-q-P-Q-KFC.
10:   SDNPPC – REAL (KIND=nag_wp) arrayOutput
On exit: the standard deviations corresponding to the parameters in the model (p autoregressive, q moving average, P seasonal autoregressive, Q seasonal moving average and c, if estimated, in that order). If the routine exits with IFAIL containing a value other than 0 or 9, or if the required number of iterations is zero, the contents of SD will be indeterminate.
11:   NPPC – INTEGERInput
On entry: the number of ϕ, θ, Φ, Θ and c parameters to be estimated. NPPC=p+q+P+Q+1 if the constant is being estimated and NPPC=p+q+P+Q if not.
Constraint: NPPC=NPAR+KFC.
12:   CMLDCMNPPC – REAL (KIND=nag_wp) arrayOutput
On exit: the correlation coefficients associated with each pair of the NPPC parameters. These are held in the first NPPC rows and the first NPPC columns of CM. These correlation coefficients are indeterminate if IFAIL contains on exit a value other than 0 or 9, or if the required number of iterations is zero.
13:   LDCM – INTEGERInput
On entry: the first dimension of the array CM as declared in the (sub)program from which G13AFF is called.
Constraint: LDCMNPPC.
14:   STNX – REAL (KIND=nag_wp) arrayOutput
On exit: the value of the state set in its first NST elements. If the routine exits with IFAIL containing a value other than 0 or 9, the contents of ST will be indeterminate.
15:   NST – INTEGEROutput
On exit: the size of the state set. NST=P×s+D×s+d+q+maxp,Q×s.
NST should be used subsequently in G13AGF and G13AHF as the dimension of ST.
16:   KPIV – INTEGERInput
On entry: must be nonzero if the progress of the optimization is to be monitored using the built-in printing facility. Otherwise KPIV must contain zero. If selected, monitoring output will be sent to the current advisory message unit defined by X04ABF. For each iteration, the heading
followed by the parameter values, and residual sum of squares, are printed.
17:   NIT – INTEGERInput
On entry: the maximum number of iterations to be performed.
Constraint: NIT0.
18:   ITC – INTEGEROutput
On exit: the number of iterations performed.
19:   ISF4 – INTEGER arrayOutput
On exit: the first four elements of ISF contain success/failure indicators, one for each of the four types of parameter in the model (autoregressive, moving average, seasonal autoregressive, seasonal moving average), in that order.
Each indicator has the interpretation:
-2 On entry parameters of this type have initial estimates which do not satisfy the stationarity or invertibility test conditions.
-1 The search procedure has failed to converge because the latest set of parameter estimates of this type is invalid.
-0 No parameter of this type is in the model.
-1 Valid final estimates for parameters of this type have been obtained.
20:   RESIRES – REAL (KIND=nag_wp) arrayOutput
On exit: the first NRES elements of RES contain the model residuals derived from the differenced series. If the routine exits with IFAIL holding a value other than 0 or 9, these elements of RES will be indeterminate. The rest of the array RES is used as workspace.
21:   IRES – INTEGERInput
On entry: the dimension of the array RES as declared in the (sub)program from which G13AFF is called.
Constraint: IRES15×Q+11n+13×NPPC+8×P+12+2× Q+NPPC 2, where P=p+P×s and Q=q+Q×s.
22:   NRES – INTEGEROutput
On exit: the number of model residuals returned in RES.
23:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to 0, -1​ or ​1. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1​ or ​1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, because for this routine the values of the output parameters may be useful even if IFAIL0 on exit, the recommended value is -1. When the value -1​ or ​1 is used it is essential to test the value of IFAIL on exit.
On exit: IFAIL=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

If on entry IFAIL=0 or -1, explanatory error messages are output on the current error message unit (as defined by X04AAF).
Note: G13AFF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
On entry,NPARp+q+P+Q,
orthe orders vector MR is invalid (check the constraints in Section 5),
orKFC0 or 1,
On entry, NX-d-D×sNPAR+KFC, i.e., the number of terms in the differenced series is not greater than the number of parameters in the model. The model is over-parameterised.
On entry,NIT<0.
On entry, the required size of the state set array ST is greater than NX. This occurs only for very unusual models with long seasonal periods or large numbers of parameters. First check that the orders vector MR has been set up as intended. If it has, change to G13AEF with ST dimensioned at least (NST), where NST is the value returned by G13AFF, or computed using the formula in Section 5 of this document.
On entry, the workspace array RES is too small. Check the value of IRES against the constraints in Section 5.
On entry,LDCM<NPPC.
The search procedure in the algorithm has failed. This may be due to a badly conditioned sum of squares function, or the default convergence criterion may be too strict. Use G13AEF with a less strict convergence criterion.
Some output parameters may contain meaningful values; see Section 5 for details.
The inversion of the Hessian matrix in the calculation of the covariance matrix of the parameter estimates has failed.
Some output parameters may contain meaningful values; see Section 5 for details.
This indicates a failure when solving the equations giving the latest estimates of the backforecasts.
Some output parameters may contain meaningful values; see Section 5 for details.
Satisfactory parameter estimates could not be obtained for all parameter types in the model. Inspect array ISF for further information on the parameter type(s) in error.
An internal error has arisen in partitioning RES for use by G13AEF. This error should not occur; report it to NAG via your site representative.
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.8 in the Essential Introduction for further information.
Your licence key may have expired or may not have been installed correctly.
See Section 3.7 in the Essential Introduction for further information.
Dynamic memory allocation failed.
See Section 3.6 in the Essential Introduction for further information.

7  Accuracy

The computations are believed to be stable.

8  Parallelism and Performance

G13AFF is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
G13AFF makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9  Further Comments

The time taken by G13AFF is approximately proportional to NX×ITC× q+Q×s+NPPC 2.

10  Example

This example reads 30 observations from a time series relating to the rate of the earth's rotation about its polar axis. Differencing of order 1 is applied, and the number of non-seasonal parameters is 3, one autoregressive (ϕ) and two moving average (θ). No seasonal effects are taken into account.
The constant is estimated. Up to 50 iterations are allowed.
The initial estimates of ϕ1, θ1, θ2 and c are zero.
Some intermediate monitoring output from G13AFZ has been omitted.

10.1  Program Text

Program Text (g13affe.f90)

10.2  Program Data

Program Data (g13affe.d)

10.3  Program Results

Program Results (g13affe.r)

G13AFF (PDF version)
G13 Chapter Contents
G13 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2015