From 358d2728f65e7cf74429d9f2e532d83e39617564 Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Wed, 25 Mar 2026 06:34:21 +1300
Subject: [PATCH 1/4] Remaining PType and TPath changes (no ExpAsset)

These are the non-ExpAsset changes needed thus far to progress a model through all the ptype and tpath things I have encountered.
---
 .../FHorz/LifeCycleProfiles_FHorz_Case1.m     |   2 +-
 ...FnOnAgentDist_AllStats_FHorz_Case1_PType.m |  83 ++-
 .../EvalFnOnTransPath_AggVars_Case1.m         |  13 +-
 ...lFnOnTransPath_AggVars_Case1_FHorz_PType.m |   7 +-
 .../HeteroAgentStationaryEqm_Case1_FHorz.m    |   1 +
 ...teroAgentStationaryEqm_Case1_FHorz_PType.m |  63 +--
 ...entStationaryEqm_Case1_FHorz_PType_subfn.m |  35 +-
 .../HeteroAgentStationaryEqm_Case1_PType.m    |   4 +-
 ...teroAgentStationaryEqm_Case1_PType_subfn.m |  16 +-
 Optimization/setupGEnewprice3_shooting.m      |  17 +-
 .../TransPath/AgentDistOnTransPath_Case1.m    |   2 +-
 .../AgentDistOnTransPath_Case1_FHorz_PType.m  |  26 +-
 .../PType/TransitionPath_Case1_FHorz_PType.m  | 519 ++++++++++++------
 ...ransitionPath_Case1_FHorz_PType_shooting.m | 100 +++-
 ...nsitionPath_FHorz_PType_singlepath_e_raw.m |   7 +-
 ...ath_FHorz_PType_singlepath_fastOLG_e_raw.m |   5 +-
 ...ransitionPath_FHorz_PType_singlepath_raw.m |   6 +-
 .../TransitionPath_InfHorz_PType_shooting.m   |   1 +
 .../TransitionPath_InfHorz_PType_singlepath.m | 168 ++++++
 .../AgentDist_InfHorz_TPath_SingleStep.m      |   2 +-
 .../ValueFnIter_InfHorz_TPath_SingleStep.m    |  10 +-
 .../PricePathParamPath_FHorz_StructToMatrix.m |  29 +-
 .../Subcodes/inputsFindtplus1tminus1.m        |  83 ++-
 .../ValueFnIter_nod_HowardGreedy_raw.m        |   9 -
 .../TransPath/ValueFnOnTransPath_Case1.m      |   2 +-
 .../ValueFnOnTransPath_Case1_FHorz_PType.m    |   6 +-
 26 files changed, 885 insertions(+), 331 deletions(-)
 create mode 100644 TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m

diff --git a/EvaluateFnOnAgentDist/FHorz/LifeCycleProfiles_FHorz_Case1.m b/EvaluateFnOnAgentDist/FHorz/LifeCycleProfiles_FHorz_Case1.m
index 38019cf3..1008bc95 100644
--- a/EvaluateFnOnAgentDist/FHorz/LifeCycleProfiles_FHorz_Case1.m
+++ b/EvaluateFnOnAgentDist/FHorz/LifeCycleProfiles_FHorz_Case1.m
@@ -463,7 +463,7 @@
             for jj=j1:jend
                 % Includes check for cases in which no parameters are actually required
                 FnToEvaluateParamsCell=CreateCellFromParams(Parameters,FnsToEvaluateParamNames(ff).Names,jj);
-                Values(:,:,jj-j1+1)=EvalFnOnAgentDist_Grid(FnsToEvaluate{ff}, FnToEvaluateParamsCell,PolicyValuesPermuteJ(:,:,:,jj),l_daprime,n_a,n_z,a_gridvals,z_gridvals_J(:,:,jj));                
+                Values(:,:,jj-j1+1)=EvalFnOnAgentDist_Grid(FnsToEvaluate{ff}, FnToEvaluateParamsCell,PolicyValuesPermuteJ(:,:,:,jj),l_daprime,n_a,n_z,a_gridvals,z_gridvals_J(:,:,jj));
             end
 
             Values=reshape(Values,[N_a*N_z*(jend-j1+1),1]);
diff --git a/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m b/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m
index 7ed93d53..1d1a782a 100644
--- a/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m
+++ b/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m
@@ -303,18 +303,35 @@
     N_a_temp=prod(n_a_temp);
     
     a_gridvals_temp=CreateGridvals(n_a_temp,a_grid_temp,1);
-    % Turn (semiz,z,e) into z_gridvals_J_temp as FnsToEvalute do not distinguish them
-    [n_z_temp,z_gridvals_J_temp,N_z_temp,l_z_temp,simoptions_temp]=CreateGridvals_FnsToEvaluate_FHorz(n_z_temp,z_grid_temp,N_j_temp,simoptions_temp,Parameters_temp);
+    if isfinite(N_j_temp)
+        % Turn (semiz,z,e) into z_gridvals_J_temp as FnsToEvalute do not distinguish them
+        [n_z_temp,z_gridvals_J_temp,N_z_temp,l_z_temp,simoptions_temp]=CreateGridvals_FnsToEvaluate_FHorz(n_z_temp,z_grid_temp,N_j_temp,simoptions_temp,Parameters_temp);
+    else
+        l_z_temp=length(n_z_temp);
+        N_z_temp=prod(n_z_temp);
+        z_gridvals_temp=CreateGridvals(n_z_temp,z_grid_temp,1);
+    end
     if N_z_temp==0
         N_z_temp=1; % Just makes things easier below
     end
 
     % Switch to PolicyVals
-    PolicyValues_temp=PolicyInd2Val_FHorz(PolicyIndexes_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,d_grid_temp,a_grid_temp,simoptions_temp,1);
-    if l_z_temp==0
-        PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,1]); % (N_a,N_j,l_daprime)
+    if isfinite(N_j_temp)
+        PolicyValues_temp=PolicyInd2Val_FHorz(PolicyIndexes_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,d_grid_temp,a_grid_temp,simoptions_temp,1);
+        if l_z_temp==0
+            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,1]); % (N_a,N_j,l_daprime)
+        else
+            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,4,1]); % (N_a,N_z,N_j,l_daprime)
+        end
+        StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
     else
-        PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,4,1]); % (N_a,N_z,N_j,l_daprime)
+        PolicyValues_temp=PolicyInd2Val_Case1(PolicyIndexes_temp,n_d_temp,n_a_temp,n_z_temp,d_grid_temp,a_grid_temp,simoptions_temp,1);
+        if l_z_temp==0
+            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,1]); % (N_a,l_daprime)
+        else
+            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,1]); % (N_a,N_z,l_daprime)
+        end
+        StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
     end
     l_daprime_temp=size(PolicyValues_temp,1);
     
@@ -322,7 +339,6 @@
     FnsAndPTypeIndicator(:,ii)=FnsAndPTypeIndicator_ii;
 
     %% Some things that don't need to go in the loop over FnsToEvalaute
-    StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
     % Eliminate all the zero-weighted points (this doesn't really save runtime for the exact calculation and often can increase it, but
     % for the createDigest it slashes the runtime. So since we want it then we may as well do it now.)    
     temp=logical(StationaryDist_ii~=0);
@@ -348,14 +364,21 @@
                 CondlRestnFnParamNames={};
             end
 
-            if l_z_temp==0
-                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,2); % j in 2nd dimension: (a,j,l_d+l_a), so we want j to be after N_a
-                RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,0,a_gridvals_temp,[]));
+            if isfinite(N_j_temp)
+                if l_z_temp==0
+                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,2); % j in 2nd dimension: (a,j,l_d+l_a), so we want j to be after N_a
+                    RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,0,a_gridvals_temp,[]));
+                else
+                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,3); % j in 3rd dimension: (a,z,j,l_d+l_a), so we want j to be after N_a and N_z
+                    RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp));
+                end
+                RestrictionValues=reshape(RestrictionValues,[N_a_temp*N_z_temp*N_j_temp,1]);
             else
-                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,3); % j in 3rd dimension: (a,z,j,l_d+l_a), so we want j to be after N_a and N_z
-                RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp));
+                CondlRestnFnParamsCell=CreateCellFromParams(Parameters_temp,CondlRestnFnParamNames);
+
+                RestrictionValues=logical(EvalFnOnAgentDist_Grid(CondlRestnFn, CondlRestnFnParamsCell,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_temp));
+                RestrictionValues=reshape(RestrictionValues,[N_a_temp*N_z_temp,1]);
             end
-            RestrictionValues=reshape(RestrictionValues,[N_a_temp*N_z_temp*N_j_temp,1]);
 
             RestrictedStationaryDistVec=StationaryDist_ii;
             RestrictedStationaryDistVec(~RestrictionValues)=0; % zero mass on all points that do not meet the restriction
@@ -383,24 +406,36 @@
         if FnsAndPTypeIndicator_ii(ff)==1 % If this function is relevant to this ptype
 
             % Get parameter names for current FnsToEvaluate functions
-            tempnames=getAnonymousFnInputNames(FnsToEvaluate.(FnsToEvalNames{ff}));
+            % Get parameter names for current FnsToEvaluate functions
+            if isstruct(FnsToEvaluate.(FnsToEvalNames{ff}))
+                tempfn=FnsToEvaluate.(FnsToEvalNames{ff}).(Names_i{ii});
+            else
+                tempfn=FnsToEvaluate.(FnsToEvalNames{ff});
+            end
+            tempnames=getAnonymousFnInputNames(tempfn);
             if length(tempnames)>(l_daprime_temp+l_a_temp+l_z_temp)
                 FnsToEvaluateParamNames={tempnames{l_daprime_temp+l_a_temp+l_z_temp+1:end}}; % the first inputs will always be (d,aprime,a,z)
             else
                 FnsToEvaluateParamNames={};
             end
-            if l_z_temp==0
-                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,2);
+            if isfinite(N_j_temp)
+                if l_z_temp==0
+                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,2);
+                else
+                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,3);
+                end
+                
+                %% We have set up the current PType, now do some calculations for it.
+                simoptions_temp.keepoutputasmatrix=1;
+                ValuesOnGrid_ii=EvalFnOnAgentDist_Grid_J(tempfn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp);
+                ValuesOnGrid_ii=reshape(ValuesOnGrid_ii,[N_a_temp*N_z_temp*N_j_temp,1]);
+    
+                % StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
             else
-                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,3);
+                FnToEvaluateParamsCell=CreateCellFromParams(Parameters_temp,FnsToEvaluateParamNames);
+                ValuesOnGrid_ii=EvalFnOnAgentDist_Grid(tempfn, FnToEvaluateParamsCell,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_temp);
+                ValuesOnGrid_ii=reshape(ValuesOnGrid_ii,[N_a_temp*N_z_temp,1]);
             end
-            
-            %% We have set up the current PType, now do some calculations for it.
-            simoptions_temp.keepoutputasmatrix=1;
-            ValuesOnGrid_ii=EvalFnOnAgentDist_Grid_J(FnsToEvaluate.(FnsToEvalNames{ff}),CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp);
-            ValuesOnGrid_ii=reshape(ValuesOnGrid_ii,[N_a_temp*N_z_temp*N_j_temp,1]);
-
-            % StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
 
             % Eliminate all the zero-weighted points (this doesn't really save runtime for the exact calculation and often can increase it, but
             % for the createDigest it slashes the runtime. So since we want it then we may as well do it now.)
diff --git a/EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_AggVars_Case1.m b/EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_AggVars_Case1.m
index 5ed6b654..74a05423 100644
--- a/EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_AggVars_Case1.m
+++ b/EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_AggVars_Case1.m
@@ -1,6 +1,17 @@
-function AggVarsPath=EvalFnOnTransPath_AggVars_Case1(FnsToEvaluate,AgentDistPath,PolicyPath,PricePath,ParamPath, Parameters, T, n_d, n_a, n_z, d_grid, a_grid,z_grid,simoptions)
+function AggVarsPath=EvalFnOnTransPath_AggVars_Case1(FnsToEvaluate,AgentDistPath,PolicyPath,PricePath,ParamPath, Parameters, T, n_d, n_a, n_z, d_grid, a_grid,z_grid,transpathoptions,simoptions)
 % AggVarsPath is T periods long (periods 0 (before the reforms are announced) & T are the initial and final values).
 
+%% Check which transpathoptions have been used, set all others to defaults 
+if exist('transpathoptions','var')==0
+    disp('No transpathoptions given, using defaults')
+    % If transpathoptions is not given, just use all the defaults
+    transpathoptions.verbose=0;
+else
+    % Check transpathoptions for missing fields, if there are some fill them with the defaults
+    if ~isfield(transpathoptions,'verbose')
+        transpathoptions.verbose=0;
+    end
+end
 if ~exist('simoptions','var')
     % If simoptions is not given, just use all the defaults
     simoptions.experienceasset=0;
diff --git a/EvaluateFnOnAgentDist/TransPathFHorz/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m b/EvaluateFnOnAgentDist/TransPathFHorz/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
index 73201f17..9f6ab038 100644
--- a/EvaluateFnOnAgentDist/TransPathFHorz/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
+++ b/EvaluateFnOnAgentDist/TransPathFHorz/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
@@ -184,8 +184,11 @@
 
     [FnsToEvaluate_temp,~, ~,~]=PType_FnsToEvaluate(FnsToEvaluate,Names_i,ii,l_d_temp,l_a_temp,l_z_temp,0);
 
