From dae1b434edec76c9802a4e1171f55f43d9568dbd Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Fri, 24 Jul 2026 05:58:05 +1200
Subject: [PATCH 1/5] Demonstrate optimize_nprobs

LifeCycleModel21 demonstrates grid interpolation
LifeCycleModel41 demonstrates ExpAsset
LifeCycleModel42 demonstrates ExpAssetu

In all cases I adjusted some grid parameters based on feedback about how well (or how poorly) the model fit the initial grid conditions.  LifeCycleModel42 so far exceeded the confines of the grid that I changed a bit more of the model to make it fit.
---
 Models20to30/LifeCycleModel21.m          | 12 ++++++++--
 Models40to44/LifeCycleModel41.m          |  9 ++++----
 Models40to44/LifeCycleModel42.m          | 28 ++++++++++++++----------
 Models40to44/LifeCycleModel42_ReturnFn.m |  6 ++---
 4 files changed, 35 insertions(+), 20 deletions(-)

diff --git a/Models20to30/LifeCycleModel21.m b/Models20to30/LifeCycleModel21.m
index f22c798e..c9557a05 100644
--- a/Models20to30/LifeCycleModel21.m
+++ b/Models20to30/LifeCycleModel21.m
@@ -77,7 +77,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=15*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % Grid for labour choice
 h_grid=linspace(0,1,n_d)'; % Notice that it is imposing the 0<=h<=1 condition implicitly
@@ -137,6 +137,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=1;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % Note: Because we have simoptions.ExogShockFn, what we input for z_grid and  pi_z will just be ignored.
 
@@ -165,7 +168,7 @@
 title('Life Cycle Profile: Fraction Time Worked (h)')
 subplot(5,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z h)')
-subplot(5,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(5,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(5,1,4); plot(1:1:Params.J,[AgeConditionalStats.fractionunemployed.Mean(1:Params.Jr-1),nan(1,Params.J-Params.Jr+1)])
 title('Life Cycle Profile: Fraction Unemployment (z==0)')
@@ -194,11 +197,16 @@
 ReturnFn=@(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,wg1,wg2,wg3,beta,sj) ...
     LifeCycleModel8_ReturnFn(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,wg1,wg2,wg3,beta,sj);
 [V_nomedical, Policy_nomedical]=ValueFnIter_Case1_FHorz(n_d,n_a,n_z,N_j, d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], vfoptions);
+simoptions.optimize_nProbs=1;
+simoptions.verbose=1;
 StationaryDist_nomedical=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy_nomedical,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 AgeConditionalStats_nomedical=LifeCycleProfiles_FHorz_Case1(StationaryDist_nomedical,Policy_nomedical,FnsToEvaluate,Params,[],n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,simoptions);
 
 figure(2)
+hold on
 plot(1:1:Params.J,AgeConditionalStats.assets.Mean,1:1:Params.J,AgeConditionalStats_nomedical.assets.Mean)
+plot(1:1:Params.J,AgeConditionalStats.assets.Mean,1:1:Params.J,AgeConditionalStats_nomedical.assets.Maximum)
+hold off
 title('Life Cycle Profile: Assets (a)')
 legend('Medical Expense Shocks','No Medical Shocks')
 % Notice that medical expense shocks late in life cause elderly households
diff --git a/Models40to44/LifeCycleModel41.m b/Models40to44/LifeCycleModel41.m
index 7c730e74..ac427ae6 100644
--- a/Models40to44/LifeCycleModel41.m
+++ b/Models40to44/LifeCycleModel41.m
@@ -30,7 +30,7 @@
 
 % Grid sizes to use
 n_d=2; % Endogenous labour choice (working or not-working)
-n_a=[201,21]; % Endogenous asset holdings, female-labor-force-participation
+n_a=[201,31]; % Endogenous asset holdings, female-labor-force-participation
 n_z=11; % Exogenous labor productivity units shock
 N_j=Params.J; % Number of periods in finite horizon
 
@@ -96,7 +96,7 @@
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
 asset_grid=10*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
-h_grid=linspace(0.1,3,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
+h_grid=linspace(0.1,4,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
 % Note: deliberately omit 0 from h_grid
 
 % First, the AR(1) process z
@@ -140,7 +140,8 @@
 % vfoptions.lowmemory=1; % default=0, set =1 to use loops (over e, or z if no e) if you get a gpu out-of-memory error, the loops reduce memory use but slow the runtimes [models with e & z can set =2]
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
-
+simoptions.optimize_nProbs=1;
+simoptions.verbose=1;
 
 %% Now, create the return function
 DiscountFactorParamNames={'beta','sj'};
@@ -252,7 +253,7 @@
 title('Life Cycle Profile: Male Earnings (y_m)')
 subplot(5,1,4); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean)
 title('Life Cycle Profile: Assets (a)')
-subplot(5,1,5); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Mean)
+subplot(5,1,5); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Maximum); hold off
 title('Life Cycle Profile: female labor force participation history (h)')
 % Remember, we added a huge childcare cost, so women choose not to work
 % during ages 27-31. You can see that their human capital (femaleLFPH) has
