V2.2: Same same but faster

robertdkirkby · November 20, 2023, 12:37am

v2.2 of VFI Toolkit has been released.

Not really any new features, but everything is faster now. Mostly due to rolling out the Tan improvement to all agent distribution computations, plus because the agent distribution was then so fast the time taken for FnsToEvaluate got embarassing so I went and overhauled this for speed too.

aledinola · November 20, 2023, 12:07pm

Great news, thanks for this update Robert! The function SimPanelValues_FHorz_Case1 still doesn’t work for the case of semi-exogenous shocks, right? I ran our code in v3 and I got the same error message as before, reproduced below:

Error using gpuArray/reshape
Number of elements must not change. Use as one of the size inputs to automatically calculate the appropriate size for that dimension.

Error in KronPolicyIndexes_FHorz_Case1 (line 22)
Policy=reshape(Policy,[size(Policy,1),N_a,N_z,N_j]);

Error in SimPanelValues_FHorz_Case1 (line 219)
PolicyIndexesKron=KronPolicyIndexes_FHorz_Case1(Policy, n_d, n_a, n_z, N_j); % Create it here as want it both here and inside SimPanelIndexes_FHorz_Case1 (which will recognise that it is already in this form)

Error in main_v3 (line 297)
SimPanelValues=SimPanelValues_FHorz_Case1(InitialDist,Policy,FnsToEvaluate,,Params,n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,pi_z,simoptions);

aledinola · November 20, 2023, 12:13pm

I checked also the simple example in examples\Finite horizon and I ran the code here:

github.com

vfitoolkit/VFItoolkit-matlab-examples/blob/master/FiniteHorizon/FiniteHorzStochConsSavings_Extended.m

% This code implements a simple ten-period stochastic consumption savings
% problem to demonstrate the finite value function command of the VFI
% Toolkit.
% This 'extended' version then goes on to simulate panel data and
% life-cycle profiles for earnings and assets. ('extension' begins on line 84)

% Main strength of the finite value function command is that you simply
% declare parameters that depend on age as row vectors, and the VFI Toolkit
% automatically realizes that these are age dependent.


% Agents live ten periods. Income consists of a deterministic function of
% age Wj plus an AR(1) shock Wz.

% Age
% One endogenous variable: assets
% One stochastic exogenous variables: income shock

J=10; % Ages 1 to 10 inclusive.

This file has been truncated. show original

However I get an error at line 102 when running

SimPanelValues=SimPanelValues_FHorz_Case1(InitialDist,Policy, FnsToEvaluate,,Params,0,n_a,n_z,N_j,0,a_grid,z_grid,pi_z, simoptions);

The error message is the following. Note that before the error the code gives a warning about a potential problem with parfor

Warning: File: SimLifeCycleProfiles_FHorz_Case1.m Line: 223 Column: 17
The temporary variable ‘d_val’ will be cleared at the beginning of each iteration of the parfor-loop. If ‘d_val’ is used before it is set, a runtime
error will occur. For more information, see Parallel for Loops in MATLAB, “Uninitialized Temporaries”.

In FiniteHorzStochConsSavings_Extended (line 124)
Error using SimLifeCycleIndexes_FHorz_Case1_raw
Index in position 3 exceeds array bounds. Index must not exceed 1.

Error in SimPanelIndexes_FHorz_Case1 (line 156)
parfor ii=1:simoptions.numbersims % This is only change from the simoptions.parallel==0

Error in SimLifeCycleProfiles_FHorz_Case1 (line 101)
SimPanelIndexes=SimPanelIndexes_FHorz_Case1(InitialDist,PolicyIndexesKron,n_d,n_a,n_z,N_j,pi_z, simoptions);

Error in FiniteHorzStochConsSavings_Extended (line 124)
LifeCycleProfiles=SimLifeCycleProfiles_FHorz_Case1(InitialDist,Policy, FnsToEvaluate,Params,,0,n_a,n_z,N_j,0,a_grid,z_grid,pi_z,simoptions);

robertdkirkby · November 21, 2023, 4:12am

Panel data simulation should work now (totally overhauled it so it is also faster)

aledinola · November 21, 2023, 2:04pm

Hi, thank you so much as always!
I perfomed a check on the panel simulation commands and I found some potential problem. I computed the life-cycle profile of a variable of interest, say, earnings, using the distribution as

x = AgeConditionalStats.earnings.Mean

Then I simulate a T*N panel with N=50000 and I average out across the individuals,

y = mean(simPanelValues.earnings,2)

For N large enough, x and y should be very close. Unfortunately, this does not seem to be the case for earnings (assets are ok-ish). I show here the results of my check:

Note that the discrepancy between the “simulation” and the “distribution” tends to disappear for older ages. This suggests that there might be a problem with the initial condition for the simulation. I don’t understand why though, because I pass the following initial distribution to both functions, so with N large there should not be any significant difference:

jequaloneDist = zeros([n_a,n_semiz],‘gpuArray’);

% All agents start with zero assets, 80% are employed (both male and
% female), there is an expansion (good state). lm and lf are statistically
% independent, so that the prob that both m and f are employed is 0.80.8,
% etc.
% jequaloneDist(a,lm,lf,za)
jequaloneDist(1,:,:,1) = [0.80.8,0.80.2;
0.20.8,0.2*0.2];

Please let me know if I am overlooking something

NOTES

(1) I have programmed my own simulation on the cpu and it confirms the life-cycles obtained using the distribution (the cpu code is commented out). This means that the potential problem is in the panel simulation codes.

(2) I have tried also N larger but the results are pretty much the same

(3) Everything in this post refers to main_v3.

aledinola · November 21, 2023, 2:59pm

I did a bit of digging in the codes…
The relevant function is SimPanelIndexes_FHorz_Case1_noz_semiz because we are simulating the panel indexes for semi-exogenous shock and no z variable. In line 64 we have

parfor ii=1:simoptions.numbersims
seedpointvec=max(cumsumInitialDistVec>rand(1));
seedpoints(ii,:)=[ind2sub_homemade([N_a,N_semiz],seedpointvec),1];
end

Here we have an initial distribution (i.e. the initial condition for the simulation) and we must draw the initial sample from this distribution. We do this for all numbersims individuals that we want to simulate. I would replace the line

seedpointvec=max(cumsumInitialDistVec>rand(1));

with

[~,seedpointvec] = max(cumsumInitialDistVec>rand(1));

because you want to find the position of rand(1) in the vector cumsumInitialDistVec so we need the argmax, not the max. Note: you can do the same using the find function in Matlab.

Note also that in other places in the code we have the argmax, so I guess it’s just a typo

aledinola · November 21, 2023, 3:11pm

Last thing: the parfor for drawing the initial conditions slows down a bit my code. Since we are using the parallel option, I would replace this:

parfor ii=1:simoptions.numbersims % TODO faster without parfor
%seedpointvec=max(cumsumInitialDistVec>rand(1)); TODO
[~,seedpointvec]=max(cumsumInitialDistVec>rand(1));
seedpoints(ii,:)=[ind2sub_homemade([N_a,N_semiz],seedpointvec),1];
end

with the following

if simoptions.parallel==0

for ii=1:simoptions.numbersims % TODO faster without parfor
%seedpointvec=max(cumsumInitialDistVec>rand(1)); TODO
[~,seedpointvec]=max(cumsumInitialDistVec>rand(1));
seedpoints(ii,:)=[ind2sub_homemade([N_a,N_semiz],seedpointvec),1];
end

else

parfor ii=1:simoptions.numbersims % TODO faster without parfor
%seedpointvec=max(cumsumInitialDistVec>rand(1)); TODO
[~,seedpointvec]=max(cumsumInitialDistVec>rand(1));
seedpoints(ii,:)=[ind2sub_homemade([N_a,N_semiz],seedpointvec),1];
end
end

robertdkirkby · November 21, 2023, 11:21pm

Fixed and vectorized (no loops at all so is also marginally faster)

You were absolutely right that problem was in the initial distribution and that it should have been argmax (which is why it was wrong, the baseline case was correct but for some reason not semiz).

Trick to vectorizing is that if cumsum_aaa is (cumulative sum of) a column vector and bbb is a row vector of rand() then [~,index]=max(cumsum_aaa>bbb) gives the indexes for each of the bbb being drawn from cumsum_aaa .
[That and I have an ind2sub_vec_homemade() for when the input is a column vector precisely to allow vectorizing this things]

Actually really I should vectorize the whole of the panel simulation, but I just don’t have the heart for it right now (spent last 2-3 months overhauling the insides of the toolkit, I need a break from overhauling, do something new )

aledinola · November 22, 2023, 11:43am

Great, thanks a lot! Don’t worry for the panel simulation, it is already faster now (for some reasons the parallel cpu was slowing things down)