amckay / ec702 Goto Github PK

using Plots, LinearAlgebra
#
PolyBasis(K)     = [ones(size(K)) K K.^2] 
PolyGetCoef(K,Y) = PolyBasis(K) \ Y
#
β = 0.99; α = 0.36; z=1.0; γ = 2; δ = 1.0;
s_ss = ( (1/β -1 +δ) / (α*z) )^(1/(α-1)); 
f(s1; z=z,α=α,δ=δ)      =   z*(s1^α)     +(1-δ)*s1
fprime(s1; z=z,α=α,δ=δ) = α*z*(s1^(α-1)) +(1-δ)
#
S1 = [range(0.75*s_ss, 1.25*s_ss,length=100);]
T1(s1,c1) = f(s1) - c1
Π10 = [ f(s1) - δ*s1  for s1 in S1] 
EE1(s1,sp1,cp1; β=β,γ=γ) = (cp1) * ((β*fprime(sp1))^(-1/γ))
#
@inline function FPI_EE(S1, Π1; k_max=100, tol=0.001)
    dist = 10.0; k=1;
    @inbounds while(dist > tol && k<=k_max) 
        println(k)
        #
        b  = PolyGetCoef(S1, Π1)
        c1 = PolyBasis(S1)*b # ≈Π1
        sp1 = T1.(S1, c1)
        cp1 = PolyBasis(sp1)*b
        Πp1 = EE1.(S1, sp1, cp1)
        dist = norm(Π1 - Πp1) 
        k += 1
        Π1 = Πp1
    end 
    return Π1
end
@time Pol_EE = FPI_EE(S1, Π10; k_max=500)
plot(S1, Pol_EE, legend=:false)