Bayesian Modelling

# Generate a sequence of binary random variables
from scipy.stats import bernoulli, beta, uniform
p = 0.3
r = bernoulli.rvs(p, size=400)
r[:10]

import matplotlib.pyplot as plt
import numpy as np

# prior distribution: uniform prior
prior_alpha = 1
prior_beta = 1

# Now we plot the distribution as we add more data points
fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(17,17))
N = np.linspace(0, len(r), 9)
for i in range(9):
  n = N[i]
  i_x = int(i/3)
  i_y = i % 3 
  r_trunc = r[:int(n)]  
  p_grid = np.linspace(0, 1, 1000)
  post_alpha = prior_alpha + np.sum(r_trunc)
  post_beta = prior_beta + len(r_trunc)- np.sum(r_trunc)
  mean = beta.mean(post_alpha, post_beta)
  std = beta.std(post_alpha, post_beta)
  ax[i_x][i_y].plot(p_grid, beta.pdf(p_grid, post_alpha, post_beta))
  ax[i_x][i_y].set_title("N = " + str(int(n))+ ", est = " + str(np.round(mean, 3)) + " +/- " + str(np.round(std, 3)))

import matplotlib.pyplot as plt
import numpy as np

# prior distribution: non - informative prior (approximation)
prior_alpha = 0.000001
prior_beta = 0.000001

# Now we plot the distribution as we add more data points
fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(17,17))
N = np.linspace(0, len(r), 9)
for i in range(9):
  n = N[i]
  i_x = int(i/3)
  i_y = i % 3 
  r_trunc = r[:int(n)]  
  post_alpha = prior_alpha + np.sum(r_trunc)
  post_beta = prior_beta + len(r_trunc)- np.sum(r_trunc)
  mean = beta.mean(post_alpha, post_beta)
  std = beta.std(post_alpha, post_beta)
  ax[i_x][i_y].plot(p_grid, beta.pdf(p_grid, post_alpha, post_beta))
  ax[i_x][i_y].set_title("N = " + str(int(n))+ ", est = " + str(np.round(mean, 3)) + " +/- " + str(np.round(std, 3)))

import matplotlib.pyplot as plt
import numpy as np

# prior distribution: wrong confident prior
prior_alpha = 30
prior_beta = 30

# Now we plot the distribution as we add more data points
fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(17,17))
N = np.linspace(0, len(r), 9)
for i in range(9):
  n = N[i]
  i_x = int(i/3)
  i_y = i % 3 
  r_trunc = r[:int(n)]  
  post_alpha = prior_alpha + np.sum(r_trunc)
  post_beta = prior_beta + len(r_trunc)- np.sum(r_trunc)
  mean = beta.mean(post_alpha, post_beta)
  std = beta.std(post_alpha, post_beta)
  ax[i_x][i_y].plot(p_grid, beta.pdf(p_grid, post_alpha, post_beta))
  ax[i_x][i_y].set_title("N = " + str(int(n))+ ", est = " + str(np.round(mean, 3)) + " +/- " + str(np.round(std, 3)))

import matplotlib.pyplot as plt
import numpy as np

# prior distribution: correct confident prior
prior_alpha = 18
prior_beta = 42

# Now we plot the distribution as we add more data points
fig, ax = plt.subplots(nrows=3, ncols=3, figsize=(17,17))
N = np.linspace(0, len(r), 9)
for i in range(9):
  n = N[i]
  i_x = int(i/3)
  i_y = i % 3 
  r_trunc = r[:int(n)]  
  post_alpha = prior_alpha + np.sum(r_trunc)
  post_beta = prior_beta + len(r_trunc)- np.sum(r_trunc)
  mean = beta.mean(post_alpha, post_beta)
  std = beta.std(post_alpha, post_beta)
  ax[i_x][i_y].plot(p_grid, beta.pdf(p_grid, post_alpha, post_beta))
  ax[i_x][i_y].set_title("N = " + str(int(n))+ ", est = " + str(np.round(mean, 3)) + " +/- " + str(np.round(std, 3)))

from scipy.stats import norm
import numpy as np

class TheGoodAndBadDataModel():
  def __init__(self, prior_p_good, prior_mean, prior_std_good, prior_std_bad):
    self.p_good = prior_p_good
    self.mean = prior_mean
    self.std_good = prior_std_good
    self.std_bad = prior_std_bad
    self.loglik_history = []  # store likelihood values

  def predict(self, X):
    p_x_bad_pbad = (1 - self.p_good) * norm.pdf(X, loc=self.mean, scale=self.std_bad)
    p_x_good_pgood = self.p_good * norm.pdf(X, loc=self.mean, scale=self.std_good)
    p_bad_x = p_x_bad_pbad / (p_x_good_pgood + p_x_bad_pbad)
    return p_bad_x

  def compute_loglik(self, X):
    mixture_pdf = (
        self.p_good * norm.pdf(X, loc=self.mean, scale=self.std_good) +
        (1 - self.p_good) * norm.pdf(X, loc=self.mean, scale=self.std_bad)
    )
    return np.sum(np.log(mixture_pdf + 1e-12))  # add epsilon to avoid log(0)

  def learn(self, X, max_iter=1000, tolerance=1e-5, print_error=False, track_likelihood=True):
    iter = 0
    while True:
      iter += 1
      # E-step
      p_bad_s = self.predict(X)
      p_good_s = 1 - p_bad_s

      # M-step
      p_good_sp1 = np.mean(p_good_s)
      std_good_sp1 = np.sqrt(np.sum(p_good_s * (X - self.mean)**2) / (len(X) * p_good_sp1))
      std_bad_sp1 = np.sqrt(np.sum(p_bad_s * (X - self.mean)**2) / (len(X) * (1 - p_good_sp1)))

      # compute change (for stopping)
      error = np.sqrt(((p_good_sp1 - self.p_good)/self.p_good)**2
                      + ((std_good_sp1 - self.std_good)/self.std_good)**2
                      + ((std_bad_sp1 - self.std_bad)/self.std_bad)**2)

      # update parameters
      self.p_good = p_good_sp1
      self.std_good = std_good_sp1
      self.std_bad = std_bad_sp1

      # track log-likelihood
      if track_likelihood:
        ll = self.compute_loglik(X)
        self.loglik_history.append(ll)
        if print_error:
          print(f"Iter {iter}: error={error:.6f}, loglik={ll:.6f}")

      # stopping condition
      if (error < tolerance or iter >= max_iter):
        break

import yfinance as yf
import numpy as np

# Define ticker
bbva_tkr = yf.Ticker("BBVA.MC")

# Get 10 years of daily adjusted close prices
end_date = "2025-07-31"
start_date = "2015-07-31"  # 10 years earlier
data = bbva_tkr.history(start=start_date, end=end_date, interval="1d")["Close"]

data = bbva_tkr.history(period="10y", interval="1d")["Close"]

# Calculate daily returns
bbva_ret = data.pct_change().dropna()

# Print mean and standard deviation of returns
print("Mean return:", np.mean(bbva_ret))
print("Std deviation:", np.std(bbva_ret))

# We learn the model over the historical data. As a prior we assume there are only 10% anomalies in the dataset
# The result is relatively robust to the choice of prior
gbdm = TheGoodAndBadDataModel(0.99, np.mean(bbva_ret), np.std(bbva_ret), 2*np.std(bbva_ret))
gbdm.learn(bbva_ret)

# We have a look at the results: according to the model, there are 15% anomalies, with roughly a 3x standard deviation
# This means the model detects that the distribution of returns is not accurately described by a single Gaussian
print(gbdm.p_good, gbdm.mean, gbdm.std_good, gbdm.std_bad)

# Plot the log-likelihood history
import matplotlib.pyplot as plt
plt.plot(gbdm.loglik_history, marker="o")
plt.xlabel("Iteration")
plt.ylabel("Log-likelihood")
plt.title("EM Log-likelihood evolution")
plt.show()

import matplotlib.pyplot as plt

bbva_ret.hist(bins=50, figsize=(8,5))
plt.xlabel("Daily Return")
plt.ylabel("Frequency")
plt.title("Histogram of BBVA Daily Returns (10y)")
plt.show()

# Let us flag anomalies as those with 50% of more probability of belonging to the "bad data" mixture component
anomalies = (gbdm.predict(bbva_ret) > 0.5)
# Add anomaly flag to DataFrame
pd_bbva_ret = bbva_ret.reset_index()
pd_bbva_ret.columns = ["Date", "Returns"]
pd_bbva_ret["anomaly"] = anomalies

# Plot histogram of returns, segmented by anomaly flag
pd_bbva_ret[pd_bbva_ret["anomaly"] == True]["Returns"].hist(
    bins=100, alpha=0.7, color="red", label="Anomalies")
