Monte Carlo simulation_Implied volatilities_Black–Scholes–Merton (BSM)_Geometric_Geometric Brownian

import numpy as np
from math import log
import matplotlib.pyplot as plt

x = []
y = []

# Equation that defines profit
def generateProfit(d):
    global s
    
    if d>=s :
        return 0.6*s
    else:
        return 0.6*d - 0.4*(s-d)

# Although y comes from uniform distribution in [80, 140]
# we are running simulation for d in [20,305]

maxprofit = 0

for s in range(20,305): # s:the number of units sold
    # Run a simulation for n =1000
    # Even if we run for n=10,000 the result would be almost the same
    for i in range(1,1000):
        ###### generate a random value of d ######
        d = np.random.randint(low=10, high=200)
        
        # for this random value of d（demand）, find profit and update maxprofit
        profit = generateProfit(d)
        if profit>maxprofit:
            maxprofit = profit
            
    # store the value of s to be plotted along X axis
    x.append(s)

    # store the value of maxprofit plotted along Y axis
    y.append(log(maxprofit))#plotted on log scale

plt.plot(x,y)

plt.title('Find MaxProfit by Monte Carlo Simulation')
plt.xlabel('Units Sold')
plt.ylabel('Max Profits (log scale)')
print('Max Profit:', round(maxprofit,1))

plt.annotate('log({})={}'.format( round(maxprofit,1), round( log(maxprofit),3) ), 
             ha='left', va='top', 
             xytext = (200, 4.), #The xytext parameter specifies the text position.
             xy=( 200, round( log(maxprofit),3) ),     #The xy parameter specifies the arrow's start position
             arrowprops={'facecolor':'blue',  'arrowstyle':'<-'})

plt.show()

# Will display this Max Profit: 119.4

"in a classroom full of 30 students, what is the probability that more than one person will have the same birthday?" We will assume that it is not a leap year and instead of using the month and day, we will just use the days in the calendar year, that is, 365 days. The logic is pretty simple. The following code shows how one can calculate the probability of more than one person having the same birthday in a room of 30 students: 1 - ( 365 364 363 ... 336 /365^30 ) = 1－0.2937=0.7063

import numpy as np

numstudents =30
numTrials = 10000
count_numWithSameBday =0

for trial in range(numTrials):
    year = [0]*365  #a list with 365 elements (or called days)
    for i in range(numstudents): #30 students' birthday
        newBDay = np.random.randint(365)#############################
        year[newBDay] = year[newBDay] + 1
    
    haveSameBday = False
    for num in year:
        if num > 1:#############################
            haveSameBday=True
    
    if haveSameBday == True:
        count_numWithSameBday += 1    #############################
    
prob = float(count_numWithSameBday)/float(numTrials)#############################
print('The probability of a shared birthday in a class of ', numstudents, ' is ', prob)

Monte Carlo simulation in basketball¶

"when down by three and left with only 30 seconds is it better to attempt a hard 3-point shot or an easy 2-point shot and get another possession?"

import numpy as np
import matplotlib.pyplot as plt

colors = [  ( 31,119,180), (174, 199, 232), (255,127, 14), (255, 187, 120), 
            ( 44,160, 44), (214,  39,  40), (148,103,189), (152, 223, 138), 
            (255,152,150), (197, 176, 213), (140, 86, 75), (196, 156, 148), 
            (227,119,194), (247, 182, 210), (127,127,127), (199, 199, 199),
            (188,189, 34), (219, 219, 141), (23, 190,207), (158, 218, 229),
            (217,217,217)
         ]

# Scale RGB values to the [0, 1] range, format matplotlib accepts.
for i in range( len(colors) ):
    r, g, b = colors[i]
    colors[i] = (r/255., g/255., b/255.)

def attemptThree():
    if np.random.randint(0, high=100) < threePtPercent:
        if np.random.randint(0, high=100) < overtimePercent:
            return True # We won!!
    return False # we either missed the 3 or lost in OT

Assuming that on an average, it takes only 5 seconds to attempt a two-point shot and the two-point scoring percent labeled twoPtPercent of the player is pretty high, then they score a two-point shot, which will be deducted from the value in the pointsDown variable

