I have checked also some inequality statistics like the Gini coefficient.
The Gini of assets conditional on health=sick is
g_toolkit =
0.4956
and
g_mycode =
0.4956
So that’s good (small discrepancies with further digits, but who cares). But if I do the gini of assets conditional on health=sick and on age, I get some differences b/w the toolkit and my own code, please see plot below:
Here comes the bad news 
I went back to my project and if I run the code without inequality statistics, with simoptions.whichstats = [1,1,1,0,0,0,0];, all is good.
But if I set simoptions.whichstats = [1,1,1,2,0,0,0];, then I get the following error:
Error using gpuArray/interp1
Interpolation requires at least two sample points in each dimension.
Error in StatsFromWeightedGrid (line 158)
temp=interp1(CumSumSortedWeights2,CumSumSortedWeightedValues(u1index),llvec(1:end-1));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in EvalFnOnAgentDist_AllStats_FHorz_Case1_PType (line 407)
AllStats.(CondlRestnFnNames{rr}).(FnsToEvalNames{ff}).(Names_i{ii})=StatsFromWeightedGrid(SortedValues,RestrictedSortedWeights,simoptions.npoints,simoptions.nquantiles,simoptions.tolerance,1,simoptions.whichstats); % 1 is presorted
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in main (line 303)
AllStats=EvalFnOnAgentDist_AllStats_FHorz_Case1_PType(StatDist,Policy,FnsToEval,Params,n_d,n_a,n_z,N_j,N_i,d_grid,a_grid,z_grid_J,simoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Surprisingly, if I do not set simoptions.whichstats at all (leaving the deafult settings on), I do not get any error
Btw, I checked the Gini and it seems OK (modulo a small difference, probably due to the fact that I use my own function mylorenz which uses a slightly different method).
I would be happy to use the defaults for simoptions.whichstats, but computational time is a bit of an issue for me and so I’d like to be able to set the last three of simoptions.whichstats to zero (I need Lorenz and Gini but not the other statistics).
@robertdkirkby Please let me know if you would like me to share the code for the project I am referring to.
The defaults are now:
simoptions.whichstats=[1,1,1,1,1,1,1] for AllStats
simoptions.whichstats=[1,1,1,2,1,2,1] for LifeCycleProfiles
(4th entry is Lorenz+Gini, 6th entry is Quantiles)
For AllStats I figure that it is better to worry more about memory and less about speed. For LifeCycleProfiles the memory is less likely to be an issue (as just doing conditional on age) and speed will be more important (as do lot’s of stats, doing them for each age).
Please share code that generates these errors, I don’t get that the thing fails on any of my examples (dropbox link is good if you don’t want to be putting beta codes on github).
Thanks for your reply, Robert! I sent to your email a Dropbox link.
Please run main.m with line 36 commented out:
%simoptions.whichstats = [1,1,1,2,0,0,0];
All is good and the screen output on my machine is:
Check pi_zh_J.low...
Check pi_zh_J.high...
Check pi_z_J.low...
Check pi_z_J.high...
Test ValueFnIter
Check size of value function
n_a, n_z(1), n_z(2), n_e, Age
201 7 2 5 80
201 7 2 5 80
Check size of policy function
2, n_a, n_z(1), n_z(2), n_e, Age
2 201 7 2 5 80
2 201 7 2 5 80
errV =
5.6843e-13
errP =
0
err_D =
8.6736e-19
ans =
'Possible Warning: Number of columns of n_z is the same as the number of permanent types.
This may just be coincidence as number of z variables is equal to number of permanent types.
If they are intended to be permanent types then n_z should have them as different rows (not columns).
'
> In EvalFnOnAgentDist_ValuesOnGrid_FHorz_Case1_PType (line 164)
In fun_model_moments (line 45)
In main (line 307)
ans =
'Possible Warning: Number of columns of n_z is the same as the number of permanent types.
This may just be coincidence as number of z variables is equal to number of permanent types.
If they are intended to be permanent types then n_z should have them as different rows (not columns).
'
> In EvalFnOnAgentDist_ValuesOnGrid_FHorz_Case1_PType (line 164)
In fun_model_moments (line 45)
In main (line 307)
ans =
logical
0
ans =
logical
0
ans =
logical
0
ans =
logical
0
RUNNING TIMES
VFI toolkit = 1.188982
VFI cpu = 3.514840
AgeStats = 17.585771
AllStats = 1.328889
Moments (my code) = 0.338687
If you uncomment line 36, on my machine I get the following error (after I run VFI and distribution…):
Check pi_zh_J.low...
Check pi_zh_J.high...
Check pi_z_J.low...
Check pi_z_J.high...
Test ValueFnIter
Check size of value function
n_a, n_z(1), n_z(2), n_e, Age
201 7 2 5 80
201 7 2 5 80
Check size of policy function
2, n_a, n_z(1), n_z(2), n_e, Age
2 201 7 2 5 80
2 201 7 2 5 80
errV =
5.6843e-13
errP =
0
err_D =
8.6736e-19
Error using gpuArray/interp1
Interpolation requires at least two sample points in each dimension.
Error in StatsFromWeightedGrid (line 158)
temp=interp1(CumSumSortedWeights2,CumSumSortedWeightedValues(u1index),llvec(1:end-1));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in EvalFnOnAgentDist_AllStats_FHorz_Case1_PType (line 407)
AllStats.(CondlRestnFnNames{rr}).(FnsToEvalNames{ff}).(Names_i{ii})=StatsFromWeightedGrid(SortedValues,RestrictedSortedWeights,simoptions.npoints,simoptions.nquantiles,simoptions.tolerance,1,simoptions.whichstats); % 1 is presorted
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in main (line 303)
AllStats=EvalFnOnAgentDist_AllStats_FHorz_Case1_PType(StatDist,Policy,FnsToEval,Params,n_d,n_a,n_z,N_j,N_i,d_grid,a_grid,z_grid_J,simoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Yeah, it was a rounding error at numerical accuracy thing.
It is possible that, e.g., Weights has two elements, but one of them is 1e-20~ish, and then you do CumWeights=cumsum(Weights) and just get that it equals [1,1]. So when I do unique(CumWeights) I am left with one element.
The main formula for LorenzCurve with whichstats(4)==6 requires at least two elements, and because it was left with just one element things were going wrong. I just put in a check so if there is only one point left then it just gives LorenzCurve=1/npoints:1/npoints:1.
Anyway, fixed now
For graph of “Gini coefficient of assets, by age, given sick”: Do both codes necessarily give same answer (does one use more points/interpolation? and if yes, then if you turn up grids on both do they get closer?)
Or are they both being calculated from toolkit solutions for Policy and StatDist?
Actually, running your main.m, I think now the graph shows both with exact same line?? So this is now also working right???
Thanks, Robert!
For the Gini, I use the distribution and the policy computed by the toolkit, but then I compute the Lorenz curve and the coefficient using my own function mylorenz.m. This function uses a slightly different method (no interpolation) so the results might be slightly different.
Another check I can do is to plot the coefficient of variation (standard deviation divided by the mean): this should give exactly the same answer.
EDIT
After running main.m, I also verified that the lines for the Gini are the same (except for the first point, which is however uninmportant since all agents start with the same level of assets, so the Gini coefficient is zero or undefined anyways). So now all is working well!