rm(list = ls()) #install.packages(c("stats","urca", "forecast")) library("stats") library("urca") library("forecast") library("xtable") getDFSizesForDGP <- function( p, type = "none", T_obs, reps, DGP = replicate( reps , arima.sim(n = T_obs, list(ar = c(0.6, 0.3)), sd = sqrt( 5 ))), alpha_s = 0.05){ # Calculate test statistics for each time series / column in the DGPs test_statistics <- apply(DGP, 2, function(x){ X <- matrix(NA, nrow=length(p), ncol=1) if(length(p) == 1){dimnames(X)[[1]] <- list(p) } else {dimnames(X)[[1]] <- p} for(i in p){X[as.character(i), ] <- ur.df(x, type = type, lags = i, selectlags = "Fixed")@teststat[1]} return(X)}) dimnames(test_statistics)[[1]] <- p # Find critical value depending on the test type (independent of p) if(alpha_s == 0.01){ index <- 1 } else if(alpha_s == 0.05){ index <- 2 } else if(alpha_s == 0.1){ index <- 3 } else{stop("Wrong alpha was chosen")} crit_value <- summary(ur.df(DGP[ ,1], type = type, lags = 0, selectlags = "Fixed"))@cval[1,index ] # How often was H_0 rejected? rejections <- rowMeans(test_statistics < crit_value) return(rejections) } # a simple test set.seed(1234) getDFSizesForDGP( p = 0:10, type = "none", T_obs = 200, reps = 100, DGP = replicate( 100 , arima.sim(n = 200, list(ar = c(0.6, 0.3)), sd = sqrt( 5 ))), alpha_s = 0.05) # now a extended simulation p <- 0:10 # a sequence of lags included in augmentation T_obs <- 50 # number of observations per sample reps <- 10000 # number of observations per sample error_variance <- 5 # variance of innovations alpha_s <- 0.05 # significance level types <- c("none", "drift", "trend") names(types) <- c("Type - 0", "Type - Mu", "Type - Tau") rejection_matrix <- matrix(NA, ncol = 3, nrow = length(p), dimnames = list(as.character(p), names(types))) set.seed(1234) for(type in types){ # type = "drift" rejection_matrix[ ,names(types[which(type == types)])] <- getDFSizesForDGP( p = p, type = type, T_obs = T_obs, reps = reps, DGP = replicate( reps , arima.sim(n = T_obs, list(ar = c(0.6, 0.3)), sd = sqrt( error_variance ))), alpha_s = alpha_s) } mat1_50 <- rejection_matrix # now a extended simulation p <- 0:10 # a sequence of lags included in augmentation T_obs <- 200 # number of observations per sample reps <- 10000 # number of observations per sample error_variance <- 5 # variance of innovations alpha_s <- 0.05 # significance level types <- c("none", "drift", "trend") names(types) <- c("Type - 0", "Type - Mu", "Type - Tau") rejection_matrix <- matrix(NA, ncol = 3, nrow = length(p), dimnames = list(as.character(p), names(types))) set.seed(1234) for(type in types){ # type = "drift" rejection_matrix[ ,names(types[which(type == types)])] <- getDFSizesForDGP( p = p, type = type, T_obs = T_obs, reps = reps, DGP = replicate( reps , arima.sim(n = T_obs, list(ar = c(0.6, 0.3)), sd = sqrt( error_variance ))), alpha_s = alpha_s) } mat1_200 <- rejection_matrix # now as non-stationary AR(2) with parameters 0.6 and 0.4; need a function ar2.gauss <- function(T, alpha , y.start = c(0,0,0), sigma = 1){ y <- rep(0, T+2) # initialize the output vector; # the length is set to T+3, because of the initial values y[1:2] <- y.start # overwrite the first three values with the specified vector of starting values y <- ts(y, start = 0) # transform vector y to a ts object # compute y iteratively; start for the fourth row of y, which is for time period 1 for(t in 3:(T+2)){ y[t] <- alpha[1] * y[t-1] + alpha[2] * y[t-2] + rnorm(1, 0 ,sigma) } return(window(y, start = 2)) # discard the initial values } ar2.gauss(T = T_obs, alpha =c(0.6,0.4), y.start = c(0,0), sigma = 1) p <- 0:10 # a sequence of lags included in augmentation T_obs <- 200 # number of observations per sample reps <- 10000 # number of observations per sample error_variance <- 5 # variance of innovations alpha_s <- 0.05 # significance level types <- c("none", "drift", "trend") names(types) <- c("Type - 0", "Type - Mu", "Type - Tau") rejection_matrix <- matrix(NA, ncol = 3, nrow = length(p), dimnames = list(as.character(p), names(types))) set.seed(1234) for(type in types){ # type = "drift" rejection_matrix[ ,names(types[which(type == types)])] <- getDFSizesForDGP( p = p, type = type, T_obs = T_obs, reps = reps, DGP = replicate( reps , ar2.gauss(T = T_obs, alpha =c(0.6,0.4), y.start = c(0,0), sigma = sqrt(error_variance))), alpha_s = alpha_s) } mat2_200 <- rejection_matrix # now a extended simulation p <- 0:10 # a sequence of lags included in augmentation T_obs <- 50 # number of observations per sample reps <- 10000 # number of observations per sample error_variance <- 5 # variance of innovations alpha_s <- 0.05 # significance level types <- c("none", "drift", "trend") names(types) <- c("Type - 0", "Type - Mu", "Type - Tau") rejection_matrix <- matrix(NA, ncol = 3, nrow = length(p), dimnames = list(as.character(p), names(types))) set.seed(1234) for(type in types){ # type = "drift" rejection_matrix[ ,names(types[which(type == types)])] <- getDFSizesForDGP( p = p, type = type, T_obs = T_obs, reps = reps, DGP = replicate( reps , ar2.gauss(T = T_obs, alpha =c(0.6,0.4), y.start = c(0,0), sigma = sqrt(error_variance))), alpha_s = alpha_s) } mat2_50 <- rejection_matrix