###cargar los datos

library(stats)
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
# Esto abrirá una ventana, busca tu archivo Excel y selecciónalo
ventas_empresa <- read_excel("C:\\Users\\MINEDUCYT\\Desktop\\EMA118\\SEGUNDO EXAMEN PARCIAL\\ventas_empresa.xlsx")

# Verificas
head(ventas_empresa, n=10)
## # A tibble: 10 × 4
##        V     C     P     M
##    <dbl> <dbl> <dbl> <dbl>
##  1   607   197   173   110
##  2   590   208   152   107
##  3   543   181   150    99
##  4   558   194   150   102
##  5   571   192   163   109
##  6   615   196   179   114
##  7   606   203   169   113
##  8   593   200   166   113
##  9   582   198   159   115
## 10   646   221   206   119

##Estimación del modelo

library(stargazer)
modelo_estimado <- lm(V ~ C + P + M, data = ventas_empresa)
stargazer(modelo_estimado,title = "MODELO_ESTIMADO",type = "text")
## 
## MODELO
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                  V             
## -----------------------------------------------
## C                            0.923***          
##                               (0.223)          
##                                                
## P                            0.950***          
##                               (0.156)          
##                                                
## M                            1.298***          
##                               (0.431)          
##                                                
## Constant                    107.444***         
##                              (18.057)          
##                                                
## -----------------------------------------------
## Observations                    24             
## R2                             0.980           
## Adjusted R2                    0.977           
## Residual Std. Error       9.506 (df = 20)      
## F Statistic           323.641*** (df = 3; 20)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

#Pruebas de normalidad ##Prueba de Jarque_Bera ##Usando la libreria tseries

options(scipen = 9999)
library(tseries)
salida_JB<-jarque.bera.test(modelo_estimado$residuals)
salida_JB
## 
##  Jarque Bera Test
## 
## data:  modelo_estimado$residuals
## X-squared = 1.4004, df = 2, p-value = 0.4965
#Resultado: Dado que el p-value=1.4004 y el nivel de signficancia es de 0.05 entonces 1.4004>0.05 por lo tanto no se se rechaza la hipotesis nula y se concluye que los datos siguen una distribución normal

##Calculo manual de la prueba Jarque-Bera

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(gt)
# 1. Extraer residuales y definir n
residuos<- modelo_estimado$residuals
n <- length(residuos)

# 2. Calcular los componentes de la prueba
# Asimetría (Skewness)
asimetria <- sum((residuos - mean(residuos))^3) / (n * (sd(residuos) * sqrt((n-1)/n))^3)

# Curtosis
curtosis <- sum((residuos - mean(residuos))^4) / (n * (sd(residuos) * sqrt((n-1)/n))^4)

# Estadístico Jarque-Bera (JB)
jb_est <- (n / 6) * (asimetria^2 + ((curtosis - 3)^2 / 4))

# Valor p (usa una distribución Chi-cuadrado con 2 grados de libertad)
p_val_jb <- pchisq(jb_est, df = 2, lower.tail = FALSE)

# 3. Crear la tabla de resultados
tabla_JB <- tibble(
  Estadistico = c("Asimetría", "Curtosis", "Jarque-Bera (JB)", "Valor p"),
  Valor = c(asimetria, curtosis, jb_est, p_val_jb),
  Interpretacion = c(
    ifelse(asimetria > 0, "Positiva", "Negativa"),
    ifelse(curtosis > 3, "Leptocúrtica", "Platicúrtica"),
    "Estadístico de prueba",
    ifelse(p_val_jb > 0.05, "No Rechaza Normalidad", "Rechaza Normalidad")
  )
)

# 4. Formatear con gt
tabla_JB %>%
  gt() %>%
  tab_header(title = "Resultados de la Prueba Jarque-Bera",
             subtitle = "Evaluación de Asimetría y Curtosis de los Residuales") %>%
  fmt_number(columns = Valor, decimals = 4) %>%
  cols_label(Estadistico = "Métrica", Valor = "Valor Calculado")
Resultados de la Prueba Jarque-Bera
Evaluación de Asimetría y Curtosis de los Residuales
Métrica Valor Calculado Interpretacion
Asimetría −0.4023 Negativa
Curtosis 2.1322 Platicúrtica
Jarque-Bera (JB) 1.4004 Estadístico de prueba
Valor p 0.4965 No Rechaza Normalidad

##Gráfico de la prueba Jarque-BERA

library(fastGraph)
## Warning: package 'fastGraph' was built under R version 4.5.3
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)))

