####################################################################################### # Charles B. Moss - cbmoss@ufl.edu April 16, 2010 # ####################################################################################### # Data Section # ####################################################################################### # Read in the draw matrix from the estimated Leontief function # dta <- read.csv("Leontief-1.csv") dta[,5:6] <- 10000*dta[,5:6] # Preliminary Regression # y <- dta[,2] n <- dta[,3] p <- dta[,4] nn <- sqrt(dta[,3]) pp <- sqrt(dta[,4]) np <- nn * pp olsorig <- lm(y ~ n + p + nn + pp + np) # Create the matrix of zeros compatable with the data matrix # o <- matrix(1,nrow=nrow(dta),1) # Combine the x data into a single matrix # xdta <- as.matrix(cbind(o,dta[,3:6])) ####################################################################################### # Generate Missing Data # ####################################################################################### # Create b1 and b2 to draw missing values based on the logit # b1 <- as.matrix(c(-0.4500, 0.01125,-0.01775,-0.00749,-0.00624)) b2 <- as.matrix(c(-0.4875,-0.02275, 0.01675,-0.00852,-0.00710)) # Create xhat or the estimated probabilities based on the logit and the above b1 and # # b2 vectors # xhat <- cbind(exp(xdta %*% b1)/(1+exp(xdta %*% b1)), exp(xdta %*% b2)/(1+exp(xdta %*% b2))) #*************************************************************************************# # Loop to generate datasets with missing variables # #*************************************************************************************# # Draw i samples with missing variables # for (i in 1:10) { # ndta create a temporary dataset to hold the 'observed' data with the missing # # values # ndta <- dta[,2:6] # phat create a temporary dataset to contain the 'trigger' probabilities for # # setting the missing. Draw rbinom(1=number of randomes,1=number of draws,prob) # phat <- xhat for (j in 1:nrow(dta)) { for (k in 1:2) phat[j,k] <- rbinom(1,1,1-xhat[j,k]) for (k in 1:2) if (phat[j,k] == 0) ndta[j,k+1] <- NA } # Based on the random draw, setup the data for the Leontief production function # y <- ndta[,1] n <- ndta[,2] p <- ndta[,3] nn <- sqrt(ndta[,2]) pp <- sqrt(ndta[,3]) np <- nn * pp # Estimation the least square regression model dropping the # # missing variables # lsq <- lm(y ~ n + p + nn + pp + np) # Collect the results of least squares regression with dropped variables into # # lsqres # if (i == 1) lsqres <- lsq$coefficients else lsqres <- cbind(lsqres,lsq$coefficients) if (i == 1) lsqcov <- sqrt(diag(vcov(lsq))) else lsqcov <- cbind(lsqcov,sqrt(diag(vcov(lsq)))) # Construct the dataset without the observations with missing variables # nndta <- na.omit(ndta) o <- matrix(1,nrow=nrow(nndta)) # Construct the x1 matrix with the data that can be used to estimate the first # # missing variable (nitrogen) and y1 the observed levels of nitrogen # x1 <- as.matrix(cbind(o,nndta[,c(3,4,5)])) y1 <- as.matrix(nndta[,2]) lsq1 <- lm(y1~x1[,2] + x1[,3] + x1[,4]) # Estimate bx1 the estimated coefficients for predicting x1 # bx1 <- solve(t(x1)%*%x1)%*%(t(x1)%*%y1) # Estimate the mean squared error for the regresion for x1 # mse1 <- t(y1 - x1%*%bx1)%*%(y1 - x1%*%bx1)/(nrow(y1)-ncol(x1)) # Compute the covariance estimate for bx1 # vcb1 <- mse1[1,1]*solve(t(x1)%*%x1) dvcb1 <- diag(sqrt(diag(vcb1))) vvcb1 <- solve(dvcb1) %*% vcb1 %*% solve(dvcb1) chvvcb1 <- chol(vvcb1) # Construct the x2 matrix with the data that can be used to estimate the second # # missing variable (phosphorous) and y2 the observed levels of phosphorous # x2 <- as.matrix(cbind(o,nndta[,c(2,4,5)])) y2 <- as.matrix(nndta[,3]) lsq2 <- lm(y2~x2[,2]+x2[,3]+x2[,4]) # Estimate bx2 the estimated coefficients for the regression for x2 # bx2 <- solve(t(x2)%*%x2)%*%(t(x2)%*%y2) # Estimate the mean squared error for the regression for x2 # mse2 <- t(y2-x2%*%bx2)%*%(y2-x2%*%bx2)/(nrow(y2)-ncol(x2)) # Compute the covariance estimate of bx2 # vcb2 <- mse2[1,1]*solve(t(x2)%*%x2) dvcb2 <- diag(sqrt(diag(vcb2))) vvcb2 <- solve(dvcb2) %*% vcb2 %*% solve(dvcb2) chvvcb2 <- chol(vvcb2) #----------------------------------------------------------------------------------# # Draw from the estimated distributions for bx1 and bx2 (the regresion values on # # the missing values # #----------------------------------------------------------------------------------# jcount <- 0 for (j in 1:1000) { odraw <- matrix(1,nrow=nrow(ndta)) # Draw from the b1 distribution # b1hat <- bx1+dvcb1%*%(chvvcb1%*%rnorm(4)) # Compute the new values of the x1 dataset # tdta1 <- as.matrix(cbind(odraw,ndta[,c(3,4,5)])) %*% b1hat # Draw from the b2 distribution # b2hat <- bx2+dvcb2%*%(chvvcb2%*%rnorm(4)) # Compute the new values of the x2 dataset # tdta2 <- as.matrix(cbind(odraw,ndta[,c(2,4,5)])) %*% b2hat # Set up the original data from the dataset with missing values # nndta <- ndta for (k in 1:nrow(nndta)) { if (is.na(nndta[k,2])) nndta[k,2] <- tdta1[k] if (is.na(nndta[k,3])) nndta[k,3] <- tdta2[k] } # Test for the appropriateness of the drawn values (the Leontief requires all # # all inputs to be positive # t1 <- min(na.omit(nndta[,2])) t2 <- min(na.omit(nndta[,3])) if ((t1 >= 0) & (t2 >= 0)) { # Omit observations that still contain missing values (missing both N and P) # yxdta <- na.omit(nndta) # Set up the new regressions - y vector, vector of ones, and x matrix # oo <- matrix(1,1,nrow=nrow(yxdta)) y <- yxdta[,1] xx <- cbind(oo,yxdta[,2],yxdta[,3],sqrt(yxdta[,2]),sqrt(yxdta[,3]), sqrt(yxdta[,2]*yxdta[,3])) # Compute the covariance matrix for the regression # cxx <- t(xx)%*%xx/nrow(xx) cxy <- t(xx)%*%y/nrow(xx) # Compute condition number rest for the invertability of the matrix # tt <- eigen(cxx)$values t3 <- max(tt)/min(tt) if (t3 > 0) { # Matrix is invertabile increment the counter by one # jcount <- jcount + 1 # Compute the estimate of beta # bb <- solve(cxx) %*% cxy msebb <- as.numeric(t(y - xx %*% bb)%*%(y - xx %*% bb)/(nrow(xx)-ncol(xx))) vcbb <- msebb * solve(t(xx)%*%xx) if (jcount == 1) bobs <- bb else bobs <- cbind(bobs,bb) if (jcount == 1) vbb <- vcbb else vbb <- vbb + vcbb } } } # End of j loop tbobs <- t(bobs) o <- matrix(1,nrow=nrow(tbobs),1) mbobs <- (bobs %*% o)/jcount v1 <- vbb/jcount v2 <- var(tbobs) vt <- v1+(1+1/jcount)*v2 std <- sqrt(diag(vt)) print(cbind(mbobs,std,mbobs/std)) if (i==1) bb1 <- mbobs else bb1 <- cbind(bb1,mbobs) if (i==1) vv1 <- vt else vv1 <- cbind(vv1,vt) } # End of i loop # Test for the condition numbers of vcbb