Oil Price Shocks and Fuel Price Pass-through in Portugal

import io
import re
import json
import requests
import numpy as np
import pandas as pd


# ============================================================
# 0. SETTINGS
# ============================================================

START_DATE = "2001-01-01"
COUNTRY = "PT"   # Portugal

WB_OIL_URL = (
    "https://thedocs.worldbank.org/en/doc/"
    "18675f1d1639c7a34d463f59263ba0a2-0050012025/"
    "related/CMO-Historical-Data-Monthly.xlsx"
)

EUROSTAT_BASE = "https://ec.europa.eu/eurostat/api/dissemination/statistics/1.0/data"


# ============================================================
# 1. HELPER: PARSE WORLD BANK MONTHS
# ============================================================

def parse_worldbank_month(x):
    """
    Handles dates like:
    - 2001M01
    - 2001-01
    - Jan-2001
    - Excel datetime
    """
    if pd.isna(x):
        return pd.NaT

    if isinstance(x, (pd.Timestamp,)):
        return x.to_period("M").to_timestamp()

    s = str(x).strip()

    m = re.match(r"^(\d{4})M(\d{1,2})$", s)
    if m:
        year = int(m.group(1))
        month = int(m.group(2))
        return pd.Timestamp(year=year, month=month, day=1)

    dt = pd.to_datetime(s, errors="coerce")
    if pd.isna(dt):
        return pd.NaT

    return dt.to_period("M").to_timestamp()


# ============================================================
# 2. DOWNLOAD WORLD BANK CRUDE OIL DATA
# ============================================================

def download_worldbank_oil():
    r = requests.get(WB_OIL_URL, timeout=60)
    r.raise_for_status()

    xls = pd.ExcelFile(io.BytesIO(r.content))

    monthly_sheet = None
    for sh in xls.sheet_names:
        if "month" in sh.lower():
            monthly_sheet = sh
            break

    if monthly_sheet is None:
        raise ValueError(f"No monthly sheet found. Available sheets: {xls.sheet_names}")

    raw = pd.read_excel(xls, sheet_name=monthly_sheet, header=None)

    # Find the header row containing "Crude oil, average"
    hit = raw.apply(
        lambda row: row.astype(str).str.contains(
            "Crude oil, average", case=False, na=False
        ).any(),
        axis=1
    )

    if not hit.any():
        raise ValueError("Could not find 'Crude oil, average' in the World Bank Excel file.")

    header_row = hit[hit].index[0]

    data = raw.iloc[header_row + 1:].copy()
    data.columns = raw.iloc[header_row]

    date_col = data.columns[0]

    oil_col = None
    for c in data.columns:
        if isinstance(c, str) and "crude oil" in c.lower() and "average" in c.lower():
            oil_col = c
            break

    if oil_col is None:
        raise ValueError("Could not identify the crude oil average column.")

    oil = data[[date_col, oil_col]].copy()
    oil.columns = ["date_raw", "oil_usd"]

    oil["date"] = oil["date_raw"].apply(parse_worldbank_month)
    oil["oil_usd"] = pd.to_numeric(oil["oil_usd"], errors="coerce")

    oil = (
        oil[["date", "oil_usd"]]
        .dropna()
        .drop_duplicates(subset="date")
        .sort_values("date")
        .reset_index(drop=True)
    )

    return oil


# ============================================================
# 3. HELPER: DOWNLOAD EUROSTAT JSON-STAT SERIES
# ============================================================

def eurostat_json_to_series(dataset, params, value_name):
    """
    Downloads one filtered Eurostat series and returns:
    date | value_name

    Works with Eurostat JSON-stat where time is one dimension.
    """
    url = f"{EUROSTAT_BASE}/{dataset}"

    query = {
        "format": "JSON",
        "lang": "en",
        **params
    }

    r = requests.get(url, params=query, timeout=60)
    r.raise_for_status()

    js = r.json()

    if "error" in js:
        raise RuntimeError(js["error"])

    ids = js["id"]
    sizes = js["size"]
    dims = js["dimension"]

    time_dim = None
    for d in ids:
        if d.lower() in ["time", "time_period"]:
            time_dim = d
            break

    if time_dim is None:
        raise ValueError(f"No time dimension found. Dimensions are: {ids}")

    time_pos = ids.index(time_dim)
    time_index = dims[time_dim]["category"]["index"]

    # Sort time labels by Eurostat internal position
    time_labels = sorted(time_index.keys(), key=lambda k: time_index[k])

    values = js.get("value", {})

    def flat_index(coords):
        idx = 0
        multiplier = 1
        for coord, size in zip(reversed(coords), reversed(sizes)):
            idx += coord * multiplier
            multiplier *= size
        return idx

    out = []

    for t in time_labels:
        coords = [0] * len(ids)
        coords[time_pos] = time_index[t]

        idx = flat_index(coords)

        if isinstance(values, dict):
            val = values.get(str(idx), np.nan)
        else:
            val = values[idx] if idx < len(values) else np.nan

        out.append((t, val))

    ser = pd.DataFrame(out, columns=["time_raw", value_name])
    ser["date"] = pd.PeriodIndex(ser["time_raw"], freq="M").to_timestamp()
    ser[value_name] = pd.to_numeric(ser[value_name], errors="coerce")

    ser = (
        ser[["date", value_name]]
        .dropna()
        .drop_duplicates(subset="date")
        .sort_values("date")
        .reset_index(drop=True)
    )

    return ser


