%% Initial setup: %Takes 26,3h to run @50x50!
tic;
clear all;filename5='spec';
%% Define specs:
saveasfilename=v25_reservevariation_specs(filename5);

%% Solve model:
v25_call_solving(filename5);

%% Simulate model:
load(['results/' saveasfilename]);
toc;tic;
[SSIM,XSIM,ESIM,STAT]=v25_run_simulation(saveasfilename,filename5);
%save(['results/' saveasfilename '_simresults']);
%% Analyze the results:
for i=1:size(STAT,1) %Delete first 20 observations from XSIMsaveasfilename
    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; ([.005*0.999 0.005*1.001]'*[1:20])'];
[moments,pricestat,pctls,freqs]=v25_calcstatistics(STAT,XSIM,pctls_values,freqs_range);
toc;

%Create table with moments:
%1=Mean; 2=Std. Dev.;3=Skewness; 4=Kurtosis; 5=Min; 6=Max; 7=%LB; 8=%UB
%Put moments 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; 0.99 1.01];
%               22       23
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.03 0.00 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=table(:,:,1,1)*0.06; 
% costs(:,end+1)=costs(:,1)+costs(:,2);
[reserveopcosts]=v25_calcreserveopcosts(XSIM);
costs=costs+reserveopcosts;

%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=costs./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'])

pctls(:,:,7,13)
pctls(:,:,7,12)
xaxisp
yaxisp=[1:20]*0.005;
figure;contourf(pctls(:,:,7,12)');
figure;contourf(yaxisp,xaxisp,pctls(:,:,7,12)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;ylabel(c,'Price in 99.9 percentile');

figure;
contourf(yaxisp(1:17),xaxisp(1:16),pctls(1:17,1:16,7,12)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;ylabel(c,'Price in 99.9 percentile');

for i=1:size(freqs,1)
    freq_nonfull(i,:)=freqs(i,:,1,i+1);
end;
freq_nonfull=1-freq_nonfull;

pos_full=[0.05 0.71 0.34 0.28];

figure;
cutyend=0;
subplot(3,2,1);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,7,12)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;ylabel(c,'Price at 99.9%');
    %pos = get(gca, 'Position');
    pos(1) = 0.07;
    pos(2) = 0.72;
    pos(3) = 0.325;
    pos(4) = 0.26;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,2);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,7,11)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;ylabel(c,'Price at 99%');
    pos(1) = 0.55;
    pos(2) = 0.72;
    pos(3) = 0.325;
    pos(4) = 0.26;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,3);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),min(pctls(1:end,1:end-cutyend,7,9),1.2399)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;ylabel(c,'Price at 90%');
    pos(1) = 0.07;
    pos(2) = 0.38;
    pos(3) = 0.325;
    pos(4) = 0.26;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,4);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),100*costs(1:end,1:end-cutyend)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;ylabel(c,'Costs');
    pos(1) = 0.55;
    pos(2) = 0.38;
    pos(3) = 0.325;
    pos(4) = 0.26;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,5);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),freqs(1:end,1:end-cutyend,1,1)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;ylabel(c,'Frequency of empty reserve');
    pos(1) = 0.07;
    pos(2) = 0.05;
    pos(3) = 0.325;
    pos(4) = 0.26;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,6);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),freq_nonfull(1:end,1:end-cutyend)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;ylabel(c,'Frequency of non-full reserve');
    pos(1) = 0.55;
    pos(2) = 0.05;
    pos(3) = 0.325;
    pos(4) = 0.26;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
% saveas(gcf,'figures/reserve_saveas.pdf');
% export_fig 'figures/reservevariation' -pdf -a4 -r100 -RGB -opengl -nocrop -transparent

PP_old=get(gcf,'PaperPosition');
PS_old=get(gcf,'PaperSize');

set(gcf, 'color', [1 1 1]);set(gcf, 'Color', 'None');
% set(gca,'FontSize',30,'fontWeight','bold')
set(gca,'FontSize',14)%,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',14)%,'fontWeight','bold')

set(gcf, 'PaperPosition', [-0.16 -0.06 10.3 12.2]); %Position the plot further to the left and down. Extend the plot to fill entire paper.
set(gcf, 'PaperSize', [10 12]); %Keep the same paper size
saveas(gcf, 'figures/reservevariation', 'pdf')








