This talk introduces the numerical solutions

of Caputo--Hadamard fractional stochastic differential equations. Firstly, we construct an Euler--Maruyama (EM) scheme for the equations, and the corresponding convergence rate is investigated. Secondly, we propose a fast EM scheme based on the sum-of-exponentials approximation to decrease the computational cost of the EM scheme. More concretely, the fast EM scheme reduces the computational cost from $O(N^2)$ to $O(N\log^2 N)$ and the storage from $O(N)$ to $O(\log^2 N)$ when the final time $T\approx e$, where $N$ is the total number of time steps. Moreover, considering the statistical errors from Monte Carlo path approximations, multilevel Monte Carlo (MLMC) techniques are utilized to reduce computational complexity. In particular, for a prescribed accuracy $\epsilon>0$, the EM scheme and the fast EM scheme, integrated with the MLMC technique, respectively reduce the standard EM scheme's computational complexity from $O(\epsilon^{-2-\frac{2}{\widetilde{\alpha}}})$ to $O(\epsilon^{-\frac{2}{\widetilde{\alpha}}})$ and the fast EM scheme's complexity to $O(\epsilon^{-\frac{1}{\widetilde{\alpha}}}\abs{\log \epsilon}^3)$, where $0<\widetilde{\alpha}=\alpha-\frac 12<\frac 12$. Finally, numerical examples are included to verify the theoretical results and demonstrate the performance of our methods.