data bivnormal;
do i = 1 to 1000;
z1 = rannor(104);
z2 = rannor(104);
z3 = rannor(104);
x = 3*z1+z2;
y = 3*z1+z3;
output;
end;
run;
ods html style=statistical;
ods graphics on;
proc kde data=bivnormal;
bivar x y / plots=contour surface;
run;
ods html close;