Trying to understand pymc3 through a simple example

jaredisathome · December 18, 2018, 5:45am

Hi Guys and Girls,

Sorry for the low level questions here - but I am working through trying to grasp pymc… I am hoping that someone can explain to me the following through an example, that helps me to deepen my understanding. At its core i am trying to work through the following question:

on the assumption we’re working with a time series, is it possible, at time t(n) to reference t(n-1) or t(n-2) etc.

To put it into a practical context - please humor me and assume that this is a real world application of pymc3 - given the following simple example. There may be a better example thought…

If want to model revenue growth through over the next ten years:

Using numpy, i create a 10x1000 matrix of revenue growth, sampled from the normal distribution with mean of .12 and variance of 0.03.

growth = np.random.normal(0.12,0.03,(10,1000))

And therefore my revenue growth is modeled as

revenue = (1+sample).cumprod(0)

Any help would be appreciated in modelling this in pymc3 so i can understand how, sampling from distributions randomly can then be applied into the above function, which may be stored as a deterministic variable…-> even if it pointing me in the right direction.

Cheers

junpenglao · December 18, 2018, 5:59am

If you write down something like

with pm.Model() as model:
  growth = pm.Normal('growth', mu=0.12, sd=0.03, shape=(10,1000))

It will generate a RandomVariable growth, which you can understand as a theano.tensor that draw randomly from a normal distribution. You can then do something like revenue = (1+growth).cumprod(0).

jaredisathome · December 18, 2018, 6:02am

thank you junpenglao, i’ll try that !!! who would have thought it was that easy…

further to that - the reason i defined the matrix in numpy was so simulate what i thought would be done naturally by pymc3 when sampling from the distributions anyway…

