MAX_ARMA

Exact maximum likelihood estimation of the parameters in a univariate ARMA (auto‑regressive, moving average) time series model.

Required Arguments

W — Vector of length NOBS containing the stationary time series. (Input)

PAR — Vector of length NPAR. On input PAR contains initial estimates for the autoregressive parameters. On output these are replaced by the exact maximum likelihood estimates for the autoregressive parameters. (Input/Output)

PMA — Vector of length NPMA. On input PMA contains initial estimates for the moving average parameters. On output these are replaced by the exact maximum likelihood estimates for the moving average parameters. (Input/Output)

Optional Arguments

NOBS — Number of values in the time series. (Input)
Default: NOBS = size(W,1).

NPAR — Number of autoregressive parameters. (Input)
Default: NPAR =  size(PAR,1).

NPMA — Number of moving average parameters. (Input)
Default: NPMA =  size(PMA,1).

WMEAN — Estimate of the mean of the time series W. (Input)
Default: WMEAN = arithmetic mean of w.

IPRINT — Printing option. (Input)

 

IPRINT

Action

0

No printing

1

Prints final results only

2

Prints intermediate and final results

Default: IPRINT = 0.

MAXIT — Maximum number of estimation iterations. (Input)
Default: MAXIT = 500.

CNST — Estimate of the constant term θ0 in the model. (Output)

AVAR — Estimate of the noise variance. (Output)

F — Value of –2*(ln(likelihood)) for fitted model. (Output)

EWS — Array of length NOBS containing the residuals of the requested ARMA fit. (Output)

FORTRAN 90 Interface

Generic: CALL MAX_ARMA (W, PAR, PMA [, …])

Specific: The specific interface names are S_MAX_ARMA and D_MAX_ARMA.

Description

Routine MAX_ARMA is derived from the maximum likelihood estimation algorithm described by Akaike, Kitagawa, Arahata and Tada (1979), and the XSARMA routine published in the TIMSAC‑78 Library.

Using the notation developed in the introduction to this chapter, the stationary time series Wt with mean μ can be represented by the nonseasonal autoregressive moving average (ARMA) model by the following relationship:

ɸ(B)(Wt − μ) = θ(B)At

where

 

B is the backward shift operator defined by

 

 

and

 

MAX_ARMA estimates the coefficients

 

and

 

using maximum likelihood estimation.

MAX_ARMA checks the initial estimates for the autoregressive coefficients to ensure that they represent a stationary series. If

 

are the initial estimates for a stationary series then all (complex) roots of the following polynomial will fall outside the unit circle:

 

MAX_ARMA computes the roots of this polynomial for the initial estimates supplied in the vector PAR. If these estimates represent a non‑stationary series, MAX_ARMA issues a warning message and replaces PAR with initial values that are stationary.

Initial estimates can be obtained using NSPE and NSLSE procedures for calculating the autoregressive and moving average parameters of a series.

MAX_ARMA also validates its final estimates to ensure that they too represent a stationary series. This is done to guard against the possibility that MAX_ARMA converged to a non‑stationary solution. If non‑stationary estimates are encountered, MAX_ARMA quits and issues a fatal message. Routines IERCD and N1RTY (see the Reference Material section of this manual) can be used to verify that the stationary condition was met.

The ARMA process

ɸ(B)(Wt − μ) = θ(B)At

can equivalently be written in the form

 

where the constant term θ0 is defined by .

MAX_ARMA estimates μ always by the sample mean of the series.

For model selection, the ARMA model with the minimum value for AIC might be preferred

AIC = F+2p

where p = NPAR + NPMA.

Comments

Informational errors

 

Type

Code

Description

3

1

Input values for autoregressive coefficients are invalid. They do not represent a stationary time series. New values have been generated.

4

1

Maximum number of iterations exceeded. Try increasing MAXIT or use double precision.

4

2

Estimation process converged to a non‑stationary solution.

Example

Consider the Wolfer Sunspot Data (Box and Jenkins, 1976, page 530) consisting of the number of sunspots observed each year from 1770 through 1869. In this example, MAX_ARMA is used to fit the following ARMA model:

 

For these data, MAX_ARMA calculated the following estimates:

 

Letting , we can obtain the following equivalent representations:

 

 

USE MAX_ARMA_INT

USE GDATA_INT

USE NSPE_INT

IMPLICIT NONE

! SPECIFICATIONS FOR LOCAL VARIABLES

INTEGER I

REAL(KIND(1E0)) PAR(2), PMA(1), AVAR, F

REAL(KIND(1E0)) X(176,2)

REAL(KIND(1E0)) CONST

INTEGER NCOL, NROW

! Get Wolfer Sunspot Data

CALL GDATA(2,X,NROW,NCOL)

! Get preliminary PAR and PMA estimates

CALL NSPE(X(22:,2),CONST, PAR, PMA, AVAR, NOBS=100)

! TEST #1: DOCUMENT EXAMPLE

CALL MAX_ARMA(x(22:,2), PAR, PMA, nobs=100, MAXIT=12000, &

AVAR=AVAR, F=F)

WRITE(*,99994) SIZE(PAR)

WRITE (*,99996) (PAR(I),I=1,SIZE(PAR))

WRITE(*,99995) SIZE(PMA)

WRITE(*,99996) (PMA(I),I=1,SIZE(PMA))

 

WRITE(*,*) "-2*LN(MAXIMUM LOG LIKELIHOOD) = ", F

WRITE(*,*) "WHITE NOISE VARIANCE = ", AVAR

99994 FORMAT(//1H ,5('-'),2X,'FINAL PAR(I)',2X,'NPAR=',I3,2X,5('-'))

99995 FORMAT(//1H ,5('-'),2X,'FINAL PMA(I)',2X,'NPMA =',I3,2X,5('-'))

99996 FORMAT(1H ,5E20.10,/(1H ,5E20.10))

 

END

Output

----- FINAL PAR(I) NPAR= 2 -----

0.1224243164E+01 -0.5600821972E+00

 

 

----- FINAL PMA(I) NPMA = 1 -----

-0.3847315013E+00

-2*LN(MAXIMUM LOG LIKELIHOOD) = 539.5841

WHITE NOISE VARIANCE = 214.50406