1function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ... 2 dr,options,oo) 3 4%@info: 5%! @deftypefn {Function File} {[@var{dr},@var{info}] =} dyn_risky_steadystate_solver (@var{ys0},@var{M},@var{dr},@var{options},@var{oo}) 6%! @anchor{dyn_risky_steadystate_solver} 7%! @sp 1 8%! Computes the second order risky steady state and first and second order reduced form of the DSGE model. 9%! @sp 2 10%! @strong{Inputs} 11%! @sp 1 12%! @table @ @var 13%! @item ys0 14%! Vector containing a guess value for the risky steady state 15%! @item M 16%! Matlab's structure describing the model (initialized by @code{dynare}). 17%! @item dr 18%! Matlab's structure describing the reduced form solution of the model. 19%! @item options 20%! Matlab's structure describing the options (initialized by @code{dynare}). 21%! @item oo 22%! Matlab's structure gathering the results (initialized by @code{dynare}). 23%! @end table 24%! @sp 2 25%! @strong{Outputs} 26%! @sp 1 27%! @table @ @var 28%! @item dr 29%! Matlab's structure describing the reduced form solution of the model. 30%! @item info 31%! Integer scalar, error code. 32%! @sp 1 33%! @table @ @code 34%! @item info==0 35%! No error. 36%! @item info==1 37%! The model doesn't determine the current variables uniquely. 38%! @item info==2 39%! MJDGGES returned an error code. 40%! @item info==3 41%! Blanchard & Kahn conditions are not satisfied: no stable equilibrium. 42%! @item info==4 43%! Blanchard & Kahn conditions are not satisfied: indeterminacy. 44%! @item info==5 45%! Blanchard & Kahn conditions are not satisfied: indeterminacy due to rank failure. 46%! @item info==6 47%! The jacobian evaluated at the deterministic steady state is complex. 48%! @item info==19 49%! The steadystate routine has thrown an exception (inconsistent deep parameters). 50%! @item info==20 51%! Cannot find the steady state, info(2) contains the sum of square residuals (of the static equations). 52%! @item info==21 53%! The steady state is complex, info(2) contains the sum of square of imaginary parts of the steady state. 54%! @item info==22 55%! The steady has NaNs. 56%! @item info==23 57%! M_.params has been updated in the steadystate routine and has complex valued scalars. 58%! @item info==24 59%! M_.params has been updated in the steadystate routine and has some NaNs. 60%! @end table 61%! @end table 62%! @end deftypefn 63%@eod: 64 65% Copyright (C) 2001-2018 Dynare Team 66% 67% This file is part of Dynare. 68% 69% Dynare is free software: you can redistribute it and/or modify 70% it under the terms of the GNU General Public License as published by 71% the Free Software Foundation, either version 3 of the License, or 72% (at your option) any later version. 73% 74% Dynare is distributed in the hope that it will be useful, 75% but WITHOUT ANY WARRANTY; without even the implied warranty of 76% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 77% GNU General Public License for more details. 78% 79% You should have received a copy of the GNU General Public License 80% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 81 82 83info = 0; 84lead_lag_incidence = M.lead_lag_incidence; 85order_var = dr.order_var; 86endo_nbr = M.endo_nbr; 87exo_nbr = M.exo_nbr; 88 89[~,dr.i_fwrd_g,i_fwrd_f] = find(lead_lag_incidence(3,order_var)); 90dr.i_fwrd_f = i_fwrd_f; 91nd = nnz(lead_lag_incidence) + M.exo_nbr; 92dr.nd = nd; 93kk = reshape(1:nd^2,nd,nd); 94kkk = reshape(1:nd^3,nd^2,nd); 95dr.i_fwrd2_f = kk(i_fwrd_f,i_fwrd_f); 96dr.i_fwrd2a_f = kk(i_fwrd_f,:); 97dr.i_fwrd3_f = kkk(dr.i_fwrd2_f,:); 98dr.i_uu = kk(end-exo_nbr+1:end,end-exo_nbr+1:end); 99if options.k_order_solver 100 func = @risky_residuals_k_order; 101else 102 func = @risky_residuals; 103end 104 105if isfield(options,'portfolio') && options.portfolio == 1 106 pm = portfolio_model_structure(M,options); 107 108 x0 = ys0(pm.v_p); 109 n = length(x0); 110 [x, info] = solve1(@risky_residuals_ds,x0,1:n,1:n,0,options.gstep, ... 111 options.solve_tolf,options.solve_tolx, ... 112 options.steady.maxit,options.debug,pm,M,dr, ... 113 options,oo); 114 if info 115 error('DS approach can''t be computed') 116 end 117 %[x, info] = csolve(@risky_residuals_ds,x0,[],1e-10,100,M,dr,options,oo); 118 % ys0(l_var) = x; 119 [resids,dr1] = risky_residuals_ds(x,pm,M,dr,options,oo); 120 ys1 = dr1.ys; 121else 122 pm = model_structure(M,options); 123end 124 125[ys, info] = solve1(func,ys0,1:endo_nbr,1:endo_nbr,0,options.gstep, ... 126 options.solve_tolf,options.solve_tolx, ... 127 options.steady.maxit,options.debug,pm,M,dr,options,oo); 128% [ys, info] = csolve(func,ys0,[],1e-10,100,M,dr,options,oo); 129if info 130 error('RSS approach can''t be computed') 131end 132dr.ys = ys; 133 134[resid,dr] = func(ys,pm,M,dr,options,oo); 135dr.ghs2 = zeros(M.endo_nbr,1); 136 137for i=1:M.endo_nbr 138 if isfield(options,'portfolio') && options.portfolio == 1 139 disp(sprintf('%16s %12.6f %12.6f', M.endo_names{i}, ys1(i), ... 140 ys(i))) 141 else 142 disp(sprintf('%16s %12.6f %12.6f', M.endo_names{i}, ys(i))) 143 end 144end 145 146end 147 148function [resid,dr] = risky_residuals(ys,pm,M,dr,options,oo) 149 150lead_lag_incidence = M.lead_lag_incidence; 151iyv = lead_lag_incidence'; 152iyv = iyv(:); 153iyr0 = find(iyv) ; 154 155if M.exo_nbr == 0 156 oo.exo_steady_state = [] ; 157end 158 159z = repmat(ys,1,3); 160z = z(iyr0) ; 161[resid1,d1,d2] = feval([M.fname '.dynamic'],z,... 162 [oo.exo_simul ... 163 oo.exo_det_simul], M.params, dr.ys, 2); 164if ~isreal(d1) || ~isreal(d2) 165 pause 166end 167 168if options.use_dll 169 % In USE_DLL mode, the hessian is in the 3-column sparse representation 170 d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ... 171 size(d1, 1), size(d1, 2)*size(d1, 2)); 172end 173 174if isfield(options,'portfolio') && options.portfolio == 1 175 pm = portfolio_model_structure(M,options); 176 x = ys(pm.v_p); 177 dr = first_step_ds(x,pm,M,dr,options,oo); 178 dr.ys = ys; 179else 180 pm = model_structure(M,options); 181 [dr,info] = dyn_first_order_solver(d1,M,dr,options,0); 182 if info 183 print_info(info,options.noprint,options); 184 end 185 dr = dyn_second_order_solver(d1,d2,dr,M,... 186 options.threads.kronecker.sparse_hessian_times_B_kronecker_C); 187end 188 189gu1 = dr.ghu(pm.i_fwrd_g,:); 190 191resid = resid1+0.5*(d1(:,pm.i_fwrd_f1)*dr.ghuu(pm.i_fwrd_g,:)+ ... 192 d2(:,pm.i_fwrd_f2)*kron(gu1,gu1))*vec(M.Sigma_e); 193end 194 195function [resid,dr] = risky_residuals_ds(x,pm,M,dr,options,oo) 196 197v_p = pm.v_p; 198v_np = pm.v_np; 199 200% computing steady state of non-portfolio variables consistent with 201% assumed portfolio 202dr.ys(v_p) = x; 203ys0 = dr.ys(v_np); 204f_h =str2func([M.fname '.static']); 205[dr.ys(v_np),info] = csolve(@ds_static_model,ys0,[],1e-10,100,f_h,x,pm.eq_np,v_np,v_p, ... 206 M.endo_nbr,M.exo_nbr,M.params); 207if info 208 error('can''t compute non-portfolio steady state') 209end 210 211dr_np = first_step_ds(x,pm,M,dr,options,oo); 212 213lead_lag_incidence = M.lead_lag_incidence; 214iyv = lead_lag_incidence'; 215iyv = iyv(:); 216iyr0 = find(iyv) ; 217 218z = repmat(dr.ys,1,3); 219z = z(iyr0) ; 220[resid1,d1,d2] = feval([M.fname '.dynamic'],z,... 221 [oo.exo_simul ... 222 oo.exo_det_simul], M.params, dr.ys, 2); 223if ~isreal(d1) || ~isreal(d2) 224 pause 225end 226 227if options.use_dll 228 % In USE_DLL mode, the hessian is in the 3-column sparse representation 229 d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ... 230 size(d1, 1), size(d1, 2)*size(d1, 2)); 231end 232 233 234gu1 = dr_np.ghu(pm.i_fwrd_g,:); 235 236resid = resid1+0.5*(d2(:,pm.i_fwrd_f2)*kron(gu1,gu1))*vec(M.Sigma_e); 237 238resid = resid(pm.eq_p) 239end 240 241function dr_np = first_step_ds(x,pm,M,dr,options,oo) 242 243lead_lag_incidence = M.lead_lag_incidence; 244iyv = lead_lag_incidence'; 245iyv = iyv(:); 246iyr0 = find(iyv) ; 247 248ys = dr.ys; 249ys(pm.v_p) = x; 250 251z = repmat(ys,1,3); 252z = z(iyr0) ; 253[resid1,d1,d2] = feval([M.fname '.dynamic'],z,... 254 [oo.exo_simul ... 255 oo.exo_det_simul], M.params, dr.ys, 2); 256if ~isreal(d1) || ~isreal(d2) 257 pause 258end 259 260if options.use_dll 261 % In USE_DLL mode, the hessian is in the 3-column sparse representation 262 d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ... 263 size(d1, 1), size(d1, 2)*size(d1, 2)); 264end 265 266d1_np = d1(pm.eq_np,pm.i_d1_np); 267d2_np = d2(pm.eq_np,pm.i_d2_np); 268 269[dr_np,info] = dyn_first_order_solver(d1_np,pm.M_np,pm.dr_np,options,0); 270if info 271 print_info(info, 0, options); 272 return 273end 274 275dr_np = dyn_second_order_solver(d1_np,d2_np,dr_np,pm.M_np,... 276 options.threads.kronecker.sparse_hessian_times_B_kronecker_C); 277end 278 279function [resid,dr] = risky_residuals_k_order(ys,pm,M,dr,options,oo) 280exo_nbr = M.exo_nbr; 281endo_nbr = M.endo_nbr; 282 283iyv = M.lead_lag_incidence'; 284iyv = iyv(:); 285iyr0 = find(iyv) ; 286 287if exo_nbr == 0 288 oo.exo_steady_state = [] ; 289end 290 291z = repmat(ys,1,3); 292z = z(iyr0) ; 293[resid1,d1,d2] = feval([M.fname '.dynamic'],z,... 294 [oo.exo_simul ... 295 oo.exo_det_simul], M.params, dr.ys, 2); 296 297if isfield(options,'portfolio') && options.portfolio == 1 298 eq_np = pm.eq_np; 299 300 d1_np = d1(eq_np,pm.i_d1_np); 301 d2_np = d2(eq_np,pm.i_d2_np); 302 303 M_np = pm.M_np; 304 dr_np = pm.dr_np; 305 306 [dr_np,info] = dyn_first_order_solver(d1_np,pm.M_np,pm.dr_np,options,0); 307 if info 308 print_info(info, 0, options); 309 return 310 end 311 312 dr_np = dyn_second_order_solver(d1_np,d2_np,dr_np,pm.M_np,... 313 options.threads.kronecker.sparse_hessian_times_B_kronecker_C); 314end 315 316i_fwrd_f1 = pm.i_fwrd_f1; 317i_fwrd_f2 = pm.i_fwrd_f2; 318i_fwrd_f3 = pm.i_fwrd_f3; 319i_fwrd_g = pm.i_fwrd_g; 320gu1 = dr_np.ghu(i_fwrd_g,:); 321ghuu = dr_np.ghuu; 322 323resid = resid1+0.5*(d1(:,i_fwrd_f1)*ghuu(i_fwrd_g,:)+d2(:,i_fwrd_f2)* ... 324 kron(gu1,gu1))*vec(M.Sigma_e); 325 326if nargout > 1 327 [resid1,d1,d2,d3] = feval([M.fname '.dynamic'],z,... 328 [oo.exo_simul ... 329 oo.exo_det_simul], M.params, dr.ys, 2); 330 331 332 [a,b,c] = find(d2(eq_np,pm.i_d2_np)); 333 d2_np = [a b c]; 334 335 [a,b,c] = find(d3(eq_np,pm.i_d3_np)); 336 d3_np = [a b c]; 337 338 options.order = 3; 339 % space holder, unused by k_order_pertrubation 340 dr_np.ys = dr.ys(pm.v_np); 341 nu2 = exo_nbr*(exo_nbr+1)/2; 342 nu3 = exo_nbr*(exo_nbr+1)*(exo_nbr+2)/3; 343 M_np.NZZDerivatives = [nnz(d1_np); nnz(d2_np); nnz(d3_np)]; 344 [err, dynpp_derivs] = k_order_perturbation(dr_np,M_np,options,d1_np,d2_np,d3_np); 345 mexErrCheck('k_order_perturbation', err); 346 g_0 = dynpp_derivs.g_0; 347 g_1 = dynpp_derivs.g_1; 348 g_2 = dynpp_derivs.g_2; 349 g_3 = dynpp_derivs.g_3; 350 351 gu1 = g_1(i_fwrd_g,end-exo_nbr+1:end); 352 ghuu = unfold2(g_2(:,end-nu2+1:end),exo_nbr); 353 ghsuu = get_ghsuu(g_3,size(g_1,2),exo_nbr); 354 355 i_fwrd1_f2 = pm.i_fwrd1_f2; 356 i_fwrd1_f3 = pm.i_fwrd1_f3; 357 n = size(d1,2); 358 d1b = d1 + 0.5*( ... 359 d1(:,i_fwrd_f1)*... 360 d2(:,i_fwrd1_f2)*kron(eye(n),dr_np.ghuu(i_fwrd_g,:)*vec(M.Sigma_e))... 361 + 0.5*d3(:,i_fwrd1_f3)*kron(eye(n),kron(gu1,gu1)*vec(M.Sigma_e))); 362 format short 363 kk1 = [nonzeros(M.lead_lag_incidence(:,1:6)'); ... 364 nnz(M.lead_lag_incidence)+[1; 2]] 365 kk2 = [nonzeros(M.lead_lag_incidence(:,1:6)'); ... 366 nnz(M.lead_lag_incidence)+[3; 4]] 367 format short 368 gu1 369 kron(gu1,gu1)*vec(M.Sigma_e) 370 disp(d1(:,:)) 371 disp(d1b(:,:)) 372 aa2=d2(:,i_fwrd1_f2)*kron(eye(n),dr_np.ghuu(i_fwrd_g,:)*vec(M.Sigma_e)); 373 aa3=d3(:,i_fwrd1_f3)*kron(eye(n),kron(gu1,gu1)*vec(M.Sigma_e)); 374 disp(d3(4,7+6*n+6*n*n)) 375 disp(d3(4,8+16*n+17*n*n)) %8,17,18 376 disp(d3(4,8+17*n+16*n*n)) %8,17,18 377 disp(d3(4,7*n+17+17*n*n)) %8,17,18 378 disp(d3(4,7*n+18+16*n*n)) %8,17,18 379 disp(d3(4,7*n*n+16*n+18)) %8,17,18 380 disp(d3(4,7*n*n+17+17*n)) %8,17,18 381 pause 382 disp(aa2(:,kk1)) 383 disp(aa2(:,kk2)) 384 disp(aa3(:,kk1)) 385 disp(aa3(:,kk2)) 386 [dr,info] = dyn_first_order_solver(d1b,M,dr,options,0); 387 if info 388 print_info(info, 0, options); 389 return 390 end 391 392 disp_dr(dr,dr.order_var,[]); 393 394end 395end 396 397function y=unfold2(x,n) 398y = zeros(size(x,1),n*n); 399k = 1; 400for i=1:n 401 for j=i:n 402 y(:,(i-1)*n+j) = x(:,k); 403 if i ~= j 404 y(:,(j-1)*n+i) = x(:,k); 405 end 406 k = k+1; 407 end 408end 409end 410 411function y=unfold3(x,n) 412y = zeros(size(x,1),n*n*n); 413k = 1; 414for i=1:n 415 for j=i:n 416 for m=j:n 417 y(:,(i-1)*n*n+(j-1)*n+m) = x(:,k); 418 y(:,(i-1)*n*n+(m-1)*n+j) = x(:,k); 419 y(:,(j-1)*n*n+(i-1)*n+m) = x(:,k); 420 y(:,(j-1)*n*n+(m-1)*n+i) = x(:,k); 421 y(:,(m-1)*n*n+(i-1)*n+j) = x(:,k); 422 y(:,(m-1)*n*n+(j-1)*n+i) = x(:,k); 423 424 k = k+1; 425 end 426 end 427end 428end 429 430function pm = model_structure(M,options) 431 432 433lead_index = M.maximum_endo_lag+2; 434lead_lag_incidence = M.lead_lag_incidence; 435dr = struct(); 436dr = set_state_space(dr,M,options); 437pm.i_fwrd_g = find(lead_lag_incidence(lead_index,dr.order_var)'); 438 439i_fwrd_f1 = nonzeros(lead_lag_incidence(lead_index,dr.order_var)); 440pm.i_fwrd_f1 = i_fwrd_f1; 441n = nnz(lead_lag_incidence)+M.exo_nbr; 442ih = reshape(1:n*n,n,n); 443i_fwrd_f2 = ih(i_fwrd_f1,i_fwrd_f1); 444pm.i_fwrd_f2 = i_fwrd_f2(:); 445i_fwrd1_f2 = ih(i_fwrd_f1,:); 446pm.i_fwrd1_f2 = i_fwrd1_f2(:); 447 448ih = reshape(1:n*n*n,n,n,n); 449i_fwrd_f3 = ih(i_fwrd_f1,i_fwrd_f1,i_fwrd_f1); 450pm.i_fwrd_f3 = i_fwrd_f3(:); 451i_fwrd1_f3 = ih(i_fwrd_f1,i_fwrd_f1,:); 452pm.i_fwrd1_f3 = i_fwrd1_f3(:); 453end 454 455function pm = portfolio_model_structure(M,options) 456 457i_d3_np = []; 458i_d3_p = []; 459 460lead_index = M.maximum_endo_lag+2; 461lead_lag_incidence = M.lead_lag_incidence; 462eq_tags = M.equations_tags; 463n_tags = size(eq_tags,1); 464eq_p = cell2mat(eq_tags(strcmp(eq_tags(:,2), ... 465 'portfolio'),1)); 466pm.eq_p = eq_p; 467pm.eq_np = setdiff(1:M.endo_nbr,eq_p); 468v_p = zeros(n_tags,1); 469for i=1:n_tags 470 v_p(i) = find(strncmp(eq_tags(i,3),M.endo_names, ... 471 length(cell2mat(eq_tags(i,3))))); 472end 473if any(lead_lag_incidence(lead_index,v_p)) 474 error(['portfolio variables appear in the model as forward ' ... 475 'variable']) 476end 477pm.v_p = v_p; 478v_np = setdiff(1:M.endo_nbr,v_p); 479pm.v_np = v_np; 480lli_np = lead_lag_incidence(:,v_np)'; 481k = find(lli_np); 482lead_lag_incidence_np = lli_np; 483lead_lag_incidence_np(k) = 1:nnz(lli_np); 484lead_lag_incidence_np = lead_lag_incidence_np'; 485pm.lead_lag_incidence_np = lead_lag_incidence_np; 486i_d1_np = [nonzeros(lli_np); nnz(lead_lag_incidence)+(1:M.exo_nbr)']; 487pm.i_d1_np = i_d1_np; 488 489n = nnz(lead_lag_incidence)+M.exo_nbr; 490ih = reshape(1:n*n,n,n); 491i_d2_np = ih(i_d1_np,i_d1_np); 492pm.i_d2_np = i_d2_np(:); 493 494ih = reshape(1:n*n*n,n,n,n); 495i_d3_np = ih(i_d1_np,i_d1_np,i_d1_np); 496pm.i_d3_np = i_d3_np(:); 497 498M_np = M; 499M_np.lead_lag_incidence = lead_lag_incidence_np; 500M_np.lead_lag_incidence = lead_lag_incidence_np; 501M_np.endo_nbr = length(v_np); 502M_np.endo_names = M.endo_names(v_np); 503dr_np = struct(); 504dr_np = set_state_space(dr_np,M_np,options); 505pm.dr_np = dr_np; 506pm.M_np = M_np; 507pm.i_fwrd_g = find(lead_lag_incidence_np(lead_index,dr_np.order_var)'); 508 509i_fwrd_f1 = nonzeros(lead_lag_incidence(lead_index,:)); 510pm.i_fwrd_f1 = i_fwrd_f1; 511n = nnz(lead_lag_incidence)+M.exo_nbr; 512ih = reshape(1:n*n,n,n); 513i_fwrd_f2 = ih(i_fwrd_f1,i_fwrd_f1); 514pm.i_fwrd_f2 = i_fwrd_f2(:); 515i_fwrd1_f2 = ih(i_fwrd_f1,:); 516pm.i_fwrd1_f2 = i_fwrd1_f2(:); 517 518ih = reshape(1:n*n*n,n,n,n); 519i_fwrd_f3 = ih(i_fwrd_f1,i_fwrd_f1,i_fwrd_f1); 520pm.i_fwrd_f3 = i_fwrd_f3(:); 521i_fwrd1_f3 = ih(i_fwrd_f1,i_fwrd_f1,:); 522pm.i_fwrd1_f3 = i_fwrd1_f3(:); 523end 524 525function r=ds_static_model(y0,f_h,p0,eq_np,v_np,v_p,endo_nbr,exo_nbr,params) 526ys = zeros(endo_nbr,1); 527ys(v_p) = p0; 528ys(v_np) = y0; 529r = f_h(ys,zeros(exo_nbr,1),params); 530r = r(eq_np); 531end 532 533function ghsuu = get_ghsuu(g,ns,nx) 534nxx = nx*(nx+1)/2; 535m1 = 0; 536m2 = ns*(ns+1)/2; 537kk = 1:(nx*nx); 538ghsuu = zeros(size(g,1),(ns*nx*nx)); 539 540for i=1:n 541 j = m1+(1:m2); 542 k = j(end-nxx+1:end); 543 ghsuu(:,kk) = unfold2(g(:,k),nx); 544 m1 = m1+m2; 545 m2 = m2 - (n-i+1); 546 kk = kk + nx*nx; 547end 548end 549