rm(list = ls())                              # Διαγράφει όλα τα αντικείμενα από το current R workspace
cat("\014")                                  # Καθαρίζει την κονσόλα της R
graphics.off()                              # Κλείνει όλα τα ανοιχτά γραφήματα


# install.packages("wooldridge")
# install.packages("lmtest")
# install.packages("estimatr")
# install.packages("orcutt")
# install.packages("sandwich")

# install.packages("remotes")
# remotes::install_version("orcutt", version = "2.3", repos = "https://cran.r-project.org")

library(wooldridge)
library(lmtest)
library(estimatr)

library(orcutt)

library(sandwich)


data("vote1")
data("phillips")
data("fish")


#*********** Chow test ***********************************
#*******************************************************
#*******************************************************

rm(list=ls())  
library(wooldridge)
data("sleep75")
#  help(sleep75)

names(sleep75)

#********** Different ols for men and women **************

ols1a=lm(sleep ~ totwrk + educ + age + agesq + yngkid ,data=sleep75, subset=male==1 )
summary(ols1a)

ols1b=lm(sleep ~ totwrk + educ + age + agesq + yngkid ,data=sleep75, subset=male==0 )
summary(ols1b)


#** plotting the quadratic functon of age for men and women ***

# δημιουργεί ένα διάνυσμα με όλες τις ακέραιες ηλικίες από τη μικρότερη μέχρι τη μεγαλύτερη ηλικία που υπάρχει στο dataset sleep75.
years = min(sleep75$age) : max(sleep75$age) 

plot( years , ols1a$coeff[4]*years + ols1a$coeff[5]*years ^2 , type="l", ylab="min of sleep", col="blue", ylim=c(-650,280))
lines(years , ols1b$coeff[4]*years + ols1b$coeff[5]*years ^2 , col="green")
legend("right", c("men","women"), fil=c("blue","green"))

-ols1b$coeff[4]/(2*ols1b$coeff[5])  # global minimum for women quadratic relation

#*********** Chow test ***********************************


# male:educ = male * educ, Το male:μεταβλητη στην R σημαίνει όρος αλληλεπίδρασης (interaction term).

ols2a=lm(sleep ~ totwrk + educ + age + agesq + yngkid + male  + male:totwrk + male:educ 
         +male:age + male:agesq + male:yngkid , data=sleep75)
summary(ols2a)

ols2b=lm(sleep ~ totwrk + educ + age + agesq + yngkid , data=sleep75)
summary(ols2b)

#   H0: b6=b7=b8=b9=b10=b11=0
#   H1: bi<>0 for some i
anova(ols2a,ols2b)  #*********** Chow test results


#********* "educ" for men ***************
ols1a$coeff[3] 
ols2a$coeff[3] + ols2a$coeff[9]

#******* and "educ" for women **********
ols1b$coeff[3] 
ols2a$coeff[3] 

# *********** Second Chow test *********************************

ols3b=lm(sleep ~ totwrk + educ + age + agesq + yngkid + male , data=sleep75)
summary(ols3b)

#   H0: b7=b8=b9=b10=b11=0
#   H1: bi<>0 for some i

anova(ols2a,ols3b)  #****** Second Chow test results

#*** considering:
anova(ols2a,ols2b)

#***  and:
anova(ols2a,ols3b)

#** the best model is:
summary(ols3b)
#*********************************************************
#*********************************************************
#*********************************************************
#
#*********************************************************
#********* HETEROSCEDASTICITY TESTS *********************
#*********************************************************
rm(list = ls())                              # Διαγράφει όλα τα αντικείμενα από το current R workspace
cat("\014")                                  # Καθαρίζει την κονσόλα της R
graphics.off()                              # Κλείνει όλα τα ανοιχτά γραφήματα


#
#*********************************************************
#*********** Breusch–Pagan and White tests *******************
#*********************************************************

library(wooldridge)
# help(vote1)

names(vote1)


#***************** simple OLS regression *****************

