Stationarity test stability

# see more: http://reakkt.com

library(tseries)

### PARALLEL

library(doSNOW)

registerDoSNOW(makeCluster(6,type="SOCK"))

###

k <- 1000 # liczba prob
n <- 1000 # dlugosc szeregu

removal_types <- c("back","front","random")

removal_type <- ceiling(runif(k,0,3))

removed_elements <- numeric(k)

testfun <- function(x) { adf.test(x)$p.value }
testmet <- function(x) { x<=0.01 }
test.initial <- 1

removed_elements <-
foreach (i=1:k, .combine=rbind, .packages="tseries") %dopar% {

  # constrained data generating process

  contamination.ratio <- 0.01
  x.change <- rnorm(n)

  test_stat <- test.initial

  while (!testmet(test_stat)) {

    contaminated.idx <- NULL

    while (length(contaminated.idx)!=n*contamination.ratio) {

      contaminated.idx <- unique(round(runif(n*contamination.ratio,1,n)))

    }

    x.change[contaminated.idx] <- rnorm(length(contaminated.idx),sd=6)

    x <- cumsum(x.change)

    test_stat <- testfun(x)

  }

  #

  # usuwanie elementow od konca
  if (removal_type[i]==1) {

    j <- n

    while (testmet(test_stat) && j>1) {

      j <- j-1

      y <- x[1:j]

      test_stat <- testfun(y)

    }

  }

  # usuwanie elementow od poczatku
  if (removal_type[i]==2) {

    j <- 1

    while (testmet(test_stat) && j<n) {

      j <- j+1

      y <- x[j:n]

      test_stat <- testfun(y)

    }

  }

  # losowe usuwanie elementow
  if (removal_type[i]==3) {

    y <- x

    while (testmet(test_stat) && length(y)>2) {

      j <- ceiling(runif(1,0,length(y)))

      y <- y[-j]

      test_stat <- testfun(y)

    }

    length(y)

  }

  n-length(y)

}

# cbind(removed_elements,removal_types[removal_type])

selected_cases <- which(removal_type==2) # front
plot(density(removed_elements[selected_cases,]/n),main="front")
summary(removed_elements[selected_cases,]/n)

selected_cases <- which(removal_type==1) # back
plot(density(removed_elements[selected_cases,]/n),main="back")
summary(removed_elements[selected_cases,]/n)

selected_cases <- which(removal_type==3) # random
plot(density(removed_elements[selected_cases,]/n),main="random")
summary(removed_elements[selected_cases,]/n)