WIP: BARK: Bellman Abstract Representation toolKit.

[sbenthall]

This PR introduces a new module, stage, which contains the implementation of a class representing a general Bellman-form stage. This stage can be solved using a generic value iteration algorithm. Stages can be composed to form more complex problems. A notebook demonstrates how the stage module can be used to create and solve a Portfolio Consumption problem. The current PR includes some basic tests.

This is the start of a much larger project that addresses some points raised in several issues, such as #1116, #997, #620, among others. It is work in progress. However, the state of this PR is at a point where it would be great to get reviews.

• Tests for new functionality/models or Tests to reproduce the bug-fix in code.
• Updated documentation of features that add new functionality.
• Update CHANGELOG.md with major/minor changes.

[sbenthall]

Looks like my current implementation uses some Python 3.9 features that are available in a different form in 3.8.

https://stackoverflow.com/questions/59955751/abcmeta-object-is-not-subscriptable-when-trying-to-annotate-a-hash-variable

[alanlujan91]

Look at estimagic for alternative optimizers https://github.com/OpenSourceEconomics/estimagic

They also have constraints https://estimagic.readthedocs.io/en/stable/how_to_guides/optimization/how_to_specify_constraints.html

Perhaps the demo notebook should be reversed; start with growth, then allocation, then consumption

[codecov-commenter]

Codecov Report

Base: 73.62% // Head: 72.92% // Decreases project coverage by -0.69% ⚠️

Coverage data is based on head (d662ae5) compared to base (e6c936b).
Patch coverage: 79.09% of modified lines in pull request are covered.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1186      +/-   ##
==========================================
- Coverage   73.62%   72.92%   -0.70%     
==========================================
  Files          72       74       +2     
  Lines       11509    12023     +514     
==========================================
+ Hits         8473     8768     +295     
- Misses       3036     3255     +219     

