Actividad 2. Modelos de pronósticos en series de tiempo

Series de Tiempo

Cuatro ejercicios. Los dos primeros comparan cinco modelos de pronóstico sobre una misma serie y eligen el mejor por MAPE. Los dos últimos son informes gerenciales con descomposición estacional.

Instalar paquetes y llamar librerías

#install.packages("forecast")
#install.packages("readxl")
library(forecast)
library(readxl)

Funciones que se usan en todo el documento

# Promedio Móvil de orden k. El pronóstico de t usa las k observaciones previas.
pm <- function(y, k){
  f <- rep(NA, length(y))
  for(t in (k+1):length(y)) f[t] <- mean(y[(t-k):(t-1)])
  f
}

# Promedio Móvil Ponderado. El primer peso del vector va al dato inmediatamente
# anterior, no al más viejo. Los pesos deben sumar 1.
pmp <- function(y, w){
  k <- length(w); f <- rep(NA, length(y))
  for(t in (k+1):length(y)) f[t] <- sum(y[(t-1):(t-k)] * w)
  f
}

# Suavizador Exponencial. Arranca con F2 = y1.
suav <- function(y, a){
  f <- rep(NA, length(y)); f[2] <- y[1]
  for(t in 3:length(y)) f[t] <- a*y[t-1] + (1-a)*f[t-1]
  f
}
suav_sig <- function(y, a){
  f <- y[1]; for(t in 2:length(y)) f <- a*y[t] + (1-a)*f; f
}

# Las seis métricas. Cada modelo se divide entre SUS propios errores.
métricas <- function(real, pron){
  ok <- !is.na(pron) & !is.na(real)
  y <- real[ok]; f <- pron[ok]; e <- y - f; n <- length(e)
  c(n=n, ME=sum(e)/n, MAE=sum(abs(e))/n, MSE=sum(e^2)/n,
    RMSE=sqrt(sum(e^2)/n), MPE=sum(e/y*100)/n, MAPE=sum(abs(e/y*100))/n)
}

# Índices estacionales por el método de razón al promedio móvil centrado.
# Se normalizan para que sumen el número de periodos del ciclo.
índices_estacionales <- function(serie){
  f <- frequency(serie)
  cma <- stats::filter(serie, c(0.5, rep(1, f-1), 0.5)/f, sides = 2)
  razón <- as.numeric(serie) / as.numeric(cma)
  crudos <- tapply(razón, cycle(serie), mean, na.rm = TRUE)
  crudos * f / sum(crudos)
}

---

Ejercicio 1. Ventas Semanales

Instrucciones: modela la serie con los cinco modelos vistos en clase, arma la tabla resumen y concluye.

semana <- c(1:12)
valor  <- c(17,21,19,23,18,16,20,18,22,20,15,22)
ts1 <- ts(valor, c(2025,1), frequency=52)
plot(ts1, ylab="Valor", xlab="Semana", main="Ventas semanales")

# Modelo 1. Naive
naive1 <- naive(ts1, h=6); summary(naive1)

## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = ts1, h = 6) 
## 
## Residual sd: 4.0339 
## 
## Error measures:
##                     ME     RMSE      MAE       MPE     MAPE MASE       ACF1
## Training set 0.4545455 4.033947 3.727273 0.1082169 19.24431  NaN -0.4553872
## 
## Forecasts:
##          Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 2025.231             22 16.830289 27.16971 14.093609 29.90639
## 2025.250             22 14.688925 29.31108 10.818675 33.18132
## 2025.269             22 13.045798 30.95420  8.305730 35.69427
## 2025.288             22 11.660578 32.33942  6.187219 37.81278
## 2025.308             22 10.440175 33.55983  4.320773 39.67923
## 2025.327             22  9.336846 34.66315  2.633377 41.36662

f1_naive <- c(NA, head(valor, -1))

# Modelo 2. Promedio Móvil (k=3)
f1_pm <- pm(valor, 3)

# Modelo 3. Promedio Móvil Ponderado (k=3, pesos 3/6, 2/6 y 1/6)
f1_pmp <- pmp(valor, c(3,2,1)/6)

# Modelo 4. Suavizador Exponencial (alfa = 0.2)
f1_ses <- suav(valor, 0.2)

# Modelo 5. ARIMA
arima1 <- auto.arima(ts1); arima1

## Series: ts1 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##          mean
##       19.2500
## s.e.   0.6985
## 
## sigma^2 = 6.386:  log likelihood = -27.63
## AIC=59.26   AICc=60.59   BIC=60.23

f1_arima <- as.numeric(fitted(arima1))

tabla1 <- data.frame(Semana=semana, Real=valor, Naive=f1_naive, PM3=f1_pm,
                     PMP3=f1_pmp, SES=f1_ses, ARIMA=f1_arima)