# ============================================================
# 4. DOWNLOAD EUROSTAT HICP DATA
# ============================================================

def download_hicp_portugal():
    """
    HICP indices, 2025 = 100.
    TOTAL = total HICP
    FUEL = liquid fuels + fuels and lubricants for personal transport equipment
    """

    common = {
        "freq": "M",
        "unit": "I25",
        "geo": COUNTRY
    }

    total_hicp = eurostat_json_to_series(
        dataset="prc_hicp_minr",
        params={**common, "coicop18": "TOTAL"},
        value_name="hicp_total"
    )

    fuel_hicp = eurostat_json_to_series(
        dataset="prc_hicp_minr",
        params={**common, "coicop18": "FUEL"},
        value_name="fuel_hicp"
    )

    return total_hicp, fuel_hicp


# ============================================================
# 5. DOWNLOAD EUROSTAT EXCHANGE RATE DATA
# ============================================================

def download_eurusd():
    """
    Eurostat gives the monthly average exchange rate as:
    USD per EUR.

    For modelling, we also construct:
    EUR per USD = 1 / USD per EUR.
    """

    eurusd = eurostat_json_to_series(
        dataset="ert_bil_eur_m",
        params={
            "freq": "M",
            "statinfo": "AVG",
            "unit": "NAC",
            "currency": "USD"
        },
        value_name="usd_per_eur"
    )

    eurusd["eur_per_usd"] = 1 / eurusd["usd_per_eur"]

    return eurusd


# ============================================================
# 6. BUILD FINAL DATAFRAME
# ============================================================

oil = download_worldbank_oil()
total_hicp, fuel_hicp = download_hicp_portugal()
fx = download_eurusd()

df = (
    oil
    .merge(fx, on="date", how="inner")
    .merge(fuel_hicp, on="date", how="inner")
    .merge(total_hicp, on="date", how="inner")
    .sort_values("date")
    .reset_index(drop=True)
)

df = df[df["date"] >= pd.Timestamp(START_DATE)].copy()

# Oil price also expressed in EUR, useful for checking transmission,
# but the main four variables remain oil_usd, eur_per_usd, fuel_hicp, hicp_total.
df["oil_eur"] = df["oil_usd"] * df["eur_per_usd"]


# ============================================================
# 7. LOGS AND 12-MONTH LOG DIFFERENCES
# ============================================================

main_vars = [
    "oil_usd",
    "eur_per_usd",
    "fuel_hicp",
    "hicp_total"
]

extra_vars = [
    "usd_per_eur",
    "oil_eur"
]

for col in main_vars + extra_vars:
    df[f"ln_{col}"] = np.log(df[col])

for col in main_vars + extra_vars:
    df[f"yoy_ln_{col}"] = df[f"ln_{col}"].diff(12)

# Monthly log change of oil price: this will be useful for the LP shock later.
df["dln_oil_usd"] = df["ln_oil_usd"].diff()


# ============================================================
# 8. FINAL MODEL DATAFRAME
# ============================================================

df_model = df[
    [
        "date",

        # Raw four variables
        "oil_usd",
        "eur_per_usd",
        "fuel_hicp",
        "hicp_total",

        # Optional checks
        "usd_per_eur",
        "oil_eur",

        # Log-levels
        "ln_oil_usd",
        "ln_eur_per_usd",
        "ln_fuel_hicp",
        "ln_hicp_total",

        # 12-month log differences
        "yoy_ln_oil_usd",
        "yoy_ln_eur_per_usd",
        "yoy_ln_fuel_hicp",
        "yoy_ln_hicp_total",

        # Later LP shock variable
        "dln_oil_usd"
    ]
].copy()

df_model = df_model.reset_index(drop=True)

# ============================================================
# DATA VISUALISATION AND UNIT ROOT TESTS
# Portugal: oil price, exchange rate, fuel HICP, total HICP
# Assumption: df_model already exists from the previous data download step
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings

from statsmodels.tsa.stattools import adfuller, kpss

warnings.filterwarnings("ignore")


# ============================================================
# 1. PREPARE DATA
# ============================================================

df = df_model.copy()
df["date"] = pd.to_datetime(df["date"])
df = df.sort_values("date").reset_index(drop=True)

raw_vars = {
    "oil_usd": "Crude oil price, USD",
    "eur_per_usd": "Exchange rate, EUR per USD",
    "fuel_hicp": "Fuel HICP, Portugal",
    "hicp_total": "Total HICP, Portugal"
}

log_vars = []
yoy_vars = []

for var in raw_vars.keys():
    log_col = f"ln_{var}"
    yoy_col = f"yoy_ln_{var}"

    df[log_col] = np.log(df[var])
    df[yoy_col] = df[log_col].diff(12)

    log_vars.append(log_col)
    yoy_vars.append(yoy_col)


# ============================================================
# 2. PLOT LOG-LEVELS
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 8), sharex=True)
axes = axes.flatten()

for ax, raw_var, log_col in zip(axes, raw_vars.keys(), log_vars):
    ax.plot(df["date"], df[log_col])
    ax.set_title(raw_vars[raw_var])
    ax.set_ylabel("Log level")
    ax.grid(True, alpha=0.3)

plt.suptitle("Log-levels of the four variables", fontsize=14)
plt.tight_layout()
plt.show()


# ============================================================
# 3. PLOT 12-MONTH LOG DIFFERENCES
# ============================================================

