# Trying to understand pymc3 through a simple example

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

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)`.

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?

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.

@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):
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):
``````

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 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)))
``````

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 = 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.

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)`

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

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