customer_zakka.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Analysis of Zakka  Canada Customers\n",
    "\n",
    "In this notebook, I take a look at customer purchase information from a real jewelery display company called Zakka Canada.    I apply the standard BG/NBD model from _lifetimes_ python package to fit and perform inference on the data.  The data details all orders from 2007 to 2015."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from lifetimes.plotting import *\n",
    "from lifetimes.utils import *\n",
    "from lifetimes.estimation import *\n",
    "from pandas import read_csv\n",
    "import pandas as pd\n",
    "from scipy.stats import beta, gamma, invgamma, norm\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import datetime as dt\n",
    "import warnings\n",
    "from random import uniform \n",
    "from statsmodels.api import OLS\n",
    "import numpy as np\n",
    "warnings.filterwarnings('ignore')\n",
    "%cd\n",
    "sns.set_palette(\"husl\")\n",
    "sns.set(rc={'image.cmap': 'coolwarm'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trans_data = read_csv('Downloads/orders.csv')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trans_data['Date'] = pd.to_datetime(trans_data['Date'])\n",
    "trans_data.head()\n",
    "print 'Number of Entries: %s' % len(trans_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The models presented in _lifetimes_ measure customer's __repeat__ purchase behaviors based on his/her recency $x_t$, frequency $x$ and analysis period $T$.  Thus, each customer can be summarized into a vector $(x, t_x, T)$.  We can take our basic customer transaction list and parse it into a usable RF(M) dataframe for the model. We should use the function already made within _lifetimes_ called __summary_data_from_transaction_data()__:\n",
    "- Recency ($t_x$) is calculated as the most recent transaction within $T$ less time of the first transaction\n",
    "- Frequency ($x$) is number of repeat transactions the customer has made\n",
    "- $T$ is equal to the end observation period less first purchase made by the customer (implied age)\n",
    "\n",
    "The unit of time is in days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = summary_data_from_transaction_data(trans_data, 'Customer ID', 'Date', monetary_value_col='Subtotal', observation_period_end='2015-12-31')\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data[data['frequency'] > 0]['frequency'].plot(kind='hist', bins=20)\n",
    "print data[data['frequency'] > 0]['frequency'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most customers are one time buyers on Zakka Canada with only ~13.6% of customers making two or more repeat purchases while 72.3% of customers only make one purchase.  Below we look at their recency distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data[data['recency'] != 0]['recency'].plot(kind='hist', bins=20)\n",
    "print sum(data['recency'] != 0)/float(len(data))\n",
    "print data[data['recency'] != 0]['recency'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Customers that previously have a frequency of zero (customers that have not made repeat purchases) also have a recency of zero so we isolate for repeat customer only.  Considering the analysis ranges within a span of eight years, having 50% of your customer churn before a year (245) is even up and 75% of your customers churn before two years (661) should be a concern"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the BG/NBD Model\n",
    "\n",
    "In quick summary, the BG/NBD model assumes that each customer makes repeat purchases at different intervals.  This interval between purchases is modeled by an exponential distribution with a transaction rate $\\lambda$.  Furthermore, after each transaction, the customer has a $p%$ chance of going inactive (or dying) which is modelled by a shifted Geometric distribution.  These two major parameters are further refined for heterogenity across customers through a Gamma distribution for $\\lambda$ and a Beta distribution for $p$.  Thus we have four total parameters $r$, $\\alpha$ (gamma distro) and $a$, $b$ (Beta distro) estimated below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# similar API to scikit-learn and lifelines.\n",
    "bgf = BetaGeoFitter(penalizer_coef=0.0)\n",
    "bgf.fit(data['frequency'], data['recency'], data['T'], )\n",
    "print bgf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since everything goes back to the two major parameters, we can use our heterogenity parameters to model the uncertainty in the two.  As shown below, the probability of a customer dropping out has a large range but mostly centered between 15 to 60% after a certain transaction.  On the other hand, our customer transaction rate falls mostly within 0 and 0.01."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "gbd = beta.rvs(bgf.params_['a'], bgf.params_['b'], size = 50000)\n",
    "ggd = gamma.rvs(bgf.params_['r'], scale=1./bgf.params_['alpha'], size = 50000)\n",
    "plt.figure(figsize=(14,4))\n",
    "plt.subplot(121)\n",
    "plt.title('Heterogenity of $p$')\n",
    "temp = plt.hist(gbd, 20, facecolor='pink', alpha=0.75)\n",
    "plt.subplot(122)\n",
    "plt.title('Heterogenity of $\\lambda$')\n",
    "temp = plt.hist(ggd, 500, facecolor='pink', alpha=0.75)\n",
    "plt.xlim(0.0,0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fact Checking\n",
    "\n",
    "Peter Fader in his [Youtube talk](https://www.youtube.com/watch?v=guj2gVEEx4s) on CLV provides four simple plots to illustrate the accuracy of the model and its associated behavioral story when applied to real data sets.  I will replicate two of them below on our data to ensure that the model is indeed a good fit for our customers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_period_transactions(bgf, max_frequency=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After fitting the BG/NBD model, we can fact check through looking at the fitted repeat purchase for the entire customer database.  Above plots the Actual repeat purchase frequency versus the Model repeat purchase frequency.  As shown, the model performs relatively well in fitting the repeat purchase pattern. It overestimates one and six repeat purchases while underestimating three repeat purchases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "summary_cal_holdout = calibration_and_holdout_data(trans_data, 'Customer ID', 'Date', \n",
    "                                        calibration_period_end='2015-01-01',\n",
    "                                        observation_period_end='2015-12-19' )   \n",
    "bgf.fit(summary_cal_holdout['frequency_cal'], summary_cal_holdout['recency_cal'], summary_cal_holdout['T_cal'])\n",
    "plot_calibration_purchases_vs_holdout_purchases(bgf, summary_cal_holdout, n=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can further validate our model by checking for presence of overfitting.  Separating our dataset into two sets: training (2007- 2014 end) and test (2015), we can fit our model with training and see how well it forecasts our customer's out of sample purchasing habit.  The plot above shows our fitted model prediction (green) versus our out of sample repeat purchase frequency (blue).  Each tick on the x-axis is composed of customers that have been grouped such that they have had $x$ repeat purchases.  the Y axis shows the mean of that group's subsequent repeat purchase.  For example, for all customers within our training period that have had 4 repeat purchases, the same group had on average about 0.42 purchases and the model predicted the same as well.  From this plot, we can see that the model performs relatively well in identifying repeat purchases within the lower ranges (less than 6) but is unable to properly forecast buying behaviors of customers with more frequent purchasing habits. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R/F Matrix Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_frequency_recency_matrix(bgf, T=365)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a given random customer with his/her observed $(x,t_x)$, we can calculate through the model the expected number of purchases within a time period which I have set equivalent to a year.  As shown, customers that purchase frequently and recently have a much higher expected future repeat purchase although this is very slim region on the chart.  For example, unless a customer has purchased recently (say within the past year and a half) and frequently (more than 5 times), don't expect them to show up and place an order within the next year."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_probability_alive_matrix(bgf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a given random customer with his/her observed $(x,t_x)$, we can calculate the probability as of today that they are still alive.  In other words, are they still interested in maintaining a purchasing relationship with our company in the future?  According to the model, it's very likely as long as they bought recently.  Even if they only bought once, they have about a 60-70% chance that they will buy again in the future.  Customers that are considered dead are usually ones that haven't made an order in a while.  \n",
    "\n",
    "It is interesting to note that for a given recency, customer's with more frequency are more likely to be considered dead.  This is a property of the model that illustrates a clear behavioral story:  If we observe customers making more frequent purchases, then having a prolonged period of inactivity makes it very likely that they've died off.  For example, Two customers $A$ and $B$ that both last purchased 1.5 years ago ($t_x\\approx 2700$) but purchased 10 and 30 times respectively.  We believe that it is more likely that customer B has died ($p=0.0$) while customer A still has a fair chance of being alive and making purchases in the future ($p\\approx 0.4$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing Customer's Historical Path\n",
    "\n",
    "One last plot that was developed in __lifetimes__ allows us to fully appreciate the model at a granular customer-base level.  Indeed, we can observe each of our customer's historical purchasing path and measure our confidence that they have not dropped out in between purchases.  To get a small and relevant subsample, I will let the _bgf_ model create projections for the next month purchasing pattern of my entire current customer-base and pick out 6 of my best customers.  As shown, our best customers that are projected to make more future repeat purchases are often our newest ones.  This makes sense since the model has revealed to us before that our customer base has a tendency to drop off after a short period of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t = 31\n",
    "data['predicted_purchases'] = bgf.conditional_expected_number_of_purchases_up_to_time(t, data['frequency'], data['recency'], data['T'])\n",
    "best_projected_cust = data.sort('predicted_purchases').tail(6)\n",
    "print data.sort('predicted_purchases').tail(6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we can plot these six individuals over the course of their lifetime to observe their purchasing behavior.  Among these top 6 customers, they are sorted by expected future purchases ascendingly.  Adding a one month cushion, we can infer the probability that these customers are alive or dead if we have not seen them purchase by the end of the month.  For example, for customer $\\tt{14263}$ with a short life-time but very recent purchases, we expect him/her to  return and buy at-leasts once during the month of January.  If however, we don't see this person make an order by the middle of Jan 2016, there is a good chance that they're dead.\n",
    "\n",
    "Looking at the overall picture, there's not much pattern that can be generalized across all 6 individuals but we can identify certain seasonal behaviors for certain people, for example: \n",
    "\n",
    "1. Customer $\\tt{13577}$ tends to buy at an equal interval with the exception between September 2015 and Nov 2015.\n",
    "2. Customers $\\tt{12780}$,  $\\tt{13527}$ and $\\tt{10920}$ tends to buy during the early stage of their life-cycle and then repeat at a much slower interval\n",
    "3. Customer $\\tt{14008}$ first bought at a slower interval but increased his/her frequency over time. We can reasonably assume that they started bying more as their business expanded and became more successful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,12))\n",
    "for ind,i in enumerate(best_projected_cust.index.tolist()):\n",
    "    ax = plt.subplot(4,2,ind+1)\n",
    "    best_T = data.ix[i]['T']+31 #add a month\n",
    "    best_trans = trans_data[trans_data['Customer ID'] == i]\n",
    "    plot_history_alive(bgf, best_T, best_trans, 'Date', freq='D', ax=ax)\n",
    "    ax.set_title('ID: '+str(i))\n",
    "    plt.xticks(rotation=25)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Customer Value in Dollars: The Gamma-Gamma Model\n",
    "In a traditional RFM matrix, each customer can be assigned a monetary value with his or her customer profile.  This value is usually calculated as the mean transaction value of his or her purchasing history.  Transaction value is a loosely based word around any metric that the business might be interested in.  Zakka Canada reports an average profit margin of 55%, I will multiply that by the revenue.  Note that we are only interested in customers that have repeat purchases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ggf = GammaGammaFitter(penalizer_coef = 0)\n",
    "ret_cust_data = data[data['frequency'] > 0]\n",
    "ret_cust_data['monetary_value'] = ret_cust_data['monetary_value']*0.55 # Set Monetary value to gross profit only\n",
    "ggf.fit(ret_cust_data['frequency'], ret_cust_data['monetary_value'])\n",
    "p,q,v = ggf._unload_params('p', 'q', 'v')\n",
    "print ggf "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.special import beta\n",
    "def pzbar(zbar,p,q,v,x):\n",
    "    return (1./(zbar*beta(p*x,q)))*((v/(v+x*zbar))**(q))*(((x*zbar)/(v+x*zbar))**(p*x))\n",
    "\n",
    "x = range(1,100)\n",
    "numtrans = [1,3,5,25,100]\n",
    "dist = pd.DataFrame(0,index=x, columns=numtrans)\n",
    "for j in range(0,5):\n",
    "    for i in x:\n",
    "        dist.loc[i,numtrans[j]] = pzbar(i, p,q,v,numtrans[j])\n",
    "dist.plot(figsize=(10,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The authors now introduce another complicated model to inferring transaction value called the __Gamma-Gamma Model__.  The GG model behaves much like the underlying concepts driving the BGF model: _heterogenity_.  What is the mean transaction value of our customer base? Well we can often estimate that with the sample mean and variance...IF our customer base was symmetrically and normal distributed.  Since this isn't the case (1. non-negative transaction value 2. right skew), we can use the Gamma distribution to approximate the distribution of the mean transaction value ($\\bar{z}=\\sum_{i} z_i / x$) for EACH customer:\n",
    "$$\\bar{z}\\sim Gamma(px,\\nu x)$$ Allowing the shape parameter ($\\nu$) to be heterogeneous over the entire customer base, here comes the second gamma term:\n",
    "$$\\nu\\sim Gamma(q,\\gamma)$$\n",
    "\n",
    "We can calculate the conditional average profit of our customer base versus the observed average profit.  In short, the conditional average profit uses the same concept of unobserved average transaction value in the distribution and provide a shrinkage mean between the unobserved mean and the customer's sample mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "avg_cond_profit = ggf.conditional_expected_average_profit(ret_cust_data['frequency'], ret_cust_data['monetary_value'])\n",
    "print \"Expected conditional Average profit: %s, Population Average profit %s, Average profit: %s\" % (\n",
    "    avg_cond_profit, \n",
    "    (p*v)/(q-1),\n",
    "    ret_cust_data['monetary_value'].mean()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In short, the shrinkage trans. value has the highest average profit followed by what we have observed and lastly the population distribution.  Although these differences are minute, it still implies that the model believes customers are currently slightly under purchasing to their true transaction value.  The shrinkage model also allows us to closely examine the distribution on a customer by customer basis.  Recall that we previously discovered our six best customers for the month of January, let's look at their transaction value we expect from them over the next month."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "customer_scale = p*(v+ret_cust_data['frequency']*ret_cust_data['monetary_value'])\n",
    "customer_shape = p*ret_cust_data['frequency']+q\n",
    "ez = customer_scale/(customer_shape-1)\n",
    "best_projected_cust['predicted_purchase_value'] = best_projected_cust['predicted_purchases']*ez\n",
    "plt.figure(figsize=(12,6))\n",
    "plt.xlabel('$ Profit')\n",
    "for ind,i in enumerate(best_projected_cust.index.tolist()):\n",
    "    ccig = invgamma.rvs(customer_shape[i], scale=customer_scale[i], size = 10000)\n",
    "    plt.text(ez[i]+5, invgamma.pdf(ez[i], customer_shape[i], scale=customer_scale[i]),\n",
    "             '$'+str(round(best_projected_cust['predicted_purchase_value'][i],2))\n",
    "            +'| '+str(int(ret_cust_data.loc[i,'frequency'])))\n",
    "    sns.distplot(ccig, hist=False, rug=False, label=str(i));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The legends are in the order from least to most expected future purchases (in quantity, not trans. value) for our top 6 customers.  After factoring in their underlying mean purchase distribution, we can gain a lot of insight:\n",
    "- The order of ranking changes, customer $\\texttt{14008}$ is now our second best customer based on his transaction value.\n",
    "- $\\texttt{13577}$ and $\\texttt{10920}$ are frequent small buyers, they make frequent purchases worth approximately \\$21 to \\$28.  Thus, we are more confident about their range of mean transaction prices.\n",
    "- Conversely, $\\texttt{14263}$ being our best customer, purchase infrequently but in large quantities, we are less sure of his mean transaction value but he still stochastically dominates other top customers. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Valuation: Estimating From the Ground Up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Potential investors interested in a promising company often need to determine how much the company should cost at the present day.  According to traditional financial theory, the cost (also referred to as worth or value) of a company should be equivalent to its discounted future cash flow net of all costs and investments (more on this on the blog).\n",
    "\n",
    "For Zakka-Canada, it's pretty simple:  The company incurs a profit margin on each order (which I have established as 55% beforehand) and pays a fixed cost each month.  Taking this, we can make a few assumptions on their operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. I took the 2015 annual cost data, divided into equal month payments and used this as the estimate for future fixed cost."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "monthly_fc = 195000/12 # Estimated from 2015 annual data\n",
    "print 'Estimated Fixed Cost per month: $%s' % monthly_fc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.  Average annual equity return ranges 6-8%, if Zakka Canada were to sell off their operations today to investors, we can assume a 15% annualized return for the new owners due to its small size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dr = 0.15 # An annualized discrete rate\n",
    "mr = (1+dr)**(1/12.)-1. # Monthly effective rate\n",
    "print 'Effective Monthly Rate: %s%%' % mr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. New customer sales have risen historically but stagnated over time with cyclical purchase trends."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "unique_cohort = trans_data[['Customer ID', 'Date']].sort('Date')\n",
    "unique_cohort = unique_cohort.drop_duplicates('Customer ID')\n",
    "unique_cohort['Date_M'] = unique_cohort['Date'].apply(lambda x: dt.datetime(day=1,month=x.month,year=x.year))\n",
    "cohort_over_time = unique_cohort[['Date_M', 'Customer ID']].groupby('Date_M').count()\n",
    "ax = cohort_over_time.plot(legend=False, figsize=(12,4))\n",
    "ax.set_xlabel('Date'); ax.set_ylabel('New Customer Cohort Size')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. Observed customers have their own individual purchasing patterns while the newly arrived customers will assume model population parameters.\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given our assumptions, the value of Zakka Canada is simply the summation of:\n",
    "\n",
    "1. The present value of all new and repeat sales from each month's new cohort of unobserved customers.\n",
    "2. The present value of all gross profit on repeat sales from our observed customers.\n",
    "3. Less the present value of fixed cost for each month.\n",
    "\n",
    "Repeat until infinity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first address (1) by forecasting out our seasonal cohort of new customers.  We fit the time-series model to only the stagnation period (2010 and beyond).\n",
    "$$A_t = C + \\sum_{i=1}^{12}{\\beta_{i}\\delta_{m(t)=i}} $$\n",
    "Where $\\delta_{m(t)}$ is equal to the indicator variable denoting which month $t$ is in.  $C$ is equal to the mean monthly  estimated new arrivals.  Example plot of the arrivals forecast is below but remember that this forecast extends indefinitely as we also believe Zakka Canada can operate to an indefinite horizon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def prepareDesignMatrix(x):\n",
    "    x['M'] = x.index.month\n",
    "    for elem in x['M'].unique():\n",
    "        x['M'+str(elem)] = (x['M'] == elem).astype(int)\n",
    "    x['c'] = 1.\n",
    "    return x.drop('M', axis=1)\n",
    "cohort_over_time = cohort_over_time[cohort_over_time.index >= dt.datetime(2010,2,1)]\n",
    "cohort_over_time = prepareDesignMatrix(cohort_over_time)\n",
    "res = OLS(cohort_over_time['Customer ID'], cohort_over_time[['M1','M2','M3','M4','M5','M6','M7','M8','M9','M10','M11','M12','c']]).fit()\n",
    "print res.summary()\n",
    "cohort_over_time['predicted'] = res.predict()\n",
    "f_cohort = pd.DataFrame(None, index= pd.date_range(cohort_over_time.index.max(), periods=5*12, freq='M'))\n",
    "f_cohort = prepareDesignMatrix(f_cohort)\n",
    "f_cohort['forecast'] = res.predict(f_cohort)\n",
    "cohort_over_time['predicted'].plot(style='--')\n",
    "cohort_over_time['Customer ID'].plot(legend=None)\n",
    "ax = f_cohort['forecast'].plot(style='--', figsize=(12,4))\n",
    "ax.set_xlabel('Date'); ax.set_ylabel('New Customer Cohort Size')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The present value calculation is a bit tricky for both (1) and (2) which we will refer to as Discounted Repeat Value (DRV) and Residual Customer Value (RCV) respectively.  For both these measures, customers will repeatedly buy over multiple periods, these periods are then each discounted back to their present value.  Fader and his colleagues have proposed a closed form calculation for their Pareto/NBD model but not one for BG/NBD.  Through my emails with him and through other authors, they confirm that the actual calculation for BG/NBD model is messy and long so I'm going to approximate it numerically.  We can first start with RCV:\n",
    "\n",
    "The RCV consists of first predicting the expected purchases of each observed customer given their $(x, t_x, T)$.  We can then infer their expected $Z$ from the GG model as shown before to obtain an expected gross profit from them at each month. Note that our periodicity is monthly while the model calculates in days.  Our __residual value__ formula becomes:\n",
    "$$RCV = \\sum_{t=1}^{\\infty}{\\frac{\\sum_{k=1}^{K}{E(Z\\mid p,q,\\gamma,\\bar{z}_{k},x_{k})(E(Y(30t)\\mid x_{k},t_{x,k},T_{k},r,\\alpha,a,b)-E(Y(30(t-1))\\mid x_{k},t_{x,k},T_{k},r,\\alpha,a,b))}}{(1+\\delta)^t}}$$\n",
    "Where $k$ is the k-th customer observed.  See [here](http://www.brucehardie.com/notes/025/gamma_gamma.pdf) and [here](http://mktg.uni-svishtov.bg/ivm/resources/Counting_Your_Customers.pdf) for notation definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def RCV(t, model_trans, d, rfm, zeta):\n",
    "    discount_factor = 1./(1+d)**t\n",
    "    dif = (model_trans.conditional_expected_number_of_purchases_up_to_time(t*30, rfm['frequency'], rfm['recency'], rfm['T'])\n",
    "        - model_trans.conditional_expected_number_of_purchases_up_to_time((t-1)*30, rfm['frequency'], rfm['recency'], rfm['T']))\n",
    "    return ((dif*zeta).sum())*discount_factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DRV (of only ONE customer) is similar in that we discount the expected gross profit of a generic customer from our two models.\n",
    "$$DRV = \\sum_{t=1}^{\\infty}{\\frac{E(Z\\mid p,q,\\gamma)(E(X(30t)\\mid r,\\alpha,a,b)-E(X(30(t-1))\\mid r,\\alpha,a,b))}{(1+\\delta)^t}}$$\n",
    "$$ = \\sum_{t=1}^{\\infty}{\\frac{\\frac{p\\gamma}{q-1}(E(X(30t)\\mid r,\\alpha,a,b)-E(X(30(t-1))\\mid r,\\alpha,a,b))}{(1+\\delta)^t}}$$\n",
    "Since the model only measures repeat transactions, we have to consider the fact that when one cohort of customers arrive, they must first make a purchase before beginning they buy until they die.  For example, if one customer enters at $t=0$ we would assume they would have first made a purchase at that time period and then repeat purchases at $t=1,2,3,...,\\infty$.\n",
    "$$ = \\sum_{t=1}^{\\infty}{\\frac{p\\gamma}{q-1}(\\frac{(E(X(30t)\\mid r,\\alpha,a,b)-E(X(30(t-1))\\mid r,\\alpha,a,b))}{(1+\\delta)^t}+\\delta_{t=1}1)}$$\n",
    "where $\\delta_{t=1}$ is the indicator variable for when if $t=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def DRV(t, model, d, zeta):\n",
    "    discount_factor = 1./(1+d)**t\n",
    "    dif = model.expected_number_of_purchases_up_to_time(t*30) - model.expected_number_of_purchases_up_to_time((t-1)*30)\n",
    "    return zeta*(dif*discount_factor+(t==1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, the DRV is only designed for one customer that arrives at some arbitrary time and buys until infinity.  At each time period $\\tau$, a cohort of size $A_t$ arrives and has a present value of $A_t(DRV)$.  If we allow $\\tau\\rightarrow\\infty$ and take the summation of each present value of cohort DRV, we end up with the __terminal value__ of the firm.\n",
    "$$TV = \\frac{A_{1}DRV}{(1+\\delta)}+\\frac{A_{2}DRV}{(1+\\delta)^2}+...+\\frac{A_{\\infty}DRV}{(1+\\delta)^{\\infty}}$$\n",
    "$$ = \\sum_{\\tau=1}^{\\infty}{\\frac{A_{\\tau}DRV}{(1+\\delta)^{\\tau}}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def TermVal(t, model, d, zeta, ts_model, drv):\n",
    "    q = dt.datetime(2016,1,1)+pd.DateOffset(months=t)\n",
    "    x = pd.DataFrame([[0,0,0,0,0,0,0,0,0,0,0,0]], index=[q], columns=['M1','M2','M3','M4','M5','M6','M7','M8','M9','M10','M11','M12'])\n",
    "    cohort_size = ts_model.predict(prepareDesignMatrix(x))[0]\n",
    "    val = drv * cohort_size/( (1+d)**t )\n",
    "    #print 'present value of cohort size %s at t = %s : %s' % (cohort_size, t,val)\n",
    "    return val\n",
    "\n",
    "def approximate(fn, t=1, eps_tol=1e-6, eps=0, **kwargs):\n",
    "    eps = 0\n",
    "    cf = 0\n",
    "    while True:\n",
    "        cf += fn(t=t, **kwargs)\n",
    "        if(cf - eps < eps_tol):\n",
    "            break\n",
    "        eps = cf; t+=1\n",
    "    return cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "terminal_value = (approximate(TermVal, model=bgf, \n",
    "                              d=mr, zeta=(p*v)/(q-1), \n",
    "                              ts_model=res, \n",
    "                              drv=approximate(DRV, model=bgf, d=mr, zeta=(p*v)/(q-1)),\n",
    "                              eps_tol=1e-8))\n",
    "print 'Discount Repeated Value: $%.2f' % terminal_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "residual_value = approximate(RCV, model_trans=bgf, d=mr, rfm=ret_cust_data, zeta=ez)\n",
    "print 'Residual Customer Value: $%.2f' % residual_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate our firm value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pv_cost = monthly_fc / mr\n",
    "print 'Present Value of Fixed Cost: $%.2f' % pv_cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "firm_value = terminal_value + residual_value - pv_cost\n",
    "print 'Zakka Canada Firm Price: $%.2f' % firm_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Miscellaneous"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Other stuff that's been tried out before.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def prepareDesignMatrix(x):\n",
    "    x['isQ1'] = [1 if d.quarter == 1 else 0 for d in x.index]; x['isQ2'] = [1 if d.quarter == 2 else 0 for d in x.index]; x['isQ3'] = [1 if d.quarter == 3 else 0 for d in x.index]; x['c'] = 1\n",
    "    return x\n",
    "cohort_over_time = cohort_over_time[cohort_over_time.index >= dt.datetime(2010,2,1)]\n",
    "cohort_over_time = prepareDesignMatrix(cohort_over_time)\n",
    "res = OLS(cohort_over_time['Customer ID'], cohort_over_time[['isQ1', 'isQ2','isQ3','c']]).fit()\n",
    "print res.summary()\n",
    "cohort_over_time['predicted'] = res.predict()\n",
    "cohort_over_time[['predicted', 'Customer ID']].plot(legend=None)\n",
    "f_cohort = pd.DataFrame(None, index= pd.date_range(cohort_over_time.index.max(), periods=5*12, freq='M'))\n",
    "f_cohort = prepareDesignMatrix(f_cohort)\n",
    "f_cohort['forecast'] = res.predict(f_cohort)\n",
    "f_cohort['forecast'].plot(style='--', figsize=(12,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "growth_value = 0\n",
    "for t,i in enumerate(f_cohort['forecast']):\n",
    "    growth_value += (i*approximate(MCV, model=bgf, d=mr, zeta=(p*v)/(q-1)))/ ( (1+mr)**(t+1) )\n",
    "print growth_value"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}