Aggvars in ptype parameters

DiscretizationMethods/discretizeAR1_FarmerToda/discretizeAR1_FarmerToda.m

@@ -33,7 +33,7 @@
 %%%%%%%%%%%%%%%
 % Original paper:
 % Farmer & Toda (2017) - Discretizing Nonlinear, Non-Gaussian Markov Processes with Exact Conditional Moments
-% They how that this method outperforms both Tauchen and Rouwenhorst for almost all discretization of Gaussian AR(1).
+% They show that this method outperforms both Tauchen and Rouwenhorst for almost all discretization of Gaussian AR(1).
 
 if rho>=0.99
     fprintf('COMMENT: When discretizing gaussian AR(1) process with autocorrelation (rho) greater than 0.99 (which you currently have), the Rouwenhorst method tends to outperform Farmer-Toda method. \n')

EvaluateFnOnAgentDist/AggShocksPath/EvalFnOnAggShocksPath_AllStats_InfHorz.m

@@ -0,0 +1,103 @@
+function AllStatsPath=EvalFnOnAggShocksPath_AllStats_InfHorz(GeneralizedTransitionFn,T,n_d,n_a,n_z,n_S,d_grid,a_grid,z_grid,FnsToEvaluate,Parameters,simoptions)
+
+if ~exist('simoptions','var')
+    simoptions.npoints=100;
+    simoptions.nquantiles=20;
+    simoptions.whichstats=ones(7,1); % See StatsFromWeightedGrid(), zeros skip some stats and can be used to reduce runtimes 
+else
+    if ~isfield(simoptions,'npoints')
+        simoptions.npoints=100;
+    end
+    if ~isfield(simoptions,'nquantiles')
+        simoptions.nquantiles=20;
+    end
+    if ~isfield(simoptions,'whichstats')
+        simoptions.whichstats=ones(7,1); % See StatsFromWeightedGrid(), zeros skip some stats and can be used to reduce runtimes 
+    end
+end
+
+%% Setup
+PricePathNames=fieldnames(GeneralizedTransitionFn.PricePath);
+AggShockNames=fieldnames(GeneralizedTransitionFn.AggShocksPath);
+
+FnsToEvalNames=fieldnames(FnsToEvaluate);
+for ff=1:length(FnsToEvalNames)
+    AllStatsPath.(FnsToEvalNames{ff}).Mean=zeros(1,T);
+    AllStatsPath.(FnsToEvalNames{ff}).Median=zeros(1,T);
+    AllStatsPath.(FnsToEvalNames{ff}).RatioMeanToMedian=zeros(1,T);
+    AllStatsPath.(FnsToEvalNames{ff}).Variance=zeros(1,T);
+    AllStatsPath.(FnsToEvalNames{ff}).StdDeviation=zeros(1,T);
+    AllStatsPath.(FnsToEvalNames{ff}).LorenzCurve=zeros(simoptions.npoints,T);
+    AllStatsPath.(FnsToEvalNames{ff}).Gini=zeros(1,T);
+    AllStatsPath.(FnsToEvalNames{ff}).QuantileCutoffs=zeros(simoptions.nquantiles+1,T);  % Includes the min and max values
+    AllStatsPath.(FnsToEvalNames{ff}).QuantileMeans=zeros(simoptions.nquantiles,T);
+end
+
+
+%% Loop over tt=1:T and get AllStats
+% We already have GeneralizedTransitionFn.AgentDistPath so just calculate from this.
+for tt=1:T
+    for pp=1:length(PricePathNames)
+        temp=GeneralizedTransitionFn.PricePath.(PricePathNames{pp});
+        Parameters.(PricePathNames{pp})=temp(tt);
+    end
+    for SS_c=1:length(n_S)
+        temp=GeneralizedTransitionFn.AggShocksPath.(AggShockNames{SS_c});
+        Parameters.(AggShockNames{SS_c})=temp(tt);
+    end
+
+    AgentDist=GeneralizedTransitionFn.AgentDistPath(:,:,tt);
+    Policy=GeneralizedTransitionFn.PolicyPath(:,:,:,tt);
+
+    % would be slightly faster if I pass FnsToEvaluateCell and FnsToEvaluateParamNames
+    AllStats_tt=EvalFnOnAgentDist_AllStats_Case1(AgentDist,Policy,FnsToEvaluate,Parameters,[],n_d,n_a,n_z,d_grid,a_grid,z_grid,simoptions);
+
+    for ff=1:length(FnsToEvalNames)
+        temp=AllStatsPath.(FnsToEvalNames{ff}).Mean;
+        temp(tt)=AllStats_tt.(FnsToEvalNames{ff}).Mean;
+        AllStatsPath.(FnsToEvalNames{ff}).Mean=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).Median;
+        temp(tt)=AllStats_tt.(FnsToEvalNames{ff}).Median;
+        AllStatsPath.(FnsToEvalNames{ff}).Median=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).RatioMeanToMedian;
+        temp(tt)=AllStats_tt.(FnsToEvalNames{ff}).RatioMeanToMedian;
+        AllStatsPath.(FnsToEvalNames{ff}).RatioMeanToMedian=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).Variance;
+        temp(tt)=AllStats_tt.(FnsToEvalNames{ff}).Variance;
+        AllStatsPath.(FnsToEvalNames{ff}).Variance=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).StdDeviation;
+        temp(tt)=AllStats_tt.(FnsToEvalNames{ff}).StdDeviation;
+        AllStatsPath.(FnsToEvalNames{ff}).StdDeviation=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).LorenzCurve;
+        temp(:,tt)=AllStats_tt.(FnsToEvalNames{ff}).LorenzCurve;
+        AllStatsPath.(FnsToEvalNames{ff}).LorenzCurve=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).Gini;
+        temp(tt)=AllStats_tt.(FnsToEvalNames{ff}).Gini;
+        AllStatsPath.(FnsToEvalNames{ff}).Gini=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).QuantileCutoffs;
+        temp(:,tt)=AllStats_tt.(FnsToEvalNames{ff}).QuantileCutoffs;
+        AllStatsPath.(FnsToEvalNames{ff}).QuantileCutoffs=temp;
+
+        temp=AllStatsPath.(FnsToEvalNames{ff}).QuantileMeans;
+        temp(:,tt)=AllStats_tt.(FnsToEvalNames{ff}).QuantileMeans;
+        AllStatsPath.(FnsToEvalNames{ff}).QuantileMeans=temp;
+    end
+
+end
+
+
+
+
+
+
+
+
+
+end

EvaluateFnOnAgentDist/EvalFnOnAgentDist_AggVars_Case1.m

@@ -77,7 +77,7 @@
 if isfield(simoptions,'outputasstructure')
     if simoptions.outputasstructure==1
         FnsToEvaluateStruct=1;
-        FnsToEvalNames=simoptions.AggVarNames;
+        % FnsToEvalNames=simoptions.AggVarNames;
     elseif simoptions.outputasstructure==0
         FnsToEvaluateStruct=0;
     end

EvaluateFnOnAgentDist/EvalFnOnAgentDist_AllStats_Case1.m

@@ -149,7 +149,7 @@
     Values=reshape(Values,[N_a*N_z,1]);
 
     AllStats.(FnsToEvalNames{ff})=StatsFromWeightedGrid(Values,StationaryDistVec,simoptions.npoints,simoptions.nquantiles,simoptions.tolerance,0,simoptions.whichstats);
-
+    
     %% If there are any conditional restrictions then deal with these
     % Evaluate AllStats, but conditional on the restriction being one.
     if useCondlRest==1

EvaluateFnOnAgentDist/EvalFnOnAgentDist_ValuesOnGrid_Case1.m

@@ -26,9 +26,13 @@
 if simoptions.gridinterplayer==1
     l_daprime=l_daprime-1;
 end
-a_gridvals=CreateGridvals(n_a,a_grid,1);
-[z_gridvals,~,simoptions]=ExogShockSetup(n_z,z_grid,[],Parameters,simoptions,1);
-
+if Parallel==2
+    a_gridvals=CreateGridvals(n_a,a_grid,1);
+    [z_gridvals,~,simoptions]=ExogShockSetup(n_z,z_grid,[],Parameters,simoptions,1);
+elseif Parallel==1
+    a_gridvals=CreateGridvals(n_a,a_grid,2);
+    z_gridvals=CreateGridvals(n_z,z_grid,2); % CPU, so must just be simple stacked column for z
+end
 
 %% Implement new way of handling FnsToEvaluate
 if isstruct(FnsToEvaluate)
@@ -76,9 +80,12 @@
         Values=reshape(Values,[N_a*N_z,1]);
         ValuesOnGrid.(AggVarNames{ff})=reshape(Values,[n_a,n_z]);
     end
