%{
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % 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/>.   %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%}


function [y L]=probit(parameters,transfers,distance,balance)
  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;
  load amounts 

  save parameters parameters price_before price_after jac amounts

  alpha=parameters(47);
  beta=parameters(57);

  mu_d=abs(parameters(49));
  sigma_d=abs(parameters(50));
  gamma=parameters(48);

  mu_b1=abs(parameters(53));
  sigma_b1=abs(parameters(54));
  mu_b2=abs(parameters(55));
  sigma_b2=abs(parameters(56));

  bound=sqrt(1/(4*cos(pi/(33+1))^2)-1e-6);
  rho1=2*atan(parameters(58))/pi*bound;

  L=[parameters(59) 0 0;parameters(60) parameters(61) 0;parameters(62) parameters(63) parameters(64)];
  omega=L*L';
  sigma_alpha=sqrt(omega(1,1));
  sigma_beta=sqrt(omega(2,2));
  sigma_alpha_beta=omega(1,2)/(sigma_alpha*sigma_beta);

  fprintf('Transfer dummies:\n');
  for i=1:33
    fprintf('%f ',parameters(i));
  end
  fprintf('\n');
  fprintf('Time dummies:\n');
  for i=1:6
    fprintf('%f ',parameters(33+i));
  end
  fprintf('\n');
  for i=7:13
    fprintf('%f ',parameters(33+i));
  end
  fprintf('\n');

  fprintf('Price:\t %f, %f, %f, %f, %f\n',parameters(47), parameters(57), sigma_alpha, sigma_beta, sigma_alpha_beta);
  fprintf('Persistence:\t %f, %f\n',parameters(51),parameters(52));
  fprintf('Rho:\t %f\n',rho1);
  fprintf('Balance:\t %f, %f, %f, %f\n',mu_b1,sigma_b1,mu_b2,sigma_b2);
  fprintf('Distance\t %f, %f, %f\n', gamma, mu_d, sigma_d);

  system('cp output.mat output_save.mat');
  system('rm output.mat');

  system('./probit dataTran.mat parameters.mat output.mat ext.mat');


  load output

  for i=1:33
    fprintf('Choice %d: %f\t(%f)\t%f\t(%f)\n',i,moments(i),momentsData(i),moments(i+33),momentsData(i+33));
  end
  for i=1:14
    fprintf('Choice %d: %f\t(%f)\n',i,moments(66+i),momentsData(66+i));
  end
  for i=1:4
    fprintf('Distance %d: %f\t(%f)\n',i,moments(80+i),momentsData(80+i));
  end
  for i=1:3
    fprintf('Joint %d: ',i-1);
    for j=1:3
      fprintf('%f\t(%f)\t',moments(84+(i-1)*3+j),momentsData(84+(i-1)*3+j));
    end
    fprintf('\n');
  end
  for i=1:4
    fprintf('Balance %d: %f\t(%f)\n',i,moments(94+i),momentsData(94+i));
  end
  fprintf('Same transfer: %f (%f)\n',moments(93),momentsData(93));

  for i=1:4
    fprintf('Interactions %d: %f\t(%f)\n',i,moments(98+i),momentsData(98+i));
  end


  deleteMoments=[67 74 93];
  deleteParam=[];

  moments(deleteMoments)=[];momentsData(deleteMoments)=[];
 
  % Uncomment only of computing derivative
  %{
  J=(momentsGradient-repmat(moments,1,size(parameters,2)))*1e3;
  J(deleteMoments,:)=[];
  J(:,deleteParam)=[];
  %}

  load Omega

  f=moments-momentsData;
  y=f'*Omega*f;
  
  L=[];
  %L=2*f'*Omega*J;
end
