It looks very tricky to make these two play together. Should I try? Or keep a safe distance?
I see that ExpAsset works with GI, and that looks hard enough.
1 Like
I plan on implementing this next month, so if you care to wait that is an option 
It will work, but obvs not a trivial one to get going.
2 Likes
Sounds like a good deal. I suspect refactoring first will ultimately make that easier, but that’s just my opinion.
I’m taking a serious run at ExpAsset in TPath with slowOLG today. I expect to produce a few more cherries to pick along the way.
You’ve probably already noticed, but for value function you largely just copy-paste from the FHorz commands to create the SingleStep commands. It is much the same thing, just the setup has to be treated very slightly differently.
Indeed. Quick (?) question:
In the transition path loop, I see this (updated 3 years ago):
%% Iterate backwards through j.
for reverse_j=1:N_j-1
jj=N_j-reverse_j;
%% [...]
% VKronNext_j=Vtemp_j; % Has been presaved before it was
% VKronNext_j=sum(VKronNext_j.*pi_e_J(1,1,:,jj+1),3); % Take expectations over e
%
% Vtemp_j=V(:,:,jj); % Grab this before it is replaced/updated
VKronNext_j=Vnext(:,:,1,jj+1);
whereas in the stand-alone I see all this extra stuff (mostly from 2 years ago, and other updates more recent):
%% Iterate backwards through j.
for reverse_j=1:N_j-1
jj=N_j-reverse_j;
%% [...]
aprimeFnParamsVec=CreateVectorFromParams(Parameters, aprimeFnParamNames,jj);
[a2primeIndex,a2primeProbs]=CreateExperienceAssetFnMatrix_Case1(aprimeFn, n_d2, n_a2, d2_grid, a2_grid, aprimeFnParamsVec,2); % Note, is actually aprime_grid (but a_grid is anyway same for all ages)
% Note: aprimeIndex is [N_d2,N_a2], whereas aprimeProbs is [N_d2,N_a2]
aprimeIndex=repelem((1:1:N_a1)',N_d2,N_a2)+N_a1*repmat(a2primeIndex-1,N_a1,1,1); % [N_d2*N_a1,N_a2]
aprimeplus1Index=repelem((1:1:N_a1)',N_d2,N_a2)+N_a1*repmat(a2primeIndex,N_a1,1,1); % [N_d2*N_a1,N_a2]
aprimeProbs=repmat(a2primeProbs,N_a1,1,N_z); % [N_d2*N_a1,N_a2,N_z]
EVpre=sum(shiftdim(pi_e_J(:,N_j),-2).*V(:,:,:,jj+1),3); % First, switch V_Jplus1 into Kron form
Vlower=reshape(EVpre(aprimeIndex(:),:),[N_d2*N_a1,N_a2,N_z]);
Vupper=reshape(EVpre(aprimeplus1Index(:),:),[N_d2*N_a1,N_a2,N_z]);
% Skip interpolation when upper and lower are equal (otherwise can cause numerical rounding errors)
skipinterp=(Vlower==Vupper);
aprimeProbs(skipinterp)=0; % effectively skips interpolation
% Switch EV from being in terps of a2prime to being in terms of d2 and a2
EV=aprimeProbs.*Vlower+(1-aprimeProbs).*Vupper; % (d2,a1prime,a2,u,zprime)
% Already applied the probabilities from interpolating onto grid
EV=EV.*shiftdim(pi_z_J(:,:,jj)',-2);
EV(isnan(EV))=0; % remove nan created where value fn is -Inf but probability is zero
EV=squeeze(sum(EV,3));
% EV is over (d2,a1prime,a2,z)
Is this just code drift with the stand-alone being updated and the TransitionPath stuff being less up-to-date? Since my first cut will be in the transition path world, I’d like to start from the correct starting point.
Since the real action is creating the return matrix and calculating V and Policy I think I am safe grafting the loop from the stand-alone version into the TransitionPath code. But if not, sing out…
All the stuff around ‘aprimeFnParamsVec’, ‘aprimeIndex’, ‘Vupper’, is the experience asset being handled.
With experience asset, we have Vnext(a1prime,a2prime,zprime), and a2prime(d2,a2), and that code is switching Vnext to instead by Vnext(d2,a1prime,a2,zprime) and then take expectations to get EV(d2,a1prime,a2,z) which can then be put together with the ReturnFn which is (d1,d2,a1prime,a1,a2,z)
[a2 is the experience asset, so a2prime cannot be chosen directly and hence is not in the ReturnFn]
I figured as much. I’m already using that as my starting point for the update. Thanks for the quick response!
Do I read correctly that there should be a UnKronPolicyIndexes_Case1_FHorz (with and without _e) for the ExpAsset case and there is none such yet?
Answering my own question: yes. There is a UnKronPolicyIndexes_InfHorz_TransPath_ExpAsset but it is not yet complete.
I copied the guts of some semiz code to kickstart experienceasset. I have code connectivity running through transition paths but am sure I don’t quite have the maths right for setting Policy (in particular, the PolicyIndexes for (d2,aprime)).
I’m at a meeting this evening but will post code for review soon. Once I have the slowOLG _e case working, the other few should come quickly.
Alas, my model won’t fit into memory with fastOLG:
The matrix multiply in AgentDist_FHorz_TPath_SingleStep_IterFast_e_raw needs to multiply a 1M row vector by a 3.1M column vector. That’s a large matrix multiplication! My N_j dimension is 17 (5 year periods).
I think I’ve been able to relax the assumption that simoptions.fastOLG=1 in Transpath_Case1_FHorz_PType.m, but I’m unsure about comments like this when setting up jequalOneDist:
jequalOneDist_temp=reshape(jequalOneDist_temp,[N_a*N_z*N_e,1]); % simoptions.fastOLG==1
How would this be different when simoptions.fastOLG=0? After all, the main difference between fastOLG and slowOLG is how N_j contributes to the shape of the array, and since jequalOneDist is itself a by-N_j entity, when would its shape be otherwise (and not automatically reshaped as needed)?
I would double-check that simoptions.fastOLG=1 is being done on cpu not gpu (put a whos just inside the agentdist command, make sure that everything that gets passed into sparse() is on the cpu), before I tried to do simoptions.fastOLG=0.
Thanks! I’ve made it to EvalFnOnAgentDist_AggVars_FHorz_fastOLG which mentions in a comment the optionality of fastOLG, but does not present a fastOLGpolicy parameter. Thoughts on that?
And I need to add l_aprime==4 logical to my l_d==2 case, which is easy enough, though I need to specialize it for ExpAsset.
If you look at a fastOLG shooting command and a slowOLG shooting command, from memory, you can see both use EvalFnOnAgentDist_AggVars_FHorz_fastOLG but one has to do a ‘permute’ of Policy before passing it as an input and the other doesn’t.
[Note, from memory, so I could be mistaken]
Looking at the code, what I see is that the permute step is done for the respective cases before calling the Eval on the AggVars, so I think the comment in the code needs to reflect that there’s just one way to do it. Here’s the snippet that preps the call from the slowOLG shooting function:
PolicyValuesPath=permute(PolicyValuesPath,[2,3,1,4]); %[N_a,N_j,l_d+l_aprime,T-1] % fastOLG ordering is needed for AggVars
So, I’m creating a hand-crafted version of the Eval on AggVars for my experienceasset case, but we’ll need your code generator to flesh it out properly.
This parameter
z_gridvals_J_fastOLG(1,:,:,1)
is reading as a multidimensional vector when passed to arrayfun, and my AggVar function is seeing that as feeding z only and not also e.
I don’t know the trick for splitting it up so that the z part of the joint grid feeds z and the e part feeds e.
This comment in EvalFnOnAgentDist_AggVars_FHorz_fastOLG tells me to send it a joint grid, but not how to unpack it:
% If you have e or semiz, just disguise them as z for this command
You should not be unpacking it. When it comes to FnsToEvaluate there is no distinction between z and e. Hence why you just disguies the e as z, and then you can use the same code that handles just z.
[Obviously the z vs e distinction is important for value fns, agent dist, etc., but it is irrelevant for FnsToEvaluate, so for AggVars, AllStats, etc.]
Here’s the error message:
Caused by:
@(labor,buyhouse,sprime,aprime,cprime,hprime,s,a,car,h,solarpv,z,e,kappa_j,Lhscale)labor*kappa_j*exp(z+e)*Lhscale takes at least 15 inputs. 14 were supplied.
From this call:
Values=arrayfun(FnsToEvaluate{ff}, PolicyValues_d(:,:,:,1), PolicyValues_d(:,:,:,2), PolicyValues_aprime(:,:,:,1), PolicyValues_aprime(:,:,:,2), PolicyValues_aprime(:,:,:,3), PolicyValues_aprime(:,:,:,4), a_gridvals(:,1), a_gridvals(:,2), a_gridvals(:,3), a_gridvals(:,4), a_gridvals(:,5), z_gridvals_J_fastOLG(1,:,:,1), ParamCell{:});
All the data is there, it’s just formatted in a way that’s not making it through. Here’s an example of a call that doesn’t error out, because the zs and es are filling the parameter buckets in the way the GPU likes:
Fmatrix=arrayfun(ReturnFn, d_gridvals(:,1),d_gridvals(:,2), a1prime1vals,a1prime2vals,a1prime3vals,a1prime4vals, a1vals,a2vals,a3vals,a4vals, expassetvals, z1vals, e1vals, ReturnFnParamsCell{:});
I am this >< close!
Aha! I went back and checked how it was possible that I was getting this error with TPath and not without. Turns out that my FnsToEvaluate and FnsToEvaluate2 functions did not replicate the household labor function signature exactly. I see that I got the signature right in one case and wrong with another. I will give it a go with the correct function signature.
I hit a new snag. In TransitionPath_Case1_Fhorz_PType.m, there’s this code at line 426:
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});
% Figure out FnsToEvaluateParamNames
temp=getAnonymousFnInputNames(FnsToEvaluate.(FnNames{kk}).(Names_i{ii}));
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)
On the last line you can see that if there is an ‘e’ in the grid (l_e>0), then the position normally occupied by the e parameter will be skipped in the list of parameters named by temp. In my case:
temp =
15×1 cell array
{'labor' }
{'buyhouse'}
{'sprime' }
{'aprime' }
{'cprime' }
{'hprime' }
{'s' }
{'a' }
{'car' }
{'h' }
{'solarpv' }
{'z' }
{'w' }
{'kappa_j' }
{'Lhscale' }
becomes (with special experience asset code temp{l_d+l_aprime+l_a+l_z+l_e+1:end}, meaning that the l_aprime is not consuming more than it should):
ans =
1×2 cell array
{'kappa_j'} {'Lhscale'}
the w is being absorbed as if it were e. So if I construct a parameter list that’s correct for later in the Transition Paths (no e), it gets botched early on. And if I construct a parameter list correct for this initial setup, it gets botched when the Transition Path function doesn’t expect to see an e.
What am I missing?