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-p40;
cards;
                                  0.77    0.35    0.08    0.48
                                  0.73    0.74    0.87    0.20
                                  0.63    0.08    0.01    0.17
                                  0.55    0.75    0.73    0.60
                                  0.05    0.64    0.93    0.40
                                  0.46    0.54    0.61    0.23
                                  0.89    0.53    0.82    0.92
                                  0.88    0.38    0.40    0.55
                                  0.40    0.76    0.91    0.23
                                  0.69    0.05    0.13    0.33
;




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

data v; set v; 
model='nfp'; 
 nsim=200; **number of simulations;
        do lamb=2.5;
		   do w=.25;
   output; 
 end; end; 


data v; set v;
 array paya{10,2,2} p1-p40; **payoff matrix;
 ARRAY sp{10,2} sp1-sp20; **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 10;

************************ 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;


  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-ob180;
cards;
          0.892    0.634    0.566    0.114    0.512    0.832    0.506    0.590    0.674
          0.328    0.490    0.588    0.106    0.196    0.174    0.262    0.436    0.278
          0.770    1.000    0.582    0.996    0.870    0.894    0.802    0.976    0.674
          0.458    0.090    0.452    0.524    0.360    0.374    0.218    0.404    0.368
          0.622    0.944    0.572    0.668    0.322    0.324    0.832    0.540    0.422
          0.284    0.510    0.000    0.220    0.096    0.108    0.500    0.096    0.184
          0.094    0.328    0.230    0.356    0.282    0.470    0.004    0.526    0.178
          0.742    0.404    0.566    0.664    0.662    0.462    0.316    0.456    0.242
          0.426    0.356    0.654    0.348    0.262    0.346    0.388    0.356    0.270
          0.508    0.224    0.098    0.484    0.294    0.194    0.380    0.368    0.334
          0.746    0.476    0.602    0.598    0.616    0.608    0.638    0.692    0.762
          0.620    0.546    0.420    0.436    0.268    0.314    0.464    0.420    0.200
          0.000    0.310    0.696    0.118    0.284    0.478    0.296    0.290    0.184
          0.498    0.612    0.252    0.560    0.424    0.454    0.516    0.548    0.830
          0.278    0.544    0.612    0.000    0.506    0.428    0.512    0.458    0.260
          0.026    0.140    0.188    0.408    0.192    0.260    0.264    0.282    0.276
          0.598    0.388    0.646    0.554    0.762    0.630    0.400    0.546    0.530
          0.424    0.468    0.530    0.446    0.726    0.312    0.482    0.170    0.484
          0.278    0.098    0.710    0.420    0.236    0.146    0.362    0.228    0.404
          0.224    0.156    0.374    0.344    0.160    0.000    0.160    0.000    0.402
;

data all; merge v obs; by zero;

 ARRAY ob{10,2,9} ob1-ob180;
 ARRAY sp{10,2} sp1-sp20;
 ARRAY avg{10,2} avg1-avg20; 
 
mse=0; 

do g=1 to 10; do pl=1 to 2; avg{g,pl}=0;
  do id=1 to 9; 
    avg{g,pl}=avg{g,pl}+ob{g,pl,id}/9; **mean over subjects;
    mse=mse+(1/180)*(sp{g,pl}-ob{g,pl,id})**2; ***mse;
  end; 
end; end;
s2=0;
do g=1 to 10; do pl=1 to 2; do id=1 to 9;
   s2=s2+((ob{g,pl,id}-avg{g,pl})**2)*(1/8)*(1/20); ****pooled variance;
end; end; end;
eno=s2/(mse-s2); ***eno;


proc print;
var lamb w mse s2 eno;

run;