Mathestate Logo

 

Zolotarev-Nolan Transform

This is a technical page to describe one of the algorithms in the StableM7 package.  Zolotarev first described this transform, but it was clarified and made useful by Nolan in 1997 (see Reference Page).  The transform can be derived step by step using Mathematica.  We start with the characteristic function in the Nolan 1-parameterization because it is algebraically simpler.  We will do the derivation only for α ≠ 1.

Zolotarev_1.gif

The Real portion of the Fourier transform of the characteristic function gives the stable density.  There is no closed formula for the general case of this integral.  Mathematica can numerically integrate some of the parameter range for this integral.

Plot[Re[NIntegrate[ Exp[-I t x] cf1/(2 Pi) /. {α -> 1.5, β -> -1}, {t, -∞, ∞}, Method -> {Automatic, "SymbolicProcessing" -> 0}]], {x, -5, 5}, PlotRange -> All, PerformanceGoal -> "Speed"]

Zolotarev_2.gif

The integral can be simplified to just use two times the portion from {t, 0, ∞} as the real portion is symmetric about t.

ig1 = 2 FullSimplify[Exp[-I t x] cf1/(2 Pi), {t > 0}]

Zolotarev_3.gif

Simplifying the real part of this, we can get a more tractable formula.  It does not perform well when x is large or α is less than about 0.7.

ig1a = FullSimplify[ComplexExpand[ig1][[1]], t > 0]
Plot[{NIntegrate[ig1a /. {α -> 1.5, β -> -1}, {t, 0, ∞}, Method -> {Automatic, "SymbolicProcessing" -> 0}], Quiet[NIntegrate[ig1a /. {α -> 0.4, β -> -0.5}, {t, 0, ∞}, Method -> {Automatic, "SymbolicProcessing" -> 0}]]}, {x, -5, 5}, PlotRange -> All, PerformanceGoal -> "Speed"]

Zolotarev_4.gif

Zolotarev_5.gif

The Zolotarev transformation transforms the parameter t of the integrand to t →z = r Zolotarev_6.gif

ig2 = ig1 /. t -> r Exp[I θ]

Zolotarev_7.gif

The above equation will need to be solved for r such that the integrand is Real.  The easiest way to do this is to remove the factor Zolotarev_8.gif then solve for r such that the exponent is real.

ig3 = PowerExpand[Log[ π ig2]]

Zolotarev_9.gif

To make some of the terms simpler we substitute -Zolotarev_10.gif→ ζ.

ig4 = ig3 /. β Tan[π α/2] -> ζ /. ζ -> -ζ

Zolotarev_11.gif

The next steps isolate the real and imaginary parts of the integrand.

ig5 = Expand[ExpToTrig[ig4]]
igi = Select[ig5, #[[1]] == i || #[[1]] == -i &]
igr = ig5 - igi

Zolotarev_12.gif

Zolotarev_13.gif

Zolotarev_14.gif

Next we solve for r so that the imaginary will be equal to zero.

rsol = Solve[-I igi == 0, r]

Zolotarev_15.gif

Zolotarev_16.gif

The exponent of the integrand then becomes

ig6 = igr /. rsol[[1]]
ig7 = FullSimplify[PowerExpand[ig6]]

Zolotarev_17.gif

Zolotarev_18.gif

Both x and r must be greater than zero for this to be real

r1 = rsol[[1]][[1]][[2]]
t0 = Solve[r1 == 0, θ];
t1 = FullSimplify[t0, {0 < α <= 2, ζ ∈ Reals}]

Zolotarev_19.gif

Zolotarev_20.gif

Zolotarev_21.gif

Using the Nolan nomenclature we will take θ0 to be Zolotarev_22.gif and ζ to be -Tan[α θ0] the domain of the integral will be  to -Pi/2 to θ0

ig8 = ig7 /. ζ -> - Tan[α θ0];
ig9 = FullSimplify[ig8]

Zolotarev_23.gif

Next we turn attention to converting Re[dt] to Re[dθ]

ExpToTrig[r Exp[I θ]]
dz1 = D[r1 Cos[θ], θ] /. ζ -> - Tan[α θ0];
dz2 = FullSimplify[dz1]

Zolotarev_24.gif

Zolotarev_25.gif

Remove the x terms

v1 = (-1 + α)/α dz2 /. x -> 1
v2 = FullSimplify[-ig9 /. x -> 1]

Zolotarev_26.gif

Zolotarev_27.gif

v1 and v2 are actually the same, but it takes a little manipulation to show this.  Removing Zolotarev_28.gif from the first part of v2 and adding it to the second term will accomplish the proof.

Zolotarev_29.gif

Zolotarev_30.gif

Zolotarev_31.gif

Zolotarev_32.gif

Zolotarev_33.gif

Zolotarev_34.gif

Now collect all the pieces and put them together.  The absolute value prevents a negative density when α <  1.  We reverse the sign of θ to put the integration in the same direction described in the Nolan paper.

Zolotarev_35.gif

Zolotarev_36.gif

Zolotarev_37.gif

The integral that we have derived is valid only for x > 0.  We get x < 0 by stable symmetry.

Zolotarev_38.gif

We haven't defined what the value is when x = 0, but it happens that we can get that from the original form.

ig1a
Integrate[ig1a /. x -> 0, {t, 0, ∞}, Assumptions -> {0 < α <= 2, -1 <= β <= 1}]

Zolotarev_39.gif

Zolotarev_40.gif

Thus

Zolotarev_41.gif

Now we have the whole thing and can show the better performance of the transformed integral, which is not troubled by infinite oscillation.

Plot[{pdf[x, 1.3, -1], pdf[x, 0.4, -0.5]}, {x, -10, 10}, PlotRange -> All, PerformanceGoal -> "Speed"]

Zolotarev_42.gif

This page is an experiment to see if there is interest in more of this level information.  If you would like to see more technical information on stable distributions, send an email.  A graphical demonstration of the integrand for the distribution function, is on the Wolfram Demonstrations Project.

Zolotarev_43.gif



© Copyright 2009 mathestate    Wed 14 Jan 2009