library(copula)
library(ks)
library(kdecopula)
library(sm)
library(VineCopula)


### Conditional PDF of rainfall at ENSO index = 0.5C f(X | Y = Yc = 0.5C)

rain=matrix(scan("http://civil.colorado.edu/~balajir/r-session-files/aismr.txt"),ncol=13,byrow=T)

#rain=matrix(scan("aismr.txt"),ncol=13,byrow=T)
N=1999-1900+1

N2=1900-1871+1
N3=N2+N-1

sesrain=apply(rain[,7:10],1,mean)
sesrain=sesrain[N2:N3]

n=length(sesrain)

#read the NINO3 index..

test=matrix(scan("http://civil.colorado.edu/~balajir/r-session-files/nino3-index.txt"), ncol=3, byrow=T)


#test=matrix(scan("nino3-index.txt"),ncol=3,byrow=T)
nino3=matrix(test[,3],ncol=12,byrow=T)

#get the nino3 of the summer (June - Sep) season
nino3ses=apply(nino3[,6:9],1,mean)

nino3ses=nino3ses[1:n]


###

X = sesrain
Y = nino3ses

Z = cbind(X,Y)

u = pobs(Z)
nd = ncol(Z)
N = nrow(Z)

## or get the CDFs from Kernel density estimators

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

### Condition at NINO3 index = 0.5C
### Yc = 0.5C
### Compute the CDF corresponding to this.

Yc = 0.5  ## condition at NINO3 index = 0.5

fhat = kde(Z[,2])
zz = pkde(Yc,fhat)

nsim=5000	## number of simulations
Ycc = rep(zz,nsim)


## Fit Frank Copula or Normal Copula and plot PDF

fc = frankCopula(dim=nd)

#fc = normalCopula(dim=nd)

ffc=fitCopula(fc,u)

params = coef(ffc)

mycopula = frankCopula(dim=2,param=params)

#mycopula = normalCopula(dim=2,param=params)

contour(mycopula, dCopula, xlab="Rain",ylab="NINO3")
points(u[,1],u[,2])

## simulate from Bivariate Copula - Standard
zz=BiCop(family=5, par=params)		#Frank Copula
#zz=BiCop(family=1, par=params)		#Normal Copula

## simulate from conditional Copula

ycsim = BiCopCondSim(nsim, cond.val=Ycc, cond.var=2,zz)
### Plot CDF

fhat = kde(Z[,1])
simY=qkde(ycsim,fhat)

hist(simY,probability=TRUE,xlab="Rain (mm) | NINO3 = 0.5")
plot(kde(simY), add=TRUE)

## or
#zz=sm.density(simY, add=TRUE)