### Monitoria de EAE 0325 - Econometria II
### Prof. Raphael Corbi/ Monitor Alan Leal

if(!require(mfx)){install.packages("mfx");require(mfx)}
### Modelos de Escolhas Qualitativas

###### A) Estimando beta para o PROBIT e o LOGIT
# MLE do PROBIT e do LOGIT
probit_my <- function(beta){
  p <- sum(y*log(pnorm(X%*%beta))+(1-y)*log(1-pnorm(X%*%beta)))
  return(-p) # Optim minimiza a função, por isso o -p como objeto retornado
  # optim(a,f): a é o chute inicial e f é a função a ser otimizada. 
}

logit_my <- function(beta){
  p <- sum(y*log(plogis(X%*%beta))+(1-y)*log(1-plogis(X%*%beta)))
  return(-p)
}
#Lendo a base de dados (item do 2 do primeiro exame prático)
bangladesh_impact <- read.csv("C:\\Users\\Alanm\\OneDrive\\PAE\\2022I\\Exames Práticos\\Primeiro Exame Prático\\Datasets\\Bangladesh_microcredit_impact.csv")
# Criando as matrizes y e X:
y <- bangladesh_impact$treat
X <- cbind(1,bangladesh_impact$sexhead,bangladesh_impact$educhead,bangladesh_impact$famsize)
# primeiramente, vamos clacular o PROBIT:
probit_coefs <- optim(rep(0,4),probit_my)
# Comparar nossos resultados com a função glm:
probit_cp <-  glm(treat~sexhead+educhead+famsize,data = bangladesh_impact, family=binomial(link="probit"))
cbind(probit_coefs$par,probit_cp$coefficients)

# Agora, para o LOGIT:
logit_coefs <- optim(rep(0,4),logit_my)
# Comparar nossos resultados com a função glm:
logit_cp <-  glm(treat~sexhead+educhead+famsize,data = bangladesh_impact, family=binomial(link="logit"))
cbind(logit_coefs$par,logit_cp$coefficients)

###### B) Calculando os efeitos parciais do PROBIT
# O cálculo dos efeitos parciais difere entre variáveis contínuas e discretas
# PEA (PROBIT)
probit_medias <- X %>% colMeans() %>% as.matrix() %>% t()
b <- probit_coefs$par %>% as.vector()
k <- length(b)
N <- dim(X)[1]
PEA <- matrix(0,ncol=1,nrow=k)
# Temos duas variáveis contínuas: educhead e famsize (var 3 e 4):
for (i in 3:4){
  PEA[i] <- dnorm(probit_medias%*%b)*b[i]
}
# Temos também uma variável discreta no modelo (var 2: sexhead):
media0 <- probit_medias
media0[2] <- 0
media1 <- probit_medias
media1[2] <- 1
PEA[2] <- pnorm(media1%*%b)-pnorm(media0%*%b)

pea_mfx <- probitmfx(y~sexhead+educhead+famsize,data=bangladesh_impact,
                     atmean=TRUE) # Efeitos marginais na média
cbind(pea_mfx$mfxest[,1],PEA[2:4,])

# APE (PROBIT):
APE <- matrix(0,ncol=1,nrow=k)
e_gb <- matrix(0,nrow=N,ncol=1)
# Vou calcular G(xbeta) para toda observação:
for (i in 1:N){
  e_gb[i] <- dnorm(X[i,]%*%b)
}
# No caso de variáveis contínuas:
e_gb <- mean(e_gb)
for (i in c(3,4)){
  APE[i] <- e_gb*b[i]
}
# No caso da variável discreta, temos:
e_gb_d <- matrix(0,ncol=1,nrow=N)
beta0 <- b[-2]
for (j in 1:N){
  e_gb_d[j] <- pnorm(X[j,-2]%*%beta0+b[2])-pnorm(X[j,-2]%*%beta0)
}
APE[2] <- mean(e_gb_d)
# Comparando com o comando próprio do R:
ape_mfx <- probitmfx(y~sexhead+educhead+famsize,data=bangladesh_impact, 
                     atmean=FALSE)
cbind(ape_mfx$mfxest[,1],APE[2:4,])

###### C) Comparando PROBIT, LOGIT e MPL
lpm <- lm(y~sexhead+educhead+famsize,data=bangladesh_impact)
cbind(probit_coefs$par,logit_coefs$par,lpm$coefficients)