stationarity <- function(T=1000, MC=100000, deg_f=Inf, lambda=0, rho=0, c=FALSE, s=0, k=FALSE, lags=0, bootstrap=FALSE, p=.1) { ############################################################################################ # # # Stationarity Function # # # # Breno Neri (http://homepages.nyu.edu/~bpn207) # # Department of Economics, New York University # # # # August, 2008 - Version 9 # # # ############################################################################################ # # # Parameter: Default: Description: # # # # T 1000 Number of Observations # # MC 100000 Number of Replications # # deg_f Inf t Distribution's Degrees of Freedom (Gaussian=Inf) # # lambda 0 Random Walk Parameter in The DGP (lambda>=0) # # rho 0 Auto-Regressive Coefficient # # c FALSE Local to Finite Variance # # s 0 Scaling Factor # # k FALSE Varying Kurtosis Component # # lags 0 Lags in Kernel: 0 for 0 # # 1 for round(4*(T/100)^(1/4)) # # 2 for data dependent # # bootstrap FALSE If Bootstrap Is to Be Used # # 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 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) n <- sqrt(2) + 8 * 1:T / T cv <- c(.463, .463, 2.07, 1.77) results <- array(0, c(MC, 4)) results2 <- array(0, c(MC, 4)) fr <- numeric() # 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) # Varying Scale if(s) scale <- sqrt(1+s*1:T) else scale <- 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))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)) } # Monte Carlo for(mc in 1:MC) { # Local to Finie Variance if(c) local1 <- rt(T, df=1) / sqrt(T) else local1 <- 0 # Varying Kurtosis if(k) kurt1 <- n/sqrt(2) * (as.numeric((e1<-runif(T))>1-1/(n*n)) - as.numeric(e1<1/(n*n))) else kurt1 <- 0 # Generating Unir Root Component if(lambda) { if(c) local2 <- rt(T, df=1) / sqrt(T) else local2 <- 0 if(k) kurt2 <- n/sqrt(2) * (as.numeric((e2<-runif(T))>1-1/(n*n)) - as.numeric(e2<1/n*n)) else kurt2 <- 0 r <- lambda * cumsum(scale*(rt(T, df=deg_f) + local2 + kurt2)) } else r <- 0 # Defining DGP e <- scale * (rt(T, df=deg_f) + local1 + kurt1) if(rho) { e.ar <- e[1] for(t in 2:T) e.ar[t] <- rho*e.ar[t-1] + e[t] x <- r + e.ar } else x <- r + e # Generating Results gamma <- n.lags(abs(ar(x, order.max=1)[[2]])) results[mc,] <- stats(x) if(bootstrap) { pos <- numeric() u <- .Internal(sample(T, T, TRUE, NULL)) l <- rgeom(T, p) i <- 1 while(length(pos)quantile(results2[,i], .95)) } else fr <- apply(t(results)>cv, 1, mean) # Printing Output cat("\n\n\n--- Parameters ---", fill=TRUE) cat("T:", T, fill=TRUE) cat("MC:", MC, fill=TRUE) cat("deg_f:", deg_f, fill=TRUE) cat("lambda:", lambda, fill=TRUE) cat("rho:", rho, fill=TRUE) cat("c:", c, fill=TRUE) cat("s:", s, fill=TRUE) cat("k:", k, fill=TRUE) cat("lags:", lags, fill=TRUE) cat("bootstrap:", bootstrap, fill=TRUE) cat("p:", p, fill=TRUE) cat("\n--- Frequency of Rejection ---", fill=TRUE) cat("KPSS:", fr[1], fill=TRUE) cat("IKPSS:", fr[2], fill=TRUE) cat("XL:", fr[3], fill=TRUE) cat("SS:", fr[4], fill=TRUE) }