def attemptTwo():
    havePossession = True
    pointsDown = 3
    timeLeft = 30
    while (timeLeft > 0):
        # What to do if we have possession
        if havePossession :
            # If we are down by 3 or more, we take 2 quickly.           
            if pointsDown >= 3:
                timeLeft -= timeToShoot2
            else: # if we are down by 2 or less,we run down the clock first
                timeLeft = 0
            
            # Do we make the shot?
            if (np.random.randint(0, high =100) < twoPtPercent):
                pointsDown -=2
                havePossession = False
        else:
            # Does the opponent team rebound?
            # If so, we lose possession.
            #This doesn't really matter when we run the clock down
            if (np.random.randint(0, high=100) >= offenseReboundPercent):
                havePossession = False
            else: #cases where we don't have possession
                if (pointsDown > 0): #foul to get back possession
                    #takes time to foul
                    timeLeft -= timeToFoul
                    
                    #opponent takes 2 free throws
                    if (np.random.randint(0, high=100) < oppFtPercent):
                        pointsDown += 1
                    
                    if (np.random.randint(0, high=100) < oppFtPercent):
                        pointsDown += 1
                        havePossession = True
                else:
                    if (np.random.randint(0, high=100) >= ftReboundPercent):
                        #you were able to rebound the missed it
                        havePossession = true
                    else:
                        #tied or up so don't want to foul;
                        #assume opponent to run out clock and take
                        if (np.random.randint(0, high=100) < oppTwoPtPercent):
                            pointsDown += 2 # They made the 2
                        timeLeft = 0
    
    if pointsDown > 0:
        return False
    else:
        if pointsDown <0:
            return True
        else:
            if np.random.randint(0, high=100) < overtimePercent:
                return True
            else:
                return False

names=['Lebron James', 'Kyrie Irving', 'Steph Curry', 'Kyle Krover', 'Dirk Nowitzki']
threePercents = [35.4, 46.8, 44.3, 49.2, 38.0]
twoPercents =   [53.6, 49.1, 52.8, 47.0, 48.6]

plt.figure(figsize=(14,14))

colind = 0
for i in range(5):
    x = []
    y1 = []
    y2 = []
    
    trials = 400 #Number of trials to run for simulation
    
    threePtPercent = threePercents[i] #Your % chance of making 3-pt shot
    twoPtPercent = twoPercents[i] #Your % chance of making a 2-pt shot
    
    oppTwoPtPercent = 40 #Opponent % chance making 2-pter
    oppFtPercent = 70 #Opponent's FT %    
    
    timeToShoot2 = 5 #How many seconds elapse to shoot a 2
    timeToFoul = 5 #How many seconds elapse to foul opponent
    
    offenseReboundPercent = 25 #% of regular offense rebound
    ftReboundPercent = 15 #% of offense rebound after missed FT
    
    overtimePercent = 50 #% chance of winning in overtime
    
    winsTakingThree = 0
    lossTakingThree = 0
    winsTakingTwo = 0
    lossTakingTwo = 0
    curTrial = 1
    
    #Each time draw() is called, we run another trial and update
    #graphics
    while curTrial < trials:
        #run a trial take 3
        if(attemptThree()):
            winsTakingThree +=1
        else:
            lossTakingThree +=1
            
        #run a trial take a 2
        if attemptTwo() == True:
            winsTakingTwo +=1
        else:
            lossTakingTwo +=1
            
        #print " winsTakingThree=",winsTakingThree," ->percentage",(100.0*winsTakingThree/curTrial),"winsTakingTwo=",winsTakingTwo," ->percentage",(100.0*winsTakingTwo/curTrial)
        x.append(curTrial)
        y1.append(winsTakingThree)
        y2.append(winsTakingTwo)
        curTrial += 1
        
    plt.plot(x,y1, color=colors[colind], label=names[i] + " Wins Taking Three Point")
    plt.plot(x,y2, color=colors[20], label=names[i] + " Wins Taking Two Point")
    colind +=2
    
legend = plt.legend(loc='upper left', shadow=True)
for legobj in legend.legendHandles:
    legobj.set_linewidth(2.4)
plt.show()

   

#import pandas.io.data as stockdata
import pandas_datareader.data as stockdata
import datetime
import numpy as np
import matplotlib.pyplot as plt

