library(readxl)
ventas_empresa <- read_excel("C:/Users/Rebeca/Downloads/ventas_empresa.xlsx")
View(ventas_empresa)
# 1. CARGA DE LIBRERÍAS Y DATOS
library(dplyr)      
library(gt)         
library(tseries)    # Jarque-Bera
library(nortest)    # Lilliefors
library(fastGraph)  # Gráficas
library(stargazer)  # Tablas de regresión
#ESTIMACIÓN DEL MODELO
modelo_ventas <- lm(V ~ C + P + M, data = ventas_empresa)

# Resultados del modelo
stargazer(modelo_ventas, title = "Modelo de Regresión de Ventas", type = "text")
## 
## Modelo de Regresión de Ventas
## ===============================================
##                         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
# Extracción de residuos
residuos <- modelo_ventas$residuals
# 3. VERIFICACIÓN DEL SUPUESTO DE NORMALIDAD
# A) PRUEBA JARQUE-BERA (JB)
options(scipen = 999)
salida_JB <- jarque.bera.test(residuos)
print(salida_JB)
## 
##  Jarque Bera Test
## 
## data:  residuos
## X-squared = 1.4004, df = 2, p-value = 0.4965
# Gráfica JB (Distribución Chi-cuadrado)
alpha_sig <- 0.05
JB_stat <- salida_JB$statistic
gl_jb <- salida_JB$parameter
VC_jb <- qchisq(1 - alpha_sig, gl_jb)

shadeDist(JB_stat, ddist = "dchisq", parm1 = gl_jb, lower.tail = FALSE, xmin = 0,
          sub = paste("VC:", round(VC_jb, 2), " ", "JB:", round(JB_stat, 2)))

# B) PRUEBA LILLIEFORS (KS)
prueba_KS <- lillie.test(residuos)
print(prueba_KS)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.13659, p-value = 0.2935
# Tabla de cálculo manual (Estructura solicitada)
tabla_KS <- residuos %>%
  as_tibble() %>%
  mutate(posicion = row_number()) %>%
  arrange(value) %>%
  mutate(dist1 = row_number() / n(),
         dist2 = (row_number() - 1) / n(),
         zi = as.vector(scale(value, center = TRUE)),
         pi = pnorm(zi),
         dif1 = abs(dist1 - pi),
         dif2 = abs(dist2 - pi))

# Mostrar tabla con destacados
tabla_KS %>%
  gt() %>%
  tab_header("Tabla para cálculo del Estadístico KS (Lilliefors)") %>%
  tab_style(style = cell_fill(color = "lightblue"),
            locations = cells_body(columns = dif1, rows = dif1 == max(dif1)))
Tabla para cálculo del Estadístico KS (Lilliefors)
value 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
 
# C) PRUEBA SHAPIRO-WILK (SW)
# ---------------------------------------------------------
# 1. Preparación de tabla manual y coeficientes ai
n_obs <- length(residuos)
residuos %>%    
  as_tibble() %>%  
  rename(residuales = value) %>%  
  arrange(residuales) %>%  
  mutate(pi = (row_number() - 0.375) / (n() + 0.25)) %>%  
  mutate(mi = qnorm(pi)) %>%   
  mutate(ai = 0) -> tabla_SW

# Cálculo de ai (Aproximación de Royston para n=24)
m_val <- sum(tabla_SW$mi^2)
theta <- 1 / sqrt(n_obs)

# Coeficientes extremos
tabla_SW$ai[n_obs] <- -2.706056*theta^5 + 4.434685*theta^4 - 2.071190*theta^3 - 0.147981*theta^2 + 0.2211570*theta + tabla_SW$mi[n_obs]/sqrt(m_val)
tabla_SW$ai[n_obs-1] <- -3.582633*theta^5 + 5.682633*theta^4 - 1.752461*theta^3 - 0.293762*theta^2 + 0.042981*theta + tabla_SW$mi[n_obs-1]/sqrt(m_val)
tabla_SW$ai[1] <- -tabla_SW$ai[n_obs]
tabla_SW$ai[2] <- -tabla_SW$ai[n_obs-1]

