anitamette.jl

β = 0.9995
Wₕ = 100
W̄ₕ = Wₕ/(1-β)
Wₗ = 50
W̄ₗ = Wₗ/(1-β)
pₕ = 0.35
pₗ = 0.65
Bₖ = 0.6*Wₗ
U(x) = sqrt(x)
Uₕ  = U(Bₖ-1)
Uₗ  = U(Bₖ-1)
Uₗₕ = U(Bₖ-2)
Uᵤ  = U(Bₖ)
vₕ(Vᵤ)   = Uₕ+β*(pₕ*W̄ₕ+(1-pₕ)*Vᵤ)
vₕ′(Vᵤ)  = β*(1-pₕ)
vₗ(Vᵤ)   = Uₗ+β*(pₗ*W̄ₕ+(1-pₗ)*Vᵤ)
vₗ′(Vᵤ)  = β*(1-pₗ)
vₗₕ(Vᵤ)  = Uₗₕ+β*(pₕ*W̄ₕ+pₗ*(1-pₕ)*W̄ₗ+(1-pₗ)*(1-pₕ)*Vᵤ)
vₗₕ′(Vᵤ) = β*(1-pₗ)*(1-pₕ)
vᵤ(Vᵤ)   = Uᵤ+β*Vᵤ
vᵤ′(Vᵤ)  = β

function Γ(Vᵤ)
    vₕVᵤ, vₗVᵤ, vₗₕVᵤ, vᵤVᵤ = vₕ(Vᵤ), vₗ(Vᵤ), vₗₕ(Vᵤ),  vᵤ(Vᵤ)

    K = max(vₕVᵤ, vₗVᵤ, vₗₕVᵤ, vᵤVᵤ)

    return K + log(exp(vₕ(Vᵤ)-K) + exp(vₗ(Vᵤ)-K) +
                   exp(vₗₕ(Vᵤ)-K) + exp(vᵤ(Vᵤ)-K))
end

tol = []
Vᵤ = Γ(0.0)
for i ∈ 1:10^6
    Vᵤold = Vᵤ
    Vᵤ = Γ(Vᵤold)
    push!(tol, abs(Vᵤold-Vᵤ))
    if tol[end] < 1e-12
        break
    end
end

# gem da den bliver overskrevet nedenfor
Vᵤ_vfi = Vᵤ

using Plots
plot(1:length(tol), log.(tol))

function Γ′(Vᵤ)
    vₕVᵤ, vₗVᵤ, vₗₕVᵤ, vᵤVᵤ = vₕ(Vᵤ), vₗ(Vᵤ), vₗₕ(Vᵤ), vᵤ(Vᵤ)

    K = max(vₕVᵤ, vₗVᵤ, vₗₕVᵤ, vᵤVᵤ)

    Σvⱼ = exp(vₕVᵤ-K) + exp(vₗVᵤ-K) + exp(vₗₕVᵤ-K) + exp(vᵤVᵤ-K)
    Pₗ  = exp(vₗVᵤ-K)/Σvⱼ
    Pₕ  = exp(vₕVᵤ-K)/Σvⱼ
    Pₗₕ = exp(vₗₕVᵤ-K)/Σvⱼ
    Pᵤ  = exp(vᵤVᵤ-K)/Σvⱼ
    Γ′Vᵤ = Pₗ*vₗ′(Vᵤ)+Pₕ*vₕ′(Vᵤ)+Pₗₕ*vₗₕ′(Vᵤ)+Pᵤ*vᵤ′(Vᵤ)
end

tol_newton = []
Vᵤ = Γ(0.0)
for i ∈ 1:10^6
    Vᵤold = Vᵤ
    Vᵤ = Vᵤ+(Γ(Vᵤ)-Vᵤ)/(1-Γ′(Vᵤ))
    push!(tol_newton, abs(Vᵤold-Vᵤ))
    if tol_newton[end] < 1e-12
        break
    end
end
tol_newton
plot(1:length(tol_newton), tol_newton)