diff --git a/doc/dynare.texi b/doc/dynare.texi
index c06b434c0c6d7207dac6e56e89b89a29e49dbd8f..ae366df3800b904f5b17e10918e1217c7f6ee2f6 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 c3e70fac4bd5153fd83fe3150142f96abc217f4a..80d9e573154d6f00c1c3273e4b67ac773a8ee997 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 13e404f226223b2b83718de27801d54e58c0a501..3fecd13c8a4062eeba7bbfd4a51251288263661e 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 6fce62bf7005308079f1fc8a992aa25930eba781..065ff8cfcbe0466916121a59cb030f7746a2ee05 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 0000000000000000000000000000000000000000..02719c4271edf6567a0867f3dcb3c4f28a113cb2
--- /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