Spectral analysis using overlapping sliding windows is among the most widely used techniques in analyzing non-stationary time series. Although sliding window analysis is convenient to implement, the resulting estimates are sensitive to the window length and overlap size. In addition, it undermines the dynamics of the time series as the estimate associated to each window uses only the data within. Finally, the overlap between consecutive windows hinders a precise statistical assessment. In this paper, we address these shortcomings by explicitly modeling the spectral dynamics through integrating the multitaper method with state-space models in a Bayesian estimation framework. The underlying states pertaining to the eigen-spectral quantities arising in multitaper analysis are estimated using instances of the Expectation-Maximization algorithm, and are used to construct spectrograms and their respective confidence intervals. We propose two spectral estimators that are robust to noise and are able to capture spectral dynamics at high spectrotemporal resolution. We provide theoretical analysis of the bias-variance trade-off, which establishes performance gains over the standard overlapping multitaper method. We apply our algorithms to synthetic data as well as real data from human EEG and electric network frequency recordings, the results of which validate our theoretical analysis.