A stochastic method for making earthquakes compatible with target response spectrum was proposed. In the meantime, the earthquakes had fully nonstationary characteristics. First, the difference spectrum was gotten by comparing the response spectra of a recorded ground motion and the design response spectrum, then generated difference earthquakes compatible with difference spectrum through spectral representation method; second, the probability density function of lognormal distribution was adopted to express intensity nonstationary characteristics; finally, obtained spectral compatible fully nonstationary earthquakes by adding different earthquakes to the recorded earthquake. The frequency nonstationary characteristics of simulated earthquakes were close to ground motion records, and the response spectrum of generated earthquakes fitted the target response spectrum well. Hence, the method can be used in structural seismic analysis.