Bugs in VFItoolkit-matlab-examples

1

I will collect in this post bugs related to the codes in VFItoolkit-matlab-examples.

  • StochasticNeoClassicalGrowthModel if Pijoan=1 (if I run Pijoan=0 and then Pijoan=1 WITHOUT clearing the workspace, no error.)

StochasticNeoClassicalGrowthModel
Unrecognized field name “sigma_epsilon”.

Error in StochasticNeoClassicalGrowthModel (line 50)
[z_grid,pi_z]=discretizeAR1_Tauchen(0,Params.rho,Params.sigma_epsilon,n_z,Tauchen_q);
^^^^^^^^^^^^^^^^^^^^

1 Like
2

I propose an updated version of BasicRealBusinessCycleModel.m
Nothing of substance, but something useful to interpret the model results:

  1. Convert policy from indices to values. Get policy for labor d(k,z) and policy for next-period capital k’(k,z).
  2. Add vfoptions=struct(); It is not really needed but is clearer that we are using the default settings.
    2b) Add vfoptions for grid interpolation BUT commented out, just to tell users that there is this option, but is more advanced.
  3. Add plots for policy functions, besides the plot for value function.
% Example using a variant of the Basic RBC model (following Aruoba, Fernandez-Villaverde, & Rubio-Ramirez, 2006)
%
% Note: On GPU this will take fractions of a second, on CPUs you will want to make the grids smaller (n_d,n_a,n_z) or it will take a while.

clear
tic;
UseAlternativeParams=1;

%% Set up

% Number of grid points
n_d=50;
n_a=250;
n_z=21;  %Note, for the figures to correctly use z=0 this must be an odd number (make it 39)

%Discounting rate
Params.beta = 0.9896;

%Parameter values
Params.alpha = 0.4; % alpha
Params.theta=0.357;
Params.rho = 0.95; % rho
Params.tau=2;
Params.delta = 0.0196; % delta
Params.sigma_epsilon=0.007;

if UseAlternativeParams==1
    Params.tau=50;
    Params.sigma_epsilon=0.035;
end

%% Compute the steady state (just use these when picking grids)
varphi=(1/Params.alpha*(1/Params.beta-1+Params.delta))^(1/(1-Params.alpha));
Psi=Params.theta/(1-Params.theta)*(1-Params.alpha)*varphi^(-Params.alpha);
Omega=varphi^(1-Params.alpha)-Params.delta;
K_ss=Psi/(Omega+varphi*Psi);


%% Create grids (it is very important that each grid is defined as a column vector)
Tauchen_q=3; %Parameter for the Tauchen method
[z_grid,pi_z]=discretizeAR1_Tauchen(0,Params.rho,Params.sigma_epsilon,n_z,Tauchen_q);

a_grid=linspace(0.01,5*K_ss,n_a)';
d_grid=linspace(0,1,n_d)';

%% Now, create the return function 
DiscountFactorParamNames={'beta'};

ReturnFn=@(d_val,aprime_val, a_val, s_val,alpha,delta,theta,tau) BasicRealBusinessCycleModel_ReturnFn(d_val,aprime_val, a_val, s_val,alpha,delta,theta,tau);

%% Solve
%Do the value function iteration. Returns both the value function itself, and the optimal policy function.
vfoptions = struct(); % Use default options

% If you want to solve the VFI using linear interpolation (GPU only), uncomment the following:
% vfoptions.gridinterplayer  = 1;
% vfoptions.ngridinterp      = 15;   % number of extra points between grid points
% simoptions.gridinterplayer = vfoptions.gridinterplayer;
% simoptions.ngridinterp     = vfoptions.ngridinterp;

[V,Policy]=ValueFnIter_Case1(n_d,n_a,n_z,d_grid,a_grid,z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames,[]);
time1=toc;

% Convert policy from indices to values
PolicyValues = PolicyInd2Val_Case1(Policy, n_d, n_a, n_z, d_grid, a_grid, vfoptions);
PolicyValues_d = squeeze(PolicyValues(1,:,:));
PolicyValues_k = squeeze(PolicyValues(2,:,:));

%% Report some output
fprintf('Run time for value function iteration: %8.2f seconds \n', time1)

%% Graph of the value function.
figure(1)
surf(a_grid*ones(1,n_z), ones(n_a,1)*z_grid', V)

%% Graph of the policy function

figure(2)
plot(a_grid,a_grid,'--')
hold on
plot(a_grid,PolicyValues_k(:,1))
hold on
plot(a_grid,PolicyValues_k(:,n_z))
legend('45 line','Given z(1)','Given z(n_z)')
title('Policy function for next period capital')
xlabel('Current period capital')

figure(3)
plot(a_grid,PolicyValues_d(:,1))
hold on
plot(a_grid,PolicyValues_d(:,n_z))
legend('Given z(1)','Given z(n_z)')
title('Policy function for labor')
xlabel('Current period capital')
1 Like
3

Another bug in BasicRealBusinessCycleModel_BusinessCycleStatistics

Running BasicRealBusinessCycleModel_BusinessCycleStatistics.m
No vfoptions given, using defaults
Run time for value function iteration:     1.53 seconds 
Unrecognized function or variable 'TimeSeries_Case1'.

Error in BasicRealBusinessCycleModel_BusinessCycleStatistics (line 18)
TimeSeries=TimeSeries_Case1(Policy, FnsToEvaluate, Params, n_d, n_a, n_z, d_grid, a_grid, z_grid,pi_z,simoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
4

Bug in BasicRealBusinessCycleModel_NumericalErrors

`Running BasicRealBusinessCycleModel_NumericalErrors.m
WARNING: Marcet-DenHaan statistics are not correctly implemented
No vfoptions given, using defaults
Run time for value function iteration:     1.19 seconds 
Unrecognized function or variable 'StationaryDist_Case1_Simulation'.

Error in BasicRealBusinessCycleModel_NumericalErrors (line 47)
    StationaryDist(:,:,ii)=StationaryDist_Case1_Simulation(Policy,n_d,n_a,n_z,pi_z,simoptions);
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 `
5

@robertdkirkby I’d be interested in simulation commands for infinite horizon models (for either replication or to compute statistics that are not easy to get from the stationary distribution). Could you please point to some example of toolkit command that does simulation for infinite horizon models?
Thanks a lot!!

P.S. Maybe instead of updating these older examples (based on representative agent models) it could be a good idea to add a simulation command to the Aiyagari example. I like that example since it shows many things the toolkit can do in the context of the textbook heterogeneous agent, infinite horizon model.

2 Likes
6

Fixed. Except I didn’t really bother with BasicRealBusinessCycleModel_NumericalErrors

There are two commands:
SimPanelValues_Case1
and
SimTimeSeriesValues_Case1
both are based around using FnToEvaluate and produce the simulations for these.

[Note, the later is actually internally just calling the former with simoptions.numbersims=1 anyway]

1 Like
7

Great, thanks!
I will use the commands for sim panel values then

8

I checked BasicRealBusinessCycleModel and you set

vfoptions.gridinterplayer=1; % for next period capital, interpolate between grid points...
vfoptions.ngridinterp=20; % ...with 20 more evenly spaced grid points

in lines 51-52. This runs fine on the GPU but quite surprisingly also on the CPU, but it shouldn’t since gridinterplayer is not implemented on CPU. More precisely, it runs the VFI using pure discretization and then it breaks later on with this error:

Error using gpuArray
Unable to find a supported GPU device.

Error in PolicyInd2Val_Case1 (line 80)
                aprime_grid=interp1(gpuArray(1:1:N_a)',a_grid,linspace(1,N_a,N_a+(N_a-1)*vfoptions.ngridinterp)');
                                    ^^^^^^^^^^^^^^^^^^
Error in BasicRealBusinessCycleModel (line 68)
PolicyVals=PolicyInd2Val_Case1(Policy,n_d,n_a,n_z,d_grid,a_grid,vfoptions);
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Proposed fix: The CPU codes should return an informative error if the user sets vfoptions.gridinterplayer=1 that says “gridinterplayer not implemented for cpu”. Lines 58-61 of ValueFnIter_InfHorz_CPU seem a good place to put this check, since there is already a check for experienceasset.

%% Anything but the basics is only for GPU
if vfoptions.experienceasset==1
    error('Cannot use vfoptions.experienceasset=1 without a GPU')
end
9

Good point.

They now check for gpu, and only do gridinterplayer if you have one

1 Like