Re: How to calculate cofidence intervals for cospectrum?

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
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk