## Bootstrapping confidence bands for spectrum

Questions and discussions on Time Series Analysis

### Bootstrapping confidence bands for spectrum

Hello.

I would like to estimate confidence bands for spectrum estimates. I have seen the discussion in a previous post on this section, but it relies on a Chi-sq assumption. Since I am somewhat short of time, I was wondering if there are routines or code readily available to obtain confidence bands obtained through bootstrapping along the lines of

J. Franke and W. Hardle (1992) “On Bootstrapping Kernel Spectral Estimates” The Annals of StatisticsVol. 20 No 1

Thanks for any help.
rod1

Posts: 5
Joined: Tue Sep 13, 2011 11:53 am

### Re: Bootstrapping confidence bands for spectrum

This looks like it does it. The smoothing that's used for doing the "source" for the bootstrap doesn't need to match the one that is used to actually smooth the spectral estimate. Some quick experiments would indicate that it isn't a good idea to make <<nwidthboot>> too large.

bootspectrum.rpf
Program file

sunspots.dat
Data file
TomDoan

Posts: 2725
Joined: Wed Nov 01, 2006 5:36 pm

### Re: Bootstrapping confidence bands for spectrum

Thank you very much, I am going to work with it.
rod1

Posts: 5
Joined: Tue Sep 13, 2011 11:53 am

### Re: Bootstrapping confidence bands for spectrum

Can't Find Match for %IF(INTEGER,COMPLEX,COMPLEX). Closest Match is %IF(REAL,REAL,REAL)
## SX27. Illegal Combination of Data Types for Operation
>>>><=n/2,cs,cs(n-t+2))<<<<
## SX21. A END or } Here is Unneeded or Unexpected

Also, it seems that I need to download cseriessymm from procedure and example list. My 7.3 does not have it.
ivory4

Posts: 149
Joined: Mon Aug 24, 2009 12:16 pm

### Re: Bootstrapping confidence bands for spectrum

To get CSERIESSYMM to work on RATS prior to v8 you would need to replace

cset cs 1 n = %if(t<=n/2,cs,cs(n-t+2))

with

cset cs 1 n/2 = cs
cset cs n/2+1 n = cs(n-t+2)
TomDoan

Posts: 2725
Joined: Wed Nov 01, 2006 5:36 pm

### Re: Bootstrapping confidence bands for spectrum

Hello, Tom.

Thanks for the code. I have used it subtituting for my data, and it works. What i didn't realize was that it computes the confidence bands and plots them in the time domain, and what I basically need is boostrapped confidence bands for the spectrum computed in the procedure @spectrum. I have tried to adapt @spectrum to produce confidence bands, to no avail. How should I proceed? Thanks again.
rod1

Posts: 5
Joined: Tue Sep 13, 2011 11:53 am

### Re: Bootstrapping confidence bands for spectrum

The SPECTRUM procedure does exactly the same thing (plots everything in the time domain). You'll have to post what you tried to do so that didn't work.
TomDoan

Posts: 2725
Joined: Wed Nov 01, 2006 5:36 pm

### Re: Bootstrapping confidence bands for spectrum

Sorry, I meant to say that I need to plot the smoothed spectrum and its confidence bands with fractions of pi on the horizontal axis, as the procedure @spectrum does.

I tried to add the boostrapping part of @bootspectrum to @spectrum, and adding some changes:

- plotting without log scale
- plotting the spectrum and confidence bands in the (0, pi) range

Below I enclose the code for this modified version of @spectrum. It does run, but the confidence band looks very wide.

I'm not sure whether what I did is right, if you could check it I would be very grateful

Regards.

Code:

