R中稳定分布参数的对数似然函数拟合

zd287kbt  于 2023-04-09  发布在  其他
关注(0)|答案(1)|浏览(132)

我试图用最大似然法来拟合稳定分布到我的数据中。我知道有类似“libstableR”的软件包,但我想自己估计参数,我不知道FFT如何帮助我获得对数似然函数。
例如,对于t-distirubtion,我估计的参数如下:

nllDax1 <- function(pars, xsD1){
  v <- pars[1]
  b <- pars[2]
  a <- pars[3] 
  mllk <- sum(log((1/a)*((gamma((v+1)/2)/(gamma(v/2)*gamma(0.5)*sqrt(v)))*
               (1+((((xsD1-b)/a)**2)/v))**(-(v+1)/2))))
  return(-mllk)}
mleD1 <- optim(c(v = 2,b = 0,a = 0.5), fn = nllDax1, method = "BFGS", 
               x = dx$DAX1$Tagesrendite[2:441])

现在我想用同样的方法来估计稳定分布的参数,但我不知道如何得到对数似然函数。
对于特征函数,我想使用这种方法:
特性功能:


有谁知道我该怎么编程吗。

eimct9ow

eimct9ow1#

您可以考虑以下方法:

library(DEoptim)
library(elliptic)

integrand_Stable <- function(x, u, m, c, alpha, beta)
{
  i <- complex(real = 0.0, imaginary = 1.0)
  term1 <- i * m * u

  if(alpha == 1)
  {
    term2 <- c * abs(u) * (1 + i * beta * sign(u) * 2 / pi * log(abs(u)))

  }else
  {
    term2 <- c ^ alpha * abs(u) ^ alpha * (1 - i * beta * sign(u) * tan(pi * alpha / 2))
  }

  term3 <- exp(term1 - term2) * exp(- i * u * x) / (2 * pi)

  return(term3)
}

my_log_Lik <- function(param, val)
{
  m <- param[1]
  c_Par <- param[2]
  alpha <- param[3]
  beta <- param[4]

  if((alpha > 2) | (alpha < 0) | (beta < -1) | (beta > 1) | (c_Par < 0))
  {
    return(10 ^ 30)

  }else if(((alpha < 1) & (beta == -1) & any(val > m)) | ((alpha < 1) & (beta == 1) & any(val < m)))
  {
    return(10 ^ 30)

  }else
  {
    nb_Data <- length(X)
    log_Lik_Val <- 0

    for(i in 1 : nb_Data)
    {
      temp_Val <- tryCatch(Re(myintegrate(f = integrand_Stable, lower = -1000, upper = 1000, x = val[i], m = m, c = c_Par, alpha = alpha, beta = beta)),
                           error = function(e) NA)
      log_Lik_Val <- log_Lik_Val + log(temp_Val)
    }

    if(is.na(log_Lik_Val) | is.nan(log_Lik_Val) | is.infinite(log_Lik_Val))
    {
      return(10 ^ 30)

    }else
    {
      return(-log_Lik_Val)
    }
  }
}

x_Val <- rnorm(100) + rexp(100)

obj_DEoptim <- DEoptim(fn = my_log_Lik, lower = c(-5, -5, 0, -1), upper = c(5, 5, 2, 1),
                       control = list(itermax = 1000),  val = x_Val)

obj_DEoptim$optim$bestmem

相关问题