Welcome to the Harris Geospatial product documentation center. Here you will find reference guides, help documents, and product libraries.
﻿

>  Docs Center  >  IDL Reference  >  Advanced Math and Stats  >  IMSL_ARMA

### IMSL_ARMA

IMSL_ARMA

The IMSL_ARMA function computes method-of-moments or least-squares estimates of parameters for a nonseasonal IMSL_ARMA model.

## Syntax

Result = IMSL_ARMA(z, p, q [, AR_LAGS=array] [, AUTOCOV=variable] [, BACKWARD_ORIGIN=value] [, CONFIDENCE=value] [, CONSTANT]
[, /DOUBLE] [, ERR_REL=value] [, FORECAST=variable] [, INIT_EST_AR=array] [, INIT_EST_MA=array] [, ITMAX=value] [, /LSQ] [, LGTH_BACKCAST=value] [, MA_LAGS=array] [, MEAN_EST=value]
[, /MOMENTS] [, N_PREDICT=value] [, /NO_CONSTANT] [, PARAM_EST_COV=variable] [, RESIDUAL=variable] [, SS_RESIDUAL=variable] [, TOL_BACKCAST=value] [, TOL_CONVERGENCE=value])

## Return Value

An array of length 1 + p + q with the estimated constant, AR, and MA parameters. If No_Constant is specified, the 0-th element of this array is 0.0.

## Arguments

#### p

Number of autoregressive parameters.

#### q

Number of moving average parameters.

#### z

One-dimensional array containing the observations.

## Keywords

#### AR_LAGS

One-dimensional array of length p containing the order of the nonzero autoregressive parameters. The elements of Ar_Lags must be greater than or equal to 1. Default: Ar_Lags = [1, 2, ..., p]

#### AUTOCOV

Named variable into which an array of length p + q + 2 containing the variance and autocovariances of the time series z is stored. Keyword Autocov(0) contains the variance of the series z. Keyword Autocov(k) contains the autocovariance of lag k, where k = 1, ..., p + q + 1.

#### BACKWARD_ORIGIN

Maximum backward origin. Keyword Backward_Origin must be greater than or equal to zero and less than or equal to N_ELEMENTS(z) – (max(maxar, maxma)), where maxar = max(Ar_Lags) and maxma = max(Ma_Lags).

Forecasts at origins N_ELEMENTS(z) – Backward_Origin through N_ELEMENTS(z) are generated. Default: Backward_Origin = 0

#### CONFIDENCE

Value in the exclusive interval (0, 100) used to specify the confidence level of the forecasts. Typical choices for Confidence are 90.0, 95.0, and 99.0. Default: Confidence = 95.0

#### CONSTANT

If present and nonzero, the time series is centered about its mean. Keywords No_Constant and Constant cannot be used together.

#### DOUBLE

If present and nonzero, double precision is used.

#### ERR_REL

Stopping criterion for use in the nonlinear equation solver used in both the method-of-moments and least-squares algorithms. Default: Err_Rel = 100 x e, where e is machine precision

#### FORECAST

Named variable into which an array of length N_Predict x (Backward_Origin + 3) containing the forecasts up to N_Predict steps ahead and the information necessary to obtain confidence intervals is stored. Keywords Forecast and N_Predict must be used together.

#### INIT_EST_AR

Array of length p containing preliminary estimates of the autoregressive parameters, internally. Keywords Init_Est_Ar and Init_Est_Ma must be used together and are only applicable if Lsq is also present and nonzero.

#### INIT_EST_MA

Array of length q containing preliminary estimates of the moving average parameters. Keywords Init_Est_Ar and Init_Est_Ma must be used together and are only applicable if Lsq is also present and nonzero.

The following keywords are used to forecast up to N_Predict steps ahead and the information necessary to obtain confidence intervals:

#### ITMAX

Maximum number of iterations allowed in the nonlinear equation solver used in both the method-of-moments and least-squares algorithms. Default: Itmax = 200

#### LSQ

If present and nonzero, the autoregressive and moving average parameters are estimated by a least-squares procedure. Keywords Moments and Lsq cannot be used together.

#### LGTH_BACKCAST

Specifies the maximum length of backcasting. Must be greater than or equal to zero. Keywords Lgth_Backcast and Tol_Backcast must be used together. Default: Lgth_Backcast = 10

#### MA_LAGS

One-dimensional array of length q containing the order of the nonzero moving average parameters. The elements of Ma_Lags must be greater than or equal to 1. Default: Ma_Lags = [1, 2, ..., q]

#### MEAN_EST

Initial estimate of the mean of the time series z.

Default:

#### MOMENTS

