Simulating Bounded Jacobi Volatility with the Milstein Scheme

# monte carlo simulation


# monte carlo simulation












I'm a graduate student working on a real options research problem suggested to me by my advisor. I'm not looking for a solution, but I'd like to know about the feasibility of numerically solving the following as I have not made progress for a very long time.

In the problem, I use a Jacobi process to model volatility. The Jacobi process follows the dynamics: $$dV(t)=\kappa(\theta-V(t))dt+ {\sigma}_v \sqrt{\frac{(V(t)-v_{min})(v_{max}-V(t))}{(\sqrt{v_{max}}-\sqrt{v_{min}})^2}}dW^V(t)$$ where $W^V(t)$ is a standard Brownian motion.

Let $C_1,C_2,C_3$ be constants. I have reduced my problem to evaluating $$E[e^{\int_t^sV(p)(C_1(e^{-a(s-p)}-1)-C_2)dp}|V(t)=v] $$ where $s>t$. I would like to find a closed form solution for this, but I don't know how to proceed in that direction. The goal is to calculate $$L_j(t)=C_3\int_t^{j+b}[\int_t^qE[e^{\int_t^sV(p)(C_1(e^{-a(s-p)}-1)-C_2)dp}|V(t)=v]ds]e^{-\frac{q}{\lambda}}dq$$ The trouble is that I am supposed to calculate this for all integers $j$ from $1$ to $20,000$ and for all $t$ from $j$ to $j+120$. If I had to calculate $L_j(t)$ one time, monte carlo simulation would be fine, but I will have to run a large number of monte carlo simulations.

My initial thought was that if I were to do a monte carlo simulation, all these conditional expectations will lead to a branching sequence where at each point in time I'm generating more and more paths. If the total time values I need to do this for is $t=20,000$, I'm not sure this is realistic computationally.

This is what I have tried in regards to calculating the expectation:

- I applied Ito's formula to the exponential, integrated, took expectations, and converted what I had to an ODE. The plan was to solve the ODE, but I got stuck here.

- I considered expanding the exponential in a Taylor series, but the integral inside the exponential is not bounded by 1 unless I pick very unfavorable parameters. The parameters arise in the constants $C_1,C_2$.

- I applied Feynman-Kac. The resulting PDE is: $$\frac{\partial u(v,t)}{\partial t}+\kappa(\theta-v)\frac{\partial u(v,t)}{\partial v}+\frac{\sigma_v^2}{2}\frac{(v-v_{min})(v_{max}-v)}{(\sqrt{v_{max}}-\sqrt{v_{min}})^2}\frac{\partial^2 u(v,t)}{\partial v^2}+v[C_2-C_1(e^{-a(s-t)}-1)]=0$$ with $u(v,s)=1$. I don't know enough PDE's to know if I can solve this analytically.

- I have tried Monte Carlo simulation, but this seems very computationally expensive. The paths of the Jacobi process are bounded. So, $v_{min}\leq V(t)\leq v_{max}$. As $V(t)$ approaches $v_{min}$ or $v_{max}$, the drift term begins taking over. Maybe at each point in time, I can partition $[v_{min},v_{max}]$ and average the conditional expectations given $V(t)=v(i)$ where $v(i)$ is a partition point.

Questions:

If I need to resort to monte carlo simulation, does this look feasible? The Jacobi process is stationary with the stationary distribution being a beta distribution. Maybe I can use this fact to make the simulation more efficient. There are formulas involving infinite summations for the conditional expectations of the form $E[V(s)|V(t)=v]$. In calculating the conditional expectations, it seems like at each point in time, I need to generate increasingly more paths to the point that it becomes infeasible. Is this correct?

## Answer by alexbougias (score 1)

https://quant.stackexchange.com/a/54410

For academic purposes, you can resort to a MC simulation. If you are using the Euler scheme for the discretization of the SDE, then it is quite reasonable to expect a pathological behaviour near the bounds $v_{min}$ and $v_{max}$ in your MC simulation. I have also encountered this problem when I had to simulate paths from the Jacobi diffusion. What solved the problem for me was using the Milstein scheme (since the Jacobi has state-dependent drift and volatility parameters). All the simulated values $\hat{V}_{i\Delta t,j}$ were inside the bounds even for $N_{paths}>100000$ and $N_{steps}>5000$.

If you are using Matlab, you can make your algorithm less time-consuming using all the available workers of the cores through the parfor loop.

Regarding the Milstein scheme, you can read the following paper:

Higham, Desmond J. "An algorithmic introduction to numerical simulation of stochastic differential equations." SIAM review 43.3 (2001): 525-546.