Interpolate then maximize

Question

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]]]]])

• Gareth Rees: Could you give us enough code (and sample data) that we can actually run it?
• FooBar: I did so. Notice that you need to set V0 = the copy paste of the last part before running it.

Answer 1

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() 

• Nobody moving away from SE: 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)