% *********************** THE MODEL and PARAMETERS to be estimated ***************************************
%
% The model is written as
%
% Pobs(t)   = f(Psto(t-1),z(t))
%
% The parameters conatained in the vector theta
%
% a and b<0 of the inverse demand function P(t) = a + b*c(t)
% k the marginal storage cost
% possibly the decay rate delta

[ComPrices,Com,Year,T,t,Residuals,Aobs,options,solveroptions,N,mina,maxa,r,delta] = OverallInitialization;
Tab      = NaN(length(Com),16);

options.ActiveParams = [1 1 0 1];
Nsample = options.solveroptions{1}.Size;

%% Run a loop to launch estimate for each commodity
for i = 1:length(Com)
  disp(Com{i})

  Pobs = table2array(ComPrices(:,Com{i}));

  %% Generate a set of first-guess parameters
  rng(options.seed)
  fgparamsrange  = [0.5   4;
                   0.03  0.5;
                   0    0.03]';
  fgparams       = fgparamsrange(ones(Nsample-1,1),:)+bsxfun(@times,fgparamsrange(2,:)-fgparamsrange(1,:),rand(Nsample-1,3));
  fgparams       = [fgparams;
                    1.2 0.2 0.01]; %#ok<*AGROW>

  thetasample    = [fgparams(:,1)*mean(Pobs) -std(Pobs)./fgparams(:,2) ...
                    zeros(Nsample,1) fgparams(:,3)*mean(Pobs)]';
  [model,interp] = InitLogLikPb(Pobs,{'params' thetasample(:,end)}, r, delta, ...
                                mina, maxa, N, 'without', options);

  %% Likelihood maximization
  options.bounds.lb = [0.5*mean(Pobs)                 -20   -inf 0]';
  options.bounds.ub = [max(10*mean(Pobs),3*max(Pobs)) -1E-3 +inf 0.1*mean(Pobs)]';

  [theta,ML,vcov,grad,H,exitflag,output] = ...
      MaxLik(@(theta,obs) LogLikVectorized(theta,obs,model,interp,options,false),...
             thetasample,Pobs,options);

  %% Postprocessing
  V                         = sqrt(diag(vcov));
  if isnan(ML)
    CondH                   = NaN;
    exitflagFP              = 0;
    omega                   = NaN(T,1);
    pstar                   = NaN;
    NBstockouts             = NaN;
  else
    [~,omega,pstar,Aobstmp,exitflagFP] = LogLikelihood(theta,Pobs,model,interp,options,true);
    Aobs(:,Com{i})          = array2table(Aobstmp);
    if all(isfinite(H(:)))
      CondH                 = max(eig(H))/min(eig(H));
    else
      CondH                 = NaN;
    end
    NBstockouts             = length(find(Pobs>=pstar));
  end
  par                       = num2cell(theta);
  [a, b, ~, k]              = par{:};
  G                         = norm(grad,'inf');
  NormG                     = grad'*(H\grad);
  Residuals(:,Com{i})       = array2table(omega);
  AIC                       = 2*(length(theta)-1-T*ML);
  BIC                       = -2*T*ML+(length(theta)-1)*log(T);

  %% Plot observed, simulated, and detrended price series, and the residuals
  figure
  plot(Year,[Pobs pstar(ones(T,1))]);
  h = legend('\fontname{Times} \fontsize{7} Pobs',...
             '\fontname{Times} \fontsize{7} p*');
  set(h, 'Box', 'off','Location','NorthEast');

  figure
  plot(Year,omega);

  Prices(i).commodity = Com{i}; %#ok<*SAGROW>
  Prices(i).Pobs      = Pobs;
  Prices(i).pstar     = pstar;

  Tab(i,:)    = [a V(1) b V(2) k V(3) ML AIC BIC NBstockouts pstar G NormG ...
                 CondH exitflag exitflagFP];
end

Tableau = array2table(Tab,'RowNames',Com,'VariableNames',{'a' 'SD_a' 'b' 'SD_b' ...
                    'k' 'SD_k' 'ML' 'AIC' 'BIC' 'NBstockouts' 'pstar' 'Gradient' ...
                    'Normalized_Grad' 'CondH' 'exitflag' 'exitflagFP'});
save(fullfile('..','Results','WithoutTrend.mat'),'Tableau','Residuals','Aobs');
save(fullfile('..','Results','WithoutPrices.mat'),'Prices','-v6')