library(wooldridge)

data("hprice1")
head(force(hprice1), n=5)

# Mostrar las primeras 5 observaciones

library(stargazer)
modelo_estimado <- lm(formula =price~lotsize+sqrft+bdrms, data = hprice1)
stargazer(modelo_estimado, title = "Modelo estimado", type = "text")

library(tseries)
salida_jb <- jarque.bera.test(modelo_estimado$residuals)
salida_jb

# Librerías necesarias
library(tseries)
library(fastGraph)

# Obtener residuos del modelo
residuos <- residuals(modelo_estimado)

# Realizar la prueba de Jarque-Bera
salida_JB <- jarque.bera.test(residuos)

# Parámetros para la gráfica
alpha_sig <- 0.05
JB <- salida_JB$statistic
gl <- salida_JB$parameter
VC <- qchisq(1 - alpha_sig, gl, lower.tail = TRUE)

# Gráfica de la región crítica
shadeDist(JB,
          ddist = "dchisq",
          parm1 = gl,
          lower.tail = FALSE,
          xmin = 0,
          sub = paste("VC:", round(VC, 2), " ", "JB:", round(JB, 2)))

library(tseries)
salida_jb <- jarque.bera.test(modelo_estimado$residuals)
salida_jb

library(dplyr)
library(gt)
library(gtExtras)
residuos<-modelo_estimado$residuals  
residuos %>%  
  as_tibble() %>%  
  mutate(posicion=row_number()) %>%  
  arrange(value) %>%  
  mutate(dist1=row_number()/n()) %>%  
  mutate(dist2=(row_number()-1)/n()) %>%  
  mutate(zi=as.vector(scale(value,center=TRUE))) %>%  
  mutate(pi=pnorm(zi,lower.tail = TRUE)) %>%  
  mutate(dif1=abs(dist1-pi)) %>%  
  mutate(dif2=abs(dist2-pi)) %>%  
  rename(residuales=value) -> tabla_KS  


 tabla_KS %>%  
  gt() %>%  
  tab_header("Tabla para calcular el Estadistico KS") %>%  
  tab_source_note(source_note = "Fuente: Elaboración propia") %>%  
  tab_style(  
    style = list(
      cell_fill(color = "#A569BD"),  
      cell_text(style = "italic")  
      ),
    locations = cells_body(  
      columns = dif1,  
      rows = dif1==max(dif1)  
    )) %>%
   tab_style(  
    style = list(
      cell_fill(color = "#3498DB"),  
      cell_text(style = "italic")  
      ),
    locations = cells_body(  
      columns = dif2,  
      rows = dif2==max(dif2)  
    ))

D <- max(max(tabla_KS$dif1),max(tabla_KS$dif2))
print(D)

matriz_X <- model.matrix(modelo_estimado)
n <- nrow(matriz_X)
print(n)

n_KS <- 88
resultado_ks <- 0.855897/sqrt(n_KS)
print(resultado_ks)

library(nortest)
prueba_KS <- lillie.test(modelo_estimado$residuals)
prueba_KS

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_SW<-nrow(hprice1)
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") %>%  
  tab_source_note(source_note = "Fuente: Elaboración propia")

W <- (sum(tabla_SW$ai_ui)^2)/sum(tabla_SW$ui2)
print(W)

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)

p.value <- pnorm(Wn, lower.tail = FALSE)
print(p.value)

library(fastGraph)
shadeDist(Wn,ddist = "dnorm",lower.tail = FALSE)

salida_SW <- shapiro.test(modelo_estimado$residuals)
print(salida_SW)

Wn_salida <- qnorm(salida_SW$p.value, lower.tail = FALSE)
print(Wn_salida)