importacion

#Enuanciado: Sea el conjunto de datos, indicados en el enlace de abajo, tomados en 24 meses correspondientes a los gastos de comercialización (C) de una empresa, el nivel de ventas (V), su coste de personal (P) y los costes de materias primas (M); se trata de estimar el nivel de ventas a partir de las restantes variables.

library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
real_state_excel <- read_excel(
  "C:/Users/gabri/Downloads/ventas_empresa.xlsx",
  col_names = TRUE
)

# eliminar columnas completamente vacías
real_state_excel <- real_state_excel[, colSums(is.na(real_state_excel)) < nrow(real_state_excel)]

head(real_state_excel, 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

Estime el modelo de regresión lineal, correspondiente y verifique el supuesto de normalidad, usando todas las pruebas vistas en clase. EL CASO DEL JB Y DEL SW REPRESENTE SUS RESULTADOS DE FORMA GRÁFICA Y COMENTE AL RESPECTO.

modelo_estimado<- lm( V~ C + P + M, data = real_state_excel )
summary(modelo_estimado)
## 
## Call:
## lm(formula = V ~ C + P + M, data = real_state_excel)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.279  -6.966   1.555   6.721  14.719 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 107.4435    18.0575   5.950 8.08e-06 ***
## C             0.9226     0.2227   4.142 0.000505 ***
## P             0.9502     0.1558   6.097 5.86e-06 ***
## M             1.2978     0.4307   3.013 0.006872 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.506 on 20 degrees of freedom
## Multiple R-squared:  0.9798, Adjusted R-squared:  0.9768 
## F-statistic: 323.6 on 3 and 20 DF,  p-value: < 2.2e-16

Ajuste de los residuos a la Distribución Normal

library(fitdistrplus)
## Warning: package 'fitdistrplus' was built under R version 4.5.3
## Cargando paquete requerido: MASS
## Cargando paquete requerido: survival
## Warning: package 'survival' was built under R version 4.5.2
fit_normal<-fitdist(data = modelo_estimado$residuals,distr = "norm")
plot(fit_normal) 

options(scipen = 999)
summary(fit_normal)
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##                      estimate Std. Error
## mean 0.0000000000000004254771   1.771258
## sd   8.6773587050283431665321   1.252469
## Loglikelihood:  -85.91174   AIC:  175.8235   BIC:  178.1796 
## Correlation matrix:
##                   mean                sd
## mean 1.000000000000000 0.000000002955564
## sd   0.000000002955564 1.000000000000000

Verique el supuesto de normalidad, a través de:

a) La prueba JB

Usando tseries

library(tseries)
## Warning: package 'tseries' was built under R version 4.5.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
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
# construyendo la regla de decisión para saber: ✔️ Si los errores son normales ❌ o no son normales
library(fastGraph)
## Warning: package 'fastGraph' was built under R version 4.5.3
options(scipen = 9 )
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,xmax = max(JB, VC) + 5,col = "red", #El estadístico Jarque-Bera (JB = 32.28) es mayor que el valor crítico (5.99), por lo que se rechaza la hipótesis nula de normalidad. Esto indica que los residuos del modelo no siguen una distribución normal.
          sub=paste("VC:",round(VC,2)," ","JB:",round(JB,2)))

cat("JB =", salida_JB$statistic, "\n")
## JB = 1.400415
cat("p-value =", salida_JB$p.value, "\n")
## p-value = 0.4964822

Según la prueba de Jarque-Bera (JB = 1.4004; p-value = 0.4965), no se rechaza la hipótesis nula de normalidad de los residuos, por lo que se concluye que estos se distribuyen normalmente.

c) La prueba SW

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(gt)
## Warning: package 'gt' was built under R version 4.5.3
# Residuales del modelo
residuos <- modelo_estimado$residuals

# Tamaño de muestra
n <- length(residuos)

tabla_SW <- residuos %>%  
  as_tibble() %>%
  rename(residuales = value) %>%
  arrange(residuales) %>%
  mutate(
    pi = (row_number() - 0.375) / (n + 0.25),
    mi = qnorm(pi),
    ai = 0
  )

