% ===========================================================================================
% Replication: Barigozzi, Lissona & Tonni (2026) - "Large datasets for the Euro Area and its
% member countries and the dynamic effects of the common monetary policy"
%
% Correlation analysis linking the country-level monetary policy IRFs (loaded from
% outputs/ccvar_monthly_iv/Results.mat, produced by ccvar_monthly_iv.m) to cross-country
% structural characteristics. Produces Table 6 of the paper.
% -------------------------------------------------------------------------------------------
% Author: Lorenzo Tonni (lorenzo.tonni@unimi.it)
% Last update: 09/Jun/2026
% Version: MATLAB 2026a
% -------------------------------------------------------------------------------------------

% Settings and data import
clear all; clear session; close all; clc
warning off all

% ===============|| Set Output Folder ||============= %
if isfolder(join([pwd,'\outputs'],'')); Folder_output = join([pwd,'\outputs'],''); % windows
else; Folder_output = join([pwd,'/outputs'],''); % mac
end
addpath(join([pwd,'/outputs/ccvar_monthly_iv'],''));

% ====================|| Set path ||================= %
if isfolder(join([pwd,'\institutional_data'],'')); addpath(join([pwd,'\institutional_data'],'')); % windows
else; addpath(join([pwd,'/institutional_data'],''));
end

if isfolder(join([pwd,'\Functions'],'')); addpath(join([pwd,'\Functions'],'')); % windows
else; addpath(join([pwd,'/Functions'],''));
end

if isfolder(join([pwd,'\Results'],'')); addpath(join([pwd,'\Results'],'')); % windows
else; addpath(join([pwd,'/Results'],''));
end

% ===============|| Set Figures Folder ||============= %
if isfolder(join([pwd,'\Figures'],'')); folder_fig = join([pwd,'\Figures'],''); % windows
else; folder_fig = join([pwd,'/Figures'],''); % mac
end

%% Import Covariates

data_temp = readtable('dataset.xlsx', 'Sheet', 'Foglio1', 'VariableNamingRule','preserve');
covariates.data = table2array(data_temp(:,2:end));
covariates.countries = data_temp.Properties.VariableNames; covariates.countries(:,1) = [];
covariates.id=table2array(data_temp(:,1));

%% Import IRFs

results_file = fullfile(Folder_output, 'ccvar_monthly_iv', 'Results.mat');
loaded_results = load(results_file); TR = 3;
if isfield(loaded_results, 'results')
    results = loaded_results.results;
elseif isfield(loaded_results, 'IRF')
    results = loaded_results.IRF;
else
    error('Results.mat must contain either a results or IRF structure.')
end

irf_temp = squeeze(results.irf_nationals_cum); irf_temp(:,20:29) = squeeze(results.irf_nationals(:,20:29));
irf_IPMN = min(irf_temp(:,1:9), [], 1); irf_IPMN = irf_IPMN';
irf_HICP = min(irf_temp(:,10:19), [], 1); irf_HICP = irf_HICP';
irf_LTIRT = max(irf_temp(:,20:29), [], 1); irf_LTIRT = irf_LTIRT';
irf_SHIX = min(irf_temp(:,30:39), [], 1); irf_SHIX = irf_SHIX';
irf_UNETOT = max(irf_temp(:,40:49), [], 1); irf_UNETOT = irf_UNETOT';

countries_default = {'AT','BE','DE','EL','ES','FR','IE','IT','NL','PT'};
countries_ipmn = {'AT','BE','DE','EL','ES','FR','IT','NL','PT'};

var_IPMN = explained_variance(results.EXP_VAR, 'IPMN', countries_ipmn);
var_HICP = explained_variance(results.EXP_VAR, 'HICPOV', countries_default);
var_LTIRT = explained_variance(results.EXP_VAR, 'LTIRT', countries_default);
var_SHIX = explained_variance(results.EXP_VAR, 'SHIX', countries_default);
var_UNETOT = explained_variance(results.EXP_VAR, 'UNETOT', countries_default);

%% Select Covariates
variables = { ...
'employment_protection_index', ...
'wage_adj_more_1_once_year_AVG_A', ...
'wage_adj_once_year_AVG_A', ...
'wage_adj_less_once_year_AVG_A', ...
'price_change_freq_baseline', ...
'leverage_to_gdp_nf_comp', ...
'homeownership_rate_tot', ...
'homeownership_rate_mortgage', ...
'homeownership_rate_no_mortgage', ...
'2022_FRM_share', ...
'ltv', ...
'htm', ...
'wealthy_htm', ...
'saving_rate_percentage_gdp' ...
};

table6_labels = { ...
'Employment Protection', ...
'Wage Adj. Frequency: more than once a year', ...
'Wage Adj. Frequency: once a year', ...
'Wage Adj. Frequency: less than once a year', ...
'Price Flexibility', ...
'NFCLeverage', ...
'Homeownership', ...
'Homeownership with mortgage', ...
'Homeownership without mortgage', ...
'Fixed-Rate Mortgages', ...
'Loan to Value', ...
'Share of HtM', ...
'Share of WHtM', ...
'Saving Rate' ...
};

[found_variables,posizioni] = ismember(variables, covariates.id);
if any(~found_variables)
    error('The following Table 6 variables are missing from dataset.xlsx: %s', strjoin(variables(~found_variables), ', '))
end

regressors = covariates.data(posizioni, :)';
ireland_idx = find(ismember(covariates.countries, {'IR','IE'}), 1);
if isempty(ireland_idx)
    error('Ireland must be identified in dataset.xlsx as either IR or IE.')
end

regressors_IPMN = regressors; regressors_IPMN(ireland_idx,:) = [];
regressors_IPMN = [regressors_IPMN, var_IPMN];
regressors_HICP = [regressors, var_HICP];
regressors_LTIRT = [regressors, var_LTIRT];
regressors_SHIX = [regressors, var_SHIX];
regressors_UNETOT = [regressors,var_UNETOT];

labels_IPMN = {'AT','BE','DE','EL','ES','FR','IT','NL','PT'};
labels_default = {'AT','BE','DE','EL','ES','FR','IE','IT','NL','PT'};

%% Correlation Analysis (IPMN)

X = regressors_IPMN;
y = irf_IPMN;

variables_f = [variables, {'var_IPMN'}];
IPMN = correlation_table(variables_f, X, y);

%% Correlation Analysis (HICP)

X = regressors_HICP;
y = irf_HICP;

variables_f = [variables, {'var_HICP'}];
HICP = correlation_table(variables_f, X, y);

%% Correlation Analysis (LTIRT)

X = regressors_LTIRT;
y = irf_LTIRT;

variables_f = [variables, {'var_LTIRT'}];
LTIRT = correlation_table(variables_f, X, y);

%% Correlation Analysis (SHIX)

X = regressors_SHIX;
y = irf_SHIX;

variables_f = [variables, {'var_SHIX'}];
SHIX = correlation_table(variables_f, X, y);

%% Correlation Analysis (UNETOT)
X = regressors_UNETOT;
y = irf_UNETOT;

variables_f = [variables, {'var_UNETOT'}];
UNETOT = correlation_table(variables_f, X, y);

table6_labels = [table6_labels, {'Explained Variance'}];
print_table6(table6_labels, IPMN, HICP, LTIRT, SHIX, UNETOT);

function out = correlation_table(variable_names, X, y)

nVars = size(X,2);
corrVals = nan(nVars,1);
pVals = nan(nVars,1);
stdErrs = nan(nVars,1);

for i = 1:nVars
    X_temp = X(:,i);
    y_temp = y;
    valid_obs = ~isnan(X_temp) & ~isnan(y_temp);
    X_temp = X_temp(valid_obs,:);
    y_temp = y_temp(valid_obs,:);
    n = length(X_temp);

    if n > 2 && std(X_temp) > 0 && std(y_temp) > 0
        [r, p] = corr(X_temp, y_temp, 'Type', 'Pearson');
        corrVals(i) = r;
        pVals(i) = p;
        stdErrs(i) = sqrt((1 - r^2) / (n - 2));
    end
end

out = table(variable_names', corrVals, pVals, stdErrs, ...
    'VariableNames', {'Variable', 'Correlation', 'pValue', 'stdErrs'});

end

function exp_var = explained_variance(EXP_VAR, series_code, countries)

row_names = EXP_VAR.Properties.RowNames;
point_row = find(strcmp(row_names, 'ExpVar'), 1);
if isempty(point_row)
    error('EXP_VAR must contain a row named ExpVar.')
end

exp_var = nan(numel(countries),1);
for i = 1:numel(countries)
    var_name = [series_code, '_', countries{i}];
    var_idx = find(strcmp(EXP_VAR.Properties.VariableNames, var_name), 1);
    if isempty(var_idx)
        error('EXP_VAR does not contain %s.', var_name)
    end
    exp_var(i) = EXP_VAR{point_row,var_idx};
end

end

function print_table6(row_labels, IPMN, HICP, LTIRT, SHIX, UNETOT)

tables = {IPMN, HICP, LTIRT, SHIX, UNETOT};
headers = {'IPMN', 'HICPOV', 'LTIRT', 'SHIX', 'UNETOT'};

disp(' ')
disp('  ==================================================================================')
disp('  = Table 6: Cross-country characteristics and IRFs dynamics: correlation analysis =')
disp('  ==================================================================================')
disp(' ')

fprintf('  %-4s  %-45s', 'ID', 'Variable');
for j = 1:numel(headers)
    fprintf('%12s', headers{j});
end
fprintf('\n');

for i = 1:numel(row_labels)
    fprintf('  %-4d  %-45s', i, row_labels{i});
    for j = 1:numel(tables)
        fprintf('%12s', formatted_correlation(tables{j}.Correlation(i), tables{j}.pValue(i)));
    end
    fprintf('\n');

    fprintf('  %-4s  %-45s', '', '');
    for j = 1:numel(tables)
        fprintf('%12s', formatted_standard_error(tables{j}.stdErrs(i)));
    end
    fprintf('\n');
end

fprintf('  Notes: standard errors in parentheses; asterisks denote statistical significance at 99%% (***), 95%% (**), and 90%% (*).\n');

end

function txt = formatted_correlation(value, pvalue)

if isnan(value)
    txt = '-';
else
    txt = sprintf('%0.2f%s', value, significance_stars(pvalue));
end

end

function txt = formatted_standard_error(value)

if isnan(value)
    txt = '(-)';
else
    txt = sprintf('(%0.2f)', value);
end

end

function stars = significance_stars(pvalue)

if isnan(pvalue)
    stars = '';
elseif pvalue < 0.01
    stars = '***';
elseif pvalue < 0.05
    stars = '**';
elseif pvalue < 0.10
    stars = '*';
else
    stars = '';
end

end