pd_bbva_ret[pd_bbva_ret["anomaly"] == False]["Returns"].hist(
    bins=50, alpha=0.7, color="blue", label="Normal")

plt.xlabel("Daily Return")
plt.ylabel("Frequency")
plt.title("Histogram of BBVA Daily Returns")
plt.legend()
plt.show()

# In terms of time-series, we flat the anomalies in the following plot
colors = {False: "blue", True: "red"}
pd_bbva_ret.reset_index().plot.scatter(x = "Date", y = "Returns", c = pd_bbva_ret["anomaly"].map(colors).values, title = "Time Series of BBVA Daily Returns")

# An interesting question is whether these anomalies tend to cluster. Looking at the auto-correlation it hints this is a possibility
# This means it could make sense to analyse these anomalies as regime changes, i.e. use a hidden markov model
from pandas.plotting import autocorrelation_plot
ax = autocorrelation_plot(pd_bbva_ret["anomaly"])
ax.set_xlim([0, 100])

# Check that relevant periods of financial stress have been correctly classified
pd_bbva_ret["Date"] = pd.to_datetime(pd_bbva_ret["Date"]).dt.date

# Define event dates
event_dates = [
    pd.to_datetime("2020-03-16").date(),  # COVID lockdown Spain
    pd.to_datetime("2016-06-24").date(),  # Brexit referendum (first trading day)
    pd.to_datetime("2016-11-09").date()   # Trump election (first trading day after)
]

# Filter rows that match event dates
events_df = pd_bbva_ret[pd_bbva_ret["Date"].isin(event_dates)][["Date", "Returns", "anomaly"]]

print(events_df)

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline

def simulate_local_level(T=200, sigma_v2=1.0, sigma_w2=0.05, seed=42):
    rng = np.random.default_rng(seed)
    y = np.zeros(T)
    x = np.zeros(T)
    y[0] = rng.normal(0.0, np.sqrt(sigma_w2))
    x[0] = y[0] + rng.normal(0.0, np.sqrt(sigma_v2))
    for t in range(1, T):
        y[t] = y[t-1] + rng.normal(0.0, np.sqrt(sigma_w2))
        x[t] = y[t] + rng.normal(0.0, np.sqrt(sigma_v2))
    return x, y



def kf_rts_identity_cross(x, sigma_v2, sigma_w2, m0=0.0, P0=1e6):
    T = len(x)
    m_pred = np.zeros(T)
    P_pred = np.zeros(T)
    m_filt = np.zeros(T)
    P_filt = np.zeros(T)
    K = np.zeros(T)
    innov = np.zeros(T)
    S = np.zeros(T)

    m_prev, P_prev = m0, P0
    loglik = 0.0
    for t in range(T):
        # predict
        m_pred[t] = m_prev
        P_pred[t] = P_prev + sigma_w2
        # update
        innov[t] = x[t] - m_pred[t]
        S[t] = P_pred[t] + sigma_v2
        K[t] = P_pred[t] / S[t]
        m_filt[t] = m_pred[t] + K[t] * innov[t]
        P_filt[t] = (1 - K[t]) * P_pred[t]
        # log-likelihood
        loglik += -0.5 * (np.log(2*np.pi*S[t]) + innov[t]**2 / S[t])
        m_prev, P_prev = m_filt[t], P_filt[t]

    # RTS smoother
    m_smooth = np.zeros(T)
    P_smooth = np.zeros(T)
    J = np.zeros(T-1)
    m_smooth[-1] = m_filt[-1]
    P_smooth[-1] = P_filt[-1]
    for t in range(T-2, -1, -1):
        J[t] = P_filt[t] / P_pred[t+1]
        m_smooth[t] = m_filt[t] + J[t] * (m_smooth[t+1] - m_pred[t+1])
        P_smooth[t] = P_filt[t] + J[t]**2 * (P_smooth[t+1] - P_pred[t+1])

    # lag-one smoothed covariance via identity: P_{t-1,t|T} = J_{t-1} P_{t|T}
    P_cross = np.zeros(T-1)
    for t in range(1, T):
        P_cross[t-1] = J[t-1] * P_smooth[t]

    return {
        "m_pred": m_pred, "P_pred": P_pred,
        "m_filt": m_filt, "P_filt": P_filt,
        "m_smooth": m_smooth, "P_smooth": P_smooth,
        "innov": innov, "S": S, "K": K, "J": J,
        "P_cross": P_cross, "loglik": loglik
    }



