Understanding code for pricing Asian Option

Question

I am looking at the code here for Asian option pricing:

https://github.com/saulwiggin/finance-with-python/blob/master/Monte%20Carlo%20and%20Pricing%20Exotic%20Options/asian-option.py.

Here is the code:

import scipy as sp
s0=40. #today stock price
x=40. #excercise price
T=0.5 #maturity in years
r=0.05 #risk-free rate
sigma=0.2 # volatility
n_simulation =100 # number of simulations
n_steps=100
dt=T/n_steps
call=sp.zeros([n_simulation],dtype=float)
for j in range(0, n_simulation):
    sT*=s0
    total=0
for i in range(0,int(n_steps)):
        e=sp.random.normal()
        sT*=sp.exp((r-0.5*sigma*sigma)*dt+sigma*e*sp.sqrt(dt))
        total+=sT
        price_average=total/n_steps
        call[j]=max(price_average-x,0)
        call_price=mean(call)*exp(-r*T)
        print 'call price = ', round(call_price,3)

What is the sT variable defined as here? I don't see a definition. What is the for loop for j in range(0, n_simulation) doing exactly? From context it seems like it should be the stock price at time T. If so how should sT be defined?

Thanks.

EDIT: I understand the code is illegal/won't compile, I'm just asking if anyone can help fill in the blanks

Answer 1

Here's an attempt at figuring out what was intended. Obviously some guesswork and assumptions are involved in this.

The changes involved are:

1. Fixing the line you questioned from sT*=s0 to sT = s0 (i.e. at the start of each simulation, reset the current stock price back to today's price)

2. Fixing some indentation so the loops make sense. (Do j simulations; in each simulation evolve the price for i timesteps and calculate the option payoff; finally average the option payoff over all simulations.)

3. Fixing the line that draws values from the normal distribution.

4. Fixing some imports (to be fair, these might have worked fine in earlier versions of scipy).

5. Fixing the print to be Python 3 compatible.

This now runs in Python 3 and prices the call at about $1.50 - $1.60. (Running more simulations would get you tighter bounds on it.)

import scipy as sp
import math
import statistics

s0=40. #today stock price
x=40. #excercise price
T=0.5 #maturity in years
r=0.05 #risk-free rate
sigma=0.2 # volatility
n_simulation =100 # number of simulations
n_steps=100
dt=T/n_steps
call=sp.zeros([n_simulation],dtype=float)
for j in range(0, n_simulation):
    sT=s0  # Reset the stock price to its t=0 value for a new simulation
    total=0  # Reset the total for the average price calculation
    for i in range(0,int(n_steps)):
        e, = sp.stats.norm().rvs(1)  # Get one value from normal dist
        sT *= math.exp((r-0.5*sigma*sigma)*dt + sigma*e*sp.sqrt(dt))
        total += sT
    price_average=total/n_steps
    call[j]=max(price_average-x,0)
call_price = statistics.mean(call)*math.exp(-r*T)
print('call price = ', round(call_price,3))

In general, I'd gently suggest that a resource isn't the best learning vehicle if the code needs this much massaging in order to run and produce sensible results :)

---

As an alternative, using numpy array operations the same algorithm could be done faster and without explicit for loops:

import numpy as np

spot = 40
strike = 40
T = 0.5
r = 0.05
sigma = 0.2
n_sims = 100000
n_steps = 100

dt = T / n_steps

# Perform geometric Brownian motion:
# Make a matrix price_factors representing the ratio of each price
# to the previous timestep. We treat rows as representing simulations,
# columns as representing timesteps
rands = np.random.normal(size=(n_sims, n_steps))
price_factors = np.exp((r - sigma * sigma / 2) * dt + sigma * rands * np.sqrt(dt))

# Now calculate the prices using the cumulative product (along the time axis)
# of price_factors
prices = spot * np.cumprod(price_factors, axis=1)

# Calculate the Asian call payoff per simulation
payoffs = np.maximum(0, np.mean(prices, axis=1) - strike)

# Final result is the average over simulations
call_price = np.mean(payoffs)
print(call_price)