Kullback-Leibler Divergence

Reference

Ullah, A. (2002) Uses of entropy and divergence measures for evaluating econometric approximations and inference, Journal of Econometrics 107 313-326.

Basic Result

The discrepacy between two normal distributions with mean zero and covariance matricies (Subscript[Σ, g]) and (Subscript[Σ, h]) is given by

\(TraditionalForm\`I(g, h) = 
      1\/2\ \((\ 
          tr(Σ\_g\ Σ\_h\%\(-1\)) - 
            log(\(|\)\(Σ\_g\ Σ\_h\%\(-1\)\)\(|\)) - n)\)\)
\(TraditionalForm\`\(\(=\)\(1\/2\ \((\ 
          tr(Σ\_g\ Σ\_h\%\(-1\)) - 
            log(\(|\)\(Σ\_g\)\(|\)\(\(/\)\(\(|\)\(Σ\_h\)\(|\)\)\)) - n)\)\)\)\)

h is the approximation to g

Mathematica Code

<< FitARMA.m

<< TrenchInverse.m

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}]

MA(1) Case

<<Graphics`MultipleListPlot`

phi = {}; theta = {0.9}; ac = TacvfARMA[phi, theta, 1000];
rs = Range[10, 50]; ns = {100, 200, 500, 1000};

Timing[tb = 
          Table[{n = ns[[i]], 
              r = rs[[j]], 
              KLDiscrepancy[ac, r, n]}, {i, Length[ns]}, {j, 
              Length[rs]}];]

tbInv = 
        Table[{n = ns[[i]], 
            r = rs[[j]], 
            KLDiscrepancyInverted[ac, r, n, phi, theta]}, {i, 
            Length[ns]}, {j, Length[rs]}];

Dimensions[tb2 = Transpose[tb, {2, 3, 1}]]

ListPlot[
        Rest /@ tb[[1]]];

Multiple List Plot

g = 
        Table[MultipleListPlot[{Rest /@ 
                tb[[i]], 
              Rest /@ tbInv[[i]]}, 
            PlotLabel -> 
              StringForm["n = ``", 
                ns[[i]]], 
            Axes -> False, Frame -> True, 
            FrameLabel -> {"r", "ℐ"}, 
            RotateLabel -> False], {i, Length[ns]}];

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

List Plot

g = 
        Table[ListPlot[Rest /@ tb[[i]], 
            PlotLabel -> 
              StringForm["n = ``", 
                ns[[i]]], 
            Axes -> False, Frame -> True, 
            PlotStyle -> {PointSize[0.03]}, 
            FrameLabel -> {"r", "ℐ"}, 
            RotateLabel -> False], {i, Length[ns]}];

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

Graphic For Paper

g = 
        ListPlot[Rest /@ tb[[2]], 
          Axes -> False, Frame -> True, 
          PlotStyle -> {PointSize[0.02]}, 
          FrameLabel -> {"r", "ℐ"}, 
          RotateLabel -> False];

Export["d:/math/2007/faster/kl/KL.eps", g]

Export["e:/tex/2007/faster/KL.eps", g]

Table

tb200 = 
      Rest /@ tb[[2]][[{1, 11, 21, 31, 41}]]

tbInv200 = 
      Rest /@ tbInv[[2]][[{1, 11, 21, 31, 41}]]

{x, y1} = Transpose[tb200]; {x, y2} = 
      Transpose[tbInv200];
\(xyz = {x, y1, y2} // Transpose; 
    TableForm[xyz, 
      TableHeadings -> {None, {"n", "MMSE Predictor", "Inverted Form"}}]\)

xyz = {x, ToString[PaddedForm[#, {3, 2}]] & /@ 
            y1, ToString[PaddedForm[#, {3, 2}]] & /@ y2} // 
        Transpose

\(TableForm[xyz, 
      TableHeadings -> {None, {"n", "MMSE Predictor", "Inverted Form"}}]\)

Export["d:/math/2007/faster/kl/KLTab.dat", xyz]