If present and nonzero, the autoregressive and moving average parameters are estimated by a method-of-moments procedure. Keywords Moments and Lsq cannot be used together. (Default)

#### N_PREDICT

Maximum lead time for forecasts. Keyword N_Predict must be greater than zero. Keywords Forecast and N_Predict must be used together.

#### NO_CONSTANT

If present and nonzero, the time series is not centered about its mean. Keywords No_Constant and Constant cannot be used together.

#### PARAM_EST_COV

Named variable into which an array, containing the covariance matrix of the final parameter estimates, is stored. The array is of size np x np, where np = p + q + 1 if z is centered about its mean and np = p + q if z is not centered. The ordering of variables in Param_Est_Cov is Mean_Est, Ar_lags, and Ma_lags.

#### RESIDUAL

Named variable into which an array of length N_ELEMENTS(z) – (max(Ar_Lags)) + Lgth_Backcast containing the residuals (including backcasts) at the final parameter estimate point in the first N_ELEMENTS(z) – (max(Ar_Lags)) + nb, where nb is the number of values backcast is stored.

#### SS_RESIDUAL

Named variable into which the sum of squares of the random error is stored.

#### TOL_BACKCAST

Specifies the tolerance level used to determine convergence of the backcast algorithm. Typically, Tol_Backcast is set to a fraction of an estimate of the standard deviation of the time series. Keywords Lgth_Backcast and Tol_Backcast must be used together. Default: Tol_Backcast = 0.01 x standard deviation of l

#### TOL_CONVERGENCE

Tolerance level used to determine convergence of the nonlinear least-squares algorithm. Keyword Tol_Convergence represents the minimum relative decrease in sum of squares between two iterations required to determine convergence. Hence, Tol_Convergence must be greater than or equal to zero. Default: max {10–10, e2 / 3} for single precision,
max {10–20, e2 / 3} for double precision, where e is machine precision.

## Discussion

The IMSL_ARMA function computes estimates of parameters for a nonseasonal IMSL_ARMA model given a sample of observations, {Zt}, for t = 1, 2, ..., n, where
n = N_ELEMENTS(z). You may choose either method of moments or least squares. The default is method of moments.

Choose the method-of-moments algorithm with the keyword Moments. The least-squares algorithm is used if Lsq is specified. If you wish to use the least-squares algorithm, the preliminary estimates are the method-of-moments estimates by default; otherwise, you can input initial estimates by specifying keywords Init_Est_Ar and Init_Est_Ma. Table 20-1 lists the appropriate keywords for both the method-of-moments and least-squares algorithm:

Table 20-1: Method-of-Moments and Least-Squares Keywords

Method of
Moments
only
Least Squares only
Both Method of
Moments and
Least Squares
Moments
Lsq
Err_Rel

Constant (or No_Constant)
Itmax

Ar_Lags
Mean_Estimate

Ma_Lags
Autocov

Lgth_Backcast
Forecast

Tol_Backcast
N_Predict

Tol_Convergence
Confidence

Init_Est_Ar
Backward_Origin

Init_Est_Ma

Residual

Param_Est_Cov

Ss_Residual

##### Method-of-moments Estimation

Suppose the time series {Zt } is generated by an IMSL_ARMA(p, q) model of the form:

for

Let

be the estimate of the mean m of the time series {Zt}, where:

equals the following:

The autocovariance function is estimated by:

for k = 0, 1, ..., K, where K = p + q + 1. Note that:

is an estimate of the sample variance.

Given the sample autocovariances, the function computes the method-of-moments estimates of the autoregressive parameters using the extended Yule-Walker equations as follows:

where:

The overall constant q0 is estimated by the following:

The moving average parameters are estimated based on a system of nonlinear equations given K = p + q + 1 autocovariances, s(k) for k = 1, ..., K, and p autoregressive parameters fi for i = 1, ..., p.

Let Z't = f(B)Zt. The autocovariances of the derived moving average process Z't = q(B)At are estimated by the following relation:

The iterative procedure for determining the moving average parameters is based on the relation:

where s(k) denotes the autocovariance function of the original Zt process.

Let t = (t0, t1, ..., tq)T, f = (f0, f1, ..., fq)T, and T be a (q + 1) x (q + 1) matrix, where tj , fj , and T are as follows:

and:

Then, the value of t at the (i + 1)-th iteration is determined by:

t i + 1 = t i – (T i)–1 f i

The estimation procedure begins with the initial value:

and terminates at iteration i when either |f i| is less than Err_Rel or i equals Itmax. The moving average parameter estimates are obtained from the final estimate of t by setting:

for j = 1, ..., q. The random error variance is estimated by the following:

See Box and Jenkins (1976, pp. 498–500) for a description of a function that performs similar computations.

##### Least-squares Estimation

