Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Como utilizar o algoritmo "value iteration" para resolver um problema determinístico de programação dinâmica, tempo discreto, horizonte infinito e estado contínuo?

0 votos
27 visitas
perguntada Jul 12 em Economia por Carlos A T Haraguchi (21 pontos)  

Resolva o problema
\[ \max\limits_{c_t} \sum\limits_{t=0}^{\infty} \beta_t u(c_t) \]
s.a. \(k_{t+1}=F(k_t)-c_t\)

onde \(\beta=0.95, u(c)=\log c\), e \(F(k)=k+0.5(2+\sin{2\pi k})k^{0.25}\),
discretizando \(k\) e usando value function iteration para \(V(k), k \in [0.1,2]\).

Esta pergunta é parte do exercício 4, capítulo 12, p. 443 do Livro "Numerical Methods in Economics", de Kenneth L. Judd (1998).

Compartilhe

1 Resposta

+1 voto
respondida Jul 12 por Carlos A T Haraguchi (21 pontos)  
editado Jul 12 por Carlos A T Haraguchi

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.

A imagem será apresentada aqui.

A imagem será apresentada aqui.

comentou Jul 13 por Pedro Gomes (16 pontos)  
Excelente resolução.
Outra forma de resolver esse problema é formar um lagrangeano com o problema de maximização, e resolver para \( c_{t} \), \( k_{t+1} \) e \( k_{t} \) , obter as condições de primeira ordem e resolver um sistema de equações, também é necessário algumas suposições para o modelo convergir nesse caso, tal como:

\[ k_{T-1} = k_{T} \]
para t = 1,..., T
Obtendo o sistema de equações, no python, pode-se utilizar a funçãoscipy.optimize.broyden1 para resolver o problema.
comentou Jul 13 por Carlos A T Haraguchi (21 pontos)  
Pedro, obrigado pelo comentário. Agradeço também a dica da solução por Lagrange e a função do python. Vou tentar resolver pela sua sugestão e comparar os resultados.
...