We formulate standard and multilevel Monte Carlo methods for the kth moment M-k (epsilon)[xi] of a Banach space valued random variable xi : Omega -> E, interpreted as an element of the k -fold injective tensor product space circle times(k)(epsilon) E. For the standard Monte Carlo estimator of M-epsilon(k) [xi], we prove the k-independent convergence rate 1- 1/p in the L-q(Omega; circle times E-k(epsilon))-norm, provided that (i) xi E L-kq(Omega; E) and (ii) q is an element of [p,infinity), where p is an element of [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; circle times E-k(epsilon))-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.(c) 2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/).