Suppose the time series {Zt } is generated by a nonseasonal IMSL_ARMA model of the form:

where B is the backward-shift operator, m is the mean of Zt , and:

with p autoregressive and q moving average parameters. Without loss of generality, the following is assumed:

so that the nonseasonal IMSL_ARMA model is of order (p', q'), where:

and

Note that the usual hierarchial model assumes the following:

Consider the sum-of-squares function:

where:

and T = Lgth_Backcast is the length of backcasting from the beginning of the series. The random errors {At } are assumed to be independent and identical distributed N(0, sA2) random variables. Hence, the log-likelihood function is given by:

where f (m, f, q) is a function of m, f, and q.

For T = 0, the log-likelihood function is conditional on the past values of both Zt and At required to initialize the model. The method of selecting these initial values usually introduces transient bias into the model (Box and Jenkins 1976, pp. 210–211). For T = infinity, this dependency vanishes, and the estimation problem concerns maximization of the unconditional log-likelihood function. Box and Jenkins (1976, p. 213) argue that:

dominates

The parameter estimates that minimize the sum-of-squares function are called least-squares estimates. For large n, the unconditional least-squares estimates are approximately equal to the maximum likelihood-estimates.

In practice, a finite value of T enables sufficient approximation of the unconditional sum-of-squares function. The values of [At] needed to compute the unconditional sum of squares are computed iteratively with initial values of Zt obtained by backcasting. The residuals (including backcasts), estimate of random error variance, and covariance matrix of the final parameter estimates also are computed. ARIMA parameters can be computed using the IMSL_DIFFERENCE, together with IMSL_ARMA.

##### Forecasting Option

The Box-Jenkins forecasts and their associated confidence intervals for a nonseasonal IMSL_ARMA model are computed given a sample of n = N_ELEMENTS(z) {Zt} for t = 1, 2, ..., n.

Suppose the time series {Zt} is generated by a nonseasonal IMSL_ARMA model of the form:

f (B) Zt = q0 + q (B) At

for

where B is the backward-shift operator, q0 is the constant, and:

with p autoregressive and q moving average parameters. Without loss of generality, the following is assumed:

so that the nonseasonal IMSL_ARMA model is of order (p', q'), where:

and

Note that the usual hierarchal model assumes the following:

The Box-Jenkins forecast at origin t for lead time l of Zt + l is defined in terms of the difference equation:

where the following is true:

The 100(1 – a)-percent confidence interval for Zt + l is given by:

where

is the 100 (1 – a / 2)-percentile of the standard normal distribution, sA is the standard deviation of the random error, and yj is defined as follows:

In this equation, fi = 0 for i > p and qj = 0 for j > q. Note that the forecasts are computed for lead times l = 1, 2, ..., L at origins t = (nb), (nb + 1), ..., n, where L = N_Predict and b = Backward_Origin.

The Box-Jenkins forecasts minimize the mean-square error:

Also, the forecasts are easily updated according to the following equation:

This approach and others are discussed in Chapter 5 of Box and Jenkins (1976).

## Examples

#### Example 1

Consider the Wolfer Sunspot Data (Anderson 1971, p. 660) consisting of the number of sunspots observed each year from 1749 through 1924. The data set for this example consists of the number of sunspots observed from 1770 through 1869 and is shown in Figure 20-1. The method-of-moments estimates:

, and

for the IMSL_ARMA(2,1) model are:

where Zt is "raw" data and the errors At are independently and identical normally distributed with mean zero and variance s2A.

`temp = IMSL_STATDATA(2)  ; Get the Wolfer Sunspot Data.  z = TEMP(21:120, 1)  ; Use only 100 observations, 1770-1869.  years = FINDGEN(100) + 1770  PLOT, years, z, XStyle = 1, Psym   = -6, \$     Title  = ';Wolfer Sunspot Data', XTitle = 'Year', \$     YTitle = ';Number of Sunspots'  ; Plot the data.  p = 2  q = 1  parameters = IMSL_ARMA(z, p, q)  ; Perform time-series analysis.  PRINT, ';AR estimates:', parameters(1), parameters(2)  PRINT, ';MA estimate :', parameters(3)    AR estimates: 1.24426 -0.575149  MA estimate : -0.124094  `

Figure 20-1: Wolfer Sunspot Data Plot

#### Example 2

The data for this example are the same as that for the initial example. Preliminary method-of-moments estimates are computed by default, and the method of least squares is used to find the final estimates.

