¿Cómo se hizo y cómo se repite?

DescargarSHA-256 b2e92a4d6693a708d14b87e701866226861e5f6caf9bdb60af4230ebb052570e

aeropuerto,fecha,internacional,nacional,total,fuente
RCH,2020-01-01,652,23902,24554,datos.gov.co gb6w-ynu4
RCH,2020-02-01,529,19016,19545,datos.gov.co gb6w-ynu4
RCH,2020-03-01,283,11747,12030,datos.gov.co gb6w-ynu4
RCH,2020-04-01,0,60,60,datos.gov.co gb6w-ynu4
RCH,2020-05-01,0,21,21,datos.gov.co gb6w-ynu4
RCH,2020-06-01,0,31,31,datos.gov.co gb6w-ynu4
RCH,2020-07-01,0,14,14,datos.gov.co gb6w-ynu4
RCH,2020-08-01,0,19,19,datos.gov.co gb6w-ynu4
…

DescargarSHA-256 458ddcf511e6238fd6e53fe494f48bdd31f773fd2250556633693fa461fed569

aeropuerto,fecha,sillas,vuelos,fuente
RCH,2020-01-01,27418,188,datos.gov.co djjf-g4q4
RCH,2020-02-01,22362,152,datos.gov.co djjf-g4q4
RCH,2020-03-01,17546,116,datos.gov.co djjf-g4q4
RCH,2020-04-01,376,14,datos.gov.co djjf-g4q4
RCH,2020-05-01,376,8,datos.gov.co djjf-g4q4
RCH,2020-06-01,0,8,datos.gov.co djjf-g4q4
RCH,2020-07-01,0,6,datos.gov.co djjf-g4q4
RCH,2020-08-01,0,6,datos.gov.co djjf-g4q4
…

DescargarSHA-256 378f31c5ed9b6fc429e1e635be4db57a68042000ebe7182126cdee949bdb6d07

ap,fecha,internacional,nacional,total,sillas
BAQ,2020-01-01,26954,201021,227975,273313
BAQ,2020-02-01,29945,202045,231990,277059
BAQ,2020-03-01,16093,114143,130236,185257
BAQ,2020-04-01,48,232,280,714
BAQ,2020-05-01,341,454,795,1884
BAQ,2020-06-01,60,819,879,1054
BAQ,2020-07-01,24,560,584,1675
BAQ,2020-08-01,129,531,660,200
…

DescargarSHA-256 f32aa3c2389446cc7528d11981f673cba62e06eb590cb330eb470923a6461893

origen,anio,destinos,destinos_internacionales,pais_destino
BAQ,2020,19,5,"['ESPANA', 'ESTADOS UNIDOS', 'PANAMA', 'PERU']"
BAQ,2021,29,17,"['ANTILLAS HOLANDESAS', 'ARGENTINA', 'CHILE', 'ESPANA', 'ESTADOS UNIDOS', 'MEXICO', 'PANAMA', 'PERU']"
BAQ,2022,38,25,"['ANTILLAS HOLANDESAS', 'ARGENTINA', 'CANADA', 'CHILE', 'ESPAÑA', 'ESTADOS UNIDOS', 'MEXICO', 'PANAMA', 'PERU', 'PUERTO RICO', 'REPUBLICA DOMINICANA']"
BAQ,2023,51,34,"['ANTILLAS HOLANDESAS', 'ARGENTINA', 'BRASIL', 'CANADA', 'CHILE', 'COSTA RICA', 'ECUADOR', 'ESPAÑA', 'ESTADOS UNIDOS', 'GUATEMALA', 'INGLATERRA', 'MEXICO', 'PANAMA', 'PERU', 'PUERTO RICO', 'REPUBLICA DOMINICANA', 'VENEZUELA']"
BAQ,2024,54,36,"['ANTILLAS HOLANDESAS', 'ARGENTINA', 'BRASIL', 'CANADA', 'CHILE', 'COSTA RICA', 'ECUADOR', 'ESPAÑA', 'ESTADOS UNIDOS', 'GUATEMALA', 'INGLATERRA', 'MEXICO', 'PANAMA', 'PERU', 'PUERTO RICO', 'REPUBLICA DOMINICANA', 'SURAMERICA']"
BAQ,2025,57,40,"['ANTILLAS HOLANDESAS', 'ARGENTINA', 'BRASIL', 'CANADA', 'CHILE', 'COSTA RICA', 'ECUADOR', 'ESPAÑA', 'ESTADOS UNIDOS', 'FRANCIA', 'GUATEMALA', 'INGLATERRA', 'MEXICO', 'PANAMA', 'PERU', 'PUERTO RICO', 'REPUBLICA DOMINICANA']"
BAQ,2026,45,31,"['ANTILLAS HOLANDESAS', 'ARGENTINA', 'BRASIL', 'CANADÁ', 'CHILE', 'ECUADOR', 'ESPAÑA', 'ESTADOS UNIDOS', 'GUATEMALA', 'MÉXICO', 'PANAMÁ', 'PERÚ', 'PUERTO RICO', 'REPUBLICA DOMINICANA']"
CTG,2020,35,25,"['ARGENTINA', 'BRASIL', 'CANADA', 'CHILE', 'EL SALVADOR', 'ESPANA', 'ESTADOS UNIDOS', 'FRANCIA', 'HOLANDA', 'INGLATERRA', 'MEXICO', 'PANAMA', 'PERU']"
…

DescargarSHA-256 85db66c96225625e3e504d2273982f3dd7e5cb2077c217d133a69d7b5cced07a

aeropuerto,sentido,consulta_url,fecha_consulta,filas
SMR,embarcados,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i%2Csum%28pasajeros%29&%24where=origen%3D%27SMR%27&%24group=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i&%24limit=50000,2026-11-11T14:32:34,150
SMR,desembarcados,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i%2Csum%28pasajeros%29&%24where=destino%3D%27SMR%27&%24group=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i&%24limit=50000,2026-11-11T14:32:38,150
RCH,embarcados,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i%2Csum%28pasajeros%29&%24where=origen%3D%27RCH%27&%24group=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i&%24limit=50000,2026-11-11T14:32:38,150
RCH,desembarcados,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i%2Csum%28pasajeros%29&%24where=destino%3D%27RCH%27&%24group=a_o%2Cn_mero_de_mes%2Ctr_fico_n_i&%24limit=50000,2026-11-11T14:32:39,113
SMR,capacidad_salidas,https://www.datos.gov.co/resource/djjf-g4q4.json?%24select=a_o%2Cmes%2Csum%28sillas_ofrecidas%29%2Csum%28n_mero_de_vuelos%29&%24where=origen%3D%27SMR%27&%24group=a_o%2Cmes&%24limit=50000,2026-11-11T14:32:40,77
SMR,capacidad_llegadas,https://www.datos.gov.co/resource/djjf-g4q4.json?%24select=a_o%2Cmes%2Csum%28sillas_ofrecidas%29%2Csum%28n_mero_de_vuelos%29&%24where=destino%3D%27SMR%27&%24group=a_o%2Cmes&%24limit=50000,2026-11-11T14:32:49,77
RCH,capacidad_salidas,https://www.datos.gov.co/resource/djjf-g4q4.json?%24select=a_o%2Cmes%2Csum%28sillas_ofrecidas%29%2Csum%28n_mero_de_vuelos%29&%24where=origen%3D%27RCH%27&%24group=a_o%2Cmes&%24limit=50000,2026-11-11T14:32:50,78
RCH,capacidad_llegadas,https://www.datos.gov.co/resource/djjf-g4q4.json?%24select=a_o%2Cmes%2Csum%28sillas_ofrecidas%29%2Csum%28n_mero_de_vuelos%29&%24where=destino%3D%27RCH%27&%24group=a_o%2Cmes&%24limit=50000,2026-11-11T14:32:53,78
…

DescargarSHA-256 18e1b574de45905b5268aee87815bfd2eff45d94e16001d875d84cce175730c0

servicio,consulta,fecha
https://www.datos.gov.co/resource/gb6w-ynu4.json,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=origen+as+ap%2Ca_o%2Cn_mero_de_mes%2Ctr_fico_n_i%2Csum%28pasajeros%29&%24where=origen+in+%28%27SMR%27%2C%27RCH%27%2C%27CZU%27%2C%27TLU%27%2C%27CVN%27%2C%27CTG%27%2C%27BAQ%27%29&%24group=origen%2Ca_o%2Cn_mero_de_mes%2Ctr_fico_n_i&%24limit=50000,2026-11-12T07:43:25
https://www.datos.gov.co/resource/gb6w-ynu4.json,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=destino+as+ap%2Ca_o%2Cn_mero_de_mes%2Ctr_fico_n_i%2Csum%28pasajeros%29&%24where=destino+in+%28%27SMR%27%2C%27RCH%27%2C%27CZU%27%2C%27TLU%27%2C%27CVN%27%2C%27CTG%27%2C%27BAQ%27%29&%24group=destino%2Ca_o%2Cn_mero_de_mes%2Ctr_fico_n_i&%24limit=50000,2026-11-12T07:43:28
https://www.datos.gov.co/resource/djjf-g4q4.json,https://www.datos.gov.co/resource/djjf-g4q4.json?%24select=origen+as+ap%2Ca_o%2Cmes%2Csum%28sillas_ofrecidas%29&%24where=origen+in+%28%27SMR%27%2C%27RCH%27%2C%27CZU%27%2C%27TLU%27%2C%27CVN%27%2C%27CTG%27%2C%27BAQ%27%29&%24group=origen%2Ca_o%2Cmes&%24limit=50000,2026-11-12T07:43:33
https://www.datos.gov.co/resource/djjf-g4q4.json,https://www.datos.gov.co/resource/djjf-g4q4.json?%24select=destino+as+ap%2Ca_o%2Cmes%2Csum%28sillas_ofrecidas%29&%24where=destino+in+%28%27SMR%27%2C%27RCH%27%2C%27CZU%27%2C%27TLU%27%2C%27CVN%27%2C%27CTG%27%2C%27BAQ%27%29&%24group=destino%2Ca_o%2Cmes&%24limit=50000,2026-11-12T07:43:34
https://www.datos.gov.co/resource/gb6w-ynu4.json,https://www.datos.gov.co/resource/gb6w-ynu4.json?%24select=origen%2Ca_o%2Cdestino%2Cpais_destino%2Csum%28pasajeros%29&%24where=origen+in+%28%27SMR%27%2C%27RCH%27%2C%27CZU%27%2C%27TLU%27%2C%27CVN%27%2C%27CTG%27%2C%27BAQ%27%29+and+pasajeros+%3E+0&%24group=origen%2Ca_o%2Cdestino%2Cpais_destino&%24limit=50000,2026-11-12T07:43:36

…

DescargarSHA-256 6030ecfd0a0ba5c8808455662966dea63748fd1444d1a74ca2d2651beea78229

aeropuerto,origen,objetivo,h,modelo,pronostico,observado,error,ape
RCH,2024-12-01,2025-01-01,1,Ingenuo estacional,53738.0,57935.0,4197.0,7.244325537240011
RCH,2024-12-01,2025-01-01,1,ETS,61276.08966299964,57935.0,-3341.0896629996423,5.766962394061694
RCH,2024-12-01,2025-01-01,1,SARIMA,66289.8706251288,57935.0,-8354.870625128795,14.421110943520835
RCH,2024-12-01,2025-01-01,1,XGBoost,60678.4566552639,57935.0,-2743.4566552639008,4.735404600438251
RCH,2024-12-01,2025-01-01,1,Capacidad × ocupación (programada),56391.03493304248,57935.0,1543.9650669575203,2.664995368874636
RCH,2024-12-01,2025-01-01,1,SARIMAX + capacidad (programada),56027.495602442614,57935.0,1907.5043975573863,3.2924905455379063
RCH,2024-12-01,2025-01-01,1,XGBoost + capacidad (programada),49999.847453832626,57935.0,7935.152546167374,13.696647184201904
RCH,2024-12-01,2025-01-01,1,Capacidad × ocupación (estimada),66198.94537585092,57935.0,-8263.945375850919,14.264167387332217
…

