function [MiddleH_cum, MiddleH, difference_national_ea, difference_national_ea_cum, intervals,  F, rhoU]=irf_iv(st_f_ccvar,vardata, n,p,c,hor,mpi, date_EU,date_iv_nonan, conf)

[C,~,~,~,~,~,~,err]=wold_irf_LP(vardata,n,p,c,hor,st_f_ccvar); % reduced-form errors 

% External instrument - Two stage regression:
[iv_index, mpi] = align_proxy_sample(date_EU,date_iv_nonan,mpi,p,err);
startsample = min(iv_index)-1;
endsample = size(err,1)-max(iv_index);

% Consider the errors overlapping with the instrument's observations
eps_p=err(iv_index,1); % reduced-form errors mp shock (mpi overlapping period only)           
eps_q=err(iv_index,2:end); % reduced-form errors other shocks (mpi overlapping period only)

%% ==========================================================================
% Autocorrelation Functions of VAR Residuals
%==========================================================================

%acf of u(t)
rhoU = acf([eps_p,eps_q],7);

%% Define the instrument:
Z = mpi; 

% === STEP 1: Normalize the instrument ===

% Z = Z - mean(Z);                
% [Q, ~] = qr(Z / sqrt(length(Z)), 0);     % QR decomposition: makes Q orthonormal with respect to (1/S)
% Z_ort = Q * sqrt(length(Z));

% Lewis et al.
Z_ort = orthonormalize_proxy(Z); 
T = size(Z_ort,1);
%(Z_ort'*Z_ort)/T


% ==== Step 1st stage regression
reg_temp = fitlm(Z_ort,eps_p); u_temp = table2array(reg_temp.Residuals(:,1));
ols_instr_ort=reg_temp.Coefficients.Estimate;
F.F_statistic = first_stage_fstat(reg_temp,T,size(Z_ort,2));
% === STEP 2: Compute Y'ZZ'Y ===
T = length(eps_p);
num = (eps_p' * (Z_ort * Z_ort') * eps_p); 

% === STEP 3: estimate u variance
Zv2 = Z_ort.*repmat(u_temp,1,size(Z_ort,2)); 

W2_nw = nw(Zv2, T); 

% === STEP 4: Compute Effective F-statistic ===
F.F_eff = num / (T * trace(W2_nw));

%%

% Compute fitted u_{t}^{p}:
u_hat_p=[ones(size(Z_ort,1),1),Z_ort]*ols_instr_ort;

% Second step - Estimate s_{q}/s_{p} and normalise s_{p}=1 in order to have
% an impact response of the policy rate of 1% after the monetary policy
% shock:

sq_sp=(u_hat_p'*u_hat_p)\u_hat_p'*eps_q; 
%[sq_sp_2]=OLS(eps_q,u_hat_p,0); % check

% Step 3 - Normalise s_p=1 and compute IRFs

% Structural coefficients for the monetary policy shock
s=[1,sq_sp]';

% Compute the IRFs:

H=zeros(n,1,hor);
for i=1:hor
    H(:,:,i)=C(:,:,i)*s; % IV wold respesentation
end

% Organise the IRFs in a hor x n*n matrix:

MiddleH=(reshape(permute(H,[3 2 1]),hor,n,[]));

cc=1; % constant for the two stages regression

Z = NaN(size(err,1),size(mpi,2));
Z(iv_index,:) = mpi;

% cumulated

MiddleH_cum=cumsum(MiddleH,1);

national_single_irf = MiddleH(:,end);

for jjj = 1:size(MiddleH,2)-1
difference_national_ea(:,jjj) = national_single_irf - MiddleH(:,jjj);
end

national_single_irf_cum = MiddleH_cum(:,end);

for jjj = 1:size(MiddleH,2)-1
difference_national_ea_cum(:,jjj) = national_single_irf_cum - MiddleH_cum(:,jjj);
end

% Bootstrap:

nboot=10000;

if conf == 1
    [intervals]=bootstrapProxy_LP_rotation_b(st_f_ccvar,vardata,hor,c,cc,nboot,68,p,n,Z,startsample,endsample,iv_index);
elseif conf == 0

    intervals.HighH = nan(hor,n);
    intervals.LowH = nan(hor,n);

    %--------------------------------

    intervals.HighH_cum = nan(hor,n);
    intervals.LowH_cum = nan(hor,n);

    %--------------------------------

    intervals.HighH_diff = nan(hor,n);
    intervals.LowH_diff = nan(hor,n);

    %--------------------------------

    intervals.HighH_diff_cum = nan(hor,n);
    intervals.LowH_diff_cum = nan(hor,n);

end

end

function [iv_index, mpi] = align_proxy_sample(date_EU,date_iv_nonan,mpi,p,err)

if size(mpi,1) ~= size(date_iv_nonan,1)
    error('The proxy instrument and proxy date vector have different numbers of observations.');
end

valid_mpi = all(~isnan(mpi),2);
date_iv_nonan = date_iv_nonan(valid_mpi,:);
mpi = mpi(valid_mpi,:);

err_dates = date_EU(p+1:end);
[in_sample, iv_index] = ismember(date_iv_nonan,err_dates);
if ~any(in_sample)
    error('No valid proxy observations overlap with VAR residual dates.');
end

iv_index = iv_index(in_sample);
mpi = mpi(in_sample,:);

if any(iv_index < 1) || any(iv_index > size(err,1))
    error('Proxy dates map outside the reduced-form residual sample.');
end

end

function Z_ort = orthonormalize_proxy(Z)

Z_centered = Z - mean(Z,1);
G = (Z_centered'*Z_centered)/size(Z_centered,1);
G = (G + G')/2;
[V,D] = eig(G);
d = real(diag(D));
tol = max(size(G))*eps(max(d));

if any(d <= tol)
    error('The proxy instrument matrix is rank deficient after demeaning.');
end

Z_ort = Z_centered*V*diag(1./sqrt(d))*V';

end

function F_stat = first_stage_fstat(model,T,k)

R2 = model.Rsquared.Ordinary;
if k <= 0 || T <= k + 1 || R2 >= 1
    F_stat = NaN;
else
    F_stat = (R2/k)/((1-R2)/(T-k-1));
end

end








