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]';
