##
    ## Fit a multivariate K-NN lag-1 model.
    ##


## Suppose the monthly streamflow matrix at Lees Ferry and Green River are
## in lf and gr, respectively.

## Lees Ferry   2006-1906
lf = matrix(scan("LeesFerry-monflows-1906-2006.txt"), ncol=13, byrow=T)
lf = lf[,2:13]
lf = lf*0.0004690502	# convert to cms

lfmean = colMeans(lf)
lfsd = apply(lf,2,sd)

#lfscale = scale(lf)
lfscale = t( (t(lf) - lfmean) / lfsd)

## Green River Ut
gr = matrix(scan("GreenRiver-UT-monflows-1906-2006.txt"), ncol=13, byrow=T)
gr = gr[,2:13]
gr = gr*0.0004690502     #convert to cms

grmean = colMeans(gr)
grsd = apply(gr,2,sd)

#grscale = scale(gr)
grscale = t( (t(gr) - grmean) / grsd)

    ## Bind the data into a 2 x 101(NumYears) X 12(months) matrix
library(abind)

    z = abind(lfscale, grscale, along=0)       # binds matrices.

nyrs = length(z[1,,1])   # number of years
nyrs1=nyrs-1

########### simulate..
nsims=250   # number of simulation..
sims = array(0,dim=c(nsims,dim(z)))

### sims is a four dimensional array..
### first dimension is the simulation number
### second dimension is the location 1, for lees ferry and 2 for green river
### third dimension is rows - i.e., the years
### fourth is the columns - i.e, months

K = sqrt(nyrs)
W=1:K
W=1/W
W=W/sum(W)
W=cumsum(W)



for(i in 1:nsims){

this.sim = array(0,dim=c((nyrs*12),2))


## Start with a previ Dec value..
sample.num = round(runif(1,1,nyrs))
Z.prev = as.matrix(z[,sample.num,12])

for(j in 1:(nyrs*12)){

imon = j%%12    ## the month we are simulating
if(imon ==0)imon = 12

 
 
     if(imon == 1) {
	imon1 = 12
	prev.vector = rbind(t(z[,1:nyrs1,imon1]))
	prev.vector = rbind(t(Z.prev),prev.vector) 
		N = nyrs1}
		else
		    {
		imon1 = imon-1
	prev.vector = rbind(t(z[,,imon1])) 
	prev.vector = rbind(t(Z.prev),prev.vector)
		N = nyrs}
		

        ordered.distance = order(as.matrix(dist(prev.vector))[1,2:(N+1)])

		xx=runif(1,0,1)
		xy=c(xx,W)
		xx=rank(xy)
		i1=ordered.distance[xx[1]]
		if(imon == 1)i1=i1+1

		Z.t = t(z[,i1,imon])

	this.sim[j,] = Z.t

## set the previous value and continue.
	Z.prev = t(Z.t)
			}
		

## Put back the mean and std dev..  and 
## Store the simulation back in the big sim array
this.simlf = matrix(this.sim[,1],ncol=12,byrow=T)
this.simgr = matrix(this.sim[,2],ncol=12,byrow=T)

this.simlf = t(t(this.simlf) * lfsd + lfmean)

this.simgr = t(t(this.simgr) * grsd + grmean)

sims[i,1,,] = this.simlf
sims[i,2,,] = this.simgr

}

#### You can grab each simulation for each location separately and compute
### the stats and do the boxplots. You can use the codes you have from 
### HW 2.