Manual calculation of ARIMA(p,d,q) forecast

I have obtained the answer on my own question Manual calculation of ARIMA(1,1,0) forecast, then I have seen the question and answer Forecasting with ARMA models.

I'd like to write the formula to the manual calculation of ARIMA(p,d,q) forecast in common case. I'm starting with ARIMA(2,1,2) which is:

$$y_t = \beta_1(y_{t-1}-y_{t-2}-\mu)+\beta_2(y_{t-2}-y_{t-3}-\mu)-(y_{t-1}+ \mu)+\varepsilon_{t}+\alpha_1\varepsilon_{t-1}+\alpha_2\varepsilon_{t-2}$$ or $$y_t = \underbrace{\beta_1(\Delta y_{t-1}-\mu)+\beta_2(\Delta y_{t-2}-\mu)-(y_{t-1}+ \mu)}_{part~AR(2)}+\underbrace{\varepsilon_{t}+\alpha_1\varepsilon_{t-1}+\alpha_2\varepsilon_{t-2}}_{part~MA(2)},$$

Lets say we have n=252 observations and we want to calculate an one-step ahead ARIMA(2,1,2) forecast in-sample. This is easy in R:

library(forecast)
set.seed(1)
n<-252
y<- runif(n, 9000, 10000)
fit1<-Arima(y[1:(n-1)], order=c(2,1,2))
coef(fit1)
#            ar1         ar2         ma1         ma2
#    -0.49265434  0.07456189 -0.46064013 -0.53934992


We have fitted the ARIMA(2,1,2) model:

$$y_t = -0.49265434(y_{t-1}-y_{t-2}-\mu)+0.07456189(y_{t-2}-y_{t-3}-\mu)-(y_{t-1}+ \mu)+\varepsilon_{t} -0.46064013\varepsilon_{t-1}-0.53934992\varepsilon_{t-2}$$ The one-step ahead ARIMA(2,1,2) forecast in-sample is

forecast(fit1, h=1)$mean[1] # [1] 9538.81  Now I'd like to do the the one-step ahead forecast in-sample manually. We can estimate the mean,$\mu$, and$\varepsilon_t$: $$\hat{\mu}=\frac{1}{n-1}\sum_{i=1}^{n-1}y_i,$$ $$\hat{\varepsilon_t}=y_{t-1} -\beta_1 y_{t-2} - \beta_2 y_{t-3},$$ But now I have a problem how to write the estimations$\hat{\varepsilon}_{t-1}$and$\hat{\varepsilon}_{t-2}$. Question. Are my representation of the ARIMA(2,1,2) model, i.e.$y_t$, and estimations$\hat{\mu}$and$\hat{\varepsilon_t}\$ correct?