Максимизация функции правдоподобия, содержащей pbivnorm

Question

У меня есть функция правдоподобия, содержащая двумерный нормальный CDF. Я продолжаю получать значения, близкие к одному для корреляции, даже если истинное значение равно нулю.

Выбор образца R-пакета максимизирует функцию правдоподобия, содержащую нормальный CDF бивараита (как в Van de Ven и Van Praag (1981)). Я попытался посмотреть исходный код пакета, но не смог найти, как они пишут вероятность. Для справки, Ван де Вен и бумаги Ван Praag в:

функция

The Demand for Deductibles in Private Health Insurance: A Probit Model with Sample Selection.

вероятность того, уравнение (19), где Н обозначает стандартное нормальное CDF и H_2 обозначает двумерное нормальное КОР.

Мой вопрос:

1. Может кто-нибудь сказать мне, как функция правдоподобия записывается в пакете sampleSelection? или

2. Может кто-нибудь сказать мне, почему я получаю значения близких к одному для корреляции в коде ниже?

Вот код, который держит меня на ночь:

######################################################## 
# 
# Trying to code Van de Ven and Van Praag (1981) 
# 
# 
######################################################## 
library(MASS) 
library(pbivnorm) 
library(mnormt) 
library(maxLik) 
library(sampleSelection) 
set.seed(1) 

# Sample size 
full_sample <- 1000 

# Parameters 
rho  <- .1 
beta0 <- 0 
beta1 <- 1 
gamma0 <- .2 
gamma1 <- .5 
gamma2 <- .5 
varcovar <- matrix(c(1,rho,rho,1), nrow = 2, ncol = 2) 
# Vectors for storing values 
    y <- rep(0,full_sample) 
    s <- rep(0,full_sample) 
outcome <- rep(0,full_sample) 
select <- rep(0,full_sample) 

####################### 
# Simulate data 
####################### 
x <- rnorm(full_sample) 
z <- rnorm(full_sample) 

errors <- mvrnorm(full_sample, rep(0,2), varcovar) 
# note: 1st element for selection eq; 2nd outcome 
s <- gamma0 + gamma1*x + gamma2*z + errors[,1] 
y <- beta0 + beta1*x + errors[,2] 

for(i in 1:full_sample){ 
    if(s[i]>=0){ 
     select[i] <- 1 
     if(y[i]>=0){ 
     outcome[i] <- 1 
     }else{ 
     outcome[i] <- 0 
     } 
    }else{ 
     outcome[i] <- NA 
     select[i] <- 0 
    } 

} 
####################################### 
# Writing the log likelihood 
## 
# Note: vega1= beta0, 
#  vega2= beta1, 
#  vega3= gamma0, 
#  vega4= gamma1, 
#  vega5= gamma2, 
#  vega6= rho 
####################################### 
first.lf <- function(vega) { 

# Transforming this parameter becuase 
# correlation is bounded between -1 aad 1 
    corr <- tanh(vega[6]) 

# Set up vectors for writing log likelihood 
    y0 <- 1-outcome 
    for(i in 1:full_sample) {  
    if(is.na(y0[i])){ y0[i]<- 0} 
    if(is.na(outcome[i])){ outcome[i]<- 0} 
    } 
    yt0 <- t(y0) 
    yt1 <- t(outcome) 
    missing <- 1 - select 
    ytmiss <- t(missing) 

# Terms in the selection and outcome equations 
    A <- vega[3]+vega[4]*x+vega[5]*z 
    B <- vega[1]+vega[2]*x 
    term1 <- pbivnorm(A,B,corr) 
    term0 <- pbivnorm(A,-B,corr) 
    term_miss <- pnorm(-A) 
    log_term1 <- log(term1) 
    log_term0 <- log(term0) 
    log_term_miss <- log(term_miss) 
# The log likelihood 
    logl <- sum(yt1%*%log_term1 + yt0%*%log_term0 + ytmiss%*%log_term_miss) 
return(logl) 
} 

startv <- c(beta0,beta1,gamma0,gamma1,gamma2,rho) 
# Maxmimizing my likelihood gives 
maxLik(logLik = first.lf, start = startv, method="BFGS") 
# tanh(7.28604) = 0.9999991, far from the true value of .1 

# Using sampleSelection package for comparison 
outcomeF<-factor(outcome) 
selectEq <- select ~ x + z 
outcomeEq <- outcomeF ~ x 
selection(selectEq, outcomeEq) 
# Notice the value of -0.2162 for rho compared to my 0.9999991 

Answer 1

Случается, что есть опечатка в статье в уравнении (19). Слагаемые от i = N_1 + 1 до N должны иметь -rho, а не rho. Таким образом, используя

term0 <- pbivnorm(A,-B,-corr) 

дает

maxLik(logLik = first.lf, start = startv, method="BFGS") 
# Maximum Likelihood estimation 
# BFGS maximization, 40 iterations 
# Return code 0: successful convergence 
# Log-Likelihood: -832.5119 (6 free parameter(s)) 
# Estimate(s): 0.3723783 0.9307454 0.1349979 0.4693686 0.4572421 -0.219618 

по мере необходимости.