figure;
cutyend=0;
subplot(3,2,1);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,7,12)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;title('Price at 99.9%');
    %pos = get(gca, 'Position');
    pos(1) = 0.07;
    pos(2) = 0.72;
    pos(3) = 0.325;
    pos(4) = 0.25;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,2);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,7,11)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;title('Price at 99%');
    pos(1) = 0.55;
    pos(2) = 0.72;
    pos(3) = 0.325;
    pos(4) = 0.25;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,3);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),min(pctls(1:end,1:end-cutyend,7,9),1.2399)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;title('Price at 90%');
    pos(1) = 0.07;
    pos(2) = 0.385;
    pos(3) = 0.325;
    pos(4) = 0.25;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,4);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),100*costs(1:end,1:end-cutyend)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;title('Costs');
    pos(1) = 0.55;
    pos(2) = 0.385;
    pos(3) = 0.325;
    pos(4) = 0.25;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,5);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),freqs(1:end,1:end-cutyend,1,1)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;title('Frequency of empty reserve');
    pos(1) = 0.07;
    pos(2) = 0.05;
    pos(3) = 0.325;
    pos(4) = 0.25;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
subplot(3,2,6);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),freq_nonfull(1:end,1:end-cutyend)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;title('Frequency of non-full reserve');
    pos(1) = 0.55;
    pos(2) = 0.05;
    pos(3) = 0.325;
    pos(4) = 0.25;
    set(gca, 'Position', pos);set(gca,'FontSize',14);
% saveas(gcf,'figures/reserve_saveas.pdf');
% export_fig 'figures/reservevariation' -pdf -a4 -r100 -RGB -opengl -nocrop -transparent

PP_old=get(gcf,'PaperPosition');
PS_old=get(gcf,'PaperSize');

set(gcf, 'color', [1 1 1]);set(gcf, 'Color', 'None');
% set(gca,'FontSize',30,'fontWeight','bold')
set(gca,'FontSize',14)%,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',14)%,'fontWeight','bold')

set(gcf, 'PaperPosition', [-0.18 -0.1 10.58 12.15]); %Position the plot further to the left and down. Extend the plot to fill entire paper.
set(gcf, 'PaperSize', [10 12]); %Keep the same paper size
saveas(gcf, 'figures/reservevariation2', 'pdf')


%% Crowding out of private stocks:
figure;
cutyend=0;
subplot(2,2,1);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,3,12)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;title('Private stocks at 99.9%');
    pos = get(gca, 'Position');
%     pos(1) = 0.14;
%     pos(2) = 0.66;
    pos(3) = 0.28;
%     pos(4) = 0.26;
    set(gca, 'Position', pos);
set(gca,'FontSize',14);
subplot(2,2,2);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,3,11)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;title('Private stocks at 99%');
    pos = get(gca, 'Position');
%     pos(1) = 0.55;
%     pos(2) = 0.72;
    pos(3) = 0.28;
%     pos(4) = 0.26;
     set(gca, 'Position', pos);
set(gca,'FontSize',14);
subplot(2,2,3);
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),pctls(1:end,1:end-cutyend,3,9)');
xlabel('Reserve capacity');ylabel('Trigger Price');
c=colorbar;title('Private stocks at 90%');
    pos = get(gca, 'Position');
%     pos(1) = 0.07;
%     pos(2) = 0.38;
     pos(3) = 0.28;
%     pos(4) = 0.26;
     set(gca, 'Position', pos);
set(gca,'FontSize',14);
subplot(2,2,4);%imagesc(moments(1:end,1:end-cutyend,5,1));
contourf(yaxisp(1:end),xaxisp(1:end-cutyend),moments(1:end,1:end-cutyend,5,1)');
xlabel('Reserve capacity');%ylabel('Trigger Price');
c=colorbar;title('Mean Private stocks');
    pos = get(gca, 'Position');
%     pos(1) = 0.55;
%     pos(2) = 0.38;
     pos(3) = 0.28;
%     pos(4) = 0.26;
     set(gca, 'Position', pos);
set(gca,'FontSize',14);

% saveas(gcf,'figures/reserve_saveas.pdf');
% export_fig 'figures/reservevariation' -pdf -a4 -r100 -RGB -opengl -nocrop -transparent

PP_old=get(gcf,'PaperPosition');
PS_old=get(gcf,'PaperSize');

set(gcf, 'color', [1 1 1]);set(gcf, 'Color', 'None');
% set(gca,'FontSize',30,'fontWeight','bold')
set(gca,'FontSize',14)%,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',14)%,'fontWeight','bold')

set(gcf, 'PaperPosition', [-0.9 -0.6 11.5 9.9]); %Position the plot further to the left and down. Extend the plot to fill entire paper.
set(gcf, 'PaperSize', [10 9]); %Keep the same paper size
saveas(gcf, 'figures/reservevariation_privstocks', 'pdf')