-    AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1_FHorz(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
-
+    if isfinite(N_j_temp)
+        AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1_FHorz(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
+    else
+        AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
+    end
 
     FnNames_temp=fieldnames(FnsToEvaluate_temp);
     for ff=1:length(FnNames_temp)
diff --git a/HeterogeneousAgent/FHorz/HeteroAgentStationaryEqm_Case1_FHorz.m b/HeterogeneousAgent/FHorz/HeteroAgentStationaryEqm_Case1_FHorz.m
index 7d87b0e8..3561886d 100644
--- a/HeterogeneousAgent/FHorz/HeteroAgentStationaryEqm_Case1_FHorz.m
+++ b/HeterogeneousAgent/FHorz/HeteroAgentStationaryEqm_Case1_FHorz.m
@@ -133,6 +133,7 @@
     heteroagentoptions.useCustomModelStats=1;
     % Stash some of the inputs so they can be passed to CustomModelStats later (only things we otherwise overright).
     % So that user gets exactly what they input, not any internally reworked things
+    % In this (non-PType) context, these assignments are arrays or other simple types
     heteroagentoptions.CustomModelStatsInputs.z_grid=z_grid;
     heteroagentoptions.CustomModelStatsInputs.pi_z=pi_z;
     % Need the following two as otherwise they would contain alreadygridvals=1
diff --git a/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m b/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m
index 64d1e5be..a5787337 100644
--- a/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m
+++ b/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m
@@ -165,8 +165,9 @@
 heteroagentoptions.useCustomModelStats=0;
 if isfield(heteroagentoptions,'CustomModelStats')
     heteroagentoptions.useCustomModelStats=1;
-    % Stash some of the inputs so they can be passed to CustomModelStats later (only things we otherwise overright).
+    % Stash some of the inputs so they can be passed to CustomModelStats later (only things we otherwise overwrite).
     % So that user gets exactly what they input, not any internally reworked things
+    % In this (PType) context, these assignments are typically structs that each have all the PType fields within them
     heteroagentoptions.CustomModelStatsInputs.FnsToEvaluate=FnsToEvaluate;
     heteroagentoptions.CustomModelStatsInputs.n_d=n_d;
     heteroagentoptions.CustomModelStatsInputs.n_a=n_a;
@@ -431,6 +432,32 @@
         PTypeStructure.(iistr).pi_z=pi_z;
     end
 
+    %% Parameter Structure
+    % Parameters are allowed to be given as structure, or as vector/matrix
+    % (in terms of their dependence on permanent type). So go through each of
+    % these in term.
+    % ie. Parameters.alpha=[0;1]; or Parameters.alpha.ptype1=0; Parameters.alpha.ptype2=1;
+    % Need to establish this in PTypeStructure before we use it below.
+    PTypeStructure.(iistr).Parameters=Parameters;
+    FullParamNames=fieldnames(Parameters); % all the different parameters
+    nFields=length(FullParamNames);
+    for kField=1:nFields
+        if isa(Parameters.(FullParamNames{kField}), 'struct') % Check the current parameter for permanent type in structure form
+            % Check if this parameter is used for the current permanent type (it may or may not be, some parameters are only used be a subset of permanent types)
+            if isfield(Parameters.(FullParamNames{kField}),Names_i{ii})
+                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=Parameters.(FullParamNames{kField}).(Names_i{ii});
+            end
+        elseif sum(size(Parameters.(FullParamNames{kField}))==PTypeStructure.N_i)>=1 % Check for permanent type in vector/matrix form.
+            temp=Parameters.(FullParamNames{kField});
+            [~,ptypedim]=max(size(Parameters.(FullParamNames{kField}))==PTypeStructure.N_i); % Parameters as vector/matrix can be at most two dimensional, figure out which relates to PType, it should be the row dimension, if it is not then give a warning.
+            if ptypedim==1
+                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=temp(ii,:);
+            elseif ptypedim==2
+                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=temp(:,ii);
+            end
+        end
+    end
+
     % Check if using ExogShockFn or EiidShockFn, and if so, do these use a parameter that is being determined in general eqm
     heteroagentoptions.gridsinGE(ii)=0;
     if isfield(PTypeStructure.(iistr).vfoptions,'ExogShockFn')
@@ -493,34 +520,8 @@
     end
     PTypeStructure.(iistr).ReturnFnParamNames=ReturnFnParamNamesFn(PTypeStructure.(iistr).ReturnFn,PTypeStructure.(iistr).n_d,PTypeStructure.(iistr).n_a,PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).vfoptions,PTypeStructure.(iistr).Parameters);
     
-    %% Parameter Structure
-    % Parameters are allowed to be given as structure, or as vector/matrix
-    % (in terms of their dependence on permanent type). So go through each of
-    % these in term.
-    % ie. Parameters.alpha=[0;1]; or Parameters.alpha.ptype1=0; Parameters.alpha.ptype2=1;
-    PTypeStructure.(iistr).Parameters=Parameters;
-    FullParamNames=fieldnames(Parameters); % all the different parameters
-    nFields=length(FullParamNames);
-    for kField=1:nFields
-        if isa(Parameters.(FullParamNames{kField}), 'struct') % Check the current parameter for permanent type in structure form
-            % Check if this parameter is used for the current permanent type (it may or may not be, some parameters are only used be a subset of permanent types)
-            if isfield(Parameters.(FullParamNames{kField}),Names_i{ii})
-                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=Parameters.(FullParamNames{kField}).(Names_i{ii});
-            end
-        elseif sum(size(Parameters.(FullParamNames{kField}))==PTypeStructure.N_i)>=1 % Check for permanent type in vector/matrix form.
-            temp=Parameters.(FullParamNames{kField});
-            [~,ptypedim]=max(size(Parameters.(FullParamNames{kField}))==PTypeStructure.N_i); % Parameters as vector/matrix can be at most two dimensional, figure out which relates to PType, it should be the row dimension, if it is not then give a warning.
-            if ptypedim==1
-                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=temp(ii,:);
-            elseif ptypedim==2
-                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=temp(:,ii);
-            end
-        end
-    end
-
     %% jequaloneDist and AgeWeightsParamNames
     if isfinite(PTypeStructure.(iistr).N_j) % FHorz
-
         if isstruct(jequaloneDist)
             if isfield(jequaloneDist,PTypeStructure.Names_i{ii})
                 if isa(jequaloneDist, 'function_handle')
@@ -529,9 +530,7 @@
                     PTypeStructure.(iistr).jequaloneDist=jequaloneDist.(PTypeStructure.Names_i{ii});
                 end
             else
-                if isfinite(PTypeStructure.(iistr).N_j)
-                    error(['You must input jequaloneDist for permanent type ', PTypeStructure.Names_i{ii}, ' \n'])
-                end
+                error(['You must input jequaloneDist for permanent type ', PTypeStructure.Names_i{ii}, ' \n'])
             end
         else
             PTypeStructure.(iistr).jequaloneDist=jequaloneDist;
@@ -542,9 +541,7 @@
             if isfield(AgeWeightsParamNames,Names_i{ii})
                 PTypeStructure.(iistr).AgeWeightParamNames=AgeWeightsParamNames.(Names_i{ii});
             else
-                if isfinite(PTypeStructure.(iistr).N_j)
-                    error(['You must input AgeWeightParamNames for permanent type ', Names_i{ii}, ' \n'])
-                end
+                error(['You must input AgeWeightParamNames for permanent type ', Names_i{ii}, ' \n'])
             end
         end
     end
diff --git a/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType_subfn.m b/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType_subfn.m
index 41e20634..2b811f76 100644
--- a/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType_subfn.m
+++ b/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType_subfn.m
@@ -3,6 +3,31 @@
 heteroagentparamsvecindex=0:1:length(GEpricesvec);
 [GEpricesvec,penalty]=ParameterConstraints_TransformParamsToOriginal(GEpricesvec,heteroagentparamsvecindex,GEPriceParamNames,heteroagentoptions);
 
+if heteroagentoptions.verbose>0
+    GEpricesvec_tminus1=zeros(nGEprices,1);
+    AggVars_tminus1=NaN(length(AggVarNames),1);
+
+    for pp=1:nGEprices
+        GEpricesvec_tminus1(pp)=Parameters.(GEPriceParamNames{pp});
+    end
+    GEpricesvec_delta=GEpricesvec-GEpricesvec_tminus1; % Compute this to show max change per round
+    for aa=1:length(AggVarNames)
+        if isfield(Parameters,AggVarNames{aa})
+            AggVars_tminus1(aa)=Parameters.(AggVarNames{aa});
+        end
+    end
+    if heteroagentoptions.useintermediateEqns==1
+        intEqnnames=fieldnames(heteroagentoptions.intermediateEqns);
+        intermediateEqns_tminus1=zeros(length(intEqnnames),1);
+        for aa=1:length(intEqnnames)
+            if isfield(Parameters,intEqnnames{aa})
+                intEqns_tminus1(aa)=Parameters.(intEqnnames{aa});
+            end
+        end
+    end
+    % We don't do anything special for CustomModelStats, which are not as easily done as others above.
+end
+
 if heteroagentoptions.verbose==2
     fprintf(' \n')
     fprintf('Current GE prices: \n')
@@ -33,7 +58,7 @@
 
 
 %%
-AggVars_ConditionalOnPType=zeros(PTypeStructure.numFnsToEvaluate,PTypeStructure.N_i,'gpuArray'); % Create AggVars conditional on ptype.
+AggVars_ConditionalOnPType=zeros(PTypeStructure.numFnsToEvaluate,PTypeStructure.N_i); % Create AggVars conditional on ptype.
 
 for ii=1:PTypeStructure.N_i
     
@@ -82,12 +107,12 @@
         StationaryDist.(iistr)=StationaryDist_ii;
     end
 end
-AggVars=gather(sum(AggVars_ConditionalOnPType.*PTypeStructure.ptweights,2));
+AggVars=sum(AggVars_ConditionalOnPType.*PTypeStructure.ptweights',2);
 % Note: AggVars is a vector
 
 
 
-%% Put GE parameters  and AggVars in structure, so they can be used for intermediateEqns and GeneralEqmEqns
+%% Put GE parameters and AggVars in structure, so they can be used for intermediateEqns and GeneralEqmEqns
 % already did the basic GE params
 % for pp=1:nGEprices
 %     Parameters.(GEPriceParamNames{pp})=GEprices(pp);
@@ -172,8 +197,8 @@
     end
     if heteroagentoptions.useCustomModelStats==1
         fprintf('Current CustomModelStats variables: \n')
-        for ii=1:length(customstatnames)
-            fprintf(heteroagentoptions.verboseaccuracy1,customstatnames{ii},CustomStats.(customstatnames{ii}))
+        for aa=1:length(customstatnames)
+            fprintf(heteroagentoptions.verboseaccuracy1,customstatnames{aa},CustomStats.(customstatnames{aa}))
         end
     end
     fprintf('Current GeneralEqmEqns: \n')
diff --git a/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType.m b/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType.m
index c9ba76c6..17c5efbe 100644
--- a/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType.m
+++ b/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType.m
@@ -441,7 +441,7 @@
     FullParamNames=fieldnames(Parameters); % all the different parameters
     nFields=length(FullParamNames);
     for kField=1:nFields
-        if isa(Parameters.(FullParamNames{kField}), 'struct') % Check the current parameter for permanent type in structure form
+        if isstruct(Parameters.(FullParamNames{kField})) % Check the current parameter for permanent type in structure form            
             % Check if this parameter is used for the current permanent type (it may or may not be, some parameters are only used be a subset of permanent types)
             if isfield(Parameters.(FullParamNames{kField}),Names_i{ii})
                 PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=Parameters.(FullParamNames{kField}).(Names_i{ii});
@@ -473,6 +473,8 @@
     PTypeStructure.(iistr).FnsToEvaluate=FnsToEvaluate_temp;
     PTypeStructure.(iistr).FnsToEvaluateParamNames=FnsToEvaluateParamNames_temp;
     PTypeStructure.(iistr).WhichFnsForCurrentPType=WhichFnsForCurrentPType;
+    % Copied from FHorz/PType version for consistency
+    PTypeStructure.(iistr).FnsAndPTypeIndicator_ii=FnsAndPTypeIndicator_ii;
     
     %% PType masses
     if isa(PTypeDistParamNames, 'array')
diff --git a/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType_subfn.m b/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType_subfn.m
index 38e1e59f..f133e22d 100644
--- a/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType_subfn.m
+++ b/HeterogeneousAgent/InfHorz/PType/HeteroAgentStationaryEqm_Case1_PType_subfn.m
@@ -32,6 +32,7 @@
 
 
 %%
+AggVars_ConditionalOnPType=zeros(PTypeStructure.numFnsToEvaluate,PTypeStructure.N_i); % Create AggVars conditional on ptype.
 AggVars=zeros(PTypeStructure.numFnsToEvaluate,1,'gpuArray'); % numFnsToEvaluate is independent of the ptype
 
 for ii=1:PTypeStructure.N_i
@@ -55,11 +56,13 @@
     StationaryDist_ii=StationaryDist_Case1(Policy_ii,PTypeStructure.(iistr).n_d,PTypeStructure.(iistr).n_a,PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).pi_z,PTypeStructure.(iistr).simoptions,PTypeStructure.(iistr).Parameters);
     % PTypeStructure.(iistr).simoptions.outputasstructure=0; % Want AggVars_ii as matrix to make it easier to add them across the PTypes (is set outside this script)
     AggVars_ii=EvalFnOnAgentDist_AggVars_Case1(StationaryDist_ii, Policy_ii, PTypeStructure.(iistr).FnsToEvaluate, PTypeStructure.(iistr).Parameters, PTypeStructure.(iistr).FnsToEvaluateParamNames, PTypeStructure.(iistr).n_d, PTypeStructure.(iistr).n_a, PTypeStructure.(iistr).n_z, PTypeStructure.(iistr).d_grid, PTypeStructure.(iistr).a_grid, PTypeStructure.(iistr).z_grid, PTypeStructure.(iistr).simoptions);
-    
-    for kk=1:PTypeStructure.numFnsToEvaluate
-        jj=PTypeStructure.(iistr).WhichFnsForCurrentPType(kk);
-        if jj>0
-            AggVars(kk)=AggVars(kk)+PTypeStructure.(iistr).PTypeWeight*AggVars_ii(jj);
+    AggVars_ConditionalOnPType(PTypeStructure.(iistr).FnsAndPTypeIndicator_ii,ii)=AggVars_ii;
+    % Put updated AggVars into subsequent PTypeStructure Parameters, so they can be used for subsequent PType evaluations
+    FnsToEvaluate_aa=fieldnames(PTypeStructure.(iistr).FnsToEvaluate);
+    for jj=ii+1:PTypeStructure.N_i
+        jjstr=PTypeStructure.iistr{jj};
+        for aa=1:length(AggVars_ii)
+            PTypeStructure.(jjstr).Parameters.(FnsToEvaluate_aa{aa})=AggVars_ii(aa);
         end
     end
 
@@ -69,7 +72,8 @@
         StationaryDist.(iistr)=StationaryDist_ii;
     end
 end
-
+AggVars=sum(AggVars_ConditionalOnPType.*PTypeStructure.ptweights',2);
+% Note: AggVars is a vector
 
 
 %% Put GE parameters  and AggVars in structure, so they can be used for intermediateEqns and GeneralEqmEqns
diff --git a/Optimization/setupGEnewprice3_shooting.m b/Optimization/setupGEnewprice3_shooting.m
index b34b5fdc..6073e08a 100644
--- a/Optimization/setupGEnewprice3_shooting.m
+++ b/Optimization/setupGEnewprice3_shooting.m
@@ -91,15 +91,16 @@
     for gg=1:nGeneralEqmEqns
         if options.GEptype(gg)==1
             for gg2=1:size(options.GEnewprice3.howtoupdate,1)
-                if strcmp(options.GEnewprice3.howtoupdate{gg2,1},GEPriceParamNames{gg})
+                if strcmp(options.GEnewprice3.howtoupdate{gg2,1},GEeqnNames{gg})
                     pricename_gg=options.GEnewprice3.howtoupdate{gg2,2};
-                end
-            end
-            for pp=1:length(GEPriceParamNames)
-                if strcmp(pricename_gg,GEPriceParamNames{pp})
-                    if PricePathSizeVec(2,pp)-PricePathSizeVec(1,pp)+1~=N_i
-                        fprintf('Following error relates to GE condition %s and to price %s \n',GEeqnNames{gg},GEPriceParamNames{pp})
-                        error('You declared a GE condition to depend on permenent type, but the price that relates to it (in options.GEnewprice3.howtoupdate) does not depend on ptype')
+                    for pp=1:length(GEPriceParamNames)
+                        if strcmp(pricename_gg,GEPriceParamNames{pp})
+                            if PricePathSizeVec(2,pp)-PricePathSizeVec(1,pp)+1~=N_i
+                                fprintf('Following error relates to GE condition %s and to price %s \n',GEeqnNames{gg},GEPriceParamNames{pp})
+                                error('You declared a GE condition to depend on permenent type, but the price that relates to it (in options.GEnewprice3.howtoupdate) does not depend on ptype')
+                            end
+                            break
+                        end
                     end
                 end
             end
diff --git a/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m b/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m
index 0ac065fd..cd8501c7 100644
--- a/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m
+++ b/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m
@@ -1,4 +1,4 @@
-function AgentDistPath=AgentDistOnTransPath_Case1(AgentDist_initial, PolicyPath,n_d,n_a,n_z,pi_z,T,simoptions,Parameters,PricePath,ParamPath)
+function AgentDistPath=AgentDistOnTransPath_Case1(AgentDist_initial,PricePath,ParamPath,PolicyPath,n_d,n_a,n_z,pi_z,T,Parameters,simoptions)
 n_e=0; % NOT YET IMPLEMENTED FOR TRANSITION PATHS
 
 N_d=prod(n_d);
diff --git a/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m b/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m
index d1851e16..37708531 100644
--- a/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m
+++ b/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m
@@ -129,15 +129,17 @@
     else
         AgentDist_initial_temp=AgentDist_initial; % NEED TO DEAL WITH THIS PROPERLY
     end
-    if isstruct(AgeWeightsParamNames)
-        AgeWeightsParamNames_temp=AgeWeightsParamNames.(Names_i{ii});
-    else
-        AgeWeightsParamNames_temp=AgeWeightsParamNames;
-    end
-    if isstruct(jequalOneDist)
-        jequalOneDist_temp=jequalOneDist.(iistr);
-    else
-        jequalOneDist_temp=jequalOneDist;
+    if isfinite(N_j_temp)
+        if isstruct(AgeWeightsParamNames)
+            AgeWeightsParamNames_temp=AgeWeightsParamNames.(Names_i{ii});
+        else
+            AgeWeightsParamNames_temp=AgeWeightsParamNames;
+        end
+        if isstruct(jequalOneDist)
+            jequalOneDist_temp=jequalOneDist.(Names_i{ii});
+        else
+            jequalOneDist_temp=jequalOneDist;
+        end
     end
     
 
@@ -177,7 +179,11 @@
     
     
     % Compute the agent distribution path for permanent type ii
-    AgentDistPath_ii=AgentDistOnTransPath_Case1_FHorz(AgentDist_initial_temp, jequalOneDist_temp, PricePath_temp, ParamPath_temp, PolicyPath_temp, AgeWeightsParamNames_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,pi_z_temp, T,Parameters_temp, transpathoptions_temp, simoptions_temp);
+    if isfinite(N_j_temp)
+        AgentDistPath_ii=AgentDistOnTransPath_Case1_FHorz(AgentDist_initial_temp, jequalOneDist_temp, PricePath_temp, ParamPath_temp, PolicyPath_temp, AgeWeightsParamNames_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,pi_z_temp, T,Parameters_temp, transpathoptions_temp, simoptions_temp);
+    else
+        AgentDistPath_ii=AgentDistOnTransPath_Case1(AgentDist_initial_temp, PricePath_temp, ParamPath_temp, PolicyPath_temp, n_d_temp,n_a_temp,n_z_temp,pi_z_temp,T,Parameters_temp,simoptions_temp);
+    end
     % Note: T cannot depend on ptype, nor can PricePath depend on ptype
 
     AgentDistPath.(Names_i{ii})=AgentDistPath_ii;
diff --git a/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType.m b/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType.m
index 5670e520..daf58de4 100644
--- a/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType.m
+++ b/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType.m
@@ -1,4 +1,4 @@
-function varargout=TransitionPath_Case1_FHorz_PType(PricePath0, ParamPath, T, V_final, AgentDist_init, jequalOneDist, n_d, n_a, n_z, N_j, Names_i, d_grid,a_grid,z_grid, pi_z, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, AgeWeightsParamNames, PTypeDistParamNames, transpathoptions, simoptions, vfoptions)
+function varargout=TransitionPath_Case1_FHorz_PType(PricePath0, ParamPath, T, V_final, AgentDist_initial, jequalOneDist, n_d, n_a, n_z, N_j, Names_i, d_grid,a_grid,z_grid, pi_z, ReturnFn, FnsToEvaluate, GeneralEqmEqns, Parameters, DiscountFactorParamNames, AgeWeightsParamNames, PTypeDistParamNames, transpathoptions, simoptions, vfoptions)
 % This code will work for all transition paths except those that involve at
 % change in the transition matrix pi_z (can handle a change in pi_z, but
 % only if it is a 'surprise', not anticipated changes) 
@@ -138,15 +138,60 @@
 
 %% Some internal commands require a few vfoptions and simoptions to be set
 vfoptions.EVpre=0; % Not actually an option that can be used here
-if isfield(vfoptions,'lowmemory')==0
-    vfoptions.lowmemory=0;
+if ~isfield(vfoptions,'verbose')
+    vfoptions.verbose=0;
+end
+if ~isfield(vfoptions,'experienceasset')
+    for ii=1:N_i
+        vfoptions.experienceasset.(Names_i{ii})=0;
+    end
+else
+    % User created vfoptions.PType.experienceasset; we have Names_i
+    for ii=1:N_i
+        if isfield(vfoptions.experienceasset, Names_i{ii})
+            if ~isfield(vfoptions, 'aprimeFn') || ~isfield(vfoptions.aprimeFn, Names_i{ii})
+                error('To use an experience asset you must define vfoptions.aprimeFn')
+            end
+        else
+            vfoptions.experienceasset.(Names_i{ii})=0;
+        end
+    end
+end
+if ~isfield(vfoptions,'lowmemory')
+    for ii=1:N_i
+        vfoptions.lowmemory.(Names_i{ii})=0;
+    end
+else
+    % User created vfoptions.PType.lowmemory; we have Names_i
+    for ii=1:N_i
+        if ~isfield(vfoptions.lowmemory, Names_i{ii})
+            vfoptions.lowmemory.(Names_i{ii})=0;
+        end
+    end
 end
 
-%% Get AgeWeights from Parameters   
-try
-    AgeWeights=Parameters.(AgeWeightsParamNames{1});
-catch
-    error(['Failed to find parameter ', AgeWeightsParamNames{1}])
+%% Get AgeWeights from Parameters
+if isstruct(AgeWeightsParamNames)
+    for ii=1:N_i
+        try
+            AgeWeights=Parameters.(AgeWeightsParamNames.(Names_i{ii}){1});
+            if length(AgeWeights)~=N_j.(Names_i{ii})
+                error('Ageweights does not have age-like length')
+            else
+                break
+            end
+        catch
+            if ii==N_i
+                error(['Failed to find parameter ', AgeWeightsParamNames.(Names_i{ii}){1}])
+            end
+        end
+    end
+else
+    try
+        AgeWeights=Parameters.(AgeWeightsParamNames{1});
+    catch
+        error(['Failed to find parameter ', AgeWeightsParamNames{1}])
+    end
 end
 % Later, when creating PTypeStructure, we get the ptype-specific versions out of this and create an AgeWeights_T structure.
 
@@ -166,11 +211,11 @@
 if isstruct(GeneralEqmEqns)
     if length(PricePathNames)~=length(fieldnames(GeneralEqmEqns))
         fprintf('length(PricePathNames)=%i and length(fieldnames(GeneralEqmEqns))=%i (relates to following error) \n', length(PricePathNames), length(fieldnames(GeneralEqmEqns)))
-        error('Initial PricePath contains less variables than GeneralEqmEqns (structure) \n')
+        error('Initial PricePath contains fewer variables than GeneralEqmEqns (structure) \n')
     end
 else
     if length(PricePathNames)~=length(GeneralEqmEqns)
-        error('Initial PricePath contains less variables than GeneralEqmEqns')
+        error('Initial PricePath contains fewer variables than GeneralEqmEqns')
     end
 end
 
@@ -232,7 +277,7 @@
         PTypeStructure.(iistr).vfoptions.divideandconquer=0; %default
     else
         if PTypeStructure.(iistr).vfoptions.divideandconquer==1
-            PTypeStructure.(iistr).vfoptions.level1n=ceil(n_a/50); % default
+            PTypeStructure.(iistr).vfoptions.level1n=ceil(n_a.(iistr)/50); % default
         end
     end
     if ~isfield(PTypeStructure.(iistr).vfoptions,'gridinterplayer')
@@ -263,7 +308,7 @@
     
     % Horizon is determined via N_j
     if isstruct(N_j)
-        PTypeStructure.(iistr).N_j=N_j.(Names_i{ii});
+        PTypeStructure.(iistr).N_j=N_j.(iistr);
     elseif isscalar(N_j)
         PTypeStructure.(iistr).N_j=N_j;
     else
@@ -271,7 +316,7 @@
     end
     
     if isa(n_d,'struct')
-        PTypeStructure.(iistr).n_d=n_d.(Names_i{ii});
+        PTypeStructure.(iistr).n_d=n_d.(iistr);
     else
         PTypeStructure.(iistr).n_d=n_d;
     end
@@ -280,37 +325,41 @@
     if N_d==0
         PTypeStructure.(iistr).l_d=0;
     else
-        PTypeStructure.(iistr).l_d=length(n_d);
+        PTypeStructure.(iistr).l_d=length(PTypeStructure.(iistr).n_d);
     end
     if isa(n_a,'struct')
-        PTypeStructure.(iistr).n_a=n_a.(Names_i{ii});
+        PTypeStructure.(iistr).n_a=n_a.(iistr);
     else
         PTypeStructure.(iistr).n_a=n_a;
     end
     N_a=prod(PTypeStructure.(iistr).n_a);
     PTypeStructure.(iistr).N_a=N_a;
     if isa(n_z,'struct')
-        PTypeStructure.(iistr).n_z=n_z.(Names_i{ii});
+        PTypeStructure.(iistr).n_z=n_z.(iistr);
     else
         PTypeStructure.(iistr).n_z=n_z;
     end
     N_z=prod(PTypeStructure.(iistr).n_z);
     PTypeStructure.(iistr).N_z=N_z;
+    if N_z==0 && ~isfinite(PTypeStructure.(iistr).N_j)
+        % This simplifies InfoHorz conditional logic for the purposes of this function
+        N_z=1;
+    end
     N_e=prod(PTypeStructure.(iistr).n_e);
     PTypeStructure.(iistr).N_e=N_e;
 
     if isa(d_grid,'struct')
-        PTypeStructure.(iistr).d_grid=gpuArray(d_grid.(Names_i{ii}));
+        PTypeStructure.(iistr).d_grid=gpuArray(d_grid.(iistr));
     else
         PTypeStructure.(iistr).d_grid=gpuArray(d_grid);
     end
     if isa(a_grid,'struct')
-        PTypeStructure.(iistr).a_grid=gpuArray(a_grid.(Names_i{ii}));
+        PTypeStructure.(iistr).a_grid=gpuArray(a_grid.(iistr));
     else
         PTypeStructure.(iistr).a_grid=gpuArray(a_grid);
     end
     if isa(z_grid,'struct')
-        PTypeStructure.(iistr).z_grid=gpuArray(z_grid.(Names_i{ii}));
+        PTypeStructure.(iistr).z_grid=gpuArray(z_grid.(iistr));
     else
         PTypeStructure.(iistr).z_grid=gpuArray(z_grid);
     end
@@ -341,14 +390,14 @@
     % end
 
     if isa(pi_z,'struct')
-        PTypeStructure.(iistr).pi_z=pi_z.(Names_i{ii}); % Different grids by permanent type, but not depending on age. (same as the case just above; this case can occour with or without the existence of vfoptions, as long as there is no vfoptions.agedependentgrids)
+        PTypeStructure.(iistr).pi_z=pi_z.(iistr); % Different grids by permanent type, but not depending on age. (same as the case just above; this case can occour with or without the existence of vfoptions, as long as there is no vfoptions.agedependentgrids)
     else
         PTypeStructure.(iistr).pi_z=pi_z;
     end
     
     PTypeStructure.(iistr).ReturnFn=ReturnFn;
     if isa(ReturnFn,'struct')
-        PTypeStructure.(iistr).ReturnFn=ReturnFn.(Names_i{ii});
+        PTypeStructure.(iistr).ReturnFn=ReturnFn.(iistr);
     end
     
     % Parameters are allowed to be given as structure, or as vector/matrix (in terms of their dependence on permanent type). 
@@ -360,8 +409,8 @@
     for kField=1:nFields
         if isa(Parameters.(FullParamNames{kField}), 'struct') % Check the current parameter for permanent type in structure form
             % Check if this parameter is used for the current permanent type (it may or may not be, some parameters are only used be a subset of permanent types)
-            if isfield(Parameters.(FullParamNames{kField}),Names_i{ii})
-                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=Parameters.(FullParamNames{kField}).(Names_i{ii});
+            if isfield(Parameters.(FullParamNames{kField}),iistr)
+                PTypeStructure.(iistr).Parameters.(FullParamNames{kField})=Parameters.(FullParamNames{kField}).(iistr);
             end
         elseif sum(size(Parameters.(FullParamNames{kField}))==PTypeStructure.N_i)>=1 % Check for permanent type in vector/matrix form.
             temp=Parameters.(FullParamNames{kField});
@@ -385,7 +434,7 @@
     % The parameter names can be made to depend on the permanent-type
     PTypeStructure.(iistr).DiscountFactorParamNames=DiscountFactorParamNames;
     if isa(DiscountFactorParamNames,'struct')
-        PTypeStructure.(iistr).DiscountFactorParamNames=DiscountFactorParamNames.(Names_i{ii});
+        PTypeStructure.(iistr).DiscountFactorParamNames=DiscountFactorParamNames.(iistr);
     end
 
     % Implement new way of handling ReturnFn inputs (note l_d, l_a, l_z are just created for this and then not used for anything else later)
@@ -408,13 +457,8 @@
             l_e=length(PTypeStructure.(iistr).vfoptions.n_e);
         end
     end
-    % Figure out ReturnFnParamNames from ReturnFn
-    temp=getAnonymousFnInputNames(PTypeStructure.(iistr).ReturnFn);
-    if length(temp)>(l_d+l_a+l_a+l_z+l_e) % This is largely pointless, the ReturnFn is always going to have some parameters
-        ReturnFnParamNames={temp{l_d+l_a+l_a+l_z+l_e+1:end}}; % the first inputs will always be (d,aprime,a,z)
-    else
-        ReturnFnParamNames={};
-    end
+    %% Implement new way of handling ReturnFn inputs
+    ReturnFnParamNames=ReturnFnParamNamesFn(PTypeStructure.(iistr).ReturnFn,PTypeStructure.(iistr).n_d,PTypeStructure.(iistr).n_a,PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).vfoptions,Parameters);
     PTypeStructure.(iistr).ReturnFnParamNames=ReturnFnParamNames;
     
 
@@ -422,17 +466,18 @@
     % Only the relevant ones need to be evaluated.
     % The dependence of FnsToEvaluateFn and FnsToEvaluateFnParamNames are necessarily the same.
     PTypeStructure.(iistr).FnsToEvaluate={};
+    PTypeStructure.(iistr).FnsToEvaluateParamNames={};
 
     FnNames=fieldnames(FnsToEvaluate);
     PTypeStructure.numFnsToEvaluate=length(fieldnames(FnsToEvaluate));
     PTypeStructure.(iistr).WhichFnsForCurrentPType=zeros(PTypeStructure.numFnsToEvaluate,1);
     jj=1; % jj indexes the FnsToEvaluate that are relevant to the current PType
     for kk=1:PTypeStructure.numFnsToEvaluate
-        if isa(FnsToEvaluate.(FnNames{kk}),'struct')
-            if isfield(FnsToEvaluate.(FnNames{kk}), Names_i{ii})
-                PTypeStructure.(iistr).FnsToEvaluate{jj}=FnsToEvaluate.(FnNames{kk}).(Names_i{ii});
+        if isstruct(FnsToEvaluate.(FnNames{kk}))
+            if isfield(FnsToEvaluate.(FnNames{kk}), iistr)
+                PTypeStructure.(iistr).FnsToEvaluate{jj}=FnsToEvaluate.(FnNames{kk}).(iistr);
                 % Figure out FnsToEvaluateParamNames
-                temp=getAnonymousFnInputNames(FnsToEvaluate.(FnNames{kk}).(Names_i{ii}));
+                temp=getAnonymousFnInputNames(FnsToEvaluate.(FnNames{kk}).(iistr));
                 PTypeStructure.(iistr).FnsToEvaluateParamNames(jj).Names={temp{l_d+l_a+l_a+l_z+l_e+1:end}}; % the first inputs will always be (d,aprime,a,z)
                 PTypeStructure.(iistr).WhichFnsForCurrentPType(kk)=jj; jj=jj+1;
                 PTypeStructure.FnsAndPTypeIndicator(kk,ii)=1;
@@ -450,17 +495,20 @@
             PTypeStructure.FnsAndPTypeIndicator(kk,ii)=1;
         end
     end
-    PTypeStructure.(iistr).AggVarNames=FnNames;
+    PTypeStructure.(iistr).AggVarNames=FnNames(logical(PTypeStructure.FnsAndPTypeIndicator(:,ii)));
 
 
     %% Set up exogenous shock processes
-    [PTypeStructure.(iistr).z_gridvals_J, PTypeStructure.(iistr).pi_z_J, PTypeStructure.(iistr).pi_z_J_sim, PTypeStructure.(iistr).e_gridvals_J, PTypeStructure.(iistr).pi_e_J, PTypeStructure.(iistr).pi_e_J_sim, PTypeStructure.(iistr).ze_gridvals_J_fastOLG, transpathoptions, simoptions]=ExogShockSetup_TPath_FHorz(PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).z_grid,PTypeStructure.(iistr).pi_z,PTypeStructure.(iistr).N_a,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).Parameters,PricePathNames,ParamPathNames,transpathoptions,PTypeStructure.(iistr).simoptions,4);
-    % Convert z and e to age-dependent joint-grids and transtion matrix
-    % output: z_gridvals_J, pi_z_J, e_gridvals_J, pi_e_J, transpathoptions,vfoptions,simoptions
-    
+    if isfinite(PTypeStructure.(iistr).N_j)
+        [PTypeStructure.(iistr).z_gridvals_J, PTypeStructure.(iistr).pi_z_J, PTypeStructure.(iistr).pi_z_J_sim, PTypeStructure.(iistr).e_gridvals_J, PTypeStructure.(iistr).pi_e_J, PTypeStructure.(iistr).pi_e_J_sim, PTypeStructure.(iistr).ze_gridvals_J_fastOLG, transpathoptions, simoptions]=ExogShockSetup_TPath_FHorz(PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).z_grid,PTypeStructure.(iistr).pi_z,PTypeStructure.(iistr).N_a,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).Parameters,PricePathNames,ParamPathNames,transpathoptions,PTypeStructure.(iistr).simoptions,4);
+        % Convert z and e to age-dependent joint-grids and transtion matrix
+        % output: z_gridvals_J, pi_z_J, e_gridvals_J, pi_e_J, transpathoptions,vfoptions,simoptions
+    end
+
     %% We can precompute some things that the fastOLG needs for the agent dist
     % But only bother with this when using fastOLG=1