def em_local_level(x, sigma_v2_init=2.0, sigma_w2_init=0.2, m0=0.0, P0=1e6, max_iter=500, tol=1e-8):
    sigma_v2, sigma_w2 = float(sigma_v2_init), float(sigma_w2_init)
    ll_hist = []
    for it in range(max_iter):
        out = kf_rts_identity_cross(x, sigma_v2, sigma_w2, m0, P0)
        mu, P, Pc = out["m_smooth"], out["P_smooth"], out["P_cross"]
        ll_hist.append(out["loglik"])

        # E-step expectations
        E1 = (x - mu)**2 + P
        diff_mu = mu[1:] - mu[:-1]
        E2 = diff_mu**2 + P[1:] + P[:-1] - 2.0 * Pc

        # M-step (MLE scaling)
        sigma_v2_new = max(np.mean(E1), 1e-12)
        sigma_w2_new = max(np.mean(E2), 1e-12)

        # convergence
        rel = max(abs(sigma_v2_new - sigma_v2) / (sigma_v2 + 1e-12),
                  abs(sigma_w2_new - sigma_w2) / (sigma_w2 + 1e-12))
        sigma_v2, sigma_w2 = sigma_v2_new, sigma_w2_new
        if rel < tol:
            ll_hist.append(kf_rts_identity_cross(x, sigma_v2, sigma_w2, m0, P0)["loglik"])
            break
    return {"sigma_v2": sigma_v2, "sigma_w2": sigma_w2, "ll_history": np.array(ll_hist)}

# --- Configuration ---
T = 200
sigma_v2_true = 1.0
sigma_w2_true = 0.05
seed = 42

# EM seeds
sigma_v2_init = 2.0
sigma_w2_init = 0.2

# --- Simulate ---
x, y = simulate_local_level(T, sigma_v2_true, sigma_w2_true, seed)

# --- Run EM ---
em = em_local_level(x, sigma_v2_init=sigma_v2_init, sigma_w2_init=sigma_w2_init, max_iter=1000, tol=1e-5)
sigma_v2_hat, sigma_w2_hat = em["sigma_v2"], em["sigma_w2"]
ll = em["ll_history"]

# --- Smooth with estimated params ---
out = kf_rts_identity_cross(x, sigma_v2_hat, sigma_w2_hat)

# --- Compare parameters ---
df = pd.DataFrame({
    "Parameter": [r"$\sigma_v^2$", r"$\sigma_w^2$"],
    "True": [sigma_v2_true, sigma_w2_true],
    "EM estimate": [sigma_v2_hat, sigma_w2_hat],
})
df

# Log-likelihood (should be non-decreasing up to numerical noise)
plt.figure()
plt.plot(ll, marker="o")
plt.title("EM log-likelihood over iterations")
plt.xlabel("iteration"); plt.ylabel("log-likelihood")
plt.tight_layout()
plt.show()

# Data and smoothed state with ±1σ band
t = np.arange(len(x))
m = out["m_smooth"]
std = np.sqrt(out["P_smooth"])

plt.figure()
plt.plot(x, label="observations $x_t$")
plt.plot(y, label="true state $y_t$")
plt.plot(m, label="smoothed mean $\hat{y}_{t|T}$")
plt.fill_between(t, m - 2*std, m + 2*std, alpha=0.3, facecolor="tab:green", 
                  edgecolor="none", zorder=0,label="smoother ±2σ")
plt.legend()
plt.title("Local level model simulation and estimation")
plt.xlabel("t"); plt.ylabel("value")
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression, Ridge, RidgeCV, Lasso, LassoCV
from sklearn.metrics import mean_squared_error as mse

FIGURES_DIR = '../markdown/figures'

# Synthetic dataset: 20 features, 10 informative
n_samples, n_features, n_informative = 1000, 20, 10
rng = np.random.RandomState(0)
X, y, coef = make_regression(
    n_samples, n_features, n_informative=n_informative,
    noise=10, random_state=rng, coef=True
)
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=rng)
ftr_labels = [str(i + 1) for i in range(n_features)]

# Baseline OLS
reg = LinearRegression().fit(X_train, y_train)
print(f"OLS  MSE: {mse(y_test, reg.predict(X_test)):.2f}")

# Ridge with cross-validation
ridge_cv = RidgeCV(cv=5).fit(X_train, y_train)
print(f"Ridge MSE: {mse(y_test, ridge_cv.predict(X_test)):.2f}")

