Pruebas de Normalidad para los residuos

Econometria UES

Importación de Datos

library(readr)
ejemplo_regresion_1_ <- read_csv("C:/Users/Genovez/Downloads/ejemplo_regresion (1).csv")
head(ejemplo_regresion_1_)
## # A tibble: 6 × 3
##      X1    X2     Y
##   <dbl> <dbl> <dbl>
## 1  3.92  7298  0.75
## 2  3.61  6855  0.71
## 3  3.32  6636  0.66
## 4  3.07  6506  0.61
## 5  3.06  6450  0.7 
## 6  3.11  6402  0.72

Modelo Estimado

library(stargazer)
modelo_estimado<-lm("Y~.",data=ejemplo_regresion_1_)
stargazer(modelo_estimado,title = "Modelo para Ejemplo", type = "html")
Modelo para Ejemplo
Dependent variable:
NA
X1 0.237***
(0.056)
X2 -0.0002***
(0.00003)
Constant 1.564***
(0.079)
Observations 25
R2 0.865
Adjusted R2 0.853
Residual Std. Error 0.053 (df = 22)
F Statistic 70.661*** (df = 2; 22)
Note: p<0.1; p<0.05; p<0.01

Ajuste de los residuos a la Distribución Normal

Verificando el ajuste de los residuos a la distribución normal, se usará la librería fitdistrplus

library(fitdistrplus)
fit_normal<-fitdist(data = modelo_estimado$residuals,distr = "norm")
plot(fit_normal) 

summary(fit_normal)
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##          estimate  Std. Error
## mean 1.804926e-18 0.010000382
## sd   5.000191e-02 0.007058615
## Loglikelihood:  39.41889   AIC:  -74.83778   BIC:  -72.40002 
## Correlation matrix:
##               mean            sd
## mean  1.000000e+00 -9.796171e-14
## sd   -9.796171e-14  1.000000e+00

1. Prueba de Normalidad de Jarque Bera

Usando tseries

library(tseries)
salida_JB<-jarque.bera.test(modelo_estimado$residuals)
salida_JB
## 
##  Jarque Bera Test
## 
## data:  modelo_estimado$residuals
## X-squared = 0.93032, df = 2, p-value = 0.628
library(fastGraph)
alpha_sig<-0.05
JB<-salida_JB$statistic
gl<-salida_JB$parameter
VC<-qchisq(1-alpha_sig,gl,lower.tail = TRUE)
shadeDist(JB,ddist = "dchisq",
          parm1 = gl,
          lower.tail = FALSE,xmin=0,
          sub=paste("VC:",round(VC,2)," ","JB:",round(JB,2)))

Usando notmtest (descatalogada)

library(normtest)
jb.norm.test(modelo_estimado$residuals)
## 
##  Jarque-Bera test for normality
## 
## data:  modelo_estimado$residuals
## JB = 0.93032, p-value = 0.4925

2. Prueba de Kolmogorov Smirnov -Lilliefors

Cálculo Manual

