# Algorithm for producing a Moving Average (as in ARIMA) model

I have a time series $$X_t$$ and I want to produce an ARMA forecast (without using any automated packages - the purpose of my project is to understand how those work).

So far, I have the AR(p) part down, I simply use p-many previous time steps as predictors and perform a linear regression for the next time step. My question is, how do I proceed with the MA(q) part?

Clearly, I don't have the errors of the previous time steps until I actually make a forecast and retrieve residuals. So, should I fit the AR(p) model, get the residuals, consider them as estimates for the error terms and use those to perform a second regression, and thus get the MA(q) part? Is there some other iterative process? I am trying to find papers or source codes demonstrating very explicitly how the numerical algorithm for calculating ARIMAs work.