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
### commands to plot Annual flow stats

#### 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.

## boxplot the May pdfs - commands at the end of this file
### Or you can write out the May PDFs from simulations and 
### use/modify the boxpdfs.R file

# This code is to model a Seasonal AR(1) or Thomas-Fiering Model
# 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-2016.txt"),ncol=13,byrow=T)
WYmonflow=test[,2:13]

nyrs = dim(test)[1]

nyrs1=nyrs-1

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

WYmonflow=WYmonflow/10^6			 # convert to MAF
  
# get the co-efficients of the Thomas Fiering Model 
coef1=1:12
coef2=1:12
for(j in 1:12){
j1=j-1
if(j == 1)j1=12 
if(j == 1)coef1[j]=(cor(WYmonflow[2:nyrs,j], WYmonflow[1:nyrs1,j1])) 
if(j > 1)coef1[j]=cor(WYmonflow[1:nyrs,j], WYmonflow[1:nyrs,j1])
coef2[j]=sqrt((var(WYmonflow[,j])) * (1. - coef1[j]*coef1[j]))
coef1[j]=coef1[j]*sqrt(var(WYmonflow[,j])/var(WYmonflow[,j1]))
}

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

nsim=250	# number of simulations 

# get the co-efficients of the Thomas Fiering Model 
coef1=1:12
coef2=1:12
for(j in 1:12){
j1=j-1
if(j == 1)j1=12 
if(j == 1)coef1[j]=(cor(WYmonflow[2:nyrs,j], WYmonflow[1:nyrs1,j1])) 
if(j > 1)coef1[j]=cor(WYmonflow[1:nyrs,j], WYmonflow[1:nyrs,j1])
coef2[j]=sqrt((var(WYmonflow[,j])) * (1. - coef1[j]*coef1[j]))
coef1[j]=coef1[j]*sqrt(var(WYmonflow[,j])/var(WYmonflow[,j1]))
}

#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


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=1:nmons

i=round(runif(1,1,nyrs))
xsim[1]=WYmonflow[i,1]
xprev=xsim[1]

for(i in 2:nmons){
j=i %% 12
if(j == 0)j=12
j1=j-1

if(j == 1)j1=12

x1=xprev-mean(WYmonflow[,j1])
x2=coef2[j]*rnorm(1,0,1)

xsim[i]=mean(WYmonflow[,j]) + x1*coef1[j] + x2

xprev=xsim[i]}

#put the array inot a matrix 
simdismon=matrix(xsim,nrow=12)
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


### light grey spaggeti plots of simulated PDFs and 
### historic PDF in solid red

plot(xeval,xdensityorig,lwd=3,col="red",ylim=range(simpdf,xdensityorig))
for(i in 1:nsim)lines(xeval,simpdf[i,],col='grey')
lines(xeval,xdensityorig,lwd=3,col="red")

#### OR

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