function y = mompai(pai,y,g,z,TM)
%y = mompai(pai,y,g,z,TM) computes second moments of pai. Specifically, it computes mean, std, acorr, corr(pai,y), corr(pai,g), and corr(pai,z). 
%Inputs: The stochastic processes pai, y, g, and z, and the transition matrix TM. Note that the process pai is given by a 4x4 matrix because both the previous and current states are relevant in determining this variable. 
%Output: y, a vector containing the moments listed above.
%(c) Stephanie Schmitt-Grohe and Martin Uribe


n=length(g); %number of states

e=ones(n,1);

for j=1:n

Epai(j,1) = mean((pai(:,:,j).*TM)*e);
Ey(j,1) = mean(y(j,:));
Eg(j,1) = mean(g(j,:));
Ez(j,1) = mean(z(j,:));

dpai = pai(:,:,j)-Epai(j);
dy = y(j,:)'-Ey(j);
dg = g(j,:)'-Eg(j);
dz = z(j,:)'-Ez(j);

Vpai(j,1) = mean((dpai.*dpai.*TM)*e); 
Vy(j,1) = mean(dy.*dy); 
Vg(j,1) = mean(dg.*dg); 
Vz(j,1) = mean(dz.*dz); 

ACOVpai(j,1) = 0;
for pa=1:n
for pr=1:n
for fu=1:n
ACOVpai(j,1) = ACOVpai(j) + 1/n * TM(pa,pr) * TM(pr,fu)*dpai(pa,pr) * dpai(pr,fu);
end
end
end
COVpaiy(j,1) = mean(((dpai*diag(dy)).*TM)*e); 
COVpaig(j,1) = mean(((dpai*diag(dg)).*TM)*e);
COVpaiz(j,1) = mean(((dpai*diag(dz)).*TM)*e);

end

EPAI = mean(Epai);

VPAI = mean(Vpai);
VY = mean(Vy);
VG = mean(Vg);
VZ = mean(Vz);

SDPAI = sqrt(VPAI);
SDY = sqrt(VY);
SDG = sqrt(VG);
SDZ = sqrt(VZ);

ACOVPAI = mean(ACOVpai);
COVPAIY = mean(COVpaiy);
COVPAIG = mean(COVpaig);
COVPAIZ = mean(COVpaiz);

ACORPAI = ACOVPAI/VPAI;
CORPAIY = COVPAIY/SDPAI/SDY;
CORPAIG = COVPAIG/SDPAI/SDG;
CORPAIZ = COVPAIZ/SDPAI/SDZ;

y = [EPAI SDPAI ACORPAI CORPAIY CORPAIG CORPAIZ];