# Importación de librería y funciones específicas:
import os
import warnings
import pickle
import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
import statsmodels.api as sm
import statsmodels as sm2
from statsmodels.tsa.api import VAR
from statsmodels.formula.api import ols
from numpy import transpose
from numpy.linalg import inv, slogdet
from scipy.stats import skew, kurtosis
from scipy.optimize import minimize
from matplotlib.ticker import FormatStrFormatter
from matplotlib.colors import to_rgba
from matplotlib.ticker import MultipleLocator
from openpyxl.styles import PatternFill, Font
from openpyxl.utils import get_column_letter

# Parámetro condicional:
optimizar = 1
archivo_resultado = 'result_var.pkl'

# Parámetros de colores (código Hex):
color_azul_principal  = '#1F3F6D'
color_azul_secundario = '#2A5593'
color_rojo_principal  = '#920000'
color_rojo_secundario = '#92172F'
color_gris_principal  = '#D9D9D9'
color_gris_secundario = '#8AABB8'

# Función para crear vector que escala variables según efecto de la pandemia
def crear_escalar(data, s , rho): 
    #
    # Creamos variable escalar (s)
    escalar = np.ones([len(data), 1]) 
    #
    # Buscamos cuándo empieza la pandemia
    inicio_covid = data.index.get_loc('2020-06').start
    #
    # Modificamos escalar tras inicio de covid
    escalar[inicio_covid:inicio_covid + 3] = s 
    #
    # Efecto decae con el tiempo
    distancia = np.arange(len(escalar)) - (inicio_covid + 2)
    #
    escalar[inicio_covid + 3:] = np.array(1 + (s[2, 0] - 1)*rho**np.maximum(distancia[inicio_covid + 3:], 0)).reshape(-1, 1)
    #
    return escalar

# Función para estimar VAR con restricciones de LP y control por volatilidad de la pandemia
def negative_ll_restricted_VAR(params): 
    #
    # Parámetros:
    gamma = params[:len(gamma_0)].reshape(-1, 1)
    theta = params[len(gamma_0):]
    s     = theta[:-1].reshape(-1, 1)
    rho   = theta[-1]
    #
    # Creamos vector para reescalar shocks de pandemia
    escalar = crear_escalar(data_var, s, rho) 
    #
    # Calculamos variables reescaladas
    data_var_s = data_var/escalar          
    #
    # Estimamos errores
    #   Preparación datos
    #     Variables endógenas
    X_s = np.array(data_var_s[1:]) 
    #
    #     Vectorización:
    x_s = X_s.flatten().reshape(-1, 1) 
    #
    #     Variables endógenas rezagadas:
    Z_s = np.concatenate((np.ones([len(data_var_s), 1]), data_var_s), axis = 1)[:-1].T 
    #
    # Residuos:
    Z_s_kron = np.kron(Z_s.T,I_L)
    z_s = x_s - Z_s_kron@r
    # Estimador MCO:
    betas = inv(Z_s@Z_s.T)@(Z_s@X_s)
    #
    # Estimación errores
    u = z_s - Z_s_kron@(R@gamma)
    U = X_s - Z_s.T@betas
    # Varianza-covarianza de los residuos:
    sigma2 = (U.T@U)/(len(data_var) - 1)
    #
    # Función LL 
    nll = -n_v*np.sum(np.log(escalar[1:])) - 0.5*len(data_var)*slogdet(sigma2)[1] - (u.T@np.kron(I_T, inv(sigma2))@u)[0,0]
    #
    # Retorna el negativo del log-likelihood (para maximizar en vez de minimizar)
    return -nll

# Ignorar advertencias de openpyxl relacionadas con formato condicional
warnings.filterwarnings("ignore", category = UserWarning, module = "openpyxl")

# Nombre del fichero Excel con las bases de datos:
file = 'Datos/base_datos.xlsx'

# Importación de base con datos trimestrales:
data_trim = pd.read_excel(file, sheet_name = 'trimestral', header = 0, index_col = 'fecha')

# Asegurarnos que el índice es frecuencia trimestral:
data_trim.index = pd.to_datetime(data_trim.index)
data_trim = data_trim.asfreq('QE')  # Año calendario

