% 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'));
quadratic = Tableau;
load(fullfile('..','Results','CubicTrend'));
cubic     = Tableau;

load(fullfile('..','Results','CompDataFeatures.mat'));
ActualDataFeatures      = rmfield(ActualDataFeatures,{'RCS3' 'RCS4'});
PredictedDataFeatures   = rmfield(PredictedDataFeatures,{'RCS3' 'RCS4'});
ActualNbStockouts       = rmfield(ActualNbStockouts,{'RCS3' 'RCS4'});
PredictedProbaStockouts = rmfield(PredictedProbaStockouts,{'RCS3' 'RCS4'});

%% Compare models performance with AIC, BIC, Loglik ratio, and AICc
BIC      = [without.BIC linear.BIC quadratic.BIC cubic.BIC];
AIC      = [without.AIC linear.AIC quadratic.AIC cubic.AIC];
NormAIC  = [without.AIC linear.AIC./AIC(:,1) quadratic.AIC./AIC(:,1)...
            cubic.AIC./AIC(:,1)];
LRStat   = 2*T*[(linear.ML-without.ML) (quadratic.ML-without.ML)...
                (cubic.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;

% Quantiles of the Chi2 distribution
Chi2Compa = NaN(3,3);
for i = 1:3
  Chi2Compa(i,:) = chi2inv([0.9 0.95 0.99],i+1);
end

%% 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','NormCompa.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);

% LR tests
fid = fopen(fullfile('..','Article','Tables','LRstat.tex'), 'w');
for i=1:13
  [~,indmin] = min(AIC(i,2:end));
  fprintf(fid,'%s',Com{i});
  for j=1:3
    if j~=indmin
      if LRStat(i,j)>=Chi2Compa(j,3)
        fprintf(fid,'&%4.2f&$^{\\ast \\ast \\ast}$',LRStat(i,j));
      elseif LRStat(i,j)>=Chi2Compa(j,2)
        fprintf(fid,'&%4.2f&$^{\\ast \\ast}$',LRStat(i,j));
      elseif LRStat(i,j)>=Chi2Compa(j,1)
        fprintf(fid,'&%4.2f&$^{\\ast}$',LRStat(i,j));
      else
        fprintf(fid,'&%4.2f&',LRStat(i,j));
      end
    end
    if j==indmin
      if LRStat(i,j)>=Chi2Compa(j,3)
        fprintf(fid,'&\\bf{%4.2f}&$^{\\ast \\ast \\ast}$',LRStat(i,j));
      elseif LRStat(i,j)>=Chi2Compa(j,2)
        fprintf(fid,'&\\bf{%4.2f}&$^{\\ast \\ast}$',LRStat(i,j));
      elseif LRStat(i,j)>=Chi2Compa(j,1)
        fprintf(fid,'&\\bf{%4.2f}&$^{\\ast}$',LRStat(i,j));
      else
        fprintf(fid,'&\\bf{%4.2f}&',LRStat(i,j));
      end
    end
  end
  fprintf(fid,'\\\\ \n');
end

%% Joint estimation results

% k in percentage, demand elasticity
Percentagek = 100*[without.k./without.a linear.k./linear.a quadratic.k./quadratic.a...
                   cubic.k./cubic.a];

Variancek = ([without.SD_k linear.SD_k quadratic.SD_k...
              cubic.SD_k]).^2;

Variancea = ([without.SD_a linear.SD_a quadratic.SD_a...
              cubic.SD_a]).^2;

Varianceb = ([without.SD_b linear.SD_b quadratic.SD_b...
              cubic.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).*(quadratic.k.^2)./(quadratic.a.^2))./(quadratic.a.^2)...
                    (Variancek(:,4)+Variancea(:,4).*(cubic.k.^2)./(cubic.a.^2))./(cubic.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

% change the values found to be null or too low with the restricted cubic splines standard errors (banana, rice, tea, and tin)
% CoefVar(CoefVar<=0.02) = Sigma_mu.CV_Rcs(CoefVar<=0.02);

ElastD        = [CoefVar.*without.a./without.b CoefVar.*linear.a./linear.b CoefVar.*quadratic.a./quadratic.b...
                 CoefVar.*cubic.a./cubic.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).*(quadratic.a.^2)./(quadratic.b.^2))./(quadratic.b.^2)...
                 (CoefVar.^2).*(Variancea(:,4)+Varianceb(:,4).*(cubic.a.^2)./(cubic.b.^2))./(cubic.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','Params1.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 quadratic.NBstockouts...
             cubic.NBstockouts];

fid = fopen(fullfile('..','Article','Tables','Params2.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={'PercentileNT.tex' 'PercentileBest.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={'PercentileNT20.tex' 'PercentileBest20.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={'ProbaStockoutsNT.tex' 'ProbaStockoutsBest.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={'ProbaStockoutsNTBest.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