Questioning AgentDist_InfHorz_TPath_SingleStep

1

Here’s a small function I’m trying to understand–in particular, the Gammatranspose line.

function AgentDist=AgentDist_InfHorz_TPath_SingleStep(AgentDist,Policy_aprimez,II1,IIones,N_a,N_z,pi_z_sparse)
% II1 and IIones are inputs just so they can be precomputed and then reused each period of transition
% II1=gpuArray(1:1:N_a*N_z);
% IIones=ones(N_a*N_z,1,'gpuArray');
% AgentDist is [N_a*N_z,1]
% Policy_aprimez is either [2,N_a,N_z] if N_d>0 or [N_a,N_z] if N_d==0
% pi_z_sparse is [N_z,N_zprime]

% AgentDist is already sparse and on cpu

Gammatranspose=sparse(gather(Policy_aprimez),II1,IIones,N_a*N_z,N_a*N_z);

% Two steps of the Tan improvement
AgentDist=reshape(Gammatranspose*AgentDist,[N_a,N_z]);
AgentDist=reshape(AgentDist*pi_z_sparse,[N_a*N_z,1]);

end

In my case, N_d>0, so my Policy_aprimez is 2x201x3. But II1 and IIones are 1x603 and 603x1, respectively. The sparse constructor wants to see either half of Policy_aprimez or double the II1. Thoughts?

2

Sorry this comment

Policy_aprimez is either [2,N_a,N_z] if N_d>0 or [N_a,N_z] if N_d==0

is incorrect. Policy_aprimez does not depend on the size of d.

But you must be looking at an outdated version of that code, because the comment is not there,

PS. Gammatranspose essentially the Policy matrix for aprime, but converted into a transition matrix from (a,z) to (a’,z) [not a typo that ’ is on a but not on z, this is what the Tan improvement is about]

1 Like
3

Thanks for the answers! I added the comment based on the construction I observed. Will update when I get back to my computer.

4

There are some models where having d changes things, like if you have experienceasset/riskyasset or semi-exogenous state. But in all these situations you would end up in a different AgentDist SingleStep command (nProbs or SemiExo, respectively). The fact you are in AgentDist_InfHorz_TPath_SingleStep tells you that d cannot be relevant.

1 Like
5

Looks like it’s my fault. The call stack is:

From TransitionPath_Case1_FHorz_PType:

%% Shooting algorithm
if transpathoptions.GEnewprice~=2
    % For permanent types, there is just one shooting command,
    % because things like z,e, and fastOLG are handled on a per-PType basis (to permit that they differ across ptype)
    [PricePath,GEcondnPath]=TransitionPath_Case1_FHorz_PType_shooting(PricePath0, PricePathNames, ParamPath, ParamPathNames, T, V_final, AgentDist_init, jequalOneDist_T, AgeWeights_T, FnsToEvaluate, GeneralEqmEqns, PricePathSizeVec, ParamPathSizeVec, PricePathSizeVec_ii, ParamPathSizeVec_ii, use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, transpathoptions, PTypeStructure);

At this point we are actually getting answers for an InfHorz calculation (mixed-horizon), and I’m following literally the comment that “there is just one shooting command [operating across PTypes]”. In other words, I’m delegating down to TransitionPath_Case1_FHorz_PType_shooting the handling of the InfHorz case. Thus, from TransitionPath_Case1_FHorz_PType_shooting we get to:

        % For current ptype, do the backward iteration of V and Policy, then forward iterate agent dist and get the AggVarsPath
        if isfinite(PTypeStructure.(iistr).N_j)
            if vfoptions.experienceasset==1
                AggVarsPath=TransitionPath_FHorz_PType_ExpAsset_singlepath(PricePathOld_ii, ParamPath_ii, PricePathNames,ParamPathNames,T,V_final.(iistr),AgentDist_init.(iistr),jequalOneDist_T.(iistr),AgeWeights_T.(iistr),l_d,N_d,n_d,N_a,n_a,N_z,n_z,N_e,n_e,N_j,d_grid,a_grid,d_gridvals,aprime_gridvals,a_gridvals,z_gridvals_J, pi_z_J,pi_z_J_sim,e_gridvals_J,pi_e_J,pi_e_J_sim,ze_gridvals_J_fastOLG,ReturnFn, FnsToEvaluate, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, AggVarNames, PricePathSizeVec_ii, ParamPathSizeVec_ii, use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, II1orII, II2, exceptlastj,exceptfirstj,justfirstj, transpathoptions, vfoptions, simoptions);
            else
                AggVarsPath=TransitionPath_FHorz_PType_singlepath(PricePathOld_ii, ParamPath_ii, PricePathNames,ParamPathNames,T,V_final.(iistr),AgentDist_init.(iistr),jequalOneDist_T.(iistr),AgeWeights_T.(iistr),l_d,N_d,n_d,N_a,n_a,N_z,n_z,N_e,n_e,N_j,d_grid,a_grid,d_gridvals,aprime_gridvals,a_gridvals,z_gridvals_J, pi_z_J,pi_z_J_sim,e_gridvals_J,pi_e_J,pi_e_J_sim,ze_gridvals_J_fastOLG,ReturnFn, FnsToEvaluate, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, AggVarNames, PricePathSizeVec_ii, ParamPathSizeVec_ii, use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, II1orII, II2, exceptlastj,exceptfirstj,justfirstj, transpathoptions, vfoptions, simoptions);
            end
        else
            AggVarsPath=TransitionPath_InfHorz_PType_singlepath(PricePathOld_ii, ParamPath_ii, PricePathNames,ParamPathNames,T,V_final.(iistr),AgentDist_init.(iistr),l_d,N_d,n_d,N_a,n_a,N_z,n_z,d_grid,a_grid,d_gridvals,aprime_gridvals,a_gridvals,z_grid, pi_z,ReturnFn, FnsToEvaluate, Parameters, DiscountFactorParamNames, ReturnFnParamNames, FnsToEvaluateParamNames, AggVarNames, transpathoptions, vfoptions, simoptions);
        end
        % AggVarsPath=zeros(length(FnsToEvaluate),T-1);

