# Euler Equation Iteration to solve an Optimization problem

**URL:** <https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135>\
**Category:** Finance and Economics\
**Tags:** optimization\
**Created:** [July 12, 2022, 11:49pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135 "2022-07-12T23:49:17Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 12, 2022, 11:49pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/1 "2022-07-12T23:49:18Z")

</div>

I would like to solve a simple/routine economics problem with [Euler Equation Iteration](https://alisdairmckay.com/Notes/NumericalCrashCourse/EulerIt.html).  
The following code works well:

```julia
using Interpolations, Plots, LinearAlgebra
#
PolyBasis(K) = [ones(size(K)) K K.^2] 
PolyGetCoef(K,Y) = PolyBasis(K) \ Y
#
β = 0.99; α = 0.36; z=1.0; γ = 2.0; δ = 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/γ))
#EE1(s1,sp1,cp1; β=β,γ=γ) = min((cp1) * ((β*fprime(sp1))^(-1/γ)), f(s1))
#
@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=5000)
plot(S1, Pol_EE, legend=:false)

```

Everything works well. But it stops converging when `γ<1.78`.

Does anyone know why? Or how to fix this?  
Here is the solution (policy function) for a grid of `γ`:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/6/46b2db347750bdbae06a08512d35d1c5e6008e4b.png)

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 13, 2022, 3:19pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/2 "2022-07-13T15:19:54Z")

</div>

I think your grid for K is too narrow. For γ=1.7, the solution converges when changing the line  
S1 = [range(0.75 \* s\_ss, 1.25 \* s\_ss, length=100);]  
to  
S1 = [range(0.75 \* s\_ss, 3 \* s\_ss, length=100);]

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 13, 2022, 3:28pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/3 "2022-07-13T15:28:32Z")

</div>

I’ll try it later.  
My goal is to solve the case gamma =1 so I can compare to closed form.  
Does that case still work?

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 13, 2022, 3:46pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/4 "2022-07-13T15:46:24Z")

</div>

γ=1 seems a hard one. I’ve tried some things including different starting guesses for Π10 and I keep getting negative values or cycles. If you actually know the solution you think you should get, you can try it as initial guess for Π10, or something very close, and see whether at least the fixed point is attracting. Otherwise you would have to solve it in a different way I’m afraid.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 13, 2022, 3:49pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/5 "2022-07-13T15:49:24Z")

</div>

It’s the usual closed form solution to the NGM w log utility:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/b/e/bed333cdfd51f1661240e017d982b54b7b41dcdf.png)

I tried using it as a guess yesterday without luck

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 13, 2022, 6:25pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/6 "2022-07-13T18:25:40Z")

</div>

I tried a different version w/ a different `approx` & starting w/ a guess equal to the closed form solution does converge:

```julia
using Interpolations, Plots, LinearAlgebra
@inline approx(gr, val0) = LinearInterpolation(gr, val0, extrapolation_bc = Line())
#
β = 0.99; α = 0.36; z=1.0; δ = 1.0; γ = 1.0;
f(s1; z=z,α=α,δ=δ) = z*(s1^α) +(1-δ)*s1
fprime(s1; z=z,α=α,δ=δ) = α*z*(s1^(α-1)) +(1-δ)
c_CF(s1; z=z,α=α,β=β) = (1-α*β)*z*(s1^α)
#
s_ss = ( (1/β -1 +δ) / (α*z) )^(1/(α-1)) 
c_ss = f(s_ss) -s_ss
#
S1 = [range(0.025*s_ss, 3*s_ss,length=1000);]
T1(s1,c1) = f(s1) - c1
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)
        #
        Πh1 = approx((S1), Π1)
        sp1 = T1.(S1, Π1)
        cp1 = Πh1.(sp1)
        Πp1 = EE1.(S1,sp1,cp1) # Π Improvement/Update
        #
        plot!(S1, Πp1) |> display
        #k % 50 ==0 ? plot!(S, Πp) |> display : nothing  
        dist = norm(Π1 - Πp1) 
        k += 1
        Π1 = Πp1
    end 
    return Π1
end
Π10 = [c_CF(s1) for s1 in S1] 
#
plot(S1, f, legend=:false)
plot!(S1, Π10)
plot!([s_ss], seriestype = :vline, lab="", color="grey")
plot!([c_ss], seriestype = :hline, lab="", color="grey")
@time Pol_EE = FPI_EE(S1, Π10; k_max=200)

```

