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)


#Nimetään muuttujat
Test_score = datam[,"testscr"]
STR = datam[,"str"]
PctEL= datam[,"el_pct"]
Expn = datam[,"expn_stu"]
Income = datam[,"avginc"]




################################
###                          ###
### 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)




# 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







######################################
###                                ###
### Epälineaarinen regressio S&W 8 ###
###                                ###
######################################

#Polynomit


plot(Income,Test_score,xlim=c(0,60),ylim=c(600,740))

model = lm(Test_score~Income+c(Income^2))
model$newse = vcovHC(model,type="HC1")
coeftest(model,model$newse)

beta0 = summary(model)$coeff[1]
beta1 = summary(model)$coeff[2]
beta2 = summary(model)$coeff[3]

abline(lm(Test_score~Income),lty=2)
fit = beta0 + beta1*c(0:55) + beta2*c(0:55)^2 
lines(fit,lty=1)


#dY/dX
dTest_score = beta0 + beta1*6 + beta2*6^2 - (beta0 + beta1*5 + beta2*5^2)
dTest_score
dTest_score = beta0 + beta1*26 + beta2*26^2 - (beta0 + beta1*25 + beta2*25^2)
dTest_score
dTest_score = beta0 + beta1*36 + beta2*36^2 - (beta0 + beta1*35 + beta2*35^2)
dTest_score




model = lm(Test_score~Income+c(Income^2)+c(Income^3))
model$newse = vcovHC(model,type="HC1")
coeftest(model,model$newse)

beta0 = summary(model)$coeff[1]
beta1 = summary(model)$coeff[2]
beta2 = summary(model)$coeff[3]
beta3 = summary(model)$coeff[4]

fit = beta0 + beta1*c(0:55) + beta2*c(0:55)^2 + beta3*c(0:55)^3
lines(fit,lty=3)



#H0: lineaarinen regression (beta2 = beta3 = 0)
betahat = c(beta0,beta1,beta2,beta3)
k = 3
q = 2
r = rep(0,q)
R = matrix(0,nc=k+1,nr=q)
R[1,3] = 1
R[2,4] = 1

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

# Tai
newse = vcovHC(model,type = "HC1")
COVbetahat = 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)