#r,g,b=(31,  119, 180)
#colornow=(r/255.,g/255.,b/255.)
ibmquotes = stockdata.DataReader( name='IBM', data_source='yahoo', start='2005-10-1' )
ibmquotes['Volatility'] = np.log( ibmquotes['Close']/ibmquotes['Close'].shift(1) )############################
# today's volatility == today close price/ previous close price

#ibmquotes[['Close', 'Volatility']].plot(figsize=(12,10),subplots=True, color=colornow) 
ibmquotes[['Close','Volatility']].plot(figsize=(12,10), subplots=True,color='Red')
plt.show()

                    

Above: The higher the variance in their values(Close), the more volatile it turns out

Volatility is the measure of variation of price, which can have various peaks. Moreover, the exponential function lets us plug in time and gives us growth; logarithm (inverse of exponential) lets us plug in growth and gives us the time measure.

import pandas as pd

url = 'http://www.stoxx.com/download/historical_values/h_vstoxx.txt'
vstoxx_index = pd.read_csv(url, index_col=0, header=2,
                          parse_dates=True, #####################
                          dayfirst=True,#DD/MM format dates, international(YYYY-MM-DD) and European format.
                          sep=',') # header on the 2nd row
vstoxx_index.head()

This data file consists of Euro Stoxx volatility indices, which can all be plotted via one simple filter mechanism of dates between Dec 31, 2011 to May 1, 2015. The following code can be used to plot the VSTOXX data:

vstoxx_short = vstoxx_index[('2011/12/31' < vstoxx_index.index) 
                           & (vstoxx_index.index<'2015/5/1')]
vstoxx_short.head()

# to plot all together
vstoxx_short.plot(figsize=(15,14))

plt.title('it creates a plot that compares Euro Stoxx volatility index:')
plt.show() 

# to plot in subplots separately
vstoxx_short.plot(subplots=True, grid=True, color='r', figsize=(20,20), linewidth=2)

plt.title("plotting in subplots separately")
plt.show()

from math import log, sqrt, exp
from scipy import stats
import pandas as pd
import matplotlib.pyplot as plt

colors = [(31, 119, 180), (174, 199, 232), (255,128,0),
(255, 15, 14), (44, 160, 44), (152, 223, 138), (214, 39, 40),
(255, 152, 150),(148, 103, 189), (197, 176, 213), (140, 86, 75),
(196, 156, 148),(227, 119, 194), (247, 182, 210), (127, 127, 127),
(199, 199, 199),(188, 189, 34), (219, 219, 141), (23, 190, 207),
(158, 218, 229)]
# Scale the RGB values to the [0, 1] range, which is the format matplotlib accepts.
for i in range(len(colors)):
    r, g, b = colors[i]
    colors[i] = (r / 255., g / 255., b / 255.)

######comment begin

######comment end

#S==S0 as t==0      K==X       initial
def black_scholes_merton(S, r, sigma, X,T):
    S = float(S)  #convert to float
    
    d1 =( log(S/X) + ( (r + 0.5* sigma**2)*T) ) / ( sigma*sqrt(T) )
    d2 =( log(S/X) + ( (r - 0.5* sigma**2)*T) ) / ( sigma*sqrt(T) )
    
    # stats.norm.cdf -> cumulative distribution function X~N(0,1)=N(u,sigma**2)
                                #    u     sigma square - variance
    value = ( S * stats.norm.cdf(d1, 0.0, 1.0) - 
             X * exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0))
    
    return value #present value of the European call option

def vega(S, r, sigma, X, T):
    S = float(S)
    
    d1 =( log(S/X) + ( (r + 0.5* sigma**2)*T) ) / ( sigma*sqrt(T) )
    
    vagaValue = S * stats.norm.pdf(d1, 0.0, 1.0) * sqrt(T)
    
    return vagaValue

                            #sigma_est: initialVal  #iteration=100
                            #X: strike price  or exercise prices
                            #T: time to expiration
def impliedVolatility(S, r, sigma_est, X, T, Cstar, it):###########################
    
    for i in range(it):
        numer = (black_scholes_merton(S, r, sigma_est, X, T) - Cstar)
        denom = vega(S,r,sigma_est, X, T)
        sigma_est -= numer/denom
        
    return sigma_est  #sigma_estimated #implied volatilty

Step1

import pandas as pd

h5 = pd.HDFStore('C:/Users/LlQ/Data Visualization/cp5/vstoxx_data_31032014.h5', 'r',encoding='utf-8')

