# ---------------------------- ECONOMETRIA -------------------------------#
#         Trabajo Práctico No. 1 - Modelo lineal y estimación MCO         #
#                           Resolucion oficial                            #
# ------------------------------------------------------------------------#

# Para borrar todo lo que pueda estar abierto tenemos:
rm(list=ls()) #abreviación de rm(list = ls(all = TRUE)) 

# Directorio de la carpeta donde se encuentra la base de datos
dataDir <- "C:\\Users\\pupih\\Desktop\\Asistente full time\\2019\\Econometría Avanzada\\Trabajo Práctico 1\\Datos"
# Directorio de la carpeta donde guardaremos los resultados
resultsDir <- "C:\\Users\\pupih\\Desktop\\Asistente full time\\2019\\Econometría Avanzada\\Trabajo Práctico 1\\Resultados"


# Cargamos los paquetes que necesitamos
install.packages("stargazer")  # paquete para tener una salida estética de las regresiones
install.packages("car")  # paquete para realizar test de hipótesis lineales
library(corrplot)
library(RColorBrewer)
library(stargazer)
library(car)                    

# Cargamos la base de datos
setwd(dataDir) 
TP1data <- read.csv("TP1_MCO.csv")


# Revisamos la base de datos
View(TP1data) # para ver la base de datos como una tabla
str(TP1data) # para ver el tipo de datos guardados: character, numeric, integer, logical, complex, factor

# Transformo en factor a sex y everbfed. ¿Qué es un factor?
# Son variables en R que toman un número limitado de valores diferentes. 
# En este caso, hombre o mujer para sex y sí o no para everbfed.
# Más información en https://www.stat.berkeley.edu/~s133/factors.html 

TP1data$sex <- factor(TP1data$sex,
                    levels = c(0,1),
                    labels = c("female", "male"))

TP1data$everbfed <- factor(TP1data$everbfed,
                      levels = c(0,1),
                      labels = c("No", "yes"))
str(TP1data) # Chequeo


# PUNTO 2) -----> ESTADÍSTICA DESCRIPTIVA

summary(TP1data) 

# Gráfico de dispersión de la variables
pairs(TP1data[,-c(5,6)], col="blue") # sin las variables cualitativas!

# Matriz de correlaciones
round(cor(TP1data[,-c(5,6)]), 2) # sin las variables cualitativas!

# Gráfico de Matriz de correlaciones
corS<-cor(TP1data[,-c(5,6)])
corrplot(corS,tl.pos = "lt", tl.col = "black" )

# Puntajes por sexo y niños que fueron amamantados:
layout(matrix(1:2, nrow = 1,ncol = 2))
for (i in 5:6){
  plot(TP1data[,i], TP1data$score_ppvt, ylab="Puntaje obtenido en PPVT en 2007", main = names(TP1data[i]), col=i)
}
# Explicación del comando: quiero armar un gráfico con los puntajes según sexo y según si fueron amamantados o no.
# Con "for (i in 5:6)" armo un loop llamando a las columnas 5 y 6 del dataset: sex y everbfed.
# ¿Cómo sé que llamo al dataset correcto? Gracias a "TP1data[,i]", donde le digo que tome todas las filas
# de las columnas i (es decir, las columnas 5 y 6). Con "TP1data$score_ppvt" le digo a R que quiero ver 
# solo los datos de score para esas columnas.
# Con "ylab="Puntaje obtenido en PPVT en 2007"" indico el nombre del eje y con "main = names(TP1data[i])" el nombre 
# de cada gráfico.
# Por último, con "col=i" indico que quiero que use un color distinto para sex y para everbfed. Fíjense que si lo 
# sacan, ambos gráficos estarán en blanco y negro.

# PUNTO 3) -----> ESTANDARIZACIÓN
# Usaremos el comando scale, función genérica que centra o escala las columnas de una matriz.
TP1data$score_ppvt_std <- scale(TP1data$score_ppvt)
TP1data$bwght_std <- scale(TP1data$bwght)
TP1data$mtweight_std <- scale(TP1data$mtweight)

# También es posible hacerlo "a mano"
# TP1data$score_ppvt_std <- (TP1data$score_ppvt - mean(TP1data$score_ppvt)) / sd(TP1data$score_ppvt)
# TP1data$bwght_std <- (TP1data$bwght - mean(TP1data$bwght)) / sd(TP1data$bwght)
# TP1data$mtweight_std <- (TP1data$mtweight - mean(TP1data$mtweight)) / sd(TP1data$mtweight)


# Estadística descriptiva incluyendo estas variables)
summary(TP1data)


# PUNTO 4) -----> MODELO LINEAL

# E(score_ppvt_std|bwght_std, mtweight_std, agechild, sex, everbfed) = b0 + b1*bwght_std + b2*mtweight_std + b3*agechild + b4*sex + b5*everbfed + ui


# PUNTOS 5 Y 6) -----> ESTIMACIÓN DEL MODELO LINEAL (1)
attach(TP1data) # Para no tener que llamar siempre a la base antes de nombrar una columna,
# existe este comando que nos permite fijar la base. Si en algún momento quiero cambiar de 
# base, debo sacar esta fijación usando detach(nombre de la base).

lm.fit=lm(score_ppvt_std ~ bwght_std + mtweight_std + agechild + sex + everbfed)
# La función lm (Linear Models) sirve para correr regresiones lineales
lm.fit # Para ver los resultados

# Testeo si los coeficientes son iguales: bwght_std = mtweight_std
linearHypothesis(lm.fit, "bwght_std = mtweight_std")

# PUNTO 7) -----> ESTIMACIÓN DEL MODELO LINEAL (2)
# La función I() facilita realizar operaciones matemáticas dado que indica a R que símbolos
# tales como "+", "-", "*" y "^" son operadores aritméticos.
# En este caso, elevo al cuadrado la variable estandarizada de peso al nacer. ¿Por qué?
TP1data$bwght_std2 <- I(bwght_std^2)

lm.fit2=lm(score_ppvt_std ~ bwght_std + mtweight_std + bwght_std2 + agechild + sex + everbfed )


# Presentamos los resultados en una sola tabla 
setwd(resultsDir) # Indico el directorio donde guardaré los resultados, definido más arriba
stargazer(lm.fit, lm.fit2, type="text",
          dep.var.labels=c("score_ppvt_std")
          , out="models1.txt")

# ------------------------------------------------------------------------#