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

Question

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}$

Answer 1

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.

Answer 2

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

Answer 3

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.