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
Published date: 03/19/2020
Last modified date: 03/19/2020