title: How to solve a system of stochastic differential equations with time dependent parameters in R?
body: <p>I am trying to solve this system of SDEs in R, dX = (-X*a + (Y-X)<em>b + h)dt + g dW  and dY = (-Y</em>a + (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. I have the following code for the system for constant h=1:</p>
  <pre><code>library(diffeqr) ##for solving ODEs SDEs DDEs..
  library(deSolve)
  library(Julia)  ##used  for the sde 
  library(JuliaCall) 
  library(data.table)
  library(minpack.lm)
  
  
  tryCatch(system2(file.path(JULIA_HOME, &quot;julia&quot;),
                   &quot;-E \&quot;try println(JULIA_HOME) catch e println(Sys.BINDIR) end;\&quot;&quot;, stdout = TRUE)[1],
           warning = function(war) {},
           error = function(err) NULL)
  
  julia_setup(JULIA_HOME =&quot;/home/tim/.local/share/R/JuliaCall/julia/v1.5.0/bin&quot;)
  
  ##calling Julia
  JuliaCall:::julia_locate(JULIA_HOME = &quot;/home/tim/.local/share/R/JuliaCall/julia/v1.5.0/bin&quot;)
  
  
  ## setup of the packages
  de &lt;- diffeqr::diffeq_setup() 
  julia &lt;- 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
  u0&lt;-c(0,0)
  
  
  ## parameters
  a=0.33;
  b=0.2;
  g=1;
  h=1;
  
  ## parameters of the sde 
  p&lt;-c(a, b, g, h)
  
  ## drift part in the SDE
  f &lt;- JuliaCall::julia_eval(&quot;
  function f(dy,y,p,t)
    dy[1] = -y[1]*p[1] + (y[2]-y[1])*p[2] + p[4]
    dy[2] = -y[2]*p[1] + (y[1]-y[2])*p[2]
  end&quot;)
  
  
  
  ## diffusion part in the SDE
  g &lt;- JuliaCall::julia_eval(&quot;
  function g(dy,y,p,t)
    dy[1,1] = p[3]
    dy[2,1] = 0
  end&quot;)
  
  
  
  ## matrix associated to the noise
  noise_rate_prototype &lt;- matrix(rep(1,2), nrow = 2, ncol = 1)
  
  tspan &lt;- c(0.,200)
  
  JuliaCall::julia_assign(&quot;u0&quot;, u0); ## initial value
  JuliaCall::julia_assign(&quot;tspan&quot;, tspan); ## grid
  JuliaCall::julia_assign(&quot;p&quot;, p); ## parameters
  JuliaCall::julia_assign(&quot;noise_rate_prototype&quot;, noise_rate_prototype); ## matrix associated to the noise
  
  ## call the SDE problem
  prob &lt;- JuliaCall::julia_eval(&quot;SDEProblem(f, g, u0,tspan,p, noise_rate_prototype=noise_rate_prototype)&quot;);
  
  sol &lt;- de$solve(prob)
  
  udf &lt;- as.data.frame(t(sapply(sol$u,identity)))
      
  #plot the 2 figures, 
  plot(sol$t,udf$V1,type=&quot;l&quot;,col=&quot;purple&quot;, ylab=&quot;Y1&quot;) #L
  lines(sol$t,udf$V2,col=&quot;red&quot;,ylab=&quot;Y2&quot;) #R
  </code></pre>
  <p>I am not sure how to include the time dependence of <code>h</code>, I tried to do it by including a <a href="https://stackoverflow.com/questions/69843063/how-to-solve-a-system-of-ode-with-time-dependent-parameters-in-r/69846444#69846444">forcing variable</a>  but I got a lengthy Julia error (stating that <code>h</code> is not defined- code below) . So perhaps thats not the way to go. I'd appreciate any help regarding this.</p>
  <p>The chunk of the code I changed to include the time dependent parameters is this (the new or changes lines have an arrow to their right):</p>
  <pre><code>## drift part in the SDE
  f &lt;- JuliaCall::julia_eval(&quot;
  function f(dy,y,p,t)
    MU &lt;- signal(t)                                    # &lt;-
    dy[1] = -y[1]*p[1] + (y[2]-y[1])*p[2] + MU         # &lt;- 
    dy[2] = -y[2]*p[1] + (y[1]-y[2])*p[2]
  end&quot;)
  
  signal &lt;- approxfun(x = c(0, 35, 125, 200),          # &lt;- 
                      y = c(0, 1,  0,  0), 
                      method = &quot;constant&quot;, rule = 2)
  
  ## diffusion part in the SDE
  g &lt;- JuliaCall::julia_eval(&quot;
  function g(dy,y,p,t)
    dy[1,1] = p[3]
    dy[2,1] = 0
  end&quot;)
  
  
  
  ## matrix associated to the noise
  noise_rate_prototype &lt;- matrix(rep(1,2), nrow = 2, ncol = 1)
  
  tspan &lt;- c(0,200)
  
  JuliaCall::julia_assign(&quot;u0&quot;, u0); ## initial value
  JuliaCall::julia_assign(&quot;tspan&quot;, tspan); ## grid
  JuliaCall::julia_assign(&quot;p&quot;, p); ## parameters
  JuliaCall::julia_assign(&quot;signal&quot;, signal);                   # &lt;-
  JuliaCall::julia_assign(&quot;noise_rate_prototype&quot;, noise_rate_prototype); ## matrix associated to the noise
  
  ## call the SDE problem
  prob &lt;- JuliaCall::julia_eval(&quot;SDEProblem(f, g, u0,tspan,p, noise_rate_prototype=noise_rate_prototype, signal=signal)&quot;);      # &lt;-
  </code></pre>
  <p>The rest of the code remained the same.</p>