diff --git a/Models40to44/LifeCycleModel42.m b/Models40to44/LifeCycleModel42.m
index 76d65185..bc76e84f 100644
--- a/Models40to44/LifeCycleModel42.m
+++ b/Models40to44/LifeCycleModel42.m
@@ -64,9 +64,9 @@
 
 % Uncertain Human Capital
 Params.alpha_h=0.6; % returns to scale of human capital production fn
-Params.delta_h=0.02; % depreciation rate of human capital
-Params.ability=3; % ability to produce human capital
-Params.sscaler=5; % scale from units of study-time into something more appropriate for the human capital production function
+Params.delta_h_j=0.02*[ones(1,31), linspace(1,2,46-31), linspace(2,4,61-46), 4*ones(1,81-61)]; % depreciation rate of human capital
+Params.ability_j=[linspace(3,1,21), ones(1,81-21)]; % ability to produce human capital
+Params.sscaler=0.5; % scale from units of study-time into something more appropriate for the human capital production function
 % i.i.d. shocks that make human capital uncertain
 Params.sigma_epsilon_u=0.01;
 
@@ -94,13 +94,17 @@
 Params.wg2=3; % degree to which bequests are a luxury good (>=1; =1 would be a normal good)
 Params.wg3=Params.sigma; % By using the same curvature as the utility of consumption it makes it much easier to guess appropriate parameter values for the warm glow
 
+% Tuition
+% When tuition is more expensive (8), agents get advanced degree later, when they can afford it
+% When tutition is less expensive (6), they hurry up and study while kappa_j is low
+Params.tuition=7;
 
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
 asset_grid=10*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