-    if transpathoptions.fastOLG==1
+    if transpathoptions.fastOLG==1 && isfinite(PTypeStructure.(iistr).N_j)
+        N_j_temp=PTypeStructure.(iistr).N_j;
         if PTypeStructure.(iistr).simoptions.gridinterplayer==1
             N_probs=2;
         else
@@ -469,57 +517,57 @@
         if N_z==0
             if N_e==0 % no z, no e
                 if PTypeStructure.(iistr).simoptions.gridinterplayer==0
-                    II1=1:1:N_a*(N_j-1);
-                    PTypeStructure.(iistr).II2=ones(N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',N_j-1,1)+repelem(N_a*(0:1:N_j-2)',N_a,1); % Note: there is one use of N_j which is because we want to index AgentDist
+                    II1=1:1:N_a*(N_j_temp-1);
+                    PTypeStructure.(iistr).II2=ones(N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',N_j_temp-1,1)+repelem(N_a*(0:1:N_j_temp-2)',N_a,1); % Note: there is one use of N_j which is because we want to index AgentDist
                     PTypeStructure.(iistr).exceptfirstj=[]; % not needed
                     PTypeStructure.(iistr).justfirstj=[]; % not needed
                 elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
-                    II=repelem((1:1:N_a*(N_j-1))',1,N_probs);
+                    II=repelem((1:1:N_a*(N_j_temp-1))',1,N_probs);
                     PTypeStructure.(iistr).exceptlastj=[]; % not needed
                     PTypeStructure.(iistr).exceptfirstj=[]; % not needed
                     PTypeStructure.(iistr).justfirstj=[]; % not needed
                 end
             else % no z, yes e
                 if PTypeStructure.(iistr).simoptions.gridinterplayer==0
-                    II1=1:1:N_a*(N_j-1)*N_e;
-                    PTypeStructure.(iistr).II2=ones(N_a*(N_j-1)*N_e,1);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j-1)*N_e,1)+repmat(repelem(N_a*(0:1:N_j-2)',N_a,1),N_e,1)+repelem(N_a*N_j*(0:1:N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j-1)*N_e,1)+repmat(repelem(N_a*(1:1:N_j-1)',N_a,1),N_e,1)+repelem(N_a*N_j*(0:1:N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_e,1)+N_a*N_j*repelem((0:1:N_e-1)',N_a,1);
+                    II1=1:1:N_a*(N_j_temp-1)*N_e;
+                    PTypeStructure.(iistr).II2=ones(N_a*(N_j_temp-1)*N_e,1);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j_temp-1)*N_e,1)+repmat(repelem(N_a*(0:1:N_j_temp-2)',N_a,1),N_e,1)+repelem(N_a*N_j_temp*(0:1:N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j_temp-1)*N_e,1)+repmat(repelem(N_a*(1:1:N_j_temp-1)',N_a,1),N_e,1)+repelem(N_a*N_j_temp*(0:1:N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_e,1)+N_a*N_j_temp*repelem((0:1:N_e-1)',N_a,1);
                 elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
-                    II=repelem((1:1:N_a*(N_j-1)*N_e)',1,N_probs);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j-1)*N_e,1)+repmat(repelem(N_a*(0:1:N_j-2)',N_a,1),N_e,1)+repelem(N_a*N_j*(0:1:N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j-1)*N_e,1)+repmat(repelem(N_a*(1:1:N_j-1)',N_a,1),N_e,1)+repelem(N_a*N_j*(0:1:N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_e,1)+N_a*N_j*repelem((0:1:N_e-1)',N_a,1);
+                    II=repelem((1:1:N_a*(N_j_temp-1)*N_e)',1,N_probs);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j_temp-1)*N_e,1)+repmat(repelem(N_a*(0:1:N_j_temp-2)',N_a,1),N_e,1)+repelem(N_a*N_j_temp*(0:1:N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j_temp-1)*N_e,1)+repmat(repelem(N_a*(1:1:N_j_temp-1)',N_a,1),N_e,1)+repelem(N_a*N_j_temp*(0:1:N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_e,1)+N_a*N_j_temp*repelem((0:1:N_e-1)',N_a,1);
                 end
             end
         else % N_z>0
             if N_e==0 % z, no e
                 if PTypeStructure.(iistr).simoptions.gridinterplayer==0
-                    II1=1:1:N_a*(N_j-1)*N_z;
-                    PTypeStructure.(iistr).II2=ones(N_a*(N_j-1)*N_z,1);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j-1)*N_z,1)+repmat(repelem(N_a*(0:1:N_j-2)',N_a,1),N_z,1)+repelem(N_a*N_j*(0:1:N_z-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j-1)*N_z,1)+repmat(repelem(N_a*(1:1:N_j-1)',N_a,1),N_z,1)+repelem(N_a*N_j*(0:1:N_z-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z,1)+N_a*N_j*repelem((0:1:N_z-1)',N_a,1);
+                    II1=1:1:N_a*(N_j_temp-1)*N_z;
+                    PTypeStructure.(iistr).II2=ones(N_a*(N_j_temp-1)*N_z,1);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j_temp-1)*N_z,1)+repmat(repelem(N_a*(0:1:N_j_temp-2)',N_a,1),N_z,1)+repelem(N_a*N_j_temp*(0:1:N_z-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j_temp-1)*N_z,1)+repmat(repelem(N_a*(1:1:N_j_temp-1)',N_a,1),N_z,1)+repelem(N_a*N_j_temp*(0:1:N_z-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z,1)+N_a*N_j_temp*repelem((0:1:N_z-1)',N_a,1);
                 elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
-                    II=repelem((1:1:N_a*(N_j-1)*N_z)',1,N_probs);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j-1)*N_z,1)+repmat(repelem(N_a*(0:1:N_j-2)',N_a,1),N_z,1)+repelem(N_a*N_j*(0:1:N_z-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j-1)*N_z,1)+repmat(repelem(N_a*(1:1:N_j-1)',N_a,1),N_z,1)+repelem(N_a*N_j*(0:1:N_z-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z,1)+N_a*N_j*repelem((0:1:N_z-1)',N_a,1);
+                    II=repelem((1:1:N_a*(N_j_temp-1)*N_z)',1,N_probs);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j_temp-1)*N_z,1)+repmat(repelem(N_a*(0:1:N_j_temp-2)',N_a,1),N_z,1)+repelem(N_a*N_j_temp*(0:1:N_z-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j_temp-1)*N_z,1)+repmat(repelem(N_a*(1:1:N_j_temp-1)',N_a,1),N_z,1)+repelem(N_a*N_j_temp*(0:1:N_z-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z,1)+N_a*N_j_temp*repelem((0:1:N_z-1)',N_a,1);
                 end
             else % z and e
                 if PTypeStructure.(iistr).simoptions.gridinterplayer==0
-                    II1=1:1:N_a*(N_j-1)*N_z*N_e;
-                    PTypeStructure.(iistr).II2=ones(N_a*(N_j-1)*N_z*N_e,1);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j-1)*N_z*N_e,1)+repmat(repelem(N_a*(0:1:N_j-2)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j*(0:1:N_z*N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j-1)*N_z*N_e,1)+repmat(repelem(N_a*(1:1:N_j-1)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j*(0:1:N_z*N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z*N_e,1)+N_a*N_j*repelem((0:1:N_z*N_e-1)',N_a,1);
+                    II1=1:1:N_a*(N_j_temp-1)*N_z*N_e;
+                    PTypeStructure.(iistr).II2=ones(N_a*(N_j_temp-1)*N_z*N_e,1);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j_temp-1)*N_z*N_e,1)+repmat(repelem(N_a*(0:1:N_j_temp-2)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j_temp*(0:1:N_z*N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j_temp-1)*N_z*N_e,1)+repmat(repelem(N_a*(1:1:N_j_temp-1)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j_temp*(0:1:N_z*N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z*N_e,1)+N_a*N_j_temp*repelem((0:1:N_z*N_e-1)',N_a,1);
                 elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
-                    II=repelem((1:1:N_a*(N_j-1)*N_z*N_e)',1,N_probs);
-                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j-1)*N_z*N_e,1)+repmat(repelem(N_a*(0:1:N_j-2)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j*(0:1:N_z*N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j-1)*N_z*N_e,1)+repmat(repelem(N_a*(1:1:N_j-1)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j*(0:1:N_z*N_e-1)',N_a*(N_j-1),1);
-                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z*N_e,1)+N_a*N_j*repelem((0:1:N_z*N_e-1)',N_a,1);
+                    II=repelem((1:1:N_a*(N_j_temp-1)*N_z*N_e)',1,N_probs);
+                    PTypeStructure.(iistr).exceptlastj=repmat((1:1:N_a)',(N_j_temp-1)*N_z*N_e,1)+repmat(repelem(N_a*(0:1:N_j_temp-2)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j_temp*(0:1:N_z*N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).exceptfirstj=repmat((1:1:N_a)',(N_j_temp-1)*N_z*N_e,1)+repmat(repelem(N_a*(1:1:N_j_temp-1)',N_a,1),N_z*N_e,1)+repelem(N_a*N_j_temp*(0:1:N_z*N_e-1)',N_a*(N_j_temp-1),1);
+                    PTypeStructure.(iistr).justfirstj=repmat((1:1:N_a)',N_z*N_e,1)+N_a*N_j_temp*repelem((0:1:N_z*N_e-1)',N_a,1);
                 end
             end
         end
@@ -530,160 +578,250 @@
             PTypeStructure.(iistr).II1orII=II;
             PTypeStructure.(iistr).II2=[];
         end
-    else
-        
+    elseif transpathoptions.fastOLG==1 % && ~isfinite(PTypeStructure.(iistr).N_j)
+        if PTypeStructure.(iistr).simoptions.gridinterplayer==1
+            N_probs=2;
+        else
+            N_probs=1;
+        end
+        if N_z==0 % Never N_e in InfHorz
+            if PTypeStructure.(iistr).simoptions.gridinterplayer==0
+                II1=1:1:N_a;
+                PTypeStructure.(iistr).II2=ones(N_a,1);
+            elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
+                II=repelem((1:1:N_a)',1,N_probs);
+            end
+        else % N_z>0, no N_e to worry about
+            if PTypeStructure.(iistr).simoptions.gridinterplayer==0
+                II1=1:1:N_a*N_z;
+                PTypeStructure.(iistr).II2=ones(N_a*N_z,1);
+            elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
+                II=repelem((1:1:N_a*N_z)',1,N_probs);
+            end
+        end
+        PTypeStructure.(iistr).exceptlastj=[]; % not needed
+        PTypeStructure.(iistr).exceptfirstj=[]; % not needed
+        PTypeStructure.(iistr).justfirstj=[]; % not needed
+        % To keep inputs simpler
+        if PTypeStructure.(iistr).simoptions.gridinterplayer==0
+            PTypeStructure.(iistr).II1orII=II1;
+        elseif PTypeStructure.(iistr).simoptions.gridinterplayer==1
+            PTypeStructure.(iistr).II1orII=II;
+            PTypeStructure.(iistr).II2=[];
+        end        
     end
 
     
-    %% Organise V_final and AgentDist_initial
+    %% Organise V_final, AgeWeights and AgentDist_initial
     % Reshape V_final
-    if transpathoptions.fastOLG==0
+    if ~isfinite(PTypeStructure.(iistr).N_j)
+        % If no z, then N_z=1 here
+        V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_z]);
+    elseif transpathoptions.fastOLG==0
+        N_j_temp=PTypeStructure.(iistr).N_j;
         if N_z==0
             if N_e==0
-                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_j]);
+                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_j_temp]);
             else
-                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_e,N_j]);
+                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_e,N_j_temp]);
             end
         else
             if N_e==0
-                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_z,N_j]);
+                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_z,N_j_temp]);
             else