library(dplyr)  # Carga la librería dplyr para manipulación de datos
library(gt)  # Carga la librería gt para crear tablas de datos
library(gtExtras)  # Carga la librería gtExtras para agregar funcionalidades a las tablas creadas con gt
residuos<-modelo_estimado$residuals  # Crea un vector con los residuos del modelo estimado
residuos %>%  # Utiliza el operador %>% para encadenar las operaciones siguientes al vector residuos
  as_tibble() %>%  # Convierte el vector residuos en una tibble (tabla) de una columna
  mutate(posicion=row_number()) %>%  # Agrega una columna llamada "posicion" con el número de fila
  arrange(value) %>%  # Ordena la tabla por los valores de residuos en orden ascendente
  mutate(dist1=row_number()/n()) %>%  # Agrega una columna "dist1" con los percentiles según su posición en la tabla (usando la función row_number() y n() para obtener el número de filas)
  mutate(dist2=(row_number()-1)/n()) %>%  # Agrega una columna "dist2" con los percentiles según su posición en la tabla, pero ajustando en una unidad para evitar problemas con los extremos de la distribución
  mutate(zi=as.vector(scale(value,center=TRUE))) %>%  # Agrega una columna "zi" con los valores de residuos escalados para tener media cero y varianza uno
  mutate(pi=pnorm(zi,lower.tail = TRUE)) %>%  # Agrega una columna "pi" con los valores de la función de distribución acumulada (CDF) de una distribución normal estándar evaluada en los valores de zi
  mutate(dif1=abs(dist1-pi)) %>%  # Agrega una columna "dif1" con las diferencias absolutas entre los percentiles según la posición y los valores de pi
  mutate(dif2=abs(dist2-pi)) %>%  # Agrega una columna "dif2" con las diferencias absolutas entre los percentiles ajustados según la posición y los valores de pi
  rename(residuales=value) -> tabla_KS  # Renombra la columna "value" como "residuales" y asigna la tabla resultante a la variable tabla_KS


#Formato
 tabla_KS %>%  # Utiliza el operador %>% para encadenar las operaciones siguientes a la tabla tabla_KS
  gt() %>%  # Crea una tabla con la función gt()
  tab_header("Tabla para calcular el Estadistico KS") %>%  # Agrega un encabezado a la tabla
  tab_source_note(source_note = "Fuente: Elaboración propia") %>%  # Agrega una nota de fuente a la tabla
  tab_style(  # Cambia el estilo de algunas celdas de la tabla
    style = list(
      cell_fill(color = "#A569BD"),  # Cambia el color de fondo de las celdas a un tono de morado
      cell_text(style = "italic")  # Cambia el estilo de texto de las celdas a itálico
      ),
    locations = cells_body(  # Aplica el estilo a las celdas del cuerpo de la tabla que cumplan las siguientes condiciones:
      columns = dif1,  # Que pertenezcan a la columna "dif1"
      rows = dif1==max(dif1)  # Que pertenezcan a la fila donde el valor de "dif1" es máximo
    )) %>%
   tab_style(  # Cambia el estilo de algunas celdas de la tabla
    style = list(
      cell_fill(color = "#3498DB"),  # Cambia el color de fondo de las celdas a un tono de azul
      cell_text(style = "italic")  # Cambia el estilo de texto de las celdas a itálico
      ),
    locations = cells_body(  # Aplica el estilo a las celdas del cuerpo de la tabla que cumplan las siguientes condiciones:
      columns = dif2,  # Que pertenezcan a la columna "dif2"
      rows = dif2==max(dif2)  # Que pertenezcan a la fila donde el valor de "dif2" es máximo
    ))
