function [f] = MDE_Timevarying_HigherOrderMom_boot6(x,iNUM_YEARS,m2htk,m3htk,m4htk,Ihtk,Contraction_index,NumObs)

    maxage=36;
    
    m2etaR = exp(x(1));
    m2etaE = exp(x(2));
    m3etaR = x(3);
    m3etaE =x(4);
    rho     = 2*exp(x(5))/(1+exp(x(5))) - 1;
    m2eps = exp(x(6));
    m3eps = x(7);
    m2fx = exp(x(8));
    m3fx = (x(9));
    b2fx=x(10);
    b3fx=x(11);
    
    m280s=m2etaR*ones(36,1);
    m290s=m2etaE*ones(17,1);
    m2eta=[m280s;m290s];
    
    m380s=m3etaR*ones(36,1);
    m390s=m3etaE*ones(17,1);
    m3eta=[m380s;m390s];
    
    m2eps=m2eps*ones(iNUM_YEARS,1);
    m3eps=m3eps*ones(iNUM_YEARS,1);
    
% MOMENTS BASED ON LEVELS%
%% Calculate the initial  theoretical moments
PAST_EZ2ht      = zeros(maxage,maxage);
PAST_EZ3ht   = zeros(maxage,maxage);

PAST_EZ2ht2      = zeros(maxage,maxage);
PAST_EZ3ht2      = zeros(maxage,maxage);


    for j = 1:maxage
         for t = maxage+(1-j):maxage
        PAST_EZ2ht(t,j) = rho.^(2*(maxage-t));
        PAST_EZ3ht(t,j) = rho.^(3*(maxage-t));
        
        PAST_EZ2ht2(t,j) = m2eta(t);
        PAST_EZ3ht2(t,j) = m3eta(t);
         end
    end 

   
MATV=PAST_EZ2ht.*PAST_EZ2ht2;
InitialV=sum(MATV,1);
MATS=PAST_EZ3ht.*PAST_EZ3ht2;
InitialS=sum(MATS,1);

EZ2ht      =zeros(iNUM_YEARS,maxage);
EZ3ht      = zeros(iNUM_YEARS,maxage);

EZ2ht(1,:) = InitialV;
EZ3ht(1,:) =InitialS;
    
    for t = 2:iNUM_YEARS
        EZ2ht(t,1) = m2eta(t+35);
        EZ2ht(t,2:end) = rho^2*EZ2ht(t-1,1:end-1) + m2eta(t+35);
        EZ3ht(t,1) = m3eta(t+35);
        EZ3ht(t,2:end) = rho^3*EZ3ht(t-1,1:end-1) + m3eta(t+35);
    end
    
%% Calculate theoretical moments for second moments and third moments
    Tm2 = zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    EZhtZhk= zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    Tm3    = zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    COSKEW = zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    
    j = 1;
    for t = 1:iNUM_YEARS
          Tm2(j,:)= EZ2ht(t,:) + m2eps(t)+(t+35-[1:36])*b2fx;
          EZhtZhk(j,:)=EZ2ht(t,:);
          Tm3(j,:)     = EZ3ht(t,:) + m3eps(t)+(t+35-[1:36])*b3fx;
          COSKEW(j,:)=EZ3ht(t,:);
          
        j = j + 1;
        for s = t+1:iNUM_YEARS
                Tm2(j,:) = rho^(s-t)*EZ2ht(t,:)+(t+35-[1:36])*b2fx;
                EZhtZhk(j,:)= rho^(s-t)*EZ2ht(t,:);
                
                Tm3(j,:) = rho^(s-t)*EZ3ht(t,:)+(t+35-[1:36])*b3fx;
                COSKEW(j,:)= rho^(s-t)*EZ3ht(t,:);
            j = j + 1;
        end
    end
    Tm2 = Tm2 + m2fx;   % fixed effects    
    Tm3 = Tm3 + m3fx;   % fixed effects  
        
    X2 = Tm2(:).*Ihtk(:);  
    X3 = Tm3(:).*Ihtk(:); 
    M2 = m2htk(:).*Ihtk(:);
    M3 = m3htk(:).*Ihtk(:);
    M=[M2;M3];
    X=[X2;X3];
    DD=[Ihtk(:);Ihtk(:)];
    ii = find(not(isnan(M)));
    M  = M(ii);
    X  = X(ii);
    DD=DD(ii);
    iii = find(DD>NumObs);
    M  = M(iii);
    X  = X(iii);
    
    
    f = 100*(M-X);

end