
%- Number of factors
r1=5; nfac.unobserved1 = r1;% number of the yield curve factors 
r2=8; nfac.unobserved2 = r2;% number of the macro factors 
est_par.fac_par.w = [];     % don't use if all the factors are treated as unobserved -- here we treat all the factors as unobserved
nfac.observed = 0;          % number of observed factors
r = r1+r2; nfac.unobserved = r;
nfac.total = nfac.unobserved+nfac.observed;
est_par.fac_par.nfac = nfac;

%- Identification here is for yield news shock, orthogonal to realized
% volatility and business cycle news shocks
% First set horizon for max FEV share problem 
k1=0; % horizon for realized volatility shock
k2=11; % horizon for business cycle news shock

% Target variable
m1=find(ismember(name_vec,'RVol'));
m2=[];
for i=1:length(leadlist)
    m2=[m2,find(ismember(name_vec,leadlist(i)))];
end

identification_par=[];
identification_par.k_yield_news_shock=11; % horizon used for identification of yield news shock
identification_par.k{1}=k1;
identification_par.m{1}=m1;
identification_par.k{2}=k2;
identification_par.m{2}=m2;
identification_par.res_ynews=1;

%- Note, scale of shocks set so that they each lead to the same peak decline of the two-year yield 
scl_par=[];
indstrvar=find(ismember(name_vec,{'YIELD24'})); 
scl_par.indstrvar=indstrvar; % This contains the index of the two-year yield

% horizon for IRFs and FEVDs
hor=121;  

%-- Construct estimates of factors
fac_est_out =factor_estimation_ls_full_simple(data,inclcode_vec,est_par);
fac_est_out.varout.coef.G0=fac_est_out.varout.coef.G;

%-- Construct estimates of IRFs and FEVDs
[Hy,H]=max_fev_share(fac_est_out,est_par,hor,tcode_vec,identification_par);
irf_fevd_out = irf_fevd_full(Hy,H,scl_par,fac_est_out,est_par,hor,tcode_vec,identification_par.res_ynews);   

%-- Compute SEs for IRFs and FEVDs
se_irf_vdecomp_out = se_irf_fevd(data,fac_est_out,est_par,...
    hor,identification_par,scl_par,n_rep,inclcode_vec,tcode_vec);
