Independence <- function(regressor, regressand, N=999, lag=1, Kernel=1, logistic=FALSE, q=c(.5, 1, 2, 4), epsilon=1e-3, trim.type=1, quantiles=c(.2, .8), msd=1) { #################################################################################### # # # Nonparametric Entropy-based Tests of Independence between Stochastic Processes # # Fernandes, Marcelo and Breno Neri, 2007 # # # # Department of Economics, Queen Mary, University of London # # m.fernandes@qmul.ac.uk http://www.econ.qmul.ac.uk/staff/fernandes.htm # # # # Department of Economics, New York University # # breno.neri@nyu.edu http://homepages.nyu.edu/~bpn207 # # # # January, 2008 - Version 3 # # # #################################################################################### # # # Parameter: Default: Description: # # # # N 999 Number of Bootstrap Samples # # lag 1 Lag in the serial independence test # # Kernel 1 (1) Gaussian, (2) Gamma, (3) Beta [only 1 for now] # # logistic FALSE Logistic transformation: TRUE or FALSE # # q c(.5, 1, 2, 4) Vector of entropic indices # # epsilon 1e-3 Avoid division by zero or bondary problems # # trim.type 1 (1) Quantiles or (2) Multiples of standard deviat. # # quantiles c(.2, .8) Vector of trimming quatiles # # msd 1 Multiples of standard deviations (for trimming) # # # #################################################################################### # # # 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. # # # #################################################################################### # Initializing library(tseries) nq <- length(q) test.statistic <- array(0, c(N+1, nq)) bds.test.statistic <- array(0, c(N+1, nq)) # Weigh Function if(trim.type==1) w <- function(x) x > quantile(x, quantiles[1]) & x < quantile(x, quantiles[2]) else if(trim.type==2) w <- function(x) abs(x-mean(x)) < msd*sd(x) else stop("trim.type must be 1 or 2!") # Bandwidths bandwidthSilverman <- function(x) bw.nrd(x) #CVgamma <- function(h) #{ # h <<- h # term1 <- integrate(gammakernel2, epsilon, 2, subdivisions=500, stop.on.error=FALSE)[[1]] # term2 <- numeric() # for(j in 1:n) term2[j] <- gammakernel(x[-j], x[j], h) # return(term1 - 2*mean(term2)) #} #CVbeta <- function(h) #{ # h <<- h # term1 <- integrate(betakernel2, epsilon, 1-epsilon, subdivisions=500, stop.on.error=FALSE)[[1]] # term2 <- numeric() # for(j in 1:n) term2[j] <- betakernel(x[-j], x[j], h) # return(term1 - 2*mean(term2)) #} #bandwidthgamma <- function(x) return(optim(c(solve(n<<-length(x))), CVgamma, method="L-BFGS-B", lower=epsilon)[[1]]) #bandwidthbeta <- function(x) return(optim(c(solve(n<<-length(x))), CVbeta, method="L-BFGS-B", lower=epsilon)[[1]]) # Gaussian Kernel Estimation GK <- function(u) exp(-u^2/2) / sqrt(2*pi) # Gaussian Kernel Gf <- function(x, x0, bw) mean(GK((x-x0)/bw)) / bw # Density at x0 Gfxy <- function(x, y, x0, y0, bwx, bwy) mean(GK((x-x0)/bwx) * GK((y-y0)/bwy)) / (bwx*bwy) # Density at (x0, y0) # Asymmetric Kernels Estimation #gammakernel <- function(x, x0, bwx) mean(dgamma(x, shape=x0/bwx+1, scale=bwx)) # Density at x0 using Gamma Kernel #betakernel <- function(x, x0, bwx) mean(dbeta(x, shape1=x0/bwx+1, shape2=(1-x0)/bwx+1)) # Density at x0 using Beta Kernel #gammakernel2 <- function(x0) #{ # y <- numeric() # for(i in 1:length(x0)) y[i] <- gammakernel(x, x0[i], h)^2 # return(y) #} #betakernel2 <- function(x0) #{ # y <- numeric() # for(i in 1:length(x0)) y[i] <- betakernel(x, x0[i], h)^2 # return(y) #} #gammakernel.xy <- function(x, y, x0, y0, bwx, bwy) mean(dgamma(x, shape=x0/bwx+1, scale=bwx) * dgamma(y, shape=y0/bwy+1, #scale=bwy)) # Density at (x0, y0) using Gamma Kernel #betakernel.xy <- function(x, y, x0, y0, bwx, bwy) mean(dbeta(x, shape1=x0/bwx+1, shape2=(1-x0)/bwx+1) * dbeta(y, #shape1=y0/bwy+1, shape2=(1-y0)/bwy+1)) # Density at (x0, y0) using Beta Kernel # Test Statistic r <- function(x, y) { W <- matrix(rep(wx*wy, nq), c(t-lag, nq)) # if(Kernel==1) # { for(i in 1:(t-lag)) W[i, q!=1] <- W[i, q!=1]*(1-(Gf(x, x[i], bwx)*Gf(y, y[i], bwy)/Gfxy(x, y, x[i], y[i], bwx, bwy))^(1-q[q!=1]))/(1-q[q!=1]) if(sum(q==1)>0) for(i in 1:(t-lag)) W[i, q==1] <- W[i, q==1]*log(Gfxy(x, y, x[i], y[i], bwx, bwy)/(Gf(x, x[i], bwx)*Gf(y, y[i], bwy))) # } else if(Kernel==2) # { # bwx <- bandwidthgamma(x) # bwy <- bandwidthgamma(y) # for(i in 1:(t-lag)) W[i, q!=1] <- W[i, q!=1]*(1-(gammakernel(x, x[i], bwx)*gammakernel(y, y[i], bwy)/gammakernel.xy(x, y, x[i], y[i], bwx, bwy))^(1-q[q!=1]))/(1-q[q!=1]) # if(sum(q==1)>0) for(i in 1:(t-lag)) W[i, q==1] <- W[i, q==1]*log(gammakernel.xy(x, y, x[i], y[i], bwx, bwy)/(gammakernel(x, x[i], bwx)*gammakernel(y, y[i], bwy))) # } else if(Kernel==3) # { ## bwx <- bandwidthbeta(x) ## bwy <- bandwidthbeta(y) # for(i in 1:(t-lag)) W[i, q!=1] <- W[i, q!=1]*(1-(betakernel(x, x[i], bwx)*betakernel(y, y[i], bwy)/betakernel.xy(x, y, x[i], y[i], bwx, bwy))^(1-q[q!=1]))/(1-q[q!=1]) # if(sum(q==1)>0) for(i in 1:(t-lag)) W[i, q==1] <- W[i, q==1]*log(betakernel.xy(x, y, x[i], y[i], bwx, bwy)/(betakernel(x, x[i], bwx)*betakernel(y, y[i], bwy))) # } else stop("Kernel must be 1, 2 or 3!") return(apply(W, 2, mean)/sqrt(var(wx)/(t-lag))) } # Hong and White Test Statistic hw <- function(x, y) { temp <- numeric() for(i in 1:(t-lag)) temp[i] <- log(Gqxy(x, y, x[i], y[i], bwx, bwy)/(Gq(x, x[i], bwx)*Gq(y, y[i], bwy))) return(mean(temp)) } # Quartic Kernel Estimation (for Hong and White) GQ <- function(u) (15/16)*(1-u^2)^2*(abs(u)<=1) Gq <- function(x, x0, bw) mean(GQ((x-x0)/bw)) / bw # Density at x0 Gqxy <- function(x, y, x0, y0, bwx, bwy) mean(GQ((x-x0)/bwx) * GQ((y-y0)/bwy)) / (bwx*bwy) # Density at (x0, y0) # Computations ols <- lm(regressand~regressor) x <- ols$resid t <- length(x) if(logistic) x <- (temp<-exp(x))/(1+temp) bwx <- bandwidthSilverman(x[1:(t-lag)]) bwy <- bandwidthSilverman(x[(1+lag):t]) wx <- w(x[1:(t-lag)]) wy <- w(x[(1+lag):t]) test.statistic[1,] <- r(x[1:(t-lag)], x[(1+lag):t]) hw.test.statistic <- hw(x[1:(t-lag)], x[(1+lag):t]) bds.test.statistic[1,] <- c(bds.test(x, 2, q*sd(x))[[1]]) # Bootstrap for(bs in 1:N+1) { u <- sample(x, replace=TRUE) y <- cbind(1, regressor)%*%ols$coef + u u <- lm(y~regressor)$resid wx <- w(u[1:(t-lag)]) wy <- w(u[(1+lag):t]) test.statistic[bs,] <- r(u[1:(t-lag)], u[(1+lag):t]) hw.test.statistic[bs] <- hw(u[1:(t-lag)], u[(1+lag):t]) bds.test.statistic[bs,] <- c(bds.test(u, 2, q*sd(u))[[1]]) } # Results print(paste("Bootstrap Samples: ", N)) print(paste("Lags: ", lag)) print(paste("Logistic Transformation: ", logistic)) if(trim.type==1) { print("Trimming Quantiles:") print(quantiles) } else print(paste("Multiples of Standard Deviation: ", msd)) print("BDS P-Value:") print(apply(bds.test.statistic>bds.test.statistic[1,], 2, mean)) print(paste("Hong & White P-Value:", mean(hw.test.statistic>hw.test.statistic[1]))) print("Fernandes & Neri P-Value:") print(apply(test.statistic>test.statistic[1,], 2, mean)) }