spec.env <- function(x, m) { #computes spectral envelope, takes multivariate or univariate sequence as input #x (should be a matrix nxk) #input m is the desired degree of smoothness of the periodogram n <- dim(x)[1] k <- dim(x)[2] x <- x[1:n,] dft <- apply(x,2,fft)/sqrt(n) dft <- dft[-1,] ift <- array(dim=c(k,k,as.integer(n/2))) for (ii in 1:(n/2)) { ift[,,ii] <- Re(matrix(dft[ii,],nrow=k,ncol=1) %*% Conj(matrix(dft[ii,],nrow=1,ncol=k))) } fw <- trismooth(ift, as.integer(n/2), m, k)$fw lambdaw <- matrix(nrow=dim(fw)[3], ncol=1) betaw <- matrix(nrow=dim(fw)[3], ncol=k) #S <- var(x)*(n-1)/n S <- (t(x-colMeans(x)) %*% (x-colMeans(x)))/n S.eig <- eigen(S) Sroot <- S.eig$vectors %*% diag(1/sqrt(S.eig$values)) %*% t(S.eig$vectors) for (ii in 1:(n/2-2*m)) { deq <- 2/n* Sroot %*% matrix(fw[,,ii],nrow=k,ncol=k) %*% Sroot deq.eig <- eigen(deq) lambdaw[ii] <- deq.eig$values[1] betaw[ii,] <- matrix(deq.eig$vectors[,1],ncol=k,nrow=1) %*% Sroot } j <- seq((m+1):(as.integer(n/2)-m)) omegaj <- (j+m)/n plot(omegaj, lambdaw, type="l") list(SpecEnvelope = lambdaw, scale = betaw, freq=omegaj, m=m) } trismooth <- function(ift, n, m, k){ #symmetric moving average smoothing with triangular weighting fw <- array(dim = c(k, k, (n - 2 * m))) for(ii in (1 + m):(n - m)) { fw[, , (ii - m)] <- (m + 1)/(m + 1)^2 * matrix(ift[, , ii], nrow = k, ncol = k) if(m >= 1) { for(jj in 1:m) fw[, , (ii - m)] <- fw[, , (ii - m)] + jj/(m + 1)^2 * matrix(ift[, , (ii - m + jj - 1)], nrow = k, ncol= k) for(kk in 1:m) fw[, , (ii - m)] <- fw[, , (ii - m)] + kk/(m + 1)^2 * matrix(ift[, , (ii + m - kk + 1)], nrow = k, ncol= k) } } list(fw = fw) }