%% Code to calibrate parameters for oligopoly model

clear; 
clc;
addpath('\Model\Tables\CompEcon2017\CEtools')

format short g

% Create parametr saving object
params = zeros(4,4,2,2);
% Import data
pre_data=xlsread('\Inputs for model.xlsx','pre_olig','d2:k5');
post_data=xlsread('\Inputs for model.xlsx','post_olig','d2:k5');

% Toggle if use from labour share regression ,or mu regresion specfciation - preffered is labour sahre - 2
for mu_spec = 1:1
% Loop over pre and post period
for pre = 1:2

if pre == 1 
data = pre_data;
else
data=post_data;
end
% Loop over the various specficiations 1,2 are main versions in paper. 
for option = 1:2


% Put parameters into relevant places

   Markuptarget = data(option,3);                      % our target for aggregate markup

 N            =  data(option,1);                      % mean number of firms in sector

p.sigma      = 1;                         %  elast between
p.gamma      = 5;                          % elast wtihin
p.xi         = 4;                           %  tail
cr4        = data(option,7); 
cr20        = data(option,8); 
selast_mu      = data(option,2);
selast_ls      = -0.15;                     % old note check note and run version Ls
sectors     = 440;

x            = [N;     p.sigma;   p.gamma;   p.xi];
 
lb           = [3400;     0.9;       2;     2.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 = 500;

        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(:);
% SAave outputs

params(:,option,pre,mu_spec) = x;
end
end
end

% export
% 
% 
 writematrix(params(:,:,1,2), '\out_oli_ls_pre.xlsx')
 writematrix(params(:,:,2,2), '\out_oli_ls_post.xlsx')