## $Id: simul_ar1.R,v 1.8 2008/12/11 11:10:50 jracine Exp jracine $

## This code simulates Table 1 along with the summaries reported in
## Table 2 (i.e., boxplots are for simulated optimal block lengths
## used to construct Table 2 that have not been divided by the optimal
## block length in Table 1, but rmse corresponds to Table 2) in
## A. Patton, D.N. Politis, and H. White, `CORRECTION TO ``Automatic
## Block-Length Selection for the Dependent Bootstrap'' by D. Politis
## and H. White'. Note that as our code is slightly different from
## that used in Patton et al which was written in Matlab. Thus our
## results differ ever so slightly due to random sampling variation
## and the use of the acf rather than correlations and covariances
## having dropped the first mmax observations as is done the Patton's
## Matlab code.

rm(list=ls())

## The file ppw.R must be in the same library as this file otherwise
## you need to specify the location in source("ppw.R").

source("ppw.R")
library(forecast)

set.seed(12345)

do.round <- FALSE

## Set the number of Monte Carlo replications.

M <- 1000

## We replicate simulations for the stationary and circular
## bootstrap. These functions return the optimal block lengths for the
## stationary and circular bootstraps.

opt.sb <- function(n) {(4*(rho^2)/((1-(rho^2))^2) )^(1/3) * n^(1/3)}
opt.cb <- function(n) {(6*(rho^2)/((1-(rho^2))^2) )^(1/3) * n^(1/3)}

## Set values for rho and n. You can change these as you wish.

rho.vals <- c(0.7,0.1,-0.4)
n.vals <- c(200,800)

Bstar <- numeric(length=M)

## First, run the simulation for the stationary bootstrap.

## Plot the six boxplots and summary information in one plot and write
## this to boxplot_sb.pdf.

pdf(file="boxplot_sb.pdf")
par(mfrow=c(3,2))

for(rho in rho.vals) {
  for(n in n.vals) {
    for(i in 1:M) {
      Bstar[i] <- b.star(arima.sim(n = n, list(ar = c(rho)),sd = 1),round=do.round)[1,1]
    }

    ## Making sure SB block length is weakly greater than 1: A. Patton

    Bstar <- ifelse(Bstar<1,1,Bstar)

    ## Write results to a file.

    write(Bstar,file=paste("Bstar_SB_",rho,"_",n,".out",sep=""),ncol=1)
    
    rmse <- sqrt(mean((Bstar/opt.sb(n)-1)^2))

    boxplot(Bstar,horizontal=TRUE,notch=TRUE,
            main=paste("rho=",rho,
              ", n=",n,
              ", median=",signif(median(Bstar),digits=3),
              ", mean=",signif(mean(Bstar),digits=3),sep=""),
            sub=paste("median.rel=",signif(median(Bstar/opt.sb(n)),digits=3),
              ", mean.rel=",signif(mean(Bstar/opt.sb(n)),digits=3),
              ", opt=",signif(opt.sb(n),digits=3),
              ", rmse=",signif(rmse,digits=3),sep=""))

  }
}

dev.off()

## Second, run the simulation for the circular bootstrap.

## Plot the six boxplots and summary information in one plot and write
## this to boxplot_cb.pdf.

pdf(file="boxplot_cb.pdf")
par(mfrow=c(3,2))

for(rho in rho.vals) {
  for(n in n.vals) {
    for(i in 1:M) {
      Bstar[i] <- b.star(arima.sim(n = n, list(ar = c(rho)),sd = 1),round=do.round)[1,2]
    }

    ##  Rounding CB block length up to nearest integer: A. Patton
    
    Bstar <- ceiling(Bstar)

    ## Write results to a file.

    write(Bstar,file=paste("Bstar_CB_",rho,"_",n,".out",sep=""),ncol=1)

    rmse <- sqrt(mean((Bstar/opt.cb(n)-1)^2))

    boxplot(Bstar,horizontal=TRUE,notch=TRUE,
            main=paste("rho=",rho,
              ", n=",n,
              ", median=",signif(median(Bstar),digits=3),
              ", mean=",signif(mean(Bstar),digits=3),sep=""),
            sub=paste("median.rel=",signif(median(Bstar/opt.cb(n)),digits=3),
              ", mean.rel=",signif(mean(Bstar/opt.cb(n)),digits=3),
              ", opt=",signif(opt.cb(n),digits=3),
              ", rmse=",signif(rmse,digits=3),sep=""))

  }
}

dev.off()