#*****************************************************************************
# Name: Peyton Capristo
# Pledge: Peyton Capristo...
#
# This program will perform a Monte-Carlo Simulation of a European call option
# given a certain set of inputs and produce a graph with different outcomes. 
#
# Inputs: So = Current Stock price 
#         K = Strike price
#         r = risk free interest rate
#         sigma = volatility of the asset 
#         mu = expected return
#         T = Time to expiration 
#
# Outputs: C = price of black scholes call option
#          S[i] = path of the option price
#
#******************************************************************************

# Import Libraries 
import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt
from scipy import stats 

# Initialize a random seed to reproduce data
np.random.seed(1)

# Now need to define the Black Scholes Metron Call option Pricing Model

def BSM_Call (So, K, r, sigma, T):
    d1 = (np.log(So/K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    C = (So * stats.norm.cdf(d1)) - (K * np.exp(-r * T) * stats.norm.cdf(d2))
    return C

def BSM_Put (So, K, r, sigma, T):
    d1 = (np.log(So/K) + (r + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    P = (K * np.exp(-r * T) * stats.norm.cdf(-d2)) - (So * stats.norm.cdf(-d1))
    return P

# Define Inputs 
So = 95
K = 100
r = 0.1
sigma = 0.2
mu = r # Make assumption that the expected return from the underlying asset is the risk-free interest rate
T = 0.5


# Number of simulation is equal to 100000
N = 100000

# Number of trading days in a year =  252
NumberOfSteps = int(T * 252)
dt = T / NumberOfSteps

# Initialize an array S and its starting point to be equal to the current stock option price
S = np.zeros((NumberOfSteps, N))
S[0] = So

# This for loop runs the simulation starting at the first point until the expiration date for an option path
# It does this for the number of simulations N desired
for i in range (1, NumberOfSteps):
    Z = np.random.standard_normal(N) 
    S[i] = S[i-1] * np.exp((mu - 0.5 * sigma ** 2) * dt + (sigma * np.sqrt(dt) * Z))

# Averages the Final S[i] values
Avg_Si = sum(S[i]) / len(S[i])

# Finds highest and lowest ending value
Max_Si = np.amax(S[i])
Min_Si = np.amin(S[i])

# Collects all the final data points and calculates the BSM value
All_Final_Call_Prices = BSM_Call(S[i], K, r, sigma, T)

# Averages the final BSM values 
Avg_Final_Call_Price = sum(All_Final_Call_Prices) / len(All_Final_Call_Prices)

# Sums and Discounts all the values
Sum_and_Discounted_Call_Price = np.exp(-r * T) * 1 / N * np.sum(np.maximum(S[-1] - K, 0))

# Collects all the final data points and calculates the BSM value
All_Final_Put_Prices = BSM_Put(S[i], K, r, sigma, T)

# Averages the final BSM values 
Avg_Final_Put_Price = sum(All_Final_Put_Prices) / len(All_Final_Put_Prices)

# Sums and Discounts all the values
Sum_and_Discounted_Put_Price = 1 / N * np.sum(np.maximum(K - S[-1], 0)) * np.exp(-r * T) 


# print("All_Final_Call_Price = ", All_Final_Call_Price) #If want list of all final uncomment out
print("Average Final Stock Values =", Avg_Si)
print("Maximum Ending Stock Value =", Max_Si)
print("Minimum Ending Stock Value =", Min_Si)
print("\nInitial BSM call price =", BSM_Call(So, K, r, sigma, T))
print("Average Final Call Price =", Avg_Final_Call_Price)
print("Sum and discounted Call Price =", Sum_and_Discounted_Call_Price)
print("\nInitial BSM put price =", BSM_Put(So, K, r, sigma, T))
print("Average Final Put Price =", Avg_Final_Put_Price)
print("Sum and discounted Put Price =", Sum_and_Discounted_Put_Price)


# Set ups and displays the graph
plt.plot(S[:,0:100000]) # Prepare a  plot of option price movements
plt.axhline(K, c="k", xmin=0,      # Plot a Horizontal line at the strike price
            xmax=NumberOfSteps,
           label="Strike Price")
plt.title('$25 Differnce Between So and K') # Give a title for the graph
plt.xlabel('Time in Days') # Give x axis label for the graph
plt.ylabel('Option Price') # Give y axis label for the graph
plt.grid(True, which='both')  # Show grid lines in both directions

plt.show() # Finally, show the plot