futures_data = h5['futures_data'] # VSTOXX futures data
options_data = h5['options_data'] # VSTOXX call option data

h5.close()

futures_data['DATE']= pd.to_datetime(futures_data['DATE'])
futures_data['MATURITY']= pd.to_datetime(futures_data['MATURITY'])
futures_data

options_data['DATE']= pd.to_datetime(options_data['DATE'])
options_data['MATURITY']= pd.to_datetime(options_data['MATURITY'])

options_data['IMP_VOL'] = 0.0  #initial implied volatility################################
V0 = 17.6639 #the closing value of the index
r=0.04       #risk free interest rate
sigma_est = 2
tol = 0.5    #tolerance level for moneyness

step 2

for option in options_data.index:
    #iterating over all option quotes
    futureVal = futures_data[ futures_data['MATURITY'] == options_data.loc[option]['MATURITY'] ]['PRICE'].values[0]
    
    # picking the right futures values
    if( futureVal * (1-tol) < options_data.loc[option]['STRIKE'] < futureVal * (1+tol) ):
        #impliedVolatility(S, r, sigma_est, X, T, Cstar, it):
        impliedVol = impliedVolatility(V0, 
                                        r, 
                                        sigma_est,      #volatility inital value or first 
                                        options_data.loc[option]['STRIKE'],
                                        options_data.loc[option]['TTM'],
                                        options_data.loc[option]['PRICE'], #option market price/option quote
                                        it=100)#iterations
        options_data['IMP_VOL'].loc[option] = impliedVol################################

options_data.loc[46170]

from datetime import datetime

plot_data = options_data[ options_data['IMP_VOL']>0 ]################################
maturities = sorted( set(options_data['MATURITY']) )
maturities

#plot_data = options_data[ options_data['IMP_VOL']>0 ]
#maturities = sorted( set(options_data['MATURITY']) )

plt.figure(figsize =(15,10))

i=0
for maturity in maturities:
    data = plot_data[ options_data.MATURITY == maturity ]
    
    # select data for this maturity
    plot_args = {'lw':3, 'markersize': 9}
    plt.plot(data['STRIKE'], data['IMP_VOL'], label=maturity.date(), marker='o', color=colors[i], **plot_args)
    i+=1
    
plt.grid(True)
plt.xlabel(r'Strike Price $X$', fontsize=18)
plt.ylabel(r'Implied volatility of $\sigma$', fontsize=18)
plt.title('Short Maturity Window (Volatility Smile)', fontsize=22)

plt.legend()
plt.show()

keep=['PRICE', 'IMP_VOL']

group_data = plot_data.groupby(['MATURITY','STRIKE'])[keep]

group_data=group_data.sum()

#tell us each maturity date and strike price will have unique option  price and implied violatility
group_data 
#The resulting DataFrame object has two index levels and two columns.

Monte Carlo Simulation

Monte Carlo simulation is one of the most important algorithms in finance and numerical science in general. Its importance stems from the fact that it is quite powerful when it comes to option pricing or risk management problems. In comparison to other numerical methods, the Monte Carlo method can easily cope with high-dimensional problems where the complexity and computational demand, respectively, generally increase in linear fashion. The downside of the Monte Carlo method is that it is per se computationally demanding and often needs huge amounts of memory even for quite simple problems. Therefore, it is necessary to implement Monte Carlo algorithms efficiently. The example that follows illustrates different implementation strategies in Python and offers three different implementation approaches for a Monte Carlo-based valuation of a European option.[15] The three approaches are:[16]

Pure Python

This example sticks with the standard library — i.e., those libraries and packages that come with a standard Python installation — and uses only built-in Python capabilities to implement the Monte Carlo valuation.

Vectorized NumPy

This implementation uses the capabilities of NumPy to make the implementation more compact and much faster.

Fully vectorized NumPy

The final example combines a different mathematical formulation with the vectorization capabilities of NumPy to get an even more compact version of the same algorithm.

 

# from the previous example, we can calculate the exact option value

# black_scholes_merton(S, r, sigma, X,T): return value # present value of the European call option (VS quote)
from bsm_functions import black_scholes_merton

S0 = 100.
K = 105.
T = 1.0
r = 0.05
sigma=0.2

