simple control variates using terminal spot

sde_mc/helpers.py

@@ -131,3 +131,8 @@ def inner(x):
  i3 = quad(lambda x: c_plus * g(x) * inner(x), 0, 1)[0]
  i4 = quad(lambda x: c_plus * q(x) * tail(x), 1, np.inf)[0]
  return i1 + i2 + i3 + i4
+
+
+def sample_cov(x, y):
+ return torch.sum((x - x.mean()) * (y - y.mean())) / (len(x) - 1)
+

sde_mc/mc.py

@@ -5,7 +5,7 @@
 from .varred import train_diffusion_control_variate, apply_diffusion_control_variate, train_adapted_control_variates, \
  apply_adapted_control_variates, EarlyStopping
 from .nets import NormalJumpsPathData, NormalPathData, AdaptedPathData, Mlp
-from .helpers import partition, mc_estimates, ceil_mult
+from .helpers import partition, mc_estimates, ceil_mult, sample_cov
 from .options import ConstantShortRate
 import gc
 
@@ -326,6 +326,55 @@ def mc_adaptive_cv(models, opt, solver, trials, steps, payoff, discounter, sim_b
  return MCStatistics(mn, sd, train_time + test_time, test_trials)
 
 
+def mc_terminal_cv(num_trials, sde_solver, payoff, discounter=None, bs=None, return_normals=False):
+ """Uses the terminal spot price as a control variate"""
+ if discounter is None:
+ discounter = ConstantShortRate(r=0.0)
+
+ if not bs:
+ start = time.time()
+ out, normals = sde_solver.solve(bs=num_trials, return_normals=return_normals)
+ spots = out[:, -1]
+ payoffs = payoff(spots) * discounter(sde_solver.time_interval)
+ cv = (discounter(sde_solver.time_interval) * spots[:, 0] - sde_solver.sde.init_value[0])
+ b = sample_cov(cv, payoffs) / cv.var()
+ cv_payoffs = payoffs - b * cv
+ mn = cv_payoffs.mean()
+ sd = cv_payoffs.std() / np.sqrt(num_trials)
+ end = time.time()
+ tt = end - start
+
+ return MCStatistics(mn, sd, tt, num_trials, out, payoffs, normals)
+ else:
+ remaining_trials = num_trials
+ sample_sum, sample_sum_sq = 0.0, 0.0
+ start = time.time()
+ first_batch = True
+ while remaining_trials:
+ if remaining_trials < bs:
+ bs = remaining_trials
+
+ remaining_trials -= bs
+ out, normals = sde_solver.solve(bs=bs, return_normals=False)
+ spots = out[:, -1]
+ payoffs = payoff(spots) * discounter(sde_solver.time_interval)
+ cv = (discounter(sde_solver.time_interval) * spots[:, 0] - sde_solver.sde.init_value[0])
+ if first_batch:
+ # estimate b on the first batch
+ b = sample_cov(cv, payoffs) / cv.var()
+ first_batch = False
+ cv_payoffs = payoffs - b * cv
+
+ sample_sum += cv_payoffs.sum()
+ sample_sum_sq += (cv_payoffs**2).sum()
+
+ mn = sample_sum / num_trials
+ sd = torch.sqrt((sample_sum_sq/num_trials - mn**2) * (num_trials / (num_trials-1))) / np.sqrt(num_trials)
+ end = time.time()
+ tt = end-start
+ return MCStatistics(mn, sd, tt, num_trials)
+
+
 def simulate_data(trials, solver, payoff, discounter, bs=1000, inference=False):
  """Simulates trajectories of an SDE and returns the trajectories, payoffs and random variables in a DataLoader
  which can be used for training or inference"""
@@ -378,6 +427,13 @@ def find_num_trials(problem, eps, models=None, init_trials=1e5, bs=1e5):
  return int(trials)
 
 
+def find_num_trials_terminal_cv(problem, eps, init_trials, bs):
+ mc_stats = mc_terminal_cv(init_trials, problem.solver, problem.payoff, problem.discounter, bs)
+ ratio = (mc_stats.sample_std * 1.96 / eps) ** 2
+ trials = np.ceil(ratio * init_trials)
+ return int(trials)
+
+
 def run_mc(problem, eps, bs=1e5, init_trials=1e5):
  trials = find_num_trials(problem, eps, None, init_trials, bs)
  payoff_time = 'adapted' if problem.solver.has_jumps else 'terminal'
@@ -409,3 +465,8 @@ def run_cv_mc(problem, models, opt, eps, train_size, step_factor=30, sim_bs=1e5,
  test_time = mc_stats.time_elapsed
  mc_stats.time_elapsed += train_time_end - train_time_start
  return mc_stats, train_time_end - train_time_start, test_time
+
+
+def run_mc_terminal_cv(problem, eps, bs=1e5, init_trials=1e5):
+ trials = find_num_trials_terminal_cv(problem, eps, init_trials, bs)
+ return mc_terminal_cv(trials, problem.solver, problem.payoff, problem.discounter, bs=bs)

tests/test_mc.py

@@ -51,4 +51,13 @@ def test_simulate_adapted(merton_1d_solver, euro_call, constant_short_rate):
 def test_mc_adaptive_cv(merton_1d_solver, mlps_1d, euro_call, constant_short_rate):
  adam = torch.optim.Adam(list(mlps_1d[0].parameters()) + list(mlps_1d[1].parameters()))
  mc_stats = mc_adaptive_cv(mlps_1d, adam, merton_1d_solver, (100, 100), (10, 20), euro_call,
- constant_short_rate, sim_bs=(100, 100), bs=(10, 10), print_losses=False)
+ constant_short_rate, sim_bs=(100, 100), bs=(10, 10), print_losses=False)
+
+
+def test_mc_terminal_cv(gbm_1d_solver, euro_call, constant_short_rate):
+ mc_stats = mc_terminal_cv(100, gbm_1d_solver, euro_call, constant_short_rate, 10)
+
+
+def test_run_mc_terminal_cv(heston_problem):
+ mc_stats = run_mc_terminal_cv(heston_problem, 0.01, bs=1e3, init_trials=1e4)
+

tests/test_simple.py

@@ -114,3 +114,9 @@ def icdf():
 @pytest.fixture
 def exp_example_levy():
  return ExpExampleLevy(1, 1, 0.5, 2, 0.02, 0.2, 0.1, 0.01)
+
+
+@pytest.fixture
+def heston_problem():
+ return HestonEuroCall.default_params(100, 'cpu')
+