
function [results] = loop_ccvar(conf, DATA, HEADER, DATES, ID, policy_instrument, mpi)

rr = 6; % Number of common factors

[header_national_rotations, nationals_rotations]=nationals_rotation_monthly(HEADER.header_AT,HEADER.header_BE,HEADER.header_DE,HEADER.header_EL,HEADER.header_ES,HEADER.header_FR,HEADER.header_IE,HEADER.header_IT,HEADER.header_NL,HEADER.header_PT,          DATA.data_AT,DATA.data_BE,DATA.data_DE,DATA.data_EL,DATA.data_ES,DATA.data_FR,DATA.data_IE,DATA.data_IT,DATA.data_NL,DATA.data_PT);

% Merge EA + nationals
data_rotations=[DATA.data, nationals_rotations];
header_rotations=[HEADER.header, header_national_rotations];
n=size(DATA.data,2);

% % Estimate r 
kmax = 20;
cmax = 3;
[data_rotations_std, ~] =standardize(data_rotations);
[rhat1_abc, ~] = ABC_crit(data_rotations, kmax, cmax, 0);          
[results]=NbFactors(data_rotations_std,kmax); 
r_bai_ng=results.khat_IC(2);
[r_ER,r_GR] = AH(data_rotations_std,kmax,0);
out=onatski(data_rotations_std,kmax);

disp(' ')
disp('  ========================================')
disp('  = Table 3: Estimated Number of Factors =')
disp('  ========================================')
disp(' ')

disp(sprintf('  According to Bai and Ng (2002):  %1.0f common factor(s)',r_bai_ng)), disp(' ')

disp(sprintf('  According to Alessi et al. (2010):  %1.0f common factor(s)',rhat1_abc)), disp(' ')

disp(sprintf('  According to Onatski (2010):  %1.0f common factor(s)',out)), disp(' ')

disp(sprintf('  According to Ahn and Horenstein (2013):  %1.0f common factor(s)',r_ER)), disp(' ')

% PC and Common Components 
[pc_rotations, common_components_rotations, finaldata, loadings,~,EXP_VAR] = pca_monthly(data_rotations,DATES.date_EU, rr, header_rotations);
print_explained_variance_table(EXP_VAR);

results.factors = pc_rotations;
results.loadings = loadings;
% Store variables' names
[id] = save_id(header_rotations); 

N=size(common_components_rotations,2);
chi =[common_components_rotations, pc_rotations];

% =================%
%   Select data    %
% =================%

psi = [id.IRT3M_EACC, id.IPMN_EA, id.HICPOV_EA, id.LTIRT_EACC, id.SHIX_EA, id.UNETOT_EA];  %   
                          
hor=36;  %  Select IRF horizon                                      

	% Last pre-COVID observation
	T0 = find(DATES.date_EU == datetime(2020,2,1), 1); 
	if isempty(T0) && DATES.date_EU(end) >= datetime(2020,2,1)
	    error('The date vector does not contain February 2020, needed for LP volatility correction.');
	end

% Figure details
VARnames={'2-Year Interest Rate';'IP'; 'HICP'; '10-Year Interest Rate'; 'SHIX' ; 'Unetot'};

% national variables in common components         
	rotation_variables=chi(:,(n+1):N);                
	
	n_rotations=size(rotation_variables,2);
	valid_iv = all(~isnan(mpi),2);
	date_iv_nonan = DATES.date_iv(valid_iv,:);
	mpi_nonan = mpi(valid_iv,:);


results.irf_ccvar=zeros(hor,size(VARnames,1),n_rotations);
results.upper_ccvar=zeros(hor,size(VARnames,1),n_rotations);
results.lower_ccvar=zeros(hor,size(VARnames,1),n_rotations);

results.irf_ccvar_cum=zeros(hor,size(VARnames,1),n_rotations);
results.upper_ccvar_cum=zeros(hor,size(VARnames,1),n_rotations);
results.lower_ccvar_cum=zeros(hor,size(VARnames,1),n_rotations);

results.irf_nationals=zeros(hor,1,n_rotations);
results.upper_nationals=zeros(hor,1,n_rotations);
results.lower_nationals=zeros(hor,1,n_rotations);

results.irf_nationals_cum=zeros(hor,1,n_rotations);
results.upper_nationals_cum=zeros(hor,1,n_rotations);
results.lower_nationals_cum=zeros(hor,1,n_rotations);

% =========== Estimate CC-VAR Covid Volatility ============= %

pp = 8; % number of lags

t=size(DATA.data,1);

