Forecasting Customer Lifetime Value Using RFM-Analysis and Markov Chain

# Import required libraries
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import squarify
from tabulate import tabulate

# Set up environment
pd.options.display.float_format = '{:,.2f}'.format
plt.rcParams["figure.figsize"] = (12, 8)

# Load the dataset
columns = ['customer_id', 'purchase_amount', 'date_of_purchase']
df = pd.read_csv('purchases.txt', header=None, sep='\t', names=columns)

# Convert to datetime and create time-based features
df['date_of_purchase'] = pd.to_datetime(df['date_of_purchase'], format='%Y-%m-%d')
df['year_of_purchase'] = df['date_of_purchase'].dt.year

# Set reference date for calculating recency
basedate = pd.Timestamp('2016-01-01')
df['days_since'] = (basedate - df['date_of_purchase']).dt.days

def calculate_rfm(dataframe, reference_date=None, days_lookback=None):
    """
    Calculate RFM metrics for customers.
    
    Parameters:
    -----------
    dataframe : DataFrame
        Transaction data with customer_id, purchase_amount, days_since
    days_lookback : int, optional
        Only consider transactions within this many days (default: None, uses all)
    
    Returns:
    --------
    DataFrame with customer_id and RFM metrics
    """
    if days_lookback is not None:
        df_filtered = dataframe[dataframe['days_since'] > days_lookback].copy()
        df_filtered['days_since'] = df_filtered['days_since'] - days_lookback
    else:
        df_filtered = dataframe.copy()
    
    rfm = df_filtered.groupby('customer_id').agg({
        'days_since': ['min', 'max'],
        'customer_id': 'count',
        'purchase_amount': ['mean', 'max']
    })
    
    rfm.columns = ['recency', 'first_purchase', 'frequency', 'avg_amount', 'max_amount']
    rfm = rfm.reset_index()
    
    return rfm


def segment_customers(rfm_data):
    """
    Segment customers based on RFM metrics using managerial rules.
    """
    customers = rfm_data.copy()
    
    conditions = [
        customers['recency'] > 365 * 3,
        (customers['recency'] <= 365 * 3) & (customers['recency'] > 365 * 2),
        (customers['recency'] <= 365 * 2) & (customers['recency'] > 365) & 
            (customers['first_purchase'] <= 365 * 2),
        (customers['recency'] <= 365 * 2) & (customers['recency'] > 365) & 
            (customers['avg_amount'] >= 100),
        (customers['recency'] <= 365 * 2) & (customers['recency'] > 365) & 
            (customers['avg_amount'] < 100),
        (customers['recency'] <= 365) & (customers['first_purchase'] <= 365),
        (customers['recency'] <= 365) & (customers['avg_amount'] >= 100),
        (customers['recency'] <= 365) & (customers['avg_amount'] < 100)
    ]
    
    choices = [
        'inactive',
        'cold',
        'new warm',
        'warm high value',
        'warm low value',
        'new active',
        'active high value',
        'active low value'
    ]
    
    customers['segment'] = np.select(conditions, choices, default='other')
    
    segment_order = [
        'inactive', 'cold', 'warm high value', 'warm low value', 'new warm',
        'active high value', 'active low value', 'new active'
    ]
    customers['segment'] = pd.Categorical(
        customers['segment'],
        categories=segment_order,
        ordered=True
    )
    
    return customers.sort_values('segment')


# Calculate segments for both years
customers_2014 = calculate_rfm(df, days_lookback=365)
customers_2014 = segment_customers(customers_2014)

customers_2015 = calculate_rfm(df)
customers_2015 = segment_customers(customers_2015)

print("2014 Segment Distribution:")
print(customers_2014['segment'].value_counts())
print("\n2015 Segment Distribution:")
print(customers_2015['segment'].value_counts())

# Merge 2014 and 2015 segments by customer_id
# Use left join to keep all 2014 customers
transitions_df = customers_2014[['customer_id', 'segment']].merge(
    customers_2015[['customer_id', 'segment']], 
    on='customer_id', 
    how='left',
    suffixes=('_2014', '_2015')
)

# Handle customers who didn't appear in 2015 data (extremely rare edge case)
# They would have transitioned to inactive
transitions_df['segment_2015'] = transitions_df['segment_2015'].fillna('inactive')

