I want to use Pymc to estimate the parameters of Dynamic Discrete Choice Model. The model consists of states and action, and the utility function is defined as

where s is state and p is action.

Likelihood Function is defined as

And the measure of model fit in Bayesian statistics is the marginal density of the

data according to the model, defined as

can I use pymc to estimate the model?

If I’m reading the setup right, the challenge I see is that you need to know the value function by solving (3) every time you change the parameters. Do you have a closed form solution for it? Otherwise you will need to do VFI (or similar) on each parameter draw. PyMC can do this, but not in a way that would give you gradients (i.e. you lose NUTS).

Once you have V_a(s | \theta_i), though, PyMC could do it no problem. The solution would be a table of size `(states, actions)`

, and (5) is just a softmax over the action dimension. Then you can use a `pm.Potential`

to insert the P(\text{data} | \theta) term into the likelihood. PyMC would handle the prior part for you.

So VFI is Value Function Iteration? do you now how to implement that in pymc?

Yes, VFI is value function iteration. There is a very nice collection of lectures about numerical solutions to Bellman equations on QuantEcon, starting from lecture 43 until lecture 48.

Once you had a function to implement VFI, you would wrap it in a Pytensor `Op`

, and call it within your model, similar to this example. In the example, the author implements a likelihood function puts it into a `pm.Potential`

, but you can follow the exact same procedure to run arbitrary python code inside your PyMC model.

do you know how to write the softmax likelihood function?

Use `pm.Categorical`

as the likelihood. The softmax computes the class (choice) probabilities, it’s `pm.math.softmax`

.

Thanks! since there are may states=, should I write multiple `pm.Categorical`

?

No, you should gather everything together. From the equations it looks like you are trying to infer actions? So the states should just be integer indexes, and your data should be a matrix with dimensions `(observation, time)`

, is that right?