###Define the tukey window and constants postscript("test.ps") ###first compute the contants auxwindow <- function(x)(1-x^3)^3 TT <- 5001 N <- (TT-1)/2 waux <- apply(as.array(c(N:0,1:N)/N),1,auxwindow) aux <- c(1:TT)/TT - .5/TT W0 <- mean(waux) W1 <- mean(aux*waux) W2 <- mean(aux^2*waux) U0 <- mean(waux^2) U1 <- mean(aux*waux^2) U2 <- mean(aux^2*waux^2) a0 <- (W0*W2-W1^2)^{-2} a1 <- (U0*U2 - U1^2) a2 <- W0^(-2)*(W0*U1 - W1*U0)^2 b0 <- (U2*W0^2 - 2*W1*U1*W0+U0*W1^2) b1 <- (U2*W1^2 - 2*W2*U1*W1+U0*W2^2) C0 <- a0*b0 C1 <- U0/W0^2 C2 <- a0*b1 C3 <- a0*W1/(W0^2)*(W0^2*W1*U2 - W1^3*U0 - 2*W0^2*W2*U1 + 2*W0*W1*W2*U0) C4 <- a0*(W0*W1*U2 + W1*W2*U0 - U1*(W1^2 + W0*W2)) rm(a0,a1,a2,b0,b1,auxwindow) Window <- function(n,h){ aux <- seq(-1,1,len=round(n*h)) Naux1 <- round(n- length(aux))/2 Naux2 <- n- length(aux) - Naux1 c(rep(0,Naux1),(1-abs(aux)^3)^3,rep(0,Naux2)) } ##now we run the program that computes all the variances cycles <- c(1:30) par(mfrow=c(1,2)) plot(cycles,rep(0,length(cycles)),ylab="percent difference",main="Comparison of asymptotic vars\nof fund freq",ylim=c(0,100),type="n",las=1) plot(cycles,rep(0,length(cycles)),ylab="percent difference",main="Comparison of asymptotic vars\nof total amp",ylim=c(0,100),type="n",las=1) Ns <- c(50,100,500,1000) A <- 2 B<-1 sigma2 <- .01 for(j in 1:4){ N <- Ns[j] tt <- 1:N w <- Window(N,1) vars <- sapply(cycles,function(i){ l<- i*pi/N ltt <- l*tt J1<-matrix(c(A^2*sum(w*tt^2*sin(ltt)^2)+B^2*sum(w*tt^2*cos(ltt)^2) -2*A*B*sum(w*tt^2*cos(ltt)*sin(ltt)), -A*sum(w*tt*cos(ltt)*sin(ltt))+B*sum(w*tt*cos(ltt)^2), -A*sum(w*tt*sin(ltt)^2)+B*sum(w*tt*cos(ltt)*sin(ltt)), -A*sum(w*tt*cos(ltt)*sin(ltt))+B*sum(w*tt*cos(ltt)^2), sum(w*cos(ltt)^2), sum(w*cos(ltt)*sin(ltt)), -A*sum(w*tt*sin(ltt)^2)+B*sum(w*tt*cos(ltt)*sin(ltt)), sum(w*cos(ltt)*sin(ltt)), sum(w*sin(ltt)^2)),3,3) J1 <- solve(J1) J2<-matrix(c(A^2*sum(w^2*tt^2*sin(ltt)^2)+B^2*sum(w^2*tt^2*cos(ltt)^2) -2*A*B*sum(w^2*tt^2*cos(ltt)*sin(ltt)), -A*sum(w^2*tt*cos(ltt)*sin(ltt))+B*sum(w^2*tt*cos(ltt)^2), -A*sum(w^2*tt*sin(ltt)^2)+B*sum(w^2*tt*cos(ltt)*sin(ltt)), -A*sum(w^2*tt*cos(ltt)*sin(ltt))+B*sum(w^2*tt*cos(ltt)^2), sum(w^2*cos(ltt)^2), sum(w^2*cos(ltt)*sin(ltt)), -A*sum(w^2*tt*sin(ltt)^2)+B*sum(w^2*tt*cos(ltt)*sin(ltt)), sum(w^2*cos(ltt)*sin(ltt)), sum(w^2*sin(ltt)^2)),3,3) aux <- (J1%*%J2%*%J1) c(aux[1,1],2/N^3*C0/(A^2+B^2),1/(A^2+B^2)*(A^2*aux[2,2]+B^2*aux[3,3]+2*A*B*aux[2,3]),2/N*(U0/W0^2)) }) par(mfg=c(1,1,1,2)) lines(cycles,100*abs(vars[1,]-vars[2,])/vars[1,],lty=j) par(mfg=c(1,2,1,2)) lines(cycles,100*abs(vars[3,]-vars[4,])/vars[3,],lty=j) } par(mfg=c(1,1,1,2)) legend(5,90,sapply(as.character(Ns),function(x) paste("N=",x,sep="")),lty=1:length(Ns)) par(mfg=c(1,2,1,2)) legend(5,90,sapply(as.character(Ns),function(x) paste("N=",x,sep="")),lty=1:length(Ns)) dev.off()