options linesize=78 pagesize=999;
libname out '';

***This sas program estimates (using grid search) the parameters of the NFP model 
(see Ert & Erev, 2006) based on the data of the 30 random constant sum games studied 
in Erev, Roth, Slonim & Barron;

************THE 30 GAMES********************;
DATA v; 
input p1-p120;
cards;
                         0.58    0.94    0.98    0.51
                         1.00    0.24    0.05    0.96
                         0.82    0.37    0.19    0.84
                         0.55    0.36    0.26    0.59
                         0.28    0.85    0.78    0.00
                         0.07    0.97    0.87    0.38
                         0.70    0.26    0.41    0.87
                         0.38    0.87    1.00    0.60
                         0.91    0.46    0.49    0.71
                         0.98    0.57    0.15    0.93
                         0.96    0.49    0.05    0.83
                         0.55    0.39    0.14    0.64
                         0.80    0.62    0.67    1.00
                         0.55    0.16    0.08    0.29
                         0.95    0.03    0.28    0.49
                         0.51    0.33    0.05    0.97
                         0.72    0.44    0.60    0.66
                         0.48    0.27    0.02    0.36
                         0.59    0.95    0.78    0.04
                         0.88    0.15    0.56    0.69
                         0.57    0.17    0.22    0.31
                         0.19    0.05    0.01    0.58
                         0.96    0.27    0.78    0.97
                         1.00    0.49    0.71    0.76
                         0.48    0.39    0.00    0.50
                         0.37    0.97    0.83    0.79
                         0.35    0.98    0.81    0.79
                         0.47    0.10    0.46    0.59
                         0.87    0.00    0.03    0.09
                         0.25    0.91    0.94    0.90
;




**************the grid of the different parameters*********************************;

data v; set v; 
model='nfp'; 
 nsim=100; **number of simulations;
 ***the grid: The grid below is large, and running it takes a hour on my PC.  
 To check the program it is better to start with a smaller grid; 
        do lamb=1 to 5 by .25; 
		   do w=.05 to .6 by .05;
   output; 
 end; end; 


data v; set v;
 array paya{30,2,2} p1-p120; **payoff matrix;
 ARRAY sp{30,2} sp1-sp60; **proprtion of A choices in the simulation;
 ARRAY pA{2,2} paa1-paa4; ** propensities;
 ARRAY tpA{2,2} tpa1-tpa4; 
 array av{2} av1-av2; ***a(1);
 array sd{2} sd1-sd2; ***normalization factor;
 array va{2} va1 va2; *obtained payoff;
 array vaf{2} vaf1 vaf2; **forgone payoff;
 array dec{2} dec1 dec2; **decisons;

r=500;
****new game****;

do g=1 to 30;

************************ new virual subject********;

do s=1 to nsim; k=1;

*****initial values ******;
sd1=0; sd2=0; av1=0; av2=0;

do aa=1 to 2; do bb=1 to 2;
 av1=av1+paya{g,aa,bb}/(4);
 av2=av2+(1-paya{g,aa,bb})/(4);
end; end;

do  aa=1 to 2; do bb=1 to 2;
 sd1=sd1+abs((1-av1)*paya{g,aa,bb}+av1*(1-paya{g,aa,bb}))/4;
 sd2=sd2+abs((1-av2)*(1-paya{g,aa,bb})+av2*paya{g,aa,bb})/4;
end; end;

do id=1 to 2;
 do z=1 to 2;
    pa{z,id}=av{id};
 end;

end;


******* the r trials ********************************************************;
do t=1 to r;
do player=1 to 2; a=player;

************ choice *******************************;
***finding the max to reduce the risk of numbers that are too large***;
maxx=.;
do z=1 to 2; if pa{z,a}>maxx then maxx=pa{z,a}; end;
do z=1 to 2; tpa{z,a}=pa{z,a}-maxx; end;
  como=0; fl=0; sa=0; sb=0;

if norma='n' then sd{a}=1;
  DO Z=1 TO 2; Sa=sum(Sa,exp(tPA{Z,a}*lamb/sd{a})); end;

  Ra=Sa*ranuni(0);
  str=.;

  do z=1 to 2;
    COMO=sum(COMO,exp(tPA{Z,a}*lamb/sd{a}));
    IF COMO>=RA AND FL=0 THEN DO; dec{a}=Z; FL=1; END;
  END;
 como=0; fl=0;

end;******of a (player);

**outcomes***;
 Va1=PAYA{g,Dec1,Dec2};
 Va2=1-PAYA{g,Dec1,Dec2};
 Vaf1=PAYA{g,3-Dec1,Dec2};
 Vaf2=1-PAYA{g,Dec1,3-Dec2};

do player =1 to 2; a=player;
****pa statistics********************;
  ca=1; if dec{a}=2 then ca=0;
  sp{g,a}=sum(sp{g,a},ca/(nsim*500));

**** learning ****************************************************************;
    PA{dec{a},a}=PA{dec{a},a}*(1-w) + va{a}*w;
    PA{3-dec{a},a}=PA{3-dec{a},a}*(1-w) + vaf{a}*w;
  
***sd**************************;
 sd{a}=sd{a}*(1-w)+w*abs(va{a}-max(va{a},vaf{a}));

end;****of player;


end; **of trials***;
end; **of sim;
end; *** of g;

zero=0;



****the observed choice proportions;
data obs; zero=0;
input ob1-ob60;
cards;
 0.554 0.474
 0.610 0.460
 0.562 0.596
 0.364 0.510
 0.498 0.494
 0.424 0.526
 0.344 0.548
 0.266 0.330
 0.616 0.362
 0.792 0.220
 0.800 0.426
 0.722 0.650
 0.328 0.650
 0.328 0.240
 0.256 0.224
 0.922 0.578
 0.312 0.168
 0.476 0.604
 0.924 0.612
 0.238 0.526
 0.278 0.174
 0.708 0.684
 0.024 0.860
 0.464 0.336
 0.984 0.014
 0.166 0.226
 0.168 0.314
 0.114 0.372
 0.252 0.018
 0.322 0.558
;

data all; merge v obs; by zero;

 ARRAY ob{30,2} ob1-ob60;
 ARRAY sp{30,2} sp1-sp60; 
 
mse=0; 

do g=1 to 30; do pl=1 to 2; 
  mse=mse+(1/60)*(sp{g,pl}-ob{g,pl})**2;
end; end;

proc print;
var lamb w mse;

run;