Monte Carlo for MultiFactor Ornstein Uhlenbeck

Question

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 1

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 Matlab