## Tue Mar 5 18:13:25 2024
## Original file Copyright © 2024 A.C. Guidoum, K. Boukhetala
## This file is part of the R package Sim.DiffProc
## Department of Mathematics and Computer Science,
## Faculty of Sciences and Technology,
## University of Tamanghasset
## and
## Department of Probabilities & Statistics
## Faculty of Mathematics
## University of Science and Technology Houari Boumediene
## BP 32 El-Alia, U.S.T.H.B, Algiers
## Algeria
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
## A copy of the GNU General Public License is available at
## http://www.r-project.org/Licenses/
## Unlimited use and distribution (see LICENCE).
###################################################################################################
#####
##### bridgesde1d
bridgesde1d <- function(N, ...) UseMethod("bridgesde1d")
bridgesde1d.default <- function(N =1000,M=1,x0=0,y=0,t0=0,T=1,Dt,drift,diffusion,
alpha=0.5,mu=0.5,type=c("ito","str"), method=c(
"euler","milstein","predcorr","smilstein","taylor",
"heun","rk1","rk2","rk3"),...)
{
if (!is.numeric(x0) || !is.numeric(y))
stop("'x0' and 'y' must be numeric")
if (any(!is.numeric(t0) || !is.numeric(T)))
stop(" 't0' and 'T' must be numeric")
if (any(!is.numeric(N) || (N - floor(N) > 0) || N <= 1))
stop(" 'N' must be a positive integer ")
if (any(!is.numeric(M) || (M - floor(M) > 0) || M <= 0))
stop(" 'M' must be a positive integer ")
if (any(!is.expression(drift) || !is.expression(diffusion) ))
stop(" coefficient of 'drift' and 'diffusion' must be expressions in 't' and 'x'")
if (missing(type)) type <- "ito"
method <- match.arg(method)
if (method =="predcorr"){
if (any(alpha > 1 || alpha < 0))
stop("please use '0 <= alpha <= 1' ")
if (any(mu > 1 || mu < 0))
stop("please use '0 <= mu <= 1' ")
}
if (t0 < 0 || T < 0)
stop(" please use positive times! (0 <= t0 < T) ")
if (missing(Dt)) {
t <- seq(t0, T, by=Dt)
} else {
t <- c(t0, t0 + cumsum(rep(Dt, N)))
T <- t[N + 1]
}
Dt <- (T - t0)/N
X1 <- snssde1d(N,M,x0,t0,T,Dt,drift,diffusion,alpha,mu,type, method,...)$X
if (M > 1){X2 <- apply(data.frame(snssde1d(N,M,x0=y,t0,T,Dt,drift,diffusion,alpha,mu,type, method,...)$X),2,rev)
}else{
X2 <- rev(snssde1d(N,M,x0=y,t0,T,Dt,drift,diffusion,alpha,mu,type, method,...)$X)}
G <- rep(NA,M)
if (M > 1){
for(j in 1:M){
if (X1[1,M] >= X2[1,M]){
if (!all(X1[,j] > X2[,j]))
G[j] <- min(which((X1[,j]-X2[,j]) <= 0)) - 1
}else{ if (!all(X1[,j] < X2[,j]))
G[j] <- min(which((X1[,j]-X2[,j]) >= 0)) - 1
}
}
}else{
if (X1[1] >= X2[1]){
if (!all(X1 > X2))
G <- min(which((X1-X2) <= 0)) - 1
}else{ if (!all(X1 < X2))
G <- min(which((X1-X2) >= 0)) - 1
}
}
G[which(G==0)]=NA
name <- "X"
name <- if(M > 1)
paste("X",1:length(which(!is.na(G))),sep="")
if (M == 1){
if (is.na(G) ){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else{X <- ts(c(X1[1:G],X2[-(1:G)]),start=t0,deltat=deltat(X1),names=name)}
}else if (M > 1){ if(length(which(is.na(G))) == M){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else if (length(which(!is.na(G))) == 1){X <- ts(c(X1[,-which(is.na(G))][1:G[which(!is.na(G))]],X2[,-which(is.na(G))][-(1:G[which(!is.na(G))])]),start=t0,deltat=deltat(X1),names=name)
}else if (length(which( is.na(G))) == 0){X <- ts(sapply(1:length(which(!is.na(G))), function(j) c(X1[,j][1:G[!is.na(G)][j]],X2[,j][-(1:G[!is.na(G)][j])])),start=t0,deltat=deltat(X1),names=name)
}else{ X1 <- X1[,-c(which(is.na(G)))]
X2 <- X2[,-c(which(is.na(G)))]
X <- ts(sapply(1:length(which(!is.na(G))), function(j) c(X1[,j][1:G[!is.na(G)][j]],X2[,j][-(1:G[!is.na(G)][j])])),
start=t0,deltat=deltat(X1),names=name)
}
}
structure(list(X=X,drift=drift[[1]], diffusion=diffusion[[1]],type=type,method=method,
x0=x0,y=y, N=N,M=M,Dt=Dt,t0=t0,T=T,C=G,dim="1d",call=match.call()),class="bridgesde1d")
}
###
print.bridgesde1d <- function(x, digits=NULL, ...)
{
class(x) <- "bridgesde1d"
Ito = "It\xf4"
Encoding(Ito) <- "latin1"
if (x$method=="euler")
sch <- "Euler scheme with order 0.5"
else if (x$method=="milstein")
sch <- "First-order Milstein scheme"
else if (x$method=="predcorr")
sch <- "Predictor-corrector method with order 1"
else if (x$method=="smilstein")
sch <- "Second-order Milstein scheme"
else if (x$method=="taylor")
sch <- "Taylor scheme with order 1.5"
else if (x$method=="heun")
sch <- "Heun scheme with order 2"
else if (x$method=="rk1")
sch <- "Runge-Kutta method with order 1"
else if (x$method=="rk2")
sch <- "Runge-Kutta method with order 2"
else if (x$method=="rk3")
sch <- "Runge-Kutta method with order 3"
##Dr <- gsub(pattern = 'x', replacement = 'X(t)', x = as.expression(x$drift), ignore.case = F,fixed = T)
##DD <- gsub(pattern = 'x', replacement = 'X(t)', x = as.expression(x$diffusion), ignore.case = F,fixed = T)
Dr <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)))), list(e = x$drift))))
DD <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)))), list(e = x$diffusion))))
if(x$type=="ito"){
cat(Ito," Bridge Sde 1D:","\n",
"\t| dX(t) = ", Dr," * dt + ", DD," * dW(t)","\n",
"Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$C))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial value","\t\t| x0 = ",format(x$x0,digits=digits),".","\n",
"\t| Ending value","\t\t| y = ",format(x$y,digits=digits),".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}else{
cat("Stratonovich Bridge Sde 1D:","\n",
"\t| dX(t) = ", Dr," * dt + ", DD," o dW(t)","\n",
"Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$C))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial value","\t\t| x0 = ",format(x$x0,digits=digits),".","\n",
"\t| Ending value","\t\t| y = ",format(x$y,digits=digits),".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}
invisible(x)
}
time.bridgesde1d <- function(x,...)
{
class(x) <- "bridgesde1d"
as.vector(time(x$X))
}
summary.bridgesde1d <- function(object, at,digits=NULL, ...)
{
class(object) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (is.null(digits)){digits = base::options()$digits}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$C)))
if (C == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)}else{X = object$X}
xx <- as.vector(X[which(time(object)==as.character(at)),])
if (length(xx) == 0){
if (C == 1){ F <- lapply(1:C,function(i) approxfun(time(object),object$X))}else{
F <- lapply(1:C,function(i) approxfun(time(object),object$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
cat("\nMonte-Carlo Statistics for X(t) at time t = ",at,"\n",
"\t| Crossing realized ",C," among ",object$M,"\n",
sep="")
res <- as.data.frame(matrix(c(mean(xx,na.rm = TRUE),var(xx,na.rm = TRUE),Median(xx),Mode(xx),
quantile(xx,0.25,na.rm = TRUE),quantile(xx,0.75,na.rm = TRUE),min(xx,na.rm=TRUE),
max(xx,na.rm=TRUE),skewness(xx),kurtosis(xx),cv(xx),moment(xx,order=3,center=TRUE),
moment(xx,order=4,center=TRUE),moment(xx,order=5,center=TRUE),moment(xx,order=6,center=TRUE)),
ncol=1))
dimnames(res) <- list(c("Mean","Variance","Median","Mode","First quartile","Third quartile","Minimum","Maximum","Skewness","Kurtosis",
"Coef-variation","3th-order moment","4th-order moment","5th-order moment","6th-order moment"),c(""))
print(round(res,digits=digits), quote = FALSE, right = TRUE,...)
invisible(object)
}
mean.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
mean(xx,na.rm = TRUE,...)
}
cv.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
cv(xx,...)
}
min.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
min(xx,na.rm = TRUE)
}
max.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
max(xx,na.rm = TRUE)
}
skewness.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
skewness(xx)
}
kurtosis.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
kurtosis(xx)
}
Median.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
Median(xx)
}
Mode.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
Mode(xx)
}
quantile.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
quantile(xx,...)
}
moment.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
moment(xx,...)
}
bconfint.bridgesde1d <- function(x,at,...)
{
class(x) <- "bridgesde1d"
if (missing(at)) {at = (as.numeric(x$T)-as.numeric(x$t0))/2}
if (any(x$T < at || x$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(x$C)))
if (C == 1){ X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
xx <- as.vector(X[which(time(x)==as.character(at)),])
if (length(xx) == 0){
if (C==1){ F <- lapply(1:C,function(i) approxfun(time(x),x$X))}else{
F <- lapply(1:C,function(i) approxfun(time(x),x$X[,i]))}
xx <- sapply(1:length(F),function(i) F[[i]](at))
}
bconfint(xx,...)
}
##
## Plot
plot.bridgesde1d <- function(x,...)
{
class(x) <- "bridgesde1d"
X <- x$X
if (x$M > 10) {plot(X,plot.type="single",...)}else{plot(X,...)}
}
lines.bridgesde1d <- function(x,...)
{
class(x) <- "bridgesde1d"
X <- x$X
if (length(which(!is.na(x$C))) >=2){
for (i in 1:dim(X)[2]){
lines(time(x),X[,i],...)}}else{
lines(time(x),X,...)}
}
points.bridgesde1d <- function(x,...)
{
class(x) <- "bridgesde1d"
X <- x$X
if (length(which(!is.na(x$C))) >=2){
for (i in 1:dim(X)[2]){
points(time(x),X[,i],...)}}else{
points(time(x),X,...)}
}
#####
##### bridgesde2d
bridgesde2d <- function(N, ...) UseMethod("bridgesde2d")
bridgesde2d.default <- function(N =1000,M=1,x0=c(0,0),y=c(0,0),t0=0,T=1,Dt,drift,diffusion,
corr = NULL,alpha=0.5,mu=0.5,type=c("ito","str"), method=c(
"euler","milstein","predcorr","smilstein","taylor",
"heun","rk1","rk2","rk3"),...)
{
if (!is.numeric(x0) || length(x0) != 2)
stop("'x0' must be numeric, and length(x0) = 2 ")
if (!is.numeric(y) || length(y) != 2)
stop("'y' must be numeric, and length(y) = 2 ")
if (any(!is.numeric(t0) || !is.numeric(T)))
stop(" 't0' and 'T' must be numeric")
if (any(!is.numeric(N) || (N - floor(N) > 0) || N <= 1))
stop(" 'N' must be a positive integer ")
if (any(!is.numeric(M) || (M - floor(M) > 0) || M <= 0))
stop(" 'M' must be a positive integer ")
if (length(drift) !=2 )
stop("drift must be expression 2d (vector of 2 expression)")
if (length(diffusion) !=2 )
stop("diffusion must be expression 2d (vector of 2 expression)")
if (any(!is.expression(drift) || !is.expression(diffusion) ))
stop(" coefficient of 'drift' and 'diffusion' must be expressions in 't', 'x' and 'y'")
if (!is.null(corr)) {
if (any(!is.matrix(corr) || det(corr) <= 0 || nrow(corr)!=2 || ncol(corr)!=2 || !isSymmetric(corr)))
stop("the correlation structure of W1(t) and W2(t) must be a real symmetric\n positive-definite square matrix of dimension 2")
}
if (missing(type)) type <- "ito"
method <- match.arg(method)
if (!is.null(corr) && method != "euler" && method != "milstein") {method="euler"}
if (method =="predcorr"){
if (any(alpha > 1 || alpha < 0))
stop("please use '0 <= alpha <= 1' ")
if (any(mu > 1 || mu < 0))
stop("please use '0 <= mu <= 1' ")
}
if (t0 < 0 || T < 0)
stop(" please use positive times! (0 <= t0 < T) ")
if (missing(Dt)) {
t <- seq(t0, T, by=Dt)
} else {
t <- c(t0, t0 + cumsum(rep(Dt, N)))
T <- t[N + 1]
}
Dt <- (T - t0)/N
X1 <- snssde2d(N,M,x0=c(x0[1],x0[2]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$X
Y1 <- snssde2d(N,M,x0=c(x0[1],x0[2]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Y
if (M > 1){X2 <- apply(data.frame(snssde2d(N,M,x0=c(y[1],y[2]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$X),2,rev)
Y2 <- apply(data.frame(snssde2d(N,M,x0=c(y[1],y[2]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Y),2,rev)
}else{X2 <- rev(snssde2d(N,M,x0=c(y[1],y[2]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$X)
Y2 <- rev(snssde2d(N,M,x0=c(y[1],y[2]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Y)
}
Gx = Gy <- rep(NA,M)
if (M > 1){
for(j in 1:M){
if (X1[1,M] >= X2[1,M]){
if (!all(X1[,j] > X2[,j]))
Gx[j] <- min(which((X1[,j]-X2[,j]) <= 0)) - 1
}else{ if (!all(X1[,j] < X2[,j]))
Gx[j] <- min(which((X1[,j]-X2[,j]) >= 0)) - 1
}
}
}else{
if (X1[1] >= X2[1]){
if (!all(X1 > X2))
Gx <- min(which((X1-X2) <= 0)) - 1
}else{ if (!all(X1 < X2))
Gx <- min(which((X1-X2) >= 0)) - 1
}
}
if (M > 1){
for(j in 1:M){
if (Y1[1,M] >= Y2[1,M]){
if (!all(Y1[,j] > Y2[,j]))
Gy[j] <- min(which((Y1[,j]-Y2[,j]) <= 0)) - 1
}else{ if (!all(Y1[,j] < Y2[,j]))
Gy[j] <- min(which((Y1[,j]-Y2[,j]) >= 0)) - 1
}
}
}else{
if (Y1[1] >= Y2[1]){
if (!all(Y1 > Y2))
Gy <- min(which((Y1-Y2) <= 0)) - 1
}else{ if (!all(Y1 < Y2))
Gy <- min(which((Y1-Y2) >= 0)) - 1
}
}
Gx[which(Gx==0)]=NA
Gy[which(Gy==0)]=NA
G <- na.omit(data.frame(Gx,Gy))
Gx <- G$Gx;Gy <- G$Gy
namex <- "X"
namey <- "Y"
namex <- if(M > 1) paste("X",1:length(which(!is.na(Gx))),sep="")
namey <- if(M > 1) paste("Y",1:length(which(!is.na(Gy))),sep="")
if (M == 1){ if (is.na(Gx) ){stop( "A crossing has been no realized,trying again (Repeat)..." )}else{X <- ts(c(X1[1:Gx],X2[-(1:Gx)]),start=t0,deltat=deltat(X1),names=namex)}
}else if (M > 1){ if(length(which(is.na(Gx))) == M){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else if (length(which(!is.na(Gx))) == 1){X <- ts(c(X1[,-which(is.na(Gx))][1:Gx[which(!is.na(Gx))]],X2[,-which(is.na(Gx))][-(1:Gx[which(!is.na(Gx))])]),start=t0,deltat=deltat(X1),names=namex)
}else if (length(which( is.na(Gx))) == 0){X <- ts(sapply(1:length(which(!is.na(Gx))), function(j) c(X1[,j][1:Gx[!is.na(Gx)][j]],X2[,j][-(1:Gx[!is.na(Gx)][j])])),start=t0,deltat=deltat(X1),names=namex)
}else{ X1 <- X1[,-c(which(is.na(Gx)))]
X2 <- X2[,-c(which(is.na(Gx)))]
X <- ts(sapply(1:length(which(!is.na(Gx))), function(j) c(X1[,j][1:Gx[!is.na(Gx)][j]],X2[,j][-(1:Gx[!is.na(Gx)][j])])),
start=t0,deltat=deltat(X1),names=namex)
}
}
if (M == 1){ if (is.na(Gy) ){stop( "A crossing has been no realized,trying again (Repeat)..." )}else{Y <- ts(c(Y1[1:Gy],Y2[-(1:Gy)]),start=t0,deltat=deltat(Y1),names=namey)}
}else if (M > 1){ if(length(which(is.na(Gy))) == M){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else if (length(which(!is.na(Gy))) == 1){Y <- ts(c(Y1[,-which(is.na(Gy))][1:Gy[which(!is.na(Gy))]],Y2[,-which(is.na(Gy))][-(1:Gy[which(!is.na(Gy))])]),start=t0,deltat=deltat(Y1),names=namey)
}else if (length(which( is.na(Gy))) == 0){Y <- ts(sapply(1:length(which(!is.na(Gy))), function(j) c(Y1[,j][1:Gy[!is.na(Gy)][j]],Y2[,j][-(1:Gy[!is.na(Gy)][j])])),start=t0,deltat=deltat(Y1),names=namey)
}else{ Y1 <- Y1[,-c(which(is.na(Gy)))]
Y2 <- Y2[,-c(which(is.na(Gy)))]
Y <- ts(sapply(1:length(which(!is.na(Gy))), function(j) c(Y1[,j][1:Gy[!is.na(Gy)][j]],Y2[,j][-(1:Gy[!is.na(Gy)][j])])),
start=t0,deltat=deltat(Y1),names=namey)
}
}
structure(list(X=X,Y=Y, driftx=drift[[1]], diffx=diffusion[[1]],drifty=drift[[2]], diffy=diffusion[[2]],corrmat=corr,type=type,method=method,
x0=x0,y=y, N=N,M=M,Dt=Dt,t0=t0,T=T,Cx=Gx,Cy=Gy,dim="2d",call=match.call()),class="bridgesde2d")
}
###
print.bridgesde2d <- function(x, digits=NULL, ...)
{
class(x) <- "bridgesde2d"
Ito = "It\xf4"
Encoding(Ito) <- "latin1"
if (x$method=="euler")
sch <- "Euler scheme with order 0.5"
else if (x$method=="milstein")
sch <- "First-order Milstein scheme"
else if (x$method=="predcorr")
sch <- "Predictor-corrector method with order 1"
else if (x$method=="smilstein")
sch <- "Second-order Milstein scheme"
else if (x$method=="taylor")
sch <- "Taylor scheme with order 1.5"
else if (x$method=="heun")
sch <- "Heun scheme with order 2"
else if (x$method=="rk1")
sch <- "Runge-Kutta method with order 1"
else if (x$method=="rk2")
sch <- "Runge-Kutta method with order 2"
else if (x$method=="rk3")
sch <- "Runge-Kutta method with order 3"
Drx <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)))), list(e = x$driftx))))
DDx <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)))), list(e = x$diffx))))
Dry <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)))), list(e = x$drifty))))
DDy <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)))), list(e = x$diffy))))
if (is.null(x$corrmat)){
if(x$type=="ito"){
cat(Ito," Bridge Sde 2D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," * dW1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," * dW2(t)","\n",
"Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}else{
cat("Stratonovich Bridge Sde 2D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," o dW1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," o dW2(t)","\n",
"Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}} else {
if(x$type=="ito"){
cat(Ito," Bridge Sde 2D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," * dB1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," * dB2(t)","\n",
"\t| Correlation structure:",sep="")
prmatrix(x$corrmat,rowlab = rep(" ", 2), collab = rep("", 2),digits=digits)
cat("Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}else{
cat("Stratonovich Bridge Sde 2D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," o dB1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," o dB2(t)","\n",
"\t| Correlation structure:",sep="")
prmatrix(x$corrmat,rowlab = rep(" ", 2), collab = rep("", 2),digits=digits)
cat("Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}
}
invisible(x)
}
##
## Plot
plot.bridgesde2d <- function(x,...) .plot.bridgesde2d(x,...)
lines.bridgesde2d <- function(x,...)
{
class(x) <- "bridgesde2d"
if (length(which(!is.na(x$Cx))) == 1){X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
if (length(which(!is.na(x$Cy))) == 1){Y = matrix(x$Y,nrow=length(x$Y),ncol=1)}else{Y = x$Y}
for (i in 1:length(which(!is.na(x$Cx)))){lines(time(x),X[,i],...)}
for (i in 1:length(which(!is.na(x$Cy)))){lines(time(x),Y[,i],...)}
}
points.bridgesde2d <- function(x,...)
{
class(x) <- "bridgesde2d"
if (length(which(!is.na(x$Cx))) == 1){X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
if (length(which(!is.na(x$Cy))) == 1){Y = matrix(x$Y,nrow=length(x$Y),ncol=1)}else{Y = x$Y}
for (i in 1:length(which(!is.na(x$Cx)))){points(time(x),X[,i],...)}
for (i in 1:length(which(!is.na(x$Cy)))){points(time(x),Y[,i],...)}
}
plot2d.bridgesde2d <- function(x,...) .plot2d.bridgesde2d(x,...)
lines2d.bridgesde2d <- function(x,...)
{
class(x) <- "bridgesde2d"
if (length(which(!is.na(x$Cx))) == 1){X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
if (length(which(!is.na(x$Cy))) == 1){Y = matrix(x$Y,nrow=length(x$Y),ncol=1)}else{Y = x$Y}
X <- X[,1]
Y <- Y[,1]
lines2d(X,Y,...)
}
points2d.bridgesde2d <- function(x,...)
{
class(x) <- "bridgesde2d"
if (length(which(!is.na(x$Cx))) == 1){X = matrix(x$X,nrow=length(x$X),ncol=1)}else{X = x$X}
if (length(which(!is.na(x$Cy))) == 1){Y = matrix(x$Y,nrow=length(x$Y),ncol=1)}else{Y = x$Y}
X <- X[,1]
Y <- Y[,1]
points2d(X,Y,...)
}
##
## summary
summary.bridgesde2d <- function(object,at,digits=NULL,...)
{
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (is.null(digits)){digits = base::options()$digits}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
cat("\nMonte-Carlo Statistics for (X(t),Y(t)) at time t = ",at,"\n",
"\t| Crossing realized ",C," among ",object$M,"\n",
sep="")
res <- as.data.frame(matrix(c(mean(x,na.rm = TRUE),var(x,na.rm = TRUE),Median(x),Mode(x),
quantile(x,0.25,na.rm = TRUE),quantile(x,0.75,na.rm = TRUE),min(x,na.rm=TRUE),
max(x,na.rm=TRUE),skewness(x),kurtosis(x),cv(x),moment(x,order=3,center=TRUE),
moment(x,order=4,center=TRUE),moment(x,order=5,center=TRUE),moment(x,order=6,center=TRUE),
mean(y,na.rm = TRUE),var(y,na.rm = TRUE),Median(y),Mode(y),
quantile(y,0.25,na.rm = TRUE),quantile(y,0.75,na.rm = TRUE),min(y,na.rm=TRUE),
max(y,na.rm=TRUE),skewness(y),kurtosis(y),cv(y),moment(y,order=3,center=TRUE),
moment(y,order=4,center=TRUE),moment(y,order=5,center=TRUE),moment(y,order=6,center=TRUE)),
ncol=2))
dimnames(res) <- list(c("Mean","Variance","Median","Mode","First quartile","Third quartile","Minimum","Maximum","Skewness","Kurtosis",
"Coef-variation","3th-order moment","4th-order moment","5th-order moment","6th-order moment"),c("X","Y"))
print(round(res,digits=digits), quote = FALSE, right = TRUE,...)
invisible(object)
}
mean.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(mean(x,na.rm = TRUE,...),mean(y,na.rm = TRUE,...)))
}
cv.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(cv(x,...),cv(y,...)))
}
max.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(max(x,na.rm = TRUE),max(y,na.rm = TRUE)))
}
min.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(min(x,na.rm = TRUE),min(y,na.rm = TRUE)))
}
skewness.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(skewness(x),skewness(y)))
}
kurtosis.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(kurtosis(x),kurtosis(y)))
}
Median.bridgesde2d <- function(x, at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(Median(x),Median(y)))
}
Mode.bridgesde2d <- function(x, at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(Mode(x),Mode(y)))
}
quantile.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(list(x=quantile(x,...),y=quantile(y,...)))
}
moment.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(c(moment(x,...),moment(y,...)))
}
bconfint.bridgesde2d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde2d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)}else{
X = object$X
Y = object$Y}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
return(list(x=bconfint(x,...),y=bconfint(y,...)))
}
time.bridgesde2d <- function(x,...)
{
class(x) <- "bridgesde2d"
as.vector(time(x$X))
}
#####
##### bridgesde3d
bridgesde3d <- function(N, ...) UseMethod("bridgesde3d")
bridgesde3d.default <- function(N =1000,M=1,x0=c(0,0,0),y=c(0,0,0),t0=0,T=1,Dt,drift,diffusion,
corr=NULL, alpha=0.5,mu=0.5,type=c("ito","str"), method=c(
"euler","milstein","predcorr","smilstein","taylor",
"heun","rk1","rk2","rk3"),...)
{
if (!is.numeric(x0) || length(x0) != 3)
stop("'x0' must be numeric, and length(x0) = 3 ")
if (!is.numeric(y) || length(y) != 3)
stop("'y' must be numeric, and length(y) = 3 ")
if (any(!is.numeric(t0) || !is.numeric(T)))
stop(" 't0' and 'T' must be numeric")
if (any(!is.numeric(N) || (N - floor(N) > 0) || N <= 1))
stop(" 'N' must be a positive integer ")
if (any(!is.numeric(M) || (M - floor(M) > 0) || M <= 0))
stop(" 'M' must be a positive integer ")
if (length(drift) !=3 )
stop("drift must be expression 3d (vector of 3 expression)")
if (length(diffusion) !=3 )
stop("diffusion must be expression 3d (vector of 3 expression)")
if (any(!is.expression(drift) || !is.expression(diffusion) ))
stop(" coefficient of 'drift' and 'diffusion' must be expressions in 't', 'x', 'y' and 'z'")
if (missing(type)) type <- "ito"
if (!is.null(corr)) {
if (any(!is.matrix(corr) || det(corr) <= 0 || nrow(corr)!=3 || ncol(corr)!=3 || !isSymmetric(corr)))
stop("the correlation structure of W1(t), W2(t) and W3(t) must be a real symmetric\n positive-definite square matrix of dimension 3")
}
method <- match.arg(method)
if (!is.null(corr) && method != "euler" && method != "milstein") {method="euler"}
if (method =="predcorr"){
if (any(alpha > 1 || alpha < 0))
stop("please use '0 <= alpha <= 1' ")
if (any(mu > 1 || mu < 0))
stop("please use '0 <= mu <= 1' ")
}
if ( t0 < 0 || T < 0 )
stop(" please use positive times! (0 <= t0 < T) ")
if (missing(Dt)) {
t <- seq(t0, T, by=Dt)
} else {
t <- c(t0, t0 + cumsum(rep(Dt, N)))
T <- t[N + 1]
}
Dt <- (T - t0)/N
X1 <- snssde3d(N,M,x0=c(x0[1],x0[2],x0[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$X
Y1 <- snssde3d(N,M,x0=c(x0[1],x0[2],x0[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Y
Z1 <- snssde3d(N,M,x0=c(x0[1],x0[2],x0[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Z
if (M > 1){X2 <- apply(data.frame(snssde3d(N,M,x0=c(y[1],y[2],y[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$X),2,rev)
Y2 <- apply(data.frame(snssde3d(N,M,x0=c(y[1],y[2],y[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Y),2,rev)
Z2 <- apply(data.frame(snssde3d(N,M,x0=c(y[1],y[2],y[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Z),2,rev)
}else{X2 <- rev(snssde3d(N,M,x0=c(y[1],y[2],y[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$X)
Y2 <- rev(snssde3d(N,M,x0=c(y[1],y[2],y[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Y)
Z2 <- rev(snssde3d(N,M,x0=c(y[1],y[2],y[3]),t0,T,Dt,drift,diffusion,corr,alpha,mu,type, method,...)$Z)
}
Gx = Gy = Gz <- rep(NA,M)
if (M > 1){
for(j in 1:M){
if (X1[1,M] >= X2[1,M]){
if (!all(X1[,j] > X2[,j]))
Gx[j] <- min(which((X1[,j]-X2[,j]) <= 0)) - 1
}else{ if (!all(X1[,j] < X2[,j]))
Gx[j] <- min(which((X1[,j]-X2[,j]) >= 0)) - 1
}
}
}else{
if (X1[1] >= X2[1]){
if (!all(X1 > X2))
Gx <- min(which((X1-X2) <= 0)) - 1
}else{ if (!all(X1 < X2))
Gx <- min(which((X1-X2) >= 0)) - 1
}
}
if (M > 1){
for(j in 1:M){
if (Y1[1,M] >= Y2[1,M]){
if (!all(Y1[,j] > Y2[,j]))
Gy[j] <- min(which((Y1[,j]-Y2[,j]) <= 0)) - 1
}else{ if (!all(Y1[,j] < Y2[,j]))
Gy[j] <- min(which((Y1[,j]-Y2[,j]) >= 0)) - 1
}
}
}else{
if (Y1[1] >= Y2[1]){
if (!all(Y1 > Y2))
Gy <- min(which((Y1-Y2) <= 0)) - 1
}else{ if (!all(Y1 < Y2))
Gy <- min(which((Y1-Y2) >= 0)) - 1
}
}
if (M > 1){
for(j in 1:M){
if (Z1[1,M] >= Z2[1,M]){
if (!all(Z1[,j] > Z2[,j]))
Gz[j] <- min(which((Z1[,j]-Z2[,j]) <= 0)) - 1
}else{ if (!all(Z1[,j] < Z2[,j]))
Gz[j] <- min(which((Z1[,j]-Z2[,j]) >= 0)) - 1
}
}
}else{
if (Z1[1] >= Z2[1]){
if (!all(Z1 > Z2))
Gz <- min(which((Z1-Z2) <= 0)) - 1
}else{ if (!all(Z1 < Z2))
Gz <- min(which((Z1-Z2) >= 0)) - 1
}
}
Gx[which(Gx==0)]=NA
Gy[which(Gy==0)]=NA
Gz[which(Gz==0)]=NA
G <- na.omit(data.frame(Gx,Gy,Gz))
Gx <- G$Gx;Gy <- G$Gy; Gz <- G$Gz
namex <- "X"
namey <- "Y"
namey <- "Z"
namex <- if(M > 1) paste("X",1:length(which(!is.na(Gx))),sep="")
namey <- if(M > 1) paste("Y",1:length(which(!is.na(Gy))),sep="")
namez <- if(M > 1) paste("Z",1:length(which(!is.na(Gz))),sep="")
if (M == 1){ if (is.na(Gx) ){stop( "A crossing has been no realized,trying again (Repeat)..." )}else{X <- ts(c(X1[1:Gx],X2[-(1:Gx)]),start=t0,deltat=deltat(X1),names=namex)}
}else if (M > 1){ if(length(which(is.na(Gx))) == M){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else if (length(which(!is.na(Gx))) == 1){X <- ts(c(X1[,-which(is.na(Gx))][1:Gx[which(!is.na(Gx))]],X2[,-which(is.na(Gx))][-(1:Gx[which(!is.na(Gx))])]),start=t0,deltat=deltat(X1),names=namex)
}else if (length(which( is.na(Gx))) == 0){X <- ts(sapply(1:length(which(!is.na(Gx))), function(j) c(X1[,j][1:Gx[!is.na(Gx)][j]],X2[,j][-(1:Gx[!is.na(Gx)][j])])),start=t0,deltat=deltat(X1),names=namex)
}else{ X1 <- X1[,-c(which(is.na(Gx)))]
X2 <- X2[,-c(which(is.na(Gx)))]
X <- ts(sapply(1:length(which(!is.na(Gx))), function(j) c(X1[,j][1:Gx[!is.na(Gx)][j]],X2[,j][-(1:Gx[!is.na(Gx)][j])])),
start=t0,deltat=deltat(X1),names=namex)
}
}
if (M == 1){ if (is.na(Gy) ){stop( "A crossing has been no realized,trying again (Repeat)..." )}else{Y <- ts(c(Y1[1:Gy],Y2[-(1:Gy)]),start=t0,deltat=deltat(Y1),names=namey)}
}else if (M > 1){ if(length(which(is.na(Gy))) == M){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else if (length(which(!is.na(Gy))) == 1){Y <- ts(c(Y1[,-which(is.na(Gy))][1:Gy[which(!is.na(Gy))]],Y2[,-which(is.na(Gy))][-(1:Gy[which(!is.na(Gy))])]),start=t0,deltat=deltat(Y1),names=namey)
}else if (length(which( is.na(Gy))) == 0){Y <- ts(sapply(1:length(which(!is.na(Gy))), function(j) c(Y1[,j][1:Gy[!is.na(Gy)][j]],Y2[,j][-(1:Gy[!is.na(Gy)][j])])),start=t0,deltat=deltat(Y1),names=namey)
}else{ Y1 <- Y1[,-c(which(is.na(Gy)))]
Y2 <- Y2[,-c(which(is.na(Gy)))]
Y <- ts(sapply(1:length(which(!is.na(Gy))), function(j) c(Y1[,j][1:Gy[!is.na(Gy)][j]],Y2[,j][-(1:Gy[!is.na(Gy)][j])])),
start=t0,deltat=deltat(Y1),names=namey)
}
}
if (M == 1){ if (is.na(Gz) ){stop( "A crossing has been no realized,trying again (Repeat)..." )}else{Z <- ts(c(Z1[1:Gz],Z2[-(1:Gz)]),start=t0,deltat=deltat(Z1),names=namez)}
}else if (M > 1){ if(length(which(is.na(Gz))) == M){stop( "A crossing has been no realized,trying again (Repeat)..." )
}else if (length(which(!is.na(Gz))) == 1){Z <- ts(c(Z1[,-which(is.na(Gz))][1:Gz[which(!is.na(Gz))]],Z2[,-which(is.na(Gz))][-(1:Gz[which(!is.na(Gz))])]),start=t0,deltat=deltat(Z1),names=namez)
}else if (length(which( is.na(Gz))) == 0){Z <- ts(sapply(1:length(which(!is.na(Gz))), function(j) c(Z1[,j][1:Gz[!is.na(Gz)][j]],Z2[,j][-(1:Gz[!is.na(Gz)][j])])),start=t0,deltat=deltat(Z1),names=namez)
}else{ Z1 <- Z1[,-c(which(is.na(Gz)))]
Z2 <- Z2[,-c(which(is.na(Gz)))]
Z <- ts(sapply(1:length(which(!is.na(Gz))), function(j) c(Z1[,j][1:Gz[!is.na(Gz)][j]],Z2[,j][-(1:Gz[!is.na(Gz)][j])])),
start=t0,deltat=deltat(Z1),names=namez)
}
}
structure(list(X=X,Y=Y,Z=Z, driftx=drift[[1]], diffx=diffusion[[1]],drifty=drift[[2]], diffy=diffusion[[2]],driftz=drift[[3]],diffz=diffusion[[3]],
type=type,corrmat=corr,method=method, x0=x0,y=y, N=N,M=M,Dt=Dt,t0=t0,T=T,Cx=Gx,Cy=Gy,Cz=Gz,dim="3d",call=match.call()),class="bridgesde3d")
}
###
print.bridgesde3d <- function(x, digits=NULL, ...)
{
class(x) <- "bridgesde3d"
Ito = "It\xf4"
Encoding(Ito) <- "latin1"
if (x$method=="euler")
sch <- "Euler scheme with order 0.5"
else if (x$method=="milstein")
sch <- "First-order Milstein scheme"
else if (x$method=="predcorr")
sch <- "Predictor-corrector method with order 1"
else if (x$method=="smilstein")
sch <- "Second-order Milstein scheme"
else if (x$method=="taylor")
sch <- "Taylor scheme with order 1.5"
else if (x$method=="heun")
sch <- "Heun scheme with order 2"
else if (x$method=="rk1")
sch <- "Runge-Kutta method with order 1"
else if (x$method=="rk2")
sch <- "Runge-Kutta method with order 2"
else if (x$method=="rk3")
sch <- "Runge-Kutta method with order 3"
Drx <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)),z=quote(Z(t)))), list(e = x$driftx))))
DDx <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)),z=quote(Z(t)))), list(e = x$diffx))))
Dry <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)),z=quote(Z(t)))), list(e = x$drifty))))
DDy <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)),z=quote(Z(t)))), list(e = x$diffy))))
Drz <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)),z=quote(Z(t)))), list(e = x$driftz))))
DDz <- deparse(eval(substitute(substitute(e, list(x=quote(X(t)),y=quote(Y(t)),z=quote(Z(t)))), list(e = x$diffz))))
if (is.null(x$corrmat)){
if(x$type=="ito"){
cat(Ito," Bridge Sde 3D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," * dW1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," * dW2(t)","\n",
"\t| dZ(t) = ", Drz," * dt + ", DDz," * dW3(t)","\n",
"Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),",",format(x$x0[3],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),",",format(x$y[3],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}else{
cat("Stratonovich Bridge Sde 3D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," o dW1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," o dW2(t)","\n",
"\t| dZ(t) = ", Drz," * dt + ", DDz," o dW3(t)","\n",
"Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),",",format(x$x0[3],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),",",format(x$y[3],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}} else {
if(x$type=="ito"){
cat(Ito," Bridge Sde 3D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," * dB1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," * dB2(t)","\n",
"\t| dZ(t) = ", Drz," * dt + ", DDz," * dB3(t)","\n",
"\t| Correlation structure:",sep="")
prmatrix(x$corrmat,rowlab = rep(" ", 3), collab = rep("", 3),digits=digits)
cat("Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),",",format(x$x0[3],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),",",format(x$y[3],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}else{
cat("Stratonovich Bridge Sde 3D:","\n",
"\t| dX(t) = ", Drx," * dt + ", DDx," o dB1(t)","\n",
"\t| dY(t) = ", Dry," * dt + ", DDy," o dB2(t)","\n",
"\t| dZ(t) = ", Drz," * dt + ", DDz," o dB3(t)","\n",
"\t| Correlation structure:",sep="")
prmatrix(x$corrmat,rowlab = rep(" ", 3), collab = rep("", 3),digits=digits)
cat("Method:","\n",
"\t| ",sch,"\n",
"Summary:","\n",
"\t| Size of process","\t| N = ",format(x$N+1,digits=digits),".","\n",
"\t| Crossing realized","\t| C = ",format(length(which(!is.na(x$Cx))),digits=digits)," among ",format(x$M,digits=digits),".","\n",
"\t| Initial values","\t| x0 = ","(",format(x$x0[1],digits=digits),",",format(x$x0[2],digits=digits),",",format(x$x0[3],digits=digits),")",".","\n",
"\t| Ending values","\t\t| y = ","(",format(x$y[1],digits=digits),",",format(x$y[2],digits=digits),",",format(x$y[3],digits=digits),")",".","\n",
"\t| Time of process","\t| t in [",format(x$t0,digits=digits),",",format(x$T,digits=digits),"].","\n",
"\t| Discretization","\t| Dt = ",format(x$Dt,digits=digits),".","\n",
sep="")}
}
invisible(x)
}
##
## summary
summary.bridgesde3d <- function(object,at,digits=NULL,...)
{
class(object) <- "bridgesde3d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (is.null(digits)){digits = base::options()$digits}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)
Z = matrix(object$Z,nrow=length(object$Z),ncol=1)}else{
X = object$X
Y = object$Y
Z = object$Z}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
z <- as.vector(Z[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
if (length(z) == 0){
if (as.numeric(C)==1){ Fz <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Z))}else{
Fz <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Z[,i]))}
z <- sapply(1:length(Fz),function(i) Fz[[i]](at))
}
cat("\nMonte-Carlo Statistics for (X(t),Y(t),Z(t)) at time t = ",at,"\n",
"\t| Crossing realized ",C," among ",object$M,"\n",
sep="")
res <- as.data.frame(matrix(c(mean(x,na.rm = TRUE),var(x,na.rm = TRUE),Median(x),Mode(x),
quantile(x,0.25,na.rm = TRUE),quantile(x,0.75,na.rm = TRUE),min(x,na.rm=TRUE),max(x,na.rm=TRUE),
skewness(x),kurtosis(x),cv(x),moment(x,order=3,center=TRUE),
moment(x,order=4,center=TRUE),moment(x,order=5,center=TRUE),moment(x,order=6,center=TRUE),
mean(y,na.rm = TRUE),var(y,na.rm = TRUE),Median(y),Mode(y),
quantile(y,0.25,na.rm = TRUE),quantile(y,0.75,na.rm = TRUE),min(y,na.rm=TRUE),max(y,na.rm=TRUE),
skewness(y),kurtosis(y),cv(y),moment(y,order=3,center=TRUE),
moment(y,order=4,center=TRUE),moment(y,order=5,center=TRUE),moment(y,order=6,center=TRUE),
mean(z,na.rm = TRUE),var(z,na.rm = TRUE),Median(z),Mode(z),
quantile(z,0.25,na.rm = TRUE),quantile(z,0.75,na.rm = TRUE),min(z,na.rm=TRUE),max(z,na.rm=TRUE),
skewness(z),kurtosis(z),cv(z),moment(z,order=3,center=TRUE),
moment(z,order=4,center=TRUE),moment(z,order=5,center=TRUE),moment(z,order=6,center=TRUE)),
ncol=3))
dimnames(res) <- list(c("Mean","Variance","Median","Mode","First quartile","Third quartile","Minimum","Maximum","Skewness","Kurtosis",
"Coef-variation","3th-order moment","4th-order moment","5th-order moment","6th-order moment"),c("X","Y","Z"))
print(round(res,digits=digits), quote = FALSE, right = TRUE,...)
invisible(object)
}
mean.bridgesde3d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde3d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)
Z = matrix(object$Z,nrow=length(object$Z),ncol=1)}else{
X = object$X
Y = object$Y
Z = object$Z}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
z <- as.vector(Z[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
if (length(z) == 0){
if (as.numeric(C)==1){ Fz <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Z))}else{
Fz <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Z[,i]))}
z <- sapply(1:length(Fz),function(i) Fz[[i]](at))
}
return(c(mean(x,na.rm = TRUE,...),mean(y,na.rm = TRUE,...),mean(z,na.rm = TRUE,...)))
}
cv.bridgesde3d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde3d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)
Z = matrix(object$Z,nrow=length(object$Z),ncol=1)}else{
X = object$X
Y = object$Y
Z = object$Z}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
z <- as.vector(Z[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X[,i]))}
x <- sapply(1:length(Fx),function(i) Fx[[i]](at))
}
if (length(y) == 0){
if (as.numeric(C)==1){ Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y))}else{
Fy <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Y[,i]))}
y <- sapply(1:length(Fy),function(i) Fy[[i]](at))
}
if (length(z) == 0){
if (as.numeric(C)==1){ Fz <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Z))}else{
Fz <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$Z[,i]))}
z <- sapply(1:length(Fz),function(i) Fz[[i]](at))
}
return(c(cv(x,...),cv(y,...),cv(z,...)))
}
max.bridgesde3d <- function(x,at,...)
{
object <- x
class(object) <- "bridgesde3d"
if (missing(at)) {at = (as.numeric(object$T)-as.numeric(object$t0))/2}
if (any(object$T < at || object$t0 > at) )
stop( " please use 't0 <= at <= T'")
C = length(which(!is.na(object$Cx)))
if (as.numeric(C) == 1){ X = matrix(object$X,nrow=length(object$X),ncol=1)
Y = matrix(object$Y,nrow=length(object$Y),ncol=1)
Z = matrix(object$Z,nrow=length(object$Z),ncol=1)}else{
X = object$X
Y = object$Y
Z = object$Z}
x <- as.vector(X[which(time(object)==as.character(at)),])
y <- as.vector(Y[which(time(object)==as.character(at)),])
z <- as.vector(Z[which(time(object)==as.character(at)),])
if (length(x) == 0){
if (as.numeric(C)==1){ Fx <- lapply(1:as.numeric(C),function(i) approxfun(time(object),object$X))}else{
Fx <- lapply(1: