Pronósticos y Series de Tiempo. Un enfoque aplicado usando R.

Juan Pablo Acuña González1*

1Universidad Autónoma de Guerrero

*22253567@uagro.mx

Contents

Se presentan los ejemplos numéricos sobre series de tiempo del texto de Bowerman, O’Connell, and Koehler (2005) implementados en R (R Core Team 2023), usando las paqueterías fpp3 (R. J. Hyndman and Athanasopoulos 2021; Rob Hyndman 2023) y forecast (R. J. Hyndman and Khandakar 2008; R. Hyndman et al. 2024).

1 Introducción a los pronósticos

2 Conceptos básicos de estadística

3 Regresión lineal simple

4 Regresión lineal múltiple

5 Construcción de modelos y análisi residual

6 Regresión de series temporales

6.1 Venta de bacalao

La Figura 1 muestra la gŕafica de datos de venta de bacalao.

library(fpp3)

## ── Attaching packages ────────────────────────────────────────────── fpp3 0.5 ──

## ✔ tibble      3.2.1     ✔ tsibble     1.1.4
## ✔ dplyr       1.1.4     ✔ tsibbledata 0.4.1
## ✔ tidyr       1.3.1     ✔ feasts      0.3.2
## ✔ lubridate   1.9.3     ✔ fable       0.3.4
## ✔ ggplot2     3.5.1     ✔ fabletools  0.4.2

## ── Conflicts ───────────────────────────────────────────────── fpp3_conflicts ──
## ✖ lubridate::date()    masks base::date()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ tsibble::intersect() masks base::intersect()
## ✖ tsibble::interval()  masks lubridate::interval()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ tsibble::setdiff()   masks base::setdiff()
## ✖ tsibble::union()     masks base::union()

datos_venta_bacalao <- tsibble(
    date = rep(yearmonth("2007 Jan") + 0:23),
    value =  c(362,381,317,297,399,402,375,349,386,328,389,343,276,334,394,334,384,314,344,337,345,362,314,365),
    index = date
    ) 
datos_venta_bacalao  |>
    autoplot(value) +
    labs(title = "Venta bacalao", y = "ventas") + 
  geom_hline(yintercept = mean(datos_venta_bacalao$value), linetype = "dashed", color = "red")

Figure 1: Gráfica de pesca de bacalao (toneladas) contra tiempo (meses) del Ejemplo 6.1

El modelo lineal es:

datos_venta_bacalao |>
  model(TSLM(value ~ 1)) |>
  report()

## Series: value 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -75.292 -18.792  -4.292  30.458  50.708 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  351.292      6.904   50.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.82 on 23 degrees of freedom

La predicción e intervalos del modelo con intercepto son:

datos_venta_bacalao |>
  model(TSLM(value ~ 1)) |>
    forecast(h = "3 years") |>
    autoplot(datos_venta_bacalao, level = 95) +
    labs(y = "ventas", title = "Venta bacalao") 

Figure 2: Pronósticos de pesca de bacalao (toneladas) contra tiempo (meses) del Ejemplo 6.1

datos_venta_bacalao |>
  model(TSLM(value ~ 1)) |>
  forecast(h = "3 years") |>
  hilo(level = c(95))

## # A tsibble: 36 x 5 [1M]
## # Key:       .model [1]
##    .model              date        value .mean                  `95%`
##    <chr>              <mth>       <dist> <dbl>                 <hilo>
##  1 TSLM(value ~ 1) 2009 ene N(351, 1192)  351. [283.6289, 418.9545]95
##  2 TSLM(value ~ 1) 2009 feb N(351, 1192)  351. [283.6289, 418.9545]95
##  3 TSLM(value ~ 1) 2009 mar N(351, 1192)  351. [283.6289, 418.9545]95
##  4 TSLM(value ~ 1) 2009 abr N(351, 1192)  351. [283.6289, 418.9545]95
##  5 TSLM(value ~ 1) 2009 may N(351, 1192)  351. [283.6289, 418.9545]95
##  6 TSLM(value ~ 1) 2009 jun N(351, 1192)  351. [283.6289, 418.9545]95
##  7 TSLM(value ~ 1) 2009 jul N(351, 1192)  351. [283.6289, 418.9545]95
##  8 TSLM(value ~ 1) 2009 ago N(351, 1192)  351. [283.6289, 418.9545]95
##  9 TSLM(value ~ 1) 2009 sep N(351, 1192)  351. [283.6289, 418.9545]95
## 10 TSLM(value ~ 1) 2009 oct N(351, 1192)  351. [283.6289, 418.9545]95
## # ℹ 26 more rows

6.2 Venta de calculadoras

La Figura 3 muestra la gŕafica de datos de ventas de calculadoras.

datos_venta_calculadoras <- tsibble(
  date = rep(yearmonth("2007 Jan") + 0:23),
  value =  c(197,211,203,247,239,269,308,262,258,256,261,288, 296,276,305,308,356,393,363,386,443,308,358,384),
  index = date
) 
datos_venta_calculadoras  |>
  autoplot(value) +
  labs(title = "Venta de calculadoras", y = "Ventas") + 
  geom_hline(yintercept = mean(datos_venta_calculadoras$value), linetype = "dashed", color = "red")

Figure 3: Gráfica de ventan de calculadoras en función del tiempo (meses) del Ejemplo 6.2

El modelo lineal es:

datos_venta_calculadoras  |>
  model(TSLM(value ~ trend())) |>
  report()

## Series: value 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -67.665 -19.227  -6.958  17.682  75.410 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 198.0290    13.3444  14.840 6.10e-13 ***
## trend()       8.0743     0.9339   8.646 1.59e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.67 on 22 degrees of freedom
## Multiple R-squared: 0.7726,  Adjusted R-squared: 0.7623
## F-statistic: 74.75 on 1 and 22 DF, p-value: 1.5893e-08

El ANOVA es:

lm(value ~ 1 + date, data = as.data.frame(datos_venta_calculadoras)) |>
  anova()

## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## date       1  74932   74932  74.561 1.624e-08 ***
## Residuals 22  22109    1005                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pronósticos:

datos_venta_calculadoras |>
  model(TSLM(value ~ trend())) |>
  forecast(h = "3 years") |>
  autoplot(datos_venta_calculadoras, level = 95) +
  labs(y = "ventas", title = "venta calculadoras") 

Figure 4: Pronósticos de ventan de calculadoras en función del tiempo (meses) del Ejemplo 6.2

datos_venta_calculadoras |>
  model(TSLM(value ~ trend())) |>
  forecast(h = "3 years") |>
  hilo(level = c(95)) |> 
  print(n=36)

## # A tsibble: 36 x 5 [1M]
## # Key:       .model [1]
##    .model                    date        value .mean                  `95%`
##    <chr>                    <mth>       <dist> <dbl>                 <hilo>
##  1 TSLM(value ~ trend()) 2009 ene N(400, 1181)  400. [332.5293, 467.2461]95
##  2 TSLM(value ~ trend()) 2009 feb N(408, 1204)  408. [339.9601, 475.9640]95
##  3 TSLM(value ~ trend()) 2009 mar N(416, 1228)  416. [347.3481, 484.7246]95
##  4 TSLM(value ~ trend()) 2009 abr N(424, 1254)  424. [354.6946, 493.5268]95
##  5 TSLM(value ~ trend()) 2009 may N(432, 1282)  432. [362.0010, 502.3692]95
##  6 TSLM(value ~ trend()) 2009 jun N(440, 1312)  440. [369.2684, 511.2505]95
##  7 TSLM(value ~ trend()) 2009 jul N(448, 1343)  448. [376.4982, 520.1693]95
##  8 TSLM(value ~ trend()) 2009 ago N(456, 1376)  456. [383.6918, 529.1245]95
##  9 TSLM(value ~ trend()) 2009 sep N(464, 1411)  464. [390.8504, 538.1146]95
## 10 TSLM(value ~ trend()) 2009 oct N(473, 1448)  473. [397.9753, 547.1384]95
## 11 TSLM(value ~ trend()) 2009 nov N(481, 1486)  481. [405.0677, 556.1946]95
## 12 TSLM(value ~ trend()) 2009 dic N(489, 1526)  489. [412.1291, 565.2820]95
## 13 TSLM(value ~ trend()) 2010 ene N(497, 1568)  497. [419.1604, 574.3993]95
## 14 TSLM(value ~ trend()) 2010 feb N(505, 1612)  505. [426.1630, 583.5454]95
## 15 TSLM(value ~ trend()) 2010 mar N(513, 1657)  513. [433.1380, 592.7191]95
## 16 TSLM(value ~ trend()) 2010 abr N(521, 1704)  521. [440.0865, 601.9193]95
## 17 TSLM(value ~ trend()) 2010 may N(529, 1753)  529. [447.0097, 611.1448]95
## 18 TSLM(value ~ trend()) 2010 jun N(537, 1804)  537. [453.9085, 620.3947]95
## 19 TSLM(value ~ trend()) 2010 jul N(545, 1856)  545. [460.7840, 629.6679]95
## 20 TSLM(value ~ trend()) 2010 ago N(553, 1910)  553. [467.6371, 638.9635]95
## 21 TSLM(value ~ trend()) 2010 sep N(561, 1966)  561. [474.4689, 648.2804]95
## 22 TSLM(value ~ trend()) 2010 oct N(569, 2024)  569. [481.2801, 657.6178]95
## 23 TSLM(value ~ trend()) 2010 nov N(578, 2083)  578. [488.0718, 666.9749]95
## 24 TSLM(value ~ trend()) 2010 dic N(586, 2144)  586. [494.8446, 676.3507]95
## 25 TSLM(value ~ trend()) 2011 ene N(594, 2207)  594. [501.5995, 685.7446]95
## 26 TSLM(value ~ trend()) 2011 feb N(602, 2271)  602. [508.3371, 695.1556]95
## 27 TSLM(value ~ trend()) 2011 mar N(610, 2338)  610. [515.0582, 704.5832]95
## 28 TSLM(value ~ trend()) 2011 abr N(618, 2406)  618. [521.7636, 714.0266]95
## 29 TSLM(value ~ trend()) 2011 may N(626, 2475)  626. [528.4537, 723.4851]95
## 30 TSLM(value ~ trend()) 2011 jun N(634, 2547)  634. [535.1294, 732.9581]95
## 31 TSLM(value ~ trend()) 2011 jul N(642, 2620)  642. [541.7912, 742.4451]95
## 32 TSLM(value ~ trend()) 2011 ago N(650, 2695)  650. [548.4396, 751.9453]95
## 33 TSLM(value ~ trend()) 2011 sep N(658, 2772)  658. [555.0753, 761.4583]95
## 34 TSLM(value ~ trend()) 2011 oct N(666, 2850)  666. [561.6988, 770.9836]95
## 35 TSLM(value ~ trend()) 2011 nov N(674, 2931)  674. [568.3105, 780.5205]95
## 36 TSLM(value ~ trend()) 2011 dic N(682, 3013)  682. [574.9109, 790.0688]95

6.3 Solicitudes de préstamo

La Figura 5 muestra la gŕafica de datos de solicitudes de préstamo de la universidad estatal.

datos_universidad_estatal <- tsibble(
  Year = 1:24,
  value =  c(297,249,340,406,464,481,549, 553, 556,642,670,712,808,809,867,855,965,921,956, 990,1019, 1021,1033,1127),
  index = Year
) 
datos_universidad_estatal  |>
  autoplot(value) +
  labs(title = "Solicitudes préstamo", y = "solicitudes") + 
  geom_hline(yintercept = mean(datos_universidad_estatal$value), linetype = "dashed", color = "red")

Figure 5: Solicitudes de préstamo UE del Ejemplo 6.3

El modelo lineal es:

datos_universidad_estatal |>
  model(TSLM(value ~ Year + I(Year^2))) |>
  report()