fig, axes = plt.subplots(2, 2, figsize=(14, 8), sharex=True)
axes = axes.flatten()

for ax, raw_var, yoy_col in zip(axes, raw_vars.keys(), yoy_vars):
    ax.plot(df["date"], df[yoy_col])
    ax.axhline(0, linewidth=1)
    ax.set_title(raw_vars[raw_var])
    ax.set_ylabel("12-month log difference")
    ax.grid(True, alpha=0.3)

plt.suptitle("Year-on-year log growth rates", fontsize=14)
plt.tight_layout()
plt.show()


# ============================================================
# 4. OPTIONAL COMPARISON PLOT: OIL VS DOMESTIC FUEL PRICES
# ============================================================

fig, ax = plt.subplots(figsize=(12, 5))

ax.plot(
    df["date"],
    df["yoy_ln_oil_usd"],
    label="Crude oil price, USD"
)

ax.plot(
    df["date"],
    df["yoy_ln_fuel_hicp"],
    label="Fuel HICP, Portugal"
)

ax.axhline(0, linewidth=1)
ax.set_title("Oil price growth and domestic fuel price inflation")
ax.set_ylabel("12-month log difference")
ax.grid(True, alpha=0.3)
ax.legend()

plt.tight_layout()
plt.show()


# ============================================================
# 5. UNIT ROOT TEST FUNCTIONS
# ============================================================

def adf_test(series, regression="ct"):
    """
    Augmented Dickey-Fuller test.
    Null hypothesis: the series has a unit root.

    regression='ct' means constant + linear trend.
    This is suitable for log-level macro price series.
    """
    x = series.dropna()

    result = adfuller(
        x,
        regression=regression,
        autolag="AIC"
    )

    return {
        "ADF statistic": result[0],
        "ADF p-value": result[1],
        "ADF lags": result[2],
        "ADF nobs": result[3],
        "ADF 1% critical value": result[4]["1%"],
        "ADF 5% critical value": result[4]["5%"],
        "ADF 10% critical value": result[4]["10%"],
        "ADF decision, 5%": (
            "Reject unit root"
            if result[1] < 0.05
            else "Do not reject unit root"
        )
    }


def kpss_test(series, regression="ct"):
    """
    KPSS test.
    Null hypothesis: the series is trend-stationary.

    regression='ct' means stationarity around constant + linear trend.
    """
    x = series.dropna()

    result = kpss(
        x,
        regression=regression,
        nlags="auto"
    )

    return {
        "KPSS statistic": result[0],
        "KPSS p-value": result[1],
        "KPSS lags": result[2],
        "KPSS 1% critical value": result[3]["1%"],
        "KPSS 5% critical value": result[3]["5%"],
        "KPSS 10% critical value": result[3]["10%"],
        "KPSS decision, 5%": (
            "Reject trend stationarity"
            if result[1] < 0.05
            else "Do not reject trend stationarity"
        )
    }


# ============================================================
# 6. RUN UNIT ROOT TESTS ON LOG-LEVELS
# ============================================================

unit_root_results = []

for raw_var, log_col in zip(raw_vars.keys(), log_vars):
    adf_res = adf_test(df[log_col], regression="ct")
    kpss_res = kpss_test(df[log_col], regression="ct")

    row = {
        "Variable": raw_vars[raw_var],
        "Log variable": log_col,
        **adf_res,
        **kpss_res
    }

    unit_root_results.append(row)

unit_root_table = pd.DataFrame(unit_root_results)


# ============================================================
# 7. PRINT RESULTS
# ============================================================

pd.set_option("display.max_columns", None)
pd.set_option("display.width", 200)

print("\n============================================================")
print("UNIT ROOT TESTS ON LOG-LEVELS")
print("ADF null: unit root")
print("KPSS null: trend stationarity")
print("============================================================\n")

print(unit_root_table.round(4))


# ============================================================
# 8. SHORTER SUMMARY TABLE
# ============================================================

summary_table = unit_root_table[
    [
        "Variable",
        "ADF statistic",
        "ADF p-value",
        "ADF lags",
        "ADF decision, 5%",
        "KPSS statistic",
        "KPSS p-value",
        "KPSS lags",
        "KPSS decision, 5%"
    ]
].copy()

print("\n============================================================")
print("SUMMARY")
print("============================================================\n")

print(summary_table.round(4))


# ============================================================
# 9. OPTIONAL SAVE OUTPUTS
# ============================================================

unit_root_table.to_csv("unit_root_tests_log_levels_full.csv", index=False)
summary_table.to_csv("unit_root_tests_log_levels_summary.csv", index=False)

import pandas as pd
from statsmodels.tsa.api import VAR

# ============================================================
# VAR LAG SELECTION TABLE: AIC AND BIC
# Assumption: df_model already exists
# ============================================================

var_cols = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_fuel_hicp",
    "ln_hicp_total"
]

var_data = (
    df_model[["date"] + var_cols]
    .dropna()
    .copy()
)

var_data["date"] = pd.to_datetime(var_data["date"])
var_data = var_data.sort_values("date").set_index("date")

var_model = VAR(var_data)

max_lags = 12

ic_rows = []

for p in range(1, max_lags + 1):
    res = var_model.fit(p)
    
    ic_rows.append({
        "lag": p,
        "AIC": res.aic,
        "BIC": res.bic
    })

lag_ic_table = pd.DataFrame(ic_rows)

lag_ic_table["AIC_min"] = lag_ic_table["AIC"] == lag_ic_table["AIC"].min()
lag_ic_table["BIC_min"] = lag_ic_table["BIC"] == lag_ic_table["BIC"].min()

print(lag_ic_table.round(4))

selected_lag_aic = lag_ic_table.loc[lag_ic_table["AIC"].idxmin(), "lag"]
selected_lag_bic = lag_ic_table.loc[lag_ic_table["BIC"].idxmin(), "lag"]

print("\nAIC-selected lag:", selected_lag_aic)
print("BIC-selected lag:", selected_lag_bic)

# Estimate baseline VAR using AIC-selected lag
var_results_aic = var_model.fit(selected_lag_aic)

# ============================================================
# CHOLESKY IRFs FROM VAR IN LOG-LEVELS
# Baseline ordering:
# oil price -> exchange rate -> total HICP -> fuel HICP
# Assumption: df_model already exists
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR

# ------------------------------------------------------------
# 1. Baseline Cholesky ordering
# ------------------------------------------------------------

ordering_baseline = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

labels = {
    "ln_oil_usd": "Oil price",
    "ln_eur_per_usd": "EUR per USD",
    "ln_hicp_total": "Total HICP",
    "ln_fuel_hicp": "Fuel HICP"
}

# ------------------------------------------------------------
# 2. Prepare VAR data
# ------------------------------------------------------------

var_data = (
    df_model[["date"] + ordering_baseline]
    .dropna()
    .copy()
)

var_data["date"] = pd.to_datetime(var_data["date"])
var_data = var_data.sort_values("date").set_index("date")

# ------------------------------------------------------------
# 3. Select lag length by AIC and estimate VAR
# ------------------------------------------------------------

var_model = VAR(var_data)

lag_selection = var_model.select_order(maxlags=12)
selected_lag_aic = lag_selection.selected_orders["aic"]

var_results = var_model.fit(selected_lag_aic)

print(f"AIC-selected lag length: {selected_lag_aic}")

# ------------------------------------------------------------
# 4. Compute Cholesky IRFs
# ------------------------------------------------------------

horizon = 24
irf = var_results.irf(horizon)

# Orthogonalized IRFs, i.e. Cholesky IRFs
# Shape: horizons x responses x impulses
orth_irfs = irf.orth_irfs

horizons = np.arange(horizon + 1)

oil_shock_index = ordering_baseline.index("ln_oil_usd")

# ------------------------------------------------------------
# 5. Plot responses to an oil price shock
# ------------------------------------------------------------

fig, axes = plt.subplots(2, 2, figsize=(14, 8), sharex=True)
axes = axes.flatten()

for i, var in enumerate(ordering_baseline):
    response = orth_irfs[:, i, oil_shock_index]

    axes[i].plot(horizons, response)
    axes[i].axhline(0, linewidth=1)
    axes[i].set_title(f"Response of {labels[var]} to oil price shock")
    axes[i].set_xlabel("Months after shock")
    axes[i].set_ylabel("Log response")
    axes[i].grid(True, alpha=0.3)

plt.suptitle("Cholesky impulse responses to an oil price shock", fontsize=14)
plt.tight_layout()
plt.show()

# ------------------------------------------------------------
# 6. Save and print fuel-price response table
# ------------------------------------------------------------

fuel_response = orth_irfs[:, ordering_baseline.index("ln_fuel_hicp"), oil_shock_index]

fuel_irf_table = pd.DataFrame({
    "horizon": horizons,
    "fuel_hicp_response_to_oil_shock": fuel_response
})


fuel_irf_table.to_csv("baseline_cholesky_irf_fuel_response.csv", index=False)

# ============================================================
# BASELINE VS ALTERNATIVE CHOLESKY IRFs
# Four figures/panels: each contains both baseline and alternative IRF
# VAR estimated with 12 lags
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR

# ------------------------------------------------------------
# 1. Settings
# ------------------------------------------------------------

lag_length = 12
horizon = 24

ordering_baseline = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

ordering_alternative = [
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp",
    "ln_oil_usd"
]

response_vars = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

labels = {
    "ln_oil_usd": "Oil price",
    "ln_eur_per_usd": "EUR per USD",
    "ln_hicp_total": "Total HICP",
    "ln_fuel_hicp": "Fuel HICP"
}

# ------------------------------------------------------------
# 2. Prepare data
# ------------------------------------------------------------

df_var = df_model.copy()
df_var["date"] = pd.to_datetime(df_var["date"])
df_var = df_var.sort_values("date")

# ------------------------------------------------------------
# 3. Estimate baseline VAR
# ------------------------------------------------------------

data_baseline = (
    df_var[["date"] + ordering_baseline]
    .dropna()
    .set_index("date")
)

res_baseline = VAR(data_baseline).fit(lag_length)
irf_baseline = res_baseline.irf(horizon).orth_irfs

# ------------------------------------------------------------
# 4. Estimate alternative-ordering VAR
# ------------------------------------------------------------

data_alternative = (
    df_var[["date"] + ordering_alternative]
    .dropna()
    .set_index("date")
)

res_alternative = VAR(data_alternative).fit(lag_length)
irf_alternative = res_alternative.irf(horizon).orth_irfs

# ------------------------------------------------------------
# 5. Extract responses to oil price shock
# ------------------------------------------------------------

