  A new implementation of stochastic many-body perturbation theory for periodic
2D systems is presented. The method is used to compute quasiparticle
excitations in twisted bilayer phosphorene. Excitation energies are studied
using stochastic $G_0W_0$ and partially self-consistent $\bar \Delta GW_0$
approaches. The approach is inexpensive; it is used to study twisted systems
with unit cells containing $>2,700$ atoms ($>13,500$ valence electrons), which
corresponds to a minimum twisting angle of $\approx 3.1^\circ$. Twisted
bilayers exhibit band splitting, increased localization and formation of
localized Moir\'e impurity states, as documented by band-structure unfolding.
Structural changes in twisted structures lift band degeneracies. Energies of
the impurity states vary with the twisting angle due to an interplay between
non-local exchange and polarization effects. The mechanisms of quasiparticle
energy (de)stabilization due to twisting are likely applicable to a wide range
of low-dimensional Moir\'{e} superstructures.
