{ "cells": [ { "cell_type": "code", "execution_count": 51, "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": 52, "id": "cfd11919-0941-400e-a516-72871881f733", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_1311/1940519970.py:1: DtypeWarning:\n", "\n", "Columns (1,2,3,4) have mixed types. Specify dtype option on import or set low_memory=False.\n", "\n", "/tmp/ipykernel_1311/1940519970.py:2: DtypeWarning:\n", "\n", "Columns (1,2,3,4) have mixed types. Specify dtype option on import or set low_memory=False.\n", "\n" ] } ], "source": [ "stocks=pd.read_csv('stocks.csv')\n", "flows = pd.read_csv('flows.csv')" ] }, { "cell_type": "code", "execution_count": 53, "id": "b99e3402-fe26-4f4e-8c1c-5f07847bce94", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_1311/3613746644.py:1: DtypeWarning:\n", "\n", "Columns (1) have mixed types. Specify dtype option on import or set low_memory=False.\n", "\n" ] } ], "source": [ "merged = pd.read_csv('merged.csv')" ] }, { "cell_type": "code", "execution_count": 54, "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": 55, "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": 56, "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": 57, "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": "iVBORw0KGgoAAAANSUhEUgAAAxoAAAFoCAYAAADQNY2xAAAQAElEQVR4AezdB5xU1d3/8d+UBZalCCigYAF779gVe3lsUSGxRRProz722BOJGiyJJhpji0nMo8a/GhM1Ro2aaB5jNxgL9oZSBAFBlt1ld2fmP98Ld7077OzO7LRbPrw4O3fuPefcc97nzsz93TITz/APAQQQQAABBBBAAAEEECizQNz4hwACPhOgOQgggAACCCCAQPAFCDSCP4b0AAEEEECg0gLUjwACCCBQtACBRtFkFEAAAQQQQAABBBCotQDr978AgYb/x4gWIoAAAggggAACCCAQOAECjcANWakNpjwCCCCAAAIIIIAAApUXINCovDFrQAABBLoXYCkCCCCAAAIhFCDQCOGg0iUEEEAAAQQQKE2A0gggULoAgUbphtSAAAIIIIAAAggggAACOQJlDjRyaucpAggggAACCCCAAAIIRFKAQCOSw06nIyVAZxFAAAEEEEAAgRoIEGjUAJ1VIoAAAghEW4DeI4AAAlEQINCIwijTRwQQQAABBBBAAIHuBFhWAQECjQqgUiUCCCCAAAIIIIAAAlEXINCI+hZQav8pjwACCCCAAAIIIIBAFwIEGl2gMAsBBBAIsgBtRwABBBBAwA8CBBp+GAXagAACCCCAAAJhFqBvCERSgEAjksNOpxFAAAEEEEAAAQQQqKyAvwONyvad2hFAAAEEEEAAAQQQQKBCAgQaFYKlWgTCKkC/EEAAAQQQQACBQgQINApRIg8CCCCAAAL+FaBlCCCAgC8FCDR8OSw0CgEEEEAAAQQQQCC4ArRcAgQaUiAhgAACCCCAAAIIIIBAWQUINMrKSWWlClAeAQQQQAABBBBAIBwCBBrhGEd6gQACCFRKgHoRQAABBBDolQCBRq/YKIQAAggggAACCNRKgPUiEAwBAo1gjBOtRAABBBBAAAEEEEAgUAKRCjQCNTI0FgEEEEAAAQQQQACBAAsQaAR48Gg6AiEQoAsIIIAAAgggEFIBAo2QDizdQgABBBBAoHcClEIAAQTKI0CgUR5HakEAAQQQQAABBBBAoDICAa2VQCOgA0ezEUAAAQQQQAABBBDwswCBhp9Hh7aVKkB5BBBAAAEEEEAAgRoJEGjUCJ7VIoAAAtEUoNcIIIAAAlERINCIykjTTwQQQAABBBBAoCsB5iFQIQECjQrBUi0CCCCAAAIIIIAAAlEWINDo/ehTEgEEEEAAAQQQQAABBPIIEGjkgWE2AggEUYA2I4AAAggggIBfBAg0/DIStAMBBBBAAIEwCtAnBBCIrACBRmSHno4jgAACCCCAAAIIRFGgWn0m0KiWNOtBAAEEEEAAAQQQQCBCAgQaERpsulqqAOURQAABBBBAAAEEChUg0ChUinwIIIAAAv4ToEUIIIAAAr4VINDw7dDQMAQQQAABBBBAIHgCtBgBV4BAw5XgEQEEEEAAAQQQQAABBMomQKBRNspSK6I8AggggAACCCCAAALhESDQCM9Y0hMEECi3APUhgAACCCCAQK8FCDR6TUdBBBBAAAEEEKi2AOtDAIHgCBBoBGesaCkCCCCAAAIIIIAAAn4TyNseAo28NCxAAAEEEEAAAQQQQACB3goQaPRWjnIIlCpAeQQQQAABBBBAIMQCBBohHly6hgACCCBQnAC5EUAAAQTKJ0CgUT5LakIAAQQQQAABBBAorwC1BViAQCPAg0fTEUAAAQQQQAABBBDwqwCBhl9HptR2UR4BBBBAAAEEEEAAgRoKEGjUEJ9VI4BAtAToLQIIIIAAAlESINCI0mjTVwQQQAABBBDwCjCNAAIVFCDQqCAuVSOAAAIIIIAAAgggEFWB3gUaUdWi3wgggAACCCCAAAIIIFCQAIFGQUxkQsD/ArQQAQQQQAABBBDwkwCBhp9Gg7YggAACCIRJgL4ggAACkRYg0Ij08NN5BBBAAAEEEEAgSgL0tZoCBBrV1GZdCCCAAAIIIIAAAghERIBAIyIDXWo3KY8AAggggAACCCCAQDECBBrFaJEXAQQQ8I8ALUEAAQQQQMDXAgQavh4eGocAAggggAACwRGgpQgg4BUg0PBqMI0AAggggAACCCCAAAJlEfBFoFGWnlAJAggggAACCCCAAAII+EaAQMM3Q0FDEPCVAI1BAAEEEEAAAQRKEiDQKImPwggggAACCFRLgPUggAACwRIg0AjWeNFaBBBAAAEEEEAAAb8I0I5uBQg0uuVhIQIIIIAAAggggAACCPRGgECjN2qUKVWA8ggggAACCCCAAAIhFyDQCPkA0z0EEECgMAFyIYAAAgggUF4BAo3yelIbAggggAACCCBQHgFqQSDgAgQaAR9Amo8AAggggAACCCCAgB8Fwhho+NGZNiGAAAIIIIAAAgggECkBAo1IDTedRaBWAqwXAQQQQAABBKImQKARtRGnvwgggAACCEiAhAACCFRYgECjwsBUjwACCCCAAAIIIIBAIQJhy0OgEbYRpT8IIIAAAggggAACCPhAgEDDB4NAE0oVoDwCCCCAAAIIIICA3wQINPw2IrQHAQQQCIMAfUAAAQQQiLwAgUbkNwEAEEAAAQQQQCAKAvQRgWoLEGhUW5z1IYAAAggggAACCCAQAQECjR4HmQwIIIAAAggggAACCCBQrACBRrFi5EcAgdoL0AIEEEAAAQQQ8L0AgYbvh4gGIoAAAggg4H8BWogAAgjkChBo5IrwHAEEEEAAAQQQQACB4AvUvAcEGjUfAhqAAAIIIIAAAggggED4BAg0wjem9KhUAcojgAACCCCAAAIIlCxAoFEyIRUggAACCFRagPoRQAABBIInQKARvDGjxQgggAACCCCAQK0FWD8CPQoQaPRIRAYEEEAAAQQQQAABBBAoVoBAo1ixUvNTHgEEEEAAAQQQQACBCAgQaERgkOkiAgh0L8BSBBBAAAEEECi/AIFG+U2pEQEEEEAAAQRKE6A0AgiEQIBAIwSDSBcQQAABBBBAAAEEEKisQPG1E2gUb0YJBBBAAAEEEEAAAQQQ6EGAQKMHIBYjUKoA5RFAAAEEEEAAgSgKEGhEcdTpMwIIIJAVSM+Zae2vv2xtz//dWp940FoevNOa777Zmm77qS3+xY+s8Sdn26KLT7KvzzzCFp5wgC2YuKMdf0abnXVxm110RZtd/rN2u/ZX7XbTb9vtjj+k7N4HU/bXJ9L2witpe/+jjM2bn12JP//TKgQQQACBKggQaFQBmVUggAACNRdobrL2N16xlgfusMVXnWcLj9/fvj5tohNMLP7FpdZ0+8+s5Q+32pKH7rbWpx7KBh//cIKQ1AdTLT3zM8ss/KqjC4sazeZ8aTbt84y9837GpryesX+9lLYnn07bn/+ast/clbJrbmi383/c5gQm501qs6uvb7fb70zZQ4+l7M23M9ba2lEdEwgggICZgRBGAQKNMI4qfUIAgWgLpFKW+uhdW/K3P1nTjVc4ZyQWHLOXNV5xlrXce7u1TXneMl8vqJrR/GyM8sHHGXvx1bT95fG0XX9ru516XptdcW273f9Qyt6YmrGWJVVrDitCAAEEEKiSAIFGlaArtRrqRQABBCSgwGHJo/dZ46Wn2oLDd7FFFx5vzb+5zlr/73HnjITy+CllMmaffpaxv/0jbTfc1m6nZQOPy65ZGngoKPFTW2kLAggggEDvBAg0eudGKQQQQCCfQPXmtzRnA4m/WePkc2zhiQdZ8x03WPs7r1dv/WVe02czlgYeuszq/Elt9seHUzZ9ZjYiKfN6qA4BBBBAoDoCBBrVcWYtCCCAQHkEUinn0qfF10+yBcfvb003Xm7t/3nJLJ0qT/0+qWXeV2aP/z1tk65utx9ObrNH/pa2ufN80jiaEUABmowAArUQINCohTrrRAABBIoUSL33pjXffm32zMWBzs3cbc89ZdYajRsbZs02e/DRlF1wWZv95Np251utiuQjOwIIIIBADQS6DTRq0B5WiQACCCDgEWh75VlbdN6xtuiH/21LnvizZRYt9CyN3uQnn2Wcb7W64Mdt9o9n09baFj0DeowAAggERYBAIygjRTsRWCrA34gItE95wRZdcLwt/umFlvr0w4j0uvBuzp1v9oc/puy8S9vs4cfTtrip8LLkRAABBBCojgCBRnWcWQsCCCBQkIB+60I/ktd41Q8s9fG7BZWJcqbGxWYPP5ayH/yozQk8vqret/Z62JlEAAEEEOhKgECjKxXmIYAAAlUWaH/3DefyKP3WhX4kr8qrD/zqdAmVLqXSjwTqtzmWROP2lcCPGx1AoGICVOwLAQINXwwDjUAAgagKpN6fao2XnW6NPzrFdMN3VB3K1e902pzf5rjoijZ7eUr2Sbkqph4EEEAAgaIFCDSKJgt1ATqHAAJVEkh9+oHp7MWiS06y9remVGmt0VnNwq/Nbvt9yq6+od1mfsFvcURn5OkpAgj4SYBAw0+jQVsQQCD8AtlD7i0P3OHc6K37MXruMDlKEfjgo4zzWxz3PZgyLqcqRZKyCCCAQPECBBrFm1ECAQQQ6JVAev6XtuhHp1jLvbeH7gf2egVSpULZ2M6eeDptl0xusw8+5nKqKrGHezX0DgEEChIg0CiIiUwIIIBAaQKtL/3Tvj7rSEu9/1ZpFVG61wL6RqprbkjZg39NmYKPXldEQQQQQACBggSqGWgU1CAyIYAAAqESaFpsi2/4sTVde7FZMz/2UOuxzWTMHnkibVf9ot3mzss+qXWDWD8CCCAQYgECjRAPLl1DoGcBclRSoP29N21h9ixG27+erORqqLsXAh9Py9ilV7XbK69xKVUv+CiCAAIIFCRAoFEQE5kQQACBIgRSKWu551bnK2szX80toiBZqymwpNXs1jtSdvudKX/dKF5NBNaFAAIIVFCAQKOCuFSNAALRE0jNnmGLLjzeWv58p5mu04keQeB6/OKrabv8Z222YCGXUgVu8GgwAlUSYDW9EyDQ6J0bpRBAAIHlBNo/fMcWnfc9029kLLeQGb4W+GKO2U+ua7dZswk2fD1QNA4BBAIlQKARqOEKWmNpLwLREdD9GI0/Pp0bvgM85PpWqsk/b7dPpnHfRoCHkaYjgICPBAg0fDQYNAUBBIIp4AQZl51htqTZ/x2ghd0KNGeH8Jpfpuz1twg2uoViIQIIIFCAAIFGAUhkQQABBPIJtE99zRoVZLS15svC/IAJtLWZ3Xh7yp5/mWAjYEMX2ObScATCKkCgEdaRpV8IIFBxgbY3XrHGn5xtRpBRcetqr0D38f/27pQ9/neCjWrbsz4EEAiPQIADjfAMAj1BAIHgCbRNed4WX3muWXv28Hfwmk+LCxT448Mpe+a5VIG5yYYAAggg4BUg0PBqMI0AAqUJRKR068v/tMXXXGiWYgc0CkN+131pe+nfnNmIwljTRwQQKK8AgUZ5PakNAQRCLtD63FPWdO0PzdIEGSEf6k7d04/6BfUG8U4d4QkCCCBQRQECjSpisyoEEAi2QPu7b1jTDZeZZTi6HeyRLL71umfjpt+mbOq7HSvZbgAAEABJREFUjH3xepRAAIEcgcg8JdCIzFDTUQQQKEUgPXe2Lb76fIKMUhADXlZXyv3y1yn78BOCjYAPJc1HAIEqCRBoVAma1ZRBgCoQqJFApnWJNV55rmUWL6pRC1itXwTa281+fnPKPp+R8UuTaAcCCCDgWwECDd8ODQ1DAAG/CDRdP8nSn3/il+b4qh1RbMySJQo22q1xcSaK3afPCCCAQMECBBoFU5ERAQSiKNDy5zut7ZVno9h1+tyNwNfZk1s3/zZlGd280U0+FiFQAwFWiYBvBAg0fDMUNAQBBPwm0Pbai9Zyz61+axbt8YnAex9m7K9PclbDJ8NBMxBAwIcCBBruoPCIAAIIeARSs2fY4ut+6JnDJALLCzz0aMo+/Jibw5eXYQ4CCCBgRqDBVoAAAr4VqFXDMi1NtnjyuWZLmmvVBNYbEAFdOaWvveV+jYAMGM1EAIGqChBoVJWblSGAgN8FdM29zmSkZ33u96bSPp8IuPdrpNORuIzKJ+o0AwEEgiBAoBGEUaKNCCBQNYHWJ/5s7f95qWrrY0XhEND9Go//nUAjHKNJLxAImoB/20ug4d+xoWUIIFBlgfTXC6z57purvFZWFxaBhx9L2ZfzCDbCMp70AwEEShcg0CjdkBoCKkCzEcgVaPnfG81auC8j14XnhQm0p8zu+EP2T2HZyYUAAgiEXoBAI/RDTAcRQKAQgfb33rTW/3u8kKzkqZxA4GvWJVTPv8y3UAV+IOkAAgiURYBAoyyMVIIAAkEWyKRT1nTT5CB3gbb7SOCeP6X41XAfjQdNKVWA8gj0XoBAo/d2lEQAgZAItP71PuNbpkIymD7oRnOz2b1/5qyGD4aCJiCAQI0FCDQqNABUiwACwRBIL5hnzf/v18FoLK0MjMALr6TtvQ8INgIzYDQUAQQqIkCgURFWKkUAAR8KdNmk5t/+wqyttctlzESgFIG77k+ZfpellDooiwACCARZgEAjyKNH2xFAoCSBtremWNuLT5dUB4URyCcwa7bZq6/xdbf5fJbO5y8CCIRZgEAjzKNL3xBAIK9Apr3dmm+5Ku9yFiBQDoE//TVl6TTBRjksqQMBBKokUMbVEGiUEZOqEEAgOAKtj91v6Tkzg9PgWrS0Tx+LDVrB4sNXtvhqYy2x7sa2+qoxGzSwFo0J5jq/nGv2Cmc1gjl4tBoBBEoWINAomZAKEHAE+BMggUw6bS0P/yFALa5CU/v0teRGW1i/w75nA350va1w1z+cNPj2R2zQjffboJ/9rw28/Gb74blJu+6KOrv9+jr72WV1dsk5STvt+KQdNTFh++wet1VWjlWhscFaxUOP8SN+wRoxWosAAuUSINAolyT1IIBAYATaX/6nZRZ+FZj2VqqhscFDrO9/TbQBP/5VNqj4ezbAuMH6TTwuG3BsaZY9m9HTelcYbLbGajHbbOOYjd8hbocdmLDLLkjaVZfW2bcPTtjaY2MWi/VUSyWX+6PuOV+a8SN+/hgLWoEAAtUVINCorjdrQwABHwi0PHq/D1pRmybEBgy0Prsf4AQVg299yOqPOd2S629a1sasONRsz13jdv4ZSfv5T+rse0ckbOMNoh1xPPQoZzXKupFRWXAFaHmkBAg0IjXcdBYBBFLTP7XUu29EDkL3WtQffaopuOh/0vnZsxZbmMUr/xEwoMFsh23idsZJSZt0ftI22TAWOXt1eF72BNq/XuJ3NWRBQgCB6AhU/lMmOpaV7Cl1I4BAmQSWRO3ejP4N1m/i8Tb4V3+0vgccblbXp0ySxVczepWYnX5i0rmvY/11ohdwPPUMZzWK32oogQACQRaIB7nxtB0BBBAoRiCzaKG1PvNoMUW6yev/RX33m+gEGP0OO9asbz/fNFj3dZxzatLOPyNhY1ePTsAxfabZp5/xVbe+2RBpCAIIVFyAQKPixKwAAQT8IrDkyYeKasr1jXU24av6jvRSa6Kj/Ly02X8v6NuxTPmOWdDPPmhfuuP8QHOyY5mm3YIqd87XfTvyufPL+RgbsqINmHSj1R97usUa/PtdtGuPjdtFZyftiMMSlviGtpwUvqvr2RezG47vWhWiBtEVBBDwlUDcV62hMQgggECFBPSVtkseK/wmcAUEasrvBzfb/UOabfLAFrujuW65AOHchlZnufL8foUWWzuZMZX9T3vCbsmWVdK05qm+Z5Ykbfu6lJNPz8ud6sbtbIN+frclN9is3FVXrL7ddorbxdmAY+gKFVuFbyp+4ZW0tbVxVsM3A0JDEECgogIKNCq6AipHAAEE/CDQ/tIzRX2l7bDsu+MZA9qsf/ZR7R8az+4cZjI2P71shmbmSfPTMetvGavPLlfStOYp2FDQMb5ve3ZJmf/X97f+p/3QGs6dbLH+DWWuvPLVrTY6ZpMuSNpG6y89I1T5NdZmDa2txg/41YaetSKAQA0Eev7ErEGjWCUCCCBQboGWx/5YUpUfZs9QNFnMhsbTner52eI+ziVS3sumhmaDEuVtzuZU0rTm6WzGZsmUKYjJLirb/9jAwTbwytutz857l63OWlTUvz5mZ56ctAkHJawKX4hViy466/wXl085DvxBAIHwC8TD30V6iAACURdIz5tjvf1KW91zoSBCAcUp/ds6LnlSsHDzCks6Lps6MHuW4rrGPs5lU1qmgOLkhfWmpGmNgXs2w3vvh/e+D+UpNulH9wZcfrMlVlmt2KK+zb/3bnE7//SkDfLb7SVlEnv/o4zpR/zKVB3VIIAAAr4VINDw7dDQMAQQKJdA+2sv9roq3XOhey90r8UdTUnLFxiMzwYausxKl0hpZYfWt3cEIZp2z2ZoeVMmZrr3Q/d93NeSdIITlSk2OTd9X3FLqIIM12DNMTE767+T1qd238brNqUij8++yFfdVgSWSiMnQIf9LUCg4e/xoXUIIFAGgbbXXyq5Fp2lWC+Ztump4u8h8N6boXs8dCmVGqTLqfoXX52KmvWttwGX3mCJEaOWPg/h31VHxezU4xIW662Rj01efa3zJXg+bipNQwABBHotQKDRa7ogF6TtCERHIJPJWNt/ig80dMnUbYuTHVB6PqU9YaMTGWeezmwoOU+yf3TGQkHDqPjS5dlZHf+1TJdPKVgZGk87N4pr4dKzG5oqPjWcc0Uoz2TkSmy4XtyOPCx8H1VfzjOb99Xy20pu/3mOAAIIBFkgfO/eQR4N2o4AAmUXSL/3ptmSlqLrVcDweTrh3Oit38i4aFE/0z0a2/RZesmLAoabmuo6lj+1JGFnNizp+JYqd4UKUN5vj9u+fdudWWsnM7ZO9szIMQvrTXVO7Ne+9OZwZ2lhf/p950Sr22ybwjKHINf4HRO2287h+7h6YyqBRgg2T7qAAALdCITvnbubzrIIAQSiJ9D2+su96rTut7h84Dc3e+t3MtwgQxUqYNC9G5qvpBvDdcZCy7xJ+S4c2NopANE9Gyqj5K3TWy7fdHLLHazfId/Ntzi087/zrbhtsG64rqF6+10un/L7Bkv7EECgNAECjdL8KI0AAj4XaPtP728E91vXYgMHW8Ppl/qtWVVpTzwes5OOSVi/vlVZXVVWMvW9jKVSmaqsi5UggAACtRCoQKBRi26wTgQQQGB5gUxTo6U+enf5BQGdU3/M6Rar7x/Q1pfe7IaGmB2wT6L0inxSg3687+NpPmkMzUAAAQQqIECgUQFUqkTAdwIRbVApX2vrN7LE2HUD/4N85TDdfZeYDRtajpr8UcfUd7h8yh8jQSsQQKASAgQalVClTgQQ8IVAWy++bcoXDe+iEf1PuaiLudGblUzEbMJB33wbWJAF1Pa33+PSKTmQEEAgnAIEGuEcV3qFAAJZgbYpL2T/Bv9/330OtcRqawa/I2XqwVabxWzNNcJxY/i06RlLpwk2yrRpUA0CpQpQvswCBBplBqU6BBDwh0Bq2oeWWbTAH40ppRV9+li/iceVUkMoyx5xWDju1UilzGZ/GY6gKZQbGp1CAIGSBOIllaYwAhIgIeBDgdQnH/iwVcU3qc8u+1lswKDiC4a8xOqrxmybLcPxETZ9JvdphHxzpXsIRFYgHO/SkR0+Oo4AAvkE0nNm5FsUqPl9Dzy8V+2NQqGD9g3HWY0Zs7h0KgrbK31EIIoCBBpRHHX6jEAEBNJzZgW+l8lNx1lixKjA96NSHRi+kpnObFSq/mrV+/kMAo1qWdd4PawegcgJEGhEbsjpMALREEjNnhn4jvbbb2Lg+1DpDmy1WfA/xqbPJNCo9HZC/QggUBsB/79D18aFtSKAQMAFAn9Go77BEpttE/BRqHzzx20Z/Bup5803W7KEYKPyWwtrQACBagsQaFRbnPUhEAIBv3chk05Z5qu5fm9mt+2r22I7i8WCvxPdbSfLsHDYkFg4Lp8K/gm4MowmVSCAQNgECDTCNqL0BwEELDNreuAV6jblbEahgxiGy6fmziv5jEahXORDAAEEqiZAoFE1alaEAALVEkiH4P6M5BbbVosr8OsJw+VTjYsDPwx0AAEElhNgBoEG2wACCIROIDUn2NehJNZYy+KDhoRuXCrVIV0+NXhQpWqvTr2NizmjUR1p1oIAAtUUINCopjbrKkiATAiUKpAOeqCx5vqlEkSu/IiVgn0/y2LOaERum6XDCERBgEAjCqNMHxGImEA64L+h4cPfzvD9FjR8xWAHGosINHy/jdFABBAoXoBAo3gzSiCAgM8F0gE/oxEfsYrPhf3XvJVKCDQuuzBpN15TZ9tt/c1H4jmnJu326+uclLvsoP0Sdst1S5f99LI6W2/tb4Ic1aWyxQot5tKpYsnIbxAg4H+Bb95V/d9WWogAAggUJJBZ3FhQPr9mio8Y5dem+bZdw1fqXdMUGKwy8ptAQbUc852EjVk9Zr+5K2XHn9Fm87/K2CEHJDoCii03jdkLL6edZc3NGdtmy6UfpQpA6upi9tcnUqqmqMTN4EVxkRkBBAIisPTdMSCNLUczqQMBBCIgEA/2W1t85OgIDFJ5u9ibS6fcMw/PvpDuaIzOTmy0Qdw+mZaxF15ZOv/fr2esvp/ZumvHnWBDwcSCr5cWWZh9XHHY0kBFAci776ft3Q+Kv7G7sbH4MktbwF8EEEDAvwLB/jT2rystQwCBwgXKnzOeKH+dVaoxNmCQxfo3VGlt4VnNiOFLd/YL7ZGCjJEjYvaHP/Z89mHOl0uDgBUGmRNEtLVlTNNal77tSr+B4Z7NeOnfS4MTLSsmtbQWk5u8CCCAQDAECDSCMU60EgEEihEI8C9qx/oPKKan5F0m0K+v2YAC4zNdGqUg4zd3tjuBw7IqnAedjfhidsa0XGc3nJk5f3SGY7txcef+jfr6mCm4cM9m6DIq994OXZaVUzTv03hxcVLeesq3gJoQQACB0gUINEo3pAYEEPCbQIDPaGRS7X7TDEx7MktPPPTYXl3qNGSw2bmnLb3he6ft4qZA5cgJCeeGcPceC3f5cUclLJk0cy+XeujRlJ18dptzj8YPftTmXNUNTaEAABAASURBVFKly6ne/yhj660Tt7/8LW0/u7HdFIQoqOmxQdkMAY6Ns63nPwIIVEUggCsh0AjgoNFkBBDoXiAW5Hs0CDS6H9xuljY1d7PQs+jaX7U7QYJu9FbSPRotS8zuvj/l3JehsxoKILRMSYHDVwvM3vug68ui3LMZWoXu5dCj6tCN4pouJCWCe7VfId0jDwIIRFSAQCOiAx+hbtPVKAoE+TqUds5o9GaTXZINFAo9o1FM/frK2z3Hxy3fTd7eezO+WpCx5palteuyK53RWPqs579Bjo177h05EEAgqgIEGlEdefqNQJgFuHTK56Nb/uYVejajkDUruNBvZ+hei2MOT9iTz6Tt9/9v+ZvGFUxsu1XcXnx16TdN6SyGpvfdI+5clqUzGl2V66oNCT6Nu2JhHgIIBFyAt7aADyDNRwCBLgSCfHi4uckybXwFURej2u2spuYCb9DoohYFA6ed1+ZcNqXF+lpbPddlU7oXQ/dkaH5uUmBx4WVt5l2uaZVR2R9dWfjZqSBvsrkuPA+oAM1GoAICBBoVQKVKBBCosUDA99rSsz6vMWDwVl/OMxq16H3AN9lakLFOBBAIgACBRmmDRGkEEPCjQCzYb23pWdP9qOrrNgU+0EjEfO1L4xBAAIHeCAT707g3PaYMAgiEXMAslqyzIP9Lf0GgUez46bcvii3jp/wN9X5qDW1BAAEEyiNAoFEeR2pBAAEfCcSHruij1hTflBSXThWN9s77XX/1bNEV1ajAsKEhP6NRI1dWiwACtRUg0KitP2tHAIEKCMRXXq0CtVavytSH71RvZSFYUyqVMf1YXpC7MnRIkFtP2xFAIIgC1WgzgUY1lFkHAghUVSCx8uiqrq/cK0t99pGlF35V7mpDW9+nn5m1BvyLujijEdrNk44hEGkBAo1IDz+dL16AEkEQiAc80JBx+2sv6IFUgMA7H/T+q20LqL4qWYYN4dKpqkCzEgQQqKoAgUZVuVkZAghUQyA+ctVqrKai62j7z0sVrT9Mlb/3QbDvz9BYcEZDCiQEEAibAIFG2EaU/iCAgMX6N1hscLAvem/PBhqZTPCP1Fd6c9T9GR98FHynYUOC34dKjzX1B0uA1iIgAQINKZAQQCB0AvGVg31WI9PUaKnXXgzduJS7Q+99aNaeKnet1a1vQINZnz5cOlVdddaGAALVECDQqIZywesgIwIIlEsgMTLYN4TLoeXR+/RA6kbg8b8HPMrI9m30KgQZWQb+I4BACAUINEI4qHQJAQTMynZGo4aY7W+8Yqnpn9SwBf5e9bTPM/b2e8G/5GjsGjF/Q9M6BBBAoJcCBBq9hKMYAgj4WyAR8EunXN0lj/7RneQxR+CRvwX/bIa6NGY1PorlUEwiLwIIBEOAd7dgjBOtRACBIgXiIfiKW3W59f8et0zTYk2SPAIzZmXstTeDfzZDXVprrP6SEEAAgUALdNl4Ao0uWZiJAAJBF4ivslrQu7C0/a1LrOWBO5ZO87dD4NEng/+VturMkMFmAwdw6ZQsSAggED4BAo3wjSk9CpIAba2YQKyujyVWX6ti9Vez4iWP3GupGdOquUpfr2ve/Iy9PCUcgcaY1QkyfL2x0TgEEChJgECjJD4KI4CAnwXqttzez80rvG2ZtDXdclXh+UOcU78t8tu7U1bJnxipJt+Y1fkYrqY360IAgeoK8A5XXW/WhgACVRSo22KHKq6tsqtKvfemtT77RGVXEoDaH30yY+99GI57M8Q9ljMaYiAh0JMAywMqQKAR0IGj2Qgg0LNAfO0NLDZwhZ4zBiRH8x3XW6a5KSCtLX8zP5mWtgcfDcc3TUmnrs5szTGaIiGAAALhFCDQCOe4Lu0VfxGIuEAsFrO6LbYLjcKt39rHDpvxT0tbeI7oFzo4TU0Z+9VvwnPJlPq90XoxSya5R0MWJAQQCKcAgUY4x5VeIYDAMgG/3aexrFlFPcxfaUWbcMHJdv7KSfvTwk/t1GnPFlU+DJlv/X3KFiwMQ0++6cNmG/MR/I0GUwggEEYB3uXCOKr0CQEEOgSSm21rFgvuW90/d93FtjvhQHsqs8jcf7fMnWrXz3nDfRr6xyefSdvUd8N3FmfTjTib4ZONl2YggECFBIL76VshEKpFAIFwCcT61Vtyw80D16m2PnU26bTj7FvjRtuc1JLl2n/W58/ZIws+XW5+2Ga8MTVj9/45PPdluOOzzpoxG9BAoOF68IgAAuEU6H2gEU4PeoUAAiEUqNsyWN8+NW2tsbbHud+36xta8t6NoeP7Ez5+wl5dPMfC+u/FV9P2y1+3h7J7m27Ex28oB5ZOIYBAJwHe6Tpx8ASBYAvQ+q4FkgEKNO6Z+C3bYcJO9kaqsevOeOa2ZFK2zwd/temtPef1FAvE5GNPpe32O8N187cXfstNOZvh9WAaAQTCKUCgEc5xpVcIIOARSIwcZfFVVvPM8d/kwsGD7ZjzTrZTxvS3xenCj+LPS7XYNu/+yV5ePNt/nepli+68L2UP/CU0l0stp7DyCLMVhxFoLAfDDAQQCJ0AgUbohpQOIYBAVwJ9d9u/q9m+mDdl3Fa246kT7OHYNzd8F9OwmW2Lbcf3Hgz8DeLpdMZu/l27/fO5dDHdD1ze7cclAtdmGoxAuAToTbUECDSqJc16EECgpgJ99jzYrE/fmrYhd+WpRMKuPfG7tteu69n0VHPu4qKet2XSdubnz9mBHz5qi1KtRZX1Q+bPpmfs8p+127//o7tP/NCiyrQhmTDbeTvOZlRGl1oRQMBvAgQafhsRH7eHpiEQZIFYfX/ru/sBvunCjNVWtf3OO8GuGJKyVN5bvotv7l8WTrNN3r7P3mn+qvjCNSixJBsT6VulFGR8PqMGDajyKrfZKm4NfNtUldVZHQII1EqAQKNW8qwXAQSqLtBnvwlmsdofTX58371s+yN3t5fTvbtUyr751+XUp62LbLN37rPTPnvWZrU1dZnHDzPfmJq2S65osyefSVsm3CcyOrj32Y2P3Q4MJhBAIPQCvOOFfojpIAIIuAKJEaOsbuud3adVf2xqaLCzzzrBDt9kJfs63VbR9bdm0varL9+ysW/eZedMf97mtrdUdH3FVL7wa7ObfttuN9yWsq9C9mvf3TmsNSZmK4+sfaDbXRtZVg4B6kAAAVeAQMOV4BEBBCIh0O+A79Skn1M33dh2OuMI+12f6p5h0FfgXjf7dVvjzTvtwhkv2vwufvyvWiCzZpvd92DKLs6exZjyekROYXhwd9+Zj1wPB5MIIBABAd+860XAmi4igIAPBBLrbmyJtTaoaktuPeY7ttu+m9nHqeoGGd5O6itzr/riNRv1+u+dG8Z/O/fdqpzlWJztsi6N+vE17fbDyW32xNNpa1nibVk0pgcNNNuC386IxmDTSwQQ6BAg0OigYAIBBHIEQvu03/7VOavx5cgRNuGCk+2CkQnTpUx+ANUZDt0wfty0p23E67+znd970HTGQ/d1lKt97SmzV/+TsV/e1m5nXdxmutn78xnRO4Ph9dx1x4QlElw25TVhGgEEwi8QD38X6SECCCDQWSC57XiLDV2p88wyP3t6911su+/tZ09lKnPDdzmaq1+reLZxlnMPx5g377KNp95rj0yZZZ8+l7a572esZUH+tejm7bnzzKa+m7Gnn007wcQN2cDikp+02ak/aLNbftdur0/NWForyV9NJJb062u2287lCjIiQUYnEUAgJAIEGiEZSLqBAAKFC8TicavUWY22PnV26anft0O3Gm3z0q2FN8oHOd9qmW9NL9Tbq79L2TPXtNsj57bZH49vs4fPbLPHL2qzf0xut6cuy84/p80evbbdLriszX5+c7vd/ceU881Rb2QDiy/mmKWyZzR80B3fNGH/vRPW0J9AwzcDQkMQKLcA9eUVINDIS8MCBBAIs0Cf3Q8w61df1i5+tO46tse537cbBiwp4y9jlLWJPVZWP3N5k9ZGs8ZsADH/44wt+Cx7pmOhWbyyX5rVYzuDkmHwILM9xseC0lzaiQACCJRVgECjrJxUVoQAWRGoqYB+wK/+0GPL1oY/fPsQ2/mQ7eyNVHavvGy1Vreiwe19bcnsAj8WCDQKGpzDDkxYknszCrIiEwIIhE+gwE+U8HWcHiGAAAJ9DviOxVcdUxLEwsGD7ZjzTrZT16i3pnR7SXXVuvAejWMLbkKGQKNHq5VHmG23NR+zPUKRAQEEQivAO2Boh5aOIYBATwKxeML6n3pxT9nyLp8ybivb8dQJ9nBsUd48QVqw5YJRBTc33RLtb5EqBOqoiYlCspEHge4FWIpAgAUINAI8eDQdAQRKF0iOXc/67vWtoipKJRL2sxO/a3vtup5NTzUXVdbPmVf/cljBzWtvKThrJDNutH7M1l2Lj9hIDj6dRgCBDoGwvgt2dJAJBBBAoCeBfkf9t8UGD+0pm7N8xmqr2n7nnWA/GZKyVGBv+Xa6styfAV/0X25evhntiq84qdElTzJpdsRhnM3oEoeZCCAQKQECjUgNN51FoJYC/l13rF9/63/cWT028PF997Ltj9zdXk6H41Ipb4f7phLWOqO4nePBy39BlbfKyE7rBvDhK8Yi2386jgACCLgCBBquBI8IIBBpgbptd7XkpuO6NGhqaLCzzjrRDt9kJfs6Hc67oHdvWsOKPUEzsE+XXJGeueYaMds9SD/OF+nRovMIIFBpAQKNSgtTPwIIBEag/8kXmNV13nueuunGttMZR9gdfRYHph+9aejWC1YtulhD36KLhLpAn+ymc+KxCYvFOJsR6oGmcwhUWCBM1cfD1Bn6ggACCJQiEB823Oq/c0JHFbcc8x3bbd/N7ONUU8e8ck5sdM/LNuYf73aqcvhbM2yPC/7Ukba65Z+WaPnmLIryu8s17Rbuu7DZtv3FUzbos/nurKIe15y7YlH5lbk/vw8hho404aCEDRtCkNEBwgQCCERegEAj8ptAWADoBwLlEej7XxNt/hbjbMIFJ9uFIxPWmkmXp2JPLQoQFCyMfH26Z+7SyYY5i+yNo7axp646xJ6edIAzc/0/v+Y8KpgY9v5se/bCfZ2kac3TwlX+Pc1mbzLavl5tqJ4WnQZ/MaDoMvVJ7gZ30dZdK2a77shHquvBIwIIICAB3hWlQEIAAQRcgXg2uDjzh/ZCTF+r5M4s7+Mnu63nBBJfbDp6uYq1bM5GS3/PItWvzuatM8L6Zc9W6KyGgor2vknzJs1TUtAxc8vVl6uvkBnxbLzQ9nmykKyd8vSN93D0vlPu8D7RJVMnfDcR3g7SMwQQQKCXAgQavYSjGAIIhFdgvX5D7N6xe/migw2zv7aWwfWmoGNJ9jG5pN28SfN0NkMBiaZ70+idmlazTKr4ktyisdTshKMTtsJggq6lGvz1uwDtQ6CaAgQa1dRmXQggEBiB/xq8ul08couatlf3awyeNs8+2G9jpx0KJBRQ7HTlY6akaS1wz2bong9dkqWkslpWSNpuwWqFZFsuT132TMhyWJXiAAAQAElEQVRyMyM2Y/+94rb5JnyURmzY6S4CCBQowLtjQVBkQgCBKApcvso4233g0suYqt1/BQob/HGKvXnENqYAw12/Lq3S/RtKmnbPZujyqWRLm3Nfx8unjLexT71jmueW6+5x7Xkrdbc477Ko36Kx0foxO2g/PkbzbiAsQACByAvwDhn5TQAABAIqUIVmx2Ixe2DNvW3VuuJvlC6leW6QMeX7O3R7c7cCCfdsRr+vm51LqrReBSbt/eo0WVAaNntQQflyM8V6cblVbh1BfT5iuNnJ30vwVbZBHUDajQACVREg0KgKMytBAIGgCgxO9LW/rr2f1ceKv1m6N31WkLH2I2/YC2ft0W2QobrdsxkKLFoG1Ts3iWu+AhCd3dB0ISn1WeFBibe+WLv3WXSm6/uZnXlS0vr15b4My/nHUwQQQMArQKDh1WAaAQQQ6EJg4/phdteY3btY0rtZ3q+3XfOJt238pL90/P7F8DdnWP2CZuceDN1roeRd7q5Rv5eh+zc+235NZ5a+1nbh6sNs12xd4256xj7eY/1Ol1w5mbr4s1XzKpZu7WJBIbOWFJIpfHlOOS5hK61IkBG+kaVHCIRSoKadItCoKT8rRwCBoAgcMmSs/XTUdmVpru6t0D0Wbnpm0gEdZy/eOnyc89W37jI9epe7DVBg8Z/v7eB8G5U7z1uv+xW57rJ8jzsu7N2N4Kov1Rq9u8EnHpyw9dfho1PjT0IAAQR6EuDdsichlkdTgF4j0IXAuSM3s0tW3rKLJcGdtd684b1ufKql10UDWfCQA+K21658bAZy8Gg0AgjURIB3zJqws1IEEAiqgL6J6tSVNgpq85dr9/A5g5ebV+iMtsWF5ixPvlrWoiBjvz0StWwC60YAAQQCJ0CgEbgho8EIIFBrgRtX28mOX3H9WjejLOuPTetTUj0DencfeUnrrHZhgoxqi7O+AAnQVAS6FSDQ6JaHhQgggEDXArettosdNXSdrhcGZO76S4ZZe3NpNzUP7BeQzvaymQQZvYSjGAIIIJAVINDIIlT9PytEAIHAC8RiMfv9GrvZQYPXCGxfdlk4puS2DyjthEjJ669kBQQZldSlbgQQiIIAgUYURpk+IoBAjwK9yRDPBhv6Qb+gXka14VcjetPtTmX615V2RqRTZT55kh1WO/aIhHFPhk8GhGYggEBgBQg0Ajt0NBwBBPwgkIjF7derj7frRm9vQXtDXXnOCiUT9q/O7xiW3M5CK+jX1+zc0xK24zZBG81CexiofDQWAQQCLsA7acAHkOYjgIA/BM4asak9tOa+1jcWnG8mqptZ+g0W/UJ0QmPYULMf/SBp667FR6M/XlW0AgEE/CdQXIt4Ny3Oi9wIIIBAXoH9V1jDnlv3WzY8WZ83j18WjGxvsNb5pUcJfUqvwhcka64Rc4KM4SuFpEO+UKURCCAQdQECjahvAfS/KgKsJDoCWzasZFPWP8zW7zfE153ee+FaZWlfGK6c2mFczM47PWEN/QkyyrJRUAkCCCCwTIBAYxkEDwgggEC5BEb1GWCvZoMNP98kvulXK5elu3WZslRTi0qcdU48OGHfOzJpiQRBhgPCHwQQQKCMAgQaZcSkKgQQQMAV6B9POjeJP7zmvjY40ced7ZvH0XOGlqUt8fayVFP1SlYZaXbpeUnba1c+BquOzwoRyCvAgrAJ8A4bthGlPwgg4CuBA1ZYw97d8HDbbeAoX7Wr36zSbwRXh2Ip/Q1O0lfX7rtHPBtk1NmqoziLEZyRo6UIIBBEAQKNII5aTpt5igAC/hYYWdff/r7OgfbzVXewvrHav+0Obu9rS+aUqR1L/G3vbd2I4WYXn5O0Qw9IWCI4Xw7m7QLTCCCAQKAE4oFqLY1FAAEEgiHQZSvPHL6Jvbb+RFuv3wpdLq/WzL0a1yzbqtKt/r9JQ2cxdInUjy+oszVW5SxG2QafihBAAIEeBAg0egBiMQIIIFBOgfXrh9gbG3zbObsxJNG3nFUXXNeW88t3GVeqpeDV1iTjxhvEbNL5SdNN30nOYtRkDPyzUlqCAALVFiDQqLY460MAgcgL1MXiprMbH218pJ2VPcvRJ/u8miirzS3PjeBqc2uj/vovjVk9ZheckbAzTkraqJU5i+G/EaJFCCAQBYEeA40oINBHBBBAoBYCOqNx3ao72HsbHWGHDRlbtSYM+KKhrOuqrytrdSVVtsrImJ12fNIuPjtpa43lI64kTAojgAACJQrwLlwiIMURqIEAqwyZwBp9Btr9Y/e2F9c7xLbsv1JFe9c3lbAlM8v71j+oNleAdXIaNsTs+0cm7McXJG2zjTmD0QmHJwgggECNBMr7aVOjTrBaBBBAIAwC2zSMcH7o76l1DrADBq9ulXiD3nNx9sxJme/fbqjhz4SMyMZlR387YVdPqrPtx8VNN37XZltgrQgggAACuQLx3Bk8RwABBBCorcDuA0fbw2vtZx9vfJTpXo5B8fJdm7T1V6PL3rmGZHXPICiY2GKTmP3gtKT95JI622V7PsrKPqhUiEAYBOhDzQV4d675ENAABBBAoGuB1fsMdL6dauamx9ovV93R1uk7uOuMRcwdM3dYEbkLy1qtezQGNJjtt2fcfvrjOjvluKStu3Z1A5zCNMiFAAIIIOAKEGi4Ejy6AjwigIDPBBriSTtt+MbOTeMvrHeInT1iU1utz4BetXLw7IG9KtddofoKf23s6qvGnPsvfnZZnR2yf8JWKD3e6q47LEMAAQQQKJMAgUaZIKkGAQQQqJzANzVv2zDCrh29vU3b+Gjn5vFigo54xqxtevmjgnLfopHIfjKtlz1bMfHghF1xcZ398Nykc/9FMvmNA1MIIIAAAv4XyL6d+7+RtBABBBBAYHkB3TzuBh1vbfhtu37VHe2gFdawwYmud/3HL17DMqnl6yl1TtdrK67WQdkTLTtsE7f//n7Sbriqzs49LWn6Ne+Rw4urh9wIVE2AFSGAQI8CBBo9EpEBAQQQ8L/Ahv2G2unDN7YH19zX5m32fXtpvUPtylHb2B4DR1t9bOlZjG0Xlv9GcMkk0/pbeEpkm7PGajHbdae4HXfU0rMW111RZ987ImFbbhqzvn0Lr4ucCCCAAAL+FYhXuWmsDgEEEECgwgIJi9m4huF2wcgt7Ml1DrCmLU60f6xzoO05epSN3ipug0fHrIxfZGWJjOX9V19vNnqVbHu2iNu3v5WwC89M2q3X1dkl5yTtyMMStt3WceOsRV4+FiCAAAKBFiDQCPTw0XgEyiFAHVEQ2HXgKBu/xVDb9uSE7TkpaYfcXGf7XV1nO2d3+DfLnklYc7e4jdgwZkPHxmzwqjEbuLJZ/xXN+g02q+tv1tXVWMm+Zg0rmQ0cFrOtNo/bvnvE7Yhs8HDmyUmbdH7Sbvppnf3yqjpn+sRjsusdH7c1x8SiwE0fEUAAAQSyAvFs4j8CCCCAQAQF+g8zG75+zNbKBhmbZ4ONnc5K2m4XJW3PS5O29+XZQCQbJOx/bZ0ddEOdfeumOjvs9s7p4F/V2b5X1tlOxyXs5GMTdugBCdttp7htlK1TZzH69Ikgarm6TD0IIIBACAQINEIwiHQBAQQQQAABBBBAoLIC1F68AIFG8WaUQAABBBBAAAEEEEAAgR4ECDR6AGJxqQKURwCB3gg8/fxrtuH4Y02PPZWf+t6nNuHESfblvAWdsv7mnkedOlSPpt2Fyqf8KufO4xEBBBBAAIFyCxBolFuU+hBAAIESBRRcnHbR9T3WooBhnyPOs4knTbJFjU2d8mvZcy+/ac888AsnaVrzlOnhJ563fXYdZxuuu4aekhBAAAEEEKiIAIFGRVipFAEEEOidgM4y3HTHQ3b/bZNs1VW6/7W6lYatYI//4Rq779ZJNnBA/04rnDN3gQ0YUG8N/fs5SdOap2BDQceBe23fKT9PEECgtgKsHYEwChBohHFU6RMCCARSQEHG2ZN+ZZPOOdYURJTSieErrmCNjc22uKnFSZrWvIezZzN2GLdxyfWX0jbKIoAAAghEQyDggUY0BoleIoBA+AV0pmHStXfYdZNOLcslTQpUFFCMP/RMU9K0FN2zGRdd+euO+zd0qZaWkRBAAAEEECinAIFGOTWpCwEEzDDolYAua5o+c45zv4Vu3lZw8Hn2ue7V6G0gcNzh+9nUZ+5wkqbdsxlaV2NTs73y2C3OZVe6VEuBTq8aTiEEEEAAAQTyCBBo5IFhNgIIIFBNAd2Y/cIjNzlBgYID3cStezRunHyG7br95iU3RYGEezZjzryvnMuqVKkupxrYUK9JUogF6BoCCCBQCwECjVqos04EEECgFwL6ilp9y5SChmKLu2czdEnV8GFDnBvFVYfObixa3KxJEgIIIIBA9QQisSYCjUgMM51EAIEwCijgUOChr7d9+/1PnXsxdO9Fbl91k/nrb39ohx+8m7NIZ0823WAt23rfk51LtU459iBuDndk+IMAAgggUE4BAo1yalJX5QVYAwIREdCZB311rfeyKd1noXlaJgY96rkutXLT5AtP0KJOSYHFDZefbv3r+3XMV11uGe86OjIwgQACCCCAQIkCBBolAlIcAQQQiLoA/UcAAQQQQKArAQKNrlSYhwACCCCAAAIIBFeAliPgCwECDV8MA41AAAEEEEAAAQQQQCBcAgQa3vFkGgEEEEAAAQQQQAABBMoiQKBRFkYqQQCBSglQLwIIIIAAAggEU4BAI5jjRqsRQAABBBColQDrRQABBAoSINAoiIlMCCCAAAIIIIAAAgj4VcCf7SLQ8Oe40CoEEEAAAQQQQAABBAItQKAR6OGj8aUKUB4BBBBAAAEEEECgMgIEGpVxpVYEEEAAgd4JUAoBBBBAICQCBBohGUi6gQACCCCAAAIIVEaAWhHonQCBRu/cKIUAAggggAACCCCAAALdCBBodINT6iLKI4AAAggggAACCCAQVQECjaiOPP1GIJoC9BoBBBBAAAEEqiRAoFElaFaDAAIIIIAAAl0JMA8BBMIqQKAR1pGlXwgggAACCCCAAAII9EagTGUINMoESTUIIIAAAggggAACCCDwjQCBxjcWTCFQqgDlEUAAAQQQQAABBJYJEGgsg+ABAQQQQCCMAvQJAQQQQKBWAgQatZJnvQgggAACCCCAQBQF6HNkBAg0IjPUdBQBBBBAAAEEEEAAgeoJEGhUz7rUNUWm/EVX/to2HH9sp/T9s662puaWshv85p5HrVJ1l72xORXKQ21XH7yLnn7+Ndtu/1Ns6nufWq3+LWo0+7/n0wWn197I1KqprBcBXwmkPvvIWp96uOCUen+qr9pPYxBAAAGvAIGGV4Np3whss/n69spjt9jUZ+5wHtWw0y66viLBhuoOU9p1+83thUdusg3XXaNm3Zo7P2P/e2+q4PTok6mKtTVfQFaeFRZWi4I/BdCF5a5sLrVlnyPOsy/nLSh6RaWULXplOQXccVRwremcxVV7Kjf5VWo82994xZpuu6bg1PriQWcR/gAAEABJREFU01XreyEr0thojHIPgBRStlp5tB3rYJYeq7XOcq1HB5B0IMnPvuXqK/WEQ4BAIxzjGOpe9K/vZ0dP2Mtmzp5ni5uWntXQB4Q+7PWh73Zeb8B7H/6DjiP5bh496o1ZHyx6VD6V0fzrbr3PXnrtHdt635OdMyh681adqlvLlU8p98NTy9w8qlN1q6zyuuU1T0kfuiqvZV0lN/8jT73gnF1xy3w2/QvTOvTcTd6dmyt+cafTdvVBy9UO9c1tm+p116d5yuMmt63u8jA85vbR9fBD357856u25y5b+aEpgW3DPQ/+w0YOH2q//fn5pveEYjui18OEEyd1vD+45fWa0bai10bua1Wvt9zXykrDVrA7f3mRffDJjOXqcusMyqP7msnto57LQ49+6ovGyvsen7dt3SzQdnDTHQ/ZfbdOMh2U8WZV/dpGlMc7v5BplVFZ1eHNr/d+bVfyVMo11Tam+UoaD7es6lF9qtedp0cdQLr92vPs8adf7tXBAtVBQqCaAgQa1dRmXb0W+HjarF6V/XzmHLvz/ifs7/df55wdOf7I/e2cH9/kvEHrQ+bskyaa9+zJcYfvV/B6cutWWX0oHP0/k+380w531qczMto5KuRszE+ygcM5J3/bKaedqfpsgLXFxms7Z3RUzzMP/MKmvPmBuR9Ul5x5tNN29UHL853F0IfXRZN/7XywKp/quf8vz3TUU3CHfZxRJt4+qp/6MH7ob/+qeau1TXwxZ75ttO6YmrdFDdg1e8br8T9cY9ph1vMgJBlqx+rIQ/bstrnKp+BcgYN21JTZ3dEbf+iZNj37fqB53nT3n560yRed4Lzu9FrVgQctV/nGpmY7/ODd9LRTkt0pxx5kKttpQQCfjF55pU47rTKUteYHsDs9NvnhJ563tceM6nTGV33WdjPxpEm2qLGpxzq8GXravnRASNuV3pPc9169J6sOdxvT2XsFPn9+7NmOs/batrSNaVtTXm9SsLHPruNMffHOZxoBPwrElzWKBwR8K6APAe0YTzhgfNE7R6uuMtyuvuSkjiOgB+61vdPPt977xHks5U9u3apLb/wKDrQzp+dK2jma8cVc++SzL/Q0b9LOjj5A3Az6gJl84QkdbddzGXz06Qw3S0GPOpq+6w6bd3ywuvU89/KbHR9qBVXk00zu9pHrJ8uLTj+qy1brg15HEJW0g6E6lFGPeq75SjoSqR0JLdNOgY6mXn/7A87ZLy3X0Ugt6y5pW9OOhtzz5VM9qs9N7nrddepRZdUWLVP79VzztVPtllM9mq+Ury8qqz5qufLlS8rn1qt1aF25ed32uPkU8Ll5lF/l3GXetimf+nH6D2/osNT63LK5jzIc2FBvY1YbmbvIee6uS0G+zjZ4g26d/VDgrp280dn3A6fAsj9qv4KJ4cOGOHPWXGOUuQc1tKP3rX136nj9ORk8fxQ4KoDsydFTxJeTo0auaDuO27hjp1XvYXqu+d4Gu8YaT42rnrvLZaBtSst0dtgN1tzlGlstU1I+5dey3O1A24iWKY/yKmk70TgpXXvLvTZ91pemgEDLVK/q0aOeK6ms6tD83KQ69L6Xe3ZRr00F39rZHzigf26xbp93t32pHTrzpc8AVaL16PNB78l6PmfeVzagf72zjQ1fcQVrbGx2ztrLVtulDoIpX1dp2y02sNff/jAU7+Fd9Y954REg0AjPWIaqJ/qg0geWPjgOPOYiu/bSU0xnDErtZEP/frbKiGEdOxOl1pdbXkHAQ397ztRuN+lD8etFi3OzFvDcnEsz9KHu1qXLpLRzow/MQipQPuXP/WDVh9SixUs/1Aqpx895tBOq9mnHT489Je2UeM9+KHj7+W33O8VUl3s2Sjumulzvngf/4SzTH43j7C/nO0e/tVxnmFSfluVL2qnI9ffm1c6W6lF9OuqpM1Te5d1Nv/Xux6YzNyqnnaSnn3vN3PaoT+qbliltttFazk5Md/W5y1SH10jrcC9bdPNo29KZOgVRql/t15Fw7SQpT3dt03K9xrUjr7I3Tj7Drr7xHudMo5blJu3875DdGdZOnXeZ7PTa0A6ozlpqZ1E7c9483U2rPu3oaYdP+fT6Hbv6ys7rrqcdPb2XqIy2GT0GOR209w6mHfBPP//CZKDn3v5oh1lnghXMa7z0ePw51zhO2g7Ov+JW0w60lunovHcHOXdb0japbdOt37sd6MCKPLt6DWqsdMZXZ1q0rWtdOqDTU/3uevTobsOFvleoTClpztwFpvcMbx0KZvWeLDcFuNrONK28AwbUm7arnoJc1afAZNbs+T0ewFJeEgK1FCDQqKU+684roA8qfWBp52XwoAGmHQm9Gect4KMF+pDWh6A3eY+wFtpU7UQpSNGHultXMTuhha4nDPkUPOoDupC+aMffe4ZHQZeOOmpnSjsuSqpHO6zaedKOl54rDRrYYEcuu3xH69N6NT9f0k636s63Y6NtWjt42vnS+vLVk2/+tw/areNMlY72r7/26k4QrXq1M+Nt++nHHVrwGcFcI50dGrf5ep2aoTN0OlPneqj9uiTlxSlvO/nytc1ZmP2j17hSdtK0wzWwmyPJ6ocCAOXtKinY0Y5oV8t6mqf2K6hSwCIztcnd0dNOsOYr6Wi7U9eyP1qf1rvsaaAftNOqQO7cy2427QjrubdD2vnXc3c7lpG2NY117nagfN6Uuy15X2/Kp7qUNK2k15+SprVN5b4GNd+beqrfm9e7M++dX8lpvWfkerrr0+tKga4Oqum9XoG3PBV8yFpnc7Tt6WCT3kvccnrU+8/KI4ZqkoSArwUINHw9PDROHzQ6m/HOB9NM17qWKqIjWjNnzzN3p8V9LLVet7w+pLWzoh09d15vH7VzpcDC/dDNraeQHR03jz6MveW1g6BLUfRh5Z0f1GmNqca20PZ7zzrpA957XbZ2KPXhrqR8hdbZVT45a+db23FXy915vd0OtZ25OyPaWdGOserUuOuSQZ0pUT+UFLhqmTdp50U7MVqu5M2jbdmbt6tpHa2Vn8oqeb3yta2reno7T2c5FYTrjJHWv08vvk1LO3s6EKB6dImVd0dP18zr6LkOeuh1raPnvW2r38spAIjFYuZeXprbXgXV7vuFti9vkNXdzrTq0Xah8VHS9uJ9vWl5bir2NVhs/bnr6+65xlztzk1qY3fl3GV6jSjAcZ/nPuosjrY9Jb3Xu0GugjsZa74ONgXpYFtuH6vxnHX4V4BAw79jQ8uWCWhHQG+0+jBxd4R0BHTh142mN2Nl006N3oj1pq7n+ZJO2esD03sETUGMdi7cMvowVR7vzrkun3F34tx8XT3qw1r1eYMitU0fSjpi3lWZ7ubpaLfKK492Cm+/+xFNdkrePJ0WLHuinTBdUqPymqV26J4XHcHUDoPmBTnpyJ/a724Lmu4pKYDTB7ib3EtuNE7aodSOpZbp7FRPdeVbrnHT2Mg/Xx53vi4NcqcLfVT93kuX1Gbvdq3gRv1SP7SzfN/DTzuXunjr12vL3clWPu24u8sV6LrT+R51z4POOqqsm1RHT23LV1+++Qp6ujPSDprWr/szdJ+Ggid3e89XZ7757o6els+aPV8PzjX02ulznoT0j7aF+2+blPes18zsARo3mNf46nXiUuh9t7ud6XyvN7e897E3r8FC69eZBfc+CO86u5t2ty1tX96kAKG7clqm9SkI07Sb9LrSttTVe6+2WZ3N0OvYu73r827AgHq3Ch4RCJQAgUaghqucjQ1WXXqz13XcukdBH0T6UFTwoR0tHWnafcLZtvf4rS33TV3fDKVvm1EeJX04qh73TV716jIaHWXTcgUyWpZ7NLhxcbPzDU89qaldD/9+svPtUKpPSUeaVU47fnosNF1y5tFOVpVXPQqkvnPw7s48989ZJ05wvvZXefLtXKmPsnL7KA9dqqMdQreeID/KVf3R5S/6oHb7ounJN9zlPu141I6/AjYtd2fe8JsHOu4PcHcCcnem3LyFPip41X0wbiDUVTltawr4FJBofUqadvNqRyUWi9mceV85sxTsKjlPlv3RTrgmtROonUFNqx7dtK5HPS82ycgbnOqorpK3Hl2qpbNiCt7d+TK996Fv7mnpqm1u3mIedcZHLj31R9uCgisFT3otFrMO5VX73R09HXBwL03RevXeoTxuUsCued2Nr5s36I9uH3WjuPqibVAHVHRgRduBbhzX2Tst03avZZpW0rbU3etNeXJTvtdg7utB5YqpX2OqMsUclFD+3iZtjzqjqeBVdWib0VlGtVnPc5Py6fIpvS9om3eX6/WvAMl9rkc561H+eiQh4FcBAg2/jkyE26UjRbqEQW+2XgbtMOuIkpZrvvtc87RjoWvC/3bPTzuuWVcefTOU94hrV/WqPtWh5O586wNCOyyap3TG8Yc639/vLte6tVz5tB5v0jwtUzk3aR3ePN5pN7/q9M5X/9Vetw5Nu+3QMuV1yyqPDLRzpXq0fi1THiXNUx43uf3QsjAk9ccbTCkw082qB+2943Ldk0VuXu0wykvX62sHW+UVvOoehOUqKHCGdrz0FZSqt7si7tenKljcPRswe49cqqxujHUDal3Ko6Odqk/bgH5fRsG32qsvTchkMlrkHIHXTeuqU8sUZF7wP0d0em04GfP8kZG+ClrlVF5BnI6qerNr/QraZac8SjLfaL2xzvrztc1bR6HT2tFV0ObuXLnlnIMO44/t9OULaoc36FaQoMvLFGC//f6nzjcW6bnmu/XoUc9vvethO+moA532q3+aVp/kqJ1fuSivknZWNU9jpOdhTuqjLmFVwCBfbQ/6ggC938hJN2m7yyZde4fpTJfrIbN8rzc3j/exu9eg2qGDCu7rQcFvMfWrrQrsvWertW4FALrsTtu7thFtK9q2tKynpO1G25PKqKzq0HPNV1kdMHJfI8qj9qvNWuZN6oueu8vcR9dbxmq/8ijp/WXTDdZytlU9JyHgVwECDb+ODO1CAIGiBPTB7AZSenQDL304K0hTMOJWmJvXDQS146RybnkFrt5leq48qsetV3XpuTdpJ0NfPakjvt75XU279bjr1M6DN5/q1zKlGy4/3Ql4NU959Kj5Smq32uf2U+3WfDcpr1smNxDV/Nyketyyqlv9Vh3est62K6+bT3Upr+Ypab63bapbY6Lyyqu6u7tsRzuYCtp0xFf53ZTbR61LSetTncqndWhdmu8mPdd8LXeTnsvXLaf5mlZdKqd1aZ6SdkxvuuOhji8G0LwgJo1RPgvN1zi5/fJayETPu1qmcVTyltV6ZOgm11J5tB7Zd1WX1qPtxs2vPCrj1qN6NU+P7jw9evNruTfpHhR9QYPOXrnztX1pu1ZZN3nrcPN19ai2qw9uOT3queYrvx71XPOV1H7Nz03qQ+469Vxl5OD1Vtv1DW/qS249PEfAbwIEGn4bEdqDQAgEBg2I2c7bxwtOm28Srrci98g7lzWUb2PWmR8dGfYeLS5f7YXXpCBD96nZJA4AAA0mSURBVIHokhjvzl/hNXSfM7HamtZnjwMLTsl1Nuy+QpZ2ElBQccqxBzlnttyzCJ0y+PyJggydZVPgrb74vLk0DwGr0Kc7sgjUXkBHiHSUijfj6o/FsKFm3/12ouC0357heivSDqiOjutoZrH6OuKpI6C9KVvsuoKUXx5yUdJ0rdqu9xO9r+hocyXakNxka+t/4nkFp7ptd61EM0Jdpz4bdKZAj0HrqN5bdIZD7xNBazvtjaZAuD7dozmG9BqBwgTIhQACCCCAAAIIVFGAQKOK2KwKAQQQQAABrwDTCCCAQJgFCDTCPLr0DQEEEEAAAQQQQKAYAfKWUYBAo4yYVIUAAggggAACCCCAAAJLBQg0ljrwt1QByiOAAAIIIIAAAggg4BEg0PBgMIkAAuUR+LK92X499+2C04MLPinPiqmlkwBPgiewcHrGPvm/dMFp/sdLf6QxeD2lxQggEAUBAo0ojDJ9RKDKAp+2LrITp/2z4HTlrCkVbaG+L1+//KvfQOhpRb+551HrzW816JeEC11HT23w03KZqV/qn5/aFda2zH47Y//+31TBafqr6bBS9NgvbZPaNrWN9pi5RhnUNrVRbV3WBB4QiJQAgUakhpvOIhA+AfeDXMGEt3f6YFfyzqvUtNatH5P7028uM/3Ognc9aoOCF+88TavMhuOPNSXtiKgfml9IUp0q5ybV5S2n9bnLcoMm7zJNu+W0/gknTjL9IJg7T4/qz52/vMg++GTGcsu0nBRsAW0D7rbiPuZuT37sodqY7zXXU3t7ev2obtci97Xp9dK0uy5eP64Ejwh0FghGoNG5zTxDAAEEihLQD3PpR9a001xUwQIyNzW32J33P2FHT9jLvD8kp50Q7aw89LfnlqtFO/NX/fIPdt+t2R37Z+6wCQeMt/OvuNVU13KZc2Zoh0azXnnsFtOPjqkO1aU6NV87Sff/5Rl75oFfOMtHDh9qV/ziTi0ylX3u5TedZVquac3TwoefeN722XWc6QfB9Nyb5HbKsQfZ3X960jub6YAL5G4r7vb0wqtTfd0zvU66es0V0mh3e8/3+tHrSK8nva7k4X1tqqxeM3rtKGla87ReXj9SICGwvACBxvImzEEAgQIEgpRFO1TeI5PaOdBzBQJuUmDg7dMFk29zzjZoucp7l3mnP/nsC1u0uNk2WneMd7Ydd/h+zo7+QXvv0Gm+nrw45W3bcpN1Onbqt91iA5vxxVxTXVreXdJO/+QLT+gIaoavuILFYjGbM+8rp9iT/3zVCVyUTzP23GUrm/LmB06QMWfuAhswoN4a+vdzkqY1Tx7aaTpwr+1VpMuk/ukIsvJ2mYGZgRPQtrLFxmt3OgunQPOi04/q6IvG2/tacV8Ler3kni3TmQIlFVYwoOV6/Sgpv+a79V1/+wO23f6nOEk791qmPMqrpLKqQ/Nzk14nXb3m3PJuG3PL6bleF929frp7beq1otcMrx9JkhAoTIBAozAnciGAQEgEtPOiswfawdIRSx3Z3Gbz9Tv17qXX3rFv7buTEyjcOPkMu/rGe5wd9U6Zlj3RDv7aY0Z12llbtijvw0efzui0LDdY6LSwhydvvfeJLfy60YYPG+KcEVEw4C2i+ZlMxrSTpPU0Njbb4qYWJ2la83Q0dodxG3fbB+1cqV6tT48kXwoU1ag11xhlOuOWb8dcQcHR/zPZCVz1WtFR/pvueMh5LSgo1c6+dvq1UuXV5XVHHrKnsx2edtH1prNpKqej/48//XKnS+9en/qh/f3+6+yFR25yAm4FCbln4u558B+qermU7zXnBvcKoBSs5OuXt0Jtz+7rR/O7e23qtaLXDK8fSZEQKEyAQKMwJ3IhgIDPBbRjo50LN2kHqqsma8dIZw+0Q9TVcs1T4KGkae2oDxzQX5Ndpo+nzTLtsHW5sJuZvSnjrU5HgXVEWP2efNEJzs6au3zs6iu7k50edTRXAcX4Q880JU0rg3s2Q0ejXb/cnTRdFqYdR+UnhUNAO+ZnnzTRtA25466zFwoa1EPthOtRQYUex6w20gY21Jvma1tSgK0zAFqmeVqmPLmvsdy8yu+91FDBv7ZBXaakvFquM3Gap2V67k09veZ0xkIHEHR5lV4jeq14y2ta87RMfc99/eR7baptes3otaOkadWldsqI1480SMsLRHsOgUa0x5/eIxAaAZ150NFTN3V1yZLb2UEDG0xHJ62G/3KPnLpNcXeA3B0/9zH3UhJd4qKjwTparDMu3sBAO2JufbmP2rl0jTTtns3QGY/GpmbTDpr3yHVueZ6HS0DbgLs9aOxXGTGs0/1Cn8+c4wSl2g633vdk09k+V8AbDOgsgna8FZBq+deLFtvEkyaZyinlC/yV103X3XpfR34FAO783jyqHb/9+fnOWZNrb7nXvK8P1dfd6yffa1PlvF6a5vUjFRIC+QUINPLbsKSGAqwagUoKaCdIO9blWIfOHnS3Y9LVOnKPmKoturxJZ0/cHSB358991E6Tdp5y69NRVl0GpuBCy3PPOugyk1gs1mVgpSPX7tFY5dNlIapfQZiOTmuaFB0BbT860zBz9jzn0jr1XGf2FIC426Eed91+cy1y7kvS5VMKPnTJno7qOwuyf0avMtz50gHld5N2zLOL8v7PPViQb5sv5DWnMyEKznefcLadc/K3zW1z7sq9rx8t6+61qeXexOvHq8E0Al0LEGh07cJcBBAIqYAu7Rg1ckVzL/nQZR7vfDCt171VcKBr07XTUWgluvn732+833HNutqiNqltPdWhMx6Tb7irI5ueP/3ca6adL83UUWZd6+62R0eaFYhoh0rLvck9Gqtl6seAAfXOYgU+2oF0niz78+W8BaadSd0UvmwWDwEX0KU+uUf6vduLxlpBh/deCeVXUte13eibyn71uwed+zH0XPO1HStQ/flt9+upk7Sd3vtQ1/dcKMDR2RCdmXO3WwUKumFcj04Fnj/aVrt7zalfOvuioEln/RS8u8XVju5eP8W8Nnn9uKo8IpBfgEAjvw1LEEAghALaqdERztvvfsS5TGPStXeYjr66O+rFdtndqdI16t6yurnVvWREl4ToenDt5CiPdnwu+J8jOi4tUWBw9SUndXyTlPLkS1rfh5/McNqu+nV5iq4xd4/Y6lHXuusaci1XcHDJmUcvV53a8vrbH9rhB+/mLFObNt1gLdMOmurU19m6O47KoP7pbIl3nuaTgiugoFSXKGk7cZN6o3sc9Kix1m+oaPt0lysYUACi5UraMZ8x60tTXXqupNeYzk5o23PLHX/ONbbRemO1uMuksx3e7Vbb4YCG+i5fE3oNKJDRNumtzH3NqS06i6LXgne5plW2u9ePXgeFvDb98/pRr0gI+FeAQMO/Y0PLEECgAAHtDOk3MnJ3KrSzpKQqtEx5lFfPtTOhI53aGbnpyjM1y/nWJk1oh8d7yYby3n/bpLzfyKSdKh051Y2n3qOvqkf1u0nrU11ah5La5C7ztk3Luktan9rnltWj6vKW8a5beVXGu1zTassNl5/eaUfOW85bp44y69uGuruBXnWSgiWgMdb2403ua8btiV4z2j7dPJrWPHe5tiNt26rLnadHbXPa9txyyqO8Kqs6cvOrjHf7Uzk91/zcpLq7e811Vbdbh8p626X15ObXc81XUlvVZre8+6i+8PpxNXhEIL9AJAON/BwsQQCBcggMT9bbiStuUHA6eEjn36AoRxu6q8M98qmjrTryr6P32nHorkx3y7RjoqP9hxz3I+erP7vLG7RlCjL0Faf6hqFSjILW71q1d4XRMRuzc7zgNHRsND/Gg/Ka4/VTq1cS6/WLQDTfofyiTzsQCKnA6n0G2q2r71JwunDkFlWV0JFSHa10k3ZaSm2AjgTnO/pZat21LK+jueqX+lfLdkRl3cM3iNmW300UnEZvFYsKzXL91DapbVPb6HILfTJDbVMb1VafNIlmIFBVAQKNqnKzMgQQQAABBIIiQDsRQACB0gQINErzozQCCCCAAAIIIIAAAtURCNhaCDQCNmA0FwEEEEAAAQQQQACBIAgQaARhlGhjqQKURwABBBBAAAEEEKiyAIFGlcFZHQIIIICABEgIIIAAAmEXINAI+wjTPwQQQAABBBBAoBAB8iBQZgECjTKDUh0CCCCAAAIIIIAAAgiYEWiUvhVQAwIIIIAAAggggAACCOQIEGjkgPAUAQTCIEAfEEAAAQQQQKDWAgQatR4B1o8AAggggEAUBOgjAghEToBAI3JDTocRQAABBBBAAAEEEDCrtAGBRqWFqR8BBBBAAAEEEEAAgQgKEGhEcNDpcqkClEcAAQQQQAABBBDoSYBAoychliOAAAII+F+AFiKAAAII+E6AQMN3Q0KDEEAAAQQQQACB4AvQAwQINNgGEEAAAQQQQAABBBBAoOwCBBplJy21QsojgAACCCCAAAIIIBB8AQKN4I8hPUAAgUoLUD8CCCCAAAIIFC1AoFE0GQUQQAABBBBAoNYCrB8BBPwv8P8BAAD//70igKwAAAAGSURBVAMAykPB0Il8Yk4AAAAASUVORK5CYII=" }, "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": 58, "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": 59, "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", "