Examples where method of moments can beat maximum likelihood in small samples? Maximum likelihood estimators (MLE) are asymptotically efficient; we see the practical upshot in that they often do better than method of moments (MoM) estimates (when they differ), even at small sample sizes 
Here 'better than' means in the sense of typically having smaller variance when both are unbiased, and typically smaller mean square error (MSE) more generally.
The question, occurs, however: 
Are there cases where the MoM can beat the MLE - on MSE, say - in small samples?
(where this isn't some odd/degenerate situation - i.e. given that conditions for ML to exist/be asymptotically efficient hold) 
A followup question would then be 'how big can small be?' - that is, if there are examples, are there some which still hold at relatively large sample sizes, perhaps even all finite sample sizes? 
[I can find an example of a biased estimator that can beat ML in finite samples, but it isn't MoM.]

Note added retrospectively: my focus here is primarily on the univariate case (which is actually where my underlying curiosity is coming from). I don't want to rule out multivariate cases, but I also don't particularly want to stray into extended discussions of James-Stein estimation. 
 A: Additional sources in favor of MOM:
Hong, H. P., and W. Ye. 2014. Analysis of extreme ground snow loads for Canada using snow depth records. Natural Hazards 73 (2):355-371.

The use of
MML could give unrealistic predictions if the sample size is small (Hosking et al. 1985;
Martin and Stedinger 2000).


Martins, E. S., and J. R. Stedinger. 2000. Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data. Water Resources Research 36 (3):737-744.
Abstract:

The three-parameter generalized extreme-value (GEV) distribution has found wide application for describing annual floods, rainfall, wind speeds, wave heights, snow depths, and other maxima. Previous studies show that small-sample maximum-likelihood estimators (MLE) of parameters are unstable and recommend L moment estimators. More recent research shows that method of moments quantile estimators have for −0.25 < κ < 0.30 smaller root-mean-square error than L moments and MLEs. Examination of the behavior of MLEs in small samples demonstrates that absurd values of the GEV-shape parameter κ can be generated. Use of a Bayesian prior distribution to restrict κ values to a statistically/physically reasonable range in a generalized maximum likelihood (GML) analysis eliminates this problem. In our examples the GML estimator did substantially better than moment and L moment quantile estimators for − 0.4 ≤ κ ≤ 0.

In the Introduction and Literature Review sections they cite additional papers which concluded that MOM in some cases outperform MLE (again extreme value modelling), e.g.

Hosking et al. [1985a] show that small-sample MLE parameter
estimators are very unstable and recommend probability-weighted moment (PWM) estimators
which are equivalent to L moment estimators [Hosking, 1990].
[...]
Hosking et al. [1985a] showed that the probability-weighted moments (PM) or equivalent L moments(LM) estimators for the GEV distribution are better
than maximum-likelihood estimators (MLE) in terms of bias
and variance for sample sizes varying from 15 to 100. More
recently, Madsen et al. [1997a] showed that the method of
moments (MOM) quantile estimators have smaller RMSE
(root-mean-squareer ror) for -0.25 < K < 0.30 than LM and
MLE when estimating the 100-year event for sample sizes of
10-50. MLEs are preferable only when K > 0.3 and the sample
sizes are modest (n >= 50).

K (kappa) is the shape parameter of GEV.
papers which appear in the quotes:
Hosking J, Wallis J, Wood E (1985) Estimation of the generalized extreme-value distribution by the method of probability-weighted moments. Technometrics 27:251–261.
Madsen, H., P. F. Rasmussen and D. Rosbjerg (1997) Comparison of annual
maximum series and partial duration series methods for modeling
extreme hydrologic events, 1, At-site modeling, Water Resour. Res.,
33(4), 747-758.
Hosking, J. R. M., L-moments: Analysis and estimation of distributions
using linear combinations of order statistics, J. R. Stat. Soc., Ser. B,
52, 105-124, 1990.

Additionally, I have the same experience as concluded in the above papers, in case of modeling extreme events with small and moderate sample size (<50-100 which is typical) MLE can give unrealistic results, simulation shows that MOM is more robust and has smaller RMSE.
A: In the process of answering this:  Estimating parameters for a binomial
I stumbled over this paper:
Ingram Olkin, A John Petkau, James V Zidek:  A comparison of N estimators for the Binomial Distribution.  Jasa 1981.
which gives an example where method of moments, at least in some cases, beats maximum likelihood.  The problem is estimation of $N$ in the binomial distribution $\text{Bin}(N,p)$ where both parameters are unknown.  It appears for example in trying to estimate animal abundance when you cannot see all the animals, and the sighting probability $p$ also is unknown.
A: This may be considered... cheating, but the OLS estimator is a MoM estimator. Consider a standard linear regression specification (with $K$ stochastic regressors, so magnitudes are conditional on the regressor matrix), and a sample of size $n$. Denote $s^2$ the OLS estimator of the variance $\sigma^2$ of the error term. It is unbiased so
$$ MSE(s^2) = \operatorname {Var}(s^2) = \frac {2\sigma^4}{n-K} $$
Consider now the MLE of $\sigma^2$. It is
$$\hat \sigma^2_{ML} = \frac {n-K}{n}s^2$$
Is it biased. Its MSE is
$$MSE (\hat \sigma^2_{ML}) = \operatorname {Var}(\hat \sigma^2_{ML}) + \Big[E(\hat \sigma^2_{ML})-\sigma^2\Big]^2$$
Expressing the MLE in terms of the OLS and using the expression for the OLS estimator variance we obtain
$$MSE (\hat \sigma^2_{ML}) = \left(\frac {n-K}{n}\right)^2\frac {2\sigma^4}{n-K} + \left(\frac {K}{n}\right)^2\sigma^4$$
$$\Rightarrow MSE (\hat \sigma^2_{ML}) = \frac {2(n-K)+K^2}{n^2}\sigma^4$$
We want the conditions (if they exist) under which
$$MSE (\hat \sigma^2_{ML}) > MSE (s^2) \Rightarrow \frac {2(n-K)+K^2}{n^2} > \frac {2}{n-K}$$
$$\Rightarrow 2(n-K)^2+K^2(n-K)> 2n^2$$
$$ 2n^2 -4nK + 2K^2 +nK^2 - K^3 > 2n^2 $$
Simplifying we obtain
$$ -4n + 2K +nK - K^2 > 0 \Rightarrow K^2 - (n+2)K + 4n < 0 $$
Is it feasible for this quadratic in $K$ to obtain negative values? We need its discriminant to be positive. We have
$$\Delta_K = (n+2)^2 -16n = n^2 + 4n + 4 - 16n = n^2 -12n + 4$$
which is another quadratic, in $n$ this time. This discriminant is
$$\Delta_n = 12^2 - 4^2 = 8\cdot 16$$
so
$$n_1,n_2 = \frac {12\pm \sqrt{8\cdot 16}}{2} = 6 \pm 4\sqrt2 \Rightarrow n_1,n_2 = \{1, 12\}$$
to take into account the fact that $n$ is an integer. If $n$ is inside this interval we have that  $\Delta_K <0$ and the quadratic in $K$ takes always positive values, so we cannot obtain the required inequality. So: we need a sample size larger than 12. 
Given this the roots for $K$-quadratic are
$$K_1, K_2 = \frac {(n+2)\pm \sqrt{n^2 -12n + 4}}{2} = \frac n2 +1 \pm \sqrt{\left(\frac n2\right)^2 +1 -3n}$$
Overall : for sample size $n>12$ and number of regressors $K$ such that $\lceil K_1\rceil <K<\lfloor K_2\rfloor $
we have
$$MSE (\hat \sigma^2_{ML}) > MSE (s^2)$$
For example, if $n=50$ then one finds that the number of regressors must be $5<K<47$ for the inequality to hold. It is interesting that for small numbers of regressors the MLE is better in MSE sense.
ADDENDUM
The equation for the roots of the $K$-quadratic can be written
$$K_1, K_2 = \left(\frac n2 +1\right) \pm \sqrt{\left(\frac n2 +1\right)^2  -4n}$$
which by a quick look I think implies that the lower root will always be $5$ (taking into account the "integer-value" restriction) -so MLE will be MSE-efficient when regressors are up to $5$ for any (finite) sample size.
A: An example that is admittedly connected with the James-Stein phenomenon, albeit in dimension one.
In the case of estimating the squared norm $\theta=||\mu||^2$ of a Gaussian mean vector, when observing $X\sim\mathcal N_p(\mu,\mathbf I_p)$, the MLE
$$\hat\theta^\text{MLE}=||x||^2$$
is doing quite poorly [in terms of squared error loss] when compared with the moment estimator
$$\hat\theta^\text{MM}=||x||^2-p$$
itself outperformed by the left-truncated version
$$\hat\theta^\text{TMM}=(||x||^2-p)^+$$
Surprisingly, the MLE of $\theta$ based on the original MLE distribution $$\hat\theta^\text{MLE}\sim\chi^2_p(\theta)$$
is different and apparently admissible, standing between $\hat\theta^\text{TMM}$ and $$\hat\theta^\text{JS}=(||x||^2-p+1)^+$$
A: "In this article, we consider a new parametrization of the two-parameter Inverse
    Gaussian distribution. We find the estimators for parameters of the Inverse Gaussian
    distribution by the method of moments and the method of maximum likelihood. Then, we
    compare the efficiency of the estimators for the two methods based on their bias and
    mean square error (MSE). For this we fix values of parameters, run simulations, and
    report MSE and bias for estimates obtained by both methods. The conclusion is that when 
    sample sizes are 10, the method of moments tends to be more efficient than the maximum
    likelihood method for estimates of both parameters (lambda and theta)...." read more
Nowadays one cannot (or should not) trust everything published, but the paper's last page appears promising. I hope this addresses your note added retrospectively.
A: The method of moments (MM) can beat the maximum likelihood (ML) approach when it is possible to specify only some population moments. If the distribution is ill-defined, the ML estimators will not be consistent. 
Assuming finite moments and i.i.d observations, the MM can provide good estimators with nice asymptotic properties. 
Example: Let $X_1, \ldots, X_n$ be an i.i.d sample of $X \sim f$, where $f: \mathbb{R} \to \mathbb{R}_+$ is an unknown probability density function. Define $\nu_k = \int_{\mathbb{R}} x^k f(x)dx$ the $k$th moment and consider that the interest is to estimate the forth moment $\nu_4$.
Let $\bar{X_k} = \frac{1}{n}\sum_{i=1}^n X_i^k$, then by assuming that $\nu_8 < \infty$, the central limit theorem guarantees that
$$
\sqrt{n}(\bar{X_4} - \nu_4) \stackrel{d}{\to} N(0, \nu_8 - \nu_4^2),
$$ where "$\stackrel{d}{\to}$" means "converges in distribution to". Moreover, by the Slutsky's theorem, 
$$
\frac{\sqrt{n}(\bar{X_4} - \nu_4)}{\sqrt{\bar{X_8} - \bar{X_4}^2}} \stackrel{d}{\to} N(0, 1)
$$ since $\bar{X_8} - \bar{X_4}^2 \stackrel{P}{\to} \nu_8 - \nu_4^2$ (convergence in probability).
That is, we can draw (approximate) inferences for $\nu_4$ by using the moment approach (for large samples), we just have to make some assumptions on the population moments of interest. Here, the maximum likelihood estimators cannot be defined without knowing the shape of $f$.
A Simulation study:
Patriota et al. (2009) conducted some simulation studies to verify the rejection rates of hypothesis testings in an errors-in-variables model. The results suggest that the MM approach produces error rates under the null hypothesis closer to the nominal level than the ML one for small samples. 
Historical note:
The method of moments was proposed by K. Pearson in 1894 "Contributions to the Mathematical Theory of Evolution". The method of maximum likelihood was proposed by R.A. Fisher in 1922 "On the Mathematical Foundations of Theoretical Statistics". Both papers where published in the Philosophical Transactions of the Royal Society of London, Series A. 
Reference:
Fisher, RA (1922). On the Mathematical Foundations of Theoretical Statistics, Philosophical Transactions of the Royal Society of London, Series A, 222, 309-368.
Patriota, AG, Bolfarine, H, de Castro, M (2009). A heteroscedastic structural errors-in-variables model with equation error, Statistical Methodology 6 (4), 408-423 (pdf)
Pearson, K (1894). Contributions to the Mathematical Theory of Evolution, Philosophical Transactions of the Royal Society of London, Series A, 185, 71-110.
A: I found one:
For the asymmetric exponential power distribution
$$f(x) = \frac{\alpha}{\sigma\Gamma(\frac{1}{\alpha})} \frac{\kappa}{1+\kappa^2}\exp\left(-\frac{\kappa^\alpha}{\sigma^\alpha}[(x-\theta)^+]^\alpha -\frac{1}{\kappa^\alpha \sigma^\alpha}[(x-\theta)^-]^\alpha\right)\,,\quad \alpha,\sigma,\kappa>0, \text{ and } x,\theta\in \mathbb R$$
the simulation results of Delicado and Goria (2008) suggest that for some of the parameters at the smaller sample sizes, method of moments can outperform MLE; for example in the known-$\theta$ case at sample size 10, when estimating $\sigma$, the MSE of MoM  is smaller than for ML.
Delicado and Goria (2008),
A small sample comparison of maximum likelihood, moments and L-moments methods for the
asymmetric exponential power distribution,
Journal Computational Statistics & Data Analysis
Volume 52 Issue 3, January, pp 1661-1673 
(also see http://www-eio.upc.es/~delicado/my-public-files/LmomAEP.pdf)
A: According to simulations run by Hosking and Wallis (1987) in "Parameter and Quantile Estimation for the Generalized Pareto Distribution", the parameters of the two-parameter generalized Pareto distribution given by the cdf
$G(y)=  \begin{cases}  1-\left(1+ \frac{\xi y}{\beta} \right)^{-\frac{1}{\xi}} & \xi \neq 0 \\
1-\exp\left(-\frac{y}{\beta}\right) & \xi=0 \end{cases}$
or the density 
$g(y)=  \begin{cases}  \frac{1}{\beta} \left( 1+\frac{\xi y}{\beta}      \right)^{-1-\frac{1}{\xi}} & \xi \neq 0 \\
 \frac{1}{\beta} \exp\left(-\frac{y}{\beta}   \right)    & \xi=0 
\end{cases}$
are more reliable if they are estimated by means of MOM as opposed to ML. This holds for samples up to size 500. The MOM estimates are given by
$\widehat\beta = \frac{\overline y \overline{y^2}}{2(\overline{y^2} - (\overline y)^2)}$
and 
$\widehat\xi = \frac{1}{2} - \frac{(\overline y)^2}{2(\overline{y^2} - (\overline y)^2)}$
with 
$\overline{y^2} = \frac{1}{n} \sum_{i=1}^n y_i^2$
The paper contains quite a few typos (at least my version does). Results for the MOM estimators given above were kindly provided by "heropup" in this thread.
