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

?TrenchInverseUpdate

<< FitARMA.m

<< Statistics`

New Functions

TrenchForecast::usage="TrenchForecast[z,r,mu,n] 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 and the mu is the mean parameter. The result is a list with (m+1) components. Each component is a matrix with 2 columns correponding to the forecasts and their variances.  Multi-step forecasts can be obtained using the option, MaxLeadTime.  TrenchForecast[z,r,mu,n,MaxLeadTime->L] computes forecasts for all lead times 1, ..., L for origins n, n+1,...,n+m.";

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

TrenchForecast::"lengerror" = "Length error";

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_, \ opts___?OptionQ] := 
      Module[{n = Length[z], zc, n1, m, z1, G, gk, GI, gkGI, vF, zF, zk, F, 
          V, maxLead}, 
If[\((n + 1)\) ≠ 
            Length[r], Message[TrenchForecast::"lengerror"]; 
          Abort[]]; 
If[\((ForecastOrigin > 
                n)\) || \((ForecastOrigin < 1)\), \ \ \ Message[
            TrenchForecast::"errOrig", n]; Abort[]]; 
        maxLead = \(MaxLeadTime /. {opts}\) /. 
            Options[TrenchForecast]; 
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 = 
            TrenchInverseUpdate[{r[[
                  n1 + 𝓂 - 
                    1]]}]; 
{F, V} = 
            FastMatrixForecast[GK, GI, 
              zk]; 
zF[[𝓂]] = 
            F; 
vF[[𝓂]] = V;
, {𝓂, 2, 
            m + 1}]; 
{zF + zMean, 
          r[[1]] - vF}]\)

Fast22 = 
        Compile[{{GI, _Real, 2}, {v, _Real}, {g, _Real, 
              1}}, 
1/(v - g . GI . g)];

Fast12 = 
        Compile[{{GI, _Real, 2}, {e, _Real}, {g, _Real, 1}}, -GI . g*
            e];

Fast11 = 
        Compile[{{GI, _Real, 2}, {e, _Real}, {g, _Real,  
              1}}, 
GI + (GI . 
                  Table[g[[i]]*
                      g[[j]], {i, 
                      Length[g]}, {j, Length[g]}] . GI)*e];

TrenchInverseUpdate[GI_, r_, rnew_] := 
      Module[{v = First[r], 
          g = {r, rnew} // Flatten // Reverse // Most, e, f, 
          A}, 
e = Fast22[GI, v, g]; 
f = 
          Fast12[GI, e, g]; 
A = 
          Fast11[GI, e, 
            g]; 
Partition[{Partition[Flatten[{A, f}], 
                  Length[f]] // Transpose, f, e} // Flatten, 
          Length[f] + 1]
]

\(\(TrenchInverseUpdate::usage = 
"TrenchInverseUpdate[GI,r,rnew] uses partitioned matrix result to obtain the order n+1 inverse given the order n inverse of G, GI; the first row of G, r; and the next element of the first row rnew.";\)\)

Save Package

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

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

AR Simulation Test Error Conditions

n = 100; Muz = 0. ; n1 = 90;
r = TacvfARMA[phi, theta, n];
phi = {0.8};
theta = {};
z = SimulateGaussianARMA[phi, theta, n];
z1 = Take[z, n1];
r1 = Rest[r];

ans = TrenchForecast[z, r1, 0.0, n1];

ans = TrenchForecast[z, r1, 0.0, n + 1];

ans = TrenchForecast[z, r1, 0.0, 0];

AR Simulation Example - n=n1=100

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

ans = TrenchForecast[z, r, 0.0, n1]

AR Simulation Example - Fast

n = 100; Muz = 0. ; n1 = 90;
phi = {0.8};
theta = {};
sig = 8.0;
r = TacvfARMA[phi, theta, n]*(sig^2);
z = 
        SimulateGaussianARMA[phi, theta, n]*sig;
z1 = Take[z, n1];

ans = TrenchForecast[z, r, 0.0, n1];

ans

{zF, vF} = ans; znew = Drop[z, n1];

{Mean[errMMSE = (((znew - Most[zF // Flatten]))^2)], 
      MeanCI[errMMSE]}

Short[vF // Flatten, 3]

ARMA Simulation Example

n = 100; Muz = 0. ; n1 = 80;
phi = {0.8};
theta = {0.4};
sig = 9. ;
r = TacvfARMA[phi, theta, n]*(sig^2);
z = SimulateGaussianARMA[phi, theta, n]*sig;

ans = TrenchForecast[z, r, 0.0, n1];

{zF, vF} = ans;
znew = Drop[z, n1];

{Mean[errMMSE = (((znew - Most[zF // Flatten]))^2)], 
      MeanCI[errMMSE]}

Short[vF // Flatten, 3]

ARMA Simulation Example

n = 400; Muz = 0. ; n1 = 200;
phi = {1.2399401603634308, -0.2653412932838738};
theta = {0.8436410920940121};
sig = 100.0;
r = TacvfARMA[phi, theta, n]*(sig^2);
z = SimulateGaussianARMA[phi, theta, n]*sig;

Testing computation from first principles

G = toe[Take[r, n1]];
GI = Inverse[G];
gk = Take[Drop[r, 1], n1];
r[[1]] - gk . GI . gk

G = toe[Take[r, n1]];
GI = TrenchInverse[G];
gk = Take[Drop[r, 1], n1];
r[[1]] - gk . GI . gk

ans = TrenchForecast[z, r, 0.0, n1];

{zF, vF} = ans;
z1 = Take[z, n1];
znew = Drop[z, n1];

{Mean[errMMSE = (((znew - Most[zF // Flatten]))^2)], 
      MeanCI[errMMSE]}

Short[vF // Flatten, 3]

FGN Simulation Example

m = 200; n1 = 200; n = n1 + m; H = 0.8; Muz = 0. ; 
    sig = 9. ;
r = Table[FGNAcf[k, H], {k, 0, n}]*(sig^2);
z = SimulateFGN[H, n];
z1 = Take[z, n1];

Variance[z]

Variance[SimulateFGN[H, (2^15)]]

?SimulateFGN

ans = TrenchForecast[z, r, 0.0, n1];

{zF, vF} = ans;
znew = Drop[z, n1];

{Mean[errMMSE = (((znew - Most[zF // Flatten]))^2)], 
      MeanCI[errMMSE]}

Short[vF // Flatten, 3]

Application to Nile Minimums

Percival & Walden comment in Sec 5.7 (p. 190), "Toussoun's record extends up to 1921.  It is complete for 622-1284.  For 1285-1470 there are scattered missing values."  We linearly interpolated to complete the record.

Initialize

z = ReadData["d:/math/2006/Trench/data/Minimum.1"];

n1 = 1284 - 622 + 1;
n = Length[z];
m = n - n1;
z1 = Take[z, n1]; znew = Drop[z, n1];

FGN Case

ansFit = FitFGN[z1, MeanMLEQ -> True]

Hz = HurstCoefficient /. Last[ansFit];
Muz = MuHat /. Last[ansFit];
r = 
        Table[FGNAcf[i - 1, Hz], {i, n + 1}]*
          Variance[Take[z, n1]];

ans = TrenchForecast[z, r, Muz, n1];

{zF, vF} = ans;

{Mean[errMMSE = (((znew - Most[zF]))^2)], MeanCI[errMMSE]}

Short[vF // Flatten, 5]

Short[zF // Flatten, 5]

zFFGN = zF;

Save Data

DeleteFile["d:/math/2006/Trench/data/FGNNile.m"]

Save["d:/math/2006/Trench/data/FGNNile.m", {zFFGN}]

ARMA Forecast

ans = FitARMA[z, {2, 0, 1}];
CoefARMA[ans]

\(\(\(\(phi = Take[ARMACoefficients /. ans, 2];\)

    \(theta = Drop[ARMACoefficients /. ans, 2];\)

    \(Muz = MeanValue /. ans;\)

    \(Sig = \ ResidualSd /. ans;\)

    \(\(r = TacvfARMA[phi, theta, n]*Sig\^2\) \); \)\(\ \)\( (*VNB : \ 
        must\ be\ \(autocovariances!\)*) \)\)\)

TrenchForecast[z, r, Muz, n]

ans = TrenchForecast[z, r, Muz, n1];

{zF, vF} = ans;

{Mean[errMMSE = (((znew - Most[zF]))^2)], MeanCI[errMMSE]}

Naturally the out-of-sample mse is larger than the residual mse.

Short[z]

Short[zF]

Short[vF]

Computation from first principles.  This agrees with the residual sd.

G = toe[Take[r, n1]];
GI = Inverse[G];
gk = Take[Drop[r, 1], n1];
r[[1]] - gk . GI . gk // 
      Sqrt