Exact Simulation of a Multivariate Ornstein–Uhlenbeck Process
Summary
The document asks how to simulate a multifactor Ornstein–Uhlenbeck process with matrix-valued mean reversion and diffusion. It contrasts the scalar exact time-step formula with the multivariate solution, whose random increment has correlated components. Diagonalizing the mean-reversion matrix can transform the system into scalar equations, but the transformed noise still requires its covariance structure to be handled correctly; scalar simulation does not by itself make all components independent.
The answer says direct simulation is possible if correlated increments are generated for the original state vector, while the transformed approach offers a way to work with decoupled dynamics. For checking an implementation, it recommends reproducing a known multivariate example and comparing against an exact solution or a suitable alternative scheme. The exchange does not give a complete algorithm, covariance factorization, or benchmark values, and the displayed equations may need careful scrutiny before implementation. It is guidance on the modeling issue and validation approach, rather than a turnkey simulation specification.
Key ideas
- The scalar Ornstein–Uhlenbeck update cannot be applied independently to correlated multivariate state variables without accounting for their joint noise.
- Direct simulation can use correlated increments for the original multivariate process.
- Diagonalization can simplify the dynamics, while the transformed noise covariance still matters.
- A known multivariate example and an exact solution can serve as benchmarks for implementation.
Tags
Full text
# Monte Carlo for MultiFactor Ornstein Uhlenbeck
# Monte Carlo for MultiFactor Ornstein Uhlenbeck
I'm following loosely the exposition given in "Monte Carlo Methods in Financial Engineering by Glasserman.
For a multifactor OU process:
$dX(t)=C(b-X(t))dt+DdW(t)$
Where C and D are d*d matrices and b and X(t) are vectors on length d, and W is a d dimensional brownian motion.
He notes that this can be used to define an "exact" discretization similar to the 1D case (Shown below)
1D case:
$dr(t)=\alpha(b-r(t))dt+\sigma dW(t)$
which has solution
$r(t)=exp^{-\alpha(t-u)}r(u)+\alpha\int_u^t \exp^{- \alpha(t-s)}b(s) ds+ \sigma \int_u^t \exp^{- \alpha(t-s)}dW(s)$
and can be simulated as
$r(t+1)=exp^{-\alpha(t_{i+1}-t_i)}r(t_i)+\mu(t_i,t_{t+1})+\sigma_r(t_i,t_{i+1})Z_{i_1} $ -- EQ1
where $\mu(u,t)=\alpha \int_u^t \exp^{-\alpha(t-s)}b(s)ds$ if b is constant $\mu=b(1-\exp^{-\alpha(t_{i+1}-t_i)})$
and
$\sigma_r^2(u,t)=\sigma^2 \int_u^t \exp^{-2\alpha(t-s)}ds=\frac{\sigma^2}{2 \alpha}(1-exp^{-2 \alpha(t-u)})$
returning to multi-factor case
The general solution is
$X(t)=exp^{-C(t-u)}X(u)+\int_u^t \exp^{- C(t-s)}b ds+ \int_u^t \exp^{- C(t-s)}DdW(s)$ -- EQ2
However he also goes on to show that when C is diagonizable using $VCV^{-1}=\Delta$
$dY(t)=VX(t)$
which after simplification becomes
$dY(t)=\Delta(\tilde{b}-Y(t))dt+d \tilde W(t)$
with $\tilde{W}$ a $BM(0,\Sigma)$ where $\Sigma=VDD'T'$
This can be modeled as
$Y_j(t+1)=exp^{-\lambda_j (t_{i+1}-t_i)}Y_j(t_i)+(\exp^{\lambda_j(t_{i+1}-t_i)}-1)\tilde{b}_j +\sqrt{\frac{1}{2 \lambda_j}(1-exp^{(-2\lambda_j(t_{i+1}-t_i)})}\xi_j(i+1)$ EQ3
where $\xi(1), \xi(2)$ are independent $N(0,\Sigma)$. Which reduces to a system of scalar simulations
You can then recover the $X(t)$ using $X(t)=V(t)^{-1}Y(t)$
My question is
1.) Could I run the multifactor OU process without doing the diagonalisation, i.e using the general solution Eq2 and the discretization scheme shown in EQ1? (if not why not?)
I realize that the second solution utilizing the diagonalisation is prolly neater (and faster) but I would like to be able to simulate the process two ways so that I can test the models against each other.
2.) If the answer to question 1 is no that how can I benchmark my model to enusre it's correct.
## Answer by user12348 (score 4, accepted)
https://quant.stackexchange.com/a/11432
EQ1 is uni-variate case. EQ2 is multivariate case, in which you have to use correlated $X_t$. His way of doing is making $Y_t$ independent so that you can simulate freely. He does so by finding PC on $\Delta$. Alternatively, you could generate correlated $X_t$ in your simulation.
To benchmark your model / code, you should first test and reproduce a given example of multivariate OU and see if you match the solution. There are ways to use exact solution Multivariate exact solution , another by Gillespie, goes up to bi-variate case and here is Gillespie implemented in MatlabShown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.