{
WYmonflow_matrix(scan("leesferry-mon.dat"),ncol=13,byrow=T)
y_WYmonflow[1:76,2:13] #76x12
WYmonflow_WYmonflow[1:76,2:13]
annlee_1:76
for(i in 1:76)annlee[i]_sum(WYmonflow[i,1:12])
x_annlee

prop_matrix(0,nrow=76,ncol=12)
for(i in 1:76)prop[i,1:12]_y[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,100,13)
arstdev_matrix(0,100,13)
arcor_matrix(0,100,13)
arskw_matrix(0,100,13)

for(k in 1:100){


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

index_runif(76,1,76)
index_round(index)

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

#***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) {
                arcor[k,j]_cor(simdismon[,j],simdismon[,j-1])
        }
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,]))
        }
        armean[k,13]_mean(arann)
        arstdev[k,13]_stdev(arann)
        arskw[k,13]_sum((arann-armean[k,13])^3)
        arskw[k,13]_arskw[k,13]/76
        arskw[k,13]_arskw[k,13]/arstdev[k,13]^3

        arcor[k,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){
        obscor[i]_cor(WYmonflow[,i], WYmonflow[,i-1])
}
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])



monbox_function(data)
{

#data is the matrix containing the monthly statistic - 100 rows and 13
columns
months_c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Ann")

xs_1:13
zz_boxplot(split(t(data),xs),plot=F,cex=1.0)
zz$names_rep("",length(zz$names))
z1_bxp(zz,ylim=range(data),xlab="",ylab="",style.bxp="old",cex=1.25)
axis(1,at=z1,labels=months,cex=1.25)
points(z1,obs,pch=17,cex=1.25)
lines(z1[1:12],obs[1:12],lty=1)
#title(main="Boxplots of monthly mean",cex=1.0)
}

#******Prints .ps files for each stastistic to the current
#directory*****
#******I use the lpr -P printer filename command to print each
#boxplot***

obs_obsmean
monbox(armean)
title(main="Boxplots of Simulated Mean of Flows for Dissaggregated
Annual AR(1) ",cex=.7)
dev.copy(postscript, file="armean.ps")
dev.off()
obs_obsstdev
monbox(arstdev)
title(main="Boxplots of Simulated Standard Deviation of Flows AR(1)
",cex=1.0)
dev.copy(pscript, file="arstdev.ps")
dev.off()
obs_obscor
monbox(arcor)
title(main="Boxplots of Simulated Lag(1) Correlation of Flows AR(1)
",cex=1.0)
dev.copy(pscript, file="arcor.ps")
dev.off()
obs_obsskw
monbox(arskw)
title(main="Boxplots of Simulated Skew of Flows AR(1) ",cex=1.0)
dev.copy(pscript, file="arskw.ps")
dev.off()

}

}