-h_grid=linspace(1,4,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
+h_grid=linspace(1,20,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
 % Note: deliberately omit 0 from h_grid
 
 % First, the AR(1) process z
@@ -130,7 +134,7 @@
 % the experience asset).
 
 % aprimeFn gives the value of hprime
-vfoptions.aprimeFn=@(s,h,u,alpha_h, delta_h, ability,sscaler) u*(ability*(h*(sscaler*s))^alpha_h+h*(1-delta_h));
+vfoptions.aprimeFn=@(s,h,u,alpha_h, delta_h_j, ability_j,sscaler) u*(ability_j*(h*(sscaler*s))^alpha_h+h*(1-delta_h_j));
 % The first three inputs must be (d,a,u) [in the sense of aprime(d,a,u)], then any parameters
 
 % We also need to tell simoptions about the experienceassetu
@@ -138,6 +142,8 @@
 simoptions.aprimeFn=vfoptions.aprimeFn;
 simoptions.d_grid=d_grid; % Needed to handle aprimeFn 
 simoptions.a_grid=a_grid; % Needed to handle aprimeFn
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 
 % And we need to define the i.i.d. shocks u
 [u_grid,pi_u]=discretizeAR1_FarmerToda(0,0,Params.sigma_epsilon_u,n_u);
@@ -162,12 +168,12 @@
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
 
-% To better understand the human capital production function, here is a graph of it
+% To better understand the human capital production function, here is a graph of it (at age 1)
 hprime_shu=zeros([n_a,n_u]);
 for hh=1:n_a(2)
     for ss=1:n_d(2)
         for uu=1:n_u
-            hprime_shu(ss,hh,uu)=vfoptions.aprimeFn(s_grid(ss),h_grid(hh),u_grid(uu),Params.alpha_h,Params.delta_h,Params.ability,Params.sscaler);
+            hprime_shu(ss,hh,uu)=vfoptions.aprimeFn(s_grid(ss),h_grid(hh),u_grid(uu),Params.alpha_h,Params.delta_h_j(1),Params.ability_j(1),Params.sscaler);
         end
     end
 end
@@ -187,8 +193,8 @@
 DiscountFactorParamNames={'beta','sj'};
 
 % Use 'LifeCycleModel42_ReturnFn'
-ReturnFn=@(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj)...
-    LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj);
+ReturnFn=@(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj,tuition)...
+    LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj,tuition);
 % Notice how we have (l,s,aprime,a,h,z,...)
 % Follow same decision-next endo-endo-exo ordering as usual, but because h
 % is an experienceassetu, we do not include hprime as it is not chosen
@@ -257,7 +263,7 @@
 %% Initial distribution of agents at birth (j=1)
 % Before we plot the life-cycle profiles we have to define how agents are at age j=1. We will give them all zero assets.
 jequaloneDist=zeros([n_a,n_z],'gpuArray'); % Put no households anywhere on grid
-jequaloneDist(1,4,floor((n_z+1)/2))=1; % All agents start with zero assets, h_grid(4) of human capital, and the median shock [h_grid(4) is roughly same as y_m(1), just my arbitrary decision]
+jequaloneDist(1,2,floor((n_z+1)/2))=1; % All agents start with zero assets, h_grid(2) of human capital, and the median shock [h_grid(2) is roughly same as y_m(1), just my arbitrary decision]
 
 %% We now compute the 'stationary distribution' of households
 % Start with a mass of one at initial age, use the conditional survival
@@ -295,7 +301,7 @@
 title('Life Cycle Profile: Study time (s)')
 subplot(5,1,3); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Earnings (w h z l)')
-subplot(5,1,4); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Mean)
+subplot(5,1,4); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Maximum); hold off
 title('Life Cycle Profile: Human Capital (h)')
 subplot(5,1,5); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean)
 title('Life Cycle Profile: Assets (a)')
diff --git a/Models40to44/LifeCycleModel42_ReturnFn.m b/Models40to44/LifeCycleModel42_ReturnFn.m
index 61fcea76..13cfbb2e 100644
--- a/Models40to44/LifeCycleModel42_ReturnFn.m
+++ b/Models40to44/LifeCycleModel42_ReturnFn.m
@@ -1,12 +1,12 @@
-function F=LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj)
+function F=LifeCycleModel42_ReturnFn(l,s,aprime,a,h,z,w,sigma,eta,psi,agej,Jr,pension,r,wg1,wg2,wg3,beta,sj,tuition)
 
 leisure=1-s-l; % one unit of time, minus time spent studying and working
 
 F=-Inf;
 if agej<Jr % If working age
-    c=w*h*z*l+(1+r)*a-aprime; % Add z here
+    c=w*h*z*l+(1+r)*a-aprime-s*tuition; % Add z here
 else % Retirement
-    c=pension+(1+r)*a-aprime;
+    c=pension+(1+r)*a-aprime-s*tuition/2;
 end
 
 if c>0 && leisure>0

From 417eea03613dc269f5a9ae49579785e76fafdcd9 Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Fri, 24 Jul 2026 07:40:32 +1200
Subject: [PATCH 2/5] Additional demonstrations using `e`

