1 Importamos las librerías

library(forecast)
library(plm)
library(dplyr)
library(maps)

2 Cargar datos

datos<-read.csv("historical_state_population_by_year.csv")
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 dataframe 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))
}

LS0tDQp0aXRsZTogIkFjdGl2aWRhZCBJbnRlZ3JhZG9yYSINCmF1dGhvcjogSm9zw6kgR2FicmllbCBVc2nDsWEgTW9ncm8gQTAwODMxNDM1ICBMb3JlbmEgVmlsbGFycmVhbCBWZWdhIEEwMTcyMDgwMiBYaW1lbmENCiAgU29sw61zIElzbGFzIEEwMDgzMTM3MQ0KZGF0ZTogIjIwMjQtMDItMjEiDQpvdXRwdXQ6DQogIGh0bWxfZG9jdW1lbnQ6DQogICAgdG9jOiB0cnVlDQogICAgdG9jX2RlcHRoOiAzDQogICAgbnVtYmVyX3NlY3Rpb25zOiBUUlVFDQogICAgdG9jX2Zsb2F0OiB0cnVlDQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KLS0tDQoNCiFbXShDOlxcVXNlcnNcXHhzaV9zXFxEb3dubG9hZHNcXHVzYSBtYXAuZ2lmKQ0KDQoNCiMgSW1wb3J0YW1vcyBsYXMgbGlicmVyw61hcw0KYGBge3IgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0NCmxpYnJhcnkoZm9yZWNhc3QpDQpsaWJyYXJ5KHBsbSkNCmxpYnJhcnkoZHBseXIpDQpsaWJyYXJ5KG1hcHMpDQpgYGANCg0KIyBDYXJnYXIgZGF0b3MNCmBgYHtyfQ0KZGF0b3M8LXJlYWQuY3N2KCJoaXN0b3JpY2FsX3N0YXRlX3BvcHVsYXRpb25fYnlfeWVhci5jc3YiKQ0KY29sbmFtZXMoZGF0b3MpPC1jKCJFc3RhZG8iLCJBw7FvIiwiUG9ibGFjaW9uIikNCmBgYA0KDQoNCiMgU2VyaWUgZGUgVGllbXBvDQojIyBDb25zdHJ1Y2Npw7NuIGRlbCBNb2RlbG8gDQpSZWFsaXphcmVtb3MgdW4gbG9vcCBwYXJhIHBvZGVyIGhhY2VyIGxhIHByZWRpY2Npw7NuIGRlIHRvZG9zIGxvcyBlc3RhZG9zDQpgYGB7cn0NCiMgT3JkZW5hciBsb3MgZGF0b3MgcG9yIGHDsW8NCmRhdG9zIDwtIGRhdG9zW29yZGVyKGRhdG9zJEHDsW8pLCBdDQojIERpdmlkaXIgbG9zIGRhdG9zIGVuIHVuYSBsaXN0YSBwb3IgZXN0YWRvDQpwYW5lbF9saXN0IDwtIHNwbGl0KGRhdG9zLCBkYXRvcyRFc3RhZG8pDQojIEluaWNpYWxpemFyIHVuYSBsaXN0YSBwYXJhIGFsbWFjZW5hciBsb3MgbW9kZWxvcyBBUklNQQ0KYXJpbWFfbW9kZWxzIDwtIGxpc3QoKQ0KDQojIEl0ZXJhciBzb2JyZSBjYWRhIGVzdGFkbyBwYXJhIGFqdXN0YXIgdW4gbW9kZWxvIEFSSU1BIHkgaGFjZXIgcHJvbsOzc3RpY29zDQpmb3IgKHN0YXRlIGluIG5hbWVzKHBhbmVsX2xpc3QpKSB7DQogIHN0YXRlX2RhdGEgPC0gcGFuZWxfbGlzdFtbc3RhdGVdXQ0KICB0c19zdGF0ZSA8LSB0cyhkYXRhID0gc3RhdGVfZGF0YSRQb2JsYWNpb24sIHN0YXJ0ID0gbWluKHN0YXRlX2RhdGEkQcOxbyksIGZyZXF1ZW5jeSA9IDEpDQogIGFyaW1hX21vZGVsIDwtIGF1dG8uYXJpbWEodHNfc3RhdGUpDQogIGZvcmVjYXN0X3ZhbHVlcyA8LSBmb3JlY2FzdChhcmltYV9tb2RlbCwgaCA9IDYwKSAgIyBIYWNlbW9zIHByb27Ds3N0aWNvcyBoYXN0YSBlbCBhw7FvIDIwNzANCiAgY2xlYW5fc3RhdGVfbmFtZSA8LSBzdGF0ZQ0KICBhcmltYV9tb2RlbHNbW2NsZWFuX3N0YXRlX25hbWVdXSA8LSBsaXN0KG1vZGVsID0gYXJpbWFfbW9kZWwsIGZvcmVjYXN0ID0gZm9yZWNhc3RfdmFsdWVzKQ0KICANCiAgcGxvdChmb3JlY2FzdF92YWx1ZXMsIG1haW4gPSBwYXN0ZSgiUHJvbsOzc3RpY28gZGUgcG9ibGFjacOzbiBlbiIsIHN0YXRlKSkNCn0NCg0KYGBgDQoNCiMjIEV4cG9ydGFjacOzbiBkZSBkYXRvcyBwcmVkaWNob3MNCg0KR2VuZXJhbW9zIHVuYSBkYXRhZnJhbWUgZGlzdGludG8gcGFyYSBwb2RlciBncmFmaWNhciBsb3MgbWFwYXMNCmBgYHtyfQ0KIyBDcmVhciB1biBkYXRhZnJhbWUgdmFjw61vIHBhcmEgYWxtYWNlbmFyIGxvcyBwcm9uw7NzdGljb3MgZW4gZm9ybWEgZGUgcGFuZWwNCmZvcmVjYXN0X3BhbmVsIDwtIGRhdGEuZnJhbWUoc3RhdGUgPSBjaGFyYWN0ZXIoKSwgeWVhciA9IG51bWVyaWMoKSwgcHJlZGljdGlvbiA9IG51bWVyaWMoKSkNCg0KIyBJdGVyYXIgc29icmUgbG9zIG5vbWJyZXMgZW4gcGFuZWxfbGlzdA0KZm9yIChzdGF0ZSBpbiBuYW1lcyhwYW5lbF9saXN0KSkgew0KICAjIEFjY2VkZXIgYWwgcHJvbsOzc3RpY28gY29ycmVzcG9uZGllbnRlIGVuIGFyaW1hX21vZGVscw0KICBmb3JlY2FzdCA8LSBhcmltYV9tb2RlbHNbW3N0YXRlXV0kZm9yZWNhc3QNCiAgIyBDcmVhciB1biBkYXRhZnJhbWUgdGVtcG9yYWwgcGFyYSBlbCBlc3RhZG8gYWN0dWFsDQogIHN0YXRlX2ZvcmVjYXN0IDwtIGRhdGEuZnJhbWUoc3RhdGUgPSByZXAoc3RhdGUsIGxlbmd0aChmb3JlY2FzdCkpLCB5ZWFyID0gc2VxKDIwMjAsIDIwNzkpLCBwcmVkaWN0aW9uID0gZm9yZWNhc3QpDQogICMgQWdyZWdhciBlbCBkYXRhZnJhbWUgdGVtcG9yYWwgYWwgZGF0YWZyYW1lIGZvcmVjYXN0X3BhbmVsDQogIGZvcmVjYXN0X3BhbmVsIDwtIHJiaW5kKGZvcmVjYXN0X3BhbmVsLCBzdGF0ZV9mb3JlY2FzdCkNCn0NCg0KIyBTZWxlY2Npb25hciB5IHJlbm9tYnJhciBjb2x1bW5hcyBwZXJ0aW5lbnRlcw0KZm9yZWNhc3RfcGFuZWwgPC0gZm9yZWNhc3RfcGFuZWwgJT4lIHNlbGVjdCgxOjMpDQpjb2xuYW1lcyhmb3JlY2FzdF9wYW5lbCkgPC0gYygiRXN0YWRvIiwgIkHDsW8iLCAiUG9ibGFjaW9uIikNCg0KYGBgDQoNCiMjIENvbnN0cnVjY2nDs24gZGUgTWFwYXMNCk51ZXZhbWVudGUgY3JlYXJlbW9zIHVuIGxvb3AgcGFyYSBsYSBjcmVhY2nDs24gZGUgbG9zIG1hcGFzIA0KYGBge3J9DQojIERlZmluaXIgYcOxb3MgZXNwZWPDrWZpY29zIHBhcmEgZWwgYW7DoWxpc2lzDQphbmlvcyA8LSBjKDIwMzAsIDIwNDAsIDIwNTAsIDIwNjAsIDIwNzApDQoNCiMgSXRlcmFyIHNvYnJlIGxvcyBhw7FvcyBwYXJhIGdlbmVyYXIgbWFwYXMNCmZvciAoaSBpbiBhbmlvcykgew0KICAjIEZpbHRyYXIgZWwgZGF0YWZyYW1lIGRlIHByb27Ds3N0aWNvcyBwYXJhIGVsIGHDsW8gYWN0dWFsDQogIGRmMiA8LSBmb3JlY2FzdF9wYW5lbCAlPiUgZmlsdGVyKEHDsW8gPT0gaSkNCiAgDQogICMgQ2FsY3VsYXIgbG9zIGN1YXJ0aWxlcyB5IGFzaWduYXIgY29sb3JlcyBjb3JyZXNwb25kaWVudGVzDQogIGRmMiA8LSB3aXRoaW4oZGYyLCBxdWFydGlsZSA8LSBhcy5pbnRlZ2VyKGN1dChQb2JsYWNpb24sIHF1YW50aWxlKFBvYmxhY2lvbiwgcHJvYnMgPSAwOjQvNCksIGluY2x1ZGUubG93ZXN0ID0gVFJVRSkpKQ0KICBkZjIkY29sW2RmMiRxdWFydGlsZSA9PSAxXSA8LSAnY2FkZXRibHVlMicNCiAgZGYyJGNvbFtkZjIkcXVhcnRpbGUgPT0gMl0gPC0gJ2NhZGV0Ymx1ZTMnDQogIGRmMiRjb2xbZGYyJHF1YXJ0aWxlID09IDNdIDwtICdjYWRldGJsdWUnDQogIGRmMiRjb2xbZGYyJHF1YXJ0aWxlID09IDRdIDwtICdjYWRldGJsdWU0Jw0KICANCiAgIyBSZW5vbWJyYXIgbG9zIGVzdGFkb3MgY29uIGPDs2RpZ29zIGEgbm9tYnJlcyBjb21wbGV0b3MNCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdUWCddIDwtICdUZXhhcycNCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdORCddIDwtICdOb3J0aCBEYWtvdGEnDQogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnU0QnXSA8LSAnU291dGggRGFrb3RhJw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ0tTJ10gPC0gJ0thbnNhcycNCiAgZGYyJEVzdGFkb1tkZjIkRXN0YWRvID09ICdJQSddIDwtICdJb3dhJw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1ZUJ10gPC0gJ1Zlcm1vbnQnDQogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnUEEnXSA8LSAnUGVubnN5bHZhbmlhJw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ0tZJ10gPC0gJ0tlbnR1Y2t5Jw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1dWJ10gPC0gJ1dlc3QgVmlyZ2luaWEnDQogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnVkEnXSA8LSAnVmlyZ2luaWEnDQogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnVE4nXSA8LSAnVGVubmVzc2VlJw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ05DJ10gPC0gJ05vcnRoIENhcm9saW5hJw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ1NDJ10gPC0gJ1NvdXRoIENhcm9saW5hJw0KICBkZjIkRXN0YWRvW2RmMiRFc3RhZG8gPT0gJ0dBJ10gPC0gJ0dlb3JnaWEnDQogIGRmMiRFc3RhZG9bZGYyJEVzdGFkbyA9PSAnTEEnXSA8LSAnTG91aXNpYW5hJw0KICANCiAgIyBHZW5lcmFyIG1hcGEgY29uIGxvcyBjb2xvcmVzIGFzaWduYWRvcw0KICBtYXAoZGF0YWJhc2UgPSAnc3RhdGUnKQ0KICBtYXAoZGF0YWJhc2UgPSAnc3RhdGUnLCByZWdpb25zID0gZGYyJHN0YXRlLCBjb2wgPSBkZjIkY29sLCBmaWxsID0gVFJVRSwgYWRkID0gVFJVRSkNCiAgdGl0bGUocGFzdGUoIk1hcGEgZGVsIGHDsW8iLCBpKSkNCn0NCg0KYGBgDQoNCg==