1 Importamos las librerías

library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(plm)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plm':
## 
##     between, lag, lead
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(maps)

2 Cargar datos

library(readxl)
datos <- read_excel("~/Desktop/historical_state_population_by_year.xlsx")
colnames(datos)<-c("Estado","Año","Poblacion")

3 Serie de Tiempo

3.1 Construcción del Modelo

Realizaremos un loop para poder hacer la predicción de todos los estados

# Ordenar los datos por año
datos <- datos[order(datos$Año), ]
# Dividir los datos en una lista por estado
panel_list <- split(datos, datos$Estado)
# Inicializar una lista para almacenar los modelos ARIMA
arima_models <- list()

# Iterar sobre cada estado para ajustar un modelo ARIMA y hacer pronósticos
for (state in names(panel_list)) {
  state_data <- panel_list[[state]]
  ts_state <- ts(data = state_data$Poblacion, start = min(state_data$Año), frequency = 1)
  arima_model <- auto.arima(ts_state)
  forecast_values <- forecast(arima_model, h = 60)  # Hacemos pronósticos hasta el año 2070
  clean_state_name <- state
  arima_models[[clean_state_name]] <- list(model = arima_model, forecast = forecast_values)
  
  plot(forecast_values, main = paste("Pronóstico de población en", state))
}

3.2 Exportación de datos predichos

Generamos una data frame distinto para poder graficar los mapas

# Crear un dataframe vacío para almacenar los pronósticos en forma de panel
forecast_panel <- data.frame(state = character(), year = numeric(), prediction = numeric())

# Iterar sobre los nombres en panel_list
for (state in names(panel_list)) {
  # Acceder al pronóstico correspondiente en arima_models
  forecast <- arima_models[[state]]$forecast
  # Crear un dataframe temporal para el estado actual
  state_forecast <- data.frame(state = rep(state, length(forecast)), year = seq(2020, 2079), prediction = forecast)
  # Agregar el dataframe temporal al dataframe forecast_panel
  forecast_panel <- rbind(forecast_panel, state_forecast)
}

# Seleccionar y renombrar columnas pertinentes
forecast_panel <- forecast_panel %>% select(1:3)
colnames(forecast_panel) <- c("Estado", "Año", "Poblacion")

3.3 Construcción de Mapas

Nuevamente crearemos un loop para la creación de los mapas

# Definir años específicos para el análisis
anios <- c(2030, 2040, 2050, 2060, 2070)

# Iterar sobre los años para generar mapas
for (i in anios) {
  # Filtrar el dataframe de pronósticos para el año actual
  df2 <- forecast_panel %>% filter(Año == i)
  
  # Calcular los cuartiles y asignar colores correspondientes
  df2 <- within(df2, quartile <- as.integer(cut(Poblacion, quantile(Poblacion, probs = 0:4/4), include.lowest = TRUE)))
  df2$col[df2$quartile == 1] <- 'cadetblue2'
  df2$col[df2$quartile == 2] <- 'cadetblue3'
  df2$col[df2$quartile == 3] <- 'cadetblue'
  df2$col[df2$quartile == 4] <- 'cadetblue4'
  
  # Renombrar los estados con códigos a nombres completos
  df2$Estado[df2$Estado == 'TX'] <- 'Texas'
  df2$Estado[df2$Estado == 'ND'] <- 'North Dakota'
  df2$Estado[df2$Estado == 'SD'] <- 'South Dakota'
  df2$Estado[df2$Estado == 'KS'] <- 'Kansas'
  df2$Estado[df2$Estado == 'IA'] <- 'Iowa'
  df2$Estado[df2$Estado == 'VT'] <- 'Vermont'
  df2$Estado[df2$Estado == 'PA'] <- 'Pennsylvania'
  df2$Estado[df2$Estado == 'KY'] <- 'Kentucky'
  df2$Estado[df2$Estado == 'WV'] <- 'West Virginia'
  df2$Estado[df2$Estado == 'VA'] <- 'Virginia'
  df2$Estado[df2$Estado == 'TN'] <- 'Tennessee'
  df2$Estado[df2$Estado == 'NC'] <- 'North Carolina'
  df2$Estado[df2$Estado == 'SC'] <- 'South Carolina'
  df2$Estado[df2$Estado == 'GA'] <- 'Georgia'
  df2$Estado[df2$Estado == 'LA'] <- 'Louisiana'
  
  # Generar mapa con los colores asignados
  map(database = 'state')
  map(database = 'state', regions = df2$state, col = df2$col, fill = TRUE, add = TRUE)
  title(paste("Mapa del año", i))
}

