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

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,sum)
sesrain=sesrain[N2:N3]
X=sesrain/10	#cm

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]

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

## Fit Frank Copula and plot PDF

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

params = coef(ffc)

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

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


### Plot CDF


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



## simulate from fitted Copula
par(mfrow=c(2,2))
plot(u[,1],u[,2])

zsim = rCopula(500,frankCopula(coef(ffc),dim=nd))
fhat = kde(Z[,1])
simX=qkde(zsim[,1],fhat)

fhat = kde(Z[,2])
simY=qkde(zsim[,2],fhat)

plot(Z[,1],Z[,2])
points(simX,simY,col='red')

hist(Z[,1],freq=FALSE)
sm.density(simX,add=TRUE,col="blue")
sm.density(Z[,1],add=TRUE,col="red")
points(Z[,1],rep(0,nrow(Z)))

hist(Z[,2],freq=FALSE)
sm.density(simY,add=TRUE,col="blue")
sm.density(Z[,2],add=TRUE,col="red")
points(Z[,2],rep(0,nrow(Z)))