Application of the Trench Algorithm to Forecasting

Introduction

Let \(TraditionalForm\`\(z\&^\_n\)(k)\) denote the minimum mean square error linear predictor of (Subscript[z, n + k]) given the data z = ((Subscript[z, 1]),   ...  ,  (Subscript[z, n])) and the autocovariances (Subscript[γ, l]),  l = 0,  ... , n - 1.  Then as shown by Hamilton (1994, §4.3) the finite-sample predictor may be written

\(TraditionalForm\`\(z\&^\_n\)(k) = μ + 
          g\_k\ \(\(Γ\_n\%\(-1\)\)(z - μ)\)\)

where \(TraditionalForm\`g\_k\  = \ \((γ\_\(n + k - 1\)\ , \  ... \ , \ γ\_k)\)\).  Note that (1), the autocovariances may be replaced by autocorrelations.  The variance of the forecast is

\(TraditionalForm\`var(\(z\&^\_n\)(k))\  = γ\_0\  - \ 
          g\_k\%\(\(\ \)\(T\)\)\ \(Γ\_n\%\(-1\)\) g\_k\).

Note that in this case the autocovariances are required.

In many forecasting applications, a model is fitted to the data to the first n observed time series values and then k-step ahead forecasting are needed for each of m additional observations, ((Subscript[z, n + 1]),   ...  ,  (Subscript[z, n + m])).  These forecasts may be written

\(TraditionalForm\`\(\(\(z\&^\_\(n + l\)\)(k) = μ + 
          g\_\(n, l\)\ \(\(Γ\_\(n + j\)\%\(-1\)\)(
              z\^\((n + l)\) - μ)\)\)\(,\)\)\)

where \(TraditionalForm\`\(\(\ \)\(g\_\(n, l\) = \ \((γ\_\(n + l + k - 2\)\ , \  ... \ , \ γ\_k)\)\)\)\) and (z^((n + l))) = ((Subscript[z, 1]),   ...  ,  
          (Subscript[z, n + l - 1])).  Rather than evaluating m matrix inverses of dimension n + 1,  ... ,  n + m we just need to evaluate (Subscript[Γ, n]^(-1)) and the remaining inverses can be obtained by updating our comptations.  This effectively reduces the matrix inverse computations from O(((m(n + m))^2)) flops to just O(n (n + m)) flops.  Hamilton (1994, §4.4) discusses the use of the Cholesky decomposition for finite sample forecast updating but the Trench algorithm provides a more computationally efficient approach.

More generally, we may be interested in computing the k-step predictor for k = 1, ⋯, L.  The vector of forecasts \(TraditionalForm\`\(z\&^\)\_n\%\((L)\) = \((\(z\_n\)(
            1), ⋯, \ \(z\_n\)(L))\)\) may be written,

\(TraditionalForm\`\(z\&^\)\_n\%\((L)\) = μ + 
          G\_L\ \(\(Γ\_n\%\(-1\)\)(z - μ)\)\)

where \(TraditionalForm\`\(\(G\_L = \((γ\_\(n + i - j\))\)\_\(n×L\)\)\(,\)\)\) where i = 1, ⋯, L and j = 1, ⋯, n.  For updating from n to n + m we can write,

\(TraditionalForm\`\(\(\(z\&^\)\_\(n + m\)\%\((L)\) = μ + 
          G\_\(L, m\)\ \(\(Γ\_\(n + m\)\%\(-1\)\)(
              z\^\((n + m)\) - μ)\)\)\(,\)\)\)

where \(TraditionalForm\`\(\(G\_\(L, m\) = \((γ\_\(n + m + i - j - 1\))\)\_\(\((n + m - 1)\)×L\)\)\(,\)\)\) where i = 1, ⋯, L and j = 1, ⋯, n + m - 1.

Hamilton, J.D. (1994).  Time Series Analysis.  New Jersey: Princeton University Press.

Initialize

<< FGN.m

<< FitARMA.m

<< Statistics`

New Functions

Clear[TrenchForecast];

TrenchForecast::usage="TrenchForecast[z,r,mu,n,L] for series of length n+m in z, the one-step ahead forecasts & their variances for origin times n, n+1, ..., n+m are produced. The vector r is the theoretical autocovariances at lags 0,...,n+m+L-1 and the mu is the mean parameter. The result is a list with 2 components corresponding to the forecasts and their variances. Each component is a matrix with L columns correponding to the lead-time and m+1 rows corresponding to the forecast origin.";

Options[TrenchForecast] = {MaxLeadTime -> 1};

\(\(TrenchForecast::"lengerror" = "Length error, Length[r] < ``";\)\)

TrenchForecast::"errOrig" = "Origin must be between 1 and \
``";

FastMatrixForecast = 
        Compile[{{GL,  _Real, 2}, {M, _Real, 2}, {z, _Real, 1}}, 
          With[{P = GL . M}, {P . z, 
              Table[P[[i]] . 
                  GL[[i]], {i, 
                  Dimensions[GL] // 
                    First}]}]
];

\(TrenchForecast[z_, \ r_, zMean_, ForecastOrigin_, \ MaxLeadTime_] := 
      Module[{n = Length[z], zc, n1, m, z1, G, gk, GI, gkGI, vF, zF, zk, F, 
          V, maxLead}, 
If[\((ForecastOrigin > 
                n)\) || \((ForecastOrigin < 1)\), \ \ \ Message[
            TrenchForecast::"errOrig", n]; 
          Abort[]]; 
maxLead = 
          MaxLeadTime; 
If[\((n + maxLead)\) > \ Length[r], 
          Message[TrenchForecast::"lengerror", n + maxLead]; 
          Abort[]]; 
zc = 
          z - zMean; 
n1 = 
          ForecastOrigin; 
m = 
          n - n1; 
z1 = 
          Take[zc, n1]; 
zF = \(vF = 
            Array[0 &, m + 1]\); 
G = 
          toe[Take[r, n1]]; 
GK = 
          Table[r[[n1 + i - j + 1]], {i, 
              maxLead}, {j, n1}]; 
GI = 
          TrenchInverse[G]; 
ans = 
          FastMatrixForecast[GK, GI, z1]; 
F = 
          First[ans]; 
V = 
          Last[ans]; 
zF[[1]] = 
          F; 
vF[[1]] = V; 
zk = 
          z1; 
Do[
GK = 
            Table[r[[n1 + 𝓂 + i - 
                  j]], {i, maxLead}, {j, 
                n1 + 𝓂 - 1}]; 
zk = 
            Append[zk, 
              zc[[n1 + 𝓂 - 
                  1]]]; 
GI = 
            TrenchInverse[
              toe[Take[r, n1 + 𝓂 - 1]]]; 
{F, V} = 
            FastMatrixForecast[GK, GI, 
              zk]; 
zF[[𝓂]] = 
            F; 
vF[[𝓂]] = V;\ 
, {𝓂, 2, 
            m + 1}]; 
{zF + zMean, 
          r[[1]] - vF}]\)

Save Package

DeleteFile["d:/math/2006/Trench/code/forecast.m"]

Save["d:/math/2006/Trench/code/forecast.m", {TrenchForecast}]

AR Simulation Example -

n = 10; Muz = 0. ; n1 = 5; maxl = 3;
phi = {0.8};
theta = {};
sig = 5. ;
r = 
        TacvfARMA[phi, theta, n + maxl - 1]*
          (sig^2);
z = 
        SimulateGaussianARMA[phi, theta, n]*sig;
z1 = Take[z, n1];

{F, V} = ans = TrenchForecast[z, r, Muz, n1, maxl]

MatrixForm[F]

MatrixForm[V]

AR(1) Direct Forecast Computations

(Subscript[z, n])(
          l) = μ + (ϕ^l)((Subscript[z, n]) - μ) and (Subscript[V, l]) = (Subscript[σ, a]^2)(
            1 - (ϕ^l))/(1 - (ϕ^2))

\(AR1[z_, \ phi_, sigsq_, zMean_, ForecastOrigin_, \ MaxLeadTime_] := 
      Module[{}, 
maxLead = 
          MaxLeadTime; 
If[\((n + maxLead)\) ≠ 
            Length[r], Message[TrenchForecast::"lengerror", n + maxLead]; 
          Abort[]]; 
zc = 
          z - zMean; 
n1 = 
          ForecastOrigin; 
m = 
          n - n1; 
z1 = 
          Take[zc, n1]; 
zF = \(vF = 
            Array[0 &, m + 1]\); 
F = 
          Table[zMean + 
              phi\^l*\((z[[n1]] - 
                    zMean)\), {l, maxLead}]; 
V = 
          sigsq*\ Table[\((1 - phi\^\(2  l\))\)/\((1 - phi\^2)\), {l, 
                maxLead}]; 
zF[[1]] = 
          F; 
vF[[1]] = V; 
zk = 
          z1; 
Do[
F = 
            Table[zMean + 
                phi\^l*\((z[[
                        n1 + 𝓂 - 1]] - 
                      zMean)\), {l, 
                maxLead}]; 
zF[[𝓂]] = 
            F; 
vF[[𝓂]] = V;\ 
, {𝓂, 2, 
            m + 1}]; 
{zF, vF}]\)

{F1, V1} = 
      AR1[z, phi[[1]], sig^2, Muz, n1, 
        maxl]

MatrixForm[F1]

MatrixForm[V1]