m <- sum(tabla_SW$mi^2)
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 <- tabla_SW %>%
  mutate(
    ai_ui = ai * residuales,
    ui2 = residuales^2
  )

tabla_SW %>%
  gt() %>%
  tab_header(title = "Tabla para calcular el Estadístico W") %>%
  tab_source_note(source_note = "Fuente: Elaboración propia")
Tabla para calcular el Estadístico W
residuales pi mi ai ai_ui ui2
-17.27950969 0.02577320 -1.94690278 -0.44751326 7.7328097356 298.5814551698
-15.50384667 0.06701031 -1.49843365 -0.31302439 4.8530821064 240.3692617022
-13.84401165 0.10824742 -1.23590240 -0.25502179 3.5305246473 191.6566585921
-9.91393148 0.14948454 -1.03864671 -0.21431914 2.1247453168 98.2860374647
-8.43534410 0.19072165 -0.87524006 -0.18060106 1.5234321207 71.1550300366
-7.90808316 0.23195876 -0.73241136 -0.15112913 1.1951417479 62.5377792373
-6.65267950 0.27319588 -0.60317579 -0.12446207 0.8280062569 44.2581445625
-3.33614422 0.31443299 -0.48332361 -0.09973122 0.3327177289 11.1298582608
-3.32226362 0.35567010 -0.37005675 -0.07635921 0.2536854377 11.0374355619
-2.43553165 0.39690722 -0.26136061 -0.05393035 0.1313490763 5.9318144083
-0.40251845 0.43814433 -0.15567569 -0.03212284 0.0129300352 0.1620211022
0.01280369 0.47938144 -0.05170609 -0.01066927 -0.0001366061 0.0001639345
3.09667909 0.52061856 0.05170609 0.01066927 0.0330393120 9.5894213810
4.04562361 0.56185567 0.15567569 0.03212284 0.1299569147 16.3670704012
4.98802235 0.60309278 0.26136061 0.05393035 0.2690057954 24.8803669567
5.88324503 0.64432990 0.37005675 0.07635921 0.4492399641 34.6125720380
6.28004614 0.68556701 0.48332361 0.09973122 0.6263166548 39.4389794795
6.69809559 0.72680412 0.60317579 0.12446207 0.8336588371 44.8644845291
6.78865323 0.76804124 0.73241136 0.15112913 1.0259632737 46.0858126562
7.37251119 0.80927835 0.87524006 0.18060106 1.3314833668 54.3539212565
8.70403920 0.85051546 1.03864671 0.21431914 1.8654422374 75.7602984046
9.77138853 0.89175258 1.23590240 0.25502179 2.4919170034 95.4800337177
10.67367778 0.93298969 1.49843365 0.31302439 3.3411214463 113.9273972668
14.71907878 0.97422680 1.94690278 0.44751326 6.5869829488 216.6512801774
Fuente: Elaboración propia

Cálculo del Estadistico W

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

Cálculo del 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
library(fastGraph)
shadeDist(Wn,ddist = "dnorm",lower.tail = FALSE)

Según la prueba de Shapiro-Wilk p-value = 0.3166 > 0.05 No se rechaza H₀(W = 0.9531; p-value = 0.3166), no se rechaza la hipótesis nula de normalidad, por lo que se concluye que los residuos del modelo siguen una distribución normal.

* Utilizando todas las herramientas vistas en clase, evalué la situación de colinealidad de los regresores del modelo. EN EL CASO DE LA LA PRUEBA DE FG, REPRESENTE SUS RESULTADOS GRÁFICAMENTE Y COMENTE AL RESPECTO.