black_scholes_merton(S0, r, sigma, K, T)

# Monte Carlo valuation of European call options with pure Python

from time import time
from math import exp, sqrt, log
from random import gauss, seed

#import numpy as np
#np.random.seed(2000)

seed(2000)
t0 =time()

# Parameters
S0 = 100.    # initial value
K = 105.     # strike price
T = 1.0      # maturity
r = 0.05     # riskless short rate
sigma = 0.2   # volatility
M = 50       # number of time steps
dt = T/M     # length of time interval #delta t
I = 250000    # number of paths

#Simulating I paths with M time steps
S=[]
for i in range(I):
    path = []
    for t in range(M+1): #0,1,2--50
        if t==0:
           path.append(S0)
        else:
           Z = gauss(0.0, 1.0) # np.random.standard_normal()#
           St = path[t-1] * exp( (r - 0.5*sigma**2) * dt + sigma * sqrt(dt) * Z )       
           path.append(St)
        
    S.append(path)

# Calculating the Monte Carlo estimator
#Determine the inner value hT of the European call option at T as hT(ST(i)) = max(ST(i) – K,0).
C0 = exp(-r * T) * sum([max(path[-1] - K, 0) for path in S]) / I
                       
# Results output
tpy = time() - t0
print("European Option Value {:7.3f}".format(C0))  #European Option Value   8.011
print("Duration in Seconds {:7.3f}".format(tpy))   #Duration in Seconds  39.034

# C0 = exp(-r * T) * sum([max(path[-1] - K, 0) for path in S]) / I
sum_val = 0.0
for path in S:
    # C-like iteration for comparison
    sum_val += max(path[-1]-K, 0)
C0 = exp(-r * T) * sum_val /I
round(C0, 3)

# Monte Carlo valuation of European call options with Numpy
# mcs_vector_numpy.py

import math
import numpy as np
from time import time

np.random.seed(20000)
t0 = time()

# Parameters
S0 = 100.
K = 105.
T = 1.0
r = 0.05
sigma = 0.2
M = 50
dt= T/M
I = 250000

# Simulating I paths with M time steps
S = np.zeros( (M+1,I) )
S[0] = S0 #the zero -row  with 250,000 columns filled with same value S0=100
for t in range(1, M+1):
    #In the line in question, 250,000 numbers are generated in a single step
    Z = np.random.standard_normal(I) #psedurandom numbers    ###############
    
         # path[t-1] * exp( (r - 0.5*sigma**2) * dt + sigma * sqrt(dt) * Z ) 
    S[t] = S[t-1] * np.exp( (r-0.5*sigma**2) * dt + sigma*math.sqrt(dt)*Z )
    
    # vectorized operation per time step over all paths
    
# Calculating the Monte Carlo estimator
C0 = math.exp(-r * T) * np.sum( np.maximum(S[-1] - K, 0 ))/I

# Results output
tnp1 = time() - t0
print("European Option Value {:7.3f}".format(C0))
print("Duration in Seconds {:7.3f}".format(tnp1))

Vectorization brings a speedup of more than 30 times in comparison to pure Python. The estimated Monte Carlo value is again quite close to the benchmark value(==black_scholes_merton(S0, r, sigma, K, T)==8.021352235143176).

# Monte Carlo valuation of European call options with Numpy (log version)
# mcs_full_vector_numpy.py

import math
from numpy import *
from time import time
# start import for shorter code

random.seed(20000)
t0 = time()

# Parameters
S0 = 100.
K = 105.
T = 1.0
r = 0.05
sigma = 0.2
M = 50
dt = T/M
I = 250000

# S=((M+1, I))=((51, 250000))

# Simulating I paths with M time steps 

#row 0: S0 * exp(...) + S0 *exp(...) + 。。。== S0 * (exp(...）+exp(...))
#S[t] = S[t-1] * np.exp( (r-0.5*sigma**2) * dt + sigma*math.sqrt(dt)*Z )

#以下先将原先的公式 两边求对:  指数移动到等号左边  右边剩下指数部分，然后累加 log S[t] -log S[0] = cumsum(...)
#减数部分再移回右边 然后  两边再求指 还原公式:  S[t] = S[0] * exp(cumsum(...) )  
#这里先用S0 代替S[0], 使得所有的column都一样
S = S0 * exp( cumsum( (r - 0.5*sigma**2)*dt 
                        + sigma * math.sqrt(dt)
                        * random.standard_normal((M+1, I)), axis=0 #sum previous value on row to next row
                    ) 
            )

