diff --git a/Project.toml b/Project.toml index 786d4eb..5aa5c67 100644 --- a/Project.toml +++ b/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" diff --git a/examples/models/consumption_savings_iid.yaml b/examples/models/consumption_savings_iid.yaml index b162031..d902161 100644 --- a/examples/models/consumption_savings_iid.yaml +++ b/examples/models/consumption_savings_iid.yaml @@ -49,6 +49,6 @@ exogenous: options: discretization: endo: - n: [100] + n: [1000] exo: n: 7 diff --git a/examples/models/rbc_mc.yaml b/examples/models/rbc_mc.yaml index c87cf4a..4be845c 100644 --- a/examples/models/rbc_mc.yaml +++ b/examples/models/rbc_mc.yaml @@ -79,7 +79,6 @@ exogenous: options: - - grid: !Cartesian - # mu: 3 - orders: [20] + discretization: + endo: + n: [200] diff --git a/examples/models/zero_to_aiyagari_iid.yaml b/examples/models/zero_to_aiyagari_iid.yaml new file mode 100644 index 0000000..759b181 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl new file mode 100644 index 0000000..396caeb --- /dev/null +++ b/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 + + + diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index 42be21a..347ae6c 100644 --- a/src/algos/ergodic.jl +++ b/src/algos/ergodic.jl @@ -6,15 +6,6 @@ # return AxisArray(μ; d...) # end - -function ergodic_distribution(model, sol) where T where Q - G = distG(model, sol) - μ = ergodic_distribution(G) - return μ - # return label_density(μ, model.symbols, sol.dr.grid_exo, sol.dr.grid_endo) -end - - # represents the way that distributions are discretized struct distG{ID, n_m, n_s, n_x, Gx, Ge} @@ -48,6 +39,23 @@ function distG(model, sol) end +""" +Computes the ergodic distribution of the states (endo & exo). +# Arguments +* `model::NumericModel`: Model object that describes the current model environment. +* `sol`: solution of the model obtained after a time iteration algorithm +# Optional Argument +* `smooth::boolean`: Indicates whether, when we discretize the transition results on the grid, we want to smooth the initial linear ponderation of nodes or not. +# Returns +* `μ1::distG`: the ergodic distribution +""" +function ergodic_distribution(model, sol; smooth=true) where T where Q + G = distG(model, sol) + μ = ergodic_distribution(G; smooth=smooth) + return μ + # return label_density(μ, model.symbols, sol.dr.grid_exo, sol.dr.grid_endo) +end + function ergodic_distribution(G::distG; x0=G.x0, kwargs...) P = transition_matrix(G.model, G.dprocess, x0, G.grid; kwargs...) M = [P'-I; ones(1, size(P,1))] @@ -58,22 +66,41 @@ function ergodic_distribution(G::distG; x0=G.x0, kwargs...) end - -function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; diff=false) where n_x - +""" +Uses the current distribution of the states (endogeneous & exogeneous) and the controls to compute the next distribution of the states on the grid. Exogeneous parameters can be included +as well as the possibility to differentiate this function with respect to the states distribution, the vector of controls and the exogeneous parameters. +# Arguments +* `μ0::AbstractVector{Float64}`: current distribution of the states (endo & exo) +* `x0::MSM{Point{n_x}}`: vector of controls +# Optional Arguments +* `exo`: nothing or (z0, z1) +* `diff::boolean`: Indicates whether we want to evaluate the derivative of G wrt μ, x and possibly z1 and z2 +* `smooth::boolean`: Indicates whether, when we discretize the transition results on the grid, we want to smooth the initial linear ponderation of nodes or not. +# Returns +* `μ1::distG`: It is the new distribution of the states (endogeneous & exogeneous) +# Optionally returns also +* `∂G_∂μ::LinearMaps.FunctionMap{Float64}`: It is the jacobian of G wrt μ. +* `∂G_∂x::LinearMaps.FunctionMap{Float64}`: It is the jacobian of G wrt x. +* `∂G_∂z1::LinearMaps.FunctionMap{Float64}`: It is the jacobian of G wrt z1. +* `∂G_∂z2::LinearMaps.FunctionMap{Float64}`: It is the jacobian of G wrt z2. +""" +function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false, smooth=true) where n_x if !diff - P = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=nothing, diff=false) + P = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=false) μ1 = P'*μ0 return μ1 end - P, P_x = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=nothing, diff=true) + P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true, smooth=smooth) - μ1 = P'μ0 + + μ1 = P'μ0 #new distribution of the states M = length(μ0) N = length(x0.data)*length(x0.data[1]) + N_z1 = length(exo[1]) + N_z2 = length(exo[2]) function fun_x(dx::AbstractVector{Float64}) d_x = MSM(copy(reinterpret(SVector{n_x, Float64}, dx)), x0.sizes) @@ -86,50 +113,37 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; diff=fals return d_μ end + function fun_z1(dz1) + P_dz1 = [(P_z1[i,j]'*dz1) for i=1:size(P,1), j=1:size(P,2)] + d_μ = P_dz1'*μ0 + return d_μ + end + + function fun_z2(dz2) + P_dz2 = [(P_z2[i,j]'*dz2) for i=1:size(P,1), j=1:size(P,2)] + d_μ = P_dz2'*μ0 + return d_μ + end + + ∂G_∂μ = LinearMap( μ -> P'*μ, M, M) ∂G_∂x = LinearMap(fun_x, M, N) - return μ1, ∂G_∂μ, ∂G_∂x + ∂G_∂z1 = LinearMap(fun_z1, M, N_z1) + ∂G_∂z2 = LinearMap(fun_z2, M, N_z2) + + return μ1, ∂G_∂μ, ∂G_∂x, ∂G_∂z1, ∂G_∂z2 end -function (G::distG)(μ0::AbstractVector{Float64}, x0::AbstractVector{Float64}; diff=false) +function (G::distG)(μ0::AbstractVector{Float64}, x0::AbstractVector{Float64}; exo=nothing, diff=false, smooth=true) n_x = length(G.x0.data[1]) x = MSM(copy(reinterpret(SVector{n_x, Float64}, x0)), G.x0.sizes) - return G(μ0, x; diff=diff) + return G(μ0, x; exo=exo, diff=diff, smooth=smooth) end -""" -Updates A. -# Arguments -* `A`: the transition matrix that will be updated. -* `x` : vector of controls. -* `w::Float64` : vector of weights. -* `exo_grid` : exogenous grid that will be used to determine the type of rescaling to do. -* `a` : SVector containing the minimum values on the UCGrids -* `b` : SVector containing the maximum values on the UCGrids - -# Optional -* `M` : future node considered when the exogenous grid is a UCGrid to rescale x - -# Modifies -* `A` : the updated transition matrix -""" -function trembling_hand_rescaled!(A, x, w::Float64, exo_grid, a, b; M=0) - if typeof(exo_grid) == Dolo.UnstructuredGrid{ndims(exo_grid)} - x = [(x[n]-a)./(b-a) for n=1:length(x)] - trembling_hand!(A, x, w) - elseif typeof(exo_grid) == Dolo.UCGrid{ndims(exo_grid)} - V = [(SVector(M..., el...)-a)./(b.-a) for el in x] - trembling_hand!(A, V, w) - else - - x = [(x[n]-a)./(b-a) for n=1:length(x)] - trembling_hand!(A, x, w) - end -end """ Calculates the new transition matrix for a given model, a given discretized exogenous process, given control values (x0) and given grids (exogenous and endogenous). @@ -137,14 +151,16 @@ Calculates the new transition matrix for a given model, a given discretized exog * `model::NumericModel`: Model object that describes the current model environment. * `dprocess::`: Discretized exogenous process. * `x0::Dolo.MSM{SVector{2, Float64}}`: Initial control values. -* `exo_grid`: Exogenous grid that can be of type either UnstructuredGrid or UCGrid or EmptyGrid (in the three following functions). -* `endo_grid::UCGrid`: Endogenous grid. +* `grid`: grid +# Optional arguments * `exo`: nothing or (z0, z1) +* `diff`: true or false. Indicates whether we want to evaluate the derivative of G wrt μ, x and possibly z1 and z2 +* `smooth`: true or false. Indicates whether, when we discretize the transition results on the grid, we want to smooth the initial linear ponderation of nodes or not. # Returns * `Π0::`: New transition matrix. """ -function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing, diff=false) where n_x +function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing, diff=false, smooth = true) where n_x # the restriction here is that dp is a descrete_process @@ -154,11 +170,17 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing endo_grid = grid.endo N_m = max(1, n_nodes(exo_grid)) - N_s = n_nodes(endo_grid) + N_s = n_nodes(endo_grid) N = N_m*N_s - Π = zeros(N_s, N_m, endo_grid.n..., N_m) - if diff - dΠ = zeros(SVector{n_x, Float64}, N_s, N_m, endo_grid.n..., N_m) + Π = zeros(N_s, N_m, endo_grid.n..., N_m) #initialization of the transition matrix + if diff #initialization of its derivatives + dΠ_x = zeros(SVector{n_x, Float64}, N_s, N_m, endo_grid.n..., N_m) + if !(exo === nothing) + n_z1 = length(exo[1]) + n_z2 = length(exo[2]) + dΠ_z1 = zeros(SVector{n_z1, Float64}, N_s, N_m, endo_grid.n..., N_m) + dΠ_z2 = zeros(SVector{n_z2, Float64}, N_s, N_m, endo_grid.n..., N_m) + end end s = nodes(endo_grid) a = SVector(endo_grid.min...) @@ -167,7 +189,7 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing x = x0.views[i_m] m = node(exo_grid, i_m) if !(exo === nothing) - m = Dolo.repsvec(exo[1], m) # z0 + m = Dolo.repsvec(exo[1], m) # completes z0 by elements of m until reaching a length equal to max(length(m),length(z0)) end for i_M in 1:n_inodes(dp, i_m) @@ -176,13 +198,14 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing else i_MM = i_M end - M = inode(Point, dp, i_m, i_M) + M = inode(Point, dp, i_m, i_M) # future node if !(exo === nothing) M = Dolo.repsvec(exo[2], M) # z1 end w = iweight(dp, i_m, i_M) if diff - S, S_x = transition(model, Val{(0,3)}, m, s, x, M, parms) + S, S_z1, S_x, S_z2 = transition(model, Val{(0,1,3,4)}, m, s, x, M, parms) # computes the next vector of states and its derivatives wrt z1, x and z2. + # The points obtained are located between points of the endogenous grid and a discretization needs therefore to be done to adjust the transition matrix. else S = transition(model, m, s, x, M, parms) end @@ -190,11 +213,22 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing ind_s = tuple((Colon() for k in 1:(ndims(Π)-3))...) Π_view = view(Π,:,i_m,ind_s..., i_MM) if !diff - trembling_hand!(Π_view, S, w) + trembling_hand!(Π_view, S, w; smooth=smooth) else + S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] - dΠ_view = view(dΠ,:,i_m,ind_s...,i_MM) - trembling_foot!(Π_view, dΠ_view, S, S_x, w) + S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] for n=1:length(S)] + S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] for n=1:length(S)] + + S_z1 = [M[:,1:length(exo[1])] for M in S_z1] + S_z2 = [M[:,1:length(exo[2])] for M in S_z2] + + + dΠ_view_x = view(dΠ_x,:,i_m,ind_s...,i_MM) + dΠ_view_z1 = view(dΠ_z1,:,i_m,ind_s...,i_MM) + dΠ_view_z2 = view(dΠ_z2,:,i_m,ind_s...,i_MM) + + trembling_foot!(Π_view, dΠ_view_x, dΠ_view_z1, dΠ_view_z2, S, S_x, S_z1, S_z2, w; smooth = smooth) end end end @@ -203,36 +237,60 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing if !diff return Π0 else - dΠ0 = reshape(dΠ,N,N) - return Π0, dΠ0 + dΠ_0_x = reshape(dΠ_x,N,N) + dΠ_0_z1 = reshape(dΠ_z1,N,N) + dΠ_0_z2 = reshape(dΠ_z2,N,N) + return Π0, dΠ_0_x, dΠ_0_z1, dΠ_0_z2 end end -function transition_matrix(model, sol; diff=false, exo=nothing) +function transition_matrix(model, sol; diff=false, exo=nothing, smooth=true) x0 = Dolo.MSM([sol.dr(i, sol.dr.grid_endo.nodes) for i=1:max(1,Dolo.n_nodes(sol.dr.grid_exo))]) grid = ProductGrid(sol.dr.grid_exo, sol.dr.grid_endo) - Dolo.transition_matrix(model, sol.dprocess, x0, grid; diff=diff, exo=nothing); + Dolo.transition_matrix(model, sol.dprocess, x0, grid; diff=diff, exo=nothing, smooth=smooth); end -function transition_matrix(G::distG; dp=G.dprocess, x0=G.x0, grid=G.grid, exo=nothing, diff=false) +function transition_matrix(G::distG; dp=G.dprocess, x0=G.x0, grid=G.grid, exo=nothing, diff=false, smooth=true) - return transition_matrix(G.model, dp, x0, grid; exo=exo, diff=diff) + return transition_matrix(G.model, dp, x0, grid; exo=exo, diff=diff,smooth=smooth) end """ -Updates A. +Evaluate a polynomial that allows to smooth the repartition done by trembling_foot. +# Arguments +* `x`: a float number. +# Returns +* the value of 2x^3 - 3x^2 +1 +""" +function smoothing(x::SVector{2}) + return 2. .* x .^3 .- 3. .* x.^2 .+1. +end + +""" +Evaluate the derivative of a polynomial that allows to smooth the repartition done by trembling_foot. +# Arguments +* `x`: a float number. +# Returns +* the value of the derivative of x |—> 2x^3 - 3x^2 +1 +""" +function ∂smoothing(x::SVector{2}) + return 6. .* x .^2 .- 6. .* x +end +""" +Updates A. # Arguments * `A`: the transition matrix that will be updated. * `x::Vector{Point{d}}` : vector of controls. * `w::Float64` : vector of weights. - +# Optional Argument +* `smooth::boolean` : indicates whether the discretization is made with a smooth ponderation or a linear one # Modifies * `A` : the updated transition matrix """ -function trembling_hand!(A, x::Vector{Point{d}}, w::Float64) where d +function trembling_hand!(A, x::Vector{Point{d}}, w::Float64; smooth=true) where d @assert ndims(A) == d+1 shape_A = size(A) @@ -250,8 +308,12 @@ function trembling_hand!(A, x::Vector{Point{d}}, w::Float64) where d λn = (xn./δ.-qn) # ∈[0,1[ by construction qn_ = round.(Int,qn) .+ 1 - λn_weight_vector = tuple( (SVector((1-λn[i]),λn[i]) for i in 1:d)... ) - + if smooth==false + λn_weight_vector = tuple( (SVector((1-λn[i]),λn[i]) for i in 1:d)... ) + else + λn_weight_vector = tuple( (smoothing(SVector((λn[i]),1-λn[i])) for i in 1:d)... ) + end + indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) # Filling transition matrix @@ -263,30 +325,58 @@ function trembling_hand!(A, x::Vector{Point{d}}, w::Float64) where d end -function trembling_foot!(Π, dΠ, S::Vector{Point{d}}, S_x::Vector{SMatrix{d,n_x,Float64,_}}, w::Float64) where d where n_x where _ +""" +Updates Π, the transition matrix, and its derivatives with respect to x, z1 and z2. This is done using the previous Π, dΠ_x, dΠ_z1 and dΠ_z2 and the derivatives +of the transition function wrt x, z1 and z2. +# Arguments +* `Π`: the transition matrix that will be updated. +* `dΠ_x` : its derivative wrt x +* `dΠ_z1` : its derivative wrt z1 +* `dΠ_z2` : its derivative wrt z2 +* `S::Vector{Point{d}}` : Vector of SVectors of states (endo & exo) with d the dimension of the state space and the size of the vector corresponding to the size of the grid +* `S_x` : Derivative of S wrt x +* `S_z1` : Derivative of S wrt z1 +* `S_z2` : Derivative of S wrt z2 +* `w::Float64` : weight +# Optional Argument +* `smooth::boolean` : indicates whether the discretization is made with a smooth ponderation or a linear one +# Modifies +* `Π`: the transition matrix that will be updated. +* `dΠ_x`: its derivative wrt x +* `dΠ_z1`: its derivative wrt z1 +* `dΠ_z2`: its derivative wrt z2 +""" + +function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Vector{SMatrix{d,n_x,Float64,_}}, S_z1, S_z2, w::Float64; smooth=true) where d where n_x where _ @assert ndims(Π) == d+1 shape_Π = size(Π) grid_dimension = d δ = SVector{d,Float64}(1.0./(shape_Π[1+i]-1) for i in 1:d ) N = shape_Π[1] + for n in 1:N Sn = S[n] Sn_x = S_x[n] - - inbound = (0. .<= Sn) .& ( Sn .<= 1.0) - + Sn_z1 = S_z1[n] + Sn_z2 = S_z2[n] + + inbound = (0. .<= Sn) .& ( Sn .<= 1.0) # detects the points that have fallen outside of the grid + Sn = min.(max.(Sn, 0.0),1.0) qn = div.(Sn, δ) qn = max.(0, qn) - # inbound = inbound .& (qn .<= shape_Π[2:d+1].-2) qn = min.(qn, shape_Π[2:d+1].-2) λn = (Sn./δ.-qn) # ∈[0,1[ by construction qn_ = round.(Int,qn) .+ 1 - λn_weight_vector_Π = tuple( (SVector((1-λn[i]),λn[i]) for i in 1:d)... ) + if smooth==false + λn_weight_vector_Π = tuple( (SVector((1-λn[i]),λn[i]) for i in 1:d)... ) # linear ponderation + else + λn_weight_vector_Π = tuple( (smoothing(SVector((λn[i]),1-λn[i])) for i in 1:d)... ) # smoothed ponderation —> allows to have continuous derivatives + end # # # Filling transition matrix rhs_Π = outer(λn_weight_vector_Π...) @@ -294,14 +384,28 @@ function trembling_foot!(Π, dΠ, S::Vector{Point{d}}, S_x::Vector{SMatrix{d,n_x indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) Π[indexes_to_be_modified...] .+= w.*rhs_Π - + + # # # Computing the associated derivatives for k=1:d - λ_vec = tuple( (i==k ? SVector( -1. /δ[k] * inbound[k], 1. / δ[k] * inbound[k]) : (SVector((1-λn[i]),λn[i])) for i in 1:d)... ) + + if smooth==false + λ_vec = tuple( (i==k ? SVector( -1. /δ[k] * inbound[k], 1. / δ[k] * inbound[k]) : (SVector((1-λn[i]),λn[i])) for i in 1:d)... ) + else + λ_vec = tuple( (i==k ? ∂smoothing(SVector((λn[i]),1-λn[i])) .* SVector{2,Float64}(1.,-1.) ./δ[k] .* inbound[k] : (smoothing(SVector((λn[i]),1-λn[i]))) for i in 1:d)... ) + end + A = outer(λ_vec...) - rhs_dΠ = outer2(A, Sn_x[k,:]) - dΠ[indexes_to_be_modified...] .+= w*rhs_dΠ + + rhs_dΠ_x = outer2(A, Sn_x[k,:]) + dΠ_x[indexes_to_be_modified...] .+= w*rhs_dΠ_x + + rhs_dΠ_z1 = outer2(A, Sn_z1[k,:]) + dΠ_z1[indexes_to_be_modified...] .+= w*rhs_dΠ_z1 + + rhs_dΠ_z2 = outer2(A, Sn_z2[k,:]) + dΠ_z2[indexes_to_be_modified...] .+= w*rhs_dΠ_z2 + end - end end @@ -310,10 +414,8 @@ end """ Computes the outer product. - # Argument * `λn_weight_vector::Vararg{Point{2}}`: tuple of Point{2} to be multiplied by outer product - # Returns * the outer product """ @@ -321,15 +423,12 @@ function outer(λn_weight_vector::Vararg{SVector{2}}) return [prod(e) for e in Iterators.product(λn_weight_vector...)] end -# TODO: define and document """ Computes an outer product between a vector and a matrix and returns a vector of matrices. - # Arguments * `A` : SVector of any size * `x`: matrix of any size - # Returns * a sized vector obtained by making the outer product of A and x and dividing the matrix obtained in a vector of matrices for which the value of the A[i] term changes """ diff --git a/test/test_G_derivative.jl b/test/test_G_derivative.jl index 3b95d6e..953adcb 100644 --- a/test/test_G_derivative.jl +++ b/test/test_G_derivative.jl @@ -1,36 +1,48 @@ using FiniteDiff using Dolo -@testset "Test the G derivative w.r.t. x and μ" begin +@testset "Test the G derivative w.r.t. x, μ, z1 and z2" begin model_rbc_mc = yaml_import("examples/models/rbc_mc.yaml") model_rbc = yaml_import("examples/models/rbc.yaml") model_rbc_iid = yaml_import("examples/models/rbc_iid.yaml") + for model in [model_rbc_mc, model_rbc, model_rbc_iid] 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 - x0_not_flat = G.x0 - x0_flat = cat(G.x0.data...; dims = 1) - μ0 = G.μ0 - for x0 in [x0_not_flat, x0_flat] + μ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) + + 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) + + - μ1, ∂G_∂μ, ∂G_∂x = G(μ0, x0, diff = true) + @assert maximum(abs, Jμ_num - Jμ_exact) < 1e-8 - Jμ_exact = convert(Matrix, ∂G_∂μ) - Jμ_num = FiniteDiff.finite_difference_jacobian(mu -> G(mu, x0), μ0) + @assert maximum(abs, Jx_num - Jx_exact) < 1e-8 - Jx_exact = convert(Matrix, ∂G_∂x) - Jx_num = FiniteDiff.finite_difference_jacobian(x -> G(μ0, x), x0) + @assert maximum(abs, Jz1_num - Jz1_exact) < 1e-8 - @assert maximum(abs, Jμ_num - Jμ_exact) < 1e-8 + @assert maximum(abs, Jz2_num - Jz2_exact) < 1e-5 - @assert maximum(abs, Jx_num - Jx_exact) < 1e-8 - end end