Implicit_FDM.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[5.17924845]\n",
      " [0.99937502]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from numpy import linalg as la\n",
    "\n",
    "def compute_transition_probabilities(M: int, dt: float, r: float, q: float, sigma: float):\n",
    "    \"\"\"Calculate the transition probabilities for the finite difference matrix.\"\"\"\n",
    "    j = np.arange(0, M, dtype=np.float64)\n",
    "    a = 0.5 * (r - q) * j * dt - 0.5 * sigma**2 * j**2 * dt\n",
    "    b = 1.0 + sigma**2 * j**2 * dt + r * dt\n",
    "    c = -0.5 * (r - q) * j * dt - 0.5 * sigma**2 * j**2 * dt\n",
    "    return a, b, c\n",
    "\n",
    "def create_tridiagonal_matrix(a: np.ndarray, b: np.ndarray, c: np.ndarray):\n",
    "    \"\"\"Construct and invert the tridiagonal matrix from transition probabilities.\"\"\"\n",
    "    M = len(b)\n",
    "    A = np.diag(b) + np.diag(a[1:], k=-1) + np.diag(c[:-1], k=1)\n",
    "    B = la.inv(A)\n",
    "    return B, la.norm(B, np.inf)\n",
    "\n",
    "def initialize_option_matrix(N: int, M: int, S_max: float, dS: float, K: float, r: float, dt: float, is_call: bool):\n",
    "    \"\"\"Set up the initial and boundary conditions for the option pricing matrix.\"\"\"\n",
    "    FV = np.zeros((N + 1, M + 1))\n",
    "    if is_call:\n",
    "        FV[:, 0] = 0.0\n",
    "        FV[:, M] = S_max * np.exp(-r * (N - np.arange(N + 1)) * dt)\n",
    "        FV[N, :] = np.maximum(np.arange(0, S_max + dS, dS) - K, 0)\n",
    "    else:\n",
    "        FV[:, 0] = K * np.exp(-r * (N - np.arange(N + 1)) * dt)\n",
    "        FV[:, M] = 0.0\n",
    "        FV[N, :] = np.maximum(K - np.arange(0, S_max + dS, dS), 0)\n",
    "    return FV\n",
    "\n",
    "def backward_iteration(FV: np.ndarray, B: np.ndarray, a: np.ndarray, c: np.ndarray, N: int, M: int, dS: float, K: float, is_call: bool, is_american: bool):\n",
    "    \"\"\"Perform backward iteration to calculate the option price at each time step.\"\"\"\n",
    "    for i in range(N - 1, -1, -1):\n",
    "        k_i = np.zeros(M)\n",
    "        k_i[0] = -a[0] * FV[i + 1, 0]\n",
    "        k_i[M - 1] = -c[M - 1] * FV[i + 1, M]\n",
    "        \n",
    "        FV[i, :M] = np.dot(FV[i + 1, :M] + k_i, B)\n",
    "\n",
    "        if is_american:\n",
    "            intrinsic_value = np.maximum(np.arange(0, M + 1) * dS - K if is_call else K - np.arange(0, M + 1) * dS, FV[i, :])\n",
    "            FV[i, :] = intrinsic_value\n",
    "    return FV\n",
    "\n",
    "def vanilla_option_implicit(S_0: float, T: float, K: float, r: float, q: float, sigma: float, N: int, M: int, is_call: bool, is_american: bool) -> np.ndarray:\n",
    "    \"\"\"Calculate the fair value of a vanilla European or American call/put option.\"\"\"\n",
    "    results = np.zeros((2, 1), dtype=float)\n",
    "    S_max = 2.0 * K\n",
    "    dS = S_max / M\n",
    "    dt = T / N\n",
    "\n",
    "    # Compute transition probabilities\n",
    "    a, b, c = compute_transition_probabilities(M, dt, r, q, sigma)\n",
    "    \n",
    "    # Create tridiagonal matrix and its inverse\n",
    "    B, infinity_norm = create_tridiagonal_matrix(a, b, c)\n",
    "    results[1] = infinity_norm\n",
    "\n",
    "    # Initialize option price matrix with boundary conditions\n",
    "    FV = initialize_option_matrix(N, M, S_max, dS, K, r, dt, is_call)\n",
    "\n",
    "    # Perform backward iteration\n",
    "    FV = backward_iteration(FV, B, a, c, N, M, dS, K, is_call, is_american)\n",
    "    \n",
    "    # Store fair value of option in results\n",
    "    results[0] = FV[0, M // 2]\n",
    "    return results\n",
    "\n",
    "# Sample test execution\n",
    "if __name__ == \"__main__\":\n",
    "    S_0 = 50.0\n",
    "    T = 1.0\n",
    "    K = 50.0\n",
    "    r = 0.1\n",
    "    q = 0.00\n",
    "    sigma = 0.25\n",
    "    N = 100\n",
    "    M = 40\n",
    "    is_call = True\n",
    "    is_american = False\n",
    "    print(vanilla_option_implicit(S_0, T, K, r, q, sigma, N, M, is_call, is_american))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import unittest\n",
    "\n",
    "# Define the utility functions and main function here (or import them if they are in a module)\n",
    "\n",
    "# Define the Test Cases\n",
    "class TestVanillaOptionImplicit(unittest.TestCase):\n",
    "    # Test methods go here as defined previously\n",
    "    def test_compute_transition_probabilities(self):\n",
    "        M, dt, r, q, sigma = 40, 0.01, 0.1, 0.02, 0.25\n",
    "        a, b, c = compute_transition_probabilities(M, dt, r, q, sigma)\n",
    "        self.assertEqual(len(a), M)\n",
    "        self.assertEqual(len(b), M)\n",
    "        self.assertEqual(len(c), M)\n",
    "        self.assertAlmostEqual(a[1], 8.750000000000001e-05, places=4)\n",
    "        self.assertAlmostEqual(b[1], 1.001625, places=4)\n",
    "        self.assertAlmostEqual(c[1], -0.0007125, places=4)\n",
    "\n",
    "    def test_create_tridiagonal_matrix(self):\n",
    "        M = 40\n",
    "        a = np.full(M, -0.5)\n",
    "        b = np.full(M, 2.0)\n",
    "        c = np.full(M, -0.5)\n",
    "        B, infinity_norm = create_tridiagonal_matrix(a, b, c)\n",
    "        self.assertEqual(B.shape, (M, M))\n",
    "        self.assertGreater(infinity_norm, 0)\n",
    "        self.assertTrue(np.allclose(B @ (np.diag(b) + np.diag(a[1:], k=-1) + np.diag(c[:-1], k=1)), np.eye(M)))\n",
    "\n",
    "    def test_initialize_option_matrix(self):\n",
    "        N, M = 100, 40\n",
    "        S_max, K, r, dt = 100, 50, 0.1, 0.01\n",
    "        dS = S_max / M\n",
    "        FV_call = initialize_option_matrix(N, M, S_max, dS, K, r, dt, is_call=True)\n",
    "        self.assertEqual(FV_call.shape, (N + 1, M + 1))\n",
    "        self.assertAlmostEqual(FV_call[N, 0], 0.0)\n",
    "        self.assertTrue(np.all(FV_call[N, :] >= 0))\n",
    "        FV_put = initialize_option_matrix(N, M, S_max, dS, K, r, dt, is_call=False)\n",
    "        self.assertEqual(FV_put.shape, (N + 1, M + 1))\n",
    "        self.assertAlmostEqual(FV_put[N, M], 0.0)\n",
    "        self.assertTrue(np.all(FV_put[N, :] >= 0))\n",
    "\n",
    "    def test_backward_iteration(self):\n",
    "        N, M = 5, 5\n",
    "        S_max, K, dS = 100, 50, 20\n",
    "        r, dt = 0.1, 0.01\n",
    "        is_call = True\n",
    "        is_american = False\n",
    "        a = np.full(M, -0.5)\n",
    "        b = np.full(M, 2.0)\n",
    "        c = np.full(M, -0.5)\n",
    "        B, _ = create_tridiagonal_matrix(a, b, c)\n",
    "        FV = initialize_option_matrix(N, M, S_max, dS, K, r, dt, is_call)\n",
    "        FV_result = backward_iteration(FV, B, a, c, N, M, dS, K, is_call, is_american)\n",
    "        self.assertEqual(FV_result.shape, (N + 1, M + 1))\n",
    "        self.assertTrue(np.all(FV_result >= 0))\n",
    "\n",
    "    def test_vanilla_option_implicit(self):\n",
    "        S_0, T, K, r, q, sigma, N, M = 50.0, 1.0, 50.0, 0.1, 0.0, 0.25, 100, 40\n",
    "        is_call, is_american = True, False\n",
    "        result = vanilla_option_implicit(S_0, T, K, r, q, sigma, N, M, is_call, is_american)\n",
    "        self.assertEqual(result.shape, (2, 1))\n",
    "        option_value = result[0, 0]\n",
    "        self.assertGreater(option_value, 0)\n",
    "        self.assertLess(option_value, S_0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      ".....\n",
      "----------------------------------------------------------------------\n",
      "Ran 5 tests in 0.010s\n",
      "\n",
      "OK\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<unittest.runner.TextTestResult run=5 errors=0 failures=0>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unittest.TextTestRunner().run(unittest.TestLoader().loadTestsFromTestCase(TestVanillaOptionImplicit))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}