print("\nSample transitions:")
print(tabulate(transitions_df.tail(10), headers='keys', tablefmt='psql', showindex=False))

# Create cross-tabulation of transitions (raw counts)
transition_counts = pd.crosstab(
    transitions_df['segment_2014'], 
    transitions_df['segment_2015'],
    dropna=False
)

print("\nTransition Counts (2014 → 2015):")
print(transition_counts)

# Convert counts to probabilities
transition_matrix = transition_counts.div(transition_counts.sum(axis=1), axis=0)

print("\nTransition Matrix (Probabilities):")
print(transition_matrix.round(3))

# Get 2014 segment counts as a vector
segment_2014_counts = customers_2014['segment'].value_counts().reindex(
    transition_matrix.index, fill_value=0
)

# Predict 2015 distribution by matrix multiplication
predicted_2015 = segment_2014_counts.dot(transition_matrix).round(0)

# Get actual 2015 distribution
actual_2015 = customers_2015['segment'].value_counts().reindex(
    transition_matrix.columns, fill_value=0
)

# Compare predictions to actuals
validation = pd.DataFrame({
    '2014_actual': segment_2014_counts,
    '2015_predicted': predicted_2015.astype(int),
    '2015_actual': actual_2015,
    'error': actual_2015 - predicted_2015.astype(int),
    'pct_error': ((actual_2015 - predicted_2015) / actual_2015 * 100).round(1)
})

print("\nTransition Matrix Validation:")
print(tabulate(validation, headers='keys', tablefmt='psql'))

# Initialize projection matrix
# Rows = segments, Columns = years
years = np.arange(2015, 2026)
segment_projection = pd.DataFrame(
    0, 
    index=customers_2015['segment'].cat.categories,
    columns=years
)

# Populate 2015 with actual counts
segment_projection[2015] = customers_2015['segment'].value_counts().reindex(
    segment_projection.index, fill_value=0
)

# Project forward using matrix multiplication
for year in range(2016, 2026):
    segment_projection[year] = segment_projection[year-1].dot(transition_matrix).round(0)
    
# Convert to integers (can't have fractional customers)
segment_projection = segment_projection.astype(int)

print("\nSegment Projection (2015-2025):")
print(segment_projection)

# Plot inactive customer trajectory
fig, ax = plt.subplots(figsize=(12, 6))
segment_projection.loc['inactive'].plot(kind='bar', ax=ax, color='#9b59b6', alpha=0.7)
ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Number of Customers', fontsize=12)
ax.set_title('Projected Inactive Customer Growth (No New Acquisition)', fontsize=14)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig('inactive_projection.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate average 2015 revenue by segment
revenue_2015 = df[df['year_of_purchase'] == 2015].groupby('customer_id')['purchase_amount'].sum()
customers_with_revenue = customers_2015.merge(
    revenue_2015.rename('revenue_2015'), 
    on='customer_id', 
    how='left'
)
customers_with_revenue['revenue_2015'] = customers_with_revenue['revenue_2015'].fillna(0)

segment_revenue = customers_with_revenue.groupby('segment', observed=True)['revenue_2015'].mean()

print("\nAverage 2015 Revenue by Segment:")
print(segment_revenue.round(2))

# Create revenue vector (aligned with segment order in our projection matrix)
revenue_vector = segment_revenue.reindex(segment_projection.index, fill_value=0)

# Calculate revenue projection for each year
revenue_projection = segment_projection.multiply(revenue_vector, axis=0)

# Sum across segments to get total yearly revenue
yearly_revenue = revenue_projection.sum(axis=0)

print("\nProjected Yearly Revenue (2015-2025):")
print(yearly_revenue.round(0))

# Plot revenue projection
fig, ax = plt.subplots(figsize=(12, 6))
yearly_revenue.plot(kind='bar', ax=ax, color='#3498db', alpha=0.7)
ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Revenue ($)', fontsize=12)
ax.set_title('Projected Annual Revenue (No Acquisition, No Discounting)', fontsize=14)
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, p: f'${x/1000:.0f}K'))
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig('revenue_projection.png', dpi=150, bbox_inches='tight')
plt.show()

# Set discount rate and calculate discount factors
discount_rate = 0.10
discount_factors = 1 / ((1 + discount_rate) ** np.arange(0, 11))