knitr::kable(tabla1, digits=2, caption="Pronósticos de cada modelo, Ejercicio 1")

Pronósticos de cada modelo, Ejercicio 1

Semana	Real	Naive	PM3	PMP3	SES	ARIMA
1	17	NA	NA	NA	NA	19.25
2	21	17	NA	NA	17.00	19.25
3	19	21	NA	NA	17.80	19.25
4	23	19	19	19.33	18.04	19.25
5	18	23	21	21.33	19.03	19.25
6	16	18	20	19.83	18.83	19.25
7	20	16	19	17.83	18.26	19.25
8	18	20	18	18.33	18.61	19.25
9	22	18	18	18.33	18.49	19.25
10	20	22	20	20.33	19.19	19.25
11	15	20	20	20.33	19.35	19.25
12	22	15	19	17.83	18.48	19.25

modelos1 <- list("1. Naive"=f1_naive, "2. PM k=3"=f1_pm, "3. PMP k=3"=f1_pmp,
                 "4. SES alfa=0.2"=f1_ses, "5. ARIMA"=f1_arima)
res1 <- t(sapply(modelos1, métricas, real=valor))
knitr::kable(res1, digits=4, caption="Tabla resumen de resultados, Ejercicio 1")

Tabla resumen de resultados, Ejercicio 1

	n	ME	MAE	MSE	RMSE	MPE	MAPE
1. Naive	11	0.4545	3.7273	16.2727	4.0339	0.1082	19.2443
2. PM k=3	9	0.0000	2.6667	10.2222	3.1972	-2.3101	14.3566
3. PMP k=3	9	0.0556	2.9815	11.4907	3.3898	-2.1299	15.9925
4. SES alfa=0.2	11	0.9932	2.5963	8.9822	2.9970	3.2600	13.4024
5. ARIMA	12	0.0000	2.0833	5.8542	2.4195	-1.6623	11.1853

sig1 <- c(Naive=tail(valor,1), PM3=mean(tail(valor,3)),
          PMP3=sum(rev(tail(valor,3))*c(3,2,1)/6), SES=suav_sig(valor,0.2),
          ARIMA=as.numeric(forecast(arima1, h=1)$mean))
round(sig1, 4)

##   Naive     PM3    PMP3     SES   ARIMA 
## 22.0000 19.0000 19.3333 19.1850 19.2500

Conclusión del Ejercicio 1

El MAPE más bajo lo da el ARIMA con 11.19%, y atrás van el suavizador exponencial con 13.40%, el promedio móvil simple con 14.36%, el ponderado con 15.99% y hasta el final el Naive con 19.24%.

Vale la pena ver qué modelo escogió el auto.arima, porque salió un ARIMA(0,0,0) con media. Eso quiere decir que decidió que la serie es puro ruido alrededor de un nivel fijo y que su pronóstico para todas las semanas es el mismo número, 19.25. O sea que el mejor modelo acabó siendo no hacerle caso a nada de lo reciente y quedarse en el promedio.

Eso encaja con el resto de la tabla. El orden va con el peso que cada modelo le da al último dato, porque el Naive le da todo y es el peor, el ponderado la mitad, el simple un tercio, el suavizador una quinta parte y el ARIMA prácticamente nada. Entre menos caso le hacen al último dato, mejor pronostican. La razón es que la serie se mueve alrededor de 19 sin subir ni bajar, entonces el brinco de una semana a otra es ruido y creerle nada más mete error.

El ARIMA se evalúa sobre las 12 semanas, mientras que los promedios móviles pierden 3 y el Naive pierde 1, y su media está estimada con los mismos datos con los que se le mide el error, así que la comparación no es del todo pareja. Aun con eso, los cinco modelos apuntan a la misma conclusión.

---

Ejercicio 2. Hershey’s

Instrucciones: mismo análisis que el ejercicio 1 pero con la base Ventas_Historicas_Lechitas.xlsx. Pronostica los siguientes 6 meses.

# El Excel trae los tres años en bloques de dos columnas uno junto al otro
# (2017 en A-B, 2018 en C-D, 2019 en E-F), por eso se apilan.
ruta <- "Ventas_Historicas_Lechitas.xlsx"
crudo <- read_excel(ruta, col_names=FALSE)

bloques <- list(crudo[,1:2], crudo[,3:4], crudo[,5:6])
df2 <- do.call(rbind, lapply(bloques, function(b){
  names(b) <- c("Mes","Ventas")
  b <- b[!is.na(suppressWarnings(as.numeric(b$Ventas))), ]
  data.frame(Mes=as.Date(as.numeric(b$Mes), origin="1899-12-30"),
             Ventas=as.numeric(b$Ventas))
}))
df2 <- df2[order(df2$Mes), ]
ventas <- df2$Ventas
ts2 <- ts(ventas, c(2017,1), frequency=12)
str(df2)

## 'data.frame':    36 obs. of  2 variables:
##  $ Mes   : Date, format: "2017-01-01" "2017-02-01" ...
##  $ Ventas: num  25521 23740 26254 25868 27073 ...

plot(ts2, ylab="Ventas en miles de USD", xlab="Año",
     main="Ventas históricas de leche saborizada Hershey México")

# Modelo 1. Naive
naive2 <- naive(ts2, h=6); summary(naive2)

## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = ts2, h = 6) 
## 
## Residual sd: 1112.6435 
## 
## Error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 266.4474 1112.644 902.8514 0.8162353 3.013136 0.2579763 -0.5648044
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2020       34846.17 33420.26 36272.08 32665.43 37026.91
## Feb 2020       34846.17 32829.63 36862.71 31762.14 37930.20
## Mar 2020       34846.17 32376.42 37315.92 31069.02 38623.32
## Apr 2020       34846.17 31994.35 37697.99 30484.69 39207.65
## May 2020       34846.17 31657.74 38034.60 29969.88 39722.46
## Jun 2020       34846.17 31353.42 38338.92 29504.47 40187.87

f2_naive <- c(NA, head(ventas, -1))

# Modelo 2. Promedio Móvil (k=3)
f2_pm <- pm(ventas, 3)

# Modelo 3. Promedio Móvil Ponderado (k=3, pesos 3/6, 2/6 y 1/6)
f2_pmp <- pmp(ventas, c(3,2,1)/6)

# Modelo 4. Suavizador Exponencial (alfa = 0.2)
f2_ses <- suav(ventas, 0.2)

# Modelo 5. ARIMA
arima2 <- auto.arima(ts2); arima2

## Series: ts2 
## ARIMA(1,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##          ar1     sar1     drift
##       0.6383  -0.5517  288.8979
## s.e.  0.1551   0.2047   14.5026
## 
## sigma^2 = 202701:  log likelihood = -181.5
## AIC=371   AICc=373.11   BIC=375.72

f2_arima <- as.numeric(fitted(arima2))

pronóstico2 <- forecast(arima2, level=c(95), h=6)
pronóstico2

##          Point Forecast    Lo 95    Hi 95
## Jan 2020       35498.90 34616.48 36381.32
## Feb 2020       34202.17 33155.28 35249.05
## Mar 2020       36703.01 35596.10 37809.92
## Apr 2020       36271.90 35141.44 37402.36
## May 2020       37121.98 35982.07 38261.90
## Jun 2020       37102.65 35958.90 38246.40

plot(pronóstico2, ylab="Ventas en miles de USD", xlab="Año",
     main="Pronóstico a 6 meses con ARIMA")

modelos2 <- list("1. Naive"=f2_naive, "2. PM k=3"=f2_pm, "3. PMP k=3"=f2_pmp,
                 "4. SES alfa=0.2"=f2_ses, "5. ARIMA"=f2_arima)
res2 <- t(sapply(modelos2, métricas, real=ventas))
knitr::kable(res2, digits=4, caption="Tabla resumen de resultados, Ejercicio 2")

Tabla resumen de resultados, Ejercicio 2

	n	ME	MAE	MSE	RMSE	MPE	MAPE
1. Naive	35	266.4474	902.8514	1237975.6	1112.6435	0.8162	3.0131
2. PM k=3	33	591.0093	870.4659	1063495.6	1031.2592	1.9011	2.8117
3. PMP k=3	33	481.5529	814.1926	919646.3	958.9819	1.5425	2.6245
4. SES alfa=0.2	35	1255.3971	1402.2860	2458027.1	1567.8096	3.9842	4.5734
5. ARIMA	36	25.2216	227.1700	118242.4	343.8640	0.0806	0.7070

# Promedios anuales, para ver la tendencia
tapply(ventas, rep(2017:2019, each=12), mean)

##     2017     2018     2019 
## 26838.29 30567.67 33837.78

Conclusión del Ejercicio 2

El mejor es el ARIMA, con el MAPE más bajo por mucho, y el pronóstico de los seis meses que pide el ejercicio es el que sale en la tabla de arriba.

Gana porque es el único que toma en cuenta que la serie va subiendo. Los promedios de cada año van de 26,838 en 2017 a 30,568 en 2018 y 33,838 en 2019, y con ese crecimiento los modelos que nada más promedian lo que ya pasó siempre se quedan cortos. Se nota en que los cuatro modelos simples tienen el ME positivo, o sea que pronostican por debajo del valor real mes tras mes.

Y aquí el orden entre los modelos simples se voltea respecto al Ejercicio 1. El suavizador con alfa de 0.2 pasa de ser el mejor a ser el peor, y el ponderado le gana al simple. Cuando la serie trae dirección el dato más reciente sí sirve para saber qué sigue, así que suavizar de más estorba en lugar de ayudar.

---

Ejercicio 3. Vintage Restaurant

Karen Payne quiere saber cómo se comportan las ventas del restaurante y qué esperar del cuarto año. Los datos son las ventas mensuales de alimentos y bebidas en miles de dólares de los tres primeros años.

año1 <- c(242,235,232,178,184,140,145,152,110,130,152,206)
año2 <- c(263,238,247,193,193,149,157,161,122,130,167,230)
año3 <- c(282,255,265,205,210,160,166,174,126,148,173,235)

vintage <- ts(c(año1, año2, año3), start=c(1,1), frequency=12)
meses <- c("Enero","Febrero","Marzo","Abril","Mayo","Junio",
           "Julio","Agosto","Septiembre","Octubre","Noviembre","Diciembre")

knitr::kable(data.frame(Mes=meses, `Año 1`=año1, `Año 2`=año2, `Año 3`=año3,
                        check.names=FALSE),
             caption="Ventas de alimentos y bebidas del restaurante Vintage, miles de dólares")

Ventas de alimentos y bebidas del restaurante Vintage, miles de dólares

Mes	Año 1	Año 2	Año 3
Enero	242	263	282
Febrero	235	238	255
Marzo	232	247	265
Abril	178	193	205
Mayo	184	193	210
Junio	140	149	160
Julio	145	157	166
Agosto	152	161	174
Septiembre	110	122	126
Octubre	130	130	148
Noviembre	152	167	173
Diciembre	206	230	235

1. Gráfica de la serie de tiempo

plot(vintage, ylab="Ventas en miles de USD", xlab="Año de operación",
     main="Ventas mensuales del Vintage Restaurant", lwd=2, col="#2f5183")
points(vintage, pch=19, cex=0.5, col="#2f5183")
grid()

En la gráfica se ven dos cosas al mismo tiempo. La primera es un patrón que se repite igual cada año, con las ventas más altas en enero, febrero y marzo y las más bajas en septiembre. La segunda es que ese mismo patrón va subiendo de nivel año con año, o sea que además del ciclo hay tendencia.

Para Captiva Island tiene todo el sentido, porque es la temporada alta de turismo en Florida cuando en el norte hace frío, y septiembre cae en plena temporada de huracanes, que es cuando se vacía la isla.

2. Índices estacionales

ie <- índices_estacionales(vintage)
knitr::kable(data.frame(Mes=meses, `Índice estacional`=round(as.numeric(ie),4),
                        check.names=FALSE),
             caption="Índices estacionales por el método de razón al promedio móvil centrado")

Índices estacionales por el método de razón al promedio móvil centrado

Mes	Índice estacional
Enero	1.4436
Febrero	1.2997
Marzo	1.3441
Abril	1.0412
Mayo	1.0494
Junio	0.8004
Julio	0.8283
Agosto	0.8530
Septiembre	0.6280
Octubre	0.7003
Noviembre	0.8528
Diciembre	1.1593

cat("Los índices suman:", sum(ie), "\n")

## Los índices suman: 12

barplot(as.numeric(ie), names.arg=substr(meses,1,3), las=2,
        col=ifelse(as.numeric(ie)>1, "#1a7f37", "#c8502f"),
        main="Índice estacional por mes", ylab="Índice")
abline(h=1, lwd=2, lty=2)

Los índices sí tienen sentido intuitivo y bastante. Enero sale en 1.4436, o sea que un enero vende 44% más que un mes promedio del año, y marzo y febrero andan cerca con 1.3441 y 1.2997. Los tres meses de invierno son la temporada alta y ahí es donde el restaurante hace su dinero.

Del otro lado septiembre sale en 0.6280, que es vender 37% menos que un mes normal, y octubre en 0.7003. Ese es el fondo del año.

La diferencia entre el mejor y el peor mes es enorme, porque enero vende más del doble que septiembre con los mismos costos fijos encima. Eso ya es una recomendación por sí sola para Karen, los meses buenos tienen que cargar con los malos, y conviene planear el personal y el inventario mes por mes y no con un promedio anual.

3. Serie desestacionalizada

desestacionalizada <- as.numeric(vintage) / as.numeric(ie[cycle(vintage)])
ts_des <- ts(desestacionalizada, start=c(1,1), frequency=12)

plot(vintage, ylab="Ventas en miles de USD", xlab="Año de operación",
     main="Serie original contra serie desestacionalizada", col="#b8c4d4", lwd=2)
lines(ts_des, col="#c8502f", lwd=2)
legend("topright", c("Original","Desestacionalizada"),
       col=c("#b8c4d4","#c8502f"), lwd=2, bty="n")
grid()

t <- 1:36
tendencia <- lm(desestacionalizada ~ t)
summary(tendencia)

## 
## Call:
## lm(formula = desestacionalizada ~ t)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1892 -2.2527 -0.4848  0.8131  9.4189 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 169.3494     1.0926  155.00   <2e-16 ***
## t             1.0213     0.0515   19.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.21 on 34 degrees of freedom
## Multiple R-squared:  0.9204, Adjusted R-squared:  0.9181 
## F-statistic: 393.3 on 1 and 34 DF,  p-value: < 2.2e-16

Sí hay tendencia y es clara. Al quitarle el efecto de los meses, la serie queda en una línea que sube parejo, y la regresión da una pendiente de 1.0213 con un R² de 0.9204.

Eso quiere decir que el restaurante crece alrededor de mil dólares de venta cada mes, o sea unos 12,255 dólares más al año. El negocio está creciendo de verdad, no es que unos meses tapen a otros.

4. Pronóstico del año 4 por descomposición

t4 <- 37:48
tend4 <- coef(tendencia)[1] + coef(tendencia)[2]*t4
pron_descomposición <- tend4 * as.numeric(ie)

knitr::kable(data.frame(Mes=meses, Tendencia=round(tend4,2),
                        `Índice`=round(as.numeric(ie),4),
                        `Pronóstico`=round(pron_descomposición,2), check.names=FALSE),
             caption="Pronóstico del año 4 por descomposición, miles de dólares")

Pronóstico del año 4 por descomposición, miles de dólares

Mes	Tendencia	Índice	Pronóstico
Enero	207.14	1.4436	299.02
Febrero	208.16	1.2997	270.54
Marzo	209.18	1.3441	281.16
Abril	210.20	1.0412	218.86
Mayo	211.22	1.0494	221.65
Junio	212.24	0.8004	169.88
Julio	213.27	0.8283	176.65
Agosto	214.29	0.8530	182.78
Septiembre	215.31	0.6280	135.21
Octubre	216.33	0.7003	151.50
Noviembre	217.35	0.8528	185.35
Diciembre	218.37	1.1593	253.15

cat("Total del año 4 estimado:", round(sum(pron_descomposición),2), "miles de dólares\n")

## Total del año 4 estimado: 2545.76 miles de dólares

5. Pronóstico del año 4 con regresión de variables ficticias

mes_factor <- factor(rep(meses, 3), levels=meses)
regresión_dummies <- lm(as.numeric(vintage) ~ t + mes_factor)
summary(regresión_dummies)

## 
## Call:
## lm(formula = as.numeric(vintage) ~ t + mes_factor)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.1250 -1.9583  0.1667  2.2292  7.4583 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           249.10764    2.74567  90.727  < 2e-16 ***
## t                       1.01736    0.07554  13.467 2.14e-12 ***
## mes_factorFebrero     -20.68403    3.62688  -5.703 8.30e-06 ***
## mes_factorMarzo       -16.36806    3.62924  -4.510 0.000158 ***
## mes_factorAbril       -73.38542    3.63317 -20.199 3.90e-16 ***
## mes_factorMayo        -70.73611    3.63866 -19.440 8.97e-16 ***
## mes_factorJunio      -117.75347    3.64571 -32.299  < 2e-16 ***
## mes_factorJulio      -112.43750    3.65431 -30.768  < 2e-16 ***
## mes_factorAgosto     -107.12153    3.66445 -29.233  < 2e-16 ***
## mes_factorSeptiembre -151.13889    3.67611 -41.114  < 2e-16 ***
## mes_factorOctubre    -135.48958    3.68928 -36.725  < 2e-16 ***
## mes_factorNoviembre  -108.50694    3.70395 -29.295  < 2e-16 ***
## mes_factorDiciembre   -49.85764    3.72009 -13.402 2.36e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.441 on 23 degrees of freedom
## Multiple R-squared:  0.9942, Adjusted R-squared:  0.9911 
## F-statistic:   327 on 12 and 23 DF,  p-value: < 2.2e-16

nuevos <- data.frame(t=t4, mes_factor=factor(meses, levels=meses))
pron_dummies <- as.numeric(predict(regresión_dummies, nuevos))

knitr::kable(data.frame(Mes=meses,
                        `Descomposición`=round(pron_descomposición,2),
                        `Dummies`=round(pron_dummies,2),
                        `Diferencia`=round(pron_dummies-pron_descomposición,2),
                        check.names=FALSE),
             caption="Los dos pronósticos del año 4 comparados")

Los dos pronósticos del año 4 comparados

Mes	Descomposición	Dummies	Diferencia
Enero	299.02	286.75	-12.27
Febrero	270.54	267.08	-3.46
Marzo	281.16	272.42	-8.75
Abril	218.86	216.42	-2.45
Mayo	221.65	220.08	-1.57
Junio	169.88	174.08	4.21
Julio	176.65	180.42	3.77
Agosto	182.78	186.75	3.97
Septiembre	135.21	143.75	8.54
Octubre	151.50	160.42	8.92
Noviembre	185.35	188.42	3.07
Diciembre	253.15	248.08	-5.07

cat("Total por descomposición:", round(sum(pron_descomposición),2), "\n")

## Total por descomposición: 2545.76

cat("Total por dummies:      ", round(sum(pron_dummies),2), "\n")

## Total por dummies:       2544.67

6. Resumen, error de pronóstico y conclusiones

plot(ts(c(as.numeric(vintage), pron_descomposición), start=c(1,1), frequency=12),
     col="#2f5183", lwd=2, ylab="Ventas en miles de USD", xlab="Año de operación",
     main="Historia y pronóstico del año 4")
lines(ts(c(rep(NA,36), pron_dummies), start=c(1,1), frequency=12), col="#d9a021", lwd=2)
abline(v=4, lty=2)
legend("topright", c("Histórico y descomposición","Dummies"),
       col=c("#2f5183","#d9a021"), lwd=2, bty="n")
grid()

real_enero <- 295
e_desc <- real_enero - pron_descomposición[1]
e_dumm <- real_enero - pron_dummies[1]

knitr::kable(data.frame(
  Método = c("Descomposición","Regresión con dummies"),
  Pronóstico = round(c(pron_descomposición[1], pron_dummies[1]),2),
  Real = real_enero,
  Error = round(c(e_desc, e_dumm),2),
  `Error %` = round(c(e_desc, e_dumm)/real_enero*100, 2),
  check.names=FALSE),
  caption="Error de pronóstico de enero del año 4")

Error de pronóstico de enero del año 4

Método	Pronóstico	Real	Error	Error %
Descomposición	299.02	295	-4.02	-1.36
Regresión con dummies	286.75	295	8.25	2.80

cat("R² de la regresión con dummies:", round(summary(regresión_dummies)$r.squared,4), "\n")

## R² de la regresión con dummies: 0.9942

cat("R² de la tendencia desestacionalizada:", round(summary(tendencia)$r.squared,4), "\n")

## R² de la tendencia desestacionalizada: 0.9204

Informe para Karen

Las ventas del Vintage traen dos cosas encima al mismo tiempo, un patrón de temporada muy marcado y una tendencia de crecimiento. El patrón es el que salta a la vista, porque enero vende 44% arriba del mes promedio y septiembre 37% abajo, y esa diferencia se repite igualita los tres años. La tendencia está tapada por ese patrón y solo se ve al desestacionalizar, y ahí queda una línea limpia que sube 1,021 dólares cada mes, o sea como 12 mil dólares más de venta cada año.

Los dos métodos pronostican casi lo mismo para el año 4. Por descomposición sale 2,545.76 mil dólares en el año y por regresión con dummies 2,544.67, o sea que se separan por mil dólares sobre dos millones y medio. Donde sí difieren es mes por mes.

Sobre enero, el pronóstico por descomposición fue de 299.02 mil y las ventas reales salieron en 295 mil, así que se pasó por 4 mil dólares, un 1.4%. El de dummies dio 286.75 mil y se quedó corto por 8.25 mil, un 2.8%. Los dos le atinaron bien, y si a Karen le preocupa esa diferencia lo que hay que decirle es que 4 mil dólares sobre 295 mil es ruido normal de un mes y no una falla del método.

Lo que sí le recomendaría para no quedarse con la duda mes a mes es dejar de ver el pronóstico como un número solo. Tres cosas concretas. Que se maneje un rango en vez de una cifra, porque el error típico del modelo se puede calcular y le da un piso y un techo para planear. Que se recalculen los índices estacionales cada año conforme entren datos nuevos, ya que con tres años apenas van tres observaciones por mes y con cinco o seis la estimación se vuelve más firme. Y que se lleve registro del error mes con mes, porque si los errores empiezan a salir todos del mismo lado quiere decir que el negocio cambió de nivel y el modelo se tiene que reajustar.

---

Ejercicio 4. Carlson Department Store

Carlson estuvo cerrada de septiembre a diciembre del año 5 por el huracán. Hay que estimar cuánto habría vendido de no haber pasado, y si le corresponde algo del aumento de actividad comercial que hubo en la zona después.

# Ventas de Carlson en los 48 meses previos al huracán, de septiembre del año 1
# a agosto del año 5, en millones de dólares
carlson <- c(1.71,1.90,2.74,4.20,
             1.45,1.80,2.03,1.99,2.32,2.20,2.13,2.43,1.90,2.13,2.56,4.16,
             2.31,1.89,2.02,2.23,2.39,2.14,2.27,2.21,1.89,2.29,2.83,4.04,
             2.31,1.99,2.42,2.45,2.57,2.42,2.40,2.50,2.09,2.54,2.97,4.35,
             2.56,2.28,2.69,2.48,2.73,2.37,2.31,2.23)

# Ventas de todas las tiendas departamentales del condado, mismos 48 meses
condado <- c(55.80,56.40,71.40,117.60,
             46.80,48.00,60.00,57.60,61.80,58.20,56.40,63.00,57.60,53.40,71.40,114.00,
             46.80,48.60,59.40,58.20,60.60,55.20,51.00,58.80,49.80,54.60,65.40,102.00,
             43.80,45.60,57.60,53.40,56.40,52.80,54.00,60.60,47.40,54.60,67.80,100.20,
             48.00,51.60,57.60,58.20,60.00,57.00,57.60,61.80)

# Ventas reales del condado en los cuatro meses del huracán
condado_real <- c(69.00, 75.00, 85.20, 121.80)

ts_car <- ts(carlson, start=c(1,9), frequency=12)
ts_con <- ts(condado, start=c(1,9), frequency=12)

cat("Meses de Carlson:", length(carlson), "  Meses del condado:", length(condado), "\n")

## Meses de Carlson: 48   Meses del condado: 48

par(mfrow=c(2,1), mar=c(4,4,3,1))
plot(ts_car, ylab="Millones de USD", xlab="Año", col="#2f5183", lwd=2,
     main="Carlson Department Store, 48 meses previos al huracán"); grid()
plot(ts_con, ylab="Millones de USD", xlab="Año", col="#c8502f", lwd=2,
     main="Tiendas departamentales del condado, mismos 48 meses"); grid()

par(mfrow=c(1,1))

1. Ventas que Carlson habría registrado sin el huracán

ie_car <- índices_estacionales(ts_car)
des_car <- as.numeric(ts_car) / as.numeric(ie_car[cycle(ts_car)])
tc <- 1:48
tend_car <- lm(des_car ~ tc)

knitr::kable(data.frame(Mes=meses, `Índice Carlson`=round(as.numeric(ie_car),4),
                        check.names=FALSE),
             caption="Índices estacionales de Carlson")

Índices estacionales de Carlson

Mes	Índice Carlson
Enero	0.9566
Febrero	0.8194
Marzo	0.9071
Abril	0.9294
Mayo	1.0113
Junio	0.9372
Julio	0.9357
Agosto	0.9744
Septiembre	0.7967
Octubre	0.9357
Noviembre	1.1191
Diciembre	1.6774

summary(tend_car)

## 
## Call:
## lm(formula = des_car ~ tc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69031 -0.06049  0.01174  0.07180  0.32586 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.149084   0.048946   43.91  < 2e-16 ***
## tc          0.011407   0.001739    6.56 4.18e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1669 on 46 degrees of freedom
## Multiple R-squared:  0.4833, Adjusted R-squared:  0.4721 
## F-statistic: 43.03 on 1 and 46 DF,  p-value: 4.181e-08

t5 <- 49:52   # septiembre a diciembre del año 5
est_car <- (coef(tend_car)[1] + coef(tend_car)[2]*t5) * as.numeric(ie_car[9:12])

knitr::kable(data.frame(Mes=meses[9:12],
                        `Ventas estimadas sin huracán`=round(est_car,3),
                        check.names=FALSE),
             caption="Carlson, septiembre a diciembre del año 5, millones de dólares")

Carlson, septiembre a diciembre del año 5, millones de dólares

Mes	Ventas estimadas sin huracán
Septiembre	2.158
Octubre	2.544
Noviembre	3.056
Diciembre	4.600

cat("Total estimado sin huracán:", round(sum(est_car),3), "millones\n")

## Total estimado sin huracán: 12.358 millones

2. Ventas que el condado habría tenido sin el huracán

ie_con <- índices_estacionales(ts_con)
des_con <- as.numeric(ts_con) / as.numeric(ie_con[cycle(ts_con)])
tend_con <- lm(des_con ~ tc)
summary(tend_con)

## 
## Call:
## lm(formula = des_con ~ tc)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.520 -2.402 -0.015  1.484  5.598 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 62.22039    0.74999   82.96  < 2e-16 ***
## tc          -0.07194    0.02665   -2.70  0.00967 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.558 on 46 degrees of freedom
## Multiple R-squared:  0.1368, Adjusted R-squared:  0.118 
## F-statistic: 7.289 on 1 and 46 DF,  p-value: 0.009675

est_con <- (coef(tend_con)[1] + coef(tend_con)[2]*t5) * as.numeric(ie_con[9:12])

knitr::kable(data.frame(Mes=meses[9:12],
                        `Estimado sin huracán`=round(est_con,3),
                        `Real con huracán`=condado_real,
                        `Exceso`=round(condado_real-est_con,3),
                        `Exceso %`=round((condado_real/est_con-1)*100,2),
                        check.names=FALSE),
             caption="Condado, septiembre a diciembre del año 5, millones de dólares")

