clc
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Codes for PLOT 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Introducing carbon tax in Italy: Is there room for a quadruple-dividend effect?
% Journal: Energy Economics 2023

%% Import data from Vensim simulations

% Base = importdata('BASE.tab');
S_base = importdata('Base_sim1000.tab'); %Baseline as PNIEC plan
S_CT0 = importdata('CT_sim1000.tab'); % Carbon tax with redistribution
S_CTbeta = importdata('CTb_sim1000.tab'); % Consumer reaction to CT
S_CTbeta_gamma = importdata('CTbg_sim1000.tab'); % Firms reaction to CT

all_variables = S_base.textdata';  
all_variables(1:2) = [];
seed = S_CTbeta.data(:,2);
variab = S_CTbeta.textdata';
variab(1:2) = [];

S_CTbeta.data(:,1:2) = []; %deleting seeds and label
S_CTbeta_gamma.data(:,1:2) = [];
S_CT0.data(:,1:2) = [];
S_base.data(:,1:2) = [];

year = 2010:2050;
nyear = length(year);
list = variab(1:nyear:end); list=string(list);

nreal = length(2010:2019); ySIM = 2020; %YEAR OF STARTING SIMULATIONS
scenario = {'Baseline', 'CT_0', 'CT_{\beta}', 'CT_{\beta + \gamma}'}';
niter = size(S_CTbeta.data,1); l = length(scenario);
clc
clo = [0.1,0.1,0.1; 0.05, 0.4, 0.05; 0.64,0.08,0.21; 0.05, 0.45, 0.75;...%system of colors
       0.3,0.3,0.3;  0.80, 0.95, 0.80; 0.95, 0.75, 0.75;0.60, 0.90, 0.90 ];

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2: Scenario analysis: environmental-energy indicator
    % 2a : Carbon emissions
    % 2b : Clean electricity
    % 2c : TPES/GDP
    % 2d : HH energy consumption (COICOP)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 2a : Carbon emissions
% CO2_reduction_1990_level
x = find (list=='CO2_reduction_1990_level');
in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;

co2red = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    co2red(i,:) = nanmean(s);
    co2red(i+l,:) = 2*nanstd(s);
    
end
%aaa =100*aaa;
data = [96.923	94.500	90.692	84.846	81.789	83.250	82.289];