if ~isempty(T0) && t > T0
st_f_all_ccvar=nan(t-pp,n_rotations);

for j=1:n_rotations

% dataset for VAR estimates
vardata_star = chi(:, psi); % fixed EA variables
vardata = [vardata_star, rotation_variables(:,j)]; % add national LTIRT and national variable of interest


% Use the policy rate without the common component
vardata(:,1)=policy_instrument; 

% IRFs
	c=1; % constant
% IRFs
[st_f, theta_f]=volatility_LP(vardata,pp,c,T0,2);
st_f_all_ccvar(:,j)=st_f;
end

st_f_ccvar=median(st_f_all_ccvar,2);

else

st_f_ccvar= ones(t-pp,1);

end

% =============== Estimate CC-VAR ============== %


F_stat = nan(n_rotations,2);
rhoU_all = cell(n_rotations,1);

% Progress bar over rotations (each rotation runs the full bootstrap, so this
% tracks the slow conf==1 path; skipped when no bands are estimated)
if conf == 1                                                            % only when bootstrapping bands
    wb = waitbar(0, 'Bootstrap: 0%');                                  % create progress bar
    wbCleanup = onCleanup(@() delete(wb));                             % safeguard: close bar on error/exit
end

for j=1:n_rotations
% dataset for VAR estimates
vardata_star = chi(:, psi);

vardata = [vardata_star, rotation_variables(:,j)];

% policy rate not in cc form
vardata(:,1)=policy_instrument; 

% IRFs
n=size(vardata,2); % # of variables
t=size(vardata,1); % # of periods
c=1; % constant

% IRFs
	[MiddleH_cum, MiddleH, difference_national_ea_temp, difference_national_ea_cum_temp, intervals, F, rhoU]=irf_iv(st_f_ccvar, vardata, n,pp,c,hor,mpi_nonan,DATES.date_EU,date_iv_nonan, conf);
	F_stat(j,1) = F.F_statistic;
	F_stat(j,2) = F.F_eff;
	rhoU_all{j} = rhoU;

difference_national_ea(:,:,j)=difference_national_ea_temp;
difference_national_ea_cum(:,:,j)=difference_national_ea_cum_temp;

% EA - National differences
if j <= 9
   index = 2;
elseif j >= 10 && j <= 19
   index = 3;
elseif j >= 20 && j <= 29
   index = 4;
elseif j >= 30 && j <= 39
   index = 5;
elseif j >= 40 && j <= 49
   index = 6;
end

% National Differences with EA

results.difference(:,j)=difference_national_ea(:,index,j);
results.upper_difference(:,j)=intervals.HighH_diff(:,index);
results.lower_difference(:,j)=intervals.LowH_diff(:,index);

%-------------------------------------------------------------
% National Differences with EA - Cumulate

results.difference_cum(:,j)=difference_national_ea_cum(:,index,j);
results.upper_difference_cum(:,j)=intervals.HighH_diff_cum(:,index);
results.lower_difference_cum(:,j)=intervals.LowH_diff_cum(:,index);

%-------------------------------------------------------------
% EA fixed variables

results.irf_ccvar(:,:,j)=MiddleH(:,1:size(VARnames,1)); 
results.upper_ccvar(:,:,j)=intervals.HighH(:,1:size(VARnames,1));
results.lower_ccvar(:,:,j)=intervals.LowH(:,1:size(VARnames,1));

%-------------------------------------------------------------
% EA fixed variables (Cumulate)

results.irf_ccvar_cum(:,:,j)=MiddleH_cum(:,1:size(VARnames,1)); 
results.upper_ccvar_cum(:,:,j)=intervals.HighH_cum(:,1:size(VARnames,1));
results.lower_ccvar_cum(:,:,j)=intervals.LowH_cum(:,1:size(VARnames,1));

%-------------------------------------------------------------
% National Variables Absolute

results.irf_nationals(:,:,j)=MiddleH(:,n);
results.upper_nationals(:,:,j)=intervals.HighH(:,n);
results.lower_nationals(:,:,j)=intervals.LowH(:,n);

%-------------------------------------------------------------
% National Varables Absolute (Cumulate)

results.irf_nationals_cum(:,:,j)=MiddleH_cum(:,n);
results.upper_nationals_cum(:,:,j)=intervals.HighH_cum(:,n);
results.lower_nationals_cum(:,:,j)=intervals.LowH_cum(:,n);

% Refresh progress bar after each completed rotation
if conf == 1                                                            % only when bootstrapping bands
    frac = j/n_rotations;                                               % fraction of rotations done
    waitbar(frac, wb, sprintf('Bootstrap: %.0f%%', 100*frac));         % update bar with percentage only
