Application of the Trench Algorithm to Forecasting

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)\)\) and \(TraditionalForm\`var \((\(\(z\&^\)\_n\)(k))\)\  = γ\_0\  - \ 
          g\_k\%\(\(\ \)\(T\)\)\ \(Γ\_n\%\(-1\)\) g\_k\).  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 + j\)\)(k) = μ + 
          g\_\(n, j\)\ \(\(Γ\_\(n + j\)\%\(-1\)\)(
              z\^\((n + j)\) - μ)\)\)\(,\)\)\)

where \(TraditionalForm\`\(\(\ \)\(g\_\(n, j\) = \ \((γ\_\(n + j + k - 2\)\ , \  ... \ , \ γ\_k)\)\)\)\) and (z^((n + j))) = ((Subscript[z, 1]),   ...  ,  
          (Subscript[z, n + j - 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.

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

Further Development Needed

Only for lead time 1.  Need to compute for l = 1,  ... , L for general L.  Result should be a pair of matricies.

Initialize

<< FGN.m

<< FitARMA.m

<< Statistics`

New Functions

TrenchForecast::usage="TrenchForecast[z,r,mean,n,L] for series of length n+m in z, the forecasts & their variances starting at origin n with lead time L are produced. The vector r are the theoretical autocovariances at lags 0,...,n+m and the mean parameter, mean. The dimension of the result is (m+1)-by-2, corresponding to (Subscript[z, n + j])(L),j=0...m";

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

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

FastIn = 
        Compile[{{a, _Real, 1}, {b,  _Real, 1}}, Total[a b]];

FastbA = 
        Compile[{{b,  _Real, 1}, {A, _Real, 2}}, b . A];

\(TrenchForecast[z_, \ r_, zMean_, ForecastOrigin_, \ LeadTime_] := 
      Module[{n = Length[z], zc, n1, m, z1, zF, G, gk, GI, gkGI, vF, 
          zk}, 
If[\((n + 1)\) ≠ Length[r], 
          Message[TrenchForecast::"lengerror"]; 
          Abort[]]; 
If[\((ForecastOrigin > 
                n)\) || \((ForecastOrigin < 1)\), 
          Message[TrenchForecast::"errOrig", n]; 
          Abort[]]; 
zc = 
          z - zMean; 
n1 = 
          ForecastOrigin; 
m = 
          n - n1; 
z1 = 
          Take[zc, n1]; 
zF = \(vF = 
            Array[0.0 &, m + 1]\); 
G = 
          toe[Take[r, n1]]; 
gk = 
          Reverse[Take[Drop[r, 1], n1]]; 
GI = 
          TrenchInverse[G]; 
gkGI = 
          FastbA[gk, GI] // 
            N; 
zF[[1]] = zMean + 
            FastIn[gkGI, 
              z1]; 
vF[[1]] = 
          r[[1]] - 
            FastIn[gkGI, gk]; 
zk = 
          z1; 
Do[
gk = 
            Prepend[gk, 
              r[[n1 + j - 
                  1]]]; 
zk = 
            Append[zk, 
              zc[[n1 + j - 
                  1]]]; 
GI = 
            TrenchInverseUpdate[{r[[
                  n1 + j]]}]; 
gkGI = 
            FastbA[gk, GI]; 
zF[[
              j]] = 
            zMean + FastIn[gkGI, 
                zk]; 
vF[[
              j]] = 
            r[[1]] - 
              FastIn[gkGI, gk];
, 
{j, 2, 
            m + 1}
]; 
{zF, 
          vF}]\)

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, 1];

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

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

AR Simulation Example - n=n1=100

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

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

AR Simulation Example - Fast

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

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

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

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

Short[vF, 3]

AR Simulation Example

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

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

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

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

Short[vF, 3]

ARMA Simulation Example

n = 400; Muz = 0. ; n1 = 200;
r = TacvfARMA[phi, theta, n];
phi = {1.2399401603634308, -0.2653412932838738};
theta = {0.8436410920940121};
z = SimulateGaussianARMA[phi, theta, n];

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

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

{Mean[errMMSE = (((znew - zF))^2)], 
      Mean[errM = (((znew - Mean[z1]))^2)], 
      MeanDifferenceCI[errMMSE, errM]}

Short[vF, 3]

FGN Simulation Example

m = 200; n1 = 200; n = n1 + m; H = 0.8; 
    Muz = 0. ;
r = Table[FGNAcf[k, H], {k, 0, n - 1}];
z = SimulateFGN[H, n];
z1 = Take[z, n1];

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

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

{Mean[errMMSE = (((znew - zF))^2)], 
      Mean[errM = (((znew - Mean[z1]))^2)], 
      MeanDifferenceCI[errMMSE, errM]}

Short[vF, 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

ans = FitFGN[z1, MeanMLEQ -> True]

Hz = HurstCoefficient /. Last[ans];
Muz = MuHat /. Last[ans];
r = Table[FGNAcf[i - 1, Hz], {i, n + 1}];

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

{zF, vF} = ans;

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

ARMA Forecast

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

phi = Take[ARMACoefficients /. ans, 2];
theta = Drop[ARMACoefficients /. ans, 2]; 
Muz = MeanValue /. ans;
r = TacvfARMA[phi, theta, n];

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

{zF, vF} = ans;

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