### Generate April and May flows jointly
### replace nd with 12 to simulate all
### the 12 monthly flows jointly


library(copula)
library(MASS)
library(ks)
library(sm)


# 1. Generate synthetic data, standard normal marginals
# Read data - Lees Ferry monthly flows

flows=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-3-2018/LeesFerry-monflows-1906-2016.txt"),ncol=13,byrow=T)
flows=flows[,2:13]/10^6

## Select Apr and May flows. May flows has bimodality
data = flows[,4:5]
nd=2

## use empirical CDF
u = pobs(data)


### or fit GAmma PDFs and use the CDFs
## fit Gamma for Apr flows
zgamma = fitdistr(data[,1], "gamma")
uu = pgamma(data[,1], shape=zgamma$estimate[1], scale=1/zgamma$estimate[2])

## Gamma CDF for April and Empirical CDF for May
u = cbind(uu,pobs(data)[,2])


# 2. Fit the copula 
fc=normalCopula(dim=nd,disp='un')
fnc = fitCopula(fc,u)

# 3. random generate from the fitted copula
cop_sim = rCopula(10000,normalCopula(fnc@estimate,dim=nd,dispstr='un'))

# 4.  convert back to data space, assuming standard normal marginals

## X comes from Gamma and Y from the sample
data_sim_x = qgamma(cop_sim[,1],shape=zgamma$estimate[1], scale=1/zgamma$estimate[2])
data_sim_y = quantile(data[,2],cop_sim[,2])

# 5. plot contours of the simulated data and the data together
contour(kde2d(data_sim_x, data_sim_y,n=100))
points(data[,1],data[,2],cex=0.5,pch=19)


contour(kde2d(data_sim_x, data_sim_y,n=100),nlevels=20)

#or specify the contour levels:

contour(kde2d(data_sim_x, data_sim_y,n=100),levels=c(0.3,0.2,0.1,0.05,0.01,0.001))

### Suppose you wish to make many
### simulations, each same length as historical data
## Plot the May histogram

hist(data[,2],probability=TRUE)
N = dim(data)[1]


### Simulate and plot the PDF of May flows

nsim=250
for(isim in 1:nsim){

cop_sim = rCopula(N,normalCopula(fnc@estimate,dim=nd,dispstr='un'))

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

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

}
### add the PDF of the historical data a solid red line
sm.density(data[,2],add=TRUE,col="red",lwd=2)

##### to simulate all the 12 months

data = flows

## or use empirical CDF
u = pobs(data)

## or use Kernel CDFs
nd=12
 u = c()
 for(i in 1:nd){
fhat = kde(data[,i])
 zz = pkde(data[,i],fhat)
 u = cbind(u,zz)
}


# Fit the copula 
fnc = fitCopula(normalCopula(dim=nd,disp='un'),u)
###
## Simulat nsim simulations each of length N
### Compute the monthly stats for each simulation
### boxplot them.
## Use the commands from the code for Problem 1