Error in asset pricing code

[jstac]

From the rules of vectorization it seems to me that we should have a reshaping of g in the definition of the operator A, as shown below. And yet maybe it's not necessary. @HumphreyYang or @shlff , could you please take a look?

If both work I prefer the version below, since it's more explicit.

Thanks to Pilar, Conor and Joaquín for bringing my attention to this.

def A(g, sv_model, shapes):
    # Set up
    P, hc_grid, Q, hd_grid, R, z_grid, β, γ, bar_σ, μ_c, μ_d = sv_model
    I, J, K = shapes
    # Reshape and broadcast over (i, j, k, i', j', k')
    hc = np.reshape(hc_grid,     (I, 1, 1, 1,  1,  1))
    hd = np.reshape(hd_grid,     (1, J, 1, 1,  1,  1))
    z = np.reshape(z_grid,       (1, 1, K, 1,  1,  1))
    P = np.reshape(P,            (I, 1, 1, I,  1,  1))
    Q = np.reshape(Q,            (1, J, 1, 1,  J,  1))
    R = np.reshape(R,            (1, 1, K, 1,  1,  K))
    g = np.reshape(g,            (1, 1, 1, I,  J,  K))
    a = μ_d - γ * μ_c
    b = bar_σ**2 * (jnp.exp(2 * hd) + γ**2 * jnp.exp(2 * hc)) / 2
    κ = jnp.exp(a + (1 - γ) * z + b)
    A = β * κ * P * Q * R
    Ag = jnp.sum(A * g, axis=(3, 4, 5))
    return Ag

[jstac]

Thanks @sgcbayli024 and @pilarbarros

[shlff]

Thanks guys. This is an interesting observation.

According to my study, the original code is correct since the arrays in the product A * gfollow the broadcasting rule, see here or here.

Here is a simple illustration: https://gist.github.com/shlff/c09def39413dbabf42dcf4fe0d41ddfd#file-investigation-ipynb

From my understanding the difference is that

• with g = np.reshape(g, (1, 1, 1, I, J, K)), it is easier to read, debug and extend,
• without it, running the code requires a bit less memory.