Plotting Several OU processes using ggplot and zoo

Question

I need assistance with adapting my ggplot and zoo code to produce plots similar to those generated by my base R code. I am working on a presentation about stochastic methods and need to plot several Ornstein-Uhlenbeck (OU) processes. These processes have various parameters influencing their behavior, and I want to display multiple OU walks on a single plot for comparison. While I can achieve this using base R, I am aiming to utilize ggplot and zoo as they result in cleaner plots.

Here's my base R code:

# Calculate the feature values according to an OU process
sim.ou <- function(theta, alpha, sigma, start_value) {
  x <- numeric(100)     # Time steps
  x[1] <- start_value   # Set the initial value
  for (i in 1:99) {
                        # OU simulation
    x[i + 1] <- x[i] - alpha * (x[i] - theta) + rnorm(n = 1, mean = 0, sd = sigma)
  }
  return(x)
}

# Plot the OU process
plotOU_diffusion <- function(sigma, theta, alpha, start_value, ...) {
  Xou <- replicate(1000, sim.ou(theta, alpha, sigma, start_value))
  matplot(Xou, type = "l", ylab = "", xlab = "", main = "", ...)
}

plotOU_diffusion(sigma=1,theta=1,alpha=1,start_value=1,col="#009E73",ylim=c(-7,7))
par(new=TRUE) 
plotOU_diffusion(sigma=0.5,theta=0.5,alpha=0.5,start_value=7,col="#56B4E9",ylim=c(-7,7))
par(new=TRUE)
plotOU_diffusion(sigma=0.1,theta=0.1,alpha=0.1,start_value=-4,col="#E69F00",ylim=c(-7,7))

This results in this plot:

And here's my ggplot and zoo code:

library(tidyverse)
library(zoo)

# Calculate the feature values according to an OU process
sim.ou <- function(theta, alpha, sigma, start_value) {
  x <- numeric(150)     # Time steps
  x[1] <- start_value   # Set the initial value
  for (i in 1:149) {
                        # OU simulation
    x[i + 1] <- x[i] - alpha * (x[i] - theta) + rnorm(n = 1, mean = 0, sd = sigma)
  }
  return(x)
}

# Plot the OU process
plotOU_diffusion <- function(sigma, theta, alpha, start_value, num_walks, col, tit = NULL, xlabs = NULL, ylabs = NULL) {
  # Create a matrix of OU values
  Xou <- replicate(num_walks, sim.ou(theta, alpha, sigma, start_value))
  # Plot Xou object
  autoplot(zoo(Xou), facet = NULL, alpha = 0.5) + 
    theme_classic() + 
    geom_hline(yintercept = theta, linetype = "longdash", color = "black") +
    theme(legend.position = "none") +
    labs(title = tit, x = xlabs, y = ylabs) +
    scale_color_manual(values = rep(col,num_walks)) +
    scale_y_continuous(limits = c(-5, 5), expand = c(0,0)) +
    scale_x_continuous(limits = c(0, 150), expand = c(0,0))
}

plotOU_diffusion(sigma=1,theta=1,alpha=1,start_value=1,num_walks=50,col="#009E73")
plotOU_diffusion(sigma=0.5,theta=0.5,alpha=0.5,start_value=7,num_walks=50,col="#56B4E9")
plotOU_diffusion(sigma=0.1,theta=0.1,alpha=0.1,start_value=-4,num_walks=50,col="#E69F00")

This results in the 3 following plots:

I'm seeking advice on how to modify the ggplot code to achieve a visual output similar to the base R plots. Any suggestions are greatly appreciated.

Answer 1

One option would be to "vectorize" your plotting function and instead of relying on autoplot uses ggplot(...) + geom_line() to draw the plot:

library(ggplot2)
library(zoo)

plotOU_diffusion <- function(
    sigma, theta, alpha, start_value, num_walks, col,
    tit = NULL, xlabs = NULL, ylabs = NULL) {
  Xou_list <- Map(
    function(sigma, theta, alpha, start_value, num_walks) {
      replicate(num_walks, sim.ou(theta, alpha, sigma, start_value))
    }, sigma, theta, alpha, start_value, num_walks
  )

  Xou_list <- Map(\(x, y) {
    df <- fortify(zoo(x), melt = TRUE)
    df$color <- y
    df
  }, Xou_list, seq_along(Xou_list))
  dat <- do.call(rbind, Xou_list)

  ggplot(dat, aes(Index, Value,
    color = factor(color),
    group = interaction(Series, color)
  )) +
    geom_line() +
    theme_classic() +
    geom_hline(yintercept = theta, linetype = "longdash", color = "black") +
    theme(legend.position = "none") +
    labs(title = tit, x = xlabs, y = ylabs) +
    scale_color_manual(values = col) +
    scale_y_continuous(limits = c(-5, 5), expand = c(0, 0)) +
    scale_x_continuous(limits = c(0, 150), expand = c(0, 0))
}

plotOU_diffusion(
  sigma = c(1, .5, .1),
  theta = c(1, .5, .1),
  alpha = c(1, .5, .1),
  start_value = c(1, 7, -4),
  num_walks = 50,
  col = c("#009E73", "#56B4E9", "#E69F00")
)
#> Warning: Removed 50 rows containing missing values (`geom_line()`).