lag1knn = function(annflow,nsims){

### This function implements K-NN lag-1 resampling
### It lags the input data vector by 1 into 
## vectors Xt1- the previous time and
### Xt the current time
## neighbors are found on Xt1 and simulated from Xt

# length of the annual flow
N=length(annflow)
N1 = N-1

Xt = annflow[2:N]
Xt1 = annflow[1:N1]


## number of simulations desired is nsim
## each simulation will be of length N

nsim = nsims	# you can change this

simflow = matrix(0,ncol=nsim, nrow=N)

## the weight metric to do the K-NN resampling
K=round(sqrt(N1))
W=1:K
W=1/W
W=W/sum(W)
W=cumsum(W)


for(i in 1:nsim){

#pick a flow at random to start
simflow[1,i] = sample(Xt1,1)
xprev = simflow[1,i]

 for(j in 2:N) {

# find the dist of xprev to all the N-1 years in the historical record
# order the distances
xdist = order(abs(xprev - Xt1))
	
#select one of the nearest neighbor with the weight function
# above.

		xx=runif(1,0,1)
		xy=c(xx,W)
		xx=rank(xy)
		index=xdist[xx[1]]
		
# the successor value of sampled index, which is in vector Xt		
simflow[j,i] = Xt[index]	
xprev = Xt[index]
}

}

# write out the simulations..



out=c()
out$sims = simflow

as.list(out)

}