pdf("rpd4_1r.pdf")

X <- matrix(c(1,0.02,1,0.07,1,0.11,1,0.15),byrow=T,ncol=2)

Y <- matrix(c(242,237,231,201),ncol=1)

XPX <- t(X) %*% X
XPY <- t(X) %*% Y
YPY <- t(Y) %*% Y
XPXINV <- solve(XPX)

BETAHAT <- XPXINV %*% XPY

P <- X %*% XPXINV %*% t(X)

YHAT <- P %*% Y

I4 <- diag(1,nrow=4,ncol=4)

e <- (I4-P) %*% Y

SSTOTALU <- YPY
SSMODEL <- t(YHAT) %*% YHAT
SSRESID <- t(e) %*% e

J4 <- matrix(rep(1,16),ncol=4,nrow=4)

SSMU <- t(Y) %*% (J4/4) %*% Y
SSTOTALC <- SSTOTALU-SSMU
SSREG <- (t(BETAHAT) %*% XPY)-SSMU

S2 <- SSRESID/2

S2BETAHAT <- S2[1,1]*XPXINV

MSREG <- SSREG/1
MSRESID <- SSRESID/2

FSTAT <- MSREG/MSRESID
FCRIT <- qf(.95,1,2)
PROBF <- 1-pf(FSTAT,1,2)

SSTOTALU
SSMU
SSTOTALC
SSREG
SSRESID

BETAHAT
S2BETAHAT

FSTAT
FCRIT
PROBF

dev.off()