# Trying to understand pymc3 through a simple example

**URL:** <https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404>\
**Category:** Questions\
**Created:** [December 18, 2018, 5:45am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404 "2018-12-18T05:45:22Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![jaredisathome](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jaredisathome/32/1170_2.png) [@jaredisathome](https://discourse.pymc.io/u/jaredisathome)\
**Post date:** [December 18, 2018, 5:45am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/1 "2018-12-18T05:45:23Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 18, 2018, 5:59am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/2 "2018-12-18T05:59:55Z")

</div>

> [@jaredisathome](#):
>
> sampling from distributions randomly can then be applied into the above function

If you write down something like

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

---

<div class="post-metadata">

**Author:** ![jaredisathome](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jaredisathome/32/1170_2.png) [@jaredisathome](https://discourse.pymc.io/u/jaredisathome)\
**Post date:** [December 18, 2018, 6:02am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/3 "2018-12-18T06:02:03Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 18, 2018, 8:02am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/4 "2018-12-18T08:02:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jaredisathome](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jaredisathome/32/1170_2.png) [@jaredisathome](https://discourse.pymc.io/u/jaredisathome)\
**Post date:** [December 21, 2018, 1:08am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/5 "2018-12-21T01:08:14Z")

</div>

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

```

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/8/89532031c10405b21be7f4e9ef85564ed756d1f4.png)

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

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

```

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/4/48cfb9a119f9636102d05134fa1c022085bc527f.png)

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 🙏

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 21, 2018, 11:27am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/6 "2018-12-21T11:27:32Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![jaredisathome](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jaredisathome/32/1170_2.png) [@jaredisathome](https://discourse.pymc.io/u/jaredisathome)\
**Post date:** [December 30, 2018, 12:17am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/7 "2018-12-30T00:17:20Z")

</div>

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 a[free cash flow model](https://en.wikipedia.org/wiki/Free_cash_flow) 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.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/0/0476e15993a21abab87d6bb19a7888aed2448adf.png)

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

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 30, 2018, 2:17am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/8 "2018-12-30T02:17:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jaredisathome](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jaredisathome/32/1170_2.png) [@jaredisathome](https://discourse.pymc.io/u/jaredisathome)\
**Post date:** [December 30, 2018, 2:55am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/9 "2018-12-30T02:55:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 30, 2018, 3:20am UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/10 "2018-12-30T03:20:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![pthaz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/pthaz/32/5557_2.png) [@pthaz](https://discourse.pymc.io/u/pthaz)\
**Post date:** [August 22, 2022, 1:06pm UTC](https://discourse.pymc.io/t/trying-to-understand-pymc3-through-a-simple-example/2404/11 "2022-08-22T13:06:45Z")

</div>

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 ?