Instrucciones generales: el examen es de carácter individual y consta de tres series, cada una con las mismas diez tareas (1 punto por tarea, 30 puntos en total). Para cada serie se proporciona únicamente la sintaxis para cargarla en R — todo el análisis debe ser generado por el estudiante. Entrega: debe adjuntarse todo el código de R utilizado, además de presentar los resultados (salidas, gráficos e interpretaciones) de cada tarea. Instrucciones generales: toda comparación de estadísticos (AIC, BIC, RMSE, MAE, u otros) debe presentarse en un cuadro (tabla), no únicamente en texto corrido.

library(forecast)
## Warning: package 'forecast' was built under R version 4.5.3
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(prophet)
## Warning: package 'prophet' was built under R version 4.5.3
## Cargando paquete requerido: Rcpp
## Cargando paquete requerido: rlang
library(knitr)

Serie 1: JohnsonJohnson

Utilidades trimestrales por acción de la empresa Johnson & Johnson, entre 1960 y 1980. Es un ejemplo clásico de serie financiera con crecimiento compuesto: en épocas de expansión, las utilidades no solo crecen, sino que tienden a crecer a una tasa cada vez mayor — por eso este tipo de series suele modelarse mejor en escala logarítmica que en su escala original. Es una serie ampliamente citada en la literatura de series de tiempo precisamente por combinar ese crecimiento exponencial con una estacionalidad trimestral persistente, ligada al calendario de reportes financieros de la empresa.

data("JohnsonJohnson")
serie <- JohnsonJohnson
serie
##       Qtr1  Qtr2  Qtr3  Qtr4
## 1960  0.71  0.63  0.85  0.44
## 1961  0.61  0.69  0.92  0.55
## 1962  0.72  0.77  0.92  0.60
## 1963  0.83  0.80  1.00  0.77
## 1964  0.92  1.00  1.24  1.00
## 1965  1.16  1.30  1.45  1.25
## 1966  1.26  1.38  1.86  1.56
## 1967  1.53  1.59  1.83  1.86
## 1968  1.53  2.07  2.34  2.25
## 1969  2.16  2.43  2.70  2.25
## 1970  2.79  3.42  3.69  3.60
## 1971  3.60  4.32  4.32  4.05
## 1972  4.86  5.04  5.04  4.41
## 1973  5.58  5.85  6.57  5.31
## 1974  6.03  6.39  6.93  5.85
## 1975  6.93  7.74  7.83  6.12
## 1976  7.74  8.91  8.28  6.84
## 1977  9.54 10.26  9.54  8.73
## 1978 11.88 12.06 12.15  8.91
## 1979 14.04 12.96 14.85  9.99
## 1980 16.20 14.67 16.02 11.61

Tarea 1.

Graficar la serie. Presentar el gráfico de la serie temporal completa e identificar visualmente sus características generales (tendencia, variabilidad, posibles cambios de comportamiento).

par(mfrow = c(2, 1))
plot(serie, type = "o", pch = 20, xlab = "Año", ylab = "USD por acción",
     main = "Utilidades trimestrales por acción de Johnson & Johnson")
plot(log(serie), type = "o", pch = 20, xlab = "Año", ylab = "log(USD)",
     main = "Serie en escala logarítmica")

par(mfrow = c(1, 1))

Interpretación:

La serie crece durante todo el período. Las utilidades pasan de valores entre 0.44 y 0.85 USD por trimestre en 1960 a valores entre 11.61 y 16.20 USD en 1980. El avance es lento en los sesenta y se acelera desde 1970. La curva se vuelve cada vez más empinada. Esto es típico del crecimiento compuesto. La variabilidad aumenta con el nivel. Al principio, las oscilaciones trimestrales casi no se notan. En 1979 y 1980, superan los 4 USD dentro de un mismo año. Las caídas del cuarto trimestre son cada vez más profundas. Por ejemplo, de 12.15 a 8.91 en 1978 y de 14.85 a 9.99 en 1979. No se observan quiebres bruscos. El cambio más visible es esa mayor amplitud al final de la muestra. En escala logarítmica, la tendencia es casi una recta. Esto indica una tasa de crecimiento relativa casi constante. Las oscilaciones dejan de crecer con el nivel. La serie tiene, por tanto, una estructura multiplicativa. Por eso, se justifica la descomposición multiplicativa y el uso del logaritmo en los modelos.

Tarea 2.

Descomposición y comentario. Aplicar decompose() para separar la serie en tendencia, estacionalidad y componente irregular, y comentar brevemente lo que muestra cada componente.

desc <- decompose(serie, type = "multiplicative")
plot(desc)

setNames(round(desc$figure, 3), paste0("Trim", 1:4))
## Trim1 Trim2 Trim3 Trim4 
## 0.993 1.033 1.114 0.860

Interpretación:

Se usó la descomposición multiplicativa porque la amplitud estacional crece con el nivel de la serie.

La tendencia sube todo el tiempo, sin bajar, desde menos de 1 USD a principios de los sesenta. Llega a casi 14 USD al final. En los sesenta, la línea casi no cambia. En los setenta, sube rápido y toma una forma exponencial.

La estacionalidad se repite igual cada año. Según los factores, el tercer trimestre está 11.4 % por encima de la tendencia. El segundo trimestre se ubica 3.3 % arriba. El primero está casi igual que la tendencia, con -0.7 %. El cuarto trimestre queda 14 % por debajo. Así, el cuarto trimestre es el más débil del año.

El componente irregular cambia alrededor de 1. En general, se mantiene entre 0.9 y 1.1, sin mostrar una tendencia. Las desviaciones más grandes, cerca de 0.8 y de 1.2, aparecen en 1960-1961, 1967-1968 y 1978-1980. Esto muestra que el patrón estacional no fue tan estable como supone el método. La tendencia y el irregular no tienen valores en los dos primeros ni en los dos últimos trimestres. Esto se debe a la media móvil centrada.

Tarea 3

Pruebas de estacionariedad y su interpretación. Aplicar dos pruebas formales de estacionariedad sobre la serie original (ADF y KPSS), e interpretar sus resultados: hipótesis nula de cada una, decisión al 5% de significancia, y conclusión sobre si la serie es o no estacionaria.

adf= adf.test(serie)
## Warning in adf.test(serie): p-value greater than printed p-value
kpss = kpss.test(serie, null = "Level")
## Warning in kpss.test(serie, null = "Level"): p-value smaller than printed
## p-value
adf
## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie
## Dickey-Fuller = 1.9321, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary
kpss
## 
##  KPSS Test for Level Stationarity
## 
## data:  serie
## KPSS Level = 1.9804, Truncation lag parameter = 3, p-value = 0.01
pruebas = data.frame(Prueba = c(
  "ADF", "KPSS"), H0 = c("Raíz unitaria (no estacionaria)", "Estacionaria en nivel"),
  Estadistico = unname(c(adf$statistic, kpss$statistic)), p_valor = c(adf$p.value, kpss$p.value))
pruebas$Decision <- ifelse(pruebas$p_valor < 0.05, "Rechazar H0", "No rechazar H0")
kable(pruebas, digits = 3, caption = "Pruebas de estacionariedad sobre la serie original")
Pruebas de estacionariedad sobre la serie original
Prueba H0 Estadistico p_valor Decision
ADF Raíz unitaria (no estacionaria) 1.932 0.99 No rechazar H0
KPSS Estacionaria en nivel 1.980 0.01 Rechazar H0

Interpretación:

ADF: la hipótesis nula dice que la serie tiene raíz unitaria. Esto significa que no es estacionaria. El estadístico es positivo (1.932). El p-valor es mayor que 0.99. R indica que el valor real es más alto que el mostrado. Por eso, al 5 % no se rechaza H0. Esto pasa aunque la regresión del ADF incluye constante y tendencia lineal.