##Calculo de la prueba de Kolmogorov Smimov-Liliefors ###libreria nortest

library(nortest)
prueba_KS<-lillie.test(modelo_estimado$residuals)
prueba_KS
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  modelo_estimado$residuals
## D = 0.13659, p-value = 0.2935
#Resultado, dado que p-value=0.2935> 0.05 entonces no se rechaza Ho, entonces los datos siguen una distribución normal

#Calculo manual de la prueba K-S

library(dplyr)  
library(gt) 
library(gtExtras)
## Warning: package 'gtExtras' was built under R version 4.5.3
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
  
  #Formato
 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)
    ))
Tabla para calcular el Estadistico KS
residuales posicion dist1 dist2 zi pi dif1 dif2
-17.27950969 11 0.04166667 0.00000000 -1.949405386 0.02562352 0.0160431507 0.025623516
-15.50384667 10 0.08333333 0.04166667 -1.749082164 0.04013841 0.0431949244 0.001528258
-13.84401165 19 0.12500000 0.08333333 -1.561826195 0.05916447 0.0658355305 0.024168864
-9.91393148 5 0.16666667 0.12500000 -1.118450221 0.13168738 0.0349792889 0.006687378
-8.43534410 9 0.20833333 0.16666667 -0.951641887 0.17063932 0.0376940180 0.003972649
-7.90808316 24 0.25000000 0.20833333 -0.892158410 0.18615402 0.0638459847 0.022179318
-6.65267950 23 0.29166667 0.25000000 -0.750528776 0.22646815 0.0651985168 0.023531850
-3.33614422 8 0.33333333 0.29166667 -0.376370489 0.35332074 0.0199874078 0.061654074
-3.32226362 4 0.37500000 0.33333333 -0.374804535 0.35390292 0.0210970793 0.020569587
-2.43553165 3 0.41666667 0.37500000 -0.274766969 0.39174764 0.0249190286 0.016747638
-0.40251845 21 0.45833333 0.41666667 -0.045410527 0.48189005 0.0235567120 0.065223379
0.01280369 14 0.50000000 0.45833333 0.001444461 0.50057626 0.0005762565 0.042242923
3.09667909 20 0.54166667 0.50000000 0.349354988 0.63658859 0.0949219225 0.136588589
4.04562361 7 0.58333333 0.54166667 0.456411125 0.67595282 0.0926194830 0.134286150
4.98802235 12 0.62500000 0.58333333 0.562728793 0.71319021 0.0881902126 0.129856879
5.88324503 22 0.66666667 0.62500000 0.663724246 0.74656659 0.0798999242 0.121566591
6.28004614 18 0.70833333 0.66666667 0.708489765 0.76067942 0.0523460837 0.094012750
6.69809559 15 0.75000000 0.70833333 0.755652438 0.77507120 0.0250711962 0.066737863
6.78865323 16 0.79166667 0.75000000 0.765868789 0.77812281 0.0135438559 0.028122811
7.37251119 2 0.83333333 0.79166667 0.831737317 0.79722138 0.0361119495 0.005554717
8.70403920 6 0.87500000 0.83333333 0.981954997 0.83693899 0.0380610096 0.003605657
9.77138853 13 0.91666667 0.87500000 1.102369093 0.86484938 0.0518172880 0.010150621
10.67367778 1 0.95833333 0.91666667 1.204161768 0.88573647 0.0725968632 0.030930197
14.71907878 17 1.00000000 0.95833333 1.660547780 0.95159785 0.0484021522 0.006735485
Fuente: Elaboración propia

##Cálculo del estadistico SK

D<-max(max(tabla_KS$dif1),max(tabla_KS$dif2))
print(D)
## [1] 0.1365886
##Resultado: dado que 0.1365<0.2693. Es decir que D<V.C entonces la hipótesis nula no se rechaza y por lo tanto, los datos siguen una distribución normal

##Cálculo de la prueba Shapiro Will ##Usando la libreria stats

salida_SW<-shapiro.test(modelo_estimado$residuals)
print(salida_SW)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo_estimado$residuals
## W = 0.95315, p-value = 0.3166
Wn_salida<-qnorm(salida_SW$p.value,lower.tail = FALSE)
print(Wn_salida)
## [1] 0.4773519
#Conclusión: Dado que 0.3166(p-value)>0.05 entonces no se rechaza la hipótesis nula y los datos siguen una distribución normal

##calculo manual de la prueba Shapiro Will