oil_shock_baseline_idx = ordering_baseline.index("ln_oil_usd")
oil_shock_alternative_idx = ordering_alternative.index("ln_oil_usd")

h = np.arange(horizon + 1)

# ------------------------------------------------------------
# 6. Plot four panels: baseline and alternative on each
# ------------------------------------------------------------

fig, axes = plt.subplots(2, 2, figsize=(14, 8), sharex=True)
axes = axes.flatten()

for ax, var in zip(axes, response_vars):

    baseline_response = irf_baseline[
        :,
        ordering_baseline.index(var),
        oil_shock_baseline_idx
    ]

    alternative_response = irf_alternative[
        :,
        ordering_alternative.index(var),
        oil_shock_alternative_idx
    ]

    ax.plot(h, baseline_response, label="Baseline ordering")
    ax.plot(h, alternative_response, label="Alternative ordering", linestyle="--")

    ax.axhline(0, linewidth=1)
    ax.set_title(f"Response of {labels[var]} to oil price shock")
    ax.set_xlabel("Months after shock")
    ax.set_ylabel("Log response")
    ax.grid(True, alpha=0.3)
    ax.legend()

plt.suptitle("Cholesky IRFs: baseline vs alternative ordering", fontsize=14)
plt.tight_layout()
plt.show()

# ------------------------------------------------------------
# 7. Optional: save comparison table
# ------------------------------------------------------------

irf_comparison = pd.DataFrame({"horizon": h})

for var in response_vars:
    irf_comparison[f"{var}_baseline"] = irf_baseline[
        :,
        ordering_baseline.index(var),
        oil_shock_baseline_idx
    ]

    irf_comparison[f"{var}_alternative"] = irf_alternative[
        :,
        ordering_alternative.index(var),
        oil_shock_alternative_idx
    ]


irf_comparison.to_csv(
    "cholesky_irfs_baseline_vs_alternative_ordering_lag12.csv",
    index=False
)

# ============================================================
# COMPARE FUEL HICP RESPONSES TO AN OIL PRICE SHOCK
# Specifications:
# 1. Baseline ordering, VAR(12)
# 2. Baseline ordering, VAR(2)
# 3. Alternative ordering, VAR(12)
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR

# ------------------------------------------------------------
# 1. Settings
# ------------------------------------------------------------

horizon = 24

ordering_baseline = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

ordering_alternative = [
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp",
    "ln_oil_usd"
]

# ------------------------------------------------------------
# 2. Prepare common data
# ------------------------------------------------------------

df_var = df_model.copy()
df_var["date"] = pd.to_datetime(df_var["date"])
df_var = df_var.sort_values("date")

all_vars = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

data_common = (
    df_var[["date"] + all_vars]
    .dropna()
    .set_index("date")
)

# ------------------------------------------------------------
# 3. Function to estimate VAR and extract fuel response
# ------------------------------------------------------------

def get_fuel_response_to_oil_shock(data, ordering, lag_length, horizon=24):
    """
    Estimates a VAR with the given ordering and lag length,
    then returns the Cholesky IRF of fuel HICP to an oil price shock.
    """

    data_ordered = data[ordering].copy()

    model = VAR(data_ordered)
    results = model.fit(lag_length)

    irf = results.irf(horizon)
    orth_irfs = irf.orth_irfs

    oil_idx = ordering.index("ln_oil_usd")
    fuel_idx = ordering.index("ln_fuel_hicp")

    fuel_response = orth_irfs[:, fuel_idx, oil_idx]

    return fuel_response, results


# ------------------------------------------------------------
# 4. Estimate the three specifications
# ------------------------------------------------------------

fuel_baseline_lag12, res_baseline_lag12 = get_fuel_response_to_oil_shock(
    data=data_common,
    ordering=ordering_baseline,
    lag_length=12,
    horizon=horizon
)

fuel_baseline_lag2, res_baseline_lag2 = get_fuel_response_to_oil_shock(
    data=data_common,
    ordering=ordering_baseline,
    lag_length=2,
    horizon=horizon
)

fuel_alternative_lag12, res_alternative_lag12 = get_fuel_response_to_oil_shock(
    data=data_common,
    ordering=ordering_alternative,
    lag_length=12,
    horizon=horizon
)

# ------------------------------------------------------------
# 5. Create comparison table
# Log responses are multiplied by 100, so they are approximately percent responses
# ------------------------------------------------------------

h = np.arange(horizon + 1)

fuel_irf_comparison = pd.DataFrame({
    "horizon": h,
    "baseline_ordering_VAR12": 100 * fuel_baseline_lag12,
    "baseline_ordering_VAR2": 100 * fuel_baseline_lag2,
    "alternative_ordering_VAR12": 100 * fuel_alternative_lag12
})



# ------------------------------------------------------------
# 6. Plot comparison
# ------------------------------------------------------------

plt.figure(figsize=(11, 6))

plt.plot(
    fuel_irf_comparison["horizon"],
    fuel_irf_comparison["baseline_ordering_VAR12"],
    label="Baseline ordering, VAR(12)"
)

plt.plot(
    fuel_irf_comparison["horizon"],
    fuel_irf_comparison["baseline_ordering_VAR2"],
    label="Baseline ordering, VAR(2)",
    linestyle="--"
)

plt.plot(
    fuel_irf_comparison["horizon"],
    fuel_irf_comparison["alternative_ordering_VAR12"],
    label="Alternative ordering, VAR(12)",
    linestyle=":"
)

