From 518e9b3212bd97dce53882a4a9f6b2df4837d71a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C2=A8ezgioz=C2=A8?= <¨ezgiozsogut@gmail.com¨> Date: Fri, 4 Oct 2019 17:23:21 +0200 Subject: [PATCH 1/3] Carroll_2017_MPC paper model with no aggregate shock --- examples/Caroll_2017_no_aggregate.yaml | 71 ++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 examples/Caroll_2017_no_aggregate.yaml diff --git a/examples/Caroll_2017_no_aggregate.yaml b/examples/Caroll_2017_no_aggregate.yaml new file mode 100644 index 0000000..207290f --- /dev/null +++ b/examples/Caroll_2017_no_aggregate.yaml @@ -0,0 +1,71 @@ +symbols: + exogenous: [w,r,ʊ,ψ,ξ] + states: [m] + controls: [c] + parameters: [ρ, β, δ, D, Đ, wbar, rbar, ʊbar, μ_ue, τ, θ, μ_θ, σ_θ, l, val] + + +equations: + + transition: + - m = ((D+r)/Đ)*(m(-1)-c(-1))*(1/ψ) + ξ + + arbitrage: + - ( β*(c(1)/c)^(-ρ) )*(ψ(1)^(-ρ))*(D+r(1))/Đ - 1 | 0.0<=c<=m + + + +calibration: + ρ: 1 + β: 0.99 + δ: 0.025 + D: 1-δ + Đ: 1-0.00625 + wbar: 1 + rbar: 0.01 + δ + ʊbar: 0.07 + μ_ue: 0.15 + τ: 0 + l: 1/0.9 + + μ_θ: 1 + σ_θ: 0.010*4 + μ_ψ: 1 + σ_ψ: (0.010*4)/11 + + val: (D+r)/Đ + m: 1.2 + c: (m*(val-1) + 1 )/val + θ: 1 + + w : wbar + r : rbar + ʊ : ʊbar + ψ : 1 + ξ : 1 + +domain: + m: [0.1, 10.0] + +exogenous: + w,r,ʊ: !ConstantProcess + μ: [wbar, rbar, ʊbar] + + ψ: !LogNormal + μ: 1 + σ: 0.00004 + + ξ: !Mixture + index: !Bernouilli + π: ʊbar + distributions: + 0: !ConstantProcess + μ: [μ_ue] + 1: !LogNormal + μ: (1-τ)*l*μ_θ + σ: (((1-τ)*l)^2)*(σ_θ)^2 + ### transformed for μ_e=(1-τ)*l*θ) + +options: + grid: !Cartesian + orders: [1000] From d164b79b71b5bfeb943db4a5a307552e5f125e1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C2=A8ezgioz=C2=A8?= <¨ezgiozsogut@gmail.com¨> Date: Fri, 11 Oct 2019 18:01:42 +0200 Subject: [PATCH 2/3] Create heterogenous agent version of CS model from Carol 2017 --- examples/bfs_2017_beta_het.yaml | 75 +++++++++++++++++++ examples/bfs_2017_het.yaml | 105 ++++++++++++++++++++++++++ examples/bfs_het_example.py | 128 ++++++++++++++++++++++++++++++++ 3 files changed, 308 insertions(+) create mode 100644 examples/bfs_2017_beta_het.yaml create mode 100644 examples/bfs_2017_het.yaml create mode 100644 examples/bfs_het_example.py diff --git a/examples/bfs_2017_beta_het.yaml b/examples/bfs_2017_beta_het.yaml new file mode 100644 index 0000000..d203aaa --- /dev/null +++ b/examples/bfs_2017_beta_het.yaml @@ -0,0 +1,75 @@ +symbols: + exogenous: [w,r,ʊ,ψ,ξ] + states: [m] + controls: [c] + parameters: [ρ, β, δ, D, Đ, wbar, rbar, ʊbar, μ_ue, τ, θ, μ_θ, σ_θ, l, val] + + +equations: + + transition: + - m = ((D+r)/Đ)*(m(-1)-c(-1))*(1/ψ) + ξ + + arbitrage: + - ( β*(c(1)/c)^(-ρ) )*(ψ(1)^(-ρ))*(D+r(1))/Đ - 1 | 0.0<=c<=m + + + +calibration: + ρ: 1 + β: 0.99 + δ: 0.025 + D: 1-δ + Đ: 1-0.00625 + wbar: 1 + rbar: 0.01 + δ + ʊbar: 0.07 + μ_ue: 0.15 + τ: 0 + l: 1/0.9 + + μ_θ: 1 + σ_θ: 0.010*4 + μ_ψ: 1 + σ_ψ: (0.010*4)/11 + + val: (D+r)/Đ + m: 1.2 + c: (m*(val-1) + 1 )/val + θ: 1 + + w : wbar + r : rbar + ʊ : ʊbar + ψ : 1 + ξ : 1 + +domain: + m: [0.1, 10.0] + +exogenous: + w,r,ʊ: !ConstantProcess + μ: [wbar, rbar, ʊbar] + + ψ: !LogNormal + μ: 1 + σ: 0.00004 + + ξ: !Mixture + index: !Bernouilli + π: 1-ʊbar + distributions: + 0: !ConstantProcess + μ: [μ_ue] + 1: !UNormal + μ: (1-τ)*l + σ: (1-τ)*l*σ_θ + + # 1: !LogNormal + # μ: (1-τ)*l*μ_θ + # σ: (((1-τ)*l)^2)*(σ_θ)^2 + ### transformed for μ_e=(1-τ)*l*θ) + +options: + grid: !Cartesian + orders: [1000] diff --git a/examples/bfs_2017_het.yaml b/examples/bfs_2017_het.yaml new file mode 100644 index 0000000..bf5d8fe --- /dev/null +++ b/examples/bfs_2017_het.yaml @@ -0,0 +1,105 @@ +name: Consumption Savings + +symbols: + exogenous: [w,r,ʊ,ψ,ξ] + states: [m] + controls: [c] + parameters: [ρ, β, δ, D, Đ, wbar, rbar, ʊbar, μ_ue, τ, μ_θ, σ_θ, μ_ψ, σ_ψ, l] + +equations: + + transition: + - m = ((D+r)/Đ)*(m(-1)-c(-1))*(1/ψ) + ξ + + arbitrage: + - ( β*(c(1)/c)^(-ρ) )*(ψ(1)^(-ρ))*(D+r(1))/Đ - 1 | 0.0<=c<=m + +calibration: + ρ: 1 + β: 0.99 + δ: 0.025 + D: 1-δ + Đ: 1-0.00625 + l: 1/0.9 + + wbar: 1 + rbar: 0.01 + δ + ʊbar: 0.07 + μ_ue: 0.15 + τ: 0 + + μ_θ: 1 + σ_θ: 0.010*4 + μ_ψ: 1 + σ_ψ: (0.010*4)/11 + + m: 1.2 + c: (m*(((D+r)/Đ)-1) + 1 )/((D+r)/Đ) + θ: 1 + + w : wbar + r : rbar + ʊ : ʊbar + ψ : 1 + ξ : 1 + +domain: + m: [0.1, 10.0] + +exogenous: + w,r,ʊ: !ConstantProcess + μ: [wbar, rbar, ʊbar] + + ψ: !UNormal + μ: 1 + σ: 0.00004 + + ξ: !Mixture + index: !Bernouilli + π: 1-ʊbar + distributions: + 0: !ConstantProcess + μ: [μ_ue] + 1: !UNormal + μ: (1-τ)*l + σ: (1-τ)*l*σ_θ + + # 1: !LogNormal + # μ: (1-τ)*l*μ_θ + # σ: (((1-τ)*l)^2)*(σ_θ)^2 + ### transformed for μ_e=(1-τ)*l*θ) + +options: + grid: !Cartesian + orders: [1000] + + +--- + +name: Caroll + +symbols: + exogenous: [Z, Ψ, Ξ] + aggregate: [K] + parameters: [α, YK, L, l] + +calibration: + α: 0.36 + YK: 10.26 + L: 1 + l: 1/0.9 + Z: 1.0 + Ψ: 1.0 + Ξ: 1.0 + K: l*(YK)**(1/(α-1)) + + +exogenous: !ConstantProcess + μ: [1.0, 1.0, 1.0] + +projection: + r: α*YK + w: (1-α)*Z*(K/(l*L))**(α) + +equilibrium: + - K = m-c diff --git a/examples/bfs_het_example.py b/examples/bfs_het_example.py new file mode 100644 index 0000000..cdec303 --- /dev/null +++ b/examples/bfs_het_example.py @@ -0,0 +1,128 @@ +# -*- coding: utf-8 -*- +# %% +from matplotlib import pyplot as plt + +# %% +# Let's import the heterogeneous agents model +from dolark import HModel +aggmodel = HModel('examples/bfs_2017_het.yaml', i_options={'to':'iid'}) +aggmodel = HModel('examples/ayiagari.yaml') + +aggmodel # TODO: find a reasonable representation of this object + +# %% [markdown] +# First we can check whether the one-agent sub-part of it works + +# %% +from dolo import time_iteration +discretization_options = {"N": 2} +model = aggmodel.model +mc = model.exogenous.discretize(to='mc', options=[{},discretization_options]) +sol0 = time_iteration(model, details=True, dprocess=mc) + +# %% +# TEMP: we need to supply a projection function, which maps aggregate variables +# into the exogenous shocks received by idiosyncratic agents. +# This should be read from the YAML file. For now, we monkey-patch + +def projection(self, m: 'n_e', y: "n_y", p: "n_p"): + + from numpy import exp + z = m[0] + K = y[0] + A = [0] + alpha = p[1] + delta = p[2] + N = 1 + r = alpha*exp(z)*(N/K)**(1-alpha) - delta + w = (1-alpha)*exp(z)*(K/N)**(alpha) + return {'r': r, "w": w} + +import types +aggmodel.projection = types.MethodType(projection, aggmodel) + +# %% +# We can now solve for the aggregate equilibrium + +eq = aggmodel.find_steady_state() +eq + +# %% +# lot's look at the aggregate equilibrium +for i in range(eq.μ.shape[0]): + s = eq.dr.endo_grid.nodes() # grid for states (temporary) + plt.plot(s, eq.μ[i,:]*(eq.μ[i,:].sum()), label=f"y={eq.dr.exo_grid.node(i)[2]: .2f}") +plt.plot(s, eq.μ.sum(axis=0), label='total', color='black') +plt.grid() +plt.legend(loc='upper right') +plt.title("Wealth Distribution by Income") + +# %% +# alternative way to plot equilibrium + +import altair as alt +df = eq.as_df() +spec = alt.Chart(df).mark_line().encode( + x = 'a', + y = 'μ', + color = 'i_m:N' +) +spec + +# %% +# alternative way to plot equilibrium (with some interactivity) +# TODO: function to generate it automatically. + +import altair as alt +single = alt.selection_single(on='mouseover', nearest=True) +df = eq.as_df() +ch = alt.Chart(df) +spec = ch.properties(title='Distribution', height=100).mark_line().encode( + x = 'a', + y = 'μ', + color = alt.condition(single, 'i_m:N', alt.value('lightgray')) +).add_selection( + single +) + ch.mark_line(color='black').encode( + x = 'a', + y = 'sum(μ)' +) & ch.properties(title='Decision Rule', height=100).mark_line().encode( + x = 'a', + y = 'i', + color = alt.condition(single, 'i_m:N', alt.value('lightgray')) +).add_selection( + single +) + +# %% + +# Resulting object can be saved to a file. (try to open this file in jupyterlab) +open('distrib.json','tw').write(spec.to_json()) + +# %% + +# %% +import xarray + +# %% +# now we compute the perturbation +peq = aggmodel.perturb(eq) + + +# %% +# and we simulate given initial value of aggregate shock +sim = peq.response([0.1]) + +# %% +plt.subplot(121) +for t, (m,μ,x,y) in enumerate(sim): + plt.plot(μ.sum(axis=0), color='red', alpha=0.01) +plt.xlabel('a') +plt.ylabel('density') +plt.grid() +plt.subplot(122) +plt.plot( [e[3][0] for e in sim]) +plt.xlabel("t") +plt.ylabel("k") +plt.grid() +plt.tight_layout() From 4e34ae504a08d66dfabfdb60d806d6207a8a7c2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C2=A8ezgioz=C2=A8?= <¨ezgiozsogut@gmail.com¨> Date: Fri, 20 Dec 2019 17:57:59 +0100 Subject: [PATCH 3/3] corrected yaml file --- examples/bfs_2017.py | 39 +++++++++++++++++++++++++++++++++------ examples/bfs_2017.yaml | 31 ++++++++++++++++--------------- 2 files changed, 49 insertions(+), 21 deletions(-) diff --git a/examples/bfs_2017.py b/examples/bfs_2017.py index e36fd43..7dc6bba 100644 --- a/examples/bfs_2017.py +++ b/examples/bfs_2017.py @@ -5,8 +5,8 @@ # text_representation: # extension: .py # format_name: light -# format_version: '1.4' -# jupytext_version: 1.2.4 +# format_version: '1.5' +# jupytext_version: 1.3.0 # kernelspec: # display_name: Python 3 # language: python @@ -15,6 +15,14 @@ # WARNING: this is not working yet +# + jupyter={"source_hidden": true} +from dolo import groot +# - + +groot() + +pwd + # + from dolo import * from matplotlib import pyplot as plt @@ -22,7 +30,7 @@ import altair as alt from dolark import HModel -hmodel = HModel("bfs_2017.yaml") +hmodel = HModel("examples/bfs_2017.yaml") hmodel.features # - hmodel.agent @@ -43,9 +51,24 @@ # Agents' decision rule is not correct yet +# + + + +from dolo.numeric.decision_rule import CustomDR + +values = { + #'c': 'min(m, 1 + 0.05*(m-1))' + 'c': 'm*0.2' +} + + +cdr = CustomDR(values, hmodel.agent) + +# - + from dolo import time_iteration, improved_time_iteration -dr = time_iteration(hmodel.agent) -# # %time dr = improved_time_iteration(hmodel.model, dr0=dr, verbose=True, details=False) +dr = time_iteration(hmodel.agent, dr0 = cdr, maxit=100) +# %time dr = improved_time_iteration(hmodel.model, dr0=dr, verbose=True, details=False) from dolo import tabulate @@ -53,16 +76,19 @@ tab = tabulate(hmodel.agent, dr, 'm') plt.plot(tab['m'], tab['c']) +#plt.ylim(0,1) # ergodic distribution (premature) Π, μ = ergodic_distribution(hmodel.model, dr) df_μ = μ.to_dataframe('μ').reset_index() +plt.plot(df_μ['m'], df_μ['μ']) + ch = alt.Chart(tab) g1 = ch.mark_line(color='black',strokeDash=[1,1]).encode(x='m', y='m') + \ ch.mark_line().encode(x='m', y='c') -g2 = alt.Chart(df_μ).mark_line().encode(x='m:Q', y= 'mu:Q') +g2 = alt.Chart(df_μ).mark_line().encode(x='m:Q', y= 'μ:Q') g2 @@ -74,5 +100,6 @@ # here are the values projected from market equilibrium, given default level of capital m0, y0, p = hmodel.calibration['exogenous','aggregate', 'parameters'] hmodel.projection(m0, y0, p) # values for r, w, ω (not the same at all) +# This should give w, r, ʊ diff --git a/examples/bfs_2017.yaml b/examples/bfs_2017.yaml index 1fa1dba..92f7a13 100644 --- a/examples/bfs_2017.yaml +++ b/examples/bfs_2017.yaml @@ -22,7 +22,7 @@ calibration: Đ: 1-0.00625 l: 1/0.9 - wbar: 1 + wbar: 2.37 rbar: 0.01 + δ ʊbar: 0.07 μ_ue: 0.15 @@ -44,14 +44,14 @@ calibration: ξ : 1 domain: - m: [0.1, 10.0] + m: [0.1, 20.0] exogenous: w,r,ʊ: !ConstantProcess μ: [wbar, rbar, ʊbar] - ψ: !UNormal - μ: 1 + ψ: !LogNormal + μ: 1.0 σ: 0.00004 ξ: !Mixture @@ -60,14 +60,15 @@ exogenous: distributions: 0: !ConstantProcess μ: [μ_ue] - 1: !UNormal - μ: (1-τ)*l - σ: (1-τ)*l*σ_θ + +# 1: !UNormal +# μ: (1-τ)*l +# σ: (1-τ)*l*σ_θ - # 1: !LogNormal - # μ: (1-τ)*l*μ_θ - # σ: (((1-τ)*l)^2)*(σ_θ)^2 - ### transformed for μ_e=(1-τ)*l*θ) + 1: !LogNormal + μ: (1-τ)*l*μ_θ + σ: (((1-τ)*l)**2)*(σ_θ)**2 + ## transformed for μ_e=(1-τ)*l*θ) distribution: β: !Uniform @@ -76,7 +77,7 @@ distribution: options: grid: !Cartesian - orders: [100] + orders: [200] --- @@ -90,7 +91,7 @@ symbols: calibration: α: 0.36 - YK: 10.26 + YK: 1/10.26 L: 1 l: 1/0.9 Z: 1.0 @@ -105,10 +106,10 @@ exogenous: !ConstantProcess μ: [1.0, 1.0, 1.0] projection: - w: (1-α)*Z*(K/(l*L))**(α) + w: (1-α)*Z*(YK)**(α/(α-1)) r: α*YK ʊ: ʊbar equilibrium: - - K = m-c + - K = m-c \ No newline at end of file