Zolotarev

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.

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"]

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}]

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"]

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

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

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 then solve for r such that the exponent is real.

ig3 = PowerExpand[Log[ π ig2]]

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

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

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

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

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

The exponent of the integrand then becomes

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

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}]

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

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

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

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

Remove the x terms

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

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

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.

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

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}]

Thus

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"]

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.