* * Example of bootstrapping a spectral density, from Franke and * Hardle(1992) "On Bootstrapping Kernel Spectral Estimates", Annals of * Statistics, vol. 20, no 1, 121-145. * * Data set is the Wolfer Sunspot series * open data sunspots.dat calendar 1770 data(format=free,org=columns) 1770:1 1869:1 sunspots * diff(center) sunspots / cspots compute nobs=%allocend() * * <> is the window width to be used in estimation. * <> is the window width to be used to standardize the * periodogram for bootstrapping. These can be different. * compute nwidth=9 compute nwidthboot=9 * freq 5 nobs rtoc # cspots # 1 * * Compute the periodogram back into complex series 1. We need this for * bootstrapping. * fft 1 cmult(scale=1.0/(2*%pi*nobs)) 1 1 * * Get the smoothed estimate with the desired width to complex series 2. * window(type=tent,width=nwidth) 1 / 2 * * Get the smoothed estimate with the width for bootstrapping to complex * series 5. * window(type=tent,width=nwidthboot) 1 / 5 compute nboot=1000 * * The reshuffled periodogram needs to be symmetrical, so we only use * half the ordinates. Because the 0 ordinate in the periodogram is * forceably zero (as a result of de-meaning), we also leave that out. * compute nshuffle=nobs/2+1 * * We need to bootstrap the periodogram divided by the estimated * spectral density using the width <>. * set source = %real(%z(t,1))/%real(%z(t,5)) * set sumsqr 1 nobs = 0.0 * do draw=1,nboot * * Draw a set of standardized periodogram ordinates. Normalize to mean * one and symmetrize. * boot shuffle 2 nshuffle cset 3 1 nshuffle = %if(t==1,0.0,source(shuffle(t))) sstats(mean) 2 nshuffle %real(%z(t,3))>>shufflesum cset 3 1 nshuffle = %z(t,3)/shufflesum @cseriessymm 3 * * Recolor it by the standardizing spectral estimate. * cset 3 = %z(t,3)*%z(t,5) * * Smooth using the desired width. * window(type=tent,width=nwidth) 3 / 4 * * Get sample statistics with respect to log ratio to the original * spectral density estimate. * set sumsqr 1 nobs = sumsqr+log(%real(%z(t,4))/%real(%z(t,2)))^2 end do draw * set logspect 1 nobs = log(%real(%z(t,2))) set lower 1 nobs = logspect-2.0*sqrt(sumsqr/nboot) set upper 1 nobs = logspect+2.0*sqrt(sumsqr/nboot) graph(footer="Confidence bands using bootstrap") 3 # logspect # lower / 2 # upper / 2 * set loweradt = logspect+log(%edf/%invchisqr(.025,%edf)) set upperadt = logspect+log(%edf/%invchisqr(.975,%edf)) graph(footer="Confidence bands using asymptotic distribution") 3 # logspect # loweradt / 2 # upperadt / 2