Barrier Option Pricing using Python

Question

This is my implementation of pricing an exotic option (in this case an up-and-in barrier option) using the Monte Carlo simulation in Python. I use NumPy where I can. Any ideas to optimize this code?

################################################
# Exotic Barrier Option (up-and-in) Pricer 
# By: DudeWah
################################################
#-----------------------------------------------
################################################
# Monte Carlo Simulation pricing an up-and-in  # 
# call option. This call option is a barrier   #
# option in which pyoffs are zero unless the   #
# asset crosses some predifned barrier at some #
# time in [0,T]. If the barrier is crossed,    #
# the payoff becomes that of a European call.  #
# Note: Monte Carlo tends to overestimate the  #
# price of an option. Same fo Barrier Options. #
################################################
#-----------------------------------------------

#import libraries
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import style
from random import gauss
from math import exp, sqrt

#-----------------------------------------------
#Class for all parameters
#-----------------------------------------------

class parameters():
    """parameters to be used for program"""

    #--------------------------
    def __init__(self):
        """initializa parameters"""

        self.S = 5                     #underlying asset price
        self.v = 0.30                  #volatility
        self.r = 0.05                  #10 year risk free rate
        self.T = 365.0/365.0           #years until maturity
        self.K = 6                     #strike price
        self.B = 8                     #barrier price
        self.delta_t = .001            #timestep
        self.N = self.T/self.delta_t   #Number of discrete time points
        self.simulations = 1000        #num of simulations

    #---------------------------
    def print_parameters(self):
        """print parameters"""

        print "---------------------------------------------"
        print "---------------------------------------------"
        print "Pricing an up-and-in option"
        print "---------------------------------------------"
        print "Parameters of Barrier Option Pricer:"
        print "---------------------------------------------"
        print "Underlying Asset Price = ", self.S
        print "Volatility = ", self.v
        print "Risk-Free 10 Year Treasury Rate =", self.r
        print "Years Until Expiration = ", self.T
        print "Time-Step = ", self.delta_t
        print "Discrete time points =", self.N
        print "Number of Simulations = ", self.simulations
        print "---------------------------------------------"
        print "---------------------------------------------"


#-----------------------------------------------
#Class for Monte Carlo
#-----------------------------------------------    

