The MA parameters are estimated in many different ways, including MLE. It involves solving a set of non-linear equations, which is why everyone eventually ends up resorting to numerical methods.
Here are a few links that should get you moving forward.
One good place to start is Prof. Hyndman's textbook: Forecasting: Principles and Practice
For this question, I'd start with Section 8.4, and then Section 8/7.
Some of the theory behind estimating the MA thetas is in this lecture:
http://www2.econ.iastate.edu/classes/econ674/bunzel/documents/Lecture4.pdf
Here's the idea: First, each error term is recursively estimated ($\epsilon_0$ is assumed to be zero.) by exploiting the fact that the $\epsilon$ 's are Normal white noise.
Even after this, you have to revert to numerical methods to get the MLE.
Specifically, look at the slides 52-55.
Wolfram does a good job of explaining this here. They assume you will be using $Mathematica$, but the examples are relevant even if you are not using it.
Specifically, look at the section on Innovations Algorithm, where they have an example. (One drawback is that the implementation details of the Innovations Estimation are not shared.)
If you believe that your noise is zero-mean and Normally distributed, then be sure to also read the section on "Maximum Likelihood Method."
Auto.Arima in R
This Journal of Statistical Computing paper is well worth reading. A practical way to get moving on finding the best $P,d,q$ for ARIMA is implemented in forecast and the pseudo-code is in the JSS paper.
Hope some of these help you move forward.