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 + 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`

ARMA Simulation Example

m = 100; n = 200; ℕ = n + m; Muz = 0. ; 
    phi1 = 0.8; phi = {phi1}; theta = {};
maxLead = 7;
SeedRandom[72177791];
r = TacvfARMA[phi, theta, ℕ];
zTOT = 
        SimulateGaussianARMA[phi, 
          theta, ℕ];
z = Take[zTOT, n];
zFF = vFF = 
          Table[0, {maxLead}, {m + 1}];
ans = Array[0 &,  maxLead];

\(Do[
znew = Drop[zTOT, n]; 
zF = \(vF = Array[0 &, m + 1]\); 

      G = toe[Take[r, n]]; 
gk = Take[Drop[r, k], n]; 

      GI = TrenchInverse[G]; 

      zF[[1]] = 
        Muz + Reverse[gk] . GI . \((z - Muz)\); 

      vF[[1]] = 
        r[[1]] - gk . GI . gk; 

      zk = z; 
VarF = Total[Table[phi1\^\(2  i\), {i, 0, k - 1}]]; 

      Do[
gk = 
          Append[gk, 
            r[[n + j - 
                1]]]; 
zk = 
          Append[zk, 
            zTOT[[n + j - 
                1]]]; 
GI = 
          TrenchInverseUpdate[{r[[
                n + j]]}]; 
zF[[j]] = 
          Muz + Reverse[gk] . 
              GI . \((zk - Muz)\); 
vF[[
            j]] = 
          r[[1]] - 
            gk . GI . gk;, 
{j, 2, m + 1}]; 

      ans[[k]] = {Mean[
            errMMSE = \((Drop[znew, k - 1] - Drop[zF, \(-k\)])\)\^2], 
          MeanCI[errMMSE], VarF}; 

      zFF[[k]] = zF; 

      vFF[[k]] = 
        vF; 
Print[" ok done. lead time: ", 
        k];
, {k, maxLead}]\)

\(TableForm[ans, TableDepth -> 2, 
      TableHeadings -> {None, {"Emp MSE", "C.I.", "MSE"}}]\)

Get the matrix, GK

First direct coding based on existing program,

maxLead = 3; n = 10; 
    GK = Array[0 &, maxLead];
r = 
        Table[Subscript[γ, i], {i, 20}];
Do[
gk = 
        Take[Drop[r, k], n]; 
GK[[
          k]] = gk;
, {k, 
        maxLead}]

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,

MatrixForm[GK]

Now for the short direct method,

Table[
        r[[n + i - j + 1]], {i, 
          maxLead}, {j, n}] // MatrixForm

Vectorized Computation

m = 100; n = 200; ℕ = n + m; Muz = 0. ; 
    phi1 = 0.8; phi = {phi1}; theta = {};
maxLead = 7;
SeedRandom[72177791];
r = 
        TacvfARMA[phi, theta, ℕ + maxLead];
zTOT = 
        SimulateGaussianARMA[phi, 
          theta, ℕ];
z = Take[zTOT, n];
ans = Array[0 &,  maxLead];

znew = Drop[zTOT, n];
zF = vF = Array[0 &, m + 1];
G = toe[Take[r, n]];
GK = 
        Table[r[[n + i - j + 1]], {i, 
            maxLead}, {j, n}];
GI = TrenchInverse[G];
zF[[1]] = 
        Muz + GK . GI . (z - Muz);
zk = z;
VarF = 
        Rest[FoldList[Plus, 0, 
            Table[(phi1^(2  i)), {i, 0, 
                maxLead - 1}]]];
Do[
GK = 
          Table[r[[n + 𝓂 + i - j + 
                1]], {i, maxLead}, {j, 
              n + 𝓂 - 1}]; 
zk = 
          Append[zk, 
            zTOT[[n + 𝓂 - 
                1]]]; 
GI = 
          TrenchInverseUpdate[{r[[
                n + 𝓂]]}]; 
zF[[𝓂]] = 
          Muz + GK . 
              GI . (zk - 
                  Muz);
, 
{𝓂, 
          2, m + 1}];

zF = Transpose[zF];
ans = 
        Table[{Mean[
              errMMSE = (((Drop[znew, k - 1] - Drop[zF[[k]], -k]))^2)], MeanCI[errMMSE], 
            VarF[[k]]}, {k, 
            maxLead}];

\(TableForm[ans, TableDepth -> 2, 
      TableHeadings -> {None, {"Emp MSE", "C.I.", "MSE"}}]\)

Special Compile Functions

Symbolic Example

We need to evaluate \(TraditionalForm\`var \((\(\(z\&^\)\_n\)(k))\)\  = γ\_0\  - \ 
          g\_k\%\(\(\ \)\(T\)\)\ \(Γ\_n\%\(-1\)\) g\_k\) for k = 1,  ⋯, L.  Letting 𝔳=\(\((var(\(\(z\&^\)\_n\)(1)), ⋯, 
                  var(\(\(z\&^\)\_n\)(L)))\)\_\(1×L\)\)we can write 𝔳 = 
        diag((Subscript[G, L]) (Subscript[Γ, n]^(-1)) (Subscript[G, L]^′)) or for computation, 𝔳=

GL = 
        Partition[Table[Subscript[a, i], {i, 9}], 
          3];
M = Partition[Table[Subscript[b, i], {i, 9}], 3];

Table[
      GL[[i]] . M . 
          GL[[i]] // Apart, {i, 
        3}]

Transpose[GL . M . Transpose[GL], {1, 1}] // Apart

P = GL . M;

Table[
      P[[i]] . 
          GL[[i]] // Apart, {i, 
        3}]

Numerical Example

GL = 
        Partition[Table[Prime[i], {i, 12}], 4];
M = Partition[Table[(2^i), {i, 16}], 4];

Dimensions[GL]

Dimensions[M]

Table[
      GL[[i]] . M . 
        GL[[i]], {i, 3}]

Transpose[GL . M . Transpose[GL], {1, 1}]

P = GL . M;

Table[
      P[[i]] . 
          GL[[i]] // Apart, {i, 
        3}]

Compile Function

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}]}]
];

GLR = GL // N; MR = M // N; z = Range[4] // N;

FastMatrixForecast[GLR, MR, z]

Faster Vectorized Computation

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}]}]
];

m = 100; n = 200; ℕ = n + m; Muz = 0. ; 
    phi1 = 0.8; phi = {phi1}; theta = {};
maxLead = 7;
SeedRandom[72177791];
r = 
        TacvfARMA[phi, theta, ℕ + maxLead];
zTOT = 
        SimulateGaussianARMA[phi, 
          theta, ℕ];
z = Take[zTOT, n];
zF = vF = Array[0 &, {maxLead + 1}];

znew = Drop[zTOT, n];
zF = vF = Array[0 &, m + 1];
G = toe[Take[r, n]];
GK = 
        Table[r[[n + i - j + 1]], {i, 
            maxLead}, {j, n}];
GI = TrenchInverse[G];
ans = 
        FastMatrixForecast[GK, GI, z - Muz];
F = First[ans];
V = Last[ans];
zF[[1]] = 
        F;
vF[[1]] = 
        V;
zk = z;
VarF = 
        Rest[FoldList[Plus, 0, 
            Table[(phi1^(2  i)), {i, 0, 
                maxLead - 1}]]];
Do[
GK = 
          Table[r[[n + 𝓂 + i - 
                j]], {i, maxLead}, {j, 
              n + 𝓂 - 1}]; 
zk = 
          Append[zk, 
            zTOT[[n + 𝓂 - 
                1]]]; 
GI = 
          TrenchInverseUpdate[{r[[
                n + 𝓂 - 
                  1]]}]; 
{F, V} = 
          FastMatrixForecast[GK, GI, 
            zk - Muz]; 
zF[[𝓂]] = 
          F; 
vF[[𝓂]] = V;
, {𝓂, 2, 
          m + 1}];

zF = Transpose[zF];
ans = 
        Table[{Mean[
              errMMSE = (((Drop[znew, k - 1] - Drop[zF[[k]], -k]))^2)], MeanCI[errMMSE], 
            VarF[[k]]}, {k, 
            maxLead}];

\(TableForm[ans, TableDepth -> 2, 
      TableHeadings -> {None, {"Emp MSE", "C.I.", "MSE"}}]\)

Dimensions /@ {zFF, zF}

eFR = Abs[zFF - zF]/Abs[zFF];

eFR // Chop // Flatten // Max