class up_and_in_mc(parameters):
    """This is the class to price the barrier option 
       defined as an up and in option. Paramaters are
       fed in as a subclass"""

    #---------------------------   
    def __init__(self):
        """initialize the class including the subclass"""
        parameters.__init__(self)
        self.payoffs = []
        self.price_trajectories = []
        self.discount_factor = exp(-self.r * self.T)


    #---------------------------
    def call_payoff(self,s):
        """use to price a call"""
        self.cp = max(s - self.K, 0.0)
        return self.cp



    #---------------------------
    def calculate_payoff_vector(self):
        """Main iterative loop to run the 
        Monte Carlo simulation. Returns a vector
        of different potential payoffs"""

        for i in xrange(0, self.simulations):
            self.stock_path = []
            self.S_j = self.S
            for j in xrange(0, int(self.N-1)):
                self.xi = gauss(0,1.0)

                self.S_j *= (exp((self.r-.5*self.v*self.v) * self.delta_t + self.v *sqrt(self.delta_t) * self.xi))

                self.stock_path.append(self.S_j)
            self.price_trajectories.append(self.stock_path)
            if max(self.stock_path) > self.B:
                self.payoffs.append(self.call_payoff(self.stock_path[-1]))
            elif max(self.stock_path) < self.B:
                self.payoffs.append(0)

        return self.payoffs           

    #---------------------------
    def compute_price(self):
        """Uses payoff vector and discount
        factor to compute the price of the
        option. Numpy used for efficiency"""
        self.np_payoffs = np.array(self.payoffs, dtype=float) 
        self.np_Vi = self.discount_factor*self.np_payoffs 
        self.price = np.average(self.np_Vi)

    #---------------------------
    def print_price(self):
        """prints the option price to terminal"""
        print str("Call Price: %.4f") % self.price
        print "---------------------------------------------"

    #---------------------------
    def calc_statistics(self):
        """uses payoffs and price to calc
        variance, standard deviation, and
        a 95% confidence interval for price"""
        self.variance = np.var(self.np_Vi)
        self.sd = np.std(self.np_Vi)
        #95% C.I. uses 1.96 z-value
        self.CI = [self.price - (1.96*self.sd/sqrt(float(self.simulations))),
                   self.price + (1.96*self.sd/sqrt(float(self.simulations)))]

    #---------------------------
    def print_statistics(self):
        """prints option statistics that were
        calculated to the terminal"""
        print "Variance: %.4f" % self.variance
        print "Standard Deviation: %.4f" % self.sd
        print "95% Confidence Interval:", self.CI
        print "---------------------------------------------"
        print "---------------------------------------------\n"

    #---------------------------
    def plot_trajectories(self):
    """uses matplotlib to plot the trajectory
    of each individual stock stored in earlier
    trajectory array"""
    print "Creating Plot..."
    #use numpy to plot 
    self.np_price_trajectories = np.array(self.price_trajectories, dtype=float)
    self.times = np.linspace(0, self.T, self.N-1)

    #style/plot/barrier line
    style.use('dark_background')

    self.fig = plt.figure()
    self.ax1 = plt.subplot2grid((1,1),(0,0))
    for sublist in self.np_price_trajectories:
        if max(sublist) > self.B:
            self.ax1.plot(self.times,sublist,color = 'cyan')
        else:
            self.ax1.plot(self.times,sublist,color = '#e2fb86')
    plt.axhline(y=8,xmin=0,xmax=self.T,linewidth=2, color = 'red', label = 'Barrier')

    #rotate and add grid
    for label in self.ax1.xaxis.get_ticklabels():
        label.set_rotation(45)
    self.ax1.grid(True)

    #plotting stuff
    plt.xticks(np.arange(0, self.T+self.delta_t, .1))
    plt.suptitle('Stock Price Trajectory', fontsize=40)
    plt.legend()
    self.leg = plt.legend(loc= 2)
    self.leg.get_frame().set_alpha(0.4)
    plt.xlabel('Time (in years)', fontsize = 30)
    plt.ylabel('Price', fontsize= 30)
    plt.show()

#-----------------------------------------------
#Main Program
#-----------------------------------------------

#Initialize and print parameters
prm = parameters()
prm.print_parameters()

#Price/print the option
ui_mc = up_and_in_mc()
ui_mc.calculate_payoff_vector()
ui_mc.compute_price()
ui_mc.print_price()

#caclulate/print stats
ui_mc.calc_statistics()
ui_mc.print_statistics()

#plot
ui_mc.plot_trajectories()

Answer 1

Numerics

There is a minor advantage to initializing S to 5. instead of 5 as static typing systems will recognize it as a float, which is closer to reality. Same for K, B.

If T is one year, then just... write 1..

Use the Numpy functions for exp, sqrt instead of their math alternatives.

The for j loop needs to go away and be vectorised. The for i loop can optionally be vectorised, though this is less important and takes more memory. If you're serious about this being a library, then you should add a partial vectorisation over the i (simulation) dimension that calculates the simulations in batches of perhaps 100 at a time.

The offset (self.r-.5*self.v*self.v) * self.delta_t and coefficient self.v *sqrt(self.delta_t) should be calculated outside of any loop.

Don't represent price_trajectories as a list; pre-allocate it as an ndarray via np.zeros. Don't allocate stock_path at all; reference it as a slice out of the trajectories matrix.

You can fully remove

        elif max(self.stock_path) < self.B:
            self.payoffs.append(0)

because you can preallocate to zeros as above.

