From b87b19b61740c1c2b71659491dc8c96fbf669cfe Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 09:51:42 +0100 Subject: [PATCH 01/22] trying first to modify the function transition matrix to derive G wrt z1, z2 --- experiments/dev_G_derivative_wrt_z.jl | 125 ++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 experiments/dev_G_derivative_wrt_z.jl diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl new file mode 100644 index 0000000..fc3f6d2 --- /dev/null +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -0,0 +1,125 @@ +function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) where n_x + + + if !diff + 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=exo, diff=true) + + μ1 = P'μ0 + + M = length(μ0) + N = length(x0.data)*length(x0.data[1]) + + function fun_x(dx::AbstractVector{Float64}) + d_x = MSM(copy(reinterpret(SVector{n_x, Float64}, dx)), x0.sizes) + N = size(P,2) + d_μ = [ sum([μ0[k]*(P_x[k,i]'*d_x.data[k]) for k=1:N]) for i=1:N ] + # alternate calculation + P_dx = [(P_x[i,j]'*d_x.data[i]) for i=1:size(P,1), j=1:size(P,2)] + d_μμ = P_dx'*μ0 + @assert (maximum(abs, d_μ-d_μμ)<1e-10) + return d_μ + end + + ∂G_∂μ = LinearMap( μ -> P'*μ, M, M) + ∂G_∂x = LinearMap(fun_x, M, N) + return μ1, ∂G_∂μ, ∂G_∂x + +end + +function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing, diff=false) where n_x + + # the restriction here is that dp is a descrete_process + + parms = SVector(model.calibration[:parameters]...) + + exo_grid = grid.exo + endo_grid = grid.endo + + N_m = max(1, n_nodes(exo_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Π_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...) + b = SVector(endo_grid.max...) + for i_m in 1:max(1, n_nodes(exo_grid)) + x = x0.views[i_m] + m = node(exo_grid, i_m) + if !(exo === nothing) + m = Dolo.repsvec(exo[1], m) # z0 + dm_dz = Dolo.repsvec((@SVector ones(n_z1)),m*0) + end + for i_M in 1:n_inodes(dp, i_m) + + if isa(grid.exo, EmptyGrid) + i_MM = 1 # used to index transition matrix + else + i_MM = i_M + end + M = inode(Point, dp, i_m, i_M) + if !(exo === nothing) + M = Dolo.repsvec(exo[2], M) # z1 + dM_dz = Dolo.repsvec((@SVector ones(n_2)),M*0) + end + w = iweight(dp, i_m, i_M) + if diff + S, S_z1, S_x, S_z2 = transition(model, Val{(0,1,3,4)}, m, s, x, M, parms) + else + S = transition(model, m, s, x, M, parms) + end + S = [(S[n]-a)./(b-a) for n=1:length(S)] + 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) + else + S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] + S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] for n=1:length(S)] * dm_dz + S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] for n=1:length(S)] * dM_dz + + 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) + end + end + end + + Π0 = (reshape(Π,N,N)) + if !diff + return Π0 + else + 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) + 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); +end + +function transition_matrix(G::distG; dp=G.dprocess, x0=G.x0, grid=G.grid, exo=nothing, diff=false) + + return transition_matrix(G.model, dp, x0, grid; exo=exo, diff=diff) +end From a2d0ead911d82e5c1b212c3bdb188302d63c45e1 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 10:11:46 +0100 Subject: [PATCH 02/22] =?UTF-8?q?modifying=20trembling=20foot=20to=20updat?= =?UTF-8?q?e=20d=CE=A0=5Fz1=20and=20d=CE=A0=5Fz2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- experiments/dev_G_derivative_wrt_z.jl | 53 +++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index fc3f6d2..30f59e1 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -123,3 +123,56 @@ function transition_matrix(G::distG; dp=G.dprocess, x0=G.x0, grid=G.grid, exo=no return transition_matrix(G.model, dp, x0, grid; exo=exo, diff=diff) end + + +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) 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] + Sn_z1 = S_z1[n] + Sn_z2 = S_z2[n] + + Sn = min.(max.(Sn, 0.0),1.0) + qn = div.(Sn, δ) + qn = max.(0, qn) + 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)... ) + + # # # Filling transition matrix + rhs_Π = outer(λn_weight_vector_Π...) + + indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) + + Π[indexes_to_be_modified...] .+= w.*rhs_Π + + @assert d==1 + + λ_vec = (SVector( -1. /δ[1], 1. / δ[1]), ) + A = outer(λ_vec...) + + rhs_dΠ_x = outer2(A, Sn_x[1,:]) + rhs_dΠ_x = [ -1. /δ[1].*Sn_x[1,:] , 1. / δ[1].*Sn_x[1,:]] + dΠ_x[indexes_to_be_modified...] .+= w*rhs_dΠ_x + + rhs_dΠ_z1 = outer2(A, Sn_z1[1,:]) + rhs_dΠ_z1 = [ -1. /δ[1].*Sn_z1[1,:] , 1. / δ[1].*Sn_z1[1,:]] + dΠ_z1[indexes_to_be_modified...] .+= w*rhs_dΠ_z1 + + rhs_dΠ_z2 = outer2(A, Sn_z2[1,:]) + rhs_dΠ_z2 = [ -1. /δ[1].*Sn_z2[1,:] , 1. / δ[1].*Sn_z2[1,:]] + dΠ_z2[indexes_to_be_modified...] .+= w*rhs_dΠ_z2 + + end + +end From e709c1a929e309d7d85fde3699282f76ec2c2259 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 10:56:43 +0100 Subject: [PATCH 03/22] modifying G to return the derivatives wrt z1, z2 --- experiments/dev_G_derivative_wrt_z.jl | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index 30f59e1..51859f0 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -7,12 +7,14 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth return μ1 end - P, P_x = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true) + P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true) μ1 = P'μ0 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) @@ -25,9 +27,23 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth return d_μ end + function fun_z1(dz1::Point{n}) where n + d_μ = (P_z1'*dz1)'*μ0 + return d_μ + end + + function fun_z2(dz2::Point{n}) where n + d_μ = (P_z2'*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 From f177e19e11359b59e1aeec9f26413e4e2571b7f3 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 10:58:00 +0100 Subject: [PATCH 04/22] removing unmodified functions from the dev file --- experiments/dev_G_derivative_wrt_z.jl | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index 51859f0..6d2b43b 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -129,18 +129,6 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing end -function transition_matrix(model, sol; diff=false) - 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); -end - -function transition_matrix(G::distG; dp=G.dprocess, x0=G.x0, grid=G.grid, exo=nothing, diff=false) - - return transition_matrix(G.model, dp, x0, grid; exo=exo, diff=diff) -end - - 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) where d where n_x where _ @assert ndims(Π) == d+1 From a967a70657c400cfc9a730a1c86b8cffe180a92c Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 11:32:45 +0100 Subject: [PATCH 05/22] correcting fun_z --- experiments/dev_G_derivative_wrt_z.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index 6d2b43b..7ebb14d 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -28,12 +28,14 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth end function fun_z1(dz1::Point{n}) where n - d_μ = (P_z1'*dz1)'*μ0 + 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::Point{n}) where n - d_μ = (P_z2'*dz2)'*μ0 + 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 From 877f21c80fb926a3b9252966952a1257a2ad8aaf Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 11:40:15 +0100 Subject: [PATCH 06/22] correcting a composite derivative in transition_matrix --- experiments/dev_G_derivative_wrt_z.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index 7ebb14d..a9a4ff6 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -106,8 +106,8 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing trembling_hand!(Π_view, S, w) else S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] - S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] for n=1:length(S)] * dm_dz - S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] for n=1:length(S)] * dM_dz + S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] for n=1:length(S)] .* dm_dz + S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] for n=1:length(S)] .* dM_dz dΠ_view_x = view(dΠ_x,:,i_m,ind_s...,i_MM) dΠ_view_z1 = view(dΠ_z1,:,i_m,ind_s...,i_MM) From 7dc6fdf561040533fb305a46b70cd6f435729d64 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 13:12:53 +0100 Subject: [PATCH 07/22] correcting a spelling --- experiments/dev_G_derivative_wrt_z.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index a9a4ff6..3077626 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -1,3 +1,6 @@ +using Dolo +import Dolo + function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) where n_x @@ -91,7 +94,7 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing M = inode(Point, dp, i_m, i_M) if !(exo === nothing) M = Dolo.repsvec(exo[2], M) # z1 - dM_dz = Dolo.repsvec((@SVector ones(n_2)),M*0) + dM_dz = Dolo.repsvec((@SVector ones(n_z2)),M*0) end w = iweight(dp, i_m, i_M) if diff From 5784d00a3a50a795f063c165244ae202dbb1f338 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 13:32:40 +0100 Subject: [PATCH 08/22] modifying a product --- experiments/dev_G_derivative_wrt_z.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index 3077626..009c688 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -109,8 +109,8 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing trembling_hand!(Π_view, S, w) else S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] - S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] for n=1:length(S)] .* dm_dz - S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] for n=1:length(S)] .* dM_dz + S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] .* dm_dz for n=1:length(S)] + S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] .* dM_dz for n=1:length(S)] dΠ_view_x = view(dΠ_x,:,i_m,ind_s...,i_MM) dΠ_view_z1 = view(dΠ_z1,:,i_m,ind_s...,i_MM) From 35704fb96b631017e801fa6bfcf3a7f345cf9257 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 16:48:34 +0100 Subject: [PATCH 09/22] test file showing that the derivatives of G wrt z are zero --- experiments/dev_ergodic.jl | 40 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 experiments/dev_ergodic.jl diff --git a/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl new file mode 100644 index 0000000..8585922 --- /dev/null +++ b/experiments/dev_ergodic.jl @@ -0,0 +1,40 @@ +using Dolo +using BenchmarkTools +using StaticArrays +using FiniteDiff + +model = model_rbc_mc +sol = Dolo.improved_time_iteration(model) + +G = Dolo.distG(model, sol) +z10 = SVector{2,Float64}( model.exogenous.values[1,:]) +z20 = SVector{2,Float64}(model.exogenous.values[2,:]) +x0 = G.x0 +μ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) + +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) < 1e-8) + +# print(maximum(abs, Jx_num - Jx_exact) < 1e-8) + +print(maximum(abs, Jz1_num - Jz1_exact) < 1e-8) + +print(maximum(abs, Jz2_num - Jz2_exact) < 1e-8) + + + +maximum(Jz1_num*100000000) + From 6a553e70d407b446f8cea97b9d11c5d13915b073 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 16:51:49 +0100 Subject: [PATCH 10/22] modification --- experiments/dev_ergodic.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl index 8585922..45502b9 100644 --- a/experiments/dev_ergodic.jl +++ b/experiments/dev_ergodic.jl @@ -3,6 +3,8 @@ using BenchmarkTools using StaticArrays using FiniteDiff + +model_rbc_mc = yaml_import("examples/models/rbc_mc.yaml") model = model_rbc_mc sol = Dolo.improved_time_iteration(model) From 22846a16acf28415c96263876dbbfe592ed2f8d8 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 10 Jan 2022 17:34:14 +0100 Subject: [PATCH 11/22] adding modifications --- experiments/dev_G_derivative_wrt_z.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index 009c688..f8117bf 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -1,6 +1,7 @@ using Dolo import Dolo + function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) where n_x @@ -30,13 +31,13 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth return d_μ end - function fun_z1(dz1::Point{n}) where n + 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::Point{n}) where n + 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_μ @@ -82,7 +83,7 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing m = node(exo_grid, i_m) if !(exo === nothing) m = Dolo.repsvec(exo[1], m) # z0 - dm_dz = Dolo.repsvec((@SVector ones(n_z1)),m*0) + dm_dz = Dolo.repsvec((@SVector ones(length(exo[1]))),m*0) end for i_M in 1:n_inodes(dp, i_m) @@ -94,7 +95,7 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing M = inode(Point, dp, i_m, i_M) if !(exo === nothing) M = Dolo.repsvec(exo[2], M) # z1 - dM_dz = Dolo.repsvec((@SVector ones(n_z2)),M*0) + dM_dz = Dolo.repsvec((@SVector ones(length(exo[2]))),M*0) end w = iweight(dp, i_m, i_M) if diff @@ -109,8 +110,8 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing trembling_hand!(Π_view, S, w) else S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] - S_z1 = [( 1.0 ./(b-a)) .* S_z1[n] .* dm_dz for n=1:length(S)] - S_z2 = [( 1.0 ./(b-a)) .* S_z2[n] .* dM_dz for n=1:length(S)] + S_z1 = [( 1.0 ./(b-a)) .* S_z1[n].* dm_dz for n=1:length(S)] + S_z2 = [( 1.0 ./(b-a)) .* S_z2[n].* dM_dz for n=1:length(S)] dΠ_view_x = view(dΠ_x,:,i_m,ind_s...,i_MM) dΠ_view_z1 = view(dΠ_z1,:,i_m,ind_s...,i_MM) @@ -165,7 +166,7 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve Π[indexes_to_be_modified...] .+= w.*rhs_Π - @assert d==1 + #@assert d==1 λ_vec = (SVector( -1. /δ[1], 1. / δ[1]), ) A = outer(λ_vec...) @@ -184,4 +185,4 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve end -end +end \ No newline at end of file From 5f3648070a434c4694ebf437eb8eb6a0f9e7b7aa Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Thu, 13 Jan 2022 10:31:56 +0100 Subject: [PATCH 12/22] =?UTF-8?q?correcting=20the=20computation=20of=20?= =?UTF-8?q?=E2=88=82g=5F=E2=88=82z1=20and=20=E2=88=82g=5F=E2=88=82z2?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- experiments/dev_G_derivative_wrt_z.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl index f8117bf..d626688 100644 --- a/experiments/dev_G_derivative_wrt_z.jl +++ b/experiments/dev_G_derivative_wrt_z.jl @@ -83,7 +83,6 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing m = node(exo_grid, i_m) if !(exo === nothing) m = Dolo.repsvec(exo[1], m) # z0 - dm_dz = Dolo.repsvec((@SVector ones(length(exo[1]))),m*0) end for i_M in 1:n_inodes(dp, i_m) @@ -95,7 +94,6 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing M = inode(Point, dp, i_m, i_M) if !(exo === nothing) M = Dolo.repsvec(exo[2], M) # z1 - dM_dz = Dolo.repsvec((@SVector ones(length(exo[2]))),M*0) end w = iweight(dp, i_m, i_M) if diff @@ -110,8 +108,11 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing trembling_hand!(Π_view, S, w) else S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] - S_z1 = [( 1.0 ./(b-a)) .* S_z1[n].* dm_dz for n=1:length(S)] - S_z2 = [( 1.0 ./(b-a)) .* S_z2[n].* dM_dz for n=1:length(S)] + 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 = S_z1[:,1:length(exo[1])] + S_z2 = S_z2[:,1:length(exo[2])] dΠ_view_x = view(dΠ_x,:,i_m,ind_s...,i_MM) dΠ_view_z1 = view(dΠ_z1,:,i_m,ind_s...,i_MM) From efea269a922037282b860dd02c45d8e93f4f07ad Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Thu, 13 Jan 2022 15:18:26 +0100 Subject: [PATCH 13/22] little modifications --- experiments/dev_ergodic.jl | 33 +-- src/algos/ergodic.jl | 444 ++++++++++++++++++++++++++----------- 2 files changed, 340 insertions(+), 137 deletions(-) diff --git a/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl index 45502b9..b37ca9e 100644 --- a/experiments/dev_ergodic.jl +++ b/experiments/dev_ergodic.jl @@ -4,39 +4,46 @@ using StaticArrays using FiniteDiff -model_rbc_mc = yaml_import("examples/models/rbc_mc.yaml") -model = model_rbc_mc +model = yaml_import("examples/models/rbc_iid.yaml") sol = Dolo.improved_time_iteration(model) + G = Dolo.distG(model, sol) -z10 = SVector{2,Float64}( model.exogenous.values[1,:]) -z20 = SVector{2,Float64}(model.exogenous.values[2,:]) +z10 = SVector(model.calibration[:exogenous]...) # for rbc_iid +z20 = z10 # for rbc_iid 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) +z10 + -# Jμ_exact = convert(Matrix, ∂G_∂μ) -# Jμ_num = FiniteDiff.finite_difference_jacobian(mu -> G(mu, x0, exo = [z10,z20]), μ0) +μ1, ∂G_∂μ, ∂G_∂x, ∂G_∂z1, ∂G_∂z2 = G(μ0, x0; exo = [z10,z20], diff = true) -# Jx_exact = convert(Matrix, ∂G_∂x) -# Jx_num = FiniteDiff.finite_difference_jacobian(x -> G(μ0, x, exo = [z10,z20]), x0) +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) +Jz2_num = FiniteDiff.finite_difference_jacobian(z2 -> G(μ0, x0; exo = [z10,z2]), z20) -# print( maximum(abs, Jμ_num - Jμ_exact) < 1e-8) +print( maximum(abs, Jμ_num - Jμ_exact) < 1e-8) -# print(maximum(abs, Jx_num - Jx_exact) < 1e-8) +print(maximum(abs, Jx_num - Jx_exact)) print(maximum(abs, Jz1_num - Jz1_exact) < 1e-8) -print(maximum(abs, Jz2_num - Jz2_exact) < 1e-8) +print(maximum(abs, Jz2_num - Jz2_exact)) maximum(Jz1_num*100000000) +print( (abs.(Jz2_exact-Jz2_num).>1e-8)) + +G(μ0, x0, exo = [z10,z20]) diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index 6ad2253..2fc51af 100644 --- a/src/algos/ergodic.jl +++ b/src/algos/ergodic.jl @@ -59,44 +59,44 @@ end -function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; diff=false) where n_x +# function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; diff=false) where n_x - if !diff - P = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=nothing, diff=false) - μ1 = P'*μ0 - return μ1 - end +# if !diff +# P = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=nothing, 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 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=nothing, diff=true) - μ1 = P'μ0 +# μ1 = P'μ0 - M = length(μ0) - N = length(x0.data)*length(x0.data[1]) - - function fun_x(dx::AbstractVector{Float64}) - d_x = MSM(copy(reinterpret(SVector{n_x, Float64}, dx)), x0.sizes) - N = size(P,2) - d_μ = [ sum([μ0[k]*(P_x[k,i]'*d_x.data[k]) for k=1:N]) for i=1:N ] - # alternate calculation - P_dx = [(P_x[i,j]'*d_x.data[i]) for i=1:size(P,1), j=1:size(P,2)] - d_μμ = P_dx'*μ0 - @assert (maximum(abs, d_μ-d_μμ)<1e-10) - return d_μ - end - - ∂G_∂μ = LinearMap( μ -> P'*μ, M, M) - ∂G_∂x = LinearMap(fun_x, M, N) - return μ1, ∂G_∂μ, ∂G_∂x +# M = length(μ0) +# N = length(x0.data)*length(x0.data[1]) + +# function fun_x(dx::AbstractVector{Float64}) +# d_x = MSM(copy(reinterpret(SVector{n_x, Float64}, dx)), x0.sizes) +# N = size(P,2) +# d_μ = [ sum([μ0[k]*(P_x[k,i]'*d_x.data[k]) for k=1:N]) for i=1:N ] +# # alternate calculation +# P_dx = [(P_x[i,j]'*d_x.data[i]) for i=1:size(P,1), j=1:size(P,2)] +# d_μμ = P_dx'*μ0 +# @assert (maximum(abs, d_μ-d_μμ)<1e-10) +# return d_μ +# end + +# ∂G_∂μ = LinearMap( μ -> P'*μ, M, M) +# ∂G_∂x = LinearMap(fun_x, M, N) +# return μ1, ∂G_∂μ, ∂G_∂x -end +# 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) 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) end @@ -144,70 +144,70 @@ Calculates the new transition matrix for a given model, a given discretized exog * `Π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) where n_x - # the restriction here is that dp is a descrete_process +# # the restriction here is that dp is a descrete_process - parms = SVector(model.calibration[:parameters]...) +# parms = SVector(model.calibration[:parameters]...) - exo_grid = grid.exo - endo_grid = grid.endo - - N_m = max(1, n_nodes(exo_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) - end - s = nodes(endo_grid) - a = SVector(endo_grid.min...) - b = SVector(endo_grid.max...) - for i_m in 1:max(1, n_nodes(exo_grid)) - x = x0.views[i_m] - m = node(exo_grid, i_m) - if !(exo === nothing) - m = Dolo.repsvec(exo[1], m) # z0 - end - for i_M in 1:n_inodes(dp, i_m) +# exo_grid = grid.exo +# endo_grid = grid.endo + +# N_m = max(1, n_nodes(exo_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) +# end +# s = nodes(endo_grid) +# a = SVector(endo_grid.min...) +# b = SVector(endo_grid.max...) +# for i_m in 1:max(1, n_nodes(exo_grid)) +# x = x0.views[i_m] +# m = node(exo_grid, i_m) +# if !(exo === nothing) +# m = Dolo.repsvec(exo[1], m) # z0 +# end +# for i_M in 1:n_inodes(dp, i_m) - if isa(grid.exo, EmptyGrid) - i_MM = 1 # used to index transition matrix - else - i_MM = i_M - end - M = inode(Point, dp, i_m, i_M) - 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) - else - S = transition(model, m, s, x, M, parms) - end - S = [(S[n]-a)./(b-a) for n=1:length(S)] - 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) - 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) - end - end - end - - Π0 = (reshape(Π,N,N)) - if !diff - return Π0 - else - dΠ0 = reshape(dΠ,N,N) - return Π0, dΠ0 - end +# if isa(grid.exo, EmptyGrid) +# i_MM = 1 # used to index transition matrix +# else +# i_MM = i_M +# end +# M = inode(Point, dp, i_m, i_M) +# 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) +# else +# S = transition(model, m, s, x, M, parms) +# end +# S = [(S[n]-a)./(b-a) for n=1:length(S)] +# 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) +# 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) +# end +# end +# end + +# Π0 = (reshape(Π,N,N)) +# if !diff +# return Π0 +# else +# dΠ0 = reshape(dΠ,N,N) +# return Π0, dΠ0 +# end -end +# end function transition_matrix(model, sol; diff=false) @@ -263,53 +263,53 @@ 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 _ +# function trembling_foot!(Π, dΠ, S::Vector{Point{d}}, S_x::Vector{SMatrix{d,n_x,Float64,_}}, w::Float64) 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] - - Sn = min.(max.(Sn, 0.0),1.0) - qn = div.(Sn, δ) - qn = max.(0, qn) - qn = min.(qn, shape_Π[2:d+1].-2) - λn = (Sn./δ.-qn) # ∈[0,1[ by construction - qn_ = round.(Int,qn) .+ 1 +# @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] + +# Sn = min.(max.(Sn, 0.0),1.0) +# qn = div.(Sn, δ) +# qn = max.(0, qn) +# 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)... ) +# λn_weight_vector_Π = tuple( (SVector((1-λn[i]),λn[i]) for i in 1:d)... ) - # # # Filling transition matrix - rhs_Π = outer(λn_weight_vector_Π...) +# # # # Filling transition matrix +# rhs_Π = outer(λn_weight_vector_Π...) - indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) +# indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) - Π[indexes_to_be_modified...] .+= w.*rhs_Π +# Π[indexes_to_be_modified...] .+= w.*rhs_Π - # for k=1:d - # λ_vec = tuple( (i==k ? SVector( -1. /δ[k], 1. / δ[k]) : (SVector((1-λn[i]),λn[i])) for i in 1:d)... ) - # A = outer(λ_vec...) - # rhs_dΠ = outer2(A, Sn_x[k,:]) - # dΠ[indexes_to_be_modified...] .+= w*rhs_dΠ - # end +# # for k=1:d +# # λ_vec = tuple( (i==k ? SVector( -1. /δ[k], 1. / δ[k]) : (SVector((1-λn[i]),λn[i])) for i in 1:d)... ) +# # A = outer(λ_vec...) +# # rhs_dΠ = outer2(A, Sn_x[k,:]) +# # dΠ[indexes_to_be_modified...] .+= w*rhs_dΠ +# # end - @assert d==1 +# @assert d==1 - λ_vec = (SVector( -1. /δ[1], 1. / δ[1]), ) - A = outer(λ_vec...) - rhs_dΠ = outer2(A, Sn_x[1,:]) - rhs_dΠ = [ -1. /δ[1].*Sn_x[1,:] , 1. / δ[1].*Sn_x[1,:]] - dΠ[indexes_to_be_modified...] .+= w*rhs_dΠ +# λ_vec = (SVector( -1. /δ[1], 1. / δ[1]), ) +# A = outer(λ_vec...) +# rhs_dΠ = outer2(A, Sn_x[1,:]) +# rhs_dΠ = [ -1. /δ[1].*Sn_x[1,:] , 1. / δ[1].*Sn_x[1,:]] +# dΠ[indexes_to_be_modified...] .+= w*rhs_dΠ - end +# end -end +# end @@ -340,3 +340,199 @@ Computes an outer product between a vector and a matrix and returns a vector of """ outer2(A, x) = [A[i]*x for i in CartesianIndices(A)] +function outer3(δ_weight_vector) + +end + + + + + + +### dev + + +function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) where n_x + + + if !diff + P = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=false) + μ1 = P'*μ0 + return μ1 + end + + P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true) + + μ1 = P'μ0 + + 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) + N = size(P,2) + d_μ = [ sum([μ0[k]*(P_x[k,i]'*d_x.data[k]) for k=1:N]) for i=1:N ] + # alternate calculation + P_dx = [(P_x[i,j]'*d_x.data[i]) for i=1:size(P,1), j=1:size(P,2)] + d_μμ = P_dx'*μ0 + @assert (maximum(abs, d_μ-d_μμ)<1e-10) + 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) + ∂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 transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing, diff=false) where n_x + + # the restriction here is that dp is a descrete_process + + parms = SVector(model.calibration[:parameters]...) + + exo_grid = grid.exo + endo_grid = grid.endo + + N_m = max(1, n_nodes(exo_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Π_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...) + b = SVector(endo_grid.max...) + for i_m in 1:max(1, n_nodes(exo_grid)) + x = x0.views[i_m] + m = node(exo_grid, i_m) + if !(exo === nothing) + m = Dolo.repsvec(exo[1], m) # z0 + end + for i_M in 1:n_inodes(dp, i_m) + + if isa(grid.exo, EmptyGrid) + i_MM = 1 # used to index transition matrix + else + i_MM = i_M + end + M = inode(Point, dp, i_m, i_M) + if !(exo === nothing) + M = Dolo.repsvec(exo[2], M) # z1 + end + w = iweight(dp, i_m, i_M) + if diff + S, S_z1, S_x, S_z2 = transition(model, Val{(0,1,3,4)}, m, s, x, M, parms) + else + S = transition(model, m, s, x, M, parms) + end + S = [(S[n]-a)./(b-a) for n=1:length(S)] + 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) + else + S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] + 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) + end + end + end + + Π0 = (reshape(Π,N,N)) + if !diff + return Π0 + else + 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 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) 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] + Sn_z1 = S_z1[n] + Sn_z2 = S_z2[n] + + Sn = min.(max.(Sn, 0.0),1.0) + qn = div.(Sn, δ) + qn = max.(0, qn) + 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)... ) + + # # # Filling transition matrix + rhs_Π = outer(λn_weight_vector_Π...) + + indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) + + Π[indexes_to_be_modified...] .+= w.*rhs_Π + + + for k=1:d + λ_vec = tuple( (i==k ? SVector( -1. /δ[k], 1. / δ[k]) : (SVector((1-λn[i]),λn[i])) for i in 1:d)... ) + A = outer(λ_vec...) + + 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 + From ea1067f11327c4427835694da0ae6687df8a242b Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 17 Jan 2022 08:16:28 +0100 Subject: [PATCH 14/22] =?UTF-8?q?pb=20with=20the=20computation=20of=20?= =?UTF-8?q?=E2=88=82G=5F=E2=88=82z=20in=20rbc=5Fiid?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- experiments/dev_ergodic.jl | 10 +++-- src/algos/ergodic.jl | 83 +++++++++++++++++++++++++++++--------- 2 files changed, 69 insertions(+), 24 deletions(-) diff --git a/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl index b37ca9e..345714b 100644 --- a/experiments/dev_ergodic.jl +++ b/experiments/dev_ergodic.jl @@ -15,7 +15,6 @@ x0 = G.x0 x0_flat = cat(G.x0.data...; dims=1) μ0 = G.μ0 -z10 μ1, ∂G_∂μ, ∂G_∂x, ∂G_∂z1, ∂G_∂z2 = G(μ0, x0; exo = [z10,z20], diff = true) @@ -41,9 +40,12 @@ print(maximum(abs, Jz1_num - Jz1_exact) < 1e-8) print(maximum(abs, Jz2_num - Jz2_exact)) +# using Plots -maximum(Jz1_num*100000000) +# spy((abs.(Jz2_exact-Jz2_num).>1e-8)) +# maximum(Jz1_num*100000000) -print( (abs.(Jz2_exact-Jz2_num).>1e-8)) +# print( (abs.(Jz2_exact-Jz2_num).>1e-8)) -G(μ0, x0, exo = [z10,z20]) +# print(Jz2_num) +# G(μ0, x0, exo = [z10,z20]) diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index 687a857..d03ad8c 100644 --- a/src/algos/ergodic.jl +++ b/src/algos/ergodic.jl @@ -59,21 +59,22 @@ end -function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; diff=false) where n_x - +function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) 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) μ1 = P'μ0 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,17 +87,33 @@ 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) 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) end @@ -156,7 +173,13 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing 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) + 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...) @@ -180,7 +203,7 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing 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) else S = transition(model, m, s, x, M, parms) end @@ -191,8 +214,17 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing trembling_hand!(Π_view, S, w) 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) end end end @@ -201,8 +233,10 @@ 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 @@ -259,7 +293,7 @@ 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 _ +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) where d where n_x where _ @assert ndims(Π) == d+1 shape_Π = size(Π) @@ -271,13 +305,14 @@ function trembling_foot!(Π, dΠ, S::Vector{Point{d}}, S_x::Vector{SMatrix{d,n_x Sn = S[n] Sn_x = S_x[n] - + Sn_z1 = S_z1[n] + Sn_z2 = S_z2[n] + inbound = (0. .<= Sn) .& ( Sn .<= 1.0) - + 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 @@ -290,14 +325,22 @@ 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_Π + 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)... ) 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 From a1dc351a23558fde0199e4d2beedf75569993777 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 31 Jan 2022 11:12:00 +0100 Subject: [PATCH 15/22] some graphs explaining that the mistake in the derivation wrt z occurs with a repeted pattern of difference with the numerical value, every 15 values on the example of rbc_iid --- experiments/dev_ergodic.jl | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl index 345714b..083c8a0 100644 --- a/experiments/dev_ergodic.jl +++ b/experiments/dev_ergodic.jl @@ -40,12 +40,24 @@ print(maximum(abs, Jz1_num - Jz1_exact) < 1e-8) print(maximum(abs, Jz2_num - Jz2_exact)) -# using Plots +using Plots -# spy((abs.(Jz2_exact-Jz2_num).>1e-8)) -# maximum(Jz1_num*100000000) +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) -# print( (abs.(Jz2_exact-Jz2_num).>1e-8)) +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)) -# print(Jz2_num) -# G(μ0, x0, exo = [z10,z20]) + +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)) + +Jz2_num \ No newline at end of file From 833845bdafe8b0f21113f02fb80873a4c99a5c96 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Thu, 3 Feb 2022 16:30:40 +0100 Subject: [PATCH 16/22] some temporary indications about the sizes --- src/algos/ergodic.jl | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index d03ad8c..f03cd8f 100644 --- a/src/algos/ergodic.jl +++ b/src/algos/ergodic.jl @@ -69,6 +69,9 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true) + #dev + #P_z2 = @SMatrix [reshape([P_z2[:,15*i+j] for j=1:15, i=0:14]',1,225)[i][j] for j=1:225, i=1:225] + μ1 = P'μ0 M = length(μ0) @@ -93,8 +96,12 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth return d_μ end + print(" typeof Pz2 ",typeof(P_z2)," ", "size of Pz2", size(P_z2)," ") + print(" typeof P ",typeof(P)," ", "size of P", size(P)," ") + function fun_z2(dz2) P_dz2 = [(P_z2[i,j]'*dz2) for i=1:size(P,1), j=1:size(P,2)] + print("typeof Pdz2 ",typeof(P_dz2)," size of P_dz2 ", size(P_dz2)," ") d_μ = P_dz2'*μ0 return d_μ end @@ -200,6 +207,7 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing M = inode(Point, dp, i_m, i_M) if !(exo === nothing) M = Dolo.repsvec(exo[2], M) # z1 + #println("M : ",M," ") end w = iweight(dp, i_m, i_M) if diff @@ -207,23 +215,34 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing else S = transition(model, m, s, x, M, parms) end + print("S_z2 = ", S_z2," ") S = [(S[n]-a)./(b-a) for n=1:length(S)] 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) else + #print("type of Sz1 0 : ",typeof(S_z1)," length of Sz1 0 : ", length(S_z1), " ") # Vector{SMatrix{2, 1, Float64, 2}} of length 225 + #print("type of Sx 0 : ",typeof(S_x)," length of Sx 0 : ", length(S_x), " ") # Vector{SMatrix{2, 2, Float64, 4}} of length 225 + S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] 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)] + #print("type of Sz1 1 : ",typeof(S_z1)," length of Sz1 1 : ", length(S_z1), " ") # Vector{SMatrix{2, 1, Float64, 2}} of length 225 + #print("type of Sx 1 : ",typeof(S_x)," length of Sx 1 : ", length(S_x), " ") # Vector{SMatrix{2, 2, Float64, 4}} of length 225 + S_z1 = [M[:,1:length(exo[1])] for M in S_z1] - S_z2 = [M[:,1:length(exo[2])] for M in S_z2] + S_z2 = [M[:,1:length(exo[2])] for M in S_z2] + + #print("type of Sz1 2 : ",typeof(S_z1)," length of Sz1 2 : ", length(S_z1), " size of Sz1[1] 2",size(S_z1[1])," ") # Vector{Matrix{Float64}} of length 225 containing matrices of size (2,1) for rbc_iid 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) + print("dΠdΠ_view_z2_z2 : ", dΠ_view_z2," ") + trembling_foot!(Π_view, dΠ_view_x, dΠ_view_z1, dΠ_view_z2, S, S_x, S_z1, S_z2, w) end end @@ -300,6 +319,9 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve grid_dimension = d δ = SVector{d,Float64}(1.0./(shape_Π[1+i]-1) for i in 1:d ) N = shape_Π[1] + + #print("typeof Snx : ", typeof(S_x[1]), " ") # SMatrix{2, 2, Float64, 4} + #print("typeof Snz2 : ", typeof(S_z2[1]), "size of Snz2 : ",size(S_z2[1])," ") # Matrix{Float64} (2,1) for n in 1:N @@ -326,21 +348,28 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve Π[indexes_to_be_modified...] .+= w.*rhs_Π - + #print("d = ",d," ") # 2 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)... ) A = outer(λ_vec...) rhs_dΠ_x = outer2(A, Sn_x[k,:]) dΠ_x[indexes_to_be_modified...] .+= w*rhs_dΠ_x + + #print("size of rhs_dΠ_x : ", size(rhs_dΠ_x)," typeof rhs_dΠ_x", typeof(rhs_dΠ_x)," ") #SizedMatrix{2, 2, SVector{2, Float64}, 2, Matrix{SVector{2, Float64}}} + 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 + #print("size of rhs_dΠ_z2 : ", size(rhs_dΠ_z2)," typeof rhs_dΠ_z2", typeof(rhs_dΠ_z2)," ") #SizedMatrix{2, 2, Vector{Float64}, 2, Matrix{Vector{Float64}}} + end + # dΠ_z2_view = view(dΠ_z2,n,1,:,:,1) + # dΠ_z2_view = view(dΠ_z2,n,1,:,:,1)' end end @@ -358,7 +387,6 @@ 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. From b7d7c814b25e7ebd46a92458cd1ae67c8ed4f1a4 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 7 Mar 2022 22:54:28 +0100 Subject: [PATCH 17/22] modifiing trembling_hand and trembling_foot to smooth the distribution --- src/algos/ergodic.jl | 148 ++++++++++++++++++++++--------------------- 1 file changed, 76 insertions(+), 72 deletions(-) diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index f03cd8f..9012201 100644 --- a/src/algos/ergodic.jl +++ b/src/algos/ergodic.jl @@ -59,7 +59,7 @@ end -function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) where n_x +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=exo, diff=false) @@ -67,10 +67,8 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth return μ1 end - P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true) + P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true, smooth=smooth) - #dev - #P_z2 = @SMatrix [reshape([P_z2[:,15*i+j] for j=1:15, i=0:14]',1,225)[i][j] for j=1:225, i=1:225] μ1 = P'μ0 @@ -96,12 +94,8 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth return d_μ end - print(" typeof Pz2 ",typeof(P_z2)," ", "size of Pz2", size(P_z2)," ") - print(" typeof P ",typeof(P)," ", "size of P", size(P)," ") - function fun_z2(dz2) P_dz2 = [(P_z2[i,j]'*dz2) for i=1:size(P,1), j=1:size(P,2)] - print("typeof Pdz2 ",typeof(P_dz2)," size of P_dz2 ", size(P_dz2)," ") d_μ = P_dz2'*μ0 return d_μ end @@ -116,42 +110,42 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth end -function (G::distG)(μ0::AbstractVector{Float64}, x0::AbstractVector{Float64}; exo=nothing, 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; exo=exo, 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 +# """ +# 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). @@ -166,7 +160,7 @@ Calculates the new transition matrix for a given model, a given discretized exog * `Π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 @@ -207,7 +201,6 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing M = inode(Point, dp, i_m, i_M) if !(exo === nothing) M = Dolo.repsvec(exo[2], M) # z1 - #println("M : ",M," ") end w = iweight(dp, i_m, i_M) if diff @@ -215,35 +208,26 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing else S = transition(model, m, s, x, M, parms) end - print("S_z2 = ", S_z2," ") S = [(S[n]-a)./(b-a) for n=1:length(S)] 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 - #print("type of Sz1 0 : ",typeof(S_z1)," length of Sz1 0 : ", length(S_z1), " ") # Vector{SMatrix{2, 1, Float64, 2}} of length 225 - #print("type of Sx 0 : ",typeof(S_x)," length of Sx 0 : ", length(S_x), " ") # Vector{SMatrix{2, 2, Float64, 4}} of length 225 - + S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] 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)] - #print("type of Sz1 1 : ",typeof(S_z1)," length of Sz1 1 : ", length(S_z1), " ") # Vector{SMatrix{2, 1, Float64, 2}} of length 225 - #print("type of Sx 1 : ",typeof(S_x)," length of Sx 1 : ", length(S_x), " ") # Vector{SMatrix{2, 2, Float64, 4}} of length 225 - S_z1 = [M[:,1:length(exo[1])] for M in S_z1] S_z2 = [M[:,1:length(exo[2])] for M in S_z2] - #print("type of Sz1 2 : ",typeof(S_z1)," length of Sz1 2 : ", length(S_z1), " size of Sz1[1] 2",size(S_z1[1])," ") # Vector{Matrix{Float64}} of length 225 containing matrices of size (2,1) for rbc_iid - + 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) - print("dΠdΠ_view_z2_z2 : ", dΠ_view_z2," ") - - trembling_foot!(Π_view, dΠ_view_x, dΠ_view_z1, dΠ_view_z2, S, S_x, S_z1, S_z2, w) + 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 @@ -261,15 +245,30 @@ function transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing 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 + +""" +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 + +function ∂smoothing(x::SVector{2}) + return 6. .* x .^2 .- 6. .* x end """ @@ -281,7 +280,7 @@ Updates A. # 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) @@ -299,8 +298,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 @@ -312,7 +315,8 @@ function trembling_hand!(A, x::Vector{Point{d}}, w::Float64) where d end -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) where d where n_x where _ + +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(Π) @@ -320,8 +324,6 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve δ = SVector{d,Float64}(1.0./(shape_Π[1+i]-1) for i in 1:d ) N = shape_Π[1] - #print("typeof Snx : ", typeof(S_x[1]), " ") # SMatrix{2, 2, Float64, 4} - #print("typeof Snz2 : ", typeof(S_z2[1]), "size of Snz2 : ",size(S_z2[1])," ") # Matrix{Float64} (2,1) for n in 1:N @@ -339,7 +341,11 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve λ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)... ) + else + λn_weight_vector_Π = tuple( (smoothing(SVector((λn[i]),1-λn[i])) for i in 1:d)... ) + end # # # Filling transition matrix rhs_Π = outer(λn_weight_vector_Π...) @@ -348,28 +354,26 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve Π[indexes_to_be_modified...] .+= w.*rhs_Π - #print("d = ",d," ") # 2 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Π_x = outer2(A, Sn_x[k,:]) dΠ_x[indexes_to_be_modified...] .+= w*rhs_dΠ_x - - #print("size of rhs_dΠ_x : ", size(rhs_dΠ_x)," typeof rhs_dΠ_x", typeof(rhs_dΠ_x)," ") #SizedMatrix{2, 2, SVector{2, Float64}, 2, Matrix{SVector{2, Float64}}} - 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 - #print("size of rhs_dΠ_z2 : ", size(rhs_dΠ_z2)," typeof rhs_dΠ_z2", typeof(rhs_dΠ_z2)," ") #SizedMatrix{2, 2, Vector{Float64}, 2, Matrix{Vector{Float64}}} - + end - # dΠ_z2_view = view(dΠ_z2,n,1,:,:,1) - # dΠ_z2_view = view(dΠ_z2,n,1,:,:,1)' end end From 75d14358b51f53b5f7ebe19c12a415f0f17ba4dd Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Thu, 10 Mar 2022 09:59:31 +0100 Subject: [PATCH 18/22] removing a file --- experiments/dev_G_derivative_wrt_z.jl | 189 -------------------------- 1 file changed, 189 deletions(-) delete mode 100644 experiments/dev_G_derivative_wrt_z.jl diff --git a/experiments/dev_G_derivative_wrt_z.jl b/experiments/dev_G_derivative_wrt_z.jl deleted file mode 100644 index d626688..0000000 --- a/experiments/dev_G_derivative_wrt_z.jl +++ /dev/null @@ -1,189 +0,0 @@ -using Dolo -import Dolo - - -function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =nothing, diff=false) where n_x - - - if !diff - P = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=false) - μ1 = P'*μ0 - return μ1 - end - - P, P_x, P_z1, P_z2 = transition_matrix(G.model, G.dprocess, x0, G.grid; exo=exo, diff=true) - - μ1 = P'μ0 - - 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) - N = size(P,2) - d_μ = [ sum([μ0[k]*(P_x[k,i]'*d_x.data[k]) for k=1:N]) for i=1:N ] - # alternate calculation - P_dx = [(P_x[i,j]'*d_x.data[i]) for i=1:size(P,1), j=1:size(P,2)] - d_μμ = P_dx'*μ0 - @assert (maximum(abs, d_μ-d_μμ)<1e-10) - 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) - ∂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 transition_matrix(model, dp, x0::MSM{<:SVector{n_x}}, grid; exo=nothing, diff=false) where n_x - - # the restriction here is that dp is a descrete_process - - parms = SVector(model.calibration[:parameters]...) - - exo_grid = grid.exo - endo_grid = grid.endo - - N_m = max(1, n_nodes(exo_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Π_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...) - b = SVector(endo_grid.max...) - for i_m in 1:max(1, n_nodes(exo_grid)) - x = x0.views[i_m] - m = node(exo_grid, i_m) - if !(exo === nothing) - m = Dolo.repsvec(exo[1], m) # z0 - end - for i_M in 1:n_inodes(dp, i_m) - - if isa(grid.exo, EmptyGrid) - i_MM = 1 # used to index transition matrix - else - i_MM = i_M - end - M = inode(Point, dp, i_m, i_M) - if !(exo === nothing) - M = Dolo.repsvec(exo[2], M) # z1 - end - w = iweight(dp, i_m, i_M) - if diff - S, S_z1, S_x, S_z2 = transition(model, Val{(0,1,3,4)}, m, s, x, M, parms) - else - S = transition(model, m, s, x, M, parms) - end - S = [(S[n]-a)./(b-a) for n=1:length(S)] - 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) - else - S_x = [( 1.0 ./(b-a)) .* S_x[n] for n=1:length(S)] - 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 = S_z1[:,1:length(exo[1])] - S_z2 = S_z2[:,1:length(exo[2])] - - 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) - end - end - end - - Π0 = (reshape(Π,N,N)) - if !diff - return Π0 - else - 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 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) 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] - Sn_z1 = S_z1[n] - Sn_z2 = S_z2[n] - - Sn = min.(max.(Sn, 0.0),1.0) - qn = div.(Sn, δ) - qn = max.(0, qn) - 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)... ) - - # # # Filling transition matrix - rhs_Π = outer(λn_weight_vector_Π...) - - indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) - - Π[indexes_to_be_modified...] .+= w.*rhs_Π - - #@assert d==1 - - λ_vec = (SVector( -1. /δ[1], 1. / δ[1]), ) - A = outer(λ_vec...) - - rhs_dΠ_x = outer2(A, Sn_x[1,:]) - rhs_dΠ_x = [ -1. /δ[1].*Sn_x[1,:] , 1. / δ[1].*Sn_x[1,:]] - dΠ_x[indexes_to_be_modified...] .+= w*rhs_dΠ_x - - rhs_dΠ_z1 = outer2(A, Sn_z1[1,:]) - rhs_dΠ_z1 = [ -1. /δ[1].*Sn_z1[1,:] , 1. / δ[1].*Sn_z1[1,:]] - dΠ_z1[indexes_to_be_modified...] .+= w*rhs_dΠ_z1 - - rhs_dΠ_z2 = outer2(A, Sn_z2[1,:]) - rhs_dΠ_z2 = [ -1. /δ[1].*Sn_z2[1,:] , 1. / δ[1].*Sn_z2[1,:]] - dΠ_z2[indexes_to_be_modified...] .+= w*rhs_dΠ_z2 - - end - -end \ No newline at end of file From 12b18d520da51532e5ff47e1e777eb3723073c3b Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Thu, 10 Mar 2022 17:01:10 +0100 Subject: [PATCH 19/22] commenting the functions of ergodic.jl + adding a test file of the derivatives + adding plots in dev_ergodic.jl about the difference in terms of ergodic distribution with or without smoothing option --- examples/models/consumption_savings_iid.yaml | 2 +- examples/models/rbc_mc.yaml | 7 +- experiments/dev_ergodic.jl | 148 ++++++++++++++++++- src/algos/ergodic.jl | 133 ++++++++++------- test/test_G_derivative.jl | 38 +++-- 5 files changed, 252 insertions(+), 76 deletions(-) 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/experiments/dev_ergodic.jl b/experiments/dev_ergodic.jl index 083c8a0..7ac2207 100644 --- a/experiments/dev_ergodic.jl +++ b/experiments/dev_ergodic.jl @@ -4,13 +4,13 @@ using StaticArrays using FiniteDiff -model = yaml_import("examples/models/rbc_iid.yaml") +model = yaml_import("examples/models/rbc.yaml") sol = Dolo.improved_time_iteration(model) G = Dolo.distG(model, sol) -z10 = SVector(model.calibration[:exogenous]...) # for rbc_iid -z20 = z10 # for rbc_iid +z10 = SVector(model.calibration[:exogenous]...) +z20 = z10 x0 = G.x0 x0_flat = cat(G.x0.data...; dims=1) μ0 = G.μ0 @@ -31,15 +31,15 @@ Jz1_num = FiniteDiff.finite_difference_jacobian(z1 -> G(μ0, x0; exo = [z1,z20]) 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) < 1e-8) +print( maximum(abs, Jμ_num - Jμ_exact) ) print(maximum(abs, Jx_num - Jx_exact)) -print(maximum(abs, Jz1_num - Jz1_exact) < 1e-8) +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") @@ -60,4 +60,138 @@ 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)) -Jz2_num \ No newline at end of file + +# 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_iid.yaml") +sol = Dolo.time_iteration(model) + +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)) + +df = DataFrame(w=vcat(w,w), μ=vcat(μ_smoothing[1:n], μ_no_smoothing[1:n]), smooth=vcat(["smoothing" for i in 1:n], ["no smoothing" for i in 1:n])) + +PlotlyJS.plot( + PlotlyJS.scatter(df, x=:w, y=:μ, group=:smooth), + Layout( + xaxis_title="w", + yaxis_title="μ" + ) + ) + + + diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index 9012201..9ab9d21 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,7 +66,24 @@ function ergodic_distribution(G::distG; x0=G.x0, kwargs...) end - +""" +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 @@ -70,7 +95,7 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::MSM{Point{n_x}}; exo =noth 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]) @@ -119,33 +144,6 @@ function (G::distG)(μ0::AbstractVector{Float64}, x0::AbstractVector{Float64}; e 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). @@ -153,9 +151,11 @@ 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. """ @@ -170,10 +170,10 @@ 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 + Π = 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]) @@ -189,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) @@ -198,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_z1, S_x, S_z2 = transition(model, Val{(0,1,3,4)}, 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 @@ -267,6 +268,13 @@ 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 @@ -277,6 +285,8 @@ Updates A. * `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 """ @@ -315,7 +325,27 @@ function trembling_hand!(A, x::Vector{Point{d}}, w::Float64; smooth=true) where end - +""" +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` : +# 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 @@ -332,7 +362,7 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve Sn_z1 = S_z1[n] Sn_z2 = S_z2[n] - inbound = (0. .<= Sn) .& ( Sn .<= 1.0) + 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, δ) @@ -342,9 +372,9 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve qn_ = round.(Int,qn) .+ 1 if smooth==false - λn_weight_vector_Π = tuple( (SVector((1-λn[i]),λn[i]) for i in 1:d)... ) + λ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)... ) + λ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 @@ -353,7 +383,8 @@ function trembling_foot!(Π, dΠ_x, dΠ_z1, dΠ_z2, S::Vector{Point{d}}, S_x::Ve indexes_to_be_modified = tuple(n, UnitRange.(qn_,qn_.+1)...) Π[indexes_to_be_modified...] .+= w.*rhs_Π - + + # # # Computing the associated derivatives for k=1:d if smooth==false diff --git a/test/test_G_derivative.jl b/test/test_G_derivative.jl index 3b95d6e..70feae0 100644 --- a/test/test_G_derivative.jl +++ b/test/test_G_derivative.jl @@ -6,31 +6,43 @@ using Dolo 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 From 348bfe0a947f1e9f623e68a54eb58a6071a269b4 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Thu, 10 Mar 2022 17:03:07 +0100 Subject: [PATCH 20/22] small modif --- test/test_G_derivative.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_G_derivative.jl b/test/test_G_derivative.jl index 70feae0..953adcb 100644 --- a/test/test_G_derivative.jl +++ b/test/test_G_derivative.jl @@ -1,7 +1,7 @@ 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") From e471fd5a6c1160e579693ae4204522c74afd77e0 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 21 Mar 2022 14:54:48 +0100 Subject: [PATCH 21/22] updating project.toml with plotlyjs --- Project.toml | 1 + examples/models/zero_to_aiyagari_iid.yaml | 52 ++++++++++++ experiments/dev_ergodic.jl | 97 ++++++++++++++++++++--- 3 files changed, 140 insertions(+), 10 deletions(-) create mode 100644 examples/models/zero_to_aiyagari_iid.yaml 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/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 index 7ac2207..396caeb 100644 --- a/experiments/dev_ergodic.jl +++ b/experiments/dev_ergodic.jl @@ -170,9 +170,12 @@ fig ## Aiyagari model -model = yaml_import("examples/models/consumption_savings_iid.yaml") +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)] @@ -181,17 +184,91 @@ w = [nodes[i][1] for i in 1:length(nodes)] μ_smoothing = Dolo.ergodic_distribution(model, sol; smooth=true) -n = Int(length(μ_smoothing)) +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 + ) -df = DataFrame(w=vcat(w,w), μ=vcat(μ_smoothing[1:n], μ_no_smoothing[1:n]), smooth=vcat(["smoothing" for i in 1:n], ["no smoothing" for i in 1:n])) +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) + ] +) -PlotlyJS.plot( - PlotlyJS.scatter(df, x=:w, y=:μ, group=:smooth), - Layout( - xaxis_title="w", - yaxis_title="μ" - ) - ) +fig From e513e646d7867aac14b81f9376ad4e517d0e2c14 Mon Sep 17 00:00:00 2001 From: Gabrielle Queran Date: Mon, 21 Mar 2022 15:44:06 +0100 Subject: [PATCH 22/22] trying to solve the error --- src/algos/ergodic.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/algos/ergodic.jl b/src/algos/ergodic.jl index 9ab9d21..347ae6c 100644 --- a/src/algos/ergodic.jl +++ b/src/algos/ergodic.jl @@ -337,15 +337,16 @@ of the transition function wrt x, z1 and z2. * `S_x` : Derivative of S wrt x * `S_z1` : Derivative of S wrt z1 * `S_z2` : Derivative of S wrt z2 -* `w::Float64` : +* `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 +* `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