Skip to content

Commit

Permalink
Merge pull request #1342 from wdu9/master
Browse files Browse the repository at this point in the history
Fix Jacobian Methods
  • Loading branch information
Mv77 authored Aug 17, 2023
2 parents 4f644bf + 0cdecbd commit b18c71f
Show file tree
Hide file tree
Showing 4 changed files with 92 additions and 163 deletions.
2 changes: 1 addition & 1 deletion Documentation/CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ Release Date: TBD
- Adds option `sim_common_Rrisky` to control whether risky-asset models draw common or idiosyncratic returns in simulation. [#1250](https://github.com/econ-ark/HARK/pull/1250),[#1253](https://github.com/econ-ark/HARK/pull/1253)
- Addresses [#1255](https://github.com/econ-ark/HARK/issues/1255). Makes age-varying stochastic returns possible and draws from their discretized version. [#1262](https://github.com/econ-ark/HARK/pull/1262)
- Fixes bug in the metric that compares dictionaries with the same keys. [#1260](https://github.com/econ-ark/HARK/pull/1260)

- Fixes bug in the calc_jacobian method. [#1342](https://github.com/econ-ark/HARK/pull/1342)
### 0.13.0

Release Date: February, 16, 2023
Expand Down
103 changes: 78 additions & 25 deletions HARK/ConsumptionSaving/ConsIndShockModel.py
Original file line number Diff line number Diff line change
Expand Up @@ -2932,38 +2932,91 @@ def calc_jacobian(self, shk_param, T):
########
# STEP4 # of the algorithm
########

# Function to compute jacobian matrix from fake news matrix
def J_from_F(F):
J = F.copy()
for t in range(1, F.shape[0]):
J[1:, t] += J[:-1, t-1]
return J

J_A = J_from_F(Curl_F_A)
J_C = J_from_F(Curl_F_C)

########
# Additional step due to compute Zeroth Column of the Jacobian
########

params = deepcopy(self.__dict__["parameters"])
params["T_cycle"] = 2 # Dimension of Jacobian Matrix

params["LivPrb"] = params["T_cycle"] * [self.LivPrb[0]]
params["PermGroFac"] = params["T_cycle"] * [self.PermGroFac[0]]
params["PermShkStd"] = params["T_cycle"] * [self.PermShkStd[0]]
params["TranShkStd"] = params["T_cycle"] * [self.TranShkStd[0]]
params["Rfree"] = params["T_cycle"] * [self.Rfree]
params["UnempPrb"] = params["T_cycle"] * [self.UnempPrb]
params["IncUnemp"] = params["T_cycle"] * [self.IncUnemp]
params['IncShkDstn'] = params['T_cycle']* [self.IncShkDstn[0]]
params['cFunc_terminal_'] = deepcopy(self.solution[0].cFunc)

# Create instance of a finite horizon agent for calculation of zeroth
ZerothColAgent = IndShockConsumerType(**params)
ZerothColAgent.cycles = 1 # required

# If parameter is in time invariant list then add it to time vary list
ZerothColAgent.del_from_time_inv(shk_param)
ZerothColAgent.add_to_time_vary(shk_param)

# Jacobian Matrices
J_A = np.zeros((T, T)) # Asset Jacobian
J_C = np.zeros((T, T)) # Consumption Jacobian
for t in range(T):
for s in range(T):
if (t == 0) or (s == 0):
J_A[t][s] = Curl_F_A[t][s]
J_C[t][s] = Curl_F_C[t][s]
else:
J_A[t][s] = J_A[t - 1][s - 1] + Curl_F_A[t][s]
J_C[t][s] = J_C[t - 1][s - 1] + Curl_F_C[t][s]
# Update income process if perturbed parameter enters the income shock distribution
ZerothColAgent.update_income_process()

# Zeroth Column of the Jacobian
dD_0_0 = np.dot(tranmat_t[-2] - tranmat_ss, D_ss)
# Solve
ZerothColAgent.solve()

# this condition is because some attributes are specified as lists while other as floats
if type(getattr(self, shk_param)) == list:
peturbed_list = (
[getattr(self, shk_param)[0] + dx]
+ (params["T_cycle"] - 1) * [getattr(self, shk_param)[0]]
) # Sequence of interest rates the agent faces
else:
peturbed_list = (
[getattr(self, shk_param) + dx]
+ (params["T_cycle"] - 1) * [getattr(self, shk_param)]
) # Sequence of interest rates the agent

setattr(ZerothColAgent, shk_param, peturbed_list) # Set attribute to agent

D_curl_0_0 = dD_0_0 / dx
# Use Harmenberg Neutral Measure
ZerothColAgent.neutral_measure = True
ZerothColAgent.update_income_process()

c_first_col_0 = []
a_first_col_0 = []
for i in range(params["T_cycle"]):
c_first_col_0.append(np.dot(exp_vecs_c[i], D_curl_0_0))
a_first_col_0.append(np.dot(exp_vecs_a[i], D_curl_0_0))
# Calculate Transition Matrices
ZerothColAgent.define_distribution_grid()
ZerothColAgent.calc_transition_matrix()

tranmat_t_zeroth_col = ZerothColAgent.tran_matrix
dstn_t_zeroth_col = self.vec_erg_dstn.T[0]

C_t_no_sim = np.zeros(T)
A_t_no_sim = np.zeros(T)

c_first_col_0 = np.array(c_first_col_0)
a_first_col_0 = np.array(a_first_col_0)
for i in range(T):
if i ==0:
dstn_t_zeroth_col = np.dot(tranmat_t_zeroth_col[i],dstn_t_zeroth_col)
else:
dstn_t_zeroth_col = np.dot(tranmat_ss,dstn_t_zeroth_col)

C_t_no_sim[i] = np.dot(self.cPol_Grid ,dstn_t_zeroth_col)
A_t_no_sim[i] = np.dot( self.aPol_Grid ,dstn_t_zeroth_col)

J_A.T[0] = (A_t_no_sim - self.A_ss)/dx
J_C.T[0] = (C_t_no_sim - self.C_ss)/dx

return J_C, J_A

# Fill zeroth column of jacobian matrix
J_A.T[0] = a_first_col_0
J_C.T[0] = c_first_col_0

return J_C, J_A

def make_euler_error_func(self, mMax=100, approx_inc_dstn=True):
"""
Expand Down
2 changes: 1 addition & 1 deletion HARK/ConsumptionSaving/tests/test_IndShockConsumerType.py
Original file line number Diff line number Diff line change
Expand Up @@ -921,4 +921,4 @@ def test_calc_jacobian(self):

self.assertAlmostEqual(CJAC_Perm.T[30][29], -0.06120, places=HARK_PRECISION)
self.assertAlmostEqual(CJAC_Perm.T[30][30], 0.05307, places=HARK_PRECISION)
self.assertAlmostEqual(CJAC_Perm.T[30][31], 0.04674, places=HARK_PRECISION)
self.assertAlmostEqual(CJAC_Perm.T[30][31], 0.04674, places=HARK_PRECISION)
148 changes: 12 additions & 136 deletions examples/ConsIndShockModel/IndShockConsumerType_Jacobian_Example.ipynb

Large diffs are not rendered by default.

0 comments on commit b18c71f

Please sign in to comment.