%% Initial setup:  (This routine takes about 1h to solve)
clear all;filename5='specs';%Define how to call the temporary files with the specifications (named spec files in other scripts)
%% Define specs:
saveasfilename=v25_defaultvalues_specs(filename5);

%% Solve model: (Remove % if the model hasn't been solved yet, otherwise keep % to skip the solving and directly access the saved solution of the model which has been solved before)
% solvepar.startat=1;
% v25_call_solving(filename5,solvepar);
%This script saves the solution at the end.

%% Simulate model:
load(['results/' saveasfilename]);%Load the solution saved by the solving-script.
[SSIM,XSIM,ESIM,STAT]=v25_run_simulation(saveasfilename,filename5);%Run simulation
%save(['results/' saveasfilename '_simresults']);
%% Analyze the results:

%CALCULATE STATISTICAL DATA:
%Define percentiles and frequences which shall be calculated
for i=1:size(STAT,1) %Delete first 20 observations from XSIM
    XSIM(i).x=XSIM(i).x(:,:,21:end);
end;
pctls_values=[0.01 0.1 1 5 10 25 50 75 90 95 99 99.9 99.99];
freqs_range=[0 0.001; 0.035*0.999 0.035*1.001; 0.02*0.999 0.02*1.001];
%Calculate moments, summary statistics, percentiles and frequencies
[moments,pricestat,pctls,freqs]=v25_calcstatistics(STAT,XSIM,pctls_values,freqs_range);


% figure;boxplot_jb(squeeze(pricestat(:,1,:))','whisker',100);
% figure;boxplot_jb(squeeze(pricestat(:,1,[1 2 3 4 4]))','whisker',100);
% 
% figure;plot(squeeze(pricestat(:,1,1:4)))
% figure;plot(squeeze(pricestat(:,1,:)))
% figure;semilogy(squeeze(pricestat(:,1,1:4)))
% figure;semilogy(squeeze(pricestat(:,1,:)))

%CREATE OUTPUT TABLES AND FIGURES:
%Create table with moments:
%1=Mean; 2=Std. Dev.;3=Skewness; 4=Kurtosis; 5=Min; 6=Max; 7=%LB; 8=%UB
%Put moments for all response variables into table:
table=moments(:,:,3:end,:);
%Add percentiles:
%pctls_values=[0.01 0.1 1  5  10 25 50 75 90 95 99 99.9 99.99];
%               9   10  11 12 13 14 15 16 17 18 19 20   21
table(:,:,:,size(table,4)+1:size(table,4)+size(pctls,4))=pctls(:,:,:,:);
%Add frequences:
%freqs_range=[0 0.001; 0.0345 0.0355; 0.0195 0.0205];
%               22       23             24
table(:,:,:,size(table,4)+1:size(table,4)+size(freqs,4))=freqs(:,:,:,:);
%Add costs:
% 24: Costs from storage subsidy; 25: Costs from public reserve; 26: Total
% costs;
par2ka=[0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.06 0.03 0.00 0.03 0.04]; % THIS NEEDS TO BE UPDATED FOR ALL MODELS!
par2ka_sub=0.06-par2ka;
%This makes only sense for variable 3 - private stocks:
costs=table(:,:,3,1).*(0.06-par2ka(:));
%This makes only sense for variable 1 - public reserve:
costs(:,end+1)=table(:,:,1,1)*par2ka(1,1); 
costs(:,end+1)=costs(:,1)+costs(:,2);

%Calculate revenues from reserve's buying and selling operations and add to total costs:
reserveopcosts=v25_calcreserveopcosts(XSIM);
costs(:,3)=costs(:,3)+reserveopcosts(:);

%Calculate revenues from trade quota or from trade tariff: Calculated as 
%exports*price differences (The tariff phi cannot be used directly because 
%it is approximated only very poorly)
for i=1:size(XSIM,2)
    if i==7 || i==6
        quotarevenue(i)=sum(sum(max(XSIM(i).x(:,8,:,:)-XSIM(i).x(:,7,:,:)-par.thetaab,0).*XSIM(i).x(:,9,:,:)))...
            /(size(XSIM(i).x,1)*size(XSIM(i).x,3));
    else
        quotarevenue(i)=0;
    end;
end;
costs(:,3)=costs(:,3)-quotarevenue(:);%Add these revenues to total costs.
quotarevenue(:)' %Display costs

%Devide costs by agricultural GDP:
for i=1:size(XSIM,2)
    agriGDP(i,:)=sum(sum(XSIM(i).x(:,7,2:end,:).*XSIM(i).x(:,5,2:end,:).*ESIM(:,1,2:end,:)))/(size(XSIM(i).x,1)*size(XSIM(i).x,3)-size(XSIM(i).x,3));
end;
% No more division by Agri-GDP: costs(:,3)=costs(:,3)./agriGDP;

%Clear unnecessary variables and save the rest:
%clear ESIM SSIM XSIM interp model
%save(['results/' saveasfilename '_stats']);
%FIRST RUN BEGINNING AND THEN: load(['results/' saveasfilename '_stats'])

%Create tables
input.data=[squeeze(table(:,:,7,[1:4 11 15 17:20])) 100*costs(:,3)];
input.header=0;
input.tableColLabels={'\textbf{Mean}','\textbf{SD}','\textbf{Skew}','\textbf{Kurt}','\textbf{1\%}','\textbf{50\%}', ...
    '\textbf{90\%}','\textbf{95\%}','\textbf{99\%}','\textbf{99.9\%}','(in \%)'};
input.tableCol2Labels={'Prices','Costs'};
input.tableCol2multi=[10,1];
input.tableRowLabels={'1','2','3','4','5','6','7','8','9','10','11','12','13'};
input.dataFormat = {'%.3f',2,'%.2f',1,'%.1f',1,'%.2f',6,'%.3f',1};
input.display=0;
input.tableBorders = 1;%input.tableCaption = 'MyTableCaption';input.tableLabel = 'MyTableLabel';
input.filename='tables/price_defaultvalues.tex';
latexTable_jb2(input);

reservedata=[moments(:,:,1,1) moments(:,:,1,2) squeeze(table(:,:,1,[1 2 22 23])) 100*costs(:,3)];
reservedata([9 13],6)=squeeze(table([9 13],:,1,24));
reservedata(:,6)=1-reservedata(:,6);
input.tableColLabels={'\textbf{Mean}','\textbf{SD}','\textbf{Mean}','\textbf{SD}','\textbf{Freq(empty)}','\textbf{Freq(nun-full)}','(in \%)'};
input.tableCol2Labels={'Availability','Reserve','Costs'};
input.tableCol2multi=[2,4,1];
input.dataFormat = {'%.3f',7};%'%.2f',1};%,'%.3f',1,'%.2f',7};
input.data=reservedata;
input.filename='tables/reserve_defaultvalues.tex';
latexTable_jb2(input);

input.tableColLabels={'\textbf{Mean}','\textbf{SD}','\textbf{50\%}','\textbf{90\%}','\textbf{Mean}', ...
    '\textbf{SD}','\textbf{Mean}','\textbf{SD}','\textbf{90\%}','\textbf{Freq(trade)}'};
input.tableCol2Labels={'Private storage','Production','Exports'};
input.tableCol2multi=[4,2,4];
freqtrade=1-squeeze(table(:,:,9,[22]));
input.data=[squeeze(table(:,:,3,[1 2 15 17])) squeeze(table(:,:,5,[1 2])) squeeze(table(:,:,9,[1 2 17])) freqtrade];
input.dataFormat = {'%.3f',10};%,'%.3f',1,'%.2f',7};
input.filename='tables/combined_defaultvalues.tex';
latexTable_jb2(input);
input=rmfield(input,'tableCol2Labels');


