# Kentaro Fukumoto # Note (or Warning): This is very preliminary (due to time limitation). # I've not checked this script worlks well. # It just conveys the sense the models are implemented. rm(list=ls()) ################### # Univariate Models # truncated normal # Eq. 9 log.likelihood.truncated.normal <- function(parameter,y,x,cutoff){ standard.deviantion <- exp(parameter[1]) coefficient <- parameter[-1] linear.predictor <- x%*%coefficient sum( dnorm(y,linear.predictor,standard.deviantion,log=T) - pnorm(cutoff,linear.predictor,standard.deviantion,lower.tail=F,log.p=T) ) } # censored normal (or Tobit) # Eq. 19 or 24 log.likelihood.censored.normal <- function(parameter,y,x,cutoff,censor){ standard.deviantion <- exp(parameter[1]) coefficient <- parameter[-1] linear.predictor <- x%*%coefficient sum( (censor*dnorm(y,linear.predictor,standard.deviantion,log=T)) + ((1-censor)*pnorm(cutoff,linear.predictor,standard.deviantion,lower.tail=T,log.p=T)) ) } # survival models as a special case ################## # Bivariate Models # SUR (p. 8) library(mvtnorm) log.likelihood.SUR <- function(parameter,y1,y2,x1,x2){ k1 <- ncol(x1); k2 <- ncol(x2) correlation <- (1-exp(-parameter[1]))/(1+exp(-parameter[1])) standard.deviantion1 <- exp(parameter[(1+1)]) standard.deviantion2 <- exp(parameter[(1+1+k1+1)]) coefficient1 <- parameter[(1+1+1):(1+1+k1)] coefficient2 <- parameter[(1+1+1+k1+1):(1+1+k2+k1+1)] linear.predictor1 <- x1%*%coefficient1 linear.predictor2 <- x1%*%coefficient2 standardozed.y1 <- (y1-linear.predictor1)/standard.deviantion1 standardozed.y2 <- (y2-linear.predictor2)/standard.deviantion2 y <- cbind(standardozed.y1,standardozed.y2) sum(dmvnorm(y,mean=c(0,0),sigma=matrix(c(1,correlation,correlation,1),2,2),log=T)) } # Heckman selection model (p. 6) Heckman.two.step.estimator <- function(parameter,y1,y2,x1,x2){ selection.out <- glm(y1 ~ x1, family = binomial(link = "probit")) linear.predictor1.hat <- x1%*%selection.out$coef lambda.hat <- dnorm(linear.predictor1.hat)/pnorm(linear.predictor1.hat) lm(y2 ~ x2 + lambda.hat) } # ML estimtor # bivariate probit (Eq. 28) log.likelihood.bivariate.probit <- function(parameter,y1,y2,x1,x2){ number.of.observartions <- length(y1) k1 <- ncol(x1); k2 <- ncol(x2) correlation <- (1-exp(-parameter[1]))/(1+exp(-parameter[1])) coefficient1 <- parameter[(1+1):(1+k1)] coefficient2 <- parameter[(1+1+k1):(1+k2+k1)] linear.predictor1 <- x1%*%coefficient1 linear.predictor2 <- x1%*%coefficient2 prob1 <- prob2 <- prob12 <- array(dim=number.of.observartions) for(i in 1:number.of.observartions){ prob12[i] <-pmvnorm(lower=c(-Inf,-Inf), upper=c(linear.predictor1[i],linear.predictor2[i]), mean=c(0,0), sigma=matrix(c(1,correlation,correlation,1),2,2)) } prob1 <- pnorm(linear.predictor1) prob2 <- pnorm(linear.predictor2) sum((y1*y2*log(prob12)) + (y1*(1-y2)*log(prob1 - prob12)) + ((1-y1)*y2*log(prob2 - prob12)) + ((1-y1)*(1-y2)*log(1 - prob1 - prob2 - prob12)) ) } ####################### # Zero Inflated Poisson log.likelihood.ZIP <- function(parameter,y,x){ no.risk.probablity <- 1/(1+exp(-parameter[1])) coefficient <- parameter[-1] expectation <- exp(x%*%coefficient) sum(log( (no.risk.probablity*as.integer(y==0)) + ((1-no.risk.probablity)*dpois(y,expectation)) )) }