function [m,S]=moments(yB,yA,y1,y2,y3,y4,xB,xA,x1,x2,x3,x4,theta,should_scale)

% if should_scale ~= 0
%     fprintf('should_scale should be 0')
% end

	
    z1=[yA yB x1]*theta;
    z2=[y1 yA x2]*theta;
    z3=[y2 y1 x3]*theta;
    z4=[y3 y2 x4]*theta;
    
    xb1=x1*theta(3:length(theta));
    xb2=x2*theta(3:length(theta));
    xb3=x3*theta(3:length(theta));
    xb4=x4*theta(3:length(theta));
    
    gamma1=theta(1);
	gamma2=theta(2);
	
	ooo=ones(length(y1),1);
    
    d1=(1-y1);
    d2=(1-y2);
    d3=(1-y3);
    d4=(1-y4);


    a1=exp(xb2+gamma2*yA-xb3)-exp(xb4-xb3+gamma1);
	a2=exp(xb2+gamma2*yA-xb4-gamma1)-1;
	a3=-ooo;
	a4=exp(xb4-xb1-gamma1*yA-gamma2*yB+gamma1);
	a5=exp(xb4-xb1-gamma1*yA-gamma2*yB+gamma1).*(1+exp(xb2+gamma1+gamma2*yA-xb3-gamma2)-exp(xb4-xb3+gamma1-gamma2));
	a6=exp(xb2+gamma1+gamma2*yA-xb1-gamma1*yA-gamma2*yB);
	
    ma=d1.*d2.*y3.*d4.*a1 ...
        +d1.*d2.*y3.*y4.*a2 ...
        +d1.*y2.*a3...
        +y1.*d2.*d3.*a4 ...
        +y1.*d2.*y3.*d4.*a5 ...
        +y1.*d2.*y3.*y4.*a6 ;
    if should_scale==1
        scal=abs(a1)+abs(a2)+abs(a3)+abs(a4)+abs(a5)+abs(a6);
        ma=ma./scal;
    end

	a1=exp(xb1+gamma1*yA+gamma2*yB-xb2-gamma2*yA);
	a2=exp(xb1+gamma1*yA+gamma2*yB-xb4-gamma2).*(1+exp(xb3+gamma1-xb2-gamma2*yA)-exp(xb3+gamma1-xb4-gamma2));
	a3=exp(xb1+gamma1*yA+gamma2*yB-xb4-gamma2);
	a4=-ooo;
	a5=exp(xb4+gamma2-xb2-gamma2*yA-gamma1)-1;
	a6=exp(xb3+gamma1+gamma2-xb2-gamma1-gamma2*yA)-exp(xb3-xb4+gamma1);
    mb=d1.*y2.*d3.*d4.*a1...
        +d1.*y2.*d3.*y4.*a2...
        +d1.*y2.*y3.*a3...
        +y1.*d2.*a4...
        +y1.*y2.*d3.*d4.*a5...
        +y1.*y2.*d3.*y4.*a6;
     if should_scale==1
        scal=abs(a1)+abs(a2)+abs(a3)+abs(a4)+abs(a5)+abs(a6);
        mb=mb./scal;
     end
     
    
    
	a1=(exp(xb2+gamma2*yA-xb4)-1).*(1-exp(xb3-xb4));
	a2=exp(xb2+gamma2*yA-xb4)-1;
	a3=-ooo;
	a4=exp(xb4-xb1-gamma1*yA-gamma2*yB);
	a5=exp(xb2+gamma2*yA-xb1-gamma1*yA-gamma2*yB).*(exp(gamma1)+exp(xb3-xb2-gamma2*yA+gamma2)-exp(xb3-xb4+gamma1+gamma2)); % CHECH
	a6=exp(xb2+gamma2*yA-xb1-gamma1*yA-gamma2*yB+gamma1);
    mc=d1.*d2.*d3.*y4.*a1...
        +d1.*d2.*y3.*a2...
        +d1.*y2.*a3...
        +y1.*d2.*d3.*d4.*a4...
        +y1.*d2.*d3.*y4.*a5...
        +y1.*d2.*y3.*a6;
     if should_scale==1
        scal=abs(a1)+abs(a2)+abs(a3)+abs(a4)+abs(a5)+abs(a6);
        mc=mc./scal;
     end
    
	a1=exp(xb1+gamma1*yA+gamma2*yB-xb2-gamma2*yA);
	a2=exp(xb1+gamma1*yA+gamma2*yB-xb2-gamma2*yA).*(1+exp(xb2+gamma2*yA-xb3-gamma1)-exp(xb4-xb3+gamma2));
	a3=exp(xb1+gamma1*yA+gamma2*yB-xb4-gamma1-gamma2);
	a4=-ooo;
	a5=exp(xb4-xb2-gamma2*yA+gamma2)-1;
	a6=(exp(xb4-xb2-gamma2*yA+gamma2)-1).*(1-exp(xb4-xb3));
    md=d1.*y2.*d3.*a1...
        +d1.*y2.*y3.*d4.*a2...
        +d1.*y2.*y3.*y4.*a3...
        +y1.*d2.*a4 ...
        +y1.*y2.*d3.*a5...
        +y1.*y2.*y3.*d4.*a6;
     if should_scale==1
        scal=abs(a1)+abs(a2)+abs(a3)+abs(a4)+abs(a5)+abs(a6);
        md=md./scal;
     end
    

    xx=[yA.*yB (1-yA).*yB yA.*(1-yB) (1-yA).*(1-yB) x2-x1 x3-x2 x4-x3];
    kk=size(xx,2);
    m0a=repmat(ma,1,kk);
    m0b=repmat(mb,1,kk);
    m0c=repmat(mc,1,kk);
    m0d=repmat(md,1,kk);
    
    mm=[m0a.*xx m0b.*xx m0c.*xx m0d.*xx];


    m=mean(mm);

        
    if nargout>1
        S=cov(mm);
%         S=mm'*mm/length(y1);
    end
 

end
    
    
            
            
                
    
