Para resolver o problema, primeiro montamos a equação de Bellman, cuja função valor é:
\[
V(k)=\max\limits_{c} u(c) + \beta V(F(k)-c)
\]
Ao discretizar os estados \(k \in [0.1,2]\), a variável de controle \(c\) deve ser limitada de modo que mantenha a variável de estado \(k\) dentro do conjunto escolhido. Por exemplo, se discretizarmos os estados em \(k \in \{0.1,0.2,...,2\}\), a variável de controle \(c\) teria que estar no conjunto \(C\) tal que para cada possível variação em \(k\), \(c=F(k)-k^+\), onde \(k^+\) é o valor do estado no período seguinte. Veja que \(C\), neste exemplo, pode chegar a \(20^2\) elementos (o quadrado da quantidade de estados escolhidos). Além disso, pode haver alguma dificuldade em tratar arredondamentos. Uma maneira de lidar com isso é usar a variável de estado do período seguinte \(k^+\) como variável de controle. Assim, tanto o conjunto de estados como de controles serão os mesmos. A equação de Bellman é modificada para:
\[
V(k)=\max\limits_{k^+} u(F(k)-k^+) + \beta V(k^+)
\]
Uma vez discretizado, podemos utilizar o seguinte teorema (Denardo, 1967):
Teorema. Se o conjunto de estado \(X\) for compacto e a função payoff \(\pi(x,u)\) (payoff quando o estado inicial é \(x \in X\) e o controle aplicado é \(u \in D(x)\), onde \(D(x)\) é o conjunto de controles possíveis quando o estado é \(x\)) for limitada, acima e abaixo, então o operador linear \(T:X \rightarrow X\)
\[ TV=\sup\limits_{u \in D(x)} \pi(x,u)+\beta E\left[V(x^+) \mid x,u\right] \]
i) é monótono em \(V\), ii) é uma contração uniforme com módulo \(\beta\) no espaço de funções limitadas, e iii) tem um único ponto fixo.
Este teorema é um caso do teorema do ponto fixo em espaços de Banach (ou teorema da contração uniforme) que, além de garantir a existência e unicidade da solução para \(T(v)=v\), quando \(T\) é uma contração, induz ao método de aproximações sucessivas, dado que, independente do "chute" inicial, o valor e as regras de controle ótimo irão convergir. Uma boa referência do teorema é Qual é o teorema do ponto fixo de Banach?.
O código a seguir, em Python, implementa o método "value iteration".
import numpy as np
import matplotlib.pyplot as plt
# Função payoff
def payoff(c):
return (np.log(c))
# Função F(k)
def F(k):
return(k+0.5*(2+np.sin(2*np.pi*k))*k**(0.25))
# Chute inicial para reduzir o número de iterações necessárias
def Vc(K_list, beta):
return np.array(list(np.log(0.5*(2+np.sin(2*np.pi*k))*k**(0.25))/(1-beta) for k in K_list))
# Matriz de payoffs estado-ação
def payoff_matrix(K_list, F_list):
n = len(K_list)
# Matriz estados iniciais x estados finais (problema transformado: ações = estados finais)
R = np.zeros([n, n])
# Iteração por todos os estados iniciais
for k_index in range(n):
F_k = F_list[k_index]
# Iteração por todas as ações (= estados finais)
for k1_index in range(n):
k1 = K_list[k1_index]
# k(t+1) deve ser menor que F(k(t)), pois c(t) não pode ser negativo
if k1 <= F_k:
R[k_index, k1_index] = payoff(F_k-k1) # Utilidade do consumo
else:
R[k_index, k1_index] = -np.inf
return(R)
# Função "value iteration"
def value_iteration(k_inferior, k_superior, n, beta, v=[], maxiterations=1000, theta=0.000000001):
# Monta um vetor discretizando os estados
K_list = list(round(k_inferior + (k_superior-k_inferior)*i/n,4) for i in range(n+1))
# Monta um vetor com os valores de F(k)
F_list = list(F(K_list[i]) for i in range(n+1))
# Default: política sem variação do estado como chute inicial
if len(v)==0: v=Vc(K_list,beta)
# Inicializa vetor de políticas
u=np.zeros([n+1])
# Calcula matriz de payoffs estado-ação
R=payoff_matrix(K_list, F_list)
# delta = variação do valor em relação ao inicial. Inicialmente maior que theta para entrar no while
delta=theta+0.001
# Menor valor observado para delta.
mindelta=+np.inf
# Contador de iterações
i=0
# Iterações até que os valores convirjam ou maxiterations iterações
while(delta>theta):
i+=1
delta=0
oldV=v.copy()
# Iterações por cada estado inicial
for state in range(n+1):
# Cálculo dos valores ao levar do estado inicial state para todos os estados finais
values = R[state,:]+beta*v
# Seleciona o valor e o argumento máximos
v[state]=max(values) # valor máximo
u[state]=K_list[np.argmax(values)] # política ótima
delta=max(delta,np.abs(v[state]-oldV[state]))
mindelta = min(delta,mindelta)
print("Iteração ",i,":",delta)
if i==maxiterations: break
return v,u,i,mindelta
Além da função que realiza as iterações value_iteration - o método de aproximações sucessivas de fato - foram desenvolvidas em separado as funções payoff e F, apenas por uma questão de clareza e para facilitar o uso do código com outras configurações. Outra função auxiliar é payoff_matrix, que mapeia os payoffs aos possíveis estado-ação.
Uma dica importante é a escolha do "chute" inicial. Uma das políticas possíveis é aquela em que não há variação na poupança, ou seja, \(k^+=k\). Neste caso, o consumo será \(c=F(k)-k\), constante, e o valor será \(V^c(k)=\frac{F(k)-k}{1-\beta}\), que chamaremos de "chute inteligente". Este valor provavelmente não estará muito longe do valor verdadeiro maximizado. No código, \(V^c\) é calculado pela função Vc e utilizado, por default, como o chute inicial.
Para testar, rodamos o código seguinte.
# Principal
if __name__ == '__main__':
# Fator de desconto
beta = 0.95
# Inicializa variáveis que receberão resultados
v = 2*[[]] # listas de valores maximizados
u = 2*[[]] # listas de políticas
iterations = 2*[None] # lista com o número de iterações
deltas = 2*[None] # lista com o valor da diferença atingida na convergência
# Aplica a função "value iteration" para encontrar os valores maximizados e a política ótima
## Chama a função usando um chute inicial "inteligente" (default)
[v[0],u[0],iterations[0],deltas[0]] = value_iteration(0.10, 2.00, 190, beta)
## Chama a função usando como chute inicial uma matriz de zeros
[v[1],u[1],iterations[1],deltas[1]] = value_iteration(0.10, 2.00, 190, beta, v=np.zeros([190+1]))
# Desenha o gráfico com as políticas ótimas encontradas
plt.figure(1)
plt.scatter(list(0.1 + (2-0.1)*i/190 for i in range(190+1)), u[0], c=u[0])
plt.title('Políticas ótimas')
plt.xlabel('$k$')
plt.ylabel('$k^+$')
# Desenha o gráfico com valores ótimos encontrados
plt.figure(2)
plt.scatter(list(0.1 + (2-0.1)*i/190 for i in range(190+1)), v[0], c=v[0])
plt.title('Valores ótimos')
plt.xlabel('$k$')
plt.ylabel('$v(k)$')
plt.show
# Diferença entre o número de iterações partindo de um valor inicial sugerido e de zeros
dif_iterations_qty=iterations[0]-iterations[1]
dif_iterations_perc=round(dif_iterations_qty/iterations[1],3)
# Verifica a diferença entre os valores médios (assumindo distribuição uniforme na probabilidade de ocorrência de qualquer estado inicial)
dif_values_mean=np.mean(v[0])-np.mean(v[1])
# Verifica qual a diferença absoluta nos valores ótimos
dif_values_abs=sum(np.abs(v[0]-v[1]))
# Verifica se quantas políticas são diferentes
dif_controls_qty=sum(np.abs(u[0]-u[1])!=0)
O método implementado computacionalmente foi utilizado para resolver o problema proposto, usando primeiro o "chute inteligente" (\(V^c\)) e, depois, o vetor de zeros. Os valores obtidos em ambos os casos foram praticamente idênticos. As políticas foram iguais. A grande diferença é que foram necessárias 289 iterações, usando o vetor de zeros, e apenas 5 (!) iterações usando o vetor "inteligente". Representa uma queda de mais de 98% ou 284 iterações a menos.
As políticas e os valores ótimos encontrados estão sintetizados nos gráficos seguintes.
