%% Date: 30/10/2025
%% Table 6
cd(ini.dirs.code)
%% Initialize all data matrix
load ccydates; 
T=size(ccydates,1);
%% intervention
load intervention; load interventiondates;
idx=ismember(ccydates,interventiondates);
int=NaN(T,1);
int(idx)=intervention;
int(isnan(int))=0;

%% dates
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));

%% RR & BF
load('lagrr'); load('lagbf');
%% CDS
load('lagcds')

%% ccys
load ilsnyccy; load dilsnyccy;
%load ilslondonccy;
load NEERny;

%% Change in NEER
load lagNEER1
load lagNEER2
load lagNEER5
load lagNEER60
load lagNEER10
load newforwards
load newbasis
%load lagNEER3

%% Flows
load laginstflows
load lagforeignflows


%%%% Exogenous
%% lagged interest rate change
%% 5-day change in interest rate
load uslibor; load telbor;
p=10;
lead=5;
dtelbor=[repmat(NaN,lead,6);telbor(lead+1:end,:)-telbor(1:end-lead,:)];
dlibor=[repmat(NaN,lead,6);uslibor(lead+1:end,:)-uslibor(1:end-lead,:)];
%% Take the three-month change
lagtelborchange=lagmatrix(dtelbor(:,2),0:p);
lagliborchange=lagmatrix(dlibor(:,2),0:p);
save('lagtelborchange','lagtelborchange');
save('lagliborchange','lagliborchange');

%1-day libor change 
lead=1;
dlibor=[repmat(NaN,lead,6);uslibor(lead+1:end,:)-uslibor(1:end-lead,:)];
dtelbor=[repmat(NaN,lead,6);telbor(lead+1:end,:)-telbor(1:end-lead,:)];
lagtelborchange1=lagmatrix(dtelbor(:,2),0:p);
%% Take the three-month change
lagliborchange1=lagmatrix(dlibor(:,2),0:p);
save('lagliborchange1','lagliborchange1');
save('lagtelborchange1','lagtelborchange1');


load lageurusd;
load intweeklag;

%%% dforward
dforwards=[repmat(NaN,1,5);100*log(newforwards(2:end,:)./newforwards(1:end-1,:))]; %% in %
save('dforwards','dforwards');
load newbasis;

%% verbal intervention
load verbalindex;
%% pub index
load pubindex;

%% verbal intervention
load verbalindex;
%% pub index
load pubindex;
%% vix
load vix; load vixdates;
vixfull=NaN(T,1);

% Find matching dates
[~, idx_vixdates, idx_ccydates] = intersect(vixdates, ccydates);

% Assign VIX values to matching dates
vixfull(idx_ccydates) = vix(idx_vixdates);
vixfull = fillmissing(vixfull,'previous');lagvix=lagmatrix(vixfull,0:p);
save('lagvix','lagvix');


%% Panel A:
%%  Regular OLS %%%
idx=find(ccydates==datenum('01/01/2013','dd/mm/yyyy')); %% beginning of sample
idx1=find(ccydates==datenum('31/12/2019','dd/mm/yyyy'));

load NSindex1; load NSindex2;
NSindex1lag=lagmatrix(NSindex1,0:20);
NSindex2lag=lagmatrix(NSindex2,0:20);

%% if interest rate
dCITI=CITIindex(2:end)-CITIindex(1:end-1);
dCITI=[NaN;dCITI];