-                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_z,N_e,N_j]);
+                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_z,N_e,N_j_temp]);
             end
         end
-    else
+    else % transpathoptions.fastOLG==1
+        N_j_temp=PTypeStructure.(iistr).N_j;
         if N_z==0
             if N_e==0
-                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_j]);
+                V_final.(iistr)=reshape(V_final.(iistr),[N_a,N_j_temp]);
             else
-                V_final.(iistr)=reshape(permute(V_final.(iistr),[1,3,2]),[N_a*N_j,N_e]);
+                V_final.(iistr)=reshape(permute(V_final.(iistr),[1,3,2]),[N_a*N_j_temp,N_e]);
             end
         else
             if N_e==0
-                V_final.(iistr)=reshape(permute(V_final.(iistr),[1,3,2]),[N_a*N_j,N_z]);
+                V_final.(iistr)=reshape(permute(V_final.(iistr),[1,3,2]),[N_a*N_j_temp,N_z]);
             else
-                V_final.(iistr)=reshape(permute(V_final.(iistr),[1,4,2,3]),[N_a*N_j,N_z,N_e]);
+                V_final.(iistr)=reshape(permute(V_final.(iistr),[1,4,2,3]),[N_a*N_j_temp,N_z,N_e]);
             end
         end
     end
-    % Reshape AgentDist_initial
-    if N_z==0
-        if N_e==0
-            AgentDist_init.(iistr)=reshape(AgentDist_init.(iistr),[N_a,N_j]); % if simoptions.fastOLG==0
-            AgeWeights_init.(iistr)=sum(AgentDist_init.(iistr),1); % [1,N_j]
-            if PTypeStructure.(iistr).simoptions.fastOLG==1
-                AgentDist_init.(iistr)=reshape(AgentDist_init.(iistr),[N_a*N_j,1]);
-                AgeWeights_init.(iistr)=repelem(AgeWeights_init.(iistr)',N_a,1);
+    % Reshape AgentDist_initial and turn AgeWeights_T into appropriate size so that we can always just do AgentDist.*AgeWeights
+    % Note when simoptions.fastOLG==1 we have shapes of [N_a*N_j_temp*whatever,1-or-N_e] intead of [N_a,whatever-and-maybe-N_e,N_j_temp]
+    AgentDist_init=AgentDist_initial.(iistr);
+    if isfinite(PTypeStructure.(iistr).N_j)
+        N_j_temp=PTypeStructure.(iistr).N_j;
+        if N_z==0
+            if N_e==0
+                AgentDist_init=reshape(AgentDist_init,[N_a,N_j_temp]);
+                AgeWeights_init=sum(AgentDist_init,1); % [1,N_j]
+                if simoptions.fastOLG
+                    AgentDist_init=reshape(AgentDist_init,[N_a*N_j_temp,1]);
+                end
+            else
+                AgentDist_init=reshape(AgentDist_init,[N_a*N_e,N_j_temp]);
+                AgeWeights_init=sum(AgentDist_init,1); % [1,N_j]
+                if simoptions.fastOLG
+                    AgentDist_init=reshape(permute(reshape(AgentDist_init,[N_a,N_e,N_j_temp]),[1,3,2]),[N_a*N_j_temp,N_e]);
+                end
             end
         else
-            AgentDist_init.(iistr)=reshape(AgentDist_init.(iistr),[N_a*N_e,N_j]); % if simoptions.fastOLG==0
-            AgeWeights_init.(iistr)=sum(AgentDist_init.(iistr),1); % [1,N_j]
-            if PTypeStructure.(iistr).simoptions.fastOLG==1 % simoptions.fastOLG==1, so AgentDist is treated as : (a,j,z)-by-1
-                AgentDist_init.(iistr)=reshape(permute(reshape(AgentDist_init.(iistr),[N_a,N_e,N_j]),[1,3,2]),[N_a*N_j*N_e,1]);
-                AgeWeights_init.(iistr)=repelem(AgeWeights_init.(iistr)',N_a,1);
+            if N_e==0
+                AgentDist_init=reshape(AgentDist_init,[N_a*N_z,N_j_temp]);
+                AgeWeights_init=sum(AgentDist_init,1); % [1,N_j]
+                if simoptions.fastOLG
+                    AgentDist_init=reshape(permute(reshape(AgentDist_init,[N_a,N_z,N_j_temp]),[1,3,2]),[N_a*N_j_temp*N_z,1]);
+                end
+            else
+                AgentDist_init=reshape(AgentDist_init,[N_a*N_z*N_e,N_j_temp]);
+                AgeWeights_init=sum(AgentDist_init,1); % [1,N_j]
+                if simoptions.fastOLG
+                    AgentDist_init=reshape(permute(reshape(AgentDist_init,[N_a,N_z,N_e,N_j_temp]),[1,4,2,3]),[N_a*N_j_temp*N_z,N_e]);
+                end
             end
         end
-    else
-        if N_e==0
-            AgentDist_init.(iistr)=reshape(AgentDist_init.(iistr),[N_a*N_z,N_j]); % if simoptions.fastOLG==0
-            AgeWeights_init.(iistr)=sum(AgentDist_init.(iistr),1); % [1,N_j]
-            if PTypeStructure.(iistr).simoptions.fastOLG==1 % simoptions.fastOLG==1, so AgentDist is treated as : (a,j,z)-by-1
-                AgentDist_init.(iistr)=reshape(permute(reshape(AgentDist_init.(iistr),[N_a,N_z,N_j]),[1,3,2]),[N_a*N_j*N_z,1]);
-                AgeWeights_init.(iistr)=repelem(AgeWeights_init.(iistr)',N_a,1);
-            end
+
+        % Get AgeWeights and switch into the transpathoptions.ageweightstrivial=0 setup (and this is what subfns hardcode when doing PTypes)
+        % It is assumed there is only one Age Weight Parameter (name))
+        % AgeWeights_T is (a,j,z)-by-T (create as N_j-by-T to start, then switch)
+        if isstruct(AgeWeights)
+            AgeWeights_ii=AgeWeights.(iistr);
         else
-            AgentDist_init.(iistr)=reshape(AgentDist_init.(iistr),[N_a*N_z*N_e,N_j]); % if simoptions.fastOLG==0
-            AgeWeights_init.(iistr)=sum(AgentDist_init.(iistr),1); % [1,N_j]
-            if PTypeStructure.(iistr).simoptions.fastOLG==1 % simoptions.fastOLG==1, so AgentDist is treated as : (a,j,z)-by-1
-                AgentDist_init.(iistr)=reshape(permute(reshape(AgentDist_init.(iistr),[N_a,N_z,N_e,N_j]),[1,4,2,3]),[N_a*N_j*N_z,N_e]);
-                AgeWeights_init.(iistr)=repelem(AgeWeights_init.(iistr)',N_a,1);
-            end
+            % not a structure, so must apply to all permanent types
+            AgeWeights_ii=AgeWeights;
         end
-    end
-    % Get AgeWeights and switch into the transpathoptions.ageweightstrivial=0 setup (and this is what subfns hardcode when doing PTypes)
-    % It is assumed there is only one Age Weight Parameter (name))
-    % AgeWeights_T is (a,j,z)-by-T (create as j-by-T to start, then switch)
-    if isstruct(AgeWeights)
-        AgeWeights_ii=AgeWeights.(iistr);
-        if all(size(AgeWeights_ii)==[N_j,1])
+        if all(size(AgeWeights_ii)==[N_j_temp,1])
             % Does not depend on transition path period
-            AgeWeights_T.(iistr)=gather(AgeWeights_ii.*ones(1,T));
-        elseif all(size(AgeWeights)==[1,N_j])
+        elseif all(size(AgeWeights_ii)==[1,N_j_temp])
             % Does not depend on transition path period
-            AgeWeights_T.(iistr)=gather(AgeWeights_ii'.*ones(1,T));
+            % Make AgeWeights a column vector
+            AgeWeights_ii=AgeWeights_ii';
         else
-            fprintf('Following error applies to agent permanent type: %s \n',Names_i{ii})
+            fprintf('Following error applies to agent permanent type: %s \n',iistr)
             error('The age weights parameter seems to be the wrong size')
         end
-    else % not a structure, so must apply to all permanent types
-        if all(size(AgeWeights)==[N_j,1])
-            % Does not depend on transition path period
-            AgeWeights_T.(iistr)=gather(AgeWeights.*ones(1,T));
-        elseif all(size(AgeWeights)==[1,N_j])
-            % Does not depend on transition path period
-            AgeWeights_T.(iistr)=gather(AgeWeights'.*ones(1,T));
+        clear AgeWeights
+        
+        % Check ParamPath to see if the AgeWeights vary over the transition
+        % (and overwrite AgeWeights_T.(iistr) if it does)
+        temp=strcmp(ParamPathNames,AgeWeightsParamNames.(iistr){1});
+        if any(temp)
+            transpathoptions.ageweightstrivial=0; % AgeWeights vary over the transition
+            [~,kk]=max(temp); % Get index for the AgeWeightsParamNames{1} in ParamPathNames
+            % Create AgeWeights_T
+            AgeWeights_T.(iistr)=ParamPath(:,ParamPathSizeVec(1,kk):ParamPathSizeVec(2,kk))'; % This will always be N_j-by-T (as transpose)
+            % Note: still leave it in ParamPath just in case it is used in AggVars or somesuch
         else
-            error('The age weights parameter seems to be the wrong size')
+            AgeWeights_T.(iistr)=repelem(AgeWeights_ii,1,T); % N_j-by-T
         end
-    end
-    % Check ParamPath to see if the AgeWeights vary over the transition
-    % (and overwrite AgeWeights_T.(iistr) if it does)
-    temp=strcmp(ParamPathNames,AgeWeightsParamNames{1});
-    if any(temp)
-        transpathoptions.ageweightstrivial=0; % AgeWeights vary over the transition
-        [~,kk]=max(temp); % Get index for the AgeWeightsParamNames{1} in ParamPathNames
-        % Create AgeWeights_T
-        AgeWeights_T.(iistr)=ParamPath(:,ParamPathSizeVec(1,kk):ParamPathSizeVec(2,kk))'; % This will always be N_j-by-T (as transpose)
-        % Note: still leave it in ParamPath just in case it is used in AggVars or somesuch
-    end
-    % Because ptypes hardcodes transpathoptions.ageweightstrivial=0 and fastOLG=1, we need
-    if N_z==0
-        AgeWeights_T.(iistr)=repelem(AgeWeights_T.(iistr),N_a,1); % simoptions.fastOLG=1 so this is (a,j)-by-1
-    else
-        AgeWeights_T.(iistr)=repmat(repelem(AgeWeights_T.(iistr),N_a,1),N_z,1); % simoptions.fastOLG=1 so this is (a,j,z)-by-1
-    end
+        clear AgeWeights_ii
 
-    %% Set up jequalOneDist_T.(iistr) [hardcodes transpathoptions.trivialjequalonedist=0 and simoptions.fastOLG=1]
-    if ~isstruct(jequalOneDist)
-        jequalOneDist_temp=gpuArray(jequalOneDist);
-    else % jequalOneDist is a structure
-        jequalOneDist_temp=gpuArray(jequalOneDist.(iistr));
-    end
-    % Check if jequalOneDistPath is a path or not (and reshape appropriately)
-    temp=size(jequalOneDist_temp);
-    if temp(end)==T % jequalOneDist depends on T
-        % transpathoptions.trivialjequalonedist=0; hardcoded for ptypes
-        if N_z==0
+        % Turn AgeWeights_T into appropriate size so that we can always just do AgentDist.*AgeWeights
+        % Currently it is N_j-by-T
+
+        % PTypes hardcodes transpathoptions.ageweightstrivial=0, we need
+        if simoptions.fastOLG==0
             if N_e==0
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a,T]);
-            else
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_e,T]); % simoptions.fastOLG==1
+                if N_z==0
+                    AgeWeights_T.(iistr)=repelem(shiftdim(AgeWeights_T.(iistr),-1),N_a,1,1); % [N_a,N_j,T]
+                else
+                    AgeWeights_T.(iistr)=repelem(shiftdim(AgeWeights_T.(iistr),-1),N_a*N_z,1,1); % [N_a*N_z,N_j,T]
+                end
+            else % N_e>0
+                if N_z==0
+                    AgeWeights_T.(iistr)=repelem(shiftdim(AgeWeights_T.(iistr),-1),N_a*N_e,1,1); % [N_a*N_e,N_j,T]
+                else
+                    AgeWeights_T.(iistr)=repelem(shiftdim(AgeWeights_T.(iistr),-1),N_a*N_z*N_e,1,1); % [N_a*N_z*N_e,N_j,T]
+                end
             end
         else
             if N_e==0
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z,T]); % simoptions.fastOLG==1
-            else
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z*N_e,T]); % simoptions.fastOLG==1
+                if N_z==0
+                    AgeWeights_T.(iistr)=repelem(AgeWeights_T.(iistr),N_a,1); % simoptions.fastOLG=1 so this is N_a*N_j-by-T
+                else
+                    AgeWeights_T.(iistr)=repmat(repelem(AgeWeights_T.(iistr),N_a,1),N_z,1); % simoptions.fastOLG=1 so this is N_a*N_j*N_z-by-T
+                end
+            else % N_e>0
+                if N_z==0
+                    AgeWeights_T.(iistr)=repelem(reshape(AgeWeights_T.(iistr),[N_j.(iistr),1,T]),N_a,N_e); % [N_a*N_j,N_e,T]
+                else
+                    AgeWeights_T.(iistr)=repmat(repelem(reshape(AgeWeights_T.(iistr),[N_j.(iistr),1,T]),N_a,1),N_z,N_e); % [N_a*N_j*N_z,N_e,T]
+                end
             end
         end
-        jequalOneDist_T.(iistr)=jequalOneDist_temp;
-    else
-        if N_z==0
+
+        %% Remove the age weights and do all the iterations. Only put the age weights back in when performing FnsToEvaluate (faster as saves putting weights in and then removing them T times)
+        % Weights are all in AgeWeights_T
+        if simoptions.fastOLG==0
+            AgentDist_init=AgentDist_init./AgeWeights_init; % AgentDist_init conveniently shaped as whatever-by-N_j
+        else
             if N_e==0
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a,1]);
+                if N_z==0
+                    AgentDist_init=AgentDist_init./repelem(AgeWeights_init',N_a,1); % remove age weights
+                else % N_e>0
+                    AgentDist_init=AgentDist_init./repmat(repelem(AgeWeights_init',N_a,1),N_z,1); % remove age weights
+                end
+            else % N_e>0
+                if N_z==0
+                    AgentDist_init=AgentDist_init./repelem(AgeWeights_init',N_a,1); % remove age weights
+                else % N_e>0
+                    AgentDist_init=AgentDist_init./repmat(repelem(AgeWeights_init',N_a,1),N_z,1); % remove age weights
+                end
+            end
+        end
+        clear AgeWeights_init
+
+    
+        %% Set up jequalOneDist_T.(iistr) [hardcodes transpathoptions.trivialjequalonedist=0]
+        if ~isstruct(jequalOneDist)
+            jequalOneDist_temp=gpuArray(jequalOneDist);
+        else % jequalOneDist is a structure
+            jequalOneDist_temp=gpuArray(jequalOneDist.(iistr));
+        end
+        % Check if jequalOneDistPath is a path or not (and reshape appropriately)
+        temp=size(jequalOneDist_temp);
+        if temp(end)==T % jequalOneDist depends on T
+            % transpathoptions.trivialjequalonedist=0; hardcoded for ptypes
+            if N_z==0
+                if N_e==0
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a,T]);
+                else
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_e,T]); % simoptions.fastOLG==1
+                end
             else
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_e,1]); % simoptions.fastOLG==1
+                if N_e==0
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z,T]); % simoptions.fastOLG==1
+                else
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z*N_e,T]); % simoptions.fastOLG==1
+                end
             end
+            jequalOneDist_T.(iistr)=jequalOneDist_temp;
         else
-            if N_e==0
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z,1]); % simoptions.fastOLG==1
+            if N_z==0
+                if N_e==0
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a,1]);
+                else
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_e,1]); % simoptions.fastOLG==1
+                end
             else
-                jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z*N_e,1]); % simoptions.fastOLG==1
+                if N_e==0
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z,1]); % simoptions.fastOLG==1
+                else
+                    jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z*N_e,1]); % simoptions.fastOLG==1 (how different than simoptions.fastOLG==0?)
+                end
             end
+            jequalOneDist_T.(iistr)=jequalOneDist_temp.*ones(1,T,'gpuArray');
         end
-        jequalOneDist_T.(iistr)=jequalOneDist_temp.*ones(1,T,'gpuArray');
+    else
+        % If no z, then N_z=1 here
+        AgentDist_init=reshape(AgentDist_init,[N_a*N_z,1]);
     end
 
+    AgentDist_initial.(iistr)=AgentDist_init;
+    clear AgentDist_init
+    
     %% Which parts of ParamPath and PricePath relate to ptype ii
     % Some ParamPath and PricePath parameters may depend on ptype
     PTypeStructure.(iistr).RelevantPricePath=ones(1,size(PricePath0,2)); % start will all relevant
@@ -741,9 +879,15 @@
 %% If using a shooting algorithm, set that up
 transpathoptions=setupGEnewprice3_shooting(transpathoptions,GeneralEqmEqns,PricePathNames,N_i,PricePathSizeVec);
 
-%% Check if using _tminus1 and/or _tplus1 variables.
+%% Check if using _tminus1 and/or _tplus1 variables, and update PTypeStructure
 if isstruct(FnsToEvaluate) && isstruct(GeneralEqmEqns)
-    [tplus1priceNames,tminus1priceNames,tminus1AggVarsNames,tminus1paramNames,tplus1pricePathkk]=inputsFindtplus1tminus1(FnsToEvaluate,GeneralEqmEqns,PricePathNames,ParamPathNames);
+    [tplus1priceNames,tminus1priceNames,tminus1AggVarsNames,tminus1paramNames,tplus1pricePathkk]=inputsFindtplus1tminus1(FnsToEvaluate,GeneralEqmEqns,PricePathNames,ParamPathNames,Names_i);
+    if isstruct(tminus1AggVarsNames)
+        AggVarsPTypes=fieldnames(tminus1AggVarsNames);
+        for ii=1:length(AggVarsPTypes)
+            PTypeStructure.(AggVarsPTypes{ii}).tminus1AggVarsNames=tminus1AggVarsNames.(AggVarsPTypes{ii});
+        end
+    end
 else
     tplus1priceNames=[];
     tminus1priceNames=[];
