Needs["StableM`"];
SetOptions[FourierTransform,FourierParameters-> {1,1}];
SetOptions[InverseFourierTransform,FourierParameters-> {1,1}];
SetOptions[Fourier, FourierParameters -> {-1,1}];
SetOptions[InverseFourier, FourierParameters -> {-1,1}];
ReOrder::usage = "ReOrder[list] reorders a frequency domain list for the time domain and vice versa.  If Length[list] is odd the output list will be truncated to the next smallest integer.";
XTable::usage="XTable[SR, N] puts out a list of x values that correspond to the stable x-axis in the time domain for a given length, N, and SR of the DFFT list.";
DFTCF::usage= "DFTCF[cf,SR,N] puts out a discrete FT of a distribution specified by the characteric function cf in the frequency domain compatible with FourierParameters -> {-1, 1}.  SR is the sampling rate; SR is greater than 0.  N is the number of points and must be even.  The output is scaled ready for inversion.";
CF1c::usage= "CF1c[u, α, β, γ, δ] outputs a compiled version of the S(α, β, γ, δ; 1) characteristic function.  All the inputs must be real numbers.";
ReOrder=With[{n=Length[#]},Join[Drop[#,Floor[n/2]],Drop[#,-Floor[n/2]]]]&;
XTable  = Compile[{SR, {N,_Integer}},Module[{n},
n=2 Floor[N/2];
Range[-n/2,n/2-1,1]/(SR)]];
DFTCF [cf_,SR_,n_Integer]:=Module[{rlist,n2},
n2=Floor[n/2];
rlist=2. Pi SR/n Join[Range[0,n2-1],Range[-n2,-1]];
Map[cf[ #]&,rlist]SR/n
	];
CF1c=Compile[{#1,#2,#3,#4,#5},If[#1==0,1,If[#2==1,E^(I*#5*#1-#4*Abs[#1]*((2*I*#3*Log[Abs[#1]]*Sign[#1])/Pi+1)),E^(I*#1*#5-Abs[#1]^#2*#4^#2*(1-I*Sign[#1]*#3*Tan[(Pi*#2)/2]))]]];

Stable Approximation with FFT

We showed in the page on differentiation with the FFT that is was a simple matter to calculate the serial derivatives of a function in the Fourier domain.  Here we use that information to enhance the accuracy of an interpolated routine.  The interpolation function in Mathematica allows the addition of derivatives of the function; the addition will improve the accuracy of the function.  The StableM package uses the FFT of the characteristic function with the interpolation function to create a stable density approximation for maximum likelihood fitting, but the approximation does not take full advantage of the interpolation as a compromise for speed.  Below a stable characteristic function {1.5, 1, 1, 0} is sampled to create  the density and its first two derivatives.

sr=5;nn=4096;
sdft=DFTCF[CF1c[#1,1.5,1,1,0]&,sr,nn];
dens=Chop[InverseFourier[sdft ]];
g1=ListPlot[Transpose[{XTable[sr,nn],ReOrder[dens]}],PlotRange->{{-10,10},All},PlotLabel-> "Stable Density",Frame->True];
differentiate = Compile[{t} , -I t];
scale = nn/sr;
df = scale DFTCF[differentiate[#1] &, sr, nn];(*this isn't scaled so it is multiplied by scale*)
dens1 = Chop[InverseFourier[sdft df]];
g2 = ListPlot[Transpose[{XTable[sr, nn], ReOrder[dens1]}], PlotRange -> {{-10, 10}, All}, PlotLabel-> "First Derivative Stable Density",Frame->True];
dens2 = Chop[InverseFourier[sdft df df]];
g3 = ListPlot[Transpose[{XTable[sr, nn], ReOrder[dens2]}], PlotRange -> {{-10, 10}, All}, PlotLabel-> "Second Derivative Stable Density",Frame->True];
Show[g1,g2,g3,PlotLabel-> "Stable Density with 1st and 2nd Derivatives\nα = 1.5\nβ = 1"]

Below is a graph created from all of that information.

grid = Transpose[{{#} & /@ XTable[sr, nn], ReOrder[dens], ReOrder[dens1], ReOrder[dens2]}];
f = Interpolation[grid];
Plot[f[x], {x, -10, 10}, PlotRange -> All, PlotLabel -> "Stable Density\nα = 1.5\nβ =1", Frame -> True]

grid=Transpose[{{#}&/@XTable[sr, nn],ReOrder[dens],ReOrder[dens1],ReOrder[dens2]}];
f=Interpolation[grid];
Plot[f[x],{x,-10,10},PlotRange->All,PlotLabel-> "Stable Density\nα = 1.5\nβ =1",Frame->True]

The plot below shows the difference between the Zolotarev integral calculation and the third order interpolated function with and without the use of the first two derivatives.  The interpolated function with the derivatives is more accurate.

g4=Quiet[Plot[SPDF[x,{1.5,1}]-f[x],{x,-50,50},PlotRange->All,PlotLabel-> "Integral Routine Minus Interpolation with Derivatives",Frame->True]]
f1=Interpolation[Transpose[{XTable[sr, nn],ReOrder[dens]}]];
Quiet[Plot[SPDF[x,{1.5,1}]-f1[x],{x,-50,50},PlotRange->All,PlotLabel-> "Integral Routine Minus Interpolation without Derivatives",Frame->True]]

It can be seen from the first of the above graphs that the interpolated function over estimates the density.  This is because the mean of the Fourier transform was set to 1, the value of the characteristic function at t = 0.  However, our sample does not run from (-∞, ∞), but rather on the interval below.  We have calculated the value of the integral over that range.

f[[1]][[1]]
mass = Quiet[SCDF[-f[[1]][[1]][[1]], {1.5, 1}] - SCDF[f[[1]][[1]][[1]], {1.5, 1}]]

f[[1]][[1]]
mass=Quiet[SCDF[-f[[1]][[1]][[1]],{1.5,1}]-SCDF[f[[1]][[1]][[1]],{1.5,1}]]

Below is a new graph changing the first term of our sampled characteristic function in the Fourier domain to that reflect mass.  We have brought the difference closer to zero, but now are underestimating the density slightly.

sdft1=ReplacePart[sdft,1->mass/scale];
dens=Chop[InverseFourier[sdft1 ]];
dens1 = Chop[InverseFourier[sdft1 df]];
dens2 = Chop[InverseFourier[sdft1 df df]];
grid=Transpose[{{#}&/@XTable[sr, nn],ReOrder[dens],ReOrder[dens1],ReOrder[dens2]}];
f=Interpolation[grid];
g5=Quiet[Plot[SPDF[x,{1.5,1}]-f[x],{x,-50,50},PlotRange->All,PlotStyle->Darker[Red], PlotLabel-> "Integral Routine Minus Interpolation with Derivatives",Frame->True]];
Show[g4,g5,PlotRange->All]

txcolor=RGBColor[181/255,18/255,27/255];
t1 = Graphics[Text[Hyperlink[Framed[Style["Financial Home", {txcolor, Bold, Larger}]], "../map/Overview.html"], {-5, 0}]];
t2 = Graphics[Text[Hyperlink[Framed[Style["Next", {txcolor, Bold, Larger}]], "../ftr/CFReferences.html"], {0, -0.5}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../ft08/FFTInterpolation.nb"], {5, 0}]];
t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]], "../ft07/StableConvolution.html"], {0, 0.5}]];
Show[ t1, t2, t3, t4, ImageSize -> 500]
CellPrint[TextCell["\n\n© Copyright " <> ToString[DateList[][[1]]] <> " mathestate    " <> DateString[DateList[], "DateShort"], "Text"]]

© Copyright 2008 mathestate    Sat 8 Mar 2008