library(dplyr) # Carga la librería para manipulación de datos (pipes, mutate, etc.)
library(gt)    # Carga la librería para la creación de tablas con formato profesional
residuos <-modelo_estimado$residuals # Extrae los residuales del objeto del modelo previamente ajustado
residuos |>   # Toma el vector de residuales y lo pasa al siguiente paso
  as_tibble() |>  # Convierte el vector en una estructura de datos tipo tabla (tibble)
  rename(residuales = value) |>   # Cambia el nombre de la columna por defecto 'value' a 'residuales'
  arrange(residuales) |>   # Ordena los residuales de menor a mayor (requisito de Shapiro-Wilk)
  mutate(pi = (row_number() - 0.375) / (n() + 0.25)) |>   # Calcula la posición de graficación (Blom's plotting position)
  mutate(mi = qnorm(pi, lower.tail = TRUE)) |>   # Calcula los cuantiles teóricos (esperados) de una normal estándar
  mutate(ai = 0) -> tabla_SW # Inicializa la columna 'ai' en cero y guarda el resultado en 'tabla_SW'

m <-sum(tabla_SW$mi^2) # Calcula la suma de los cuadrados de los cuantiles teóricos (m)
n <-nrow(tabla_SW) # Obtiene el tamaño de la muestra original
theta<-1/sqrt(n) # Define una variable auxiliar basada en la raíz inversa de n para las aproximaciones

# Calcula el coeficiente 'ai' para el último dato usando el polinomio de aproximación de Royston
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)

# Calcula el coeficiente 'ai' para el penúltimo dato usando la fórmula polinómica correspondiente
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] # El primer coeficiente es el negativo del último (por simetría)
tabla_SW$ai[2] <- -tabla_SW$ai[n-1] # El segundo coeficiente es el negativo del penúltimo

# Calcula un factor de escala (omega) para los coeficientes restantes basándose en m y los extremos calculados
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)

# Asigna los coeficientes 'ai' intermedios escalándolos por la raíz de omega
tabla_SW$ai[3:(n-2)] <- tabla_SW$mi[3:(n-2)] / sqrt(omega)

# Calcula los productos (ai * residual) y los residuales al cuadrado para el estadístico W final
tabla_SW |>   
  mutate(ai_ui = ai * residuales, ui2 = residuales^2) -> tabla_SW # Actualiza la tabla con los cálculos de los productos y cuadrados
print(tabla_SW)
## # A tibble: 24 × 6
##    residuales     pi     mi      ai ai_ui    ui2
##         <dbl>  <dbl>  <dbl>   <dbl> <dbl>  <dbl>
##  1     -17.3  0.0258 -1.95  -0.448  7.73  299.  
##  2     -15.5  0.0670 -1.50  -0.313  4.85  240.  
##  3     -13.8  0.108  -1.24  -0.255  3.53  192.  
##  4      -9.91 0.149  -1.04  -0.214  2.12   98.3 
##  5      -8.44 0.191  -0.875 -0.181  1.52   71.2 
##  6      -7.91 0.232  -0.732 -0.151  1.20   62.5 
##  7      -6.65 0.273  -0.603 -0.124  0.828  44.3 
##  8      -3.34 0.314  -0.483 -0.0997 0.333  11.1 
##  9      -3.32 0.356  -0.370 -0.0764 0.254  11.0 
## 10      -2.44 0.397  -0.261 -0.0539 0.131   5.93
## # ℹ 14 more rows

##calclo del estadistico w de la prueba shapiro

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

##calculo del estadistico Wn y 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.4773519
p.value<-pnorm(Wn,lower.tail = FALSE)
print(p.value)
## [1] 0.3165558
#Dado que el pvalue=0.3166>0.05 entonces no se rechaza la hipotesis nula y los datos siguen una distribución normal

##Grafico de la prueba Shapiro-will

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

#PRUEBAS DE MULTICOLANIEDAD ##CALCULO DEL INDICE DE CONDICIÓN CON LA LIBRERIA MCTEST

library(mctest)
X_mat<-model.matrix(modelo_estimado)
mctest(mod = modelo_estimado)
## 
## Call:
## omcdiag(mod = mod, Inter = TRUE, detr = detr, red = red, conf = conf, 
##     theil = theil, cn = cn)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0346         0
## Farrar Chi-Square:        71.2080         1
## Red Indicator:             0.8711         1
## Sum of Lambda Inverse:    20.9196         1
## Theil's Method:            0.5430         1
## Condition Number:         71.1635         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test
#dado que el indice de condición k(x)=71.2080 y es superior a 30. Entonces como k(x)> 30 se considera colinealidad severa

