function [Hy,H]=max_fev_share(fac_est_out,est_par,hor,tcode,identification_par)

res_ynews  = identification_par.res_ynews;
G0 = fac_est_out.varout.coef.G0; % Cholesky factor
n_fac      = size(fac_est_out.fac,2);

%-- Construct estimates of the IRFs and FEVs from a recursive VAR 
irf_fevd_out = irf_fevd(fac_est_out,est_par,hor,tcode);  
imp_fac_mat=irf_fevd_out.imp_fac_mat; 
vtotal_fac_mat=irf_fevd_out.vtotal_fac_mat;
imp_y_fac_mat=irf_fevd_out.imp_y_fac_mat;
vtotal_y_mat=irf_fevd_out.vtotal_y_mat;
        
%-- Identification for the level, slope and yield news shocks
% Note, level and slope are the first two factors in the VAR.

m=1:2; % indexes of level and slope 
k=identification_par.k_yield_news_shock; % horizon for the max FEV share problem

d=vtotal_fac_mat(m,k+1);
D=diag((1./d));
S1=zeros(n_fac);
for l=0:k
    tmp=imp_fac_mat(:,:,l+1);
    S1=S1+tmp(m,:)'*D*tmp(m,:);    
end

Gy=eye(n_fac);
[G22,~,~]=svd(S1(m(2)+1:n_fac,m(2)+1:n_fac));
Gy(m(2)+1:n_fac,m(2)+1:n_fac)=G22;

Hy=G0*Gy;
Sy=S1;


%-- Identification for measures of implied and realized stock market 
% volatility, leading economic indicators, etc. 
Hs=[];
S=[];
if ~isempty(identification_par.m) 
    Gs=[];
    for i=1:length(identification_par.m) 
        m=identification_par.m{i}; % indexes of target variables
        k=identification_par.k{i}; % horizon for max FEV share problem 
        
        d=vtotal_y_mat(m,k+1);
        D=diag((1./d));
        
        S1=zeros(n_fac);
        for l=0:k
            tmp=imp_y_fac_mat(:,:,l+1);
            S1=S1+tmp(m,:)'*D*tmp(m,:);
        end
        S{i}=S1;
        
        [Gtmp,~,~]=svd(S1);
        Gs{i}=Gtmp;
        
        Hs{i}=G0*Gs{i};
        
    end
    
    H=Hs{1};
    
    
    %-- Identification for a residual yield news shock.
    % Note, shocks to realized volatility and business cycle news
    % are identified first. Then, a residual yield news shock is
    % identified by targeting the combined FEVs of level and slope, 
    % subject to the constraints that the shock is orthogonal to
    % realized volatility and business cycle news shocks and has zero 
    % effect on level and slope contemporaneously.
    if res_ynews==1 

        G1=Gs{1}(:,1); %-- G1 is the column of Gs{1} corresponding to the realized volatility shock
        G2=Gs{2}(:,1); %-- G2 is the column of Gs{2} corresponding to the business cycle news shock
        
        % convert to a minimization problem
        Q_m=-Sy; 
        
        % normalize to be orthonormal 
        H_m=[];k_v=[];d_s=[]; 
        H_m{1}=eye(size(Q_m,1));
        k_v{1}=zeros(size(Q_m,1),1);
        d_s{1}=-1/2;
        
        % impose orthogonality to G1 
        H_m{2}=zeros(size(Q_m,1)); 
        k_v{2}=G1;
        d_s{2}=0;
        
        % impose orthogonality to G2
        H_m{3}=zeros(size(Q_m,1)); 
        k_v{3}=G2;
        d_s{3}=0;        
        
        % impose contemporaneous zero restriction for level 
        H_m{4}=zeros(size(Q_m,1)); 
        k_v{4}=zeros(size(G1)); 
        k_v{4}(1)=1;
        d_s{4}=0;        
        
        % impose contemporaneous zero restriction for slope 
        H_m{5}=zeros(size(Q_m,1)); 
        k_v{5}=zeros(size(G1)); 
        k_v{5}(2)=1;
        d_s{5}=0;           
        
        options = optimset('Algorithm','interior-point',...
            'Display','off',...
            'GradObj','on','GradConstr','on',...
            'HessFcn',@(x,lambda)quadhess(x,lambda,Q_m,H_m));
        
        fun = @(x)quadobj(x,Q_m);
        nonlconstr = @(x)quadconstr(x,H_m,k_v,d_s);
        x0 = ones(size(Q_m,1),1)/size(Q_m,1); % column vector
        [x,fval,eflag,output,lambda] = fmincon(fun,x0,...
            [],[],[],[],[],[],nonlconstr,options);
                
        Gseq=x;        
        
        [Q,R]=qr(Gseq);
        Q(:,1)=Q(:,1)*sign(R(1));
        H=G0*Q;        

    end
else
    H=[];
end

end

function [y,grady] = quadobj(x,Q)
y = 1/2*x'*Q*x ;
if nargout > 1
    grady = Q*x;
end
end

function [y,yeq,grady,gradyeq] = quadconstr(x,H,k,d)
jj = length(H); % jj is the number of equality constraints

yeq = zeros(1,jj);
for i = 1:jj
    yeq(i) = 1/2*x'*H{i}*x + k{i}'*x + d{i};
end
y = [];

if nargout > 2
    gradyeq = zeros(length(x),jj);
    for i = 1:jj
        gradyeq(:,i) = H{i}*x + k{i};
    end
end
grady = [];
end

function hess = quadhess(x,lambda,Q,H)
hess = Q;
jj = length(H); % jj is the number of equality constraints
for i = 1:jj
    hess = hess + lambda.eqnonlin(i)*H{i};
end
end
