We develop a theoretical approach to study the transient dynamics and the time-dependent statistics for the Anderson-Holstein model in the regime of strong electron-phonon coupling. For this purpose we adapt a recently introduced diagrammatic approach to the time domain. The generating function for the time-dependent charge transfer probabilities is evaluated numerically by discretizing the Keldysh contour. The method allows us to analyze the system evolution to the steady state after a sudden connection of the dot to the leads, starting from different initial conditions. Simple analytical results are obtained in the regime of very short times. We study in particular the apparent bistable behavior occurring for strong electron-phonon coupling, small bias voltages, and a detuned dot level. The results obtained are in remarkably good agreement with numerically exact results obtained by quantum Monte Carlo methods. We analyze the waiting time distribution and charge transfer probabilities, showing that only a single electron transfer is responsible for the rich structure found in the short-time regime. A universal scaling (independent of the model parameters) is found for the relative amplitude of the higher order current cumulants in the short-time regime, starting from an initially empty dot. We finally analyze the convergence to the steady state of the differential conductance and of the differential Fano factor at the inelastic threshold, which exhibits a peculiar oscillatory behavior.