% *********************** THE MODEL to be estimated ***************************************

% The model is written as
% Pobs(t)   = Ptrend(t)*Psto(t)
% Psto(t)   = f(Psto(t-1),z(t)) with z(t) normally iid with zero mean and unit variance
% Ptrend(t) =  exp(g1*t +g2*(t-k1)^3*(t>k1))

% *********************** THE PARAMETERS to be estimated ***************************************
% 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
% g1 and g2 from the restricted cubic spline trend : Ptrend(t) = exp(g1*t +g2*(t-k1)^3*(t>k1))

[ComPrices,Com,Year,T,t,Residuals,Aobs,options,solveroptions,N,mina,maxa,r,delta] = OverallInitialization;
Tab      = NaN(length(Com),20);
Lin      = load(fullfile('..','Results','LinearTrend.mat'));
fgLin    = Lin.Tableau{:,[1 3 5 7]};
Quadra   = load(fullfile('..','Results','QuadraticTrend.mat'));
fgQuadra = Quadra.Tableau{:,[1 3 5 7 9]};
Cubic    = load(fullfile('..','Results','CubicTrend.mat'));
fgCubic  = Cubic.Tableau{:,[1 3 5 7 9 11]};

options.ActiveParams = [1 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}));
  basis3 = importdata(fullfile('Basis','PricensBSM3.mat'));
  poly   = Legendre(t);

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

  %% First guess for the trend
  lPobslm        = fitlm(basis3,log(Pobs));
  g              = lPobslm.Coefficients.Estimate(2:end);
  Sigma2         = lPobslm.CoefficientCovariance;
  R              = chol(Sigma2(2:end,2:end));
  Psto           = Pobs.*exp(-basis3*g);
  trendparams    = repmat(g',Nsample-3,1)+5*fgparams(:,1:2)*R; % fgparams(1:2): normal distrib

  %% Best Linear solution
  Pstolin        = Pobs.*exp(-poly(:,1)*fgLin(i,1));
  MeanPtrendLin  = mean(exp(-poly(:,1)*fgLin(i,1)));
  lPobslmlin     = fitlm(basis3,poly(:,1)*fgLin(i,1));
  MeanPtrend     = mean(exp(-basis3*lPobslmlin.Coefficients.Estimate(2:end)));
  MeanCorrection = MeanPtrend/MeanPtrendLin;
  thetaLin       = [lPobslmlin.Coefficients.Estimate(2:end)' fgLin(i,2:3)*MeanCorrection delta fgLin(i,end)*MeanCorrection];

  %% Best Quadratic solution
  PstoQuadra       = Pobs.*exp(-poly(:,1:2)*fgQuadra(i,1:2)');
  lPobslmQuadra    = fitlm(basis3,poly(:,1:2)*fgQuadra(i,1:2)');
  ToOptim          = @(BETA) sum((PstoQuadra-Pobs.*exp(-[ones(T,1) basis3]*BETA)).^2);
  beta             = fminsearch(ToOptim,lPobslmQuadra.Coefficients.Estimate);
  [model,interp]   = InitLogLikPb(Pobs,{'params' [beta(2:end)' fgQuadra(i,3:4) ...
                      delta fgQuadra(i,end)]'}, r, delta, mina, maxa, N, 'RCS3', ...
                                  options);
  ToSolve          = @(X) Quadra.Tableau.ML(i)*T- sum(LogLikelihood([beta(2:end)' ...
                      fgQuadra(i,(end-2):(end-1))*X delta fgQuadra(i,end)*X]', ...
                                                    Pobs,model,interp, ...
                                                    options));
  correct          = fminsearch(ToSolve,1);
  thetaQuadra      = [beta(2:end)' fgQuadra(i,3:4)*correct delta fgQuadra(i,end)*correct];

  %% Best Cubic solution
  PstoCubic        = Pobs.*exp(-poly*fgCubic(i,1:3)');
  lPobslmCubic     = fitlm(basis3,poly*fgCubic(i,1:3)');
  ToOptim          = @(BETA) sum((PstoCubic-Pobs.*exp(-[ones(T,1) basis3]*BETA)).^2);
  beta             = fminsearch(ToOptim,lPobslmCubic.Coefficients.Estimate);
  [model,interp]   = InitLogLikPb(Pobs,{'params' [beta(2:end)' fgCubic(i,4:5) ...
                      delta fgCubic(i,end)]'}, r, delta, mina, maxa, N, 'RCS3', ...
                                  options);
  ToSolve          = @(X) Cubic.Tableau.ML(i)*T- sum(LogLikelihood([beta(2:end)' ...
                      fgCubic(i,(end-2):(end-1))*X delta fgCubic(i,end)*X]', ...
                                                    Pobs,model,interp, ...
                                                    options));
  correct          = fminsearch(ToSolve,1);
  thetaCubic       = [beta(2:end)' fgCubic(i,4:5)*correct delta fgCubic(i,end)*correct];

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

  %% Likelihood maximization
  options.bounds.lb     = [-5*ones(1,2) 0.5*mean(Psto) -20 -inf 0]';
  options.bounds.ub     = [+5*ones(1,2) max(10*mean(Psto),3*max(Pobs)) -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);
    NBstockouts         = length(find(Psto>=pstar));
  end
  par                   = num2cell(theta);
  [g1, g2, a, b, ~, k]  = par{:};
  G                     = norm(grad,'inf');
  NormG                 = grad'*(H\grad);
  Psto                  = Pobs.*exp(-grate);
  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) g2 V(2) a V(3) b V(4) k V(5) ML AIC BIC NBstockouts pstar ...
              G NormG CondH exitflag exitflagFP];
end

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