Math623Homework3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d3f5d89-b649-42ea-9427-91cbdadd4442",
   "metadata": {},
   "source": [
    "# Problem 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 306,
   "id": "a2a4ca89-a1f6-4360-a537-eb9721603844",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "id": "25725630-5592-4a02-b2d4-4affb453e74e",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = 0.9\n",
    "lambd = 0.25\n",
    "mew = 0.5\n",
    "n = 100000\n",
    "m = int(n/2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "id": "321e79ca-253d-4b82-a7b6-ef7125830613",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = list()\n",
    "Y = list()\n",
    "\n",
    "for i in range(n):\n",
    "    V = np.random.uniform()\n",
    "    W = np.random.uniform()\n",
    "    \n",
    "    U = fsolve(lambda u: theta*(u-1)*u*(2*V-1) + u - W,0)[0]\n",
    "    \n",
    "    X.append(-np.log(1-U)/lambd)\n",
    "    Y.append(-np.log(1-V)/mew)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 346,
   "id": "534aa2db-d164-4812-8351-611c9ce9e4be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.22194739],\n",
       "       [0.22194739, 1.        ]])"
      ]
     },
     "execution_count": 346,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.corrcoef(X,Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "id": "2a7dae2a-34f6-45ff-8039-e4ee49cf987f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The correlation coefficient between X and Y is:  0.2219473929907042\n"
     ]
    }
   ],
   "source": [
    "print(\"The correlation coefficient between X and Y is: \", np.corrcoef(X,Y)[0][1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f716e98-062d-4d8b-a953-b34dfcd07033",
   "metadata": {},
   "source": [
    "# Problem 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "05d40e67-a69d-4372-a9c5-89da2668ca9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 1\n",
    "B = np.sqrt(2/np.pi) * 1/(1-math.erf(A/np.sqrt(2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "e1b75d03-8b3b-42f1-a479-ecd97fe2c031",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 0\n",
    "c = 2\n",
    "X = list()\n",
    "numSamples = 1000\n",
    "count=0\n",
    "\n",
    "while n < numSamples:\n",
    "    Y = -np.log(1-np.random.uniform())/A + A\n",
    "    U = np.random.uniform()\n",
    "    \n",
    "    ratio = B*np.exp(-Y**2 / 2) / (c*(A + np.exp(-Y)/A))\n",
    "    \n",
    "    if U <= ratio:\n",
    "        X.append(Y)\n",
    "        n+=1\n",
    "    \n",
    "    count+=1\n",
    "    \n",
    "    if ratio > 1:\n",
    "        print(\"Make c bigger\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "e5809701-e1a1-4308-b602-4a2b9a09789f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/sethkelley/opt/anaconda3/lib/python3.9/site-packages/seaborn/distributions.py:2619: FutureWarning: `distplot` is a deprecated function and will be removed in a future version. Please adapt your code to use either `displot` (a figure-level function with similar flexibility) or `histplot` (an axes-level function for histograms).\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-0.25, 2.0)"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD8CAYAAABZ/vJZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5i0lEQVR4nO3dd3hc9ZXw8e+Z0ahLVrVsNfduXGUbCKHGxpQAAZYYk1AS1puEbLIhBbJkQ5Ykz7KbN7t505YWEpKXkkKJARtCDZDgCtjGDRtXybIlS1axujTn/eM3wmNZsmVJozuSzud55pmZe3/33nNnpHvmV+69oqoYY4wxp8vndQDGGGMGJksgxhhjesQSiDHGmB6xBGKMMaZHLIEYY4zpEUsgxhhjesSzBCIiBSLymohsEZHNIvLVTsqIiPxURHaKyEYRmRM27yYR2RF63NS/0RtjjBGvzgMRkZHASFV9R0RSgPXAVaq6JazMpcA/A5cCC4D/q6oLRCQDWAcUARpadq6qHunv/TDGmKHKsxqIqpaq6juh17XAViCvQ7Ergd+qswpICyWei4GXVLUylDReAhb3Y/jGGDPkxXgdAICIjAZmA6s7zMoD9oe9Lw5N62p6Z+teBiwDSEpKmjt58uS+CdoYY4aI9evXH1bV7I7TPU8gIpIMPAn8i6rW9PX6VfUB4AGAoqIiXbduXV9vwhhjBjUR2dvZdE9HYYlIAJc8HlXVpzopUgIUhL3PD03raroxxph+4uUoLAF+BWxV1f/uothy4MbQaKwzgWpVLQVeBBaJSLqIpAOLQtOMMcb0Ey+bsD4GfBbYJCLvhab9K1AIoKr3AStwI7B2AvXALaF5lSLyfWBtaLl7VLWy/0I3xhjjWQJR1bcAOUUZBW7rYt7DwMMRCM0YY/pUS0sLxcXFNDY2eh3KScXHx5Ofn08gEOhWec870U3nHlu975Rlli4o7IdIjDG9VVxcTEpKCqNHj8a13kcfVaWiooLi4mLGjBnTrWXsUibGGBNhjY2NZGZmRm3yABARMjMzT6uWZAnEGGP6QTQnj3anG6MlEGOMMT1iCcQYY4aIF154gUmTJjF+/HjuvffeXq/PEogxxgwBbW1t3HbbbaxcuZItW7bw+OOPs2XLllMveBKWQIwxJto8+iiMHg0+n3t+9NFer3LNmjWMHz+esWPHEhsby5IlS/jzn//cq3VaAjHGmGjy6KOwbBns3Quq7nnZsl4nkZKSEgoKjl0BKj8/n5KS3l0ByhKIMcZEk7vugvr646fV17vpUcYSiDHGRJN9XZxE3NX0bsrLy2P//mN3wSguLiYvr9O7YHSbJRBjjIkmhV1cYaKr6d00b948duzYwe7du2lubuaJJ57giiuu6NU6LYEYY0w0+eEPITHx+GmJiW56L8TExPDzn/+ciy++mClTpnDdddcxbdq03q2zV0sbY4zpWzfc4J7vuss1WxUWuuTRPr0XLr30Ui699NJer6edJRBjjIk2N9zQJwkj0qwJyxhjTI94WgMRkYeBy4EyVZ3eyfxvAu1pOAaYAmSHbii1B6gF2oBWVS3qn6iNMeb0qWrUX1DR3YKp+7yugfwGWNzVTFX9karOUtVZwLeBv3a48+AFofmWPIwxUSs+Pp6KiorTPkD3p/b7gcTHx3d7GU9rIKr6hoiM7mbx64HHIxiOMcZERH5+PsXFxZSXl3sdykm135GwuwZEJ7qIJOJqKl8Om6zAX0REgftV9QFPgjPGmFMIBALdvsvfQDIgEgjwSeBvHZqvzlHVEhEZDrwkIttU9Y2OC4rIMmAZQGEvT8QxxhhzjNd9IN21hA7NV6paEnouA54G5ne2oKo+oKpFqlqUnZ0d8UCNMWaoiPoEIiLDgPOAP4dNSxKRlPbXwCLgfW8iNMaYocnrYbyPA+cDWSJSDNwNBABU9b5QsU8Bf1HVurBFc4CnQ0PiYoDHVPWF/oo7Wjy2+tQXV1u6wJrtjDGR4fUorOu7UeY3uOG+4dN2ATMjE5UxxpjuiPomLGOMMdHJEogxxpgesQRijDGmRyyBGGOM6RFLIMYYY3rEEogxxpgesQRijDGmRyyBGGOM6RFLIMYYY3rEEogxxpgesQRijDGmRyyBGGOM6ZGBckMp0w2twSCrdlVypL6ZlLgYFozJ9DokY8wgZglkkKisa+b/rdrLwZpG4mJ8NLUGeXtXBZNGpHDOhCyvwzPGDELWhDUItAWVJ9buo6qhmc+eOYq7PzmNL50/joSAn88/spb1e494HaIxZhCyBDIIvLmjnOIjDVw1K48pI1MByE9P5NaPj2XksHhufWQtJVUNHkdpjBlsPE0gIvKwiJSJSKe3oxWR80WkWkTeCz2+GzZvsYhsF5GdInJn/0UdXeqaWnl1WxnT84YxIz/tuHnJcTH8+pb5NLUG+fZTm1BVb4I0xgxKXtdAfgMsPkWZN1V1VuhxD4CI+IFfAJcAU4HrRWRqRCONUuv2VNIaVC6aPLzT+WOykrhj8WTe+KCcP64v7ufojDGDmacJRFXfACp7sOh8YKeq7lLVZuAJ4Mo+DW4AaAsqq3ZXMi47iZzU+C7LffbMUcwbnc69K7dR09jSjxEaYwYzr2sg3XGWiGwQkZUiMi00LQ/YH1amODTtBCKyTETWici68vLySMfar7YdrKG6oYWzxp58uK7PJ3z38mlU1jVz3+sf9lN0xpjBLtoTyDvAKFWdCfwMeOZ0V6CqD6hqkaoWZWdn93V8ntpYXE1yXAyTQx3nJ3NG/jCumpXLr97aTWm1dagbY3ovqhOIqtao6tHQ6xVAQESygBKgIKxofmjakNEaDPLBoVomj0jBJ9KtZb6+aBJtQeX+v+6KcHTGmKEgqhOIiIwQcUdHEZmPi7cCWAtMEJExIhILLAGWexdp/9tzuJ6m1uBHw3a7oyAjkatm5/HE2n0cPtoUweiMMUOBp2eii8jjwPlAlogUA3cDAQBVvQ+4FviiiLQCDcASdWNRW0Xky8CLgB94WFU3e7ALntlaWkPAL4zLTj5pucdW7zvufUF6Ik0tQb7xhw0smjYCgKULCiMWpzFm8PI0gajq9aeY/3Pg513MWwGsiERc0U5V2XqwhvHZycTGnF4lMjsljml5w3h7VwXnTswmPuCPUJTGmMEuqpuwTOcq65qpqm9h4oiUHi1//sRsmlqDrNpV0ceRGWOGEksgA9CeinoARmcm9Wj53LQEJuWk8NbOwzS3BvsyNGPMEGIJZADaU1FHQsBPdkpcj9dx/qRs6pvbWL+3J+dxGmOMJZABac/hOkZnJnZ7+G5nRmUmUZCewN8+rKAtaNfIMsacPksgA0xtYwsVdc2MzupZ81W4cyZkU1nXzMtbD/VBZMaYocYSyADT2/6PcFNHppKeGOBXb+7u9bqMMUOPJZABZl9FHQG/kJuW0Ot1+X3C2eOyWLOnkg37q3ofnDFmSLEEMsCUVDUyclgCfl/P+z/CFY1KJyUuhofeslqIMeb0WAIZQIKqlFY3kJvW9aXbT1dcwM/1CwpZsanU7lpojDktlkAGkMqjzTS1Bskd1vvmq3A3nT0agN/8zWohxpjuswQygLTXEPLS+zaB5KUlcNkZI3lizX5q7YZTxphusgQygByoasDvE4an9F0TVrtbPz6G2qZWfr92/6kLG2MMHl9M0ZyekuoGRqTG91kHerv2K/aOzkzi56/tJC7Gf8I27Iq9xpiOrAYyQKgqB6oayOuD4btdOWd8FlX1LWw+UB2xbRhjBg9LIANEVUMLjS1BRvbhCKyOJo9MITMplrd2HsbddsUYY7rmaQIRkYdFpExE3u9i/g0islFENonI30VkZti8PaHp74nIuv6L2huHahoBGJEauQTiE+Fj47MoPtLAvsr6iG3HGDM4eF0D+Q2w+CTzdwPnqeoZwPeBBzrMv0BVZ6lqUYTiixqHatwtaCPRgR5uTmE6CQE/b+08HNHtGGMGPk8TiKq+AXR5PXFV/buqHgm9XQXk90tgUaisppHU+BgSYiN7B8HYGB/zx2Sw5UANFXbfdGPMSXhdAzkdnwdWhr1X4C8isl5ElnW1kIgsE5F1IrKuvLw84kFGyqGaRnIi2HwV7qyxmfhE+PuHdsdCY0zXBkQCEZELcAnkjrDJ56jqHOAS4DYRObezZVX1AVUtUtWi7Ozsfoi27wVVKatt6rcEkpoQYEb+MNbvPUJDc1u/bNMYM/BEfQIRkRnAQ8CVqvrRT2JVLQk9lwFPA/O9iTDyKuuaaQ0qOak9vwPh6TpnQhbNbUHW7LE7FhpjOhfVCURECoGngM+q6gdh05NEJKX9NbAI6HQk12DQPgKrv2ogACOHJTAuO4m3PzxMa9Dum26MOZGnZ6KLyOPA+UCWiBQDdwMBAFW9D/gukAn8UtztW1tDI65ygKdD02KAx1T1hX7fgX7SnkAiPQKro3PGZ/HI23vZVGwnFhpjTuRpAlHV608x/1bg1k6m7wJmnrjE4FRW20R6YoDYmP6tME7ISWFEajyvbS+jtS1IjD+qK6zGmH5mR4QB4HBtE9kp/df/0c4nwkVThnP4aDPPvHeg37dvjIlulkCiXFCV8qNNZCX3fwIBd9/03LR4fvrKDlrarC/EGHOMJZAoV9PQQkubelIDARARPjElh32V9fxpfbEnMRhjopMlkCh3+GgzgGc1EIBJOSnMKkjjZ6/soKnVzgsxxjiWQKJcea0bgeVVDQRcLeTriyZyoLqRJ9bYDaeMMY4lkChXfrSZuBgfKXHe3vvrnPFZnDU2k5+8/AFV9c2exmKMiQ6WQKLc4VrXgR4658UzIsJ3PzmV6oYWfvLyDk9jMcZEB0sgUa78qDdDeDszZWQqSxcU8rtVe9lxqNbrcIwxHrMEEsWaW4NUN7R42oHe0e0LJ5EU6+ee57bYXQuNGeIsgUSxijp3P46s5FiPIzkmIymWry2cyJs7DvPSlkNeh2OM8ZAlkChWEQVDeDvzmTNHMTEnme8t38zRplavwzHGeMQSSBSrqHMJJCMpemogAAG/j3uvmUFpTSP/uXKb1+EYYzxiCSSKVRxtIikuhvhAZG9j2xNzCtO55ewx/G7VXtbstnuGGDMUWQKJYhV1zWRGWe0j3DcunkhBRgJ3PLmRxhY7Q92YoaZbCUREnhKRy0TEEk4/qozyBJIYG8O9V89g9+E6/uelD069gDFmUOluQvglsBTYISL3isikvti4iDwsImUi0undBMX5qYjsFJGNIjInbN5NIrIj9LipL+KJJi1tbghvZhSNwOrMx8Zncf38Ah58cxfr91pTljFDSbcSiKq+rKo3AHOAPcDLIvJ3EblFRAK92P5vgMUnmX8JMCH0WAb8L4CIZODuXrgAdy/0u0UkvRdxRJ3KUAd6ZlJ0jcDqzF2XTSU3LYGv/2ED9c02KsuYoaLbF1gSkUzgM8BngXeBR4FzgJtwt6U9bar6hoiMPkmRK4HfqjtjbZWIpInIyND2XlLVylBsL+ES0eM9iSMafZRAorwGApAcF8OPrp3J9Q+u4vO/WccnZ+Z2WXbpgsJ+jMwYE0nd7QN5GngTSAQ+qapXqOrvVfWfgeQIxpcHhF/+tTg0ravpJxCRZSKyTkTWlZeXRyzQvlZx1J1EOBBqIABnjcvk7HGZvL2rgg/Lj3odjjGmH3S3BvKgqq4InyAicarapKpFEYirz6jqA8ADAEVFRQPm2hsVdc0kBPwkxEbHEN7HVu87ZZlFU0fwwaFanlxfzFcumhCVw4+NMX2nu53oP+hk2tt9GUgXSoCCsPf5oWldTR80KuqaB0TzVbjYGB/Xzi2guqGFFZtKvQ7HGBNhJ00gIjJCROYCCSIyW0TmhB7n45qzIm05cGNoNNaZQLWqlgIvAotEJD3Ueb4oNG3QqDjaFHVnoHdHYUYi507MZt3eI2w7WON1OMaYCDpVE9bFwM24X/j/HTa9FvjX3m5cRB7HdYhniUgxbmRVAEBV7wNWAJcCO4F64JbQvEoR+T6wNrSqe9o71AeD5tYgVfUtzCoYGP0fHV00eTjbD9by9DslfPWiRBI9vhmWMSYyTvqfraqPAI+IyDWq+mRfb1xVrz/FfAVu62Lew8DDfR1TNCg+Uo8yMEZgdSbG7+Paufn88vWdLN94gCXzbOSVMYPRSROIiHxGVf8fMFpEbu84X1X/u5PFTC/tragHiOqz0E8lNy2BCyfn8PLWQ0wdWcWM/DSvQzLG9LFTdaInhZ6TgZROHiYC9lTUAZAZZZdxP13nTcwmPz2B5RsOUNvY4nU4xpg+dqomrPtDz//eP+EYcDWQuBgfSVEyhLen/D7h2jn5/Oy1nax8/yDXFRWceiFjzIDR3RMJ/0tEUkUkICKviEi5iHwm0sENVXsq6shMikVEvA6l14anxnPuhGze219lJxgaM8h09zyQRapaA1yOuxbWeOCbkQpqqNtbUU/GAG++Cnf+pGwykmJZ/t4BmluDXodjjOkj3U0g7U1dlwF/VNXqCMUz5LW2BdlfWT+gO9A7Cvh9fHJGLuVHm3jorV1eh2OM6SPdTSDPicg2YC7wiohkA42RC2voOlDVSGtQB1UCAZg0IoVpuan89JUd7K+s9zocY0wf6O7l3O8EzgaKVLUFqMNdKdf0scEyAqszl50xEkH4zxfsPurGDAanc4rwZNz5IOHL/LaP4xny9rYnkEFWAwFIS4zlrHGZPLexlPy0rRRmJp1Qxi73bszA0d1RWL8D/g/u/h/zQo+ovgrvQLWnop74gI+U+MF5+Y+PT8giJT6G5zeV4i40YIwZqLp7lCoCpqr9x0fc3oo6RmcmDYohvJ2Ji/GzcEoOT71bwqaSajtD3ZgBrLud6O8DIyIZiHH2VNQzKrM/LnTsnTmj0hk5LJ4XNx+kpc2G9RozUHU3gWQBW0TkRRFZ3v6IZGBDUVtQ2VdRz+hO+gYGE58Il0wfyZH6Ft7+sMLrcIwxPdTdJqzvRTII4xysaaS5LcioQZ5AAMYPT2byiBRe217GnFHpJNsl340ZcLo7jPevuDPQA6HXa4F3IhjXkLT3sBuBNXqQN2G1Wzx9BC1tQV7ZesjrUIwxPdDdUVj/CPwJuD80KQ94prcbF5HFIrJdRHaKyJ2dzP8fEXkv9PhARKrC5rWFzRsUzWl7QpdxH5U1+GsgAMNT4pk/JoO1eyopq7XzUo0ZaLrbB3Ib8DGgBkBVdwDDe7NhEfEDvwAuAaYC14vI1PAyqvo1VZ2lqrOAnwFPhc1uaJ+nqlf0JpZosbeijtgYHyNT470Opd9cODmHgN/Hi5utFmLMQNPdBNKkqs3tb0InE/Z2SO98YKeq7gqt+wlOfnb79cDjvdxmVNtTUUdBegI+3+AcwtuZ5LgYzpuYzdbSGnaHmvCMMQNDdxPIX0XkX4EEEVkI/BF4tpfbzgP2h70vDk07gYiMAsYAr4ZNjheRdSKySkSu6mojIrIsVG5deXl5L0OOrD2H6xkzRJqvwp09LovU+BhWvm8nFxozkHQ3gdwJlAObgH8CVgDfiVRQnVgC/ElV28KmjVLVImAp8BMRGdfZgqr6gKoWqWpRdnZ2f8TaI8GgsqeibkgmkNgYHwun5lB8pIHnN5V6HY4xppu6OworiOs0/5KqXquqD/bBWeklQPgt6vJD0zqzhA7NV6paEnreBbwOzO5lPJ4qrWmkqTXI6CGYQABmF6YzIjWe/3phO02tbadewBjjuZMmEHG+JyKHge3A9tDdCL/bB9teC0wQkTEiEotLEieMphKRyUA68HbYtHQRiQu9zsJ18G/pg5g8syfU/j8UayDgTi5cPH0E+yrreXTVPq/DMcZ0w6lqIF/DHZznqWqGqmYAC4CPicjXerNhVW0Fvgy8CGwF/qCqm0XkHhEJH1W1BHiiQ41nCrBORDYArwH3quqATiC7hngCAZgwPJlzxmfx01d3UN3Q4nU4xphTONXpv58FFqrq4fYJqrordD/0vwD/05uNq+oKXH9K+LTvdnj/vU6W+ztwRm+2HW32HK4jPuAjJ2XoDOHtSES485LJfPLnb/HL13fy7UumeB2SMeYkTlUDCYQnj3aqWg4EIhPS0LT7sLsK71AawtuZ6XnD+NSsPH79tz2UVDV4HY4x5iROlUCaezjPnKY9h4fmCKzOfP3iSQD8+MXtHkdijDmZUyWQmSJS08mjlkHWhOSl1rYg+yqH5jkgnclLS+BzHxvDU++WsGF/ldfhGGO6cNIEoqp+VU3t5JGiqtaE1UeKjzTQGtQhO4S3M7ddMI6s5Di+9+xmgkE7udCYaNTdEwlNBO2usBFYHaXEB7hj8STe3VfFM+91dXqQMcZLlkCiwO5ySyCduWZOPjML0viPlds42tTqdTjGmA4sgUSBPRV1pMTFkJkU63UoUcXnE/79immU1zbxs1d3eB2OMaYDSyBRYPfhOkZnJSEytIfwdmZWQRrXzs3n4bd2s7PsqNfhGGPCWAKJArttCO9J3bF4MgkBP3c9vcmu1mtMFLEE4rGm1jZKqhpsBNZJZKfE8e1Lp7B6dyV/Wl/sdTjGmJBTXcrERNi+inpUYawlEAAeW935hRSDqozKSOTu5Zu5aEoOGdZfZIznrAbisfa78FkN5OR8Ilw1O4/GljZ++PxWr8MxxmAJxHPtCWRMpiWQU8lJjefcCdk8+U4xb+044RJtxph+ZgnEY3sq6shIimVYop3Y3x0XTB7O2Owk7nhyI7WNdsl3Y7xkCcRju8rrGJ2Z6HUYA0bA7+NH186ktLqB/1i5zetwjBnSPE0gIrJYRLaLyE4RubOT+TeH7oD4Xuhxa9i8m0RkR+hxU/9G3nc+LD/K+OHJXocxoMwdlc6tHx/LY6v3WVOWMR7yLIGIiB/4BXAJMBW4XkSmdlL096o6K/R4KLRsBnA37u6I84G7RSS9n0LvM0fqmjl8tJkJw1O8DmXAuX3hRGvKMsZjXg7jnQ/sVNVdACLyBHAl3bu3+cXAS6paGVr2JWAx8HiEYo2IneXuzGqrgZye9qG+i6aO4P6/fsiNv1rDPxQVHFdm6YJCL0IzZkjxsgkrD9gf9r44NK2ja0Rko4j8SUTajxLdXRYRWSYi60RkXXl5eV/E3Wd2HLIE0huFGYlcMHk47+6vYkNxldfhGDPkRHsn+rPAaFWdAbwEPHK6K1DVB1S1SFWLsrOz+zzA3thRVktCwE9eWoLXoQxYF0waTmFGIn9+r4Qj9XaTTGP6k5cJpAQIb3fID037iKpWqGpT6O1DwNzuLjsQ7CxzHehD/T7oveH3CdcVFaAKf1y3n6BdK8uYfuNlAlkLTBCRMSISCywBlocXEJGRYW+vANpPQX4RWCQi6aHO80WhaQPKjkNHmWDNV72WkRTLFTNz2VNRzytby7wOx5ghw7NOdFVtFZEv4w78fuBhVd0sIvcA61R1OfAVEbkCaAUqgZtDy1aKyPdxSQjgnvYO9YGitrGFgzWNjLME0idmFaTxYXkdr28vY5SdV2NMv/D0YoqqugJY0WHad8Nefxv4dhfLPgw8HNEAI6j93hZWA+kbIsIVM3M5UNXAH9btZ9m5Y8m1viVjIiraO9EHre0HawGYmGPngPSV2BgfS+cX0hZUbnvsHZpbg16HZMygZgnEI9sO1pIU66cww5pb+lJWShxXz8nn3X1VfO/ZzXYDKmMiyBKIR7aU1jBpRIqNwIqAM/KG8U/nuUud/G7VXq/DMWbQsgTiAVVla2kNk0emeh3KoPWtiyfziSnD+fdnt/Dmjug6gdSYwcISiAcOVDdS29jKFEsgEeP3CT9ZMpvx2cnc9ug7fBi6bIwxpu9YAvHA1gM1AEwdaR3okZQcF8NDNxUR8Pu48VdrKK1u8DokYwYVSyAe2FrqEsikEVYDibSCjEQe+dx8qhtauPFXazhSZ5c7MaavWALxwLaDtRRmJJIc5+lpOEPG9LxhPHhjEXsr67n5N2upa2r1OiRjBgVLIB7YUlrDFGu+6ldnjcvk59fP5v2Sam56eA01dg8RY3rNEkg/q65vYffhOmbkp3kdypCzaNoIfrpkNu/tr+KGB1dTac1ZxvSKJZB+trGkCoCZlkA8cdmMkTx4YxEfHKrl0/e/zaGaRq9DMmbAskb4frZhfxUAZ+QP8zaQQa79roVd+eyZo3h8zT6u/PnfuO+zc5lVkNY/gRkziFgNpJ+9t7+asdlJDEsIeB3KkDY2O5k/fuFsYvzCdfe/zR/X7T/1QsaY41gC6UeqyobiKmu+ihJTc1N59svnUDQqnW/+aSPffmojR22EljHdZk1Y/ehgTSPltU3MtOarqNDezHXJ9JHE+Hw8sWY/Kzcd5Oo5+R/dp37pgkIvQzQmqnlaAxGRxSKyXUR2isidncy/XUS2iMhGEXlFREaFzWsTkfdCj+Udl41G7+2rAmCGtbdHFb9PWDx9BP907lhi/D4e/ttu/rR+P1V2j3VjTsqzGoiI+IFfAAuBYmCtiCxX1S1hxd4FilS1XkS+CPwX8OnQvAZVndWfMffW6t2VxAd8TMu1M9CjUWFmEv984Xhe3VbGWzsPs7G4murGFr503niGJVqflTEdeVkDmQ/sVNVdqtoMPAFcGV5AVV9T1frQ21VAfj/H2KdW765kTmE6cTF+r0MxXQj4fVw8bQS3L5zIGXnDeOCNXZx17yt8b/lm9lbUeR2eMVHFywSSB4QPfSkOTevK54GVYe/jRWSdiKwSkau6WkhEloXKrSsv9+6y3tX1LWw7WMOZYzM9i8F0X3piLP9QVMCKr3ycxdNG8OjqvZz/f17n879Zy8pNpTS1tnkdojGeGxCd6CLyGaAIOC9s8ihVLRGRscCrIrJJVT/suKyqPgA8AFBUVOTZ7enW7KlEFRaMyfAqBNMD7+6romh0BhNzUli1q4I1eyp5ZVsZCQE/M/KHMbswnTsWT0LEbgxmhh4vE0gJUBD2Pj807Tgi8gngLuA8VW1qn66qJaHnXSLyOjAbOCGBRIvVuyqIjfEx0zrQB6TUhACLpo3gE1Nz2Fl2lHf3HeGdfUdYvbuSP6zbz8z8YczITyMnNf6EZW0klxmsvEwga4EJIjIGlziWAEvDC4jIbOB+YLGqloVNTwfqVbVJRLKAj+E62KPWqt0VzClMIz5g/R8DmU+EiTkpTMxJobGljc0Hqtmwv5rXt5fz2vZyRqTGf5RM0pNivQ7XmIjyLIGoaquIfBl4EfADD6vqZhG5B1inqsuBHwHJwB9DTQT7VPUKYApwv4gEcf0493YYvRVVymobeb+khq8vnOh1KKYPxQf8zB2VwdxRGdQ2trCppJqNxdW8uOUQL245RGFGIjPyh7Fwag7ZKXFeh2tMn/O0D0RVVwArOkz7btjrT3Sx3N+BMyIbXd95bZurPF00JcfjSEykpMQHOHtcFmePy+JIXTMbi6vYUFzNcxtLWbGplLPHZXHFzFwunj7CLmNjBo0B0Yk+0L2ytYzcYfF2D5AhIj0plvMmDee8ScM5VNNIW1BZvuEA33pyI9955n0WTs3hhgWFnDUu0zrfzYBmCSTCGlvaeHPHYa6Zm2cHiyEoJzWepQsK+fqiiWworuaZd0t4+t0Snt9UytisJJYuKOSaOfnWX2IGJEsgEfb2rgoaWtqs+WoIC7+0/MScFG5fOJH3S6pZvbuSHzy/lXtXbuPT8wpYdu5YRmUmeRipMafHEkiEPfveAVLiYjjLTiA0IQG/j9mF6cwuTOdgdSNv7zrMH9cV8/iafVw2I5cvnDeWabl2wU0T/exy7hF0tKmVle8f5PKZuTZ813RqxLB4PjU7n7fuuIB/PHcsr20r47KfvsXNv17D+r1HvA7PmJOyBBJBKzaV0tDSxrVzT3aFFmNgeGo8375kCn+780K+efEkNhZXc83//p0bHlrF6l0VXodnTKesCSuCnlxfzJisJOYUpnsdioly4f0k6YmxfOXCCazeXcGbOw7z6QdWMSYriR9cNZ2zbeSWiSJWA4mQ9k7S64oK7B/enLbYGB8fn5DNNy+exOUzRlJxtIkbHlrNNf/7d17fXoaqZ5d1M+YjVgOJkP99/UNS4mK44Uy7DpLpuYDfx9njspg3OgOfT/jf13Zy86/XMn54MkvmFXD1nHwybAiw8YglkAj4sPwoK94v5UvnjyM13s46Nr0X8LvGgi+cP44N+6tZs7uCHzy/lf9YuY3x2clMHpnCpJwUvnTBeI8jNUOJJZA+pqr8x4ptxMf4ueVjY7wOxwwyMT4fc0elM3eUGwK8fm8lW0pr2H6oFoDH1+7jjLxhTMsdxvjhyeSlJZCfnsCwhMBHTanh/S1dsSsIm+6wBNLHXtx8iJe3HuJfL51MVrJdQM9Ezohh8Vw2I5dLzxhJeW0T2w/V4hPh/QPVrNh08LiySbF+ctMSyE1LoL65lWEJsaQlBBiWGCAtIUBaYix+n/XVmdNjCaQPldU2cvfy95kyMtVqH6bfiAjDU+MZHrpsCkB1Qwv7Kuopqaqn+EgDJVUNlBxpoLS6kV2H66hraj1uHX6fkJkUy/CUOLJT4kmK8zN+eDLjspPtHCbTJUsgfaSxpY1lv11PTUMrv7pp3kdt1sb0p86apxJjY5gwPIUJw49dzLOlLUhNQwtVDS1U1TdTXttMeW0jpdWNbD5Qw2vb3RWkRSAvLYHRmUmMzkp0z5lJjEyLJyc1nozEWHxWcxmyLIH0gZrGFm579B3e21/FfZ+Zy/Q8uwyFiW4Bv4/M5DgyO2lmbWkLsmBsBjvLjrKz7Ci7yuvYW1HH8vcOUNN4fM0lxidkp8R9VHMZnhpHdnIcWSlxbCutISUuhuT4AMlxMcTGnPijyvpaBjZPE4iILAb+L+6GUg+p6r0d5scBvwXmAhXAp1V1T2jet4HPA23AV1T1xX4M/SPr9x7hjic3sudwHT+6dgaLp4/o/Uo3bYKd5ZCcBEePgj8G4uJCr/0QFwutbdDUCAkJEBsLdfXQ1uZexwagvh5iAq58Wyu0trrXMTHQ0OjKBWLcMo1NEB8PzU3g87mHiCvTvk2AQCz4BOrq3OvYWDhaC7FxkJICTU1QdxR8fhd7fX1oPXGgQWhocNtpbHRlEuLdftTXQUoqtLZASwskJrpt+Hxu2YZ6t59tbRAMQmoqVFUdW3dDg9teczOounL+GDc/xu/iSElx+xIX5+JMSnbPLc1uXkODWz4tzU1vbXWfleLKJCa6zzg29Bk0N7vPxR/jPsfaWveZoC7uYNAt39zipsXHu32KiwPxuWlZWXDkiFs2Lj60P4FjywSD7nNpbnLP1dXHPivU7VdGJmRkwMGD7rtIT4fMLPd9tbbA4QoXrwhUVLj9aP8O/DHub+noUfc9ZGRAdTUBESbPGcHkEbnH/VmqKlUHK9i94QMO1bdyKDGNspJyylqEQ34/xUfqeXffESrrm+nsNJU4H6TFCukJAdLSkslIiiUjKUB+eiIF6YmkJsS4jv7Dh6G01O1Lfn7v/5+qq2H/fkhKgjGdNC03NsLu3e77GjPGPUdK+7Z8Phg7FgIDe5SmZwlERPzAL4CFQDGwVkSWd7iz4OeBI6o6XkSWAP8JfFpEpuJugTsNyAVeFpGJqtrWH7G3tAVZtauCR1ft44XNBxmRGs9vPzefs8dn9X7lzz0Ht90Gd/4EnnsWtm51B51Fi2DbdnfAuWQxPPNnqK2BkSPhyith3Tr3GD0GzjvXLdPWBofK4EglrFoFQYUFC9xBcPcuuPQyePllmDUT/vLSsfWddx7k5sIf/wg7drgDzsKFMGIEPPF7Vy43Fy64AH7/B0hOhuv+Adaug3ffcQfDhQthymTYuRPy8uHZ5VBQAE3NsH6dSyCfuwVWvgDDs90/99/+5mIsmgvJKfDX12H+fHege/dduOxSdyDc/oFb7+WXwV/+4vY5bRgMHw5vvQUlJS6mhQvdP2zhKPjtb6G83CWIT10NmzfDK6+4RDB1KkycCM88A/kF8MnLXVwzZ8BLL7uDdmGhi6W+HkoPwnvvurguXuRiev559/4TF7mEXHoARo+G115z25g8BcaNhdVrXNyBUJJ/6mnYu8clyCuvcq8LCuDNN+Gqq+DXv3bru/RSeGGlSz6LFsLKlVBTAzkj4NPXwaOPugSRMwJuWArDc+D559x3cuUV8Oqrbv/nz3ef04svumQ9c6b7Luvq3A+LdetAhMfePtftS/bwY3+bB0pg9Wp45VVoaSZ22jTyC0eR//prcNllMHs2ZI2mLajU7S/h6J+fo/ZQBUeTh3F01lxqahuoKq2lKjGV3UcyaMLHivePdfanxMWQGwe5O98nt2QXeW315C29mhHnnkl6SjxpiQHSEmI7rcl0acsW+MIX4M03eWzBVe6zmjcvlOyBygr483JYvYqlW16Fb30LvvIVyM7u1b9xp/btgx/8AH71K5dAvvY1uP129381QIlXZ7SKyFnA91T14tD7bwOo6n+ElXkxVOZtEYkBDgLZwJ3hZcPLnWybRUVFum7dutOO9e0PK9h8oJpdh+vYXV7HppJqjja1kpYY4IYFhXzx/PEkx/VBLl63zh30pk3jsemfcAf9cFdd5X4prlgBwbBcOTwHrrgCHnrQvZ88BeYVwZ49MGyYS0rhFi92v55ffx0uvBCe77C+nBw45xx48snjl/vUp+Dpp4/f7tixsOptl7AuvNAd2NotXeoO3n/4vftFee657sAFMGeO+6VZUgKXXOqSZbgLL3QHv9oauPhi+OsbrsZ1443wyCPHpjU3u8RSetAdtEtKjq1DfHDzzfD0U67G0tV+ABQVwcFDULzfJZGLF8HDD3PcT+mx41xy+8Mfjl/2yith+bOultX+HoE/P3N8uXnz4MABqKmF274ETz0FH3xwfJnPf94lo1Gj3Y+HuXNccnnpZVcTuvxy9xm3hX1fWdkwaaJLwOBqNtdcC/ff55L8+negptoluiuuPDGuCy5w398LLxw//TOfgY+f6143N8Pbf4fHHju+zIIzYd9e9/3e+o9uH5ub4cEHYeOG48tedZVL0riKXeOXv8KRgnFU1jVTVd/Mkeo6qjZto1oCVCWkUB+bQGeSYv0kxsUQ6/cRG+Mj4BcCfh9BhWBQaQ0GCSq0tbXRdqictqYm2sRHfSA+tIJExB+DANLcBA0N+DRIWsNRYoKtxBTkE5uZToxPiPG79cf4fATaX/t9BHxumzGhbcf4hECMm76ltOajfUTdsyqwZze6ZQuKMLVsN0ER+OQniZk+Fb/PrcPvEwJ+Oe79R8+hOMLfd2e5gN+HP9TU2NO+WRFZr6pFHad72YSVB+wPe18MLOiqTOge6tVAZmj6qg7LRuyKhQ+9uYtXtpWRnhhgTFYSV8zK5dwJWZw/aXjfjlDZt88d6C67DJavOXF+Q4P75RrsUNEqO+T+aZOTXXPE9m2uFlFa6g5WHW3bDud+3P3qbGs7cX2HDrlf7h01NZ243aK5x+Z1vGTLkSMwMhd27YIJE+HDD4/NGzEC3nnH/WI+UMIJduyEUaPg/U2wazcU5LvaUGWlm+/3Q2OD+6VYetDVYtatPX4dGnTNX+HJIybgPqOOtmx1Sbd4v3uIcEI7zK4PXW2to/Jyl6irjhzb77hOhnBv3eoS5xtvuFpFx+QB7vs6eNCVq6l2yQNc8gDXrNXW4fs6XA7z54W9PwwVh93rhAS3HnB/H5WdXJjxaJ37sdHR5i3HEkhtjft76mjbVpg+3c0rO+Sa22pqYdPGE8s2N3/0UoCE/ftIOGM6uWmhfdzfAA8+day4P4aq+BSKll3PkVlFrsO/rpkj9S00tLTS3Ko0twVpaQ3S0hZE5NhB0+cT/HV1+N9Yg1/b8AeDfJgROkRMmoTmZEBbG7phJ1pfR1B8FFQfotUXQ0tqgJbCHFrblJa2II0tQVraWmlpU1rb3LZa2lyiam1zMbSG3re0df2DXDQRGTMXUN4eNRNQKPMT/OsuWoOR/yH/8u3nMX54cp+uc9B3oovIMmAZQGFhzzrs7rlqOj+O9ZOWGOFLRmRluQPjtm0s9ZfDhr8dP/+aj7nnDR26e9LSIPdmWLPcNUtMmgRJNdCyHzLHnVh+7ufAXw5bXnXr7Dg/PR2Szj1xeseyaWnQONVNE4Grzz5+/tILYXwSVG6BtdvgoouOzT97LJRthD2r4NZbT9zWlCXwzlvuoFuUDy+95BLs5y9zZa86Eza95Nr0Z38RanfAgXfcwTzc166HD95wyRdcnFfOP3F7Z58NO95003NyYHjwxDIFBZC08MTpl86B1c8cS7qXF0Fs64nlzjkHtv8Vtv0VJvybi3nXruPLfPXTsH8tNExxZS+b65oQN/7FJbTOvq/UVGg549j01FQoXOreXzbXracx1O91/uQTl5+TC5k+eKbD9KvPhvZO7ookeK/txGUvuAC2vQ4bXoUvfsp9r1VV0LTXfXfhrj3n+OXvuOnY+gFygrB39fEJH+DHX4WZx/fHdEtFBfx4JWzbdvz0lSth8VkuGX/xd/DoA8fP/8Uv4Jb5p789XD9RezLwibhajrih1tx+O/zP/xy/wI9/DLffjqoSVNc83hZ062gL1abaghpKUEpbMEhr6H17udYTljm+XPv74al9f16al2NNS4CCsPf5oWmdlgk1YQ3DdaZ3Z1kAVPUBVS1S1aLsHrZr5qUlRD55gGtG+f73XXv2t77l+gXaffzjsHeva8v+3OeOTff54DvfcR3vLS3ul+9tt7kksHSpO5iMD7u8xciR7v1998Hdd7v13Xjj8ev71rdgxgx34Gq3aJFbZ3i5b37zWJPGXXcda54CVwOaPds1KX3nO67pZswYyAv9Cnz8cbjjjlCHfMAlvXY5OTBrljsA5ee7voR9++D6612STUtzzT9f/aqrTYjA+vXwjW8c6/AHWLIE1qyB//zPY7UjVVf7OeecY+WSk+G661z8fj/827/B8uWuyaVdIOC2FxfnOpvbzZ3r1tmePObMca9LSuDMM4+VS0mBa691fU7f+pY7eP3wh+6g3u7qq+H9993n+uijcM897nN69ln44hddmZdeOvH7/+EPXXlw+/mjH7mmpNxct/w3vuHmNTe7fpN5YbWVtDQ44wwYN879bbQbN871u7TLzHR/g/PDDqypqa7p9LXXXNl589znl5npDsLhtbArrzw+WZ5zjkva4UaPhvvvd/vU7pvfhGnT6JHMTPjlL4+P4x/+wf1dgtvOl77k+oTanXGGa0buIREJNXX5PqoJfXQx1Vtucd9Ju8mTXWtDaDm/T4gP+EmKi2FYQoCMpFiGp8QzclgCBRmJjMlKYvzwFCaPSGV63jBmFqQxd1Q6C8Zmcvb4LM6dmM0Fk4ezcGoOi6eP5PIZuVw1O49r5+bz6XmFEbmskpd9IDHAB8BFuIP/WmCpqm4OK3MbcIaqfiHUiX61ql4nItOAx4D5uE70V4AJp+pE72kfSL8qKXEdvDU17kC6d69LJDk5rskjPt4dVI8ccU1NeXnuH2XfPrfMyJHuoHDwoDto+f3ugFZR4V7n5bl1DhvmytXVHRsddOSImx4T45JFS4tbb1ycawpJS3MH7NJSdxBOTHQHhYwM9+u8rs41U8XHu+anlhbX5FRQcCze3Fy3jzEx7iBVWenejxvnYqyvd/0qe/a4g2F+vhu1khm6o2Mw6P7hP/jAxRoXB2VlbpmqKncwr6x0B7dAwD0fOuRiP3TI7Vd9vVvfkSNuf8aMcSN1Dh50cQQCUFzslm1tPfa5lpW5bSYkuFE9CQmuCS0mxn0OiYluflub25eRI902Ghrc515S4uJISXHNTJMnu23u3u0+w2Gh4d9VVW6/YkIjyerrj63rwAEXY12di2f8ePfYsMF9L2PHugN5UpL7DN9/38Wp6mLMzXWfQUmJm5af7z6Xujq3ng8/dNudOROmTDn+b1MVNm50v+hbW90Bf+9eF+u0aW5/wstu2XLsbzYtzY02Ky11n/28ee5vuqPmZrfcrl1u/hlnuO+hp8LjSE9362v/W2q3a5crEwi45rj2HzmR0P6d+P1uWwUFp1wkGnTVB+JZAgEQkUuBn+CG8T6sqj8UkXuAdaq6XETigd8Bs4FKYImq7gotexfwOaAV+BdVXdnZNsINiARijDFRJioTSH+zBGKMMaevqwRi19swxhjTI5ZAjDHG9IglEGOMMT1iCcQYY0yPWAIxxhjTI5ZAjDHG9IglEGOMMT1iCcQYY0yPWAIxxhjTI5ZAjDHG9IglEGOMMT1iCcQYY0yPWAIxxhjTI5ZAjDHG9IglEGOMMT1iCcQYY0yPeJJARCRDRF4SkR2h5/ROyswSkbdFZLOIbBSRT4fN+42I7BaR90KPWf26A8YYYzyrgdwJvKKqE3D3M7+zkzL1wI2qOg1YDPxERNLC5n9TVWeFHu9FOmBjjDHH8yqBXAk8Enr9CHBVxwKq+oGq7gi9PgCUAdn9FaAxxpiTi/FouzmqWhp6fRDIOVlhEZkPxAIfhk3+oYh8l1ANRlWbulh2GbAs9PaoiGzvVeT9Jws47HUQEWD7NbDYfg0ckdynUZ1NFFWNyNZE5GVgRCez7gIeUdW0sLJHVPWEfpDQvJHA68BNqroqbNpBXFJ5APhQVe/p0x3wmIis6+wm9gOd7dfAYvs1cHixTxGrgajqJ7qaJyKHRGSkqpaGkkFZF+VSgeeBu9qTR2jd7bWXJhH5NfCNPgzdGGNMN3jVB7IcuCn0+ibgzx0LiEgs8DTwW1X9U4d5I0PPgus/eT+SwRpjjDmRVwnkXmChiOwAPhF6j4gUichDoTLXAecCN3cyXPdREdkEbMK1+/2gX6PvHw94HUCE2H4NLLZfA0e/71PE+kCMMcYMbnYmujHGmB6xBGKMMaZHLIF4TEQeFpEyEel0IIA4PxWRnaFLuszp7xhPVzf26XwRqQ7r2/puf8fYEyJSICKviciW0CV2vtpJmQH1fXVznwbc9yUi8SKyRkQ2hPbr3zspEycivw99V6tFZLQHoZ6Wbu7XzSJSHvZ93RqxgFTVHh4+cAMF5gDvdzH/UmAlIMCZwGqvY+6DfTofeM7rOHuwXyOBOaHXKcAHwNSB/H11c58G3PcV+vyTQ68DwGrgzA5lvgTcF3q9BPi913H30X7dDPy8P+KxGojHVPUNoPIkRa7EDWVWdefCpLUPY45W3dinAUlVS1X1ndDrWmArkNeh2ID6vrq5TwNO6PM/GnobCD06jhgKv6TSn4CLQqcGRK1u7le/sQQS/fKA/WHvixkE/+DAWaFq+EoRmeZ1MKcr1NwxG/cLMNyA/b5Osk8wAL8vEfGLyHu4E5VfUtUuvytVbQWqgcx+DbIHurFfANeEmlD/JCIFkYrFEojxwjvAKFWdCfwMeMbbcE6PiCQDTwL/oqo1XsfTF06xTwPy+1LVNlWdBeQD80Vkusch9Ylu7NezwGhVnQG8xLFaVp+zBBL9SoDwXxD5oWkDlqrWtFfDVXUFEBCRLI/D6hYRCeAOtI+q6lOdFBlw39ep9mkgf18AqloFvIa7LUS4j74rEYkBhgEV/RpcL3S1X6paoccuLvsQMDdSMVgCiX7LgRtDo3vOBKr12LXABiQRGdHe1izuSss+BsA/bijmXwFbVfW/uyg2oL6v7uzTQPy+RCRbQvcPEpEEYCGwrUOx8EsqXQu8qqFe6GjVnf3q0Od2Ba5fKyK8upy7CRGRx3GjXLJEpBi4G9cxhqreB6zAjezZibvJ1i3eRNp93dina4Evikgr0AAsifZ/3JCPAZ8FNoXaoAH+FSiEAft9dWefBuL3NRJ4RET8uIT3B1V9TkTuAdap6nJc4vydiOzEDfpY4l243dad/fqKiFwBtOL26+ZIBWOXMjHGGNMj1oRljDGmRyyBGGOM6RFLIMYYY3rEEogxxpgesQRijDGmRyyBGGOM6RFLIMYYY3rk/wPY+YEbutqAowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "zero = [0 for i in range(len(X))]\n",
    "df = pd.DataFrame(zero,X)\n",
    "fig2 = sns.scatterplot(data = df,marker = '^',palette = ['r'])\n",
    "fig2 = sns.distplot(X)\n",
    "plt.ylim(-.25,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5301920-c96e-41dd-a7c6-b9cbd9d5ca83",
   "metadata": {},
   "source": [
    "# Problem 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "id": "0a0dad17-7cc9-4e6d-bfd3-be8b7bb7cca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 0.01\n",
    "theta = 0.08\n",
    "k = 2\n",
    "Beta = 0.2\n",
    "rho = -0.7\n",
    "Y0 = 0.04\n",
    "S0 = 50 \n",
    "A = 80\n",
    "B = 30\n",
    "K = 50\n",
    "n = 10000\n",
    "dt = (1/3)/360\n",
    "days = 90/360"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "id": "1c35afb1-f66f-4dc6-bebe-37cab67ad43b",
   "metadata": {},
   "outputs": [],
   "source": [
    "Pairs = list()\n",
    "Ypaths = list()\n",
    "Spaths = list()\n",
    "\n",
    "for i in range(n):\n",
    "    Y = list()\n",
    "    S = list()\n",
    "    Y.append(Y0)\n",
    "    S.append(S0)\n",
    "    \n",
    "    for _ in range(int(days/dt)):\n",
    "        rand1 = np.random.normal(0,1)\n",
    "        Ynext = Y[-1] + k*(theta/k-Y[-1])*dt + Beta*np.sqrt(Y[-1])*(rho*np.sqrt(dt)*rand1 + np.sqrt(1-rho**2)*np.sqrt(dt)*np.random.normal(0,1))\n",
    "        Snext = S[-1] + r*S[-1]*dt + np.sqrt(Y[-1])*S[-1]*np.sqrt(dt)*rand1\n",
    "        \n",
    "        if Ynext < 0:\n",
    "            Ynext = 0\n",
    "            Y.append(Ynext)\n",
    "        else:\n",
    "            Y.append(Ynext)\n",
    "            \n",
    "        S.append(Snext)\n",
    "    \n",
    "    Ypaths.append(Y)\n",
    "    Spaths.append(S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "id": "149673fe-e870-4730-acd0-d60ac481e17f",
   "metadata": {},
   "outputs": [],
   "source": [
    "terminalY = [Path[-1] for Path in Ypaths]\n",
    "terminalS = [Path[-1] for Path in Spaths]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "id": "e2e18bc7-03bf-4428-8efc-89203bd48562",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(50.10907720968314, 0.040115502019253305)"
      ]
     },
     "execution_count": 288,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expectedY = sum(terminalY)/len(terminalY)\n",
    "expectedS = sum(terminalS)/len(terminalS)\n",
    "expectedS,expectedY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "id": "edb1cedf-3d71-42fb-bc70-903a4d050072",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7fa4f33d48b0>"
      ]
     },
     "execution_count": 289,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYRUlEQVR4nO3df5Bd5X3f8ffXElooYEngDSgLGcmzSlubynZmjZ2sp25gnGDyA5o62E0aKy6BP8AdPHYccNyZOpl2ijNpiFoVZxhILTxJgFAbSMLgypjYtGNjX7B8bTCO1yDBbgAJkIQNowXBt3/cg7h7taDV7p773B/v18zOPec551x9987uR2ef85znRGYiSeq+15UuQJKGlQEsSYUYwJJUiAEsSYUYwJJUyMrSBSzFOeeck3fccUfpMiTpSGK+xr4+A37yySdLlyBJi9bXASxJ/cwAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKqSvp6NU983OztJoNOa0TUxMMDIyUqgiqX8ZwDoqjUaDy66+hdVj4wDsn5liyyUwOTlZuDKp/xjAOmqrx8YZHd8EwEsHX6DZbB62j2fF0pEZwDqks3vh+eefB2DVqlWH2prNJpmvTO7/zOO72LrzAKc+/MrlBM+KpYWpNYAjYifwI+BF4GBmTkTEScCNwHpgJ3BBZu6NiAC2AOcCzwG/nZn31Vmf5ursXpjZ8RVWnHAyp46fcWifmR13s2bjxJzjTly34dAZsaSF68YZ8M9nZvuzg64A7szMKyPiimr9cuC9wMbq6x3AZ6pXdVF798K+6SmOWbtuTrjum54qVZo0cEoMQzsP2FYtbwPOb2u/Plu+DqyJiHUF6pOkrqg7gBP4PxFxb0RcXLWdkpmPVcuPA6dUy2PAo23HTldtc0TExRHRiIjGnj176qpbkmpXdxfEuzJzJiJ+AtgeEQ+2b8zMjIg8mjfMzGuAawAmJiaO6lhJ6iW1ngFn5kz1uhv4AnAm8MTLXQvV6+5q9xng9LbDT6vaJGkg1XYGHBHHA6/LzB9Vy78A/CFwG7AZuLJ6vbU65DbgwxFxA62Lb/vbuiq0zOa7o61ziJmketXZBXEK8IXW6DJWAn+ZmXdExDeBmyLiQmAXcEG1/+20hqBN0RqG9qEaaxt6nUPOYP4hZpLqU1sAZ+ZDwFvmaX8KOHue9gQuraseHa59yBk4xEzqNmdDk6RCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsT5gIdE540X3nQhlWcAD4nD5/r1pgupNAN4iHTO9SupLPuAJakQA1iSCjGAJakQ+4C17OZ7VL2PqZcOZwBr2XU+qt7H1EvzM4BVCx9VLx2ZfcCSVIhnwAPIxw1J/cEAHkA+bkjqDwbwgPJxQ1Lvsw9YkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgpxLog+NN9sZz5xQuo/BnAf6pztbO8jD3LRu5ts2tSafMepJ6X+YAD3qfbZzvZNT7F1+wOHHgHk1JNSfzCAB0T7I4CcelLqD16Ek6RCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCnA1NtXvp4As0m83D2p1EXsPOAFbtnnl8F1t3Hjg0XzHA/pkptlwCk5OTBSuTyjKA1RXt8xVLarEPWJIKqT2AI2JFRHwrIv62Wt8QEfdExFRE3BgRq6r2kWp9qtq+vu7aJKmkbpwBXwZ8r23908BVmTkO7AUurNovBPZW7VdV+0nSwKo1gCPiNOCXgGur9QDOAm6udtkGnF8tn1etU20/u9pfkgZS3WfAfwr8HvBStX4ysC8zD1br08BYtTwGPApQbd9f7T9HRFwcEY2IaOzZs6fG0iWpXrUFcET8MrA7M+9dzvfNzGsycyIzJ0ZHR5fzrSWpq+ochjYJ/GpEnAscC7we2AKsiYiV1VnuacBMtf8McDowHRErgdXAUzXWJ0lF1XYGnJmfyMzTMnM98AHgy5n5m8BdwPuq3TYDt1bLt1XrVNu/nJlZV32SVFqJccCXAx+NiClafbzXVe3XASdX7R8FrihQmyR1TVfuhMvMvwf+vlp+CDhznn0OAL/ejXokqRd4J5wkFWIAS1IhBrAkFWIAS1IhBrAkFWIAS1IhBrAkFWIAS1IhBrAkFWIAS1IhPpSzx83OztJoNOa0NZtNMp2rXup3BnCPazQaXHb1LaweGz/UNrPjbtZsnChYlaTlYAD3gdVj43Me6b5veqpgNZKWiwGsIl46+ALNZnNO28TEBCMjI4UqkrrPAFYRzzy+i607D3Dqw63rwPtnpthyCUxOThauTOoeA1jFnLhuw5yuFWnYOAxNkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpEANYkgoxgCWpkJWlC5AAXjr4As1mc07bxMQEIyMjhSqS6mcAqyc88/gutu48wKkPt/4o2/vIg1z07iabNm2as5+hrEFiAKtnnLhuA6PjrcDdNz3F1u0PHApkgP0zU2y5BCYnJ0uVKC0rA1g9qz2QpUHkRThJKsQAlqRC7ILoMbOzszQajUPrzWaTzChYkaS6GMA9ptFocNnVt7B6bByAmR13s2bjROGqJNXBAO5Bq8fG54wGkDSYausDjohjI+IbEfHtiLg/Iv6gat8QEfdExFRE3BgRq6r2kWp9qtq+vq7aJKkX1HkRbhY4KzPfArwVOCci3gl8GrgqM8eBvcCF1f4XAnur9quq/SRpYNUWwNny42r1mOorgbOAm6v2bcD51fJ51TrV9rMjwqtPkgZWrcPQImJFROwAdgPbgR8C+zLzYLXLNDBWLY8BjwJU2/cDJ8/znhdHRCMiGnv27KmzfEmqVa0BnJkvZuZbgdOAM4F/tgzveU1mTmTmxOjo6FLfTpKK6cooiMzcFxF3AT8LrImIldVZ7mnATLXbDHA6MB0RK4HVwFPdqE/9wRnTNGhqC+CIGAVeqML3OOA9tC6s3QW8D7gB2AzcWh1yW7X+tWr7lzMz66pP/adzxjQn51G/q/MMeB2wLSJW0OrquCkz/zYiHgBuiIj/DHwLuK7a/zrgcxExBTwNfKDG2tSnnKBHg6S2AM7MJvC2edofotUf3Nl+APj1uuqRpF7jZDySVMiCAjgiDutkm69NkrRwCz0D/h8LbJMkLdBr9gFHxM8CPweMRsRH2za9HlhRZ2GSNOiOdBFuFXBCtd+Jbe3P0BoqJklapNcM4Mz8CvCViPhsZu7qUk2SNBQWOgxtJCKuAda3H5OZZ9VRlCQNg4UG8F8DfwZcC7xYXzmSNDwWGsAHM/MztVYiSUNmocPQ/iYiLomIdRFx0stftVYmSQNuoWfAm6vXj7e1JfDG5S1HkobHggI4MzfUXYgkDZsFBXBEfHC+9sy8fnnLkaThsdAuiLe3LR8LnA3cBxjAkrRIC+2C+A/t6xGxhtaE6pKkRVrsdJTPAvYLS9ISLLQP+G9ojXqA1iQ8/xy4qa6iJGkYLLQP+I/blg8CuzJzuoZ6JGloLKgLopqU50FaM6KtBZ6vsyhJGgYLfSLGBcA3aD2z7QLgnohwOkpJWoKFdkF8Enh7Zu6GQ4+c/xJwc12FSdKgW+goiNe9HL6Vp47iWEnSPBZ6BnxHRHwR+Ktq/f3A7fWUNDxmZ2dpNBpz2prNJplRqCJJ3XSkZ8KNA6dk5scj4teAd1Wbvgb8Rd3FDbpGo8FlV9/C6rHxQ20zO+5mzcaJglVJ6pYjnQH/KfAJgMz8PPB5gIj4F9W2X6mxtqGwemyc0fFNh9b3TU8VrEZSNx2pH/eUzPxOZ2PVtr6WiiRpSBwpgNe8xrbjlrEOSRo6RwrgRkRc1NkYEb8D3FtPSZI0HI7UB/wR4AsR8Zu8ErgTwCrgX9dYl3RELx18gWazeVj7xMQEIyMjBSqSjs5rBnBmPgH8XET8PHBG1fx3mfnl2isbQJ3DzhxytjTPPL6LrTsPcOrDr/wht39mii2XwOTkZMHKpIVZ6HzAdwF31VzLwOscduaQs6U7cd2GOaNIpH6y0BsxtEzah5055Ewabt5OLEmFGMCSVIgBLEmFGMCSVIgBLEmFOApCA2W+mzO8MUO9ygDWQOm8OcMbM9TLDGANHG/OUL+wD1iSCjGAJakQA1iSCjGAJakQA1iSCjGAJakQh6FpoPnUDPUyA1gDzadmqJcZwBp43pihXmUfsCQVYgBLUiG1BXBEnB4Rd0XEAxFxf0RcVrWfFBHbI+IH1evaqj0i4r9HxFRENCPiZ+qqTZJ6QZ1nwAeBj2Xmm4B3ApdGxJuAK4A7M3MjcGe1DvBeYGP1dTHwmRprk6TiagvgzHwsM++rln8EfA8YA84DtlW7bQPOr5bPA67Plq8DayJiXV31SVJpXekDjoj1wNuAe4BTMvOxatPjwCnV8hjwaNth01Vb53tdHBGNiGjs2bOnvqIlqWa1B3BEnAD8b+AjmflM+7bMTCCP5v0y85rMnMjMidHR0WWsVJK6q9YAjohjaIXvX2Tm56vmJ17uWqhed1ftM8DpbYefVrVJ0kCqcxREANcB38vMP2nbdBuwuVreDNza1v7BajTEO4H9bV0VkjRw6rwTbhL4LeA7EbGjavt94Ergpoi4ENgFXFBtux04F5gCngM+VGNtklRcbQGcmf8XiFfZfPY8+ydwaV31SFKv8U44SSrEAJakQgxgSSrEAJakQpwPWENnvqdk+IQMlWAAa+h0PiVj7yMPctG7m2zaNHfSdkNZdTOANZTan5Kxb3qKrdsf8LFF6joDWMLHFqkML8JJUiEGsCQVYgBLUiEGsCQVYgBLUiEGsCQVYgBLUiEGsCQVYgBLUiEGsCQVYgBLUiEGsCQVYgBLUiEGsCQVYgBLUiHOB1yT2dlZGo3GnLZms0lmFKpIUq8xgGvSaDS47OpbWD02fqhtZsfdrNk4UbAqSb3EAK7R6rHxOU9Z2Dc9VbAaSb3GPmBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCDGBJKsQAlqRCvBFDmsdLB1+g2WzOaZuYmGBkZKRQRRpEBrA0j2ce38XWnQc49eHWH4n7Z6bYcglMTk4WrkyDxACWXsWJ6zbMuZVcWm72AUtSIQawJBViAEtSIQawJBViAEtSIQawJBViAEtSIY4DlhbAO+NUBwNYWgDvjFMdDGBpgdrvjJvvjBg8K9bRMYClReg8IwbPinX0DGBpkZwrQkvlKAhJKsQAlqRCauuCiIg/B34Z2J2ZZ1RtJwE3AuuBncAFmbk3IgLYApwLPAf8dmbeV1dtdZidnaXRaBxabzabZEbBiiT1ujr7gD8LbAWub2u7ArgzM6+MiCuq9cuB9wIbq693AJ+pXvtGo9HgsqtvYfXYOAAzO+5mzcaJwlVJ6mW1dUFk5leBpzuazwO2VcvbgPPb2q/Plq8DayJiXV211WX12Dij45sYHd/E8W8YK12OpB7X7T7gUzLzsWr5ceCUankMeLRtv+mq7TARcXFENCKisWfPnvoqlaSaFbsIl5kJ5CKOuyYzJzJzYnR0tIbKJKk7uh3AT7zctVC97q7aZ4DT2/Y7rWqTpIHV7QC+DdhcLW8Gbm1r/2C0vBPY39ZVIUkDqc5haH8F/CvgDRExDfwn4Ergpoi4ENgFXFDtfjutIWhTtIahfaiuuiSpV9QWwJn5b19l09nz7JvApXXVInWDU1bqaDkXhLRMnLJSR8sAlpaRU1bqaBjAUk2cslJHYgBLNXLKSr0WZ0OTpEIMYEkqxACWpEIMYEkqxACWpEIMYEkqxACWpEIMYEkqxACWpEK8E24ROp+ADD4FWdLRM4AXofMJyOBTkLUwTlmpdgbwIr38BOSX7ZueKliN+oVTVqqdASx1mRP06GVehJOkQgxgSSrEAJakQgxgSSrEAJakQgxgSSrEAJakQgxgSSrEAJakQgxgSSrEW5GlHjPfbHtO2DOYDGCpoPlmR2s2m1z71YdYc3prtj0n7BlcBvACdJ6ROPevlkvn7GjwytSmTtgz+AzgBeic/9e5f7WcOmdHc2rT4WEAL1D7/L/+gkhaDgaw1OPm6ycGL8wNAgNY6nHz9RPvfeRBLnp3k02bXum6MJD7jwEs9YH5+om3bn/ARxv1OQNY6lPtoWw3RX8ygKUBMF83hWfFvc8AlgaED/vsP84FIUmFGMCSVIhdEPPw1mMNgvkuzHlRrrcYwPPw1mMNgs4Lc16U6z0G8Kvw1mMNAi/M9TYDWBpS8807DHZTdJMBLA2pzq42OPwW5+effx6AVatWzTnWkF4eQx/A850FeNFNw6K9qw0Ov8V5ZsdXWHHCyZw6fsahfexLXj5DH8DznQV40U3DrL3feN/0FMesXWc/ck2GPoBh/rMASaqbASwNic5xwXa1lTd0AexNFhpWneOC7Worb+gC2JssNMw6+3dVVk8FcEScA2wBVgDXZuaVdfw73mQhLd5CbnGeb3SRQ9cO1zMBHBErgP8JvAeYBr4ZEbdl5gNlK5PUrrMrY77HIzWbTa796kOsOX38VffpHGM835jj5Qj2xd5w0o3/RHomgIEzganMfAggIm4AzgOWPYD3z7xy1vvskzOsOHCAPccf/6ptR1rvtX16vb5B+B56vb7av4cTTj60/tzTT/BfPzfF2nXfPdT25EPfZc2GTUfcZ8Vxr2ftup+ad/3Zpx7jY+9/z2HB/t9u3M7xJ6971X06dR6z2OOefeoxtv2Xjy3r+OfIzGV7s6WIiPcB52Tm71TrvwW8IzM/3LHfxcDF1eo/Bb7/Km/5BuDJmsqtUz/W3Y81g3V3Wz/WvVw1P5mZ53Q29tIZ8IJk5jXANUfaLyIamdl3V9f6se5+rBmsu9v6se66a+6lCdlngNPb1k+r2iRpIPVSAH8T2BgRGyJiFfAB4LbCNUlSbXqmCyIzD0bEh4Ev0hqG9ueZef8S3vKI3RQ9qh/r7seawbq7rR/rrrXmnrkIJ0nDppe6ICRpqBjAklRIXwZwRJwTEd+PiKmIuGKe7SMRcWO1/Z6IWN+27RNV+/cj4hd7veaIeE9E3BsR36lez+pWzUupu237T0XEjyPid7tWNEv+GdkUEV+LiPurz/3YXq45Io6JiG1Vrd+LiE90o96jqPtfRsR9EXGwGu/fvm1zRPyg+trcvaoXX3dEvLXt56MZEe9fdBGZ2VdftC7Q/RB4I7AK+Dbwpo59LgH+rFr+AHBjtfymav8RYEP1Pit6vOa3AT9ZLZ8BzPTDZ922/Wbgr4Hf7Ye6aV2YbgJvqdZP7oOfkd8AbqiW/wmwE1jfQ5/1emATcD3wvrb2k4CHqte11fLaPqj7p4GN1fJPAo8BaxZTRz+eAR+6ZTkznwdevmW53XnAtmr5ZuDsiIiq/YbMnM3Mh4Gp6v16tubM/FZm/mPVfj9wXER0a0aTpXzWRMT5wMO06u6mpdT9C0AzM78NkJlPZeaLPV5zAsdHxErgOOB54Jku1AwLqDszd2ZmE3ip49hfBLZn5tOZuRfYDhx2t1hNFl13Zv5DZv6gWv5HYDcwupgi+jGAx4BH29anq7Z598nMg8B+WmcyCzm2Dkupud2/Ae7LzNma6uy06Loj4gTgcuAPulBnp6V83j8NZER8sfrz8/e6UO+ceipHU/PNwLO0zsQeAf44M5+uu+DOmipH8ztV6vdx2f7tiDiT1hn0DxdTRM+MA9Zri4g3A5+mdYbWDz4FXJWZP65OiPvFSuBdwNuB54A7I+LezLyzbFmv6UzgRVp/Dq8F7o6IL2U1sZXqERHrgM8BmzOz8+x+QfrxDHghtywf2qf6s2w18NQCj63DUmomIk4DvgB8MDMX9T/tIi2l7ncAfxQRO4GPAL9f3WjTDUupexr4amY+mZnPAbcDP1N7xUur+TeAOzLzhczcDfw/oFtzLizld6rk9ANL+rcj4vXA3wGfzMyvL7qKbnR4L3Pn+UpanfUbeKXz/M0d+1zK3IsVN1XLb2buRbiH6M4FlqXUvKba/9f66bPu2OdTdPci3FI+77XAfbQuZq0EvgT8Uo/XfDnwv6rl42lN4bqpVz7rtn0/y+EX4R6uPvO11fJJfVD3KuBO4CNLrqMb32wNH965wD/Q6nf5ZNX2h8CvVsvH0rryPgV8A3hj27GfrI77PvDeXq8Z+I+0+vd2tH39RK/X3fEen6KLAbwMPyP/jtaFw+8Cf9TrNQMnVO330wrfj/fYZ/12Wn9ZPEvrjP3+tmP/ffX9TAEf6oe6q5+PFzp+J9+6mBq8FVmSCunHPmBJGggGsCQVYgBLUiEGsCQVYgBLUiEGsCQVYgBLUiH/H8Opdkf4FpElAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.displot(terminalY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "id": "03224ce8-d390-4105-86e5-f52b6c9a858d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<seaborn.axisgrid.FacetGrid at 0x7fa4f0c83fd0>"
      ]
     },
     "execution_count": 297,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX9klEQVR4nO3df5Bd5X3f8fc3yJI9yNbyQ10kxFTIYuwmtEZUpuvYk4lhnBEkNaTjUHvcoPGQaqZgxq47dnA90zYzyYzdaYsxU8gQcCyS1DahdkVcBocATmknxha+CC6GlLVqjbSAJAjI2MYLsr794z4r312vpN3Vnvvcvff9mtm55zzn3Hu/PMt+dO5zz3lOZCaSpN77hdoFSNKwMoAlqRIDWJIqMYAlqRIDWJIqWVa7gJOxZcuWvPfee2uXIUlTYj47L+kj4Oeff752CZK0YEs6gCVpKTOAJakSA1iSKjGAJakSA1iSKjGAJakSA1iSKjGAJakSA1iSKjGAJakSA1iSKjGAJakSA1iSKlnS01FK/WBycpJWqzWtbdOmTaxYsaJSRVoqDGDpJLVaLa67eQer1m4A4NAzu7npGhgbG6tcmfqdASwtglVrN3DmhvNrl6ElxjFgSarEAJakSgxgSarEAJakSgxgSaqk0QCOiO9HxOMR8WhE7Cxtp0fEfRHxdHk8rbRHRHwuIsYj4rGIuLDJ2iSptl4cAb87My/IzM1l/Xrg/sw8D7i/rANcCpxXfrYBt/SgNkmqpsYQxOXA9rK8Hbiiq/2O7PgmMBIRayrUJ0k90XQAJ/CXEfFIRGwrbaOZ+WxZfg4YLctnA3u7nruvtE0TEdsiYmdE7Dx48GBTdUtS45q+Eu5dmTkREX8PuC8inuremJkZETmfF8zMW4FbATZv3jyv50pSP2n0CDgzJ8rjAeCrwEXA/qmhhfJ4oOw+AZzT9fR1pU2SBlJjARwRp0bEG6eWgV8D2sDdwNay21ZgR1m+G7iqnA0xBhzqGqqQpIHT5BDEKPDViJh6n/+WmfdGxLeBOyPiamAPcGXZ/x7gMmAc+DHwoQZrk6TqGgvgzNwNvG2W9heAS2ZpT+DapuqRFsr5ftUUp6OUTsD5ftUUA1iag/nM93vk8Gu02+1pbR4xazYGsLTIXt6/lxv3vMLoeOcsSY+YdSwGsNSAlWet9w4ZOiFnQ5OkSgxgSarEIQipYbN9KQd+MScDWGrczC/lwC/m1GEAS/M084i23W6TR44/L5Rfymk2BrA0TzOPaCd2PcTIRm/govkzgKUF6D6iPTSxu3I1Wqo8C0KSKvEIWEPPyXZUiwGsoedkO6rFAJaY32Q70mJxDFiSKjGAJakSA1iSKjGAJakSv4STZljIpcbSQhjA0gxeaqxeMYClWXipsXrBMWBJqsQAlqRKDGBJqsQAlqRKDGBJqsQAlqRKDGBJqsQAlqRKDGBJqsQAlqRKvBRZQ2W2+7852Y5qMYA1VGbe/w2cbEf1GMAaOjPv/+ZkO6rFMWBJqsQAlqRKDGBJqsQAlqRK/BJOA23maWeecqZ+YgBroM087cxTztRPDGANvO7TzjzlTP3EAJYqOHL4Ndrt9rS2TZs2sWLFikoVqQYDWKrg5f17uXHPK4yOd8ajDz2zm5uugbGxscqVqZcMYKmSlWetn3ZFnoaPp6FJUiUGsCRVYgBLUiUGsCRVYgBLUiWNB3BEnBIRrYj4Wlk/NyIejojxiPhyRCwv7SvK+njZvr7p2iSppl4cAX8EeLJr/TPADZm5EXgRuLq0Xw28WNpvKPtJ0sBqNIAjYh3w68BtZT2Ai4G7yi7bgSvK8uVlnbL9krK/JA2kpo+APwt8AjhS1s8AXsrMw2V9H3B2WT4b2AtQth8q+08TEdsiYmdE7Dx48GCDpUtSsxoL4Ij4DeBAZj6ymK+bmbdm5ubM3Lx69erFfGlJ6qkmL0V+J/DeiLgMeD3wJuBGYCQilpWj3HXARNl/AjgH2BcRy4BVwAsN1idJVTV2BJyZn8zMdZm5Hng/8EBmfhB4EHhf2W0rsKMs313WKdsfyExnzpY0sGqcB/y7wMciYpzOGO/tpf124IzS/jHg+gq1SVLP9GQ2tMz8BvCNsrwbuGiWfX4C/FYv6pGkfuCVcJJUiQEsSZUYwJJUiQEsSZUYwJJUiQEsSZUYwJJUiXdF1kCZnJyk1WodXW+32+QRL6hUfzKANVBarRbX3byDVWs3ADCx6yFGNl5YuSppdgawBs6qtRs4c8P5ABya2F25GunYHAOWpEo8AtaS5pivljIDWEuaY75aygxgLXmDMOZ75PBrtNvtaW2bNm1ixYoVlSpSLxjAUh94ef9ebtzzCqPjneGTQ8/s5qZrYGxsrHJlapIBLPWJlWetP3okr+HgWRCSVIkBLEmVGMCSVIljwFIf8qyI4WAAS33IsyKGgwEs9SnPihh8jgFLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRVYgBLUiUGsCRV4i2JtKRMTk7SarWOrrfbbfJIVqxIWjgDWEtKq9Xiupt3sGrtBgAmdj3EyMYLK1clLYwBrCVn1doNR29WeWhid+VqesPb1A8mA1haArxN/WAygKUlwtvUDx7PgpCkSgxgSarEAJakSgxgSarEAJakShoL4Ih4fUR8KyJ2RcQTEfF7pf3ciHg4IsYj4ssRsby0ryjr42X7+qZqk6R+0OQR8CRwcWa+DbgA2BIRY8BngBsycyPwInB12f9q4MXSfkPZT5IGVmMBnB0/LKuvKz8JXAzcVdq3A1eU5cvLOmX7JRERTdUnSbU1OgYcEadExKPAAeA+4HvAS5l5uOyyDzi7LJ8N7AUo2w8BZ8zymtsiYmdE7Dx48GCT5UtSoxoN4Mz8aWZeAKwDLgLeugiveWtmbs7MzatXrz7Zl5OkanpyKXJmvhQRDwLvAEYiYlk5yl0HTJTdJoBzgH0RsQxYBbzQi/rUv5x+UoOssQCOiNXAayV83wC8h84Xaw8C7wO+BGwFdpSn3F3W/6ZsfyAz/Usbck4/qUHW5BHwGmB7RJxCZ6jjzsz8WkR8F/hSRPw+0AJuL/vfDvxJRIwDfwe8v8HatIQM4/STGg5zCuCIeGdm/p8TtXXLzMeATbO076YzHjyz/SfAb82lHkkaBHP9Eu6mObZJkubouEfAEfEO4JeB1RHxsa5NbwJOabIwSRp0JxqCWA6sLPu9sav9B3S+KJMkLdBxAzgz/xr464j4Qmbu6VFNkjQU5noWxIqIuBVY3/2czLy4iaIkaRjMNYD/HPhD4Dbgp82VI0nDY64BfDgzb2m0EkkaMnM9De0vIuKaiFgTEadP/TRamSQNuLkeAW8tjx/vaktgw+KWI0nDY04BnJnnNl2IpLk7cvg12u32tLZNmzaxYsWKShVpIeZ6KfJVs7Vn5h2LW46kuXh5/15u3PMKo+Od+aoOPbObm66BsbGxypVpPuY6BPH2ruXXA5cA3wEMYKmSlWetPzpJkZamuQ5BXNe9HhEjdKaTlCQt0ELviPEjwHFhSToJcx0D/gs6Zz1AZxKefwDc2VRRkjQM5joG/J+6lg8DezJzXwP1SNLQmNMQRJmU5yk6M6KdBrzaZFGSNAzmFMARcSXwLTp3rLgSeDginI5Skk7CXIcgPgW8PTMPwNEbbv4VcFdThUnSoJvrWRC/MBW+xQvzeK4kaRZzPQK+NyK+DnyxrP9z4J5mSpKk4XCie8JtBEYz8+MR8c+Ad5VNfwP8WdPFSdIgO9ER8GeBTwJk5leArwBExD8s2/5pg7VpyExOTtJqtaa1tdtt8kge4xnS0naiAB7NzMdnNmbm4xGxvpmSNKxarRbX3byDVWt/NsvpxK6HGNl4YcWqpOacKIBHjrPtDYtYhwTAqrUbpk0wc2hid8VqpGadKIB3RsS/zMw/6m6MiN8BHmmuLEnzMdv8wOAcwf3uRAH8UeCrEfFBfha4m4HlwG82WJekeZg5PzA4R/BScNwAzsz9wC9HxLuBqc+F/zMzH2i8Mknz4vzAS89c5wN+EHiw4Vokaah4NZskVWIAS1IlBrAkVWIAS1IlBrAkVWIAS1IlBrAkVTLX+YAlLTGzXZ7spcn9xQCWBtTMy5O9NLn/GMDSAPPy5P7mGLAkVeIRsKqZeQcM736hYWMAq5qZd8Dw7hcaNgawquq+A4Z3v9CwcQxYkioxgCWpEgNYkioxgCWpEgNYkioxgCWpEgNYkippLIAj4pyIeDAivhsRT0TER0r76RFxX0Q8XR5PK+0REZ+LiPGIeCwiPCNf0kBr8gj4MPBvMvMXgTHg2oj4ReB64P7MPA+4v6wDXAqcV362Abc0WJskVddYAGfms5n5nbL8MvAkcDZwObC97LYduKIsXw7ckR3fBEYiYk1T9UlSbT0ZA46I9cAm4GFgNDOfLZueA0bL8tnA3q6n7SttM19rW0TsjIidBw8ebK5oSWpY4wEcESuB/w58NDN/0L0tMxOY1/RXmXlrZm7OzM2rV69exEolqbcaDeCIeB2d8P2zzPxKad4/NbRQHg+U9gngnK6nryttkjSQmjwLIoDbgScz8790bbob2FqWtwI7utqvKmdDjAGHuoYqJGngNDkd5TuB3wYej4hHS9u/BT4N3BkRVwN7gCvLtnuAy4Bx4MfAhxqsTZKqayyAM/N/A3GMzZfMsn8C1zZVjyT1G6+Ek6RKDGBJqsQAlqRKDGBJqsQAlqRKDGBJqsTb0qsxk5OTtFqtaW2bNm1ixYoVlSqS+osBrMa0Wi2uu3kHq9ZuAODQM7u56RoYGxurXNlwOnL4Ndrt9rQ2/0GsywBWo1at3cCZG86vXYaAl/fv5cY9rzA63pn/yn8Q6zOA1TMzj8Da7TZ5ZF6T4ekkrTxrvf8g9hEDWD0z8whsYtdDjGz0zlMaXgaweqr7COzQxO7K1Uh1eRqaJFViAEtSJQawJFViAEtSJQawJFViAEtSJQawJFViAEtSJQawJFViAEtSJQawJFViAEtSJQawJFViAEtSJQawJFXifMDSkPIecfUZwNKQ8h5x9RnA0hDzHnF1OQYsSZUYwJJUiQEsSZUYwJJUiV/CacEmJydptVpH11999VUAli9fDkC73SaPZJXapKXAANaCtVotrrt5B6vWbgBgYtdDLFt5OqNv/qWj6yMbL6xZotTXDGCdlFVrNxw9jenQxG6WjYxOW5d0bI4BS1IlHgFLArw0uQYDWBLgpck1GMCSjvLS5N5yDFiSKjGAJakSA1iSKjGAJakSA1iSKjGAJakSA1iSKjGAJakSA1iSKmksgCPi8xFxICLaXW2nR8R9EfF0eTyttEdEfC4ixiPisYhwDkNJA6/JI+AvAFtmtF0P3J+Z5wH3l3WAS4Hzys824JYG65KkvtBYAGfm/wL+bkbz5cD2srwduKKr/Y7s+CYwEhFrmqpNkvpBr8eARzPz2bL8HDBals8G9nbtt6+0SdLAqvYlXGYmMO8bhkXEtojYGRE7Dx482EBlktQbvQ7g/VNDC+XxQGmfAM7p2m9dafs5mXlrZm7OzM2rV69utFhJalKvA/huYGtZ3grs6Gq/qpwNMQYc6hqqkKSB1NiE7BHxReBXgTMjYh/w74FPA3dGxNXAHuDKsvs9wGXAOPBj4ENN1SVJ/aKxAM7MDxxj0yWz7JvAtU3VImn+vEdc87wlkeZkcnKSVqs1ra3dbpNH5v09qpYI7xHXPANYc9Jqtbju5h2sWrvhaNvErocY2ehFi4PMe8Q1ywDWnK1au2HaH+Ohid0Vq5GWPifjkaRKDGBJqsQAlqRKDGBJqsQAlqRKDGBJqsQAlqRKDGBJqsQAlqRKDGBJqsQAlqRKnAtCs5o5+5kzn0mLzwDWrGbOfubMZ9LiM4B1TN2znznzmWaboB2cpP1kGMCS5mTmBO3gJO0nywCWNGdO0L64PAtCkioxgCWpEgNYkioxgCWpEgNYkioxgCWpEgNYkirxPGBJCzbb1XFeGTd3BrCkBZt5dZxXxs2PATykZs52Bh65aGG8Om7hDOAhNXO2sxf3/l/+1bvbnH9+5w/J6Sel5hnAQ2zmbGc3fr199KOk009KzTOAh8RcJljv/ijp9JNS8wzgATVb4P7hN8YZWfdmwCNcqR8YwAPqWHe08AhXTfK0tPkxgAeYd7RQr3la2vwYwJIWVfd3CR4RH58BLKkxHhEfnwEsqVFeqHFsTsYjSZUYwJJUiUMQA2C2eR28lFjqfwbwAJh5zi94oYW0FBjAA6L7nF/wvF/1J09Lm84AltQznpY2nQEsqac8Le1nDGBJ1Qz7kIQBLKmaYR+SMICXoLnM7SstFcM8JGEALwHO7SsNJgN4CXBuXw2L2caEYXDHhQ3gPjDzCPfVV18FYPny5UDniPdNZ51r4GrgzRwThp+/YSwMTiD3VQBHxBbgRuAU4LbM/HTlkhoxlyGFZStPZ/TNv3R03SEGDYuZY8Izbxg784u62S7FXyoB3TcBHBGnAP8VeA+wD/h2RNydmd9drPeY7RcFx/9lnejo9ETPn+01jjWG232Eu2xk1CNeqTjeJO8z/57mG9ALyYXF0jcBDFwEjGfmboCI+BJwObBoAdxqtbjqUzdw6hlrjrb96IVn+cQH3jPt4023drvNf/zifUef8/zuNqe84Y2ctubvz+n5x3qNkXOn7//D577P86ee2ll+foJlP3llweuL8Rq+x9J+z0F+j2ef+Ca///DLnLamE8Kz/T3NDOjuv7+Zf7Mzt0/tc8cf/OvGT4eLzP44fSki3gdsyczfKeu/DfyTzPzwjP22AdvK6luAv13kUs4Enl/k11wI65jOOqbrhzr6oQborzqeyswtc31CPx0Bz0lm3grc2tTrR8TOzNzc1Otbh3UMSh39UEMf1jHn8IX+mpB9Ajina31daZOkgdRPAfxt4LyIODcilgPvB+6uXJMkNaZvhiAy83BEfBj4Op3T0D6fmU9UKKWx4Y15so7prGO6fqijH2qAJVxH33wJJ0nDpp+GICRpqBjAklTJ0AZwRLw+Ir4VEbsi4omI+L3Sfm5EPBwR4xHx5fKFYI06vhAR/y8iHi0/FzRZR1c9p0REKyK+VtZ72h/HqaPn/RER34+Ix8v77Sxtp0fEfRHxdHk8rVId/yEiJrr647Ie1DESEXdFxFMR8WREvKNSf8xWR0/7IyLe0vVej0bEDyLio/Ptj6ENYGASuDgz3wZcAGyJiDHgM8ANmbkReBG4ulIdAB/PzAvKz6MN1zHlI8CTXeu97o9j1QF1+uPd5f2mzjO9Hrg/M88D7i/rNeqAzu9lqj/u6UENNwL3ZuZbgbfR+f3U6I/Z6oAe9kdm/u3UewH/GPgx8FXm2R9DG8DZ8cOy+rryk8DFwF2lfTtwRaU6ei4i1gG/DtxW1oMe98dsdfSZy+n0A/SoP/pBRKwCfgW4HSAzX83Ml+hxfxynjpouAb6XmXuYZ38MbQDD0Y+5jwIHgPuA7wEvZebhsss+4Oxe15GZD5dNfxARj0XEDRHRi6mdPgt8AjhS1s+gQn/MUseUXvdHAn8ZEY+US+ABRjPz2bL8HDBaqQ6AD5f++HwPPvqfCxwE/rgMDd0WEafS+/44Vh3Q2/7o9n7gi2V5Xv0x1AGcmT8tHyHW0ZkM6K39UEdEnA98stTzduB04HebrCEifgM4kJmPNPk+J1FHT/ujeFdmXghcClwbEb/SvTE753D24tPKbHXcAryZzrDVs8B/briGZcCFwC2ZuQn4ETM+XveoP45VR6/7A4Dynch7gT+fuW0u/THUATylfIR5EHgHMBIRUxeo9PRy6K46tmTms2V4YhL4Yzr/QDTpncB7I+L7wJfoDD3cSO/74+fqiIg/rdAfZOZEeTxAZ3zvImB/RKwBKI8HatSRmfvLP9xHgD+i+f7YB+zr+nR2F50g7HV/zFpHhf6YcinwnczcX9bn1R9DG8ARsToiRsryG+jMQ/wknQB8X9ltK7CjQh1Pdf0Sg8440s/fp2URZeYnM3NdZq6n85Hqgcz8ID3uj2PU8S963R8RcWpEvHFqGfi18p530+kH6M3/H7PWMdUfxW/S/P8fzwF7I+ItpekSOlPF9rQ/jlVHr/ujywf42fADzLc/MnMof4B/BLSAx+j8sv5dad8AfAsYp/OxYkWlOh4AHi9tfwqs7GHf/CrwtRr9cZw6etof5b97V/l5AvhUaT+DzrfbTwN/BZxeqY4/Kf3xWPmjX9OD38cFwM7ynv8DOK3X/XGcOmr0x6nAC8CqrrZ59YeXIktSJUM7BCFJtRnAklSJASxJlRjAklSJASxJlRjAklSJASxJlfx/XOOjCDfu9cMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.displot(terminalS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "id": "31aac4cf-552e-4e9c-8738-23fb5b25265e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/bv/j678dvtd34bc5mtymgsd2c6w0000gn/T/ipykernel_74024/1735570090.py:10: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  convergeS.append(ES/i)\n",
      "/var/folders/bv/j678dvtd34bc5mtymgsd2c6w0000gn/T/ipykernel_74024/1735570090.py:11: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  convergeY.append(EY/i)\n"
     ]
    }
   ],
   "source": [
    "convergeS = list()\n",
    "convergeY = list()\n",
    "ES = 0\n",
    "EY = 0\n",
    "\n",
    "for i in range(n):\n",
    "    ES += terminalS[i]\n",
    "    EY += terminalY[i]\n",
    "    \n",
    "    convergeS.append(ES/i)\n",
    "    convergeY.append(EY/i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 299,
   "id": "271acc67-4bfd-43a5-b5ed-36c5aed566c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of paths used is:  10000\n",
      "The final estimate of expected S is:  50.10907720968314\n",
      "The final estimate of expected Y is:  0.040115502019253305\n",
      "The standard error of expected S is:  0.6000622771919076\n",
      "The standard error of expected Y is:  0.000624858146651417\n"
     ]
    }
   ],
   "source": [
    "print(\"The number of paths used is: \",n)\n",
    "print(\"The final estimate of expected S is: \",expectedS)\n",
    "print(\"The final estimate of expected Y is: \",expectedY)\n",
    "print(\"The standard error of expected S is: \",np.std(np.array(convergeS[1:])))\n",
    "print(\"The standard error of expected Y is: \",np.std(np.array(convergeY[1:])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 300,
   "id": "fe6a0701-9000-4da5-a2db-3dd44d305dc3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 300,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUvklEQVR4nO3df4xdZ33n8ff3zvU4sUMTO5kak8Q4qFHaiBUhjFAQFFWYH0laNVGFoqBKmG621m6rtrR/UCP+qFaqVKjQ7hbtCpot7ZrdkgIp1BGU0tRlu1qJTRlTNjEJiU1+EGf9Y/LLITGxZ+Z+94/zzPje4Yw9M3ecmXPyfklX99zn/Hqee8Yfn/Occ+6JzESS1C6d1a6AJGnlGe6S1EKGuyS1kOEuSS1kuEtSC3VXuwIAl112WW7fvn21qyFJjbJ///6nM3OsbtyaCPft27czMTGx2tWQpEaJiCcWGnfObpmI+POIOB4RB/rKNkfEvRFxsLxvKuUREZ+KiEMRcX9EXL8yTZAkLcVi+tz/G3DjvLLdwL7MvBrYVz4D3ARcXV67gE+vTDUlSUtxznDPzP8FPDuv+BZgTxneA9zaV/65rPwf4JKI2LpCdZUkLdJyr5bZkplHyvBRYEsZvhx4sm+6w6XsJ0TEroiYiIiJycnJZVZDklRn6Eshs/pxmiX/QE1m3pmZ45k5PjZWe7JXkrRMyw33Y7PdLeX9eCl/Criyb7orSpkk6RW03HC/B9hZhncCe/vKP1iumrkBONHXfSNJeoUs5lLIu4BvAddExOGIuAP4OPCeiDgIvLt8Bvhb4FHgEPBfgd84L7Uu/ujrD3HbZ77F/ieeO5+rkaTGOedNTJn5gQVG7aiZNoHfHLZSi/XI0R/xz48/y4kfn36lVilJjeBvy0hSCxnuktRCrQh3nxQoSYMaHe4RsdpVkKQ1qdHhLkmqZ7hLUgu1Itztc5ekQY0Od3vcJaleo8NdklTPcJekFjLcJamFWhHunk+VpEGNDnfvYZKkeo0Od0lSPcNdklqoFeGe3sUkSQMaHu52uktSnYaHuySpjuEuSS3UinC3x12SBjU63L3OXZLqNTrcJUn1DHdJaqFWhLuXuUvSoEaHu13uklSv0eEuSapnuEtSCxnuktRCLQl3z6hKUr9Gh7s3MUlSvUaHuySp3lDhHhG/ExEHIuJ7EfHhUrY5Iu6NiIPlfdOK1FSStGjLDveIeCPw68BbgTcBvxQRPwPsBvZl5tXAvvL5vPImJkkaNMye+88B92XmycycBv4J+BXgFmBPmWYPcOtQNTyL8DYmSao1TLgfAH4+Ii6NiA3AzcCVwJbMPFKmOQpsGbKOkqQl6i53xsx8KCI+Afw98BLwXWBm3jQZEbWdJhGxC9gFsG3btuVWQ5JUY6gTqpn52cx8S2a+E3gOeAQ4FhFbAcr78QXmvTMzxzNzfGxsbJhqeJW7JM0z7NUyP13et1H1t38euAfYWSbZCewdZh1nX//5WrIkNduyu2WKv46IS4Ep4Dcz8/mI+DjwxYi4A3gCuG3YSkqSlmaocM/Mn68pewbYMcxyJUnD8Q5VSWqhVoS7NzFJ0qBGh7snVCWpXqPDXZJUz3CXpBZqRbintzFJ0oBGh7s/HCZJ9Rod7pKkeoa7JLVQK8Ld69wlaVCzw90ud0mq1exwlyTVMtwlqYVaEe52uUvSoEaHu13uklSv0eEuSapnuEtSCxnuktRCrQj39C4mSRrQ6HAPn9YhSbUaHe6SpHqGuyS1kOEuSS3U6HC3x12S6jU63CVJ9Qx3SWqhVoS7l7lL0qBGh7uXuUtSvUaHuySpnuEuSS3UinBPH9chSQMaHe52uUtSvaHCPSJ+NyK+FxEHIuKuiLggIq6KiPsi4lBEfCEiRleqspKkxVl2uEfE5cBvA+OZ+UZgBLgd+ATwHzPzZ4DngDtWoqKSpMUbtlumC1wYEV1gA3AEeBdwdxm/B7h1yHVIkpZo2eGemU8BnwR+SBXqJ4D9wPOZOV0mOwxcXjd/ROyKiImImJicnFxuNUpdhppdklpnmG6ZTcAtwFXA64CNwI2LnT8z78zM8cwcHxsbW24dljWfJLXdMN0y7wYey8zJzJwCvgy8HbikdNMAXAE8NWQdJUlLNEy4/xC4ISI2RLULvQN4EPgm8P4yzU5g73BVlCQt1TB97vdRnTj9DvBAWdadwO8DvxcRh4BLgc+uQD3PUZfzvQZJapbuuSdZWGb+AfAH84ofBd46zHIXyx53SarX6DtUJUn1DHdJaqFWhLtd7pI0qNnhbqe7JNVqdrhLkmoZ7pLUQoa7JLVQK8I9vYtJkgY0OtzDM6qSVKvR4S5Jqme4S1ILtSLc7XGXpEGNDnef1SFJ9Rod7pKkeoa7JLVQO8LdTndJGtDocLfLXZLqNTrcJUn1DHdJaqFWhHva6S5JAxod7l7nLkn1Gh3ukqR6hrsktZDhLkkt1Ipw91kdkjSo0eHuwzokqV6jw12SVM9wl6QWakW42+UuSYMaHe7exCRJ9Rod7pKkeoa7JLXQssM9Iq6JiO/2vV6IiA9HxOaIuDciDpb3TStZ4Tpe5y5Jg5Yd7pn5cGZel5nXAW8BTgJfAXYD+zLzamBf+Xxe2OcuSfVWqltmB/CDzHwCuAXYU8r3ALeu0DokSYu0UuF+O3BXGd6SmUfK8FFgS90MEbErIiYiYmJycnKFqiFJghUI94gYBX4Z+NL8cZmZLHAZembemZnjmTk+NjY2bDUkSX1WYs/9JuA7mXmsfD4WEVsByvvxFVjHWfkkJkkatBLh/gHOdMkA3APsLMM7gb0rsI4FeEZVkuoMFe4RsRF4D/DlvuKPA++JiIPAu8tnSdIrqDvMzJn5EnDpvLJnqK6ekSStklbcoepNTJI0qNHh7k1MklSv0eEuSapnuEtSC7Ui3O1yl6RBjQ53u9wlqV6jw12SVM9wl6QWake4e6G7JA1odLh7nbsk1Wt0uEuS6hnuktRChrsktVArwt3TqZI0qNHhHt7GJEm1Gh3ukqR6hrsktVArwt17mCRpUKPD3ZuYJKleo8NdklTPcJekFmpFuKed7pI0oNHhbpe7JNVrdLhLkuoZ7pLUQq0Id3vcJWlQo8M9vNBdkmo1OtwlSfUMd0lqIcNdklqoFeHuPUySNGiocI+ISyLi7oj4fkQ8FBFvi4jNEXFvRBws75tWqrKSpMUZds/9T4C/y8yfBd4EPATsBvZl5tXAvvJZkvQKWna4R8TFwDuBzwJk5unMfB64BdhTJtsD3DpcFSVJSzXMnvtVwCTwFxHxLxHxZxGxEdiSmUfKNEeBLXUzR8SuiJiIiInJyckhquFNTJI03zDh3gWuBz6dmW8GXmJeF0xWP9dYm72ZeWdmjmfm+NjY2LIq4D1MklRvmHA/DBzOzPvK57upwv5YRGwFKO/Hh6uiJGmplh3umXkUeDIirilFO4AHgXuAnaVsJ7B3qBpKkpasO+T8vwX8ZUSMAo8Cv0b1H8YXI+IO4AngtiHXcU4+rEOSBg0V7pn5XWC8ZtSOYZa7WOHjOiSpVivuUJUkDTLcJamFDHdJaqFGh/vsde6eT5WkQY0O95FOle49012SBjQ63Gf33HtmuyQNaHS4d8I9d0mq0/Bwr969iUmSBjU83Kt0n+mtckUkaY1pRbjbLSNJg1oR7nbLSNKghod79e7VMpI0qNnhXtJ9xj13SRrQ7HC3z12SajU83Kt3s12SBjU83Mueu53ukjSg0eHuzw9IUr1Gh7t97pJUr9HhPhvq33782VWuiSStLY0O9z/82kMAfO//vbDKNZGktaXR4S5JqtfocH/TlZesdhUkaU1qdLi/7uILVrsKkrQmNTrcd9/0s3PDM14PKUlzGh3ur790IxtGRwA4eXp6lWsjSWtHo8Md4OTpGQBO/HhqlWsiSWtH48N91rs++U+rXQVJWjNaE+6nfdaeJM1pTbgDTBvwkgS0INz/atcNc8Of+9YTq1gTSVo7Gh/ul24cnRv+w689uIo1kaS1o/Hhvv2yjXPDvYTtu7/GvQ8es4tG0qtad5iZI+Jx4EfADDCdmeMRsRn4ArAdeBy4LTOfG66aC1s30uFrv/0OfvFT/3uu7Nc/NwHAY390M5MvnmLsovXE7I+/L0Fm8ujTL/Hksyc5/qNTvPDjKaZmkis2XcjWiy8gIrjmta/hovVDfY2StOIih/gt9BLu45n5dF/ZHwPPZubHI2I3sCkzf/9syxkfH8+JiYmlV+CRR+CBBwDY99AxvrT/8Dlned+1r+XIiR9z/1MnAPhXl1/M0y+e4siJlwHojgTj2zYx1UsOP3uSYz86dc5ljl20ns0bR7lofZcLR0eIOPPov14vWb9uhJ+6sAsJF6zr0Mvq+vzpXhJU03Y68OKpaU6cnGa0G2RWD/7uBHQ7HbojAQkRwUgHuiOduXER1W/bdyLodKBDcHJqhlNT1ToAuuWZhCOdYLRbHbBN95LsJdEJpmeSqZkep6Z7ZCbdkQ6j3eDCbpduB3pArwc9quVFWVYH6HY7rCt1jLLc6V6e+Z39rI6qpns9pmZ6cw9X6URV3uslM5n0epT36nNmEgQRVbujzBMRc22mr2yml2TC+nUdRkc69Mr83Q4QQS+TDlVbMsuLM3//vcxSxtx3NlJe3U6QQJYDws68Y94gSJKZHvR6Sz9qXO6/wlf6vuxhH50w+29j/nf/E9OV73MpOuXvYlbfn19fYf/gwsufP2aml4t64ltE0J+pOfd+pjJzZWVg89vfys23vuOcy15gffszc7xu3PnY5bwF+IUyvAf4n8BZw33Z9u6Fj3wEgB3ldU5/c15qIknLcvfzH4NlhvvZDBvuCfx9RCTwp5l5J7AlM4+U8UeBLXUzRsQuYBfAtm3blrf2D30IbryxdtRjT7/Iv/0f3+HW617Ho0+/xP2HT9ROd9MbXwvA1w8cBeC6Ky5mx7VbODXV4/WXbuDnXvtTrOsufGri5akZfvjMS7zw8jQvvDzFyVMzzGSybqSaJ6j2yJ958TTruh1Onp5iJDpsKHv43U6HTifIXnLh6AiXXbSeqV6v2mOMYLqXzPR6nJrOuT3d6V6PU1PVPvTsHsXsnu5Mr9oD3bBuhA3rR+iWXczZ8dMzcHp6hqTa+x+JoEfS7QSjIx0uHO0QEUzNJKemZ3jp1AzTvR4dgk7Zi4Vqr6PaC0+mZpKp6R5TvRkyq66y6sjizJ51J6j27rvVOqGqZ2f2SKTTqY4EOsG68t6Jau+tl1UXWWYyM7vXl2fKe2UvsEMQEbw8NcPp6Rk6nWrvfaYcRXQjmCEZIYjSjtk9vc7ckUHM7V3OZNWu6ogjiYi5h7LP3zmf3TMbHamWXS1taZY+x5DzvdIzliOj/u9+9m/kzF52dbTVy6yO2Ija1WUycIScVEeh88fXzdxfNH87DY47Y/bvdBFNnDvSnL/MvpK54QjYds1V517wMgzbLXN5Zj4VET8N3Av8FnBPZl7SN81zmbnpbMtZdreMJL2Kna1bZqirZTLzqfJ+HPgK8FbgWERsLSveChwfZh2SpKVbdrhHxMaIeM3sMPBe4ABwD7CzTLYT2DtsJSVJSzNMn/sW4Culb6kLfD4z/y4ivg18MSLuAJ4Abhu+mpKkpVh2uGfmo8CbasqfYZEXrkiSzo/G36EqSfpJhrsktZDhLkktZLhLUgsNdRPTilUiYpLqyprluAx4+pxTtYttfnWwza8Ow7T59Zk5VjdiTYT7MCJiYqE7tNrKNr862OZXh/PVZrtlJKmFDHdJaqE2hPudq12BVWCbXx1s86vDeWlz4/vcJUk/qQ177pKkeQx3SWqhRod7RNwYEQ9HxKHyvNZGiogrI+KbEfFgRHwvIn6nlG+OiHsj4mB531TKIyI+Vdp9f0Rc37esnWX6gxGxc6F1rhURMRIR/xIRXy2fr4qI+0rbvhARo6V8ffl8qIzf3reMj5byhyPifavUlEWJiEsi4u6I+H5EPBQRb2v7do6I3y1/1wci4q6IuKBt2zki/jwijkfEgb6yFduuEfGWiHigzPOpiEU8Dmv28WVNewEjwA+ANwCjwP8Frl3tei2zLVuB68vwa4BHgGuBPwZ2l/LdwCfK8M3A16me3HUDcF8p3ww8Wt43leFNq92+c7T994DPA18tn78I3F6GPwP8uzL8G8BnyvDtwBfK8LVl268Hrip/EyOr3a6ztHcP8G/K8ChwSZu3M3A58BhwYd/2/VDbtjPwTuB64EBf2YptV+Cfy7RR5r3pnHVa7S9liC/zbcA3+j5/FPjoatdrhdq2F3gP8DCwtZRtBR4uw38KfKBv+ofL+A9QPcuWuunW2gu4AtgHvAv4avnDfRrozt/GwDeAt5Xhbpku5m/3/unW2gu4uARdzCtv7XYu4f5kCaxu2c7va+N2BrbPC/cV2a5l3Pf7ygemW+jV5G6Z2T+aWYdLWaOVw9A3A/ex8MPGF2p7076T/wR8BJh93PSlwPOZOV0+99d/rm1l/IkyfZPafBUwCfxF6Yr6s/IUs9Zu56wexflJ4IfAEarttp92b+dZK7VdLy/D88vPqsnh3joRcRHw18CHM/OF/nFZ/ZfdmutWI+KXgOOZuX+16/IK6lIdun86M98MvER1uD6nhdt5E3AL1X9srwM2AjeuaqVWwWps1yaH+1PAlX2fryhljRQR66iC/S8z88uleKGHjS/U9iZ9J28HfjkiHgf+iqpr5k+ASyJi9glh/fWfa1sZfzHwDM1q82HgcGbeVz7fTRX2bd7O7wYey8zJzJwCvky17du8nWet1HZ9qgzPLz+rJof7t4Gry1n3UaqTL/escp2WpZz5/izwUGb+h75RCz1s/B7gg+Ws+w3AiXL49w3gvRGxqewxvbeUrTmZ+dHMvCIzt1Ntu3/MzF8Fvgm8v0w2v82z38X7y/RZym8vV1lcBVxNdfJpzcnMo8CTEXFNKdoBPEiLtzNVd8wNEbGh/J3Ptrm127nPimzXMu6FiLihfIcf7FvWwlb7JMSQJzBuprqy5AfAx1a7PkO04x1Uh2z3A98tr5up+hr3AQeBfwA2l+kD+C+l3Q8A433L+tfAofL6tdVu2yLb/wucuVrmDVT/aA8BXwLWl/ILyudDZfwb+ub/WPkuHmYRVxGscluvAybKtv4bqqsiWr2dgX8PfB84APx3qiteWrWdgbuozilMUR2h3bGS2xUYL9/fD4D/zLyT8nUvf35Aklqoyd0ykqQFGO6S1EKGuyS1kOEuSS1kuEtSCxnuktRChrsktdD/B6YcLBknDPB0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.lineplot(data = convergeS[1:],linewidth=2)\n",
    "sns.lineplot(data = np.ones(n)*expectedS,color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "id": "191b9bc8-b26a-4332-8285-1ecd53eb8098",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 301,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfWUlEQVR4nO3dfZQc1X3m8e/T3fOid0AMoCBhCZDtVQzBMMtLbOxsWNsi9lk550AQS2ySsGEdH3bj9ck6ONlwCH8kxps1sTcc2zhgY2wMWRzHk1iOEsAvya6taLCFQQTZg4BIAgchCQlJ89bdv/2j7gytnpKmZzSjkaaezzl9VHXrVnHv1DBP173V1YoIzMysmEoz3QAzM5s5DgEzswJzCJiZFZhDwMyswBwCZmYFVpnpBkzEqaeeGsuXL5/pZpiZnVAee+yxlyOiK29bSyEgaTXwSaAM/HlEfKxpewfwReAiYBdwTUQ8J+k64L83VD0fuDAiNkn6NrAE6E/b3hkRLx2pHcuXL6e3t7eVJpuZWSLp+cNtG3c4SFIZuBO4ElgFXCtpVVO1G4A9EXEucAdwO0BEfDkiLoiIC4D3Ac9GxKaG/a4b2T5eAJiZ2dRrZU7gYqAvIrZGxBDwALCmqc4a4N60/BBwhSQ11bk27WtmZseJVkLgTGBbw/r2VJZbJyKqwF5gcVOda4CvNJV9XtImSX+QExoASLpRUq+k3p07d7bQXDMza9UxuTtI0iXAwYh4sqH4uog4D7g8vd6Xt29E3BUR3RHR3dWVO69hZmaT1EoI7ACWNawvTWW5dSRVgEVkE8Qj1tJ0FRARO9K/rwL3kw07mZnZMdRKCGwEVkpaIamd7A96T1OdHuD6tHwV8GikJ9NJKgG/QsN8gKSKpFPTchvwHuBJzMzsmBr3FtGIqEq6CVhPdovoPRGxWdJtQG9E9AB3A/dJ6gN2kwXFiLcB2yJia0NZB7A+BUAZeBj43JT0yMzMWqYT6VHS3d3dMZnPCXzp+88zv6PCe9/cPJ9tZjb7SXosIrrzthXisREPbPwX/vrxF2a6GWZmx51ChICZmeUrTAicOINeZmbHTiFCQOR+Ds3MrPAKEQJmZpavMCFwIt0FZWZ2rBQiBPKfSmRmZoUIATMzy+cQMDMrsMKEgGcEzMzGKkQIeErAzCxfIULAzMzyFSYEfIeomdlYxQgB3yNqZparGCFgZma5ChMCHg0yMxurECHgwSAzs3yFCAEzM8tXmBDwA+TMzMZqKQQkrZa0RVKfpJtztndIejBt3yBpeSq/TtKmhldd0gVp20WSnkj7fEqavlt4fHOQmVm+cUNAUhm4E7gSWAVcK2lVU7UbgD0RcS5wB3A7QER8OSIuiIgLgPcBz0bEprTPp4HfBFam1+qj7o2ZmU1IK1cCFwN9EbE1IoaAB4A1TXXWAPem5YeAK3Le2V+b9kXSEmBhRHw/snGaLwLvnVwXzMxssloJgTOBbQ3r21NZbp2IqAJ7gcVNda4BvtJQf/s4xwRA0o2SeiX17ty5s4Xm5hxjUnuZmc1+x2RiWNIlwMGIeHKi+0bEXRHRHRHdXV1d09A6M7PiaiUEdgDLGtaXprLcOpIqwCJgV8P2tbx2FTBSf+k4xzQzs2nWSghsBFZKWiGpnewPek9TnR7g+rR8FfBoGutHUgn4FdJ8AEBEvAjsk3Rpmjt4P/D1o+rJOHyHqJnZWJXxKkREVdJNwHqgDNwTEZsl3Qb0RkQPcDdwn6Q+YDdZUIx4G7AtIrY2HfqDwBeAOcA302taTOPdp2ZmJ7RxQwAgItYB65rKbmlYHgCuPsy+3wYuzSnvBd40gbaamdkUK84nhv0IOTOzMQoRAh4MMjPLV4gQMDOzfIUJAd8dZGY2ViFCwDcHmZnlK0QImJlZvsKEgIeDzMzGKkQIyPcHmZnlKkQImJlZvsKEgD8sZmY2VjFCwKNBZma5ihECZmaWyyFgZlZghQkB3yJqZjZWIULAUwJmZvkKEQJmZpavMCHg0SAzs7EKEQJ+gJyZWb6WQkDSaklbJPVJujlne4ekB9P2DZKWN2w7X9L3JG2W9ISkzlT+7XTMTel12pT1yszMWjLudwxLKgN3Au8AtgMbJfVExFMN1W4A9kTEuZLWArcD10iqAF8C3hcRj0taDAw37Hdd+q7h6efxIDOzMVq5ErgY6IuIrRExBDwArGmqswa4Ny0/BFwhScA7gR9FxOMAEbErImpT0/TW+QFyZmb5WgmBM4FtDevbU1lunYioAnuBxcDrgZC0XtIPJH2kab/Pp6GgP0ihMYakGyX1SurduXNnC801M7NWTffEcAV4K3Bd+veXJV2Rtl0XEecBl6fX+/IOEBF3RUR3RHR3dXVNuiF+gJyZ2VithMAOYFnD+tJUllsnzQMsAnaRXTV8NyJejoiDwDrgQoCI2JH+fRW4n2zYaVr47iAzs3ythMBGYKWkFZLagbVAT1OdHuD6tHwV8GhEBLAeOE/S3BQObweeklSRdCqApDbgPcCTR98dMzObiHHvDoqIqqSbyP6gl4F7ImKzpNuA3ojoAe4G7pPUB+wmCwoiYo+kT5AFSQDrIuIbkuYB61MAlIGHgc9NQ/8a+jGdRzczOzGNGwIAEbGObCinseyWhuUB4OrD7PslsttEG8sOABdNtLGT5eEgM7N8hfjEsJmZ5XMImJkVWGFCwFMCZmZjFSIE/IlhM7N8hQgBMzPLV5gQCN8jamY2RiFCwLeImpnlK0QImJlZvsKEgAeDzMzGKkwImJnZWA4BM7MCK0wI+OYgM7OxChECh/nSMjOzwitECJiZWb7ChIBHg8zMxipECHgwyMwsXyFCwMzM8jkEzMwKrDgh4HtEzczGaCkEJK2WtEVSn6Sbc7Z3SHowbd8gaXnDtvMlfU/SZklPSOpM5Rel9T5Jn9I03sfpO0TNzPKNGwKSysCdwJXAKuBaSauaqt0A7ImIc4E7gNvTvhWyL5n/QET8LPALwHDa59PAbwIr02v10XbGzMwmppUrgYuBvojYGhFDwAPAmqY6a4B70/JDwBXpnf07gR9FxOMAEbErImqSlgALI+L7kT3o/4vAe4++O4fnwSAzs7FaCYEzgW0N69tTWW6diKgCe4HFwOuBkLRe0g8kfaSh/vZxjgmApBsl9Urq3blzZwvNzTnGpPYyM5v9Ksfg+G8F/i1wEHhE0mNkIdGSiLgLuAugu7vbb+jNzKZQK1cCO4BlDetLU1lunTQPsAjYRfYO/7sR8XJEHATWARem+kvHOeaU8s1BZmZjtRICG4GVklZIagfWAj1NdXqA69PyVcCjaax/PXCepLkpHN4OPBURLwL7JF2a5g7eD3x9CvqTyw+QMzPLN+5wUERUJd1E9ge9DNwTEZsl3Qb0RkQPcDdwn6Q+YDdZUBAReyR9gixIAlgXEd9Ih/4g8AVgDvDN9DIzs2OopTmBiFhHNpTTWHZLw/IAcPVh9v0S2W2izeW9wJsm0tijEb4/yMxsjEJ8YtiDQWZm+QoRAmZmlq8wIeC7g8zMxipECPjmIDOzfIUIATMzy+cQMDMrsMKEgOcEzMzGKkgIeFLAzCxPQULAzMzyFCYEPBpkZjZWIULAt4iameUrRAiYmVm+woRA+PYgM7MxChECHg0yM8tXiBAwM7N8DgEzswIrRAj47iAzs3yFCAEzM8vXUghIWi1pi6Q+STfnbO+Q9GDavkHS8lS+XFK/pE3p9ZmGfb6djjmy7bQp61UO3xxkZjbWuN8xLKkM3Am8A9gObJTUExFPNVS7AdgTEedKWgvcDlyTtj0TERcc5vDXpe8anlby/UFmZrlauRK4GOiLiK0RMQQ8AKxpqrMGuDctPwRcIXkk3szseNdKCJwJbGtY357KcutERBXYCyxO21ZI+qGk70i6vGm/z6ehoD84XGhIulFSr6TenTt3ttDcfOGnB5mZjTHdE8MvAmdFxJuBDwP3S1qYtl0XEecBl6fX+/IOEBF3RUR3RHR3dXVNqhG+JjEzy9dKCOwAljWsL01luXUkVYBFwK6IGIyIXQAR8RjwDPD6tL4j/fsqcD/ZsJOZmR1DrYTARmClpBWS2oG1QE9TnR7g+rR8FfBoRISkrjSxjKSzgZXAVkkVSaem8jbgPcCTR98dMzObiHHvDoqIqqSbgPVAGbgnIjZLug3ojYge4G7gPkl9wG6yoAB4G3CbpGGgDnwgInZLmgesTwFQBh4GPjfVnTu0H9N5dDOzE9O4IQAQEeuAdU1ltzQsDwBX5+z3VeCrOeUHgIsm2tjJ8pyAmVk+f2LYzKzAChMCHg0yMxurECHgTwybmeUrRAiYmVm+woSAv17SzGysYoSAR4PMzHIVIwTMzCxXYULAg0FmZmMVIgQ8GmRmlq8QIWBmZvmKEwIeDzIzG6MQIXBgsMr2V/pnuhlmZsedlh4gd6L71pbJfyOZmdlsVogrATMzy+cQMDMrsEKEwK9eehaL57XPdDPMzI47hQgB8M1BZmZ5ChECfpS0mVm+QoQA+CmiZmZ5WgoBSaslbZHUJ+nmnO0dkh5M2zdIWp7Kl0vql7QpvT7TsM9Fkp5I+3xKmr5vAvZ3DJuZ5Rs3BCSVgTuBK4FVwLWSVjVVuwHYExHnAncAtzdseyYiLkivDzSUfxr4TWBleq2efDfG5+sAM7OxWrkSuBjoi4itETEEPACsaaqzBrg3LT8EXHGkd/aSlgALI+L7kY3TfBF470Qb3ypfCJiZ5WslBM4EtjWsb09luXUiogrsBRanbSsk/VDSdyRd3lB/+zjHBEDSjZJ6JfXu3Dn5T/56SsDMbKzpnhh+ETgrIt4MfBi4X9LCiRwgIu6KiO6I6O7q6ppUI6ZxusHM7ITWSgjsAJY1rC9NZbl1JFWARcCuiBiMiF0AEfEY8Azw+lR/6TjHnFK+O8jMbKxWQmAjsFLSCkntwFqgp6lOD3B9Wr4KeDQiQlJXmlhG0tlkE8BbI+JFYJ+kS9PcwfuBr09Bf8zMbALGfYpoRFQl3QSsB8rAPRGxWdJtQG9E9AB3A/dJ6gN2kwUFwNuA2yQNA3XgAxGxO237IPAFYA7wzfSaNr4OMDMbq6VHSUfEOmBdU9ktDcsDwNU5+30V+OphjtkLvGkijZ0sTwmYmeUrzCeGfSlgZjZWIULAzw4yM8tXiBAAXwiYmeUpRAh4TsDMLF8hQgD8OQEzszyFCAFfCJiZ5StECIDnBMzM8hQiBDwnYGaWrxAhAH6KqJlZnkKEgJ8iamaWrxAhABCeFTAzG6MQIeDrADOzfIUIgcFqnYHh+kw3w8zsuFOIEPjC/3sOgD0Hhma2IWZmx5lChMCIoZqvBszMGhUqBColzw6YmTUqVAiYmdmhChECv/bzywE/OsLMrFlLISBptaQtkvok3ZyzvUPSg2n7BknLm7afJWm/pN9pKHtO0hOSNknqPeqeHME5p80HoO6PDZuZHWLcEJBUBu4ErgRWAddKWtVU7QZgT0ScC9wB3N60/RPkf5H8v4uICyKie8Itn4DRqQBngJnZIVq5ErgY6IuIrRExBDwArGmqswa4Ny0/BFyh9KwGSe8FngU2T0mLJ6GUHhtRdwiYmR2ilRA4E9jWsL49leXWiYgqsBdYLGk+8LvAH+YcN4C/k/SYpBsn2vCJGLkQ8HCQmdmhKtN8/FuBOyJif85D3N4aETsknQb8vaSnI+K7zZVSQNwIcNZZZ02qESNXAo4AM7NDtXIlsANY1rC+NJXl1pFUARYBu4BLgI9Leg74EPB7km4CiIgd6d+XgK+RDTuNERF3RUR3RHR3dXW11qtmKX/qHg8yMztEKyGwEVgpaYWkdmAt0NNUpwe4Pi1fBTwamcsjYnlELAf+FPijiPgzSfMkLQCQNA94J/Dk0XcnX8mPkjYzyzXucFBEVNO79/VAGbgnIjZLug3ojYge4G7gPkl9wG6yoDiS04GvpSGiCnB/RPztUfTjiDwnYGaWr6U5gYhYB6xrKrulYXkAuHqcY9zasLwV+LmJNPRojFwIOAPMzA5ViE8Me2LYzCxfIUJg5ErAw0FmZocqRAiMcAaYmR2qECHw2t1BTgEzs0aFCIHXhoNmth1mZsebYoRAuknUw0FmZocqRAiMPEU0PBxkZnaIQoTA6HCQv2LYzOwQhQiBkc8M+0rAzOxQhQiBbbsPAtCz6YUZbomZ2fGlECHwk5deBeCz3906wy0xMzu+FCIE2sqF6KaZ2YQV4q9jpfTao6S/+cSLM9gSM7PjSzFCoOFK4Le+/AN/uYyZWVKMECgd+qUyB4dro8sHBqv89eOeMDazYpru7xg+LpSbQuDAYJX5HVnX3/2pf+C5XQfZc3CI91+2fAZaZ2Y2cwp5JXDJHz0yOiR06vwOAD73D6/dObT85m+Mvm74wkZq9eDLG57nlYND7D04zIHB6mjdkdtPzcxORIW4EiDnO4Y3v7CPMxZ10vv8HgC27e7nL3+wnQ//xeOH1Hvk6Zc45/eyL1X7/a8d/muQn/vYu6ewwWZmx0YhrgTynhy34dldXPnJ7x5S1hgAv/6W5RP6T3zjRxO76ygiqNeD7XsO0j9UY6ha5zPfeYY/ffjH/ONPXuZj33ya8BPvzGyaqZU/NJJWA58k+6L5P4+IjzVt7wC+CFwE7AKuiYjnGrafBTwF3BoRf9LKMfN0d3dHb29vaz1r8PVNO/jtBzZNaJ/nPvZuXh0Y5pWDw1z+8W/xxjMW8L+vfTNzOyr8z799mr/a9AISfOiK13PHwz8e3W/zH76LeWm+Yccr/dzas5l/enY33a87mUeefmnCbW/2xjMW8PRPXx1d/8Dbz+GCZYv4X3/3Y37x35zGgcEqX/vBDkolcWCwyjld81l2ylzOPW0+23YfZG57had/uo/NL+wD4OqLlvLTfQN0ze/g/z7zMpVSidMWdvCG0xfwysFhnnxhL6ct6GD7nn5KEp1tJS5YdhJDtexBTJ2VMo9ueYn9A1VOmtvOxStOJgKGa3XecMYC9vVXmdtepmtBB51tZXqf282Sk+Ywt63MG5cs5IVX+tm1f5BzTpvPKfPaOWlOOyu65o3O2ZjZ0ZP0WER0524bLwQklYEfA+8AtgMbgWsj4qmGOh8Ezo+ID0haC/xyRFzTsP0hsm902RARf9LKMfNMNgRi3Tr+y2e/c9jtt615E0/u2Mtf9G4D4H+8exVnLOp8bf/0M1LOsBLAAxv/hX/8ycuj6+ecNp9nXto/4XYCdC8/md7n9kxq3wkR437HjpRdRJVLotZ0W+3COW0cGKyOlp+/dBG7Dwzx6kCVvf3Do/U6KiUGqxN/ct+iOW0smtPGgs425raXqdaD9nKJxfPbOWNhJ6WSOLtrHgs72w7Zr1qvs/fgMAEsntd+yDmLCAaGawSwf6BKe6XE3PYK7ZXsgrheD0ql/HN8JNVanf7hGrv3D1GLYE5bmVcHq+x8dZDBap1qLRiu1dnbP0T/cJ0lizopCeZ3VJDEy/sHgSw4K6XS6DHmtpeZ11EZ/Rm2l0scGKoyMFSjs71MuSQ6K2XqEZQk5nVUsn3ay3S0lemolMb0X9KE+hkRDNeCwWqN4VowVK3zr/sG2D9YJQIGqjVKEpWSKKdXR6VEe6VEvR4MVLOfDXHo71NbWbSVS1RrwYGhbI6tJFFS9r0fEYz2Kwiq9aCtVKJar1Otx+h/Y+TXOAJq9eCV/iEGh+vUA4aqdeZ3VpjfUaGtPLa/9YCh4RqlUon2imivlGgrl5jbXqazrUy9HnS2lUf7NpB+jweHa9QJdu0fYiiVlUsabX+lVKJUGumPqEd2nOb/5er1rF+1elCLoFp7rd2VcvZz7Ejnd2C4xqUf+c+oMrk3R0cbApeRvYN/V1r/KEBE/HFDnfWpzvckVYCfAl0REZLeC7wFOADsTyEw7jHzTDYE+NmfhaeOmC9mZse1/r37mbNw3qT2PVIItBIrZwLbGta3A5ccrk5EVCXtBRZLGgB+l+wd/+9M8Jgjjb8RuBHgrLPOaqG5OXp6uOKP/x6A33jL6zh/2Ul86IFs/P9rH7yMhXPaJ3fcHC++0s+v3v1PvPv8M/i1n1/BKfOm7thFNzBUY9ueg+w+MEStHjy/+wD7B2q8uHeASjm7uDlpbjtz2rIrh739w7xycIjndx2krVLiDacvYPH8dtrLJboWdDBcC/b1D/FK/zDVWjC/s8Le/mF2vjqYrkIqvPzqEMO1Oqcv7GRguMa+wSrD1TpzOyoMDtcYrNZZcepcliyaQ9eCDtorJfYeHGbhnAqvO2Ue8zsrtKV3hu2VMpC946/VglcHhxmuBacvzK46S3rtarNaq7N/oMq+wWEODtZoq5QYHK4xr6PC/M4K/UM1BBwcqlGPIAL2D1ZHr8T2Dwyzf6hKtZa9m24rl6iUle0neGHvAGWUvWMtZVcHc9vL1OowXK8zp63Mws422iolOipiTluF9oroqJQ5Y1EHp87vpCTRXhalkqjWsne11VqdA0M1BqtVRHZ1MnLFM3IlUq3VGajWqNaCckksmtOGgFrDu/+Ssp9HrZ5dQbRVSgwP1ymXRXu5RH/62Y8QUC6L+e0VOtrKo+URwcHhGtWcq1EJ5rRXsiuWdLzBao29/cP0D9cRcGAw+/lW6/XRd/Pz2isg6FrQwYI0bFlNc3zVWvauvlZ/7TUyNDuipOyKoDxyXkolymVoK5VGr6Sq9aB/uJbOl5jXXub0uZ1j+jAVpnvg9VbgjojYf7ihlPFExF3AXZBdCUzqIOecwzOnPg3Az1x2Eee98TQ+93Pn8fA//ysLLzxnUoc8nCXAI5ddOKXHtEwnsLJh/ednqiFHqS29jvS/dAU4Kb1srI4W6wlo5b1z4/GWTbw5J7RWQmAHh/5clqayvDrb03DQIrIJ4kuAqyR9nOz3uZ6uDh5r4ZjTYuSDY2d3zefGrvnH4j9pZnbcaiUENgIrJa0g+0O9FviPTXV6gOuB7wFXAY9GNtlw+UgFSbeSzQn8WQqK8Y45LZo/OGZmVmTjhkAa478JWE92O+c9EbFZ0m1Ab0T0AHcD90nqA3aT/VGf8DGPsi8taX6EhJlZkbU0JxAR64B1TWW3NCwPAFePc4xbxzvmseAQMDN7TTE+MdxgshPUZmazUQFDYKZbYGZ2/CheCMx0A8zMjiOFCYGRK4DJPMLAzGy2KkwInHXKXABOW9Dqx0zMzGa/wjyq8Us3XMLXfriDFadO7tkbZmazUWFCYNkpc/mvV6wcv6KZWYEUZjjIzMzGcgiYmRWYQ8DMrMAcAmZmBeYQMDMrMIeAmVmBOQTMzArMIWBmVmDKvgDsxCBpJ/D8JHc/FXh5CptzInCfi6FofS5af+Ho+/y6iOjK23BChcDRkNQbEd0z3Y5jyX0uhqL1uWj9hents4eDzMwKzCFgZlZgRQqBu2a6ATPAfS6GovW5aP2FaexzYeYEzMxsrCJdCZiZWROHgJlZgc36EJC0WtIWSX2Sbp7p9hwNScskfUvSU5I2S/rtVH6KpL+X9JP078mpXJI+lfr+I0kXNhzr+lT/J5Kun6k+tUpSWdIPJf1NWl8haUPq24OS2lN5R1rvS9uXNxzjo6l8i6R3zVBXWiLpJEkPSXpa0j9Lumy2n2dJ/y39Xj8p6SuSOmfbeZZ0j6SXJD3ZUDZl51XSRZKeSPt8Shr5dvUjiIhZ+wLKwDPA2UA78DiwaqbbdRT9WQJcmJYXAD8GVgEfB25O5TcDt6flXwK+CQi4FNiQyk8BtqZ/T07LJ890/8bp+4eB+4G/Set/AaxNy58BfistfxD4TFpeCzyYllel898BrEi/F+WZ7tcR+nsv8J/Scjtw0mw+z8CZwLPAnIbz+2uz7TwDbwMuBJ5sKJuy8wr8U6qrtO+V47Zppn8o0/wDvwxY37D+UeCjM92uKezf14F3AFuAJalsCbAlLX8WuLah/pa0/Vrgsw3lh9Q73l7AUuAR4BeBv0m/4C8DlebzDKwHLkvLlVRPzee+sd7x9gIWpT+Iaiqftec5hcC29Ietks7zu2bjeQaWN4XAlJzXtO3phvJD6h3uNduHg0Z+sUZsT2UnvHT5+2ZgA3B6RLyYNv0UOD0tH67/J9rP5U+BjwD1tL4YeCUiqmm9sf2jfUvb96b6J1KfVwA7gc+nIbA/lzSPWXyeI2IH8CfAvwAvkp23x5jd53nEVJ3XM9Nyc/kRzfYQmJUkzQe+CnwoIvY1bovsLcCsue9X0nuAlyLisZluyzFUIRsy+HREvBk4QDZMMGoWnueTgTVkAfgzwDxg9Yw2agbMxHmd7SGwA1jWsL40lZ2wJLWRBcCXI+IvU/G/SlqSti8BXkrlh+v/ifRzeQvwHyQ9BzxANiT0SeAkSZVUp7H9o31L2xcBuzix+rwd2B4RG9L6Q2ShMJvP878Hno2InRExDPwl2bmfzed5xFSd1x1pubn8iGZ7CGwEVqY7DNrJJpB6ZrhNk5Zm+u8G/jkiPtGwqQcYuUPgerK5gpHy96e7DC4F9qbLzvXAOyWdnN6BvTOVHXci4qMRsTQilpOdv0cj4jrgW8BVqVpzn0d+Flel+pHK16a7SlYAK8km0Y47EfFTYJukN6SiK4CnmMXnmWwY6FJJc9Pv+UifZ+15bjAl5zVt2yfp0vQzfH/DsQ5vpidJjsEkzC+R3UXzDPD7M92eo+zLW8kuFX8EbEqvXyIbC30E+AnwMHBKqi/gztT3J4DuhmP9BtCXXr8+031rsf+/wGt3B51N9j93H/B/gI5U3pnW+9L2sxv2//30s9hCC3dNzHBfLwB607n+K7K7QGb1eQb+EHgaeBK4j+wOn1l1noGvkM15DJNd8d0wlecV6E4/v2eAP6Pp5oK8lx8bYWZWYLN9OMjMzI7AIWBmVmAOATOzAnMImJkVmEPAzKzAHAJmZgXmEDAzK7D/D3gPDO09xtbRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.lineplot(data = convergeY)\n",
    "sns.lineplot(data = np.ones(n)*expectedY,color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "id": "d16c0098-95dd-4771-83d9-dc0bb67db08b",
   "metadata": {},
   "outputs": [],
   "source": [
    "optionPrice = list()\n",
    "\n",
    "for i in range(len(Spaths)):\n",
    "    \n",
    "    if (max(Spaths[i]) > A) or (min(Spaths[i]) < B):\n",
    "        optionPrice.append(0)\n",
    "    else:\n",
    "        optionPrice.append(max(Spaths[i][-1] - K,0))\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "id": "abdcbe10-7930-4bb4-b6fa-2acbeda5752b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0311188574285244"
      ]
     },
     "execution_count": 303,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optionPriceInitial = sum(optionPrice)/len(optionPrice)\n",
    "optionPriceInitial"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}