# sum instead of cumsum would also do
# if only the final values are of interest
S[0] = S0 # the zero row filled with same value S0 then do previous step

# Calculating the Monte Carlo estimator
C0 = math.exp(-r * T) * sum(maximum(S[-1]-K,0))/I

#Results output
tnp2 = time() - t0
print("European Option Value {:7.3f}".format(C0))
print("Duration in Seconds {:7.3f}".format(tnp2))

plt.hist(S[-1], bins=50)

plt.grid(True)
plt.xlabel('index level')
plt.ylabel('frequency')

plt.show()

the output, this time illustrating the (approximately) log-normal distribution of the end-of-period index level values

 

plt.hist(np.maximum(S[-1]-K, 0), bins=50)

plt.grid(True)
plt.xlabel('option inner value')  # S[-1]-K  == option price - strike price
plt.ylabel('frequency')
plt.ylim(0,50000)  # if option price < K, then inner value==0

In this case, the majority of the simluated values are zero, indicating that the European call option expires worthless in a significant amount of cases.

Technical Analysis

In finance, technical analysis is a security analysis methodology for forecasting the direction of prices through the study of past market data, primarily price and volume.

import numpy as np
import pandas as pd
import pandas_datareader.data as web

sp500 = web.DataReader('^GSPC', data_source='yahoo', start='1/1/2000', end='4/14/2014')
sp500.info()

sp500['Close'].plot(grid=True, figsize=(8,5))

sp500.head()

# window: int offset
# Size of the moving window. This is the number of observations 
# used for calculating the statistic. Each window will be a fixed size.

sp500['42d'] = np.round(sp500['Close'].rolling(window=42, center=False).mean(), 2)
sp500['252d'] = np.round(sp500['Close'].rolling(window=252, center=False).mean(), 2)
print(sp500[['Close', '42d', '252d']].tail())

sp500[40:50]

sp500[['Close', '42d', '252d']].plot(grid=True, figsize=(8,5))

As the market Close price went down, the new Close value is lower than original average value, or the average curve is over new Close curve(pull down action)

As the market Close price went up, the new Close value is higher than original average value, or the average curve is under the new Close curve(pull up action)

sp500['42-252'] = sp500['42d'] - sp500['252d']
sp500['42-252'].tail()

SD = 50
sp500['Regime'] = np.where(sp500['42-252'] >  SD,  1, 0)   #1 buy 0 wait
sp500['Regime'] = np.where(sp500['42-252'] < -SD, -1, sp500['Regime'])   # -1 sell; 0 wait
sp500['Regime'].value_counts()

sp500['Regime'].plot(lw=1.5)
plt.ylim([-1.1, 1.1])

Based on the respective regime, the investor either is long or short in the market (index) or parks his wealth in cash, which does not bear any interest. This simplified strategy allows us to work with market returns only. The investor makes the market return when he is long (1), makes the negative market returns when he is short (–1), and makes no returns (0) when he parks his wealth in cash.

Note that the shift method shifts a time series by as many index entries as desired — in our case by one trading day, such that we get daily log returns:

sp500['Market'] = np.log(sp500['Close']/sp500['Close'].shift(1))

Monte Carlo valuation of European call options with pure Python

from time import time from math import exp, sqrt, log from random import gauss, seed

seed(2000) t0 =time()

Parameters

S0 = 100. # initial value K = 105. # strike price T = 1.0 # maturity r = 0.05 # riskless short rate sigma = 0.02 # volatility M = 50 # number of time steps dt = T/M # length of time interval #delta t I = 250000 # number of paths

#Simulating I paths with M time steps S=[] for i in range(I): path = [] for t in range(M+1): #0,1,2--50 if t==0: path.append(S0) else:# #pseudorandom numbers Z = gauss(0.0, 1.0) St = path[t-1] exp((r - 0.5sigma*2) dt + sigma sqrt(dt)Z)
path.append(St)

S.append(path)

Calculating the Monte Carlo estimator

#Determine the inner value hT of the European call option at T as hT(ST(i)) = max(ST(i) – K,0). C0 = exp(-r T) sum([max(path[-1] - K, 0) for path in S]) / I

