days<- read.table("C:/CLIM/Days5.txt",header=T)  
attach(days)
postscript(file="C:/CLIM/DaysCros.ps",height=6,width=16,horiz=F)
#par(mfrow=c(2,2))                            #for further plots


#-----Supplement-----------------------------------------------
ChiSqu<- function(mat,k,m)  #Pearson's chi^2 of a kxm matrix mat 
{n<- sum(mat); ch2<- 0
for(i in 1:k){ni<- sum(mat[i,])                        #row sums
for(j in 1:m){nj<- sum(mat[,j])                     #column sums
eij<- ni*nj/n                              #expected frequencies
ch2<- ch2+(mat[i,j] - eij)^2/eij}}                       
return(ch2)
}
#-----End-of-Supplement-------------------------------------------

Vau<- function(Xa,Xb,m,n,l){   #Cramer's V of an mxm table YaxYb
Ya<- Xa[1:(n-l)]; Yb<- Xb[(1+l):n]                   #time lag l
chi2<-chisq.test(Ya,Yb)$statistic   #Pearson's $chi^2 statistic
##or with the user function ChiSqu:
#tab<- table(Ya,Yb);  chi2<- ChiSqu(tab,m,m)
Vau<- sqrt(chi2/((n-l)*(m-1)))
return(Vau)
}

Vabba<- function(Pra,Prb,m,n,brk,lx){        #cross corr fuction
Preca<- pmin(Pra/10,100); Precb<- pmin(Prb/10,100)    
Prca<-  cut(Preca,brk)                     #Categorical variable
Prcb<-  cut(Precb,brk)                     #Categorical variable
Vab<- 1:(lx+1); Vba<- 1:(lx+1)
for(l in 0:lx){
Vab[l+1]<- Vau(Prca,Prcb,m,n,l)
Vba[l+1]<- Vau(Prcb,Prca,m,n,l)
}
return(cbind(Vab,Vba))        #cbind produces an (lx+1)x2 matrix
}

plotCrosV<- function(mx,cra,crb,yc,xtxt,ytxt,ttxt,Stab,Stba){
plot(0:mx,cra,type="l",lty=1,xlim=c(0,mx),ylim=yc,
            xlab="day lag",ylab="Cramers V")
points(0:mx,cra,pch=16);  title(main=ttxt,cex=0.6)
lines(0:mx,crb,lty=2);  points(0:mx,crb,pch=4)
legend(xtxt,ytxt,legend=c(Stab,Stba),lty=c(1,2)) 
}

#---------------------------------------------------------------
ttxt<- "Daily Precipitation 2004-2010"; ttxt
 
n<- length(Year)   
m<- 6; brk<-c(-1,0,1,2.5,5,10,100)                 #6 categories
"6 Classes: [0],(0,1],(1,2.5],(2.5,5],(5,10],(10,100])"
yc<- c(0.0,0.26); xtxt<- 5; ytxt<- 0.25     #plotting parameters 
lx<- 8                                         #maximal time lag

#---------------------------------------------------------------
Pra<- PrAa; Prb<- PrHo
Stab<- "Aa -> Ho"; Stba<- "Ho -> Aa"
V<- Vabba(Pra,Prb,m,n,brk,lx)
Vab<- V[,1]; Vba<- V[,2]           #1. and 2. column of matrix V
Stab; Vab; Stba; Vba
plotCrosV(lx,Vab,Vba,yc,xtxt,ytxt,ttxt,Stab,Stba)

#Continue with: Pra<- PrAa; Prb<- PrPo  etc.
 
dev.off() 

detach(days)                                            ###
rm(list=objects())                                      ###