  We formulate standard and multilevel Monte Carlo methods for the $k$th moment
$\mathbb{M}^k_\varepsilon[\xi]$ of a Banach space valued random variable
$\xi\colon\Omega\to E$, interpreted as an element of the $k$-fold injective
tensor product space $\otimes^k_\varepsilon E$. For the standard Monte Carlo
estimator of $\mathbb{M}^k_\varepsilon[\xi]$, we prove the $k$-independent
convergence rate $1-\frac{1}{p}$ in the $L_q(\Omega;\otimes^k_\varepsilon
E)$-norm, provided that (i) $\xi\in L_{kq}(\Omega;E)$ and (ii)
$q\in[p,\infty)$, where $p\in[1,2]$ is the Rademacher type of $E$. By using the
fact that Rademacher averages are dominated by Gaussian sums combined with a
version of Slepian's inequality for Gaussian processes due to Fernique, we
moreover derive corresponding results for multilevel Monte Carlo methods,
including a rigorous error estimate in the $L_q(\Omega;\otimes^k_\varepsilon
E)$-norm and the optimization of the computational cost for a given accuracy.
Whenever the type of the Banach space $E$ is $p=2$, our findings coincide with
known results for Hilbert space valued random variables.
  We illustrate the abstract results by three model problems: second-order
elliptic PDEs with random forcing or random coefficient, and stochastic
evolution equations. In these cases, the solution processes naturally take
values in non-Hilbertian Banach spaces. Further applications, where physical
modeling constraints impose a setting in Banach spaces of type $p<2$, are
indicated.
