diff --git a/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb b/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb index c9de4c1a8..5d24aa7c8 100644 --- a/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb +++ b/examples/ConsNewKeynesianModel/KS-HARK-presentation.ipynb @@ -60,6 +60,20 @@ "from HARK.ConsumptionSaving.ConsNewKeynesianModel import NewKeynesianConsumerType" ] }, + { + "cell_type": "markdown", + "id": "28d24bc5", + "metadata": {}, + "source": [ + "# Calibration and setup\n", + "\n", + "This notebook uses the HARK toolkit to solve [Krusell and Smith (1998)](https://www.journals.uchicago.edu/doi/abs/10.1086/250034), applying the [Sequence Space Jacobian](https://onlinelibrary.wiley.com/doi/abs/10.3982/ECTA17434) and [SSJ toolkit](https://github.com/shade-econ/sequence-jacobian). \n", + "\n", + "## Firm setup\n", + "\n", + "Collect calibration for the production economy in a dictionary." + ] + }, { "cell_type": "code", "execution_count": 3, @@ -75,43 +89,23 @@ "outputs": [], "source": [ "calibration = {\n", - " \"eis\": 1,\n", - " \"delta\": 0.025,\n", - " \"alpha\": 0.11,\n", - " \"L\": 1.0,\n", - " \"K\": 1.0,\n", - " \"Y\": 1.0,\n", - " \"r\": 0.01,\n", + " \"eis\": 1, # Elasticity of intertemporal substitution\n", + " \"delta\": 0.025, # Depreciation rate\n", + " \"alpha\": 0.11, # Capital share of income\n", + " \"L_ss\": 1.0, # Steady state labor\n", + " \"Y_ss\": 1.0, # Steady state output\n", + " \"r_ss\": 0.01, # Steady state real interest rate\n", "}" ] }, - { - "cell_type": "code", - "execution_count": 4, - "id": "2a97701f", - "metadata": { - "execution": { - "iopub.execute_input": "2024-07-11T15:18:16.176291Z", - "iopub.status.busy": "2024-07-11T15:18:16.176051Z", - "iopub.status.idle": "2024-07-11T15:18:16.178499Z", - "shell.execute_reply": "2024-07-11T15:18:16.178018Z" - } - }, - "outputs": [], - "source": [ - "L_ss = 1.0 # Steady state labor\n", - "r_ss = 0.01 # steady state interest rate\n", - "Y_ss = 1.0 # steady state output" - ] - }, { "cell_type": "markdown", "id": "1bac99bf", "metadata": {}, "source": [ - "## \n", + "### Steady State Capital \n", "\n", - "Given these steady state choices, we will need to find $K, Z$ to clear the firm first order conditions.\n" + "Find steady state capital (`K_ss`) and productivity (`Z`) implied by the firm's first order condition, aggregate resource constraint and the above steady state values of labor (`L_ss`), the real interest rate (`r_ss`) and output (`Y_ss`). " ] }, { @@ -132,15 +126,18 @@ "\n", "\n", "def your_funcs(X):\n", - " L = calibration[\"L\"]\n", + " L_ss = calibration[\"L_ss\"]\n", " alpha = calibration[\"alpha\"]\n", " delta = calibration[\"delta\"]\n", + " r_ss = calibration[\"r_ss\"]\n", + " Y_ss = calibration[\"Y_ss\"]\n", "\n", " K, Z = X\n", - " # all RHS have to be 0\n", + "\n", + " # Firms first order condition and aggregate resource constraint\n", " f = [\n", - " alpha * Z * (K / L) ** (alpha - 1) - delta - r_ss, # r = MPK\n", - " Z * K**alpha * L ** (1 - alpha) - Y_ss, # Y = Z*F(K,L)\n", + " alpha * Z * (K / L_ss) ** (alpha - 1) - delta - r_ss, # r_ss = MPK\n", + " Z * K**alpha * L_ss ** (1 - alpha) - Y_ss, # Y = Z*F(K,L)\n", " ]\n", "\n", " return f\n", @@ -148,7 +145,10 @@ "\n", "sol = root(your_funcs, [1.0, 1.0]) # find roots\n", "\n", - "K_ss, Z_ss = sol.x" + "K_ss, Z_ss = sol.x\n", + "\n", + "calibration[\"K_ss\"] = K_ss\n", + "calibration[\"Z_ss\"] = Z_ss" ] }, { @@ -191,7 +191,7 @@ "id": "922dce5d", "metadata": {}, "source": [ - " Let's double check the roots we find produce our chosen steady state values for $ r, Y , L$." + "Double check the roots we find produce our chosen steady state values. " ] }, { @@ -211,17 +211,19 @@ "def firm(\n", " K,\n", " Z,\n", - " L=calibration[\"L\"],\n", + " L_ss=calibration[\"L_ss\"],\n", " alpha=calibration[\"alpha\"],\n", " delta=calibration[\"delta\"],\n", "):\n", - " r = alpha * Z * (K / L) ** (alpha - 1) - delta\n", - " w = (1 - alpha) * Z * (K / L) ** alpha\n", - " Y = Z * K**alpha * L ** (1 - alpha)\n", + " r = alpha * Z * (K / L_ss) ** (alpha - 1) - delta\n", + " w = (1 - alpha) * Z * (K / L_ss) ** alpha\n", + " Y = Z * K**alpha * L_ss ** (1 - alpha)\n", " return r, w, Y\n", "\n", "\n", - "r_ss, w_ss, Y_ss = firm(sol.x[0], sol.x[1])" + "r_ss, w_ss, Y_ss = firm(sol.x[0], sol.x[1])\n", + "\n", + "calibration[\"w_ss\"] = w_ss" ] }, { @@ -246,7 +248,30 @@ } ], "source": [ - "print(r_ss, w_ss, Y_ss)" + "print(\"--------------------------------------+------------\")\n", + "print(f\"Steady state capital | {K_ss:.4f}\")\n", + "print(f\"Steady state output | {Y_ss:.4f}\")\n", + "print(f\"Steady state real interest rate | {r_ss:.4f}\")\n", + "print(f\"Steady state wage | {w_ss:.4f}\")\n", + "print(\"--------------------------------------+------------\")" + ] + }, + { + "cell_type": "markdown", + "id": "7ca88d16", + "metadata": {}, + "source": [ + "## HARK agent\n", + "\n", + "HARK represents agent problems where aggregate variables can affect an individual income process as instances of `NewKeynesianConsumerType`, a subclass of `AgentType` (see [A Gentle Introduction to HARK](https://docs.econ-ark.org/examples/Gentle-Intro/Gentle-Intro-To-HARK.html)).\n", + "\n", + "Instances of the `NewKeynesianConsumerType` class contain functions that solve policy functions of optimizing agents and compute Jacobian matrices in response to a policy shocks using the [SSJ toolkit](https://github.com/shade-econ/sequence-jacobian). \n", + "\n", + "**_NOTE:_** `NewKeynesianConsumerType` is still a microeconomic agent solving an income fluctuation problem. The only difference from [`IndShockConsumerType`](https://docs.econ-ark.org/examples/HowWeSolveIndShockConsumerType/HowWeSolveIndShockConsumerType.html) is that additional aggregate labour income variables are passed to the agent's income process. This allows the agent's problem to be defined within a general equilibrium framework such as (but not restricted to) HANK. \n", + "\n", + "**_NOTE:_** Researchers can also create their own agent types by subclassing `AgentType`. See [HARK's documentation](https://docs.econ-ark.org/examples/HowWeSolveIndShockConsumerType/HowWeSolveIndShockConsumerType.html) for more information.\n", + "\n", + "To create an instance of `NewKeynesianConsumerType`, first specify parameters in a dictionary. To start, we use a discount factor of 0.98." ] }, { @@ -263,40 +288,42 @@ }, "outputs": [], "source": [ + "L_ss = calibration[\"L_ss\"]\n", + "w_ss = calibration[\"w_ss\"]\n", + "r_ss = calibration[\"r_ss\"]\n", + "\n", "HANK_Dict = {\n", - " # Parameters shared with the perfect foresight model\n", + " # Individual agent 'preferences' (shared with perfect foresight model)\n", " \"CRRA\": calibration[\"eis\"], # Coefficient of relative risk aversion\n", - " \"Rfree\": (1 + r_ss), # Interest factor on assets\n", " \"DiscFac\": 0.98, # Intertemporal discount factor\n", " \"LivPrb\": [0.99375], # Survival probability\n", + " # Individual lifcycle income process parameters\n", " \"PermGroFac\": [1.00], # Permanent income growth factor\n", - " # Parameters that specify the income distribution over the lifecycle\n", - " # Standard deviation of log permanent shocks to income\n", - " \"PermShkStd\": [0.06],\n", + " \"PermShkStd\": [0.06], # Standard deviation of log permanent shocks to income\n", " \"PermShkCount\": 5, # Number of points in discrete approximation to permanent income shocks\n", - " # Standard deviation of log transitory shocks to income\n", - " \"TranShkStd\": [0.2],\n", - " \"TranShkCount\": 5,\n", - " # HANK params\n", - " \"tax_rate\": [\n", - " 0,\n", - " ], # set to 0.0 because we are going to assume that labor here is actually after tax income\n", - " \"labor\": [L_ss],\n", - " \"wage\": [w_ss],\n", - " # Number of points in discrete approximation to transitory income shocks\n", + " \"TranShkStd\": [0.2], # Standard deviation of log transitory shocks to income\n", + " \"TranShkCount\": 5, # Number of points in discrete approximation to transitory income shocks\n", + " # Parameters related to unemployment and retirement\n", " \"UnempPrb\": 0.0, # Probability of unemployment while working\n", " \"IncUnemp\": 0.0, # Unemployment benefits replacement rate\n", " \"UnempPrbRet\": 0.0000, # Probability of \"unemployment\" while retired\n", " \"IncUnempRet\": 0.0, # \"Unemployment\" benefits when retired\n", " \"T_retire\": 0.0, # Period of retirement (0 --> no retirement)\n", - " # Parameters for constructing the \"assets above minimum\" grid\n", + " # Aggregates affecting the agent's decision\n", + " \"Rfree\": (1 + r_ss), # Interest factor for assets faced by agents\n", + " \"wage\": [w_ss], # Wage rate faced by agents\n", + " \"tax_rate\": [\n", + " 0\n", + " ], # set to 0.0 because we are going to assume that labor here is actually after tax income\n", + " \"labor\": [L_ss], # Aggregate (mean) labor supply\n", + " # Parameters for constructing \"assets above minimum\" grid\n", " \"aXtraMin\": 0.0001, # Minimum end-of-period \"assets above minimum\" value\n", " \"aXtraMax\": 2000, # Maximum end-of-period \"assets above minimum\" value\n", " \"aXtraCount\": 200, # Number of points in the base grid of \"assets above minimum\"\n", " # Exponential nesting factor when constructing \"assets above minimum\" grid\n", " \"aXtraNestFac\": 3,\n", " \"aXtraExtra\": None, # Additional values to add to aXtraGrid\n", - " # Transition Matrix simulation parameters\n", + " # Transition matrix simulation parameters\n", " \"mCount\": 200,\n", " \"mMax\": 2000,\n", " \"mMin\": 0.0001,\n", @@ -309,7 +336,7 @@ "id": "3be9593e", "metadata": {}, "source": [ - "# Create HARK agent" + "Let's create a temporary instance of a `NewKeynesianConsumerType` agent." ] }, { @@ -326,7 +353,7 @@ }, "outputs": [], "source": [ - "Agent = NewKeynesianConsumerType(**HANK_Dict)" + "tempAgent = NewKeynesianConsumerType(**HANK_Dict)" ] }, { @@ -334,10 +361,28 @@ "id": "fa266888", "metadata": {}, "source": [ - "# Find Steady state discount factor clear asset market\n", + "Given a discount factor, the agent will supply a steady state capital stock `A_ss`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ebe1792e", + "metadata": {}, + "outputs": [], + "source": [ + "A_ss = tempAgent.compute_steady_state()[0]\n", + "print(f\"Steady state agent asset supply for beta = 0.98 is: {A_ss:.3f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "911805ab", + "metadata": {}, + "source": [ + "### Steady State Assets \n", "\n", - "We will estimate the discount factor to ensure that asset supply equals the steady state capital we found earlier. \n", - "\n" + "Since we are interested in computing a steady state equilibrium, we find `discFac` such that `A_ss` clears the asset market given that `K_ss` is the steady state firm capital demand. " ] }, { @@ -365,7 +410,10 @@ "def A_ss_func(beta):\n", " HANK_Dict[\"DiscFac\"] = beta\n", "\n", + " # For a given value of beta, we need to re-create the Agent\n", " Agent_func = NewKeynesianConsumerType(**HANK_Dict, verbose=False)\n", + "\n", + " # And then solve for the steady state supply\n", " A_ss = Agent_func.compute_steady_state()[0]\n", "\n", " return A_ss\n", @@ -380,7 +428,15 @@ "\n", "DiscFac = optimize.brentq(ss_dif, 0.8, 0.9999)\n", "\n", - "print(\"Time taken to solve for steady state\", time.time() - start)" + "print(f\"Time taken to solve for steady state {time.time() - start:.3f} secs.\")" + ] + }, + { + "cell_type": "markdown", + "id": "01be6a8a", + "metadata": {}, + "source": [ + "We can now create the steady state agent using the general equilibrium discount factor." ] }, { @@ -399,8 +455,8 @@ "source": [ "# Create a new agent\n", "HANK_Dict[\"DiscFac\"] = DiscFac\n", - "Agent_GE = NewKeynesianConsumerType(**HANK_Dict, verbose=False)\n", - "A_ss, C_ss = Agent_GE.compute_steady_state()" + "steadyHANK = NewKeynesianConsumerType(**HANK_Dict, verbose=False)\n", + "A_ss, C_ss = steadyHANK.compute_steady_state()" ] }, { @@ -426,9 +482,12 @@ } ], "source": [ - "# to make sure goods and asset markets clear\n", - "print(\"goods_clearing\", Y_ss - C_ss - calibration[\"delta\"] * K_ss)\n", - "print(\"asset_clearing\", A_ss - K_ss)" + "# To make sure goods and asset markets clear\n", + "print(\n", + " \"Final goods clearing:\",\n", + " calibration[\"Y_ss\"] - C_ss - calibration[\"delta\"] * calibration[\"K_ss\"],\n", + ")\n", + "print(\"Asset clearing:\", A_ss - calibration[\"K_ss\"])" ] }, { @@ -436,7 +495,15 @@ "id": "a8167383", "metadata": {}, "source": [ - "# Computing Heterogenous Agent Jacobians" + "# Computing Jacobians using SSJ\n", + "\n", + "With the steady state agent in hand, we can compute the Jacobians of the steady state HANK agent's policy with respect to the aggregate state variables.\n", + "\n", + "The `calc_jacobian` method of our `NewKeynesianConsumerType` agent computes the Jacobians of the steady state aggregate consumption (`CJACW`) and assets (`AJACW`) with respect to a pertubation of a specified variable.\n", + "\n", + "Recall `CJACW[s,t]` is the time $t$ response to a shock at time $s$. \n", + "\n", + "Here is an example where we shock the wage rate. " ] }, { @@ -463,9 +530,9 @@ "source": [ "start = time.time()\n", "\n", - "CJACW, AJACW = Agent_GE.calc_jacobian(\"wage\", 300) # Wage jacobians\n", + "CJACW, AJACW = steadyHANK.calc_jacobian(\"wage\", 300) # Wage jacobians\n", "\n", - "print(\"Time taken to compute jacobians\", time.time() - start)" + "print(f\"Time taken to compute wage Jacobians: {time.time() - start:.3f} seconds\")" ] }, { @@ -493,15 +560,16 @@ } ], "source": [ - "plt.plot(CJACW.T[0])\n", - "plt.plot(CJACW.T[20])\n", - "plt.plot(CJACW.T[50])\n", - "plt.plot(CJACW.T[100])\n", + "plt.plot(CJACW.T[0], label=\"s=0\")\n", + "plt.plot(CJACW.T[20], label=\"s=20\")\n", + "plt.plot(CJACW.T[50], label=\"s=50\")\n", + "plt.plot(CJACW.T[100], label=\"s=100\")\n", "plt.xlim(-2, 300)\n", "plt.plot(np.arange(300), np.zeros(300), color=\"k\")\n", - "plt.title(\"Consumption Jacobian Wage\")\n", - "plt.xlabel(\"quarters\")\n", - "plt.ylabel(\"C response\")\n", + "plt.title(\"Consumption Jacobian (Wage) \")\n", + "plt.xlabel(\"Quarters\")\n", + "plt.ylabel(\"Consumption response\")\n", + "plt.legend()\n", "plt.show()" ] }, @@ -529,9 +597,19 @@ "source": [ "start = time.time()\n", "\n", - "CJACR, AJACR = Agent_GE.calc_jacobian(\"Rfree\", 300) # Rfree jacobians\n", + "CJACR, AJACR = steadyHANK.calc_jacobian(\"Rfree\", 300) # Rfree jacobians\n", "\n", - "print(\"Time taken to compute jacobians\", time.time() - start)" + "print(\n", + " f\"Time taken to compute return factor Jacobians: {time.time() - start:.3f} seconds\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "311a3a60", + "metadata": {}, + "source": [ + "Here is an example where we shock the return factor. " ] }, { @@ -559,16 +637,25 @@ } ], "source": [ - "plt.plot(CJACR.T[0])\n", - "plt.plot(CJACR.T[30])\n", - "plt.plot(CJACR.T[50])\n", + "plt.plot(CJACR.T[0], label=\"s =0\")\n", + "plt.plot(CJACR.T[20], label=\"s=20\")\n", + "plt.plot(CJACR.T[50], label=\"s=50\")\n", + "plt.plot(CJACR.T[100], label=\"s=100\")\n", "plt.plot(np.arange(300), np.zeros(300), color=\"k\")\n", - "plt.title(\"Consumption Jacobian interest rate\")\n", - "plt.xlabel(\"quarters\")\n", - "plt.ylabel(\"C response\")\n", + "plt.title(\"Consumption Jacobian (Return Factor)\")\n", + "plt.xlabel(\"Quarters\")\n", + "plt.ylabel(\"Consumption response\")\n", "plt.show()" ] }, + { + "cell_type": "markdown", + "id": "c5d0b70e", + "metadata": {}, + "source": [ + "Let's store the Jacobians in a dictionary for later use." + ] + }, { "cell_type": "code", "execution_count": 18, @@ -621,7 +708,7 @@ "id": "d3669fab", "metadata": {}, "source": [ - "## Other Blocks of the Model" + "## Impulse Response Functions and Simulations " ] }, { @@ -675,7 +762,7 @@ "id": "f3d569b5", "metadata": {}, "source": [ - "# Solving for Impulse Responses" + "### Solving for Impulse Responses" ] }, { @@ -692,8 +779,8 @@ }, "outputs": [], "source": [ - "T = 300 # <-- the length of the IRF\n", - "rho_Z = 0.8 # persistence of IRF shock\n", + "T = 300 # Length of the IRF\n", + "rho_Z = 0.8 # Persistence of IRF shock\n", "dZ = 0.001 * Z_ss * rho_Z ** np.arange(T)\n", "shocks = {\"Z\": dZ}\n", "\n", @@ -701,7 +788,6 @@ "unknowns = [\"K\"]\n", "targets = [\"asset_mkt\"]\n", "\n", - "\n", "irfs_Z = ks.solve_impulse_linear(SteadyState_Dict, unknowns, targets, shocks)" ] }, @@ -723,7 +809,7 @@ " irfs_list,\n", " variables,\n", " labels=[\" \"],\n", - " ylabel=r\"Percentage points (dev. from ss)\",\n", + " ylabel=r\"Percentage point (dev. from ss)\",\n", " T_plot=50,\n", " figsize=(18, 6),\n", "):\n", @@ -734,13 +820,17 @@ " for i in range(n_var):\n", " # plot all irfs\n", " for j, irf in enumerate(irfs_list):\n", - " ax[i].plot(irf[variables[i]][:50], label=labels[j])\n", + " ax[i].plot(irf[variables[i]][:50] * 1e2, label=labels[j])\n", " ax[i].set_title(variables[i])\n", - " ax[i].set_xlabel(r\"$t$\")\n", + " ax[i].set_xlabel(r\"Quarters\")\n", " if i == 0:\n", " ax[i].set_ylabel(ylabel)\n", - " ax[i].legend()\n", - " plt.show()" + " # ax[i].legend()\n", + " # ax[i].x\n", + "\n", + " plt.tight_layout()\n", + " plt.show()\n", + " # plt.savefig(\"irf.png\")" ] }, { @@ -777,7 +867,7 @@ "id": "a0130fc3", "metadata": {}, "source": [ - "# Simulating the model" + "### Simulating the model" ] }, { @@ -829,7 +919,6 @@ "rhos = {\"Z\": 0.8}\n", "impulses = {}\n", "\n", - "\n", "for i in inputs:\n", " own_shock = {i: sigmas[i] * rhos[i] ** np.arange(T)}\n", " impulses[i] = ks.solve_impulse_linear(\n", @@ -844,14 +933,6 @@ "data_simul = simulate(list(impulses.values()), outputs, T_sim)\n", "plot_timeseries(data_simul, (1, 3), figsize=(12, 4))" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "27732edd", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": {