Don't iteratively self.S_j *=. Instead call into the ufunc np.multiply.accumulate.

Don't self.np_payoffs = np.array(self.payoffs, dtype=float). Start with an array and don't use a list.

Prefer self.vi.mean() over average(). Same for var and std.

You write:

# 95% C.I. uses 1.96 z-value

but... that's hard-coded, when it shouldn't be; and it's inaccurate. Call into scipy.stats.norm.interval instead.

Data vis

Remove all of the "---------------------------------------------"; they're just noise.

For the Risk-Free 10 Year Treasury Rate, show it as a percentage via {:%}.

You need to separate out your semantic plot calls from your styling plot calls. Also, I just don't think most of the style calls are a good idea, particularly the giant fonts.

Add the payout and non-payout entries to the legend.

Don't y=8; instead use your class constant self.B.

You definitely shouldn't plot all of your trajectories. That creates a giant, uninformative smear. Because this is a Monte Carlo process, sort the trajectories by some indicator like the final price, and then sample from them evenly to display a small (20?) number of representative curves.

Class structure

It's not at all helpful that parameters (which should be written Parameters) is a parent class. If you really need the parameters to be, well, parametric, then you should pass in an instance. I demonstrate a simpler version that just has one class with these parameters set as class (not instance) variables.

call_payout should be a @classmethod.

You have a habit of tossing nearly every variable that should be local-only into class state (e.g. xi). Don't do this; it's healthy for variables to be forgotten.

When you instantiate your class, do so in a main() protected by a __main__ guard rather than in the global namespace.

Demonstration

"""
Exotic Barrier Option (up-and-in) Pricer

Monte Carlo Simulation pricing an up-and-in call option. This call option is a barrier
option in which payoffs are zero unless the asset crosses some predefined barrier at some
time in [0,T]. If the barrier is crossed, the payoff becomes that of a European call.
Note: Monte Carlo tends to overestimate the price of an option. Same fo Barrier Options.
"""

import numpy as np
from matplotlib import pyplot as plt
import scipy.stats


