0
$\begingroup$

It's winter break (happy new year!), and I'm trying implement a few options pricing models (bapm, tapm, monte carlo, Fast Fourier etc.) for practice.

The issue: My BAPM CRR model converges to 8.45544504853379 for "ec", which is consistent with online results. My Monte Carlo implementation converges to 4.748831589623564 :(. I can't find the error in my code...I was hoping someone could tell me if such a large error is possible when using ~100k GBM simulations.

type = ["ec", "ep", "ac", "ap"] (Euro call, euro put, etc.)

Binomial Asset Pricing Model

# Initial parameters
S0 = 100      # initial stock price
K = 100       # strike price
T = 1         # time to maturity in years
r = 0.06      # annual risk-free rate
N = 100         # number of time steps
sig = .13      # sigma - std of stock annualized 

import numpy as np
import matplotlib.pyplot as plt

class BAPM():
    def __init__(self, S0, K, T, N, r, c, sig) -> None:
        self.S0 = S0
        self.K = K
        self.T = T
        self.N = N
        self.r = r
        self.c = c
        self.sig = sig
        self.dt = self.T/self.N

    def JarrowRudd(self, type):
        u = np.exp((self.r - self.sig**2 / 2)*self.dt + self.sig*np.sqrt(self.dt))
        d = np.exp((self.r - self.sig**2 / 2)*self.dt - self.sig*np.sqrt(self.dt))
        p = (np.exp(self.r*self.dt) - d)/(u - d)
        q = 1-p
        disc = np.exp(-self.r*self.dt)
        return self.optionPrice(type, p, q, u, d, disc)

    def CRR(self, type, drift = 0):
        u = np.exp(drift*self.dt + self.sig * np.sqrt(self.dt))
        d = np.exp(drift*self.dt - self.sig * np.sqrt(self.dt))
        p = (np.exp(self.r*self.dt) - d)/(u - d)
        q = 1-p
        disc = np.exp(-self.r*self.dt)
        return self.optionPrice(type, p, q, u, d, disc)

    def Tian(self, type):
        v = np.exp(self.sig **2 * self.dt)
        u = .5 * np.exp(self.r * self.dt) * v * (v + 1 + np.sqrt(v**2 + 2*v -3))
        d = .5 * np.exp(self.r * self.dt) * v * (v + 1 - np.sqrt(v**2 + 2*v -3))
        p = (np.exp(self.r * self.dt) - d)/(u - d)
        q = 1-p
        disc = np.exp(-self.r * self.dt)
        return self.optionPrice(type, p, q, u, d, disc)

    def optionPrice(self, type, p, q, u, d, disc):
        # European Option
        if type == "ec" or type == "ep":
            # terminal value of stock
            S = self.S0 * (d ** (np.arange(self.N,-1,-1))) * (u ** (np.arange(0,self.N+1,1)))
            V = 0
            # value of option at terminal time
            if type == "ec":
                V = np.maximum(S - self.K, np.zeros(self.N+1))
            elif type == "ep":
                V = np.maximum(self.K - S, np.zeros(self.N+1))
            # take 2 at a time and get weighted-discounted value
            for i in np.arange(self.N,0,-1):
                V = disc * ( p * V[1:i+1] + q * V[0:i])
            # final result is option price
            return V[0]

        elif type == "ac" or type == "ap":
            print("test successful")
            # terminal value of stock
            S = self.S0 * (d**(np.arange(self.N,-1,-1))) * (u ** (np.arange(0,self.N+1,1)))

            V = 0
            # value of payoff at terminal time
            if type == "ac":
                V = np.maximum(0, S - self.K)
            elif type == "ap":
                V = np.maximum(0, self.K - S)

            # backtrack through options tree
            for i in np.arange(self.N-1, -1, -1):
                # recalculate prices at current level
                S = self.S0 * d**(np.arange(i,-1,-1)) * u**(np.arange(0,i+1,1))

                V[:i+1] = disc * ( p*V[1:i+2] + q*V[0:i+1] )                
                V = V[:-1]

                if type == 'ap':
                    V = np.maximum(V, self.K - S)
                else:
                    V = np.maximum(V, S - self.K)
  
            return V[0]

bapm = BAPM(S0, K, T, N, r, c, sig)
print(bapm.CRR(type = "ec"))

8.45544504853379

Monte Carlo Pricing (vanilla only...American not implemented yet)

import numpy as np
import matplotlib.pyplot as plt

class MCPricing():
    def __init__(self, S0, drift, sig, T, tIncrement, r) -> None:
        self.S0 = S0
        self.mu = drift
        self.sig = sig
        self.T = T
        self.dt = tIncrement
        self.r = r

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

    def simulateGBM(self):
        W = np.random.standard_normal(size = self.N)
        W = np.cumsum(W)*np.sqrt(self.dt)
        X = (self.mu - .5*self.sig**2)*self.t + self.sig*W
        S = self.S0*np.exp(X)
        return S.tolist()

    def calculateOptionPrice(self, K, nsim, type):
        if type == "ec" or type == "ep":
            sims = [self.simulateGBM() for i in range(nsim)]
            payoffs = []
            for path in sims:
                if type == "ec":
                    payoffs.append(np.maximum(path[-1] - K, 0))
                elif type == "ep":
                    payoffs.append(np.maximum(K - path[-1], 0))
            payoffs = np.array(payoffs)
            return np.average(np.exp(-1*self.r * self.T) * payoffs)

mcpm = MCPricing(S0 = 100, drift=0, sig =.13, T=1, tIncrement=.01, r=.06)
path1 = mcpm.simulateGBM()
plt.plot(mcpm.t, path1)

mcpm.calculateOptionPrice(100, 100000, "ec")

4.748831589623564

I'd appreciate any help! Thanks

$\endgroup$
3
  • $\begingroup$ Your mc drift is scaled by self.t, whereas I think dt would be more appropriate, given the sqrt dt scaling for the Brownian $\endgroup$ Dec 31, 2021 at 21:43
  • $\begingroup$ @JamesSpencer-Lavan I made that change, and the MC is now consistent at $5.30. This is much better but is still far off from the BAPM. I tried decreasing the dt but get the same result. $\endgroup$ Dec 31, 2021 at 22:23
  • 1
    $\begingroup$ Your drift term is wrong. We're pricing the option in a risk-neutral world meaning that the drift needs to be equal to r, not 0. Try setting it to .06 and you will probably get the correct answer. $\endgroup$
    – Oscar
    Jan 1 at 21:19

1 Answer 1

0
$\begingroup$

Really wish I knew why my previous version didn't work...but I found some modifications to GBM that got me the answer.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

class MCPricing():
    def __init__(self, S0, drift, sig, T, tIncrement, r) -> None:
        self.S0 = S0
        self.mu = drift
        self.sig = sig
        self.T = T
        self.dt = tIncrement
        self.r = r

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

    def simulateGBM(self, nsims):
        Xt = np.log(self.S0) +\
             np.cumsum(((self.mu - self.sig**2/2)*self.dt + \
             self.sig*np.sqrt(self.dt) * \
             np.random.normal(size=(self.N,nsims))),axis=0)
             
        return np.exp(Xt)

    def calculateOptionPrice(self, K, nsims, type):
        if type == "ec" or type == "ep":
            paths = self.simulateGBM(nsims)
            if type == "ec":
                payoffs = np.maximum(paths[-1]-K, 0)
            elif type == "ep":
                payoffs = np.maximum(K - paths[-1], 0)
                
            option_price = np.mean(payoffs)*np.exp(-self.r * self.T) #discounting back to present value
            
            return option_price

mcpm = MCPricing(S0 = 100, drift=0.06, sig =.13, T=1, tIncrement=.001, r=.06)
paths = mcpm.simulateGBM(1000)
plt.plot(paths)
mcpm.calculateOptionPrice(100, 10000, "ec")

8.373535799654391

$\endgroup$
2
  • 1
    $\begingroup$ It works now because you set drift and r to be equal, which is exactly what I told you to do in a comment a few days ago. $\endgroup$
    – Oscar
    Jan 5 at 21:31
  • $\begingroup$ @Oscar yeah I made that change in the previous implementation but was still getting an incorrect answer. But yeah, the mismatched rates was a stupid mistake. Thank you for the help. $\endgroup$ Jan 7 at 7:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.