Demonstrates noz_e and e cases.
---
 Models20to30/LifeCycleModel23.m    | 3 +++
 ModelsAppendixA/LifeCycleModelA6.m | 6 ++++--
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/Models20to30/LifeCycleModel23.m b/Models20to30/LifeCycleModel23.m
index d537505b..3f99841d 100644
--- a/Models20to30/LifeCycleModel23.m
+++ b/Models20to30/LifeCycleModel23.m
@@ -145,6 +145,9 @@
 vfoptions.ngridinterp=20; % 20 evenly-spaced points between each pair of consecutive a_grid points
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
+simoptions.optimize_nProbs=1;
+simoptions.verbose=1;
+simoptions.tolerance=1e-7;
 
 
 % Grid for labour choice
diff --git a/ModelsAppendixA/LifeCycleModelA6.m b/ModelsAppendixA/LifeCycleModelA6.m
index ac3cbb71..230c810a 100644
--- a/ModelsAppendixA/LifeCycleModelA6.m
+++ b/ModelsAppendixA/LifeCycleModelA6.m
@@ -139,7 +139,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=16*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % Grid for labour choice
 h_grid=linspace(0,1,n_d)'; % Notice that it is imposing the 0<=h<=1 condition implicitly
@@ -183,6 +183,8 @@
 vfoptions.ngridinterp=20; % 20 evenly-spaced points between each pair of consecutive a_grid points
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
+simoptions.optimize_nProbs=1;
+simoptions.verbose=1;
 
 %% Now, create the return function
 DiscountFactorParamNames={'beta','sj'};
@@ -260,7 +262,7 @@
 title('Life Cycle Profile: Fraction Time Worked (h)')
 subplot(3,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j h e1 e2)')
-subplot(3,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(3,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); subplot(3,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); plot(1:1:Params.J,AgeConditionalStats.assets.Minimum); hold off
 title('Life Cycle Profile: Assets (a)')
 
 

From 45d891ef3583bf108a46cc293c728bb59309c33a Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Mon, 27 Jul 2026 06:00:40 +1200
Subject: [PATCH 3/5] Fix Maximum plots

---
 Models20to30/LifeCycleModel21.m | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Models20to30/LifeCycleModel21.m b/Models20to30/LifeCycleModel21.m
index c9557a05..1b6c01dc 100644
--- a/Models20to30/LifeCycleModel21.m
+++ b/Models20to30/LifeCycleModel21.m
@@ -205,10 +205,10 @@
 figure(2)
 hold on
 plot(1:1:Params.J,AgeConditionalStats.assets.Mean,1:1:Params.J,AgeConditionalStats_nomedical.assets.Mean)
-plot(1:1:Params.J,AgeConditionalStats.assets.Mean,1:1:Params.J,AgeConditionalStats_nomedical.assets.Maximum)
+plot(1:1:Params.J,AgeConditionalStats.assets.Maximum,1:1:Params.J,AgeConditionalStats_nomedical.assets.Maximum)
 hold off
 title('Life Cycle Profile: Assets (a)')
-legend('Medical Expense Shocks','No Medical Shocks')
+legend('Medical Expense Shocks','No Medical Shocks','Medical Expense Shocks (Maximum)','No Medical (Maximum)')
 % Notice that medical expense shocks late in life cause elderly households
 % to hold more assets (as self-insurance against medical expense shocks)
 

From e65ea7fe4a82a9986f9e3895350a71b30de424ff Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Mon, 27 Jul 2026 09:48:53 +1200
Subject: [PATCH 4/5] Add semiExo cases; cleaning

---
 Models20to30/LifeCycleModel21.m    |  4 ++--
 Models20to30/LifeCycleModel23.m    |  7 ++++---
 Models20to30/LifeCycleModel28.m    |  3 +++
 Models31to35/LifeCycleModel31.m    |  5 ++++-
 Models31to35/LifeCycleModel35.m    | 12 +++++++-----
 Models40to44/LifeCycleModel41.m    |  6 ++++--
 Models40to44/LifeCycleModel42.m    |  5 +++--
 ModelsAppendixA/LifeCycleModelA6.m |  6 ++++--
 8 files changed, 31 insertions(+), 17 deletions(-)

diff --git a/Models20to30/LifeCycleModel21.m b/Models20to30/LifeCycleModel21.m
index 1b6c01dc..ab4b80fd 100644
--- a/Models20to30/LifeCycleModel21.m
+++ b/Models20to30/LifeCycleModel21.m
@@ -139,7 +139,7 @@
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
 
 simoptions.optimize_nProbs=1;
-simoptions.verbose=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % Note: Because we have simoptions.ExogShockFn, what we input for z_grid and  pi_z will just be ignored.
 
@@ -198,7 +198,7 @@
     LifeCycleModel8_ReturnFn(h,aprime,a,z,w,sigma,psi,eta,agej,Jr,pension,r,kappa_j,wg1,wg2,wg3,beta,sj);
 [V_nomedical, Policy_nomedical]=ValueFnIter_Case1_FHorz(n_d,n_a,n_z,N_j, d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], vfoptions);
 simoptions.optimize_nProbs=1;
