Skip to content
Snippets Groups Projects
Verified Commit 14b68ee3 authored by Stéphane Adjemian's avatar Stéphane Adjemian
Browse files

Rename routines.

parent 5d430700
No related branches found
No related tags found
No related merge requests found
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);
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);
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);
......@@ -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);
......@@ -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);
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment