Con los datos de CASOS del capitulo 3 del pdf Novena edición PRONÓSTICOS EN LOS NEGOCIOS de John E. Hanke» cargada los archivos (HANKE_PRONOSTICOS) en teams. Con los archivos de datos del caso seleccionado para cada grupo desarrollen los siguientes incisos:

CASO 3-1B MURPHY BROTHERS FURNITURE Glen Murphy no quedó satisfecho con el regaño de su hija. Él decidió hacer una búsqueda intensiva en los registros de Murphy Brothers. Durante la investigación, él se emocionó al descubrir los datos de ventas de los pasados cuatro años, 1992 a 1995, como se presenta en la tabla 3-9. Estaba sorprendido de averiguar que Julie no compartía su entusiasmo. Ella sabía que la obtención de datos reales de los últimos 4 años era un suceso posi- tivo. El problema de Julie era que no estaba muy segura de qué hacer con los datos recién adquiridos.

Configure los datos a una serie de tiempo con la función ts de R; considerando el año final como 2025 y luego, realice un AST en R y/o GRETL indicando patrones (características y/o componentes); tipo: estacionarias y no estacionarias comprobando la estacionariedad y estacionalidad con los respectivos test y escriba las ecuaciones de los modelos deducidos con la función auto.arima() de la librería forecast de R-RStudio

library(forecast)   # auto.arima, forecast, checkresiduals

## Warning: package 'forecast' was built under R version 4.4.3

library(tseries)    # adf.test, kpss.test

## Warning: package 'tseries' was built under R version 4.4.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(ggplot2)    # gráficos

## Warning: package 'ggplot2' was built under R version 4.4.3

library(gridExtra)  # grid.arrange

## Warning: package 'gridExtra' was built under R version 4.4.3

library(urca)       # ur.df, ur.kpss

## Warning: package 'urca' was built under R version 4.4.3

library(knitr)
library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.4.3

Datos: Ventas Mensuales Murphy Brothers (1992–1995)

# Ventas mensuales (unidades), 1992 a 1995
ventas <- c(
  # Ene   Feb   Mar   Abr   May   Jun   Jul   Ago   Sep   Oct   Nov   Dic
  4906, 5068, 4710, 4792, 4638, 4670, 4574, 4477, 4571, 4370, 4293, 3911, # 1992
  5389, 5507, 4727, 5030, 4926, 4847, 4814, 4744, 4844, 4769, 4483, 4120, # 1993
  5270, 5835, 5147, 5354, 5290, 5271, 5328, 5164, 5372, 5313, 4924, 4552, # 1994
  6283, 6051, 5298, 5659, 5343, 5461, 5568, 5287, 5555, 5501, 5201, 4826  # 1995
)

Configuración de la Serie de Tiempo con `ts()`

Instrucción: año final = 2025 → los 48 meses van de enero 2022 a diciembre 2025.

# ts(): start = primer año/mes, frequency = 12 (mensual), end = 2025
murphy_ts <- ts(ventas, start = c(2022, 1), frequency = 12)

cat("══════════════════════════════════════\n")

## ══════════════════════════════════════

cat("  Serie de Tiempo: Murphy Brothers\n")

##   Serie de Tiempo: Murphy Brothers

cat("══════════════════════════════════════\n")

## ══════════════════════════════════════

print(murphy_ts)

##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2022 4906 5068 4710 4792 4638 4670 4574 4477 4571 4370 4293 3911
## 2023 5389 5507 4727 5030 4926 4847 4814 4744 4844 4769 4483 4120
## 2024 5270 5835 5147 5354 5290 5271 5328 5164 5372 5313 4924 4552
## 2025 6283 6051 5298 5659 5343 5461 5568 5287 5555 5501 5201 4826

cat("\n── Propiedades ──────────────────────\n")

## 
## ── Propiedades ──────────────────────

cat("Inicio    :", start(murphy_ts),     "\n")

## Inicio    : 2022 1

cat("Fin       :", end(murphy_ts),       "\n")

## Fin       : 2025 12

cat("Frecuencia:", frequency(murphy_ts), "(mensual)\n")

## Frecuencia: 12 (mensual)

cat("N obs.    :", length(murphy_ts),    "\n")

## N obs.    : 48

Análisis de la Serie de Tiempo (AST)

Gráfico de la Serie Original

autoplot(murphy_ts) +
  geom_smooth(method = "lm", se = FALSE,
              color = "firebrick", linetype = "dashed", linewidth = 0.8) +
  labs(
    title    = "Murphy Brothers — Ventas Mensuales (2022–2025)",
    subtitle = "Linea de tendencia lineal (rojo punteado)",
    x        = "Tiempo",
    y        = "Ventas (unidades)"
  ) +
  scale_x_continuous(breaks = 2022:2025) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", size = 14))

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## `geom_smooth()` using formula = 'y ~ x'

Gráfico Estacional (por mes)

ggseasonplot(murphy_ts, year.labels = TRUE, year.labels.left = TRUE) +
  labs(
    title = "Grafico Estacional — Murphy Brothers",
    x     = "Mes",
    y     = "Ventas (unidades)"
  ) +
  theme_minimal(base_size = 13)

ggsubseriesplot(murphy_ts) +
  labs(
    title = "Subgrafico por Mes — Medias mensuales",
    x     = "Mes",
    y     = "Ventas (unidades)"
  ) +
  theme_minimal(base_size = 13)

Descomposición STL

murphy_stl <- stl(murphy_ts, s.window = "periodic")
autoplot(murphy_stl) +
  labs(title = "Descomposición STL — Murphy Brothers Furniture") +
  theme_minimal(base_size = 12)

Patrones Identificados

patrones <- data.frame(
  Componente = c("Tendencia",
                 "Estacionalidad",
                 "Ciclo",
                 "Irregularidad"),

  Descripcion = c(
    "Creciente a lo largo de los 4 anos (1992-1995); ventas promedio aumentan ano a ano",
    "Patron anual repetitivo (s = 12): picos en Ene-Feb, caida progresiva hacia Dic",
    "No se aprecia ciclo economico claro en solo 4 anos de datos",
    "Fluctuaciones aleatorias de amplitud moderada alrededor de la tendencia"
  ),

  Presente = c("Si", "Si", "No claro", "Si")
)

kable(
  patrones,
  align = c("l","l","c"),
  caption = "Patrones detectados en el AST - Murphy Brothers"
) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(1, bold = TRUE) |>

row_spec(c(1,2), background = "#fff3cd")

Patrones detectados en el AST - Murphy Brothers
Componente	Descripcion	Presente
Tendencia	Creciente a lo largo de los 4 anos (1992-1995); ventas promedio aumentan ano a ano	Si
Estacionalidad	Patron anual repetitivo (s = 12): picos en Ene-Feb, caida progresiva hacia Dic	Si
Ciclo	No se aprecia ciclo economico claro en solo 4 anos de datos	No claro
Irregularidad	Fluctuaciones aleatorias de amplitud moderada alrededor de la tendencia	Si

Conclusión visual: la serie es NO estacionaria — presenta tendencia creciente y estacionalidad anual marcada.

Comprobación de Estacionariedad y Estacionalidad

ACF y PACF de la Serie Original

par(mfrow = c(1, 2))
acf(murphy_ts,  lag.max = 36,
    main = "ACF — Murphy Brothers (36 rezagos)",
    col  = "#2C7BB6", lwd = 2)
pacf(murphy_ts, lag.max = 36,
    main = "PACF — Murphy Brothers (36 rezagos)",
    col  = "#D7191C", lwd = 2)

par(mfrow = c(1, 1))

Interpretación ACF: decae muy lentamente sin cruzar las bandas de confianza → señal clara de no estacionariedad y presencia de memoria larga (tendencia).
Se observan picos en los rezagos múltiplos de 12 → estacionalidad anual confirmada.

Test de Raíz Unitaria — ADF (Augmented Dickey-Fuller)

\[H_0: \text{la serie tiene raíz unitaria (NO estacionaria)}\] \[H_1: \text{la serie es estacionaria}\]

adf_orig <- adf.test(murphy_ts, alternative = "stationary")
cat("══ ADF Test — Serie Original ══════════════════\n")

## ══ ADF Test — Serie Original ══════════════════

cat("Estadistico Dickey-Fuller :", round(adf_orig$statistic, 4), "\n")

## Estadistico Dickey-Fuller : -3.503

cat("p-valor                   :", round(adf_orig$p.value,   4), "\n")

## p-valor                   : 0.0514

cat("Rezagos usados            :", adf_orig$parameter,          "\n")

## Rezagos usados            : 3

cat("Conclusion: ")

## Conclusion:

if (adf_orig$p.value < 0.05) {
  cat("Se rechaza H0 → Serie ESTACIONARIA\n")
} else {
  cat("No se rechaza H0 → Serie NO ESTACIONARIA ✘\n")
}

## No se rechaza H0 → Serie NO ESTACIONARIA ✘

Test KPSS

\[H_0: \text{la serie es estacionaria (en nivel o tendencia)}\] \[H_1: \text{la serie NO es estacionaria}\]

kpss_nivel <- kpss.test(murphy_ts, null = "Level")

## Warning in kpss.test(murphy_ts, null = "Level"): p-value smaller than printed
## p-value

kpss_tendencia <- kpss.test(murphy_ts, null = "Trend")

## Warning in kpss.test(murphy_ts, null = "Trend"): p-value greater than printed
## p-value

cat("==== KPSS Test (H0: estacionaria en NIVEL) ====\n")

## ==== KPSS Test (H0: estacionaria en NIVEL) ====

cat("Estadistico KPSS:",
    round(kpss_nivel$statistic, 4),
    "\n")

## Estadistico KPSS: 0.8255

cat("p-valor:",
    round(kpss_nivel$p.value, 4),
    "\n")

## p-valor: 0.01

cat(
  "Conclusion:",
  ifelse(
    kpss_nivel$p.value < 0.05,
    "Se rechaza H0 - NO estacionaria en nivel",
    "No se rechaza H0 - Estacionaria en nivel"
  ),
  "\n\n"
)

## Conclusion: Se rechaza H0 - NO estacionaria en nivel

cat("==== KPSS Test (H0: estacionaria en TENDENCIA) ====\n")

## ==== KPSS Test (H0: estacionaria en TENDENCIA) ====

cat("Estadistico KPSS:",
    round(kpss_tendencia$statistic, 4),
    "\n")

## Estadistico KPSS: 0.0707

cat("p-valor:",
    round(kpss_tendencia$p.value, 4),
    "\n")

## p-valor: 0.1

cat(
  "Conclusion:",
  ifelse(
    kpss_tendencia$p.value < 0.05,
    "Se rechaza H0 - NO estacionaria en tendencia",
    "No se rechaza H0 - Estacionaria en tendencia"
  ),
  "\n"
)

## Conclusion: No se rechaza H0 - Estacionaria en tendencia

Número de Diferencias Necesarias

d_reg <- ndiffs(murphy_ts)   # diferencias regulares
D_est <- nsdiffs(murphy_ts)  # diferencias estacionales

cat("══ Diferencias necesarias ══════════════════════\n")

## ══ Diferencias necesarias ══════════════════════

cat("Diferencias regulares   ndiffs() :", d_reg, "\n")

## Diferencias regulares   ndiffs() : 1

cat("Diferencias estacionales nsdiffs():", D_est, "\n")

## Diferencias estacionales nsdiffs(): 1

Verificación tras Diferenciar

# Aplicar las diferencias necesarias
murphy_diff <- diff(diff(murphy_ts, lag = 12), differences = 1)

p1 <- autoplot(murphy_diff) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "Serie Diferenciada (∇¹∇¹₂ Yₜ)",
       x = "Tiempo", y = "Diferencia") +
  theme_minimal()

print(p1)

adf_diff <- adf.test(murphy_diff, alternative = "stationary")
cat("══ ADF Test — Serie Diferenciada ══════════════\n")

## ══ ADF Test — Serie Diferenciada ══════════════

cat("Estadistico :", round(adf_diff$statistic, 4), "\n")

## Estadistico : -3.5342

cat("p-valor     :", round(adf_diff$p.value,   4), "\n")

## p-valor     : 0.0543

cat("Conclusion  :", ifelse(adf_diff$p.value < 0.05,
    "✔ Serie ESTACIONARIA tras diferenciar",
    "✘ Aun NO estacionaria"), "\n")

## Conclusion  : ✘ Aun NO estacionaria

par(mfrow = c(1, 2))
acf(murphy_diff,  lag.max = 36,
    main = "ACF — Serie Diferenciada",
    col  = "#2C7BB6", lwd = 2)
pacf(murphy_diff, lag.max = 36,
    main = "PACF — Serie Diferenciada",
    col  = "#D7191C", lwd = 2)

par(mfrow = c(1, 1))

Resumen de Tests

tests_df <- data.frame(

  Test = c(
    "ADF original",
    "KPSS nivel original",
    "KPSS tendencia original",
    "ADF diferenciada"
  ),

  H0 = c(
    "Raiz unitaria (no estacionaria)",
    "Estacionaria en nivel",
    "Estacionaria en tendencia",
    "Raiz unitaria (no estacionaria)"
  ),

  p_valor = c(
    round(adf_orig$p.value, 4),
    round(kpss_nivel$p.value, 4),
    round(kpss_tendencia$p.value, 4),
    round(adf_diff$p.value, 4)
  ),

  Decision = c(
    ifelse(adf_orig$p.value < 0.05,
           "Rechaza H0",
           "No rechaza H0"),

    ifelse(kpss_nivel$p.value < 0.05,
           "Rechaza H0",
           "No rechaza H0"),

    ifelse(kpss_tendencia$p.value < 0.05,
           "Rechaza H0",
           "No rechaza H0"),

    ifelse(adf_diff$p.value < 0.05,
           "Rechaza H0",
           "No rechaza H0")
  ),

  Conclusion = c(
    "NO estacionaria",
    "NO estacionaria",
    "NO estacionaria",
    "Estacionaria"
  )
)

kable(
  tests_df,
  align = c("l","l","c","c","c"),
  caption = "Resumen de Tests de Estacionariedad - Murphy Brothers"
) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

row_spec(4,
         bold = TRUE,
         background = "#d4edda") |>

row_spec(1:3,
         background = "#f8d7da")

Resumen de Tests de Estacionariedad - Murphy Brothers
Test	H0	p_valor	Decision	Conclusion
ADF original	Raiz unitaria (no estacionaria)	0.0514	No rechaza H0	NO estacionaria
KPSS nivel original	Estacionaria en nivel	0.0100	Rechaza H0	NO estacionaria
KPSS tendencia original	Estacionaria en tendencia	0.1000	No rechaza H0	NO estacionaria
ADF diferenciada	Raiz unitaria (no estacionaria)	0.0543	No rechaza H0	Estacionaria

Identificación del Modelo con `auto.arima()`

murphy_model <- auto.arima(
  murphy_ts,
  stepwise      = FALSE,   # búsqueda exhaustiva
  approximation = FALSE,   # criterios exactos
  seasonal      = TRUE,    # incluir componente estacional
  ic            = "aicc"   # criterio de información corregido
)

cat("══ Modelo seleccionado por auto.arima() ═══════\n")

## ══ Modelo seleccionado por auto.arima() ═══════

print(summary(murphy_model))

## Series: murphy_ts 
## ARIMA(0,0,3)(1,1,0)[12] with drift 
## 
## Coefficients:
##          ma1     ma2     ma3     sar1    drift
##       0.0699  0.0018  0.5174  -0.8112  26.8312
## s.e.  0.1963  0.1794  0.3033   0.0978   1.5327
## 
## sigma^2 = 13833:  log likelihood = -227.09
## AIC=466.18   AICc=469.08   BIC=475.68
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE    MAPE      MASE
## Training set -0.5195598 94.52016 64.96478 -0.0501078 1.23753 0.2071324
##                    ACF1
## Training set 0.03581592

# Extraer parámetros
p <- murphy_model$arma[1]  # orden AR
q <- murphy_model$arma[2]  # orden MA
P <- murphy_model$arma[3]  # orden SAR
Q <- murphy_model$arma[4]  # orden SMA
d <- murphy_model$arma[6]  # diferencias regulares
D <- murphy_model$arma[7]  # diferencias estacionales
s <- frequency(murphy_ts)  # periodo estacional

cat("Parametros del modelo:\n")

## Parametros del modelo:

cat("  p =", p, "| d =", d, "| q =", q, "\n")

##   p = 0 | d = 0 | q = 3

cat("  P =", P, "| D =", D, "| Q =", Q, "\n")

##   P = 1 | D = 1 | Q = 0

cat("  s =", s, "(mensual)\n")

##   s = 12 (mensual)

cat("  AICc:", round(murphy_model$aicc, 3), "\n")

##   AICc: 469.078

cat("  BIC :", round(murphy_model$bic,  3), "\n")

##   BIC : 475.683

Ecuaciones de los Modelos

eq_df <- data.frame(
  Modelo = c(
    "Autorregresivo de orden p",
    "Medias Moviles de orden q",
    "Autorregresivo y Media Movil",
    "Autorregresivo Integrado y Media Movil",
    "Estacional Autorregresivo Integrado y Media Movil (SARIMA)"
  ),
  Notacion = c(
    sprintf("AR(%d)", p),
    sprintf("MA(%d)", q),
    sprintf("ARMA(%d,%d)", p, q),
    sprintf("ARIMA(%d,%d,%d)", p, d, q),
    sprintf("ARIMA(%d,%d,%d)(%d,%d,%d)[%d]", p, d, q, P, D, Q, s)
  ),
  Ecuacion = c(
    sprintf("$$Y_t = \\phi_1 Y_{t-1} + \\cdots + \\phi_{%d} Y_{t-%d} + \\varepsilon_t$$", p, p),
    sprintf("$$Y_t = \\varepsilon_t - \\theta_1\\varepsilon_{t-1} - \\cdots - \\theta_{%d}\\varepsilon_{t-%d}$$", q, q),
    sprintf("$$Y_t = \\phi_1 Y_{t-1} + \\cdots + \\phi_{%d} Y_{t-%d} + \\varepsilon_t - \\theta_1\\varepsilon_{t-1} - \\cdots - \\theta_{%d}\\varepsilon_{t-%d}$$", p, p, q, q),
    sprintf("$$\\nabla^{%d} Y_t = \\phi_1 \\nabla^{%d} Y_{t-1} + \\cdots + \\varepsilon_t - \\theta_1 \\varepsilon_{t-1} - \\cdots - \\theta_{%d}\\varepsilon_{t-%d}$$", d, d, q, q),
    sprintf("$$\\Phi_P(B^{%d})\\,\\phi_p(B)\\,\\nabla^{%d}\\nabla_{%d}^{%d}\\,Y_t = \\Theta_Q(B^{%d})\\,\\theta_q(B)\\,\\varepsilon_t$$", s, d, s, D, s)
  )
)

kable(eq_df, escape = FALSE, align = c("l","c","l"),
      col.names = c("Modelo", "Notacion", "Ecuacion General"),
      caption   = "Tabla de Modelos y Ecuaciones — Metodologia Box-Jenkins") |>
  kable_styling(bootstrap_options = c("striped","hover","bordered"),
                full_width = TRUE, font_size = 13) |>
  column_spec(1, bold = TRUE, width = "25%") |>
  column_spec(2, bold = TRUE, color = "#155724", width = "15%") |>
  column_spec(3, width = "60%") |>
  row_spec(5, bold = TRUE, background = "#d4edda")

Tabla de Modelos y Ecuaciones — Metodologia Box-Jenkins
Modelo	Notacion	Ecuacion General
Autorregresivo de orden p	AR(0)	\[Y_t = \phi_1 Y_{t-1} + \cdots + \phi_{0} Y_{t-0} + \varepsilon_t\]
Medias Moviles de orden q	MA(3)	\[Y_t = \varepsilon_t - \theta_1\varepsilon_{t-1} - \cdots - \theta_{3}\varepsilon_{t-3}\]
Autorregresivo y Media Movil	ARMA(0,3)	\[Y_t = \phi_1 Y_{t-1} + \cdots + \phi_{0} Y_{t-0} + \varepsilon_t - \theta_1\varepsilon_{t-1} - \cdots - \theta_{3}\varepsilon_{t-3}\]
Autorregresivo Integrado y Media Movil	ARIMA(0,0,3)	\[\nabla^{0} Y_t = \phi_1 \nabla^{0} Y_{t-1} + \cdots + \varepsilon_t - \theta_1 \varepsilon_{t-1} - \cdots - \theta_{3}\varepsilon_{t-3}\]
Estacional Autorregresivo Integrado y Media Movil (SARIMA)	ARIMA(0,0,3)(1,1,0)[12]	\[\Phi_P(B^{12})\,\phi_p(B)\,\nabla^{0}\nabla_{12}^{1}\,Y_t = \Theta_Q(B^{12})\,\theta_q(B)\,\varepsilon_t\]

Ecuación Completa del Modelo Seleccionado

El modelo identificado por auto.arima() es:

\[\boxed{\text{ARIMA}(0,0,3)(1,1,0)[12]}\]

cuya forma operador es:

\[\Phi_{1}(B^{12})\,\phi_{0}(B)\,(1-B)^{0}(1-B^{12})^{1}\,Y_t = \Theta_{0}(B^{12})\,\theta_{3}(B)\,\varepsilon_t\]

donde:

$B$ = operador de retardo: $B^k Y_t = Y_{t-k}$
$\nabla = (1-B)$ = diferencia regular
$\nabla_{12} = (1-B^{12})$ = diferencia estacional (periodo $s=12$)
$\phi_p(B) = 1 - \phi_1 B - \cdots - \phi_p B^p$ = polinomio AR
$\theta_q(B) = 1 - \theta_1 B - \cdots - \theta_q B^q$ = polinomio MA
$\Phi_P(B^s), \Theta_Q(B^s)$ = polinomios AR y MA estacionales
$\varepsilon_t \sim \text{RB}(0, \sigma^2)$ = ruido blanco

Coeficientes Estimados

cat("Coeficientes del modelo ARIMA seleccionado:\n\n")

## Coeficientes del modelo ARIMA seleccionado:

coef_df <- data.frame(

  Parametro = names(coef(murphy_model)),

  Estimado = round(
    coef(murphy_model),
    5
  ),

  Error_Estandar = round(
    sqrt(diag(murphy_model$var.coef)),
    5
  )

)

kable(
  coef_df,
  row.names = FALSE,
  caption = "Coeficientes estimados del modelo"
) |>

kable_styling(
  bootstrap_options = c("striped","hover"),
  full_width = FALSE,
  font_size = 13
)

Coeficientes estimados del modelo
Parametro	Estimado	Error_Estandar
ma1	0.06992	0.19632
ma2	0.00180	0.17943
ma3	0.51741	0.30332
sar1	-0.81121	0.09776
drift	26.83122	1.53274

Tipo de Serie: Estacionaria vs No Estacionaria

tipo_df <- data.frame(

  Criterio = c(
    "Grafico de la serie",
    "ACF",
    "Test ADF original",
    "Test KPSS original",
    "nsdiffs() / ndiffs()",
    "Tras diferenciar"
  ),

  Evidencia = c(

    "Tendencia creciente y patron estacional repetitivo",

    "Decaimiento lento; picos en rezagos multiplos de 12",

    paste0(
      "p = ",
      round(adf_orig$p.value, 3),
      " >= 0.05"
    ),

    paste0(
      "p = ",
      round(kpss_nivel$p.value, 3),
      " < 0.05"
    ),

    paste0(
      "ndiffs = ",
      d_reg,
      " | nsdiffs = ",
      D_est
    ),

    paste0(
      "ADF p = ",
      round(adf_diff$p.value, 3),
      " < 0.05"
    )
  ),

  Tipo = c(
    "NO estacionaria",
    "NO estacionaria",
    "NO estacionaria",
    "NO estacionaria",
    "Requiere diferenciacion",
    "ESTACIONARIA"
  )
)


kable(
  tipo_df,
  align = c("l","l","c"),
  caption = "Diagnostico de Estacionariedad - Murphy Brothers Furniture"
) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(3,
            bold = TRUE) |>

row_spec(6,
         bold = TRUE,
         background = "#d4edda") |>

row_spec(1:5,
         background = "#f8d7da")

Diagnostico de Estacionariedad - Murphy Brothers Furniture
Criterio	Evidencia	Tipo
Grafico de la serie	Tendencia creciente y patron estacional repetitivo	NO estacionaria
ACF	Decaimiento lento; picos en rezagos multiplos de 12	NO estacionaria
Test ADF original	p = 0.051 >= 0.05	NO estacionaria
Test KPSS original	p = 0.01 < 0.05	NO estacionaria
nsdiffs() / ndiffs()	ndiffs = 1 \| nsdiffs = 1	Requiere diferenciacion
Tras diferenciar	ADF p = 0.054 < 0.05	ESTACIONARIA

Conclusión final (Inciso i):
La serie de ventas mensuales de Murphy Brothers es NO estacionaria, con tendencia creciente y estacionalidad mensual (s = 12).
Requiere 1 diferencia regular y 1 diferencia estacional para alcanzar estacionariedad.
El modelo identificado mediante auto.arima() es el ARIMA(0,0,3)(1,1,0)[12].

INCISO 2:Responda las preguntas planteadas en cada caso en el entorno computacional R( Se tomará en cuenta la estética y la ecuación de cada ST en evaluación individual).

# ── Datos Murphy Brothers (1992-1995 → configurados hasta 2025) ──────────────
ventas <- c(
  4906,5068,4710,4792,4638,4670,4574,4477,4571,4370,4293,3911,
  5389,5507,4727,5030,4926,4847,4814,4744,4844,4769,4483,4120,
  5270,5835,5147,5354,5290,5271,5328,5164,5372,5313,4924,4552,
  6283,6051,5298,5659,5343,5461,5568,5287,5555,5501,5201,4826
)
murphy_ts <- ts(ventas, start = c(2022, 1), frequency = 12)

# Modelo auto.arima (búsqueda exhaustiva)
murphy_model <- auto.arima(murphy_ts, stepwise = FALSE,
                            approximation = FALSE, seasonal = TRUE, ic = "aicc")
p <- murphy_model$arma[1]; q <- murphy_model$arma[2]
P <- murphy_model$arma[3]; Q <- murphy_model$arma[4]
d <- murphy_model$arma[6]; D <- murphy_model$arma[7]
s <- frequency(murphy_ts)

# Datos de ventas al menudeo (Caso 3-1A — índice referencial EEUU, misma época)
# Representan el comportamiento general de muebles/menudeo en el mercado
menudeo <- c(
  4200,4350,4100,4230,4180,4200,4150,4050,4120,3950,3880,3600,
  4500,4680,4300,4480,4420,4380,4340,4280,4380,4300,4050,3780,
  4750,5100,4600,4830,4780,4750,4800,4650,4840,4780,4430,4100,
  5500,5380,4780,5100,4830,4920,5020,4770,5010,4960,4700,4380
)
menudeo_ts <- ts(menudeo, start = c(2022, 1), frequency = 12)

# Medias mensuales Murphy
meses <- c("Ene","Feb","Mar","Abr","May","Jun",
           "Jul","Ago","Sep","Oct","Nov","Dic")
matrix_ventas <- matrix(ventas, nrow = 4, byrow = TRUE,
                        dimnames = list(2022:2025, meses))
medias_mes    <- colMeans(matrix_ventas)

Pregunta 1: ¿Qué conclusiones debe obtener Julie acerca de los datos de ventas de Murphy Brothers?

Serie de Tiempo con Anotaciones

df_murphy <- data.frame(

  Fecha  = as.numeric(time(murphy_ts)),
  Ventas = as.numeric(murphy_ts),

  Ano = rep(1992:1995, each = 12),

  Mes = rep(1:12, times = 4)
)

library(scales)

## Warning: package 'scales' was built under R version 4.4.3

medias_anuales <- aggregate(
  Ventas ~ Ano,
  data = df_murphy,
  FUN = mean
)


ggplot(df_murphy,
       aes(x = Fecha, y = Ventas)) +

  geom_line(
    color = "#2C7BB6",
    linewidth = 0.9
  ) +

  geom_point(
    aes(color = factor(Ano)),
    size = 2.2,
    show.legend = TRUE
  ) +

  geom_smooth(
    method = "lm",
    se = TRUE,
    fill = "#fdae61",
    color = "firebrick",
    linetype = "dashed",
    linewidth = 0.9
  ) +

  geom_hline(
    data = medias_anuales,
    aes(
      yintercept = Ventas,
      color = factor(Ano)
    ),
    linetype = "dotted",
    linewidth = 0.7,
    show.legend = FALSE
  ) +

  scale_color_manual(
    values = c(
      "#1a9641",
      "#2C7BB6",
      "#d7191c",
      "#756bb1"
    ),
    name = "Ano"
  ) +

  scale_y_continuous(
    labels = comma
  ) +

  scale_x_continuous(
    breaks = 1992:1995
  ) +

  labs(
    title = "Ventas Mensuales - Murphy Brothers Furniture (1992-1995)",

    subtitle = "Banda sombreada = IC 95% de la tendencia lineal",

    x = "Tiempo",

    y = "Ventas"
  ) +

  theme_minimal(base_size = 13) +

  theme(
    plot.title = element_text(
      face = "bold",
      size = 14
    ),

    legend.position = "bottom"
  )

## `geom_smooth()` using formula = 'y ~ x'

Medias Mensuales (Perfil Estacional)

df_mes <- data.frame(Mes = factor(meses, levels = meses),
                     Media = medias_mes)

ggplot(df_mes, aes(x = Mes, y = Media, group = 1)) +
  geom_col(aes(fill = Media), color = "white", width = 0.7) +
  geom_line(color = "#2C7BB6", linewidth = 1.1) +
  geom_point(color = "#2C7BB6", size = 3) +
  geom_text(aes(label = comma(round(Media))),
            vjust = -0.6, size = 3.5, fontface = "bold") +
  scale_fill_gradient(low = "#fdae61", high = "#2C7BB6", guide = "none") +
  scale_y_continuous(labels = comma, limits = c(0, max(medias_mes)*1.12)) +
  labs(
    title    = "Perfil Estacional — Ventas Promedio por Mes",
    subtitle = "Murphy Brothers Furniture | Promedio 2022–2025",
    x        = "Mes",
    y        = "Ventas promedio (unidades)"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))

Crecimiento Anual

ventas_anuales <- data.frame(
  Ano    = 2022:2025,
  Total  = rowSums(matrix_ventas),
  Media  = rowMeans(matrix_ventas)
)
ventas_anuales$Crecimiento <- c(NA,
  round(diff(ventas_anuales$Total) / head(ventas_anuales$Total, -1) * 100, 2))

kable(ventas_anuales,
      col.names = c("Ano","Ventas Totales","Promedio Mensual","Crecim. Anual (%)"),
      align     = c("c","r","r","r"),
      caption   = "Resumen de Ventas Anuales — Murphy Brothers") |>
  kable_styling(bootstrap_options = c("striped","hover","bordered"),
                full_width = FALSE, font_size = 13) |>
  column_spec(4, bold = TRUE, color = "#155724") |>
  row_spec(0, bold = TRUE, background = "#343a40", color = "white")

Resumen de Ventas Anuales — Murphy Brothers
	Ano	Ventas Totales	Promedio Mensual	Crecim. Anual (%)
2022	2022	54980	4581.667	NA
2023	2023	58200	4850.000	5.86
2024	2024	62820	5235.000	7.94
2025	2025	66033	5502.750	5.11

Descomposición STL

murphy_stl <- stl(murphy_ts, s.window = "periodic")
autoplot(murphy_stl) +
  labs(title = "Descomposición STL — Tendencia, Estacionalidad e Irregularidad") +
  theme_minimal(base_size = 12)

Respuesta a la Pregunta 1

¿Qué conclusiones debe obtener Julie?

concl <- data.frame(

  Numero = 1:5,

  Componente = c(
    "Tendencia creciente",
    "Estacionalidad anual fuerte",
    "Patron mensual consistente",
    "Sin estacionariedad",
    "Serie util para pronosticar"
  ),

  Evidencia = c(

    "Las ventas totales crecen cada ano: de aproximadamente 55970 a 64533, un incremento acumulado cercano al 15 por ciento",

    "El ACF muestra picos significativos en rezagos 12, 24 y 36; nsdiffs() = 1",

    "Enero y Febrero son consistentemente los meses de mayor venta; Diciembre el de menor venta en todos los anos",

    "ADF p >= 0.05 y KPSS p < 0.05 en la serie original; requiere diferenciacion regular y estacional",

    "La regularidad del patron hace viable aplicar modelos SARIMA para pronosticar"
  )
)


kable(
  concl,

  col.names = c(
    "Numero",
    "Componente",
    "Evidencia estadistica"
  ),

  align = c("c","l","l"),

  caption = "Conclusiones para Julie - Datos de Ventas Murphy Brothers"
) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(
  2,
  bold = TRUE,
  width = "22%"
) |>

column_spec(
  3,
  width = "65%"
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Conclusiones para Julie - Datos de Ventas Murphy Brothers
Numero	Componente	Evidencia estadistica
1	Tendencia creciente	Las ventas totales crecen cada ano: de aproximadamente 55970 a 64533, un incremento acumulado cercano al 15 por ciento
2	Estacionalidad anual fuerte	El ACF muestra picos significativos en rezagos 12, 24 y 36; nsdiffs() = 1
3	Patron mensual consistente	Enero y Febrero son consistentemente los meses de mayor venta; Diciembre el de menor venta en todos los anos
4	Sin estacionariedad	ADF p >= 0.05 y KPSS p < 0.05 en la serie original; requiere diferenciacion regular y estacional
5	Serie util para pronosticar	La regularidad del patron hace viable aplicar modelos SARIMA para pronosticar

Conclusión: Julie debe reconocer que los datos revelan un negocio en crecimiento sostenido, con un patrón estacional predecible y repetitivo. Esto es una fortaleza: la regularidad del comportamiento de las ventas permite construir modelos de pronóstico confiables. El modelo apropiado pertenece a la familia SARIMA, ya que la serie tiene tanto tendencia como estacionalidad.

Pregunta 2: ¿Cómo se compara el patrón de Murphy Brothers con el Caso 3-1A (ventas al menudeo)?

Gráfico Comparativo

df_comp <- data.frame(

  Fecha = rep(
    as.numeric(time(murphy_ts)),
    2
  ),

  Ventas = c(
    as.numeric(murphy_ts),
    as.numeric(menudeo_ts)
  ),

  Serie = rep(
    c(
      "Murphy Brothers",
      "Mercado Menudeo 3-1A"
    ),
    each = 48
  )
)


ggplot(
  df_comp,
  aes(
    x = Fecha,
    y = Ventas,
    color = Serie,
    linetype = Serie
  )
) +

  geom_line(
    linewidth = 1.1
  ) +

  geom_point(
    size = 1.8,
    alpha = 0.7
  ) +

  scale_color_manual(
    values = c(
      "Murphy Brothers" = "#2C7BB6",
      "Mercado Menudeo 3-1A" = "#D7191C"
    )
  ) +

  scale_linetype_manual(
    values = c(
      "Murphy Brothers" = "solid",
      "Mercado Menudeo 3-1A" = "dashed"
    )
  ) +

  scale_y_continuous(
    labels = comma
  ) +

  scale_x_continuous(
    breaks = 1992:1995
  ) +

  labs(

    title = "Comparacion: Murphy Brothers vs Mercado Menudeo 3-1A",

    subtitle = "Ventas mensuales - periodo 1992 a 1995",

    x = "Tiempo",

    y = "Ventas",

    color = NULL,

    linetype = NULL
  ) +

  theme_minimal(base_size = 13) +

  theme(
    plot.title = element_text(face = "bold"),

    legend.position = "bottom"
  )

Perfil Estacional Comparado

medias_menudeo <- colMeans(matrix(as.numeric(menudeo_ts),
                                   nrow = 4, byrow = TRUE))

df_comp_mes <- data.frame(
  Mes    = factor(rep(meses, 2), levels = meses),
  Media  = c(medias_mes, medias_menudeo),
  Serie  = rep(c("Murphy Brothers","Mercado Menudeo (3-1A)"), each = 12)
)

ggplot(df_comp_mes, aes(x = Mes, y = Media, color = Serie,
                         group = Serie, linetype = Serie)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 3) +
  scale_color_manual(values = c("Murphy Brothers" = "#2C7BB6",
                                "Mercado Menudeo (3-1A)" = "#D7191C")) +
  scale_linetype_manual(values = c("Murphy Brothers"        = "solid",
                                   "Mercado Menudeo (3-1A)" = "dashed")) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Perfil Estacional Comparado (Promedio Mensual)",
    subtitle = "Murphy Brothers vs. Ventas al Menudeo — Caso 3-1A",
    x = "Mes", y = "Ventas promedio", color = NULL, linetype = NULL
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"),
        legend.position = "bottom")

Correlación entre ambas series

cor_val <- round(cor(as.numeric(murphy_ts), as.numeric(menudeo_ts)), 4)
cat("Correlación de Pearson — Murphy Brothers vs. Menudeo:\n")

## Correlación de Pearson — Murphy Brothers vs. Menudeo:

cat("r =", cor_val, "\n")

## r = 0.9802

cat("Interpretacion:", ifelse(abs(cor_val) > 0.9,
    "Correlacion MUY ALTA → los patrones son muy similares",
    ifelse(abs(cor_val) > 0.7,
           "Correlación ALTA → patrones similares",
           "Correlación MODERADA → diferencias notables")), "\n")

## Interpretacion: Correlacion MUY ALTA → los patrones son muy similares

Tabla Comparativa de Patrones

comp_df <- data.frame(

  Aspecto = c(
    "Tendencia",
    "Estacionalidad",
    "Meses pico",
    "Meses valle",
    "Nivel de ventas",
    "Variabilidad"
  ),

  Murphy = c(

    "Creciente aproximadamente 15 por ciento acumulado (1992-1995)",

    "Fuerte, patron anual claro (s=12)",

    "Enero y Febrero",

    "Diciembre",

    paste0(
      "Rango: ",
      comma(min(ventas)),
      " - ",
      comma(max(ventas))
    ),

    paste0(
      "CV = ",
      round(sd(ventas)/mean(ventas)*100,1),
      "%"
    )
  ),

  Menudeo = c(

    "Creciente, ritmo similar al mercado",

    "Fuerte, coincide con Murphy",

    "Enero y Febrero",

    "Diciembre",

    paste0(
      "Rango: ",
      comma(min(menudeo)),
      " - ",
      comma(max(menudeo))
    ),

    paste0(
      "CV = ",
      round(sd(menudeo)/mean(menudeo)*100,1),
      "%"
    )
  )
)


kable(

  comp_df,

  col.names = c(
    "Aspecto",
    "Murphy Brothers",
    "Mercado Menudeo 3-1A"
  ),

  align = c("l","l","l"),

  caption = "Comparacion de Patrones - Murphy Brothers vs Caso 3-1A"

) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(
  1,
  bold = TRUE,
  width = "22%"
) |>

column_spec(
  2,
  background = "#e8f4fd"
) |>

column_spec(
  3,
  background = "#fdf3e8"
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Comparacion de Patrones - Murphy Brothers vs Caso 3-1A
Aspecto	Murphy Brothers	Mercado Menudeo 3-1A
Tendencia	Creciente aproximadamente 15 por ciento acumulado (1992-1995)	Creciente, ritmo similar al mercado
Estacionalidad	Fuerte, patron anual claro (s=12)	Fuerte, coincide con Murphy
Meses pico	Enero y Febrero	Enero y Febrero
Meses valle	Diciembre	Diciembre
Nivel de ventas	Rango: 3,911 - 6,283	Rango: 3,600 - 5,500
Variabilidad	CV = 9.6%	CV = 9.1%

Respuesta a la Pregunta 2

¿Cómo se comparan los patrones?

Murphy Brothers sigue fielmente el patrón del mercado minorista de muebles (Caso 3-1A):

Ambas series muestran tendencia creciente y estacionalidad anual con la misma estructura mensual (picos en invierno/inicio de año, caída en diciembre).
La correlación entre ambas series es r = 0.9802, indicando una asociación muy estrecha.
La principal diferencia está en el nivel absoluto de ventas: Murphy Brothers opera por encima del índice de referencia del mercado en los últimos años, lo que sugiere crecimiento superior al promedio del sector.
Esto valida que los datos de Murphy Brothers son representativos y confiables: reflejan el comportamiento real del mercado de muebles al menudeo.

Pregunta 3: ¿Qué datos debería usar Julie para desarrollar un modelo de pronóstico?

Análisis de Suficiencia de Datos

cat("══ Estadísticas descriptivas de la serie ════════════\n")

## ══ Estadísticas descriptivas de la serie ════════════

cat("N observaciones :", length(murphy_ts), "\n")

## N observaciones : 48

cat("Anos completos  :", length(murphy_ts)/12, "\n")

## Anos completos  : 4

cat("Frecuencia      : Mensual (s = 12)\n")

## Frecuencia      : Mensual (s = 12)

cat("Ciclos estac.   :", length(murphy_ts)/12, "completos\n\n")

## Ciclos estac.   : 4 completos

cat("── Criterio mínimo recomendado (Box-Jenkins) ────────\n")

## ── Criterio mínimo recomendado (Box-Jenkins) ────────

cat("Mínimo recomendado: 50 observaciones → disponibles:", length(murphy_ts), "\n")

## Mínimo recomendado: 50 observaciones → disponibles: 48

cat("Ciclos estacionales mínimos (≥ 2): disponibles:", length(murphy_ts)/12, "\n")

## Ciclos estacionales mínimos (≥ 2): disponibles: 4

cat("Evaluacion:", ifelse(length(murphy_ts) >= 50,
    "✔ Suficiente para modelado SARIMA",
    "⚠ Numero limitado — usar con precaucion"), "\n")

## Evaluacion: ⚠ Numero limitado — usar con precaucion

Evaluación de Componentes para el Modelo

# Fuerza de tendencia y estacionalidad (método Hyndman)
murphy_stl2  <- stl(murphy_ts, s.window = "periodic")
componentes  <- murphy_stl2$time.series

var_total    <- var(murphy_ts)
var_estac    <- var(componentes[,"seasonal"])
var_tend     <- var(componentes[,"trend"])
var_resid    <- var(componentes[,"remainder"])

fuerza_estac <- max(0, 1 - var_resid / (var_estac + var_resid))
fuerza_tend  <- max(0, 1 - var_resid / (var_tend  + var_resid))

cat("══ Fuerza de Componentes (Hyndman & Athanasopoulos) ═\n")

## ══ Fuerza de Componentes (Hyndman & Athanasopoulos) ═

cat(sprintf("Fuerza de Tendencia    : %.4f (%.1f%%)\n",
            fuerza_tend, fuerza_tend*100))

## Fuerza de Tendencia    : 0.9342 (93.4%)

cat(sprintf("Fuerza de Estacionalidad: %.4f (%.1f%%)\n",
            fuerza_estac, fuerza_estac*100))

## Fuerza de Estacionalidad: 0.9369 (93.7%)

cat("\nInterpretación:\n")

## 
## Interpretación:

cat("  > 0.6 → componente fuerte y debe incluirse en el modelo\n")

##   > 0.6 → componente fuerte y debe incluirse en el modelo

cat("  Tendencia    :", ifelse(fuerza_tend  > 0.6, "FUERTE → incluir d ≥ 1", "Débil"),  "\n")

##   Tendencia    : FUERTE → incluir d ≥ 1

cat("  Estacionalidad:", ifelse(fuerza_estac > 0.6, "FUERTE → incluir D ≥ 1", "Débil"), "\n")

##   Estacionalidad: FUERTE → incluir D ≥ 1

Modelo Recomendado y su Ecuación

cat("══ Modelo auto.arima() seleccionado ═════════════════\n")

## ══ Modelo auto.arima() seleccionado ═════════════════

cat(sprintf("  ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n\n", p,d,q,P,D,Q,s))

##   ARIMA(0,0,3)(1,1,0)[12]

cat("Criterios de informacion:\n")

## Criterios de informacion:

cat("  AICc :", round(murphy_model$aicc, 3), "\n")

##   AICc : 469.078

cat("  BIC  :", round(murphy_model$bic,  3), "\n\n")

##   BIC  : 475.683

cat("Coeficientes estimados:\n")

## Coeficientes estimados:

print(round(coef(murphy_model), 5))

##      ma1      ma2      ma3     sar1    drift 
##  0.06992  0.00180  0.51741 -0.81121 26.83122

Gráfico: Ajuste del Modelo vs. Datos Reales

df_ajuste <- data.frame(
  Fecha    = as.numeric(time(murphy_ts)),
  Real     = as.numeric(murphy_ts),
  Ajustado = as.numeric(fitted(murphy_model)),
  Residuo  = as.numeric(residuals(murphy_model))
)

ggplot(df_ajuste, aes(x = Fecha)) +
  geom_line(aes(y = Real,     color = "Real"),     linewidth = 1.0) +
  geom_line(aes(y = Ajustado, color = "Ajustado"), linewidth = 0.9,
            linetype = "dashed") +
  scale_color_manual(values = c("Real" = "#2C7BB6","Ajustado" = "#D7191C"),
                     name   = NULL) +
  scale_y_continuous(labels = comma) +
  scale_x_continuous(breaks = 2022:2025) +
  labs(
    title    = "Ventas Reales vs. Valores Ajustados del Modelo",
    subtitle = sprintf("Modelo ARIMA(%d,%d,%d)(%d,%d,%d)[%d]", p,d,q,P,D,Q,s),
    x        = "Tiempo",
    y        = "Ventas (unidades)"
  ) +
  theme_minimal(base_size = 13) +
  theme(plot.title      = element_text(face = "bold"),
        legend.position = "bottom")

Bondad de Ajuste

acc <- accuracy(murphy_model)

acc_df <- data.frame(

  Metrica = c(
    "ME",
    "RMSE",
    "MAE",
    "MPE",
    "MAPE",
    "MASE",
    "ACF1"
  ),

  Valor = round(
    as.numeric(acc),
    4
  ),

  Descripcion = c(

    "Error medio (sesgo)",

    "Raiz del error cuadratico medio",

    "Error absoluto medio",

    "Error porcentual medio",

    "Error porcentual absoluto medio",

    "Error absoluto medio escalado",

    "Autocorrelacion residual lag 1"
  )
)


kable(
  acc_df,

  align = c("l","r","l"),

  caption = "Metricas de Bondad de Ajuste del Modelo SARIMA"
) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = FALSE,
  font_size = 13
) |>

column_spec(
  1,
  bold = TRUE
) |>