@@ -759,9 +903,20 @@
 use_tminus1price=0;
 if ~isempty(tminus1priceNames)
     use_tminus1price=1;
-    for ii=1:length(tminus1priceNames)
-        if ~isfield(transpathoptions.initialvalues,tminus1priceNames{ii})
-            error('Using %s as an input (to FnsToEvaluate or GeneralEqmEqns) but it is not in transpathoptions.initialvalues \n',tminus1priceNames{ii})
+    if isstruct(tminus1AggVarsNames)
+        AggVarsPTypes=fieldnames(tminus1AggVarsNames);
+        for nn=1:length(AggVarsPTypes)
+            for ii=1:length(tminus1AggVarsNames.(AggVarsPTypes{nn}))
+                if ~isfield(transpathoptions.initialvalues,tminus1AggVarsNames.(AggVarsPTypes{nn}){ii})
+                    error('Using %s as an input (to FnsToEvaluate or GeneralEqmEqns) but it is not in transpathoptions.initialvalues \n',tminus1AggVarsNames{ii})
+                end
+            end
+        end
+    else
+        for ii=1:length(tminus1AggVarsNames)
+            if ~isfield(transpathoptions.initialvalues,tminus1AggVarsNames{ii})
+                error('Using %s as an input (to FnsToEvaluate or GeneralEqmEqns) but it is not in transpathoptions.initialvalues \n',tminus1AggVarsNames{ii})
+            end
         end
     end
 end
@@ -800,9 +955,9 @@
 
 %% Shooting algorithm
 if transpathoptions.GEnewprice~=2
-    % For permanent type, 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);
+    % For permanent types, there is just one shooting command,
+    % because things like z,e, and fastOLG, as well as ExpAsset 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_initial, 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);
 
     % Switch the solution into structure for output.
     pp_indexinpricepath=zeros(1,length(PricePathNames));
diff --git a/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m b/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m
index 46b31937..3b0898d0 100644
--- a/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m
+++ b/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m
@@ -157,21 +157,37 @@
         a_gridvals=PTypeStructure.(iistr).a_gridvals;
         d_gridvals=PTypeStructure.(iistr).d_gridvals;
         aprime_gridvals=PTypeStructure.(iistr).aprime_gridvals;
-        if N_z>0
-            z_gridvals_J=PTypeStructure.(iistr).z_gridvals_J;
-            pi_z_J=PTypeStructure.(iistr).pi_z_J;
-            pi_z_J_sim=PTypeStructure.(iistr).pi_z_J_sim;
-        else
-            z_gridvals_J=[]; pi_z_J=[]; pi_z_J_sim=[];
-        end
-        if N_e>0
-            e_gridvals_J=PTypeStructure.(iistr).e_gridvals_J;
-            pi_e_J=PTypeStructure.(iistr).pi_e_J;
-            pi_e_J_sim=PTypeStructure.(iistr).pi_e_J_sim;
+        if isfinite(N_j)
+            if N_z>0
+                z_gridvals_J=PTypeStructure.(iistr).z_gridvals_J;
+                pi_z_J=PTypeStructure.(iistr).pi_z_J;
+                pi_z_J_sim=PTypeStructure.(iistr).pi_z_J_sim;
+            else
+                z_gridvals_J=[]; pi_z_J=[]; pi_z_J_sim=[];
+            end
+            if N_e>0
+                e_gridvals_J=PTypeStructure.(iistr).e_gridvals_J;
+                pi_e_J=PTypeStructure.(iistr).pi_e_J;
+                pi_e_J_sim=PTypeStructure.(iistr).pi_e_J_sim;
+            else
+                e_gridvals_J=[]; pi_e_J=[]; pi_e_J_sim=[];
+            end
+            ze_gridvals_J_fastOLG=PTypeStructure.(iistr).ze_gridvals_J_fastOLG;
         else
-            e_gridvals_J=[]; pi_e_J=[]; pi_e_J_sim=[];
+            if N_z>0
+                z_grid=PTypeStructure.(iistr).z_grid;
+                pi_z=PTypeStructure.(iistr).pi_z;
+            else
+                z_grid=[]; pi_z=[];
+            end
+            if N_e>0
+                e_gridvals=PTypeStructure.(iistr).e_gridvals;
+                pi_e=PTypeStructure.(iistr).pi_e;
+            else
+                e_gridvals=[]; pi_e=[];
+            end
+            % ze_gridvals_fastOLG=PTypeStructure.(iistr).ze_gridvals_fastOLG;
         end
-        ze_gridvals_J_fastOLG=PTypeStructure.(iistr).ze_gridvals_J_fastOLG;
         if transpathoptions.fastOLG==1
             exceptlastj=PTypeStructure.(iistr).exceptlastj;
             exceptfirstj=PTypeStructure.(iistr).exceptfirstj;
@@ -191,10 +207,16 @@
         FnsToEvaluate=PTypeStructure.(iistr).FnsToEvaluate;
         FnsToEvaluateParamNames=PTypeStructure.(iistr).FnsToEvaluateParamNames;
         AggVarNames=PTypeStructure.(iistr).AggVarNames;
-
+        if isfield(PTypeStructure.(iistr), 'tminus1AggVarsNames')
+            use_tminus1AggVars=1;
+            tminus1AggVarsNames=PTypeStructure.(iistr).tminus1AggVarsNames;
+        else
+            use_tminus1AggVars=0;
+        end
 
         % Following few lines I would normally do outside of the while loop, but have to set them for each ptype
         % AgentDist=AgentDist_initial.(iistr);
+        % WARNING: The following would overwrite themselves next iteration
         % V_final=V_final.(iistr);
         % AgeWeights_T=AgeWeights_T.(iistr);
         % jequalOneDist_T=jequalOneDist_T.(iistr);
@@ -206,10 +228,23 @@
         % PricePathSizeVec_ii, ParamPathSizeVec_ii
         
         % For current ptype, do the backward iteration of V and Policy, then forward iterate agent dist and get the AggVarsPath
-        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);
+        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_d1,n_d2,n_a1,n_a2,N_z,n_z,N_e,n_e,N_j,d2_grid,a1_gridvals,a2_grid,d_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, aprimeFn, FnsToEvaluate, Parameters, DiscountFactorParamNames, ReturnFnParamNames, aprimeFnParamNames, 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
+            if vfoptions.experienceasset==1
+                error("not yet implemented")
+                AggVarsPath=TransitionPath_InfHorz_PType_ExpAsset_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, PricePathSizeVec_ii, ParamPathSizeVec_ii, use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, transpathoptions, vfoptions, simoptions);
+            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, PricePathSizeVec_ii, ParamPathSizeVec_ii, use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, transpathoptions, vfoptions, simoptions);
+            end
+        end
         % AggVarsPath=zeros(length(FnsToEvaluate),T-1);
         
-        AggVarsFullPath(PTypeStructure.(iistr).WhichFnsForCurrentPType,:,ii)=AggVarsPath;
+        AggVarsFullPath(logical(PTypeStructure.(iistr).WhichFnsForCurrentPType),:,ii)=AggVarsPath;
 
     end % done loop over ii
     
@@ -222,7 +257,22 @@
     %% Do the general eqm conditions and create PricePathNew based on these
     if all(transpathoptions.GEptype==0)
         GECondnPath=zeros(T,length(GEeqnNames));
-            
+
+        % Restore all AggVarNames and tminus1AggVarsNames for GEeqns
+        AggVarNames=cell(1,N_i);
+        tminus1AggVarsNames=cell(1,N_i);
+        use_tminus1AggVars=0;
+        for ii=1:N_i
+            iistr=PTypeStructure.Names_i{ii};
+            AggVarNames{ii}=PTypeStructure.(iistr).AggVarNames;
+            if isfield(PTypeStructure.(iistr), 'tminus1AggVarsNames')
+                use_tminus1AggVars=1;
+                tminus1AggVarsNames{ii}=PTypeStructure.(iistr).tminus1AggVarsNames;
+            end
+        end
+        AggVarNames=vertcat(AggVarNames{:});
+        tminus1AggVarsNames=vertcat(tminus1AggVarsNames{:});
+        
         % Parameters that may be relevant to General Eqm
         Parameters=PTypeStructure.ParametersRaw;
 
@@ -289,7 +339,23 @@
         end % Done loop over tt, evaluating the GE conditions
     else % Some GE conditions depend on PType
         GECondnPath=zeros(T,nGeneralEqmEqns_acrossptypes);
+        error("need to fit these loops together")
 
+        % Restore AggVarNames and tminus1AggVarsNames by PType for GEeqns
+        AggVarNames=cell(1,N_i);
+        tminus1AggVarsNames=cell(1,N_i);
+        use_tminus1AggVars=0;
+        for ii=1:N_i
+            iistr=PTypeStructure.Names_i{ii};
+            AggVarNames{ii}=PTypeStructure.(iistr).AggVarNames;
+            tminus1AggVarsNames{ii}=PTypeStructure.(iistr).tminus1AggVarsNames;
+            if ~isempty(tminus1AggVarsNames{ii})
+                use_tminus1AggVars=1;
+            end
+        end
+        AggVarNames=vertcat(AggVarNames{:});
+        tminus1AggVarsNames=vertcat(tminus1AggVarsNames{:});
+        
         % Parameters that may be relevant to General Eqm
         Parameters=PTypeStructure.ParametersRaw;
         for ii=1:N_i
diff --git a/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_e_raw.m b/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_e_raw.m
index b3b454ca..371e8529 100644
--- a/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_e_raw.m
+++ b/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_e_raw.m
@@ -9,7 +9,10 @@
 else
     l_d=length(n_d);
 end
-l_aprime=length(n_a);
+l_a=length(n_a);
+l_aprime=l_a;
+N_ze=prod(n_z)*prod(n_e);
+l_ze=length(n_z)+length(n_e);
 
 %%
 % Shapes:
@@ -192,7 +195,7 @@
     end
 
     %% AggVars
-    AggVars=EvalFnOnAgentDist_AggVars_FHorz_fastOLG(AgentDist.*AgeWeights, PolicyValuesPath(:,:,:,1:l_d,tt), PolicyValuesPath(:,:,:,l_d+1:end,tt), FnsToEvaluateCell,FnsToEvaluateParamNames,AggVarNames,Parameters,N_j,l_d,l_a,l_a,l_ze,N_a,N_ze,a_gridvals,ze_gridvals_J_fastOLG,1);
+    AggVars=EvalFnOnAgentDist_AggVars_FHorz_fastOLG(AgentDist.*AgeWeights, PolicyValuesPath(:,:,:,1:l_d,tt), PolicyValuesPath(:,:,:,l_d+1:end,tt), FnsToEvaluateCell,FnsToEvaluateParamNames,AggVarNames,Parameters,N_j,l_d,l_aprime,l_a,l_ze,N_a,N_ze,a_gridvals,ze_gridvals_J_fastOLG,1);
     for ff=1:length(AggVarNames)
         Parameters.(AggVarNames{ff})=AggVars.(AggVarNames{ff}).Mean;
     end
diff --git a/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_fastOLG_e_raw.m b/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_fastOLG_e_raw.m
index 0f197e6c..6d1461a3 100644
--- a/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_fastOLG_e_raw.m
+++ b/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_fastOLG_e_raw.m
@@ -12,6 +12,9 @@
 l_a=length(n_a);
 l_aprime=length(n_a);
 l_z=length(n_z);
+l_e=length(n_e);
+l_ze=l_z+l_e;
+N_ze=N_z*N_e;
 
 %%
 % fastOLG so everything is (a,j,z,e)
@@ -158,7 +161,7 @@
     end
     
     %% AggVars
-    AggVars=EvalFnOnAgentDist_AggVars_FHorz_fastOLG(AgentDist.*AgeWeights, PolicyValuesPath(:,:,:,1:l_d,tt), PolicyValuesPath(:,:,:,l_d+1:end,tt), FnsToEvaluateCell,FnsToEvaluateParamNames,AggVarNames,Parameters,N_j,l_d,l_a,l_a,l_ze,N_a,N_ze,a_gridvals,ze_gridvals_J_fastOLG,1);
+    AggVars=EvalFnOnAgentDist_AggVars_FHorz_fastOLG(AgentDist.*AgeWeights, PolicyValuesPath(:,:,:,1:l_d,tt), PolicyValuesPath(:,:,:,l_d+1:end,tt), FnsToEvaluateCell,FnsToEvaluateParamNames,AggVarNames,Parameters,N_j,l_d,l_aprime,l_a,l_ze,N_a,N_ze,a_gridvals,ze_gridvals_J_fastOLG,1);
     for ff=1:length(AggVarNames)
         Parameters.(AggVarNames{ff})=AggVars.(AggVarNames{ff}).Mean;
     end
diff --git a/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_raw.m b/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_raw.m
index d8cc9222..7c39a925 100644
--- a/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_raw.m
+++ b/TransitionPaths/FHorz/PType/singlepath/TransitionPath_FHorz_PType_singlepath_raw.m
@@ -9,7 +9,9 @@
 else
     l_d=length(n_d);
 end
-l_aprime=length(n_a);
+l_a=length(n_a);
+l_aprime=l_a;
+l_z=length(n_z);
 
 %%
 % Shapes:
@@ -185,7 +187,7 @@
     end
 
     %% AggVars
-    AggVars=EvalFnOnAgentDist_AggVars_FHorz_fastOLG(AgentDist.*AgeWeights, PolicyValuesPath(:,:,:,1:l_d,tt), PolicyValuesPath(:,:,:,l_d+1:end,tt), FnsToEvaluateCell,FnsToEvaluateParamNames,AggVarNames,Parameters,N_j,l_d,l_a,l_a,l_z,N_a,N_z,a_gridvals,z_gridvals_J_fastOLG,1);
+    AggVars=EvalFnOnAgentDist_AggVars_FHorz_fastOLG(AgentDist.*AgeWeights, PolicyValuesPath(:,:,:,1:l_d,tt), PolicyValuesPath(:,:,:,l_d+1:end,tt), FnsToEvaluateCell,FnsToEvaluateParamNames,AggVarNames,Parameters,N_j,l_d,l_aprime,l_a,l_z,N_a,N_z,a_gridvals,z_gridvals_J_fastOLG,1);
     for ff=1:length(AggVarNames)
         Parameters.(AggVarNames{ff})=AggVars.(AggVarNames{ff}).Mean;
     end
diff --git a/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_shooting.m b/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_shooting.m
index ffa4987c..c727ec3d 100644
--- a/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_shooting.m
+++ b/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_shooting.m
@@ -31,6 +31,7 @@
     AgentDist_initial.(iistr)=reshape(StationaryDist_init.(iistr),[N_a*N_z,1]);
 end
 PricePathNew=zeros(size(PricePathOld),'gpuArray'); PricePathNew(T,:)=PricePathOld(T,:);
+GEeqnNames=fieldnames(GeneralEqmEqns);
 GEcondnPath=zeros(T-1,length(GEeqnNames),'gpuArray');
 
 