#Cálculo “manual”:
library(stargazer)
## Warning: package 'stargazer' was built under R version 4.5.2
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
X_mat<-model.matrix(modelo_estimado)
stargazer(head(X_mat,n=6),type="text")
## 
## =========================
##   (Intercept)  C   P   M 
## -------------------------
## 1      1      197 173 110
## 2      1      208 152 107
## 3      1      181 150 99 
## 4      1      194 150 102
## 5      1      192 163 109
## 6      1      196 179 114
## -------------------------
XX_matrix<-t(X_mat)%*%X_mat
stargazer(XX_matrix,type = "text")
## 
## ===================================================
##             (Intercept)     C         P        M   
## ---------------------------------------------------
## (Intercept)     24        5,308     4,503    2,971 
## C              5,308    1,187,852 1,007,473 664,534
## P              4,503    1,007,473  859,157  564,389
## M              2,971     664,534   564,389  372,387
## ---------------------------------------------------

La matriz XtX debe ser normalizada, para evitar sesgo de la escala de las variables. Cálculo de la matriz de normalización:

library(stargazer)
options(scipen = 999)
Sn<-solve(diag(sqrt(diag(XX_matrix))))
stargazer(Sn,type = "text")
## 
## =======================
## 0.204   0     0     0  
## 0     0.001   0     0  
## 0       0   0.001   0  
## 0       0     0   0.002
## -----------------------

XtX normalizada

library(knitr)
## Warning: package 'knitr' was built under R version 4.5.2
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
options(scipen = 99)

XX_norm <- (Sn %*% XX_matrix) %*% Sn

kable(XX_norm,
      digits = 4,
      format = "html",
      align = "c",
      caption = "Matriz X'X normalizada") %>%
  kable_styling(bootstrap_options = c("bordered", "striped"),
                full_width = FALSE,
                position = "center")
Matriz X’X normalizada
1.0000 0.9941 0.9917 0.9938
0.9941 1.0000 0.9973 0.9992
0.9917 0.9973 1.0000 0.9978
0.9938 0.9992 0.9978 1.0000 
# Autovalores de XtX Normalizada:
library(knitr)
library(kableExtra)

lambdas <- eigen(XX_norm, symmetric = TRUE)

valores <- data.frame(
  i = 1:length(lambdas$values),
  Autovalores = lambdas$values
)

kable(valores,
      digits = 4,
      format = "html",
      align = "c") %>%   #  CENTRA TODO
  kable_styling(bootstrap_options = c("bordered", "striped"))
i Autovalores
1 3.9869
2 0.0095
3 0.0028
4 0.0008

Cálculo de κ(x)= √λmax/λmin

K<-sqrt(max(lambdas$values)/min(lambdas$values))
print(K)
## [1] 71.16349

El número de condición obtenido (κ = 71.16) indica la presencia de multicolinealidad fuerte entre las variables explicativas del modelo, lo cual puede afectar la estabilidad y precisión de los estimadores.

Prueba de Farrar-Glaubar

calculo manual Calculo de |R| c

library(knitr)
library(kableExtra)

Zn <- scale(X_mat[,-1])

kable(head(Zn, 6), digits = 3, format = "html") %>%
  kable_styling(bootstrap_options = c("bordered", "striped"))
C P M
-0.983 -0.587 -0.975
-0.536 -1.430 -1.187
-1.634 -1.510 -1.753
-1.105 -1.510 -1.541
-1.186 -0.988 -1.046
-1.024 -0.346 -0.692

Calcular la matriz R

library(knitr)
library(kableExtra)

n <- nrow(Zn)
R <- (t(Zn) %*% Zn) * (1/(n-1))

kable(R,
      digits = 4,
      format = "html",
      align = "c",
      caption = "Matriz  R") %>%
  kable_styling(bootstrap_options = c("bordered", "striped"),
                full_width = FALSE,
                position = "center")
Matriz R
C P M
C 1.0000 0.8205 0.9312
P 0.8205 1.0000 0.8579
M 0.9312 0.8579 1.0000

Calcular |R|

determinante_R<-det(R)
print(determinante_R)
## [1] 0.03459107

Aplicando la prueba de Farrer Glaubar (Bartlett)

