Actividad 2

Ejercicio en clase: Población

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.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

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

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

Importar la base de datos

#file.choose()
poblacion <- read.csv("/Users/danielemilianonajeraotero/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="Población de Texas", x ="Año", y="Población")

ts_texas <- ts(poblacion_texas$population, start = 1950, frequency=1) #Serie de tiempo anual
#ts_texas <- ts(poblacion_texas$population, start = c(1900, 4), frequency=4) #Serie de tiempo trimestral
#ts_texas <- ts(poblacion_texas$population, start = c(1900, 8), frequency=) #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
## 2070       29398472 29199487 29597457
## 2071       29806827 29463665 30149990
## 2072       30215183 29742956 30687410
## 2073       30623538 30024100 31222977
## 2074       31031894 30303359 31760429
## 2075       31440249 30579246 32301253
## 2076       31848605 30851090 32846119
## 2077       32256960 31118581 33395339
## 2078       32665316 31381587 33949044
## 2079       33073671 31640070 34507272

plot(pronostico_texas, main = "Población en Texas")

Crear un Mapa

## Instalar y cargar las librerías necesarias para análisis y visualización


#install.packages(c("dplyr", "ggplot2", "forecast", "maps"))
library(dplyr)
library(ggplot2)
library(forecast)
library(maps)

## 
## Attaching package: 'maps'

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

## 1) Importar la base de datos de población 

poblacion <- read.csv("/Users/danielemilianonajeraotero/Downloads/population.csv")

## 2) Proyectar la población estatal hasta 2050 y combinar con datos originales

# Integración de pronósticos con la base de datos inicial
# Crear un dataset ampliado que parte de la base de datos original
poblacion_extendida <- poblacion

# Extraer la lista de estados únicos presentes en los datos
estados <- unique(poblacion$state)

# Modelar y pronosticar la población de cada estado hasta el año 2050
for(st in estados){
  
  # Filtrar y organizar los datos por año para cada estado
  datos_st <- poblacion %>%
    filter(state == st) %>%
    arrange(year)
  
  # Identificar el último año disponible en la serie histórica
  ultimo_anio <- max(datos_st$year)
  
  # Generar la serie de tiempo anual para la población del estado
  ts_st <- ts(datos_st$population,
              start = min(datos_st$year),
              end   = ultimo_anio,
              frequency = 1)  # Datos anuales
  
  # Ajustar un modelo ARIMA de manera automática para el estado
  modelo_st <- auto.arima(ts_st)
  
  # Calcular el horizonte de pronóstico necesario
  # (se genera forecast solo si faltan años para llegar a 2050)
  h_years <- 2050 - ultimo_anio
  
  if(h_years > 0){
    # Generar el pronóstico para el periodo faltante
    pronostico <- forecast(modelo_st, h = h_years)
    
    # Crear un dataframe con las proyecciones generadas
    anios_pronostico <- (ultimo_anio + 1):2050
    poblacion_pronosticada <- as.numeric(pronostico$mean)
    
    df_forecast <- data.frame(
      state = st,
      year  = anios_pronostico,
      population = poblacion_pronosticada
    )
    
    # Incorporar las proyecciones al dataset extendido
    poblacion_extendida <- rbind(poblacion_extendida, df_forecast)
  }
}

## 3) Definir una función para visualizar el mapa de población por año 

plot_map <- function(year) {
  
  # Filtrar el dataset para el año especificado
  data_year <- poblacion_extendida %>%
    filter(year == !!year)
  
  # Cargar la información geográfica de los estados de EE.UU.
  states_map <- map_data("state")
  
  # Relacionar las abreviaturas estatales con los nombres completos en minúsculas
  # Utilizando los vectores auxiliares state.abb y state.name
  data_year <- data_year %>%
    mutate(region = tolower(state.name[match(state, state.abb)])) %>%
    right_join(states_map, by = "region")
  
  # Generar el mapa temático por población
  ggplot(data_year, aes(x = long, y = lat, group = group, fill = population)) +
    geom_polygon(color = "black") +
    # Aplicar un gradiente de color verde (bajo) a rojo (alto)
    scale_fill_gradient(
      low = "green",   # Representa la menor población
      high = "red",    # Indica la mayor población
      name = "Población"
    ) +
    labs(
      title = paste("Población por Estado en", year)
    ) +
    theme_void() +
    theme(
      legend.position = "right",
      plot.title = element_text(size = 16, face = "bold")
    )
}

## 4) Visualizar la evolución demográfica cada década (1950 - 2050) -----------

for(year in seq(1950, 2050, by = 10)) {
  print(plot_map(year))
}

Actividad 2. Leche saborizada Hershey’s

##Instalar paquetes y llamar librerias

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

##Importar base de datos

#file.choose()
ventas <- read.csv("/Users/danielemilianonajeraotero/Downloads/Ventas_Históricas_Lechitas.csv")

##1. Modelo AUTO.ARIMA

str(ventas)

## 'data.frame':    12 obs. of  6 variables:
##  $ Mes     : chr  "ene-17" "feb-17" "mar-17" "abr-17" ...
##  $ Ventas  : chr  " 25,520.51 " " 23,740.11 " " 26,253.58 " " 25,868.43 " ...
##  $ Mes.1   : chr  "ene-18" "feb-18" "mar-18" "abr-18" ...
##  $ Ventas.1: chr  " 28,463.69 " " 26,996.11 " " 29,768.20 " " 29,292.51 " ...
##  $ Mes.2   : chr  "ene-19" "feb-19" "mar-19" "abr-19" ...
##  $ Ventas.2: chr  " 32,496.44 " " 31,287.28 " " 33,376.02 " " 32,949.77 " ...

head(ventas)

##      Mes      Ventas  Mes.1    Ventas.1  Mes.2    Ventas.2
## 1 ene-17  25,520.51  ene-18  28,463.69  ene-19  32,496.44 
## 2 feb-17  23,740.11  feb-18  26,996.11  feb-19  31,287.28 
## 3 mar-17  26,253.58  mar-18  29,768.20  mar-19  33,376.02 
## 4 abr-17  25,868.43  abr-18  29,292.51  abr-19  32,949.77 
## 5 may-17  27,072.87  may-18  29,950.68  may-19  34,004.11 
## 6 jun-17  27,150.50  jun-18  30,099.17  jun-19  33,757.89

ventas$Ventas <- gsub(",", "", ventas$Ventas)  # Elimina comas
ventas$Ventas <- as.numeric(as.character(ventas$Ventas))

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,1,0) with drift 
## 
## Coefficients:
##           ar1     drift
##       -0.8190  292.0190
## s.e.   0.1727  122.4236
## 
## sigma^2 = 608020:  log likelihood = -88.31
## AIC=182.62   AICc=186.05   BIC=183.81
## 
## Training set error measures:
##                     ME    RMSE     MAE        MPE     MAPE MASE      ACF1
## Training set -65.80225 675.289 575.958 -0.2638862 2.179855  NaN 0.1402595

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