Results output

tpy = time() - t0 print("European Option Value %7.3f" % C0) print("Duration in Seconds %7.3f" % tpy)

European Option Value 0.858 Duration in Seconds 40.877

The portfolio valuation

The common sense of portfolio valuation of an entity is to estimate its current worth. Valuations are usually applicable on a financial asset or liability and can be performed on stocks, options, business enterprises, or intangible assets.

import pandas as pd # gets numpy as pd.np
from pandas_datareader.data import get_data_yahoo
import matplotlib.pyplot as plt

three funds from the Vanguard, such as Vanguard US Total (vus.to), Vanguard Canadian Capped (vre.to), and Vanguard Emerging Markets (vee.to).

# get data
data = get_data_yahoo(['vus.to', 'vre.to', 'vee.to'], start='2014-01-01') ['Adj Close']
data.head()

data.plot(figsize=(10,10), lw=2)
plt.ylabel('Adj Close')

plt.show()

# convert prices to log returns
retn = data.apply(pd.np.log, axis=0).diff()
retn.head()

# make scatterplot  to show correlation
pd.tools.plotting.scatter_matrix(retn, figsize=(10,10))
plt.show()

# some more stats
retn.skew()
retn.kurt()

import matplotlib.pyplot as plt
import numpy as np

'''
Geometric Brownian Motion with drift!
u=drift factor
sigma: volatility
T: time span
dt: length of steps
S0: Stock Price in t=0
W: Brownian Motion with Drift N[0,1]
'''
rect = [0.1, 5.0, 0.1, 0.1]
fig = plt.figure( figsize=(10,10) )

T = 2                       # T: time span
mu = 0.1
sigma = 0.04                # sigma: volatility
S0 = 20                     # S0: Stock Price in t=0
dt = 0.01                   # dt: length of steps

N = round(T/dt)
t = np.linspace(0, T, N)

# Standare normal distrib
W = np.random.standard_normal( size=N )
W = np.cumsum(W) * np.sqrt(dt)

# exponent
X = (mu-0.5*sigma**2)*t + sigma*W

# Brownian Motion
S = S0 * np.exp(X)

plt.plot(t,S, lw=2)

plt.xlabel("Time t", fontsize=16)
plt.ylabel("S", fontsize=16)
plt.title("Geometric Brownian Motion (Simulation)", fontsize=18)

plt.show()

The simulating stock prices using Brownian motion

import pylab, random

class Stock(object):
    def __init__(self, price, distribution):
        self.price = price
        self.history = [price]
        self.distribution = distribution
        self.lastChange = 0
        
    def setPrice(self, price):
        self.price = price
        self.history.append(price)
    
    def getPrice(self):
        return self.price
    
    def walkIt(self, marketBias, mo):
        oldPrice = self.price
        baseMove = self.distribution() + marketBias
        self.price = self.price * (1.0 + baseMove)
        
        if mo:
            self.price = self.price + random.gauss(.5, .5) * self.lastChange
            
        if self.price < 0.01:
            self.price = 0.0
        
        self.history.append(self.price)
        self.lastChange = oldPrice - self.price
        
    def plotIt(self, figNum):
        pylab.figure(figNum)
        
        pylab.plot(self.history)
        
        pylab.title( 'Closing Price Simulation Run -' + str(figNum) )
        pylab.xlabel('Day')
        pylab.ylabel('Price') 

def testStockSimulation():
    def runSimulation(stocks, fig, mo):
        mean = 0.0
        
        for s in stocks:
            for d in range(numDays):
                s.walkIt(bias, mo)
                
            s.plotIt(fig)
            mean += s.getPrice()
            
        mean = mean/float(numStocks)
        pylab.axhline(mean)
        
    pylab.figure( figsize=(12,12) )
    numStocks = 20
    numDays = 400
    stocks = []
    bias = 0.0
    mo = False
    startvalues = [100,500,200,300,100,100,100,200,200, 300,300,400,500,00,300,100,100,100,200,200,300]
    
    for i in range(numStocks):
        volatility = random.uniform(0,0.2)
        
        dl = lambda: random.uniform(-volatility, volatility)
        stocks.append( Stock(startvalues[i], dl))
    
    runSimulation(stocks, 1, mo)

testStockSimulation()
pylab.show()