Stable PRNG

(*Algorithms for the stable random number generator.*)
Z=Function[{α, β, u1, u2}, Module[{W,θ, t0},
t0 = ArcTan[β Tan[((π α)/2)]]/α;
θ=  Pi(u1-1/2);
W=- Log[u2];
Sin[α(t0+θ)](Cos[α t0] Cos[θ])^(-1/α) (Cos[α t0 + θ(α-1)]/W)^((1-α)/α)]];

Zc=Compile[{α, β, u1, u2}, Module[{W,θ, at0},
at0 = ArcTan[β Tan[((π α)/2)]];
θ=  Pi(u1-1/2);
W=- Log[u2];
Sin[at0+α θ](Cos[at0] Cos[θ])^(-1/α) (Cos[at0 + θ(α-1)]/W)^((1-α)/α)],CompilationTarget->"C",Parallelization->True,RuntimeAttributes->{Listable}];

Iz=Compile[{x,α, β, u1, u2}, Module[{W,θ, at0},
at0 = ArcTan[β Tan[((π α)/2)]];
θ=  Pi(u1-1/2);
W=- Log[u2];
Boole[Sin[at0+α θ](Cos[at0] Cos[θ])^(-1/α) (Cos[at0 + θ(α-1)]/W)^((1-α)/α)≤ x]],CompilationTarget->"C",Parallelization->True,RuntimeAttributes->{Listable}];

IZ[x_,a_,b_,u1_?NumericQ,u2_?NumericQ]:=Iz[x,a,b,u1,u2];