end

end

if conf == 1, clear wbCleanup; end                                     % trigger onCleanup to close the bar

results.MiddleH_Big_cum=squeeze(results.irf_nationals_cum);
results.MiddleH_Big=squeeze(results.irf_nationals);

results.rhoU = rhoU_all{end};
results.rhoU_by_rotation = rhoU_all;

col_names = {'F', 'F-Olea'};
results.F_statistics = array2table(F_stat, 'VariableNames', col_names);
results.EXP_VAR = EXP_VAR;
results.r_factors.BN = r_bai_ng;
results.r_factors.ABC1 = rhat1_abc;
results.r_factors.AH1 = r_ER;
results.r_factors.AH2 = r_GR;
results.r_factors.onatski = out;


end

function print_first_stage_f_statistics(F_stat, col_names, rotation_names)

disp(' ')
disp('  ===============================================')
disp('  = First-stage F-statistics for the instrument =')
disp('  ===============================================')
disp(' ')

f_tol = 1e-10;
same_f_stat = all(abs(F_stat(:,1) - F_stat(1,1)) < f_tol) && ...
              all(abs(F_stat(:,2) - F_stat(1,2)) < f_tol);

if same_f_stat
    disp(array2table(F_stat(1,:), 'VariableNames', col_names, 'RowNames', {'All_Rotations'}));
else
    disp(array2table(F_stat, 'VariableNames', col_names, 'RowNames', rotation_names(:)));
end

end

function print_explained_variance_table(EXP_VAR)

countries = {'EA','AT','BE','DE','EL','ES','FR','IE','IT','NL','PT'};
series_codes = {'IPMN','HICPOV','LTIRT','SHIX','UNETOT'};
headers = {'IP: Manufacturing','HICP: Overall','10-years Interest Rate','Stock Price Index','Unemployment Rate'};
sub_headers = {'(IPMN)','(HICPOV)','(LTIRT)','(SHIX)','(UNETOT)'};

row_names = EXP_VAR.Properties.RowNames;
upper_row = find(strcmp(row_names, 'Upper68'), 1);
point_row = find(strcmp(row_names, 'ExpVar'), 1);
lower_row = find(strcmp(row_names, 'Lower68'), 1);

if isempty(upper_row) || isempty(point_row) || isempty(lower_row)
    warning('EXP_VAR does not contain the expected Upper68, ExpVar, and Lower68 rows. Skipping explained variance table printout.')
    return
end

disp(' ')
disp('  ==============================================================')
disp('  = Table 4: Share of Explained Variance by the common factors =')
disp('  ==============================================================')
disp(' ')

fprintf('  %-8s', 'Country');
for j = 1:numel(headers)
    fprintf('%24s', headers{j});
end
fprintf('\n');

fprintf('  %-8s', '');
for j = 1:numel(sub_headers)
    fprintf('%24s', sub_headers{j});
end
fprintf('\n');

for i = 1:numel(countries)
    point_values = cell(1,numel(series_codes));
    ci_values = cell(1,numel(series_codes));

    for j = 1:numel(series_codes)
        var_name = explained_variance_variable_name(series_codes{j}, countries{i});
        var_idx = find(strcmp(EXP_VAR.Properties.VariableNames, var_name), 1);

        if isempty(var_idx)
            point_values{j} = '-';
            ci_values{j} = '(-)';
        else
            point = EXP_VAR{point_row,var_idx};
            lower = EXP_VAR{lower_row,var_idx};
            upper = EXP_VAR{upper_row,var_idx};
            point_values{j} = sprintf('%0.2f', point);
            ci_values{j} = sprintf('(%0.2f-%0.2f)', lower, upper);
        end
    end

    fprintf('  %-8s', countries{i});
    for j = 1:numel(point_values)
        fprintf('%24s', point_values{j});
    end
    fprintf('\n');

    fprintf('  %-8s', '');
    for j = 1:numel(ci_values)
        fprintf('%24s', ci_values{j});
    end
    fprintf('\n');
end

fprintf('  Notes: entries report the share of variance explained by the common component; parentheses contain 68%% confidence bands.\n');

end

function var_name = explained_variance_variable_name(series_code, country)

if strcmp(country, 'EA') && strcmp(series_code, 'LTIRT')
    var_name = 'LTIRT_EACC';
elseif strcmp(country, 'EA')
    var_name = [series_code, '_EA'];
else
    var_name = [series_code, '_', country];
end

end
