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]]]]])
V0 =
the copy paste of the last part before running it. \$\endgroup\$