% This function transcript the model parameters a, b and k in percentage of
% storage cost and in price elasticity of demand

Com = {'Banana' 'Cocoa' 'Coffee'  'Copper'  'Cotton' 'Jute' 'Maize' 'Palm oil' 'Rice'  'Sugar'  'Tea' 'Tin' 'Wheat'};
T = 112;

%% Load results
load(fullfile('..','Results','WithoutTrend.mat'));
without   = Tableau;
load(fullfile('..','Results','LinearTrend.mat'));
linear    = Tableau;
load(fullfile('..','Results','QuadraticTrend.mat'));
quadratic = Tableau;
load(fullfile('..','Results','CubicTrend.mat'));
cubic     = Tableau;
load(fullfile('..','Results','RCS3Trend.mat'));
RCS3Trend = Tableau;
load(fullfile('..','Results','RCS4Trend.mat'));
RCS4Trend = Tableau;
load(fullfile('..','Results','CompDataFeatures.mat'));

ActualDataFeatures      = rmfield(ActualDataFeatures,{'Quadratic' 'Cubic'});
PredictedDataFeatures   = rmfield(PredictedDataFeatures,{'Quadratic' 'Cubic'});
ActualNbStockouts       = rmfield(ActualNbStockouts,{'Quadratic' 'Cubic'});
PredictedProbaStockouts = rmfield(PredictedProbaStockouts,{'Quadratic' 'Cubic'});

%% Compare models performance with AIC, BIC, Loglik ratio, and AICc
BIC      = [without.BIC linear.BIC RCS3Trend.BIC RCS4Trend.BIC];
AIC      = [without.AIC linear.AIC RCS3Trend.AIC RCS4Trend.AIC];
NormAIC  = [without.AIC linear.AIC./AIC(:,1) RCS3Trend.AIC./AIC(:,1)...
            RCS4Trend.AIC./AIC(:,1)];
LRStat   = 2*T*[(linear.ML-without.ML) (RCS3Trend.ML-without.ML)...
                (RCS4Trend.ML-without.ML)];
Penalize = repmat([2*3*4/(T-3-1) 2*4*5/(T-4-1) 2*5*6/(T-5-1) 2*6*7/(T-6-1)],13,1);
AICc     = AIC + Penalize;

