source("http://civil.colorado.edu/~balajir/CVEN6833/R-sessions/session3/files-4HW3/skew.R")
library(sm)


##### Things you have to add/develop to this######
### Annual flow PDF

##### you can also compute drought/surplus statistics if you wish...
#### Drought and Surplus Stats for the annual flow
### First compute the annual flow - i.e. sum the monthly flows
### Select Median of the observed annual flow as the threshold
### Annual flows above this is surplus and below is drought
### Sum the excess from the median for surplus and deficit from the median
### for the deficit.
## you can also compute the average/max length of drought and surplus


### The matrices of statistics are written in separate files
#### You can use/modify the boxplots.R to create boxplots of the stats.
### Again feel free to modify/change these codes.


# This code is to simulate from an ARMA fitted to the data
#
# for monthly streamflow time series.

# you should HAVE THE DATE IN THIS MATRIX WYmonflow before you 
# run this code..

#the data is in WYmonflow - that contains number of years as
# rows and 12 columns (Jan through Dec)


test=matrix(scan("LeesFerry-monflows-1906-2006.txt"),ncol=13,byrow=T)
WYmonflow=test[,2:13]

#WYmonflow=WYmonflow*0.0004690502      # conv ert to cms or

WYmonflow=WYmonflow/10^6			 # convert to MAF
  	
### Scale the data..
fmean = apply(WYmonflow,2,mean)		#monthly mean
fsdev = apply(WYmonflow,2,sd)		#monthly standard deviation

# get the data into one long array...
# scale the flow data.. - i.e. remove the monthly mean and divide #by monthly standard deviation

#X = t(t(WYmonflow) - fmean ) / fsdev

X = t((t(WYmonflow) - fmean) / fsdev)

X=array(t(X))		#standardized anomalies..

## Fit various ARMA models
p = 2    #order of AR component
q = 1    # order of MA component
d = 0
zz=arima0(X,order=c(p,0,q))    

modelaic = zz$aic  #this gives the AIC
#### Repeat the above for a number of models Do p = 1,2,3 and
### q = 1,2,3 and their combinations. The best model is the one
### with the least AIC

################### Simulate
#say ARMA(2,1)

nsim=250	# number of simulations 




#initialize the matrices that store the statistics
# the first row contains the statistic of the historical
#data. 

armean=matrix(0,nsim,12)	
arstdev=matrix(0,nsim,12)
arcor=matrix(0,nsim,12)
arskw=matrix(0,nsim,12)
armax=matrix(0,nsim,12)
armin=matrix(0,nsim,12)

## Store the May simulations and array for the PDF
maysim=1:nyrs
simpdf=matrix(0,nrow=nsim,ncol=nevals)

## Points to evaluate the May pdf

May = WYmonflow[,5]
xeval=seq(min(May)-0.25*sd(May),max(May)+0.25*sd(May),length=100)
nevals=length(xeval)



nyrs=length(WYmonflow[,1])	#number of years of data for each month
nyrs1=nyrs-1


# get the data into one long array...
# scale the flow data.. - i.e. remove the monthly mean and divide #by monthly standard deviation


#X = t(t(WYmonflow) - fmean ) / fsdev

X = t((t(WYmonflow) - fmean) / fsdev)

X=array(t(X))		#standardized anomalies..
zz=arima0(X,order=c(p,0,q))    


for(k in 1:nsim){
nmons=nyrs*12	#number of values to be generated
		#nyrs is the number of years and 12 is the number of
		#months.
xsim=arima.sim(n=nmons,list(ar=c(zz$coef[1:p]),ma=c(zz$coef[(p+1):(p+q)])), sd=sqrt(zz$sigma2))+mean(X)
#xsim=arima.sim(n=nmons,list(order=c(p,d,q),ar=c(zz$coef[1:2]),ma=c(zz$coef[3])), sd=sqrt(zz$sigma2))+mean(X)

#put the array inot a matrix 
simdismon=matrix(xsim,nrow=12)



# back standardize
simdismon = simdismon *fsdev + fmean

simdismon=t(simdismon)

	maysim = simdismon[,5]
	simpdf[k,]=sm.density(maysim,eval.points=xeval,display="none")$estimate

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

for(j in 1:12){

                armean[k,j]=mean(simdismon[,j])
                armax[k,j]=max(simdismon[,j])
                armin[k,j]=min(simdismon[,j])
                arstdev[k,j]=sd(simdismon[,j])
		arskw[k,j]=skew(simdismon[,j])

        }

# correlation between one month to another..

for(j in 2:12) {
j1=j-1
                arcor[k,j]=cor(simdismon[,j],simdismon[,j1])
        }
arcor[k,1]=cor(simdismon[1:nyrs1,12],simdismon[2:nyrs,1])

}

#*************************
#Compute statistics from the Historical DAta.

obsmean=1:12
obsstdev=1:12
obscor=1:12
obsskw=1:12
obsmax=1:12
obsmin=1:12

for(i in 1:12){

        obsmax[i]=max(WYmonflow[,i])
        obsmin[i]=min(WYmonflow[,i])
        obsmean[i]=mean(WYmonflow[,i])
        obsstdev[i]=sd(WYmonflow[,i])
	obsskw[i]=skew(WYmonflow[,i])

}
for(i in 2:12){
i1=i-1
        obscor[i]=cor(WYmonflow[,i], WYmonflow[,i1])
}
obscor[1]= cor(WYmonflow[1:nyrs1,12], WYmonflow[2:nyrs,1])



### bind the stats of the observed data at the top..
        armean=rbind(obsmean,armean)
        arstdev=rbind(obsstdev,arstdev)
        arskw=rbind(obsskw,arskw)
        arcor=rbind(obscor,arcor)
        armax=rbind(obsmax,armax)
        armin=rbind(obsmin,armin)

#write the stats in separate files..

write(t(arcor),file="arcorrs1",ncol=12)
write(t(armean),file="armeans",ncol=12)
write(t(arstdev),file="arstdevs",ncol=12)
write(t(arskw),file="arskews",ncol=12)
write(t(armax),file="armax",ncol=12)
write(t(armin),file="armin",ncol=12)



#boxplots of the May PDF....
## Obtain the PDF of the observed May flows

zz=sm.density(May,eval.points=xeval,display="none")
xdensityorig=zz$estimate


par(mfrow=c(1,1))
xs = 1:length(xeval)
zz=boxplot(split(t(simpdf),xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(simpdf,xdensityorig),xlab="May flow MAF",ylab="PDF",cex=1.25)

z2=1:6
n1=1:6
z2[1]=z1[1]
z2[2]=z1[20]
z2[3]=z1[40]
z2[4]=z1[60]
z2[5]=z1[80]
z2[6]=z1[100]

n1[1]=xeval[1]
n1[2]=xeval[20]
n1[3]=xeval[40]
n1[4]=xeval[60]
n1[5]=xeval[80]
n1[6]=xeval[100]
n1=round(n1,dig=0)
n1=as.character(n1)

axis(1,at=z2,labels=n1,cex=1.00)
lines(z1,xdensityorig,lty=2,lwd=2,col="red")