NAG Library Routine Document
G13BEF
1 Purpose
G13BEF fits a multi-input model relating one output series to the input series with a choice of three different estimation criteria: nonlinear least squares, exact likelihood and marginal likelihood. When no input series are present, G13BEF fits a univariate ARIMA model.
2 Specification
SUBROUTINE G13BEF ( |
MR, NSER, MT, PARA, NPARA, KFC, NXXY, XXY, LDXXY, KEF, NIT, KZSP, ZSP, ITC, SD, CM, LDCM, S, D, NDF, KZEF, RES, STTF, ISTTF, NSTTF, WA, IWA, MWA, IMWA, KPRIV, IFAIL) |
INTEGER |
MR(7), NSER, MT(4,NSER), NPARA, KFC, NXXY, LDXXY, KEF, NIT, KZSP, ITC, LDCM, NDF, KZEF, ISTTF, NSTTF, IWA, MWA(IMWA), IMWA, KPRIV, IFAIL |
REAL (KIND=nag_wp) |
PARA(NPARA), XXY(LDXXY,NSER), ZSP(4), SD(NPARA), CM(LDCM,NPARA), S, D, RES(NXXY), STTF(ISTTF), WA(IWA) |
|
3 Description
The output series , for , is assumed to be the sum of (unobserved) components which are due respectively to the inputs , for .
Thus where is the error, or output noise component.
A typical component
may be either
(a) |
a simple regression component, (here is called a simple input), or |
(b) |
a transfer function model component which allows for the effect of lagged values of the variable, related to by
|
The noise
is assumed to follow a (possibly seasonal) ARIMA model, i.e., may be represented in terms of an uncorrelated series,
, by the hierarchy of equations
(i) |
|
(ii) |
|
(iii) |
|
as outlined in
Section 3 in G13AEF.
Note: the orders appearing in each of the transfer function models and the ARIMA model are not necessarily the same; is the result of applying non-seasonal differencing of order and seasonal differencing of seasonality and order to the series : the differenced series is then of length ; the constant term parameter may optionally be held fixed at its initial value (usually, but not necessarily zero) rather than being estimated.
For the purpose of defining an estimation criterion it is assumed that the series is a sequence of independent Normal variates having mean and variance . An allowance has to be made for the effects of unobserved data prior to the observation period. For the noise component an allowance is always made using a form of backforecasting.
For each transfer function input, you have to decide what values are to be assumed for the pre-period terms and which are in theory necessary to re-create the component series , during the estimation procedure.
The first choice is to assume that all these values are zero. In this case, in order to avoid undesirable transient distortion of the early values , you are advised first to correct the input series by subtracting from all the terms a suitable constant to make the early values , close to zero. The series mean is one possibility, but for a series with strong trend the constant might be simply .
The second choice is to treat the unknown pre-period terms as nuisance parameters and estimate them along with the other parameters. This choice should be used with caution. For example, if and , it is equivalent to fitting to the data a decaying geometric curve of the form , for , along with the other inputs, this being the form of the transient. If the output contains a strong trend of this form, which is not otherwise represented in the model, it will have a tendency to influence the estimate of away from the value appropriate to the transfer function model.
In most applications the first choice should be adequate, with the option possibly being used as a refinement at the end of the modelling process. The number of nuisance parameters is then , with a corresponding loss of degrees of freedom in the residuals. If you align the input with the output by using in its place the shifted series , then setting in the transfer function model, there is some improvement in efficiency. On some occasions when the model contains two or more inputs, each with estimation of pre-period nuisance parameters, these parameters may be co-linear and lead to failure of the routine. The option must then be ‘switched off’ for one or more inputs.
3.2 The Estimation Criterion
This is a measure of how well a proposed set of parameters in the transfer function and noise ARIMA models matches the data. The estimation routine searches for parameter values which minimize this criterion. For a proposed set of parameter values it is derived by calculating
(i) |
the components as the responses to the input series using the equations (a) or (b) above, |
(ii) |
the discrepancy between the output and the sum of these components, as the noise
|
(iii) |
the residual series from by reversing the recursive equations (i), (ii) and (iii) above. |
This last step again requires treatment of the effect of unknown pre-period values of
and other terms in the equations regenerating
.
This is identical to the treatment given in
Section 3 in G13AEF, and leads to a criterion which is a sum of squares function
, of the residuals
.
It may be shown that the finite algorithm presented there is equivalent to taking the infinite set of past values
, as (linear) nuisance parameters. The pre-period nuisance parameters for the input series are included in the reduction of
, as is the constant if it is estimated.
The covariance matrix of the vector of model parameter estimates is given by
where
is the linearized least squares matrix taken from the final iteration of the algorithm of Marquardt. From this expression are derived the vector of standard deviations, and the correlation matrix of parameter estimates. These are approximations which are only valid asymptotically, and must be treated with great caution when the parameter estimates are close to their constraint boundaries.
The residual series is available upon completion of the iterations over the range corresponding to the differenced noise series .
Because of the algorithm used for backforecasting, these are only true residuals for , provided this is positive. Estimation of pre-period terms for the inputs will also tend to reduce the magnitude of the early residuals, sometimes severely.
The model component series and may optionally be returned in place of the supplied series values, in order to assess the effects of the various inputs on the output.
3.3 Forecasting Information
For the purpose of constructing forecasts of the output series at future time points
using
G13BHF, it is not necessary to use the whole set of observations
and
, for
. It is sufficient to retain a limited set of quantities constituting the ‘state set’ as follows: for each series which appears with lagged subscripts in equations
(a),
(b),
(i),
(ii) and
(iii) above, include the values at times
for
up to the maximum lag associated with that series in the equations. Note that
(i) implicitly includes past values of
and intermediate differences of
such as
.
If later observations of the series become available, it is possible to update the state set (without re-estimating the model) using
G13BGF. If time series data is supplied with a previously estimated model, it is possible to construct the state set (and forecasts) using
G13BJF.
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: MR() – INTEGER arrayInput
On entry: the orders vector
of the ARIMA model for the output noise component.
, , and refer respectively to the number of autoregressive , moving average , seasonal autoregressive and seasonal moving average parameters.
, and refer respectively to the order of non-seasonal differencing, the order of seasonal differencing and the seasonal period.
Constraints:
- , , , , , , ;
- ;
- ;
- if , ;
- if , ;
- ;
- .
- 2: NSER – INTEGERInput
On entry: the total number of input and output series. There may be any number of input series (including none), but always one output series.
Constraints:
- ;
- if there are no parameters in the model (that is, and ), .
- 3: MT(,NSER) – INTEGER arrayInput
On entry: the transfer function model orders
,
and
of each of the input series. The order parameters for input series
are held in column
. Row
holds the value
, row 2 holds the value
and row 3 holds the value
. For a simple input,
.
Row 4 holds the value , where for a simple input, for a transfer function input for which no allowance is to be made for pre-observation period effects, and for a transfer function input for which pre-observation period effects will be treated by estimation of appropriate nuisance parameters.
When , any nonzero contents of rows 1, 2, and 3 of column are ignored.
Constraint:
, for .
- 4: PARA(NPARA) – REAL (KIND=nag_wp) arrayInput/Output
On entry: initial values of the multi-input model parameters. These are in order, firstly the ARIMA model parameters:
values of
parameters,
values of
parameters,
values of
parameters and
values of
parameters. These are followed by initial values of the transfer function model parameters
,
for the first of any input series and similarly for each subsequent input series. The final component of
PARA is the initial value of the constant
, whether it is fixed or is to be estimated.
On exit: the latest values of the estimates of these parameters.
- 5: NPARA – INTEGERInput
On entry: the exact number of , and parameters.
Constraint:
, the summation being over all the
supplied in
MT.
must be included, whether fixed or estimated.
- 6: KFC – INTEGERInput
On entry: must be set to if the constant is to remain fixed at its initial value, and if it is to be estimated.
Constraint:
or .
- 7: NXXY – INTEGERInput
On entry: the (common) length of the original, undifferenced input and output time series.
- 8: XXY(LDXXY,NSER) – REAL (KIND=nag_wp) arrayInput/Output
On entry: the columns of
XXY must contain the
NXXY original, undifferenced values of each of the input series and the output series
in that order.
On exit: if
,
XXY remains unchanged on exit.
If
, the columns of
XXY hold the corresponding values of the input component series
in place of
and the output noise component
in place of
, in that order.
- 9: LDXXY – INTEGERInput
On entry: the first dimension of the array
XXY as declared in the (sub)program from which G13BEF is called.
Constraint:
.
- 10: KEF – INTEGERInput
On entry: indicates the likelihood option.
- Gives least squares.
- Gives exact likelihood.
- Gives marginal likelihood.
Constraint:
, or .
- 11: NIT – INTEGERInput
On entry: the maximum required number of iterations.
- No change is made to any of the model parameters in array PARA except that the constant (if ) and any relating to simple input series are estimated. (Apart from these, estimates are always derived for the nuisance parameters relating to any backforecasts and any pre-observation period effects for transfer function inputs.)
Constraint:
.
- 12: KZSP – INTEGERInput
On entry: must be set to
if the routine is to use the input values of
ZSP in the minimization procedure, and to any other value if the default values of
ZSP are to be used.
- 13: ZSP() – REAL (KIND=nag_wp) arrayInput/Output
On entry: if
, then
ZSP must contain the four values used to control the strategy of the search procedure.
- Contains , the value used to constrain the magnitude of the search procedure steps.
- Contains , the multiplier which regulates the value of .
- Contains , the value of the stationarity and invertibility test tolerance factor.
- Contains , the value of the convergence criterion.
If
before entry, default values of
ZSP are supplied by the routine. These are
,
,
and
, respectively.
On exit: contains the values, default or otherwise, used by the routine.
Constraint:
if , , , , .
- 14: ITC – INTEGEROutput
On exit: the number of iterations carried out.
- Indicates that the only estimates obtained up to this point have been for the nuisance parameters relating to backforecasts, unless the marginal likelihood option is used, in which case estimates have also been obtained for simple input coefficients and for the constant (if ). This value of ITC usually indicates a failure in a consequent step of estimating transfer function input pre-observation period nuisance parameters.
- Indicates that estimates have been obtained up to this point for the constant (if ), for simple input coefficients and for the nuisance parameters relating to the backforecasts and to transfer function input pre-observation period effects.
- 15: SD(NPARA) – REAL (KIND=nag_wp) arrayOutput
On exit: the
NPARA values of the standard deviations corresponding to each of the parameters in
PARA. When the constant is fixed its standard deviation is returned as zero. When the values of
PARA are valid, the values of
SD are usually also valid. However, if an exit value of
,
or
, then the contents of
SD will be indeterminate.
- 16: CM(LDCM,NPARA) – REAL (KIND=nag_wp) arrayOutput
On exit: the first
NPARA rows and columns of
CM contain the correlation coefficients relating to each pair of parameters in
PARA. All coefficients relating to the constant will be zero if the constant is fixed. The contents of
CM will be indeterminate under the same conditions as
SD.
- 17: LDCM – INTEGERInput
On entry: the first dimension of the array
CM as declared in the (sub)program from which G13BEF is called.
Constraint:
.
- 18: S – REAL (KIND=nag_wp)Output
On exit: the residual sum of squares, , at the latest set of valid parameter estimates.
- 19: D – REAL (KIND=nag_wp)Output
On exit: the objective function, , at the latest set of valid parameter estimates.
- 20: NDF – INTEGEROutput
On exit: the number of degrees of freedom associated with .
- 21: KZEF – INTEGERInput
On entry: must not be set to
, if the values of the input component series
and the values of the output noise component
are to overwrite the contents of
XXY on exit, and must be set to
if
XXY is to remain unchanged.
- 22: RES(NXXY) – REAL (KIND=nag_wp) arrayOutput
On exit: the values of the residuals relating to the differenced values of the output series. The remainder of the first
NXXY terms in the array will be zero.
- 23: STTF(ISTTF) – REAL (KIND=nag_wp) arrayOutput
On exit: the
NSTTF values of the state set array.
- 24: ISTTF – INTEGERInput
On entry: the dimension of the array
STTF as declared in the (sub)program from which G13BEF is called.
Constraint:
, where over all input series for which .
- 25: NSTTF – INTEGEROutput
On exit: the number of values in the state set array
STTF.
- 26: WA(IWA) – REAL (KIND=nag_wp) arrayWorkspace
- 27: IWA – INTEGERInput
On entry: the dimension of the array
WA as declared in the (sub)program from which G13BEF is called.
It is not practical to outline a method for deriving the exact minimum permissible value of
IWA, but the following gives a reasonably good conservative approximation. (It should be noted that if
IWA is too small (but not grossly so) then the exact minimum is returned in
and is also printed if
.)
Let and where the orders of the output noise model are , , , , , , .
Let there be input series, where .
Let
where the transfer function model orders for input
are given by
,
,
,
.
Let .
Let and .
Finally, let , and then increment by every time any of the following conditions is satisfied. (The last six conditions should be applied separately to each input series, so that, for example, if we have two input series and if and then is incremented by .)
Then .
- 28: MWA(IMWA) – INTEGER arrayWorkspace
- 29: IMWA – INTEGERInput
On entry: the dimension of the array
MWA as declared in the (sub)program from which G13BEF is called.
Constraint:
, where the derivation of
is shown under
IWA.
If
IMWA is too small then the exact minimum needed is returned in
IMWA and if
it is also printed.
- 30: KPRIV – INTEGERInput
On entry: must not be set to
, if it is required to monitor the course of the optimization or to print out the requisite minimum values of
IWA or
IMWA in the event of an error of the type
or
. The course of the optimization is monitored by printing out at each iteration the iteration count (
ITC), the residual sum of squares (
S), the objective function (
D) and a description and value for each of the parameters in the
PARA array. The descriptions are PHI for
, THETA for
, SPHI for
, STHETA for
, OMEGA/SI for
in a simple input, OMEGA for
in a transfer function input, DELTA for
and CONSTANT for
. In addition SERIES 1, SERIES 2, etc. indicate the input series relevant to the OMEGA and DELTA parameters.
KPRIV must be set to
if the print-out of the above information is not required.
- 31: IFAIL – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. 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
is recommended. If the output of error messages is undesirable, then the value
is recommended. Otherwise, because for this routine the values of the output parameters may be useful even if
on exit, the recommended value is
.
When the value is used it is essential to test the value of IFAIL on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Note: G13BEF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
On entry, | , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | , |
or | and there are no parameters in the model ( and ). |
On entry, there is inconsistency between
NPARA and
KFC on the one hand and the orders in arrays
MR and
MT on the other, or one of the
, stored in
,
or
.
-
On entry or during execution, one or more sets of parameters do not satisfy the stationarity or invertibility test conditions.
On entry, | when , , |
or | , |
or | , |
or | , |
or | . |
On entry,
IWA is too small by a considerable margin. No information is supplied about the requisite minimum size.
On entry,
IWA is too small, but the requisite minimum size is returned in
, which is printed if
.
On entry,
IMWA is too small, but the requisite minimum size is returned in
, which is printed if
.
This indicates a failure in
F04ASF which is used to solve the equations giving the latest estimates of the parameters.
This indicates a failure in the inversion of the second derivative matrix. This is needed in the calculation of the correlation matrix and the standard deviations of the parameter estimates.
On entry or during execution, one or more sets of the ARIMA (, , or ) parameters do not satisfy the stationarity or invertibility test conditions.
On entry,
ISTTF is too small. The state set information will not be produced and if
array
XXY will remain unchanged. All other parameters will be produced correctly.
The routine has failed to converge after
NIT iterations. If steady decreases in the objective function,
, were monitored up to the point where this exit occurred, then the exit probably occurred because
NIT was set too small, so the calculations should be restarted from the final point held in
PARA.
On entry,
ISTTF is too small (see
) and
NIT iterations were carried out without the convergence conditions being satisfied (see
).
7 Accuracy
The computation used is believed to be stable.
The time taken by G13BEF is approximately proportional to .
9 Example
After the full
iterations, the following are computed and printed out: the final values of the
PARA parameters and their standard errors, the correlation matrix, the residuals for the
differenced values, the values of
and
, the values of the state set and the number of degrees of freedom.
9.1 Program Text
Program Text (g13befe.f90)
9.2 Program Data
Program Data (g13befe.d)
9.3 Program Results
Program Results (g13befe.r)