Trench Forecast - R code development

Initialization

<< TrenchInverse.m

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

z = 
        ReadList["d:/r/2006/trench/ar1sim.dat", 
          Real];
z = Take[z, 10];
zMean = Mean[z];
n = Length[z];
r  =  
        With[{phi = 0.8}, (1/(1 - phi^2))*
            phi^Range[0, n]];
ForecastOrigin  =  5;
maxL =  1;

TrenchForecast Algorithm - Lead 1 Case

?TrenchForecast

\(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 = 
            TrenchInverse[
              toe[Take[r, n1 + 𝓂 - 1]]]; 
{F, V} = 
            FastMatrixForecast[GK, GI, 
              zk]; 
zF[[𝓂]] = 
            F; 
vF[[𝓂]] = V;
, {𝓂, 2, 
            m + 1}]; 
{zF + zMean, 
          r[[1]] - vF}]\)

TrenchForecast[z,  r, zMean, ForecastOrigin, 
      MaxLeadTime ->  maxL]

Deconstructed

Part 1

maxLead = maxL;
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];

GK

GI

GI

gkGI

z

z1

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

ans = FastMatrixForecast[GK, GI, z1]

With[{P = GK . GI}, {P . z1, 
        Table[P[[i]] . 
            GK[[i]], {i, 
            Dimensions[GK] // First}]}]

F = First[ans];
V = Last[ans];
zF[[1]] = 
        F;
vF[[1]] = V;

Part 2

zk = z1;
\(\(Do[
GK = 
          Table[r[[n1 + 𝓂 + i - 
                j]], {i, maxLead}, {j, 
              n1 + 𝓂 - 1}]; 
Print["GK=", 
          GK]; 
zk = 
          Append[zk, 
            zc[[n1 + 𝓂 - 
                1]]]; 
GI = 
          TrenchInverse[
            toe[Take[r, n1 + 𝓂 - 1]]]; 
{F, V} = 
          With[{P = GK . GI}, {P . zk, 
              Table[P[[i]] . 
                  GK[[i]], {i, 
                  Dimensions[GK] // 
                    First}]}]; 
zF[[𝓂]] = 
          F; 
vF[[𝓂]] = V;
, {𝓂, 2, 
          m + 1}];\)\)
{zF + zMean, 
      r[[1]] - vF}

TrenchForecast Algorithm - Lead 5 Case

Result

?TrenchForecast

\(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}]\)

z = 
        ReadList["d:/r/2006/trench/ar1sim.dat", 
          Real];
z = Take[z, 10];
zMean = Mean[z];
n = Length[z];
maxL = 5;
phi = 0.8;
r  =  (1/(1 - phi^2))*
          phi^Range[0, n + maxL - 1];
ForecastOrigin  =  5;

TrenchForecast[z,  r, zMean, ForecastOrigin,  maxL]

Deconstructed

maxLead = maxL;
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;

F

(GK . GI) . z1

V

Table[(GK . GI)[[i]] . 
        GK[[i]], {i, 
        Dimensions[GK] // First}]

z1

GK // MatrixForm

MatrixForm[GK . GI] // Chop

F

zMean + F

\(\(Do[
GK = 
          Table[r[[n1 + 𝓂 + i - 
                j]], {i, maxLead}, {j, 
              n1 + 𝓂 - 1}]; 
Print[" GK = ", 
          MatrixForm[GK]]; 
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}];\)\)
{F, 
        V} = ans = {zF + zMean, 
          r[[1]] - vF}

Dimensions[ans]

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[{F, V}, 
maxLead = 
          MaxLeadTime; 
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}]

sig = 1;
{F1, V1} = ans = AR1[z, phi, sig^2, zMean, n1, maxL]

MatrixForm[F1]

MatrixForm[V1]

Abs[F1 - F] // Flatten // Total

Abs[V1 - V] // Flatten // Total

Second Example

z = ReadList["d:/r/2006/trench/z.dat", Real];

AR1[z,  0.8, 30^2, 0, 5, 1]

r = 
        900*(0.8^Range)[0, Length[z] + 1]/(1 - 
                (0.8^2));
ans = TrenchForecast[z,  r, 0, 5,  1]

First[ans] // Flatten

\(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]; 
Print["GK=", 
          GK]; 
Print["z1=", 
          z1]; 
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}]\)