`f.hf1` <- function(cex=1){ cat("Haar Fisz Sim\n") op <- par(mfrow=c(2,2), cex=cex) x <- rpois(1000, lambda=2) y <- rpois(1000, lambda=2) zeta <- fisz(x,y) - fisz(2,4) A <- sqrt(x+ 3/8) bw="bcv" plot(density(A, bw=bw, from=0), main="a.", xlab="A") qqnorm(A, main="b.") plot(density(zeta, bw=bw, from=-3, to=3), main="c.", xlab="zeta") qqnorm(zeta, main="d.") cat("Var of A is ", var(2*A), "\n") cat("Var of zeta is ", var(zeta), "\n") par(op) } `f.hf2` <- function(cex=1){ op <- par(mfrow=c(2,2), cex=cex) #lambda <- DJ.EX(n=512)$blocks #lambda <- (lambda - min(lambda) + 1)/8 lambda <- c(rep(1,256), rep(10,256)) n <- length(lambda) X <- rpois(n, lambda=lambda) plot(1:n, lambda, type="l", xlab="i", ylab="lambdai") plot(1:n, X, type="l", xlab="i", ylab="Xi") XHF <- hft(X) XAn <- 2*sqrt(X+3/8) plot(1:n, XHF, type="l", xlab="i", ylab="FXi") plot(1:n, XAn, type="l", xlab="i", ylab="Anscombei") vhf1 <- var(XHF[1:256]) vhf2 <- var(XHF[257:512]) va1 <- var(XAn[1:256]) va2 <- var(XAn[257:512]) cat("Variance of HF on first section is ", vhf1, "\n") cat("Variance of HF on second section is ", vhf2, "\n") cat("Variance of Anscombe on first section is ", va1, "\n") cat("Variance of Anscome on second section is ", va2, "\n") l <- list(lambda=lambda, X=X, XHF=XHF, XAn=XAn) par(op) l } `f.hf3` <- function(obj.hf3, cex=1){ op <- par(mfrow=c(2,2), cex=cex) if (missing(obj.hf3)) { lambda <- DJ.EX(n=512)$blocks lambda <- (lambda - min(lambda) + 1)/8 #lambda <- c(rep(1,256), rep(10,256)) n <- length(lambda) X <- rpois(n, lambda=lambda) } else { X <- obj.hf3$X lambda <- obj.hf3$lambda n <- length(lambda) } plot(1:n, lambda, type="l") plot(1:n, X, type="l") XHF <- hft(X) XAn <- 2*sqrt(X+3/8) plot(1:n, XHF, type="l") plot(1:n, XAn, type="l") l <- list(lambda=lambda, X=X, XHF=XHF, XAn=XAn) par(op) l } `f.hf4` <- function(obj.hf3, ll=3, cex=1){ op <- par(mfrow=c(2,2), cex=cex) XHF <- obj.hf3$XHF XAn <- obj.hf3$XAn n <- length(XAn) XHFwd <- wd(XHF, filter.number=1, family="DaubExPhase") XAnwd <- wd(XAn, filter.number=1, family="DaubExPhase") XHFwdT <- ebayesthresh.wavelet(XHFwd) XAnwdT <- ebayesthresh.wavelet(XAnwd) XHFdn <- wr(XHFwdT) XAndn <- wr(XAnwdT) plot(1:n, XHFdn, type="l", xlab="i.") plot(1:n, XAndn, type="l", xlab="i.") XHFest <- hft.inv(XHFdn) XAnest <- ((XAndn)/2)^2 - 3/8 lambda <- obj.hf3$lambda plot(1:n, XHFest, type="l", xlab="i.") lines(1:n, lambda, lty=2) errHF <- sqrt(sum((XHFest-lambda)^2)) plot(1:n, XAnest, type="l", xlab="i.") lines(1:n, lambda, lty=2) errAn <- sqrt(sum((XAnest-lambda)^2)) cat("HF error is ", errHF, "\n") cat("An error is ", errAn, "\n") par(op) } `f.hf5` <- function(cex=1){ op <- par(cex=cex) ts.plot(AirPassengers, xlab="Year") par(op) } `f.hf6` <- function(cex=1){ op <- par(cex=cex) ts.plot(AirPassengers, xlab="Year") par(op) } `f.hf7` <- function(cex=1){ op <- par(cex=cex) AP128 <- AirPassengers[1:128] AP128 <- ts(AP128, start=c(1949,1), frequency=12) APhft <- ddhft.np.2(AP128) plot(APhft$mu, sqrt(APhft$sigma2), xlab="MU", ylab="SIGMA") lines(APhft$mu, APhft$sigma) lines(APhft$mu, sqrt(APhft$mu), lty=2) par(op) } `f.hf8` <- function(cex=1){ op <- par(cex=cex) AP128 <- AirPassengers[1:128] AP128 <- ts(AP128, start=c(1949,1), frequency=12) APhft <- ddhft.np.2(AP128) APhftTS <- ts(APhft$hft, start=c(1949,1), frequency=12) ts.plot(APhftTS, xlab="Year") AP128sqrt <- sqrt(AP128) # # Do linear transform to put AP128sqrt on same scale as APhftTS # AP128sqrt <- (AP128sqrt - min(AP128sqrt))/(max(AP128sqrt)-min(AP128sqrt)) AP128sqrt <- AP128sqrt*(max(APhftTS)-min(APhftTS)) + min(APhftTS) lines(AP128sqrt, lty=2) cat("var First 32 HF: ", var(APhftTS[1:32]), "\n") cat("var Last 32 HF: ", var(APhftTS[97:128]), "\n") cat("var First 32 sqrt: ", var(AP128sqrt[1:32]), "\n") cat("var Last 32 sqrt: ", var(AP128sqrt[97:128]), "\n") AP128log <- log(AP128) # # Do linear transform to put AP128log on same scale as APhftTS # AP128log <- (AP128log - min(AP128log))/(max(AP128log)-min(AP128log)) AP128log <- AP128log*(max(APhftTS)-min(APhftTS)) + min(APhftTS) lines(AP128log, lty=4) cat("var First 32 log: ", var(AP128log[1:32]), "\n") cat("var Last 32 log: ", var(AP128log[97:128]), "\n") par(op) }