Stable Characteristic Function Fitting

Below is the stable characteristic function in the 1-parameterization.

This can be rewritten in polar coordinate form:

where

and

If we take the absolute value of a result for the stable cf, we have the value of r.

At the point t = ± 1, and when γ = 1,  r = 1/e .

This fact can be used with the sample characteristic function to find γ from a stable sample.  The primary assumption is that the sample characteristic function converges at t = ± 1.

The empirical characteristic function is defined below.

In Mathematica the empirical characteristic function is conveniently given by the formula below, where s is a list of identically distributed random variables:

Here is the algorithm we use to estimate γ from a single point on the characteristic function is below.  We need a starting point for the search and the QuartileDeviation is a good one when α > 1.

s = SRandom[1000, {1.3, 0.5, 0.345, 1}];
c0 = QuartileDeviation[s]; cc /. FindRoot[Abs[Mean[Exp[I s/cc]]] - 1/E, {cc, c0}]

The empirical characteristic function can be used to estimate all four stable parameters.  Below is the algorithm that is used in StableM to compute stable parameters.  The routine is a little different than published algorithms.  It is in the 0-paramteterization where all four parameters are continuous and uses the above algorithm and also takes advantage of the fact the δ0 can similarly be estimated at t = 1 when γ = 1.  This gives starting points for a more detailed parameter search without other methods such as the quantile method.

SConvert[{a_, b_, c_, d_}, input_Integer, output_Integer] := Which[input == output, {a, b, c, d},
   input == 0 && output == 1, {a, b, c, If[a == 1, d - b 2/Pi c Log[c], d - b c Tan[Pi a/2]]},
   input == 1 && output == 0, {a, b, c, If[a == 1, d + b 2/Pi c Log[c], d + b c Tan[Pi a/2]]}];

(*sample characteristic function compiled version*)
cfs = Compile[{t, {s, _Real, 1}}, Mean[Exp[I t s]]];
(*sample cf with test to prevent symbolic processiing*)
CFs[(t_)?NumericQ, s_] := cfs[t, s];

CFFit[s_, type_Integer: 1] := Module[{soln, cc, s1, c, c0, c1, c2, i, t, a, d, d0, d1, d2, y, b, tval, cfval},
   (*get c0*)
   c0 = QuartileDeviation[s];
   c = Abs[FindRoot[Log[Abs[CFs[1/cc, s] ]] + 1, {cc, c0}][[1]][[2]]];
   (*get d0*)
   d = Median[s];(*use median to move d closer to zero*)
   d1 = Arg[CFs[1, (s - d)/c]];(*doesn't work unless d close to zero*)
   d = d + c  d1;
   c0 = c;
   tval = Range[32]/32;
   
   c1 = 1; c2 = 2;
   i = 0;
   While[Abs[c2 - 1] > 10^-8 && i < 20, i++;
    c = c c1;
    cfval = {Log[#], Log[-Log[Chop[cfs[#, (s - d)/c] cfs[-#, (s - d)/c]]]]} & /@ tval;
    y = Fit[cfval, {t, 1}, t];
    a = y[[2]][[1]];
    
    c1 = Exp[(y[[1]] - Log[2])/a];
    c2 = c1;
    ];
   a = If[a >= 2, 2, a];
   
   d0 = d;
   
   d1 = 0; d2 = 1;(*allows entry into while and preserves correct relation of b and d*)
   While[Abs[d2] > 10.^-8,
    d = d + c  d1;
    cfval = {#, Arg[cfs[#, (s - d)/c]]/(# )} & /@ tval;
    y = Fit[cfval, {t^(a - 1), 1}, t];
    If[Abs[ y[[2]][[1]]/ Tan[Pi a/2]] > 1, Break[]];
    d1 = y[[2]][[1]] + y[[1]];
    d2 = d1];
   
   b = y[[2]][[1]]/ Tan[Pi a/2];
   b = If[a >= 2, 0, b];
   b = If[Abs[b] > 1, Sign[b], b];
   If[type == 1, SConvert[{a, b, c, d}, 0, 1], {a, b, c, d}]
   ];

CFFit[s]

Characteristic function fits improve in accuracy as sample size increases but not so quickly as does the maximum likelihood fit.  The fit for γ based on only one point at t = 1, is less accurate than the full characteristic function fit, but it is much faster.  The full characteristic function fit is likewise much faster than the maximum likelihood fit.

© Copyright 2008 mathestate    Fri 29 Feb 2008