We present theory for calculating the third-order non-linear response function of a molecular aggregate in the weak inter-chromophore coupling regime. This approach is based on the perturbative expansion of the system evolution with respect to the resonance coupling, while the system-bath interaction is treated non-perturbatively by means of cumulant expansion. An explicit expression for the time-resolved fluorescence signal is then obtained. This allows us to investigate the ultrafast time-dependent Stokes shift, signatures of coherent dynamics, and the excitonic polaron formation in the excited state of the aggregate. Numerical simulations of the time-resolved fluorescence spectra of a pair of coupled molecules demonstrate these effects.