Two_Factor_Model-SimpleVersion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Two-Factor Forward Price Model for Energy Commodities: the maths, calibration, implementation and applications to energy transport contracts\n",
    "* * * * * "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##About this notebook\n",
    "Documention prepared by **Jesus Perez Colino**.\n",
    "Version 0.1, Released 25/11/2014, Alpha\n",
    "\n",
    "\n",
    "- This work is licensed under a [Creative Commons Attribution-ShareAlike 3.0 Unported License](http://creativecommons.org/licenses/by-sa/3.0/deed.en_US). This work is offered for free, with the hope that it will be useful.\n",
    "\n",
    "\n",
    "- **Summary**: This technical note has as main goal not only to show the details of the **Clewlow and Strickland (1999)** multivariate model implementation in Python, but also propose an implementation of an adapted version for tradeable products in Python, that is consistent with Black-76 model. \n",
    "\n",
    "\n",
    "- **Index**:\n",
    "\n",
    "    0. About this notebook\n",
    "    1. Basics about the Two-Factor Forward Model\n",
    "    2. Monte-Carlo Simulation of the Forward Curve based in a Two-Factor Model\n",
    "    3. The Volatility Term-Structure and the Equivalent Volatility\n",
    "    4. Fitting the Model to the Market: a Convex Calibration procedure\n",
    "    5. Log-Normal approximation for the spread option valuation using the Two-Factor model\n",
    "\n",
    "\n",
    "- **Reference**: [Clewlow L. and C. Strickland, 1999b, \"A Multi-Factor Model for Energy Derivatives and Risk Management\", Working Paper, University of Tecnology, Sydney.](https://www.google.de/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&sqi=2&ved=0CCEQFjAAahUKEwjqgbyF0LnHAhVEWBQKHftCBvQ&url=http%3A%2F%2Fwww.erasmusenergy.com%2Fdownloadattachment.php%3FaId%3Db68bca6b17e223fab5bfdaff87404006%26articleId%3D103&ei=kM_WVer5H8SwUfuFmaAP&usg=AFQjCNFYYDeAudvzk_PhdZ9Fm-woZJBUpg&sig2=7yTr3iKyMJI8b-J1SJ17sA)\n",
    "\n",
    "\n",
    "- **Reproducibility conditions**: the content in this technical note can be reproduce in your computer under the following versions of the Python's packages: \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------------- Reproducibility conditions for the Two-factor model--------------\n",
      " \n",
      " Python Version: 2.7.10 |Anaconda 2.3.0 (x86_64)| (default, May 28 2015, 17:04:42) \n",
      "[GCC 4.2.1 (Apple Inc. build 5577)]\n",
      " Numpy Version: 1.9.2\n",
      " \n",
      "--------------------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import datetime, timedelta\n",
    "import scipy.stats as st\n",
    "\n",
    "print ' Reproducibility conditions for the Two-factor model'.center(80, '-')\n",
    "print ' '\n",
    "print ' Python Version: ' + sys.version\n",
    "print ' Numpy Version: ' + np.__version__\n",
    "print ' ' \n",
    "print '-'*80"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basics about Two-Factor Forward-Price Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The 'Instantaneous' Forward-Prices dynamic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "All over this technical note, we will consider a two-dimensional Brownian motion $B_{S,t},B_{L,t}$ with correlation matrix $\\rho$ on the probability space $(\\Omega,\\mathcal{F},\\mathbb{P} )$ endowed with the natural (completed and right-continuous) filtration $\\mathbb{F}=(\\mathcal{F_t})_{t \\geq T}$ generated by ${B_{S,t},B_{L,t}}$ up to some fixed time horizon $T>0$. \n",
    "\n",
    "The main goal of this section is to define the dynamic of a commodity **forward curve model** based in the following stochastic differential equation (SDE), where we suppose that the uncertainties in the forward price will be explained by the defined Brownian motions ${B_{S,t},B_{L,t}}$ (**Clewlow and Strickland, 1999**): \n",
    "\n",
    "\\begin{equation} \n",
    "\\frac{dF(t,T)}{F(t,T)}=\\sigma_{S}(t,T)e^{-a(T-t)}dB_{S,t}+\\sigma_{L}(t,T)dB_{L,t}\n",
    "\\end{equation}\n",
    "\n",
    "where the following are used:\n",
    "\n",
    "* $ F(t,T) $: forward price at time t for delivery of one quantity of the underlying commodity at time T\n",
    "* $ \\sigma_{S}(t,T) $: short-term volatility,\n",
    "* $ \\sigma_{L}(t,T) $: long-term volatility,\n",
    "* $ a $: a mean-reversion speed\n",
    "* $ dB_{S,t} $ and $ dB_{L,t} $: brownian motions with correlation \u001a$\\rho = \\rho(dB_{S,t}, dB_{L,t})$ \n",
    "\n",
    "The two-factors model here introduced assumes the dynamic of forward market prices is described by two distinct sources of risk (random shocks here introduced as gaussian processes or Brownian motions). In this model, the first Brownian motion describes swift changes in the future curve, the second one describes structural changes in the gas market and deals with long term changes in the curve. Due to the decay effect given by the exponential term, the first risk factor affects particularly the behaviour of prices with short-time to delivery (front-end of the forward curve) and hence is usually referred to as **”short-term factor”** and denoted by the sub-index $S$. The second risk factor has the greatest influence on the behaviour of the back-end of the curve (long time-to-delivery) and will be referred to as **”long-term factor”** denoted by the sub-index $L$. \n",
    "\n",
    "Therefore, the infinitesimal behaviour of the Forward Prices in time $t$,for delivery at time $T$, is clearly random and driven by the following constituents:\n",
    "* the sum of two gaussian process, possibly correlated, here represented by the Brownian motions $B_L$ and $B_S$\n",
    "* the value in $t$ of the forward price, or $F(t,T)$ \n",
    "\n",
    "Another way of viewing the two factors comes from the economic interpretation: The first factor captures the increased trading activity as knowledge about weather, unexpected outages etc. becomes available. The second factor models a long-term uncertainty, that is common to all products in the market. This uncertainty can be explained by technological advances, political changes, price developments in other commodity markets and many more.\n",
    "\n",
    "The main reasons because this is the most extended model in the commodities market to price derivatives and structured products is, not only his **simplicity** and the **flexibility** to represent volatility term structures with exponential-decay pattern, but most importantly, this model is able to **aggregate** the position and hedging of complex derivatives' portfolio with many different maturities in *the two most liquid forward contracts*: calendar contracts representing the ”long-term factor”, and front-month (or weekahead) for the \"short-term\" factor. Or in other words, *the factors are not only tradeable, but the most liquid tradeble products in the market*. \n",
    "\n",
    "In order to get the solution of the previous SDE, we just need to apply Itô’s formula (Kiyosi Itô, 1950) to the $\\log F(t,T)$ taking into account the previous SDE, and integrate the expression to obtain the forward prices:\n",
    "\n",
    "\\begin{equation}\n",
    "d\\log{F(t,T)} = \\frac{dF(t,T)}{F(t,T)} - \\frac{1}{2}\\left(\\sigma_{S}e^{-2a(T-t)}+\\sigma_{L}+2\\rho\\sigma_{L}\\sigma_{S}e^{-a(T-t)}\\right)dt\n",
    "\\end{equation}\n",
    "\n",
    "Integrating the previous expression and using the following definitions:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{lll}\n",
    "V(t_0,t,T) & =& \\int_{t_0}^{t}\\sigma_{S}^2e^{-2a(T-u)}du+\\int_{t_0}^t\\sigma^2_Ldu+\\int_{t_0}^t2\\rho\\sigma_L\\sigma_Se^{-a(T-u)}du\\\\\n",
    "W_S(t_0,t,T)& =& \\int_{t_0}^{t}\\sigma_{S}e^{au}dB_{S,t}\\\\\n",
    "W_L(t_0,t,T)& =& \\int_{t_0}^{t}\\sigma_{L}dB_{L,t}\\\\\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "we have a much easier expression to understand of the solution of the original SDE: \n",
    "\n",
    "\\begin{equation}\n",
    "F(t,T) = F(t_0,T)e^{-\\frac{1}{2}V(t_0,t,T)+e^{-aT}W_S(t_0,t,T)+W_L(t_0,t,T)}\n",
    "\\end{equation}\n",
    "\n",
    "The SDE that define the dynamics of the spot price can be obtained directly from the previous equation, simply defining the spot price $S_t=F(t,t)$ as the price at time $t$ for instant delivery, such as\n",
    "\n",
    "\\begin{equation}\n",
    "S_t = F(t,t)e^{-\\frac{1}{2}V(t_0,t,t)+e^{-at}W_S(t_0,t,t)+W_L(t_0,t,t)}\n",
    "\\end{equation}\n",
    "\n",
    "and, following similar procedure than in **Clewlow and Strickland (1999)**, we obtain the SDE for the dynamic of the spot prices:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d S_t}{S_t} = \\frac{\\partial lnF(t_0,t)}{\\partial t} - \\frac{1}{2} \\frac{\\partial V(t_0,t,t)}{\\partial t}\n",
    "-a\\left\\{lnS_t-lnF(t_0,t)+\\frac{1}{2}V(t_0,t,t)-\\sigma_L(B_{L,t}-B_{L,t_0} )\\right\\}+\\sigma_SdB_{S,t}+\\sigma_LdB_{S,t}\n",
    "\\end{equation}\n",
    "\n",
    "However, the special properties of gas and power markets, especially that of the limited or total non-storability of the underlying commodity, lead to non-validity of the spot-forward relationship here showed. Traditional theory utilises no-arbitrage and cash-and-carry arguments to derive this well-known spot-forward relationship, however the limited (gas) or total non-storability (power) of the underlying commodity implies the break down of traditional cash-and-carry arguments as well as the well-known spot- forward relationship. Additionally, any estimation of the parameters of the SDE for any forward prices model using statistical information from the spot market is wrong."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 'Tradeable' Forward-prices dynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the mathematical development written above is only valid for intantaneous forward prices, or in other words, forward prices for instantaneous product delivery. In reality, there is no such thing in comodities markets. Forwards contracts in commodities are for a specified delivery period; days, weeks, months, seasons or years. The simulations produced for products with long delivery periods will have volatilies that are too high in reality, because the actual price path within the product is in effect averaged to give the actual product price, and this will affect the volatility of the product itself. It is quite intuitive to expect the volatility of yearly products to be much less than daily products, for example.\n",
    "\n",
    "In this subsection we explain how the instanateous model can be adjusted to describe with accuracy the dynamics of tradeable forward products with longer delivery periods. Given an instanateous forward-price curve (e.g. daily for gas or hourly curve for power), from a simple no arbitrage opportunity argument, it is possible to derive the price for any kind of forward product specified by a delivery time $T$ and a delivery length $\\theta$ where $F_p(t,T,\\theta)$ is defined as\n",
    "\n",
    "\\begin{equation}\n",
    "F_p(t,T,\\theta) = \\frac{1}{\\theta}\\int_0^{\\theta} F(t,T + u)du\n",
    "\\end{equation}\n",
    "\n",
    "It is not unusual to find that most of the practicioners in the commodities markets are achiving this relationship estimating statisitically some scaling factors usually defined as **shaping factors** and most commonly defined as \n",
    "\n",
    "\\begin{equation}\n",
    "\\lambda_{t,T,u} = \\frac{F(t,T+u)}{F_p(t,T,\\theta)}  \\forall u \\in [0,\\theta)\n",
    "\\end{equation}\n",
    "\n",
    "Each of these shaping factors can be interpreted as a weight of the hour or day $u \\in [0, \\theta]$ compared with the tradeable forward price of a contract that has 'continuous' delivery inside the interval $[T, T + \\theta]$. They are intuitive in the energy markets as they can be related to the various seasonalities observable in the market (for instance, forward price for delivery on a Monday at 8am is more expensive than for delivery on Sunday at midnight for instance).\n",
    "\n",
    "Taking into account no-arbitrage conditions, it is straightforward to show that the average of the shaping factors over the delivery period is equals to:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{1}{\\theta}\\int_0^{\\theta}\\lambda_{t,T,u} du=\\frac{1}{\\theta}\\sum_{u=1}^{\\theta}\\lambda_{t,T,u} = 1\n",
    "\\end{equation}\n",
    "\n",
    "Therefore, we can re-define the previous relationship between instantaneous and tradeable forward prices, including the shape factors, as:\n",
    "\n",
    "\\begin{equation}\n",
    "F_p(t,T,\\theta) = \\frac{1}{\\theta}\\int_0^{\\theta} \\lambda_{t,T,u}F(t,T+u)du\n",
    "\\end{equation}\n",
    "\n",
    "Using the results of the previous section, we will get the expression:\n",
    "\n",
    "\\begin{equation}\n",
    "dF_p(t,T,\\theta) = \\frac{F_p(t,T,\\theta)}{\\theta} \\int_0^{\\theta}\\lambda_{t,T,u} \\left[\\sigma_S e^{-a(T+u-t)}dB_{S,t}+\\sigma_LdB_{L,t} \\right]du  \\forall u \\in [0,\\theta)\n",
    "\\end{equation}\n",
    "\n",
    "and, assuming that shaping factors are constant throught time, we can define \n",
    "\n",
    "\\begin{equation}\n",
    "\\Psi_{t,T,\\theta_u}=\\frac{1}{\\theta}\\int_0^{\\theta}\\lambda_{t,T,u}e^{-au}du\n",
    "\\end{equation}\n",
    "\n",
    "and obtain an SDE similar to the used for the instantaneous forward prices in the previous section:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dF_p(t,T,\\theta)}{F_p(t,T,\\theta)}= \\sigma_S \\Psi_{t,T,\\theta_u} e^{-a(T-t)}dB_{S,t}+\\sigma_LdB_{L,t} \n",
    "\\end{equation}\n",
    "\n",
    "**As a particular case**, and just for the sake of simplicity, we will assume for the rest of this technical note that $\\lambda_{t,T,\\theta_u} = 1$, $\\forall u \\in [0,\\theta)$ which leads to \n",
    "\n",
    "\\begin{equation}\n",
    "\\Psi_{t,T,\\theta_u}=\\Psi_{\\theta}=\\frac{1-e^{-a\\theta}}{a\\theta}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def PsiFunction(meanReverting, \n",
    "                deliveryStart, deliveryEnd, \n",
    "                blockSize  = 0, granularity = 1 / 365.):\n",
    "    \n",
    "    \"\"\"\n",
    "    Correction term to apply on short-term volatility \n",
    "    when instaneous model is adapted for higher granularity products.\n",
    "    \"\"\"\n",
    "    try: \n",
    "        meanReverting = float(meanReverting)\n",
    "        if meanReverting < 0 or deliveryEnd < deliveryStart:\n",
    "            raise ValueError('Negative inputs.')\n",
    "    except: \n",
    "        print('Error passing vol inputs to PsiFunction')\n",
    "        \n",
    "    if (blockSize == 0):\n",
    "        blockSize = (deliveryEnd - deliveryStart) / 365.\n",
    "    if (blockSize == granularity) or (meanReverting == 0):\n",
    "        PsiFunc = 1\n",
    "    else:\n",
    "        PsiFunc = (1 - np.exp(- meanReverting * blockSize)) / float(meanReverting * blockSize)\n",
    "        \n",
    "    return PsiFunc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be focus in the implementation of this version of the model. Simulating forward products with any delivery length is simply done using the techniques that will show in the next sections after adjusting the short-term volatility by the appropriate factor \t$\\Psi_{t,T,\\theta_u}$. Therefore, the model that we are interested now reads as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dF_p(t,T,\\theta)}{F_p(t,T,\\theta)}= \\sigma_S \\Psi_{\\theta} e^{-a(T-t)}dB_{S,t}+\\sigma_LdB_{L,t} \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And using the same argument that in the previous section, we have\n",
    "\n",
    "\\begin{equation}\n",
    "F_p(t,T,\\theta) = F_p(t_0,T,\\theta)e^{-\\frac{1}{2}V(t_0,t,T)+e^{-aT}W_S(t_0,t,T)+W_L(t_0,t,T)}\n",
    "\\end{equation}\n",
    "\n",
    "where $V(t_0,t,T)= V_{S}(t_0,t) + V_{L}(t_0,t) + V_{SL}(t_0,t)$ and each component is defined as \n",
    "\n",
    "1.- \n",
    "\n",
    "\\begin{equation}\n",
    "V_{S}(t_0,t)= \\int_{t_0}^t \\sigma_{S}^2 \\Psi_{\\theta}^2 e^{-2a(T-u)}du = \\frac{\\sigma_{S}^2 \\Psi_{\\theta}^2}{2a} [e^{-2a(T-t)}-e^{-2a(T-t_0)}]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def V_S(STvol, meanReverting, \n",
    "        valuationdate , maturitydate, \n",
    "        deliveryStart, deliveryEnd,\n",
    "        blockSize = 0):\n",
    "    try: \n",
    "        STvol = float(STvol)\n",
    "        meanReverting = float(meanReverting)                                                                                            \n",
    "        if STvol < 0 or meanReverting < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if valuationdate > maturitydate or deliveryStart > deliveryEnd:\n",
    "            raise ValueError('Wrong dates')\n",
    "    except: \n",
    "        print('Error passing inputs to V_S')\n",
    "\n",
    "    TimeToMaturity = (maturitydate - valuationdate) / 365.\n",
    "    psi = PsiFunction(meanReverting, deliveryStart, deliveryEnd, blockSize)\n",
    "        \n",
    "    if meanReverting == 0:\n",
    "        V_S = STvol * psi * np.sqrt(TimeToMaturity)\n",
    "    else:\n",
    "        V_S = np.exp(-2. * meanReverting * (deliveryStart - maturitydate) / 365.) \\\n",
    "            - np.exp(-2. * meanReverting * (deliveryStart - valuationdate) / 365.)\n",
    "        V_S = (STvol ** 2 * psi ** 2 * V_S) / float(2. * meanReverting)\n",
    "\n",
    "    return V_S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.- \n",
    "\n",
    "\\begin{equation}\n",
    "V_{L}(t_0,t)= \\int_{t_0}^t\\sigma^2_Ldu = \\sigma_{L}^2 (t-t_0)\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def V_L(LTvol, valuationdate, maturitydate):\n",
    "    try: \n",
    "        LTvol = float(LTvol)                                                                                           \n",
    "        if LTvol < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if valuationdate > maturitydate:\n",
    "            raise ValueError('Wrong dates')\n",
    "    except: \n",
    "        print('Error passing inputs to V_L')\n",
    "        \n",
    "    V_L = LTvol * ((maturitydate - valuationdate) / 365.)\n",
    "    return V_L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3.- \n",
    "\n",
    "\\begin{equation}\n",
    "V_{SL}(t_0,t)= \\int_{t_0}^t2\\rho\\sigma_L\\sigma_S\\Psi_{\\theta}e^{-a(T-u)}du =  \\frac{2\\rho\\sigma_L\\sigma_S\\Psi_{\\theta}}{a}[e^{-a(T-t)}-e^{-a(T-t_0)}]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def V_SL(LTvol, STvol, meanReverting,\n",
    "         valuationdate, maturitydate, \n",
    "         deliveryStart, deliveryEnd, \n",
    "         corr, blockSize = 0, granularity = 1/365.): \n",
    "\n",
    "    try: \n",
    "        LTvol = float(LTvol)\n",
    "        STvol = float(STvol)\n",
    "        meanReverting = float (meanReverting)\n",
    "        corr = float(corr)\n",
    "        if LTvol < 0 or STvol < 0 or meanReverting < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if deliveryStart > deliveryEnd or valuationdate > maturitydate :\n",
    "            raise ValueError('Negative time periods')\n",
    "        if corr > 1 or corr < -1 :\n",
    "            raise ValueError('Incorrect Correlation')\n",
    "    except: \n",
    "        print('Error passing vol inputs to V_SL')\n",
    "        \n",
    "    TimeToMaturity = (maturitydate - valuationdate) / 365.\n",
    "    psi = PsiFunction(meanReverting, deliveryStart, deliveryEnd, blockSize, granularity)\n",
    "    if meanReverting == 0 :\n",
    "        V_SL = 0\n",
    "    else:\n",
    "        V_SL = (np.exp(-meanReverting * (deliveryStart - maturitydate) / 365.)\n",
    "                - np.exp(-meanReverting * (deliveryStart - valuationdate) / 365.))\n",
    "        V_SL = (2.* corr * LTvol * STvol * psi  * V_SL / float(meanReverting))     \n",
    "        #V_SL = (2.* corr * LTvol * STvol * psi  * V_SL / float(meanReverting * TimeToMaturity)) \n",
    "                                                               #<--*TimetoMaturity?????\n",
    "        \n",
    "    return V_SL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the stochastic noise is defined by these two stochastic integrals, quite similar as the ones defined in the previous model, \n",
    "\n",
    "\\begin{equation} \n",
    "\\begin{array}{lll}\n",
    "W_S(t_0,t,T)& = & \\int_{t_0}^{t}\\sigma_{S}e^{au}\\Psi_{\\theta}dB_{S,t}\\\\\n",
    "W_L(t_0,t,T)& = & \\int_{t_0}^{t}\\sigma_{L}dB_{L,t}\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "that can be treated as stochastic variables.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte-Carlo Price Simulation with the Two-Factor Forward Price Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following section describes how the model can be used to generate random samples of forward prices using Monte Carlo techniques. Here, we are particularly interested in the generation of samples of prices for an specific tradeable forward contract, in a particular point of time $t+\\delta$ where $\\delta > 0$, in the following sense: \n",
    "\n",
    "\\begin{equation}\n",
    "F_p(t+\\delta,T,\\theta)=F_p(t,T,\\theta)e^{-\\frac{1}{2}V(t,t+\\delta,T)+e^{-aT}W_{S,t+\\delta}+W_{L,t+\\delta}}\n",
    "\\end{equation}\n",
    "\n",
    "Notice that, according with the notation showed in previous sections, we defined:\n",
    "\n",
    "\\begin{equation}\n",
    "V(t,t+\\delta,T) =  V_{S}(t,t+\\delta)+V_{L}(t,t+\\delta)+V_{SL}(t,t+\\delta)\n",
    "\\end{equation}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{lllll}\n",
    "V_{S}(t,t+\\delta) & = & \\int_{t}^{t+\\delta}\\Psi_{\\theta}^2 \\sigma_{S}^2 e^{2au}du &= & e^{-2a(T-t)} \\sigma_{S}^2 \\Psi_{\\theta}^2 \\frac{(e^{2a\\delta}-1)}{2a} \\\\\n",
    "V_{L}(t,t+\\delta) & = & \\int_{t}^{t+\\delta}\\sigma^2_Ldu &=& \\sigma_{L}^2 \\delta \\\\\n",
    "V_{SL}(t,t+\\delta) & = & \\int_{t}^{t+\\delta}2\\Psi_{\\theta}\\rho\\sigma_L\\sigma_Se^{au}du & = & e^{-a(T-t)}2\\rho\\sigma_L\\sigma_S\\Psi_{\\theta}\\frac{(e^{a\\delta}-1)}{a} \\\\\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "and, therefore, we obtain finally\n",
    " \n",
    "\\begin{equation}\n",
    "V(t,t+\\delta,T) = e^{-2a(T-t)} \\sigma_{S}^2 \\Psi_{\\theta}^2 \\frac{(e^{2a\\delta}-1)}{2a} + \\sigma_{L}^2 \\delta + e^{-a(T-t)}2\\rho\\sigma_L\\sigma_S\\Psi_{\\theta}\\frac{(e^{a\\delta}-1)}{a}\n",
    "\\end{equation}\n",
    "\n",
    "and notice that, in the particular case where we are interested in the forward price simulation at maturity (e.g. european options whith forward contract as underlying ), the last expression will become \n",
    "\n",
    "\\begin{equation}\n",
    "V(t,t+\\delta,t+\\delta) = \\Psi_{\\theta}^2\\sigma_{S}^2 \\frac{1-e^{-2a(\\delta)}}{2a} + \\sigma_{L}^2 \\delta + 2\\Psi_{\\theta}\\rho\\sigma_L\\sigma_S \\frac{1-e^{-a\\delta}}{a}\n",
    "\\end{equation}\n",
    "\n",
    "See here below the Python implementation of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def V_Forward(LTvol, STvol, meanReverting, \n",
    "              valuationdate, maturitydate, \n",
    "              deliveryStart, deliveryEnd, \n",
    "              corr = 0):\n",
    "    \"\"\"\n",
    "    Compute the drift term resulting from Ito in the forward price formula\n",
    "    for a tradeable product assuming a 2 factors model over forward-prices\n",
    "    --> Assumming that t + delta = T = maturitydate\n",
    "    \"\"\"\n",
    "    try: \n",
    "        LTvol = float(LTvol)\n",
    "        STvol = float(STvol)\n",
    "        if meanReverting != 0:\n",
    "            meanReverting = float(meanReverting)\n",
    "        if corr != 0 :\n",
    "            corr = float(corr)                                                                                                \n",
    "        if LTvol < 0 or STvol < 0 or meanReverting < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if deliveryStart > deliveryEnd or valuationdate > maturitydate :\n",
    "            raise ValueError('Incorrect dates')\n",
    "        if corr > 1 or corr < -1 :\n",
    "            raise ValueError('Wrong correlation')\n",
    "    except: \n",
    "        print('Error passing inputs to V_Forward')\n",
    "    \n",
    "    delta_time = (maturitydate - valuationdate) / 365.\n",
    "    \n",
    "    psi = PsiFunction(meanReverting, \n",
    "                      deliveryStart, deliveryEnd, \n",
    "                      blockSize = 0, granularity = 1 / 365.)\n",
    "\n",
    "    Vs = V_S(STvol, meanReverting, \n",
    "             valuationdate , maturitydate, \n",
    "             deliveryStart, deliveryEnd, \n",
    "             blockSize = 0)\n",
    "    \n",
    "    Vl = V_L(LTvol, valuationdate, maturitydate)\n",
    "    \n",
    "    Vsl = V_SL(LTvol, STvol, meanReverting,\n",
    "                valuationdate, maturitydate, \n",
    "                deliveryStart, deliveryEnd, \n",
    "                corr, blockSize = 0, granularity = 1/365.)\n",
    "\n",
    "    if meanReverting == 0:\n",
    "        V_Forward = V_L\n",
    "    else:\n",
    "        V_Forward = Vs + Vl + Vsl\n",
    "        \n",
    "    return V_Forward"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and, additionally, notice that using Ito isometry property to sample $W_S$ and $W_L$ \n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{llll}\n",
    "W_S(t,t+\\delta) & =  \\int_{t}^{t+\\delta}\\sigma_{S}e^{au}dB_{S,u} & \\rightarrow & N ( 0, e^{2at}\\sigma_S^2\\frac{e^{2a\\delta}-1}{2a}) \\\\\n",
    "W_L(t,t+\\delta) & =  \\int_{t}^{t+\\delta}\\sigma_{L}dB_{L,u} & \\rightarrow & N ( 0, \\sigma_L^2\\delta) \n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "will allow us to obtain the initial expression, that can be easily implemented for simulation where $\\delta = $ time to maturity in daily terms over a year. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def GenerateBrownian(corr, size = 1, seed = None ):\n",
    "    try:\n",
    "        if corr != 0 :\n",
    "            corr = float(corr)     \n",
    "        if size <= 0:\n",
    "            raise ValueError('Wrong sample size')\n",
    "        else:\n",
    "            assert isinstance(size, int)\n",
    "        if seed is not None:\n",
    "            assert isinstance(seed, int)\n",
    "    except:\n",
    "        print 'Error passing inputs to GenerateBrownian' \n",
    "    np.random.seed(seed)\n",
    "    Bs = np.random.randn(size)\n",
    "    Bl = np.random.randn(size)\n",
    "    Bl = corr * Bs + np.sqrt(1 - corr ** 2) * Bl\n",
    "    \n",
    "    return Bs, Bl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def GenerateTwoFactors(STvol, LTvol, meanReverting, corr,\n",
    "                       maturitydate, valuationdate,\n",
    "                       deliveryStart, deliveryEnd,\n",
    "                       size = 1, seed = None):\n",
    "    try:\n",
    "        if STvol < 0 or LTvol < 0 or corr < -1 or corr > 1 or meanReverting < 0 : \n",
    "            raise ValueError('Wrong vol/corr/meanrev. inputs')\n",
    "        if corr != 0 :\n",
    "            corr = float(corr)     \n",
    "        if size <= 0 :\n",
    "            raise ValueError('Wrong sample size')\n",
    "        if valuationdate > maturitydate or deliveryStart > deliveryEnd: \n",
    "            raise ValueError('Wrong dates')\n",
    "        else :\n",
    "            assert isinstance(size, int)\n",
    "    except:\n",
    "        print 'Error passing inputs to GenerateTwoFactors' \n",
    "    \n",
    "    delta_time = (maturitydate - valuationdate) / 365.\n",
    "    \n",
    "    [Bs, Bl] = GenerateBrownian(corr, size, seed)\n",
    "    psi = PsiFunction(meanReverting, deliveryStart, deliveryEnd)\n",
    "    \n",
    "    Ws = Bs * np.sqrt (np.exp(2 * meanReverting * valuationdate / 365.) * psi**2 * STvol**2 \\\n",
    "                    * (np.exp(2 * meanReverting * delta_time)-1) \\\n",
    "                    / 2 * float(meanReverting)) * np.exp(-meanReverting*maturitydate/365.)\n",
    "    \n",
    "    Wl = Bl * LTvol * np.sqrt(delta_time)\n",
    "    \n",
    "    return Ws, Wl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def MCSimulation_2F(ForwardPrice, \n",
    "                    STvol, LTvol, meanReverting, corr, \n",
    "                    maturitydate, valuationdate, \n",
    "                    deliveryStart, deliveryEnd, \n",
    "                    size = 1, seed = None):\n",
    "    try: \n",
    "        LTvol = float(LTvol)\n",
    "        STvol = float(STvol)\n",
    "        if meanReverting != 0:\n",
    "            meanReverting = float(meanReverting)\n",
    "        if corr != 0 :\n",
    "            corr = float(corr)                                                                                                \n",
    "        if LTvol < 0 or STvol < 0 or meanReverting < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if deliveryStart > deliveryEnd or valuationdate > maturitydate :\n",
    "            raise ValueError('Incorrect dates')\n",
    "        if corr > 1 or corr < -1 :\n",
    "            raise ValueError('Wrong correlation')\n",
    "    except: \n",
    "        print('Error passing inputs to MCSimulation_2F')\n",
    "    \n",
    "    V = V_Forward(LTvol, STvol, meanReverting, \n",
    "                  valuationdate, maturitydate, \n",
    "                  deliveryStart, deliveryEnd, corr)\n",
    "    \n",
    "    [Ws, Wl] = GenerateTwoFactors(STvol, LTvol, meanReverting, corr, \n",
    "                                  maturitydate, valuationdate, \n",
    "                                  deliveryStart, deliveryEnd, \n",
    "                                  size, seed)\n",
    "    \n",
    "    F_simul = ForwardPrice * np.exp(-.5*V + Ws + Wl)\n",
    "    \n",
    "    return F_simul"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+oAAAI6CAYAAACw17JJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VOXd//HPTEJWAiHBECACYTeBChgQi2yBgAoWaikg\nWNmUFhCppbbgvjS2PohSKNY+BQx1q1QrtfK4BGQVZU1AAgoBAZGESEIgErLN3L8//DEykzAZIBkO\n5P26Li5zznzPfb5ncieXn5xlbMYYIwAAAAAAYAn2y90AAAAAAAD4AUEdAAAAAAALIagDAAAAAGAh\nBHUAAAAAACyEoA4AAAAAgIUQ1AEAAAAAsBCCOgDUsPj4eD3zzDOXu40rSkVFhSZOnKjGjRvLbrdr\n3bp1l7ulq9748eOVkpLil3098cQTateuXY2MZbfb9frrr9fIWDVh9uzZatKkiex2u/7xj39c7nYs\nr7a+fwcPHpTdbtfGjRtrfGwAuBwI6gDgA2+hxvN/PLdu3aoHHnjAp3E3bNggu92uw4cP10ifV6q3\n335bb7zxht577z3l5ubqpptuqrIuNTVVvXv3VoMGDWS323X06FGfxt+7d68GDx6s8PBwXXPNNZoy\nZYqKi4vdanJycjRy5Eg1bNhQDRs21J133qlvv/3W9frZIFDVv+eff/6ij72oqEj33nuvGjdurPr1\n6+u2227TgQMHKtWtW7dOycnJioiIUEREhHr06KGDBw9e9H4XLFigt95666K3v1A2m+2C6gcOHKgJ\nEyZUWp+bm6uf/exnNdXWJdm0aZOeffZZLV68WLm5uRo5cmSNjf3EE09UOdf+53/+p8b2YUVpaWlu\nxxsbG6vbb79du3bt8rpdixYtlJubqx49evipUwCoXYGXuwEAuBLYbDafg0Z0dPQFj2+MueBtLkZZ\nWZmCgoL8sq8LsW/fPjVv3lw9e/b0WldWVqbhw4fr9ttv16xZs3wa+7vvvtOAAQPUpUsXffrpp8rP\nz9fEiRNVWFioN954Q5LkdDo1dOhQBQYGauXKlXI6nZo6daqGDx+uTz75xG28d999t1IYiIiIqHLf\nBw8eVOvWreV0Os/b3y9+8Qvt2rVLb7/9tho2bKjZs2crJSVFWVlZCgkJkSR98MEHGjFihGbNmqX5\n8+crLCxMe/bsUVhYmE/vQVXO13Ntqak5HhMTUyPj1IR9+/bJbrdr6NChtTJ+fHy8Pv30U7d19evX\nv+jxKioqFBhY+//rd6m/ZwICAvTNN99Ikg4cOKAZM2bolltu0Z49e6qct+Xl5apXr56l5gYAXDID\nAKjWuHHjTEpKSpWv2Ww289prr7mWW7Zsaf7whz+4lpcvX266dOliwsLCTGRkpOnRo4fJyMgwX331\nlbHZbG7/+vfvb4wxxul0mjlz5pj4+HgTFBRk2rRpY+bNm+e23+PHj5sRI0aY8PBwExsba5588kkz\nbtw4M3DgQFdN3759zaRJk8wjjzxiYmNjTdOmTY0xxrz22mumR48epmHDhqZx48ZmyJAhZu/eva7t\nzvb2+uuvm0GDBpmwsDBz3XXXmfXr15tDhw6ZwYMHm/DwcJOQkGDWr19f7fvn7Vj69u3r9h7Ex8dX\nO97q1auNzWYz33zzTbW1f/vb30xoaKg5deqUa92KFSuMzWYzBw8eNMYY8+GHHxqbzeb2HmRlZRmb\nzWbWrFnj9p588skn1e7zrLPbnM+XX35pbDabSU9Pd607ceKECQ4ONmlpacYYYxwOh4mPjzcPPfSQ\nz/s1xpiTJ0+a8ePHm9jYWBMcHGyuvfZa85vf/Mb1uudcObs8f/5807x5c1O/fn3zy1/+0lRUVJgF\nCxaYFi1amEaNGpnJkyebsrIy13Z9+/Y199xzj9u+n376adOqVSvX8uOPP27atm3rWj5w4ID56U9/\napo1a2bCwsJM586dzSuvvOLWi+fPxtq1a40xlX/ejh49akaNGmUiIyNNaGio6devn9m6davr9bNz\nJT093fTu3duEhYWZhIQE8/7771f7HqalpZnrrrvOBAUFmbi4OPPII4+YioqKKnu02+3nHWfevHmm\nS5cupn79+iY2NtaMHj3a5OTkeN2353vmydfjXrFihenVq5cJCQkxf/3rX01cXJz5+9//7qq7++67\njc1mM9nZ2a51zZs3N//7v/9rjDFm27Zt5pZbbjExMTGmfv36pnv37uaDDz5w66Vly5bmkUceMVOm\nTDHR0dGmZ8+exhhjPv74Y9O5c2cTEhJifvSjH5mPP/640vfP08svv2wCAwPd1n3yySeu7+HZn6nX\nXnvN3HrrrSY8PNzMmjWryp/PY8eOmfHjx5smTZqYkJAQ06FDB7NkyRLX6/v27TN33HGHiYyMNI0a\nNTKDBg0yn3/++Xl7AwB/4tJ3APCR8fGM4Lln33Nzc/Xzn/9cY8eO1e7du/XZZ5/pgQceUGBgoFq0\naKH//Oc/kqQtW7YoNzdX//73vyVJL774oh577DE99NBD2r17tx588EHNmjVLS5Ysce1nwoQJ+vzz\nz7VixQqtWrVKBw8e1H/+859KZ/6XLVum/Px8rV69Wunp6ZK+P+P12GOPKSMjQytXrlRAQICGDBmi\n8vJyt20fffRRTZs2TZmZmerYsaNGjx6tcePGacqUKcrIyFBCQoLGjBmjioqK874fCxcu9Hos77zz\njmbOnKlWrVopNzdXW7Zs8el99tUnn3yiH//4x25n4lJSUmS3211nyz/55BO1bt3a7T7qhIQExcXF\nacOGDW7j+ToPfO2tXr16GjBggGtdZGSkevTo4drv9u3bdfDgQTVv3lx9+vRRkyZNdNNNN2n58uVe\nx37kkUeUkZGhd999V9nZ2XrzzTeVkJDger2qq0Q2b96s7du3a9WqVXrjjTe0dOlSDRkyRFu3btVH\nH32kV199Va+88ooWL17sdZzqnD59WgMHDtQHH3ygXbt2afLkyZowYYLWrFkjSZo/f7569+6tUaNG\nKTc397y3QxhjNHz4cO3du1crVqzQ5s2b1aRJE6WkpCg/P9+t9re//a0eeeQR7dy5UzfeeKNGjRql\nwsLC8/a4YsUKTZo0SePGjVNWVpbmzp2rhQsX6sknn3T1OG/ePAUEBCg3N1c5OTnnHctms2nu3Lna\ntWuX3nnnHR0+fFijR4++oPfsYo975syZmj17tr744gvdfvvtSk5O1scff+x6ffXq1YqJiXGt+/LL\nL3X06FElJydL+v7WjDvvvFNr1qxRRkaGBg8erJ/85Cfat2+f237mz5+v2NhYffbZZ3r55Zd19OhR\nDR06VN27d1dGRobmzp2rGTNmXNTxnr2ypKyszLXu97//vX7xi18oKytLv/rVryptc+bMGfXt21ef\nf/65Xn/9dX3xxRd68cUXFR4eLkk6duyYbr75ZsXGxmrDhg3atGmTOnTooH79+un48eMX1ScA1KjL\n+3cCALgyjBs3zgQGBpr69etX+ud5hqhVq1YmNTXVGGPM9u3b3c7celq/fr2x2Wzm0KFDbuvj4uLM\n73//e7d1DzzwgGndurUxxpi9e/cam81mPv74Y9fr5eXl5tprr3U789+3b1/ToUOHao8vPz/f2Gw2\ns3HjRmPMD2eC//znP7tqtmzZYmw2m3n++edd6zIyMozNZjNZWVnnHbu6YzGm+rOHni7kjHpKSooZ\nO3ZspfXXXHONee6554wxxtx7772mV69elWq6d+9u7rvvPmPMD+9JWFhYpTnw2WefVbnv6s6op6am\nmmbNmlVaP2LECDN06FBjjDH//Oc/jc1mM40aNTJ///vfzY4dO8wf/vAHY7fb3c7Eexo2bJgZP378\neV+v6ox6kyZNTHl5uWvdkCFDzDXXXON2Bn3YsGFmxIgRruV+/fqZe++9123s6s6on6/fc8cZOHCg\nmTBhQqW6c3/eVq5caWw2m9mzZ4/r9dLSUtO0aVPz1FNPGWN+mCvvvPOOq+bYsWPGZrOZjz766Lz9\n3HzzzWbUqFFu6/785z+b0NBQ13tU1dlfX5z9vXD06NHz1jz++OPGbre7zbOIiAjjdDov6LhfffVV\nt3HT0tJMkyZNjDHf/x4JCwszf/jDH8zo0aONMca8+OKLpmXLll77v/76612/44z5/oz6uXPJGGMe\nfvhh06pVK+NwOFzr3nvvvQs+o56Xl2eGDh1qGjZsaL799lvXz9S5Vy0ZU/mKl0WLFpmQkJDz/o54\n/PHHXWf+z3I6nVVevQQAlwP3qAOAj3r27KmlS5e6rTPGeH2a9fXXX6/BgwerU6dOSklJUb9+/XTH\nHXcoLi7uvNucOnVK33zzjfr06eO2vk+fPvrzn/+skpIS7d6929XTWYGBgUpKStJ3333ntt0NN9xQ\naR+ZmZl68skntWPHDh0/ftx1lvjQoUNuZy6vv/5619dNmjSRJP3oRz+qtC4vL8/tbO2FHMvZs2W1\nxdezvcbHM+VpaWmV3tPmzZu7vq5fv75rn2fHPPdsvs1m06lTp3zalyTX/e333nuv7rnnHknffw8+\n++wzLViwQAMHDqxyu6lTp+pnP/uZtm7dqgEDBuiWW27R4MGDvb4f1113nds9zE2aNFGHDh1Ur149\nt3VffPGFz/1Xpbi4WE899ZTee+895eTkqKysTKWlpa6zuL7KyspSdHS0Onbs6FoXFBSkG2+8UVlZ\nWW61Xbp0cX0dExOjgIAAHTt27Lxj7969W3feeafbuj59+qikpET79+9Xhw4dfO5zzZo1+uMf/6g9\ne/aosLDQ9T09dOiQmjZtet7trr32Wrez39L38+dCjtvzeQr9+vVTXl6esrKytGHDBvXu3VuDBw/W\nggULJEkff/yx+vfv76r/9ttv9fjjj2v16tXKzc1VRUWFSkpK3B6AabPZKu1n9+7d6tGjh+z2Hy7e\n7NWrl9f36SyHw+H6mTl9+rSuu+46vf3222rcuLHr91t1D43btm2bEhMT1axZsypf37Jli7Zt21bp\nnveSkhJlZ2f71CcA1CaCOgD4KCQkRK1bt76gbex2u95//31t2bJFK1eu1Ntvv61Zs2bpX//6l4YM\nGXLJPXmGLs+wabPZXJd6nlVcXKxBgwapT58+SktLU5MmTWSMUWJiotulpZLcAtrZfVW1ztvD0i63\npk2b6uuvv3ZbV15eroKCAldIatq0qVatWlVp29zc3EpBqnnz5l7nwc6dO11fHzlyRP369dOOHTvO\n29vZP5Sc+708duyYK4Sd3X9iYqLbtgkJCfrwww/P28egQYN0+PBhffjhh1qzZo3uuusude7cWatW\nrXILT+fyfNCYzWarct2532+73V5p3nneQuHpwQcf1LvvvqsXXnhBHTp0UFhYmGbOnKmTJ0963c5X\nnu+npCofbuaPeXv48GHddtttGjdunJ544gk1btxYX3/9tQYOHFjp581TvXr1Luh3TlXH7fnz37Jl\nS7Vu3VqrVq3Sxo0blZycrG7duqm0tFS7du3S2rVrNXfuXFf9+PHjdeTIEc2ZM0fx8fEKCQnR6NGj\nK/XuuR+bzXbRt4kEBARox44dstlsiomJqTR2Vfurirf9G2M0cOBA/eUvf6n0WsOGDS+sYQCoBdyj\nDgA+utD7cM/VvXt3zZ49W2vXrlXfvn318ssvS/ohPDgcDldtgwYNFBcXp7Vr17qNsXbtWrVu3Voh\nISGus9fnfmZwRUWFtm3bVm0ve/bs0fHjx5Wamqo+ffqoQ4cOKigoqJUnz/tyLLWtV69e+vTTT1VU\nVORal56eLqfT6TrDd/PNN+urr75yO5O2e/duHTlyRDfffPMF7a9169aufy1atKi07tzg1atXL5WX\nl7v9kaCwsFCbN2927TcpKUmhoaHas2eP236+/PJLxcfHe+2lUaNGGj16tF566SWtWLFCa9eurTTO\nuS5mjsfExLie0H3W9u3bvY61fv163XXXXRoxYoQ6d+6s+Ph4ffnll27bBAUFeX32gfT9Hy/y8/Pd\njqm0tFSbNm1Sp06dLvhYPMeuat6GhYWpTZs2Po+zZcsWlZSUaN68ebrpppvUrl075ebmXnJvl3Lc\n/fv318qVK7VmzRoNGDBAdrtdffv21QsvvKDjx4+7Xdmwfv16TZ06VUOHDlViYqJiY2O1f//+aveR\nkJCgzZs3u/0xxPMTFLxp3bq14uPjfQrkVUlKStLu3bsrzc1zX9+1a5frD2/n/ruYT+4AgJpGUAcA\nH/kaZM+t27hxo55++mlt3rxZhw8f1qpVq7Rz507X2dGWLVvKbrdrxYoVysvLc51RnD17thYsWKBF\nixZp3759+tvf/qaXXnpJDz30kCSpXbt2uv322zVt2jStW7dOu3fv1i9/+UudOnXKLewYYyr13bJl\nSwUHB2v+/Pnav3+/Vq1apRkzZlzSHyK8qe5YLsThw4eVmZnpCtRZWVnKzMzUiRMnXDV33323xo0b\n51oeM2aMGjdurDFjxmjnzp1avXq1pk2bptGjR6tly5aSvv/M7m7duumuu+7Sli1btGnTJt199926\n6aabKl22n5+f73rA2dl/nrcb+Kp9+/YaNmyYpkyZonXr1ikzM1NjxoxRXFycRo0aJen7S+mnT5+u\nhQsX6o033lB2drbmzZun9957z+vDuR5++GG98847+vLLL7Vv3z69+uqrioiIcP3xoCoX88eagQMH\nauXKlXrrrbeUnZ2tP/3pT9qwYYPXsTp06KDly5dry5Yt2r17tyZPnqycnBy3beLj47Vt2zYdOHBA\nx48frzK0DxgwQD169NCYMWO0ceNG7dq1S3fffbfKyso0ZcqUCz6Wc82ePVtvv/22nn32We3du1fL\nli3Tk08+qZkzZ17QR5y1b99eNptNzz33nL766istX75cTz/99CX1dqnHnZycrPfff19lZWXq1q2b\na90//vEPtWvXzu1Wjg4dOujVV1/Vrl27lJmZqTvvvFNOp9Pte1XV93rKlCn69ttvNXnyZO3Zs0er\nVq3Sww8/fEnHfSHuvPNOtWzZUj/5yU+0atUqffXVV1q1apWWLVsmSbrvvvvkcDg0bNgwbdiwQQcP\nHtSGDRv08MMPV/pIPAC4HAjqAOCDC3my9bl1kZGR+uyzzzRs2DC1b99ekyZN0l133aVHH31U0vf3\n+/7xj3/Un/70JzVr1kw//elPJX3/P7lPPfWUnnnmGSUmJmrOnDl69tlnNWHCBNfYL7/8sjp16qRb\nb71VycnJiouL06BBgxQcHOy178aNG+vVV19Venq6OnXqpN/97neaO3dupcuhqzpeX9edy5dj8fX9\nfeyxx9StWzdNnjxZNptNgwcP1g033KD//ve/rpqvv/7a7VL38PBwrVy5UmVlZbrpppv085//XLfc\nckulJ5e/9957atGihQYMGKBBgwapXbt2rqfyn2vYsGFq1qyZ2z9vf3So7rheeeUV9e/fXz/96U9d\nZ/g/+ugjt+/jM888owceeEC/+93v1KVLF73++utavny5+vXrd95xQ0ND9dhjjykpKUndu3fXrl27\n9P7777vuyfV8z6v6Hviybty4cZo2bZqmTZum7t2765tvvtH999/vdewXXnhBLVu2VP/+/TVw4EBd\ne+21GjFihFvNzJkz1bhxY11//fWKiYlxu3rkXMuXL1fHjh01ZMgQ9ejRQ3l5eUpPT1dUVJTb/i/U\nrbfeqiVLlmjp0qXq3LmzfvOb32jatGl6/PHHK70f3nTu3FkLFizQ3/72NyUmJur555/XvHnzqt2u\nup+JSznu/v37y+FwqG/fvq6a5ORkORyOSs8JePnll+V0OtWjRw/dcccduu2229S9e/dK319PzZo1\n03//+19t3rxZXbt21QMPPKAXXnjB6zFX13d1r5+7PjQ0VGvXrlWnTp00evRoJSQkaPr06SopKZH0\n/ZUgn376qRo3bqw77rhDHTt21F133aWvv/76vPe1A4A/2UxtXOsIAPA7h8Ohjh07avjw4ZozZ87l\nbgcAAAAX6ao9o+751FOAOYGqXMnzYv369Xrrrbe0f/9+ZWZmauLEiTp8+LDGjx9/uVu7ol3JcwK1\ngzkBT8wJeGJOwNOlzgmCOuoM5gSqciXPC4fDodTUVHXp0kXJyck6ePCgVq9eXenp4LgwV/KcQO1g\nTsATcwKemBPwdKlzgo9nA4ArVL9+/ZSRkXG52wAAAEANu2rPqAMAAAAAcCXiYXIAAAAAAFjIVX3p\n+9GjRy93C7CQiIgIFRUVXe42YDHMC3hiTsATcwKemBPwxJyAp0v9qEcufQcAAAAAwEII6gAAAAAA\nWAhBHQAAAAAACyGoAwAAAABgIQR1AAAAAAAshKAOAAAAAICFENQBAAAAALAQgjoAAAAAABZCUAcA\nAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1AAAAAAAshKAOAAAAAICFENQBAAAAALAQ\ngjoAAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1AAAAAAAshKAOAAAA\nAICFENQBAAAAALAQgjoAAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1\nAAAAAAAshKAOAAAAAICFENQBAAAAALAQgjoAAAAAABYS6K8dZWZmKi0tTU6nU8nJyRo+fHilmiVL\nligzM1PBwcGaOnWq4uPjJUmnT5/WSy+9pCNHjkiSpkyZovbt2/urdQAAAAAA/MYvQd3pdGrx4sV6\n9NFHFRUVpdmzZyspKUlxcXGumu3bt+vYsWOaP3++9u3bp0WLFik1NVWS9PLLL6tr166aOXOmHA6H\nSktL/dE2AAAAAAB+55dL37OzsxUbG6uYmBgFBgaqV69e2rp1q1vN1q1b1bdvX0lSu3btdPr0aRUW\nFqq4uFhffPGFkpOTJUkBAQEKCwvzR9sAAAAAAPidX86oFxQUKDo62rUcFRWl7OxsrzXR0dEqKCiQ\n3W5XgwYN9OKLL+rQoUOKj4/XhAkTFBwc7I/WAQAAAADwK7/do+4LY0yldQ6HQ1999ZUmTpyotm3b\nKi0tTcuXL9eoUaPc6rKyspSVleVaHjlypCIiImq9Z1w5goKCmBOohHkBT8wJeGJOwBNzAp6YE6jK\nsmXLXF8nJiYqMTHR5239EtSjoqKUn5/vWs7Pz1dUVJTPNVFRUWrbtq0kqWfPnlq+fHmlfVR14EVF\nRTV2DLjyRUREMCdQCfMCnpgT8MScgCfmBDwxJ+ApIiJCI0eOvOjt/XKPeps2bZSbm6u8vDxVVFRo\n48aNSkpKcqtJSkrSunXrJEl79+5VeHi4IiMjFRkZqcaNG+vo0aOSpJ07d7o9hA4AAAAAgKuJX86o\nBwQEaOLEiUpNTXV9PFtcXJzS09MlSSkpKerWrZsyMjI0ffp0hYSEaMqUKa7tJ0yYoAULFqiiokJN\nmjTR1KlT/dE2AAAAAAB+ZzNV3Rh+lTh7Fh6QuCQJVWNewBNzAp6YE/DEnIAn5gQ8NWvW7JK298ul\n7wAAAAAAwDcEdQAAAAAALISgDgAAAACAhRDUAQAAAACwEII6AAAAAAAWQlAHAAAAAMBCCOoAAAAA\nAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAAAACwEII6AAAAAAAWQlAH\nAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAAAACw\nEII6AAAAAAAWQlAHAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCGBl7sBAKgpeWfylFOc\n47WmaVhTxYTG+KkjAAAA4MIR1AFcNXKKc/Ta7te81oxNGEtQBwAAgKVx6TsAAAAAABZCUAcAAAAA\nwEII6gAAAAAAWAhBHQAAAAAAC+FhcgCuCL480b24vNhP3QAAAAC1h6AO4IrgyxPdU1ql+KkbAAAA\noPZw6TsAAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAAC+FhcgDqFIdxaEf+DtdycFGwSstK\n3WqahjVVTGiMv1sDAAAAJBHUAdQx+WfylX4w3bUcVC9IZeVlbjVjE8YS1AEAAHDZcOk7AAAAAAAW\nQlAHAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAA\nAACwEII6AAAAAAAWQlAHAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISg\nDgAAAACAhRDUAQAAAACwkMDL3QAAWI3DOLQjf4fXmqZhTRUTGuOnjgAAAFCXENQBwEP+mXylH0z3\nWjM2YSxBHQAAALWCS98BAAAAALAQzqgDuOzyzuQppzjHa01xebGfugEAAAAuL4I6gMsupzhHr+1+\nzWtNSqsUP3UDAAAAXF5c+g4AAAAAgIUQ1AEAAAAAsBCCOgAAAAAAFkJQBwAAAADAQgjqAAAAAABY\nCEEdAAAAAAALIagDAAAAAGAhBHUAAAAAACyEoA4AAAAAgIUE+mtHmZmZSktLk9PpVHJysoYPH16p\nZsmSJcrMzFRwcLCmTp2q+Ph4SdK0adMUGhoqu92ugIAA/fGPf/RX2wAAAAAA+JVfgrrT6dTixYv1\n6KOPKioqSrNnz1ZSUpLi4uJcNdu3b9exY8c0f/587du3T4sWLVJqaqrr9SeeeEL169f3R7sAAAAA\nAFw2frn0PTs7W7GxsYqJiVFgYKB69eqlrVu3utVs3bpVffv2lSS1a9dOp0+fVmFhoet1Y4w/WgUA\nAAAA4LLyyxn1goICRUdHu5ajoqKUnZ3ttSY6OloFBQWKjIyUzWbT008/LbvdroEDB2rgwIH+aBsA\nAAAAAL/z2z3qvjjfWfOnnnpKUVFROnXqlJ5++mk1b95c1113nVtNVlaWsrKyXMsjR45URERErfaL\nK0tQUBBzwqKCi4IVVC/Ia01AQECt1NgD7ApSkNeaKnsOCmY+XaX4XQFPzAl4Yk7AE3MCVVm2bJnr\n68TERCUmJvq8rV+CelRUlPLz813L+fn5ioqK8rnm7H8bNGigHj16KDs7u1JQr+rAi4qKavQ4cGWL\niIhgTlgr+uoOAAAgAElEQVRUaVmpysrLvNY4HI5aqQlSUKVtfBmntKyU+XSV4ncFPDEn4Ik5AU/M\nCXiKiIjQyJEjL3p7v9yj3qZNG+Xm5iovL08VFRXauHGjkpKS3GqSkpK0bt06SdLevXsVHh6uyMhI\nlZaW6syZM5KkkpIS7dy5Uy1atPBH2wAAAAAA+J1fzqgHBARo4sSJSk1NdX08W1xcnNLT0yVJKSkp\n6tatmzIyMjR9+nSFhIRoypQpkqTCwkI999xzkr5/evzNN9+s66+/3h9tAwAAAADgd367R71r167q\n2rWr27qUlBS35UmTJlXarkmTJpozZ06t9gYAAAAAgFX45dJ3AAAAAADgG4I6AAAAAAAWQlAHAAAA\nAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAAAACwEII6\nAAAAAAAWQlAHAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCGBl7sBALgSOYxDO/J3eK1p\nGtZUMaExfuoIAAAAVwuCOgBchPwz+Uo/mO61ZmzCWII6AAAALhiXvgMAAAAAYCEEdQAAAAAALISg\nDgAAAACAhRDUAQAAAACwEB4mB6DW5J3JU05xTrV1xeXFfugGAAAAuDIQ1AHUmpziHL22+7Vq61Ja\npfihGwAAAODKwKXvAAAAAABYCEEdAAAAAAALIagDAAAAAGAhBHUAAAAAACyEoA4AAAAAgIUQ1AEA\nAAAAsBCCOgAAAAAAFkJQBwAAAADAQgjqAAAAAABYCEEdAAAAAAALIagDAAAAAGAhBHUAAAAAACyE\noA4AAAAAgIUQ1AEAAAAAsBCCOgAAAAAAFkJQBwAAAADAQgjqAAAAAABYCEEdAAAAAAALIagDAAAA\nAGAhBHUAAAAAACyEoA4AAAAAgIUQ1AEAAAAAsBCCOgAAAAAAFkJQBwAAAADAQgjqAAAAAABYCEEd\nAAAAAAALIagDAAAAAGAhBHUAAAAAACyEoA4AAAAAgIUQ1AEAAAAAsJDAy90AAFytHMahHfk7vNY0\nDWuqmNAYP3UEAACAKwFBHQBqSf6ZfKUfTPdaMzZhLEEdAAAAbrj0HQAAAAAACyGoAwAAAABgIQR1\nAAAAAAAshKAOAAAAAICFENQBAAAAALAQgjoAAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAA\nCyGoAwAAAABgIQR1AAAAAAAshKAOAAAAAICFENQBAAAAALAQgjoAAAAAABYS6K8dZWZmKi0tTU6n\nU8nJyRo+fHilmiVLligzM1PBwcGaOnWq4uPjXa85nU7NmjVLUVFRmjVrlr/aBgAAAADAr/xyRt3p\ndGrx4sV66KGH9Pzzz+uTTz7RkSNH3Gq2b9+uY8eOaf78+Zo8ebIWLVrk9vr//d//KS4uTjabzR8t\nAwAAAABwWfglqGdnZys2NlYxMTEKDAxUr169tHXrVrearVu3qm/fvpKkdu3a6fTp0yosLJQk5efn\nKyMjQ8nJyTLG+KNlAAAAAAAuC79c+l5QUKDo6GjXclRUlLKzs73WREdHq6CgQJGRkVq6dKnuuusu\nnTlzxh/tAvBB3pk85RTneK0pLi/2UzcAAADA1cNv96j7oqqz5du2bVODBg0UHx+vrKys826blZXl\n9vrIkSMVERFRK33iyhQUFMScqEFfFH2hf+37l9ea5BbJCqoXVO1YAQEB1dbVVo09wK4gBXmtqc1+\ngoOCmZcWw+8KeGJOwBNzAp6YE6jKsmXLXF8nJiYqMTHR5239EtSjoqKUn5/vWs7Pz1dUVJRPNZ99\n9pm2bdumjIwMlZeX68yZM/rLX/6i++67z237qg68qKioFo4GV6qIiAjmRA0qLStVWXmZ1xqHw1Ft\nja91tVUTpKBK2/izn9KyUualxfC7Ap6YE/DEnIAn5gQ8RUREaOTIkRe9vV+Ceps2bZSbm6u8vDxF\nRUVp48aNmjFjhltNUlKSPvzwQ/Xq1Ut79+5VeHi4IiMjNWbMGI0ZM0aStHv3br377ruVQjoAAAAA\nAFcLvwT1gIAATZw4Uampqa6PZ4uLi1N6erokKSUlRd26dVNGRoamT5+ukJAQTZkypcqxeOo7AAAA\nAOBq5rd71Lt27aquXbu6rUtJSXFbnjRpktcxEhISlJCQUOO9AQAAAABgFX75eDYAAAAAAOAbgjoA\nAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1AAAAAAAshKAOAAAAAICF\nENQBAAAAALAQgjoAAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1AAAA\nAAAshKAOAAAAAICFENQBAAAAALAQgjoAAAAAABZCUAcAAAAAwEICL3cDAFCXOYxDO/J3eK1pGtZU\nMaExfuoIAAAAlxtBHQAuo/wz+Uo/mO61ZmzCWII6AABAHcKl7wAAAAAAWAhBHQAAAAAACyGoAwAA\nAABgIQR1AAAAAAAshKAOAAAAAICF+BTUt2zZIofDUdu9AAAAAABQ5/kU1N98801NnjxZixcv1r59\n+2q7JwAAAAAA6iyfPkf9ueee08GDB7Vu3TrNnTtXwcHB6tOnj3r37q2YGD7bFwAAAACAmuJTUJek\nVq1aqVWrVvrFL36hzz//XK+88oqWLVumjh07asCAAbr55ptlt3PLOwAAAAAAl8LnoC5Jubm5Wrdu\nnTZs2CCbzaaRI0fqmmuu0YcffqhNmzbpwQcfrK0+AQAAAACoE3wK6h988IHWr1+vo0eP6qabbtJ9\n992n9u3bu16/8cYbdc8999RakwAAAAAA1BU+BfWMjAwNHTpUN9xwg4KCgiq9HhwcrJkzZ9Z4cwAA\nAAAA1DU+BfWZM2fKbrcrMPCH8oqKCjmdTldw79KlS+10CAAAAABAHeLT099SU1N14MABt3UHDhzQ\nM888UytNAQAAAABQV/kU1A8dOqS2bdu6rWvbtq0OHjxYGz0BAAAAAFBn+RTUw8PDdfLkSbd1J0+e\nVEhISK00BQAAAABAXeVTUL/xxhs1f/58HT58WKWlpTp06JD+8pe/qGfPnrXdHwAAAAAAdYpPD5Mb\nPXq0/vGPf+ihhx5SeXm56tWrp/79+2vMmDG13R8AAAAAAHWKT0E9KChI99xzjyZOnKiioiJFRETI\nbvfpZDwAAAAAALgAPgV1SSouLtbRo0dVUlLitr5Tp0413hQAAAAAAHWVT0F9zZo1Wrx4sUJCQlyf\nm37WwoULa6UxAAAAAADqIp+C+htvvKHf/OY36tq1a233AwAAAABAnebTjeZOp1PXX399bfcCAAAA\nAECd51NQHzZsmN566y05nc7a7gcAAAAAgDrNp0vf33vvPZ08eVLvvvuuIiIi3F7761//WiuNAQAA\nAABQF/kU1KdPn17bfQAAAAAAAPkY1BMTE2u7DwAAAAAAIB+DellZmd566y1t3LhRRUVFWrp0qXbs\n2KGcnBzdcssttd0jgDqiue0GlRRGVvnad1+3VfRJ78/JqKomJLJQ35htNdYjAAAAUNt8CupLly5V\nQUGB7r//fj3zzDOSpGuvvVZpaWkEdQA1pqQwUqs/aljla4ciw3WgsOrXvNX0HyTJ+2YAAACApfgU\n1Ddv3qwFCxYoJCRENptNkhQVFaWCgoJabQ4AAAAAgLrGp49nq1evnhwOh9u6U6dOqUGDBrXSFAAA\nAAAAdZVPZ9R79uyphQsXaty4cZKkEydOKC0tTT/+8Y9rtTkA1ubtnnLJt/vKz61raLtW0qka7BAA\nAAC48vgU1O+880699tpr+u1vf6uysjLdf//9GjBggEaMGFHb/QGwMG/3lEu+3Vd+bt3QPqEiqAMA\nAKCu8ymo16tXT+PHj9e4ceN06tQpRUREyG736ap5AAAAAABwAXwK6seOHXNbLikpcX3dpEmTmu0I\nAAAAAIA6zKegfv/995/3tTfffLPGmgEAAAAAoK7zKah7hvHCwkItW7ZM1113Xa00BQAAAABAXeVT\nUPcUGRmp8ePH69e//rV69+5d0z0BuMzyzuQppzjHa01xebGfugEAAADqlosK6pJ09OhRlZaW1mQv\nACwipzhHr+1+zWtNSqsUP3UDAAAA1C0+BfXHHnvMbbm0tFRHjhzRz372s1ppCgAAAACAusqnoJ6c\nnOy2HBwcrJYtW6pZs2a10hQAAAAAAHWVT0G9X79+tdwGAOB8HMahHfk7vNY0DWuqmNAYP3UEAACA\n2uRTUP/nP/8pm81Wab0xxvW1zWbTqFGjaq4zAIAkKf9MvtIPpnutGZswlqAOAABwlfApqOfm5mrT\npk1q27atGjdurOPHjys7O1s33nijgoKCZIypMsgDAAAAAIAL4/NT32fMmKGePXu6ljdt2qRPP/1U\nU6dOrZXGAAAAAACoi3wK6hkZGbr//vvd1t1www168cUXa6UpAKgp1wRdK50c4Fr+7uu2ij7pdC0H\nBATI4XC4beNZc66QyEJ9Y7bVTrMAAACAfAzqsbGx+uCDD3Tbbbe51n300UeKjY2ttcYAoCYUnwrV\n6nUNXcuHIsN1oPCH5aqCumfNufoPklT1SwAAAECN8Cmo/+pXv9KcOXP0n//8R1FRUSooKFBAQIB+\n+9vf1nZ/AAAAAADUKT4F9fj4eM2fP1979+7ViRMn1KhRI7Vv316BgT7f4q7MzEylpaXJ6XQqOTlZ\nw4cPr1SzZMkSZWZmKjg4WFOnTlV8fLzKysr0xBNPqLy8XBUVFerevbvGjBnj+xECAAAAAHAF8T1p\n/382m00JCQkqKSlRRUWFQkJCqt3G6XRq8eLFevTRRxUVFaXZs2crKSlJcXFxrprt27fr2LFjmj9/\nvvbt26dFixYpNTVVQUFBevzxxxUcHCyHw6HHHntMX3zxhTp27HihrQMAAAAAYHk+BfXDhw/r2Wef\nVb169ZSfn68f//jH2r17t9auXasHHnig2u2zs7MVGxurmJjvP+O3V69e2rp1q1tQ37p1q/r27StJ\nateunU6fPq3CwkJFRkYqODhYklRRUSGn06n69etf8IECAAAAAHAlsPtS9Pe//10jR47UvHnzXJe7\nJyQk6IsvvvBpJwUFBYqOjnYtn73P3VtNdHS0q8bpdOrBBx/Uvffeq8TERLeADwAAAADA1cSnM+pH\njhxRnz593NYFBwerrKysRpsxxlS53m63a86cOSouLlZqaqqysrKUmJjoVpOVlaWsrCzX8siRIxUR\nEVGj/eHKFhQUVOfnRE6OUV5e9X+fO3GypWK+G+S15sw37dXIXqyAgKLz1tjtdgUEBFS7v7N1Npvt\nvPW+jFVVjeeYnjVV7dPbvgICAhRUL8j1X2/8WRMcFFzn53dN4XcFPDEn4Ik5AU/MCVRl2bJlrq8T\nExMrZVhvfArqjRs31v79+9W2bVvXuv379/v88WxRUVHKz893Lefn5ysqKuqCa8LCwtS1a1ft37+/\n0kFWdeBFRecPEKh7IiIi6vycOHw4SO+/H1xtXW5xPX3+rfdbTA5EhijhRw45HIXnrXE6nZU++sxb\nnTHmvPW+jFVVjeeYnjVVfTybt305HA6VlZe5/uuNP2tKy0rr/PyuKfyugCfmBDwxJ+CJOQFPERER\nGjly5EVv79Ol76NHj9azzz6rN998UxUVFfr3v/+tuXPnatSoUT7tpE2bNsrNzVVeXp4qKiq0ceNG\nJSUludUkJSVp3bp1kqS9e/cqPDxckZGROnXqlE6fPi1JKisr0+eff674+PgLOUYAAAAAAK4YPp1R\nv+GGG/TQQw9p5cqVSkhI0PHjx/Xggw+qdevWPu0kICBAEydOVGpqquvj2eLi4pSeni5JSklJUbdu\n3ZSRkaHp06crJCREU6ZMkSQVFhZq4cKFcjqdMsaoT58+6ty580UeLgAAAAAA1lZtUHc4HPr1r3+t\n559/Xvfee+9F76hr167q2rWr27qUlBS35UmTJlXarkWLFnr22Wcver8AAAAAAFxJqr30/ewDnmr6\nwXEAAAAAAKAyny59HzJkiObNm6fhw4crOjpaNpvN9VqTJk1qrTkAAAAAAOoar0G9sLBQkZGRWrJk\niSRp586dlWrefPPN2ukMAAAAAIA6yGtQnzFjhpYuXeoK43PmzNGDDz7ol8YAAAAAAKiLvN6jboxx\nW969e3etNgMAAAAAQF3n0+eoAwAAAAAA//B66bvT6dSuXbskfX923eFwuJbP6tSpU+11BwAAAABA\nHeM1qDds2FB//etfXcsRERFuy5K0cOHC2ukMAAAAAIA6yGtQJ4QDAAAAAOBf3KMOAAAAAICFeD2j\nDgBwd03QtdLJAfru67aKPun0WutLTUBx85psDwAAAFcBgjoAXIDiU6Fava6hDkWG60BhQ6+1vtS0\njwqtyfYAAABwFeDSdwAAAAAALISgDgAAAACAhRDUAQAAAACwEII6AAAAAAAWQlAHAAAAAMBCCOoA\nAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAAAACwkMDL3QAA4NI5\njEM78nd4rWka1lQxoTF+6ggAAAAXi6AOAFeB/DP5Sj+Y7rVmbMJYgjoAAMAVgEvfAQAAAACwEII6\nAAAAAAAWQlAHAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACA\nhRDUAQAAAACwEII6AAAAAAAWQlAHAAAAAMBCCOoAAAAAAFhI4OVuAMD5HT8eqLy8mvt72unTATU2\nFgAAAIDaQVAHLCwvz6733w+usfG6dXPU2FgAAAAAageXvgMAAAAAYCGcUQfqmNPlp1VUXuS1ptxZ\n7qduAAAAAHgiqAN1TFF5kT7/9nOvNa0jW/upG4SooaJPDvBa893XbRV90ulzTUhkob4x22qsRwAA\nAPgXQR0ALqNThUFa/VFDrzWHIsN1oND3mv6DJHkvBwAAgIVxjzoAAAAAABbCGXUA6tHEqFHZSddy\nTNgRtbMVet0mJuyImpeXqrzRSa811Y1zbl2ncnPe8TzHOhHUUJuP2aodGwAAALjSENQBqFHZSdXf\nttm1HB52jeoXf+t1m/CwaxR6rFT1957yWlPdOOfWhRY1OO94lca6oYekyGrHBgAAAK40XPoOAAAA\nAICFENQBAAAAALAQgjoAAAAAABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIXyO\nOlALAk6ckD0//5LHafxtoK4Pqrkf02u/M7ox9JSaNSp0W9+mgXSsxvYCAAAA4FIQ1IFaYM/PV73P\nPrvkcYJzAhW6q+Z+TIPznArLLlL9om/c1of16VRj+wAAAABwabj0HQAAAAAACyGoAwAAAABgIQR1\nAAAAAAAshKAOAAAAAICFENQBAAAAALAQgjoAAAAAABbCx7MBQB3hMA7tyN9RbV3TsKaKCY3xQ0cA\nAACoCkEdAOqI/DP5Sj+YXm3d2ISxBHUAAIDLiEvfAQAAAACwEII6AAAAAAAWQlAHAAAAAMBCuEcd\nwBWpdXSABqvQbV1M2BG1s7mv61RuVN7o5Hlr7Da7nMZZ7Tie43mrqWqcE0ENtfmYrfoDAwAAQJ1H\nUAdwRQovKVL9bbvc14Vdo/rF37qtCy1qoPp7T523xm63y+l0D+pVjeM5nreaKse5oYekyOoOCwAA\nAPBvUM/MzFRaWpqcTqeSk5M1fPjwSjVLlixRZmamgoODNXXqVMXHx+v48eNauHChTp48KZvNpgED\nBui2227zZ+sAAAAAAPiF34K60+nU4sWL9eijjyoqKkqzZ89WUlKS4uLiXDXbt2/XsWPHNH/+fO3b\nt0+LFi1SamqqAgMDNW7cOLVq1UolJSX6/e9/rx/96Edu2wIAAAAAcDXw28PksrOzFRsbq5iYGAUG\nBqpXr17aunWrW83WrVvVt29fSVK7du10+vRpFRYWKjIyUq1atZIkhYSEqHnz5jpx4oS/WgcAAAAA\nwG/8FtQLCgoUHR3tWo6KilJBQYHXmujo6Eo1eXl5OnjwoNq1a1e7DQMAAAAAcBlY7uPZjDHnfa2k\npETPP/+8xo8fr5CQED92BQAAAACAf/jtHvWoqCjl5+e7lvPz8xUVFeVzTUVFhebOnavevXurR48e\nlcbPyspSVlaWa3nkyJGKiIio6cPAFSwoKMh/cyIkRPVq4I9JgfWkgMCAGmjoewEBUoA9QHa7+9/o\n7LK7rbPZbJVqPNlsNtnkvc6Xcc6t8zae51iePZ9vf55jetZUuY2Xvs+O5+t7dLbGbrMrIKDy9/J8\n691q7BdW0ySkhQK+G1Sp5sw37RXzXfUfEXficEsdONnAtRwT41TTpnXno+X8+rsCVwTmBDwxJ+CJ\nOYGqLFu2zPV1YmKiEhMTfd7Wb0G9TZs2ys3NVV5enqKiorRx40bNmDHDrSYpKUkffvihevXqpb17\n9yo8PFyRkZEyxuill15S8+bNNWTIkCrHr+rAi4qKau14cOWJiIjw25yoV1IiR0nJJY9TUR4oR0UN\nNPT/ORxOOZyOSh9H5pTTbZ0xplKNJ2OMjLzX+TLOuXXexvMcy7Pn8+3Pc0zPmqo+ns1b32fH8/U9\nOlvjNE45HI5KNedb71bjvLCa7wpDtHJd/Uo1ByJDdKCw8npPx66pp9iwH65uuvXWMtWvX1btdlcL\nf/6uwJWBOQFPzAl4Yk7AU0REhEaOHHnR2/stqAcEBGjixIlKTU11fTxbXFyc0tPTJUkpKSnq1q2b\nMjIyNH36dIWEhGjKlCmSpC+//FLr169XixYt9Lvf/U6SNGbMGHXp0sVf7QMAAAAA4Bd+/Rz1rl27\nqmvXrm7rUlJS3JYnTZpUabuOHTvqzTffrNXeAAAAAACwAss9TA4AAAAAgLqMoA4AAAAAgIUQ1AEA\nAAAAsBCCOgAAAAAAFkJQBwAAAADAQgjqAAAAAABYCEEdAAAAAAAL8evnqAMX6tSpUzp58mSNjBUa\nGqozZ87UyFjVCSwoUKDDccnjFARKRQ1NDXT0veMhTp1pZGRCgtzWnwksk7nmh3WlwZIJD/Lc3E1A\neIik0hrrDQAAAMD3COqwtJMnT2rPnj01MlZYWJiKi4trZKzqBObkKODLLy95nBMnAnToUEANdPS9\ngpNOHco5o/wz+W7rTwaX6lj2Iddyw+CGOlnq/Q8kDRMbSt6zPAAAAICLwKXvAAAAAABYCEEdAAAA\nAAALIagDAAAAAGAhBHUAAPD/2rv/4Kiq+//jr9272WwSsgm7IfIjRYOgBiwURKGDlYKD/VCmVavG\njh9nCur0W1Fa/UNFZhw702FaiqBioTgFpdqfjAUsQ2ccW38h2KkIjJYfYlAEhJCQhE1INvvr7vcP\nPizJZrPZJJvdm+zzMdMZ995zz74Dh2Zfe889BwAAWAhBHQAAAAAACyGoAwAAAABgIWzPBqRJc7Mh\nn+/Cd19un6G8pv5vqxYK2frdBwAAAIDBhaAOpInPZ9dHH+VJkiodhpxp2P98xAiz330AAAAAGFwI\n6gCATqKKqratNvb6VJtfoYavOrUZVThK5QXlmS4NAAAgJxDUAQCdtIXa9Pm5z2Ovy2p9avD/q1Ob\n/534vwR1AACAAcJicgAAAAAAWAhBHQAAAAAAC2HqOzCEhMyQQmaw2/PBiCEzygJ12TDOa+g7Otfl\n+Hj/MX1neNfjHZUXntQEW+ptrg1FFRru61M/Hds1OUv0nzPsPAAAAJBpBHVgCAmZQTX4G7o9Xxwq\nlBnNy2BFuKiovUXDPvpvl+MFga807L9fJbiiw7WFIzSsrT7lNgUtbg070tynfjq1u+4GSaU9tgcA\nAEB6MfUdAAAAAAAL4Y46ACCpEc6vSb6bOx2r/2KMDp5x9qm/8nJTZWXhdJQGAAAwJBHUAQBJtTUX\n6O33SjodOzuiQCML8/vU3/z5AZWVpaMyAACAoYmp7wAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1\nAAAAAAAshMXkAPSJw+NSixFSdET3K38H8qVoUc8rg19s11Jk77a/+L78jmCXtoneL77P+DZRu11R\n00y57ov9pfKzdWyTqF5JaiqIJv0z7KmeRG3M0kJJXfdRBwAAwOBAUIflGa2tsjWnIXQ4nXIEg/3v\npxtun6FKhyFJKsn3yz9g72QN/lC7Th89omOn2rptU5JfIl/A12NfF9sF2wq77S++L19+QGdqvuzx\n/eL7jG9js9kVjXYO6snqvthfKj9bxzaJ6pUkI9igI8cbUu4nlTbjKr6RtC0AAACsjaAOy7M1N8v4\n9NN+92M4nTIGMKjnNRlyfnkhqOddPWLIB3UAAAAAA4Nn1AEAAAAAsBCCOgAAAAAAFkJQBwAAAADA\nQgjqAAAAAABYCEEdAAAAAAALIagDAAAAAGAhBHUAAAAAACyEoA4AAAAAgIUQ1AEAAAAAsBCCOgAA\nAAAAFkJQBwAAAADAQgjqAAAAAABYiCPbBQAABp+ooqptq03apjivWEV5RRmqCAAAYOggqAMAeq0t\n1KbPz32etM3XR3ydoA4AANAHTH0HAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAvhGXUAQEYFg3YdPOhM\nW3/l5abKysJp6w8AACDbCOoAgIxqarJr7968tPU3f35AZWVp6w4AACDrmPoOAAAAAICFcEcdAJDQ\nOK+h7+icrg1FFRru63SuvPCkJtjOJb3+awXHVOps6nr8fFQRpy1tdZbVh5VXE5bp9SoyfHja+gUA\nAMgWgjoAIKGi9hYN++i/Kmhxa9iR5s7nCkdoWFt90usLi79SgbOky/H8OlMFx9M3oSu/Pay8UWGF\nZs4kqAMAgCGBqe8AAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAA\nAACwEII6AAAAAAAWQlAHAAAAAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALMSR\n7Q0AWtIAABywSURBVAIApCZkhhQyg0nbmFEzQ9UAAAAAGCgEdWCQCJlBNfgbkrYpyS/JUDVAanxB\nX5dj/rBDvmBYkuQyXMo38jNdFgAAgKVlNKjv379fmzZtkmmamjt3rm677bYubV566SXt379f+fn5\nWrx4sSorKyVJ69at0759++R2u7Vq1apMlg0A6INgJKj6tvoux70Bt75qaZYkjSkeQ1AHAACIk7Fn\n1E3T1MaNG7Vs2TKtXr1au3bt0smTJzu12bt3r86cOaM1a9boxz/+sTZs2BA7N2fOHC1btixT5QIA\nAAAAkBUZC+o1NTUaOXKkysvL5XA4NGvWLO3Zs6dTmz179mj27NmSpAkTJqi1tVXnzp2TJFVVVamo\nqChT5QIAAAAAkBUZm/re2Ngor9cbe+3xeFRTU5O0jdfrVWNjo0pLSzNVJnJIe7tdfr8tbf2FQunr\nC0DqwmGbTp92KPCVQ2eDzn7353LZ5HY7VFYWTkN1AAAAvWe5xeSi0Wifrjtw4IAOHDgQe11dXa3i\n4uJ0lYUsKSgokNPplOHs/4dvu90uZ4d+mpulEyfTN6lkRJkpu3GhP7vNHvvv/rDbdanPiF02W099\n2pK2scmWsI1N8X0n70eSbDabbD2267mfju2S99f5XNeaE79f1z7j+rHZ1HVyUfd1XOovlZ/tUpvE\n9SaqL3k/qbWxyW5P8F62xMe7a2eXXXb7hb+X+OtS6au7Nh37M+yGDIfRY03JBALSZzWG/MNc+m+4\n//+/bxiG/ud/pMrKvv0+wtDjdDr5TIFOGBOIx5hAIps3b47996RJkzRp0qSUr81YUPd4PGpouLRi\ndUNDgzweT6/bdCfRD97S0tKPimEFfr9fwWBQRjD5tmSpcDqdCnboJxIxZEb63W2MaZqx/syoKTPS\n/63S4vuM9rj9WjRpm6iiCdtEFd938n6kC1+qRXts13M/Hdsl76/zua41J36/rn3Gv7an1E/X/lL5\n2S61SVxvovqS95Nam6hMM8F7RRMf766dKVOmeeHvJf66VPrqrk3H/iJmRJFw//4hRiKmImEpEg53\n+jfeV06nU+3t7Wpp6X9fGBqKi4v5TIFOGBOIx5hAvOLiYlVXV/f5+owF9SuvvFK1tbWqq6uTx+PR\n7t279bOf/axTm+nTp+uNN97QrFmzdOTIERUVFTHtHQB6yeYuVHREW5fjgXwpWtTz7JSL7fyOoKIj\nnGopsis6wpmwTSr9xOvYn78gKrsj1GNNyZx1mWopsStkr5dhuPrVlyTZ7U41NbXpxIn+1ZVISUmJ\n3G532vsFAABDS8aCumEYuu+++7R8+fLY9mwVFRV68803JUnz5s3TtGnTtG/fPi1ZskQul0sPPvhg\n7PrnnntOhw4dUktLix588EFVV1drzpw5mSofAAaNtmC7jtQc6XK8JL9EvkDXfc27a+fLD+hMzZcK\nthXq2Km2hG1S6Sdex/68BQ0qdPRvodBGn6n6eruCkr4It/arL0lyOBxyufxqbU1/UK+qqiKoAwCA\nHmX0GfWpU6dq6tSpnY7Nmzev0+v7778/4bWPPPLIgNUFAAAAAIBVZGx7NgAAAAAA0DPLrfoOAEBf\nlBQEVek/0e9+7IYhty8ohyONq03+H0dpqfICgdhr0+tVZPjwtL8PAAAY3AjqAIAhIS94Xs6a+n73\nYzfsymsPyTg7AEFdUp5xaTu60MyZBHUAANAFU98BAAAAALAQgjoAAAAAABZCUAcAAAAAwEII6gAA\nAAAAWAhBHQAAAAAAC2HVdwBAVrWFW5Oez7M7lWfPy1A1AAAA2UdQBwBkTdgMyxfwJW3jLfAS1AEA\nQE5h6jsAAAAAABZCUAcAAAAAwEKY+g4AQJxIxKamJiNt/RUUROVymWnrDwAADG0EdQAA4gQCUn19\n+oL65ZdH5HKlrTsAADDEMfUdAAAAAAALIagDAAAAAGAhBHUAAAAAACyEoA4AAAAAgIUQ1AEAAAAA\nsBCCOgAAAAAAFkJQBwAAAADAQthHHQCALLEHAsqrqcl2GSkzvV5Fhg/PdhkAAAx5BHUAALLEfu6c\njC++yHYZKQvNnElQBwAgA5j6DgAAAACAhXBHHbCAkBlSyAx2OR6MGGoLRyRJZtTMdFkAAAAAsoCg\nLqmtrU2RSCTbZSCHhcygGvwNXY4XhwrV4G+TJJXkl2S6LMAy2sKt3Z4LRgyFTCOD1QAAAAwsgrqk\no0eP6vTp09kuA3EcDoeuuOKKbJcBIMvCZli+gK/b88WhQoXMggxW1HuRiE1NTYYaGw05w5d+9TpL\nDQVP9/5XcXGxqWHDmGUDAMBQRVCXFI1Gs10CEuDvBcBQEQhI9fWGXHmGzvou/eotG2bX2f/2/lfx\ntdeGCeoAAAxhBHUAADLENjxPLQrFXtvVrJaSUJIrEqt1RBTMwiNb4cZGhU+cyPj7ZkJJSYncbne2\nywAAQBJBHQCAjGkL+XW85mjstdtoVHNNfa/7aQ9FNHx45oN6RFL43LmMv28mVFVVEdQBAJbB9mwA\nAAAAAFgId9RzkNHaKltzc7bL6JHhcMhRVCRbIJDtUgAAAAAgYwjqOcjW3Czj00+zXUaPDKdThiS7\n35/tUgAAAAAgY5j6DgAAAACAhRDUAQAAAACwEII6AAAAAAAWwjPqAIAhIRAJqC3c2u35PLtTefa8\nDFY09NhCITlOn852GSmLut2KFBVluwwAAHqNoA4AGPTCZlgtwRY1+Bu6beMt8A6ZoB6J2NTUZKSt\nv4KCqFwus8d29tZW6cyZtL3vQItcfbVEUAcADEIEdQwa7e12+f22Pl9vGFIkcumDbSjU974AIJsC\nAam+Pn1B/fLLI3K50tYdAADoJ4I6Bg2/36Yvv+z7B1O7YZcZufR6xIie7x4BAAAAQKaxmBwAAAAA\nABZCUAcAAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1AAAAAAAshO3ZAAA5oy3cmvR8\nnt2pfCM/Q9VYRyRiU1NTz9tfOorsCqfQTpIKCqJyuQbPNpjt7e06ceJEwnMFBQXy+/0ZrgiSVFJS\nIrfbne0yACDjCOrAAAuZIYXMYNI2ZnTwfJgFBquwGZYv4EvaxlvgVb5yL6gHAlJ9fc8B3O2yq/nL\n1IL65ZdH5HL1t7L+sYVCcpw+nVLbNrtdJ776KuE5p9OpYDD5/4+nS9TtVqSoKCPvNRhUVVUR1AHk\nJII6MMBCZlAN/oakbUrySzJUDQDkDntrq3TmTEptHaYp47PPEp4znE4ZGQrqkauvlgjqAJDzeEYd\nAAAAAAALIagDAAAAAGAhBHUAAAAAACyEZ9QBAOigLdiq4ohNbeFIwvN5dqfy7HkZrgoAAOQSgjoA\nAP8nbIbVHGxRUcilBn9bwjbeAi9BPQWpbvmWqsG23RsAAP1BUAcAAGmX6pZvqbLCdm+Z0Jst5ayA\n7eQAYGAQ1AEAACyiN1vKWQHbyQHAwGAxOQAAAAAALIQ76gAAALCk9vZ2nThxIttldFFQUCC/35/t\nMrKqpKREbrc722UAQxZBHQAAWF5fFqdzFNkV7uYaFqcbHM6fP68vvvgi22V0UVhYqLa2xAtO5oqq\nqiqCOjCACOoAAPRSW7g16Xm2cEu/vixO53bZ1fxl4mtyZXE6AMDgRFAHAKAXwmZYvoAvaRu2cEOu\nGOhV6u12e9r7Z6V6AIMBQR3oh5AZUsgMSpKCEUNt4YgkKT/ijt1xM6NMrQQADE0DvUq9wzRlfPZZ\nWvtkpXoAgwFBHQOmvd0uv9/W5+sNp6Hh5+xqOX/hGcNQqO99DZSQGVSDv0GSVBwqVIP/wvNqjuCw\n2PGS/JKs1QcgezpOj7/4RV57pD12PM/uzFZpUNdn3p3N9m6fgTcMKRLpedo9z70DANKFoI4B4/fb\n9GU3zwamIs9ll7PIUL3/wjOGI0bw4QfA4BA/Pf7iF3negC/2JZ63wJut8qCuz7xHh9t1vJvfWXbD\nLjPSc588945cYtUV+bPFSjsBsCL/0EBQBwAAQM5Iy3P1TqccwWB6CuqBVZ+pt+qK/NlipZ0AWJF/\naCCoAwCQJYFIIOkK8qweP7j0ZQu5eB23lGMq/cBIx3P1htMpI0NBnWfqgdxEUAcAIAvCZlgtwZbY\nVPhEWD1+cOnLFnLxOm4px1R6SAO/sn5fJVuR36qzAIDBhKAOAABgQem4Q98Rd+gHp4FeWb+vkq3I\nzyyA7GL9AGsYPXp0v67PWFDfv3+/Nm3aJNM0NXfuXN12221d2rz00kvav3+/8vPztXjxYlVWVqZ8\nLZBuHbde6w5brwEYaImmxnfcDlJiivxQlY479B1VVJj92o0lXkFBVNzwRyJWnQXQnaE2A4D1A6xh\nxowZ/bo+I0HdNE1t3LhRTz31lDwej5588klNnz5dFRUVsTZ79+7VmTNntGbNGn322WfasGGDli9f\nntK1WRWNyhaNZrsKDICOW691h63XAAyk+NXjL+q4HaR0YYp8T18sEuaR7uB/+eURgjoSsuosgO4w\nAwBWlJGgXlNTo5EjR6q8vFySNGvWLO3Zs6dT2N6zZ49mz54tSZowYYJaW1t17tw51dXV9XhtNjlO\nn5YxmL6xKi5WZMSIbFcBAEij7gJ9R8nC/MU79IR59EYkYlPL+UuL36UD0/ORDYNtJwBJUl6eFAol\nPJVs/YBsGWqzFjIhI0G9sbFRXu+l/WI9Ho9qamqStvF6vWpsbEzp2qwyzW7/kVhSOJzwcHu7Pa3T\n4SQpFEpvf+kUMkMKRiKdpo7GY1o7gKEkWZiP7fPew535/IhbgUi7ItHkm4rn2Z2S0hfeYE2BgNR8\n9tLid+kQPz3f2Wzv13P6hiFF4oZrx5X1e4svEoamwbYTgCTpssu6rTnZ+gHZYo4bJ1tzc7bL6JVs\nf7lgqcXkolmaQl5cXNzn9zYcDtmz8BdoRu19u9DlUtTlktF8vtNhn88u3/H0fqgqKzcVtadYZ4JM\nb+QZstts3Z7vTnfPltsj9lj4NqOmzocCnaaOxmNaO4Bc09OdeUdwmAKRQEp37+Ofo49n2IweA3/H\nNt31xyyAoSV+en50uF3H+/FFwIgRpurrO38W6biyfm9d/CLhwhcA/f/clOiLhHi9/WIhlT57o7v+\nkn2JkqzmdNfHlyeDw2B7HELK/pcLGQnqHo9HDQ2XnvVtaGiQx+NJqU04HO7xWkk6cOCADhw4EHtd\nXV2d8kp7/V2RDwCs5v9luwAAAIAct3nz5th/T5o0SZMmTUr52j7elu2dK6+8UrW1taqrq1M4HNbu\n3bs1ffr0Tm2mT5+u9957T5J05MgRFRUVqbS0NKVrpQs/eHV1dex/Hf9QAEmMCSTEuEA8xgTiMSYQ\njzGBeIwJxNu8eXOnfNqbkC5l6I66YRi67777Yqu4z507VxUVFXrzzTclSfPmzdO0adO0b98+LVmy\nRC6XSw8++GDSawEAAAAAGIoy9oz61KlTNXXq1E7H5s2b1+n1/fffn/K1AAAAAAAMRRmZ+p4NvZ1a\ngKGPMYFEGBeIx5hAPMYE4jEmEI8xgXj9HRO2aLaWWgcAAAAAAF0M2TvqAAAAAAAMRgR1AAAAAAAs\nJGOLyWXS/v37tWnTptgq8bfddlu2S0KGnT17VmvXrpXP55PNZtPNN9+s7373uzp//ryeffZZnT17\nViNGjNCjjz6qoqKibJeLDDJNU0uXLpXH49HSpUsZEzmutbVV69ev18mTJyVJixcv1qhRoxgTOWzr\n1q3auXOnbDabxo4dq8WLFysQCDAmcsy6deu0b98+ud1urVq1SpKS/r7YunWr3n77bdntdi1atEhT\npkzJZvkYAInGxKuvvqq9e/fK4XDosssu0+LFi1VYWCiJMZELEo2Ji7Zv364//OEP2rhxo4YNGyap\n92NiyN1RN01TGzdu1LJly7R69Wrt2rUr9gEMucPhcOhHP/qRVq9ereXLl+uNN97QyZMntW3bNk2e\nPFnPP/+8rr32Wm3bti3bpSLD/vGPf6iiokI2m02SGBM57uWXX9bUqVP17LPP6plnntGYMWMYEzms\nrq5O//rXv7RixQqtWrVKpmlq165djIkcNGfOHC1btqzTse7GwcmTJ7V7926tXr1ay5Yt04YNG2Sa\nZjbKxgBKNCamTJmiVatWaeXKlRo1apS2bt0qiTGRKxKNCenCDcOPP/5YZWVlsWN9GRNDLqjX1NRo\n5MiRKi8vl8Ph0KxZs7Rnz55sl4UMKy0t1RVXXCFJcrlcGjNmjBobG7Vnzx7Nnj1bkvTtb39bH374\nYRarRKY1NDRo3759mjt3ri6uo8mYyF1tbW06fPiw5s6dK0kyDEOFhYWMiRxWWFgowzAUCAQUiUQU\nCATk8XgYEzmoqqqqy6yJ7sbBhx9+qFmzZsnhcKi8vFwjR45UTU1NxmvGwEo0JiZPniy7/UKcmjBh\nghoaGiQxJnJFojEhSa+88oruvffeTsf6MiaG3NT3xsZGeb3e2GuPx8M/jBxXV1enY8eOacKECfL5\nfCotLZUklZSUyOfzZbk6ZNLvf/973XvvvfL7/bFjjIncVVdXJ7fbrXXr1unLL79UZWWlFi5cyJjI\nYcOGDdP3vvc9LV68WE6nU1OmTNHkyZMZE5DU/e+LpqYmTZgwIdbO6/WqsbExKzUie9566y3deOON\nkhgTuezDDz+Ux+PR5Zdf3ul4X8bEkLujDnTU3t6uVatWaeHChSooKOh07uLUZ+SGjz76SG63W5WV\nlepuV0rGRG6JRCL64osvdMstt2jFihVyuVxdpjQzJnJLbW2tduzYobVr1+rFF19Ue3u73nvvvU5t\nGBOQeh4HjJPcsmXLFjkcjlhQT4QxMfQFAgFt3bpV1dXVsWPJdkLvaUwMuTvqHo8nNu1EujDV1ePx\nZLEiZEs4HNaqVat000036YYbbpB04Rvwc+fOqbS0VE1NTSopKclylciUTz/9VB999JH27dunUCgk\nv9+vF154gTGRw7xerzwej8aPHy9JmjlzprZu3arS0lLGRI76/PPPdfXVV6u4uFiSNGPGDB05coQx\nAUndf4bgs2due+edd7Rv3z499dRTsWOMidx05swZ1dfX67HHHpN0Yab30qVLtXz58j6NiSF3R/3K\nK69UbW2t6urqFA6HtXv3bk2fPj3bZSHDotGo1q9frzFjxmjBggWx49OnT9c777wjSXr33Xd1/fXX\nZ6lCZNo999yj3/72t1q7dq0eeeQRTZo0SUuWLGFM5LDS0lKVlZXp1KlTkqSPP/5YX/va13Tdddcx\nJnLU6NGj9dlnnykYDCoajerjjz9WRUUFYwKSuv8MMX36dO3atUvhcFh1dXWqra2NfQGIoW3//v36\n+9//rscee0xOpzN2nDGRm8aOHavf/e53Wrt2rdauXSuPx6MVK1aotLS0T2PCFk12P36Q2rdvX6ft\n2W6//fZsl4QMO3z4sJ5++mmNHTs2Nq3knnvu0fjx49liBzp48KC2b9+uJ554gu3ZctyxY8f04osv\nKhwOx7bWMU2TMZHDXn/9db377ruy2WyqrKzUT37yE7W3tzMmcsxzzz2nQ4cOqbm5WaWlpaqurtb1\n11/f7TjYsmWL3n77bRmGoYULF+ob3/hGln8CpFv8mLjrrru0bds2hcPh2PZbV111lR544AFJjIlc\ncHFMtLS0qKSkRNXV1ZozZ07s/MMPP6xf/epXsfHR2zExJIM6AAAAAACD1ZCb+g4AAAAAwGBGUAcA\nAAAAwEII6gAAAAAAWAhBHQAAAAAACyGoAwAAAABgIQR1AAAAAAAshKAOAABS9tBDD+mTTz5JS1+/\n/OUv9d5776WlLwAAhhJHtgsAACBXPfTQQ/L5fLLbL3xvbrPZ9Pzzz6u0tDTLlSVns9kSHl+7dq12\n7dolh8Mhh8OhcePG6b777tPo0aMTtn/yyScHskwAAAYtgjoAAFm0dOlSXXvttX2+PhKJyDCMNFbU\n935tNptuvfVW3X333QoGg1q/fr3Wrl2r5cuXd2oXjUZj7QEAQFcEdQAALCYUCumPf/yjPvjgA0nS\nN7/5Td17771yOBw6cOCAXnjhBc2fP187duzQ5MmTdfbsWc2fP18zZszQ4cOH9fTTT+uJJ57QtGnT\n9Mknn+jVV1/Vr3/9a9XW1urFF1/U8ePHJUlTpkzRAw88oMLCQkkX7vDfcsst2rlzp06fPq1XX31V\n77//vv7yl78oEAhowYIFKf8MTqdTN954o5577jlJ0s9//nNdc801OnDggI4dO6aVK1dq/fr1uumm\nmzR37lxJ0j//+U/t2LFDjY2N8nq9WrJkiSorK9XY2KiXXnpJhw8flsvl0oIFCzR//vx0/pEDAGAp\nBHUAALLo4t3ljrZs2aKamhqtXLlSkrRy5Ur97W9/09133y1J8vl8am1t1bp162SaprZt26YDBw5o\nxowZOnjwoMrLy3Xo0CFNmzZNBw8e1MSJE2N9/+AHP1BVVZXa2tq0atUqbd68WQsXLoyd3717t5Yt\nW6bi4mKdOnVKGzZs0LJlyzR+/Hj96U9/UmNjY0o/T3t7u3bu3KnKysrYuZ07d+rJJ5/U6NGjFY1G\nO91R/+CDD/Taa6/p8ccf17hx41RbWyuHwyHTNLVixQrdcMMNevTRR9XQ0KBf/OIXGj16tKZMmdL7\nP3AAAAYBFpMDACCLVq5cqUWLFmnRokV65plnJEnvv/++7rzzTrndbrndbt15552dFl2z2Wyqrq6W\nw+GQ0+nUxIkTdejQIUnS4cOHdfvtt+vgwYOS1Cmojxw5Ul//+tflcDjkdru1YMGC2HUXzZ8/Xx6P\nR3l5efr3v/+t6667Ttdcc40cDofuvvvupNPVo9Gotm/frkWLFumnP/2pAoGAHnrooVjNs2fPVkVF\nhex2e5dp9W+99ZZuvfVWjRs3LlZrWVmZjh49qpaWFt1xxx0yDEPl5eWaO3eudu3a1Z8/dgAALI07\n6gAAZNHjjz/e5Rn1pqYmlZWVxV6XlZWpqakp9trtdsvhuPQrfMKECTp16pR8Pp+OHTumxx9/XJs3\nb1ZLS4uOHj0aC+rnzp3Tpk2bdPjwYfn9fkWjUQ0bNqzTe3u93k51dHydn5+v4uLibn8Wm82m73//\n+7E7//E6/kzxGhoadNlll3U5Xl9fr6amJi1atCh2zDRNVVVVddsXAACDHUEdAACLGT58uOrr61VR\nUSFJOnv2rIYPHx47H39XOz8/X+PGjdOOHTs0duxYORwOXXXVVdq+fbtGjhwZC+N//vOfZbfbtWrV\nKhUVFek///mPXn755U59dex7+PDhOnnyZOx1IBBQS0tL0toTTeVPhdfrVW1tbZfjZWVlKi8v1/PP\nP9+nfgEAGIyY+g4AgMXMmjVLW7ZsUXNzs5qbm/Xaa6/ppptuSnrNxIkT9cYbb8Tunk+aNKnTa+nC\nc+P5+fkqKChQY2Ojtm/fnrTPmTNnau/evTp8+LDC4bD++te/Jg3ifQ3pknTzzTdr+/bt+vzzzxWN\nRlVbW6uzZ89q/Pjxcrlcev311xUMBmWapo4fP66jR4/2+b0AALA67qgDAGAxd9xxh/x+vx577DFJ\nF1Z9v+OOO5JeM3HiRG3bti02Jbyqqkrt7e2dpojfdddd+s1vfqOFCxdq1KhR+ta3vqUdO3Z022dF\nRYXuv/9+rVmzJrbqe8ep8PFsNluft1ybOXOmWlpatGbNGjU2Nqq8vFwPP/ywysrKtHTpUr3yyit6\n+OGHFQqFNGbMGP3whz/s0/sAADAY2KL9+fobAAAAAACkFVPfAQAAAACwEII6AAAAAAAWQlAHAAAA\nAMBCCOoAAAAAAFgIQR0AAAAAAAshqAMAAAAAYCEEdQAAAAAALISgDgAAAACAhRDUAQAAAACwkP8P\nzJ0IKaCGvjgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10efc34d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Example of simulation and ploting of 'tradeable' Forward-Prices \n",
    "from matplotlib import rcParams\n",
    "from matplotlib.pyplot import *\n",
    "import matplotlib.pyplot as plt\n",
    "rcParams['figure.figsize'] = (14, 8)\n",
    "rcParams['figure.dpi'] = 150\n",
    "plt.style.use('ggplot')\n",
    "\n",
    "num_simul = 1000000\n",
    "seed = 12345\n",
    "num_bins = 60\n",
    "\n",
    "forwardsimul1 = MCSimulation_2F(40, .15, .05, 50, .05, 30, 0, 30, 60, num_simul, seed)\n",
    "forwardsimul2 = MCSimulation_2F(40, .25, .15, 50, .35, 30, 0, 30, 60, num_simul, seed)\n",
    "forwardsimul3 = MCSimulation_2F(40, .35, .25, 50, .65, 30, 0, 30, 60, num_simul, seed)\n",
    "forwardsimul4 = MCSimulation_2F(40, .45, .35, 50, .95, 30, 0, 30, 60, num_simul, seed)\n",
    "\n",
    "plt.hist(forwardsimul1, num_bins, normed=1, facecolor='green', alpha=0.6)\n",
    "plt.hist(forwardsimul2, num_bins, normed=1, facecolor='blue', alpha=0.5)\n",
    "plt.hist(forwardsimul3, num_bins, normed=1, facecolor='red', alpha=0.4)\n",
    "plt.hist(forwardsimul4, num_bins, normed=1, facecolor='black', alpha=0.3)\n",
    "\n",
    "plt.xlabel('Forward Price')\n",
    "plt.ylabel('Frequency')\n",
    "plt.axis([0, 140, 0, 0.065])\n",
    "plt.title(r'Histogram of %.2E simulation of a Forward Price' %num_simul)\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Example of European Call Option valuation using Monte Carlo - 2 factor model\n",
    "\n",
    "def Call_PriceMC(ForwardPrice, strike, interest_rate,\n",
    "                STvol, LTvol, meanReverting, corr, \n",
    "                maturitydate, valuationdate, \n",
    "                deliveryStart, deliveryEnd, \n",
    "                size = 1000, seed = None):\n",
    "    try: \n",
    "        ForwardPrice = float(ForwardPrice)\n",
    "        LTvol = float(LTvol)\n",
    "        STvol = float(STvol)\n",
    "        if meanReverting != 0:\n",
    "            meanReverting = float(meanReverting)\n",
    "        if corr != 0 :\n",
    "            corr = float(corr)                                                                                                \n",
    "        if LTvol < 0 or STvol < 0 or meanReverting < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if deliveryStart > deliveryEnd or valuationdate > maturitydate :\n",
    "            raise ValueError('Incorrect dates')\n",
    "        if corr > 1 or corr < -1 :\n",
    "            raise ValueError('Wrong correlation')\n",
    "    except: \n",
    "        print('Error passing inputs to Call_PriceMC')\n",
    "        \n",
    "    TimeToMaturity = (maturitydate - valuationdate) / 365.\n",
    "    \n",
    "    forwardsimul = MCSimulation_2F (ForwardPrice, STvol, LTvol, \n",
    "                                   meanReverting, corr, \n",
    "                                   maturitydate, valuationdate, \n",
    "                                   deliveryStart, deliveryEnd, \n",
    "                                   size, seed)\n",
    "    \n",
    "    payoff = np.maximum(forwardsimul - strike, np.zeros(size, dtype=np.float64))   \n",
    "    return (np.exp(-1.0 * interest_rate * TimeToMaturity) * np.sum(payoff) / float(size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Call price: 0.747535349135\n"
     ]
    }
   ],
   "source": [
    "print 'Call price:', Call_PriceMC(40, 40, 0.00, .25, .15, 40, .75, 30/365., 0, 30, 60, 1000000, 12345)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Term-Structure, Intantaneous and Equivalent Volatility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As opposition with the simpler Geometric Brownian Motion widely used in simple financial models, the 2-factor model exhibits a volatility term-structure: this mean that the actual volatility of a particular point in the forward curve depends on the valuation date $t$ and the delivery time $T$ . The volatility term-structure gives how the price behaves depending on how far we are from delivery. It can be seen in different ways."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section describe how the volatility term-structure given by the instantaneous model specified in the previous section, behaves. As starting point, we can redefine the model to obtain a new expression for the volatilities grouped in a new definition, the **instataneous total volatility**\n",
    "\n",
    "\\begin{equation} \n",
    "\\frac{dF(t,T)}{F(t,T)}=\\sigma_{S}(t,T)e^{-a(T-t)}dB_{S,t}+\\sigma_{L}(t,T)dB_{L,t} \\approx \\Sigma_{inst}(t,T)dZ_t\n",
    "\\end{equation}\n",
    "\n",
    "where $dZ_t$ is a new Brownian motion, and the instanteous total volatility expression can be obtained from the variance of the previous SDE,\n",
    "\n",
    "\\begin{equation} \n",
    "\\Sigma_{inst}(t,T)=\\sqrt{Var\\left\\{\\frac{dF(t,T)}{F(t,T)}\\right\\}}=\\sqrt{\\sigma_S^2e^{-2a(T-t)}+\\sigma_L^2+2\\rho\\sigma_S\\sigma_Le^{-a(T-t)}}\n",
    "\\end{equation}\n",
    "\n",
    "Noting that for energy markets, where the short-term volatility is usually much higher that the long-term\n",
    "one, previous equation exhibits price characteristics actually observable on the market: The closer we get to delivery, the more volatile the price is, and if our product is far from delivery it follows a geometric brownian motion.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def VolInst2F( STvol, LTvol, meanReverting, \n",
    "              valuationdate, maturitydate, deliveryStart, deliveryEnd, corr): \n",
    "\n",
    "    try: \n",
    "        LTvol = float(LTvol)\n",
    "        STvol = float(STvol)\n",
    "        meanReverting = float (meanReverting)\n",
    "        corr = float(corr)\n",
    "        if LTvol < 0 or STvol < 0 or meanReverting < 0 or valuationdate > maturitydate:\n",
    "            raise ValueError('Negative inputs.')\n",
    "        if corr > 1 or corr < -1 :\n",
    "            raise ValueError('Incorrect Correlation')\n",
    "    except: \n",
    "        print('Error passing vol inputs')\n",
    "     \n",
    "    TimeToMaturity = (maturitydate - valuationdate) / 365.   \n",
    "    \n",
    "    if TimeToMaturity == 0 :\n",
    "        SquareVolInst = 0\n",
    "    elif meanReverting == 0 : \n",
    "        SquareVolInst = LTvol ** 2 + STvol ** 2 + 2 * corr * LTvol * STvol\n",
    "    else:   \n",
    "        SquareVolInst = (LTvol ** 2 \n",
    "                         + STvol ** 2 * \n",
    "                         (np.exp(- 2. * meanReverting * (deliveryStart-valuationdate)/365.))\n",
    "                         /(2. * meanReverting * TimeToMaturity) \n",
    "                         + 2. * STvol * LTvol * corr * \n",
    "                         (np.exp(-meanReverting*(deliveryStart-maturitydate)/365.)\n",
    "                         - np.exp(-meanReverting * (deliveryStart-valuationdate)/365.))\n",
    "                         /float(meanReverting * TimeToMaturity))\n",
    "        \n",
    "    return np.sqrt(SquareVolInst)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, we will call the **equivalent volatility** of the forward contract for delivery at time $T$ to the average volatility $\\Sigma_{eqv}(t_0, t_m, T )$ of the log-return of $F(t,T)$ between 2 dates $t_0$ (the quotation date) and $t_m$ ($m$\n",
    "for maturity).\n",
    "\n",
    "\\begin{equation} \n",
    "\\Sigma_{eq}(t_0,t_m,T)=\\sqrt{\\frac{V(t_0,t_m,T)}{t_m-t_0}}=\\sqrt{ \\frac{1}{t_m-t_0} \\int_{t_0}^{t_m} \\left( \\sigma_S^2e^{-2a(T-t)}+\\sigma_L^2+2\\rho\\sigma_S\\sigma_Le^{-a(T-t)} \\right) dt}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def volEqM2F( STvol, LTvol, meanReverting, \n",
    "             valuationdate, maturitydate, \n",
    "             deliveryStart, deliveryEnd, corr, \n",
    "             blockSize = 0, granularity = 1 / 365.): \n",
    "    try: \n",
    "        LTvol = float(LTvol)\n",
    "        STvol = float(STvol)\n",
    "        meanReverting = float (meanReverting)\n",
    "        corr = float(corr)\n",
    "        if LTvol < 0 or STvol < 0 or meanReverting < 0 :\n",
    "            raise ValueError('Negative inputs')\n",
    "        if deliveryStart > deliveryEnd or valuationdate > maturitydate:\n",
    "            raise ValueError('Wrong dates')\n",
    "        if corr > 1 or corr < -1 :\n",
    "            raise ValueError('Incorrect Correlation')\n",
    "    except: \n",
    "        print('Error passing vol inputs')\n",
    "        \n",
    "    TimeToMaturity = float((maturitydate - valuationdate) / 365.)\n",
    "    psi = PsiFunction(meanReverting, deliveryStart, deliveryEnd)\n",
    "    \n",
    "    if TimeToMaturity == 0 :\n",
    "        SquareVolEq = 0\n",
    "    elif meanReverting == 0 : \n",
    "        SquareVolEq = LTvol ** 2 + STvol ** 2 + 2 * corr * LTvol * STvol\n",
    "    else:\n",
    "        SquareVolEq = (LTvol ** 2 + \n",
    "                      STvol ** 2 * psi ** 2 * \n",
    "                       (np.exp(- 2 * meanReverting * (deliveryStart - maturitydate) / 365.) \n",
    "                        - np.exp(- 2 * meanReverting * (deliveryStart - valuationdate) / 365.)) \n",
    "                       / ( 2 * meanReverting * TimeToMaturity) + \n",
    "                       2 * STvol * LTvol * corr * psi * \n",
    "                       (np.exp( - meanReverting * (deliveryStart - maturitydate)/365.) - \n",
    "                        np.exp( - meanReverting * (deliveryStart - valuationdate)/365.)) \n",
    "                       / (meanReverting * TimeToMaturity))\n",
    "                            \n",
    "        VForward = np.sqrt(SquareVolEq)\n",
    "        \n",
    "    return VForward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Equivalent volatility: 0.167739441383\n"
     ]
    }
   ],
   "source": [
    "volEq = volEqM2F( .25, .15, 40, 0, 30, 30, 60, .75)\n",
    "print 'Equivalent volatility:', volEq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Call price using Black76 Model: 0.767321908898\n",
      "Call price using Monte-Carlo Two-Factor Model: 0.747535349135\n"
     ]
    }
   ],
   "source": [
    "Call_Black76 = Black(40, 40, 30/365., 0., volEq, 1)\n",
    "Call_MC_2_Factors = Call_PriceMC(40, 40, 0., .25, .15, 40, .75, 30/365., 0, 30, 60, 1000000, 12345)\n",
    "\n",
    "print 'Call price using Black76 Model:', Call_Black76.Price\n",
    "print 'Call price using Monte-Carlo Two-Factor Model:', Call_MC_2_Factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Load this snippet before to run the previous cell\n",
    "\n",
    "\n",
    "class Black(object):\n",
    "    \n",
    "    ''' Class for valuation of European Options (Calls/P.uts) using Black-76 model.\n",
    "    F0 : float : initial forward price/index level\n",
    "    K : float : strike price\n",
    "    t : float : maturity (in year fractions)\n",
    "    r : float : constant risk-free short rate\n",
    "    vol : float : volatility factor in diffusion term (std)\n",
    "    CallPut : integer : 1 for a Call, and -1 for a Put\n",
    "    '''\n",
    "    \n",
    "    def __init__(self, F0, K, t, r, vol, CallPut):\n",
    "        try: \n",
    "            self.F0 = float(F0)\n",
    "            self.K = float(K)\n",
    "            self.t = float (t)\n",
    "            self.r = float(r)\n",
    "            self.vol = float(vol)\n",
    "            self.CallPut = int(CallPut)           \n",
    "            if t < 0 or r < 0 or vol < 0:\n",
    "                raise ValueError('Negative inputs.')\n",
    "            if CallPut != 1 and CallPut != -1:\n",
    "                raise ValueError('For a Call: CallPut=1, or -1 for a Put')\n",
    "        except: \n",
    "            print('Error passing option inputs to Black76 function')\n",
    "        \n",
    "    @property    \n",
    "    def Price(self):\n",
    "        d1 = np.log(self.F0 / self.K) / (self.vol * np.sqrt(self.t)) \\\n",
    "            + 0.5 * self.vol * np.sqrt(self.t)\n",
    "        d2 = d1 - self.vol * np.sqrt(self.t)\n",
    "        Price = np.exp(-self.r * self.t) * self.CallPut * \\\n",
    "                    (self.F0 * st.norm(0, 1).cdf(self.CallPut * d1) \\\n",
    "                    - self.K * st.norm(0, 1).cdf(self.CallPut * d2))\n",
    "        return Price\n",
    "    \n",
    "    @property    \n",
    "    def Delta(self):\n",
    "        d1 = np.log(self.F0 / self.K) / (self.vol * np.sqrt(self.t))\\\n",
    "            + 0.5 * self.vol * np.sqrt(self.t)\n",
    "        Delta = np.exp(-self.r * self.t) * self.CallPut * st.norm(0, 1).cdf(self.CallPut * d1)\n",
    "        return Delta\n",
    "    \n",
    "    @property\n",
    "    def Gamma(self):\n",
    "        d1 = np.log(self.F0 / self.K) / (self.vol * np.sqrt(self.t))\\\n",
    "            + 0.5 * self.vol * np.sqrt(self.t)\n",
    "        Gamma = np.exp(-self.r * self.t) * st.norm(0, 1).pdf(d1)\\\n",
    "            / (self.F0 * self.vol * np.sqrt(self.t))\n",
    "        return Gamma\n",
    "\n",
    "    @property\n",
    "    def Vega(self):\n",
    "        d1 = np.log(self.F0 / self.K) / (self.vol * np.sqrt(self.t)) \\\n",
    "            + 0.5 * self.vol * np.sqrt(self.t)\n",
    "        Vega = np.exp(-self.r * self.t) * self.F0 \\\n",
    "            * np.sqrt(self.t) * st.norm(0, 1).pdf(d1)\n",
    "        return Vega\n",
    "    \n",
    "    @property\n",
    "    def Rho(self):\n",
    "        Rho = -self.t * self.Price()\n",
    "        return Rho\n",
    "\n",
    "    @property\n",
    "    def Theta(self):\n",
    "        d1 = np.log(self.F0 / self.K) / (self.vol * np.sqrt(self.t))\\\n",
    "            + 0.5 * self.vol * np.sqrt(self.t)\n",
    "        d2 = d1 - self.vol * np.sqrt(self.t)\n",
    "        Theta = np.exp(-self.r * self.t) * \\\n",
    "                (-self.F0 * st.norm(0, 1).pdf(d1) * self.vol / 2. / np.sqrt(self.t) + \\\n",
    "                self.CallPut * self.r * (self.F0 * st.norm(0, 1).cdf(self.CallPut * d1) - \\\n",
    "                               self.K * st.norm(0, 1).cdf(self.CallPut * d2)))\n",
    "        return Theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}