{
WYmonflow_matrix(scan("leesferry-mon.dat"),ncol=13,byrow=T)
WYmonflow_WYmonflow[1:76,2:13]
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_annlee[1:75]
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)
index_1:76

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

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

#	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_simann[1]


for(i in 1:38){
j1_(i-1)*2+1
j2_j1+1

xprev_simann[j1]
x1_abs(xprev-database[,1])
x2_abs(simann[j2]-database[,2])
xx_sqrt(x1^2 + x2^2)

#xx_x2
xz_order(xx)
xz_xz[1:kk]
xx_runif(1,0,1)
xy_c(xx,W)
xx_rank(xy)
i1_xz[xx[1]]
index[j1]_i1
index[j2]_i1+1
#if(index[i] > 76)index[i]_76}

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[,6]
febsim[1:76,k]_simdismon[,2]
aprsim[1:76,k]_simdismon[,4]
augsim[1:76,k]_simdismon[,8]
octsim[1:76,k]_simdismon[,10]

#***calculate basic statistics****************

for(j in 1:12){

                armean[k,j]_mean(simdismon[,j])
                arstdev[k,j]_stdev(simdismon[,j])
                arskw[k,j]_sum((simdismon[,j]-armean[k,j])^3)
                arskw[k,j]_arskw[k,j]/76
                arskw[k,j]_arskw[k,j]/arstdev[k,j]^3
        }
for(j in 2:12) {
j1_j-1
                arcor[k,j]_cor(simdismon[,j],simdismon[,j1])
        }
arcor[k,1]_cor(simdismon[1:75,12],simdismon[2:76,1])

arann_1:76
        for(i in 1:76){
                arann[i]_(sum(simdismon[i,]))
        }
k1_k+1
        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])
#cat(k)
}

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




        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]

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

}