# Example R code to calculate the SPECTRAL ENVELOPE for a continuous time series # *** Must source mvspec.R first *** # u=scan("gnpres.dat") # data (use the residuals from example 3.35 on gnp) x=cbind(u, abs(u), u^2) # possible transformations (identity, absolute value, square) # Var = var(x) # var-cov matrix xspec = mvspec(x, spans=c(7,7)) # spectral matrices are an array called fxx fxxr = Re(xspec\$fxx) # fxxr is real(fxx) # #----- compute Q = Var^-1/2 -----# ev = eigen(Var) Q = ev\$vectors%*%diag(1/sqrt(ev\$values))%*%t(ev\$vectors) # #--- compute spec env and scale vectors ---# num = xspec\$n.used # sample size used for FFT nfreq = length(xspec\$freq) # number of freqs used specenv = matrix(0,nfreq,1) # initialize the spec envelope beta=matrix(0,nfreq,3) # initialize the scale vectors for (k in 1:nfreq){ ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num) # get evalues of normalized spectral matrix at freq k/n specenv[k] = ev\$values[1] # spec env at freq k/n is max evalue b = Q%*%ev\$vectors[,1] # beta at freq k/n beta[k,] = b/b[1] # first coef is always 1 } # #--- output and graphics ---# frequency = (0:(nfreq-1))/num plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)") ## add significance threshold to plot ## m=xspec\$kernel\$m etainv=sqrt(sum(xspec\$kernel[-m:m]^2)) thresh=100*(2/num)*exp(qnorm(.9999)*etainv)*matrix(1,nfreq,1) lines(frequency,thresh, lty="dashed", col="blue") #-- details --# output = cbind(frequency, specenv, beta) colnames(output)=c("freq","specenv","x", "|x|", "x^2") # output # uncomment to write to screen # write.table(output, file="output.txt", quote=FALSE, row.names=FALSE) # uncomment to write to file