Following recent advances in terahertz (THz) technology, there is a consensus on the crucial role of THz communications in the next generation of wireless systems. Aiming at catalyzing THz communications research, we propose TeraMIMO, an accurate stochastic MATLAB simulator of statistical THz channels. We simulate ultra-massive multiple-input multiple-output antenna configurations as critical infrastructure enablers that overcome the limitation in THz communications distances. We consider both line-of-sight and multipath components and propose frequency- and delay-domain implementations for single- and multi-carrier paradigms in both time-invariant and time-variant scenarios. We implement exhaustive molecular absorption computations based on radiative transfer theory alongside alternative sub-THz approximations. We further model THz-specific constraints, including wideband beam split effects, spherical wave propagation, and misalignment fading. We verify TeraMIMO by analogy with measurement-based channel models in the literature and ergodic capacity analysis. We introduce a graphical user interface and a guide for using TeraMIMO in THz channel generation and analyses.