theta0 = ArcTan[#2 Tan[Pi #1/2]]/#1 &;

Z1=Function[{b,u1,u2},Block[{w, theta},
theta=π (-(1/2)+u1);
 w = -Log[u2]; 
   (1/π)(-2 b Log[((π w Cos[theta])/(π+2 b theta))]+(π+2 b theta) Tan[theta]) ]];

The stable pseudorandom number generator provides some insight into the behavior and structure of stable random variables and contains the necessary formulas to compute stable probabilities.  This page will demonstrate the algorithm when α ≠ 1, for the standardized 1-parameterization.  The basic algorithm in Mathematica is shown below.

Z=Function[{α,β,u1,u2},Module[{W,θ,t0},t0=ArcTan[β Tan[(π α)/2]]/α;
θ=Pi(u1-1/2);
W=-Log[u2];
Sin[α(t0+θ)](Cos[α t0] Cos[θ])^(-1/α) (Cos[α t0+θ(α-1)]/W)^((1-α)/α)]];

u1 and u2 are independent uniform random variables with a range (0, 1).  As such this function can also be thought of as a two dimensional inverse quantile function for stable distributions spread over a unit square.  It also shows that stable random variables can arise as the product of two random variables.  This idea will be explored in more detail below.

The 3D plot below shows the function Z for α = 1.5, β = 0.5 in green.  A plane set at Z = 1 is shown in blue.  The x and y axes of the graph range from zero to one, thus the blue plane has an area of one.  The plane cuts the green plot in such a manner that the area of the square to the left of the cut is the probability of Z ≤ 1.  The contour of the cut is the Zolotarev integrand for the stable distribution function scaled to have its domain constrained by the unit square.

{aa,bb}={1.5,-0.5};
Plot3D[{1,Zc[aa,bb,x,y]},{x,0,1},{y,0,1},PlotStyle->{{Blue,Opacity[0.5]},Green},Exclusions->{x==0,y==0,x==1,y==1},ClippingStyle->None, PlotLabel-> "Z( 1.5, 0.5, x, y)"]

The curves from the integrand are shown as contours on the unit square in the graph below.  The stable probability function could be calculated as the area to the left of the curve using an indicator function, (Subscript[I, z]), which returns 1 if Z ≤ z and 0 otherwise.

(∫_0^1)(∫_0^1)(Subscript[I, z])(α,β,x,y)dydx

But usually it is easier and more accurate to calculate the area under the curve and subtract from 1.

Note also that there is one vertical curve in the group, this is at Z = 0, where there is a known function for the distribution function, for α ≠ 1.

F(0)=(1/2)-(((tan^(-1))(β tan(((π α)/2))))/(π α))

The lower portions of the curves converge to this vertical point, when α > 1.

(lim_(y->0)) Z(α,β,(1/2)-(((tan^(-1))(β tan(((π α)/2))))/(π α)),y)=0

ContourPlot[Zc[aa,bb,x,y],{x,0,1},{y,0,1},Exclusions->{x==0,y==0,x==1,y==1},Contours->50,ContourShading->False,PlotLabel-> "Contour Plot\nZ( 1.5, 0.5, x, y)"]

When α < 1, the convergence to the line is in the other direction.

(lim_(y->1)) Z(α,β,(1/2)-(((tan^(-1))(β tan(((π α)/2))))/(π α)),y)=0

{aa,bb}={0.5,0.5};
ContourPlot[Zc[aa,bb,x,y],{x,0,1},{y,0,1},Exclusions->{x==0,y==0,x==1,y==1},Contours->50,ContourShading->False,PlotLabel-> "Contour Plot\nZ( 0.5, 0.5, x, y)"]

The formula can be rewritten so that Z is the product of two random variables, each based on an independent uniform random variable, q and p.  The separate formulas are quantile functions for each random variable.

Z=(((-(1/(log(q)))))^((1-α)/α)) (sin(α ((((tan^(-1))(β tan(((π α)/2))))/α)+π (p-(1/2)))) (cos^((1-α)/α))((tan^(-1))(β tan(((π α)/2)))+π (α-1) (p-(1/2))) (((((cos(π (p-(1/2))))/(Sqrt[(β^2) (tan^2)(((π α)/2))+1]))))^(-1/α)))

The formula for the q component can be solved to a distribution function.

W(x,α)=({{{, ({{1-(e^(-(x^(α/(-1+α))))), x>0&&α>1}, {(e^(-(x^(α/(-1+α))))), x>0&&α<1}, {0, True}})}})

This distribution function is a variant of the Weibull and Fréchet extreme value distributions.  A plot is shown below for curves where α = {0.4, 0.8, 1.2, 1.6, 2} from bottom to top.

W[x_,α_]:=Piecewise[{{1-(E^(-(x^(α/(-1+α))))),x>0&&α>1},{(E^(-(x^(α/(-1+α))))),x>0&&α<1}},0]
Quiet@Plot[Evaluate[W[x,#]&/@(Rest[Range[0,20,4]]/20)],{x,0,5},PlotLabel-> "W(x, α)\nα = {0.4, 0.8, 1.2, 1.6, 2.0}"]

The distribution function for the other component, p, cannot be solved, but it can be plotted as an inverse quantile function.

(V^(-1))(p,α,β) = sin((tan^(-1))(β tan(((π α)/2)))+π α (p-(1/2))) (cos^((1-α)/α))((tan^(-1))(β tan(((π α)/2)))+π (α-1) (p-(1/2))) (((((cos(π (p-(1/2))))/(Sqrt[(β^2) (tan^2)(((π α)/2))+1]))))^(-1/α))

Plots for this distribution function are shown below for the symmetric case and the same α’s.  This component of the random variable determines the sign of Z and most of the tail behavior.  For the normal stable case when α = 2 the distribution is finite and has a formula.

InverseV[p_,α_,β_]:=((Cos[(-(1/2)+p) π (-1+α)+ArcTan[β Tan[((π α)/2)]]])^((1-α)/α)) Sin[(-(1/2)+p) π α+ArcTan[β Tan[((π α)/2)]]] (((((Cos[(-(1/2)+p) π])/(Sqrt[1+(β^2) ((Tan[((π α)/2)])^2)]))))^(-1/α));
g1=ParametricPlot[{InverseV[p,0.4,0],p},{p,0.001,0.999},PlotRange->{{-5,5},All},AspectRatio->1/GoldenRatio];
g2=ParametricPlot[{InverseV[p,0.8,0],p},{p,0.001,0.999},PlotRange->{{-5,5},All},AspectRatio->1/GoldenRatio];
g3=ParametricPlot[{InverseV[p,1.2,0],p},{p,0.001,0.999},PlotRange->{{-5,5},All},AspectRatio->1/GoldenRatio];
g4=ParametricPlot[{InverseV[p,1.6,0],p},{p,0.001,0.999},PlotRange->{{-5,5},All},AspectRatio->1/GoldenRatio];
g5=ParametricPlot[{InverseV[p,2,0],p},{p,0.001,0.999},PlotRange->{{-5,5},All},AspectRatio->1/GoldenRatio];
Show[g1,g2,g3,g4,g5,PlotLabel -> "V(x, α, 0)\nα = {0.4, 0.8, 1.2, 1.6, 2.0}"]

For the normal stable case when α = 2, the distribution is finite and has a formula.

V(x,2,0)=(((cos^(-1))(-(x/2)))/π)

In that case the tail behavior is determined by W (x, 2).

W(x,2)= ({{{, ({{1-(e^(-(x^2))), x>0}, {0, True}})}})

The loglog plot below shows the Pareto tail behavior of the function, V, when α < 2.  Calculations of this show that the negative slope of the loglog plot matches α.

FV[x_,a_,b_]:=FindRoot[InverseV[p,a,b]==x,{p,0.5,0,1}][[1]][[2]];
g1=Quiet@LogLogPlot[{FV[-x,0.5,0.5],1-FV[x,0.5,0.5]},{x,0.1,1000}];
g2=Quiet@LogLogPlot[{FV[-x,1.5,0.5],1-FV[x,1.5,0.5]},{x,0.1,500}];
Show[g1,g2, PlotLabel-> "LogLog Distribution plots showing tail behavior\nα = 0.5, α = 1.5 with β = 0.5\nRight Tail Red"]

We have shown the operation of the stable pseudo random generator and that stable random variables are generated as a product of two random variables, one belonging to a light tailed extreme value type of distribution with a tail exponent of (α/(α-1)) and the other of a peculiar distribution derived from the Zolotarev transform of the Fourier transform of the stable characteristic function.  This latter distribution determines the sign and tail exponent of the stable random variable.

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

© Copyright 2011 mathestate    Fri 17 Jun 2011