/* revcons1.gau  for Revenue-Constrained STIP Note 21-22/5/2003:
  Simulate welfare and profits as functions of R&D and output subsidies in R&D+Cournot game and print Fig. 1*/
new;   /*************************************************************************/
/*      USER-SPECIFIED VARIABLES                                     */
eta = 0.2;   /* Relative Effectiveness of R&D */
stepsize=.005;   /* stepsize in s1 (R&D Subsidy) loops: .005 for publication; insufficient workspace memory for 0.0025 */
stepsizeB=stepsize;    /* stepsize in s2 (Output Subsidy) loops  */
s1min=-0.5;
s1max=1.0;
n=1+(s1max-s1min)/stepsize;  /* n: # of values of s1 considered (=31 for stepsize=.05)*/
s2min=s1min;
m=n;                             /* m: # of values of s2 considered */
/*************************************************************************/
/*              INITIALISE VECTORS */
s1=zeros(1,n);  /* A row, not a column, vector: a mistake, but too tedious to rectify it 22.5.03 */
s2=zeros(m,1);
W=zeros(m,n);
Ppi=zeros(m,n);
s2PCW=zeros(1,n);
/*************************************************************************/
/*  Preliminary calculations (independent of s1 and s2)  */
Delt=(1-(4/3)*eta)*(3-(4/3)*eta);  /* Det. of LHS Matrix of Generalised Reac. Funcs */
CC=0.5/(1-(2/3)*eta) -(4/9)*eta;
DD=(1-(4/3)*eta)/Delt;
EE=2*(1-(2/3)*eta)/Delt;
WFT=(1-(8/9)*eta)*DD^2;  /* Welfare in free trade; q=DD, x=(4/3)*eta*q */
sig0=(4/3)*DD;
sig1=2*(1+2*EE*eta^2);
sig2=(7/3)*EE;
phi0=DD;phi1=(7/3)*eta*EE;
phi2=2*EE;
ome0=(1/3-CC*EE)*DD;
ome1=1+(5/3-CC*EE)*eta*EE;
ome2=(4/3-CC*EE)*EE; 
psi0=-CC*DD*EE;
psi1=(4/3-CC*EE)*eta*EE;
psi2=(1-CC*EE)*EE;
AAA=ome2*phi2-psi2*sig2;
/*************************************************************************/
i=1;                                  /*      BEGIN i - s1 (R&D Subsidy) DO LOOP    */
ss1=s1min-stepsize;  /* ss1: Current value of s1 */
DO until i>n;
ss1=ss1+stepsize;
s1[1,i]=ss1;		
/*************************************************************************/
/*  Calculate profit-constrained welfare (PCW) locus:  W_1.S_2=W_2.S_1 */
sig3=sig0+sig1*ss1;
phi3=phi0+phi1*ss1;
ome3=ome0+ome1*ss1;
psi3=psi0+psi1*ss1;
BBB=(ome3*phi2+ome2*phi3)-(psi3*sig2+psi2*sig3);
CCC=ome3*phi3-psi3*sig3;
DDD=BBB^2-4*AAA*CCC;  
IF DDD<=0;
	s2PCW[1,i]=-0.6;  /* WAS 1 */
	GOTO ENDPCW;
ENDIF;  
s2PCW[1,i]=(-BBB+DDD^0.5)/(2*AAA); /* Calculate the Profit-Constrained Welfare Locus */
ENDPCW:
/*************************************************************************/
/*              BEGIN j - s2 (Output Subsidy) DO LOOP    */
j=1;
ss2=s2min-stepsizeB;  /* ss2: Current value of s2 */
DO until j>m;
ss2=ss2+stepsizeB;
s2[j,1]=ss2;		
/*************************************************************************/
q= (1-(4/3)*eta +(2-(4/3)*eta)*(eta*ss1+ss2)) / Delt;  /* Home Output */
qstar= (1-(4/3)*eta - (eta*ss1+ss2)) / Delt;           /* Foreign Output */
x=eta*((4/3)*q+ss1);   /* Home R&D */
w[j,i]=((q-ss2)*q - 0.5*eta*((4/3)*q+ss1)^2)/WFT; /* Welfare at delta=1, normalised by FT value */
ppi[j,i]=w[j,i]+(ss1*x+ss2*q)/WFT;
/***  Check for negative output, profit or R&D levels and set contours accordingly ***/
  IF qstar < 0;
        ppi[j,i]=8;   /* set at this level to avoid unnecessary contour lines */
      GOTO W_EQ_0;
  ENDIF;                                   
  IF q <0;				
  ppi[j,i]=0;   W[j,i]=0;
  ENDIF;					
   IF x <0;
  ppi[j,i]=0;   W[j,i]=0;
  ENDIF;								
IF ppi[j,i] <0;   GOTO W_EQ_0;   ENDIF;  
  GOTO ENDLOOP;
  W_EQ_0:
  W[j,i]=0;
/*************************************************************************/
ENDLOOP:
        j=j+1;		/*  END OF j LOOP */
    endo;
    i=i+1;			/*  END OF i LOOP */
endo;
/*************************************************************************/
library pgraph;		/* COMPUTE AND PRINT FIGURES */
BEGWIND;
MAKEWIND(7.69,6.27,.6,0,0);  /*Create a transparent window that fills the page */
graphset;

/********************     PRINT FIG. 1     **************************/
graphset;
title("Fig. 1:  Revenue-Constrained Optimal Subsidy Locus in R&D-Cournot Game");     /* Home Profits (Contour) */
fonts("complex simgrma");  /* load fonts: \201: standard; \202: Greek/maths */
xlabel("s]1[");
ylabel("s]2[");
_paxes=1;
_ptitlht=0;
_pcross=1;      /* axes intersect at zero */
xtics(-.5,0.5,0.5,1); 
ytics(-.5,1,0.5,1); 
ztics(0,1,1,0);  
xy ( s1', 1-(4/3)*eta - eta*s1');  /* Exact s2(s1) boundary corresponding to q*=0 */
s2_x0 = - (  (3/4)*delt*s1'  +  1  -  (4/3)*eta  )   /   (2-(4/3)*eta)   -   eta*s1' ;  /* Exact s2(s1) boundary corresponding to x (R&D)=0; only for s1<0 */
s2_q0 =  (-(1-(4/3)*eta)/(2-(4/3)*eta) - eta*s1') ;  /* Exact s2(s1) boundary corresponding to q=0; only for s1>0 */
xy ( s1' ,  maxc (  (s2_x0 ~  s2_q0)' )  ) ;  /* Display the maximum of the two preceding boundaries */
	_pmsgctl= { -0.25 0.38  .15  0  1  2  0};  /* message control */
	_pmsgstr = "A" ;                   /* legend string */
ztics(1,8,1,0);  	contour(s1,s2,Ppi); 
ztics(0,.5,0.5,0);  	contour(s1,s2,Ppi);   /* Add a contour line at pi=0.5 */
ztics(1,1.20,0.1,0);  	contour(s1,s2,W); 
ztics(0,0.8,0.2,0);  contour(s1,s2,W); 
xy(s1',s2PCW');   /* Draw Profit-Constrained Welfare Locus  */
PRINTEND:
ENDWIND;
finish:
end