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