Suppose I have a $p$'th order vector auto regression
$$\vec Z_t = F_1\vec Z_{t-1}+F_2\vec Z_{t-2} + \cdots +F_p \vec Z_{t - p} + \vec \epsilon_t,\qquad \vec\epsilon_t\sim N_q(\vec0,Q)$$
where $Z_t\in\mathbb{R}^q$. Then we can put this into the state-space form
\begin{aligned} \vec X_t &= (\vec Z_t,\vec Z_{t-1},\dots, \vec Z_{t - p + 1}) \\ \vec X_t &= F\vec X_{t-1} + R\vec\epsilon_t\\ F &= \begin{pmatrix} F_1 & \cdots & \cdots & F_{p-1} & F_p \\ I_q & 0 & \cdots & 0 & 0 \\ 0 & I_q & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots &0 & 0 \\ 0 & \cdots & 0 & I_q & 0 \end{pmatrix} & R &= \begin{pmatrix} I_q \\ 0 \\ \vdots \\ 0 \end{pmatrix} \end{aligned}
It follows that the conditional mean is
\begin{aligned} E(\vec X_{t + p} \mid \vec X_t = \vec x)&= E(\vec E(\vec X_{t + p} \mid \vec X_{t + p - 1}) \mid \vec X_t = \vec x) \\ &= E(F\vec X_{t + p -1} \mid \vec X_t = \vec x) = \cdots = F^p\vec x \end{aligned}
and for the conditional covariance, we use that
\begin{aligned} \vec Z_{t+h} &\propto \sum_{i = 1}^h G(h-i)\vec \epsilon_{t + i} \\ G(k) &=\left\{\begin{matrix} I_q & k = 0 \\ \sum_{i = 1}^{\min (k, p)}F_iG(k-i) & k > 0 \end{matrix}\right. \end{aligned}
Thus,
\begin{aligned} \text{Var}(\vec Z_{t+h}\mid \vec X_t) &= \sum_{i=1}^h G(h-i)QG(h-i)^\top \\ \text{cov}(\vec Z_{t+h}, \vec Z_{t+l}\mid \vec X_t) &= \text{cov}\left( \sum_{i = 1}^h G(h-i)\vec \epsilon_{t + i}, \sum_{i = 1}^l G(l-i)\vec \epsilon_{t + i} \mid \vec X_t\right) \\ &= \text{cov}\left( \sum_{i = 1}^l G(h-i)\vec \epsilon_{t + i}, \sum_{i = 1}^l G(l-i)\vec \epsilon_{t + i} \mid \vec X_t\right) \\ &= \sum_{i = 1}^l G(h-i)QG(l-i)^\top \end{aligned}
where I assume that $h > l$. I hope the above is correct though this is not the question. My question is what are similar expressions if we only condition on $Z_t$ and not $X_t= (\vec Z_t,\vec Z_{t-1}, \vec Z_{t - p + 1})$? I.e., what are
$$ E(\vec X_{t+p}\mid \vec Z_t),\qquad \text{Var}(\vec Z_{t+h}\mid \vec Z_t),\qquad \text{Cov}(\vec Z_{t+h},\vec Z_{t+l}\mid \vec Z_t),\qquad h > l > 0 $$
Update
Say
$$ \vec X_0 \sim N(\vec \mu_0, Q_0) $$
which we either set to be the stationary distribution or a distribution selected for convenience. Then
$$ \vec Z_k = \sum_{i = 1}^k G(k-i)\vec \epsilon_i + R^\top F^k\vec X_0 $$
and thus,
\begin{align*} E(\vec Z_k) &= R^\top F^k\vec\mu_0 \\ \text{Var}(\vec Z_k) &= \sum_{i = 1}^k G(k-i)QG(k-i)^\top + R^\top F^kQ_0F^{k\top} R \\ \text{Cov}(\vec Z_k, \vec Z_l) &= \sum_{i = 1}^l G(k-i)QG(l-i)^\top + R^\top F^kQ_0F^{l\top} R, & k &> l > 0 \end{align*}
Using the above, we can compute the joint mean and covariance matrix and find that
\begin{aligned} \begin{pmatrix} \vec X_{t+p} \\ \vec Z_t \end{pmatrix} = \begin{pmatrix} \vec Z_{t + p} \\ \vec Z_{t + p - 1} \\ \vdots \\ \vec Z_t \end{pmatrix} & \sim N\left(\vec\mu, \Sigma \right) \\ \vec\mu &= (E(\vec X_{t+p})^\top, E(\vec Z_t)^\top)^\top\\ \Sigma &= \begin{pmatrix} \text{Var}(\vec X_{t+p}) & \text{Cov}(\vec X_{t+p}, \vec Z_t) \\ \text{Cov}(\vec Z_t,\vec X_{t+p}) & \text{Var}(\vec Z_t) \end{pmatrix} \end{aligned}
From which it is follows that
\begin{aligned} \vec X_{t+p} \mid \vec Z_t = \vec z &\sim N(k_\vec z, K_\vec z) \\ k_\vec z &= E(\vec X_{t+p}) + \text{Cov}(\vec X_{t+p}, \vec Z_t) \text{Var}(\vec Z_t)^{-1}(\vec z - E(\vec Z_t)) \\ K_\vec z &= \text{Var}(\vec X_{t+p}) - \text{Cov}(\vec X_{t+p}, \vec Z_t) \text{Var}(\vec Z_t)^{-1} \text{Cov}(\vec Z_t, \vec X_{t+p}) \end{aligned}
Is the above correct? Can I simplify the final expression further to something that is faster to compute or neater? I gather the latter is simple when we have stationary model since the unconditional means and covariances are independent of time. What about in the non-stationary case where we assume some distribution for $\vec X_0$?
