/* DEQUIL.G: procedure to compute equilibria of the prisoner's dilemman game

             and the derivative of the equilibria with respect to the
	     (identified) parameters of the players' payoff functions.

	     This is basically the same as equil.g except for the 
	     extra 2 returns, the derivative of the equilibrium
	     probability of confessing by player a with respect to the
	     parameters (t1_a,t2_a,t1_b,t2_b) where t1_a is the
	     coefficient of x_a and t2_a is the coefficient of x_a*p_b
	     in player a's best response function, and t2_a is the
	     coefficient of x_b and t2_b is the coefficient of x_b*p_a
	     in player 2's best response function.

	     Fixed point (equilibrium condition) is given by

	     p_a(x_a,x_b)=br2_a(x_a,x_b,p_a(x_a,x_b))  where

             br2_a(x_a,x_b,p_a)=br_a(br_b(p_a,x_b),x_a)

	     and 
	     
	     br_a(p_b,x_a)=1/(1+exp(t1_a*x_a+t2_a*x_a*p_b))
	     br_b(p_a,x_b)=1/(1+exp(t1_b*x_b+t2_b*x_b*p_a))

	     By the Implicit function theorem we have

	     d(p_a)/d(t1_a) = [dbr2_a/dt1_a]/[1-dbr2_a(x_a,x_b,p_a)/dp_a]
	     d(p_a)/d(t2_a) = [dbr2_a/dt2_a]/[1-dbr2_a(x_a,x_b,p_a)/dp_a]
	     d(p_a)/d(t1_b) = [dbr2_a/dt1_b]/[1-dbr2_a(x_a,x_b,p_a)/dp_a]
	     d(p_a)/d(t2_b) = [dbr2_a/dt2_b]/[1-dbr2_a(x_a,x_b,p_a)/dp_a]

	     where

	     dbr2_a/dt1_a = -p_a(1-p_a)*x_a
	     dbr2_a/dt2_a = -p_a(1-p_a)*x_a*p_b
	     dbr2_a/dt1_b = p_a(1-p_a)*t2_a*x_a*p_b(1-p_b)*x_b
	     dbr2_a/dt2_b = p_a(1-p_a)*t2_a*x_a*p_b(1-p_b)*x_b*p_a

	     Under the symmetry condition that 

	     t1_a = t1_b
	     t2_a = t2_b

	     then there are effectively only 2 parameters (t1,t2) and the
	     derivatives are

	     d(p_a)/d(t1) = d(p_a)/d(t1_a) + d(p_a)/d(t1_b)
	     d(p_a)/d(t2) = d(p_a)/d(t2_a) + d(p_a)/d(t2_b)

	     The derivatives of p_b with respect to (t1_a,t2_a,t1_b,t2_b) 
	     (using the chain rule and the condition p_b=br_b(p_a,x_b), where
	     p_a is the equilibrium probability a confesses) is:

	     d(p_b)/d(t1_a) = [d(br_b)/d(p_a)] d(p_a)/d(t1_a)
	     d(p_b)/d(t2_a) = [d(br_b)/d(p_a)] d(p_a)/d(t2_a)
	     d(p_b)/d(t1_b) = [d(br_b)/d(p_a)] d(p_a)/d(t1_a) + d(br_b)/d(t1_b)
	     d(p_b)/d(t2_b) = [d(br_b)/d(p_a)] d(p_a)/d(t2_a) + d(br_b)/d(t2_b)

	     The derivatives are returned as two 4x1 vectors,

	     dpa=(dp_a/dt1_a,dp_a/dt2_a,dp_a/dt1_b,dp_a/dt2_b)'
	     dpb=(dp_b/dt1_a,dp_b/dt2_a,dp_b/dt1_b,dp_b/dt2_b)'

	     John Rust, University of Maryland, November 2005 */

    proc (4)=dequil(xa,xb);

       local pt,pt1,df,i;
       local dpa,dpb,dbr2dt1a,dbr2dt2a,dbr2dt1b,dbr2dt2b;
       local dbradpb,dbr2dpa,dbrbdpa,dbrbdt1b,dbrbdt2b;

       pt1=1;
       pt=0;

       i=1;

       if (sa_steps > 0);

          if debg_e;

	     "Equil.g: Initial successive approximation steps ";

          endif;

          do until (abs(pt-pt1) < sa_tol);
             
	     pt=pt1;
	     {pt1,df}=br2_a(xa,xb,pt1); 

	     if debg_e;
                "successive approximation "$+ftos(i,"*.*lf",3,0);; 
		pt;; pt1;; abs(pt1-pt);
	     endif;

	     i=i+1;

             if (i > sa_steps); 
               "warning: terminating successive approximation iterations at ";;
	       ""$+ftos(i,"*.*lf",3,0)$+" iterations.";
               break; 
             endif;

          endo;
	  
       endif;

       if debg_e;
         "starting equil with p=";; pt1;
       endif;

       i=0;

       do until (abs(pt-pt1) < ctol);

       i=i+1;

       if (i > maxnewts); 
         "warning: terminating Newton iterations at ";;
	 ""$+ftos(maxnewts,"*.*lf",3,0)$+" iterations.";
	 "xa="$+ftos(xa,"*.*lf",6,3)$+" xb="$+ftos(xb,"*.*lf",6,3);
         break; 
       endif;

       pt=pt1;

       {pt1,df}=br2_a(xa,xb,pt);


       if debg_e;
         "newton "$+ftos(i,"*.*lf",1,0);; 
	 pt;; (pt-pt1);; (1-df);; pt-(pt-pt1)/(1-df); 
       endif;

       pt1=pt-(pt-pt1)/(1-df);

       if (pt1 < 0);
          pt1=pt/2;
       endif;

       endo;

       dbr2dpa=(1-df);
       dbradpb=df;

       {pt,df}=br_b(pt1,xb);

       dpa=zeros(4,1);
       dpb=dpa;

       if (df > 0);

         dbradpb=dbradpb/df;
     
       else;
        
	 dbradpb=0;

       endif;

       dbr2dt1a=-pt1*(1-pt1)*xa;
       dbr2dt2a=dbr2dt1a*pt;
       /* dbr2dt1b=pt1*(1-pt1)*(dc_a-cc_a+cd_a-dd_a)*pt*(1-pt)*xb;
       dbr2dt2b=dbr2dt1b*pt1; */
       dbr2dt1b=-pt*(1-pt)*xb*dbradpb;
       dbr2dt2b=dbr2dt1b*pt1;

       dpa[1]=dbr2dt1a/dbr2dpa;
       dpa[2]=dbr2dt2a/dbr2dpa;
       dpa[3]=dbr2dt1b/dbr2dpa;
       dpa[4]=dbr2dt2b/dbr2dpa;

       dbrbdt1b=-pt*(1-pt)*xb;
       dbrbdt2b=-pt*(1-pt)*xb*pt1;

       dpb[1]=df*dpa[1];
       dpb[2]=df*dpa[2];
       dpb[3]=df*dpa[3]+dbrbdt1b;
       dpb[4]=df*dpa[4]+dbrbdt2b;

       retp(pt1,pt,dpa,dpb);

    endp;