its likely i would misunderstand pymc3 but i was guessing that you dont need to dim multiiple dimensions (ie, forecast 10 years, 1000 times as per the above, as thats what pymc3 does regardles…

so in this instance, defining shape means you’ll generate samples x matrix(10,1000) worth of samples, correct?

junpenglao · December 18, 2018, 8:02am

Well, it depends - most likely, if you considering each year to be a repeated observation, you would generate a shape=(1, 1000) and account for the 10 years in some random variable downstream.

jaredisathome · December 21, 2018, 1:08am

@junpenglao

thanks… im still a bit confused here, so i’m hoping this illustrates my problem further. So, i’m wondering if you could help me with some code… this should help me understand this problem. or point me to the correct docs…

As an example, I want to model growth in revenues that looks like this over a 20 year period. (see fade function below)

a = []
for i in range(20):
    a.append(fade2(.15,.03, 20, i))
plt.plot(a)

This revenue growth rate, converts to actual revenue as follows:

b = np.array(a)
plt.plot((1+b).cumprod())

So, im trying to understand how to use pymc3 to do this kind of timeseries stuff, because of the ability to do the beysian inference… the examples that i can find are a bit complex for me at this stage, im tying to get basic building blocks out…

so, at this stage i have the following model defined as follows:

with pm.Model() as revModel:
    growth_st = pm.Normal('growth_st',15,1)
    growth_lt = pm.Normal('growth_lt',3,.05)
    t_halfway = pm.Normal('t_halfway',5,1)

the fade function (used above) is defined as follows:

def fade(start, end, periods, t):
    halfway = periods / 2
    return start + (end - start) / ( 1 + np.exp(-(t - halfway)))

how do i go from the above snippets of code, to generate the same thing, but convert it to the model in pymc3… i tried the following:

rev_growth = []
for i in range(10):
    rev_growth.append(fade(growth_st, growth_lt, t_halfway, i))

and got the following here after printing it…

[Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0, Elemwise{add,no_inplace}.0]

am i on the right track?

thanks so much

junpenglao · December 21, 2018, 11:27am

Appending a python list and using a python for loop is usually pretty slow. In this case, you can vectorized your operation:

t = np.arrange(10)
halfway = t_halfway / 2
growth_st + (growth_lt - growth_st) / ( 1 + pm.math.exp(-(t - halfway)))

jaredisathome · December 30, 2018, 12:17am

thanks for your advice… vectorizing the function solved my issue with building out the model.

I’ve now built out the model that I want, and have two further questions to help my understanding,.

For background, the model is describing afree cash flow model to value a company. The fade function is used to help fade revenue and ebit etc, to change through time… the starting and ending points are determined by the priors rev_g_st and rev_g_lt and ebit_m_st and ebit_m_lt…

def fade(start, end, periods):
    #vectorised function.. ie, its passing i as a vector, not as a loop.
    t = 1+np.arange(periods)
    halfway = periods / 2
    return start + (end - start) / ( 1 + np.exp(-(t - halfway)))

with pm.Model() as fcf_model:    
    #Priors
    #revenue and margins
    rev_g_st = pm.Normal('rev_g_st',.25,.01)
    rev_g_lt = pm.Normal('rev_g_lt',.12,.04)
    ebit_m_st = pm.Normal('ebit_m_st', -0.2, 0.01)
    ebit_m_lt = pm.Normal('ebit_m_lt', 0.22, 0.03)
    
    #working capital turn, and depreciation and amortization as a percentage of cashflow
    wc = pm.Uniform('wc', 0.05, 0.09)
    da = pm.Uniform('da', 0.07, 0.10)
    
    #other inputs
    wacc = pm.Normal('wacc',0.12,0.02)
    tax = pm.Normal('tax',0.27,0.01)
    lt_g = 0.05
    period = 10
    shs_on_issue = 160
    BASE_REVENUE = 75
    
    #calculate free cash flows for the interval described by the period - ie years 1 to 10 - at each sample.
    rev_growth = pm.Deterministic('rev_growth',fade(rev_g_st, rev_g_lt, period))
    rev = pm.Deterministic('rev', (1+rev_growth).cumprod(0) * BASE_REVENUE)
    ebit = pm.Deterministic('ebit', rev * fade(ebit_m_st, ebit_m_lt, period))
    fcf = pm.Deterministic('fcf', ebit*(1-tax) - wc - da)
    
    #terminal cash flow - at each sample. 
    term_fcf = pm.Deterministic('term_fcf', fcf[period-1]/(wacc-lt_g))
    
    #calculate present value - for all years described by period, at each sample.
    disc_factor = 1/((1+wacc)**(np.arange(period)))
    pv_fcf = pm.Deterministic('pv_fcf',fcf * disc_factor)
    pv_term_fcf = pm.Deterministic('pv_term_fcf', term_fcf / ((1+wacc)**(period+1)))
    
    #calculate value per share
    value = pm.Deterministic('value', (pv_fcf.sum() + pv_term_fcf)/shs_on_issue)

If I plot the values sampled from the model for the deterministic functions, they follow the general shape of that I would expect.

Question 1:

In the first instance, I know I am doing a sort of monte carlo type of simulation, and so am only really interested in the deterministic variables. Note there is no observed data in the model yet.

For the deterministic variable - ‘value’ - plotted as a histogram at the bottom left of the figure, I want to work out the probability of a specific value occurring, and so tried the following:

value.logp({'value':1.445})

but got the following error:

AttributeError: 'TensorVariable' object has no attribute 'logp'

Do we have to fit a distribution to this to calculate the pdf or can pymc3 do this automatically?

Question 2:

The histogram at the bottom right of the figure is actual observed prices for the security I am modelling out above using the theoretical FCF model.

I am interested to know, that given my model above, if we are able to work backwards, given the observed data and work out the rev_g_st, rev_g_lt, ebit_m_st, ebit_m_lt etc using beysian inference. (I am aware that these are transformed using the fade function to get the deterministic variables rev_g and ebit…)

I tried setting the observed parameter in the deterministic variable as follows:

    value = pm.Deterministic('value', (pv_fcf.sum() + pv_term_fcf)/shs_on_issue, observed=obs_data)

but that didnt work…

Is what I am trying to achieve here possible?? It feels like to me that it is, based on my very limited knowledge of pymc3… thanks for all suggestions so far, its really helping me deepen my knowledge of pymc3.

junpenglao · December 30, 2018, 2:17am

Your two questions are related, and both do not have a good solution unfortunately: you want to associate logp and observed to a deterministic variable. This could not be done because deterministic variable does not have a logp function. The current work around is to model it using a Gaussian with small as:

value = pm.Normal('value', (pv_fcf.sum() + pv_term_fcf)/shs_on_issue, .01, observed=obs_data)

jaredisathome · December 30, 2018, 2:55am

Thanks! I’m assuming we can use any distribution there… like gamma for example???

junpenglao · December 30, 2018, 3:20am

That’s even better if your observed is from some distribution.

pthaz · August 22, 2022, 1:06pm

Hey Jared. I am trying to recreate your steps for a similar DCF model. Can you help with the following?
a. How did you sample from the priors (since you dont have a likelihood defined)? Was it using sample_prior_predictive
b. How did you plot those charts for each of the variables. Was it using something built into pymc3 or custom plots based on the trace of each variable ?

Topic		Replies	Views
Help with pm.Potential Questions	10	1130	January 13, 2021
Help with Hierarchical Model Example Questions	3	498	October 5, 2021
How to interpret a very basic PyMC3 model? Questions	2	1025	November 15, 2021
How to model time-dependent variables in pymc3 Questions	5	4037	November 6, 2017
New to pymc3 - learning through example - financial valuation question Questions	3	878	October 11, 2018

Trying to understand pymc3 through a simple example

Related topics