This notebook demonstrates two methods of fitting stable parameters to a sample.  Stock market price data is used as an example with the data contained in the file wmt.csv.  The quantile method is shown, but the main use of this is to determine starting parameters for a maximum likelihood algorithm based on a FFT of the characteristic function of the stable distribution to quickly determine an approximate PDF function for given parameters.  The algorithms are found in the FFTFit.m file, comments and suggestions for improving the algorithm would be appreciated.
Bob Rimmer
rrimmer@npgcable.com

Change file paths appropriate to your computer via Input/Get File Path.

<< FFTFit.m

StableFunctions.m available at: http://library.wolfram.com/infocenter/MathSource/4377/

<< StableFunctions.m

wmt is a list of four years of daily closing prices for the stock Walmart, ending 5/24/2004, loaded from the file wmt.csv.  Similar more comprehensive data can be downloaded from: "http://finance.yahoo.com/q/hp?s=GE"

wmt = Flatten[Import["C:........wmt.csv"]];

Some routines:  LogReturn generates a seqence of daily logarithmic returns from a daily price time series.  ScaleS1 scales a set of S1 stable parameters by n.

<<  "Graphics`Graphics`"
Off[General::"spell1"];
Difference[
          data_] :=   Drop[data, 1] -  
          Drop[data, -1];
LogReturn[data_]  :=  
        Difference[Log[N[data]]];
ScaleS1[n_, 
          s1params_List]  :=  {s1params[[1]], s1params[[2]], 
          n^(1/s1params[[1]]) s1params[[3]], 
          n s1params[[4]]};
\(EmpiricCDFList[
          data_]\  := \ 
        Block[{prob, sel, 
            del}, 
					
					prob = \ Table[\((i - 0.5)\)/Length[data], {i, 
                Length[data]}]; 
					Transpose[
            Join[{Sort[data]}, {prob}]]];\)

A plot of WMT stock price.

ListPlot[wmt, PlotJoined -> True];

A plot of daily log returns for WMT.

lr =  LogReturn[wmt];
ListPlot[lr, PlotJoined -> True];

A quantile fit for stable parameters: {α, β, γ, δ} in the Nolan 1 parameterization.  See StableFunctionsInformation.nb from above Wolfram website for more information.

qft = SQFit[lr, 1]

The maximum likelihood fit for the same data.

mft = fMLFit[lr, 1]

Loglikelihood for the quanitle fit.

Plus @@ Log[Map[SPDF[#, qft, 1] &, lr]]

Loglikelihood for the FFT Maxlikelihood fit. It is greater than the quantile fit: the fMLFit is better.

Plus @@ Log[Map[SPDF[#, mft, 1] &, lr]]

A plot of the fit to the empirical distribution.

\(\(\(
\)\(\(p1\  = \ 
        ListPlot[EmpiricCDFList[lr], PlotStyle -> Hue[0.15], 
          PlotRange -> All, DisplayFunction -> Identity];\)

    \(p2\  = \ 
        Plot[SCDF[x, mft, 1], {x, Min[lr], Max[lr]}, 
          PlotStyle -> {Hue[0.67], Thickness[0.003]}, 
          DisplayFunction -> Identity];\)

    \(Show[p1, p2, DisplayFunction -> $DisplayFunction, 
        PlotLabel -> \ "Empirical Distribution Fit - WMT Close Data"];\)\)\)\)

A log - log plot of the fit.

lrlist = EmpiricCDFList[lr];
ltlr  =  Select[lrlist, #[[1]] < 0 &];
rtlr =  Select[lrlist, #[[1]] > 0 &];
lp1 = 
        LogLogListPlot[Abs[ltlr], PlotStyle -> RGBColor[0, 0, 1], 
          PlotRange -> All, 
          DisplayFunction -> Identity];
trt = Transpose[rtlr];
lp2 =  
        LogLogListPlot[Transpose[{trt[[1]], 1 -  trt[[2]]}], 
          PlotStyle -> RGBColor[1, 0, 0], PlotRange -> All, 
          DisplayFunction -> Identity];
lp3 =  
        LogLogPlot[SCDF[-x, mft, 1], {x, 0.00001, -4 Min[lr]}, 
          PlotStyle -> {RGBColor[ .5,  .5, 1], Thickness[0.003]}, 
          PlotRange -> All, 
          DisplayFunction -> Identity];
lp4 =  
        LogLogPlot[1 - SCDF[x, mft, 1], {x, 0.00001, 4 Max[lr]}, 
          PlotStyle -> {RGBColor[1,  .5, 0.5], Thickness[0.003]}, 
          PlotRange -> All, 
          DisplayFunction -> Identity];
q1 = 
        Show[lp1, lp2, lp3, lp4, 
          PlotLabel -> "Log Log Left Tail Blue, Right Red", 
          AspectRatio -> Automatic, PlotRange -> All, 
          DisplayFunction -> $DisplayFunction];

How it works: 
The algorithm is modified from Rachev & Mittnick, Stable Paretian Models in Finance,  p. 121 - 122.
 twoexp is the power of 2 that determines the number of sampling points used.
 Mathematica can do Fourier transforms of lists of varying lengths, but it is still more efficient to use powers of 2.
 xrange is the x - axis interval over which the function is transformed, scaled to γ = 1. 
 The length of the period of the stable distribution, of course, is infinite, but if a range is selected that contains nearly all of the probability density, the calculations can be very accurate.  
 In this algorithm the xrange or x--period length is roughly centered about zero. 
 It corresponds to 2Pi/h in Rachev's algorithm, where h is the fundamental angular frequency.  xrange must be large enough to capture essentially all the probability of the distribution for the parameters used.
  When α ~ 1, both the xrange and the twoexp need to be increased progressively as α -> 0 to maintain accuracy. 
  Accuracy goal is the parameter fed to FindMaximum.

Options[fMLFit]

AccuracyGoal -> 4 usually strikes a good compromise between speed and accuracy.

mft = fMLFit[lr, 1]; // AbsoluteTiming
mft
Plus @@ Log[Map[SPDF[#, mft, 1] &, lr]]

mft1 = fMLFit[lr, 1, AccuracyGoal -> 6]; // 
      AbsoluteTiming
mft1
Plus @@ Log[Map[SPDF[#, mft1, 1] &, lr]]

twoexp determines the number of data points used in the FFT: (2^twoexp).  Increasing it will take more time but increase accuracy.

mft = fMLFit[lr, 1, twoexp -> 14]; // 
      AbsoluteTiming
2^14
mft
Plus @@ Log[Map[SPDF[#, mft, 1] &, lr]]

xrange is the range on the x-scale of a standardized stable distribution (γ=1, δ=0).  The actual range is {-xrange/2, xrange (1/2 -   (2. ^(-twoexp)))}, because an even number of points are used.  In general xrange should not be more than  (2^(twoexp - 1)).

{-xrange/2, 
        xrange (1/2 -   (2. ^(-twoexp)))} /. {xrange -> 1024, 
        twoexp -> 13}

mft = fMLFit[lr, 1, twoexp -> 10, xrange -> 200]; // 
      AbsoluteTiming
2^10
mft
Plus @@ Log[Map[SPDF[#, mft, 1] &, lr]]

fMLFit1 is an adaptive routine under development that attempts to strike a compromise between speed and accuracy based on the sample.

mft = fMLFit1[lr, 1]; // AbsoluteTiming
mft
Plus @@ Log[Map[SPDF[#, mft, 1] &, lr]]