However if I start slightly below/above the correct solution:  
`Π10 = [.99*c_CF(s1) for s1 in S1] `  
it blows up.

BTW, part of the reason is the for Kp \geq 0 we need c \leq f(K).  
I can change the EulerEquation function to incorporate this:  
`EE1(s1,sp1,cp1; β=β,γ=γ) = min((cp1) * ((β*fprime(sp1))^(-1/γ)), f(s1))`  
Now I don’t have complex number problems, but it still doesn’t converge…

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 13, 2022, 10:49pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/7 "2022-07-13T22:49:08Z")

</div>

Ok, so your map is clearly not a contraction. This is required to guarantee convergence. Luckily SpeedMapping can still work for expanding maps as long as it does not diverge too fast to infeasible points. So your change to the Euler equation is helpful.

```julia
function FPI_EE_map(Πp1, S1, Π1)
    Πh1 = approx((S1), Π1)
    sp1 = T1.(S1, Π1)
    cp1 = Πh1.(sp1)
    Πp1 .= EE1.(S1, sp1, cp1) # Π Improvement/Update
end

badΠ10 = [0.9c_CF(s1) for s1 in S1]
trueΠ10 = [c_CF(s1) for s1 in S1]

using SpeedMapping
res = speedmapping(badΠ10; (m!)=(Πp1, Π1) -> FPI_EE_map(Πp1, S1, Π1)).minimizer

norm(res - trueΠ10)

```

Which gives

```julia
julia> norm(res - trueΠ10)
1.1706510230617607e-6

```

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 13, 2022, 11:07pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/8 "2022-07-13T23:07:26Z")

</div>

It works (robustly for \gamma\>.75), I wish I understood why  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/8/f/8f51ac71ea31d2166286fc488d0eee0403c64c72.png)

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 13, 2022, 11:29pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/9 "2022-07-13T23:29:31Z")

</div>

PS: I’ve been interested in creating a generic/composable package for Fixed Point Iteration (FPI) methods & it looks like `SpeedMapping.jl` can be super-helpful.

Here is a generic 1-state/1-choice discrete-time/deterministic problem:  
(I can make it more generic…)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/6/36a755e1efcd43c56f773d7f20d83a8422e45758.png)

Here is my interpretation of the Algorithms in Coleman, Lyon, Maliar, Maliar:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/8/68743adf89e1540058d4129faa9e63964f8be9d8.png)

Here are the functional forms that need to be plugged in:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/c/2c6ebe4deb9a1916b2dce727e62307896a9707b2.png)

Are you interested in also working on creating generic FPI solvers for dynamic optimization problems in Julia?  
The world would be a much better place…

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 14, 2022, 12:55pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/10 "2022-07-14T12:55:11Z")

</div>

I’m happy it works. It does for many cases, although I can’t guarantee it always will. Some arbitrary choices must be made as to what the right direction is (thinking about it, I could give more freedom to users on this. I will think some more).

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 14, 2022, 2:12pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/11 "2022-07-14T14:12:22Z")

</div>

Initially I thought SpeedMapping.jl only speeds up function iteration. It turns out it also increases the probability the iteration will converge (based on my usage).

This is important for economics bc there is a large class of Fixed Point Iteration methods for solving dynamic optimization problems.

A few of these FPI methods (like VFI & extensions) iterate on a contraction, so they are guaranteed to converge  
1: from any guess  
2: to the unique fixed point

Many other FPI methods in general are not guaranteed to converge, from any guess, to a unique fixed point. These other methods (Euler iteration etc) don’t always work, but when they do, they often outperform VFI by a lot.  
Your package increases the probability these other FPI methods will converge which the world needs to know.

To clarify: your package doesn’t just do FP acceleration, it increases the probability of convergence.

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 14, 2022, 11:34pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/12 "2022-07-14T23:34:21Z")

</div>

I agree. This could be very helpful for many domains (for instance, I’ve had a recent discussion with a quantum physicist who needed the package adapted to his needs).

I don’t think much exists in Julia except Anderson Acceleration in NLsolve.jl and N-GMRES and O-ACCEL in Optim.jl.

While developing SpeedMapping, I tried comparing its performance to the best competition I could find, so I already wrote in Julia:

Zhou, H., Alexander, D. and Lange, K. (2011). A quasi-Newton acceleration for high-  
dimensional optimization algorithms. Statistics and computing 21: 261-273.

and

Henderson, N. and Varadhan, R. (2019). Damped Anderson acceleration with restarts  
and monotonicity control for accelerating EM and EM-like algorithms. Journal of  
Computational and Graphical Statistics 28: 834-846.

