{ "cells": [ { "cell_type": "markdown", "id": "be27992f", "metadata": {}, "source": [ "# Optimal Transport in linear ICA applied to Econometrics\n", "\n", "#### In this application we test our methodology to find ICs in a simulated price discovery case between distinct markets" ] }, { "cell_type": "code", "execution_count": 1, "id": "ce700a45", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import torch\n", "import scipy.stats\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import linear_sum_assignment\n", "\n", "from wasserstein_ica import WassersteinICA" ] }, { "cell_type": "code", "execution_count": 2, "id": "efd6fc0a", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "\n", "# Define a consistent Thesis Theme\n", "def set_thesis_theme():\n", " # Academic, colorblind-friendly palette\n", " # Blue, Orange, Green, Red, Purple, Brown\n", " thesis_colors = ['#0173B2', '#DE8F05', '#029E73', '#D55E00', '#CC78BC', '#CA9161']\n", " \n", " mpl.rcParams.update({\n", " # Figure and Layout\n", " 'figure.figsize': (8, 5),\n", " 'figure.dpi': 300, # High resolution for print\n", " 'axes.prop_cycle': mpl.cycler(color=thesis_colors),\n", " \n", " # Grid lines (light and unobtrusive)\n", " 'axes.grid': True,\n", " 'grid.alpha': 0.3,\n", " 'grid.linestyle': '--',\n", " 'axes.axisbelow': True, # Grid goes behind data\n", " \n", " # Spines (remove top and right borders for a cleaner look)\n", " 'axes.spines.top': False,\n", " 'axes.spines.right': False,\n", " \n", " # Fonts and Text\n", " 'font.size': 11,\n", " 'axes.titlesize': 13,\n", " 'axes.labelsize': 12,\n", " 'xtick.labelsize': 10,\n", " 'ytick.labelsize': 10,\n", " \n", " # Legends\n", " 'legend.frameon': False, # No box around the legend\n", " 'legend.fontsize': 10,\n", " \n", " # Lines\n", " 'lines.linewidth': 2.0\n", " })\n", "\n", "# Run this before plotting\n", "set_thesis_theme()" ] }, { "cell_type": "code", "execution_count": 3, "id": "03a48eff", "metadata": {}, "outputs": [], "source": [ "def apply_economic_constraints(W_est, W_white):\n", " \"\"\"\n", " Post-process the statistical ICA output to enforce economic meaning.\n", " \n", " Addresses the critique that standard ICA results are permutation/sign indefinite.\n", " We apply:\n", " 1. Dominant Diagonal: The shock impacting Market i the most is labeled Shock i.\n", " 2. Positive Spillover: Structural shocks are normalized to have positive effects.\n", " \n", " Args:\n", " W_est (Torch Tensor): Unmixing matrix from ICA (Whitened space).\n", " W_white (Torch Tensor): Whitening matrix used.\n", " \n", " Returns:\n", " B_identified (Numpy Array): The economically identified Mixing Matrix.\n", " \"\"\"\n", " # 1. Recover the Raw Mixing Matrix (B_raw)\n", " # The total estimated unmixing is W_tot = W_est @ W_white\n", " # The mixing matrix B is the inverse: X = B * S\n", " W_total = torch.matmul(W_est, W_white)\n", " B_raw = torch.linalg.inv(W_total).cpu().numpy()\n", " \n", " # 2. Permutation Constraint (Dominant Diagonal)\n", " # We use the Hungarian Algorithm (linear_sum_assignment) to match \n", " # Shocks (columns) to Markets (rows) such that the diagonal magnitude is maximized.\n", " # We pass negative absolute values because the algorithm minimizes cost.\n", " row_ind, col_ind = linear_sum_assignment(-np.abs(B_raw))\n", " \n", " # Reorder the columns (shocks) of B\n", " B_permuted = B_raw[:, col_ind]\n", " \n", " # 3. Sign Constraint (Positive Impact)\n", " # Ensure the diagonal elements are positive (B_ii > 0).\n", " # If B_ii is negative, we assume we found the \"negative\" shock and flip the column.\n", " diagonals = np.diag(B_permuted)\n", " sign_corrections = np.sign(diagonals)\n", " \n", " # Broadcast multiplication to flip columns where sign is negative\n", " B_identified = B_permuted * sign_corrections[np.newaxis, :]\n", " \n", " return B_identified" ] }, { "cell_type": "code", "execution_count": 4, "id": "c8662488", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Running Experiment on cpu...\n", "Target True IS: [0.12 0.24 0.64]\n", "Constraint Mode: Dominant Diagonal & Positive Signs\n" ] } ], "source": [ "# ==============================================================================\n", "# EXPERIMENT SETUP (Replicating Zema & Cordoni DGP)\n", "# ==============================================================================\n", "\n", "# Simulation Settings\n", "N_SIMULATIONS = 500 # Number of Monte Carlo runs\n", "T_SAMPLES = 5000 # Observations per run (Daily returns)\n", "N_MARKETS = 3\n", "DEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n", "\n", "# Define the \"True\" Structure\n", "# We target Information Shares: [0.12, 0.24, 0.64]\n", "# Assuming a standard VECM where the orthogonal complement to beta is [1, 1, 1],\n", "# The IS is proportional to the squared column sums of B.\n", "target_is = np.array([0.12, 0.24, 0.64])\n", "psi = np.array([1.0, 1.0, 1.0]) # Long-run impact vector\n", "\n", "# Construct a True Mixing Matrix B that yields these shares.\n", "# For simplicity, we use a diagonal B scaled by sqrt(target).\n", "# (Note: Even with a diagonal B_true, the ICA has to work because we whiten the data,\n", "# effectively destroying the structure, and ICA must rotate it back).\n", "B_true = np.diag(np.sqrt(target_is)) \n", "\n", "print(f\"Running Experiment on {DEVICE}...\")\n", "print(f\"Target True IS: {target_is}\")\n", "print(f\"Constraint Mode: Dominant Diagonal & Positive Signs\")\n", "\n", "results_is = []" ] }, { "cell_type": "code", "execution_count": 5, "id": "6873a940", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Simulation 10/500 completed.\n", "Simulation 20/500 completed.\n", "Simulation 30/500 completed.\n", "Simulation 40/500 completed.\n", "Simulation 50/500 completed.\n", "Simulation 60/500 completed.\n", "Simulation 70/500 completed.\n", "Simulation 80/500 completed.\n", "Simulation 90/500 completed.\n", "Simulation 100/500 completed.\n", "Simulation 110/500 completed.\n", "Simulation 120/500 completed.\n", "Simulation 130/500 completed.\n", "Simulation 140/500 completed.\n", "Simulation 150/500 completed.\n", "Simulation 160/500 completed.\n", "Simulation 170/500 completed.\n", "Simulation 180/500 completed.\n", "Simulation 190/500 completed.\n", "Simulation 200/500 completed.\n", "Simulation 210/500 completed.\n", "Simulation 220/500 completed.\n", "Simulation 230/500 completed.\n", "Simulation 240/500 completed.\n", "Simulation 250/500 completed.\n", "Simulation 260/500 completed.\n", "Simulation 270/500 completed.\n", "Simulation 280/500 completed.\n", "Simulation 290/500 completed.\n", "Simulation 300/500 completed.\n", "Simulation 310/500 completed.\n", "Simulation 320/500 completed.\n", "Simulation 330/500 completed.\n", "Simulation 340/500 completed.\n", "Simulation 350/500 completed.\n", "Simulation 360/500 completed.\n", "Simulation 370/500 completed.\n", "Simulation 380/500 completed.\n", "Simulation 390/500 completed.\n", "Simulation 400/500 completed.\n", "Simulation 410/500 completed.\n", "Simulation 420/500 completed.\n", "Simulation 430/500 completed.\n", "Simulation 440/500 completed.\n", "Simulation 450/500 completed.\n", "Simulation 460/500 completed.\n", "Simulation 470/500 completed.\n", "Simulation 480/500 completed.\n", "Simulation 490/500 completed.\n", "Simulation 500/500 completed.\n" ] } ], "source": [ "for sim in range(N_SIMULATIONS):\n", " \n", " # --- A. Data Generation (Student-t Innovations) ---\n", " # Independent shocks with degrees of freedom 5, 6, 7 (Heavy tails)\n", " s1 = np.random.standard_t(df=5, size=T_SAMPLES)\n", " s2 = np.random.standard_t(df=6, size=T_SAMPLES)\n", " s3 = np.random.standard_t(df=7, size=T_SAMPLES)\n", " \n", " # Stack and Normalize to unit variance (Standard ICA assumption)\n", " S = np.vstack([s1, s2, s3])\n", " S = S / S.std(axis=1, keepdims=True)\n", " \n", " # Mix the signals\n", " X_np = B_true @ S\n", " X_torch = torch.tensor(X_np, dtype=torch.float32, device=DEVICE)\n", " \n", " # --- B. Wasserstein ICA ---\n", " model = WassersteinICA(X_torch)\n", " model.whiten()\n", " \n", " dim = X_torch.shape[0]\n", " num_restarts = min(dim * 4, 150)\n", " \n", " # \n", " # OPTIMIZATION PHASE 1: Deflationary SGD\n", " # Sequentially extract components to form a robust initialization matrix\n", " W_init_list = []\n", " prev_components = torch.empty((0, dim), device=DEVICE)\n", " \n", " for k in range(dim):\n", " w_k, _ = model.optimize_wasserstein2(\n", " prev_components=prev_components,\n", " continuous=True,\n", " max_iter=200,\n", " lr=0.1,\n", " n_restarts=num_restarts\n", " )\n", " W_init_list.append(w_k)\n", " prev_components = torch.cat([prev_components, w_k.unsqueeze(0)], dim=0)\n", " \n", " W_init = torch.stack(W_init_list)\n", " \n", " # OPTIMIZATION PHASE 2: Symmetric Stiefel Optimization\n", " # Refine the initial matrix simultaneously on the Stiefel manifold\n", " W_opt = model.optimize_symmetric(\n", " n_components=dim,\n", " optimizer='stiefel',\n", " init_w=W_init,\n", " use_sinkhorn=False,\n", " max_iter=300,\n", " lr=1.0\n", " )\n", " \n", " # --- C. Economic Identification ---\n", " B_final = apply_economic_constraints(W_opt, model.W_white)\n", " \n", " # --- D. Calculate Information Share ---\n", " # Calculate the long-run impact of each shock: Psi * B_column\n", " # Since Psi=[1,1,1], this is just the column sum.\n", " impacts = psi @ B_final\n", " \n", " # Variance contribution = impact^2\n", " variances = impacts ** 2\n", " \n", " # Information Share = Variance / Total Variance\n", " estimated_is = variances / np.sum(variances)\n", " results_is.append(estimated_is)\n", " \n", " if (sim + 1) % 10 == 0:\n", " print(f\"Simulation {sim+1}/{N_SIMULATIONS} completed.\")" ] }, { "cell_type": "code", "execution_count": 6, "id": "fc361326", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "============================================================\n", "Market / Source | True IS | Est. IS (Mean) | Std Dev \n", "------------------------------------------------------------\n", "Market 1 | 0.1200 | 0.1199 | 0.0169 \n", "Market 2 | 0.2400 | 0.2380 | 0.0241 \n", "Market 3 | 0.6400 | 0.6421 | 0.0268 \n", "============================================================\n" ] } ], "source": [ "#==============================================================================\n", "# RESULTS TABLE\n", "# ==============================================================================\n", "results_is = np.array(results_is)\n", "mean_is = np.mean(results_is, axis=0)\n", "std_is = np.std(results_is, axis=0)\n", "\n", "print(\"\\n\" + \"=\"*60)\n", "print(f\"{'Market / Source':<15} | {'True IS':<10} | {'Est. IS (Mean)':<15} | {'Std Dev':<10}\")\n", "print(\"-\" * 60)\n", "\n", "for i in range(N_MARKETS):\n", " print(f\"Market {i+1:<14} | {target_is[i]:<10.4f} | {mean_is[i]:<15.4f} | {std_is[i]:<10.4f}\")\n", "\n", "print(\"=\"*60)" ] }, { "cell_type": "code", "execution_count": 10, "id": "1759bbd6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Successfully saved huge figure to shock_recovery_scatter.pdf\n" ] } ], "source": [ "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "# Colorblind-friendly palette EXCLUDING Blue and Orange (reserved for FastICA/OT-ICA)\n", "# Market 1: Teal/Green, Market 2: Purple/Pink, Market 3: Brown\n", "market_colors = ['#029E73', '#CC78BC', '#CA9161']\n", "market_labels = ['Market 1 (Target 0.12)', 'Market 2 (Target 0.24)', 'Market 3 (Target 0.64)']\n", "targets = [0.12, 0.24, 0.64]\n", "\n", "# ---------------------------------------------------------\n", "# Plot 1: Information Share Distribution (Boxplot)\n", "# ---------------------------------------------------------\n", "plt.figure(figsize=(14, 10))\n", "\n", "# Use seaborn for a cleaner boxplot that natively accepts our color palette\n", "sns.boxplot(data=results_is, palette=market_colors, width=0.5, \n", " boxprops=dict(alpha=0.6), medianprops=dict(color=\"black\", linewidth=2.5))\n", "\n", "plt.title(f\"Information Share Distribution (N={N_SIMULATIONS})\\nMetric: Wasserstein-ICA\", \n", " fontsize=32, pad=20, fontweight='bold')\n", "plt.ylabel(\"Information Share\", fontsize=30, fontweight='bold')\n", "plt.xticks([0, 1, 2], ['Market 1', 'Market 2', 'Market 3'], fontsize=28)\n", "plt.yticks(fontsize=28)\n", "\n", "# Draw True IS lines with matching colors\n", "for i in range(3):\n", " plt.axhline(targets[i], color=market_colors[i], linestyle='--', \n", " linewidth=3, alpha=0.9, label=f'True {market_labels[i]}')\n", "\n", "plt.legend(fontsize=28, loc='center left', frameon=False)\n", "plt.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", "\n", "file_name_1 = \"shock_recovery_scatter.pdf\"\n", "plt.savefig(file_name_1, format='pdf', bbox_inches='tight')\n", "print(f\"Successfully saved huge figure to {file_name_1}\")\n", "plt.close()" ] }, { "cell_type": "code", "execution_count": 8, "id": "c2a9aae7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Successfully saved huge figure to econ_applicatn_is_distribution.pdf\n" ] } ], "source": [ "# ---------------------------------------------------------\n", "# Plot 2: Distribution of Estimated Information Shares (KDE)\n", "# ---------------------------------------------------------\n", "plt.figure(figsize=(14, 10))\n", "\n", "for i in range(3):\n", " # Extract the column for Market i across all simulations\n", " data = results_is[:, i]\n", " \n", " # Plot KDE with matching colors\n", " sns.kdeplot(data, color=market_colors[i], fill=True, alpha=0.4, linewidth=3.5, label=market_labels[i])\n", " # Add vertical line for the true target\n", " plt.axvline(targets[i], color=market_colors[i], linestyle='--', linewidth=3, alpha=0.9)\n", "\n", "plt.title(f'Empirical Distribution of Estimated Information Shares\\n(Wasserstein-ICA, N={N_SIMULATIONS}, T={T_SAMPLES})', \n", " fontsize=30, pad=20, fontweight='bold')\n", "plt.xlabel('Information Share', fontsize=28, fontweight='bold')\n", "plt.ylabel('Probability Density', fontsize=28, fontweight='bold')\n", "plt.xticks(fontsize=24)\n", "plt.yticks(fontsize=24)\n", "\n", "plt.legend(fontsize=28, loc='upper center', frameon=False)\n", "plt.grid(True, alpha=0.3)\n", "plt.tight_layout()\n", "\n", "file_name_2 = \"econ_applicatn_is_distribution.pdf\"\n", "plt.savefig(file_name_2, format='pdf', bbox_inches='tight')\n", "print(f\"Successfully saved huge figure to {file_name_2}\")\n", "plt.close()" ] }, { "cell_type": "markdown", "id": "d51b5ef0", "metadata": {}, "source": [ "### Economic Identification Strategy\n", "Unlike the statistical output of raw ICA, which is permutation and sign indeterminate, we impose two strict economic constraints to recover the structural parameters:\n", "\n", "1. **Dominant Diagonal (Permutation):** We strictly enforce that the structural shock assigned to Market $i$ must be the one that contributes the maximum variance to that specific market. This is implemented via the Hungarian Algorithm to solve the assignment problem.\n", "2. **Positive Impact (Sign):** We assume positive spillovers, meaning a structural innovation (good news) must have a positive initial impact on its own price. Negative coefficients are sign-flipped." ] }, { "cell_type": "code", "execution_count": null, "id": "1115e192", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "ot_ica", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.4" } }, "nbformat": 4, "nbformat_minor": 5 }