KPSS: la hipótesis nula dice que la serie es estacionaria en nivel. El estadístico es 1.980 y el p-valor es menor que 0.01. Por eso, al 5 %, se rechaza H0.

Ambas pruebas llevan a la misma conclusión. La serie original no es estacionaria. Esto concuerda con la tendencia creciente y la varianza que aumenta con el nivel, vistas en las tareas 1 y 2.

Tarea 4.

Correlogramas y su interpretación. Presentar el ACF y el PACF de la serie, y comentar los patrones observados (decaimiento, cortes abruptos, etc.).

par(mfrow = c(1, 2))
acf_s <- Acf(serie, lag.max = 24, main = "ACF")
pacf_s <- Pacf(serie, lag.max = 24, main = "PACF")

par(mfrow = c(1, 1))
acf_s
## 
## Autocorrelations of series 'serie', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.925 0.888 0.833 0.824 0.764 0.718 0.675 0.654 0.608 0.564 0.526 0.500 
##    13    14    15    16    17    18    19    20    21    22    23    24 
## 0.456 0.423 0.390 0.370 0.334 0.307 0.278 0.260 0.226 0.199 0.172 0.153
pacf_s
## 
## Partial autocorrelations of series 'serie', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.925  0.225 -0.089  0.265 -0.257 -0.085  0.139  0.004 -0.104 -0.011  0.021 
##     12     13     14     15     16     17     18     19     20     21     22 
## -0.015 -0.089  0.078  0.005 -0.036 -0.008  0.011 -0.034  0.011 -0.052 -0.014 
##     23     24 
## -0.009 -0.016

Interpretación:

La ACF baja muy despacio. Empieza en 0.925 en el rezago 1 y sigue fuera de las bandas de confianza (±0.214) hasta el rezago 21. En el rezago 24 todavía es 0.153. Ese descenso lento, casi recto, es típico de una serie no estacionaria dominada por la tendencia. La estacionalidad casi no se nota. En los rezagos 4, 8, 12 y 16 el descenso se detiene un poco, pero no hay picos estacionales claros porque la tendencia los tapa.

La PACF muestra un corte abrupto después del rezago 1 (0.925). Solo los rezagos 2 (0.225, apenas), 4 (0.265) y 5 (-0.257) salen de las bandas. El par 4-5 muestra la relación con el mismo trimestre del año anterior. Los otros valores son pequeños. Una ACF que no desaparece junto con una PACF cuyo primer valor está cerca de 1 confirma lo que mostraron ADF y KPSS. La serie necesita diferenciación, regular y probablemente estacional, antes de buscar términos AR o MA.

Tarea 5.

Ajuste SARIMA/ARIMA. Ajustar un modelo mediante auto.arima() y reportar el orden (p,d,q)(P,D,Q) de cada componente junto con los coeficientes estimados.

modelo_arima= auto.arima(serie, lambda = 0, stepwise = FALSE, approximation = FALSE)
modelo_arima
## Series: serie 
## ARIMA(0,0,2)(0,1,0)[4] with drift 
## Box Cox transformation: lambda= 0 
## 
## Coefficients:
##          ma1     ma2   drift
##       0.2398  0.3616  0.0384
## s.e.  0.1076  0.1089  0.0039
## 
## sigma^2 = 0.00787:  log likelihood = 81.65
## AIC=-155.29   AICc=-154.76   BIC=-145.76
arimaorder(modelo_arima)
##         p         d         q         P         D         Q Frequency 
##         0         0         2         0         1         0         4

Interpretación:

Con la serie transformada a logaritmos (lambda = 0) y una búsqueda completa, auto.arima() eligió por AICc un modelo SARIMA(0,0,2)(0,1,0)[4] con deriva. La parte regular tiene p = 0, d = 0 y q = 2. Incluye dos términos de media móvil y ninguna diferencia regular. La parte estacional tiene P = 0, D = 1 y Q = 0. Incluye una diferencia estacional de período 4, sin términos AR ni MA estacionales. Esta diferencia compara cada trimestre con el mismo trimestre del año anterior. Como (1 - B^4) contiene el factor (1 - B), también aplica una diferencia regular. El modelo no añade otra diferencia regular.

Los coeficientes son ma1 = 0.2398 (error estándar = 0.1076), ma2 = 0.3616 (error estándar = 0.1089) y deriva = 0.0384 (error estándar = 0.0039). Los tres resultan significativos al 5 %, ya que sus valores absolutos son mayores que el doble de sus errores estándar.

La ecuación estimada es:

\[(1-B^4)\log Y_t = 0.1536 + \varepsilon_t + 0.2398\,\varepsilon_{t-1} + 0.3616\,\varepsilon_{t-2}\]

donde \(B^4\) es el rezago de cuatro trimestres y \(\varepsilon_t\) es el error del modelo.

El valor 0.1536 es igual a 4 × 0.0384 y muestra el crecimiento interanual medio en escala logarítmica. Esto equivale a cerca de 16.6 % anual en la escala original.

La varianza residual es 0.00787 en escala logarítmica. Esto equivale a un error típico aproximado del 9 %. Los criterios de información son AIC = -155.29, AICc = -154.76 y BIC = -145.76.

Tarea 6.

Ajuste de Prophet. Ajustar Prophet sobre la misma serie, reportando sus parámetros o componentes principales.

set.seed(2026)
fechas= seq(as.Date("1960-01-01"), by = "quarter", 
            length.out = length(serie))
datos = data.frame(ds = fechas, y = log(as.numeric(serie)))

ajustar_prophet = function(d) {capture.output(m <- prophet(d, yearly.seasonality = 2,
                                                            weekly.seasonality = FALSE,
                                                            daily.seasonality = FALSE)) 
  m}
modelo_prophet = ajustar_prophet(datos)
pred_prophet = predict(modelo_prophet, datos)
prophet_plot_components(modelo_prophet, pred_prophet)

round(exp(tapply(pred_prophet$yearly, cycle(serie), mean)), 3)
##     1     2     3     4 
## 1.009 1.040 1.117 0.855

Interpretación:

Prophet se ajustó al logaritmo de las utilidades con una tendencia lineal por tramos y estacionalidad anual de Fourier de orden 2. Esto es suficiente para mostrar el patrón trimestral. Se desactivaron las estacionalidades semanal y diaria porque los datos son trimestrales. Se usaron los valores por defecto: 25 posibles puntos de cambio en el primer 80 % de la muestra, changepoint.prior.scale = 0.05 y seasonality.prior.scale = 10. Al modelar en logaritmos, la estacionalidad aditiva se convierte en un efecto multiplicativo sobre las utilidades originales.

La tendencia es casi lineal. Muestra cambios de pendiente alrededor de 1964 y 1973. El crecimiento se acelera a mediados de los sesenta. Se modera un poco durante los setenta. La estacionalidad anual solo se interpreta en las fechas observadas: 1 de enero, abril, julio y octubre.

Los factores estacionales son 1.009 para el primer trimestre, 1.040 para el segundo, 1.117 para el tercero y 0.855 para el cuarto. Si los comparamos con la tendencia, equivalen a utilidades que son aproximadamente 0.9 %, 4.0 % y 11.7 % más altas en los primeros tres trimestres, y 14.5 % más bajas en el cuarto. El efecto del cuarto trimestre es parecido al que obtuvo decompose() en la Tarea 2.

Tarea 7.

