datamat = t(matrix(scan("C:/HY-Data/JPLUOTO/OPETUS/EKONOMETRIAN_PERUSKURSSI/test_data.txt"),nrow=13))
datam = matrix(datamat,nc=13)
colnames(datam) = c("enrl_tot","teachers","calw_pct","meal_pct","computer","testscr","comp_stu","expn_stu","str","avginc","el_pct","read_scr","math_scr")


#Ladataan tarvittavat paketit
library(sandwich)
library(lmtest)
library(car)



Test_score = datam[,"testscr"]
STR = datam[,"str"]
PctEL= datam[,"el_pct"]
Expn = datam[,"expn_stu"]



################################
###                          ###
### Hypoteesin testaus S&W 7 ###
###                          ###
################################


model = lm(Test_score~STR+PctEL)
model$newse = vcovHC(model,type="HC1")
coeftest(model,model$newse)
2*pnorm(-abs(-1.101294/0.432847)) # Vert.
summary(model)



# H0: beta_STR = 0, beta_PctEL = 0 ja beta_Expn = 0

model_unrestricted = lm(Test_score~STR+PctEL+Expn)
betahat = model_unrestricted$coeff
model_restricted = lm(Test_score~1)
n = dim(datam)[1]
q = 3
k = length(betahat)-1
R2unrestricted = summary(model_unrestricted)$r.squared
R2restricted = summary(model_restricted)$r.squared
F = ((R2unrestricted-R2restricted)/q)/((1-R2unrestricted)/(n-k-1) )
#p-arvo
pf(F, df1 = q, df2 = Inf, lower.tail = FALSE, log.p = FALSE)
#Vertaa
summary(model_unrestricted)
F


# H0: beta_STR = 0 ja beta_Expn = 0 vert. S&W sivu 267-268

model_unrestricted = lm(Test_score~STR+PctEL+Expn)
model_restricted = lm(Test_score~PctEL)
n = dim(datam)[1]
q = 2
R2unrestricted = summary(model_unrestricted)$r.squared
R2restricted = summary(model_restricted)$r.squared
F = ((R2unrestricted-R2restricted)/q)/((1-R2unrestricted)/(n-k-1) )
#p-arvo
pf(F, df1 = q, df2 = Inf, lower.tail = FALSE, log.p = FALSE)
F







# Lineaariset rajoitteet
# H0: R*beta = r 
# H0 -> R*betahat ~ N(r,R*COVbetahat*R')



# Esimerkki 1
# H0: beta_STR = 0 ja beta_Expn = 0

model_unrestricted = lm(Test_score~STR+PctEL+Expn)
betahat = model_unrestricted$coeff
k = length(betahat)-1
q = 2
r = rep(0,q)
R = matrix(0,nc=k+1,nr=q)
R[1,2] = 1
R[2,4] = 1

model$newse = vcovHC(model_unrestricted,type = "const")                     # Kokeile type = "HC1" ja vertaa S&W sivu 268
COVbetahat = model$newse
F = t(R%*%betahat-r)%*%solve(R%*%COVbetahat%*%t(R))%*%(R%*%betahat-r)/q
pf(F, df1 = q, df2 = Inf, lower.tail = FALSE, log.p = FALSE)              
F

# Voit käyttää myös suoraan car-pakettia
linearHypothesis(model_unrestricted,R,r,white.adjust="hc1")



# Esimerkki 2
# H0: beta_STR = beta_Expn

model_unrestricted = lm(Test_score~STR+PctEL+Expn)
betahat = model_unrestricted$coeff
k = length(betahat)-1
q = 1
r = rep(0,q)
R = matrix(0,nc=k+1,nr=q)
R[1,2] = 1
R[1,4] = -1

linearHypothesis(model_unrestricted,R,r,white.adjust="hc1")

# Tai
model$newse = vcovHC(model_unrestricted,type = "HC1")
COVbetahat = model$newse
F = t(R%*%betahat-r)%*%solve(R%*%COVbetahat%*%t(R))%*%(R%*%betahat-r)/q
pf(F, df1 = q, df2 = Inf, lower.tail = FALSE, log.p = FALSE)

#Vert. edellä laskettua p-arvoa seuraavaan 
model = lm(Test_score~STR+c(STR+Expn)+PctEL) #Uudelleen järjestellään regressio (S&W sivu 269)
model$newse = vcovHC(model,type = "HC1")
coeftest(model,model$newse)
2*pnorm(-abs(-0.2902382/0.4812497)) #Poimitaan STR:n kerroin ja sen keskivirhe ja lasketaan p-arvo

#vertaa myös
F
(-0.2902382/0.4812497)^2