function[parameters,K1Q_X,rinfQ,sigma_cP,L0,L1,hess,K0Q_cP, K1Q_cP,sigma_e,para1]=...
    ML_estimation(yields,L,mature,dt)

options  =  optimset('fminunc');
options  =  optimset(options , 'TolFun'      , 1e-20);
options  =  optimset(options , 'TolX'        , 1e-15);
options  =  optimset(options , 'Display'     , 'on');
options  =  optimset(options , 'Diagnostics' , 'off');
options  =  optimset(options , 'LargeScale'  , 'off');
options  =  optimset(options , 'MaxFunEvals' , 20000) ;
options  =  optimset(options , 'MaxIter' , 10000) ;



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:3)';
pc1 = yields*W';


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


K0P_cP= para1(end,:)';
K1P_cP= para1(1:end-1,:)';
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));
K1Q_X = [.99,.98,.9];
parameters0 = convert_ind2v(sort(K1Q_X,'descend'),rinfQ,sigma_cP0);
% optimization
disp('optimization');
% [parameters,FVAL,EXITFLAG,OUTPUT] =  fminunc('jsz_like', parameters, options, yields, W, mature, dt);
[parameters1,FVAL1,EXITFLAG,OUTPUT] =  fminsearch('jsz_like_L_macro', parameters0, options, yields, W, mature, dt,u1t_star',L);
%CSMINWEL
 [FVAL,parameters,gh,hess,itct,fcount,retcodeh] = ...
     csminwel('jsz_like_L_macro', parameters1,eye(length(parameters1))*.5,[],1e-15,10000,yields, W, mature, dt,u1t_star',L);

 [LLF,llk,penalty,sigma_e,K0Q_cP, K1Q_cP] = ....
     jsz_like_L_macro(parameters,yields, W, mature, dt, u1t_star',L);
[K1Q_X,rinfQ,sigma_cP] = convert_v2ind(parameters);


L0=(K0P_cP-K0Q_cP);
L1=(K1P_cP-(K1Q_cP+eye(N)));
