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: 
Functions for Fractional Gaussian Noise Modelling. For examples see the
notebook FGN-Documentation.nb.
*)

BeginPackage["FGN`", "TrenchInverse`"]

Needs["Statistics`NormalDistribution`"]; 

LoglikelihoodFGN;RegressionCoefficients;HurstCoefficient;
IterationsCount;MuHat;

FGNLL /: FGNLL::usage = "FGNLL[H,z] evaluates loglikelihood for FGN model for \
time series z. Mean zero assumed";
 
FGNAcf /: FGNAcf::usage = "FGNAcf[k,H] autocorrelation function of FGN with \
Hurst coefficient H at lag k.";
 
GetFitFGN /: GetFitFGN::usage = "GetFitFGN[z] exact mle for Hurst coefficient \
H for mean zero time series z. GetFitFGN is not normally used since it is a \
helper function for FitFGN.";

Options[FitFGN] = {MeanZeroQ -> False, MeanMLEQ -> False};
 
FitFGN /: FitFGN::usage = "FitFGN[z] exact mle for Hurst coefficient H for \
time series z. By default the series is corrected by the sample mean. If the \
mean is known to be zero, use the option MeanZeroQ->True.  If you wish to \
obtain the MLE estimate of the mean, use the option MeanMLEQ->True.";

FitRegressionFGN /: FitRegressionFGN::usage = "FitRegressionFGN[X,y] fits \
regression with FGN noise inputs are X the design matrix and y the observed \
time series";

FitRegressionFGN /: FitRegressionFGN::SingularMatrix = 
     "Matrix is singular or not p.d.";
 
FitRegressionFGN /: FitRegressionFGN::TrenchM = "TrenchM.exe returned $Failed";
 
StepIAXMatrix /: StepIAXMatrix::usage = "StepIAMatrix[n,T] sets up design \
matrix for simple step intervention with series length n and intervention =0 \
t<T and =1, t>=T.";
 
HurstK /: HurstK::usage = 
     "HurstK[z] computes Hurst K, a robust estimate of the Hurst coefficient \
H.";

SimulateGLP /: SimulateGLP::usage = "SimulateGP[c] simulates a Gaussian time \
series with autocovariance function c. If c is not a power of 2 it will be \
padded with zeros.";

SimulateFGN /: SimulateFGN::usage = "SimulateFGN[H,n] simulate a series of \
length n from FGN with parameter H. For maximum efficiency choose n to be a \
power of 2.";

Options[Acvf] = {MaxLag -> 10, MeanZero -> False};
 
Acvf /: Acvf::usage = "Acvf[z] returns the sample autocovariances for time \
series z at lags 0, 1, ... , MaxLag. It takes options MaxLag->1 and \
MeanZero->False.  If MeanZero->False, then the sample mean of z is subtracted \
and the sacvf is computed for the corrected series.";
 
Options[Acf] = {MaxLag -> 1, MeanZero -> False};
 
Acf /: Acf::usage = "Acf[z] returns the sample autocorrelations for time \
series z at lags 1, ... , MaxLag. It takes options with default values: \
MaxLag->1 and MeanZero->False.  If MeanZero->False, then the sample mean of z \
is subtracted and the sacf is computed for the corrected series. If \
MeanZero->True, no sample mean correction is made.";

ProfileLogLRegressionFGN /: ProfileLogLRegressionFGN::usage = "ProfileLogLReg\
ressionFGN[\[Alpha],X,y] computes the profile loglikelihood function, that is,
given \[Alpha] the loglikelihood maximized for H is obtained.";

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

FGNAcf=Compile[{{k, _Integer},{H, _Real}},
      0.5*(Abs[k+1.0]^(2.0*H)-2*Abs[k]^(2.0*H)+Abs[k-1.0]^(2.0*H))];

StepIAXMatrix[n_, T_] := Module[{}, Transpose[{Table[1, {n}], 
       Table[If[i < T, 0, 1], {i, n}]}]];

GetFitFGN[(z_)?NumberVectorQ] := Module[{}, 
 Check[
   FindMaximum[FGNLL[HurstCoefficient, z], {HurstCoefficient, 0.4, 
        0.8}], 
   NMaximize[{FGNLL[HurstCoefficient, z], 0 < HurstCoefficient < 
         1}, {HurstCoefficient}, Method -> "NelderMead"]
 ]
];
 
FitFGN[(z_)?NumberVectorQ, (opts___)?OptionQ] := 
 Module[{x, MQ, \[ScriptL]0 = -10^20, n = Length[z], iter = 0, Maxit = 10, 
      mn = 0, zav = 0, \[ScriptL]1, ans, h, r, y}, 
  MQ = MeanZeroQ /. {opts} /. Options[FitFGN]; 
  MMLEQ = MeanMLEQ /. {opts} /. Options[FitFGN]; 
  x = If[ !MQ || MMLEQ, z - (zav = Mean[z]), z]; 
  ans = GetFitFGN[x]; 
  h = HurstCoefficient /. Last[ans]; 
  \[ScriptL]1 = First[ans]; 
  If[MMLEQ, 
    While[iter++ < Maxit && \[ScriptL]1 - \[ScriptL]0 > 0.001, 
     \[ScriptL]0 = \[ScriptL]1; 
     r = Table[FGNAcf[k, h], {k, 0, n - 1}]; 
     mn = TrenchMean[r, x]; 
     y = x - mn; 
     ans = GetFitFGN[y]; 
     h = HurstCoefficient /. Last[ans]; 
     \[ScriptL]1 = First[ans]; 
    ]; 
  ]; 
  {\[ScriptL]1, {HurstCoefficient -> h, MuHat -> mn + zav, 
        IterationsCount -> iter}}
];

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