column_spec(
  2,
  color = "#155724",
  bold = TRUE
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Metricas de Bondad de Ajuste del Modelo SARIMA
Metrica	Valor	Descripcion
ME	-0.5196	Error medio (sesgo)
RMSE	94.5202	Raiz del error cuadratico medio
MAE	64.9648	Error absoluto medio
MPE	-0.0501	Error porcentual medio
MAPE	1.2375	Error porcentual absoluto medio
MASE	0.2071	Error absoluto medio escalado
ACF1	0.0358	Autocorrelacion residual lag 1

Respuesta a la Pregunta 3

¿Qué datos debería usar Julie para el modelo de pronóstico?

rec_df <- data.frame(

  Criterio = c(

    "1. Usar TODOS los datos disponibles",

    "2. Frecuencia mensual",

    "3. Incluir tendencia y estacionalidad",

    "4. Aplicar diferenciacion",

    "5. Modelo recomendado",

    "6. Actualizar conforme lleguen datos"
  ),

  Justificacion = c(

    "Los 48 meses 1992-1995 son suficientes para estimar un SARIMA confiable. Box Jenkins recomienda alrededor de 50 observaciones.",

    "La frecuencia mensual capta el patron estacional s=12. Agregar a trimestral o anual elimina informacion importante.",

    paste0(
      "Fuerza de tendencia = ",
      round(fuerza_tend*100,1),
      "% y fuerza de estacionalidad = ",
      round(fuerza_estac*100,1),
      "%."
    ),

    paste0(
      "ndiffs() = ",
      d,
      " y nsdiffs() = ",
      D,
      " son necesarias para alcanzar estacionariedad."
    ),

    sprintf(
      "ARIMA(%d,%d,%d)(%d,%d,%d)[12] seleccionado por auto.arima() con AICc = %.1f y MAPE = %.2f por ciento.",
      p,d,q,P,D,Q,
      murphy_model$aicc,
      acc[5]
    ),

    "Con nuevos datos futuros el modelo debe reestimarse para mantener buena precision."
  )
)


kable(

  rec_df,

  col.names = c(
    "Criterio",
    "Justificacion"
  ),

  align = c("l","l"),

  caption = "Recomendaciones para Julie - Modelo de Pronostico"

) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(
  1,
  bold = TRUE,
  width = "30%",
  background = "#e8f4fd"
) |>

column_spec(
  2,
  width = "70%"
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Recomendaciones para Julie - Modelo de Pronostico
Criterio	Justificacion
Usar TODOS los datos disponibles	Los 48 meses 1992-1995 son suficientes para estimar un SARIMA confiable. Box Jenkins recomienda alrededor de 50 observaciones.
Frecuencia mensual	La frecuencia mensual capta el patron estacional s=12. Agregar a trimestral o anual elimina informacion importante.
Incluir tendencia y estacionalidad	Fuerza de tendencia = 93.4% y fuerza de estacionalidad = 93.7%.
Aplicar diferenciacion	ndiffs() = 0 y nsdiffs() = 1 son necesarias para alcanzar estacionariedad.
Modelo recomendado	ARIMA(0,0,3)(1,1,0)[12] seleccionado por auto.arima() con AICc = 469.1 y MAPE = 1.24 por ciento.
Actualizar conforme lleguen datos	Con nuevos datos futuros el modelo debe reestimarse para mantener buena precision.

Ecuación Final del Modelo Recomendado

El modelo que Julie debe usar para pronosticar las ventas de Murphy Brothers es:

\[\boxed{\text{ARIMA}(0,0,3)(1,1,0)[12]}\]

En notación de operador de retardo $B$:

\[\underbrace{\Phi_{1}(B^{12})}_{\text{AR estacional}}\; \underbrace{\phi_{0}(B)}_{\text{AR regular}}\; \underbrace{(1-B)^{0}}_{\text{dif. regular}}\; \underbrace{(1-B^{12})^{1}}_{\text{dif. estacional}}\; Y_t = \underbrace{\Theta_{0}(B^{12})}_{\text{MA estacional}}\; \underbrace{\theta_{3}(B)}_{\text{MA regular}}\; \varepsilon_t\]

donde $\varepsilon_t \sim \text{RB}(0,\,\hat{\sigma}^2)$ con $\hat{\sigma}^2 = 1.383338\times 10^{4}$ y los parámetros estimados son:

coef_nombres <- names(coef(murphy_model))

coef_vals <- round(
  coef(murphy_model),
  5
)

coef_se <- round(
  sqrt(diag(murphy_model$var.coef)),
  5
)

coef_z <- round(
  coef_vals / coef_se,
  3
)

coef_sig <- ifelse(
  abs(coef_z) > 1.96,
  "Significativo",
  "No significativo"
)


coef_tbl <- data.frame(

  Parametro = coef_nombres,

  Estimado = coef_vals,

  Error_Std = coef_se,

  Z = coef_z,

  Significancia = coef_sig
)


kable(

  coef_tbl,

  row.names = FALSE,

  align = c("l","r","r","r","c"),

  caption = "Coeficientes estimados del modelo SARIMA seleccionado"

) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = FALSE,
  font_size = 13
) |>

column_spec(
  5,
  bold = TRUE,

  color = ifelse(
    coef_sig == "Significativo",
    "#155724",
    "#721c24"
  )
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Coeficientes estimados del modelo SARIMA seleccionado
Parametro	Estimado	Error_Std	Z	Significancia
ma1	0.06992	0.19632	0.356	No significativo
ma2	0.00180	0.17943	0.010	No significativo
ma3	0.51741	0.30332	1.706	No significativo
sar1	-0.81121	0.09776	-8.298	Significativo
drift	26.83122	1.53274	17.505	Significativo

Resumen de Respuestas

resumen_df <- data.frame(

  Pregunta = c(

    "1. Que conclusiones debe obtener Julie?",

    "2. Como se compara con el Caso 3-1A?",

    "3. Que datos usar para el modelo?"
  ),

  Respuesta = c(

    "Las ventas muestran tendencia creciente cercana al 15 por ciento y fuerte estacionalidad anual. Los mayores valores ocurren en Enero y Febrero y los menores en Diciembre. La serie NO es estacionaria y requiere diferenciacion. El patron es regular y adecuado para pronosticos SARIMA.",

    paste0(
      "El patron de Murphy Brothers es similar al mercado minorista con correlacion r = ",
      round(cor_val,3),
      ". Ambos presentan la misma estacionalidad y los mismos meses pico y valle."
    ),

    sprintf(
      "Julie debe usar los 48 meses disponibles con frecuencia mensual y aplicar el modelo ARIMA(%d,%d,%d)(%d,%d,%d)[12]. El modelo captura tendencia y estacionalidad y debe actualizarse con nuevos datos.",
      p,d,q,P,D,Q
    )
  )
)



kable(

  resumen_df,

  col.names = c(
    "Pregunta",
    "Respuesta Sintetizada"
  ),

  align = c("l","l"),

  caption = "Resumen Ejecutivo - Inciso ii Caso 3-1B"

) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(
  1,
  bold = TRUE,
  width = "30%",
  background = "#e8f4fd"
) |>

column_spec(
  2,
  width = "70%"
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Resumen Ejecutivo - Inciso ii Caso 3-1B
Pregunta	Respuesta Sintetizada
Que conclusiones debe obtener Julie?	Las ventas muestran tendencia creciente cercana al 15 por ciento y fuerte estacionalidad anual. Los mayores valores ocurren en Enero y Febrero y los menores en Diciembre. La serie NO es estacionaria y requiere diferenciacion. El patron es regular y adecuado para pronosticos SARIMA.
Como se compara con el Caso 3-1A?	El patron de Murphy Brothers es similar al mercado minorista con correlacion r = 0.98. Ambos presentan la misma estacionalidad y los mismos meses pico y valle.
Que datos usar para el modelo?	Julie debe usar los 48 meses disponibles con frecuencia mensual y aplicar el modelo ARIMA(0,0,3)(1,1,0)[12]. El modelo captura tendencia y estacionalidad y debe actualizarse con nuevos datos.

INCISO3:Calcule los pronósticos para los periodos correspondientes con IC al 94% y utilizando la metodología de B-J.

Serie de Tiempo y Modelo

# ── Datos Murphy Brothers ────────────────────────────────────────────────────
ventas <- c(
  4906,5068,4710,4792,4638,4670,4574,4477,4571,4370,4293,3911,
  5389,5507,4727,5030,4926,4847,4814,4744,4844,4769,4483,4120,
  5270,5835,5147,5354,5290,5271,5328,5164,5372,5313,4924,4552,
  6283,6051,5298,5659,5343,5461,5568,5287,5555,5501,5201,4826
)

murphy_ts <- ts(ventas, start = c(2022, 1), frequency = 12)

# ── Modelo Box-Jenkins (búsqueda exhaustiva) ─────────────────────────────────
murphy_model <- auto.arima(
  murphy_ts,
  stepwise      = FALSE,
  approximation = FALSE,
  seasonal      = TRUE,
  ic            = "aicc"
)

# Extraer parámetros
p <- murphy_model$arma[1]; q <- murphy_model$arma[2]
P <- murphy_model$arma[3]; Q <- murphy_model$arma[4]
d <- murphy_model$arma[6]; D <- murphy_model$arma[7]
s <- frequency(murphy_ts)

cat("Modelo identificado por auto.arima():\n")

## Modelo identificado por auto.arima():

cat(sprintf("  ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n\n", p,d,q,P,D,Q,s))

##   ARIMA(0,0,3)(1,1,0)[12]

cat("Coeficientes:\n")

## Coeficientes:

print(round(coef(murphy_model), 5))

##      ma1      ma2      ma3     sar1    drift 
##  0.06992  0.00180  0.51741 -0.81121 26.83122

cat("\nVarianza residual (σ²):", round(murphy_model$sigma2, 4), "\n")

## 
## Varianza residual (σ²): 13833.38

cat("AICc:", round(murphy_model$aicc, 3), "\n")

## AICc: 469.078

Ecuación del Modelo

El modelo seleccionado por la metodología Box-Jenkins es:

\[\boxed{\text{ARIMA}(0,0,3)(1,1,0)[12]}\]

En notación de operador de retardo $B$:

\[\Phi_{1}(B^{12})\;\phi_{0}(B)\; (1-B)^{0}(1-B^{12})^{1}\;Y_t \;=\; \Theta_{0}(B^{12})\;\theta_{3}(B)\;\varepsilon_t, \quad \varepsilon_t \sim \text{RB}(0,\,1.383338\times 10^{4})\]

Metodología Box-Jenkins — Etapas

bj_df <- data.frame(
  Etapa = c("1. Identificacion","2. Estimacion","3. Diagnostico","4. Pronostico"),
  Procedimiento = c(
    "ACF/PACF + tests ADF y KPSS → serie NO estacionaria; ndiffs=1, nsdiffs=1",
    sprintf("auto.arima() → ARIMA(%d,%d,%d)(%d,%d,%d)[%d] mínimo AICc = %.2f",
            p,d,q,P,D,Q,s, murphy_model$aicc),
    "Ljung-Box sobre residuos → residuos = ruido blanco (modelo válido)",
    "forecast() con h=12 meses y level=94 → IC al 94%"
  ),
  Estado = c("✔ Completada","✔ Completada","✔ Completada","✔ En curso")
)

kable(bj_df, col.names = c("Etapa","Procedimiento","Estado"),
      align   = c("l","l","c"),
      caption = "Etapas de la Metodología Box-Jenkins — Murphy Brothers") |>
  kable_styling(bootstrap_options = c("striped","hover","bordered"),
                full_width = TRUE, font_size = 13) |>
  column_spec(1, bold = TRUE, width = "18%") |>
  column_spec(3, bold = TRUE, color = "#155724") |>
  row_spec(4, background = "#d4edda") |>
  row_spec(0, bold = TRUE, background = "#343a40", color = "white")

Etapas de la Metodología Box-Jenkins — Murphy Brothers
Etapa	Procedimiento	Estado
Identificacion	ACF/PACF + tests ADF y KPSS → serie NO estacionaria; ndiffs=1, nsdiffs=1	✔ Completada
Estimacion	auto.arima() → ARIMA(0,0,3)(1,1,0)[12] mínimo AICc = 469.08	✔ Completada
Diagnostico	Ljung-Box sobre residuos → residuos = ruido blanco (modelo válido)	✔ Completada
Pronostico	forecast() con h=12 meses y level=94 → IC al 94%	✔ En curso

Diagnóstico Final del Modelo

checkresiduals(murphy_model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,3)(1,1,0)[12] with drift
## Q* = 8.7316, df = 6, p-value = 0.1892
## 
## Model df: 4.   Total lags used: 10

lb <- Box.test(residuals(murphy_model), lag = 12, type = "Ljung-Box")
cat("Ljung-Box Test (lag = 12):\n")

## Ljung-Box Test (lag = 12):

cat("  Estadistico Q :", round(lb$statistic, 4), "\n")

##   Estadistico Q : 9.9963

cat("  p-valor       :", round(lb$p.value,   4), "\n")

##   p-valor       : 0.6163

cat("  Conclusion    :", ifelse(lb$p.value > 0.05,
    "✔ Residuos = Ruido Blanco → Modelo válido para pronosticar",
    "✘ Residuos con estructura → Revisar modelo"), "\n")

##   Conclusion    : ✔ Residuos = Ruido Blanco → Modelo válido para pronosticar

Los residuos se comportan como ruido blanco → el modelo está bien especificado y es apto para generar pronósticos confiables.

Pronósticos — 12 Periodos con IC al 94%

# h = 12 meses (enero a diciembre 2026)
# level = 94  → Intervalo de Confianza al 94%
murphy_fc <- forecast(murphy_model, h = 12, level = 94)
print(murphy_fc)

##          Point Forecast    Lo 94    Hi 94
## Jan 2026       6031.557 5810.346 6252.768
## Feb 2026       6496.620 6274.869 6718.371
## Mar 2026       5768.379 5546.627 5990.130
## Apr 2026       5994.745 5745.198 6244.292
## May 2026       5883.170 5633.623 6132.717
## Jun 2026       5890.034 5640.487 6139.581
## Jul 2026       5956.473 5706.926 6206.021
## Aug 2026       5770.385 5520.838 6019.932
## Sep 2026       5989.713 5740.165 6239.260
## Oct 2026       5931.656 5682.109 6181.204
## Nov 2026       5559.459 5309.912 5809.006
## Dec 2026       5186.892 4937.345 5436.439

Tabla de Pronósticos

meses_nombres <- c(
  "Enero","Febrero","Marzo","Abril",
  "Mayo","Junio","Julio","Agosto",
  "Septiembre","Octubre","Noviembre","Diciembre"
)

fc_df <- data.frame(

  Periodo = paste(meses_nombres, 2026),

  Pronostico = round(murphy_fc$mean, 2),

  LI_94 = round(murphy_fc$lower, 2),

  LS_94 = round(murphy_fc$upper, 2),

  Amplitud = round(
    murphy_fc$upper - murphy_fc$lower,
    2
  )
)



kable(

  fc_df,

  col.names = c(
    "Periodo",
    "Pronostico",
    "Limite Inferior 94%",
    "Limite Superior 94%",
    "Amplitud IC"
  ),

  align = c("l","r","r","r","r"),

  caption = "Pronosticos Murphy Brothers 2026 - IC 94%"

) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(
  2,
  bold = TRUE,
  color = "#155724"
) |>

column_spec(
  3,
  color = "#721c24"
) |>

column_spec(
  4,
  color = "#0c5460"
) |>

column_spec(
  5,
  color = "#856404"
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
) |>

row_spec(
  c(1,2),
  background = "#d4edda"
)

Pronosticos Murphy Brothers 2026 - IC 94%
Periodo	Pronostico	Limite Inferior 94%	Limite Superior 94%	Amplitud IC
Enero 2026	6031.56	5810.35	6252.77	442.42
Febrero 2026	6496.62	6274.87	6718.37	443.50
Marzo 2026	5768.38	5546.63	5990.13	443.50
Abril 2026	5994.74	5745.20	6244.29	499.09
Mayo 2026	5883.17	5633.62	6132.72	499.09
Junio 2026	5890.03	5640.49	6139.58	499.09
Julio 2026	5956.47	5706.93	6206.02	499.09
Agosto 2026	5770.39	5520.84	6019.93	499.09
Septiembre 2026	5989.71	5740.17	6239.26	499.09
Octubre 2026	5931.66	5682.11	6181.20	499.09
Noviembre 2026	5559.46	5309.91	5809.01	499.09
Diciembre 2026	5186.89	4937.35	5436.44	499.09

cat("── Resumen de pronósticos 2026 ──────────────────────\n")

## ── Resumen de pronósticos 2026 ──────────────────────

cat("Mes de mayor venta pronosticada :",
    meses_nombres[which.max(murphy_fc$mean)], "→",
    comma(round(max(murphy_fc$mean))), "unidades\n")

## Mes de mayor venta pronosticada : Febrero → 6,497 unidades

cat("Mes de menor venta pronosticada :",
    meses_nombres[which.min(murphy_fc$mean)], "→",
    comma(round(min(murphy_fc$mean))), "unidades\n")

## Mes de menor venta pronosticada : Diciembre → 5,187 unidades

cat("Total anual pronosticado 2026   :",
    comma(round(sum(murphy_fc$mean))), "unidades\n")

## Total anual pronosticado 2026   : 70,459 unidades

Gráfico de Pronósticos

Gráfico Principal con IC 94%

autoplot(murphy_fc) +
  autolayer(murphy_ts, series = "Ventas históricas", color = "#2C7BB6") +
  scale_color_manual(values = c("Ventas históricas" = "#2C7BB6"),
                     name = NULL) +
  scale_y_continuous(labels = comma) +
  scale_x_continuous(breaks = 2022:2026) +
  labs(
    title    = "Murphy Brothers — Pronósticos 2026 con IC al 94%",
    subtitle = sprintf("Modelo Box-Jenkins: ARIMA(%d,%d,%d)(%d,%d,%d)[%d]",
                       p,d,q,P,D,Q,s),
    x        = "Tiempo",
    y        = "Ventas (unidades)",
    caption  = "Banda sombreada = Intervalo de Confianza al 94%"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(color = "#555555"),
    plot.caption  = element_text(color = "#777777", size = 10),
    legend.position = "bottom"
  )

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Zoom: Solo Pronósticos 2026

# Construir data frame de pronósticos para ggplot personalizado
fechas_fc <- seq(as.Date("2026-01-01"), by = "month", length.out = 12)

df_fc <- data.frame(
  Mes  = factor(meses_nombres, levels = meses_nombres),
  Yhat = as.numeric(murphy_fc$mean),
  LI   = as.numeric(murphy_fc$lower),
  LS   = as.numeric(murphy_fc$upper)
)

ggplot(df_fc, aes(x = Mes, y = Yhat, group = 1)) +
  geom_ribbon(aes(ymin = LI, ymax = LS), fill = "#fdae61", alpha = 0.5) +
  geom_line(color = "#D7191C", linewidth = 1.1) +
  geom_point(color = "#D7191C", size = 3.5) +
  geom_text(aes(label = comma(round(Yhat))),
            vjust = -1.0, size = 3.4, fontface = "bold", color = "#D7191C") +
  geom_text(aes(y = LI, label = comma(round(LI))),
            vjust = 1.4, size = 2.8, color = "#721c24", alpha = 0.8) +
  geom_text(aes(y = LS, label = comma(round(LS))),
            vjust = -0.5, size = 2.8, color = "#0c5460", alpha = 0.8) +
  scale_y_continuous(labels = comma,
                     limits = c(min(df_fc$LI)*0.97, max(df_fc$LS)*1.05)) +
  labs(
    title    = "Pronosticos Mensuales — Murphy Brothers 2026",
    subtitle = "IC al 94% | Banda naranja = Intervalo de Confianza",
    x        = "Mes",
    y        = "Ventas pronosticadas (unidades)",
    caption  = "Rojo = Pronostico puntual | Naranja claro = LI | Azul oscuro = LS"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(color = "#555555"),
    axis.text.x   = element_text(angle = 30, hjust = 1)
  )

Interpretación del Intervalo de Confianza al 94%

\[\hat{Y}_{t+h} \pm z_{0.03}\;\hat{\sigma}_h\]

donde $z_{0.03} = 1.8808$ (percentil 97° de la normal estándar, puesto que el IC es bilateral al 94%) y $\hat{\sigma}_h$ es el error estándar del pronóstico al horizonte $h$.

z_94 <- qnorm(0.97)

cat(
  "z critico para IC al 94% (bilateral):",
  round(z_94, 4),
  "\n"
)

## z critico para IC al 94% (bilateral): 1.8808

cat(
  "P(-1.8808 < Z < 1.8808) = 0.94\n\n"
)

## P(-1.8808 < Z < 1.8808) = 0.94

# Error estandar del pronostico
se_fc <- (
  as.numeric(murphy_fc$upper) -
  as.numeric(murphy_fc$mean)
) / z_94



se_df <- data.frame(

  Mes = meses_nombres,

  Pronostico = round(
    as.numeric(murphy_fc$mean),
    2
  ),

  Error_Std = round(
    se_fc,
    2
  ),

  LI_94 = round(
    as.numeric(murphy_fc$lower),
    2
  ),

  LS_94 = round(
    as.numeric(murphy_fc$upper),
    2
  )
)



kable(

  se_df,

  col.names = c(
    "Mes",
    "Pronostico",
    "Error Std",
    "LI 94%",
    "LS 94%"
  ),

  align = c("l","r","r","r","r"),

  caption = "Calculo del IC 94% - Metodo Box Jenkins"

) |>

kable_styling(
  bootstrap_options = c("striped","hover","bordered"),
  full_width = TRUE,
  font_size = 13
) |>

column_spec(
  2,
  bold = TRUE,
  color = "#155724"
) |>

column_spec(
  3,
  color = "#6c757d"
) |>

row_spec(
  0,
  bold = TRUE,
  background = "#343a40",
  color = "white"
)

Calculo del IC 94% - Metodo Box Jenkins
Mes	Pronostico	Error Std	LI 94%	LS 94%
Enero	6031.56	117.62	5810.35	6252.77
Febrero	6496.62	117.90	6274.87	6718.37
Marzo	5768.38	117.90	5546.63	5990.13
Abril	5994.74	132.68	5745.20	6244.29
Mayo	5883.17	132.68	5633.62	6132.72
Junio	5890.03	132.68	5640.49	6139.58
Julio	5956.47	132.68	5706.93	6206.02
Agosto	5770.39	132.68	5520.84	6019.93
Septiembre	5989.71	132.68	5740.17	6239.26
Octubre	5931.66	132.68	5682.11	6181.20
Noviembre	5559.46	132.68	5309.91	5809.01
Diciembre	5186.89	132.68	4937.35	5436.44

Nota: El error estándar del pronóstico $\hat{\sigma}_h$ crece con el horizonte $h$ porque la incertidumbre se acumula; por eso la amplitud del IC aumenta mes a mes.

Conclusión

conc_df <- data.frame(
  Item = c("Modelo B-J","IC utilizado","Horizonte","Pico pronosticado",
           "Valle pronosticado","Total anual 2026"),
  Resultado = c(
    sprintf("ARIMA(%d,%d,%d)(%d,%d,%d)[12]", p,d,q,P,D,Q),
    "94%  (z = 1.8808)",
    "12 meses (enero – diciembre 2026)",
    paste0(meses_nombres[which.max(murphy_fc$mean)], " 2026 → ",
           comma(round(max(murphy_fc$mean))), " unidades"),
    paste0(meses_nombres[which.min(murphy_fc$mean)], " 2026 → ",
           comma(round(min(murphy_fc$mean))), " unidades"),
    paste0(comma(round(sum(murphy_fc$mean))), " unidades")
  )
)

kable(conc_df, col.names = c("Elemento","Resultado"),
      align   = c("l","l"),
      caption = "Resumen del Pronóstico — Inciso iii") |>
  kable_styling(bootstrap_options = c("striped","hover","bordered"),
                full_width = FALSE, font_size = 13) |>
  column_spec(1, bold = TRUE, width = "35%") |>
  column_spec(2, width = "65%", color = "#155724") |>
  row_spec(0, bold = TRUE, background = "#343a40", color = "white")

Resumen del Pronóstico — Inciso iii
Elemento	Resultado
Modelo B-J	ARIMA(0,0,3)(1,1,0)[12]
IC utilizado	94% (z = 1.8808)
Horizonte	12 meses (enero – diciembre 2026)
Pico pronosticado	Febrero 2026 → 6,497 unidades
Valle pronosticado	Diciembre 2026 → 5,187 unidades
Total anual 2026	70,459 unidades

Resuelva los siguientes ejercicios en R, realizando lo mismo que los incisos i y iii: (actualizando y configurando,…) y tomando en cuenta que las ST se encuentran en la carpeta: CASOS PRÁCTICOS en teans de ST

En la tabla siguiente, se presenta el número de matrimonios en Estados Unidos. Calcule las primeras diferencias para tales datos. Grafique los datos originales y los datos de diferencia como una serie de tiempo. ¿Hay tendencia en cualquiera de estas series? Explique su respuesta

Datos y Primeras Diferencias

Carga de datos

# ── Datos originales ──────────────────────────────────────────────────────────
anio <- 1985:2004

matrimonios <- c(
  2413, 2407, 2403, 2396, 2403, 2443, 2371, 2362, 2334, 2362,
  2336, 2344, 2384, 2244, 2358, 2329, 2345, 2254, 2245, 2279
)

df <- data.frame(Anio = anio, Matrimonios = matrimonios)

knitr::kable(
  df,
  caption = "Matrimonios en Estados Unidos 1985–2004 (miles)",
  col.names = c("Anio", "Matrimonios (miles)"),
  align = "cc"
)

Matrimonios en Estados Unidos 1985–2004 (miles)
Anio	Matrimonios (miles)
1985	2413
1986	2407
1987	2403
1988	2396
1989	2403
1990	2443
1991	2371
1992	2362
1993	2334
1994	2362
1995	2336
1996	2344
1997	2384
1998	2244
1999	2358
2000	2329
2001	2345
2002	2254
2003	2245
2004	2279

Primeras diferencias

\[\Delta Y_t = Y_t - Y_{t-1}\]

primera_diff <- diff(matrimonios)

df_diff <- data.frame(
  Anio       = anio[-1],
  Matrimonios = matrimonios[-1],
  Primera_Diff = primera_diff
)

knitr::kable(
  df_diff,
  caption  = "Serie original y primeras diferencias",
  col.names = c("Anio", "Matrimonios (miles)", "Primera Diferencia"),
  align    = "ccc"
)

Serie original y primeras diferencias
Anio	Matrimonios (miles)	Primera Diferencia
1986	2407	-6
1987	2403	-4
1988	2396	-7
1989	2403	7
1990	2443	40
1991	2371	-72
1992	2362	-9
1993	2334	-28
1994	2362	28
1995	2336	-26
1996	2344	8
1997	2384	40
1998	2244	-140
1999	2358	114
2000	2329	-29
2001	2345	16
2002	2254	-91
2003	2245	-9
2004	2279	34

Gráfica de la serie original y diferenciada

par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))

# Serie original
plot(anio, matrimonios,
     type = "b", col = "steelblue", pch = 16, lwd = 2,
     xlab = "Anio", ylab = "Matrimonios (miles)",
     main = "Serie Original: Matrimonios en EE.UU. (1985–2004)")
grid(col = "gray85")

# Primeras diferencias
plot(anio[-1], primera_diff,
     type = "b", col = "firebrick", pch = 16, lwd = 2,
     xlab = "Anio", ylab = "Δ Matrimonios (miles)",
     main = "Primeras Diferencias de Matrimonios en EE.UU.")
abline(h = 0, lty = 2, col = "gray50")
grid(col = "gray85")

par(mfrow = c(1, 1))

Interpretación:

La serie original muestra una leve tendencia decreciente con algunas fluctuaciones irregulares a lo largo del período.
Las primeras diferencias oscilan alrededor de cero sin patrón sistemático claro, lo que sugiere que la diferenciación de primer orden puede ser suficiente para inducir estacionariedad.

Análisis de la Serie de Tiempo (AST)

Configuración con `ts()` — año final 2025

library(forecast)
library(tseries)
library(ggplot2)
library(gridExtra)

# Serie de tiempo 1985–2004
matrimonios_ts <- ts(matrimonios, start = 1985, frequency = 1)

# Extensión de la serie hasta 2025 con NA para pronóstico posterior
# (la serie observada va de 1985 a 2004; los años 2005-2025 serán pronósticos)
cat("Serie de tiempo configurada:\n")

## Serie de tiempo configurada:

print(matrimonios_ts)

## Time Series:
## Start = 1985 
## End = 2004 
## Frequency = 1 
##  [1] 2413 2407 2403 2396 2403 2443 2371 2362 2334 2362 2336 2344 2384 2244 2358
## [16] 2329 2345 2254 2245 2279

cat(sprintf("\nInicio: %d | Fin: %d | Frecuencia: %d\n",
            start(matrimonios_ts)[1],
            end(matrimonios_ts)[1],
            frequency(matrimonios_ts)))

## 
## Inicio: 1985 | Fin: 2004 | Frecuencia: 1

Gráfica de la serie con `autoplot`

autoplot(matrimonios_ts) +
  labs(
    title  = "Matrimonios en EE.UU. — Serie de Tiempo (1985–2004)",
    x      = "Anio",
    y      = "Matrimonios (miles)"
  ) +
  theme_minimal(base_size = 13) +
  geom_point(color = "steelblue", size = 2)

Descomposición y patrones (AST)

Como la serie es anual (frecuencia = 1), no existe componente estacional en sentido estricto. Se analiza tendencia y componente irregular.

# Media móvil centrada (suavizado) para identificar tendencia
trend_ma <- stats::filter(matrimonios_ts, filter = rep(1/5, 5), sides = 2)

plot(matrimonios_ts,
     col = "steelblue", lwd = 2,
     main = "Serie Original y Tendencia (MA-5)",
     ylab = "Matrimonios (miles)", xlab = "Anio")
lines(trend_ma, col = "firebrick", lwd = 2, lty = 2)
legend("topright",
       legend = c("Serie original", "Tendencia MA(5)"),
       col    = c("steelblue", "firebrick"),
       lty    = c(1, 2), lwd = 2, bty = "n")
grid(col = "gray85")

Patrones identificados:

Componente	Descripción
Tendencia	Leve tendencia decreciente a lo largo del período 1985–2004
Estacionalidad	No aplica (datos anuales, frecuencia = 1)
Irregularidad	Fluctuaciones irregulares presentes, especialmente en 1990, 1998 y 2002
Tipo de serie	Aparentemente no estacionaria en media (tendencia decreciente)

Prueba de Estacionariedad

Prueba Aumentada de Dickey-Fuller (ADF)

\[H_0: \text{La serie tiene raíz unitaria (no estacionaria)}\] \[H_1: \text{La serie es estacionaria}\]

adf_orig <- adf.test(matrimonios_ts, alternative = "stationary")
cat("=== ADF — Serie Original ===\n")

## === ADF — Serie Original ===

print(adf_orig)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  matrimonios_ts
## Dickey-Fuller = -3.5916, Lag order = 2, p-value = 0.05116
## alternative hypothesis: stationary

if (adf_orig$p.value < 0.05) {
  cat("\n→ p-valor < 0.05: Se RECHAZA H₀ → La serie ES estacionaria.\n")
} else {
  cat("\n→ p-valor ≥ 0.05: NO se rechaza H₀ → La serie NO es estacionaria.\n")
}

## 
## → p-valor ≥ 0.05: NO se rechaza H₀ → La serie NO es estacionaria.

Prueba KPSS

\[H_0: \text{La serie es estacionaria}\] \[H_1: \text{La serie no es estacionaria (tiene raíz unitaria)}\]

kpss_orig <- kpss.test(matrimonios_ts, null = "Level")
cat("=== KPSS — Serie Original ===\n")

## === KPSS — Serie Original ===

print(kpss_orig)

## 
##  KPSS Test for Level Stationarity
## 
## data:  matrimonios_ts
## KPSS Level = 0.68692, Truncation lag parameter = 2, p-value = 0.01473

if (kpss_orig$p.value < 0.05) {
  cat("\n→ p-valor < 0.05: Se RECHAZA H₀ → La serie NO es estacionaria.\n")
} else {
  cat("\n→ p-valor ≥ 0.05: NO se rechaza H₀ → La serie ES estacionaria.\n")
}

## 
## → p-valor < 0.05: Se RECHAZA H₀ → La serie NO es estacionaria.

ADF sobre primeras diferencias

diff_ts <- diff(matrimonios_ts)
adf_diff <- adf.test(diff_ts, alternative = "stationary")
cat("=== ADF — Primera Diferencia ===\n")

## === ADF — Primera Diferencia ===

print(adf_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff_ts
## Dickey-Fuller = -3.9089, Lag order = 2, p-value = 0.02793
## alternative hypothesis: stationary

if (adf_diff$p.value < 0.05) {
  cat("\n→ p-valor < 0.05: Se RECHAZA H₀ → La primera diferencia ES estacionaria.\n")
} else {
  cat("\n→ p-valor ≥ 0.05: NO se rechaza H₀ → La primera diferencia NO es estacionaria.\n")
}

## 
## → p-valor < 0.05: Se RECHAZA H₀ → La primera diferencia ES estacionaria.

Función de Autocorrelación (ACF) y Autocorrelación Parcial (PACF)

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))

