# AER: Applied Econometrics with R
#
# Functions, data sets, examples, demos, and vignettes for the book 
# Christian Kleiber and Achim Zeileis (2008), Applied Econometrics with R, 
# Springer-Verlag, New York. ISBN 978-0-387-77316-2. 
#
# https://cran.r-project.org/web/packages/AER/vignettes/AER.pdf
# https://www.rdocumentation.org/packages/AER/versions/1.2-9/topics/HMDA

library(AER)

data(HMDA)

# Excluding the 4 largst pirat (pirat>1.0)
HMDA = HMDA[HMDA[,2]<1,]

attach(HMDA)

# Response variable
#deny: Was the mortgage denied?

# Dependent variables
#pirat: Payments to income ratio.
#afam:   Is the individual African-American?


boxplot(pirat~deny,outline=FALSE)

table(afam,deny)


n = nrow(HMDA)
y = rep(0,n)
y[deny=="yes"]=1
x = rep(0,n)
x[afam=="yes"]=1

plot(density(pirat[(x==0)&(y==1)]),xlim=c(0,1),ylim=c(0,7),xlab="Payments to income ratio",main="")
lines(density(pirat[(x==1)&(y==1)]),col=2)
lines(density(pirat[(x==0)&(y==0)]),col=1,lty=2)
lines(density(pirat[(x==1)&(y==0)]),col=2,lty=2)
legend("topright",legend=c("Not black, denied","Black, denied","Not black, not denied","Black, not denied"),col=c(1,2,1,2),lty=c(1,1,2,2),bty="n")


fit.probit <- glm(deny ~ pirat + afam, family = binomial(link = "probit"), data = HMDA)
fit.logit  <- glm(deny ~ pirat + afam, family = binomial(link = "logit"), data = HMDA)

summary(fit.probit)
summary(fit.logit)


N      = 1000
pirat1 = seq(0,1,length=N)
fit1   = fit.probit$coef[1]+fit.probit$coef[2]*pirat1
fit2   = fit.probit$coef[1]+fit.probit$coef[2]*pirat1+fit.probit$coef[3]
fit3   = fit.logit$coef[1]+fit.logit$coef[2]*pirat1
fit4   = fit.logit$coef[1]+fit.logit$coef[2]*pirat1+fit.logit$coef[3]
fit1   = pnorm(fit1)
fit2   = pnorm(fit2)
fit3   = 1/(1+exp(-fit3))
fit4   = 1/(1+exp(-fit4))

cutoff = c(
-fit.probit$coef[1]/fit.probit$coef[2],
-(fit.probit$coef[1]+fit.probit$coef[3])/fit.probit$coef[2],
-fit.logit$coef[1]/fit.logit$coef[2],
-(fit.logit$coef[1]+fit.logit$coef[3])/fit.logit$coef[2])

cutoff = round(cutoff,3)

plot(pirat1,fit1,col=x+1,ylim=c(-0.1,1.1),xlim=c(0,1.5),type="l",xlab="Payments to income ratio",
     ylab="Probability mortgage was denied")
abline(h=0.5,col=grey(0.5),lty=3)
lines(pirat1,fit2,col=2)
lines(pirat1,fit3,col=1,lty=2)
lines(pirat1,fit4,col=2,lty=2)

abline(v=cutoff[1],col=grey(0.5),lty=3)
abline(v=cutoff[2],col=grey(0.5),lty=3)
abline(v=cutoff[3],col=grey(0.5),lty=3)
abline(v=cutoff[4],col=grey(0.5),lty=3)

legend("topright",legend=c(paste("Probit link (AIC=",round(summary(fit.probit)$aic,1),")",sep=""),
paste("Not black - ",cutoff[1],sep=""),
paste("Black - ",cutoff[2],sep="")),col=0:2,lty=1,bty="n")

legend("bottomright",legend=c(paste("Logit link (AIC=",round(summary(fit.logit)$aic,1),")",sep=""),
paste("Not black - ",cutoff[3],sep=""),
paste("Black - ",cutoff[4],sep="")),col=0:2,lty=2,bty="n")

points(pirat[x==0],y[x==0]+rnorm(sum(x==0),0,0.01),col=1)
points(pirat[x==1],y[x==1]+rnorm(sum(x==1),0.05,0.01),col=2)

title(paste("Black = ",round(100*mean(x==1),1),"%",sep=""))