1
$\begingroup$

I am using the Chambers-Mallows-Stuck method for simulation of skewed stable random variables http://www.sciencedirect.com/science/article/pii/0167715295001131.

There are two equations, one for $\alpha=1$ and one for leftover values of $\alpha$.

I would expect that both equation will return similar output for similar values of $\alpha$, but it is not.

I double checked the code against the paper. And I have not found any problem.

I attached the resulting histograms and the Python code at bottom.

According to the: https://en.wikipedia.org/wiki/Stable_distribution#/media/File:Levy_distributionPDF.png

It seems that my implementation (or algorithm by itself) is really wrong. Can somebody verify if the described equations in paper makes sense? (eq 3.8 and 3.9)

Case 1 $\alpha=0.99$ enter image description here

Case 2 $\alpha=1.00$ enter image description here

Case 3 $\alpha=1.01$ enter image description here

def levy_noise(n, alpha=2., beta=1., sigma=1., position=0.):
    # correct the inputs or throw error
    alpha = float(alpha)
    beta = float(beta)
    check_type_or_raise(n, int, "n")
    if not 0. <= alpha <= 2.:
        raise ValueError("Alpha must be between 0 and 2")
    if not -1. <= beta <= 1.:
        raise ValueError("Beta must be between -1 and 1")
    # generate random series
    v = np.random.uniform(-0.5*np.pi, 0.5*np.pi, n)
    w = np.random.exponential(1, n)
    if alpha == 1:
        arg1 = (0.5 * np.pi) + (beta * v)
        arg2 = w * np.cos(v)
        arg3 = 0.5 * np.pi * beta * sigma + np.log10(sigma)
        x = 0.5 * np.pi * (arg1 * np.tan(v) - beta * np.log10(arg2 / arg1))
        return (sigma * x) + arg3 + position
    else:
        arg1 = 0.5*np.pi*alpha
        b_ab = np.arctan(beta*np.tan(arg1)) / alpha
        s_ab = (1 + (beta**2) * np.tan(arg1)**2)**(1/(2.*alpha))
        arg2 = alpha * (v + b_ab)
        n1 = np.sin(arg2)
        d1 = np.cos(v)**(1/alpha)
        n2 = np.cos(v - arg2)
        x = s_ab * (n1/d1) * (n2/w)**((1-alpha)/alpha)
        return (sigma * x) + position
$\endgroup$
0

1 Answer 1

1
$\begingroup$

Typo in my code. The line:

x = 0.5 * np.pi * (arg1 * np.tan(v) - beta * np.log10(arg2 / arg1))

should be:

x = 2 / np.pi ...

Now the results seems to be ok.

The algorithm in the paper (and similar papers) is ok, the problem was on my side (this is the better case).

$\endgroup$
3
  • $\begingroup$ +1, with thanks for posting a resolution to your question -- but I'm also voting to close the thread because this solution makes the problem outside the scope of our interests on this site. $\endgroup$
    – whuber
    Commented May 16, 2017 at 13:21
  • $\begingroup$ @whuber Yes I agree with the reason. But I would like to defend the question - My answer is a prove that the algorithm is correct - that is useful information in the scope of this community. It can be really useful for somebody else in future. $\endgroup$
    – matousc
    Commented May 16, 2017 at 13:24
  • $\begingroup$ Yes--which is why it will remain visible and searchable on this site. (Indeed, I upvoted the question when it first appeared. I appreciate clear, well-formulated, well-researched posts.) $\endgroup$
    – whuber
    Commented May 16, 2017 at 13:31

Not the answer you're looking for? Browse other questions tagged or ask your own question.