[ComPrices,Com,Year,T,t,Residuals,Aobs,options,solveroptions,N,mina,maxa,r,delta] = OverallInitialization;
Tab      = NaN(length(Com),18);
Wout   = load(fullfile('..','Results','WithoutTrend.mat'));
fgWout = Wout.Tableau{:,[1 3 5]};

options.ActiveParams = [1 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})
  rng(options.seed)
  Pobs = table2array(ComPrices(:,Com{i}));

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

  %% First guess for the trend
  poly            = Legendre(t);
  lPobslm         = fitlm(poly(:,1),log(Pobs));
  g               = lPobslm.Coefficients.Estimate(2:end);
  Sigma2          = lPobslm.CoefficientCovariance;
  Psto            = Pobs.*exp(-poly(:,1)*g);
  trendparams     = g+3*fgparams(:,1)*sqrt(Sigma2(2,2)); % fgparams(1): normal distrib

  %% Best solution without trend
  thetaWout     = [0 fgWout(i,1:2) delta fgWout(i,end)];

  %% Matrix of first guess
  thetasample     = [trendparams(:,1) fgparams(:,2)*mean(Psto,1) -std(Psto)./fgparams(:,3) ...
                     zeros(Nsample-2,1) fgparams(:,4)*mean(Psto);
                     thetaWout]';

  %% Likelihood maximization
  [model,interp]           = InitLogLikPb(Pobs,{'params' thetasample(:,end)}, ...
                                          r, delta, mina, maxa, N, 'linear', ...
                                          options);

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

  [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;
    exitflag                = 0;
    omega                   = NaN(T,1);
    pstar                   = NaN;
    NBstockouts             = NaN;
  else
   [~,omega,pstar,Aobstmp,exitflagFP,grate] = 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
    Residuals(:,Com{i})     = array2table(omega);
  end
  par                       = num2cell(theta);
  [g1, a, b, ~, k]          = par{:};
  G                         = norm(grad,'inf');
  NormG                     = grad'*(H\grad);
  Psto                      = Pobs.*exp(-grate);
  NBstockouts               = length(find(Psto>=pstar));
  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
  Ptrend = mean(Psto)*exp(grate);
  plot(Year,[Psto pstar(ones(T,1)) Pobs Ptrend]);
  h = legend ('\fontname{Times} \fontsize{7} Psto','\fontname{Times} \fontsize{7} p*',...
              '\fontname{Times} \fontsize{7} Pobs','\fontname{Times} \fontsize{7} Trend');
  set(h, 'Box', 'off','Location','NorthEast');

  figure
  plot(Year,omega);

  Prices(i).commodity = Com{i};
  Prices(i).Pobs      = Pobs;
  Prices(i).Psto      = Psto;
  Prices(i).Ptrend    = Ptrend;
  Prices(i).pstar     = pstar;

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

end

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