Stable Characteristic Function Fitting

If you are not well versed in complex numbers and characteristic functions, try this page first.

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 above algorithm is accurate enough that we use it for our rapid method to calculate stable γ on our volatility page.

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.