%****************************** All trend specifications ******************************%
%                                                                                      %
% Calculate the predicted data features (mean, the two first order coefficents of      %
% autocorrelation, coefficent of variation, skewness and kurtosis) using the parameter %
% estimates of the storage model when trend is ignored and for all trend specifications%
%                                                                                      %
% Perform the percentiles of the distribution of the predicted data features           %
% on all consecutive subsamples of length T                                            %
%                                                                                      %
%**************************************************************************************%

[ComPrices,Com,Year,T,~,~,~,options,~,N,mina,maxa,r,delta] = OverallInitialization;

trend = {'Without' 'Linear' 'Quadratic' 'Cubic' 'RCS3' 'RCS4'};

ActualDataFeatures      = struct();
PredictedDataFeatures   = struct();
ActualNbStockouts       = struct();
PredictedProbaStockouts = struct();
VariableNames = {'Mean' 'AC_1','AC_2' 'CV' 'Skewness' 'Kurtosis' };

%% Define the grid of interpolation (equally space points)
interp.s           = linspace(mina,maxa,N)';
interp.PriceInterp = [];

%% Define the shocks structure by a discretization of the gaussian density
truncated_normal_rule(3,10,0,1,-5,5,fullfile('..','Results','GaussianQuadrature'))
model.shocks.e  = load(fullfile('..','Results','GaussianQuadrature_x.txt'));
model.shocks.w  = load(fullfile('..','Results','GaussianQuadrature_w.txt'));
TruncatedNormal = truncate(makedist('Normal'),-5,5);
model.funrand   = @(nrep,nper) TruncatedNormal.random(nrep,nper);

%% Number of periods and replications
Nsample         = 2E4;                % Number of replications in the MC experiment
burninper       = 30;                 % Number of burn-in periods
nper            = 49+T+burninper;     % Number of periods of simulation
A0              = zeros(Nsample,1);   % Initial vector of availability

disp('Number of samples of length T on the asymptotic distribution')
disp(Nsample*(nper+1-burninper-T))

%% Allocate memory
Tab           = NaN(length(Com),6);
Res           = NaN(length(Com),6);
Tabstouckouts = NaN(length(Com),2);
Resstockouts  = NaN(length(Com),6);
statdes       = NaN(Nsample*(nper+1-burninper-T),6);
NBstockouts   = NaN(Nsample*(nper+1-burninper-T),1);

for itrend=trend
  disp('Trend:')
  disp(itrend{1})

  %% Load the estimated detrended prices Psto and the estimated parameters
  load(fullfile('..','Results',[itrend{1} 'Prices.mat']))
  load(fullfile('..','Results',[itrend{1} 'Trend.mat']))
  switch itrend{1}
    %% Parametrization used to simulate the model
    case 'Without'
      target  = Tableau{:,[1 3 5]};
      Pcutoff = Tableau{:,11};
    case 'Linear'
      target  = Tableau{:,[3 5 7]};
      Pcutoff = Tableau{:,13};
    case 'Quadratic'
      target  = Tableau{:,[5 7 9]};
      Pcutoff = Tableau{:,15};
    case 'Cubic'
      target  = Tableau{:,[7 9 11]};
      Pcutoff = Tableau{:,17};
    case 'RCS3'
      target  = Tableau{:,[5 7 9]};
      Pcutoff = Tableau{:,15};
    case 'RCS4'
      target  = Tableau{:,[7 9 11]};
      Pcutoff = Tableau{:,17};
  end

  for i = 1:length(Com)
    disp(Com{i})

    if strcmp(itrend{1},'Without')
      Psto = Prices(i).Pobs; % When trend is ignored Psto = Pobs
    else
      Psto = Prices(i).Psto; % Load the cyclical data Psto
    end
    theta = target(i,:);
    pstar = Pcutoff(i);

    %% Simulate the model
    model.params  = [theta(1:2) delta theta(3) r];
    interp        = SolveStorageDL(model,interp,options);
    rng(options.seed)             % fix the seed of the random number generator
    [~,Xsim]      = SimulStorage(model,interp,A0,nper);
    P             = squeeze(Xsim(:,2,(burninper+1):end))';

    %% Statistics for all the possible subsamples of length T
    for t=1:(nper+1-burninper-T)
      statdes((1+(t-1)*Nsample):(t*Nsample),:) = [...
          mean(P(t:(t+T-1),:))' ...
          autocor(P(t:(t+T-1),:),2) ...
          (std(P(t:(t+T-1),:))./mean(P(t:(t+T-1),:)))' ...
          skewness(P(t:(t+T-1),:))' ...
          kurtosis(P(t:(t+T-1),:))'-3 ];
      NBstockouts((1+(t-1)*Nsample):(t*Nsample)) = sum(P(t:(t+T-1),:)>=pstar)';
    end

    %% Actual data features
    Tab(i,:) = [mean(Psto) autocor(Psto) std(Psto)/mean(Psto) ...
                skewness(Psto) kurtosis(Psto)-3];

    %% Percentiles
    Res(i,:) = diag(invprctile(statdes,Tab(i,:)));

    %% Predicted and actual number of stockouts
    Tabstouckouts(i,:) = [sum(Psto>=pstar) pstar];
    Resstockouts(i,:)  = sum([NBstockouts>=1 NBstockouts>=2 NBstockouts>=3 ...
                              NBstockouts>=4 NBstockouts>=5 NBstockouts>=10])/length(NBstockouts);
  end

  ActualFeatures        = array2table(Tab,'RowNames',Com,...
                                      'VariableNames',VariableNames);
  DataFeatures          = array2table(Res,'RowNames',Com,...
                                      'VariableNames',VariableNames);
  ActualDataFeatures    = setfield(ActualDataFeatures,itrend{1},ActualFeatures); %#ok<*SFLD>
  PredictedDataFeatures = setfield(PredictedDataFeatures,itrend{1}, ...
                                                 DataFeatures);

  ActualStockouts = array2table(Tabstouckouts,'RowNames',Com,...
                                'VariableNames',{'NbStockouts' 'Pstar' });
  ProbaStockouts  = array2table(Resstockouts,'RowNames',Com,...
                                'VariableNames',{'P1Stockouts' 'P2Stockouts' 'P3Stockouts'...
                      'P4Stockouts' 'P5Stockouts' 'P10Stockouts'});
  ActualNbStockouts       = setfield(ActualNbStockouts,itrend{1},ActualStockouts);
  PredictedProbaStockouts = setfield(PredictedProbaStockouts,itrend{1},ProbaStockouts);

end

save(fullfile('..','Results','CompDataFeatures.mat'),...
     'PredictedDataFeatures','ActualDataFeatures',...
     'ActualNbStockouts','PredictedProbaStockouts');