Bug in Aiyagari1994 example on CPU

aledinola · January 28, 2026, 3:47pm

Just tested the Aiyagari example on a laptop without a supported GPU. I got the error message shown below. I like testing this example once in a while

Aiyagari1994

Grid sizes are: 512 points for assets, and 11 points for exogenous shock
Calculating price vector corresponding to the stationary general eqm

Error using [zeros](matlab:matlab.lang.internal.introspective.errorDocCallback('zeros'))
Unable to find a supported GPU device.

Error in [UnKronPolicyIndexes_Case1](matlab:matlab.lang.internal.introspective.errorDocCallback('UnKronPolicyIndexes_Case1', 'C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\SubCodes\UnKronPolicyIndexes\UnKronPolicyIndexes_Case1.m', 16)) ([line 16](matlab: opentoline('C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\SubCodes\UnKronPolicyIndexes\UnKronPolicyIndexes_Case1.m',16,0)))
Policy=zeros(l_aprime+extra,N_a,N_z,'gpuArray');
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [ValueFnIter_Case1](matlab:matlab.lang.internal.introspective.errorDocCallback('ValueFnIter_Case1', 'C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\ValueFnIter\InfHorz\ValueFnIter_Case1.m', 649)) ([line 649](matlab: opentoline('C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\ValueFnIter\InfHorz\ValueFnIter_Case1.m',649,0)))
Policy=UnKronPolicyIndexes_Case1(Policy, n_d, n_a, n_z,vfoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [HeteroAgentStationaryEqm_Case1_subfn](matlab:matlab.lang.internal.introspective.errorDocCallback('HeteroAgentStationaryEqm_Case1_subfn', 'C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\HeterogeneousAgent\InfHorz\HeteroAgentStationaryEqm_Case1_subfn.m', 29)) ([line 29](matlab: opentoline('C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\HeterogeneousAgent\InfHorz\HeteroAgentStationaryEqm_Case1_subfn.m',29,0)))
[V,Policy]=ValueFnIter_Case1(n_d,n_a,n_z,d_grid,a_grid,z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames,ReturnFnParamNames,vfoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [HeteroAgentStationaryEqm_Case1>@(p)HeteroAgentStationaryEqm_Case1_subfn(p,n_d,n_a,n_z,d_grid,a_grid,z_gridvals,pi_z,ReturnFn,FnsToEvaluate,FnsToEvaluateCell,GeneralEqmEqnsCell,Parameters,DiscountFactorParamNames,ReturnFnParamNames,FnsToEvaluateParamNames,GeneralEqmEqnParamNames,GEPriceParamNames,GEeqnNames,AggVarNames,nGEprices,heteroagentoptions,simoptions,vfoptions)](matlab:matlab.lang.internal.introspective.errorDocCallback('HeteroAgentStationaryEqm_Case1>@(p)HeteroAgentStationaryEqm_Case1_subfn(p,n_d,n_a,n_z,d_grid,a_grid,z_gridvals,pi_z,ReturnFn,FnsToEvaluate,FnsToEvaluateCell,GeneralEqmEqnsCell,Parameters,DiscountFactorParamNames,ReturnFnParamNames,FnsToEvaluateParamNames,GeneralEqmEqnParamNames,GEPriceParamNames,GEeqnNames,AggVarNames,nGEprices,heteroagentoptions,simoptions,vfoptions)', 'C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\HeterogeneousAgent\InfHorz\HeteroAgentStationaryEqm_Case1.m', 381)) ([line 381](matlab: opentoline('C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\HeterogeneousAgent\InfHorz\HeteroAgentStationaryEqm_Case1.m',381,0)))
GeneralEqmConditionsFnOpt=@(p) HeteroAgentStationaryEqm_Case1_subfn(p, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, FnsToEvaluate, FnsToEvaluateCell, GeneralEqmEqnsCell, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames, GEeqnNames, AggVarNames, nGEprices, heteroagentoptions, simoptions, vfoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [fminsearch](matlab:matlab.lang.internal.introspective.errorDocCallback('fminsearch', 'C:\Program Files\MATLAB\R2025b\toolbox\matlab\optimfun\optimfun\fminsearch.m', 209)) ([line 209](matlab: opentoline('C:\Program Files\MATLAB\R2025b\toolbox\matlab\optimfun\optimfun\fminsearch.m',209,0)))
f = funfcn(x,varargin{:});
^^^^^^^^^^^^^^^^^^^^^
Error in [HeteroAgentStationaryEqm_Case1](matlab:matlab.lang.internal.introspective.errorDocCallback('HeteroAgentStationaryEqm_Case1', 'C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\HeterogeneousAgent\InfHorz\HeteroAgentStationaryEqm_Case1.m', 411)) ([line 411](matlab: opentoline('C:\Users\dinolaa\Documents\GitHub\VFIToolkit-matlab\HeterogeneousAgent\InfHorz\HeteroAgentStationaryEqm_Case1.m',411,0)))
[p_eqm_vec,GeneralEqmConditions]=fminsearch(GeneralEqmConditionsFnOpt,GEparamsvec0,minoptions);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in [Aiyagari1994](matlab:matlab.lang.internal.introspective.errorDocCallback('Aiyagari1994', 'C:\Users\dinolaa\Dropbox\1 - RESEARCH PROJECTS\VFI_TOOLKIT\Aiyagari_example\Aiyagari1994.m', 101)) ([line 101](matlab: opentoline('C:\Users\dinolaa\Dropbox\1 - RESEARCH PROJECTS\VFI_TOOLKIT\Aiyagari_example\Aiyagari1994.m',101,0)))
[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);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[Related documentation](matlab:helpview('parallel-computing','error_gpu_device_NoCUDADevice'))

robertdkirkby · January 29, 2026, 2:06am

Should be fixed. I’ve separated out the CPU version to be branched out early on, so hopefully it is permanently fixed (I’d done this before for FHorz, but now also for InfHorz).

yechen · February 3, 2026, 9:34am

Hi, Rob

I ran it on the CPU and received the following messages.

>> Aiyagari1994
Grid sizes are: 512 points for assets, and 11 points for exogenous shock 
Calculating price vector corresponding to the stationary general eqm 
Starting parallel pool (parpool) using the 'Processes' profile ...
2026-02-03 17:22:12: Connected to 7 of 112 parallel pool workers.
Connected to parallel pool with 112 workers.
Unrecognized function or variable 'n_aprime'.

Error in UnKronPolicyIndexes_InfHorz_CPU (line 18)
            optA=ind2sub_homemade(n_aprime',optaindexKron);
                                  ^^^^^^^^^
Error in ValueFnIter_InfHorz_CPU (line 168)
    Policy=UnKronPolicyIndexes_InfHorz_CPU(Policy, n_d, n_a, n_z);
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in ValueFnIter_Case1 (line 168)
    [V,Policy]=ValueFnIter_InfHorz_CPU(n_d,n_a,n_z,d_grid,a_grid,z_grid, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in HeteroAgentStationaryEqm_Case1_subfn (line 29)
[V,Policy]=ValueFnIter_Case1(n_d,n_a,n_z,d_grid,a_grid,z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames,ReturnFnParamNames,vfoptions);
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in HeteroAgentStationaryEqm_Case1>@(p)HeteroAgentStationaryEqm_Case1_subfn(p,n_d,n_a,n_z,d_grid,a_grid,z_gridvals,pi_z,ReturnFn,FnsToEvaluate,FnsToEvaluateCell,GeneralEqmEqnsCell,Parameters,DiscountFactorParamNames,ReturnFnParamNames,FnsToEvaluateParamNames,GeneralEqmEqnParamNames,GEPriceParamNames,GEeqnNames,AggVarNames,nGEprices,heteroagentoptions,simoptions,vfoptions) (line 392)
        GeneralEqmConditionsFnOpt=@(p) HeteroAgentStationaryEqm_Case1_subfn(p, n_d, n_a, n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, FnsToEvaluate, FnsToEvaluateCell, GeneralEqmEqnsCell, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, GeneralEqmEqnParamNames, GEPriceParamNames, GEeqnNames, AggVarNames, nGEprices, heteroagentoptions, simoptions, vfoptions);
                                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in fminsearch (line 209)
f = funfcn(x,varargin{:});
    ^^^^^^^^^^^^^^^^^^^^^
Error in HeteroAgentStationaryEqm_Case1 (line 422)
        [p_eqm_vec,GeneralEqmConditions]=fminsearch(GeneralEqmConditionsFnOpt,GEparamsvec0,minoptions);
                                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Aiyagari1994 (line 103)
[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);
                          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

aledinola · February 3, 2026, 12:45pm

Thanks for testing the code. I ran the same example on CPU and I got the same mistake you reported.

The bug is in the function UnKronPolicyIndexes_InfHorz_CPU, line 18:

optA=ind2sub_homemade(n_aprime',optaindexKron);

The variable n_aprime does not exist, it should be the size of A (or the number of elements of A, not clear) . By the way, why is it transposed?

@robertdkirkby You can easily reproduce the error message even on a computer with GPU by setting

tauchenoptions.parallel=0;
vfoptions.parallel=0;

This reminds me of a previous “bug” I had found: if you want to test the CPU codes on a computer with GPU, it is not enough to set vfoptions.parallel=0. If you do this, z_grid and pi_z are still generated on the GPU and then the CPU codes are super-slow on GPU arrays (since they use loops).

aledinola · February 3, 2026, 1:00pm

@yechen Reading the message you reported, I see that you are using the Processes profile. If you are using a single computer, or a single node on a cluster, you might find ‘Threads’ to be faster. ‘Threads’ uses a shared memory parallelism (similar to OpenMP for C or Fortran), while ‘Processes’ uses distributed memory (similar to MPI).

https://www.mathworks.com/help/parallel-computing/choose-between-thread-based-and-process-based-environments.html

robertdkirkby · February 3, 2026, 7:06pm

I ran it on the CPU and received the following messages

Fixed. Somehow the copy of Aiyagari code I tested it on was a version that had been modified to have endogenous labor while I was playing round with something else Thanks for letting me know.

Still errors when it gets to AllStats, but this is a deliberate error message because AllStats is not supported for CPU. By this point you have already solved the Aiyagari model, just don’t get the same detailed model stats as you would with GPU. I will think about doing a version of this just for CPU just so the Aiyagari example can be run all the way through.

This reminds me of a previous “bug” I had found: if you want to test the CPU codes on a computer with GPU, it is not enough to set vfoptions.parallel=0. If you do this, z_grid and pi_z are still generated on the GPU and then the CPU codes are super-slow on GPU arrays (since they use loops).

Pretty sure this no longer happens due to rewrite that separates the CPU version off early. But not 100% certain (confident it doesn’t if you run value fn command, not confident if you run stationary general eqm command).

aledinola · February 4, 2026, 12:27am

I have tried to improve the Aiyagari example and sent a PR.

Among other fixes, I added a commented block that explains some vfoptions that users might find useful:

vfoptions  = struct();  % Default options for value function iteration
simoptions = struct();  % Default options for stationary distribution

% If you run out of memory, you can enable low-memory mode:
% vfoptions.lowmemory = 1;

% 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;

I have also added the command to convert policy from indexes to values, which is often useful to look at policy functions:

% Convert policy from indices to values
PolicyValues = PolicyInd2Val_Case1(Policy, n_d, n_a, n_z, d_grid, a_grid, vfoptions);

Another thing is a function for consumption defined in an M-file and values on grid:

FnsToEvaluate2.Consumption = @(aprime, a, z, alpha, delta, r) ...
    Aiyagari1994_ConsumptionFn(aprime, a, z, alpha, delta, r);

% Values on the grid
ValuesOnGrid = EvalFnOnAgentDist_ValuesOnGrid_Case1( Policy, FnsToEvaluate2, Params, [], n_d, n_a, n_z, d_grid, a_grid, z_grid, simoptions);

The purpose of FnsToEvaluate2.Consumption is also to show that you can set up FnsToEvaluate as M-files and not necessarily as one-line functions.

Finally, I’ve wrapped AllStats in

if gpuDeviceCount>0
    simoptions.npoints = 1000;
    AllStats = EvalFnOnAgentDist_AllStats_Case1(StationaryDist, Policy, FnsToEvaluate2, ...
        Params,[],n_d, n_a, n_z, d_grid, a_grid, z_grid, simoptions);
    disp(' ');
    disp('Distributions of Earnings and Wealth');
    fprintf('Gini (Earnings): %8.4f\n', AllStats.Earnings.Gini);
    fprintf('Gini (Income):   %8.4f\n', AllStats.Income.Gini);
    fprintf('Gini (Wealth):   %8.4f\n', AllStats.Wealth.Gini);
end

Overall, the main logic behind these additions is that

new users might struggle to find all these features since there isn’t a document collecting all options/commands/features, especially for infinite horizon models and
even more experienced users (like me!) would appreciate a place where they can look and immediately find what they need.
Users w/out supported GPU can run the entire file without getting any error message.

I hope this is useful

yechen · February 4, 2026, 3:05am

A big thanks! I’ve learned a lot.

aledinola · February 4, 2026, 1:58pm

@robertdkirkby Thanks a lot for merging my pull request!

Regarding the Guerrieri and Lorenzoni example, you write that it is hard-coded to run on GPU only. However, if you set vfoptions.gridinterplayer=0 the code runs also on CPU. I have tested now on my low-specs low-speed low-everything university laptop and it manages to compute the general equilibrium, stationary distribution and aggregate variables, until line 268.

Some information that might be useful to other users: with the exception of the transition, the Guerrieri and Lorenzoni (2017) model is relatively simple: it is a Huggett (1993) model with endogenous labor supply, so the action space is simply one d variable (labor supply), one a’ and one a variable and one z variable (idiosyncratic labor productivity), and this is coded on the CPU as well.