Gqueran/derivative wrt z

Project.toml

@@ -22,6 +22,7 @@ LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuantEcon = "fcd29c91-0bd7-5a09-975d-7ac3f643a60c"

examples/models/consumption_savings_iid.yaml

@@ -49,6 +49,6 @@ exogenous:
 options:
     discretization:
         endo:
-            n: [100]
+            n: [1000]
         exo:
             n: 7

examples/models/rbc_mc.yaml

@@ -79,7 +79,6 @@ exogenous:
 
 options:
 
-
-    grid: !Cartesian
-        # mu: 3
-        orders: [20]
+    discretization:
+        endo:
+            n: [200]

examples/models/zero_to_aiyagari_iid.yaml

@@ -0,0 +1,52 @@
+name: Consumption Savings
+
+symbols:
+    states: [m]
+    exogenous: [e]
+    parameters: [beta, B, a_max, r, w]
+    controls: [c]
+
+
+# definitions: |
+#     i[t] = m[t] - c[t]
+
+equations:
+
+    arbitrage: |
+        beta*(1+r)*(c[t+1]/c[t])^(-4.0)-1 ⟂ 0 <= c[t] <= m[t]
+
+    transition: |
+        m[t] = w*exp(e[t]) + (m[t-1]-c[t-1])*(1+r)
+
+calibration:
+
+    beta: 0.99
+    B: -1
+    a_max: 10
+    m: 1.0
+    c: 2
+    a: 1
+    i: a
+    K: 40.
+    alpha: 0.36
+    A: 1
+    N: 1
+    delta: 0.045
+    r: alpha*(N/K)^(1-alpha) - delta
+    w: (1-alpha)*(K/N)^(alpha)
+    e: 0
+
+
+exogenous:  
+    e: !UNormal
+        σ: 0.1
+
+
+domain:
+    m: [0.1, a_max]
+
+
+options:
+    discretization:
+        endo:
+            n: 100

experiments/dev_ergodic.jl

