Molecular excitons in large extended systems are often not well described by local time-dependent density functional theory (TDDFT) due to highly delocalized states with long range electronic coupling. The issue of long-range coupling is made exceptionally more difficult when we consider excitons delocalized over many large molecules in aggregates ranging up to micron scale. In this thesis, we develop a series of electronic structure theory methods leveraging stochastic techniques that enable us to perform higher quality calculations on molecular excitons, and enable us to study extremely large systems in the context of molecular aggregates. We have developed a linear scaling method that can study spectroscopic observables such as the density of states and participation ratio in systems of millions of coupled dye dipoles. For the study of single excitons in large molecular complexes, we have developed a stochastic formalism of the Bethe-Salpeter equation, the linear response formalism that arises from the GW approximation of many-body perturbation theory. Through a series of algorithmic improvements to the method, we have developed new approximations to capture the screened Coulomb interaction at lower computational cost, leading to the study of systems with several thousand electrons.