function [Tableau_CML, nb_discarded] = Simul_CML(type,MC)

%% Initialization of the Monte-Carlo experiment
switch MC
  case 1
    target = [1 -1 0 0.02];
  case 2
    target = [1 -2 0 0.02];
end
[model,interp,Psample,target,Inittheta,Nsample,options] = InitMonteCarlo(target,type);
options.estimator = 'CML';
discard           = 0;   % To count the number of discarded replications

%% Allocate memory
Tab             = NaN(Nsample,9);
MLsample        = NaN(Nsample,1);
thetasample     = NaN(Nsample,4);
gradsample      = NaN(Nsample,3);
NormGsample     = NaN(Nsample,1);
CondHsample     = NaN(Nsample,1);
exitflagSample  = NaN(Nsample,1);

%% Run a loop for each replication of the MC experiment
parfor iSample=1:Nsample

  Pobs  = Psample(iSample,:)';

  [theta,ML,~,grad,H,exitflag] = MaxLik(@(theta,obs) LogLikelihood(theta,obs,model,interp,options,false),...
                                        Inittheta(iSample,1:4)', ...
                                        Pobs,options);

  if isnan(ML)
    MLsample(iSample)       = NaN;
    thetasample(iSample,:)  = NaN(1,4);
    gradsample(iSample,:)   = NaN(1,3);
    NormGsample(iSample)    = NaN;
    CondHsample(iSample)    = NaN;
    exitflagSample(iSample) = 0;

  else
    MLsample(iSample)       = ML;
    thetasample(iSample,:)  = theta;
    gradsample(iSample,:)   = grad;
    NormGsample(iSample)    = grad'*(H\grad);
    if all(isfinite(H(:))), CondHsample(iSample) = max(eig(H))/min(eig(H));
    else                    CondHsample(iSample) = NaN;
    end
    exitflagSample(iSample) = exitflag;
  end
end

%% Postprocessing of results
for j = 1:Nsample
  Pobs        = Psample(j,:)';
  ML          = MLsample(j);
  theta       = thetasample(j,:)';
  grad        = gradsample(j,:);
  NormG       = NormGsample(j);
  CondH       = CondHsample(j);
  exitflag    = exitflagSample(j);


  if isnan(ML)
    pstar       = NaN;
    NBstockouts = NaN;
    G           = NaN;
  else
    model.params(1:4)    = theta;
    [~,~,pstar]          = SolveStorageDL(model,interp,options);
    par                  = num2cell(model.params);
    [a, b, k]            = par{[1:2 4]};
    G                    = norm(grad,'inf');

    NBstockouts = length(find(Pobs>=pstar)); % Number of stockouts
  end

  %% Discard the replications that failed
  if exitflag<=0
    Tab(j,:)    = NaN(1,9);
    discard     = discard+1;
  else
    Tab(j,:)    = [a  b  k  pstar NBstockouts G NormG CondH exitflag];
  end
end

target(:,3)      = [];
MeanTab          = nanmean(Tab(:,1:4));
SDTab            = nanstd(Tab(:,1:4));
Bias             = MeanTab-target;
RMSE             = sqrt(Bias.*Bias+SDTab.*SDTab);
nb_discarded     = discard;

Tableau_CML = array2table([MeanTab;SDTab;Bias;RMSE],...
                          'RowNames',{'Mean','St.dev.','Bias','RMSE'},...
                          'VariableNames',{'a' 'b' 'k' 'pstar'});

save(fullfile('..','..','Results','MonteCarlo',...
              ['Results' num2str(MC) '_CML' num2str(size(Psample,2)) '.mat']),...
     'Tableau_CML','Tab','nb_discarded')
