{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "### Προσομοίωση κίνησης Brown\n",
        "\n",
        "Θεωρούμε μια κίνηση Brown $(W_t)_{t \\ge 0}$, δηλαδή μια στοχαστική διαδικασία με τις ιδιότητες:\n",
        "\n",
        "- $W_0 = 0$ σχεδόν βέβαια,\n",
        "- $W_t - W_s \\sim \\mathcal{N}(0, t - s)$ για $0 \\le s < t$ (αύξηση με κανονική κατανομή),\n",
        "- ανεξάρτητες αυξήσεις,\n",
        "- σχεδόν σίγουρα συνεχείς τροχιές.\n",
        "\n",
        "Στο διακριτό πλέγμα χρόνου $0 = t_0 < t_1 < \\dots < t_N = T$ με βήμα\n",
        "$$\n",
        "\\Delta t = \\frac{T}{N},\n",
        "$$\n",
        "προσομοιώνουμε την κίνηση Brown με\n",
        "$$\n",
        "W_{t_{k+1}} = W_{t_k} + \\sqrt{\\Delta t}\\, Z_k,\n",
        "$$\n",
        "όπου $Z_k \\sim \\mathcal{N}(0,1)$ ανεξάρτητες τ.μ.\n",
        "\n",
        "Στο επόμενο κελί γράφουμε μια συνάρτηση Python που:\n",
        "\n",
        "- δημιουργεί το πλέγμα χρόνου,\n",
        "- προσομοιώνει $n_{\\text{paths}}$ τροχιές,\n",
        "- τις απεικονίζει γραφικά.\n"
      ],
      "metadata": {
        "id": "Y6QZA5MtSWtb"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "def simulate_brownian_motion(T=1.0, N=1000, n_paths=5, seed=None):\n",
        "    \"\"\"\n",
        "    Προσομοίωση κίνησης Brown σε [0, T] με βήμα dt = T/N.\n",
        "    Επιστρέφει:\n",
        "      t: πίνακας χρόνου σχήματος (N+1,)\n",
        "      W: πίνακας τροχιών σχήματος (n_paths, N+1)\n",
        "    \"\"\"\n",
        "    if seed is not None:\n",
        "        np.random.seed(seed)\n",
        "\n",
        "    dt = T / N\n",
        "    t = np.linspace(0.0, T, N + 1)\n",
        "\n",
        "    # Κανονικές αυξήσεις ~ N(0, dt)\n",
        "    dW = np.sqrt(dt) * np.random.randn(n_paths, N)\n",
        "\n",
        "    # Συσσώρευση για να πάρουμε W_t\n",
        "    W = np.zeros((n_paths, N + 1))\n",
        "    W[:, 1:] = np.cumsum(dW, axis=1)\n",
        "\n",
        "    return t, W\n",
        "\n",
        "# Παράδειγμα χρήσης\n",
        "T = 1.0\n",
        "N = 1000\n",
        "n_paths = 5\n",
        "\n",
        "t, W = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=123)\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "for i in range(n_paths):\n",
        "    plt.plot(t, W[i], alpha=0.8)\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"W_t\")\n",
        "plt.title(\"Προσομοιωμένες τροχιές κίνησης Brown\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 353
        },
        "id": "v3XB8aoOSYgS",
        "outputId": "0745e68b-2ab4-47fc-b08a-c3ad214291b2"
      },
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 110840 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Έλεγχος της κατανομής του $W_t$\n",
        "\n",
        "Για κάθε $t \\ge 0$ ισχύει θεωρητικά:\n",
        "$$\n",
        "W_t \\sim \\mathcal{N}(0, t).\n",
        "$$\n",
        "\n",
        "Θα επιλέξουμε ένα συγκεκριμένο χρόνο, π.χ. $t = T$, και θα:\n",
        "\n",
        "1. Προσομοιώσουμε πολλές ανεξάρτητες τροχιές κίνησης Brown.\n",
        "2. Θα κρατήσουμε τις τιμές $W_T$.\n",
        "3. Θα σχεδιάσουμε το ιστόγραμμα των $W_T$.\n",
        "4. Θα το συγκρίνουμε με την πυκνότητα της κανονικής $\\mathcal{N}(0, T)$.\n",
        "\n",
        "Αυτό οπτικοποιεί τον ορισμό της κίνησης Brown.\n"
      ],
      "metadata": {
        "id": "1XEpX6u1Wd48"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from math import sqrt, pi, exp\n",
        "\n",
        "# Παράμετροι\n",
        "T = 1.0\n",
        "N = 1000\n",
        "n_paths = 10000\n",
        "\n",
        "t, W = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=42)\n",
        "W_T = W[:, -1]  # τιμές στο χρόνο T\n",
        "\n",
        "# Θεωρητική πυκνότητα N(0, T)\n",
        "def normal_pdf(x, mean=0.0, var=1.0):\n",
        "    return (1.0 / np.sqrt(2 * np.pi * var)) * np.exp(-(x - mean)**2 / (2 * var))\n",
        "\n",
        "x_vals = np.linspace(W_T.min(), W_T.max(), 200)\n",
        "pdf_vals = normal_pdf(x_vals, mean=0.0, var=T)\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "# Ιστόγραμμα κανονικοποιημένο\n",
        "plt.hist(W_T, bins=50, density=True, alpha=0.6, label=\"Εμπειρική πυκνότητα\")\n",
        "plt.plot(x_vals, pdf_vals, linewidth=2, label=\"Θεωρητική N(0, T)\")\n",
        "plt.xlabel(\"w\")\n",
        "plt.ylabel(\"πυκνότητα\")\n",
        "plt.title(f\"Κατανομή του W_T για T = {T}\")\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 357
        },
        "id": "FE-ns36XWhB4",
        "outputId": "62b5a4b4-372b-4269-aaec-cf7115fde665"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 41540 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Brownian Bridge από κίνηση Brown\n",
        "\n",
        "Στο βιβλίο ορίζεται η **Brownian Bridge** στο $[0,T]$ ως στοχαστική διαδικασία $(X_t)_{0 \\le t \\le T}$ με\n",
        "\n",
        "- $X_0 = X_T = 0$,\n",
        "- κατάλληλες κανονικές κατανομές και συνδιακυμάνσεις.\n",
        "\n",
        "Ένα κλασικό παράδειγμα είναι:\n",
        "$$\n",
        "X_t = W_t - \\frac{t}{T} W_T, \\qquad 0 \\le t \\le T,\n",
        "$$\n",
        "όπου $(W_t)$ είναι κίνηση Brown. Τότε $X_0 = X_T = 0$ και $X_t$ είναι Brownian Bridge.\n",
        "\n",
        "Στο επόμενο κελί:\n",
        "- Προσομοιώνουμε μια κίνηση Brown $W$,\n",
        "- Ορίζουμε την Brownian Bridge $X_t$,\n",
        "- Τη σχεδιάζουμε γραφικά.\n"
      ],
      "metadata": {
        "id": "v78q__3uWlDE"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def simulate_brownian_bridge(T=1.0, N=1000, n_paths=5, seed=None):\n",
        "    t, W = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=seed)\n",
        "    # X_t = W_t - (t/T)*W_T\n",
        "    W_T = W[:, -1][:, None]      # shape (n_paths,1)\n",
        "    t_over_T = (t / T)[None, :]  # shape (1, N+1)\n",
        "    X = W - t_over_T * W_T\n",
        "    return t, X, W\n",
        "\n",
        "T = 1.0\n",
        "N = 1000\n",
        "n_paths = 5\n",
        "\n",
        "t, X, W = simulate_brownian_bridge(T=T, N=N, n_paths=n_paths, seed=123)\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "for i in range(n_paths):\n",
        "    plt.plot(t, X[i], alpha=0.8)\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"X_t\")\n",
        "plt.title(\"Προσομοιωμένες τροχιές Brownian Bridge\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 348
        },
        "id": "9dXxCA8GWoUk",
        "outputId": "32896380-2057-4021-a723-303900d69261"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 136572 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Αριθμητική προσέγγιση στοχαστικού ολοκληρώματος με ντετερμινιστικό ολοκληρωτέο\n",
        "\n",
        "Στην ενότητα 10.8.1 ορίζεται το στοχαστικό ολοκλήρωμα τύπου Itô\n",
        "$$\n",
        "\\int_0^T f(s) \\, dW_s\n",
        "$$\n",
        "αρχικά για απλές (προβλέψιμες) διαδικασίες και στη συνέχεια επεκτείνεται σε $L^2$.\n",
        "\n",
        "Αν $f$ είναι ντετερμινιστική και τετραγωνικά ολοκληρώσιμη στο $[0,T]$, τότε\n",
        "$$\n",
        "\\int_0^T f(s) \\, dW_s \\sim \\mathcal{N}\\!\\left(0, \\int_0^T f(s)^2\\,ds\\right),\n",
        "$$\n",
        "και\n",
        "$$\n",
        "\\mathbb{E}\\left[\\left(\\int_0^T f(s)\\, dW_s\\right)^2\\right]\n",
        "= \\int_0^T f(s)^2\\, ds\n",
        "$$\n",
        "(ισομετρία του Itô).\n",
        "\n",
        "Σε διακριτό χρόνο, προσεγγίζουμε:\n",
        "$$\n",
        "\\int_0^T f(s)\\, dW_s \\approx \\sum_{k=0}^{N-1} f(t_k)\\, \\Delta W_k,\n",
        "$$\n",
        "με $t_k = k \\Delta t$ και $\\Delta W_k = W_{t_{k+1}} - W_{t_k}$.\n",
        "\n",
        "Θα πάρουμε $f(s) = \\sin(s)$, θα προσομοιώσουμε πολλές φορές το ολοκλήρωμα, και θα συγκρίνουμε την εμπειρική διασπορά με την θεωρητική\n",
        "$$\n",
        "\\int_0^T \\sin^2(s)\\, ds = \\frac{T}{2} - \\frac{\\sin(2T)}{4}.\n",
        "$$\n"
      ],
      "metadata": {
        "id": "p_JOyLKXWs4Q"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def ito_integral_deterministic(f, T=1.0, N=1000, n_paths=5000, seed=None):\n",
        "    \"\"\"\n",
        "    Προσομοίωση του ∫_0^T f(s) dW_s με απλή προσεγγιστική άθροιση.\n",
        "    Επιστρέφει έναν πίνακα (n_paths,) με τις προσομοιωμένες τιμές.\n",
        "    \"\"\"\n",
        "    if seed is not None:\n",
        "        np.random.seed(seed)\n",
        "\n",
        "    dt = T / N\n",
        "    t = np.linspace(0.0, T, N + 1)\n",
        "\n",
        "    # Προσομοίωση αυξήσεων Brown\n",
        "    dW = np.sqrt(dt) * np.random.randn(n_paths, N)\n",
        "\n",
        "    # f(t_k) για k = 0,...,N-1\n",
        "    f_vals = f(t[:-1])  # shape (N,)\n",
        "\n",
        "    # Στοχαστικό άθροισμα\n",
        "    integrals = np.sum(f_vals[None, :] * dW, axis=1)\n",
        "    return integrals\n",
        "\n",
        "# Παράδειγμα με f(s) = sin(s)\n",
        "T = 1.0\n",
        "N = 1000\n",
        "n_paths = 5000\n",
        "\n",
        "f = np.sin\n",
        "integrals = ito_integral_deterministic(f, T=T, N=N, n_paths=n_paths, seed=123)\n",
        "\n",
        "emp_mean = np.mean(integrals)\n",
        "emp_var = np.var(integrals)\n",
        "\n",
        "# Θεωρητική διασπορά ∫_0^T sin^2(s) ds\n",
        "theo_var = T/2 - np.sin(2*T)/4\n",
        "\n",
        "print(f\"Εμπειρική μέση τιμή ≈ {emp_mean:.4f} (θεωρητικά 0)\")\n",
        "print(f\"Εμπειρική διασπορά ≈ {emp_var:.4f}\")\n",
        "print(f\"Θεωρητική διασπορά   = {theo_var:.4f}\")\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.hist(integrals, bins=50, density=True, alpha=0.6)\n",
        "plt.title(r\"Ιστόγραμμα $\\int_0^T \\sin(s)\\, dW_s$\")\n",
        "plt.xlabel(\"τιμή ολοκληρώματος\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 432
        },
        "id": "1njgQAWzWvr8",
        "outputId": "72158aa5-ace5-4fdd-e8fd-d67720e47317"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Εμπειρική μέση τιμή ≈ 0.0080 (θεωρητικά 0)\n",
            "Εμπειρική διασπορά ≈ 0.2724\n",
            "Θεωρητική διασπορά   = 0.2727\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 21884 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Παράδειγμα Itô: $\\int_0^T W_s\\, dW_s = \\frac{1}{2}(W_T^2 - T)$\n",
        "\n",
        "Με τη φόρμουλα του Itô για $f(x) = x^2$ παίρνουμε\n",
        "$$\n",
        "d(W_t^2) = 2 W_t\\, dW_t + dt,\n",
        "$$\n",
        "άρα\n",
        "$$\n",
        "\\int_0^T W_s\\, dW_s\n",
        "= \\frac{1}{2}\\left(W_T^2 - T\\right).\n",
        "$$\n",
        "\n",
        "Θα προσομοιώσουμε:\n",
        "\n",
        "1. Τροχιές κίνησης Brown $(W_t)$.\n",
        "2. Το στοχαστικό άθροισμα\n",
        "   $$\n",
        "   \\sum_{k=0}^{N-1} W_{t_k}\\, \\Delta W_k\n",
        "   $$\n",
        "   ως προσέγγιση του $\\int_0^T W_s\\, dW_s$.\n",
        "3. Τη διαφορά με το θεωρητικό $\\frac{1}{2}(W_T^2 - T)$.\n",
        "\n",
        "Θα ελέγξουμε ότι η διαφορά είναι μικρή (κατά μέσο όρο).\n"
      ],
      "metadata": {
        "id": "x_O2la2dW0kJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def ito_integral_W_dW(T=1.0, N=2000, n_paths=2000, seed=None):\n",
        "    if seed is not None:\n",
        "        np.random.seed(seed)\n",
        "\n",
        "    dt = T / N\n",
        "    t = np.linspace(0.0, T, N + 1)\n",
        "\n",
        "    # Προσομοίωση Brown\n",
        "    dW = np.sqrt(dt) * np.random.randn(n_paths, N)\n",
        "    W = np.zeros((n_paths, N + 1))\n",
        "    W[:, 1:] = np.cumsum(dW, axis=1)\n",
        "\n",
        "    # Στοχαστικό άθροισμα ∑ W_{t_k} dW_k\n",
        "    integral_approx = np.sum(W[:, :-1] * dW, axis=1)\n",
        "\n",
        "    # Θεωρητικό: 0.5 * (W_T^2 - T)\n",
        "    W_T = W[:, -1]\n",
        "    integral_theo = 0.5 * (W_T**2 - T)\n",
        "\n",
        "    return integral_approx, integral_theo\n",
        "\n",
        "T = 1.0\n",
        "N = 2000\n",
        "n_paths = 5000\n",
        "\n",
        "approx_vals, theo_vals = ito_integral_W_dW(T=T, N=N, n_paths=n_paths, seed=42)\n",
        "diff = approx_vals - theo_vals\n",
        "\n",
        "print(f\"Μέση διαφορά (προσομοίωση - θεωρία): {np.mean(diff):.4e}\")\n",
        "print(f\"Τυπική απόκλιση διαφοράς: {np.std(diff):.4e}\")\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.hist(diff, bins=50, density=True, alpha=0.6)\n",
        "plt.title(r\"Διαφορά $\\sum W_{t_k}\\Delta W_k - \\frac{1}{2}(W_T^2 - T)$\")\n",
        "plt.xlabel(\"διαφορά\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 410
        },
        "id": "18oqbXMWW3nK",
        "outputId": "6462b8c5-dded-4c59-ad45-0cb07a123028"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Μέση διαφορά (προσομοίωση - θεωρία): -1.3584e-05\n",
            "Τυπική απόκλιση διαφοράς: 1.5680e-02\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 18908 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Φόρμουλα Itô για μονοδιάστατη διαδικασία Itô\n",
        "\n",
        "Αν $(X_t)$ είναι διαδικασία Itô\n",
        "$$\n",
        "dX_t = a(t)\\, dt + b(t)\\, dW_t,\n",
        "$$\n",
        "και $f(t,x)$ έχει επαρκώς ομαλές μερικές παραγώγους, τότε (φόρμουλα Itô)\n",
        "$$\n",
        "df(t,X_t) =\n",
        "\\left(f_t(t,X_t) + a(t) f_x(t,X_t) + \\tfrac{1}{2} b(t)^2 f_{xx}(t,X_t)\\right) dt\n",
        "+ b(t) f_x(t,X_t)\\, dW_t.\n",
        "$$\n",
        "\n",
        "Για $X_t = W_t$ (καθαρή κίνηση Brown) και $f(x) = x^2$, παίρνουμε:\n",
        "$$\n",
        "df(W_t) = d(W_t^2) = 2 W_t\\, dW_t + dt,\n",
        "$$\n",
        "άρα\n",
        "$$\n",
        "W_t^2 - t\n",
        "$$\n",
        "είναι martingale, και ειδικότερα\n",
        "$$\n",
        "\\mathbb{E}[W_t^2] = t.\n",
        "$$\n",
        "\n",
        "Στο επόμενο κελί:\n",
        "- Προσομοιώνουμε πολλές τροχιές $(W_t)$,\n",
        "- Υπολογίζουμε εμπειρικά το $\\mathbb{E}[W_t^2]$ σε διάφορους χρόνους $t$,\n",
        "- Το συγκρίνουμε με την ευθεία $t$.\n"
      ],
      "metadata": {
        "id": "4eNYhr49W-6y"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Παράμετροι\n",
        "T = 1.0\n",
        "N = 1000\n",
        "n_paths = 5000\n",
        "\n",
        "t, W = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=1)\n",
        "\n",
        "# Επιλέγουμε μερικά σημεία χρόνου\n",
        "times_idx = np.linspace(0, N, 11, dtype=int)  # 11 σημεία\n",
        "times = t[times_idx]\n",
        "\n",
        "E_W2 = np.mean(W[:, times_idx]**2, axis=0)\n",
        "\n",
        "for ti, ew2 in zip(times, E_W2):\n",
        "    print(f\"t = {ti:.2f}, E[W_t^2] ≈ {ew2:.4f}\")\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.plot(times, E_W2, \"o-\", label=r\"Εμπειρικό $E[W_t^2]$\")\n",
        "plt.plot(times, times, \"--\", label=\"t (θεωρητικό)\")\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(r\"$E[W_t^2]$\")\n",
        "plt.title(r\"Έλεγχος ότι $E[W_t^2] = t$\")\n",
        "plt.grid(True)\n",
        "plt.legend()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 548
        },
        "id": "oCUwDPoqXBkI",
        "outputId": "81db70c1-d9f8-4507-c614-28fa1cf75a76"
      },
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "t = 0.00, E[W_t^2] ≈ 0.0000\n",
            "t = 0.10, E[W_t^2] ≈ 0.1000\n",
            "t = 0.20, E[W_t^2] ≈ 0.2026\n",
            "t = 0.30, E[W_t^2] ≈ 0.2972\n",
            "t = 0.40, E[W_t^2] ≈ 0.3966\n",
            "t = 0.50, E[W_t^2] ≈ 0.4933\n",
            "t = 0.60, E[W_t^2] ≈ 0.5988\n",
            "t = 0.70, E[W_t^2] ≈ 0.6976\n",
            "t = 0.80, E[W_t^2] ≈ 0.7982\n",
            "t = 0.90, E[W_t^2] ≈ 0.9079\n",
            "t = 1.00, E[W_t^2] ≈ 1.0147\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 40384 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Εκθετική martingale: $X_t = \\exp(W_t - \\tfrac{1}{2} t)$\n",
        "\n",
        "Στο Παράδειγμα 551 ορίζεται\n",
        "$$\n",
        "X_t = \\exp\\!\\left(W_t - \\frac{t}{2}\\right).\n",
        "$$\n",
        "Η φόρμουλα Itô δείχνει ότι $X_t$ είναι martingale και ειδικά\n",
        "$$\n",
        "\\mathbb{E}[X_t] = \\mathbb{E}\\left[e^{W_t - t/2}\\right] = 1.\n",
        "$$\n",
        "\n",
        "Θα το ελέγξουμε αριθμητικά:\n",
        "\n",
        "- Προσομοιώνουμε πολλές τροχιές $(W_t)$,\n",
        "- Υπολογίζουμε $X_t$ σε διάφορους χρόνους,\n",
        "- Βλέπουμε ότι η εμπειρική μέση τιμή είναι περίπου $1$.\n"
      ],
      "metadata": {
        "id": "Jg9AfIrKXF-T"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "T = 2.0\n",
        "N = 2000\n",
        "n_paths = 5000\n",
        "\n",
        "t, W = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=7)\n",
        "X = np.exp(W - 0.5 * t[None, :])  # X_t = exp(W_t - t/2)\n",
        "\n",
        "times_idx = np.linspace(0, N, 9, dtype=int)\n",
        "times = t[times_idx]\n",
        "E_X = np.mean(X[:, times_idx], axis=0)\n",
        "\n",
        "for ti, ex in zip(times, E_X):\n",
        "    print(f\"t = {ti:.2f}, E[X_t] ≈ {ex:.4f}\")\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.plot(times, E_X, \"o-\", label=r\"Εμπειρικό $E[X_t]$\")\n",
        "plt.axhline(1.0, color=\"k\", linestyle=\"--\", label=\"1 (θεωρητικό)\")\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(r\"$E[X_t]$\")\n",
        "plt.title(r\"Martingale $X_t = e^{W_t - t/2}$\")\n",
        "plt.grid(True)\n",
        "plt.legend()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 511
        },
        "id": "Ne7sdinuXI8h",
        "outputId": "1f2de582-f418-44ec-9565-c945136c19a6"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "t = 0.00, E[X_t] ≈ 1.0000\n",
            "t = 0.25, E[X_t] ≈ 0.9933\n",
            "t = 0.50, E[X_t] ≈ 0.9831\n",
            "t = 0.75, E[X_t] ≈ 0.9755\n",
            "t = 1.00, E[X_t] ≈ 0.9748\n",
            "t = 1.25, E[X_t] ≈ 0.9757\n",
            "t = 1.50, E[X_t] ≈ 0.9690\n",
            "t = 1.75, E[X_t] ≈ 0.9618\n",
            "t = 2.00, E[X_t] ≈ 0.9508\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 42140 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Γενική προσομοίωση SDE με σχήμα Euler–Maruyama\n",
        "\n",
        "Μια στοχαστική διαφορική εξίσωση Itô έχει τη μορφή\n",
        "$$\n",
        "dX_t = a(t, X_t)\\, dt + b(t, X_t)\\, dW_t, \\qquad X_0 = x_0.\n",
        "$$\n",
        "\n",
        "Σε διακριτό πλέγμα $0 = t_0 < t_1 < \\dots < t_N = T$ με $\\Delta t = T/N$ χρησιμοποιούμε το σχήμα Euler–Maruyama:\n",
        "$$\n",
        "X_{k+1} = X_k + a(t_k, X_k)\\, \\Delta t + b(t_k, X_k)\\, \\Delta W_k,\n",
        "$$\n",
        "όπου $\\Delta W_k \\sim \\mathcal{N}(0, \\Delta t)$.\n",
        "\n",
        "Στο επόμενο κελί γράφουμε μια γενική συνάρτηση Python που προσομοιώνει λύσεις SDE με αυτό το σχήμα.\n"
      ],
      "metadata": {
        "id": "MUJOIl6zXM0q"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def euler_maruyama(a, b, x0, T=1.0, N=1000, n_paths=1, seed=None):\n",
        "    \"\"\"\n",
        "    Euler–Maruyama για SDE:\n",
        "      dX_t = a(t, X_t) dt + b(t, X_t) dW_t,  X_0 = x0\n",
        "    Επιστρέφει:\n",
        "      t: (N+1,)\n",
        "      X: (n_paths, N+1)\n",
        "    \"\"\"\n",
        "    if seed is not None:\n",
        "        np.random.seed(seed)\n",
        "\n",
        "    dt = T / N\n",
        "    t = np.linspace(0.0, T, N + 1)\n",
        "\n",
        "    X = np.zeros((n_paths, N + 1))\n",
        "    X[:, 0] = x0\n",
        "\n",
        "    dW = np.sqrt(dt) * np.random.randn(n_paths, N)\n",
        "\n",
        "    for k in range(N):\n",
        "        tk = t[k]\n",
        "        Xk = X[:, k]\n",
        "        X[:, k+1] = Xk + a(tk, Xk) * dt + b(tk, Xk) * dW[:, k]\n",
        "\n",
        "    return t, X\n",
        "\n",
        "# Μικρό test με a=0, b=1 (πρέπει να έχουμε κίνηση Brown)\n",
        "a = lambda t, x: 0.0\n",
        "b = lambda t, x: 1.0\n",
        "\n",
        "t_sde, X_sde = euler_maruyama(a, b, x0=0.0, T=1.0, N=1000, n_paths=5, seed=0)\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "for i in range(5):\n",
        "    plt.plot(t_sde, X_sde[i], alpha=0.8)\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"X_t\")\n",
        "plt.title(\"SDE με a=0, b=1 (κίνηση Brown)\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 353
        },
        "id": "aNJ5pR2-XPkj",
        "outputId": "2c16c1e2-f0bb-404e-b055-ba7904f76ab2"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 116680 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Διαδικασία Ornstein–Uhlenbeck\n",
        "\n",
        "Στο Παράδειγμα 552 (10.8.2) εμφανίζεται η διαδικασία Ornstein–Uhlenbeck. Η αντίστοιχη SDE είναι\n",
        "$$\n",
        "dY_t = -a Y_t\\, dt + b\\, dW_t, \\qquad a > 0, \\; b \\in \\mathbb{R}.\n",
        "$$\n",
        "\n",
        "Η διαδικασία αυτή είναι μέση-επιστρεφόμενη (mean-reverting) προς το $0$.\n",
        "\n",
        "Θα χρησιμοποιήσουμε τον αλγόριθμο Euler–Maruyama για να προσομοιώσουμε μερικές τροχιές.\n"
      ],
      "metadata": {
        "id": "2g4BejFtXT6Z"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "a_param = 2.0\n",
        "b_param = 1.0\n",
        "\n",
        "def a_OU(t, x):\n",
        "    return -a_param * x\n",
        "\n",
        "def b_OU(t, x):\n",
        "    return b_param * np.ones_like(x)\n",
        "\n",
        "T = 2.0\n",
        "N = 2000\n",
        "n_paths = 5\n",
        "x0 = 0.0\n",
        "\n",
        "t_OU, Y_OU = euler_maruyama(a_OU, b_OU, x0=x0, T=T, N=N, n_paths=n_paths, seed=1)\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "for i in range(n_paths):\n",
        "    plt.plot(t_OU, Y_OU[i], alpha=0.8)\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"Y_t\")\n",
        "plt.title(\"Διαδικασία Ornstein–Uhlenbeck (Euler–Maruyama)\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 353
        },
        "id": "jKXl7iYFXWhm",
        "outputId": "e31b837a-b03d-44da-9bc6-649e7ab802cf"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 139500 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Γεωμετρική κίνηση Brown (GBM)\n",
        "\n",
        "Η γεωμετρική κίνηση Brown $S_t$ ικανοποιεί την SDE\n",
        "$$\n",
        "dS_t = \\mu S_t\\, dt + \\sigma S_t\\, dW_t, \\qquad S_0 > 0.\n",
        "$$\n",
        "\n",
        "Η κλειστή μορφή της λύσης είναι\n",
        "$$\n",
        "S_t = S_0 \\exp\\!\\left( \\left(\\mu - \\tfrac{1}{2}\\sigma^2\\right) t + \\sigma W_t \\right).\n",
        "$$\n",
        "\n",
        "Στο επόμενο κελί:\n",
        "\n",
        "1. Προσομοιώνουμε $S_t$ με Euler–Maruyama.\n",
        "2. Προσομοιώνουμε και με την κλειστή μορφή χρησιμοποιώντας κίνηση Brown.\n",
        "3. Συγκρίνουμε τις τροχιές.\n"
      ],
      "metadata": {
        "id": "0Az00BtFXdxX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mu = 0.1\n",
        "sigma = 0.3\n",
        "\n",
        "def a_gbm(t, S):\n",
        "    return mu * S\n",
        "\n",
        "def b_gbm(t, S):\n",
        "    return sigma * S\n",
        "\n",
        "S0 = 1.0\n",
        "T = 2.0\n",
        "N = 2000\n",
        "n_paths = 3\n",
        "\n",
        "t_gbm, S_euler = euler_maruyama(a_gbm, b_gbm, x0=S0, T=T, N=N, n_paths=n_paths, seed=2)\n",
        "\n",
        "# Κλειστή μορφή με την ίδια Brownian κίνηση για μια τροχιά\n",
        "t_BM, W_BM = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=2)\n",
        "S_closed = S0 * np.exp((mu - 0.5 * sigma**2) * t_BM[None, :] + sigma * W_BM)\n",
        "\n",
        "plt.figure(figsize=(10, 4))\n",
        "\n",
        "for i in range(n_paths):\n",
        "    plt.plot(t_gbm, S_euler[i], alpha=0.8, label=f\"Euler path {i+1}\" if i == 0 else None)\n",
        "    plt.plot(t_BM, S_closed[i], \"--\", alpha=0.8, label=f\"Closed form {i+1}\" if i == 0 else None)\n",
        "\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"S_t\")\n",
        "plt.title(\"Γεωμετρική κίνηση Brown: Euler–Maruyama vs κλειστή λύση\")\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 296
        },
        "id": "-vm-CeAEXgS7",
        "outputId": "99473ae1-ecc9-412f-bd7a-2edea7091f86"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 140732 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Παράδειγμα στοχαστικού Fubini\n",
        "\n",
        "Στην ενότητα 10.8.4 δίνεται το παράδειγμα\n",
        "$$\n",
        "Y_t = \\int_0^t W_r \\, dr.\n",
        "$$\n",
        "\n",
        "Χρησιμοποιώντας το στοχαστικό θεώρημα Fubini, μπορούμε να γράψουμε\n",
        "$$\n",
        "Y_t = \\int_0^t (t - s)\\, dW_s.\n",
        "$$\n",
        "\n",
        "Θα ελέγξουμε αριθμητικά ότι οι δύο εκφράσεις δίνουν (σχεδόν) την ίδια ποσότητα:\n",
        "\n",
        "1. Προσομοιώνουμε μια τροχιά $W_s$.\n",
        "2. Υπολογίζουμε προσεγγιστικά\n",
        "   $$\n",
        "   Y_t^{(1)} \\approx \\sum W_{r_k}\\, \\Delta r_k\n",
        "   $$\n",
        "   ως κλασικό ολοκλήρωμα Riemann.\n",
        "3. Υπολογίζουμε προσεγγιστικά\n",
        "   $$\n",
        "   Y_t^{(2)} \\approx \\sum (t - s_k)\\, \\Delta W_k\n",
        "   $$\n",
        "   ως στοχαστικό ολοκλήρωμα Itô.\n",
        "4. Συγκρίνουμε $Y_t^{(1)}$ και $Y_t^{(2)}$.\n"
      ],
      "metadata": {
        "id": "hMuDXi3cXkEb"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "T = 1.0\n",
        "N = 4000\n",
        "dt = T / N\n",
        "t, W = simulate_brownian_motion(T=T, N=N, n_paths=1, seed=123)\n",
        "W = W[0]  # μια τροχιά\n",
        "\n",
        "# 1) Προσέγγιση Y_t = ∫_0^t W_r dr με Riemann\n",
        "Y1 = np.cumsum(W[:-1] * dt)  # μέθοδος \"left Riemann\"\n",
        "Y1 = np.concatenate([[0.0], Y1])  # Y_0 = 0\n",
        "\n",
        "# 2) Προσέγγιση Y_t = ∫_0^t (t-s) dW_s\n",
        "dW = np.diff(W)  # ΔW_k\n",
        "# (t - s_k) όπου s_k = t_k\n",
        "kernel = (T - t[:-1])  # αν θέλουμε συγκεκριμένο t = T\n",
        "Y2_T = np.sum(kernel * dW)  # τιμή στο χρόνο T\n",
        "\n",
        "print(f\"Προσέγγιση ∫_0^T W_r dr με Riemann: {Y1[-1]:.4f}\")\n",
        "print(f\"Προσέγγιση ∫_0^T (T-s) dW_s με Itô: {Y2_T:.4f}\")\n",
        "\n",
        "# Γραφική παράσταση της τροχιάς (για οπτικοποίηση)\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.plot(t, W)\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"W_t\")\n",
        "plt.title(\"Τροχιά κίνησης Brown για το παράδειγμα Fubini\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 387
        },
        "id": "R4F6q6psXnnH",
        "outputId": "ead8b398-103b-4768-be34-e312de6ee3b0"
      },
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Προσέγγιση ∫_0^T W_r dr με Riemann: -0.1621\n",
            "Προσέγγιση ∫_0^T (T-s) dW_s με Itô: -0.1619\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 54876 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Απλό αριθμητικό παράδειγμα Girsanov\n",
        "\n",
        "Στο Θεώρημα Girsanov εξετάζουμε αλλαγή μέτρου έτσι ώστε μια κίνηση Brown με drift\n",
        "να γίνει «καθαρή» κίνηση Brown κάτω από νέο μέτρο.\n",
        "\n",
        "Αν $(W_t)$ είναι κίνηση Brown υπό $\\mathbb{P}$ και $\\theta \\in \\mathbb{R}$, ορίζουμε\n",
        "τον εκθετικό μαρτινγάλη\n",
        "$$\n",
        "L_t = \\exp\\left( -\\theta W_t - \\frac{1}{2}\\theta^2 t \\right).\n",
        "$$\n",
        "\n",
        "Τότε, αν ορίσουμε νέο μέτρο $\\mathbb{Q}$ με\n",
        "$$\n",
        "\\frac{d\\mathbb{Q}}{d\\mathbb{P}}\\Big|_{\\mathcal{F}_t} = L_t,\n",
        "$$\n",
        "η διαδικασία\n",
        "$$\n",
        "\\widetilde{W}_t = W_t + \\theta t\n",
        "$$\n",
        "είναι κίνηση Brown υπό $\\mathbb{Q}$.\n",
        "\n",
        "Αριθμητικά, θα ελέγξουμε ότι\n",
        "$$\n",
        "\\mathbb{E}_{\\mathbb{Q}}[\\widetilde{W}_T] = 0\n",
        "\\quad\\Leftrightarrow\\quad\n",
        "\\mathbb{E}_{\\mathbb{P}}[L_T \\widetilde{W}_T] \\approx 0.\n",
        "$$\n"
      ],
      "metadata": {
        "id": "rmOHi8GMXrkS"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "theta = 1.0\n",
        "T = 1.0\n",
        "N = 1000\n",
        "n_paths = 100000  # μεγάλο για να φαίνεται καλά το αποτέλεσμα\n",
        "\n",
        "t, W = simulate_brownian_motion(T=T, N=N, n_paths=n_paths, seed=123)\n",
        "W_T = W[:, -1]\n",
        "\n",
        "# Εκθετικός μαρτινγάλη L_T\n",
        "L_T = np.exp(-theta * W_T - 0.5 * theta**2 * T)\n",
        "\n",
        "# Νέα διαδικασία W-tilde\n",
        "W_tilde_T = W_T + theta * T\n",
        "\n",
        "# Εμπειρικές προσδοκίες\n",
        "E_P_Wtilde = np.mean(W_tilde_T)\n",
        "E_P_L_Wtilde = np.mean(L_T * W_tilde_T)\n",
        "\n",
        "print(f\"E_P[W_tilde_T] ≈ {E_P_Wtilde:.4f} (δεν είναι 0)\")\n",
        "print(f\"E_P[L_T * W_tilde_T] ≈ {E_P_L_Wtilde:.4f} (πρέπει να ≈ 0)\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-_-AqvOyXuz5",
        "outputId": "47d9dadd-11e5-4a8a-dcd2-d2e559478d55"
      },
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "E_P[W_tilde_T] ≈ 1.0019 (δεν είναι 0)\n",
            "E_P[L_T * W_tilde_T] ≈ 0.0022 (πρέπει να ≈ 0)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Ευρωπαϊκό call option σε γεωμετρική κίνηση Brown\n",
        "\n",
        "Στην ενότητα 10.8.6 εξετάζεται η εφαρμογή των SDE και της φόρμουλας Itô στα δικαιώματα προαίρεσης.\n",
        "\n",
        "Υποθέτουμε ότι η τιμή μιας μετοχής $(S_t)$ υπό το **risk–neutral** μέτρο ικανοποιεί\n",
        "$$\n",
        "dS_t = r S_t\\, dt + \\sigma S_t\\, dW_t, \\qquad S_0 > 0,\n",
        "$$\n",
        "με λύση\n",
        "$$\n",
        "S_T = S_0 \\exp\\!\\left((r - \\tfrac{1}{2}\\sigma^2)T + \\sigma W_T\\right).\n",
        "$$\n",
        "\n",
        "Η τιμή ενός **Ευρωπαϊκού call** με τιμή εξάσκησης $K$ και λήξη $T$ είναι\n",
        "$$\n",
        "C_0 = e^{-rT} \\mathbb{E}^{\\mathbb{Q}}\\big[(S_T - K)^+\\big].\n",
        "$$\n",
        "\n",
        "Θα υπολογίσουμε τον $C_0$ με Monte Carlo:\n",
        "\n",
        "1. Προσομοιώνουμε πολλές τιμές $S_T$ από τη γεωμετρική κίνηση Brown.\n",
        "2. Υπολογίζουμε το μέσο των προεξοφλημένων payoffs $(S_T - K)^+$.\n"
      ],
      "metadata": {
        "id": "jQoQSjVHXz5W"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from math import log\n",
        "from math import sqrt\n",
        "\n",
        "# Παράμετροι\n",
        "S0 = 100.0\n",
        "K = 100.0\n",
        "r = 0.05\n",
        "sigma = 0.2\n",
        "T = 1.0\n",
        "n_paths = 200000\n",
        "\n",
        "# Προσομοίωση W_T ~ N(0, T)\n",
        "np.random.seed(0)\n",
        "W_T = np.sqrt(T) * np.random.randn(n_paths)\n",
        "\n",
        "# S_T under risk-neutral measure\n",
        "S_T = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * W_T)\n",
        "\n",
        "# Payoff και προεξόφληση\n",
        "payoffs = np.maximum(S_T - K, 0.0)\n",
        "C0_MC = np.exp(-r * T) * np.mean(payoffs)\n",
        "std_err = np.exp(-r * T) * np.std(payoffs) / np.sqrt(n_paths)\n",
        "\n",
        "print(f\"Monte Carlo τιμή call ≈ {C0_MC:.4f}\")\n",
        "print(f\"Εκτίμηση σφάλματος StdErr ≈ {std_err:.4f}\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "95yj7690X2ft",
        "outputId": "173fc895-48b0-401d-8ed4-0c9c49d7232e"
      },
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Monte Carlo τιμή call ≈ 10.4709\n",
            "Εκτίμηση σφάλματος StdErr ≈ 0.0329\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Τύπος Black–Scholes για την τιμή του Ευρωπαϊκού call\n",
        "\n",
        "Ο τύπος Black–Scholes δίνει κλειστή μορφή για την τιμή του Ευρωπαϊκού call:\n",
        "$$\n",
        "C_0 = S_0 N(d_1) - K e^{-rT} N(d_2),\n",
        "$$\n",
        "όπου\n",
        "$$\n",
        "d_1 = \\frac{\\ln(S_0/K) + (r + \\tfrac{1}{2}\\sigma^2)T}{\\sigma \\sqrt{T}},\n",
        "\\qquad\n",
        "d_2 = d_1 - \\sigma \\sqrt{T},\n",
        "$$\n",
        "και $N(\\cdot)$ είναι η αθροιστική συνάρτηση κατανομής της $\\mathcal{N}(0,1)$.\n",
        "\n",
        "Στο επόμενο κελί:\n",
        "- Υλοποιούμε τον τύπο Black–Scholes.\n",
        "- Τον συγκρίνουμε με την Monte Carlo εκτίμηση του προηγούμενου κελιού.\n"
      ],
      "metadata": {
        "id": "BPfJuf64X8CE"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from math import erf\n",
        "\n",
        "def N_cdf(x):\n",
        "    \"\"\"CDF της N(0,1) χωρίς scipy.\"\"\"\n",
        "    return 0.5 * (1.0 + erf(x / np.sqrt(2.0)))\n",
        "\n",
        "def black_scholes_call(S0, K, r, sigma, T):\n",
        "    d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))\n",
        "    d2 = d1 - sigma * np.sqrt(T)\n",
        "    return S0 * N_cdf(d1) - K * np.exp(-r * T) * N_cdf(d2)\n",
        "\n",
        "C0_BS = black_scholes_call(S0, K, r, sigma, T)\n",
        "\n",
        "print(f\"Τιμή Black–Scholes C0 = {C0_BS:.4f}\")\n",
        "print(f\"Τιμή Monte Carlo      ≈ {C0_MC:.4f} ± {1.96*std_err:.4f} (95% CI)\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "oTEFh_YzX9TH",
        "outputId": "cc189175-e444-4d07-d9a5-b9639484c313"
      },
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Τιμή Black–Scholes C0 = 10.4506\n",
            "Τιμή Monte Carlo      ≈ 10.4709 ± 0.0645 (95% CI)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Αριθμητική επίλυση στοχαστικής διαφορικής εξίσωσης με Euler–Maruyama\n",
        "\n",
        "Θεωρούμε τη στοχαστική διαφορική εξίσωση (SDE) στο $\\,[0,T]$:\n",
        "$$\n",
        "dX_t = f(t, X_t)\\, dt + g(t, X_t)\\, dW_t, \\quad X_0 = x_0,\n",
        "$$\n",
        "όπου $W_t$ είναι κίνηση Brown, $f$ είναι ο όρος \\emph{drift} και $g$ ο όρος \\emph{diffusion}.\n",
        "\n",
        "Σε διακριτό πλέγμα χρόνου\n",
        "$$\n",
        "0 = t_0 < t_1 < \\dots < t_N = T, \\quad \\Delta t = \\frac{T}{N},\n",
        "$$\n",
        "ο σχήμα Euler–Maruyama δίνεται από:\n",
        "$$\n",
        "X_{k+1} = X_k + f(t_k, X_k)\\, \\Delta t + g(t_k, X_k)\\, \\Delta W_k,\n",
        "$$\n",
        "όπου\n",
        "$$\n",
        "\\Delta W_k = W_{t_{k+1}} - W_{t_k} \\sim \\mathcal{N}(0,\\Delta t).\n",
        "$$\n",
        "\n",
        "Στον κώδικα που ακολουθεί:\n",
        "\n",
        "- Ο **χρήστης ορίζει** τις συναρτήσεις\n",
        "\n",
        "  - $f(t,x)$ για το drift\n",
        "  - $g(t,x)$ για το diffusion\n",
        "\n",
        "- Η συνάρτηση `euler_maruyama`:\n",
        "  - προσομοιώνει $n_{\\text{paths}}$ τροχιές,\n",
        "  - επιστρέφει το διάνυσμα χρόνου $t$ και τον πίνακα $X$ με τις προσομοιωμένες λύσεις.\n",
        "\n",
        "**Σημαντικό για χρήση σε NumPy:**\n",
        "\n",
        "- Οι συναρτήσεις `f(t, x)` και `g(t, x)` πρέπει να δουλεύουν με διανύσματα `x` (NumPy arrays).\n",
        "- Να χρησιμοποιούνται συναρτήσεις της `numpy` (π.χ. `np.sin`, `np.exp`, κτλ.) και όχι της `math`.\n",
        "\n",
        "Παράδειγμα SDE (Ornstein–Uhlenbeck):\n",
        "$$\n",
        "dX_t = -a X_t\\, dt + \\sigma\\, dW_t.\n",
        "$$\n",
        "Σε αυτή την περίπτωση\n",
        "$$\n",
        "f(t,x) = -a x, \\quad g(t,x) = \\sigma.\n",
        "$$\n"
      ],
      "metadata": {
        "id": "DlAM2j0eYAmA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "def euler_maruyama(f, g, x0, T=1.0, N=1000, n_paths=1, seed=None):\n",
        "    \"\"\"\n",
        "    Euler–Maruyama για SDE:\n",
        "        dX_t = f(t, X_t) dt + g(t, X_t) dW_t,  X_0 = x0\n",
        "\n",
        "    Παραμέτρους:\n",
        "    - f : συνάρτηση f(t, x) (drift)\n",
        "    - g : συνάρτηση g(t, x) (diffusion)\n",
        "    - x0 : αρχική τιμή X_0 (scalar)\n",
        "    - T : τελικός χρόνος\n",
        "    - N : αριθμός βημάτων (Δt = T/N)\n",
        "    - n_paths : πόσες ανεξάρτητες τροχιές θα προσομοιώσουμε\n",
        "    - seed : (optional) ακέραιος για αναπαραγωγιμότητα\n",
        "\n",
        "    Επιστρέφει:\n",
        "    - t : διάνυσμα χρόνου σχήματος (N+1,)\n",
        "    - X : πίνακας λύσεων σχήματος (n_paths, N+1)\n",
        "    \"\"\"\n",
        "    if seed is not None:\n",
        "        np.random.seed(seed)\n",
        "\n",
        "    dt = T / N\n",
        "    t = np.linspace(0.0, T, N + 1)\n",
        "\n",
        "    # Πίνακας για τις τροχιές\n",
        "    X = np.zeros((n_paths, N + 1))\n",
        "    X[:, 0] = x0\n",
        "\n",
        "    # Ανεξάρτητες αυξήσεις Brownian: ΔW_k ~ N(0, dt)\n",
        "    dW = np.sqrt(dt) * np.random.randn(n_paths, N)\n",
        "\n",
        "    # Βρόχος Euler–Maruyama\n",
        "    for k in range(N):\n",
        "        tk = t[k]\n",
        "        Xk = X[:, k]          # σχήμα (n_paths,)\n",
        "        drift = f(tk, Xk)     # f(t_k, X_k)\n",
        "        diff  = g(tk, Xk)     # g(t_k, X_k)\n",
        "        X[:, k+1] = Xk + drift * dt + diff * dW[:, k]\n",
        "\n",
        "    return t, X\n",
        "\n",
        "\n",
        "# ============================================\n",
        "# Παράδειγμα χρήσης: Ornstein–Uhlenbeck SDE\n",
        "# dX_t = -a X_t dt + sigma dW_t\n",
        "# ============================================\n",
        "\n",
        "a = 2.0\n",
        "sigma = 0.5\n",
        "\n",
        "def f_example(t, x):\n",
        "    # drift f(t,x) = -a x\n",
        "    return -a * x\n",
        "\n",
        "def g_example(t, x):\n",
        "    # diffusion g(t,x) = sigma\n",
        "    # ΠΡΟΣΟΧΗ: πρέπει να επιστρέψει array ίδιου σχήματος με το x\n",
        "    return sigma * np.ones_like(x)\n",
        "\n",
        "# Παράμετροι προσομοίωσης\n",
        "x0 = 0.0\n",
        "T = 2.0\n",
        "N = 2000\n",
        "n_paths = 5\n",
        "\n",
        "t, X = euler_maruyama(f_example, g_example, x0, T=T, N=N, n_paths=n_paths, seed=123)\n",
        "\n",
        "plt.figure(figsize=(8, 4))\n",
        "for i in range(n_paths):\n",
        "    plt.plot(t, X[i], alpha=0.8)\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"X_t\")\n",
        "plt.title(\"Παράδειγμα SDE με Euler–Maruyama (Ornstein–Uhlenbeck)\")\n",
        "plt.grid(True)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 353
        },
        "id": "0ZpaI_XeZK1p",
        "outputId": "ff57a37a-0348-469b-8205-6f36344820ef"
      },
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 137076 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Πιο πολύπλοκο παράδειγμα: Στοχαστική λογιστική εξίσωση (logistic SDE)\n",
        "\n",
        "Θεωρούμε τη στοχαστική διαφορική εξίσωση (SDE):\n",
        "$$\n",
        "dX_t = r X_t \\left(1 - \\frac{X_t}{K}\\right)\\, dt + \\sigma X_t\\, dW_t,\n",
        "\\quad X_0 = x_0,\n",
        "$$\n",
        "όπου:\n",
        "\n",
        "- $r > 0$ είναι ο ρυθμός ανάπτυξης,\n",
        "- $K > 0$ είναι η χωρητικότητα (carrying capacity),\n",
        "- $\\sigma \\ge 0$ η ένταση του θορύβου.\n",
        "\n",
        "Ο όρος\n",
        "$$\n",
        "f(t,x) = r x \\left(1 - \\frac{x}{K}\\right)\n",
        "$$\n",
        "είναι μη γραμμικός drift (λογιστική ανάπτυξη) και\n",
        "$$\n",
        "g(t,x) = \\sigma x\n",
        "$$\n",
        "είναι diffusion **ανάλογο της κατάστασης** (πολλαπλασιαστικός θόρυβος).\n",
        "\n",
        "Αν $\\sigma = 0$, παίρνουμε την κλασική ντετερμινιστική λογιστική εξίσωση:\n",
        "$$\n",
        "\\frac{dX_t}{dt} = r X_t \\left(1 - \\frac{X_t}{K}\\right),\n",
        "$$\n",
        "η οποία έχει ρητή λύση\n",
        "$$\n",
        "X_t^{\\text{det}} =\n",
        "\\frac{K x_0 e^{rt}}{K + x_0 (e^{rt} - 1)}.\n",
        "$$\n",
        "\n",
        "Σχημα Euler–Maruyama σε διακριτό πλέγμα:\n",
        "$$\n",
        "X_{k+1} = X_k\n",
        "+ r X_k\\left(1 - \\frac{X_k}{K}\\right)\\Delta t\n",
        "+ \\sigma X_k \\Delta W_k,\n",
        "$$\n",
        "με\n",
        "$$\n",
        "\\Delta W_k \\sim \\mathcal{N}(0,\\Delta t).\n",
        "$$\n",
        "\n",
        "Στον κώδικα που ακολουθεί:\n",
        "\n",
        "- Ορίζουμε τις συναρτήσεις\n",
        "  - $f(t,x) = r x (1 - x/K)$,\n",
        "  - $g(t,x) = \\sigma x$.\n",
        "- Χρησιμοποιούμε τη γενική ρουτίνα `euler_maruyama(f, g, x0, ...)`.\n",
        "- Σχεδιάζουμε:\n",
        "  - πολλές στοχαστικές τροχιές $X_t$,\n",
        "  - μαζί με τη ντετερμινιστική λύση $X_t^{\\text{det}}$ για σύγκριση.\n"
      ],
      "metadata": {
        "id": "CwrsGXxfaEWW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Παράμετροι της logistic SDE\n",
        "r = 1.0      # ρυθμός ανάπτυξης\n",
        "K = 2.0      # χωρητικότητα\n",
        "sigma = 0.2  # ένταση θορύβου\n",
        "\n",
        "def f_logistic(t, x):\n",
        "    \"\"\"\n",
        "    Drift: f(t,x) = r x (1 - x/K)\n",
        "    Δουλεύει για x ως numpy array (σχήμα (n_paths,))\n",
        "    \"\"\"\n",
        "    return r * x * (1 - x / K)\n",
        "\n",
        "def g_logistic(t, x):\n",
        "    \"\"\"\n",
        "    Diffusion: g(t,x) = sigma x\n",
        "    Επιστρέφει array ίδιου σχήματος με το x.\n",
        "    \"\"\"\n",
        "    return sigma * x\n",
        "\n",
        "# Παράμετροι προσομοίωσης\n",
        "x0 = 0.2\n",
        "T = 5.0\n",
        "N = 5000\n",
        "n_paths = 10\n",
        "\n",
        "# Προσομοίωση με Euler–Maruyama\n",
        "t, X = euler_maruyama(f_logistic, g_logistic, x0, T=T, N=N, n_paths=n_paths, seed=42)\n",
        "\n",
        "# Ντετερμινιστική λογιστική λύση (για σύγκριση, χωρίς θόρυβο)\n",
        "def logistic_deterministic(t, r, K, x0):\n",
        "    # κλασική ρητή λύση της logistic ODE\n",
        "    return (K * x0 * np.exp(r * t)) / (K + x0 * (np.exp(r * t) - 1))\n",
        "\n",
        "X_det = logistic_deterministic(t, r, K, x0)\n",
        "\n",
        "# Γραφική παράσταση\n",
        "plt.figure(figsize=(9, 4))\n",
        "\n",
        "# Στοχαστικές τροχιές\n",
        "for i in range(n_paths):\n",
        "    plt.plot(t, X[i], alpha=0.7)\n",
        "\n",
        "# Ντετερμινιστική λύση\n",
        "plt.plot(t, X_det, \"k\", linewidth=2.5, label=\"Ντετερμινιστική λύση (σ = 0)\")\n",
        "\n",
        "plt.xlabel(\"t\")\n",
        "plt.ylabel(\"X_t\")\n",
        "plt.title(\"Στοχαστική λογιστική εξίσωση (Euler–Maruyama)\")\n",
        "plt.axhline(K, linestyle=\"--\", label=\"Χωρητικότητα K\")\n",
        "plt.grid(True)\n",
        "plt.legend()\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 324
        },
        "id": "jWQyQukIaFPn",
        "outputId": "21eb1661-50c5-4fd1-de50-e85ee82a29f6"
      },
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 900x400 with 1 Axes>"
            ],
            "image/png": "[encoded data, 153548 characters]\n"
          },
          "metadata": {}
        }
      ]
    }
  ]
}