function [m2htk,m3htk,m4htk,Ihtk]=empirical_moments_double_sq(AGEi,Ui,CLEAN,iNUM_YEARS,num_entries);


lifecycle= 36;
iota     = ones(iNUM_YEARS,num_entries);
i        = (Ui == -99999); 
iota(i)  = 0;
Ui(i)    = 0;
AGEi(i)  = 0;


AIhtk  = zeros(iNUM_YEARS,iNUM_YEARS);
Am2htk = zeros(iNUM_YEARS,iNUM_YEARS);
Am3htk = zeros(iNUM_YEARS,iNUM_YEARS);
Am4htk = zeros(iNUM_YEARS,iNUM_YEARS);

H     = iNUM_YEARS*(iNUM_YEARS+1)/2;
Ihtk  = zeros(H,lifecycle);
m2htk = zeros(H,lifecycle);
m3htk = zeros(H,lifecycle);
m4htk = zeros(H,lifecycle);

for age = 1:36
  for t = 1:iNUM_YEARS
      log_obs  = find(AGEi(t,:) >= age+24-2 & AGEi(t,:) <=age+24+2);
      jj = length(log_obs);
      if jj > 0
        Uhatbart = sum(Ui(t,log_obs))./sum(iota(t,log_obs));
      else
          continue
      end

      for s = t:iNUM_YEARS
          Ihts=sum(iota(s,log_obs));
          if Ihts >0
          Uhatbartplusk = sum(Ui(s,log_obs))./Ihts;
          m2htki =  (Ui(t,log_obs)-Uhatbart)    .*iota(t,log_obs).*(Ui(s,log_obs) -Uhatbartplusk)    .*iota(s,log_obs);
          m3htki = ((Ui(t,log_obs)-Uhatbart).^2).*iota(t,log_obs).*(Ui(s,log_obs) -Uhatbartplusk)    .*iota(s,log_obs);
          m4htki = ((Ui(t,log_obs)-Uhatbart).^2).*iota(t,log_obs).*((Ui(s,log_obs)-Uhatbartplusk).^2).*iota(s,log_obs);
          log_obs_ij   = iota(t,log_obs) .* iota(s,log_obs);
          AIhtk(t,s)   = sum(log_obs_ij);
          Am2htk(t,s)  = sum(m2htki) / AIhtk(t,s);
          Am3htk(t,s)  = sum(m3htki) / AIhtk(t,s);
          Am4htk(t,s)  = sum(m4htki) / AIhtk(t,s);
          else
               continue
          end
      end
  end
  
          Ihtk(:,age)   = vector_upper_tr(AIhtk);
          m2htk(:,age)  = vector_upper_tr(Am2htk);
          m3htk(:,age)  = vector_upper_tr(Am3htk);
          m4htk(:,age)  = vector_upper_tr(Am4htk);  
end

m2htk=m2htk.*CLEAN;
m3htk=m3htk.*CLEAN;
m4htk=m4htk.*CLEAN;
Ihtk =Ihtk.*CLEAN;