## S-Plus code (coded by Eric Van den Berg and modified by JDG)
## File: Nonlinear.prg
##
## NOTE: rename file name "Resnick-vdBergNonlinear.txt"   
## 
## INPUT:
##   x     =  data file
##   alpha =  (>0) tail index of Pareto distribution  
##   h     =  time lag
##
## OUTPUT:
##   p-value of the test statistic 
##     T(h) = (N/log(N))^{1/\alpha} * |\hat{\rho}_{1}(h)-\hat{\rho}_{2}(h)|
##     for each value of h. Here, \hat{\rho}_{i}(h) (i=1,2) denotes the uncentered
##     heavy tailed sample ACF of two disjoint subsamples of length N, and defined 
##     by 
##        
##     \hat{\rho}_{i}(h) = \sum^{N-h} Y_{t}Y_{t+h}/\sum^{N}_{t=1} Y^{2}_{t},
##
##     as an estimator of 
##
##     \rho(h) = \sum^{\infty}_{j=0} c_{j}c_{j+h}/\sum^{\infty}_{j=0} c^{2}_{j},
##
##     if alpha<1, and the classical mean-centered sample ACF if alpha>1.
##
##     Assume Y_{t}=sum^{\infty}_{j=0} c_{j}Z_{t-j} (null hypothesis), where 
##     c_{j} satisfies a summability condition and Z_{t} are i.i.d. Under
##     the asumption of linearity, and Z_{t} having infinite variance, 
##     \hat{\rho}_{h) converges to \rho(h). However, for many nonlinear models 
##     the sample ACF converges to a nondegenerate random variable. Thm. 4.1 of
##     Resnick - Van den Berg (2000b) gives the asymptotic null distribution of
##     |\hat{\rho}_{1}(h)-\hat{\rho}_{2}(h)|.
##
## NOTE:  rstab returns random deviates from the stable family
##        of probability functions.
##
## References:
## Resnick, S. and Van den Berg, E. (2000a).
##   Sample correlation behavior for the heavy tailed general bilinear process.
##   Communication in Statistics: Stochastic Models, 16(2), 233-258.
##   DOI: 10.1080/15326340008807586
## Resnick, S. and Van den Berg, E. (2000b).
##   A test for nonlinearity of time series with infinite variance.
##   Extremes, 3(2), 145-172.
#############################################################
nonlin <- function(x,alpha,h) 
{
c <- sdisp(x,h,alpha) 
r <- dim(c) 
A <- stabmat(alpha,10000.,r[1.]) 
B <- A %*% c 
Y <- apply(abs(B),1.,max) 
m <- rhotest(x, alpha, h) 
p <- length(Y[Y>m])/10000. 
p 
}

sdisp <- function(x, h, alpha) 
{ 
n <- length(x) 
if(alpha < 4./3.) 
  {rn <- h + floor((n/2.)^0.499) } 
else 
  {rn <- h + floor((n/2.)^(2./alpha - 1.)) } 
if(alpha > 1) {rho <- c(1, acf(xx,rn+h,plot = F)$acf[, , 1])} 
  else { rho <- c(1., acf.htsid(xx,rn+h, 0.,1., F)$corr)} 
if(h > 1.) {rho <- c(rev(rho[2.:h]),rho)} 
sapply(1.:h, cdisp,rho,rn,h,alpha) 
}

cdisp <- function(l,rho,rn,h,alpha) 
{ 
y <- abs(rho[(h+l+1.):(h+l+rn)]+rho[(1.+(h-l)):(rn+(h-l))]-2.*rho[(h+1.):(h+rn)]*rho[h])^(alpha) 
y 
}

rhotest <- function(x,alpha,h) 
{ 
n <- floor(length(x)/2.) 
if(alpha > 1) {
  rho1 <- acf(x[1:n], h, plot = T)$acf[, , 1] 
  rho2 <- acf(x[(n + 1):(2*n)],h,plot=T)$acf[, , 1]} 
else { 
  rho1 <- acf.htsid(x[1.:n], h, 1., 1.,T)$corr 
  rho2 <- acf.htsid(x[(n + 1.):(2.*n)],h,1.,1.,T)$corr 
} 
rd <- (n/log(n))^(1./alpha) * abs(rho1-rho2) 
m  <- max(rd) 
# plot(rd,type = "h", ylim = c(0.,ceiling(1.1*m))) 
# abline(h = 0.,lty = 1.)
m 
}

acf.htsid <- function(x,p,lower,upper,plot= T) 
{ 
len  <- length(x) 
den  <- sum(x*x) 
max  <- 0. 
corr <- matrix(0.,p,1)
for(j in 1.:p) { 
   corr[j] <- sum(x[1.:(len-j)]*x[(1.+j):len])/den
} 
l <- max(abs(corr)) 
if(plot) {plot(1.:p,corr,type="h",ylim=c(-lower,upper),ylab="Heavy Tailed ACF",xlab="Lag") 
  abline(h = 0.,lty = 1.) 
} 
return(list(corr = corr,maxcorr=l)) 
}

stabmat <- function(alpha,N,rn) 
#  u = S_{j} and S_{j}^{'} variables
#  v = S_{0} and S_{0}^{'} variables
#  Output: (4.1) in Resnick and Van den Berg (2000b)
{ 
u <- (matrix(rstab(N*2.*rn, alpha,1.),nrow=N)+tan((pi*alpha)/2.))*(cos((pi*alpha)/2.)*gamma(1.-alpha))^(1./alpha) 
v <- (matrix(rstab(N*2.,alpha/2.,1.),ncol=2.)+tan((pi*alpha)/4.))*(cos((pi*alpha)/4.)*gamma(1.-alpha/2.))^(2./alpha)
u[,1.:rn]/v[,1.]-u[,(rn+1.):(2.*rn)]/v[,2.] 
}