Bug in LifeCycleProfiles_FHorz_Case1_PType

[aledinola]

The call to function LifeCycleProfiles_FHorz_Case1_PType returned this error:

Unrecognized function or variable 'tempStats2'.

Error in LifeCycleProfiles_FHorz_Case1_PType (line 1127)
                                AgeConditionalStats.(CondlRestnFnNames{rr}).(FnsToEvalNames{ff}).Mean(jj)=tempStats2.Mean;

If I remove the field conditionalrestrictions from simoptions, then the function runs fine.

By the way, EvalFnOnAgentDist_AllStats_FHorz_Case1_PType works fine even with conditional restrictions, but there is a weird display on the screen:

Parameter values for the current permanent type 

Parameters_temp = 

  struct with fields:

         inp_dir: 'C:\Users\aledi\Dropbox\1 - RESEARCH PROJECTS\18_health_search\shared_all\Codes\v15\inputs'
               J: 76
       age_entry: 25
        age_exit: 100
            agej: [76×1 double]
      age_actual: [76×1 double]
              Jr: 41
            mewj: [76×1 double]
            beta: 0.9800
               r: 0.0400
           sigma: 2
          kappa1: 1.2500
        uc_level: 0.0100
            tau0: 0.1700
            tau1: 0.0620
             pen: 0.2500
           pen_j: 0.2500
             n_a: 200
             n_d: 5
             n_z: 7
           a_min: 0
           a_max: 600
     a_curvature: 3
       n_theta_e: 2
       n_theta_h: 2
       n_theta_y: 2
      lambda0_e0: 0.8500
      lambda1_e0: 0.0450
      lambda2_e0: -5.5000e-04
      lambda0_e1: 1.0500
      lambda1_e1: 0.0500
      lambda2_e1: -6.0000e-04
        wp_h0_e0: 0.2000
        wp_h0_e1: 0.1500
    theta_y_vals: [2×1 double]
    omega_bar_e0: 0.7000
    omega_bar_e1: 0.8500
     delta_h0_e0: [76×1 double]
     delta_h1_e0: [76×1 double]
     delta_h0_e1: [76×1 double]
     delta_h1_e1: [76×1 double]
      kappa_h0_j: [76×1 double]
      kappa_h1_j: [76×1 double]
           rho_z: 0.9850
     sigma_eps_z: 0.1483
              p0: 0.1500
     ptype_probs: 0.2973
     theta_e_val: 1
     theta_h_val: 1
     theta_y_val: 0.1200
         N_semiz: 42

[aledinola]

I have sent a PR that should fix the bug. But the code was complicated, so please @robertdkirkby have a look

[robertdkirkby]

@aledinola

By the way, EvalFnOnAgentDist_AllStats_FHorz_Case1_PType works fine even with conditional restrictions, but there is a weird display on the screen

Was a typo in setting simoptions defaults. Fixed now.

[robertdkirkby]

I ran LifeCycleProfiles_FHorz_Case1_PType with simoptions.conditionalrestrictions and didn’t get any error about Unrecognized function or variable 'tempStats2'.

Are you able to email me a code that generates this email?

[robertdkirkby]

I added LifeCycleProfiles and simoptions.conditionalrestrictions to the example in

Semi-exogenous shocks with permanent types

I’m done with the semi-exogenous states for permanent types. You must set up vfoptions.semiz_grid, and then one of vfoptions.SemiExoShockFn or vfoptions.pi_semiz. You can include age and ptype as dimensions of semiz_grid and pi_semiz, or not as you wish. Dependence on ptype can be done as a (the last) dimension, or as fieldnames. Here is the code that I was checking all the combos. % Age % One endogenous variable: assets % One stochastic exogenous variables: income shock J=10; % Ages 1 to 1…

works just fine there.

[robertdkirkby]

Copy-pasting the following from a different topic of the same title, no need for two topics on essentially same issue.

---

I think there is a bug in

EvaluateFnOnAgentDist/FHorz/PType/LifeCycleProfiles_FHorz_Case1_PType.m

when all of these are used together:

• PType

• simoptions.conditionalrestrictions

• simoptions.groupptypesforstats = 1

I get:

Error using accumarray

Second input must be a vector with one element for each row in the first input, or a scalar.

Error in LifeCycleProfiles_FHorz_Case1_PType (line 1117)

AllRestrictedWeights_rrffjj=accumarray(sortindex,AllRestrictedWeights_rrffjj,[],@sum);

The issue seems to be this:

Earlier, in the grouped unconditional-stats block, the code does

[AllValues.(FnsToEvalNames{ff}).(jgroupstr{jj}),~,sortindex]=unique(AllValues.(FnsToEvalNames{ff}).(jgroupstr{jj}));

AllWeights.(FnsToEvalNames{ff}).(jgroupstr{jj})=accumarray(sortindex,AllWeights.(FnsToEvalNames{ff}).(jgroupstr{jj}),[],@sum);

so AllValues.(...) gets overwritten by its unique support points.

Later, in the grouped conditional-restriction block, the code does

AllValues_rrffjj=gather(AllValues.(FnsToEvalNames{ff}).(jgroupstr{jj})(:));

AllRestrictedWeights_rrffjj=gather(AllRestrictedWeights.(CondlRestnFnNames{rr}).(FnsToEvalNames{ff}).(jgroupstr{jj})(:));