%% With intermediatary constraints (from Zhuigo He's Website)
load hkm;
temphkm=hkm(:,4);
temphkm(temphkm==0)=NaN;
temphkm=fillmissing(temphkm,'previous');

% Load proprietary HKM files for robustness
load('..\_data\hkm\ilshkm_aggregated.mat');
load('..\_data\hkm\global_hkm_clean.mat');
load('..\_data\hkm\foreign_hkm_true');

%% Align it with ccydates

% Filter for specific dates in ccydates
% Convert ccydates from MATLAB serial numbers to datetime
ccydates_datetime = datetime(ccydates, 'ConvertFrom', 'datenum');

%% Initialize aligned data arrays
n_ccy_dates = length(ccydates_datetime);

% Initialize arrays for aligned data
ils_hkm_aligned = NaN(n_ccy_dates, 1);
foreign_hkm_aligned = NaN(n_ccy_dates, 1);
global_hkm_aligned = NaN(n_ccy_dates, 1);

%% Align Israeli HKM data
fprintf('\nAligning Israeli HKM data...\n');

for i = 1:n_ccy_dates
    target_date = ccydates_datetime(i);
    
    % Find the most recent Israeli HKM observation on or before the target date
    valid_idx = ilshkm_aggregated.validDatetime <= target_date;
    
    if any(valid_idx)
        % Get the most recent observation
        recent_dates = ilshkm_aggregated.validDatetime(valid_idx);
        recent_values = ilshkm_aggregated.validCapitalRatio(valid_idx);
        
        [~, most_recent_idx] = max(recent_dates);
        ils_hkm_aligned(i) = recent_values(most_recent_idx);
    end
end

ils_valid_count = sum(~isnan(ils_hkm_aligned));
fprintf('Israeli HKM: %d/%d dates filled (%.1f%%)\n', ...
        ils_valid_count, n_ccy_dates, 100*ils_valid_count/n_ccy_dates);

%% Align Foreign HKM data
fprintf('Aligning Foreign HKM data...\n');

for i = 1:n_ccy_dates
    target_date = ccydates_datetime(i);
    
    % Find the most recent foreign HKM observation on or before the target date
    valid_idx = foreign_hkm_true.Date <= target_date;
    
    if any(valid_idx)
        % Get the most recent observation (using filled data if available)
        recent_dates = foreign_hkm_true.Date(valid_idx);
        
        % Use filled data if available, otherwise use original
        if ismember('HKM_Ratio_Filled', foreign_hkm_true.Properties.VariableNames)
            recent_values = foreign_hkm_true.HKM_Ratio_Filled(valid_idx);
        else
            recent_values = foreign_hkm_true.HKM_Ratio(valid_idx);
        end
        
        [~, most_recent_idx] = max(recent_dates);
        if ~isnan(recent_values(most_recent_idx))
            foreign_hkm_aligned(i) = recent_values(most_recent_idx);
        end
    end
end

foreign_valid_count = sum(~isnan(foreign_hkm_aligned));
fprintf('Foreign HKM: %d/%d dates filled (%.1f%%)\n', ...
        foreign_valid_count, n_ccy_dates, 100*foreign_valid_count/n_ccy_dates);

%% Align Global HKM data (if available)
fprintf('Aligning Global HKM data...\n');

if exist('global_hkm_clean', 'var') && istable(global_hkm_clean)
    for i = 1:n_ccy_dates
        target_date = ccydates_datetime(i);
        
        % Find the most recent global HKM observation on or before the target date
        valid_idx = global_hkm_clean.Date <= target_date;
        
        if any(valid_idx)
            % Get the most recent observation
            recent_dates = global_hkm_clean.Date(valid_idx);
            recent_values = global_hkm_clean.Combined_HKM_Ratio(valid_idx);
            
            [~, most_recent_idx] = max(recent_dates);
            if ~isnan(recent_values(most_recent_idx))
                global_hkm_aligned(i) = recent_values(most_recent_idx);
            end
        end
    end
    
    global_valid_count = sum(~isnan(global_hkm_aligned));
    fprintf('Global HKM: %d/%d dates filled (%.1f%%)\n', ...
            global_valid_count, n_ccy_dates, 100*global_valid_count/n_ccy_dates);
else
    fprintf('Global HKM data not available or not in expected format\n');
end


%%%%%%%%%%%% Some correlation analysis for HKM measures %%%%%%%%%%%%

%% Create aligned results table
hkm_aligned = table();
hkm_aligned.Date = ccydates_datetime;
hkm_aligned.DateNum = ccydates; % Keep original datenum format
hkm_aligned.ILS_HKM = ils_hkm_aligned;
hkm_aligned.Foreign_HKM = foreign_hkm_aligned;
hkm_aligned.Global_HKM = global_hkm_aligned;

load hkm;
temphkm=hkm(:,4);
temphkm(temphkm==0)=NaN;

temphkm=fillmissing(temphkm,'previous');

plot(temphkm(idx:idx1)); hold on; plot(1./hkm_aligned.Global_HKM(idx:idx1).^2,'g')

%% Correlation Analysis: temphkm vs 1/HKM^2
fprintf('Correlation Analysis: temphkm vs 1/Global_HKM^2\n');
fprintf('===============================================\n');

% Extract the data for the specified range
series1 = temphkm(idx:idx1);
series2 = 1 ./ hkm_aligned.Global_HKM(idx:idx1).^2;

% Check data availability
fprintf('Data range: index %d to %d (%d observations)\n', idx, idx1, length(series1));

% Remove NaN values for correlation analysis
valid_data = ~isnan(series1) & ~isnan(series2) & isfinite(series1) & isfinite(series2);
series1_clean = series1(valid_data);
series2_clean = series2(valid_data);

fprintf('Valid observations (no NaN/Inf): %d/%d (%.1f%%)\n', ...
        sum(valid_data), length(series1), 100*sum(valid_data)/length(series1));

if sum(valid_data) < 2
    error('Insufficient valid data points for correlation analysis');
end

%% Levels Correlation
fprintf('\n=== LEVELS CORRELATION ===\n');
corr_levels = corr(series1_clean, series2_clean, 'Type', 'Pearson');
corr_levels_spearman = corr(series1_clean, series2_clean, 'Type', 'Spearman');

fprintf('Pearson correlation (levels): %.4f\n', corr_levels);
fprintf('Spearman correlation (levels): %.4f\n', corr_levels_spearman);

% Test significance
[~, p_value_levels] = corrcoef(series1_clean, series2_clean);
fprintf('P-value (levels): %.4f\n', p_value_levels(1,2));

if p_value_levels(1,2) < 0.01
    significance_levels = '***';
elseif p_value_levels(1,2) < 0.05
    significance_levels = '**';
elseif p_value_levels(1,2) < 0.10
    significance_levels = '*';
else
    significance_levels = '';
end

fprintf('Significance (levels): %s\n', significance_levels);

%% First Differences Correlation
fprintf('\n=== FIRST DIFFERENCES CORRELATION ===\n');

% Calculate first differences
diff1 = diff(series1_clean);
diff2 = diff(series2_clean);

% Remove any remaining NaN/Inf from differences
valid_diff = ~isnan(diff1) & ~isnan(diff2) & isfinite(diff1) & isfinite(diff2);
diff1_clean = diff1(valid_diff);
diff2_clean = diff2(valid_diff);

fprintf('Valid first differences: %d/%d (%.1f%%)\n', ...
        sum(valid_diff), length(diff1), 100*sum(valid_diff)/length(diff1));

if sum(valid_diff) < 2
    warning('Insufficient valid data points for first differences correlation');
    corr_diff = NaN;
    corr_diff_spearman = NaN;
    p_value_diff = NaN;
else
    corr_diff = corr(diff1_clean, diff2_clean, 'Type', 'Pearson');
    corr_diff_spearman = corr(diff1_clean, diff2_clean, 'Type', 'Spearman');
    
    fprintf('Pearson correlation (first differences): %.4f\n', corr_diff);
    fprintf('Spearman correlation (first differences): %.4f\n', corr_diff_spearman);
    
    % Test significance
    [~, p_value_diff] = corrcoef(diff1_clean, diff2_clean);
    fprintf('P-value (first differences): %.4f\n', p_value_diff(1,2));
    
    if p_value_diff(1,2) < 0.01
        significance_diff = '***';
    elseif p_value_diff(1,2) < 0.05
        significance_diff = '**';
    elseif p_value_diff(1,2) < 0.10
        significance_diff = '*';
    else
        significance_diff = '';
    end
    
    fprintf('Significance (first differences): %s\n', significance_diff);
end

%% Summary Statistics
fprintf('\n=== SUMMARY STATISTICS ===\n');
fprintf('Series 1 (temphkm):\n');
fprintf('  Mean: %.4f, Std: %.4f\n', mean(series1_clean), sqrt(var(series1_clean)));
fprintf('  Min: %.4f, Max: %.4f\n', min(series1_clean), max(series1_clean));

fprintf('Series 2 (1/HKM^2):\n');
fprintf('  Mean: %.4f, Std: %.4f\n', mean(series2_clean), sqrt(var(series2_clean)));
fprintf('  Min: %.4f, Max: %.4f\n', min(series2_clean), max(series2_clean));

%% Simple Plot of Two HKM Series with LaTeX Legend and Date Index

% Extract the data
series1 = temphkm(idx:idx1);
series2 = 1 ./ hkm_aligned.Global_HKM(idx:idx1).^2;

% Create date vector if available
if exist('ccydates', 'var') && length(ccydates) >= idx1
    x_axis = datetime(ccydates(idx:idx1), 'ConvertFrom', 'datenum');
else
    x_axis = idx:idx1;
end

% Create the plot
figure('Position', [100, 100, 1000, 600]);
plot(x_axis, series1, 'b-', 'LineWidth', 2.5);
hold on;
plot(x_axis, series2, 'r-', 'LineWidth', 2.5);

% Add labels with large fonts
xlabel('Date', 'FontSize', 20, 'FontWeight', 'bold');
ylabel('Value', 'FontSize', 20, 'FontWeight', 'bold');

% Create legend with LaTeX
legend({'HKM', '$\overline{HKM}$'}, 'Interpreter', 'latex', 'FontSize', 18, 'Location', 'best');
% Set axis font size
set(gca, 'FontSize', 16);

grid on;

%%%%%%%%%%%%%%%%%%%%%%

%% OLS %%

data=[dilsnyccy,ones(T,1),int/1000,temphkm,int/1000.*temphkm,lageurusd(:,1),lagvix(:,1)-lagvix(:,6),ILshockindex,NSindex1,NSindex2,CPIindex,CITIindex/1000]; %% dataset used for the estimation
[coeff,tstat,r,n]=fastnw(data(idx:idx1,1),data(idx:idx1,2:end),0);
x1=coeff(1),x2=coeff(2),x3=coeff(3);x4=coeff(4);x5=coeff(5);

%% GMM was ran in EVIEWS