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 + \[ScriptM] + i - j]], {i, maxLead}, 
          {j, n1 + \[ScriptM] - 1}]; zk = Append[zk, 
          zc[[n1 + \[ScriptM] - 1]]]; GI = TrenchInverse[
          toe[Take[r, n1 + \[ScriptM] - 1]]]; {F, V} = FastMatrixForecast[GK, 
          GI, zk]; zF[[\[ScriptM]]] = F; vF[[\[ScriptM]]] = V; , 
       {\[ScriptM], 2, m + 1}]; {zF + zMean, r[[1]] - vF}]
 
Options[TrenchForecast] = {MaxLeadTime -> 1}
 
TrenchForecast /: TrenchForecast::errOrig = "Origin must be between 1 and ``"
 
TrenchForecast /: TrenchForecast::lengerror = "Length error, Length[r] < ``"
 
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.\!\(\*Cell[\"\"]\)"
 
toe[TrenchInverse`Private`m_] := Table[TrenchInverse`Private`m[[
      1 + Abs[TrenchInverse`Private`i - TrenchInverse`Private`j]]], 
     {TrenchInverse`Private`i, Length[TrenchInverse`Private`m]}, 
     {TrenchInverse`Private`j, Length[TrenchInverse`Private`m]}]
 
toe /: toe::usage = "toe[x] forms the Toeplitz matrix from vector x."
 
TrenchInverse[(TrenchInverse`Private`x_)?MatrixQ] := 
    Module[{TrenchInverse`Private`ncol, TrenchInverse`Private`invx, 
      TrenchInverse`Private`z, TrenchInverse`Private`scz, 
      TrenchInverse`Private`fault, TrenchInverse`Private`EPS = 
       $MachineEpsilon}, TrenchInverse`Private`ncol = 
       Last[Dimensions[TrenchInverse`Private`x]]; TrenchInverse`Private`scz = 
       N[TrenchInverse`Private`x[[1,1]]]; If[TrenchInverse`Private`scz <= 0, 
       Message[TrenchInverse::SingularMatrix]; Abort[]]; 
      TrenchInverse`Private`z = First[TrenchInverse`Private`x]/
        TrenchInverse`Private`scz; 
      TrenchInverse`Private`$TrenchInverseState = {TrenchInverse`Private`scz, 
        TrenchInverse`Private`ncol}; TrenchInverse`Private`invx = 
       TrenchM[TrenchInverse`Private`z, TrenchInverse`Private`EPS]; 
      If[TrenchInverse`Private`invx === $Failed, 
       Message[TrenchInverse::TrenchM]; Abort[]]; 
      TrenchInverse`Private`fault = Round[Last[TrenchInverse`Private`invx]]; 
      If[TrenchInverse`Private`fault != 0, 
       Message[TrenchInverse::SingularMatrix]; Abort[]]; 
      TrenchInverse`Private`invx = Chop[Drop[TrenchInverse`Private`invx, -1]/
         TrenchInverse`Private`scz]; 
      If[VectorQ[TrenchInverse`Private`invx] && 
        And @@ NumericQ /@ TrenchInverse`Private`invx, 
       Partition[TrenchInverse`Private`invx, TrenchInverse`Private`ncol], 
       TrenchInverse`Private`invx]]
 
TrenchInverse /: TrenchInverse::SingularMatrix = 
     "Matrix is singular or not p.d."
 
TrenchInverse /: TrenchInverse::TrenchM = "TrenchM.exe returned $Failed"
 
TrenchInverse /: TrenchInverse::usage = "TrenchInverse[x] uses the Trench \
algorithm to compute the inverse of a symmetric p.d. Toeplitz matrix"
 
TrenchM[r_List, EPS_Real] := ExternalCall[LinkObject["trenchM.exe", 2, 2], 
     CallPacket[0, {r, EPS}]]
 
TrenchM /: TrenchM::usage = "TrenchM[r, eps, mode] MathLink function: invert \
a Toeplitz matrix whose first row is r. Assumed scaled so that r[[1]]=1. \
Returns vector of length square of the length of r plus 1.  The last element \
is a fault indicator. The rest is the matrix stored by row, eps is machine \
epsilon.  For usual inverse  set mode=1.  For recursive updating r is the \
additional elements and mode=2."
 
Attributes[ExternalCall] = {ReadProtected}
 
LinkObject["trenchM.exe"] = LinkObject["trenchM.exe", 2, 2]
 
LinkObject["trenchQFM.exe"] = LinkObject["trenchQFM.exe", 3, 3]
 
FastMatrixForecast = CompiledFunction[{{_Real, 2}, {_Real, 2}, {_Real, 1}}, 
     {{3, 2, 0}, {3, 2, 1}, {3, 1, 2}, {3, 2, 6}}, {0, 5, 2, 0, 8}, 
     {{1, 5}, {86, 0, 1, 3}, {86, 3, 2, 4}, {4, 2, 0}, {56, 0, 0, 2, 5}, 
      {4, 1, 0}, {64, 5, 0, 0, 0, 2}, {4, 0, 3}, {58, 2, 3, 5}, {4, 0, 4}, 
      {9, 4, 0}, {79, 0, 2, 6}, {64, 3, 0, 0, 1, 6}, {64, 0, 0, 0, 1, 7}, 
      {86, 6, 7, 1}, {59, 3, 1, 5}, {42, -5}, {62, 4, 5, 6}, {2}}, 
     Function[{GL, M, z}, With[{P = GL . M}, {P . z, Table[P[[i]] . GL[[i]], 
         {i, First[Dimensions[GL]]}]}]], Evaluate]
