# Simulate from a SDE where drift and diffusion terms are matrices using Yuima in R

I'm trying to implement an SDE in R using Yuima. Most of examples and literature show how to implement and how the math works for SDE where drift and diffusion terms are scalar. What if I want to implement an SDE like this: $$dX_t = (b + BX(t))dt + \Lambda dW(t)$$ where $$b$$ is a vector of length $$3$$, $$B$$ is a $$3\times 3$$ matrix, $$\Lambda$$ is $$3\times 18$$, $$W_t$$ is $$R^{18}$$-valued.

The simulation of this process should return a $$t\times 3$$ matrix wich columns are realizations of $$X_t = X_{t1}, X_{t2}, X_{t3}$$

• How is $\Lambda$ $3 \times 18$? Commented Jan 22, 2023 at 11:14
• First of all, I think a $d$ is missing from the beginning of your equation...If those coefficients are scalars, then you can simply solve the SDE and implement the solution (if any exists)...however if those are not constants, then there is the chanche that you can't even solve them explicit. In this case I would just simply discretize the SDE: calculate $\Delta X_{t}$ with the discretized $\Delta X_{t}=\left(b+BX_{t}\right)\Delta t+\Lambda\Delta W_{t}$ SDE. Use small $\Delta t=t_{i+1}-t_{i}$ differences and simulate $\Delta X_{t}=X_{t_{i+1}}-X_{t_{i}}$ with a for cycle ... Commented Jan 22, 2023 at 12:56
• ... Here $\Delta W_{t}=W_{t_{i+1}}-W_{t_{i}}$ is a vector with correlated normal variables with $\mu=0$ and $\sigma^{2}=\Delta t$ parameters. How to correlate these variables is an other topic (see e.g.: Cholesky decomposition) ... about discretized SDE you should check Euler-Maruyama mthod (basically this is almost what I've described above without some material details) or Milstein method... Commented Jan 22, 2023 at 13:01
• @KapesMate thank you, i’ll give a try with a for loop. Do you see any problem with the dimensions of the matrices? I’m implementing a paper and I’m not sure about the derivation of b,B and Lambda
– Nic
Commented Jan 22, 2023 at 22:24
• @BobJansen do you think there is something wrong? I derived Lambda following a procedure described in a paper, and it looks correct
– Nic
Commented Jan 22, 2023 at 22:28

### Producing a working example:

We can produce a working example, following section 3.3 on p. 8 in the Yuima article. For the sake of simplicity, we reduce the dimension of $$\Lambda \in \mathbb{R}^{3\times 18}$$ to $$\Lambda \in \mathbb{R}^{3\times 4}$$. Then, we define the 3-dimensional Ito process with arbitrarily chosen values for $$b$$, $$B$$ and $$\Lambda$$ as:

\begin{align} X_t &= (b + B X_t) \: dt + \Lambda \: dW_t\\ &=\left(\begin{bmatrix} 0.5\\ 0.45\\ 0.4 \end{bmatrix} + \begin{bmatrix} 0.1 & 0.3 & 0.2\\ 0.3 & 0.2 & 0.3\\ 0.2 & 0.3 & 0.3\\ \end{bmatrix} \begin{bmatrix} X_{1t}\\ X_{2t}\\ X_{3t} \end{bmatrix} \right) \: dt + \begin{bmatrix} 0.1 & 0.3 & 0.2 & 0.4\\ 0.3 & 0.2 & 0.3 & 0.2\\ 0.2 & 0.3 & 0.3 & 0.5\\ \end{bmatrix} \begin{bmatrix} dW_{1t}\\ dW_{2t}\\ dW_{3t} \end{bmatrix}\\ &=\left(\begin{bmatrix} 0.5 + 0.1X_{1t}+0.3X_{2t}+0.2X_{3t}\\ 0.45+0.3X_{1t}+0.2X_{2t}+0.3X_{3t}\\ 0.4+0.2X_{1t}+0.3X_{2t}+0.3X_{3t} \end{bmatrix} \right) \: dt + \begin{bmatrix} 0.1 & 0.3 & 0.2 & 0.4\\ 0.3 & 0.2 & 0.3 & 0.2\\ 0.2 & 0.3 & 0.3 & 0.5\\ \end{bmatrix} \begin{bmatrix} dW_{1t}\\ dW_{2t}\\ dW_{3t} \end{bmatrix}.\\ \end{align}

The last equation follows the same matrix-form as provided in section 3.3. Therefore, we stringently follow the code setup from the example in the aforementioned section:

library(yuima)

sol <- c("x1", "x2", "x3")

a <- c("0.5+0.1*x1+0.3*x2+0.2*x3",
"0.45+0.3*x1+0.2*x2+0.3*x3",
"0.4+0.2*x1+0.3*x2+0.3*x3")

b <- matrix(c("0.1", "0.3", "0.2", "0.4",
"0.3", "0.2", "0.3", "0.2",
"0.2", "0.3", "0.3", "0.5"), 3, 4, byrow = T)

mod3 <- setModel(drift = a, diffusion = b, solve.variable = sol)

set.seed(123)
X <- yuima::simulate(mod3, xinit = 0.1)
yuima::plot(X, plot.type = "single", lty = 1:3)


The last code-snippet produces the following graph of the simulated 3-dimensional SDE:

You can further access the raw simulation data with X@[email protected]. With the example provided above, you should be able to extend it by your own means. I hope my answer provides some insight.

• this is completely ok, I'll add my contribution extending this example on cases when matrices are big and specifying the Lambda matrix element one by one is tedious.
– Nic
Commented Jan 24, 2023 at 10:53
• If you only need to simulate the SDE with a big $\Lambda$ once, then doing the labour intensive work might be the easiest solution. However, if you have multiple multivariate Ito processes you need to simulate, then I agree with you, it will become tedious to manually specify all the elements one by one. If you have $\Lambda$ in a numeric matrix and want an easy conversion, you can do class(lambda) <- "character" and all of the entries will become strings (just to help).
– Pleb
Commented Jan 24, 2023 at 11:23
• By doing it, I found that actually the most annoying part is specifying the drift, also because coefficients have generally lots of digits.
– Nic
Commented Jan 24, 2023 at 11:29
• @Nic Yes, I also realize that now. However, your own paste0 solution seems to do the trick, which is more than fine :-)
– Pleb
Commented Jan 24, 2023 at 11:35