Recently, I’ve found many others. The interest in Anderson Acceleration has been growing steadily in recent years (to non-smooth maps, to domain-specific problems, etc). Choosing the right ones is the interesting part.

SpeedMapping is fairly new; I presented it at JuliaCon last year. Some problems probably still need to be ironed out, but the only way to find them at this point is for people to try it and complain back to me. So I’d happily include it in a generic FPI package. Let’s talk some more.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 15, 2022, 9:58am UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/13 "2022-07-15T09:58:31Z")

</div>

> [@nicolas](#):
>
> I don’t think much exists in Julia except Anderson Acceleration in NLsolve.jl and N-GMRES and O-ACCEL in Optim.jl.

And one baked into the ODE solvers 😅

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 16, 2022, 2:55am UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/14 "2022-07-16T02:55:41Z")

</div>

Are you aware of this package:

> **[GitHub - s-baumann/FixedPointAcceleration.jl: Fixed Point Acceleration for Julia](https://github.com/s-baumann/FixedPointAcceleration.jl)**
>
> Fixed Point Acceleration for Julia. Contribute to s-baumann/FixedPointAcceleration.jl development by creating an account on GitHub.

It has 8 algorithms for fixed point acceleration.

Would be awesome if all the algorithms were together in one good well maintained package

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 17, 2022, 12:49pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/15 "2022-07-17T12:49:43Z")

</div>

Ah yes I think I saw that one too a while ago but forgot about it. I think that’s a great start and new algorithms can surely be incorporated.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 21, 2022, 6:16pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/16 "2022-07-21T18:16:39Z")

</div>

I’d like to apply your package to a problem with 2-state variables & 2 -choice variables.  
\max\_{c\_t, b\_{t+1}, i\_{t}, h\_{t+1}} \sum\_{t=0}^{\infty} \beta^t u(c\_t, h\_t)  
b\_{t+1}= Rb\_t -c\_t -i\_t -\phi(i\_{t})  
h\_{t+1} = (1-\delta)h\_{t} + i\_{t}   
Let: u(c,h) = \alpha ln(c) + \gamma ln(h) and \phi(i) = \frac{\theta}{2} (i)^2

Two state variables: b\_{t+1}, h\_{t+1} we are given b\_0, h\_0  
Two choice variables: c\_t, i\_t  
The equations that need to be solved: (4 equations in 4 variables)  
c\_{t} = (R\beta)^{-1} c\_{t+1}  
i\_{t} = \frac{\gamma}{\alpha R \theta} \frac{c\_{t+1}}{h\_{t+1}} +\frac{1-\delta}{R} i\_{t+1} - \frac{r+\delta}{R \theta}  
b\_{t+1}= Rb\_t -c\_t -i\_t -\phi(i\_{t})  
h\_{t+1} = (1-\delta)h\_{t} + i\_{t}

Note:  
h\_{t+1} \geq 0 \Rightarrow i\_{t} \geq -(1-\delta)h\_{t}  
b\_{t+1} \geq -\eta \Rightarrow c\_t +i\_t +\phi(i\_{t}) \leq Rb\_t +\eta {natural borrowing constraint}  
c\_{t} \geq 0 \Rightarrow i\_t +\phi(i\_{t}) \leq Rb\_t - b\_{t+1} \leq Rb\_t +\eta

Maximum consumption if you sell all capital & borrow to limit:  
c\_{t} = Rb\_t -b\_{t+1} -i\_t -\phi(i\_{t})  
c\_{t} = Rb\_t +\eta +(1-\delta)h\_{t}- \frac{\theta}{2} ((1-\delta)h\_{t})^2

Thus:  
 -(1-\delta)h\_{t} \leq i\_{t} \leq \frac{1}{\theta} (-1+ \sqrt{1+2\theta (Rb\_t +\eta)} )   
 0 \leq c\_{t} \leq Rb\_t +\eta +(1-\delta)h\_{t}- \frac{\theta}{2} ((1-\delta)h\_{t})^2

Now Euler Iteration code: c\_{1t}\equiv c\_{t}, c\_{2t}\equiv i\_{t}, s\_{1t}\equiv b\_{t}, s\_{2t}\equiv h\_{t}