diff --git a/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m b/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m
new file mode 100644
index 00000000..647fcac6
--- /dev/null
+++ b/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m
@@ -0,0 +1,168 @@
+function AggVarsPath=TransitionPath_InfHorz_PType_singlepath(PricePathOld, ParamPath, PricePathNames, ParamPathNames, T, V_final, AgentDist_initial, ...
+    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, PricePathSizeVec, ParamPathSizeVec, ...
+    use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, ...
+    transpathoptions, vfoptions, simoptions)
+% PricePathOld is matrix of size T-by-'number of prices'
+% ParamPath is matrix of size T-by-'number of parameters that change over path'
+
+% Remark to self: No real need for T as input, as this is anyway the length of PricePathOld
+
+% For this agent type, first go back through the value & policy fns.
+% Then forwards through agent dist and agg vars.
+
+l_a=length(n_a);
+PolicyIndexesPath=zeros(l_d+l_a,N_a,N_z,T-1,'gpuArray'); %Periods 1 to T-1
+
+% This and much of the rest of this code borrowed from TransitionPath_InfHorz_shooting.m
+if simoptions.gridinterplayer==0
+    II1=(1:1:N_a*N_z); % Index for this period (a,z)
+    IIones=ones(N_a*N_z,1); % Next period 'probabilities'
+elseif simoptions.gridinterplayer==1
+    PolicyProbsPath=zeros(N_a*N_z,2,T-1,'gpuArray'); % preallocate
+    II2=([1:1:N_a*N_z; 1:1:N_a*N_z]'); % Index for this period (a,z), note the 2 copies
+end
+
+if size(V_final)~=[N_a,N_z]
+    error("V_final wrong shape")
+    V_final=reshape(V_final,[N_a,N_z]);
+end
+
+%First, go from T-1 to 1 calculating the Value function and Optimal
+%policy function at each step. Since we won't need to keep the value
+%functions for anything later we just store the next period one in
+%Vnext, and the current period one to be calculated in V
+Vnext=V_final;
+for ttr=1:T-1 %so t=T-i
+    % The following digs deeper into PricePathOld and ParamPath in
+    % FHorz case--check it
+    for kk=1:length(PricePathNames)
+        Parameters.(PricePathNames{kk})=PricePathOld(T-ttr,kk);
+    end
+    for kk=1:length(ParamPathNames)
+        Parameters.(ParamPathNames{kk})=ParamPath(T-ttr,kk);
+    end
+    
+    [V, Policy]=ValueFnIter_InfHorz_TPath_SingleStep(Vnext,n_d,n_a,n_z,d_grid, a_grid, z_grid, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
+    % The VKron input is next period value fn, the VKron output is this period.
+    % Policy is kept in the form where it is just a single-value in (d,a')
+
+    PolicyIndexesPath(:,:,:,T-ttr)=Policy;
+    Vnext=V;
+    
+end
+% Free up space on GPU by deleting things no longer needed
+clear V Vnext
+
+%% Modify PolicyIndexesPath into forms needed for forward iteration
+% Create version of PolicyIndexesPath in form we want for the agent distribution iteration
+% Creates PolicyaprimezPath, and when using grid interpolation layer also PolicyProbsPath 
+if isscalar(n_a)
+    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:),[N_a*N_z,T-1]); % aprime index
+elseif length(n_a)==2
+    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:)+n_a(1)*(PolicyIndexesPath(l_d+2,:,:,:)-1),[N_a*N_z,T-1]);
+elseif length(n_a)==3
+    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:)+n_a(1)*(PolicyIndexesPath(l_d+2,:,:,:)-1)+n_a(1)*n_a(2)*(PolicyIndexesPath(l_d+3,:,:,:)-1),[N_a*N_z,T-1]);
+elseif length(n_a)==4
+    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:)+n_a(1)*(PolicyIndexesPath(l_d+2,:,:,:)-1)+n_a(1)*n_a(2)*(PolicyIndexesPath(l_d+3,:,:,:)-1)+n_a(1)*n_a(2)*n_a(3)*(PolicyIndexesPath(l_d+4,:,:,:)-1),[N_a*N_z,T-1]);
+end
+PolicyaprimezPath=PolicyaprimePath+repelem(N_a*gpuArray(0:1:N_z-1)',N_a,1);
+if simoptions.gridinterplayer==1
+    PolicyaprimezPath=reshape(PolicyaprimezPath,[N_a*N_z,1,T-1]); % reinterpret this as lower grid index
+    PolicyaprimezPath=repelem(PolicyaprimezPath,1,2,1); % create copy that will be the upper grid index
+    PolicyaprimezPath(:,2,:)=PolicyaprimezPath(:,2,:)+1; % upper grid index
+    PolicyProbsPath(:,2,:)=reshape(PolicyIndexesPath(l_d+l_aprime+1,:,:),[N_a*N_z,1,T-1]); % L2 index
+    PolicyProbsPath(:,2,:)=(PolicyProbsPath(:,2,:)-1)/(1+simoptions.ngridinterp); % probability of upper grid point
+    PolicyProbsPath(:,1,:)=1-PolicyProbsPath(:,2,:); % probability of lower grid point
+end
+% Create PolicyValuesPath from PolicyIndexesPath for use in calculating model stats
+PolicyValuesPath=PolicyInd2Val_InfHorz_TPath(PolicyIndexesPath,n_d,n_a,n_z,T-1,d_grid,a_grid,vfoptions,1);
+PolicyValuesPath=permute(reshape(PolicyValuesPath,[size(PolicyValuesPath,1),N_a,N_z,T-1]),[2,3,1,4]); %[N_a,N_z,l_d+l_a,T-1]
+
+%Now we have the full PolicyIndexesPath, we go forward in time from 1
+%to T using the policies to update the agents distribution generating a
+%new price path
+%Call AgentDist the current periods distn
+
+AgentDist=sparse(gather(reshape(AgentDist_initial,[N_a*N_z,1])));
+AggVarsPath=zeros(length(FnsToEvaluate),T-1);
+pi_z_sparse=sparse(gather(pi_z)); % Need full pi_z for value fn, and sparse for agent dist
+
+for tt=1:T-1
+    %% Setup the Parameters for period tt
+
+    % Get t-1 PricePath and ParamPath before we update them
+    if use_tminus1price==1
+        for pp=1:length(tminus1priceNames)
+            if tt>1
+                Parameters.([tminus1priceNames{pp},'_tminus1'])=Parameters.(tminus1priceNames{pp});
+            else
+                Parameters.([tminus1priceNames{pp},'_tminus1'])=transpathoptions.initialvalues.(tminus1priceNames{pp});
+            end
+        end
+    end
+    if use_tminus1params==1
+        for pp=1:length(tminus1paramNames)
+            if tt>1
+                Parameters.([tminus1paramNames{pp},'_tminus1'])=Parameters.(tminus1paramNames{pp});
+            else
+                Parameters.([tminus1paramNames{pp},'_tminus1'])=transpathoptions.initialvalues.(tminus1paramNames{pp});
+            end
+        end
+    end
+    % Get t-1 AggVars before we update them
+    if use_tminus1AggVars==1
+        for pp=1:length(tminus1AggVarsNames)
+            if tt>1
+                % The AggVars have not yet been updated, so they still contain previous period values
+                Parameters.([tminus1AggVarsNames{pp},'_tminus1'])=Parameters.(tminus1AggVarsNames{pp});
+            else
+                Parameters.([tminus1AggVarsNames{pp},'_tminus1'])=transpathoptions.initialvalues.(tminus1AggVarsNames{pp});
+            end
+        end
+    end
+
+    % Update current PricePath and ParamPath
+    for kk=1:length(PricePathNames)
+        Parameters.(PricePathNames{kk})=PricePathOld(tt,PricePathSizeVec(1,kk):PricePathSizeVec(2,kk));
+    end
+    for kk=1:length(ParamPathNames)
+        Parameters.(ParamPathNames{kk})=ParamPath(tt,ParamPathSizeVec(1,kk):ParamPathSizeVec(2,kk));
+    end
+
+    % Get t+1 PricePath
+    if use_tplus1price==1
+        for pp=1:length(tplus1priceNames)
+            kk=tplus1pricePathkk(pp);
+            Parameters.([tplus1priceNames{pp},'_tplus1'])=PricePathOld(tt+1,PricePathSizeVec(1,kk):PricePathSizeVec(2,kk)); % Make is so that the time t+1 variables can be used
+        end
+    end
+    
+    %% Get the current optimal policy, and iterate the agent dist
+    if simoptions.gridinterplayer==0
+        AgentDistnext=AgentDist_InfHorz_TPath_SingleStep(AgentDist,PolicyaprimezPath(:,tt),II1,IIones,N_a,N_z,pi_z_sparse);
+    elseif simoptions.gridinterplayer==1
+        AgentDistnext=AgentDist_InfHorz_TPath_SingleStep_nProbs(AgentDist,PolicyaprimezPath(:,:,tt),II2,PolicyProbsPath(:,:,tt),N_a,N_z,pi_z_sparse);
+    end
+
+    %% AggVars
+    AggVars_Means=EvalFnOnAgentDist_InfHorz_TPath_SingleStep_AggVars(gpuArray(full(AgentDist)), PolicyValuesPath(:,:,:,tt), FnsToEvaluate, Parameters, FnsToEvaluateParamNames, AggVarNames, n_a, n_z, a_gridvals, z_grid,1);
+    AggVars=zeros(length(AggVars_Means),1);
+    if length(fieldnames(AggVars_Means))~=length(AggVarNames)
+        error(["AggVar length disparity:";"---------";AggVarNames;"--- vs ---";fieldnames(AggVars_Means)]);
+    end
+    for ii=1:length(AggVarNames)
+        AggVars(ii)=AggVars_Means.(AggVarNames{ii}).Mean;
+        Parameters.(AggVarNames{ii})=AggVars(ii);
+    end
+
+    % Do nothing with IntermediateEqns and GeneralEqmEqns as they are outside PType scope (they are where all the PTypes meet).
+
+    AgentDist=AgentDistnext;
+
+    AggVarsPath(:,tt)=AggVars;
+end
+
+
+end
diff --git a/TransitionPaths/InfHorz/subcodes/AgentDistSingleStep/AgentDist_InfHorz_TPath_SingleStep.m b/TransitionPaths/InfHorz/subcodes/AgentDistSingleStep/AgentDist_InfHorz_TPath_SingleStep.m
index 17809915..27998ce7 100644
--- a/TransitionPaths/InfHorz/subcodes/AgentDistSingleStep/AgentDist_InfHorz_TPath_SingleStep.m
+++ b/TransitionPaths/InfHorz/subcodes/AgentDistSingleStep/AgentDist_InfHorz_TPath_SingleStep.m
@@ -3,7 +3,7 @@
 % 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 
+% 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
diff --git a/TransitionPaths/InfHorz/subcodes/ValueFnSingleStep/ValueFnIter_InfHorz_TPath_SingleStep.m b/TransitionPaths/InfHorz/subcodes/ValueFnSingleStep/ValueFnIter_InfHorz_TPath_SingleStep.m
index 3427ee74..e9722b22 100644
--- a/TransitionPaths/InfHorz/subcodes/ValueFnSingleStep/ValueFnIter_InfHorz_TPath_SingleStep.m
+++ b/TransitionPaths/InfHorz/subcodes/ValueFnSingleStep/ValueFnIter_InfHorz_TPath_SingleStep.m
@@ -4,8 +4,8 @@
 % vfoptions must be already fully set up (this command is for internal use only so it should be)
 
 N_d=prod(n_d);
-% N_a=prod(n_a);
-% N_z=prod(n_z);
+N_a=prod(n_a);
+N_z=prod(n_z);
 
 %% 
 if strcmp(vfoptions.exoticpreferences,'QuasiHyperbolic')
@@ -21,9 +21,11 @@
 if vfoptions.gridinterplayer==0
     if vfoptions.divideandconquer==0
         if N_d==0
-            [VKron,PolicyKron]=ValueFnIter_InfHorz_TPath_SingleStep_nod_raw(VKron,n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
+            [VKron,Policy]=ValueFnIter_InfHorz_TPath_SingleStep_nod_raw(VKron,n_a, n_z, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
+            PolicyKron=reshape(Policy,[size(Policy,1),N_a,N_z]);
         else
-            [VKron, PolicyKron]=ValueFnIter_InfHorz_TPath_SingleStep_raw(VKron,n_d,n_a,n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
+            [VKron, Policy2]=ValueFnIter_InfHorz_TPath_SingleStep_raw(VKron,n_d,n_a,n_z, d_grid, a_grid, z_gridvals, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
+            PolicyKron=reshape(Policy2,[size(Policy2,1),N_a,N_z]);
         end
     elseif vfoptions.divideandconquer==1
         if isscalar(n_a)
diff --git a/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m b/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m
index 52035383..4eb92486 100644
--- a/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m
+++ b/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m
@@ -71,6 +71,7 @@
     % Note: Internally PricePathOld is matrix of size T-by-'number of prices'.
     % ParamPath is matrix of size T-by-'number of parameters that change over the transition path'.
     % Actually, some of those prices are 1-by-N_j or N_i or both, so is more subtle than this.
+    Names_i=fieldnames(N_j);
     PricePathNames=fieldnames(PricePathStruct);
     PricePathSizeVec=zeros(1,length(PricePathNames)); % Allows for a given price param to depend on age (or permanent type)
     for pp=1:length(PricePathNames)
@@ -79,15 +80,21 @@
             temp=PricePathStruct.(PricePathNames{pp}).(tempptypenames{1});
             tempsize=size(temp);
             PricePathSizeVec(pp)=length(tempptypenames)*tempsize(tempsize~=T); % Get the dimension which is not T
-            if ~any(PricePathSizeVec(pp)==[1,N_i,N_j,N_i*N_j])
-                error(['PricePath for ', PricePathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+            for ii=1:N_i
+                N_j_temp=N_j.(Names_i{ii});
+                if isfinite(N_j_temp) && ~any(PricePathSizeVec(pp)==[1,N_i,N_j_temp,N_i*N_j_temp])
+                    error(['PricePath for ', PricePathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+                end
             end
         else
             temp=PricePathStruct.(PricePathNames{pp});
             tempsize=size(temp);
             PricePathSizeVec(pp)=tempsize(tempsize~=T); % Get the dimension which is not T
-            if ~any(PricePathSizeVec(pp)==[1,N_i,N_j])
-                error(['PricePath for ', PricePathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+            for ii=1:N_i
+                N_j_temp=N_j.(Names_i{ii});
+                if isfinite(N_j_temp) && ~any(PricePathSizeVec(pp)==[1,N_i,N_j_temp])
+                    error(['PricePath for ', PricePathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+                end
             end
         end
     end
@@ -136,16 +143,22 @@
             temp=ParamPathStruct.(ParamPathNames{pp}).*(tempptypenames{1});
             tempsize=size(temp);
             ParamPathSizeVec(pp)=length(tempptypenames)*tempsize(tempsize~=T); % Get the dimension which is not T
-            if ~any(ParamPathSizeVec(pp)==[1,N_i,N_j,N_i*N_j])
-                error(['ParamPath for ', ParamPathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+            for ii=1:N_i
+                N_j_temp=N_j.(Names_i{ii});
+                if isfinite(N_j_temp) && ~any(ParamPathSizeVec(pp)==[1,N_i,N_j_temp,N_i*N_j_temp])
+                    error(['ParamPath for ', ParamPathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+                end
             end
 
         else
             temp=ParamPathStruct.(ParamPathNames{pp});
             tempsize=size(temp);
             ParamPathSizeVec(pp)=tempsize(tempsize~=T); % Get the dimension which is not T
-            if ~any(ParamPathSizeVec(pp)==[1,N_i,N_j])
-                error(['ParamPath for ', ParamPathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+            for ii=1:N_i
+                N_j_temp=N_j.(Names_i{ii});
+                if isfinite(N_j_temp) && ~any(ParamPathSizeVec(pp)==[1,N_i,N_j_temp])
+                    error(['ParamPath for ', ParamPathNames{pp}, ' appears to be the wrong size (should be 1-by-T or N_j-by-T or N_i-by-T)'])
+                end
             end
         end
     end
diff --git a/TransitionPaths/Subcodes/inputsFindtplus1tminus1.m b/TransitionPaths/Subcodes/inputsFindtplus1tminus1.m
index 1a744443..5865ebd6 100644
--- a/TransitionPaths/Subcodes/inputsFindtplus1tminus1.m
+++ b/TransitionPaths/Subcodes/inputsFindtplus1tminus1.m
@@ -1,20 +1,50 @@
-function [tplus1priceNames,tminus1priceNames,tminus1AggVarsNames,tminus1paramNames,tplus1pricePathkk]=inputsFindtplus1tminus1(FnsToEvaluate,GeneralEqmEqns,PricePathNames,ParamPathNames)
+function [tplus1priceNames,tminus1priceNames,tminus1AggVarsNames,tminus1paramNames,tplus1pricePathkk]=inputsFindtplus1tminus1(FnsToEvaluate,GeneralEqmEqns,PricePathNames,ParamPathNames,Names_i)
 % Subscript that is used to determine if there are any '_tminus1' or
 % '_tplus1' variables used as inputs to FnsToEvaluate or GeneralEqmEqns
 % (Used as part of transition path codes)
 % Look for the _tplus1 in PricePath, and for _tminus1 in AggVars, PricePath, and ParamPath
 
-% Check for using _tminus1 (for price or AggVars) or _tplus1
-FnInputNames={};
+if ~exist('Names_i','var')
+    Names_i={};
+end
+
+% Check for using _tminus1 (for price or AggVars) or _tplus1, honoring PType
+FnInputNames={}; % FnNames beloning to no PType
+PTypeFnInputNames=struct(); % Built as we find them
 ninputs=0;
 % Get all the input names for FnsToEvaluate
 AggVarNames=fieldnames(FnsToEvaluate);
+PTypeAggVarNames=struct();
 for ff=1:length(AggVarNames)
-    temp=getAnonymousFnInputNames(FnsToEvaluate.(AggVarNames{ff}));
+    if ~isempty(Names_i)
+        temp={};
+        for ii=1:length(Names_i)
+            if isfield(FnsToEvaluate.(AggVarNames{ff}), Names_i{ii})
+                if ~isfield(PTypeAggVarNames, Names_i{ii})
+                    PTypeAggVarNames.(Names_i{ii})={};
+                    PTypeFnInputNames.(Names_i{ii})={};
+                end
+                PTypeAggVarNames.(Names_i{ii})=[PTypeAggVarNames.(Names_i{ii});AggVarNames{ff}];
+                temp=getAnonymousFnInputNames(FnsToEvaluate.(AggVarNames{ff}).(Names_i{ii}));
+                PTypeFnInputNames.(Names_i{ii})=[PTypeFnInputNames.(Names_i{ii});temp]; % Note, this will include the (d,aprime,a,z), but that is irrelevant to our current purposes
+                break
+            end
+        end
+        if isempty(temp)
+            error(['Could not find', AggVarNames{ff}, 'with PType in {', Names_i, '}'])
+        end
+    else
+        temp=getAnonymousFnInputNames(FnsToEvaluate.(AggVarNames{ff}));
+        FnInputNames={FnInputNames{:},temp{:}}; % Note, this will include the (d,aprime,a,z), but that is irrelevant to our current purposes
+    end
     tempninputs=length(temp);
-    FnInputNames={FnInputNames{:},temp{:}}; % Note, this will include the (d,aprime,a,z), but that is irrelevant to our current purposes
     ninputs=ninputs+tempninputs;
 end
+for ii=1:length(Names_i) % If Names_i is empty this loop is not executed
+    name_matches=ismember(AggVarNames, PTypeAggVarNames.(Names_i{ii}));
+    AggVarNames(name_matches)=[];
+end
+
 % Get all the input names for GeneralEqmEqns
 GEeqnNames=fieldnames(GeneralEqmEqns);
 for gg=1:length(GEeqnNames)
@@ -25,6 +55,9 @@
 end
 % Now, get rid of all the duplicates
 FnInputNames=unique(FnInputNames);
+for ii=1:length(Names_i) % If Names_i is empty this loop is not executed
+    PTypeFnInputNames.(Names_i{ii})=unique(PTypeFnInputNames.(Names_i{ii}));
+end
 
 % Find any which are _tplus1 or _tminus1
 tplus1Names={};
@@ -53,6 +86,10 @@
 tminus1priceNames={}; ntminus1prices=0;
 tminus1paramNames={}; ntminus1params=0;
 tminus1AggVarsNames={}; ntminus1AggVars=0;
+for ii=1:length(Names_i)
+    PType_tminus1AggVarsNames.(Names_i{ii})={};
+    PType_ntminus1AggVars.(Names_i{ii})=0;
+end
 tplus1pricePathkk=[];
 
 % Check that they are prices, otherwise error
@@ -87,14 +124,38 @@
         end
     end
 end
-for ii=1:ntminus1
-    for kk=1:length(AggVarNames)
-        if strcmp(tminus1Names{ii},AggVarNames{kk})
-            ntminus1AggVars=ntminus1AggVars+1;
-            tminus1AggVarsNames{ntminus1AggVars}=tminus1Names{ii};
-            tminus1UsedAsPriceOrAggVar(ii)=1;
+if isempty(Names_i)
+    for ii=1:ntminus1
+        for kk=1:length(AggVarNames)
+            if strcmp(tminus1Names{ii},AggVarNames{kk})
+                ntminus1AggVars=ntminus1AggVars+1;
+                tminus1AggVarsNames{ntminus1AggVars}=tminus1Names{ii};
+                tminus1UsedAsPriceOrAggVar(ii)=1;
+            end
+        end
+    end
+    % Return the array of tminus1AggVarsNames
+else
+    for ii=1:ntminus1
+        for nn=1:length(Names_i)
+            AggVarNames_nn=PTypeAggVarNames.(Names_i{nn});
+            for kk=1:length(AggVarNames_nn)
+                if strcmp(tminus1Names{ii},AggVarNames_nn{kk})
+                    PType_ntminus1AggVars.(Names_i{nn})=PType_ntminus1AggVars.(Names_i{nn})+1;
+                    PType_tminus1AggVarsNames.(Names_i{nn}){PType_ntminus1AggVars.(Names_i{nn})}=tminus1Names{ii};
+                    tminus1UsedAsPriceOrAggVar(ii)=1;
+                end
+            end
         end
     end
+    % It is easier to leave empty values than remove/retest the fields
+    % for ii=1:length(Names_i)
+    %     if isempty(PType_tminus1AggVarsNames.(Names_i{ii}))
+    %         PType_tminus1AggVarsNames=rmfield(PType_tminus1AggVarsNames, Names_i{ii});
+    %     end
+    % end
+    % Return the PType structure of the arrays of tminus1AggVarsNames
+    tminus1AggVarsNames=PType_tminus1AggVarsNames;
 end
 % Check that they have all been used, otherwise error
 if prod(tplus1UsedAsPriceOrAggVar)==0
diff --git a/ValueFnIter/InfHorz/ValueFnIter_nod_HowardGreedy_raw.m b/ValueFnIter/InfHorz/ValueFnIter_nod_HowardGreedy_raw.m
index 6c1fefa4..62c78ed9 100644
--- a/ValueFnIter/InfHorz/ValueFnIter_nod_HowardGreedy_raw.m
+++ b/ValueFnIter/InfHorz/ValueFnIter_nod_HowardGreedy_raw.m
@@ -56,14 +56,5 @@
     warning('Value fn iteration has stopped due to reaching the maximum number of iterations (not due to convergence); can be set by vfoptions.maxiter.')
 end
 
-%% Howards greedy cannot solve models where V contains values of -Inf. Can kind of test for this by looking for -Inf in Ftemp
-if any(~isfinite(Ftemp))
-    warning('Howards-greedy cannot be used for models where V contains values of -Inf. This model looks like it may be one where V takes a value of -Inf at some points in the state-space. Consider checking solution against that with vfoptions.howardsgreedy=0')
-end
-
-if tempcounter>=maxiter
-    warning('Value fn iteration has stopped due to reaching the maximum number of iterations (not due to convergence); can be set by vfoptions.maxiter.')
-end
-
 
 end
diff --git a/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m b/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m
index 766a563c..5a756eae 100644
--- a/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m
+++ b/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m
@@ -1,4 +1,4 @@
-function [VPath,PolicyPath]=ValueFnOnTransPath_Case1(PricePath, ParamPath, T, V_final, Policy_final, Parameters, n_d, n_a, n_z, pi_z, d_grid, a_grid,z_grid, DiscountFactorParamNames, ReturnFn, transpathoptions, vfoptions)
+function [VPath,PolicyPath]=ValueFnOnTransPath_Case1(PricePath, ParamPath, T, V_final, Policy_final, Parameters, n_d, n_a, n_z, d_grid, a_grid,z_grid, pi_z, DiscountFactorParamNames, ReturnFn, transpathoptions, vfoptions)
 % transpathoptions, vfoptions and simoptions are optional inputs
 
 %% Check which transpathoptions have been used, set all others to defaults 
diff --git a/ValueFnIter/TransPathFHorz/ValueFnOnTransPath_Case1_FHorz_PType.m b/ValueFnIter/TransPathFHorz/ValueFnOnTransPath_Case1_FHorz_PType.m
index f0ae812a..4c02800b 100644
--- a/ValueFnIter/TransPathFHorz/ValueFnOnTransPath_Case1_FHorz_PType.m
+++ b/ValueFnIter/TransPathFHorz/ValueFnOnTransPath_Case1_FHorz_PType.m
@@ -190,7 +190,11 @@
         end
     end
 
-    [VPath_ii,PolicyPath_ii]=ValueFnOnTransPath_Case1_FHorz(PricePath_temp, ParamPath_temp, T, V_final_temp, Policy_final_temp, Parameters_temp, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp,pi_z_temp, DiscountFactorParamNames_temp, ReturnFn_temp, transpathoptions_temp, vfoptions_temp);
+    if isfinite(N_j_temp)
+        [VPath_ii,PolicyPath_ii]=ValueFnOnTransPath_Case1_FHorz(PricePath_temp, ParamPath_temp, T, V_final_temp, Policy_final_temp, Parameters_temp, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp,pi_z_temp, DiscountFactorParamNames_temp, ReturnFn_temp, transpathoptions_temp, vfoptions_temp);
+    else
+        [VPath_ii,PolicyPath_ii]=ValueFnOnTransPath_Case1(PricePath_temp, ParamPath_temp, T, V_final_temp, Policy_final_temp, Parameters_temp, n_d_temp, n_a_temp, n_z_temp, d_grid_temp, a_grid_temp,z_grid_temp, pi_z_temp, DiscountFactorParamNames_temp, ReturnFn_temp, transpathoptions_temp, vfoptions_temp);
+    end
     % Note: T cannot depend on ptype, nor can PricePath depend on ptype
 
     if vfoptions_temp.ptypestorecpu==1

From ce0f30f33a41367817ba809351a2355c5bd1a4fb Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Mon, 13 Apr 2026 12:33:11 +1200
Subject: [PATCH 2/4] Revert function signature changes

Revert changes to AgentDistOnTransPath_Case1.m and ValueFnOnTransPath_Case1.m function signatures.
---
 StationaryDist/TransPath/AgentDistOnTransPath_Case1.m | 2 +-
 ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m      | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m b/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m
index cd8501c7..0ac065fd 100644
--- a/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m
+++ b/StationaryDist/TransPath/AgentDistOnTransPath_Case1.m
@@ -1,4 +1,4 @@
-function AgentDistPath=AgentDistOnTransPath_Case1(AgentDist_initial,PricePath,ParamPath,PolicyPath,n_d,n_a,n_z,pi_z,T,Parameters,simoptions)
+function AgentDistPath=AgentDistOnTransPath_Case1(AgentDist_initial, PolicyPath,n_d,n_a,n_z,pi_z,T,simoptions,Parameters,PricePath,ParamPath)
 n_e=0; % NOT YET IMPLEMENTED FOR TRANSITION PATHS
 
 N_d=prod(n_d);
diff --git a/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m b/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m
index 5a756eae..766a563c 100644
--- a/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m
+++ b/ValueFnIter/TransPath/ValueFnOnTransPath_Case1.m
@@ -1,4 +1,4 @@
-function [VPath,PolicyPath]=ValueFnOnTransPath_Case1(PricePath, ParamPath, T, V_final, Policy_final, Parameters, n_d, n_a, n_z, d_grid, a_grid,z_grid, pi_z, DiscountFactorParamNames, ReturnFn, transpathoptions, vfoptions)
+function [VPath,PolicyPath]=ValueFnOnTransPath_Case1(PricePath, ParamPath, T, V_final, Policy_final, Parameters, n_d, n_a, n_z, pi_z, d_grid, a_grid,z_grid, DiscountFactorParamNames, ReturnFn, transpathoptions, vfoptions)
 % transpathoptions, vfoptions and simoptions are optional inputs
 
 %% Check which transpathoptions have been used, set all others to defaults 

From a5474f3d302120e2d78cb8f004c092f102d2e8e1 Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Fri, 8 May 2026 10:24:28 +1200
Subject: [PATCH 3/4] Cleanups for cleaner Pull Request

---
 ...lFnOnTransPath_AggVars_Case1_FHorz_PType.m |   6 +-
 ...teroAgentStationaryEqm_Case1_FHorz_PType.m | 100 ++++-------
 .../AgentDistOnTransPath_Case1_FHorz_PType.m  |  26 ++-
 ...ransitionPath_Case1_FHorz_PType_shooting.m |   2 -
 .../TransitionPath_InfHorz_PType_singlepath.m | 168 ------------------
 .../PricePathParamPath_FHorz_StructToMatrix.m |   1 -
 6 files changed, 50 insertions(+), 253 deletions(-)
 delete mode 100644 TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m

diff --git a/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m b/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
index 9f6ab038..840d0432 100644
--- a/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
+++ b/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
@@ -184,11 +184,7 @@
 
     [FnsToEvaluate_temp,~, ~,~]=PType_FnsToEvaluate(FnsToEvaluate,Names_i,ii,l_d_temp,l_a_temp,l_z_temp,0);
 
-    if isfinite(N_j_temp)
-        AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1_FHorz(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
-    else
-        AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
-    end
+    AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1_FHorz(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
 
     FnNames_temp=fieldnames(FnsToEvaluate_temp);
     for ff=1:length(FnNames_temp)
diff --git a/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m b/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m
index 46fdf133..2fbebcb3 100644
--- a/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m
+++ b/HeterogeneousAgent/FHorz/PType/HeteroAgentStationaryEqm_Case1_FHorz_PType.m
@@ -444,6 +444,7 @@
         end
     end
 
+
     %% Set up exogenous shock grids now (so they can then just be reused every time)
 
     if isstruct(z_grid)
@@ -489,45 +490,27 @@
     end
 
     % Switch over to joint-grids
-    if isfinite(PTypeStructure.(iistr).N_j) % FHorz
-        % If z (and e) are not determined in GE, then compute z_gridvals_J and pi_z_J now (and e_gridvals_J and pi_e_J)
-        if heteroagentoptions.gridsinGE(ii)==0
-            [PTypeStructure.(iistr).z_gridvals_J, PTypeStructure.(iistr).pi_z_J, PTypeStructure.(iistr).vfoptions]=ExogShockSetup_FHorz(PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).z_grid,PTypeStructure.(iistr).pi_z,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).vfoptions,3);
-            % Note: these are actually z_gridvals_J and pi_z_J
-            PTypeStructure.(iistr).simoptions.e_gridvals_J=PTypeStructure.(iistr).vfoptions.e_gridvals_J; % Note, will be [] if no e
-            PTypeStructure.(iistr).simoptions.pi_e_J=PTypeStructure.(iistr).vfoptions.pi_e_J; % Note, will be [] if no e
-        else
-            % % Create placeholders, as these will need to be created in general eqm since they depend on General eqm parameters
-            % PTypeStructure.(iistr).z_gridvals_J=[];
-            % PTypeStructure.(iistr).pi_z_J=[];
-        end
-        PTypeStructure.(iistr)=rmfield(PTypeStructure.(iistr),'z_grid'); % Should not be used, as now have z_gridvals_J
-        PTypeStructure.(iistr)=rmfield(PTypeStructure.(iistr),'pi_z'); % Should not be used, as now have pi_z_J
-        if isfield(PTypeStructure.(iistr).simoptions,'ExogShockFn') % Note: ExogShockSetup_FHorz() removed ExogShockFn from vfoptions but not from simoptions
-            if heteroagentoptions.useCustomModelStats==1
-                heteroagentoptions.CustomModelStatsInputs.z_grid=PTypeStructure.(iistr).z_gridvals_J;
-                heteroagentoptions.CustomModelStatsInputs.pi_z=PTypeStructure.(iistr).pi_z_J;
-            end
-            PTypeStructure.(iistr).simoptions=rmfield(simoptions,'ExogShockFn');
-        end
-    else % InfHorz
-        if heteroagentoptions.gridsinGE(ii)==0
-            [PTypeStructure.(iistr).z_gridvals, PTypeStructure.(iistr).pi_z, PTypeStructure.(iistr).vfoptions]=ExogShockSetup_InfHorz(PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).z_grid,PTypeStructure.(iistr).pi_z,PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).vfoptions,3);
-            PTypeStructure.(iistr).simoptions.e_gridvals=PTypeStructure.(iistr).vfoptions.e_gridvals; % Note, will be [] if no e
-            PTypeStructure.(iistr).simoptions.pi_e=PTypeStructure.(iistr).vfoptions.pi_e; % Note, will be [] if no e
-        else
-            % % Create placeholders, as these will need to be created in general eqm since they depend on General eqm parameters
-            % PTypeStructure.(iistr).z_gridvals=[];
-            % PTypeStructure.(iistr).pi_z=[];
-        end
-        if isfield(PTypeStructure.(iistr).simoptions,'ExogShockFn') % Note: ExogShockSetup_InfHorz() removed ExogShockFn from vfoptions but not from simoptions
-            if heteroagentoptions.useCustomModelStats==1
-                heteroagentoptions.CustomModelStatsInputs.z_grid=PTypeStructure.(iistr).z_gridvals_J;
-                heteroagentoptions.CustomModelStatsInputs.pi_z=PTypeStructure.(iistr).pi_z_J;
-            end
-            PTypeStructure.(iistr).simoptions=rmfield(PTypeStructure.(iistr).simoptions,'ExogShockFn');
+    % If z (and e) are not determined in GE, then compute z_gridvals_J and pi_z_J now (and e_gridvals_J and pi_e_J)
+    if heteroagentoptions.gridsinGE(ii)==0
+        [PTypeStructure.(iistr).z_gridvals_J, PTypeStructure.(iistr).pi_z_J, PTypeStructure.(iistr).vfoptions]=ExogShockSetup_FHorz(PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).z_grid,PTypeStructure.(iistr).pi_z,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).vfoptions,3);
+        % Note: these are actually z_gridvals_J and pi_z_J
+        PTypeStructure.(iistr).simoptions.e_gridvals_J=PTypeStructure.(iistr).vfoptions.e_gridvals_J; % Note, will be [] if no e
+        PTypeStructure.(iistr).simoptions.pi_e_J=PTypeStructure.(iistr).vfoptions.pi_e_J; % Note, will be [] if no e
+    else
+        % % Create placeholders, as these will need to be created in general eqm since they depend on General eqm parameters
+        % PTypeStructure.(iistr).z_gridvals_J=[];
+        % PTypeStructure.(iistr).pi_z_J=[];
+    end
+    PTypeStructure.(iistr)=rmfield(PTypeStructure.(iistr),'z_grid'); % Should not be used, as now have z_gridvals_J
+    PTypeStructure.(iistr)=rmfield(PTypeStructure.(iistr),'pi_z'); % Should not be used, as now have pi_z_J
+    if isfield(PTypeStructure.(iistr).simoptions,'ExogShockFn') % Note: ExogShockSetup_FHorz() removed ExogShockFn from vfoptions but not from simoptions
+        if heteroagentoptions.useCustomModelStats==1
+            heteroagentoptions.CustomModelStatsInputs.z_grid=PTypeStructure.(iistr).z_gridvals_J;
+            heteroagentoptions.CustomModelStatsInputs.pi_z=PTypeStructure.(iistr).pi_z_J;
         end
+        PTypeStructure.(iistr).simoptions=rmfield(simoptions,'ExogShockFn');
     end
+
     % Regardless of whether they are done here of in _subfn, they will be precomputed by the time we get to the value fn, staty dist, etc. So
     PTypeStructure.(iistr).vfoptions.alreadygridvals=1;
     PTypeStructure.(iistr).simoptions.alreadygridvals=1;
@@ -556,13 +539,10 @@
             end
         end
         
-        if isfinite(PTypeStructure.(iistr).N_j) % FHorz
-            PTypeStructure.(iistr).vfoptions=SemiExogShockSetup_FHorz(PTypeStructure.(iistr).n_d,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).d_grid,PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).vfoptions,2,3);
-            PTypeStructure.(iistr).simoptions.semiz_gridvals_J=PTypeStructure.(iistr).vfoptions.semiz_gridvals_J;
-            PTypeStructure.(iistr).simoptions.pi_semiz_J=PTypeStructure.(iistr).vfoptions.pi_semiz_J;
-        else
-            error('Semi-exogenous state not yet implemented for InfHorz')
-        end
+        PTypeStructure.(iistr).vfoptions=SemiExogShockSetup_FHorz(PTypeStructure.(iistr).n_d,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).d_grid,PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).vfoptions,2,3);
+        PTypeStructure.(iistr).simoptions.semiz_gridvals_J=PTypeStructure.(iistr).vfoptions.semiz_gridvals_J;
+        PTypeStructure.(iistr).simoptions.pi_semiz_J=PTypeStructure.(iistr).vfoptions.pi_semiz_J;
+
         % Regardless of whether they are done here of in _subfn, they will be precomputed by the time we get to the value fn, staty dist, etc. So
         PTypeStructure.(iistr).vfoptions.alreadygridvals_semiexo=1;
         PTypeStructure.(iistr).simoptions.alreadygridvals_semiexo=1;        
@@ -584,28 +564,26 @@
     PTypeStructure.(iistr).ReturnFnParamNames=ReturnFnParamNamesFn(PTypeStructure.(iistr).ReturnFn,PTypeStructure.(iistr).n_d,PTypeStructure.(iistr).n_a,PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).N_j,PTypeStructure.(iistr).vfoptions,PTypeStructure.(iistr).Parameters);
     
     %% jequaloneDist and AgeWeightsParamNames
-    if isfinite(PTypeStructure.(iistr).N_j) % FHorz
-        if isstruct(jequaloneDist)
-            if isfield(jequaloneDist,PTypeStructure.Names_i{ii})
-                if isa(jequaloneDist, 'function_handle')
-                    [PTypeStructure.(iistr).jequaloneDist,~,PTypeStructure.(iistr).Parameters]=jequaloneDist_PType(jequaloneDist.(iistr),PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).simoptions,PTypeStructure.(iistr).n_a,PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).N_i,PTypeStructure.(iistr).PTypeDistParamNames,0);
-                else
-                    PTypeStructure.(iistr).jequaloneDist=jequaloneDist.(PTypeStructure.Names_i{ii});
-                end
+    if isstruct(jequaloneDist)
+        if isfield(jequaloneDist,PTypeStructure.Names_i{ii})
+            if isa(jequaloneDist, 'function_handle')
+                [PTypeStructure.(iistr).jequaloneDist,~,PTypeStructure.(iistr).Parameters]=jequaloneDist_PType(jequaloneDist.(iistr),PTypeStructure.(iistr).Parameters,PTypeStructure.(iistr).simoptions,PTypeStructure.(iistr).n_a,PTypeStructure.(iistr).n_z,PTypeStructure.(iistr).N_i,PTypeStructure.(iistr).PTypeDistParamNames,0);
             else
-                error(['You must input jequaloneDist for permanent type ', PTypeStructure.Names_i{ii}, ' \n'])
+                PTypeStructure.(iistr).jequaloneDist=jequaloneDist.(PTypeStructure.Names_i{ii});
             end
         else
-            PTypeStructure.(iistr).jequaloneDist=jequaloneDist;
+            error(['You must input jequaloneDist for permanent type ', PTypeStructure.Names_i{ii}, ' \n'])
         end
+    else
+        PTypeStructure.(iistr).jequaloneDist=jequaloneDist;
+    end
 
-        PTypeStructure.(iistr).AgeWeightParamNames=AgeWeightParamNames;
-        if isstruct(AgeWeightParamNames)
-            if isfield(AgeWeightParamNames,Names_i{ii})
-                PTypeStructure.(iistr).AgeWeightParamNames=AgeWeightParamNames.(Names_i{ii});
-            else
-                error(['You must input AgeWeightParamNames for permanent type ', Names_i{ii}, ' \n'])
-            end
+    PTypeStructure.(iistr).AgeWeightParamNames=AgeWeightParamNames;
+    if isstruct(AgeWeightParamNames)
+        if isfield(AgeWeightParamNames,Names_i{ii})
+            PTypeStructure.(iistr).AgeWeightParamNames=AgeWeightParamNames.(Names_i{ii});
+        else
+            error(['You must input AgeWeightParamNames for permanent type ', Names_i{ii}, ' \n'])
         end
     end
 
diff --git a/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m b/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m
index 37708531..4513e910 100644
--- a/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m
+++ b/StationaryDist/TransPathFHorz/AgentDistOnTransPath_Case1_FHorz_PType.m
@@ -129,17 +129,15 @@
     else
         AgentDist_initial_temp=AgentDist_initial; % NEED TO DEAL WITH THIS PROPERLY
     end
-    if isfinite(N_j_temp)
-        if isstruct(AgeWeightsParamNames)
-            AgeWeightsParamNames_temp=AgeWeightsParamNames.(Names_i{ii});
-        else
-            AgeWeightsParamNames_temp=AgeWeightsParamNames;
-        end
-        if isstruct(jequalOneDist)
-            jequalOneDist_temp=jequalOneDist.(Names_i{ii});
-        else
-            jequalOneDist_temp=jequalOneDist;
-        end
+    if isstruct(AgeWeightsParamNames)
+        AgeWeightsParamNames_temp=AgeWeightsParamNames.(Names_i{ii});
+    else
+        AgeWeightsParamNames_temp=AgeWeightsParamNames;
+    end
+    if isstruct(jequalOneDist)
+        jequalOneDist_temp=jequalOneDist.(Names_i{ii});
+    else
+        jequalOneDist_temp=jequalOneDist;
     end
     
 
@@ -179,11 +177,7 @@
     
     
     % Compute the agent distribution path for permanent type ii
-    if isfinite(N_j_temp)
-        AgentDistPath_ii=AgentDistOnTransPath_Case1_FHorz(AgentDist_initial_temp, jequalOneDist_temp, PricePath_temp, ParamPath_temp, PolicyPath_temp, AgeWeightsParamNames_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,pi_z_temp, T,Parameters_temp, transpathoptions_temp, simoptions_temp);
-    else
-        AgentDistPath_ii=AgentDistOnTransPath_Case1(AgentDist_initial_temp, PricePath_temp, ParamPath_temp, PolicyPath_temp, n_d_temp,n_a_temp,n_z_temp,pi_z_temp,T,Parameters_temp,simoptions_temp);
-    end
+    AgentDistPath_ii=AgentDistOnTransPath_Case1_FHorz(AgentDist_initial_temp, jequalOneDist_temp, PricePath_temp, ParamPath_temp, PolicyPath_temp, AgeWeightsParamNames_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,pi_z_temp, T,Parameters_temp, transpathoptions_temp, simoptions_temp);
     % Note: T cannot depend on ptype, nor can PricePath depend on ptype
 
     AgentDistPath.(Names_i{ii})=AgentDistPath_ii;
diff --git a/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m b/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m
index c02e68b6..870fdd7e 100644
--- a/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m
+++ b/TransitionPaths/FHorz/PType/TransitionPath_Case1_FHorz_PType_shooting.m
@@ -129,7 +129,6 @@
     AggVarsFullPath=zeros(PTypeStructure.numFnsToEvaluate,T-1,N_i); % Does not include period T
     for ii=1:N_i
         iistr=PTypeStructure.Names_i{ii};
-        
 
         % Following few lines I would normally do outside of the while loop, but have to set them for each ptype
         % AgentDist=AgentDist_initial.(iistr);
@@ -142,7 +141,6 @@
         % Get just the values that correspond to the current ptype
         PricePathOld_ii=PricePathOld(:,PTypeStructure.(iistr).RelevantPricePath);
         ParamPath_ii=ParamPath(:,PTypeStructure.(iistr).RelevantParamPath);
-
         % Have not yet set up the following to allow dependence on ptype (should do this at some point)
         PricePathNames_ii=PricePathNames;
         ParamPathNames_ii=ParamPathNames;
diff --git a/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m b/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m
deleted file mode 100644
index 647fcac6..00000000
--- a/TransitionPaths/InfHorz/PType/TransitionPath_InfHorz_PType_singlepath.m
+++ /dev/null
@@ -1,168 +0,0 @@
-function AggVarsPath=TransitionPath_InfHorz_PType_singlepath(PricePathOld, ParamPath, PricePathNames, ParamPathNames, T, V_final, AgentDist_initial, ...
-    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, PricePathSizeVec, ParamPathSizeVec, ...
-    use_tminus1price, use_tminus1params, use_tplus1price, use_tminus1AggVars, tminus1priceNames, tminus1paramNames, tplus1priceNames, tminus1AggVarsNames, ...
-    transpathoptions, vfoptions, simoptions)
-% PricePathOld is matrix of size T-by-'number of prices'
-% ParamPath is matrix of size T-by-'number of parameters that change over path'
-
-% Remark to self: No real need for T as input, as this is anyway the length of PricePathOld
-
-% For this agent type, first go back through the value & policy fns.
-% Then forwards through agent dist and agg vars.
-
-l_a=length(n_a);
-PolicyIndexesPath=zeros(l_d+l_a,N_a,N_z,T-1,'gpuArray'); %Periods 1 to T-1
-
-% This and much of the rest of this code borrowed from TransitionPath_InfHorz_shooting.m
-if simoptions.gridinterplayer==0
-    II1=(1:1:N_a*N_z); % Index for this period (a,z)
-    IIones=ones(N_a*N_z,1); % Next period 'probabilities'
-elseif simoptions.gridinterplayer==1
-    PolicyProbsPath=zeros(N_a*N_z,2,T-1,'gpuArray'); % preallocate
-    II2=([1:1:N_a*N_z; 1:1:N_a*N_z]'); % Index for this period (a,z), note the 2 copies
-end
-
-if size(V_final)~=[N_a,N_z]
-    error("V_final wrong shape")
-    V_final=reshape(V_final,[N_a,N_z]);
-end
-
-%First, go from T-1 to 1 calculating the Value function and Optimal
-%policy function at each step. Since we won't need to keep the value
-%functions for anything later we just store the next period one in
-%Vnext, and the current period one to be calculated in V
-Vnext=V_final;
-for ttr=1:T-1 %so t=T-i
-    % The following digs deeper into PricePathOld and ParamPath in
-    % FHorz case--check it
-    for kk=1:length(PricePathNames)
-        Parameters.(PricePathNames{kk})=PricePathOld(T-ttr,kk);
-    end
-    for kk=1:length(ParamPathNames)
-        Parameters.(ParamPathNames{kk})=ParamPath(T-ttr,kk);
-    end
-    
-    [V, Policy]=ValueFnIter_InfHorz_TPath_SingleStep(Vnext,n_d,n_a,n_z,d_grid, a_grid, z_grid, pi_z, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
-    % The VKron input is next period value fn, the VKron output is this period.
-    % Policy is kept in the form where it is just a single-value in (d,a')
-
-    PolicyIndexesPath(:,:,:,T-ttr)=Policy;
-    Vnext=V;
-    
-end
-% Free up space on GPU by deleting things no longer needed
-clear V Vnext
-
-%% Modify PolicyIndexesPath into forms needed for forward iteration
-% Create version of PolicyIndexesPath in form we want for the agent distribution iteration
-% Creates PolicyaprimezPath, and when using grid interpolation layer also PolicyProbsPath 
-if isscalar(n_a)
-    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:),[N_a*N_z,T-1]); % aprime index
-elseif length(n_a)==2
-    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:)+n_a(1)*(PolicyIndexesPath(l_d+2,:,:,:)-1),[N_a*N_z,T-1]);
-elseif length(n_a)==3
-    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:)+n_a(1)*(PolicyIndexesPath(l_d+2,:,:,:)-1)+n_a(1)*n_a(2)*(PolicyIndexesPath(l_d+3,:,:,:)-1),[N_a*N_z,T-1]);
-elseif length(n_a)==4
-    PolicyaprimePath=reshape(PolicyIndexesPath(l_d+1,:,:,:)+n_a(1)*(PolicyIndexesPath(l_d+2,:,:,:)-1)+n_a(1)*n_a(2)*(PolicyIndexesPath(l_d+3,:,:,:)-1)+n_a(1)*n_a(2)*n_a(3)*(PolicyIndexesPath(l_d+4,:,:,:)-1),[N_a*N_z,T-1]);
-end
-PolicyaprimezPath=PolicyaprimePath+repelem(N_a*gpuArray(0:1:N_z-1)',N_a,1);
-if simoptions.gridinterplayer==1
-    PolicyaprimezPath=reshape(PolicyaprimezPath,[N_a*N_z,1,T-1]); % reinterpret this as lower grid index
-    PolicyaprimezPath=repelem(PolicyaprimezPath,1,2,1); % create copy that will be the upper grid index
-    PolicyaprimezPath(:,2,:)=PolicyaprimezPath(:,2,:)+1; % upper grid index
-    PolicyProbsPath(:,2,:)=reshape(PolicyIndexesPath(l_d+l_aprime+1,:,:),[N_a*N_z,1,T-1]); % L2 index
-    PolicyProbsPath(:,2,:)=(PolicyProbsPath(:,2,:)-1)/(1+simoptions.ngridinterp); % probability of upper grid point
-    PolicyProbsPath(:,1,:)=1-PolicyProbsPath(:,2,:); % probability of lower grid point
-end
-% Create PolicyValuesPath from PolicyIndexesPath for use in calculating model stats
-PolicyValuesPath=PolicyInd2Val_InfHorz_TPath(PolicyIndexesPath,n_d,n_a,n_z,T-1,d_grid,a_grid,vfoptions,1);
-PolicyValuesPath=permute(reshape(PolicyValuesPath,[size(PolicyValuesPath,1),N_a,N_z,T-1]),[2,3,1,4]); %[N_a,N_z,l_d+l_a,T-1]
-
-%Now we have the full PolicyIndexesPath, we go forward in time from 1
-%to T using the policies to update the agents distribution generating a
-%new price path
-%Call AgentDist the current periods distn
-
-AgentDist=sparse(gather(reshape(AgentDist_initial,[N_a*N_z,1])));
-AggVarsPath=zeros(length(FnsToEvaluate),T-1);
-pi_z_sparse=sparse(gather(pi_z)); % Need full pi_z for value fn, and sparse for agent dist
-
-for tt=1:T-1
-    %% Setup the Parameters for period tt
-
-    % Get t-1 PricePath and ParamPath before we update them
-    if use_tminus1price==1
-        for pp=1:length(tminus1priceNames)
-            if tt>1
-                Parameters.([tminus1priceNames{pp},'_tminus1'])=Parameters.(tminus1priceNames{pp});
-            else
-                Parameters.([tminus1priceNames{pp},'_tminus1'])=transpathoptions.initialvalues.(tminus1priceNames{pp});
-            end
-        end
-    end
-    if use_tminus1params==1
-        for pp=1:length(tminus1paramNames)
-            if tt>1
-                Parameters.([tminus1paramNames{pp},'_tminus1'])=Parameters.(tminus1paramNames{pp});
-            else
-                Parameters.([tminus1paramNames{pp},'_tminus1'])=transpathoptions.initialvalues.(tminus1paramNames{pp});
-            end
-        end
-    end
-    % Get t-1 AggVars before we update them
-    if use_tminus1AggVars==1
-        for pp=1:length(tminus1AggVarsNames)
-            if tt>1
-                % The AggVars have not yet been updated, so they still contain previous period values
-                Parameters.([tminus1AggVarsNames{pp},'_tminus1'])=Parameters.(tminus1AggVarsNames{pp});
-            else
-                Parameters.([tminus1AggVarsNames{pp},'_tminus1'])=transpathoptions.initialvalues.(tminus1AggVarsNames{pp});
-            end
-        end
-    end
-
-    % Update current PricePath and ParamPath
-    for kk=1:length(PricePathNames)
-        Parameters.(PricePathNames{kk})=PricePathOld(tt,PricePathSizeVec(1,kk):PricePathSizeVec(2,kk));
-    end
-    for kk=1:length(ParamPathNames)
-        Parameters.(ParamPathNames{kk})=ParamPath(tt,ParamPathSizeVec(1,kk):ParamPathSizeVec(2,kk));
-    end
-
-    % Get t+1 PricePath
-    if use_tplus1price==1
-        for pp=1:length(tplus1priceNames)
-            kk=tplus1pricePathkk(pp);
-            Parameters.([tplus1priceNames{pp},'_tplus1'])=PricePathOld(tt+1,PricePathSizeVec(1,kk):PricePathSizeVec(2,kk)); % Make is so that the time t+1 variables can be used
-        end
-    end
-    
-    %% Get the current optimal policy, and iterate the agent dist
-    if simoptions.gridinterplayer==0
-        AgentDistnext=AgentDist_InfHorz_TPath_SingleStep(AgentDist,PolicyaprimezPath(:,tt),II1,IIones,N_a,N_z,pi_z_sparse);
-    elseif simoptions.gridinterplayer==1
-        AgentDistnext=AgentDist_InfHorz_TPath_SingleStep_nProbs(AgentDist,PolicyaprimezPath(:,:,tt),II2,PolicyProbsPath(:,:,tt),N_a,N_z,pi_z_sparse);
-    end
-
-    %% AggVars
-    AggVars_Means=EvalFnOnAgentDist_InfHorz_TPath_SingleStep_AggVars(gpuArray(full(AgentDist)), PolicyValuesPath(:,:,:,tt), FnsToEvaluate, Parameters, FnsToEvaluateParamNames, AggVarNames, n_a, n_z, a_gridvals, z_grid,1);
-    AggVars=zeros(length(AggVars_Means),1);
-    if length(fieldnames(AggVars_Means))~=length(AggVarNames)
-        error(["AggVar length disparity:";"---------";AggVarNames;"--- vs ---";fieldnames(AggVars_Means)]);
-    end
-    for ii=1:length(AggVarNames)
-        AggVars(ii)=AggVars_Means.(AggVarNames{ii}).Mean;
-        Parameters.(AggVarNames{ii})=AggVars(ii);
-    end
-
-    % Do nothing with IntermediateEqns and GeneralEqmEqns as they are outside PType scope (they are where all the PTypes meet).
-
-    AgentDist=AgentDistnext;
-
-    AggVarsPath(:,tt)=AggVars;
-end
-
-
-end
diff --git a/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m b/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m
index bb637d0a..723a6f3a 100644
--- a/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m
+++ b/TransitionPaths/Subcodes/PricePathParamPath_FHorz_StructToMatrix.m
@@ -71,7 +71,6 @@
     % Note: Internally PricePathOld is matrix of size T-by-'number of prices'.
     % ParamPath is matrix of size T-by-'number of parameters that change over the transition path'.
     % Actually, some of those prices are 1-by-N_j or N_i or both, so is more subtle than this.
-    Names_i=fieldnames(N_j);
     PricePathNames=fieldnames(PricePathStruct);
     PricePathSizeVec=zeros(1,length(PricePathNames)); % Allows for a given price param to depend on age (or permanent type)
     if isstruct(N_j)

From 9333f543a8d07ca5a944cbaca7286460a4f6854f Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Fri, 8 May 2026 10:35:05 +1200
Subject: [PATCH 4/4] More cleaning

---
 ...FnOnAgentDist_AllStats_FHorz_Case1_PType.m | 75 ++++++-------------
 ...lFnOnTransPath_AggVars_Case1_FHorz_PType.m |  1 +
 2 files changed, 24 insertions(+), 52 deletions(-)

diff --git a/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m b/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m
index 31c4b14f..a738be35 100644
--- a/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m
+++ b/EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m
@@ -335,35 +335,18 @@
     N_a_temp=prod(n_a_temp);
     
     a_gridvals_temp=CreateGridvals(n_a_temp,a_grid_temp,1);
-    if isfinite(N_j_temp)
-        % Turn (semiz,z,e) into z_gridvals_J_temp as FnsToEvalute do not distinguish them
-        [n_z_temp,z_gridvals_J_temp,N_z_temp,l_z_temp,simoptions_temp]=CreateGridvals_FnsToEvaluate_FHorz(n_z_temp,z_grid_temp,N_j_temp,simoptions_temp,Parameters_temp);
-    else
-        l_z_temp=length(n_z_temp);
-        N_z_temp=prod(n_z_temp);
-        z_gridvals_temp=CreateGridvals(n_z_temp,z_grid_temp,1);
-    end
+    % Turn (semiz,z,e) into z_gridvals_J_temp as FnsToEvalute do not distinguish them
+    [n_z_temp,z_gridvals_J_temp,N_z_temp,l_z_temp,simoptions_temp]=CreateGridvals_FnsToEvaluate_FHorz(n_z_temp,z_grid_temp,N_j_temp,simoptions_temp,Parameters_temp);
     if N_z_temp==0
         N_z_temp=1; % Just makes things easier below
     end
     
     % Switch to PolicyVals
-    if isfinite(N_j_temp)
-        PolicyValues_temp=PolicyInd2Val_FHorz(PolicyIndexes_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,d_grid_temp,a_grid_temp,simoptions_temp,1);
-        if l_z_temp==0
-            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,1]); % (N_a,N_j,l_daprime)
-        else
-            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,4,1]); % (N_a,N_z,N_j,l_daprime)
-        end
-        StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
+    PolicyValues_temp=PolicyInd2Val_FHorz(PolicyIndexes_temp,n_d_temp,n_a_temp,n_z_temp,N_j_temp,d_grid_temp,a_grid_temp,simoptions_temp,1);
+    if l_z_temp==0
+        PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,1]); % (N_a,N_j,l_daprime)
     else