%% Matrix of critical values for the inverse of the chi-square cdf with degrees of freedom (1 for linear, 2 for RCS3, and 3 for RCS4)
% and levels alpha= 1,5, and 10%
DoF     = repmat([1 2 3],3,1);
levels  = repmat([0.9 0.95 0.99]',1,3);
CritVal = chi2inv(levels,DoF);
%% Export Table Results

% Without
without.ML = without.ML*T;

ExportResults(without,'Without',Com,0,true)

% Linear
linear.ML    = linear.ML*T;

% Express g1 in annual growth rates
t            = (1:1:T)'/T;
lt           = min(t);
ut           = max(t);

g1           = linear.g1;
sdg1         = linear.SD_g1;
gnew         = (2/(ut-lt)).*g1;
sdgnew       = (2/(ut-lt)).*sdg1;
G            = gnew/T;
SD_G         = sdgnew/T;

G            = table(G);
SD_G         = table(SD_G);
linear       = [G SD_G linear];


ExportResults(linear,'Linear',Com,1,false)

% Quadratic
quadratic.ML = quadratic.ML*T;


ExportResults(quadratic,'Quadratic',Com,2,false)

% Cubic
cubic.ML = cubic.ML*T;

ExportResults(cubic,'Cubic',Com,3,false)

% RCS3
RCS3Trend.ML = RCS3Trend.ML*T;

ExportResults(RCS3Trend,'RCS3',Com,2,false)

% RCS4
RCS4Trend.ML = RCS4Trend.ML*T;

ExportResults(RCS4Trend,'RCS4',Com,3,false)

%% Export AIC

AIC = AIC(:,1:end);
BestModels = NaN(13,1);
for i=1:13, [~,BestModels(i)] = min(AIC(i,:)); end

fid = fopen(fullfile('..','Article','Tables','RCSNormCompa.tex'), 'w');
for i=1:13
  [~,indmin] = min(AIC(i,:));
  fprintf(fid,'%s',Com{i});
  for j=1:4
    if j~=indmin
      fprintf(fid,'&%4.3f',NormAIC(i,j));
    else
      fprintf(fid,'&\\mathbf{%4.3f}',NormAIC(i,j));
    end
  end
  fprintf(fid,'\\\\ \n');
end
fclose(fid);

%% Export LRStat

LRStat = LRStat(:,1:end);

fid = fopen(fullfile('..','Article','Tables','LRStat.tex'), 'w');
for i=1:13
  indk       = [without.k(i) linear.k(i) RCS3Trend.k(i) RCS4Trend.k(i)];
  [~,indmin] = min(AIC(i,:));
  fprintf(fid,'%s',Com{i});
  for j=1:3
      if all(indk)~=0 && i~=12
        if j~=(indmin-1) % because LR stat has only three columns while AIC has 4;
           if CritVal(1,j)>LRStat(i,j)
              fprintf(fid,'&%4.2f',LRStat(i,j));
           elseif CritVal(2,j)>LRStat(i,j)
              fprintf(fid,'&%4.2f$^{\\ast }$',LRStat(i,j));
           elseif CritVal(3,j)>LRStat(i,j)
              fprintf(fid,'&%4.2f$^{\\ast \\ast }$',LRStat(i,j));
           else
              fprintf(fid,'&%4.2f$^{\\ast \\ast \\ast }$',LRStat(i,j));
          end
        else
            if CritVal(1,j)>LRStat(i,j)
              fprintf(fid,'&\\bf{%4.2f}',LRStat(i,j));
            elseif CritVal(2,j)>LRStat(i,j)
              fprintf(fid,'&\\bf{%4.2f$^{\\ast }$}',LRStat(i,j));
            elseif CritVal(3,j)>LRStat(i,j)
              fprintf(fid,'&\\bf{%4.2f$^{\\ast \\ast }$}',LRStat(i,j));
            else
              fprintf(fid,'&\\bf{%4.2f$^{\\ast \\ast \\ast }$}',LRStat(i,j));
            end
        end
      elseif i ==12&& j==1
           if CritVal(1,j)>LRStat(i,j)
              fprintf(fid,'&%4.2f',LRStat(i,j));
           elseif CritVal(2,j)>LRStat(i,j)
              fprintf(fid,'&%4.2f$^{\\ast }$',LRStat(i,j));
           elseif CritVal(3,j)>LRStat(i,j)
              fprintf(fid,'&%4.2f$^{\\ast \\ast }$',LRStat(i,j));
           else
              fprintf(fid,'&%4.2f$^{\\ast \\ast \\ast }$',LRStat(i,j));
          end
      else
          fprintf(fid,'&--');
      end
  end
  fprintf(fid,'\\\\ \n');
end
fclose(fid);
%% Export k in percentage, demand elasticity, number of stockouts and cutoff price

% Joint estimation results

% k in percentage, demand elasticity
Percentagek = 100*[without.k./without.a linear.k./linear.a ...
                   RCS3Trend.k./RCS3Trend.a RCS4Trend.k./RCS4Trend.a];
Variancek = [without.SD_k linear.SD_k RCS3Trend.SD_k...
             RCS4Trend.SD_k].^2;
Variancea = [without.SD_a linear.SD_a RCS3Trend.SD_a...
             RCS4Trend.SD_a].^2;
Varianceb = [without.SD_b linear.SD_b RCS3Trend.SD_b...
             RCS4Trend.SD_b].^2;

% Variance of k in percentage of the mean price: (100/a)^2*[V(k)+(k/a)^2*V(a)]
VPercentagek  = (100.^2).*[(Variancek(:,1)+Variancea(:,1).*(without.k.^2)./(without.a.^2))./(without.a.^2)...
                           (Variancek(:,2)+Variancea(:,2).*(linear.k.^2)./(linear.a.^2))./(linear.a.^2)...
                           (Variancek(:,3)+Variancea(:,3).*(RCS3Trend.k.^2)./(RCS3Trend.a.^2))./(RCS3Trend.a.^2)...
                           (Variancek(:,4)+Variancea(:,4).*(RCS4Trend.k.^2)./(RCS4Trend.a.^2))./(RCS4Trend.a.^2)];

SDPercentagek = sqrt(VPercentagek);

% load the standard error of supply shocks taken from the detrended data quantity (1961-2011)
Sigma_mu      = readtable(fullfile('..','Results','CVQuantities.csv'),...
                          'Delimiter','\t','ReadRowNames',true);
CoefVar       = Sigma_mu.Production; % get the coefficient of variation (sigma/mu) detrending with the restricted cubic splines

ElastD        = [CoefVar.*without.a./without.b CoefVar.*linear.a./linear.b ...
                 CoefVar.*RCS3Trend.a./RCS3Trend.b CoefVar.*RCS4Trend.a./RCS4Trend.b];

% Variance of demand elasticity eta:(sigma/(b*mu))^2*[V(a)+(a/b)^2*V(b)]
VElastD       = [(CoefVar.^2).*(Variancea(:,1)+Varianceb(:,1).*(without.a.^2)./(without.b.^2))./(without.b.^2)...
                           (CoefVar.^2).*(Variancea(:,2)+Varianceb(:,2).*(linear.a.^2)./(linear.b.^2))./(linear.b.^2)...
                           (CoefVar.^2).*(Variancea(:,3)+Varianceb(:,3).*(RCS3Trend.a.^2)./(RCS3Trend.b.^2))./(RCS3Trend.b.^2)...
                           (CoefVar.^2).*(Variancea(:,4)+Varianceb(:,4).*(RCS4Trend.a.^2)./(RCS4Trend.b.^2))./(RCS4Trend.b.^2)];

SDElastD = sqrt(VElastD);
ind = sub2ind([13 4],(1:13)',BestModels);
% Percentagek = [Percentagek(:,1) Percentagek(ind)]
Result   = [Percentagek ElastD];
SDResult = [SDPercentagek SDElastD];

fid = fopen(fullfile('..','Article','Tables','RCSParams1.tex'), 'w');
for i=1:13
  [~,indmin] = min(AIC(i,:));
  fprintf(fid,'%s',Com{i});
  for j=1:4
    if j~=indmin
      if Result(i,j)==0
        fprintf(fid,'&%1g',Result(i,j));
      else
        fprintf(fid,'&%4.2f',Result(i,j));
      end
    else
      if Result(i,j)==0
        fprintf(fid,'&\\bf{%1g}',Result(i,j));
      else
        fprintf(fid,'&\\bf{%4.2f}',Result(i,j));
      end
    end
  end
  for j= 5:8
    if j~=indmin+4
      if Result(i,j)==0
        fprintf(fid,'&$%1g$',Result(i,j));
      else
        fprintf(fid,'&$%4.3f$',Result(i,j));
      end
    else
      if Result(i,j)==0
        fprintf(fid,'&$\\bf{%1g}$',Result(i,j));
      else
        fprintf(fid,'&$\\bf{%4.3f}$',Result(i,j));
      end
    end
  end
  fprintf(fid,'\\\\ \n');
  for j =1:4
    if isnan(SDResult(i,j))
      fprintf(fid,'&--');
    else
      fprintf(fid,'&(%4.2f)',SDResult(i,j));
    end
  end
  for j= 5:8
    fprintf(fid,'&$(%4.3f)$',SDResult(i,j));
  end
  fprintf(fid,'\\\\ \n');
end
fclose(fid);

% Effect of storage on the price standard deviation and Coefficient of variation
% of the Consumption and number of stockouts
Stockouts = [without.NBstockouts linear.NBstockouts RCS3Trend.NBstockouts...
             RCS4Trend.NBstockouts];

fid = fopen(fullfile('..','Article','Tables','RCSParams2.tex'), 'w');
for i=1:13
  [~,indmin] = min(AIC(i,:));
  fprintf(fid,'%s',Com{i});
  for j=1:4
    if j~=indmin
      fprintf(fid,'&$%1g$',Stockouts(i,j));
    else
      fprintf(fid,'&$\\mathbf{%1g}$',Stockouts(i,j));
    end
  end
  fprintf(fid,'\\\\ \n');
end
fclose(fid);

%% Comparison of statistics
ActualDataFeatures    = structfun(@table2array,ActualDataFeatures, ...
                               'UniformOutput',false);
ActualDataFeatures    = struct2cell(ActualDataFeatures);
ActualDataFeatures    = cell2mat(ActualDataFeatures);
ActualDataFeatures    = reshape(ActualDataFeatures,[13 4 6]);

PredictedDataFeatures = structfun(@table2array,PredictedDataFeatures, ...
                               'UniformOutput',false);
PredictedDataFeatures = struct2cell(PredictedDataFeatures);
PredictedDataFeatures = cell2mat(PredictedDataFeatures);
PredictedDataFeatures = reshape(PredictedDataFeatures,[13 4 6]);

Actual  = NaN(13,2,6);
Actual(:,1,:)           = ActualDataFeatures(:,1,:);
for com = 1:13
  Actual(com,2,:)       = ActualDataFeatures(com,BestModels(com),:);
end

Predicted = NaN(13,2,6);
Predicted(:,1,:)        = PredictedDataFeatures(:,1,:);
for com   = 1:13
  Predicted(com,2,:)    = PredictedDataFeatures(com,BestModels(com),:);
end

% iter = 0;
% for file={'RCSPercentileNT.tex' 'RCSPercentileBest.tex'}
%   iter = iter+1;
%   statistics = squeeze(Actual(:,iter,:));
%   percentile = abs(squeeze(Predicted(:,iter,:))-50);
%   fid = fopen(fullfile('..','Article','Tables',file{1}), 'w');
%   for i=1:13
%     fprintf(fid,'%s',Com{i});
%     for j=1:6
%       if percentile(i,j)<40
%         fprintf(fid,'&%4.2f',statistics(i,j));
%       elseif percentile(i,j)<45
%         fprintf(fid,'&%4.2f^{\\bullet} ',statistics(i,j));
%       elseif percentile(i,j)<47.5
%         fprintf(fid,'&%4.2f^{\\ast} ',statistics(i,j));
%       elseif percentile(i,j)<49.5
%         fprintf(fid,'&%4.2f^{\\ast \\ast} ',statistics(i,j));
%       else
%         fprintf(fid,'&%4.2f^{\\ast \\ast \\ast}',statistics(i,j));
%       end
%     end
%     fprintf(fid,'\\\\ \n');
%   end
%   fclose(fid);
% end

iter = 0;
for file={'RCSPercentileNT.tex' 'RCSPercentileBest.tex'}
  iter = iter+1;
  statistics = squeeze(Actual(:,iter,:));
  percentile = squeeze(Predicted(:,iter,:));
  fid = fopen(fullfile('..','Article','Tables',file{1}), 'w');
  for i=1:13
    fprintf(fid,'%s',Com{i});
    fprintf(fid,'&%s','Observed moments');
    for j=1:6
        fprintf(fid,'&%4.2f',statistics(i,j));
    end
    fprintf(fid,'\\\\ \n');
    fprintf(fid,'&%s','Model percentiles');
    for j =1:6
      fprintf(fid,'&%4.2f',percentile(i,j));
    end
  fprintf(fid,'\\\\ \n');
  end
  fclose(fid);
end

iter = 0;
for file={'RCSPercentileNT20.tex' 'RCSPercentileBest20.tex'}
  iter = iter+1;
  statistics = squeeze(Actual(:,iter,:));
  percentile = abs(squeeze(Predicted(:,iter,:))-50);
  fid = fopen(fullfile('..','Article/Tables',file{1}), 'w');
  for i=1:13
    fprintf(fid,'%s',Com{i});
    for j=1:6
      if percentile(i,j)<40
        fprintf(fid,'&%4.2f',statistics(i,j));
      elseif percentile(i,j)<47.5
        fprintf(fid,'&%4.2f^{\\ast}',statistics(i,j));
      elseif percentile(i,j)<49.5
        fprintf(fid,'&%4.2f^{\\ast \\ast} ',statistics(i,j));
      else
        fprintf(fid,'&%4.2f^{\\ast \\ast \\ast}',statistics(i,j));
      end
    end
    fprintf(fid,'\\\\ \n');
  end
  fclose(fid);
end

%% Comparison of implied number of stockouts
ActualNbStockouts = structfun(@table2array,ActualNbStockouts, ...
                               'UniformOutput',false);
ActualNbStockouts = struct2cell(ActualNbStockouts);
ActualNbStockouts = cell2mat(ActualNbStockouts);
ActualNbStockouts = reshape(ActualNbStockouts,[13 4 2]);

PredictedProbaStockouts = structfun(@table2array,PredictedProbaStockouts, ...
                               'UniformOutput',false);
PredictedProbaStockouts = struct2cell(PredictedProbaStockouts);
PredictedProbaStockouts = cell2mat(PredictedProbaStockouts);
PredictedProbaStockouts = reshape(PredictedProbaStockouts,[13 4 6]);

Actual = NaN(13,2,2);
Actual(:,1,:) = ActualNbStockouts(:,1,:);
for com=1:13
  Actual(com,2,:) = ActualNbStockouts(com,BestModels(com),:);
end

Predicted = NaN(13,2,6);
Predicted(:,1,:) = PredictedProbaStockouts(:,1,:);
for com=1:13
  Predicted(com,2,:) = PredictedProbaStockouts(com,BestModels(com),:);
end

iter = 0;
for file={'RCSProbaStockoutsNT.tex' 'RCSProbaStockoutsBest.tex'}
  iter = iter+1;
  statistics = squeeze(Predicted(:,iter,:));
  fid = fopen(fullfile('..','Article/Tables',file{1}), 'w');
  for i=1:13
    fprintf(fid,'%s',Com{i});
    for j=1:6
      fprintf(fid,'&$%4.2f$',statistics(i,j));
    end
    fprintf(fid,'\\\\ \n');
  end
  fclose(fid);
end

for file={'RCSProbaStockoutsNTBest.tex'}
  statistics = [squeeze(Predicted(:,1,:)) squeeze(Predicted(:,2,:))] ;
  fid = fopen(fullfile('..','Article','Tables',file{1}), 'w');
  for i=1:13
    fprintf(fid,'%s',Com{i});
    for j=[1:3 5:9 11:12]
      fprintf(fid,'&$%4.2f$',statistics(i,j));
    end
    fprintf(fid,'\\\\ \n');
  end
  fclose(fid);
end

%% Combining nb of stockouts and impled number of stockouts
% Nb of stockouts for no trend and preferred models
NBStockouts = [Stockouts(:,1) Stockouts(sub2ind(size(Stockouts),...
                                                (1:size(Stockouts,1))', ...
                                                BestModels))];
% Probabilities of 1, 2, 5, and 10 stockouts
ProbImplied = reshape(permute(Predicted(:,:,[1:2 5:6]),[1 3 2]),13,8);
% Export to LaTeX
fid = fopen(fullfile('..','Article','Tables','RCSStockouts.tex'), 'w');
for i=1:13
  fprintf(fid,'%s',Com{i});
  fprintf(fid,'&$%1g$',NBStockouts(i,1));
  for j=1:4
    fprintf(fid,'&$%4.2f$',ProbImplied(i,j));
  end
  fprintf(fid,'&$%1g$',NBStockouts(i,2));
  for j=5:8
    fprintf(fid,'&$%4.2f$',ProbImplied(i,j));
  end
  fprintf(fid,'\\\\ \n');
end
fclose(fid);

