Estimación por \(MCO\)

Para la simulación vamos a generar un conjunto de datos basandonos en un modelo con parametros previamente definidos: Yi=5+(3Xi)+εi, en donde podemos notar que β0=5, β1=3 y εi ~ Normal(0,5).

Paso 1

#Importar la base de datos 
require(dplyr)
library(pacman)
p_load(data.table, fixest, lattice, magrittr, ggplot2, kableExtra)

data_sm <- read.csv("~/OneDrive - PUJ Cali/Simulación MCO.txt", sep=";")
data_sm = select(data_sm, -X)

x=runif(30,1,30) #Valores iniciales arbitrarios para X
y = 5+(3*x) + rnorm(n =30, mean = 0, sd = 10) #Valores de y de acuerdo con el modelo

plot(x,y, xlim = c(0,30), ylim = c(0,100), las=1, pch=19, col="#b0394a")

datos=data.frame(x,y)
head(datos,10)
##            x        y
## 1   6.103840 30.59239
## 2  20.247299 55.39141
## 3   3.867035 37.97723
## 4  26.869554 84.31660
## 5  10.061096 46.21643
## 6  23.491288 70.61461
## 7  20.522702 74.98180
## 8  18.501761 57.28081
## 9   5.377932 21.75875
## 10 15.641313 60.15640
b0=4  # intercepto estimado
b1=2.5 # pendiente estimada

x1=1:30
y1 = b0 + b1 * x1

plot(x,y, xlim = c(0,30), ylim = c(0,100), las=1,pch=19 ,col="red")
lines(x1, y1 , col = 2, type = "l")
grid()

Paso 3

Se realiza el Paso 2 para otros posibles valores de los parametros y se guarda el valor de SCE (Suma de Cuadrados del Error)

beta0_est= seq(3,7,0.5) # rango de valores al rededor de beta0 = 5
beta1_est= seq(1,5,0.3)  # rango de valores al rededor de beta1 = 3 

betas=expand.grid(beta0_est,beta1_est)
names(betas)=c("beta0_est","beta1_est")
SCE=array(NA,dim(betas)[1])
plot(x,y, ylim = c(0,100))

for(i in 1:dim(betas)[1]){
y_est=betas$beta0_est[i] + (betas$beta1_est[i]*x)
lines(x,y_est, type = "l")
SCE[i]=(y_est-y)^2
}

x=data.frame(betas,SCE)

Paso 4

Graficar e interpretar la relación entre los posibles Parametros (betas) y la SCE.

library(plotly)

resultados=data.frame(betas,SCE)
# 
# plot_ly(x=resultados$beta0_est,
#         y=resultados$beta1_est,
#         z=resultados$SCE,
#         size=.5)

En la Grafica se puede observar que existe un punto mínimo en la SCE para una combinación de los valores β0ˆ y β1ˆ, esto en general lo que nos indica es que el método de MCO busca justamente cuales son esos dos valores tal que nos garantizan una mínima suma de cuadrados del error.

Finalmente se realiza la estimación por MCO utilizando la función lm de R, sabiendo que estos valores corresponden a la combinacion de valores de los estiadores que minimizan SCE

Estimación por \(MLE\)

library(pacman)
p_load(data.table, fixest, lattice, magrittr, ggplot2, kableExtra)

Paso 1

Se procede a generar los datos con los cuales se pretende trabajar la estación por maxima verisimilitud

N = 30        # Sample size
beta_0 = 5    # Beta cero
beta_1 = 3    # Beta uno
sigma_2 = 5  # sigma^2 (Distribution of U)
DT <- data.table(X = runif(N,1,N),
                 U = rnorm(N, mean = 0, sd = 10)) %>%
    .[, Y := beta_0 + beta_1*X + U]

# Y = Modelo
# X = Variable regresora o predictora
# U = Error

Paso 2

Se procede a identificar las derivadas de la función de maxima verisimitud

# Función de verosimilitud
log_like <- function(theta, Y, X){
  X <- as.matrix(X); Y <- as.matrix(Y);
  N       <- nrow(X)
  beta_0    <- theta[1]
  beta_1 <- theta[2]
  sigma_2 <- theta[3]
  e <- Y - beta_0 - beta_1*X
  loglik  <- -.5*N*log(2*pi)-.5*N*log(sigma_2) - ( (t(e) %*% e)/ (2*sigma_2) )
  return(-loglik)
}


# Graficar la función de verosimilitud

log_like_graph <- function(beta_0, beta_1, sigma_2){
  X <- as.matrix(DT$X); Y <- as.matrix(DT$Y);
  N       <- nrow(X)
  e       <- Y - beta_0 - beta_1*X
  loglik  <- -.5*N*log(2*pi)-.5*N*log(sigma_2) - ( (t(e) %*% e)/ (2*sigma_2) )
  return(loglik)
}
log_like_graph <- Vectorize(log_like_graph)

#Establecer cuadrícula de valores beta y sigma2

beta_0_vals <- seq(-30,30, by =1)
beta_1_vals <- seq(-30,30, by =1)
sigma2_vals <- seq(-30,30, by =1)
# log_vals <- outer(beta_0_vals, beta_1_vals, sigma2_vals, log_like_graph)
# 
# persp(beta_0_vals, beta_1_vals, sigma2_vals, log_vals, theta = 7, phi =8 , r= 30)