# Homework Assignment 3: Old Heteroskedastic # 36-402, Data Analysis, Spring 2011 # Solutions by Z. Kurtz # Setup setwd("C:/Documents and Settings/zach/My Documents/My Dropbox/402") library(MASS) data(geyser) summary(geyser) ############################################################ # Problem 1 #pdf("hw3/graphics/p1.pdf") # Plot the data points plot(geyser\$duration, geyser\$waiting, cex = 0.5, pch = 16, main = "Waiting time as function of geyser duration", xlab = "Duration (minutes)", ylab = "Waiting time (minutes)") mtext("Black line = LS regression line") # Build linear model: lm.1 = lm(waiting ~ duration, data = geyser) # Add the regression line to the plot abline(lm.1) #dev.off() ############################################################ # Problem 2 #pdf("hw3/graphics/p2and3.pdf") # Plot the squared residuals against duration plot(geyser\$duration, lm.1\$residuals^2, cex = 0.5, pch = 16, main = "Squared residuals versus geyser duration", xlab = "Geyser duration", ylab = "Squared residuals from linear model") ############################################################ # Problem 3 # To estimate the variance, I will use the same kernel regression function # from HW Assignment 2: require(np) reg.var3 = npreg(lm.1\$residuals^2 ~ geyser\$duration) # Pair each observation, indexed by duration, with its fitted value fits = cbind(geyser\$duration, reg.var3\$mean) # Order the observation-fit pairs by duration (column 1) fits3 = fits[order(fits[,1]),] # Add the fits to the plot lines(fits3, col=2, lwd = 2.5) mtext("Red curve: kernel regression function", col = 2) #dev.off() ############################################################ # Problem 4 # Build a weighted linear model: lm.4 = lm(waiting ~ duration, data = geyser, weights = 1/reg.var3\$mean) # Consider the coefficients and standard errors summary(lm.1)\$coefficients summary(lm.4)\$coefficients ############################################################ # Problem 5 #pdf("hw3/graphics/p5.pdf") # Do a nonparametric kernel regression of waiting on duration regression5 = npreg(waiting ~ duration, data = geyser, residuals = T) # Plot the data and then plot the regression results plot(geyser\$duration, geyser\$waiting, cex = 0.5, pch = 16, main = "Waiting time as function of geyser duration with a kernel regression curve", xlab = "Duration (minutes)", ylab = "Waiting time (minutes)") # Pair each observation, indexed by duration, with its fitted value fits = cbind(geyser\$duration, regression5\$mean) # Order the observation-fit pairs by duration (column 1) fits5 = fits[order(fits[,1]),] # Add the fits to the plot lines(fits, col=2, lwd = 2.5) #dev.off() ############################################################ # Problem 6 # Do a kernel regression of the squared residuals (from part 5) on duration reg.var6 = npreg(regression5\$resid^2 ~ geyser\$duration) # Pair each observation, indexed by duration, with its fitted value fits = cbind(geyser\$duration, reg.var6\$mean) # Order the observation-fit pairs by duration (column 1) fits6 = fits[order(fits[,1]),] #pdf("hw3/graphics/p6.pdf") # Plot the fits as a blue dashed curve plot(fits6, type = "l", col = 4, lty = 2, lwd = 2.5, xlab = "Duration (minutes)", ylab = "Estimated variance of waiting given duration", main = "Estimated variance of waiting given duration: Solid black is for unweighted linear regression, and dashed blue is for kernel regression" ) # Add to the plot the fits from part 3 lines(fits3, lwd = 2) #dev.off() ############################################################ # Problem 7 require(np) #pdf("hw3/graphics/p7.pdf") # Produce a bandwidth specification bw = npcdensbw(waiting~duration, data = geyser, tol=.1, ftol=.1) # Produce a conditional density estimate and plot it npc.dens7 = npcdens(bws = bw) npplot(bws=bw, theta = -30, phi = 30, view = "fixed", main = "Conditional density estimate for waiting given duration") #dev.off()