rm(list=ls())
library(rio)
acs <- import("https://github.com/marctrussler/IIS-Data/raw/main/ACSCountyData.csv")

cor(acs$population.density, acs$percent.transit.commute, use="pairwise.complete")

acs$density[acs$population.density<50] <- 1
acs$density[acs$population.density>=50 & acs$population.density<1000] <- 2
acs$density[acs$population.density>=1000] <- 3

table(acs$density)

plot(acs$density, acs$percent.transit.commute)

boxplot(acs$percent.transit.commute ~ acs$density)

boxplot(acs$percent.transit.commute ~ acs$density, outline=F)

m <- lm(percent.transit.commute ~ density, data=acs)
summary(m)

plot(acs$density, acs$percent.transit.commute, xlim=c(0,4))
abline(m, col="firebrick", lwd=2)

m1 <- lm(median.income ~ percent.transit.commute, data=acs)
summary(m1)
plot(acs$percent.transit.commute, acs$median.income)
abline(lm(median.income ~ percent.transit.commute, data=acs), lwd=2, col="firebrick")

table(acs$census.region)

acs$midwest[acs$census.region=="midwest"]<- 1
acs$midwest[acs$census.region!="midwest"]<- 0
acs$northeast[acs$census.region=="northeast"]<- 1
acs$northeast[acs$census.region!="northeast"]<- 0
acs$south[acs$census.region=="south"]<- 1
acs$south[acs$census.region!="south"]<- 0
acs$west[acs$census.region=="west"]<- 1
acs$west[acs$census.region!="west"]<- 0

#Check if this worked:
table(acs$census.region, acs$midwest)
table(acs$census.region, acs$northeast)
table(acs$census.region, acs$south)
table(acs$census.region, acs$west)

#acs[,29:32] <- NULL
#
##install.packages("fastDummies")
#library(fastDummies)
#acs <- dummy_cols(acs, select_columns = "census.region")
#
##Doesn't give them great names, but that's easy to fix:
#names(acs)[29:32] <- c("west","midwest", "northeast","south")

m <- lm(median.income ~ midwest + northeast + south + west, data=acs)
summary(m)

counties <- c(1055,4007,9001,17173)

c <- acs[acs$county.fips %in% counties,]
keep <- c("county.name","state.full","census.region","midwest","northeast","south","west")
c <- c[keep]
c

m <- lm(median.income ~  northeast + south + west, data=acs)
summary(m)

mean(acs$median.income[acs$census.region=="midwest"])

coef(m)["(Intercept)"] + coef(m)["northeast"]
mean(acs$median.income[acs$census.region=="northeast"])

m <- lm(median.income ~  midwest + south + west, data=acs)
summary(m)

m1 <- lm(median.income ~ percent.transit.commute, data=acs)
summary(m1)

m2 <- lm(median.income ~ percent.transit.commute + northeast + south + west, data=acs)
summary(m2)

acs$census.region <- as.factor(acs$census.region)

m <- lm(median.income ~ percent.transit.commute + census.region, data=acs)
summary(m)

acs$census.region <- relevel(acs$census.region, ref="northeast")


m <- lm(median.income ~ percent.transit.commute + census.region, data=acs)
summary(m)


m1 <- lm(median.income ~ percent.transit.commute, data=acs)
summary(m1)

m2 <- lm(median.income ~ percent.transit.commute + northeast + south + west, data=acs)
summary(m2)

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

coefplot(m2)

library(stargazer)

#The "text" type will output into R in a nice looking format:
stargazer(m1, m2, type="text")


stargazer(m1, m2, type="html", style="apsr",
          dep.var.labels = "Median Income",
          covariate.labels = c("% Transit Commuters",
                               "Region: NorthEast",
                               "Region: South",
                               "Region: West"))



set.seed(19104)
#X exists:
x <- rnorm(100, mean=50, sd=3)
#Y data generating process:
y <- 4 + 3*x + rnorm(100,mean=0, sd=2)
plot(x,y)
#Everytime we runw we get something new

y <- rnorm(100,50,25)
x <- rep(0,100)
plot(x,y)
lm(y ~ x)

