%% Settings and data import
% ==========================================================================
% Reproduces the three main-text figures of Barigozzi, Lissona & Tonni (2026)
% from the saved CC-VAR/IV results:
%   Figure 1  - EA-wide IRFs (2x3)
%   Figure 2  - Country-level IRFs, countries overlaid per variable (3+2)
%   Figure 3  - Country-minus-EA IRF differences, split into one subfigure (3a-3e)
%               per variable for readability
% --------------------------------------------------------------------------
% SEE ALSO: ccvar_monthly_iv.m, correlation_analysis_IRFs.m
% --------------------------------------------------------------------------
% Author: Claudio Lissona (claudio.lissona@gmail.com)
% Last update: 09/Jun/2026
% Version: MATLAB R2026a
% Required Toolboxes: /
% ==========================================================================
clear; clc; close all;                                                   % reset

root = fileparts(mfilename('fullpath'));                                 % script dir
if isempty(root); root = pwd; end                                        % cell-mode fallback

folder_fig = fullfile(root,'Figures');                                   % output folder
if ~isfolder(folder_fig); mkdir(folder_fig); end                         % create if missing

load(fullfile(root,'outputs','ccvar_monthly_iv','Results.mat'));         % load results

H = 36;                                                                  % IRF horizon
hor = (0:H-1)';                                                          % horizon axis

has_bands = ~all(isnan(results.lower_ccvar(:)));                         % bands present
if ~has_bands
    warning(['Band fields are NaN (conf=0): plotting point estimates only. ' ...
             'Re-run ccvar_monthly_iv.m with conf=1 for bands.']);
end

% Country colours (RGB/255)
colmap = [204 37 41; 57 106 177; 62 150 81; 218 124 48; 107 76 154; ...
          44 160 178; 146 36 40; 148 139 61; 216 82 24; 80 122 156]/255;
ctry_all = {'AT','BE','DE','EL','ES','FR','IE','IT','NL','PT'};

% Variable blocks: column range into the rotation dimension and plot options
blk(1) = struct('title','IP: Manufacturing','start', 1,'ctry',{ctry_all([1:6 8:10])},'scale',100,'lev',false);
blk(2) = struct('title','HICP: Overall',    'start',10,'ctry',{ctry_all},'scale',100,'lev',false);
blk(3) = struct('title','10-Year Rate',     'start',20,'ctry',{ctry_all},'scale',  1,'lev',true );
blk(4) = struct('title','Stock Price Index','start',30,'ctry',{ctry_all},'scale',100,'lev',false);
blk(5) = struct('title','Unemployment Rate','start',40,'ctry',{ctry_all},'scale',  1,'lev',false);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Figure 1 - EA-wide IRFs %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ea_title = {'2-Year Rate','IP: Manufacturing','HICP: Overall', ...
            '10-Year Rate','Stock Price Index','Unemployment Rate'};
ea_scale = [1 100 100 1 100 1];                                          % scale to % pts
ea_lev   = [true false false true false false];                          % level vs cumulated
ea_ord   = [1 2 3 5 4 6];                                                % panel order

fig1 = figure('Position',[100 100 1080 560],'Color','w');
t1 = tiledlayout(fig1,2,3,'TileSpacing','compact','Padding','compact');

for nn = 1:6
    col = ea_ord(nn);                                                    % data column
    ax = nexttile(t1);
    if ea_lev(col)
        pt = ea_scale(col)*results.irf_ccvar(1:H,col,1);                 % point estimate
        lo = ea_scale(col)*median(results.lower_ccvar(1:H,col,:),3);     % 68% lower
        hi = ea_scale(col)*median(results.upper_ccvar(1:H,col,:),3);     % 68% upper
    else
        pt = ea_scale(col)*results.irf_ccvar_cum(1:H,col,1);             % point estimate
        lo = ea_scale(col)*median(results.lower_ccvar_cum(1:H,col,:),3); % 68% lower
        hi = ea_scale(col)*median(results.upper_ccvar_cum(1:H,col,:),3); % 68% upper
    end
    shadeband(hor,lo,hi,[.6 .6 .6]);                                     % 68% shade
    hold(ax,'on');
    plot(ax,hor,pt,'k-','LineWidth',2.5);                                % point estimate
    yline(ax,0);                                                         % zero line
    title(ax,ea_title{col},'FontWeight','normal');
    xlim(ax,[0 H-1]); xticks(ax,0:6:H-1);                                % horizon ticks
    setylim(ax,[pt;lo;hi],[0 0]);                                        % margins
    grid(ax,'on'); box(ax,'on');                                         % grid + box
    hold(ax,'off');