acf(matrimonios_ts,  main = "ACF — Serie Original",         lag.max = 15)
pacf(matrimonios_ts, main = "PACF — Serie Original",        lag.max = 15)
acf(diff_ts,         main = "ACF — Primera Diferencia",     lag.max = 15)
pacf(diff_ts,        main = "PACF — Primera Diferencia",    lag.max = 15)

par(mfrow = c(1, 1))

Identificación automática del modelo con `auto.arima()`

modelo_arima <- auto.arima(
  matrimonios_ts,
  seasonal     = FALSE,   # serie anual, sin estacionalidad
  stepwise     = FALSE,
  approximation = FALSE,
  ic           = "aicc",
  trace        = TRUE
)

## 
##  ARIMA(0,1,0)                    : 207.2958
##  ARIMA(0,1,0) with drift         : 209.4758
##  ARIMA(0,1,1)                    : 202.7636
##  ARIMA(0,1,1) with drift         : Inf
##  ARIMA(0,1,2)                    : 205.0621
##  ARIMA(0,1,2) with drift         : Inf
##  ARIMA(0,1,3)                    : 208.0176
##  ARIMA(0,1,3) with drift         : Inf
##  ARIMA(0,1,4)                    : Inf
##  ARIMA(0,1,4) with drift         : Inf
##  ARIMA(0,1,5)                    : Inf
##  ARIMA(0,1,5) with drift         : Inf
##  ARIMA(1,1,0)                    : 203.3528
##  ARIMA(1,1,0) with drift         : 204.7933
##  ARIMA(1,1,1)                    : 205.1195
##  ARIMA(1,1,1) with drift         : Inf
##  ARIMA(1,1,2)                    : 208.3346
##  ARIMA(1,1,2) with drift         : Inf
##  ARIMA(1,1,3)                    : Inf
##  ARIMA(1,1,3) with drift         : Inf
##  ARIMA(1,1,4)                    : Inf
##  ARIMA(1,1,4) with drift         : Inf
##  ARIMA(2,1,0)                    : 205.6182
##  ARIMA(2,1,0) with drift         : 206.5008
##  ARIMA(2,1,1)                    : 208.3757
##  ARIMA(2,1,1) with drift         : Inf
##  ARIMA(2,1,2)                    : Inf
##  ARIMA(2,1,2) with drift         : Inf
##  ARIMA(2,1,3)                    : Inf
##  ARIMA(2,1,3) with drift         : Inf
##  ARIMA(3,1,0)                    : 207.9021
##  ARIMA(3,1,0) with drift         : 207.6583
##  ARIMA(3,1,1)                    : 211.5416
##  ARIMA(3,1,1) with drift         : Inf
##  ARIMA(3,1,2)                    : Inf
##  ARIMA(3,1,2) with drift         : Inf
##  ARIMA(4,1,0)                    : 211.2476
##  ARIMA(4,1,0) with drift         : 208.901
##  ARIMA(4,1,1)                    : 213.9713
##  ARIMA(4,1,1) with drift         : Inf
##  ARIMA(5,1,0)                    : 214.2257
##  ARIMA(5,1,0) with drift         : 213.7678
## 
## 
## 
##  Best model: ARIMA(0,1,1)

cat("\n=== Modelo seleccionado por auto.arima() ===\n")

## 
## === Modelo seleccionado por auto.arima() ===

summary(modelo_arima)

## Series: matrimonios_ts 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -0.5708
## s.e.   0.1530
## 
## sigma^2 = 2033:  log likelihood = -99.01
## AIC=202.01   AICc=202.76   BIC=203.9
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -15.41511 42.77714 29.65635 -0.6846123 1.281196 0.7958626
##                    ACF1
## Training set -0.3132728

Ecuación del modelo deducido

p <- modelo_arima$arma[1]   # AR
d <- modelo_arima$arma[6]   # diferenciación
q <- modelo_arima$arma[2]   # MA

cat(sprintf("Modelo identificado: ARIMA(%d, %d, %d)\n\n", p, d, q))

## Modelo identificado: ARIMA(0, 1, 1)

# Coeficientes
coef_modelo <- coef(modelo_arima)
print(coef_modelo)

##       ma1 
## -0.570791

La ecuación general del modelo $\text{ARIMA}(p, d, q)$ se expresa como:

\[\phi(B)(1-B)^d Y_t = \theta(B)\varepsilon_t\]

Donde:

$B$ es el operador de retardo ($BY_t = Y_{t-1}$)
$\phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p$ es el polinomio autorregresivo
$\theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q$ es el polinomio de medias móviles
$(1-B)^d$ es el operador de diferenciación de orden $d$
$\varepsilon_t \sim \text{WN}(0, \sigma^2)$ es el ruido blanco

# Presentar la ecuación específica con los coeficientes estimados
cat("╔══════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════╗

cat(sprintf("  Modelo: ARIMA(%d, %d, %d)\n", p, d, q))

##   Modelo: ARIMA(0, 1, 1)

cat("──────────────────────────────────────────────────────\n")

## ──────────────────────────────────────────────────────

if (length(coef_modelo) > 0) {
  for (nm in names(coef_modelo)) {
    cat(sprintf("  %s = %.4f\n", nm, coef_modelo[nm]))
  }
}

##   ma1 = -0.5708

cat("╚══════════════════════════════════════════════════════╝\n")

## ╚══════════════════════════════════════════════════════╝

Diagnóstico de residuos

checkresiduals(modelo_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)
## Q* = 3.6231, df = 3, p-value = 0.3051
## 
## Model df: 1.   Total lags used: 4

lb_test <- Box.test(residuals(modelo_arima), lag = 10, type = "Ljung-Box")
cat("=== Prueba Ljung-Box sobre residuos ===\n")

## === Prueba Ljung-Box sobre residuos ===

print(lb_test)

## 
##  Box-Ljung test
## 
## data:  residuals(modelo_arima)
## X-squared = 10.881, df = 10, p-value = 0.3668

if (lb_test$p.value > 0.05) {
  cat("\n→ p-valor > 0.05: Los residuos son RUIDO BLANCO (modelo adecuado).\n")
} else {
  cat("\n→ p-valor ≤ 0.05: Los residuos NO son ruido blanco (posible mejora del modelo).\n")
}

## 
## → p-valor > 0.05: Los residuos son RUIDO BLANCO (modelo adecuado).

Pronósticos con Metodología Box-Jenkins — IC al 94%

Pronósticos para el período 2005–2025

Usando la metodología Box-Jenkins (B-J), se generan pronósticos para los 21 períodos restantes (2005–2025) con un intervalo de confianza del 94%.

# Número de períodos a pronosticar: 2005 a 2025 = 21 años
h_pronostico <- 2025 - 2004

pronosticos <- forecast(modelo_arima,
                        h     = h_pronostico,
                        level = 94)   # IC al 94%

cat("=== Pronósticos ARIMA — 2005 a 2025 (IC 94%) ===\n")

## === Pronósticos ARIMA — 2005 a 2025 (IC 94%) ===

print(pronosticos)

##      Point Forecast    Lo 94    Hi 94
## 2005       2277.921 2193.114 2362.728
## 2006       2277.921 2185.633 2370.210
## 2007       2277.921 2178.714 2377.129
## 2008       2277.921 2172.247 2383.596
## 2009       2277.921 2166.153 2389.689
## 2010       2277.921 2160.375 2395.467
## 2011       2277.921 2154.869 2400.974
## 2012       2277.921 2149.598 2406.245
## 2013       2277.921 2144.535 2411.307
## 2014       2277.921 2139.658 2416.185
## 2015       2277.921 2134.946 2420.896
## 2016       2277.921 2130.386 2425.457
## 2017       2277.921 2125.962 2429.881
## 2018       2277.921 2121.663 2434.179
## 2019       2277.921 2117.479 2438.363
## 2020       2277.921 2113.402 2442.440
## 2021       2277.921 2109.423 2446.419
## 2022       2277.921 2105.537 2450.306
## 2023       2277.921 2101.735 2454.107
## 2024       2277.921 2098.015 2457.828
## 2025       2277.921 2094.369 2461.473

Tabla de pronósticos

anios_pron <- 2005:2025


tabla_pron <- data.frame(

  Anio = anios_pron,

  Pronostico = round(
    as.numeric(pronosticos$mean),
    2
  ),

  LI_94 = round(
    as.numeric(pronosticos$lower),
    2
  ),

  LS_94 = round(
    as.numeric(pronosticos$upper),
    2
  )
)



knitr::kable(

  tabla_pron,

  caption = "Pronosticos de Matrimonios en EEUU 2005-2025 con IC 94%",

  col.names = c(
    "Anio",
    "Pronostico miles",
    "Limite Inferior 94%",
    "Limite Superior 94%"
  ),

  align = "cccc"
)

Pronosticos de Matrimonios en EEUU 2005-2025 con IC 94%
Anio	Pronostico miles	Limite Inferior 94%	Limite Superior 94%
2005	2277.92	2193.11	2362.73
2006	2277.92	2185.63	2370.21
2007	2277.92	2178.71	2377.13
2008	2277.92	2172.25	2383.60
2009	2277.92	2166.15	2389.69
2010	2277.92	2160.38	2395.47
2011	2277.92	2154.87	2400.97
2012	2277.92	2149.60	2406.24
2013	2277.92	2144.54	2411.31
2014	2277.92	2139.66	2416.18
2015	2277.92	2134.95	2420.90
2016	2277.92	2130.39	2425.46
2017	2277.92	2125.96	2429.88
2018	2277.92	2121.66	2434.18
2019	2277.92	2117.48	2438.36
2020	2277.92	2113.40	2442.44
2021	2277.92	2109.42	2446.42
2022	2277.92	2105.54	2450.31
2023	2277.92	2101.74	2454.11
2024	2277.92	2098.01	2457.83
2025	2277.92	2094.37	2461.47

Gráfica de pronósticos

autoplot(pronosticos) +
  autolayer(matrimonios_ts, series = "Observado", color = "steelblue", size = 1) +
  labs(
    title    = "Pronósticos de Matrimonios en EE.UU. (2005–2025)",
    subtitle = "Metodología Box-Jenkins — ARIMA | IC al 94%",
    x        = "Anio",
    y        = "Matrimonios (miles)",
    color    = "Serie"
  ) +
  theme_minimal(base_size = 13) +
  scale_color_manual(values = c("Observado" = "steelblue")) +
  theme(legend.position = "bottom")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the forecast package.
##   Please report the issue at <https://github.com/robjhyndman/forecast/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Ignoring unknown labels:
## • colour : "Serie"

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Gráfica completa: observado + pronóstico

# Serie completa con extensión a 2025
anios_todos  <- c(anio, anios_pron)
valores_obs  <- c(matrimonios, rep(NA, h_pronostico))
valores_pron <- c(rep(NA, length(matrimonios)), as.numeric(pronosticos$mean))
li_94        <- c(rep(NA, length(matrimonios)), as.numeric(pronosticos$lower))
ls_94        <- c(rep(NA, length(matrimonios)), as.numeric(pronosticos$upper))

plot(anios_todos, valores_obs,
     type = "b", col = "steelblue", lwd = 2, pch = 16,
     ylim = range(c(li_94, ls_94, matrimonios), na.rm = TRUE),
     xlab = "Anio", ylab = "Matrimonios (miles)",
     main = "Serie Observada y Pronosticos de Matrimonios en EE.UU.")

# Banda de confianza
polygon(c(anios_pron, rev(anios_pron)),
        c(ls_94[!is.na(ls_94)], rev(li_94[!is.na(li_94)])),
        col = adjustcolor("salmon", alpha.f = 0.25), border = NA)

# Pronósticos puntuales
lines(anios_pron, valores_pron[!is.na(valores_pron)],
      col = "firebrick", lwd = 2, lty = 2)
points(anios_pron, valores_pron[!is.na(valores_pron)],
       col = "firebrick", pch = 17, cex = 0.8)

# Límites IC
lines(anios_pron, ls_94[!is.na(ls_94)], col = "salmon", lty = 3, lwd = 1.5)
lines(anios_pron, li_94[!is.na(li_94)], col = "salmon", lty = 3, lwd = 1.5)

legend("topright",
       legend = c("Observado", "Pronostico", "IC 94%"),
       col    = c("steelblue", "firebrick", "salmon"),
       lty    = c(1, 2, 1), lwd = c(2, 2, 8), pch = c(16, 17, NA),
       bty    = "n", pt.cex = 0.9)

grid(col = "gray88")

Resumen y Conclusiones

cat("╔════════════════════════════════════════════════════════════════╗\n")

## ╔════════════════════════════════════════════════════════════════╗

cat("             RESUMEN DEL ANALISIS DE SERIES DE TIEMPO             \n")

##              RESUMEN DEL ANALISIS DE SERIES DE TIEMPO

cat("╚════════════════════════════════════════════════════════════════╝\n\n")

## ╚════════════════════════════════════════════════════════════════╝

cat("► INCISO 1 — Primeras Diferencias\n")

## ► INCISO 1 — Primeras Diferencias

cat("  • Se calcularon las primeras diferencias ΔYt = Yt - Yt-1.\n")

##   • Se calcularon las primeras diferencias ΔYt = Yt - Yt-1.

cat("  • La serie original muestra tendencia decreciente leve.\n")

##   • La serie original muestra tendencia decreciente leve.

cat("  • La serie diferenciada oscila alrededor de 0, sin tendencia.\n\n")

##   • La serie diferenciada oscila alrededor de 0, sin tendencia.

cat("► INCISO 2 — Análisis AST\n")

## ► INCISO 2 — Análisis AST

cat(sprintf("  • Modelo identificado: ARIMA(%d, %d, %d)\n", p, d, q))

##   • Modelo identificado: ARIMA(0, 1, 1)

cat("  • Serie: NO estacionaria en la original (ADF/KPSS).\n")

##   • Serie: NO estacionaria en la original (ADF/KPSS).

cat("  • Con una diferenciación (d=1) se logra estacionariedad.\n")

##   • Con una diferenciación (d=1) se logra estacionariedad.

cat("  • Sin componente estacional (datos anuales).\n")

##   • Sin componente estacional (datos anuales).

cat("  • Residuos tipo ruido blanco (Ljung-Box).\n\n")

##   • Residuos tipo ruido blanco (Ljung-Box).

cat("► INCISO 3 — Pronósticos B-J (2005–2025, IC 94%)\n")

## ► INCISO 3 — Pronósticos B-J (2005–2025, IC 94%)

cat("  • Se proyectaron 21 períodos con intervalos al 94%.\n")

##   • Se proyectaron 21 períodos con intervalos al 94%.

cat(sprintf("  • Pronóstico 2025: %.2f miles de matrimonios.\n",
            as.numeric(pronosticos$mean)[h_pronostico]))

##   • Pronóstico 2025: 2277.92 miles de matrimonios.

cat(sprintf("  • IC 94%% para 2025: [%.2f, %.2f] miles.\n",
            as.numeric(pronosticos$lower)[h_pronostico],
            as.numeric(pronosticos$upper)[h_pronostico]))

##   • IC 94% para 2025: [2094.37, 2461.47] miles.

Tendencia en la serie original

La serie de matrimonios en EE.UU. (1985–2004) presenta una tendencia decreciente moderada: partiendo de ~2,413 miles en 1985 y llegando a ~2,279 miles en 2004. Esta tendencia no es perfectamente lineal; se observan repuntes en 1990 (2,443) y 1999 (2,358). En términos generales, sí existe tendencia en la serie original, lo cual la hace no estacionaria en media. La serie diferenciada, en cambio, no presenta tendencia evidente.

A Allie White, la jefa de la oficina de préstamos del Dominion Bank, le gustaría analizar el portafolios de préstamos del banco para los años 2001 a 2006.Los datos se presentan en la tabla siguiente.

Calcule las autocorrelaciones para los retrasos de tiempo.
Use R para graficar los datos y calcule las autocorrelaciones para los primeros seis retrasos de tiempo. ¿Esta serie de tiempo es estacionaria?
Calcule las primeras diferencias de los datos de préstamos trimestrales del Dominion Bank.
Calcule los coeficientes de autocorrelación para los retrasos de tiempo 1-6, usando los datos diferenciados.
Graficar los datos diferenciados y para calcular sus autocorrelaciones de los primeros seis retrasos de tiempo. ¿Esta serie es estacionaria? Luego aplique el test correspondiente para confirmar.

Datos del Dominion Bank

Carga y organización de datos

# ── Datos trimestrales de préstamos (millones de dólares) ──────────────────
trimestres <- c(
  "2001-Q1","2001-Q2","2001-Q3","2001-Q4",
  "2002-Q1","2002-Q2","2002-Q3","2002-Q4",
  "2003-Q1","2003-Q2","2003-Q3","2003-Q4",
  "2004-Q1","2004-Q2","2004-Q3","2004-Q4",
  "2005-Q1","2005-Q2","2005-Q3","2005-Q4",
  "2006-Q1","2006-Q2","2006-Q3","2006-Q4"
)

prestamos <- c(
  2313, 2495, 2609, 2792,
  2860, 3099, 3202, 3161,
  3399, 3471, 3545, 3851,
  4458, 4850, 5093, 5318,
  5756, 6013, 6158, 6289,
  6369, 6568, 6646, 6861
)

df <- data.frame(
  Trimestre  = trimestres,
  Prestamos  = prestamos
)

knitr::kable(
  matrix(prestamos, nrow = 6, byrow = TRUE,
         dimnames = list(2001:2006, c("Mar.31","Jun.30","Sep.30","Dic.31"))),
  caption = "Portafolio de Préstamos del Dominion Bank (millones USD)",
  align   = "cccc"
)

Portafolio de Préstamos del Dominion Bank (millones USD)
	Mar.31	Jun.30	Sep.30	Dic.31
2001	2313	2495	2609	2792
2002	2860	3099	3202	3161
2003	3399	3471	3545	3851
2004	4458	4850	5093	5318
2005	5756	6013	6158	6289
2006	6369	6568	6646	6861

Configuración como serie de tiempo trimestral

library(forecast)
library(tseries)
library(ggplot2)

# Serie trimestral: inicio 2001, trimestre 1, frecuencia 4
prestamos_ts <- ts(prestamos, start = c(2001, 1), frequency = 4)

cat("Serie de tiempo configurada:\n")

## Serie de tiempo configurada:

print(prestamos_ts)

##      Qtr1 Qtr2 Qtr3 Qtr4
## 2001 2313 2495 2609 2792
## 2002 2860 3099 3202 3161
## 2003 3399 3471 3545 3851
## 2004 4458 4850 5093 5318
## 2005 5756 6013 6158 6289
## 2006 6369 6568 6646 6861

cat(sprintf("\nInicio: %d-Q%d | Fin: %d-Q%d | Frecuencia: %d (trimestral)\n",
            start(prestamos_ts)[1], start(prestamos_ts)[2],
            end(prestamos_ts)[1],   end(prestamos_ts)[2],
            frequency(prestamos_ts)))

## 
## Inicio: 2001-Q1 | Fin: 2006-Q4 | Frecuencia: 4 (trimestral)

Inciso a) — Autocorrelaciones para los primeros 6 retrasos

Cálculo manual de las autocorrelaciones

La autocorrelación de orden $k$ se define como:

\[r_k = \frac{\sum_{t=k+1}^{n}(Y_t - \bar{Y})(Y_{t-k} - \bar{Y})}{\sum_{t=1}^{n}(Y_t - \bar{Y})^2}\]

n <- length(prestamos)

Ybar <- mean(prestamos)



# Funcion para calcular rk manualmente
rk_manual <- function(serie, k) {

  n <- length(serie)

  Ybar <- mean(serie)

  num <- sum(
    (serie[(k+1):n] - Ybar) *
    (serie[1:(n-k)] - Ybar)
  )

  den <- sum(
    (serie - Ybar)^2
  )

  return(num / den)
}



lags <- 1:6


autocorr_vals <- sapply(
  lags,
  function(k) rk_manual(prestamos, k)
)



tabla_autocorr <- data.frame(

  Retraso = lags,

  rk = round(
    autocorr_vals,
    4
  )
)



knitr::kable(

  tabla_autocorr,

  caption = "Autocorrelaciones para los primeros 6 retrasos - Serie Original",

  col.names = c(
    "Retraso k",
    "rk"
  ),

  align = "cc"
)

Autocorrelaciones para los primeros 6 retrasos - Serie Original
Retraso k	rk
1	0.8953
2	0.7884
3	0.6733
4	0.5582
5	0.4331
6	0.3094

Interpretación del Inciso a):

Los coeficientes de autocorrelación $r_1$ a $r_6$ son todos muy cercanos a 1 y decaen lentamente, lo cual es un indicador clásico de que la serie tiene tendencia fuerte y no es estacionaria. En una serie estacionaria, las autocorrelaciones decaerían rápidamente a cero para retrasos moderados.

Inciso b) — Gráfica de la serie y diagnóstico de estacionariedad

Gráfica de la serie original

autoplot(prestamos_ts) +
  geom_point(color = "steelblue", size = 2.5) +
  labs(
    title    = "Portafolio de Préstamos — Dominion Bank (2001–2006)",
    subtitle = "Datos trimestrales en millones de dolares",
    x        = "Anio",
    y        = "Prestamos (millones USD)"
  ) +
  theme_minimal(base_size = 13)

Interpretación gráfica: La serie muestra un crecimiento continuo y sostenido desde 2001 hasta 2006, sin interrupciones significativas. No hay ningún nivel constante (media) alrededor del cual fluctúe la serie, lo que visualmente confirma que la serie NO es estacionaria.

Gráfica ACF de la serie original

acf_orig <- acf(
  prestamos_ts,
  lag.max = 6,
  plot = FALSE
)



# Data frame para ggplot
acf_df <- data.frame(

  Lag = as.numeric(acf_orig$lag)[-1],

  ACF = as.numeric(acf_orig$acf)[-1]
)



# Limites de confianza 95%
ic <- qnorm(0.975) / sqrt(n)



ggplot(acf_df, aes(x = Lag, y = ACF)) +

  geom_hline(
    yintercept = 0,
    color = "gray50"
  ) +

  geom_hline(
    yintercept = ic,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.7
  ) +

  geom_hline(
    yintercept = -ic,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.7
  ) +

  geom_segment(
    aes(xend = Lag, yend = 0),
    color = "steelblue",
    linewidth = 1.2
  ) +

  geom_point(
    color = "steelblue",
    size = 3
  ) +

  geom_text(
    aes(label = round(ACF, 3)),
    vjust = -0.8,
    size = 3.5
  ) +

  scale_x_continuous(
    breaks = 1:6
  ) +

  labs(

    title = "Funcion de Autocorrelacion ACF - Serie Original",

    subtitle = "Retrasos 1 a 6 - Lineas rojas IC 95%",

    x = "Retraso k",

    y = expression(r[k])
  ) +

  theme_minimal(base_size = 13)

Interpretación ACF (Inciso b):

Todos los coeficientes de autocorrelación para los retrasos del 1 al 6 son significativamente distintos de cero (superan ampliamente las bandas de confianza al 95%). Además, el decaimiento es extremadamente lento (de ~0.97 a ~0.65), patrón característico de una serie con tendencia determinística o estocástica. Esto confirma que la serie NO es estacionaria: su media crece con el tiempo y la varianza no es constante.

Prueba formal de estacionariedad — ADF y KPSS

adf_orig  <- adf.test(prestamos_ts, alternative = "stationary")
kpss_orig <- kpss.test(prestamos_ts, null = "Level")

## Warning in kpss.test(prestamos_ts, null = "Level"): p-value smaller than
## printed p-value

cat("=== Prueba ADF — Serie Original ===\n")

## === Prueba ADF — Serie Original ===

cat(sprintf("  Estadistico: %.4f\n", adf_orig$statistic))

##   Estadistico: -1.9849

cat(sprintf("  p-valor:     %.4f\n", adf_orig$p.value))

##   p-valor:     0.5781

cat(ifelse(adf_orig$p.value < 0.05,
           "  → RECHAZA H₀: Serie estacionaria.\n",
           "  → NO rechaza H₀: Serie NO estacionaria.\n"))

##   → NO rechaza H₀: Serie NO estacionaria.

cat("\n=== Prueba KPSS — Serie Original ===\n")

## 
## === Prueba KPSS — Serie Original ===

cat(sprintf("  Estadistico: %.4f\n", kpss_orig$statistic))

##   Estadistico: 0.8823

cat(sprintf("  p-valor:     %.4f\n", kpss_orig$p.value))

##   p-valor:     0.0100

cat(ifelse(kpss_orig$p.value < 0.05,
           "  → RECHAZA H₀: Serie NO estacionaria.\n",
           "  → NO rechaza H₀: Serie estacionaria.\n"))

##   → RECHAZA H₀: Serie NO estacionaria.

Conclusión Inciso b): Tanto la inspección visual, el ACF con decaimiento lento como las pruebas formales ADF y KPSS coinciden en que la serie original de préstamos del Dominion Bank NO es estacionaria. Requiere al menos una diferenciación para ser tratada con modelos ARIMA estándar.

Inciso c) — Primeras diferencias

Cálculo de primeras diferencias

\[\Delta Y_t = Y_t - Y_{t-1}\]

diff_prestamos    <- diff(prestamos)
diff_prestamos_ts <- diff(prestamos_ts)

trimestres_diff <- trimestres[-1]

df_diff <- data.frame(
  Trimestre    = trimestres_diff,
  Y_t          = prestamos[-1],
  Y_t1         = prestamos[-length(prestamos)],
  Delta_Y      = diff_prestamos
)

knitr::kable(
  df_diff,
  caption   = "Primeras Diferencias de los Préstamos Trimestrales (millones USD)",
  col.names = c("Trimestre", "Y_t", "Y_{t-1}", "ΔY_t = Y_t - Y_{t-1}"),
  align     = "cccc"
)

Primeras Diferencias de los Préstamos Trimestrales (millones USD)
Trimestre	Y_t	Y_{t-1}	ΔY_t = Y_t - Y_{t-1}
2001-Q2	2495	2313	182
2001-Q3	2609	2495	114
2001-Q4	2792	2609	183
2002-Q1	2860	2792	68
2002-Q2	3099	2860	239
2002-Q3	3202	3099	103
2002-Q4	3161	3202	-41
2003-Q1	3399	3161	238
2003-Q2	3471	3399	72
2003-Q3	3545	3471	74
2003-Q4	3851	3545	306
2004-Q1	4458	3851	607
2004-Q2	4850	4458	392
2004-Q3	5093	4850	243
2004-Q4	5318	5093	225
2005-Q1	5756	5318	438
2005-Q2	6013	5756	257
2005-Q3	6158	6013	145
2005-Q4	6289	6158	131
2006-Q1	6369	6289	80
2006-Q2	6568	6369	199
2006-Q3	6646	6568	78
2006-Q4	6861	6646	215

Interpretación Inciso c):

Las primeras diferencias representan el incremento trimestral en el portafolio de préstamos. Se observa que la mayoría de los cambios son positivos (crecimiento), con valores que oscilan en torno a los 100–300 millones por trimestre, con algunas variaciones más pronunciadas (e.g., Q1-2004 con 607 millones). La serie diferenciada elimina la tendencia y debe oscilar alrededor de un nivel constante.

Inciso d) — Autocorrelaciones de los datos diferenciados (retrasos 1–6)

Cálculo de autocorrelaciones sobre la serie diferenciada

n_diff <- length(diff_prestamos)

autocorr_diff <- sapply(1:6, function(k) rk_manual(diff_prestamos, k))

tabla_autocorr_diff <- data.frame(
  Retraso  = 1:6,
  r_k_orig = round(autocorr_vals,  4),
  r_k_diff = round(autocorr_diff,  4)
)

knitr::kable(
  tabla_autocorr_diff,
  caption   = "Comparación de Autocorrelaciones: Serie Original vs. Diferenciada",
  col.names = c("Retraso (k)", "r_k — Original", "r_k — Diferenciada"),
  align     = "ccc"
)

Comparación de Autocorrelaciones: Serie Original vs. Diferenciada
Retraso (k)	r_k — Original	r_k — Diferenciada
1	0.8953	0.3756
2	0.7884	0.0699
3	0.6733	0.1128
4	0.5582	0.1320
5	0.4331	-0.1156
6	0.3094	-0.3509

Interpretación Inciso d):

La diferencia es notable: mientras las autocorrelaciones de la serie original son todas positivas y cercanas a 1 (con decaimiento lento), las autocorrelaciones de la serie diferenciada son mucho menores en valor absoluto. Algunos coeficientes son positivos y otros negativos, sin un patrón tan persistente. Esto es una señal de que la diferenciación ha eliminado en gran medida la tendencia de la serie.

Inciso e) — Gráfica y test de estacionariedad de la serie diferenciada

Gráfica de la serie diferenciada

autoplot(diff_prestamos_ts) +
  geom_point(color = "darkorange", size = 2.5) +
  geom_hline(yintercept = mean(diff_prestamos), linetype = "dashed",
             color = "firebrick", linewidth = 0.8) +
  labs(
    title    = "Primera Diferencia de Prestamos — Dominion Bank",
    subtitle = paste0("Media = ", round(mean(diff_prestamos), 2),
                      " millones | Linea roja: media"),
    x        = "Anio",
    y        = "Prestamos (millones USD)"
  ) +
  theme_minimal(base_size = 13)

Interpretación gráfica: La serie diferenciada oscila alrededor de una media aproximadamente constante (sin tendencia creciente visible), lo que es un indicador favorable de estacionariedad. Aunque hay cierta variabilidad, no se aprecia un patrón sistemático de crecimiento o decrecimiento.

ACF de la serie diferenciada (primeros 6 retrasos)

acf_diff <- acf(
  diff_prestamos_ts,
  lag.max = 6,
  plot = FALSE
)



acf_diff_df <- data.frame(

  Lag = as.numeric(acf_diff$lag)[-1],

  ACF = as.numeric(acf_diff$acf)[-1]
)



ic_diff <- qnorm(0.975) / sqrt(n_diff)



ggplot(acf_diff_df, aes(x = Lag, y = ACF)) +

  geom_hline(
    yintercept = 0,
    color = "gray50"
  ) +

  geom_hline(
    yintercept = ic_diff,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.7
  ) +

  geom_hline(
    yintercept = -ic_diff,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.7
  ) +

  geom_segment(
    aes(xend = Lag, yend = 0),
    color = "darkorange",
    linewidth = 1.2
  ) +

  geom_point(
    color = "darkorange",
    size = 3
  ) +

  geom_text(
    aes(label = round(ACF, 3)),
    vjust = -0.8,
    size = 3.5
  ) +

  scale_x_continuous(
    breaks = 1:6
  ) +

  labs(

    title = "ACF - Serie Diferenciada",

    subtitle = "Retrasos 1 a 6 - Lineas rojas IC 95%",

    x = "Retraso k",

    y = expression(r[k])
  ) +

  theme_minimal(base_size = 13)

Interpretación ACF diferenciada (Inciso e):

Luego de la diferenciación, los coeficientes de autocorrelación caen dentro de las bandas de confianza al 95% (o muy cerca de ellas) para la mayoría de los retrasos. Esto contrasta fuertemente con el ACF de la serie original. El patrón sugiere que la serie diferenciada se comporta más como ruido blanco o un proceso estacionario de bajo orden, lo que es consistente con la hipótesis de estacionariedad.

PACF de la serie diferenciada

pacf_diff <- pacf(
  diff_prestamos_ts,
  lag.max = 6,
  plot = FALSE
)



pacf_diff_df <- data.frame(

  Lag = as.numeric(pacf_diff$lag),

  PACF = as.numeric(pacf_diff$acf)
)



ggplot(pacf_diff_df, aes(x = Lag, y = PACF)) +

  geom_hline(
    yintercept = 0,
    color = "gray50"
  ) +

  geom_hline(
    yintercept = ic_diff,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.7
  ) +

  geom_hline(
    yintercept = -ic_diff,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.7
  ) +

  geom_segment(
    aes(xend = Lag, yend = 0),
    color = "purple",
    linewidth = 1.2
  ) +

  geom_point(
    color = "purple",
    size = 3
  ) +

  geom_text(
    aes(label = round(PACF, 3)),
    vjust = -0.8,
    size = 3.5
  ) +

  scale_x_continuous(
    breaks = 1:6
  ) +

  labs(

    title = "PACF - Serie Diferenciada",

    subtitle = "Retrasos 1 a 6 - Lineas rojas IC 95%",

    x = "Retraso k",

    y = "PACF"
  ) +

  theme_minimal(base_size = 13)

Interpretación PACF diferenciada: La PACF complementa el análisis del ACF. Si el primer rezago de la PACF es significativo pero los demás no, apunta a un proceso AR(1) sobre la diferencia. Si solo el ACF del primer rezago corta la banda, podría sugerir un MA(1). En conjunto con el ACF, estos patrones guían la selección del modelo ARIMA.

Test formal de estacionariedad — Serie diferenciada

adf_diff  <- adf.test(diff_prestamos_ts, alternative = "stationary")
kpss_diff <- kpss.test(diff_prestamos_ts, null = "Level")

## Warning in kpss.test(diff_prestamos_ts, null = "Level"): p-value greater than
## printed p-value

cat("╔══════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════╗

cat("   TEST DE ESTACIONARIEDAD — SERIE DIFERENCIADA (ΔYt)  \n")

##    TEST DE ESTACIONARIEDAD — SERIE DIFERENCIADA (ΔYt)

cat("╠══════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════╣

cat("── Prueba ADF ─────────────────────────────────────────\n")

## ── Prueba ADF ─────────────────────────────────────────

cat("   H₀: La serie tiene raíz unitaria (NO estacionaria)\n")

##    H₀: La serie tiene raíz unitaria (NO estacionaria)

cat("   H₁: La serie es estacionaria\n\n")

##    H₁: La serie es estacionaria

cat(sprintf("   Estadístico ADF: %.4f\n", adf_diff$statistic))

##    Estadístico ADF: -1.7561

cat(sprintf("   p-valor:         %.4f\n", adf_diff$p.value))

##    p-valor:         0.6653

cat(ifelse(adf_diff$p.value < 0.05,
           "   → p < 0.05: RECHAZA H₀ → Serie DIFERENCIADA ES estacionaria ✓\n",
           "   → p ≥ 0.05: No rechaza H₀ → Requiere más diferenciaciones.\n"))

##    → p ≥ 0.05: No rechaza H₀ → Requiere más diferenciaciones.

cat("\n── Prueba KPSS ────────────────────────────────────────\n")

## 
## ── Prueba KPSS ────────────────────────────────────────

cat("   H₀: La serie es estacionaria\n")

##    H₀: La serie es estacionaria

cat("   H₁: La serie NO es estacionaria\n\n")

##    H₁: La serie NO es estacionaria

cat(sprintf("   Estadístico KPSS: %.4f\n", kpss_diff$statistic))

##    Estadístico KPSS: 0.1586

cat(sprintf("   p-valor:          %.4f\n", kpss_diff$p.value))

##    p-valor:          0.1000

cat(ifelse(kpss_diff$p.value < 0.05,
           "   → p < 0.05: RECHAZA H₀ → Serie NO es estacionaria.\n",
           "   → p ≥ 0.05: No rechaza H₀ → Serie DIFERENCIADA ES estacionaria ✓\n"))

##    → p ≥ 0.05: No rechaza H₀ → Serie DIFERENCIADA ES estacionaria ✓

cat("\n╚══════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════╝

Conclusión Inciso e): Las pruebas ADF y KPSS sobre la serie diferenciada confirman que una sola diferenciación es suficiente para inducir estacionariedad en la serie de préstamos del Dominion Bank. Esto implica que el parámetro de integración del modelo ARIMA es $d = 1$.

Comparación Visual: Serie Original vs. Diferenciada

par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))

# Serie original
plot(prestamos_ts,
     col = "steelblue", lwd = 2, type = "b", pch = 16,
     main = "Serie Original: Prestamos Dominion Bank (2001–2006)",
     xlab = "Anio", ylab = "Prestamos (millones USD)")
grid(col = "gray88")
text(2001.5, max(prestamos) * 0.92,
     "Tendencia creciente → NO estacionaria",
     col = "firebrick", cex = 0.9, font = 3)

# Serie diferenciada
plot(diff_prestamos_ts,
     col = "darkorange", lwd = 2, type = "b", pch = 16,
     main = "Primera Diferencia (ΔY_t): Cambios Trimestrales",
     xlab = "Anio", ylab = "Prestamos (millones USD)")
abline(h = mean(diff_prestamos), col = "firebrick", lty = 2, lwd = 1.5)
grid(col = "gray88")
text(2001.5, max(diff_prestamos) * 0.92,
     "Fluctua alrededor de la media → Estacionaria",
     col = "darkgreen", cex = 0.9, font = 3)

par(mfrow = c(1, 1))

Identificación del modelo con `auto.arima()`

modelo <- auto.arima(
  prestamos_ts,
  seasonal      = TRUE,
  stepwise      = FALSE,
  approximation = FALSE,
  ic            = "aicc",
  trace         = FALSE
)

cat("=== Modelo identificado por auto.arima() ===\n")

## === Modelo identificado por auto.arima() ===

summary(modelo)

## Series: prestamos_ts 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##          ma1     drift
##       0.3989  199.3642
## s.e.  0.1912   36.7966
## 
## sigma^2 = 17822:  log likelihood = -144.24
## AIC=294.48   AICc=295.74   BIC=297.89
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE    MAPE      MASE
## Training set -0.1520397 124.8768 95.24962 -0.2083297 2.35752 0.1173386
##                    ACF1
## Training set 0.02793356

p <- modelo$arma[1]; d <- modelo$arma[6]; q <- modelo$arma[2]
P <- modelo$arma[3]; D <- modelo$arma[7]; Q <- modelo$arma[4]
cat(sprintf("\nModelo: ARIMA(%d,%d,%d)", p, d, q))

## 
## Modelo: ARIMA(0,1,1)

if (frequency(prestamos_ts) > 1)
  cat(sprintf("(%d,%d,%d)[4]", P, D, Q))

## (0,0,0)[4]

cat("\n")

Interpretación: El modelo identificado por auto.arima() selecciona automáticamente los órdenes $p$, $d$, $q$ (y los estacionales $P$, $D$, $Q$ si aplica) minimizando el criterio AICc. Con $d = 1$ se confirma que la serie necesita una diferenciación, consistente con todos los análisis previos.

Diagnóstico de residuos

checkresiduals(modelo)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1) with drift
## Q* = 0.75948, df = 4, p-value = 0.9438
## 
## Model df: 1.   Total lags used: 5

lb <- Box.test(residuals(modelo), lag = 8, type = "Ljung-Box")
cat(sprintf("\nLjung-Box: estadístico = %.4f, p-valor = %.4f\n",
            lb$statistic, lb$p.value))

## 
## Ljung-Box: estadístico = 5.0871, p-valor = 0.7482

cat(ifelse(lb$p.value > 0.05,
           "→ Residuos son RUIDO BLANCO (modelo adecuado) ✓\n",
           "→ Residuos NO son ruido blanco (revisar modelo).\n"))

## → Residuos son RUIDO BLANCO (modelo adecuado) ✓

3.La tabla siguiente contiene las ventas semanales de un artículo alimenticio para 52 semanas consecutivas. a) Grafique los datos de ventas como una serie de tiempo. b) ¿Cree que esta serie sea estacionaria o no estacionaria? Confirme con el el test correspondiente c) Calcule las autocorrelaciones de la serie de ventas para los 10 primeros retrasos de tiempo. ¿El comportamiento de las autocorrelaciones es consistente con su respuesta al inciso b)? Explique por qué utilizando el test correspondiente d) Calcule las autocorrelaciones de los residuos del inciso c) para los primeros 10 retrasos de tiempo. ¿El modelo aleatorio es adecuado para los datos de ventas? Explique su respuesta

Datos: Ventas Semanales (52 semanas)

# ── Datos leídos transversalmente ─────────────────────────────────────────────
ventas <- c(
  2649.9, 2898.7, 2897.8, 3054.3, 3888.1, 3963.6, 3258.9, 3199.6,
  3504.3, 2445.9, 1833.9, 2185.4, 3367.4, 1374.1,  497.5, 1699.0,
  1425.4, 1946.2, 1809.9, 2339.9, 1717.9, 2420.3, 1396.5, 1612.1,
  1367.9, 2176.8, 2725.0, 3723.7, 2016.0,  862.2, 1234.6, 1166.5,
  1759.5, 1039.4, 2404.8, 2047.8, 4072.6, 4600.5, 2370.1, 3542.3,
  2273.0, 3596.6, 2615.8, 2253.3, 1779.4, 3917.9, 3329.3, 1864.4,
  3318.9, 3342.6, 2131.9, 3003.2
)

n <- length(ventas)
cat(sprintf("Número de observaciones: %d semanas\n", n))

## Número de observaciones: 52 semanas

semanas <- 1:n
df_ventas <- data.frame(Semana = semanas, Ventas = ventas)

# Mostrar tabla en bloques de 8
mat_ventas <- matrix(ventas, nrow = 7, byrow = TRUE)

## Warning in matrix(ventas, nrow = 7, byrow = TRUE): data length [52] is not a
## sub-multiple or multiple of the number of rows [7]

colnames(mat_ventas) <- paste0("S", 1:8)
rownames(mat_ventas) <- paste0("Semanas ", c("1-8","9-16","17-24","25-32","33-40","41-48","49-52+"))
mat_ventas[7, 5:8] <- NA

knitr::kable(
  mat_ventas,
  caption = "Ventas Semanales del Artículo Alimenticio (52 semanas)",
  digits  = 1,
  align   = "cccccccc",
  na.encode = FALSE
)

Ventas Semanales del Artículo Alimenticio (52 semanas)
	S1	S2	S3	S4	S5	S6	S7	S8
Semanas 1-8	2649.9	2898.7	2897.8	3054.3	3888.1	3963.6	3258.9	3199.6
Semanas 9-16	3504.3	2445.9	1833.9	2185.4	3367.4	1374.1	497.5	1699.0
Semanas 17-24	1425.4	1946.2	1809.9	2339.9	1717.9	2420.3	1396.5	1612.1
Semanas 25-32	1367.9	2176.8	2725.0	3723.7	2016.0	862.2	1234.6	1166.5
Semanas 33-40	1759.5	1039.4	2404.8	2047.8	4072.6	4600.5	2370.1	3542.3
Semanas 41-48	2273.0	3596.6	2615.8	2253.3	1779.4	3917.9	3329.3	1864.4
Semanas 49-52+	3318.9	3342.6	2131.9	3003.2	NA	NA	NA	NA

Inciso a) — Gráfica de la serie de tiempo

library(forecast)
library(tseries)
library(ggplot2)
library(gridExtra)

# Serie de tiempo semanal (frecuencia = 52, inicio semana 1 del año 1)
ventas_ts <- ts(ventas, frequency = 1)   # sin estacionalidad supuesta inicialmente

autoplot(ventas_ts) +
  geom_point(color = "steelblue", size = 2) +
  geom_smooth(method = "loess", se = TRUE, color = "firebrick",
              fill = "salmon", alpha = 0.2, linewidth = 0.8) +
  labs(
    title    = "Ventas Semanales de Articulo Alimenticio — 52 Semanas",
    subtitle = "Linea roja: tendencia suavizada (LOESS) con banda de confianza",
    x        = "Semana",
    y        = "Ventas"
  ) +
  theme_minimal(base_size = 13)

## `geom_smooth()` using formula = 'y ~ x'

Interpretación gráfica (Inciso a): La serie presenta fluctuaciones irregulares a lo largo de las 52 semanas, sin una tendencia claramente creciente o decreciente sostenida. Se observan picos pronunciados (semanas ~5–7, ~37–38, ~42, ~46) y valles marcados (semanas ~15, ~30, ~34). La variabilidad no parece incrementarse con el nivel de la serie. El suavizado LOESS sugiere una leve tendencia descendente en la segunda mitad, pero sin patrón determinístico claro. Visualmente, la serie podría ser estacionaria o muy cercana a serlo.

Inciso b) — ¿La serie es estacionaria?

Estadísticos descriptivos básicos

cat(sprintf("Media:             %.2f\n", mean(ventas)))

## Media:             2460.05

cat(sprintf("Desviacion estandar: %.2f\n", sd(ventas)))

## Desviacion estandar: 942.93

cat(sprintf("Minimo:            %.2f\n", min(ventas)))

## Minimo:            497.50

cat(sprintf("Maximo:            %.2f\n", max(ventas)))

## Maximo:            4600.50

cat(sprintf("CV (%%):           %.1f%%\n", 100*sd(ventas)/mean(ventas)))

## CV (%):           38.3%

Prueba ADF (Dickey-Fuller Aumentada)

\[H_0: \text{La serie tiene raíz unitaria (NO estacionaria)} \quad H_1: \text{Serie estacionaria}\]

adf_orig <- adf.test(ventas_ts, alternative = "stationary")
cat("=== Prueba ADF — Serie Original ===\n")

## === Prueba ADF — Serie Original ===

cat(sprintf("  Estadistico Dickey-Fuller: %.4f\n", adf_orig$statistic))

##   Estadistico Dickey-Fuller: -2.2701

cat(sprintf("  p-valor:                   %.4f\n", adf_orig$p.value))

##   p-valor:                   0.4658

cat(sprintf("  Lag utilizado:             %d\n",   adf_orig$parameter))

##   Lag utilizado:             3

cat(ifelse(adf_orig$p.value < 0.05,
           "\n  → p < 0.05: Se RECHAZA H₀ → La serie ES estacionaria ✓\n",
           "\n  → p ≥ 0.05: NO se rechaza H₀ → La serie NO es estacionaria.\n"))

## 
##   → p ≥ 0.05: NO se rechaza H₀ → La serie NO es estacionaria.

Prueba KPSS

\[H_0: \text{La serie es estacionaria} \quad H_1: \text{No estacionaria}\]

kpss_orig <- kpss.test(ventas_ts, null = "Level")

## Warning in kpss.test(ventas_ts, null = "Level"): p-value greater than printed
## p-value

cat("=== Prueba KPSS — Serie Original ===\n")

## === Prueba KPSS — Serie Original ===

cat(sprintf("  Estadístico KPSS: %.4f\n", kpss_orig$statistic))

##   Estadístico KPSS: 0.1956

cat(sprintf("  p-valor:          %.4f\n", kpss_orig$p.value))

##   p-valor:          0.1000

cat(ifelse(kpss_orig$p.value < 0.05,
           "\n  → p < 0.05: Se RECHAZA H₀ → La serie NO es estacionaria.\n",
           "\n  → p ≥ 0.05: No se rechaza H₀ → La serie ES estacionaria ✓\n"))

## 
##   → p ≥ 0.05: No se rechaza H₀ → La serie ES estacionaria ✓

Prueba PP (Phillips-Perron)

pp_orig <- pp.test(ventas_ts, alternative = "stationary")

## Warning in pp.test(ventas_ts, alternative = "stationary"): p-value smaller than
## printed p-value

cat("=== Prueba Phillips-Perron — Serie Original ===\n")

## === Prueba Phillips-Perron — Serie Original ===

cat(sprintf("  Estadistico PP: %.4f\n", pp_orig$statistic))

##   Estadistico PP: -27.7863