Condado, septiembre a diciembre del año 5, millones de dólares

Mes	Estimado sin huracán	Real con huracán	Exceso	Exceso %
Septiembre	50.549	69.0	18.451	36.50
Octubre	53.195	75.0	21.805	40.99
Noviembre	66.783	85.2	18.417	27.58
Diciembre	103.115	121.8	18.685	18.12

cat("Total estimado sin huracán:", round(sum(est_con),3), "millones\n")

## Total estimado sin huracán: 273.642 millones

cat("Total real con huracán:    ", round(sum(condado_real),3), "millones\n")

## Total real con huracán:     351 millones

3. Pérdida de Carlson y el argumento del exceso

factor_exceso <- sum(condado_real) / sum(est_con)
pérdida_base  <- sum(est_car)
pérdida_ajust <- pérdida_base * factor_exceso

knitr::kable(data.frame(
  Concepto = c("Ventas normales estimadas de Carlson",
               "Factor de exceso del condado",
               "Ventas estimadas con el efecto del huracán",
               "Pérdida total reclamable"),
  Valor = c(round(pérdida_base,3), round(factor_exceso,4),
            round(pérdida_ajust,3), round(pérdida_ajust,3))),
  caption="Estimación de la pérdida de Carlson, millones de dólares")

Estimación de la pérdida de Carlson, millones de dólares

Concepto	Valor
Ventas normales estimadas de Carlson	12.3580
Factor de exceso del condado	1.2827
Ventas estimadas con el efecto del huracán	15.8510
Pérdida total reclamable	15.8510

cat("El condado vendió", round((factor_exceso-1)*100,2), "% más de lo esperado\n")

## El condado vendió 28.27 % más de lo esperado

cat("Pérdida sin ajustar:", round(pérdida_base,3), "millones\n")

## Pérdida sin ajustar: 12.358 millones

cat("Pérdida con el ajuste por huracán:", round(pérdida_ajust,3), "millones\n")

## Pérdida con el ajuste por huracán: 15.851 millones

cat("Diferencia entre las dos cifras:", round(pérdida_ajust-pérdida_base,3), "millones\n")

## Diferencia entre las dos cifras: 3.494 millones

Informe para los directivos de Carlson

Cuánto se perdió. Con el patrón de los 48 meses previos, Carlson habría vendido alrededor de 12.36 millones de dólares entre septiembre y diciembre del año 5. Esa es la pérdida de ventas por los cuatro meses que estuvo cerrada, calculada solo con su propia historia y sin meter nada del huracán.

Diciembre es el mes que más pesa de los cuatro, con 4.6 millones estimados, porque su índice estacional es de 1.6774 y vende casi 68% arriba de un mes normal. Cerrar en diciembre es lo que más caro salió.

El argumento del exceso, y sí se sostiene. El condado debió haber vendido cerca de 273.64 millones en esos cuatro meses según su propia tendencia, y vendió 351. Son 77.36 millones de más, un 28.27% arriba de lo normal, y eso no se explica con el comportamiento histórico del condado sino con los más de 8 mil millones de ayuda federal y pagos de seguros que entraron a la zona.

El exceso no está parejo. Octubre fue el mes más fuerte con 41% arriba de lo esperado y septiembre el más flojo con 36%, mientras que diciembre apenas subió 18%. Se ve como gente reponiendo lo que perdió en los meses inmediatos al huracán, y no como un cambio permanente del mercado.

Si Carlson hubiera estado abierta, no hay razón para pensar que se habría quedado fuera de ese movimiento, porque el dinero se repartió entre las tiendas departamentales del condado y Carlson es una de ellas. Aplicando el mismo factor de 1.2827, las ventas que habría hecho suben de 12.36 a 15.85 millones, o sea 3.49 millones más de reclamación.

Lo que le recomendaría a la tienda. Reclamar los 12.36 millones como pérdida base, que esa es la parte sólida y sale de su propia historia, y presentar los 3.49 millones adicionales por separado, con el análisis del condado como respaldo. Conviene presentarlos aparte y no revueltos, porque la primera cifra es difícil de discutir y la segunda depende de aceptar que Carlson habría capturado el mismo porcentaje de aumento que las demás tiendas.

Ahí está el punto débil que hay que reconocer antes de que lo saquen del otro lado. El supuesto es que Carlson habría crecido igual que el promedio del condado, y eso no está demostrado, solo es razonable. Además la regresión de tendencia del condado tiene un R² bajo, así que su estimación es menos firme que la de Carlson. Vale la pena decirlo de una vez y ofrecer la comparación mes por mes, que es donde el argumento se ve más convincente.