armaForecast¶

Computes forecasts and their associated probability limits for an ARMA model.

Synopsis¶

armaForecast (armaInfo, nPredict)

Required Arguments¶

structure armaInfo (Input): A structure that is passed from the arma function.
int nPredict (Input): Maximum lead time for forecasts. Argument nPredict must be greater than 0.

Return Value¶

An array of length nPredict × (backwardOrigin + 3) containing the forecasts up to nPredict steps ahead and the information necessary to obtain pairwise confidence intervals.

Optional Arguments¶

confidence, float (Input)

Value in the exclusive interval (0, 100) used to specify the confidence percent probability limits of the forecasts. Typical choices for confidence are 90.0, 95.0, and 99.0.

Default: confidence = 95.0.

backwardOrigin, int (Input)

If specified, the maximum backward origin. Argument backwardOrigin must be greater than or equal to 0 and less than or equal to nObservations - max(maxar, maxma), where maxar = max(arLags[i]), maxma = max (maLags[j]), and nObservations = the number of observations in the series, as input in function arma. nPredict forecasts beginning at origins nObservations - backwardOrigin +1 through nObservations are generated.

Default: backwardOrigin = 0.

Description¶

The Box-Jenkins forecasts and their associated probability limits for a nonseasonal ARMA model are computed given a sample of n = nObservations \(\{Z_t\}\) for \(t=1,2,\ldots,n\), where nObservations = the number of observations in the series, as input in function arma.

Suppose the time series \(\{Z_t\}\) is generated by a nonseasonal ARMA model of the form

\[φ(B)Z_t = θ_0 + θ(B)A_t\]

for \(t\in\{0,\pm 1,\pm 2,\ldots\}\), where B is the backward shift operator, \(\theta_0\) is the constant, and

\[\begin{split}\begin{array}{l} \phi(B) = 1 - \phi_1B^{l_\phi(1)} - \phi_2B^{l_\phi(2)} - \ldots - \phi_pB^{l_\phi(p)} \\ \theta(B) = 1 - \theta_1B^{l_\theta(1)} - \theta_2B^{l_\theta(2)} - \ldots - \theta_qB^{l_\theta(q)} \end{array}\end{split}\]

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

\[1 ≤ l_φ (1) ≤ l_φ (2) ≤ … ≤ l_φ (p)\]

\[1 ≤ l_θ (1) ≤ l_θ (2) ≤ … ≤ l_θ (q)\]

so that the nonseasonal ARMA model is of order \((p',q')\), where \(p'=l_\varphi(p)\) and \(q'=l_\theta(q)\). Note that the usual hierarchical model assumes the following:

\[l_φ (i) = i, 1 ≤ i ≤ p\]

\[l_θ (j) = j, 1 ≤ j ≤ q\]

The Box-Jenkins forecast at origin t for lead time l of \(Z_{t+l}\) is defined in terms of the difference equation

\[\begin{split}\begin{array}{l} \hat{Z}_t(l) = \theta_0 + \phi_1\left[Z_{t+l-l_\phi(1)}\right] + \ldots + \phi_p\left[Z_{t+l-l_\phi(p)}\right] \\ + \left[A_{t+l}\right] - \theta_1\left[A_{t+l-l_\theta(1)}\right] - \ldots - \theta_q\left[A_{t+l-l_\theta(q)}\right] \end{array}\end{split}\]

where the following is true:

\[\begin{split}\left[Z_{t+k}\right] = \begin{cases} Z_{t+k} & \text{for } k = 0,-1,-2, \ldots \\ \hat{Z}_t(k) & \text{for } k=1,2, \ldots \end{cases}\end{split}\]

\[\begin{split}\left[A_{t+k}\right] = \begin{cases} Z_{t+k} - \hat{Z}_{t+k-1}(1) & \text{for } k = 0, -1, -2, \ldots \\ 0 & \text{for } k = 1, 2, \ldots \\ \end{cases}\end{split}\]

The 100(1 - α) percent probability limits for \(Z_{t+l}\) are given by

\[\hat{Z}_t(l) \pm z_{\alpha/2} \left\{1 + \sum_{j=1}^{l-1} \mathit{\Psi}_j^2\right\}^{1/2} \sigma_A\]

where \(z_{(a/2)}\) is the \(100(1-\alpha/2)\) percentile of the standard normal distribution

\[\sigma_A^2\]

(returned from arma) and

\[\left\{\mathit{\Psi}_j\right\}\]

are the parameters of the random shock form of the difference equation. Note that the forecasts are computed for lead times \(l=1,2,\ldots,L\) at origins \(t=(n-b),(n-b+1),\ldots,n\), where L = nPredict and b = backwardOrigin.