print("\nDiscount Factors (10% rate):")
discount_df = pd.DataFrame({
    'year': years,
    'discount_factor': discount_factors,
    'dollar_value': discount_factors
}).round(3)
print(tabulate(discount_df, headers='keys', tablefmt='psql', showindex=False))

# Calculate discounted yearly revenue
discounted_revenue = yearly_revenue * discount_factors

print("\nRevenue Comparison (Nominal vs. Present Value):")
comparison = pd.DataFrame({
    'year': years,
    'nominal_revenue': yearly_revenue,
    'discount_factor': discount_factors,
    'present_value': discounted_revenue
}).round(0)
print(tabulate(comparison, headers='keys', tablefmt='psql', showindex=False))

# Plot nominal vs. discounted revenue
fig, ax = plt.subplots(figsize=(12, 6))

x = np.arange(len(years))
width = 0.35

bars1 = ax.bar(x - width/2, yearly_revenue, width, label='Nominal Revenue', 
               color='#3498db', alpha=0.7)
bars2 = ax.bar(x + width/2, discounted_revenue, width, label='Present Value', 
               color='#e74c3c', alpha=0.7)

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Revenue ($)', fontsize=12)
ax.set_title('Nominal vs. Discounted Revenue Projections', fontsize=14)
ax.set_xticks(x)
ax.set_xticklabels(years)
ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, p: f'${x/1000:.0f}K'))
ax.legend()
ax.grid(axis='y', alpha=0.3)

plt.tight_layout()
plt.savefig('revenue_discounted_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate cumulative metrics
cumulative_nominal = yearly_revenue.cumsum()
cumulative_pv = discounted_revenue.cumsum()

# Calculate total CLV (sum of all discounted cash flows)
total_clv = discounted_revenue.sum()
total_customers = len(customers_2015)
clv_per_customer = total_clv / total_customers

print(f"\nCustomer Lifetime Value Analysis:")
print(f"Total customer base (2015): {total_customers:,}")
print(f"Total CLV (10-year, discounted): ${total_clv:,.2f}")
print(f"Average CLV per customer: ${clv_per_customer:,.2f}")
print(f"\nComparison:")
print(f"Total nominal revenue (10-year): ${yearly_revenue.sum():,.2f}")
print(f"Discount impact: {(1 - total_clv/yearly_revenue.sum())*100:.1f}% reduction")

# Plot cumulative CLV over time
fig, ax = plt.subplots(figsize=(12, 6))

ax.plot(years, cumulative_pv/1000, marker='o', linewidth=2.5, 
        markersize=8, color='#2ecc71', label='Cumulative CLV (Present Value)')
ax.fill_between(years, 0, cumulative_pv/1000, alpha=0.2, color='#2ecc71')

ax.set_xlabel('Year', fontsize=12)
ax.set_ylabel('Cumulative Value ($K)', fontsize=12)
ax.set_title('Cumulative Customer Lifetime Value (10-year, Discounted)', fontsize=14)
ax.grid(alpha=0.3)
ax.legend(fontsize=11)

# Add annotation for total CLV
ax.annotate(f'Total CLV: ${total_clv/1000:.0f}K', 
            xy=(2025, cumulative_pv.iloc[-1]/1000),
            xytext=(2022, cumulative_pv.iloc[-1]/1000 + 200),
            fontsize=11,
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5),
            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.3'))

plt.tight_layout()
plt.savefig('cumulative_clv.png', dpi=150, bbox_inches='tight')
plt.show()

# Calculate segment-specific CLV
segment_clv = revenue_projection.multiply(discount_factors, axis=1).sum(axis=1) / segment_projection[2015]
segment_clv = segment_clv.replace([np.inf, -np.inf], 0)  # Handle divide by zero

# Create comprehensive segment analysis
segment_analysis = pd.DataFrame({
    '2015_count': segment_projection[2015],
    'avg_2015_revenue': revenue_vector,
    'total_segment_clv': revenue_projection.multiply(discount_factors, axis=1).sum(axis=1),
    'clv_per_customer': segment_clv
}).round(2)

segment_analysis['pct_of_total_clv'] = (
    segment_analysis['total_segment_clv'] / segment_analysis['total_segment_clv'].sum() * 100
).round(1)

print("\nSegment-Level CLV Analysis:")
print(tabulate(segment_analysis, headers='keys', tablefmt='psql'))