## Series: value 
## Model: TSLM 
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.063 -17.274  -5.134  21.455  63.531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 199.6196    20.8480   9.575 4.12e-09 ***
## Year         50.9366     3.8424  13.256 1.14e-11 ***
## I(Year^2)    -0.5677     0.1492  -3.805  0.00104 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.25 on 21 degrees of freedom
## Multiple R-squared: 0.9871,  Adjusted R-squared: 0.9859
## F-statistic: 802.3 on 2 and 21 DF, p-value: < 2.22e-16

El ANOVA es:

lm(value ~ 1 + Year + I(Year^2), data = as.data.frame(datos_universidad_estatal)) |>
  anova()

## Analysis of Variance Table
## 
## Response: value
##           Df  Sum Sq Mean Sq  F value    Pr(>F)    
## Year       1 1552596 1552596 1590.178 < 2.2e-16 ***
## I(Year^2)  1   14134   14134   14.477  0.001035 ** 
## Residuals 21   20504     976                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pronósticos:

datos_universidad_estatal |>
    model(TSLM(value ~ Year + I(Year^2))) |>
    refit(datos_universidad_estatal) |>
    forecast(h = 3) |>
    autoplot(datos_universidad_estatal, level = 95) +
    labs(y = "solicitudes", title = "Solicitudes préstamos") 

Figure 6: Pronóstico de solicitudes de préstamo UE del Ejemplo 6.3

datos_universidad_estatal |>
  model(TSLM(value ~ Year + I(Year^2))) |>
  refit(datos_universidad_estatal) |>
  forecast(h = 3) |>
  hilo(level = c(95))

## # A tsibble: 3 x 5 [1]
## # Key:       .model [1]
##   .model                         Year         value .mean                  `95%`
##   <chr>                         <dbl>        <dist> <dbl>                 <hilo>
## 1 TSLM(value ~ Year + I(Year^2…    25 N(1118, 1411) 1118. [1044.584, 1191.829]95
## 2 TSLM(value ~ Year + I(Year^2…    26 N(1140, 1574) 1140. [1062.441, 1217.937]95
## 3 TSLM(value ~ Year + I(Year^2…    27 N(1161, 1782) 1161. [1078.293, 1243.780]95

6.4 Residuales

Los residuales con respecto al tiempo de venta de bacalao son:

datos_venta_bacalao |>
  model(trend_model = TSLM(value ~ 1)) |>
  augment() |>
  autoplot(.innov) +
  labs(y = "ventas",
       title = "Residuales")

Figure 7: Residuales del modelo con intercepto de la venta de bacalao

Los residuales con respecto al tiempo de venta de calculadoras son:

datos_venta_calculadoras |>
  model(TSLM(value ~ trend())) |>
  augment() |>
  autoplot(.innov) +
  labs(y = "ventas",
       title = "Residuales")

Figure 8: Residuales del modelo con tendencia de la venta de calculadoras

Los residuales con respecto al tiempo (meses) de las solicitudes de préstamo son:

datos_universidad_estatal |>
  model(trend_model = TSLM(value ~ 1 + Year + I(Year^2))) |>
  augment() |>
  autoplot(.innov) +
  labs(y = "solicitudes",
       title = "Residuales")

Figure 9: Residuales del modelo exponencial de las solicitudes de préstamo

6.5 Durbin-Watson

La estadística de Durbin-Watson para el ejemplo de venta de bacalao es:

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

lm(value ~ 1, data = datos_venta_bacalao) |>
  durbinWatsonTest(simulate = TRUE)

##  lag Autocorrelation D-W Statistic p-value
##    1     -0.06797788      2.124457   0.782
##  Alternative hypothesis: rho != 0

La estadística de Durbin-Watson del ejemplo de venta de calculadoras es:

lm(value ~ 1 + date, data = datos_venta_calculadoras) |>
  durbinWatsonTest(simulate = TRUE)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.1584812      1.676269   0.244
##  Alternative hypothesis: rho != 0

La estadística de Durbin-Watson del ejemplo de solicitudes de préstamo es:

lm(value ~ 1 + Year + I(Year^2), data = datos_universidad_estatal) |>
  durbinWatsonTest(simulate = TRUE)

##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1287283           2.1    0.85
##  Alternative hypothesis: rho != 0

6.6 Promedios habitaciones

promedio_cuartos <- tsibble(
  date = rep(yearmonth("2007 Jan") + 0:167),
  value = c(501,488,504,578,545,632,728,725,585,542,480,530,518,489,528,599,572,659,739,758,602,587,497,558,555,523,532,623,598,683,774,780,609,604,531,592,578,543,565,648,615,697,785,830,645,643,551,606,585,553,576,665,656,720,826,838,652,661,584,664,623,553,599,657,680,759,878,881,705,684,577,656,645,593,617,686,679,773,906,934,713,710,600,676,645,602,601,709,706,817,930,983,745,735,620,698,665,626,649,740,729,824,937,994,781,759,643,728,691,649,656,735,748,837,995,1040,809,793,692,763,723,655,658,761,768,885,1067,1038,812,790,692,782,758,709,715,784,794,893,1046,1075,812,822,714,802,748,731,748,827,788,937,1076,1125,840,864,717,813,811,732,745,844,833,935,1110,1124,868,860,762,877), 
  index = date
) 
promedio_cuartos |>
  autoplot(value) +
  labs(title = "Promedios Habs", y = "Promedios")

Figure 10: Gráfica de los promedios de habitaciones del Ejemplo 6.6

promedio_cuartos |>
  autoplot(sqrt(value)) +
  labs(title = "Raíces Promedios Habs", y = "Raíces Promedios")

Figure 11: Raíz cuadrada de los promedios de habitaciones del Ejemplo 6.6

promedio_cuartos |>
  autoplot((value)^(1/4)) +
  labs(title = "Raíz Cuarta Promedios", y = "Raíz Cuarta Promedios")

Figure 12: Raíz cuarta de los promedios de habitaciones del Ejemplo 6.6

promedio_cuartos |>
  autoplot(log(value)) +
  labs(title = "Logaritmo Promedios", y = "Logaritmo Promedios")

Figure 13: Logaritmos de los promedios de habitaciones del Ejemplo 6.6

6.7 Promedios habitaciones estacional

promedio_cuartos <- promedio_cuartos |>
  mutate(
    Mes = month(date),
    M1 = as.numeric(Mes == 1),
    M2 = as.numeric(Mes == 2),
    M3 = as.numeric(Mes == 3),
    M4 = as.numeric(Mes == 4),
    M5 = as.numeric(Mes == 5),
    M6 = as.numeric(Mes == 6),
    M7 = as.numeric(Mes == 7),
    M8 = as.numeric(Mes == 8),
    M9 = as.numeric(Mes == 9),
    M10 = as.numeric(Mes == 10),
    M11 = as.numeric(Mes == 11)
  )
promedio_cuartos |>
  model(TSLM(log(value) ~ trend() + M1 + M2 + M3 + M4 + M5 + M6 + M7 + M8 + M9 + M10 + M11)) |> 
  report()

## Series: value 
## Model: TSLM 
## Transformation: log(value) 
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0593245 -0.0138500  0.0008501  0.0138988  0.0500432 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.290e+00  6.528e-03 963.555  < 2e-16 ***
## trend()      2.722e-03  3.432e-05  79.321  < 2e-16 ***
## M1          -4.382e-02  8.142e-03  -5.383 2.67e-07 ***
## M2          -1.143e-01  8.140e-03 -14.041  < 2e-16 ***
## M3          -8.667e-02  8.139e-03 -10.649  < 2e-16 ***
## M4           3.726e-02  8.138e-03   4.579 9.54e-06 ***
## M5           1.819e-02  8.137e-03   2.235   0.0268 *  
## M6           1.447e-01  8.136e-03  17.787  < 2e-16 ***
## M7           2.868e-01  8.135e-03  35.258  < 2e-16 ***
## M8           3.090e-01  8.134e-03  37.987  < 2e-16 ***
## M9           5.379e-02  8.134e-03   6.614 5.74e-10 ***
## M10          3.735e-02  8.133e-03   4.593 9.00e-06 ***
## M11         -1.144e-01  8.133e-03 -14.066  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02152 on 155 degrees of freedom
## Multiple R-squared: 0.9883,  Adjusted R-squared: 0.9874
## F-statistic:  1091 on 12 and 155 DF, p-value: < 2.22e-16

lm(log(value) ~ 1 + date + M1 + M2 + M3 + M4 + M5 + M6 + M7 + M8 + M9 + M10 + M11, data = as.data.frame(promedio_cuartos)) |>
  anova()

## Analysis of Variance Table
## 
## Response: log(value)
##            Df  Sum Sq Mean Sq  F value    Pr(>F)    
## date        1 3.02103 3.02103 6525.493 < 2.2e-16 ***
## M1          1 0.11745 0.11745  253.700 < 2.2e-16 ***
## M2          1 0.42614 0.42614  920.471 < 2.2e-16 ***
## M3          1 0.37521 0.37521  810.473 < 2.2e-16 ***
## M4          1 0.03741 0.03741   80.797 8.175e-16 ***
## M5          1 0.08726 0.08726  188.490 < 2.2e-16 ***
## M6          1 0.02896 0.02896   62.554 4.564e-13 ***
## M7          1 0.61522 0.61522 1328.894 < 2.2e-16 ***
## M8          1 1.10995 1.10995 2397.519 < 2.2e-16 ***
## M9          1 0.06628 0.06628  143.169 < 2.2e-16 ***
## M10         1 0.08350 0.08350  180.363 < 2.2e-16 ***
## M11         1 0.09168 0.09168  198.025 < 2.2e-16 ***
## Residuals 155 0.07176 0.00046                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lm(log(value) ~ date + M1 + M2 + M3 + M4 + M5 + M6 + M7 + M8 + M9 + M10 + M11, data = as.data.frame(promedio_cuartos)) |>
  durbinWatsonTest(simulate = TRUE)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.4015944      1.170453       0
##  Alternative hypothesis: rho != 0

promedio_cuartos |>
  model(TSLM(log(value) ~ trend() + season())) |>
  forecast(h = "2 years") |>
  autoplot(promedio_cuartos) +
  labs(y = "promedios", title = "Promedios habitaciones")

Figure 14: Pronósticos con tendencia y estacionalidad

promedio_cuartos |>
  model(TSLM(value ~ trend() + season())) |>
  forecast(h = "2 years") |>
  hilo(level = c(97.5))

## # A tsibble: 24 x 5 [1M]
## # Key:       .model [1]
##    .model                     date        value .mean                    `97.5%`
##    <chr>                     <mth>       <dist> <dbl>                     <hilo>
##  1 TSLM(value ~ trend() … 2021 ene  N(822, 510)  822. [ 771.0949,  872.3721]97.5
##  2 TSLM(value ~ trend() … 2021 feb  N(779, 510)  779. [ 728.2378,  829.5150]97.5
##  3 TSLM(value ~ trend() … 2021 mar  N(797, 510)  797. [ 745.8806,  847.1578]97.5
##  4 TSLM(value ~ trend() … 2021 abr  N(880, 510)  880. [ 828.9521,  930.2293]97.5
##  5 TSLM(value ~ trend() … 2021 may  N(869, 510)  869. [ 818.5949,  919.8721]97.5
##  6 TSLM(value ~ trend() … 2021 jun  N(965, 510)  965. [ 914.3092, 1015.5864]97.5
##  7 TSLM(value ~ trend() … 2021 jul N(1090, 510) 1090. [1039.0235, 1140.3007]97.5
##  8 TSLM(value ~ trend() … 2021 ago N(1113, 510) 1113. [1062.4521, 1163.7293]97.5
##  9 TSLM(value ~ trend() … 2021 sep  N(903, 510)  903. [ 851.9521,  953.2293]97.5
## 10 TSLM(value ~ trend() … 2021 oct  N(894, 510)  894. [ 843.0949,  944.3721]97.5
## # ℹ 14 more rows

6.8 Promedios habitaciones trigonométrico

