# Kentaro Fukumoto # simulation n <- 100 seed <- 10072013 set.seed(seed) X <- cbind(1,rnorm(n)) beta <- c(-1,1) Xbeta <- X%*%beta iteration <- 1000 sigma <- 2 library(quantreg) beta.OLS <- se.OLS <- beta.LAD <- matrix(NA,iteration,length(beta)) date() for(i in 1:iteration){ set.seed(seed+(i*10000)) Y <- rnorm(n,Xbeta,sigma) OLS.out <- summary(lm(Y~X-1))$coef beta.OLS[i,] <- OLS.out[,1] se.OLS[i,] <- OLS.out[,2] beta.LAD[i,] <- rq(Y~X-1)$coef } date() apply(beta.OLS,2,mean) apply(beta.LAD,2,mean) apply(beta.OLS,2,sd) apply(se.OLS,2,mean) sqrt(diag(solve(t(X)%*%X)))*sigma apply(beta.LAD,2,sd) mean(( beta.OLS[,2]+(se.OLS[,2]*qnorm(0.025))-beta[2])* (beta.OLS[,2]+(se.OLS[,2]*qnorm(0.975))-beta[2])<0) # bootstrap library(faraway) attach(gala) lm.out <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent) covariates <- model.matrix(lm.out)[,-1] B <- 1000 N <- length(Species) boot.out <- matrix(NA,B,length(lm.out$coef)) date() for(b in 1:B){ set.seed(seed+b) resample <- sample(N,rep=T) boot.out[b,] <- lm(Species[resample] ~ covariates[resample,])$coef } date() round( cbind( summary(lm.out)$coef[,1], apply(boot.out,2,mean), summary(lm.out)$coef[,2], apply(boot.out,2,sd) ) ,3)