end

title(t1,'Figure 1: EA IRFs');
exportgraphics(fig1,fullfile(folder_fig,'Figure1_EA.pdf'),'ContentType','vector');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2. Figure 2 - Country-level IRFs %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ord  = [1 2 4 3 5];                                                      % panel order
ypad = {[.1 0],[.1 0],[.25 0],[.05 0],[.8 0]};                           % y-pad per variable

fig2 = figure('Position',[100 100 980 620],'Color','w');
t2 = tiledlayout(fig2,2,6,'TileSpacing','compact','Padding','compact');
tile_pos = {1,3,5,8,10};                                                 % centred 3+2 tiles
leg_h = gobjects(1,numel(ctry_all)); leg_txt = cell(1,numel(ctry_all));

for nn = 1:5
    v = ord(nn);                                                         % block index
    cols = blk(v).start:(blk(v).start+numel(blk(v).ctry)-1);             % result columns
    nC = numel(cols);                                                    % country count
    if blk(v).lev
        pt = blk(v).scale*results.MiddleH_Big(1:H,cols);                 % point estimates
        lo = blk(v).scale*results.lower_nationals(1:H,cols);             % 68% lower
        hi = blk(v).scale*results.upper_nationals(1:H,cols);             % 68% upper
    else
        pt = blk(v).scale*results.MiddleH_Big_cum(1:H,cols);             % point estimates
        lo = blk(v).scale*results.lower_nationals_cum(1:H,cols);         % 68% lower
        hi = blk(v).scale*results.upper_nationals_cum(1:H,cols);         % 68% upper
    end
    sig = sign(lo)==sign(hi);                                            % significant if band excludes 0
    cpos = cellfun(@(c) find(strcmp(ctry_all,c)), blk(v).ctry);          % colour index
    cmap = colmap(cpos,:);                                               % fixed colour per country
    ax = nexttile(t2,tile_pos{nn},[1 2]);
    hold(ax,'on');
    for cc = 1:nC
        psig = pt(:,cc); psig(~sig(:,cc)) = NaN;                         % significant horizons
        hs = plot(ax,hor,psig,'-','Color',cmap(cc,:),'LineWidth',2.5);   % solid: significant
        plot(ax,hor,pt(:,cc),':','Color',cmap(cc,:),'LineWidth',2.5);    % dotted: full
        if v==2; leg_h(cc) = hs; leg_txt{cc} = blk(v).ctry{cc}; end      % legend from HICP
    end
    yline(ax,0);                                                         % zero line
    title(ax,blk(v).title,'FontWeight','normal');
    xlim(ax,[0 H-1]); xticks(ax,0:6:H-1);                                % horizon ticks
    setylim(ax,pt(:),ypad{v});                                           % margins from point est. (no bands drawn)
    grid(ax,'on'); box(ax,'on');                                         % grid + box
    hold(ax,'off');
end

lgd = legend(leg_h,leg_txt,'Orientation','horizontal','NumColumns',10, ...
             'Interpreter','latex','IconColumnWidth',15);
lgd.Layout.Tile = 'south';
title(t2,'Figure 2: Country-level IRFs');
exportgraphics(fig2,fullfile(folder_fig,'Figure2_countries.pdf'),'ContentType','vector');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3.   Figure 3 - Country - EA IRFs (3a-3e)   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ord  = [1 2 3 4 5];                                                      % block order
core = {'AT','BE','DE','FR','NL'};                                       % core (upper row)
peri = {'EL','ES','IE','IT','PT'};                                       % periphery (lower row)
grp  = [core, peri];                                                     % 10 cells per block
lett = {'a','b','c','d','e'};                                            % subfigure letters
tag  = {'IP','HICP','LTIRT','SHIX','UNETOT'};                            % filename tags

