# demonstration of weighted least squares

# First we specify the way that the variance depends on the mean

meanvarfunc <- function(mu,delta){

return(mu^2 * delta)

}

# Set the number of observations in the simulation

numobs <- 101

# Set the number of repliactions in the simulation

numreps <- 2000

# We will fake data from a simple linear regression with parameters

# beta0 and beta1, but the noise wil be heteroscedastic.

# The x-values will be between 10 and 20.

x <- seq(10,90,length=numobs)

# Now a function to generate the data from the regression

fakedata <- function(x,beta0,beta1,delta){

goupline <- beta0 + beta1 * x

noise <- sqrt(meanvarfunc(x,delta)) * rnorm(numobs)

return(goupline + noise)

}

# Time for a simulation. Plan:

# 1. Run a regression.

# 2. Read off estimate and standard error

# 3. Calculate mean and variance of the estimates

# 4. Check to see if CI contains true value

# The first is the OLS simulation.

# Make a matrix to hold all the results

results.ols <- matrix(rep(0,3*numreps),ncol=3)

for( i in 1:numreps){

simdata <- fakedata(x,10,2,.25)

lm.ols <- lm(simdata ~ x)

results.ols[i,1] <- coef(lm.ols)[2]

results.ols[i,2] <- summary(lm.ols)$coefficients[2,2]

if((2 < (results.ols[i,1] + 1.96 * results.ols[i,2])) &&

(2 > (results.ols[i,1] - 1.96 * results.ols[i,2]))){ results.ols[i,3] <- 1}

print(paste("OLS Iteration",i))

}

mean(results.ols[,1]) # This finds the average estimate

mean(results.ols[,2]^2) # This finds the average estimated s.e.

var (results.ols[,1]) # This finds the "true" s.e.

# Now the simulation using the bootstrap standard error.

#

results.boot <- matrix(rep(0,3*numreps),ncol=3)

for( i in 1:numreps){

simdata <- fakedata(x,10,2,.25)

simdataframe <- as.data.frame(cbind(simdata,x))

lm.boot <- lm(simdata ~ x)

results.boot[i,1] <- coef(lm.boot)[2]

boot.out <- bootstrap(simdataframe,

coef(lm(simdata ~ x,data=simdataframe))[2],

B=100)

results.boot[i,2] <- boot.out$estimate[[3]]

if((2 < (limits.emp(boot.out)[4])) &&

(2 > (limits.emp(boot.out)[1])))

{ results.boot[i,3] <- 1}

print(paste("BOOT Iteration",i))

}

mean(results.boot[,1]) # This finds the average estimate

mean(results.boot[,2]^2) # This finds the average estimated s.e.

var (results.boot[,1]) # This finds the "true" s.e.

# Finally the WLS simulation.

results.wls <- matrix(rep(0,3*numreps),ncol=3)

for( i in 1:numreps){

simdata <- fakedata(x,10,2,.25)

w <- 1/meanvarfunc(x,.25)

lm.wls <- lm(simdata ~ x,weights=w)

results.wls[i,1] <- coef(lm.wls)[2]

results.wls[i,2] <- summary(lm.wls)$coefficients[2,2]

if((2 < (results.wls[i,1] + 1.96 * results.wls[i,2])) &&

(2 > (results.wls[i,1] - 1.96 * results.wls[i,2]))){ results.wls[i,3] <- 1}

print(paste("WLS Iteration",i))

}

mean(results.wls[,1]) # This finds the average estimate

mean(results.wls[,2]^2) # This finds the average estimated s.e.

var (results.wls[,1]) # This finds the "true" s.e.