function out = se_irf_fevd(datain,fac_est_out,est_par,...
    hor,identification_par,scl_par,n_rep,inclcode,tcode)
% Compute standard errors for IRFs and FEVDs using (Gaussian) parametric bootstrap.
% This is adapated from James H. Stock and Mark W. Watson (2016).

% Bootstrap Esimation Parameters
% ... Factors (including observed factors are simulated
bs_est_par = est_par;

bs_est_par.var_par.selectlagorder=0;

bs_est_par.fac_par.w = 1;

% ------ Easy names for paramters ----
[dnobs, n_fac] = size(fac_est_out.fac);
n_series      = size(datain,2);
uar_coef_mat = fac_est_out.uar_coef_mat;
uar_ser_mat = fac_est_out.uar_ser_mat;
M = fac_est_out.varout.coef.M;
G = fac_est_out.varout.coef.G0;
nshocks=size(fac_est_out.varout.coef.G,2);
n_uarlag = est_par.n_uarlag;
n_varlag = fac_est_out.p;
n_initial = 300;      % number of discarded initial values for VAR and AR simulations;

% Matrices Used below;
bpdata_NaN = isnan(datain);
Mu = zeros(n_uarlag,n_uarlag);
if n_uarlag > 1;
    Mu(2:end,1:end-1) = eye(n_uarlag-1);
end;
Gu = zeros(n_uarlag,1);

% --- Matrices for saving sum and sum^2 of draws
sum_btstrp_imp_y_fac_mat = zeros(n_series,nshocks,hor);
sum_squared_btstrp_imp_y_fac_mat = zeros(n_series,nshocks,hor);
sum_btstrp_imp_y_fac_mat_scl = zeros(n_series,nshocks,hor);
sum_squared_btstrp_imp_y_fac_mat_scl = zeros(n_series,nshocks,hor);
sum_btstrp_vfrac_y_comp_mat = zeros(n_series,nshocks,hor);
sum_squared_btstrp_vfrac_y_comp_mat = zeros(n_series,nshocks,hor);
sum_btstrp_vfrac_y_comp123_mat = zeros(n_series,hor);
sum_squared_btstrp_vfrac_y_comp123_mat = zeros(n_series,hor);
for irep = 1:n_rep
    t1 = cputime;
    % ---- Generate Bootstap Data
    % Generate Factor;
    f_btstrp = zeros(n_fac*n_varlag,n_initial+dnobs);
    
    for t = 2:dnobs+n_initial       % bootstrap sample of {f(t)}
        f_btstrp(:,t) =  M*f_btstrp(:,t-1) + G*randn(size(G,2),1);
    end
    f_btstrp = f_btstrp(1:n_fac,n_initial+1:end);  % Note this is nfac X dnobs;
    
    % Generate u
    u_btstrp = NaN(dnobs,n_series);
    for i = 1:n_series
        Mu(1,:) = uar_coef_mat(i,:);
        Gu(1,1) = uar_ser_mat(i);
        u = zeros(max([n_uarlag,1]),n_initial+dnobs);
        for t = 2:dnobs+n_initial
            u(:,t) =  Mu*u(:,t-1) + Gu*randn(1,1);
        end
        u_btstrp(:,i) = u(1,n_initial+1:end)';
    end
    % Generate Y
    y_btstrp = (fac_est_out.lam_mat*f_btstrp)' + u_btstrp;
    % Replace values with missing values as in original dataset
    y_btstrp(bpdata_NaN==1) = NaN;
    % Save oberved factors
    if bs_est_par.fac_par.nfac.observed > 0
        bs_est_par.fac_par.w = f_btstrp(bs_est_par.fac_par.nfac.unobserved1+1:...
            bs_est_par.fac_par.nfac.unobserved1+bs_est_par.fac_par.nfac.observed,:)';        
%         bs_est_par.fac_par.w = f_btstrp(1:bs_est_par.fac_par.nfac.observed,:)';        
    end
    
    % Compute VD and IRF
    btstrp_fac_est_out = factor_estimation_ls_full_simple(y_btstrp,inclcode,bs_est_par);
    btstrp_fac_est_out.varout.coef.G0=btstrp_fac_est_out.varout.coef.G;
    [Hy_btstrp,H_btstrp]=max_fev_share(btstrp_fac_est_out,bs_est_par,hor,tcode,identification_par);
    irf_fevd_out_btstrp = irf_fevd_full(Hy_btstrp,H_btstrp,scl_par,btstrp_fac_est_out,bs_est_par,hor,tcode,identification_par.res_ynews);  

    sum_btstrp_imp_y_fac_mat = sum_btstrp_imp_y_fac_mat + irf_fevd_out_btstrp.imp_y_fac_mat;
    sum_squared_btstrp_imp_y_fac_mat = sum_squared_btstrp_imp_y_fac_mat+ irf_fevd_out_btstrp.imp_y_fac_mat.^2;
    sum_btstrp_imp_y_fac_mat_scl = sum_btstrp_imp_y_fac_mat_scl + irf_fevd_out_btstrp.imp_y_fac_mat_scl;
    sum_squared_btstrp_imp_y_fac_mat_scl = sum_squared_btstrp_imp_y_fac_mat_scl+ irf_fevd_out_btstrp.imp_y_fac_mat_scl.^2;
    sum_btstrp_vfrac_y_comp_mat = sum_btstrp_vfrac_y_comp_mat + irf_fevd_out_btstrp.vfrac_y_comp_mat;
    sum_squared_btstrp_vfrac_y_comp_mat =  sum_squared_btstrp_vfrac_y_comp_mat + irf_fevd_out_btstrp.vfrac_y_comp_mat.^2;

    tmp=irf_fevd_out_btstrp.vfrac_y_comp_mat(:,[1,2,3],:);
    tmp=reshape(sum(tmp,2),n_series,size(irf_fevd_out_btstrp.vfrac_y_comp_mat,3));
    sum_btstrp_vfrac_y_comp123_mat = sum_btstrp_vfrac_y_comp123_mat + tmp;
    sum_squared_btstrp_vfrac_y_comp123_mat =  sum_squared_btstrp_vfrac_y_comp123_mat + tmp.^2;
    
    
    if irep == 5
        t2 = cputime-t1;
        fprintf('Carry out simulations to compute standard errors \n');
        fprintf('  Number of simulations: %6i \n',n_rep);
        fprintf('  Seconds per simulation: %6.3f \n',t2);
        fprintf('  Total time (in seconds) required %10.2f \n',t2*n_rep);
    end
end

mean_btstrp_imp_y_fac_mat = sum_btstrp_imp_y_fac_mat/n_rep;
var_btstrp_imp_y_fac_mat = sum_squared_btstrp_imp_y_fac_mat/n_rep - mean_btstrp_imp_y_fac_mat.^2;
mean_btstrp_imp_y_fac_mat_scl = sum_btstrp_imp_y_fac_mat_scl/n_rep;
var_btstrp_imp_y_fac_mat_scl = sum_squared_btstrp_imp_y_fac_mat_scl/n_rep - mean_btstrp_imp_y_fac_mat_scl.^2;
mean_btstrp_vfrac_y_comp_mat = sum_btstrp_vfrac_y_comp_mat/n_rep;
var_btstrp_vfrac_y_comp_mat =  sum_squared_btstrp_vfrac_y_comp_mat/n_rep - mean_btstrp_vfrac_y_comp_mat.^2;
mean_btstrp_vfrac_y_comp123_mat = sum_btstrp_vfrac_y_comp123_mat/n_rep;
var_btstrp_vfrac_y_comp123_mat =  sum_squared_btstrp_vfrac_y_comp123_mat/n_rep - mean_btstrp_vfrac_y_comp123_mat.^2;

se_btstrp_imp_y_fac_mat = sqrt(var_btstrp_imp_y_fac_mat);
se_btstrp_imp_y_fac_mat_scl = sqrt(var_btstrp_imp_y_fac_mat_scl);
se_btstrp_vfrac_y_comp_mat = sqrt(var_btstrp_vfrac_y_comp_mat);
se_btstrp_vfrac_y_comp123_mat = sqrt(var_btstrp_vfrac_y_comp123_mat);

%output
out.mean_imp_y_fac_mat = mean_btstrp_imp_y_fac_mat;
out.mean_imp_y_fac_mat_scl = mean_btstrp_imp_y_fac_mat_scl;
out.mean_vfrac_y_comp_mat = mean_btstrp_vfrac_y_comp_mat;
out.mean_vfrac_y_comp123_mat = mean_btstrp_vfrac_y_comp123_mat;
out.se_imp_y_fac_mat = se_btstrp_imp_y_fac_mat;
out.se_imp_y_fac_mat_scl = se_btstrp_imp_y_fac_mat_scl;
out.se_vfrac_y_comp_mat = se_btstrp_vfrac_y_comp_mat;
out.se_vfrac_y_comp123_mat = se_btstrp_vfrac_y_comp123_mat;

end
