This morning I discovered a divergence between TransitionPath_Case1_FHorz and TransitionPath_Case1_FHorz_PType introduced by “major overhaul of FHorz TPath” (commit 79fa959). As I contemplated what it would take to port that across to the PType version, I contemplated how simple it might be to refactor the PType code so that it would operate like this (from a top-down perspective):
Step 1: for each PType in the incoming parameter list, collect the PType-specific values into a fresh list of parameters. This is the bulk of what is already done at the start of TransitionPath_Case1_FHorzPType.
Step 2: Call the relevant FHorz/InfHorz transition path functions that have no knowledge of PTypes whatsover, but are gracious enough to return all values that might later be interesting from a PType-specific perspective. This already happens.
Step 3: Assemble PType-aware results from Step 2. Which we have to do anyway, the only question is whether to do all of the heavy lifting by duplicating the functionality of the non-PType-aware functions, or simply collecting their results and organizing it into PType structures.
From a bottom-up perspective, we see that the PType-specific code calls things like TransitionPath_FHorz_PType_singlepath_fastOLG_e_raw whereas the non-PType code calls things like TransitionPath_FHorz_shooting_fastOLG_e. The differences between the two are:
The shooting code runs a loop with these major steps
- First, go from T-1 to 1 calculating the Value function and Optimal policy function at each steps
- Modify PolicyIndexesPath into forms needed for forward iteration
- Iterate forward over t: iterate agent dist, calculate aggvars, evaluate general eqm, which loops over:
- Get the current optimal policy, and iterate the agent dist (conditional logic!)
- AggVars
- General Eqm Eqns (missing from singlepath!)
- Update prices, give verbose feedback, and check for convergence
There are 16 versions of this code owing to the combinatorial explosion of d/nod, z/noz, e/noe, fastOLG/no fastOLG, each of which works on different Policy shapes for each permutation. Which is very understandable for steps 1 and 2, which are thick the the shaping specializations.
Step 3.1 has a small bit of shaping conditions, but in this case (where we have d, z, and e), there’s yet conditional logic:
%% Get the current optimal policy, and iterate the agent dist
if transpathoptions.zepathtrivial==0
ze_gridvals_J_fastOLG=transpathoptions.ze_gridvals_J_T_fastOLG(:,:,:,:,tt);
end
if transpathoptions.zpathtrivial==0
pi_z_J_sim=transpathoptions.pi_z_J_sim_T(:,:,:,tt);
end
if transpathoptions.epathtrivial==0
pi_e_J_sim=transpathoptions.pi_e_J_sim_T(:,:,tt); % (a,j,z)-by-e
end
if transpathoptions.trivialjequalonedist==0
jequalOneDist=jequalOneDist_T(:,tt+1); % Note: t+1 as we are about to create the next period AgentDist
end
It’s not nearly as complicated as seen in 1 and 2, and many of the conditions tested (such as transpathoptions.trivialjequalonedist) are tested across the many cases.
The singlepath code is not a loop, but executes steps 1, 2, and 3 (without using the General Eqm Eqns to generate a new Price Path that the shooting code does in 3.3). And it does not do Step 4, which is a step that appears perfectly replicated across all shooting versions with no differences between the 16 variations. That step, which must be PType-aware, is done by the enclosing logic of TransitionPath_Case1_FHorz_PType. That function does call a PType shooting dispatcher that does call down to specialized singlepath versions, but the shooting dispatcher does not itself run the loop that the shooting functions implement in the non-PType world.
I would argue that there would be no loss in performance, and a great improvement in maintainability, if the shooting code for the non-PType case were refactored to match the factorization of the PType case. This would allow both cases could use the identical singlepath logic for steps 1-3. By removing the complex shape-aware logic from the shooting functions to the singlepath functions, all the non-PType shooting code can be simplified to calling a singlepath dispatcher (as the PType code already does). This would collapse the shooting codes to the same three versions that PType has: FHorz, FHorz_ExpAsset, and InfHorz.
Then the PType and non-PType shooting codes can be further harmonized into a single case because the outer path loop does not depend on PType, either:
%% Now we just check for convergence, update prices, and give some feedback on progress
% See how far apart the price paths are
PricePathDist=max(abs(reshape(PricePathNew(1:T-1,:)-PricePathOld(1:T-1,:),[numel(PricePathOld(1:T-1,:)),1])));
% Notice that the distance is always calculated ignoring the time t=T periods, as these needn't ever converges
if transpathoptions.verbose==1
% Would be nice to have a way to get the iteration count without having the whole printout of path values (I think that would be useful?)
pathnametitles
[PricePathOld,PricePathNew]
end
if transpathoptions.graphpricepath==1
% Do a graph of the GE prices
figure(1);
for pp=1:length(PricePathNames)
if PTypeStructure.PricePath_Idependsonptype(pp)==0
subplot(nrows_pricepath,ncolumns_pricepath,pp); plot(1:1:T,PricePathOld(:,pp_indexinpricepath(pp)))
else
subplot(nrows_pricepath,ncolumns_pricepath,pp); plot(1:1:T,PricePathOld(:,pp_indexinpricepath(pp)))
hold on
for ii=2:N_i
subplot(nrows_pricepath,ncolumns_pricepath,pp); plot(1:1:T,PricePathOld(:,pp_indexinpricepath(pp)+ii-1))
end
hold off
end
title(PricePathNames{pp})
end
end
if transpathoptions.graphaggvarspath==1
% Do a graph of the AggVar paths
figure(2);
for aa=1:length(FullAggVarNames)
subplot(nrows_aggvars,ncolumns_aggvars,aa); plot(1:1:T-1,AggVarsPooledPath(aa,:))
title(FullAggVarNames{aa})
if aa_aggvarbyptype(aa)==1 % We care about the value conditional on ptype
hold on
for ii=1:N_i
subplot(nrows_aggvars,ncolumns_aggvars,aa); plot(1:1:T-1,AggVarsFullPath(aa,:,ii),'-')
end
hold off
end
end
end
if transpathoptions.graphGEcondns==1
% Do a graph of the General eqm conditions
figure(3);
for gg=1:length(GEeqnNames)
if transpathoptions.GEptype(gg)==0
subplot(nrows_GEcondns,ncolumns_GEcondns,gg); plot(1:1:T,GECondnPath(:,gg_indexinGEcondns(gg)))
else
subplot(nrows_GEcondns,ncolumns_GEcondns,gg); plot(1:1:T,GECondnPath(:,gg_indexinGEcondns(gg)))
hold on
for ii=2:N_i
subplot(nrows_GEcondns,ncolumns_GEcondns,gg); plot(1:1:T,GECondnPath(:,gg_indexinGEcondns(gg)+ii-1))
end
hold off
end
title(GEeqnNames{gg})
end
end
% Set price path to be 9/10ths the old path and 1/10th the new path (but making sure to leave prices in periods 1 & T unchanged).
if transpathoptions.weightscheme==0
PricePathOld=PricePathNew; % The update weights are already in GEnewprice setup
elseif transpathoptions.weightscheme==1 % Just a constant weighting
PricePathOld(1:T-1,:)=transpathoptions.oldpathweight.*PricePathOld(1:T-1,:)+(1-transpathoptions.oldpathweight).*PricePathNew(1:T-1,:);
elseif transpathoptions.weightscheme==2 % A exponentially decreasing weighting on new path from (1-oldpathweight) in first period, down to 0.1*(1-oldpathweight) in T-1 period.
% I should precalculate these weighting vectors
Ttheta=transpathoptions.Ttheta;
PricePathOld(1:Ttheta,:)=transpathoptions.oldpathweight*PricePathOld(1:Ttheta,:)+(1-transpathoptions.oldpathweight)*PricePathNew(1:Ttheta,:);
PricePathOld(Ttheta:T-1,:)=((transpathoptions.oldpathweight+(1-exp(linspace(0,log(0.2),T-Ttheta)))*(1-transpathoptions.oldpathweight))'*ones(1,l_p)).*PricePathOld(Ttheta:T-1,:)+((exp(linspace(0,log(0.2),T-Ttheta)).*(1-transpathoptions.oldpathweight))'*ones(1,l_p)).*PricePathNew(Ttheta:T-1,:);
elseif transpathoptions.weightscheme==3 % A gradually opening window.
if (pathcounter*3)<T-1
PricePathOld(1:(pathcounter*3),:)=transpathoptions.oldpathweight*PricePathOld(1:(pathcounter*3),:)+(1-transpathoptions.oldpathweight)*PricePathNew(1:(pathcounter*3),:);
else
PricePathOld(1:T-1,:)=transpathoptions.oldpathweight*PricePathOld(1:T-1,:)+(1-transpathoptions.oldpathweight)*PricePathNew(1:T-1,:);
end
elseif transpathoptions.weightscheme==4 % Combines weightscheme 2 & 3
if (pathcounter*3)<T-1
PricePathOld(1:(pathcounter*3),:)=((transpathoptions.oldpathweight+(1-exp(linspace(0,log(0.2),pathcounter*3)))*(1-transpathoptions.oldpathweight))'*ones(1,l_p)).*PricePathOld(1:(pathcounter*3),:)+((exp(linspace(0,log(0.2),pathcounter*3)).*(1-transpathoptions.oldpathweight))'*ones(1,l_p)).*PricePathNew(1:(pathcounter*3),:);
else
PricePathOld(1:T-1,:)=((transpathoptions.oldpathweight+(1-exp(linspace(0,log(0.2),T-1)))*(1-transpathoptions.oldpathweight))'*ones(1,l_p)).*PricePathOld(1:T-1,:)+((exp(linspace(0,log(0.2),T-1)).*(1-transpathoptions.oldpathweight))'*ones(1,l_p)).*PricePathNew(1:T-1,:);
end
end
TransPathConvergence=PricePathDist/transpathoptions.tolerance; %So when this gets to 1 we have convergence
if transpathoptions.verbose==1
fprintf('Number of iterations on transition path: %i \n',pathcounter)
fprintf('Current distance between old and new price path (in L-Infinity norm): %8.6f \n', PricePathDist)
fprintf('Ratio of current distance to the convergence tolerance: %.2f (convergence when reaches 1) \n',TransPathConvergence)
end
if transpathoptions.historyofpricepath==1
% Store the whole history of the price path and save it every ten iterations
PricePathHistory{pathcounter,1}=PricePathDist;
PricePathHistory{pathcounter,2}=PricePathOld;
if rem(pathcounter,10)==1
save ./SavedOutput/TransPath_Internal.mat PricePathHistory
end
end
pathcounter=pathcounter+1;
A pricepath function implementing the above logic would save 16 replications in the non-PType case, and could re-use the code in the PType case.
Then in all cases, the principal varieties (FHorz, FHorzExpAsset, and InfHorz) can simply return the PricePathOld and GECondnPath as they should, to be used directly in the non-PType case, or sorted into PType structures in the PType case.
Otherwise we create a large amount of work trying to keep the PType code up to date with the non-PType code (and vice-versa).