Comparación de métricas y selección del modelo. Construir un cuadro comparativo de ambos modelos con AIC, BIC, RMSE y MAE, señalando explícitamente cuáles de estas métricas son comparables entre las dos familias de modelo y cuáles no, y seleccionar el modelo preferido justificando la elección.

train = window(serie, end = c(1979, 4))
test = window(serie, start = c(1980, 1))

arima_train = auto.arima(train, lambda = 0, stepwise = FALSE, approximation = FALSE)
arimaorder(arima_train)
##         p         d         q         P         D         Q Frequency 
##         0         0         2         0         1         0         4
pred_arima = forecast(arima_train, h = 4)$mean

prophet_train= ajustar_prophet(datos[1:80, ])
pred_prophet_test = exp(predict(prophet_train, datos[81:84, ])$yhat)

e_arima= as.numeric(test - pred_arima)
e_prophet = as.numeric(test) - pred_prophet_test

cuadro = data.frame(Modelo = c("auto.arima (log)", "Prophet (log)"),
                    AIC = c(modelo_arima$aic, NA), BIC = c(modelo_arima$bic, NA),
                    RMSE = c(sqrt(mean(e_arima^2)), sqrt(mean(e_prophet^2))), 
                    MAE = c(mean(abs(e_arima)), mean(abs(e_prophet))))
kable(cuadro, digits = 3, caption = "Comparación de modelos (RMSE y MAE: pronóstico de 1980)")
Comparación de modelos (RMSE y MAE: pronóstico de 1980)
Modelo AIC BIC RMSE MAE
auto.arima (log) -155.291 -145.763 0.697 0.495
Prophet (log) NA NA 1.558 1.452

Interpretacion. Modelo seleccionado

Ambos modelos se ajustaron otra vez con los datos de 1960 a 1979. Pronosticaron los cuatro trimestres de 1980. El SARIMA mantuvo el orden (0,0,2)(0,1,0)[4]. El AIC y el BIC corresponden al SARIMA de la Tarea 5. El RMSE y el MAE se calcularon con los pronósticos de 1980 y están expresados en USD por acción.

Comparabilidad de las métricas entre SARIMA y Prophet
Métrica ¿Comparable? Motivo
AIC No Prophet no lo reporta: se estima por máximo a posteriori y no tiene un número fijo de parámetros.
BIC No Mismo motivo que el AIC.
RMSE Sí Mismos cuatro trimestres de 1980, en USD por acción y con los mismos datos de entrenamiento.
MAE Sí Igual que el RMSE.

La comparación de los errores favorece al SARIMA en esos cuatro trimestres. Obtuvo un RMSE de 0.697 frente a 1.558 de Prophet. También logró un MAE de 0.495 frente a 1.452 USD por acción. El SARIMA proyecta desde lo observado en 1979 e incluye la deriva. Prophet extrapola la tendencia y el patrón estacional estimados con toda la historia.

Se elige el SARIMA(0,0,2)(0,1,0)[4] con deriva, en logaritmos. Es un modelo simple con tres coeficientes importantes. Sin embargo, la evaluación solo incluye cuatro observaciones. Por eso, esta diferencia no es suficiente para decir que el SARIMA será mejor en general. La Tarea 8 revisará la elección usando validación por origen móvil.

Tarea 8.

Validación del modelo seleccionado. Evaluar el modelo seleccionado mediante una partición train-test o, preferiblemente, validación cruzada por origen móvil.

f_sarima = function(x, h) {modelo = Arima(x, order = c(0, 0, 2), seasonal = c(0, 1, 0), include.drift = TRUE, lambda = 0)
forecast(modelo, h = h)}

e_cv = tsCV(train, f_sarima, h = 4, initial = 39)

validacion = data.frame(Horizonte = 1:4, RMSE = sqrt(colMeans(e_cv^2, na.rm = TRUE)),
                        MAE = colMeans(abs(e_cv), na.rm = TRUE), Pronosticos = colSums(!is.na(e_cv)))
kable(validacion, digits = 3, row.names = FALSE, caption = "Validación por origen móvil del SARIMA (1970-1979)")
Validación por origen móvil del SARIMA (1970-1979)
Horizonte RMSE MAE Pronosticos
1 0.540 0.430 40
2 0.537 0.440 39
3 0.558 0.449 38
4 0.558 0.446 37

Interpretacion:

La validación reestimó el SARIMA(0,0,2)(0,1,0)[4] con deriva en 40 orígenes sucesivos. Usó una ventana expansiva desde 40 trimestres, de 1960 a 1969. Pronosticó de 1 a 4 trimestres adelante dentro de 1970 a 1979. El año 1980 quedó fuera porque ya se había usado para elegir el modelo. El RMSE se mantiene entre 0.537 y 0.558 USD, y el MAE entre 0.430 y 0.449 USD en los cuatro horizontes. El error casi no crece al pronosticar más lejos. Estos errores son del mismo orden que los de la prueba de 1980, con RMSE de 0.697 y MAE de 0.495. En 1980 son algo mayores porque el nivel de la serie es más alto. El desempeño del modelo es estable y no depende de un solo año.

Tarea 9.

Evaluación de residuales. Presentar la evaluación de residuales del modelo seleccionado, incluyendo una prueba de normalidad y la prueba de Ljung-Box, interpretando ambos resultados.

checkresiduals(modelo_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,2)(0,1,0)[4] with drift
## Q* = 6.9236, df = 6, p-value = 0.328
## 
## Model df: 2.   Total lags used: 8
shapiro.test(residuals(modelo_arima))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modelo_arima)
## W = 0.98065, p-value = 0.2383

Interpretacion:

Los residuales, en escala logarítmica, fluctúan alrededor de cero sin mostrar un patrón visible. Su dispersión es un poco mayor antes de 1973. El valor más extremo, cerca de -0.28, aparece en el primer trimestre de 1961. Los cuatro primeros valores son casi cero porque coinciden con el inicio de la diferencia estacional. En la ACF, solo el rezago 7 sale apenas de las bandas. Esto puede pasar por azar al revisar muchos rezagos.

Ljung-Box: la hipótesis nula dice que no hay autocorrelación en los residuales hasta el rezago 8. Con Q* igual a 6.924, 6 grados de libertad y p igual a 0.328, no se rechaza H0 al 5 %.

Shapiro-Wilk: la hipótesis nula dice que los residuales tienen una distribución normal. Con W = 0.981 y p = 0.238, no se rechaza H0 al 5 %. El histograma es casi simétrico. El valor de 1961 está en la cola izquierda.

Los residuales actúan como ruido blanco con una distribución que es casi normal. Así, el modelo recoge bien la dependencia de la serie.

Tarea 10.

Pronóstico y comentario final. Generar el pronóstico del modelo seleccionado y comentar brevemente qué tan razonable resulta, considerando lo observado en la serie.

pronostico = forecast(modelo_arima, h = 4, level = c(80, 95))
pronostico
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1981 Q1       18.50781 16.51885 20.73625 15.55400 22.02257
## 1981 Q2       17.21349 15.31420 19.34832 14.39514 20.58362
## 1981 Q3       18.68227 16.50472 21.14713 15.45669 22.58098
## 1981 Q4       13.53940 11.96128 15.32573 11.20176 16.36487
plot(pronostico, include = 24, xlab = "Año", ylab = "USD por acción")

Interpretacion:

El modelo predice para 1981 utilidades de 18.51, 17.21, 18.68 y 13.54 USD por acción en los trimestres 1 a 4. Cada valor supera al del mismo trimestre de 1980 entre 14 % y 17 %. Esto sigue la tendencia estimada, cerca de 16.6 % anual. El patrón estacional reciente se mantiene: primer y tercer trimestre altos, segundo algo menor, y una caída marcada en el cuarto. Los intervalos del 95 % van de unos 17 % por debajo a unos 20 % por encima de cada pronóstico. Son asimétricos porque el modelo se ajustó en logaritmo. El pronóstico tiene sentido porque mantiene el crecimiento compuesto y la estacionalidad de la serie, pero solo a corto plazo, ya que supone que el ritmo de crecimiento de 1960-1980 sigue igual.

