I am trying to use pymc to estimate parameters for a model. As I am not familiar with the methodology I first generated the data myself and tried to derive the posterior distribution using pymc.

The data consists of a (time-) series of the number of missed payments $D$ for each quarter. For each quarter, the probability of a client missing a payment is a random variable, $p(x,\alpha)$, that depends on two parameters, $\alpha$ which is an unknown constant (in time) that I'm trying to estimate, and $x \sim \mathcal{N}(0,1)$, which is a random variable and drawn for each quarter.

Edit: for clarity, the function $p(x,\alpha)$ below describes the probability of a client to miss a payment. This probability is normally 3% per quarter but fluctuates a bit however, based on a correlation between clients ($\alpha$) and a random number. Please note each quarter is independent from the other. Only $x$ makes this a random variable. I have coded this as a deterministic in PyMC as its value it determined completely by $\alpha$ and $x$.

The actual probability of missing a payment in a quarter is described by the following function:

$$p(x,\alpha) = \Phi \left( \frac{ \Phi^{-1}(rate) - \alpha x }{\sqrt{1-\alpha^2}}\right) $$

where $\Phi$ is the normal cumulative distribution function and $\Phi^{-1}$ its inverse.

To generate a data sample for 50 quarters I have used this code (with the sample $\alpha = 0.37$ ).

sample_alpha = 0.37
const_z = 0.03 # basis prob of missing payment (3%)
nr_clients = 10000
nr_quarters = 50

Xs = np.random.normal(loc=0.0, scale=1.0,size=nr_quarters)

def p(x,a):
    return sp.norm.cdf( (sp.norm.ppf(const_z) - sample_alpha*Xs)/(np.sqrt(1-sample_alpha*sample_alpha)) )

sample_mp = np.random.binomial(n=nr_clients, p=p(Xs,sample_alpha))

py.ylabel('nr of missed payments')

Below is shown a random realization of this process. To give more insight: every quarter the random variable $x$ and the actual $\alpha$ produce a certain probability of missing a payment. What is shown below for each quarter is then a drawing from a random binomial variable where the probability of success is this probability of missing a payment.

enter image description here

Now I'm trying to estimate the posterior distribution of $\alpha$ from this realization. I am using the following pymc code:

xfactor = pm.Normal('xfactor', mu=0, tau=1/(1*1), size = sample_mp.shape[0])

alpha = pm.Uniform("alpha", 0.01,0.99)

def P(x=xfactor,a=alpha):
    return sp.norm.cdf( (sp.norm.ppf(const_z) - a*x)/(np.sqrt(1-a*a)) )

missed_payments = pm.Binomial('defaults',n=10000,p=P, value=sample_mp,observed=True)
model = pm.Model([xfactor, alpha, P, missed_payments])

And then try to run it to get some feeling of what is going on (I haven't set burn-in and thinning at this point):

mcmc = pm.MCMC(model)

The results look like this:

enter image description here

When I increase the nr of samples the plot is similar. I understand from this that $\alpha$ doesn't seem to converge, and the autocorrelation plot looks problematic too, but I am uncertain how to resolve the issue.

  • $\begingroup$ Can you explain your model (esp. sp.norm.cdf part)? What is the parameter alpha and what it governs? $\endgroup$ – Vladislavs Dovgalecs Feb 12 '17 at 23:31
  • $\begingroup$ Also, since your data is time-series, RandomWalk would be perhaps a better fit (for real-world data). But you're right, this model should fit the synthetic dataset where observations are i.i.d. $\endgroup$ – Vladislavs Dovgalecs Feb 12 '17 at 23:33
  • $\begingroup$ I think Negative Binomial distribution would fit better the problem. $\endgroup$ – Vladislavs Dovgalecs Feb 13 '17 at 2:33
  • $\begingroup$ @xeon: I have added the information in the post. I will have a look at RandomWalk; although since the realizations of p are independent from eachother I'm not sure how this would work out. $\endgroup$ – Letranger2014 Feb 13 '17 at 8:03
  • $\begingroup$ As for the negative binomial: the number of clients is fixed (10,000) and the nr of clients that miss a payment fluctuates per quarter (base rate = 3%). For extra info: alpha here governs an intrinsic correlation between clients, so a higher alpha results in a higher p(x,a), and thus for the same x, a higher alpha means a higher p(x,a). I was thinking this (pymc-devs.github.io/pymc/…) might be of interest but I've been stuck on it for a day or 2 now. Any input is very welcome.. $\endgroup$ – Letranger2014 Feb 13 '17 at 8:10

Here is my shot at the problem in PyMC3. I can be wrong how the model is built, so please correct me where I am wrong.

The data are 50 observations (50 binomial draws) that are i.i.d. This assumptions is strong one.

I set the true parameter value (p_true=0.37) and set number of Bernoulli trials to 10,000. With these parameters set, I draw (K=50) samples from the Binomial distribution.

import pymc3 as pm
import matplotlib.pyplot as plt

from scipy.stats import binom

p_true = 0.37
n = 10000
K = 50

X = binom.rvs( n=n, p=p_true, size=K )
print( X )

model = pm.Model()

with model:
    p = pm.Beta( 'p', alpha=2, beta=2 )
    y_obs = pm.Binomial( 'y_obs', p=p, n=n, observed=X )
    step = pm.Metropolis()
    trace = pm.sample( 10000, step=step, progressbar=True )

pm.traceplot( trace )
| cite | improve this answer | |
  • $\begingroup$ Thanks a lot for your help. I will take some time to study your version. There seems to be 1 misunderstanding which is that the p_true as in your model is not actually fixed: it is also a random variable. In my generation code above you can see that as I first produce the x, that lead to a different p_true for each drawing. Please note that the drawings indeed are independent as the x's are independent. $\endgroup$ – Letranger2014 Feb 13 '17 at 8:31
  • $\begingroup$ @Letranger2014 Thanks for your feedback! It looks like you DO want to have time-series as your "p" changes from quarter to quarter. Here you would want to use RandomWalk. $\endgroup$ – Vladislavs Dovgalecs Feb 13 '17 at 9:46
  • $\begingroup$ Indeed, P changes from quarter to quarter. In fact, p is in [0,1] and its pdf can be derived analytically -- it's just somewhat of an ugly function involving Erf and inverse Erf). I therefore opted for sampling x and producing p as a result. I havent succeeded in trying RandomWalk as the feature doesn't seem to be documented very well. $\endgroup$ – Letranger2014 Feb 15 '17 at 14:22
  • $\begingroup$ Take a look at this example of RandomWalk github.com/pymc-devs/pymc3/blob/… and pymc-devs.github.io/pymc3/notebooks/survival_analysis.html $\endgroup$ – Vladislavs Dovgalecs Feb 15 '17 at 22:37
  • $\begingroup$ @Letranger2014 I would try replacing the "p = pm.Beta(...)" with "p = pm.GaussianRandomWalk(...)". This would model p at time t conditioned on p at time t-1: prob(p_t | p_t-1) $\endgroup$ – Vladislavs Dovgalecs Feb 15 '17 at 22:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.