[AllValues_rrffjj,~,sortindex]=unique(AllValues_rrffjj);

AllRestrictedWeights_rrffjj=accumarray(sortindex,AllRestrictedWeights_rrffjj,[],@sum);

But at that point:

• AllValues.(...) has already been collapsed by unique(...)

• AllRestrictedWeights.(...) is still at the old pre-collapsed length

So sortindex and AllRestrictedWeights_rrffjj no longer have matching lengths, which seems to explain the accumarray error.

My guess is that the grouped unconditional-statistics block should use temporary variables rather than overwrite AllValues.(...), or else the grouped conditional-restriction block needs to rebuild matching values/weights from the same pre-collapsed objects.

This looks like a generic bookkeeping bug in the grouped-PType conditional-restriction path…

[robertdkirkby]

I ran this combo:

• PType
• simoptions.conditionalrestrictions
• simoptions.groupptypesforstats = 1
  and didn’t get any error.

@aledinola if you have a code that creates this error please email me a copy.

[aledinola]

Thanks for looking into this. I will send you the code tomorrow.

[aledinola]

I synced my github and the error disappeared. However, the call to LifeCycleProfiles_FHorz_Case1_PType takes a huge amount of time relative to the rest of the model. See AgeStats which is the output of LifeCycleProfiles_FHorz_Case1_PType

Solve model timings (seconds):
  ValueFnIter        2.530
  PolicyInd2Val      0.029
  StationaryDist     0.489
  ValuesOnGrid       0.229
  AllStats           0.543
  AgeStats          49.428
  Total             53.249

I have 8 permanent types and a few conditional restrictions; on the other hand I chose to calculate only the means, so I am a bit surprised the code is still so slow.

simoptions = 

  struct with fields:

                    verbose: 0
                 whichstats: [7×1 double]
            gridinterplayer: 0
                ngridinterp: 0
    conditionalrestrictions: [1×1 struct]

>> simoptions.whichstats

ans =

     1
     0
     0
     0
     0
     0
     0

I shared the code with you on Dropbox. You can see that I have replaced LifeCycleProfiles_FHorz_Case1_PType with my own function fun_agestats which calculated only the means and it takes only 0.6 seconds.
I think LifeCycleProfiles_FHorz_Case1_PType does something redundant:

[SortedValues_jj,~,sortindex]=unique(Values_jj);
                    SortedWeights_jj=accumarray(sortindex,StationaryDistVec_jj,[],@sum);

                    SortedWeights_jj=SortedWeights_jj/sum(SortedWeights_jj(:)); % Normalize conditional on jj (is later renormalized ii weight before storing for groupstats)

                    %% Use the full ValuesOnGrid_ii and StationaryDist_ii to calculate various statistics for the current PType-FnsToEvaluate (current ii and ff)
                    tempStats=StatsFromWeightedGrid(SortedValues_jj,SortedWeights_jj,simoptions.npoints,simoptions.nquantiles,simoptions.tolerance,1,simoptions.whichstats); % 1 is presorted

If I want to compute only the means (i.e. the aggregate values, E(X) = sum(x.*statdist)) I do not need to sort the big arrays (right?), but the code above calls the unique function which is expensive.

For me improving the performance of LifeCycleProfiles_FHorz_Case1_PType is NOT a priority however: most of the time I replace it with my own function(s). There are also other tricks I haven’t used: in practice I need age stats on only a subset of functions to evaluate, so I could create a smaller FntoEval2 and feed that to LifeCycleProfiles_FHorz_Case1_PType. Same goes for conditional restrictions.

Solve model timings (seconds):
  ValueFnIter        2.552
  PolicyInd2Val      0.031
  StationaryDist     0.460
  ValuesOnGrid       0.249
  AllStats           0.544
  AgeStats           0.073
  Total              3.908

Thanks for fixing the code!!

[robertdkirkby]

By default LifeCycleProfiles with PType is setting simoptions.whichstats=[1,1,1,2,1,2,1], but if you use it in a calibration command and only target means, the calibration command should automatically set it to simoptions.whichstats=[1,0,0,0,0,0,0] which will substantially cut runtimes. [=0 is ‘do not calculate’, =2 is ‘calculate in the faster but more memory intensive manner’]

From what I understand though you are saying even if you set simoptions.whichstats=[1,0,0,0,0,0,0] it is still taking a long time, versus hardcoding it yourself.

As you also mention, the other trick is to pass version of FnsToEvaluate that contains only you need rather age-condtional stats for rather than all the FnsToEvaluate, this slashes the runtimes.

Toolkit commands will of course never be as fast as hardcoding things yourself. But they are very easy to use

I know how I could substantially cut the runtimes of LifeCycleProfiles by parallel over j rather than loop over j, but would be substantially more memory intensive. I haven’t got round to doing it as only in last year or two are GPUs getting powerful enough that this would be an appealing approach (without just running out of memory).

[aledinola]

Thanks for the nice explanation. I think in my case what slows things down is that I have many ptypes and many conditional restrictions. I cannot avoid conditional restrictions because I have to condition on being alive (I have health defined as h=0,1,2, where 0-1 refer to bad and good health and h=2 is dead)