How to simulate the average trajectories of the SDE variables?


I have a code to solve to SDEs I am working on, in R. The equations in question are, dX = (-X*a + (Y-X)b + h)dt + g dW and dY = (-Ya + (X-Y)*b)dt for time [0,200], a=0.30, b=0.2, g=1 and h is 1 for time [50,70] and 0 otherwise.

library(diffeqr) ##for solving ODEs SDEs DDEs..
library(Julia)  ##used  for the sde 

tryCatch(system2(file.path(JULIA_HOME, "julia"),
                 "-E \"try println(JULIA_HOME) catch e println(Sys.BINDIR) end;\"", stdout = TRUE)[1],
         warning = function(war) {},
         error = function(err) NULL)

julia_setup(JULIA_HOME ="/home/tim/.local/share/R/JuliaCall/julia/v1.5.0/bin")

##calling Julia
JuliaCall:::julia_locate(JULIA_HOME = "/home/tim/.local/share/R/JuliaCall/julia/v1.5.0/bin")

## setup of the packages
de <- diffeqr::diffeq_setup() 
julia <- julia_setup()

#u0 is the starting value for the vector
#func1 is the ODE TO SOLVE
#f is the drift part in the sde
#g is the diffusion part in the sde
#p are the parameters associated to the sde
# noise_rate_prototype is a sparse matrix by making a dense matrix and setting some values as not zero

#starting value for the three variables

## parameters

## parameters of the sde 
p<-c(a, b, g, h)

## drift part in the SDE
f <- JuliaCall::julia_eval("
function f(dy,y,p,t)
  MU = signal(t)                                    # <-
  dy[1] = -y[1]*p[1] + (y[2]-y[1])*p[2] + MU         # <- 
  dy[2] = -y[2]*p[1] + (y[1]-y[2])*p[2]

signal <- approxfun(x = c(0, 35, 125, 200),          # <- 
                    y = c(0, 1,  0,  0), 
                    method = "constant", rule = 2)

## diffusion part in the SDE
g <- JuliaCall::julia_eval("
function g(dy,y,p,t)
  dy[1,1] = p[3]
  dy[2,1] = 0

## matrix associated to the noise
noise_rate_prototype <- matrix(rep(1,2), nrow = 2, ncol = 1)

tspan <- c(0,200)

JuliaCall::julia_assign("u0", u0); ## initial value
JuliaCall::julia_assign("tspan", tspan); ## grid
JuliaCall::julia_assign("p", p); ## parameters
JuliaCall::julia_assign("signal", signal);                   # <-
JuliaCall::julia_assign("noise_rate_prototype", noise_rate_prototype); ## matrix associated to the noise

## call the SDE problem
prob <- JuliaCall::julia_eval("SDEProblem(f, g, u0,tspan,p, noise_rate_prototype=noise_rate_prototype, signal=signal)");      # <-
sol <- de$solve(prob)

udf <-$u,identity)))
#plot the 2 figures, 
plot(sol$t,udf$V1,type="l",col="purple", ylab="Y1") #L
lines(sol$t,udf$V2,col="red",ylab="Y2") #R

I want to plot the average behavior of V1 and V2 over the course of, say 100 trials. I have tried using a for loop to just run this code a 100 times and stack the results in a data frame side by side, but I get an error because not all of the udf$V1 will have the same length in the loop. How do I go about it?

I can not speak to the particulars of calling from R, but generally if you want average properties of the solution of an SDE, you wrap the problem in an Ensemble Problem and then you can use the EnsembleSummary. Hopefully that helps?

1 Like

And EnsembleProblem can be used from R as well. It’s showcased in the GPU acceleration examples.