Remaining PType and TPath changes (no ExpAsset)

EvaluateFnOnAgentDist/FHorz/LifeCycleProfiles_FHorz_Case1.m

@@ -467,7 +467,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]);

EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AllStats_FHorz_Case1_PType.m

@@ -335,26 +335,42 @@
     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);
 
     [~,~,~,FnsAndPTypeIndicator_ii]=PType_FnsToEvaluate(FnsToEvaluate,Names_i,ii,l_d_temp,l_a_temp,l_z_temp,0);
     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);
@@ -381,14 +397,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
@@ -416,24 +439,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.)

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;

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)

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

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;
@@ -443,7 +444,6 @@
         end
     end
 
-
     %% Set up exogenous shock grids now (so they can then just be reused every time)
 
     if isstruct(z_grid)
@@ -571,7 +571,6 @@
 
     %% 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')
@@ -580,9 +579,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;
@@ -593,9 +590,7 @@
             if isfield(AgeWeightParamNames,Names_i{ii})
                 PTypeStructure.(iistr).AgeWeightParamNames=AgeWeightParamNames.(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

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')
@@ -18,7 +43,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
 
@@ -78,12 +103,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);
@@ -168,8 +193,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')

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')