# This routine is for doing Trend Analysis on a time series based on moving windows of a fixed width
# and Plotting relative to overall median (standardized to 1). The confidence interval
# is based on bootstrapping (with replacement) the entire series and then moving window through
# the resampled series and recalculating median for each window segment.  Did not simply
# bootstrap the series median because that would not account for differences in numbers
# of samples in the sliding window. This script does not move by equal steps but instead uses 
# existing time(txa) values.

# STEP 1:
# Read in data and options
# txa = times for a series(first column), tya=y values for series, 
# zz = which column in dataframe to use for y
# the total length of time series
#NOTE: Script assumes dates are in ascending order regardless if AD or BP!!

trender3 <-function(data, win=300,zz=2,plt=F,title="",nruns=1,p=.95) {
# Forces dataset to be a matrix and then takes memory and reads in data
	dataa <- as.matrix(data)
	options(object.size = 1000000000)
		txa_as.vector(dataa[,1][!is.na(dataa[,1])])
		ta_as.vector(dataa[,zz][!is.na(dataa[,1])])
tya_ta
# STEP 2:  
# Calculate the number of rows in output matrix depending on wind and steps
winhalf_win/2
rowss_length(txa[txa>(min(txa)+winhalf) & txa<(max(txa)-winhalf)])
TREout_matrix(NA, nrow=rowss, ncol=3)
ciout <- matrix(NA, nrow = rowss, ncol =3)

## Calculate the Global Median
GlobMed_median(tya)

## Identify first sample point that in range (e.g. minimum time value + 1/2 window)
TREout[,3]_1
mincut_min(txa)+winhalf
maxcut_max(txa)-winhalf
row_length(txa)
k_0
for (i in 1:row) {	
if (txa[i]>=mincut & txa[i]<=maxcut){
	k_k+1
	winmed_txa[i]
	winmax_winmed+winhalf
	winmin_winmed-winhalf
	TREout[k,1]_winmed
	TREout[k,2]_median(tya[txa>winmin & txa<winmax])/GlobMed
}
}
#Starts the bootstrapping if nruns>1
if (nruns> 1) {
medconf_matrix(NA, nrow=rowss, ncol=nruns)
for (i in 1:nruns){
tya_as.matrix(sample(ta,replace=T))
k_0
for (j in 1:row) {	
if (txa[j]>=mincut & txa[j]<=maxcut){
	k_k+1
	winmed_txa[j]
	winmax_winmed+winhalf
	winmin_winmed-winhalf
	medconf[k,i]_median(tya[txa>winmin & txa<winmax])/GlobMed
}
#Closes rowss loop
}
#Bracket below closes the actual bootstrapping loop
}
ciout[,1]_TREout[,1]
for (i in 1:rowss)  {
   	 ciout[i,2] <- quantile(medconf[i,-1], probs=((1-p)/2), na.rm = T) 
     ciout[i,3] <- quantile(medconf[i,-1], probs=(p+((1-p)/2)),na.rm=T)
					}
#Bracket below closes the if nruns>1 loop
}
if(plt) {
if (nruns>1) {
plot(TREout[,1], TREout[,2],type='b',xlab = "Years BP", ylab = "Median",pch=16, ylim=c((min(min(ciout[,2]),min(TREout[,2]))),(max(max(ciout[,3]),max(TREout[,2])))))
lines(TREout[,1], TREout[,3],lwd=3)
lines(ciout[,1],ciout[,2],lty=2)
lines(ciout[,1],ciout[,3],lty=2)
}
if (nruns < 2) {
plot(TREout[,1], TREout[,2],type='b',xlab = "Years BP", ylab = "Median",pch=16)
lines(TREout[,1], TREout[,3],lwd=3)
}
}
#optional
#res1<-list(Time=TREout[,1], WinMed = TREout[,2],GlobalMedian=TREout[,3])
res1<-list(Time=TREout[,1], WinMed = TREout[,2],lower=ciout[,2], upper=ciout[,3])
#exportData(res1,"out.xls",type="EXCEL")	
#Closes Function
}