Impacted Files	Coverage Δ
HARK/stage.py	73.41% <73.41%> (ø)
HARK/tests/test_stage.py	100.00% <100.00%> (ø)
HARK/ConsumptionSaving/ConsIndShockModel.py	79.47% <0.00%> (-5.78%)	⬇️
HARK/utilities.py	29.63% <0.00%> (-5.23%)	⬇️
HARK/distribution.py	82.47% <0.00%> (-3.75%)	⬇️
HARK/tests/test_distribution.py	100.00% <0.00%> (ø)
HARK/numba.py
HARK/numba_tools.py	18.97% <0.00%> (ø)
HARK/core.py	91.48% <0.00%> (+0.01%)	⬆️
...nsumptionSaving/tests/test_IndShockConsumerType.py	78.55% <0.00%> (+3.39%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[sbenthall]

If we use estimagic in one part of the library, we probably should use it throughout.
I just looked it up and found that estimagic is an OpenSourceEconomics project, which is on brand for us!
But the development team seems very small: https://github.com/OpenSourceEconomics/estimagic/graphs/contributors

I'd like to talk more during a meeting about what it would mean to tie HARK to these tools.

[sbenthall]

Things seem to be working now.

There is one remaining issue I'd like to discuss, which involves the case where a stage is defined with a pre-given value for pi* for some state value. This is used in the consumption problems for the point where m = 0, so c = 0, but because of a binding constraint. The 'optimal' utility in this case is negative infinity, and this complicates both numerical optimization and interpolation.

The question is: in generality, is it safe to dispense with the role of shocks in these cases? It is straightforward to say, "if (input state) x has value a, then fix pi*(a) = b. But in this Bellman stage framework, the policy function typically takes both input states and realized shocks: pi : X, K -> A.

I hope we can discuss this in a meeting.

[alanlujan91]

If we use estimagic in one part of the library, we probably should use it throughout. I just looked it up and found that estimagic is an OpenSourceEconomics project, which is on brand for us! But the development team seems very small: https://github.com/OpenSourceEconomics/estimagic/graphs/contributors

I'd like to talk more during a meeting about what it would mean to tie HARK to these tools.

I've meet their lead contributor, and they seem to be on their way to being part of NumFocus so I am sure they would be happy to see us use their tools and maybe even future collaboration.

[alanlujan91]

Things seem to be working now.

There is one remaining issue I'd like to discuss, which involves the case where a stage is defined with a pre-given value for pi* for some state value. This is used in the consumption problems for the point where m = 0, so c = 0, but because of a binding constraint. The 'optimal' utility in this case is negative infinity, and this complicates both numerical optimization and interpolation.

The question is: in generality, is it safe to dispense with the role of shocks in these cases? It is straightforward to say, "if (input state) x has value a, then fix pi*(a) = b. But in this Bellman stage framework, the policy function typically takes both input states and realized shocks: pi : X, K -> A.

I hope we can discuss this in a meeting.

I think it's safe to dispense with the shocks, as well as anything else in the dynamics. At $m=0, c=0, u(c) = - \infty, v(m) = - \infty, v'(m) = + \infty$, BUT $v^{-1}(m) \equiv u^{-1}(v(m)) = 0$ and $v'^{-1}(m) = 0$. This is important not for the $m=0$ case per se, but for the $m = \varepsilon$ small which falls between 0 and the next smallest point $m = \delta$ in our approximation. Interpolation between $- \infty$ and a finite point is undefined, but interpolation in the inverse space is well behaved.

So for $m = \delta$ small we use the appropriate dynamics and shocks to evaluate all variables, but for $m= 0$ (or a different borrowing constraint) we need a good anchor for interpolation when an agent is unlucky and happens to fall into low wealth $m = \varepsilon$ between 0 and $\delta$.

[sbenthall]

Thanks @alanlujan91 this is all very helpful.

I've been digging a bit deeper into this and have a new take it.

How it works in this PR.
It's possible to define a $\pi^*(\vec{x}_p, \cdot) = \vec{a}_p$ in order to set the optimal action at a specific state $\vec{x}_a$.
From this, it is easy to compute $v_x(0) = u(\pi^*(0)) + \beta E[v_y(m - c)] = -\infty$.
But this only works for the consumption stage.

Whatever $\pi^*(0))$ is in the allocation stage, $v_x(0)$ will just be $v_y(0, \alpha)$.

And during the growth stage, there are no actions, so setting $\pi^*(0))$ doesn't make much sense.

Does this mean that for generic Bellman stages, it should be possible to add points to the input value function $v_x(0) = -\infty$ directly?

[alanlujan91]

Thanks @alanlujan91 this is all very helpful.

I've been digging a bit deeper into this and have a new take it.

How it works in this PR. It's possible to define a π∗(x→p,⋅)=a→p in order to set the optimal action at a specific state x→a. From this, it is easy to compute vx(0)=u(π∗(0))+βE[vy(m−c)]=−∞. But this only works for the consumption stage.

Whatever π∗(0)) is in the allocation stage, vx(0) will just be vy(0,α).

Correct. Presumably, $v_y(0, \alpha) = - \infty \ \forall \alpha \in [0, 1]$. This would probably also have to be added manually.

And during the growth stage, there are no actions, so setting π∗(0)) doesn't make much sense.

Does this mean that for generic Bellman stages, it should be possible to add points to the input value function vx(0)=−∞ directly?

I think it would be useful to be able to add points to the input value function. This however could be a generic treatment of the borrowing constraint. Perhaps not true universally, but for a few models if we set BorroCnstArt = None then the natural borrowing constraint will generally have this property of tending toward $- \infty$.

[sbenthall]

Ok, I am adding in the capability to add arbitrary points to both $\pi^*$ and $v_x$ for any stage.
This will come in the next few commits.

I am now wondering how to handle this for the growth stage, which takes both $a$ and $alpha$ as inputs.
I am tempted to say as you've suggested that $v_x^{\text{growth}}(0, \alpha) = - \infty \ \forall \alpha \in [0, 1]$.

