Summary: This paper proposes the discrete analytical method (DAM) to determine exactly and efficiently the fully nonstationary random responses of rectangular Kirchhoff plates under temporally and spectrally nonstationary acceleration excitation of earthquake ground motions. First, the fully nonstationary power spectral density (PSD) model is suggested by replacing the filtered frequency and damping of Gaussian filtered white-noise model with the time-variant ones. The exact solutions of free vibration of thin plates with two opposite edges simply supported boundary conditions are introduced. Then, the full analytical procedure for random vibration analysis of the plate is established by using a pseudo excitation method (PEM) that can consider all modal auto-correlation and cross-correlation terms. Owing to involving a series of Duhamel time integrals of single degree of freedom systems, it is difficult to fully analytically evaluate the PSD of time-variant responses such as the transverse deflection, velocity, acceleration and stress components. Thus, DAM that combines the PEM with precise integration technique is developed to enhance the computational efficiency. Finally, comparison of the results by the DAM with Monte Carlo simulations and the analytical stationary random vibration analysis demonstrates the high efficiency and accuracy of DAM. Moreover, the fully nonstationary excitation imposes a remarkable effect on the response PSD of rectangular Kirchhoff plates.