```julia
#
r = 0.05; R = (1+r); θ = 1.00; δ = 0.10; α = 0.50; γ = 0.50; β=1/R; 
ϕ(c2;θ=θ) = .5*θ*(c2^2)
#
S1 = [(.01 : .01 : 2.5);] 
S2 = [(.01 : .01 : 2.5);] #S = [S1 S2] 
#
T1(s1,s2,c1,c2; R=R) = R*s1 -c1 -c2 -ϕ(c2)
T2(s1,s2,c1,c2; δ=δ) = (1-δ)*s2 +c2
#
η=0.0
c1_EE1(s1,s2,s1p,s2p,c1p,c2p; R=R,β=β) = c1p*((R*β)^(-1))
c1_min(s1,s2) = 0.00001 # c1 >= 0.0
c1_max(s1,s2;R=R,δ=δ,η=η) = R*s1 +η +(1-δ)*s2 - ϕ((1-δ)*s2) # c1 <= R*(s1) + z - η => s1p>=0
#
c2_EE2(s1,s2,s1p,s2p,c1p,c2p; R=R,γ=γ,α=α,θ=θ,δ=δ) = (γ/(α*R*θ))*(c1p/s2p) +((1-δ)/R)*c2p -((R-1+δ)/(R*θ))
c2_min(s1,s2) = -(1-δ)*s2
c2_max(s1,s2;R=R,θ=θ,η=η) = (1/θ)*(-1 + sqrt(1+2*θ*(R*s1+η)))
#
EE1(s1,s2,s1p,s2p,c1p,c2p) = clamp(c1_EE1(s1,s2,s1p,s2p,c1p,c2p), c1_min(s1,s2), c1_max(s1,s2)) # clamp(x, lo, hi) 
EE2(s1,s2,s1p,s2p,c1p,c2p) = clamp(c2_EE2(s1,s2,s1p,s2p,c1p,c2p), c2_min(s1,s2), c2_max(s1,s2)) # clamp(x, lo, hi) 
#
function FPI_EE_map!(Π1p, Π2p, Π1, Π2, S1, S2)
    #
    Π1h = approx((S1, S2), Π1)
    Π2h = approx((S1, S2), Π2)
    S1p = [T1(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    S2p = [T2(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    C1p = [Π1h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    C2p = [Π2h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    Π1p .= [EE1(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    Π2p .= [EE2(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
end
#
Π10 = [.1 for s1 in S1, s2 in S2] #
Π20 = [.1 for s1 in S1, s2 in S2] #
#
res = speedmapping(Π10,Π20; (m!)=(Π1p, Π2p, Π1, Π2) -> FPI_EE_map!(Π1p, Π2p, Π1, Π2, S1, S2))

```

Is it possible to do this with your package?

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 27, 2022, 12:07am UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/17 "2022-07-27T00:07:33Z")

</div>

A cleaner & simpler example of a 2-state/2-choice problem is the Uzawa-Lucas model.  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/4/f4e7e03fd81fb8a06439e4027d0b6bc521ac3669.png)

A special case (full depreciation & log utility) has a closed form solution:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/a/3a335038b454e567e68e2dfa02cf8d1f91be02a7.png)

The Euler equations if \gamma =0 (no externality)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/7/d76f3dec95ce0d78beedf98225d0b0a4af168789.png)

Now in Julia

