function [MOMENTS, SE]=my_stat(data) T=size(data,1); m1=mean(data); m2=var(data); m3=mean((data-m1).^3); m4=mean((data-m1).^4); MEAN=m1; STD=sqrt(m2); SKEW=m3/m2^(3/2); KURT=m4/m2^2; h1=data-m1; h2=h1.^2-m2; h3=h1.^3-m3; h4=h1.^4-m4; h=[h1 h2 h3 h4]; T=length(h); R=h'*h/T; n_moving=5; for i=1:n_moving R_temp=h(i+1:T,:)'*h(1:T-i,:)/T; R=R+(R_temp'+R_temp)*(1-i/(n_moving+1)); end W=inv(R); D=[-1 0 0 0; 0 -1 0 0; 3*m2 0 -1 0; 4*m3 0 0 -1]; COV=inv(D'*W*D); SE=sqrt(diag(COV)/T); D2=1/2/sqrt(m2); C2=COV(2,2); SE(2)=sqrt(D2*C2*D2/T); D23=[-1.5*m3/m2^(5/2) 1/m2^(3/2)]; C23=COV(2:3,2:3); SE(3)=sqrt(D23*C23*D23'/T); D24=[-2*m4/m2^3 1/m2^2]; C24=COV([2 4],[2 4]); SE(4)=sqrt(D24*C24*D24'/T); MOMENTS=[MEAN STD SKEW KURT]';