sde.sim: Simulation of stochastic differential equation

Description

Generic interface to different methods of simulation of solutions to stochastic differential equations.

Usage

sde.sim(t0 = 0, T = 1, X0 = 1, N = 100, delta, drift, sigma, 
   drift.x, sigma.x, drift.xx, sigma.xx, drift.t, 
   method = c("euler", "milstein", "KPS", "milstein2", 
   "cdist","ozaki","shoji","EA"), 
   alpha = 0.5, eta = 0.5, pred.corr = T, rcdist = NULL, 
   theta = NULL, model = c("CIR", "VAS", "OU", "BS"),
   k1, k2, phi, max.psi = 1000, rh, A, M=1)

Value

x: returns an invisible ts object

Arguments

t0: time origin.
T: horizon of simulation.
X0: initial value of the process.
N: number of simulation steps.
M: number of trajectories.
delta: time step of the simulation.
drift: drift coefficient: an expression of two variables t and x.
sigma: diffusion coefficient: an expression of two variables t and x.
drift.x: partial derivative of the drift coefficient w.r.t. x: a function of two variables t and x.
sigma.x: partial derivative of the diffusion coefficient w.r.t. x: a function of two variables t and x.
drift.xx: second partial derivative of the drift coefficient w.r.t. x: a function of two variables t and x.
sigma.xx: second partial derivative of the diffusion coefficient w.r.t. x: a function of two variables t and x.
drift.t: partial derivative of the drift coefficient w.r.t. t: a function of two variables t and x.
method: method of simulation; see details.
alpha: weight alpha of the predictor-corrector scheme.
eta: weight eta of the predictor-corrector scheme.
pred.corr: boolean: whether to apply the predictor-correct adjustment; see details.
rcdist: a function that is a random number generator from the conditional distribution of the process; see details.
theta: vector of parameters for cdist; see details.
model: model from which to simulate; see details.
k1: lower bound for psi(x); see details.
k2: upper bound for psi(x); see details.
phi: the function psi(x) - k1.
max.psi: upper value of the support of psi to search for its maximum.
rh: the rejection function; see details.
A: A(x) is the integral of the drift between 0 and x.

Author

Stefano Maria Iacus

Details

The function returns a ts object of length N+1; i.e., X0 and the new N simulated values if M=1. For M>1, an mts (multidimensional ts object) is returned, which means that M independent trajectories are simulated. If the initial value X0 is not of the length M, the values are recycled in order to have an initial vector of the correct length. If delta is not specified, then delta = (T-t0)/N. If delta is specified, then N values of the solution of the sde are generated and the time horizon T is adjusted to be N * delta.

The function psi is psi(x) = 0.5*drift(x)^2 + 0.5*drift.x(x).

If any of drift.x, drift.xx, drift.t, sigma.x, and sigma.xx are not specified, then numerical derivation is attempted when needed.

If sigma is not specified, it is assumed to be the constant function 1.

The method of simulation can be one among: euler, KPS, milstein, milstein2, cdist, EA, ozaki, and shoji. No assumption on the coefficients or on cdist is checked: the user is responsible for using the right method for the process object of simulation.

The model is one among: CIR: Cox-Ingersoll-Ross, VAS: Vasicek, OU Ornstein-Uhlenbeck, BS: Black and Scholes. No assumption on the coefficient theta is checked: the user is responsible for using the right ones.

If the method is cdist, then the process is simulated according to its known conditional distribution. The random generator rcdist must be a function of n, the number of random numbers; dt, the time lag; x, the value of the process at time t - dt; and the vector of parameters theta.

For the exact algorithm method EA: if missing k1 and k2 as well as A, rh and phi are calculated numerically by the function.

References

See Chapter 2 of the text.

Examples

Run this code

# Ornstein-Uhlenbeck process
set.seed(123)
d <- expression(-5 * x)
s <- expression(3.5) 
sde.sim(X0=10,drift=d, sigma=s) -> X
plot(X,main="Ornstein-Uhlenbeck")

# Multiple trajectories of the O-U process
set.seed(123)
sde.sim(X0=10,drift=d, sigma=s, M=3) -> X
plot(X,main="Multiple trajectories of O-U")

# Cox-Ingersoll-Ross process
# dXt = (6-3*Xt)*dt + 2*sqrt(Xt)*dWt
set.seed(123)
d <- expression( 6-3*x ) 
s <- expression( 2*sqrt(x) ) 
sde.sim(X0=10,drift=d, sigma=s) -> X
plot(X,main="Cox-Ingersoll-Ross")

# Cox-Ingersoll-Ross using the conditional distribution "rcCIR"

set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), rcdist=rcCIR, 
        method="cdist") -> X
plot(X, main="Cox-Ingersoll-Ross")

set.seed(123)
sde.sim(X0=10, theta=c(6, 3, 2), model="CIR") -> X
plot(X, main="Cox-Ingersoll-Ross")

# Exact simulation
set.seed(123)
d <- expression(sin(x))
d.x <- expression(cos(x)) 
A <- function(x) 1-cos(x)
sde.sim(method="EA", delta=1/20, X0=0, N=500, 
        drift=d, drift.x = d.x, A=A) -> X
plot(X, main="Periodic drift")

Run the code above in your browser using DataLab