#!/usr/bin/env python
# coding: utf-8

# <img src="https://raw.githubusercontent.com/israeldi/quantlab/master/assets/images/Program-Logo.png" width="400px" align="right">
# 
# # QuantLab: Stochastics
# ### [(Go to Quant Lab)](https://israeldi.github.io/quantlab/)
# 
# &copy; Dr. Yves J. Hilpisch | The Python Quants GmbH
# 
# <img src="http://hilpisch.com/images/py4fi_2nd_shadow.png" width="200px" align="left">

# Initially import all the modules we will be using for our notebook

# In[1]:


import math
import numpy as np
import numpy.random as npr 

# COMMAND PROMPT: pip install matplotlib
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import os


# ### 4.3 Pricing European Options by Monte Carlo
# 
# The payoff of a European call option on an index at maturity is given by $h(S_T)\equiv max\{S_T-K, 0\}$, where $S_T$ is the index level at maturity date $T$ and $K$ is the strike price. Given a risk-neutral measure for the relevant stochastic process (e.g., geometric Brownian motion), the price of such an option is given by the formula:
# 
# $$C_{0}=e^{-rT}\mathbb{E}_{0}^{Q}[h(S_{T})]=e^{-rT}\intop_{0}^{\infty}h(s)q(s)ds$$
# 
# The equation below provides the respective Monte Carlo estimator for the European option, where $\tilde{S}_{T}^{i}$ is the $T$-th simulated index level at maturity.
# 
# $$\tilde{C}_{0}=e^{-rT}\frac{1}{I}\sum_{i=1}^{I}h(\tilde{S}_{T}^{i})$$
# 
# Consider the following parameterization for the geometric Brownian motion and the valuation function `gbm_mcs_stat()`, taking as a parameter only the strike price. Here, only the index level at maturity is simulated. As a reference, consider the case with a strike price of $K = 105$:

# In[35]:


def bsm_call_value(S0, K, T, r, sigma):
    ''' Valuation of European call option in BSM model.
    Analytical formula.

    Parameters
    ==========
    S0: float
        initial stock/index level
    K: float
        strike price
    T: float
        maturity date (in year fractions)
    r: float
        constant risk-free short rate
    sigma: float
        volatility factor in diffusion term

    Returns
    =======
    value: float
        present value of the European call option
    '''
    from math import log, sqrt, exp
    from scipy import stats

    S0 = float(S0)
    d1 = (log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
    d2 = (log(S0 / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
    # stats.norm.cdf --> cumulative distribution function
    #                    for normal distribution
    value = (S0 * stats.norm.cdf(d1, 0.0, 1.0) -
             K * exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0))
    return value


# In[36]:


def gen_sn(M, I):
    ''' Function to generate random numbers for simulation.
    
    Parameters
    ==========
    M: int
        number of time intervals for discretization
    I: int
        number of paths to be simulated
    '''
    sn = npr.standard_normal((M + 1, I))
    return sn


# In[58]:


S0 = 100.
r = 0.05
sigma = 0.25
T = 1.0
I = 50000
K = 105


# In[59]:


def gbm_mcs_stat(K):
    ''' Valuation of European call option in Black-Scholes-Merton
    by Monte Carlo simulation (of index level at maturity)
    
    Parameters
    ==========
    K: float
        (positive) strike price of the option
    
    Returns
    =======
    C0: float
        estimated present value of European call option
    '''
    sn = gen_sn(1, I)
    
    # Simulate Prices at Maturity
    
    # calculate payoff at maturity
    
    # calculate discounted price

    return C0


# In[60]:


# The Monte Carlo estimator value for the European call option.
gbm_mcs_stat(K)


# In[61]:


# The analytical solution is
bsm_call_value(S0, K, T, r, sigma)


# Next, consider the dynamic simulation approach and allow for European put options in addition to the call option. The function `gbm_mcs_dyna()` implements the algorithm. The code also compares option price estimates for a call and a put stroke at the same level:

# In[63]:


def gbm_mcs_dyna(K):
    ''' Valuation of European options in Black-Scholes-Merton
    by Monte Carlo simulation (of index level paths)
    
    Parameters
    ==========
    K: float
        (positive) strike price of the option
    option : string
        type of the option to be valued ('call', 'put')
    
    Returns
    =======
    C0: float
        estimated present value of European call option
    '''
    dt = T / M
    
    # simulation of index level paths
    S = np.zeros((M + 1, I))
    S[0] = S0
    sn = gen_sn(M, I)
    
    for t in range(1, M + 1):
        S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt 
                + sigma * math.sqrt(dt) * sn[t])
    
    # calculation of payoff
    hT = np.maximum(S[-1] - K, 0)
    
    # calculation of MCS estimator
    C0 = math.exp(-r * T) * np.mean(hT)
    
    return C0


# In[49]:


# The Monte Carlo estimator value for the European call option (Dynamic Simulation Approach)
gbm_mcs_dyna(K)  


# ### Euler Discretization
# 
# Consider the stochastic differential equation $\mathrm{d} X_{t}=a\left(X_{t}\right) \mathrm{d} t+b\left(X_{t}\right) \mathrm{d} W_{t}$ with initial condition $X_{0}=x_{0},$ where $W_{t}$ stands for the Wiener process, and suppose that we wish to solve this SDE on some internal of time $[0, T]$. Then the Euler-Maruyama approximation to the true solution $X$ is the Markov chain $Y$ defined as follows:
# 
# - partition the interval $[0,7 \text { into } N \text { equal subintervals of width } \Delta t>0$ :
# $$
# 0=\tau_{0}<\tau_{1}<\cdots<\tau_{N}=T \text { and } \Delta t=T / N
# $$
# - set $Y_{0}=x_{0}$
# - recursively define $Y_{n}$ for $1 \leq n \leq N$ by
# $$
# \begin{array}{l}
# {Y_{n+1}=Y_{n}+a\left(Y_{n}\right) \Delta t+b\left(Y_{n}\right) \Delta W_{n}} \\
# {\Delta W_{n}=W_{\tau_{n+1}}-W_{\tau_{n}}}
# \end{array}
# $$

# In[108]:


# Calculate Price with these parameters again
S0 = 100.
r = 0.05
sigma = 0.25
T = 1.0
I = 5000
K = 105


# In[109]:





# In[110]:





# In[111]:





# ### Finite Difference Scheme
# (Next Time!)