function y=moments(x,y,g,z,TM)
%y=moments(x,y,g,z,TM) computes second moments of the variable x. Specifically, it computes mean, std, acorr, corr(x,y), corr(x,g), and corr(x,z). 
%Inputs: The stochastic processes x, y, g, and z, and the transition matrix TM. Note that the stochastic process of each variable is given by just 4 numbers, one for each possible state of the economy. 
%Output: y, a vector containing the moments listed above.
%(c) Stephanie Schmitt-Grohe and Martin Uribe

e=ones(size(x(:,1)));

for j=1:size(x,1)
Ex(j,1) = mean(x(j,:));
Ey(j,1) = mean(y(j,:));
Eg(j,1) = mean(g(j,:));
Ez(j,1) = mean(z(j,:));

dx = x(j,:)'-Ex(j);
dy = y(j,:)'-Ey(j);
dg = g(j,:)'-Eg(j);
dz = z(j,:)'-Ez(j);

Vx(j,1) = mean(dx.*dx); 
Vy(j,1) = mean(dy.*dy); 
Vg(j,1) = mean(dg.*dg); 
Vz(j,1) = mean(dz.*dz); 

ACOVx(j,1) = mean(( dx*dx' .* TM )*e);
COVxy(j,1) = mean(dx.*dy);
COVxg(j,1) = mean(dx.*dg);
COVxz(j,1) = mean(dx.*dz);

end

EX = mean(Ex);

VX = mean(Vx);
VY = mean(Vy);
VG = mean(Vg);
VZ = mean(Vz);

SDX = sqrt(VX);
SDY = sqrt(VY);
SDG = sqrt(VG);
SDZ = sqrt(VZ);

ACOVX = mean(ACOVx);
COVXY = mean(COVxy);
COVXG = mean(COVxg);
COVXZ = mean(COVxz);

ACORX = ACOVX/VX;
CORXY = COVXY/SDX/SDY;
CORXG = COVXG/SDX/SDG;
CORXZ = COVXZ/SDX/SDZ;

y = [EX SDX ACORX CORXY CORXG CORXZ];