R Graphical Manual

Browse All

Last data update: 2014.03.03

R: First Passage Time in 3-Dim SDE

fptsde3d

R Documentation

First Passage Time in 3-Dim SDE

Description

The (S3) generic function fptsde3d for simulate first-passage-time (f.p.t) in 3-dim stochastic differential equations.

Usage

fptsde3d(N, ...)
## Default S3 method:
fptsde3d(N = 1000, M = 100, x0 = 0, y0 = 0, z0 = 0, t0 = 0, T = 1, Dt, 
   boundary, driftx, diffx, drifty, diffy, driftz, diffz, alpha = 0.5, mu = 0.5, 
   type = c("ito", "str"), method = c("euler", "milstein", "predcorr", 
   "smilstein", "taylor", "heun", "rk1", "rk2", "rk3"), ...)
						   
						   
## S3 method for class 'fptsde3d'
summary(object, ...)	
## S3 method for class 'fptsde3d'
mean(x, ...)
## S3 method for class 'fptsde3d'
median(x, ...)
## S3 method for class 'fptsde3d'
quantile(x, ...)
## S3 method for class 'fptsde3d'
kurtosis(x, ...)
## S3 method for class 'fptsde3d'
skewness(x, ...)
## S3 method for class 'fptsde3d'
moment(x, order = 2, ...)
## S3 method for class 'fptsde3d'
bconfint(x, level=0.95, ...)	
## S3 method for class 'fptsde3d'
plot(x, ...)

Arguments

`N`	size of sde.
`M`	size of fpt.
`x0, y0, z0`	initial value of the process X(t), Y(t) and Z(t) at time `t0`.
`t0`	initial time.
`T`	final time.
`Dt`	time step of the simulation (discretization). If it is `missing` a default Dt = (T-t0)/N.
`boundary`	an `expression` of a constant or time-dependent boundary.
`driftx, drifty, driftz`	drift coefficient: an `expression` of four variables `t`, `x`, `y` and `z` for process X(t), Y(t) and Z(t).
`diffx, diffy, diffz`	diffusion coefficient: an `expression` of four variables `t`, `x`, `y` and `z` for process X(t), Y(t) and Z(t).
`alpha, mu`	weight of the predictor-corrector scheme; the default `alpha = 0.5` and `mu = 0.5`.
`type`	sde of the type Ito or Stratonovich.
`method`	numerical methods of simulation, the default `method = "euler"`; see `snssde3d`.
`x, object`	an object inheriting from class `"fptsde3d"`.
`order`	order of moment.
`level`	the confidence level required.
`...`	further arguments for (non-default) methods.

Details

The function fptsde3d returns a random variable (tau(X(t),S(t)),tau(Y(t),S(t)),tau(Z(t),S(t))) "first passage time", is defined as :

tau(X(t),S(t))={t>=0; X(t) >= S(t)}, if X(t0) < S(t0)

tau(Y(t),S(t))={t>=0; Y(t) >= S(t)}, if Y(t0) < S(t0)

tau(Z(t),S(t))={t>=0; Z(t) >= S(t)}, if Z(t0) < S(t0)

and:

tau(X(t),S(t))={t>=0; X(t) <= S(t)}, if X(t0) > S(t0)

tau(Y(t),S(t))={t>=0; Y(t) <= S(t)}, if Y(t0) > S(t0)

tau(Z(t),S(t))={t>=0; Z(t) <= S(t)}, if Z(t0) > S(t0)

with S(t) is through a continuous boundary (barrier).

Value

fptsde3d returns an object inheriting from class "fptsde3d".

fptx, fpty, fptz

a vector of triplet 'fpt' (tau(X(t),S(t)),tau(Y(t),S(t)),tau(Z(t),S(t))).

Author(s)

A.C. Guidoum, K. Boukhetala.

References

Argyrakisa, P. and G.H. Weiss (2006). A first-passage time problem for many random walkers. Physica A. 363, 343–347.

Aytug H., G. J. Koehler (2000). New stopping criterion for genetic algorithms. European Journal of Operational Research, 126, 662–674.

Boukhetala, K. (1996) Modelling and simulation of a dispersion pollutant with attractive centre. ed by Computational Mechanics Publications, Southampton ,U.K and Computational Mechanics Inc, Boston, USA, 245–252.

Boukhetala, K. (1998a). Estimation of the first passage time distribution for a simulated diffusion process. Maghreb Math.Rev, 7(1), 1–25.

Boukhetala, K. (1998b). Kernel density of the exit time in a simulated diffusion. les Annales Maghrebines De L ingenieur, 12, 587–589.

Ding, M. and G. Rangarajan. (2004). First Passage Time Problem: A Fokker-Planck Approach. New Directions in Statistical Physics. ed by L. T. Wille. Springer. 31–46.

Roman, R.P., Serrano, J. J., Torres, F. (2008). First-passage-time location function: Application to determine first-passage-time densities in diffusion processes. Computational Statistics and Data Analysis. 52, 4132–4146.

Roman, R.P., Serrano, J. J., Torres, F. (2012). An R package for an efficient approximation of first-passage-time densities for diffusion processes based on the FPTL function. Applied Mathematics and Computation, 218, 8408–8428.

Gardiner, C. W. (1997). Handbook of Stochastic Methods. Springer-Verlag, New York.

Examples


## dX(t) = 4*(-1-X(t))*Y(t) dt + 0.2 * dW1(t) 
## dY(t) = 4*(1-Y(t)) *X(t) dt + 0.2 * dW2(t) 
## dZ(t) = 4*(1-Z(t)) *Y(t) dt + 0.2 * dW3(t) 
## x0 = 0, y0 = -2, z0 = 0, and barrier -3+5*t.       
## W1(t), W2(t) and W3(t) three independent Brownian motion      
set.seed(1234)

fx <- expression(4*(-1-x)*y)
gx <- expression(0.2)
fy <- expression(4*(1-y)*x)
gy <- expression(0.2)
fz <- expression(4*(1-z)*y)
gz <- expression(0.2)

St <- expression(-3+5*t)

res <- fptsde3d(driftx=fx,diffx=gx,drifty=fy,diffy=gy,driftz=fz,diffz=gz,boundary=St,
                x0=2,y0=-2,z0=0,M=50)
res
summary(res)
plot(res,union=TRUE)
dev.new()
plot(res,union=FALSE)
##

fptx <- res$fptx
fpty <- res$fpty
fptz <- res$fptz
X1 <- cbind(fptx,fpty,fptz)
## library(sm)
## sm.density(X1,display="rgl")

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Sim.DiffProc)
Package 'Sim.DiffProc' version 3.2 loaded.
help(Sim.DiffProc) for summary information.
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Sim.DiffProc/fptsde3d.Rd_%03d_medium.png", width=480, height=480)
> ### Name: fptsde3d
> ### Title: First Passage Time in 3-Dim SDE
> ### Aliases: fptsde3d fptsde3d.default summary.fptsde3d mean.fptsde3d
> ###   median.fptsde3d quantile.fptsde3d kurtosis.fptsde3d skewness.fptsde3d
> ###   moment.fptsde3d bconfint.fptsde3d plot.fptsde3d
> ### Keywords: fpt sde ts mts
> 
> ### ** Examples
> 
> 
> ## dX(t) = 4*(-1-X(t))*Y(t) dt + 0.2 * dW1(t) 
> ## dY(t) = 4*(1-Y(t)) *X(t) dt + 0.2 * dW2(t) 
> ## dZ(t) = 4*(1-Z(t)) *Y(t) dt + 0.2 * dW3(t) 
> ## x0 = 0, y0 = -2, z0 = 0, and barrier -3+5*t.       
> ## W1(t), W2(t) and W3(t) three independent Brownian motion      
> set.seed(1234)
> 
> fx <- expression(4*(-1-x)*y)
> gx <- expression(0.2)
> fy <- expression(4*(1-y)*x)
> gy <- expression(0.2)
> fz <- expression(4*(1-z)*y)
> gz <- expression(0.2)
> 
> St <- expression(-3+5*t)
> 
> res <- fptsde3d(driftx=fx,diffx=gx,drifty=fy,diffy=gy,driftz=fz,diffz=gz,boundary=St,
+                 x0=2,y0=-2,z0=0,M=50)
> res
$SDE
Ito Sde 3D:
	| dX(t) = 4 * (-1 - X(t)) * Y(t) * dt + 0.2 * dW1(t)
	| dY(t) = 4 * (1 - Y(t)) * X(t) * dt + 0.2 * dW2(t)
	| dZ(t) = 4 * (1 - Z(t)) * Y(t) * dt + 0.2 * dW3(t)
Method:
	| Euler scheme of order 0.5
Summary:
	| Size of process	| N  = 1000.
	| Number of simulation	| M  = 50.
	| Initial values	| (x0,y0,z0) = (2,-2,0).
	| Time of process	| t in [0,1].
	| Discretization	| Dt = 0.001.

$boundary
-3 + 5 * t

$fptx
 [1] 0.5760063 0.5832781 0.5795617 0.5931902 0.5791569 0.5725338 0.5762656
 [8] 0.5878879 0.5679412 0.5874620 0.5781744 0.5718304 0.5827081 0.5722382
[15] 0.5736284 0.5814066 0.5752476 0.5732806 0.5885270 0.5769123 0.5797888
[22] 0.5698742 0.5997994 0.5853098 0.5834944 0.5790055 0.5773006 0.5910663
[29] 0.5921050 0.5846778 0.5817661 0.5902033 0.5753017 0.5814413 0.5942903
[36] 0.5757611 0.5882927 0.5773588 0.5804145 0.5868764 0.5779249 0.5785095
[43] 0.5887054 0.5850905 0.5842715 0.5868127 0.5682496 0.5871819 0.5863559
[50] 0.5770388

$fpty
 [1] 0.7994094 0.8226149 0.7468067 0.7998463 0.7522640 0.7917487 0.7860098
 [8] 0.7864210 0.7175875 0.8348322 0.7857938 0.8550556 0.7894133 0.8169060
[15] 0.8655540 0.7847123 0.8178045 0.9250828 0.7887015 0.7801047 0.7352194
[22] 0.7425373 0.8022543 0.7466123 0.7723577 0.7676493 0.8106421 0.7359943
[29] 0.8069275 0.7889495 0.8233341 0.7755564 0.8140864 0.7373236 0.7908124
[36] 0.7547203 0.7902368 0.7726842 0.7478743 0.7457673 0.7989163 0.7725229
[43] 0.7885825 0.7404833 0.7895216 0.7766209 0.7539435 0.8036853 0.8253030
[50] 0.7890954

$fptz
 [1] 0.7861733 0.7956446 0.7414119 0.7785966 0.7476102 0.7843249 0.7773501
 [8] 0.7595979 0.7350067 0.7966396 0.7819250 0.8029107 0.7654809 0.7923767
[15] 0.8033168 0.7689548 0.7871758 0.8213037 0.7660745 0.7668914 0.7354981
[22] 0.7460324 0.7802632 0.7442508 0.7655039 0.7509376 0.7836662 0.7207627
[29] 0.7691037 0.7638728 0.7935643 0.7559995 0.7927764 0.7277230 0.7775181
[36] 0.7551167 0.7709414 0.7624840 0.7448931 0.7368079 0.7759685 0.7651244
[43] 0.7788589 0.7298758 0.7488274 0.7532690 0.7640534 0.7851165 0.7916897
[50] 0.7671621

attr(,"class")
[1] "fptsde3d"
> summary(res)

	Monte-Carlo Statistics for the F.P.T of (X(t),Y(t),Z(t))
| T(S,X) = inf{t >=  0  : X(t) <=  -3 + 5 * t }
| T(S,Y) = inf{t >=  0  : Y(t) <=  -3 + 5 * t }
| T(S,Z) = inf{t >=  0  : Z(t) <=  -3 + 5 * t }

                       fpt(x)   fpt(y)    fpt(z)
NA's                        0        0         0
Mean                 0.581430 0.786938  0.767929
Variance             0.000051 0.001400  0.000492
Median               0.580911 0.788642  0.767027
First quartile       0.576427 0.757953  0.751520
Third quartile       0.586860 0.803328  0.784160
Skewness             0.223249 0.956670 -0.050127
Kurtosis             2.495032 5.137701  2.410876
Moment of order 2    0.000050 0.001372  0.000482
Moment of order 3    0.000000 0.000050 -0.000001
Moment of order 4    0.000000 0.000010  0.000001
Moment of order 5    0.000000 0.000001  0.000000
Bound conf Inf (95%) 0.568615 0.735394  0.728207
Bound conf Sup (95%) 0.594043 0.863192  0.803225
> plot(res,union=TRUE)
> dev.new()
Error in dev.new() : no suitable unused file name for pdf()
Execution halted