Preliminary estimates of ARIMA in R?

We know that dealing with model involving MA factors is not easy to estimate, since there are past values of errors to be computed recursively. And this recursive estimation requires preliminary (initial) estimates of the parameters. For example, an ARMA(1,2) $$z_t=\phi z_{t-1}-\theta_1 \varepsilon_{t-1}-\theta_2 \varepsilon_{t-2}+\varepsilon_t$$ To estimate the parameters, we need to compute first the values of $\varepsilon_{t-1}$ and $\varepsilon_{t-2}$, since these are not available yet. And they are computed using $$\varepsilon_t=z_t-\phi z_{t-1}+\theta_1 \varepsilon_{t-1}+\theta_2 \varepsilon_{t-2}$$ Procedures for obtaining preliminary estimates of the parameters is available in Box and Jenkins, Time Series: Forecasting and Control. And this estimation is already available in much statistical software. My question is, "Is there a function for obtaining a preliminary estimate of the parameters in R?"

I need this to obtain a preliminary estimate for my Space-Time ARIMA. Another question is, "How does arima function of R compute preliminary estimates of the parameters when there are MA factors involved?"