Act 2 - GEFA

Instalar paquetes y llamar libreria

library(forecast)

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

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

#file.choose()
poblacion <- read.csv("/Users/jenaromtzg/Desktop/Generación de Escenarios Futuros/R/Semana 2/population.csv")

summary(poblacion)

##     state                year        population      
##  Length:6020        Min.   :1900   Min.   :   43000  
##  Class :character   1st Qu.:1930   1st Qu.:  901483  
##  Mode  :character   Median :1960   Median : 2359000  
##                     Mean   :1960   Mean   : 3726003  
##                     3rd Qu.:1990   3rd Qu.: 4541883  
##                     Max.   :2019   Max.   :39512223

str(poblacion)

## 'data.frame':    6020 obs. of  3 variables:
##  $ state     : chr  "AK" "AK" "AK" "AK" ...
##  $ year      : int  1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 ...
##  $ population: int  135000 158000 189000 205000 215000 222000 224000 231000 224000 224000 ...

head(poblacion)

##   state year population
## 1    AK 1950     135000
## 2    AK 1951     158000
## 3    AK 1952     189000
## 4    AK 1953     205000
## 5    AK 1954     215000
## 6    AK 1955     222000

poblacion_texas <- poblacion %>% filter(state=="TX")
ggplot(poblacion_texas, aes(x=year, y=population)) + geom_line() +
  labs(title="Población de Texas", x = "Año", 
       y ="Población")

ts_texas <- ts(poblacion_texas$population, start = 1900, frequency=1)

# Serie de tiempo anual
## ts_texas <- ts(poblacion_texas$population start = c(1900, 4), frequency=4) - TRIMESTRAL
## ts_texas <- ts(poblacion_texas$population start = c(1900, 8), frequency=12) - MENSUAL


arima_texas <- auto.arima(ts_texas)
summary(arima_texas)

## Series: ts_texas 
## ARIMA(0,2,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.5950  -0.1798
## s.e.   0.0913   0.0951
## 
## sigma^2 = 1.031e+10:  log likelihood = -1527.14
## AIC=3060.28   AICc=3060.5   BIC=3068.6
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE      MAPE      MASE
## Training set 12147.62 99818.31 59257.39 0.1046163 0.5686743 0.2672197
##                     ACF1
## Training set -0.02136734

pronostico_texas <- forecast(arima_texas, level=95, h=50)
pronostico_texas

##      Point Forecast    Lo 95    Hi 95
## 2020       29398472 29199487 29597457
## 2021       29806827 29463665 30149990
## 2022       30215183 29742956 30687410
## 2023       30623538 30024100 31222977
## 2024       31031894 30303359 31760429
## 2025       31440249 30579246 32301253
## 2026       31848605 30851090 32846119
## 2027       32256960 31118581 33395339
## 2028       32665316 31381587 33949044
## 2029       33073671 31640070 34507272
## 2030       33482027 31894047 35070007
## 2031       33890382 32143561 35637204
## 2032       34298738 32388674 36208801
## 2033       34707093 32629456 36784730
## 2034       35115449 32865983 37364914
## 2035       35523804 33098330 37949278
## 2036       35932160 33326573 38537746
## 2037       36340515 33550788 39130242
## 2038       36748871 33771046 39726695
## 2039       37157226 33987418 40327034
## 2040       37565581 34199972 40931191
## 2041       37973937 34408774 41539100
## 2042       38382292 34613887 42150698
## 2043       38790648 34815371 42765925
## 2044       39199003 35013284 43384723
## 2045       39607359 35207682 44007036
## 2046       40015714 35398618 44632810
## 2047       40424070 35586145 45261995
## 2048       40832425 35770311 45894540
## 2049       41240781 35951163 46530399
## 2050       41649136 36128748 47169524
## 2051       42057492 36303110 47811874
## 2052       42465847 36474290 48457405
## 2053       42874203 36642330 49106076
## 2054       43282558 36807269 49757848
## 2055       43690914 36969145 50412683
## 2056       44099269 37127994 51070544
## 2057       44507625 37283853 51731396
## 2058       44915980 37436755 52395205
## 2059       45324336 37586734 53061937
## 2060       45732691 37733822 53731560
## 2061       46141047 37878050 54404044
## 2062       46549402 38019447 55079357
## 2063       46957758 38158044 55757471
## 2064       47366113 38293868 56438358
## 2065       47774469 38426948 57121989
## 2066       48182824 38557310 57808338
## 2067       48591180 38684979 58497380
## 2068       48999535 38809982 59189088
## 2069       49407891 38932343 59883438

plot(pronostico_texas, main="Poblacion en Texas")

Ejercicio en Clase Lunes 17: Mapa

library(forecast)
library(tidyverse)
library(ggplot2)
library(maps)

## 
## Attaching package: 'maps'

## The following object is masked from 'package:purrr':
## 
##     map

#file.choose()
poblacion <- read.csv("/Users/jenaromtzg/Desktop/Generación de Escenarios Futuros/R/Semana 2/population.csv")

summary(poblacion)

##     state                year        population      
##  Length:6020        Min.   :1900   Min.   :   43000  
##  Class :character   1st Qu.:1930   1st Qu.:  901483  
##  Mode  :character   Median :1960   Median : 2359000  
##                     Mean   :1960   Mean   : 3726003  
##                     3rd Qu.:1990   3rd Qu.: 4541883  
##                     Max.   :2019   Max.   :39512223

str(poblacion)

## 'data.frame':    6020 obs. of  3 variables:
##  $ state     : chr  "AK" "AK" "AK" "AK" ...
##  $ year      : int  1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 ...
##  $ population: int  135000 158000 189000 205000 215000 222000 224000 231000 224000 224000 ...

head(poblacion)

##   state year population
## 1    AK 1950     135000
## 2    AK 1951     158000
## 3    AK 1952     189000
## 4    AK 1953     205000
## 5    AK 1954     215000
## 6    AK 1955     222000

# Crear un mapa de EUA por década, con un gradiente verde-rojo de la poblacion por estado, desde 1900, hasta el 2100
map(database="state")
map(database="state", regions = "New York", col = "green", fill = TRUE, add = TRUE)
map(database="state", regions = "Texas", col = "red", fill = TRUE, add = TRUE)
title(main = "Pronóstico de Población en EE.UU.")

Generar Pronostico por Cada Estado

# Suponiendo que en 'poblacion' tienes columnas:
# state (abreviatura, ej. "TX"), year (año), population (valor)
poblacion_extended <- poblacion %>%
  group_by(state) %>%
  arrange(year) %>%
  # Para cada estado, creamos un df con valores observados + pronóstico
  do({
    df <- .
    # Años mínimos y máximos en tus datos
    min_year <- min(df$year)
    max_year <- max(df$year)
    
    # Creamos la serie de tiempo anual
    ts_pop <- ts(df$population, start = min_year, frequency = 1)
    
    # Determinamos cuántos años faltan para llegar a 2050
    horizon <- 2050 - max_year
    
    # Si el dataset llega hasta antes de 2050, hacemos forecast
    if(horizon > 0){
      fit <- auto.arima(ts_pop)
      fc <- forecast(fit, h = horizon)
      
      # Data frame con los datos pronosticados
      years_forecast <- (max_year + 1):2050
      df_forecast <- data.frame(
        state      = unique(df$state),
        year       = years_forecast,
        population = as.numeric(fc$mean)
      )
      
      # Unimos histórico + forecast
      df_all <- bind_rows(
        # Histórico (columnas relevantes)
        df %>% select(state, year, population),
        # Futuro
        df_forecast
      )
    } else {
      # Si ya tenemos datos hasta 2050 o más, no pronosticamos
      df_all <- df
    }
    
    df_all
  }) %>%
  ungroup()

Convertir abreviaturas de estado (p.ej., “TX”) a nombres completos

df_state_names <- data.frame(
  state_abb  = state.abb,
  state_name = tolower(state.name),
  stringsAsFactors = FALSE
)

poblacion_full <- poblacion_extended %>%
  left_join(df_state_names, by = c("state" = "state_abb"))

Unir la geometría de los estados con los datos de población

# 4.1) Obtenemos las coordenadas de polígonos con 'map_data("state")'
states_map <- map_data("state")

# 4.2) Hacemos left_join para heredar la columna 'population' (por estado y año)
map_data_joined <- states_map %>%
  left_join(poblacion_full, by = c("region" = "state_name"))

## Warning in left_join(., poblacion_full, by = c(region = "state_name")): Each row in `x` is expected to match at most 1 row in `y`.
## ℹ Row 1 of `x` matches multiple rows.
## ℹ If multiple matches are expected, set `multiple = "all"` to silence this
##   warning.

Filtrar décadas entre 1950 y 2050

map_data_decadas <- map_data_joined %>%
  mutate(year = as.numeric(year)) %>%
  filter(year >= 1950, year <= 2050, year %% 10 == 0)

Graficar con ggplot2 (rojo = población baja, verde = población alta)Filtrar décadas entre 1950 y 2050

ggplot(map_data_decadas, aes(x = long, 
                             y = lat, 
                             group = group, 
                             fill = population)) +
  geom_polygon(color = "black", size = 0.1) +
  # Gradiente de rojo (baja) a verde (alta)
  scale_fill_gradient(low = "red", high = "green", 
                      na.value = "grey90") +
  facet_wrap(~ year) +
  coord_fixed(1.3) +
  labs(title = "Población de EUA por Estado (1950 - 2050, por década)",
       fill = "Población estimada") +
  theme_void() +
  theme(legend.position = "right",
        strip.text = element_text(face = "bold"))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Actividad 2. Leche saborizada Hershey’s

library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

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

library(forecast)
library(tidyverse)
library(ggplot2)

#file.choose()
ventas <- read.csv("/Users/jenaromtzg/Desktop/Generación de Escenarios Futuros/R/Semana 2/Ventas_Históricas_Lechitas.csv")

Modelo Arima

ts_ventas <- ts(ventas$Ventas,start=c(2017,1), frequency = 12)
autoplot(ts_ventas) + labs(title = "Ventas de leche saborizada Hershey's", x = "Tiempo" , y = "Miles de Dólares")

arima_ventas <- auto.arima(ts_ventas)
summary(arima_ventas)

## Series: ts_ventas 
## ARIMA(1,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##          ar1     sar1     drift
##       0.6383  -0.5517  288.8980
## s.e.  0.1551   0.2047   14.5026
## 
## sigma^2 = 202700:  log likelihood = -181.5
## AIC=371   AICc=373.11   BIC=375.72
## 
## Training set error measures:
##                    ME    RMSE      MAE        MPE      MAPE       MASE
## Training set 25.22163 343.863 227.1699 0.08059942 0.7069541 0.06491041
##                   ACF1
## Training set 0.2081043

pronostico_ventas <- forecast(arima_ventas, level=95, h = 12)
pronostico_ventas

##          Point Forecast    Lo 95    Hi 95
## Jan 2020       35498.90 34616.48 36381.32
## Feb 2020       34202.17 33155.29 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.91 38246.40
## Jul 2020       37151.04 36005.74 38296.35
## Aug 2020       38564.65 37418.71 39710.59
## Sep 2020       38755.23 37609.03 39901.42
## Oct 2020       39779.03 38632.73 40925.33
## Nov 2020       38741.63 37595.29 39887.97
## Dec 2020       38645.86 37499.50 39792.22

autoplot(pronostico_ventas) + labs(title = "Pronostico de ventas 2020 de lecehe saborizada Hershey's", x="Tiempo", y = "Miles de Dólares")

Modelo de regresión Lineal

ventas$mes <- 1:36
regresion_ventas <- lm(Ventas ~ mes, data = ventas)
summary(regresion_ventas)

## 
## Call:
## lm(formula = Ventas ~ mes, data = ventas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2075.79  -326.41    33.74   458.41  1537.04 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 24894.67     275.03   90.52   <2e-16 ***
## mes           298.37      12.96   23.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 808 on 34 degrees of freedom
## Multiple R-squared:  0.9397, Adjusted R-squared:  0.9379 
## F-statistic: 529.8 on 1 and 34 DF,  p-value: < 2.2e-16

siguiente_anio <- data.frame(mes=37:48)
prediccion_regresion <- predict(regresion_ventas, siguiente_anio)
prediccion_regresion

##        1        2        3        4        5        6        7        8 
## 35934.49 36232.86 36531.23 36829.61 37127.98 37426.35 37724.73 38023.10 
##        9       10       11       12 
## 38321.47 38619.85 38918.22 39216.59

plot(ventas$mes, ventas$Ventas, main = "Pronostico de Ventas 20220 de leche saborizadas Hershey's", xlab = "Tiempo", ylab = "Miles de Dólares")
abline(regresion_ventas, col = "blue")
points(siguiente_anio$mes, prediccion_regresion, col = "red")

prediccion_reales <- predict(regresion_ventas, ventas)

MAPE <- mean(abs((ventas$Ventas - prediccion_reales)/ventas$Ventas))*100
MAPE

## [1] 2.011298

## El mejor modelo que se adapta a la serie es el SARIMA con un MAPE DE 0.71%, comparado con la de regresion lineal que su MAPE es de 2.01%. Para el siguiente año, la proyección de ventas es la siguiente:

tabla <- data.frame(
  `Mes y Año` = c("Jan 2020", "Feb 2020", "Mar 2020", "Apr 2020", "May 2020",
                  "Jun 2020", "Jul 2020", "Aug 2020", "Sep 2020", "Oct 2020"),
  Esperado = c(35498.90, 34202.17, 36703.01, 36271.90, 37121.98, 
               37102.65, 37151.04, 38564.64, 38755.22, 39779.02),
  Pesimista = c(34616.48, 33155.28, 35596.10, 35141.44, 35982.07, 
                35958.90, 36005.73, 37418.70, 37609.03, 38632.72),
  Optimista = c(36381.32, 35249.05, 37809.92, 37402.36, 38261.90, 
                38246.40, 38296.34, 39710.58, 39901.42, 40925.32)
)

# Generar tabla con kable
tabla %>%
  kable("html", caption = "Proyección Financiera 2020") %>%
  kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover"))

Proyección Financiera 2020

Mes.y.Año	Esperado	Pesimista	Optimista
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
Jul 2020	37151.04	36005.73	38296.34
Aug 2020	38564.64	37418.70	39710.58
Sep 2020	38755.22	37609.03	39901.42
Oct 2020	39779.02	38632.72	40925.32

#file.choose()
ventas_por_anio <- read.csv("/Users/jenaromtzg/Desktop/Generación de Escenarios Futuros/R/Semana 2/ventas_por_anio.csv")

ggplot(ventas_por_anio, aes(x=mes, y=ventas, col=as.factor(anio),group=anio)) + geom_line() + 
  labs(title = "Ventas de leche saborizada Hershey's por año", x="Ventas",y="Miles de Dólares")

## Nuestra recomendación seria realizar campañas publicitarias para aumentar el consumo de leches saborizadas Hershey’s en el primer semestre del año.