m<-ncol(X_mat[,-1])
n<-nrow(X_mat[,-1])
chi_FG<--(n-1-(2*m+5)/6)*log(determinante_R)
print(chi_FG)
## [1] 71.20805

El estadístico de Farrar-Glauber (χ² = 71.21) es mayor que el valor crítico, por lo que se rechaza la hipótesis nula de ausencia de multicolinealidad. Se concluye que existe multicolinealidad significativa entre las variables explicativas del modelo.

Cálculo de los VIF’s usando “performance”

library(performance)
## Warning: package 'performance' was built under R version 4.5.3
VIFs<-multicollinearity(x = modelo_estimado,verbose = FALSE)
VIFs
## # Check for Multicollinearity
## 
## Low Correlation
## 
##  Term  VIF    VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##     P 3.84 [2.41,  6.73]     1.96      0.26     [0.15, 0.42]
## 
## Moderate Correlation
## 
##  Term  VIF    VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
##     C 7.63 [4.49, 13.61]     2.76      0.13     [0.07, 0.22]
##     M 9.45 [5.49, 16.91]     3.07      0.11     [0.06, 0.18]
library(performance)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
VIFs <- multicollinearity(modelo_estimado)
df_vif <- VIFs[, c("Term", "VIF")]

ggplot(df_vif, aes(x = Term, y = VIF)) +
  geom_col() +
  labs(title = "VIF por variable",
       x = "Variables",
       y = "VIF")