m <- lm(percent.transit.commute ~ median.income, data=acs)

#Homoskedastic Errors
x <- runif(100, min=0, max=100)
y <- 2 + 3*x + rnorm(100, mean=0, sd=15)

plot(x,y)
abline(lm(y~x), col="firebrick" ,lwd=2)

#Heteroskedastic Errors
x <- runif(100, min=0, max=100)
y <- 2 + 3*x + rnorm(100, mean=0, sd=x*2)

plot(x,y)
abline(lm(y~x), col="firebrick" ,lwd=2)


betas <- rep(NA,10000)
ses <- rep(NA,10000)

for(i in 1:10000){
x <- runif(100, min=0, max=100)
y <- 2 + 3*x + rnorm(100, mean=0, sd=x*2)
m <- lm(y ~ x)

betas[i] <- coef(m)["x"]
ses[i] <- sqrt(vcov(m)["x","x"])
}

#Empirical standard error
sd(betas)

#Average calculated standard errors in samples
mean(ses)


library(AER)
data("CASchools")
head(CASchools)

#Creating Variables
CASchools$STR <- CASchools$students/CASchools$teachers
summary(CASchools$STR)

CASchools$score <- (CASchools$read + CASchools$math)/2
summary(CASchools$score)
dat <- CASchools

dat$esl <- dat$english

#lm(math ~ esl)

plot( dat$esl, dat$math)
abline(lm(dat$math ~ dat$esl), col="firebrick", lwd=2)

library(lmtest)
m <- lm(math ~ esl, data=dat)
bptest(m)

m <- lm(math ~ esl + lunch + calworks + income, data=dat)
summary(m)
bptest(m)

library(sandwich)
library(lmtest)
m <- lm(math ~ esl, data=dat)
robust.vcov <- vcovHC(m, type="HC1")
coeftest(m, vcov=robust.vcov)

gdp <- import("https://github.com/marctrussler/IIS-Data/raw/main/GDP.csv")
gdp <- gdp[213:307,]

plot(gdp$DATE, gdp$GDP)
abline(lm(gdp$GDP ~ gdp$DATE), col="firebrick", lwd=2)

gdp$residuals <- residuals(lm(gdp$GDP ~ gdp$DATE))
plot(gdp$DATE, gdp$residuals)
abline(h=0)

x <- runif(100, min=0, max=100)
y <- 2 + 3*x + rnorm(100, mean=0, sd=15)
resid <- residuals(lm(y~x))
plot(x,resid)
abline(h=0)

head(acs)
acs.reduc <- acs[acs$state.abbr %in% c("MS","NY"),]
m <- lm(median.income ~ percent.college, data=acs.reduc)
plot(acs.reduc$percent.college, acs.reduc$median.income, pch=16)
abline(m, col="black")

plot(acs.reduc$percent.college[acs.reduc$state.abbr=="MS"], acs.reduc$median.income[acs.reduc$state.abbr=="MS"], col="firebrick", pch=16)
points(acs.reduc$percent.college[acs.reduc$state.abbr=="NY"], acs.reduc$median.income[acs.reduc$state.abbr=="NY"], col="darkblue", pch=16)
abline(m, col="black")

acs.reduc$resid <- residuals(m)

plot(acs.reduc$percent.college[acs.reduc$state.abbr=="MS"], acs.reduc$resid[acs.reduc$state.abbr=="MS"], col="firebrick", pch=16)
points(acs.reduc$percent.college[acs.reduc$state.abbr=="NY"], acs.reduc$resid[acs.reduc$state.abbr=="NY"], col="darkblue", pch=16)
abline(h=0)



m <- lm(median.income ~ percent.college, data=acs.reduc)
summary(m)
library(sandwich)
cluster.vcov <- vcovCL(m, cluster= ~ state.abbr, type="HC1")
cluster.vcov2 <- vcovCL(m, cluster= ~ state.abbr)

#Compare to original var-covar matrix
vcov(m)
#New hypothesis test
coeftest(m, vcov = cluster.vcov )


#You can also put the new variance covariance matrix into stargazer:

