# 静态T-Copula
kappa1 <- cor(qnorm(u_afe_t), qnorm(u_gp_t))
tcopula_CL <- function(theta, Zdata) {
  # 解析参数
  kappa <- theta[1]
  nu <- theta[2]
  # 计算 x 和 y
  if (nu >= 100) {
    x <- qnorm(Zdata[, 1])
    y <- qnorm(Zdata[, 2])
  } else {
    x <- qt(Zdata[, 1], nu)
    y <- qt(Zdata[, 2], nu)
  }
  # 计算 CL
  CL <- lgamma((nu+2)/2) + lgamma(nu/2) - 2*lgamma((nu+1)/2) - 0.5*log(1-kappa^2)
  CL <- CL - (nu+2)/2*log(1+(x^2 + y^2 - 2*kappa*x*y)/(nu*(1-kappa^2)))
  CL <- CL + (nu+1)/2*log(1+x^2/nu) + (nu+1)/2*log(1+y^2/nu)
  CL <- sum(CL)
  CL <- -CL
  return(CL)
}
lower_tcopula <- c(-0.9, 2.1)  # 约束下限，包括nu的范围
upper_tcopula  <- c(0.9, 100)   # 约束上限，包括nu的范围
theta0_tcopula  <- c(kappa1, 20)  # 初始参数 [theta1, theta2, nu]
result_tcopula <- optim(theta0_tcopula, tcopula_CL, Zdata = urets, method = "L-BFGS-B", lower = lower_tcopula, upper = upper_tcopula)
result_tcopula$par[1]

# 动态Normal-Copula
tvbinormal_CL <- function(theta, Zdata, rhobar) {
  # 解析参数
  c <- theta[1]
  a <- theta[2]
  b <- theta[3]
  T <- nrow(Zdata)
  x <- qnorm(Zdata[, 1])
  y <- qnorm(Zdata[, 2])
  # 初始化 kappa
  kappa <- rep(-999.99, T)
  kappa[1] <- rhobar
  # 计算 kappa
  for (jj in 2:T) {
    if (jj <= 10) {
      psi <- c + a * mean(x[1:(jj-1)] * y[1:(jj-1)]) + b * kappa[jj - 1]
    } else {
      psi <- c + a * mean(x[(jj - 10):(jj - 1)] * y[(jj - 10):(jj - 1)]) + b * kappa[jj - 1]
    }
    kappa[jj] <- 1.998 / (1 + exp(-psi)) - 0.999 #归一化处理
  }
  rhohat <- kappa
  # 计算 CL
  CL <- -1 * (2 * (1 - kappa^2))^(-1) * (x^2 + y^2 - 2 * kappa * x * y)
  CL <- CL + 0.5 * (x^2 + y^2)  
  CL <- CL - 0.5 * log(1 - kappa^2)
  CL <- sum(CL)
  CL <- -CL
  return(list(CL = CL, rhohat = rhohat))
}
lower_tvbinormal <- c(-5, -5, -5)  # 约束下限
upper_tvbinormal <- c(5, 5, 5)   # 约束上限
theta0_tvbinormal <- c(0, 0, 0)  # 初始参数 [theta1, theta2, theta3]
opt_fun_tvbinormal <- function(theta, Zdata, rhobar) {
  tvbinormal_CL(theta[1:3], Zdata, rhobar)$CL
}
result_tvbinormal <- optim(theta0_tvbinormal, opt_fun_tvbinormal, Zdata = urets, rhobar = result_tcopula$par[1], method = "L-BFGS-B", lower = lower_tvbinormal, upper = upper_tvbinormal)
tvbinormal_rho <- tvbinormal_CL(result_tvbinormal$par, Zdata = urets, result_tcopula$par[1])$rhohat

# 动态T-Copula
tvt_CL <- function(theta, Zdata, rhobar) {
  # 解析参数
  c <- theta[1]
  a <- theta[2]
  b <- theta[3]
  nu <- theta[4]
  T <- nrow(Zdata)
  x <- qt(Zdata[, 1],nu)
  y <- qt(Zdata[, 2],nu)
  # 初始化 kappa
  kappa <- rep(-999.99, T)
  kappa[1] <- rhobar
  # 计算 kappa
  for (jj in 2:T) {
    if (jj <= 10) {
      psi <- c + a * mean(x[1:(jj-1)] * y[1:(jj-1)]) + b * kappa[jj - 1]
    } else {
      psi <- c + a * mean(x[(jj - 10):(jj - 1)] * y[(jj - 10):(jj - 1)]) + b * kappa[jj - 1]
    }
    kappa[jj] <- 1.998 / (1 + exp(-psi)) - 0.999 #归一化处理
  }
  rhohat <- kappa
  # 计算 CL
  CL <- lgamma((nu+2)/2) + lgamma(nu/2) - 2*lgamma((nu+1)/2) - 0.5*log(1-kappa^2)
  CL <- CL - (nu+2)/2*log(1+(x^2 + y^2 - 2*kappa*x*y)/(nu*(1-kappa^2)))
  CL <- CL + (nu+1)/2*log(1+x^2/nu) + (nu+1)/2*log(1+y^2/nu)
  CL <- sum(CL)
  CL <- -CL
  return(list(CL = CL, rhohat = rhohat))
}
lower_tvt <- c(-5, -5, -5, 2.1)  # 约束下限
upper_tvt <- c(5, 5, 5, 100)   # 约束上限
theta0_tvt <- c(0, 0, 0, 0)  # 初始参数 [theta1, theta2, theta3, nu]
opt_fun_tvt <- function(theta, Zdata, rhobar) {
  tvt_CL(theta[1:4], Zdata, rhobar)$CL
}
result_tvt <- optim(theta0_tvt, opt_fun_tvt, Zdata = urets, rhobar = result_tcopula$par[1], method = "L-BFGS-B", lower = lower_tvt, upper = upper_tvt)
tvt_rho <- tvt_CL(result_tvt$par, Zdata = urets, result_tcopula$par[1])$rhohat

为什么我动态T-Copula和动态Normal-Copula跑出来的结果完全相反？理论上动态Normal-Copula更正确？所以我的动态T-Copula的CL写错了？还是qt()返回的问题？