function [Est,f] = MDE_Timevarying_main_boot6(iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,init_points,Contraction_index,NumObs)  

XG=zeros(27,14);
opt = optimset('Display','off','GradObj','on','TolX',1e-7,'TolFun',1e-7,'MaxFunEvals',10000,'MaxIter',3000);

parfor  s=1:27
            X0=init_points(s,:);
            
            UB=[0               0       1  1 (log(1+0.9999)-log(1-0.9999))    0           0       1  1     0       1   0.05   0.05];
            LB=[log(1e-04)  log(1e-04) -1 -1  (log(1+0.3)-log(1-0.3))      log(1e-04) log(1e-04) -1 -1 log(1e-04) -1  -0.05  -0.05];       
            
   [params,Fval]= lsqnonlin('MDE_Timevarying_HigherOrderMom_boot6',X0,LB,UB,opt,iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,Contraction_index,NumObs);
   
    m2etaR = exp(params(1));
    m2etaE = exp(params(2));
    m3etaR = params(3);
    m3etaE =params(4);
    rho     = 2*exp(params(5))/(1+exp(params(5))) - 1;
    m2epsR = exp(params(6));
    m2epsE = exp(params(7));
    m3epsR = params(8);
    m3epsE = params(9);
    m2fx = exp(params(10));
    m3fx = (params(11));
    b2fx=params(12);
    b3fx=params(13);
     XG(s,:)=[m2etaR,m2etaE,m3etaR,m3etaE,rho,m2epsR,m2epsE,m3epsR,m3epsE,m2fx,m3fx,b2fx,b3fx,Fval];
end 

 XX=sortrows(XG,14);
 
 Est=XX(1,1:end-1)';
 f=XX(1,end);

