FixedPointAcceleration.jl/ConsumptionSmoothingProblem.jl at master · s-baumann/FixedPointAcceleration.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
using SchumakerSpline
using HCubature
using Distributions
using Random
using Optim
using FixedPointAcceleration
delta = 0.2
beta = 0.95
periodic_income = 1.0
shock_var = 1.0
shock_process = LogNormal(0.0, shock_var)
BudgetStateSpace = vcat( collect(0:0.015:periodic_income), collect(1.05:0.05:(3*periodic_income)))
InitialGuess = sqrt.(BudgetStateSpace)

function ValueGivenShock(Budget::Float64, epsilon::Float64, NextValueFunction::Schumaker)
    opt = optimize(x ->  -1.0*(epsilon*(x^delta) + beta*evaluate(NextValueFunction, Budget - x + periodic_income)), 0.0, Budget)
    return -1.0 * opt.minimum
end

function ExpectedUtility(Budget::Float64, NextValueFunction::Schumaker)
    if Budget > 0.00001
        integ = hcubature(epsilon-> ValueGivenShock(Budget, epsilon[1], NextValueFunction)* pdf(shock_process, epsilon[1]), [quantile(shock_process,0.0001)], [quantile(shock_process, 0.9999)])
        return integ[1]
    else
        return beta * evaluate(NextValueFunction, periodic_income)
    end
end

function OneIterateBudgetValues(BudgetValues::Array{Float64,1})
    NextValueFunction = Schumaker(BudgetStateSpace, BudgetValues)
    new_budget_values = zeros(length(BudgetStateSpace))
    for i in 1:length(BudgetStateSpace)
        new_budget_values[i] = ExpectedUtility(BudgetStateSpace[i], NextValueFunction)
    end
    return new_budget_values
end

function new_point(level::Float64, stepp::Float64, max_grad::Float64, grad_frac::Float64)
    expp = exp(grad_frac)
    frac = expp / (expp+1)
    if isinf(max_grad)
        return level + expp
    else
        return level + stepp * max_grad * frac
    end
end

function invert_new_point(level::Float64, previous_level::Float64, previous_stepp::Float64, previous_max_grad::Float64)
    rise = level - previous_level
    if isinf(previous_max_grad)
        return log(rise)
    else
        R = rise/(previous_stepp * previous_max_grad)
        return log(R/(1-R))
    end
end

function ShapeToBud(ShapeVec)
    lenlen = length(ShapeVec)
    BudgetValues = Array{Float64,1}(undef, lenlen)
    if lenlen > 0 BudgetValues[1] = ShapeVec[1] end
    if lenlen > 1 BudgetValues[2] = new_point(BudgetValues[1], 1.0,Inf,ShapeVec[2]) end
    if lenlen > 2
        for i in 3:lenlen
            stepp = BudgetStateSpace[i] - BudgetStateSpace[i-1]
            previous_grad = (BudgetValues[i-1]-BudgetValues[i-2])/(BudgetStateSpace[i-1] - BudgetStateSpace[i-2])
            BudgetValues[i] = new_point(BudgetValues[i-1], stepp, previous_grad, ShapeVec[i])
        end
    end
    return BudgetValues
end

function BudToShape(BudgetValues)
    lenlen = length(BudgetValues)
    new_shape_vec = Array{Float64,1}(undef, lenlen)
    if lenlen > 0 new_shape_vec[1] = BudgetValues[1] end
    if lenlen > 1 new_shape_vec[2] = invert_new_point(BudgetValues[2],BudgetValues[1], 1.0,Inf) end
    if lenlen > 2
        for i in 3:lenlen
            stepp = BudgetStateSpace[i] - BudgetStateSpace[i-1]
            previous_grad = (BudgetValues[i-1]-BudgetValues[i-2])/(BudgetStateSpace[i-1] - BudgetStateSpace[i-2])
            new_shape_vec[i] = invert_new_point(BudgetValues[i],BudgetValues[i-1], stepp, previous_grad)
        end
    end
    return new_shape_vec
end

function Reparameterised_FP_Vector(ShapeVec)
    BudgetValues = ShapeToBud(ShapeVec)
    new_Budget_Values = OneIterateBudgetValues(BudgetValues)
    new_shape_vec = BudToShape(new_Budget_Values)
    return new_shape_vec
end
function shapeconvergence(Inputs, Outputs)
    return maximum(abs.(ShapeToBud(Inputs) .- ShapeToBud(Outputs)))
end

# Testing that the conversions work properly.
fp = fixed_point(OneIterateBudgetValues, InitialGuess; PrintReports = true, MaxIter = Integer(1))
shape_guess = BudToShape(InitialGuess)
fp_reparam = fixed_point(Reparameterised_FP_Vector, shape_guess; PrintReports = true, MaxIter = Integer(1), ConvergenceMetric = shapeconvergence)
iterated = ShapeToBud(fp_reparam.Outputs_[:,1])
sum(abs.(fp.Outputs_[:,1] .- iterated) .> 1e-10)  == 0

# fixed point acceleration with the reparameterised version.
fp = fixed_point(OneIterateBudgetValues, InitialGuess; PrintReports = true, ConvergenceMetricThreshold = 1e-06)
fpfp = fp.FixedPoint_
fp_ander = fixed_point(Reparameterised_FP_Vector, shape_guess; PrintReports = true, ReportingSigFig = Integer(10), ConvergenceMetricThreshold = 1e-06, ConvergenceMetric = shapeconvergence)
iterated = ShapeToBud(fp_ander.FixedPoint_)
sum(abs.(fpfp .- iterated) .> 1e-4)  == 0
fp_aitken = fixed_point(Reparameterised_FP_Vector, shape_guess; PrintReports = true, ReportingSigFig = Integer(10), ConvergenceMetricThreshold = 1e-06, Algorithm = Aitken, ConvergenceMetric = shapeconvergence)
iterated = ShapeToBud(fp_aitken.FixedPoint_)
sum(abs.(fpfp .- iterated) .> 1e-4)  == 0
fp_newton = fixed_point(Reparameterised_FP_Vector, shape_guess; PrintReports = true, ReportingSigFig = Integer(10), ConvergenceMetricThreshold = 1e-06, Algorithm = FixedPointAcceleration.Newton, ConvergenceMetric = shapeconvergence)
iterated = ShapeToBud(fp_newton.FixedPoint_)
sum(abs.(fpfp .- iterated) .> 1e-4)  == 0
fp_sea = fixed_point(Reparameterised_FP_Vector, shape_guess; PrintReports = true, ReportingSigFig = Integer(10), ConvergenceMetricThreshold = 1e-06, Algorithm = SEA,  ConvergenceMetric = shapeconvergence)
iterated = ShapeToBud(fp_sea.FixedPoint_)
sum(abs.(fpfp .- iterated) .> 1e-4)  == 0