Dear Robert,
I am interested in a recent paper by Boar and Midrigan (2022), Efficient Redistribution, published on the Journal of Monetary Economics.
Having the replication of this paper in the toolkit would greatly help me since the part om transition is quite complicated to code.
Interesting! Just to let you know that you can find the replication code in Matlab for their AER Insights paper which uses a similar model. The code is available on Data and Code for: Should We Tax Capital Income or Wealth?
Thanks @aledinola
I tried to use their code but it is based on a toolkit by Miranda and Fackler that requires C mex and it does not work on my computer
At a brief glance it seems do-able. I will read the paper properly at the end of the week and see how it goes.
Hi @aledinola, I stand corrected: it turns out that you can use the files of Miranda and Fackler even without mex, just that they are very slow. They double coded everything in plain matlab
Thanks for clarifying this point, @javier_fernan. I wanted to add that you can certainly run the âplain matlabâ files as you mentioned, but maybe you can also just run the .mex files without the need to recompile them. MF provide the executables on their website. It depends on your OS.
That said, I would strongly recommend using the VFI Toolkit instead. The MF package is still useful, but it has become somewhat outdated. In my experience, the mex approach is conceptually appealing (Iâve used with Fortran), but it requires dealing with many low-level system details (for example, ensuring that MATLAB can properly detect and use your C compiler). By contrast, the VFI Toolkit is much more âbatteries includedâ and works smoothly right out of the box.
A good starting point for Boar and Midrigan might be my repo on GitHub - aledinola/PijoanMasTaxes: Aiyagari model with endogenous labor (as in Pijoan-Mas 2006), with VFI toolkit
Here are some codes implementing the core exercise of Boar & Midrigan (2022).
This shows how to solve for the initial stationary general eqm, then for one specific tax policy reform it solves for the final stationary general eqm and the transition path, as well as then evaluating the social welfare function.
The full paper of Boar & Midrigan (2022) is essentially to repeat this for lots of different tax policy reforms to find the optimal one (defined as the tax policy reform that maximizes the social welfare function).
Be warned, the codes are currently very slow (I did the bare minimum to solve). But there are no end of ways to speed them up, and I describe a few here.
i) the initial stationary general eqm could be sped up, but you only ever want to solve it once, so would be a waste of your time to do so. Instead, I would add a binary variable near start of code SkipInitialGeneralEqm=1 and then around the code that solves the initial eqm I would put
if SkipInitialGeneralEqm=1
existing code that solves the initial stationary general eqm
save BM2022initeqm.mat
else
load BM2022initeqm.mat
end
so you just run once with =1, and then you can set it to =0 from then on.
ii) the final stationary eqm could be massively sped up by changing to heteroagentoptions.fminalgo=5 or fminalgo=8 (5 requires you to specify update rules, see the Appendix to the Intro to OLG Models; 8 is Matlabâs lsqnonlin(), the only reason it is not the default is that it requires Matlabâs Optimization Toolbox). You could also massively speed up with a good guess. If I was doing the BM2022 full exercise, I would first put a 10x10 grid on tau and tau_a (the main tax rates, as the paper points out the xi and xi_a that control progressivity are secondary in the model results), and then I would use the âclosestâ result from this grid as my initial guess for the general eqm prices while trying to find the optimal tax policy.
iii) The transition path is currently using a bad initial guess, so takes way way longer. As described in second half of (ii) I would repeat this same process to generate initial guesses for the transition PricePath0. You should also play with the âdampingâ hyperparameters that determine how big the updates are (bigger is faster, but too big and it becomes unstable; a good initial guess for PricePath0 reduces the chance of instability and allows you to set larger values for the damping hyperparameters).
You can also play around with the size of n_d, n_a, and vfoptions.ngridinterp (and for the transition path with vfoptions.level1n). And maybe adjust
heteroagentoptions.toleranceGEprices=10^(-4); % Accuracy of general eqm prices
heteroagentoptions.toleranceGEcondns=10^(-4); % Accuracy of general eqm eqns
which control the tolerance for the convergence of the stationary eqm, and for the tolerance on the transition path
transpathoptions.tolerance=10^(-5);
transpathoptions.updateaccuracycutoff=10^(-9); % If the suggested update is less than this then don't bother; 10^(-9) is decent odds to be numerical error anyway (currently only works for transpathoptions.GEnewprice=3)
As always the easiest way to speed them up is just to find a more powerful GPU 
More generally, I need to implement some alternative algorithms for the transition path. VFI Toolkit currently only has shooting algorithm. BM2022 say that they started out using a Jacobian-based update [If I implement Matlabâs lsqnonlin() it will be a much better version of this that also uses Hessian], but that they switched to Anderson Type-I acceleration [Jacobian-based, but you keep a memory of the Jacobianâs from the last few updates, and update based on an average of the current and lagged Jacobians, rather than just the current Jacobian]
Whether or not all of this is enough to get you fast enough for a full replication I am not sure. Will certainly be close. Unfortunately BM2022 makes no mention of their hardware or their runtimes. They used equivalent of n_a=500. They used same n_z setup, but they discretize with Rouwenhourst and I used Farmer-Toda (which is known to perform better for the specific AR(1) in this model).
I think that the only thing in the paper I donât presently do is to calculate the âchange in welfareâ (just calculate the level). But given how they calculate welfare this is probably trivial. I may or may not come back and do this later.
If anyone wants to submit an improved version of the code that does some of the above, please do and I will replace it. Also, I did not do much testing, so possible there is a typo.
One comment on the BM2022 paper itself.
They use a social welfare function based on âan approach simliar to that of Benabou (2002) and Bakis, Kaymak & Poschke (2015)â. This massively simplifies the calculation of the social welfare, or more accurately of the calculation of the compensating equivalent variation which is normally heinous in a model with consumption and leisure (see, e.g., the calculations in Conesa, Kitao & Krueger (2009) for what happens if you do the compensating equivalent variation in full). What I am not sure of is if this simplification is in any way misleading? Does anyone out there know of a paper that compares the two in a model with endogenous labor? It seems like a calculation that will give much the same answer in a model where there is only consumption, but this is not clear once you also have endogenous labor. Iâd be interested so see results on if the difference matters or not.
After writing the first paragraph in this post I went and did a bit more reading of Benabou (2002). It is not really a SWF. It is to evaluate âaggregate efficiencyâ. This makes sense when you do the Delta=0 version. As he describes it in the intro, âTo evaluate more generally the extent to which market distortions and imperfections are worsened or improved by alternatives policies, the paper proposes a new measure of overall economic efficiency. This criterion properly reflects (dynamic) variations in the aggregate consumption of goods and leisure and in the idiosyncratic risks that agents face; but, unlike a standard social welfare function, it does not reward interpersonal equality per se. The underlying idea is straightforward. Instead of aggregating individual incomes or consumptions (which washes out all idiosyncratic uncertainty), or individual utilities (which introduces a bias towards egalitarian allocations), one sums up consumption certainty-equivalents, so as to obtain a kind of risk-adjusted GDP.â
This makes sense of the Delta=0 version, where it is about âaggregate efficiencyâ, rather than SWF. But I still find the Delta=1 version to be odd, BM2022 call theirs âutilitarianâ, but it is not really (as that should be the full calculation like in CKK2009). As Benabou (2002) puts it in his Section 5 about âA Criterion of Aggregate Economic Efficiencyâ âMore generally, when aggregating over individuals the relevant parameter should be societyâs degree of inequality-aversionâ. The BM2022 exercise with Delta=1 seems better understood as an easy-to-compute back-of-envelope approximation of utilitarian.
Looking at BKP2015, they actually do the full utilitarian calculation (their eqn 6), and point out that once you have consumption and leisure, even just getting the certainty-equivalent for use in the aggregate efficiency formula requires making a decision about how exactly to calculate it (their eqn 7 and the paragraph before it). Note that calculating welfare (and hence the optimal tax) is no harder with the full utilitarian calculation than with the SWF-approximation of BM2022, it is only calculating the compensating equivalent variation that is hard with the full utilitarian calculation.
This looks quite interesting, Iâll have a look and see if I can further improve it performance-wise.
I have noticed that grid interp layer with d variable is quite slow in infinite horizon, see VFI in InfHorz: Howards, PFI, and relation to implicit finite differences - #35 by aledinola
[quote=âaledinola, post:35, topic:408â]
As far as I understand, with refinement we precompute F(d,a',a,z) and then maximize with respect to d to obtain F^\*(a',a,z) and d^\*(a',a,z). But we need the optimal d for all possible points on the fine grid over a'. This means that F will be a huge matrix. Even when using lowmemory, where we loop over z, we still have to store F(d,aâ,a) which is big.
[/quote]
A possible solution would be not precomputing F(d,aâ,a,z) for all the points on the fine grid aâ, i.e. not do refinement. So for all combinations of a and z, evaluate the RHS of the Bellman equation on (d,aâ) where aâ belongs to the original small grid. Find the maximum and only then build the fine grid and recalculate the RHS of the Bellman equation for each (d,aâ) where now aâ in on the fine grid. In this way you avoid precomputing the big array F(d,aâ,a,z) before the proper value function iteration.
I donât know if this makes sense 
Furthermore, I can try to run the original code of BM on my computer and report the runtimes. But I am not sure if I can make the C-mex work: I might run into the same problems encountered by @javier_fernan
Yeah, I think you are right that there is definitely some better way to implement the interaction between ârefineâ and âgrid interpolation layerâ. The challenge is just finding what the right combo is while still always giving the correct solution (without having to assume anything beyond conditional monotonicity).
Yes, agree: a better way of combining ârefineâ with âgridinterplayerâ would potentially speed up the computation of the steady state.
I tried to suggest an improvement for âgridinterplayerâ in this post, but it does not consider any d variable. Still, maybe the idea is useful: Faster Aiyagari example with interpolation
The performance of the transition is another issue, though. There we do the value function by backward induction over time, like in a life-cycle model, and we donât use refinement (please correct me if I am wrong). Hence refinement if not responsible for the slow performance of the transition in infinite horizon models.
Refinement is not useful in finite horizon problems (including transition paths in an infinite horizon model) because you would only ever use the ârefinedâ solution once, so it is a waste to refine it to being with.
[Exception is portfolio-choice with endogenous labor, where there is one decision variable that appears in return fn but not expectations, and another decision variable that appears in expectation but not return fn, so by refining before combining you avoid having to consider the interactions.]
In any case, the tookit is pretty good in finite horizon thanks to the combination of divide-and-conquer with grid interpolation layer. The BM2022 example code does a single iteration of the transition path in a decent runtime.
Maybe I should think about doing a âvfoptions.maxaprimediffâ based on markov chain approximation approach for InfHorz problems. You would only consider aprime indexes that are the current a index ± vfoption.maxaprimediff; markov chain approximation, including implicity finite difference, is to use vfoptions.maxaprimediff=1, but no reason it cannot be set to any integer, which will also allow you to use it in more situations, and to still work with larger n_a grids. This could be done together with grid interpolation layer and refine.
Is what you call âmarkov chain approximation approach for InfHorz problemsâ equivalent to the method used by Pontus R. in the discrete time model in his JEDC paper?
As far as I understand, he iterates on the value function using linear interpolation and Newtonâs method. To compute the derivative of the value function, he uses finite differences.
But the choice of aâ is not restricted to anything: if you see his code discrete_time.m, the choice variable is consumption. Once he has consumption, he computes the policy functions for aprime, called ape and apu in the code. Then he does the Howard acceleration like you do in the gridinterplayer code, just that you use interp1 and he computes manually the left grid point and the correspinding weight
Yeah. But when is the âleft grid point and corresponding weightâ a valid way to implement linear interpolation? Answer: when you are moving at most one grid point (down)! [You need the left for âdownwindâ and the right for âupwindâ]. So this will only work to do the linear interpolation if you move at most one grid point down (or up). If you move more than one grid point the linear interpolation is incorrect, and we know that value fn is concave, so we know this will be incorrect (not just that it might be).
If the true solution involved moving at most two grid points, the error introduced is probably pretty minor, but if you move five the errors likely start getting substantial.
You can kind of see this in Section 3.3.1 of Rendahl (2022) when he shows a matrix âDâ that is the continuous time implicit finite differences algorithm he is explaining how to implement a discrete time equivalent of. Notice how âDâ is only non-zero for up/down one point.
PS. The value fn is generally concave (it inherits it from the return fn, which inherits it from the utility fn which is concave), but not at a few points like when a binding constraint switches to non-binding, or with non-convex choice sets. Theorems establishing this are somewhere in SLP1989 (from memory).
PPS. There are two other important difference between grid interpolation layer and implict finite differences (or markov approximation methods in general) The first is that while both linearly interpolates the expectations, the grid interpolation layer explicitly evaluates the return fn while the later linearly interpolates the return fn (which also means the former cannot violate constraints, while the later can unless you check for them elsewhere). The second is that the former is still stuck on a âfine gridâ while the later is not.
Thank you Robert, your codes run on my computer! The steady-state is slow but it works. Havenât checked the transition yet
Hi @aledinola, I tried your function VFI_interp_fast and it is indeed fast, with interpolation! It took me a bit to make it work because the main file calls many other functions that I had to comment out.
One question: I was able to run Boar and Midrigan with fixed labor supply using your function. Can I have variable labor?
Hi @javier_fernan thanks for your interest in my repo. I will polish it so that it will be easy to run Aiyagari with my fast interpolation method.
Unfortunately you cannot have variable labor out of the box since I coded it for a model with (aâ,a,z) and not (d,aâ,a,z).
Thanks for the very nice explanation Robert. I will have to read seriously these papers at some point