class UpAndInMc:
    __slots__ = (
        'payoffs', 'price_trajectories', 'vi', 'rand',
        'price', 'variance', 'sd', 'CI',
    )

    S = 5.    # underlying asset price
    v = 0.30  # volatility
    r = 0.05  # 10 year risk-free rate
    T = 1.    # years until maturity
    K = 6.    # strike price
    B = 8.    # barrier price
    delta_t = .001       # timestep
    confidence = 0.95    # confidence interval around price
    simulations = 1_000  # num of simulations
    N = int(T/delta_t)   # Number of discrete time points
    discount_factor = np.exp(-r*T)

    def __init__(self, seed: None | int = None) -> None:
        self.payoffs = np.zeros(self.simulations)
        self.price_trajectories = np.empty((self.simulations, self.N - 1))
        self.vi = np.empty(self.simulations)
        self.CI = np.empty(2)
        self.rand = np.random.default_rng(seed=seed)
        self.price: float | None = None
        self.variance: float | None = None
        self.sd: float | None = None

    @classmethod
    def format_params(cls) -> str:
        return (
f'''Underlying asset price ${cls.S:.2f}
Volatility = {cls.v}
Risk-free 10 year treasury rate = {cls.r:.1%}
Years until expiration = {cls.T}
Time step = {cls.delta_t}
Discrete time points = {cls.N}
Number of simulations = {cls.simulations}''')

    @classmethod
    def call_payoff(cls, s: float) -> float:
        """price a call"""
        return max(s - cls.K, 0.)

    def calculate_payoff_vector(self) -> None:
        """
        Main iterative loop to run the Monte Carlo simulation. Returns a vector
        of different potential payoffs
        """
        offset = (self.r - 0.5*self.v**2)*self.delta_t
        coef = self.v*np.sqrt(self.delta_t)

        for i in range(self.simulations):
            xi = self.rand.normal(loc=0., scale=1., size=self.N - 1)
            prod_series = np.multiply.accumulate(np.exp(offset + coef*xi))
            self.price_trajectories[i] = self.S*prod_series
            if self.price_trajectories[i].max() > self.B:
                self.payoffs[i] = self.call_payoff(self.price_trajectories[i,-1])

    def compute_price(self) -> None:
        """
        Uses payoff vector and discount factor to compute the price of the
        option. Numpy used for efficiency
        """
        self.vi[:] = self.discount_factor*self.payoffs
        self.price = self.vi.mean()

    def print_price(self) -> None:
        """prints the option price to terminal"""
        print(f'Call price: ${self.price:.4f}')

    def calc_statistics(self) -> None:
        """
        uses payoffs and price to calc
        variance, standard deviation, and
        a 95% confidence interval for price
        """
        self.variance = self.vi.var()
        self.sd = self.vi.std()
        interval = scipy.stats.norm.interval(confidence=self.confidence, loc=self.price, scale=self.sd)
        self.CI[:] = (interval - self.price)/np.sqrt(self.simulations) + self.price

    def print_statistics(self) -> None:
        print(f'Variance: {self.variance:.4f}')
        print(f'Standard deviation: {self.sd:.4f}')
        a, b = self.CI
        print(f'{self.confidence:.0%} confidence interval: {a:.4f}, {b:.4f}')

    def plot_trajectories(self) -> None:
        """
        uses matplotlib to plot the trajectory of each individual stock stored in earlier
        trajectory array
        """
        times = np.linspace(start=0, stop=self.T, num=self.N - 1)

        fig, ax = plt.subplots()
        n_curves = 20
        order = self.price_trajectories[:, -1].argsort()
        prices = self.price_trajectories[order]
        idx = np.concatenate((
            np.arange(0, self.simulations, self.simulations//n_curves),
            (-1,),
        ))
        labelled = set()
        for price_series in prices[idx]:
            if price_series.max() > self.B:
                color = 'orange'
                name = 'payoff'
            else:
                color = 'blue'
                name = 'no payoff'
            if name in labelled:
                label = None
            else:
                label = name
                labelled.add(name)
            ax.plot(times, price_series, color=color, linewidth=0.5, label=label)

        ax.axhline(y=self.B, xmin=0, xmax=self.T, color='red', label='Barrier')
        ax.grid(True)
        ax.legend()
        ax.set_title('Stock Price Trajectory')
        ax.set_xlabel('Time (years)')
        ax.set_ylabel('Price')


def main() -> None:
    ui_mc = UpAndInMc(seed=0)
    print(ui_mc.format_params())
    ui_mc.calculate_payoff_vector()
    ui_mc.compute_price()
    ui_mc.print_price()
    ui_mc.calc_statistics()
    ui_mc.print_statistics()

    # Regression test for seed=0
    assert np.isclose(a=0.25416602693109610, atol=0, rtol=1e-15, b=ui_mc.price)
    assert np.isclose(a=0.66033430179464100, atol=0, rtol=1e-15, b=ui_mc.variance)
    assert np.isclose(a=0.81260956294806240, atol=0, rtol=1e-15, b=ui_mc.sd)
    assert np.isclose(a=0.20380088989925702, atol=0, rtol=1e-15, b=ui_mc.CI[0])
    assert np.isclose(a=0.30453116396293500, atol=0, rtol=1e-15, b=ui_mc.CI[1])

    ui_mc.plot_trajectories()
    plt.show()


if __name__ == '__main__':
    main()

Underlying asset price $5.00
Volatility = 0.3
Risk-free 10 year treasury rate = 5.0%
Years until expiration = 1.0
Time step = 0.001
Discrete time points = 1000
Number of simulations = 1000
Call price: $0.2542
Variance: 0.6603
Standard deviation: 0.8126
95% confidence interval: 0.2038, 0.3045