Mathematica Code

<< FitARMA.m

<< TrenchInverse.m

?TacvfFDWN

KLDiscrepancy::usage = "KLDiscrepancy[c,r,n] KL-discrepancy for \
series of length n with autocovariance function c ";

KLDiscrepancy::"invalid"  =  "n must be < Length[ac]. n reset \
to ``.";

KLDiscrepancy[ac_, r_, n_] := 
        Module[{c, 𝓃, g, phiApprox, 
            h}, 
𝓃 = 
            If[n > Length[ac], 
              Message[KLDiscrepancy::"invalid", Length[ac]]; Length[ac], 
              n]; 
g = 
            Toeplitz[
              c = Take[ac, 𝓃]]; 
phiApprox = 
            ARCoefficients /. 
              PacfDL[Take[c, r],  
                LinearPredictor -> True]; 
h = 
            Inverse[Toeplitz[
                TacvfAR[
                  phiApprox, 𝓃 - 1]]]; 
0.5*(Tr[
                  g . h] - Log[Det[g]*Det[h]] - 
                Length[g])
];

KLDiscrepancyInverted[ac_, r_, n_,  phi_,  theta_] := 
        Module[{c, 𝓃, g, phiapprox, 
            h}, 
𝓃 = 
            If[n > Length[ac], 
              Message[KLDiscrepancy::"invalid", Length[ac]]; Length[ac], 
              n]; 
g = 
            Toeplitz[
              c = Take[ac, 𝓃]]; 
phiapprox = 
            Rest[InvertedCoefficientsARMA[phi, theta, 
                r]]; 
h = 
            Inverse[Toeplitz[
                TacvfAR[
                  phiapprox, 𝓃 - 1]]]; 
0.5*(Tr[
                  g . h] - Log[Det[g]*Det[h]] - 
                Length[g])
];

ac = TacvfARMA[{}, {0.95}, 1000];

KLDiscrepancy[ac, 20, 200]

KLDiscrepancyInverted[ac, 20, 200, {}, {0.95}]

FDWN Case

ac = TacvfFDWN[0.4, 1000];

KLDiscrepancy[ac, 20, 200]

n = 200; rs = Range[10, 50]; ds = {0.1, 0.2, 0.3, 0.4};

Timing[tb = 
          Table[{d = ds[[i]], 
              r = rs[[j]], 
              ac = TacvfFDWN[ds[[i]], 
                  1000]; KLDiscrepancy[ac, r, n]}, {i, Length[ds]}, {j, 
              Length[rs]}];]

Round[
        Max[Transpose[tb[[1]]] // 
              Last]*100]/100.0

g = 
        Table[ymax = 
            Round[0.8*
                  Max[Transpose[
                        tb[[i]]] // 
                      Last]*100]/100.0; 
ymin = 
            Round[1.2*
                  Min[Transpose[
                        tb[[i]]] // 
                      Last]*100]/100.0; 
ListPlot[
            Rest /@ tb[[i]], 
            PlotLabel -> 
              StringForm["d = ``", 
                PaddedForm[ds[[i]]], {3, 
                  1}], Axes -> False, Frame -> True, 
            PlotStyle -> {PointSize[0.03]}, 
            FrameLabel -> {"r", "ℐ"}, 
            RotateLabel -> False, PlotRange -> {Automatic, {0, 1}}, 
            FrameTicks -> {{10, 30, 50}, {0.1, 0.9}, None, None}], {i, 
            Length[ds]}];

ga = Show[GraphicsArray[Partition[g, 2]]];

Export["d:/math/2007/faster/kl/KLFDWN.eps", ga]

Export["e:/tex/2007/faster/KLFDWN.eps", ga]