# Mostrar las últimas filas con formato
data_trim.tail().style\
    .format(precision = 1, thousands = ".", decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Importación de base de datos anuales:
data_anual = pd.read_excel(file, sheet_name = 'anual', header = 0, index_col = 'fecha')

# Asegurarnos que el índice es frecuencia anual:
data_anual.index = pd.to_datetime(data_anual.index)
data_anual = data_anual.asfreq('YE')  # Año calendario

# Mostrar las últimas filas con formato
data_anual.tail().style\
    .format(precision = 1, thousands = ".", decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Importación de base con proyecciones del MDD:
deuda_det = pd.read_excel(file, sheet_name = 'Deuda modelo deterministico', header = 0, index_col = 'fecha')*100

# Asegurarnos que el índice es frecuencia anual:
deuda_det.index = pd.to_datetime(deuda_det.index)
deuda_det = deuda_det.asfreq('YE')  # Año calendario

# Mostrar las últimas filas con formato
deuda_det.tail().style\
    .format(precision = 1, decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Seed para números aleatorios
np.random.seed(123)

# Fecha de calibración: 04 de junio de 2025.

# Parámetros
trim_final = '2025-03' # Último mes de último trimestre usado en cálculos trimestrales (formato yyyy-mm). /data_trim.iloc[:,0].astype(str).values[-1][:7]
y_final    = '2024-12-31' # Último mes de último año usado en cálculos anuales. / data_anual[data_anual['bal'].notna()]['fecha'].astype(str).values[-1][:7]

n_h        = 20        # Trimestres a proyectar (4*años a estimar).
n_b        = 20_000    # Número de simulaciones (replicaciones).

prop       = 0.362     # Proporción de deuda en moneda extranjera (a fin de último año disponible).
prop_uf    = 0.296     # Proporción de deuda en UF (a fin de último año disponible).
dur_p      = 9.9       # Duración promedio de los bonos en pesos chilenos últimos 10 años.
dur_uf     = 12.0      # Duración promedio de los bonos en UF últimos 10 años.
dur_f      = 12.7      # Duración promedio de los bonos en moneda extranjera últimos 10 años.
j_d_0      = 0.0281    # Tasa de interés efectiva de deuda en moneda local al inicio de la proyección.
j_e_0      = 0.0261    # Tasa de interés efectiva de deuda en moneda extranjera al inicio de la proyección.

p_ref_t    = 4.09      # Precio de referencia del cobre del año en curso.
infl_usa0  = 0.030     # Proyección de inflación de EE.UU para el año en curso.
infl_usa1  = 0.025     # Proyección de inflación de EE.UU para el año siguiente.
infl_usa2  = 0.021     # Proyección de inflación de EE.UU para dos años más.
infl_usa3  = 0.022     # Proyección de inflación de EE.UU para tres años más.
infl_usaLR = 0.020     # Proyección de inflación de EE.UU para cuatro años más.

bpmf       = 0.04      # Balance Primario Máximo Factible.
riesgo     = 0.1       # Proporción máxima de veces que se puede superar el bpmf.

# Valores de largo plazo (LP) para las variables en el VAR
g_ss       = ((0.02 + 1)**(1/4) - 1)*100                                               # Crecimiento tendencial del PIB
def_ss     = ((0.03 + 1)**(1/4) - 1)*100                                               # Deflactor del PIB en el LP
e_ss       = ((0.01 + 1)**(1/4) - 1)*100                                               # Depreciación del TCN trimestral en el LP
c_ss       = 4.09/(1 + infl_usa0) if float(trim_final[5:7]) in range(9, 12) else 4.09  # precio del cobre en el LP   /(1+infl_usa0)
i_ss       = 5.6 - 3.0                                                                 # Tasa de interés local real en el LP
iu_ss      = 3.35 - 2.0 + 1.25                                                         # Tasa de interés extranjera real en el LP (tasa USA a LP - meta inf. + spread)
inf_ss     = ((0.03 + 1)**(1/4) - 1)*100                                               # Inflación trimestral según IPC en el LP

# Opción de añadir directamente estimación de balance primario para el año en curso.
usar_bp_0  = 0       # 1 para usar pb_0, 0 para no usarlo
bp_0       = -0.005

# Opción de añadir Precio de referencia del cobre del próximo año (de haber).
usar_pr_1  = 0       # 1 para usar p_ref_1, 0 para no usarlo
p_ref_1    = 4.09    # Precio de referencia a usar si usar_pr_1 = 1

# Percentiles para fancharts
p1         = 5
p2         = 20
p3         = 30
p4         = 40

# Añadir variables en variación trimestral porcentual
data2 = data_trim.pct_change(fill_method = None)*100
data2.drop(columns = ['btu10', 'treasury_r10'], errors = 'ignore', inplace = True)
data2 = data2.add_suffix('_chg')
data2.rename(columns = {'gdpr_cl_chg': 'g_real',
                        'gdpr_usa_chg': 'g_real_usa',
                        'ipc_cl_chg': 'inf_cl',
                        'tc_clp_usd_chg': 'dep_tcn',
                        'cpi_usa_chg': 'inf_usa'}, inplace = True)

# Añadir estimación de tasa de interés real de deuda en moneda extranjera de Chile
data_trim['i_usa10_mas_EMBI'] = data_trim['treasury_r10'] + data_trim['embi_cl']/100

# Añadir precio real del cobre
ipc_usa_base = data_trim.loc['2023-12-31', 'cpi_usa']
data_trim['copper_r'] = data_trim['copper_nom']*ipc_usa_base/data_trim['cpi_usa']

# Unir bases y seleccionar periodo para preparación de base de datos
data_trim = pd.concat([data_trim, data2], axis = 1).loc['2003-03':trim_final]

# Identificar variables duplicadas
duplicadas = data_trim.columns[data_trim.columns.duplicated()]
if not duplicadas.empty:
    print("Variables duplicadas:", list(duplicadas))

# Mostrar las últimas filas con formato
data_trim.tail().style\
    .format(precision = 1, thousands = ".", decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Transformar y crear variables anuales

# Variables básicas:
data_anual['g']           = data_anual['gdpr'].pct_change(fill_method = None)          # Crecimiento real
data_anual['inf_ipc']     = data_anual['ipc'].pct_change(fill_method = None)*100       # Inflación IPC
data_anual['chg_c_n']     = data_anual['copper_n'].diff()                              # Var. precio cobre nominal
data_anual['chg_p_ref']   = data_anual['copper_ref'].diff()                            # Var. precio de referencia

# PIB nominal en millones
gdpn_m                    = data_anual['gdpn']*1_000                                   # PIB nominal

# Cálculo de razones fiscales sobre PIB
for var in ['bp', 'bal', 'gint', 'orc', 'deuda']:
    data_anual[f'{var}_pib'] = data_anual[var]/gdpn_m

data_anual.rename(columns = {'bp_pib': 'bp_pib',                                       # Balance primario (% del PIB)
                             'bal_pib': 'b_pib',                                       # Balance total (% del PIB)
                             'gint_pib': 'gi_pib',                                     # Gastos por intereses (% del PIB)
                             'orc_pib': 'orc_pib',                                     # ORC (% del PIB)
                             'deuda_pib': 'd_p'},                                      # Deuda (% del PIB)
                  inplace = True)    

# Rezagos y potencias de la deuda
data_anual['d_p_lag']         = data_anual['d_p'].shift()                              # Rezago de la deuda (% del PIB)
data_anual['d_p_lag_squared'] = data_anual['d_p_lag']**2                               # Rezago de la deuda al cuadrado
data_anual['d_p_lag_cubed']   = data_anual['d_p_lag']**3                               # Rezago de la deuda al cubo

# Fondos del Tesoro Público como porcentaje del PIB
for fondo in ['fees',                                                                  # FEES
              'frp',                                                                   # FRP
              'oatp',                                                                  # OATP
              'oatp2',                                                                 # OATP sin Fondo por la Educación
              'ret_fees']:                                                             # Retiros del FEES
    data_anual[f'{fondo}_p'] = data_anual[fondo]*data_anual['tcn_a']/gdpn_m

# Deuda neta = deuda bruta - activos del Tesoro Público
data_anual['deuda_neta']     = data_anual['d_p'] - (data_anual['fees_p'] + data_anual['frp_p'] + data_anual['oatp_p'])

# Variables dummy
data_anual['crisis']         = (data_anual['g'] < 0).astype(int)
data_anual['regla_fiscal']   = (data_anual.index >= '2001-12-01').astype(int)

# Brecha del cobre (rellenando NaN por 0)
data_anual['brecha_cobre'] = (data_anual['copper_n'] - data_anual['copper_ref']).fillna(0)

# Verificar variables duplicadas
duplicadas = data_anual.columns[data_anual.columns.duplicated()]
if not duplicadas.empty:
    print("Variables duplicadas:", list(duplicadas))

# Mostrar las últimas filas con formato
data_anual.tail().style\
    .format(precision = 1, thousands = ".", decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Calcular brecha de producto con filtro HP sobre el PIB real desde 1980
gap_aux = data_anual.loc['1980-12':y_final, ['gdpr']].copy()
gap_aux['cycle'], gap_aux['trend'] = sm.tsa.filters.hpfilter(gap_aux['gdpr'], lamb=6.25)

# Brecha del PIB: desviación porcentual del PIB respecto de su tendencia
data_anual['gdpr_gap'] = gap_aux['gdpr'] / gap_aux['trend'] - 1

# Brecha positiva y negativa como variables separadas
data_anual['brecha_pos']      = (data_anual['gdpr_gap'] > 0).astype(int)
data_anual['brecha_pos_int']  = data_anual['gdpr_gap'].where(data_anual['gdpr_gap'] > 0, 0)
data_anual['brecha_neg_int']  = data_anual['gdpr_gap'].where(data_anual['gdpr_gap'] <= 0, 0)

# Limitar periodo de análisis
data_anual = data_anual.loc['1990-12':y_final]

# Detectar y reportar variables duplicadas
duplicadas = data_anual.columns[data_anual.columns.duplicated()]
if not duplicadas.empty:
    print("Variables duplicadas:", list(duplicadas))

# Mostrar las últimas filas con formato
data_anual.tail().style\
    .format(precision = 1, thousands = ".", decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Selección de variables para graficar:
data_plot = data_trim[['g_real','def_chile_chg','dep_tcn','copper_r','btu10','i_usa10_mas_EMBI','inf_cl']].copy()
data_plot.rename(columns = {'g_real': 'Crecimiento real',
                            'def_chile_chg': 'Variación deflactor del PIB',
                            'dep_tcn': 'Depreciación TCN',
                            'copper_r': 'Precio real del cobre',
                            'btu10': 'BTU a 10 años',
                            'i_usa10_mas_EMBI': 'US treasury a 10 años indexado + EMBI',
                            'inf_cl': 'Variación del IPC'}, inplace = True)

# Número de variables a graficar
n_vars = len(data_plot.columns)

# Gráficos:
fig, axes = plt.subplots(nrows = 4, ncols = 2, figsize = (8, 10), dpi = 180)
axes = axes.flatten()

for i, ax in enumerate(axes[:len(data_plot.columns)]):
    # Selección de variable:
    data_aux = data_plot[data_plot.columns[i]]

    # Formateo de estilo:
    ax.plot(data_aux, color = color_azul_principal, linewidth = 1.5, label = data_plot.columns[i])
    ax.set_title(data_plot.columns[i], fontsize = 12, fontweight = 'bold', fontfamily = 'Tw Cen MT')
    ax.spines[['top', 'right']].set_visible(False)
    ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
    ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
    ax.axhline(y = 0, color = 'black', alpha = 0.3, linestyle = '-', linewidth = 0.5)
    ax.set_ylabel('US$ de 2023' if data_plot.columns[i] == 'Precio real del cobre' else 'Porcentaje', fontsize = 12, fontfamily = 'Tw Cen MT')
        
# Si hay menos variables que subplots, apagar los ejes sobrantes
for j in range(n_vars, len(axes)):
    axes[j].axis('off')

plt.tight_layout(pad = 2.0)
plt.subplots_adjust(top = 0.92)
plt.suptitle("Variables económicas usadas en el modelo VAR para MED", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
plt.show()

# Selección y limpieza de datos para el VAR
data_var = (data_trim[['g_real', 'def_chile_chg', 'dep_tcn', 'copper_r', 'btu10', 'i_usa10_mas_EMBI', 'inf_cl']].copy().dropna().loc[lambda df: ~df.index.duplicated()].sort_index())
data_var.index = pd.to_datetime(data_var.index)
data_var = data_var.asfreq('QE') 

# Estimación del modelo VAR con 1 rezago y constante
var_model = VAR(data_var)
var_model_fitted = var_model.fit(maxlags = 1, trend = 'c', verbose = True)

# Visualización de los coeficientes estimados
var_model_fitted.params

# Diccionario de nombres descriptivos
nombres = {
    'g_real': 'Crecimiento real',
    'def_chile_chg': 'Variación deflactor del PIB',
    'dep_tcn': 'Depreciación TCN',
    'copper_r': 'Precio real del cobre',
    'btu10': 'BTU a 10 años',
    'i_usa10_mas_EMBI': 'US treasury a 10 años indexado + EMBI',
    'inf_cl': 'Variación del IPC'}

# Obtener residuos
residuos = var_model_fitted.resid

# Crear gráfico
fig, ax = plt.subplots(dpi=180)
color_base = np.array(to_rgba(color_azul_secundario))
n_vars = residuos.shape[1]

# Graficar residuos con gradiente de color
for i, col in enumerate(residuos.columns):
    color = color_base.copy()
    factor = 1 - i/(1.5*n_vars)
    color[:2] = np.minimum(color[:2] + 0.2*i/n_vars, 1.0) 
    color[3] = factor  # opacidad

    ax.plot(residuos.index, residuos[col], linewidth = 1.5, label = nombres.get(col, col), color = color)

# Formato gráfico
ax.set_title("Residuos del modelo VAR", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
ax.set_yticks(ax.get_yticks())
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out')
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)
ax.legend(loc = 'upper right', fontsize = 8, frameon = False)
plt.tight_layout()
plt.show()

# Construcción de matriz R de tamaño L(L + 1)xV = 56x49
R = np.zeros((56, 49))

# Valores negativos del estado estacionario
valores_ss = np.array([-g_ss, -def_ss, -e_ss, -c_ss, -i_ss, -iu_ss, -inf_ss])

# Restricciones de estado estacionario en las primeras 7 filas
for i in range(7):
    R[i, i::7] = valores_ss

# Coeficientes sin restricciones (identidad desplazada hacia abajo)
for i in range(7, 56):
    R[i, i - 7] = 1

# Print:
R

array([[-0.49629316,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        , -0.49629316,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        , -0.49629316, ...,  0.        ,
         0.        ,  0.        ],
       ...,
       [ 0.        ,  0.        ,  0.        , ...,  1.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         1.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        ,  1.        ]])

# Construcción vector r de tamaño L(L + 1) = 56x1:
r = np.concatenate([-valores_ss, np.zeros(49)]).reshape(56, 1)

# Print:
r

array([[0.49629316],
       [0.74170718],
       [0.24906793],
       [4.09      ],
       [2.6       ],
       [2.6       ],
       [0.74170718],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ],
       [0.        ]])

# Número de variables endógenas y observaciones
n_v, TT = data_var.shape[1], len(data_var)

# Coeficientes iniciales (sin constante) vectorizados
gamma_0 = var_model_fitted.params.iloc[1:].to_numpy().flatten()

# Matriz de variables endógenas (sin primera fila por rezago implícito)
X = data_var.iloc[1:].to_numpy()
x = X.flatten().reshape(-1, 1)

# Matriz de regresores Z (constante + variables rezagadas)
Z = np.column_stack([np.ones(TT), data_var.to_numpy()])[:-1].T  # Shape: (n_v + 1, T - 1)

# Matrices identidad
I_L = np.eye(n_v)
I_T = np.eye(TT - 1)

# Vector r de largo L(L+1) ya definido antes
# Cálculo de residuos del modelo lineal bajo estado estacionario
z = x - np.kron(Z.T, I_L) @ r

# Parámetros iniciales de volatilidad (relacionados con pandemia)
theta = np.array([4.0, 3.0, 2.0, 0.5])  # [s0, s1, s2, rho]

# Concatenar todo en el vector inicial de parámetros
initial_guess = np.concatenate([gamma_0, theta])

# Restricciones (bounds) para optimización
eps = 1e-8
bounds = [(None, None)] * len(gamma_0) + [(eps, None)] * 3 + [(0, 0.7)]

## Realizamos maximización
if optimizar == 1:
    print("Estimando modelo VAR con restricciones de Largo Plazo...")
    result = minimize(negative_ll_restricted_VAR,
                      initial_guess,
                      method = 'L-BFGS-B',
                      bounds = bounds)
    # Guardando los resultados:
    with open(archivo_resultado, 'wb') as f:
        pickle.dump(result, f)
    
    print(f"Optimización completada y guardada en el archivo: {archivo_resultado}")
    
elif optimizar == 2:
    if os.path.exists(archivo_resultado):
        with open(archivo_resultado, 'rb') as f:
            result = pickle.load(f)
        print(f"Resultado cargado desde archivo: {archivo_resultado}")
    else:
        raise FileNotFoundError(f"No se encontró el archivo {archivo_resultado}. Correr con optimizar = 1 primero.")
        
else:
    raise ValueError("La variable 'optimizar' debe ser 1 (para estimar) o 2 (para cargar resultado).")

Estimando modelo VAR con restricciones de Largo Plazo...
Optimización completada y guardada en el archivo: result_var.pkl

# Separar parámetros estimados
gamma_tilda    = result.x[:len(gamma_0)]
theta          = result.x[len(gamma_0):]
s_est, rho_est = theta[:-1].reshape(-1, 1), theta[-1]

# Coeficientes estructurales estimados del VAR con restricciones
beta_tilda = (R @ gamma_tilda.reshape(-1, 1)) + r
B_tilda    = beta_tilda.reshape(n_v + 1, n_v)

# Escalamiento por efectos de pandemia
escalar = crear_escalar(data_var, s_est, rho_est)

# Ajustar base de datos por escalamiento
data_var_s = data_var / escalar

# Matriz de variables endógenas y sus rezagos
X_s = data_var_s.iloc[1:].values
Z_s = np.column_stack([np.ones(len(data_var_s)), data_var_s.values])[:-1]

# Cálculo de residuos
U = X_s - Z_s @ B_tilda

# Estimación de matriz de varianzas y covarianzas
sigma_tilda = U.T @ U / (len(data_var_s) - 1)

# Separar constante y coeficientes autoregresivos del VAR
constantes = B_tilda[0, :]       # Vector de constantes
alphas     = B_tilda[1:, :].T    # Coeficientes por ecuación (forma más natural para simulaciones)

# Alias para matriz de varianzas
sigma2 = sigma_tilda

# Obtener residuos
residuos = var_model_fitted.resid

# Crear gráfico
fig, ax = plt.subplots(dpi=180)
color_base = np.array(to_rgba(color_azul_secundario))
n_vars = residuos.shape[1]

# Graficar residuos con gradiente de color
for i in range(n_vars):
    color = color_base.copy()
    factor = 1 - i/(1.5*n_vars)
    color[:2] = np.minimum(color[:2] + 0.2*i/n_vars, 1.0) 
    color[3] = factor  # opacidad

    ax.plot(residuos.index, U[:,i], linewidth = 1.5, label = nombres.get(residuos.columns[i], residuos.columns[i]), color = color)

# Formato gráfico
ax.set_title("Residuos del modelo VAR restringido y escalado en periodo de pandemia Covid-19", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
ax.set_yticks(ax.get_yticks())
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out')
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)
ax.legend(loc = 'upper right', fontsize = 8, frameon = False)
plt.tight_layout()
plt.show()

# Proyección de shocks con distribución normal y trayectorias VAR

# Seed para números aleatorios
np.random.seed(123)

# Simular shocks VAR: (n_h, n_v, n_b)
shocks_var = np.random.multivariate_normal(mean = np.zeros(n_v), cov = sigma2, size = (n_b, n_h))
shocks_var = np.transpose(shocks_var, (1, 2, 0))  # Reordenar a (n_h, n_v, n_b)

# Inicializar trayectorias: (n_h + 1, n_v, n_b)
tray = np.empty((n_h + 1, n_v, n_b))

# Última observación como punto de partida
tray[0, :, :] = data_var.iloc[-1, :].to_numpy()[:, None]

# Generar proyecciones iterativas
for j in range(1, n_h + 1):
    tray[j, :, :] = constantes[:, None] + alphas@tray[j - 1, :, :] + shocks_var[j - 1, :, :]

# Remover primer periodo (observado), nos quedamos con proyecciones
tray = tray[1:, :, :]

# Gráficos de fancharts de simulaciones de variables económicas
# títulos para cada fanchart
titles = {'g_real': 'Crecimiento real',
          'def_chile_chg': 'Variación deflactor del PIB',
          'dep_tcn': 'Depreciación TCN',
          'copper_r': 'Precio real del cobre',
          'btu10': 'BTU a 10 años',
          'i_usa10_mas_EMBI': 'US Treasury a 10 años indexado + EMBI',
          'inf_cl': 'Variación del IPC'}

timestamps = data_var.index[-40:]
trim_sim = pd.date_range(start = timestamps[-1], periods = n_h + 1, freq = 'QE')

fig, axes = plt.subplots(nrows = 4, ncols = 2, figsize = (10, 12), dpi = 180)
axes = axes.flatten()

for i, var in enumerate(data_var.columns):
    ax = axes[i]
    perc = np.percentile(tray[:, i, :], np.arange(10, 100, 10), axis = 1)
    ultimo_efectivo = data_var.iloc[-1, i]
    mediana_simulacion = np.concatenate([[ultimo_efectivo], perc[4]])

    ax.plot(timestamps, data_var.iloc[-40:, i], color = 'black', linewidth = 1.5, label = "Datos históricos")
    ax.plot(trim_sim, mediana_simulacion, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")
    
    for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
        lower = np.concatenate([[ultimo_efectivo], perc[j]])
        upper = np.concatenate([[ultimo_efectivo], perc[8 - j]])
        ax.fill_between(trim_sim, lower, upper, color = color_rojo_secundario, alpha = alpha, label = f"Rango P{10*(j + 1)} - P{90 - 10*j}")

    # Línea punteada vertical y texto
    ax.axvline(x = trim_sim[0], color = 'gray', linestyle = '--', linewidth = 1)
    ax.text(trim_sim[0], ax.get_ylim()[1]*0.95, "Proyecciones", rotation = 90, fontsize = 9,
            verticalalignment = 'top', color = 'gray', fontfamily = 'Tw Cen MT')
    
    ax.set_title(titles[var], fontsize = 12, fontweight = 'bold', fontfamily = 'Tw Cen MT')
    ax.set_xlabel("Trimestres", fontsize = 10, fontfamily = 'Tw Cen MT')
    ax.set_ylabel('US$ de 2023' if data_plot.columns[i] == 'Precio real del cobre' else 'Porcentaje', fontsize = 12, fontfamily = 'Tw Cen MT')
    ax.spines[['top', 'right']].set_visible(False)
    ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
    ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
    ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)

        
# Si hay menos variables que subplots, apagar los ejes sobrantes
for j in range(len(titles), len(axes)):
    axes[j].axis('off')

# Definir la leyenda en el espacio vacío del subplot:
handles, labels = axes[0].get_legend_handles_labels()
axes[-1].legend(handles, labels, loc = 'center left', fontsize = 9, frameon = False)

plt.tight_layout()
plt.subplots_adjust(top = 0.95)
plt.suptitle("Fan Charts de Variables Económicas", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
plt.show()

# Añadimos a las trayectorias las observaciones de la base de datos que son del año en curso.
# Obtener número de trimestres del año en curso ya observados
trimestres_actual = int(data_trim.index[-1].month / 3)
obs_actuales = data_var.tail(trimestres_actual).to_numpy() if trimestres_actual <= 3 else None

# Inicializar matriz extendida
tray_aux = np.empty((n_h + trimestres_actual, n_v, n_b))

for i in range(n_b):
    if trimestres_actual <= 3 and obs_actuales is not None:
        tray_aux[:, :, i] = np.vstack([obs_actuales, tray[:, :, i]])
    else:
        tray_aux[:n_h, :, i] = tray[:, :, i]

# Cálculo de trayectorias anuales (tasa de crecimiento o promedio)
n_años = int(n_h / 4)
tray_anual = np.empty((n_años, n_v, n_b))

# Índices de diciembre de cada año proyectado
diciembre_idx = np.arange(3, 3 + 4*n_años, 4)

# Cálculo vectorizado por variable, muestra y año
for i in range(n_v):
    if i < 3 or i == 6:  # crecimiento compuesto
        for k in range(n_b):
            for t, j in enumerate(diciembre_idx):
                vals = tray_aux[j - 3:j + 1, i, k]/100 + 1
                tray_anual[t, i, k] = 100*(np.prod(vals) - 1)
    else:  # promedio simple
        for k in range(n_b):
            for t, j in enumerate(diciembre_idx):
                tray_anual[t, i, k] = tray_aux[j - 3:j + 1, i, k].mean()

# Gráficos de fancharts de simulaciones de variables macroeconómicas en años
año_inicio = data_var.index.year[-1:] + 1 if data_var.index.month[-1:] == 12 else data_var.index.year[-1:]
años_sim = np.arange(año_inicio[0], año_inicio[0] + int(n_h/4)).astype(str)

# Diccionario de títulos
titles = {'g_real': 'Crecimiento real',
          'def_chile_chg': 'Variación deflactor del PIB',
          'dep_tcn': 'Depreciación TCN',
          'copper_r': 'Precio real del cobre',
          'btu10': 'BTU a 10 años',
          'i_usa10_mas_EMBI': 'US Treasury a 10 años indexado + EMBI',
          'inf_cl': 'Variación del IPC'}

fig, axes = plt.subplots(nrows = 4, ncols = 2, figsize = (10, 12), dpi = 180)
axes = axes.flatten()

for i, var in enumerate(data_var.columns):
    ax = axes[i]
    
    # Extraer percentiles desde trayectorias anuales
    perc = np.percentile(tray_anual[:, i, :], np.arange(10, 100, 10), axis = 1)
    mediana_simulacion = perc[4]  # No se concatena nada extra

    # Plot principal
    ax.plot(años_sim, mediana_simulacion, color = color_rojo_principal, linestyle = '--',
            linewidth = 2, label = "Mediana simulación")
    
    # Rellenar bandas de confianza
    for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
        ax.fill_between(años_sim, perc[j], perc[8 - j], color = color_rojo_secundario, alpha = alpha, label = f"Rango P{10*(j + 1)} - P{90 - 10*j}")

    # Estética 
    ax.set_title(titles[var], fontsize = 12, fontweight = 'bold', fontfamily = 'Tw Cen MT')
    ax.set_xlabel("Años", fontsize = 10, fontfamily = 'Tw Cen MT')
    ax.set_ylabel('US$ de 2023' if var == 'copper_r' else 'Porcentaje', fontsize = 12, fontfamily = 'Tw Cen MT')
    ax.spines[['top', 'right']].set_visible(False)
    ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
    ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
    if (perc < 0).any():
        ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)

# Si hay menos variables que subplots, apagar los ejes sobrantes
for j in range(len(titles), len(axes)):
    axes[j].axis('off')

# Definir la leyenda en el espacio vacío del subplot:
handles, labels = axes[0].get_legend_handles_labels()
axes[-1].legend(handles, labels, loc = 'center left', fontsize = 9, frameon = False)
plt.tight_layout()
plt.subplots_adjust(top = 0.95)
plt.suptitle("Fan Charts de Variables Económicas (anual)", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
plt.show()

# Gráfico de variables usadas en FRF
# Variables a graficar y nombres descriptivos
variables = ['bp_pib', 'd_p', 'brecha_cobre', 'gdpr_gap', 
             'crisis', 'regla_fiscal', 'brecha_pos', 'ret_fees_p']

nombres = ['Balance primnario',
           'Deuda bruta del GC',
           'Diferencia precio del cobre y su referencial',
           'Brecha del producto',
           'Dummy de crisis económica',
           'Dummy Regla fiscal en uso',
           'Brecha del producto positiva',
           'Retiros del FEES (% del PIB)']

y_labels = ['% del PIB', '% del PIB', 'USc/lb', '% desviación', 
            '1 = Crisis económica', '1 = Regla fiscal en uso', 
            '1 = Brecha positiva', '% del PIB']

# Preparar datos
data_plot = data_anual[variables].copy()
data_plot[['bp_pib', 'd_p']]*= 100  # Escalar a porcentaje del PIB
data_plot.columns = nombres

# Graficar
fig, axes = plt.subplots(nrows = 4, ncols = 2, figsize = (10, 12), dpi = 180)
axes = axes.flatten()

for i, ax in enumerate(axes[:len(nombres)]):
    col = nombres[i]
    ax.plot(data_plot[col], color = color_azul_principal, linewidth = 1.5, label = col)
    ax.set_title(col, fontsize = 12, fontweight = 'bold', fontfamily = 'Tw Cen MT')
    ax.set_ylabel(y_labels[i], fontsize = 9, fontfamily = 'Tw Cen MT')
    ax.spines[['top', 'right']].set_visible(False)
    ax.tick_params(axis = 'both', direction = 'out')
    ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
    ax.tick_params(labelsize = 8)
    if (data_plot[col] < 0).any():
        ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)

plt.tight_layout(pad=2.0)
plt.subplots_adjust(top=0.92)
plt.suptitle("Variables usadas en Función de Reacción Fiscal (FRF)", 
             fontsize=14, fontweight='bold', fontfamily='Tw Cen MT')
plt.show()

# Imputar valores faltantes
data_anual['ret_fees_p'] = data_anual['ret_fees_p'].fillna(0)

# Variables explicativas
exog_vars = ['d_p_lag', 'd_p_lag_squared', 'd_p_lag_cubed',
             'brecha_cobre', 'regla_fiscal', 'ret_fees_p',
             'brecha_pos_int', 'brecha_neg_int', 'crisis']

# Estimación ARIMA(1,0,0) con términos exógenos y constante
frf = sm2.tsa.arima.model.ARIMA(endog = data_anual['bp_pib'],
                                exog  = data_anual[exog_vars],
                                order = (1, 0, 0),
                                trend = 'c').fit(cov_type = "robust", method_kwargs = {"maxiter": 500})
# Guardamos parámetros obtenidos.
frf_params = frf.params 

# Resultados
print(frf.summary())

# Diagnóstico del modelo
fig = frf.plot_diagnostics(figsize = (10, 8), lags = 10)
plt.tight_layout()
plt.show()

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                 bp_pib   No. Observations:                   35
Model:                 ARIMA(1, 0, 0)   Log Likelihood                  94.695
Date:                Tue, 01 Jul 2025   AIC                           -165.390
Time:                        10:04:20   BIC                           -146.725
Sample:                    12-31-1990   HQIC                          -158.947
                         - 12-31-2024                                         
Covariance Type:               robust                                         
===================================================================================
                      coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------
const               0.0511      0.070      0.733      0.463      -0.086       0.188
d_p_lag            -0.1912      0.887     -0.216      0.829      -1.929       1.546
d_p_lag_squared    -0.5199      3.704     -0.140      0.888      -7.780       6.740
d_p_lag_cubed       2.2320      4.792      0.466      0.641      -7.160      11.625
brecha_cobre       -0.0007      0.009     -0.073      0.942      -0.019       0.017
regla_fiscal       -0.0131      0.014     -0.951      0.342      -0.040       0.014
ret_fees_p          1.4622      0.404      3.624      0.000       0.671       2.253
brecha_pos_int      0.5362      0.396      1.356      0.175      -0.239       1.311
brecha_neg_int      0.2714      0.307      0.885      0.376      -0.329       0.872
crisis              0.0050      0.014      0.358      0.720      -0.022       0.032
ar.L1               0.7344      0.186      3.949      0.000       0.370       1.099
sigma2              0.0003   5.49e-05      4.657      0.000       0.000       0.000
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):                 0.43
Prob(Q):                              0.99   Prob(JB):                         0.81
Heteroskedasticity (H):               1.72   Skew:                             0.04
Prob(H) (two-sided):                  0.36   Kurtosis:                         2.47
===================================================================================

Warnings:
[1] Quasi-maximum likelihood covariance matrix used for robustness to some misspecifications; calculated using the observed information matrix (complex-step) described in Harvey (1989).

# Calcular diferencias
# Asegurar alineación temporal entre valores ajustados y efectivos
bp_estimado = frf.fittedvalues
bp_efectivo = data_anual['bp_pib'].iloc[0:(1 + len(bp_estimado))]
diferencia = bp_efectivo.values - bp_estimado.values

# Calcular ticks simétricos para ambos ejes
tick_interval = 0.03  # puedes ajustar esto según necesidad
tick_range = np.arange(-0.09, 0.10, tick_interval)  # Por ejemplo, -9% a 9% del PIB.

# Gráfico con eje secundario
fig, ax1 = plt.subplots(dpi = 180)

# Eje primario: líneas efectivas y estimadas
ax1.plot(bp_estimado.index, bp_estimado, color = color_rojo_principal, label = 'Estimado (ajustado)')
ax1.plot(bp_efectivo.index, bp_efectivo, color = color_azul_principal, linestyle = '--', label = 'Efectivo (observado)')
ax1.set_title("Balance primario: efectivo vs estimado", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
ax1.set_xlabel("Año", fontsize = 12, fontfamily = 'Tw Cen MT')

# Estetica eje primario:
ax1.set_ylabel("Balance primario (% del PIB)", fontsize = 12, fontfamily = 'Tw Cen MT')
ax1.set_ylim(-0.09, 0.09)
ax1.set_yticks(tick_range)
ax1.spines[['top', 'right']].set_visible(False)
ax1.tick_params(axis = 'both', direction = 'out')
ax1.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax1.tick_params(labelsize = 8)
ax1.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)

# Eje secundario: diferencia como barras
ax2 = ax1.twinx()
ax2.bar(bp_estimado.index, diferencia, width = 250, color = color_gris_secundario, alpha = 0.3, edgecolor = 'black', linewidth = 0.5, label = 'Diferencia (efectivo - estimado)')
ax2.set_ylabel("Diferencia (pp del PIB)", fontsize = 12, fontfamily = 'Tw Cen MT')
ax2.set_ylim(-0.09, 0.09)
ax2.set_yticks(tick_range)
ax2.tick_params(axis = 'y', labelsize = 8)
ax2.spines[['top', 'left']].set_visible(False)

# Combinar leyendas
handles1, labels1 = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles1 + handles2, labels1 + labels2, loc = 'lower left', fontsize=9, frameon=False)

plt.tight_layout()
plt.show()

# Proyección de shocks fiscales con distribución normal

# Seed para números aleatorios
np.random.seed(123)

# Guardar parámetros estimados del modelo FRF
psi = frf.params['ar.L1']
sigma = np.sqrt((1 - psi**2)*frf.params['sigma2'])

# Creamos espacio para almacenar los shocks
shocks_FRF = np.empty([int(n_h/4), 1, n_b]) 

# Simulación de shocks FRF (normal, media cero, desviación estándar ajustada)
for j in range(0,n_b):            
    shocks_FRF[:, 0, j] = np.random.normal(loc = 0, 
                                           scale = sigma, 
                                           size = int(n_h/4))


# Agregar shocks FRF como nueva variable a la trayectoria anual
tray_anual = np.concatenate([tray_anual, shocks_FRF], axis = 1)

# Ajustar modelo AR(1) para ORC desde 2008 en adelante
orc = sm2.tsa.arima.model.ARIMA(data_anual['orc_pib']['2008':y_final], order = (1, 0, 0), trend = 'c').fit()
orc_params = orc.params
orc_shock_sd = np.sqrt((1 - orc_params['ar.L1']**2)*orc_params['sigma2'])

# Simulación de trayectorias ORC
horizonte = int(n_h / 4)
sim_orc = np.empty((horizonte, 1, n_b))
orc_const = orc_params['const']
orc_phi = orc_params['ar.L1']
orc_0 = data_anual['orc_pib'].iloc[-1]

# Simular en bloque para mejorar eficiencia
for i in range(n_b):
    sim_orc[0, 0, i] = orc_const + orc_phi*(orc_0 - orc_const) + np.random.normal(0, orc_shock_sd)
    for j in range(1, horizonte):
        sim_orc[j, 0, i] = orc_const + orc_phi*(sim_orc[j - 1, 0, i] - orc_const) + np.random.normal(0, orc_shock_sd)

# Agregar a trayectorias anuales
tray_anual = np.concatenate((tray_anual, sim_orc), axis = 1)

# Seed para números aleatorios
np.random.seed(123)

# Preparar datos para el modelo de corrección de errores
data_oatp = data_anual.loc['2007-12':y_final, ['oatp_p']]
media_oatp2 = data_anual['oatp2_p'].iloc[-10:].mean()

# Crear variables: término de error y cambio en OATP
data_oatp['oatp_mod'] = media_oatp2 - data_oatp['oatp_p'].shift(1)
data_oatp['chg_oatp'] = data_oatp['oatp_p'].diff()

# Estimar modelo sin constante (OLS sin intercepto)
modelo_oatp = sm2.regression.linear_model.OLS(data_oatp.loc['2008-12':y_final, 'chg_oatp'],
                                              data_oatp.loc['2008-12':y_final, 'oatp_mod'],
                                              hasconst = False)
oatp_result = modelo_oatp.fit()

# Guardamos parámetros necesarios para simulación
oatp_params = oatp_result.params

# Imprimir resultados
print("Coeficiente estimado:", oatp_result.params.iloc[0])
print("MSE de residuales:", oatp_result.mse_resid)

# Parámetros para simulación
phi_hat = oatp_result.params.iloc[0]
shock_sd = np.sqrt(oatp_result.mse_resid)

# Simulación de shocks
sim_shock_oatp = np.empty([int(n_h/4), 1, n_b])
for i in range(0, n_b):
    for j in range(0, int(n_h/4)):
        sim_shock_oatp[j, 0, i] = np.random.normal(0, shock_sd)

Coeficiente estimado: 0.2771481071147569
MSE de residuales: 9.481007230586045e-05

# Preparar base de datos con cambios porcentuales
data_copper = data_anual[['chg_c_n', 'chg_p_ref']].copy()

# Agregar observación para año en curso y si se da el caso, para el año siguiente
if usar_pr_1 == 0: # Agregamos variación del precio de referencia del año en curso
    data_copper.loc[pd.Timestamp(f"{data_copper.index[-1].year + 1}-12-31")] = [np.nan, p_ref_t - data_anual['copper_ref'].iloc[-1]]
else:              # Si ya existe un precio de referencia para el próximo año también agregamos variación del próximo año
    data_copper.loc[pd.Timestamp(f"{data_copper.index[-1].year  +1}-12-31")] = [np.nan, p_ref_t - data_anual['copper_ref'].iloc[-1]]
    data_copper.loc[pd.Timestamp(f"{data_copper.index[-1].year + 2}-12-31")] = [np.nan, p_ref_1 - p_ref_t.iloc[-1]]

# Crear rezagos de chg_c_n
for lag in range(1, 5):
    data_copper[f'chg_c_n_l{lag}'] = data_copper['chg_c_n'].shift(lag)

# Estimar modelo OLS con rezagos
formula = 'chg_p_ref ~ ' + ' + '.join([f'chg_c_n_l{lag}' for lag in range(1, 5)])
p_ref_ols = ols(formula, data = data_copper)

# Guardar resultados
p_ref_ols_result = p_ref_ols.fit()

# Mostrar resumen de regresión
print(p_ref_ols_result.summary())

# Guardar parámetros estimados
pref_params = p_ref_ols_result.params

                            OLS Regression Results                            
==============================================================================
Dep. Variable:              chg_p_ref   R-squared:                       0.449
Model:                            OLS   Adj. R-squared:                  0.333
Method:                 Least Squares   F-statistic:                     3.872
Date:                Tue, 01 Jul 2025   Prob (F-statistic):             0.0183
Time:                        10:04:35   Log-Likelihood:                 9.2722
No. Observations:                  24   AIC:                            -8.544
Df Residuals:                      19   BIC:                            -2.654
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.0581      0.043      1.347      0.194      -0.032       0.148
chg_c_n_l1     0.2125      0.070      3.023      0.007       0.065       0.360
chg_c_n_l2     0.1576      0.070      2.246      0.037       0.011       0.305
chg_c_n_l3     0.0689      0.071      0.975      0.342      -0.079       0.217
chg_c_n_l4     0.1101      0.070      1.567      0.134      -0.037       0.257
==============================================================================
Omnibus:                        1.721   Durbin-Watson:                   1.654
Prob(Omnibus):                  0.423   Jarque-Bera (JB):                1.078
Skew:                           0.519   Prob(JB):                        0.583
Kurtosis:                       2.955   Cond. No.                         2.25
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

# Gráfico precio de referencia efectivo vs estimado dentro de la muestra

# Calcular diferencia entre observado y ajustado
valores_ajustados = p_ref_ols_result.fittedvalues
valores_efectivos = data_copper['chg_p_ref'].reindex(valores_ajustados.index)
mask = valores_efectivos.notna() & valores_ajustados.notna()
valores_ajustados = valores_ajustados[mask]
valores_efectivos = valores_efectivos[mask]
diferencia = valores_efectivos - valores_ajustados

# Establecer límites y ticks comunes
ymin = min(valores_efectivos.min(), valores_ajustados.min(), -abs(diferencia).max())
ymax = max(valores_efectivos.max(), valores_ajustados.max(), abs(diferencia).max())
ylim = max(abs(ymin), abs(ymax))
yticks = np.round(np.linspace(-ylim, ylim, num = 9), 2)  # Mismos valores, redondeados

# Crear gráfico
fig, ax1 = plt.subplots(dpi = 180)

# Curvas principales
ax1.plot(valores_ajustados, color = color_rojo_principal, label = 'Estimado (ajustado)')
ax1.plot(valores_efectivos, color = color_azul_principal, linestyle = '--', label = 'Efectivo (observado)')
ax1.set_title("Precio de referencia del cobre: efectivo vs estimado", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 15)
ax1.set_xlabel("Año", fontsize=12, fontfamily='Tw Cen MT')
ax1.set_ylabel("Variación en precio de referencia", fontsize=12, fontfamily='Tw Cen MT')
ax1.set_ylim(-ylim, ylim)
ax1.set_yticks(yticks)
ax1.spines[['top', 'right']].set_visible(False)
ax1.tick_params(axis = 'both', direction = 'out')
ax1.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax1.tick_params(labelsize = 8)
ax1.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)

# Eje secundario: diferencia como barras
ax2 = ax1.twinx()
ax2.bar(valores_ajustados.index, diferencia, width = 250, color = color_gris_secundario, alpha = 0.3, edgecolor = 'black', linewidth = 0.5, label = 'Diferencia (efectivo - estimado)')
ax2.set_ylabel("Diferencia", fontsize = 12, fontfamily = 'Tw Cen MT')
ax2.set_ylim(-ylim, ylim)
ax2.set_yticks(yticks)
ax2.tick_params(axis = 'y', labelsize = 8)
ax2.spines[['top', 'left']].set_visible(False)

# Combinar leyendas
handles1, labels1 = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(handles1 + handles2, labels1 + labels2, loc = 'lower left', fontsize=9, frameon=False)

plt.tight_layout()
plt.show()

# Simulación optimizada de brecha del precio del cobre
sim_cobre = np.empty([int(n_h/4), 3, n_b])

for i in range(0, n_b):
    for j in range(0, int(n_h/4)):
        # Calcula precio nominal (0) y de referencia (1)
        if j == 0:
            sim_cobre[j, 0, i] = tray_anual[j, 3, i]*(1 + infl_usa0)
            sim_cobre[j, 1, i] = p_ref_t
        elif j == 1:
            sim_cobre[j, 0, i] = tray_anual[j, 3, i]*(1 + infl_usa0)*(1 + infl_usa1)
            if usar_pr_1 == 1:
                sim_cobre[j, 1, i] = p_ref_1
            elif usar_pr_1 == 0:
                sim_cobre[j, 1, i] = sim_cobre[j - 1, 1, i] + \
                    (pref_params['Intercept'] + \
                     pref_params['chg_c_n_l1']*(sim_cobre[j - 1, 0, i] - data_anual['copper_n'].iloc[-1]) + \
                     pref_params['chg_c_n_l2']*data_anual['chg_c_n'].iloc[-1] + \
                     pref_params['chg_c_n_l3']*data_anual['chg_c_n'].iloc[-2] + \
                     pref_params['chg_c_n_l4']*data_anual['chg_c_n'].iloc[-3])  
        elif j == 2:
            sim_cobre[j, 0, i] = tray_anual[j, 3, i]*(1 + infl_usa0)*(1 + infl_usa1)*(1 + infl_usa2)
            sim_cobre[j, 1, i] = sim_cobre[j - 1, 1, i] + \
                (pref_params['Intercept'] + \
                 pref_params['chg_c_n_l1']*(sim_cobre[j - 1, 0, i] - sim_cobre[j - 2, 0, i]) + \
                 pref_params['chg_c_n_l2']*(sim_cobre[j - 2, 0, i] - data_anual['copper_n'].iloc[-1]) + \
                 pref_params['chg_c_n_l3']*data_anual['chg_c_n'].iloc[-1] + \
                 pref_params['chg_c_n_l4']*data_anual['chg_c_n'].iloc[-2])
        elif j == 3:
            sim_cobre[j, 0, i] = tray_anual[j, 3, i]*(1 + infl_usa0)*(1 + infl_usa1)*(1 + infl_usa2)*(1 + infl_usa3)
            sim_cobre[j, 1, i] = sim_cobre[j - 1, 1, i] + \
                (pref_params['Intercept'] + \
                 pref_params['chg_c_n_l1']*(sim_cobre[j - 1, 0, i] - sim_cobre[j - 2, 0, i]) + \
                 pref_params['chg_c_n_l2']*(sim_cobre[j - 2, 0, i] - sim_cobre[j - 3, 0, i]) + \
                 pref_params['chg_c_n_l3']*(sim_cobre[j - 3, 0, i] - data_anual['copper_n'].iloc[-1]) + \
                 pref_params['chg_c_n_l4']*data_anual['chg_c_n'].iloc[-1])
        elif j == 4:
            sim_cobre[j, 0, i] = tray_anual[j, 3, i]*(1 + infl_usa0)*(1 + infl_usa1)*(1 + infl_usa2)*(1 + infl_usa3)*(1 + infl_usaLR)
            sim_cobre[j, 1, i] = sim_cobre[j - 1, 1, i] + \
                (pref_params['Intercept'] + \
                 pref_params['chg_c_n_l1']*(sim_cobre[j - 1, 0, i] - sim_cobre[j - 2, 0, i]) + \
                 pref_params['chg_c_n_l2']*(sim_cobre[j - 2, 0, i] - sim_cobre[j - 3, 0, i]) + \
                 pref_params['chg_c_n_l3']*(sim_cobre[j - 3, 0, i] - sim_cobre[j - 4, 0, i]) + \
                 pref_params['chg_c_n_l4']*(sim_cobre[j - 4, 0, i] - data_anual['copper_n'].iloc[-1]))
        else:
            sim_cobre[j, 0, i] = tray_anual[j, 3, i]*(1 + infl_usa0)*(1 + infl_usa1)*(1 + infl_usa2)*(1 + infl_usa3)*(1 + infl_usaLR)**(j - 4)
            sim_cobre[j, 1, i] = sim_cobre[j - 1, 1, i] + \
                (pref_params['Intercept'] + \
                 pref_params['chg_c_n_l1']*(sim_cobre[j - 1, 0, i] - sim_cobre[j - 2, 0, i]) + \
                 pref_params['chg_c_n_l2']*(sim_cobre[j - 2, 0, i] - sim_cobre[j - 3, 0, i]) + \
                 pref_params['chg_c_n_l3']*(sim_cobre[j - 3, 0, i] - sim_cobre[j - 4, 0, i]) + \
                 pref_params['chg_c_n_l4']*(sim_cobre[j - 4, 0, i] - sim_cobre[j - 5, 0, i]))
        
        # Calcula brecha del precio del cobre
        sim_cobre[j, 2, i] = sim_cobre[j, 0, i] - sim_cobre[j, 1, i]
        
tray_anual = np.concatenate((tray_anual, sim_cobre[:, 2:3, :]), axis = 1)

# Extraer percentiles desde trayectorias anuales
perc = np.percentile(tray_anual[:, 9, :], np.arange(10, 100, 10), axis = 1)
mediana_simulacion = perc[4]

# Plot principal
fig, ax = plt.subplots(dpi = 180)
ax.plot(años_sim, mediana_simulacion, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")
    
# Rellenar bandas de confianza
for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
    ax.fill_between(años_sim, perc[j], perc[8 - j], color = color_rojo_secundario, alpha = alpha, label = f"Rango P{10*(j + 1)} - P{90 - 10*j}")

ax.set_title("Brecha del precio del cobre ($USc/lb)", fontsize = 16, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 20)
ax.set_xlabel("Años", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_xticks(años_sim)
ax.set_ylabel("Centavos de dólar", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)
ax.legend(loc = "lower left", fontsize = 10, frameon = False)

plt.tight_layout()
plt.show()

# Inicialización de matriz de simulación
sim_pibr = np.empty((int(n_h/4), 7, n_b))
horizonte = int(n_h/4)
gdpr_hist = data_anual['gdpr'].to_numpy()
ultimo_gdpr = gdpr_hist[-1]

for i in range(n_b):
    crecimiento = tray_anual[:, 0, i]/100
    # Calcular PIB real proyectado
    sim_pibr[0, 0, i] = ultimo_gdpr*(1 + crecimiento[0])
    for j in range(1, horizonte):
        sim_pibr[j, 0, i] = sim_pibr[j - 1, 0, i]*(1 + crecimiento[j])

    # Calcular tendencia con filtro HP
    serie_total = np.concatenate((gdpr_hist, sim_pibr[:, 0, i]))
    _, tendencia = sm.tsa.filters.hpfilter(serie_total, lamb=6.25)
    tendencia_sim = tendencia[-horizonte:]

    # Guardar en matriz
    sim_pibr[:, 1, i] = tendencia_sim
    brecha = sim_pibr[:, 0, i]/tendencia_sim - 1
    sim_pibr[:, 2, i] = brecha
    sim_pibr[:, 3, i] = np.where(brecha > 0, brecha, 0)
    sim_pibr[:, 4, i] = np.where(brecha <= 0, brecha, 0)
    sim_pibr[:, 5, i] = (crecimiento < 0).astype(int)

# Concatenar brechas positivas, negativas y crisis
tray_anual = np.concatenate((tray_anual, sim_pibr[:, 3:6, :]), axis = 1)

# Extraer percentiles desde trayectorias anuales
perc = np.percentile(sim_pibr[:, 2, :]*100, np.arange(10, 100, 10), axis = 1)
mediana_simulacion = perc[4]

# Plot principal
fig, ax = plt.subplots(dpi = 180)
ax.plot(años_sim, mediana_simulacion, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")
    
# Rellenar bandas de confianza
for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
    ax.fill_between(años_sim, perc[j], perc[8 - j], color = color_rojo_secundario, alpha = alpha, label = f"Rango P{10*(j + 1)} - P{90 - 10*j}")

ax.set_title("Brecha del PIB (% desviación)", fontsize = 16, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 20)
ax.set_xlabel("Años", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_xticks(años_sim)
ax.set_ylabel("Porcentaje", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out')
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.tick_params(labelsize = 8)
ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)
ax.legend(loc = "lower left", fontsize = 10, frameon = False)
plt.tight_layout()
plt.show()

# Inflación promedio últimos 10 años
inf_prom = np.mean(data_anual["inf_ipc"][-10:])

# Tasa de interés UF efectiva del último año
j_uf_0 = j_d_0 - ((1 - prop - prop_uf)/(1 - prop))*inf_prom/100  

# Tasa de interés pesos efectiva del último año
j_p_0 = j_d_0 + (prop_uf/(1 - prop))*inf_prom/100              

# Simulación:
sim_tasas = np.empty([int(n_h/4), 6, n_b])
for i in range(0, n_b):
    for j in range(0, int(n_h/4)):
        # Tasa int. nom. en pesos chilenos
        sim_tasas[j, 0, i] = np.where(tray_anual[j, 4, i] + tray_anual[j, 6, i] > 0, tray_anual[j, 4, i] + tray_anual[j, 6, i], 0) 

        # Tasa int. en UF
        sim_tasas[j, 1, i] = tray_anual[j, 4, i]                                         

        # Tasa int. nom. mon. extranjera
        sim_tasas[j, 2, i] = np.where(tray_anual[j, 5, i] + tray_anual[j, 6, i] > 0, tray_anual[j, 5, i] + tray_anual[j, 6, i], 0) 
        
        if j == 0:
            # Tasa de interés de la deuda en mon. local
            sim_tasas[j, 3, i] = (1 - 1/dur_p)*j_p_0*100 + (1/dur_p)*sim_tasas[j, 0, i] 

            # Tasa de interés de la deuda en UF
            sim_tasas[j, 4, i] = (1 - 1/dur_uf)*j_uf_0*100 + (1/dur_uf)*sim_tasas[j, 1, i] 

            # Tasa de interés de la deuda en mon. extranjera
            sim_tasas[j, 5, i] = (1 - 1/dur_f)*j_e_0*100 + (1/dur_f)*sim_tasas[j, 2, i]    
        else:
            sim_tasas[j, 3, i] = (1 - 1/dur_p)*sim_tasas[j - 1, 3, i] + (1/dur_p)*sim_tasas[j, 0, i]
            sim_tasas[j, 4, i] = (1 - 1/dur_uf)*sim_tasas[j - 1, 4, i] + (1/dur_uf)*sim_tasas[j, 1, i]
            sim_tasas[j, 5, i] = (1 - 1/dur_f)*sim_tasas[j - 1, 5, i] + (1/dur_f)*sim_tasas[j, 2, i]
            
tray_anual = np.concatenate((tray_anual, sim_tasas[:, 3:6, :]), axis = 1)

# Simulación de trayectorias de balance primario, deuda, activos y gastos por intereses

# Para guardar resultados deuda, balance primario e intereses
sim_bp_d = np.empty([int(n_h/4), 4, n_b]) 

# Para guardar resultados activos
sim_act = np.empty([int(n_h/4), 10, n_b])

for i in range(0, n_b):
    for j in range(0, int(n_h/4)):
        if j == 0:
            # Simular activos:
            
            # FRP: ----
            sim_act[j, 0, i] = max(0, min(data_anual['b_pib'].iloc[-1], 0.005))                                  # Aporte
            sim_act[j, 1, i] = 0.001 if data_anual['frp_p'].iloc[-1] > 0.001 else data_anual['frp_p'].iloc[-1]   # Retiro
            sim_act[j, 2, i] = data_anual['frp_p'].iloc[-1] + sim_act[j, 0, i] - sim_act[j, 1, i]                # Saldo final
            
            # FEES: ----
            # Aporte
            if data_anual['b_pib'].iloc[-1] > max(0, min(data_anual['b_pib'].iloc[-1], 0.005)):
                sim_act[j, 3, i] = data_anual['b_pib'].iloc[-1] - max(0, min(data_anual['b_pib'].iloc[-1], 0.005)) \
                if data_anual['fees_p'].iloc[-1] + data_anual['b_pib'].iloc[-1] - max(0,min(data_anual['b_pib'].iloc[-1], 0.005)) < 0.07 \
                else 0.07 - data_anual['fees_p'].iloc[-1]
            else:
                sim_act[j, 3, i] = 0

            # Saldo final antes de retiros
            sim_act[j, 4, i] = data_anual['fees_p'].iloc[-1] + sim_act[j, 3, i]   
            
            # Retiro
            if tray_anual[j, 12, i] == 1 and sim_pibr[j, 2, i] <= -0.03:
                sim_act[j, 5, i] = 2*tray_anual[j, 0, i] if sim_act[j, 4, i] >= -2*tray_anual[j, 0, i] else -sim_act[j, 4, i]
            else:
                sim_act[j, 5, i] = 0

            # Saldo final
            sim_act[j, 6, i] = sim_act[j, 4, i] + sim_act[j, 5, i]   
            
            # OATP: ----
            sim_act[j, 7, i] = data_anual['oatp_p'].iloc[-1] + oatp_params['oatp_mod']*(np.mean(data_anual['oatp2_p'][-10:]) - data_anual['oatp_p'].iloc[-1]) +\
                           sim_shock_oatp[j, 0, i]
            sim_act[j, 7, i] = sim_act[j, 7, i] if sim_act[j, 7, i] > 0 else 0
            
            # Total activos: ----
            sim_act[j, 8, i] = sim_act[j, 2, i] + sim_act[j, 6, i] + sim_act[j, 7, i]
            
            # Simular balance primario: ----
            if usar_bp_0 == 0:
                sim_bp_d[j, 0, i] = frf_params['const'] + frf_params['ret_fees_p']*sim_act[j, 5, i] +\
                                    frf_params[['d_p_lag', 'd_p_lag_squared', 'd_p_lag_cubed']]@[data_anual['d_p'].iloc[-1], data_anual['d_p'].iloc[-1]**2, data_anual['d_p'].iloc[-1]**3]+ \
                                    frf_params['brecha_cobre']*tray_anual[j, 9, i] +\
                                    frf_params['brecha_pos_int']*tray_anual[j, 10, i] +\
                                    frf_params['brecha_neg_int']*tray_anual[j, 11, i] +\
                                    frf_params['regla_fiscal'] +\
                                    frf_params['crisis']*tray_anual[j, 12, i] +\
                                    frf_params['ar.L1']*frf.resid.iloc[-1] + tray_anual[j, 7, i]
            elif usar_bp_0 == 1:
                sim_bp_d[j, 0, i] = bp_0 + (frf_params['ar.L1']*frf.resid.iloc[-1] + tray_anual[j, 7, i])*(1 - data_trim.tail().index[0].month/12)
          
            mu = frf_params['ar.L1']*frf.resid.iloc[-1] + tray_anual[j, 7, i]
            
            # Simular deuda: ----
            sim_bp_d[j, 1, i] =((prop*(1 + tray_anual[j, 15, i]/100)*(1 + tray_anual[j, 2, i]/100)+\
                                prop_uf*(1 + tray_anual[j, 14, i]/100)*(1 + tray_anual[j, 6, i]/100)+\
                                (1 - prop - prop_uf)*(1 + tray_anual[j, 13, i]/100))/ \
                                ((1 + tray_anual[j, 0, i]/100)*(1 + tray_anual[j, 1, i]/100)))* \
                                data_anual['d_p'].iloc[-1] -\
                                sim_bp_d[j, 0, i] + tray_anual[j, 8, i] +\
                                (sim_act[j, 8, i] - (data_anual['fees_p'].iloc[-1] + data_anual['frp_p'].iloc[-1] + data_anual['oatp_p'].iloc[-1])/((1 + tray_anual[j, 0, i]/100)*(1 + tray_anual[j, 1, i]/100)))
                            
            # Simular gastos por intereses: ----
            sim_bp_d[j, 2, i] = ((prop*tray_anual[j, 15, i]/100*(1 + tray_anual[j, 2, i]/100) + prop_uf*tray_anual[j, 14, i]/100*(1 + tray_anual[j, 6, i]/100) + (1 - prop - prop_uf)*tray_anual[j, 13, i]/100)/ \
                                ((1 + tray_anual[j, 0, i]/100)*(1 + tray_anual[j, 1, i]/100)))* \
                                data_anual['d_p'].iloc[-1]
        else:
            # Simular activos
            
            # FRP: ----
            sim_act[j, 0, i] = max(0, min(sim_bp_d[j - 1, 3, i], 0.005))                       # Aporte
            sim_act[j, 1, i] = 0.001 if sim_act[j - 1, 1, i] > 0.001 else sim_act[j - 1, 1, i] # Retiro
            sim_act[j, 2, i] = sim_act[j - 1, 2, i] + sim_act[j, 0, i] - sim_act[j, 1, i]      # Saldo final
            
            # FEES: ----
            # Aporte
            if sim_bp_d[j - 1, 3, i] > max(0, min(sim_bp_d[j - 1, 3, i], 0.005)):
                sim_act[j, 3, i] = sim_bp_d[j - 1, 3, i] - max(0, min(sim_bp_d[j - 1, 3, i], 0.005)) \
                if sim_act[j - 1, 4, i] + sim_bp_d[j - 1, 3, i] - max(0, min(sim_bp_d[j - 1, 3, i], 0.005)) < 0.07 else 0.07 - sim_act[j - 1, 4, i]
            else:
                sim_act[j, 3, i] = 0

            # Saldo final antes de retiros
            sim_act[j, 4, i] = sim_act[j - 1, 6, i] + sim_act[j, 3, i]
            
            # Retiro
            if tray_anual[j, 12, i] == 1 and sim_pibr[j, 2, i] <= -0.03:
                sim_act[j, 5, i ]= 2*tray_anual[j, 0, i] if sim_act[j, 4, i] >= -2*tray_anual[j, 0, i] else -sim_act[j, 4, i]
            else:
                sim_act[j, 5, i] = 0

            # Saldo final
            sim_act[j, 6, i] = sim_act[j, 4, i] + sim_act[j, 5, i]   
            
            # OATP: ----
            sim_act[j, 7, i] = sim_act[j - 1, 7, i] + oatp_params['oatp_mod']*(np.mean(data_anual['oatp2_p'][-10:]) - sim_act[j - 1, 7, i]) +\
                           sim_shock_oatp[j, 0, i]
            sim_act[j, 7, i] = sim_act[j, 7, i] if sim_act[j, 7, i] > 0 else 0
            
            # Total activos: ----
            sim_act[j, 8, i] = sim_act[j, 2, i] + sim_act[j, 6, i] + sim_act[j, 7, i]
            
            # Simular balance primario: ----
            sim_bp_d[j, 0, i] = frf_params['const'] + frf_params['ret_fees_p']*sim_act[j, 5, i] +\
                                frf_params[['d_p_lag', 'd_p_lag_squared', 'd_p_lag_cubed']]@[sim_bp_d[j - 1, 1, i], sim_bp_d[j - 1, 1, i]**2, sim_bp_d[j - 1, 1, i]**3]+ \
                                frf_params['brecha_cobre']*tray_anual[j, 9, i] +\
                                frf_params['brecha_pos_int']*tray_anual[j, 10, i] +\
                                frf_params['brecha_neg_int']*tray_anual[j, 11, i] +\
                                frf_params['regla_fiscal'] +\
                                frf_params['crisis']*tray_anual[j, 12, i] +\
                                frf_params['ar.L1']*mu+tray_anual[j, 7, i]
            
            mu = frf_params['ar.L1']*mu + tray_anual[j, 7, i]
            
            # Simular deuda: ----
            sim_bp_d[j, 1, i] = ((prop*(1 + tray_anual[j, 15, i]/100)*(1 + tray_anual[j, 2, i]/100)+\
                                prop_uf*(1 + tray_anual[j, 14, i]/100)*(1 + tray_anual[j, 6, i]/100)+\
                                (1 - prop - prop_uf)*(1 + tray_anual[j, 13, i]/100))/ \
                                ((1 + tray_anual[j, 0, i]/100)*(1 + tray_anual[j, 1, i]/100)))* \
                                sim_bp_d[j - 1, 1, i] -\
                                sim_bp_d[j, 0, i] + tray_anual[j, 8, i] +\
                                (sim_act[j, 8, i] - sim_act[j - 1, 8, i]/((1 + tray_anual[j, 0, i]/100)*(1 + tray_anual[j, 1, i]/100)))
                            
            # Gastos por intereses: ---
            sim_bp_d[j, 2, i] = ((prop*tray_anual[j, 15, i]/100*(1 + tray_anual[j, 2, i]/100) + prop_uf*tray_anual[j, 14, i]/100*(1 + tray_anual[j, 6, i]/100) + (1 - prop - prop_uf)*tray_anual[j, 13, i]/100)/ \
                                ((1 + tray_anual[j, 0, i]/100)*(1 + tray_anual[j, 1, i]/100)))* \
                                sim_bp_d[j - 1, 1, i]
        # Balance total: ----
        sim_bp_d[j, 3, i] = sim_bp_d[j, 0, i] - sim_bp_d[j, 2, i]

# Simulación como porcentaje del PIB:
sim_bp_d = sim_bp_d*100

# Gráficos de fancharts de simulaciones de balance primario y deuda bruta
titles = ['Balance primario (% del PIB)',
          'Deuda bruta del Gobierno Central (% del PIB)',
          'Gastos por intereses (% del PIB)',
          'Balance total (% del PIB)']

# Escalamos a porcentaje
data_plot = data_anual[['bp_pib', 'd_p', 'gi_pib', 'b_pib']]*100

# Años históricos y de simulación
años = data_plot.index.year[-20:]
años_sim = np.arange(años[-1] + 1, años[-1] + 1 + int(n_h/4))
años_con = np.concatenate([[años[-1]], años_sim])

# Percentiles
percentiles = [p1, p2, p3, p4]
alphas = [0.1, 0.2, 0.3, 0.4]

fig, axes = plt.subplots(nrows = 2, ncols = 2, figsize = (12, 10), dpi = 180)
axes = axes.flatten()

for i in range(4):
    ax = axes[i]
    perc = np.percentile(sim_bp_d[:, i, :],
                         [p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1],
                         axis = 1)

    ax.plot(años, data_plot.iloc[-20:, i], color = 'black', linewidth = 1.5, label = 'Datos efectivos')
    ultimo_efectivo = data_plot.iloc[-1, i]
    mediana_simulacion = np.concatenate([[ultimo_efectivo], perc[4]])
    ax.plot(años_con, mediana_simulacion, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")

    for j, alpha in enumerate(alphas[::-1]):
        lower = np.concatenate([[ultimo_efectivo], perc[j]])
        upper = np.concatenate([[ultimo_efectivo], perc[8 - j]])
        ax.fill_between(años_con, lower, upper, color = color_rojo_secundario, alpha = alpha,
                        label = f"Rango {percentiles[::-1][j]} - {100 - percentiles[::-1][j]}")

    # Línea punteada vertical y texto
    ax.axvline(x = años_con[0], color = 'gray', linestyle = '--', linewidth = 1)
    ax.text(años_con[0], ax.get_ylim()[1]*0.95, "Proyecciones", rotation = 90, fontsize = 9,
            verticalalignment = 'top', color = 'gray', fontfamily = 'Tw Cen MT')

    ax.set_title(titles[i], fontsize = 12, fontweight = 'bold', fontfamily = 'Tw Cen MT')
    ax.set_xlabel("Años", fontsize = 10, fontfamily = 'Tw Cen MT')
    ax.set_xticks(np.arange(años[0], años_sim[-1] + 1, 2))
    ax.set_ylabel("Porcentaje del PIB", fontsize = 10, fontfamily = 'Tw Cen MT')
    ax.spines[['top', 'right']].set_visible(False)
    ax.tick_params(axis = 'both', direction = 'out', labelsize=8)
    ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
    if (perc < 0).any():
        ax.axhline(y = 0, color = 'black', linestyle = '-', linewidth = 0.5, alpha = 0.3)

# Leyenda en la última subfigura vacía
handles, labels = axes[0].get_legend_handles_labels()
axes[2].legend(handles, labels, loc = 'center left', fontsize = 9, frameon = False)

plt.tight_layout()
plt.subplots_adjust(top = 0.92)
plt.suptitle("Fan Charts de Balance y Deuda del Gobierno Central", fontsize = 14, fontweight = 'bold', fontfamily = 'Tw Cen MT')
plt.show()

# Preparar años simulados con frecuencia correcta
años_sim2 = pd.date_range(start = data_plot.index[-1] + pd.DateOffset(years = 1), periods = int(n_h / 4), freq = "YE-DEC")
años_sim_labels = años_sim2.year
años_con = np.concatenate([[data_plot.index.year[-1]], años_sim_labels])

# Calcular percentiles de deuda simulada y diferencia con la mediana
percentiles = [p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1]
perc_db = np.percentile(sim_bp_d[:, 1, :], percentiles, axis = 1)

# centrar respecto a la mediana
dif_mediana = perc_db - perc_db[4]  

# Asegurar índice datetime para match correcto
deuda_det.index = pd.to_datetime(deuda_det.index)
col_deuda = deuda_det.columns[0]
db_det = deuda_det.loc[deuda_det.index.isin(años_sim2), col_deuda]

# Asegurar dimensiones compatibles
if dif_mediana.shape[1] != len(db_det):
    dif_mediana = dif_mediana.T
per_db_det = dif_mediana + db_det.to_numpy()

# Preparar datos para gráfico
data_db = data_anual[['d_p']]*100
años = data_db.index.year[-20:]
ultimo_efectivo = data_db.iloc[-1, 0]

fig, ax = plt.subplots(dpi = 180)

# Línea de datos efectivos
ax.plot(años, data_db.iloc[-20:, 0], color = 'black', linewidth = 1.5, label = 'Datos efectivos')

# Mediana simulada
mediana = np.concatenate([[ultimo_efectivo], per_db_det[4]])
ax.plot(años_con, mediana, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")

# Fanchart con bucle
for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
    lower = np.concatenate([[ultimo_efectivo], per_db_det[j]])
    upper = np.concatenate([[ultimo_efectivo], per_db_det[8 - j]])
    ax.fill_between(años_con, lower, upper, color = color_rojo_secundario, alpha = alpha, label = f"Rango {percentiles[j]}-{percentiles[8 - j]}")

# Línea horizontal de referencia
nivel_prudente = 45
ax.axhline(nivel_prudente, color = color_azul_secundario, linestyle = '--', linewidth = 1)
ax.text(años[0], nivel_prudente + 0.8, "Nivel prudente de deuda (45% del PIB)", color = color_azul_secundario, fontsize = 10, va = 'bottom')

# Línea punteada vertical y texto
ax.axvline(x = años_con[0], color = 'gray', linestyle = '--', linewidth = 1)
ax.text(años_con[0], ax.get_ylim()[1]*0.5, "Proyecciones", rotation = 90, fontsize = 9,
            verticalalignment = 'top', color = 'gray', fontfamily = 'Tw Cen MT')

# Estética del gráfico
ax.set_title('Deuda bruta del Gobierno Central ("baseline centered") (% del PIB)',
             fontsize = 16, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 20)
ax.set_xlabel("Años", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_ylabel("% del PIB", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.legend(loc = "center left", fontsize = 10, frameon = False)

plt.tight_layout()
plt.savefig('fanchart_deuda_bruta_baseline_centered.jpg', transparent = True)
plt.show()

# Calcular percentiles de simulación de activos
percentiles = [p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1]
perc = np.percentile(sim_act[:, 8, :]*100, percentiles, axis = 1) 

# Transponer una sola vez para eficiencia
perc_T = perc.T  

# Preparar gráfico
fig, ax = plt.subplots(dpi = 180)

# Línea de la mediana
ax.plot(años_sim, perc_T[:, 4], color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")

# Relleno de fancharts en bucle
for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
    ax.fill_between(años_sim, perc_T[:, j], perc_T[:, 8 - j],
                    color = color_rojo_secundario, alpha = alpha,
                    label = f"Rango {percentiles[j]}-{percentiles[8 - j]}")

# Estética del gráfico
ax.set_title("Total activos del Tesoro Público (% del PIB)", fontsize = 16, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 20)
ax.set_xlabel("Años", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_ylabel("% del PIB", fontsize = 12, fontfamily = 'Tw Cen MT')  
ax.set_xticks(años_sim)
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.legend(loc = "upper left", fontsize = 8, frameon = False)

plt.tight_layout()
plt.show()

# Calcular deuda neta simulada (vectorizado)
sim_act[:, 9, :] = sim_bp_d[:, 1, :] - sim_act[:, 8, :]*100

# Percentiles
percentiles = [p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1]
perc = np.percentile(sim_act[:, 9, :], percentiles, axis = 1)
perc_T = perc.T 

# Preparar gráfico
fig, ax = plt.subplots(dpi = 180)

# Datos efectivos
ax.plot(años, data_anual['deuda_neta'].iloc[-20:]*100, color = 'black', linewidth = 1.5, label = 'Datos efectivos')

# Mediana simulada
ultimo_efectivo = data_anual['deuda_neta'].iloc[-1]*100
mediana_sim = np.concatenate([[ultimo_efectivo], perc_T[:, 4]])
ax.plot(años_con, mediana_sim, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")

# Fancharts (rango de percentiles)
for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
    lower = np.concatenate([[ultimo_efectivo], perc_T[:, j]])
    upper = np.concatenate([[ultimo_efectivo], perc_T[:, 8 - j]])
    ax.fill_between(años_con, lower, upper, color = color_rojo_secundario, alpha = alpha, label = f"Rango {percentiles[j]}-{percentiles[8 - j]}")

# Línea punteada vertical y texto
ax.axvline(x = años_con[0], color = 'gray', linestyle = '--', linewidth = 1)
ax.text(años_con[0], ax.get_ylim()[1]*0.5, "Proyecciones", rotation = 90, fontsize = 9,
            verticalalignment = 'top', color = 'gray', fontfamily = 'Tw Cen MT')

# Estética del gráfico
ax.set_title('Deuda neta del Gobierno Central (% del PIB)', fontsize = 16, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 20)
ax.set_xlabel("Años", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_ylabel("% del PIB", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_xticks(np.arange(años[0], años_sim[-1] + 1, 2))  
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.legend(loc = "upper left", fontsize = 10, frameon = False)

plt.tight_layout()
plt.savefig('fanchart_deuda_neta.jpg', transparent=True)
plt.show()

# Preparar años simulados con frecuencia correcta
años_sim2 = pd.date_range(start = data_plot.index[-1] + pd.DateOffset(years = 1), periods = int(n_h / 4), freq = "YE-DEC")
años_sim_labels = años_sim2.year
años_con = np.concatenate([[data_plot.index.year[-1]], años_sim_labels])

# Calcular percentiles de deuda simulada y diferencia con la mediana
percentiles = [p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1]
perc_dn = np.percentile(sim_act[:, 9, :], percentiles, axis = 1)

# centrar respecto a la mediana
dif_mediana = perc_dn - perc_dn[4]  

# Asegurar índice datetime para match correcto
deuda_det.index = pd.to_datetime(deuda_det.index)
col_deuda = deuda_det.columns[1]
dn_det = deuda_det.loc[deuda_det.index.isin(años_sim2), col_deuda]

# Asegurar dimensiones compatibles
if dif_mediana.shape[1] != len(dn_det):
    dif_mediana = dif_mediana.T
per_dn_det = dif_mediana + dn_det.to_numpy()

# Preparar datos para gráfico
data_dn = data_anual[['deuda_neta']]*100
años = data_dn.index.year[-20:]
ultimo_efectivo = data_dn.iloc[-1, 0]

fig, ax = plt.subplots(dpi = 180)

# Línea de datos efectivos
ax.plot(años, data_dn.iloc[-20:, 0], color = 'black', linewidth = 1.5, label = 'Datos efectivos')

# Mediana simulada
mediana = np.concatenate([[ultimo_efectivo], per_dn_det[4]])
ax.plot(años_con, mediana, color = color_rojo_principal, linestyle = '--', linewidth = 2, label = "Mediana simulación")

# Fanchart con bucle
for j, alpha in zip(range(3, -1, -1), [0.4, 0.3, 0.2, 0.1]):
    lower = np.concatenate([[ultimo_efectivo], per_dn_det[j]])
    upper = np.concatenate([[ultimo_efectivo], per_dn_det[8 - j]])
    ax.fill_between(años_con, lower, upper, color = color_rojo_secundario, alpha = alpha, label = f"Rango {percentiles[j]}-{percentiles[8 - j]}")

# Línea punteada vertical y texto
ax.axvline(x = años_con[0], color = 'gray', linestyle = '--', linewidth = 1)
ax.text(años_con[0], ax.get_ylim()[1]*0.5, "Proyecciones", rotation = 90, fontsize = 9,
            verticalalignment = 'top', color = 'gray', fontfamily = 'Tw Cen MT')

# Estética del gráfico
ax.set_title('Deuda neta del Gobierno Central ("baseline centered") (% del PIB)',
             fontsize = 16, fontweight = 'bold', fontfamily = 'Tw Cen MT', pad = 20)
ax.set_xlabel("Años", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_ylabel("% del PIB", fontsize = 12, fontfamily = 'Tw Cen MT')
ax.set_xticks(np.arange(años[0], años_sim_labels[-1] + 1, 2))
ax.spines[['top', 'right']].set_visible(False)
ax.tick_params(axis = 'both', direction = 'out', labelsize = 8)
ax.grid(axis = 'y', color = color_gris_principal, alpha = 0.3, linestyle = '--')
ax.legend(loc = "upper left", fontsize = 10, frameon = False)

plt.tight_layout()
plt.savefig('fanchart_deuda_neta_MED_baseline_centered.jpg', transparent = True)
plt.show()

# Media y mediana anual de la proyección de deuda, balance primario, gastos por intereses y balance total para horizonte de proyección.
resultados = pd.DataFrame()
resultados['mediana deuda'] = np.percentile(sim_bp_d[:, 1, :], 50, axis = 1).reshape(-1)
resultados['media deuda'] = np.mean(sim_bp_d[:, 1, :], axis = 1)
resultados['mediana bp'] = np.percentile(sim_bp_d[:, 0, :], 50, axis = 1).reshape(-1)
resultados['media bp'] = np.mean(sim_bp_d[:, 0, :], axis = 1)
resultados['mediana gint'] = np.percentile(sim_bp_d[:, 2, :], 50, axis = 1).reshape(-1)
resultados['media gint'] = np.mean(sim_bp_d[:, 2, :], axis = 1)
resultados['mediana bal'] = np.percentile(sim_bp_d[:, 3, :], 50, axis = 1).reshape(-1)
resultados['media bal'] = np.mean(sim_bp_d[:, 3, :], axis = 1)
resultados.index = años_sim2

# Probabilidad de superar el nivel prudente de deuda al final del horizonte de proyección.
sup_prud = (sim_bp_d[:, 1, :] > 45).astype(int)
resultados['prob>n.prud'] = np.round(np.sum(sup_prud, axis = 1)/n_b*100, 1)

# Estadísticos de la deuda
resultados['desv.est. deuda'] = np.std(sim_bp_d[:, 1, :], axis = 1)
resultados['skewness deuda'] = skew(sim_bp_d[:, 1, :], axis = 1)
resultados['kurtosis deuda'] = kurtosis(sim_bp_d[:, 1, :], axis = 1)

# Percentiles deuda modelo MED
perc_db_med = pd.DataFrame(perc_db.T)
perc_db_med.index = años_sim2
perc_db_med.columns = p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1

# Percentiles deuda modelo MED con mediana de modelo deterministico
perc_db_det = pd.DataFrame(per_db_det.T)
perc_db_det.index = años_sim2
perc_db_det.columns = p1, p2, p3, p4, 50, 100 - p4, 100 - p3, 100 - p2, 100 - p1

# Mostrar las últimas filas con formato
resultados.style\
    .format(precision = 1, thousands = ".", decimal = ",")\
    .set_table_styles([{'selector': 'caption',
                        'props': [('caption-side', 'top'), ('font-size', '14px'), ('font-weight', 'bold')]}])\
    .set_properties(**{'text-align': 'center'})\
    .set_table_attributes('class="table table-striped"')

# Exportamos resultados a Excel
with pd.ExcelWriter('resultados modelo MED.xlsx') as writer:
    resultados.to_excel(excel_writer = writer, sheet_name = 'Resumen resultados')
    perc_db_med.to_excel(excel_writer = writer, sheet_name = 'Percentiles deuda MED')
    perc_db_det.to_excel(excel_writer = writer, sheet_name = 'Percentiles deuda mod det')
    
    workbook = writer.book
    for sheet_name in writer.sheets:
        sheet = writer.sheets[sheet_name]
        
        # Customizar los headers
        header_fill = PatternFill(start_color = "07434B", end_color = "07434B", fill_type = "solid")
        header_font = Font(bold = True, color = "FFFFFF")  # Texto en negrita
        
        for col_num, cell in enumerate(sheet[1], start = 1):  # Primera fila (títulos)
            cell.fill = header_fill
            cell.font = header_font
            
        # Customizar el resto de las celdas
        for row_idx, row in enumerate(sheet.iter_rows(min_row = 2), start = 2):  # Excluye fila con títulos
            for col_idx, cell in enumerate(row, start = 1):
                if col_idx == 1:  # Primrera columna (Index)
                    cell.fill = PatternFill(start_color = "CFDDE2", end_color = "CFDDE2", fill_type = "solid")
                else:  # Resto de las columnas
                    cell.fill = PatternFill(start_color = "E8EEF0", end_color = "E8EEF0", fill_type = "solid")
        
        for col_idx, cell in enumerate(sheet[1], start = 1):  # Itera sólo por los títulos (Headers)
            if col_idx > 1:  # Saltarse la primera columna
                header_length = len(str(cell.value))  # Obtener el largo del header
                adjusted_width = header_length + 2  # Ajustar el ancho
                sheet.column_dimensions[get_column_letter(col_idx)].width = adjusted_width

	gdpr_cl	def_chile	ipc_cl	tc_clp_usd	btu10	embi_cl	treasury_r10	copper_nom	gdpr_usa	cpi_usa
fecha
2024-03-31 00:00:00	52.099,9	145,9	102,7	945,5	2,5	135,1	1,9	3,8	23.053,5	312,1
2024-06-30 00:00:00	51.915,3	147,3	103,4	935,4	2,8	120,8	2,1	4,4	23.223,9	313,1
2024-09-30 00:00:00	52.734,1	148,8	104,5	931,6	2,6	128,0	1,8	4,2	23.400,3	314,9
2024-12-31 00:00:00	52.979,5	151,8	105,6	961,6	2,3	116,9	2,0	4,2	23.542,3	317,6
2025-03-31 00:00:00	53.342,2	153,5	107,7	964,0	2,6	125,0	2,1	4,2	23.526,1	319,6

	bal	bp	gint	deuda	orc	gdpn	gdpr	g_tend	copper_n	copper_ref	ipc	tcn_a	fees	frp	oatp	oatp2	ret_fees
fecha
2021-12-31 00:00:00	-18.497.666,2	-16.439.782,4	2.057.883,8	87.262.776,5	5.985.833,3	239.418,1	199.170,1	0,0	4,2	2,9	86,7	850,2	2.457,2	7.472,9	5.905,8	5.703,4	-6.196,8
2022-12-31 00:00:00	2.958.526,7	5.598.134,0	2.639.607,3	99.721.087,2	7.842.082,2	263.065,4	203.460,2	0,0	4,0	3,3	97,7	859,5	7.514,2	6.475,3	5.630,5	5.425,2	-0,1
2023-12-31 00:00:00	-6.718.708,5	-3.746.108,1	2.972.600,4	111.094.685,1	4.916.125,7	281.857,5	204.521,0	0,0	3,8	3,7	101,0	884,6	6.030,1	8.638,6	1.236,9	1.236,9	0,0
2024-12-31 00:00:00	-8.880.648,3	-5.110.858,6	3.769.789,8	129.794.268,2	6.324.650,1	311.630,9	209.929,1	0,0	4,1	3,9	105,6	992,1	3.618,2	9.378,3	1.106,2	1.106,2	-1.800,0
2025-12-31 00:00:00	nan	nan	nan	nan	nan	nan	nan	0,0	nan	4,1	nan	nan	nan	nan	nan	nan	nan

	db_pib_det	dn_pib_det
fecha
2070-12-31 00:00:00	44,4	40,8
2071-12-31 00:00:00	44,4	40,8
2072-12-31 00:00:00	44,4	40,8
2073-12-31 00:00:00	44,4	40,8
2074-12-31 00:00:00	44,4	40,8

	gdpr_cl	def_chile	ipc_cl	tc_clp_usd	btu10	embi_cl	treasury_r10	copper_nom	gdpr_usa	cpi_usa	i_usa10_mas_EMBI	copper_r	g_real	def_chile_chg	inf_cl	dep_tcn	embi_cl_chg	copper_nom_chg	g_real_usa	inf_usa
fecha
2024-03-31 00:00:00	52.099,9	145,9	102,7	945,5	2,5	135,1	1,9	3,8	23.053,5	312,1	3,2	3,8	1,4	3,8	1,6	5,5	-6,2	3,4	0,4	1,1
2024-06-30 00:00:00	51.915,3	147,3	103,4	935,4	2,8	120,8	2,1	4,4	23.223,9	313,1	3,3	4,4	-0,4	1,0	0,7	-1,1	-10,5	15,5	0,7	0,3
2024-09-30 00:00:00	52.734,1	148,8	104,5	931,6	2,6	128,0	1,8	4,2	23.400,3	314,9	3,1	4,1	1,6	1,0	1,1	-0,4	5,9	-5,6	0,8	0,5
2024-12-31 00:00:00	52.979,5	151,8	105,6	961,6	2,3	116,9	2,0	4,2	23.542,3	317,6	3,1	4,0	0,5	2,0	1,0	3,2	-8,7	-0,3	0,6	0,9
2025-03-31 00:00:00	53.342,2	153,5	107,7	964,0	2,6	125,0	2,1	4,2	23.526,1	319,6	3,3	4,1	0,7	1,1	2,0	0,2	6,9	1,8	-0,1	0,6

	bal	bp	gint	deuda	orc	gdpn	gdpr	g_tend	copper_n	copper_ref	ipc	tcn_a	fees	frp	oatp	oatp2	ret_fees	g	inf_ipc	chg_c_n	chg_p_ref	bp_pib	b_pib	gi_pib	orc_pib	d_p	d_p_lag	d_p_lag_squared	d_p_lag_cubed	fees_p	frp_p	oatp_p	oatp2_p	ret_fees_p	deuda_neta	crisis	regla_fiscal	brecha_cobre
fecha
2021-12-31 00:00:00	-18.497.666,2	-16.439.782,4	2.057.883,8	87.262.776,5	5.985.833,3	239.418,1	199.170,1	0,0	4,2	2,9	86,7	850,2	2.457,2	7.472,9	5.905,8	5.703,4	-6.196,8	0,1	7,2	1,4	0,0	-0,1	-0,1	0,0	0,0	0,4	0,3	0,1	0,0	0,0	0,0	0,0	0,0	-0,0	0,3	0	1	1,3
2022-12-31 00:00:00	2.958.526,7	5.598.134,0	2.639.607,3	99.721.087,2	7.842.082,2	263.065,4	203.460,2	0,0	4,0	3,3	97,7	859,5	7.514,2	6.475,3	5.630,5	5.425,2	-0,1	0,0	12,8	-0,2	0,4	0,0	0,0	0,0	0,0	0,4	0,4	0,1	0,0	0,0	0,0	0,0	0,0	-0,0	0,3	0	1	0,7
2023-12-31 00:00:00	-6.718.708,5	-3.746.108,1	2.972.600,4	111.094.685,1	4.916.125,7	281.857,5	204.521,0	0,0	3,8	3,7	101,0	884,6	6.030,1	8.638,6	1.236,9	1.236,9	0,0	0,0	3,4	-0,2	0,4	-0,0	-0,0	0,0	0,0	0,4	0,4	0,1	0,1	0,0	0,0	0,0	0,0	0,0	0,3	0	1	0,1
2024-12-31 00:00:00	-8.880.648,3	-5.110.858,6	3.769.789,8	129.794.268,2	6.324.650,1	311.630,9	209.929,1	0,0	4,1	3,9	105,6	992,1	3.618,2	9.378,3	1.106,2	1.106,2	-1.800,0	0,0	4,5	0,3	0,1	-0,0	-0,0	0,0	0,0	0,4	0,4	0,2	0,1	0,0	0,0	0,0	0,0	-0,0	0,4	0	1	0,3
2025-12-31 00:00:00	nan	nan	nan	nan	nan	nan	nan	0,0	nan	4,1	nan	nan	nan	nan	nan	nan	nan	nan	nan	nan	0,2	nan	nan	nan	nan	nan	0,4	0,2	0,1	nan	nan	nan	nan	nan	nan	0	1	0,0

Consejo Fiscal Autónomo de Chile (CFA)¶

Modelo Estocástico de la Deuda bruta del Gobierno Central de Chile (MED)¶

Versión 1.0¶

I. PREÁMBULO:¶

II. CARGA DE DATOS: importar bases de datos trimestrales y anuales¶

III. CALIBRACIÓN: parámetros del modelo¶

IV. TRANSFORMACIONES: preparación de bases de datos trimestrales y anuales¶

V. VISUALIZACIÓN DE DATOS: gráficos de las variables usadas en el modelo VAR¶

VI. ESTIMACIÓN: estimación del modelo VAR¶

VII. SIMULACIONES: análisis de modelo VAR¶

VIII. PROYECCIÓN DEUDA¶

	bal	bp	gint	deuda	orc	gdpn	gdpr	g_tend	copper_n	copper_ref	ipc	tcn_a	fees	frp	oatp	oatp2	ret_fees	g	inf_ipc	chg_c_n	chg_p_ref	bp_pib	b_pib	gi_pib	orc_pib	d_p	d_p_lag	d_p_lag_squared	d_p_lag_cubed	fees_p	frp_p	oatp_p	oatp2_p	ret_fees_p	deuda_neta	crisis	regla_fiscal	brecha_cobre	gdpr_gap	brecha_pos	brecha_pos_int	brecha_neg_int
fecha
2020-12-31 00:00:00	-14.642.921,5	-12.705.786,7	1.937.134,8	65.167.461,5	-1.756.088,0	201.257,7	178.924,9	0,0	2,8	2,9	80,9	711,2	8.955,2	10.156,8	4.391,5	4.189,4	-4.090,0	-0,1	3,0	0,1	-0,1	-0,1	-0,1	0,0	-0,0	0,3	0,3	0,1	0,0	0,0	0,0	0,0	0,0	-0,0	0,2	1	1	-0,1	-0,1	0	0,0	-0,1
2021-12-31 00:00:00	-18.497.666,2	-16.439.782,4	2.057.883,8	87.262.776,5	5.985.833,3	239.418,1	199.170,1	0,0	4,2	2,9	86,7	850,2	2.457,2	7.472,9	5.905,8	5.703,4	-6.196,8	0,1	7,2	1,4	0,0	-0,1	-0,1	0,0	0,0	0,4	0,3	0,1	0,0	0,0	0,0	0,0	0,0	-0,0	0,3	0	1	1,3	0,0	1	0,0	0,0
2022-12-31 00:00:00	2.958.526,7	5.598.134,0	2.639.607,3	99.721.087,2	7.842.082,2	263.065,4	203.460,2	0,0	4,0	3,3	97,7	859,5	7.514,2	6.475,3	5.630,5	5.425,2	-0,1	0,0	12,8	-0,2	0,4	0,0	0,0	0,0	0,0	0,4	0,4	0,1	0,0	0,0	0,0	0,0	0,0	-0,0	0,3	0	1	0,7	0,0	1	0,0	0,0
2023-12-31 00:00:00	-6.718.708,5	-3.746.108,1	2.972.600,4	111.094.685,1	4.916.125,7	281.857,5	204.521,0	0,0	3,8	3,7	101,0	884,6	6.030,1	8.638,6	1.236,9	1.236,9	0,0	0,0	3,4	-0,2	0,4	-0,0	-0,0	0,0	0,0	0,4	0,4	0,1	0,1	0,0	0,0	0,0	0,0	0,0	0,3	0	1	0,1	-0,0	0	0,0	-0,0
2024-12-31 00:00:00	-8.880.648,3	-5.110.858,6	3.769.789,8	129.794.268,2	6.324.650,1	311.630,9	209.929,1	0,0	4,1	3,9	105,6	992,1	3.618,2	9.378,3	1.106,2	1.106,2	-1.800,0	0,0	4,5	0,3	0,1	-0,0	-0,0	0,0	0,0	0,4	0,4	0,2	0,1	0,0	0,0	0,0	0,0	-0,0	0,4	0	1	0,3	0,0	1	0,0	0,0

	g_real	def_chile_chg	dep_tcn	copper_r	btu10	i_usa10_mas_EMBI	inf_cl
const	0.416911	2.707989	0.799521	-0.034313	-0.085440	0.022978	-0.276300
L1.g_real	-0.081005	0.126322	-0.444411	0.054891	0.036170	-0.002687	-0.019540
L1.def_chile_chg	0.140598	0.176823	-1.142962	0.225075	0.039275	-0.054476	0.178769
L1.dep_tcn	-0.067954	0.040302	0.207128	-0.020636	-0.014670	0.014171	-0.013442
L1.copper_r	0.171701	-0.406509	0.759525	0.965120	0.024476	-0.001724	0.242074
L1.btu10	0.304462	-0.081682	0.067258	-0.054454	0.914638	0.125146	0.007538
L1.i_usa10_mas_EMBI	-0.326529	0.026332	-1.006222	0.043188	0.021853	0.825124	-0.126405
L1.inf_cl	-0.250536	-0.070410	0.776165	-0.116622	0.044284	0.186217	0.399545

	mediana deuda	media deuda	mediana bp	media bp	mediana gint	media gint	mediana bal	media bal	prob>n.prud	desv.est. deuda	skewness deuda	kurtosis deuda
2025-12-31 00:00:00	40,6	40,6	2,6	2,6	1,2	1,2	1,4	1,4	4,1	2,5	0,0	0,1
2026-12-31 00:00:00	41,2	41,3	2,1	2,1	1,2	1,2	0,9	0,8	12,0	3,2	0,2	0,3
2027-12-31 00:00:00	41,1	41,2	2,6	2,6	1,3	1,3	1,3	1,3	13,2	3,5	0,2	0,2
2028-12-31 00:00:00	41,2	41,3	2,8	2,8	1,3	1,4	1,4	1,4	14,8	3,6	0,2	0,2
2029-12-31 00:00:00	41,1	41,1	3,0	3,1	1,4	1,4	1,6	1,7	14,4	3,7	0,2	0,2