```julia
using Interpolations, Plots, LinearAlgebra, SpeedMapping
@inline approx(gr, val0) = LinearInterpolation(gr, val0, extrapolation_bc = Line())
#
α=.33; β=.97; #σ=1;δ=0;A=1.0; 
F(K,H,u;α=α) = (K^α)*(u^(1-α))*(H^(1-α))
#
γ=0.0;
ū=((1-α)*(1-β))/(1-α+β*γ)
B= 1/(1.0 - ū) + .02
#B= 1/(1.0 - ū) - .02
#B= 1/(1.0 - ū) 
#
S1 = [(.01 : .01 : 1.0);] 
S2 = [(.01 : .01 : 1.0);] #S = [S1 S2] PolyBasis(S1, S2)
#
T1(s1,s2,c1,c2) = F(s1,s2,c2) -c1
T2(s1,s2,c1,c2; B=B) = B*(1-c2)*s2
#
c1_EE1(s1,s2,s1p,s2p,c1p,c2p; α=α,β=β) = c1p*((β*α*(F(s1p,s2p,c2p)/s1p))^(-1))
c1_min(s1,s2) = 0.00001 # c1 >= 0.0
c1_max(s1,s2) = F(s1,s2,1.0) # c1 <=
#
c2_EE2(s1,s2,s1p,s2p,c1p,c2p;B=B,α=α) = ((s1p*B/α)*(1/s2p +1/c2p -1)/((s1^α)*(s2^(1-α))))^(-1/α)
c2_min(s1,s2) = 0.0
c2_max(s1,s2) = 1.0
#
EE1(s1,s2,s1p,s2p,c1p,c2p) = clamp(c1_EE1(s1,s2,s1p,s2p,c1p,c2p), c1_min(s1,s2), c1_max(s1,s2))
EE2(s1,s2,s1p,s2p,c1p,c2p) = clamp(c2_EE2(s1,s2,s1p,s2p,c1p,c2p), c2_min(s1,s2), c2_max(s1,s2))
#
Π10 = [(1-α*β)*F(s1,s2,ū) for s1 in S1, s2 in S2] #
Π20 = [ū for s1 in S1, s2 in S2] #
#
function FPI_EE_map!(Π1p, Π2p, Π1, Π2, S1, S2)
    #
    Π1h = approx((S1, S2), Π1)
    Π2h = approx((S1, S2), Π2)
    S1p = [T1(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    S2p = [T2(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    C1p = [Π1h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    C2p = [Π2h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    Π1p .= [EE1(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
    Π2p .= [EE2(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
end
#res = speedmapping(Π10,Π20; (m!)=(Π1p, Π2p, Π1, Π2) -> FPI_EE_map!(Π1p, Π2p, Π1, Π2, S1, S2)) #, tol=1e-10
@inline function FPI_EE(S1, S2, Π1, Π2; k_max=100, tol=0.001)
    dist = 10.0; k=1;
    Π1p = similar(Π1)
    Π2p = similar(Π2)
    @inbounds while(dist > tol && k<=k_max) 
        println(k)
        #
        Π1h = approx((S1, S2), Π1)
        Π2h = approx((S1, S2), Π2)
        S1p = [T1(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        S2p = [T2(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        C1p = [Π1h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        C2p = [Π2h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        Π1p .= [EE1(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        Π2p .= [EE2(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        #
        dist = 0.5*norm(Π1 - Π1p) + 0.5*norm(Π2 - Π2p)
        k += 1
        Π1 = Π1p
        Π2 = Π2p
    end 
    return Π1, Π2 #, va
end
@time Pol1_EE, Pol2_EE = FPI_EE(S1, S2, Π10, Π20; k_max=600, tol=0.000001)

```

Strangely, plugging in the closed form solution as my guess doesn’t seem to converge…

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [July 27, 2022, 6:05am UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/18 "2022-07-27T06:05:37Z")

</div>

I am not sure I understand the details but are you solving for the roots of a specific fuctional F(X,\gamma) = 0 as function of \gamma?

---

<div class="post-metadata">