-else
+elseif Parallel==1
+    % CPU
+    simoptions.experienceasset=0; % needs to be set so can use CreateGridvals_Policy()
+    simoptions.experienceassetu=0; % needs to be set so can use CreateGridvals_Policy()
     [d_gridvals, aprime_gridvals]=CreateGridvals_Policy(Policy,n_d,n_a,n_a,n_z,d_grid,a_grid,simoptions,1, 2);
-        
+
     if l_d>0
 
         for ff=1:length(FnsToEvaluate)

EvaluateFnOnAgentDist/FHorz/EvalFnOnAgentDist_ValuesOnGrid_FHorz_Case1.m

@@ -32,7 +32,11 @@
 
 %% Implement new way of handling FnsToEvaluate
 % Figure out l_daprime from Policy
-l_daprime=size(Policy,1);
+if simoptions.gridinterplayer==0
+    l_daprime=size(Policy,1);
+else
+    l_daprime=size(Policy,1)-1;
+end
 
 % Note: l_z includes e and semiz (when appropriate)
 if isstruct(FnsToEvaluate)

EvaluateFnOnAgentDist/FHorz/PType/EvalFnOnAgentDist_AggVars_FHorz_Case1_PType.m

@@ -12,7 +12,7 @@
 % AgeWeightParamNames is either same for all permanent types, or must be passed as a structure.
 %
 % The stationary distribution be a structure and will contain both the
-% weights/distribution across the permenant types, as well as a pdf for the
+% weights/distribution across the permanent types, as well as a pdf for the
 % stationary distribution of each specific permanent type.
 %
 % How exactly to handle these differences between permanent (fixed) types

EvaluateFnOnAgentDist/FHorz/PType/LifeCycleProfiles_FHorz_Case1_PType.m

@@ -11,13 +11,13 @@
 % jequaloneDist can either be same for all permanent types, or must be passed as a structure.
 % AgeWeightParamNames is either same for all permanent types, or must be passed as a structure.
 %
-% The stationary distribution be a structure and will contain both the
+% The stationary distribution must be a structure and will contain both the
 % weights/distribution across the permenant types, as well as a pdf for the
 % stationary distribution of each specific permanent type.
 %
 % How exactly to handle these differences between permanent (fixed) types
 % is to some extent left to the user. You can, for example, input
-% parameters that differ by permanent type as a vector with different rows f
+% parameters that differ by permanent type as a vector with different rows
 % for each type, or as a structure with different fields for each type.
 %
 % Any input that does not depend on the permanent type is just passed in
@@ -176,7 +176,12 @@
             ii=ii+1;
             Names_i{ii}=Names_i2{ii2};
         else % tell the user about it
-            fprintf(['LifeCycleProfiles_FHorz_Case1_PType: Ignoring the ', num2str(ii2), '-th PType, ',Names_i2{ii2}, ' because it is infinite horizon \n']);
+            suffix='th';
+            if ii2<4 % Good up to 20 PTypes
+                suffix_cells={'st', 'nd', 'rd'};
+                suffix=suffix_cells{ii2};
+            end
+            fprintf('LifeCycleProfiles_FHorz_Case1_PType: Ignoring the %d%s PType { %s } because it is infinite horizon \n', ii2, suffix, Names_i2{ii2});
         end
     end   
     % Eliminate any no longer relevant functions from FnsToEvaluate (those which are only used for infinite horizon)
@@ -228,7 +233,7 @@
             end
         end
     end
-elseif length(simoptions.agejshifter)==1 % not using agejshifter
+elseif isscalar(simoptions.agejshifter) % not using agejshifter
     simoptions.agejshifter=zeros(N_i,1);
     N_j_max2=N_j_max;
 else % have inputed as a vector
@@ -563,7 +568,12 @@
             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}));
+                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
@@ -577,7 +587,7 @@
 
                 %% We have set up the current PType, now do some calculations for it.
                 simoptions_temp.keepoutputasmatrix=2;
-                ValuesOnGrid_ffii=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_ffii=EvalFnOnAgentDist_Grid_J(tempfn,CellOverAgeOfParamValues,PolicyValuesPermute_temp,l_daprime_temp,n_a_temp,n_z_temp,a_gridvals_temp,z_gridvals_J_temp);
 
                 ValuesOnGrid_ffii=reshape(ValuesOnGrid_ffii,[N_a_temp*N_z_temp,N_j_temp]);
                 % StationaryDist_ii=reshape(StationaryDist.(Names_i{ii}),[N_a_temp*N_z_temp,N_j_temp]); % Note: does not impose *StationaryDist.ptweights(ii)
