####################
#MATH 5305 Examples#
#Chapters 2 and 4  #
####################

#Generating Example Data
math=rnorm(1000,mean=500,sd=100)
hist(math)

epsilon = rnorm(1000,0,86.6)

verbal = 250 + 0.5*math + epsilon
beta1=250
beta2=0.5
sigma = 86.6


#Basic Exploratory Data Analysis
mean(math)
sd(math)
hist(math)

mean(verbal)
sd(verbal)
hist(verbal)

cor(math,verbal)
plot(math,verbal)


#Fitting regression line using chapter 2 methods

beta2hat = cor(math,verbal)*sd(verbal)/sd(math)
beta2hat


beta1hat= mean(verbal)-beta2hat*mean(math)
beta1hat

plot(math,verbal)
lines(0:1000,beta1hat+beta2hat*(0:1000),col="red")


#Fitting regression line using chapter 4 methods

X=cbind(rep(1,1000),math)
X

head(X)
head(math)

Y=verbal

betahat=solve(t(X)%*%X)%*%t(X)%*%Y
betahat

beta1hat
beta2hat


#Fitting regression line using lm command

model=lm(verbal~math)
summary(model)

coef(model)



#Yhat and residuals e
Yhat = X%*%betahat
e = verbal - Yhat



#All variables for first student.
math[1]
verbal[1]
epsilon[1]

#True model for first student
250+0.5*math[1]
epsilon[1]

250+0.5*math[1]+epsilon[1]
verbal[1]


#We can't see the true model.  What does the model that we fit predict?
beta1hat
beta2hat

beta1hat+beta2hat*math[1]     #Model prediction for 1st student (fitted value)
Yhat[1]

verbal[1]                     #Actual value of Y for 1st student

verbal[1]-yhat.firststudent   #Residual for 1st student
e[1]
epsilon[1]

#True Model
beta1
beta2
epsilon[1]

beta1+beta2*math[1]+epsilon[1]
verbal[1]

#Fitted Model
beta1hat
beta2hat
e[1]

beta1hat+beta2hat*math[1]+e[1]
verbal[1]


#Predictions and residuals using lm
head(Yhat)
head(e)

summary(model)
names(model)

head(cbind(e,
           model$residuals))

head(cbind(Yhat,
           model$fitted.values))



#Estimating Sigma

sigma

sigmahat = sqrt(sum(e^2)/(1000-2))
sigmahat

summary(model)


#Estimating cov(betahat)

covhat=sigmahat^2*solve(t(X)%*%X)
covhat

#Standard Error for beta1hat
sqrt(covhat[1,1])
sqrt(covhat[2,2])