##calculo del indice usando la libreria olsrr

library(olsrr)
## 
## Adjuntando el paquete: 'olsrr'
## The following object is masked from 'package:datasets':
## 
##     rivers
ols_eigen_cindex(model = modelo_estimado)
##     Eigenvalue Condition Index    intercept             C            P
## 1 3.9869237681         1.00000 0.0007213226 0.00009632167 0.0002714936
## 2 0.0095007154        20.48523 0.8775733849 0.00494748750 0.0877003286
## 3 0.0027882470        37.81406 0.1182991896 0.15940219484 0.8478049007
## 4 0.0007872695        71.16349 0.0034061029 0.83555399599 0.0642232771
##               M
## 1 0.00008216652
## 2 0.00754342318
## 3 0.06362314121
## 4 0.92875126909
#dado que el indice de condición k(x)=71.2080 y es superior a 30. Entonces como k(x)> 30 se considera colinealidad severa

##PRUEBA DE Farrar-Glauber ##usando la libreria physis

library(psych)
FG_test<-cortest.bartlett(X_mat[,-1])
## R was not square, finding R from data
print(FG_test)
## $chisq
## [1] 71.20805
## 
## $p.value
## [1] 0.000000000000002352605
## 
## $df
## [1] 3

##REGLA DE DESCISIÓN: DADO QUE EL pvalue=0.000000<0.05 se rechaza la hipotesis nula y por lo tanto hay evidencia de colinealidad entre las variables explicativas. Otra forma de comprobarlo es mediante la regla de decisión utilizando el valor critico. si FG>= V.C se rechaza H0. Dado que FG= 71.81 y con tres grados de libertad de y un 5% de significancia en la tabla chi cuadrado el V.C= 7.81— 71.81>7.81 entonces se rechaza H0.

calculo de la prueba Farrar-Glauber usando la libreria mctest

library(mctest)
mctest::omcdiag(mod = modelo_estimado)
## 
## Call:
## mctest::omcdiag(mod = modelo_estimado)
## 
## 
## Overall Multicollinearity Diagnostics
## 
##                        MC Results detection
## Determinant |X'X|:         0.0346         0
## Farrar Chi-Square:        71.2080         1
## Red Indicator:             0.8711         1
## Sum of Lambda Inverse:    20.9196         1
## Theil's Method:            0.5430         1
## Condition Number:         71.1635         1
## 
## 1 --> COLLINEARITY is detected by the test 
## 0 --> COLLINEARITY is not detected by the test

##GRAFICANDO LA PRUEBA DE FARRAR

library(fastGraph)
shadeDist(xshade = 71.208, 
          ddist = "dchisq", 
          parm1 = 3, 
          lower.tail = FALSE, 
          main = "Prueba Global de Farrar-Glauber",
          sub = "Distribución Chi-cuadrado con df=3")

alpha_sig <- 0.05
farrar_stat <- 71.20805  
gl_farrar <- 3          
VC <- qchisq(1 - alpha_sig, gl_farrar)
shadeDist(xshade = farrar_stat, 
          ddist = "dchisq", 
          parm1 = gl_farrar, 
          lower.tail = FALSE, 
          xmin = 0, 
          xmax = max(20, farrar_stat + 5), 
          main = "Prueba Global de Multicolinealidad (Farrar-Glauber)",
          sub = paste("VC:", round(VC, 2), " | Estadístico Farrar:", round(farrar_stat, 2)))

#VERFICACIÓN DEL SUPUESTO DE HOMOCEDASTICIDAD

library(lmtest)
prueba_white <- bptest(modelo_estimado, 
                       ~ C + P + M + 
                         I(C^2) + I(P^2) + I(M^2) + 
                         C*P + C*M + M*P, 
                       data = ventas_empresa)
