Quasi-periodic orbital cycles may be expressed as peaks in the power spectrum of seismic images of sedimentary deposits, if the bandwidth includes the relevant spatial frequencies and if cyclicity is preserved in the impedance. Several previous studies have shown apparent Milankovic cycles in seismic data1,2. In addition to the problem of spurious (non-orbital) peaks inherent to spectral analysis on core or borehole logs, seismic data suffers from being i) a representation of the reflectivity (a derivative of the impedance) ii) filtered by the source spectrum iii) subject to attenuation with propagation due to spreading and absorption and iv) contaminated by internal multiples. Here we present a method to assign significance levels to spectral peaks in seismic data. We create an ensemble of models for a given scenario, based on geological parameters (sedimentation rate), geophysical parameters (velocity and density logs, the source wavelet) and the relative contribution of background sedimentary noise versus orbital forcing (an astronomical solution). We model the borehole log and the seismic response, using a visco-acoustic scheme to include absorption and multiples. The output is an ensemble of synthetic logs and traces that do and do not contain cyclicity. We apply a spectral whitening procedure (adjusted to cope with an arbitrary seismic wavelet) and estimate the spectral baseline and corresponding significance thresholds. We compare the true and false positive detection rates across the spectrum at each significance to estimate the sensitivity and specificity for each orbital period, for both the borehole and seismic cases. We demonstrate the method on a synthetic scenario with constant sedimentation rate and a 150 Hz Ricker source. We also re-evaluate the real data case study of ODP Site 10821. Close to the dominant source frequency the statistical power of seismic cyclostratigraphy is comparable to - but always lower than - conventional (borehole) cyclostratigraphy. Our results suggest that for sedimentation rates above ~10 cm ka-1 and common seismic bandwidths this method may be suitable for detecting 100 kyr eccentricity and 41 kyr obliquity cycles directly from seismic data.