data = read.table("hprice1.txt",header=TRUE)

attach(data)


# ----------
# Full model
# ----------
X = cbind(1,lassess,lotsize,sqrft,bdrms)

y = lprice

p = ncol(X)

n = nrow(X)

iXtX    = solve(t(X)%*%X)

bhat    = iXtX%*%t(X)%*%y

ehat    = y-X%*%bhat

s2      = t(ehat)%*%ehat/(n-p)

s       = sqrt(s2)

se.bhat = s*sqrt(diag(iXtX))

round(cbind(bhat,se.bhat,bhat/se.bhat,2-2*pnorm(bhat/se.bhat)),8)

# -------------------------------------------------
# Restricted model 2: H0: beta1=beta2=beta3=beta4=0
# -------------------------------------------------
R = cbind(0,diag(1,4))

r = matrix(0,4,1)

num = t(R%*%bhat - r)%*%solve(R%*%iXtX%*%t(R))%*%(R%*%bhat - r)

den = sum((y-mean(y))^2)/n

LM  = num/den

LM

criticalvalue = qgamma(0.95,4/2,1/2)

criticalvalue

# Alternative solution
# --------------------
res      = y-mean(y)
bhat1    = iXtX%*%t(X)%*%res
ehat1    = res-X%*%bhat1
s21      = t(ehat1)%*%ehat1/(n-p)
s1       = sqrt(s21)
se.bhat1 = s*sqrt(diag(iXtX))

round(cbind(bhat1,se.bhat1,bhat1/se.bhat1,2-2*pnorm(bhat1/se.bhat1)),8)

R2  = 1-sum(ehat1^2)/sum(res^2)

LM1 = n*R2



# -------------------------------
# Restricted model 2: H0: beta1=1
# -------------------------------
R = matrix(c(0,1,0,0,0),1,p)
r = 1

num = t(R%*%bhat - r)%*%solve(R%*%iXtX%*%t(R))%*%(R%*%bhat - r)

y1  = y - assess
X1 = cbind(1,lotsize,sqrft,bdrms)
p1 = ncol(X1)
iXtX1    = solve(t(X1)%*%X1)
bhat1    = iXtX1%*%t(X1)%*%y1
ehat1    = y1-X1%*%bhat1
s21      = t(ehat1)%*%ehat1/(n-p1)

LM  = num/s21

LM

criticalvalue = qgamma(0.95,1/2,1/2)

criticalvalue