This runs the following simulations.
Each process is simulated first using the Euler scheme and then repeated using the Exact method.
In [1]:
source('one_dimension/general_euler.r')
source('one_dimension/driftless_exact.r')
source('one_dimension/drifted_exact.r')
source('multiple_dimensions/general_euler_multi_d.r')
source('multiple_dimensions/drifted_exact_multi_d.r')
source('run_simulation.r')
source('coefficients.r')
source('payoffs.r')
source('utility_functions.r')
In [3]:
num_paths <- 100000
beta <- 0.2
Texp <- 1
X0 <- 1
strikes <- seq(0.6, 1.5, by=0.1)
one_d_from_paper_vols <- data.frame(K=strikes)
In [4]:
title <- 'Euler - 10 timesteps'
num_steps_euler <- 10
euler51_process <- General_Euler(num_paths, num_steps_euler,
mu=make_constant_coefficient(0),
sigma=sigma_from_paper,
Texp=Texp, X0=X0)
one_d_from_paper_vols[title] <- run_simulation(euler51_process, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE)
In [5]:
title <- 'Euler - 50 timesteps'
num_steps_euler <- 50
euler51_process <- General_Euler(num_paths, num_steps_euler,
mu=make_constant_coefficient(0),
sigma=sigma_from_paper,
Texp=Texp, X0=X0)
one_d_from_paper_vols[title] <- run_simulation(euler51_process, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE)
In [6]:
title <- 'Exact'
exact51_process <- Driftless_Exact(num_paths, beta,
sigma_from_paper,
sigma_deriv_from_paper,
Texp=Texp, X0=X0)
one_d_from_paper_vols[title] <- run_simulation(exact51_process, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE)
In [7]:
title <- 'Euler'
num_steps_euler <- 10
euler52_process <- General_Euler(num_paths, num_steps_euler,
mu=mu_1d_from_paper,
sigma=make_constant_coefficient(1),
Texp=Texp, X0=X0,
convert_y_to_x=convert_y_to_x_1d,
convert_x_to_y=convert_x_to_y_1d)
one_d_from_paper_vols[title] <- run_simulation(euler52_process, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE)
In [8]:
title <- 'Exact - Lamperti'
exact52_process <- Drifted_Exact(num_paths, beta,
mu_1d_from_paper,
sigma0=1,
Texp=Texp, X0=X0,
convert_y_to_x=convert_y_to_x_1d,
convert_x_to_y=convert_x_to_y_1d)
one_d_from_paper_vols[title] <- run_simulation(exact52_process, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE)
In [9]:
plot_values(one_d_from_paper_vols, 'One-Dimensional Process from Reference Paper')
Out[9]:
In [5]:
d <- 4
num_paths <- 10000*d
num_steps_euler <- 100
beta <- 0.2
Texp <- 1
X0 <- rep(1, d)
strikes <- seq(0.6, 1.5, by=0.1)
multi_d_from_paper_vols <- data.frame(K=strikes)
In [6]:
title <- 'Euler Section 5.2'
euler_section_52_process <- General_Euler_Multi_D(num_paths, num_steps_euler,
mu=sapply(1:d, make_mu_section_52_euler),
sigma=sapply(1:d, make_sigma_section_52_euler),
cov_matrix=matrix(.5, nrow=d, ncol=d) + diag(d)*.5,
Texp=Texp, X0=X0)
multi_d_from_paper_vols[title] <- run_simulation(euler_section_52_process, strikes, Texp, X0[1], make_basket_call_payoff, plot_title=title)
In [7]:
title <- 'Exact Section 5.2'
exact_section_52_process <- Drifted_Exact_Multi_D(num_paths, beta,
mu=sapply(1:d, make_mu_from_paper_multi_dimension_lamperti),
sigma0=diag(d),
cov_matrix=matrix(.5, nrow=d, ncol=d) + diag(d)*.5,
Texp=Texp, X0=X0,
convert_y_to_x=convert_y_to_x_multi_d,
convert_x_to_y=convert_x_to_y_multi_d)
multi_d_from_paper_vols[title] <- run_simulation(exact_section_52_process, strikes, Texp, X0[1], make_basket_call_payoff_exact, plot_title=title)
In [8]:
plot_values(multi_d_from_paper_vols, 'PHL Multi-Dimensional')
Out[8]:
In [10]:
paramsBCC <- list(lambda = 1.15,rho = -0.64,eta = 0.39, vbar = 0.04,v0 = 0.04)
paramsBCC1 <- list(lambda = 1.15,rho = -0.64,eta = 0.39/2,vbar = 0.04,v0 = 0.04)
paramsBCC2 <- list(lambda = 1.15,rho = -0.64,eta = 0.39*2,vbar = 0.04,v0 = 0.04)
c(feller(paramsBCC), feller(paramsBCC1), feller(paramsBCC2))
Out[10]:
In [11]:
params <- paramsBCC1
cat('Is the Feller condition satisfied: ', feller(params))
y_to_x_hv <- make_convert_y_to_x_heston_variance(params$eta, params$v0)
x_to_y_hv <- make_convert_x_to_y_heston_variance(params$eta, params$v0)
num_paths <- 1000000
num_steps_euler <- 100
beta <- 0.2
Texp <- 1
X0 <- params$v0
strikes <- seq(0.05, 0.12, by=0.01)
heston_variance_vols <- data.frame(K=strikes)
In [12]:
title <- 'Euler'
euler_heston_variance <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_heston_variance(params$lambda, params$vbar),
sigma=make_sigma_cir(params$eta),
Texp=Texp, X0=X0)
heston_variance_vols[title] <- run_simulation(euler_heston_variance, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
The BCC and BCC2 parameters fail the Feller condition. Using those parameters would make 0 an attainable variance. We therefore use the BCC1 parameters.
Using parameters from regular Heston variance
In [13]:
title <- 'Euler - Lamperti'
euler_heston_variance_lamperti <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_heston_variance_lamperti(params$lambda, params$vbar,
params$eta, params$v0),
sigma=make_constant_coefficient(1),
Texp=Texp, X0=X0,
convert_y_to_x=y_to_x_hv,
convert_x_to_y=x_to_y_hv)
heston_variance_vols[title] <- run_simulation(euler_heston_variance_lamperti, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [14]:
title <- 'Exact - Lamperti'
exact_heston_variance_lamperti <- Drifted_Exact(num_paths, beta,
mu=make_mu_heston_variance_lamperti(params$lambda, params$vbar,
params$eta, params$v0),
sigma0=1,
Texp=Texp, X0=X0,
convert_y_to_x=y_to_x_hv,
convert_x_to_y=x_to_y_hv)
heston_variance_vols[title] <- run_simulation(exact_heston_variance_lamperti, strikes, Texp, X0,
make_call_payoff, plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [15]:
plot_values(heston_variance_vols, 'Heston Variance Process', ylab='Prices')
Out[15]:
In [16]:
params <- list(lambda = 1.15,rho = -0.64,eta = 0.39/2,vbar = 1.0,v0 = 1.0)
cat('Is the Feller condition satisfied: ', feller(params))
y_to_x_hv2 <- make_convert_y_to_x_heston_variance(params$eta, params$v0)
x_to_y_hv2 <- make_convert_x_to_y_heston_variance(params$eta, params$v0)
num_paths <- 1000000
num_steps_euler <- 100
beta <- 0.2
Texp <- 1
X0 <- params$v0
strikes <- seq(0.5, 1.5, by=0.1)
heston_variance_vols_2 <- data.frame(K=strikes)
In [17]:
title <- 'Euler'
euler_heston_variance_2 <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_heston_variance(params$lambda, params$vbar),
sigma=make_sigma_cir(params$eta),
Texp=Texp, X0=X0)
heston_variance_vols_2[title] <- run_simulation(euler_heston_variance_2, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [18]:
title <- 'Euler - Lamperti'
euler_heston_variance_lamperti_2 <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_heston_variance_lamperti(params$lambda, params$vbar,
params$eta, params$v0),
sigma=make_constant_coefficient(1),
Texp=Texp, X0=X0,
convert_y_to_x=y_to_x_hv2,
convert_x_to_y=x_to_y_hv2)
heston_variance_vols_2[title] <- run_simulation(euler_heston_variance_lamperti_2, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [19]:
title <- 'Exact - Lamperti'
exact_heston_variance_lamperti_2 <- Drifted_Exact(num_paths, beta,
mu=make_mu_heston_variance_lamperti(params$lambda, params$vbar,
params$eta, params$v0),
sigma0=1,
Texp=Texp, X0=X0,
convert_y_to_x=y_to_x_hv2,
convert_x_to_y=x_to_y_hv2)
heston_variance_vols_2[title] <- run_simulation(exact_heston_variance_lamperti_2, strikes, Texp, X0,
make_call_payoff, plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [20]:
plot_values(heston_variance_vols_2, 'Heston Variance Process - Avoiding Truncation', ylab='Prices')
Out[20]:
In [21]:
num_paths <- 1000000
num_steps_euler <- 100
beta <- 0.2
Texp <- 1
X0 <- 1
mu0 <- 0.1
sigma0 <- 0.3
y_to_x_bs <- make_convert_y_to_x_black_scholes(sigma0, X0)
x_to_y_bs <- make_convert_x_to_y_black_scholes(sigma0, X0)
strikes <- seq(0.5, 1.5, by=0.1)
black_scholes_vols <- data.frame(K=strikes)
In [22]:
title <- 'Euler'
euler_black_scholes <- General_Euler(num_paths, num_steps_euler,
mu=make_affine_coefficient(mu0, 0),
sigma=make_affine_coefficient(sigma0, 0),
Texp=Texp, X0=X0)
black_scholes_vols[title] <- run_simulation(euler_black_scholes, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [23]:
title <- 'Euler - Lamperti'
euler_black_scholes_lamperti <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_black_scholes_lamperti(mu0, sigma0),
sigma=make_constant_coefficient(1),
Texp=Texp, X0=X0,
convert_y_to_x=y_to_x_bs,
convert_x_to_y=x_to_y_bs)
black_scholes_vols[title] <- run_simulation(euler_black_scholes_lamperti, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [24]:
title <- 'Exact - Lamperti'
exact_black_scholes_lamperti <- Drifted_Exact(num_paths, beta,
mu=make_mu_black_scholes_lamperti(mu0, sigma0),
sigma0=1,
Texp=Texp, X0=X0,
convert_y_to_x=y_to_x_bs,
convert_x_to_y=x_to_y_bs)
black_scholes_vols[title] <- run_simulation(exact_black_scholes_lamperti, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [25]:
plot_values(black_scholes_vols, 'Black-Scholes Process', ylab='Prices')
Out[25]:
where $\sigma_0 = 0.3 X_0$
Note that price under this model is $$ \mathbb{E}[(X_T-K)^+|X_0] = (\mu-K) N(a) + \frac{\sigma_0\sqrt{T}}{\sqrt{2\pi}}\exp\left(-\frac{a^2}{2}\right) $$ where \begin{align*} \mu &= X_0-\frac{1}{2\pi n}\cos\left(2\pi nT - \frac{\pi}{2}\right)+T \\ a &= \frac{\mu-K}{\sigma_0\sqrt{T}} \end{align*} and $N(\cdot)$ is the standard normal CDF.
In [26]:
num_paths <- 1000000
beta <- 0.2
Texp <- 1
X0 <- 10
n <- 100
sigma0 <- 0.3*X0
strikes <- seq(9.5, 10.5, by=0.1)
sin_vols <- data.frame(K=strikes)
In [27]:
true_price_sine <- function(K){
mu <- X0 - cos(2*pi*n*Texp - 0.5*pi)/(2*pi*n) + Texp
a <- (mu-K)/(sigma0*sqrt(Texp))
(mu-K)*pnorm(a) + sigma0*sqrt(Texp)/sqrt(2*pi)*exp(-0.5*a^2)
}
sin_vols['True Value'] <- true_price_sine(strikes)
In [28]:
title <- 'Euler - 100 Timesteps'
num_steps_euler <- 100
euler_sin <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_sin(a=1, k=2*pi*n, phi=-0.5*pi, b=1),
sigma=make_constant_coefficient(sigma0),
Texp=Texp, X0=X0)
sin_vols[title] <- run_simulation(euler_sin, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [29]:
title <- 'Euler - 99 Timesteps'
num_steps_euler <- 99
euler_sin <- General_Euler(num_paths, num_steps_euler,
mu=make_mu_sin(a=1, k=2*pi*n, phi=-0.5*pi, b=1),
sigma=make_constant_coefficient(sigma0),
Texp=Texp, X0=X0)
sin_vols[title] <- run_simulation(euler_sin, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [30]:
title <- 'Exact'
exact_sin <- Drifted_Exact(num_paths, beta,
mu=make_mu_sin(a=1, k=2*pi*n, phi=-0.5*pi, b=1),
sigma0=sigma0,
Texp=Texp, X0=X0)
sin_vols[title] <- run_simulation(exact_sin, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [31]:
plot_values(sin_vols, 'Sine Process', ylab='Prices')
Out[31]:
In [32]:
num_paths <- 1000000
num_steps_euler <- 100
beta <- 0.2
Texp <- 1
X0 <- 1
lambda <- 0.1
xbar <- 1
sigma0 <- 0.3
strikes <- seq(0.5, 1.5, by=0.1)
ou_vols <- data.frame(K=strikes)
In [33]:
title <- 'Euler'
euler_ou <- General_Euler(num_paths, num_steps_euler,
mu=make_mean_reverting_coefficient(lambda, xbar),
sigma=make_constant_coefficient(sigma0),
Texp=Texp, X0=X0)
ou_vols[title] <- run_simulation(euler_ou, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [34]:
title <- 'Exact'
euler_ou <- Drifted_Exact(num_paths, beta,
mu=make_mean_reverting_coefficient(lambda, xbar),
sigma=sigma0,
Texp=Texp, X0=X0)
ou_vols[title] <- run_simulation(euler_ou, strikes, Texp, X0, make_call_payoff,
plot_title=title, num_paths_to_plot=FALSE, return_ivols=FALSE)
In [35]:
plot_values(ou_vols, 'Ornstein-Uhlenbeck Process', ylab='Prices')
Out[35]: