Note: 1) This is a simulation(Monte Carlo), the results are close to the theoretical results.
          2) The theoretical results in the paper is valid for the assumption that offerings is the same within a given chioce -- which I implemented as,
    x1=    0.85+0.15*1;
    x2=    0.85+0.15*2;
    x3=    0.85+0.15*3;
    x4=    0.85+0.15*4;

         3) The numeric definition should be valid at any time.
         4) adjusted chg=0.01 to have a 'better' numerica derivative.
         5) when calculate elasticity for x1, everything else needs to be fixed including the errors.

Hope this is clear.





options notes;
%let size=80000;   ***simulation size;
%let Eulerconstant=0.5772156649015328606;


***error of extreme distribution;
data a;
  seed=7813430;
  
  do i=1 to &size;
    e1=-log(-log(ranuni(seed)) ) - &Eulerconstant;
    e2=-log(-log(ranuni(seed)) ) - &Eulerconstant;  
    e3=-log(-log(ranuni(seed)) ) - &Eulerconstant;  
    e4=-log(-log(ranuni(seed)) ) - &Eulerconstant;  
    ****alternative offering values;
    x1=    0.85+0.15*1;
    x2=    0.85+0.15*2;
    x3=    0.85+0.15*3;
    x4=    0.85+0.15*4;

    output;
  end;
  keep e: x:;
run;

%macro e(ovar=, factor=, dsn=);
data tmp;
   set a;
   ***4 choices**;

   array x(*)  x1-x4 ;
   array u(*) u1-u4;

   array e(*) e1-e4;
   do i=1 to 4;
       if upcase(vname(x(i)))= %upcase("&ovar") then  x(i)=&factor+x(i);
   end;
   ***unitility function, the parameters estimated from the model***;
  *** may have many covariates****;
  ***for simplicity, coef=1 and only one variable;
   do i=1 to 4;
       u(i)=0 + 1*x(i)  + e(i) ;
   end;

   max=max(u1,u2,u3,u4);

   do i = 1 to 4;
      if u(i)=max then choice=i;     
   end;
run;

****calculate means on each choice****;
proc sql;
  create table &dsn as
  select choice, count( choice ) as n, calculated n/&size as p,
           mean(x1) as x1m , mean(x2) as x2m,
           mean(x3) as x3m , mean(x4) as x4m
from tmp
group by  choice
;
quit;

proc print data=&dsn;run;
%mend;

%macro e2(dsn1=x1,dsn2=xx1, ovar=x1,chg=0.01) ;
%e(ovar=&ovar, factor=0, dsn=&dsn1);
%e(ovar=&ovar, factor=&chg, dsn=&dsn2);

data all;
  set &dsn1(keep=n choice p &ovar.m  rename=(p=p1  &ovar.m=&ovar));
  set &dsn2(keep=n p &ovar.m rename=(p=p2 n=n2 &ovar.m=&ovar._2));
  e1= ( (p1-p2)/(&ovar-&ovar._2) ) *( &ovar/p1);    ***elasticity evaluate at first point x(i);
  e2= ( (p1-p2)/(&ovar-&ovar._2) ) *(&ovar._2/p2);  ***elasticity evaluate at second point define as x(i)=&factor+x(i);
  &ovar.Chg=put( 100*(&ovar-&ovar._2) /&ovar, 8.2); ***change in x(i) in %;
run;

proc print; run;
%mend;

%e2(dsn1=x1,dsn2=xx1, ovar=x1,chg=0.01) ;

%e2(dsn1=x2,dsn2=xx2, ovar=x2,chg=0.01) ;

%e2(dsn1=x3,dsn2=xx3, ovar=x3,chg=0.01) ;

%e2(dsn1=x4,dsn2=xx4, ovar=x4,chg=0.01) ;