{
 "cells": [
  {"cell_type": "markdown", "metadata": {}, "source": [
    "# Volatilité rugueuse avec trajectoires dépendantes\n",
    "## Projet de fin d'études -- DUFE 2025--2026 (sujet 5, encadré par G. Pagès)\n",
    "\n",
    "**Auteur :** Lionel Feliho  \n",
    "**Date :** 2026\n",
    "\n",
    "Ce notebook accompagne le rapport LaTeX. Il reprend pas à pas la construction du modèle, la simulation et l'étude numérique. On suppose le lecteur familier avec le **calcul stochastique d'Itô** (Brownien, intégrale stochastique, formule d'Itô) mais on n'exige *aucun* pré-requis sur la volatilité rugueuse.\n",
    "\n",
    "**Plan :**\n",
    "1. Rappels motivants -- pourquoi une volatilité *rugueuse* ? Pourquoi path-dependente ?\n",
    "2. Concepts clés : mouvement Brownien fractionnaire et équations de Volterra.\n",
    "3. Le pont elephant $\\longleftrightarrow$ goldfish de [BCGP23].\n",
    "4. Schéma d'Euler hybride et étude numérique de la convergence forte.\n",
    "5. Modèle quadratic rough Heston : trajectoires, feedback et smile.\n",
    "6. Structure par terme du skew ATM : la signature rough.\n",
    "7. Discussion, limites et propositions de recherche.\n",
    "\n",
    "**Références principales :**\n",
    "- [BCGP23] O. Bonesini, G. Callegaro, M. Grasselli, G. Pagès. *From elephant to goldfish (and back) : memory in stochastic Volterra processes*. arXiv:2306.02708.\n",
    "- [GJR20] J. Gatheral, P. Jusselin, M. Rosenbaum. *The quadratic rough Heston model and the joint S&P500/VIX smile calibration problem*. Risk Magazine, arXiv:2001.01789.\n",
    "- [GJR14] J. Gatheral, T. Jaisson, M. Rosenbaum. *Volatility is rough*. QF, 2018."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "%matplotlib inline\n",
    "import os, sys\n",
    "sys.path.insert(0, os.path.abspath('..'))\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import gamma as Gamma\n",
    "\n",
    "from src import (QuadRoughHestonParams, simulate_elephant_goldfish,\n",
    "                 simulate_asset, fractional_kernel, phi_burst,\n",
    "                 smile_from_paths, strong_error_study)\n",
    "\n",
    "plt.rcParams.update({'figure.dpi': 110, 'font.size': 10, 'axes.grid': True, 'grid.alpha': 0.3})\n",
    "rng = np.random.default_rng(2026)"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "---\n",
    "## 1. Pourquoi parle-t-on de *rough* volatility ?\n",
    "\n",
    "Dans les modèles classiques (Black-Scholes, Heston, SABR, Bergomi standard) la volatilité stochastique est une **diffusion d'Itô** gouvernée par un mouvement Brownien $W$ : sa régularité est donc *höldérienne d'indice strictement inférieur à $1/2$*, mais \"proche\" de $1/2$. Depuis [GJR14], il est empiriquement admis que les séries historiques de volatilité (construites à partir de données intraday haute fréquence) ont en fait un indice de régularité **$H \\approx 0{,}05$--$0{,}15$**, bien inférieur à $1/2$. On dit que la volatilité est *rough*.\n",
    "\n",
    "La conséquence théorique est brutale : aucune diffusion d'Itô classique ne peut avoir de trajectoires aussi peu régulières. Il faut sortir du cadre markovien et introduire une *mémoire du passé* dans l'intégrande. Le modèle jouet est le **mouvement Brownien fractionnaire** $W^H$ d'indice de Hurst $H \\in (0,1)$, caractérisé par\n",
    "\n",
    "$$\\mathbb{E}[W^H_t W^H_s] = \\tfrac{1}{2}\\bigl(|t|^{2H} + |s|^{2H} - |t-s|^{2H}\\bigr).$$\n",
    "\n",
    "Pour $H = 1/2$ c'est un Brownien standard ; pour $H<1/2$ les accroissements sont **anti-corrélés** et les trajectoires plus rugueuses qu'un Brownien ; pour $H>1/2$, au contraire, elles sont *persistantes* (mémoire longue). Dans la littérature rough volatility, on prend $H<1/2$ -- la volatilité est donc à la fois **irrégulière** et **non-markovienne**.\n",
    "\n",
    "Une représentation très utile est la **représentation de Mandelbrot-Van Ness** (version Riemann-Liouville, tronquée à $t\\geq 0$) :\n",
    "\n",
    "$$W^H_t \\;=\\; \\frac{1}{\\Gamma(H+1/2)} \\int_0^t (t-s)^{H-1/2} \\, dW_s,$$\n",
    "\n",
    "avec $W$ Brownien standard. Reconnaissons un **noyau fractionnaire** $K_\\alpha(u) = u^{\\alpha-1}/\\Gamma(\\alpha)$ avec $\\alpha = H+1/2 \\in (0,1)$. Ce noyau est singulier en 0 (exposant $\\alpha-1<0$) et à décroissance polynomiale lente : c'est exactement ce qui introduit la *mémoire longue* et la *rugosité*."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "# Visualisation de K_alpha pour plusieurs H\n",
    "t = np.linspace(1e-3, 2.0, 400)\n",
    "fig, ax = plt.subplots(1, 2, figsize=(11, 3.5))\n",
    "for H in [0.05, 0.1, 0.25, 0.5, 0.75]:\n",
    "    ax[0].plot(t, fractional_kernel(t, H + 0.5), label=f'H={H}')\n",
    "ax[0].set(xlabel='t', ylabel=r'$K_\\alpha(t)$', title='noyau fractionnaire, $\\\\alpha=H+1/2$')\n",
    "ax[0].set_xscale('log'); ax[0].set_yscale('log'); ax[0].legend()\n",
    "# Terme de 'burst' phi(t) = t^{1-alpha}/Gamma(2-alpha)\n",
    "for H in [0.05, 0.1, 0.25, 0.5]:\n",
    "    ax[1].plot(t, phi_burst(t, H + 0.5), label=f'H={H}')\n",
    "ax[1].set(xlabel='t', ylabel=r'$\\phi(t) = t^{1-\\alpha}/\\Gamma(2-\\alpha)$',\n",
    "          title='terme de burst initial')\n",
    "ax[1].legend()\n",
    "fig.tight_layout(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "**Observations.** Le noyau fractionnaire explose en 0 et décroît comme une loi de puissance à l'infini. Plus $H$ est petit, plus la singularité à 0 est marquée (l'*incrément présent* a un impact énorme) et plus la décroissance à l'infini est lente : les modèles rough ont une **mémoire longue** (corrélations polynomiales plutôt qu'exponentielles).\n",
    "\n",
    "---\n",
    "## 2. Équations de Volterra stochastiques\n",
    "\n",
    "Une **équation de Volterra stochastique** est la généralisation des EDS d'Itô où le noyau $K$ est convolu avec le drift et la diffusion :\n",
    "\n",
    "$$X_t \\;=\\; \\xi_0 \\;+\\; \\int_0^t K(t-s)\\, b(s, X_s)\\, ds \\;+\\; \\int_0^t K(t-s)\\, \\sigma(s, X_s)\\, dW_s.$$\n",
    "\n",
    "Lorsque $K = \\mathbf{1}$ (constant), on retrouve une EDS usuelle. Dès que $K$ dépend de $t-s$ de façon non-triviale, le processus $X$ est **non-markovien** et n'est pas une semi-martingale en général. C'est le prix à payer pour capturer la mémoire.\n",
    "\n",
    "Dans [BCGP23], les auteurs considèrent une classe **path-dépendante** encore plus riche : le drift et la diffusion dépendent de la convolution $(\\bar K \\star X)_t$ où $\\bar K$ est un *co-noyau* satisfaisant $K\\star \\bar K = e^{-\\lambda\\cdot}$. C'est la clef du théorème principal qui suit.\n",
    "\n",
    "Nous nous concentrons sur le cas $\\lambda=0$, **kernel fractionnaire** $K(u) = u^{\\alpha-1}/\\Gamma(\\alpha)$ avec $\\alpha = H+1/2 \\in (1/2, 1)$. Alors $\\bar K(u) = u^{-\\alpha}/\\Gamma(1-\\alpha)$ et $(K\\star\\bar K)(t) = 1$ (identité de Riemann-Liouville).\n",
    "\n",
    "---\n",
    "## 3. Le pont éléphant / poisson-rouge (BCGP23)\n",
    "\n",
    "**Idée.** On souhaite simuler un processus $X$ rough (*mémoire d'éléphant*) à un coût comparable à celui d'une diffusion markovienne. Les auteurs montrent qu'il existe un processus auxiliaire markovien $\\eta$ (*mémoire de poisson-rouge*) dont la dynamique est une EDS **standard**, telle que :\n",
    "\n",
    "$$\\eta_t \\;=\\; \\xi_0\\, \\underbrace{\\frac{t^{1-\\alpha}}{\\Gamma(2-\\alpha)}}_{\\text{burst de mémoire}} \\;+\\; \\int_0^t b(s,\\eta_s)\\, ds \\;+\\; \\int_0^t \\sigma(s,\\eta_s)\\, dW_s,$$\n",
    "\n",
    "et le processus de Volterra d'origine se reconstruit simplement par **intégration fractionnaire** :\n",
    "\n",
    "$$Z_t \\;=\\; \\int_0^t \\frac{(t-s)^{\\alpha-1}}{\\Gamma(\\alpha)} \\bigl( b(s,\\eta_s)\\, ds + \\sigma(s,\\eta_s)\\, dW_s \\bigr).$$\n",
    "\n",
    "**Pourquoi c'est important.** Le schéma d'Euler standard appliqué à $\\eta$ converge fortement à la vitesse $h^{1/2}$ (comme n'importe quelle EDS d'Itô). On peut ensuite reconstruire $Z$ par combinaison linéaire des mêmes increments, et [BCGP23, Thm 5.3] montre que cette reconstruction garde la **même vitesse $h^{1/2}$**, *indépendamment* de $H$. C'est une amélioration considérable sur les schémas d'Euler directs pour le modèle rough Heston, dont le taux se dégrade à $H$ (typiquement $\\approx 0.1$).\n",
    "\n",
    "Le terme $\\xi_0\\, t^{1-\\alpha}/\\Gamma(2-\\alpha)$ s'appelle **burst de mémoire** : il encode (dès que $\\xi_0\\neq 0$) l'influence d'un \"passé avant 0\" fictif et produit un bond initial très rapide, visible sur les figures."
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "### 3.1 Illustration : trajectoires de $(\\eta, Z, V)$\n",
    "\n",
    "On choisit le modèle inspiré de [GJR20, BCGP23] :\n",
    "\n",
    "$$\\begin{cases}\n",
    "\\eta_t = \\eta_0\\, \\dfrac{t^{1-\\alpha}}{\\Gamma(2-\\alpha)} + \\displaystyle\\int_0^t \\nu(\\mu-\\eta_s)\\,ds + \\int_0^t \\theta\\,\\sigma(\\eta_s)\\,dW_s,\\\\[2mm]\n",
    "Z_t = \\displaystyle\\int_0^t \\dfrac{(t-s)^{\\alpha-1}}{\\Gamma(\\alpha)} \\bigl(\\nu(\\mu-\\eta_s)\\,ds + \\theta\\,\\sigma(\\eta_s)\\,dW_s\\bigr),\\\\[2mm]\n",
    "V_t = a(\\eta_t - b)^2 + c \\qquad \\text{(variance spot quadratique)},\\\\[2mm]\n",
    "dS_t = S_t \\sqrt{V_t}\\, dB_t, \\qquad d\\langle W,B\\rangle_t = \\rho\\, dt.\n",
    "\\end{cases}$$\n",
    "\n",
    "La dépendance quadratique $V = a(\\eta - b)^2 + c$ est l'ingrédient-clé de [GJR20] : elle introduit un *effet feedback fort* (Zumbach), crucial pour calibrer simultanément smiles SPX et VIX."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "T, N = 5.0, 8000\n",
    "cases = [\n",
    "    ('A : eta0 = 0',      dict(eta0=0.0,      mu=2.0, nu=1.2, theta=0.1)),\n",
    "    ('B : eta0 = mu/nu',  dict(eta0=2.0/1.2,  mu=2.0, nu=1.2, theta=0.1)),\n",
    "    ('C : rev. rapide',   dict(eta0=2.0/20.0, mu=2.0, nu=20.0, theta=0.01)),\n",
    "]\n",
    "fig, axes = plt.subplots(3, 3, figsize=(13, 8), sharex=True)\n",
    "for i, (title, kw) in enumerate(cases):\n",
    "    p = QuadRoughHestonParams(H=0.1, a=0.384, b=0.095, c=0.0025, **kw)\n",
    "    sim = simulate_elephant_goldfish(p, T, N, n_paths=1,\n",
    "                                     rng=np.random.default_rng(i+7))\n",
    "    for j, key in enumerate(['eta', 'Z', 'V']):\n",
    "        axes[i,j].plot(sim['t'], sim[key][0], lw=0.7)\n",
    "        axes[i,j].set_title(f'{title} -- {key}')\n",
    "fig.tight_layout(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "**Lecture des graphiques.**\n",
    "- **Ligne A (`eta0=0`).** Le burst est *nul* ($\\eta_0 \\phi(t) = 0$) et $\\eta$ part doucement. Comme dans le Lemme 5.4 de [BCGP23], $\\mathbb{E}[\\eta_t] \\to \\mu$ et $\\mathbb{E}[Z_t] \\to 0$ quand $t\\to\\infty$. L'éléphant $Z$ est significativement plus rugueux que $\\eta$ (on voit bien l'effet du noyau fractionnaire).\n",
    "- **Ligne B (`eta0 = mu/nu`).** Le burst est *visible* : $\\eta$ fait un bond quasi-instantané dès $t=0^+$. C'est l'interprétation \"mémoire d'avant zéro\" : on démarre comme si le processus avait déjà une histoire.\n",
    "- **Ligne C (forte vitesse de réversion).** Retour quasi-instantané vers la limite stationnaire. Ce régime est intéressant pour calibrer des maturités longues (c'est l'approche développée dans [P24] et sujet 6 du DUFE).\n",
    "\n",
    "**Remarque visuelle.** $\\eta$ est *globalement markovien* (une vraie EDS) : ses trajectoires sont $1/2^-$-höldériennes. En revanche $Z$ hérite de la rugosité du noyau et est $\\alpha-1/2 = H^-$-höldérien."
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "### 3.2 Effet du paramètre de Hurst $H$"
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "T, N = 1.0, 4000\n",
    "fig, axes = plt.subplots(1, 4, figsize=(14, 3))\n",
    "for i, H in enumerate([0.05, 0.1, 0.3, 0.49]):\n",
    "    p = QuadRoughHestonParams(H=H, mu=0.04, nu=1.0, theta=0.3, eta0=0.04)\n",
    "    sim = simulate_elephant_goldfish(p, T, N, n_paths=3,\n",
    "                                     rng=np.random.default_rng(i))\n",
    "    for j in range(3):\n",
    "        axes[i].plot(sim['t'], sim['Z'][j], lw=0.7)\n",
    "    axes[i].set(title=f'H={H}', xlabel='t')\n",
    "axes[0].set_ylabel('$Z_t$'); fig.tight_layout(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "On voit nettement la transition : pour $H$ petit, $Z$ a des oscillations \"en dents de scie\" incessantes ; pour $H\\to 1/2$, les trajectoires redeviennent \"diffusives\".\n",
    "\n",
    "---\n",
    "## 4. Schéma numérique et convergence forte\n",
    "\n",
    "### 4.1 Le schéma hybride\n",
    "\n",
    "On pose $t_k = kh$, $h = T/n$, et l'incrément\n",
    "\n",
    "$$\\Delta F_k \\;=\\; \\nu(\\mu - \\eta_{t_k}) h \\;+\\; \\theta\\,\\sigma(\\eta_{t_k})\\,\\Delta W_k.$$\n",
    "\n",
    "**Schéma pour $\\eta$** (Euler standard avec burst) :\n",
    "$$\\eta_{t_{k+1}} \\;=\\; \\eta_0\\, \\phi(t_{k+1}) + \\sum_{\\ell=0}^{k} \\Delta F_\\ell.$$\n",
    "\n",
    "**Schéma pour $Z$** -- deux variantes :\n",
    "\n",
    "(a) *Stepwise kernel* (éq. (5.5) de [BCGP23]) :\n",
    "$$Z_{t_{k+1}}^{(a)} \\;=\\; \\sum_{\\ell=0}^{k} \\frac{(t_{k+1}-t_\\ell)^{\\alpha-1}}{\\Gamma(\\alpha)}\\, \\Delta F_\\ell.$$\n",
    "\n",
    "(b) *Semi-integrated* (éq. (4.6)) : on intègre exactement le noyau sur chaque pas pour la partie drift :\n",
    "$$Z_{t_{k+1}}^{(b)} \\;=\\; \\sum_{\\ell=0}^{k} \\frac{(t_{k+1}-t_\\ell)^{\\alpha} - (t_{k+1}-t_{\\ell+1})^{\\alpha}}{\\Gamma(\\alpha+1)}\\, \\nu(\\mu-\\eta_{t_\\ell}) \\;+\\; \\sum_{\\ell=0}^{k} \\frac{(t_{k+1}-t_\\ell)^{\\alpha-1}}{\\Gamma(\\alpha)}\\, \\theta\\sigma(\\eta_{t_\\ell})\\,\\Delta W_\\ell.$$\n",
    "\n",
    "La complexité par trajectoire est $O(n^2)$ (calcul de tous les points), mais seulement $O(n)$ si l'on ne veut que $Z_T$.\n",
    "\n",
    "### 4.2 Étude numérique de la vitesse de convergence\n",
    "\n",
    "On utilise la méthode classique dite *\"de Cameron-Clark\"* : on simule un Brownien très fin $(N_{\\mathrm{ref}}\\gg 1)$, on sous-échantillonne les mêmes incréments pour obtenir des grilles $n \\ll N_{\\mathrm{ref}}$, et on mesure l'erreur forte\n",
    "$$\\mathrm{err}(h) \\;=\\; \\left\\|\\sup_{k}\\,|X^{h}_{t_k} - X^{h_{\\mathrm{ref}}}_{t_k}|\\right\\|_{L^2}.$$\n",
    "La théorie prédit $\\mathrm{err}(h) \\leq C\\, h^{1/2}$ (Thm 5.3 de [BCGP23]), *indépendamment de $H$*."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "T, N_ref = 1.0, 8192\n",
    "n_list = [32, 64, 128, 256, 512, 1024, 2048]\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 4))\n",
    "for H in [0.1, 0.25, 0.4]:\n",
    "    p = QuadRoughHestonParams(H=H, mu=0.04, nu=1.5, theta=0.3, eta0=0.04)\n",
    "    res = strong_error_study(p, T, n_list, N_ref, n_paths=800, p=2.0,\n",
    "                             rng=np.random.default_rng(int(1000*H)))\n",
    "    axes[0].loglog(res['h_list'], res['err_eta'], 'o-', label=f'H={H}')\n",
    "    axes[1].loglog(res['h_list'], res['err_Z'],   'o-', label=f'H={H}')\n",
    "hh = np.array([T/max(n_list), T/min(n_list)])\n",
    "axes[0].loglog(hh, 0.02*hh**0.5, 'k--', label='pente 1/2')\n",
    "axes[1].loglog(hh, 0.2*hh**0.5,  'k--', label='pente 1/2')\n",
    "axes[0].set(xlabel='h', ylabel='err L2 sur eta', title='Goldfish'); axes[0].legend()\n",
    "axes[1].set(xlabel='h', ylabel='err L2 sur Z',   title='Elephant'); axes[1].legend()\n",
    "fig.tight_layout(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "**Interprétation.** Les droites sont parallèles à la pente théorique $1/2$ et *ne se séparent pas* lorsque $H$ varie : c'est la confirmation empirique du résultat principal de [BCGP23]. Pour comparaison, un schéma d'Euler direct sur l'équation de Volterra (sans passer par le goldfish) aurait une pente $H \\approx 0.1$ : à pas $h = 10^{-3}$ cela signifie une erreur ~$10^{-0{,}3}\\approx 0.5$ vs. ~$10^{-1{,}5}\\approx 0.03$ pour notre schéma, soit un gain d'un facteur 15.\n",
    "\n",
    "---\n",
    "## 5. Quadratic rough Heston : trajectoires de l'actif et effet feedback\n",
    "\n",
    "Passons à l'actif sous-jacent $S$. On choisit $\\rho=-0{,}7$ (corrélation négative, effet leverage)."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "p = QuadRoughHestonParams(H=0.1, mu=0.04, nu=1.0, theta=0.3,\n",
    "                          a=0.384, b=0.095, c=0.0025, eta0=0.04)\n",
    "T, N = 1.0, 4000\n",
    "sim = simulate_asset(p, T, N, n_paths=5, S0=1.0, rho=-0.7,\n",
    "                     rng=np.random.default_rng(11))\n",
    "fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)\n",
    "for i in range(5):\n",
    "    axes[0].plot(sim['t'], sim['S'][i], lw=0.8)\n",
    "    axes[1].plot(sim['t'], np.sqrt(sim['V'][i]), lw=0.8)\n",
    "axes[0].set_ylabel('$S_t$'); axes[1].set(xlabel='t', ylabel='$\\\\sqrt{V_t}$')\n",
    "fig.tight_layout(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "**Effet feedback.** On observe nettement que les pics de volatilité suivent les baisses de $S$ (corrélation négative + noyau rough = mémoire longue des chocs récents). C'est un *fait stylisé* des marchés et une motivation empirique centrale du quadratic rough Heston.\n",
    "\n",
    "---\n",
    "## 6. Smile de volatilité implicite et term-structure du skew ATM\n",
    "\n",
    "Le modèle produit-il des smiles réalistes ? Pricing Monte-Carlo d'un call européen, puis inversion de Black-Scholes."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "S0 = 1.0\n",
    "strikes = np.array([0.85, 0.9, 0.95, 0.975, 1.0, 1.025, 1.05, 1.1, 1.15])\n",
    "fig, ax = plt.subplots(figsize=(8, 4.5))\n",
    "for T in [1/12, 3/12, 6/12, 1.0]:\n",
    "    N = max(200, int(400*T))\n",
    "    sim = simulate_asset(p, T, N, n_paths=30000, S0=S0, rho=-0.7,\n",
    "                         rng=np.random.default_rng(int(1000*T)))\n",
    "    sm = smile_from_paths(sim['S'], strikes, T, S0)\n",
    "    ax.plot(strikes/S0, 100*sm['iv'], 'o-', label=f'T={T:.3f}')\n",
    "ax.set(xlabel='K/S0', ylabel='IV (%)', title='Smile quadratic rough Heston')\n",
    "ax.legend(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "**Ce qu'on observe.**\n",
    "- Fort **skew négatif** (ATM plus haut que le call OTM) : signature du levier $\\rho<0$ combiné au feedback quadratique.\n",
    "- Le skew *explose* à courte maturité (cf. Bergomi [BE16]). C'est la signature du noyau rough : le skew ATM se comporte comme $T^{H-1/2}$ quand $T\\to 0$.\n",
    "- À maturité longue, le smile s'aplatit mais reste asymétrique.\n",
    "\n",
    "### 6.1 Structure par terme de l'ATM skew\n",
    "On définit le skew ATM via différence finie : $\\text{skew}(T) \\approx \\dfrac{\\mathrm{IV}(0{,}97) - \\mathrm{IV}(1{,}03)}{\\log(1{,}03) - \\log(0{,}97)}$. La théorie rough prévoit $|\\text{skew}(T)| \\sim T^{H-1/2}$."
  ]},
  {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [
    "Ts = np.array([1/52, 2/52, 1/12, 2/12, 3/12, 6/12, 9/12, 1.0, 1.5, 2.0])\n",
    "strikes_rel = np.array([0.97, 1.0, 1.03])\n",
    "skew = np.zeros_like(Ts)\n",
    "for i, T in enumerate(Ts):\n",
    "    N = max(200, int(400*T))\n",
    "    sim = simulate_asset(p, T, N, n_paths=40000, S0=S0, rho=-0.7,\n",
    "                         rng=np.random.default_rng(100+i))\n",
    "    sm = smile_from_paths(sim['S'], strikes_rel, T, S0)\n",
    "    skew[i] = (sm['iv'][0]-sm['iv'][2])/(np.log(strikes_rel[2])-np.log(strikes_rel[0]))\n",
    "mask = (Ts > 1/12) & (Ts < 1.5)\n",
    "coef = np.polyfit(np.log(Ts[mask]), np.log(np.abs(skew[mask])), 1)\n",
    "print(f'Pente empirique : {coef[0]:.3f}   /   theorie H-1/2 = {p.H-0.5:.3f}')\n",
    "fig, ax = plt.subplots(figsize=(7,4))\n",
    "ax.loglog(Ts, np.abs(skew), 'o', label='simulation')\n",
    "ax.loglog(Ts, np.exp(coef[1])*Ts**coef[0], 'r--', label=f'fit: pente {coef[0]:.2f}')\n",
    "ax.set(xlabel='T', ylabel='|ATM skew|'); ax.legend(); plt.show()"
  ]},
  {"cell_type": "markdown", "metadata": {}, "source": [
    "**Discussion.** La pente empirique sera typiquement de l'ordre de $-0{,}35$ à $-0{,}45$ (proche mais pas exactement égale à $H-1/2 = -0{,}4$). L'écart s'explique par :\n",
    "\n",
    "1. **Bruit Monte-Carlo** : la dispersion des estimateurs de IV à courte maturité est élevée.\n",
    "2. **Feedback quadratique non-linéaire** : strictement parlant le résultat théorique $T^{H-1/2}$ est établi pour le rough Bergomi linéaire ; le quadratic rough Heston modifie légèrement l'exposant.\n",
    "3. **Effet de bord** du burst initial : aux très courtes maturités, l'initial burst distord le smile.\n",
    "\n",
    "Ce résultat est néanmoins cohérent avec la littérature [BFG16, GE22] et valide qualitativement l'usage du modèle.\n",
    "\n",
    "---\n",
    "## 7. Discussion, limites et idées de recherche\n",
    "\n",
    "### Avantages du cadre BCGP / quadratic rough Heston\n",
    "1. **Vitesse $h^{1/2}$ indépendante de $H$** : un gain de temps de calcul majeur pour la calibration (10-100x plus rapide qu'un Euler direct).\n",
    "2. **Deux dynamiques couplées accessibles** : $\\eta$ markovien pour le pricing de vanille et $Z$ rough pour les faits stylisés historiques.\n",
    "3. **Burst de mémoire** modélise physiquement \"la mémoire d'avant 0\", sans introduire d'état initial fictif.\n",
    "4. **Effet feedback quadratique** : résout la conjecture de Guyon sur la calibration jointe SPX/VIX.\n",
    "\n",
    "### Limites observées\n",
    "1. Le **coût mémoire** du schéma reste $O(n^2)$ (matrice des poids fractionnaires). Pour $n \\geq 10^4$, il faut envisager une décomposition du noyau en somme d'exponentielles (*Multi-factor* de Abi Jaber, Callegaro-Fiorin-Grasselli, etc.).\n",
    "2. Le **choix du burst** $\\eta_0 \\phi(t)$ est un artefact : il dépend de la convention choisie pour \"l'information avant 0\". Pour des applications de calibration, il faudrait le voir comme un paramètre à ajuster (ou le remplacer par le régime stationnaire décrit dans le sujet 6 du DUFE).\n",
    "3. La **formule semi-explicite** de la fonction caractéristique du rough Heston (équation de Riccati fractionnaire) n'existe **pas** pour la version quadratique : pas de pricing par Fourier. Monte-Carlo est obligatoire.\n",
    "4. Le noyau fractionnaire **diverge en 0** : tout schéma perd en stabilité quand $H\\to 0$. Dans la pratique, $H<0{,}05$ est difficile à simuler.\n",
    "\n",
    "### Propositions de recherche\n",
    "\n",
    "**Idée 1 -- Décomposition multi-facteurs du noyau.**  \n",
    "Approximer $K_\\alpha(u) \\approx \\sum_{i=1}^{M} c_i e^{-\\lambda_i u}$ permettrait de ramener la simulation à un système d'EDS markoviennes en dimension $M$, avec complexité $O(n \\cdot M)$ (linéaire en $n$). Il existe déjà des approximations par quadrature de Gauss-Jacobi (Abi Jaber-El Euch 2019) ; il s'agirait de les combiner avec le pont BCGP pour une **convergence controllée**.\n",
    "\n",
    "**Idée 2 -- Quantification fonctionnelle.**  \n",
    "[BCGP23] mentionne explicitement comme piste future la quantification. On pourrait construire une grille quantifiée du *goldfish* $\\eta$ (processus markovien 1D) et utiliser [LMP22] pour le pricing d'options exotiques. Le gain : évaluer des prix avec quelques dizaines de points au lieu de dizaines de milliers de trajectoires MC.\n",
    "\n",
    "**Idée 3 -- Calibration hybride SPX/VIX.**  \n",
    "Tester empiriquement la capacité du modèle à calibrer simultanément les smiles SPX et VIX (problème historique [GJR20]). Le schéma rapide BCGP rend la *calibration* faisable (chaque évaluation du loss est $O(n)$).\n",
    "\n",
    "**Idée 4 -- Hedging par *pathwise* Malliavin.**  \n",
    "Les Greeks dans un modèle rough sont notoirement difficiles à calculer. Via le lien BCGP, on peut *réduire* leur calcul à des Greeks markoviens sur $\\eta$, puis transporter via la structure de l'intégration fractionnaire. C'est une piste peu explorée.\n",
    "\n",
    "**Idée 5 -- Extension path-dépendante complète.**  \n",
    "Ce projet a traité le cas où $b, \\sigma$ dépendent de $\\eta_s = (\\bar K \\star Z)_s$. Le cadre [BCGP23] permet plus général : $b(s, (\\bar K \\star Z)_s, Z_s)$, intégrant le *spot* et la *mémoire* conjointement. Cela se rapproche du cadre *path-dependent volatility* à la Guyon [G22], moyennant le changement de noyau.\n",
    "\n",
    "**Idée 6 -- Schéma variance-reduction pour le smile.**  \n",
    "Le Monte-Carlo des smiles est bruité aux extrémités de la moneyness. Utiliser comme variable de contrôle le prix Black-Scholes avec volatilité de forward (connue analytiquement pour rough Bergomi) devrait réduire la variance d'un facteur 10-100.\n",
    "\n",
    "---\n",
    "## Références\n",
    "\n",
    "- **[BCGP23]** O. Bonesini, G. Callegaro, M. Grasselli, G. Pagès. *From elephant to goldfish (and back)*. arXiv:2306.02708, 2023.\n",
    "- **[GJR20]** J. Gatheral, P. Jusselin, M. Rosenbaum. *The quadratic rough Heston model*. arXiv:2001.01789.\n",
    "- **[GJR14]** J. Gatheral, T. Jaisson, M. Rosenbaum. *Volatility is rough*. QF, 2018.\n",
    "- **[BFG16]** C. Bayer, P. Friz, J. Gatheral. *Pricing under rough volatility*. QF, 2016.\n",
    "- **[BE16]** L. Bergomi. *Stochastic Volatility Modeling*. CRC, 2016.\n",
    "- **[ER19]** O. El Euch, M. Rosenbaum. *The characteristic function of rough Heston models*. MF, 2019.\n",
    "- **[AJE19]** E. Abi Jaber, O. El Euch. *Multifactor approximation of rough volatility models*. SIAM JFM, 2019.\n",
    "- **[LMP22]** V. Lemaire, T. Montes, G. Pagès. *Stationary Heston model: Calibration and pricing of exotics using Product Recursive Quantization*. QF, 2022.\n",
    "- **[P24]** G. Pagès. *Volterra equations with affine drift: looking for stationarity*. arXiv:2401.15021.\n",
    "- **[G22]** J. Guyon. *The VIX future in Bergomi models*. QF, 2022.\n",
    "- **[GE22]** J. Guyon, M. El Amrani. *Does the term-structure of equity ATM skew follow a power law?*. SSRN 4174538."
  ]}
 ],
 "metadata": {
  "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"},
  "language_info": {"name": "python", "version": "3.13"}
 },
 "nbformat": 4,
 "nbformat_minor": 5
}