DescargarSHA-256 5a8c20d01a55a94a6afc59ebe7ad54918ad319dc4691c3425c555a54fe01aa14

aeropuerto,bloque,modelo,MAPE_%,MASE,n
RCH,h1-3,SARIMAX + capacidad (programada),4.559147234955397,0.3086059400540041,51.0
RCH,h1-3,Ingenuo estacional,4.926794803624814,0.34691220339885825,51.0
RCH,h1-3,Capacidad × ocupación (programada),5.183184038576818,0.3487991485483154,51.0
RCH,h1-3,ETS,5.703852492804424,0.388309002939068,51.0
RCH,h1-3,SARIMAX + capacidad (estimada),6.101745092223645,0.4177104635356704,51.0
RCH,h1-3,XGBoost + capacidad (programada),6.1394847328184206,0.43269722520490594,51.0
RCH,h1-3,SARIMA,6.804970976013591,0.4636568071173935,51.0
RCH,h1-3,Capacidad × ocupación (estimada),7.973893335746071,0.5458926048249619,51.0
…

DescargarSHA-256 35b91c875f8d37c91fa3b46be03cd375864e9a92ef17105afeaaa7581fa04e31

aeropuerto,modelo,h,MAPE_%,MASE,n
RCH,Capacidad × ocupación (estimada),1,7.04259706289692,0.4965031951743589,18.0
RCH,Capacidad × ocupación (estimada),2,8.446002637061994,0.5699693119001379,17.0
RCH,Capacidad × ocupación (estimada),3,8.519985510053202,0.5758741894145163,16.0
RCH,Capacidad × ocupación (estimada),4,8.592690354529193,0.5915270614193673,15.0
RCH,Capacidad × ocupación (estimada),5,9.868681812955495,0.6873155026399538,14.0
RCH,Capacidad × ocupación (estimada),6,9.381177521052003,0.6752830419973817,13.0
RCH,Capacidad × ocupación (programada),1,5.237866031000732,0.35306085909728935,18.0
RCH,Capacidad × ocupación (programada),2,5.389211364066973,0.3603107755995968,17.0
…

DescargarSHA-256 33399afdfccd9a33f54daf9aa5031bd3f375a2f940f970a87a9400584b4fcff0

aeropuerto,h,modelo,n,DM_HLN,p_valor,lectura
RCH,1,ETS,18,0.06030118449250476,0.9526189918914969,peor que la referencia (no significativo)
RCH,1,SARIMA,18,1.5212965821711464,0.14656768869675457,peor que la referencia (no significativo)
RCH,1,XGBoost,18,2.168614539327769,0.04458804606525601,peor que la referencia (significativo al 5 %)
RCH,1,Capacidad × ocupación (programada),18,0.36771703912858794,0.7176223878004655,peor que la referencia (no significativo)
RCH,1,SARIMAX + capacidad (programada),18,-0.05987311381472571,0.9529549170188881,mejor que la referencia (no significativo)
RCH,1,XGBoost + capacidad (programada),18,1.6617227229027787,0.11489627862586566,peor que la referencia (no significativo)
RCH,1,Capacidad × ocupación (estimada),18,1.9305961117670012,0.07038642554401445,peor que la referencia (no significativo)
RCH,1,SARIMAX + capacidad (estimada),18,0.7248778238767013,0.47838649221831314,peor que la referencia (no significativo)
…

DescargarSHA-256 703b37664a3c75c9de4e2e17f39e370e1a26a7c45adebca8468f009b07cca228

,XGBoost,XGBoost + capacidad (programada)
es_RCH,0.004141056,0.005398220382630825
g_capacidad,,0.09674830314132475
g_media3,0.028982593,0.013329492881894112
g_objetivo_hace12,0.09641151,0.03002931736409664
g_origen,0.035365473,0.008871021680533886
h,0.006903906,0.0008095098892226815
mes_objetivo,0.024607923,0.012652028352022171
nivel_rel,0.036783423,0.012929890304803848
…

DescargarSHA-256 a4e401bb5adee4ecf898d1e5afa643a9fc1db808c8c3dcdedc0904243bd112dc

aeropuerto,fecha,internacional,nacional,total,fuente_x,sillas,vuelos,fuente_y,ocupacion,media_12m_previa,exceso,alta_documentada
RCH,2022-01-01,361,33632,33993,datos.gov.co gb6w-ynu4,37248,212,datos.gov.co djjf-g4q4,0.9126127577319587,0.8323227526722308,0.08029000505972794,True
RCH,2022-02-01,260,29328,29588,datos.gov.co gb6w-ynu4,33420,206,datos.gov.co djjf-g4q4,0.8853381208856972,0.8375052887842359,0.047832832101461276,False
RCH,2022-03-01,274,32481,32755,datos.gov.co gb6w-ynu4,39160,221,datos.gov.co djjf-g4q4,0.8364402451481103,0.8449267995917734,-0.008486554443663086,False
RCH,2022-04-01,515,37250,37765,datos.gov.co gb6w-ynu4,47058,266,datos.gov.co djjf-g4q4,0.8025202941051468,0.8440441758720282,-0.041523881766881376,True
RCH,2022-05-01,423,35424,35847,datos.gov.co gb6w-ynu4,46408,265,datos.gov.co djjf-g4q4,0.7724314773314945,0.8578119397880889,-0.08538046245659436,False
RCH,2022-06-01,306,36248,36554,datos.gov.co gb6w-ynu4,44122,257,datos.gov.co djjf-g4q4,0.828475590408413,0.8678648865236421,-0.03938929611522912,True
RCH,2022-07-01,363,43396,43759,datos.gov.co gb6w-ynu4,49682,284,datos.gov.co djjf-g4q4,0.8807817720703676,0.86742311642034,0.013358655650027584,True
RCH,2022-08-01,421,39302,39723,datos.gov.co gb6w-ynu4,45076,256,datos.gov.co djjf-g4q4,0.8812450084302068,0.867311749938822,0.013933258491384826,False
…

DescargarSHA-256 9535387115e7ef6aac86dc9c7144dd7d8383eecd7a8c6ca5a6d72549a2c5d4b3

region,fecha,ocupacion,media_12m_previa,exceso,alta_documentada
HCA,2022-01-01,0.615586282834778,0.46103711949616694,0.15454916333861107,True
HCA,2022-02-01,0.518966007970916,0.4770366537171901,0.04192935425372585,False
HCA,2022-03-01,0.558842844810114,0.4905253129841094,0.06831753182600464,False
HCA,2022-04-01,0.576627116524489,0.5047627141897978,0.07186440233469116,True
HCA,2022-05-01,0.49616777240606197,0.5315210091829693,-0.035353236776907326,False
HCA,2022-06-01,0.569883696146666,0.5493429273886701,0.02054076875799593,True
HCA,2022-07-01,0.605571553365246,0.5602211721405967,0.045350381224649294,True
HCA,2022-08-01,0.5791353649930879,0.5674773392390341,0.011658025754053836,False
…

DescargarSHA-256 04fdd44925e2f26ebc00a9f6750421336f7f549dfd3d77a4c6d6cef1c32aeace

anio,mes,visitantes_mes,dias_abiertos,dias_temporada_alta,fuente_total,media_diaria,alta_documentada
2022,1,55814.400313348626,31,15,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1800.4645262370525,True
2022,2,20849.201219632116,13,0,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1603.7847092024704,False
2022,3,51199.045784878035,31,0,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1651.5821220928399,False
2022,4,51962.100387644474,30,8,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1732.0700129214824,True
2022,5,49777.846162982336,31,0,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1605.7369729994302,False
2022,6,26868.695351401035,15,15,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1791.246356760069,True
2022,7,60054.81566214294,31,15,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1937.252118133643,True
2022,8,55670.29557525241,31,0,"PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1",1795.8159862984646,False
…

DescargarSHA-256 7ef1b2c5a778a7fa2979e3d28ce5f62dca03ed805113652cdde119b06cd28586

"""P1 · Serie mensual de pasajeros de Santa Marta (SMR) y Riohacha (RCH).

Fuente: Aerocivil, «Transporte Aéreo Comercial – Tráfico Origen-Destino (Colombia)»,
datos.gov.co, conjunto gb6w-ynu4 (2020 en adelante), API SODA abierta.
Pasajeros del aeropuerto = embarcados (origen = aeropuerto) + desembarcados (destino = aeropuerto).
Guarda datos/serie_pasajeros_datosgov.csv y registra la consulta en datos/_procedencia.csv.
"""
import csv, datetime as dt, sys, time
from pathlib import Path
import pandas as pd, requests

sys.stdout.reconfigure(encoding="utf-8")
DATOS = Path(__file__).resolve().parents[1] / "datos"
URL = "https://www.datos.gov.co/resource/gb6w-ynu4.json"
AEROPUERTOS = {"SMR": "Santa Marta (Simón Bolívar)", "RCH": "Riohacha (Almirante Padilla)"}

def consulta(campo, codigo):
    params = {"$select": "a_o,n_mero_de_mes,tr_fico_n_i,sum(pasajeros)",
              "$where": f"{campo}='{codigo}'",
              "$group": "a_o,n_mero_de_mes,tr_fico_n_i", "$limit": 50000}
    for intento in range(5):
        r = requests.get(URL, params=params, timeout=180)
        if r.ok:
            break
        time.sleep(5 * (intento + 1))
    r.raise_for_status()
    return pd.DataFrame(r.json()).rename(columns={"sum_pasajeros": "pax"}), r.url

filas, procedencia = [], []
for cod in AEROPUERTOS:
    for campo, sentido in (("origen", "embarcados"), ("destino", "desembarcados")):
        df, url = consulta(campo, cod)
        df["aeropuerto"], df["sentido"] = cod, sentido
        filas.append(df)
        procedencia.append([cod, sentido, url, dt.datetime.now().isoformat(timespec="seconds"), len(df)])

d = pd.concat(filas)
d["pax"] = pd.to_numeric(d["pax"]); d["anio"] = d["a_o"].astype(int); d["mes"] = d["n_mero_de_mes"].astype(int)
d["fecha"] = pd.to_datetime(dict(year=d.anio, month=d.mes, day=1))
d["trafico"] = d["tr_fico_n_i"].map({"N": "nacional", "I": "internacional"}).fillna(d["tr_fico_n_i"])
serie = (d.pivot_table(index=["aeropuerto", "fecha"], columns="trafico", values="pax", aggfunc="sum", fill_value=0)
          .reset_index())
serie["total"] = serie.drop(columns=["aeropuerto", "fecha"]).sum(axis=1)
serie["fuente"] = "datos.gov.co gb6w-ynu4"
serie.to_csv(DATOS / "serie_pasajeros_datosgov.csv", index=False)

# Capacidad: sillas ofrecidas y vuelos (Aerocivil, «Transporte Aéreo – Tráfico por Equipo», djjf-g4q4)
URL_CAP = "https://www.datos.gov.co/resource/djjf-g4q4.json"
cap_filas = []
for cod in AEROPUERTOS:
    for campo, sentido in (("origen", "salidas"), ("destino", "llegadas")):
        params = {"$select": "a_o,mes,sum(sillas_ofrecidas),sum(n_mero_de_vuelos)", "$where": f"{campo}='{cod}'",
                  "$group": "a_o,mes", "$limit": 50000}
        for intento in range(5):
            r = requests.get(URL_CAP, params=params, timeout=180)
            if r.ok:
                break
            time.sleep(5 * (intento + 1))
        r.raise_for_status()
        c = pd.DataFrame(r.json()); c["aeropuerto"], c["sentido"] = cod, sentido
        cap_filas.append(c)
        procedencia.append([cod, f"capacidad_{sentido}", r.url, dt.datetime.now().isoformat(timespec="seconds"), len(c)])
c = pd.concat(cap_filas)
c["fecha"] = pd.to_datetime(dict(year=c.a_o.astype(int), month=c.mes.astype(int), day=1))
c["sillas"] = pd.to_numeric(c["sum_sillas_ofrecidas"]); c["vuelos"] = pd.to_numeric(c["sum_n_mero_de_vuelos"])
cap = c.groupby(["aeropuerto", "fecha"])[["sillas", "vuelos"]].sum().reset_index()
cap["fuente"] = "datos.gov.co djjf-g4q4"
cap.to_csv(DATOS / "serie_capacidad_datosgov.csv", index=False)
print(cap.groupby("aeropuerto").agg(desde=("fecha", "min"), hasta=("fecha", "max"), meses=("fecha", "count")))

with open(DATOS / "_procedencia.csv", "w", newline="", encoding="utf-8") as f:
    w = csv.writer(f); w.writerow(["aeropuerto", "sentido", "consulta_url", "fecha_consulta", "filas"]); w.writerows(procedencia)

print(serie.groupby("aeropuerto").agg(desde=("fecha", "min"), hasta=("fecha", "max"), meses=("fecha", "count"),
                                      pax_2025=("total", lambda s: s[serie.loc[s.index, "fecha"].dt.year == 2025].sum())))
print(serie.tail(8).to_string())

DescargarSHA-256 bea9ecfe19939c2090f2c3b8b4085d085c56c4703fcfddf5678860d30ddc4e6e

# %% [markdown]
# # P1 · Prueba mínima reproducible de pronóstico de demanda (E-13)
# ATARRAYA AZUL · Nota de concepto TRL de entrada · CIDS, Universidad del Magdalena
#
# **Pregunta.** ¿Un modelo estadístico o de aprendizaje automático pronostica los pasajeros
# mensuales de los aeropuertos de Santa Marta (SMR) y Riohacha (RCH) con menor error fuera de
# muestra que la referencia estacional ingenua, a horizontes de 3 y 6 meses?
#
# **Datos.** Aerocivil, «Transporte Aéreo Comercial – Tráfico Origen-Destino (Colombia)»,
# datos.gov.co (gb6w-ynu4). Pasajeros = embarcados + desembarcados. Se descargan con
# `01_descargar_datosgov.py`.
#
# **Diseño.** Entrenamiento desde ene-2022 (se excluyen 2020–2021 por la pandemia). Evaluación con
# origen móvil: en cada mes de origen entre dic-2024 y may-2026 se reentrena con los datos
# disponibles hasta ese mes y se pronostican h = 1…6 meses. Ningún modelo ve datos posteriores
# a su origen.
#
# **Modelos.** (0) Ingenuo estacional: ŷ(t) = y(t−12), referencia. (1) ETS aditivo sobre log con
# tendencia amortiguada. (2) SARIMA(0,1,1)(0,1,1)₁₂ sobre log (modelo «aerolínea»). (3) XGBoost
# global (ambos aeropuertos), estrategia directa sobre el crecimiento interanual en log, con
# explicación SHAP.
#
# **Criterio (nota de concepto, 7.3).** Error fuera de muestra (MAPE o MASE) inferior al de la
# referencia. Se complementa con la prueba de Diebold-Mariano (corrección HLN).

# %% Configuración
import hashlib, json, platform, sys, warnings
from pathlib import Path
import numpy as np, pandas as pd
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
import statsmodels, xgboost as xgb
from scipy import stats
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.statespace.sarimax import SARIMAX

sys.stdout.reconfigure(encoding="utf-8")
warnings.filterwarnings("ignore")
RAIZ = Path(__file__).resolve().parents[1]
DATOS, RES = RAIZ / "datos", RAIZ / "resultados"
RES.mkdir(exist_ok=True)
SEMILLA = 2026
INICIO_TRAIN = pd.Timestamp("2022-01-01")
PRIMER_ORIGEN = pd.Timestamp("2024-12-01")
H_MAX = 6
NOMBRES = {"SMR": "Santa Marta (SMR)", "RCH": "Riohacha (RCH)"}
BASE = ["Ingenuo estacional", "ETS", "SARIMA", "XGBoost"]
CAPM = ["Capacidad × ocupación", "SARIMAX + capacidad", "XGBoost + capacidad"]
MODOS = ["programada", "estimada"]
MODELOS = BASE + [f"{m} ({md})" for md in MODOS for m in CAPM]

# %% Datos
archivo = DATOS / "serie_pasajeros_datosgov.csv"
sha = hashlib.sha256(archivo.read_bytes()).hexdigest()
s = pd.read_csv(archivo, parse_dates=["fecha"])
Y = (s.pivot(index="fecha", columns="aeropuerto", values="total")
       .asfreq("MS").fillna(0.0).loc[INICIO_TRAIN:])
ULTIMO = Y.index.max()
print(f"Serie {Y.index.min():%Y-%m} a {ULTIMO:%Y-%m}; {len(Y)} meses; sha256={sha[:16]}…")

# Capacidad: sillas ofrecidas (salidas + llegadas). Se usa desde 2020 solo como base de comparación interanual.
archivo_cap = DATOS / "serie_capacidad_datosgov.csv"
sha_cap = hashlib.sha256(archivo_cap.read_bytes()).hexdigest()
CAPT = (pd.read_csv(archivo_cap, parse_dates=["fecha"]).pivot(index="fecha", columns="aeropuerto", values="sillas")
          .asfreq("MS").fillna(0.0))
MAS = lambda f, k: f + pd.DateOffset(months=k)

def cap_hat(a, o, t, modo):
    """Sillas del mes t vistas desde el origen o.
    programada: las sillas efectivamente ofrecidas (aproxima la programación que las aerolíneas publican
                con meses de anticipación; cota superior).
    estimada:   sin información futura: sillas del mismo mes del año anterior × crecimiento interanual
                medio de la capacidad en los tres meses previos al origen."""
    if modo == "programada" and t in CAPT.index:
        return CAPT.at[t, a]
    # «estimada», o meses sin dato publicado (fuera de la ventana evaluada)
    g = np.mean([np.log(CAPT.at[MAS(o, -k), a] / CAPT.at[MAS(o, -k - 12), a]) for k in range(3)])
    return CAPT.at[MAS(t, -12), a] * np.exp(g)

# %% Modelos estadísticos por aeropuerto
def f_ingenuo(y, h):
    return np.array([y.iloc[len(y) - 12 + (k - 1) % 12] if k <= 12 else np.nan for k in range(1, h + 1)])

def f_ets(y, h):
    m = ExponentialSmoothing(np.log(y.values), trend="add", damped_trend=True, seasonal="add",
                             seasonal_periods=12, initialization_method="estimated").fit()
    return np.exp(m.forecast(h))

def f_sarima(y, h):
    m = SARIMAX(np.log(y.values), order=(0, 1, 1), seasonal_order=(0, 1, 1, 12),
                enforce_stationarity=False, enforce_invertibility=False).fit(disp=False)
    return np.exp(m.forecast(h))

# %% Modelos con capacidad
def f_cxo(a, y, o, modo, h):
    """Capacidad × ocupación: sillas del mes × ocupación (pasajeros/sillas) del mismo mes del año anterior."""
    out = []
    for k in range(1, h + 1):
        t = MAS(o, k); b = MAS(t, -12)
        out.append(cap_hat(a, o, t, modo) * y.at[b] / CAPT.at[b, a])
    return np.array(out)

def f_sarimax_cap(a, y, o, modo, h):
    """Regresión del log de pasajeros sobre el log de sillas, con errores AR(1) y AR estacional(1)."""
    ex = np.log(CAPT.loc[y.index, a].values)
    m = SARIMAX(np.log(y.values), exog=ex, order=(1, 0, 0), seasonal_order=(1, 0, 0, 12), trend="c",
                enforce_stationarity=False, enforce_invertibility=False).fit(disp=False)
    ex_f = np.log([cap_hat(a, o, MAS(o, k), modo) for k in range(1, h + 1)])
    return np.exp(m.forecast(h, exog=ex_f.reshape(-1, 1)))

# %% XGBoost global, directo, sobre crecimiento interanual en log
FEATS = ["h", "mes_objetivo", "es_RCH", "g_origen", "g_media3", "g_objetivo_hace12", "nivel_rel"]
FEATS_CAP = FEATS + ["g_capacidad"]

def filas_supervisadas(Yo, origenes, con_objetivo=True, modo=None):
    L = np.log(Yo.where(Yo > 0))
    G = L - L.shift(12)                        # crecimiento interanual en log
    out = []
    for a in Yo.columns:
        for o in origenes:
            if o not in G.index or pd.isna(G.at[o, a]):
                continue
            g3 = G[a].loc[:o].tail(3)
            if g3.isna().any():
                continue
            for h in range(1, H_MAX + 1):
                t = o + pd.DateOffset(months=h)
                base = t - pd.DateOffset(months=12)
                if base not in L.index:
                    continue
                fila = {"aeropuerto": a, "origen": o, "objetivo": t, "h": h, "mes_objetivo": t.month,
                        "es_RCH": int(a == "RCH"), "g_origen": G.at[o, a], "g_media3": g3.mean(),
                        "g_objetivo_hace12": G.at[base, a] if not pd.isna(G.at[base, a]) else 0.0,
                        "nivel_rel": L.at[o, a] - L[a].loc[:o].tail(12).mean(), "base": Yo.at[base, a]}
                if modo:
                    fila["g_capacidad"] = np.log(cap_hat(a, o, t, modo) / CAPT.at[base, a])
                if con_objetivo:
                    if t not in G.index or pd.isna(G.at[t, a]):
                        continue
                    fila["y_obj"] = G.at[t, a]
                out.append(fila)
    return pd.DataFrame(out)

def ajustar_xgb(Yo, modo=None):
    feats = FEATS_CAP if modo else FEATS
    entren = filas_supervisadas(Yo, Yo.index, modo=modo)
    entren = entren[entren["objetivo"] <= Yo.index.max()]
    modelo = xgb.XGBRegressor(n_estimators=300, max_depth=2, learning_rate=0.05, subsample=0.9,
                              colsample_bytree=0.9, min_child_weight=3, reg_lambda=1.0,
                              random_state=SEMILLA, n_jobs=1)
    modelo.fit(entren[feats], entren["y_obj"])
    return modelo, len(entren)

def pronosticar_xgb(Yo, o, modo=None):
    """Devuelve {aeropuerto: array h=1..H_MAX} y las filas con contribuciones SHAP."""
    feats = FEATS_CAP if modo else FEATS
    modelo, _ = ajustar_xgb(Yo, modo)
    X = filas_supervisadas(Yo, [o], con_objetivo=False, modo=modo)
    X["pred"] = X["base"].values * np.exp(modelo.predict(X[feats]))
    contrib = modelo.get_booster().predict(xgb.DMatrix(X[feats]), pred_contribs=True)
    sal = {a: X.loc[X.aeropuerto == a].sort_values("h")["pred"].values for a in Yo.columns}
    sh = [dict(zip(feats + ["sesgo"], contrib[i]), aeropuerto=f.aeropuerto, horizonte=f.h, objetivo=f.objetivo)
          for i, f in enumerate(X.itertuples())]
    return sal, sh

def todos_los_pronosticos(Yo, o, modos=MODOS):
    xg_base, sh_base = pronosticar_xgb(Yo, o)
    xg_cap = {md: pronosticar_xgb(Yo, o, md) for md in modos}
    prev = {}
    for a in Yo.columns:
        y = Yo[a]
        p = {"Ingenuo estacional": f_ingenuo(y, H_MAX), "ETS": f_ets(y, H_MAX), "SARIMA": f_sarima(y, H_MAX),
             "XGBoost": xg_base[a]}
        for md in modos:
            p[f"Capacidad × ocupación ({md})"] = f_cxo(a, y, o, md, H_MAX)
            p[f"SARIMAX + capacidad ({md})"] = f_sarimax_cap(a, y, o, md, H_MAX)
            p[f"XGBoost + capacidad ({md})"] = xg_cap[md][0][a]
        prev[a] = p
    shap = {"XGBoost": sh_base}
    if "programada" in xg_cap:
        shap["XGBoost + capacidad (programada)"] = xg_cap["programada"][1]
    return prev, shap

