From 4b83c1bf763fe5029c94b7e55aba257135426de1 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sun, 6 Sep 2015 11:41:49 +0200 Subject: [PATCH 1/2] Integrates Tom Holden's robust linear solver Supersedes #984 --- doc/dynare.texi | 3 + matlab/global_initialization.m | 1 + matlab/perfect-foresight-models/sim1.m | 85 +++++++++++++++++- tests/Makefile.am | 3 +- tests/simul/Irreversible_investment.mod | 112 ++++++++++++++++++++++++ 5 files changed, 201 insertions(+), 3 deletions(-) create mode 100644 tests/simul/Irreversible_investment.mod diff --git a/doc/dynare.texi b/doc/dynare.texi index c06b434c0..ae366df38 100644 --- a/doc/dynare.texi +++ b/doc/dynare.texi @@ -2879,6 +2879,9 @@ add new variables one by one. Determines the maximum number of iterations used in the non-linear solver. The default value of @code{maxit} is 50. +@item robust_lin_solve +Triggers the use of a robust linear solver for the default @code{solve_algo=4}. + @item solve_algo = @var{INTEGER} @anchor{solve_algo} Determines the non-linear solver to use. Possible values for the option are: diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index c3e70fac4..80d9e5731 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -63,6 +63,7 @@ options_.minimal_workspace = 0; options_.dp.maxit = 3000; options_.steady.maxit = 50; options_.simul.maxit = 50; +options_.simul.robust_lin_solve = 0; options_.mode_check.status = 0; options_.mode_check.neighbourhood_size = .5; diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m index 13e404f22..3fecd13c8 100644 --- a/matlab/perfect-foresight-models/sim1.m +++ b/matlab/perfect-foresight-models/sim1.m @@ -162,9 +162,17 @@ for iter = 1:options.simul.maxit if endogenous_terminal_period && iter>1 dy = ZERO; - dy(1:i_rows(end)) = -A(1:i_rows(end),1:i_rows(end))\res(1:i_rows(end)); + if options.simul.robust_lin_solve + dy(1:i_rows(end)) = -lin_solve_robust( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)) ); + else + dy(1:i_rows(end)) = -lin_solve( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)) ); + end else - dy = -A\res; + if options.simul.robust_lin_solve + dy = -lin_solve_robust( A, res ); + else + dy = -lin_solve( A, res ); + end end Y(i_upd) = Y(i_upd) + dy; @@ -223,3 +231,76 @@ end if verbose skipline(); end + +function x = lin_solve( A, b ) + if norm( b ) < sqrt( eps ) % then x = 0 is a solution + x = 0; + return + end + + x = A\b; + x( ~isfinite( x ) ) = 0; + relres = norm( b - A * x ) / norm( b ); + if relres > 1e-6 + fprintf( 'WARNING : Failed to find a solution to the linear system.\n' ); + end + +function [ x, flag, relres ] = lin_solve_robust( A, b ) + if norm( b ) < sqrt( eps ) % then x = 0 is a solution + x = 0; + flag = 0; + relres = 0; + return + end + + x = A\b; + x( ~isfinite( x ) ) = 0; + [ x, flag, relres ] = bicgstab( A, b, [], [], [], [], x ); % returns immediately if x is a solution + if flag == 0 + return + end + + disp( relres ); + + fprintf( 'Initial bicgstab failed, trying alternative start point.\n' ); + old_x = x; + old_relres = relres; + [ x, flag, relres ] = bicgstab( A, b ); + if flag == 0 + return + end + + fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' ); + if old_relres < relres + x = old_x; + end + [ x, flag, relres ] = gmres( A, b, [], [], [], [], [], x ); + if flag == 0 + return + end + + fprintf( 'Initial gmres failed, trying alternative start point.\n' ); + old_x = x; + old_relres = relres; + [ x, flag, relres ] = gmres( A, b ); + if flag == 0 + return + end + + fprintf( 'Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' ); + if old_relres < relres + x = old_x; + relres = old_relres; + end + old_x = x; + old_relres = relres; + x = pinv( full( A ) ) * b; + relres = norm( b - A * x ) / norm( b ); + if old_relres < relres + x = old_x; + relres = old_relres; + end + flag = relres > 1e-6; + if flag ~= 0 + fprintf( 'WARNING : Failed to find a solution to the linear system\n' ); + end \ No newline at end of file diff --git a/tests/Makefile.am b/tests/Makefile.am index 6fce62bf7..065ff8cfc 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -149,7 +149,8 @@ MODFILES = \ simul/Solow_no_varexo.mod \ simul/simul_ZLB_purely_forward.mod \ simul/simul_ZLB_purely_forward_no_solution.mod \ - conditional_forecasts/1/fs2000_cal.mod \ + + simul/Irreversible_investment.mod \ conditional_forecasts/2/fs2000_est.mod \ conditional_forecasts/3/fs2000_conditional_forecast_initval.mod \ conditional_forecasts/4/fs2000_conditional_forecast_histval.mod \ diff --git a/tests/simul/Irreversible_investment.mod b/tests/simul/Irreversible_investment.mod new file mode 100644 index 000000000..02719c427 --- /dev/null +++ b/tests/simul/Irreversible_investment.mod @@ -0,0 +1,112 @@ +@#define NumberOfCountries = 2 + +var log_kappa; +@#for Country in 1:NumberOfCountries + var mu@{Country}, logit_l@{Country}, k@{Country}, a@{Country}; +@#endfor + +parameters alpha beta varrho delta theta rhoA sigmaA; + +alpha = 0.3; +beta = 0.99; +varrho = 1.5; +delta = 0.025; +theta = 0.9; +rhoA = 0.95; +sigmaA = 0.05; + +@#for Country in 1:NumberOfCountries + varexo epsilonA@{Country}; +@#endfor + +model; + + #kappa = exp( log_kappa ); + #LEAD_kappa = exp( log_kappa(+1) ); + + #min_A = theta ^ ( 1/alpha ); + #mean_a = log( 1 - min_A ); + + @#for Country in 1:NumberOfCountries + + #K@{Country} = exp( k@{Country} ); + #LAG_K@{Country} = exp( k@{Country}(-1) ); + #LEAD_A@{Country} = min_A + exp( a@{Country}(+1) ); + #A@{Country} = min_A + exp( a@{Country} ); + #LAG_A@{Country} = min_A + exp( a@{Country}(-1) ); + #L@{Country} = 1 / ( 1 + exp( -logit_l@{Country} ) ); + #LEAD_L@{Country} = 1 / ( 1 + exp( -logit_l@{Country}(+1) ) ); + + #C@{Country} = kappa ^ ( -1/varrho ) - theta / ( 1-alpha ) * ( A@{Country}*(1-L@{Country}) ) ^ ( 1-alpha ); + #LEAD_phi@{Country} = ( 1 - theta * ( LEAD_A@{Country}*(1-LEAD_L@{Country}) ) ^ ( -alpha ) ) * LEAD_kappa; + #I@{Country} = K@{Country} - ( 1-delta ) * LAG_K@{Country}; + + @#endfor + + @#for Country in 1:NumberOfCountries + + ( a@{Country} - mean_a ) = rhoA * ( a@{Country}(-1) - mean_a ) - sigmaA * epsilonA@{Country}; + L@{Country} = min( LAG_K@{Country} / A@{Country}, 1 - theta ^ ( 1/alpha ) / A@{Country} ); + kappa - mu@{Country} = beta * ( ( 1-delta ) * ( LEAD_kappa - mu@{Country}(+1) ) + LEAD_phi@{Country} ); + kappa = max( beta * ( ( 1-delta ) * ( LEAD_kappa - mu@{Country}(+1) ) + LEAD_phi@{Country} ), ( theta / ( 1-alpha ) * ( A@{Country}*(1-L@{Country}) ) ^ ( 1-alpha ) + @#for OtherCountry in 1:NumberOfCountries + + A@{OtherCountry} * L@{OtherCountry} + @#if OtherCountry != Country + - C@{OtherCountry} - I@{OtherCountry} + @#endif + @#endfor + ) ^ ( -varrho ) ); + + @#endfor + + 0 = 0 + @#for OtherCountry in 1:NumberOfCountries + + A@{OtherCountry} * L@{OtherCountry} + - C@{OtherCountry} - I@{OtherCountry} + @#endfor + ; + +end; + +steady_state_model; + + a = log( 1 - theta ^ ( 1/alpha ) ); + mu = 0; + L = 1 - ( 1/theta * ( 2 - 1/beta - delta ) ) ^ ( -1/alpha ); + K = L; + I = delta * K; + C = ( 1 - delta ) * L; + + kappa_ = ( C + theta / ( 1-alpha ) * ( 1-L ) ^ ( 1-alpha ) ) ^ ( -varrho ); + + log_kappa = log( kappa_ ); + + @#for Country in 1:NumberOfCountries + + mu@{Country} = mu; + logit_l@{Country} = log( L / ( 1 - L ) ); + k@{Country} = log( K ); + a@{Country} = a; + + @#endfor + +end; + +steady; +check; + +shocks; + +var epsilonA1; periods 1; values 2; +@#for Country in 2:NumberOfCountries + var epsilonA@{Country}; periods 1; values 0; +@#endfor + +end; + +options_.simul.robust_lin_solve=1; +simul( periods = 400 ); + +if ~oo_.deterministic_simulation.status + error('Model did not solve') +end \ No newline at end of file From 0669ac4c1c3d2d7c90cc64460cb75e376f5e2e34 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Fri, 25 Mar 2016 20:14:25 +0100 Subject: [PATCH 2/2] Adjust verbosity of solvers --- .../perfect_foresight_solver.m | 2 +- matlab/perfect-foresight-models/sim1.m | 32 ++++++++++++------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver.m b/matlab/perfect-foresight-models/perfect_foresight_solver.m index 464bfff24..09b3bf435 100644 --- a/matlab/perfect-foresight-models/perfect_foresight_solver.m +++ b/matlab/perfect-foresight-models/perfect_foresight_solver.m @@ -176,7 +176,7 @@ end if oo_.deterministic_simulation.status == 1 disp('Perfect foresight solution found.') else - warning('Failed to solve perfect foresight model') + disp('Failed to solve perfect foresight model') end skipline() diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m index 3fecd13c8..0f92ca8ed 100644 --- a/matlab/perfect-foresight-models/sim1.m +++ b/matlab/perfect-foresight-models/sim1.m @@ -163,15 +163,15 @@ for iter = 1:options.simul.maxit if endogenous_terminal_period && iter>1 dy = ZERO; if options.simul.robust_lin_solve - dy(1:i_rows(end)) = -lin_solve_robust( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)) ); + dy(1:i_rows(end)) = -lin_solve_robust( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)),verbose ); else - dy(1:i_rows(end)) = -lin_solve( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)) ); + dy(1:i_rows(end)) = -lin_solve( A(1:i_rows(end),1:i_rows(end)), res(1:i_rows(end)), verbose ); end else if options.simul.robust_lin_solve - dy = -lin_solve_robust( A, res ); + dy = -lin_solve_robust( A, res, verbose ); else - dy = -lin_solve( A, res ); + dy = -lin_solve( A, res, verbose ); end end @@ -232,7 +232,7 @@ if verbose skipline(); end -function x = lin_solve( A, b ) +function x = lin_solve( A, b,verbose) if norm( b ) < sqrt( eps ) % then x = 0 is a solution x = 0; return @@ -241,11 +241,11 @@ function x = lin_solve( A, b ) x = A\b; x( ~isfinite( x ) ) = 0; relres = norm( b - A * x ) / norm( b ); - if relres > 1e-6 + if relres > 1e-6 && verbose fprintf( 'WARNING : Failed to find a solution to the linear system.\n' ); end -function [ x, flag, relres ] = lin_solve_robust( A, b ) +function [ x, flag, relres ] = lin_solve_robust( A, b , verbose) if norm( b ) < sqrt( eps ) % then x = 0 is a solution x = 0; flag = 0; @@ -262,7 +262,9 @@ function [ x, flag, relres ] = lin_solve_robust( A, b ) disp( relres ); - fprintf( 'Initial bicgstab failed, trying alternative start point.\n' ); + if verbose + fprintf( 'Initial bicgstab failed, trying alternative start point.\n' ); + end old_x = x; old_relres = relres; [ x, flag, relres ] = bicgstab( A, b ); @@ -270,7 +272,9 @@ function [ x, flag, relres ] = lin_solve_robust( A, b ) return end - fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' ); + if verbose + fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' ); + end if old_relres < relres x = old_x; end @@ -279,7 +283,9 @@ function [ x, flag, relres ] = lin_solve_robust( A, b ) return end - fprintf( 'Initial gmres failed, trying alternative start point.\n' ); + if verbose + fprintf( 'Initial gmres failed, trying alternative start point.\n' ); + end old_x = x; old_relres = relres; [ x, flag, relres ] = gmres( A, b ); @@ -287,7 +293,9 @@ function [ x, flag, relres ] = lin_solve_robust( A, b ) return end - fprintf( 'Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' ); + if verbose + fprintf( 'Alternative start point also failed with gmres, using the (SLOW) Moore-Penrose Pseudo-Inverse.\n' ); + end if old_relres < relres x = old_x; relres = old_relres; @@ -301,6 +309,6 @@ function [ x, flag, relres ] = lin_solve_robust( A, b ) relres = old_relres; end flag = relres > 1e-6; - if flag ~= 0 + if flag ~= 0 && verbose fprintf( 'WARNING : Failed to find a solution to the linear system\n' ); end \ No newline at end of file