promedio_cuartos |>
  model(TSLM(log(value)~ trend() + fourier(K = 2))) |> 
  report()

## Series: value 
## Model: TSLM 
## Transformation: log(value) 
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.100542 -0.045572 -0.002805  0.046947  0.122897 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          6.333e+00  8.740e-03 724.597  < 2e-16 ***
## trend()              2.733e-03  8.975e-05  30.445  < 2e-16 ***
## fourier(K = 2)C1_12 -1.598e-01  6.144e-03 -26.015  < 2e-16 ***
## fourier(K = 2)S1_12 -2.433e-02  6.152e-03  -3.954 0.000115 ***
## fourier(K = 2)C2_12  6.712e-02  6.144e-03  10.924  < 2e-16 ***
## fourier(K = 2)S2_12  1.638e-02  6.145e-03   2.665 0.008483 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0563 on 162 degrees of freedom
## Multiple R-squared: 0.9162,  Adjusted R-squared: 0.9137
## F-statistic: 354.4 on 5 and 162 DF, p-value: < 2.22e-16

promedio_cuartos |>
  model(TSLM(log(value) ~ trend() + fourier(K = 2))) |>
  augment() |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "Data")) +
  geom_line(aes(y = .fitted, colour = "Fitted"))  +
  labs(y = "promedios", title = "Promedios habitaciones")

Figure 15: Pronósticos con estacioonalidad trigonométrica

promedio_cuartos |>
  model(TSLM(log(value)~ trend() + fourier(K = 2))) |>
  forecast(h = "3 years") |>
  autoplot(promedio_cuartos) +
  labs(y = "promedios", title = "Promedios habitaciones")  

Figure 16: Valores ajustados vs serie temporal

promedio_cuartos |>
  model(TSLM(value ~ trend() + season())) |>
  forecast(h = "3 years") |>
  hilo(level = c(97.5))

## # A tsibble: 36 x 5 [1M]
## # Key:       .model [1]
##    .model                     date        value .mean                    `97.5%`
##    <chr>                     <mth>       <dist> <dbl>                     <hilo>
##  1 TSLM(value ~ trend() … 2021 ene  N(822, 510)  822. [ 771.0949,  872.3721]97.5
##  2 TSLM(value ~ trend() … 2021 feb  N(779, 510)  779. [ 728.2378,  829.5150]97.5
##  3 TSLM(value ~ trend() … 2021 mar  N(797, 510)  797. [ 745.8806,  847.1578]97.5
##  4 TSLM(value ~ trend() … 2021 abr  N(880, 510)  880. [ 828.9521,  930.2293]97.5
##  5 TSLM(value ~ trend() … 2021 may  N(869, 510)  869. [ 818.5949,  919.8721]97.5
##  6 TSLM(value ~ trend() … 2021 jun  N(965, 510)  965. [ 914.3092, 1015.5864]97.5
##  7 TSLM(value ~ trend() … 2021 jul N(1090, 510) 1090. [1039.0235, 1140.3007]97.5
##  8 TSLM(value ~ trend() … 2021 ago N(1113, 510) 1113. [1062.4521, 1163.7293]97.5
##  9 TSLM(value ~ trend() … 2021 sep  N(903, 510)  903. [ 851.9521,  953.2293]97.5
## 10 TSLM(value ~ trend() … 2021 oct  N(894, 510)  894. [ 843.0949,  944.3721]97.5
## # ℹ 26 more rows

A continuación se muestra el código en base R para realizar estos cálculos:

# Datos
t <- 1:nrow(promedio_cuartos)
sin_2pi_t_12 <- sin(2 * pi * t / 12)
cos_2pi_t_12 <- cos(2 * pi * t / 12)
sin_4pi_t_12 <- sin(4 * pi * t / 12)
cos_4pi_t_12 <- cos(4 * pi * t / 12)

# Matriz de diseño X
X <- cbind(1, t, sin_2pi_t_12, cos_2pi_t_12, sin_4pi_t_12, cos_4pi_t_12)

# Variable dependiente (logaritmo de promedio de cuartos ocupados)
ln_y <- log(promedio_cuartos$value)

# Ajustar regresión lineal
modelo_manual <- lm(ln_y ~ X - 1)  # -1 para evitar el intercepto adicional

# Coeficientes del modelo manual
coeficientes_manual <- coef(modelo_manual)
print(coeficientes_manual)

##             X            Xt Xsin_2pi_t_12 Xcos_2pi_t_12 Xsin_4pi_t_12 
##   6.333133882   0.002732605  -0.100986973  -0.126257621   0.066312449 
## Xcos_4pi_t_12 
##   0.019375927

# Predicción puntual y16
t_16 <- 169
ln_y_16 <- coeficientes_manual[1] +
            coeficientes_manual[2] * t_16 +
            coeficientes_manual[3] * sin(2 * pi * t_16 / 12) +
            coeficientes_manual[4] * cos(2 * pi * t_16 / 12) +
            coeficientes_manual[5] * sin(4 * pi * t_16 / 12) +
            coeficientes_manual[6] * cos(4 * pi * t_16 / 12)
y_16 <- exp(ln_y_16)
print(y_16)

##        X 
## 814.2151

# Intervalo de predicción al 95% para y16
# Error estándar residual
sigma <- summary(modelo_manual)$sigma
# Cuantil t para alpha = 0.05 y n - p grados de libertad
t_quantil <- qt(0.975, nrow(promedio_cuartos) - ncol(X))
# Intervalo de predicción
intervalo_prediccion <- exp(ln_y_16 + c(-1, 1) * t_quantil * sigma)
print(intervalo_prediccion)

## [1] 728.5370 909.9692

6.9 Cadena de restaurantes

western <- tsibble(
  date = rep(yearmonth("2007 Jan") + 0:14),
  value = c(11,14,16,22,28,36,46,67,82,99,119,156,257,284,403),
  index = date
) 
modelo <- western |>
  model(TSLM(log(value) ~ trend()))

# Reporte del modelo
modelo |>
  report()

## Series: value 
## Model: TSLM 
## Transformation: log(value) 
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11668 -0.04172 -0.01747  0.06304  0.13951 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.070120   0.041032   50.45 2.67e-16 ***
## trend()     0.256880   0.004513   56.92  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07552 on 13 degrees of freedom
## Multiple R-squared: 0.996,   Adjusted R-squared: 0.9957
## F-statistic:  3240 on 1 and 13 DF, p-value: < 2.22e-16

lm(log(value) ~ date, data = as.data.frame(western)) |>
  anova()

## Analysis of Variance Table
## 
## Response: log(value)
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## date       1 18.4802 18.4802  3410.2 < 2.2e-16 ***
## Residuals 13  0.0704  0.0054                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lm(log(value) ~ date, data = as.data.frame(western)) |>
  durbinWatsonTest(simulate = TRUE)

##  lag Autocorrelation D-W Statistic p-value
##    1     -0.03401941       1.92181   0.668
##  Alternative hypothesis: rho != 0

western |>
  model(TSLM(log(value) ~ trend())) |>
  forecast(h = "2 years") |>
  autoplot(western) +
  labs(y = "restaurantes", title = "Cadena de restaurantes") 

western |>
  model(TSLM(log(value) ~ trend())) |>
  forecast(h = "2 years") |>
  hilo(level = c(97.5))

## # A tsibble: 24 x 5 [1M]
## # Key:       .model [1]
##    .model                date             value .mean                    `97.5%`
##    <chr>                <mth>            <dist> <dbl>                     <hilo>
##  1 TSLM(log(value) … 2008 abr t(N(6.2, 0.0074))  485. [ 398.4460,  585.7196]97.5
##  2 TSLM(log(value) … 2008 may t(N(6.4, 0.0077))  627. [ 512.8529,  760.6586]97.5
##  3 TSLM(log(value) … 2008 jun t(N(6.7, 0.0081))  811. [ 659.8404,  988.2510]97.5
##  4 TSLM(log(value) … 2008 jul   t(N(7, 0.0085)) 1048. [ 848.6334, 1284.4273]97.5
##  5 TSLM(log(value) … 2008 ago  t(N(7.2, 0.009)) 1356. [1091.0600, 1669.9540]97.5
##  6 TSLM(log(value) … 2008 sep t(N(7.5, 0.0095)) 1753. [1402.2849, 2171.9028]97.5
##  7 TSLM(log(value) … 2008 oct   t(N(7.7, 0.01)) 2268. [1801.7482, 2825.5698]97.5
##  8 TSLM(log(value) … 2008 nov    t(N(8, 0.011)) 2933. [2314.3691, 3676.9780]97.5
##  9 TSLM(log(value) … 2008 dic  t(N(8.2, 0.011)) 3793. [2972.0875, 4786.1418]97.5
## 10 TSLM(log(value) … 2009 ene  t(N(8.5, 0.012)) 4905. [3815.8391, 6231.3277]97.5
## # ℹ 14 more rows

6.10 Promedios habitaciones autocorrelación

# Ajustar un modelo ARIMA a los residuos para capturar la autocorrelación
lm(log(value) ~ 1 + date + M1 + M2 + M3 + M4 + M5 + M6 + M7 + M8 + M9 + M10 + M11, data = as.data.frame(promedio_cuartos)) |> 
  resid() |> 
  arima(order = c(1, 0, 0)) |> 
  coef()

##           ar1     intercept 
##  4.099543e-01 -1.333176e-05

7 Métodos de descomposición

7.1 Venta de Tasty cola

La Figura 17 muestra la gŕafica de datos de venta de Tasty cola.

datos_venta_cola <- tsibble(
    date = rep(yearmonth("2007 Jan") + 0:35),
    value =  c(189,229,249,289,260,431,660,777,915,613,485,277,244,296,319,370,313,556,831,960,1152,759,607,371,298,378,373,443,374,660,1004,1153,1388,904,715,441),
    index = date
    ) 
datos_venta_cola  |>
    autoplot(value) +
    labs(title = "Venta cola", y = "ventas") + 
  geom_hline(yintercept = mean(datos_venta_cola$value), linetype = "dashed", color = "red")  +
  labs(y = "ventas", title = "Ventas Tasty cola")

Figure 17: Ventas de Tasty cola contra tiempo (meses) del Ejemplo 7.1

Los componentes son:

datos_venta_cola |>
   model(
    classical_decomposition(value, type = "multiplicative")
  )  |>
  components()

## # A dable: 36 x 7 [1M]
## # Key:     .model [1]
## # :        value = trend * seasonal * random
##    .model                         date value trend seasonal random season_adjust
##    <chr>                         <mth> <dbl> <dbl>    <dbl>  <dbl>         <dbl>
##  1 "classical_decomposition(… 2007 ene   189   NA     0.493 NA              383.
##  2 "classical_decomposition(… 2007 feb   229   NA     0.596 NA              384.
##  3 "classical_decomposition(… 2007 mar   249   NA     0.595 NA              418.
##  4 "classical_decomposition(… 2007 abr   289   NA     0.680 NA              425.
##  5 "classical_decomposition(… 2007 may   260   NA     0.564 NA              461.
##  6 "classical_decomposition(… 2007 jun   431   NA     0.986 NA              437.
##  7 "classical_decomposition(… 2007 jul   660  450.    1.47   0.999          450.
##  8 "classical_decomposition(… 2007 ago   777  455.    1.69   1.01           459.
##  9 "classical_decomposition(… 2007 sep   915  461.    1.99   0.998          460.
## 10 "classical_decomposition(… 2007 oct   613  467.    1.31   1.00           469.
## # ℹ 26 more rows

La visualización de la tendencia es:

datos_venta_cola |>
  model(
    classical_decomposition(value, type = "multiplicative")
  ) |>
  components() |>
  as_tsibble() |>
  autoplot(value, colour="gray") +
  geom_line(aes(y=trend), colour = "#D55E00") +
  labs(
    y = "Ventas (mensuales)",
    title = "Ventas mensuales Tasty cola"
  )

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

