Calculando a função densidade de probabilidade da distribuição estável usando FFT.

Distribuições estáveis ​​paretianas têm propriedades atraentes para modelagem empírica em finanças, porque incluem a distribuição normal como um caso especial, mas também pode permitir caudas mais pesadas e assimetria.

Uma razão principal para a pouca utilização dessa distribuição em trabalhos acadêmicos aplicados é devido ao fato de que, em geral, não há expressão de "forma fechada" para a a função de densidade de probabilidade, e que as aproximações numéricas computacionais são não-triviais e computacionalmente extensivas.

Nesse post vou mostrar como é possível calcular a função densidade de probabilidade via Fast-Fourier Transform (FFT).

O trabalho original sobre esse assunto foi produzido por Mittnik, Doganoglu e Chenyao (1999).

A distribuição alfa-estável.

A distribuição alfa-estável, em geral, não possui expressão analítica para sua função densidade de probabilidade (f.d.p) ou ainda para a sua função distribuição acumulada (f.d.a), mas pode ser escrita por meio de sua função característica (Rachev e Mittnik, 2000 ):

onde $0<\alpha<2$ é o expoente da distribuição ou índice de cauda, $-1\le \beta\le 1$ é o parâmetro de assimetria, $\sigma>0$ é o parâmetro de escala e $\delta \in \Re$ é o parâmetro de locação, a função $\omega(.,.)$ é dada por:

A distribuição alfa-estável representada acima e definida pela notação $S_{\alpha}^{0}(\beta,\delta,\sigma)$ é denominada parametrização $S_{0}$ segundo Nolan (2010).

A função densidade de probabilidade pode ser aproximada utilizando o método FFT (Fast Fourier Transform) o qual é computacionalmente eficiente e permite um processo de aproximação mais rápido do que expansão por séries (Bergström, 1952) ou integração direta (Nolan, J. P., 2001. Maximum likelihood estimation of stable parameters. Manuscrito não publicado.).

Segundo Durrett (2010) página 106 uma função densidade de probabilidade pode ser escrita pela Transformada de Fourier da função característica, em outras palavras:

A integral acima pode ser calculada para $n$ pontos igualmente espaçados com distância $h$ e soma resultante pode ser computada por meio do método FFT (Fast Fourier Transform). Mittnik e Doganoglu (1999) sugerem que os valores de $n$ e $h$ devem ser respectivamente $2^{13}$ e $h=0.01$ para que uma boa aproximação seja possível.

Implementação no R.

Podemos então implementar o método descrito anteriormente no R.

Inicialmente, definimos os valores dos parâmetros para os quais desejamos calcular a função densidade de probabilidade:

#Limpa o console
rm(list=ls(all=TRUE))

#Parâmetros para o calculo da PDF
alpha<-0.5
beta<-0.3
sigma<-1
delta<-0

O limite superior para o suporte da distribuição $xmax$ e o número de pontos utilizados para a aproximação $n$ deve ser definido:

#Valor máximo em módulo da PDF
xmax<-20

#2^N espaços iguais entre [-xmax,xmax]
n<-13

Devido à natureza do método FFT, valores distantes do centro da distribuição podem ser subestimados. Por esta razão o método aqui sugerido calcula a função densidade de probabilidade da distribuição alfa-estável para o intervalo $[-xmax, xmax] \times 2^{mult}$ onde o multiplicador $mult$ é utilizado para aumentar (temporariamente) o suporte da distribuição, em seguida, esse intervalo aumentado é truncado em seu intervalo original. Seja o valor $mult=4$ como uma primeira tentativa, no entanto, para melhor precisão deve-se utilizar $mult>4$.

#Define o multiplicador
mult<-4

#Calcula o intervalo do suporte aumentado
xmax <- xmax*(2^mult)

#Aumenta n para cobrir o intervalo maior com os mesmos pontos do grid
n <- n + mult
M <- 2^n
R <- pi/xmax
dt <- 1/(R*M)

Como o método FFT admite a utilização de números complexos e devido a função característica da distribuição alfa-estável possuir número imaginário na sua estrutura precisamos definí-lo no R:

#Número imaginário
i<-complex(real = 0, imaginary = 1)

#Define o grid de valores com pontos uniformemente espaçados 
xx <-seq(from = -2^(n-1)+0.5, to = (2^(n-1)-0.5))/(2^n*dt)

Em seguida, a função característica é definida:

#Função característica da distribuição alfa-estável
piby2 <- pi/2
if(abs(alpha-1)<0.001)
  {
    yy <- exp( -sigma*(abs(xx))*( 1+i*beta*sign(xx)/piby2*log(sigma*abs(xx))) + i*delta*xx )
  }
else
{
    yy <- exp( -(sigma*abs(xx))^alpha*( 1+i*beta*sign(xx)*tan(alpha*piby2)*( (sigma*abs(xx))^(1-alpha)-1 ) ) + i*delta*xx )
}

Note que caso $\alpha$ esteja muito próximo de um, admitimos que $\alpha=1$ e utilizamos a estrutura da função característica nessa condição. O próximo passo é aplicar o método FFT:

#Utiliza o método FFT
yy1 <- c(yy[seq(from=(2^(n-1)+1),to=2^n)], yy[seq(from=1,to=2^(n-1))])
z <- Re( fft(yy1) )/(2*pi)*R

E após a utilização do método FFT pela função fft do R, a função densidade de probabilidade é calculada:

#Computa a função densidade de probabilidade
x <- (2*pi)*(seq(from=0,to=(M-1),by=1)/(M*R)-1/(2*R))
y <- c(z[seq(from=(2^(n-1)+1),to=2^n)], z[seq(from=1,to=2^(n-1))])

Por fim, o intervalo original do suporte da distribuição é obtido e o esboço de sua função densidade de probabilidade é apresentado:

#Encontra os pontos da f.d.p no intervalo original
T <- which((x<=xmax/(2^mult)) & (x>=-xmax/(2^mult)))
x <- x[T]
y <- y[T]

#Plota a função densidade de probabilidade
plot(x,y,type="l", main ="Distribuição alfa-estável"
+,xlab="Suporte",ylab="f.d.p",col="blue")

O valor da função densidade de probabilidade para um ponto arbitrário qualquer pode ser obtido por meio de interpolação entre os pontos que contém o valor desejado. Já a aproximação via FFT para outras parametrizações segue o mesmo processo, exceto pela definição da expressão da função característica na sintaxe apresentada anteriormente que deve ser alterada para a parametrização desejada.