-simoptions.verbose=1;
+simoptions.verbose=2;
 StationaryDist_nomedical=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy_nomedical,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 AgeConditionalStats_nomedical=LifeCycleProfiles_FHorz_Case1(StationaryDist_nomedical,Policy_nomedical,FnsToEvaluate,Params,[],n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,simoptions);
 
diff --git a/Models20to30/LifeCycleModel23.m b/Models20to30/LifeCycleModel23.m
index 3f99841d..af5ba0bd 100644
--- a/Models20to30/LifeCycleModel23.m
+++ b/Models20to30/LifeCycleModel23.m
@@ -145,9 +145,6 @@
 vfoptions.ngridinterp=20; % 20 evenly-spaced points between each pair of consecutive a_grid points
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
-simoptions.optimize_nProbs=1;
-simoptions.verbose=1;
-simoptions.tolerance=1e-7;
 
 
 % Grid for labour choice
@@ -209,6 +206,10 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
+simoptions.tolerance=1e-7;
 StationaryDist=StationaryDist_Case1_FHorz_PType(jequaloneDist,AgeWeightsParamNames,PTypeDistParamNames,Policy,n_d,n_a,n_z,N_j,N_i,pi_z,Params,simoptions);
 
 %% FnsToEvaluate are how we say what we want to graph the life-cycles of
diff --git a/Models20to30/LifeCycleModel28.m b/Models20to30/LifeCycleModel28.m
index 7cb2dcaf..89c4b4ae 100644
--- a/Models20to30/LifeCycleModel28.m
+++ b/Models20to30/LifeCycleModel28.m
@@ -79,6 +79,7 @@
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
 
+
 % At the bottom of this code/script there are some lines showing you what
 % pi_semiz_J which is created internally looks like.
 
@@ -223,6 +224,8 @@
 % Because evaluating pi_semiz_J requires the d_grid we also have to provide
 simoptions.d_grid=d_grid;
 
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 %% FnsToEvaluate are how we say what we want to graph the life-cycles of
diff --git a/Models31to35/LifeCycleModel31.m b/Models31to35/LifeCycleModel31.m
index cccbce32..6933a9aa 100644
--- a/Models31to35/LifeCycleModel31.m
+++ b/Models31to35/LifeCycleModel31.m
@@ -111,7 +111,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=13*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=16*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % First, the AR(1) process z
 [z_grid,pi_z]=discretizeAR1_FarmerToda(0,Params.rho_z,Params.sigma_epsilon_z,n_z);
@@ -266,6 +266,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % riskyasset requires the grids when simulating the agent distribution to be able to handle aprime(d,u). The grids are passed in simoptions.
 
diff --git a/Models31to35/LifeCycleModel35.m b/Models31to35/LifeCycleModel35.m
index 12903fe1..076fc593 100644
--- a/Models31to35/LifeCycleModel35.m
+++ b/Models31to35/LifeCycleModel35.m
@@ -34,8 +34,8 @@
 Params.J=100-Params.agejshifter; % =81, Number of period in life-cycle
 
 % Grid sizes to use
