This is the third post in a series about Python for Finance. We will work again with the example of the first and second post of the series to illustrate now how to automate the valuation of a complete set of American put options.
In the first two posts, we presented a rather concise and nevertheless quite fast implementation of the Longstaff-Schwartz algorithm (LSM) to value American options by Monte Carlo simulation (MCS). In this post, we will build on the previous implementations and will value all American put options from Table 1 in the seminal paper by Longstaff-Schwartz (2001); see publisher’s Web site or the recent book “Derivatives Analytics with Python” by Dr. Yves J. Hilpisch, MD of Visixion GmbH, for details of the algorithm.
The Python script presented below is still quite compact, even if not as compact as the ones before. This script now also takes care of the PEP 8 Style Guide for Python (cf. PEP 8 Style Guide).
When running the script, you might get results that look like this:
run LSM_Complete_Bench.py
S0 36 | VOL 0.20 | T 1 | Value 4.479 | Benchmark 4.478 | Abs 0.001 | Rel 0.016
S0 36 | VOL 0.20 | T 2 | Value 4.827 | Benchmark 4.840 | Abs -0.013 | Rel -0.265
S0 36 | VOL 0.40 | T 1 | Value 7.098 | Benchmark 7.101 | Abs -0.003 | Rel -0.040
S0 36 | VOL 0.40 | T 2 | Value 8.502 | Benchmark 8.508 | Abs -0.006 | Rel -0.069
S0 38 | VOL 0.20 | T 1 | Value 3.240 | Benchmark 3.250 | Abs -0.010 | Rel -0.306
S0 38 | VOL 0.20 | T 2 | Value 3.735 | Benchmark 3.745 | Abs -0.010 | Rel -0.273
S0 38 | VOL 0.40 | T 1 | Value 6.166 | Benchmark 6.148 | Abs 0.018 | Rel 0.296
S0 38 | VOL 0.40 | T 2 | Value 7.665 | Benchmark 7.670 | Abs -0.005 | Rel -0.068
S0 40 | VOL 0.20 | T 1 | Value 2.316 | Benchmark 2.314 | Abs 0.002 | Rel 0.086
S0 40 | VOL 0.20 | T 2 | Value 2.878 | Benchmark 2.885 | Abs -0.007 | Rel -0.239
S0 40 | VOL 0.40 | T 1 | Value 5.318 | Benchmark 5.312 | Abs 0.006 | Rel 0.108
S0 40 | VOL 0.40 | T 2 | Value 6.946 | Benchmark 6.920 | Abs 0.026 | Rel 0.381
S0 42 | VOL 0.20 | T 1 | Value 1.626 | Benchmark 1.617 | Abs 0.009 | Rel 0.554
S0 42 | VOL 0.20 | T 2 | Value 2.207 | Benchmark 2.212 | Abs -0.005 | Rel -0.208
S0 42 | VOL 0.40 | T 1 | Value 4.570 | Benchmark 4.582 | Abs -0.012 | Rel -0.263
S0 42 | VOL 0.40 | T 2 | Value 6.264 | Benchmark 6.248 | Abs 0.016 | Rel 0.258
S0 44 | VOL 0.20 | T 1 | Value 1.104 | Benchmark 1.110 | Abs -0.006 | Rel -0.511
S0 44 | VOL 0.20 | T 2 | Value 1.689 | Benchmark 1.690 | Abs -0.001 | Rel -0.084
S0 44 | VOL 0.40 | T 1 | Value 3.946 | Benchmark 3.948 | Abs -0.002 | Rel -0.050
S0 44 | VOL 0.40 | T 2 | Value 5.647 | Benchmark 5.647 | Abs -0.000 | Rel -0.000
Time in Seconds 24.598
Time per Option 1.282
Abs stands for the absolute error in currency units and Rel for the relative error in percent. Although the script simulates 50,000 paths with 50 time intervals each, the execution time is quite fast (a bit more than a second per option on standard notebook under Linux).
Here, the full Python script:
#
# Least-Squares Monte Carlo for American Put Options
# with Antithetic Paths and Moment Matching
# Using In-The-Money Paths Only for Regression
#
# (c) Visixion GmbH
# Script For Illustration Purposes Only.
# January 2013
#
from pylab import *
from time import time
import warnings
warnings.simplefilter('ignore', np.RankWarning)
t0 = time()
seed(100000)
#
# Parameters
#
# Option Parameters
s0 = [36., 38., 40., 42., 44.] # Initial Index Levels
voL = [0.2, 0.4] # Constant Volatilities
tL = [1.0, 2.0] # Times-to-maturity
K = 40. # Strike Price
r = 0.06 # Risk-Free Short Rate
# Simulation Parameters
M0 = 50 # Time Steps p.a.
I = 50000 # Simulation Paths
# Variance Reduction Techniques
antiPaths = True # Antithetic Paths
moMatch = True # Moment Matching
# Benchmark Values
bL = (4.478, 4.840, 7.101, 8.508,
3.250, 3.745, 6.148, 7.670,
2.314, 2.885, 5.312, 6.920,
1.617, 2.212, 4.582, 6.248,
1.110, 1.690, 3.948, 5.647)
#
# Random Numbers
#
def RNG(M, I):
''' Function to generate pseudo-random numbers with variance reduction.
M: number of discrete time intervals
I: number of simulated paths '''
if antiPaths:
ran = standard_normal((M, I / 2))
rand = concatenate((ran, -ran), 1) # Antithetic Variates
else:
rand = standard_normal((M, I))
if moMatch:
rand = rand - mean(rand) # Matching of 1st Moment
rand = rand / std(rand) # Matching of 2nd Moment
return rand
#
# LSM Algorithm
#
z = 0
for S0 in s0:
for vol in voL:
for T in tL:
# Index Level Paths
M = int(T * M0) # Total Time Steps
dt = T / M # Length of Time Interval
df = exp(-r * dt) # Discount Factor
rand = RNG(M, I) # Random Numbers
S = ones((M + 1, I), 'd')
S[1:, :] = exp(cumsum((r - vol ** 2 / 2) * dt +
vol * rand * sqrt(dt),
axis=0))
V = maximum(K - S * S0, 0) # Value Array
# Dynamic Optimization/Exercise Decisions
for t in range(M - 1, 0, -1):
ITM = V[t, :] > 0 # ITM or not
relS = compress(ITM, S[t, :]) # Use only ITM-Paths
relV = compress(ITM, V[t + 1, :] * df)
regval = polyval(polyfit(relS, relV, 3), relS) # Regression
C = zeros(I, 'd')
# Put Regression Values in New Array
put(C, nonzero(ITM), regval)
# Exercise Decision
V[t, :] = where(V[t, :] > C, V[t, :], V[t + 1, :] * df)
V0 = sum(V[1, :] * df) / I # LSM Estimator
# Output
print ("S0 %3.0f | VOL %2.2f | T %2.0f | Value %6.3f |"
+ " Benchmark %6.3f | Abs %6.3f | Rel %6.3f") \
% (S0, vol, T, V0, bL[z], V0 - bL[z],
(V0 - bL[z]) / bL[z] * 100) # Output Values
z += 1
print "Time in Seconds %8.3f" % (time() - t0)
print "Time per Option %8.3f" % ((time() - t0) / len(bL))
Go to Python for Finance (I)
Go to Python for Finance (II)
