########### COPULA #######################################
install.packages('mvtnorm')
require(mvtnorm)
S <- matrix(c(1,.8,.8,1),2,2) #Correlation matrix
AB <- rmvnorm(mean=c(0,0),sig=S,n=1000) #Our gaussian variables
U <- pnorm(AB) #Now U is uniform - check using hist(U[,1]) or hist(U[,2])
x <- qgamma(U[,1],2) #x is gamma distributed
y <- qbeta(U[,2],1,2) #y is beta distributed
plot(x,y) #They correlate!

########## tCopula #######################################
# Here we will show two ways to simulate a bivariate dataset with
# Normal(0,1) marginals but different dependency structures.

# (a) (Bivariate normal distribution)

 rho=0.8
 z1=rnorm(10000)
 z2=rho*z1+sqrt(1-rho^2)*rnorm(10000)
	# the last two lines mean that both z1 and z2 are N(0,1)
	# (z1,z2) is bivariate normal with cor(z1,z2)=rho
 plot(z1,z2,pch=20,cex=0.5)	# pch=20 means that we get solid dots
				# cex=0.5 means character expansion 0.5 times normal

# next do a QQnorm plot to verify graphically that the marginals are normal

 qqnorm(z1)
 qqnorm(z2)

# (b) Normal marginals using the t-copula
#     Keep z1 and z2 as above, but now use these in a different way
#     in order to create a standard bivariate t

 nu=3
 y=rchisq(10000,df=nu)	# generates 10000 i.i.d. chi-sqaured r.v.'s
			# with nu degrees of freedom
 x1=z1/sqrt(y/nu)	# this defines a standard t_nu random variable 
 x2=z2/sqrt(y/nu)
 u1=pt(x1,df=nu)	# u1 will be i.i.d. ~ U[0,1]
 u2=pt(x2,df=nu)

# first plot the empirical t-copula 

 plot(u1,u2,pch=20,cex=0.5)

# now generate the normal random variables

 w1=qnorm(u1)
 w2=qnorm(u2)
 plot(w1,w2,pch=20,cex=0.5)

# It helps to be able to compare the results of (a) and (b) side
# by side as follows:

 par(mfrow=c(1,2))	# mfrow=c(1,2) means set up a plotting region that
			# will have one row and two plots per row
 plot(w1,w2,pch=20,cex=0.5)
 plot(z1,z2,pch=20,cex=0.5)
 
###### Empirical Distribution Function ###################
# set the seed for consistent replication of random numbers
set.seed(1)

# generate 150 random numbers from the standard normal distribution
normal.numbers = rnorm(150)

# empirical normal CDF of the 150 normal random numbers
normal.ecdf = ecdf(normal.numbers)

# plot normal.ecdf (notice that the only argument needed is normal.ecdf)
# use png() and dev.off() to print this plot to your chosen folder
# png('INSERT YOUR DIRECTORY PATH HERE/ecdf standard normal.png')

plot(normal.ecdf, xlab = 'Quantiles of Random Standard Normal Numbers', ylab = '', main = 'Empirical Cumluative Distribution\nStandard Normal Quantiles')

# add label to y-axis with mtext()
# side = 2 denotes the left veritical axis
# line = 1.8 sets the position of the label
mtext(text = expression(hat(F)[n](x)), side = 2, line = 1.8)