The Box-Jenkins forecasts minimize the mean-square error

\[E \left[ Z_{t+l} - \hat{Z}_t(l) \right]^2\]

Also, the forecasts can be easily updated according to the following equation:

\[\hat{Z}_{t+1}(l) = \hat{Z}_t(l+1) + \mathit{\Psi}_l A_{t+1}\]

This approach and others are discussed in Chapter 5, Categorical and Discrete Data Analysis of Box and Jenkins (1976).

Example¶

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. Function armaForecast computes forecasts and 95-percent probability limits for the forecasts for an ARMA(2, 1) model fit using function arma with the method of moments option. With backwardOrigin = 3, columns zero through three of forecasts provide forecasts starting with 1867, 1868, 1869, and 1870, respectively. Note that the values in the first row are the one-step ahead forecasts for 1867, 1868, 1869, and 1870; the values in the second row are the two-step ahead forecasts for 1868, 1869, 1870, and 1871; etc. Column four gives the deviations for computing probability limits, and column five gives the psi weights, which can be used to update forecasts when more data is available. For example, the forecast for the 102nd observation (year 1871) given the data through the 100th observation (year 1869) is 77.21; and 95-percent probability limits are given by \(77.21\mp 56.30\). After observation 101 (\(Z_{101}\) for year 1870) is available, the forecast can be updated by using

\[\hat{Z}_t(l) \pm z_{\alpha/2} \left\{1 + \sum_{j=1}^{l-1} \mathit{\Psi}_j^2\right\}^{1/2} \sigma_A\]

with the psi weight (\(\psi_1=1.37\)) and the one-step-ahead forecast error for observation 101 (\(Z_{101}-83.72\)) to give the following:

\[77.21 + 1.37 × (Z_{101} - 83.72)\]

Since this updated forecast is one step ahead, the 95-percent probability limits are now given by the forecast \(\mp 33.22\).

from __future__ import print_function
from numpy import *
from pyimsl.stat.arma import arma
from pyimsl.stat.armaForecast import armaForecast
from pyimsl.stat.dataSets import dataSets
from pyimsl.stat.writeMatrix import writeMatrix

p = 2
q = 1
n_observations = 100
max_iterations = 0
n_predict = 12
backward_origin = 3
z = empty(100)
rel_error = 0.0
arma_info = []

col_labels = [
    "Lead Time",
    "Forecast From 1866",
    "Forecast From 1867",
    "Forecast From 1868",
    "Forecast From 1869",
    "Dev. for Prob. Limits",
    "Psi"]

w = dataSets(2)
for i in range(0, n_observations):
    z[i] = w[21 + i][1]

parameters = arma(z, p, q,
                  relativeError=rel_error,
                  maxIterations=max_iterations,
                  armaInfo=arma_info)

print("Method of Moments initial estimates:")
print("AR estimates are %11.4f and %11.4f" % (parameters[1], parameters[2]))
print("MA estimate is %11.4f." % parameters[3])

forecasts = armaForecast(arma_info[0], n_predict,
                         backwardOrigin=backward_origin)

writeMatrix("* * * Forecast Table * * *\n",
            forecasts,
            colLabels=col_labels,
            writeFormat="%11.4f")

Output¶

Method of Moments initial estimates:
AR estimates are      1.2443 and     -0.5751
MA estimate is     -0.1241.
 
                     * * * Forecast Table * * *

Lead Time  Forecast From  Forecast From  Forecast From  Forecast From
                    1866           1867           1868           1869
        1        18.2833        16.6150        55.1894        83.7198
        2        28.9182        32.0188        62.7608        77.2094
        3        41.0100        45.8275        61.8924        63.4609
        4        49.9387        54.1496        56.4572        50.0988
        5        54.0937        56.5623        50.1939        41.3803
        6        54.1283        54.7780        45.5268        38.2174
        7        51.7815        51.1701        43.3220        39.2964
        8        48.8417        47.7073        43.2630        42.4581
        9        46.5335        45.4736        44.4577        45.7715
       10        45.3524        44.6861        45.9781        48.0758
       11        45.2103        44.9908        47.1827        49.0372
       12        45.7128        45.8230        47.8072        48.9081
 
Lead Time  Dev. for Prob.          Psi
                   Limits             
        1         33.2179       1.3683
        2         56.2980       1.1274
        3         67.6168       0.6158
        4         70.6432       0.1178
        5         70.7515      -0.2076
        6         71.0869      -0.3261
        7         71.9074      -0.2863
        8         72.5336      -0.1687
        9         72.7498      -0.0452
       10         72.7653       0.0407
       11         72.7779       0.0767
       12         72.8225       0.0720