removed globals from solve1.m
parent
267b8fd995
commit
70e162c736
|
@ -83,6 +83,7 @@ function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
|
||||||
info = 0;
|
info = 0;
|
||||||
lead_lag_incidence = M.lead_lag_incidence;
|
lead_lag_incidence = M.lead_lag_incidence;
|
||||||
order_var = dr.order_var;
|
order_var = dr.order_var;
|
||||||
|
endo_nbr = M.endo_nbr;
|
||||||
exo_nbr = M.exo_nbr;
|
exo_nbr = M.exo_nbr;
|
||||||
|
|
||||||
M.var_order_endo_names = M.endo_names(dr.order_var,:);
|
M.var_order_endo_names = M.endo_names(dr.order_var,:);
|
||||||
|
@ -104,420 +105,283 @@ function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
|
||||||
end
|
end
|
||||||
|
|
||||||
if isfield(options,'portfolio') && options.portfolio == 1
|
if isfield(options,'portfolio') && options.portfolio == 1
|
||||||
eq_tags = M.equations_tags;
|
pm = portfolio_model_structure(M,options);
|
||||||
n_tags = size(eq_tags,1);
|
|
||||||
l_var = zeros(n_tags,1);
|
|
||||||
for i=1:n_tags
|
|
||||||
l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
|
|
||||||
length(cell2mat(eq_tags(i,3)))));
|
|
||||||
end
|
|
||||||
dr.ys = ys0;
|
|
||||||
x0 = ys0(l_var);
|
|
||||||
% dr = first_step_ds(x0,M,dr,options,oo);
|
|
||||||
n = size(ys0);
|
|
||||||
%x0 = ys0;
|
|
||||||
[x, info] = solve1(@risky_residuals_ds,x0,1:n_tags,1:n_tags,0,1,M,dr,options,oo);
|
|
||||||
%[x, info] = solve1(@risky_residuals,x0,1:n,1:n,0,1,M,dr,options,oo);
|
|
||||||
% ys0(l_var) = x;
|
|
||||||
ys0(l_var) = x;
|
|
||||||
dr.ys = ys0;
|
|
||||||
oo.dr = dr;
|
|
||||||
oo.steady_state = ys0;
|
|
||||||
disp_steady_state(M,oo);
|
|
||||||
end
|
|
||||||
|
|
||||||
[ys, info] = csolve(func,ys0,[],1e-10,100,M,dr,options,oo);
|
x0 = ys0(pm.v_p);
|
||||||
|
n = length(x0);
|
||||||
if options.k_order_solver
|
[x, info] = solve1(@risky_residuals_ds,x0,1:n,1:n,0,1, options.gstep, ...
|
||||||
[resid,dr] = risky_residuals_k_order(ys,M,dr,options,oo);
|
options.solve_tolf,options.solve_tolx, ...
|
||||||
else
|
options.solve_maxit,options.debug,pm,M,dr, ...
|
||||||
[resid,dr] = risky_residuals(ys,M,dr,options,oo);
|
options,oo);
|
||||||
end
|
|
||||||
|
|
||||||
dr.ys = ys;
|
|
||||||
for i=1:M.endo_nbr
|
|
||||||
disp(sprintf('%16s %12.6f %12.6f',M.endo_names(i,:),ys0(i), ys(i)))
|
|
||||||
end
|
|
||||||
|
|
||||||
dr.ghs2 = zeros(size(dr.ghs2));
|
|
||||||
|
|
||||||
k_var = setdiff(1:M.endo_nbr,l_var);
|
|
||||||
dr.ghx(k_var,:) = dr.ghx;
|
|
||||||
dr.ghu(k_var,:) = dr.ghu;
|
|
||||||
dr.ghxx(k_var,:) = dr.ghxx;
|
|
||||||
dr.ghxu(k_var,:) = dr.ghxu;
|
|
||||||
dr.ghuu(k_var,:) = dr.ghuu;
|
|
||||||
dr.ghs2(k_var,:) = dr.ghs2;
|
|
||||||
end
|
|
||||||
|
|
||||||
function [resid,dr] = risky_residuals(ys,M,dr,options,oo)
|
|
||||||
persistent old_ys old_resid
|
|
||||||
|
|
||||||
lead_lag_incidence = M.lead_lag_incidence;
|
|
||||||
iyv = lead_lag_incidence';
|
|
||||||
iyv = iyv(:);
|
|
||||||
iyr0 = find(iyv) ;
|
|
||||||
|
|
||||||
if M.exo_nbr == 0
|
|
||||||
oo.exo_steady_state = [] ;
|
|
||||||
end
|
|
||||||
|
|
||||||
z = repmat(ys,1,3);
|
|
||||||
z = z(iyr0) ;
|
|
||||||
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
|
||||||
[oo.exo_simul ...
|
|
||||||
oo.exo_det_simul], M.params, dr.ys, 2);
|
|
||||||
if ~isreal(d1) || ~isreal(d2)
|
|
||||||
pause
|
|
||||||
end
|
|
||||||
|
|
||||||
if options.use_dll
|
|
||||||
% In USE_DLL mode, the hessian is in the 3-column sparse representation
|
|
||||||
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
|
||||||
size(d1, 1), size(d1, 2)*size(d1, 2));
|
|
||||||
end
|
|
||||||
|
|
||||||
if isfield(options,'portfolio') && options.portfolio == 1
|
|
||||||
eq_tags = M.equations_tags;
|
|
||||||
n_tags = size(eq_tags,1);
|
|
||||||
portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
|
|
||||||
'portfolio'),1));
|
|
||||||
eq = setdiff(1:M.endo_nbr,portfolios_eq);
|
|
||||||
l_var = zeros(n_tags,1);
|
|
||||||
for i=1:n_tags
|
|
||||||
l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
|
|
||||||
length(cell2mat(eq_tags(i,3)))));
|
|
||||||
end
|
|
||||||
k_var = setdiff(1:M.endo_nbr,l_var);
|
|
||||||
lli1 = lead_lag_incidence(:,k_var);
|
|
||||||
lead_incidence = lli1(3,:)';
|
|
||||||
k = find(lli1');
|
|
||||||
lli2 = lli1';
|
|
||||||
lli2(k) = 1:nnz(lli1);
|
|
||||||
lead_lag_incidence = lli2';
|
|
||||||
x = ys(l_var);
|
|
||||||
dr = first_step_ds(x,M,dr,options,oo);
|
|
||||||
|
|
||||||
|
|
||||||
M.lead_lag_incidence = lead_lag_incidence;
|
|
||||||
lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
|
|
||||||
d1a = d1(eq,lli1a);
|
|
||||||
ih = 1:size(d2,2);
|
|
||||||
ih = reshape(ih,size(d1,2),size(d1,2));
|
|
||||||
ih1 = ih(lli1a,lli1a);
|
|
||||||
d2a = d2(eq,ih1);
|
|
||||||
|
|
||||||
M.endo_nbr = M.endo_nbr-n_tags;
|
|
||||||
dr = set_state_space(dr,M,options);
|
|
||||||
|
|
||||||
[junk,dr.i_fwrd_g] = find(lead_lag_incidence(3,dr.order_var));
|
|
||||||
i_fwrd_f = nonzeros(lead_incidence(dr.order_var));
|
|
||||||
i_fwrd2_f = ih(i_fwrd_f,i_fwrd_f);
|
|
||||||
dr.i_fwrd_f = i_fwrd_f;
|
|
||||||
dr.i_fwrd2_f = i_fwrd2_f;
|
|
||||||
nd = nnz(lead_lag_incidence) + M.exo_nbr;
|
|
||||||
dr.nd = nd;
|
|
||||||
kk = reshape(1:nd^2,nd,nd);
|
|
||||||
kkk = reshape(1:nd^3,nd^2,nd);
|
|
||||||
dr.i_fwrd2a_f = kk(i_fwrd_f,:);
|
|
||||||
% dr.i_fwrd3_f = kkk(i_fwrd2_f,:);
|
|
||||||
dr.i_uu = kk(end-M.exo_nbr+1:end,end-M.exo_nbr+1:end);
|
|
||||||
else
|
|
||||||
d1a = d1;
|
|
||||||
d2a = d2;
|
|
||||||
end
|
|
||||||
|
|
||||||
% $$$ [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
|
|
||||||
% $$$ b = zeros(M.endo_nbr,M.endo_nbr);
|
|
||||||
% $$$ b(:,cols_b) = d1a(:,cols_j);
|
|
||||||
% $$$
|
|
||||||
% $$$ [dr,info] = dyn_first_order_solver(d1a,b,M,dr,options,0);
|
|
||||||
% $$$ if info
|
|
||||||
% $$$ [m1,m2]=max(abs(ys-old_ys));
|
|
||||||
% $$$ disp([m1 m2])
|
|
||||||
% $$$ % print_info(info,options.noprint);
|
|
||||||
% $$$ resid = old_resid+info(2)/40;
|
|
||||||
% $$$ return
|
|
||||||
% $$$ end
|
|
||||||
% $$$
|
|
||||||
% $$$ dr = dyn_second_order_solver(d1a,d2a,dr,M);
|
|
||||||
|
|
||||||
gu1 = dr.ghu(dr.i_fwrd_g,:);
|
|
||||||
|
|
||||||
resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)+ ...
|
|
||||||
d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
|
|
||||||
disp(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)*vec(M.Sigma_e));
|
|
||||||
old_ys = ys;
|
|
||||||
disp(max(abs(resid)))
|
|
||||||
old_resid = resid;
|
|
||||||
end
|
|
||||||
|
|
||||||
function [resid,dr] = risky_residuals_ds(x,M,dr,options,oo)
|
|
||||||
persistent old_ys old_resid old_resid1 old_d1 old_d2
|
|
||||||
|
|
||||||
dr = first_step_ds(x,M,dr,options,oo);
|
|
||||||
|
|
||||||
lead_lag_incidence = M.lead_lag_incidence;
|
|
||||||
iyv = lead_lag_incidence';
|
|
||||||
iyv = iyv(:);
|
|
||||||
iyr0 = find(iyv) ;
|
|
||||||
|
|
||||||
if M.exo_nbr == 0
|
|
||||||
oo.exo_steady_state = [] ;
|
|
||||||
end
|
|
||||||
|
|
||||||
eq_tags = M.equations_tags;
|
|
||||||
n_tags = size(eq_tags,1);
|
|
||||||
portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
|
|
||||||
'portfolio'),1));
|
|
||||||
eq = setdiff(1:M.endo_nbr,portfolios_eq);
|
|
||||||
l_var = zeros(n_tags,1);
|
|
||||||
for i=1:n_tags
|
|
||||||
l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
|
|
||||||
length(cell2mat(eq_tags(i,3)))));
|
|
||||||
end
|
|
||||||
k_var = setdiff(1:M.endo_nbr,l_var);
|
|
||||||
lli1 = lead_lag_incidence(:,k_var);
|
|
||||||
k = find(lli1');
|
|
||||||
lli2 = lli1';
|
|
||||||
lli2(k) = 1:nnz(lli1);
|
|
||||||
lead_lag_incidence = lli2';
|
|
||||||
|
|
||||||
ys = dr.ys;
|
|
||||||
ys(l_var) = x;
|
|
||||||
|
|
||||||
z = repmat(ys,1,3);
|
|
||||||
z = z(iyr0) ;
|
|
||||||
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
|
||||||
[oo.exo_simul ...
|
|
||||||
oo.exo_det_simul], M.params, dr.ys, 2);
|
|
||||||
% $$$ if isempty(old_resid)
|
|
||||||
% $$$ old_resid1 = resid1;
|
|
||||||
% $$$ old_d1 = d1;
|
|
||||||
% $$$ old_d2 = d2;
|
|
||||||
% $$$ old_ys = ys;
|
|
||||||
% $$$ else
|
|
||||||
% $$$ if ~isequal(resid1,old_resid)
|
|
||||||
% $$$ disp('ys')
|
|
||||||
% $$$ disp((ys-old_ys)');
|
|
||||||
% $$$ disp('resids1')
|
|
||||||
% $$$ disp((resid1-old_resid1)')
|
|
||||||
% $$$ old_resid1 = resid1;
|
|
||||||
% $$$ pause
|
|
||||||
% $$$ end
|
|
||||||
% $$$ if ~isequal(d1,old_d1)
|
|
||||||
% $$$ disp('d1')
|
|
||||||
% $$$ disp(d1-old_d1);
|
|
||||||
% $$$ old_d1 = d1;
|
|
||||||
% $$$ pause
|
|
||||||
% $$$ end
|
|
||||||
% $$$ if ~isequal(d2,old_d2)
|
|
||||||
% $$$ disp('d2')
|
|
||||||
% $$$ disp(d2-old_d2);
|
|
||||||
% $$$ old_d2 = d2;
|
|
||||||
% $$$ pause
|
|
||||||
% $$$ end
|
|
||||||
% $$$ end
|
|
||||||
if ~isreal(d1) || ~isreal(d2)
|
|
||||||
pause
|
|
||||||
end
|
|
||||||
|
|
||||||
if options.use_dll
|
|
||||||
% In USE_DLL mode, the hessian is in the 3-column sparse representation
|
|
||||||
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
|
||||||
size(d1, 1), size(d1, 2)*size(d1, 2));
|
|
||||||
end
|
|
||||||
|
|
||||||
% $$$ if isfield(options,'portfolio') && options.portfolio == 1
|
|
||||||
% $$$ lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
|
|
||||||
% $$$ d1a = d1(eq,lli1a);
|
|
||||||
% $$$ ih = 1:size(d2,2);
|
|
||||||
% $$$ ih = reshape(ih,size(d1,2),size(d1,2));
|
|
||||||
% $$$ ih1 = ih(lli1a,lli1a);
|
|
||||||
% $$$ d2a = d2(eq,ih1);
|
|
||||||
% $$$
|
|
||||||
% $$$ M.endo_nbr = M.endo_nbr-n_tags;
|
|
||||||
% $$$ dr = set_state_space(dr,M);
|
|
||||||
% $$$
|
|
||||||
% $$$ dr.i_fwrd_g = find(lead_lag_incidence(3,dr.order_var)');
|
|
||||||
% $$$ else
|
|
||||||
% $$$ d1a = d1;
|
|
||||||
% $$$ d2a = d2;
|
|
||||||
% $$$ end
|
|
||||||
% $$$
|
|
||||||
% $$$ [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
|
|
||||||
% $$$ b = zeros(M.endo_nbr,M.endo_nbr);
|
|
||||||
% $$$ b(:,cols_b) = d1a(:,cols_j);
|
|
||||||
% $$$
|
|
||||||
% $$$ [dr,info] = dyn_first_order_solver(d1a,b,M,dr,options,0);
|
|
||||||
% $$$ if info
|
|
||||||
% $$$ [m1,m2]=max(abs(ys-old_ys));
|
|
||||||
% $$$ disp([m1 m2])
|
|
||||||
% $$$ % print_info(info,options.noprint);
|
|
||||||
% $$$ resid = old_resid+info(2)/40;
|
|
||||||
% $$$ return
|
|
||||||
% $$$ end
|
|
||||||
% $$$
|
|
||||||
% $$$ dr = dyn_second_order_solver(d1a,d2a,dr,M);
|
|
||||||
|
|
||||||
gu1 = dr.ghu(dr.i_fwrd_g,:);
|
|
||||||
|
|
||||||
% resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*dr.ghuu(dr.i_fwrd_g,:)+ ...
|
|
||||||
% d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
|
|
||||||
resid = resid1+0.5*(d2(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
|
|
||||||
|
|
||||||
% $$$ if isempty(old_resid)
|
|
||||||
% $$$ old_resid = resid;
|
|
||||||
% $$$ else
|
|
||||||
% $$$ disp('resid')
|
|
||||||
% $$$ dr = (resid-old_resid)';
|
|
||||||
% $$$ % disp(dr)
|
|
||||||
% $$$ % disp(dr(portfolios_eq))
|
|
||||||
% $$$ old_resid = resid;
|
|
||||||
% $$$ end
|
|
||||||
resid = resid(portfolios_eq)
|
|
||||||
end
|
|
||||||
|
|
||||||
function [dr] = first_step_ds(x,M,dr,options,oo)
|
|
||||||
|
|
||||||
lead_lag_incidence = M.lead_lag_incidence;
|
|
||||||
iyv = lead_lag_incidence';
|
|
||||||
iyv = iyv(:);
|
|
||||||
iyr0 = find(iyv) ;
|
|
||||||
|
|
||||||
if M.exo_nbr == 0
|
|
||||||
oo.exo_steady_state = [] ;
|
|
||||||
end
|
|
||||||
|
|
||||||
eq_tags = M.equations_tags;
|
|
||||||
n_tags = size(eq_tags,1);
|
|
||||||
portfolios_eq = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
|
|
||||||
'portfolio'),1));
|
|
||||||
eq = setdiff(1:M.endo_nbr,portfolios_eq);
|
|
||||||
l_var = zeros(n_tags,1);
|
|
||||||
for i=1:n_tags
|
|
||||||
l_var(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
|
|
||||||
length(cell2mat(eq_tags(i,3)))));
|
|
||||||
end
|
|
||||||
k_var = setdiff(1:M.endo_nbr,l_var);
|
|
||||||
lli1 = lead_lag_incidence(:,k_var);
|
|
||||||
k = find(lli1');
|
|
||||||
lli2 = lli1';
|
|
||||||
lli2(k) = 1:nnz(lli1);
|
|
||||||
lead_lag_incidence = lli2';
|
|
||||||
M.lead_lag_incidence = lead_lag_incidence;
|
|
||||||
|
|
||||||
ys = dr.ys;
|
|
||||||
ys(l_var) = x;
|
|
||||||
|
|
||||||
z = repmat(ys,1,3);
|
|
||||||
z = z(iyr0) ;
|
|
||||||
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
|
||||||
[oo.exo_simul ...
|
|
||||||
oo.exo_det_simul], M.params, dr.ys, 2);
|
|
||||||
if ~isreal(d1) || ~isreal(d2)
|
|
||||||
pause
|
|
||||||
end
|
|
||||||
|
|
||||||
if options.use_dll
|
|
||||||
% In USE_DLL mode, the hessian is in the 3-column sparse representation
|
|
||||||
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
|
||||||
size(d1, 1), size(d1, 2)*size(d1, 2));
|
|
||||||
end
|
|
||||||
|
|
||||||
if isfield(options,'portfolio') && options.portfolio == 1
|
|
||||||
lli1a = [nonzeros(lli1'); size(d1,2)+(-M.exo_nbr+1:0)'];
|
|
||||||
d1a = d1(eq,lli1a);
|
|
||||||
ih = 1:size(d2,2);
|
|
||||||
ih = reshape(ih,size(d1,2),size(d1,2));
|
|
||||||
ih1 = ih(lli1a,lli1a);
|
|
||||||
d2a = d2(eq,ih1);
|
|
||||||
|
|
||||||
M.endo_nbr = M.endo_nbr-n_tags;
|
|
||||||
dr = set_state_space(dr,M,options);
|
|
||||||
|
|
||||||
dr.i_fwrd_g = find(lead_lag_incidence(3,dr.order_var)');
|
|
||||||
else
|
|
||||||
d1a = d1;
|
|
||||||
d2a = d2;
|
|
||||||
end
|
|
||||||
|
|
||||||
[junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
|
|
||||||
b = zeros(M.endo_nbr,M.endo_nbr);
|
|
||||||
b(:,cols_b) = d1a(:,cols_j);
|
|
||||||
|
|
||||||
[dr,info] = dyn_first_order_solver(d1a,M,dr,options,0);
|
|
||||||
if info
|
if info
|
||||||
[m1,m2]=max(abs(ys-old_ys));
|
error('DS approach can''t be computed')
|
||||||
disp([m1 m2])
|
end
|
||||||
% print_info(info,options.noprint);
|
%[x, info] = csolve(@risky_residuals_ds,x0,[],1e-10,100,M,dr,options,oo);
|
||||||
resid = old_resid+info(2)/40;
|
% ys0(l_var) = x;
|
||||||
|
[resids,dr1] = risky_residuals_ds(x,pm,M,dr,options,oo);
|
||||||
|
ys1 = dr1.ys;
|
||||||
|
end
|
||||||
|
|
||||||
|
[ys, info] = solve1(func,ys0,1:endo_nbr,1:endo_nbr,0,1, options.gstep, ...
|
||||||
|
options.solve_tolf,options.solve_tolx, ...
|
||||||
|
options.solve_maxit,options.debug,pm,M,dr,options,oo);
|
||||||
|
% [ys, info] = csolve(func,ys0,[],1e-10,100,M,dr,options,oo);
|
||||||
|
if info
|
||||||
|
error('RSS approach can''t be computed')
|
||||||
|
end
|
||||||
|
dr.ys = ys;
|
||||||
|
|
||||||
|
[resid,dr] = func(ys,pm,M,dr,options,oo);
|
||||||
|
|
||||||
|
for i=1:M.endo_nbr
|
||||||
|
disp(sprintf('%16s %12.6f %12.6f',M.endo_names(i,:),ys1(i), ys(i)))
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
function [resid,dr] = risky_residuals(ys,pm,M,dr,options,oo)
|
||||||
|
|
||||||
|
lead_lag_incidence = M.lead_lag_incidence;
|
||||||
|
iyv = lead_lag_incidence';
|
||||||
|
iyv = iyv(:);
|
||||||
|
iyr0 = find(iyv) ;
|
||||||
|
|
||||||
|
if M.exo_nbr == 0
|
||||||
|
oo.exo_steady_state = [] ;
|
||||||
|
end
|
||||||
|
|
||||||
|
z = repmat(ys,1,3);
|
||||||
|
z = z(iyr0) ;
|
||||||
|
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
||||||
|
[oo.exo_simul ...
|
||||||
|
oo.exo_det_simul], M.params, dr.ys, 2);
|
||||||
|
if ~isreal(d1) || ~isreal(d2)
|
||||||
|
pause
|
||||||
|
end
|
||||||
|
|
||||||
|
if options.use_dll
|
||||||
|
% In USE_DLL mode, the hessian is in the 3-column sparse representation
|
||||||
|
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
||||||
|
size(d1, 1), size(d1, 2)*size(d1, 2));
|
||||||
|
end
|
||||||
|
|
||||||
|
if isfield(options,'portfolio') && options.portfolio == 1
|
||||||
|
pm = portfolio_model_structure(M,options);
|
||||||
|
x = ys(pm.v_p);
|
||||||
|
dr = first_step_ds(x,pm,M,dr,options,oo);
|
||||||
|
dr.ys = ys;
|
||||||
|
else
|
||||||
|
d1a = d1;
|
||||||
|
d2a = d2;
|
||||||
|
end
|
||||||
|
|
||||||
|
gu1 = dr.ghu(pm.i_fwrd_g,:);
|
||||||
|
|
||||||
|
resid = resid1+0.5*(d1(:,pm.i_fwrd_f1)*dr.ghuu(pm.i_fwrd_g,:)+ ...
|
||||||
|
d2(:,pm.i_fwrd_f2)*kron(gu1,gu1))*vec(M.Sigma_e);
|
||||||
|
end
|
||||||
|
|
||||||
|
function [resid,dr] = risky_residuals_ds(x,pm,M,dr,options,oo)
|
||||||
|
|
||||||
|
v_p = pm.v_p;
|
||||||
|
v_np = pm.v_np;
|
||||||
|
|
||||||
|
% computing steady state of non-portfolio variables consistent with
|
||||||
|
% assumed portfolio
|
||||||
|
dr.ys(v_p) = x;
|
||||||
|
ys0 = dr.ys(v_np);
|
||||||
|
f_h =str2func([M.fname '_static']);
|
||||||
|
[dr.ys(v_np),info] = csolve(@ds_static_model,ys0,[],1e-10,100,f_h,x,pm.eq_np,v_np,v_p, ...
|
||||||
|
M.endo_nbr,M.exo_nbr,M.params);
|
||||||
|
if info
|
||||||
|
error('can''t compute non-portfolio steady state')
|
||||||
|
end
|
||||||
|
|
||||||
|
dr_np = first_step_ds(x,pm,M,dr,options,oo);
|
||||||
|
|
||||||
|
lead_lag_incidence = M.lead_lag_incidence;
|
||||||
|
iyv = lead_lag_incidence';
|
||||||
|
iyv = iyv(:);
|
||||||
|
iyr0 = find(iyv) ;
|
||||||
|
|
||||||
|
z = repmat(dr.ys,1,3);
|
||||||
|
z = z(iyr0) ;
|
||||||
|
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
||||||
|
[oo.exo_simul ...
|
||||||
|
oo.exo_det_simul], M.params, dr.ys, 2);
|
||||||
|
if ~isreal(d1) || ~isreal(d2)
|
||||||
|
pause
|
||||||
|
end
|
||||||
|
|
||||||
|
if options.use_dll
|
||||||
|
% In USE_DLL mode, the hessian is in the 3-column sparse representation
|
||||||
|
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
||||||
|
size(d1, 1), size(d1, 2)*size(d1, 2));
|
||||||
|
end
|
||||||
|
|
||||||
|
|
||||||
|
gu1 = dr_np.ghu(pm.i_fwrd_g,:);
|
||||||
|
|
||||||
|
resid = resid1+0.5*(d2(:,pm.i_fwrd_f2)*kron(gu1,gu1))*vec(M.Sigma_e);
|
||||||
|
|
||||||
|
resid = resid(pm.eq_p)
|
||||||
|
end
|
||||||
|
|
||||||
|
function dr_np = first_step_ds(x,pm,M,dr,options,oo)
|
||||||
|
|
||||||
|
lead_lag_incidence = M.lead_lag_incidence;
|
||||||
|
iyv = lead_lag_incidence';
|
||||||
|
iyv = iyv(:);
|
||||||
|
iyr0 = find(iyv) ;
|
||||||
|
|
||||||
|
ys = dr.ys;
|
||||||
|
ys(pm.v_p) = x;
|
||||||
|
|
||||||
|
z = repmat(ys,1,3);
|
||||||
|
z = z(iyr0) ;
|
||||||
|
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
||||||
|
[oo.exo_simul ...
|
||||||
|
oo.exo_det_simul], M.params, dr.ys, 2);
|
||||||
|
if ~isreal(d1) || ~isreal(d2)
|
||||||
|
pause
|
||||||
|
end
|
||||||
|
|
||||||
|
if options.use_dll
|
||||||
|
% In USE_DLL mode, the hessian is in the 3-column sparse representation
|
||||||
|
d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
||||||
|
size(d1, 1), size(d1, 2)*size(d1, 2));
|
||||||
|
end
|
||||||
|
|
||||||
|
d1_np = d1(pm.eq_np,pm.i_d1_np);
|
||||||
|
d2_np = d2(pm.eq_np,pm.i_d2_np);
|
||||||
|
|
||||||
|
[dr_np,info] = dyn_first_order_solver(d1_np,pm.M_np,pm.dr_np,options,0);
|
||||||
|
if info
|
||||||
|
print_info(info,0);
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
dr = dyn_second_order_solver(d1a,d2a,dr,M,...
|
dr_np = dyn_second_order_solver(d1_np,d2_np,dr_np,pm.M_np,...
|
||||||
options.threads.kronecker.A_times_B_kronecker_C,...
|
options.threads.kronecker.A_times_B_kronecker_C,...
|
||||||
options.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
options.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||||
end
|
end
|
||||||
|
|
||||||
function [resid,dr] = risky_residuals_k_order(ys,M,dr,options,oo)
|
function [resid,dr] = risky_residuals_k_order(ys,pm,M,dr,options,oo)
|
||||||
|
|
||||||
lead_lag_incidence = M.lead_lag_incidence;
|
|
||||||
npred = dr.npred;
|
|
||||||
exo_nbr = M.exo_nbr;
|
exo_nbr = M.exo_nbr;
|
||||||
vSigma_e = vec(M.Sigma_e);
|
endo_nbr = M.endo_nbr;
|
||||||
|
|
||||||
iyv = lead_lag_incidence';
|
iyv = M.lead_lag_incidence';
|
||||||
iyv = iyv(:);
|
iyv = iyv(:);
|
||||||
iyr0 = find(iyv) ;
|
iyr0 = find(iyv) ;
|
||||||
|
|
||||||
if M.exo_nbr == 0
|
if exo_nbr == 0
|
||||||
oo.exo_steady_state = [] ;
|
oo.exo_steady_state = [] ;
|
||||||
end
|
end
|
||||||
|
|
||||||
z = repmat(ys,1,3);
|
z = repmat(ys,1,3);
|
||||||
z = z(iyr0) ;
|
z = z(iyr0) ;
|
||||||
|
[resid1,d1,d2] = feval([M.fname '_dynamic'],z,...
|
||||||
|
[oo.exo_simul ...
|
||||||
|
oo.exo_det_simul], M.params, dr.ys, 2);
|
||||||
|
|
||||||
|
if isfield(options,'portfolio') && options.portfolio == 1
|
||||||
|
eq_np = pm.eq_np;
|
||||||
|
|
||||||
|
d1_np = d1(eq_np,pm.i_d1_np);
|
||||||
|
d2_np = d2(eq_np,pm.i_d2_np);
|
||||||
|
|
||||||
|
M_np = pm.M_np;
|
||||||
|
dr_np = pm.dr_np;
|
||||||
|
|
||||||
|
[dr_np,info] = dyn_first_order_solver(d1_np,pm.M_np,pm.dr_np,options,0);
|
||||||
|
if info
|
||||||
|
print_info(info,0);
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
dr_np = dyn_second_order_solver(d1_np,d2_np,dr_np,pm.M_np,...
|
||||||
|
options.threads.kronecker.A_times_B_kronecker_C,...
|
||||||
|
options.threads.kronecker.sparse_hessian_times_B_kronecker_C);
|
||||||
|
end
|
||||||
|
|
||||||
|
i_fwrd_f1 = pm.i_fwrd_f1;
|
||||||
|
i_fwrd_f2 = pm.i_fwrd_f2;
|
||||||
|
i_fwrd_f3 = pm.i_fwrd_f3;
|
||||||
|
i_fwrd_g = pm.i_fwrd_g;
|
||||||
|
gu1 = dr_np.ghu(i_fwrd_g,:);
|
||||||
|
ghuu = dr_np.ghuu;
|
||||||
|
|
||||||
|
resid = resid1+0.5*(d1(:,i_fwrd_f1)*ghuu(i_fwrd_g,:)+d2(:,i_fwrd_f2)* ...
|
||||||
|
kron(gu1,gu1))*vec(M.Sigma_e);
|
||||||
|
|
||||||
|
if nargout > 1
|
||||||
[resid1,d1,d2,d3] = feval([M.fname '_dynamic'],z,...
|
[resid1,d1,d2,d3] = feval([M.fname '_dynamic'],z,...
|
||||||
[oo.exo_simul ...
|
[oo.exo_simul ...
|
||||||
oo.exo_det_simul], M.params, dr.ys, 2);
|
oo.exo_det_simul], M.params, dr.ys, 2);
|
||||||
|
|
||||||
% $$$ hessian = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
|
||||||
% $$$ size(d1, 1), size(d1, 2)*size(d1, 2));
|
[a,b,c] = find(d2(eq_np,pm.i_d2_np));
|
||||||
% $$$ fy3 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
|
d2_np = [a b c];
|
||||||
% $$$ size(d1, 1), size(d1, 2)^3);
|
|
||||||
|
[a,b,c] = find(d3(eq_np,pm.i_d3_np));
|
||||||
|
d3_np = [a b c];
|
||||||
|
|
||||||
options.order = 3;
|
options.order = 3;
|
||||||
|
% space holder, unused by k_order_pertrubation
|
||||||
|
dr_np.ys = dr.ys(pm.v_np);
|
||||||
nu2 = exo_nbr*(exo_nbr+1)/2;
|
nu2 = exo_nbr*(exo_nbr+1)/2;
|
||||||
% $$$ d1_0 = d1;
|
nu3 = exo_nbr*(exo_nbr+1)*(exo_nbr+2)/3;
|
||||||
% $$$ gu1 = dr.ghu(dr.i_fwrd_g,:);
|
M_np.NZZDerivatives = [nnz(d1_np); nnz(d2_np); nnz(d3_np)];
|
||||||
% $$$ guu = dr.ghuu;
|
[err,g_0, g_1, g_2, g_3] = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np);
|
||||||
% $$$ for i=1:2
|
|
||||||
% $$$ d1 = d1_0 + 0.5*(hessian(:,dr.i_fwrd2a_f)*kron(eye(dr.nd),guu(dr.i_fwrd_g,:)*vSigma_e)+ ...
|
|
||||||
% $$$ fy3(:,dr.i_fwrd3_f)*kron(eye(dr.nd),kron(gu1,gu1)*vSigma_e));
|
|
||||||
% $$$ [junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
|
|
||||||
% $$$ b = zeros(M.endo_nbr,M.endo_nbr);
|
|
||||||
% $$$ b(:,cols_b) = d1(:,cols_j);
|
|
||||||
|
|
||||||
% $$$ [dr,info] = dyn_first_order_solver(d1,b,M,dr,options,0);
|
|
||||||
[err,g_0, g_1, g_2, g_3] = k_order_perturbation(dr,M,options);
|
|
||||||
mexErrCheck('k_order_perturbation', err);
|
mexErrCheck('k_order_perturbation', err);
|
||||||
gu1 = g_1(dr.i_fwrd_g,end-exo_nbr+1:end);
|
|
||||||
guu = unfold(g_2(:,end-nu2+1:end),exo_nbr);
|
|
||||||
d1old = d1;
|
|
||||||
% disp(max(max(abs(d1-d1old))));
|
|
||||||
% end
|
|
||||||
|
|
||||||
[junk,cols_b,cols_j] = find(lead_lag_incidence(2,dr.order_var));
|
gu1 = g_1(i_fwrd_g,end-exo_nbr+1:end);
|
||||||
|
ghuu = unfold2(g_2(:,end-nu2+1:end),exo_nbr);
|
||||||
|
ghsuu = get_ghsuu(g_3,size(g_1,2),exo_nbr);
|
||||||
|
|
||||||
resid = resid1+0.5*(d1(:,dr.i_fwrd_f)*guu(dr.i_fwrd_g,:)+hessian(:,dr.i_fwrd2_f)*kron(gu1,gu1))*vec(M.Sigma_e);
|
i_fwrd1_f2 = pm.i_fwrd1_f2;
|
||||||
|
i_fwrd1_f3 = pm.i_fwrd1_f3;
|
||||||
|
n = size(d1,2);
|
||||||
|
d1b = d1 + 0.5*( ...
|
||||||
|
d1(:,i_fwrd_f1)*...
|
||||||
|
d2(:,i_fwrd1_f2)*kron(eye(n),dr_np.ghuu(i_fwrd_g,:)*vec(M.Sigma_e))...
|
||||||
|
+ 0.5*d3(:,i_fwrd1_f3)*kron(eye(n),kron(gu1,gu1)*vec(M.Sigma_e)));
|
||||||
|
format short
|
||||||
|
kk1 = [nonzeros(M.lead_lag_incidence(:,1:6)'); ...
|
||||||
|
nnz(M.lead_lag_incidence)+[1; 2]]
|
||||||
|
kk2 = [nonzeros(M.lead_lag_incidence(:,1:6)'); ...
|
||||||
|
nnz(M.lead_lag_incidence)+[3; 4]]
|
||||||
|
format short
|
||||||
|
gu1
|
||||||
|
kron(gu1,gu1)*vec(M.Sigma_e)
|
||||||
|
disp(d1(:,:))
|
||||||
|
disp(d1b(:,:))
|
||||||
|
aa2=d2(:,i_fwrd1_f2)*kron(eye(n),dr_np.ghuu(i_fwrd_g,:)*vec(M.Sigma_e));
|
||||||
|
aa3=d3(:,i_fwrd1_f3)*kron(eye(n),kron(gu1,gu1)*vec(M.Sigma_e));
|
||||||
|
disp(d3(4,7+6*n+6*n*n))
|
||||||
|
disp(d3(4,8+16*n+17*n*n)) %8,17,18
|
||||||
|
disp(d3(4,8+17*n+16*n*n)) %8,17,18
|
||||||
|
disp(d3(4,7*n+17+17*n*n)) %8,17,18
|
||||||
|
disp(d3(4,7*n+18+16*n*n)) %8,17,18
|
||||||
|
disp(d3(4,7*n*n+16*n+18)) %8,17,18
|
||||||
|
disp(d3(4,7*n*n+17+17*n)) %8,17,18
|
||||||
|
pause
|
||||||
|
disp(aa2(:,kk1))
|
||||||
|
disp(aa2(:,kk2))
|
||||||
|
disp(aa3(:,kk1))
|
||||||
|
disp(aa3(:,kk2))
|
||||||
|
[dr,info] = dyn_first_order_solver(d1b,M,dr,options,0);
|
||||||
|
if info
|
||||||
|
print_info(info,0);
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
disp_dr(dr,dr.order_var,[]);
|
||||||
|
|
||||||
if nargout > 1
|
|
||||||
[dr,info] = k_order_pert(dr,M,options,oo);
|
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
function y=unfold(x,n)
|
function y=unfold2(x,n)
|
||||||
y = zeros(size(x,1),n*n);
|
y = zeros(size(x,1),n*n);
|
||||||
k = 1;
|
k = 1;
|
||||||
for i=1:n
|
for i=1:n
|
||||||
|
@ -526,6 +390,122 @@ function y=unfold(x,n)
|
||||||
if i ~= j
|
if i ~= j
|
||||||
y(:,(j-1)*n+i) = x(:,k);
|
y(:,(j-1)*n+i) = x(:,k);
|
||||||
end
|
end
|
||||||
|
k = k+1;
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
|
function y=unfold3(x,n)
|
||||||
|
y = zeros(size(x,1),n*n*n);
|
||||||
|
k = 1;
|
||||||
|
for i=1:n
|
||||||
|
for j=i:n
|
||||||
|
for m=j:n
|
||||||
|
y(:,(i-1)*n*n+(j-1)*n+m) = x(:,k);
|
||||||
|
y(:,(i-1)*n*n+(m-1)*n+j) = x(:,k);
|
||||||
|
y(:,(j-1)*n*n+(i-1)*n+m) = x(:,k);
|
||||||
|
y(:,(j-1)*n*n+(m-1)*n+i) = x(:,k);
|
||||||
|
y(:,(m-1)*n*n+(i-1)*n+j) = x(:,k);
|
||||||
|
y(:,(m-1)*n*n+(j-1)*n+i) = x(:,k);
|
||||||
|
|
||||||
|
k = k+1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function pm = portfolio_model_structure(M,options)
|
||||||
|
|
||||||
|
i_d3_np = [];
|
||||||
|
i_d3_p = [];
|
||||||
|
|
||||||
|
lead_index = M.maximum_endo_lag+2;
|
||||||
|
lead_lag_incidence = M.lead_lag_incidence;
|
||||||
|
eq_tags = M.equations_tags;
|
||||||
|
n_tags = size(eq_tags,1);
|
||||||
|
eq_p = cell2mat(eq_tags(strcmp(eq_tags(:,2), ...
|
||||||
|
'portfolio'),1));
|
||||||
|
pm.eq_p = eq_p;
|
||||||
|
pm.eq_np = setdiff(1:M.endo_nbr,eq_p);
|
||||||
|
v_p = zeros(n_tags,1);
|
||||||
|
for i=1:n_tags
|
||||||
|
v_p(i) = find(strncmp(eq_tags(i,3),M.endo_names, ...
|
||||||
|
length(cell2mat(eq_tags(i,3)))));
|
||||||
|
end
|
||||||
|
if any(lead_lag_incidence(lead_index,v_p))
|
||||||
|
error(['portfolio variables appear in the model as forward ' ...
|
||||||
|
'variable'])
|
||||||
|
end
|
||||||
|
pm.v_p = v_p;
|
||||||
|
v_np = setdiff(1:M.endo_nbr,v_p);
|
||||||
|
pm.v_np = v_np;
|
||||||
|
lli_np = lead_lag_incidence(:,v_np)';
|
||||||
|
k = find(lli_np);
|
||||||
|
lead_lag_incidence_np = lli_np;
|
||||||
|
lead_lag_incidence_np(k) = 1:nnz(lli_np);
|
||||||
|
lead_lag_incidence_np = lead_lag_incidence_np';
|
||||||
|
pm.lead_lag_incidence_np = lead_lag_incidence_np;
|
||||||
|
i_d1_np = [nonzeros(lli_np); nnz(lead_lag_incidence)+(1:M.exo_nbr)'];
|
||||||
|
pm.i_d1_np = i_d1_np;
|
||||||
|
|
||||||
|
n = nnz(lead_lag_incidence)+M.exo_nbr;
|
||||||
|
ih = reshape(1:n*n,n,n);
|
||||||
|
i_d2_np = ih(i_d1_np,i_d1_np);
|
||||||
|
pm.i_d2_np = i_d2_np(:);
|
||||||
|
|
||||||
|
ih = reshape(1:n*n*n,n,n,n);
|
||||||
|
i_d3_np = ih(i_d1_np,i_d1_np,i_d1_np);
|
||||||
|
pm.i_d3_np = i_d3_np(:);
|
||||||
|
|
||||||
|
M_np = M;
|
||||||
|
M_np.lead_lag_incidence = lead_lag_incidence_np;
|
||||||
|
M_np.lead_lag_incidence = lead_lag_incidence_np;
|
||||||
|
M_np.endo_nbr = length(v_np);
|
||||||
|
M_np.endo_names = M.endo_names(v_np,:);
|
||||||
|
dr_np = struct();
|
||||||
|
dr_np = set_state_space(dr_np,M_np,options);
|
||||||
|
pm.dr_np = dr_np;
|
||||||
|
M_np.var_order_endo_names = M_np.endo_names(dr_np.order_var,:);
|
||||||
|
pm.M_np = M_np;
|
||||||
|
pm.i_fwrd_g = find(lead_lag_incidence_np(lead_index,dr_np.order_var)');
|
||||||
|
|
||||||
|
i_fwrd_f1 = nonzeros(lead_lag_incidence(lead_index,:));
|
||||||
|
pm.i_fwrd_f1 = i_fwrd_f1;
|
||||||
|
n = nnz(lead_lag_incidence)+M.exo_nbr;
|
||||||
|
ih = reshape(1:n*n,n,n);
|
||||||
|
i_fwrd_f2 = ih(i_fwrd_f1,i_fwrd_f1);
|
||||||
|
pm.i_fwrd_f2 = i_fwrd_f2(:);
|
||||||
|
i_fwrd1_f2 = ih(i_fwrd_f1,:);
|
||||||
|
pm.i_fwrd1_f2 = i_fwrd1_f2(:);
|
||||||
|
|
||||||
|
ih = reshape(1:n*n*n,n,n,n);
|
||||||
|
i_fwrd_f3 = ih(i_fwrd_f1,i_fwrd_f1,i_fwrd_f1);
|
||||||
|
pm.i_fwrd_f3 = i_fwrd_f3(:);
|
||||||
|
i_fwrd1_f3 = ih(i_fwrd_f1,i_fwrd_f1,:);
|
||||||
|
pm.i_fwrd1_f3 = i_fwrd1_f3(:);
|
||||||
|
end
|
||||||
|
|
||||||
|
function r=ds_static_model(y0,f_h,p0,eq_np,v_np,v_p,endo_nbr,exo_nbr,params)
|
||||||
|
ys = zeros(endo_nbr,1);
|
||||||
|
ys(v_p) = p0;
|
||||||
|
ys(v_np) = y0;
|
||||||
|
r = f_h(ys,zeros(exo_nbr,1),params);
|
||||||
|
r = r(eq_np);
|
||||||
|
end
|
||||||
|
|
||||||
|
function ghsuu = get_ghsuu(g,ns,nx)
|
||||||
|
nxx = nx*(nx+1)/2;
|
||||||
|
m1 = 0;
|
||||||
|
m2 = ns*(ns+1)/2;
|
||||||
|
kk = 1:(nx*nx);
|
||||||
|
ghsuu = zeros(size(g,1),(ns*nx*nx));
|
||||||
|
|
||||||
|
for i=1:n
|
||||||
|
j = m1+(1:m2);
|
||||||
|
k = j(end-nxx+1:end);
|
||||||
|
ghsuu(:,kk) = unfold2(g(:,k),nx);
|
||||||
|
m1 = m1+m2;
|
||||||
|
m2 = m2 - (n-i+1);
|
||||||
|
kk = kk + nx*nx;
|
||||||
|
end
|
||||||
|
end
|
|
@ -74,7 +74,9 @@ if options_.solve_algo == 0
|
||||||
end
|
end
|
||||||
elseif options_.solve_algo == 1
|
elseif options_.solve_algo == 1
|
||||||
nn = size(x,1);
|
nn = size(x,1);
|
||||||
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,varargin{:});
|
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,1,options_.gstep, ...
|
||||||
|
options_.solve_tolf,options_.solve_tolx, ...
|
||||||
|
options_.solve_maxit,options_.debug,varargin{:});
|
||||||
elseif options_.solve_algo == 2 || options_.solve_algo == 4
|
elseif options_.solve_algo == 2 || options_.solve_algo == 4
|
||||||
nn = size(x,1) ;
|
nn = size(x,1) ;
|
||||||
tolf = options_.solve_tolf ;
|
tolf = options_.solve_tolf ;
|
||||||
|
@ -125,14 +127,19 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
|
||||||
if options_.debug
|
if options_.debug
|
||||||
disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
|
disp(['DYNARE_SOLVE (solve_algo=2|4): solving block ' num2str(i) ', of size ' num2str(r(i+1)-r(i)) ]);
|
||||||
end
|
end
|
||||||
[x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, bad_cond_flag, varargin{:});
|
[x,info]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag, ...
|
||||||
|
bad_cond_flag, options_.gstep, ...
|
||||||
|
options_.solve_tolf,options_.solve_tolx, ...
|
||||||
|
options_.solve_maxit,options_.debug,varargin{:});
|
||||||
if info
|
if info
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
fvec = feval(func,x,varargin{:});
|
fvec = feval(func,x,varargin{:});
|
||||||
if max(abs(fvec)) > tolf
|
if max(abs(fvec)) > tolf
|
||||||
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, varargin{:});
|
[x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag, bad_cond_flag, ...
|
||||||
|
options_.gstep, options_.solve_tolf,options_.solve_tolx, ...
|
||||||
|
options_.solve_maxit,options_.debug,varargin{:});
|
||||||
end
|
end
|
||||||
elseif options_.solve_algo == 3
|
elseif options_.solve_algo == 3
|
||||||
if jacobian_flag
|
if jacobian_flag
|
||||||
|
|
|
@ -34,7 +34,11 @@ function sim1_purely_backward()
|
||||||
yb = oo_.endo_simul(:,it-1); % Values at previous period, also used as guess value for current period
|
yb = oo_.endo_simul(:,it-1); % Values at previous period, also used as guess value for current period
|
||||||
yb1 = yb(iyb);
|
yb1 = yb(iyb);
|
||||||
|
|
||||||
tmp = solve1(model_dynamic, [yb1; yb], 1:M_.endo_nbr, nyb+1:nyb+M_.endo_nbr, 1, 1, oo_.exo_simul, M_.params, oo_.steady_state, it);
|
tmp = solve1(model_dynamic, [yb1; yb], 1:M_.endo_nbr, nyb+1:nyb+ ...
|
||||||
|
M_.endo_nbr, 1, 1, options_.gstep, ...
|
||||||
|
options_.solve_tolf,options_.solve_tolx, ...
|
||||||
|
options_.solve_maxit,options_.debug,oo_.exo_simul, ...
|
||||||
|
M_.params, oo_.steady_state, it);
|
||||||
oo_.endo_simul(:,it) = tmp(nyb+1:nyb+M_.endo_nbr);
|
oo_.endo_simul(:,it) = tmp(nyb+1:nyb+M_.endo_nbr);
|
||||||
end
|
end
|
||||||
|
|
|
@ -33,7 +33,11 @@ function sim1_purely_forward()
|
||||||
yf = oo_.endo_simul(:,it+1); % Values at next period, also used as guess value for current period
|
yf = oo_.endo_simul(:,it+1); % Values at next period, also used as guess value for current period
|
||||||
yf1 = yf(iyf);
|
yf1 = yf(iyf);
|
||||||
|
|
||||||
tmp = solve1(model_dynamic, [yf; yf1], 1:M_.endo_nbr, 1:M_.endo_nbr, 1, 1, oo_.exo_simul, M_.params, oo_.steady_state, it);
|
tmp = solve1(model_dynamic, [yf; yf1], 1:M_.endo_nbr, 1:M_.endo_nbr, ...
|
||||||
|
1, 1, options_.gstep, options_.solve_tolf, ...
|
||||||
|
options_.solve_tolx, options_.solve_maxit, ...
|
||||||
|
options_.debug,oo_.exo_simul, M_.params, oo_.steady_state, ...
|
||||||
|
it);
|
||||||
oo_.endo_simul(:,it) = tmp(1:M_.endo_nbr);
|
oo_.endo_simul(:,it) = tmp(1:M_.endo_nbr);
|
||||||
end
|
end
|
||||||
|
|
|
@ -104,6 +104,10 @@ DynareOutput.endo_simul(:,1) = DynareOutput.steady_state;
|
||||||
for it = 2:sample_size+1
|
for it = 2:sample_size+1
|
||||||
y(jdx) = DynareOutput.endo_simul(:,it-1); % A good guess for the initial conditions is the previous values for the endogenous variables.
|
y(jdx) = DynareOutput.endo_simul(:,it-1); % A good guess for the initial conditions is the previous values for the endogenous variables.
|
||||||
y(hdx) = y(jdx(iy1)); % Set lagged variables.
|
y(hdx) = y(jdx(iy1)); % Set lagged variables.
|
||||||
y(jdx) = solve1(model_dynamic, y, idx, jdx, 1, 1, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
|
y(jdx) = solve1(model_dynamic, y, idx, jdx, 1, 1, DynareOptions.gstep, ...
|
||||||
|
DynareOptions.solve_tolf,DynareOptions.solve_tolx, ...
|
||||||
|
DynareOptions.solve_maxit,DynareOptions.debug, ...
|
||||||
|
DynareOutput.exo_simul, DynareModel.params, ...
|
||||||
|
DynareOutput.steady_state, it);
|
||||||
DynareOutput.endo_simul(:,it) = y(jdx);
|
DynareOutput.endo_simul(:,it) = y(jdx);
|
||||||
end
|
end
|
|
@ -1,4 +1,4 @@
|
||||||
function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
|
function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,gstep,tolf,tolx,maxit,debug,varargin)
|
||||||
% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
|
% function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
|
||||||
% Solves systems of non linear equations of several variables
|
% Solves systems of non linear equations of several variables
|
||||||
%
|
%
|
||||||
|
@ -11,6 +11,12 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
|
||||||
% jacobian_flag=0: jacobian obtained numerically
|
% jacobian_flag=0: jacobian obtained numerically
|
||||||
% bad_cond_flag=1: when Jacobian is badly conditionned, use an
|
% bad_cond_flag=1: when Jacobian is badly conditionned, use an
|
||||||
% alternative formula to Newton step
|
% alternative formula to Newton step
|
||||||
|
% gstep increment multiplier in numercial derivative
|
||||||
|
% computation
|
||||||
|
% tolf tolerance for residuals
|
||||||
|
% tolx tolerance for solution variation
|
||||||
|
% maxit maximum number of iterations
|
||||||
|
% debug debug flag
|
||||||
% varargin: list of arguments following bad_cond_flag
|
% varargin: list of arguments following bad_cond_flag
|
||||||
%
|
%
|
||||||
% OUTPUTS
|
% OUTPUTS
|
||||||
|
@ -37,18 +43,13 @@ function [x,check] = solve1(func,x,j1,j2,jacobian_flag,bad_cond_flag,varargin)
|
||||||
% 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/>.
|
||||||
|
|
||||||
global M_ options_ fjac
|
|
||||||
|
|
||||||
nn = length(j1);
|
nn = length(j1);
|
||||||
fjac = zeros(nn,nn) ;
|
fjac = zeros(nn,nn) ;
|
||||||
g = zeros(nn,1) ;
|
g = zeros(nn,1) ;
|
||||||
|
|
||||||
tolf = options_.solve_tolf ;
|
|
||||||
tolx = options_.solve_tolx;
|
|
||||||
tolmin = tolx ;
|
tolmin = tolx ;
|
||||||
|
|
||||||
stpmx = 100 ;
|
stpmx = 100 ;
|
||||||
maxit = options_.solve_maxit ;
|
|
||||||
|
|
||||||
check = 0 ;
|
check = 0 ;
|
||||||
|
|
||||||
|
@ -77,7 +78,7 @@ for its = 1:maxit
|
||||||
fvec = fvec(j1);
|
fvec = fvec(j1);
|
||||||
fjac = fjac(j1,j2);
|
fjac = fjac(j1,j2);
|
||||||
else
|
else
|
||||||
dh = max(abs(x(j2)),options_.gstep(1)*ones(nn,1))*eps^(1/3);
|
dh = max(abs(x(j2)),gstep(1)*ones(nn,1))*eps^(1/3);
|
||||||
|
|
||||||
for j = 1:nn
|
for j = 1:nn
|
||||||
xdh = x ;
|
xdh = x ;
|
||||||
|
@ -89,33 +90,10 @@ for its = 1:maxit
|
||||||
end
|
end
|
||||||
|
|
||||||
g = (fvec'*fjac)';
|
g = (fvec'*fjac)';
|
||||||
if options_.debug
|
if debug
|
||||||
disp(['cond(fjac) ' num2str(cond(fjac))])
|
disp(['cond(fjac) ' num2str(cond(fjac))])
|
||||||
end
|
end
|
||||||
M_.unit_root = 0;
|
if bad_cond_flag && rcond(fjac) < sqrt(eps)
|
||||||
if M_.unit_root
|
|
||||||
if first_time
|
|
||||||
first_time = 0;
|
|
||||||
[q,r,e]=qr(fjac);
|
|
||||||
n = sum(abs(diag(r)) < 1e-12);
|
|
||||||
fvec = q'*fvec;
|
|
||||||
p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
|
|
||||||
disp(' ')
|
|
||||||
disp('STEADY with unit roots:')
|
|
||||||
disp(' ')
|
|
||||||
if n > 0
|
|
||||||
disp([' The following variable(s) kept their value given in INITVAL' ...
|
|
||||||
' or ENDVAL'])
|
|
||||||
disp(char(e(:,end-n+1:end)'*M_.endo_names))
|
|
||||||
else
|
|
||||||
disp(' STEADY can''t find any unit root!')
|
|
||||||
end
|
|
||||||
else
|
|
||||||
[q,r]=qr(fjac*e);
|
|
||||||
fvec = q'*fvec;
|
|
||||||
p = e*[-r(1:end-n,1:end-n)\fvec(1:end-n);zeros(n,1)];
|
|
||||||
end
|
|
||||||
elseif bad_cond_flag && cond(fjac) > 1/sqrt(eps)
|
|
||||||
fjac2=fjac'*fjac;
|
fjac2=fjac'*fjac;
|
||||||
p=-(fjac2+1e6*sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
|
p=-(fjac2+1e6*sqrt(nn*eps)*max(sum(abs(fjac2)))*eye(nn))\(fjac'*fvec);
|
||||||
else
|
else
|
||||||
|
@ -126,7 +104,7 @@ for its = 1:maxit
|
||||||
|
|
||||||
[x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:});
|
[x,f,fvec,check]=lnsrch1(xold,fold,g,p,stpmax,func,j1,j2,varargin{:});
|
||||||
|
|
||||||
if options_.debug
|
if debug
|
||||||
disp([its f])
|
disp([its f])
|
||||||
disp([xold x])
|
disp([xold x])
|
||||||
end
|
end
|
||||||
|
|
Loading…
Reference in New Issue