function [pc, common_components, finaldata, A,common_components0, EXP_VAR] = pca_monthly(data,date_EU, r, header_rotations)

% Standardize data
[finaldata, ~] =standardize(data);

% plot standardized data
%figure; plot(date_EU, finaldata); sgtitle('Standardized EU Aggregate data');


% Estimate Factors and Common Components
Sigma = cov(finaldata);
opts.disp = 0;
[ A, L ] = eigs(Sigma,r,'LM',opts); L = diag(L); % A --> E-Vectors / L --> E-Values

% fi
pc = finaldata*A; % Factors EA

% Variance Factors
Exp_r_static = zeros(1,r);
for s=1:r
    spazio = L(1:s);
    ratio = sum(spazio)/trace(Sigma);
    Exp_r_static(1,s) = ratio;
end

% Common Components EA
common_components=pc*A'; common_components0=pc*A';
common_components=common_components.*std(data); % de-standardize
common_components = common_components + mean(data);% de-demean data
idio = data - common_components;

variance_cc = var(common_components);
variance_tot = var(data);

exp_var = variance_cc./variance_tot;

Exp_var = exp_var;

T = size(pc,1);

%% Confidence Intervals

n_boot = 1000;
exp_var_boot = nan(n_boot,size(idio,2));

for boot_iter = 1:n_boot

    % Factors (indipendent resampling)
    f_bs = NaN(T, r);
    w = 12;
    for f_j = 1:r
        [f_bs(:,f_j), ~]=stationary_bootstrap(pc(:,f_j),w);
    end

    common_components_bs = f_bs*A';
    common_components_bs = common_components_bs.*std(data); % de-standardize
    common_components_bs = common_components_bs + mean(data);% de-demean data

    % New data
    data_new = common_components_bs + idio;

    [finaldata_bs, ~] =standardize(data_new);
    Sigma_bs = cov(finaldata_bs);
    opts.disp = 0;
    [ A_bs, L_bs ] = eigs(Sigma_bs,r,'LM',opts); L_bs = diag(L_bs); % A --> E-Vectors / L --> E-Values
    
    % fi
    pc_bs = finaldata_bs*A_bs; % Factors EA

    cc_bs=pc_bs*A_bs';
    cc_bs=cc_bs.*std(data_new); % de-standardize
    cc_bs = cc_bs + mean(data_new);% de-demean data

    % Explained Variance CC
    variance_cc_new = var(cc_bs);
    variance_tot_new = var(data_new);
    exp_var_boot(boot_iter,:) = variance_cc_new./variance_tot_new;


end

Exp_var_boot = exp_var_boot;

conf_l = 68;

upper_68 = prctile(Exp_var_boot, (100+conf_l)/2, 1);
lower_68 = prctile(Exp_var_boot, (100-conf_l)/2, 1);


%%

table = [upper_68; Exp_var; lower_68];

col_names_v = header_rotations;

for i = 1:r
    factors_names{i} = ['R' num2str(i)];
end

col_names = [col_names_v];

row_names = {'Upper68', 'ExpVar', 'Lower68'};
EXP_VAR = array2table(table, 'VariableNames', col_names, 'RowNames', row_names);



end
