stationarity    <- function(x, R=1000, bootstrap=TRUE, lags=1, p=.1)
{
############################################################################################
#                                                                                          #
# Stationarity Test                                                                        #
#                                                                                          #
# Breno Neri <breno.neri@nyu.edu> (http://homepages.nyu.edu/~bpn207)                       #
# Department of Economics, New York University                                             #
#                                                                                          #
#                                                                 August, 2008 - Version 3 #
#                                                                                          #
############################################################################################
#                                                                                          #
# Parameter:    Default:                Description:                                       #
#                                                                                          #
# R             1000                    Number of Replications                             #
# bootstrap     TRUE                    Bootstrap (recommended)                            #
# lags          1                       Lags in Kernel: 0 for 0                            #
#                                                       1 for round(4*(T/100)^(1/4))       #
#                                                       2 for data dependent               #
# p             .1                      Probability of The Geometric Distribution          #
#                                                                                          #
############################################################################################
#                                                                                          #
# You may use this code only if you accept the conditions:                                 #
# (1) I am not liable for any problem caused by this code (i.e. you use it at your own     #
# risk);                                                                                   #
# (2) You must give me credit in your papers where this code has been used.                #
#                                                                                          #
############################################################################################

    # Inicializing
    T                       <- length(x)
    tau                     <- seq(.1, .9, .01)
    qs                      <- length(tau)
    Mtau                    <- t(array(tau, c(qs, T)))
    A                       <- round(4*(T/100)^(1/4))
    B                       <- (6*T)^(1/3)
    C                       <- 8*(T/100)^(1/3)
    results2                <- array(0, c(R, 4))
    cv                      <- array(0, c(3, 4))
    cv[1,]                  <- c(0.347, 0.347, 1.89, 1.65)
    cv[2,]                  <- c(0.463, 0.463, 2.07, 1.77)
    cv[3,]                  <- c(0.739, 0.739, 2.40, 2.01)

    # Defining Kernel
    if(lags==0)      n.lags     <- function(r) 0
    else if(lags==1) n.lags     <- function(r) A
    else if(lags==2) n.lags     <- function(r) round(min(B*(r/(1-r*r))^(2/3), C))
    else stop("lags must be 0, 1 or 2!")
    K                       <- function(j) 1 - j/(gamma+1)

    # Defining Tests Statistics Functions
    stats                   <- function(x)
    {
        # Creating Artificial Data
        y.KPSS                  <- x - mean(x)
        y.IKPSS                 <- sign(x-median(x))
        y.XL                    <- y.KPSS*y.KPSS - (s11 <- y.KPSS%*%y.KPSS) / T
        psi                     <- (array(x, c(T, qs))<t(array(quantile(x, tau), c(qs, T)))) - Mtau

        # Computing Long Run Variances
        sigma2.KPSS             <- y.KPSS %*% y.KPSS
        sigma2.IKPSS            <- T
        s12                     <- y.KPSS %*% y.XL
        s22                     <- y.XL %*% y.XL
        pi2                     <- apply(psi*psi, 2, sum)
        if(gamma>0) for(j in 1:gamma)
        {
            sigma2.KPSS             <- sigma2.KPSS + 2*K(j)*y.KPSS[1:(T-j)]%*%y.KPSS[(1+j):T]
            sigma2.IKPSS            <- sigma2.IKPSS + 2*K(j)*y.IKPSS[1:(T-j)]%*%y.IKPSS[(1+j):T]
            s11                     <- s11 + 2*K(j)*y.KPSS[1:(T-j)]%*%y.KPSS[(1+j):T]
            s12                     <- s12 + K(j) * (y.KPSS[1:(T-j)]%*%y.XL[(1+j):T] + y.XL[1:(T-j)]%*%y.KPSS[(1+j):T])
            s22                     <- s22 + 2*K(j)*y.XL[1:(T-j)]%*%y.XL[(1+j):T]
            pi2                     <- pi2 + 2*K(j)*apply(psi[1:(T-j),]*psi[(1+j):T,], 2, sum)      
        }
        sigma2.XL               <- matrix(c(s11, s12, s12, s22), nrow=2)

        # Computing Test Statistic
        stat.KPSS               <- sum((c.y.KPSS<-cumsum(y.KPSS)) * c.y.KPSS) / (sigma2.KPSS*T)
        stat.IKPSS              <- sum((c.y.IKPSS<-cumsum(y.IKPSS)) * c.y.IKPSS) / (sigma2.IKPSS*T)
        stat.XL                 <- max(apply(abs(solve(t(chol(sigma2.XL))) %*% rbind(cumsum(y.KPSS), cumsum(y.XL))), 2, sum))
        stat                    <- max(abs(apply(psi, 2, cumsum) - outer(1:T, apply(psi, 2, mean))) / t(array(sqrt(pi2), c(qs, T))))

        return(c(stat.KPSS, stat.IKPSS, stat.XL, stat))
    }

    gamma                   <- n.lags(abs(ar(x, order.max=1)[[2]]))
    results                 <- stats(x)

    # Bootstrap
    if(bootstrap)
    {
        for(r in 1:R)
        {
            pos                     <- numeric()
            u                       <- .Internal(sample(T, T, TRUE, NULL))
            l                       <- rgeom(T, p)
            i                       <- 1
            while(length(pos)<T)
            {
                 pos                    <- c(pos, u[i]:min(u[i]+l[i], T))
                 i                      <- i + 1
            }
            x.resampled             <- x[pos[1:T]]
            gamma                   <- n.lags(abs(ar(x.resampled, order.max=1)[[2]]))
            results2[r,]            <- stats(x.resampled)
        }
        pv                      <- apply(t(results2)>results, 1, mean)

        # Printing Output
        cat("--- P-Values ---", fill=TRUE)
        cat("KPSS:", pv[1], fill=TRUE)
        cat("IKPSS:", pv[2], fill=TRUE)
        cat("XL:", pv[3], fill=TRUE)
        cat("SS:", pv[4], fill=TRUE)
    } else
    {
        cat("--- Test Statistics ---", fill=TRUE)
        cat("KPSS:", results[1], fill=TRUE)
        cat("IKPSS:", results[2], fill=TRUE)
        cat("XL:", results[3], fill=TRUE)
        cat("SS:", results[4], fill=TRUE)
        cat("--- Critical Values at 10% Significance Level ---", fill=TRUE)
        cat("KPSS:", cv[1, 1], fill=TRUE)
        cat("IKPSS:", cv[1, 2], fill=TRUE)
        cat("XL:", cv[1, 3], fill=TRUE)
        cat("SS:", cv[1, 4], fill=TRUE)
        cat("--- Critical Values at 5% Significance Level ---", fill=TRUE)
        cat("KPSS:", cv[2, 1], fill=TRUE)
        cat("IKPSS:", cv[2, 2], fill=TRUE)
        cat("XL:", cv[2, 3], fill=TRUE)
        cat("SS:", cv[2, 4], fill=TRUE)
        cat("--- Critical Values at 1% Significance Level ---", fill=TRUE)
        cat("KPSS:", cv[3, 1], fill=TRUE)
        cat("IKPSS:", cv[3, 2], fill=TRUE)
        cat("XL:", cv[3, 3], fill=TRUE)
        cat("SS:", cv[3, 4], fill=TRUE)
    }
}