If you have to deal with big matrices like $$\Lambda$$ or with complex drift like $$b + BX_t$$ and specifying drift and diffusion elements one-by-one is annoying this solution uses paste0 to create the character needed by drift parameters and a simple transformation from a numeric matrix to a character matrix of $$\Lambda$$ to create the diffusion parameter.

library(yuima)
n = 3
b = rnorm(n, 0, 1)
B = matrix(rnorm(n*n, 0, 1), n, n)
Lambda = matrix(rnorm(n*18, 0, 1), n, 18)

x_helper_X = paste0('*x', 1:n)
drift_character_X = character(n)

for(i in 1:n){
drift_character_X[i] = paste0(b[i], '+', paste0(paste0(B[i,], x_helper_X), collapse = '+'))
}

diffusion_character_X = apply(Lambda, 2, as.character)

sol = paste('x', 1:n)
drift_X = drift_character_X
diffusion_X = diffusion_character_X
mod_X = setModel(drift = drift_X, diffusion = diffusion_X, solve.variable = sol)

set.seed(123)
X_sim = simulate(mod_X)
plot(X_sim)


The result is the plot of the 3-dimension SDE

Here is my take of how you might implement the given SDE in R using base functions and Euler - Maruyama method :


# Define the drift and diffusion functions
b <- function(x) {c(3,2,1) + x}
Lambda <- function(x) {matrix(c(1,0,0,0,2,0,0,0,3), nrow = 3)}

# Set initial condition and simulation parameters
x0 <- c(0,0,0)
tmax <- 5
n <- 1000
dt <- tmax/n

# Initialize the solution matrix
X <- matrix(nrow = n, ncol = 3)
X[1,] <- x0

# Use the Euler-Maruyama method to simulate the solution
set.seed(123) # for reproducibility
for (i in 2:n) {
W_i <- rnorm(3)
X[i,] <- X[i-1,] + b(X[i-1,])*dt + Lambda(X[i-1,]) %*% W_i * sqrt(dt)
}



In this example, the variable X will be a matrix of size n rows and 3 columns, which are the realizations of $$X_{t1},X_{t2},X_{t3}.$$

Note that the Euler-Maruyama method is only one of the many methods that can be used to simulate a solution for a SDE, and it may not be the most accurate method for this particular SDE. In addition, it is usually a good idea to check the moments of the simulated solution to ensure that they match the moments of the true solution.

• It looks like a familiar style of answering!
– Nic
Commented Jan 22, 2023 at 22:57
• This doesn’t use the Yuima package Commented Jan 23, 2023 at 9:17