rm(list=ls())

# Simulating data
set.seed(13579)
n  = 100
x = rnorm(n)
z = rnorm(n)
e = rnorm(n,0,0.5)
a = 1
b = 2
y = 2 + (a*x+b*z)^2 + e

par(mfrow=c(1,3))
plot(x,y)
plot(z,y)
plot(x,z)

reg = lm(y~x+z)
summary(reg)

# RESET test
yfit2 = reg$fit^2
reg1  = lm(y~x+z+yfit2)
summary(reg1)

e      = reg1$res
SSR    = t(e)%*%e
bhat   = reg1$coef
X      = cbind(1,x,z,yfit2)
p      = ncol(X)
iXtX   = solve(t(X)%*%X)

g      = 1
R      = matrix(c(0,0,0,1),g,p,byrow=TRUE)
r      = matrix(0,g,1)
Rbr    = R%*%bhat-r
num   = t(Rbr)%*%solve(R%*%iXtX%*%t(R))%*%(Rbr)/g
den   = SSR/(n-p)
Ftest = num/den
criticalvalue = qf(0.95,g,n-p)
pvalue = 1-pf(Ftest,g,n-p)
c(Ftest,criticalvalue,pvalue)

x2 = x^2
z2 = z^2
xz = x*z
lm(y~x2+z2+xz)