function IV_F= my_first_stage_F_jae(E_opt, T_eff, Params, Impact_IRF,V_opt)
 

Params.N=size(E_opt,3);
PREDICTORS_Y=ones((T_eff)*Params.N,Params.p+1);
PREDICTORS_Z=ones((T_eff)*Params.N,Params.pz+1);
for pp=1:Params.p
    PREDICTORS_Y(:,pp)=reshape(E_opt(pp,1:T_eff,:), (T_eff)*Params.N,1);
    if pp<=Params.pz
        PREDICTORS_Z(:,pp)=reshape(E_opt(Params.p+pp,1:T_eff,:), (T_eff)*Params.N,1);
    end;
end;

Q_IV=chol(V_opt(1:Params.p,1:Params.p),'lower')\Impact_IRF;
Q=[Q_IV,null(Q_IV')];
Full_IRF=chol(V_opt(1:Params.p,1:Params.p),'lower')*Q; 
%transform errors into structural shocks
PREDICTORS_Y(:,1:Params.p)=PREDICTORS_Y(:,1:Params.p)/Full_IRF;

Winsor_Cutoff = 0;
for pp=1:Params.p
    PREDICTORS_Y(PREDICTORS_Y(:,pp)<prctile(PREDICTORS_Y(:,pp), Winsor_Cutoff) ,pp) = prctile(PREDICTORS_Y(:,pp), Winsor_Cutoff);
    PREDICTORS_Y(PREDICTORS_Y(:,pp)>prctile(PREDICTORS_Y(:,pp), 100-Winsor_Cutoff) ,pp) = prctile(PREDICTORS_Y(:,pp), 100-Winsor_Cutoff);
end;

WinsorE=NaN((T_eff)*Params.N,Params.pz);


for ppzz=1:Params.pz
    WinsorE(:,ppzz)=reshape(E_opt(Params.p+ppzz,1:(T_eff),:), (T_eff)*Params.N,1);
    WinsorE(WinsorE(:,ppzz)<prctile(WinsorE(:,ppzz),Winsor_Cutoff),ppzz) = prctile(WinsorE(:,ppzz),Winsor_Cutoff);
    WinsorE(WinsorE(:,ppzz)>prctile(WinsorE(:,ppzz),100-Winsor_Cutoff),ppzz) = prctile(WinsorE(:,ppzz),100-Winsor_Cutoff);
    if ppzz<=Params.pv1
        PREDICTORS_Z(PREDICTORS_Z(:,ppzz)<prctile(PREDICTORS_Z(:,ppzz),Winsor_Cutoff),ppzz)=prctile(PREDICTORS_Z(:,ppzz),Winsor_Cutoff);
        PREDICTORS_Z(PREDICTORS_Z(:,ppzz)>prctile(PREDICTORS_Z(:,ppzz),100-Winsor_Cutoff),ppzz)=prctile(PREDICTORS_Z(:,ppzz),100-Winsor_Cutoff);
    end;
end;

IV_R2=NaN(1,Params.pz); % -> regress each structural shocks on instruments. 

for ppzz=1:Params.pz
    PredictedB=PREDICTORS_Y(:,[1:Params.pv1, Params.pz+1])\WinsorE(:,ppzz);
    PredictedVal=PREDICTORS_Y(:,[1:Params.pv1, Params.pz+1])*PredictedB;
    IV_R2(1,ppzz)=corr(PredictedVal,reshape(E_opt(Params.p+ppzz,1:(T_eff),:), (T_eff)*Params.N,1))^2;
end;

IV_F= ( IV_R2' ./(1-IV_R2') ) * ((T_eff)*Params.N - Params.p - 1 ) / Params.p;
