A distribuição estável positiva em R

As distribuições estáveis ​​positivas são descritas por quatro parâmetros: o parâmetro skewness , o parâmetro de escala , o parâmetro de localização , e assim parâmetro de índice chamado . Quando é zero, a distribuição é simétrica em torno de \ mu , quando é positivo (resp. negativo) a distribuição é inclinada para a direita (resp. para a esquerda) Distribuições estáveis ​​permitem caudas de gordura quando \ alpha diminui.σ > 0 μ ∈ ( - ∞ , ∞ ) α ∈ ( 0 , 2 ] β μ $\beta\in[-1,1]$$\sigma>0$$\mu\in(-\infty,\infty)$$\alpha\in(0,2]$$\beta$$\mu$$\alpha$

Quando $\alpha$ é estritamente menor que um e $\beta=1$ o suporte da distribuição se restringe a $(\mu,\infty)$ .

A função de densidade possui apenas uma expressão de forma fechada para algumas combinações particulares de valores para os parâmetros. Quando $\mu=0$ , $\alpha<1$ , $\beta=1$ e $\sigma=\alpha$ for (consulte a fórmula (4.4) aqui ):

$f(y) = -\frac{1}{\pi y} \sum_{k=1}^{\infty} \frac{\Gamma(k\alpha+1)}{k!} (-y^{-\alpha})^k \sin(\alpha k \pi)$

Tem infinita média e variância.

Questão

Eu gostaria de usar essa densidade em R. eu uso

> alpha <- ...
> dstable(y, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)

onde a função dstable vem com o pacote fBasics.

Você pode confirmar que este é o caminho certo para calcular essa densidade em R?

Agradeço antecipadamente!

EDITAR

Uma razão pela qual desconfio é que, na saída, o valor do delta é diferente daquele na entrada. Exemplo:

> library(fBasics)
> alpha <- 0.4
> dstable(4, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)
[1] 0.02700602
attr(,"control")
   dist alpha beta gamma    delta pm
stable   0.4    1   0.4 0.290617  1

A resposta curta é que seu é bom, mas o seu é errado. Para obter a distribuição estável positiva dada pela sua fórmula em R, você precisa definir γ γ = | 1 - i tan ( π α / 2 ) | - 1 / a$\delta$$\gamma$

$$\gamma = |1 - i \tan \left(\pi \alpha / 2\right)|^{-1/\alpha}.$$

O exemplo mais antigo que pude encontrar da fórmula que você deu foi (Feller, 1971), mas só encontrei esse livro em forma física. Contudo (Hougaard, 1986) fornece a mesma fórmula, juntamente com a transformação de Laplace No manual ( usado em ), a parametrização é de (Samorodnitsky e Taqqu, 1994), outro recurso cuja reprodução on-line me escapou. Entretanto (Weron, 2001) fornece a função característica na parametrização de Samorodnitsky e Taqqu para como

$$\mathrm{L}(s) = \mathrm{E}\left[\exp(-sX)\right] = \exp\left(-s^\alpha\right).$$ stablediststabledistfBasicspm=1$\alpha \neq 1$ $$\varphi(t) = \mathrm{E}\left[\exp(i t X) \right] = \exp\left[i \delta t - \gamma^\alpha |t|^\alpha \left(1 - i \beta \mathrm{sign}(t) \tan{\frac{\pi \alpha}{2}} \right) \right].$$ Renomeei alguns parâmetros do artigo de Weron para coincidir com a notação que estamos usando. Ele usa para e para . De qualquer forma, ao conectar e , obtemos $\mu$$\delta$$\sigma$$\gamma$$\beta=1$$\delta=0$ $$\varphi(t) = \exp\left[-\gamma^\alpha |t|^\alpha \left(1 - i \mathrm{sign}(t) \tan \frac{\pi \alpha}{2} \right) \right].$$

Observe que para e que . Formalmente, , portanto, definindo em obtemos Um ponto interessante a ser observado é que que corresponde a também é , portanto, se você tentar ou$(1 - i \tan (\pi \alpha / 2)) / |1 - i \tan(\pi \alpha / 2)| = \exp(-i \pi \alpha / 2)$$\alpha \in (0, 1)$$i^\alpha = \exp(i \pi \alpha / 2)$γ = | 1 - i tan ( π α / 2 ) | - 1 / α φ ( t ) φ ( i s ) = exp ( - s α ) = L ( s ) . γ ct = 1 / 2 1 / 2 γ = ct γ = 1 - ct $\mathrm{L}(s)=\varphi(is)$$\gamma = |1 - i \tan \left(\pi \alpha / 2\right)|^{-1/\alpha}$$\varphi(t)$

