A novel spectral incremental dynamic analysis methodology for analysing structural response in nonlinear systems with fractional derivative elements is presented, aligning with modern seismic design codes, like Eurocode 8. Drawing inspiration from the concept of fully non-stationary stochastic processes, the vector of the imposed seismic excitations is characterised by time and frequency evolving power spectra stochastically compatible with elastic response spectra of specified damping ratio and ground acceleration. The proposed method efficiently determines the nonlinear system time-dependent probability density functions for the non-stationary system response amplitude by employing potent nonlinear stochastic dynamics concepts, such as stochastic averaging and statistical linearization. Unlike traditional incremental dynamic analysis curves found in the literature, the herein proposed method introduces a three-dimensional alternative counterpart, that of stochastic engineering demand parameter surfaces, providing with higher-order statistics of the system response. An additional noteworthy aspect involves the derivation of response evolutionary power spectra as function of spectral acceleration, offering a deeper insight into the underlying system dynamics. Besides its capabilities, the method maintains the coveted element of a particularly low associated computational cost, increasing its attractiveness and practicality among diverse applications of engineering interest. Numerical examples comprising the bilinear hysteretic model endowed with fractional derivative elements subject to an Eurocode 8 elastic design spectrum demonstrate the capabilities and reliability of the proposed methodology. Its accuracy is assessed by juxtaposing the derived results with germane Monte Carlo Simulation data.