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.