Revenue & Cost Forecasting using Python

7.1 Why Forecasting Is the Heart of Valuation

In a DCF valuation, the explicit forecast period — typically 5 years — is where the analyst's judgment has the greatest impact. The terminal value, which we will study in Chapter 13, often contributes 60–80% of the total DCF result. But those first 5 years matter enormously: they establish the base from which terminal growth is projected, and they are the years where company-specific insight (rather than generic perpetuity assumptions) carries the most weight.

Forecasting is not about predicting the future with precision — that is impossible. It is about building reasonable, internally consistent, assumption-transparent projections that can be stress-tested. A good forecast:

• Is grounded in historical data — the trend, volatility, and patterns in past performance inform the range of plausible futures.
• Respects economic logic — revenue cannot grow at 25% forever while the industry grows at 5%. Margins cannot expand indefinitely without attracting competition.
• Separates signal from noise — a single exceptional year (pandemic boom, commodity spike) should not drive a 5-year forecast.
• Is internally consistent — revenue growth, capex, working capital, and margins must be linked. If revenue doubles, working capital must also increase.

7.2 Preparing the Forecasting Environment

Install the forecasting libraries we need:

pip install yfinance pandas numpy matplotlib seaborn statsmodels prophet scikit-learn

Let us start by fetching the historical data we will use throughout this chapter:

import yfinance as yf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')

sns.set_style("darkgrid")
plt.rcParams['figure.figsize'] = (14, 5)

# Fetch historical financials for a sample company
ticker = yf.Ticker("ASIANPAINT.NS")
info = ticker.info

# Extract annual financials
is_df = ticker.financials       # Income Statement
bs_df = ticker.balance_sheet    # Balance Sheet
cf_df = ticker.cashflow         # Cash Flow

# Sort chronologically (oldest first) for time series work
is_df = is_df.sort_index(axis=1)
bs_df = bs_df.sort_index(axis=1)
cf_df = cf_df.sort_index(axis=1)

# Extract key series
def safe_get(df, candidates):
    for c in ([candidates] if isinstance(candidates, str) else candidates):
        if c in df.index:
            return df.loc[c]
    return None

revenue = safe_get(is_df, ['Total Revenue', 'Revenue'])
cogs = safe_get(is_df, ['Cost Of Goods Sold', 'Cost of Revenue'])
ebit = safe_get(is_df, ['Operating Income', 'EBIT'])
net_income = safe_get(is_df, ['Net Income', 'Net Income Common Stockholders'])

# Convert to crores and create a clean DataFrame
years = [d.year for d in revenue.index]  # Extract just the year

df = pd.DataFrame({
    'Year': years,
    'Revenue': revenue.values / 1e7,
    'COGS': cogs.values / 1e7 if cogs is not None else None,
    'EBIT': ebit.values / 1e7 if ebit is not None else None,
    'Net_Income': net_income.values / 1e7 if net_income is not None else None,
})

df = df.set_index('Year').sort_index()
df['Gross_Margin'] = ((df['Revenue'] - df['COGS']) / df['Revenue']) * 100
df['Operating_Margin'] = (df['EBIT'] / df['Revenue']) * 100

print(f"Company: {info.get('longName')}")
print(f"Data range: {df.index.min()} – {df.index.max()}")
print(f"\n{df.round(2)}")

# Quick visualization of the historical revenue trend
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

axes[0].plot(df.index, df['Revenue'], marker='o', linewidth=2, color='#6c8cff', markersize=8)
axes[0].set_title('Revenue (Rs. Crore)', fontweight='bold')
axes[0].set_ylabel('Rs. Crore')
axes[0].grid(True, alpha=0.3)

axes[1].bar(df.index, df['Revenue'].pct_change() * 100, color='#00c9a7', edgecolor='white')
axes[1].axhline(y=0, color='red', linestyle='--', alpha=0.5)
axes[1].set_title('Revenue Growth Rate (%)', fontweight='bold')

axes[2].plot(df.index, df['Gross_Margin'], marker='s', linewidth=2, color='#f0a040', label='Gross Margin')
axes[2].plot(df.index, df['Operating_Margin'], marker='^', linewidth=2, color='#6c8cff', label='Operating Margin')
axes[2].set_title('Margin Trends (%)', fontweight='bold')
axes[2].legend()

