How to approximate stochastic volatility model with finite-state Markov chain?

Question

A common practice when computing solutions to stochastic dynamic optimization problems is to approximate an exogenous forcing process $z_{t+1} = \rho z_t + \sigma \epsilon_{t+1}$ with a finite-state Markov chain, e.g. by Tauchen's or Rouwenhorst's method. What would be a good way to discretize AR(1) process with stochastic volatility?

That is, if the original AR(1)-SV process looks something like this: $$ \begin{split} z_{t+1} &= \rho_z z_t + \mathrm{e}^{v_t} \sigma_z \epsilon_{t+1} \\ v_{t+1} &= \rho_v v_t + \sigma_v \eta_{t+1} \end{split} $$ with $\epsilon, \eta$ being independent standard gaussian shocks, my goal is to obtain values $z_i$ and transition probabilities $p_{ij}$ ($i,j=1,\dots,n$) such that the corresponding Markov chain approximates the original continuous-valued process.

Answer 1

You could do a binomial-tree approximation to the process for $z_t$ and then have a different process control the number of steps you take on the tree. This preserves recombining property and it is essentially the method explored in On the Computation of Continuous Time Option Prices Using Discrete Approximations (Amin (1991))

We develop a class of discrete, path-independent models to compute prices of American options within the Black-Scholes (1973) framework, including models in which state variables have time-varying volatility functions and models with multiple state variables. Time-varying volatility functions are illustrated with applications to term structure models developed by Vasicek (1977) and Heath, Jarrow, and Morton (1988), (1990). Distinct from previous work in the literature, the multivariate models suggested in this paper are consistent with arbitrarily large, though constant, covariance functions. Finally, we compare and contrast the numerical accuracy of a large number of models with simulation results.

You could then use a binomial AR(1) process for the number of steps.