Serie 2: USAccDeaths

Número mensual de muertes accidentales en Estados Unidos, registradas entre 1973 y 1978. Es una serie de vigilancia: los accidentes (de tránsito, laborales, domésticos, etc.) no ocurren de manera uniforme a lo largo del año — factores como el clima y la mayor actividad recreativa al aire libre generan un componente estacional marcado, con picos característicos en los meses de verano. Es una serie corta y relativamente estable, sin una tendencia dominante, útil para practicar el flujo completo de modelado sobre un caso más sencillo que el de las otras dos series.

data("USAccDeaths")
serie = USAccDeaths
serie
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 1973  9007  8106  8928  9137 10017 10826 11317 10744  9713  9938  9161  8927
## 1974  7750  6981  8038  8422  8714  9512 10120  9823  8743  9129  8710  8680
## 1975  8162  7306  8124  7870  9387  9556 10093  9620  8285  8466  8160  8034
## 1976  7717  7461  7767  7925  8623  8945 10078  9179  8037  8488  7874  8647
## 1977  7792  6957  7726  8106  8890  9299 10625  9302  8314  8850  8265  8796
## 1978  7836  6892  7791  8192  9115  9434 10484  9827  9110  9070  8633  9240

Tarea 11.

Graficar la serie. Presentar el gráfico de la serie temporal completa e identificar visualmente sus características generales (tendencia, variabilidad, posibles cambios de comportamiento).

plot(serie, type = "o", pch = 20, xlab = "Año", ylab = "Número de muertes", main = "Muertes accidentales mensuales en EE. UU.")

Interpretación:

La serie no tiene una tendencia dominante: los valores van de 6892 muertes (febrero de 1978) a 11317 (julio de 1973). El nivel baja entre 1973 y 1974 (julio, por ejemplo, pasa de 11317 a 10120 muertes), llega a su punto más bajo en 1976 y se recupera un poco en 1977 y 1978 sin volver a los valores de 1973. Lo que domina es una estacionalidad muy marcada y regular: cada año el máximo ocurre en julio y el mínimo en febrero. La diferencia entre ambos meses se mantiene entre unas 2600 y 3700 muertes y no crece con el nivel, lo que indica una estructura aditiva. El único cambio de comportamiento visible es el descenso de nivel entre 1973 y 1974.

Tarea 12.

Descomposición y comentario. Aplicar decompose() para separar la serie en tendencia, estacionalidad y componente irregular, y comentar brevemente lo que muestra cada componente.

desc = decompose(serie, type = "additive")
plot(desc)

setNames(round(desc$figure), month.abb)
##   Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec 
##  -806 -1523  -741  -515   340   745  1679   986  -109   264  -261   -59

Interpretación:

Se usó la descomposición aditiva porque la amplitud estacional es estable y no depende del nivel de la serie.

La tendencia baja de unas 9600 muertes a mediados de 1973 a unas 8700 a mediados de 1974, sigue descendiendo con suavidad hasta un mínimo cercano a 8350 en 1976 y luego sube hasta unas 8800 a mediados de 1978. Sus cambios son menores que las oscilaciones estacionales.

La estacionalidad deja las muertes por encima de la tendencia de mayo a agosto, con el máximo en julio (+1679), y por debajo de enero a abril, con el mínimo en febrero (-1523); entre ambos meses hay unas 3200 muertes de diferencia. Octubre muestra un pequeño repunte (+264), y septiembre, noviembre y diciembre quedan cerca de la tendencia.

El componente irregular oscila alrededor de cero, casi siempre entre -400 y 400 muertes, sin patrón claro. Los mayores desvíos ocurren a inicios de 1974, cerca de -500, y en febrero de 1976, cerca de 580. La tendencia y el irregular no tienen valores en los seis primeros ni en los seis últimos meses por la media móvil centrada.

Tarea 13.

Pruebas de estacionariedad y su interpretación. Aplicar dos pruebas formales de estacionariedad sobre la serie original (ADF y KPSS), e interpretar sus resultados: hipótesis nula de cada una, decisión al 5% de significancia, y conclusión sobre si la serie es o no estacionaria.

adf = adf.test(serie)
kpss = kpss.test(serie, null = "Level")
## Warning in kpss.test(serie, null = "Level"): p-value greater than printed
## p-value
adf
## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie
## Dickey-Fuller = -3.8221, Lag order = 4, p-value = 0.02268
## alternative hypothesis: stationary
kpss
## 
##  KPSS Test for Level Stationarity
## 
## data:  serie
## KPSS Level = 0.19799, Truncation lag parameter = 3, p-value = 0.1
pruebas = data.frame(Prueba = c("ADF", "KPSS"), H0 = c("Raíz unitaria (no estacionaria)", "Estacionaria en nivel"), 
                     Estadistico = unname(c(adf$statistic, kpss$statistic)), p_valor = c(adf$p.value, kpss$p.value))

pruebas$Decision = ifelse(pruebas$p_valor < 0.05, "Rechazar H0", "No rechazar H0")
kable(pruebas, digits = 3, caption = "Pruebas de estacionariedad sobre la serie original")
Pruebas de estacionariedad sobre la serie original
Prueba H0 Estadistico p_valor Decision
ADF Raíz unitaria (no estacionaria) -3.822 0.023 Rechazar H0
KPSS Estacionaria en nivel 0.198 0.100 No rechazar H0

Interpretación:

ADF: la hipótesis nula es que la serie tiene raíz unitaria. El estadístico es -3.822 con p = 0.023, menor que 0.05, así que se rechaza H0.

KPSS: la hipótesis nula es que la serie es estacionaria en nivel. El estadístico es 0.198 y el p-valor es mayor que 0.10 (R advierte que el valor real supera al impreso), por lo que no se rechaza H0.

Las dos pruebas coinciden: la serie es estacionaria en nivel, sin raíz unitaria, lo que concuerda con la ausencia de una tendencia dominante en la Tarea 11. Estas pruebas no evalúan la estacionalidad, que sigue presente y deberá modelarse.

Tarea 14.

Correlogramas y su interpretación. Presentar el ACF y el PACF de la serie, y comentar los patrones observados (decaimiento, cortes abruptos, etc.).

par(mfrow = c(1, 2))
acf_s = Acf(serie, lag.max = 36, main = "ACF")
pacf_s = Pacf(serie, lag.max = 36, main = "PACF")

par(mfrow = c(1, 1))
acf_s
## 
## Autocorrelations of series 'serie', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.707  0.409  0.084 -0.182 -0.294 -0.423 -0.346 -0.285 -0.065  0.162 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.414  0.629  0.429  0.221 -0.038 -0.216 -0.278 -0.362 -0.268 -0.227 -0.054 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.116  0.307  0.450  0.285  0.119 -0.123 -0.239 -0.311 -0.359 -0.281 -0.242 
##     33     34     35     36 
## -0.095  0.023  0.179  0.302
pacf_s
## 
## Partial autocorrelations of series 'serie', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.707 -0.184 -0.266 -0.167  0.024 -0.284  0.148 -0.149  0.277  0.043  0.409 
##     12     13     14     15     16     17     18     19     20     21     22 
##  0.129 -0.488 -0.042  0.101 -0.011  0.093 -0.055  0.126 -0.155  0.037  0.032 
##     23     24     25     26     27     28     29     30     31     32     33 
##  0.002 -0.090  0.025 -0.076 -0.036  0.014 -0.121  0.025 -0.112  0.056 -0.089 
##     34     35     36 
## -0.083 -0.030  0.103

Interpretación:

La ACF tiene forma de onda con periodo de 12 meses: es positiva en los rezagos 12 (0.629), 24 (0.450) y 36 (0.302) y negativa en 6 (-0.423), 18 (-0.362) y 30 (-0.359), todos fuera de las bandas (±0.231). Los picos estacionales disminuyen lentamente, señal de una estacionalidad persistente que conviene tratar con una diferencia estacional. En los primeros rezagos la ACF cae rápido (0.707, 0.409 y 0.084 en los rezagos 1 a 3), lo que confirma que no hay una tendencia fuerte.

La PACF tiene un pico grande en el rezago 1 (0.707) y otros significativos en los rezagos 3, 6, 9, 11 y 13, los mayores alrededor del ciclo anual (0.409 en el 11 y -0.488 en el 13). Desde el rezago 14 todos los valores quedan dentro de las bandas. Esto apunta a una dependencia de corto plazo de tipo autorregresivo junto con una estructura estacional de periodo 12.

Tarea 15.

Ajuste SARIMA/ARIMA. Ajustar un modelo mediante auto.arima() y reportar el orden (p,d,q)(P,D,Q) de cada componente junto con los coeficientes estimados.

modelo_arima = auto.arima(serie, stepwise = FALSE, approximation = FALSE)
modelo_arima
## Series: serie 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4303  -0.5528
## s.e.   0.1228   0.1784
## 
## sigma^2 = 102860:  log likelihood = -425.44
## AIC=856.88   AICc=857.32   BIC=863.11
arimaorder(modelo_arima)
##         p         d         q         P         D         Q Frequency 
##         0         1         1         0         1         1        12

Interpretación:

auto.arima() eligió un SARIMA(0,1,1)(0,1,1)[12]. En la parte regular p = 0, d = 1 y q = 1: una diferencia regular y un término de media móvil. En la parte estacional P = 0, D = 1 y Q = 1: una diferencia estacional de periodo 12 y un término de media móvil estacional. La diferencia estacional recoge el patrón anual que mostró la ACF. La regular la decidió auto.arima con su prueba KPSS sobre la serie ya diferenciada estacionalmente, por lo que no contradice la Tarea 13, que evaluó la serie original.

Los coeficientes son ma1 = -0.4303 (e.e. 0.1228) y sma1 = -0.5528 (e.e. 0.1784); ambos son significativos al 5 %, porque superan el doble de su error estándar. El modelo estimado es

\[(1-B)(1-B^{12})\,Y_t = (1-0.4303\,B)(1-0.5528\,B^{12})\,\varepsilon_t\]

donde \(B\) es el operador de rezago y \(\varepsilon_t\) el error del modelo. La varianza residual es 102860, un error típico de unas 321 muertes por mes, y los criterios de información son AIC = 856.88, AICc = 857.32 y BIC = 863.11.

Tarea 16.

Ajuste de Suavizamiento Exponencial (Holt-Winters). Ajustar Suavizamiento Exponencial (Holt-Winters) sobre la misma serie, reportando sus parámetros o componentes principales.

modelo_hw = hw(serie, seasonal = "additive")
modelo_hw$model
## Holt-Winters' additive method 
## 
## Call:
## hw(y = serie, seasonal = "additive")
## 
##   Smoothing parameters:
##     alpha = 0.5378 
##     beta  = 0.0012 
##     gamma = 0.0037 
## 
##   Initial states:
##     l = 9933.1305 
##     b = -20.0469 
##     s = 58.2216 -260.4927 230.8796 -47.9817 988.7754 1698.957
##            751.926 333.9133 -514.4812 -741.2456 -1510.742 -987.7303
## 
##   sigma:  301.4197
## 
##      AIC     AICc      BIC 
## 1145.850 1157.183 1184.553
plot(modelo_hw$model)

Interpretación:

Se ajustó el Holt-Winters aditivo, coherente con la amplitud estacional estable de la Tarea 12. Los parámetros de suavizamiento son α = 0.5378, β = 0.0012 y γ = 0.0037. El nivel se actualiza con rapidez moderada, con algo más de la mitad del peso en la observación más reciente, mientras que la pendiente y la estacionalidad casi no cambian en el tiempo. El nivel inicial es 9933 muertes y la pendiente inicial, -20.05 muertes por mes. Los estados estacionales iniciales, que se leen de diciembre hacia enero, van de -1511 en febrero a +1699 en julio, casi los mismos efectos que dio decompose(). La desviación estándar residual es 301.42 muertes, con AIC = 1145.85, AICc = 1157.18 y BIC = 1184.55.

En el gráfico de componentes, el nivel baja de unas 9900 muertes en 1973 a unas 8700 en 1974, llega a su mínimo en 1976 y sube hasta unas 9000 al final de 1978. La pendiente se mantiene cerca de -20 y la estacionalidad se repite prácticamente igual cada año.

Tarea 17.

Comparación de métricas y selección del modelo. Construir un cuadro comparativo de ambos modelos con AIC, BIC, RMSE y MAE, señalando explícitamente cuáles de estas métricas son comparables entre las dos familias de modelo y cuáles no, y seleccionar el modelo preferido justificando la elección.

train = window(serie, end = c(1977, 12))
test = window(serie, start = c(1978, 1))

arima_train = auto.arima(train, stepwise = FALSE, approximation = FALSE)
arimaorder(arima_train)
##         p         d         q         P         D         Q Frequency 
##         0         1         1         0         1         1        12
pred_arima = forecast(arima_train, h = 12)$mean
pred_hw = hw(train, seasonal = "additive", h = 12)$mean

e_arima = as.numeric(test - pred_arima)
e_hw = as.numeric(test - pred_hw)

cuadro = data.frame(Modelo = c("auto.arima", "Holt-Winters aditivo"),
                    AIC = c(modelo_arima$aic, modelo_hw$model$aic),
                    BIC = c(modelo_arima$bic, modelo_hw$model$bic),
                    RMSE = c(sqrt(mean(e_arima^2)), sqrt(mean(e_hw^2))),
                    MAE = c(mean(abs(e_arima)), mean(abs(e_hw))))
kable(cuadro, digits = 2,
      caption = "Comparación de modelos (RMSE y MAE: pronóstico de 1978)")
Comparación de modelos (RMSE y MAE: pronóstico de 1978)
Modelo AIC BIC RMSE MAE
auto.arima 856.88 863.11 288.83 231.61
Holt-Winters aditivo 1145.85 1184.55 414.41 346.63

Interpretación:

Ambos modelos se reajustaron con los datos de 1973 a 1977 y pronosticaron los 12 meses de 1978; el SARIMA reajustado conserva el orden (0,1,1)(0,1,1)[12]. El AIC y el BIC corresponden a los ajustes con toda la serie (tareas 15 y 16), y el RMSE y el MAE, a los pronósticos de 1978, en número de muertes.

Comparabilidad de las métricas entre SARIMA y Holt-Winters
Métrica ¿Comparable? Motivo
AIC No El del SARIMA usa la verosimilitud de la serie diferenciada (59 observaciones) y el del Holt-Winters la de las 72 observaciones con otra formulación, así que la diferencia no indica cuál es mejor.
BIC No Mismo motivo que el AIC.
RMSE Sí Mismos 12 meses de 1978, en número de muertes y con los mismos datos de entrenamiento.
MAE Sí Igual que el RMSE.

