function out = factor_estimation_ls_full_simple(data, inclcode, est_par)
% This program estimates the factors, their loadings, a VAR on the 
% estimated factors, ARs for error terms, and the R-squared of the 
% common components. 
% This is adapated from James H. Stock and Mark W. Watson (2016).

% PRELIMINARIES
n_series = size(data,2);                    % number of series
nfirst   = est_par.smpl_par.nfirst;         % start date
nlast    = est_par.smpl_par.nlast;          % end date
calvec   = est_par.smpl_par.calvec;         % calendar
nper     = est_par.smpl_par.nper;           % number of periods a year
n_uarlag = est_par.n_uarlag;                % number of AR lags
ntmin    = est_par.lambda.nt_min;           % minimum number of Obs

% USE SUBSET OF DATA TO ESTIMATE FACTORS
est_data = data(:,inclcode==1);
catcode=est_par.fac_par.catcode(inclcode==1);
lsout = factor_estimation_ls_simple(est_data,est_par,catcode);

% Construct estimates of factor loadings;
lam_mat = NaN(n_series,est_par.fac_par.nfac.total);
cons_vec = NaN(n_series,1);
ismpl = smpl(calvec,nfirst,nlast,nper);
uar_coef_mat = NaN(n_series,n_uarlag);
uar_ser_mat = NaN(n_series,1);
r2_mat = NaN(n_series,1);                              % R-squared value
for is = 1:n_series;
    tmp = packr([data(ismpl==1,is) lsout.fac(ismpl==1,:)]);
    if size(tmp,1) >= ntmin;
       y = tmp(:,1);
       x = [tmp(:,2:end), ones(size(tmp,1),1)];
       xxi = inv(x'*x);
       bols = xxi*(x'*y);
       b = bols;
       lam_mat(is,:) = b(1:end-1)';
       cons_vec(is)=b(end);
       u = y - x*b;
       % Compute R-squared 
       ssr = sum(u.^2);
       ym = y - mean(y);
       tss = sum(ym.^2);
       r2_mat(is) = 1-(ssr/tss);
       % Compute AR model for errors
       if r2_mat(is) < 0.9999;
           if n_uarlag>0
               [arcoef, ser] = uar(u,n_uarlag);    % AR Coefficients and ser
           else
               arcoef=0;
               ser=std(u);
           end
       else
           arcoef = zeros(max([n_uarlag,1]),1);
           ser = 0.0;
       end;
       uar_coef_mat(is,:) = arcoef'; 
       uar_ser_mat(is,1) = ser;
    end;
end;

% % Enforce average positive loadings
% negs=mean(lam_mat(est_par.fac_par.cacode==1,:),1)<0;
% lam_mat(:,negs)=-lam_mat(:,negs);
% lsout.fac(:,negs)=-lsout.fac(:,negs);

varout = varest(lsout.fac,est_par.var_par,est_par.smpl_par);

% SAVE OUTPUT
out.est_data     = est_data;
out.fac          = lsout.fac;
out.lam_mat      = lam_mat;
out.cons_vec     = cons_vec;
out.uar_coef_mat = uar_coef_mat;
out.uar_ser_mat  = uar_ser_mat;
out.varout       = varout;
out.r2           = r2_mat;
out.p            = est_par.var_par.nlag;
end