-n_d=[51,201]; % Decisions: riskyshare, savings
-n_a=[5,201]; % Endogenous housing and asset holdings
+n_d=[51,301]; % Decisions: riskyshare, savings
+n_a=[5,301]; % Endogenous housing and asset holdings
 n_z=7; % Exogenous labor productivity units shock
 n_u=5; % Between period i.i.d. shock
 N_j=Params.J; % Number of periods in finite horizon
@@ -44,7 +44,6 @@
 simoptions.riskyasset=1;
 % When there is more than one endogenous state, the riskyasset is the last one
 
-
 % Specify Epstein-Zin preferences
 vfoptions.exoticpreferences='EpsteinZin';
 vfoptions.EZpositiveutility=0; % Epstein-Zin preferences in utility-units have to be handled differently depending on whether the utility function is positive or negative valued (this is all done internally, you just need to use vfoptions to specify which)
@@ -122,8 +121,8 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=-3+13*(linspace(0,1,n_a(2)))'; % Note, I use equal spacing (normally would put most points near zero)
-% note: will go from -3 to 13-3
+asset_grid=-3+30*(linspace(0,1,n_a(2)))'; % Note, I use equal spacing (normally would put most points near zero)
+% note: will go from -3 to 30-3
 % Make it so that there is a zero assets
 % Find closest to zero assets
 [~,zeroassetindex]=min(abs(asset_grid));
@@ -294,6 +293,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 % riskyasset requires the grids when simulating the agent distribution to be able to handle aprime(d,u). The grids are passed in simoptions.
 
diff --git a/Models40to44/LifeCycleModel41.m b/Models40to44/LifeCycleModel41.m
index ac427ae6..60b996ea 100644
--- a/Models40to44/LifeCycleModel41.m
+++ b/Models40to44/LifeCycleModel41.m
@@ -140,8 +140,7 @@
 % vfoptions.lowmemory=1; % default=0, set =1 to use loops (over e, or z if no e) if you get a gpu out-of-memory error, the loops reduce memory use but slow the runtimes [models with e & z can set =2]
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
-simoptions.optimize_nProbs=1;
-simoptions.verbose=1;
+
 
 %% Now, create the return function
 DiscountFactorParamNames={'beta','sj'};
@@ -226,6 +225,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 %% FnsToEvaluate are how we say what we want to graph the life-cycles of
diff --git a/Models40to44/LifeCycleModel42.m b/Models40to44/LifeCycleModel42.m
index bc76e84f..1bcca56d 100644
--- a/Models40to44/LifeCycleModel42.m
+++ b/Models40to44/LifeCycleModel42.m
@@ -142,8 +142,6 @@
 simoptions.aprimeFn=vfoptions.aprimeFn;
 simoptions.d_grid=d_grid; % Needed to handle aprimeFn 
 simoptions.a_grid=a_grid; % Needed to handle aprimeFn
-simoptions.optimize_nProbs=1;
-simoptions.verbose=2;
 
 % And we need to define the i.i.d. shocks u
 [u_grid,pi_u]=discretizeAR1_FarmerToda(0,0,Params.sigma_epsilon_u,n_u);
@@ -275,6 +273,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 
diff --git a/ModelsAppendixA/LifeCycleModelA6.m b/ModelsAppendixA/LifeCycleModelA6.m
index 230c810a..ece4142d 100644
--- a/ModelsAppendixA/LifeCycleModelA6.m
+++ b/ModelsAppendixA/LifeCycleModelA6.m
@@ -183,8 +183,7 @@
 vfoptions.ngridinterp=20; % 20 evenly-spaced points between each pair of consecutive a_grid points
 simoptions.gridinterplayer=vfoptions.gridinterplayer; % grid interpolation layer must also be set in simoptions (because it changes Policy size/interpretation)
 simoptions.ngridinterp=vfoptions.ngridinterp;
-simoptions.optimize_nProbs=1;
-simoptions.verbose=1;
+
 
 %% Now, create the return function
 DiscountFactorParamNames={'beta','sj'};
@@ -236,6 +235,9 @@
 end
 Params.mewj=Params.mewj./sum(Params.mewj); % Normalize to one
 AgeWeightsParamNames={'mewj'}; % So VFI Toolkit knows which parameter is the mass of agents of each age
+
+simoptions.optimize_nProbs=1;
+simoptions.verbose=2;
 StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightsParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
 
 

From 052661d9ea06fddaa60a6dcc3fd9a5944539dc67 Mon Sep 17 00:00:00 2001
From: Michael Tiemann <mdtiemann@gmail.com>
Date: Mon, 27 Jul 2026 12:07:10 +1200
Subject: [PATCH 5/5] Improve grid limits and plots

---
 Models20to30/LifeCycleModel28.m    |  4 ++--
 Models31to35/LifeCycleModel31.m    |  2 +-
 Models31to35/LifeCycleModel35.m    | 10 +++++-----
 Models40to44/LifeCycleModel41.m    |  8 ++++----
 Models40to44/LifeCycleModel42.m    |  4 ++--
 ModelsAppendixA/LifeCycleModelA6.m |  2 +-
 6 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/Models20to30/LifeCycleModel28.m b/Models20to30/LifeCycleModel28.m
index 89c4b4ae..e263eb5e 100644
--- a/Models20to30/LifeCycleModel28.m
+++ b/Models20to30/LifeCycleModel28.m
@@ -132,7 +132,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=12*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % First, the AR(1) process z
 [z_grid,pi_z]=discretizeAR1_FarmerToda(0,Params.rho_z,Params.sigma_epsilon_z,n_z);
@@ -253,7 +253,7 @@
 title('Life Cycle Profile: Fraction Time Worked (h)')
 subplot(5,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z h)')