`temp = IMSL_STATDATA(2)  ; Get the Wolfer Sunspot Data.  z = TEMP(21:120, 1)  ; Use only 100 observations, 1770-1869.  parameters = IMSL_ARMA(z, 2, 1, /Lsq, Tol_Convergence = .125)   ; Perform time-series analysis using method of moments. The  ; warning error can be ignored in this case.  PRINT, ';AR estimates:', parameters(1), parameters(2)  PRINT, ';MA estimate :', parameters(3)    AR estimates: 1.39257 -0.732948  MA estimate : -0.137512  `

#### Example 3

Consider the Wolfer Sunspot Data (Anderson 1971, p. 660) consisting of the number of sunspots observed each year from 1749 through 1924. The data set for this example consists of the number of sunspots observed from 1770 through 1869. IMSL_ARMA computes forecasts and 95-percent confidence limits for the forecasts for an IMSL_ARMA(2, 1) model fit using IMSL_ARMA with the method-of-moments option. With Backward_Origin = 3, columns zero through three of Forecast provide forecasts given the data through 1866, 1867, 1868, and 1869. Column five gives the deviations from the forecast for computing confidence limits, and column six gives the psi weights, which can be used to update forecasts when more data is available. For example, the forecast for the 102-nd observation (year 1871) given the data through the 100-th observation (year 1869) is 77.21; 95-percent confidence limits are given by:

After observation 101 (Z101 for year 1870) is available, the forecast can be updated by using:

with the psi weight (y1 = 1.37) and the one-step-ahead forecast error for observation 101 (Z101 – 83.72) to give the following:

77.21 + 1.37 x (Z101 – 83.72)

Since this updated forecast is one step ahead, the 95-percent confidence limits are now given by the forecast:

First, define a procedure to output the results:

`.RUN  PRO print_results, parameters, forecast     PRINT, ';Method-of-moments initial estimates:'     PRINT, ';AR estimates:', parameters(1), parameters(2)     PRINT, ';MA estimate :', parameters(3)     PRINT     lead_time = INDGEN(12) + 1     forecast = [[lead_time], [forecast]]     PRINT, ';Forecasts from ...'     PRINT, ';lead time', ' 1866', ' 1867', \$        '; 1868', ' 1869', ' Deviat.', ' Psi'     PM, forecast, FORMAT = ';(i6, 3x, 6f9.4)'  END    temp = IMSL_STATDATA(2)  ; Get the Wolfer Sunspot Data.  z = TEMP(21:120, 1)  ; Use only 100 observations, 1770-1869.  parameters = IMSL_ARMA(z, 2, 1, Itmax = 0, Err_Rel = 0.0, \$     Forecast = forecast, N_Predict = 12, Backward_Origin = 3)  ; Perform time-series analysis using method-of-moments.  print_results, parameters, forecast  years = INDGEN(100) + 1770  PLOT, years, z, \$     Psym   = -6, Symsize = .5, XStyle = 1, XRange = [1770, 1885], \$     YRange = [-50, 175], Title  = ';Wolfer Sunspot Data', \$     XTitle = ';Year', YTitle = 'Number of Sunspots'  ; Plot the data along with the forecasted values with confidence  ; intervals.  OPLOT, INDGEN(10) + 1870, forecast(*, 3), Psym   = 4, Symsize = .5  ERRPLOT, indgen(10) + 1870, forecast(*, 3) - forecast(*, 4), \$     forecast(*, 3) + forecast(*, 4), Width  = .005    Method-of-moments initial estimates:     AR estimates:      1.24426    -0.575149     MA estimate :    -0.124094  Forecasts from ...     lead time  1866     1867     1868     1869     Deviat.   Psi       1     18.2833  16.6151  55.1893  83.7196  33.2179   1.3684       2     28.9182  32.0189  62.7606  77.2092  56.2980   1.1274       3     41.0101  45.8275  61.8922  63.4608  67.6168   0.6158       4     49.9387  54.1496  56.4571  50.0987  70.6432   0.1178       5     54.0937  56.5623  50.1939  41.3803  70.7515  -0.2076       6     54.1282  54.7780  45.5268  38.2174  71.0869  -0.3261       7     51.7815  51.1701  43.3221  39.2965  71.9074  -0.2863       8     48.8417  47.7072  43.2631  42.4582  72.5337  -0.1687       9     46.5335  45.4736  44.4577  45.7715  72.7498  -0.0452      10     45.3524  44.6861  45.9781  48.0758  72.7653   0.0407      11     45.2103  44.9909  47.1827  49.0371  72.7779   0.0767      12     45.7128  45.8230  47.8072  48.9080  72.8225   0.0720  `

The plot of the forecasts and the confidence limits from year 1869 are shown in Figure 20-2.

Figure 20-2: Sunspot Data with Predicted Values and Confidence Bands

## Version History

 6.4 Introduced

© 2018 Harris Geospatial Solutions, Inc. |  Legal