{
#WYmonflow_matrix(scan("leesferry-mon.dat"),ncol=13,byrow=T)
#WYmonflow_WYmonflow[1:76,2:13]

WYmonflow_matrix(scan("water-year-flow.dat"),ncol=12,byrow=T)
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

x1_WYmonflow[1:75,12]
x2_annlee[2:76]
database_cbind(x1,x2)
xm1_mean(x1)
xm2_mean(x2)
xs1_sqrt(var(x1))
xs2_sqrt(var(x2))
database_scale(database)

prop_matrix(0,nrow=76,ncol=12)
for(i in 1:76)prop[i,1:12]_WYmonflow[i,1:12]/x[i]

#******Fit the historic data to an AR(1) model***********
obsmean_mean(x)
obsstdev_sqrt(var(x))

y_scale(x)
mody_arima.mle(y, model=list(order=c(1,0,0)))

#******Simulates AR(1) annually then performs disaggregation to
#monthly*******

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)

### drought stats..
mxsp_1:101
mxdef_1:101
maxs_1:101
maxd_1:101

index_1:76

kk_sqrt(76)
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)

for(k in 1:100){


#zy_rnorm(76,0,sqrt(mody$sigma2))
#arsim_arima.sim(76,model=list(ar=mody$model$ar), innov=zy)
#simann_arsim*obsstdev+obsmean

#do a K-NN on the annuals..
xx_round(runif(1,1,75))
simann[1]_annlee[xx]
for(i in 2:76){
i1_i-1
xx_abs(simann[i1]-annlee[1:75])
xz_order(xx)
xz_xz[1:kk]
xx_runif(1,0,1)
xy_c(xx,W)
xx_rank(xy)
i1_xz[xx[1]]
i1_i1+1
simann[i]_annlee[i1]}



#	Get the first year..

xx_abs(simann[1]-annlee)
xz_order(xx)
xz_xz[1:kk]
xx_runif(1,0,1)
xy_c(xx,W)
xx_rank(xy)
i1_xz[xx[1]]
index[1]_i1
xprev_prop[index[1],12]*simann[1]

kk_sqrt(76)/2
kk_round(kk)
W_1:kk
W_1/W
W_W/sum(W)
W_cumsum(W)

for(i in 2:76){
#x1_abs(simann[i]-annlee)
#x2_abs(xprev-WYmonflow[,12])

xx_(simann[i]-xm2)/xs2
x1_abs(xx-database[,2])
xx_(xprev-xm1)/xs1
x2_abs(xx-database[,1])
xx_sqrt(x1^2 + x2^2)
#xx_sqrt(0.25*x1^2 + 0.75*x2^2)

#xx_x1
xz_order(xx)
xz_xz[1:kk]
xx_runif(1,0,1)
xy_c(xx,W)
xx_rank(xy)
i1_xz[xx[1]]
index[i]_i1
index[i]_i1+1
xprev_prop[index[i],12]*simann[i]}

simdismon_matrix(99,76,12)
for(i in 1:76){
simdismon[i,1:12]_prop[index[i],1:12]*simann[i]
}

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])
                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)
        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


for(i in 1:12){

        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,]))
        }
        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]

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

}