Error with gpuArray/arrayfun and POWER operation in GourinchasParker2002.m

1

Hello Robert,

I am a user of the VFI Toolkit and am trying to replicate the main results from the GourinchasParker2002.m example.

The script runs for a while, but eventually terminates with a gpuArray/arrayfun error related to a POWER operation. The full error message and call stack are as follows (matlab2004b)

可能存在局部最小值。
lsqnonlin 已停止，因为当前步长小于
步长容差的值。

错误使用 gpuArray/arrayfun
POWER: needs to return a complex result, but this is not supported for real input X and Y on the GPU. Use POWER(COMPLEX(X),COMPLEX(Y,0)) instead.

出错 CreateReturnFnMatrix_Case1_Disc_DC1_nod_Par2e (第 33 行)
Fmatrix=arrayfun(ReturnFn, aprime_grid, shiftdim(a_grid,-1), shiftdim(z_gridvals(:,1),-2), shiftdim(e_gridvals(:,1),-3), ParamCell{:});
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
出错 ValueFnIter_FHorz_DC1_nod_e_raw (第 38 行)
ReturnMatrix_ii=CreateReturnFnMatrix_Case1_Disc_DC1_nod_e_raw(ReturnFn, n_z, n_e, a_grid, a_grid(level1ii), z_gridvals_J(:,:,N_j), e_gridvals_J(:,:,N_j), ReturnFnParamsVec,1);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
… (rest of the call stack)

Analysis:

Based on the error message, it appears that a POWER operation (e.g., x^y) is being applied on the GPU, where the base x is a negative real number and the exponent y is a non-integer. This would require a complex result, which the GPU does not support in this context.

The traceback points to CreateReturnFnMatrix_Case1_Disc_DC1_nod_Par2e calling ReturnFn. My hypothesis is that this error originates from the GourinchasParker2002_ReturnFn.m function, likely when the consumption c = R*a + z*e - aprime becomes a very small or negative value, causing the utility calculation c^(1-rho) to fail, as 1-rho is a non-integer.

The lsqnonlin warning about a potential local minimum also suggests that the optimization process is exploring regions of the parameter space that may lead to these numerical issues.

Question:

Could you please advise on how to handle this? Is there a common fix for this issue within the GP2002 model? I suspect the problem is that some combination of parameters or grid points is leading to negative consumption.

Is there a recommended way to:

  1. Adjust the grids (a_grid, z_grid) to prevent negative consumption?
  2. Add a small value (e.g., epsilon) to consumption, such as (c + 1e-6)^(1-rho), to prevent it from going to zero or negative?

Any guidance on how to resolve this numerical instability would be greatly appreciated.

Thank you very much for your help!

2 Likes
2

You are correct on the cause of the error, that there is a negative number to a power, and on the GPU this causes an error.

Less clear to me is why this happened. Last time I ran it was early this year and it worked fine at the time. I will run this myself to check and fix, but I am on holiday all next week so will take near two weeks before I can reply to this.

PS. If you look in the return fn, you will see that is sets F=-Inf and then only evaluates F=utility(c) if c>0. So it won’t be because of negative c (that just gives negative infinity utility, so no-one chooses it as a policy).

PPS. I anyway need to revisit this code to turn on grid interpolation layer so will do that while I am there.

2 Likes
3

There is also the term Fret=betakappaupsilon*(XplusH^(1-rho))/(1-rho)

Could it be that XplusH is negative for some values of the inputs to ReturnFn?

2 Likes
4

@Zhang_Lu Does the error message happen while running this line of GourinchasParker2022.m:

[EstimParams_robust, EstimParamsConfInts_robust, estsummary_robust]=EstimateLifeCycleModel_MethodOfMoments(EstimParamNames,TargetMoments, RobustWeightingMatrix,CovarMatrixDataMoments, n_d,n_a,n_z,N_j,d_grid, a_grid, z_grid_J, pi_z_J, ReturnFn, Params, DiscountFactorParamNames, jequaloneDist,AgeWeightParamNames, ParametrizeParamsFn, FnsToEvaluate, estimoptions, vfoptions,simoptions);

?

I ran the code on my laptop (also R2024b) and initially I did not get any error. I stopped the estimation code after a while though. It could be that the minimization routine lsqnonline selected some strange parameter values that generate negative numbers raised to a power. My last set of param values before terminating the program is:

Current (sum-of-squares of) objective fn value is 1.545497742840 
 
Current parameter values: 
    beta= 0.963002 
    rho= 0.313389 
    h= 470.923784 
    kappa= 2.537381 
    rhoJ= 0.305273
1 Like
5

Hi @Zhang_Lu, I also got the same mistake. My script runs for a while ok, but then it breaks and display the error that you reported

@robertdkirkby How do I turn on grid interpolation? Does it speed up the code?

1 Like
6

You should just be able to add the following four lines to your code:
vfoptions.gridinterplayer=1;
vfoptions.ngridinterp=10;
simoptions.gridinterplayer=vfoptions.gridinterplayer;
simoptions.ngridinterp=vfoptions.ngridinterp;

Until now, the choice of next period assets, a’, was restricted to be on the same grid as this period, so on the a_grid. These options turn on ‘grid interpolation layer’ and use 10 points (the 10 above can be any number). The grid interpolation layer means that now a’ considers an additional 10 evenly-spaced points between every two consecutive points in a_grid.

Technically, turning grid interpolation on for any given size of a_grid (for any given n_a) will be slower. But of course the point is that now it is interpolating extra points and will be more accurate, so you can choose a smaller n_a. With both a smaller n_a and grid interpolation layer at the same time, you should get something that is both faster and more accurate.

The interpolation of the ReturnFn is exact (so model still satisfies all constraints), and the interpolation of next period value fn is done as linear interpolation.

PS. Grid interpolation layer currently works for one endogenous state models (except transition paths), and also for models with two endogenous states where the second state is an experienceasset. It is slowly rolling out to everything but will take time.

1 Like
7

I have figured out the cause of the error (interaction between taking derivatives and the edges of the constraints). Hope to have it fixed later this week.

2 Likes
8

I think grid interpolation works also for models with two generic endogenous states (steady state).

9

For two standard endogenous states. Grid interpolation layer works for InfHorz. It works for some some FHorz setups, not all (I haven’t finished doing all shock combos).

2 Likes
10

I fixed what was causing this error. But the command itself is not yet fixed. I realised that there is an issue I had not previously considered, namely “how to compute the Jacobian (and hence standard errors and confidence intervals) when the parameter is up against the constraint”. I will need to go read up on what I am supposed to do in this case so will take some time, hopefully fixed next week.

In the meantime, the GMM estimation command is working, except if you end up with the parameter constraints binding on your final estimated parameter values, in which case it will still error.

2 Likes