Conversation
Got a bee in my bonnet, cranked this out quickly based on notes I wrote for DCL. Completely untested because I didn't make an AgentType subclass for it, but all the math for the solver is there.
Kind of embarrassing typo.
I added the remaining bits to the solver and added a skeletal AgentType. It can initialize and begin to run (solves T-1), but there's some problem with the lower boundary of the state space that I'm still working out. I'm sure it's something terribly dumb.
Problem at the bottom of the state space was due to tiny (1e-16) violations of the lower boundary, easily fixed. Consumption function is still "wiggly" from the perspective of (m,n) space. Might need to implement vPPfunc and cubic spline interpolation, and *also* might need to conduct a more precise search for the bounds of the region of inaction, taking expectations repeatedly. Setting aside for now.
Deleting method I added to interpolation.py during debugging.
Made these at NumFocus summit?
Alan asked me about this, so I'm checking on its status.
There was a problem hiding this comment.
Pull Request Overview
This PR introduces a work-in-progress basic consumption-saving model with liquid and illiquid assets, and reduces the resolution of shock discretizations in the Euler error function.
- Reduced the number of quadrature points (N) from 200 to 100 for both Transitory and Permanent shocks.
- Halved the tail_N parameter from 50 to 25 in both shock discretizations.
Comments suppressed due to low confidence (1)
HARK/ConsumptionSaving/ConsIndShockModel.py:2214
- Add targeted tests to verify that the reduced shock discretization still delivers acceptable Euler error accuracy under typical model settings.
N=100,
| N=100, | ||
| method="equiprobable", | ||
| tail_N=50, | ||
| tail_N=25, |
There was a problem hiding this comment.
[nitpick] Consider documenting the rationale for reducing N from 200 to 100 (and tail_N from 50 to 25), or expose these as configurable parameters so that users can adjust discretization resolution as needed.
There was a problem hiding this comment.
These are probably just temporary changes I made while working on the model. Or they were supposed to be in a different branch that I was working on at the same time-- this PR doesn't do anything with the Euler error calculator.
|
bugbot run |
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 2 out of 2 changed files in this pull request and generated 21 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| alpha = (0.0 - x0) / (x1 - x0) | ||
| OptZeroDeposit[j] = (1.0 - alpha) * k0 + alpha * k1 | ||
|
|
||
| idx = np.searchsorted(LogNvrsRatio_temp[:cut], RatioTarg) | ||
| k0 = kNrm_temp[idx - 1] | ||
| k1 = kNrm_temp[idx] | ||
| x0 = LogNvrsRatio_temp[idx - 1] | ||
| x1 = LogNvrsRatio_temp[idx] | ||
| # print(cut,x0,x1,k0,k1) | ||
| alpha = (RatioTarg - x0) / (x1 - x0) | ||
| OptZeroWithdraw[j] = (1.0 - alpha) * k0 + alpha * k1 |
There was a problem hiding this comment.
Potential division by zero on lines 658 and 667. If x1 equals x0 (i.e., the log ratio is constant between two grid points), the division will fail or produce Inf. Add a check to ensure x1 != x0 before dividing.
| for j in range(bCount): | ||
| plt.plot(kNrmNow[:, j], cNrmNow[:, j]) | ||
| plt.show() | ||
|
|
There was a problem hiding this comment.
This debugging plot should be removed before merging. It interrupts the solver execution and is not appropriate for production code.
| for j in range(bCount): | |
| plt.plot(kNrmNow[:, j], cNrmNow[:, j]) | |
| plt.show() |
| Marginal value of middle-of-period liquid kash-on-hand, as a function of | ||
| middle-of-period kash-on-hand kNrm and middle-of-period illiquid assets bNrm. | ||
| dvdb : function | ||
| Marginal value of middle-of-period illiquid assets, as a function of | ||
| middle-of-period kash-on-hand kNrm and middle-of-period illiquid assets bNrm. |
There was a problem hiding this comment.
Inconsistent naming: 'kash' appears to be intentionally spelled differently from 'cash', but the usage is inconsistent. The variable is called 'kNrm' throughout but described as 'kash-on-hand' in comments. Consider standardizing on either 'cash' or providing a clear explanation for the 'k' notation in the module docstring.
| Marginal value of middle-of-period liquid kash-on-hand, as a function of | |
| middle-of-period kash-on-hand kNrm and middle-of-period illiquid assets bNrm. | |
| dvdb : function | |
| Marginal value of middle-of-period illiquid assets, as a function of | |
| middle-of-period kash-on-hand kNrm and middle-of-period illiquid assets bNrm. | |
| Marginal value of middle-of-period liquid cash-on-hand, as a function of | |
| middle-of-period cash-on-hand kNrm and middle-of-period illiquid assets bNrm. | |
| dvdb : function | |
| Marginal value of middle-of-period illiquid assets, as a function of | |
| middle-of-period cash-on-hand kNrm and middle-of-period illiquid assets bNrm. |
|
|
||
| if __name__ == "__main__": | ||
| MyType = BasicIlliquidConsumerType() | ||
| MyType.cycles = 3 | ||
| MyType.solve() | ||
|
|
||
| B = 6.0 | ||
| f = lambda x: MyType.MPCfuncSimple(0, x, B * np.ones_like(x)) | ||
| g = lambda x: MyType.cFuncSimple(0, x, B * np.ones_like(x)) | ||
| h = lambda x: MyType.solution[0].dFunc(x, B * np.ones_like(x)) | ||
| z = lambda x: (x + MyType.solution[0].hNrm + B) * MyType.solution[0].MPCmin | ||
| c = lambda x: MyType.solution[0].cFunc(x, B * np.ones_like(x)) | ||
|
|
||
| plot_funcs([g], -7.5, 12) |
There was a problem hiding this comment.
This temporary test code in the main block should be removed before merging. Test code belongs in a proper test file in the tests directory, not in the module's main block.
| if __name__ == "__main__": | |
| MyType = BasicIlliquidConsumerType() | |
| MyType.cycles = 3 | |
| MyType.solve() | |
| B = 6.0 | |
| f = lambda x: MyType.MPCfuncSimple(0, x, B * np.ones_like(x)) | |
| g = lambda x: MyType.cFuncSimple(0, x, B * np.ones_like(x)) | |
| h = lambda x: MyType.solution[0].dFunc(x, B * np.ones_like(x)) | |
| z = lambda x: (x + MyType.solution[0].hNrm + B) * MyType.solution[0].MPCmin | |
| c = lambda x: MyType.solution[0].cFunc(x, B * np.ones_like(x)) | |
| plot_funcs([g], -7.5, 12) |
| A consumer type that faces idiosyncratic shocks to income and has a different | ||
| interest factor on saving vs borrowing. Extends IndShockConsumerType, with | ||
| very small changes. Solver for this class is currently only compatible with | ||
| linear spline interpolation. | ||
|
|
||
| Same parameters as AgentType. | ||
|
|
||
|
|
||
| Parameters | ||
| ---------- |
There was a problem hiding this comment.
The class docstring is incomplete and inaccurate. It says "Same parameters as AgentType" but doesn't document the specific parameters for this model (e.g., Rilqd, IlqdPenalty, bNrmGrid). The docstring also says it extends IndShockConsumerType when it actually extends KinkedRconsumerType.
| A consumer type that faces idiosyncratic shocks to income and has a different | |
| interest factor on saving vs borrowing. Extends IndShockConsumerType, with | |
| very small changes. Solver for this class is currently only compatible with | |
| linear spline interpolation. | |
| Same parameters as AgentType. | |
| Parameters | |
| ---------- | |
| A consumer type with two financial assets: a liquid asset and an illiquid | |
| asset. The illiquid asset earns a (potentially) higher return than the liquid | |
| asset but is costly to access: withdrawals from the illiquid asset incur a | |
| proportional penalty. Income is subject to idiosyncratic shocks, and the | |
| liquid asset faces different interest factors for saving vs borrowing, as in | |
| :class:`KinkedRconsumerType`, which this class extends. | |
| This type uses all parameters of :class:`KinkedRconsumerType`, plus | |
| additional parameters governing the illiquid asset. | |
| Parameters | |
| ---------- | |
| Rilqd : float | |
| Gross return factor on the illiquid asset between periods. | |
| IlqdPenalty : float | |
| Proportional penalty rate applied to withdrawals from the illiquid asset. | |
| A value of zero implies costless access; higher values make withdrawals | |
| more expensive. | |
| bNrmGrid : array_like | |
| Grid for normalized balances in the illiquid asset used by the solver. |
| def derivativeX(self, mNrm, nNrm, eps=1e-8): | ||
| dLo = self.__call__(mNrm - eps, nNrm) | ||
| dHi = self.__call__(mNrm + eps, nNrm) | ||
| return (dHi - dLo) / (2 * eps) | ||
|
|
||
| def derivativeY(self, mNrm, nNrm, eps=1e-8): |
There was a problem hiding this comment.
Using numerical differentiation with finite differences in derivative methods can be numerically unstable and inaccurate, especially near discontinuities in the DepositFunction (at the boundaries of the inaction region). Consider implementing analytical derivatives or at least document the potential numerical issues and the choice of eps=1e-8.
| def derivativeX(self, mNrm, nNrm, eps=1e-8): | |
| dLo = self.__call__(mNrm - eps, nNrm) | |
| dHi = self.__call__(mNrm + eps, nNrm) | |
| return (dHi - dLo) / (2 * eps) | |
| def derivativeY(self, mNrm, nNrm, eps=1e-8): | |
| def derivativeX(self, mNrm, nNrm, eps=1e-8): | |
| """ | |
| Numerical derivative of the optimal net deposit with respect to liquid | |
| market resources mNrm, computed via a central finite difference. | |
| Parameters | |
| ---------- | |
| mNrm : float or array_like | |
| Liquid market resources. | |
| nNrm : float or array_like | |
| Illiquid market resources. | |
| eps : float, optional | |
| Step size used for the central difference approximation. The default | |
| value eps=1e-8 is chosen as a small perturbation relative to typical | |
| state-space scales, but callers may wish to adjust it depending on | |
| the application. | |
| Notes | |
| ----- | |
| This method approximates the derivative using: | |
| (d(mNrm + eps, nNrm) - d(mNrm - eps, nNrm)) / (2 * eps), | |
| where d(·, ·) is the `DepositFunction` itself. Because the deposit policy | |
| generally has kinks/discontinuities at the boundaries of the inaction | |
| region, this numerical derivative can be noisy or unstable near those | |
| boundaries and should be interpreted as an approximation only. | |
| """ | |
| dLo = self.__call__(mNrm - eps, nNrm) | |
| dHi = self.__call__(mNrm + eps, nNrm) | |
| return (dHi - dLo) / (2 * eps) | |
| def derivativeY(self, mNrm, nNrm, eps=1e-8): | |
| """ | |
| Numerical derivative of the optimal net deposit with respect to illiquid | |
| market resources nNrm, computed via a central finite difference. | |
| Parameters | |
| ---------- | |
| mNrm : float or array_like | |
| Liquid market resources. | |
| nNrm : float or array_like | |
| Illiquid market resources. | |
| eps : float, optional | |
| Step size used for the central difference approximation. The default | |
| value eps=1e-8 is chosen as a small perturbation relative to typical | |
| state-space scales, but callers may wish to adjust it depending on | |
| the application. | |
| Notes | |
| ----- | |
| This method approximates the derivative using: | |
| (d(mNrm, nNrm + eps) - d(mNrm, nNrm - eps)) / (2 * eps), | |
| where d(·, ·) is the `DepositFunction` itself. As with `derivativeX`, | |
| the presence of kinks/discontinuities at inaction-region boundaries | |
| can cause this numerical derivative to be inaccurate or unstable near | |
| those points, so it should be treated as an approximate marginal effect. | |
| """ |
| # Define the functions to be searched for FOC solutions | ||
| # dvdkFunc_temp = MargValueFuncCRRA(cFunc_by_bNrm_list[j], CRRA) | ||
| # dvdbFunc_temp = MargValueFuncCRRA(dvdbNvrsFunc_by_bNrm_list[j], CRRA) | ||
| # LogRatioFunc = lambda x: np.log(dvdkFunc_temp(x) / dvdbFunc_temp(x)) | ||
| # LogRatioFuncAlt = lambda x: LogRatioFunc(x) - RatioTarg | ||
|
|
||
| # # Define the bounds of the search in k | ||
| # if j < 10: | ||
| # kBotD = kNrmNow[0, j] - BoroCnstNat[j] | ||
| # kTopD = kBotD + 3.0 #kNrmNow[-1, j] - BoroCnstNat[j] | ||
| # kBotW = kNrmNow[0, j] - BoroCnstNat[j] | ||
| # kTopW = kNrmNow[-1, j] - BoroCnstNat[j] | ||
| # else: | ||
| # kBotD = OptZeroDeposit[j - 1] - BoroCnstNat[j - 1] | ||
| # kTopD = kBotD + 1.0 | ||
| # kBotW = OptZeroWithdraw[j - 1] - BoroCnstNat[j - 1] | ||
|
|
||
| # # Perform a bounded search for optimal zero withdrawal and deposit | ||
| # print(kBotD, kTopD, LogRatioFunc(kBotD), LogRatioFunc(kTopD)) | ||
| # #plot_funcs(LogRatioFunc, kNrmNow[0, j] - BoroCnstNat[j], kNrmNow[-1, j] - BoroCnstNat[j]) | ||
| # OptZeroDeposit[j] = ( | ||
| # root_scalar( | ||
| # LogRatioFunc, | ||
| # bracket=(kBotD, kTopD), | ||
| # xtol=1e-8, | ||
| # rtol=1e-8, | ||
| # method="brentq", | ||
| # ).root | ||
| # + BoroCnstNat[j] | ||
| # ) |
There was a problem hiding this comment.
This comment appears to contain commented-out code.
| # Define the functions to be searched for FOC solutions | |
| # dvdkFunc_temp = MargValueFuncCRRA(cFunc_by_bNrm_list[j], CRRA) | |
| # dvdbFunc_temp = MargValueFuncCRRA(dvdbNvrsFunc_by_bNrm_list[j], CRRA) | |
| # LogRatioFunc = lambda x: np.log(dvdkFunc_temp(x) / dvdbFunc_temp(x)) | |
| # LogRatioFuncAlt = lambda x: LogRatioFunc(x) - RatioTarg | |
| # # Define the bounds of the search in k | |
| # if j < 10: | |
| # kBotD = kNrmNow[0, j] - BoroCnstNat[j] | |
| # kTopD = kBotD + 3.0 #kNrmNow[-1, j] - BoroCnstNat[j] | |
| # kBotW = kNrmNow[0, j] - BoroCnstNat[j] | |
| # kTopW = kNrmNow[-1, j] - BoroCnstNat[j] | |
| # else: | |
| # kBotD = OptZeroDeposit[j - 1] - BoroCnstNat[j - 1] | |
| # kTopD = kBotD + 1.0 | |
| # kBotW = OptZeroWithdraw[j - 1] - BoroCnstNat[j - 1] | |
| # # Perform a bounded search for optimal zero withdrawal and deposit | |
| # print(kBotD, kTopD, LogRatioFunc(kBotD), LogRatioFunc(kTopD)) | |
| # #plot_funcs(LogRatioFunc, kNrmNow[0, j] - BoroCnstNat[j], kNrmNow[-1, j] - BoroCnstNat[j]) | |
| # OptZeroDeposit[j] = ( | |
| # root_scalar( | |
| # LogRatioFunc, | |
| # bracket=(kBotD, kTopD), | |
| # xtol=1e-8, | |
| # rtol=1e-8, | |
| # method="brentq", | |
| # ).root | |
| # + BoroCnstNat[j] | |
| # ) | |
| # Legacy search-based implementation for determining zero-withdrawal and | |
| # zero-deposit points in k has been removed; current behavior is defined | |
| # by the active code in this method. |
| # LogRatioFunc, | ||
| # bracket=(kBotD, kTopD), | ||
| # xtol=1e-8, | ||
| # rtol=1e-8, | ||
| # method="brentq", | ||
| # ).root | ||
| # + BoroCnstNat[j] | ||
| # ) | ||
| # # print(bNrmGrid[j], 'deposit', OptZeroDeposit[j]) | ||
| # if j > 1: | ||
| # kTopW = OptZeroDeposit[j] - BoroCnstNat[j] | ||
| # try: | ||
| # OptZeroWithdraw[j] = ( | ||
| # root_scalar( | ||
| # LogRatioFuncAlt, | ||
| # bracket=(kBotW, kTopW), | ||
| # xtol=1e-8, | ||
| # rtol=1e-8, | ||
| # method="brentq", | ||
| # ).root | ||
| # + BoroCnstNat[j] | ||
| # ) | ||
| # except: | ||
| # OptZeroWithdraw[j] = BoroCnstNat[j] | ||
| # print(bNrmGrid[j], "deposit", OptZeroDeposit[j], "withdraw", OptZeroWithdraw[j]) | ||
|
|
There was a problem hiding this comment.
This comment appears to contain commented-out code.
| # LogRatioFunc, | |
| # bracket=(kBotD, kTopD), | |
| # xtol=1e-8, | |
| # rtol=1e-8, | |
| # method="brentq", | |
| # ).root | |
| # + BoroCnstNat[j] | |
| # ) | |
| # # print(bNrmGrid[j], 'deposit', OptZeroDeposit[j]) | |
| # if j > 1: | |
| # kTopW = OptZeroDeposit[j] - BoroCnstNat[j] | |
| # try: | |
| # OptZeroWithdraw[j] = ( | |
| # root_scalar( | |
| # LogRatioFuncAlt, | |
| # bracket=(kBotW, kTopW), | |
| # xtol=1e-8, | |
| # rtol=1e-8, | |
| # method="brentq", | |
| # ).root | |
| # + BoroCnstNat[j] | |
| # ) | |
| # except: | |
| # OptZeroWithdraw[j] = BoroCnstNat[j] | |
| # print(bNrmGrid[j], "deposit", OptZeroDeposit[j], "withdraw", OptZeroWithdraw[j]) |
| # + BoroCnstNat[j] | ||
| # ) | ||
| # # print(bNrmGrid[j], 'deposit', OptZeroDeposit[j]) | ||
| # if j > 1: | ||
| # kTopW = OptZeroDeposit[j] - BoroCnstNat[j] | ||
| # try: | ||
| # OptZeroWithdraw[j] = ( | ||
| # root_scalar( | ||
| # LogRatioFuncAlt, | ||
| # bracket=(kBotW, kTopW), | ||
| # xtol=1e-8, | ||
| # rtol=1e-8, | ||
| # method="brentq", | ||
| # ).root | ||
| # + BoroCnstNat[j] | ||
| # ) | ||
| # except: | ||
| # OptZeroWithdraw[j] = BoroCnstNat[j] | ||
| # print(bNrmGrid[j], "deposit", OptZeroDeposit[j], "withdraw", OptZeroWithdraw[j]) | ||
|
|
There was a problem hiding this comment.
This comment appears to contain commented-out code.
| # + BoroCnstNat[j] | |
| # ) | |
| # # print(bNrmGrid[j], 'deposit', OptZeroDeposit[j]) | |
| # if j > 1: | |
| # kTopW = OptZeroDeposit[j] - BoroCnstNat[j] | |
| # try: | |
| # OptZeroWithdraw[j] = ( | |
| # root_scalar( | |
| # LogRatioFuncAlt, | |
| # bracket=(kBotW, kTopW), | |
| # xtol=1e-8, | |
| # rtol=1e-8, | |
| # method="brentq", | |
| # ).root | |
| # + BoroCnstNat[j] | |
| # ) | |
| # except: | |
| # OptZeroWithdraw[j] = BoroCnstNat[j] | |
| # print(bNrmGrid[j], "deposit", OptZeroDeposit[j], "withdraw", OptZeroWithdraw[j]) |
| def update(self): | ||
| super().update() |
There was a problem hiding this comment.
This method requires 1 positional argument, whereas overridden Model.update may be called with arbitrarily many. This call correctly calls the base method, but does not match the signature of the overriding method.
| def update(self): | |
| super().update() | |
| def update(self, *args, **kwargs): | |
| super().update(*args, **kwargs) |
This is another new (to HARK) micro model: a basic consumption-saving problem with a liquid and an illiquid asset. I started work on this back in the fall, then forgot about it. It's definitely not complete, but I don't fully recall how incomplete it is. It looks like the solver mostly works, but has bug(s) left-- the temporary test code only has it run for 3 periods, and I left a bunch of debugging code in.
Definitely a work in progress, but figured I'd issue the PR because Alan asked about it.