Calculate $E [ Z_{t-1}| X_{t-1}]$ in an ARMA process Suppose I have an ARMA(1,1) model:
$$X_t = \phi X_{t-1} + Z_t + \theta Z_{t-1}, \quad Z_t \sim WN(0,\sigma^2)$$
Indeed, I want calculate $E[X_t | X_{t-1}]$. For this:
\begin{align}
E[X_t | X_{t-1}] &= \phi X_{t-1} + E[Z_t | X_{t-1}]  + \theta E [ Z_{t-1}| X_{t-1}]\\
&= \phi X_{t-1} + \theta E [ Z_{t-1}| X_{t-1}]
\end{align}
I used the fact that for $j>0$, $X_t$ and $Z_{t+j}$ are independent. Thus, they are mean independent: $E[Z_t|X_{t-1}]= E[Z_t]=0$. But I don't know how to calculate $E [ Z_{t-1}| X_{t-1}]$.
I don't know if this the best way to calculate $E[X_t | X_{t-1}]$. If there is another way to calculate, I would also like to know!

For the Gaussian ARMA model we can derive the conditional mean and variance using standard results for conditional moments for the multivariate normal distribution.  This requires use of the auto-correlation function for the ARMA model.  You can find details of how to compute the autocorrelation function for an ARMA model in Section 3.3 of Brockwell and Davis (1991).  Details of how to compute the conditional distribution and moments can be found in O'Neill (2021), which describes the functions in the ts.extend package.  First we derive the autocorrelation function $\gamma$ at the relevant lag values (zero and one in this case):
$$\gamma(k) \equiv \mathbb{Corr}(X_t,X_{t-k}).$$
Using these values we can then write the conditional expectation of interest as:
$$\mathbb{E}(X_t | X_{t-1} = x)
= \mu + \gamma(1) (x - \mu).$$
(In this equation I have included a mean term $\mu$ for the Gaussian ARMA model, even though this mean is zero in your question.)  The conditional expectation of interest to you can be computed using the functions in the ts.extend package.  Here is an example using the parameters $\mu = 0$, $\phi = 0.8$ and $\theta = -0.3$.
#Set model parameters
mu    <- 0
phi   <- 0.8
theta <- -0.3

#Compute autocorrelation
AUTOCORR  <- ts.extend::ARMA.autocov(n = 2, ar = phi, ma = theta, corr = TRUE)
CORR.LAG1 <- AUTOCORR[2]

#Compute and plot conditional expectation
x <- (-30:30)/10
E <- mu + CORR.LAG1*(x - mu)
plot(x, E, type = 'l', xlim = c(-3, 3), ylim = c(-3, 3),
     main = 'Conditional Expectation in Gaussian ARMA Model',
     xlab = 'X[t-1]', ylab = 'E(X[t] | X[t-1])')