Figure 18: VIsualización de la tendencia Ventas de Tasty cola

La visualización de los componentes:

datos_venta_cola |>
  model(
    classical_decomposition(value, type = "multiplicative")
  ) |>
  components() |>
  autoplot()

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).

Figure 19: Descomposición de la serie temporal de ventas de Tasty cola

Los datos ajustados a la estacionalidad:

datos_venta_cola |>
  model(
    classical_decomposition(value, type = "multiplicative")
  ) |>
  components() |>
  as_tsibble() |>
  autoplot(value, colour = "gray") +
  geom_line(aes(y=season_adjust), colour = "#0072B2") +
  labs(y = "Ventas (meses)",
       title = "Ventas mensuales Tasty cola")

Figure 20: Serie temporal de ventas de Tasty cola ajustada a la estacionalidad

Las medias móviles son:

datos_venta_cola |>
  mutate(
    `12-MA` = slider::slide_dbl(value, mean,
                .before = 5, .after = 6, .complete = TRUE)
  ) |>
  autoplot(value) +
  geom_line(aes(y = `12-MA`), colour = "#D55E00") +
  labs(y = "% ventas",
       title = "Ventas Tasty cola")

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_line()`).

Figure 21: Medias móviles ventas de Tasty cola

Las medias móviles de las medias móviles:

datos_venta_cola |>
  mutate(
    `12-MA` = slider::slide_dbl(value, mean,
                .before = 5, .after = 6, .complete = TRUE),
    `2x12-MA` = slider::slide_dbl(`12-MA`, mean,
                .before = 1, .after = 0, .complete = TRUE)
  ) |>
  autoplot(value, colour = "gray") +
  geom_line(aes(y = `2x12-MA`), colour = "#D55E00") +
  labs(y = "Ventas (mensuales)",
       title = "Ventas Tasty cola")

## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_line()`).

Figure 22: Medias móviles de medias móviles anuales ventas de Tasty cola

Las predicciones son:

datos_venta_cola |>
  model(TSLM(value ~ I(trend()) * I(season()))) |>
  augment()

## # A tsibble: 36 x 6 [1M]
## # Key:       .model [1]
##    .model                                     date value .fitted .resid .innov
##    <chr>                                     <mth> <dbl>   <dbl>  <dbl>  <dbl>
##  1 TSLM(value ~ I(trend()) * I(season())) 2007 ene   189    189. -0.167 -0.167
##  2 TSLM(value ~ I(trend()) * I(season())) 2007 feb   229    226.  2.5    2.50 
##  3 TSLM(value ~ I(trend()) * I(season())) 2007 mar   249    252. -2.67  -2.67 
##  4 TSLM(value ~ I(trend()) * I(season())) 2007 abr   289    290. -1.33  -1.33 
##  5 TSLM(value ~ I(trend()) * I(season())) 2007 may   260    259.  1.33   1.33 
##  6 TSLM(value ~ I(trend()) * I(season())) 2007 jun   431    434. -3.5   -3.50 
##  7 TSLM(value ~ I(trend()) * I(season())) 2007 jul   660    660.  0.333  0.333
##  8 TSLM(value ~ I(trend()) * I(season())) 2007 ago   777    775.  1.67   1.67 
##  9 TSLM(value ~ I(trend()) * I(season())) 2007 sep   915    915. -0.167 -0.167
## 10 TSLM(value ~ I(trend()) * I(season())) 2007 oct   613    613. -0.167 -0.167
## # ℹ 26 more rows

La visualización de datos vs ajuste:

datos_venta_cola |>
  model(TSLM(value ~ I(trend()) * I(season()))) |>
  augment() |>
  ggplot(aes(x = date)) +
  geom_line(aes(y = value, colour = "Data")) +
  geom_line(aes(y = .fitted, colour = "Fitted")) +
  labs(y = "Ventas (mensuales)",
       title = "Ventas Tasty cola")

Figure 23: Predicciones de las Ventas de Tasty cola

8 Suavización exponencial

8.1 Suavización exponencial simple

Los datos de la Figura 1 de la venta de bacalao muestran la siguiente suavización exponencial simple:

fit <- datos_venta_bacalao |>
  model(ETS(value ~ error("M") + trend("N") + season("N")))
fc <- fit |>
  forecast(h = 5)
fc

## # A fable: 5 x 4 [1M]
## # Key:     .model [1]
##   .model                                                 date        value .mean
##   <chr>                                                 <mth>       <dist> <dbl>
## 1 "ETS(value ~ error(\"M\") + trend(\"N\") + season… 2009 ene N(351, 1196)  351.
## 2 "ETS(value ~ error(\"M\") + trend(\"N\") + season… 2009 feb N(351, 1196)  351.
## 3 "ETS(value ~ error(\"M\") + trend(\"N\") + season… 2009 mar N(351, 1196)  351.
## 4 "ETS(value ~ error(\"M\") + trend(\"N\") + season… 2009 abr N(351, 1196)  351.
## 5 "ETS(value ~ error(\"M\") + trend(\"N\") + season… 2009 may N(351, 1196)  351.

8.2 Visualización suavización simple

La visualización se muestra en la siguiente gráfica:

fc |>
  autoplot(datos_venta_bacalao) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(fit)) +
  labs(y="% ventas", title="Ventas bacalao") +
  guides(colour = "none")

Figure 24: Venta de bacalao

8.3 Método de Holt

Datos de ventas de termostatos semanales por suavización exponencial con tendencia:

datos_venta_termos <- tsibble(
    date = rep(yearweek("2007 W1") + 0:51),
    value =  c(206,245,185,169,162,177,207,216,193,230,212,192,162,189,244,209,207,211,210,173,194,234,156,206,188,162,172,210,205,244,218,182,206,211,273,248,262,258,233,255,303,282,291,280,255,312,296,307,281,308,280,345),
    index = date
    ) 
git <- datos_venta_termos |>
  model(ETS(value ~ error("M") + trend("A") + season("N")))
gc <- git |>
  forecast(h = 5)
gc

## # A fable: 5 x 4 [1W]
## # Key:     .model [1]
##   .model                                                 date        value .mean
##   <chr>                                                <week>       <dist> <dbl>
## 1 "ETS(value ~ error(\"M\") + trend(\"A\") + season… 2008 W01 N(311, 1668)  311.
## 2 "ETS(value ~ error(\"M\") + trend(\"A\") + season… 2008 W02 N(312, 1832)  312.
## 3 "ETS(value ~ error(\"M\") + trend(\"A\") + season… 2008 W03 N(313, 1998)  313.
## 4 "ETS(value ~ error(\"M\") + trend(\"A\") + season… 2008 W04 N(314, 2165)  314.
## 5 "ETS(value ~ error(\"M\") + trend(\"A\") + season… 2008 W05 N(316, 2334)  316.

gc |>
  autoplot(datos_venta_termos) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(git)) +
  labs(y="% ventas", title="Ventas termostatos") +
  guides(colour = "none")

Figure 25: Venta de termostatos

8.4 Holt-Winters aditivo

Venta de bicicletas:

datos_venta_bici <- tsibble(
    date = rep(yearquarter("2007 Q1") + 0:15),
    value =  c(10,31,43,16,11,33,45,17,13,34,48,19,15,37,51,21),
    index = date
    ) 
hit <- datos_venta_bici |>
  model(ETS(value ~ error("A") + trend("A") + season("A")))
hc <- hit |>
  forecast(h = 5)
hc

## # A fable: 5 x 4 [1Q]
## # Key:     .model [1]
##   .model                                                  date       value .mean
##   <chr>                                                  <qtr>      <dist> <dbl>
## 1 "ETS(value ~ error(\"A\") + trend(\"A\") + season(\… 2011 Q1 N(18, 0.73)  17.8
## 2 "ETS(value ~ error(\"A\") + trend(\"A\") + season(\… 2011 Q2 N(40, 0.78)  39.6
## 3 "ETS(value ~ error(\"A\") + trend(\"A\") + season(\… 2011 Q3 N(53, 0.89)  52.9
## 4 "ETS(value ~ error(\"A\") + trend(\"A\") + season(\… 2011 Q4  N(24, 1.1)  24.5
## 5 "ETS(value ~ error(\"A\") + trend(\"A\") + season(\… 2012 Q1  N(21, 1.4)  20.7

hc |>
  autoplot(datos_venta_bici) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(hit)) +
  labs(y="% ventas", title="Ventas bicicletas") +
  guides(colour = "none")

Figure 26: Venta de bicicletas contra tiempo (trimestral) del Ejemplo 8.4

8.5 Holt-Winters multiplicativo

Bebidas deportivas:

datos_venta_tiger <- tsibble(
    date = rep(yearquarter("2007 Q1") + 0:31),
    value =  c(72,116,136,96,77,123,146,101,81,131,158,109,67,140,167,120,94,147,177,128,102,162,191,134,106,170,200,142,115,177,218,149),
    index = date
    ) 
jit <- datos_venta_tiger |>
  model(ETS(value ~ error("M") + trend("A") + season("M")))
jc <- jit |>
  forecast(h = 5)
jc

## # A fable: 5 x 4 [1Q]
## # Key:     .model [1]
##   .model                                                  date       value .mean
##   <chr>                                                  <qtr>      <dist> <dbl>
## 1 "ETS(value ~ error(\"M\") + trend(\"A\") + season(\… 2015 Q1  N(115, 38)  115.
## 2 "ETS(value ~ error(\"M\") + trend(\"A\") + season(\… 2015 Q2 N(186, 100)  186.
## 3 "ETS(value ~ error(\"M\") + trend(\"A\") + season(\… 2015 Q3 N(222, 143)  222.
## 4 "ETS(value ~ error(\"M\") + trend(\"A\") + season(\… 2015 Q4  N(155, 70)  155.
## 5 "ETS(value ~ error(\"M\") + trend(\"A\") + season(\… 2016 Q1  N(120, 42)  120.

jc |>
  autoplot(datos_venta_tiger) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(jit)) +
  labs(y="% ventas", title="Ventas bebidas") +
  guides(colour = "none")

Figure 27: Venta de bebidas (miles de envases) contra tiempo (trimestral) del Ejemplo 8.5

Tasty cola:

kit <- datos_venta_cola |>
  model(ETS(value ~ error("M") + trend("A") + season("M")))
kc <- kit |>
  forecast(h = 15)
kc

## # A fable: 15 x 4 [1M]
## # Key:     .model [1]
##    .model                                               date         value .mean
##    <chr>                                               <mth>        <dist> <dbl>
##  1 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 ene    N(355, 67)  355.
##  2 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 feb   N(433, 100)  433.
##  3 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 mar   N(452, 110)  452.
##  4 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 abr   N(525, 148)  525.
##  5 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 may   N(453, 110)  453.
##  6 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 jun   N(778, 325)  778.
##  7 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 jul  N(1170, 734) 1170.
##  8 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 ago  N(1356, 987) 1356.
##  9 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 sep N(1606, 1384) 1606.
## 10 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 oct  N(1057, 600) 1057.
## 11 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 nov   N(830, 370)  830.
## 12 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2010 dic   N(496, 132)  496.
## 13 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2011 ene    N(410, 91)  410.
## 14 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2011 feb   N(500, 134)  500.
## 15 "ETS(value ~ error(\"M\") + trend(\"A\") + seas… 2011 mar   N(521, 146)  521.

kc |>
  autoplot(datos_venta_cola) +
  geom_line(aes(y = .fitted), col="#D55E00",
            data = augment(kit)) +
  labs(y="% ventas", title="Ventas Tasty cola") +
  guides(colour = "none")

Figure 28: Venta de Tasty cola

8.6 Intervalos de predicción

Modelo simple venta termostatos:

gc |>
  hilo(level = c(95))

## # A tsibble: 5 x 5 [1W]
## # Key:       .model [1]
##   .model                          date        value .mean                  `95%`
##   <chr>                         <week>       <dist> <dbl>                 <hilo>
## 1 "ETS(value ~ error(\"M\") … 2008 W01 N(311, 1668)  311. [230.7041, 390.8223]95
## 2 "ETS(value ~ error(\"M\") … 2008 W02 N(312, 1832)  312. [228.0873, 395.8856]95
## 3 "ETS(value ~ error(\"M\") … 2008 W03 N(313, 1998)  313. [225.6041, 400.8152]95
## 4 "ETS(value ~ error(\"M\") … 2008 W04 N(314, 2165)  314. [223.2383, 405.6276]95
## 5 "ETS(value ~ error(\"M\") … 2008 W05 N(316, 2334)  316. [220.9762, 410.3361]95

Modelo con tendencia venta bicicletas:

hc |>
  hilo(level = c(95))

## # A tsibble: 5 x 5 [1Q]
## # Key:       .model [1]
##   .model                           date       value .mean                  `95%`
##   <chr>                           <qtr>      <dist> <dbl>                 <hilo>
## 1 "ETS(value ~ error(\"A\") + … 2011 Q1 N(18, 0.73)  17.8 [16.17050, 19.51339]95
## 2 "ETS(value ~ error(\"A\") + … 2011 Q2 N(40, 0.78)  39.6 [37.88139, 41.33974]95
## 3 "ETS(value ~ error(\"A\") + … 2011 Q3 N(53, 0.89)  52.9 [50.99881, 54.70379]95
## 4 "ETS(value ~ error(\"A\") + … 2011 Q4  N(24, 1.1)  24.5 [22.42861, 26.53566]95
## 5 "ETS(value ~ error(\"A\") + … 2012 Q1  N(21, 1.4)  20.7 [18.41347, 23.08005]95

Método Holt-Winters aditivo venta bebidas:

jc |>
  hilo(level = c(95))

## # A tsibble: 5 x 5 [1Q]
## # Key:       .model [1]
##   .model                           date       value .mean                  `95%`
##   <chr>                           <qtr>      <dist> <dbl>                 <hilo>
## 1 "ETS(value ~ error(\"M\") + … 2015 Q1  N(115, 38)  115. [102.5602, 126.7557]95
## 2 "ETS(value ~ error(\"M\") + … 2015 Q2 N(186, 100)  186. [166.2985, 205.5761]95
## 3 "ETS(value ~ error(\"M\") + … 2015 Q3 N(222, 143)  222. [198.5267, 245.4693]95
## 4 "ETS(value ~ error(\"M\") + … 2015 Q4  N(155, 70)  155. [138.5937, 171.4012]95
## 5 "ETS(value ~ error(\"M\") + … 2016 Q1  N(120, 42)  120. [107.3787, 132.8249]95

Método Holt-Winters multiplicativo venta Tasty cola:

kc |>
  hilo(level = c(95))

## # A tsibble: 15 x 5 [1M]
## # Key:       .model [1]
##    .model                      date         value .mean                    `95%`
##    <chr>                      <mth>        <dist> <dbl>                   <hilo>
##  1 "ETS(value ~ error(\"M… 2010 ene    N(355, 67)  355. [ 338.8269,  371.0261]95
##  2 "ETS(value ~ error(\"M… 2010 feb   N(433, 100)  433. [ 413.2184,  452.4944]95
##  3 "ETS(value ~ error(\"M… 2010 mar   N(452, 110)  452. [ 431.9552,  473.0199]95
##  4 "ETS(value ~ error(\"M… 2010 abr   N(525, 148)  525. [ 501.2895,  548.9543]95
##  5 "ETS(value ~ error(\"M… 2010 may   N(453, 110)  453. [ 432.3058,  473.4187]95
##  6 "ETS(value ~ error(\"M… 2010 jun   N(778, 325)  778. [ 743.0303,  813.7062]95
##  7 "ETS(value ~ error(\"M… 2010 jul  N(1170, 734) 1170. [1116.6585, 1222.8920]95
##  8 "ETS(value ~ error(\"M… 2010 ago  N(1356, 987) 1356. [1294.1070, 1417.2436]95
##  9 "ETS(value ~ error(\"M… 2010 sep N(1606, 1384) 1606. [1532.5763, 1678.4287]95
## 10 "ETS(value ~ error(\"M… 2010 oct  N(1057, 600) 1057. [1008.6760, 1104.6861]95
## 11 "ETS(value ~ error(\"M… 2010 nov   N(830, 370)  830. [ 792.3114,  867.7398]95
## 12 "ETS(value ~ error(\"M… 2010 dic   N(496, 132)  496. [ 473.3096,  518.3764]95
## 13 "ETS(value ~ error(\"M… 2011 ene    N(410, 91)  410. [ 391.6537,  428.9516]95
## 14 "ETS(value ~ error(\"M… 2011 feb   N(500, 134)  500. [ 476.8167,  522.2323]95
## 15 "ETS(value ~ error(\"M… 2011 mar   N(521, 146)  521. [ 497.5949,  544.9973]95

8.7 Simulación intervalos de predicción

Mil realizaciones por bootstrap de las trayectorias de los intervalos de predicción para los siguientes 10 periodos de la venta de bacalao:

fit |>
  forecast(h = 10, bootstrap = TRUE, times = 1000) |>
  hilo(level = c(95))

## # A tsibble: 10 x 5 [1M]
## # Key:       .model [1]
##    .model                         date        value .mean                  `95%`
##    <chr>                         <mth>       <dist> <dbl>                 <hilo>
##  1 "ETS(value ~ error(\"M\")… 2009 ene sample[1000]  349. [275.8961, 401.9051]95
##  2 "ETS(value ~ error(\"M\")… 2009 feb sample[1000]  350. [275.8959, 401.9024]95
##  3 "ETS(value ~ error(\"M\")… 2009 mar sample[1000]  351. [275.8971, 401.8938]95
##  4 "ETS(value ~ error(\"M\")… 2009 abr sample[1000]  352. [275.8974, 401.9055]95
##  5 "ETS(value ~ error(\"M\")… 2009 may sample[1000]  352. [275.8962, 401.9000]95
##  6 "ETS(value ~ error(\"M\")… 2009 jun sample[1000]  352. [275.9003, 401.9056]95
##  7 "ETS(value ~ error(\"M\")… 2009 jul sample[1000]  352. [275.8974, 401.9020]95
##  8 "ETS(value ~ error(\"M\")… 2009 ago sample[1000]  351. [275.9002, 401.9058]95
##  9 "ETS(value ~ error(\"M\")… 2009 sep sample[1000]  351. [275.9004, 401.8996]95
## 10 "ETS(value ~ error(\"M\")… 2009 oct sample[1000]  350. [275.8967, 401.9056]95

fit |>
  forecast(h = 10, bootstrap = TRUE, times = 1000) |> 
  autoplot(datos_venta_bacalao, level = 95) +
  labs(y = "ventas", title = "Venta bacalao") 

Figure 29: Simulación de intervalos venta de bacalao

9 Modelos de Box-Jenkins y su identificacion tentativa

9.1 Primeras diferencias

datos_venta_rollos <- tsibble(
    date = rep(yearweek("2007 W1") + 0:119),
    value =  c(15, 14.4064, 14.9383, 16.0734, 15.6320, 14.3975, 13.8959, 14.0765,
               16.3750, 16.5342, 16.3839, 17.1006, 17.7876, 17.7354, 17.0010, 17.7485,
               18.1888, 18.5997, 17.5859, 15.7389, 13.6971, 15.0059, 16.2574, 14.3506,
               11.9515, 12.0328, 11.2142, 11.7023, 12.5905, 12.1991, 10.7752, 10.1123,
               9.9330, 11.7435, 12.2590, 12.5009, 11.5378, 9.6649, 10.1043, 10.3452,    
               9.2835,7.7219, 6.8300, 8.2046, 8.5289, 8.8733, 8.7948, 8.1577, 7.9128,                  8.7978, 9.0775, 9.3234, 10.4739, 10.6943, 9.8367, 8.1803, 7.2509, 5.0814,                1.8313, -0.9127, -1.3173, -0.6021, 0.1400, 1.4030, 1.9280, 2.5626, 1.9615, 4.8463, 6.5454, 8.0141, 7.9746, 8.4959, 8.4538, 8.7114, 7.3780, 8.1905, 9.9720, 9.6930, 9.4506, 11.2088, 11.4986, 13.2778, 13.5910, 13.4927, 13.3125, 12.7445, 11.7979, 11.7319, 11.6523, 11.3718, 10.5502, 11.4741, 11.5568, 11.7986, 11.8867, 11.2951,                  12.7847, 13.9435, 13.6859, 14.1136, 13.8949, 14.2853, 16.3867, 17.0884,                 15.8861, 14.8227, 15.9479, 15.0982, 13.8770, 14.2746, 15.1682, 15.3818,                 14.1863, 13.9996, 15.2463, 17.0179, 17.2929, 16.6366, 15.3410, 15.6453),
    index = date
    ) 
datos_venta_rollos |> autoplot(value) |> #ACF(value) |>
  labs(subtitle = "Venta semanal de rollos de papel")#labs(subtitle = "ACF venta semanal de rollos de papel")

## [[1]]

Figure 30: Ventas de rollos de papel del ejemplo 9.1

## 
## $subtitle
## [1] "Venta semanal de rollos de papel"
## 
## attr(,"class")
## [1] "labels"

datos_venta_rollos |> autoplot(difference(value)) + 
  labs(subtitle = "Primeras diferencias venta semanal de rollos de papel")

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

Figure 31: Primeras diferencias ventas de rollos de papel

9.2 Función de autocorrelación muestral

datos_venta_rollos |> ACF(value) |>
  autoplot() + labs(subtitle = "ACF en ventas de rollos de papel")

Figure 32: Autocorrelación venta de rollos de papel

datos_venta_rollos |> ACF(difference(value)) |>
  autoplot() + labs(subtitle = "ACF diferencias en ventas de rollos de papel")

Figure 33: Autocorrelación primeras diferencias

9.3 Función de autocorrelación parcial

datos_venta_rollos |> PACF(difference(value)) |>
  autoplot() + labs(subtitle = "PACF diferencias en ventas de rollos de papel")

Figure 34: Autocorrelación parcial primeras diferencias

9.4 Venta de rollos de papel

Estimación:

fit_rol <- datos_venta_rollos |>
  model(ARIMA(value ~ pdq(1,1,0)))
report(fit_rol)

## Series: value 
## Model: ARIMA(1,1,0) 
## 
## Coefficients:
##          ar1
##       0.3479
## s.e.  0.0855
## 
## sigma^2 estimated as 1.037:  log likelihood=-170.58
## AIC=345.16   AICc=345.27   BIC=350.72

Pronósticos:

fit_rol |> forecast(h=10) |>
  autoplot(datos_venta_rollos) +
  labs(y = "Ventas", title = "Ventas rollos de papel")

Figure 35: Pronósticos

9.5 Viscodad químico XB-55-7

datos_visc <- tsibble(
    Day = 1:95,
    value =  c(25, 27, 33.5142, 35.4962, 36.9029, 37.8359, 34.2654, 31.8978, 33.7567,
               36.6298, 36.3518, 40.0762, 38.0928, 34.5412, 34.8567, 34.5316, 32.3851,
               32.6058, 34.8913,
               38.2418, 36.8926, 33.8942, 34.1710, 35.4268, 38.5831, 34.6184, 33.9741,
               30.2072, 30.5429, 34.8686, 35.8892, 35.2035, 34.4337, 35.4844, 33.2381,
               36.1684, 34.4116, 33.7668, 
               33.4246, 33.5719, 35.9222, 33.2125, 37.1668, 35.8138, 33.6847, 33.2761,
               38.8162, 42.0838, 40.0069, 33.4514, 30.8413, 30.0655, 37.9075, 39.0982,
               37.9075, 36.2393, 34.9535, 
               33.2061, 34.4261, 37.4511, 37.3335, 38.4679, 33.0976, 32.9285, 32.2754,
               33.2214, 34.5786, 32.3448, 31.5136, 37.0844, 36.0536, 35.7297, 36.79912,
               34.9502, 33.5246, 35.1012, 
               35.9574, 38.0977, 33.4598, 32.9278, 36.5121, 37.4243, 35.1550, 34.4797,
               33.2898, 33.9252, 36.1036, 36.7351, 35.4576, 37.5924, 34.4895, 39.1692,
               35.8242, 32.3875, 31.2846),
    index = Day
    ) 
datos_visc |>
  autoplot(value) +
  labs(title = "Viscocidad XB-77-5", y = "viscocidad") + 
  geom_hline(yintercept = mean(datos_visc$value), linetype = "dashed", color = "red")

Figure 36: Viscocidad del químico XB-77-5

Estimación:

fit_qui <- datos_visc |>
  model(ARIMA(value ~ pdq(2,0,0)))
report(fit_qui)

## Series: value 
## Model: ARIMA(2,0,0) w/ mean 
## 
## Coefficients:
##          ar1      ar2  constant
##       0.6838  -0.4339   26.2141
## s.e.  0.0979   0.1037    0.2204
## 
## sigma^2 estimated as 4.689:  log likelihood=-207.01
## AIC=422.02   AICc=422.46   BIC=432.23

Pronósticos:

fit_qui |> forecast(h=10) |>
  autoplot(datos_visc) +
  labs(y = "Viscocidad", title = "Viscocidad diaria XB-77-5")

Figure 37: Pronósticos viscocidad

10 Estimación, chequeo de diagnóstico, y pronóstico de modelos de Box-Jenkins no estacionales

10.1 Estimación

Modelo autoregresivo(1,1,0):

datos_visc |> 
  ACF(value) |>
  autoplot() + labs(subtitle = "ACF venta de bacalao")

Figure 38: Modelo autoregresivo(1,1,0) viscocidad XB-55-7

Modelo autoregresivo(2,0,0):

datos_visc |>
  model(ARIMA(value ~ pdq(2,0,0))) |>
  report()

## Series: value 
## Model: ARIMA(2,0,0) w/ mean 
## 
## Coefficients:
##          ar1      ar2  constant
##       0.6838  -0.4339   26.2141
## s.e.  0.0979   0.1037    0.2204
## 
## sigma^2 estimated as 4.689:  log likelihood=-207.01
## AIC=422.02   AICc=422.46   BIC=432.23

10.2 Autocorrelación de los residuales

Venta de rollos de papel:

datos_venta_rollos |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  augment() |>
  ACF(.innov) |>
  autoplot() +
  labs(title = "ACF residuales autoregresivo(1,1,0)")

Figure 39: Autocorrelación residuales autoregresivo(1,1,0)

datos_venta_rollos |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  augment() |>
  PACF(.innov) |>
  autoplot() +
  labs(title = "PACF residuales autoregresivo(1,1,0)")

Figure 40: Autocorrelación parcial residuales autoregresivo(1,1,0)

Viscocidad XB-55-7:

datos_visc |>
  model(ARIMA(value ~ pdq(2,0,0))) |>
  augment() |>
  ACF(.innov) |>
  autoplot() +
  labs(title = "ACF residuales autoregresivo(2,0,0)")

Figure 41: Autocorrelación residuales autoregresivo(2,0,0)

datos_visc |>
  model(ARIMA(value ~ pdq(2,0,0))) |>
  augment() |>
  PACF(.innov) |>
  autoplot() +
  labs(title = "PACF residuales autoregresivo(2,0,0)")

Figure 42: Autocorrelación parcial residuales autoregresivo(2,0,0)

10.3 Intervalo de confianza

Modelo autoregresivo(1,1,0):

datos_visc |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  forecast(h=10) |>
  hilo(level = c(95))

## # A tsibble: 10 x 5 [1]
## # Key:       .model [1]
##    .model                        Day      value .mean                  `95%`
##    <chr>                       <dbl>     <dist> <dbl>                 <hilo>
##  1 ARIMA(value ~ pdq(1, 1, 0))    96 N(31, 6.8)  31.3 [26.13054, 36.38022]95
##  2 ARIMA(value ~ pdq(1, 1, 0))    97  N(31, 14)  31.3 [23.91035, 38.59886]95
##  3 ARIMA(value ~ pdq(1, 1, 0))    98  N(31, 21)  31.3 [22.21855, 40.29062]95
##  4 ARIMA(value ~ pdq(1, 1, 0))    99  N(31, 28)  31.3 [20.79691, 41.71226]95
##  5 ARIMA(value ~ pdq(1, 1, 0))   100  N(31, 36)  31.3 [19.54664, 42.96253]95
##  6 ARIMA(value ~ pdq(1, 1, 0))   101  N(31, 43)  31.3 [18.41757, 44.09160]95
##  7 ARIMA(value ~ pdq(1, 1, 0))   102  N(31, 50)  31.3 [17.38008, 45.12909]95
##  8 ARIMA(value ~ pdq(1, 1, 0))   103  N(31, 57)  31.3 [16.41495, 46.09422]95
##  9 ARIMA(value ~ pdq(1, 1, 0))   104  N(31, 65)  31.3 [15.50886, 47.00031]95
## 10 ARIMA(value ~ pdq(1, 1, 0))   105  N(31, 72)  31.3 [14.65215, 47.85702]95

Modelo autoregresivo(2,0,0):

datos_visc |>
  model(ARIMA(value ~ pdq(2,0,0))) |>
  forecast(h=10) |>
  hilo(level = c(95))

## # A tsibble: 10 x 5 [1]
## # Key:       .model [1]
##    .model                        Day      value .mean                  `95%`
##    <chr>                       <dbl>     <dist> <dbl>                 <hilo>
##  1 ARIMA(value ~ pdq(2, 0, 0))    96 N(34, 4.7)  33.6 [29.30936, 37.79748]95
##  2 ARIMA(value ~ pdq(2, 0, 0))    97 N(36, 6.9)  35.6 [30.44202, 40.72506]95
##  3 ARIMA(value ~ pdq(2, 0, 0))    98 N(36, 6.9)  36.0 [30.84377, 41.13078]95
##  4 ARIMA(value ~ pdq(2, 0, 0))    99 N(35, 7.2)  35.4 [30.10935, 40.65546]95
##  5 ARIMA(value ~ pdq(2, 0, 0))   100 N(35, 7.4)  34.8 [29.45142, 40.13572]95
##  6 ARIMA(value ~ pdq(2, 0, 0))   101 N(35, 7.4)  34.7 [29.31061, 39.99615]95
##  7 ARIMA(value ~ pdq(2, 0, 0))   102 N(35, 7.5)  34.8 [29.46093, 40.16513]95
##  8 ARIMA(value ~ pdq(2, 0, 0))   103 N(35, 7.5)  35.0 [29.62505, 40.34105]95
##  9 ARIMA(value ~ pdq(2, 0, 0))   104 N(35, 7.5)  35.0 [29.67191, 40.38814]95
## 10 ARIMA(value ~ pdq(2, 0, 0))   105 N(35, 7.5)  35.0 [29.62959, 40.34717]95

10.4 Modelo venta de bacalao

Estimación autoregresivo(1,0,0):

datos_venta_bacalao |>
  model(ARIMA(value ~ pdq(1,0,0))) |>
  report()

## Series: value 
## Model: ARIMA(1,0,0) w/ mean 
## 
## Coefficients:
##           ar1  constant
##       -0.0660  374.4266
## s.e.   0.2007    6.7636
## 
## sigma^2 estimated as 1191:  log likelihood=-118
## AIC=242   AICc=243.2   BIC=245.53

ACF:

datos_venta_bacalao |>
  ACF(value) |>
  autoplot() +
  labs(title = "ACF venta de bacalao")

Figure 43: Autocorrelación venta de bacalao

PACF:

datos_venta_bacalao |>
  PACF(value) |>
  autoplot() +
  labs(title = "PACF venta de bacalao")

Figure 44: Autocorrelación parcial venta bacalao

10.5 Modelo de termostatos

Estimación autoregresivo(1,1,0):

datos_venta_termos |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  report()

## Series: value 
## Model: ARIMA(1,1,0) 
## 
## Coefficients:
##           ar1
##       -0.3954
## s.e.   0.1351
## 
## sigma^2 estimated as 909.2:  log likelihood=-245.67
## AIC=495.33   AICc=495.58   BIC=499.2

ACF residuales:

datos_venta_termos |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  augment() |>
  ACF(.innov) |>
  autoplot() +
  labs(title = "ACF residuales autoregresivo(1,1,0)")

Figure 45: Autocorrelación residuales autoregresivo(1,1,0)

PACF residuales:

datos_venta_termos |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  augment() |>
  PACF(.innov) |>
  autoplot() +
  labs(title = "PACF residuales autoregresivo(1,1,0)")

Figure 46: Autocorrelación parcial residuales autoregresivo(1,1,0)

Pronósticos:

datos_venta_termos |>
  model(ARIMA(value ~ pdq(1,1,0))) |>
  forecast(h=10) |>
  autoplot(datos_venta_termos) +
  labs(y = "Ventas", title = "Venta de termos")

Figure 47: Pronósticos autoregresivo(1,1,0) venta de termos

11 Modelos de Box-Jenkins estacionales

11.1 Transformaciones estacionarias

ACF raíces cuartas promedios habitaciones:

promedio_cuartos |> ACF(value**.25) |>
  autoplot() + labs(subtitle = "ACF raíces cuartas")

Figure 48: Autocorrelación raíces cuartas prmedios habitaciones

PACF raíces cuartas promedios habitaciones:

promedio_cuartos |> PACF(value**.25) |>
  autoplot() + labs(subtitle = "PACF raíces cuartas")

Figure 49: Autocorrelación parcial raíces cuartas prmedios habitaciones

ACF primeras diferencias rc promedios habitaciones:

promedio_cuartos |> ACF(difference(value**.25)) |>
  autoplot() + labs(subtitle = "ACF diferencias promedios habitaciones")

Figure 50: Autocorrelación primeras diferencias promedios habitaciones

PACF primeras diferencias promedios habitaciones:

promedio_cuartos |> PACF(difference(value**.25)) |>
  autoplot() + labs(subtitle = "PACF diferencias promedios habitaciones")

Figure 51: Autocorrelación parcial primeras diferencias promedios habitaciones

ACF primera diferencia estacional promedios habitaciones:

promedio_cuartos |> ACF(difference(value**.25, lag=12)) |>
  autoplot() + labs(subtitle = "ACF diferencia estacional promedios habitaciones")

Figure 52: Autocorrelación primera diferencia estacional promedios habitaciones

PACF primera diferencia estacional promedios habitaciones:

promedio_cuartos |> PACF(difference(value**.25, lag=12)) |>
  autoplot() + labs(subtitle = "PACF diferencia estacional promedios habitaciones")

Figure 53: Autocorrelación prcial primera diferencia estacional promedios habitaciones

ACF primeras diferencias combinadas promedios habitaciones:

promedio_cuartos |>
  mutate(
    diff_regular = difference(value**.25),
    diff_estacional = difference(value**.25, lag = 12),
    diff_combined = difference(value**.25) - lag(difference(value**.25, lag = 12), 1)
  ) |>
  ACF(diff_combined) |>
  autoplot() + 
  labs(subtitle = "ACF de las diferencias combinadas (regulares y estacionales)")

Figure 54: Autocorrelación diferencias combinadas promedios habitaciones

PACF primeras diferencias combinadas promedios habitaciones:

promedio_cuartos |>
  mutate(
    diff_regular = difference(value**.25),
    diff_estacional = difference(value**.25, lag = 12),
    diff_combined = difference(value**.25) - lag(difference(value**.25, lag = 12), 1)
  ) |>
  PACF(diff_combined) |>
  autoplot() + 
  labs(subtitle = "PACF de las diferencias combinadas")