However... the transition function for the growth stage is of the form $m = (f(\alpha, R, \eta)a + \theta) / (G\psi)$.
So while $v_y^{\text{growth}}(m=0) = -\infty$, starting the stage with $a = 0$ ends up with positive $m$ almost 100% of the time.

[alanlujan91]

Ok, I am adding in the capability to add arbitrary points to both π∗ and vx for any stage. This will come in the next few commits.

I am now wondering how to handle this for the growth stage, which takes both a and alpha as inputs. I am tempted to say as you've suggested that vxgrowth(0,α)=−∞ ∀α∈[0,1].

However... the transition function for the growth stage is of the form m=(f(α,R,η)a+θ)/(Gψ). So while vygrowth(m=0)=−∞, starting the stage with a=0 ends up with positive m almost 100% of the time.

You are correct. It depends on the exact income process; specifically on $\underline{\theta}$. when we add the possibility of unemployment with no UI, $\underline{\theta} = 0$ so there is some likelihood to end with $m=0$ next period, otherwise as you say the agent will end up with positive $m$.

This is why the calculation for the natural borrowing constraint is recursive. Assume last period; minimum $\underline{m}_T = 0$. So, in the second to last period $R a + \theta > \underline{m}_T$ which implies $\underline{a} = ( \underline{m}_T - \underline{\theta})/R$. If $\underline{\theta} = 0$ then $\underline{a} = 0$.

So the treatment at the boundary depends on both the income process and the imposition of a borrowing constraint. In the Portfolio choice model, the artificial borrowing constraint is 0, so it only depends on $\underline{\theta}$.

[sbenthall]

OK, thank you.

Since I am trying to build magic solver technology, I think I should be expecting the natural borrowing constraint to show up in the output of the model. But the artificial borrowing constraint must be an input to the model.

I do understand that much of the (different kind of) magic of the SolvingMicroDSOPs methodology is to use analytical model results to improve the computational solution.

What I have not yet figured out is if there is a way to use an analytically derived natural borrowing constraint to improve a magic solver.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[sbenthall]

I would like to revisit this PR, which has been on hold for several months. I thought I would add some context, as well as some ideas for next steps, here, to inform any review. Moving forward, I'm going to refer to the whole approach here, based on Bellman blocks, as BARK, to disambiguate it from the many other HARK redesign proposals floating around.

This PR is original research supported by NSF #2105301 and the Econ-ARK project. The goals of my efforts here have been to:

• Create a system for defining and solving problems that has feature parity with the scope of @llorracc's SolvingMicroDSOPs procedure, but not built around a specific problem. Rather, the user can configure the system model definition and solve it using general algorithms.
• Demonstrate a way to solve a number of outstanding issues in HARK, and resolve many ongoing frustrations with the original design, including consider explicit representation of Model state/action spaces #1116, age-varying transition equations? #997, define distributions directly in initial parameters for models #620,
• Build a functioning tool that can support future research, including my own research moving forward.

The core ideas represented in this PR are that:

• dynamic stochastic control problems can be broken up into smaller components that can be solved in sequence. These components are in Bellman form, can be composed like building blocks, and solved using backwards induction.
• these blocks can be defined as a library of components and reused in different problems, to enable easier comparison between structural models
• the structures/block can vary over the lifecycle of a finite horizon agent.
• the configuration of the model, the algorithms uses to solve models, and additional information needed by the solver (such as discretization grids) are all logically distinct.

I have been using the term 'stage' for the most basic of these blocks. I now regret this choice. Moving forward, I will take a line out of Pablo's book. There are two structural forms in BARK:

• b-blocks, representing a single Bellman update. I was calling these 'stages' before, but @llorracc seems to want to use stage to mean something else he has in mind. So be it.
• p-blocks, which can be composed of b-blocks and finish with a tick, which is an indication that the model should increment in time. (This is a convention borrowed from NetLogo, where it works quite well.)