Code: Select all
`*** @spectrum3( options )   series   start  end* Computes and optionally graphs the estimated spectrum of a series.** Options:*   ORDINATES=Number of ordinates [depends upon length,seasonal]*   WINDOW=[FLAT]/TENT/QUADRATIC*   WIDTH=Window width [depends upon ordinates]*   TAPER=TRAPEZOIDAL/COSINE/[NONE]*   WTAPER=Taper width [.25 of length]*   [GRAPH]/NOGRAPH*   HEADER=Header for high-resolution graph*   FOOTER=Footer for high-resolution graph*   [LOGSCALE]/NOLOGSCALE Graph spectrum in log scale or not.*   PERIODOGRAM/[NOPERIODOGRAM]*   SPECTRUM=series for [0,pi] spectrum (0 frequency at entry 1)*   FREQ= frequencies [0,pi] (0 frequency at entry 1)** Revision schedule*   04/1992 COMPUTE HALF added*   10/1999 Update to Version 5.00*   08/2003 Add LOGSCALE and PERIODOGRAM options*   06/2007 Add FOOTER option*   07/2007 Add WINDOW=QUADRATIC*procedure spectrum3 series start endtype series seriestype integer start endoption integer  ordinate                           ;* Number of ordinatesoption choice   window    1  flat  tent quadratic  ;* Window typeoption integer  width                              ;* Window widthoption choice   taper     3  trap cosine none      ;* Taper typeoption real     wtaper                             ;* Taper width (fraction)option switch   graph     1option string   headeroption string   footeroption switch   logscale  1option switch   periodogram 0option series  *spectrumoption series  *freqlocal integer nobs startl endl freqord halflocal series  spect frequenciesif .not.%defined(series) {   DISPLAY "Syntax: @SPECTRUM3 series  start end"   return}**  Determine the range of SERIES to use*inquire(series=series) startl<<start endl<<endcompute nobs=endl-startl+1********************************** <<nwidth>> is the window width to be used in estimation.* <<nwidthboot>> is the window width to be used to standardize the* periodogram for bootstrapping. These can be different.*compute nwidth=9compute nwidthboot=9*********************************** Determine the ordinates. If the user provides a number, use it. If* not, use the number of observations, if the user has requested the* periodogram, or the recommended number, if a more general spectral* density estimate is requested.*if %defined(ordinate)   compute freqord=ordinateelse if periodogram   compute freqord=nobselse   compute freqord=%freqsize(nobs)compute half = freqord/2+1frequency 5 freqord** Transfer the series to the frequency domain - shift it to begin at* entry 1*rtoc startl endl 1# series#   1** If TAPER is requested, taper the transferred part of series 1*if taper<>3   {   taper(type=taper,fraction=wtaper) 1 1 nobs   }else   {   compute %scaletap=nobs   compute %kappa   =1.0   }fft 1cmult(scale=1./(2.*%pi*%scaletap)) 1 1if periodogram   cset 2 = %z(t,1)else   window(type=window,width=width) 1 / 2ctor 1 half#    2#  spectif %defined(spectrum)   set spectrum 1 half = spectif .not.graph   return********************************************************************** Get the smoothed estimate with the width for bootstrapping to complex* series 5.*window(type=tent,width=nwidthboot) 1 / 5compute 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=freqord/2+1** We need to bootstrap the periodogram divided by the estimated* spectral density using the width <<nwidthboot>>.*set source = %real(%z(t,1))/%real(%z(t,5))*set sumsqr 1 freqord = 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 freqord = sumsqr+log(%real(%z(t,4))/%real(%z(t,2)))^2end do drawset lower    1 half = spect-2.0*sqrt(sumsqr/nboot)set upper    1 half = spect+2.0*sqrt(sumsqr/nboot)*******************************************************************set frequencies 1 half = (t-1.0)/halfif logscale {   scatter(style=lines,header=header,footer=footer,hlabel="Fractions of pi",vlog=10.0) 1   # frequencies spect 1 half}else {   scatter(style=lines,header=header,footer="Confidence bands using bootstrap",hlabel="Fractions of pi") 3   # frequencies spect  1 half   # frequencies lower  1 half   # frequencies upper  1 half}if %defined(freq)   set freq 1 half = frequenciesif .not.graph   returnend spectrum3`
rod1

Posts: 5
Joined: Tue Sep 13, 2011 11:53 am

### Re: Bootstrapping confidence bands for spectrum

The mistake that you're making is in your final processing. If you look at what I wrote:

Code: Select all
`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)`

you'll see that the error bands are computed for the log spectrum. Error bands for the spectral density are largely proportional to value, similar to the behavior of the variance in a regression. What is accumulated in the bootstrap loop are the squared log ratios between the bootstrapped values and the original estimates, and you need to allow for that in preparing your graph. If you want the upper and lower bounds in a non-log scale, you'll need to do:

Code: Select all
`set lower    1 half = exp(log(spect)-2.0*sqrt(sumsqr/nboot))set upper    1 half = exp(log(spect)+2.0*sqrt(sumsqr/nboot))`

They'll be asymmetrical (wider on the high side), which is what you want.
TomDoan

Posts: 2725
Joined: Wed Nov 01, 2006 5:36 pm

### Re: Bootstrapping confidence bands for spectrum

It works great, thanks a lot!
rod1

Posts: 5
Joined: Tue Sep 13, 2011 11:53 am