% *********************** 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)+g3*(t-k2)^3*(t>k2))

% *********************** 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, g2 and g3 from the restricted cubic spline trend : Ptrend(t) = exp(g1*t +g2*(t-k1)^3*(t>k1)+g3*(t-k2)^3*(t>k2))

[ComPrices,Com,Year,T,t,Residuals,Aobs,options,solveroptions,N,mina,maxa,r,delta] = OverallInitialization;
Tab        = NaN(length(Com),22);
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]};
RCS3       = load(fullfile('..','Results','RCS3Trend.mat'));
fgRCS3     = RCS3.Tableau{:,[1 3 5 7 9]};

options.ActiveParams = [1 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'));
  basis4            = importdata(fullfile('Basis','PricensBSM4.mat'));
  poly              = Legendre(t);
  ToOptim           = @(Target,BETA) sum((Target-Pobs.*exp(-[ones(T,1) basis4]*BETA)).^2);

  %% 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,3) fgparams;
                   0 0 0 1 0.3 0.02]; %#ok<*AGROW>

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

  %% Best Quadratic solution
  PstoQuadra        = Pobs.*exp(-poly(:,1:2)*fgQuadra(i,1:2)');
  lPobslmQuadra     = fitlm(basis4,poly(:,1:2)*fgQuadra(i,1:2)');
  beta              = fminsearch(@(BETA) ToOptim(PstoQuadra,BETA),lPobslmQuadra.Coefficients.Estimate);

  [model,interp]    = InitLogLikPb(Pobs,{'params' [beta(2:end)' fgQuadra(i,3:4) ...
                      delta fgQuadra(i,end)]'}, r, delta, mina, maxa, N, ...
                                   'RCS4',options);

  ToSolve = @(ML,BETA,FG,X) ML(i)*T-sum(LogLikelihood([BETA(2:end)' FG(i,(end-2):(end-1))*X delta FG(i,end)*X]',...
                                                    Pobs,model,interp,options));
  correct           = fminsearch(@(X) ToSolve(Quadra.Tableau.ML,beta,fgQuadra,X),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(basis4,poly*fgCubic(i,1:3)');
  beta              = fminsearch(@(BETA) ToOptim(PstoCubic,BETA),lPobslmCubic.Coefficients.Estimate);
  correct           = fminsearch(@(X) ToSolve(Cubic.Tableau.ML,beta,fgCubic,X),1);
  thetaCubic        = [beta(2:end)' fgCubic(i,4:5)*correct delta fgCubic(i,end)*correct];

  %% Best 3-knots RCS solution
  PstoRCS3          = Pobs.*exp(-basis3*fgRCS3(i,1:2)');
  lPobslmRCS3       = fitlm(basis4,basis3*fgRCS3(i,1:2)');
  beta              = fminsearch(@(BETA) ToOptim(PstoRCS3,BETA),lPobslmRCS3.Coefficients.Estimate);
  correct           = fminsearch(@(X) ToSolve(RCS3.Tableau.ML,beta,fgRCS3,X),1);
  thetaRCS3         = [beta(2:end)' fgRCS3(i,3:4)*correct delta fgRCS3(i,end)*correct];

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

  %% Likelihood maximization
  options.bounds.lb = [-5*ones(1,3)                0.5*mean(Psto) -20   -inf              0]';
  options.bounds.ub = [5*ones(1,3) 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
    NBstockouts             = length(find(Psto>=pstar));
    Residuals(:,Com{i})     = array2table(omega);
  end
  par                       = num2cell(theta);
  [g1, g2, g3, 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) g3 V(3) a V(4) b V(5) k V(6) ML AIC BIC NBstockouts ...
              pstar G NormG CondH exitflag exitflagFP ];
end

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