plt.suptitle(f"{info.get('longName', 'Company')} — Historical Financial Trends",
             fontsize=14, fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

7.3 Time Series Decomposition: Understanding What Drives the Data

Before you forecast, you must understand the structure of your data. Every time series can be decomposed into three components:

The decomposition can be additive (Y = T + S + R) or multiplicative (Y = T × S × R). Use additive when the seasonal swings are roughly constant in absolute terms. Use multiplicative when seasonal swings grow proportionally with the trend — revenue data is typically multiplicative.

7.3.1 Decomposing Revenue with statsmodels

from statsmodels.tsa.seasonal import seasonal_decompose

# For decomposition, we need quarterly data for meaningful seasonality.
# Let's fetch quarterly data instead of annual.
ticker_q = yf.Ticker("ASIANPAINT.NS")

# Get quarterly revenue (use quarterly financials)
quarterly_is = ticker_q.quarterly_financials.sort_index(axis=1)
q_revenue = safe_get(quarterly_is, ['Total Revenue', 'Revenue'])

# Create a proper time series
q_rev_series = pd.Series(
    q_revenue.values / 1e7,
    index=pd.to_datetime(q_revenue.index),
    name='Quarterly Revenue'
).sort_index()

# Run decomposition (multiplicative for financial data)
decomp = seasonal_decompose(q_rev_series, model='multiplicative', period=4)

# Plot
fig, axes = plt.subplots(4, 1, figsize=(14, 10))

axes[0].plot(q_rev_series.index, q_rev_series.values, color='#6c8cff', linewidth=2)
axes[0].set_title('Observed — Quarterly Revenue (Rs. Crore)', fontweight='bold')

axes[1].plot(decomp.trend.index, decomp.trend.values, color='#00c9a7', linewidth=2)
axes[1].set_title('Trend — The Long-Term Direction', fontweight='bold')

axes[2].plot(decomp.seasonal.index, decomp.seasonal.values, color='#f0a040', linewidth=1.5)
axes[2].set_title('Seasonality — Repeating Quarterly Pattern', fontweight='bold')
axes[2].axhline(y=1.0, color='gray', linestyle='--', alpha=0.5)

axes[3].bar(decomp.resid.index, decomp.resid.values, color='#e0556a' if decomp.resid.mean() < 1 else '#00c9a7', alpha=0.7)
axes[3].set_title('Residual — Unexplained / Random Component', fontweight='bold')
axes[3].axhline(y=1.0, color='gray', linestyle='--', alpha=0.5)

plt.suptitle('Time Series Decomposition — Asian Paints Quarterly Revenue',
             fontsize=14, fontweight='bold', y=1.01)
plt.tight_layout()
plt.show()

7.4 Simple Forecasting Methods: Baselines You Must Beat

Before deploying sophisticated models, establish a baseline. If your complex ARIMA or Prophet model cannot beat a simple moving average, the complex model is not adding value. Every forecasting project should start with these three baselines:

7.4.1 Naive Forecast

The simplest possible model: next year's revenue = this year's revenue. No change.

naive_forecast = df['Revenue'].iloc[-1]  # Last known value
print(f"Naive Forecast (Next Year Revenue): Rs. {naive_forecast:.0f} Cr")
print(f"Implied Growth: 0%")

7.4.2 Moving Average

Forecast = average of the last N years. Smooths out year-to-year noise.

def moving_average_forecast(series, window=3, periods=1):
    """
    Forecast using a trailing moving average.

    Parameters
    ----------
    series : pd.Series — historical values (chronological)
    window : int — number of years to average
    periods : int — number of years to forecast forward
    """
    forecasts = []
    last_value = series.iloc[-window:].mean()
    # For financial forecasting, we typically use the MA of recent growth rates
    # rather than the absolute values, then apply that growth to the last value.
    growth_rates = series.pct_change().dropna()
    avg_growth = growth_rates.iloc[-window:].mean()

    for i in range(1, periods + 1):
        forecasts.append(series.iloc[-1] * (1 + avg_growth) ** i)

    return forecasts

ma_forecast = moving_average_forecast(df['Revenue'], window=3, periods=3)
print("3-Year MA (Growth Rate) Forecast:")
for i, f in enumerate(ma_forecast, 1):
    print(f"  Year +{i}: Rs. {f:.0f} Cr (Growth: {((f/df['Revenue'].iloc[-1])**(1/i)-1)*100:.1f}% CAGR)")

7.4.3 Exponential Smoothing

Exponential smoothing gives more weight to recent observations and less to distant ones. The alpha parameter (0 to 1) controls the decay: alpha near 1 gives almost all weight to the most recent observation; alpha near 0 spreads weight evenly.

from statsmodels.tsa.holtwinters import SimpleExpSmoothing

def exp_smoothing_forecast(series, alpha=0.6, periods=3):
    """Forecast using Simple Exponential Smoothing."""
    model = SimpleExpSmoothing(series).fit(smoothing_level=alpha, optimized=False)
    forecast = model.forecast(periods)
    return forecast

es_forecast = exp_smoothing_forecast(df['Revenue'], alpha=0.6, periods=3)
print("Exponential Smoothing Forecast (alpha=0.6):")
for i, f in enumerate(es_forecast, 1):
    print(f"  Year +{i}: Rs. {f:.0f} Cr")

# Compare baselines
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(df.index, df['Revenue'], 'o-', color='#6c8cff', linewidth=3, markersize=8, label='Historical')
# Extend x-axis for forecasts
forecast_years = [df.index[-1] + i for i in range(1, 4)]
ax.plot(forecast_years, ma_forecast, 's--', color='#00c9a7', linewidth=2, label='Moving Average (3Y Growth)')
ax.plot(forecast_years, es_forecast, '^--', color='#f0a040', linewidth=2, label='Exponential Smoothing')
ax.set_title('Baseline Forecast Comparison', fontweight='bold')
ax.set_ylabel('Revenue (Rs. Crore)')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

7.5 ARIMA Models: The Workhorse of Statistical Forecasting

ARIMA (AutoRegressive Integrated Moving Average) is the most widely used statistical forecasting framework. It models a time series as a function of its own past values (autoregression) and past forecast errors (moving average), after differencing to make the series stationary.

7.5.1 The Three Parameters: p, d, q

7.5.2 Stationarity and the ADF Test

ARIMA requires the time series to be stationary — its statistical properties (mean, variance) must be constant over time. Most financial series are non-stationary (they trend upward), so we difference them.

from statsmodels.tsa.stattools import adfuller

def test_stationarity(series, name="Series"):
    """Augmented Dickey-Fuller test for stationarity."""
    result = adfuller(series.dropna(), autolag='AIC')
    print(f"\n{name} — ADF Test:")
    print(f"  ADF Statistic: {result[0]:.4f}")
    print(f"  p-value:       {result[1]:.4f}")
    print(f"  Critical Values: 1%: {result[4]['1%']:.3f}, 5%: {result[4]['5%']:.3f}")

    if result[1] <= 0.05:
        print(f"  ✓ Series IS stationary (reject H0 at 95% confidence)")
    else:
        print(f"  ✗ Series is NON-stationary — differencing required")

# Test original revenue series
test_stationarity(df['Revenue'], "Revenue (Levels)")

# Test differenced revenue
revenue_diff = df['Revenue'].diff().dropna()
test_stationarity(revenue_diff, "Revenue (First Difference)")

7.5.3 Auto-ARIMA: Let the Algorithm Choose

Manually selecting (p, d, q) by inspecting ACF/PACF plots is educational but time-consuming. The pmdarima library automates this:

pip install pmdarima

from pmdarima import auto_arima

# Fit auto-ARIMA on revenue data
auto_model = auto_arima(
    df['Revenue'],
    start_p=0, start_q=0, max_p=3, max_q=3,
    seasonal=False,          # Annual data — no seasonality
    trace=False,             # Set True to see the search
    error_action='ignore',
    suppress_warnings=True,
    stepwise=True,
    n_fits=50
)

print(f"Best ARIMA model: {auto_model}")
print(f"AIC: {auto_model.aic():.1f}")

# Forecast
arima_forecast, arima_ci = auto_model.predict(n_periods=5, return_conf_int=True)

print("\nARIMA 5-Year Forecast:")
for i, (f, (lo, hi)) in enumerate(zip(arima_forecast, arima_ci), 1):
    growth = ((f / df['Revenue'].iloc[-1]) ** (1/i) - 1) * 100
    print(f"  Year +{i}: Rs. {f:.0f} Cr (CAGR: {growth:.1f}%)  [95% CI: {lo:.0f} – {hi:.0f}]")

# Plot
fig, ax = plt.subplots(figsize=(14, 6))
ax.plot(df.index, df['Revenue'], 'o-', color='#6c8cff', linewidth=3, markersize=8, label='Historical')
forecast_years = [df.index[-1] + i for i in range(1, 6)]
ax.plot(forecast_years, arima_forecast, 's-', color='#00c9a7', linewidth=2, label='ARIMA Forecast')
ax.fill_between(forecast_years, arima_ci[:, 0], arima_ci[:, 1],
                alpha=0.2, color='#00c9a7', label='95% Confidence Interval')
ax.set_title(f'ARIMA({auto_model.order[0]},{auto_model.order[1]},{auto_model.order[2]}) Revenue Forecast',
             fontweight='bold')
ax.set_ylabel('Revenue (Rs. Crore)')
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

7.6 Facebook Prophet: Automated Forecasting with Business Intuition

Prophet (developed by Meta/Facebook) is designed for business forecasting. Unlike ARIMA, which is purely statistical, Prophet is a decomposable additive model with three built-in components that align with how business analysts think about data:

Prophet's key advantage is that it handles the realities of business data: trend changes, missing data, outliers, and seasonal patterns — all with minimal parameter tuning.

7.6.1 Forecasting Revenue with Prophet

from prophet import Prophet

# Prophet requires columns: 'ds' (date) and 'y' (value)
# For annual data, we use the fiscal year-end date
prophet_df = pd.DataFrame({
    'ds': pd.to_datetime([f'{y}-03-31' for y in df.index]),  # March 31 fiscal year-end
    'y': df['Revenue'].values
})

# Initialize and fit
prophet_model = Prophet(
    growth='linear',            # 'logistic' if you want to cap growth
    yearly_seasonality=False,   # Annual data has no within-year seasonality
    changepoint_prior_scale=0.05,  # Flexibility of trend changes (lower = less flexible)
    n_changepoints=len(prophet_df) // 2,  # Potential trend changepoints
    interval_width=0.80         # 80% prediction interval
)

prophet_model.fit(prophet_df)

# Create future DataFrame: 5 years forward
future = prophet_model.make_future_dataframe(periods=5, freq='Y')  # 'Y' = year-end
forecast = prophet_model.predict(future)

# Extract forecasts
prophet_forecast = forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail(5)

print("Prophet 5-Year Revenue Forecast:")
for _, row in prophet_forecast.iterrows():
    year = row['ds'].year
    growth = ((row['yhat'] / df['Revenue'].iloc[-1]) ** (1/5) - 1) * 100  # Approx
    print(f"  FY{year}: Rs. {row['yhat']:.0f} Cr  [80% CI: {row['yhat_lower']:.0f} – {row['yhat_upper']:.0f}]")

# Plot Prophet's components
fig = prophet_model.plot(forecast, figsize=(14, 6))
ax = fig.gca()
ax.set_title('Prophet Revenue Forecast — Asian Paints', fontweight='bold')
ax.set_ylabel('Revenue (Rs. Crore)')
plt.show()

# Plot changepoints (where Prophet detected trend shifts)
fig2 = prophet_model.plot_components(forecast, figsize=(14, 5))
plt.show()

7.6.2 Tuning Prophet with Changepoints

The changepoint_prior_scale is Prophet's most important tuning parameter. It controls how flexible the trend is:

• Low (0.001–0.01): Very rigid — nearly a straight line. Use when the business has stable, predictable growth.
• Medium (0.05–0.1): Default flexibility — captures moderate trend shifts. Use for most companies.
• High (0.5+): Highly flexible — fits almost every data point. Use with caution — this is overfitting territory.

# Compare different changepoint sensitivities
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

for ax, cps, label in zip(axes,
                           [0.01, 0.05, 0.5],
                           ['Rigid (cps=0.01)', 'Balanced (cps=0.05)', 'Flexible (cps=0.5)']):
    m = Prophet(growth='linear', changepoint_prior_scale=cps, yearly_seasonality=False)
    m.fit(prophet_df)
    f = m.predict(m.make_future_dataframe(periods=5, freq='Y'))
    ax.plot(df.index, df['Revenue'], 'o-', color='#6c8cff', linewidth=2, markersize=6, label='Historical')
    ax.plot([r.year for _, r in f.tail(5).iterrows()], f.tail(5)['yhat'], 's--', color='#00c9a7', linewidth=2)
    ax.fill_between([r.year for _, r in f.tail(5).iterrows()],
                    f.tail(5)['yhat_lower'], f.tail(5)['yhat_upper'],
                    alpha=0.2, color='#00c9a7')
    ax.set_title(label, fontweight='bold')
    ax.set_ylabel('Revenue (Rs. Cr)')
    ax.grid(True, alpha=0.3)

plt.suptitle('Prophet Sensitivity: Effect of changepoint_prior_scale on Forecast',
             fontweight='bold', y=1.02)
plt.tight_layout()
plt.show()

7.7 Regression-Based Forecasting: Linking Revenue to Economic Drivers

Time series models (ARIMA, Prophet) forecast revenue using only its own history. But revenue does not exist in a vacuum — it is driven by macroeconomic factors and company-specific investments. Regression forecasting models revenue as a function of these drivers.

7.7.1 Simple Linear Regression: Revenue vs Time

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_absolute_error, mean_squared_error

# Revenue as a function of time (year index)
X = np.arange(len(df)).reshape(-1, 1)  # Year indices: 0, 1, 2, ...
y = df['Revenue'].values

model_lr = LinearRegression()
model_lr.fit(X, y)

# Forecast next 5 years
X_future = np.arange(len(df), len(df) + 5).reshape(-1, 1)
lr_forecast = model_lr.predict(X_future)

print(f"Linear Regression: Revenue = {model_lr.intercept_:.0f} + {model_lr.coef_[0]:.0f} × YearIndex")
print(f"R² = {model_lr.score(X, y):.3f}")
print("\n5-Year Forecast:")
for i, f in enumerate(lr_forecast, 1):
    growth = ((f / y[-1]) ** (1/i) - 1) * 100
    print(f"  Year +{i}: Rs. {f:.0f} Cr (CAGR: {growth:.1f}%)")

7.7.2 Multiple Regression: Revenue vs Multiple Drivers

For a richer model, we could regress revenue against:

• GDP growth — for cyclical businesses (autos, cement, steel)
• Past capex — capacity expansion drives future revenue
• Store count / capacity — for retail, restaurants, telecom
• Industry revenue — to capture sector-wide trends

# Example: Revenue as a function of Time + Lagged Capex
# (Simplified — real GDP data would come from RBI/World Bank APIs)

# Fetch capex data
capex = safe_get(cf_df, ['Capital Expenditure', 'Capital Expenditure Reported'])
if capex is not None:
    capex_values = capex.values / 1e7  # Convert to crores
    # Lag capex by 1 year (this year's capex enables next year's revenue)
    capex_lagged = np.roll(capex_values, 1)
    capex_lagged[0] = capex_values[0]  # Fill the first year with unlagged value

    # Build feature matrix
    X_multi = np.column_stack([
        np.arange(len(df)),    # Time trend
        capex_lagged            # Lagged capex
    ])

    model_multi = LinearRegression()
    model_multi.fit(X_multi, y)

    print(f"Multiple Regression R²: {model_multi.score(X_multi, y):.3f}")
    print(f"Coefficients: Time={model_multi.coef_[0]:.0f}, LaggedCapex={model_multi.coef_[1]:.1f}")

    # To forecast, you need future capex assumptions — this links revenue to your investment forecast

7.8 Forecasting Costs and Margins

Revenue is only half the equation. To get to free cash flow, you must forecast costs. The approach depends on the cost category:

7.8.1 Cost of Goods Sold (COGS): Margin-Based Forecasting

The most common approach is to forecast gross margin rather than COGS directly. Margins are more stable than absolute costs and reflect the company's competitive position.

# Forecast gross margin using its historical average and trend
historical_gm = df['Gross_Margin'].dropna()
avg_gm = historical_gm.mean()
trend_gm = np.polyfit(range(len(historical_gm)), historical_gm.values, 1)[0]

print(f"Historical Gross Margin: Mean = {avg_gm:.1f}%, Trend = {trend_gm:.1f}pp/year")

# Three margin scenarios for forecasting
# 1. Constant: keep margin at the historical average
# 2. Trend: continue the historical trend (positive or negative)
# 3. Mean-reverting: margins revert to industry average over time

def forecast_cogs_from_margin(revenue_forecast, gross_margin):
    """Given a revenue forecast and a gross margin, compute COGS."""
    cogs = revenue_forecast * (1 - gross_margin / 100)
    return cogs

# Example: Use Prophet revenue forecast with three margin scenarios
rev_forecast_values = prophet_forecast['yhat'].values[:5]

scenarios = {
    'Constant Margin': avg_gm,
    'Trend Margin': avg_gm + trend_gm * len(df),  # Extrapolate trend
    'Industry Avg Margin': 38.0  # Hypothetical industry average
}

for scenario, margin in scenarios.items():
    cogs = forecast_cogs_from_margin(rev_forecast_values, margin)
    print(f"\n{scenario} ({margin:.1f}%):")
    for i, (r, c) in enumerate(zip(rev_forecast_values, cogs), 1):
        print(f"  Year +{i}: Revenue Rs.{r:.0f}Cr, COGS Rs.{c:.0f}Cr, "
              f"Gross Profit Rs.{r-c:.0f}Cr")

7.8.2 Operating Expenses: Fixed vs Variable Split

Operating expenses are a mix of fixed costs (rent, base salaries, depreciation) and variable costs (commissions, freight, marketing). The split determines how margins evolve as revenue grows:

• High fixed-cost businesses (IT services, telecom): Margins expand as revenue grows because fixed costs are spread over more revenue — operating leverage.
• High variable-cost businesses (trading, commodity manufacturing): Margins are stable because costs scale proportionally with revenue.

def forecast_opex(revenue_forecast, fixed_opex_pct=0.40, variable_opex_pct=0.15):
    """
    Forecast operating expenses assuming a fixed/variable split.

    Parameters
    ----------
    revenue_forecast : array of future revenue values
    fixed_opex_pct : float — fixed opex as % of CURRENT (last year) revenue
    variable_opex_pct : float — variable opex as % of revenue
    """
    last_revenue = df['Revenue'].iloc[-1]
    fixed_opex = last_revenue * fixed_opex_pct

    opex_forecast = []
    for rev in revenue_forecast:
        total_opex = fixed_opex + (rev * variable_opex_pct)
        opex_forecast.append(total_opex)

    return np.array(opex_forecast)

# Example
rev_future = prophet_forecast['yhat'].values[:5]
opex = forecast_opex(rev_future, fixed_opex_pct=0.12, variable_opex_pct=0.10)

last_op_margin = df['Operating_Margin'].iloc[-1]
future_op_margin = ((rev_future - opex) / rev_future) * 100

print("Operating Margin Evolution (with Operating Leverage):")
for i, (r, o, m) in enumerate(zip(rev_future, opex, future_op_margin), 1):
    print(f"  Year +{i}: Revenue Rs.{r:.0f}Cr, Opex Rs.{o:.0f}Cr, "
          f"Op Margin {m:.1f}% (Δ from last: {m-last_op_margin:+.1f}pp)")

7.8.3 Depreciation: A Function of Capex History

Depreciation is forecast based on the existing asset base plus future capex. A simple approach: depreciation = (historical depreciation rate) × (net PPE + new capex).

def forecast_depreciation(balance_sheet, income_stmt, future_capex_assumptions, depr_rate=None):
    """
    Forecast depreciation based on asset base and future capex.

    future_capex_assumptions : list — projected capex for each forecast year
    depr_rate : float — depreciation as % of net PPE (if None, uses historical avg)
    """
    net_ppe = safe_get(balance_sheet, ['Net PPE', 'Property Plant Equipment Net'])
    d_a = safe_get(income_stmt, ['Depreciation And Amortization',
                                  'Depreciation Amortization Depletion'])

    if net_ppe is None or d_a is None:
        print("Cannot forecast depreciation: missing data")
        return None

    if depr_rate is None:
        depr_rate = (d_a / net_ppe).mean()

    depr_forecast = []
    current_ppe = net_ppe.iloc[-1]

    for capex in future_capex_assumptions:
        # Simplified: D&A = rate × (beginning PPE + capex)
        depr = depr_rate * (current_ppe + capex)
        depr_forecast.append(depr)
        current_ppe = current_ppe + capex - depr  # Net PPE evolves

    print(f"Historical Depreciation Rate: {depr_rate*100:.1f}% of Net PPE")
    print("\nDepreciation Forecast:")
    for i, (c, d) in enumerate(zip(future_capex_assumptions, depr_forecast), 1):
        print(f"  Year +{i}: Capex Rs.{c:.0f}Cr → Depreciation Rs.{d:.0f}Cr")

    return np.array(depr_forecast)

# Example
depr_forecast = forecast_depreciation(
    bs_df, is_df,
    future_capex_assumptions=[1200, 1300, 1400, 1500, 1600]  # 5 years of assumed capex
)

7.9 Evaluating Forecast Accuracy: Choosing the Best Model

How do you know which forecasting model to trust? You test them against data they have not seen. The standard approach is backtesting: hold out the last 2–3 years of data, train models on the earlier years, forecast the held-out period, and measure the error.

7.9.1 Error Metrics

7.9.2 Backtesting Framework

from sklearn.metrics import mean_absolute_error, mean_squared_error

def backtest_models(series, test_years=3):
    """
    Backtest multiple forecasting models against held-out data.

    Parameters
    ----------
    series : pd.Series — historical values (chronological)
    test_years : int — number of years to hold out for testing

    Returns
    -------
    DataFrame of error metrics for each model.
    """
    train = series.iloc[:-test_years]
    test = series.iloc[-test_years:]

    results = {}

    # 1. Naive (last value repeated)
    naive_preds = np.repeat(train.iloc[-1], test_years)
    results['Naive'] = {
        'MAE': mean_absolute_error(test, naive_preds),
        'RMSE': np.sqrt(mean_squared_error(test, naive_preds)),
        'MAPE': np.mean(np.abs((test.values - naive_preds) / test.values)) * 100
    }

    # 2. Moving Average (3-year growth rate)
    growth = train.pct_change().dropna().iloc[-3:].mean()
    ma_preds = np.array([train.iloc[-1] * (1 + growth) ** i for i in range(1, test_years+1)])
    results['Moving Avg (3Y)'] = {
        'MAE': mean_absolute_error(test, ma_preds),
        'RMSE': np.sqrt(mean_squared_error(test, ma_preds)),
        'MAPE': np.mean(np.abs((test.values - ma_preds) / test.values)) * 100
    }

    # 3. Exponential Smoothing
    try:
        es_model = SimpleExpSmoothing(train).fit(smoothing_level=0.6, optimized=False)
        es_preds = es_model.forecast(test_years)
        results['Exp Smoothing'] = {
            'MAE': mean_absolute_error(test, es_preds),
            'RMSE': np.sqrt(mean_squared_error(test, es_preds)),
            'MAPE': np.mean(np.abs((test.values - es_preds) / test.values)) * 100
        }
    except:
        pass

    # 4. Linear Regression
    X_train = np.arange(len(train)).reshape(-1, 1)
    X_test = np.arange(len(train), len(train) + test_years).reshape(-1, 1)
    lr = LinearRegression().fit(X_train, train.values)
    lr_preds = lr.predict(X_test)
    results['Linear Regression'] = {
        'MAE': mean_absolute_error(test, lr_preds),
        'RMSE': np.sqrt(mean_squared_error(test, lr_preds)),
        'MAPE': np.mean(np.abs((test.values - lr_preds) / test.values)) * 100
    }

    # 5. ARIMA (if enough data)
    if len(train) >= 6:
        try:
            arima_m = auto_arima(train, seasonal=False, suppress_warnings=True,
                                 error_action='ignore', stepwise=True)
            arima_preds = arima_m.predict(n_periods=test_years)
            results['ARIMA'] = {
                'MAE': mean_absolute_error(test, arima_preds),
                'RMSE': np.sqrt(mean_squared_error(test, arima_preds)),
                'MAPE': np.mean(np.abs((test.values - arima_preds) / test.values)) * 100
            }
        except:
            pass

    # 6. Prophet
    try:
        prophet_train = pd.DataFrame({
            'ds': pd.to_datetime([f'{y}-03-31' for y in train.index]),
            'y': train.values
        })
        p_model = Prophet(yearly_seasonality=False)
        p_model.fit(prophet_train)
        p_future = p_model.make_future_dataframe(periods=test_years, freq='Y')
        p_forecast = p_model.predict(p_future)
        prophet_preds = p_forecast.tail(test_years)['yhat'].values
        results['Prophet'] = {
            'MAE': mean_absolute_error(test, prophet_preds),
            'RMSE': np.sqrt(mean_squared_error(test, prophet_preds)),
            'MAPE': np.mean(np.abs((test.values - prophet_preds) / test.values)) * 100
        }
    except:
        pass

    return pd.DataFrame(results).T.sort_values('RMSE')

# Run backtest
backtest_results = backtest_models(df['Revenue'], test_years=3)
print("=== Backtest Results: Revenue Forecasting Models ===")
print(f"(Lower RMSE = better. Models ranked by RMSE.)\n")
print(backtest_results.round(2))

# Highlight the winner
best_model = backtest_results.index[0]
print(f"\nBest model: {best_model} (RMSE = {backtest_results.loc[best_model, 'RMSE']:.1f})")

7.10 The Complete Forecasting Pipeline

Now we combine everything into a single, reusable forecasting pipeline that takes a ticker symbol and produces a 5-year revenue and margin forecast:

from prophet import Prophet
from pmdarima import auto_arima
from sklearn.linear_model import LinearRegression

def complete_forecast_pipeline(ticker_symbol, forecast_years=5):
    """
    End-to-end revenue and margin forecasting pipeline.

    Returns a dict with:
      - historical: DataFrame of historical revenue and margins
      - prophet_forecast: Prophet model and forecast for revenue
      - margin_scenarios: 3 margin scenarios (constant, trend, industry)
      - recommended: recommended forecast for DCF use
    """
    ticker = yf.Ticker(ticker_symbol)
    info = ticker.info
    is_df = ticker.financials.sort_index(axis=1)
    bs_df = ticker.balance_sheet.sort_index(axis=1)

    # --- Extract Historical Data ---
    def safe_get(df, candidates):
        for c in ([candidates] if isinstance(candidates, str) else candidates):
            if c in df.index:
                return df.loc[c]
        return None

    rev = safe_get(is_df, ['Total Revenue', 'Revenue']) / 1e7
    cogs = safe_get(is_df, ['Cost Of Goods Sold', 'Cost of Revenue'])
    cogs = cogs / 1e7 if cogs is not None else None
    ebit = safe_get(is_df, ['Operating Income', 'EBIT'])
    ebit = ebit / 1e7 if ebit is not None else None
    ni = safe_get(is_df, ['Net Income', 'Net Income Common Stockholders'])
    ni = ni / 1e7 if ni is not None else None

    years = [d.year for d in rev.index]

    historical = pd.DataFrame({
        'Year': years,
        'Revenue': rev.values,
    })
    if cogs is not None:
        historical['COGS'] = cogs.values
        historical['Gross_Margin'] = (historical['Revenue'] - historical['COGS']) / historical['Revenue'] * 100
    if ebit is not None:
        historical['EBIT'] = ebit.values
        historical['Operating_Margin'] = historical['EBIT'] / historical['Revenue'] * 100
    if ni is not None:
        historical['Net_Income'] = ni.values

    historical = historical.set_index('Year').sort_index()

    # --- Forecast Revenue with Prophet ---
    prophet_df = pd.DataFrame({
        'ds': pd.to_datetime([f'{y}-03-31' for y in historical.index]),
        'y': historical['Revenue'].values
    })

    p_model = Prophet(growth='linear', changepoint_prior_scale=0.05,
                      yearly_seasonality=False, interval_width=0.80)
    p_model.fit(prophet_df)
    p_future = p_model.make_future_dataframe(periods=forecast_years, freq='Y')
    p_forecast = p_model.predict(p_future)

    rev_forecast = p_forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail(forecast_years)
    rev_forecast['Year'] = rev_forecast['ds'].dt.year
    rev_forecast = rev_forecast.set_index('Year')

    # --- Margin Scenarios ---
    margin_scenarios = {}
    if 'Gross_Margin' in historical.columns:
        gm_mean = historical['Gross_Margin'].mean()
        gm_trend = np.polyfit(range(len(historical)), historical['Gross_Margin'].values, 1)[0]
        margin_scenarios['Gross Margin'] = {
            'Historical Avg': gm_mean,
            'Recent Trend': historical['Gross_Margin'].iloc[-1],
            'Trend Extrapolation': historical['Gross_Margin'].iloc[-1] + gm_trend * forecast_years
        }

    if 'Operating_Margin' in historical.columns:
        om_mean = historical['Operating_Margin'].mean()
        om_trend = np.polyfit(range(len(historical)), historical['Operating_Margin'].values, 1)[0]
        margin_scenarios['Operating Margin'] = {
            'Historical Avg': om_mean,
            'Recent Trend': historical['Operating_Margin'].iloc[-1],
            'Trend Extrapolation': historical['Operating_Margin'].iloc[-1] + om_trend * forecast_years
        }

    # --- Print Summary ---
    print(f"\n{'='*60}")
    print(f"  FORECAST SUMMARY: {info.get('longName', ticker_symbol)}")
    print(f"{'='*60}")
    print(f"\n  Revenue Forecast (Prophet):")
    for year, row in rev_forecast.iterrows():
        growth = ((row['yhat'] / historical['Revenue'].iloc[-1]) ** (1/(year-historical.index[-1])) - 1) * 100
        print(f"    FY{year}: Rs. {row['yhat']:.0f} Cr [80% CI: {row['yhat_lower']:.0f} – {row['yhat_upper']:.0f}] "
              f"→ CAGR: {growth:.1f}%")

    print(f"\n  Margin Scenarios:")
    for margin_type, scenarios in margin_scenarios.items():
        print(f"    {margin_type}:")
        for scenario_name, value in scenarios.items():
            print(f"      {scenario_name}: {value:.1f}%")

    return {
        'historical': historical,
        'prophet_model': p_model,
        'prophet_forecast_full': p_forecast,
        'revenue_forecast': rev_forecast,
        'margin_scenarios': margin_scenarios,
        'company_info': info,
        'recommendation': {
            'use_revenue_model': 'Prophet',
            'use_margin': 'Historical Average (conservative) unless trend is strongly supported by competitive analysis'
        }
    }


# --- Run the complete pipeline ---
result = complete_forecast_pipeline("ASIANPAINT.NS")

# Access components for further analysis
historical = result['historical']
rev_forecast = result['revenue_forecast']
margins = result['margin_scenarios']