#Copyright Marius Hofert

#how to get help?
?mean #get help for the function "mean"

#remove previously used objects and load libs
remove(list=objects()) 
library(copula) #load package "copula" (has to be installed on your system including dependencies)

#initialise MT rng
set.seed(1)

#====user specifications====

#set the path to your working directory
mypath="/Users/mhofert/mhofert/latex/projects/job/2009\ cop/lecture/r"

n=1000 #number of data vectors
d=2 #number of dimensions

#====example 1====

#create a d-dimensional copula object
myindcop=indepCopula(d)

#create n random variates from this copula 
myindcopdata=rcopula(myindcop,n)

#export the random variates to a file
myindcoppath=paste(mypath,"/myindcopdata.dat",sep="") #append a file name to your path
write.table(myindcopdata,myindcoppath,row.names=F,col.names=F) #or use write.matrix() from the library MASS

#read the data in a data frame
mycopdata=read.table(myindcoppath)

#visualize the data in a two-dimensional plot
plot(mycopdata[,1],mycopdata[,2])

#with axis labels
plot(mycopdata[,1],mycopdata[,2],xlab="u",ylab="v")

#with title
plot(mycopdata[,1],mycopdata[,2],xlab="u",ylab="v",main="My first copula data")

#write the plot to a .pdf file
pdf(paste(mypath,"/myfirstcopdata.pdf",sep="")) #open pdf device for plotting
plot(mycopdata[,1],mycopdata[,2],xlab="u",ylab="v",main="My first copula data") #plot
dev.off() #close plotting device

#====example 2====

#define a function
myfunc=function(u,v){
	return(u*v)
}

#define a grid on which the function is evaluated
mylength=10 #grid length
grid=expand.grid(u=seq(0.01,0.99,length=mylength),v=seq(0.01,0.99,length=mylength))

#evaluate the function at the grid points
y=numeric(mylength*mylength) #define a vector of function values
for(i in 1:(mylength*mylength)){
	y[i]=myfunc(grid[[1]][i],grid[[2]][i])
}

#create a data frame of pairs (u,v,y)
data=data.frame(u=grid[[1]],v=grid[[2]],y=y)

#wireframe plot (a similar and easier to use interface is provided by the function persp())
library(lattice) #load library "lattice" (for wireframe plot)
wireframe(data[[3]]~data[[1]]*data[[2]],drape=T,xlab="u",ylab="v",zlab="C(u,v)",main="The independence copula C(u,v)=uv")

#write plot to a .ps file
trellis.device(postscript,color=T,horizontal=F,onefile=F,file=paste(mypath,"/myfirstcopplot.ps",sep=""))
wireframe(data[[3]]~data[[1]]*data[[2]],drape=T,xlab="u",ylab="v",zlab="C(u,v)",main="The independence copula C(u,v)=uv")
dev.off()