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

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.

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))