Calculation of nonlinear spectra of chromophore aggregates using response function theory when the number of contributing chromophores is large, and the level of excitation is high is extremely complicated. The main limitation is due to the exponential growth of computational time due to the aggregate size and number of excitations when considering an arbitrary excitation intensity. Non-perturbative calculation of spectra in this case becomes advantageous. We revisit our proposed model with exciton - exciton annihilation terms and apply it to large aggregates. We generalize the equations for both paulions and bosons with a parameter that allows smooth transition from one description to another. Intermediate statistics may also be valuable as molecular electronic excitations do not strictly obey either boson or paulion statistics. Specific approximations allow efficient calculation of pump-probe spectra for a large J aggregate.