We present a stochastic optimization method that uses a fourth-order regularized model to find local minima of smooth and potentially non-convex objective functions. This algorithm uses sub-sampled derivatives instead of exact quantities and its implementation relies on tensor-vector products only. The proposed approach is shown to find an $(\epsilon_1,\epsilon_2)$-second-order critical point in at most $\bigO\left(\max\left(\epsilon_1^{-4/3}, \epsilon_2^{-2}\right)\right)$ iterations, thereby matching the rate of deterministic approaches. Furthermore, we discuss a practical implementation of this approach for objective functions with a finite-sum structure, as well as characterize the total computational complexity, for both sampling with and without replacement. Finally, we identify promising directions of future research to further improve the complexity of the discussed algorithm.