# %% Evaluación con origen móvil
registros, shap_filas = [], []
origenes = pd.date_range(PRIMER_ORIGEN, ULTIMO - pd.DateOffset(months=1), freq="MS")
for o in origenes:
    Yo = Y.loc[:o]
    prev, shap_o = todos_los_pronosticos(Yo, o)
    for a in Y.columns:
        for h in range(1, H_MAX + 1):
            t = MAS(o, h)
            if t > ULTIMO:
                continue
            for m in MODELOS:
                registros.append({"aeropuerto": a, "origen": o, "objetivo": t, "h": h, "modelo": m,
                                  "pronostico": float(prev[a][m][h - 1]), "observado": float(Y.at[t, a])})
    for nombre_m, filas in shap_o.items():
        shap_filas += [dict(f, modelo=nombre_m) for f in filas if f["objetivo"] <= ULTIMO]

P = pd.DataFrame(registros)
P["error"] = P["observado"] - P["pronostico"]
P["ape"] = (P["error"].abs() / P["observado"]) * 100
P.to_csv(RES / "pronosticos_origen_movil.csv", index=False)

# %% Métricas: MAPE y MASE (escala = MAE del ingenuo estacional dentro de muestra, 2023-01 a 2024-12)
escala = {a: (Y[a] - Y[a].shift(12)).loc["2023-01":"2024-12"].abs().mean() for a in Y.columns}
P["ase"] = P["error"].abs() / P["aeropuerto"].map(escala)

def resumen(df):
    return pd.Series({"MAPE_%": df["ape"].mean(), "MASE": df["ase"].mean(), "n": len(df)})

tab_h = P.groupby(["aeropuerto", "modelo", "h"]).apply(resumen).reset_index()
P["bloque"] = np.where(P["h"] <= 3, "h1-3", "h4-6")
tab = (P.groupby(["aeropuerto", "bloque", "modelo"]).apply(resumen).reset_index())
tab_total = P.groupby(["aeropuerto", "modelo"]).apply(resumen).reset_index().assign(bloque="h1-6")
tabla = pd.concat([tab, tab_total]).sort_values(["aeropuerto", "bloque", "MASE"])
tabla.to_csv(RES / "tabla_errores.csv", index=False)
tab_h.to_csv(RES / "tabla_errores_por_horizonte.csv", index=False)

# %% Prueba de Diebold-Mariano (pérdida cuadrática, corrección Harvey-Leybourne-Newbold)
def dm_hln(e1, e2, h):
    d = e1 ** 2 - e2 ** 2
    n = len(d); dbar = d.mean()
    gam = [np.sum((d[k:] - dbar) * (d[:n - k] - dbar)) / n for k in range(h)]
    var = (gam[0] + 2 * sum(gam[1:])) / n
    if var <= 0:
        return np.nan, np.nan
    dm = dbar / np.sqrt(var)
    hln = dm * np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)
    return hln, 2 * stats.t.sf(abs(hln), df=n - 1)

dm_filas = []
for a in Y.columns:
    for h in (1, 3, 6):
        base = P[(P.aeropuerto == a) & (P.h == h) & (P.modelo == "Ingenuo estacional")].set_index("objetivo")["error"]
        for m in MODELOS[1:]:
            e = P[(P.aeropuerto == a) & (P.h == h) & (P.modelo == m)].set_index("objetivo")["error"].reindex(base.index)
            est, p = dm_hln(e.values, base.values, h)
            if np.isnan(est):
                lectura = "no calculable (varianza no positiva con n pequeño)"
            else:
                lectura = (("mejor que la referencia" if est < 0 else "peor que la referencia")
                           + (" (significativo al 5 %)" if p < 0.05 else " (no significativo)"))
            dm_filas.append({"aeropuerto": a, "h": h, "modelo": m, "n": len(base), "DM_HLN": est,
                             "p_valor": p, "lectura": lectura})
dmt = pd.DataFrame(dm_filas); dmt.to_csv(RES / "diebold_mariano.csv", index=False)

# %% SHAP (contribuciones TreeSHAP nativas de XGBoost)
SH = pd.DataFrame(shap_filas)
ETIQ = {"g_origen": "Crecimiento interanual en el origen", "g_media3": "Crecimiento interanual (media 3 meses)",
        "g_objetivo_hace12": "Crecimiento del mismo mes el año anterior", "nivel_rel": "Nivel reciente vs. media 12 m",
        "mes_objetivo": "Mes del año pronosticado", "h": "Horizonte (meses)", "es_RCH": "Aeropuerto (RCH vs. SMR)",
        "g_capacidad": "Crecimiento interanual de las sillas ofrecidas"}
imp_por = {}
for nombre_m, feats in (("XGBoost", FEATS), ("XGBoost + capacidad (programada)", FEATS_CAP)):
    imp_por[nombre_m] = SH[SH.modelo == nombre_m][feats].abs().mean().sort_values(ascending=True)
imp = imp_por["XGBoost + capacidad (programada)"]
pd.DataFrame(imp_por).to_csv(RES / "shap_importancia.csv")
fig, ax = plt.subplots(figsize=(7.5, 3.8))
ax.barh([ETIQ[i] for i in imp.index], imp.values, color="#2a7f8f")
ax.set_xlabel("Media |SHAP| sobre el crecimiento interanual en log")
ax.set_title("P1 · Qué explica el pronóstico de XGBoost + capacidad (SHAP)")
fig.tight_layout(); fig.savefig(RES / "fig_shap.png", dpi=160); plt.close(fig)

# %% Gráficos de pronóstico (h = 3) por aeropuerto
COL = {"Ingenuo estacional": "#9a9a9a", "ETS": "#d9822b",
       "Capacidad × ocupación (programada)": "#2a7f8f", "Capacidad × ocupación (estimada)": "#7a4fa3"}
fig, axes = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
for ax, a in zip(axes, ["SMR", "RCH"]):
    ax.plot(Y.index, Y[a] / 1000, color="black", lw=1.6, label="Observado")
    for m in COL:
        q = P[(P.aeropuerto == a) & (P.h == 3) & (P.modelo == m)].sort_values("objetivo")
        ax.plot(q["objetivo"], q["pronostico"] / 1000, color=COL[m], lw=1.3, marker="o", ms=3,
                ls="--" if m == "Ingenuo estacional" else "-", label=f"{m} (h=3)")
    ax.axvspan(PRIMER_ORIGEN, ULTIMO, color="#2a7f8f", alpha=0.06)
    ax.set_ylabel("Miles de pasajeros"); ax.set_title(NOMBRES[a]); ax.grid(alpha=0.25)
axes[0].legend(ncol=2, fontsize=8, loc="upper left")
fig.suptitle("P1 · Pasajeros mensuales y pronósticos a 3 meses (origen móvil)")
fig.tight_layout(); fig.savefig(RES / "fig_pronostico_h3.png", dpi=160); plt.close(fig)

# %% Registro de reproducibilidad
meta = {"datos": str(archivo.name), "sha256_datos": sha, "datos_capacidad": archivo_cap.name, "sha256_capacidad": sha_cap,
        "serie": f"{Y.index.min():%Y-%m} a {ULTIMO:%Y-%m}",
        "origenes": f"{origenes.min():%Y-%m} a {origenes.max():%Y-%m} ({len(origenes)})", "h_max": H_MAX,
        "semilla": SEMILLA, "python": platform.python_version(), "pandas": pd.__version__, "numpy": np.__version__,
        "statsmodels": statsmodels.__version__, "xgboost": xgb.__version__}
(RES / "reproducibilidad.json").write_text(json.dumps(meta, indent=2, ensure_ascii=False), encoding="utf-8")

# %% Pronóstico hacia adelante (todos los datos) y paquete para la web
# La programación de vuelos de jul–dic 2026 no está en datos abiertos: hacia adelante los modelos con
# capacidad usan la capacidad estimada (tendencia de los tres últimos meses), también en la variante «programada».
FUT = pd.date_range(MAS(ULTIMO, 1), periods=H_MAX, freq="MS")
prev_f, _ = todos_los_pronosticos(Y, ULTIMO, modos=["estimada"])
P["log_ratio"] = np.log(P["observado"] / P["pronostico"])
Q = P.groupby(["aeropuerto", "modelo", "h"])["log_ratio"].quantile([0.1, 0.9]).unstack()
adelante, cap_futura = {}, {}
for a in Y.columns:
    adelante[a] = {"fechas": [f"{d:%Y-%m}" for d in FUT]}
    for m in MODELOS:
        m_real = m.replace("(programada)", "(estimada)")
        v = prev_f[a][m_real]
        q = Q.loc[(a, m_real)].reindex(range(1, H_MAX + 1))
        adelante[a][m] = {"punto": np.round(v).tolist(),
                          "inf80": np.round(v * np.exp(q[0.1].values)).tolist(),
                          "sup80": np.round(v * np.exp(q[0.9].values)).tolist()}
    cap_futura[a] = [round(cap_hat(a, ULTIMO, t, "estimada")) for t in FUT]

completa = (s.pivot(index="fecha", columns="aeropuerto", values="total").asfreq("MS").fillna(0.0))
web = {
    "generado": pd.Timestamp.now().strftime("%Y-%m-%d %H:%M"),
    "meta": meta,
    "aeropuertos": {a: NOMBRES[a] for a in ["SMR", "RCH"]},
    "modelos_base": BASE, "modelos_cap": CAPM, "modos": MODOS,
    "observado": {a: {"fechas": [f"{d:%Y-%m}" for d in completa.index], "pax": completa[a].round(0).tolist(),
                      "sillas": CAPT[a].reindex(completa.index).fillna(0).round(0).tolist()} for a in ["SMR", "RCH"]},
    "capacidad_futura": cap_futura,
    "inicio_entrenamiento": f"{INICIO_TRAIN:%Y-%m}", "primer_origen": f"{PRIMER_ORIGEN:%Y-%m}",
    "prueba": {a: {m: {str(h): [[f"{r.objetivo:%Y-%m}", round(r.pronostico), round(r.observado), f"{r.origen:%Y-%m}"]
                                for r in P[(P.aeropuerto == a) & (P.modelo == m) & (P.h == h)].sort_values("objetivo").itertuples()]
                       for h in range(1, H_MAX + 1)} for m in MODELOS} for a in ["SMR", "RCH"]},
    "errores": tabla.round(3).to_dict(orient="records"),
    "errores_h": tab_h.round(3).to_dict(orient="records"),
    "dm": dmt.round(3).replace({np.nan: None}).to_dict(orient="records"),
    "shap": {k: [{"variable": ETIQ[f], "valor": round(float(v), 4)} for f, v in s_.sort_values(ascending=False).items()]
             for k, s_ in imp_por.items()},
    "adelante": adelante,
}
(RES / "web_data.json").write_text(json.dumps(web, ensure_ascii=False), encoding="utf-8")

# %% Informe automático (resultados/informe_P1.md)
def md(df):
    cols = list(df.columns)
    lineas = ["| " + " | ".join(cols) + " |", "|" + "---|" * len(cols)]
    for _, r in df.iterrows():
        lineas.append("| " + " | ".join(f"{v:.2f}" if isinstance(v, float) else str(v) for v in r.values) + " |")
    return "\n".join(lineas)

SIN_FUTURO = BASE[1:] + [f"{m} (estimada)" for m in CAPM]
CON_PROG = [f"{m} (programada)" for m in CAPM]

def mejor(a, bloque, candidatos):
    t = tabla[(tabla.aeropuerto == a) & (tabla.bloque == bloque)]
    ref = t[t.modelo == "Ingenuo estacional"].iloc[0]
    top = t[t.modelo.isin(candidatos)].sort_values("MASE").iloc[0]
    return top, ref, top["MASE"] < ref["MASE"]

