%MACRO bivnorm(rho=);
/*************************************************************************
Macro File: bivnorm.sas created on 2/6/07
Macro variable: rho(correlation)
Macro out: 	3D surface and contour plots of the standard bivariate normal 
			density with a given correlation
**************************************************************************/

PROC IML;
*Initialization;
 
i=1; pi = 22/7;
mat1 = shape(0, 6561, 3); * 81*81=6561;

*Do roop to evaluate density for square grid between -4 and 4 (increment=0.1);
  do y1 = -4 to 4 by 0.1;
    do y2 = -4 to 4 by 0.1;

      fy = (1/(2*pi*sqrt(1-&rho**2))) * exp((-1/(2*(1-&rho**2))) * (y1**2 + y2**2 - 2*&rho*y1*y2)); 

	  mat1[i,] = y1 || y2 || fy;  
	  i=i+1;
    end;
  end;

  col={ "y1" "y2" "fy"};
  create data1 from mat1 [colname=col];
  append from mat1;

quit;

*Create 3D surface and contour plot;
PROC G3D data=data1;
title1 '3D bivariate normal density plot';
plot y2*y1=fy / grid rotate=-20 tilt=40;
FOOTNOTE J=L H=1 "rho=&rho";
run;

PROC GCONTOUR data=data1;
title2 'contour plot';
	plot y2*y1=fy / grid;
run;

%MEND bivnorm; 
***************end of MACRO*******************************;

DM "output; clear; log; clear";

%bivnorm(rho=0.7); 
run;  *run the macro with your choice of rho;