Update Jacobian Methods For ConsMarkovType and IndShockConsumerType [WIP]#1248
Update Jacobian Methods For ConsMarkovType and IndShockConsumerType [WIP]#1248wdu9 wants to merge 7 commits intoecon-ark:mainfrom
Conversation
Codecov ReportPatch coverage:
Additional details and impacted files@@ Coverage Diff @@
## master #1248 +/- ##
==========================================
- Coverage 73.31% 72.20% -1.11%
==========================================
Files 76 76
Lines 12548 12751 +203
==========================================
+ Hits 9199 9207 +8
- Misses 3349 3544 +195
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 in Codecov by Sentry. |
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
Codecov ReportPatch coverage:
Additional details and impacted files@@ Coverage Diff @@
## master #1248 +/- ##
==========================================
- Coverage 72.55% 71.50% -1.06%
==========================================
Files 78 78
Lines 13009 13212 +203
==========================================
+ Hits 9439 9447 +8
- Misses 3570 3765 +195
☔ View full report in Codecov by Sentry. |
There was a problem hiding this comment.
Pull request overview
This PR introduces sequence-space Jacobian utilities and extends MarkovConsumerType with transition-matrix and steady-state calculations, alongside relaxing shape checks for Markov-specific interest rate inputs.
Changes:
- Added
Fake_News_JACinHARK/utilities.pyto compute asset and consumption Jacobians using the “fake news” algorithm, factoring out logic previously inIndShockConsumerType.calc_jacobian. - Extended
MarkovConsumerTypewithcalc_transition_matrix,calc_ergodic_dist, andcompute_steady_state, plus new imports for transition-matrix utilities and sparse eigen computations. - Relaxed
MarkovConsumerType.check_markov_inputsto allowRfreeto be provided as a list of arrays as well as a single NumPy array.
Reviewed changes
Copilot reviewed 2 out of 3 changed files in this pull request and generated 18 comments.
| File | Description |
|---|---|
| HARK/utilities.py | Adds the Fake_News_JAC helper implementing the fake-news Jacobian algorithm in a reusable form (with some docstring issues to fix). |
| HARK/ConsumptionSaving/ConsMarkovModel.py | Introduces Markov-specific transition-matrix, ergodic distribution, and steady-state computation logic and broadens Rfree validation, but contains several shape, indexing, and API-consistency problems that need correction and tests. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
|
|
||
|
|
||
| LivPrb = self.LivPrb[0][0] # Update probability of staying alive |
There was a problem hiding this comment.
In the infinite-horizon branch, the survival probability used in the transition matrices is hard-coded as self.LivPrb[0][0], i.e., the survival probability of the first Markov state only. Since self.LivPrb is validated above as a vector of length StateCount, this effectively ignores any heterogeneity in survival probabilities across Markov states; consider indexing by i when constructing the transition matrix if state-specific survival is intended.
| def compute_steady_state(self,IncShkDstn_ntrl_msr): | ||
|
|
||
|
|
||
|
|
||
|
|
||
| #solve the consumer's problem | ||
|
|
||
| self.solve() | ||
|
|
||
| self.define_distribution_grid(dist_pGrid = np.array([1])) | ||
| self.calc_transition_matrix(IncShkDstn_ntrl_msr) | ||
| self.calc_ergodic_dist() | ||
|
|
||
|
|
||
| C = 0 | ||
| A = 0 | ||
| for i in range(self.state_num): | ||
|
|
||
| C += self.prb_dstn[i]*np.dot(self.cPol_Grid[i],self.vec_erg_dstns[i]) | ||
| A += self.prb_dstn[i]*np.dot(self.aPol_Grid[i],self.vec_erg_dstns[i]) | ||
|
|
||
|
|
||
| self.A_ss = A | ||
| self.C_ss = C | ||
|
|
||
| return C , A |
There was a problem hiding this comment.
compute_steady_state for MarkovConsumerType returns (C, A) while the analogous compute_steady_state method for IndShockConsumerType returns (A_ss, C_ss). This inconsistency in return order between closely related types makes it easy to misinterpret results when swapping models; consider standardizing the return tuple (and documenting it) so both classes follow the same convention.
| Parameters | ||
| ---------- | ||
| transition_matrix: List | ||
| transition matrix whose ergordic distribution is to be solved |
There was a problem hiding this comment.
The docstring for calc_ergodic_dist describes the parameter as a "transition matrix whose ergordic distribution is to be solved"; "ergordic" is misspelled and should be "ergodic".
| transition matrix whose ergordic distribution is to be solved | |
| transition matrix whose ergodic distribution is to be solved |
| if isinstance(self.Rfree,list): | ||
|
|
||
| for rfree_t in self.Rfree: | ||
|
|
||
| # Check that arrays are the right shape | ||
| if not isinstance(rfree_t, np.ndarray) or rfree_t.shape != (StateCount,): | ||
| raise ValueError( | ||
| "Rfree not the right shape, it should an array of Rfree of all the states." | ||
| ) | ||
| else: | ||
|
|
||
| # Check that arrays are the right shape | ||
| if not isinstance(self.Rfree, np.ndarray) or self.Rfree.shape != (StateCount,): | ||
| raise ValueError( | ||
| "Rfree not the right shape, it should an array of Rfree of all the states." | ||
| ) |
There was a problem hiding this comment.
check_markov_inputs now accepts self.Rfree as a list of per-period arrays, but other methods (e.g., get_Rfree) still assume self.Rfree is a single NumPy array indexed directly by the Markov state. As a result, passing a list for Rfree will pass validation here but later fail or behave incorrectly when used, so either the rest of the class needs to be generalized to handle list-valued Rfree or this method should continue to enforce an array-only contract.
| if isinstance(self.Rfree,list): | |
| for rfree_t in self.Rfree: | |
| # Check that arrays are the right shape | |
| if not isinstance(rfree_t, np.ndarray) or rfree_t.shape != (StateCount,): | |
| raise ValueError( | |
| "Rfree not the right shape, it should an array of Rfree of all the states." | |
| ) | |
| else: | |
| # Check that arrays are the right shape | |
| if not isinstance(self.Rfree, np.ndarray) or self.Rfree.shape != (StateCount,): | |
| raise ValueError( | |
| "Rfree not the right shape, it should an array of Rfree of all the states." | |
| ) | |
| if not isinstance(self.Rfree, np.ndarray) or self.Rfree.shape != (StateCount,): | |
| raise ValueError( | |
| "Rfree not the right shape, it should be a NumPy array of Rfree for all the states." | |
| ) |
| if type(self.Rfree)==list: | ||
| b_next_i = self.Rfree[k][0]*a_next_i | ||
| bNext.append(b_next_i) | ||
| else: | ||
| b_next_i = self.Rfree*a_next_i | ||
| bNext.append(b_next_i) |
There was a problem hiding this comment.
In the finite-horizon branch of calc_transition_matrix, b_next_i is computed using either self.Rfree[k][0] * a_next_i or self.Rfree * a_next_i. This both ignores the dependence of the interest factor on the current Markov state (only the first entry is ever used from each period’s vector) and, when self.Rfree is a NumPy array of per-state rates, produces a 2D array that is inconsistent with the 1D bNext expected by gen_tran_matrix_1D/gen_tran_matrix_2D. The computation of b_next_i should use the scalar interest factor corresponding to state i in period k so that bNext[i] has shape matching a_next_i.
| if type(self.Rfree)==list: | |
| b_next_i = self.Rfree[k][0]*a_next_i | |
| bNext.append(b_next_i) | |
| else: | |
| b_next_i = self.Rfree*a_next_i | |
| bNext.append(b_next_i) | |
| if type(self.Rfree) == list: | |
| R_k = self.Rfree[k] | |
| # Allow for either scalar or per-state array at period k | |
| if np.ndim(R_k) == 0: | |
| R_ki = R_k | |
| else: | |
| R_ki = R_k[i] | |
| else: | |
| # Allow for scalar or 1-D array of per-state rates | |
| if np.ndim(self.Rfree) == 0: | |
| R_ki = self.Rfree | |
| else: | |
| R_ki = self.Rfree[i] | |
| b_next_i = R_ki * a_next_i | |
| bNext.append(b_next_i) |
|
|
||
|
|
||
|
|
||
| if shk_dstn == None: |
There was a problem hiding this comment.
Testing for None should use the 'is' operator.
| if shk_dstn == None: | |
| if shk_dstn is None: |
| elif self.T_cycle!= 0: | ||
|
|
||
| # for finite horizon, we can account for changing levels of prb_unemp because of endogenous job finding probability by imposing a list of these values, so for q'th period, the probability is slightly different | ||
| if shk_dstn == None: |
There was a problem hiding this comment.
Testing for None should use the 'is' operator.
| if shk_dstn == None: | |
| if shk_dstn is None: |
| c_first_col_0 = np.zeros(T) | ||
| a_first_col_0 = np.zeros(T) | ||
|
|
| c_first_col_0 = np.zeros(T) | ||
| a_first_col_0 = np.zeros(T) | ||
|
|
| def compute_steady_state(self,IncShkDstn_ntrl_msr): | ||
|
|
||
|
|
||
|
|
There was a problem hiding this comment.
This method requires 2 positional arguments, whereas overridden IndShockConsumerType.compute_steady_state requires 1.
| def compute_steady_state(self,IncShkDstn_ntrl_msr): | |
| def compute_steady_state(self, IncShkDstn_ntrl_msr=None): | |
| # Preserve previous behavior when the argument is omitted by | |
| # explicitly raising a TypeError, similar to Python's default message. | |
| if IncShkDstn_ntrl_msr is None: | |
| raise TypeError( | |
| "compute_steady_state() missing required positional argument " | |
| "'IncShkDstn_ntrl_msr'" | |
| ) |
Please ensure your pull request adheres to the following guidelines: