Merge pull request #1052 from JohannesPfeifer/lin_solve_robust
Integrate robust linear solver for stack_solve_algo=0time-shift
commit
aa1a0b5e25
|
@ -2879,6 +2879,9 @@ add new variables one by one.
|
||||||
Determines the maximum number of iterations used in the non-linear solver. The
|
Determines the maximum number of iterations used in the non-linear solver. The
|
||||||
default value of @code{maxit} is 50.
|
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}
|
@item solve_algo = @var{INTEGER}
|
||||||
@anchor{solve_algo}
|
@anchor{solve_algo}
|
||||||
Determines the non-linear solver to use. Possible values for the option are:
|
Determines the non-linear solver to use. Possible values for the option are:
|
||||||
|
|
|
@ -63,6 +63,7 @@ options_.minimal_workspace = 0;
|
||||||
options_.dp.maxit = 3000;
|
options_.dp.maxit = 3000;
|
||||||
options_.steady.maxit = 50;
|
options_.steady.maxit = 50;
|
||||||
options_.simul.maxit = 50;
|
options_.simul.maxit = 50;
|
||||||
|
options_.simul.robust_lin_solve = 0;
|
||||||
|
|
||||||
options_.mode_check.status = 0;
|
options_.mode_check.status = 0;
|
||||||
options_.mode_check.neighbourhood_size = .5;
|
options_.mode_check.neighbourhood_size = .5;
|
||||||
|
|
|
@ -176,7 +176,7 @@ end
|
||||||
if oo_.deterministic_simulation.status == 1
|
if oo_.deterministic_simulation.status == 1
|
||||||
disp('Perfect foresight solution found.')
|
disp('Perfect foresight solution found.')
|
||||||
else
|
else
|
||||||
warning('Failed to solve perfect foresight model')
|
disp('Failed to solve perfect foresight model')
|
||||||
end
|
end
|
||||||
|
|
||||||
skipline()
|
skipline()
|
||||||
|
|
|
@ -162,9 +162,17 @@ for iter = 1:options.simul.maxit
|
||||||
|
|
||||||
if endogenous_terminal_period && iter>1
|
if endogenous_terminal_period && iter>1
|
||||||
dy = ZERO;
|
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)),verbose );
|
||||||
else
|
else
|
||||||
dy = -A\res;
|
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, verbose );
|
||||||
|
else
|
||||||
|
dy = -lin_solve( A, res, verbose );
|
||||||
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
Y(i_upd) = Y(i_upd) + dy;
|
Y(i_upd) = Y(i_upd) + dy;
|
||||||
|
@ -223,3 +231,84 @@ end
|
||||||
if verbose
|
if verbose
|
||||||
skipline();
|
skipline();
|
||||||
end
|
end
|
||||||
|
|
||||||
|
function x = lin_solve( A, b,verbose)
|
||||||
|
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 && verbose
|
||||||
|
fprintf( 'WARNING : Failed to find a solution to the linear system.\n' );
|
||||||
|
end
|
||||||
|
|
||||||
|
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;
|
||||||
|
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 );
|
||||||
|
|
||||||
|
if verbose
|
||||||
|
fprintf( 'Initial bicgstab failed, trying alternative start point.\n' );
|
||||||
|
end
|
||||||
|
old_x = x;
|
||||||
|
old_relres = relres;
|
||||||
|
[ x, flag, relres ] = bicgstab( A, b );
|
||||||
|
if flag == 0
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
if verbose
|
||||||
|
fprintf( 'Alternative start point also failed with bicgstab, trying gmres.\n' );
|
||||||
|
end
|
||||||
|
if old_relres < relres
|
||||||
|
x = old_x;
|
||||||
|
end
|
||||||
|
[ x, flag, relres ] = gmres( A, b, [], [], [], [], [], x );
|
||||||
|
if flag == 0
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
if verbose
|
||||||
|
fprintf( 'Initial gmres failed, trying alternative start point.\n' );
|
||||||
|
end
|
||||||
|
old_x = x;
|
||||||
|
old_relres = relres;
|
||||||
|
[ x, flag, relres ] = gmres( A, b );
|
||||||
|
if flag == 0
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
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;
|
||||||
|
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 && verbose
|
||||||
|
fprintf( 'WARNING : Failed to find a solution to the linear system\n' );
|
||||||
|
end
|
|
@ -149,7 +149,8 @@ MODFILES = \
|
||||||
simul/Solow_no_varexo.mod \
|
simul/Solow_no_varexo.mod \
|
||||||
simul/simul_ZLB_purely_forward.mod \
|
simul/simul_ZLB_purely_forward.mod \
|
||||||
simul/simul_ZLB_purely_forward_no_solution.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/2/fs2000_est.mod \
|
||||||
conditional_forecasts/3/fs2000_conditional_forecast_initval.mod \
|
conditional_forecasts/3/fs2000_conditional_forecast_initval.mod \
|
||||||
conditional_forecasts/4/fs2000_conditional_forecast_histval.mod \
|
conditional_forecasts/4/fs2000_conditional_forecast_histval.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
|
Loading…
Reference in New Issue