10
\$\begingroup\$

I need to optimize the following function:

$$ \max_{a',\ m'}\ f(a, m, e, m', a') $$

I have approximated \$f\$ with a grid \$F\$, which has a shape \$(nA, nM, nE, nM, nA)\$. Now, I want to interpolate over the last two dimensions (the ones I need to maximize over), and then maximize.

The following code block is the one I want comments on. It looks to be too complicated and to inefficient to be the best way to do this. For example, there should be a way such that I do not need to iterate over all the given states (dimensions 0–2). I'm open for any input to improve it.

# takes grid indices (first three dimensions) idx and interpolates on V
def interpolateV(idx, V, Grid, Param):
    from scipy.interpolate import interp2d

    f = interp2d(Grid.mGrid, Grid.aGrid, V[idx])
    return f

(...)
    # the code in a different function
    v1Max = empty(s2)
    v1ArgMaxA = empty(s2)
    v1ArgMaxM = empty(s2)
    from scipy import optimize
    for idx in np.ndindex(V0[..., 0,0]):
        V0i = interpolateV(idx, V0, Grid, Param)
        x, f, d = optimize.fmin_l_bfgs_b(lambda x: -V0i(x[0], x[1]), array([1, 1]), bounds=[(Grid.aMin, Grid.aMax), (Grid.mMin, Grid.mMax)], approx_grad=True)
        v1Max[idx] = f
        v1ArgMaxA[idx], v1ArgMaxM[idx] = x

This is part of a way more complex project. The next code block is added in order to make the code reproducable. Please refrain from commenting structure and code in this block, as it is messy on purpose (in order to make it as short as possible).

I reduced state size here to the minimum that makes sense, and pasted the matrix V0 in the last code block, as it would be too much code to add creation code here. Maximizer does for now return (0,1) for all states, this should be correct — I think (do not judge the following, poor code already sad about having been conglomerated).

from scipy.stats import norm
from numpy import *
# i know this is bad, I'm working with code from several instances here - will clean up the old numpy = * later, promised!
import numpy as np

def getMarcovChain (lambbda, sigma, m, N):
    """create transition matrix and and state space for y_t = lambda*y_t-1 + u_t""" 

    " , where ut distributed by a gaussian with std sigma"
    " , N is the number of states"
    " , m = Ymax/Vary = - Ymin/Vary"

    sdY = (sigma**2/(1-lambbda**2))**(0.5)
    yMax = m*sdY
    yMin = -yMax
    w = (yMax-yMin)/(N-1) # length of each discretized state
    foo = sigma**2/(1-lambbda**2)
    s = linspace(yMin, yMax, N)
    Tran = zeros((N,N))
    for j in arange(0, N):
        for k in arange(1, N-1):
            C1 = norm.cdf(s[k]- lambbda*s[j] + w/2, scale=sigma)
            C2 = norm.cdf(s[k]- lambbda*s[j] - w/2, scale=sigma)


            Tran[j,k] = C1-C2
            #import pdb; pdb.set_trace()
        Tran[j,0] = norm.cdf(s[0]- lambbda*s[j] + w/2, scale=sigma)
        Tran[j,N-1] = 1 - norm.cdf(s[N-1]- lambbda*s[j] - w/2, scale=sigma)   


    if any(abs(Tran.sum(axis=1) -1) > exp(-10)):        
        print Tran
        print Tran.sum(axis=1)

        raise Exception ("some axis=1 does not sum up")
    return [s, Tran]


class Grids(object):
    nE = 2
    nA = 5
    nM = 5

    M = 3
    A = 7

    mMin = 0
    mMax = M
    aMin = 0
    aMax = A

    def __init__(self, Param):
        self.reset(Param);

    def reset(self, Param):
        self.mGrid = linspace(self.mMin, self.mMax, self.nM)
        self.aGrid = linspace(self.aMin, self.aMax, self.nA)

        import marcov # for transition matrix and shock grid

        [eGridLog, transitionE] = marcov.getMarcovChain(Param.rho, Param.sigma,
                                                        2, self.nE)

        self.eGrid = exp(eGridLog)
        self.transitionE = transitionE

class Parameters(object):
    # utility and parameters
    gamma = 1
    beta = 0.99
    sigma = 0.25
    rho = 0.9
    kappa = 0.5 # consumption loss adjustment cost of m