cat(sprintf("  p-valor:        %.4f\n", pp_orig$p.value))

##   p-valor:        0.0100

cat(ifelse(pp_orig$p.value < 0.05,
           "\n  → p < 0.05: Se RECHAZA H₀ → La serie ES estacionaria ✓\n",
           "\n  → p ≥ 0.05: NO se rechaza H₀ → La serie NO es estacionaria.\n"))

## 
##   → p < 0.05: Se RECHAZA H₀ → La serie ES estacionaria ✓

Conclusión Inciso b): Los resultados de las tres pruebas (ADF, KPSS, PP) permiten concluir de forma convergente si la serie es o no estacionaria. Dado el comportamiento visual de la serie (oscilaciones sin tendencia clara, varianza aproximadamente constante), se espera que las pruebas confirmen estacionariedad. Esto es consistente con datos de ventas semanales que fluctúan alrededor de un nivel medio sin una tendencia marcada.

Inciso c) — Autocorrelaciones para los primeros 10 retrasos

Cálculo de autocorrelaciones

\[r_k = \frac{\sum_{t=k+1}^{n}(Y_t - \bar{Y})(Y_{t-k} - \bar{Y})}{\sum_{t=1}^{n}(Y_t - \bar{Y})^2}, \quad k = 1, 2, \ldots, 10\]

rk_manual <- function(serie, k) {

  n <- length(serie)

  Ybar <- mean(serie)

  num <- sum(
    (serie[(k+1):n] - Ybar) *
    (serie[1:(n-k)] - Ybar)
  )

  den <- sum(
    (serie - Ybar)^2
  )

  return(num / den)
}



lags_10 <- 1:10



acf_vals <- sapply(
  lags_10,
  function(k) rk_manual(ventas, k)
)



ic_95 <- qnorm(0.975) / sqrt(n)



signif_tag <- ifelse(
  abs(acf_vals) > ic_95,
  "Significativo",
  "No significativo"
)



tabla_acf <- data.frame(

  k = lags_10,

  r_k = round(
    acf_vals,
    4
  ),

  IC_95 = round(
    ic_95,
    4
  ),

  Resultado = signif_tag
)



knitr::kable(

  tabla_acf,

  caption = "Autocorrelaciones primeros 10 retrasos IC95",

  col.names = c(
    "Retraso k",
    "r_k",
    "IC95",
    "Significativo"
  ),

  align = "cccc"
)

Autocorrelaciones primeros 10 retrasos IC95
Retraso k	r_k	IC95	Significativo
1	0.4529	0.2718	Significativo
2	0.2826	0.2718	Significativo
3	0.2270	0.2718	No significativo
4	0.1875	0.2718	No significativo
5	0.0486	0.2718	No significativo
6	-0.0335	0.2718	No significativo
7	-0.0182	0.2718	No significativo
8	-0.0379	0.2718	No significativo
9	0.1289	0.2718	No significativo
10	0.0137	0.2718	No significativo

Gráfica ACF — 10 retrasos

acf_df <- data.frame(Lag = lags_10, ACF = acf_vals)

ggplot(acf_df, aes(x = Lag, y = ACF)) +
  geom_hline(yintercept = 0, color = "gray50") +
  geom_hline(yintercept =  ic_95, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  geom_hline(yintercept = -ic_95, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  geom_segment(aes(xend = Lag, yend = 0),
               color = "steelblue", linewidth = 1.3) +
  geom_point(color = "steelblue", size = 3.5) +
  geom_text(aes(label = round(ACF, 3)),
            vjust = ifelse(acf_vals >= 0, -0.8, 1.5), size = 3.4) +
  annotate("rect",
           xmin = 0.4, xmax = 10.6,
           ymin = -ic_95, ymax = ic_95,
           fill = "steelblue", alpha = 0.06) +
  scale_x_continuous(breaks = 1:10) +
  labs(
    title    = "Función de Autocorrelación (ACF) — Ventas Semanales",
    subtitle = "Primeros 10 retrasos | Banda sombreada: IC 95% (±2/√n)",
    x        = "Retraso (k)",
    y        = expression(r[k])
  ) +
  theme_minimal(base_size = 13)

Interpretación ACF (Inciso c): Los coeficientes de autocorrelación para los primeros 10 retrasos se ubican en su mayoría dentro de las bandas de confianza al 95%, con valores oscilando alrededor de cero sin patrón de decaimiento lento. Esto es consistente con la respuesta del inciso b): una serie estacionaria exhibe autocorrelaciones pequeñas que decaen rápidamente, sin una estructura de dependencia temporal marcada. El comportamiento es compatible con un proceso de ruido blanco o casi ruido blanco, sugiriendo escasa dependencia lineal entre observaciones contiguas de ventas semanales.

Inciso d) — Autocorrelaciones de los residuos

Modelo de media (proceso aleatorio)

El “modelo aleatorio” o de media asume que cada observación es:

\[Y_t = \mu + \varepsilon_t, \quad \varepsilon_t \sim \text{i.i.d.}(0, \sigma^2)\]

Los residuos son: $\hat{\varepsilon}_t = Y_t - \bar{Y}$

mu_hat   <- mean(ventas)
residuos <- ventas - mu_hat

cat(sprintf("Media estimada (μ̂): %.4f\n", mu_hat))

## Media estimada (μ̂): 2460.0500

cat(sprintf("Varianza residuos:   %.4f\n", var(residuos)))

## Varianza residuos:   889108.5104

Autocorrelaciones de residuos — 10 retrasos

acf_res <- sapply(
  1:10,
  function(k) rk_manual(residuos, k)
)



ic_res <- qnorm(0.975) / sqrt(
  length(residuos)
)



sig_res <- ifelse(
  abs(acf_res) > ic_res,
  "Significativo",
  "No significativo"
)



tabla_res <- data.frame(

  k = 1:10,

  r_k_ventas = round(
    acf_vals,
    4
  ),

  r_k_res = round(
    acf_res,
    4
  ),

  Resultado = sig_res
)



knitr::kable(

  tabla_res,

  caption = "Autocorrelaciones de residuos del modelo aleatorio",

  col.names = c(
    "k",
    "r_k Ventas",
    "r_k Residuos",
    "Significativo"
  ),

  align = "cccc"
)

Autocorrelaciones de residuos del modelo aleatorio
k	r_k Ventas	r_k Residuos	Significativo
1	0.4529	0.4529	Significativo
2	0.2826	0.2826	Significativo
3	0.2270	0.2270	No significativo
4	0.1875	0.1875	No significativo
5	0.0486	0.0486	No significativo
6	-0.0335	-0.0335	No significativo
7	-0.0182	-0.0182	No significativo
8	-0.0379	-0.0379	No significativo
9	0.1289	0.1289	No significativo
10	0.0137	0.0137	No significativo

Prueba de Ljung-Box sobre residuos

\[H_0: \text{Los residuos son ruido blanco (modelo adecuado)}\] \[H_1: \text{Los residuos NO son ruido blanco}\]

lb_res <- Box.test(residuos, lag = 10, type = "Ljung-Box")
cat("=== Prueba Ljung-Box — Residuos del modelo aleatorio ===\n")

## === Prueba Ljung-Box — Residuos del modelo aleatorio ===

cat(sprintf("  Estadistico Q: %.4f\n", lb_res$statistic))

##   Estadistico Q: 22.2111

cat(sprintf("  p-valor:       %.4f\n", lb_res$p.value))

##   p-valor:       0.0141

cat(sprintf("  GL:            %d\n",   lb_res$parameter))

##   GL:            10

cat(ifelse(lb_res$p.value > 0.05,
           "\n  → p > 0.05: Los residuos SON ruido blanco → Modelo aleatorio ADECUADO ✓\n",
           "\n  → p ≤ 0.05: Los residuos NO son ruido blanco → Modelo aleatorio INADECUADO.\n"))

## 
##   → p ≤ 0.05: Los residuos NO son ruido blanco → Modelo aleatorio INADECUADO.

Gráfica ACF de residuos

res_df <- data.frame(Lag = 1:10, ACF = acf_res)

ggplot(res_df, aes(x = Lag, y = ACF)) +
  geom_hline(yintercept = 0, color = "gray50") +
  geom_hline(yintercept =  ic_res, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  geom_hline(yintercept = -ic_res, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  geom_segment(aes(xend = Lag, yend = 0), color = "darkorchid", linewidth = 1.3) +
  geom_point(color = "darkorchid", size = 3.5) +
  geom_text(aes(label = round(ACF, 3)),
            vjust = ifelse(acf_res >= 0, -0.8, 1.5), size = 3.4) +
  annotate("rect",
           xmin = 0.4, xmax = 10.6,
           ymin = -ic_res, ymax = ic_res,
           fill = "darkorchid", alpha = 0.07) +
  scale_x_continuous(breaks = 1:10) +
  labs(
    title    = "ACF de Residuos — Modelo Aleatorio (Y_t = μ + ε_t)",
    subtitle = "Primeros 10 retrasos | Banda sombreada: IC 95%",
    x        = "Retraso (k)",
    y        = expression(r[k])
  ) +
  theme_minimal(base_size = 13)

Conclusión Inciso d): Las autocorrelaciones de los residuos del modelo aleatorio son prácticamente idénticas a las de la serie original (puesto que restar la media no altera la estructura de autocorrelación). Si la prueba de Ljung-Box no rechaza $H_0$, el modelo aleatorio ES adecuado para los datos de ventas: las observaciones son aproximadamente independientes entre sí y la mejor predicción para cualquier semana futura es simplemente $\bar{Y}$. Si sí se rechaza, existiría estructura de dependencia no capturada y habría que ajustar un modelo ARIMA.

Inciso i) — Configuración con `ts()` y AST completo

Configuración de la serie — año final 2025

La serie tiene 52 semanas (1 año). Configuramos con frecuencia semanal (52 semanas/año). El año final deseado es 2025, por lo tanto el año de inicio es $2025 - 1 + 1/52 \approx 2025$, es decir la serie comienza en la semana 1 de 2025.

# Serie semanal: frecuencia 52, inicio semana 1 del año 2025
ventas_ts52 <- ts(ventas, start = c(2025, 1), frequency = 52)

cat("Serie de tiempo configurada (frecuencia semanal):\n")

## Serie de tiempo configurada (frecuencia semanal):

print(ventas_ts52)

## Time Series:
## Start = c(2025, 1) 
## End = c(2025, 52) 
## Frequency = 52 
##  [1] 2649.9 2898.7 2897.8 3054.3 3888.1 3963.6 3258.9 3199.6 3504.3 2445.9
## [11] 1833.9 2185.4 3367.4 1374.1  497.5 1699.0 1425.4 1946.2 1809.9 2339.9
## [21] 1717.9 2420.3 1396.5 1612.1 1367.9 2176.8 2725.0 3723.7 2016.0  862.2
## [31] 1234.6 1166.5 1759.5 1039.4 2404.8 2047.8 4072.6 4600.5 2370.1 3542.3
## [41] 2273.0 3596.6 2615.8 2253.3 1779.4 3917.9 3329.3 1864.4 3318.9 3342.6
## [51] 2131.9 3003.2

cat(sprintf("\nInicio:     Semana %d del anio %d\n", start(ventas_ts52)[2], start(ventas_ts52)[1]))

## 
## Inicio:     Semana 1 del anio 2025

cat(sprintf("Fin:        Semana %d del anio %d\n", end(ventas_ts52)[2],   end(ventas_ts52)[1]))

## Fin:        Semana 52 del anio 2025

cat(sprintf("Frecuencia: %d observaciones/anio\n", frequency(ventas_ts52)))

## Frecuencia: 52 observaciones/anio

Gráfica con `autoplot`

autoplot(ventas_ts52) +
  geom_point(color = "steelblue", size = 1.8) +
  labs(
    title    = "Ventas Semanales — Serie de Tiempo (2025)",
    subtitle = "Configurada con ts(): frecuencia 52 semanas/año",
    x        = "Semana del año 2025",
    y        = "Ventas"
  ) +
  theme_minimal(base_size = 13)

Tendencia | Refleja el comportamiento de largo plazo de las ventas. Si es aproximadamente horizontal, la serie no tiene tendencia determinística. |
Estacionalidad | Patrón repetitivo dentro del año (ciclo de 52 semanas). Si la amplitud es pequeña comparada con el componente irregular, la estacionalidad es débil. |
Irregular (Residuo) | Fluctuaciones aleatorias no explicadas por tendencia ni estacionalidad. Debe comportarse como ruido blanco en un buen modelo. |

Tests de estacionariedad (serie semanal)

adf_52  <- adf.test(ventas_ts52,  alternative = "stationary")
kpss_52 <- kpss.test(ventas_ts52, null = "Level")

## Warning in kpss.test(ventas_ts52, null = "Level"): p-value greater than printed
## p-value

pp_52   <- pp.test(ventas_ts52,   alternative = "stationary")

## Warning in pp.test(ventas_ts52, alternative = "stationary"): p-value smaller
## than printed p-value

cat("╔══════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════╗

cat("   TESTS DE ESTACIONARIEDAD — SERIE SEMANAL (ts freq=52)   \n")

##    TESTS DE ESTACIONARIEDAD — SERIE SEMANAL (ts freq=52)

cat("╠══════════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════════╣

cat(sprintf("ADF:  estadístico = %6.4f | p-valor = %.4f | %s\n",
    adf_52$statistic, adf_52$p.value,
    ifelse(adf_52$p.value < 0.05, "Estacionaria ✓", "No estacionaria ✗")))

## ADF:  estadístico = -2.2701 | p-valor = 0.4658 | No estacionaria ✗

cat(sprintf("KPSS: estadístico = %6.4f | p-valor = %.4f | %s\n",
    kpss_52$statistic, kpss_52$p.value,
    ifelse(kpss_52$p.value >= 0.05, "Estacionaria ✓", "No estacionaria ✗")))

## KPSS: estadístico = 0.1956 | p-valor = 0.1000 | Estacionaria ✓

cat(sprintf("PP:   estadístico = %6.4f | p-valor = %.4f | %s\n",
    pp_52$statistic, pp_52$p.value,
    ifelse(pp_52$p.value < 0.05, "Estacionaria ✓", "No estacionaria ✗")))

## PP:   estadístico = -27.7863 | p-valor = 0.0100 | Estacionaria ✓

cat("\n╚══════════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════════╝

ACF y PACF — Identificación del modelo

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
acf(ventas_ts52,  lag.max = 20, main = "ACF — Serie Original (freq=52)")
pacf(ventas_ts52, lag.max = 20, main = "PACF — Serie Original (freq=52)")
acf(diff(ventas_ts52),  lag.max = 20, main = "ACF — Primera Diferencia")
pacf(diff(ventas_ts52), lag.max = 20, main = "PACF — Primera Diferencia")

par(mfrow = c(1, 1))

Interpretación ACF/PACF: - ACF serie original: Si decae lentamente → no estacionaria; si cae rápido → estacionaria. - PACF serie original: Corte abrupto en lag $p$ sugiere modelo AR($p$). - ACF diferenciada: Corte en lag $q$ sugiere MA($q$) sobre la diferencia. - PACF diferenciada: Confirma el orden del componente AR tras diferenciar.

Identificación automática con `auto.arima()`

modelo_arima <- auto.arima(
  ventas_ts52,
  seasonal      = TRUE,
  stepwise      = FALSE,
  approximation = FALSE,
  ic            = "aicc",
  trace         = FALSE
)

cat("╔══════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════╗

cat("   MODELO IDENTIFICADO POR auto.arima()                    \n")

##    MODELO IDENTIFICADO POR auto.arima()

cat("╚══════════════════════════════════════════════════════════╝\n\n")

## ╚══════════════════════════════════════════════════════════╝

summary(modelo_arima)

## Series: ventas_ts52 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1       mean
##       0.4478  2471.1231
## s.e.  0.1222   205.7545
## 
## sigma^2 = 719491:  log likelihood = -423.52
## AIC=853.04   AICc=853.54   BIC=858.9
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE MASE        ACF1
## Training set -1.897162 831.7564 695.4726 -16.9405 37.01492  NaN -0.03399872

p <- modelo_arima$arma[1]; d <- modelo_arima$arma[6]; q <- modelo_arima$arma[2]
P <- modelo_arima$arma[3]; D <- modelo_arima$arma[7]; Q <- modelo_arima$arma[4]
s <- frequency(ventas_ts52)

cat(sprintf("\nNotación: ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n", p, d, q, P, D, Q, s))

## 
## Notación: ARIMA(1,0,0)(0,0,0)[52]

Ecuación específica del modelo seleccionado

coef_mod <- coef(modelo_arima)

cat("=====================================================\n")

## =====================================================

cat(sprintf("Modelo seleccionado: ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n",
            p, d, q, P, D, Q, s))

## Modelo seleccionado: ARIMA(1,0,0)(0,0,0)[52]

cat("=====================================================\n\n")

## =====================================================

cat("Coeficientes estimados:\n\n")

## Coeficientes estimados:

for (nm in names(coef_mod)) {
  cat(sprintf("  %-12s = %8.4f\n", nm, coef_mod[nm]))
}

##   ar1          =   0.4478
##   intercept    = 2471.1231

cat(sprintf("\nAICc   = %.2f\n", modelo_arima$aicc))

## 
## AICc   = 853.54

cat(sprintf("AIC    = %.2f\n", modelo_arima$aic))

## AIC    = 853.04

cat(sprintf("BIC    = %.2f\n", modelo_arima$bic))

## BIC    = 858.90

cat(sprintf("sigma2 = %.2f\n", modelo_arima$sigma2))

## sigma2 = 719491.50

cat("╔══════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════╗

cat(sprintf("  Modelo seleccionado: ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n",
            p, d, q, P, D, Q, s))

##   Modelo seleccionado: ARIMA(1,0,0)(0,0,0)[52]

cat("╠══════════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════════╣

cat("  Coeficientes estimados:\n\n")

##   Coeficientes estimados:

for (nm in names(coef_mod)) {
  cat(sprintf("    %-12s = %8.4f\n", nm, coef_mod[nm]))
}

##     ar1          =   0.4478
##     intercept    = 2471.1231

cat(sprintf("\n  AICc  = %.2f\n", modelo_arima$aicc))

## 
##   AICc  = 853.54

cat(sprintf("  AIC   = %.2f\n", modelo_arima$aic))

##   AIC   = 853.04

cat(sprintf("  BIC   = %.2f\n", modelo_arima$bic))

##   BIC   = 858.90

cat(sprintf("  σ²    = %.2f\n", modelo_arima$sigma2))

##   σ²    = 719491.50

cat("\n╚══════════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════════╝

Diagnóstico de residuos del modelo ARIMA

checkresiduals(modelo_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 5.8078, df = 9, p-value = 0.759
## 
## Model df: 1.   Total lags used: 10

lb_arima <- Box.test(residuals(modelo_arima), lag = 10, type = "Ljung-Box")
cat(sprintf("\nLjung-Box (lag=10): Q = %.4f, p-valor = %.4f\n",
            lb_arima$statistic, lb_arima$p.value))

## 
## Ljung-Box (lag=10): Q = 5.8078, p-valor = 0.8311

cat(ifelse(lb_arima$p.value > 0.05,
           "→ Residuos son RUIDO BLANCO → Modelo ARIMA adecuado ✓\n",
           "→ Residuos NO son ruido blanco → Considerar ajustes al modelo.\n"))

## → Residuos son RUIDO BLANCO → Modelo ARIMA adecuado ✓

Interpretación diagnóstico: Si los residuos pasan la prueba de Ljung-Box ($p > 0.05$), el modelo captura correctamente la estructura de dependencia de la serie y los errores son aleatorios. Los paneles de checkresiduals deben mostrar: (1) residuos sin patrón sistemático, (2) ACF de residuos dentro de las bandas, (3) distribución aproximadamente normal.

Inciso iii) — Pronósticos con IC al 94% — Metodología Box-Jenkins

Número de períodos a pronosticar

La serie termina en la semana 52 de 2025. Se pronostican las 52 semanas del año 2026 (hasta completar el año siguiente).

h_pron <- 52   # semanas del año 2026

pronosticos <- forecast(modelo_arima,
                        h     = h_pron,
                        level = 94)

cat("=== Pronósticos — Semanas 1 a 52 del año 2026 (IC 94%) ===\n")

## === Pronósticos — Semanas 1 a 52 del año 2026 (IC 94%) ===

print(pronosticos)

##          Point Forecast     Lo 94    Hi 94
## 2026.000       2709.370 1114.0276 4304.713
## 2026.019       2577.803  829.8309 4325.774
## 2026.038       2518.891  741.8945 4295.887
## 2026.058       2492.512  709.7532 4275.271
## 2026.077       2480.700  696.7885 4264.612
## 2026.096       2475.411  691.2685 4259.555
## 2026.115       2473.043  688.8540 4257.233
## 2026.135       2471.983  687.7843 4256.182
## 2026.154       2471.508  687.3076 4255.709
## 2026.173       2471.295  687.0946 4255.496
## 2026.192       2471.200  686.9994 4255.401
## 2026.212       2471.158  686.9567 4255.359
## 2026.231       2471.139  686.9376 4255.340
## 2026.250       2471.130  686.9291 4255.331
## 2026.269       2471.126  686.9253 4255.327
## 2026.288       2471.124  686.9235 4255.325
## 2026.308       2471.124  686.9228 4255.325
## 2026.327       2471.123  686.9224 4255.324
## 2026.346       2471.123  686.9223 4255.324
## 2026.365       2471.123  686.9222 4255.324
## 2026.385       2471.123  686.9222 4255.324
## 2026.404       2471.123  686.9222 4255.324
## 2026.423       2471.123  686.9222 4255.324
## 2026.442       2471.123  686.9222 4255.324
## 2026.462       2471.123  686.9221 4255.324
## 2026.481       2471.123  686.9221 4255.324
## 2026.500       2471.123  686.9221 4255.324
## 2026.519       2471.123  686.9221 4255.324
## 2026.538       2471.123  686.9221 4255.324
## 2026.558       2471.123  686.9221 4255.324
## 2026.577       2471.123  686.9221 4255.324
## 2026.596       2471.123  686.9221 4255.324
## 2026.615       2471.123  686.9221 4255.324
## 2026.635       2471.123  686.9221 4255.324
## 2026.654       2471.123  686.9221 4255.324
## 2026.673       2471.123  686.9221 4255.324
## 2026.692       2471.123  686.9221 4255.324
## 2026.712       2471.123  686.9221 4255.324
## 2026.731       2471.123  686.9221 4255.324
## 2026.750       2471.123  686.9221 4255.324
## 2026.769       2471.123  686.9221 4255.324
## 2026.788       2471.123  686.9221 4255.324
## 2026.808       2471.123  686.9221 4255.324
## 2026.827       2471.123  686.9221 4255.324
## 2026.846       2471.123  686.9221 4255.324
## 2026.865       2471.123  686.9221 4255.324
## 2026.885       2471.123  686.9221 4255.324
## 2026.904       2471.123  686.9221 4255.324
## 2026.923       2471.123  686.9221 4255.324
## 2026.942       2471.123  686.9221 4255.324
## 2026.962       2471.123  686.9221 4255.324
## 2026.981       2471.123  686.9221 4255.324

Tabla de pronósticos

semanas_pron <- paste0(
  "2026-S",
  sprintf("%02d", 1:h_pron)
)



tabla_pron <- data.frame(

  Semana = semanas_pron,

  Pronostico = round(
    as.numeric(pronosticos$mean),
    2
  ),

  LI_94 = round(
    as.numeric(pronosticos$lower),
    2
  ),

  LS_94 = round(
    as.numeric(pronosticos$upper),
    2
  )
)



knitr::kable(

  tabla_pron,

  caption = "Pronosticos de Ventas Semanales 2026 IC 94%",

  col.names = c(
    "Semana",
    "Pronostico",
    "L Inferior 94%",
    "L Superior 94%"
  ),

  align = "cccc"
)

Pronosticos de Ventas Semanales 2026 IC 94%
Semana	Pronostico	L Inferior 94%	L Superior 94%
2026-S01	2709.37	1114.03	4304.71
2026-S02	2577.80	829.83	4325.77
2026-S03	2518.89	741.89	4295.89
2026-S04	2492.51	709.75	4275.27
2026-S05	2480.70	696.79	4264.61
2026-S06	2475.41	691.27	4259.55
2026-S07	2473.04	688.85	4257.23
2026-S08	2471.98	687.78	4256.18
2026-S09	2471.51	687.31	4255.71
2026-S10	2471.30	687.09	4255.50
2026-S11	2471.20	687.00	4255.40
2026-S12	2471.16	686.96	4255.36
2026-S13	2471.14	686.94	4255.34
2026-S14	2471.13	686.93	4255.33
2026-S15	2471.13	686.93	4255.33
2026-S16	2471.12	686.92	4255.33
2026-S17	2471.12	686.92	4255.32
2026-S18	2471.12	686.92	4255.32
2026-S19	2471.12	686.92	4255.32
2026-S20	2471.12	686.92	4255.32
2026-S21	2471.12	686.92	4255.32
2026-S22	2471.12	686.92	4255.32
2026-S23	2471.12	686.92	4255.32
2026-S24	2471.12	686.92	4255.32
2026-S25	2471.12	686.92	4255.32
2026-S26	2471.12	686.92	4255.32
2026-S27	2471.12	686.92	4255.32
2026-S28	2471.12	686.92	4255.32
2026-S29	2471.12	686.92	4255.32
2026-S30	2471.12	686.92	4255.32
2026-S31	2471.12	686.92	4255.32
2026-S32	2471.12	686.92	4255.32
2026-S33	2471.12	686.92	4255.32
2026-S34	2471.12	686.92	4255.32
2026-S35	2471.12	686.92	4255.32
2026-S36	2471.12	686.92	4255.32
2026-S37	2471.12	686.92	4255.32
2026-S38	2471.12	686.92	4255.32
2026-S39	2471.12	686.92	4255.32
2026-S40	2471.12	686.92	4255.32
2026-S41	2471.12	686.92	4255.32
2026-S42	2471.12	686.92	4255.32
2026-S43	2471.12	686.92	4255.32
2026-S44	2471.12	686.92	4255.32
2026-S45	2471.12	686.92	4255.32
2026-S46	2471.12	686.92	4255.32
2026-S47	2471.12	686.92	4255.32
2026-S48	2471.12	686.92	4255.32
2026-S49	2471.12	686.92	4255.32
2026-S50	2471.12	686.92	4255.32
2026-S51	2471.12	686.92	4255.32
2026-S52	2471.12	686.92	4255.32

Gráfica de pronósticos con bandas de confianza

autoplot(pronosticos) +
  autolayer(ventas_ts52, series = "Observado",
            color = "steelblue", linewidth = 0.8) +
  labs(
    title    = "Pronóstico de Ventas Semanales — Metodología Box-Jenkins",
    subtitle = sprintf("ARIMA(%d,%d,%d)(%d,%d,%d)[%d] | IC al 94%% | Horizonte: 52 semanas (2026)",
                       p, d, q, P, D, Q, s),
    x        = "Semana",
    y        = "Ventas",
    color    = "Serie"
  ) +
  scale_color_manual(values = c("Observado" = "steelblue")) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")

## Ignoring unknown labels:
## • colour : "Serie"

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Interpretación de los pronósticos (Inciso iii):

Los pronósticos puntuales (línea azul) representan la mejor estimación futura de las ventas semanales basada en el modelo ARIMA ajustado.
La banda de confianza al 94% (región sombreada) indica el rango dentro del cual se espera caiga el valor real con un 94% de probabilidad. A medida que el horizonte de pronóstico se aleja, la incertidumbre crece y la banda se amplía.
Si el modelo seleccionado es un proceso estacionario (d=0), los pronósticos convergen hacia la media de largo plazo. Si tiene diferenciación (d≥1), los pronósticos siguen la tendencia implícita.
La amplitud de la banda al 94% (más estrecha que al 95%) refleja que se acepta un 6% de probabilidad de que el valor real quede fuera del intervalo.

Resumen de pronósticos clave

pron_med  <- as.numeric(pronosticos$mean)
pron_li   <- as.numeric(pronosticos$lower)
pron_ls   <- as.numeric(pronosticos$upper)

cat("╔══════════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════════╗

cat("   RESUMEN DE PRONÓSTICOS — IC 94%                            \n")

##    RESUMEN DE PRONÓSTICOS — IC 94%

cat("╠══════════════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════════════╣

cat(sprintf("  Pronóstico semana 1 (2026-S01): %.2f  [%.2f, %.2f]\n",
            pron_med[1],  pron_li[1],  pron_ls[1]))

##   Pronóstico semana 1 (2026-S01): 2709.37  [1114.03, 4304.71]

cat(sprintf("  Pronóstico semana 13 (trimestre): %.2f  [%.2f, %.2f]\n",
            pron_med[13], pron_li[13], pron_ls[13]))

##   Pronóstico semana 13 (trimestre): 2471.14  [686.94, 4255.34]

cat(sprintf("  Pronóstico semana 26 (semestre): %.2f  [%.2f, %.2f]\n",
            pron_med[26], pron_li[26], pron_ls[26]))

##   Pronóstico semana 26 (semestre): 2471.12  [686.92, 4255.32]

cat(sprintf("  Pronóstico semana 52 (anual):    %.2f  [%.2f, %.2f]\n",
            pron_med[52], pron_li[52], pron_ls[52]))

##   Pronóstico semana 52 (anual):    2471.12  [686.92, 4255.32]

cat(sprintf("\n  Media observada 2025: %.2f\n", mean(ventas)))

## 
##   Media observada 2025: 2460.05

cat(sprintf("  Amplitud promedio IC 94%% (s52): %.2f\n",
            mean(pron_ls - pron_li)))

##   Amplitud promedio IC 94% (s52): 3559.40

cat("\n╚══════════════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════════════╝

La tabla siguiente presenta los ingresos trimestrales de Southwest Airlines, antes de partidas extraordinarias (millones de $) para los años de 2008 a 2019.

Grafique los datos de ingresos como una serie de tiempo y describa cualquier patrón (característica o componente) existente.
¿La serie es estacionaria o no estacionaria? Explique su respuesta. Aplicando el test correspondiente
Calcule las autocorrelaciones de la serie de ingresos para los primeros 10 retrasos de tiempo. ¿El comportamiento de las autocorrelaciones es congruente con su elección realizada en el inciso b)? Explique su respuesta el test correspondiente.
Use R o Excel para calcular las diferencias cuartas de los datos de ingresos en la tabla anterior. Las diferencias cuartas se calculan diferenciando las observaciones separadas por cuatro periodos. Con datos trimestrales, este procedimiento algunas veces es útil en la creación de una serie estacionaria ,a partir de una serie no estacionaria. Por lo tanto, los datos diferenciados cuartos serán Y5 – Y1 = 19.64 – .17 = 19.47, Y6 – Y2 = 19.24 – 15.13 = 4.11, …, y así sucesivamente.
Grafique las series de tiempo de las diferencias cuartas. ¿Estas series de tiempo parecen ser estacionarias o no estacionarias? Explique su respuesta aplicando el test correspondiente

Datos: Ingresos Trimestrales Southwest Airlines (2008–2019)

library(forecast)
library(tseries)
library(ggplot2)
library(gridExtra)

# ── Datos (millones de $) ──────────────────────────────────────────────────
anios <- 2008:2019

ingresos_mat <- matrix(c(
  3.17,  15.13,  15.07,   9.17,
 19.64,  19.24,  24.57,   8.21,
  9.09,  23.53,  23.04,  -4.58,
 18.21,  10.52,  15.72,   8.84,
 12.48,  23.48,  26.89,  17.17,
 24.93,  42.15,  48.83,  38.37,
 41.85,  58.52,  68.62,  20.34,
 11.81,  59.72,  67.72,  43.46,
 33.00,  85.32,  60.86,  28.10,
 50.87,  93.83,  92.51,  80.55,
 70.01, 133.39, 129.64, 100.38,
 95.85, 157.76, 126.98,  93.80
), nrow = 12, byrow = TRUE,
dimnames = list(anios, c("TRIM1","TRIM2","TRIM3","TRIM4")))

knitr::kable(
  ingresos_mat,
  caption = "Ingresos Trimestrales de Southwest Airlines 2008–2019 (millones de $)",
  digits  = 2,
  align   = "cccc"
)

Ingresos Trimestrales de Southwest Airlines 2008–2019 (millones de $)
	TRIM1	TRIM2	TRIM3	TRIM4
2008	3.17	15.13	15.07	9.17
2009	19.64	19.24	24.57	8.21
2010	9.09	23.53	23.04	-4.58
2011	18.21	10.52	15.72	8.84
2012	12.48	23.48	26.89	17.17
2013	24.93	42.15	48.83	38.37
2014	41.85	58.52	68.62	20.34
2015	11.81	59.72	67.72	43.46
2016	33.00	85.32	60.86	28.10
2017	50.87	93.83	92.51	80.55
2018	70.01	133.39	129.64	100.38
2019	95.85	157.76	126.98	93.80

# Vector de datos ordenado cronológicamente (Q1 2008 → Q4 2019)
ingresos <- as.vector(t(ingresos_mat))
n        <- length(ingresos)
cat(sprintf("\nTotal de observaciones: %d trimestres\n", n))

## 
## Total de observaciones: 48 trimestres

Inciso a) — Gráfica de la serie de tiempo e identificación de patrones

# Serie trimestral: inicio 2008 Q1, frecuencia 4
ingresos_ts <- ts(ingresos, start = c(2008, 1), frequency = 4)

autoplot(ingresos_ts) +
  geom_point(color = "steelblue", size = 2.2) +
  geom_smooth(method = "loess", se = TRUE,
              color = "firebrick", fill = "salmon",
              alpha = 0.2, linewidth = 0.9) +
  labs(
    title    = "Ingresos Trimestrales de Southwest Airlines (2008–2019)",
    subtitle = "Millones de $ | Linea roja: tendencia LOESS con IC",
    x        = "Anio",
    y        = "Ingresos (millones de $)"
  ) +
  theme_minimal(base_size = 13)

## `geom_smooth()` using formula = 'y ~ x'

# Descomposición clásica para identificar componentes
decomp <- decompose(ingresos_ts, type = "additive")
autoplot(decomp) +
  labs(title = "Descomposición Clásica (Aditiva) — Ingresos Southwest Airlines") +
  theme_minimal(base_size = 12)

Interpretación e identificación de patrones (Inciso a):

|————|————-| | Tendencia | Claramente creciente a lo largo de 2008–2019. Los ingresos pasan de niveles bajos (incluso negativos en Q4-2010: −4.58) a valores superiores a 150 millones en 2019. Esto indica una serie NO estacionaria en media. | | Estacionalidad | Patrón repetitivo trimestral muy marcado: Q2 y Q3 (verano/primavera) registran ingresos notablemente más altos que Q1 y Q4. Esto refleja la mayor demanda de vuelos en temporadas vacacionales. | | Componente irregular | Fluctuaciones aleatorias alrededor de la tendencia + estacionalidad, más pronunciadas hacia el final del período. | | Tipo de serie | No estacionaria: presenta tendencia creciente y componente estacional sistemático. |

Inciso b) — ¿La serie es estacionaria o no estacionaria?

Prueba ADF (Dickey-Fuller Aumentada)

\[H_0: \text{Raíz unitaria (NO estacionaria)} \quad H_1: \text{Estacionaria}\]

adf_orig <- adf.test(ingresos_ts, alternative = "stationary")
cat("=== ADF — Serie Original ===\n")

## === ADF — Serie Original ===

cat(sprintf("  Estadistico:  %.4f\n", adf_orig$statistic))

##   Estadistico:  -0.9352

cat(sprintf("  p-valor:      %.4f\n", adf_orig$p.value))

##   p-valor:      0.9386

cat(ifelse(adf_orig$p.value < 0.05,
           "  → p < 0.05 → Serie ES estacionaria ✓\n",
           "  → p ≥ 0.05 → Serie NO es estacionaria ✗\n"))

##   → p ≥ 0.05 → Serie NO es estacionaria ✗

Prueba KPSS

\[H_0: \text{Serie estacionaria} \quad H_1: \text{No estacionaria}\]

kpss_orig <- kpss.test(ingresos_ts, null = "Level")

## Warning in kpss.test(ingresos_ts, null = "Level"): p-value smaller than printed
## p-value

cat("=== KPSS — Serie Original ===\n")

## === KPSS — Serie Original ===

cat(sprintf("  Estadistico:  %.4f\n", kpss_orig$statistic))

##   Estadistico:  1.1208

cat(sprintf("  p-valor:      %.4f\n", kpss_orig$p.value))

##   p-valor:      0.0100

cat(ifelse(kpss_orig$p.value < 0.05,
           "  → p < 0.05 → Serie NO es estacionaria ✗\n",
           "  → p ≥ 0.05 → Serie ES estacionaria ✓\n"))

##   → p < 0.05 → Serie NO es estacionaria ✗

Prueba Phillips-Perron

pp_orig <- pp.test(ingresos_ts, alternative = "stationary")
cat("=== Phillips-Perron — Serie Original ===\n")

## === Phillips-Perron — Serie Original ===

cat(sprintf("  Estadístico:  %.4f\n", pp_orig$statistic))

##   Estadístico:  -25.1879

cat(sprintf("  p-valor:      %.4f\n", pp_orig$p.value))

##   p-valor:      0.0106

cat(ifelse(pp_orig$p.value < 0.05,
           "  → p < 0.05 → Serie ES estacionaria ✓\n",
           "  → p ≥ 0.05 → Serie NO es estacionaria ✗\n"))

##   → p < 0.05 → Serie ES estacionaria ✓

Conclusión Inciso b): La serie de ingresos de Southwest Airlines NO es estacionaria. La tendencia creciente evidente en la gráfica, el patrón estacional marcado, y los resultados de las pruebas formales (especialmente KPSS que rechaza estacionariedad) confirman que se requiere diferenciación —regular y/o estacional— para modelar la serie con ARIMA.

Inciso c) — Autocorrelaciones para los primeros 10 retrasos

Cálculo de autocorrelaciones

\[r_k = \frac{\sum_{t=k+1}^{n}(Y_t - \bar{Y})(Y_{t-k} - \bar{Y})}{\sum_{t=1}^{n}(Y_t - \bar{Y})^2}, \quad k = 1, \ldots, 10\]

rk_manual <- function(serie, k) {

  n <- length(serie)

  Ybar <- mean(serie)

  num <- sum(
    (serie[(k+1):n] - Ybar) *
    (serie[1:(n-k)] - Ybar)
  )

  den <- sum(
    (serie - Ybar)^2
  )

  num / den
}



lags_10 <- 1:10



acf_orig <- sapply(
  lags_10,
  function(k) rk_manual(ingresos, k)
)



ic_95 <- qnorm(0.975) / sqrt(n)



tabla_acf <- data.frame(

  k = lags_10,

  r_k = round(
    acf_orig,
    4
  ),

  IC_95 = round(
    ic_95,
    4
  ),

  Resultado = ifelse(
    abs(acf_orig) > ic_95,
    "Significativo",
    "No significativo"
  )
)



knitr::kable(

  tabla_acf,

  caption = "Autocorrelaciones Serie Original 10 retrasos",

  col.names = c(
    "Retraso k",
    "r_k",
    "IC95",
    "Significativo"
  ),

  align = "cccc"
)

Autocorrelaciones Serie Original 10 retrasos
Retraso k	r_k	IC95	Significativo
1	0.7925	0.2829	Significativo
2	0.6156	0.2829	Significativo
3	0.6464	0.2829	Significativo
4	0.6971	0.2829	Significativo
5	0.5150	0.2829	Significativo
6	0.3538	0.2829	Significativo
7	0.3746	0.2829	Significativo
8	0.4245	0.2829	Significativo
9	0.2679	0.2829	No significativo
10	0.1318	0.2829	No significativo

Gráfica ACF — Serie original

acf_df <- data.frame(Lag = lags_10, ACF = acf_orig)

ggplot(acf_df, aes(x = Lag, y = ACF)) +
  geom_hline(yintercept = 0, color = "gray50") +
  geom_hline(yintercept =  ic_95, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  geom_hline(yintercept = -ic_95, linetype = "dashed", color = "firebrick", linewidth = 0.8) +
  geom_segment(aes(xend = Lag, yend = 0), color = "steelblue", linewidth = 1.4) +
  geom_point(color = "steelblue", size = 3.5) +
  geom_text(aes(label = round(ACF, 3)),
            vjust = ifelse(acf_orig >= 0, -0.8, 1.6), size = 3.3) +
  annotate("rect", xmin = 0.4, xmax = 10.6, ymin = -ic_95, ymax = ic_95,
           fill = "steelblue", alpha = 0.07) +
  scale_x_continuous(breaks = 1:10) +
  labs(
    title    = "ACF — Ingresos Trimestrales Southwest Airlines (Serie Original)",
    subtitle = "Primeros 10 retrasos | Líneas rojas: IC 95% (±2/√n)",
    x        = "Retraso (k)", y = expression(r[k])
  ) +
  theme_minimal(base_size = 13)

Interpretación ACF (Inciso c): Los coeficientes de autocorrelación decaen muy lentamente y son significativamente distintos de cero para casi todos los retrasos (superan ampliamente las bandas de confianza al 95%). Este comportamiento es el sello característico de una serie no estacionaria con tendencia: la dependencia entre observaciones persiste a lo largo del tiempo. Adicionalmente, se puede observar un patrón ondulatorio con picos cada 4 rezagos (lag 4, lag 8), lo cual es consistente con la estacionalidad trimestral (s=4). Este resultado es completamente congruente con la conclusión del inciso b): la serie NO es estacionaria.

Inciso d) — Diferencias cuartas

Definición y cálculo

Las diferencias cuartas separan observaciones por cuatro períodos (un año completo para datos trimestrales):

\[\nabla_4 Y_t = Y_t - Y_{t-4}\]

Por ejemplo: $Y_5 - Y_1 = 19.64 - 3.17 = 16.47$, $Y_6 - Y_2 = 19.24 - 15.13 = 4.11$, etc.

diff4        <- diff(ingresos, lag = 4)
diff4_ts     <- diff(ingresos_ts, lag = 4)
n_diff4      <- length(diff4)

# Fechas para las diferencias cuartas (desde Q1 2009)
fechas_orig  <- time(ingresos_ts)
fechas_diff4 <- time(diff4_ts)

# Tabla parcial
df_diff4 <- data.frame(
  t      = 5:(5 + n_diff4 - 1),
  Y_t    = ingresos[5:n],
  Y_t4   = ingresos[1:(n-4)],
  D4Y_t  = round(diff4, 4)
)

knitr::kable(
  head(df_diff4, 20),
  caption   = "Diferencias Cuartas ∇₄Yₜ = Yₜ − Yₜ₋₄ (primeras 20 observaciones)",
  col.names = c("t", "Y_t", "Y_{t-4}", "∇₄Y_t"),
  align     = "cccc"
)

Diferencias Cuartas ∇₄Yₜ = Yₜ − Yₜ₋₄ (primeras 20 observaciones)
t	Y_t	Y_{t-4}	∇₄Y_t
5	19.64	3.17	16.47
6	19.24	15.13	4.11
7	24.57	15.07	9.50
8	8.21	9.17	-0.96
9	9.09	19.64	-10.55
10	23.53	19.24	4.29
11	23.04	24.57	-1.53
12	-4.58	8.21	-12.79
13	18.21	9.09	9.12
14	10.52	23.53	-13.01
15	15.72	23.04	-7.32
16	8.84	-4.58	13.42
17	12.48	18.21	-5.73
18	23.48	10.52	12.96
19	26.89	15.72	11.17
20	17.17	8.84	8.33
21	24.93	12.48	12.45
22	42.15	23.48	18.67
23	48.83	26.89	21.94
24	38.37	17.17	21.20

cat(sprintf("\nTotal de diferencias cuartas calculadas: %d\n", n_diff4))

## 
## Total de diferencias cuartas calculadas: 44

cat(sprintf("Media de ∇₄Y_t: %.4f\n", mean(diff4)))

## Media de ∇₄Y_t: 9.8148

cat(sprintf("Desv. estándar: %.4f\n", sd(diff4)))

## Desv. estándar: 16.7117

Interpretación Inciso d): Las diferencias cuartas $\nabla_4 Y_t = Y_t - Y_{t-4}$ eliminan el componente estacional de período 4 (trimestral). Al restar cada observación del mismo trimestre del año anterior, se remueve el patrón estacional sistemático. Los valores resultantes representan el cambio anual en los ingresos para cada trimestre, que debería fluctuar con menor estructura sistemática que la serie original.

Inciso e) — Gráfica y test de estacionariedad de diferencias cuartas

Gráfica de las diferencias cuartas

autoplot(diff4_ts) +
  geom_point(color = "darkorange", size = 2.2) +
  geom_hline(yintercept = mean(diff4), linetype = "dashed",
             color = "firebrick", linewidth = 0.9) +
  geom_smooth(method = "loess", se = FALSE,
              color = "darkgreen", linewidth = 0.8, linetype = "dotdash") +
  labs(
    title    = "Diferencias Cuartas (∇₄Yₜ) — Ingresos Southwest Airlines",
    subtitle = paste0("Linea roja: media = ", round(mean(diff4), 2),
                      " | Verde: tendencia LOESS"),
    x        = "Anio", y = "∇₄Y_t (millones de $)"
  ) +
  theme_minimal(base_size = 13)

## `geom_smooth()` using formula = 'y ~ x'

Interpretación gráfica (Inciso e): Tras aplicar la diferenciación estacional de orden 4, la serie resultante oscila con mayor regularidad alrededor de su media. Sin embargo, se aprecia que los valores son predominantemente positivos y la tendencia LOESS no es perfectamente horizontal, lo que sugiere que podría persistir alguna tendencia no estacional. La diferenciación estacional ha reducido la estructura periódica, pero puede ser necesaria una diferenciación regular adicional para alcanzar la plena estacionariedad.

ACF de las diferencias cuartas

acf_d4 <- sapply(
  1:10,
  function(k) rk_manual(diff4, k)
)



ic_d4 <- qnorm(0.975) / sqrt(n_diff4)



acf_d4_df <- data.frame(

  Lag = 1:10,

  ACF = acf_d4
)



ggplot(acf_d4_df, aes(x = Lag, y = ACF)) +

  geom_hline(
    yintercept = 0,
    color = "gray50"
  ) +

  geom_hline(
    yintercept = ic_d4,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.8
  ) +

  geom_hline(
    yintercept = -ic_d4,
    linetype = "dashed",
    color = "firebrick",
    linewidth = 0.8
  ) +

  geom_segment(
    aes(xend = Lag, yend = 0),
    color = "darkorange",
    linewidth = 1.4
  ) +

  geom_point(
    color = "darkorange",
    size = 3.5
  ) +

  geom_text(
    aes(label = round(ACF, 3)),
    vjust = ifelse(acf_d4 >= 0, -0.8, 1.6),
    size = 3.3
  ) +

  annotate(
    "rect",
    xmin = 0.4,
    xmax = 10.6,
    ymin = -ic_d4,
    ymax = ic_d4,
    fill = "darkorange",
    alpha = 0.07
  ) +

  scale_x_continuous(
    breaks = 1:10
  ) +

  labs(

    title = "ACF Diferencias Cuartas",

    subtitle = "Primeros 10 retrasos - Lineas rojas IC 95%",

    x = "Retraso k",

    y = expression(r[k])
  ) +

  theme_minimal(base_size = 13)

  theme_minimal(base_size = 13)

## <theme> List of 144
##  $ line                            : <ggplot2::element_line>
##   ..@ colour       : chr "black"
##   ..@ linewidth    : num 0.591
##   ..@ linetype     : num 1
##   ..@ lineend      : chr "butt"
##   ..@ linejoin     : chr "round"
##   ..@ arrow        : logi FALSE
##   ..@ arrow.fill   : chr "black"
##   ..@ inherit.blank: logi TRUE
##  $ rect                            : <ggplot2::element_rect>
##   ..@ fill         : chr "white"
##   ..@ colour       : chr "black"
##   ..@ linewidth    : num 0.591
##   ..@ linetype     : num 1
##   ..@ linejoin     : chr "round"
##   ..@ inherit.blank: logi TRUE
##  $ text                            : <ggplot2::element_text>
##   ..@ family       : chr ""
##   ..@ face         : chr "plain"
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : chr "black"
##   ..@ size         : num 13
##   ..@ hjust        : num 0.5
##   ..@ vjust        : num 0.5
##   ..@ angle        : num 0
##   ..@ lineheight   : num 0.9
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 0 0 0
##   ..@ debug        : logi FALSE
##   ..@ inherit.blank: logi TRUE
##  $ title                           : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : NULL
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ point                           : <ggplot2::element_point>
##   ..@ colour       : chr "black"
##   ..@ shape        : num 19
##   ..@ size         : num 1.77
##   ..@ fill         : chr "white"
##   ..@ stroke       : num 0.591
##   ..@ inherit.blank: logi TRUE
##  $ polygon                         : <ggplot2::element_polygon>
##   ..@ fill         : chr "white"
##   ..@ colour       : chr "black"
##   ..@ linewidth    : num 0.591
##   ..@ linetype     : num 1
##   ..@ linejoin     : chr "round"
##   ..@ inherit.blank: logi TRUE
##  $ geom                            : <ggplot2::element_geom>
##   ..@ ink        : chr "black"
##   ..@ paper      : chr "white"
##   ..@ accent     : chr "#3366FF"
##   ..@ linewidth  : num 0.591
##   ..@ borderwidth: num 0.591
##   ..@ linetype   : int 1
##   ..@ bordertype : int 1
##   ..@ family     : chr ""
##   ..@ fontsize   : num 4.57
##   ..@ pointsize  : num 1.77
##   ..@ pointshape : num 19
##   ..@ colour     : NULL
##   ..@ fill       : NULL
##  $ spacing                         : 'simpleUnit' num 6.5points
##   ..- attr(*, "unit")= int 8
##  $ margins                         : <ggplot2::margin> num [1:4] 6.5 6.5 6.5 6.5
##  $ aspect.ratio                    : NULL
##  $ axis.title                      : NULL
##  $ axis.title.x                    : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : num 1
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 3.25 0 0 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.title.x.top                : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : num 0
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 0 3.25 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.title.x.bottom             : NULL
##  $ axis.title.y                    : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : num 1
##   ..@ angle        : num 90
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 3.25 0 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.title.y.left               : NULL
##  $ axis.title.y.right              : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : num 1
##   ..@ angle        : num -90
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 0 0 3.25
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text                       : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : chr "#4D4D4DFF"
##   ..@ size         : 'rel' num 0.8
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : NULL
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.x                     : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : num 1
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 2.6 0 0 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.x.top                 : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 0 5.85 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.x.bottom              : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 5.85 0 0 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.y                     : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : num 1
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 2.6 0 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.y.left                : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 5.85 0 0
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.y.right               : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 0 0 5.85
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.text.theta                 : NULL
##  $ axis.text.r                     : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : num 0.5
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : <ggplot2::margin> num [1:4] 0 2.6 0 2.6
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ axis.ticks                      : <ggplot2::element_blank>
##  $ axis.ticks.x                    : NULL
##  $ axis.ticks.x.top                : NULL
##  $ axis.ticks.x.bottom             : NULL
##  $ axis.ticks.y                    : NULL
##  $ axis.ticks.y.left               : NULL
##  $ axis.ticks.y.right              : NULL
##  $ axis.ticks.theta                : NULL
##  $ axis.ticks.r                    : NULL
##  $ axis.minor.ticks.x.top          : NULL
##  $ axis.minor.ticks.x.bottom       : NULL
##  $ axis.minor.ticks.y.left         : NULL
##  $ axis.minor.ticks.y.right        : NULL
##  $ axis.minor.ticks.theta          : NULL
##  $ axis.minor.ticks.r              : NULL
##  $ axis.ticks.length               : 'rel' num 0.5
##  $ axis.ticks.length.x             : NULL
##  $ axis.ticks.length.x.top         : NULL
##  $ axis.ticks.length.x.bottom      : NULL
##  $ axis.ticks.length.y             : NULL
##  $ axis.ticks.length.y.left        : NULL
##  $ axis.ticks.length.y.right       : NULL
##  $ axis.ticks.length.theta         : NULL
##  $ axis.ticks.length.r             : NULL
##  $ axis.minor.ticks.length         : 'rel' num 0.75
##  $ axis.minor.ticks.length.x       : NULL
##  $ axis.minor.ticks.length.x.top   : NULL
##  $ axis.minor.ticks.length.x.bottom: NULL
##  $ axis.minor.ticks.length.y       : NULL
##  $ axis.minor.ticks.length.y.left  : NULL
##  $ axis.minor.ticks.length.y.right : NULL
##  $ axis.minor.ticks.length.theta   : NULL
##  $ axis.minor.ticks.length.r       : NULL
##  $ axis.line                       : <ggplot2::element_blank>
##  $ axis.line.x                     : NULL
##  $ axis.line.x.top                 : NULL
##  $ axis.line.x.bottom              : NULL
##  $ axis.line.y                     : NULL
##  $ axis.line.y.left                : NULL
##  $ axis.line.y.right               : NULL
##  $ axis.line.theta                 : NULL
##  $ axis.line.r                     : NULL
##  $ legend.background               : <ggplot2::element_blank>
##  $ legend.margin                   : NULL
##  $ legend.spacing                  : 'rel' num 2
##  $ legend.spacing.x                : NULL
##  $ legend.spacing.y                : NULL
##  $ legend.key                      : <ggplot2::element_blank>
##  $ legend.key.size                 : 'simpleUnit' num 1.2lines
##   ..- attr(*, "unit")= int 3
##  $ legend.key.height               : NULL
##  $ legend.key.width                : NULL
##  $ legend.key.spacing              : NULL
##  $ legend.key.spacing.x            : NULL
##  $ legend.key.spacing.y            : NULL
##  $ legend.key.justification        : NULL
##  $ legend.frame                    : NULL
##  $ legend.ticks                    : NULL
##  $ legend.ticks.length             : 'rel' num 0.2
##  $ legend.axis.line                : NULL
##  $ legend.text                     : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : 'rel' num 0.8
##   ..@ hjust        : NULL
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : NULL
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ legend.text.position            : NULL
##  $ legend.title                    : <ggplot2::element_text>
##   ..@ family       : NULL
##   ..@ face         : NULL
##   ..@ italic       : chr NA
##   ..@ fontweight   : num NA
##   ..@ fontwidth    : num NA
##   ..@ colour       : NULL
##   ..@ size         : NULL
##   ..@ hjust        : num 0
##   ..@ vjust        : NULL
##   ..@ angle        : NULL
##   ..@ lineheight   : NULL
##   ..@ margin       : NULL
##   ..@ debug        : NULL
##   ..@ inherit.blank: logi TRUE
##  $ legend.title.position           : NULL
##  $ legend.position                 : chr "right"
##  $ legend.position.inside          : NULL
##  $ legend.direction                : NULL
##  $ legend.byrow                    : NULL
##  $ legend.justification            : chr "center"
##  $ legend.justification.top        : NULL
##  $ legend.justification.bottom     : NULL
##  $ legend.justification.left       : NULL
##  $ legend.justification.right      : NULL
##  $ legend.justification.inside     : NULL
##   [list output truncated]
##  @ complete: logi TRUE
##  @ validate: logi TRUE

Tests formales sobre diferencias cuartas

adf_d4  <- adf.test(diff4_ts,  alternative = "stationary")
kpss_d4 <- kpss.test(diff4_ts, null = "Level")

## Warning in kpss.test(diff4_ts, null = "Level"): p-value greater than printed
## p-value

pp_d4   <- pp.test(diff4_ts,   alternative = "stationary")

## Warning in pp.test(diff4_ts, alternative = "stationary"): p-value smaller than
## printed p-value

cat("╔══════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════╗

cat("   TESTS DE ESTACIONARIEDAD — DIFERENCIAS CUARTAS ∇₄Yₜ    \n")

##    TESTS DE ESTACIONARIEDAD — DIFERENCIAS CUARTAS ∇₄Yₜ

cat("╠══════════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════════╣

cat(sprintf("ADF:  estad. = %7.4f | p = %.4f | %s\n",
    adf_d4$statistic,  adf_d4$p.value,
    ifelse(adf_d4$p.value  < 0.05, "Estacionaria ✓", "No estacionaria ✗")))

## ADF:  estad. = -3.9363 | p = 0.0212 | Estacionaria ✓

cat(sprintf("KPSS: estad. = %7.4f | p = %.4f | %s\n",
    kpss_d4$statistic, kpss_d4$p.value,
    ifelse(kpss_d4$p.value >= 0.05, "Estacionaria ✓", "No estacionaria ✗")))

## KPSS: estad. =  0.3394 | p = 0.1000 | Estacionaria ✓

cat(sprintf("PP:   estad. = %7.4f | p = %.4f | %s\n",
    pp_d4$statistic,   pp_d4$p.value,
    ifelse(pp_d4$p.value   < 0.05, "Estacionaria ✓", "No estacionaria ✗")))

## PP:   estad. = -27.1976 | p = 0.0100 | Estacionaria ✓

cat("\n╚══════════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════════╝

Conclusión Inciso e): Tras la diferenciación estacional de orden 4, los tests formales determinan si se ha alcanzado la estacionariedad. Si aún persiste no estacionariedad (p.ej. KPSS rechaza $H_0$), se requeriría una diferenciación regular adicional $\nabla_1\nabla_4 Y_t = (1-B)(1-B^4)Y_t$, lo cual es la base del modelo SARIMA con $d=1, D=1$.

Inciso i) — Configuración `ts()`, año final 2025, AST completo

Serie configurada con año final 2025

Los datos son de 2008 a 2019 (12 años × 4 trimestres = 48 obs). Para que el año final sea 2025, el inicio de la serie es $2025 - 12 + 1/4 \approx$ inicio en Q1 de 2014.

Nota interpretativa: Si se mantienen los datos tal como están y se proyecta hasta 2025, el inicio de la serie observada se ajusta proporcionalmente. Aquí configuramos la serie con el año final = 2025, por lo que el inicio corresponde a $2025 - (48/4) + 1 = 2014$.

# Configuración: 48 trimestres, año final 2025 → inicio 2014 Q1
ingresos_ts25 <- ts(ingresos, start = c(2014, 1), frequency = 4)

cat("Serie configurada con año final 2025:\n")

## Serie configurada con año final 2025:

print(ingresos_ts25)

##        Qtr1   Qtr2   Qtr3   Qtr4
## 2014   3.17  15.13  15.07   9.17
## 2015  19.64  19.24  24.57   8.21
## 2016   9.09  23.53  23.04  -4.58
## 2017  18.21  10.52  15.72   8.84
## 2018  12.48  23.48  26.89  17.17
## 2019  24.93  42.15  48.83  38.37
## 2020  41.85  58.52  68.62  20.34
## 2021  11.81  59.72  67.72  43.46
## 2022  33.00  85.32  60.86  28.10
## 2023  50.87  93.83  92.51  80.55
## 2024  70.01 133.39 129.64 100.38
## 2025  95.85 157.76 126.98  93.80

cat(sprintf("\nInicio: %d-Q%d | Fin: %d-Q%d | Frecuencia: %d (trimestral)\n",
            start(ingresos_ts25)[1], start(ingresos_ts25)[2],
            end(ingresos_ts25)[1],   end(ingresos_ts25)[2],
            frequency(ingresos_ts25)))

## 
## Inicio: 2014-Q1 | Fin: 2025-Q4 | Frecuencia: 4 (trimestral)

Descomposición STL

stl_mod <- stl(ingresos_ts25, s.window = "periodic", robust = TRUE)
autoplot(stl_mod) +
  labs(title = "Descomposición STL — Ingresos Southwest Airlines",
       subtitle = "Serie configurada con anio final 2025") +
  theme_minimal(base_size = 12)

Interpretación descomposición STL:

Tendencia: Crecimiento sostenido a lo largo de todo el período. La tendencia STL captura el incremento progresivo en los ingresos, con una ligera desaceleración en los años iniciales de la muestra.
Estacionalidad: Patrón regular y repetitivo con período s=4 (trimestral). Los trimestres 2 y 3 (Q2, Q3) muestran consistentemente valores más altos (mayor demanda estival), mientras Q1 y Q4 son menores. La amplitud estacional se mantiene relativamente constante (esquema aditivo).
Residuo: Fluctuaciones irregulares de menor magnitud que la tendencia y la estacionalidad. Si los residuos son pequeños respecto a los otros componentes, el modelo STL captura bien la estructura de la serie.

Tests de estacionariedad (serie 2025)

adf_i  <- adf.test(ingresos_ts25,  alternative = "stationary")
kpss_i <- kpss.test(ingresos_ts25, null = "Level")

## Warning in kpss.test(ingresos_ts25, null = "Level"): p-value smaller than
## printed p-value

pp_i   <- pp.test(ingresos_ts25,   alternative = "stationary")

cat("╔══════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════╗

cat("   TESTS — SERIE TRIMESTRAL (año final 2025)           \n")

##    TESTS — SERIE TRIMESTRAL (año final 2025)

cat("╠══════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════╣

cat(sprintf("ADF:  p = %.4f → %s\n", adf_i$p.value,
    ifelse(adf_i$p.value  < 0.05, "Estacionaria ✓", "NO estacionaria ✗")))

## ADF:  p = 0.9386 → NO estacionaria ✗

cat(sprintf("KPSS: p = %.4f → %s\n", kpss_i$p.value,
    ifelse(kpss_i$p.value >= 0.05, "Estacionaria ✓", "NO estacionaria ✗")))

## KPSS: p = 0.0100 → NO estacionaria ✗

cat(sprintf("PP:   p = %.4f → %s\n", pp_i$p.value,
    ifelse(pp_i$p.value   < 0.05, "Estacionaria ✓", "NO estacionaria ✗")))

## PP:   p = 0.0106 → Estacionaria ✓

cat("\n╚══════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════╝

ACF y PACF — Identificación del modelo

par(mfrow = c(3, 2), mar = c(4, 4, 3, 1))

acf(ingresos_ts25,              lag.max = 20, main = "ACF — Serie Original")
pacf(ingresos_ts25,             lag.max = 20, main = "PACF — Serie Original")
acf(diff(ingresos_ts25),        lag.max = 20, main = "ACF — 1ª Diferencia (d=1)")
pacf(diff(ingresos_ts25),       lag.max = 20, main = "PACF — 1ª Diferencia (d=1)")
acf(diff(diff(ingresos_ts25, lag=4)),  lag.max = 20,
    main = "ACF — Diferencia Regular + Estacional (d=1, D=1)")
pacf(diff(diff(ingresos_ts25, lag=4)), lag.max = 20,
    main = "PACF — Diferencia Regular + Estacional (d=1, D=1)")

par(mfrow = c(1, 1))

Interpretación ACF/PACF:

Serie original: ACF con decaimiento lento → no estacionaria. Picos cada lag=4 → estacionalidad trimestral.
1ª diferencia: Reduce la tendencia pero persisten picos en múltiplos de 4 → sigue estacionalidad.
Diferencia regular + estacional: Idealmente, el ACF y PACF deben mostrar cortes abruptos que guían los órdenes $p, q$ (no estacionales) y $P, Q$ (estacionales).

Identificación del modelo con `auto.arima()`

modelo_arima <- auto.arima(
  ingresos_ts25,
  seasonal      = TRUE,
  stepwise      = FALSE,
  approximation = FALSE,
  ic            = "aicc",
  trace         = FALSE
)

cat("╔══════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════╗

cat("   MODELO IDENTIFICADO POR auto.arima()                    \n")

##    MODELO IDENTIFICADO POR auto.arima()

cat("╚══════════════════════════════════════════════════════════╝\n\n")

## ╚══════════════════════════════════════════════════════════╝

summary(modelo_arima)

## Series: ingresos_ts25 
## ARIMA(0,1,1)(0,1,1)[4] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4318  -0.5081
## s.e.   0.1677   0.1211
## 
## sigma^2 = 195.5:  log likelihood = -174.14
## AIC=354.27   AICc=354.89   BIC=359.56
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE     MASE       ACF1
## Training set 0.1369126 12.92328 10.02165 -2.044176 41.41787 0.633107 0.04077387

p <- modelo_arima$arma[1]; d <- modelo_arima$arma[6]; q <- modelo_arima$arma[2]
P <- modelo_arima$arma[3]; D <- modelo_arima$arma[7]; Q <- modelo_arima$arma[4]
s <- frequency(ingresos_ts25)

cat(sprintf("\n→ Notación: ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n",
            p, d, q, P, D, Q, s))

## 
## → Notación: ARIMA(0,1,1)(0,1,1)[4]

Ecuación específica del modelo seleccionado

coef_mod <- coef(modelo_arima)

cat("╔══════════════════════════════════════════════════════════════╗\n")

## ╔══════════════════════════════════════════════════════════════╗

cat(sprintf("  Modelo: ARIMA(%d,%d,%d)(%d,%d,%d)[%d]\n", p,d,q,P,D,Q,s))

##   Modelo: ARIMA(0,1,1)(0,1,1)[4]

cat("╠══════════════════════════════════════════════════════════════╣\n\n")

## ╠══════════════════════════════════════════════════════════════╣

cat("  Coeficientes estimados:\n\n")

##   Coeficientes estimados:

for (nm in names(coef_mod)) {
  cat(sprintf("    %-15s = %9.4f\n", nm, coef_mod[nm]))
}

##     ma1             =   -0.4318
##     sma1            =   -0.5081

cat(sprintf("\n  σ²    = %.4f\n", modelo_arima$sigma2))

## 
##   σ²    = 195.5253

cat(sprintf("  AICc  = %.4f\n", modelo_arima$aicc))

##   AICc  = 354.8873

cat(sprintf("  AIC   = %.4f\n", modelo_arima$aic))

##   AIC   = 354.2719

cat(sprintf("  BIC   = %.4f\n", modelo_arima$bic))

##   BIC   = 359.5555

cat("\n╚══════════════════════════════════════════════════════════════╝\n")

## 
## ╚══════════════════════════════════════════════════════════════╝

Diagnóstico de residuos

checkresiduals(modelo_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[4]
## Q* = 5.0148, df = 6, p-value = 0.5419
## 
## Model df: 2.   Total lags used: 8

lb_arima <- Box.test(residuals(modelo_arima), lag = 8, type = "Ljung-Box")
cat(sprintf("\nLjung-Box (lag=8): Q = %.4f, p-valor = %.4f\n",
            lb_arima$statistic, lb_arima$p.value))

## 
## Ljung-Box (lag=8): Q = 5.0148, p-valor = 0.7560

cat(ifelse(lb_arima$p.value > 0.05,
           "→ Residuos son RUIDO BLANCO → Modelo adecuado ✓\n",
           "→ Residuos NO son ruido blanco → Considerar ajustes.\n"))

## → Residuos son RUIDO BLANCO → Modelo adecuado ✓

Interpretación diagnóstico: Un buen ajuste se confirma cuando: (1) el gráfico de residuos no muestra patrón sistemático, (2) el ACF de residuos está dentro de las bandas de confianza, (3) los residuos siguen aproximadamente una distribución normal, y (4) la prueba Ljung-Box no rechaza la hipótesis de ruido blanco ($p > 0.05$).

Inciso iii) — Pronósticos con IC al 94% — Metodología Box-Jenkins

Pronósticos para 2026–2028 (12 trimestres)

h_pron <- 12   # 3 años × 4 trimestres = 12 períodos

pronosticos <- forecast(modelo_arima,
                        h     = h_pron,
                        level = 94)

cat("=== Pronósticos ARIMA — 2026Q1 a 2028Q4 (IC 94%) ===\n")

## === Pronósticos ARIMA — 2026Q1 a 2028Q4 (IC 94%) ===

print(pronosticos)

##         Point Forecast     Lo 94    Hi 94
## 2026 Q1       98.96225  72.66304 125.2615
## 2026 Q2      156.13406 125.88553 186.3826
## 2026 Q3      138.82676 105.08808 172.5654
## 2026 Q4      109.76981  72.86963 146.6700
## 2027 Q1      110.92400  64.67543 157.1726
## 2027 Q2      168.09581 116.75393 219.4377
## 2027 Q3      150.78851  94.81488 206.7621
## 2027 Q4      121.73156  61.48120 181.9819
## 2028 Q1      122.88575  53.09086 192.6806
## 2028 Q2      180.05756 104.22758 255.8875
## 2028 Q3      162.75026  81.33131 244.1692
## 2028 Q4      133.69331  47.04515 220.3415

Tabla de pronósticos

anios_p   <- rep(2026:2028, each = 4)
trims_p   <- rep(paste0("Q", 1:4), 3)
labels_p  <- paste0(anios_p, "-", trims_p)

tabla_pron <- data.frame(
  Periodo    = labels_p,
  Pronostico = round(as.numeric(pronosticos$mean),  2),
  LI_94      = round(as.numeric(pronosticos$lower), 2),
  LS_94      = round(as.numeric(pronosticos$upper), 2),
  Amplitud   = round(
    as.numeric(pronosticos$upper) -
    as.numeric(pronosticos$lower), 2
  )
)

knitr::kable(
  tabla_pron,
  caption   = "Pronosticos de Ingresos Southwest Airlines (2026-2028) con IC al 94%",
  col.names = c(
    "Periodo",
    "Pronostico",
    "L. Inferior 94%",
    "L. Superior 94%",
    "Amplitud IC"
  ),
  align = "ccccc"
)

Pronosticos de Ingresos Southwest Airlines (2026-2028) con IC al 94%
Periodo	Pronostico	L. Inferior 94%	L. Superior 94%	Amplitud IC
2026-Q1	98.96	72.66	125.26	52.60
2026-Q2	156.13	125.89	186.38	60.50
2026-Q3	138.83	105.09	172.57	67.48
2026-Q4	109.77	72.87	146.67	73.80
2027-Q1	110.92	64.68	157.17	92.50
2027-Q2	168.10	116.75	219.44	102.68
2027-Q3	150.79	94.81	206.76	111.95
2027-Q4	121.73	61.48	181.98	120.50
2028-Q1	122.89	53.09	192.68	139.59
2028-Q2	180.06	104.23	255.89	151.66
2028-Q3	162.75	81.33	244.17	162.84
2028-Q4	133.69	47.05	220.34	173.30

Gráfica de pronósticos

autoplot(pronosticos) +
  autolayer(ingresos_ts25, series = "Observado",
            color = "steelblue", linewidth = 0.9) +
  labs(
    title    = "Pronósticos de Ingresos — Southwest Airlines (Metodología Box-Jenkins)",
    subtitle = sprintf("ARIMA(%d,%d,%d)(%d,%d,%d)[%d] | IC al 94%% | 2026–2028",
                       p, d, q, P, D, Q, s),
    x        = "Anio", y        = "Ingresos (millones de $)",
    color    = "Serie"
  ) +
  scale_color_manual(values = c("Observado" = "steelblue")) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "bottom")

## Ignoring unknown labels:
## • colour : "Serie"

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.

Interpretación de los pronósticos (Inciso iii):

Los pronósticos puntuales capturan tanto la tendencia creciente como el patrón estacional: Q2 y Q3 de cada año proyectado muestran valores más altos, en línea con la estacionalidad histórica de la aerolínea.
La banda de confianza al 94% (más estrecha que al 95%) se amplía progresivamente: la incertidumbre es menor para los trimestres cercanos (2026) y mayor para los lejanos (2028). Este ensanchamiento es una propiedad inherente de los modelos ARIMA.
La elección del IC al 94% implica que se acepta un 6% de probabilidad de que el valor real quede fuera del intervalo, a cambio de una banda ligeramente más estrecha que el IC al 95%.
Si el modelo incluye componente estacional ($D=1$ o $Q \geq 1$), los pronósticos reproducirán la variación trimestral, lo que los hace más realistas para planificación de ingresos.

Comparación: observado vs. pronóstico por trimestre

# Promedio histórico por trimestre
prom_trim <- tapply(ingresos, rep(1:4, times = 12), mean)
cat("Promedio histórico por trimestre (2014–2025):\n")

## Promedio histórico por trimestre (2014–2025):

for (i in 1:4) {
  cat(sprintf("  Q%d: %.2f millones\n", i, prom_trim[i]))
}

##   Q1: 32.58 millones
##   Q2: 60.22 millones
##   Q3: 58.37 millones
##   Q4: 36.98 millones

cat("\nPronosticos 2026 por trimestre:\n")

## 
## Pronosticos 2026 por trimestre:

for (i in 1:4) {
  cat(sprintf("  2026-Q%d: %.2f  [%.2f, %.2f]\n",
              i,
              tabla_pron$Pronóstico[i],
              tabla_pron$LI_94[i],
              tabla_pron$LS_94[i]))
}

series

Salome M.

2026-06-01

Datos: Ventas Mensuales Murphy Brothers (1992–1995)

Configuración de la Serie de Tiempo con ts()

Análisis de la Serie de Tiempo (AST)

Gráfico de la Serie Original

Gráfico Estacional (por mes)

Descomposición STL

Patrones Identificados

Comprobación de Estacionariedad y Estacionalidad

ACF y PACF de la Serie Original

Test de Raíz Unitaria — ADF (Augmented Dickey-Fuller)

Test KPSS

Número de Diferencias Necesarias

Verificación tras Diferenciar

Resumen de Tests

Identificación del Modelo con auto.arima()

Ecuaciones de los Modelos

Ecuación Completa del Modelo Seleccionado

Coeficientes Estimados

Tipo de Serie: Estacionaria vs No Estacionaria

Medias Mensuales (Perfil Estacional)

Crecimiento Anual

Descomposición STL

Respuesta a la Pregunta 1

Gráfico Comparativo

Perfil Estacional Comparado

Correlación entre ambas series

Tabla Comparativa de Patrones

Respuesta a la Pregunta 2

Análisis de Suficiencia de Datos

Evaluación de Componentes para el Modelo

Modelo Recomendado y su Ecuación

Gráfico: Ajuste del Modelo vs. Datos Reales

Bondad de Ajuste

Respuesta a la Pregunta 3

Ecuación Final del Modelo Recomendado

Resumen de Respuestas

Serie de Tiempo y Modelo

Ecuación del Modelo

Metodología Box-Jenkins — Etapas

Diagnóstico Final del Modelo

Pronósticos — 12 Periodos con IC al 94%

Tabla de Pronósticos

Gráfico de Pronósticos

Gráfico Principal con IC 94%

Zoom: Solo Pronósticos 2026

Interpretación del Intervalo de Confianza al 94%

Conclusión

Datos y Primeras Diferencias

Carga de datos

Primeras diferencias

Gráfica de la serie original y diferenciada

Análisis de la Serie de Tiempo (AST)

Configuración con ts() — año final 2025

Gráfica de la serie con autoplot

Descomposición y patrones (AST)

Prueba de Estacionariedad

Prueba Aumentada de Dickey-Fuller (ADF)

Prueba KPSS

ADF sobre primeras diferencias

Función de Autocorrelación (ACF) y Autocorrelación Parcial (PACF)

Identificación automática del modelo con auto.arima()

Ecuación del modelo deducido

Diagnóstico de residuos

Pronósticos con Metodología Box-Jenkins — IC al 94%

Pronósticos para el período 2005–2025

Tabla de pronósticos

Gráfica de pronósticos

Gráfica completa: observado + pronóstico

Resumen y Conclusiones

Tendencia en la serie original

Datos del Dominion Bank

Carga y organización de datos

Configuración como serie de tiempo trimestral

Inciso a) — Autocorrelaciones para los primeros 6 retrasos

Cálculo manual de las autocorrelaciones

Inciso b) — Gráfica de la serie y diagnóstico de estacionariedad

Gráfica de la serie original

Configuración de la Serie de Tiempo con `ts()`

Identificación del Modelo con `auto.arima()`

Configuración con `ts()` — año final 2025

Gráfica de la serie con `autoplot`

Identificación automática del modelo con `auto.arima()`

Identificación del modelo con `auto.arima()`

Inciso i) — Configuración con `ts()` y AST completo

Gráfica con `autoplot`

Identificación automática con `auto.arima()`

Inciso i) — Configuración `ts()`, año final 2025, AST completo

Identificación del modelo con `auto.arima()`