Modelo seleccionado: SARIMA(0,1,1)(0,1,1)[12]. Sus errores de pronóstico en 1978 son entre 30 % y 33 % menores que los del Holt-Winters: RMSE de 288.83 frente a 414.41 y MAE de 231.61 frente a 346.63 muertes. Además es más parsimonioso, con 2 coeficientes frente a los 16 parámetros del Holt-Winters. El Holt-Winters proyecta una pendiente negativa casi fija (cerca de -20 muertes por mes, Tarea 16), mientras que en 1978 la serie fue, en promedio, más alta que en 1977. La Tarea 18 verifica la elección con validación por origen móvil.

Tarea 18.

Validación del modelo seleccionado. Evaluar el modelo seleccionado mediante una partición train-test o, preferiblemente, validación cruzada por origen móvil.

f_sarima = function(x, h) {
  modelo = Arima(x, order = c(0, 1, 1), seasonal = c(0, 1, 1))
  forecast(modelo, h = h)}
e_cv = tsCV(train, f_sarima, h = 12, initial = 35)

validacion = data.frame(Horizonte = 1:12,
                        RMSE = sqrt(colMeans(e_cv^2, na.rm = TRUE)),
                        MAE = colMeans(abs(e_cv), na.rm = TRUE),
                        Pronosticos = colSums(!is.na(e_cv)))
kable(validacion, digits = 1, row.names = FALSE,
      caption = "Validación por origen móvil del SARIMA (1976-1977)")
Validación por origen móvil del SARIMA (1976-1977)
Horizonte RMSE MAE Pronosticos
1 367.4 280.0 24
2 399.5 307.6 23
3 361.0 279.3 22
4 408.1 330.1 21
5 368.3 285.7 20
6 447.3 357.2 19
7 498.7 426.3 18
8 517.6 436.7 17
9 543.4 482.8 16
10 557.3 474.1 15
11 612.0 518.2 14
12 702.8 647.8 13

Interpretación:

La validación reestimó el SARIMA(0,1,1)(0,1,1)[12] en 24 orígenes sucesivos, con ventana expansiva desde 36 meses (1973-1975), y pronosticó de 1 a 12 meses adelante dentro de 1976-1977. El año 1978 quedó fuera porque ya se usó para elegir el modelo. Hasta cinco meses adelante el error se mantiene estable, con RMSE entre 361 y 408 y MAE entre 279 y 330 muertes; desde el sexto mes crece de forma sostenida hasta un RMSE de 702.8 y un MAE de 647.8 a doce meses. Estos errores son mayores que los de la prueba de 1978 (RMSE 288.83 y MAE 231.61), en parte porque los primeros orígenes se estimaron con solo tres años de datos. El modelo pronostica bien a corto plazo y pierde precisión al acercarse al año completo.

Tarea 19.

Evaluación de residuales. Presentar la evaluación de residuales del modelo seleccionado, incluyendo una prueba de normalidad y la prueba de Ljung-Box, interpretando ambos resultados.