Figure 55: Autocorrelación parcial diferencias combinadas promedios habitaciones

11.2 Modelo promedios habitaciones

Estimación:

pro_cuartos <- promedio_cuartos |>
  mutate(y_star = value^(1/4),
         z_t = y_star - lag(y_star, 12)) |>
   filter(!is.na(z_t)) 
fit_pro <- pro_cuartos |>
  model(
        ARIMA(z_t ~ pdq(5, 0, 0) + PDQ(0, 0, 1))
    ) 
fit_pro |>
  report()

## Series: z_t 
## Model: ARIMA(5,0,0)(0,0,1)[12] w/ mean 
## 
## Coefficients:
##          ar1     ar2      ar3      ar4      ar5     sma1  constant
##       0.2454  0.0899  -0.2299  -0.0239  -0.1112  -0.5594    0.0436
## s.e.  0.0836  0.0845   0.0813   0.0830   0.0821   0.0856    0.0010
## 
## sigma^2 estimated as 0.0006619:  log likelihood=350.78
## AIC=-685.56   AICc=-684.58   BIC=-661.16

ACF residuales:

fit_pro |> 
  augment() |>
  ACF(.innov) |>
  autoplot() +
  labs(title = "ACF residuales ARIMA estacional(5,0,0)(0,0,1)")

Figure 56: Autocorrelación residuales ARIMA estacional(5,0,0)(0,0,1)

PACF residuales:

fit_pro |> 
  augment() |>
  PACF(.innov) |>
  autoplot() +
  labs(title = "PACF residuales ARIMA estacional(5,0,0)(0,0,1)")

Figure 57: Autocorrelación parcial residuales ARIMA estacional(5,0,0)(0,0,1)

Pronósticos:

forecast(fit_pro, h=36) |>
  #filter(.model=='auto') |>
  autoplot(pro_cuartos) +
  labs(title = "Promedios ARIMA estacional",
       y="promedios")

Figure 58: Pronósticos ARIMA estacional promedios habitaciones

11.3 Modelo pasajeros vuelos

Transformación LN:

datos_pasajeros <- tsibble(
    date = rep(yearmonth("1983 Jan") + 0:131),
    values = c(
        112, 118, 132, 129, 121, 135, 148, 148, 136, 119, 104, 118, # Año 1
        115, 126, 141, 135, 125, 149, 170, 170, 158, 133, 114, 140, # Año 2
        145, 150, 178, 163, 172, 178, 199, 199, 184, 162, 146, 166, # Año 3
        171, 180, 193, 181, 183, 218, 230, 242, 209, 191, 172, 194, # Año 4
        196, 196, 236, 235, 229, 243, 264, 272, 237, 211, 180, 201, # Año 5
        204, 188, 235, 227, 234, 264, 302, 293, 259, 229, 203, 229, # Año 6
        242, 233, 267, 269, 270, 315, 364, 347, 312, 274, 237, 278, # Año 7
        284, 277, 317, 313, 318, 374, 413, 405, 355, 306, 271, 306, # Año 8
        315, 301, 356, 348, 355, 422, 465, 467, 404, 347, 305, 336, # Año 9
        340, 318, 362, 348, 363, 435, 491, 505, 404, 359, 310, 337, # Año 10
        360, 342, 406, 396, 420, 472, 548, 559, 463, 407, 362, 405  # Año 11
    ),
    index = date
) 
datos_pasajeros <- datos_pasajeros |> 
  mutate(ln_value = log(values))
datos_pasajeros |> 
  autoplot(ln_value) +
  labs(y = "LN pasajeros (mensuales)",
       title = "Pasajeors vuelos internacionales")

Figure 59: Modelo ARIMA estacional(5,0,0)(0,0,1)

ACF LN pasajeros:

datos_pasajeros |> ACF(ln_value) |>
  autoplot() + labs(subtitle = "ACF LN pasajeors")

Figure 60: Autocorrelación LN pasajeros vuelos

ACF primeras diferencias LN pasajeros:

datos_pasajeros |> ACF(difference(ln_value)) |>
  autoplot() + labs(subtitle = "ACF diferencias LN pasajeros")

Figure 61: Autocorrelación primeras diferencias LN pasajeros vuelos

ACF primera diferencia estacional LN pasajeros:

datos_pasajeros |> ACF(difference(ln_value, lag = 12)) |>
  autoplot() + labs(subtitle = "ACF diferencia estacional LN pasajeros")

Figure 62: Autocorrelación primera diferencia estacional LN pasajeros vuelos

Estimación:

fit_pas <- datos_pasajeros |>
  model(
        ARIMA(ln_value ~ pdq(0, 1, 1) + PDQ(0, 1, 1))
    ) 
fit_pas |>
  report()

## Series: ln_value 
## Model: ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.3484  -0.5623
## s.e.   0.0943   0.0774
## 
## sigma^2 estimated as 0.001338:  log likelihood=223.63
## AIC=-441.26   AICc=-441.05   BIC=-432.92

Pronósticos:

forecast(fit_pas, h=36) |>
  #filter(.model=='auto') |>
  autoplot(datos_pasajeros) +
  labs(title = "Pasajeros ARIMA estacional",
       y="pasajeros")

Figure 63: Pronósticos ARIMA estacional pasajeros

11.4 Modelo de variables ficticias y residuales ARIMA

Transformación variables ficticias:

# Transformar la serie usando la raíz cuarta
pro_habs <- promedio_cuartos |>
  mutate(transformed_value = value^(1/4),
  month = month(date, label = TRUE),
  t = row_number()
  ) 
pro_habs <- pro_habs |>
  pivot_wider(names_from = month, values_from = month, values_fill = 0, names_prefix = "M_", values_fn = list(month = ~ 1))

# Ajustar el modelo de regresión con variables ficticias y tendencia lineal
modelo_regresion <- lm(transformed_value ~ t + M_ene + M_feb + M_mar + M_abr + M_may + M_jun + M_jul + M_ago + M_sep + M_oct + M_nov, data = pro_habs)

# Resumen del modelo de regresión
summary(modelo_regresion)

## 
## Call:
## lm(formula = transformed_value ~ t + M_ene + M_feb + M_mar + 
##     M_abr + M_may + M_jun + M_jul + M_ago + M_sep + M_oct + M_nov, 
##     data = pro_habs)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.068385 -0.018929  0.001583  0.018820  0.079295 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.810e+00  8.582e-03 560.547  < 2e-16 ***
## t            3.511e-03  4.512e-05  77.826  < 2e-16 ***
## M_ene       -5.527e-02  1.070e-02  -5.164 7.33e-07 ***
## M_feb       -1.436e-01  1.070e-02 -13.418  < 2e-16 ***
## M_mar       -1.099e-01  1.070e-02 -10.271  < 2e-16 ***
## M_abr        4.661e-02  1.070e-02   4.357 2.39e-05 ***
## M_may        2.263e-02  1.070e-02   2.115    0.036 *  
## M_jun        1.874e-01  1.070e-02  17.520  < 2e-16 ***
## M_jul        3.797e-01  1.069e-02  35.503  < 2e-16 ***
## M_ago        4.106e-01  1.069e-02  38.397  < 2e-16 ***
## M_sep        6.864e-02  1.069e-02   6.420 1.58e-09 ***
## M_oct        4.787e-02  1.069e-02   4.477 1.46e-05 ***
## M_nov       -1.447e-01  1.069e-02 -13.534  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02829 on 155 degrees of freedom
## Multiple R-squared:  0.988,  Adjusted R-squared:  0.9871 
## F-statistic:  1066 on 12 and 155 DF,  p-value: < 2.2e-16

Pŕonósticos:

pronosticos <- predict(modelo_regresion, newdata = pro_habs)
# Crear un data frame con los pronósticos y la fecha correspondiente
pronosticos_df <- data.frame(
  fecha = pro_habs$date,  # Utilizamos la fecha original
  pronostico = pronosticos^(4)  # Revertimos la transformación a la raíz cuarta
)

# Graficar los pronósticos
plot(pronosticos_df$fecha, pronosticos_df$pronostico, type = "l", 
     xlab = "Fecha", ylab = "Pronóstico", main = "Pronósticos del modelo de regresión",
     col = "blue", ylim = c(min(pronosticos_df$pronostico), max(pronosticos_df$pronostico)))

# Agregar la línea de los datos originales para comparación si se desea
lines(pro_habs$date, pro_habs$transformed_value, type = "l", col = "red")

# Agregar leyendas y etiquetas
legend("topleft", legend = c("Pronóstico", "Datos originales"), col = c("blue", "red"), lty = 1)

Figure 64: Pronósticos variables ficticias

# Mostrar el gráfico

ARIMA residuales:

# Extraer los residuos del modelo de regresión
residuos <- residuals(modelo_regresion)

# Ajustar el modelo ARIMA a los residuos
residuos_tsibble <- tsibble(date = pro_habs$date, residuos = residuos)

## Using `date` as index variable.

modelo_arima_residuos <- residuos_tsibble |>
  model(
    arima = ARIMA(residuos ~ pdq(5, 0, 0) + PDQ(1, 0, 0, period = 12))
  )

# Mostrar el informe del modelo ARIMA ajustado a los residuos
report(modelo_arima_residuos)

## Series: residuos 
## Model: ARIMA(5,0,0)(1,0,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ar3      ar4      ar5    sar1
##       0.2904  0.0841  -0.2317  -0.0455  -0.0836  0.2678
## s.e.  0.0807  0.0824   0.0795   0.0810   0.0785  0.0885
## 
## sigma^2 estimated as 0.0005673:  log likelihood=391.92
## AIC=-769.83   AICc=-769.13   BIC=-747.97

ACF residuales:

modelo_arima_residuos |>
  augment() |>
  ACF(.innov) |>
autoplot() +
labs(title = "ACF residuales variables ficticias")

Figure 65: ACF residuales ARIMA variables fictcias

PACF residuales:

modelo_arima_residuos |>
  augment() |>
  PACF(.innov) |>
autoplot() +
labs(title = "PACF residuales variables ficticias")

Figure 66: PACF residuales ARIMA variables fictcias

Pronósticos:

forecast(modelo_arima_residuos, h=36) |>
  #filter(.model=='auto') |>
  autoplot(pro_habs) +
  labs(title = "Pronósticos ARIMA variables ficticias",
       y="promedios")

Figure 67: Pronósticos ARIMA variables ficticias

Pronósticos combinadoos:

library(forecast)
# Ajustar el modelo ARIMA a los residuos del modelo de regresión
modelo_arima_residuos <- Arima(resid(modelo_regresion), order = c(1,0,0))

# Generar pronósticos para el modelo de regresión
pronosticos_regresion <- predict(modelo_regresion, newdata = pro_habs)

# Generar pronósticos para el modelo ARIMA sobre los residuos
pronosticos_arima_residuos <- forecast(modelo_arima_residuos, h = 36)

# Sumar los pronósticos de ambos modelos para obtener los pronósticos combinados
pronosticos_combinados <- pronosticos_regresion + fitted(pronosticos_arima_residuos)

combinados_df <- data.frame(
  fecha = pro_habs$date,  # Utilizamos la fecha original
  pronostico = pronosticos_combinados^(4)  # Revertimos la transformación a la raíz cuarta
  )

# Graficar los pronósticos
plot(combinados_df$fecha, pronosticos_df$pronostico, type = "l", 
     xlab = "Fecha", ylab = "Pronóstico combinado", main = "Pronósticos combinados del modelo de regresión",
     col = "blue", ylim = c(min(pronosticos_df$pronostico), max(pronosticos_df$pronostico)))
 
# Agregar la línea de los datos originales para comparación si se desea
lines(pro_habs$date, pro_habs$transformed_value, type = "l", col = "red")
 
# Agregar leyendas y etiquetas
legend("topleft", legend = c("Pronóstico combinado", "Datos originales"), col = c("blue", "red"), lty = 1)

Figure 68: Pronósticos combinados ARIMA variables ficticias