Cálculo de los VIF’s usando “car”```{r}

library(car)
## Warning: package 'car' was built under R version 4.5.3
## 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:dplyr':
## 
##     recode
VIFs_car<-vif(modelo_estimado)
print(VIFs_car)
##        C        P        M 
## 7.631451 3.838911 9.449210

Cálculo de los VIF’s usando “mctest”

library(mctest)
## Warning: package 'mctest' was built under R version 4.5.2
mc.plot(mod = modelo_estimado,vif = 2)

### * Verifique el supuesto de Homocedasticidad. REPRESENTE GRÁFICAMENTE SUS RESULTADOS Y COMENTE SUS RESULTADOS.

VERIFICAR SUPUESTO DE HOMOCEDASTICIDAD Prueba de white calculo manual

library(stargazer)
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.5.2
## Cargando paquete requerido: zoo
## Warning: package 'zoo' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# 1. Modelo original
modelo_ejemplo <- lm(V ~ C + P + M, data = real_state_excel)

# 2. Residuos
u_i <- modelo_ejemplo$residuals

# 3. Dataset auxiliar
data_prueba_white <- as.data.frame(cbind(u_i, real_state_excel))

# 4. Regresión auxiliar (White)
reg_aux <- lm(I(u_i^2) ~ C + P + M + 
                I(C^2) + I(P^2) + I(M^2) + 
                C:P + C:M + P:M,
              data = data_prueba_white)

resumen <- summary(reg_aux)

# 5. Estadístico LM
n   <- nrow(data_prueba_white)
R_2 <- resumen$r.squared
LM_w <- n * R_2

# Grados de libertad: 3 originales + 3 cuadrados + 3 interacciones = 9
gl <- 9 

# p-value y valor crítico
p_value <- 1 - pchisq(LM_w, df = gl)
VC <- qchisq(0.95, df = gl)

# 6. Tabla bonita
tabla_resultados <- data.frame(
  Estadistico_LM = LM_w,
  Valor_Critico  = VC,
  P_Value        = p_value
)

stargazer(tabla_resultados, 
          type = "text", 
          summary = FALSE, 
          title = "Prueba de White (Cálculo Manual)",
          digits = 6)
## 
## Prueba de White (Cálculo Manual)
## =======================================
##   Estadistico_LM Valor_Critico P_Value 
## ---------------------------------------
## 1    7.122650      16.918980   0.624351
## ---------------------------------------
# 7. Prueba rápida con lmtest (versión automática)
prueba_white_fast <- bptest(modelo_ejemplo, 
                            varformula = ~ C + P + M + 
                                           I(C^2) + I(P^2) + I(M^2) + 
                                           C:P + C:M + P:M, 
                            data = real_state_excel)

print(prueba_white_fast)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ejemplo
## BP = 7.1227, df = 9, p-value = 0.6244

calculo con libreria

library(lmtest)

prueba_white <- bptest(
  modelo_ejemplo,
  varformula = ~ C + P + M + 
                 I(C^2) + I(P^2) + I(M^2) + 
                 C:P + C:M + P:M,
  data = real_state_excel
)

print(prueba_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ejemplo
## BP = 7.1227, df = 9, p-value = 0.6244
# Residuos y ajustados
residuos <- residuals(modelo_estimado)
ajustados <- fitted(modelo_estimado)

# Pantalla 2x2
par(mfrow = c(2,2))

# 1. Residuos vs Ajustados (CLAVE)
plot(ajustados, residuos,
     main = "Residuos vs Ajustados",
     xlab = "Valores Ajustados",
     ylab = "Residuos",
     pch = 19, col = "blue")
abline(h = 0, lty = 2, col = "red")

# 2. Residuos vs C
plot(real_state_excel$C, residuos,
     main = "Residuos vs C",
     xlab = "C",
     ylab = "Residuos",
     pch = 19, col = "darkgreen")

# 3. Residuos vs P
plot(real_state_excel$P, residuos,
     main = "Residuos vs P",
     xlab = "P",
     ylab = "Residuos",
     pch = 19, col = "purple")

# 4. Residuos vs M
plot(real_state_excel$M, residuos,
     main = "Residuos vs M",
     xlab = "M",
     ylab = "Residuos",
     pch = 19, col = "brown")

# Residuos y valores ajustados
residuos <- residuals(modelo_estimado)
ajustados <- fitted(modelo_estimado)

# Gráfico
plot(ajustados, residuos,
     main = "Verificacion de homocedasticidad",
     xlab = "Valores Ajustados",
     ylab = "Residuos",
     pch = 19, col = "blue")

# Línea horizontal en 0
abline(h = 0, col = "red", lty = 2)

Según la prueba de Breusch-Pagan (White), no se rechaza la hipótesis nula de homocedasticidad (p-value = 0.6244), por lo que se concluye que no existe evidencia de heterocedasticidad en el modelo.

Verifique el supuesto de No Autocorrelación, (verifique primero y segundo orden). USE LA PRUEBA DEL MULTIPLICADOR DE LAGRANGE, REPRESENTE GRÁFICAMENTE SUS RESULTADOS Y COMENTE AL RESPECTO.

Verifique si los residuos del modelo son independientes entre sí (no autocorrelación), a través de: Prueba de Durbin Watson.

# Opción 1: Usando la librería "lmtest"
library(lmtest)
# alternative = "two.sided" busca autocorrelación tanto positiva como negativa
test_dw_1 <- dwtest(modelo_estimado, alternative = "two.sided", iterations = 1000)
print(test_dw_1)
## 
##  Durbin-Watson test
## 
## data:  modelo_estimado
## DW = 1.2996, p-value = 0.05074
## alternative hypothesis: true autocorrelation is not 0

autocorrelacion de orden superior (orden m)

# 1. Modelo original
modelo <- lm(V ~ C + P + M, data = real_state_excel)

# 2. Residuos (u_i)
u_i <- resid(modelo)

# 3. Crear dataset para BG
data_prueba_BG <- real_state_excel
data_prueba_BG$u_i <- u_i
data_prueba_BG$Lag_1 <- dplyr::lag(u_i, 1)
data_prueba_BG$Lag_2 <- dplyr::lag(u_i, 2)

# 4. Modelo BG
bg_2 <- lm(u_i ~ C + P + M + Lag_1 + Lag_2,
           data = data_prueba_BG)

summary(bg_2)
## 
## Call:
## lm(formula = u_i ~ C + P + M + Lag_1 + Lag_2, data = data_prueba_BG)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.1290  -5.5254  -0.4024   6.4035  14.5180 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.51115   18.62219  -0.511    0.617
## C            0.05816    0.25231   0.231    0.821
## P            0.06608    0.16489   0.401    0.694
## M           -0.13272    0.48232  -0.275    0.787
## Lag_1        0.39436    0.26327   1.498    0.154
## Lag_2       -0.23690    0.25067  -0.945    0.359
## 
## Residual standard error: 9.294 on 16 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.1491, Adjusted R-squared:  -0.1169 
## F-statistic: 0.5605 on 5 and 16 DF,  p-value: 0.7287
library(dplyr)
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.5.2
library(kableExtra)
# Usamos los residuos u_i que ya extrajiste de tu modelo_ejemplo
data_prueba_BG <- cbind(u_i, real_state_excel) %>%
  as.data.frame() %>%
  mutate(Lag_1 = dplyr::lag(u_i, 1),
         Lag_2 = dplyr::lag(u_i, 2)) %>%
  replace_na(list(Lag_1 = 0, Lag_2 = 0))

# Ver las primeras filas para verificar los rezagos
kable(head(data_prueba_BG, 6))
u_i V C P M Lag_1 Lag_2
10.673678 607 197 173 110 0.000000 0.000000
7.372511 590 208 152 107 10.673678 0.000000
-2.435532 543 181 150 99 7.372511 10.673678
-3.322264 558 194 150 102 -2.435532 7.372511
-9.913932 571 192 163 109 -3.322264 -2.435532
8.704039 615 196 179 114 -9.913932 -3.322264 
# 2. Calculando la regresión auxiliar
# Se regresan los residuos contra las variables originales + los rezagos
regresion_auxiliar_BG <- lm(u_i ~ C + P + M + Lag_1 + Lag_2, 
                             data = data_prueba_BG)

sumario_BG <- summary(regresion_auxiliar_BG)

# 3. Cálculo del estadístico LM (Breusch-Godfrey)
n_bg <- nrow(data_prueba_BG)
R2_bg <- sumario_BG$r.squared
LM_bg <- n_bg * R2_bg

# Grados de libertad = número de rezagos (en este caso 2)
gl_bg <- 2
p_valor_bg <- 1 - pchisq(LM_bg, df = gl_bg)

# Mostrar resultado rápido
cat("Estadístico LM:", LM_bg, "\nP-valor:", p_valor_bg)
## Estadístico LM: 3.840869 
## P-valor: 0.1465433

b) Prueba del Multiplicador de Lagrange (verifique autocorrelación de primer y segundo orden).

autocorrelación de primer orden

bg_1 <- lm(u_i ~ C + P + M + Lag_1, 
           data = data_prueba_BG)

# LM
n1 <- nrow(data_prueba_BG)
R2_1 <- summary(bg_1)$r.squared
LM_1 <- n1 * R2_1

# p-valor
p_1 <- 1 - pchisq(LM_1, df = 1)

cat("Orden 1 -> LM:", LM_1, "p-valor:", p_1)
## Orden 1 -> LM: 2.596321 p-valor: 0.1071121

autocorrelación de segundo orden

# Regresión auxiliar (orden 2)
bg_2 <- lm(u_i ~ C + P + M + Lag_1 + Lag_2, 
           data = data_prueba_BG)

# LM
n2 <- nrow(data_prueba_BG)
R2_2 <- summary(bg_2)$r.squared
LM_2 <- n2 * R2_2

# p-valor
p_2 <- 1 - pchisq(LM_2, df = 2)

cat("Orden 2 -> LM:", LM_2, "p-valor:", p_2)
## Orden 2 -> LM: 3.840869 p-valor: 0.1465433

Según la prueba del Multiplicador de Lagrange (Breusch-Godfrey), no se rechaza la hipótesis nula de ausencia de autocorrelación tanto de primer como de segundo orden (p-values mayores a 0.05). Por lo tanto, se concluye que no existe evidencia de autocorrelación en los errores del modelo.