R语言中Gumbel Copula极大似然估计失败问题求助

问题背景 #

使用R的copula包对GM、AAPL、MSFT三只资产的对数收益率进行极大似然估计（MLE）拟合Gumbel Copula时，二维资产组合拟合正常，但三维组合拟合出现Hessian矩阵奇异的警告。已尝试调整初始参数start、确认样本量（超3000）、通过VineCopula验证Gumbel Copula适配性，且资产间线性相关系数为0.38、0.40、0.60，问题仍未解决。

原代码 #

library(quantmod)
library(copula)

# 获取数据
getSymbols(c("GM", "AAPL","MSFT"), src = "yahoo",  from = "2009-01-01")
mydata <- data.frame(na.omit(diff(log(merge(Ad(GM),Ad(AAPL),Ad(MSFT))))))
cor(mydata)
# 样本量
(n <- nrow(mydata))
   
# 转换为均匀边缘分布
u <- sapply(mydata, function(x){
  ecx <- ecdf(x)
  ecx(x)
})
plot(data.frame(u))

# 二维Gumbel Copula拟合（正常运行）
fit.ml2 <- fitCopula(gumbelCopula(dim=2), u[,1:2], method="ml", start = 3)
summary(fit.ml2)

fit.ml2a <- fitCopula(gumbelCopula(dim=2), u[,c(2:3)], method="ml", start = 3)
summary(fit.ml2a)

# 三维Gumbel Copula拟合（触发警告）
fit.ml3 <- fitCopula(gumbelCopula(dim=3), u, method="ml", start = 3)
summary(fit.ml3)

报错信息 #

Warning message:
In fitCopula.ml(copula, u = data, method = method, start = start, :
Hessian matrix not invertible: Lapack routine dgesv: system is exactly singular: U[1,1] = 0

解决方案 #

1. 修正边缘分布转换的边界值 #

经验分布函数（ECDF）转换得到的均匀值可能出现0或1，这会导致Gumbel Copula的密度函数计算出现数值奇点，进而引发Hessian矩阵奇异。对边界值做微小调整：

u <- sapply(mydata, function(x){
  ecx <- ecdf(x)
  p <- ecx(x)
  # 替换0和1为极小/极大非边界值
  p[p == 0] <- 1e-6
  p[p == 1] <- 1 - 1e-6
  p
})

2. 基于二维拟合结果优化初始参数 #

三维Gumbel Copula为单参数Archimedean Copula，可先拟合两两二维Gumbel模型，取参数均值作为三维拟合的初始值：

# 拟合两两二维Gumbel Copula
fit_gm_aapl <- fitCopula(gumbelCopula(dim=2), u[,1:2], method="ml")
fit_aapl_msft <- fitCopula(gumbelCopula(dim=2), u[,2:3], method="ml")
fit_gm_msft <- fitCopula(gumbelCopula(dim=2), u[,c(1,3)], method="ml")

# 计算初始参数
start_theta <- mean(c(coef(fit_gm_aapl), coef(fit_aapl_msft), coef(fit_gm_msft)))

# 重新拟合三维模型
fit.ml3 <- fitCopula(gumbelCopula(dim=3, param=start_theta), u, method="ml")
summary(fit.ml3)

3. 基于Kendall tau转换初始参数 #

Gumbel Copula的参数theta与Kendall tau满足tau = 1 - 1/theta，可先计算样本Kendall tau均值，转换为初始参数：

# 计算样本Kendall tau矩阵
tau_mat <- cor(mydata, method="kendall")
# 转换为Gumbel theta初始值
theta_init <- 1 / (1 - mean(tau_mat[upper.tri(tau_mat)]))

# 拟合三维模型
fit.ml3 <- fitCopula(gumbelCopula(dim=3, param=theta_init), u, method="ml")
summary(fit.ml3)

4. 指定带约束的优化方法 #

Gumbel Copula的theta参数下界为1，使用带边界约束的优化方法（如L-BFGS-B）可提升数值稳定性：

fit.ml3 <- fitCopula(gumbelCopula(dim=3, param=start_theta), u, method="ml",
                     optim.control = list(method = "L-BFGS-B", lower = 1, upper = 10))
summary(fit.ml3)

内容的提问来源于stack exchange，提问作者fython