%Create figures:

% 
%pctls_values=[0.01 0.1 1  5  10 25 50 75 90 95 99 99.9 99.99];
%               9   10  11 12 13 14 15 16 17 18 19 20   21
boxdata=squeeze(table(:,1,7,[11 13 15 17 19]));%95% and 99%
%boxdata(:,2)=table(:,2)-table(:,3);
%boxdata(:,3)=table(:,2);
figure;
boxplot_jb2015a(boxdata','whisker',100);
ylim([0.77 1.88]);
xlabel('Scenario');ylabel('Price');
ax1 = gca; % current axes
set(ax1,'XColor','k');set(ax1,'YColor','k');
set(ax1,'XTickLabel',{'1','2','3','4 (A)','4 (B)','5','6 (A)','6 (B)','7','8','11','12','13'});

% axes1 = axes('Parent',figure1,'YColor',[0 0 0],...
%     'XTickLabel',{'1','2','3','4 (A)','4 (B)','5','6 (A)','6 (B)','7','8','11','12','13'},...
%     'XTick',[1 2 3 4 5 6 7 8 9 10 11 12 13],...
%     'XColor',[0 0 0],...
%     'TickLabelInterpreter','none');

ax1_pos = get(ax1,'Position'); % position of first axes
ax2 = axes('Position',ax1_pos,'XAxisLocation','top','YAxisLocation','right','Color','none');
axes(ax2);hold on;
plot((1:13),100*costs(:,3),'r','LineStyle','none','Marker','X', ...
    'MarkerSize',15);ylabel('Costs (%)','Color','r');
xlim([0.5 13.5])
hold off;

export_fig 'figures/scenarios' -pdf -nocrop -transparent

% axes1=gca; 'YScale','log'
% 
% 
% axes1 = axes('Parent',gcf,'YScale','log');
% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[0.5 9.5]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[0.685753191312441 6.11795106813368]);
% box(axes1,'on');





% figure;pop=[1 1.9 5 5.5 15];
% boxplot_jb(pop,'whisker',100);