Renaming of variables to better fit pre-defined blocks together is accomplished through a simple twist operation.

This PR includes an almost complete implementation of BARK. I believe my next step before putting this aside a few months ago was to verify that it can produce results equivalent to HARK for the ConsIndShockModel.

A few other notes:

• There has been some notional friction around the variable names used to define a b-block. I'm not wedded to any notation and the friction seems to be due to how b-blocks are defined at an entirely different level of abstraction than the models currently described in HARK.
• The current implementation supports the use of EGM for each b-block.
• The current implementation has the limitation that all shocks in b-blocks are IID.
• I have not been trying to optimize the Python code in this version. I have been focusing on the interface and form. Where possible, I've used attractive Python libraries like xarray. I know there's a lot to do to improve performance, such as: replacing for-loops with matrix operations; JITing code, and so on. I wanted to wait to work on this until I had automated tests in place for the main functionality.

The main reason I put BARK aside earlier is that @llorracc has been keen on the next HARK version being integrated with Dolang and I was aiming to support that effort. However, in recent conversations, it sounds like @albop is not convinced that Dolang is necessary to support @llorracc 's use case, and proposed building a pure Python API for structural lifecycle models instead. In that case, I propose that BARK might be a solution to @llorracc 's problems. I would be immensely grateful for any feedback from @albop on the BARK design, if only for my own edification. In any case, I don't have any immediate plans to integrate BARK with Dolang; I don't have a need for Dolang, personally -- I'm fine with defining models in software code, as long as it can be done modularly.

I see a few paths forward for BARK and I'm wondering which one is best:

• Prepare BARK as an incremental improvement to HARK, replacing the current FramedAgentType, to providing an alternative approach to building and solving models that can be compared with HARK's more hand-crafted models. I think this would take no more than a month of work to wrap up, given its current scope.
• Split off BARK into a separate repository that builds on HARK.

I currently don't see any tools out there that are better than BARK for what I need for my research, and I've recently gotten excited about ways to extend BARK to solve asset pricing problems (this would involve defining new kinds of blocks). I'd be happy to pick up a different tool if somebody pointed me to it.

[sbenthall]

I'll add that one reason prefer the 'block' terminology is that it lines up with Bence's work on the Sequence Space Jacobian. While the current implementation of BARK uses only a linear sequence of blocks, I could imagine a system where blocks form a DAG for the individual agent's problem. I don't want to add to the scope of this PR at all, which is minimal. I'm just here suggesting where BARK would go, if it evolves where I think it should go.

[llorracc]

@sbenthall thanks for this cogent post.

@dedwar65 is right now working on a revision to SolvingMicroDSOPs, with the particular ambition to remove all references to the antique Mathematica code.

If BARK is ready to be used to implement all of the elements in the SolvingMicroDSOPs notebook, using it to do that would seem to be a natural next step for demonstrating how BARK can be used in practice. Having you and Decory work together on that in the near term seems like a good path forward, since he is more familiar with the math and economics and you have more expertise on coding.

[sbenthall]

I'd be happy to work with @dedwar65 on that.

[codecov]

Codecov Report

Patch coverage: 47.11% and project coverage change: -0.96 ⚠️

Comparison is base (0b108f5) 73.43% compared to head (5aa1d97) 72.47%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1186      +/-   ##
==========================================
- Coverage   73.43%   72.47%   -0.96%     
==========================================
  Files          77       79       +2     
  Lines       12845    13331     +486     
==========================================
+ Hits         9433     9662     +229     
- Misses       3412     3669     +257     

Impacted Files	Coverage Δ
HARK/stage.py	41.45% <41.45%> (ø)
HARK/tests/test_stage.py	100.00% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[sbenthall]

This PR has been superseded by more recent work on 'blocks', which have made it into the repository.

I still draw on this work a great deal in my more recent implementations, but this code won't make it into HARK, so closing the PR.