In [5]:

    

# pep8修正済み
#  pulp
# そのほか時間測る系もまとめる
# sigmaをcompute_greedyで求める

% matplotlib inline
import matplotlib.pyplot as plt
import pulp
import quantecon as qe
import numpy as np
import time
from gurobipy import *

class SimpleOG(object):

    def __init__(self, B=10, M=5, alpha=0.5, beta=0.9):

        self.B, self.M, self.alpha, self.beta  = B, M, alpha, beta
        self.n = B + M + 1
        self.m = M + 1

        self.R = np.empty((self.n, self.m))
        self.Q = np.zeros((self.n, self.m, self.n))

        self.populate_Q()
        self.populate_R()

    def u(self, c):
        return c**self.alpha

    def populate_R(self):

        for s in range(self.n):
            for a in range(self.m):
                self.R[s, a] = self.u(s - a) if a <= s else -np.inf

    def populate_Q(self):

        for a in range(self.m):
            self.Q[:, a, a:(a + self.B + 1)] = 1.0 / (self.B + 1)


# solve by linear programming
def LP_approach(setting):
    num_state = range(1, setting.num_states + 1)
    num_action = range(1, setting.num_actions + 1)
    LP = pulp.LpProblem('LP_Approach', pulp.LpMinimize)

    value_functions = pulp.LpVariable.dicts('value', num_state, 0, 50000, 'Continuous')
    
    # set constraints
    dimention = setting.num_states * setting.num_actions
    c_s = [None] * dimention
    count = 0
    for i in num_state:
        for j in num_action:
            c = value_functions[i] - \
                   setting.beta * (pulp.lpSum(setting.Q[i-1, j-1, k-1] * value_functions[k] for k in num_state))\
                   >= setting.R[i-1,j-1]
            LP += c
            c_s[count] = c
            count += 1
    
    # set objective
    LP += pulp.lpSum(value_functions[i] for i in num_state)
    
    # solve
    LP.solve(pulp.GUROBI(msg = False))
    
    # values
    v = np.empty(setting.num_states)
    for i, value in enumerate(value_functions.values()):
        v[i] = value.value()
    
    # sigma
    sigma = setting.compute_greedy(v)
    
    # result
    res = qe.markov.ddp.DPSolveResult(v=v, sigma = sigma,
                            mc=setting.controlled_mc(sigma),
                            method = 'Exact Linear Programming',
                            )
    return res
    
# linear programming time
def elapse_LP(setting):
    start = time.time()
    LP_approach(setting)
    return time.time() - start

# value function time
def elapse_value(setting):
    start = time.time()
    setting.solve(method='value_iteration')
    elapsed_time = time.time() - start
    return elapsed_time

# plot time until num_state = n
def plot_graph(n, discount):
    a = [None] * n
    b = [None] * n
    for i in range(n):
        h = SimpleOG(M = i, beta = discount)
        d = qe.markov.DiscreteDP(h.R, h.Q, h.beta)
        a[i] = elapse_LP(d)
        b[i] = elapse_value(d)

    plt.plot(a)
    plt.plot(b)
    
# time changing beta LP
def plot_beta_LP(shock, state, betas):
    n = [None] * len(betas)
    for i, value in enumerate(betas):
        model = SimpleOG(B = shock, M = state, beta = value)
        ddp = qe.markov.DiscreteDP(model.R, model.Q, model.beta)
        n[i] = elapse_LP(ddp)
    plt.plot(n)
    
#time changing beta value function
def plot_beta_value(shock, state, betas):
    n = [None] * len(betas)
    for i, value in enumerate(betas):
        model = SimpleOG(B = shock, M = state, beta = value)
        ddp = qe.markov.DiscreteDP(model.R, model.Q, model.beta)
        n[i] = elapse_value(ddp)
    plt.plot(n)


    

In [2]:

    

# お試し
h = SimpleOG()
d = qe.markov.DiscreteDP(h.R, h.Q, h.beta)
res = LP_approach(d)
res


    

    
Out[2]:





  sigma: array([0, 0, 0, 0, 1, 1, 1, 2, 2, 3, 3, 4, 5, 5, 5, 5])
 method: 'Exact Linear Programming'
     mc: Markov chain with transition matrix 
P = 
[[ 0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.          0.
   0.          0.          0.        ]
 [ 0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.          0.
   0.          0.          0.        ]
 [ 0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.          0.
   0.          0.          0.        ]
 [ 0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.          0.
   0.          0.          0.        ]
 [ 0.          0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.          0.          0.          0.        ]
 [ 0.          0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.          0.          0.          0.        ]
 [ 0.          0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.          0.          0.          0.        ]
 [ 0.          0.          0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.          0.          0.        ]
 [ 0.          0.          0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.          0.          0.        ]
 [ 0.          0.          0.          0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.          0.        ]
 [ 0.          0.          0.          0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.          0.        ]
 [ 0.          0.          0.          0.          0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.        ]
 [ 0.          0.          0.          0.          0.          0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909]
 [ 0.          0.          0.          0.          0.          0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909]
 [ 0.          0.          0.          0.          0.          0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909]
 [ 0.          0.          0.          0.          0.          0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909  0.09090909  0.09090909
   0.09090909  0.09090909  0.09090909  0.09090909]]
      v: array([ 19.01740222,  20.01740222,  20.43161578,  20.74945302,
        21.04078099,  21.30873018,  21.54479816,  21.76928181,
        21.98270358,  22.18824323,  22.3845048 ,  22.57807736,
        22.76109127,  22.94376708,  23.11533996,  23.27761762])

In [6]:

    

plot_graph(50,0.99)


    

In [4]:

    

h = SimpleOG(M = 50, beta = 0.99)
d = qe.markov.DiscreteDP(h.R, h.Q, h.beta)
elapse_LP(d)


    

    
Out[4]:





8.353132009506226

In [ ]: