Series de Tiempo

Instalar paquetes y llamar librerías

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

Ejercicio 1. Ventas Semanales

INSTRUCCIONES: Modela la sigueinte Serie de Tiempo y elige la mejor opción. Pronostica las sigueintes 6 semanas.

semana <- c(1:12)
valor <- c(17,21,19,23,18,16,20,18,22,20,15,22)
df1 <- data.frame(semana,valor)

ts1 <- ts(valor, c(2025,1), frequency = 52)

# 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

mape_naive1 <- accuracy(naive1)[1,"MAPE"]
# Modelo con menor MAPE es el mas adecuado.

# Función genérica de Promedio Móvil: promedia los últimos k valores para
# obtener el ajuste dentro de la muestra, y de forma recursiva usa los
# pronósticos ya calculados para proyectar h periodos hacia adelante.
ma_forecast <- function(serie, k, h){
  n <- length(serie)
  fitted <- rep(NA, n)
  for (i in (k+1):n) fitted[i] <- mean(serie[(i-k):(i-1)])
  mape <- mean(abs((serie[(k+1):n] - fitted[(k+1):n]) / serie[(k+1):n])) * 100
  ext <- c(serie, rep(NA, h))
  for (i in (n+1):(n+h)) ext[i] <- mean(ext[(i-k):(i-1)])
  list(fitted = fitted, mape = mape, forecast = ext[(n+1):(n+h)])
}

# Modelo 2. Promedio Móvil (k = 3)
ma1 <- ma_forecast(valor, k = 3, h = 6)
ma1$mape       # MAPE del modelo

## [1] 14.35661

ma1$forecast   # Pronóstico de las siguientes 6 semanas

## [1] 19.00000 18.66667 19.88889 19.18519 19.24691 19.44033

# Función genérica de Promedio Móvil Ponderado: los pesos se aplican del más
# antiguo al más reciente, por lo que el valor más reciente pesa más (3/6).
wma_forecast <- function(serie, pesos, h){
  k <- length(pesos)
  n <- length(serie)
  fitted <- rep(NA, n)
  for (i in (k+1):n) fitted[i] <- sum(serie[(i-k):(i-1)] * pesos)
  mape <- mean(abs((serie[(k+1):n] - fitted[(k+1):n]) / serie[(k+1):n])) * 100
  ext <- c(serie, rep(NA, h))
  for (i in (n+1):(n+h)) ext[i] <- sum(ext[(i-k):(i-1)] * pesos)
  list(fitted = fitted, mape = mape, forecast = ext[(n+1):(n+h)])
}

pesos <- c(1,2,3)/6  # 1/6 al más antiguo, 2/6 al intermedio, 3/6 al más reciente

# Modelo 3. Promedio Móvil Ponderado (k =3 Pesos: 3/6, 2/6, y 1/6)
wma1 <- wma_forecast(valor, pesos, h = 6)
wma1$mape

## [1] 15.99248

wma1$forecast

## [1] 19.33333 19.50000 19.86111 19.65278 19.69676 19.70949

# Modelo 4. Suavizador Exponencial (alfa = 0.2). Menor alfa = mejor modelo
ses1 <- ses(ts1, alpha = 0.2, h = 6)
summary(ses1)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
## ses(y = ts1, h = 6, alpha = 0.2)
## 
##   Smoothing parameters:
##     alpha = 0.2 
## 
##   Initial states:
##     l = 19.1774 
## 
##   sigma:  2.9273
## 
##      AIC     AICc      BIC 
## 57.40932 58.74266 58.37914 
## 
## Error measures:
##                      ME     RMSE      MAE       MPE     MAPE MASE       ACF1
## Training set 0.06549665 2.672288 2.312859 -1.479963 12.37329  NaN -0.3669082
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2025.231       19.33458 15.58304 23.08613 13.59709 25.07208
## 2025.250       19.33458 15.50875 23.16042 13.48347 25.18570
## 2025.269       19.33458 15.43587 23.23330 13.37201 25.29716
## 2025.288       19.33458 15.36432 23.30485 13.26259 25.40657
## 2025.308       19.33458 15.29405 23.37512 13.15512 25.51405
## 2025.327       19.33458 15.22497 23.44419 13.04948 25.61969

mape_ses1 <- accuracy(ses1)[1,"MAPE"]

# 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

pronostico1 <- forecast(arima1, level = c(95), h = 6)
plot(pronostico1)

mape_arima1 <- accuracy(arima1)[1,"MAPE"]

# Tabla con el Resumen de los Resultados.
resumen1 <- data.frame(
  Modelo = c("Naive", "Promedio Móvil (k=3)", "Promedio Móvil Ponderado",
             "Suavizador Exponencial (alfa=0.2)", "ARIMA"),
  MAPE = round(c(mape_naive1, ma1$mape, wma1$mape, mape_ses1, mape_arima1), 2)
)
resumen1 <- resumen1[order(resumen1$MAPE), ]
resumen1

##                              Modelo  MAPE
## 5                             ARIMA 11.19
## 4 Suavizador Exponencial (alfa=0.2) 12.37
## 2              Promedio Móvil (k=3) 14.36
## 3          Promedio Móvil Ponderado 15.99
## 1                             Naive 19.24

pronostico_final1 <- data.frame(
  Semana = 13:18,
  Naive = as.numeric(naive1$mean),
  Prom_Movil = ma1$forecast,
  Prom_Movil_Pond = wma1$forecast,
  Suavizador_Exp = as.numeric(ses1$mean),
  ARIMA = as.numeric(pronostico1$mean)
)
pronostico_final1

##   Semana Naive Prom_Movil Prom_Movil_Pond Suavizador_Exp ARIMA
## 1     13    22   19.00000        19.33333       19.33458 19.25
## 2     14    22   18.66667        19.50000       19.33458 19.25
## 3     15    22   19.88889        19.86111       19.33458 19.25
## 4     16    22   19.18519        19.65278       19.33458 19.25
## 5     17    22   19.24691        19.69676       19.33458 19.25
## 6     18    22   19.44033        19.70949       19.33458 19.25

# Conclusión Método ARE.
mejor1 <- resumen1$Modelo[1]
cat("El modelo con el MAPE más bajo (", round(resumen1$MAPE[1],2),
    "%) es:", mejor1,
    ". Por lo tanto es el modelo más adecuado para pronosticar las siguientes 6 semanas.\n")

## El modelo con el MAPE más bajo ( 11.19 %) es: ARIMA . Por lo tanto es el modelo más adecuado para pronosticar las siguientes 6 semanas.

Ejercicio 2. Hershey’s

INSTRUCCIONES: Modela la sigueinte Serie de Tiempo y elige la mejor opción. Pronostica las sigueintes 6 semanas.

ruta <- "/Users/adriana/Downloads/Ventas_Históricas_Lechitas.xlsx"

datos <- read_excel(ruta, sheet = 1, skip = 2)

## New names:
## • `Mes` -> `Mes...1`
## • `Ventas` -> `Ventas...2`
## • `Mes` -> `Mes...3`
## • `Ventas` -> `Ventas...4`
## • `Mes` -> `Mes...5`
## • `Ventas` -> `Ventas...6`

ventas <- as.numeric(c(datos[[2]], datos[[4]], datos[[6]]))  # concatena 2017-2018-2019
mes <- 1:length(ventas)
df2 <- data.frame(mes, ventas)

ts2 <- ts(ventas, start = c(2017,1), frequency = 12)
plot(ts2, main = "Ventas históricas de leche saborizada Hershey (Lechitas)",
     ylab = "Ventas (miles de dólares)", xlab = "Año")

# 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

mape_naive2 <- accuracy(naive2)[1,"MAPE"]
# Modelo con menor MAPE es el mas adecuado.

# Modelo 2. Promedio Móvil (k = 3)
ma2 <- ma_forecast(ventas, k = 3, h = 6)
ma2$mape

## [1] 2.811654

ma2$forecast

## [1] 35259.72 34968.60 35024.83 35084.38 35025.94 35045.05

# Modelo 3. Promedio Móvil Ponderado (k =3 Pesos: 3/6, 2/6, y 1/6)
wma2 <- wma_forecast(ventas, pesos, h = 6)
wma2$mape

## [1] 2.624516

wma2$forecast

## [1] 35045.23 34937.99 34958.44 34966.09 34958.85 34961.20

# Modelo 4. Suavizador Exponencial (alfa = 0.2). Menor alfa = mejor modelo
ses2 <- ses(ts2, alpha = 0.2, h = 6)
summary(ses2)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
## ses(y = ts2, h = 6, alpha = 0.2)
## 
##   Smoothing parameters:
##     alpha = 0.2 
## 
##   Initial states:
##     l = 26303.9678 
## 
##   sigma:  1574.86
## 
##      AIC     AICc      BIC 
## 661.0074 661.3710 664.1744 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 1111.747 1530.489 1335.458 3.458378 4.357771 0.3815873 0.3247295
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2020       34308.55 32290.28 36326.81 31221.88 37395.22
## Feb 2020       34308.55 32250.31 36366.78 31160.75 37456.34
## Mar 2020       34308.55 32211.10 36405.99 31100.78 37516.31
## Apr 2020       34308.55 32172.61 36444.48 31041.92 37575.17
## May 2020       34308.55 32134.81 36482.28 30984.10 37632.99
## Jun 2020       34308.55 32097.65 36519.44 30927.27 37689.82

mape_ses2 <- accuracy(ses2)[1,"MAPE"]

# 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

pronostico2 <- forecast(arima2, level = c(95), h = 6)
plot(pronostico2)

mape_arima2 <- accuracy(arima2)[1,"MAPE"]

# Tabla con el Resumen de los Resultados.
resumen2 <- data.frame(
  Modelo = c("Naive", "Promedio Móvil (k=3)", "Promedio Móvil Ponderado",
             "Suavizador Exponencial (alfa=0.2)", "ARIMA"),
  MAPE = round(c(mape_naive2, ma2$mape, wma2$mape, mape_ses2, mape_arima2), 2)
)
resumen2 <- resumen2[order(resumen2$MAPE), ]
resumen2

##                              Modelo MAPE
## 5                             ARIMA 0.71
## 3          Promedio Móvil Ponderado 2.62
## 2              Promedio Móvil (k=3) 2.81
## 1                             Naive 3.01
## 4 Suavizador Exponencial (alfa=0.2) 4.36

pronostico_final2 <- data.frame(
  Mes = 37:42,
  Naive = as.numeric(naive2$mean),
  Prom_Movil = ma2$forecast,
  Prom_Movil_Pond = wma2$forecast,
  Suavizador_Exp = as.numeric(ses2$mean),
  ARIMA = as.numeric(pronostico2$mean)
)
pronostico_final2

##   Mes    Naive Prom_Movil Prom_Movil_Pond Suavizador_Exp    ARIMA
## 1  37 34846.17   35259.72        35045.23       34308.55 35498.90
## 2  38 34846.17   34968.60        34937.99       34308.55 34202.17
## 3  39 34846.17   35024.83        34958.44       34308.55 36703.01
## 4  40 34846.17   35084.38        34966.09       34308.55 36271.90
## 5  41 34846.17   35025.94        34958.85       34308.55 37121.98
## 6  42 34846.17   35045.05        34961.20       34308.55 37102.65

# Conclusión Método ARE.
mejor2 <- resumen2$Modelo[1]
cat("El modelo con el MAPE más bajo (", round(resumen2$MAPE[1],2),
    "%) es:", mejor2,
    ". Por lo tanto es el modelo más adecuado para pronosticar las ventas de los próximos 6 meses de Lechitas.\n")

## El modelo con el MAPE más bajo ( 0.71 %) es: ARIMA . Por lo tanto es el modelo más adecuado para pronosticar las ventas de los próximos 6 meses de Lechitas.

Ejercicio 3. Vintage Restaurant

INSTRUCCIONES: Analiza las ventas de alimentos y bebidas del Vintage Restaurant (Tabla 18.26) e incluye: gráfica de la serie, índices estacionales, serie desestacionalizada, pronóstico por descomposición y por regresión con variables ficticias para el año 4, tablas resumen y conclusiones.

# Ventas de alimentos y bebidas del restaurante Vintage ($ miles)
vintage <- c(242,235,232,178,184,140,145,152,110,130,152,206,   # Año 1
             263,238,247,193,193,149,157,161,122,130,167,230,   # Año 2
             282,255,265,205,210,160,166,174,126,148,173,235)   # Año 3

ts_vintage <- ts(vintage, start = c(1,1), frequency = 12)

# 1. Gráfica de la serie de tiempo
plot(ts_vintage, main = "Ventas de alimentos y bebidas - Vintage Restaurant",
     ylab = "Ventas ($ miles)", xlab = "Año", type = "o")

Patrón principal: se observa un componente estacional muy marcado que se repite cada 12 meses (ventas altas en la temporada de invierno: dic-ene-feb-mar, y bajas en el verano/otoño, con el mínimo en septiembre), junto con una tendencia creciente de un año a otro conforme crece la reputación del restaurante.

# 2. Descomposición multiplicativa e índices estacionales
decomp_vintage <- decompose(ts_vintage, type = "multiplicative")
plot(decomp_vintage)

indices_estacionales <- data.frame(
  Mes = month.name,
  Indice = round(as.numeric(decomp_vintage$figure), 3)
)
indices_estacionales

##          Mes Indice
## 1    January  1.444
## 2   February  1.300
## 3      March  1.344
## 4      April  1.041
## 5        May  1.049
## 6       June  0.800
## 7       July  0.828
## 8     August  0.853
## 9  September  0.628
## 10   October  0.700
## 11  November  0.853
## 12  December  1.159

Los índices estacionales muestran temporada alta en enero-marzo y diciembre (índices arriba de 1, ventas por encima del promedio) y temporada baja en septiembre y octubre (índices muy por debajo de 1). Esto tiene sentido intuitivo: Captiva Island es un destino turístico de invierno (temporada de “snowbirds”), mientras que septiembre-octubre coincide con la temporada baja de turismo y de huracanes en Florida.

# 3. Serie desestacionalizada y análisis de tendencia
deseasonalized <- ts_vintage / decomp_vintage$seasonal
plot(deseasonalized, main = "Serie desestacionalizada - Vintage Restaurant",
     ylab = "Ventas ($ miles)", xlab = "Año", type = "o")

t <- 1:36
modelo_tendencia <- lm(as.numeric(deseasonalized) ~ t)
summary(modelo_tendencia)

## 
## Call:
## lm(formula = as.numeric(deseasonalized) ~ 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

La serie desestacionalizada sí muestra una tendencia lineal creciente: la pendiente estimada es positiva y estadísticamente significativa, lo que confirma que, quitando el efecto estacional, las ventas del restaurante crecen mes con mes.

# 4. Pronóstico por el método de descomposición para el año 4
pred_tendencia_a4 <- predict(modelo_tendencia, newdata = data.frame(t = 37:48))
estacional_a4 <- as.numeric(decomp_vintage$figure)  # 12 índices, Ene-Dic
pronostico_descomp <- pred_tendencia_a4 * estacional_a4

pronostico_descomp_tabla <- data.frame(
  Mes = month.name,
  Pronostico_Descomposicion = round(pronostico_descomp, 1)
)
pronostico_descomp_tabla

##          Mes Pronostico_Descomposicion
## 1    January                     299.0
## 2   February                     270.5
## 3      March                     281.2
## 4      April                     218.9
## 5        May                     221.7
## 6       June                     169.9
## 7       July                     176.6
## 8     August                     182.8
## 9  September                     135.2
## 10   October                     151.5
## 11  November                     185.3
## 12  December                     253.2

# 5. Pronóstico por regresión con variables ficticias (dummy variables)
mes_factor <- factor(rep(1:12, 3))
modelo_dummy <- lm(vintage ~ t + mes_factor)
summary(modelo_dummy)

## 
## Call:
## lm(formula = 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_factor2   -20.68403    3.62688  -5.703 8.30e-06 ***
## mes_factor3   -16.36806    3.62924  -4.510 0.000158 ***
## mes_factor4   -73.38542    3.63317 -20.199 3.90e-16 ***
## mes_factor5   -70.73611    3.63866 -19.440 8.97e-16 ***
## mes_factor6  -117.75347    3.64571 -32.299  < 2e-16 ***
## mes_factor7  -112.43750    3.65431 -30.768  < 2e-16 ***
## mes_factor8  -107.12153    3.66445 -29.233  < 2e-16 ***
## mes_factor9  -151.13889    3.67611 -41.114  < 2e-16 ***
## mes_factor10 -135.48958    3.68928 -36.725  < 2e-16 ***
## mes_factor11 -108.50694    3.70395 -29.295  < 2e-16 ***
## mes_factor12  -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

pronostico_regresion <- predict(
  modelo_dummy,
  newdata = data.frame(t = 37:48, mes_factor = factor(1:12, levels = levels(mes_factor)))
)

pronostico_regresion_tabla <- data.frame(
  Mes = month.name,
  Pronostico_Regresion = round(as.numeric(pronostico_regresion), 1)
)
pronostico_regresion_tabla

##          Mes Pronostico_Regresion
## 1    January                286.8
## 2   February                267.1
## 3      March                272.4
## 4      April                216.4
## 5        May                220.1
## 6       June                174.1
## 7       July                180.4
## 8     August                186.8
## 9  September                143.8
## 10   October                160.4
## 11  November                188.4
## 12  December                248.1

# 6. Tabla resumen con ambos métodos y gráfica comparativa
resumen_vintage <- data.frame(
  Mes = month.name,
  Descomposicion = round(pronostico_descomp, 1),
  Regresion_Dummy = round(as.numeric(pronostico_regresion), 1)
)
resumen_vintage

##          Mes Descomposicion Regresion_Dummy
## 1    January          299.0           286.8
## 2   February          270.5           267.1
## 3      March          281.2           272.4
## 4      April          218.9           216.4
## 5        May          221.7           220.1
## 6       June          169.9           174.1
## 7       July          176.6           180.4
## 8     August          182.8           186.8
## 9  September          135.2           143.8
## 10   October          151.5           160.4
## 11  November          185.3           188.4
## 12  December          253.2           248.1

matplot(1:12, resumen_vintage[, c("Descomposicion","Regresion_Dummy")], type = "o", pch = 16,
        col = c("blue","red"), xaxt = "n",
        main = "Pronóstico Año 4 - Vintage Restaurant",
        xlab = "Mes", ylab = "Ventas ($ miles)")
axis(1, at = 1:12, labels = month.abb)
legend("topright", legend = c("Descomposición","Regresión Dummy"), col = c("blue","red"), lty = 1, pch = 16)

# Error de pronóstico de enero del año 4 (venta real reportada = $295 miles)
venta_real_enero_a4 <- 295
error_descomp_ene <- venta_real_enero_a4 - resumen_vintage$Descomposicion[1]
error_regresion_ene <- venta_real_enero_a4 - resumen_vintage$Regresion_Dummy[1]

cat("Error de pronóstico (Descomposición) para enero del año 4:", round(error_descomp_ene,1), "\n")

## Error de pronóstico (Descomposición) para enero del año 4: -4

cat("Error de pronóstico (Regresión con dummies) para enero del año 4:", round(error_regresion_ene,1), "\n")

## Error de pronóstico (Regresión con dummies) para enero del año 4: 8.2

Conclusiones Ejercicio 3: ambos métodos arrojan pronósticos similares para enero del año 4, y el error frente a la venta real ($295 mil) es relativamente pequeño en comparación con el nivel de ventas, lo que indica que los modelos capturan bien el patrón estacional y la tendencia. Si a Karen le preocupa esa diferencia, puede reducir la incertidumbre generando un intervalo de pronóstico (en vez de un solo número), actualizando el modelo cada mes con los datos reales más recientes (pronóstico móvil / recalibración continua), y monitoreando si el error se mantiene dentro de un rango razonable o si crece de forma sistemática, lo cual señalaría que el patrón de la serie está cambiando.

Ejercicio 4. Carlson Department Store

INSTRUCCIONES: Estima las ventas que Carlson y el condado habrían tenido de septiembre a diciembre del año 5 si no hubiera ocurrido el huracán, estima la pérdida de ventas de Carlson, y evalúa si hay evidencia de un exceso de ventas en el condado relacionado con el huracán.

# Ventas de Carlson Department Store ($ millones)
# 48 meses: de septiembre del Año 1 a agosto del Año 5 (mes previo al huracán)
carlson <- c(1.71,1.90,2.74,4.20,                                          # Sep-Dic Año 1
             1.45,1.80,2.03,1.99,2.32,2.20,2.13,2.43,1.90,2.13,2.56,4.16,  # Año 2 completo
             2.31,1.89,2.02,2.23,2.39,2.14,2.27,2.21,1.89,2.29,2.83,4.04,  # Año 3 completo
             2.31,1.99,2.42,2.45,2.57,2.42,2.40,2.50,2.09,2.54,2.97,4.35,  # Año 4 completo
             2.56,2.28,2.69,2.48,2.73,2.37,2.31,2.23)                     # Ene-Ago Año 5

# Ventas de las tiendas departamentales del condado ($ millones)
# Mismos 48 meses que la serie de Carlson
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 durante los 4 meses en que Carlson permaneció cerrada
# (Sep-Dic Año 5, después del huracán)
condado_real_cierre <- c(69.00, 75.00, 85.20, 121.80)

length(carlson)  # debe ser 48

## [1] 48

length(condado)  # debe ser 48

## [1] 48

# Secuencia de mes calendario (1=Ene ... 12=Dic) para los 48 meses, empezando en septiembre
meses_seq <- c(9,10,11,12, rep(1:12, 3), 1,2,3,4,5,6,7,8)
factor_meses <- factor(meses_seq)
t48 <- 1:48

# 1. Estimación de las ventas que Carlson habría tenido de Sep a Dic del Año 5
modelo_carlson <- lm(carlson ~ t48 + factor_meses)
summary(modelo_carlson)

## 
## Call:
## lm(formula = carlson ~ t48 + factor_meses)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50750 -0.06604  0.00875  0.07458  0.28750 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.901944   0.090906  20.922  < 2e-16 ***
## t48             0.011111   0.001753   6.338 2.78e-07 ***
## factor_meses2  -0.178611   0.115236  -1.550   0.1301    
## factor_meses3   0.110278   0.115276   0.957   0.3453    
## factor_meses4   0.096667   0.115342   0.838   0.4077    
## factor_meses5   0.300556   0.115436   2.604   0.0134 *  
## factor_meses6   0.069444   0.115555   0.601   0.5517    
## factor_meses7   0.053333   0.115701   0.461   0.6477    
## factor_meses8   0.107222   0.115874   0.925   0.3611    
## factor_meses9  -0.215556   0.115436  -1.867   0.0702 .  
## factor_meses10  0.090833   0.115342   0.788   0.4363    
## factor_meses11  0.639722   0.115276   5.549 3.03e-06 ***
## factor_meses12  2.041111   0.115236  17.712  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1629 on 35 degrees of freedom
## Multiple R-squared:  0.9472, Adjusted R-squared:  0.9291 
## F-statistic: 52.35 on 12 and 35 DF,  p-value: < 2.2e-16

pronostico_carlson_cierre <- predict(
  modelo_carlson,
  newdata = data.frame(t48 = 49:52, factor_meses = factor(c(9,10,11,12), levels = levels(factor_meses)))
)
pronostico_carlson_cierre

##        1        2        3        4 
## 2.230833 2.548333 3.108333 4.520833

# 2. Estimación de las ventas que el condado habría tenido de Sep a Dic del Año 5
modelo_condado <- lm(condado ~ t48 + factor_meses)
summary(modelo_condado)

## 
## Call:
## lm(formula = condado ~ t48 + factor_meses)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5925 -2.0250  0.0475  1.4963  7.4925 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    48.46792    1.80898  26.793  < 2e-16 ***
## t48            -0.09208    0.03488  -2.640 0.012311 *  
## factor_meses2   2.19208    2.29314   0.956 0.345663    
## factor_meses3  12.48417    2.29393   5.442 4.20e-06 ***
## factor_meses4  10.77625    2.29526   4.695 4.02e-05 ***
## factor_meses5  13.71833    2.29711   5.972 8.41e-07 ***
## factor_meses6   9.91042    2.29950   4.310 0.000126 ***
## factor_meses7   8.95250    2.30241   3.888 0.000431 ***
## factor_meses8  15.34458    2.30584   6.655 1.07e-07 ***
## factor_meses9   5.93167    2.29711   2.582 0.014160 *  
## factor_meses10  8.12375    2.29526   3.539 0.001155 ** 
## factor_meses11 22.46583    2.29393   9.794 1.46e-11 ***
## factor_meses12 62.00792    2.29314  27.041  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.243 on 35 degrees of freedom
## Multiple R-squared:  0.9693, Adjusted R-squared:  0.9587 
## F-statistic: 92.02 on 12 and 35 DF,  p-value: < 2.2e-16

