text: """
  Created on Fri Feb 18 09:50:30 2022
  
  @author: udaytalwar
  """
  
  '''
  References:
  
  ASX Portfolio on YouTube
  Paul Wilmott's "The Mathematics of Financial Derivatives: A Student Introduction"
  
  '''
  
  import scipy.stats as stats
  import math as m 
  import numpy as np
  
  #calculating option price using Black-Scholes 
  
  def d(S, K, sigma, r, t):
      
      '''
      S = Current Price
      K = Strike Price
      sigma = Volatility
      r = annualized risk-free rate
      t = time to expiration
      
      returns d1, d2 for option price calculation using Black Scholes
      '''
      d1 = (m.log(S/K) + (r + (sigma**2/2))*t) * (1/(sigma*m.sqrt(t)))
      
      d2 = d1 - sigma*m.sqrt(t)
      
      return d1, d2
  
  def option_price(S, K, sigma, r, t, flag, d1 = 0, d2 = 0):
      
      '''
      S = Current Price
      K = Strike Price
      sigma = Volatility
      r = annualized risk-free rate
      t = time to expiration
      flag = 'Call' or 'Put'
      
      returns option price according to Black Scholes
      '''
      
      if d1 == 0 and d2 == 0:
      
          d1, d2 = d(S, K, sigma, r, t)
          
          if flag == 'Call':
              
              price = stats.norm.cdf(d1)*S - stats.norm.cdf(d2)*K*m.exp(-r*t)
              
          elif flag == 'Put':
              
              price = stats.norm.cdf(-d2)*K*m.exp(-r*t) - stats.norm.cdf(-d1)*S
              
          return price 
  
      else: 
          
          if flag == 'Call':
              
              price = stats.norm.cdf(d1)*S - stats.norm.cdf(d2)*K*m.exp(-r*t)
              
          elif flag == 'Put':
              
              price = stats.norm.cdf(-d2)*K*m.exp(-r*t) - stats.norm.cdf(-d1)*S
              
          return price  
      
  #calculating option price using Monte Carlo Approximation
  
  def option_price_MC(S, K, sigma, r, t, flag, steps, paths, error = False):
  
      '''
      S = Current Price
      K = Strike Price
      sigma = Volatility
      r = annualized risk-free rate
      t = time to expiration
      flag = 'Call' or 'Put'
      steps = number of steps per path 
      paths = number of paths to simulate
      error = False, if true returns (price, approximation error)
      
      returns expected option price from simulation
  
      '''
          
      timestep = t/steps 
      
      nudt = (r - 0.5*sigma**2)*timestep
      
      volsdt = vol*m.sqrt(timestep)
      
      lnS = m.log(S)
      
      random_matr = np.random.normal(size=(steps,paths))
  
      delta_lnSt = nudt + volsdt*random_matr
      
      lnSt = lnS + np.cumsum(delta_lnSt, axis = 0)
      
      lnSt = np.concatenate( (np.full(shape = (1,paths), fill_value = lnS), lnSt))
      
      ST = np.exp(lnSt)
      
      if flag == 'Call':    
          price = np.maximum(0, ST - K)
      
      elif flag == 'Put':
          price = np.maximum(0, K - ST)
      
      initial_price = np.exp(-r*t)*np.sum(price[-1])/paths
      
      err = np.sqrt(np.sum((price[-1]-initial_price)**2)/(paths-1))
      
      SE = err/m.sqrt(paths)
      
      if error == True: 
          
          return initial_price, SE
  
      else:
          return initial_price
  
  #calculating option price using Binomial Trees
  
  def option_price_BT(S, K, sigma, r, t, flag, steps):
  
      '''
      S = Current Price
      K = Strike Price
      sigma = Volatility
      r = annualized risk-free rate
      t = time to expiration
      flag = 'Call' or 'Put'
      steps = number of steps per path 
      
      returns option price using Binomial Trees
      '''
  
      
      timestep = t/steps #equivalent to dt, a discrete time step to iterate through
      
      
      #the formulas for A, d, u and p are described in Wilmott's book 
      
      A = (1.0/2.0)*(m.exp(-r*timestep)+m.exp((r+sigma**2)*timestep))
      
      d = A - m.sqrt(A**2 - 1)
      
      u = A + m.sqrt(A**2 - 1)
      
      p = (m.exp(r*timestep)-d)/(u-d)
  
      #discount function to discount the final value back to our initial value ie current option price
      discount = m.exp(-r*timestep)
  
      #price at node (i,j) at timestep i is = start price * d^(i-j)*u^(i) where j is the number of the node at the 
      # given time step, where we count from the bottom of the tree up to the top
  
      price = S * d ** (np.arange(steps,-1,-1)) * u **(np.arange(0, steps+1, 1))
      
          #price - K for Call, K - price for put
      if flag == 'Call':
      
          price = np.maximum(price - K, np.zeros(steps+1))
      
      elif flag == 'Put':
          
          price = np.maximum(K - price, np.zeros(steps+1))
      
          #apply the discount function to iteratively arrive at initial value
          
      for j in np.arange(steps, 0, -1):
          
          # formula below also given in Wilmott's book 
          
          price = discount * (p * price[1:j+1] + (1-p) * price[0:j])
      
      return price[0]
      
  
  t = 30.0/365.0
  S0 = 100
  K = 100
  flag = 'Put'
  vol = (0.4/30.0)*np.sqrt(365)
  r = 0.015
  
  steps = 20
  paths = 10000
  steps_bt = 500
  
  print('Price according to B-S = ',str(round(option_price(S0, K, vol, r, t, flag), 3)))
  
  print('Price according to Monte Carlo = ',str(round(option_price_MC(S0, K, vol, r, t, flag, steps, paths), 3)))
  
  print('Price according to Binomial Tree= ',str(round(option_price_BT(S0, K, vol, r, t, flag, steps_bt), 3)))
  