Re: How to calculate cofidence intervals for cospectrum?

From: Hyacinth Nnamchi <hyacinth.1_at_nyahnyahspammersnyahnyah>
Date: Sun Jul 15 2012 - 11:31:15 MDT

Thanks Dennis for the reply and explanations. I'll look for alternative means, and report back if successful.
Btw, spec@cospc does not run. Says spcx,xlag1,frq,spcx are undefined even when I manually specified them. And if it does, I suppose it'd give the wrong result: it is only of length N/2 but an entire series is needed to calculate the confidence intervals.

Thanks again,

Hyacinth

> ------------------------------
>
> Message: 2
> Date: Sat, 14 Jul 2012 06:52:32 -0600
> From: Dennis Shea <shea@ucar.edu>
> Subject: Re: How to calculate cofidence intervals for
> cospectrum?
> To: Hyacinth Nnamchi <hyacinth.1@hotmail.com>
> Cc: ncl-talk@ucar.edu
> Message-ID: <50016B90.1040105@ucar.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> ncl-talk is a language forum and questions on the NCL
> language are always answered. Questions on other subjects
> like statistics or physics may be posted. However, answers
> are not guaranteed. Perhaps, people are not confident
> about answering, particularly in a public forum.
>
> ---
>
> Generally, a 'geophysical' time series, x(t), has some degree
> of auto-correlation. Hence, the spectrum of x(t), SX(f), is often
> tested against a 'red noise' model. The red noise spectrum has more
> variance (power) at lower frequencies. This is discussed in numerous
> statistical text books. eg: Jenkins and Watts, (1968), Spectral
> Analysis and Its Applications
>
> http://www.ncl.ucar.edu/Applications/spec.shtml
> Example 3
>
> A less 'powerful' test would be to test against a 'white-noise'
> spectrum. This is one in which the variance is equally spread
> throughout the spectrum.
>
> The 'white-noise' approach could also be used to test the co-spectrum.
> In this case, the covariance (escovc) is equally spread throughout the
> cospectrum.
>
> For a cospectrum derived from two autocorrelated series
> [x(t),y(t) ==> SXY(f)], the cross correlation
> could be negative which would yield a 'blue spectrum'
> (more power at high frequencies). *Perhaps* the same function
> used for the univariate power spectrum could be used for the
> cospectrum! Use spec@cospc rather than spec@specx
> ****Don't know the answer****
>
> Please report back to ncl-talk if you find an appropriate test.
>
> ---
> You can see the code used for the red noise process at
>
> %> less $NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl
>
> search for: specx_ci
>
> This uses the 'chiinv' function to get the upper and lower limis
> of a red noise spectrum.
>
>
> On 7/13/12 10:45 AM, Hyacinth Nnamchi wrote:
> > Dear all;
> >
> > I didn't get any response on this and have not got it right myself. The
> > question is: How do I calculate the confidence intervals for co-spectrum
> > i.e.
> >
> > spec@cospc?
> >
> > Thanks,
> >
> > Hyacinth
> >
> >
> >
> > ------------------------------------------------------------------------
> > From: hyacinth.1@hotmail.com
> > To: ncl-talk@ucar.edu
> > Subject: Cospectrum: How about specxy_ci?
> > Date: Fri, 29 Jun 2012 22:29:32 +0800
> >
> > Hi;
> >
> > The specx_ci functiom yields the confidence intervals for one time
> > series. However, I have calculated the cospectrum line of 2 time series
> > and would like estimate the Markov confidence intervals. But, there is
> > no specxy_ci function (?). Does anyone know how I can estimate the
> > Markov confidence intervals of cospectrum plot? Relevant parts of the
> > code and the created plot are attached.
> >
> > Thanks in advance,
> >
> > Hyacinth
> >
> > ;************************************************
> > ; set function arguments for spectrum analysis
> > ;************************************************
> > d = 0 ; detrending opt: 0=>remove mean 1=>remove mean + detrend
> > sm = 5 ; smooth: should be at least 3 and odd
> > pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common.
> > ;************************************************
> > ;
> > ; Co-spectrum
> > spec = specxy_anal(NEP,SWP,d,sm,pct)*splt = specx_ci(spec,0.05,0.95)
> >
> > ;************************************************
> > ; plotting parameters
> > ;************************************************
> > wks = gsn_open_wks("ps","specxy_ci") ; Opens a ps file
> > plot = new(4,graphic) ; create graphic array
> >
> > r = True ; plot mods desired
> > r@gsnDraw = False &nbs p; ; do not draw
> > r@gsnFrame = False ; do not advance frame
> > r@trYLog = True ; log scale
> > r@trXLog = True
> > r@trXReverse = True ; Reverse X axis values
> > r@trXMinF = 2. ; manually set lower limit
> > r@trXMaxF = 112. &n bsp; ; " upper
> > r@tmXBMode = "Explicit"
> > r@tmXBValues = (/"100","50","20","10","5","4","3","2","1"/)
> > r@tmXBLabels = (/"100","50","20","10","5","4","3","2","1"/)
> > r@tmXBMinorValues = ispan(1,20,1)
> > r@xyLineThicknesses = (/4/) ; Define line thicknesses
> > r@xyDashPatterns = (/1/) ; Dash patterns
> > r@xyLineColors = (/"blue"/)
> > r@tiMainString = "" ; title
> > r@tiXAxisString = "" ; xaxis
> > r@tiYAxisString = "" &n bsp; ; yaxis
> >
> >
> > ; ************************************************
> > ; create plot of cospectrum
> > ;************************************************
> > r@trYMinF = -.01 ; manually set lower limit
> > r@trYMaxF = 5. ; " upper
> >
> > r@tiYAxisString = "Cospectrum: Confidence intervals (CI)?" ; yaxis
> > plot(0)=gsn_csm_xy(wks,1/(spec@frq),spec@cospc,r); create plot
> > &n bsp;
> >
> > r@tiYAxisString = "NEP spectrum & CI"
> > r@trYMinF = .1 ; manually set lower limit
> > r@trYMaxF = 10. ; " upper
> > plot(1)=gsn_csm_xy(wks,1/(spec@frq),splt,r) ;plot the specx create plot
> >
> >
> > r@tiYAxisString = "Spectrum of NEP" ; yaxis
> > plot(2)=gsn_csm_xy(wks,1/(spec@frq),spec@spcx,r) ; create plot&nb sp;
> >
> > r@tiYAxisString = "Spectrum of SWP" ; yaxis
> > plot(3)=gsn_csm_xy(wks,1/(spec@frq),spec@spcy,r) ; create plot
> > ;==============================================================
> > ; Panel the plots created
> > ;
> > ;==============================================================
> > resP = True ; modify the panel plot
> > resP@gsnPanelYWhiteSpacePercent = 0.0 ; default is 1.0
> > ; resP@gsnPanelFigureStrings= (/"LTREND EXPT(raw)","LTREND
> > EXPT(dtrend)"/) ; add strings to panel
> > resP@amJust = "BottomLeft"
> > resP@gsnMaximize = True ; large format
> > ; resP@gsnPaperOrientation = "portrait"
> > gsn_panel(wks,plot,(/2,2/),resP)
> > *
> > **
> > *
> > *
> >
> > *
> > *
> >
> > *_______________________________________________
> > ncl-talk mailing list
> > List instructions, subscriber options, unsubscribe:
> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk*
>
>
> ------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk@ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
> End of ncl-talk Digest, Vol 104, Issue 32
> *****************************************

                                               

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Jul 15 11:31:28 2012

This archive was generated by hypermail 2.1.8 : Wed Jul 18 2012 - 14:33:00 MDT