Multinomial_normal_pd.sas

options ls=78;
title "Bivariate Normal Density";
%let r=0.7;
data a;
  pi=3.1416;
  do x1=-4 to 4 by 0.1;
    do x2=-4 to 4 by 0.1;
      phi=exp(-(x1*x1-2*&r*x1*x2+x2*x2)/2/(1-&r*&r))/2/pi/sqrt(1-&r*&r);
      output;
    end;
  end;
  run;
proc g3d;
  plot x1*x2=phi / rotate=-20;
  run;