Iteration of third order solution

Hi,

I’ve compared a third order simulation obtained with Dynare (oo_.endo_simul) with the simulation obtained iterating the
non-linear state space returned by Dynare (i.e. oo_.g_0 etc…) for a ad hoc “cubic” model.
The results seem to differ if (e.g.) there are two mixed variables (i.e. M_.nboth=2), as opposed to 1 mixed variable (or when the loading – coefficient – of the mixed variable is essentially zero).
Here is the iteration code that I use. Can somebody spot any mistake in what I do? (the model is attached).

%third order folding matrix and weighting vector
S111=zeros(20,64);
w111=zeros(1,20);
S111(1,1)=1;  
S111(2,2)=1;  
S111(3,3)=1;  
S111(4,4)=1;  
S111(5,6)=1;  
S111(6,7)=1;  
S111(7,8)=1;  
S111(8,11)=1; 
S111(9,12)=1; 
S111(10,16)=1;
S111(11,22)=1;
S111(12,23)=1;
S111(13,24)=1;
S111(14,27)=1;
S111(15,28)=1;
S111(16,32)=1;
S111(17,43)=1;
S111(18,44)=1;
S111(19,48)=1;
S111(20,64)=1;
w111(1,1)=1; 
w111(1,2)=3; 
w111(1,3)=3; 
w111(1,4)=3; 
w111(1,5)=3; 
w111(1,6)=6; 
w111(1,7)=6; 
w111(1,8)=3; 
w111(1,9)=6; 
w111(1,10)=3;
w111(1,11)=1;
w111(1,12)=3;
w111(1,13)=3;
w111(1,14)=3;
w111(1,15)=6;
w111(1,16)=3;
w111(1,17)=1;
w111(1,18)=3;
w111(1,19)=3;
w111(1,20)=1;
% second order folding matrix and weighting vector
S11=zeros(10,16);
w11=zeros(1,10);
S11(1,1)=1;  
S11(2,2)=1;  
S11(3,3)=1;  
S11(4,4)=1;  
S11(5,6)=1;  
S11(6,7)=1;  
S11(7,8)=1;  
S11(8,11)=1; 
S11(9,12)=1; 
S11(10,16)=1;
w11(1,1)=1; 
w11(1,2)=2; 
w11(1,3)=2; 
w11(1,4)=2; 
w11(1,5)=1; 
w11(1,6)=2; 
w11(1,7)=2; 
w11(1,8)=1; 
w11(1,9)=2; 
w11(1,10)=1;
% unfold (fold) coefficient matrices (kronecker)
g1=oo_.dr.g_1;
g2=oo_.dr.g_2*diag(w11)*S11;
g3=oo_.dr.g_3*diag(w111)*S111;
shks=oo_.exo_simul';% shock series used also in Dynare
ystart=oo_.steady_state;% steady state !! Note that constants are zero in this model for the parametrization used
zNP=[ystart(M_.nstatic+1:M_.nstatic+nspred,1,end)*0;shks(:,1)];% set initial conditions to zero for all states
per = size(shks,2);
for jj=2:per 
zNP_=g1*zNP+g2*kron(zNP,zNP)+g3*kron(zNP,kron(zNP,zNP));
zNP=[zNP_(M_.nstatic+1:M_.nstatic+nspred,1);shks(:,jj)];
yNP_(oo_.dr.order_var,jj)=zNP_;
end

Then I compare

[oo_.endo_simul(:,1:10)-yNP_(:,2:11)]

where the indexes differ as my recursion does not throw away the first column of zeros.
The discrepancy exists only if I “switch on” (parameter fwd>0) an extra “mixed” variable.

I suspect that the folding of the kronecker products might be wrong, although what I think I have is a Lexicographic non-decreasing ordering.

Hope somebody can help.

Cheers

Gianni
cubic.mod (844 Bytes)

Found the problem: g_0 was missing in the iteration :blush: