## cosine bell
cbell <- function(x, p=.2)
{
    ## MM: fast version of ifelse()  {not working with NAs !}
    r <- x
    if(any(i1 <- x  <=   p/2)) r[i1] <- 1/2*(1-cos(2*pi* x[i1] / p ))
    if(any(i3 <- x > 1 - p/2)) r[i3] <- 1/2*(1-cos(2*pi*(1-x[i3])/p))
    if(any(i2 <- !(i1 | i3)))  r[i2] <- 1
    r
}

approximations <- function(N0,N,gg,ht)
{
  ## Purpose: Calculates new and usual version of the approximate variance
  ## ----------------------------------------------------------------------
  ## Arguments: N0 : length of times series
  ##            N  : number of grid points for smoothing
  ##            gg : vector containing smoothing weights
  ##            ht : vector of length N0 containing the data taper
  ## ----------------------------------------------------------------------
  ## Author: Michael Amrein, Date: 17 Jan 2011, 16:46
  tapsq <- ht^2 ## squared taper
  tapsqpad <-c(tapsq,rep(0,N-N0)) ## padded squared taper
  fftt <- fft(tapsqpad) ## fourier-transform
  ffttsq <- fftt[1:(2*M+1)]
  ffttsq <- Re(ffttsq*Conj(ffttsq)/(sum(tapsq))^2)
  ## Re is used to get rid of the 0i-terms
  CC <- toeplitz(ffttsq) ## "covariance" matrix
  
  ## approximated variance of Var(S^ds(f)/S(f)) (more accurate version)
  ## quadratic form
  varSds.mav <- c(t(gg) %*% (CC %*% gg))
  
  ## approximated variance of Var(S^ds(f)/S(f)) (old version)
  varSds.ov <- sum(gg^2)*N*(sum(tapsq^2))/((sum(tapsq))^2)

  c(varSds.mav,varSds.ov)
}


exactgauss <- function(arcoef,N0,gg,ht)
{
  ## Purpose: Calculates exact relative variances for Gaussian AR processes
  ## ----------------------------------------------------------------------
  ## Arguments: arcoef: coefficients of AR-model
  ##            N0, N, gg, ht: see function approximations
  ## ----------------------------------------------------------------------
  ## Author: Michael Amrein, Date: 17 Jan 2011, 17:00

  ## packages needed
  stopifnot(require(FitAR),
            require(TSA))

  LAG <- TacvfAR(phi=arcoef, lag.max = N0-1)
  R <- toeplitz(LAG) # = Autocovariances  s_{j-k}
  
  ## t in formula at the beginning of section 5
  tt <- matrix(rep(1:N0,N0),N0,N0)
  ## k in formula at the beginning of section 5
  kk <- matrix(rep(1:N0,N0),N0,N0,byrow = TRUE) 
  hh <- outer(ht,ht) ## contains h_j h_k

  
  ## matrix of covariances of specral estimates at two frequencies
  xx <- 0:(N-1)
  C <- matrix(NA,length(xx),length(xx))
  
  ## true spectrum at frequencies f/N
  DENS <- ARMAspec(model=list(ar=arcoef),freq=xx/N,plot=FALSE)
  SS <- as.vector(DENS$spec) 

  nu <- -2*pi*1i/N
  hR <- hh*R

  ## calculation of true covariance of the estimates at frequencies (k-v)/N
  ## for -M<=v<= M and (k-s)/N for -M<=k<= M for 0<=k<=N-1 needed in the exact
  ## formula for the variance (6) according to the formula at the beginning of
  ## section 5 whithout normalization (sum(ht^2))^2
  
  ## for (l in xx) {
  ##   ## The following double loop is  O( M^2 ) ---
  ##   ## but should be writable as one loop in  O ( M )  [M.A.]
  ##   for (u in l+(-M:M)) {
  ##     for (v in l+(-M:M)) {
  ##       c1 <- sum(hR * exp(nu*(u*tt-v*kk)))
  ##       c2 <- sum(hR * exp(nu*(u*tt+v*kk)))
  ##       C[u%%N+1,v%%N+1] <- (Re(c1)^2+Im(c1)^2+Re(c2)^2+Im(c2)^2)/SS[u%%N+1]/SS[v%%N+1]
  ##     }
  ##   }
  ##   if(l%%10==0){
  ##     cat(l)
  ##     cat(",")
  ##   }
  ## }
  
  for (l in xx) {
    for (u in l+((-2*M):(2*M))) {
      c1 <- sum(hR * exp(nu*(u*tt-l*kk)))
      c2 <- sum(hR * exp(nu*(u*tt+l*kk)))
      C[u%%N+1,l%%N+1] <- (Re(c1)^2+Im(c1)^2+Re(c2)^2+Im(c2)^2)/SS[u%%N+1]/SS[l%%N+1]
    }
    if(l%%10==0){
      cat(l)
      cat(",")
    }
  }

  ## calculation of the exact variances according to formula (6)
  x <- 0:ceiling(N/2)
  varx <- rep(NA,length(x))
  for (ix in x){
    varx[ix+1] <- t(gg)%*%C[(ix+(-M:M))%%N+1,(ix+(-M:M))%%N+1,drop=FALSE]%*%gg
  }
  cbind(x/N,varx/(sum(ht^2))^2)
}

variances.pic<- function(AP,EX)
{
  ## Purpose: Plotting exact variances and approximations
  ## ----------------------------------------------------------------------
  ## Arguments: AP: "c(new approximation, old approximation)"
  ##            EX: matric containing frequencies in 1st column and
  ##                exact variances in the 2nd column
  ## ----------------------------------------------------------------------
  ## Author: Michael Amrein, Date: 17 Jan 2011, 17:47
  plot(EX[,1],EX[,2],type="l",lty=1,col="green",lwd=2,xlab="f",ylab="relative variance",
       ylim=c(min(EX[,2]),max(max(EX[,2]),AP[2])))
  abline(h=AP[1],lwd=1,lty=1)
  abline(h=AP[2],lwd=1,lty=2)
}


## Test

## AR-model
arcoef <- c(2.7607,-3.8106,2.6535,-0.9238)
## no of obs
N0 <- 1024
## number of grid-frequencies
N <- 2*N0 
## weights
M <- 1
gg <- rep(1/(2*M+1),2*M+1)
## taper
ht <- cbell((2*(1:N0)-1)/(2*N0),p=.2)
## approximations
AP <- approximations(N0,N,gg,ht) # 0.5631923 0.7442177
## exact (takes some time ...)
EX <- exactgauss(arcoef,N0,gg,ht)
## picture
variances.pic(AP,EX)