pgram <- function(x, rm.mean=T) #--------------------------------------------------------- # # Compute the periodogram # #--------------------------------------------------------- { if (rm.mean) x <- x - mean(x) per <- Mod(fft(x))^2 / (n <- length(x)) per[1:((n / 2) + 1)] } pgram.test <- function(x,exact=T) #--------------------------------------------------------- # # Fisher's exact and approximated test # #--------------------------------------------------------- { n <- length(x) q <- (n-1) %/% 2 # number of freq's q1 <- q+1 per <- pgram(x) # peridogram a <- q*max(per[2:q1])/sum(per[2:q1]) # Fisher's test if (exact) p <- fisher.p(a,n) # exact P else p <- 1-(1-exp(-a))^q # apprx P list(test=a, pval=p) } fisher.p <- function(a,n) #--------------------------------------------------------- # # Compute the P-value for Fisher's test # Brockwell and Davis (1991, p.339) # #--------------------------------------------------------- { q <- (n-1) %/% 2 bj <- 1 # first binomial coef is 1 pr <- 1 # first term for (j in 1:q) { bj <- bj*(q-j+1)/j if (1-j*a/q <= 0) return (1-pr) else pr <- pr + (-1)^j*bj*(1-j*a/q)^(q-1) } 1-pr } minmax0 <- function(y) #--------------------------------------------------------- # # Find the max and the locations in periodogram # #--------------------------------------------------------- { n <- length(y) yo <- order(y) list(ymax = y[yo[n]], imax = yo[n]) } minmax <- function(y) #--------------------------------------------------------- # # Find the max, min, and the locations # #--------------------------------------------------------- { n <- length(y) yo <- order(y) ymin.f <- (yo[1]-1)/(2*n) ymax.f <- (yo[n]-1)/(2*n) list(ymin = y[yo[1]], ymin.f = ymin.f, ymax = y[yo[n]], ymax.f = ymax.f) } auspe <- function(x, ioptw=4, M=15, logs=F, ci=T) #--------------------------------------------------------- # # Function to estimate and plot the auto spectral density # of a time series x using window ioptw and M # # Arguments: # x1, x2: time series # ioptw: window option 1 ... 8 # M: scale parameter for window # logs: flag for logarithmic scale # ci: confidence intervals # # Value: None # #--------------------------------------------------------- { x.spe <- spe(x, ioptw, M) f11 <- x.spe$f crv <- x.spe$c n <- length(f11) if ( logs ) { f11 <- log(f11) fl <- f11-crv fu <- f11+crv if (ci) ci.plot(f11,fl,fu,y.lab="Log-spectrum") else sp.plot(f11,y.lab="Log-spectrum") } else { crv <- exp(crv) fl <- f11/crv fu <- f11*crv if (ci) ci.plot(f11,fl,fu,y.lab="Spectrum") else sp.plot(f11,y.lab="Spectrum") } } sp.plot <- function(f, y.lab="Spectrum", p.type="l") # ----------------------------------------------------- # # Function for plotting spectral quantities # # ----------------------------------------------------- { q <- length(f) plot(c(0:(q-1))/(2*(q-1)),f,xlab="frequency", ylab=y.lab,type=p.type) } ci.plot <- function(f,fl,fu,y.lab="Spectrum") # ----------------------------------------------------- # # Function for plotting spectral quantities and CI # # ----------------------------------------------------- { n <- length(f) plot(c(0,.5), c(min(fl), max(fu)), xlab = "frequency", ylab = y.lab, type = "n") lines(c(0:(n - 1))/(2 * (n - 1)), f) lines(c(0:(n - 1))/(2 * (n - 1)), fl, lty=2) lines(c(0:(n - 1))/(2 * (n - 1)), fu, lty=2) } coher <- function(x1,x2,ioptw,M) #--------------------------------------------------------- # # Function to estimate the spectra, coherence and phase # of bivariate series x1 and x2 using window ioptw and M # # Arguments: # x1, x2: time series # ioptw: window option 1 ... 8 # M: scale parameter for window # # Value: None # #--------------------------------------------------------- { f11 <- spe(x1,ioptw,M)$f f22 <- spe(x2,ioptw,M)$f f12 <- xspe(x1,x2,ioptw,M) c12 <- f12$re q12 <- f12$im amp <- polar(c12,q12)$amp phi <- polar(c12,q12)$phase coh <- amp^2/(f11*f22) par(mfrow=c(2,2)) sp.plot(f11, y.lab = "Spectrum f11") sp.plot(f22, y.lab = "Spectrum f22") ilam2 <- c(2,2./3.,.795,.539,.586,1,1.2,1.66) n <- length(x1) w12 <- sqrt(coh) se <- sqrt(abs(1-coh)*M*ilam2[ioptw]/(2*n)) lo <- w12-2*se up <- w12+2*se ci.plot(w12, lo, up, y.lab="Coherence") se <- amp*sqrt(abs((1/coh)-1)*M*ilam2[ioptw]/(2*n)) lo <- phi-2*se up <- phi+2*se ci.plot(phi, lo, up, y.lab="Phase") par(mfrow=c(1,1)) invisible() } cohtest <- function(x1,x2,ioptw,M) #--------------------------------------------------------- # # Function to estimate the spectra, coherence and phase # of bivariate series x1 and x2 using window ioptw and M # # Arguments: # x1, x2: time series # ioptw: window option 1 ... 8 # M: scale parameter for window # # Value: None # #--------------------------------------------------------- { f11 <- spe(x1,ioptw,M)$f f22 <- spe(x2,ioptw,M)$f f12 <- xspe(x1,x2,ioptw,M) c12 <- f12$re q12 <- f12$im amp <- polar(c12,q12)$amp phi <- polar(c12,q12)$phase coh <- amp^2/(f11*f22) n <- length(f11) par(mfrow=c(2,2)) sp.plot(f11, y.lab = "Spectrum f11") sp.plot(f22, y.lab = "Spectrum f22") ilam2 <- c(2,2./3.,.795,.539,.586,1,1.2,1.66) const <- (n - M*ilam2[ioptw])/(M*ilam2[ioptw]) fstat <- const*coh/(1-coh) nc <- 2*const ft <- qf(.95,2,nc) sp.plot(fstat, y.lab = "F-stat") lines(c(0:(n - 1))/(2 * (n - 1)), rep(ft,n),lty=2) sp.plot(phi, y.lab = "Phase") par(mfrow=c(1,1)) invisible() } crossp <- function(x1,x2,ioptw,M, ci=T) #--------------------------------------------------------- # # Function to estimate, plot the co and quad spectra # of bivariate series x1 and x2 using window ioptw and M # # Arguments: # x1, x2: time series # ioptw: window option 1 ... 8 # M: scale parameter for window # # Value: None # # TIMESLAB (p.338) #--------------------------------------------------------- { f11 <- spe(x1,ioptw,M)$f f22 <- spe(x2,ioptw,M)$f f12 <- xspe(x1,x2,ioptw,M) crv <- f12$c c12 <- f12$re q12 <- -f12$im # Compute the avar v1 <- f11*f22 v2 <- c12^2-q12^2 sd <- matrix(c(sqrt((v1+v2)/2),sqrt((v1-v2)/2)),ncol=2) f <- matrix(c(c12,q12),ncol=2) n <- length(f[,1]) par(mfrow=c(1,2)) for (j in (1:2)) { s <- f[,j] se <- sd[,j]*crv cil <- s - se ciu <- s + se if (j==1) { y.lab <- "Co-spectrum" } else { y.lab <- "Quad-spectrum" } if (ci) ci.plot(s, cil, ciu, y.lab = y.lab) else sp.plot(s,y.lab=y.lab) } par(mfrow=c(1,1)) invisible() } divpoly <- function(a,b,n) # ------------------------------------------------ # # Find the first n coefficients of # (1+sum_1^q b(j)z^j) / (1+sum_1^p a(i)z^i) # # ------------------------------------------------ { p <- length(a) q <- length(b) if (p+q == 0) return(rep(0,p+q)) g <- c(b,rep(0,n-q)) if (p > 0) { for (i in 1:n) { c0 <- g[i] for (j in 1:min(i,p)) { c1 <- 1 if (i-j>0) c1 <- g[i-j] c0 <- c0 - a[j]*c1 } g[i] <- c0 } } g } estfilt <- function(x1,x2,ioptw,M,kn) #--------------------------------------------------------- # # Function to estimate the best filter coefficients # of bivariate series x1 and x2 using window ioptw and M. # It uses the inverse Fourier transform on A=f21/f11. # # Arguments: # x1, x2: time series # ioptw: window option 1 ... 8 # M: scale parameter for window # kn+1: number of beta's # # Value: beta[1:(kn+1)] (coefficients of the transfer function) # # See Theorem 4.1.6 of TIMESLAB (p.321) # #--------------------------------------------------------- { f11 <- spe(x1,ioptw,M)$f f22 <- spe(x2,ioptw,M)$f f21 <- xspe(x2,x1,ioptw,M) c21 <- f21$re q21 <- f21$im c211 <- c21/f11 q211 <- q21/f11 q <- length(c21) c211 <- c(c211,c211[(q-1):2]) # extend to [0,1] via sym q211 <- c(q211,-q211[(q-1):2]) # extend to [0,1] via antisym z <- c211 + q211*1i # Form complex array for FFT beta <- Re(fft(z,inv=T))/(2*q-2) beta[1:(kn+1)] } filtest <- function(x1,x2,ioptw,M,kn) #--------------------------------------------------------- # # Function to estimate the transfer function of a filter # of bivariate series x1 and x2 using estfilt and prewhitening # # Arguments: # x1, x2: time series # ioptw: window option 1 ... 8 # M: scale parameter for window # kn+1: number of beta's # # Value: coef (coefficients of the transfer function) # # See TIMESLAB (p.344) # #--------------------------------------------------------- { z <- prewhite(x1,x2) x1 <- z$x x2 <- z$y ax <- -z$ax ay <- -z$ay beta <- estfilt(x1,x2,ioptw,M,kn) beta0 <- beta[1] beta <- beta/beta0 kn1 <- kn+1 gama <- multpoly(beta[2:kn1],ax) betay <- divpoly(a=ay,b=c(),40) gama <- multpoly(betay,gama) gama <- gama*beta0 list(coef=gama[1:kn], b0=beta0) } indtest <- function(x1, x2, nlag=20, alpha=.05) # ------------------------------------------------ # # To test the independence of two time series # using cross-correlations # # Argument: x1, x2, nlag, alpha # # Value: upper confidence limit # # ------------------------------------------------ { z <- prewhite(x1,x2) zz <- matrix(c(z$x,z$y),ncol=2) zacf <- acf(zz, lag.max=nlag) # plot acf n <- length(z$x) # for CI qnorm( (1+(1-alpha)^(1/(2*nlag+1)))/2) / (n^.5) # upper CI } multpoly <- function(a,b) # ------------------------------------------------ # # Multiply 1+sum_1^p a(i)z^i and 1+sum_1^q b(j)z^j # # ------------------------------------------------ { p <- length(a) q <- length(b) g <- rep(0,p+q) g <- g + c(a,rep(0,q)) g <- g + c(b,rep(0,p)) for (i in 1:p) { for (j in 1:q) g[i+j] <- g[i+j] + a[i]*b[j] } g } polar <- function(zr,zi) #---------------------------------------------------- # # Function for converting to polar representation # #---------------------------------------------------- { return(list(amp=sqrt(zr^2+zi^2), phase=atan(zi,zr))) } prewhite <- function(x,y) # ----------------------------------------------------- # # Use AR filters to prewhiten and align two time series # # ----------------------------------------------------- { n <- length(x) x <- x - mean(x) y <- y - mean(y) x.ar <- ar(x) y.ar <- ar(y) px <- x.ar$order py <- y.ar$order if (px == 0) xx <- x if (px > 0) { xx <- x.ar$resid xx <- xx[(px+1):n] } if (py == 0) yy <- y if (py > 0) { yy <- y.ar$resid yy <- yy[(py+1):n] } nx <- n - px ny <- n - py nn <- nx if (nx > ny) { nn <- ny n1 <- nx-ny n1 <- n1 + 1 xx <- xx[n1:nx] } if (nx < ny) { nn <- nx n1 <- ny-nx n1 <- n1 + 1 yy <- yy[n1:ny] } list(n=nn, x=xx, y=yy, ax=x.ar$ar, ay=y.ar$ar) } spe <- function(x, ioptw, M, alpha=.05) #----------------------------------------------------------- # # Function to find auto-spectra using scale parameter M, # window number ioptw: # # 1 Truncated # 2 Bartlett # 3 Tukey # 4 Parzen # 5 Bohman # 6 Daniell # 7 Bartlett-Priestley # 8 Parzen-Cogburn-Davis # # The spectral estimator and the constant c used in confidence # intervals are returned in a list. # # TIMESLAB (p.185) # #------------------------------------------------------------- { n <- length(x) if(M >= n || ioptw < 1 || ioptw > 8 || alpha <= 0. || alpha >= 1.) stop("Illegal Input to spe()") x.acf <- acf(x,lag.max=n-1,type="cov",plot=F) z <- x.acf$acf[,1,1] # c11 z1 <- rep(0,n) if(ioptw == 1) z1[1:(M+1)] <- rep(1, M+1) else if(ioptw == 2) z1[1:(M+1)] <- 1 - (c(0:M) / M) else if(ioptw == 3) z1[1:(M+1)] <- .54+.46*cos(pi*c(0:M)/M) else if(ioptw == 4) { u <- c(1:M)/M u1 <- u[u<=.5] u2 <- u[u>.5] z1[1:(M+1)] <- c(1,1-6*u1^2+6*u1^3,2*(1-u2)^3) } else if(ioptw == 5) { u <- c(0:M) / M piu <- pi * u z1[1:(M+1)] <- (1 - u) * cos(piu) + sin(piu) / pi } else if(ioptw == 6) { piu <- pi * (c(1:(n-1)) / M) z1[1:n] <- c(1,sin(piu)/piu) } else if(ioptw == 7) { piu <- pi * (c(1:(n-1)) / M) z1[1:n] <- c(1,3*((sin(piu)/piu) - cos(piu)) / (piu*piu)) } else if(ioptw == 8) { u <- c(0:(n-1))/M z1[1:n] <- 1 /(1+u^4) } z <- c(z1[1] * z[1], 2 * z1[2:n] * z[2:n]) z <- Re(fft(z)) z <- z[1:((n / 2) + 1)] ilam2 <- c(2,2./3.,.795,.539,.586,1,1.2,1.66) fac <- qnorm(1-(alpha/2)) * sqrt(M*ilam2[ioptw]/n) return(list(f=z,c=fac)) } xspe <- function(x1, x2, ioptw, M, alpha=.05) #------------------------------------------------------------ # # Function to find cross-spectra using scale parameter M, # window number ioptw: # # 1 Truncated # 2 Bartlett # 3 Tukey # 4 Parzen # 5 Bohman # 6 Daniell # 7 Bartlett-Priestley # 8 Parzen-Cogburn-Davis # # The spectral estimator and the constant c used in confidence # intervals are returned in a list. # # TIMESLAB (p.335) # #------------------------------------------------------------ { if (length(x1) != length(x2)) stop("Series with unequal length!") n <- length(x1) if(M >= n || ioptw < 1 || ioptw > 8 || alpha <= 0. || alpha >= 1.) stop("Illegal Input to xspe()") x <- matrix(c(x1,x2),n,2) x.acf <- acf(x,lag.max=n-1,type="cov",plot=F) c12 <- x.acf$acf[,1,2] c21 <- x.acf$acf[,2,1] z1 <- rep(0,n) if(ioptw == 1) z1[1:(M+1)] <- rep(1, M+1) else if(ioptw == 2) z1[1:(M+1)] <- 1 - (c(0:M) / M) else if(ioptw == 3) z1[1:(M+1)] <- .54+.46*cos(pi*c(0:M)/M) else if(ioptw == 4) { u <- c(1:M)/M u1 <- u[u<=.5] u2 <- u[u>.5] z1[1:(M+1)] <- c(1,1-6*u1^2+6*u1^3,2*(1-u2)^3) } else if(ioptw == 5) { u <- c(0:M) / M piu <- pi * u z1[1:(M+1)] <- (1 - u) * cos(piu) + sin(piu) / pi } else if(ioptw == 6) { piu <- pi * (c(1:(n-1)) / M) z1[1:n] <- c(1,sin(piu)/piu) } else if(ioptw == 7) { piu <- pi * (c(1:(n-1)) / M) z1[1:n] <- c(1,3*((sin(piu)/piu) - cos(piu)) / (piu*piu)) } else if(ioptw == 8) { u <- c(0:(n-1))/M z1[1:n] <- 1 /(1+u^4) } z <- fft(z1[1:n]*c12[1:n])+Conj(fft(c(0,z1[2:n]*c21[2:n]))) z <- z[1:((n / 2) + 1)] ilam2 <- c(2,2./3.,.795,.539,.586,1,1.2,1.66) fac <- qnorm(1-(alpha/2)) * sqrt(M*ilam2[ioptw]/n) return(list(re=Re(z), im=Im(z), c=fac)) }