It may be due to superparamagnetic behaviour which you solve by increasing the size of the sample. Another thing you could try is to look at the instantaneous magnetization rather than a mean -- it may jump up and down between the top and bottom of the hysteresis but hopefully you'll see only the values for these two well defined states rather than intermediate values resulting from taking averages of them.
Hysteresis loops at finite temperature are notoriously difficult however, because of how small the real-time of the simulations are:
"Incorporating the effects of thermal fluctuations through the introduction of a stochastic field term in the LLG equation is used to study a wide range of applications in magnetic media. However, there are limitations imposed by the requirement of picosecond integration time steps to obtain reliable solutions to the LLG equation. This makes the direct study of dynamical magnetization reversal at experimental times scales simply not feasible in many cases. Even with the fastest processors, the overall time period for a simulation is in the order of μs, whereas in the case of experimental M−H loops, for example, simulation times of many minutes are required. "
--
https://ieeexplore.ieee.org/document/6612685If you look at different articles on/with spin-dynamics simulations it seems like people try to avoid doing hysteresis loop simulations at finite temperatures.