{ "cells": [ { "cell_type": "code", "execution_count": 23, "id": "338730e2-a6de-4d4f-b438-efe3feb139ab", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import plotly.graph_objects as go\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "cfd11919-0941-400e-a516-72871881f733", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_1311/1940519970.py:1: DtypeWarning: Columns (1,2,3,4) have mixed types. Specify dtype option on import or set low_memory=False.\n", " stocks=pd.read_csv('stocks.csv')\n", "/tmp/ipykernel_1311/1940519970.py:2: DtypeWarning: Columns (1,2,3,4) have mixed types. Specify dtype option on import or set low_memory=False.\n", " flows = pd.read_csv('flows.csv')\n" ] } ], "source": [ "stocks=pd.read_csv('stocks.csv')\n", "flows = pd.read_csv('flows.csv')" ] }, { "cell_type": "code", "execution_count": 3, "id": "b99e3402-fe26-4f4e-8c1c-5f07847bce94", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_1311/3613746644.py:1: DtypeWarning: Columns (1) have mixed types. Specify dtype option on import or set low_memory=False.\n", " merged = pd.read_csv('merged.csv')\n" ] } ], "source": [ "merged = pd.read_csv('merged.csv')" ] }, { "cell_type": "code", "execution_count": 4, "id": "34e5a815-7269-4312-bfe6-e2cd12595e57", "metadata": {}, "outputs": [], "source": [ "# 1. Prepare stock dataset ISIN-by-ISIN\n", "stocks_isin = stocks[[\n", " \"Registrar Account - ID\",\n", " \"Product - Isin\",\n", " \"Centralisation Date\",\n", " \"Quantity - AUM\"\n", "]].copy()\n", "\n", "stocks_isin[\"Centralisation Date\"] = pd.to_datetime(stocks_isin[\"Centralisation Date\"])\n", "\n", "stocks_isin = stocks_isin.sort_values(\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"]\n", ")\n", "\n", "# 2. Prepare flows dataset ISIN-by-ISIN\n", "flows_isin = flows[[\n", " \"Registrar Account - ID\",\n", " \"Product - Isin\",\n", " \"Centralisation Date\",\n", " \"Quantity - NetFlows\"\n", "]].copy()\n", "\n", "flows_isin[\"Centralisation Date\"] = pd.to_datetime(flows_isin[\"Centralisation Date\"])\n", "\n", "flows_isin = (\n", " flows_isin\n", " .groupby(\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"]\n", " )[\"Quantity - NetFlows\"]\n", " .sum()\n", " .reset_index()\n", ")\n", "\n", "# 3. Merge stocks & flows ISIN-by-ISIN\n", "merged_isin = stocks_isin.merge(\n", " flows_isin,\n", " on=[\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"],\n", " how=\"left\"\n", ")\n", "\n", "merged_isin[\"Quantity - NetFlows\"] = merged_isin[\"Quantity - NetFlows\"].fillna(0)\n", "\n", "# 4. Compute expected stock per ISIN for each account\n", "merged_isin[\"prev_stock\"] = (\n", " merged_isin\n", " .groupby([\"Registrar Account - ID\", \"Product - Isin\"])[\"Quantity - AUM\"]\n", " .shift(1)\n", ")\n", "\n", "merged_isin[\"prev_netflows\"] = (\n", " merged_isin\n", " .groupby([\"Registrar Account - ID\", \"Product - Isin\"])[\"Quantity - NetFlows\"]\n", " .shift(1)\n", " .fillna(0)\n", ")\n", "\n", "merged_isin[\"expected_stock\"] = (\n", " merged_isin[\"prev_stock\"] + merged_isin[\"prev_netflows\"]\n", ")\n", "\n", "# 5. Detect ruptures ISIN-by-ISIN (no aggregation)\n", "TOL = 1e-6\n", "\n", "merged_isin[\"gap\"] = (\n", " merged_isin[\"Quantity - AUM\"] - merged_isin[\"expected_stock\"]\n", ")\n", "\n", "merged_isin[\"rupture_flag\"] = (\n", " merged_isin[\"prev_stock\"].notna()\n", " & (merged_isin[\"gap\"].abs() > TOL)\n", ")\n", "\n", "# 6. Summarize ruptures per (Account, ISIN)\n", "rupture_isin_summary = (\n", " merged_isin\n", " .groupby([\"Registrar Account - ID\", \"Product - Isin\"])\n", " .agg(\n", " n_ruptures=(\"rupture_flag\", \"sum\"),\n", " obs=(\"rupture_flag\", \"count\"),\n", " rupture_ratio=(\"rupture_flag\", \"mean\"),\n", " max_gap=(\"gap\", lambda x: x.abs().max())\n", " )\n", " .reset_index()\n", ")\n", "\n", "# Sort by worst ISIN trajectories\n", "rupture_isin_summary = rupture_isin_summary.sort_values(\n", " \"rupture_ratio\",\n", " ascending=False\n", ")" ] }, { "cell_type": "markdown", "id": "16213cb2-07d8-4e82-b9bb-252554ec47b9", "metadata": {}, "source": [ "# Détection des ruptures" ] }, { "cell_type": "code", "execution_count": 5, "id": "78c3db70-e0b6-4de2-92ca-e29cf5bf6bd1", "metadata": {}, "outputs": [], "source": [ "# ============================================================\n", "# AUM–FLOW CONSISTENCY & RUPTURE DETECTION (FINAL VERSION)\n", "# ============================================================\n", "# ------------------------------------------------------------\n", "# 1. Keep relevant columns\n", "# ------------------------------------------------------------\n", "stocks_clean = stocks[[\n", " \"Registrar Account - ID\",\n", " \"Product - Isin\",\n", " \"Centralisation Date\",\n", " \"Quantity - AUM\"\n", "]].copy()\n", "\n", "flows_clean = flows[[\n", " \"Registrar Account - ID\",\n", " \"Product - Isin\",\n", " \"Centralisation Date\",\n", " \"Quantity - NetFlows\"\n", "]].copy()\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 2. Date formatting\n", "# ------------------------------------------------------------\n", "stocks_clean[\"Centralisation Date\"] = pd.to_datetime(stocks_clean[\"Centralisation Date\"])\n", "flows_clean[\"Centralisation Date\"] = pd.to_datetime(flows_clean[\"Centralisation Date\"])\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 3. Aggregate flows per day\n", "# ------------------------------------------------------------\n", "flows_clean = (\n", " flows_clean\n", " .groupby(\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"],\n", " as_index=False\n", " )[\"Quantity - NetFlows\"]\n", " .sum()\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# 4. Merge stocks and flows\n", "# ------------------------------------------------------------\n", "df = stocks_clean.merge(\n", " flows_clean,\n", " on=[\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"],\n", " how=\"left\"\n", ")\n", "\n", "df[\"Quantity - NetFlows\"] = df[\"Quantity - NetFlows\"].fillna(0)\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 5. Sort and compute expected stock\n", "# ------------------------------------------------------------\n", "df = df.sort_values(\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"]\n", ")\n", "\n", "df[\"prev_stock\"] = df.groupby(\n", " [\"Registrar Account - ID\", \"Product - Isin\"]\n", ")[\"Quantity - AUM\"].shift(1)\n", "\n", "df[\"prev_flows\"] = df.groupby(\n", " [\"Registrar Account - ID\", \"Product - Isin\"]\n", ")[\"Quantity - NetFlows\"].shift(1).fillna(0)\n", "\n", "df[\"expected_stock\"] = df[\"prev_stock\"] + df[\"prev_flows\"]\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 6. Compute gaps\n", "# ------------------------------------------------------------\n", "df[\"gap\"] = df[\"Quantity - AUM\"] - df[\"expected_stock\"]\n", "df[\"gap_abs\"] = df[\"gap\"].abs()\n", "df[\"gap_rel\"] = df[\"gap_abs\"] / df[\"expected_stock\"].abs().clip(lower=1)\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 7. Detect ruptures (economic rule)\n", "# ------------------------------------------------------------\n", "TAU_ABS = 10.0 # minimum absolute gap (shares)\n", "TAU_REL = 0.005 # minimum relative gap (0.5%)\n", "\n", "df[\"rupture_flag\"] = (\n", " df[\"prev_stock\"].notna()\n", " & (df[\"gap_abs\"] > TAU_ABS)\n", " & (df[\"gap_rel\"] > TAU_REL)\n", ")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 8. Remove end-of-sample false positives (edge effects)\n", "# ------------------------------------------------------------\n", "last_date = df[\"Centralisation Date\"].max()\n", "\n", "df[\"rupture_flag\"] = np.where(\n", " (df[\"rupture_flag\"]) & (df[\"Centralisation Date\"] == last_date),\n", " False,\n", " df[\"rupture_flag\"]\n", ")\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "a9783dc1-e225-4142-8b6f-6f9e620b4b3d", "metadata": {}, "outputs": [], "source": [ "# ------------------------------------------------------------\n", "# 9. ISIN-level summary (AFTER CLEANING)\n", "# ------------------------------------------------------------\n", "rupture_isin_summary = (\n", " df.groupby([\"Registrar Account - ID\", \"Product - Isin\"])\n", " .agg(\n", " n_ruptures=(\"rupture_flag\", \"sum\"),\n", " total_obs=(\"rupture_flag\", \"count\"),\n", " rupture_ratio=(\"rupture_flag\", \"mean\"),\n", " max_gap=(\"gap_abs\", \"max\")\n", " )\n", " .reset_index()\n", ")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 10. Account-level summary (AFTER CLEANING)\n", "# ------------------------------------------------------------\n", "rupture_summary = (\n", " df.groupby(\"Registrar Account - ID\")\n", " .agg(\n", " n_ruptures=(\"rupture_flag\", \"sum\"),\n", " total_obs=(\"rupture_flag\", \"count\"),\n", " rupture_ratio=(\"rupture_flag\", \"mean\"),\n", " max_gap=(\"gap_abs\", \"max\")\n", " )\n", " .reset_index()\n", ")\n", "\n", "\n", "# ------------------------------------------------------------\n", "# 11. Outputs\n", "# ------------------------------------------------------------\n", "df.to_csv(\"aum_flow_gaps.csv\", index=False)\n", "rupture_isin_summary.to_csv(\"rupture_isin_summary.csv\", index=False)\n", "rupture_summary.to_csv(\"rupture_summary.csv\", index=False)" ] }, { "cell_type": "code", "execution_count": 7, "id": "f5b62558-c27a-4428-a193-8b97e0ce6b6a", "metadata": {}, "outputs": [ { "data": { "application/vnd.plotly.v1+json": { "config": { "plotlyServerURL": "https://plot.ly" }, "data": [ { "hole": 0.45, "hoverinfo": "label+percent", "labels": [ "Clean / quasi-clean (≤1%)", "Moderate (1–10%)", "High (10–30%)", "Severe (>30%)" ], "textinfo": "percent", "type": "pie", "values": { "bdata": "AAAAAACASEAAAAAAAIBBQAAAAAAAAChAZmZmZmZmEEA=", "dtype": "f8" } } ], "layout": { "legend": { "orientation": "h", "title": { "text": "Rupture ratio" }, "x": 0.5, "xanchor": "center", "y": -0.15, "yanchor": "top" }, "template": { "data": { "bar": [ { "error_x": { "color": "#2a3f5f" }, "error_y": { "color": "#2a3f5f" }, "marker": { "line": { "color": "#E5ECF6", "width": 0.5 }, "pattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 } }, "type": "bar" } ], "barpolar": [ { "marker": { "line": { "color": "#E5ECF6", "width": 0.5 }, "pattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 } }, "type": "barpolar" } ], "carpet": [ { "aaxis": { "endlinecolor": "#2a3f5f", "gridcolor": "white", "linecolor": "white", "minorgridcolor": "white", "startlinecolor": "#2a3f5f" }, "baxis": { "endlinecolor": "#2a3f5f", "gridcolor": "white", "linecolor": "white", "minorgridcolor": "white", "startlinecolor": "#2a3f5f" }, "type": "carpet" } ], "choropleth": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "type": "choropleth" } ], "contour": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "contour" } ], "contourcarpet": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "type": "contourcarpet" } ], "heatmap": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "heatmap" } ], "histogram": [ { "marker": { "pattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 } }, "type": "histogram" } ], "histogram2d": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "histogram2d" } ], "histogram2dcontour": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "histogram2dcontour" } ], "mesh3d": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "type": "mesh3d" } ], "parcoords": [ { "line": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "parcoords" } ], "pie": [ { "automargin": true, "type": "pie" } ], "scatter": [ { "fillpattern": { "fillmode": "overlay", "size": 10, "solidity": 0.2 }, "type": "scatter" } ], "scatter3d": [ { "line": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatter3d" } ], "scattercarpet": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattercarpet" } ], "scattergeo": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattergeo" } ], "scattergl": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattergl" } ], "scattermap": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattermap" } ], "scattermapbox": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scattermapbox" } ], "scatterpolar": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatterpolar" } ], "scatterpolargl": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatterpolargl" } ], "scatterternary": [ { "marker": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "type": "scatterternary" } ], "surface": [ { "colorbar": { "outlinewidth": 0, "ticks": "" }, "colorscale": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "type": "surface" } ], "table": [ { "cells": { "fill": { "color": "#EBF0F8" }, "line": { "color": "white" } }, "header": { "fill": { "color": "#C8D4E3" }, "line": { "color": "white" } }, "type": "table" } ] }, "layout": { "annotationdefaults": { "arrowcolor": "#2a3f5f", "arrowhead": 0, "arrowwidth": 1 }, "autotypenumbers": "strict", "coloraxis": { "colorbar": { "outlinewidth": 0, "ticks": "" } }, "colorscale": { "diverging": [ [ 0, "#8e0152" ], [ 0.1, "#c51b7d" ], [ 0.2, "#de77ae" ], [ 0.3, "#f1b6da" ], [ 0.4, "#fde0ef" ], [ 0.5, "#f7f7f7" ], [ 0.6, "#e6f5d0" ], [ 0.7, "#b8e186" ], [ 0.8, "#7fbc41" ], [ 0.9, "#4d9221" ], [ 1, "#276419" ] ], "sequential": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ], "sequentialminus": [ [ 0, "#0d0887" ], [ 0.1111111111111111, "#46039f" ], [ 0.2222222222222222, "#7201a8" ], [ 0.3333333333333333, "#9c179e" ], [ 0.4444444444444444, "#bd3786" ], [ 0.5555555555555556, "#d8576b" ], [ 0.6666666666666666, "#ed7953" ], [ 0.7777777777777778, "#fb9f3a" ], [ 0.8888888888888888, "#fdca26" ], [ 1, "#f0f921" ] ] }, "colorway": [ "#636efa", "#EF553B", "#00cc96", "#ab63fa", "#FFA15A", "#19d3f3", "#FF6692", "#B6E880", "#FF97FF", "#FECB52" ], "font": { "color": "#2a3f5f" }, "geo": { "bgcolor": "white", "lakecolor": "white", "landcolor": "#E5ECF6", "showlakes": true, "showland": true, "subunitcolor": "white" }, "hoverlabel": { "align": "left" }, "hovermode": "closest", "mapbox": { "style": "light" }, "paper_bgcolor": "white", "plot_bgcolor": "#E5ECF6", "polar": { "angularaxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" }, "bgcolor": "#E5ECF6", "radialaxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" } }, "scene": { "xaxis": { "backgroundcolor": "#E5ECF6", "gridcolor": "white", "gridwidth": 2, "linecolor": "white", "showbackground": true, "ticks": "", "zerolinecolor": "white" }, "yaxis": { "backgroundcolor": "#E5ECF6", "gridcolor": "white", "gridwidth": 2, "linecolor": "white", "showbackground": true, "ticks": "", "zerolinecolor": "white" }, "zaxis": { "backgroundcolor": "#E5ECF6", "gridcolor": "white", "gridwidth": 2, "linecolor": "white", "showbackground": true, "ticks": "", "zerolinecolor": "white" } }, "shapedefaults": { "line": { "color": "#2a3f5f" } }, "ternary": { "aaxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" }, "baxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" }, "bgcolor": "#E5ECF6", "caxis": { "gridcolor": "white", "linecolor": "white", "ticks": "" } }, "title": { "x": 0.05 }, "xaxis": { "automargin": true, "gridcolor": "white", "linecolor": "white", "ticks": "", "title": { "standoff": 15 }, "zerolinecolor": "white", "zerolinewidth": 2 }, "yaxis": { "automargin": true, "gridcolor": "white", "linecolor": "white", "ticks": "", "title": { "standoff": 15 }, "zerolinecolor": "white", "zerolinewidth": 2 } } } } }, "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzkAAAFoCAYAAAB0XzViAAAQAElEQVR4AezdB5wcdf3/8c+Wy6WTTggtARJ670QhoFQTqqEjKCCISFVKUIyKoXdEEFD8I/ATREFKaAIWeu+B0CEhCQTSLneX2/Lf94Q5Jpu9u929LVNePPhmdme+853v9/md3Z3PfGfm4ln+QwABBBBAAAEEEEAAAQRCJBA3/kMAgQICzEIAAQQQQAABBBAIqgBBTlB7jnojgAAC9RBgmwgggAACCARAgCAnAJ1EFRFAAAEEEEDA3wLUDgEE/CVAkOOv/qA2CCCAAAIIIIAAAgiERaBu7SDIqRs9G0YAAQQQQAABBBBAAIFqCBDkVEOVMisnQEkIIIAAAggggAACCJQoQJBTIhjZEUAAAT8IUAcEEEAAAQQQ6FiAIKdjG5YggAACCCCAQLAEqC0CCCDgCBDkOAz8gwACCCCAAAIIIIBAWAWi1y6CnOj1OS1GAAEEEEAAAQQQQCDUAgQ5oe7eyjWOkhBAAAEEEEAAAQQQCIoAQU5Qeop6IoCAHwWoEwIIIIAAAgj4UIAgx4edQpUQQAABBBAItgC1RwABBOorQJBTX3+2jgACCCCAAAIIIBAVAdpZMwGCnJpRsyEEEEAAAQQQQAABBBCohQBBTi2UK7cNSkIAAQQQQAABBBBAAIEuBAhyugBiMQIIBEGAOiKAAAIIIIAAAl8LEOR8bcErBBBAAAEEwiVAaxBAAIGIChDkRLTjaTYCCCCAAAIIIBBVAdodfgGCnPD3MS1EAAEEEEAAAQQQQCBSAgQ5ZXU3KyGAAAIIIIAAAggggIBfBQhy/Noz1AuBIApQZwQQQAABBBBAwAcCBDk+6ASqgAACCCAQbgFahwACCCBQWwGCnNp6szUEEEAAAQQQQACBpQL8i0DVBAhyqkZLwQgggAACCCCAAAIIIFAPgWAHOfUQY5sIIIAAAggggAACCCDgawGCHF93D5VDoDwB1kIAAQQQQAABBKIsQJAT5d6n7QgggEC0BGgtAggggEBEBAhyItLRNBMBBBBAAAEEECgswFwEwidAkBO+PqVFCCCAAAIIIIAAAghEWqAiQU6kBWk8AggggAACCCCAAAII+EqAIMdX3UFlQiZAcxBAAAEEEEAAAQTqIECQUwd0NokAAghEW4DWI4AAAgggUF0Bgpzq+lI6AggggAACCCBQnAC5EECgYgIEORWjpCAEEEAAAQQQQAABBBCotEA55RHklKPGOggggAACCCCAAAIIIOBbAYIc33YNFaucACUhgAACCCCAAAIIREmAICdKvU1bEUAAAa8ArxFAAAEEEAipAEFOSDuWZiGAAAIIIIBAeQKshQACwRcgyAl+H9ICBBBAAAEEEEAAAQSqLRCo8glyAtVdVBYBBBBAAAEEEEAAAQS6EiDI6UqI5ZUToCQEEEAAAQQQQAABBGogQJBTA2Q2gQACCHQmwDIEEEAAAQQQqKwAQU5lPSkNAQQQQAABBCojQCkIIIBA2QIEOWXTsSICCCCAAAIIIIAAArUWYHvFCBDkFKNEHgQQQAABBBBAAAEEEAiMAEFOYLqqchWlJAQQQAABBBBAAAEEwixAkBPm3qVtCCBQigB5EUAAAQQQQCAkAgQ5IelImoEAAggggEB1BCgVAQQQCJ4AQU7w+owaI4AAAggggAACCNRbgO37WoAgx9fdQ+UQQAABBBBAAAEEEECgVAGCnFLFKpefkhBAAAEEEEAAAQQQQKAKAgQ5VUClSAQQ6I4A6yKAAAIIIIAAAt0TIMjpnh9rI4AAAgggUBsBtoIAAgggULQAQU7RVGREAAEEEEAAAQQQ8JsA9UGgkABBTiEV5iGAAAIIIIAAAggggEBgBQhyLLB9R8URQAABBBBAAAEEEECggABBTgEUZiGAgJmBgAACCCCAAAIIBFSAICegHUe1EUAAAQTqI8BWEUAAAQT8L0CQ4/8+ooYIIIAAAggggIDfBagfAr4SIMjxVXdQGQQQQAABBBBAAAEEEOiugH+CnO62hPURQAABBBBAAAEEEEAAgZwAQU4Ogf8R8LMAdUMAAQQQQAABBBAoTYAgpzQvciOAAAII+EOAWiCAAAIIINChAEFOhzQsQAABBBBAAAEEgiZAfRFAQAIEOVIgIYAAAggggAACCCCAQGgElgtyQtMyGoIAAggggAACCCCAAAKRFCDIiWS30+gyBFgFAQQQQAABBBBAICACBDkB6SiqiQACCPhTgFohgAACCCDgPwGCHP/1CTVCAAEEEEAAgaALUH8EEKirAEFOXfnZOAIIIIAAAggggAAC0RGoVUsJcmolzXYQQAABBBBAAAEEEECgJgIEOTVhZiOVE6AkBBBAAAEEEEAAAQQ6FyDI6dyHpQgggEAwBKglAggggAACCLQLEOS0U/ACAQQQQAABBMImQHsQQCCaAgQ50ex3Wo0AAggggAACCCAQXYHQt5wgJ/RdTAMRQAABBBBAAAEEEIiWAEFOtPq7cq2lJAQQQAABBBBAAAEEfCpAkOPTjqFaCCAQTAFqjQACCCCAAAL1FyDIqX8fUAMEEEAAAQTCLkD7EEAAgZoKEOTUlJuNIYAAAggggAACCCDgCjCtlgBBTrVkKRcBBBBAAAEEEEAAAQTqIkCQUxf2ym2UkhBAAAEEEEAAAQQQQGBZAYKcZT14hwAC4RCgFQgggAACCCAQYQGCnAh3Pk1HAAEEEIiaAO1FAAEEoiFAkBONfqaVCCCAAAIIIIAAAh0JMD90AgQ5oetSGoQAAggggAACCCCAQLQFCHIq0/+UggACCCCAAAIIIIAAAj4RIMjxSUdQDQTCKUCrEEAAAQQQQACB2gsQ5NTenC0igAACCERdgPYjgAACCFRVgCCnqrwUjgACCCCAAAIIIFCsAPkQqJQAQU6lJCkHAQQQQAABBBBAAAEEfCEQsiDHF6ZUAgEEEEAAAQQQQAABBOooQJBTR3w2jUDNBNgQAggggAACCCAQIQGCnAh1Nk1FAAEEEFhWgHcIIIAAAuEUIMgJZ7/SKgQQQAABBBBAoFwB1kMg8AIEOYHvQhqAAAIIIIAAAggggAACXoHqBDneLfAaAQQQQAABBBBAAAEEEKihAEFODbHZFAIIIIAAAggggAACCFRfgCCn+sZsAQEEEAilQPqtV63t2f/aksfus9Z7b7OW226w5hsvt8W/+601XXCGLZr8E1v4syNswXH72bwjdrWXzv2jHX9am/3sl232y/NSdt7lKbv82pRd9//SdvPtabvzvrQ98p+Mvfxa1j6ZmbXWJaFko1EIIIAAAjUQiNdgG2wCAQQQQCDIAtmsZT5535Y8eq8tvu4iW3ja923egTvYwl/8yJouPNMWXz3Fmv98hbX87U/Wet/ttuTfU63tuf9Z6o0XLf3hO5b5fLbZ4iaLZTPW0mr25TyzGZ9m7Z33svbqG1l7+vmMPfq/jN3zQMZuuSNtV16Xssnnp+zHP2uzn5ze5ry+8g8pu+VvaXvs8Yx9misuyJzUHYHiBciJAALlCsTLXZH1EEAAAQTCKZBdMM8ZoWm59Vpb9OsTnVGYBaccZot/f64teehOS38w3SyTrknjm1vMGdV5+fWsPfLfjP3ltrT9YkqbnXxWm13zp5Qzb2YuYKpJZdgIAggggIA/BIqoBUFOEUhkQQABBMIukHrtBVt8zXk2/0f72vyjxjsjNC3/uMlSrz1v1rzYd81fuMjsuZeyzujO2eel7MQz2+x3N6ScEaEm/1XXd35UCAEEEAi7AEFO2HuY9hUSYB4CCOQEdClZ81+utvnH7pMbsTnBljxyj2XnzsktCd7/CmxefCXr3NtzSm6UR/f6PPVcxlpbg9cWaowAAggg0H0BgpzuG1ICAgggEBiB7OezreXOm2zBqYc5DwVo/ectlv3is6/qH45JOmPOvT7X35S2k3IBzzU3pk0BUKo2V9iFA5FWIIAAAgEXIMgJeAdSfQQQQKArgeyihc69NIvOPs7mH7eftdxyrWU+fr+r1UKxvK3N7LkXM86lbCed2WY3/y1tc78IRdNoRC0F2BYCCAROgCAncF1GhRFAAIHiBLLzvzQ90nn+MXs5T0VLTXuluBVDmktPdnv0vxk749dtdu2f084DDULaVJqFAAII1ETAzxshyPFz71A3BBBAoAwBPR2t+abf2fwff9d5pLO1LSmjlPCuks2aPftCxnk09SVXp+yNt3IzwttcWoYAAghEUoAgJ5Ld7pdGUw8EEKikQLZpoXMp2vzjJ1rr3beaLeGu+658FeAo0Pn1BSkn8FEA1NU6LEcAAQQQ8L8AQY7/+4gaIoBA1ARKbe/iJmu57Xpz7re58yazluZSS4h8/o9mZJ1L2M4+t83eficTeQ8AEEAAgaALEOQEvQepPwIIRFegZbG13HGj6bK0lr/d6Mu/ZxO0zvl0ttkFV6btuv+XtvkLglb78NeXFiKAAALFChDkFCtFPgQQQMBHAnr08/zjvmstf73edJmaj6oWiqo8/XzGJp3TZg8+mjE9kjoUjaIRCCAQVgHaVUCAIKcACrMQQAABvwqkP/3YFpz6PdMf8cwuYqihmv2kPyR6251p4xK2aipTNgIIIFAdAYKc6rgGq1RqiwACgRBovf8OW/izwy3z8XuBqG9YKjl7ztJL2G66LW36uzthaRftQAABBMIsQJAT5t6lbQgg0C0Bv6ycWTjfFv32VGv+46VmS3gcdL365d+PZ2zyBW0241MeOV2vPmC7CCCAQLECBDnFSpEPAQQQqINA20tP2cKTD7HUy0/XYetsMl9Aozq/uTBlD/870k9gy2fhPQIIIOA7AYIc33UJFUIAAQRyArkRm+brL7amKT81/XHP3Bz+94lAKm32f39Pm/6+zsJFjOr4pFuoBgI+EKAKfhIgyPFTb1AXBBBAICeQ/uhdW/DT71nrg//IveN/vwroD4n+8ryUvT6NUR2/9hH1QgCB6AoQ5Pio76kKAghEXCCbNT0aeuEZR1pm1icRxwhG8xcsNLv092l75D+54Z1gVJlaIoAAApEQIMiJRDfTSAQCLRCJymczaWu6fLLzaGhLpSLR5jA18pY7Mnb7XWnLxalhahZtQQABBAIrQJAT2K6j4gggEBaBbDpliy/+hbU98S/jv+AKPPBIxv7w55Sl07W6Tye4VtQcAQQQqLYAQU61hSkfAQQQ6ERAAU7TBWda27P/6SQXi4Ii8OyLWbv82rS1thLoBKXPqGcIBWgSAjkBgpwcAv8jgAAC9RDIplKmACf14pP12DzbrJKAHkhw3uUpW9REoFMlYopFAAEEuhQgyFmeiDkIIIBA1QWybUus6dyfGgFO1anrsoGPZ5hNuTRl8+bXZfNsFAEEEIi8AEFO5HcBABAoVoB8lRLILmm1pt+eaqlXn6tUkZTjQ4E5n5ldeGUbIzo+7BuqhAAC4RcgyAl/H9NCBBDwkUC2ZbEt+vWJlnrjRR/Viqp0S6CTlWfnAp2LrkpZSwuXrnXCxCIEEECg4gIEORUnpUAEEECgIqMg7AAAEABJREFUsEB2cZMt+tUJln77tcIZmBtKgU9mLv1bOm1tBDqh7GAa1aEACxCopwBBTj312TYCCERGIJvJWNNFkyz97rTItJmGfi3w7gdZu/K6tKVSBDpfq/AKAQQQqJ6Aj4Oc6jWakhFAAIFaC7T85XeWeu35Wm+W7flIQE9du/qPactkCHR81C1UBQEEQipAkBPSjqVZIRagaYETWPLEv6z1nr8Grt5UuPICr7yetetvSls2S6BTeV1KRAABBL4WIMj52oJXCCCAQMUFUu+/bYuvOqfi5VLg8gJBmfPMC1mb+jBBTlD6i3oigEAwBQhygtlv1BoBBAIgkFnwpTWd+zOzVFsAaksVaynwj3vT9s57mVpukm1FV4CWIxBJAYKcSHY7jUYAgWoLZFOpXIBzmmXnza32pig/gAK6Wk335yxqYkQngN1HlRFAIAACXQc5AWgEVUQAAQT8JtD8p8ss/e6bfqsW9fGRwIKFZr//I/fn+KhLqAoCCIRIgCAnRJ1JU2orwNYQ6EhgyWNTbclDd3a0mPkItAu89U7W7nmA0Zx2EF4ggAACFRIgyKkQJMUggAACEki984YtvvY8vYxqot0lCvzz/rRNm879OSWykR0BBBDoVIAgp1MeFiKAAALFC2SXtFrTxT83S6eLX4mckRfQ/Tl/+HPaFjdHniLkADQPAQRqKUCQU0tttoUAAqEWaL3jRsvOnRPqNtK46gjo/pzb7iQ4ro4upSKAgK8FqlQ5gpwqwVIsAghESyA9e4a1/POWaDWa1lZU4H9PZezd97lsraKoFIYAApEVIMiJbNeHpuE0BAFfCDRfcx6XqfmiJ4JdiT/ekra2FA8iCHYvUnsEEPCDAEGOH3qBOiCAQKAF2h5/2FKvv+izNlCdIArMnmN230OM5gSx76gzAgj4S4Agx1/9QW0QQCBgAtnmxbb4xisCVmuq62eBex/M2KzZjOZUrY8oGAEEIiFAkBOJbqaRCCBQLYGWv15v2flfVKt4yo2gQCY3kKPL1iLYdJqMAAJ1FAjbpglywtajtAcBBGomkP7oXWud+reabY8NRUfgvQ+y9tyLuWgnOk2mpQgggEBFBQhyKsoZ5cJoOwLRE1h89RSzLAei0ev52rT4rqlpy2S4bK022mwFAQTCJkCQE7YepT0IIFATgSWP3mvp997qelvkQKBMgU9nmz37IkFOmXyshgACERcgyIn4DkDzEUCgDIElrdZ88+/LWDFiq/TuY7FBQyw+YjVLrLG2JYYOt2FDzXr0iJhDN5p79/3h/QOh3WBhVQQQQKBLAYKcLonIgAACCCwrsOTfUy27YN6yMyP+Lr7iCOsxbg/r/aMzrf8V/2cDbvufDbjxAVvhmjut/2W3WL/zbrANjh5vU37eYFdf2GBXnt9g55zVYD89PmlHHZawiXslbPONY9ZIALTMnjRrjtmTz3JJ5DIovEEg3AK0rkIC8QqVQzEIIIBAZARaeNiA09fJDTa33seeYStcd7f1v/I2633cJOux43csPnwVZ3ln//TqaTZ8mNk6o2O2zRZx23WnuP3oB0m7/LwGO/GYpG2/Xdz69+ushOgs+2duNCed5rK16PQ4LUUAgUoIEORUQtFPZVAXBBCoqkDq1ecs88kHVd2GnwtPjFnfen3/JCew6Xv25dZjp/EWW2FgxaqcTJhtuF7MvndAwi45p8HOOCnpBEADVqjYJgJX0Gefmz31HEFO4DqOCiOAQF0FCHLqys/GEUCgVgKV2k5UHxmd3HQb63fRn63fOdda4+7frWhg01nfrDUq5lzKdtGvG+zg/RLWr29nucO77J9TuTcnvL1LyxBAoBoCBDnVUKVMBBAIpUBm7hxre/7xULato0Yl1t7Q+p57vfU98yJLrLZmR9lqMn+n7eN2/i8bbN/xCevdqyab9M1G5n5p9sIrVRnN8U0bqQgCCCBQSQGCnEpqUhYCCIRaoPW+282y0TjQjA9byfpOusj6/eb3llxzHd/0q57MtsfOuWBncoON3zVujY2+qVrVK/LfJxnNqToyG0CgXYAXQRcgyAl6D1J/BBCoiUC2tcVaH7qz6G1NT8Xs8Hk9beKXvZz0i4WNttjzkKw7mpPOfHe5ppcvanDKn5vL96N5jc5yTfXeWZD7R+u5+XJvq/J/j2/vaf0vvsmSm2xTlfIrUageXLD3Hgk7Z1LSRq4Wq0SRvi/jtTezNm9+NIJs33cGFUQAAd8LEORUqYsoFgEEwiXQ9u/7zVqai27UK20JO653m90+sNlJQ2IZu27x0iDGLWSdZMb+vMLS5cp3Yt82Z9FjrUn7dmPaWU9TvdcCBTsvpRJ2aK5cva90ivXtb31yoze9f3iaWWPPShdflfIGDojZmScnbedx4f850yDi/54myKnKjkShCCAQOoHw/yqErstoEAKBFghs5Vum3l5S3ffrlbKte3x9edEqiax9no0vM5rTUYGfpGOm/Fquqd7rtYKdTZJpG1yFb+7kxltbv0tvtgYfj97IoFBK5DwO2Cdhxx+VsLBfvvY/Llkz/kMAAQSKEcj9NBSTjTwIIIBAdAWcx0bP+LBsAF2mphEYjeb09nzrTkvF7fD5Sy9n816C5g1sFODovTuKM64xVXY9Olqxx7cmWN+zLrZ4BR8F3dG2qjl/kw3jNvn0pK26cjW3Uq2yiyv38y/Mpk1nNKc4LXIhgECUBTw/t1FmoO0IIIBAxwJL/vtgxwu7WKLgRYGMsh3tucxMIz26RE3pmhWaTQGP7rdRPgUyD7cmnHtyNNV7dxTni0zM3Ht98u/z0bqlph677mu9jzm91NV8m3/o4Jj98rQG54+J+raS3azY/57KdLMEVkcgQAJUFYEyBQhyyoRjNQQQiI5A2/NPlN1Y3WejQEaXmZ3b1FjwcjVdfqZ7bzRqow3p/e8HtDr35GiqeRoJUrBzX0uy/V4fjQxNbU1qcVmpccKB1vvIU8pa1+8r6Y+JbrVZOB9I8NxLGVu8mNEcv++D1A8BBOorEPYgp766bB0BBAIvkPnoPcsunNftdmzUkDYdlxb/6IKvN+mO4uhPw+i+HneJLmNzX5c6bfjmLtbrsONLXS1Q+X9waMLWGBm+QCeVMnv5dYKcQO2MVBYBBGouQJBTc3I2iIAfBKhDsQJtLz1VbNZl8v2hKWl6jLQ7UyMwvXPH2wpUdI/OTYuTpqmW634bXZa2TY/lL0PSMncUR/fzaPRG6yi5Iz96XUpKrLmO9f7RpFJWCWTeZCJmJx6TsGFDAln9Tiv9+jSCnE6BWIgAApEXIMiJ/C4AAAIIdCbQ9vIznS3ucNnGDVmbtLCnc1+N/gaORmDO7NNqClSUnHtrvnrowLG5qS5X8z6NzS34L4sbbHxjqv2Janp89I25AMktc/fcMjdvMdPYgMHW58wLLZYs/zK3Yrbjlzx9cpHlKcclrbeiy2IqFZA8r09bPiAOSNWpJgIIIFATAYKcmjCzEQQQCKJAtm2Jpd58uayqK2DRvThu+k2/pQGOW5h7r467XA8icJd5p8qnstx53vt18st083Q4TTZY30kXWbz/wA6zhHHBkMFLR3T0qOmwtG/hIrOPZzCaE5b+DEI7qCMCQRMgyAlaj1FfBBComUD61efNUkv/QGfNNlrFDfU68GhLjBxdxS34t+g1R8Vt953D9ZP3Opes+XeHo2YIIFB3gRp949e9nVQAAQQQKFmg7eWnS17HryvEh42wHnvs79fq1aReu38rbv361mRTNdnI61yyVhNnNoIAAsEUIMgJZr9R67AI0A5fC5R7P44fG9Xr2DMicx9OR/6NjTHbd3yio8WBmz/93ay1pbhkLXAdR4URQKAmAgQ5NWFmIwggEDSBzNw5lpn5UdCqXbC+DVvvYA0bbFZwmV9nVqte39w2biNWilWr+JqWm0qbvTW9pptkYwgggEBgBAhyAtNVVBQBBGopkHqh/D8AWst6drmthh7W6/ATuswWpQwH7xcPTXM//JiRnNB0ZnENIRcCCBQpEJ5v+iIbTDYEEECgGIGwXKrW+J39LT5kxWKaHJk864yO24brhWM055OZBDmR2XFpKAIIdCKw/CKCnOVNmIMAAghY+v23g68QT1hjxB82YB38d8A+CYuFIM75ZCZ/L6eDLmY2AghEXIAgJ+I7AM1fKsC/CHgFspm0ZT6b5Z0VyNcN2+xg8QGDAln3ald6+LCYbbVZ8H8CP51dbSnKRwABBIIpEPxv+GC6U2sEEPCxQPazcBw59tzjgO4qh3r9rTYPx0/gR59wyVqod1QahwACZQmE4xu+rKazEgIIIFBYIDN7RuEFAZqbWH0tS4xZP0A1rn1V11/brLGx9tut9BY/4b6cSpMWUR5ZEEDA7wIEOX7vIeqHAAI1F8jM+bTm26z0Bnvs+J1KFxm68pLJmG22UfB/Bj8myAndvkmDEAisgI8qHvxvdx9hUhUEEAiHQHrOzMA3pGGTrQPfhlo0YItNg/8zOIMgpxa7CttAAIGACQT/2z1g4FS3UwEWIuALgUzAg5zYoKEWH7GaLyz9XokwXLL2+VzuyfH7fkb9EECg9gIEObU3Z4sIIOBzgcxsv43klAbWsPnY0laIcO4wXLK2qCnCHUjTEUAAgQ4ECHI6gGE2AghEVyAT8HtyGrhUraSdd/Qawf6DOYubzTKZiI7mlNTTZEYAgSgJEOREqbdpKwIIdCmQbVti2YXzu8zn5wzJNdf1c/V8V7dhQ4Md5Ai0aXHw26B2kBBAoDIClGJGkMNegAACCHgEMjM/8rwL4Mtkg8UGDQlgxetX5WFDYvXbeIW2zCVrFYKkGAQQCI0AQU5ourKSDaEsBKIrEPRL1eLDV45u55XZ8oEDshYv89dwrz0Sds0lDXbqj5PtW992y7hddUGDXX/50uRdpky/PjPZvuzwAxOa5SSVddmUBtP6zowS/lnUxOVqJXCRFQEEIiBQ5td6BGRoIgIIRFIg27So43YHYEliRYKcUrspFovZikNLXctMQcnu345b8us4xdYZHbN9JyTs/Q+zdtSJbXb3AxkbvWbM3GBG6zQ0xOyiq1LOsnXGxJ11LPff5hvH7MVXMvbks5ncu9L+ZySnNC9yI4BA+AUIcsLfx7QQAQRKESj3lH4p26hi3jhBTlm6Q0u8ZE2jLTt+I26P/S9jX3pu4dp687j16mn2xDNLA5W77kvbnM+ytuaopZfEDehv1taWtWnTs878ZG4AaOCAmBMwKfh5+vml65XaiEWL6juSU2p9yY8AAghUW4Agp9rClI8AAoESiAU9yOFytbL2t1Luy1GAc8jEhDPq8tKrXQcl8xeY9eoVc0Zs5uVeK5jRiI8eeJBKmX05L2saxZn2dsYJfsppQOuSctZiHQQQqLIAxddRgCCnjvhsGgEEfCgQ8CAn1ruPD1H9XyUFHMXUUsGJLkd79oWM/fn/0sut8va7S0dUxqy5dOQmP4NGdjSS89PjkzZh17gpsFl7dNwU+GjdC3+99LqmP2MAABAASURBVD4e3eejS9usyP9ihTdX5NpkQwABBMInQJDj5z6lbgggUHuBgAc5pqGB2qsFf4tLY5Mu26FLy3Q52je3jTsPD1CwMnAFs3XHxJyHD+h+GgVA7nI9fEDLmpuXXqKmDZx9bsq5X0f37ChQckdxFBgpn+Y/+UzGvrHN1/fraL3OUpwgpzMeliGAQAQFCHIi2Ok0GYGgC1S1/kEPctKpqvKEtfDFzcW1TEHM8ae1tQcpeoCA7sl58+2sXfy7pfYKXBSouGnmrKy9+37hKEqjNRrF0b04QwZ/HanosrZS4tW45+EHxbWEXAgggEC4BQhywt2/tA4BBEoViAX7azGbXv4SqlIJoph/cW6kpRrt1uOjdT+OgphC5bujOHoQwedzvw6E9IACPZSg0DqF5iWW7raFFjEPAQQQiKQAX4uR7HYajQACHQoEfSSnlNP/HSJEb0GxIznFyCiw0WVqSiv0N/vZ2W2mICZ/Xfex0hr50TJ3qvW23Spu/3uq+AcRBH23VftJCFRPgJKjKECQE8Vep80IINChQCzg1/1k5s3tsG0s6Fig3JEcBS8KYtxL1bQFvXYvVdP9N5pXKCmoyV+u91r32FPaTA8pKLReoXkEOYVUmIcAAlEWIMgpovfJggACERII+NFiZtaMCHVW5ZraXOQ9OZXbYmVL4nK1ynpSGgIIBF+AICf4fUgLEKiXQDi3G/Bn8aZnfRzOfqlyqyp5uVqVq1qw+Hji64cWFMzATAQQQCBiAgQ5EetwmosAAl0IJBu6yODvxZlPP/F3BX1Yu2w2a5/O/vqm/+5XsfYl9O5V+22yRQQQQMDPAgQ5fu4d6oYAAjUXiA8aWvNtVnSDbUss+8VnFS0y7IV98JHZkiXBbuXgQYzkBLsHI1J7molADQUIcmqIzaYQQMD/ArFhw80C/hjp1DtvGP8VL6C/cVN8bn/mHDww2CNR/lSlVgggEGSBIAU5QXam7gggEBABPV0tvtIqAalt4Wq2vfR04QXMLSgwbXqm4PygzNRjqhPckxOU7qKeCCBQIwGCnBpBsxkEqidAyZUWSAwPeJDz4lOVJglteel01t5+J9ijIIMGcqlaaHdQGoYAAmULEOSUTceKCCAQVoGgj+Rk586xzMyPwto9xberiJzvfWiWSheR0cdZuB/Hx51D1RBAoG4CBDl1o2fDCCDgV4HESqv5tWpF16vtJUZzisGa9nawR3HUxsED9S8JgeIFyIlAFAQIcqLQy7QRAQRKEogH/J4cNXbJY1M1IXUh8Nq0YN+Po+YxkiMFEgIIILCsQBlBzrIF8A4BBBAIm0B8pVUD36T0B9Mt/fZrgW9HNRvw4cdZe/f94I/krDKCe3KquZ9QNgIIBFOAICeY/Uat/ShAnUIjEB+yolmyIfDtabnvtsC3oZoNuPO+gN+Mk8OJ5eKbkcGPyXMt4X8EEECgsgLxyhZHaQgggEA4BOIhuGSt7clHLTPvi7p3iB8rMOPTrL36RvBHcVZeKWY9euQiHT8iUycEEECgjgLxOm6bTSOAAAK+FUiE4JI1y2atldGcgvvY3fcHfxRHDRu1OgGOHAKaqDYCCFRRgCCnirgUjQACwRWIB/xv5bjyrfffYZkvP3ffMs0JzP4sa8+9FPxRnFxTbA2CHDGQEEAgVAKVaQxBTmUcKQUBBEImkByzQTha1NJszf/vqnC0pUKtuPfB4D9RzaVgJMeVYIoAAggsK0CQs6wH70IgQBMQqIRAcqMtLQwPH5BF2+MPW2raK3oZ+fTOexl78tlwBDkNDWYjhke+SwFAAAEECgrEC85lJgIIIBBxgVjPXpZcf9PQKCy+5jzLplKhaU85DVnUlLWr/5jWrUrlrO67dUatFrN4nHtyjP8QQACBAgIEOQVQmIUAAghIoGHzsZqEImVmfmRLIv4Qgt/nApwFC0PRnU4j1h3DT7gD0e1/KAABBMIowDdkGHuVNiGAQEUEGrb8ZkXK8UMhH48aaTuuZvZk0yyL4n/3/ytjb70TjocNuP23yYaM4rgWTBFAoAoCAS+SICfgHUj1EUCgegLxwcMsvsqo6m2gRiXfPWF32/bAHezJ9HzbY/q9Nr1lXo227I/NvP9hxv5+TzgeGe2KrtDfbNWVCXJcD6YIIIBAvgBBTr4I7yslQDkIhEKgYfPtAtuORf362fGnHGPfW2+QNWWW3o8zL73Edp5+j32Zag1su0qpeEtL1n53Q9oy4XjWQHvTt9iEn+92DF4ggAACBQT4liyAwiwEEEDAFWio+H05bsnVnb6y2ab2jZ8caDc3LFpuQx8uWWh7vHOvtWbCNbqR39C2NnMCnHnz85cE//2mG/HzHfxepAUIIFBNAb4lq6lL2QggEHiBxNobmvXqE5h2ZGMxu/IHh9i3d9nAPkwv7rDeTzXNtoPff9iy2XDdp+I2uKkpa+ddnrI33w5f+3o2mo1Z021piKY0BQEEEKigAEFOBTEpCgEEwicQi8WsR0BGc2aPGGH7nH6MnT3UrC3b9fVZf5/3nk2a+XToOu2LeVk75+KUffhx+AIcddbGG8R5dLQgSAhERIBmlicQL2811kIAAQSiI9AQgPtyHt55J9vm8F3t39nSnpF83qwX7biP/mNLigiKgtDjs2bnApyLUvbZ3CDUtrw6bspT1cqDYy0EEIiUAEFO6LubBiKAQHcFklttb7E+/bpbTFXWb+3Z08444SibuNlKNi+zpKxt/P6z123rN++wj5aUFiCVtbEqrvTBRxn77SUpC9Pfwsnnamw025ggJ5+F9wgggMByAgQ5y5EwAwEEIiFQQiNjDT2scdd9S1ijNlnfWn9d2+mUw+3aXs3d3uBLzZ/bhq/fZlPnf9TtsmpdQCptdtfUtJ13WdqaW2q99dpub9zYuDUkeXS08R8CCCDQhQBBThdALEYAAQQk0GO3/cziCb30RfrzwRNt3IQt7Y308k9PK7eCC3IjQd955147a8bTlrFg3M/y9rsZO3tKm919f8YU7JTb9iCsF4uZfWuH6v9sB8GCOiKAAAJdCfBt2ZUQyxFAAIGcQHzAIOux/S65V/X9f/6gQXbw6cfaSav2sJZsbgijwtVRaDNl1gu20et/tb/Pe8+3oU7TYrMbb0nbBVekbc7nFUbwaXEbrR+zQQNykY5P60e1EAi5AM0LmABBTsA6jOoigED9BBr3PKR+G89t+Yntx9o2x+5jU63698683vKl7ffuA7b5G7fbvfM/zG3dP/8/+WzGzjqnzf73dMY/lapBTb61vX9GEmvQXDaBAAIIdEuAIKdbfCWsTFYEEAi8QGKVkZbcaMuatyOVSNpvf3SETdh2pM1K1/amkxebP7fx79xn20y7wx5e8EnN2+5ucEmbOUHNuZel7Ia/pG1Rk7skGtPhw8zWW5tRnGj0Nq1EAIFKCBDkVEKRMhBAoGyBoK3Yc/wBNa3yx6NG2q6nH2kX9W+zeo5bPN00x3aefreNee0WO+2TJ+3JpllVv5QtmzWbNj1rf7w5bSdPanMuT3v3/dzMmvaAPza2E6M4/ugIaoEAAoERIMgJTFdRUQQQ8INAcpNtLL7y6jWpyh37TrCxB46zFyr4cIHuVnx663y7cPZLtt20f9hKL99oP/zwMbtv/ocV/Ts7n881u/PetJ3xqza76KqUPfFMxlqXdLfmgVu/vcJ6bPTYrRnFaQfhBQIIIFCEAEFOEUhkQQABBLwCPccf6H1b8deL+va1H/3sR3bU6P62MNNW8fIrVeDsVLNd9/mb9p137rNBL91gE997wKY/lLGZL2VtwYyspTsJTFpazT74KGvPvJCxex7ImC5BO/eylJ18Vpud8es2u+fBjM39slI1DXY5u4yLW2OPWLAbQe0RqJgABSFQnABBTnFO5EIAAQTaBRq239Wq9cdBX9lsU/vGCQfZ/8UXtG8vCC+aMin7eGarvfzXtD2RG3158Jcp+8dxbU6692dt9tDklD16Xsrun9Rmd/64zX6TC2TOuThlf/hz2u68L216mIAuRVtYuSdiB4Gtyzr27WO227f5qe4SigwIIIBAnkDkvjnz2s9bBBBAoGSBmP446B4TS16vsxUy8ZhdeeSh9q2d17cP04s7y+rbZd+Yt9pyddNoTnNuRGb+J1mb+07WFs0xS7Wa9emxXFZmFBDY+zsJRnGM/xBAAIHSBQhySjdjDQTCKECbShRonHCgxVYYVOJahbPPHjHCJpx+jJ09JGupqt/OX7gOlZi79txhRRfTt4HLr7rCGjrEbPttcerKieUIIIBAIQGCnEIqzEMAAQS6EIj17G29v39SF7m6XvzwzjvZNofvak9kFnad2ec5hsxeoega9kwWnbXOGeu3+f33Tlo8N8JXvxqwZQQQQCC4AgQ5we07ao4AAnUWaNhuJ0tuvFVZtWjt1dNOP/Fom7jZSjYv08kd+mWVXp+VYh83FL3hXomis0Yy4xqrx2zTDRnFiWTnB6XR1BMBnwsQ5Pi8g6geAgj4W6D3sWeYNZR2g8lb669rO518uP2hZzDvvSnUIxu2DLNUc/EH5aWJFdpiuOcduj9RYLh7mNYhgEC1BeLV3kAH5TMbAQQQCIVAfPAw63XAUUW35cZD9rdxE7a0N3z0t2+KrnwnGbdfsHonS5dfRJCzvIk75xvbxG21VYoPGN31mCKAAAIIfC1AkPO1Ba8Q8IEAVQiiQI/xB1h81VGdVn3+oEF28OnH2smrNFhLNt1p3iAuXP+LFUuqdjJTUvbIZF6hv9lB+/LTHJkOp6EIIFA1Ab5Jq0ZLwQggEBWBWDxhvX98VofNfWL7sbbNsfvYVAv+wwU6auSKc4p/6IDKiJca5GilCKQffi9hjY2M4hj/IYAAAt0UIMjpJiCrI4AAAhJIrrGONe62n162p1Qyab897vs2YduRNivd0j4/jC8SM3qW1Kx4qqTskcg8bmzc1h7Nz3IkOruCjaQoBBAoLMC3aWEX5iKAAAIlC/Q8+Bhz/3bOx6NG2q6nHWkX9Vti1Ri06P/RF7bNZQ9b4/zm9nomWtpsi2v+bd8+4+/tadhrM9qXK+/Y86Y6yzTVe3fhqEem2Qa3PuO+LWk6PNXH2uaVOPrQVtImQp950ACzA/bhJzn0HU0DEUCgVgLGN2rNqNkQAgiEXcD52zlHnmx/23dPG3vgOHuhCg8XUGCiAGWrqx8zBTVe02RrylpW6GWPTp5gD5+3r71y6Na23t9eMAVEyjfi+Q9t5lajnGWa6r3mq8zBb8+26XtsqLclp+/MG1PyOpnWbMnrhHmFHx6RsAb+QGqYu5i2IYBAjQUIcmoMzuYCIkA1EShToGGbHe2ZbbeyhZnqDFW05oKYx8/Y3Z45bpyley77d2m07LWDtmqfP3/VQdbWM2k9Fywd7ekze4E1Detn+k9TvddrBTtzx6xoWl/vS00bfDm81FUs01ryKqFd4ds7xG2tUfwch7aDaRgCCNRFgG/VurCzUQQQCLPA71b7pm3aa0jdm6jooUz1AAAQAElEQVQRmpjFrKV/L6cuTSv2tz5zlj78QFO9Vx6N4szcvLhHQDsF5f2zypwBeXO6ftsWnj8R1HVjO8mx6spm++3JT3EnRCxCAAEEyhLgm7UsNlZCAAEEOhZoiCXsnrX2sKHJpcFFxzmrt0SXso2571WbsdVIW7DaIGdDCmRGPPO+c0+Ops775z80jeIo2Bk3+W5nme7r0frOSkX802Nm6e3MpMwa4kUUHuIsffuYnXhM0hqSsRC3MjJNo6EIIOAzgYj/xPisN6gOAgiERmBEjz72zzV3t3huJKXWjVKAsumNT5juz3l/p3XaN6/L0XSpm+7X0VQL3FGc1R5/x9747mbO/Tpab7Un3tXiLtMKqUZb8nm8y3yFMvQv7YFshYoI7LxEjuzEYxI2YAUCnMB2IhVHAIEiBOqXJfc1W7+Ns2UEEEAgzALb9F3RLll1u5o20Rvg6P6czjbu3ouTakxaT89T2nQZmxX53x4LRheZc/ls/RuXnxeVOYcflLBRq/MTHJX+pp0IIFB7Ab5ha2/OFksQICsCQRc4cdhGdtDAtWrSDDfA0eVnXQU4ujzNHcXRAww0euNW0n0ggfu+s+lm80Z0trjTZb0j+jSxHcbGbbut4p3asBABBBBAoHsCfMt2z4+1EUAAgS4F/jRyJ9ug59L7YrrM3EUGBSdjz5tqeoR031kL7JvnTm3/+zZ6mIDmrfngG869Ne7fyyn0929G3/eqffSNtdqfqKbHR4++5xXTOhrV+Wi7NbuoydLFq3xWfrt6JZeWEaV/R64Ws0O+y09vlPqctiKAQH0E+KatjztbRQCBCAk0xhP24JjxNrLH0sc3d6fp3vtqdG+NkjtqowcMPPbV38jRfDe5y73b1bw5G6zcPstb7nPH7tD+GOr2DB286PNp7w6WdD27d6LrPGHKMWTw0gcNxOPch1O9fqVkBBBAYKkAQc5SB/5FAAEEqiqwUkMf+8/ae9vwZPlBQVUrWEbhjemEtc4q/2ekR4SO9QfnBrxOOyFp/fqWAc0qCCCAQHcFIrh+PIJtpskIIIBAXQRW7dHX/rP2XjYkGY7Hiu2+aC2zbPmUPcpfNVBrDs4FOKefmLRBAyIU1QWqh6gsAgiEUYAgJ4y9Wvk2USICCFRIYHTPAfafMXvbwETwHy22+ZdfX+5WDk+yGwFSOdurxzqDB5kR4NRDnm0igEDUBQhyor4H0H4EEOiGQHmrrttroD02Zi/rHw/2WMaozweXB/DVWon0Vy9COiHACWnH0iwEEAiEAEFOILqJSiKAQNgENuo92B4Zs6f1igX3EWP9ZnXvBpN4iIOcyAc4YfvA0h4EEAicAEFO4LqMCiOAQFgENu8z1P679t42LNkrcE2KZ83aZnbz8WipwDW7qAqvurLZWadwD05RWGRCIGICNLd2AgQ5tbNmSwgggMByApvnAp0X1v2urdtz4HLL/Dzj24tGWbabIzGZVj+3sLy6bbJBzCadnLT+/XjIQHmCrIUAAghURoAgpzKONSqFzSCAQBgFVu7R155ZZz8b13dEYJq39bzVul3XTEtuOKjbpfingPG7xu34o5PW0ECA459eoSYIIBBVAYKcqPY87UYgTAIhaEvfRIM9PGaCHTVk3UC0Zs25Q7pdz1RLt4vwRQGJhNmxRyRs7z1yL3xRIyqBAAIIIECQwz6AAAII+EQgEYvbdauPs4tX2c78/uU8YFa/bquldblawAdz9Mc9zzgxYVts6s8e63YnUQACCCAQUAG+lQPacVQbAQTCK3DKihvbo2P2shWT/n0gQXpGZZ4KN6B3cPtxlRExO/tnDTZqdX5Kg9uL1DyiAjQ7AgJ8M0egk2kiAggET2D7fiNs2voH2fgVVvdd5cc2rWqZJZWpVv8A/k3UeO6Xc4+d4/aLnyZt4IDKOFAKAggggEBlBXJf1ZUtMBKl0UgEEECgBgIDko1291p72J9G7mh945UZOalEtcfOX7USxThl9AnYTforDjMnuNl3fMJ0L47TCP5BAAEEEPCdQNx3NaJCCCAQWAEqXh2BIwavY6+td6Bt0XtodTZQYqmj5+aO9Etcp6Psvf0Tu3VURWd+LGa2605x+9UZDbbqyrk3zlz+QQABBBDwqwBBjl97hnohgAACHoHVG/vZU+vuZ79caQtLWH0PsofM6e+pWfde9gxAkDN0iNlZpyZt4l4JSybKai8rIYAAAgjUWIAgp8bgbA4BBBAoVyCRC24mj9jSXl3vANu+70rlFtP99T5u6H4ZX5XQ86upHyd9+5gdsE/CfjOpwUauGvNjFakTAgEXoPoIVE+AIKd6tpSMAAIIVEVg3V4D7d9r720PjB5vG/YaVJVtdFToJi0rWqq5cgf8PSpXVEdVLnl+Yw+zCbvG7bxfNtjO4+KM3hj/IYAAAsETCHSQEzxuaowAAghUTmCX/qvaS+vtb9evPs6G1+hx09+cP7JyDciVVLkxoVxh3fxfDxLYafulwc1eeySsZ2M3C2R1BBBAAIG6CRDk1I2eDSNQNQEKjpBA3GJ25JB17Z0ND7WzV9rCelf5KWzrfjG0orrxdEWLK6swPVRgmy3i9tufN9jB+yWsX9+yimElBBBAAAEfCcR9VBeqggACCCBQpkCfXHDzqxFb2scbfc/OGbGVrdzQp8ySOl9txdmV/cMwiZoGOcu37Zvbxu1XpyftqMMSNqS2V/4tXxnmIIAAAghUTIAgp2KUFIQAAgjUX2BQotHOWmlz+yA3snPrqJ1t2z4rVrRSyZkVvoarraLVK6qwwQPN9puQsMvPbbDDD0zYiJV8eGNQUS0hEwIVEqAYBEIoQJATwk6lSQgggEAyFrcDB61lT6yzr7243kQ7dNAY65Gb1x2ZNZcMtCXzKxsQZGsY5Kw7JmY/PjLpPFBg92/HrU/v7miwLgIIIICAnwUqEeT4uX3UDQEEEIi8wCa9hthNo75lszY+wnlIwa79V7WklR6s7DR/VMUtMy3ZipfpLbCxh9mO34jbOWc12Kk/TtqmG8VM9+B48/AaAQQQQCB8AgQ54etTWuQbASqCgL8EBiYanYcU3D96vM3Z5Pt2w+o72m65gKehyBGeDb4cXvEGpVoqXqQNHGC2w9i4/eTopF06pcEOmZiw4cMqvx1KRAABBBDwrwBBjn/7hpohgAACVRNQwPODIevY1FzA88UmP7C719rDThq2kfN3d2IdbHXEnFz00MGykmZ7Mqeac2+6OZijkZm11ojZvuMTNvn0pF34qwY7bP+EbbxBzHo05MrnfwQQQACByAnEI9diGowAAgggsIxA33iDjV9hdbt01bH2ynoH2OyNj7D/W2NnO3rIurZGj/7teXt82qv9dSVfrFBisSvkqrTphkuDGl2CdtUFDXbGiUnbY+e4rTKioxCtkjWmLASqI0CpCCBQOQGCnMpZUhICCCAQCoGhyV52wMC17A+rj7N3NzzE9KQ2Xdq2ygYJG7pOzHpWeECnX2PHbHoS2jqjY7brTnE77sikXTC5wS7+TYP9+KilQY0eJqD7bjougSUIIIAAAgEXKKv6BDllsbESAgggEB2B1Xv0M13attmhCdvhp0kbf1GD7XN1g33r50nb+uiErTchbqtuGbchY2I2cPWY9c+NpvQdZtZroFmPvmbJnoWtGnMjMiusErO114rb9tvFncvN9PdqNCpzfi6Yuf7yBtP0p8cnbeJeCdtso5gNypVZuDTmIoAAAggg8LVA/OuXvEIgpAI0CwEEKi6Q6GE2cGTMVt06buvlApCtj0nYuNOS9q1fJG2XXydttykN9p0LG2zPyxps76sa7LvXL58mXNJgO09O2kHfS9j3Dkg4l5tts0Xc1lojZhrBqXilKRABBBBAIDICBDmR6WoaigACCCwrwDsEEEAAAQTCKkCQE9aepV0IIIAAAgggUI4A6yCAQAgECHJC0Ik0AQEEEPAKLG5usR+cfL6T9Nq7rNDrSedeZzfcet8yiz6bO892O/g0W3/cEc5U790Myqt13PdMEUAAAQSiIBCsNhLkBKu/qC0CCCDQqYCCmuMnXW5Pv/hmp/m0UMGKgpi7Hnhcb5dJ/3zwCZs4YZy9/tiNzlTvlUHBzuPPvGon/3Ci3pIQQAABBBDwpQBBji+7JZyVolUIIFB9gXMuu8nGbrWhnXLM/l1u7MiD9nCCmL12Hbtc3nc/mGFrrL6SM19TvdcbBTsqf+jgAXpLQgABBBBAwJcCBDm+7BYqhQACERKoWFPdS8gUvHS30DVHrmzvffipU4ymeu+O4uy5y3bOfP5BAAEEEEDArwIEOX7tGeqFAAIIlCCgS8+UfcqZR2vS7aRA5va7HzNdzqap3rujOHM+n2fbjj/OWaZ7f3SJXLc3SAEILCfADAQQQKB8AYKc8u1YEwEEEPCNgC4n0701CkqULrn2Nue+HN2fU04QosvR7r/lAudyNk3VUN2Lo2Dn5r8/ZFMmHe0sGz5skN165yNaTEIAAQQQqIUA2yhKgCCnKCYyIYAAAv4W0AiOHhLgJt2Ts/Wm69pVU0603r16drvy7ihOn949bdacL9rL02Vs7W94gQACCCCAgE8ECHJ80hE1rAabQgCBCArofho9Etq9rK0UAq3rjuIoYNLojbu+RpDc10wRQAABBBDwiwBBjl96gnoggECdBaK3eQU8urRNl7np8jbdZ/P6Wx8sB3HpH263wybuYrqETQv1+Ojzr7rVuSdHozoH7b2TZpMQQAABBBDwjQBBjm+6googgAAClRPQE9b+eOnp7ZeqKUDRvTWa725Fr19/7Ebn3hpNn7znalt/7ZHu4vapLoXbcbtN29+7ZWkd7zbaM/AiXAK0BgEEEAigAEFOADuNKiOAAAIIIIAAAgjUV4Ct+1uAIMff/UPtEEAAAQQQQAABBBBAoEQBgpwSwSqXnZIQQAABBBBAAAEEEECgGgIEOdVQpUwEEChfgDURQAABBBBAAIFuChDkdBOQ1RFAAAEEEKiFANtAAAEEEChegCCneCtyIoAAAggggAACCPhLgNogUFCAIKcgCzMRQAABBBBAAAEEEEAgqAIEOUHtOeqNAAIIIIAAAggggAACBQUIcgqyMBMBBBBAAAEEEEAAAQSCKkCQE9Seo94IIIAAAvUQYJsIIIAAAgEQIMgJQCdRRQQQQAABBBBAwN8C1A4BfwkQ5PirP6gNAggggAACCCCAAAIIdFPAN0FON9vB6ggggAACCCCAAAIIIICAI0CQ4zDwDwK+FaBiCCCAAAIIIIAAAiUKEOSUCEZ2BBBAAAE/CFAHBBBAAAEEOhYgyOnYhiUIIIAAAggggECwBKgtAgg4AgQ5DgP/IIAAAggggAACCCCAQFgE8oOcsLSLdiCAAAIIIIAAAggggEBEBQhyItrxNLtUAfIjgAACCCCAAAIIBEWAICcoPUU9EUAAAT8KUCcEEEAAAQR8KECQ48NOoUoIIIAAAgggEGwBao8AAvUVIMiprz9br7DA6299YNuOatlA2wAAEABJREFUP87WH3fEMunRJ16s8JbMPps7z3Y7+DSrRtkVr2yBAlVv1V/tcBcvbm6xH5x8vk069zp3Vl2mTz2XsbumpotOsz+rSzXZKAK+E2j5x03WctsNRSffNYAKIYBA2AVq1j6CnJpRs6FaClw15UR7/bEbnaTXx0+6PLDBSC3devfqaX+89HSbcubRtdzsctt6+vmM3X1/8WnOZ9nlyqjUjELBYKXKLqYcBZ4n/OIKUwBfTP5q5lFdFATfcOt9JW+mO+uWvLECK6gfdfJD0wKLazZLdjoRU63+dIKcv/3JWopM1ry4Zm0vZkPqn/yTL8WsV6s87n6sz4Fe12q7ldqOTmD52bdS7aQcBCRAkCMFUnAEyqjp1puua0rvffips7Z+mPQDpYMNZ8ZX/+jLX0lvvXk0TwdHSnrtLj/9nGvt45lzTAGUlrk/HCpX5asM5VXy/nBrvpYrn8rzrqu8mq95btK6mt9RUn6Vd8/DT7aPXmkdzXfL0NR7YKUDrElTrnPqP26/k5z1VBdv3dztufNUhpLbTnd5GKb5bVQ75eGHtr3/0SynGqNWG+5M+ad0AY1WXn3jXXbbtZNtx+02Lb2A3Br6PBXaJzRP+4uSPne5rM7/+oxN/OFkZ8TXmfHVP0cetIcddch4u/nvD301J5gT9zOT/30ga81T0ms/tU59pdSdOt165yM2fNgg52SQTgp5y1LZ2k+884p9rfW0fn5+7VPat5TyTbWP6Xtdy/QboD5x11dZKtN97051AmuzDUfbPx98wp3FFIHQChDkhLZraZgr0LS4xWbOnuu+LWl6ybW32c47bOGMCD12x2X2wqvTTT8c+nE7/+fH2KojhplGijRqdP8tF9jQwQOKLt9btruuyn78mVft2anXONvUQZmCEf3QdVbw0y++aU88+5qzjuqy41cHcm7dNE8HVqf+6mrnoGv9tUfalElHO/VXu7RcP37529CPpoI4/agrj9LECePssJ9MccrJzx/E9zoQ2/fIs50DF7XPTX379PLF6MlTL7xhG6+3lmmf68y3FstUB4306UC9Ftur1DZ0QDd61Mqm/b6zMvX50wGjDhDdfPrsaZ4+r+48d6qDzEWLm53Pqz6r/5j6X9NnRssVxBx3xF4FvxMO2nsn03paX3mDnPQ5ka/bBr3WPPd9mKb6rrj/0WfskH13XqZZ7n5z1wOPLzO/mDdd7V/nXXmLE5zre0nfvTq55t3H9D2uZfqO1u+Atqn9SvuX9jO9z08n/3Ci6XdG7clfxnsEwiRAkBOm3qQtBQX0o6sFe+6ynSYlpVOO2b/9zK8CGP3I6MfB/ZEpqbC8zN6ytUg/OPoBPfXYA9oPaHVQtuPYTe2hfz+nLB0mjVT9/KTDllmuA1E32NECtb9f39425/N5eltU0ijCjFmfL/OjrnK08mtvva9J4JP2jxErDrZ8v0knHFrwoFh9r7OmOvBV0gGOi6DXmucmHcC4y3TgrMvO3HV1BlYHI+7yQlNt6+U33rFtNluv0GJnnspQWe42NXW3q20qORlz/2i+tq9yc2+de6+UX0llqCzNVyrUFq2n9bVMeTpKbj6Vq+Stg3cd1UfLlfLPUmsdzVfy1k2fE+W9/Po72u+/U520TW/Z7mvN12dWJyvceflTd1uarwNGb8Cvz5Dm6fOq5d40Z+6X1rd3L+fzOmzIAFu0qNl0UkWOOsjU59Kb332tYFGBqwJYd15QpxPH79B+wKy+kbXm5bfHNVZ/6rV3ufYnzVfSSRXvMvWf+lfLlJRXy7Ut737g7iNarnxuevSJF5XduVxZQYiSlqlMla2k15qnpPWdFQr8o++8frmTH/mjqvqu1T6y165jC6zV+azO9i/tH5tvNKb9e0jfA/o+1vey6q19bNjggc4G1hy5srlXKyjA3mf3bzr7pbMw7x/9likoUnvyFvEWgVAJEOSEqjtpjCugH0r9YCm9+8EMc0dK3OXlTtdYfSVnVEgHMuWW0dF6Cj4+mTnH9j9msqnebtKPckfrdDVfBxNuObos7Y23PzAdmHW1nrtceVcePsS8P+r6gdTlDu4Pqps3iFMdKOigbOxWG3Z4QOBtl/Jr39IBgg5qNAqmwFQHtVq2qGnpWX0t00Hx+VfdusyI17MvTjMFsVqukbWLr/lr+5l/73bc1zqY0Wuvv967SQd6Gp1TWSpT9dHooru8s6nWdduhdRVMu/VRe9QuladlGqWY89mXnRXXvkwOXiOtv/F6ay7XTgU4GqVU2cqjEwiX/uF2p5zO6uZkyP3z8uvv2L9uv8QZRcm9NV1GpGl+cj+rG6w9aplF2oYOkvX5UACkOuhgdZlMXbzRAaYONNVmfX779u1lfXr3dC5F6+wgU8Xqu0T7ntbV+6CmYUMHmj4/Olmgg2a91jxvexQ4zJrzhdNXGqXWa303KY/2g+tvvqd9tEKjz5qvJBvvvqT9Uful9k8tV3L3gyfvudr5nuroM6hgQkGIkvpaI5Jav6vylcdN+s5T+xSkuvOqOdVvl7d8BdKxWMz5DlcdFGDrO1p5lFf7lGy0T3YUYCuvkoKirk6eKR8pzALhbxtBTvj7OJIt1A+lfsh0oKkgQT+kQYDo369P+4+96u8m75nlYtrhHsDpYEIHFSpHBwjFHgAXs40w5dHBQTHtUdChM6nu5SoK+HQZlM646qDjxKP2aw+WdNY1m82aDn7dshVIaHRO74vZps7I6oy/ytY6+UkHlZrnjq7pdbFJdT/hyP3as+tAX5d1KijQgZMCbrfuqvMBe+3UnrezF/lGyqt189ugAyyvh7ymvz/DCQo7q5vKUzps4i6OtcpVsKZ5hZLa4AYfhZb3z41uKlgptKyreXLRgeaWux/rnJxQYKP26yBTQZU7QuCOMnjL0zZVL++8oL7W/qfg46bbHzS99rZDgYqCOTc4UH+p73Tpr76n8vcD77qy7Ojz5uZTWSpT7zXt6jOofG4qpnw3r6ZuIKHXtUoKRjralr6HdKJAgbq+6xXY6DtD+6EuXdN8JTeg9JZTzPePNz+vEQiiAEFOEHvNB3UOShV0ZlZn7vRDoDNc3a23zuTp0iadrVXS6+6W6a7vPUvnzit3qgO7frmDN903pB/+QuXoIEt5Ci1z5ymPDjJ0MODO04GJDlDC9COpfnXb19V0wcIm54BWBw9KCqLddbSP6YBW8zUip7zuslKnctZBvw7+O1tX+6D2xc7ydLRMwb/qqqQz2m4+nfXW6JDaoGUa8VB93OXuVAdPWq7kzaNgXfuzm6+jqey0rpK2tXDR10/66qhuHZVV6nwFUhrh1WiQRrBUB404lFqOTkDoJIKS3NyDTAWgCr40X/dNaBs64C+1/CDkl6WCfQUyel2ozt7vC32veL97OjuQ12dI+4b6R0n7TKHy3XmlfgZLLd/dTjFT9bcb6KrubtJ3hOpZTBkKrDrKpyBbI1jaxzQype9pN8DW/WEaJdVJLgVA+jx1VA7zEQirAEFOWHuWdrUL6F6LdUevbrqsRwdqOujXwYfOLupHSBn1A1DMj6cuq/CeOdS6OhOpqZv0Y/7m9A9NPziap23qsiW97izp4ECXgeUHZKpbOQdfOhOvAy13m7oUSE+Dc99rmp9H87xJl0npcjUduLnzdVmKXutMtaZBTtoXdGDm3Re6as8qI4aZRsV0YOEmBdM6aDnq1AucBzpovg4wdLDfVXkdLVffFbr+Pz+/O/qSP7+r99qvtK+pnqqvRj+96zhtemzpY9i1X2r/8S7Xa+8BvgIG7cOarwNHBdp63VnSSKu27Sa3jK7q1lmZ+csUbLn3yuQv03vtAzpAVB30XgeiCt70utSkfUAHmTqj7g2cdVDft2+vUosLVH7tC9pnOqq010Mjhd6AtrMD+Y4+b4W2I/+jTr2gpM9gKeUrGPO2o1AdvPO8+5b2LzcpMFGA4s1b6LW2552vz5RGh7U/eee7r/U9rVEcvf909heaOKOd+r1z3vAPAhETIMiJWIdHsbn6odEBnM5473n4JOeJWQp8ZKHLTHRQo0BFIz6a5016opKWK+lsos7I6myt8qhc3V/x6OMvOvfQ7Hbwac6lNlruPQuuJ5F9/4DdtEqXSQcK3nW1XR2IdnU2P79g/YCqrjo7rzKUdF+E93I15dG23DyFDuzURtnpTKDKULr97sfspisnFXxqVH49gvBel9fMnD3XzrnspmWqO+WKvzj7inemgj4FHt4Dfh1Y/fWuR5xsCmrcAxAdyOlg31lQxj/aJxWAqQ86Wt0NNBUQKY+m3kBWB0nqOwXzSrqcSPnctEL/vqYgQO+9B28KMpQ0v9QkI29grO3qIQGaesvS5XE6aSA/d/4VN9zhfIb0vqO6aVkpyR3lkk1X6+kgXQei+hx2lbfQcvcgU32mkx1unjlzv3QeSuC+11SXOHZ2KaLyhCHJQvvx7Xc/5vSt9gPthwqcFRRrP9DIsE4Gqb3a7zVV0r7U2edNefJTZ59B7+dB65Vavvq0lBMi2kZ3kr73n3/l7fbvIe0z+myp3vnl6nPkBtja51dacZCTRd76DnDeeP6Rs+w9s3iJQOgECHJC16XRbpAO3HWWTIGGV0I/tDpbq2XK477XAY2SDmrc5F0v/0xzfrkqS2WqDPcstNZ3D5bc+bonwV3ublt5lDc/ab7Wc5PK13by87nvlV9tU7nuPE1VV7cMTd06aL6WK2ldLVNS+1WGytJ8LVdy5ymPktsOLQtD0oHW32/4telAQEGcm3QDc767LPKDPp053mCdNZwnIOlJSAqGVYYe/aoDrnKMdMCn+igA62x91f304w9q/1tNOnj0BrLuI2QVzH9r4im28fprtRen0QYF/nogheo79ZGn25cpUFNwrflKqot7YqA9UwcvZKTLJHXgqnW17dmffeGcUfauov1QgbjrpbzajtrUWd28ZRTzWvXRQbYO6rz5dVCoy4a03fzkDfgV7Gm5TnhotFevNc9bll6789QuvXenyi9LnRBRXbRMB55dPTVP+cKS9H2ioEb7mvYHjSzo+0btk5O7TFaa5yZ5dfR5c/N4p/q8dvYZ1OdJJzRUB11GpnVLKV8nFRY2NbeP0mt9JY20q+7aP7SfaL/S/qVlXSXtN1pX62l9vdY8raf2nPGTg9svj9UJJn225KLlbtL+dO1f/mnHHLqn8znTcr3Wd5PaKm85u/nd7xe1x53HtI4CbLpqAgQ5VaOlYAQQCIqADgoU3CmIc5P3IMwb2OXn9QahWsdd/4FbLzQlHajIQcuU9FpJBx3apsrTe2/SqIMOTHTA751f6LXKcbepAyBvHpWtbWi56qmbsvVe85X0WsuU/nbdr9qfQqg6K7/mKymf8ivptQ5avdvJf616y0zrKqndhdb11t3Np7LcvJqn5K2bW7bWVV4lld9ZnXRwq/ubvAee+W3Udtyk8lSukrbjznenmqdl3qR53vW0TO+1jiy1Pc1T0pPg9MAC7zzND1Jy+0jtzq+35qn/1VfuMnApTvcAAA4XSURBVNdCHnrtztdU7zVfSa+967rb0TIl11JlK5+2pTLcpPWVT0mfPyXX2V1Hy7Qfq2wlvdY8Jbd8tzzvVOvvtuNWztPzvPO172ldN3VWhnc9vVb93fXcqeZpmZJeu/PVXtVB871JbbjiNyc4J1rc+Wqz6qF1ZeLO11Qj0Qr8C5Wl5SQEwiJAkBPsnqT2CIRSYOvN4zZht+LTsKGxUDlo1IFLSSrXpTqY0x/m1KiRe5a8cqWXVpLO+usyPT0Zq7Q1i8vdc5/DrOd3v190sl69iyuYXI6ARkc14viDk89f7rHoTgaf/6NRSo2yKvD3eVWpHgLdFiDI6TYhBYRRQGfGdHZPZ+jC2D6/t2mbLeK21+6JotOKQ/NbFOz3OvOqM7iltkIH8zrbW866pW4raPllorPamtaz7vpO0Rl2nWmvRj16KsjZ/0jrWWSqRh3CXKb726DfB70OWlv13aLvCH1XBK3u1BeBUgUIckoVIz8CCCCAAAJBFaDeCCCAQEQECHIi0tE0EwEEEEAAAQQQQKCwAHPDJ0CQE74+pUUIIIAAAggggAACCERagCCnIt1PIQgggAACCCCAAAIIIOAXAYIcv/QE9UAgjAK0CQEEEEAAAQQQqIMAQU4d0NkkAgh0LnDzF2/b5JnPFp2mt87vvECWIuAzgWpVZ9p9GXvjrnTRqVr1oFwEEECg3gIEOfXuAbaPAALLCdzyxXT71afPFZ2mt8xbroxKzVjc3GI/OPl809836arMz+bOs90OPs1K/Vss+iOV244/rqhtdFUHvy2Xm9qmNvqtbmGsz7R7cwHO3blAp8iUagmjQtdt0v6o/VL7Z9e565dD9VM9Vd/61aKmW2ZjCFRMgCCnYpQUhAACQRTQQYSCGAUzbv11QLHrQT8zTd151Zpquxdf81ebMulo099Q8W5H25/4w8mm4Mk7X+uozuuPO8KU1Abv8s5eq0wdNGk9JZWj8tx1tC0Falqm5A3YvMuUR+/d9VQH/aFB9707VZuOOmS83fz3h9xZTEMioP7XfqD9xE35+5Mfm6r9vaPPXFf17erzo7Jl4Hroc+GW6fWSm967y5SPz4+rwRSBygiEK8ipjAmlIIAAAu0C+oN/+sN/Olhvn1nBF0+/+KZT2tabrutM9Y8OfnQQtP8xk23hosWatUw657KbbPiwQaY/bvnYHZfZ7Xc/VvTo0VMvvOEEVFpXSeWoPG1AB2inn3OtTZwwzin7tmsn23lX3tIe7P3zwSfalymP3ms91ffxZ161k384UW+XSwftvZMtWtzcXs5yGZgROIH8fUX7ktJao1a29z+a5ev2FPrMFVvhzj4/KkOfJX2mZJH/2dTnRZ8bLdNU77UOnx8pkBCovABBTuVNKREB3wlQofIFdDCnM7M60+qWotfumVpNFZDoQMVd/sob75k7WqJ1VYa7LH/60L+fs7FbbWgKptxlQwcPsPtvucAUZPTr29ud7Uy1nenvz7BD9t3Zea+8m2042lSOM6OLfxSs7bjdpu251hy5ss2a84Wpjjo4XdjUbHvusp2zfNRqw23l4UNMB3aa8e4HM2yN1VfSS2eq93qjgzW1QXXR+/yktm283lrt5eQv533wBLSvzJj1uW2z2XrLVH7SCYfa+muPbJ/n/ay4nwXtw/rMeEcJNULiHT3VMn22lJRX66hQladyTvjFFaZl7uiHliuf5ilpfeUvlPRZ0f6q/dJd7q6vsvVZcOfnTzv7/KiMzj6b+rzw+ckX5T0C1RMgyKmeLSUjgEAIBXTwdP3N9zgBiM7IXjXlxOVa+fLr79i/br/Enp16jbPs1jsfcab5/+hgSiMc+QeK+fm87+d8Ps8WLGzyzjJvoLLMgi7eaPsagdGZZx3wzZn75TIjR5qnZTo4U1HaznsffqqXpqne68BOZbiBkbOwwD86uFM+bbPA4nrNYrtlCgwbMsBisZjpsq+O+lQBiUYZNaKhz4r2JX0WFAxrJEPBhrt5BdKbbzTGCZD0GZs05br2z5jyXvqH292sppGYfXb/pjPaOOXMo53LOQ/7yZT2UUadHLj6xruc+e0rffVCdS30mVOddGLh1GMPsG9NPMW5D095v1qt4ETLtU+rXfqsdPXZ1OdFnxsVpqne8/mRBgmB6ggQ5FTHlVIRQCBAAjpo2nL3Y50zwzoLrMvE8gMJtzk6MNtx7KbOwZg7L3962MRdnJEZHfjoACh/ufu+aXGLLVrUbDpgdOcVM+3fr0/J61jefzoDrjZr9s9POkwTJ41YcbD16d3TeZ3/jwIZHbTKSFO9d0dxdIDX2ejVsMEDrW/fXvlF8j6gAgoKbrpyks2cPde0H2mfUFJgoya5AYACFOXVvJ132MIUFGiZAnuNeuggX+81X8uVL/8z5s2r5bq0U0mvlV57631NlhmB7Nenl7nznYVf/dPVZ06jUE/ec7XpM6x26XPy1arLTDRfyzXT+/np7LOpz4s+N3LSVO/9+/lRy0gIBFuAICfY/UftEUCgAgI6YNKoi842K+lMsA5WOipaZ2A7WlaL+QrAFFQU2pYOvnQQlZ90dtybX2fA1VZdtnP8pMudy9W0XAetOhDU6/ykg1Wd7dZ6mmq5Dk51sKYHC+jhCVqmwE5n7LWcFF4B7/6gfteo5iXX3rbM/WF67+6L2s9cDV0K6QYiuvRNl0lusPYod7Hd9cDjy5x0KHRvWnvm3IuPZ86xcfud5Kyj4EMnLnKzy/5fl3SqTQq8Cl3C1tHnp7PPpteLz0/ZXcOKCBQtUJUgp+itkxEBBBAIoIB7+VZ3q64RE41udBSwFCpfoz75AZjqo8BCI0fuwZcO0LxJB22FytNZch1gKrDRaEs/zz1AOsOu+3U6Curcs9Bqh/K55XeU313ONJwCOlmgpEux3BYq8PHuh3qIh/ZTJQXYGrXRpWq77biVKQhw1zvlmP2dy9HcdRUUeJe7+dyptus9UaH1Cu3z2lf75kYUu/rM6aSAgjPVz62zuy3vdJnPz5AB1tln07ueXvP5kQIJgeoJEORUz5aSEcgX4H0IBHRm94VXp7df76+DoHKbpQO9vr17lXRDvg70Ro9auf2RzLrcR/VRvYqpx5Qr/rLMU840AqMz6jr4c8+u6+BLZekMe6Gby7VM23VHcdQOBVmar6SgS1Nv0oGsHj6gvN75vA6mgB4UoJv/FQi7LfDuL+pnBTHnX3Vr+2dFeS+//o72UUONAOqStamPPL3MAwy0L+u+N23DLfuKG+5oL8ed5041AqQRSO/ooYIUJTePO1W9OvvMaZu67PKm2x907qnTSQN3XU07+/yU8tnk8yNNEgLVFSDIqa4vpSOAQMgEdHZYTzNzL41R8zq7j0XLO0s6oFOwoANAN58OgPSkKN0b9MbbHziX4egyNHe57gHQyInONKseuu9B9XKXdzbddov1TeVqXSWVo7PtOvhTOv/nx5juF9Ay5TvjJwcXvP9IN4LrvgUd2FnuPz0+Wge0Wk9lHrT3Trm5S/9X215+451lDmSXLnH/ZRo0AQXEup9Ml4apz5WOOvUCu2Tyj9v3Fz2JTPum9lEtV96+fXo596upvdp3FLD3z40eqjzNU9K+rEsftf9pPSXtU8qv5flJ83V/kLvfKr/2RQU/+Xn1vrPPnB6koIeGdDR609nnR2UX+9nk8yMtEgLVFSDIqa4vpSOAgM8FdCCWf0Cjm48fuPVC52BNB/5arnxuU3R2V5fDKOnSLI1iKJ8OtnRZjQ7SvHm967rz3akus9Fr7z0Ebjkq303apvIpaVuqk7uss/KV35tUN3c9TVWOynPz5G9b+d1l3qnq413mXS+/TJ1h19lzuXrL4HVwBbTPqJ+1D7lJN+zn97H2TXe5pnrvbbX2I5Wj8rzztW8pv5uUT8u1fqH83v1P6+hzqHlaJz919pkrVLZ3/fx65edXOzRPdVBSfb3ru6/VHpXlvlddVWeto/VVjrMs9w+fnxwC/yNQhgBBThlorIIAAtUVOHjQaPvlSlsUnUb3HFDdCnlKd0dZdLZYSaMwOnvryVLSSx3M6LG1emSu+2SqkgrweWa1SZceHfLV3/XxeXUDX711vpOw9SbEi07JnoFvcskNCNJnjs9Pyd3LChERKKaZBDnFKJEHAQRqKnDIoDE2ecSWRafRjSvUrH7eM66FzrqWUxGd/dZZ8I7O+pZTpl/WUZvUNrXRL3UKcz3W2SMX4OyVC3SKTGG26Kxt2h+1X2r/7CxfvZepfqqn6lvvurB9BIImQJATtB6jvhUQoAgEEEAAAQQQQACBMAsQ5IS5d2kbAgggUIoAeRFAAAEEEAiJAEFOSDqSZiCAAAIIIIBAdQQoFQEEgidAkBO8PqPGCCCAAAIIIIAAAgjUW8DX2yfI8XX3UDkEEEAAAQQQQAABBBAoVYAgp1Qx8ldOgJIQQAABBBBAAAEEEKiCAEFOFVApEgEEEOiOAOsigAACCCCAQPcECHK658faCCCAAAIIIFAbAbaCAAIIFC1AkFM0FRkRQAABBBBAAAEEEPCbAPUpJECQU0iFeQgggAACCCCAAAIIIBBYAYKcwHZd5SpOSQgggAACCCCAAAIIhEmAICdMvUlbEECgkgKUhQACCCCAAAIBFSDICWjHUW0EEEAAAQTqI8BWEUAAAf8LEOT4v4+oIQIIIIAAAggggIDfBaifrwQIcnzVHVQGAQQQQAABBBBAAAEEuitAkNNdwcqtT0kIIIAAAggggAACCCBQAQGCnAogUgQCCFRTgLIRQAABBBBAAIHSBAhySvMiNwIIIIAAAv4QoBYIIIAAAh0KEOR0SMMCBBBAAAEEEEAAgaAJUF8EJPD/AQAA//+tpOrrAAAABklEQVQDAGIsz1Q2frUOAAAAAElFTkSuQmCC" }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Base\n", "rs = rupture_summary.copy()\n", "\n", "# Classes simplifiées\n", "bins = [0, 0.01, 0.10, 0.30, 1.01]\n", "labels = [\n", " \"Clean / quasi-clean (≤1%)\",\n", " \"Moderate (1–10%)\",\n", " \"High (10–30%)\",\n", " \"Severe (>30%)\"\n", "]\n", "\n", "rs[\"rupture_class\"] = pd.cut(\n", " rs[\"rupture_ratio\"],\n", " bins=bins,\n", " labels=labels,\n", " include_lowest=True\n", ")\n", "\n", "# Distribution en %\n", "dist = (\n", " rs[\"rupture_class\"]\n", " .value_counts(normalize=True)\n", " .sort_index()\n", " * 100\n", ").round(1)\n", "\n", "# Donut chart\n", "fig = go.Figure(\n", " data=[go.Pie(\n", " labels=dist.index,\n", " values=dist.values,\n", " hole=0.45,\n", " textinfo=\"percent\",\n", " hoverinfo=\"label+percent\"\n", " )]\n", ")\n", "\n", "fig.update_layout(\n", " legend=dict(\n", " orientation=\"h\", # horizontale\n", " yanchor=\"top\",\n", " y=-0.15, # en dessous du graphe\n", " xanchor=\"center\",\n", " x=0.5\n", " ),\n", " legend_title_text=\"Rupture ratio\"\n", ")\n", "\n", "fig.show()\n" ] }, { "cell_type": "markdown", "id": "e52cd650-df05-490d-af59-e66c058f955d", "metadata": {}, "source": [ "## AUM–FLOW CONSISTENCY & DISCONTINUITY DETECTION" ] }, { "cell_type": "code", "execution_count": 8, "id": "a7efe494-f5fa-43f8-8446-942fc2d3bd4c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Detection threshold epsilon (trimmed 99th percentile): 40.03%\n" ] } ], "source": [ "# ------------------------------------------------------------\n", "# 1. Keep relevant columns\n", "# ------------------------------------------------------------\n", "stocks_clean = stocks[\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\", \"Quantity - AUM\"]\n", "].copy()\n", "\n", "flows_clean = flows[\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\", \"Quantity - NetFlows\"]\n", "].copy()\n", "\n", "# ------------------------------------------------------------\n", "# 2. Date formatting\n", "# ------------------------------------------------------------\n", "stocks_clean[\"Centralisation Date\"] = pd.to_datetime(stocks_clean[\"Centralisation Date\"])\n", "flows_clean[\"Centralisation Date\"] = pd.to_datetime(flows_clean[\"Centralisation Date\"])\n", "\n", "# ------------------------------------------------------------\n", "# 3. Aggregate flows per day\n", "# ------------------------------------------------------------\n", "flows_clean = (\n", " flows_clean\n", " .groupby(\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"],\n", " as_index=False\n", " )[\"Quantity - NetFlows\"]\n", " .sum()\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# 4. Merge stocks and flows\n", "# ------------------------------------------------------------\n", "df = stocks_clean.merge(\n", " flows_clean,\n", " on=[\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"],\n", " how=\"left\"\n", ")\n", "\n", "df[\"Quantity - NetFlows\"] = df[\"Quantity - NetFlows\"].fillna(0)\n", "\n", "# ------------------------------------------------------------\n", "# 5. Sort and reconstruct expected stock\n", "# ------------------------------------------------------------\n", "df = df.sort_values(\n", " [\"Registrar Account - ID\", \"Product - Isin\", \"Centralisation Date\"]\n", ")\n", "\n", "df[\"prev_stock\"] = (\n", " df.groupby([\"Registrar Account - ID\", \"Product - Isin\"])\n", " [\"Quantity - AUM\"]\n", " .shift(1)\n", ")\n", "\n", "df[\"prev_flows\"] = (\n", " df.groupby([\"Registrar Account - ID\", \"Product - Isin\"])\n", " [\"Quantity - NetFlows\"]\n", " .shift(1)\n", " .fillna(0)\n", ")\n", "\n", "df[\"expected_stock\"] = df[\"prev_stock\"] + df[\"prev_flows\"]\n", "\n", "# ------------------------------------------------------------\n", "# 6. Compute accounting gaps\n", "# ------------------------------------------------------------\n", "df[\"gap\"] = df[\"Quantity - AUM\"] - df[\"expected_stock\"]\n", "df[\"gap_abs\"] = df[\"gap\"].abs()\n", "\n", "# Relative gap normalised by previous stock\n", "df[\"gap_rel\"] = (\n", " df[\"gap_abs\"] /\n", " df[\"prev_stock\"].abs().replace(0, np.nan)\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# 7. Calibration sample (valid regime)\n", "# ------------------------------------------------------------\n", "valid_gaps = df.loc[\n", " df[\"gap_rel\"].notna() & (df[\"prev_stock\"] > 0),\n", " \"gap_rel\"\n", "]\n", "\n", "# ------------------------------------------------------------\n", "# 8. Robust, data-driven threshold (epsilon)\n", "# ------------------------------------------------------------\n", "# Step 1 — trim extreme breaks to avoid calibrating on resets\n", "gap_rel_trimmed = valid_gaps[\n", " valid_gaps <= valid_gaps.quantile(0.90)\n", "]\n", "\n", "# Step 2 — define epsilon on the upper tail of the trimmed distribution\n", "EPSILON = gap_rel_trimmed.quantile(0.99)\n", "\n", "# ------------------------------------------------------------\n", "# 9. Detect discontinuities (diagnostic rule)\n", "# ------------------------------------------------------------\n", "df[\"rupture_flag\"] = (\n", " df[\"prev_stock\"].notna()\n", " & (df[\"prev_stock\"] > 0)\n", " & (df[\"gap_rel\"] > EPSILON)\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# 10. Remove end-of-sample edge effects\n", "# ------------------------------------------------------------\n", "last_date = df[\"Centralisation Date\"].max()\n", "\n", "df.loc[\n", " (df[\"rupture_flag\"]) &\n", " (df[\"Centralisation Date\"] == last_date),\n", " \"rupture_flag\"\n", "] = False\n", "\n", "# ------------------------------------------------------------\n", "# 11. ISIN-level summary\n", "# ------------------------------------------------------------\n", "rupture_isin_summary = (\n", " df.groupby([\"Registrar Account - ID\", \"Product - Isin\"])\n", " .agg(\n", " n_ruptures=(\"rupture_flag\", \"sum\"),\n", " total_obs=(\"rupture_flag\", \"count\"),\n", " rupture_ratio=(\"rupture_flag\", \"mean\"),\n", " max_gap_abs=(\"gap_abs\", \"max\"),\n", " max_gap_rel=(\"gap_rel\", \"max\")\n", " )\n", " .reset_index()\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# 12. Account-level summary\n", "# ------------------------------------------------------------\n", "rupture_summary = (\n", " df.groupby(\"Registrar Account - ID\")\n", " .agg(\n", " n_ruptures=(\"rupture_flag\", \"sum\"),\n", " total_obs=(\"rupture_flag\", \"count\"),\n", " rupture_ratio=(\"rupture_flag\", \"mean\"),\n", " max_gap_abs=(\"gap_abs\", \"max\"),\n", " max_gap_rel=(\"gap_rel\", \"max\")\n", " )\n", " .reset_index()\n", ")\n", "\n", "# ------------------------------------------------------------\n", "# 13. Outputs\n", "# ------------------------------------------------------------\n", "df.to_csv(\"aum_flow_gaps.csv\", index=False)\n", "rupture_isin_summary.to_csv(\"rupture_isin_summary.csv\", index=False)\n", "rupture_summary.to_csv(\"rupture_summary.csv\", index=False)\n", "\n", "print(f\"Detection threshold epsilon (trimmed 99th percentile): {EPSILON:.2%}\")\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "d7454212-1493-4715-a436-c331931f92fa", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
| \n", " | Registrar Account - ID | \n", "Product - Isin | \n", "n_ruptures | \n", "total_obs | \n", "rupture_ratio | \n", "max_gap_abs | \n", "max_gap_rel | \n", "
|---|---|---|---|---|---|---|---|
| 59545 | \n", "200127410 | \n", "FR0010135103 | \n", "384 | \n", "436 | \n", "0.880734 | \n", "295985.42 | \n", "3371.158214 | \n", "