Box-Jenkins Forecasting With ARIMA(p,d,q) models

Question

I want to check that I understand the general theme of forecasting with ARIMA models using box-jenkins, so I am going to take an example and then proceed from there.

We will use $B$ notation for the backshift operator. Suppose that we have the time series given by, with $Z_n\sim WN(0,\sigma^2)$:

$$(1-B)^2X_t-\alpha (1-B)^2X_{t-1}=Z_t+\theta Z_{t-1}$$

Then we can see that this is an ARIMA(1,2,1) model. Now if we denote the forecast of $X_{n+1}$ at time $n$ by $X_{n+1}^n$ then in order to make this forecast using the Box-Jenkins methodology we simply use the above equation and set $Z_{n+1}=0$. So writing in terms of $X_{n+1}$ we have:

$$X_{n+1}=(2+\alpha)X_{n}-(2\alpha+1)X_{n-1}+\alpha X_{n-2}+Z_{n+1}+\theta Z_{n}$$

Now to make the forecast we set $Z_{n+1}=0$ to obtain:

$$X_{n+1}^n=(2+\alpha)X_{n}-(2\alpha+1)X_{n-1}+\alpha X_{n-2}+Z_{n+1}+\theta Z_{1}$$

We now notice that $X_n-X_n^{n-1}=Z_n$ (the difference between a forecast value and an observed value will be the WN observed), so we may write this as:

$$X_{n+1}^n=(2+\alpha)X_{n}-(2\alpha+1)X_{n-1}+\alpha X_{n-2}+Z_{n+1}+\theta (X_n-X_n^{n-1})$$

So I am wondering if we have an iterative process to find the 1-step forecasts, how do we find an initial value of this recurrence?

Now if this was an ARMA process we could use the Wold-Decomposition to write $X_t=\vartheta(B)Z_t$, seT future values of $Z_n=0$ and forecast that way. However, as this is an ARIMA model and so not causal do we just do this iteratively using the previous forecast. That is, suppose that we want to find $X_{n+2}^n$. Do we just have:

$$X_{n+2}= (2+\alpha)X_{n+1}-(2\alpha+1)X_{n}+\alpha X_{n-1}+Z_{n+2}+\theta Z_{n+1}$$

Then set $Z_{n+1}=Z_{n+2}=0$ to get a forecast of:

$$X_{n+2}^n= (2+\alpha)X_{n+1}-(2\alpha+1)X_{n}+\alpha X_{n-1},$$

then as we have no observed value $X_{n+1}$ do we just use $X_{n+1}^n$ instead to get a final forecast value of:

$$X_{n+2}^n= (2+\alpha)X_{n+1}^n-(2\alpha+1)X_{n}+\alpha X_{n-1}$$

Then we could obtain any $X_{n+h}^n$ iteratively using this process?

Finally, if I want to find the variance of the errors using this method I can define the error in the $h$-step forecast to be:

$$e_n(h)=X_{n+h}-X^n_{n+h}$$

So clearly $e_n(1)=X_{n+1}-X_{n+1}^n=Z_{n+1}$, which gives $Var(e_n(1))=\sigma^2$, now how do I find $e_n(h)$ in general? If we consider $e_n(2)$ do we have:

$$e_n(2)=X_{n+2}-X_{n+2}^2=Z_{n+2}+\theta Z_{n+1}?$$

I am quite confused by this last bit with the variance of the errors.