# Valores intermedios
omega <- (m_val - 2*tabla_SW$mi[n_obs]^2 - 2*tabla_SW$mi[n_obs-1]^2) / (1 - 2*tabla_SW$ai[n_obs]^2 - 2*tabla_SW$ai[n_obs-1]^2)
tabla_SW$ai[3:(n_obs-2)] <- tabla_SW$mi[3:(n_obs-2)] / sqrt(omega)

# Estadístico W manual
tabla_SW %>% mutate(ai_ui = ai * residuales, ui2 = residuales^2) -> tabla_SW
W_manual <- (sum(tabla_SW$ai_ui)^2) / sum(tabla_SW$ui2)

# Normalización Wn para la gráfica
mu_sw <- 0.0038915*log(n_obs)^3 - 0.083751*log(n_obs)^2 - 0.31082*log(n_obs) - 1.5861
sigma_sw <- exp(0.0030302*log(n_obs)^2 - 0.082676*log(n_obs) - 0.4803)
Wn_stat <- (log(1 - W_manual) - mu_sw) / sigma_sw

# Gráfica Shapiro-Wilk (Distribución Normal Estándar)
shadeDist(Wn_stat, ddist = "dnorm", lower.tail = FALSE, 
          sub = paste("Wn:", round(Wn_stat, 3), "W:", round(W_manual, 4)))

# Validación automática
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.95315, p-value = 0.3166

Interpretación de la Prueba Jarque-Bera (JB) La prueba Jarque-Bera se utilizó para evaluar si la asimetría y la curtosis de los errores del modelo se desvían de los valores esperados en una distribución normal. El estadístico calculado fue de 1.4004 con un p-value de 0.4965. Dado que este valor de probabilidad es considerablemente mayor al nivel de significancia de 0.05, no existe evidencia estadística para rechazar la hipótesis nula de normalidad. Desde un enfoque gráfico, el estadístico se ubica dentro de la región de aceptación de la distribución Chi-cuadrado con 2 grados de libertad, sugiriendo que la forma de la distribución de los residuos es consistente con una campana de Gauss.

Interpretación de la Prueba Lilliefors (KS) A diferencia de la prueba Kolmogorov-Smirnov convencional, la prueba de Lilliefors se aplicó para corregir la estimación de los parámetros a partir de la muestra de 24 meses. El resultado arrojó un estadístico D de 0.13659 y un p-value de 0.2935. Al ser el p-value mayor a 0.05, se determina que no es posible rechazar la hipótesis nula. Esto indica que la discrepancia máxima entre la función de distribución acumulada de los residuos y la distribución normal teórica no es estadísticamente significativa, confirmando que los errores se comportan de manera aleatoria y normal.

Interpretación de la Prueba Shapiro-Wilk (SW)Se ejecutó la prueba de Shapiro-Wilk debido a su alta sensibilidad y eficacia en muestras pequeñas. El estadístico obtenido fue W = 0.95315, con un p-value asociado de 0.3166. Al contrastar este resultado con un nivel de significancia del 5%, se concluye que no hay fundamentos para rechazar la normalidad de los datos. En términos gráficos, al observar la distribución normal del estadístico normalizado Wn, se percibe que este cae dentro del cuerpo central de la distribución y no en las colas de rechazo, lo que ratifica la linealidad y el buen comportamiento de los términos de error.

# Evaluación de Multicolinealidad
# ---

# 1. Carga de librerías necesarias
library(stargazer)
library(mctest)
library(psych)
library(car)
library(fastGraph)
# --- 1. Índice de Condición ---
X_mat <- model.matrix(modelo_ventas) # Matriz de diseño
XX_matrix <- t(X_mat) %*% X_mat

# Normalización
options(scipen = 999)
Sn <- solve(diag(sqrt(diag(XX_matrix))))
XX_norm <- (Sn %*% XX_matrix) %*% Sn