**Author:** ![nicolas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nicolas/32/26439_2.png) [@nicolas](https://discourse.julialang.org/u/nicolas)\
**Post date:** [July 27, 2022, 8:27pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/19 "2022-07-27T20:27:22Z")

</div>

> [@Albert\_Zevelev](#):
>
> `res = speedmapping(Π10,Π20; (m!)=(Π1p, Π2p, Π1, Π2) -> FPI_EE_map!(Π1p, Π2p, Π1, Π2, S1, S2))`

Hi Albert,

In SpeedMapping, m! only takes two arguments, x\_out and x\_in. Plus I have realized that it does not accept arrays of arrays. I plan to change this in the next version. For now, I’m afraid you would need to stack the two arrays on top of each other and do a bunch of copying. It’s not that bad, but not the most elegant. Here is how you could adapt your code.

```julia
Π0 = vcat(Π10, Π20)

temp1 = similar(Π10)
temp2 = similar(Π20)
temp3 = similar(Π10)
temp4 = similar(Π20)

function FPI_EE_map!(Π_out, Π_in, Π1p, Π2p, Π1, Π2, S1, S2)
    s = size(Π1p, 1)
    Π1 .= Π_in[begin:begin+s-1, :]
    Π2 .= Π_in[begin+s:end, :]
    FPI_EE_map!(Π1p, Π2p, Π1, Π2, S1, S2)
    Π_out[begin:begin+s-1, :] .= Π1p
    Π_out[begin+s:end, :] .= Π2p
end

res = speedmapping(Π0; (m!)=(Π_out, Π_in) -> FPI_EE_map!(Π_out, Π_in, temp1, temp2, temp3, temp4, S1, S2))

```

That said, your mapping needs a fixed point. If you start at the analytical solution and it diverges, that’s bad news. I tried running the code and encountered negative values for consumption, capital, etc. I tried using max(x,1e-50) for every x that couldn’t be negative, but without great success. I’m not familiar with the Uzawa-Lucas model so it’s hard for me to see if there is a problem in the equations.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [July 27, 2022, 8:38pm UTC](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135/20 "2022-07-27T20:38:39Z")

</div>

Sorry there was a small typo in my 2nd Euler equation. Now it works w/ my FPI solver:

```julia
using Interpolations, Plots, LinearAlgebra, SpeedMapping
@inline approx(gr, val0) = LinearInterpolation(gr, val0, extrapolation_bc = Line())
#
α=.33; β=.97; #σ=1;δ=0;A=1.0; 
F(K,H,u;α=α) = (K^α)*(u^(1-α))*(H^(1-α))
#
ū= (1-β)
B= 1/(1.0 - ū) - .02
#
c2_CF(s1,s2; ū=ū) = ū
c1_CF(s1,s2; α=α,β=β) = (1-α*β)*F(s1,s2,c2_CF(s1,s2))
#
S1 = [(.01 : .01 : 1.0);] 
S2 = [(.01 : .01 : 1.0);] #S = [S1 S2] PolyBasis(S1, S2)
#
T1(s1,s2,c1,c2) = F(s1,s2,c2) -c1
T2(s1,s2,c1,c2; B=B) = B*(1-c2)*s2
#
c1_EE1(s1,s2,s1p,s2p,c1p,c2p; α=α,β=β) = c1p*((β*α*(F(s1p,s2p,c2p)/s1p))^(-1))
c1_min(s1,s2) = 0.00001      
c1_max(s1,s2) = F(s1,s2,1.0) 
#c1_max(s1,s2) = F(s1,s2,ū) 
#
c2_EE2(s1,s2,s1p,s2p,c1p,c2p;B=B,α=α) = (s1/s2) * ((α*c2p*s2p)/(B*s1p))
c2_min(s1,s2) = 0.0
c2_max(s1,s2) = 1.0
#
EE1(s1,s2,s1p,s2p,c1p,c2p) = clamp(c1_EE1(s1,s2,s1p,s2p,c1p,c2p), c1_min(s1,s2), c1_max(s1,s2))
EE2(s1,s2,s1p,s2p,c1p,c2p) = clamp(c2_EE2(s1,s2,s1p,s2p,c1p,c2p), c2_min(s1,s2), c2_max(s1,s2))
#
@inline function FPI_EE(S1, S2, Π1, Π2; k_max=100, tol=0.001)
    dist = 10.0; k=1;
    Π1p = similar(Π1)
    Π2p = similar(Π2)
    @inbounds while(dist > tol && k<=k_max) 
        println(k)
        #
        Π1h = approx((S1, S2), Π1)
        Π2h = approx((S1, S2), Π2)
        S1p = [T1(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        S2p = [T2(S1[j], S2[k], Π1[j,k], Π2[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        C1p = [Π1h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        C2p = [Π2h(S1p[j,k], S2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        Π1p .= [EE1(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        Π2p .= [EE2(S1[j], S2[j], S1p[j,k], S2p[j,k], C1p[j,k], C2p[j,k]) for j in eachindex(S1), k in eachindex(S2)]
        #
        dist = 0.5*norm(Π1 - Π1p) + 0.5*norm(Π2 - Π2p)
        k += 1
        Π1 = Π1p
        Π2 = Π2p
    end 
    return Π1, Π2 #, va
end
#

# Guess Closed Form. 
Π10 = [c1_CF(s1,s2) for s1 in S1, s2 in S2] #
Π20 = [c2_CF(s1,s2) for s1 in S1, s2 in S2] #
#
@time Pol1_EE, Pol2_EE = FPI_EE(S1, S2, Π10, Π20; k_max=600, tol=0.000001)
# WORKS!!!!!!!!!!!!!!!!!!

#
Π10 = [0.95*c1_CF(s1,s2) for s1 in S1, s2 in S2] #
Π20 = [0.95*c2_CF(s1,s2) for s1 in S1, s2 in S2] #
@time Pol1_EE, Pol2_EE = FPI_EE(S1, S2, Π10, Π20; k_max=600, tol=0.000001)

```

How can I use your package, Speedmapping.jl for this problem?

[Next page](https://discourse.julialang.org/t/euler-equation-iteration-to-solve-an-optimization-problem/84135.md?page=2)