library(copula)
library(ks)
library(kdecopula)



######### this code will simulate all the 12-monthly flows simultaneously
#### using a Copula
### It computes the PDF of May flow and overlays the historical PDF
### compute other stats and do boxpolots, as with problems 1,2


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

## MAF
WYmonflow=WYmonflow/10^6	


N = nrow(WYmonflow)

#################################

### All 12 months together..

u = pobs(WYmonflow)

## or get the CDFs from Kernel density estimators

u = c()
for(i in 1:12){
fhat = kde(WYmonflow[,i])
zz = pkde(WYmonflow[,i],fhat)
u = cbind(u,zz)
}

## Fit Frank Copula

fc = frankCopula(dim=12)
ffc=fitCopula(fc,u)

## Plot the May histogram

hist(WYmonflow[,5],probability=TRUE)
### Simulate
for(isim in 1:250){


zsim = rCopula(N,frankCopula(coef(ffc),dim=12))

## Get the Kernel quantile..
fhat = kde(WYmonflow[,5])
simmay=qkde(zsim[,2],fhat)

sm.density(simmay,add=TRUE,col="grey")

}

sm.density(WYmonflow[,5],add=TRUE,col="red",lwd=2)