Unexpected output from `OU`

[lioumens]

OU returns unexpected sample paths; I am getting sample paths that don't follow the theoretical properties of conditional variance.

library(Sim.DiffProc)
#> Package 'Sim.DiffProc', version 4.8
#> browseVignettes('Sim.DiffProc') for more informations.

# parameters
x0 <- 5
t0 <- 0
mu <- 3
sigma = .5
N <- 100
M <- 2000

set.seed(1)
ou_ts <- OU(N = N, M = M, x0 = x0, t0 = t0, mu = 3, sigma = sigma)

x <- seq(-2, 2, .01)
theoretical_sd <- sqrt(sigma^2 / 2/mu *( 1- exp(-2 * mu))) # at T = 1
theoretical_mean <- x0 * exp(-mu) # at T = 1
theory_dens <- dnorm(x, theoretical_mean, theoretical_sd) # should be normally distributed
hist(ou_ts[N+1,], freq = F, breaks = 30)
lines(x, theory_dens, col = "red")

^{Created on 2022-11-16 by the reprex package (v2.0.1)}

I believe the coded solution path in HWV.default for Ito.sum needs to have additional dependence on $t$ since we have

$x(t) = x_0e^{-\mu t} + \sigma \int_0^t e^{\mu (s - t)} dW(s)$

    hwv <- function() {
        w = c(0, cumsum(rnorm(N, mean = 0, sd = sqrt(Dt))))
        dw <- diff(w)
        X <- numeric()
        X[1] <- x0
        Ito.sum <- c(0, sapply(1:(N + 1), function(i) {
            exp(-mu * Dt) * dw[i]
        }))
        X <- sapply(1:(N + 1), function(i) {
            theta + (x0 - theta) * exp(-mu * t[i]) + sigma * 
                sum(Ito.sum[1:i])
        })
        X
    }