for nn = 1:numel(ord)
    v = ord(nn);                                                         % block index
    fig3 = figure('Position',[80 80 1150 500],'Color','w');
    t3 = tiledlayout(fig3,2,5,'TileSpacing','compact','Padding','compact');
    for cc = 1:numel(grp)
        ax = nexttile(t3);
        loc = find(strcmp(blk(v).ctry,grp{cc}),1);                       % country col in block
        if ~isempty(loc)
            col = blk(v).start+loc-1;                                    % result column
            if blk(v).lev
                pt = blk(v).scale*results.difference(1:H,col);           % difference
                lo = blk(v).scale*results.lower_difference(1:H,col);     % 68% lower
                hi = blk(v).scale*results.upper_difference(1:H,col);     % 68% upper
            else
                pt = blk(v).scale*results.difference_cum(1:H,col);       % difference
                lo = blk(v).scale*results.lower_difference_cum(1:H,col); % 68% lower
                hi = blk(v).scale*results.upper_difference_cum(1:H,col); % 68% upper
            end
            shadeband(hor,lo,hi,[.6 .6 .6]);                             % 68% shade
            hold(ax,'on');
            plot(ax,hor,pt,'k-','LineWidth',1.5);                        % difference
            yline(ax,0);                                                 % zero line
            grid(ax,'on'); box(ax,'on');                                 % grid + box
            hold(ax,'off');
            xlim(ax,[0 H-1]); xticks(ax,0:6:H-1);                        % horizon ticks
            setylim(ax,[pt;lo;hi],[0 0]);                                % per-panel margins
        else
            set(ax,'XTick',[],'YTick',[]); xlim(ax,[0 H-1]);             % blank tile
        end
        title(ax,grp{cc},'FontWeight','normal','FontSize',9);            % country code
        if cc<=5; set(ax,'XTickLabel',[]); end                           % x labels on periphery row
    end
    title(t3,sprintf('Figure 3%s - %s (country minus EA)',lett{nn},blk(v).title));
    exportgraphics(fig3,fullfile(folder_fig,sprintf('Figure3%s_%s.pdf',lett{nn},tag{nn})),'ContentType','vector');
end

fprintf('Figures saved to %s\n',folder_fig);                             % done
close all;                                                               % close all figures when done

%%%%%%%%%%%%%%%%%%%%%%%%%
%% AUXILIARY FUNCTIONS %%
%%%%%%%%%%%%%%%%%%%%%%%%%

function shadeband(x,lo,hi,col)
% ==========================================================================
% Shade the [lo, hi] band over x on the current axes; none if lo/hi all NaN.
% ==========================================================================
if all(isnan(lo)) || all(isnan(hi)); return; end                     % nothing to draw
   fill([x; flipud(x)],[lo; flipud(hi)],col,'FaceAlpha',0.5,'EdgeColor','none');
end

function setylim(ax,M,yscale)
% ============================
%     Set y-axis limits
% ============================
vals = M(isfinite(M));                                               % finite only
if isempty(vals); return; end
limy = [min(vals) max(vals)];                                        % data range
if limy(1)==limy(2); return; end                                     % flat: skip
if floor(log(abs(limy(1)))/log(10))==3
    pc = 0.0008;                                                     % large scale
elseif abs(floor(log(abs(limy(1)))/log(10)))<=1
    pc = 0.08;                                                       % small scale
else
    pc = 0.008;                                                      % mid scale
end
if sign(limy(1))==1 && sign(limy(2))==1
    limy = limy + [-pc*limy(1) pc*limy(2)];                          % both positive
elseif sign(limy(1))~=1 && sign(limy(2))~=1
    limy = limy + [pc*limy(1) -pc*limy(2)];                          % both negative
else
    limy = limy + limy.*pc;                                          % opposite signs
end
limy = limy + yscale.*limy;                                          % extra padding
if limy(2)>limy(1); ylim(ax,limy); end                               % apply if valid
end


