Solution (source code)

= Solution

For the causal case $|\phi|<1$, expand
$$
\frac{1+\theta B}{1-\phi B}
=1+\sum_{j\ge1}(\phi+\theta)\phi^{j-1}B^j.
$$
Hence the <white noise> coefficients are $\psi_0=1$ and $\psi_j=(\phi+\theta)\phi^{j-1}$ for $j\ge1$. Orthogonality of <white noise> gives $\gamma_k=\sigma^2\sum_{j\ge0}\psi_j\psi_{j+k}$ for $k\ge0$. Summing the <geometric series> yields the <autocovariance of a causal ARMA(1,1) process>:
$$
\boxed{\gamma_0=\frac{\sigma^2(1+\theta^2+2\phi\theta)}{1-\phi^2},
\qquad
\gamma_k=\frac{\sigma^2(\phi+\theta)(1+\phi\theta)}{1-\phi^2}
\phi^{k-1}\quad(k\ge1),\qquad\gamma_{-k}=\gamma_k.}
$$
For $\phi=0$, this means $\gamma_1=\sigma^2\theta$ and $\gamma_k=0$ for $k\ge2$, together with $\gamma_0=\sigma^2(1+\theta^2)$. For $\theta=-\phi$, all positive-lag <covariances> vanish and $\gamma_0=\sigma^2$, as expected from cancellation of the two filter polynomials.