-        PolicyValues_temp=PolicyInd2Val_Case1(PolicyIndexes_temp,n_d_temp,n_a_temp,n_z_temp,d_grid_temp,a_grid_temp,simoptions_temp,1);
-        if l_z_temp==0
-            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,1]); % (N_a,l_daprime)
-        else
-            PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,1]); % (N_a,N_z,l_daprime)
-        end
-        StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
+        PolicyValuesPermute_temp=permute(PolicyValues_temp,[2,3,4,1]); % (N_a,N_z,N_j,l_daprime)
     end
     l_daprime_temp=size(PolicyValues_temp,1);
 
@@ -371,6 +354,7 @@
     FnsAndPTypeIndicator(:,ii)=FnsAndPTypeIndicator_ii;
 
     %% Some things that don't need to go in the loop over FnsToEvalaute
+    StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
     % Eliminate all the zero-weighted points (this doesn't really save runtime for the exact calculation and often can increase it, but
     % for the createDigest it slashes the runtime. So since we want it then we may as well do it now.)    
     temp=logical(StationaryDist_ii~=0);
@@ -392,21 +376,14 @@
                 CondlRestnFnParamNames={};
             end
 
-            if isfinite(N_j_temp)
-                if l_z_temp==0
-                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,2); % j in 2nd dimension: (a,j,l_d+l_a), so we want j to be after N_a
-                    RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,0,a_gridvals_temp,[]));
-                else
-                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,3); % j in 3rd dimension: (a,z,j,l_d+l_a), so we want j to be after N_a and N_z
-                    RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp));
-                end
-                RestrictionValues=reshape(RestrictionValues,[N_a_temp*N_z_temp*N_j_temp,1]);
+            if l_z_temp==0
+                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,2); % j in 2nd dimension: (a,j,l_d+l_a), so we want j to be after N_a
+                RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,0,a_gridvals_temp,[]));
             else
