I'm looking for methods to derive an analytical formula for the variance of an ARMA(2,1) model.
Given $$(1-\phi_1B-\phi_2B^2)Z_t=(1-\theta_1B)a_t, $$
I've tried the usual "trick" $$\text{Var}[Z_t]=\text{E}[Z_t(\phi_1Z_{t-1}+\phi_2Z_{t-2} + a_t - \theta_1a_{t-1})], $$ which naturally gives me a relationship between the variance and the covariances $$\text{Var}[Z_t]=\gamma_0 \sim \gamma_1 + \gamma_2 + \text{const.}$$
I've tried solving this recursive relationship for $\gamma_0$ by generalizing the relationship to e.g. $$\gamma_k \sim \gamma_{k+1} + \gamma_{k+2} + \text{const.},$$ and writing out formulas for $k=\pm1$. But this introduces another unknown covariance $\gamma_3$.
I've also tried to expand $Z_t$ down to a known time value, say $t_0$. But this approach creates a messy formula, which I'm not able to evaluate for the variance.
Are there any other approaches that results in an analytical formula for the variance of an ARMA(2,1) model?