/* EQUILPLOT.GPR: program to plot equilibriun of prisoner's dilemma game

		  John Rust, University of Maryland, November 2005 */
   library pgraph;
   pqgwin many;

   #include setup.gpr;

   nplot=50;

   x=seqa(0,1/nplot,nplot+1);
   iloc=0;
   q_a=0;
   q_b=0;


   nobs=500;

   x_a=rndu(nobs,1);
   x_b=rndu(nobs,1);

   p_a=zeros(nobs,1);
   p_b=zeros(nobs,1);

   y=zeros(nobs,2);
   seeds=rndu(nobs,2);

   qa=q;
   qb=q;

   i=1; do until i > nobs;

      {p_a[i],p_b[i]}=equil(x_a[i],x_b[i]);

      y[i,1]=(seeds[i,1] < p_a[i]);
      y[i,2]=(seeds[i,2] < p_b[i]);

   i=i+1; endo;


  /* first plot exact, or "true" equilibrium probabilities 
     for a grid of possible (x_a,x_b) values */

  z=zeros(rows(x),rows(x));
  z1=zeros(rows(x),rows(x));
  i=1; do until i > rows(x);
    j=1; do until j > rows(x);

       {pa,pb}=equil(x[i],x[j]);

       z[i,j]=pa;
       z1[i,j]=pb;

    j=j+1; endo;
  i=i+1; endo;


   graphset;
   _pdate="";
   _pgrid=0|0;
   _pltype=6|6;
   _plwidth=5|5;
   xtics(0,1,.1,0);
   ytics(0,1,.1,0);
   ztics(.5,1,.1,0); 
   title("Equilibrium Probability A confesses\LPrisoner's Dilemma Game");
   xlabel("Prisoner a's x]a[ value");
   ylabel("Prisoner b's x]b[ value");

#IFUNIX
  let cv = 100 100 640 480 0 0 1 6 15 0 0 2 2;
  wxy = WinOpenPQG(cv,"Expected Payoffs in Prisoner's Dilemma Game","");
  call WinSetActive(wxy);
#ENDIF

  surface(x',x,z');

#IFUNIX
  call WinSetActive(1);
#ENDIF

   graphset;
   _pdate="";
   _pgrid=0|0;
   _pltype=6|6;
   _plwidth=5|5;
   xtics(0,1,.1,0);
   ytics(0,1,.1,0);
   ztics(.5,1,.1,0); 
   title("Equilibrium Probability B confesses\LPrisoner's Dilemma Game");
   xlabel("Prisoner a's x]a[ value");
   ylabel("Prisoner b's x]b[ value");

#IFUNIX
  let cv = 100 100 640 480 0 0 1 6 15 0 0 2 2;
  wxy = WinOpenPQG(cv,"Expected Payoffs in Prisoner's Dilemma Game","");
  call WinSetActive(wxy);
#ENDIF

  surface(x',x,z1');

#IFUNIX
  call WinSetActive(1);
#ENDIF

  /* now estimate these surfaces nonparametrically and plot them 
     for comparison */


  if (npemeth $== "series");

  ep_a=series(0~0,y[.,1],x_a~x_b,"q_a");
  ep_b=series(0~0,y[.,2],x_a~x_b,"q_b");

  load q_a;
  load q_b;
  load iloc;

  endif;

  z=zeros(rows(x),rows(x));
  z1=zeros(rows(x),rows(x));
  i=1; do until i > rows(x);
    j=1; do until j > rows(x);

      if (npemeth $== "kernel");
      
       ep_a=mvks(x[i]~x[j],y[.,1],x_a~x_b);
       ep_b=mvks(x[i]~x[j],y[.,2],x_a~x_b); 

      elseif (npemeth $== "llr");

       ep_a=llr(x[i]~x[j],y[.,1],x_a~x_b);
       ep_b=llr(x[i]~x[j],y[.,2],x_a~x_b); 

      else;

       dat=seriesdt(0,x[i],x[j]);
       ep_a=1/(1+exp(dat[iloc]*q_a));
       ep_b=1/(1+exp(dat[iloc]*q_b));

      endif;

       z[i,j]=ep_a;
       z1[i,j]=ep_b;

    j=j+1; endo;
  i=i+1; endo;


   graphset;
   _pdate="";
   _pgrid=0|0;
   _pltype=6|6;
   _plwidth=5|5;
   xtics(0,1,.1,0);
   ytics(0,1,.1,0);
   ztics(.5,1,.1,0); 
   title("Estimated Probability A confesses\LPrisoner's Dilemma Game");
   xlabel("Prisoner a's x]a[ value");
   ylabel("Prisoner b's x]b[ value");

#IFUNIX
  let cv = 100 100 640 480 0 0 1 6 15 0 0 2 2;
  wxy = WinOpenPQG(cv,"Expected Payoffs in Prisoner's Dilemma Game","");
  call WinSetActive(wxy);
#ENDIF

  surface(x',x,z');

#IFUNIX
  call WinSetActive(1);
#ENDIF

   graphset;
   _pdate="";
   _pgrid=0|0;
   _pltype=6|6;
   _plwidth=5|5;
   xtics(0,1,.1,0);
   ytics(0,1,.1,0);
   ztics(.5,1,.1,0); 
   title("Estimated Probability B confesses\LPrisoner's Dilemma Game");
   xlabel("Prisoner a's x]a[ value");
   ylabel("Prisoner b's x]b[ value");

#IFUNIX
  let cv = 100 100 640 480 0 0 1 6 15 0 0 2 2;
  wxy = WinOpenPQG(cv,"Expected Payoffs in Prisoner's Dilemma Game","");
  call WinSetActive(wxy);
#ENDIF

  surface(x',x,z1');

#IFUNIX
  call WinSetActive(1);
#ENDIF