function[parameters,K1Q_X,rinfQ,sigma_cP...
    ,hess,K0Q_cP, K1Q_cP,sigma_e,para2,u1t_star]=...
    ML_estimation_5_RR(yields,L,mature,dt,macro,NFACT,varargin)

options=optimset('MaxFunEvals',500000,'MaxIter',500000,'Display','Off','Algorithm','active-set');

if isempty(varargin)
    rnk=1;
else
    rnk=varargin{1};
end
cov_mat = cov(yields);
[eig_vec,eig_val] = eig(cov_mat);
eig_vec = eig_vec(:,end:-1:1);
% eig_vec(:,1) = -eig_vec(:,1);

W = eig_vec(:,1:NFACT)';
pc1 = yields*W';
N1=cols(pc1);

dat=[pc1 macro];
y1 = dat(L+1:end,:);
N=cols(y1);
x1 = prepare(dat,L);
x1=x1(L+1:end,:);
T=rows(y1);
para1 = x1\y1;

u1t_star = y1 - x1*para1;
var_1 = (u1t_star'*u1t_star)/(rows(y1)-cols(x1));
sigma_cP0 = chol(var_1)';

% Q-dynamics
% starting values
rinfQ = mean(yields(:,end));
if NFACT==3
K1Q_X = [.99,.98,.9];
elseif NFACT==4
K1Q_X = [.99,.98,.95,.85];
elseif NFACT==5
 K1Q_X = [.99,.98,.95,.85,.75];
else
  K1Q_X = [.99,.98,.95,ones(1,NFACT-3).*0.75];
end
parameters0 = convert_ind2v(sort(K1Q_X,'descend'),rinfQ,sigma_cP0);
disp('optimization');
[parameters,FVAL1,EXITFLAG,OUTPUT] =  ...
    fminsearch('jsz_like_L_macroRR', parameters0, options, yields, W, mature, dt,NFACT,rnk);
%  [FVAL,parameters1,gh,hess,itct,fcount,retcodeh] = ...
%      csminwel('jsz_like_L_macroRR', parameters,...
%      eye(length(parameters))*.0000001,[],1e-10,10000,yields, W, mature, dt,NFACT);
%MCMC to get variance
[LLF,llk,sigma_e,K0Q_cP, K1Q_cP,K0P_cP, K1P_cP] = ...
    jsz_like_L_macroRR(parameters,yields, W, mature, dt,NFACT,rnk);
likold=-LLF;
K=length(parameters);
naccept=0;
PP=diag(abs(parameters))*0.0000000001;
PP=chol(PP);
bold=parameters;
out=[];
for j=1:10000
    bnew=bold+(randn(1,K)*PP)';
    liknew=...
   jsz_like_L_macroRR(bnew,yields, W, mature, dt,NFACT,rnk);
    liknew=-liknew;
    accept=exp(liknew-likold);
    if accept>rand
        bold=bnew;
        naccept=naccept+1;
        likold=liknew;
    end
    out=[out;bold'];
    num2str(naccept/j)
end
hess=cov(out(1001:end,:));

 [K1Q_X,rinfQ,sigma_cP] = convert_v5ind(parameters,NFACT);
 omega=sigma_cP*sigma_cP';
cP = yields*W'; % (T+L)*N, cP stands for math caligraphic P.
[BcP, AcP, AX, BX, DUMMY,penalty] = ...
    jszLoadings(W, diag(K1Q_X-1), rinfQ, omega(1:NFACT,1:NFACT), mature, dt);
[K0Q_cP, K1Q_cP, rho0_cP, rho1_cP] = ...
    jszRotation(W, diag(K1Q_X-1), rinfQ, dt, [], [], BX, AX);
  
para2=[K1P_cP K0P_cP];



