function [sigyJ,sigxJ]=mom(gx,hx,varshock,J, method) %[sigyJ,sigxJ]=mom(gx,hx,varshock,J, method) % Computes the unconditional variance-covariance matrix of x(t) with x(t+J), that is sigxJ=E[x(t)*x(t+J)'], %and the unconditional variance covariaance matrix of y(t) with y(t+J), that is sigyJ=E[y(t)*y(t+J)'] % where x(t) evolves as % x(t+1) = hx x(t) + e(t+1) %and y(t) evolves according to % y(t) = gx x(t) %where Ee(t)e(t)'=varshock %The parameter J can be any integer %method =1 : use doubling algorithm %method neq 1 : use algebraic method %(c) Stephanie Schmitt-Grohe and Martin Uribe, April 18, 1990, renewed January 24, 2000 and August 18, 2001. if nargin<4 J=0; end if nargin<5 method =1; end if method == 1 %disp('method=doubling') %Doubling algorithm hx_old=hx; sig_old=varshock; sigx_old=eye(size(hx)); diferenz=.1; while diferenz>1e-25; sigx=hx_old*sigx_old*hx_old'+sig_old; diferenz = max(max(abs(sigx-sigx_old))); sig_old=hx_old*sig_old*hx_old'+sig_old; hx_old=hx_old*hx_old; sigx_old=sigx; end %while diferenz else %Algebraic method %Get the variance of x disp('method=kronecker') sigx = zeros(size(hx)); F = kron(hx,hx); sigx(:) = (eye(size(F))-F)\varshock(:); end %if method %Get E{x(t)*x(t+J)'} sigxJ=hx^(-min(0,J))*sigx*(hx')^(max(0,J)); %Get E{y(t)*y(t+J)'} sigyJ=real(gx*sigxJ*gx'); sigxJ=real(sigxJ);