Conditional variance of ARMA(p,q) process $X_t$ is an ARMA(p,q) process. How do you derive the conditional variance on $X_s$ $\forall s<t$, i.e
$$Var({X_t|X_s, s<t})=?$$

I don't know if this is a standard procedure, but I think this should work:
Starting from the representation $\theta (B) X = \phi (B) \epsilon$ for $X \sim \operatorname{ARMA}(p,q)$, let's assume we know the roots of $\phi(z)$ (these are usually computed anyway), then we can represent $\phi(B) = c \prod_{k=1}^q (B - \lambda_k )$, and then using the Neumann series $$ (B - \lambda)^{-1} = \sum_{k=0}^\infty (-\lambda)^{-k-1} B^k , $$
we get the $\operatorname{AR}(\infty)$-representation
$$ \frac{\theta (B)}{\phi (B)} X=\epsilon \\
\left( c^{-1} \theta(B)\prod_{k = 1}^q \sum_{n_k=0}^\infty  (- \lambda_k)^{-n_k - 1} B^{n_k} \right) X = \epsilon \\
c^{-1}  \sum_{n = 0}^\infty \left( \sum_{n_1 + \cdots + n_q = n} \prod_{k=1}^q (- \lambda_k )^{-n_k-1} \right) B^{n} \theta(B)  X = \epsilon .
$$
This should be no problem to evaluate this numerically for some orders of $n$ up to $m$. Then, given such an $\operatorname{AR}(m)$ approximation $\widetilde \theta (B) X \approx \epsilon$, you can compute the conditional variance in a straight-forward manner:
$$X_{t+h} \vert \mathcal F_t = \dots + \epsilon_{t+h} + \widetilde \theta_1 \epsilon_{t+h-1} + (\widetilde \theta_1^2 + \widetilde \theta_2) \epsilon_{t+h - 2} + (\widetilde \theta_3 + \widetilde \theta_2 \widetilde \theta_1 + \widetilde \theta_1 \widetilde \theta_2 + \widetilde \theta_1^3) \epsilon_{t+h-3} + \dots ,  $$
yielding the conditional variance
$$\operatorname{var}(X_{t+h} \vert \mathcal F_t ) = \sigma_\epsilon^2 \left\lbrace \widetilde \theta_1^2 + (\widetilde \theta_1^2 + \widetilde \theta_2)^2 + (\widetilde \theta_3 + \widetilde \theta_2 \widetilde \theta_1 + \widetilde \theta_1 \widetilde \theta_2^2 + \widetilde \theta_1^3)^2 + \dots \right\rbrace $$
(terms up to $(\widetilde \theta_h + \widetilde \theta_{h-1} \widetilde  \theta_1 + \dots)^2$)
