Economía y econometría

import numpy as np, pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

raw = sm.datasets.macrodata.load_pandas().data
raw.index = pd.PeriodIndex.from_fields(year=raw["year"].astype(int), quarter=raw["quarter"].astype(int), freq="Q")
print(sm.datasets.macrodata.NOTE.split("Number of Variables")[0].strip()[:400])
raw[["realgdp", "unemp", "infl", "tbilrate"]].tail()

df = pd.DataFrame({
    "gdp_growth": 400 * np.log(raw["realgdp"]).diff(),     # crecimiento anualizado en % (diferencia logarítmica)
    "unemp": raw["unemp"],
    "d_unemp": raw["unemp"].diff(),                        # cambio de la tasa de desempleo, en puntos
    "infl": raw["infl"],
}).dropna()

fig, ax = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
df["gdp_growth"].plot(ax=ax[0], color="#4f46e5"); ax[0].axhline(0, color="grey", lw=.6); ax[0].set_ylabel("crecimiento del PIB, % anual")
df["unemp"].plot(ax=ax[1], color="#ef4444"); ax[1].set_ylabel("desempleo, %")
plt.tight_layout(); plt.show()
df.describe().round(2)

X = sm.add_constant(df["gdp_growth"])
okun = sm.OLS(df["d_unemp"], X).fit(cov_type="HAC", cov_kwds={"maxlags": 4})   # robusto a la autocorrelación
print(okun.summary().tables[1])
alpha, beta = okun.params
print(f"\nR² = {okun.rsquared:.2f}. Cada punto adicional de crecimiento anualizado reduce el desempleo en "
      f"{-beta:.3f} puntos por trimestre (~{-4 * beta:.2f} al año).")
print(f"Crecimiento necesario solo para mantener estable el desempleo: {-alpha / beta:.1f}% al año.")

fig, ax = plt.subplots(1, 2, figsize=(11, 3.8))
ax[0].scatter(df["gdp_growth"], df["d_unemp"], s=10, alpha=.6)
xs = np.linspace(df["gdp_growth"].min(), df["gdp_growth"].max(), 50)
ax[0].plot(xs, alpha + beta * xs, color="#ef4444"); ax[0].set_xlabel("crecimiento del PIB, % anual"); ax[0].set_ylabel("cambio del desempleo")
ax[0].set_title("Ley de Okun, 1959-2009")
resid = okun.resid
sm.graphics.tsa.plot_acf(resid, lags=12, ax=ax[1]); ax[1].set_title("autocorrelación de los residuos")
plt.tight_layout(); plt.show()

rows = []
for start in range(1960, 2010, 10):
    d = df[(df.index.year >= start) & (df.index.year < start + 10)]
    fit = sm.OLS(d["d_unemp"], sm.add_constant(d["gdp_growth"])).fit()
    lo, hi = fit.conf_int().loc["gdp_growth"]
    rows.append(dict(decade=f"{start}s", beta=fit.params["gdp_growth"], ci_low=lo, ci_high=hi, quarters=len(d)))
pd.DataFrame(rows).set_index("decade").round(3)

import warnings
from statsmodels.tsa.arima.model import ARIMA
warnings.filterwarnings("ignore")

y = df["gdp_growth"].astype(float).reset_index(drop=True)
dates = df.index
start = int(np.where(dates.year == 1990)[0][0])
preds = {"ARIMA(1,0,1)": [], "media histórica": [], "último valor": []}
actual = []
for t in range(start, len(y)):
    past = y.iloc[:t]
    preds["ARIMA(1,0,1)"].append(ARIMA(past, order=(1, 0, 1)).fit().forecast(1).iloc[0])
    preds["media histórica"].append(past.mean())
    preds["último valor"].append(past.iloc[-1])
    actual.append(y.iloc[t])
actual = np.array(actual)
scores = pd.DataFrame({k: {"MAE": np.mean(np.abs(np.array(v) - actual)), "RMSE": np.sqrt(np.mean((np.array(v) - actual) ** 2))}
                       for k, v in preds.items()}).T.round(3)
print(f"{len(actual)} previsiones a un trimestre vista, 1990-2009")
scores

fit = ARIMA(y, order=(1, 0, 1)).fit()
fc = fit.get_forecast(8)
ci = fc.conf_int(alpha=0.2)
hist = y.iloc[-24:]
plt.figure(figsize=(10, 3.6))
plt.plot(range(len(y) - 24, len(y)), hist, label="crecimiento del PIB")
fx = range(len(y), len(y) + 8)
plt.plot(fx, fc.predicted_mean, color="#ef4444", label="previsión ARIMA")
plt.fill_between(fx, ci.iloc[:, 0], ci.iloc[:, 1], color="#ef4444", alpha=.2, label="intervalo del 80%")
plt.axhline(0, color="grey", lw=.6); plt.legend(); plt.title("Previsión a dos años tras la recesión de 2008-09"); plt.tight_layout(); plt.show()