#Extract standard errors from new vcov matrix:
cluster.se <- sqrt(diag(cluster.vcov))

#(I just have to copy and paster that list(NULL) part of it as I always forget how to do it.)
stargazer(m, type="text", se=list(NULL,cluster.se))

acs.reduc$ny <- acs.reduc$state.abbr=="NY"

m2<- lm(median.income ~ percent.college + ny, data=acs.reduc)
summary(m2)


plot(acs.reduc$percent.college[acs.reduc$state.abbr=="MS"], acs.reduc$median.income[acs.reduc$state.abbr=="MS"], col="firebrick", pch=16)
abline(a=coef(m2)["(Intercept)"], b=coef(m2)["percent.college"], col="firebrick")
points(acs.reduc$percent.college[acs.reduc$state.abbr=="NY"], acs.reduc$median.income[acs.reduc$state.abbr=="NY"], col="darkblue", pch=16)
abline(a=coef(m2)["(Intercept)"] + coef(m2)["nyTRUE"], b=coef(m2)["percent.college"], col="darkblue")

acs.reduc$resid <- residuals(m2)

plot(acs.reduc$percent.college[acs.reduc$state.abbr=="MS"], acs.reduc$resid[acs.reduc$state.abbr=="MS"], col="firebrick", pch=16)
points(acs.reduc$percent.college[acs.reduc$state.abbr=="NY"], acs.reduc$resid[acs.reduc$state.abbr=="NY"], col="darkblue", pch=16)
abline(h=0)

plot(acs$percent.college, acs$median.income, col="black", pch=16)

#Standard OLS
m <- lm(median.income ~ percent.college, data=acs)
summary(m)
#Cluster standard errors
cluster.vcov <- vcovCL(m, cluster = ~ state.abbr)

coeftest(m, cluster.vcov)


acs$state <- as.factor(acs$state.abbr)
m2 <- lm(median.income ~ percent.college + state, data=acs)
summary(m2)

#install.packages("lfe")
library(lfe)
m.fe <- felm(median.income ~ percent.college | state.abbr, data=acs)
summary(m.fe)

getfe(m.fe)


m.fe.cl <- felm(median.income ~ percent.college | state.abbr | 0 | state.abbr, data=acs)
summary(m.fe.cl)
#Much higher standard error!

library(AER)
data("CASchools")
head(CASchools)

#Creating Variables
CASchools$STR <- CASchools$students/CASchools$teachers
summary(CASchools$STR)

CASchools$score <- (CASchools$read + CASchools$math)/2
summary(CASchools$score)

#Changing names
dat <- CASchools
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 1), se = FALSE)
p
summary(lm(score ~ income, data=dat))

#YOU MUST USE THE RAW=T OPTION WHEN YOU USE THE POLY FUNCTION
#OR IT WILL NOT WORK

m <- lm(score ~ poly(income, 2, raw=T), data=dat )
summary(m)

#This is equivalent to: 

dat$income2  <- dat$income^2

m <- lm(score ~ income + income2, data=dat )
summary(m)

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 2), se = FALSE)
p

coef(m)

coef(m)[2] + 2*coef(m)[3]*20

#2.15
#and at 30
coef(m)[2] + 2*coef(m)[3]*30

#and at 50
coef(m)[2] + 2*coef(m)[3]*50

#Plotting this relationship:
library(ggplot2)
p <- ggplot(dat, aes(x=income, y=score)) + geom_point() +
  labs(x="Avereage Household Income", y = "School Test Scores") +
  geom_smooth(method = lm, formula = y ~ poly(x, 4), se = FALSE)
p

m <- lm(score ~ poly(income, 4, raw=T), data=dat )
coef(m)

#20
coef(m)[2] + 2*coef(m)[3]*20 + 3*coef(m)[4]*20^2 + 4*coef(m)[5]*20^3

#and at 30
coef(m)[2] + 2*coef(m)[3]*30 + 3*coef(m)[4]*30^2 + 4*coef(m)[5]*30^3

#and at 50
coef(m)[2] + 2*coef(m)[3]*50 + 3*coef(m)[4]*50^2 + 4*coef(m)[5]*50^3