# EvalFnOnAgentDist\_AggVars\_Case1 and EvalFnOnAgentDist\_AllStats\_Case1 bug

**URL:** <https://discourse.vfitoolkit.com/t/evalfnonagentdist-aggvars-case1-and-evalfnonagentdist-allstats-case1-bug/438>\
**Category:** Uncategorized\
**Created:** [September 26, 2025, 3:02pm UTC](https://discourse.vfitoolkit.com/t/evalfnonagentdist-aggvars-case1-and-evalfnonagentdist-allstats-case1-bug/438 "2025-09-26T15:02:59Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![aledinola](https://discourse.vfitoolkit.com/letter_avatar_proxy/v4/letter/a/3ec8ea/32.png) [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)\
**Post date:** [September 26, 2025, 3:02pm UTC](https://discourse.vfitoolkit.com/t/evalfnonagentdist-aggvars-case1-and-evalfnonagentdist-allstats-case1-bug/438/1 "2025-09-26T15:02:59Z")

</div>

I found that `EvalFnOnAgentDist_AggVars_Case1` and `EvalFnOnAgentDist_AllStats_Case1` give different results for the aggregate variables, i.e. field “Mean”. Please see minimum working example at the end of this post that reproduces the issue.

Suppose I want to calculate average income. I define

`FnsToEvaluate.income = @(aprime,a,z,P) f_income(aprime,a,z,P);`

Then I compute policy functions and stationary distributions.  
Then I call the functions to compute averages and other moments.

```matlab
AggVars=EvalFnOnAgentDist_AggVars_Case1(StationaryDist, Policy, FnsToEvaluate, Params, FnsToEvaluateParamNames, n_d, n_a, n_z, d_grid, a_grid, z_grid, simoptions);

AllStats=EvalFnOnAgentDist_AllStats_Case1(StationaryDist,Policy,FnsToEvaluate,Params,[],n_d,n_a,n_z,d_grid,a_grid,z_grid,simoptions);

```

The moment average income should be equal to

AggVars.income.Mean = 0.13  
AllStats.income.Mean = NaN

The reason for this discrepancy is that the function for income, when evaluated on the grid, contains at least one NaN value. Now, this is not a problem in `EvalFnOnAgentDist_AggVars_Case1` because this function computes the average by excluding NaN values, as it should be:

```matlab
temp=Values.*StationaryDistVec;
AggVars(ff)=sum(temp(~isnan(temp)));

```

However, the function `EvalFnOnAgentDist_AllStats_Case1` does not implement this check and therefore returns NaN as output for the mean of the variable. Indeed, in the subfunction stats from weighted grid we have:

```matlab
% Eliminate all the zero-weights from these (trivial increase in runtime, but makes it easier to spot when there is no variance)
temp=logical(Weights==0);
Weights=Weights(~temp);
Values=Values(~temp);
% code omitted
WeightedSortedValues=SortedValues.*SortedWeights;
AllStats.Mean=sum(WeightedSortedValues);

```

I think this should be fixed by simply replacing the last line with

```matlab
valid = ~isnan(WeightedSortedValues);
AllStats.Mean=sum(WeightedSortedValues(valid));

```

or something similar.

To provide some more context on why this happens in the code I am working on:  
In the model I am working on this turned out to be very important. In some cases the income of entrepreneurs is negative because of a fixed operating cost. This gives NaN at some point. For some reason, these points do receive a positive mass in equilibrium, so they are not eliminated by `StatsFromWeightedGrid`.

There could be something wrong in terms of economics in my model of course. However, I think the functions `EvalFnOnAgentDist_AggVars_Case1` and `EvalFnOnAgentDist_AllStats_Case1` should return the same result for shared statistics (like the mean) so this is why I dared to claim there is a bug 🙂

---

<div class="post-metadata">

**Author:** ![aledinola](https://discourse.vfitoolkit.com/letter_avatar_proxy/v4/letter/a/3ec8ea/32.png) [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)\
**Post date:** [September 28, 2025, 12:36pm UTC](https://discourse.vfitoolkit.com/t/evalfnonagentdist-aggvars-case1-and-evalfnonagentdist-allstats-case1-bug/438/2 "2025-09-28T12:36:21Z")

</div>

**Minimum working example that shows the bug:**

```matlab
% Example based on Aiyagari (1994).
clear,clc,close all
addpath(genpath('C:\Users\aledi\Documents\GitHub\VFIToolkit-matlab'))
% These codes set up and solve the Aiyagari (1994) model for a given
% parametrization. After solving the model they then show how some of the
% vfitoolkit commands to easily calculate things like the Gini coefficient
% for income, and how to plot the distribution of asset holdings.
%
% VFI Toolkit automatically detects hardware (GPU? Number of CPUs?) and
% sets defaults accordingly. It will run without a GPU, but slowly. It is
% indended for use with GPU.

%% Set some basic variables

% VFI Toolkit thinks of there as being:
% k: an endogenous state variable (assets)
% z: an exogenous state variable (exogenous labor supply)

% Size of the grids
n_k=2^9;%2^9;
n_z=11; %21;

% Parameters
Params.beta=0.96; %Model period is one-sixth of a year
Params.alpha=0.36;
Params.delta=0.08;
Params.mu=3;
Params.sigma=0.2;
Params.rho=0.6;

%% Set up the exogenous shock process
% Create markov process for the exogenous labour productivity, l.
Tauchen_q=3; % Footnote 33 of Aiyagari(1993WP, pg 25) implicitly says that he uses q=3
[z_grid,pi_z]=discretizeAR1_Tauchen(0,Params.rho,sqrt((1-Params.rho^2)*Params.sigma^2),n_z,Tauchen_q);
% Note: sigma is standard deviations of s, input needs to be standard deviation of the innovations
% Because s is AR(1), the variance of the innovations is (1-rho^2)*sigma^2
% Note: Nowadays you would be better using Farmer-Toda than Tauchen, but leaving it here to match the original

[z_mean,z_variance,z_corr,~]=MarkovChainMoments(z_grid,pi_z);
z_grid=exp(z_grid);
% Get some info on the markov process
[Expectation_l,~,~,~]=MarkovChainMoments(z_grid,pi_z); %Since l is exogenous, this will be it's eqm value 
% Note: Aiyagari (1994) actually then normalizes l by dividing it by Expectation_l (so that the resulting process has expectation equal to 1)
z_grid=z_grid./Expectation_l;
[Expectation_l,~,~,~]=MarkovChainMoments(z_grid,pi_z);
% If you look at Expectation_l you will see it is now equal to 1
Params.Expectation_l=Expectation_l;

%% Grids

% In the absence of idiosyncratic risk, the steady state equilibrium is given by
r_ss=1/Params.beta-1;
K_ss=((r_ss+Params.delta)/Params.alpha)^(1/(Params.alpha-1)); %The steady state capital in the absence of aggregate uncertainty.

% Set grid for asset holdings
k_grid=10*K_ss*(linspace(0,1,n_k).^3)'; % linspace ^3 puts more points near zero, where the curvature of value and policy functions is higher and where model spends more time

% Bring model into the notational conventions used by the toolkit
d_grid=0; %There is no d variable
a_grid=k_grid;
% z_grid
% pi_z;

n_d=0;
n_a=n_k;
% n_z

% Create functions to be evaluated
FnsToEvaluate.K = @(aprime,a,z) a; %We just want the aggregate assets (which is this periods state)

% Now define the functions for the General Equilibrium conditions
    % Should be written as LHS of general eqm eqn minus RHS, so that the closer the value given by the function is to 
    % zero, the closer the general eqm condition is to holding.
GeneralEqmEqns.CapitalMarket = @(r,K,alpha,delta,Expectation_l) r-(alpha*(K^(alpha-1))*(Expectation_l^(1-alpha))-delta); %The requirement that the interest rate corresponds to the agg capital level
% Inputs can be any parameter, price, or aggregate of the FnsToEvaluate

fprintf('Grid sizes are: %i points for assets, and %i points for exogenous shock \n', n_a,n_z)

%%
DiscountFactorParamNames={'beta'};

ReturnFn=@(aprime, a, z, alpha,delta,mu,r) Aiyagari1994_ReturnFn(aprime, a, z,alpha,delta,mu,r);
% The first inputs must be: next period endogenous state, endogenous state, exogenous state. Followed by any parameters

%%

% Use the toolkit to find the equilibrium price index
GEPriceParamNames={'r'};
% Set initial value for interest rates (Aiyagari proves that with idiosyncratic
% uncertainty, the eqm interest rate is limited above by it's steady state value
% without idiosyncratic uncertainty, that is that r<r_ss).
Params.r=0.038;

%%
% Solve for the stationary general equilbirium
vfoptions=struct(); % Use default options for solving the value function (and policy fn)
vfoptions.howardsgreedy=0;
simoptions=struct(); % Use default options for solving for stationary distribution
heteroagentoptions.verbose=1; % verbose means that you want it to give you feedback on what is going on

%fprintf('Calculating price vector corresponding to the stationary general eqm \n')
%[p_eqm,~,GeneralEqmCondn]=HeteroAgentStationaryEqm_Case1(n_d, n_a, n_z, 0, pi_z, d_grid, a_grid, z_grid, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Params, DiscountFactorParamNames, [], [], [], GEPriceParamNames,heteroagentoptions, simoptions, vfoptions);

%p_eqm % The equilibrium values of the GE prices
p_eqm.r = Params.r;

%% Now that we have the GE, let's calculate a bunch of related objects
%Params.r=p_eqm.r; % Put the equilibrium interest rate into Params so we can use it to calculate things based on equilibrium parameters

% Equilibrium wage
Params.w=(1-Params.alpha)*((Params.r+Params.delta)/Params.alpha)^(Params.alpha/(Params.alpha-1));

fprintf('Calculating various equilibrium objects \n')
[V,Policy]=ValueFnIter_Case1(n_d,n_a,n_z,d_grid,a_grid,z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], vfoptions);
% V is value function
% Policy is policy function (but as an index of k_grid, not the actual values)

% PolicyValues=PolicyInd2Val_Case1(Policy,n_d,n_a,n_s,d_grid,a_grid,vfoptions); % This will give you the policy in terms of values rather than index
StationaryDist=StationaryDist_Case1(Policy,n_d,n_a,n_z,pi_z, simoptions);

FnsToEvaluate.income = @(aprime,a,z,r,w) r*a+w*z;

FnsToEvaluate.income = @(aprime,a,z,r,w) f_income(aprime,a,z,r,w);

ValuesOnGrid=EvalFnOnAgentDist_ValuesOnGrid_Case1(Policy,FnsToEvaluate,Params,[],n_d,n_a,n_z,d_grid,a_grid,z_grid,simoptions,[],StationaryDist);

incomeOnGrid = gather(ValuesOnGrid.income);
disp('Are there NaN values in IncomeOnGrid? 1=yes,0=no')
any(isnan(incomeOnGrid),"all")
sum(isnan(incomeOnGrid(:)))/numel(incomeOnGrid)

simoptions.whichstats = zeros(7,1);
simoptions.whichstats(1) = 1;

AggVars=EvalFnOnAgentDist_AggVars_Case1(StationaryDist, Policy, FnsToEvaluate,Params, [],n_d, n_a, n_z, d_grid, a_grid,z_grid,simoptions);
AllStats=EvalFnOnAgentDist_AllStats_Case1(StationaryDist, Policy, FnsToEvaluate,Params, [],n_d, n_a, n_z, d_grid, a_grid,z_grid,simoptions);

% AggVars contains the aggregate values of the 'FnsToEvaluate' (in this model aggregates are equal to the mean expectation over the agent distribution)
% Currently the only FnsToEvaluate is assets, so we get aggregate capital stock
AggVars.income.Mean
AllStats.income.Mean

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function income = f_income(aprime,a,z,r,w)

income_temp = r*a+w*z;

if z>1.6
    income = NaN;
else
    income = income_temp;
end

end

function F=Aiyagari1994_ReturnFn(aprime, a, z,alpha,delta,mu,r)
% The return function is essentially the combination of the utility
% function and the constraints.

F=-Inf;
w=(1-alpha)*((r+delta)/alpha)^(alpha/(alpha-1));
c=w*z+(1+r)*a-aprime; % Budget Constraint
% c=w l_t+(1+r)a_t-a_{t+1}
if c>0
    if mu==1
        F=log(c);
    else
        F=(c^(1-mu) -1)/(1-mu);
    end
end

end

```

---

<div class="post-metadata">

**Author:** ![aledinola](https://discourse.vfitoolkit.com/letter_avatar_proxy/v4/letter/a/3ec8ea/32.png) [@aledinola](https://discourse.vfitoolkit.com/u/aledinola)\
**Post date:** [September 28, 2025, 12:38pm UTC](https://discourse.vfitoolkit.com/t/evalfnonagentdist-aggvars-case1-and-evalfnonagentdist-allstats-case1-bug/438/3 "2025-09-28T12:38:48Z")

</div>

Unrelated to the bug: if I run the VFI in the Aiyagari code with `vfoptions.howardsgreedy=0;`, it runs 5 times faster.  
Same number of grid points and same parameters as in the vfi toolkit examples. May be worth it to revisit the default options 🙂

---

<div class="post-metadata">

**Author:** ![robertdkirkby](https://discourse.vfitoolkit.com/user_avatar/discourse.vfitoolkit.com/robertdkirkby/32/287_2.png) [@robertdkirkby](https://discourse.vfitoolkit.com/u/robertdkirkby)\
**Post date:** [October 29, 2026, 11:33pm UTC](https://discourse.vfitoolkit.com/t/evalfnonagentdist-aggvars-case1-and-evalfnonagentdist-allstats-case1-bug/438/4 "2026-10-29T23:33:58Z")

</div>

Following up, this all got fixed a while back but I didn’t put a reply here.

Both now report NaN, which is the intented/correct answer.

And the howardsgreedy default settings are now based on a large test bank that aims to look at what is typically the fastest across various situations.