# R code to calculate the SPECTRAL ENVELOPE for a categorical time series 
# *** Must source mvspec.R first ***           
# The data set is a column of INTEGERS representing the categories 



u = factor(scan("/mydata/bnrf1ebv.dat"))   # first, input the data as factors and then 
x = model.matrix(~u-1)[,1:3]               # make an indicator matrix
#x = x[1:1000,]                  # select subsequence if desired
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/sqrt(sum(b^2))          # helps to normalize beta
  }
#
#--- 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","A", "C", "G")
# output                                                                # uncomment to write to screen
# write.table(output, file="output.txt", quote=FALSE, row.names=FALSE)  # uncomment to write to file