def linea(a, b, candidatos):
    top, ref, gana = mejor(a, b, candidatos)
    hb = b.replace("h1-3", "1 a 3 meses").replace("h4-6", "4 a 6 meses")
    if gana:
        return (f"- {NOMBRES[a]}, {hb}: **cumple**. {top.modelo}: MAPE {top['MAPE_%']:.1f} % frente a "
                f"{ref['MAPE_%']:.1f} % de la referencia (MASE {top.MASE:.2f} frente a {ref.MASE:.2f}).")
    return (f"- {NOMBRES[a]}, {hb}: **no cumple**. El mejor es {top.modelo} con MAPE {top['MAPE_%']:.1f} % "
            f"frente a {ref['MAPE_%']:.1f} % de la referencia.")

ver_a = [linea(a, b, SIN_FUTURO) for a in ["SMR", "RCH"] for b in ["h1-3", "h4-6"]]
ver_b = [linea(a, b, CON_PROG) for a in ["SMR", "RCH"] for b in ["h1-3", "h4-6"]]
veredictos = ["**A. Sin información del futuro** (solo historia; capacidad proyectada con su tendencia):"] + ver_a +              ["", "**B. Con la capacidad aérea programada** (sillas que las aerolíneas ofrecen en el mes pronosticado):"] + ver_b

t_pub = tabla[tabla.bloque != "h1-6"][["aeropuerto", "bloque", "modelo", "MAPE_%", "MASE", "n"]].copy()
t_pub["n"] = t_pub["n"].astype(int)
d_pub = dmt[dmt.h.isin([3, 6])][["aeropuerto", "h", "modelo", "n", "DM_HLN", "p_valor", "lectura"]]
imp_pub = imp.sort_values(ascending=False).rename(index=ETIQ).reset_index()
imp_pub.columns = ["Variable", "Media |SHAP|"]

informe = f"""# P1 · Pronóstico de demanda — resultados

**Proyecto:** ATARRAYA AZUL · SIGP 143652 · Nota de concepto TRL de entrada, soporte E-13
**Ejecutado por:** CIDS, Universidad del Magdalena · **Código:** `scripts/02_p1_pronostico.py`

## Datos
- Pasajeros: Aerocivil, «Transporte Aéreo Comercial – Tráfico Origen-Destino (Colombia)», datos.gov.co `gb6w-ynu4`. Embarcados + desembarcados por mes en SMR y RCH.
- Capacidad: Aerocivil, «Transporte Aéreo – Tráfico por Equipo (Colombia)», datos.gov.co `djjf-g4q4`. Sillas ofrecidas en vuelos que salen o llegan a cada aeropuerto.
- Consultas registradas en `datos/_procedencia.csv`. Serie usada: {meta['serie']} ({len(Y)} meses); se excluyen 2020–2021 por la pandemia.
- SHA-256 pasajeros: `{sha}` · capacidad: `{sha_cap}`.

## Diseño
- Evaluación fuera de muestra con origen móvil: {meta['origenes']} orígenes. En cada uno se reentrena con datos hasta ese mes y se pronostican h = 1…{H_MAX} meses.
- Referencia: ingenuo estacional, ŷ(t) = y(t−12).
- Modelos sin capacidad: ETS (log, tendencia amortiguada), SARIMA(0,1,1)(0,1,1)₁₂ sobre log, XGBoost global directo sobre el crecimiento interanual en log.
- Modelos con capacidad: Capacidad × ocupación (sillas del mes × ocupación del mismo mes del año anterior); SARIMAX del log de pasajeros sobre el log de sillas con errores AR(1)(1)₁₂; XGBoost con el crecimiento interanual de las sillas como señal adicional.
- Dos variantes de la capacidad del mes pronosticado: **programada** (sillas efectivamente ofrecidas, como aproximación a la programación que las aerolíneas publican con meses de anticipación) y **estimada** (sin información futura: mismo mes del año anterior × crecimiento reciente de la capacidad).
- Métricas: MAPE y MASE (escala: MAE del ingenuo estacional en 2023–2024). Diebold-Mariano con corrección HLN.
- Criterio (nota de concepto, 7.3): error fuera de muestra inferior al de la referencia.

## Veredicto
{chr(10).join(veredictos)}

**Lectura.** La historia de pasajeros por sí sola apenas mejora la referencia, y solo en Santa Marta a corto plazo. Cuando se conoce la capacidad aérea del mes pronosticado, el error baja a la mitad en Santa Marta y la referencia queda superada en los cuatro casos. La capacidad proyectada con su propia tendencia no sirve. El valor está en integrar la programación real de las aerolíneas como fuente, que es lo que hace el espacio de datos del proyecto (A3 y A4).

## Errores fuera de muestra
{md(t_pub.round(2))}

## Prueba de Diebold-Mariano frente a la referencia
Valores negativos: el modelo tiene menor error cuadrático que la referencia.

{md(d_pub.round(3))}

## Explicabilidad (SHAP, XGBoost + capacidad programada)
{md(imp_pub.round(4))}

![Pronósticos a 3 meses](fig_pronostico_h3.png)

![SHAP](fig_shap.png)

## Limitaciones declaradas
1. La serie oficial abierta inicia en 2020; descontada la pandemia quedan {len(Y)} meses. Con 13 a 18 orígenes ninguna diferencia alcanza significancia estadística al 5 %.
2. La variante «programada» usa las sillas efectivamente ofrecidas. La programación publicada por las aerolíneas con anticipación puede diferir por cancelaciones o cambios de equipo, así que este resultado es una cota superior. La prueba con programación real queda para el TRL 4, cuando el espacio de datos integre esa fuente.
3. Santa Marta perdió cerca del 30 % de su tráfico entre marzo y noviembre de 2023 por una caída de la capacidad (la ocupación se mantuvo). Los modelos con capacidad lo capturan; los univariados no.
4. Las cifras recientes de Aerocivil son provisionales; el último mes publicado es {ULTIMO:%m-%Y}.
5. Es una prueba de un componente aislado en laboratorio (TRL 3); no valida el sistema integrado.

## Reproducir
```
python scripts/01_descargar_datosgov.py
python scripts/02_p1_pronostico.py
python scripts/03_construir_web.py
```
Versiones: Python {meta['python']}, pandas {meta['pandas']}, numpy {meta['numpy']}, statsmodels {meta['statsmodels']}, xgboost {meta['xgboost']}; semilla {SEMILLA}.
"""
(RES / "informe_P1.md").write_text(informe, encoding="utf-8")

pd.set_option("display.width", 200)
print("\n".join(veredictos))
print(tabla.round(2).to_string(index=False))
print(dmt.round(3).to_string(index=False))
print(imp.round(4).to_string())

DescargarSHA-256 84cdfb1f46502343c0d6a8196544d96001d03137484743472f2eecdcae7868b7

# %% [markdown]
# # P2 · Prueba mínima reproducible del umbral de saturación (E-13)
# ATARRAYA AZUL · Nota de concepto TRL de entrada · CIDS, Universidad del Magdalena
#
# **Pregunta.** ¿Una regla simple de alerta (ocupación media diaria estimada ≥ umbral × capacidad de
# carga) se dispara en los mismos meses en que hay alta afluencia documentada en el PNN Tayrona?
#
# **Limitación de partida.** Parques Nacionales no publica la serie mensual ni diaria del Tayrona;
# solo totales anuales. Esta versión reparte el total anual por meses con un indicador indirecto
# (pasajeros aéreos de Santa Marta) y lo declara. El código acepta la serie diaria de taquilla
# (datos/pnn_tayrona_diario.csv con columnas fecha,visitantes) y entonces la usa directamente.
#
# **Diseño.** Calibración del umbral con 2022–2023; validación con 2024 (fuera de muestra).

# %% Configuración
import calendar, hashlib, json, sys
from pathlib import Path
import numpy as np, pandas as pd
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt

sys.stdout.reconfigure(encoding="utf-8")
RAIZ = Path(__file__).resolve().parents[1]
DATOS, RES = RAIZ / "datos", RAIZ / "resultados"

# Visitantes anuales registrados en el PNN Tayrona (ingresos terrestres)
ANUAL = {
    2022: (545_760, "PNN, Plan de Necesidades PNN Tayrona (2024), gráfico 1"),
    2023: (637_534, "PNN, Informe anual de visitantes 2023"),
    2024: (615_599, "PNN, Informe anual de visitantes 2024"),
    2025: (830_326, "PNN, Informe anual de visitantes 2025 (incluye ingresos marino-costeros desde el 2.º semestre; no comparable)"),
}
CAPACIDADES = {
    6900: "Capacidad de carga comunicada por Parques Nacionales («hasta 6.900 visitantes diarios»)",
    5876: "Capacidad de carga real, ejercicio 2015, Plan de Ordenamiento Ecoturístico PNN Tayrona (tabla 30)",
}
# Cierres anuales acordados con los cuatro pueblos de la Sierra Nevada (Informe PNN 2025)
CIERRES = [((2, 1), (2, 15)), ((6, 1), (6, 15)), ((10, 19), (11, 2))]
# Domingo de Ramos a Domingo de Resurrección
SEMANA_SANTA = {2022: ("2022-04-10", "2022-04-17"), 2023: ("2023-04-02", "2023-04-09"),
                2024: ("2024-03-24", "2024-03-31"), 2025: ("2025-04-13", "2025-04-20")}
ANIOS_CAL, ANIO_VAL = [2022, 2023], 2024

# %% Días abiertos y temporadas altas documentadas por día
def dias_del_anio(y):
    d = pd.DataFrame({"fecha": pd.date_range(f"{y}-01-01", f"{y}-12-31", freq="D")})
    md = list(zip(d.fecha.dt.month, d.fecha.dt.day))
    cerrado = np.zeros(len(d), bool)
    for (m1, d1), (m2, d2) in CIERRES:
        cerrado |= np.array([(m1, d1) <= x <= (m2, d2) for x in md])
    d["abierto"] = ~cerrado
    ss0, ss1 = SEMANA_SANTA[y]
    alta = ((d.fecha.dt.month == 12) & (d.fecha.dt.day >= 15)) | ((d.fecha.dt.month == 1) & (d.fecha.dt.day <= 15)) \
         | ((d.fecha.dt.month == 6) & (d.fecha.dt.day >= 15)) | ((d.fecha.dt.month == 7) & (d.fecha.dt.day <= 15)) \
         | d.fecha.between(ss0, ss1)
    d["temporada_alta"] = alta & d["abierto"]
    return d

# %% Perfil mensual: pasajeros aéreos de Santa Marta (indicador indirecto)
pax = pd.read_csv(DATOS / "serie_pasajeros_datosgov.csv", parse_dates=["fecha"])
smr = pax[pax.aeropuerto == "SMR"].set_index("fecha")["total"]
diario = DATOS / "pnn_tayrona_diario.csv"
usa_diario = diario.exists()

filas = []
for y, (total, fuente) in ANUAL.items():
    d = dias_del_anio(y)
    if usa_diario:
        obs = pd.read_csv(diario, parse_dates=["fecha"]).set_index("fecha")["visitantes"]
        d["visitantes"] = d.fecha.map(obs).fillna(0)
        mens = d.groupby(d.fecha.dt.month)["visitantes"].sum()
    g = d.groupby(d.fecha.dt.month).agg(dias_abiertos=("abierto", "sum"), dias_alta=("temporada_alta", "sum"),
                                        dias=("abierto", "size"))
    if not usa_diario:
        # Intensidad diaria del indicador indirecto (pasajeros/día del mes) × días abiertos: el parque no recibe
        # visitantes en los días de cierre. Se normaliza para que la suma anual sea el total oficial.
        s = smr[smr.index.year == y]
        intensidad = pd.Series(s.values, index=s.index.month) / g["dias"]
        peso = intensidad * g["dias_abiertos"]
        mens = total * peso / peso.sum()
    for m in range(1, 13):
        filas.append({"anio": y, "mes": m, "visitantes_mes": float(mens[m]), "dias_abiertos": int(g.at[m, "dias_abiertos"]),
                      "dias_temporada_alta": int(g.at[m, "dias_alta"]), "fuente_total": fuente})
