Bootstrapping confidence bands for spectrum

[rod1]

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.

[TomDoan]

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

(2.63 KiB) Downloaded 82 times

sunspots.dat

Data file

(393 Bytes) Downloaded 85 times

[rod1]

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

[ivory4]

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.

[TomDoan]

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)

[rod1]

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.

[TomDoan]

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.

[rod1]

Hello. Thanks for your reply.

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 end
type series series
type integer start end

option integer  ordinate                           ;* Number of ordinates
option choice   window    1  flat  tent quadratic  ;* Window type
option integer  width                              ;* Window width
option choice   taper     3  trap cosine none      ;* Taper type
option real     wtaper                             ;* Taper width (fraction)
option switch   graph     1
option string   header
option string   footer
option switch   logscale  1
option switch   periodogram 0
option series  *spectrum
option series  *freq

local integer nobs startl endl freqord half
local series  spect frequencies

if .not.%defined(series) {
   DISPLAY "Syntax: @SPECTRUM3 series  start end"
   return
}

*
*  Determine the range of SERIES to use
*

inquire(series=series) startl<<start endl<<end
compute 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=9
compute 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=ordinate
else if periodogram
   compute freqord=nobs
else
   compute freqord=%freqsize(nobs)

compute half = freqord/2+1
frequency 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 1
cmult(scale=1./(2.*%pi*%scaletap)) 1 1
if periodogram
   cset 2 = %z(t,1)
else
   window(type=window,width=width) 1 / 2



ctor 1 half
#    2
#  spect

if %defined(spectrum)
   set spectrum 1 half = spect

if .not.graph
   return



*********************************************************************
* 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=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)))^2
end do draw



set 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)/half
if 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 = frequencies

if .not.graph
   return


end spectrum3


[TomDoan]

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.

[rod1]

It works great, thanks a lot!