# Lasso with cross-validation
lasso_cv = LassoCV(cv=5, max_iter=10000).fit(X_train, y_train)
print(f"Lasso MSE: {mse(y_test, lasso_cv.predict(X_test)):.2f}")

# Ridge coefficient paths
n_alphas = 200
alphas_ridge = np.logspace(-1, 5, n_alphas)
ridge_coefs = []
for a in alphas_ridge:
    r = Ridge(alpha=a).fit(X_train, y_train)
    ridge_coefs.append(r.coef_)

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(alphas_ridge, ridge_coefs)
ax.set_xscale('log')
ax.set_xlabel('Regularization strength $\\lambda$')
ax.set_ylabel('Coefficient value')
ax.set_title('Ridge — coefficients as a function of regularization')
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/BLR_ridge_path.png', dpi=150)
plt.close(fig)
print('Saved BLR_ridge_path.png')

# Lasso coefficient paths
n_alphas = 200
alphas_lasso = np.logspace(-2, 4, n_alphas)
lasso_coefs = []
for a in alphas_lasso:
    l = Lasso(alpha=a, max_iter=10000).fit(X_train, y_train)
    lasso_coefs.append(l.coef_)

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(alphas_lasso, lasso_coefs)
ax.set_xscale('log')
ax.set_xlabel('Regularization strength $\\lambda$')
ax.set_ylabel('Coefficient value')
ax.set_title('Lasso — coefficients as a function of regularization')
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/BLR_lasso_path.png', dpi=150)
plt.close(fig)
print('Saved BLR_lasso_path.png')

from scipy import stats

class BayesLinRegPrior:
    """Prior parameters for Bayesian linear regression with NIG conjugate prior."""
    def __init__(self, n_features, lam, alpha, beta, weights):
        self.n_features = n_features
        self.lam = lam        # prior covariance (V0)
        self.alpha = alpha    # IG shape
        self.beta = beta      # IG scale
        self.weights = weights  # prior mean (b0)


class BayesLinReg:
    """Bayesian linear regression with Normal-Inverse-Gamma conjugate prior."""
    def __init__(self, prior):
        self.n_features = prior.n_features
        self.lam = prior.lam.copy()
        self.alpha = prior.alpha
        self.beta = float(prior.beta)
        self.weights = prior.weights.copy()

    def learn(self, x, y):
        x = np.atleast_2d(x)
        y = np.atleast_1d(y)
        lam0_inv = np.linalg.inv(self.lam)
        lamN_inv = lam0_inv + x.T @ x
        lamN = np.linalg.inv(lamN_inv)
        wN = lamN @ (lam0_inv @ self.weights + x.T @ y)
        alphaN = self.alpha + len(y) / 2
        betaN = float(self.beta + 0.5 * (
            self.weights.T @ lam0_inv @ self.weights
            + y @ y - wN.T @ lamN_inv @ wN
        ))
        self.lam = lamN
        self.weights = wN
        self.beta = betaN
        self.alpha = alphaN
        return self

    def sample_weights(self):
        sigma2 = stats.invgamma.rvs(self.alpha, scale=self.beta)
        return stats.multivariate_normal.rvs(mean=self.weights, cov=sigma2 * self.lam)

    @property
    def sigma(self):
        return float(np.sqrt(stats.invgamma.mean(self.alpha, scale=self.beta)))

# Single-feature dataset for posterior visualisation
bias = 100
rng2 = np.random.RandomState(10)
X1, y1, coef1 = make_regression(
    1000, 1, n_informative=1, bias=bias, noise=20, random_state=rng2, coef=True
)

fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(18, 18))
n_plt = 0
for n_obs in range(100, 1000, 100):
    ix, iy = divmod(n_plt, 3)
    # Fresh model for each N (avoids double-counting data)
    prior = BayesLinRegPrior(
        n_features=1,
        lam=0.1 * np.identity(1),
        alpha=1, beta=1,
        weights=np.zeros(1)
    )
    blr = BayesLinReg(prior)
    blr.learn(X1[:n_obs], y1[:n_obs])

    ax = axes[ix, iy]
    ax.scatter(X1[:n_obs, 0], y1[:n_obs], c='k', s=5, alpha=0.3, zorder=10, label='Data')
    x_plot = np.linspace(float(X1[:n_obs].min()), float(X1[:n_obs].max()), 200)
    for _ in range(10):
        b = blr.sample_weights()
        ax.plot(x_plot, bias + float(b) * x_plot, alpha=0.6, linewidth=1)
    ax.set_title(f'N = {n_obs}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    n_plt += 1

fig.suptitle('Posterior samples of regression lines — Bayesian Linear Regression', fontsize=14)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/BLR_posterior_samples.png', dpi=150)
plt.close(fig)
print('Saved BLR_posterior_samples.png')

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import networkx as nx
import yfinance as yf
from sklearn.feature_selection import f_regression

FIGURES_DIR = '../markdown/figures'

# Download 1 year of daily prices for 10 IBEX 35 stocks
ibex35_tkrs = ["MEL.MC","ENG.MC","ELE.MC","BBVA.MC","SAN.MC",
               "ACX.MC","AMS.MC","SAB.MC","MTS.MC","IAG.MC"]
raw = yf.download(tickers=" ".join(ibex35_tkrs), period="1y",
                  interval="1d", auto_adjust=True, progress=False)

# Handle multi-level columns from multi-ticker download
if isinstance(raw.columns, pd.MultiIndex):
    data = raw["Close"]
else:
    data = raw
data = data.dropna(axis=1, how="all").pct_change().dropna()

# Build undirected graph from pairwise F-test on correlations (1% level)
G = nx.Graph()
tickers = list(data.columns)
for tkr in tickers:
    F, pv = f_regression(data, data[tkr])
    for i, tkr2 in enumerate(tickers):
        if tkr2 != tkr and pv[i] < 0.01:
            G.add_edge(tkr, tkr2)

pos = nx.spring_layout(G, scale=2, seed=42)
fig, ax = plt.subplots(figsize=(8, 8))
nx.draw_networkx(G, pos=pos, ax=ax,
                 font_size=9, node_size=2200,
                 node_color="white", edgecolors="black",
                 linewidths=2, width=2)
ax.axis("off")
ax.set_title("IBEX 35 — undirected correlation graph (1-year daily returns)")
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/PGM_ibex35_correlation.png', dpi=150)
plt.close(fig)
print("Saved PGM_ibex35_correlation.png")

banks        = ["SAN.MC","BBVA.MC","SAB.MC"]
construction = ["MTS.MC","ACX.MC"]
travel       = ["IAG.MC","AMS.MC","MEL.MC"]
energy       = ["ENG.MC","ELE.MC"]

G2 = nx.DiGraph()
for tkr in tickers:
    G2.add_edge("SP Macro", tkr)
    if tkr in banks:        G2.add_edge("Banks", tkr)
    if tkr in construction: G2.add_edge("Construction", tkr)
    if tkr in travel:       G2.add_edge("Travel", tkr)
    if tkr in energy:       G2.add_edge("Energy", tkr)

pos2 = nx.spring_layout(G2, scale=5, seed=7)
fig, ax = plt.subplots(figsize=(10, 10))
nx.draw_networkx(G2, pos=pos2, ax=ax,
                 node_size=3000, node_color="white",
                 edgecolors="black", font_size=8,
                 verticalalignment="center_baseline")
ax.axis("off")
ax.set_title("IBEX 35 — directed causal graph (macro + sector factors)")
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/PGM_ibex35_causal.png', dpi=150)
plt.close(fig)
print("Saved PGM_ibex35_causal.png")

# Build a simple 6-node graph to illustrate MRF clique factorisation
G3 = nx.Graph()
G3.add_edges_from(nx.path_graph(6).edges())
G3.add_edge(0, 2)
G3.add_edge(5, 3)

pos3 = nx.spring_layout(G3, seed=0)
fig, ax = plt.subplots(figsize=(6, 4))
nx.draw_networkx(G3, pos=pos3, ax=ax,
                 font_size=12, node_size=1800,
                 node_color="white", edgecolors="black",
                 linewidths=2, width=2)
ax.axis("off")
ax.set_title("Undirected graph — three maximum cliques")
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/PGM_mrf_example.png', dpi=150)
plt.close(fig)

cliques = list(nx.find_cliques(G3))
print("Maximum cliques:", cliques)
print("Saved PGM_mrf_example.png")

# Demand-driven pricing Bayesian Network — with wealth confounder
G4 = nx.DiGraph()
G4.add_edge("p", "d")
G4.add_edge("W", "d")
G4.add_edge("W", "p")   # wealth also influences observed pricing policy

pos4 = {"W": (0, 1), "p": (-1, 0), "d": (1, 0)}
fig, ax = plt.subplots(figsize=(5, 4))
nx.draw_networkx(G4, pos=pos4, ax=ax,
                 node_size=2500, node_color="white",
                 edgecolors="black", font_size=12,
                 verticalalignment="center_baseline")
ax.axis("off")
ax.set_title("Bayesian Network: demand-driven pricing")
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/PGM_bn_pricing.png', dpi=150)
plt.close(fig)
print("Saved PGM_bn_pricing.png")

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Left: chain — conditioning on middle node blocks path
G_chain = nx.path_graph(4, create_using=nx.DiGraph)
colors_chain = ["white", "white", "lightgrey", "white"]
pos_chain = {0: (0,0), 1: (1,0), 2: (2,0), 3: (3,0)}
nx.draw_networkx(G_chain, pos=pos_chain, ax=axes[0],
                 node_color=colors_chain, edgecolors="black",
                 node_size=1800, font_size=14, linewidths=2, width=2)
axes[0].axis("off")
axes[0].set_title("Chain: node 2 conditioned (grey)\n"
                  "→ nodes 0 and 3 become d-separated")

# Right: confounder — conditioning on common cause blocks indirect path
G_conf = nx.DiGraph()
G_conf.add_edge("temp", "ice-cream")
G_conf.add_edge("temp", "crime")
pos_conf = {"temp": (1, 1), "ice-cream": (0, 0), "crime": (2, 0)}
colors_conf = ["lightgrey", "white", "white"]
nx.draw_networkx(G_conf, pos=pos_conf, ax=axes[1],
                 node_color=colors_conf, edgecolors="black",
                 node_size=3000, font_size=10, linewidths=2, width=2)
axes[1].axis("off")
axes[1].set_title("Confounder: conditioning on \"temp\" (grey)\n"
                  "→ ice-cream and crime become d-separated")

fig.suptitle("d-separation examples", fontsize=13)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/PGM_dsep.png', dpi=150)
plt.close(fig)
print("Saved PGM_dsep.png")

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Before intervention: P → L
G_before = nx.DiGraph()
G_before.add_edge("Pressure (P)", "Barometer (L)")
pos_b = {"Pressure (P)": (0, 0), "Barometer (L)": (1, 0)}
nx.draw_networkx(G_before, pos=pos_b, ax=axes[0],
                 node_color="white", edgecolors="black",
                 node_size=4000, font_size=9, linewidths=2, width=2)
axes[0].axis("off")
axes[0].set_title("Original graph:\nP → L")

# After do(L): incoming edge to L removed
G_after = nx.DiGraph()
G_after.add_node("Pressure (P)")
G_after.add_node("Barometer (L)")
nx.draw_networkx(G_after, pos=pos_b, ax=axes[1],
                 node_color="white", edgecolors="black",
                 node_size=4000, font_size=9, linewidths=2, width=2)
axes[1].axis("off")
axes[1].set_title("After do(L): edge P → L removed\n"
                  "→  P(P | do(L)) = P(P)")

fig.suptitle("Do-operator: intervening on the barometer level", fontsize=12)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/PGM_causal_barometer.png', dpi=150)
plt.close(fig)
print("Saved PGM_causal_barometer.png")

from scipy.stats import bernoulli, beta as beta_dist, norm
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

FIGURES_DIR = '../markdown/figures'

np.random.seed(42)
prior_alpha, prior_beta = 2, 2
p_true = 0.3
r = bernoulli.rvs(p_true, size=100, random_state=42)

def posterior_laplace(x, alpha_p, beta_p, r):
    N, n_H = len(r), sum(r)
    if N <= 1:
        return np.zeros_like(x)
    p_map = (n_H + alpha_p - 1) / (N + alpha_p + beta_p - 2)
    q_map = 1 - p_map
    std = np.sqrt(p_map * q_map / (N + alpha_p + beta_p - 2))
    return norm.pdf(x, loc=p_map, scale=std)

x_grid = np.linspace(0, 1, 1000)
fig, axes = plt.subplots(2, 2, figsize=(11, 9))
for idx, n in enumerate([0, 33, 66, 100]):
    ax = axes[idx // 2, idx % 2]
    r_trunc = r[:n]
    post_a = prior_alpha + sum(r_trunc)
    post_b = prior_beta + len(r_trunc) - sum(r_trunc)
    ax.plot(x_grid, beta_dist.pdf(x_grid, post_a, post_b),
            label="Exact (Beta)", lw=2)
    ax.plot(x_grid, posterior_laplace(x_grid, prior_alpha, prior_beta, r_trunc),
            label="Laplace (Gaussian)", lw=2, linestyle="--")
    ax.axvline(p_true, color="grey", linestyle=":", label=f"True p={p_true}")
    ax.set_title(f"N = {n}")
    ax.set_xlabel("p")
    ax.legend(fontsize=8)

fig.suptitle("Coin toss posterior: exact Beta vs Laplace approximation\n"
             f"(prior α={prior_alpha}, β={prior_beta})", fontsize=12)
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/BML_laplace_approximation.png', dpi=150)
plt.close(fig)
print("Saved BML_laplace_approximation.png")

from scipy.special import beta as beta_fn

np.random.seed(0)
prior_alpha, prior_beta = 2, 2
p_single = 0.3
r_h1 = bernoulli.rvs(p_single, size=500, random_state=0)

def bayes_factor(alpha_p, beta_p, r):
    N, n_H = len(r), sum(r)
    r_odd, r_even = r[::2], r[1::2]
    N_odd, N_even = len(r_odd), len(r_even)
    num = beta_fn(alpha_p + N, beta_p + N - n_H) * beta_fn(alpha_p, beta_p)
    den = (beta_fn(alpha_p + N_odd,  beta_p + N_odd  - sum(r_odd))
         * beta_fn(alpha_p + N_even, beta_p + N_even - sum(r_even)))
    return num / den

def bic_ll(alpha_p, beta_p, r):
    N, n_H = len(r), sum(r)
    if N < 2: return 0.0
    p_map = (n_H + alpha_p - 1) / (N + alpha_p + beta_p - 2)
    q_map = 1 - p_map
    return (n_H * np.log(max(p_map, 1e-15))
            + (N - n_H) * np.log(max(q_map, 1e-15)))

def bayes_factor_bic(alpha_p, beta_p, r):
    N = len(r)
    r_odd, r_even = r[::2], r[1::2]
    bic1 = bic_ll(alpha_p, beta_p, r) - 0.5 * np.log(max(N, 1))
    bic2 = (bic_ll(alpha_p, beta_p, r_odd) + bic_ll(alpha_p, beta_p, r_even)
            - np.log(max(N, 1)))
    return np.exp(bic1 - bic2)

ns = range(3, len(r_h1))
bfs_exact = [bayes_factor(prior_alpha, prior_beta, r_h1[:i]) for i in ns]
bfs_bic   = [bayes_factor_bic(prior_alpha, prior_beta, r_h1[:i]) for i in ns]

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(ns, bfs_exact, label="Exact")
ax.plot(ns, bfs_bic,   label="BIC approx", alpha=0.7)
ax.axhline(1, color="grey", linestyle=":")
ax.axhline(10, color="green", linestyle="--", alpha=0.5, label="Strong evidence (BF=10)")
ax.set_xlabel("N (observations)")
ax.set_ylabel("Bayes Factor $BF_{1,2}$")
ax.set_title("Bayes Factor evolution — data from H1 (single coin, p=0.3)")
ax.legend()
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/BML_bayes_factor_h1.png', dpi=150)
plt.close(fig)
print("Saved BML_bayes_factor_h1.png")

np.random.seed(1)
p1, p2 = 0.3, 0.7
r1 = bernoulli.rvs(p1, size=250, random_state=1)
r2 = bernoulli.rvs(p2, size=250, random_state=2)
r_h2 = np.array([x for pair in zip(r1, r2) for x in pair])

ns2 = range(3, len(r_h2))
bfs2_exact = [bayes_factor(prior_alpha, prior_beta, r_h2[:i]) for i in ns2]
bfs2_bic   = [bayes_factor_bic(prior_alpha, prior_beta, r_h2[:i]) for i in ns2]

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(ns2, bfs2_exact, label="Exact")
ax.plot(ns2, bfs2_bic,   label="BIC approx", alpha=0.7)
ax.axhline(1, color="grey", linestyle=":")
ax.axhline(0.1, color="red", linestyle="--", alpha=0.5, label="Strong evidence for H2 (BF=0.1)")
ax.set_xlabel("N (observations)")
ax.set_ylabel("Bayes Factor $BF_{1,2}$")
ax.set_title("Bayes Factor evolution — data from H2 (two alternating coins, p₁=0.3, p₂=0.7)")
ax.legend()
fig.tight_layout()
fig.savefig(f'{FIGURES_DIR}/BML_bayes_factor_h2.png', dpi=150)
plt.close(fig)
print("Saved BML_bayes_factor_h2.png")