EX1 P86 3.6 Although the two models’ R^2 are the same,their ESS&TSS are not.
Use your data to verify it.

X <- sort(runif(20))
y <- sort(rnorm(20,5))
beta2 <- sum((X-mean(X))*(y-mean(y)))/(sum((X-mean(X))^2))
beta1 <- mean(y)-beta2*mean(X)
cbind(beta1,beta2)
##      beta1 beta2
## [1,] 3.318 3.155
yhat <- beta1 + beta2*X

beta4 <- sum((y-mean(y))*(X-mean(X)))/(sum((y-mean(y))^2))
beta3 <- mean(X)-beta4*mean(y)
cbind(beta3,beta4)
##        beta3  beta4
## [1,] -0.9915 0.3047
Xhat <- beta3 + beta4*y

ESS1 <- sum((y-mean(y))^2)
TSS1 <- sum((yhat-mean(y))^2)

ESS2 <- sum((X-mean(X))^2)
TSS2 <- sum((Xhat-mean(X))^2)

identical(ESS1,ESS2)
## [1] FALSE
identical(TSS1,TSS2)
## [1] FALSE

Conclusion: ESS1 does not equal to ESS2 TSS1 does not equal to TSS2

EX2 P87 3.14 Suppose each X are mutiplied by 2,will the error term become different? Use your data to verify it.

X1 <- sort(runif(20))
y <- sort(rnorm(20,5))
X2 <- 2*X1
LM1 <- lm(y~X1)
LM2 <- lm(y~X2)
LM1
## 
## Call:
## lm(formula = y ~ X1)
## 
## Coefficients:
## (Intercept)           X1  
##        3.29         3.18
LM2
## 
## Call:
## lm(formula = y ~ X2)
## 
## Coefficients:
## (Intercept)           X2  
##        3.29         1.59
Re1 <- LM1$re
Re2 <- LM2$re
identical(Re1,Re2)
## [1] TRUE

Conclusion: Re1 equals Re2

EX3 P863 C.10 Use specified math methods to calculate instead of using function “lm”

Data <- read.table("C.10.txt", header = TRUE)

prepare the data

Y <- as.matrix(Data[1])
X2 <- Data[2]
X3 <- Data[3]
X <- as.matrix(cbind(1, X2, X3))

coefficients

betaHat <- solve(t(X)%*%X)%*%t(X)%*%Y
betaHat
##          Y
## 1  300.286
## X2   0.742
## X3   8.044

fitted values

YHat <- X %*% betaHat
YHat
##          Y
##  [1,] 1673
##  [2,] 1685
##  [3,] 1683
##  [4,] 1728
##  [5,] 1738
##  [6,] 1766
##  [7,] 1818
##  [8,] 1860
##  [9,] 1950
## [10,] 2042
## [11,] 2122
## [12,] 2181
## [13,] 2250
## [14,] 2294
## [15,] 2346

No.of observations

n <- length(Y) 
n
## [1] 15

No.of parameters

k <- length(betaHat) 
df <- n-k
uHat <- Y-YHat
df
## [1] 12
uHat
##              Y
##  [1,]   0.1674
##  [2,]   3.4140
##  [3,] -16.9838
##  [4,]   6.8735
##  [5,]  11.3460
##  [6,]  -9.7310
##  [7,]  -2.5515
##  [8,]   6.5319
##  [9,]  -2.1296
## [10,]   5.9830
## [11,]   5.9673
## [12,] -15.5309
## [13,]   6.8411
## [14,]  22.1825
## [15,] -22.3800

variance of residuals

sigma2Hat <- sum(uHat^2)/df
sigma2Hat
## [1] 164.7

variance-Covariance of coefficient

vcovbetaHat <- sigma2Hat*solve(t(X) %*% X)
vcovbetaHat
##           1        X2       X3
## 1  6133.650 -3.707941 220.2063
## X2   -3.708  0.002259  -0.1371
## X3  220.206 -0.137052   8.9015

se of coefficients

sebetaHat <- sqrt(diag(vcovbetaHat))
sebetaHat
##        1       X2       X3 
## 78.31763  0.04753  2.98355

t statistic

t <- betaHat/sqrt(diag(vcovbetaHat))
t
##         Y
## 1   3.834
## X2 15.610
## X3  2.696

RSS

RSS <- sum(uHat^2)
RSS
## [1] 1977

ESS

ESS <- sum((YHat-mean(Y))^2)
ESS
## [1] 828144

TSS

TSS <- RSS + ESS
TSS
## [1] 830121

R Squared

R2 <- ESS/TSS
R2
## [1] 0.9976

Adjusted R Squared

adjR2 <- 1- (1-R2)*(n-1)/(n-k)
adjR2
## [1] 0.9972

F statistic

F <- (ESS/(k-1))/(RSS/(n-k)) 
F
## [1] 2514

prediction for 1971 with X2=2610

XNew <- c(1, 2610, 16)
YNew <- XNew %*% betaHat
YNew
##         Y
## [1,] 2366

variance of mean prediction

varMP <- sigma2Hat*t(XNew) %*% solve(t(X) %*% X) %*% XNew
varMP
##       [,1]
## [1,] 48.64

variance of individual prediction

varIP <- sigma2Hat*(1 + t(XNew) %*% solve(t(X) %*% X) %*% XNew)
varIP
##       [,1]
## [1,] 213.4

correlation matrix

cor(cbind(Y, X2, X3))
##         Y     X2     X3
## Y  1.0000 0.9981 0.9743
## X2 0.9981 1.0000 0.9664
## X3 0.9743 0.9664 1.0000