# Autovalores y Cálculo de K
lambdas <- eigen(XX_norm, symmetric = TRUE)
K <- sqrt(max(lambdas$values) / min(lambdas$values))
print(paste("Índice de Condición (K):", round(K, 4)))
## [1] "Índice de Condición (K): 71.1635"
# --- 2. Prueba de Farrar-Glauber (Cálculo Manual) ---
Zn <- scale(X_mat[,-1]) # Escalamiento de regresores
n <- nrow(Zn)
m <- ncol(Zn)
R <- (t(Zn) %*% Zn) * (1 / (n - 1)) # Matriz de Correlación R
det_R <- det(R)

# Estadístico Chi-cuadrado FG
chi_FG <- -(n - 1 - (2 * m + 5) / 6) * log(det_R)
gl_fg <- m * (m - 1) / 2
VC_fg <- qchisq(p = 0.95, df = gl_fg)
p_val_fg <- pchisq(chi_FG, df = gl_fg, lower.tail = FALSE)

# Gráfica FG
shadeDist(chi_FG, ddist = "dchisq", parm1 = gl_fg, lower.tail = FALSE, 
          sub = paste("VC:", round(VC_fg, 2), " | FG:", round(chi_FG, 2)))

# --- 3. Factores Inflacionarios de la Varianza (VIF) ---
VIFs_car <- vif(modelo_ventas)
print(VIFs_car)
##        C        P        M 
## 7.631451 3.838911 9.449210

Se observa que el modelo presenta una situación de multicolinealidad estructural importante. En primera instancia, el Índice de Condición arrojó un valor de 71.1635, el cual, al superar ampliamente el umbral crítico de 30, clasifica la multicolinealidad como severa, indicando una fuerte interdependencia global entre los regresores. Esta conclusión se ve respaldada por la Prueba de Farrar-Glauber, donde el estadístico obtenido se situó en la zona de rechazo de la distribución Chi-cuadrado, confirmando estadísticamente que las variables no son independientes entre sí. No obstante, al analizar los Factores Inflacionarios de la Varianza (VIF), se observa que los valores para Comercialización (7.63), Personal (3.84) y Materias Primas (9.45) se mantienen por debajo de 10, lo que sugiere que, aunque existe una colinealidad grave a nivel de conjunto que podría inestabilizar el modelo, la precisión individual de los coeficientes aún no ha sido degradada al grado de invalidar las estimaciones. En conclusión, el modelo posee una alta capacidad explicativa, pero se debe tener cautela al interpretar el efecto aislado de cada variable debido a la estrecha relación detectada entre los gastos operativos y el nivel de ventas.

# VERIFICACIÓN DEL SUPUESTO DE HOMOCEDASTICIDAD
# ==========================================================
library(lmtest)
library(stargazer)
library(fastGraph)
# --- 1. Análisis Gráfico de los Residuos ---
# Graficamos residuos frente a valores ajustados
plot(modelo_ventas$fitted.values, modelo_ventas$residuals,
     main = "Gráfico de Dispersión de Residuos",
     xlab = "Valores Ajustados",
     ylab = "Residuos",
     pch = 19, col = "steelblue")
abline(h = 0, col = "red", lwd = 2)

# --- 2. Prueba de Breusch-Pagan (BP) ---
# H0: Varianza constante (Homocedasticidad)
# H1: Varianza no constante (Heterocedasticidad)
prueba_bp <- bptest(modelo_ventas)
print(prueba_bp)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ventas
## BP = 0.087869, df = 3, p-value = 0.9933
# Gráfica de la prueba BP
BP_stat <- prueba_bp$statistic
gl_bp <- prueba_bp$parameter
VC_bp <- qchisq(0.95, gl_bp)

shadeDist(BP_stat, ddist = "dchisq", parm1 = gl_bp, lower.tail = FALSE,
          sub = paste("VC:", round(VC_bp, 2), " | BP:", round(BP_stat, 2)))

# --- 3. Prueba de White ---
# En tu clase usan bptest con términos cuadráticos y cruzados
prueba_white <- bptest(modelo_ventas, ~ C*P + C*M + P*M + I(C^2) + I(P^2) + I(M^2), data = ventas_empresa)
print(prueba_white)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ventas
## BP = 7.1227, df = 9, p-value = 0.6244