%500 simulazioni
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,co2red(i,:) + co2red(i+l,:),co2red(i,:) - co2red(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, co2red(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, co2red(i,:) + co2red(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, co2red(i,:) - co2red(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, co2red(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthEast','NumColumns',2);

h(2) = plot(year, co2red(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[20 100],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' CO_2 emissions (1990 = 100) '); set(gca,'DataAspectRatio',[1 80*1.5/40 1]);
axis([2010 2050 20 100]);
saveas(gcf,'Fig_2a_CO2red_std.jpg');
saveas(gcf,'Fig_2a_CO2red_std.eps');


%% 2b : Renew_share_Electricity

x = find (list=='Renew_share_Electricity');
in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;

ren_elect = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    ren_elect(i,:) = nanmean(s);
    ren_elect(i+l,:) = 2*nanstd(s);
    
end
ren_elect =100*ren_elect;

data = [0.235	0.274	0.313	0.334	0.335	0.340	0.341	0.349	0.356	0.364]*100;
%500 simulazioni
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,ren_elect(i,:) + ren_elect(i+l,:),ren_elect(i,:) - ren_elect(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, ren_elect(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, ren_elect(i,:) + ren_elect(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, ren_elect(i,:) - ren_elect(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, ren_elect(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthWest','NumColumns',2);

h(2) = plot(year, ren_elect(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[20 100],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Share renew. elect. (%) '); set(gca,'DataAspectRatio',[1 80*1.5/40 1]);
axis([2010 2050 20 100]);
saveas(gcf,'Fig_2b_Renew_share_Electricity_std.jpg');

clear A aaa i h 

%% 2c: TPES on GDP 

x = find (list=='En_Intensity');
in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn);
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn)  ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn);

enrg_intens = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    enrg_intens(i,:) = nanmean(s);
    enrg_intens(i+l,:) = 2*nanstd(s);
    
end
enrg_intens(1:4,:) = 100*enrg_intens(1:4,:)/67.12; % set 66.8611 to get 2015=100
data = [110	108	106	104	102.000	100.000 97]; %PNIEC data 2015=100
    %(see *https://www.eea.europa.eu/data-and-maps/indicators/total-primary-energy-intensity-4/assessment-1)
%emrg_intens(1:4,:) = 100*emrg_intens(1:4,:)./repmat(emrg_intens(1:4,6),1,41);%
%FIGURE
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,enrg_intens(i,:) + enrg_intens(i+l,:),enrg_intens(i,:) - enrg_intens(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, enrg_intens(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, enrg_intens(i,:) + enrg_intens(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, enrg_intens(i,:) - enrg_intens(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, enrg_intens(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthEast','NumColumns',2);

h(2) = plot(year, enrg_intens(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[60 120],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Energy Intensity (2015 = 100) '); set(gca,'DataAspectRatio',[1 60*1.45/40 1]);
axis([2010 2050 60 111]);
saveas(gcf,'Fig_2C_TPES_std.jpg');

clear A aaa i h data

%% 2d: HH energy consumption (COICOP)

x = find (list=='share_elect_HH_real');
in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;

HHenergy = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    HHenergy(i,:) = nanmean(s);
    HHenergy(i+l,:) = 2*nanstd(s);
    
end
HHenergy = HHenergy*100;

% NOTE: \footnote{\textcolor{blue}{Note that due to the lack of data we do 
    % not report the real values for Households' energy consumption. However 
    % from Eurostat we recover the initial value (2010) from the COICOP classification 
    % by subtracting from item 045 the amount related to rentals and water use, 
    % in order to get the expenditures on energy which was around 5.7\% of total 
    % expenditures, close to our simulation (5.54\%)}}

%500 simulazioni
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,HHenergy(i,:) + HHenergy(i+l,:),HHenergy(i,:) - HHenergy(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, HHenergy(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, HHenergy(i,:) + HHenergy(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, HHenergy(i,:) - HHenergy(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, HHenergy(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;

end
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthEast','NumColumns',2);

h(2) = plot(year, HHenergy(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[0 7],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Share Elect. COICOP (%) '); set(gca,'DataAspectRatio',[1 7*1.5/40 1]);
axis([2010 2050 0 7]);
saveas(gcf,'Fig_2d_share_Elect_COICOP_std.jpg');

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 3: Scenario analysis: economic and fiscal indicators
    % 3a : GDP real
    % 3b : Tot export
    % 3c : Output/GDP
    % 3d : def/gdp
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 3a : GDP real
x = find (list=='GDP_real');%
in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ; 

GDP = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    GDP(i,:) = nanmean(s);
    GDP(i+l,:) = 2*nanstd(s);
end
GDP = GDP/1e12; %in T€
data = [1604514500000	1637462900000	1613265000000	1604599100000	...
    1621827200000	1652085400000	1689824000000]/1e12; %in T€

%FIGURE
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,GDP(i,:) + GDP(i+l,:),GDP(i,:) - GDP(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, GDP(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, GDP(i,:) + GDP(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, GDP(i,:) - GDP(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, GDP(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthWest','NumColumns',2);

h(2) = plot(year, GDP(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[1.5 2.6],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' GDP real (Trillion €) '); set(gca,'DataAspectRatio',[1 1.1*1.5/40 1]);
axis([2010 2050 1.5 2.6]);
saveas(gcf,'Fig_3A_GDP_std.jpg');


%% 3b : Tot export

% clear A
% x = find (list=='Tot EXP real');%
% in = (x-1)*nyear+1; fn = x*nyear;
% A = zeros(niter, nyear, l);
%     A(:,:, 1) = S_basex.data(:, in:fn) ;
%     A(:,:, 2) = S_CT0x.data(:, in:fn) ;
%     A(:,:, 3) = S_CTbetax.data(:, in:fn) ;
%     A(:,:, 4) = S_CTbeta_gammax.data(:, in:fn) ;
% 
% expreal100 = zeros(l,nyear);
% for i = 1:l
%     s = A(:,:, i);
%     expreal100(i,:) = nanmean(s);
%     expreal100(i+l,:) = 2*nanstd(s);
% end
%real data from EUROSTAT as Exp/GDP rearranged in real terms given the CPI
 %because of lack of data (	26.9	28.4	28.6	29.1	29.7	29.3	30.7	31.4	31.6)
data  = [100,105.614215613383,107.558327137546,110.231665427509,113.384442379182,...
         112.432486988848,118.292007434944,121.178319702602,122.032386617100];

%FIGURE confid interval
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,expreal100(i,:) + expreal100(i+l,:),expreal100(i,:) - expreal100(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, expreal100(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, expreal100(i,:) + expreal100(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, expreal100(i,:) - expreal100(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, expreal100(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthWest','NumColumns',2);

h(2) = plot(year, expreal100(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[100 180],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Export (2010 = 100) '); set(gca,'DataAspectRatio',[1 80*1.5/40 1]);
axis([2010 2050 100 180]);
saveas(gcf,'Fig_3b_Export100_std.jpg');

%% 3c Economic efficiency gdp/output
% x = find (list=='Output_on_GDP_ratio'); in = (x-1)*nyear+1; fn = x*nyear;
% A = zeros(niter, nyear, l);
%     A(:,:, 1) = S_basex.data(:, in:fn) ;
%     A(:,:, 2) = S_CT0x.data(:, in:fn) ;
%     A(:,:, 3) = S_CTbetax.data(:, in:fn) ;
%     A(:,:, 4) = S_CTbeta_gammax.data(:, in:fn) ;
%     
% econ_effic = zeros(l,nyear);
% for i = 1:l
%     s = (A(:,:, i)).^(-1);
%     econ_effic(i,:) = nanmean(s);
%     econ_effic(i+l,:) = 2*nanstd(s);
% end
% econ_effic = econ_effic*100;

%DATA
% output = [4133.509798	4475.042942	4024.546706	4098.209012	4075.402119]
%gdp = [2136 2295 2087 2142 2162]
data = [4133.509798	4475.042942	4024.546706	4098.209012	4075.402119]./...
    [2136 2295 2087 2142 2162];
%inflation b = [1	1.00036	1.01165	1.01898	1.02695];	1.03223	1.0365	1.03812	1.03882	];
data = 100*data./b;

%FIGURE 
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,econ_effic(i,:) + econ_effic(i+l,:),econ_effic(i,:) - econ_effic(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, econ_effic(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, econ_effic(i,:) + econ_effic(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, econ_effic(i,:) - econ_effic(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, econ_effic(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthEast','NumColumns',2);

h(2) = plot(year, econ_effic(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[160 200],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Output to GDP (%) '); set(gca,'DataAspectRatio',[1 40*1.5/40 1]);
axis([2010 2050 160 200]);
saveas(gcf,'Fig_3c_EcoEff_std.jpg');

%% 3d: deficit / GDP

x = find (list=='Deficit_to_GDP'); in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn);
    A(:,:, 2) = S_CT0.data(:, in:fn);
    A(:,:, 3) = S_CTbeta.data(:, in:fn);
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ; 

deficitGDP = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    deficitGDP(i,:) = nanmean(s);
    deficitGDP(i+l,:) = 2*nanstd(s);
end
deficitGDP = deficitGDP*100;
%DATA
data = [0.042	0.037	0.029	0.029	0.030	0.026	0.025	0.024	0.021]*100;

%figure
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,deficitGDP(i,:) + deficitGDP(i+l,:),deficitGDP(i,:) - deficitGDP(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, deficitGDP(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, deficitGDP(i,:) + deficitGDP(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, deficitGDP(i,:) - deficitGDP(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, deficitGDP(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthEast','NumColumns',2);

h(2) = plot(year, deficitGDP(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[-2 6],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Deficit / GDP (%) '); set(gca,'DataAspectRatio',[1 8*1.5/40 1]);
axis([2010 2050 -2 6]);
saveas(gcf,'Fig_3d_DefGDP_std.jpg');

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 4: Scenario analysis: Inequality and labour market
    % 4a : Gini
    % 4b : Lab prod
    % 4c : Occupation
    % 4d : Unempl
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 4a: GINI

x = find (list=='GINI_Coeff'); in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;
%  GINI AVG
gini = zeros(l,nyear);
for i = 1:l
    s = A(:,:, i);
    gini(i,:) = nanmean(s);
    gini(i+l,:) = 2*nanstd(s);
end
gini = gini*100;
data = [0.347	0.351	0.352	0.349	0.347	0.350]*100; %REAL DATA GINI 2010-2015

%FIGURE confid interval
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,gini(i,:) + gini(i+l,:),gini(i,:) - gini(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, gini(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, gini(i,:) + gini(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, gini(i,:) - gini(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, gini(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthWest','NumColumns',2);

h(2) = plot(year, gini(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[33 39],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Gini Coefficient (%) '); set(gca,'DataAspectRatio',[1 6*1.5/40 1]);
axis([2010 2050 33 39]);
saveas(gcf,'Fig_4a_Gini_std.jpg');

%% 4b: labor productivity

x = find (list=='avg_lambda_index');
in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;

lambda = zeros(l,nyear);
for i = 1:l
    s = A(:,:, i);
    lambda(i,:) = nanmean(s);
    lambda(i+l,:) = 2*nanstd(s);
end
%lambda = lambda*100;
data = [100.000	100.500	100.200	101.100	101.500	101.500	101.100	101.700	101.500];

%FIGURE confid interval
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,lambda(i,:) + lambda(i+l,:),lambda(i,:) - lambda(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, lambda(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, lambda(i,:) + lambda(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, lambda(i,:) - lambda(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, lambda(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthWest','NumColumns',2);

h(2) = plot(year, lambda(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[99 130],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Labour prod. index (2010 = 100)'); set(gca,'DataAspectRatio',[1 31*1.5/40 1]);
axis([2010 2050 99 130]);
saveas(gcf,'Fig_4b_lambda_std.jpg')

%% 4c: Tot employment

x = find (list=='Tot_N_E'); in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;

EMPL = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    EMPL(i,:) = nanmean(s);
    EMPL(i+l,:) = 2*nanstd(s);
    
end
EMPL =1e-6*EMPL;
data = [22.341, 22.416, 22.282, 21.86275, 21.9225]; %from ISTAT as average over trimestal data

%500 simulazioni
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,EMPL(i,:) + EMPL(i+l,:),EMPL(i,:) - EMPL(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, EMPL(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, EMPL(i,:) + EMPL(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, EMPL(i,:) - EMPL(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, EMPL(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data'); set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthWest','NumColumns',2);

h(2) = plot(year, EMPL(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[20 28],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Employment (M pers.) '); set(gca,'DataAspectRatio',[1 8*1.5/40 1]);
axis([2010 2050 20 28]);
 saveas(gcf,'Fig_4c_Empl_std.jpg');

%% 4d: Tot U rate

x = find (list=='Tot_U_rate'); in = (x-1)*nyear+1; fn = x*nyear;
A = zeros(niter, nyear, l);
    A(:,:, 1) = S_base.data(:, in:fn) ;
    A(:,:, 2) = S_CT0.data(:, in:fn) ;
    A(:,:, 3) = S_CTbeta.data(:, in:fn) ;
    A(:,:, 4) = S_CTbeta_gamma.data(:, in:fn) ;

Urate = zeros(l*2,nyear);
for i = 1:l
    s = A(:,:, i);
    Urate(i,:) = nanmean(s);
    Urate(i+l,:) = 2*nanstd(s);
    
end
Urate =100*Urate;
data = [0.084	0.084	0.107	0.121	0.127	0.119	0.117	0.112	0.106]*100;

%figure
figure('units', 'normalized', 'outerposition', [-0.0046875 0.09375 1.009375 0.91375])%[0 0.05 1 0.95]
axes('Position',[0.0837239583333332 0.100160256410257 0.885416666666665 0.871794871794872]);
h = zeros(1,10);
for i = 1:l
h(i) = jbfill(year,Urate(i,:) + Urate(i+l,:),Urate(i,:) - Urate(i+l,:),clo(i+l,:),clo(i+l,:),1,0.3); hold on;
set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, Urate(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, Urate(i,:) + Urate(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
 h(i) = plot(year, Urate(i,:) - Urate(i+l,:), 'LineStyle',':','LineWidth', 0.2,'Color',clo(i,:)); hold on; grid on;
 set(get(get(h(i),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
end
set(gca,'FontSize',20,'FontWeight','bold'); 

for i = 1:l
 plot(year, Urate(i,:), 'LineStyle','-','LineWidth', 2,'Color',clo(i,:)); hold on; grid on;
 box on;
end
plot(year(1:length(data)), data, 'Marker','*','LineWidth',2,'LineStyle','-.',...
    'Color',[0.5 0.2 0.5]); hold on; %DATA REAL
legend('\color{black} Baseline','\color[rgb]{0.05, 0.4, 0.05} CT_0',...
       '\color[rgb]{0.64,0.08,0.21} CT_{\beta}','\color[rgb]{0.05, 0.45, 0.75} CT_{\beta + \gamma}',...
       '\color[rgb]{0.5,0.2,0.5} data');
set(legend,'Orientation','horizontal','FontSize',15,'Location','NorthEast','NumColumns',2);

h(2) = plot(year, Urate(1,:), 'LineStyle','-','LineWidth', 2,'Color',clo(1,:)); hold on; grid on;
set(get(get(h(2),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
h(3) = line([ySIM ySIM],[6 14],'LineStyle','-.','LineWidth', 1.5,'Color',[0.8,0.7,0.6]);
set(get(get(h(3),'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
ylabel(' Unemployment rate (%) '); set(gca,'DataAspectRatio',[1 8*1.5/40 1]);
axis([2010 2050 6 14]);
saveas(gcf,'Fig_4d_Urate_std.jpg');

clear a aa ans b fn xn i in h r s
