を使用して対数尤度を推定していますoptim()
。固有値に問題があり、有効なヘッセ行列を見つけることができないため、標準誤差を計算できません。
私が受け取る「警告」メッセージは次のとおりです。
Warning messages:
1: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
NaNs produced
2: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
NaNs produced
3: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
NaNs produced
4: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
NaNs produced
5: In log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values) :
NaNs produced
誰かが問題を見つけるのを手伝ってくれますか?
これが私のコードです:
tobit.ll <- function(theta,x,y1,y2,y3)
{
library("copula")
library(mvtnorm)
p <- nrow(theta)
n<-nrow(x)
x<-as.matrix(cbind( 1,x))
km<- as.numeric(ncol(x))
beta1 <- theta[1:km]
beta2 <- theta[(km+1):(2*km)]
beta3 <- theta[(2*km+1):(km*3)]
s1 <- theta[(3*km+1)] ### If sigma^2 is given, take sqrt()
if (s1 < 0) return(1e-10)
s2 <- theta[(3*km+2)] ### If sigma^2 is given, take sqrt()
if (s2 < 0) return(1e-10)
s3 <- theta[(3*km+3)]
if (s3 < 0) return(1e-10)
rho12 <- theta[(3*km+4)]
if (rho12< -1 || rho12 > 1) return(NA)
rho13 <- theta[(3*km+5)]
if (rho13< -1 || rho13> 1) return(NA)
rho23 <- theta[(3*km+6)]
if (rho23< -1 || rho23> 1) return(NA)
fy1 <- ifelse(y1 >0, dnorm(y1,x%*%beta1, s1), (1-pnorm((x%*%beta1)/s1)))
fy2 <- ifelse(y2 >0, dnorm(y2,x%*%beta2, s2), (1-pnorm((x%*%beta2)/s2)))
fy3 <- ifelse(y3 >0, dnorm(y3,x%*%beta3, s3), (1-pnorm((x%*%beta3)/s3)))
Fy1 <- ifelse(y1>0, pnorm((y1-x%*%beta1)/s1),(1-pnorm((x%*%beta1)/s1)))
Fy2 <- ifelse(y2>0, pnorm((y2-x%*%beta2)/s2),(1-pnorm((x%*%beta2)/s2)))
Fy3 <- ifelse(y3>0, pnorm((y3-x%*%beta3)/s3),(1-pnorm((x%*%beta3)/s3)))
norm.cop= ellipCopula("normal", param = c(rho12, rho13, rho23),dim = 3, dispstr = "un")
u <- as.matrix(cbind(Fy1,Fy2,Fy3))
dcop<- dCopula(u, norm.cop)
fy1.lg <- ifelse(fy1==0, log(0.0001), log(fy1))
fy2.lg <- ifelse(fy2==0, log(0.0001), log(fy2))
fy3.lg <- ifelse(fy3==0, log(0.0001), log(fy3))
dcop.lg <- ifelse(dcop==0, log(0.0001), log(dcop))
lglik <- sum(fy1.lg+fy2.lg+fy3.lg) +sum(dcop.lg)
return(-(lglik))
}