29/05/2023 (worked with SageMath 9.7)

I have tried to partially solve your problem. When i say partially it's not really correct since I have solved it in finite time (hopping I have make no mystakes). I need some help to go to infinite horizon. In what concerns expectations, I have not open my Sargent, Lucas.... since at least 15 years so it needs reflexion and certainly the definition of the properties of the expectation operator since if you know the distribution there is no great difficulties. But to construct the attached notebook, I use at least 3 days of work. It was interesting since I encounter some unexpected difficulties.

Hope I have not make mistakes in the resolution but I think the code is full functional or at least shows "a" path you must follow to solve such a problem. I will try to resolve it in continuous time as soon as I will have some time.

The attached code is in french. Hope you can understand it. Else you can translate some sentences. Normaly there is no difficulties in inserting it in a Jupyter notebook. I cannot upload her a notebook (I do not know why here the two first latex centered equation are not computed by the site, but copied in a notebook they are correctly displayed.

c1###Put this cell in a markup cell

soit à résoudre le programme suivant :

\[ \max_{{c_t}} \left{\left.\left(\sum_{t=0}^\infty\beta^t \log(C_t)\right) + B(x_T)\right| K_{t+1} = Y_t - C_t +r_t K_t, K_0 = \overline{K},\text{conditions}\right} \]

On sait que pour résoudre un tel programme, il faut définir un hamiltonien et des fonctions adjointes. Commençons par résoudre le programme en temps fini :

\[ \max_{{c_t}} \left{\left.\left(\sum_{t=0}^{T-1}\beta^t \log(C_t)\right)+B(x_T)\right| K_{t+1} = Y_t - C_t +r_t K_t,K_0 = \overline{K},\text{conditions}\right} \]

On sait que pour résoudre un tel programme, il faut définir un hamiltonien et des fonctions adjointes :

\[ \overline{H}_t = \beta^t f(x_t,u_t,t) + \mu_{t+1} g(x_t,u_t,t) \,\,\,\text{pour}\,\,\,t=0,1,\ldots, T-1 \]

avec ici $f(x_t,u_t,t)=\log(C_t)$, $x_t= K_T$, $u_t = C_t$, $K_{t+1}=g(x_t,u_t,t)=Y_t - C_t +r_t K_t$, $Y_t$ et $r_t$ étant des fonctions données et $\overline{K}$ une constante et $B_T= B(x_T)$ une fonction d'héritage (qui détermine l'état du système en fin d'horizon -- $T$ est ce que l'on appelle l'horizon).

On sait qu'il est possible de réécrire l'hamiltonien en valeur courantes en posant $\mu(t) = \beta^t \lambda(t)$.

\[ \overline{H}_t = \beta^tf(x_t,u_t,t) + \beta^{t+1}\lambda_{t+1} g(x_t,u_t,t) \,\,\,\text{pour}\,\,\,t=0,1,\ldots, T-1 \]

ce qui donne

\[ \beta^{-t}\overline{H}_t= H_t = f(x_t,u_t,t) + \beta\lambda_{t+1} g(x_t,u_t,t) \,\,\,\text{pour}\,\,\,t=0,1,\ldots, T-1 \]

Les conditions nécessaires d'optimalité affirment que l'Hamiltonien $H_t$ est maximisé par rapport à la variable de commande à chaque instant, soit :

\[ \frac{\partial H_t}{\partial u_t}=0 \,\,\,\text{soit ici}\,\,\, \frac{\partial H_t}{\partial C_t}=0 \]

À ces conditions d'optimalité, on doit rajouter les conditions adjointes :

\[ \lambda_t = \frac{\partial H_t}{\partial x_t}\,\,\,\text{soit ici}\,\,\, \lambda_t = \frac{\partial H_t}{\partial K_t}\,\,\,\text{pour}\,\,\, t= 0, 1,\ldots, T-1 \]

On doit alors rajouter la condition de transversalité :

\[ \lambda_T = B^\prime(x_T) \]

Les conditions signifient qu'il peut aussi y avoir des contraintes sur les commandes $u_t$ ou sur l'état $x_t$ (mais pour l'instant nous ne les prendrons pas en considération). Dans un premier temps, nous considèrerons que $B(x_T) \equiv 0$.

Commençons par définir les fonction $U$ et $g$ et les variables. Comme <it>SageMath</it> n'accepte pas les indices de la forme $t+1$, nous aurons $t\rightarrow s$ et $t+1 \rightarrow t$.

#

c2 ###code cell

β, s, λ_t, Y_s, C_s, K_s, K_t, r_s, s = var('β s λ_t Y_s C_s K_s K_t r_s s', domain='real')
U=function('U')(x)
U(x) = log(x)
H(C_s,K_s)= log(C_s)+β*λ_t*(Y_s- C_s + r_s*K_s)
show(LatexExpr(r'\text{Le hamiltonien est définit comme la fonction suivante : }'),H)
sol = solve(diff(H(C_s,K_s),C_s)==0, C_s)
CC=sol[0].rhs()
show(LatexExpr(r'''\text{A l'optimum, le long du sentier de croissance optimale : } C_s ='''),CC)

#

c3###Markup cell

On doit encore évaluer :

#

c4###code cell

A= diff(H(C_s,K_s),K_s)
show(LatexExpr(r'\frac{\partial H_t}{\partial K_t} = '),A)

#

c4###Markup cell

Après substitution, on doit donc résoudre les deux équations de récurrence :

\[ \begin{array}{cc} K_{t+1} = Y_t - (1/\lambda_{t+1}) +r_t K_t, & K_0 = \overline{K}\ \lambda_{t} = r_t \lambda_{t+1} & \lambda_T = B^\prime(x_T) \end{array}
\]

La seconde équation est une équation de récurrence linéaire.

\[ \lambda_{t+1} = \frac{1}{r_t} \lambda_{t} \]

avec une condition terminale $\lambda_T$. Nous supposerons, qu'il n'y a pas de motif d'héritage. Donc $\lambda_T=0$. Et nous supposerons que $r_t = r$ constante.

#

c5###Code cell

t, r, K_0, T, K_T=var('t r K_0 T K_T')
from sympy import Function, Symbol, rsolve
from sympy.abc import n
λ = Function('λ')
rs=rsolve(λ(n+1)-r*λ(n),λ(n),{λ(0):K_0})
λ(n)=rs
λ(t)
show(LatexExpr(r'λ(t) = '),λ(t))

#

c6###Markup cell

On doit donc trouver la valeur de $K_0$ telle que $\lambda(t) = K_T$

#

c7###Code cell

ksol=solve(λ(T)==K_T,K_0)
show(LatexExpr(r'K_0 = '),ksol[0].rhs())
λ(t)=λ(t).substitute(K_0=ksol[0].rhs())
show(LatexExpr(r'λ(t) = '),λ(t))

#

c8###Code cell

ksol=solve(λ(T)==K_T,K_0)
show(LatexExpr(r'K_0 = '),ksol[0].rhs())
λ(t)=λ(t).substitute(K_0=ksol[0].rhs())
show(LatexExpr(r'λ(t) = '),λ(t))
c(t) = CC.subs(λ_t==λ(t))
show(LatexExpr(r'\text{Le long du sentier de croissance optimale : }C(t) = '),c(t))

#

c9###Code cell

var('t r Y_0 T γ')
from sympy import Function, Symbol, rsolve
from sympy.abc import n
K = Function('K')
rs1=rsolve(K(n+1)-r*K(n)+ Y_0*γ^n-c(n),K(n),{K(0):Y_0})
K(n)=rs1
K(t)
show(LatexExpr(r'K(t) = '),K(t).full_simplify())

#

On more time I hope ther is no mystakes.