Print[StyleForm[StringJoin[
        "Updated: August 14, 2006"
        ]
  ,FontSize->9,FontFamily->"Times", FontColor->
  RGBColor[0., 0.25,0.25]
   ]];

(* Author:  A.I. McLeod, aimcleod@uwo*)

(* Summary:
Inverse of p.d. Toeplitz matrix with applications to time series analysis.
For examples see the notebook TrenchInverse-Documentation.nb.
*)

BeginPackage["TrenchInverse`"]

(*==================Install MathLink Executable: TrenchM==============*)

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.";

linkM = Install[LinkLaunch["trenchM.exe"]];

Print[StyleForm[ToString[linkM],FontSize->12,FontFamily->"Courier", 
       FontColor->RGBColor[0,0.25,0.25]]]; 

(*==================Install MathLink Executable: TrenchQFM=============*)

TrenchQFM /: TrenchQFM::usage="TrenchQFM[r,z,EPS] MathLink function: quadratic form and determinant. \
r and z containing acf and data. EPS is machine epsilon.
Assumes that r[1]=1. Returns vector of length 3 containing quadratic form, determinant \
and fault indicator, = 0, ok; =1, singular matrix; =2, r[1] != 1;  =3, 
length(r)!=length(z).";

linkQFM = Install[LinkLaunch["trenchQFM.exe"]];

Print[StyleForm[ToString[linkQFM],FontSize->12,FontFamily->"Courier", 
       FontColor->RGBColor[0,0.25,0.25]]]; 

(*==================Trench Inverse===============================*)
 
toe /: toe::usage = "toe[x] forms the Toeplitz matrix from vector x.";

toeQ /: toeQ::usage = "toeQ[x] test if x is a Toeplitz matrix";

TrenchInverse /: TrenchInverse::usage = "TrenchInverse[x] uses the Trench \
algorithm to compute the inverse of a symmetric p.d. Toeplitz matrix";

(*=======Forecasting=====================*)

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.)";
 
FastInnovationVarianceApproximation::usage = 
    "FastInnovationVarianceApproximation[r] computes the 1-step ahead \
variance for forecast based on n=Length[r]-1 observations. For large n, this \
is approximately equal to the innovation variance.";
   
(*=======Loglikelihood, Mean, Innovation Variance=====================*)

TrenchLoglikelihood /: TrenchLoglikelihood::usage = "TrenchLoglikelihood[r,z] \
loglikelihood function given theoretical autocorrelation function r and data \
z. Assumes z is mean-corrected.  Assumes \
r\[LeftDoubleBracket]1\[RightDoubleBracket]=1."

TrenchMean /: TrenchMean::usage = "TrenchMean[c, z] computes the exact mle \
for the mean given a series z and autocovariance function c. Also referred to \
the BLUE estimate.";

TrenchInnovationVariance /: TrenchInnovationVariance::usage = "TrenchInnovati\
onVariance[c, z] computes the exact mle for the innovation variance given a \
series z and standardized autocovariance function c. The standardized \
autocovariance function is the autocovariance function divided by \
the innovation variance.";

NumericVectorQ /: NumericVectorQ::usage="NumericVectorQ[x] tests if \
x is a numeric vector.";

(*==================Error Messages=======================*)
TrenchInverse /: TrenchInverse::SingularMatrix = 
     "Matrix is singular or not p.d.";

TrenchInverse /: TrenchInverse::TrenchM = "TrenchM.exe returned $Failed";

TrenchLoglikelihood /: TrenchLoglikelihood::TrenchQFM = 
     "TrenchQFM.exe returned $Failed";
 
TrenchLoglikelihood /: TrenchLoglikelihood::UnequalLengths = 
     "Invalid input, length of both arguments must be equal";
  
TrenchInnovationVariance /: TrenchInnovationVariance::UnequalLengths = 
     "Invalid input, length of both arguments must be equal";

TrenchInnovationVariance /: TrenchInnovationVariance::SingularMatrix = 
     "Matrix is singular or not p.d."
 
TrenchInnovationVariance /: TrenchInnovationVariance::TrenchQFM = 
     "TrenchQFM.exe returned $Failed"
 
TrenchMean /: TrenchMean::UnequalLengths = 
     "Invalid input, length of both arguments must be equal";

TrenchMean /: TrenchMean::SingularMatrix = "Matrix is singular or not p.d.";
 
(*==================DLLoglikelihood=======================*)

DLLoglikelihood::usage=
"DLLoglikelihood[c, z] 
computes the exact loglikelihood for a time series with 
data z and autocovariance function c. 
The mean is assumed to be zero.
If c has length less than z, it is padded with zeros. 
The loglikelihood returned is -1/2*Log[g]-n/2 Log[S/n], where n=Length[z] 
g is the covariance determinant and S
is the sum of the standardized prediction errors.";

(*==================Compile Functions====================*)
FastIn::usage="fast inner product";
FastbA::usage="fast b.A";

(********************************************************************) 
Begin["`Private`"]
(********************************************************************) 

toe[m_] := Table[m[[1 + Abs[i - j]]], {i, Length[m]}, {j, Length[m]}];

toeQ[m_] := MatrixQ[m] && $MachineEpsilon > Max[Abs[m - toe[First[m]]]]