plt.axhline(0, linewidth=1)
plt.title("Fuel HICP response to an oil price shock across specifications")
plt.xlabel("Months after shock")
plt.ylabel("Response of log fuel HICP, percent")
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

# ------------------------------------------------------------
# 7. Optional: compact robustness summary
# ------------------------------------------------------------

summary_rows = []

for col in [
    "baseline_ordering_VAR12",
    "baseline_ordering_VAR2",
    "alternative_ordering_VAR12"
]:
    response = fuel_irf_comparison[col]

    summary_rows.append({
        "specification": col,
        "impact_response_h0": response.iloc[0],
        "response_h6": response.iloc[6],
        "response_h12": response.iloc[12],
        "response_h24": response.iloc[24],
        "peak_response": response.max(),
        "peak_horizon": response.idxmax(),
        "min_response": response.min(),
        "min_horizon": response.idxmin()
    })

fuel_irf_summary = pd.DataFrame(summary_rows)

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

# ------------------------------------------------------------
# 1. Settings
# ------------------------------------------------------------

H = 24
P = 12

shock_var = "dln_oil_usd"
response_var = "ln_fuel_hicp"

# Contemporaneous controls
# Do NOT include ln_fuel_hicp_t, because at h = 0 it is mechanically related to the dependent variable.
# Do NOT include ln_oil_usd_t, because dln_oil_usd is the shock variable.
contemporaneous_controls = [
    "ln_eur_per_usd",
    "ln_hicp_total"
]