M = pd.DataFrame(filas)
M["media_diaria"] = M["visitantes_mes"] / M["dias_abiertos"].replace(0, np.nan)
M["alta_documentada"] = M["dias_temporada_alta"] >= 7     # el mes contiene al menos una semana de temporada alta
M.to_csv(RES / "p2_mensual.csv", index=False)

# %% Regla de alerta y métricas
def evaluar(df, cap, umbral):
    a = (df["media_diaria"] / cap) >= umbral
    doc = df["alta_documentada"]
    vp, fp, fn = int((a & doc).sum()), int((a & ~doc).sum()), int((~a & doc).sum())
    prec = vp / (vp + fp) if vp + fp else np.nan
    rec = vp / (vp + fn) if vp + fn else np.nan
    f1 = 2 * prec * rec / (prec + rec) if prec and rec and not np.isnan(prec + rec) else 0.0
    return {"aciertos": vp, "falsas_alarmas": fp, "perdidas": fn, "precision": prec, "sensibilidad": rec, "F1": f1}

UMBRALES = np.round(np.arange(0.20, 0.71, 0.01), 2)
cal = M[M.anio.isin(ANIOS_CAL)]
val = M[M.anio == ANIO_VAL]
resultados = {}
for cap in CAPACIDADES:
    barrido = pd.DataFrame([{"umbral": u, **evaluar(cal, cap, u)} for u in UMBRALES])
    mejor = barrido.sort_values(["F1", "umbral"], ascending=[False, True]).iloc[0]
    u = float(mejor.umbral)
    r_val = evaluar(val, cap, u)
    # referencia sin información: alertar siempre en los meses de temporada fija (dic, ene, jul)
    ref = val["mes"].isin([12, 1, 7])
    doc = val["alta_documentada"]
    r_ref = {"aciertos": int((ref & doc).sum()), "falsas_alarmas": int((ref & ~doc).sum()), "perdidas": int((~ref & doc).sum())}
    resultados[cap] = {"umbral": u, "calibracion": {k: (None if pd.isna(v) else float(v)) for k, v in mejor.items() if k != "umbral"},
                       "validacion": {k: (None if pd.isna(v) else float(v)) for k, v in r_val.items()},
                       "barrido": barrido.round(3).replace({np.nan: None}).to_dict(orient="records"), "referencia_calendario": r_ref}
    barrido.to_csv(RES / f"p2_barrido_cap{cap}.csv", index=False)

# %% Gráfico
cap0 = 6900; u0 = resultados[cap0]["umbral"]
fig, ax = plt.subplots(figsize=(10, 4))
for i, y in enumerate([2022, 2023, 2024]):
    q = M[M.anio == y]
    x = np.arange(12) + i * 12
    ax.bar(x, q.media_diaria / cap0 * 100, color=np.where(q.alta_documentada, "#0f8b8d", "#b9c7cf"), width=0.8)
    ax.text(i * 12 + 5.5, 2, str(y), ha="center", fontsize=9, color="#13232e")
ax.axhline(u0 * 100, color="#b4462f", lw=1.5, ls="--", label=f"Umbral calibrado {u0*100:.0f} %")
ax.set_ylabel("Ocupación media diaria (% de 6.900)"); ax.set_xticks([])
ax.set_title("P2 · Ocupación media estimada del PNN Tayrona por mes (verde: temporada alta documentada)")
ax.legend(loc="upper right"); fig.tight_layout(); fig.savefig(RES / "fig_p2.png", dpi=160); plt.close(fig)

# %% P2-B · Saturación del acceso aéreo (datos mensuales reales de Aerocivil)
# Ocupación de sillas = pasajeros (embarcados + desembarcados) / sillas ofrecidas (salidas + llegadas).
SEMANA_SANTA[2026] = ("2026-03-29", "2026-04-05")
JUEVES_SANTO = {2021: "2021-04-01", 2022: "2022-04-14", 2023: "2023-04-06", 2024: "2024-03-28", 2025: "2025-04-17", 2026: "2026-04-02"}
cap_air = pd.read_csv(DATOS / "serie_capacidad_datosgov.csv", parse_dates=["fecha"])
A = pax.merge(cap_air, on=["aeropuerto", "fecha"])
A = A[A.fecha >= "2021-01-01"].sort_values(["aeropuerto", "fecha"]).copy()
A["ocupacion"] = A["total"] / A["sillas"]
A["media_12m_previa"] = A.groupby("aeropuerto")["ocupacion"].transform(lambda s: s.shift(1).rolling(12).mean())
A["exceso"] = A["ocupacion"] - A["media_12m_previa"]
def alta_doc_mes(f):
    y, m = f.year, f.month
    dias = 31 if m in (1, 7, 12) else 30
    alta = {1: 15, 12: 17, 6: 16, 7: 15}.get(m, 0)
    return alta >= 7 or (pd.Timestamp(JUEVES_SANTO[y]).month == m)
A["alta_documentada"] = A["fecha"].apply(alta_doc_mes)
A = A[A.fecha >= "2022-01-01"]

def eval_regla(df, col, u):
    a = df[col] >= u; doc = df["alta_documentada"]
    vp, fp, fn = int((a & doc).sum()), int((a & ~doc).sum()), int((~a & doc).sum())
    prec = vp / (vp + fp) if vp + fp else np.nan; rec = vp / (vp + fn) if vp + fn else np.nan
    f1 = 2 * prec * rec / (prec + rec) if (vp and not np.isnan(prec + rec)) else 0.0
    return {"aciertos": vp, "falsas_alarmas": fp, "perdidas": fn, "precision": prec, "sensibilidad": rec, "F1": f1}

REGLAS = {"absoluta": ("ocupacion", np.round(np.arange(0.70, 0.98, 0.005), 3)),
          "relativa": ("exceso", np.round(np.arange(-0.05, 0.12, 0.005), 3))}
cal_b, val_b = A[A.fecha.dt.year.isin(ANIOS_CAL)], A[A.fecha.dt.year >= 2024]
acceso = {}
for ap in ["SMR", "RCH"]:
    acceso[ap] = {}
    for regla, (col, grid) in REGLAS.items():
        c_ = cal_b[cal_b.aeropuerto == ap]; v_ = val_b[val_b.aeropuerto == ap]
        barr = pd.DataFrame([{"umbral": u, **eval_regla(c_, col, u)} for u in grid])
        best = barr.sort_values(["F1", "umbral"], ascending=[False, False]).iloc[0]
        u = float(best.umbral); rv = eval_regla(v_, col, u)
        ref = v_.fecha.dt.month.isin([12, 1, 7]); doc = v_["alta_documentada"]
        rr = eval_regla(v_.assign(_r=ref.astype(float)), "_r", 0.5)
        acceso[ap][regla] = {"umbral": u, "calibracion": {k: (None if pd.isna(x) else float(x)) for k, x in best.items() if k != "umbral"},
                             "validacion": {k: (None if pd.isna(x) else float(x)) for k, x in rv.items()},
                             "referencia_calendario": {k: (None if pd.isna(x) else float(x)) for k, x in rr.items()}}
A.to_csv(RES / "p2_acceso_aereo.csv", index=False)
for ap in acceso:
    for regla, r in acceso[ap].items():
        print(ap, regla, "umbral", r["umbral"], "cal", {k: round(v, 2) for k, v in r["calibracion"].items()},
              "| val", {k: round(v, 2) if v is not None else None for k, v in r["validacion"].items()},
              "| ref", {k: round(v, 2) if v is not None else None for k, v in r["referencia_calendario"].items()})

# %% P2-C · Saturación del alojamiento (DANE, Encuesta Mensual de Alojamiento, cuadro 4.2)
import openpyxl, requests
URL_EMA = "https://www.dane.gov.co/files/operaciones/EMA/anex-EMA-ago2026.xlsx"
ema_arch = DATOS / "dane_raw" / Path(URL_EMA).name
if not ema_arch.exists():
    ema_arch.parent.mkdir(exist_ok=True)
    ema_arch.write_bytes(requests.get(URL_EMA, headers={"User-Agent": "Mozilla/5.0"}, timeout=120).content)
sha_ema = hashlib.sha256(ema_arch.read_bytes()).hexdigest()
REGIONES_EMA = {"HCA": ("Caribe", "Hoteles · Región Costa Caribe (Atlántico, Cesar, La Guajira y Magdalena)"),
                "HGM": ("Golfo Morrosquillo y Sabana", "Hoteles · Región Golfo de Morrosquillo y Sabana (Bolívar, Córdoba y Sucre)")}
MESES_ES = ["enero", "febrero", "marzo", "abril", "mayo", "junio", "julio", "agosto", "septiembre", "octubre", "noviembre", "diciembre"]
ws = openpyxl.load_workbook(ema_arch, read_only=True, data_only=True)["4.2 Porc Mens Ocupación.reg"]
filas_ema, cab, anio = [], None, None
for r in ws.iter_rows(values_only=True):
    if r and r[0] == "Año":
        cab = [str(c).strip() if c else "" for c in r]; continue
    if cab and r and r[1] and str(r[1]).strip().lower() in MESES_ES:
        if r[0]:
            anio = int(str(r[0])[:4])
        for cod, (reg, _) in REGIONES_EMA.items():
            filas_ema.append({"region": cod, "fecha": pd.Timestamp(anio, MESES_ES.index(str(r[1]).strip().lower()) + 1, 1),
                              "ocupacion": float(r[cab.index(reg)]) / 100})
H = pd.DataFrame(filas_ema).sort_values(["region", "fecha"])
H["media_12m_previa"] = H.groupby("region")["ocupacion"].transform(lambda s: s.shift(1).rolling(12).mean())
H["exceso"] = H["ocupacion"] - H["media_12m_previa"]
JUEVES_SANTO.update({2019: "2019-04-18", 2020: "2020-04-09"})
H["alta_documentada"] = H["fecha"].apply(alta_doc_mes)
H = H[H.fecha >= "2022-01-01"]
REGLAS_H = {"absoluta": ("ocupacion", np.round(np.arange(0.30, 0.75, 0.005), 3)),
            "relativa": ("exceso", np.round(np.arange(-0.05, 0.12, 0.005), 3))}
hotel = {}
for cod in REGIONES_EMA:
    hotel[cod] = {}
    c_ = H[(H.region == cod) & H.fecha.dt.year.isin(ANIOS_CAL)]; v_ = H[(H.region == cod) & (H.fecha.dt.year >= 2024)]
    for regla, (col, grid) in REGLAS_H.items():
        barr = pd.DataFrame([{"umbral": u, **eval_regla(c_, col, u)} for u in grid])
        best = barr.sort_values(["F1", "umbral"], ascending=[False, False]).iloc[0]
        u = float(best.umbral); rv = eval_regla(v_, col, u)
        rr = eval_regla(v_.assign(_r=v_.fecha.dt.month.isin([12, 1, 7]).astype(float)), "_r", 0.5)
        hotel[cod][regla] = {"umbral": u, "calibracion": {k: (None if pd.isna(x) else float(x)) for k, x in best.items() if k != "umbral"},
                             "validacion": {k: (None if pd.isna(x) else float(x)) for k, x in rv.items()},
                             "referencia_calendario": {k: (None if pd.isna(x) else float(x)) for k, x in rr.items()}}
        print(cod, regla, "umbral", u, "| val", {k: round(x, 2) if x is not None else None for k, x in hotel[cod][regla]["validacion"].items()},
              "| ref", {k: round(x, 2) if x is not None else None for k, x in hotel[cod][regla]["referencia_calendario"].items()})
