  We propose two novel unbiased estimators of the integral
$\int_{[0,1]^{s}}f(u) du$ for a function $f$, which depend on a smoothness
parameter $r\in\mathbb{N}$. The first estimator integrates exactly the
polynomials of degrees $p<r$ and achieves the optimal error $n^{-1/2-r/s}$
(where $n$ is the number of evaluations of $f$) when $f$ is $r$ times
continuously differentiable. The second estimator is computationally cheaper
but it is restricted to functions that vanish on the boundary of $[0,1]^s$. The
construction of the two estimators relies on a combination of cubic
stratification and control ariates based on numerical derivatives. We provide
numerical evidence that they show good performance even for moderate values of
$n$.
