How is the MA part of ARMA solved for? In an AR model the coefficients on the lags can be solved for using least squares.  How is the MA part of ARMA solved for?  Since the MA part is a sum of white noise terms I imagine that it is not solved using least squares.  
 A: 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.
A: *

*Maximum likelihood estimation is perhaps the most popular method. Even if you noise is not Normally distributed, you will still obtain the Quasi-MLE which is consistent for the parameter of interest. You can check the Hamilton's Time Series Analysis to see how the likelihood function is specified.

*Another method mentioned in the links provided by Ram is Yule-Walker method which is based on the method of moments estimation (GMM). The idea is to match moments such as mean, variance and autocovariances.
