Nothing
## $Id: ppw.R,v 1.47 2008/12/12 14:52:17 jracine Exp jracine $
## Original code in Matlab by A. Patton, R translation and
## modifications by C. Parmeter and J. Racine.
##
## We are grateful to Andrew Patton and Dimitris Politis for their
## assistance and feedback. Kindly report features, deficiencies, and
## improvements to racinej@mcmaster.ca.
##
## The citation is A. Patton, D.N. Politis, and H. White (2008,
## forthcoming), "CORRECTION TO `Automatic Block-Length Selection for
## the Dependent Bootstrap' by D.N. Politis and H. White". This is
## based on the article by Politis, D.N., and H. White (2004),
## "Automatic block-length selection for the dependent bootstrap."
## Econometric Reviews, vol. 23.
##
## INPUTS: data, an n x k matrix.
##
## OUTPUTS: b.star, a 2 x k vector of optimal bootstrap block lengths
## for the stationary bootstrap and circular bootstrap (BstarSB,
## BstarCB).
## The function lam() is used to construct a "flat-top" lag window for
## spectral estimation based on Politis, D.N. and J.P. Romano (1995),
## "Bias-Corrected Nonparametric Spectral Estimation", Journal of Time
## Series Analysis, vol. 16, No. 1.
lam <- function(s){
return((abs(s)>=0)*(abs(s)<0.5)+2*(1-abs(s))*(abs(s)>=0.5)*(abs(s)<=1))
}
## The function b.star() returns the optimal bootstrap block
## lengths. Note that an example for usage appears at the bottom of
## this file. If you use this function as input into a routine such as
## tsboot() in the boot library (Angelo Canty and Brian Ripley
## (2008). boot: Bootstrap R (S-Plus) Functions. R package version
## 1.2-34.) you ought to use the option round=TRUE.
b.star <- function(data,
Kn = NULL,
mmax= NULL,
Bmax = NULL,
c = NULL,
round = FALSE){
## Convert the data object to a data frame to handle both vectors
## and matrices.
data <- data.frame(data)
n <- nrow(data)
k <- ncol(data)
## Set Defaults. Note that in footnote c, page 59, for Kn Politis
## and White (2004) use max(5,log10(n)). Since this must be an
## integer we use ceiling(log10(n)).
if (is.null(Kn)) Kn <- max(5,ceiling(log10(n)))
if (is.null(mmax)) mmax <- ceiling(sqrt(n))+Kn
if (is.null(Bmax)) Bmax <- ceiling(min(3*sqrt(n),n/3))
if (is.null(c)) c <- qnorm(0.975)
## Create two vectors of length k in which we store results.
BstarSB <- numeric(length=k)
BstarCB <- numeric(length=k)
## Now we loop through each variable in data (i.e., column,
## data[,i]).
for(i in 1:k) {
## We first obtain the autocorrelations rho(1),...,rho(mmax) (we
## need to drop the first autocorrelation as it is rho(0), hence
## acf[-1]). This is the default in acf [type="correlation"]. Note
## that Patton uses sample correlations after dropping the first
## mmax observations, while we instead use the acf to obtain
## rho(k).
rho.k <- acf(data[,i],
lag.max = mmax,
type = "correlation",
plot = FALSE)$acf[-1]
## Next we compute mhat. The use of c*sqrt(log10(n)/n) for
## critical values is given in footnote c of Politis and White
## (2004, page 59), and the approach for determining mhat is
## described in footnote c.
rho.k.crit <- c*sqrt(log10(n)/n)
## Compute the number of insignificant runs following each rho(k),
## k=1,...,mmax.
num.insignificant <- sapply(1:(mmax-Kn+1),
function(j){
sum((abs(rho.k) < rho.k.crit)[j:(j+Kn-1)])
})
## If there are any values of rho(k) for which the Kn proceeding
## values of rho(k+j), j=1,...,Kn are all insignificant, take the
## smallest rho(k) such that this holds (see footnote c for
## further details).
if(any(num.insignificant==Kn)) {
mhat <- which(num.insignificant==Kn)[1]
} else {
## If no runs of length Kn are insignificant, take the smallest
## value of rho(k) that is significant.
if(any(abs(rho.k) > rho.k.crit)) {
lag.sig <- which(abs(rho.k) > rho.k.crit)
k.sig <- length(lag.sig)
if(k.sig == 1) {
## When only one lag is significant, mhat is the sole
## significant rho(k).
mhat <- lag.sig
} else {
## If there are more than one significant lags but no runs
## of length Kn, take the largest value of rho(k) that is
## significant.
mhat <- max(lag.sig)
}
} else {
## When there are no significant lags, mhat must be the
## smallest positive integer (footnote c), hence mhat is set
## to one.
mhat <- 1
}
}
## Compute M (mhat is at least one).
M <- ifelse(2*mhat > mmax, mmax, 2*mhat)
## We compute BstarSB and BstarCB using the formulas in the above
## references. Now we require the autocovariance R(k) (hence
## type="covariance" in the acf call). Note that Patton uses
## sample covariances after dropping the first mmax observations,
## while we instead use the acf with type="covariance" to obtain
## R(k). Note also that we require R(0) hence we do not drop it as
## we did for rho(k) via acf(...)$acf[-1].
kk <- seq(-M,M)
R.k <- ccf(data[,i], data[,i],
lag.max = M,
type = "covariance",
plot = FALSE)$acf
Ghat <- sum(lam(kk/M)*abs(kk)*R.k)
DCBhat <- 4/3*sum(lam(kk/M)*R.k)^2
DSBhat <- 2*sum(lam(kk/M)*R.k)^2
BstarSB[i] <- ((2*Ghat^2)/DSBhat)^(1/3)*n^(1/3)
BstarCB[i] <- ((2*(Ghat^2)/DCBhat)^(1/3))*(n^(1/3))
}
## The user can choose whether they want rounded values returned or
## not. BstarCB is rounded up, BstarSB simply rounded but both must
## be positive integers.
if(round == FALSE) {
BstarSB <- ifelse(BstarSB > Bmax, Bmax, BstarSB)
BstarCB <- ifelse(BstarCB > Bmax, Bmax, BstarCB)
} else {
BstarSB <- ifelse(BstarSB > Bmax, Bmax, ifelse(BstarSB < 1, 1, round(BstarSB)))
BstarCB <- ifelse(BstarCB > Bmax, Bmax, ifelse(BstarCB < 1, 1, max(1,round(BstarCB))))
}
return(cbind(BstarSB,BstarCB))
}
## Here is a simple example with an n x 2 matrix containing n=10^5
## observations, where column 2 of x is more persistent than column
## 1. This requires that you first install the forecast library (i.e.,
## install.packages("forecast")).
##
## library(forecast)
## set.seed(123)
## x <- cbind(arima.sim(n = 100000, list(ar = c(.5,.0), ma = c(0,0)),sd = 1),
## arima.sim(n = 100000, list(ar = c(.5,.4), ma = c(0,0)),sd = 1))
## b.star(x)
## b.star(x,round=TRUE)
##> b.star(x)
## BstarSB BstarCB
##[1,] 50.39272 57.68526
##[2,] 251.62894 288.04323
##> b.star(x,round=TRUE)
## BstarSB BstarCB
##[1,] 50 58
##[2,] 252 289
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.