# Drawing bivariate random variables
 # Written for Advanced Statistics, Econ. NTU
 # Jin-Lung Lin, April 2007
library(ellipse);
library(fMultivar);
set.seed(12345678);
sig1=1
sig2=2
rho1=0.5
rho2=0
rho3=-0.5
mu=c(3,3)
x=rep(3,10)
y=x


# Graphics for contour of bivariate Normal ( (3,3), sig1,sig2,rho) 
op=par(mfrow=c(3,3))
plot(ellipse(rho1,scale=c(sig1,sig2),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",0.5,"    ",sigma[1]>sigma[2])),side=3)
plot(ellipse(rho1,scale=c(sig1,sig1),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",0.5,"    ",sigma[1]==sigma[2])),side=3)
plot(ellipse(rho1,scale=c(sig2,sig1),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",0.5,"    ",sigma[1]<sigma[2])),side=3)

plot(ellipse(rho2,scale=c(sig1,sig2),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",0,"    ",sigma[1]>sigma[2])),side=3)
plot(ellipse(rho2,scale=c(sig1,sig1),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",0,"    ",sigma[1]==sigma[2])),side=3)
plot(ellipse(rho2,scale=c(sig2,sig1),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",0,"    ",sigma[1]<sigma[2])),side=3)

plot(ellipse(rho3,scale=c(sig1,sig2),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",-0.5,"    ",sigma[1]>sigma[2])),side=3)
plot(ellipse(rho3,scale=c(sig1,sig1),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",-0.5,"    ",sigma[1]==sigma[2])),side=3)
plot(ellipse(rho3,scale=c(sig2,sig1),centre=mu),type="l",col="blue")
abline(h=3)
abline(v=3)
abline(a=0,b=1,col="magenta")
mtext(expression(paste(rho,"=",-0.5,"    ",sigma[1]<sigma[2])),side=3)
par(op)


 # Bivariate Normal plot for varying correlation, 0.5,0,-0.5
   O1=matrix(c(1,0,0,1),2,2)
   O2=matrix(c(1,0,0,4),2,2)
   O3=matrix(c(1,1,1,4),2,2)
   O4=matrix(c(1,-1,-1,4),2,2)
   O5=matrix(c(1,.99,.99,1),2,2)
   O6=matrix(c(1,-.99,-.99,1),2,2)
   n = 51
   x = seq(-3, 3, length = n)
   xoy = cbind(rep(x, n), as.vector(matrix(x, n, n, byrow = TRUE)))
   X = matrix(xoy, n * n, 2, byrow = FALSE)
   head(X)
   # The Bivariate Normal Case:
   Z1 = matrix(dmvsnorm(X, dim = 2,Omega=O1,mu=c(0,0)), length(x))
   Z2 = matrix(dmvsnorm(X, dim = 2,Omega=O2,mu=c(0,0)), length(x))
   Z3 = matrix(dmvsnorm(X, dim = 2,Omega=O3,mu=c(0,0)), length(x))
   Z4 = matrix(dmvsnorm(X, dim = 2,Omega=O4,mu=c(0,0)), length(x))
   Z5 = matrix(dmvsnorm(X, dim = 2,Omega=O5,mu=c(0,0)), length(x))
   Z6 = matrix(dmvsnorm(X, dim = 2,Omega=O6,mu=c(0,0)), length(x))
   
   nobs=1000;
   X1 = matrix(rmvsnorm(nobs, dim = 2,Omega=O1,mu=c(0,0)),nrow=nobs)
   X2 = matrix(rmvsnorm(nobs, dim = 2,Omega=O2,mu=c(0,0)),nrow=nobs)
   X3 = matrix(rmvsnorm(nobs, dim = 2,Omega=O3,mu=c(0,0)),nrow=nobs)
   X4 = matrix(rmvsnorm(nobs, dim = 2,Omega=O4,mu=c(0,0)),nrow=nobs)
   X5 = matrix(rmvsnorm(nobs, dim = 2,Omega=O5,mu=c(0,0)),nrow=nobs)
   X6 = matrix(rmvsnorm(nobs, dim = 2,Omega=O6,mu=c(0,0)),nrow=nobs)
   
   op5=par(mfrow=c(3,2))
   plot(X1,type="p")
   abline(v=0)
   abline(h=0)
   mtext(expression(paste(rho==0,"  ", sigma[1] ==1,"  ", sigma[2]==1),side=3))
   plot(X2,type="p")                                                     
   abline(v=0)                                                           
   abline(h=0)                                                         
   mtext(expression(paste(rho==0,"    ",sigma[1]==1,"  ", sigma[2]==2),side=3)) 
   plot(X3,type="p")                                                        
   abline(v=0)                                                              
   abline(h=0)                                                             
   mtext(expression(paste(rho==0.5,"    ",sigma[1]==1,"  ", sigma[2]==2),side=3))
   plot(X4,type="p")                                                            
   abline(v=0)                                                              
   abline(h=0)                                                             
   mtext(expression(paste(rho==-0.5,"    ",sigma[1]==1,"  ", sigma[2]==2),side=3))
   plot(X5,type="p")                                                            
   abline(v=0)                                                              
   abline(h=0)                                                             
   mtext(expression(paste(rho==0.99,"    ",sigma[1]==1,"  ", sigma[2]==1),side=3))  
   plot(X6,type="p")                                                          
   abline(v=0)                                                              
   abline(h=0)                                                            
   mtext(expression(paste(rho==-0.99,"    ",sigma[1]==1,"  ",sigma[2]==1),side=3))  
   par(op5) 
   
     # Effect of variance   
  op1=par(mfrow=c(1,2)) 
  persp(x, x, Z1, theta = -40, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", 
      main=expression(paste(rho,"=",0,"   ", sigma[1],"=",1,"  ",sigma[2],"=",1)) )
      
  persp(x, x, Z2, theta = -40, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", 
      main=expression(paste(rho,"=",0,"   ", sigma[1],"=",1,"  ",sigma[2],"=",2)) )    
  par(op1)
 
 
  # Effects of correlation       
   op2=par(mfrow=c(2,1))
   persp(x, x, Z3, theta = -40, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Z",
   main=expression(paste(rho,"=",0.5,"   ", sigma[1],"=",1,"  ",sigma[2],"=",1)) ) 
   persp(x, x, Z4, theta = -40, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Z",
   main=expression(paste(rho,"=",-.5,"   ", sigma[1],"=",1,"  ",sigma[2],"=",1)) ) 
   par(op2)
   
   op3=par(mfrow=c(2,1))
   persp(x, x, Z5, theta = -40, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Z",
       main=expression(paste(rho,"=",.99,"   ", sigma[1],"=",1,"  ",sigma[2],"=",1)) ) 
   persp(x, x, Z6, theta = -40, phi = 30, expand = 0.5, col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "X", ylab = "Y", zlab = "Z",
       main=expression(paste(rho,"=",-.99,"   ", sigma[1],"=",1,"  ",sigma[2],"=",1)) )
   par(op3)
   
   
# Graphics exampels from fMultivar
 # Generate Grid Points:
   n = 51
   x = seq(-3, 3, length = n)
   xoy = cbind(rep(x, n), as.vector(matrix(x, n, n, byrow = TRUE)))
   X = matrix(xoy, n * n, 2, byrow = FALSE)
   head(X)
   # The Bivariate Normal Case:
   Z = matrix(dmvsnorm(X, dim = 2), length(x))
   par (mfrow = c(2, 2), cex = 0.7)
   persp(x, x, Z, theta = -40, phi = 30, col = "steelblue") 
   title(main = "Bivariate Normal Plot")
   image(x, x, Z)
   title(main = "Bivariate Normal Contours")
   contour(x, x, Z, add = TRUE)
   # The Bivariate Skew-Student-t Case:
   mu = c(-0.1, 0.1)
   Omega = matrix(c(1, 0.5, 0.5, 1), 2)
   alpha = c(-1, 1)
   Z = matrix(dmvst(X, 2, mu, Omega, alpha, df = 3), length(x))
   persp(x, x, Z, theta = -40, phi = 30, col = "steelblue") 
   title(main = "Bivariate Student-t Plot")
   image(x, x, Z)
   contour(x, x, Z, add = TRUE)
   title(main = "Bivariate Student-t Contours")