## Examples of Applied Time Series Analysis, By He Shuyuan ## Chpater 4, Estimation of the mean and the ACV ## estimation of mean of an AR(2) demo.tsmean <- function(){ roots <- complex(mod=c(1/1.1, 1/4, 1/4, 1/1.5), arg=c(pi*2/3, pi*2/3, pi*5/6, pi/6)) coefs <- cbind(2*Re(roots), -Mod(roots)^2) print(apply(coefs, 1, function(a) 1/(1 - a[1] - a[2]))) L <- 1000 M <- 500 sigma <- 1 ## first simu: effect of sample size n n0 <- 200 n2 <- n0 + L eps <- rnorm(n2, 0, sigma) xx <- filter(eps, coefs[1,,drop=T], method="recursive", side=1) x1 <- xx[(n0+1):n2] xx <- filter(eps, coefs[2,,drop=T], method="recursive", side=1) x2 <- xx[(n0+1):n2] barx1 <- cumsum(x1)/seq(L) barx2 <- cumsum(x2)/seq(L) bareps <- cumsum(eps[(n0+1):n2])/seq(L) yr <- range(c(barx1+0.1, bareps, barx2-0.1)) plot(bareps, type="l", col="red", lwd=1.5, ylim=c(-0.3, 0.3), main="Effect of sample size N", xlab="N", ylab="Avgs") lines(barx1+0.1) lines(barx2-0.1) abline(h=c(0, 0.1, -0.1)) locator(1) ## sampling distributions ns <- c(10, 20, 50, 100, 200, 1000) nn <- length(ns) bx1s <- matrix(0, nrow=M, ncol=nn) bx2s <- matrix(0, nrow=M, ncol=nn) for(mm in seq(M)){ eps <- rnorm(n2, 0, sigma) xx <- filter(eps, coefs[3,,drop=T], method="recursive", side=1) x1 <- xx[(n0+1):n2] xx <- filter(eps, coefs[4,,drop=T], method="recursive", side=1) x2 <- xx[(n0+1):n2] barx1 <- cumsum(x1)/seq(L) barx2 <- cumsum(x2)/seq(L) bx1s[mm,] <- barx1[ns] bx2s[mm,] <- barx2[ns] } avgs1 <- apply(bx1s, 2, mean) stds1 <- apply(bx1s, 2, sd) avgs2 <- apply(bx2s, 2, mean) stds2 <- apply(bx2s, 2, sd) res <- rbind(avgs1, stds1, avgs2, stds2) rownames(res) <- c("Avg1", "Std1", "Avg2", "Std2") colnames(res) <- ns cat("\nSampling distribution by repeated sampling. ===\n") print(res) for(ii in seq(ns)){ hist(bx1s[,ii], main=paste("Sampling distribution: N=", ns[ii]), xlim=c(-0.7,0.7)) locator(1) } } f <- function(){ eps <- rnorm(1000,0) bare <- cumsum(eps)/seq(1000) plot(eps, type="l", col="green") lines(bare, col="red") abline(h=0) } ## Estimation of AR(2) \gamma_k demo.tsacv <- function(){ roots <- complex(mod=c(1/1.8, 1/1.1), arg=c(2.13, 2.13)) coefs <- cbind(2*Re(roots), -Mod(roots)^2) sigma <- 1 ngam <- 31 ## first simu: effect of sample size n L <- 1600 x <- ar.gen(L, a=coefs[1,]) true.gam <- arma.gamma.by.Wold(ngam, a=coefs[1,], sigma=sigma) est.gam <- acf(x[1:100], lag.max=ngam-1, type="covariance", plot=F)$acf plot(0:(ngam-1), true.gam, type="h", main="ACV of AR(2): n=100", ylim=range(c(true.gam,est.gam)), xlab="k", ylab="") abline(h=0) points(0:(ngam-1), est.gam) locator(1) est.gam <- acf(x[1:400], lag.max=ngam-1, type="covariance", plot=F)$acf plot(0:(ngam-1), true.gam, type="h", main="ACV of AR(2): n=400", ylim=range(c(true.gam,est.gam)), xlab="k", ylab="") abline(h=0) points(0:(ngam-1), est.gam) locator(1) est.gam <- acf(x[1:1600], lag.max=ngam-1, type="covariance", plot=F)$acf plot(0:(ngam-1), true.gam, type="h", main="ACV of AR(2): n=1600", ylim=range(c(true.gam,est.gam)), xlab="k", ylab="") abline(h=0) points(0:(ngam-1), est.gam) ## Next use repeated simulation to assess variation a <- coefs[1,] M <- 400 ns <- c(1600, 400, 100, 25) L <- max(ns) ngam <- 6 stderrs <- matrix(0, nrow=length(ns), ncol=ngam) colnames(stderrs) <- 0:(ngam-1) rownames(stderrs) <- ns gams <- as.list(ns) for(ii in seq(along=gams)){ gams[[ii]] <- matrix(0, nrow=M, ncol=ngam) } true.gam <- arma.gamma.by.Wold(ngam, a=a, sigma=sigma) for(mm in seq(M)){ x <- ar.gen(L, a=a) for(ig in seq(along=ns)){ n <- ns[ig] est.gam <- acf(x[1:n], lag.max=ngam-1, type="covariance", plot=F)$acf gams[[ig]][mm,] <- est.gam } } rats <- numeric(length(ns)) rats[1] <- NA for(ig in seq(along=ns)){ stderrs[ig,] <- apply(gams[[ig]], 2, sd) if(ig>1) rats[ig] <- mean(stderrs[ig]/stderrs[ig-1]) } cat("\nSimulation of a more stable AR, repetitions=", M, " ===\n") cat("True gamma's\n") print(true.gam[0:(ngam-1)]) cat("Standard deviation of \\hat\\gamma_k ===\n") print(cbind(stderrs, "ratio"=rats)) ## Use AR(2) with roots near unit circle a <- coefs[2,] stderrs <- matrix(0, nrow=length(ns), ncol=ngam) colnames(stderrs) <- 0:(ngam-1) rownames(stderrs) <- ns gams <- as.list(ns) for(ii in seq(along=gams)){ gams[[ii]] <- matrix(0, nrow=M, ncol=ngam) } true.gam <- arma.gamma.by.Wold(ngam, a=a, sigma=sigma) for(mm in seq(M)){ x <- ar.gen(L, a=a) for(ig in seq(along=ns)){ n <- ns[ig] est.gam <- acf(x[1:n], lag.max=ngam-1, type="covariance", plot=F)$acf gams[[ig]][mm,] <- est.gam } } rats <- numeric(length(ns)) rats[1] <- NA for(ig in seq(along=ns)){ stderrs[ig,] <- apply(gams[[ig]], 2, sd) if(ig>1) rats[ig] <- mean(stderrs[ig]/stderrs[ig-1]) } cat("\nSimulation of a less stable AR, repetitions=",M, " ===\n") cat("True gamma's\n") print(true.gam[0:(ngam-1)]) cat("Standard deviation of \\hat\\gamma_k ===\n") print(cbind(stderrs, "ratio"=rats)) invisible() } ## Chi-square White Noise Test ## Input: ## x -- time series ## m -- number of terms in chi-square statistic ## Output: p-value of test wn.chisq <- function(x, m=5){ n <- length(x) rhos <- c(acf(x, lag.max=m+1, plot=F)$acf)[-1] stat <- sum(rhos^2)*n pr <- 1 - pchisq(stat, m) pr } ## chi-square white noise test for AR(2) demo.wn.chisq.ar2 <- function(){ cat("==== Chi-square White Noise Test for AR(2) ====\n") mods <- c("1000", "10", "6", "4", "3", "2", "1.5") revroots <- complex(mod=1/c(1000, 10, 6, 4, 3, 2, 1.5), arg=1.13) coefs <- cbind(2*Re(revroots), -Mod(revroots)^2) nmodels <- length(revroots) sigma <- 1 N <- 400 ## sample size nrep <- 500 ## number repetitions alpha <- 0.05 ## test level m <- 5 ## number of terms used in chi-square statistic power <- numeric(nmodels) names(power) <- mods for(jj in seq(nmodels)){ pvs <- numeric(nrep) ## store p-values for(ii in seq(nrep)){ x <- ar.gen(N, a=coefs[jj,], sigma=sigma) pvs[ii] <- wn.chisq(x, m) } ## correct rate power[jj] <- mean(pvs