Print[StyleForm[StringJoin[
        "Updated: Oct 13, 2005"
        ]
  ,FontSize->9,FontFamily->"Times", FontColor->
  RGBColor[0., 0.25,0.25]
   ]];

(* Author:  A.I. McLeod, aimcleod@uwo and Z. Krougly krougly@stats.uwo.ca*)

(* 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] 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.";

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] MathLink function: quadratic form and determinant. \
Assumes that r[1]=1. Returns vector of length 3 containing quadratic form, determinant \
and fault indicator.";

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

(*=======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.";


(********************************************************************) 
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; 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]];

(****************************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]];

End[];
EndPackage[];
