r - 如何在 r 中实现分数驱动的 Copula 模型

问题描述

我想实现一个基于Creal et all (2013)的 GAS 框架的动态 copula 。我从一个双变量学生的 $t$ copula 开始，并考虑到 copula 参数的更新方程是：

$\hat\rho_{t} = \delta + \gamma\hat\rho_{t-1} + \kappa\upsilon_{\rho_{t-1}}$

copula 的对数密度分数在哪里
$\upsilon_{\rho} = \frac{\partial{\ln{c}}}{\partial{\rho}}$
，我为 copula 得出这个分数 $t$ （如Ayala 和 Blazsek https://www.tandfonline.com/doi/full/10.1080/1351847X.2018.1464488?src=recsys，这里有也是copula密度函数）：

$\upsilon_{\rho} = \frac{\rho}{1-\rho^{2}}+\frac{\nu+2}{\rho^{2}-1}\cdot\frac{\rho[ T_{\nu}^{-1}(u)]^{2} +[T_{\nu}^{-1}(v)]^{2} - (\rho^{2}+1)T_ {\nu}^{-1}(u)T_{\nu}^{-1}(v)}{[T_{\nu}^{-1}(u)]^{2}+[T_{ \nu}^{-1}(v)]^{2}- 2\rho T_{\nu}^{-1}(u)T_{\nu}^{-1}(v)-\nu( \rho^{2}-1)}$

我的问题是，当我尝试优化最大似然参数估计时，不会以任何方式收敛。不知道是不是参数约束相关的问题，还是我在函数中插入的循环判断错误 $\帽子\rho$ 。这是代码

loglik_t = function(param) {
  dfactor = numeric()
  rho = numeric()
  score = numeric()
  tau =  cor(u_sp, u_eu, method = "kendall")
  dfactor[1] = tau
  tu = qt(u_sp, nu)  
  tv = qt(u_eu, nu)   
  omega = param[1]
  dA = exp(param[2])
  dB = exp(param[3])/(1 + exp(param[3]))
  x2 = log(gamma((nu + 2)/2)) +  log(gamma(nu/2)) - log(digamma((nu + 1)/2))
  LL = 0
  for (i in 1:length(u_sp)){
    rho[i] = sin((pi/2)*dfactor[i])
    #log-likelihood for observation i
    x1 = -0.5*log(1 - rho[i]^2)
    x3 =((-nu + 2)/2)*log(1 + (tu[i]^2+tv[i]^2-2*rho[i]*tu[i]*tv[i])/(nu*(1-rho[i]^2)))
    x4 = ((-nu + 1)/2)*log(1 + tu[i]^2/nu) + ((-nu + 1)/2)*log(1 + tv[i]^2/nu)
    c_t = x1 + x2 + x3 - x4
    LL = LL + c_t

    #GAS dynamic
    x1score = rho[i]/(1 - rho[i]^2)
    x2score = (nu + 2)/(rho[i]^2 - 1)
    x3score = rho[i]*(tu[i]^2 + tv[i]^2) - (rho[i]^2 + 1)*tu[i]*tv[i]
    x4score = tu[i]^2 + tv[i]^2 - 2*rho[i]*tu[i]*tv[i] - nu*(rho[i]^2 - 1)
    score[i] = x1score + x2score*(x3score/x4score)

    dfactor[i + 1] = omega + dA*score[i] + dB*rho[i]
  }

  return(list(LL, rho, score))
}

n = length(u_sp)

qlogL.hat = -10^10
for (i in (1:10)) {

    # Initial values

    start.theta = runif(3, min = 0.01, max = 0.9999)



    qml.fit = optim(start.theta, fn = function(x) loglik_t(x)[[1]], 
                    hessian = TRUE, method = "L-BFGS-B", 
                    lower = c(-2,0.0001,0.0001),
                    upper = c(0.9999, 0.999, 0.9999),
                    control = list(trace = 1, fnscale = -1))




    if (qml.fit$value > qlogL.hat) {
        theta.hat = qml.fit$par
        inv.J.hat = solve(-qml.fit$hessian/(n-1))
        jacobian.hat = pracma:: jacobian(loglik_t1,theta.hat)   
        I.hat = (t(jacobian.hat) %*% jacobian.hat)/(n-1)
        sd.hat = sqrt(diag(inv.J.hat %*% I.hat %*% inv.J.hat)/(n-1))
        qlogL.hat = qml.fit$value
        i
    }

}

标签： roptimization

r - 如何在 r 中实现分数驱动的 Copula 模型

问题描述

解决方案

推荐阅读