{
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)                         # MD                 

}

#*************************
#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)

}