3D plot of invariant distribution of assets and consumption

1

Author: Fridmund K. Shunsaku

I refer to the http://quant-econ.net/py/ifp.html and its corresponding solution notebook.
I have managed to obtain the invariant asset distribution and the invariant consumption distribution, but I am wondering if it is possible to plot a 3D surface of x=asset, y=consumption, z=density.

I wasn’t sure how to even start about going that, so I gave myself an easier task, that is to plot the optimal consumption plan obtained straight from the coleman, Kc.
Specifically, plot a 3D surface of x=asset (a), y=shock (z), z=consumption ©

I have seen many different ways of plotting a 3D surface on python (including the one available on quantecons), so I picked one and tried it.

def f(X, Y):
    # 0 is the borrowing constraint, 4 is gridmax, 50 is gridsize
    asset_grid = np.linspace(0, 4, 50)
    
    for i in range(len(X[0,:])):
        for j in range(len(Y[:,0])):
            cf = lambda a, i_z: interp(a, asset_grid, c[:, i_z])
            z = cf(X,Y)
    return z

fig = plt.figure(figsize=(7, 7))
ax = fig.gca(projection='3d')

# This is range of x (set to 4 since gridmax=4)
x = np.arange(0, 4)
# This is range of y (set to 2 since there were two shocks)
y = np.arange(0, 2)

# I have read around and this seems to form the 2D base
X,Y = np.meshgrid(x, y)
Z = f(X,Y)

surf = ax.plot_surface(X, Y, Z, linewidth=0, antialiased=False)
ax.set_zlim('tight')

# Not sure what this part does, but I suppose it is for formatting the axis
ax.zaxis.set_major_locator(plt.LinearLocator(10))
ax.zaxis.set_major_formatter(plt.FormatStrFormatter('%.02f'))

plt.show()

However, I met with the following error:

ValueError: object too deep for desired array

Is anyone able to see what is the error here? Also I am keen on the original project/task (and if it is even possible)!

2

Author: Fridmund K. Shunsaku

If I may just add; the original task should be possible! Since “any deterministic function of a stationary process is stationary” got this from Prof. John and I found an article that managed a 3D plot of invariant asset and income, http://rsta.royalsocietypublishing.org/content/372/2028/20130397#disp-formula-15.

3

Hi Fridmund,

Could you please post all your code including import statements, etc.? If your code is getting long you can use a gist like this one.

While I couldn’t replicate your error message these are guesses from viewing your code:

  • The function f has to be vectorized for meshgrid to work. try f = np.vectorize(f) after the function definition.

  • for the x axis np.arrange(0, 4) returns (0, 1, 2, 3). Is that what you are after?

  • show quoted text -
4

Author: Fridmund K. Shunsaku

Hi John,

I am grateful for the reply. I posted my code here:

  • After adding np.vectorize I am still experiencing the same “Object too deep” error; from this error, I can’t really see the nature of the output that the vect. func will give and hence its intuition.

  • The x-axis was tricky for me since I wanted the actual asset values but in the context of extracting c from the f(X,Y), I needed it to be an index as well to access every asset point.
    Perhaps I could change it to np.arrange(0,50) since there are 50 points in the gridsize, but that would not be what I have described the x-axis to be.

Fridmund

5

Hi Fridmund,

Just to be totally clear, your aim is to produce a 3D plot of the stationary distribution of the pair (a_t, c_t) under the optimal policy?

6

Hi Fridmund,

Just to follow up, if you want to plot a density on a continuous state space from simulated observations then you probably want to use a nonparametric kernel density estimator.

See, for example,


http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html

7

Author: Fridmund K. Shunsaku

Hi John,

Thanks. That is an amazing advice. I have since manage to do my original task of a stationary distri of the pair (a_t, c_t) under the optimal policy.

Fridmund

figure_1.png844×608 32.4 KB

8

Hi Fridmund,

I know this is not exactly the question you asked, but since your distribution is defined over a 2d space, you can get good visual feedback by looking at a heatmap with level curves. There is a library that does this fairly well among other things: seaborn (with some bits of kde inside).

Best,

Pablo

9

HI Fridmund,

Well done, nice work.

It looks to me from your plot that the bandwidth could be increased a little bit — this will make the density estimate less spiky — if you are interested please google about choice of bandwidth for nonparametric density estimators and experiment with different values.

I second Pablo’s comment.

Regards, John.