checkresiduals(modelo_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 13.467, df = 12, p-value = 0.336
## 
## Model df: 2.   Total lags used: 14
res = residuals(modelo_arima)[-(1:13)]
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.97707, p-value = 0.3279

Interpretación:

Los residuales fluctúan alrededor de cero sin patrón visible, entre unas -550 y 850 muertes. Los 13 primeros son prácticamente cero porque corresponden al arranque de las dos diferencias; por eso se excluyeron de la prueba de normalidad, y también explican la barra alta en cero del histograma. En la ACF ningún rezago sale claramente de las bandas; el mayor, en el rezago 7, apenas las toca.

Ljung-Box: la hipótesis nula es que no hay autocorrelación en los residuales hasta el rezago 14. Con Q* = 13.467, 12 grados de libertad y p = 0.336, no se rechaza H0 al 5 %.

Shapiro-Wilk: la hipótesis nula es que los residuales siguen una distribución normal. Con W = 0.977 y p = 0.328, calculado sobre los 59 residuales restantes, no se rechaza H0 al 5 %.

Los residuales se comportan como ruido blanco con distribución aproximadamente normal, así que el modelo recoge adecuadamente la dependencia de la serie.

Tarea 20.

Pronóstico y comentario final. Generar el pronóstico del modelo seleccionado y comentar brevemente qué tan razonable resulta, considerando lo observado en la serie.

pronostico = forecast(modelo_arima, h = 12, level = c(80, 95))
pronostico
##          Point Forecast     Lo 80     Hi 80    Lo 95     Hi 95
## Jan 1979       8336.061  7924.712  8747.410 7706.957  8965.166
## Feb 1979       7531.829  7058.464  8005.194 6807.880  8255.778
## Mar 1979       8314.644  7786.496  8842.792 7506.911  9122.377
## Apr 1979       8616.869  8039.109  9194.629 7733.261  9500.477
## May 1979       9488.913  8865.476 10112.349 8535.449 10442.376
## Jun 1979       9859.757  9193.770 10525.745 8841.218 10878.297
## Jul 1979      10907.470 10201.492 11613.448 9827.770 11987.171
## Aug 1979      10086.508  9342.686 10830.331 8948.930 11224.086
## Sep 1979       9164.959  8385.127  9944.791 7972.309 10357.609
## Oct 1979       9384.259  8570.009 10198.510 8138.971 10629.548
## Nov 1979       8884.974  8037.702  9732.246 7589.183 10180.765
## Dec 1979       9376.574  8497.519 10255.628 8032.176 10720.971
plot(pronostico, include = 36, xlab = "Año", ylab = "Número de muertes")

Interpretación:

El modelo pronostica para 1979 un mínimo en febrero (7532 muertes) y un máximo en julio (10907), el mismo patrón estacional de los años anteriores. Los valores quedan, en promedio, unas 360 muertes por mes por encima de 1978, porque el modelo prolonga la recuperación de 1977-1978, y se mantienen dentro del rango histórico de la serie (6892 a 11317). Los intervalos del 95 % se ensanchan con el horizonte: miden unas 1260 muertes en enero y unas 2690 en diciembre. El pronóstico es razonable, sobre todo para los primeros meses; hacia el final del año pierde precisión, como mostró la validación, y supone que la recuperación reciente continúa.

Serie 3: fdeaths

Número mensual de muertes por enfermedades respiratorias crónicas en mujeres, en el Reino Unido, entre 1974 y 1979. Es una serie de salud pública con un patrón estacional muy marcado: las enfermedades respiratorias se agravan considerablemente en los meses fríos del invierno, generando picos de mortalidad recurrentes cada año. A diferencia de las otras dos series, su tendencia no es tan “de manual” (ni claramente creciente ni claramente decreciente), lo que la convierte en un buen caso para evaluar qué tan bien un modelo distingue la estacionalidad dominante de los matices más sutiles del comportamiento de largo plazo.

library(MASS)
data("fdeaths")
serie = fdeaths
serie
##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 1974  901  689  827  677  522  406  441  393  387  582  578  666
## 1975  830  752  785  664  467  438  421  412  343  440  531  771
## 1976  767 1141  896  532  447  420  376  330  357  445  546  764
## 1977  862  660  663  643  502  392  411  348  387  385  411  638
## 1978  796  853  737  546  530  446  431  362  387  430  425  679
## 1979  821  785  727  612  478  429  405  379  393  411  487  574

Tarea 21.

Graficar la serie. Presentar el gráfico de la serie temporal completa e identificar visualmente sus características generales (tendencia, variabilidad, posibles cambios de comportamiento).

plot(serie, type = "o", pch = 20,
     xlab = "A\u00f1o",
     ylab = "N\u00famero de muertes",
     main = "Muertes respiratorias mensuales en mujeres (Reino Unido)")

Interpretación: La serie no tiene una tendencia clara: el nivel se mantiene estable, con una leve caída. Domina una estacionalidad regular, con picos en enero o febrero y mínimos en agosto o septiembre, de amplitud similar cada año. Lo único atípico es febrero de 1976 (1141 muertes), muy por encima de los demás inviernos.

Tarea 22.

Descomposición y comentario. Aplicar decompose() para separar la serie en tendencia, estacionalidad y componente irregular, y comentar brevemente lo que muestra cada componente.

desc = decompose(serie, type = "additive")
plot(desc)

setNames(round(desc$figure), month.abb)
##  Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec 
##  253  277  200   39  -73 -131 -148 -195 -192 -107  -64  142

Interpretación: Se usó la versión aditiva porque la amplitud estacional es estable. La tendencia varía poco (entre 520 y 605 muertes) y baja levemente; su alza en 1975-1976 se debe al pico de febrero de 1976. La estacionalidad domina: febrero está 277 muertes sobre la tendencia y agosto 195 por debajo. El irregular oscila casi siempre dentro de ±100, salvo febrero de 1976 (+278) y febrero de 1977 (-165).

Tarea 23.

Pruebas de estacionariedad y su interpretación. Aplicar dos pruebas formales de estacionariedad sobre la serie original (ADF y KPSS), e interpretar sus resultados: hipótesis nula de cada una, decisión al 5% de significancia, y conclusión sobre si la serie es o no estacionaria.

adf  = adf.test(serie)
## Warning in adf.test(serie): p-value smaller than printed p-value
kpss = kpss.test(serie, null = "Level")
## Warning in kpss.test(serie, null = "Level"): p-value greater than printed
## p-value
adf
## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie
## Dickey-Fuller = -6.7742, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
kpss
## 
##  KPSS Test for Level Stationarity
## 
## data:  serie
## KPSS Level = 0.11246, Truncation lag parameter = 3, p-value = 0.1
p_num  = c(adf$p.value, kpss$p.value)

pruebas  = data.frame(Prueba = c("ADF", "KPSS"),
  H0 = c(paste0("Ra", intToUtf8(237), "z unitaria (no estacionaria)"),
         "Estacionaria en nivel"),
  Estadistico = as.numeric(c(adf$statistic, kpss$statistic)),
  p_valor = c("< 0.01", "> 0.10"),
  Decision = ifelse(p_num < 0.05, "Rechazar H0", "No rechazar H0"))

knitr::kable(pruebas, digits = 3,
  caption = "Pruebas de estacionariedad sobre la serie original")
Pruebas de estacionariedad sobre la serie original
Prueba H0 Estadistico p_valor Decision
ADF Raíz unitaria (no estacionaria) -6.774 < 0.01 Rechazar H0
KPSS Estacionaria en nivel 0.112 > 0.10 No rechazar H0

Interpretación: ADF (H0: raíz unitaria): estadístico -6.774 y p < 0.01, se rechaza H0. KPSS (H0: estacionaria en nivel): estadístico 0.112 y p > 0.10, no se rechaza H0. Ambas pruebas indican que la serie es estacionaria en nivel; la estacionalidad, que estas pruebas no evalúan, sigue presente.

Tarea 24.

Correlogramas y su interpretación. Presentar el ACF y el PACF de la serie, y comentar los patrones observados (decaimiento, cortes abruptos, etc.).

par(mfrow = c(1, 2))
acf_s = Acf(serie, lag.max = 36, main = "ACF")
pacf_s = Pacf(serie, lag.max = 36, main = "PACF")

par(mfrow = c(1, 1))
acf_s
## 
## Autocorrelations of series 'serie', by lag
## 
##      0      1      2      3      4      5      6      7      8      9     10 
##  1.000  0.730  0.381 -0.018 -0.389 -0.631 -0.696 -0.618 -0.396 -0.035  0.349 
##     11     12     13     14     15     16     17     18     19     20     21 
##  0.628  0.716  0.619  0.358  0.000 -0.295 -0.497 -0.599 -0.536 -0.377 -0.072 
##     22     23     24     25     26     27     28     29     30     31     32 
##  0.241  0.524  0.603  0.523  0.296  0.028 -0.219 -0.369 -0.438 -0.417 -0.296 
##     33     34     35     36 
## -0.081  0.161  0.349  0.420
pacf_s
## 
## Partial autocorrelations of series 'serie', by lag
## 
##      1      2      3      4      5      6      7      8      9     10     11 
##  0.730 -0.323 -0.358 -0.330 -0.231 -0.161 -0.220 -0.142  0.116  0.199  0.146 
##     12     13     14     15     16     17     18     19     20     21     22 
##  0.032  0.069  0.037 -0.069  0.101  0.116 -0.013  0.006 -0.172  0.098  0.010 
##     23     24     25     26     27     28     29     30     31     32     33 
##  0.138 -0.107 -0.013 -0.065  0.017  0.028  0.185  0.079  0.030 -0.077  0.001 
##     34     35     36 
##  0.023 -0.021 -0.051

Interpretación: La ACF forma una onda de periodo 12, con picos en los rezagos 12, 24 y 36 (0.716, 0.603 y 0.420) que decaen lentamente: la estacionalidad es persistente y requiere diferencia estacional. En los primeros rezagos no hay decaimiento lento, así que no hay tendencia fuerte. La PACF tiene un pico en el rezago 1 (0.730), valores negativos significativos en los rezagos 2 a 4, y luego se corta.

Tarea 25.

Ajuste SARIMA/ARIMA. Ajustar un modelo mediante auto.arima() y reportar el orden (p,d,q)(P,D,Q) de cada componente junto con los coeficientes estimados.

modelo_arima = auto.arima(serie, stepwise = FALSE, approximation = FALSE)
modelo_arima
## Series: serie 
## ARIMA(0,0,0)(2,1,0)[12] with drift 
## 
## Coefficients:
##          sar1     sar2    drift
##       -0.8721  -0.4950  -0.9642
## s.e.   0.1266   0.1298   0.4139
## 
## sigma^2 = 5887:  log likelihood = -349.88
## AIC=707.76   AICc=708.49   BIC=716.14
arimaorder(modelo_arima)
##         p         d         q         P         D         Q Frequency 
##         0         0         0         2         1         0        12

Interpretación: auto.arima() eligió un SARIMA(0,0,0)(2,1,0)[12] con deriva: sin términos regulares, con una diferencia estacional y dos términos AR estacionales. Coeficientes: sar1 = -0.8721 (e.e. 0.1266), sar2 = -0.4950 (e.e. 0.1298) y deriva = -0.9642 (e.e. 0.4139), los tres significativos al 5 %. La deriva implica una caída de unas 11.6 muertes por año. AIC = 707.76 y BIC = 716.14.

Tarea 26.

Ajuste de un Modelo Aditivo Generalizado (GAM). Ajustar un Modelo Aditivo Generalizado (GAM) sobre la misma serie, reportando sus parámetros o componentes principales.

library(mgcv)
datos = data.frame(y = as.numeric(serie),
                   t = as.numeric(time(serie)),
                   mes = as.numeric(cycle(serie)))
formula_gam = y ~ s(t, k = 6) + s(mes, bs = "cc", k = 12)
nudos = list(mes = c(0.5, 12.5))

modelo_gam = gam(formula_gam, data = datos, knots = nudos, method = "REML")
summary(modelo_gam)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(t, k = 6) + s(mes, bs = "cc", k = 12)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   560.68       8.07   69.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df      F p-value    
## s(t)   1.00  1.001  4.885  0.0307 *  
## s(mes) 5.76 10.000 40.160  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.855   Deviance explained = 86.9%
## -REML = 407.91  Scale est. = 4689.5    n = 72
plot(modelo_gam, pages = 1, shade = TRUE)

Interpretación: El GAM (gaussiano, estimado por REML) tiene un intercepto de 560.68, el nivel medio de la serie. La tendencia s(t) es lineal (edf = 1.00) y significativa (p = 0.031), con una caída de unas 10 muertes por año. La estacionalidad s(mes) (edf = 5.76, p < 2e-16) tiene su máximo en enero-febrero, unas 270 muertes sobre el promedio, y su mínimo en agosto, unas 185 por debajo. El modelo explica el 86.9 % de la devianza.

Tarea 27.

Comparación de métricas y selección del modelo. Construir un cuadro comparativo de ambos modelos con AIC, BIC, RMSE y MAE, señalando explícitamente cuáles de estas métricas son comparables entre las dos familias de modelo y cuáles no, y seleccionar el modelo preferido justificando la elección.

train = window(serie, end = c(1978, 12))
test = window(serie, start = c(1979, 1))

arima_train = auto.arima(train, stepwise = FALSE, approximation = FALSE)
arimaorder(arima_train)
##         p         d         q         P         D         Q Frequency 
##         0         0         0         1         1         0        12
pred_arima = forecast(arima_train, h = 12)$mean

gam_train = gam(formula_gam, data = datos[1:60, ], knots = nudos,
                method = "REML")
pred_gam = predict(gam_train, newdata = datos[61:72, ])

e_arima = as.numeric(test - pred_arima)
e_gam = as.numeric(test) - as.numeric(pred_gam)

cuadro = data.frame(Modelo = c("auto.arima", "GAM"),
  AIC = c(modelo_arima$aic, AIC(modelo_gam)),
  BIC = c(modelo_arima$bic, BIC(modelo_gam)),
  RMSE = c(sqrt(mean(e_arima^2)), sqrt(mean(e_gam^2))),
  MAE = c(mean(abs(e_arima)), mean(abs(e_gam))))
kable(cuadro, digits = 2,
      caption = "Comparación de modelos (RMSE y MAE: pronóstico de 1979)")
Comparación de modelos (RMSE y MAE: pronóstico de 1979)
Modelo AIC BIC RMSE MAE
auto.arima 707.76 716.14 36.82 27.77
GAM 822.77 843.30 36.86 30.85

Interpretación: Ambos modelos se reajustaron con 1974-1978 y pronosticaron 1979.

Comparabilidad de las métricas entre SARIMA y GAM
Métrica ¿Comparable? Motivo
AIC No El SARIMA se calcula sobre la serie diferenciada (60 obs.) y el GAM sobre 72 obs. con un ajuste penalizado.
BIC No Mismo motivo que el AIC.
RMSE Sí Mismos 12 meses de 1979 y misma unidad.
MAE Sí Mismos 12 meses de 1979 y misma unidad.

Modelo seleccionado: SARIMA(0,0,0)(2,1,0)[12] con deriva. Empata con el GAM en RMSE (36.82 frente a 36.86) y tiene menor MAE (27.77 frente a 30.85), con solo tres coeficientes.

Tarea 28.

Validación del modelo seleccionado. Evaluar el modelo seleccionado mediante una partición train-test o, preferiblemente, validación cruzada por origen móvil.

o = modelo_arima$arma  # p, q, P, Q, periodo, d, D
f_sarima = function(x, h) {
  m = Arima(x, order = o[c(1, 6, 2)], seasonal = o[c(3, 7, 4)],
            include.drift = "drift" %in% names(coef(modelo_arima)),
            method = "ML")
  forecast(m, h = h)}

e_cv = tsCV(train, f_sarima, h = 12, initial = 35)

validacion = data.frame(Horizonte = 1:12,
  RMSE = sqrt(colMeans(e_cv^2, na.rm = TRUE)),
  MAE = colMeans(abs(e_cv), na.rm = TRUE),
  Pronosticos = colSums(!is.na(e_cv)))

kable(validacion, digits = 1, row.names = FALSE,
      caption = "Validación por origen móvil del SARIMA (1977-1978)")
Validación por origen móvil del SARIMA (1977-1978)
Horizonte RMSE MAE Pronosticos
1 206.6 96.5 24
2 167.2 93.3 23
3 89.8 67.2 22
4 73.2 59.9 21
5 69.9 58.2 20
6 68.7 56.8 19
7 67.6 58.4 18
8 77.8 66.1 17
9 100.7 74.2 16
10 97.3 69.7 15
11 78.3 63.3 14
12 117.5 84.0 13

Interpretación: En 24 orígenes, con pronósticos dentro de 1977-1978, el error es estable entre 3 y 8 meses adelante (RMSE de 68 a 90) y sube a 78-118 hacia los 12 meses. Los horizontes 1 y 2 tienen RMSE altos (206.6 y 167.2) por unos pocos errores grandes, probablemente de febrero de 1977, tras el pico de 1976. En general, el modelo pronostica de forma aceptable a lo largo del año.

29. Evaluación de residuales. Presentar la evaluación de residuales del modelo seleccionado, incluyendo una prueba de normalidad y la prueba de Ljung-Box, interpretando ambos resultados.

checkresiduals(modelo_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0)(2,1,0)[12] with drift
## Q* = 6.8378, df = 12, p-value = 0.8681
## 
## Model df: 2.   Total lags used: 14
k = o[6] + 12 * o[7]
res = residuals(modelo_arima)[(k + 1):length(serie)]
shapiro.test(res)
## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.83716, p-value = 1.3e-06

Interpretación: Ljung-Box (H0: sin autocorrelación hasta el rezago 14): Q* = 6.838 y p = 0.868, no se rechaza H0, así que los residuales no están autocorrelacionados. Shapiro-Wilk (H0: normalidad), sin los 12 residuales iniciales: W = 0.837 y p = 1.3e-06, se rechaza H0. La falta de normalidad se debe sobre todo al residual atípico de febrero de 1976 (cerca de +385).

30. Pronóstico y comentario final. Generar el pronóstico del modelo seleccionado y comentar brevemente qué tan razonable resulta, considerando lo observado en la serie.

pronostico = forecast(modelo_arima, h = 12, level = c(80, 95))
pronostico
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1980       804.4767 706.1498 902.8036 654.0987 954.8547
## Feb 1980       721.3858 623.0590 819.7127 571.0079 871.7638
## Mar 1980       671.7043 573.3774 770.0312 521.3263 822.0823
## Apr 1980       575.0632 476.7363 673.3901 424.6852 725.4412
## May 1980       482.1035 383.7766 580.4304 331.7255 632.4815
## Jun 1980       389.7089 291.3821 488.0358 239.3310 540.0869
## Jul 1980       390.3876 292.0607 488.7145 240.0096 540.7656
## Aug 1980       329.8553 231.5284 428.1822 179.4773 480.2333
## Sep 1980       360.3786 262.0518 458.7055 210.0007 510.7566
## Oct 1980       377.9081 279.5812 476.2349 227.5301 528.2860
## Nov 1980       398.6089 300.2820 496.9357 248.2309 548.9868
## Dec 1980       617.8923 519.5655 716.2192 467.5143 768.2703
plot(pronostico, include = 36, xlab = "Anio", ylab = "Numero de muertes")

Interpretación: El pronóstico de 1980 conserva el patrón estacional, con el máximo en enero (804) y el mínimo en agosto (330). Queda unas 32 muertes por mes por debajo de 1979 por la leve caída de la serie, y algunos meses tocan su mínimo histórico. Es razonable a corto plazo; los intervalos del 95 % (±150 muertes) son aproximados por la falta de normalidad de los residuales.