##          Point Forecast    Lo 95    Hi 95
## Jan 2018       28047.46 26519.16 29575.75
## Feb 2018       28418.95 26865.82 29972.07
## Mar 2018       28645.88 26619.36 30672.39
## Apr 2018       28991.21 26912.66 31069.75
## May 2018       29239.57 26864.20 31614.93
## Jun 2018       29567.34 27120.62 32014.07
## Jul 2018       29830.07 27168.41 32491.74
## Aug 2018       30146.08 27401.36 32890.80
## Sep 2018       30418.45 27504.20 33332.71
## Oct 2018       30726.56 27723.22 33729.91
## Nov 2018       31005.40 27860.29 34150.52
## Dec 2018       31308.22 28071.72 34544.71

autoplot (pronostico_ventas) + labs (title= "pronóstico de Ventas 2020 de Leche Saborizada Hershey's", x="Tiempo", y = "Miles de Dólares")

##2. Modelo de regresión Lineal

nrow(ventas)

## [1] 12

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

## 
## Call:
## lm(formula = Ventas ~ mes, data = ventas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1752.6  -324.3   203.1   475.3   858.4 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  24894.7      479.9  51.878 1.71e-13 ***
## mes            299.0       65.2   4.586    0.001 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 779.7 on 10 degrees of freedom
## Multiple R-squared:  0.6777, Adjusted R-squared:  0.6455 
## F-statistic: 21.03 on 1 and 10 DF,  p-value: 0.001001

siguiente_anio <- data.frame(mes=37.48)
prediccion_regresion <- predict(regresion_ventas, siguiente_anio)
prediccion_regresion

##        1 
## 36101.85

plot(ventas$mes, ventas$Ventas, main="Pronostico de ventas 2020 de leche saborizada Hershey", xlab= "Tiempo", ylab= "Miles de Dolares")
abline(regresion_ventas, col="blue")
points(siguiente_anio$mes,  prediccion_regresion, col = "red")

predicciones_reales <- predict(regresion_ventas, ventas)
MAPE <- mean(abs((ventas$Ventas - 
                    predicciones_reales)/ventas$Ventas))*100
MAPE

## [1] 2.147841

Conclusiones

El mejor modelo que se adapta a la serie es el SARIMA con un MAPE de 0.71%, comparada con la regresión 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 | Esnearios 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 |

#file.choose()
ventas_por_anio <- read.csv("/Users/danielemilianonajeraotero/Downloads/ventas_por_anio (2).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 por Año",
       x="Mes", y="Miles de dolares")

Nuestra recoemdación seria realizar campañas para aumentar el consumo de leche saborizada Hersheys en el primer semestre del año.