-                CondlRestnFnParamsCell=CreateCellFromParams(Parameters_temp,CondlRestnFnParamNames);
-
-                RestrictionValues=logical(EvalFnOnAgentDist_Grid(CondlRestnFn, CondlRestnFnParamsCell,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_temp));
-                RestrictionValues=reshape(RestrictionValues,[N_a_temp*N_z_temp,1]);
+                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,CondlRestnFnParamNames,N_j_temp,3); % j in 3rd dimension: (a,z,j,l_d+l_a), so we want j to be after N_a and N_z
+                RestrictionValues=logical(EvalFnOnAgentDist_Grid_J(CondlRestnFn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp));
             end
+            RestrictionValues=reshape(RestrictionValues,[N_a_temp*N_z_temp*N_j_temp,1]);
 
             RestrictedStationaryDistVec=StationaryDist_ii;
             RestrictedStationaryDistVec(~RestrictionValues)=0; % zero mass on all points that do not meet the restriction
@@ -446,24 +423,18 @@
             else
                 FnsToEvaluateParamNames={};
             end
-            if isfinite(N_j_temp)
-                if l_z_temp==0
-                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,2);
-                else
-                    CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,3);
-                end
-                
-                %% We have set up the current PType, now do some calculations for it.
-                simoptions_temp.keepoutputasmatrix=1;
-                ValuesOnGrid_ii=EvalFnOnAgentDist_Grid_J(tempfn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp);
-                ValuesOnGrid_ii=reshape(ValuesOnGrid_ii,[N_a_temp*N_z_temp*N_j_temp,1]);
-    
-                % StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
+            if l_z_temp==0
+                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,2);
             else
-                FnToEvaluateParamsCell=CreateCellFromParams(Parameters_temp,FnsToEvaluateParamNames);
-                ValuesOnGrid_ii=EvalFnOnAgentDist_Grid(tempfn, FnToEvaluateParamsCell,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_temp);
-                ValuesOnGrid_ii=reshape(ValuesOnGrid_ii,[N_a_temp*N_z_temp,1]);
+                CellOverAgeOfParamValues=CreateCellOverAgeFromParams(Parameters_temp,FnsToEvaluateParamNames,N_j_temp,3);
             end
+            
+            %% We have set up the current PType, now do some calculations for it.
+            simoptions_temp.keepoutputasmatrix=1;
+            ValuesOnGrid_ii=EvalFnOnAgentDist_Grid_J(tempfn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp);
+            ValuesOnGrid_ii=reshape(ValuesOnGrid_ii,[N_a_temp*N_z_temp*N_j_temp,1]);
+
+            % StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp*N_j_temp,1]); % Note: does not impose *StationaryDist.ptweights(ii)
 
             % Eliminate all the zero-weighted points (this doesn't really save runtime for the exact calculation and often can increase it, but
             % for the createDigest it slashes the runtime. So since we want it then we may as well do it now.)
diff --git a/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m b/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
index 840d0432..73201f17 100644
--- a/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
+++ b/EvaluateFnOnAgentDist/TransPathFHorz/PType/EvalFnOnTransPath_AggVars_Case1_FHorz_PType.m
@@ -186,6 +186,7 @@
 
     AggVarsPath_ii=EvalFnOnTransPath_AggVars_Case1_FHorz(FnsToEvaluate_temp, AgentDistPath_temp, PolicyPath_temp, PricePath_temp, ParamPath_temp, Parameters_temp, T, n_d_temp, n_a_temp, n_z_temp, N_j_temp, d_grid_temp, a_grid_temp,z_grid_temp, transpathoptions_temp, simoptions_temp);
 
+
     FnNames_temp=fieldnames(FnsToEvaluate_temp);
     for ff=1:length(FnNames_temp)
         % Keep the ptype-conditional values