We present multivariate unbiased estimators for second, third, and fourth order cumulants $C_2(x,y)$, $C_3(x,y,z)$, and $C_4(x,y,z,w)$. Many relevant new estimators are derived for cases where some variables are average-free or pairs of variables have a vanishing second order cumulant. The well-know Fisher k-statistics is recovered for the single variable case. The variances of several estimators are explicitly given in terms of higher order cumulants and discussed with respect to random processes that are predominately Gaussian. We surprisingly find that the frequently used third order estimator $\overline{x^3}$ for $C_3(x,x,x)$ of a process $x$ with zero average is outperformed by alternative estimators. The new (Gauss-optimal) estimator $\overline{x^3} - 3 \overline{x^2}\overline{x}(m-1)/(m+1)$ improves the variance by a factor of up to $5/2$. Similarly, the estimator $\overline{x^2 z}$ for $C_3(x,x,z)$ can be replaced by another Gauss-optimal estimator. The known estimator $\overline{xyz}$ for $C_3(x,y,z)$ as well as previously known estimators for $C_2$ and $C_4$ of one average-free variable are shown to be Gauss-optimal. As a side result of our work we present two simple recursive formulas for finding multivariate cumulants from moments and vice versa.