calc_expectation returns additional dimension

[alanlujan91]

HARK/HARK/ConsumptionSaving/ConsPortfolioModel.py

Lines 657 to 663 in 1d37908

	# Calculate intermediate marginal value of bank balances by taking expectations over income shocks
	dvdb_intermed = calc_expectation(
	self.IncShkDstn, dvdb_dist, self.bNrmNext, self.ShareNext
	)
	# calc_expectation returns one additional "empty" dimension, remove it
	# this line can be deleted when calc_expectation is fixed
	dvdb_intermed = dvdb_intermed[:, :, 0]

In the above code, calc_expectation returns an additional dimension (dvdb_intermed is returned with shape [b_N, s_N, 1]). This might have to do with passing 2 arguments to the function since this is not an issue that I'm aware of when only 1 argument is passed.

[sbenthall]

This bit of code is handling the dimensionality of the calc_expectations output.

HARK/HARK/distribution.py

Lines 1341 to 1345 in 84ee03e

	# a hack.
	if f_exp.size == 1:
	f_exp = f_exp.flat[0]
	elif f_exp.shape[0] == f_exp.size:
	f_exp = f_exp.flatten()

It appears that it does not handle this use case properly. It should be outputing an [b_N, s_N] sized array.

The fix for this should be good for any number of arguments. Test-driven development would be good here.

[sbenthall]

@nicksawhney

[sbenthall]

@llorracc did some work on this in his HARK 2.0 pre-alpha PR and is interested in talking to @nicksawhney about this.

[nicksawhney]

I would be happy to discuss whenever convenient, feel free to shoot me an email @llorracc

I’ll make sure to use test driven development @sbenthall

[llorracc]

Nick and Alan and I had a conv about this, in which I suggested that Nick do a deep dive into how Mathematica handles expectations (sympy does it in very much the same way, and I'm pretty well convinced it's because there's basically one pretty compellingly right way to do it; I recommended Mathematica only because there's great documentation and tutorial material for Mma's handing of prob and stat, and of Expectations and related topics). The deep dive will take some time -- a few days, maybe -- and I wanted to make sure you knew I'd proposed this to him because it's up to you as his real boss to decide how he allocates his time.

…

On Wed, Feb 9, 2022 at 1:13 PM nicksawhney ***@***.***> wrote: I would be happy to discuss whenever convenient, feel free to shoot me an email @llorracc <https://github.com/llorracc> I’ll make sure to use test driven development @sbenthall <https://github.com/sbenthall> — Reply to this email directly, view it on GitHub <#1098 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAKCK7ZONIFSHKTXH2AYLQDU2KVE3ANCNFSM5LRAXTVA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- - Chris Carroll

[sbenthall]

@llorracc @alanlujan91 Could either of you point to what aspect of Mathematica/Sympy you are interested in reproducing in HARK?

[alanlujan91]

I am not very familiar with Mathematica either, but I was at the meeting yesterday and I think @llorracc wants an Expectation function like the one described below

https://docs.sympy.org/latest/modules/stats.html
https://reference.wolfram.com/language/ref/Expectation.html

Expectation(Distribution(mean = mu, std = sigma, *args)) = mu
Expectation(DiscreteDistribution( **from Distribution** )) = sum( pmf * x )

The DiscreteDistribution should come from the Distribution object depending on a discretization method and a number of discrete points. It should also keep the information of how it was created (Distribution parameters, discretization method, and number of discrete points.

[llorracc]

Almost.

1. There should be Distribution objects with well defined properties and sufficient information for an Expectation operation to extract what it needs to compute the expectation.
2. The most immediate specializations of the generic "Distribution" class would be something like:
   • DistributionDiscrete
   • DistributionContinuous
3. It should be possible to construct a TransformedDistribution object which reflects the results of a function applied to a n object that has a distribution.

A Mathematica (Mma) example would be something like:

Expectation[x , x \[Distributed] LogNormalDistribution[mu,sigma]]

which would return $e^(mu + sigma^2/2)$

xLogMinusMuSq = TransformedDistribution[ (Log[x]-mu)^2, x \[Distributed] LogNormalDistribution[mu,sigma]]
Expectation[z, z \[Distributed] xLogMinusMuSq]

which would return $\sigma^{2}$. See MathFacts

(Wolfram has a free online version of Mathematica in which you can paste small snippets of Mma code to try them out; it should be possible to copy and paste these there and get the promised results).

sympy does exactly the same thing (with slightly different syntax), and from a brief examination of it I believe TensorFlow permits the same thing (though for some reason they call a TransformedDistribution a Bijection). Probably pytorch can do the same thing with some syntax or other.

Both Mathematica and sympy will allow us to leave the objects in symbolic form, and are capable of doing symbolic manipulations on them while they remain in symbolic form. That is, if sigma and mu are undefined at the time the expressions above are evaluated, the output from Mma will literally be the string formula $e^(mu + sigma^2/2)$ while if before the first expression is evaluated we set something like sigma = 0.1 and mu = -(sigma^2)/2 we will get 1.0 as the answer.

Probably the best way to do this would just be to make sympy a dependency and use sympy's syntax. Or TensorFlow or maybe pytorch.

Once we adopt some foundational architecture, we can build on top of it the things we've been wanting to build, like a system for generating and keeping track of discretizations.

cc'ing @Mv77

[sbenthall]

I think the main thing blocking practical exploration of this space is #1105

[sbenthall]

From discussion with @llorracc today:

• Look at: https://www.tensorflow.org/probability/examples/TensorFlow_Distributions_Tutorial
• Bijections.
• TransformedDistribution in Mathematica

[llorracc]

Just wanted to make sure @nicksawhney saw this

[Smit-create]

Probably the best way to do this would just be to make sympy a dependency and use sympy's syntax.

Yeah, sympy can handle the custom distribution classes(the one which we can define symbolically according to our need) too which we can use to create our own distribution. I'll be happy to provide help in migrating the distribution classes under sympy dependency. See: https://docs.sympy.org/latest/modules/stats.html#examples

[sbenthall]

@Mv77 knows how to fix this, has comment in #1146