maxArma¶

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

Synopsis¶

maxArma(w, p, q)

Required Arguments¶

float w[] (Input): Array of length nObs containing the time series.
int p (Input): Non-negative number of autoregressive parameters.
int q (Input): Non-negative number of moving average parameters.

Return Value¶

An array of length 1+p+q with the estimated constant, AR and MA parameters. If no value can be computed, None is returned.

Optional Arguments¶

initialEstimates, float initAr[], float initMa[] (Input)

If specified, initAr is an array of length p containing preliminary estimates of the autoregressive parameters, and initMa is an array of length q containing preliminary estimates of the moving average parameters; otherwise, they are computed internally. If \(p=0\) or \(q=0\), then the corresponding arguments are ignored.

printLevel, int (Input)

Printing options:

`printLevel`	Action
0	No Printing.
1	Prints final results only.
2	Prints intermediate and final results.

Default: printLevel = 0.

maxIterations, int (Input)

Maximum number of estimation iterations.

Default: maxIterations = 300

varNoise (Output)

Estimate of innovation variance.

logLikelihood (Output)

Value of -2 × (ln(likelihood)) for the fitted model.

armaInfo (Output)

A structure that contains information necessary in the call to armaForecast.

meanEstimate, float (Input/Output)

Estimate of the mean of the time series w. On return, meanEstimate contains an update of the mean.

Default: Time series w is centered about its sample mean.

residual (Output)

An array of length nObs containing the residuals of the requested ARMA fit.

Description¶

The function maxArma 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 Time Domain Methodology at the beginning of this chapter, the stationary time series \(W_t\) with mean \(\mu\) can be represented by the nonseasonal autoregressive moving average (ARMA) model by the following relationship:

\[\phi (B)\left(W_t - \mu\right) = \theta(B) a_t\]

where

\[t \in ZZ = \{ \cdots, -2, -1,0,1,2, \cdots \},\]

B is the backward shift operator defined by \(B^k W_t=W_{t-k}\),

\[\phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p, \phantom{...} p \geq 0,\]

and

\[\theta(B) = 1 - \theta_1 B - \theta_2 B^2 - \cdots - \theta_q B^q, \phantom{...} q \geq 0.\]

Function maxArma estimates the AR coefficients \(\phi_1,\phi_2,\cdots,\phi_p\) and the MA coefficients \(\theta_1,\theta_2,\cdots,\theta_q\) using maximum likelihood estimation.

Function maxArma checks the initial estimates for both the autoregressive and moving average coefficients to ensure that they represent a stationary and invertible series respectively.

If

\[\phi_1, \phi_2, \cdots, \phi_p\]

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

\[1 - \phi_1 z - \phi_2 z^2 - \cdots - \phi_p z^p.\]

Initial values for the AR and MA coefficients can be supplied by vectors initAr and initMa. Otherwise, estimates are computed internally by the method of moments. maxArma computes the roots of the associated polynomials. If the AR estimates represent a non-stationary series, maxArma issues a warning message and replaces initAr with initial values that are stationary. If the MA estimates represent a non-invertible series, maxArma issues a terminal error, and new initial values have to be sought.

maxArma also validates the final estimates of the AR coefficients to ensure that they too represent a stationary series. This is done to guard against the possibility that the internal log-likelihood optimizer converged to a non-stationary solution. If non-stationary estimates are encountered, maxArma issues a fatal error message.

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

\[\mathit{AIC} = \mathrm{logLikelihood} + 2(p+q)\]

Function maxArma can also handle white noise processes, i.e. ARMA(0,0) Processes.

Examples¶

Example 1¶

Consider the Wolfer Sunspot data (Anderson 1971, p. 660) consisting of the number of sunspots observed each year from 1770 through 1869. In this example, maxArma is used to fit the following ARMA(2, 1) model:

\[\tilde{w}_t = \phi_1 \tilde{w}_{t-1} + \phi_2 \tilde{w}_{t-2} + a_t - \theta_1 a_{t-1}\]

with \(\tilde{w}_t :=w_t-\mu\), \(\mu\) the sample mean of the time series \(\left\{ w_t \right\}\).

For these data, maxArma calculated the following model:

\[\tilde{w}_t = 1.22 \tilde{w}_{t-1} - 0.56 \tilde{w}_{t-2} + a_t + 0.38 a_{t-1}\]

Defining the overall constant \(\phi_0\) by \(\phi_0 :=\mu\left( 1-\textstyle\sum_{i=1}^{p} \phi_i \right)\), we obtain the following equivalent representations:

\[w_t = \phi_0 + \phi_1 w_{t-1} + \phi_2 w_{t-2} + a_t - \theta_1 a_{t-1},\]

and

\[w_t = 15.76 + 1.22 w_{t-1} - 0.56 w_{t-2} + a_t + 0.38 a_{t-1}.\]

from __future__ import print_function
from numpy import *
from pyimsl.stat.dataSets import dataSets
from pyimsl.stat.maxArma import maxArma

n_obs = 100
p = 2
q = 1
w = empty(100)
avar = []
log_likeli = []

# get wolfer sunspot data
z = dataSets(2)
for i in range(0, n_obs):
    w[i] = z[21 + i][1]

parameters = maxArma(w, p, q,
                     maxIterations=12000,
                     varNoise=avar,
                     logLikelihood=log_likeli)

print("AR estimates are %11.4f and %11.4f" % (parameters[1], parameters[2]))
print("MA estimate is %11.4f" % parameters[3])
print("Constant estimate is %11.4f" % parameters[0])
print("-2*ln(Maximum Log Likelihood) = %11.4f" % log_likeli[0])
print("White noise variance = %11.4f" % avar[0])

Output¶

AR estimates are      1.2250 and     -0.5605
MA estimate is     -0.3830
Constant estimate is     15.7609
-2*ln(Maximum Log Likelihood) =    539.5838
White noise variance =    214.5088

Example 2¶

This example is the same as Example 1, but now initial values for the AR and MA parameters are explicitly given.

from __future__ import print_function
from numpy import *
from pyimsl.stat.dataSets import dataSets
from pyimsl.stat.maxArma import maxArma

n_obs = 100
p = 2
q = 1
w = empty(100)
avar = []
log_likeli = []
init_ar = (1.244e0, -0.575e0)
init_ma = (-0.1241e0)

# get wolfer sunspot data
z = dataSets(2)
for i in range(0, n_obs):
    w[i] = z[21 + i][1]

parameters = maxArma(w, p, q,
                     maxIterations=12000,
                     varNoise=avar,
                     logLikelihood=log_likeli,
                     initialEstimates={'initAr': init_ar, 'initMa': init_ma})

print("AR estimates are %11.4f and %11.4f" % (parameters[1], parameters[2]))
print("MA estimate is %11.4f" % parameters[3])
print("Constant estimate is %11.4f" % parameters[0])
print("-2*ln(Maximum Log Likelihood) = %11.4f" % log_likeli[0])
print("White noise variance = %11.4f" % avar[0])

Output¶

AR estimates are      1.2250 and     -0.5605
MA estimate is     -0.3830
Constant estimate is     15.7609
-2*ln(Maximum Log Likelihood) =    539.5838
White noise variance =    214.5088