# takes grid indices (first three dimensions) idx and interpolates on V
def interpolateV(idx, V, Grid, Param):
    from scipy.interpolate import interp2d

    f = interp2d(Grid.mGrid, Grid.aGrid, V[idx])
    return f

Param = Parameters()
Grid = Grids(Param)
s2 = (Grid.nM, Grid.nA, Grid.nE)
v1Max = empty(s2)
v1ArgMaxA = empty(s2)
v1ArgMaxM = empty(s2)
from scipy import optimize
for idx in np.ndindex(V0[..., 0,0].shape):
    V0i = interpolateV(idx, V0, Grid, Param)
    x, f, d = optimize.fmin_l_bfgs_b(lambda x: -V0i(x[0], x[1]), array([1, 1]), bounds=[(Grid.aMin, Grid.aMax), (Grid.mMin, Grid.mMax)], approx_grad=True)
    v1Max[idx] = f
    v1ArgMaxA[idx], v1ArgMaxM[idx] = x

And finally, the matrix V0 with shape (5, 5, 2, 5, 5), can be set with V0 = the following before running code:

array([[[[[  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195]],

         [[  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ]]],


        [[[  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195]],

         [[  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ]]],


        [[[  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195]],

         [[  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ]]],


        [[[  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195]],

         [[  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ]]],


        [[[  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195],
          [  0.423, -10.187, -11.398, -13.734, -17.195]],

         [[  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ],
          [  2.364,   1.928,   1.455,   0.923,   0.24 ]]]],



       [[[[  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734]],

         [[  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923]]],


        [[[  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734]],

         [[  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923]]],


        [[[  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734]],

         [[  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923]]],


        [[[  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734]],

         [[  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923]]],


        [[[  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734],
          [  1.05 ,   0.423, -10.187, -11.398, -13.734]],

         [[  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923],
          [  2.775,   2.364,   1.928,   1.455,   0.923]]]],



       [[[[  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398]],

         [[  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455]]],


        [[[  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398]],

         [[  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455]]],


        [[[  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398]],

         [[  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455]]],


        [[[  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398]],

         [[  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455]]],


        [[[  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398],
          [  1.565,   1.05 ,   0.423, -10.187, -11.398]],

         [[  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455],
          [  3.166,   2.775,   2.364,   1.928,   1.455]]]],



       [[[[  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187]],

         [[  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928]]],


        [[[  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187]],

         [[  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928]]],


        [[[  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187]],

         [[  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928]]],


        [[[  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187]],

         [[  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928]]],


        [[[  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187],
          [  2.028,   1.565,   1.05 ,   0.423, -10.187]],

         [[  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928],
          [  3.542,   3.166,   2.775,   2.364,   1.928]]]],



       [[[[  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423]],

         [[  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364]]],


        [[[  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423]],

         [[  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364]]],


        [[[  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423]],

         [[  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364]]],


        [[[  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423]],

         [[  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364]]],


        [[[  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423],
          [  2.458,   2.028,   1.565,   1.05 ,   0.423]],

         [[  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364],
          [  3.905,   3.542,   3.166,   2.775,   2.364]]]]])
Improve this question
\$\endgroup\$
2
  • \$\begingroup\$ Could you give us enough code (and sample data) that we can actually run it? \$\endgroup\$
    – Gareth Rees
    Jul 21 '14 at 14:35
  • \$\begingroup\$ I did so. Notice that you need to set V0 = the copy paste of the last part before running it. \$\endgroup\$
    – FooBar
    Jul 21 '14 at 14:53
2
\$\begingroup\$

It's not useful to have class Parameters just store constants...simply place those constants at the top of the code, after the import statements.

Place all of your import statements at the beginning of the code. They can be accessed by any functions or methods below them.

I'd combine the lines 

yMax = m*sdY
yMin = -yMax

to

yMax, yMin = m*sdY, -1*(m*sdY)

which seems a bit more readable to me.

Run the script with a __main__ statement at the end. 

if __name__=='__main__':
    grid = Grids()
Improve this answer
\$\endgroup\$
1
  • \$\begingroup\$ I would disagree with your readability argument. It might actually be harder to see what is assigned where. These tupled assignments are good for avoiding temporary variables though (if the two variables need the value before the other was changed) \$\endgroup\$
    – Nobody moving away from SE
    Jul 24 '14 at 18:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.