library(wooldridge)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tseries)    
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(fastGraph) 
library(knitr)
data(hprice1)
head(force(hprice1),n=5)
##   price assess bdrms lotsize sqrft colonial   lprice  lassess llotsize   lsqrft
## 1   300  349.1     4    6126  2438        1 5.703783 5.855359 8.720297 7.798934
## 2   370  351.5     3    9903  2076        1 5.913503 5.862210 9.200593 7.638198
## 3   191  217.7     3    5200  1374        0 5.252274 5.383118 8.556414 7.225482
## 4   195  231.8     3    4600  1448        1 5.273000 5.445875 8.433811 7.277938
## 5   373  319.1     4    6095  2514        1 5.921578 5.765504 8.715224 7.829630

1. Estimación del Modelo.

modelo <- lm(price ~ lotsize + sqrft + bdrms, data = hprice1)
residuos <- residuals(modelo)
n <- length(residuos)
df_res <- data.frame(res = residuos)

2. Supuestos de Normalidad:

Prueba de Jarque-Bera

df_jb_manual <- df_res %>%
  mutate(
    media = mean(res),
    m2 = (res - media)^2,
    m3 = (res - media)^3,
    m4 = (res - media)^4
  )
jb_calc <- df_jb_manual %>%
  summarise(
    S = (sum(m3) / n) / ((sum(m2) / n)^(3/2)), # Asimetría
    K = (sum(m4) / n) / ((sum(m2) / n)^2),     # Curtosis
    JB_stat = (n / 6) * (S^2 + ((K - 3)^2) / 4),
    p_valor = pchisq(JB_stat, df = 2, lower.tail = FALSE)
  )

print(kable(jb_calc, caption = "Cálculo Manual - Estadístico Jarque-Bera"))
## 
## 
## Table: Cálculo Manual - Estadístico Jarque-Bera
## 
## |         S|        K|  JB_stat| p_valor|
## |---------:|--------:|--------:|-------:|
## | 0.9606833| 5.260844| 32.27791|   1e-07|

Prueba de Kolmogorov-Smirnov

df_ks_manual <- df_res %>%
  mutate(
    media = mean(res),
    desv = sd(res)
  ) %>%
  arrange(res) %>% # 1. Ordenar los residuos
  mutate(
    i = row_number(),
    # 2. Estandarizar
    z = (res - media) / desv, 
    # 3. Función de Distribución Acumulada (FDA) Teórica Normal
    F_teorica = pnorm(z), 
    # 4. FDA Empírica
    F_emp_sup = i / n,
    F_emp_inf = (i - 1) / n,
    # 5. Calcular las diferencias absolutas máximas
    D_plus = abs(F_emp_sup - F_teorica),
    D_minus = abs(F_emp_inf - F_teorica),
    D_max_local = pmax(D_plus, D_minus)
  )
ks_calc <- data.frame(Estadistico_D = max(df_ks_manual$D_max_local))

print(kable(head(df_ks_manual, 5), caption = "Primeras 5 obs. - Proceso Manual KS"))
## 
## 
## Table: Primeras 5 obs. - Proceso Manual KS
## 
## |   |        res| media|     desv|  i|         z| F_teorica| F_emp_sup| F_emp_inf|    D_plus|   D_minus| D_max_local|
## |:--|----------:|-----:|--------:|--:|---------:|---------:|---------:|---------:|---------:|---------:|-----------:|
## |81 | -120.02645|     0| 58.79282|  1| -2.041516| 0.0205998| 0.0113636| 0.0000000| 0.0092362| 0.0205998|   0.0205998|
## |77 | -115.50870|     0| 58.79282|  2| -1.964674| 0.0247260| 0.0227273| 0.0113636| 0.0019987| 0.0133624|   0.0133624|
## |24 | -107.08089|     0| 58.79282|  3| -1.821326| 0.0342787| 0.0340909| 0.0227273| 0.0001877| 0.0115514|   0.0115514|
## |48 |  -91.24398|     0| 58.79282|  4| -1.551958| 0.0603361| 0.0454545| 0.0340909| 0.0148816| 0.0262452|   0.0262452|
## |12 |  -85.46117|     0| 58.79282|  5| -1.453599| 0.0730288| 0.0568182| 0.0454545| 0.0162106| 0.0275742|   0.0275742|
print(kable(ks_calc, caption = "Cálculo Manual - Estadístico Kolmogorov-Smirnov"))
## 
## 
## Table: Cálculo Manual - Estadístico Kolmogorov-Smirnov
## 
## | Estadistico_D|
## |-------------:|
## |     0.0754392|

Prueba Shapiro-Wilk

df_sw_manual <- df_res %>%
  arrange(res) %>%
  mutate(
    i = row_number(),
    # Cuantiles esperados de una normal estándar usando la aproximación de Blom
    m_i = qnorm((i - 0.375) / (n + 0.25)) 
  )

sw_calc <- df_sw_manual %>%
  summarise(
    W_aprox = cor(res, m_i)^2
  )

print(kable(head(df_sw_manual, 5), caption = "Primeras 5 obs. - Proceso Manual SW"))
## 
## 
## Table: Primeras 5 obs. - Proceso Manual SW
## 
## |   |        res|  i|       m_i|
## |:--|----------:|--:|---------:|
## |81 | -120.02645|  1| -2.453069|
## |77 | -115.50870|  2| -2.087675|
## |24 | -107.08089|  3| -1.884554|
## |48 |  -91.24398|  4| -1.738328|
## |12 |  -85.46117|  5| -1.621941|
print(kable(sw_calc, caption = "Cálculo Manual - Estadístico W (Aproximado)"))
## 
## 
## Table: Cálculo Manual - Estadístico W (Aproximado)
## 
## |   W_aprox|
## |---------:|
## | 0.9374775|

Resultados en Forma Tabular:

jb_test <- jarque.bera.test(residuos)
ks_test <- ks.test(residuos, "pnorm", mean(residuos), sd(residuos)) 
sw_test <- shapiro.test(residuos)
resultados_oficiales <- data.frame(
  Prueba = c("Jarque-Bera (JB)", "Kolmogorov-Smirnov (KS)", "Shapiro-Wilk (SW)"),
  Estadistico = c(jb_test$statistic, ks_test$statistic, sw_test$statistic),
  P_Valor = c(jb_test$p.value, ks_test$p.value, sw_test$p.value)
)

print(kable(resultados_oficiales, caption = "Resultados de Pruebas Usando Librerías"))
## 
## 
## Table: Resultados de Pruebas Usando Librerías
## 
## |          |Prueba                  | Estadistico|   P_Valor|
## |:---------|:-----------------------|-----------:|---------:|
## |X-squared |Jarque-Bera (JB)        |  32.2779112| 0.0000001|
## |D         |Kolmogorov-Smirnov (KS) |   0.0754392| 0.6700386|
## |W         |Shapiro-Wilk (SW)       |   0.9413208| 0.0005937|

Gráficos para Jarque-Bera y Shapiro-Wilk

par(mfrow = c(1, 2))
shadeDist(
  xstat = jb_test$statistic, 
  ddist = "dchisq", 
  parm1 = 2, 
  lower.tail = FALSE, 
  main = "Prueba Jarque-Bera (Distribución Chi-cuadrado)"
)
z_sw <- qnorm(sw_test$p.value, lower.tail = FALSE)
shadeDist(
  xstat = z_sw, 
  ddist = "dnorm", 
  lower.tail = FALSE, 
  main = "Prueba Shapiro-Wilk (P-valor mapeado en Normal)"
)

par(mfrow = c(1, 1))