## base code from Yanto

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

library(sm)

leesflow = read.table("LeesFerry-monflows-1906-2006.txt")
leesflow = leesflow[,-1]
leesflow = leesflow/10^6
nyrs = length(leesflow[,1])

### Scale the data..
fmean = apply(leesflow,2,mean)    #monthly mean
fsdev = apply(leesflow,2,sd)	#monthly standard deviation

leesflow.std  = t((t(leesflow) - fmean) / fsdev)

#get variance matrix..
zs.arp=var(leesflow.std)

#do an Eigen decomposition..
zsvd.arp=svd(zs.arp)

#Principal Components...
#pcs.arp=t(t(zsvd.arp$u) %*% t(leesflow.std))
pcs.arp = leesflow.std %*% zsvd.arp$u

#Eigen Values.. - fraction variance 
lambdas.arp=(zsvd.arp$d/sum(zsvd.arp$d))

par(mfrow=c(1,1))
plot(1:7, lambdas.arp[1:7], type="l", xlab="Modes", ylab="Frac. Var. explained")
points(1:7, lambdas.arp[1:7], col="red")

 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)	
arvar=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)
annsim=1:nyrs
annsimpdf=matrix(0,nrow=nsim,ncol=nevals)


## Points to evaluate the May pdf

May = leesflow[,5]
Annual = apply(leesflow,1,sum)
xeval=seq(min(May)-0.25*sd(May),max(May)+0.25*sd(May),length=100)
annxeval=seq(min(Annual)-0.25*sd(Annual),max(Annual)+0.25*sd(Annual),length=100)
nevals=length(xeval)

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

for(k in 1:nsim){
  
  #simulation matrix - 12 rows and N = # of years as columns 
  simdismon=matrix(0, ncol=N,nrow=12)
  
for(i in 1:12){
  if (i <= 3){
    arbest = ar(pcs.arp[,i],order.max=25)
	p = order(arbest$aic)[1]
	
#    p=armabest$order[1]
#    q=armabest$order[3]
#    zz = arima(pcs.arp[,i], order = c(p,0,q), include.mean=TRUE, method="ML")
#    xsim=arima.sim(n=N,list(ar=c(zz$coef[1:p]),ma=c(zz$coef[(p+1):(p+q)])), sd=sqrt(zz$sigma2))+ mean(as.matrix(leesflow.std))

zz = arima(pcs.arp[,i], order = c(p,0,0), include.mean=TRUE, method="ML")
    xsim=arima.sim(n=N,list(ar=c(zz$coef[1:p])), sd=sqrt(zz$sigma2))

	
	}
  else if (i >3){
   xsim = rnorm(N, mean=mean(pcs.arp[,i]), sd=sd(pcs.arp[,i]))
  }
  simdismon[i,]=xsim
}
  
   simdismon1 = t(simdismon)%*%t(zsvd.arp$u)
  simdismon = t(t(simdismon1) *fsdev + fmean)
  
    maysim = simdismon[,5]
  annsim = apply(simdismon,1,sum)
  simpdf[k,]=sm.density(maysim,eval.points=xeval,display="none")$estimate
  annsimpdf[k,]=sm.density(annsim,eval.points=annxeval,display="none")$estimate
  
  #***calculate basic statistics****************
  
  for(j in 1:12){
    
    armean[k,j]=mean(simdismon[,j])
    arvar[k,j]=var(simdismon[,j])
    arskw[k,j]=skew(simdismon[,j])
    armax[k,j]=max(simdismon[,j])
    armin[k,j]=min(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.
#obsstats = stats(leesflow)

### bind the stats of the observed data at the top..
armean=rbind(apply(leesflow,2,mean),armean)
arvar=rbind(apply(leesflow,2,var),arvar)
arskw=rbind(apply(leesflow,2,skew),arskw)  ##skew is the function
arc=1:12
  for(j in 2:12) {
    j1=j-1
    arc[j]=cor(leesflow[,j],leesflow[,j1])
  }
  arc[1]=cor(leesflow[1:nyrs1,12],leesflow[2:nyrs,1])
  
}

arcor=rbind(arc,arcor)
armax=rbind(apply(leesflow,2,max),armax)
armin=rbind(apply(leesflow,2,min),armin)

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)



## Boxplot of statistics
par(mfrow=c(2,3))
boxplot(armean, main="Monthly Mean")
lines(apply(leesflow,2,mean),col="red")
boxplot(arvar, main="Monthly Var")
lines(apply(leesflow,2,var),col="red")
boxplot(arskw, main="Monthly Skew")
lines(apply(leesflow,2,skew),col="red")
boxplot(arcor, main="Monthly Lag-1")
boxplot(armax, main="Monthly Max")
lines(apply(leesflow,2,max),col="red")
boxplot(armin, main="Monthly Min")
lines(apply(leesflow,2,min),col="red")

## boxplot of annual mean, variance, skew and lag-1 correlation
ymean = apply(armean,1,sum)
yvar = apply(arvar,1,sum)
yskw = apply(arskw,1,sum)
ycor = apply(arcor,1,sum)
ymax = apply(armax,1,sum)
ymin = apply(armin,1,sum)

par(mfrow=c(2,3))
boxplot(ymean, main="Ann Mean")
boxplot(yvar, main="Ann Var")
boxplot(yskw, main="Ann Skew")
boxplot(ycor, main="Ann Lag-1")
boxplot(ymax, main="Ann Max")
boxplot(ymin, main="Ann Min")


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


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

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


par(mfrow=c(1,1))
xs = 1:length(xeval)
zz=boxplot(split(t(annsimpdf),xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(annsimpdf,xdensityorig),xlab="Annual 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")