TransitionPath_InfHorz_PType_singlepath is a function I wrote based on extracting the guts of TransitionPath_InfHorz_PType_shooting. Its purpose is to deliver the AggVarsPath result for TransitionPath_Case1_FHorz_PType to weave into the unified mixed horizon results.

This is the inner loop taken from that, with all mistakes mine (especially the pre-computation of II1 and IIones, which I inferred by reading other codes):

% Precompute for later
II1=gpuArray(1:1:N_a*N_z);
IIones=ones(N_a*N_z,1,'gpuArray');

for tt=1:T-1
    
    %Get the current optimal policy
    if N_d>0
        Policy=PolicyIndexesPath(:,:,:,tt);
    else
        Policy=PolicyIndexesPath(:,:,tt);
    end
    
    GEprices=PricePathOld(tt,:);
    
    % Again, since we don't do tminus1, don't dig as deep as FHorz
    for nn=1:length(ParamPathNames)
        Parameters.(ParamPathNames{nn})=ParamPath(tt,nn);
    end
    for nn=1:length(PricePathNames)
        Parameters.(PricePathNames{nn})=PricePathOld(tt,nn);
    end
    
    PolicyUnKron=UnKronPolicyIndexes_Case1(Policy, n_d, n_a, n_z,vfoptions);
    AggVars=EvalFnOnAgentDist_AggVars_Case1(AgentDist, PolicyUnKron, FnsToEvaluate, Parameters, FnsToEvaluateParamNames, n_d, n_a, n_z, d_grid, a_grid, z_grid, simoptions); % The 2 is for Parallel (use GPU)
    
    AgentDist=AgentDist_InfHorz_TPath_SingleStep(AgentDist,Policy,II1,IIones,N_a,N_z,sparse(pi_z));

    AggVarsPath(:,tt)=AggVars;
end

It is in that loop that we have the N_d context arises.

6

Hmm, I feel like things are going to get messy.

Maybe toolkit needs a ‘MixHorz’ alongside ‘FHorz’ and ‘InfHorz’? What do you think @MichaelTiemann? I feel like that will be easier/cleaner/more intuitive than the current approach of hiding ‘MixHorz’ inside ‘FHorz PType’. Obviously ‘MixHorz’ will only have ‘PType’ commands (as you cannot mix horizons without having at least two permanent types).

This has two main advantages to my mind. First is that it will be more obvious/intuitive to the user. Second is that it keeps all the InfHorz stuff out of the FHorz PType commands. Should be easy enough to construct as will basically just be a matter of creating the MixHorz commands as a copy-merge of the InfHorz and FHorz commands.

[question is about how to handle models when there are both FHorz agents and InfHorz agents in the same model, e.g., an OLG model of households together with infinitely-lived firms.]

PS. Obviously we can do this by just copy-pasting what you have done so far and renaming it MixHorz PType, and then just leaving the FHorz PType commands as they are/were. I would be a lot more comfortable with this as then we know the FHorz is not being impacted, and I am happy to have MixHorz PType as somewhat buyer-beware for now [I think of all features as buyer-beware for the first 6-12 months they exist, until they get decent amount of use that tests them out.]

1 Like
7

There’s no free lunch. I can see that you have leaned into the idea of replicating code for each interesting case, which means that each interesting case has the simplest possible code to implement it. On the other hand, across all the replication forks, there’s a lot of code that needs to be updated when something upstream in the logic changes. For pure math that’s provably true, we don’t expect upstream logic changes, but when it comes to code…it’s a lot of maintenance work.

On the other hand, I’m pretty sure there was an original case of FHorz being the principal when there’s an FHorz/InfHorz split (because FHorz has the J terms and InfHorz does not). A mixed horizon can be supported with logic that looks at whether then J in question is finite or not, and then do something special/delegate as needed.

It’s not that hard to create a test early on that will delegate to FHorz, InfHorz, or MixedHorz. But maintaining it would require that a test suite is updated to ensure that regardless of where the logic gets changed, all the cases are either properly affected or can properly ignore the change.

No rush on the decision, but something to think about. Personally, I think the code would benefit from a few cases of gathering more common things together and making the logic more complex internally, but given how many special cases there are, I can imagine that bringing that control logic into one place will quickly obscure everything else. So it is a balance.

My changes are an attempt to find the natural splitting/grafting point that is a happy middle ground. But it’s ultimately your call.

On the original topic, I’m going to change the construction of II1 and IIones to suit the dimensions of the Policy shape. Nothing to do with aprimez per se.

1 Like
8

No rush on the decision, but something to think about. Personally, I think the code would benefit from a few cases of gathering more common things together and making the logic more complex internally, but given how many special cases there are, I can imagine that bringing that control logic into one place will quickly obscure everything else. So it is a balance.

This tends to happen gradually in the toolkit as I discover what the ‘general case’ is. You can see a lot of it in the FHorz commands, where things like setting up the shocks or figuring out the ReturnFn inputs are all sent off to a function (so then just edit that function once to change things everywhere). This hasn’t yet happened in TPath as that is still evolving and I am still finding and understanding newer and better ways to set things up, but it will gradually happen.

For example, the way that exogenous shocks are now handled took a few years to discover. I kept adding things —stacked-columns, joint-grids, as a function— while trying to do new things and then realized that this could all be handled together in one function and just sets it all up and creates a final joint-grid that all codes can then just assume and use. Point is, creating those ‘general functions’ takes quite a lot of time, not in the actual coding, but in understanding what they should do.

Also, nowadays the only function that really has to be coded for each individual case is the value function commands. The agent dist and model stats commands are all pretty flexible in terms of handling a wide variety of cases.

Given the desire for each case to be as fast as possible, conditional on not making too strong assumptions about the model, I think hard-coding the individual cases is just unavoidable. Once you understand the structures copy-pasting with minor edits to create all the different combinations is actually a lot easier and faster than it appears at first glance.

PS. Policy_aprimez is a name to remind me that it has to go to (a’,z), so it is not just the part of Policy that corresponds to aprime, which I would call Policy_aprime.

PPS. The code that handles PType and creates all the _temp (for the ‘current ptype’) is definitely one of the bits that could do with being replaced by a general function or three.

1 Like
9

Very interesting discussion. I remember in the early days of the toolkit (back in 2019?) there was only VFI case 1 and case 2, based on a classification in Stokey-Lucas :blush:

2 Likes
10

Yeah, and nowadays I am trying to eliminate the Case 1 vs Case 2 distinction :rofl:
[turns out when doing the compute, instead of the math theory, it is better to think about it as a property of an endogenous state, rather than as a property of the model as a whole]

1 Like
11

I’m curious about the difference between VFI Case 1 and Case 2. In my limited experience, I’ve never come across an example of Case 2 being used.

1 Like
12

In Stokey Lucas and Prescott (1989) they distinguish value fn problems into Case 1 and Case 2. In Case 1 you can directly choose next period endogenous state. In Case 2 you cannot directly choose next period endogenous state, it is a function of some things you can choose. In SLP1989, you then make assumptions on the function (next period in terms of some decision) and can do a bunch of theory results from this. Originally VFI Toolkit had Case1 and Case2 commands, that handled these (I had recently read SLP1989 so was just building on their theory when thinking how to construct models in the code).

But it became clear as VFI Toolkit progressed that many models contain a combination of these two kinds of endogenous state.

Nowadays Case 1 is just what VFI Toolkit calls a standard endogenous state. There are many different Case 2 instances, including experienceasset, riskyasset, and more.

It is common for models to have, e.g., one standard endogenous state and one experienceasset. So the toolkit now works by sending all these to the ‘Case1’ command, and just using vfoptions to note that the second endogenous state is an experience asset. Conceptually the change is that the ‘case’ is not a property of the model, it is a property of each endogenous state.

So nowadays all of these get sent to a command that for legacy reasons is called ‘Case1’, but is actually handling everything. There are still a few commands in the toolkit called Case2 but they are rarely if ever used.

My long term plan is to just delete all the ‘Case2’ commands, and remove the word ‘Case1’ from the description of the ‘Case1’ commands.

[Originally I called them Case2-TypeA, Case2-TypeB, etc., which was always impossible to remember but now I use more easy to remember names like experienceasset, riskyasset, based on their main applications. The Castaneda, Rios-Rull & Diaz-Gimenez (2003) replication used to be the most prominent use of Case2 but I recently rewrote it introducing ‘inheritanceasset’ as part of readying the toolkit to eliminate Case2 entirely.]

3 Likes
13

I see. Thanks for taking the time to explain so kindly.

2 Likes