Para verificar el supuesto de varianza constante en los errores, se realizó inicialmente un análisis gráfico donde la dispersión de los residuos frente a los valores ajustados no mostró patrones sistemáticos o formas de “embudo”, sugiriendo preliminarmente la ausencia de heterocedasticidad. Este hallazgo fue validado mediante la Prueba de Breusch-Pagan, la cual arrojó un estadístico de 0.0878 y un p-value de 0.9933. Al ser este valor de probabilidad significativamente superior al nivel de significancia de 0.05, no se rechaza la hipótesis nula , confirmando que la varianza de los residuos es constante. Complementariamente, se aplicó la Prueba de White para detectar formas más complejas de heterocedasticidad, obteniendo un p-value de 0.6244, el cual refuerza la decisión de no rechazar la homocedasticidad. Gráficamente, el estadístico de la prueba se sitúa dentro de la zona de aceptación de la distribución Chi-cuadrado. En conclusión, los errores del modelo presentan un comportamiento homocedástico, lo que asegura que los estimadores son eficientes y que las pruebas de inferencia estadística realizadas sobre las ventas son válidas.

# ==========================================================
# VERIFICACIÓN DE NO AUTOCORRELACIÓN (BREUSCH-GODFREY)
# ==========================================================
library(lmtest)
library(fastGraph)
# 1. Prueba de Breusch-Godfrey: Primer Orden [AR(1)]
# H0: No existe autocorrelación de primer orden
bg_1 <- bgtest(modelo_ventas, order = 1)
print(bg_1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  modelo_ventas
## LM test = 2.5963, df = 1, p-value = 0.1071
# 2. Prueba de Breusch-Godfrey: Segundo Orden [AR(2)]
# H0: No existe autocorrelación de segundo orden
bg_2 <- bgtest(modelo_ventas, order = 2)
print(bg_2)
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  modelo_ventas
## LM test = 3.8409, df = 2, p-value = 0.1465
# --- REPRESENTACIÓN GRÁFICA (Primer Orden) ---
LM_stat1 <- bg_1$statistic
gl1 <- bg_1$parameter
p_val1 <- bg_1$p.value
VC1 <- qchisq(0.95, gl1)

shadeDist(LM_stat1, ddist = "dchisq", parm1 = gl1, lower.tail = FALSE,
          sub = paste("Orden 1 - VC:", round(VC1, 2), " | LM:", round(LM_stat1, 2)))

# --- REPRESENTACIÓN GRÁFICA (Segundo Orden) ---
LM_stat2 <- bg_2$statistic
gl2 <- bg_2$parameter
p_val2 <- bg_2$p.value
VC2 <- qchisq(0.95, gl2)

shadeDist(LM_stat2, ddist = "dchisq", parm1 = gl2, lower.tail = FALSE,
          sub = paste("Orden 2 - VC:", round(VC2, 2), " | LM:", round(LM_stat2, 2)))

Al evaluar el supuesto de independencia en los residuos mediante la Prueba del Multiplicador de Lagrange (Breusch-Godfrey), se analizaron los patrones de dependencia para el primero y segundo orden. Para la autocorrelación de primer orden, se obtuvo un estadístico LM de 2.5963 con un p-value de 0.1071, mientras que para el segundo orden, el estadístico fue de 3.8409 con un p-value de 0.1465. Dado que en ambos casos los valores de probabilidad resultaron mayores al nivel de significancia de 0.05, no se rechaza la hipótesis nula de ausencia de autocorrelación. Gráficamente, se observa que los estadísticos calculados para ambos órdenes se sitúan dentro de la zona de aceptación de la distribución Chi-cuadrado, lo que indica que no existe evidencia de una estructura de dependencia temporal en los errores. Es por ello que, el modelo cumple satisfactoriamente con el supuesto de no autocorrelación, confirmando que los residuos son independientes y que las estimaciones de los errores estándar y las pruebas de significancia resultan consistentes y válidas.