S-PLUS : Copyright (c) 1988, 1999 MathSoft, Inc.
S : Copyright Lucent Technologies, Inc.
Version 5.1 Release 1 for Sun SPARC, SunOS 5.5 : 1999 
Working data will be in /bechtel/users1/faculty/balajir/MySwork 
Hello, Welcome to S-Plus
Library SM Attached
> source("knn-locfit.s")
> {
	WYmonflow <- matrix(scan("cisco-water-year.dat"), ncol = 12, byrow = T)
	febsim <- matrix(0, ncol = 100, nrow = 76)
	aprsim <- matrix(0, ncol = 100, nrow = 76)
	junsim <- matrix(0, ncol = 100, nrow = 76)
	augsim <- matrix(0, ncol = 100, nrow = 76)
	octsim <- matrix(0, ncol = 100, nrow = 76)
	annlee <- 1:76
	for(i in 1:76)
		annlee[i] <- sum(WYmonflow[i, 1:12])
	x <- annlee
	#fit lowess and get the residuals series..
	resids <- matrix(0, nrow = 75, ncol = 12)
	for(i in 1:12) {
		i1 <- i - 1
		if(i == 1) {
			zz <- locfit(WYmonflow[2:76, 1] ~ WYmonflow[1:75, 12])
			x1l <- zz
			resids[, i] <- residuals(zz)
		}
		else {
			zz <- locfit(WYmonflow[1:75, i] ~ WYmonflow[1:75, i1])
			if(i == 2)
				x2l <- zz
			if(i == 3)
				x3l <- zz
			if(i == 4)
				x4l <- zz
			if(i == 5)
				x5l <- zz
			if(i == 6)
				x6l <- zz
			if(i == 7)
				x7l <- zz
			if(i == 8)
				x8l <- zz
			if(i == 9)
				x9l <- zz
			if(i == 10)
				x10l <- zz
			if(i == 11)
				x11l <- zz
			if(i == 12)
				x12l <- zz
			resids[, i] <- residuals(zz)
		}
	}
	armean <- matrix(0, 101, 13)
	arstdev <- matrix(0, 101, 13)
	arcor <- matrix(0, 101, 13)
	arskw <- matrix(0, 101, 13)
	mondiffs <- matrix(0, 101, 911)
	armax <- matrix(0, 101, 13)
	armin <- matrix(0, 101, 13)
	### drought stats..
	mxsp <- 1:101
	mxdef <- 1:101
	maxs <- 1:101
	maxd <- 1:101
	index <- 1:76
	kk <- sqrt(75)
	kk <- round(kk)
	W <- 1:kk
	W <- 1/W
	W <- W/sum(W)
	W <- cumsum(W)
	frac <- 1
	th <- frac * mean(array(t(WYmonflow[1:76, 1:12])))
	th <- quantile(array(t(WYmonflow[1:76, 1:12])), 0.5)
	#	Get the first year..
	for(k in 1:100) {
		xsim <- 1:912
		i <- round(runif(1, 1, 75))
		xsim[1] <- WYmonflow[i, 1]
		xprev <- xsim[1]
		for(i in 2:912) {
			j <- i %% 12
			if(j == 0)
				j <- 12
			j1 <- j - 1
			if(j == 1) {
				xx <- abs(xprev - WYmonflow[1:75, 12])
			}
			else {
				xx <- abs(xprev - WYmonflow[1:75, j1])
			}
			xz <- order(xx)
			xz <- xz[1:kk]
			xx <- runif(1, 0, 1)
			xy <- c(xx, W)
			xx <- rank(xy)
			i1 <- xz[xx[1]]
			if(j == 1)
				xm <- predict(x1l, xprev)
			if(j == 2)
				xm <- predict(x2l, xprev)
			if(j == 3)
				xm <- predict(x3l, xprev)
			if(j == 4)
				xm <- predict(x4l, xprev)
			if(j == 5)
				xm <- predict(x5l, xprev)
			if(j == 6)
				xm <- predict(x6l, xprev)
			if(j == 7)
				xm <- predict(x7l, xprev)
			if(j == 8)
				xm <- predict(x8l, xprev)
			if(j == 9)
				xm <- predict(x9l, xprev)
			if(j == 10)
				xm <- predict(x10l, xprev)
			if(j == 11)
				xm <- predict(x11l, xprev)
			if(j == 12)
				xm <- predict(x12l, xprev)
			xsim[i] <- xm + resids[i1, j]
			#print(c(i,j,xprev,xm,resids[i1,j]))
			xprev <- xsim[i]
		}
		simdismon <- matrix(xsim, nrow = 12)
		simdismon <- t(simdismon)
		junsim[1:76, k] <- simdismon[, 9]
		febsim[1:76, k] <- simdismon[, 5]
		aprsim[1:76, k] <- simdismon[, 7]
		augsim[1:76, k] <- simdismon[, 11]
		octsim[1:76, k] <- simdismon[, 1]
		#***calculate basic statistics****************
		k1 <- k + 1
		for(j in 1:12) {
			armean[k1, j] <- mean(simdismon[, j])
			armax[k1, j] <- max(simdismon[, j])
			armin[k1, j] <- min(simdismon[, j])
			arstdev[k1, j] <- stdev(simdismon[, j])
			arskw[k1, j] <- sum((simdismon[, j] - armean[k1, j])^
				3)
			arskw[k1, j] <- arskw[k1, j]/76
			arskw[k1, j] <- arskw[k1, j]/arstdev[k1, j]^3
		}
		for(j in 2:12) {
			j1 <- j - 1
			arcor[k1, j] <- cor(simdismon[, j], simdismon[, j1])
		}
		arcor[k1, 1] <- cor(simdismon[1:75, 12], simdismon[2:76, 1])
		arann <- 1:76
		for(i in 1:76) {
			arann[i] <- (sum(simdismon[i,  ]))
		}
		armean[k1, 13] <- mean(arann)
		arstdev[k1, 13] <- stdev(arann)
		armax[k1, 13] <- max(arann)
		armin[k1, 13] <- min(arann)
		arskw[k1, 13] <- sum((arann - armean[k1, 13])^3)
		arskw[k1, 13] <- arskw[k1, 13]/76
		arskw[k1, 13] <- arskw[k1, 13]/arstdev[k1, 13]^3
		arcor[k1, 13] <- cor(arann[1:75], arann[2:76])
		## mon diffs..
		xx <- array(t(simdismon[1:76, 1:12]))
		xy <- diff(xx)
		mondiffs[k1,  ] <- xy
		###drought statistics
		#frac_1
		#th_frac*mean(array(t(WYmonflow[1:76,1:12])))
		#th_quantile(array(t(WYmonflow[1:76,1:12])),0.5)
		A <- array(t(simdismon))
		xx <- A
		Y <- A
		xx[xx <= th] <- 0
		xx[xx > th] <- 1
		n <- length(xx)
		x1 <- xx
		if(xx[1] == 1)
			x1 <- c(0, x1)
		if(xx[1] == 1)
			Y <- c(0, Y)
		if(xx[n] == 1)
			x1 <- c(x1, 0)
		if(xx[n] == 1)
			Y <- c(Y, 0)
		n <- length(x1[x1 == 1])
		xy <- order(x1)[1:(length(x1) - n)]
		y3 <- xy
		x3 <- diff(xy) - 1
		x3 <- x3[x3 > 0]
		mxsp[k1] <- max(x3)
		#Longest Surplus        LS              
		#Maximum surplus
		nsp <- length(x3)
		surp <- 1:nsp
		nsp <- length(y3)
		#print(c(k1,mxsp[k1],mxdef[k1],i1,i2,nsp)) 
		isp <- 0
		for(i in 1:(nsp - 1)) {
			i1 <- y3[i] + 1
			i2 <- y3[i + 1] - 1
			if(i2 >= i1) {
				isp <- isp + 1
				surp[isp] <- sum(Y[i1:i2])
			}
		}
		#print(c(k1,isp,nsp,length(x3)))
		maxs[k1] <- max(surp)
		# MS
		# for longest drought
		xx <- A
		Y <- A
		xx[xx <= th] <- 1
		xx[xx > th] <- 0
		n <- length(xx)
		x1 <- xx
		if(xx[1] == 1)
			x1 <- c(0, x1)
		if(xx[1] == 1)
			Y <- c(0, Y)
		if(xx[n] == 1)
			x1 <- c(x1, 0)
		if(xx[n] == 1)
			Y <- c(Y, 0)
		n <- length(x1[x1 == 1])
		xy <- order(x1)[1:(length(x1) - n)]
		y3 <- xy
		x2 <- diff(xy) - 1
		x2 <- x2[x2 > 0]
		mxdef[k1] <- max(x2)
		#Longest Drought        LD
		#Maximum Deficit
		nsp <- length(x2)
		surp <- 1:nsp
		nsp <- length(y3)
		isp <- 0
		for(i in 1:(nsp - 1)) {
			i1 <- y3[i] + 1
			i2 <- y3[i + 1] - 1
			if(i2 >= i1) {
				isp <- isp + 1
				surp[isp] <- sum(Y[i1:i2])
			}
		}
		#print(c(k1,isp,nsp,length(x3)))
		maxd[k1] <- max(surp)
	}
	#*************************
	#Basic Data
	obsmean <- 1:13
	obsstdev <- 1:13
	obscor <- 1:13
	obsskw <- 1:13
	obsmax <- 1:13
	obsmin <- 1:13
	for(i in 1:12) {
		obsmax[i] <- max(WYmonflow[, i])
		obsmin[i] <- min(WYmonflow[, i])
		obsmean[i] <- mean(WYmonflow[, i])
		obsstdev[i] <- stdev(WYmonflow[, i])
		obsskw[i] <- sum((WYmonflow[, i] - obsmean[i])^3)
		obsskw[i] <- obsskw[i]/76
		obsskw[i] <- obsskw[i]/obsstdev[i]^3
	}
	for(i in 2:12) {
		i1 <- i - 1
		obscor[i] <- cor(WYmonflow[, i], WYmonflow[, i1])
	}
	obscor[1] <- cor(WYmonflow[1:75, 12], WYmonflow[2:76, 1])
	xx <- array(t(WYmonflow[1:76, 1:12]))
	xy <- diff(xx)
	mondiffs[1,  ] <- xy
	obsann <- 1:76
	for(i in 1:76) {
		obsann[i] <- (sum(WYmonflow[i,  ]))
	}
	obsmax[13] <- max(obsann)
	obsmin[13] <- min(obsann)
	obsmean[13] <- mean(obsann)
	obsstdev[13] <- stdev(obsann)
	obsskw[13] <- sum((obsann - obsmean[13])^3)
	obsskw[13] <- obsskw[13]/76
	obsskw[13] <- obsskw[13]/obsstdev[13]^3
	obscor[13] <- cor(obsann[1:75], obsann[2:76])
	armean[1, 1:13] <- obsmean[1:13]
	arstdev[1, 1:13] <- obsstdev[1:13]
	arskw[1, 1:13] <- obsskw[1:13]
	arcor[1, 1:13] <- obscor[1:13]
	armax[1, 1:13] <- obsmax[1:13]
	armin[1, 1:13] <- obsmin[1:13]
	#frac_1
	#th_frac*mean(array(t(WYmonflow[1:76,1:12])))
	#th_quantile(array(t(WYmonflow[1:76,1:12])),0.5)
	A <- array(t(WYmonflow[1:76, 1:12]))
	xx <- A
	Y <- A
	xx[xx <= th] <- 0
	xx[xx > th] <- 1
	n <- length(xx)
	x1 <- xx
	if(xx[1] == 1)
		x1 <- c(0, x1)
	if(xx[1] == 1)
		Y <- c(0, Y)
	if(xx[n] == 1)
		x1 <- c(x1, 0)
	if(xx[n] == 1)
		Y <- c(Y, 0)
	n <- length(x1[x1 == 1])
	xy <- order(x1)[1:(length(x1) - n)]
	y3 <- xy
	x3 <- diff(xy) - 1
	x3 <- x3[x3 > 0]
	mxsp[1] <- max(x3)
	#Longest Surplus        LS              
	#Maximum Surplus
	nsp <- length(x3)
	surp <- 1:nsp
	nsp <- length(y3)
	isp <- 0
	for(i in 1:(nsp - 1)) {
		i1 <- y3[i] + 1
		i2 <- y3[i + 1] - 1
		if(i2 >= i1) {
			isp <- isp + 1
			surp[isp] <- sum(Y[i1:i2])
		}
	}
	maxs[1] <- max(surp)
	# MS
	# for longest drought
	xx <- A
	Y <- A
	xx[xx <= th] <- 1
	xx[xx > th] <- 0
	n <- length(xx)
	x1 <- xx
	if(xx[1] == 1)
		x1 <- c(0, x1)
	if(xx[1] == 1)
		Y <- c(0, Y)
	if(xx[n] == 1)
		x1 <- c(x1, 0)
	if(xx[n] == 1)
		Y <- c(Y, 0)
	n <- length(x1[x1 == 1])
	xy <- order(x1)[1:(length(x1) - n)]
	y3 <- xy
	x2 <- diff(xy) - 1
	x2 <- x2[x2 > 0]
	mxdef[1] <- max(x2)
	#Longest Drought        LD
	#Maximum deficit
	nsp <- length(x2)
	surp <- 1:nsp
	nsp <- length(y3)
	isp <- 0
	for(i in 1:(nsp - 1)) {
		i1 <- y3[i] + 1
		i2 <- y3[i + 1] - 1
		if(i2 >= i1) {
			isp <- isp + 1
			surp[isp] <- sum(Y[i1:i2])
		}
	}
	maxd[1] <- max(surp)
	# MD                 
	write(t(arcor), file = "arcorrs1", ncol = 13)
	write(t(armean), file = "armeans", ncol = 13)
	write(t(arstdev), file = "arstdevs", ncol = 13)
	write(t(arskw), file = "arskews", ncol = 13)
	write(t(armax), file = "armax", ncol = 13)
	write(t(armin), file = "armin", ncol = 13)
	write(t(mondiffs), file = "mondiffs", ncol = 911)
	droutstat <- cbind(mxsp, mxdef, maxs, maxd)
	write(t(febsim), file = "febsim", ncol = 100)
	write(t(aprsim), file = "aprsim", ncol = 100)
	write(t(junsim), file = "junsim", ncol = 100)
	write(t(augsim), file = "augsim", ncol = 100)
	write(t(octsim), file = "octsim", ncol = 100)
	write(t(droutstat), file = "drought-stats", ncol = 4)
}