你好,我想用msl方法估计一个面板probit,源代码如下,问题在于最终提示为“ initial value in 'vmmin' is not finite”,请问错误出在哪里,谢谢!
----------------------- copy starting from the next line -----------------------
* Example generated by -dataex-. To install: ssc install dataex
clear
## Code to estimate xtprobit model using MSL
##
library(haven)
library(summarytools)
library(matrixStats)
# --- prepare data ---
setwd("C:\\Users\\ThinkPad\\Desktop\\MSL for R")
xtprobit_data <- read_dta("xtprobit_data.dta")
descr(xtprobit_data)
descr(xtprobit_data, stats = c("n.valid", "mean", "med", "sd", "min", "max"), transpose = TRUE)
y <- unlist(xtprobit_data["union"])
X <- cbind(1, as.matrix(xtprobit_data[c("age", "grade", "not_smsa", "south")]))
nt <- 8
##
## 0. initial value comes from logit
##
ans_probit <- glm(union ~ age + grade + not_smsa + south,
family = gaussian(link = "identity"),
data = xtprobit_data)
initval <- c(coef(ans_probit), lnsig = 0)
## 1. xtlogit using Halton Sequence
# --- (1) set up halton sequence ---
library(randtoolbox)
nh <- 200
hgrid <- halton(nh * length(y) / nt, dim = 1, normal = TRUE)
htmat <- matrix(hgrid, length(y) / nt, nh, byrow = TRUE)
htmat <- htmat %x% rep(1, nt)
# --- (2) define likilihood function ---
xtprobit_llk_hal <- function(parm, y, X, nt, htmat) {
nx <- ncol(X)
nr <- ncol(htmat)
sig <- exp(parm[nx + 1])
Xb <- c(X %*% parm[1:nx])
mprh <- pnorm(Xb + sig * htmat)
mprh1<- log10(mprh)
mprh2<- log10(1-mprh)
mprh <- y * mprh1 + (1 - y) * mprh2
mpri <- rowProds(matrix(mprh, ncol = nt, byrow = TRUE)) #<--(N*R)*1
llki <- log(rowMeans(matrix(mpri, ncol = nr)))
return(-sum(llki))
}
# --- (3) to the opimization ---
ans_xtprobitt_hal <- optim(par = initval, fn = xtprobit_llk_hal, method = "BFGS",
control = list(trace = 1, maxit = 500), y = y,
X = X, nt = nt, htmat = htmat, hessian = TRUE)
------------------ copy up to and including the previous line ------------------