sde.sim {sde}

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)

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.

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.

Value

x

returns an invisible ts object

Author(s)

Stefano Maria Iacus

References

See Chapter 2 of the text.

Examples

# 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")