ACT 1

Instalar paquetes y llamar librerias

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

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

#install.packages("tidyverse")
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4

## ── 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

Instalar base de datos

#file.choose()
poblacion <- read.csv("C:\\Users\\Luis\\Downloads\\population.csv")

Entender la base de datos

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

Serie de Tiempo en Texas

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

ts_texas <- ts(poblacion_texas$population, start=1900, frequency = 1) # Serie de Tiempo Anual
# ts_texas <- ts(poblacion_texas$population, start=1900, frequency = 4) # Serie de Tiempo Trimestral
# ts_texas <- ts(poblacion_texas$population, start=1900, frequency = 12) # Serie de Tiempo 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 = 10)
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

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

Ejercicio en Clase Lunes 17: Mapa

Crear un Mapa

# Crear un mapa de EUA por decada, con un gradiente verde-rojo de la poblacion por estado, desde 1950 hasta 2050
maps::map(database = "state")
maps::map(database = "state", regions = "Texas", col = "red", fill = TRUE, add = TRUE)
maps::map(database = "state", regions = "New York", col = "green", fill = TRUE, add = TRUE)

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")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 102 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

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.

#ACT. 2 LECHE

#install package ("tidyverse")
library(tidyverse)
#install package ("forecast)
library(forecast)
#install package ("ggplot2")
library(ggplot2)
#install.packages("readxl")
library(readxl)

Importar la base de datos

ventas= read_excel("C:\\Users\\Luis\\Downloads\\Ventas_Históricas_Lechitas.xlsx")

1.pronostico de ventas

ts_ventas = ts (ventas$Ventas, start = c(2017,1), frequency = 12)
autoplot(ts_ventas) + labs(title = "Ventas de leche saborizadas Hershey´s", x= "Tiempo", y="Miles dlls")

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 leche saborizada", x ="Tiempo", y="Miles de dlls")

2.Modelo regresion 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_ano = data.frame(mes=37:48)
prediccion_regresion = predict(regresion_ventas, siguiente_ano)
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 2020", xlab= "Tiempo", ylab = "Miles de dlls")
abline(regresion_ventas, col = "blue")
points(siguiente_ano$mes, prediccion_regresion, col="red")

predicciones_reales = predict(regresion_ventas, ventas)
predicciones_reales

##        1        2        3        4        5        6        7        8 
## 25193.04 25491.42 25789.79 26088.16 26386.54 26684.91 26983.28 27281.66 
##        9       10       11       12       13       14       15       16 
## 27580.03 27878.40 28176.78 28475.15 28773.52 29071.90 29370.27 29668.64 
##       17       18       19       20       21       22       23       24 
## 29967.02 30265.39 30563.76 30862.14 31160.51 31458.89 31757.26 32055.63 
##       25       26       27       28       29       30       31       32 
## 32354.01 32652.38 32950.75 33249.13 33547.50 33845.87 34144.25 34442.62 
##       33       34       35       36 
## 34740.99 35039.37 35337.74 35636.11

MAPE = mean(abs((ventas$Ventas - predicciones_reales)/ventas$Ventas))*100
MAPE

## [1] 2.011298

3.conclusiones

El mejor modelo que se adapta a la serie es el SARIMA cpn un MAPE de 0.71% comparado con la regersion lineal que su MAPE es de 2.01%

para el siguiente año, la proyección de ventas es la siguiente:
| Mes y Año | Escenario Esperado | Escenario Pesimista | Escenario 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 | | Nov 2020 | 38741.63 | 37595.28 | 39887.97 | | Dec 2020 | 38645.86 | 37499.50 | 39792.22 |

ventas_por_anio = read.csv("C:\\Users\\Luis\\Downloads\\ventas_por_anio (1).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="Mes", y="Miles de dlls")

Nuestra recomendacion seria realizar campañas publicitarias para aumentar el consumo de leche saborizada en el primer semestre del año.