$$\varphi(is) = \exp\left(-s^\alpha\right) = \mathrm{L}(s).$$ $\gamma$$\alpha=1/2$$1/2$$\gamma=\alpha$$\gamma=1-\alpha$, que na verdade não é uma aproximação ruim, você acaba exatamente correto para .$\alpha=1/2$

Aqui está um exemplo em R para verificar a correção:

library(stabledist)

# Series representation of the density
PSf <- function(x, alpha, K) {
  k <- 1:K
  return(
    -1 / (pi * x) * sum(
      gamma(k * alpha + 1) / factorial(k) * 
        (-x ^ (-alpha)) ^ k * sin(alpha * k * pi)
    )
  )
}

# Derived expression for gamma
g <- function(a) {
  iu <- complex(real=0, imaginary=1)
  return(abs(1 - iu * tan(pi * a / 2)) ^ (-1 / a))
}

x=(1:100)/100
plot(0, xlim=c(0, 1), ylim=c(0, 2), pch='', 
     xlab='x', ylab='f(x)', main="Density Comparison")
legend('topright', legend=c('Series', 'gamma=g(alpha)'),
       lty=c(1, 2), col=c('gray', 'black'),
       lwd=c(5, 2))
text(x=c(0.1, 0.25, 0.7), y=c(1.4, 1.1, 0.7), 
     labels=c(expression(paste(alpha, " = 0.4")),
              expression(paste(alpha, " = 0.5")),
              expression(paste(alpha, " = 0.6"))))

for(a in seq(0.4, 0.6, by=0.1)) {
  y <- vapply(x, PSf, FUN.VALUE=1, alpha=a, K=100)
  lines(x, y, col="gray", lwd=5, lty=1)
  lines(x, dstable(x, alpha=a, beta=1, gamma=g(a), delta=0, pm=1), 
        col="black", lwd=2, lty=2)
}

$\hskip1in$

1. Feller, W. (1971). Uma Introdução à Teoria da Probabilidade e Suas Aplicações , 2 , 2ª ed. Nova York: Wiley.
2. Hougaard, P. (1986). Modelos de sobrevivência para populações heterogêneas derivadas de distribuições estáveis , Biometrika 73 , 387-396.
3. Samorodnitsky, G., Taqqu, MS (1994). Processos aleatórios não gaussianos estáveis , Chapman & Hall, Nova York, 1994.
4. Weron, R. (2001). Distribuições de Levy-estável revisitadas: índice de cauda> 2 não exclui o regime de Levy-estável , International Journal of Modern Physics C, 2001, 12 (2), 209-223.

O que eu acho que está acontecendo é que na saída deltapode estar relatando um valor interno de localização, enquanto na entrada deltaestá descrevendo a mudança. [Parece haver um problema semelhante com gammaquando pm=2.] Portanto, se você tentar aumentar o turno para 2

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1)
[1] 0.06569375
attr(,"control")
   dist alpha beta gamma    delta pm
 stable   0.4    1   0.4 2.290617  1

então você adiciona 2 ao valor do local.

Com beta=1e pm=1você tem uma variável aleatória positiva com um limite inferior de distribuição em 0.

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=1))
[1] 0.002666507

Alterne 2 e o limite inferior aumenta na mesma quantidade

> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=2, pm=1))
[1] 2.003286

Mas se você deseja que a deltaentrada seja o valor do local interno em vez do deslocamento ou limite inferior, será necessário usar uma especificação diferente para os parâmetros. Por exemplo, se você tentar o seguinte (com pm=3e tentando delta=0e o delta=0.290617que encontrou anteriormente), parece que entra deltae sai da mesma forma . Com pm=3e delta=0.290617você obtém a mesma densidade de 0,02700602 encontrada anteriormente e um limite inferior em 0. Com pm=3e delta=0você obtém um limite inferior negativo (de fato -0,290617).

> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3)
[1] 0.02464434
attr(,"control")
   dist alpha beta gamma delta pm
 stable   0.4    1   0.4     0  3
> dstable(4, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3)
[1] 0.02700602
attr(,"control")
   dist alpha beta gamma    delta pm
 stable   0.4    1   0.4 0.290617  3
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0, pm=3))
[1] -0.2876658
> min(rstable(100000, alpha=0.4, beta=1, gamma=0.4, delta=0.290617, pm=3))
[1] 0.004303485

Você pode achar mais fácil simplesmente ignorar deltaa saída e, contanto que continue beta=1usando os pm=1meios deltana entrada, é o limite inferior da distribuição, que parece que você deseja 0.

Também digno de nota: Martin Maechler apenas refatorou o código para o stable distribuído e adicionou algumas melhorias.

Seu novo pacote stabledist também será usado pelo fBasics, então você pode dar uma olhada também.