README.md
In shill1729/FeynmanKacSolver: Feynman-Kac PDE solver for Ito diffusions and processes
FeynmanKacSolver
The FeynmanKacSolver provides finite-difference solvers for PDE problems and PIDE problems arising from the Feynman-Kac connection of stochastic processes.
Note: 5/3/2020 temporarily in development overhaul, some examples may be broken
Table of contents
Installation
You can install the github version of FeynmanKacSolver with:
devtools::install_github("shill1729/FeynmanKacSolver")
Simulating sample paths
We can simulate sample paths of continuous Ito processes and diffusions and Ito-Levy processes, that is the class containing jump-diffusions and processes with continuous dynamics from a Brownian motion and jump dynamics from a jump process, here specifically using (compoound) Poisson processes. First we demonstrate the continuous only models.
Continuous paths
Here is a list of examples of stochastic processes with continuous paths with respect to time and the state-variable.
1. Brownian motion
A Brownian motion can easily be simulated by just specifiying the length of time, number of points, and its name.
library(FeynmanKacSolver)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "bm"
# Simulate a sample-path; for Brownian motion, no other arguments need to be specified or passed
x <- sample_path(t = tn, n, continuous.model)
# Compare the density of the increment of the simulated path vs the theoretical density
# The increment of Brownian motion is normally distributed with mean zero and variance equal to the interval of time
# that the increment is over
epdf <- density(diff(x$X))
exact <- dnorm(x = epdf$x, mean = 0, sd = sqrt(tn/n))
# Plotting the sample path and the two densities
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
2. Planar Brownian motion
Two independent Brownian motions determining the coordinates of a point in the plane
library(FeynmanKacSolver)
library(plot3D)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "pbm"
# Simulate a sample-path
x <- sample_path(t = tn, n, continuous.model)
# X_t-X_s ~ N(0, t-s) (in probabilist' notation of mean-variance parameterization of Normal distributions)
epdf.x <- density(diff(x$x))
epdf.y <- density(diff(x$y))
msh <- plot3D::mesh(x = epdf.x$x, y = epdf.y$x)
exact.x <- dnorm(x = msh$x, mean = 0, sd = sqrt(tn/n))
exact.y <- dnorm(x = msh$y, mean = 0, sd = sqrt(tn/n))
par(mfrow = c(2, 1))
# Sample path plot
plot(x$x, x$y, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
surf3D(x = msh$x, y = msh$y, z = exact.x*exact.y, main = "bivariate independent normal density")
3. Spatial Brownian motion
The same as the last example but in three dimensions now: a spatial Brownian motion. Recall that Brownian motion is recurrent in dimensions one and two (i.e. visits every point infinitely often) but transient in dimensions three and greater.
library(FeynmanKacSolver)
library(plot3D)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "sbm"
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model)
# Sample path plot
par(mfrow = c(1, 1))
lines3D(x$x, x$y, x$z, main = paste("sample path of", continuous.model))
4. Pure Ito diffusion
An Ito diffusion with constant drift and volatiltiy coefficients
library(FeynmanKacSolver)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "ito"
# Specify parameters
parameters <- list(mu = 0.2, volat = 0.25)
list2env(x = parameters, envir = .GlobalEnv)
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model, parameters = parameters)
# The density of the log-increment is normally distributed with specified mean and volatility above.
epdf <- density(diff(x$X))
exact <- dnorm(x = epdf$x, mean = mu*tn/n, sd = volat*sqrt(tn/n))
# Plotting the sample path and PDFs
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
5. Geometric Brownian motion
A geometric Brownian motion is like a Ito process that is exponentiated and has its drift compensated by a convexity adjustment.
library(FeynmanKacSolver)
tn <- 1
n <- 20000
# Specify continuous dynamics
continuous.model <- "gbm"
# Specify parameters
parameters <- list(mu = 0.15,
volat = 0.1
)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(spot = 25)
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model, parameters = parameters, IC = IC)
# The log-increment of a GBM has a normal distribution with an adjusted mean
epdf <- density(diff(log(x$X)))
exact <- dnorm(x = epdf$x, mean = (mu-0.5*volat^2)*tn/n, sd = volat*sqrt(tn/n))
# Plotting the sample path and PDFs of log-increments
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of log-increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
6. Vasicek short-rate model
The Vasicek short-rate model, AKA, the Ornstein–Uhlenbeck process but with a financial interpretation. I am not too sure about the exact density yet, but it is on e.g. wikipedia.
library(FeynmanKacSolver)
# TODO: verify distribution of short-rate
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "vasicek"
# Specify parameters
parameters <- list(theta = 0.9, mu = 0.05, volat = 0.25)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(x0 = 0.1)
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model, parameters = parameters, IC = IC)
# Increments are normally distributed with specific mean/sd
epdf <- density(diff(x$X))
exact <- dnorm(x = epdf$x, mean = mu*(1-exp(-theta*tn/n)), sd = sqrt((1-exp(-2*theta*tn/n))*(volat^2)/(2*theta)))
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
7. The Cox-Ingersoll-Ross short-rate model
This is identical to the Vasicek model except for a square root of the state-variable in the Brownian driver. The exact PDF is on wikipedia and will be eventually added to this example and the package as the standard "r,d,p" family for distributions in R. Until then, this example just plots a sample-path.
library(FeynmanKacSolver)
# TODO: verify distribution of short-rate
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "cir"
# Specify parameters
parameters <- list(theta = 0.9, mu = 0.05, volat = 0.25)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(x0 = 0.5)
# Simulate a sample-path
x <- sample_path(t = tn, n, continuous.model, jump.model, parameters, IC)
par(mfrow = c(1, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
Paths with jumps
The paths with jumps involve adding a compound Poisson process of IIDRVs that represent the jump-size of each random jump. There is a MGF condition of finiteness that must be imposed for these to have martingale properties when the drifts are compensated. The current models implemented are "kou", "norm", "unif", and "dkou" which stands for the displaced Kou double-exponential distribution.
1. Merton's jump diffusion model
Merton's jump diffusion model is simply a geometric Brownian motion with jumps driven by a compound Poisson process of Gaussian RVs.
library(FeynmanKacSolver)
tn <- 1
n <- 5000
# Specify continuous dynamics
continuous.model <- "gbm"
# Specify jump dynamics
jump.model <- "norm"
# Specify parameters
parameters <- list(mu = 0.15,
volat = 0.5,
lambda = function(x, t) 50,
mean = -0.08,
sd = 0.002
)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(spot = 25)
# Simulate a sample-path and compute log-returns
x <- sample_path(t = tn, n, continuous.model, jump.model, parameters, IC)
# Computing log increments
log_returns <- diff(log(x$X))
# Empirical PDF
epdf <- density(log_returns)
# Exact PDF
eta <- exp(mean+0.5*sd^2)-1
exact_pdf <- dmerton(epdf$x, t = tn/n, drift = mu-lambda(epdf$x, tn/n)*eta-0.5*volat^2, volat = volat, lambda = lambda(epdf$x, tn/n), a = mean, b = sd)
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model, "-", jump.model))
# PDF of increments plot
plot(epdf$x, exact_pdf, main = "PDF of log-increment", type = "l", ylim = c(0, max(exact_pdf, epdf$y, estimated_pdf)))
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend("topleft", legend = c("exact", "empirical"), col = c("black", "blue"), lty = 1, cex = 0.8)
As an extra bonus, we have provided a basic (and more for informational purposes rather than utility) scheme to estimate the parameters for the Merton model via MLE. Unfortunately the PDF of the Merton model does not necessarily possess a unique parameter set to maximize the likelihood for a given set of observed variates. In other words, the optimization region is possibly flat, and a handful of choices of parameters can be used to fit the shape of the density of empirical data. This toy example demonstrates that by specifiyng a parameter set, fitting a set to it using the MLE routine, and comparing them. While the volatility, mean-rate of jumps, mean jump-size and standard deviation jump-size often gets estimated correctly, within up to acceptable tolerances, the mean-drift of the diffusion component is almost never estimated correctly and highly dependent on initial guesses.
library(FeynmanKacSolver)
thresh_hold <- 0.08
tn <- 1
n <- 1000
# Specify continuous dynamics
continuous.model <- "gbm"
# Specify jump dynamics
jump.model <- "norm"
# Specify parameters
parameters <- list(mu = 0.15,
volat = 0.5,
lambda = function(x, t) 10,
mean = -0.08,
sd = 0.002
)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(spot = 25)
# Simulate a sample-path and compute log-returns
x <- sample_path(t = tn, n, continuous.model, jump.model, parameters, IC)
log_returns <- diff(log(x$X))
# Empirical PDF
epdf <- density(log_returns)
# Exact PDF
eta <- exp(mean+0.5*sd^2)-1
exact_pdf <- dmerton(epdf$x, t = tn/n, drift = mu-lambda(1, 0)*eta-0.5*volat^2, volat = volat, lambda = lambda(1, 0), a = mean, b = sd)
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model, "-", jump.model))
# MLE estimate for parameters from simulation
mle <- merton_mle(log_returns, thresh_hold = thresh_hold, tn = tn, time_step = tn/n, compensate = TRUE, direction = "both")
v <- mle$par
# Model-fit PDF with estimated parameters
eta_est <- exp(v[4]+0.5*v[5]^2)-1
estimated_pdf <- dmerton(epdf$x, t = tn/n, drift = v[1]-v[3]*eta_est-0.5*v[2]^2, volat = v[2], lambda = v[3], a = v[4], b = v[5])
# Finish plotting PDFs
# PDF of increments plot
plot(epdf$x, exact_pdf, main = "PDF of log-increment", type = "l", ylim = c(0, max(exact_pdf, epdf$y, estimated_pdf)))
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
lines(epdf$x, estimated_pdf, col = "green", lty = "dashed")
legend("topleft", legend = c("exact", "empirical", "estimated"), col = c("black", "blue", "green"), lty = 1, cex = 0.8)
exact_par <- unlist(c(parameters[1:2], lambda = parameters$lambda(0, 1), parameters[4:5]))
# Print estimated parameters vs exact
print(rbind(v, exact_par))
Computing CDFs of processes
We can numerically compute the CDF of a process at a future time given the process' value at the current time by solving the PDE with appropriate payoff function. The following examples demonstrate the coinciding of the exact CDFs with the approximate CDFs for three basic examples, Brownian motion, Ito diffusion, and geometric Brownian motion:
1. Numerical CDF vs exact CDF of Brownian motion at future time.
Work in progress. Eventually this will look closer to the syntax for supplying models to the sample_path. Right now, it's pretty ugly, but it works.
library(FeynmanKacSolver)
# Time span to simulate over
maturity <- 1
B <- 3
contract_specs <- list(spot = 1, div = 0, payoff = "cdf", maturity = maturity, payoff_param = list(K = 0))
list2env(x = contract_specs, envir = .GlobalEnv)
numeric_specs <- list(B = B, N = 100, M = 100)
# Domain of CDF
K <- seq(-B, B, length.out = numeric_specs$M)
# # Diffusion coeffcients:
f <- list(f.mu = function(x, t) {0},
f.vo = function(x, t) {1},
f.rate = function(x, t) {0}
)
# Set up space-time grid and discretized coefficient functions
pde.setup <- pde_setup(f = f, contract_specs, numeric_specs)
# Numerical CDF computed via PDE corresponding to conditional expectation of indicator function of Borel sets (-\infty, K]
ncdf <- feynman_kac(K, payoff = payoff, pde.setup)
plot(K, pnorm(K, mean = 0, sd = sqrt(maturity)), type = "l", main = "CDF of Brownian motion at time T", xlab = "x", ylab = "F(x)")
lines(ncdf, col = "blue", lty = "dashed")
legend("topleft", legend = c("Exact", "Approx"), col = c("black", "blue"), lty = 1)
2. Numerical CDF of GBM-Displaced Kou dynamics
Work in progress in terms of the function-argument syntax etc.
library(FeynmanKacSolver)
# Takes about a minute to run...
contract_specs <- list(spot = 100,
div = 0,
maturity = 30/252,
payoff = "straddle",
variational = TRUE,
payoff_param = list(K = 100)
)
numeric_specs <- list(N = 100,
M = 100,
L = 51,
B = 0.5
)
jump_specs <- list(distr = "dkou",
param = list(p = 0.5,
alpha = 0.08,
beta = 0.08,
ku = 0.05,
kd = -0.05
)
)
continuous_specs <- list(
volat = 0.2,
rate = 0.05,
lambda = 5
)
list2env(continuous_specs, envir = .GlobalEnv)
f <- list(f.mu = function(x, t) rate-0.5*volat^2,
f.vo = function(x, t) volat,
f.rate = function(x, t) rate,
f.lambda = function(x, t) lambda
)
K <- seq(0.01, contract_specs$spot*2, length.out = 100)
pide.setup <- pide_setup(f = f, contract_specs = contract_specs, jump_specs = jump_specs, numeric_specs = numeric_specs)
ncdf <- feynman_kac(K, payoff = "cdf", setup = pide.setup, FALSE, TRUE)
plot(ncdf, type = "l")
PDE/PIDE solution grids
Here we provide an example of a straddle.
library(FeynmanKacSolver)
contract_specs <- list(spot = 100,
div = 0,
maturity = 30/252,
payoff = "straddle",
variational = TRUE,
payoff_param = list(K = 100)
)
numeric_specs <- list(N = 100,
M = 100,
L = 51,
B = 0.5
)
jump_specs <- list(distr = "dkou",
param = list(p = 0.5,
alpha = 0.05,
beta = 0.05,
ku = 0.03,
kd = -0.03
)
)
continuous_specs <- list(
volat = 0.2,
rate = 0.05,
lambda = 5
)
list2env(continuous_specs, envir = .GlobalEnv)
f <- list(f.mu = function(x, t) rate-0.5*volat^2,
f.vo = function(x, t) volat,
f.rate = function(x, t) rate,
f.lambda = function(x, t) lambda
)
pide.setup <- pide_setup(f = f, contract_specs = contract_specs, jump_specs = jump_specs, numeric_specs = numeric_specs)
ugrid <- pide_solve(pide.setup, contract_specs$variational, "grid")
greeks <- pide_solve(pide.setup, contract_specs$variational, "greeks")
plot3D::persp3D(z = ugrid)
print(greeks)
shill1729/FeynmanKacSolver documentation built on May 19, 2020, 8:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
FeynmanKacSolver
The FeynmanKacSolver provides finite-difference solvers for PDE problems and PIDE problems arising from the Feynman-Kac connection of stochastic processes.
Note: 5/3/2020 temporarily in development overhaul, some examples may be broken
Table of contents
Installation
You can install the github version of FeynmanKacSolver with:
devtools::install_github("shill1729/FeynmanKacSolver")
Simulating sample paths
We can simulate sample paths of continuous Ito processes and diffusions and Ito-Levy processes, that is the class containing jump-diffusions and processes with continuous dynamics from a Brownian motion and jump dynamics from a jump process, here specifically using (compoound) Poisson processes. First we demonstrate the continuous only models.
Continuous paths
Here is a list of examples of stochastic processes with continuous paths with respect to time and the state-variable.
1. Brownian motion
A Brownian motion can easily be simulated by just specifiying the length of time, number of points, and its name.
library(FeynmanKacSolver)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "bm"
# Simulate a sample-path; for Brownian motion, no other arguments need to be specified or passed
x <- sample_path(t = tn, n, continuous.model)
# Compare the density of the increment of the simulated path vs the theoretical density
# The increment of Brownian motion is normally distributed with mean zero and variance equal to the interval of time
# that the increment is over
epdf <- density(diff(x$X))
exact <- dnorm(x = epdf$x, mean = 0, sd = sqrt(tn/n))
# Plotting the sample path and the two densities
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
2. Planar Brownian motion
Two independent Brownian motions determining the coordinates of a point in the plane
library(FeynmanKacSolver)
library(plot3D)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "pbm"
# Simulate a sample-path
x <- sample_path(t = tn, n, continuous.model)
# X_t-X_s ~ N(0, t-s) (in probabilist' notation of mean-variance parameterization of Normal distributions)
epdf.x <- density(diff(x$x))
epdf.y <- density(diff(x$y))
msh <- plot3D::mesh(x = epdf.x$x, y = epdf.y$x)
exact.x <- dnorm(x = msh$x, mean = 0, sd = sqrt(tn/n))
exact.y <- dnorm(x = msh$y, mean = 0, sd = sqrt(tn/n))
par(mfrow = c(2, 1))
# Sample path plot
plot(x$x, x$y, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
surf3D(x = msh$x, y = msh$y, z = exact.x*exact.y, main = "bivariate independent normal density")
3. Spatial Brownian motion
The same as the last example but in three dimensions now: a spatial Brownian motion. Recall that Brownian motion is recurrent in dimensions one and two (i.e. visits every point infinitely often) but transient in dimensions three and greater.
library(FeynmanKacSolver)
library(plot3D)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "sbm"
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model)
# Sample path plot
par(mfrow = c(1, 1))
lines3D(x$x, x$y, x$z, main = paste("sample path of", continuous.model))
4. Pure Ito diffusion
An Ito diffusion with constant drift and volatiltiy coefficients
library(FeynmanKacSolver)
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "ito"
# Specify parameters
parameters <- list(mu = 0.2, volat = 0.25)
list2env(x = parameters, envir = .GlobalEnv)
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model, parameters = parameters)
# The density of the log-increment is normally distributed with specified mean and volatility above.
epdf <- density(diff(x$X))
exact <- dnorm(x = epdf$x, mean = mu*tn/n, sd = volat*sqrt(tn/n))
# Plotting the sample path and PDFs
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
5. Geometric Brownian motion
A geometric Brownian motion is like a Ito process that is exponentiated and has its drift compensated by a convexity adjustment.
library(FeynmanKacSolver)
tn <- 1
n <- 20000
# Specify continuous dynamics
continuous.model <- "gbm"
# Specify parameters
parameters <- list(mu = 0.15,
volat = 0.1
)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(spot = 25)
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model, parameters = parameters, IC = IC)
# The log-increment of a GBM has a normal distribution with an adjusted mean
epdf <- density(diff(log(x$X)))
exact <- dnorm(x = epdf$x, mean = (mu-0.5*volat^2)*tn/n, sd = volat*sqrt(tn/n))
# Plotting the sample path and PDFs of log-increments
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of log-increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
6. Vasicek short-rate model
The Vasicek short-rate model, AKA, the Ornstein–Uhlenbeck process but with a financial interpretation. I am not too sure about the exact density yet, but it is on e.g. wikipedia.
library(FeynmanKacSolver)
# TODO: verify distribution of short-rate
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "vasicek"
# Specify parameters
parameters <- list(theta = 0.9, mu = 0.05, volat = 0.25)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(x0 = 0.1)
# Simulate a sample-path
x <- sample_path(t = tn, n = n, continuous.model = continuous.model, parameters = parameters, IC = IC)
# Increments are normally distributed with specific mean/sd
epdf <- density(diff(x$X))
exact <- dnorm(x = epdf$x, mean = mu*(1-exp(-theta*tn/n)), sd = sqrt((1-exp(-2*theta*tn/n))*(volat^2)/(2*theta)))
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
# PDF of increments plot
plot(epdf$x, exact, main = "PDF of increment", type = "l")
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend(x = "topleft", legend = c("exact", "approx"), col = c("black", "blue"), lty = 1, cex = 0.8)
7. The Cox-Ingersoll-Ross short-rate model
This is identical to the Vasicek model except for a square root of the state-variable in the Brownian driver. The exact PDF is on wikipedia and will be eventually added to this example and the package as the standard "r,d,p" family for distributions in R. Until then, this example just plots a sample-path.
library(FeynmanKacSolver)
# TODO: verify distribution of short-rate
tn <- 1
n <- 10000
# Specify continuous dynamics
continuous.model <- "cir"
# Specify parameters
parameters <- list(theta = 0.9, mu = 0.05, volat = 0.25)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(x0 = 0.5)
# Simulate a sample-path
x <- sample_path(t = tn, n, continuous.model, jump.model, parameters, IC)
par(mfrow = c(1, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model))
Paths with jumps
The paths with jumps involve adding a compound Poisson process of IIDRVs that represent the jump-size of each random jump. There is a MGF condition of finiteness that must be imposed for these to have martingale properties when the drifts are compensated. The current models implemented are "kou", "norm", "unif", and "dkou" which stands for the displaced Kou double-exponential distribution.
1. Merton's jump diffusion model
Merton's jump diffusion model is simply a geometric Brownian motion with jumps driven by a compound Poisson process of Gaussian RVs.
library(FeynmanKacSolver)
tn <- 1
n <- 5000
# Specify continuous dynamics
continuous.model <- "gbm"
# Specify jump dynamics
jump.model <- "norm"
# Specify parameters
parameters <- list(mu = 0.15,
volat = 0.5,
lambda = function(x, t) 50,
mean = -0.08,
sd = 0.002
)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(spot = 25)
# Simulate a sample-path and compute log-returns
x <- sample_path(t = tn, n, continuous.model, jump.model, parameters, IC)
# Computing log increments
log_returns <- diff(log(x$X))
# Empirical PDF
epdf <- density(log_returns)
# Exact PDF
eta <- exp(mean+0.5*sd^2)-1
exact_pdf <- dmerton(epdf$x, t = tn/n, drift = mu-lambda(epdf$x, tn/n)*eta-0.5*volat^2, volat = volat, lambda = lambda(epdf$x, tn/n), a = mean, b = sd)
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model, "-", jump.model))
# PDF of increments plot
plot(epdf$x, exact_pdf, main = "PDF of log-increment", type = "l", ylim = c(0, max(exact_pdf, epdf$y, estimated_pdf)))
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
legend("topleft", legend = c("exact", "empirical"), col = c("black", "blue"), lty = 1, cex = 0.8)
As an extra bonus, we have provided a basic (and more for informational purposes rather than utility) scheme to estimate the parameters for the Merton model via MLE. Unfortunately the PDF of the Merton model does not necessarily possess a unique parameter set to maximize the likelihood for a given set of observed variates. In other words, the optimization region is possibly flat, and a handful of choices of parameters can be used to fit the shape of the density of empirical data. This toy example demonstrates that by specifiyng a parameter set, fitting a set to it using the MLE routine, and comparing them. While the volatility, mean-rate of jumps, mean jump-size and standard deviation jump-size often gets estimated correctly, within up to acceptable tolerances, the mean-drift of the diffusion component is almost never estimated correctly and highly dependent on initial guesses.
library(FeynmanKacSolver)
thresh_hold <- 0.08
tn <- 1
n <- 1000
# Specify continuous dynamics
continuous.model <- "gbm"
# Specify jump dynamics
jump.model <- "norm"
# Specify parameters
parameters <- list(mu = 0.15,
volat = 0.5,
lambda = function(x, t) 10,
mean = -0.08,
sd = 0.002
)
list2env(x = parameters, envir = .GlobalEnv)
# Specify initial values
IC <- list(spot = 25)
# Simulate a sample-path and compute log-returns
x <- sample_path(t = tn, n, continuous.model, jump.model, parameters, IC)
log_returns <- diff(log(x$X))
# Empirical PDF
epdf <- density(log_returns)
# Exact PDF
eta <- exp(mean+0.5*sd^2)-1
exact_pdf <- dmerton(epdf$x, t = tn/n, drift = mu-lambda(1, 0)*eta-0.5*volat^2, volat = volat, lambda = lambda(1, 0), a = mean, b = sd)
par(mfrow = c(2, 1))
# Sample path plot
plot(x, type = "l", main = paste("sample path of", continuous.model, "-", jump.model))
# MLE estimate for parameters from simulation
mle <- merton_mle(log_returns, thresh_hold = thresh_hold, tn = tn, time_step = tn/n, compensate = TRUE, direction = "both")
v <- mle$par
# Model-fit PDF with estimated parameters
eta_est <- exp(v[4]+0.5*v[5]^2)-1
estimated_pdf <- dmerton(epdf$x, t = tn/n, drift = v[1]-v[3]*eta_est-0.5*v[2]^2, volat = v[2], lambda = v[3], a = v[4], b = v[5])
# Finish plotting PDFs
# PDF of increments plot
plot(epdf$x, exact_pdf, main = "PDF of log-increment", type = "l", ylim = c(0, max(exact_pdf, epdf$y, estimated_pdf)))
lines(epdf$x, epdf$y, col = "blue", lty = "dashed")
lines(epdf$x, estimated_pdf, col = "green", lty = "dashed")
legend("topleft", legend = c("exact", "empirical", "estimated"), col = c("black", "blue", "green"), lty = 1, cex = 0.8)
exact_par <- unlist(c(parameters[1:2], lambda = parameters$lambda(0, 1), parameters[4:5]))
# Print estimated parameters vs exact
print(rbind(v, exact_par))
Computing CDFs of processes
We can numerically compute the CDF of a process at a future time given the process' value at the current time by solving the PDE with appropriate payoff function. The following examples demonstrate the coinciding of the exact CDFs with the approximate CDFs for three basic examples, Brownian motion, Ito diffusion, and geometric Brownian motion:
1. Numerical CDF vs exact CDF of Brownian motion at future time.
Work in progress. Eventually this will look closer to the syntax for supplying models to the sample_path. Right now, it's pretty ugly, but it works.
library(FeynmanKacSolver)
# Time span to simulate over
maturity <- 1
B <- 3
contract_specs <- list(spot = 1, div = 0, payoff = "cdf", maturity = maturity, payoff_param = list(K = 0))
list2env(x = contract_specs, envir = .GlobalEnv)
numeric_specs <- list(B = B, N = 100, M = 100)
# Domain of CDF
K <- seq(-B, B, length.out = numeric_specs$M)
# # Diffusion coeffcients:
f <- list(f.mu = function(x, t) {0},
f.vo = function(x, t) {1},
f.rate = function(x, t) {0}
)
# Set up space-time grid and discretized coefficient functions
pde.setup <- pde_setup(f = f, contract_specs, numeric_specs)
# Numerical CDF computed via PDE corresponding to conditional expectation of indicator function of Borel sets (-\infty, K]
ncdf <- feynman_kac(K, payoff = payoff, pde.setup)
plot(K, pnorm(K, mean = 0, sd = sqrt(maturity)), type = "l", main = "CDF of Brownian motion at time T", xlab = "x", ylab = "F(x)")
lines(ncdf, col = "blue", lty = "dashed")
legend("topleft", legend = c("Exact", "Approx"), col = c("black", "blue"), lty = 1)
2. Numerical CDF of GBM-Displaced Kou dynamics
Work in progress in terms of the function-argument syntax etc.
library(FeynmanKacSolver)
# Takes about a minute to run...
contract_specs <- list(spot = 100,
div = 0,
maturity = 30/252,
payoff = "straddle",
variational = TRUE,
payoff_param = list(K = 100)
)
numeric_specs <- list(N = 100,
M = 100,
L = 51,
B = 0.5
)
jump_specs <- list(distr = "dkou",
param = list(p = 0.5,
alpha = 0.08,
beta = 0.08,
ku = 0.05,
kd = -0.05
)
)
continuous_specs <- list(
volat = 0.2,
rate = 0.05,
lambda = 5
)
list2env(continuous_specs, envir = .GlobalEnv)
f <- list(f.mu = function(x, t) rate-0.5*volat^2,
f.vo = function(x, t) volat,
f.rate = function(x, t) rate,
f.lambda = function(x, t) lambda
)
K <- seq(0.01, contract_specs$spot*2, length.out = 100)
pide.setup <- pide_setup(f = f, contract_specs = contract_specs, jump_specs = jump_specs, numeric_specs = numeric_specs)
ncdf <- feynman_kac(K, payoff = "cdf", setup = pide.setup, FALSE, TRUE)
plot(ncdf, type = "l")
PDE/PIDE solution grids
Here we provide an example of a straddle.
library(FeynmanKacSolver)
contract_specs <- list(spot = 100,
div = 0,
maturity = 30/252,
payoff = "straddle",
variational = TRUE,
payoff_param = list(K = 100)
)
numeric_specs <- list(N = 100,
M = 100,
L = 51,
B = 0.5
)
jump_specs <- list(distr = "dkou",
param = list(p = 0.5,
alpha = 0.05,
beta = 0.05,
ku = 0.03,
kd = -0.03
)
)
continuous_specs <- list(
volat = 0.2,
rate = 0.05,
lambda = 5
)
list2env(continuous_specs, envir = .GlobalEnv)
f <- list(f.mu = function(x, t) rate-0.5*volat^2,
f.vo = function(x, t) volat,
f.rate = function(x, t) rate,
f.lambda = function(x, t) lambda
)
pide.setup <- pide_setup(f = f, contract_specs = contract_specs, jump_specs = jump_specs, numeric_specs = numeric_specs)
ugrid <- pide_solve(pide.setup, contract_specs$variational, "grid")
greeks <- pide_solve(pide.setup, contract_specs$variational, "greeks")
plot3D::persp3D(z = ugrid)
print(greeks)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.