-subplot(5,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(5,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(5,1,4); plot(1:1:Params.J,AgeConditionalStats.ninfants.Mean)
 title('Life Cycle Profile: Number of Infants (n1)')
diff --git a/Models31to35/LifeCycleModel31.m b/Models31to35/LifeCycleModel31.m
index 6933a9aa..ec24cf28 100644
--- a/Models31to35/LifeCycleModel31.m
+++ b/Models31to35/LifeCycleModel31.m
@@ -296,7 +296,7 @@
 title('Life Cycle Profile: Share of savings invested in risky asset (riskyshare)')
 subplot(3,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z)')
-subplot(3,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(3,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 
 
diff --git a/Models31to35/LifeCycleModel35.m b/Models31to35/LifeCycleModel35.m
index 076fc593..0cfabf74 100644
--- a/Models31to35/LifeCycleModel35.m
+++ b/Models31to35/LifeCycleModel35.m
@@ -34,8 +34,8 @@
 Params.J=100-Params.agejshifter; % =81, Number of period in life-cycle
 
 % Grid sizes to use
-n_d=[51,301]; % Decisions: riskyshare, savings
-n_a=[5,301]; % Endogenous housing and asset holdings
+n_d=[51,201]; % Decisions: riskyshare, savings
+n_a=[5,201]; % Endogenous housing and asset holdings
 n_z=7; % Exogenous labor productivity units shock
 n_u=5; % Between period i.i.d. shock
 N_j=Params.J; % Number of periods in finite horizon
@@ -121,8 +121,8 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=-3+30*(linspace(0,1,n_a(2)))'; % Note, I use equal spacing (normally would put most points near zero)
-% note: will go from -3 to 30-3
+asset_grid=-3+23*(linspace(0,1,n_a(2)))'; % Note, I use equal spacing (normally would put most points near zero)
+% note: will go from -3 to 23-3
 % Make it so that there is a zero assets
 % Find closest to zero assets
 [~,zeroassetindex]=min(abs(asset_grid));
@@ -324,7 +324,7 @@
 title('Life Cycle Profile: Share of savings invested in risky asset (riskyshare)')
 subplot(4,1,2); plot(1:1:Params.J,AgeConditionalStats.earnings.Mean)
 title('Life Cycle Profile: Labor Earnings (w kappa_j z)')
-subplot(4,1,3); plot(1:1:Params.J,AgeConditionalStats.assets.Mean)
+subplot(4,1,3); hold on; plot(1:1:Params.J,AgeConditionalStats.assets.Mean); plot(1:1:Params.J,AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(4,1,4); plot(1:1:Params.J,AgeConditionalStats.housing.Mean)
 title('Life Cycle Profile: Housing (h)')
diff --git a/Models40to44/LifeCycleModel41.m b/Models40to44/LifeCycleModel41.m
index 60b996ea..ca300f08 100644
--- a/Models40to44/LifeCycleModel41.m
+++ b/Models40to44/LifeCycleModel41.m
@@ -30,7 +30,7 @@
 
 % Grid sizes to use
 n_d=2; % Endogenous labour choice (working or not-working)
-n_a=[201,31]; % Endogenous asset holdings, female-labor-force-participation
+n_a=[201,21]; % Endogenous asset holdings, female-labor-force-participation
 n_z=11; % Exogenous labor productivity units shock
 N_j=Params.J; % Number of periods in finite horizon
 
@@ -94,9 +94,9 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=10*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+asset_grid=30*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
-h_grid=linspace(0.1,4,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
+h_grid=linspace(0.1,2,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
 % Note: deliberately omit 0 from h_grid
 
 % First, the AR(1) process z
@@ -253,7 +253,7 @@
 title('Life Cycle Profile: Female Earnings (w h z p)')
 subplot(5,1,3); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.maleearnings.Mean)
 title('Life Cycle Profile: Male Earnings (y_m)')
-subplot(5,1,4); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean)
+subplot(5,1,4); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 subplot(5,1,5); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.femaleLFPH.Maximum); hold off
 title('Life Cycle Profile: female labor force participation history (h)')
diff --git a/Models40to44/LifeCycleModel42.m b/Models40to44/LifeCycleModel42.m
index 1bcca56d..16e7d775 100644
--- a/Models40to44/LifeCycleModel42.m
+++ b/Models40to44/LifeCycleModel42.m
@@ -102,7 +102,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-asset_grid=10*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+asset_grid=30*(linspace(0,1,n_a(1)).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 h_grid=linspace(1,20,n_a(2))'; % Because h is an experienceasset, it will be interpolated onto this grid and so we need less grid points than usual
 % Note: deliberately omit 0 from h_grid
@@ -304,7 +304,7 @@
 title('Life Cycle Profile: Earnings (w h z l)')
 subplot(5,1,4); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.humancapital.Maximum); hold off
 title('Life Cycle Profile: Human Capital (h)')
-subplot(5,1,5); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean)
+subplot(5,1,5); hold on; plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Mean); plot(Params.agejshifter+(1:1:Params.J),AgeConditionalStats.assets.Maximum); hold off
 title('Life Cycle Profile: Assets (a)')
 % Because the calibration is a bit silly, assets are an inferior investment
 % to human capital. Hence throughout working life households just invest in
diff --git a/ModelsAppendixA/LifeCycleModelA6.m b/ModelsAppendixA/LifeCycleModelA6.m
index ece4142d..deb527bd 100644
--- a/ModelsAppendixA/LifeCycleModelA6.m
+++ b/ModelsAppendixA/LifeCycleModelA6.m
@@ -139,7 +139,7 @@
 %% Grids
 % The ^3 means that there are more points near 0 and near 10. We know from theory that the value function will be more 'curved' near zero assets,
 % and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
-a_grid=16*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
+a_grid=10*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
 
 % Grid for labour choice
 h_grid=linspace(0,1,n_d)'; % Notice that it is imposing the 0<=h<=1 condition implicitly