# How to compute standard errors of an estimator with antithetic variates?

I'm pricing American options using Longstaff and Schwartz Least square method.

When using the following Python code, I obtain nearly the same prices and standard errors as in the Valuing American Options by Simulation: A Simple Least-Squares Approach by Longstaff and Schwartz (2001):

import numpy as np
import numpy.random as npr
import warnings
warnings.simplefilter('ignore')
from numpy.polynomial.laguerre import lagfit, lagval

# Same parameters as in the original paper
class par: pass
par.S0 = 36
par.K = 40
par.r = 0.06
par.sigma = 0.2
par.T = 1.0
par.I = 100000
par.M = 50

def gen_sn(par,anti_path):
''' Function to generate random numbers for simulation.
Parameters
==========
M : int
number of time intervals for discretization
I : int
number of paths to be simulated
'''
if anti_path is True:
sn = npr.standard_normal((par.M + 1, par.I//2))
else:
sn = npr.standard_normal((par.M + 1, par.I))
return sn

def gbm_mcs_amer(par):
''' Valuation of American option in Black-Scholes-Merton
by Monte Carlo simulation by LS algorithm.
Parameters
==========
S0 : spot price
K : strike float
r : riskless interest rate
I : int, number of paths to be simulated
T : time to maturity, in years
M : int, number of time intervals for discretization
sigma : vol
Returns
=======
C0 : float
estimated present value of American call option
'''
dt = par.T / par.M
df = np.exp(-par.r * dt) # discount function

# Generation of underlying asset process
# Stock Price Paths
S = par.S0 * np.exp(np.cumsum((par.r - 0.5 * par.sigma ** 2) * dt
+ par.sigma * np.sqrt(dt) * sn, axis=0)) # by exponentiating the Brownian motion
S[0] = par.S0 # Initiliazing underlying path

# put option pay-off
h = np.maximum(par.K - S, 0)
# LS algorithm
V = np.copy(h)
for t in range(par.M - 1, 0, -1):
reg = lagfit(S[t], V[t + 1] * df, 10)
C = lagval(S[t], reg)
V[t] = np.where(C > h[t], V[t + 1] * df, h[t])

# MCS estimator
y_i = df * V[1]
C0 = np.mean(y_i)
SE = np.std(y_i, ddof=1) / np.sqrt(par.I) # ddof = 1, bc it's sample std. dev

return C0, SE

# Regular Estimate loop to check the estimates
sn = gen_sn(par,False)

print("Reg","T:",par.T,"sigma:",par.sigma)
for par.S0 in range(36,44+1,2):
print("S0:",par.S0,"Price,SE:",gbm_mcs_amer_reg(par)[0],gbm_mcs_amer_reg(par)[1])


I then attempt implementing antithetic variates price as Glasserman proposes Antithetic paths estimator and the standard error as Boyle and Glasserman proposes Antithetic paths standard error:

def gbm_mcs_amer_AP(par):
dt = par.T / par.M
df = np.exp(-par.r * dt) # discount function

# Generation of underlying asset process
# Stock Price Paths
S = par.S0 * np.exp(np.cumsum((par.r - 0.5 * par.sigma ** 2) * dt
+ par.sigma * np.sqrt(dt) * sn , axis=0)) # by exponentiating the Brownian motion
S[0] = par.S0
S1 = par.S0 * np.exp(np.cumsum((par.r - 0.5 * par.sigma ** 2) * dt
+ par.sigma * np.sqrt(dt) * -sn , axis=0)) # Antithetic paths
S1[0] = par.S0

# put option pay-off
h = np.maximum(par.K - S, 0)
h1 = np.maximum(par.K - S1, 0)
# LS algorithm
V = np.copy(h)
V1 = np.copy(h1)
for t in range(par.M - 1, 0, -1):
reg = lagfit(S[t], V[t + 1] * df, 10)
C = lagval(S[t], reg)
V[t] = np.where(C > h[t], V[t + 1] * df, h[t])
reg1 = lagfit(S1[t], V1[t + 1] * df, 10)
C1 = lagval(S1[t], reg1)
V1[t] = np.where(C1 > h1[t], V1[t + 1] * df, h1[t])

# MCS estimator
y_i = df * (V[1]+V1[1])/2 # avg. pairs
C0 = np.mean(y_i)
SE = np.std(y_i, ddof=1) / np.sqrt(par.I) # Sample std. dev. of avg. pairs

return C0, SE

# AP Estimate loop to check the estimates
sn = gen_sn(par,True)

print("AP","T:",par.T,"sigma:",par.sigma)
for par.S0 in range(36,44+1,2):
print("S0:",par.S0,"Price,SE:",gbm_mcs_amer_AP(par)[0],gbm_mcs_amer_AP(par)[1])


The antithetic variate prices are close to the prices in the Longstaff and Schwartz paper, but the standard errors seem to be much smaller.

New contributor
AverageJoe is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.