%{
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % This package contains code for the paper "Mobile Money in Tanzania" by  %
 %                                                                         %
 % Nicholas Economides                                                     %
 % Stern School of Business, NYU; NET Institute; economides@stern.nyu.edu  %
 %                                                                         %
 % and                                                                     %
 %                                                                         %
 % Przemyslaw Jeziorski                                                    %
 % Haas School of Business, UC Berkeley; przemekj@haas.berkeley.edu        %
 %                                                                         %
 % Copyright by Przemyslaw Jeziorski, 2016                                 %
 %                                                                         %
 % This program is free software: you can redistribute it and/or modify    %
 % it under the terms of the GNU General Public License as published by    %
 % the Free Software Foundation, either version 3 of the License, or       %
 % (at your option) any later version.                                     %
 %                                                                         %
 % This program is distributed in the hope that it will be useful,         %
 % but WITHOUT ANY WARRANTY; without even the implied warranty of          %
 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the           %
 % GNU General Public License for more details.                            %
 %                                                                         %
 % You should have received a copy of the GNU General Public License       %
 % along with this program.  If not, see <http://www.gnu.org/licenses/>.   %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%}


load amounts

load parametersFinal

price_before=[0, 0.020, 0.030, 0.040, 0.050, 0.050, 0.050, 0.050, 0.050, .050, .050, .200, .200, .200, .200, .200, .200, .200, .200, .500, .500, .500, .500,...
             .500, .500, .500, 1.000, 1.000, 1.500, 1.500, 1.500, 1.500, 1.500];
price_after=[0, .020, .030, .040, .050, .050, .100, .100, .100, .100, .100, .200, .200, .250, .250, .250, .250, .250, .250, .500, .500, .500, .500,...
           .500, .500, .500, 1.000, 1.000, 1.500, 1.500, 1.500, 1.500, 1.500];

jac=0;

% In the code there is minus before 59,60,61, thus we need to put minus to get conditional mean and
% then put - again when we feed it to the code
L=[-parameters(59) 0 0;-parameters(60) -parameters(61) 0;parameters(62) parameters(63) parameters(64)];
omega2=L*L';
printMatrix(omega2);

mu1=[parameters(47);parameters(57)];
mu2=0;
sigma11=omega2(1:2,1:2);
sigma21=omega2(3,1:2);
sigma12=omega2(1:2,3);
sigma22=omega2(3,3);

D=zeros(33,5);
mud=zeros(2,5);
i=1;
for z=[-3 -1 0 1 3]
  load parametersFinal

  jac=0;

  price_before=[0, 0.020, 0.030, 0.040, 0.050, 0.050, 0.050, 0.050, 0.050, .050, .050, .200, .200, .200, .200, .200, .200, .200, .200, .500, .500, .500, .500,...
               .500, .500, .500, 1.000, 1.000, 1.500, 1.500, 1.500, 1.500, 1.500];
  price_after=[0, .020, .030, .040, .050, .050, .100, .100, .100, .100, .100, .200, .200, .250, .250, .250, .250, .250, .250, .500, .500, .500, .500,...
             .500, .500, .500, 1.000, 1.000, 1.500, 1.500, 1.500, 1.500, 1.500];

  a=z*sqrt(omega2(3,3));

  mu=mu1+sigma12*inv(sigma22)*(a-mu2)
  sigma=sigma11-sigma12*inv(sigma22)*sigma21
  L2=chol(sigma)';

  parameters(47)=mu(1);
  parameters(57)=mu(2);

  mud(:,i)=mu;

  parameters(59)=-L2(1,1);parameters(60)=-L2(2,1);parameters(61)=-L2(2,2);parameters(62:64)=0;
  parameters(1:33)=parameters(1:33)+a*amounts;

  save parametersEl parameters price_before price_after jac amounts

  system('rm outputEl.mat');
  system('./probit dataTran.mat parametersEl.mat outputEl.mat ../ext.mat');

  load outputEl
  D(:,i)=moments(1:33)/sum(moments(1:33));
  d1=sum(moments(1:33));

  price_before=price_before*1.01;

  save parametersEl parameters price_before price_after jac amounts

  system('rm outputEl.mat');
  system('./probit dataTran.mat parametersEl.mat outputEl.mat ../ext.mat');
  load outputEl
  d2=sum(moments(1:33));

  el_t(i)=(d2-d1)/d1*100;

  price_before=[0, 0.020, 0.030, 0.040, 0.050, 0.050, 0.050, 0.050, 0.050, .050, .050, .200, .200, .200, .200, .200, .200, .200, .200, .500, .500, .500, .500,...
    .500, .500, .500, 1.000, 1.000, 1.500, 1.500, 1.500, 1.500, 1.500];

  save parametersEl parameters price_before price_after jac amounts

  load ../ext.mat
  ext_before=ext_before*1.01;
  save ext_el ext_before ext_after

  system('rm outputEl.mat');
  system('./probit dataTran.mat parametersEl.mat outputEl.mat ext_el.mat');
  load outputEl
  d3=sum(moments(1:33));

  el_c(i)=(d3-d1)/d1*100;

  i=i+1;
end

save heterogeneity D el_c el_t mud
