#=============================================================== # Econ 21410 -- A final example putting together pieces -- May 28, 2014 # From code by Alex MacKay with changes made by John Eric Humphries # A Single-Factor Random Coefficient Model very similar to Aakvik, Heckman, and Vytlacil (2005) #=============================================================== rm(list=ls()) library(mvtnorm) library(statmod) library(doMC) library(RcppArmadillo) #library(doParallel) # this package and the next 3 lines if on a windows machine #cl <- makeCluster(8) #registerDoParallel(cl) #=================================== # Defining some functions #=================================== #------------------------------------------------------ # Function: Return likelihood with known theta #------------------------------------------------------ f.value1 = function(par=par,data=data){ # calculates the log likelihood assuming theta is known. beta1 = par[1] beta0 = par[2] gamma = par[3:4] alpha1 = par[5] alpha0 = par[6] likelihood = pnorm(colSums(c(gamma,1)*t(data[,3:5])))^d* dnorm(y - colSums(c(beta1,alpha1)*t(data[,c(3,5)])))^d* (1-pnorm(colSums(c(gamma,1)*t(data[,c(3:5)]))))^(1-d)* dnorm(y - colSums(c(beta0,alpha0)*t(data[,c(3,5)])))^(1-d) ll = sum(log(likelihood)) return(-ll) } #------------------------------------------------ # Function: Return likelihood with unknown theta #------------------------------------------------ f.value2 = function(par=par,data=data){ # Calculates the log likelihood omitting theta beta1 = par[1] beta0 = par[2] gamma = par[3:4] likelihood = pnorm(1/sqrt(2)*colSums(c(gamma)*t(data[,3:4])))^d* dnorm(1/sqrt(5)*(y-colSums(c(beta1)*t(data[,c(3)]))))^d* (1-pnorm(1/sqrt(2)*colSums(c(gamma)*t(data[,3:4]))))^(1-d)* dnorm(1/sqrt(2)*(y-colSums(c(beta0)*t(data[,c(3)]))))^(1-d) ll = sum(log(likelihood)) return(-ll) } f.value2.alt = function(par=par,data=data){ # Calculates the log likelihood omitting theta beta1 = par[1] beta0 = par[2] gamma = par[3:4] s = par[5] s1 = par[6] s2 = par[7] likelihood = pnorm(colSums(c(gamma)*t(data[,3:4])),sd=s)^d* dnorm((y-colSums(c(beta1)*t(data[,c(3)]))),sd=s1)^d* (1-pnorm(colSums(c(gamma)*t(data[,3:4])),sd=s))^(1-d)* dnorm((y-colSums(c(beta0)*t(data[,c(3)]))),sd=s2)^(1-d) ll = sum(log(likelihood)) return(-ll) } #---------------------------------------------------------------- # Function: Return likelihood with theta as a random coefficient #---------------------------------------------------------------- f.value3 = function(par=par,data=data){ # Calculates the log likelihood integrating out theta data_orig = data n = dim(data)[1] beta1 = par[1] beta0 = par[2] gamma = par[3:4] alpha1 = par[5] alpha0 = par[6] q_num = 16 quad = gauss.quad.prob(q_num,dist="normal",mu=0,sigma=1) values = matrix(0,q_num,n) for(i in 1:q_num){ theta = rep(quad$nodes[i],n) data = cbind(data_orig,theta) likelihood = pnorm(colSums(c(gamma,1)*t(data[,3:5])))^d* dnorm(y - colSums(c(beta1,alpha1)*t(data[,c(3,5)])))^d* (1-pnorm(colSums(c(gamma,1)*t(data[,c(3:5)]))))^(1-d)* dnorm(y - colSums(c(beta0,alpha0)*t(data[,c(3,5)])))^(1-d) values[i,] = likelihood } combined=colSums(quad$weights*values) return(-sum(log(combined))) } #================================================ # Some Faster Functions: #================================================ f.value3.cpp = function(par=par,data_orig=data){ # Calculates the log likelihood integrating out theta q_num = 16 quad = gauss.quad.prob(q_num,dist="normal",mu=0,sigma=1) likloop(quad$nodes,quad$weights,data_orig,par) } # FOR WINDOWS: # library(doParallel) # cl <- makeCluster(2) # registerDoParallel(cl) registerDoMC(8) foreach(i=1:3) %dopar% sqrt(i) f.value3.par = function(par=par,data=data){ # let it know to use 8 cores # Calculates the log likelihood integrating out theta d <- data[,1] y <- data[,2] data_orig = data n = dim(data)[1] beta1 = par[1] beta0 = par[2] gamma = par[3:4] alpha1 = par[5] alpha0 = par[6] q_num = 16 quad = gauss.quad.prob(q_num,dist="normal",mu=0,sigma=1) combined <- foreach(i = 1:q_num, .combine = '+', .inorder = F) %dopar% { theta = rep(quad$nodes[i],n) data = cbind(data_orig,theta) library(statmod) likelihood = pnorm(colSums(c(gamma,1)*t(data[,3:5])))^d* dnorm(y - colSums(c(beta1,alpha1)*t(data[,c(3,5)])))^d* (1-pnorm(colSums(c(gamma,1)*t(data[,c(3:5)]))))^(1-d)* dnorm(y - colSums(c(beta0,alpha0)*t(data[,c(3,5)])))^(1-d) (likelihood * quad$weights[i]) } return(-sum(log(combined))) } #===== Now a similar function using a lapply statement instead of a foreach loop. Lik = function(i) { require(statmod) theta = rep(quad$nodes[i],n) data = cbind(data_orig,theta) likelihood = pnorm(colSums(c(gamma,1)*t(data[,3:5])))^d* dnorm(y - colSums(c(beta1,alpha1)*t(data[,c(3,5)])))^d* (1-pnorm(colSums(c(gamma,1)*t(data[,c(3:5)]))))^(1-d)* dnorm(y - colSums(c(beta0,alpha0)*t(data[,c(3,5)])))^(1-d) return(likelihood * quad$weights[i]) } f.value3.par2 = function(par=par,data=data){ require(parallel) d <- data[,1] y <- data[,2] data_orig = data n = dim(data)[1] beta1 = par[1] beta0 = par[2] gamma = par[3:4] alpha1 = par[5] alpha0 = par[6] q_num = 16 quad = gauss.quad.prob(q_num,dist="normal",mu=0,sigma=1) combined = simplify2array(mclapply(1:q_num,Lik,mc.cores=4)) Like=rowSums(combined) return(-sum(log(Like))) } #==================================================== # Simulating data #===================================================== # Variables n = 1000 x = runif(n)*5 z = runif(n)*3 theta = rnorm(n) # Parameters beta1 = 1.5 beta0 = 1.4 alpha1 = 1.2 alpha0 = -2 gamma = c(.5,-1) true_par = c(beta1,beta0,gamma,alpha1,alpha0) mu = c(0,0,0) sigma = diag(c(1,1,1)) # Simulating outcomes. errors = rmvnorm(n,mu,sigma) matrix = cbind(x,z,theta,errors) colnames(matrix) = c("x","z","theta","e1","e0","eD") d_star = colSums(c(gamma,1,1)*t(matrix[,c(1:3,6)])) y1_star = colSums(c(beta1,alpha1,1)*t(matrix[,c(1,3,4)])) y0_star = colSums(c(beta0,alpha0,1)*t(matrix[,c(1,3,5)])) d = d_star > 0 y1 = y1_star y0 = y0_star y = d*y1+(1-d)*y0 #=============================================== # Solving the model under various assumptions. #=============================================== #---------------------------- # Data observed by researcher #---------------------------- matrix2 = cbind(d,y,y1,y0,d_star,matrix) write.csv(matrix2,file="dataset_a.csv",row.names=FALSE) #---------------------------- # Suppose we observe everything #----------------------------- data1 = matrix2[,c(1,2,6,7,8)] param1 = true_par result.all = optim(par=param1, fn=f.value1,data=data1) true_par result.all$par #---------------------------- # What happens if we just omit the factor? #--------------------------- data2 = matrix2[,c(1,2,6,7)] param2 = true_par[1:4] result.omit = optim(par=param2, fn=f.value2,data=data2) true_par result.omit$par result.omit = optim(par=c(param2,1,1,1), fn=f.value2.alt,data=data2) #------------------------------------------- # What happens if we integrate out the factor? #------------------------------------------- #result.int = optim(par=param1, fn=f.value3,data=data2) #true_par #result.all$par #result.omit$par #result.int$par # show timing for various sample sized # show how bias increases when alphas are larger. #------------------------------------------ # Timing our new functions #------------------------------------------ system.time(result.int <- optim(par=param1, fn=f.value3,data=data2,method="BFGS")) system.time(result.int.par <- optim(par=param1, fn=f.value3.par,data=data2,method="BFGS")) system.time(result.int.par <- optim(par=param1, fn=f.value3.par2,data=data2,method="BFGS")) system.time(result.int.cpp <- optim(par=param1, fn=f.value3.cpp,data=data2,method="BFGS")) #stopCluster(cl)