The Krylov subspace projection approach is a well-established tool for the reduced order modeling of dynamical systems in the time domain. In this paper, we address the main issues obstructing the application of this powerful approach to the time-domain solution of exterior wave problems. We use frequency independent perfectly matched layers to simulate the extension to infinity. Pure imaginary stretching functions based on Zolotarev's optimal rational approximation of the square root are implemented leading to perfectly matched layers with a controlled accuracy over a complete spectral interval of interest. A new Krylov-based solution method via stability-corrected operator exponents is presented which allows us to construct reduced-order models (ROMs) that respect the delicate spectral properties of the original scattering problem. The ROMs are unconditionally stable and are based on a renormalized bi-Lanczos algorithm. We give a theoretical foundation of our method and illustrate its performance through a number of numerical examples in which we simulate 2D electromagnetic wave propagation in unbounded domains, including a photonic waveguide example. The new algorithm outperforms the conventional finite-difference time domain method for problems on large time intervals.