@@ -1193,37 +1203,37 @@
             % a structure is there a need to take just a specific part and send
             % only that to the 'non-PType' version of the command.
 
-            if isa(n_d,'struct')
+            if isstruct(n_d)
                 n_d_temp=n_d.(Names_i{ii});
             else
                 n_d_temp=n_d;
             end
-            if isa(n_a,'struct')
+            if isstruct(n_a)
                 n_a_temp=n_a.(Names_i{ii});
             else
                 n_a_temp=n_a;
             end
-            if isa(n_z,'struct')
+            if isstruct(n_z)
                 n_z_temp=n_z.(Names_i{ii});
             else
                 n_z_temp=n_z;
             end
-            if isa(N_j,'struct')
+            if isstruct(N_j)
                 N_j_temp=N_j.(Names_i{ii});
             else
                 N_j_temp=N_j;
             end
-            if isa(d_grid,'struct')
+            if isstruct(d_grid)
                 d_grid_temp=d_grid.(Names_i{ii});
             else
                 d_grid_temp=d_grid;
             end
-            if isa(a_grid,'struct')
+            if isstruct(a_grid)
                 a_grid_temp=a_grid.(Names_i{ii});
             else
                 a_grid_temp=a_grid;
             end
-            if isa(z_grid,'struct')
+            if isstruct(z_grid)
                 z_grid_temp=z_grid.(Names_i{ii});
             else
                 z_grid_temp=z_grid;
@@ -1236,7 +1246,7 @@
             FullParamNames=fieldnames(Parameters);
             nFields=length(FullParamNames);
             for kField=1:nFields
-                if isa(Parameters.(FullParamNames{kField}), 'struct') % Check for permanent type in structure form
+                if isstruct(Parameters.(FullParamNames{kField})) % Check for permanent type in structure form
                     names=fieldnames(Parameters.(FullParamNames{kField}));
                     for jj=1:length(names)
                         if strcmp(names{jj},Names_i{ii})

EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_AggVars_Case1.m

@@ -83,13 +83,20 @@
 
 %%
 a_gridvals=CreateGridvals(n_a,a_grid,1);
-z_gridvals=CreateGridvals(n_z,z_grid,1);
 
 PolicyPath=reshape(PolicyPath,[size(PolicyPath,1),N_a,N_z,T]);
 % Create PolicyValuesPath from PolicyIndexesPath for use in calculating model stats
 PolicyValuesPath=PolicyInd2Val_InfHorz_TPath(PolicyPath,n_d,n_a,n_z,T,d_grid,a_grid,simoptions,1);
 PolicyValuesPath=permute(reshape(PolicyValuesPath,[size(PolicyValuesPath,1),N_a,N_z,T]),[2,3,1,4]); %[N_a,N_z,l_d+l_a,T-1]
 
+%% Switch to z_gridvals
+l_z=length(n_z);
+if all(size(z_grid)==[sum(n_z),1])
+    z_gridvals=CreateGridvals(n_z,z_grid,1); % The 1 at end indicates want output in form of matrix.
+elseif all(size(z_grid)==[prod(n_z),l_z])
+    z_gridvals=z_grid;
+end
+
 %%
 AgentDistPath=reshape(AgentDistPath,[N_a,N_z,T]);
 

EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_CorrTransProbs_InfHorz.m → EvaluateFnOnAgentDist/TransPath/EvalFnOnTransPath_AutoCorrTransProbs_InfHorz.m

@@ -1,4 +1,4 @@
-function CorrTransProbsPath=EvalFnOnTransPath_CorrTransProbs_InfHorz(FnsToEvaluate,AgentDistPath,PolicyPath,PricePath,ParamPath, Parameters, T, n_d, n_a, n_z, d_grid, a_grid,z_grid, pi_z,simoptions)
+function CorrTransProbsPath=EvalFnOnTransPath_AutoCorrTransProbs_InfHorz(FnsToEvaluate,AgentDistPath,PolicyPath,PricePath,ParamPath, Parameters, T, n_d, n_a, n_z, d_grid, a_grid,z_grid, pi_z,simoptions)
 % Returns stats on (auto) correlation and transition probabilities
 % You must input the names for the FnsToEvaluate that you want the transition probabilities for (by default it won't do any)
 % Done as simoptions.transprobs
@@ -291,4 +291,4 @@
             end
         end
     end
-end
+end