Tabla para calcular el Estadistico KS
residuales posicion dist1 dist2 zi pi dif1 dif2
-0.085089833 19 0.04 0.00 -1.66734963 0.04772245 0.007722452 0.047722452
-0.069637086 22 0.08 0.04 -1.36455045 0.08619719 0.006197193 0.046197193
-0.062183196 4 0.12 0.08 -1.21849022 0.11151887 0.008481131 0.031518869
-0.057379932 25 0.16 0.12 -1.12436944 0.13042817 0.029571834 0.010428166
-0.057277771 12 0.20 0.16 -1.12236757 0.13085309 0.069146907 0.029146907
-0.048434733 11 0.24 0.20 -0.94908676 0.17128824 0.068711757 0.028711757
-0.039102258 3 0.28 0.24 -0.76621535 0.22177409 0.058225910 0.018225910
-0.033971640 16 0.32 0.28 -0.66568001 0.25280783 0.067192174 0.027192174
-0.024886579 17 0.36 0.32 -0.48765671 0.31289651 0.047103492 0.007103492
-0.022166535 13 0.40 0.36 -0.43435700 0.33201461 0.067985390 0.027985390
-0.018597878 8 0.44 0.40 -0.36442855 0.35776901 0.082230992 0.042230992
-0.010121525 21 0.48 0.44 -0.19833299 0.42139227 0.058607729 0.018607729
-0.003341163 2 0.52 0.48 -0.06547065 0.47389964 0.046100362 0.006100362
0.012424659 6 0.56 0.52 0.24346329 0.59617674 0.036176739 0.076176739
0.016240338 5 0.60 0.56 0.31823217 0.62484558 0.024845584 0.064845584
0.017281530 20 0.64 0.60 0.33863450 0.63255746 0.007442545 0.032557455
0.026439478 10 0.68 0.64 0.51808602 0.69780087 0.017800874 0.057800874
0.026966239 18 0.72 0.68 0.52840800 0.70139191 0.018608095 0.021391905
0.030236216 7 0.76 0.72 0.59248375 0.72323665 0.036763350 0.003236650
0.035565142 15 0.80 0.76 0.69690497 0.75706887 0.042931135 0.002931135
0.040753758 14 0.84 0.80 0.79857675 0.78773206 0.052267936 0.012267936
0.072518915 23 0.88 0.84 1.42102038 0.92234458 0.042344582 0.082344582
0.073469743 1 0.92 0.88 1.43965200 0.92501706 0.005017061 0.045017061
0.074601871 24 0.96 0.92 1.46183625 0.92810696 0.031893044 0.008106956
0.105692240 9 1.00 0.96 2.07105726 0.98082328 0.019176722 0.020823278
Fuente: Elaboración propia

Cálculo del estadístico D:

D<-max(max(tabla_KS$dif1),max(tabla_KS$dif2))
print(D)
## [1] 0.08234458

Usando nortest

library(nortest)
prueba_KS<-lillie.test(modelo_estimado$residuals)
prueba_KS
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_estimado$residuals
## D = 0.082345, p-value = 0.9328
p.value<-prueba_KS$p.value

En este caso dado que 0.9327912 >0.05 No se rechaza la Hipótesis Nula: ϵ∼N(0,σ2) , por lo que los residuos siguen una distribución normal.

3. Prueba de Shapiro - Wilk

Cálculo Manual

library(dplyr)
library(gt)
residuos<-modelo_estimado$residuals
residuos %>%  
  as_tibble() %>%
  rename(residuales=value) %>%
  arrange(residuales) %>%
  mutate(pi=(row_number()-0.375)/(n()+0.25)) %>%
  mutate(mi=qnorm(pi,lower.tail = TRUE)) %>% 
  mutate(ai=0)->tabla_SW

m<-sum(tabla_SW$mi^2)
n<-nrow(ejemplo_regresion_1_)
theta<-1/sqrt(n)
tabla_SW$ai[n]<- -2.706056*theta^5+4.434685*theta^4-2.071190*theta^3-0.147981*theta^2+0.2211570*theta+tabla_SW$mi[n]/sqrt(m)
tabla_SW$ai[n-1]<- -3.582633*theta^5+5.682633*theta^4-1.752461*theta^3-0.293762*theta^2+0.042981*theta+tabla_SW$mi[n-1]/sqrt(m)
tabla_SW$ai[1]<- -tabla_SW$ai[n]
tabla_SW$ai[2]<- -tabla_SW$ai[n-1]
omega<-(m-2*tabla_SW$mi[n]^2-2*tabla_SW$mi[n-1]^2)/(1-2*tabla_SW$ai[n]^2-2*tabla_SW$ai[n-1]^2)
tabla_SW$ai[3:(n-2)]<-tabla_SW$mi[3:(n-2)]/sqrt(omega)

tabla_SW %>% 
  mutate(ai_ui=ai*residuales,ui2=residuales^2) ->tabla_SW

tabla_SW %>%
  gt() %>% tab_header("Tabla para calcular el Estadistico W") %>%  # Agrega un encabezado a la tabla
  tab_source_note(source_note = "Fuente: Elaboración propia")
Tabla para calcular el Estadistico W
residuales pi mi ai ai_ui ui2
-0.085089833 0.02475248 -1.96421684 -0.44179982 0.0375926725 7.240280e-03
-0.069637086 0.06435644 -1.51919718 -0.31084122 0.0216460772 4.849324e-03
-0.062183196 0.10396040 -1.25930333 -0.25447584 0.0158241211 3.866750e-03
-0.057379932 0.14356436 -1.06444156 -0.21509882 0.0123423556 3.292457e-03
-0.057277771 0.18316832 -0.90335666 -0.18254732 0.0104559034 3.280743e-03
-0.048434733 0.22277228 -0.76286392 -0.15415701 0.0074665537 2.345923e-03
-0.039102258 0.26237624 -0.63603674 -0.12852819 0.0050257427 1.528987e-03
-0.033971640 0.30198020 -0.51871372 -0.10481995 0.0035609055 1.154072e-03
-0.024886579 0.34158416 -0.40814351 -0.08247629 0.0020525526 6.193418e-04
-0.022166535 0.38118812 -0.30236184 -0.06110028 0.0013543815 4.913553e-04
-0.018597878 0.42079208 -0.19986756 -0.04038857 0.0007511418 3.458811e-04
-0.010121525 0.46039604 -0.09943603 -0.02009370 0.0002033789 1.024453e-04
-0.003341163 0.50000000 0.00000000 0.00000000 0.0000000000 1.116337e-05
0.012424659 0.53960396 0.09943603 0.02009370 0.0002496574 1.543722e-04
0.016240338 0.57920792 0.19986756 0.04038857 0.0006559241 2.637486e-04
0.017281530 0.61881188 0.30236184 0.06110028 0.0010559063 2.986513e-04
0.026439478 0.65841584 0.40814351 0.08247629 0.0021806299 6.990460e-04
0.026966239 0.69801980 0.51871372 0.10481995 0.0028265997 7.271780e-04
0.030236216 0.73762376 0.63603674 0.12852819 0.0038862062 9.142287e-04
0.035565142 0.77722772 0.76286392 0.15415701 0.0054826161 1.264879e-03
0.040753758 0.81683168 0.90335666 0.18254732 0.0074394891 1.660869e-03
0.072518915 0.85643564 1.06444156 0.21509882 0.0155987331 5.258993e-03
0.073469743 0.89603960 1.25930333 0.25447584 0.0186962745 5.397803e-03
0.074601871 0.93564356 1.51919718 0.31084122 0.0231893370 5.565439e-03
0.105692240 0.97524752 1.96421684 0.44179982 0.0466948120 1.117085e-02
Fuente: Elaboración propia

Cálculo del Estadistico W

W<-(sum(tabla_SW$ai_ui)^2)/sum(tabla_SW$ui2)
print(W)
## [1] 0.9700088

Cálculo del Wny su p value

mu<-0.0038915*log(n)^3-0.083751*log(n)^2-0.31082*log(n)-1.5861
sigma<-exp(0.0030302*log(n)^2-0.082676*log(n)-0.4803)
Wn<-(log(1-W)-mu)/sigma
print(Wn)
## [1] -0.3726406
p.value<-pnorm(Wn,lower.tail = FALSE)
print(p.value)
## [1] 0.645292
library(fastGraph)
shadeDist(Wn,ddist = "dnorm",lower.tail = FALSE)

En este caso dado que 0.645292 >0.05 No se rechaza la Hipótesis Nula: ϵ∼N(0,σ2) , por lo que los residuos siguen una distribución normal.

Usando la librería stats (precargada al iniciar R)

salida_SW<-shapiro.test(modelo_estimado$residuals)
print(salida_SW)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_estimado$residuals
## W = 0.97001, p-value = 0.6453

Mismos resultados que en el cálculo manual.

Importante, a partir de esta salida se puede calcular el Wn si se llegará a necesitar:

Wn_salida<-qnorm(salida_SW$p.value,lower.tail = FALSE)
print(Wn_salida)
## [1] -0.3726406