pdf("rpd4_8r.pdf")


alldata <- matrix( 
c(44,89.47,44.609,11.37,62,178,182,
40,75.07,45.313,10.07,62,185,185,
44,85.84,54.297,8.65,45,156,184,
42,68.15,59.571,8.17,40,166,172,
38,89.02,49.874,9.22,55,178,180,
47,77.45,44.811,11.63,58,176,176,
40,75.98,45.681,11.95,70,176,180,
43,81.19,49.091,10.85,64,162,170,
44,81.42,39.442,13.08,63,174,176,
38,81.87,60.055,8.63,48,170,186,
44,73.03,50.541,10.13,45,168,168,
45,87.66,37.388,14.03,56,186,192,
45,66.45,44.754,11.12,51,176,176,
47,79.15,47.273,10.6,47,162,164,
54,83.12,51.855,10.33,50,166,170,
49,81.42,49.156,8.95,44,180,185,
51,69.63,40.836,10.95,57,168,172,
51,77.91,46.672,10,48,162,168,
48,91.63,46.774,10.25,48,162,164,
49,73.37,50.388,10.08,67,168,168,
57,73.37,39.407,12.63,58,174,176,
54,79.38,46.08,11.17,62,156,176,
52,76.32,45.441,9.63,48,164,166,
50,70.87,54.625,8.92,48,146,186,
51,67.25,45.118,11.08,48,172,172,
54,91.63,39.203,12.88,44,168,172,
51,73.71,45.79,10.47,59,186,188,
57,59.08,50.545,9.93,49,148,160,
49,76.32,48.673,9.4,56,186,188,
48,61.24,47.92,11.5,52,170,176,
52,82.78,47.467,10.5,53,170,172),byrow=T,ncol=7)

Y <- alldata[,3]
X <- alldata[,4:7]

X <- cbind(rep(1,nrow(X)),X)

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

BETAHAT <- XPXINV %*% t(X) %*% Y

BETAHAT
XPXINV

ix <- diag(1,ncol=nrow(X),nrow=nrow(X))
jx <- matrix(rep(1,nrow(X)*nrow(X)),ncol=nrow(X)) 

SSRESID <- t(Y) %*% (ix - X %*% XPXINV %*% t(X)) %*% Y
DFRESID <- nrow(X)-ncol(X)
S2 <- SSRESID/DFRESID

K1 <- matrix(c(0,0,1,0,0,0,0,0,0,1),byrow=F,ncol=2)

m1 <- matrix(c(0,0),ncol=1)

Q1 <- t(t(K1)%*%BETAHAT-m1) %*% solve(t(K1) %*% XPXINV %*% K1) %*% (t(K1)%*%BETAHAT-m1)

F1 <- (Q1/ncol(K1))/S2[1,1]

FCRIT1 <- qf(.95,ncol(K1),nrow(X)-ncol(X))

PVAL1 <- 1- pf(F1,ncol(K1),nrow(X)-ncol(X))

Q1
F1
FCRIT1
PVAL1

SSTOTALC <- t(Y) %*% (ix - (1/nrow(X))*jx) %*% Y
DFTOTALC <- nrow(X)-1
MSTOTALC <- NA

SSREG <- t(Y) %*% (X%*%XPXINV%*%t(X) - (1/nrow(X))*jx) %*% Y
DFREG <- ncol(X)-1
MSREG <- SSREG/DFREG

print(cbind(DFTOTALC,SSTOTALC,MSTOTALC))
print(cbind(DFREG,SSREG,MSREG))
print(cbind(DFRESID,SSRESID,S2))



K2 <- matrix(c(1,0,0,0,0,0,0,1,0,0,0,0,0,0,1),byrow=F,ncol=3)

m2 <- matrix(c(90,0,0),ncol=1)

K2BM <- t(K2)%*%BETAHAT-m2
K2XPXIK2I <- solve(t(K2) %*% XPXINV %*% K2)


Q2 <- t(K2BM) %*% K2XPXIK2I %*% K2BM

F2 <- (Q2/ncol(K2))/S2[1,1]

FCRIT2 <- qf(.95,ncol(K2),nrow(X)-ncol(X))

PVAL2 <- 1- pf(F2,ncol(K2),nrow(X)-ncol(X))

K2BM
K2XPXIK2I
Q2
F2
FCRIT2
PVAL2


dev.off()