TrenchInverse[(x_)?MatrixQ] := Module[{ncol, invx, z, scz, fault, 
      EPS = $MachineEpsilon}, ncol = Last[Dimensions[x]]; scz = N[x[[1,1]]]; 
      If[scz <= 0, Message[TrenchInverse::SingularMatrix]; Abort[]]; 
      z = First[x]/scz; $TrenchInverseState = {scz, ncol}; 
      invx = TrenchM[z, EPS]; If[invx === $Failed, 
       Message[TrenchInverse::TrenchM]; Abort[]]; fault = Round[Last[invx]]; 
      If[fault != 0, Message[TrenchInverse::SingularMatrix]; Abort[]]; 
      invx = Chop[Drop[invx, -1]/scz]; 
      If[VectorQ[invx] && And @@ NumericQ /@ invx, Partition[invx, ncol], 
       invx]];

(****************************Forecasting*****************************) 

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}];
 
FastMatrixForecast=
    Compile[{{GL, _Real,2},{M,_Real,2},{z,_Real,1}},
      With[{P=GL.M},{P.z,
          Table[P\[LeftDoubleBracket]
                i\[RightDoubleBracket].GL\[LeftDoubleBracket]
                i\[RightDoubleBracket],{i,Dimensions[GL]//First}]}]
      ];

FastInnovationVarianceApproximation[r_] := Module[{n1 = Length[r] - 1, GI},
    GI = TrenchInverse[toe[Take[r, n1]]];
    FIVA[r, GI]
    ];

FIVA = Compile[{{r, _Real, 1}, {GI, _Real, 2}},
      r\[LeftDoubleBracket]1\[RightDoubleBracket] - 
        With[{gk = Take[Drop[r, 1], Length[r] - 1]}, Total[gk*(GI.gk)]]
      ];
 
  
(****************************Loglikelihood*****************************) 

TrenchLoglikelihood[(c_)?NumericVectorQ, (y_)?NumericVectorQ] := 
    Module[{r, z, n = Length[y], d, s, ans, EPS = $MachineEpsilon}, 
     If[n != Length[c], Message[TrenchLoglikelihood::UnequalLengths]; 
        Abort[]]; If[c[[1]] <= 0, 
       Message[TrenchLoglikelihood::SingularMatrix]; Abort[]]; r = c/c[[1]]; 
      ans = TrenchQFM[r, y, EPS]; 
      If[invx === $Failed,Message[TrenchLoglikelihood::TrenchQFM]; Abort[]]; 
      {s, ld, fault} = ans; 
      If[fault != 0,Message[TrenchLoglikelihood::SingularMatrix]; Abort[]]; 
      (-(n/2))*Log[s/n] - 0.5*ld];

TrenchMean[(c_)?NumericVectorQ, (z_)?NumericVectorQ] := 
    Module[{G1, r, c0 = First[c], Null}, 
     If[Length[z] != Length[c], Message[TrenchMean::UnequalLengths]; 
        Abort[]]; If[c0 <= 0, Message[TrenchMean::SingularMatrix]; Abort[]]; 
      r = N[c/c0]; G1 = Total[Transpose[TrenchInverse[toe[r]]]]; 
      muhat = Total[G1*z]/Total[G1]];
 
TrenchInnovationVariance[(c_)?NumericVectorQ, (z_)?NumericVectorQ] := 
    Module[{n = Length[z], c0 = First[c], r, EPS = $MachineEpsilon, ans}, 
     If[n != Length[c], Message[TrenchInnovationVariance::UnequalLengths]; 
        Abort[]]; If[c0 <= 0, 
       Message[TrenchInnovationVariance::SingularMatrix]; Abort[]]; r = c/c0; 
      ans = TrenchQFM[N[r], N[z], EPS]; If[ans === $Failed, 
       Message[TrenchInnovationVariance::TrenchQFM]; Abort[]]; 
      If[Last[ans] != 0, Message[TrenchInnovationVariance::SingularMatrix]; 
        Abort[]]; First[ans]/(c0*n)];
 
DLLoglikelihood[c_List, z_List, opts___?OptionQ] := 
      Module[{w, a, phik = {}, vk, phikm1, vkm1,d = Drop[c, 1], n = Length[z]}, 
      w = z;
      vk = c[[1]]; phik = {c[[2]]/c[[1]]}; 
      vk = c[[1]]*(1 - First[phik]^2); 
      S = w[[1]]^2/c[[1]] + (w[[2]] - First[phik]*w[[1]])^2/vk; 
      g = c[[1]]*vk; 
      Do[vkm1 = vk; phikm1 = phik; 
         a = Prepend[-phikm1, 1] . Reverse[Take[d, k]]/vk; 
         phik = Append[phikm1 - a*Reverse[phikm1], a]; 
         vk = vkm1*(1 - a^2); g *= vk; 
         S += (Prepend[-phik, 1] . Reverse[Take[w, k + 1]])^2/vk, 
       {k, 2, n - 1}]; 
      -2^(-1)*Log[g] - n/2*Log[S/n]
];

NumericVectorQ[x_] := VectorQ[x] && And @@ NumericQ /@ x;

ToeplitzQ[x_] := SquareMatrixQ[x] && x == toe[x[[1]]];

SquareMatrixQ[x_] := MatrixQ[x] && Apply[Equal, Dimensions[x]];

(*==================Compile Functions====================*)

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

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

End[];
EndPackage[];
