[Python] Least-Squares Monte Carlo implementation

I am implementing the least squares monte carlo method to price american put (and maybe call) options in the 4/2 stochastic volatility process. For this I generate paths first and then I want to use the least squares monte carlo method introduced in "Valuing American Options by Simulation: A Simple Least-Squares Approach (2001)". I tested my code with the example given in the paper and got the correct results, but for further testing I looked up some values for the Heston model and got the wrong results for those (by setting a=1, b=0 in the 4/2 model I get the Heston model). I have been looking at my code for days to find the mistake but I just can´t, I would really appreciate some help...

I simulate the prices for lists of strikes K and maturities T. Initially my code was really slow so I tried to speed up the calculation by copying values from previous maturities, if the values are the same at the previous maturity. For this I use boolean indexing (at which I am really bad at and my code is probably a mess), so I can skip rows that can use previously calculated values.

def pricer(self, K, type = "Put"):
if type not in ["Put", "Call"]:
raise ValueError("Invalid option type. Choose either 'Put' or 'Call'.")
if self.S_paths is None: # Check if paths already exist
self.simulate()
df = np.exp(-self.r * self.dt)
s = np.zeros((self.S_paths.shape[0], self.S_paths.shape[1], len(self.T), len(K)))


for k in range(len(K)): # Loop over strikes
payoff_fun = lambda x: np.maximum(K[k] - x, 0) if type == "Put" else np.maximum(x - K[k], 0) # Payoff function

for t in range(len(self.T)): # Loop over maturities
s[:, :self.N[t], t, k] = payoff_fun(self.S_paths[:, :self.N[t]]) # Calculate payoffs

chosen_row_idx = np.full(self.S_paths.shape[0], True, dtype=bool) # All rows are initially "chosen"
for i in range(self.N[t] - 2, 0, -1): # Iterate backwards in time
# Get indices of currently chosen rows
current_rows = np.where(chosen_row_idx)[0]
if current_rows.size == 0:
break

# Identify in-the-money paths among the currently chosen rows
in_the_money_flag = s[current_rows, i, t, k] > 0
if not np.any(in_the_money_flag):
continue

# Redundancy check for earlier maturities
if t > 0 and i in self.N[:t - 1]: # Check if i corresponds to a prior maturity
t_temp = self.N.index(i) # Find the earlier maturity index
temp_done_idx = np.isclose(
s[chosen_row_idx, i, t, k], # Current payoff values
s[chosen_row_idx, i, t_temp, k], # Payoff values for prior maturity
atol=1e-5
)
rows_done_idx = np.full(chosen_row_idx.shape, False, dtype=bool)
rows_done_idx[chosen_row_idx] = temp_done_idx # Update redundant rows
s[rows_done_idx, :i, t, k] = s[rows_done_idx, :i, t_temp, k] # Copy earlier values for redundant rows
chosen_row_idx[rows_done_idx] = False # Exclude redundant rows from further calculations

# Regular Least Squares Monte Carlo calculations
in_money_rows = current_rows[in_the_money_flag]
not_in_money_rows = current_rows[~in_the_money_flag] # For discounting only

# Discount the out-of-money paths:
s[not_in_money_rows, i, t, k] = s[not_in_money_rows, i + 1, t, k] * df

# Compute continuation values via regression for in-money rows:
holding_value = s[in_money_rows, i + 1, t, k] * df
X = self.S_paths[in_money_rows, i]
X2 = X * X
Xs = np.column_stack([X, X2])
Y = holding_value

beta, residuals, rank, sing_vals = np.linalg.lstsq(Xs, Y, rcond=None)
conditional_exp = Xs @ beta
#model_sklearn = LinearRegression()
#model = model_sklearn.fit(Xs, Y) # Fit regression model
#conditional_exp = model.predict(Xs) # Predict continuation value

# Early exercise decision:
immediate_payoff = s[in_money_rows, i, t, k]
early_exercise = immediate_payoff > conditional_exp
s[in_money_rows, i, t, k] = np.where(early_exercise, immediate_payoff, holding_value)

final_prices = np.zeros((len(self.T), len(K))) # Maturities x Strikes

path_values_at_time_1 = s[:, 1, :, :]
discounted_values = np.sum(path_values_at_time_1, axis=0) * df / self.paths

K = np.array(K)
intrinsic_values = np.maximum(K - self.S0, 0) if type == "Put" else np.maximum(self.S0 - K, 0)
final_prices = np.maximum(discounted_values, intrinsic_values)
return final_prices


I would also really appreciate some tips how I could make my code faster, the runtime is really an issue for me. Thanks!

Continue reading...