LS0tCnRpdGxlOiAiQWN0aXZpZGFkIEludGVncmFkb3JhIgphdXRob3I6IEpvc8OpIEdhYnJpZWwgVXNpw7FhIE1vZ3JvIEEwMDgzMTQzNSAgTG9yZW5hIFZpbGxhcnJlYWwgVmVnYSBBMDE3MjA4MDIgWGltZW5hCiAgU29sw61zIElzbGFzIEEwMDgzMTM3MQpkYXRlOiAiMjAyNC0wMi0yMSIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0aGVtZTogdW5pdGVkCiAgICBoaWdobGlnaHQ6IHRhbmdvCiAgICB0b2M6IHRydWUKICAgIHRvY19kZXB0aDogMwogICAgbnVtYmVyX3NlY3Rpb25zOiBUUlVFCiAgICB0b2NfZmxvYXQ6IHRydWUKICAgIGNvZGVfZG93bmxvYWQ6IHRydWUKLS0tCgojIEltcG9ydGFtb3MgbGFzIGxpYnJlcsOtYXMKYGBge3Igd2FybmluZz1GQUxTRX0KbGlicmFyeShmb3JlY2FzdCkKbGlicmFyeShwbG0pCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkobWFwcykKYGBgCgojIENhcmdhciBkYXRvcwpgYGB7cn0KbGlicmFyeShyZWFkeGwpCmRhdG9zIDwtIHJlYWRfZXhjZWwoIn4vRGVza3RvcC9oaXN0b3JpY2FsX3N0YXRlX3BvcHVsYXRpb25fYnlfeWVhci54bHN4IikKY29sbmFtZXMoZGF0b3MpPC1jKCJFc3RhZG8iLCJBw7FvIiwiUG9ibGFjaW9uIikKYGBgCgoKIyBTZXJpZSBkZSBUaWVtcG8KIyMgQ29uc3RydWNjacOzbiBkZWwgTW9kZWxvIApSZWFsaXphcmVtb3MgdW4gbG9vcCBwYXJhIHBvZGVyIGhhY2VyIGxhIHByZWRpY2Npw7NuIGRlIHRvZG9zIGxvcyBlc3RhZG9zCmBgYHtyfQojIE9yZGVuYXIgbG9zIGRhdG9zIHBvciBhw7FvCmRhdG9zIDwtIGRhdG9zW29yZGVyKGRhdG9zJEHDsW8pLCBdCiMgRGl2aWRpciBsb3MgZGF0b3MgZW4gdW5hIGxpc3RhIHBvciBlc3RhZG8KcGFuZWxfbGlzdCA8LSBzcGxpdChkYXRvcywgZGF0b3MkRXN0YWRvKQojIEluaWNpYWxpemFyIHVuYSBsaXN0YSBwYXJhIGFsbWFjZW5hciBsb3MgbW9kZWxvcyBBUklNQQphcmltYV9tb2RlbHMgPC0gbGlzdCgpCgojIEl0ZXJhciBzb2JyZSBjYWRhIGVzdGFkbyBwYXJhIGFqdXN0YXIgdW4gbW9kZWxvIEFSSU1BIHkgaGFjZXIgcHJvbsOzc3RpY29zCmZvciAoc3RhdGUgaW4gbmFtZXMocGFuZWxfbGlzdCkpIHsKICBzdGF0ZV9kYXRhIDwtIHBhbmVsX2xpc3RbW3N0YXRlXV0KICB0c19zdGF0ZSA8LSB0cyhkYXRhID0gc3RhdGVfZGF0YSRQb2JsYWNpb24sIHN0YXJ0ID0gbWluKHN0YXRlX2RhdGEkQcOxbyksIGZyZXF1ZW5jeSA9IDEpCiAgYXJpbWFfbW9kZWwgPC0gYXV0by5hcmltYSh0c19zdGF0ZSkKICBmb3JlY2FzdF92YWx1ZXMgPC0gZm9yZWNhc3QoYXJpbWFfbW9kZWwsIGggPSA2MCkgICMgSGFjZW1vcyBwcm9uw7NzdGljb3MgaGFzdGEgZWwgYcOxbyAyMDcwCiAgY2xlYW5fc3RhdGVfbmFtZSA8LSBzdGF0ZQogIGFyaW1hX21vZGVsc1tbY2xlYW5fc3RhdGVfbmFtZV1dIDwtIGxpc3QobW9kZWwgPSBhcmltYV9tb2RlbCwgZm9yZWNhc3QgPSBmb3JlY2FzdF92YWx1ZXMpCiAgCiAgcGxvdChmb3JlY2FzdF92YWx1ZXMsIG1haW4gPSBwYXN0ZSgiUHJvbsOzc3RpY28gZGUgcG9ibGFjacOzbiBlbiIsIHN0YXRlKSkKfQoKYGBgCgojIyBFeHBvcnRhY2nDs24gZGUgZGF0b3MgcHJlZGljaG9zCgpHZW5lcmFtb3MgdW5hIGRhdGEgZnJhbWUgZGlzdGludG8gcGFyYSBwb2RlciBncmFmaWNhciBsb3MgbWFwYXMKYGBge3J9CiMgQ3JlYXIgdW4gZGF0YWZyYW1lIHZhY8OtbyBwYXJhIGFsbWFjZW5hciBsb3MgcHJvbsOzc3RpY29zIGVuIGZvcm1hIGRlIHBhbmVsCmZvcmVjYXN0X3BhbmVsIDwtIGRhdGEuZnJhbWUoc3RhdGUgPSBjaGFyYWN0ZXIoKSwgeWVhciA9IG51bWVyaWMoKSwgcHJlZGljdGlvbiA9IG51bWVyaWMoKSkKCiMgSXRlcmFyIHNvYnJlIGxvcyBub21icmVzIGVuIHBhbmVsX2xpc3QKZm9yIChzdGF0ZSBpbiBuYW1lcyhwYW5lbF9saXN0KSkgewogICMgQWNjZWRlciBhbCBwcm9uw7NzdGljbyBjb3JyZXNwb25kaWVudGUgZW4gYXJpbWFfbW9kZWxzCiAgZm9yZWNhc3QgPC0gYXJpbWFfbW9kZWxzW1tzdGF0ZV1dJGZvcmVjYXN0CiAgIyBDcmVhciB1biBkYXRhZnJhbWUgdGVtcG9yYWwgcGFyYSBlbCBlc3RhZG8gYWN0dWFsCiAgc3RhdGVfZm9yZWNhc3QgPC0gZGF0YS5mcmFtZShzdGF0ZSA9IHJlcChzdGF0ZSwgbGVuZ3RoKGZvcmVjYXN0KSksIHllYXIgPSBzZXEoMjAyMCwgMjA3OSksIHByZWRpY3Rpb24gPSBmb3JlY2FzdCkKICAjIEFncmVnYXIgZWwgZGF0YWZyYW1lIHRlbXBvcmFsIGFsIGRhdGFmcmFtZSBmb3JlY2FzdF9wYW5lbAogIGZvcmVjYXN0X3BhbmVsIDwtIHJiaW5kKGZvcmVjYXN0X3BhbmVsLCBzdGF0ZV9mb3JlY2FzdCkKfQoKIyBTZWxlY2Npb25hciB5IHJlbm9tYnJhciBjb2x1bW5hcyBwZXJ0aW5lbnRlcwpmb3JlY2FzdF9wYW5lbCA8LSBmb3JlY2FzdF9wYW5lbCAlPiUgc2VsZWN0KDE6MykKY29sbmFtZXMoZm9yZWNhc3RfcGFuZWwpIDwtIGMoIkVzdGFkbyIsICJBw7FvIiwgIlBvYmxhY2lvbiIpCgpgYGAKCiMjIENvbnN0cnVjY2nDs24gZGUgTWFwYXMKTnVldmFtZW50ZSBjcmVhcmVtb3MgdW4gbG9vcCBwYXJhIGxhIGNyZWFjacOzbiBkZSBsb3MgbWFwYXMgCmBgYHtyfQojIERlZmluaXIgYcOxb3MgZXNwZWPDrWZpY29zIHBhcmEgZWwgYW7DoWxpc2lzCmFuaW9zIDwtIGMoMjAzMCwgMjA0MCwgMjA1MCwgMjA2MCwgMjA3MCkKCiMgSXRlcmFyIHNvYnJlIGxvcyBhw7FvcyBwYXJhIGdlbmVyYXIgbWFwYXMKZm9yIChpIGluIGFuaW9zKSB7CiAgIyBGaWx0cmFyIGVsIGRhdGFmcmFtZSBkZSBwcm9uw7NzdGljb3MgcGFyYSBlbCBhw7FvIGFjdHVhbAogIGRmMiA8LSBmb3JlY2FzdF9wYW5lbCAlPiUgZmlsdGVyKEHDsW8gPT0gaSkKICAKICAjIENhbGN1bGFyIGxvcyBjdWFydGlsZXMgeSBhc2lnbmFyIGNvbG9yZXMgY29ycmVzcG9uZGllbnRlcwogIGRmMiA8LSB3aXRoaW4oZGYyLCBxdWFydGlsZSA8LSBhcy5pbnRlZ2VyKGN1dChQb2JsYWNpb24sIHF1YW50aWxlKFBvYmxhY2lvbiwgcHJvYnMgPSAwOjQvNCksIGluY2x1ZGUubG93ZXN0ID0gVFJVRSkpKQogIGRmMiRjb2xbZGYyJHF1YXJ0aWxlID09IDFdIDwtICdjYWRldGJsdWUyJwogIGRmMiRjb2xbZGYyJHF1YXJ0aWxlID09IDJdIDwtICdjYWRldGJsdWUzJwogIGRmMiRjb2xbZGYyJHF1YXJ0aWxlID09IDNdIDwtICdjYWRldGJsdWUnCiAgZGYyJGNvbFtkZjIkcXVhcnRpbGUgPT0gNF0gPC0gJ2NhZGV0Ymx1ZTQnCiAgCiAgIyBSZW5vbWJyYXIgbG9zIGVzdGFkb3MgY29uIGPDs2RpZ29zIGEgbm9tYnJlcyBjb21wbGV0b3MKICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1RYJ10gPC0gJ1RleGFzJwogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnTkQnXSA8LSAnTm9ydGggRGFrb3RhJwogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnU0QnXSA8LSAnU291dGggRGFrb3RhJwogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnS1MnXSA8LSAnS2Fuc2FzJwogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnSUEnXSA8LSAnSW93YScKICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1ZUJ10gPC0gJ1Zlcm1vbnQnCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdQQSddIDwtICdQZW5uc3lsdmFuaWEnCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdLWSddIDwtICdLZW50dWNreScKICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1dWJ10gPC0gJ1dlc3QgVmlyZ2luaWEnCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdWQSddIDwtICdWaXJnaW5pYScKICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1ROJ10gPC0gJ1Rlbm5lc3NlZScKICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ05DJ10gPC0gJ05vcnRoIENhcm9saW5hJwogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnU0MnXSA8LSAnU291dGggQ2Fyb2xpbmEnCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdHQSddIDwtICdHZW9yZ2lhJwogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnTEEnXSA8LSAnTG91aXNpYW5hJwogIAogICMgR2VuZXJhciBtYXBhIGNvbiBsb3MgY29sb3JlcyBhc2lnbmFkb3MKICBtYXAoZGF0YWJhc2UgPSAnc3RhdGUnKQogIG1hcChkYXRhYmFzZSA9ICdzdGF0ZScsIHJlZ2lvbnMgPSBkZjIkc3RhdGUsIGNvbCA9IGRmMiRjb2wsIGZpbGwgPSBUUlVFLCBhZGQgPSBUUlVFKQogIHRpdGxlKHBhc3RlKCJNYXBhIGRlbCBhw7FvIiwgaSkpCn0KCmBgYAoK