clear; 
clc;
addpath('G:\Research\JonathanH\EMX 2022 Replication Package\Model\Tables\CompEcon2017\CEtools')

format short g

params = zeros(2,14);

Markuptarget = 1.15;                        % our target for aggregate markup

% N            = 414.6473;                      % mean number of firms in sector
% p.sigma      = 1.6220;
% p.gamma      = 59.6930; 
% p.xi         = 28.0803; 
N            = 3425;                      % mean number of firms in sector
p.sigma      = 5.2;                         %  ls elast : 1.79
p.gamma      = 8.75;                          %  ls elast : 10.83
p.xi         = 7.54;                           %  ls elast : 8.128
cr4        = 0.40; 
cr20        = 0.59; 
selast_mu      = 0.10;
selast_ls      = -0.15;
sectors     = 440;
mu_spec     =0;

% Markuptarget = 1.25;                         % our target for aggregate markup
% 
% N            = 143.24;                       % mean number of firms in sector
% p.sigma      = 1.1485;
% p.gamma      = 7.1573; 
% p.xi         = 5.1451; 

% Markuptarget = 1.30; 
% 
% 
% N            = 109.39;                     % mean number of firms in sector
% p.sigma      = 1.0664;
% p.gamma      = 5.9662; 
% p.xi         = 4.3576; 
% 
% 
% Markuptarget = 1.35;                       % our target for aggregate markup
% 
% N            = 112.44;                     % mean number of firms in sector
% p.sigma      = 0.99424;
% p.gamma      = 5.2127; 
% p.xi         = 3.7245; 


x            = [N;     p.sigma;   p.gamma;   p.xi];
 
lb           = [3300;     1.1;       5.0;     5.0];  
ub           = [3500;     6.1;       80.0;     50.0];

ftarget      = @(x) objective(x, Markuptarget, cr4, cr20, selast_mu, selast_ls, mu_spec, sectors);

 
% ftarget(x)
% return

 
switch 'fminsearch' 

    case 'fminsearch'
        
        disp('fminsearch')
        
        options             = optimset('fminsearch');
        options.Display     = 'iter';
        options.TolX        = 1e-4;
        options.MaxFunEvals = 250;

        x    = fminsearchbnd(ftarget, x, lb, ub, options);
   
    
    case 'ga'

        disp('ga')
        gaoptions = gaoptimset('Display', 'off','UseParallel', 'always', 'InitialPopulation', x');
        x         = ga(@(x)ftarget(x), size(x, 1), [], [], [], [], lb, ub, [], gaoptions); 

  
    case 'particleswarm'
       
       disp('particleswarm')

       options = optimoptions('particleswarm', 'Display', 'off', 'MaxTime',  10000, 'UseParallel', true, 'InitialSwarm', x',  'SwarmSize', 200);

       x   = particleswarm(ftarget, numel(x), lb', ub', options);   %this function complains  if I give it a structure as input
        
       
    case 'patternsearch'
        
       disp('patternsearch')

       options = optimoptions('patternsearch','Display','off', 'UseParallel', true);

       x  = patternsearch(ftarget, x, [], [], [], [], lb, ub, [], options);
       
      
end

x  = x(:);
