?>

library(AER)

library(stargazer)

 

data(HMDA)

head(HMDA)

 

str(HMDA)

 

HMDA$deny <- as.numeric(HMDA$deny)

str(HMDA)

 

HMDA$deny <- HMDA$deny -1

str(HMDA)

 

denymod1 <- lm(deny ~ pirat, data = HMDA)

summary(denymod1)

 

windows()

# plot the data

plot(x = HMDA$pirat,

y = HMDA$deny,

main = “Scatterplot Mortgage Application Denial and the Payment-to-Income Ratio”,

xlab = “P/I ratio”,

ylab = “Deny”,

pch = 20,

ylim = c(-0.4, 1.4),

cex.main = 0.8)

 

# add horizontal dashed lines and text

abline(h = 1, lty = 2, col = “darkred”)

abline(h = 0, lty = 2, col = “darkred”)

text(2.5, 0.9, cex = 0.8, “Mortgage denied”)

text(2.5, -0.1, cex= 0.8, “Mortgage approved”)

 

abline(denymod1,

lwd = 1.8,

col = “steelblue”)

 

windows()

plot(denymod1,1)

 

windows()

plot(denymod1,2)

 

# estimate the simple logit model

denylogit <- glm(deny ~ pirat,

family = binomial(link = “logit”),

data = HMDA)                 #glm是廣義的意思

 

summary(denylogit)

 

windows()

# plot data

plot(x = HMDA$pirat,

y = HMDA$deny,

main = “Logit Model of the Probability of Denial, Given P/I Ratio”,

xlab = “P/I ratio”,

ylab = “Deny”,

pch = 20,

ylim = c(-0.4, 1.4),

cex.main = 0.85)

 

# add horizontal dashed lines and text

abline(h = 1, lty = 2, col = “darkred”)

abline(h = 0, lty = 2, col = “darkred”)

text(2.5, 0.9, cex = 0.8, “Mortgage denied”)

text(2.5, -0.1, cex= 0.8, “Mortgage approved”)

 

# add estimated regression line

x <- seq(0, 3, 0.01)

y <- predict(denylogit, list(pirat = x), type = “response”)

 

lines(x, y, lwd = 1.5, col = “red”)

 

predictions <- predict(denylogit,

newdata = data.frame(“pirat” = c(0.5, 1.5)),

type = “response”)

 

predictions

diff(predictions)

 

str(HMDA)

colnames(HMDA)[colnames(HMDA) == “afam”] <- “black” #colnames的名稱由afam改為black

str(HMDA)

 

denylogit2 <- glm(deny ~ pirat + black,

family = binomial(link = “logit”),

data = HMDA)                 #2個解釋變數

 

summary(denylogit2)

 

predictions <- predict(denylogit2,

newdata = data.frame(“black” = c(“no”, “yes”),

“pirat” = c(0.3, 0.3)),

type = “response”)

 

predictions

 

denyprobit <- glm(deny ~ pirat,

family = binomial(link = “probit”),

data = HMDA)

 

summary(denyprobit)

 

windows()

# plot data

plot(x = HMDA$pirat,

y = HMDA$deny,

main = “Probit Model of the Probability of Denial, Given P/I Ratio”,

xlab = “P/I ratio”,

ylab = “Deny”,

pch = 20,

ylim = c(-0.4, 1.4),

cex.main = 0.85)

 

# add horizontal dashed lines and text

abline(h = 1, lty = 2, col = “darkred”)

abline(h = 0, lty = 2, col = “darkred”)

text(2.5, 0.9, cex = 0.8, “Mortgage denied”)

text(2.5, -0.1, cex= 0.8, “Mortgage approved”)

 

# add estimated regression line

x <- seq(0, 3, 0.01)

y <- predict(denyprobit, list(pirat = x), type = “response”)

 

lines(x, y, lwd = 1.5, col = “blue”)

 

denyprobit2 <- glm(deny ~ pirat + black,

family = binomial(link = “probit”),

data = HMDA)

 

summary(denyprobit2)

 

library(stargazer)

stargazer(denylogit,denylogit2,denyprobit,denyprobit2, type = “html”, out = “binary data analysis.html” ,title = “My Regression models”)

#修改這段試試

stargazer(denylogit,denylogit2,denyprobit,denyprobit2, type = “text“, out = “binary data analysis.text ” ,title = “My Regression models”)

 

windows()

# plot data

plot(x = HMDA$pirat,

y = HMDA$deny,

main = “Probit and Logit Models Model of the Probability of Denial, Given P/I Ratio”,

xlab = “P/I ratio”,

ylab = “Deny”,

pch = 20,

ylim = c(-0.4, 1.4),

cex.main = 0.9)

 

# add horizontal dashed lines and text

abline(h = 1, lty = 2, col = “darkred”)

abline(h = 0, lty = 2, col = “darkred”)

text(2.5, 0.9, cex = 0.8, “Mortgage denied”)

text(2.5, -0.1, cex= 0.8, “Mortgage approved”)

 

# add estimated regression line of Probit and Logit models

x <- seq(0, 3, 0.01)

y_probit <- predict(denyprobit, list(pirat = x), type = “response”)

y_logit <- predict(denylogit, list(pirat = x), type = “response”)

 

lines(x, y_probit, lwd = 1.5, col = “blue”)

#lines(x, y_logit, lwd = 1.5, col = “red”, lty = 2)

lines(x, y_logit, lwd = 1.5, col = “red”)

 

# add a legend

legend(“topleft”,

horiz = TRUE,

legend = c(“Probit”, “Logit”),

col = c(“blue”, “red”),

lty = c(1, 2))

 

221020 計量經濟學_作業 3
?>