# 
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

library(alr3)
data(longley)
attach(longley)

# 4.10.1 a scatter plot matrix:
#
postscript("../../WriteUp/Graphics/Chapter4/prob_10.eps", onefile=FALSE, horizontal=FALSE)
pairs(Def~GNP+Unemp+Mil+Pop+Emp+Year)
dev.off()

# 4.10.2 fit a regresion on these variables (not Year)
#
m = lm( Def~GNP+Unemp+Mil+Pop+Emp, data=longley )
s <- summary(m) 
print(coef(m),digits=4)

# get the standard errors for this error free regression:
print( s$coefficients[,2], digits=4 ) 

# 4.10.3 add random noise and compute the new regression coefficients 
#
longley.sim <- function(B=999){
 ans <- NULL
 for (i in 1:B) {
   Def   <- longley$Def   + 0.1 * runif(16,-1,+1) 
   GNP   <- longley$GNP   + 500 * runif(16,-1,+1) 
   Unemp <- longley$Unemp + 5   * runif(16,-1,+1) 
   Mil   <- longley$Mil   + 5   * runif(16,-1,+1) 
   Pop   <- longley$Pop   + 500 * runif(16,-1,+1) 
   Emp   <- longley$Emp   + 50  * runif(16,-1,+1) 
   m0    <- lm(Def~GNP+Unemp+Mil+Pop+Emp)
   ans   <- rbind(ans,coef(m0))
  }
  ans}
cl  <- function(x) quantile(x,c(.025,.975))
std <- function(x) sqrt(var(x)) 

b0 <- longley.sim(999)

apply(b0,2,mean)
#apply(b0,2,cl)
print( apply(b0,2,std), digits=4 )