Paso 1

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,100)

n = 30
x=runif(30,1,n)  ## Valores iniciales arbitrarios para X
y = 5+(3*x) + rnorm(n = 30,mean = 0,sd = 10 ) ##Valores de Y de acuerdo al 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   9.388830 17.916102
## 2   9.905156 31.203731
## 3   9.928494 45.243480
## 4  12.548330 36.854158
## 5   2.982914 24.703429
## 6   5.567430 13.797628
## 7  15.076827 63.220206
## 8  10.444872 49.886358
## 9  28.020525 74.177945
## 10  2.797044  6.424065

Paso 2

Se realiza el ejercicio 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.1) # rango de valores al rededor de beta0 = 5
beta1_est=seq(1,5,0.1)  # 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])
L=array(NA,dim(betas)[1])
plot(x,y, ylim = c(0,100))

varianza = 100
sigma = sqrt(varianza)

for(i in 1:dim(betas)[1]){

  #y_est=betas$beta0_est[i] + (betas$beta1_est[i]*x)

  
  L[i] = -n*log(2*pi)/2 - 
          n*log(varianza)/2 - 
          (y - betas$beta0_est[i] - betas$beta1_est[i]*x)^2/(2*varianza)
            

  #SCE[i]=(y_est-y)^2
}
resultados=data.frame(betas,L)

Paso 3

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

library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
plot_ly(x=resultados$beta0_est,
        y=resultados$beta1_est,
        z=resultados$L,
        size=.5)
## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Estimación MCO

modelo1=lm(y~x)
summary(modelo1)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.790  -7.742  -1.894   6.195  17.207 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.5368     3.7072   1.494    0.146    
## x             2.8938     0.2381  12.155  1.1e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.541 on 28 degrees of freedom
## Multiple R-squared:  0.8407, Adjusted R-squared:  0.835 
## F-statistic: 147.7 on 1 and 28 DF,  p-value: 1.1e-12

Conclusión: Al utilizar el método de máxima verosimilitud logramos ver que lo que se busca es maximizar la función a diferencia de los mínimos cuadrados que busca minimizar la suma de cuadrados del error, por lo que las gráficas son diferentes; mínimos cuadrados muestra una U y máxima verosimilitud muestra una U al revés. Por otro lado, cabe resaltar que al usar los dos métodos se llega a un resultado similar en cuanto a los valores del intercepto y la pendiente.