@@ -0,0 +1,274 @@
+using Dolo 
+using BenchmarkTools
+using StaticArrays
+using FiniteDiff
+
+
+model = yaml_import("examples/models/rbc.yaml")
+sol = Dolo.improved_time_iteration(model)
+
+
+G = Dolo.distG(model, sol)
+z10 = SVector(model.calibration[:exogenous]...) 
+z20 = z10 
+x0 = G.x0
+x0_flat = cat(G.x0.data...; dims=1)
+μ0 = G.μ0
+
+
+
+μ1, ∂G_∂μ, ∂G_∂x, ∂G_∂z1, ∂G_∂z2 = G(μ0, x0; exo = [z10,z20], diff = true)
+
+Jμ_exact = convert(Matrix, ∂G_∂μ)
+Jμ_num = FiniteDiff.finite_difference_jacobian(mu -> G(mu, x0, exo = [z10,z20]), μ0)
+
+Jx_exact = convert(Matrix, ∂G_∂x)
+Jx_num = FiniteDiff.finite_difference_jacobian(x -> G(μ0, x; exo = [z10,z20]), x0_flat)
+
+Jz1_exact = convert(Matrix, ∂G_∂z1)
+Jz1_num = FiniteDiff.finite_difference_jacobian(z1 -> G(μ0, x0; exo = [z1,z20]), z10)
+
+Jz2_exact = convert(Matrix, ∂G_∂z2)
+Jz2_num = FiniteDiff.finite_difference_jacobian(z2 -> G(μ0, x0; exo = [z10,z2]), z20)
+
+print( maximum(abs, Jμ_num - Jμ_exact) )
+
+print(maximum(abs, Jx_num - Jx_exact))
+
+print(maximum(abs, Jz1_num - Jz1_exact) )
+
+print(maximum(abs, Jz2_num - Jz2_exact))
+
+# Plotting the differences between exact values calculated in ergordic.jl and numerical values calculated by FiniteDiff
+using Plots
+
+pl1 = spy(abs.(Jx_num).>1e-10, title="Numerical")
+pl2 = spy(abs.(Jx_exact).>1e-10, title="Exact")
+pl3 = spy(abs.(Jx_exact - Jx_num).>1e-10, title="Diff")
+plot(pl1,pl2,pl3)
+
+p1 = plot([Jz1_num,Jz1_exact], label=["Jz1_num" "Jz1_exact"])
+p2 = plot([Jz1_num[1:50],Jz1_exact[1:50]], label=["Jz1_num_zoom" "Jz1_exact_zoom"])
+p3 = plot((Jz1_num-Jz1_exact), label = "diff")
+p4 = scatter((Jz1_num-Jz1_exact)[abs.(Jz1_num-Jz1_exact).>1e-8], label = "diff>1E-8", linestyle = :dot)
+plot(p1,p2,p3,p4, layout = (4,1))
+
+
+p1 = plot([Jz2_num,Jz2_exact], label=["Jz2_num" "Jz2_exact"])
+p2 = plot([Jz2_num[31:60],Jz2_exact[31:60]], label=["num_zoom" "exact_zoom"])
+p3 = plot((Jz2_num-Jz2_exact), label = "diff")
+p4 = scatter((Jz2_num-Jz2_exact)[abs.(Jz2_num-Jz2_exact).>1e-8], label = "diff>1E-8", linestyle = :dot)
+plot(p1,p2,p3,p4, layout = (4,1))
+
+
+# Comparing the ergodic distributions when there is smoothing and when there is not
+
+using PlotlyJS
+using DataFrames
+
+## rbc markov chain model
+
+model = yaml_import("examples/models/rbc_mc.yaml")
+sol = Dolo.improved_time_iteration(model)
+
+nodes = sol.dr.grid.endo.nodes
+
+k = [nodes[i][1] for i in 1:length(nodes)]
+
+μ_no_smoothing = Dolo.ergodic_distribution(model, sol; smooth=false)
+μ_smoothing = Dolo.ergodic_distribution(model, sol;  smooth=true)
+
+n = Int(length(μ_smoothing)/2)
+df1 = DataFrame(k=vcat(k,k), μ=vcat(2*μ_smoothing[1:n], 2*μ_no_smoothing[1:n]), smooth=vcat(["smoothing" for i in 1:n], ["no smoothing" for i in 1:n]))
+df2 = DataFrame(k=vcat(k,k), μ=vcat(2*μ_smoothing[n+1:2*n], 2*μ_no_smoothing[n+1:2*n]), smooth=vcat(["smoothing" for i in 1:n], ["no smoothing" for i in 1:n]))
+
+
+fig = PlotlyJS.make_subplots(
+    rows=2, cols=1,
+    column_widths=[1.],
+    row_heights=[0.5, 0.5]
+)
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df1[1:n,:],  x=:k, y=:μ, name  = "smoothing"),
+    row=1, col = 1
+    )
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df1[n+1:2*n,:],  x=:k, y=:μ, name = "no smoothing"),
+    row=1, col = 1
+    )
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df2[1:n,:],  x=:k, y=:μ, name = "smoothing"),
+    row=2, col = 1
+    )
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df2[n+1:2*n,:],  x=:k, y=:μ, name = "no smoothing"),
+    row=2, col = 1
+    )
+
+relayout!(
+    fig,
+    margin=attr(r=10, t=25, b=40, l=60),
+    annotations=[
+        attr(
+            text="k",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=-0.05),
+        attr(
+            text="k",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=0.55),
+        attr(
+            text="Exogeneous shock 2",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=0.45),
+        attr(
+            text="Exogeneous shock 1",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=1.05),
+        attr(
+            text="μ",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=-0.05,
+            y=0.8),
+        attr(
+            text="μ",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=-0.05,
+            y=0.2)
+    ]
+)
+
+fig
+
+
+
+
+
+
+## Aiyagari model
+model = yaml_import("examples/models/consumption_savings_mc.yaml")
+sol = Dolo.time_iteration(model)
+
+sim = Dolo.tabulate(model, sol.dr, :w)
+
+plot(sim[:w],sim[:c], ylabel="c", xlabel = "w",title = "Consumption")
+nodes = sol.dr.grid.endo.nodes
+
+w = [nodes[i][1] for i in 1:length(nodes)]
+
+μ_no_smoothing = Dolo.ergodic_distribution(model, sol; smooth=false)
+μ_smoothing = Dolo.ergodic_distribution(model, sol;  smooth=true)
+
+
+n = Int(length(μ_smoothing)/2)
+df1 = DataFrame(w=vcat(w,w), μ=vcat(2*μ_smoothing[1:n], 2*μ_no_smoothing[1:n]), smooth=vcat(["smoothing" for i in 1:n], ["no smoothing" for i in 1:n]))
+df2 = DataFrame(w=vcat(w,w), μ=vcat(2*μ_smoothing[n+1:2*n], 2*μ_no_smoothing[n+1:2*n]), smooth=vcat(["smoothing" for i in 1:n], ["no smoothing" for i in 1:n]))
+
+
+fig = PlotlyJS.make_subplots(
+    rows=2, cols=1,
+    column_widths=[1.],
+    row_heights=[0.5, 0.5]
+)
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df1[1:n,:],  x=:w, y=:μ, name  = "smoothing"),
+    row=1, col = 1
+    )
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df1[n+1:2*n,:],  x=:w, y=:μ, name = "no smoothing"),
+    row=1, col = 1
+    )
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df2[1:n,:],  x=:w, y=:μ, name = "smoothing"),
+    row=2, col = 1
+    )
+
+PlotlyJS.add_trace!(
+    fig,
+    PlotlyJS.scatter(df2[n+1:2*n,:],  x=:w, y=:μ, name = "no smoothing"),
+    row=2, col = 1
+    )
+
+relayout!(
+    fig,
+    margin=attr(r=10, t=25, b=40, l=60),
+    annotations=[
+        attr(
+            text="w",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=-0.05),
+        attr(
+            text="w",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=0.55),
+        attr(
+            text="Exogeneous shock 2",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=0.45),
+        attr(
+            text="Exogeneous shock 1",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=0.5,
+            y=1.05),
+        attr(
+            text="μ",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=-0.05,
+            y=0.8),
+        attr(
+            text="μ",
+            showarrow=false,
+            xref="paper",
+            yref="paper",
+            x=-0.05,
+            y=0.2)
+    ]
+)
+
+fig
+
+
+