I am trying to simulate the following OU process.
$$ dX_t = -\theta X_t dt + \sigma dW_t $$
I then simulated the process using the Eulers-Maruyama discretization approach i.e
$$ X_{t+\Delta t} = X_{t} -\theta X_t \Delta t + \sigma \sqrt{\Delta t}\epsilon $$
where $\epsilon \sim N(0,1)$
Because I want to simulate hourly data and my parameters have already been calibrated at that same frequency I've assumed that $\Delta t = 1$
I've been doing some tests to ensure that the variance of the simulations matches the closed-form solution for the variance and that it also converges to its long-term variance $$\frac{\sigma^2}{2\theta}(1-e^{-2\theta (T-t)})$$
What I've noticed is that the variance tends to be overestimated the larger the value of $\Delta t$. I think this is due to the fact that compared to the actual solution for $X_T = X_te^{-\theta (T-t)} + \int^T_t\sigma e^{-\theta(T-s)} dW_s$ the Euler-Maruyama is a first-order approximation.
It is easy to see that if we do a first-order approximation of the first term after the equal sign
$$ X_te^{-\theta \Delta t} \approx X_t(1-\theta \Delta t) $$
Something that I've noticed is that the difference between the variance of the simulations when I use $\Delta t = 1$ and the true variance of the process is proportional. Heuristically I've found that I can get the simulated variance to match the true variance from the process if I multiply the standard deviation of the Euler-Maruyama process by the scalar
$$e^{-\frac{1}{4} \theta \Delta t}$$
and
$$ X_{t+\Delta t} = X_{t} -\theta X_t \Delta t + \sigma e^{-\frac{1}{4} \theta \Delta t} \sqrt{\Delta t}\epsilon $$
I have used different values of $\sigma$ and $\theta$ to test that it is the adjustment holds when $\Delta t = 1$ however I've been unable to prove mathematically why this is the correct adjustment term.
Has anyone come across this before? I would love to find a way to explain this formally