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

    maxage=36;

    m2etaR = (x(1));
    m2etaE = (x(2));
    m3etaR = x(3);
    m3etaE =x(4);
    rho     = (x(5));
    m2eps = (x(6));
    m3eps = x(7);
    m2fx = (x(8));
    m3fx = (x(9));
    b2fx=x(10);
    b3fx=x(11);
    
   
    m280s=m2etaR*ones(42,1);
    m290s=m2etaE*ones(11,1);
    m2eta=[m280s;m290s];
    
    m380s=m3etaR*ones(42,1);
    m390s=m3etaE*ones(11,1);
    m3eta=[m380s;m390s];
    
    m2eps=m2eps*ones(iNUM_YEARS,1);
    m3eps=m3eps*ones(iNUM_YEARS,1);
    
   
% MOMENTS BASED ON LEVELS%

    m4etaR = exp(xx(1));
    m4etaE = exp(xx(2));
    m4eps = exp(xx(3));
    m4fx= exp(xx(4));
    b4fx= (xx(5));
    
    m480s=m4etaR*ones(42,1);
    m490s=m4etaE*ones(11,1);
    m4eta=[m480s;m490s];
    
    m4eps=m4eps*ones(iNUM_YEARS,1);
    
% MOMENTS BASED ON LEVELS%
%% Calculate the initial  theoretical moments
PAST_EZ2ht      = zeros(maxage,maxage);
PAST_EZ3ht    = zeros(maxage,maxage);
PAST_EZ4ht    = zeros(maxage,maxage);

PAST_EZ2ht2      = zeros(maxage,maxage);
PAST_EZ3ht2      = zeros(maxage,maxage);
PAST_EZ4ht2      = 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_EZ4ht(t,j) = rho.^(4*(maxage-t));
        
        PAST_EZ2ht2(t,j) = m2eta(t);
        PAST_EZ3ht2(t,j) = m3eta(t);
        PAST_EZ4ht2(t,j) = m4eta(t);
         end
    end
    
PAST_EZ4ht3A      =zeros(1,36);
PAST_EZ4ht3B      =zeros(1,36);
    for j = 1:maxage
         for t = maxage+(1-j):maxage
             PAST_EZ4ht3B(j) =PAST_EZ4ht3B(j)+ (PAST_EZ2ht2(t,j)^2)*rho.^(4*(maxage-t));
                for s = maxage+(1-j):maxage
                       PAST_EZ4ht3A(j) = PAST_EZ4ht3A(j)+PAST_EZ2ht2(t,j)*PAST_EZ2ht2(s,j).*rho.^(2*(maxage-t)+2*(maxage-s));
                end
         end
    end    
    
MATV=PAST_EZ2ht.*PAST_EZ2ht2;
InitialV=sum(MATV,1);
MATS=PAST_EZ3ht.*PAST_EZ3ht2;
InitialS=sum(MATS,1);

InitialKA=   3*(PAST_EZ4ht3A- PAST_EZ4ht3B);
MATK=PAST_EZ4ht.*PAST_EZ4ht2;
InitialKB=sum(MATK,1);
InitialK=InitialKA+InitialKB;

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

EZ2ht(1,:) = InitialV;
EZ3ht(1,:) =InitialS;
EZ4ht(1,:) =InitialK;
    
    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

    for t = 2:iNUM_YEARS
        GGG=zeros(1,maxage);
        EZ4ht(t,1) = m4eta(t+35);
        for j=2:maxage
                for tt = 1:j-1
                       GGG(j)= GGG(j)+m2eta(t+35-tt).*m2eta(t+35).*rho.^(2*tt);
                end   
        EZ4ht(t,j) = rho^4*EZ4ht(t-1,j-1) + m4eta(t+35)+6*GGG(j);
        end
    end

    
%% Calculate theoretical moments for second moments and third moments
    EZhtZhk= zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    EZ2htZhk = zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    Tm4   =zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    EZ2htZ2hk =zeros(iNUM_YEARS*(iNUM_YEARS+1)/2,maxage);
    
    j = 1;
    for t = 1:iNUM_YEARS
      EZhtZhk(j,:)=EZ2ht(t,:);
      EZ2htZhk(j,:)=EZ3ht(t,:);
      
      Tm4(j,:)     = (m4fx+(t+35-[1:36]).*b4fx)+m4eps(t)+EZ4ht(t,:) +6*((m2fx+(t+35-[1:36]).*b2fx)+m2eps(t)).*EZ2ht(t,:) +6*(m2fx+(t+35-[1:36]).*b2fx)*m2eps(t);
      EZ2htZ2hk(j,:)=EZ4ht(t,:);
        j = j + 1;
        for s = t+1:iNUM_YEARS
                EZhtZhk(j,:)= rho^(s-t)*EZ2ht(t,:);

                EZ2htZhk(j,:)= rho^(s-t)*EZ3ht(t,:);
                
                GGG=zeros(1,maxage);    
                for jj=1:maxage
                        for tt = 0:jj-1
                               GGG(jj)= GGG(jj)+m2eta(t+35-tt).*m2eta(s+35).*rho.^(2*tt);
                        end   
                end
                
                EZ2htZ2hk(j,:) = rho^2*EZ2htZ2hk(j-1,:)+GGG;
                Tm4(j,:) = (m4fx+(t+35-[1:36]).*b4fx)+m2eps(t)*m2eps(s)+(m2fx+(t+35-[1:36]).*b2fx)*(m2eps(t)+m2eps(s))+((m2fx+(t+35-[1:36]).*b2fx)).*(EZ2ht(t,:)+EZ2ht(s,:)+4*EZhtZhk(j,:))...
                            + EZ2htZ2hk(j,:) + m2eps(t)*EZ2ht(s,:)  +  m2eps(s)*EZ2ht(t,:);

            j = j + 1;
        end
    end
 
    X = Tm4(:).*Ihtk(:); 
    M = m4htk(:).*Ihtk(:);
    DD=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 = M-X;
    f=f'*f;

end