mixedsde.sim {mixedsde}

Simulation Of A Mixed Stochastic Differential Equation

Description

Simulation of M independent trajectories of a mixed stochastic differential equation (SDE) with linear drift and two random effects (\alpha_j, \beta_j) dX_j(t)= (\alpha_j- \beta_j X_(t))dt + \sigma a(X_j(t)) dW_j(t), for j=1, ..., M.

Usage

mixedsde.sim(M, T, N = 100, model, random, fixed = 0, density.phi,
  param, sigma, t0 = 0, X0 = 0.01, invariant = 0, delta = T/N,
  op.plot = 0, add.plot = FALSE)

Arguments

M

number of trajectories

T

horizon of simulation.

N

number of simulation steps, default Tx100.

model

name of the SDE: 'OU' (Ornstein-Uhlenbeck) or 'CIR' (Cox-Ingersoll-Ross).

random

random effects in the drift: 1 if one additive random effect, 2 if one multiplicative random effect or c(1,2) if 2 random effects.

fixed

fixed effects in the drift: value of the fixed effect when there is only one random effect, 0 otherwise. If random =2, fixed can be 0 but \beta has to be a non negative random variable for the estimation.

density.phi

name of the density of the random effects.

param

vector of parameters of the distribution of the two random effects.

sigma

diffusion parameter

t0

time origin, default 0.

X0

initial value of the process, default X0=0.

invariant

1 if the initial value is simulated from the invariant distribution, default 0.01 and X0 is fixed.

delta

time step of the simulation (T/N).

op.plot

1 if a plot of the trajectories is required, default 0.

add.plot

1 for add trajectories to an existing plot

Details

Simulation of M independent trajectories of the SDE (the Brownian motions Wj are independent), with linear drift. Two diffusions are implemented, with one or two random effects:

Ornstein-Uhlenbeck model (OU)

If random = 1, \beta is a fixed effect: dX_j(t)= (\alpha_j- \beta X_j(t))dt + \sigma dW_j(t)

If random = 2, \alpha is a fixed effect: dX_j(t)= (\alpha - \beta_j X_j(t))dt + \sigma dW_j(t)

If random = c(1,2), dX_j(t)= (\alpha_j- \beta_j X_j(t))dt + \sigma dW_j(t)

Cox-Ingersoll-Ross model (CIR)

If random = 1, \beta is a fixed effect: dX_j(t)= (\alpha_j- \beta X_j(t))dt + \sigma \sqrt{X_j(t)} dW_j(t)

If random = 2, \alpha is a fixed effect: dX_j(t)= (\alpha - \beta_j X_j(t))dt + \sigma \sqrt{X_j(t)} dW_j(t)

If random = c(1,2), dX_j(t)= (\alpha_j- \beta_j X_j(t))dt + \sigma \sqrt{X_j(t)} dW_j(t)

The initial value of each trajectory can be simulated from the invariant distribution of the process: Normal distribution with mean \alpha/\beta and variance \sigma^2/(2 \beta) for the OU, a gamma distribution \Gamma(2\alpha/\sigma^2, \sigma^2/(2\beta)) for the C-I-R model.

Density of the random effects

Several densities are implemented for the random effects, depending on the number of random effects.

If two random effects, choice between

'normalnormal': Normal distributions for both \alpha \beta and param=c(mean_\alpha, sd_\alpha, mean_\beta, sd_\beta)

'gammagamma': Gamma distributions for both \alpha \beta and param=c(shape_\alpha, scale_\alpha, shape_\beta, scale_\beta)

'gammainvgamma': Gamma for \alpha, Inverse Gamma for \beta and param=c(shape_\alpha, scale_\alpha, shape_\beta, scale_\beta)

'normalgamma': Normal for \alpha, Gamma for \beta and param=c(mean_\alpha, sd_\alpha, shape_\beta, scale_\beta)

'normalinvgamma': Normal for \alpha, Inverse Gamma for \beta and param=c(mean_\alpha, sd_\alpha, shape_\beta, scale_\beta)

'gammagamma2': Gamma +2 * \sigma^2 for \alpha, Gamma + 1 for \beta and param=c(shape_\alpha, scale_\alpha, shape_\beta, scale_\beta)

'gammainvgamma2': Gamma +2 * \sigma^2 for \alpha, Inverse Gamma for \beta and param=c(shape_\alpha, scale_\alpha, shape_\beta, scale_\beta)

If only \alpha is random, choice between

'normal': Normal distribution with param=c(mean, sd)

lognormal': logNormal distribution with param=c(mean, sd)

'mixture.normal': mixture of normal distributions p N(\mu1,\sigma1^2) + (1-p)N(\mu2, \sigma2^2) with param=c(p, \mu1, \sigma1, \mu2, \sigma2)

'gamma': Gamma distribution with param=c(shape, scale)

'mixture.gamma': mixture of Gamma distribution p \Gamma(shape1,scale1) + (1-p)\Gamma(shape2,scale2) with param=c(p, shape1, scale1, shape2, scale2)

'gamma2': Gamma distribution +2 * \sigma^2 with param=c(shape, scale)

'mixed.gamma2': mixture of Gamma distribution p \Gamma(shape1,scale1) + (1-p) \Gamma(shape2,scale2) + +2 * \sigma^2 with param=c(p, shape1, scale1, shape2, scale2)

If only \beta is random, choice between 'normal': Normal distribution with param=c(mean, sd)

'gamma': Gamma distribution with param=c(shape, scale)

'mixture.gamma': mixture of Gamma distribution p \Gamma(shape1,scale1) + (1-p) \Gamma(shape2,scale2) with param=c(p, shape1, scale1, shape2, scale2)

Value

X

matrix (M x (N+1)) of the M trajectories.

phi

vector (or matrix) of the M simulated random effects.

References

This function mixedsde.sim is based on the package sde, function sde.sim. See Simulation and Inference for stochastic differential equation, S.Iacus, Springer Series in Statistics 2008 Chapter 2

See Also

https://cran.r-project.org/package=sde

Examples

#Simulation of 5 trajectories of the OU SDE with random =1, and a Gamma distribution.

simuOU <- mixedsde.sim(M=5, T=10,N=1000,model='OU', random=1,fixed=0.5,
density.phi='gamma', param=c(1.8, 0.8) , sigma=0.1,op.plot=1)
X <- simuOU$X ; 
phi <- simuOU$phi
hist(phi)