print(prueba_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_estimado
## BP = 7.1227, df = 9, p-value = 0.6244

##Dado que el p-value=0.6244>0.05 entonces no se rechaza H0 ENTONCES hay evidencia de que la varianza de los residuos es homocedástica. Si requiere hacerlo por medio de V.C entonces daber los gl o df, luego revisar la tabla y aplicar la regla de descisión. RECHAZAR H0 SI: LMn> V.C—- BP=LMn

library(fastGraph)

# Extrae los datos de tu objeto 'prueba_white'
est_white <- prueba_white$statistic
gl_white <- prueba_white$parameter # Debería ser 9

shadeDist(xshade = est_white, 
          ddist = "dchisq", 
          parm1 = gl_white, 
          lower.tail = FALSE, 
          main = "Prueba de Heterocedasticidad de White",
          sub = paste("Estadístico:", round(est_white, 2), " | p-valor:", round(prueba_white$p.value, 4)))

##Supuestos de no autocorrelación #SUPUESTO DE AUTOCORRELACIÓN DE PRIMER ORDEN ##PRUEBA DE DURBIN WATSON

library(lmtest)
dwtest(modelo_estimado,alternative="two.sided",iterations = 1000)
## 
##  Durbin-Watson test
## 
## data:  modelo_estimado
## DW = 1.2996, p-value = 0.05074
## alternative hypothesis: true autocorrelation is not 0

###Dado que el p-valor (0.0507) es mayor al nivel de significancia (0.05), no existe evidencia estadística suficiente para rechazar la hipótesis nula (H_0). Por lo tanto, con un nivel de confianza del 95%, se concluye que no hay problemas de autocorrelación en el modelo

##Prueba con la libreria Car

library(car)
## Warning: package 'car' was built under R version 4.5.2
## Cargando paquete requerido: carData
## Warning: package 'carData' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
## The following object is masked from 'package:dplyr':
## 
##     recode
durbinWatsonTest(modelo_estimado,simulate= TRUE,reps=1000)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.3013888      1.299572    0.06
##  Alternative hypothesis: rho != 0

##Dado que en los dos casos el p-valu>0.05 entonces se rechaza H0. ##Tomando en cuenta dL=1.101 y du=1.656 entonces nuetro estadistico cae en zona no conluyente, sin embargo, nuevo p-value nos hace rechazar h0.

#GRAFICO DE LA PRUEBA DURVIN WATSON

# Script para visualizar tu posición en la prueba DW
plot(0, 0, type="n", xlim=c(0,4), ylim=c(0,1), xlab="Estadístico DW", ylab="", 
     main="Zonas de Decisión Durbin-Watson", yaxt="n")

# Dibujar las zonas
rect(0, 0, 1.10, 1, col="red", border=NA)      # Rechazo (Positiva)
rect(1.10, 0, 1.65, 1, col="gray", border=NA)   # Indecisión
rect(1.65, 0, 2.35, 1, col="green", border=NA)  # No Autocorrelación
rect(2.35, 0, 2.90, 1, col="gray", border=NA)   # Indecisión
rect(2.90, 0, 4, 1, col="red", border=NA)       # Rechazo (Negativa)

# Marcar tu valor (1.29)
abline(v = 1.2996, col = "blue", lwd = 3, lty = 2)
text(1.2996, 0.5, "Tu DW = 1.29", pos = 4, col="blue", font=2)

#P-VALUE

library(fastGraph)

# Parámetros basados en tu resultado
# El p-valor es 0.0507. 
# En una distribución normal estándar, esto corresponde a un valor Z de aprox 1.64

shadeDist(xshade = 1.64, 
          ddist = "dnorm", 
          parm1 = 0, 
          parm2 = 1, 
          lower.tail = FALSE, 
          main = "Representación del p-valor (Durbin-Watson)",
          sub = "p-valor = 0.0507 | Decisión: No se rechaza H0")

#Autocorrelación de orden superior ##Prueba del Multiplicador de Lagrange [Breusch-Godfrey] ##usando la libreria lmtest

library(lmtest)
bgtest(modelo_estimado,order = 2)
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  modelo_estimado
## LM test = 3.8409, df = 2, p-value = 0.1465

#GRAFICO DE PRUEBA DE AUTOCORRELACIÓN ORDEN 2

library(fastGraph)

# Usamos tus resultados: LM = 3.84 y df = 2
shadeDist(xshade = 3.8409, 
          ddist = "dchisq", 
          parm1 = 2, 
          lower.tail = FALSE, 
          main = "Prueba de Breusch-Godfrey (Orden 2)",
          sub = paste("Estadístico LM:", 3.84, "| p-valor:", 0.1465))

#Dado que el p-value>0.05 entonces puede concluirse que, los residuos del modelo no siguen autocorrelación de orden 2

##orden 1

library(lmtest)
bgtest(modelo_estimado,order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo_estimado
## LM test = 2.5963, df = 1, p-value = 0.1071

#Dado que el p-value>0.05 entonces puede concluirse que, los residuos del modelo no siguen autocorrelación de orden 1