12 Modelos de Box-Jenkins avanzados

12.1 Consulta al directorio de Cincinatti

Datos:

datos_consultas <- tsibble(
    date = rep(yearmonth("1962 Jan") + 0:179),
     value = c(
    350, 339, 351, 364, 369, 331, 331, 340, 346, 341, 357, 398, # Año 1
    381, 367, 383, 375, 353, 361, 375, 371, 373, 366, 382, 429, # Año 2
    406, 403, 429, 425, 427, 409, 402, 409, 419, 404, 429, 463, # Año 3
    428, 449, 444, 467, 474, 463, 432, 453, 462, 456, 474, 514, # Año 4
    489, 475, 492, 525, 527, 533, 527, 522, 526, 513, 564, 599, # Año 5
    572, 587, 599, 601, 611, 620, 579, 582, 592, 581, 630, 633, # Año 6
    638, 631, 645, 682, 601, 595, 521, 521, 516, 496, 538, 575, # Año 7
    537, 534, 542, 538, 547, 540, 526, 548, 555, 545, 594, 643, # Año 8
    625, 616, 640, 625, 637, 634, 621, 641, 654, 649, 662, 699, # Año 9
    672, 704, 700, 711, 715, 718, 652, 664, 695, 704, 733, 772, # Año 10
    716, 712, 732, 755, 761, 748, 748, 750, 744, 731, 782, 810, # Año 11
    777, 816, 840, 868, 872, 811, 810, 762, 634, 626, 649, 697, # Año 12
    657, 549, 162, 177, 175, 162, 161, 165, 170, 172, 178, 186, # Año 13
    178, 178, 189, 205, 202, 185, 193, 200, 196, 204, 206, 227, # Año 14 
    225, 217, 219, 236, 253, 213, 205, 210, 216, 218, 235, 241  # Año 15
  ),
    index = date
) 
datos_consultas |> 
  autoplot(value) +
  labs(y = "Consultas (mensuales)",
       title = "Consultas directorio Cincinatti")

Figure 69: Llamadas de consulta al directorio de Cincinnati; promedio mensual por día (en unidades de 100 llamadas), desde enero de 1962 hasta diciembre de 1976

Estimación:

# Crear la variable ficticia S_t
datos_consultas <- datos_consultas |>
  mutate(S_t = if_else(date >= yearmonth("1974 Mar"), 1, 0))
# Diferenciar la serie
diff_yt <- diff(diff(datos_consultas$value, lag = 1), lag = 12)
diff_St <- diff(diff(datos_consultas$S_t, lag = 1), lag = 12)
# Crear un data frame con las series diferenciadas
data_diff <- data.frame(z_t = diff_yt, z_t_S = diff_St)
# Modelo de regresión con la variable ficticia diferenciada
reg_model_diff <- lm(z_t ~ z_t_S, data = data_diff)
summary(reg_model_diff)

## 
## Call:
## lm(formula = z_t ~ z_t_S, data = data_diff)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -146.898  -10.898    0.102   11.102  133.102 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -0.1018     2.3564  -0.043    0.966    
## z_t_S       -404.5000    21.5327 -18.785   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.45 on 165 degrees of freedom
## Multiple R-squared:  0.6814, Adjusted R-squared:  0.6795 
## F-statistic: 352.9 on 1 and 165 DF,  p-value: < 2.2e-16

# Residuos del modelo de regresión diferenciado
residuos_diff <- resid(reg_model_diff)
# Ajustar un modelo ARIMA a los residuos
arima_residuos_diff <- Arima(residuos_diff, order = c(0, 0, 1), seasonal = list(order = c(0, 0, 1), period = 12))
summary(arima_residuos_diff)

## Series: residuos_diff 
## ARIMA(0,0,1)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##           ma1     sma1     mean
##       -0.0390  -1.0000  -0.2941
## s.e.   0.0778   0.1247   0.3550
## 
## sigma^2 = 458.9:  log likelihood = -763.42
## AIC=1534.83   AICc=1535.08   BIC=1547.31
## 
## Training set error measures:
##                     ME    RMSE      MAE      MPE     MAPE      MASE
## Training set 0.3822045 21.2282 13.61204 173.8626 262.6699 0.4502985
##                       ACF1
## Training set -0.0006298984

Pronósticos:

# Pronosticar los residuos
forecast_residuos_diff <- forecast(arima_residuos_diff, h = 34)
# Obtener la predicción de la regresión
pred_reg_diff <- predict(reg_model_diff, newdata = data_diff[147:180, ])
# Combinar las predicciones
pred_total_diff <- pred_reg_diff + forecast_residuos_diff$mean
# Graficar los resultados
autoplot(pred_total_diff) + autolayer(forecast_residuos_diff, series = "Residuos pronosticados")

Figure 70: Pronóstico residuos modelo de intervención

12.2 Pasta dental ventas semanales

# Load necessary libraries
library(tidyverse)
library(tsibble)
library(feasts)
library(fable)

# Prepare the data
data <- tibble(
  t = c(1:90),
  sales = c(235.000, 239.000, 244.090, 252.731, 264.377, 277.934, 286.687, 295.629, 
            310.444, 325.112, 336.291, 344.459, 355.399, 367.691, 384.003, 398.042, 
            412.969, 422.901, 434.960, 445.853, 455.929, 465.584, 477.894, 491.408, 
            507.712, 517.237, 524.349, 532.104, 538.097, 544.948, 551.925, 557.929, 
            564.285, 572.164, 582.926, 595.295, 607.028, 617.541, 622.941, 633.436, 
            647.371, 658.230, 670.777, 685.457, 690.992, 693.557, 700.675, 712.710, 
            726.513, 736.429, 743.203, 751.227, 764.265, 777.852, 791.070, 805.844, 
            815.122, 822.905, 830.663, 839.600, 846.962, 853.830, 860.840, 871.075, 
            877.792, 881.143, 884.226, 890.208, 894.966, 901.288, 913.138, 922.511, 
            930.786, 941.306, 950.305, 952.373, 960.042, 968.100, 972.477, 977.408, 
            977.602, 979.505, 982.934, 985.833, 991.350, 996.291, 1003.100, 1010.320, 
            1018.420, 1029.48)
)
# Convert to tsibble
sales_tsibble <- data %>%
  as_tsibble(index = t)

# Plot the data to check for stationarity
autoplot(sales_tsibble, sales) +
  labs(title = "Weekly Sales of Ultra Shine Toothpaste", y = "Sales (in 1000 tubes)", x = "Week")

Figure 71: Ventas semanales de la pasta dental Ultra Shine

La visualización inicial de las ventas semanales de la pasta dental Ultra Shine muestra una tendencia creciente a lo largo del tiempo. Esto sugiere que los datos pueden no ser estacionarios, ya que la media y la varianza cambian con el tiempo.

sales_tsibble %>%
  ACF(sales) %>%
  autoplot() +
  labs(title = "ACF of Weekly Sales")

Figure 72: ACF y PACF ventas semanales Ultra Shine

sales_tsibble %>%
  PACF(sales) %>%
  autoplot() +
  labs(title = "PACF of Weekly Sales")

Figure 73: ACF y PACF ventas semanales Ultra Shine

Las gráficas ACF y PACF muestran autocorrelaciones significativas en varios rezagos, lo que confirma que los datos no son estacionarios. La ACF disminuye lentamente y la PACF muestra algunos cortes significativos, lo que indica que los datos tienen una tendencia.

sales_diff <- sales_tsibble %>%
  mutate(sales_diff = difference(sales)) %>%
  filter(!is.na(sales_diff))

autoplot(sales_diff, sales_diff) +
  labs(title = "Differenced Weekly Sales of Ultra Shine Toothpaste", y = "Differenced Sales (in 1000 tubes)", x = "Week")

Figure 74: Primeras diferencias ventas semanales Ultra Shine

La serie diferenciada muestra que la tendencia ha sido eliminada y los datos parecen más estacionarios. Esto se verifica mediante las gráficas ACF y PACF de los datos diferenciados.

sales_diff %>%
  ACF(sales_diff) %>%
  autoplot() +
  labs(title = "ACF of Differenced Weekly Sales")

Figure 75: ACF y PACF primeras diferencias ventas semanales Ultra Shine

sales_diff %>%
  PACF(sales_diff) %>%
  autoplot() +
  labs(title = "PACF of Differenced Weekly Sales")

Figure 76: ACF y PACF primeras diferencias ventas semanales Ultra Shine

La ACF de los datos diferenciados muestra un corte rápido después de unos pocos rezagos, y la PACF muestra algunos picos significativos. Esto sugiere que un modelo ARIMA puede ser adecuado.

fit <- sales_diff %>%
  model(ARIMA(sales_diff))

report(fit)

## Series: sales_diff 
## Model: ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1      ma2
##       0.2945  -0.4972  -0.4009
## s.e.  0.1826   0.1658   0.1316
## 
## sigma^2 estimated as 7.514:  log likelihood=-212.68
## AIC=433.36   AICc=433.84   BIC=443.27

El modelo ajustado es un ARIMA(1,1,2). Los coeficientes del modelo son:

• ar1 = 0.2945

• ma1 = -0.4972

• ma2 = -0.4009

El logaritmo de verosimilitud es -212.68, y los valores de AIC, AICc y BIC son 433.36, 433.84 y 443.27 respectivamente.

fit %>%
  gg_tsresiduals()

Figure 77: Validación ARIMA primeras diferencias ventas semanales Ultra Shine

fit %>%
  augment() %>%
  features(.resid, ljung_box, lag = 10)

## # A tibble: 1 × 3
##   .model            lb_stat lb_pvalue
##   <chr>               <dbl>     <dbl>
## 1 ARIMA(sales_diff)    3.65     0.962

La gráfica de residuos muestra que los residuos del modelo no presentan patrones obvios y parecen ser ruido blanco.

La prueba de Ljung-Box para los residuos da un valor de lb_stat=3.65lb_stat=3.65 con p-value=0.962p-value=0.962. Esto indica que no hay autocorrelación significativa en los residuos, confirmando que el modelo es adecuado.

forecast <- fit %>%
  forecast(h = 10)

autoplot(forecast) +
  labs(title = "Forecast of Weekly Sales of Ultra Shine Toothpaste", y = "Sales (in 1000 tubes)", x = "Week")

Figure 78: Pronósticos ARIMA primeras diferencias ventas semanales Ultra Shine

El pronóstico para los próximos 10 periodos (semanas 91 a 100) muestra un aumento continuo en las ventas semanales de la pasta dental Ultra Shine. La visualización incluye intervalos de predicción al 95%.

El modelo ARIMA(1,1,2) es adecuado para modelar las ventas semanales de la pasta dental Ultra Shine después de diferenciar los datos para hacerlos estacionarios. Los análisis de ACF y PACF, así como la validación mediante residuos y la prueba de Ljung-Box, indican que el modelo ajustado es apropiado. El pronóstico sugiere que las ventas seguirán aumentando en el futuro cercano.

Referencias

Bowerman, Bruce L, Richard T O’Connell, and Anne B Koehler. 2005. Forecasting, Time Series, and Regression: An Applied Approach. South-Western College.

Hyndman, R, G Athanasopoulos, C Bergmeir, G Caceres, L Chhay, M O’Hara-Wild, F Petropoulos, S Razbash, E Wang, and F Yasmeen. 2024. forecast: Forecasting Functions for Time Series and Linear Models. https://pkg.robjhyndman.com/forecast/.

Hyndman, Rob. 2023. fpp3: Data for "Forecasting: Principles and Practice" (3rd Edition). https://CRAN.R-project.org/package=fpp3.

Hyndman, Rob J, and George Athanasopoulos. 2021. Forecasting: Principles and Practice. 3rd ed. Australia: OTexts.

Hyndman, Rob J, and Yeasmin Khandakar. 2008. “Automatic Time Series Forecasting: The forecast Package for R.” J. Stat. Softw. 27 (3).

R Core Team. 2023. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.