# Lagged controls
# Include fuel HICP lags to control for persistence.
# I also include lagged exchange rate and total HICP.
lagged_controls = [
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

# ------------------------------------------------------------
# 2. Prepare data
# ------------------------------------------------------------

df_lp = df_model.copy()
df_lp["date"] = pd.to_datetime(df_lp["date"])
df_lp = df_lp.sort_values("date").reset_index(drop=True)

# Monthly oil price change
df_lp["dln_oil_usd"] = df_lp["ln_oil_usd"].diff()

# Lagged fuel price used as pre-shock baseline
df_lp["ln_fuel_hicp_lag1"] = df_lp["ln_fuel_hicp"].shift(1)

# ------------------------------------------------------------
# 3. Add controls
# ------------------------------------------------------------

x_vars = [shock_var]

# Contemporaneous controls
for var in contemporaneous_controls:
    x_vars.append(var)

# Lagged controls: 1 to 12
for lag in range(1, P + 1):
    for var in lagged_controls:
        lag_name = f"{var}_lag{lag}"
        df_lp[lag_name] = df_lp[var].shift(lag)
        x_vars.append(lag_name)

# ------------------------------------------------------------
# 4. Estimate local projections
# ------------------------------------------------------------

lp_results = []

for h in range(H + 1):

    # Cumulative fuel-price response:
    # ln(Fuel_{t+h}) - ln(Fuel_{t-1})
    y_name = f"fuel_response_h{h}"
    df_lp[y_name] = df_lp[response_var].shift(-h) - df_lp["ln_fuel_hicp_lag1"]

    reg_data = df_lp[[y_name] + x_vars].dropna().copy()

    y = reg_data[y_name]
    X = sm.add_constant(reg_data[x_vars])

    model = sm.OLS(y, X)

    # HAC standard errors for overlapping horizons
    res = model.fit(cov_type="HAC", cov_kwds={"maxlags": h + 1})

    beta = res.params[shock_var]
    se = res.bse[shock_var]

    lp_results.append({
        "horizon": h,
        "beta": beta,
        "se": se,
        "ci_lower": beta - 1.96 * se,
        "ci_upper": beta + 1.96 * se,
        "nobs": int(res.nobs)
    })

lp_irf = pd.DataFrame(lp_results)

# ------------------------------------------------------------
# 5. Scale LP response to one-standard-deviation oil price shock
# ------------------------------------------------------------

shock_std = df_lp[shock_var].dropna().std()

lp_irf["beta_1sd"] = lp_irf["beta"] * shock_std
lp_irf["ci_lower_1sd"] = lp_irf["ci_lower"] * shock_std
lp_irf["ci_upper_1sd"] = lp_irf["ci_upper"] * shock_std

# Convert log response to approximate percent response
lp_irf["beta_percent"] = 100 * lp_irf["beta_1sd"]
lp_irf["ci_lower_percent"] = 100 * lp_irf["ci_lower_1sd"]
lp_irf["ci_upper_percent"] = 100 * lp_irf["ci_upper_1sd"]

# ------------------------------------------------------------
# 6. Plot local projection IRF
# ------------------------------------------------------------

plt.figure(figsize=(10, 5))

plt.plot(
    lp_irf["horizon"],
    lp_irf["beta_percent"],
    label="LP response to one-s.d. oil price shock"
)

plt.fill_between(
    lp_irf["horizon"],
    lp_irf["ci_lower_percent"],
    lp_irf["ci_upper_percent"],
    alpha=0.2,
    label="95% confidence interval"
)

plt.axhline(0, linewidth=1)
plt.title("Local projection response of fuel HICP to an oil price shock")
plt.xlabel("Months after shock")
plt.ylabel("Cumulative response of log fuel HICP, percent")
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

# ============================================================
# ALTERNATIVE OIL PRICE SHOCK:
# unpredictable component of monthly oil price changes
# ============================================================

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

# ------------------------------------------------------------
# 1. Settings
# ------------------------------------------------------------

P = 12

df_shock = df_model.copy()
df_shock["date"] = pd.to_datetime(df_shock["date"])
df_shock = df_shock.sort_values("date").reset_index(drop=True)

# Monthly oil price change
df_shock["dln_oil_usd"] = df_shock["ln_oil_usd"].diff()

# Other indicators in logs
other_indicators = [
    "ln_eur_per_usd",
    "ln_fuel_hicp",
    "ln_hicp_total"
]

# ------------------------------------------------------------
# 2. Construct lagged explanatory variables
# ------------------------------------------------------------

x_vars = []

# Own lagged oil price changes
for lag in range(1, P + 1):
    lag_name = f"dln_oil_usd_lag{lag}"
    df_shock[lag_name] = df_shock["dln_oil_usd"].shift(lag)
    x_vars.append(lag_name)

# Lagged values of the other indicators
for lag in range(1, P + 1):
    for var in other_indicators:
        lag_name = f"{var}_lag{lag}"
        df_shock[lag_name] = df_shock[var].shift(lag)
        x_vars.append(lag_name)

# ------------------------------------------------------------
# 3. Estimate predictability regression
# ------------------------------------------------------------

reg_data = df_shock[["date", "dln_oil_usd"] + x_vars].dropna().copy()

y = reg_data["dln_oil_usd"]
X = sm.add_constant(reg_data[x_vars])

oil_predictability_model = sm.OLS(y, X).fit()

# Fitted predictable component and residual shock
reg_data["dln_oil_predicted"] = oil_predictability_model.fittedvalues
reg_data["oil_shock_residual"] = oil_predictability_model.resid

# Merge residual shock back into main dataframe
df_shock = df_shock.merge(
    reg_data[["date", "dln_oil_predicted", "oil_shock_residual"]],
    on="date",
    how="left"
)

# ------------------------------------------------------------
# 4. Compare actual oil price change and residual shock
# ------------------------------------------------------------

comparison_data = df_shock[
    ["date", "dln_oil_usd", "dln_oil_predicted", "oil_shock_residual"]
].dropna().copy()

r2 = oil_predictability_model.rsquared
adj_r2 = oil_predictability_model.rsquared_adj

corr_actual_resid = comparison_data["dln_oil_usd"].corr(
    comparison_data["oil_shock_residual"]
)

std_actual = comparison_data["dln_oil_usd"].std()
std_resid = comparison_data["oil_shock_residual"].std()
std_ratio = std_resid / std_actual

summary = pd.DataFrame({
    "Statistic": [
        "R-squared",
        "Adjusted R-squared",
        "Correlation: actual oil change vs residual shock",
        "Std. dev. of actual oil change",
        "Std. dev. of residual shock",
        "Residual std. dev. / actual std. dev."
    ],
    "Value": [
        r2,
        adj_r2,
        corr_actual_resid,
        std_actual,
        std_resid,
        std_ratio
    ]
})

print(summary.round(4))

# ------------------------------------------------------------
# 5. Plot actual oil price change and residual shock
# Multiplied by 100, so values are approximately monthly percent changes
# ------------------------------------------------------------

plt.figure(figsize=(12, 5))

plt.plot(
    comparison_data["date"],
    100 * comparison_data["dln_oil_usd"],
    label="Actual monthly oil price change"
)

plt.plot(
    comparison_data["date"],
    100 * comparison_data["oil_shock_residual"],
    label="Residual oil price shock",
    linestyle="--"
)

plt.axhline(0, linewidth=1)
plt.title("Actual monthly oil price change and unpredictable oil price shock")
plt.xlabel("Date")
plt.ylabel("Monthly log change, percent")
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

# ============================================================
# COMPARE LP RESPONSES (TWO SHOCK DEFINITIONS) WITH VAR IRF
# Output: one figure only
# Assumption: df_model already exists
#
# Important:
#   - LP responses are cumulative fuel-price responses:
#       ln_fuel_hicp_{t+h} - ln_fuel_hicp_{t-1}
#   - LP coefficients are scaled to one-standard-deviation shocks,
#     so they are comparable with Cholesky VAR IRFs.
# ============================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.api import VAR

# ------------------------------------------------------------
# 1. Settings
# ------------------------------------------------------------

H = 24
P = 12

ordering_baseline = [
    "ln_oil_usd",
    "ln_eur_per_usd",
    "ln_hicp_total",
    "ln_fuel_hicp"
]

# ------------------------------------------------------------
# 2. Prepare data
# ------------------------------------------------------------

df = df_model.copy()
df["date"] = pd.to_datetime(df["date"])
df = df.sort_values("date").reset_index(drop=True)

# Monthly oil-price change
df["dln_oil_usd"] = df["ln_oil_usd"].diff()

# Pre-shock fuel-price baseline
df["ln_fuel_hicp_lag1"] = df["ln_fuel_hicp"].shift(1)

# ------------------------------------------------------------
# 3. Construct alternative shock:
#    residual from regression of oil-price change on its own lags
#    and lagged values of the other indicators
# ------------------------------------------------------------

other_indicators = [
    "ln_eur_per_usd",
    "ln_fuel_hicp",
    "ln_hicp_total"
]

shock_regressors = []

# Own lagged oil-price changes
for lag in range(1, P + 1):
    col = f"dln_oil_usd_lag{lag}"
    df[col] = df["dln_oil_usd"].shift(lag)
    shock_regressors.append(col)

# Lagged values of the other indicators
for lag in range(1, P + 1):
    for var in other_indicators:
        col = f"{var}_shock_lag{lag}"
        df[col] = df[var].shift(lag)
        shock_regressors.append(col)

shock_data = df[["date", "dln_oil_usd"] + shock_regressors].dropna().copy()

y_shock = shock_data["dln_oil_usd"]
X_shock = sm.add_constant(shock_data[shock_regressors])

shock_model = sm.OLS(y_shock, X_shock).fit()

shock_data["oil_shock_residual"] = shock_model.resid

df = df.merge(
    shock_data[["date", "oil_shock_residual"]],
    on="date",
    how="left"
)

# ------------------------------------------------------------
# 4. Function to estimate local projections
# ------------------------------------------------------------

def estimate_lp(df_in, shock_var, horizons=24, lags=12):
    df_lp = df_in.copy()

    # Contemporaneous controls
    # Do NOT include ln_fuel_hicp_t, because at h = 0 it is mechanically related to the dependent variable.
    # Do NOT include ln_oil_usd_t, because dln_oil_usd / residual oil shock is the shock variable.
    contemporaneous_controls = [
        "ln_eur_per_usd",
        "ln_hicp_total"
    ]

    # Lagged controls
    lagged_controls = [
        "ln_eur_per_usd",
        "ln_hicp_total",
        "ln_fuel_hicp"
    ]

    x_vars = [shock_var]

    # Contemporaneous controls
    for var in contemporaneous_controls:
        x_vars.append(var)

    # Lagged controls
    for lag in range(1, lags + 1):
        for var in lagged_controls:
            col = f"{var}_LP_lag{lag}"
            df_lp[col] = df_lp[var].shift(lag)
            x_vars.append(col)

    results = []

    for h in range(horizons + 1):

        # Cumulative fuel-price response:
        # ln(Fuel_{t+h}) - ln(Fuel_{t-1})
        y_col = f"fuel_response_h{h}"
        df_lp[y_col] = df_lp["ln_fuel_hicp"].shift(-h) - df_lp["ln_fuel_hicp_lag1"]

        reg_data = df_lp[[y_col] + x_vars].dropna().copy()

        y = reg_data[y_col]
        X = sm.add_constant(reg_data[x_vars])

        res = sm.OLS(y, X).fit(
            cov_type="HAC",
            cov_kwds={"maxlags": h + 1}
        )

        beta = res.params[shock_var]
        se = res.bse[shock_var]

        results.append({
            "horizon": h,
            "beta": beta,
            "se": se,
            "lower": beta - 1.96 * se,
            "upper": beta + 1.96 * se
        })

    out = pd.DataFrame(results)

    # Scale LP response to one-standard-deviation shock
    shock_std = df_lp[shock_var].dropna().std()

    out["resp_1sd"] = out["beta"] * shock_std
    out["lower_1sd"] = out["lower"] * shock_std
    out["upper_1sd"] = out["upper"] * shock_std

    return out, shock_std

# LP with raw monthly oil-price change
lp_raw, shock_std_raw = estimate_lp(
    df,
    shock_var="dln_oil_usd",
    horizons=H,
    lags=P
)

# LP with residual oil shock
lp_resid, shock_std_resid = estimate_lp(
    df,
    shock_var="oil_shock_residual",
    horizons=H,
    lags=P
)

# ------------------------------------------------------------
# 5. Estimate baseline VAR and obtain Cholesky IRF
# ------------------------------------------------------------

var_data = (
    df[["date"] + ordering_baseline]
    .dropna()
    .copy()
)

var_data = var_data.sort_values("date").set_index("date")

var_model = VAR(var_data)
selected_lag = var_model.select_order(maxlags=12).selected_orders["aic"]
var_res = var_model.fit(selected_lag)

irf = var_res.irf(H).orth_irfs

oil_idx = ordering_baseline.index("ln_oil_usd")
fuel_idx = ordering_baseline.index("ln_fuel_hicp")

var_irf = irf[:, fuel_idx, oil_idx]

# ------------------------------------------------------------
# 6. Plot comparison
# ------------------------------------------------------------

h = np.arange(H + 1)

plt.figure(figsize=(11, 6))

# VAR IRF
plt.plot(
    h,
    100 * var_irf,
    label=f"VAR Cholesky IRF, AIC lag = {selected_lag}",
    linewidth=2
)

# LP with raw oil-price change
plt.plot(
    h,
    100 * lp_raw["resp_1sd"],
    label="LP: monthly oil-price change shock",
    linestyle="--"
)

plt.fill_between(
    h,
    100 * lp_raw["lower_1sd"],
    100 * lp_raw["upper_1sd"],
    alpha=0.18
)

# LP with residual oil shock
plt.plot(
    h,
    100 * lp_resid["resp_1sd"],
    label="LP: residual oil shock",
    linestyle=":"
)

plt.fill_between(
    h,
    100 * lp_resid["lower_1sd"],
    100 * lp_resid["upper_1sd"],
    alpha=0.18
)

plt.axhline(0, linewidth=1)
plt.xlabel("Months after shock")
plt.ylabel("Cumulative response of log fuel HICP, percent")
plt.title("Fuel-price response to one-s.d. oil shocks: LP vs VAR")
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

print("Standard deviation of raw oil-price shock:", round(shock_std_raw, 6))
print("Standard deviation of residual oil shock:", round(shock_std_resid, 6))