make extended path algorithm 1 as a self contained problem usable by dynare_solve
parent
c3efb214ef
commit
26f2b301b0
|
@ -0,0 +1,194 @@
|
||||||
|
function [res,A,info] = ep_problem_2(y,x,pm)
|
||||||
|
|
||||||
|
info = 0;
|
||||||
|
res = [];
|
||||||
|
A = [];
|
||||||
|
|
||||||
|
dynamic_model = pm.dynamic_model;
|
||||||
|
ny = pm.ny;
|
||||||
|
params = pm.params;
|
||||||
|
steady_state = pm.steady_state;
|
||||||
|
order = pm.order;
|
||||||
|
nodes = pm.nodes;
|
||||||
|
nnodes = pm.nnodes;
|
||||||
|
weights = pm.weights;
|
||||||
|
h_correction = pm.h_correction;
|
||||||
|
dimension = pm.dimension;
|
||||||
|
world_nbr = pm.world_nbr;
|
||||||
|
nnzA = pm.nnzA;
|
||||||
|
periods = pm.periods;
|
||||||
|
i_rows = pm.i_rows';
|
||||||
|
i_cols = pm.i_cols;
|
||||||
|
nyp = pm.nyp;
|
||||||
|
nyf = pm.nyf;
|
||||||
|
hybrid_order = pm.hybrid_order;
|
||||||
|
i_cols_1 = pm.i_cols_1;
|
||||||
|
i_cols_j = pm.i_cols_j;
|
||||||
|
icA = pm.icA;
|
||||||
|
i_cols_T = pm.i_cols_T;
|
||||||
|
|
||||||
|
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_A = i_cols;
|
||||||
|
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) = pm.y0;
|
||||||
|
Y(end-ny+1:end,:) = repmat(steady_state,1,world_nbr);
|
||||||
|
Y(pm.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 nargin > 1
|
||||||
|
A1 = sparse([],[],[],i*ny,dimension,nnzA*world_nbr);
|
||||||
|
end
|
||||||
|
for k=1:nnodes
|
||||||
|
if nargin > 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 == 2 && (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 nargin > 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(:,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
|
||||||
|
else
|
||||||
|
d1 = dynamic_model(z,innovation,params,steady_state,i+1);
|
||||||
|
end
|
||||||
|
res(:,i,1) = res(:,i,1)+weights(k)*d1;
|
||||||
|
end
|
||||||
|
if nargin > 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 nargin > 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 nargin > 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(:,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;
|
||||||
|
if nargin > 1
|
||||||
|
i_cols_Af = i_cols_Af + ny;
|
||||||
|
end
|
||||||
|
else
|
||||||
|
% existing worlds other than 1
|
||||||
|
if nargin > 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 nargin > 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(:,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;
|
||||||
|
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 nargin > 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 nargin > 1
|
||||||
|
[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_j));
|
||||||
|
else
|
||||||
|
[ir,ic,v] = find(jacobian(:,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;
|
||||||
|
i_rows_y = i_rows_y + ny;
|
||||||
|
offset_c = offset_c + world_nbr*ny;
|
||||||
|
offset_r = offset_r + world_nbr*ny;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
if nargin > 1
|
||||||
|
iA = [nzA{:}]';
|
||||||
|
A = sparse(iA(:,1),iA(:,2),iA(:,3),dimension,dimension);
|
||||||
|
end
|
||||||
|
res = res(pm.i_upd_r);
|
|
@ -1,4 +1,4 @@
|
||||||
function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,EpOptions,pfm,order,varargin)
|
function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,Options,pfm,order,varargin)
|
||||||
|
|
||||||
% Copyright (C) 2012-2013 Dynare Team
|
% Copyright (C) 2012-2013 Dynare Team
|
||||||
%
|
%
|
||||||
|
@ -17,281 +17,135 @@ function [flag,endo_simul,err] = solve_stochastic_perfect_foresight_model_1(endo
|
||||||
% You should have received a copy of the GNU General Public License
|
% You should have received a copy of the GNU General Public License
|
||||||
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
|
||||||
if nargin < 6
|
global options_
|
||||||
homotopy_parameter = 1;
|
|
||||||
else
|
if nargin < 6
|
||||||
homotopy_parameter = varargin{1};
|
homotopy_parameter = 1;
|
||||||
|
else
|
||||||
|
homotopy_parameter = varargin{1};
|
||||||
|
end
|
||||||
|
flag = 0;
|
||||||
|
err = 0;
|
||||||
|
stop = 0;
|
||||||
|
|
||||||
|
EpOptions = Options.ep;
|
||||||
|
|
||||||
|
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;
|
||||||
|
nyp = pfm.nyp;
|
||||||
|
nyf = pfm.nyf;
|
||||||
|
i_cols_1 = pfm.i_cols_1;
|
||||||
|
i_cols_A1 = pfm.i_cols_A1;
|
||||||
|
i_cols_j = pfm.i_cols_j;
|
||||||
|
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
|
||||||
|
hybrid_order = pfm.hybrid_order;
|
||||||
|
dr = pfm.dr;
|
||||||
|
nodes = pfm.nodes;
|
||||||
|
weights = pfm.weights;
|
||||||
|
nnodes = pfm.nnodes;
|
||||||
|
|
||||||
|
maxit = pfm.maxit_;
|
||||||
|
tolerance = pfm.tolerance;
|
||||||
|
verbose = pfm.verbose;
|
||||||
|
|
||||||
|
number_of_shocks = size(exo_simul,2);
|
||||||
|
|
||||||
|
% make sure that there is a node equal to zero
|
||||||
|
% and permute nodes and weights to have zero first
|
||||||
|
k = find(sum(abs(nodes),2) < 1e-12);
|
||||||
|
if ~isempty(k)
|
||||||
|
nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
|
||||||
|
weights = [weights(k); weights(1:k-1); weights(k+1:end)];
|
||||||
|
else
|
||||||
|
error('there is no nodes equal to zero')
|
||||||
|
end
|
||||||
|
|
||||||
|
if hybrid_order > 0
|
||||||
|
if hybrid_order == 2
|
||||||
|
h_correction = 0.5*dr.ghs2(dr.inv_order_var);
|
||||||
end
|
end
|
||||||
flag = 0;
|
else
|
||||||
err = 0;
|
h_correction = 0;
|
||||||
stop = 0;
|
end
|
||||||
|
|
||||||
params = pfm.params;
|
if verbose
|
||||||
steady_state = pfm.steady_state;
|
disp ([' -----------------------------------------------------']);
|
||||||
ny = pfm.ny;
|
disp (['MODEL SIMULATION :']);
|
||||||
periods = pfm.periods;
|
fprintf('\n');
|
||||||
dynamic_model = pfm.dynamic_model;
|
end
|
||||||
lead_lag_incidence = pfm.lead_lag_incidence;
|
|
||||||
nyp = pfm.nyp;
|
|
||||||
nyf = pfm.nyf;
|
|
||||||
i_cols_1 = pfm.i_cols_1;
|
|
||||||
i_cols_A1 = pfm.i_cols_A1;
|
|
||||||
i_cols_j = pfm.i_cols_j;
|
|
||||||
i_cols_T = nonzeros(lead_lag_incidence(1:2,:)');
|
|
||||||
hybrid_order = pfm.hybrid_order;
|
|
||||||
dr = pfm.dr;
|
|
||||||
[nodes,weights,nnodes] = setup_integration_nodes(EpOptions,pfm);
|
|
||||||
|
|
||||||
maxit = pfm.maxit_;
|
|
||||||
tolerance = pfm.tolerance;
|
|
||||||
verbose = pfm.verbose;
|
|
||||||
|
|
||||||
number_of_shocks = size(exo_simul,2);
|
% Each column of Y represents a different world
|
||||||
|
% The upper right cells are unused
|
||||||
|
% The first row block is ny x 1
|
||||||
|
% The second row block is ny x nnodes
|
||||||
|
% The third row block is ny x nnodes^2
|
||||||
|
% and so on until size ny x nnodes^order
|
||||||
|
world_nbr = 1+(nnodes-1)*order;
|
||||||
|
Y = endo_simul(:,2:end-1);
|
||||||
|
Y = repmat(Y,1,world_nbr);
|
||||||
|
pfm.y0 = endo_simul(:,1);
|
||||||
|
|
||||||
% make sure that there is a node equal to zero
|
% The columns of A map the elements of Y such that
|
||||||
% and permute nodes and weights to have zero first
|
% each block of Y with ny rows are unfolded column wise
|
||||||
k = find(sum(abs(nodes),2) < 1e-12);
|
% number of blocks
|
||||||
if ~isempty(k)
|
block_nbr = (order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr);
|
||||||
nodes = [nodes(k,:); nodes(1:k-1,:); nodes(k+1:end,:)];
|
% dimension of the problem
|
||||||
weights = [weights(k); weights(1:k-1); weights(k+1:end)];
|
dimension = ny*block_nbr;
|
||||||
else
|
pfm.block_nbr = block_nbr;
|
||||||
error('there is no nodes equal to zero')
|
pfm.dimension = dimension;
|
||||||
end
|
if order == 0
|
||||||
|
i_upd_r = (1:ny*periods);
|
||||||
if hybrid_order > 0
|
i_upd_y = i_upd_r + ny;
|
||||||
if hybrid_order == 2
|
else
|
||||||
h_correction = 0.5*dr.ghs2(dr.inv_order_var);
|
i_upd_r = zeros(dimension,1);
|
||||||
end
|
i_upd_y = i_upd_r;
|
||||||
|
i_upd_r(1:ny) = (1:ny);
|
||||||
|
i_upd_y(1:ny) = ny+(1:ny);
|
||||||
|
i1 = ny+1;
|
||||||
|
i2 = 2*ny;
|
||||||
|
n1 = ny+1;
|
||||||
|
n2 = 2*ny;
|
||||||
|
for i=2:periods
|
||||||
|
k = n1:n2;
|
||||||
|
for j=1:(1+(nnodes-1)*min(i-1,order))
|
||||||
|
i_upd_r(i1:i2) = k+(j-1)*ny*periods;
|
||||||
|
i_upd_y(i1:i2) = k+ny+(j-1)*ny*(periods+2);
|
||||||
|
i1 = i2+1;
|
||||||
|
i2 = i2+ny;
|
||||||
|
end
|
||||||
|
n1 = n2+1;
|
||||||
|
n2 = n2+ny;
|
||||||
end
|
end
|
||||||
|
end
|
||||||
|
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
|
||||||
|
find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
|
||||||
|
h1 = clock;
|
||||||
|
pfm.order = order;
|
||||||
|
pfm.world_nbr = world_nbr;
|
||||||
|
pfm.nodes = nodes;
|
||||||
|
pfm.nnodes = nnodes;
|
||||||
|
pfm.weights = weights;
|
||||||
|
pfm.h_correction = h_correction;
|
||||||
|
pfm.i_rows = 1:ny;
|
||||||
|
i_cols = find(lead_lag_incidence');
|
||||||
|
pfm.i_cols = i_cols;
|
||||||
|
pfm.nyp = nyp;
|
||||||
|
pfm.nyf = nyf;
|
||||||
|
pfm.hybrid_order = hybrid_order;
|
||||||
|
pfm.i_cols_1 = i_cols_1;
|
||||||
|
pfm.i_cols_h = i_cols_j;
|
||||||
|
pfm.icA = icA;
|
||||||
|
pfm.i_cols_T = i_cols_T;
|
||||||
|
pfm.i_upd_r = i_upd_r;
|
||||||
|
pfm.i_upd_y = i_upd_y;
|
||||||
|
|
||||||
if verbose
|
|
||||||
disp ([' -----------------------------------------------------']);
|
|
||||||
disp (['MODEL SIMULATION :']);
|
|
||||||
fprintf('\n');
|
|
||||||
end
|
|
||||||
|
|
||||||
z = endo_simul(find(lead_lag_incidence'));
|
options_.solve_algo = 9;
|
||||||
[d1,jacobian] = dynamic_model(z,exo_simul,params,steady_state,2);
|
options_.steady.maxit = 100;
|
||||||
|
y = repmat(steady_state,block_nbr,1);
|
||||||
% Each column of Y represents a different world
|
y = dynare_solve(@ep_problem_2,y,1,exo_simul,pfm);
|
||||||
% The upper right cells are unused
|
endo_simul(:,2) = y(1:ny);
|
||||||
% The first row block is ny x 1
|
|
||||||
% The second row block is ny x nnodes
|
|
||||||
% The third row block is ny x nnodes^2
|
|
||||||
% and so on until size ny x nnodes^order
|
|
||||||
world_nbr = 1+(nnodes-1)*order;
|
|
||||||
Y = repmat(endo_simul(:),1,world_nbr);
|
|
||||||
|
|
||||||
% The columns of A map the elements of Y such that
|
|
||||||
% each block of Y with ny rows are unfolded column wise
|
|
||||||
dimension = ny*(order+(nnodes-1)*(order-1)*order/2+(periods-order)*world_nbr);
|
|
||||||
if order == 0
|
|
||||||
i_upd_r = (1:ny*periods);
|
|
||||||
i_upd_y = i_upd_r + ny;
|
|
||||||
else
|
|
||||||
i_upd_r = zeros(dimension,1);
|
|
||||||
i_upd_y = i_upd_r;
|
|
||||||
i_upd_r(1:ny) = (1:ny);
|
|
||||||
i_upd_y(1:ny) = ny+(1:ny);
|
|
||||||
i1 = ny+1;
|
|
||||||
i2 = 2*ny;
|
|
||||||
n1 = ny+1;
|
|
||||||
n2 = 2*ny;
|
|
||||||
for i=2:periods
|
|
||||||
k = n1:n2;
|
|
||||||
for j=1:(1+(nnodes-1)*min(i-1,order))
|
|
||||||
i_upd_r(i1:i2) = k+(j-1)*ny*periods;
|
|
||||||
i_upd_y(i1:i2) = k+ny+(j-1)*ny*(periods+2);
|
|
||||||
i1 = i2+1;
|
|
||||||
i2 = i2+ny;
|
|
||||||
end
|
|
||||||
n1 = n2+1;
|
|
||||||
n2 = n2+ny;
|
|
||||||
end
|
|
||||||
end
|
|
||||||
icA = [find(lead_lag_incidence(1,:)) find(lead_lag_incidence(2,:))+world_nbr*ny ...
|
|
||||||
find(lead_lag_incidence(3,:))+2*world_nbr*ny]';
|
|
||||||
h1 = clock;
|
|
||||||
for iter = 1:maxit
|
|
||||||
h2 = clock;
|
|
||||||
A1 = sparse([],[],[],ny*(order+(nnodes-1)*(order-1)*order/2),dimension,(order+1)*world_nbr*nnz(jacobian));
|
|
||||||
res = zeros(ny,periods,world_nbr);
|
|
||||||
i_rows = 1:ny;
|
|
||||||
i_cols = find(lead_lag_incidence');
|
|
||||||
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_A = i_cols;
|
|
||||||
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;
|
|
||||||
for i = 1:order+1
|
|
||||||
i_w_p = 1;
|
|
||||||
for j = 1:(1+(nnodes-1)*(i-1))
|
|
||||||
innovation = exo_simul;
|
|
||||||
if i <= order && j == 1
|
|
||||||
% first world, integrating future shocks
|
|
||||||
for k=1:nnodes
|
|
||||||
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
|
|
||||||
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 == 2 && (k > 1 || i == order)
|
|
||||||
y = [Y(i_cols_p,1);
|
|
||||||
Y(i_cols_s,1);
|
|
||||||
Y(i_cols_f,k1)+h_correction(i_hc)];
|
|
||||||
else
|
|
||||||
y = [Y(i_cols_p,1);
|
|
||||||
Y(i_cols_s,1);
|
|
||||||
Y(i_cols_f,k1)];
|
|
||||||
end
|
|
||||||
[d1,jacobian] = dynamic_model(y,homotopy_parameter*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,1) = res(:,i,1)+weights(k)*d1;
|
|
||||||
end
|
|
||||||
elseif j > 1 + (nnodes-1)*(i-2)
|
|
||||||
% new world, using previous state of world 1 and hit
|
|
||||||
% by shocks from nodes
|
|
||||||
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);
|
|
||||||
k = j - (nnodes-1)*(i-2);
|
|
||||||
innovation(i+1,:) = nodes(k,:);
|
|
||||||
y = [Y(i_cols_p,1);
|
|
||||||
Y(i_cols_s,j);
|
|
||||||
Y(i_cols_f,j)];
|
|
||||||
[d1,jacobian] = dynamic_model(y,homotopy_parameter*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);
|
|
||||||
res(:,i,j) = d1;
|
|
||||||
i_cols_Af = i_cols_Af + ny;
|
|
||||||
else
|
|
||||||
% existing worlds other than 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);
|
|
||||||
y = [Y(i_cols_p,j);
|
|
||||||
Y(i_cols_s,j);
|
|
||||||
Y(i_cols_f,j)];
|
|
||||||
[d1,jacobian] = dynamic_model(y,homotopy_parameter*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);
|
|
||||||
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
|
|
||||||
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,world_nbr);
|
|
||||||
parfor j=1:world_nbr
|
|
||||||
i_rows_y = find(lead_lag_incidence')+(order+1)*ny;
|
|
||||||
offset_c = ny*(order+(nnodes-1)*(order-1)*order/2+j-1);
|
|
||||||
offset_r = (j-1)*ny;
|
|
||||||
for i=order+2:periods
|
|
||||||
[d1,jacobian] = dynamic_model(Y(i_rows_y,j), ...
|
|
||||||
exo_simul,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+icA(ic), v]';
|
|
||||||
res(:,i,j) = d1;
|
|
||||||
i_rows_y = i_rows_y + ny;
|
|
||||||
offset_c = offset_c + world_nbr*ny;
|
|
||||||
offset_r = offset_r + world_nbr*ny;
|
|
||||||
end
|
|
||||||
end
|
|
||||||
err = max(abs(res(i_upd_r)));
|
|
||||||
if verbose
|
|
||||||
[err1, k1] = max(abs(res));
|
|
||||||
[err2, k2] = max(abs(err1));
|
|
||||||
[err3, k3] = max(abs(err2));
|
|
||||||
disp([iter err k1(:,k2(:,:,k3),k3) k2(:,:,k3) k3])
|
|
||||||
end
|
|
||||||
if err < tolerance
|
|
||||||
stop = 1;
|
|
||||||
flag = 0;% Convergency obtained.
|
|
||||||
endo_simul = reshape(Y(:,1),ny,periods+2);%Y(ny+(1:ny),1);
|
|
||||||
if verbose
|
|
||||||
save ep_test_s1 exo_simul endo_simul Y
|
|
||||||
fprintf('\n') ;
|
|
||||||
disp([' Total time of simulation :' num2str(etime(clock,h1))]) ;
|
|
||||||
fprintf('\n') ;
|
|
||||||
disp([' Convergency obtained.']) ;
|
|
||||||
fprintf('\n') ;
|
|
||||||
end
|
|
||||||
break
|
|
||||||
end
|
|
||||||
A2 = [nzA{:}]';
|
|
||||||
if any(isnan(A2(:,3))) || any(any(any(isnan(res))))
|
|
||||||
if verbose
|
|
||||||
disp(['solve_stochastic_foresight_model_1 encountered ' ...
|
|
||||||
'NaN'])
|
|
||||||
save ep_test_s2 exo_simul endo_simul
|
|
||||||
pause
|
|
||||||
end
|
|
||||||
flag = 1;
|
|
||||||
return
|
|
||||||
end
|
|
||||||
A = [A1; sparse(A2(:,1),A2(:,2),A2(:,3),ny*(periods-order-1)*world_nbr,dimension)];
|
|
||||||
if verbose
|
|
||||||
disp(sprintf('condest %g',condest(A)))
|
|
||||||
end
|
|
||||||
dy = -A\res(i_upd_r);
|
|
||||||
Y(i_upd_y) = Y(i_upd_y) + dy;
|
|
||||||
end
|
|
||||||
|
|
||||||
if ~stop
|
|
||||||
if verbose
|
|
||||||
fprintf('\n') ;
|
|
||||||
disp([' Total time of simulation :' num2str(etime(clock,h1))]) ;
|
|
||||||
fprintf('\n') ;
|
|
||||||
disp(['WARNING : maximum number of iterations is reached (modify options_.simul.maxit).']) ;
|
|
||||||
fprintf('\n') ;
|
|
||||||
disp(sprintf('err: %f',err));
|
|
||||||
save ep_test_s2 exo_simul endo_simul
|
|
||||||
pause
|
|
||||||
end
|
|
||||||
flag = 1;% more iterations are needed.
|
|
||||||
endo_simul = 1;
|
|
||||||
end
|
|
||||||
if verbose
|
|
||||||
disp (['-----------------------------------------------------']) ;
|
|
||||||
end
|
|
|
@ -0,0 +1,9 @@
|
||||||
|
function [X,w]=stroud_judd_7.5.8(d)
|
||||||
|
|
||||||
|
E = eye(d);
|
||||||
|
X = cell(2*d,1);
|
||||||
|
m = 1;
|
||||||
|
for i=1:d
|
||||||
|
X{m} = E(:,i);
|
||||||
|
m = m+1;
|
||||||
|
X{m} = -E(:,i);
|
Loading…
Reference in New Issue