H.to_csv(RES / "p2_ocupacion_hotelera.csv", index=False)

# %% P2-D · Alerta anticipada de saturación del acceso aéreo (P1 → P2)
# Se usan los pronósticos a h meses de la prueba P1 (origen móvil, sin ver el futuro de pasajeros).
# Ocupación pronosticada = pasajeros pronosticados / sillas del mes. Umbral fijado de antemano: 85 %.
# Evento de referencia: meses en que la ocupación MEDIDA por Aerocivil fue ≥ umbral (no es un calendario).
U_ANT, H_ANT = 0.85, 3
MOD_ANT = ["Ingenuo estacional", "ETS", "Capacidad × ocupación (programada)", "SARIMAX + capacidad (programada)",
           "SARIMAX + capacidad (estimada)"]
PB = pd.read_csv(RES / "pronosticos_origen_movil.csv", parse_dates=["objetivo", "origen"])
PB = PB.merge(cap_air.rename(columns={"fecha": "objetivo"})[["aeropuerto", "objetivo", "sillas"]], on=["aeropuerto", "objetivo"])
PB["ocup_obs"], PB["ocup_pron"] = PB.observado / PB.sillas, PB.pronostico / PB.sillas
anticipada = {"umbral": U_ANT, "h": H_ANT, "modelos": MOD_ANT, "serie": {}, "metricas": {}, "sensibilidad": {}}
for ap in ["SMR", "RCH"]:
    q = PB[(PB.aeropuerto == ap) & (PB.h == H_ANT)]
    anticipada["serie"][ap] = {m: {"fechas": q[q.modelo == m].sort_values("objetivo").objetivo.dt.strftime("%Y-%m").tolist(),
                                   "pron": q[q.modelo == m].sort_values("objetivo").ocup_pron.round(4).tolist(),
                                   "obs": q[q.modelo == m].sort_values("objetivo").ocup_obs.round(4).tolist()} for m in MOD_ANT}
    anticipada["metricas"][ap], anticipada["sensibilidad"][ap] = {}, {}
    for m in MOD_ANT:
        z = q[q.modelo == m].assign(alerta=lambda d: d.ocup_pron, evento=lambda d: (d.ocup_obs >= U_ANT))
        r = eval_regla(z.rename(columns={"evento": "alta_documentada"}), "ocup_pron", U_ANT)
        r["exactitud"] = float(((z.ocup_pron >= U_ANT) == (z.ocup_obs >= U_ANT)).mean())
        anticipada["metricas"][ap][m] = {k: (None if pd.isna(x) else float(x)) for k, x in r.items()}
        sens = []
        for u in (0.80, 0.83, 0.85, 0.87):
            zz = q[q.modelo == m].assign(alta_documentada=lambda d: d.ocup_obs >= u)
            rr = eval_regla(zz, "ocup_pron", u)
            sens.append({"umbral": u, "eventos": int(zz.alta_documentada.sum()), **{k: (None if pd.isna(x) else float(x)) for k, x in rr.items()}})
        anticipada["sensibilidad"][ap][m] = sens
    print(ap, {m: (v["aciertos"], v["falsas_alarmas"], v["perdidas"], round(v["exactitud"], 2)) for m, v in anticipada["metricas"][ap].items()})

# %% Paquete web e informe
web = {"meses": M.round(1).replace({np.nan: None}).to_dict(orient="records"),
       "acceso": {"serie": {ap: {"fechas": A[A.aeropuerto == ap].fecha.dt.strftime("%Y-%m").tolist(),
                                 "ocupacion": A[A.aeropuerto == ap].ocupacion.round(4).tolist(),
                                 "media12": A[A.aeropuerto == ap].media_12m_previa.round(4).tolist(),
                                 "pax": A[A.aeropuerto == ap].total.round(0).tolist(),
                                 "sillas": A[A.aeropuerto == ap].sillas.round(0).tolist(),
                                 "alta": A[A.aeropuerto == ap].alta_documentada.tolist()} for ap in ["SMR", "RCH"]},
                  "resultados": acceso},
       "anticipada": anticipada,
       "hotel": {"serie": {cod: {"fechas": H[H.region == cod].fecha.dt.strftime("%Y-%m").tolist(),
                                 "ocupacion": H[H.region == cod].ocupacion.round(4).tolist(),
                                 "alta": H[H.region == cod].alta_documentada.tolist()} for cod in REGIONES_EMA},
                 "nombres": {cod: v[1] for cod, v in REGIONES_EMA.items()}, "resultados": hotel, "sha256": sha_ema,
                 "fuente": URL_EMA},
       "capacidades": {str(k): v for k, v in CAPACIDADES.items()},
       "anios_calibracion": ANIOS_CAL, "anio_validacion": ANIO_VAL, "usa_diario": usa_diario,
       "resultados": {str(k): v for k, v in resultados.items()},
       "cierres": "1–15 feb, 1–15 jun, 19 oct–2 nov",
       "semana_santa": SEMANA_SANTA,
       "sha256_pasajeros": hashlib.sha256((DATOS / "serie_pasajeros_datosgov.csv").read_bytes()).hexdigest()}
(RES / "web_p2.json").write_text(json.dumps(web, ensure_ascii=False), encoding="utf-8")

lin = []
for cap, r in resultados.items():
    c, v, rf = r["calibracion"], r["validacion"], r["referencia_calendario"]
    lin.append(f"| {cap:,} | {r['umbral']*100:.0f} % | {c['aciertos']:.0f} / {c['falsas_alarmas']:.0f} / {c['perdidas']:.0f} | "
               f"{v['aciertos']:.0f} / {v['falsas_alarmas']:.0f} / {v['perdidas']:.0f} | "
               f"{(v['precision'] or 0)*100:.0f} % | {(v['sensibilidad'] or 0)*100:.0f} % | "
               f"{rf['aciertos']} / {rf['falsas_alarmas']} / {rf['perdidas']} |")
informe = f"""# P2 · Umbral de saturación — resultados

**Proyecto:** ATARRAYA AZUL · SIGP 143652 · Nota de concepto TRL de entrada, soporte E-13
**Código:** `scripts/04_p2_saturacion.py` · Sitio piloto: PNN Tayrona (Santa Marta)

## Datos
- Visitantes anuales: {", ".join(f"{y}: {v[0]:,} ({v[1]})" for y, v in ANUAL.items())}.
- Capacidad de carga: {"; ".join(f"{k:,}/día — {v}" for k, v in CAPACIDADES.items())}.
- Días abiertos: se descuentan los tres cierres anuales ({web['cierres']}).
- Reparto mensual: {"serie diaria de taquilla (datos/pnn_tayrona_diario.csv)" if usa_diario else "Parques Nacionales no publica la serie mensual; el total anual se reparte con la estacionalidad de los pasajeros aéreos de Santa Marta (Aerocivil), como indicador indirecto"}.
- Alta afluencia documentada (referencia independiente de la serie): meses con al menos 7 días de temporada alta abiertos (15 dic–15 ene, Semana Santa, 15 jun–15 jul).

## Regla y protocolo
1. Ocupación media diaria del mes = visitantes del mes / días abiertos.
2. Alerta si la ocupación media diaria ≥ umbral × capacidad de carga.
3. El umbral se calibra con {ANIOS_CAL[0]}–{ANIOS_CAL[1]} (máximo F1) y se valida, sin cambiarlo, con {ANIO_VAL}.
4. Se reporta aciertos / falsas alarmas / alertas perdidas, precisión y sensibilidad, frente a una referencia de calendario fijo (alertar siempre en dic, ene y jul).

## Resultados
| Capacidad (visit./día) | Umbral calibrado | Calibración A/FA/P | Validación {ANIO_VAL} A/FA/P | Precisión {ANIO_VAL} | Sensibilidad {ANIO_VAL} | Referencia calendario A/FA/P |
|---|---|---|---|---|---|---|
{chr(10).join(lin)}

![P2](fig_p2.png)

## Limitaciones declaradas
1. Sin la serie mensual o diaria oficial, el reparto mensual usa un indicador indirecto (pasajeros aéreos), que no capta excursionistas que llegan por tierra ni los picos de fin de semana. Por eso la regla trabaja con ocupación media, no con el día pico.
2. La ocupación media mensual queda muy por debajo de la capacidad (6.900/día); la saturación real ocurre en días puntuales. La prueba valida la regla y su protocolo, no el nivel de riesgo diario.
3. El total de 2025 no es comparable (incluye ingresos marino-costeros desde el segundo semestre) y no se usa para validar.
4. Las fechas de temporada alta y de cierres deben confirmarse con las resoluciones de Parques Nacionales de cada año.
5. Siguiente paso: solicitar a Parques Nacionales la serie diaria de taquilla por sector (Zaíno, Palangana, Bahía Concha, Calabazo). El script la usa sin cambios si se guarda en `datos/pnn_tayrona_diario.csv`.

## P2-B · Saturación del acceso aéreo (datos reales)
Misma regla aplicada al canal de llegada: ocupación de sillas = pasajeros / sillas ofrecidas por mes (Aerocivil, gb6w-ynu4 y djjf-g4q4). Calibración 2022–2023; validación enero 2024 a junio 2026.

| Aeropuerto | Regla | Umbral | Validación A/FA/P | Precisión | Sensibilidad | Calendario A/FA/P |
|---|---|---|---|---|---|---|
{{tabla_b}}

**Lectura.** La regla detecta la mayoría de los meses de temporada alta, con falsas alarmas. No supera a un calendario fijo porque la referencia de «alta afluencia documentada» es en sí misma un calendario. Lo que la regla aporta y el calendario no, es detectar saturaciones fuera de temporada: en 2023 la ocupación de sillas de Santa Marta se mantuvo en 85–90 % de marzo a noviembre tras la reducción de capacidad aérea. Con datos diarios del Tayrona la comparación pasaría a hacerse contra días de cierre por aforo, que no dependen del calendario.

## P2-C · Saturación del alojamiento (DANE, Encuesta Mensual de Alojamiento)
Porcentaje de ocupación hotelera mensual (cuadro 4.2) de las regiones Costa Caribe (Atlántico, Cesar, La Guajira y Magdalena) y Golfo de Morrosquillo y Sabana (Bolívar, Córdoba y Sucre). Calibración 2022–2023; validación enero 2024 a agosto 2026.

| Región | Regla | Umbral | Validación A/FA/P | Precisión | Sensibilidad | Calendario A/FA/P |
|---|---|---|---|---|---|---|
{{tabla_c}}

{{seccion_d}}

**Lectura C.** La ocupación hotelera regional casi no tiene estacionalidad turística: en la Costa Caribe pesan los viajes de negocios de Barranquilla y Valledupar, y en el Golfo el umbral calibrado queda tan bajo que alerta casi todos los meses. La agregación regional diluye la señal de los destinos piloto; serviría con la ocupación por municipio, que el DANE no publica.
"""
filas_b = []
for ap in acceso:
    for regla, r in acceso[ap].items():
        v, rf = r["validacion"], r["referencia_calendario"]
        u = f"{r['umbral']*100:.1f} %" if regla == "absoluta" else f"{r['umbral']*100:+.1f} pp sobre media 12 m"
        filas_b.append(f"| {ap} | {regla} | {u} | {v['aciertos']:.0f}/{v['falsas_alarmas']:.0f}/{v['perdidas']:.0f} | "
                       f"{(v['precision'] or 0)*100:.0f} % | {(v['sensibilidad'] or 0)*100:.0f} % | "
                       f"{rf['aciertos']:.0f}/{rf['falsas_alarmas']:.0f}/{rf['perdidas']:.0f} |")
informe = informe.replace("{tabla_b}", chr(10).join(filas_b))
filas_c = []
for cod in hotel:
    for regla, r in hotel[cod].items():
        v, rf = r["validacion"], r["referencia_calendario"]
        u = f"{r['umbral']*100:.1f} %" if regla == "absoluta" else f"{r['umbral']*100:+.1f} pp sobre media 12 m"
        filas_c.append(f"| {REGIONES_EMA[cod][0]} | {regla} | {u} | {v['aciertos']:.0f}/{v['falsas_alarmas']:.0f}/{v['perdidas']:.0f} | "
                       f"{(v['precision'] or 0)*100:.0f} % | {(v['sensibilidad'] or 0)*100:.0f} % | "
                       f"{rf['aciertos']:.0f}/{rf['falsas_alarmas']:.0f}/{rf['perdidas']:.0f} |")
informe = informe.replace("{tabla_c}", chr(10).join(filas_c))
filas_d = []
for ap in ["SMR", "RCH"]:
    for m in MOD_ANT:
        v = anticipada["metricas"][ap][m]
        filas_d.append(f"| {ap} | {m} | {v['aciertos']:.0f}/{v['falsas_alarmas']:.0f}/{v['perdidas']:.0f} | "
                       f"{'—' if v['precision'] is None else f'{v[chr(112)+chr(114)+chr(101)+chr(99)+chr(105)+chr(115)+chr(105)+chr(111)+chr(110)]*100:.0f} %'} | "
                       f"{'—' if v['sensibilidad'] is None else f'{v[chr(115)+chr(101)+chr(110)+chr(115)+chr(105)+chr(98)+chr(105)+chr(108)+chr(105)+chr(100)+chr(97)+chr(100)]*100:.0f} %'} | {v['exactitud']*100:.0f} % |")
sens_smr = []
for m in ["Ingenuo estacional", "ETS", "SARIMAX + capacidad (programada)"]:
    s_ = anticipada["sensibilidad"]["SMR"][m]
    sens_smr.append(f"| {m} | " + " | ".join(f"{x['aciertos']:.0f}/{x['falsas_alarmas']:.0f}/{x['perdidas']:.0f} ({x['eventos']} ev.)" for x in s_) + " |")
seccion_d = f"""## P2-D · Alerta anticipada de saturación del acceso aéreo (P1 → P2)
Pregunta: con el pronóstico de pasajeros a {H_ANT} meses de la prueba P1, ¿se anticipan los meses en que la ocupación de sillas medida por Aerocivil superó el {U_ANT*100:.0f} %? El umbral se fija de antemano y no se calibra; el evento de referencia es una medición oficial, no un calendario. Ventana: los mismos orígenes móviles de P1 (pronósticos para mar-2025 a jun-2026).

| Aeropuerto | Modelo | Aciertos/falsas alarmas/perdidas | Precisión | Sensibilidad | Exactitud |
|---|---|---|---|---|---|
{chr(10).join(filas_d)}

Sensibilidad al umbral en Santa Marta (aciertos/falsas alarmas/perdidas y número de eventos medidos):

| Modelo | 80 % | 83 % | 85 % | 87 % |
|---|---|---|---|---|
{chr(10).join(sens_smr)}

**Lectura D.** En Santa Marta las alertas anticipadas a 3 meses de ETS y de la referencia ingenua coinciden con cinco o seis de los seis meses de saturación medidos al 85 %, con una exactitud de 75 %; un calendario fijo (dic, ene, jul) acierta solo dos. Es la variante de P2 que cumple el criterio de coincidencia con datos abiertos y la que conecta P1 con P2. Los modelos con capacidad programada, mejores en P1, anticipan menos aquí porque subestiman la ocupación en los meses más llenos. Con 16 meses evaluados ninguna variante es claramente superior y el resultado depende del umbral; se reporta el umbral fijado de antemano. En Riohacha hubo un solo mes con ocupación ≥ 85 %, insuficiente para concluir.
"""
informe = informe.replace("{seccion_d}", seccion_d)
(RES / "informe_P2.md").write_text(informe, encoding="utf-8")
print(informe.split("## Resultados")[1].split("![P2]")[0])
print(M[M.anio == 2024][["mes", "visitantes_mes", "dias_abiertos", "media_diaria", "alta_documentada"]].round(0).to_string(index=False))

DescargarSHA-256 a36a41f50a07f6a8d17aba7f0b847ba9f4bda0f6f9ee92a7330df5ecf476fa02

"""Contexto · Conectividad aérea del Caribe: corredor del proyecto frente a Cartagena y Barranquilla.

Ilustra la situación de partida del numeral 5.2 (Anexo 2): los aeropuertos del corredor (Santa Marta,
Riohacha, Corozal, Tolú y Coveñas) movilizan muchos menos pasajeros y tienen menos destinos directos que
Cartagena y Barranquilla, por donde entra buena parte de la demanda hacia el Caribe. No mide la fuga de
tráfico (eso requiere la encuesta de ruta de acceso, A1/A4); muestra la asimetría de la conectividad.

Fuente: Aerocivil, datos.gov.co gb6w-ynu4 (Origen-Destino) y djjf-g4q4 (Tráfico por Equipo).
Salidas: datos/conectividad_*.csv, resultados/web_contexto.json
"""
import datetime as dt, hashlib, json, sys, time
from pathlib import Path
import pandas as pd, requests

sys.stdout.reconfigure(encoding="utf-8")
RAIZ = Path(__file__).resolve().parents[1]
DATOS, RES = RAIZ / "datos", RAIZ / "resultados"
OD, EQ = "https://www.datos.gov.co/resource/gb6w-ynu4.json", "https://www.datos.gov.co/resource/djjf-g4q4.json"
AEROPUERTOS = {
    "SMR": ("Santa Marta", "corredor"), "RCH": ("Riohacha", "corredor"), "CZU": ("Corozal", "corredor"),
    "TLU": ("Tolú", "corredor"), "CVN": ("Coveñas", "corredor"),
    "CTG": ("Cartagena", "puerta"), "BAQ": ("Barranquilla", "puerta"),
}
procedencia = []

def soda(url, params):
    for intento in range(5):
        r = requests.get(url, params=params, timeout=180)
        if r.ok:
            procedencia.append([url, r.url, dt.datetime.now().isoformat(timespec="seconds")])
            return pd.DataFrame(r.json())
        time.sleep(5 * (intento + 1))
    r.raise_for_status()

codigos = ",".join(f"'{c}'" for c in AEROPUERTOS)
# Pasajeros por mes y tipo de tráfico: embarcados (origen) + desembarcados (destino)
partes = []
for campo in ("origen", "destino"):
    d = soda(OD, {"$select": f"{campo} as ap,a_o,n_mero_de_mes,tr_fico_n_i,sum(pasajeros)",
                  "$where": f"{campo} in ({codigos})", "$group": f"{campo},a_o,n_mero_de_mes,tr_fico_n_i", "$limit": 50000})
    partes.append(d)
P = pd.concat(partes)
P["pax"] = pd.to_numeric(P["sum_pasajeros"])
P["fecha"] = pd.to_datetime(dict(year=P.a_o.astype(int), month=P.n_mero_de_mes.astype(int), day=1))
P["trafico"] = P["tr_fico_n_i"].map({"N": "nacional", "I": "internacional"})
M = P.pivot_table(index=["ap", "fecha"], columns="trafico", values="pax", aggfunc="sum", fill_value=0).reset_index()
M["total"] = M["nacional"] + M["internacional"]

# Sillas ofrecidas
partes = []
for campo in ("origen", "destino"):
    partes.append(soda(EQ, {"$select": f"{campo} as ap,a_o,mes,sum(sillas_ofrecidas)", "$where": f"{campo} in ({codigos})",
                            "$group": f"{campo},a_o,mes", "$limit": 50000}))
S = pd.concat(partes)
S["sillas"] = pd.to_numeric(S["sum_sillas_ofrecidas"])
S["fecha"] = pd.to_datetime(dict(year=S.a_o.astype(int), month=S.mes.astype(int), day=1))
S = S.groupby(["ap", "fecha"], as_index=False)["sillas"].sum()
M = M.merge(S, on=["ap", "fecha"], how="left")
M.to_csv(DATOS / "conectividad_mensual.csv", index=False)

# Destinos con vuelo directo (regular o no) por año, desde cada aeropuerto
R = soda(OD, {"$select": "origen,a_o,destino,pais_destino,sum(pasajeros)", "$where": f"origen in ({codigos}) and pasajeros > 0",
              "$group": "origen,a_o,destino,pais_destino", "$limit": 50000})
R["pax"] = pd.to_numeric(R["sum_pasajeros"]); R["anio"] = R.a_o.astype(int)
R = R[R.pax >= 1000]     # rutas con al menos 1.000 pasajeros al año (descarta vuelos esporádicos)
R["internacional"] = R["pais_destino"].str.upper() != "COLOMBIA"
rutas = R.groupby(["origen", "anio"]).agg(destinos=("destino", "nunique"),
                                          destinos_internacionales=("internacional", "sum")).reset_index()
paises = R[R.internacional].groupby(["origen", "anio"])["pais_destino"].apply(lambda s: sorted(set(s))).reset_index()
rutas = rutas.merge(paises, on=["origen", "anio"], how="left")
rutas.to_csv(DATOS / "conectividad_rutas.csv", index=False)

# Resumen anual
A = M.assign(anio=M.fecha.dt.year).groupby(["ap", "anio"]).agg(total=("total", "sum"), internacional=("internacional", "sum"),
                                                                sillas=("sillas", "sum"), meses=("fecha", "count")).reset_index()
A["grupo"] = A.ap.map(lambda a: AEROPUERTOS[a][1]); A["ciudad"] = A.ap.map(lambda a: AEROPUERTOS[a][0])
A["pct_internacional"] = A.internacional / A.total
G = A[A.meses == 12].groupby(["grupo", "anio"])[["total", "internacional"]].sum().reset_index()
print(A[A.anio.between(2022, 2025)].pivot(index="ap", columns="anio", values="total").div(1e6).round(2))
print(G.pivot(index="grupo", columns="anio", values="total").div(1e6).round(2))
print(rutas[rutas.anio == 2025][["origen", "destinos", "destinos_internacionales", "pais_destino"]].to_string(index=False))

ultimo = M.fecha.max()
web = {
    "generado": dt.datetime.now().strftime("%Y-%m-%d %H:%M"), "ultimo": f"{ultimo:%Y-%m}",
    "aeropuertos": {k: {"ciudad": v[0], "grupo": v[1]} for k, v in AEROPUERTOS.items()},
    "anual": A.round(4).to_dict(orient="records"),
    "grupos": G.to_dict(orient="records"),
    "mensual": {a: {"fechas": M[M.ap == a].sort_values("fecha").fecha.dt.strftime("%Y-%m").tolist(),
                    "total": M[M.ap == a].sort_values("fecha").total.tolist(),
                    "internacional": M[M.ap == a].sort_values("fecha").internacional.tolist()} for a in AEROPUERTOS if a in set(M.ap)},
    "rutas": [{**r, "pais_destino": r["pais_destino"] if isinstance(r["pais_destino"], list) else []}
              for r in rutas.to_dict(orient="records")],
    "documento": {"corredor": [3.3, 4.3], "puertas": [9.4, 10.3], "barranquilla": [3.0, 3.2], "cartagena": [6.4, 7.1]},
    "sha256": hashlib.sha256((DATOS / "conectividad_mensual.csv").read_bytes()).hexdigest(),
}
(RES / "web_contexto.json").write_text(json.dumps(web, ensure_ascii=False, default=int), encoding="utf-8")
pd.DataFrame(procedencia, columns=["servicio", "consulta", "fecha"]).to_csv(DATOS / "_procedencia_conectividad.csv", index=False)

DescargarSHA-256 b8331659d763c20c022c068ab41c9ce3d6883ea002c77b709586c7733da45fd7

# ATARRAYA AZUL · pruebas de pronóstico, alerta y conectividad
# Python 3.12
pandas==2.2.3
numpy==2.1.2
statsmodels==0.14.6
xgboost==3.3.0
scipy==1.18.0
matplotlib==3.10.9
openpyxl==3.1.5
requests==2.32.3
pillow==11.1.0