diff --git a/matlab/ep/ep_problem_0.m b/matlab/ep/ep_problem_0.m
new file mode 100644
index 0000000000000000000000000000000000000000..5ec81b475fd6899f67e666b24aefe01a06c37e2b
--- /dev/null
+++ b/matlab/ep/ep_problem_0.m
@@ -0,0 +1,122 @@
+function [res, A, info] = ep_problem_0(y, x, pfm)
+
+info = false;
+
+params = pfm.params;
+steady_state = pfm.steady_state;
+ny = pfm.ny;
+periods = pfm.periods;
+dynamic_model = pfm.dynamic_model;
+lead_lag_incidence = pfm.lead_lag_incidence;
+i_cols_1 = pfm.i_cols_1;
+i_cols_j = pfm.i_cols_j;
+i_cols_T = pfm.i_cols_T;
+order = pfm.order;
+hybrid_order = pfm.hybrid_order;
+h_correction = pfm.h_correction;
+nodes = pfm.nodes;
+weights = pfm.weights;
+nnodes = pfm.nnodes;
+
+i_cols_p = pfm.i_cols_p;
+i_cols_s = pfm.i_cols_s;
+i_cols_f = pfm.i_cols_f;
+i_rows = pfm.i_rows;
+
+i_cols_Ap = pfm.i_cols_Ap;
+i_cols_As = pfm.i_cols_As;
+i_cols_Af = pfm.i_cols_Af;
+i_hc = pfm.i_hc;
+
+Y = pfm.Y;
+Y(pfm.i_upd_y) = y;
+
+A1 = pfm.A1;
+res = pfm.res;
+
+for i = 1:order+1
+ i_w_p = 1;
+ for j = 1:nnodes^(i-1)
+ innovation = x;
+ if i > 1
+ innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
+ end
+ if i <= order
+ for k=1:nnodes
+ if hybrid_order && i==order
+ z = [Y(i_cols_p,i_w_p);
+ Y(i_cols_s,j);
+ Y(i_cols_f,(j-1)*nnodes+k)+h_correction(i_hc)];
+ else
+ z = [Y(i_cols_p,i_w_p);
+ Y(i_cols_s,j);
+ Y(i_cols_f,(j-1)*nnodes+k)];
+ end
+
+ [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
+ if i == 1
+ % in first period we don't keep track of
+ % predetermined variables
+ i_cols_A = [i_cols_As - ny; i_cols_Af];
+ A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
+ else
+ i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
+ A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
+ end
+ res(:,i,j) = res(:,i,j)+weights(k)*d1;
+ i_cols_Af = i_cols_Af + ny;
+ end
+ else
+ z = [Y(i_cols_p,i_w_p);
+ Y(i_cols_s,j);
+ Y(i_cols_f,j)];
+ [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
+ if i == 1
+ % in first period we don't keep track of
+ % predetermined variables
+ i_cols_A = [i_cols_As - ny; i_cols_Af];
+ A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
+ else
+ i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
+ A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
+ end
+ res(:,i,j) = d1;
+ i_cols_Af = i_cols_Af + ny;
+ end
+ i_rows = i_rows + ny;
+ if mod(j,nnodes) == 0
+ i_w_p = i_w_p + 1;
+ end
+ if i > 1
+ if mod(j,nnodes) == 0
+ i_cols_Ap = i_cols_Ap + ny;
+ end
+ i_cols_As = i_cols_As + ny;
+ end
+ end
+ i_cols_p = i_cols_p + ny;
+ i_cols_s = i_cols_s + ny;
+ i_cols_f = i_cols_f + ny;
+end
+nzA = cell(periods,pfm.world_nbr);
+for j=1:pfm.world_nbr
+ i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
+ offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
+ offset_r = (j-1)*ny;
+ for i=order+2:periods
+ [d1,jacobian] = dynamic_model(Y(i_rows_y,j), x, params, steady_state, i+1);
+ if i == periods
+ [ir,ic,v] = find(jacobian(:,i_cols_T));
+ else
+ [ir,ic,v] = find(jacobian(:,i_cols_j));
+ end
+ nzA{i,j} = [offset_r+ir,offset_c+pfm.icA(ic), v]';
+ res(:,i,j) = d1;
+ i_rows_y = i_rows_y + ny;
+ offset_c = offset_c + pfm.world_nbr*ny;
+ offset_r = offset_r + pfm.world_nbr*ny;
+ end
+end
+A2 = [nzA{:}]';
+A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*pfm.world_nbr,pfm.dimension)];
+res = res(pfm.i_upd_r);
diff --git a/matlab/ep/ep_problem_1.m b/matlab/ep/ep_problem_1.m
index 35757c83af6c05b28e4c52345bf54e3cc5781cf3..9579e7dbf1aca75b9f5bd8c4290b386b6a2138cd 100644
--- a/matlab/ep/ep_problem_1.m
+++ b/matlab/ep/ep_problem_1.m
@@ -1,122 +1,196 @@
-function [res, A, info] = ep_problem_1(y, x, pfm)
+function [res,A,info] = ep_problem_1(y, x, pfm)
info = false;
+A = [];
+dynamic_model = pfm.dynamic_model;
+ny = pfm.ny;
params = pfm.params;
steady_state = pfm.steady_state;
-ny = pfm.ny;
+order = pfm.stochastic_order;
+nodes = pfm.nodes;
+nnodes = pfm.nnodes;
+weights = pfm.weights;
+h_correction = pfm.h_correction;
+dimension = pfm.dimension;
+world_nbr = pfm.world_nbr;
+nnzA = pfm.nnzA;
periods = pfm.periods;
-dynamic_model = pfm.dynamic_model;
-lead_lag_incidence = pfm.lead_lag_incidence;
+i_rows = pfm.i_rows';
+i_cols = pfm.i_cols;
+nyp = pfm.nyp;
+nyf = pfm.nyf;
+hybrid_order = pfm.hybrid_order;
i_cols_1 = pfm.i_cols_1;
i_cols_j = pfm.i_cols_j;
+icA = pfm.icA;
i_cols_T = pfm.i_cols_T;
-order = pfm.order;
-hybrid_order = pfm.hybrid_order;
-h_correction = pfm.h_correction;
-nodes = pfm.nodes;
-weights = pfm.weights;
-nnodes = pfm.nnodes;
+eq_index = pfm.eq_index;
-i_cols_p = pfm.i_cols_p;
-i_cols_s = pfm.i_cols_s;
-i_cols_f = pfm.i_cols_f;
-i_rows = pfm.i_rows;
-i_cols_Ap = pfm.i_cols_Ap;
-i_cols_As = pfm.i_cols_As;
-i_cols_Af = pfm.i_cols_Af;
-i_hc = pfm.i_hc;
-Y = pfm.Y;
-Y(pfm.i_upd_y) = y;
+i_cols_p = i_cols(1:nyp);
+i_cols_s = i_cols(nyp+(1:ny));
+i_cols_f = i_cols(nyp+ny+(1:nyf));
+i_cols_Ap0 = i_cols_p;
+i_cols_As = i_cols_s;
+i_cols_Af0 = i_cols_f - ny;
+i_hc = i_cols_f - 2*ny;
-A1 = pfm.A1;
-res = pfm.res;
+nzA = cell(periods,world_nbr);
+res = zeros(ny,periods,world_nbr);
+Y = zeros(ny*(periods+2),world_nbr);
+Y(1:ny,1) = pfm.y0;
+Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr);
+Y(pfm.i_upd_y) = y;
+offset_r0 = 0;
for i = 1:order+1
i_w_p = 1;
- for j = 1:nnodes^(i-1)
+ for j = 1:(1+(nnodes-1)*(i-1))
innovation = x;
- if i > 1
- innovation(i+1,:) = nodes(mod(j-1,nnodes)+1,:);
- end
- if i <= order
+ if i <= order && j == 1
+ % first world, integrating future shocks
+ if nargout > 1
+ A1 = sparse([],[],[],i*(1+(nnodes-1)*(i-1))*ny,dimension,nnzA*world_nbr);
+ end
for k=1:nnodes
- if hybrid_order && i==order
- z = [Y(i_cols_p,i_w_p);
- Y(i_cols_s,j);
- Y(i_cols_f,(j-1)*nnodes+k)+h_correction(i_hc)];
+ if nargout > 1
+ if i == 2
+ i_cols_Ap = i_cols_Ap0;
+ elseif i > 2
+ i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes- ...
+ 1)*(i-2)*(i-3)/2);
+ end
+ if k == 1
+ i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ...
+ 2);
+ else
+ i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ...
+ 2+(i-1)*(nnodes-1) ...
+ + k - 1);
+ end
+ end
+ if i > 1
+ innovation(i+1,:) = nodes(k,:);
+ end
+ if k == 1
+ k1 = 1;
else
- z = [Y(i_cols_p,i_w_p);
- Y(i_cols_s,j);
- Y(i_cols_f,(j-1)*nnodes+k)];
+ k1 = (nnodes-1)*(i-1)+k;
end
-
- [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
- if i == 1
- % in first period we don't keep track of
- % predetermined variables
- i_cols_A = [i_cols_As - ny; i_cols_Af];
- A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_1);
+ if hybrid_order && (k > 1 || i == order)
+ z = [Y(i_cols_p,1);
+ Y(i_cols_s,1);
+ Y(i_cols_f,k1)+h_correction(i_hc)];
+ else
+ z = [Y(i_cols_p,1);
+ Y(i_cols_s,1);
+ Y(i_cols_f,k1)];
+ end
+ if nargout > 1
+ [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
+ if i == 1
+ % in first period we don't keep track of
+ % predetermined variables
+ i_cols_A = [i_cols_As - ny; i_cols_Af];
+ A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(eq_index,i_cols_1);
+ else
+ i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
+ A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(eq_index,i_cols_j);
+ end
else
- i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
- A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(:,i_cols_j);
+ d1 = dynamic_model(z,innovation,params,steady_state,i+1);
end
- res(:,i,j) = res(:,i,j)+weights(k)*d1;
+ res(:,i,1) = res(:,i,1)+weights(k)*d1(eq_index);
+ end
+ if nargout > 1
+ [ir,ic,v] = find(A1);
+ nzA{i,j} = [ir,ic,v]';
+ end
+ elseif j > 1 + (nnodes-1)*(i-2)
+ % new world, using previous state of world 1 and hit
+ % by shocks from nodes
+ if nargout > 1
+ i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2);
+ i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2);
+ end
+ k = j - (nnodes-1)*(i-2);
+ innovation(i+1,:) = nodes(k,:);
+ z = [Y(i_cols_p,1);
+ Y(i_cols_s,j);
+ Y(i_cols_f,j)];
+ if nargout > 1
+ [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
+ i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
+ [ir,ic,v] = find(jacobian(eq_index,i_cols_j));
+ nzA{i,j} = [i_rows(ir),i_cols_A(ic), v]';
+ else
+ d1 = dynamic_model(z,innovation,params,steady_state,i+1);
+ end
+ res(:,i,j) = d1(eq_index);
+ if nargout > 1
i_cols_Af = i_cols_Af + ny;
end
else
- z = [Y(i_cols_p,i_w_p);
+ % existing worlds other than 1
+ if nargout > 1
+ i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1);
+ i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2);
+ end
+ z = [Y(i_cols_p,j);
Y(i_cols_s,j);
Y(i_cols_f,j)];
- [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
- if i == 1
- % in first period we don't keep track of
- % predetermined variables
- i_cols_A = [i_cols_As - ny; i_cols_Af];
- A1(i_rows,i_cols_A) = jacobian(:,i_cols_1);
- else
+ if nargout > 1
+ [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
- A1(i_rows,i_cols_A) = jacobian(:,i_cols_j);
+ [ir,ic,v] = find(jacobian(eq_index,i_cols_j));
+ nzA{i,j} = [i_rows(ir),i_cols_A(ic),v]';
+ i_cols_Af = i_cols_Af + ny;
+ else
+ d1 = dynamic_model(z,innovation,params,steady_state,i+1);
end
- res(:,i,j) = d1;
- i_cols_Af = i_cols_Af + ny;
+ res(:,i,j) = d1(eq_index);
end
i_rows = i_rows + ny;
if mod(j,nnodes) == 0
i_w_p = i_w_p + 1;
end
- if i > 1
- if mod(j,nnodes) == 0
- i_cols_Ap = i_cols_Ap + ny;
- end
+ if nargout > 1 && i > 1
i_cols_As = i_cols_As + ny;
end
+ offset_r0 = offset_r0 + ny;
end
i_cols_p = i_cols_p + ny;
i_cols_s = i_cols_s + ny;
i_cols_f = i_cols_f + ny;
end
-nzA = cell(periods,pfm.world_nbr);
-for j=1:pfm.world_nbr
- i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
- offset_c = ny*(sum(nnodes.^(0:order-1),2)+j-1);
- offset_r = (j-1)*ny;
+for j=1:world_nbr
+ i_rows_y = i_cols+(order+1)*ny;
+ offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1);
+ offset_r = offset_r0+(j-1)*ny;
for i=order+2:periods
- [d1,jacobian] = dynamic_model(Y(i_rows_y,j), x, params, steady_state, i+1);
- if i == periods
- [ir,ic,v] = find(jacobian(:,i_cols_T));
+ if nargout > 1
+ [d1,jacobian] = dynamic_model(Y(i_rows_y,j),x,params, ...
+ steady_state,i+1);
+ if i < periods
+ [ir,ic,v] = find(jacobian(eq_index,i_cols_j));
+ else
+ [ir,ic,v] = find(jacobian(eq_index,i_cols_T));
+ end
+ nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]';
else
- [ir,ic,v] = find(jacobian(:,i_cols_j));
+ d1 = dynamic_model(Y(i_rows_y,j),x,params, ...
+ steady_state,i+1);
end
- nzA{i,j} = [offset_r+ir,offset_c+pfm.icA(ic), v]';
- res(:,i,j) = d1;
+ res(:,i,j) = d1(eq_index);
i_rows_y = i_rows_y + ny;
- offset_c = offset_c + pfm.world_nbr*ny;
- offset_r = offset_r + pfm.world_nbr*ny;
+ offset_c = offset_c + world_nbr*ny;
+ offset_r = offset_r + world_nbr*ny;
end
end
-A2 = [nzA{:}]';
-A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*pfm.world_nbr,pfm.dimension)];
+if nargout > 1
+ iA = [nzA{:}]';
+ A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension);
+end
res = res(pfm.i_upd_r);
diff --git a/matlab/ep/ep_problem_2.m b/matlab/ep/ep_problem_2.m
deleted file mode 100644
index 26ba78895cc1de7b74448d0c5cb0e45b94422ca3..0000000000000000000000000000000000000000
--- a/matlab/ep/ep_problem_2.m
+++ /dev/null
@@ -1,196 +0,0 @@
-function [res,A,info] = ep_problem_2(y, x, pfm)
-
-info = false;
-A = [];
-
-dynamic_model = pfm.dynamic_model;
-ny = pfm.ny;
-params = pfm.params;
-steady_state = pfm.steady_state;
-order = pfm.stochastic_order;
-nodes = pfm.nodes;
-nnodes = pfm.nnodes;
-weights = pfm.weights;
-h_correction = pfm.h_correction;
-dimension = pfm.dimension;
-world_nbr = pfm.world_nbr;
-nnzA = pfm.nnzA;
-periods = pfm.periods;
-i_rows = pfm.i_rows';
-i_cols = pfm.i_cols;
-nyp = pfm.nyp;
-nyf = pfm.nyf;
-hybrid_order = pfm.hybrid_order;
-i_cols_1 = pfm.i_cols_1;
-i_cols_j = pfm.i_cols_j;
-icA = pfm.icA;
-i_cols_T = pfm.i_cols_T;
-eq_index = pfm.eq_index;
-
-
-
-i_cols_p = i_cols(1:nyp);
-i_cols_s = i_cols(nyp+(1:ny));
-i_cols_f = i_cols(nyp+ny+(1:nyf));
-i_cols_Ap0 = i_cols_p;
-i_cols_As = i_cols_s;
-i_cols_Af0 = i_cols_f - ny;
-i_hc = i_cols_f - 2*ny;
-
-nzA = cell(periods,world_nbr);
-res = zeros(ny,periods,world_nbr);
-Y = zeros(ny*(periods+2),world_nbr);
-Y(1:ny,1) = pfm.y0;
-Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr);
-Y(pfm.i_upd_y) = y;
-
-offset_r0 = 0;
-for i = 1:order+1
- i_w_p = 1;
- for j = 1:(1+(nnodes-1)*(i-1))
- innovation = x;
- if i <= order && j == 1
- % first world, integrating future shocks
- if nargout > 1
- A1 = sparse([],[],[],i*(1+(nnodes-1)*(i-1))*ny,dimension,nnzA*world_nbr);
- end
- for k=1:nnodes
- if nargout > 1
- if i == 2
- i_cols_Ap = i_cols_Ap0;
- elseif i > 2
- i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes- ...
- 1)*(i-2)*(i-3)/2);
- end
- if k == 1
- i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ...
- 2);
- else
- i_cols_Af = i_cols_Af0 + ny*(i-1+(nnodes-1)*i*(i-1)/ ...
- 2+(i-1)*(nnodes-1) ...
- + k - 1);
- end
- end
- if i > 1
- innovation(i+1,:) = nodes(k,:);
- end
- if k == 1
- k1 = 1;
- else
- k1 = (nnodes-1)*(i-1)+k;
- end
- if hybrid_order && (k > 1 || i == order)
- z = [Y(i_cols_p,1);
- Y(i_cols_s,1);
- Y(i_cols_f,k1)+h_correction(i_hc)];
- else
- z = [Y(i_cols_p,1);
- Y(i_cols_s,1);
- Y(i_cols_f,k1)];
- end
- if nargout > 1
- [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
- if i == 1
- % in first period we don't keep track of
- % predetermined variables
- i_cols_A = [i_cols_As - ny; i_cols_Af];
- A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(eq_index,i_cols_1);
- else
- i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
- A1(i_rows,i_cols_A) = A1(i_rows,i_cols_A) + weights(k)*jacobian(eq_index,i_cols_j);
- end
- else
- d1 = dynamic_model(z,innovation,params,steady_state,i+1);
- end
- res(:,i,1) = res(:,i,1)+weights(k)*d1(eq_index);
- end
- if nargout > 1
- [ir,ic,v] = find(A1);
- nzA{i,j} = [ir,ic,v]';
- end
- elseif j > 1 + (nnodes-1)*(i-2)
- % new world, using previous state of world 1 and hit
- % by shocks from nodes
- if nargout > 1
- i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2);
- i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2);
- end
- k = j - (nnodes-1)*(i-2);
- innovation(i+1,:) = nodes(k,:);
- z = [Y(i_cols_p,1);
- Y(i_cols_s,j);
- Y(i_cols_f,j)];
- if nargout > 1
- [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
- i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
- [ir,ic,v] = find(jacobian(eq_index,i_cols_j));
- nzA{i,j} = [i_rows(ir),i_cols_A(ic), v]';
- else
- d1 = dynamic_model(z,innovation,params,steady_state,i+1);
- end
- res(:,i,j) = d1(eq_index);
- if nargout > 1
- i_cols_Af = i_cols_Af + ny;
- end
- else
- % existing worlds other than 1
- if nargout > 1
- i_cols_Ap = i_cols_Ap0 + ny*(i-2+(nnodes-1)*(i-2)*(i-3)/2+j-1);
- i_cols_Af = i_cols_Af0 + ny*(i+(nnodes-1)*i*(i-1)/2+j-2);
- end
- z = [Y(i_cols_p,j);
- Y(i_cols_s,j);
- Y(i_cols_f,j)];
- if nargout > 1
- [d1,jacobian] = dynamic_model(z,innovation,params,steady_state,i+1);
- i_cols_A = [i_cols_Ap; i_cols_As; i_cols_Af];
- [ir,ic,v] = find(jacobian(eq_index,i_cols_j));
- nzA{i,j} = [i_rows(ir),i_cols_A(ic),v]';
- i_cols_Af = i_cols_Af + ny;
- else
- d1 = dynamic_model(z,innovation,params,steady_state,i+1);
- end
- res(:,i,j) = d1(eq_index);
- end
- i_rows = i_rows + ny;
- if mod(j,nnodes) == 0
- i_w_p = i_w_p + 1;
- end
- if nargout > 1 && i > 1
- i_cols_As = i_cols_As + ny;
- end
- offset_r0 = offset_r0 + ny;
- end
- i_cols_p = i_cols_p + ny;
- i_cols_s = i_cols_s + ny;
- i_cols_f = i_cols_f + ny;
-end
-for j=1:world_nbr
- i_rows_y = i_cols+(order+1)*ny;
- offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1);
- offset_r = offset_r0+(j-1)*ny;
- for i=order+2:periods
- if nargout > 1
- [d1,jacobian] = dynamic_model(Y(i_rows_y,j),x,params, ...
- steady_state,i+1);
- if i < periods
- [ir,ic,v] = find(jacobian(eq_index,i_cols_j));
- else
- [ir,ic,v] = find(jacobian(eq_index,i_cols_T));
- end
- nzA{i,j} = [offset_r+ir,offset_c+icA(ic), v]';
- else
- d1 = dynamic_model(Y(i_rows_y,j),x,params, ...
- steady_state,i+1);
- end
- res(:,i,j) = d1(eq_index);
- i_rows_y = i_rows_y + ny;
- offset_c = offset_c + world_nbr*ny;
- offset_r = offset_r + world_nbr*ny;
- end
-end
-if nargout > 1
- iA = [nzA{:}]';
- A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension);
-end
-res = res(pfm.i_upd_r);
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
index 9196ea79bceff70a92f1d866823c964cafefa689..2773b527313f2006f036fd8d4a9d6c175ced3786 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
@@ -166,6 +166,6 @@ if update_options_struct
end
end
-[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_1, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
+[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_0, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
endo_simul(:,2) = y(1:ny);
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
index d9f0deda766f7a5045265b7a0140de804e13b9ec..0958d28d88607961c0c159b086ff933c0827d872 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
@@ -138,5 +138,5 @@ end
pfm.y0 = endo_simul(:,1);
-[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_2, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
+[y, errorflag, ~, ~, errorcode] = dynare_solve(@ep_problem_1, y, options_.simul.maxit, options_.dynatol.f, options_.dynatol.x, options_, exo_simul, pfm);
endo_simul(:,2) = y(1:pfm.ny);