function [sigxJ,sigx2J]=mom(MX,pai,varshock,J)
%[sigxJ,sigx2J]=mom(MX,pai,varshock,J)
% Computes E{x(t)*x(t+J)'}
% where x(t) evolves as
% x2(t+1) = pai x2(t) + e(t+1)
% x(t) = MX x2(t)
%where Ee(t)e(t)'=varshock
% (c) Stephanie Schmitt-Grohe and Martin Uribe
% September 19, 2001


if nargin<4
J=0;
end

[X,D]=eig(pai);
if rank(X)<size(X,1)

%Doubling algorithm
pai_old=pai;
sig_old=varshock;
sigx2_old=eye(size(pai));
diferenz=.1;
while diferenz>1e-25;
sigx2=pai_old*sigx2_old*pai_old'+sig_old;
diferenz = max(max(abs(sigx2-sigx2_old)));
sig_old=pai_old*sig_old*pai_old'+sig_old;
pai_old=pai_old*pai_old;
sigx2_old=sigx2;
end    %while diferenz

else

%Algebraic method
%Get the variance of inv(X)*x2
vare=X\varshock *(X\eye(size(X)))';

for i=1:size(pai,1)
for j=1:size(pai,1)
sig(i,j) = 1/(1-D(i,i)*D(j,j)')*vare(i,j);
end   %for j
end   %for i

%Get the variance of x2
sigx2=X*sig*X';

end   %if rank


%Get E{x2(t)*x2(t+J)'}
sigx2J=pai^(-min(0,J))*sigx2*(pai')^(max(0,J));

%Get E{x(t)*x(t+J)'}
sigxJ=real(MX*sigx2J*MX');
sigx2J=real(sigx2J);