#generates the independent variable and records its length
x <- 1:100
N <- length(x)

#creates the dependent variable with a weak trend
ymat <- rbind(.01*x+2,.01*x+2,.01*x+2)
M <- 3

#sets a seed so that results are reproducible and generates normal noise 
set.seed(15)
shared_noise <- rnorm(N)*sqrt(.5)

#adds randomly distributed noise and shared noise
for(i in 1:3){
    ymat[i,] <- ymat[i,] + rnorm(N)*sqrt(.5) + shared_noise
}

#examine the data
plot(x, ymat[1,])
points(x, ymat[2,],col='red')
points(x, ymat[3,],col='blue')

KendallS <- function(y){
    #Calculate Kendall's S statistic, S, from a series, y
    #Variables:
    #   y   input data series
    #   S   Kendall's S statisic
    #   N   length of y
    S <- 0
    N <- length(y)
    for(i in 1:(N-1)){
        S <- S + sum(sign(y[i]-y[(i+1):N]))
    }
    return(S)
}

#Svec is a list of Kendall S values for each dataset
Svec <- rep(0,M)
for(i in 1:M){
    Svec[[i]] <- KendallS(ymat[i,])
}

#store the average Kendall S value in Sm
Sm <- mean(Svec)

#records the variance for a single Kendall's S value
VarS <- (N*(N-1)*(2*N+5))/18

#calculate Z-scores and p-values for each Kendall’s S value
Zvec <- rep(0,M)
pvec <- rep(0,M)
for(i in 1:M){
    Zvec[[i]] <- Svec[[1]]/sqrt(VarS)
    pvec[[i]] <- 2*(1-pnorm(abs(Zvec[[i]])))
}

#store the variance assuming no covariance between datasets
VarSm_nocov <- N*(N-1)*(2*N+5)/18/M

#store the Z-score for Sm
Zfield_nocov <- Sm/(sqrt(VarSm_nocov)) 

#calculate p value
pfield_nocov <- 2*(1-pnorm(abs(Zfield_nocov)))

print(pfield_nocov)

#define crosscorr fuction
crosscorr <- function(y1,y2){
    #Variables: 
    #   y1          first input dataset
    #   y2          second input dataset
    #   mu1         mean of first dataset
    #   mu2         mean of second dataset
    #   sigma1      standard deviation of first dataset
    #   sigma2      standard deviation of second dataset
    mu1 <- mean(y1)
    mu2 <- mean(y2)
    sigma1 <- sqrt(var(y1))
    sigma2 <- sqrt(var(y2))
    return(mean((y1 - mu1) * (y2 - mu2)) / sigma1 / sigma2)
}

#compute the average cross-correlation coefficient
rhosum <- 0
for(k in 1:(M-1)){
    for(l in 1:(M-k)){
        rhosum <- rhosum + crosscorr(ymat[k,],ymat[l,])
    }
}
rhobar <- 2*rhosum/M/(M-1)

#compute variance, Z score, and field significance
VarSm <- (VarS/M)*(1+(M-1)*rhobar)

Zfield <- Sm/sqrt(VarSm)

pfield <- 2*(1-pnorm(abs(Zfield)))
print(pfield)















##### Repeated analysis using real data

#read Flow data from 3 rivers
x1 <- read.table('RockyRiver.txt',skip=36)[,5]
y1 <- read.table('RockyRiver.txt',skip=36)[,6]

x2 <- read.table('Cuyahoga.txt',skip=36)[,5]
y2 <- read.table('Cuyahoga.txt',skip=36)[,6]

x3 <- read.table('Vermillion.txt',skip=36)[,5]
y3 <- read.table('Vermillion.txt',skip=36)[,6]

#Record and count only years where all sites have data
x <- x1[is.element(x1,x2) & is.element(x1,x3)]
N <- length(x)

#create a matrix of flows for only shared years 
ymat = matrix(NA, 3, N)
ymat[1,] <- y1[is.element(x1,x) & is.element(x1,x)]
ymat[2,] <- y2[is.element(x2,x) & is.element(x2,x)]
ymat[3,] <- y3[is.element(x3,x) & is.element(x3,x)]


#examine the data
plot(x, ymat[1,])
points(x, ymat[2,],col='red')
points(x, ymat[3,],col='blue')

KendallS <- function(y){
    #Calculate Kendall's S statistic, S, from a series, y
    #Variables:
    #   y   input data series
    #   S   Kendall's S statistic
    #   N   length of y
    S <- 0
    N <- length(y)
    for(i in 1:(N-1)){
        S <- S + sum(sign(y[i]-y[(i+1):N]))
    }
    return(S)
}

#Svec is a list of Kendall S values for each dataset
Svec <- rep(0,M)
for(i in 1:M){
    Svec[[i]] <- KendallS(ymat[i,])
}

#store the average Kendall S value in Sm
Sm <- mean(Svec)

#records the variance for a single Kendall's S value
VarS <- (N*(N-1)*(2*N+5))/18

#calculate Z-scores and p-values for each Kendall’s S value
Zvec <- rep(0,M)
pvec <- rep(0,M)
for(i in 1:M){
    Zvec[[i]] <- Svec[[1]]/sqrt(VarS)
    pvec[[i]] <- 2*(1-pnorm(abs(Zvec[[i]])))
}

#store the variance assuming no covariance between datasets
VarSm_nocov <- N*(N-1)*(2*N+5)/18/M

#store the Z-score for Sm
Zfield_nocov <- Sm/(sqrt(VarSm_nocov)) 

#calculate p value
pfield_nocov <- 2*(1-pnorm(abs(Zfield_nocov)))

print(pfield_nocov)

#define crosscorr fuction
crosscorr <- function(y1,y2){
    #Variables: 
    #   y1          first input dataset
    #   y2          second input dataset
    #   mu1         mean of first dataset
    #   mu2         mean of second dataset
    #   sigma1      standard deviation of first dataset
    #   sigma2      standard deviation of second dataset
    mu1 <- mean(y1)
    mu2 <- mean(y2)
    sigma1 <- sqrt(var(y1))
    sigma2 <- sqrt(var(y2))
    return(mean((y1 - mu1) * (y2 - mu2)) / sigma1 / sigma2)
}

#compute the average cross-correlation coefficient
rhosum <- 0
for(k in 1:(M-1)){
    for(l in 1:(M-k)){
        rhosum <- rhosum + crosscorr(ymat[k,],ymat[l,])
    }
}
rhobar <- 2*rhosum/M/(M-1)

#compute variance, Z score, and field significance
VarSm <- (VarS/M)*(1+(M-1)*rhobar)

Zfield <- Sm/sqrt(VarSm)

pfield <- 2*(1-pnorm(abs(Zfield)))
print(pfield)