NumberVectorQ[x_] := VectorQ[x] && And @@ NumberQ /@ x;

NumericMatrixQ[x_] := MatrixQ[x] && And @@ NumericQ /@ Flatten[x];
 
FGNLL[(H_)?NumericQ, z_?NumericVectorQ] := Module[{r}, If[H <= 0.01 || H > 0.99, -10^10, 
      r = Table[FGNAcf[k, H], {k, 0, Length[z] - 1}]; TrenchLoglikelihood[r, 
        z]]];
 
FitRegressionFGN[(X_)?NumericMatrixQ, (y_)?NumericVectorQ] := 
  Module[{LLPrevious = -Infinity, LL = 0, MaxIter = 5, iter = 0, 
      n = Length[y], Xy, ans, H0, \[Alpha], res, r, Ginv, ac, XpXinv, 
      ETOL = 0.001, EPS = $MachineEpsilon}, 
  Xy = Transpose[Append[Transpose[X], y]]; 
  ans = Fit[Xy, {x1, x2}, {x1, x2}]; 
  \[Alpha] =  First[Transpose[({First[#1], Last[#1]} & ) /@ List @@ ans]]; 
  While[Abs[LL - LLPrevious] > ETOL && iter++ <= MaxIter, 
     LLPrevious = LL; 
     res = y - X . \[Alpha]; 
     ans = FitFGN[res, MeanZeroQ->True]; 
     LL = First[ans]; 
     H0 = HurstCoefficient /. Last[ans]; 
     r = Table[FGNAcf[k, H0], {k, 0, Length[y] - 1}]; 
     ans = TrenchM[r, EPS]; 
     If[ans === $Failed, 
         Message[FitRegressionFGN::TrenchM]; Abort[]]; 
     If[Last[ans] != 0, Message[FitRegressionFGN::SingularMatrix]; 
          Abort[]]; 
     Ginv = Partition[Most[ans], Length[y]]; 
     XpXinv = Inverse[Transpose[X] . Ginv . X]; 
     \[Alpha] = XpXinv . Transpose[X] . Ginv . y; 
   ]; 
   {LoglikelihoodFGN -> LL, RegressionCoefficients -> \[Alpha], 
       HurstCoefficient -> H0, IterationsCount -> iter}
];
 
HurstK[z_?NumericVectorQ] := Module[{n = Length[z], sd, S, R}, y = z - Mean[z]; 
      sd = Sqrt[Total[y^2]/n]; S = FoldList[Plus, 0, y]; 
      R = (Max[S] - Min[S])/sd; Log[R]/Log[n/2]];

SimulateGLP[(\[Gamma]_)?NumberVectorQ] := 
    Module[{m = Length[\[Gamma]], n, c, g, Z, Ncap}, 
     n = 2^Ceiling[Log[2, m - 1]]; acvf = If[n == m - 1, \[Gamma], 
        PadRight[\[Gamma], n + 1]]; Ncap = 2*n; 
      c = Join[acvf, Rest[Reverse[Rest[acvf]]]]; 
      g = Re[Fourier[c, FourierParameters -> {1, -1}]]; 
      Z = Table[Random[NormalDistribution[0, 1]], {Ncap - 2}]; 
      Z = (Complex[Sequence @@ #1] & ) /@ Partition[Z, 2]; 
      Z = Flatten[{Random[NormalDistribution[0, Sqrt[2]]], Z, 
         Random[NormalDistribution[0, Sqrt[2]]], Reverse[Conjugate[Z]]}]; 
      Take[Re[InverseFourier[Sqrt[g]*Z, FourierParameters -> {0, -1}]], m]/
       Sqrt[2]];

SimulateFGN[H_?NumericQ, n_?NumberQ] := Module[{\[ScriptCapitalN], ac}, 
     \[ScriptCapitalN] = 2^Ceiling[Log[2, n - 1]]; 
      ac = Table[FGNAcf[k, H], {k, 0, \[ScriptCapitalN]}]; 
      Take[SimulateGLP[ac], n]];

Acvf[(z_)?VectorQ, (opts___)?OptionQ] := Module[{ifmean0, maxlag, x, X}, 
     ifmean0 = MeanZero /. {opts} /. Options[Acvf]; 
      maxlag = MaxLag /. {opts} /. Options[Acvf]; If[maxlag <= 0, 
       Print["MaxLag<0"]; Return[{}]]; 
      X = {x0 = x = z - If[ifmean0, 0, Mean[z]]}; If[Plus @@ Abs[x] == 0, 
       Print["constant series!"]; Return[{}]]; 
      Do[AppendTo[X, x = Drop[Prepend[x, 0], -1]], {maxlag}]; 
      X . x0/Length[z]];

Acf[(z_)?VectorQ, (opts___)?OptionQ] := Module[{ifmean0, maxlag, x, X, c}, 
     ifmean0 = MeanZero /. {opts} /. Options[Acf]; 
      maxlag = MaxLag /. {opts} /. Options[Acf]; If[maxlag <= 0, 
       Print["MaxLag<0"]; Return[{}]]; 
      X = {x0 = x = z - If[ifmean0, 0, Mean[z]]}; If[Plus @@ Abs[x] == 0, 
       Print["constant series!"]; Return[{}]]; 
      Do[AppendTo[X, x = Drop[Prepend[x, 0], -1]], {maxlag}]; c = X . x0; 
      Drop[c, 1]/First[c]];

ProfileLogLRegressionFGN[\[Alpha]_?NumericVectorQ, X_?NumericMatrixQ, y_?NumericVectorQ] := 
Module[{z, ans}, 
     z = y - X . \[Alpha]; 
     First[FitFGN[z,MeanZeroQ->True]]
];

End[];
EndPackage[];
