Molecular spectroscopy is an essential experimental technique in both fundamental physical sciences and applied research. For example, electronic spectroscopy, in which light induces a change in the electronic state of a molecule, is a powerful tool for studying light-matter interactions appearing in solar energy conversion and light-emitting devices. To explain experimental findings and compare them with theoretical predictions, scientists often compute properties that can only be inferred indirectly, rather than simulating the observed spectroscopic signals. This is partly due to challenges associated with the simulation of electronic spectra, which requires both an accurate description of the electronic structure and the quantum-dynamical treatment of the motion of atomic nuclei. As the exact methods to solve the underlying quantum-mechanical problem scale exponentially with the number of atoms, this task must be performed approximately for all but the smallest molecules. In this thesis, we develop a set of methods to systematically study the effects of different approximations on the computed spectra. As our starting point, we consider a semiclassical method called the thawed Gaussian approximation, which has resurfaced in recent years as an efficient and robust alternative to simple but crude models or highly accurate but costly quantum dynamics methods. First, we present several practical improvements of the method, such as generalization to non-Condon effects or the efficient single-Hessian version, which enable a broader set of molecules to be studied. Second, we introduce a more fundamental modification of the methods - the inclusion of non-zero temperature effects. The underlying theory, which is inspired by the concept of so-called thermo-field dynamics, is exact and applies to any quantum dynamics method. Remarkably, within the thawed Gaussian approximation, the inclusion of non-zero temperature comes at nearly no additional computational cost. These methods, developed originally for steady-state linear spectra, are then employed to compute state-of-the-art, time-resolved spectra, which are typically measured with ultrashort laser pulses. More precisely, we study how the anharmonicity of the potential energy surfaces and mode-mode (Duschinsky) couplings affect pump-probe and two-dimensional electronic spectra. This approach, based on the thawed Gaussian wavepacket propagation, is the first method that can exactly evaluate the nonlinear electronic spectra of many-dimensional, shifted, distorted, and Duschinsky-rotated harmonic potentials. As in the case of linear absorption and emission spectra, we again use the concept of thermo-field dynamics to include non-zero temperature effects in the simulation of two-dimensional electronic spectra. The resulting wavepacket picture of two-dimensional spectroscopy at arbitrary temperature could strongly impact how we interpret and simulate multidimensional spectra in the future.