pronostico_condado_cierre <- predict(
  modelo_condado,
  newdata = data.frame(t48 = 49:52, factor_meses = factor(c(9,10,11,12), levels = levels(factor_meses)))
)
pronostico_condado_cierre

##        1        2        3        4 
##  49.8875  51.9875  66.2375 105.6875

# 3. Estimación de la pérdida de ventas de Carlson (Sep-Dic Año 5)
resumen_carlson <- data.frame(
  Mes = c("Septiembre","Octubre","Noviembre","Diciembre"),
  Ventas_Estimadas = round(as.numeric(pronostico_carlson_cierre), 3),
  Ventas_Reales = c(0, 0, 0, 0)  # la tienda estuvo cerrada, ventas reales = 0
)
resumen_carlson$Perdida <- resumen_carlson$Ventas_Estimadas - resumen_carlson$Ventas_Reales
resumen_carlson

##          Mes Ventas_Estimadas Ventas_Reales Perdida
## 1 Septiembre            2.231             0   2.231
## 2    Octubre            2.548             0   2.548
## 3  Noviembre            3.108             0   3.108
## 4  Diciembre            4.521             0   4.521

perdida_total_carlson <- sum(resumen_carlson$Perdida)
cat("Pérdida total estimada de ventas de Carlson (Sep-Dic Año 5): $",
    round(perdida_total_carlson,2), "millones\n")

## Pérdida total estimada de ventas de Carlson (Sep-Dic Año 5): $ 12.41 millones

# 4. Comparación de ventas reales del condado vs. estimadas 
resumen_condado <- data.frame(
  Mes = c("Septiembre","Octubre","Noviembre","Diciembre"),
  Ventas_Estimadas = round(as.numeric(pronostico_condado_cierre), 2),
  Ventas_Reales = condado_real_cierre
)
resumen_condado$Exceso <- resumen_condado$Ventas_Reales - resumen_condado$Ventas_Estimadas
resumen_condado

##          Mes Ventas_Estimadas Ventas_Reales Exceso
## 1 Septiembre            49.89          69.0  19.11
## 2    Octubre            51.99          75.0  23.01
## 3  Noviembre            66.24          85.2  18.96
## 4  Diciembre           105.69         121.8  16.11

exceso_total_condado <- sum(resumen_condado$Exceso)
cat("Exceso total de ventas reales del condado vs. lo estimado (Sep-Dic Año 5): $",
    round(exceso_total_condado,2), "millones\n")

## Exceso total de ventas reales del condado vs. lo estimado (Sep-Dic Año 5): $ 77.19 millones

Conclusiones Ejercicio 4:

1. Se estima que, de no haber ocurrido el huracán, Carlson habría vendido aproximadamente 12.41 millones de dólares entre septiembre y diciembre del año 5; como la tienda permaneció cerrada (ventas reales = 0), esa cifra representa la pérdida de ventas atribuible al huracán.
2. Las ventas reales del condado en esos mismos 4 meses superaron en 77.19 millones de dólares lo que el modelo estimaba que habrían sido sin el huracán.
3. Ese exceso de ventas en el condado (probablemente impulsado por los más de $8,000 millones en ayuda federal y pagos de seguros que se inyectaron a la zona) es evidencia a favor del argumento de Carlson de que hubo un incremento en la actividad comercial de la zona tras el huracán — lo cual respalda que, de haber permanecido abierta, la tienda muy probablemente también se habría beneficiado de ese repunte y no solo habría perdido sus ventas “normales”, sino también una parte de ese exceso de demanda.