ols1=lm(voteA ~ prtystrA + democA + log(expendA) + log(expendB) , data=vote1)
summary(ols1)


#** OLS of residuals with independent variables **********

ols2=lm(ols1$res ~ prtystrA + democA + log(expendA) + log(expendB) , data=vote1)
summary(ols2)


#** Breusch–Pagan heteroskedasticity  test ***************

plot(ols1$fit,ols1$res,xlab="Fitted",ylab="Residuals", col="blue")   


#install.packages("lmtest")
library(lmtest)

bptest(ols1)   


#**** Manually estimating the Breusch-Pagan test *********

BP.ols=lm(ols1$res^2 ~ prtystrA + democA + log(expendA) + log(expendB) , data=vote1)
summary(BP.ols)

BP_stat= nrow(vote1) * summary(BP.ols)$r.squared
BP_stat

BP.p.value = 1-pchisq(BP_stat, 4)  #*** 4=number of regressors in BP.ols ****
BP.p.value

#**** White's  heteroskedasticity  test  *****************

bptest(ols1,~( prtystrA + democA +  log(expendA) + log(expendB) )^2 + I(prtystrA^2) + 
         I((log(expendA))^2) + I((log(expendB))^2),data=vote1)


#**** Manually estimating the White test *****************

White.ols=lm(ols1$res^2 ~ (prtystrA + democA +  log(expendA) + log(expendB))^2 + I(prtystrA^2) + 
               I(lexpendA^2) + I(lexpendB^2),data=vote1)

summary(White.ols)

White.stat= nrow(vote1)*summary(White.ols)$r.squared
White.stat

White.p.value = 1-pchisq(White.stat, 13)  #*** 13=number of regressors in White.ols ****
White.p.value

#****************** Robust regression (for obtaining robust standard errors) **************

#install.packages("estimatr")
library(estimatr)

ols1_robust = lm_robust(voteA~prtystrA + democA + log(expendA) + log(expendB) , data=vote1)
summary(ols1_robust)
#*********************************************************
#*********************************************************
#*********************************************************
#
#
#*************** AUTOCORELLATION TESTS *******************
#*
# H0: ρ1 = 0
# Η1: ρ1 =/ 0 ή ρ1>0 ή ρ1<0
#
#*********************************************************
#*********** Darbin Watson test *****************************
#*********************************************************
rm(list=ls())  
library(wooldridge)
# help(phillips)

names(phillips)

#***************** simple OLS regression *****************

ols1=lm(inf ~ unem , data=phillips)
summary(ols1)

acf(ols1$res, col="red")

#******************  DW test *****************************
library(lmtest)

dwtest(ols1)

# ΚΡΙΤΙΚΕΣ ΤΙΜΕΣ ΓΙΑ α=5%, n=56, dL=1.528 KAI dU=1.601 ==> DW = 0.80148<  dL=1.528 => Η0 απορρίπτεται

#********  Cochrane-Orcutt estimation ********************
#install.packages("orcutt")
library(orcutt)

summary(cochrane.orcutt(ols1))


#********* Newey-West HAC estimators *********************
#*********************************************************
rm(list=ls())  
library(wooldridge)

#help(fish)
names(fish)


ols1=lm(avgprc ~ wave2 , data=fish)
summary(ols1)

plot(ols1$fit,ols1$res,xlab="Fitted",ylab="Residuals", col="blue")   

#**** Breusch–Pagan  heteroskedasticity  test  ***********
library(lmtest)
bptest(ols1)

#**** White's  heteroskedasticity  test  *****************
bptest(ols1,~ wave2 + I(wave2^2) , data=fish)

acf(ols1$res, col="red")

#**** Durbin-Watson  test  *******************************
dwtest(ols1)


#install.packages("sandwich")
library(sandwich)

# COVARIANCE MATRIX

NW.estimators=NeweyWest(ols1, lag = 1)
NW.estimators

sqrt(NW.estimators[1,1])
sqrt(NW.estimators[2,2])

coeftest(ols1, vcov=NW.estimators)

coefci(ols1, vcov = NW.estimators)