"""
Real Options Valuation Template
================================
A reusable class for valuing real options using both
Binomial Tree and Black-Scholes methods.

Usage:
    rov = RealOptionsValuer(S=400, K=500, T=2, sigma=0.35, r=0.07)
    
    # Option to Delay (Call)
    value_bs = rov.black_scholes(option_type='call', dividend_yield=0.04)
    value_bt = rov.binomial_tree(option_type='call', n_steps=50)
    
    # Option to Abandon (Put)
    value_abandon = rov.black_scholes(option_type='put')
    
    # Sensitivity analysis
    rov.sensitivity(parameter='sigma', range=(0.1, 0.8))
"""

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt


class RealOptionsValuer:
    """Real Options Valuation using Black-Scholes and Binomial Tree."""
    
    def __init__(self, S: float, K: float, T: float, sigma: float, r: float = 0.07):
        """
        Parameters:
            S: PV of project's future cash flows (underlying asset value)
            K: Investment cost / strike price
            T: Time to decision deadline (years)
            sigma: Volatility of project value
            r: Risk-free rate
        """
        self.S = S
        self.K = K
        self.T = T
        self.sigma = sigma
        self.r = r
    
    def black_scholes(self, option_type: str = 'call', dividend_yield: float = 0.0) -> dict:
        """
        Black-Scholes-Merton valuation for European-style real options.
        
        Args:
            option_type: 'call' (delay/expand) or 'put' (abandon)
            dividend_yield: Annual cash flows lost by waiting (q)
        
        Returns dict with value, d1, d2, and all intermediate calculations.
        """
        S, K, T, sigma, r, q = self.S, self.K, self.T, self.sigma, self.r, dividend_yield
        
        d1 = (np.log(S / K) + (r - q + sigma**2 / 2) * T) / (sigma * np.sqrt(T))
        d2 = d1 - sigma * np.sqrt(T)
        
        if option_type == 'call':
            value = S * np.exp(-q * T) * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
            name = "Option to Delay / Expand (Call)"
        else:
            value = K * np.exp(-r * T) * norm.cdf(-d2) - S * np.exp(-q * T) * norm.cdf(-d1)
            name = "Option to Abandon (Put)"
        
        npv = S - K if option_type == 'call' else K - S
        
        result = {
            'option_type': option_type,
            'option_name': name,
            'value': value,
            'd1': d1, 'd2': d2,
            'N_d1': norm.cdf(d1), 'N_d2': norm.cdf(d2),
            'npv': npv,
            'enpv': npv + value if option_type == 'call' else self.S + value,
            'premium_over_npv': value - max(npv, 0)
        }
        
        self._print_result(result)
        return result
    
    def binomial_tree(self, option_type: str = 'call', n_steps: int = 50,
                      dividend_yield: float = 0.0, american: bool = True) -> dict:
        """
        Binomial tree valuation. Supports both American and European exercise.
        
        Args:
            option_type: 'call' or 'put'
            n_steps: Number of time steps
            dividend_yield: Annual dividend yield
            american: If True, allows early exercise at any node
        
        Returns dict with option value and tree parameters.
        """
        S, K, T, sigma, r, q = self.S, self.K, self.T, self.sigma, self.r, dividend_yield
        n = n_steps
        dt = T / n
        
        u = np.exp(sigma * np.sqrt(dt))
        d = 1 / u
        p = (np.exp((r - q) * dt) - d) / (u - d)
        
        # Build terminal values
        opt_vals = np.zeros(n + 1)
        for i in range(n + 1):
            ST = S * (u ** (n - i)) * (d ** i)
            if option_type == 'call':
                opt_vals[i] = max(ST - K, 0)
            else:
                opt_vals[i] = max(K - ST, 0)
        
        # Backward induction
        for j in range(n - 1, -1, -1):
            for i in range(j + 1):
                cont = np.exp(-r * dt) * (p * opt_vals[i] + (1 - p) * opt_vals[i + 1])
                if american:
                    Sj = S * (u ** (j - i)) * (d ** i)
                    if option_type == 'call':
                        exercise = max(Sj - K, 0)
                    else:
                        exercise = max(K - Sj, 0)
                    opt_vals[i] = max(cont, exercise)
                else:
                    opt_vals[i] = cont
        
        value = opt_vals[0]
        exercise_style = "American" if american else "European"
        name = f"{'Option to Delay/Expand (Call)' if option_type == 'call' else 'Option to Abandon (Put)'}"
        
        result = {
            'option_type': option_type,
            'exercise': exercise_style,
            'option_name': name,
            'value': value,
            'n_steps': n_steps,
            'u': u, 'd': d, 'p': p,
            'npv': S - K if option_type == 'call' else K - S
        }
        
        self._print_result(result)
        return result
    
    def sensitivity(self, parameter: str = 'sigma', 
                    value_range: tuple = (0.1, 0.8),
                    n_points: int = 20,
                    option_type: str = 'call',
                    dividend_yield: float = 0.0,
                    method: str = 'black_scholes') -> None:
        """
        Run sensitivity analysis: vary one parameter and plot option value.
        
        Args:
            parameter: 'sigma', 'S', 'K', 'T', or 'r'
            value_range: (min, max) of parameter
            n_points: Number of points to calculate
            option_type: 'call' or 'put'
            dividend_yield: q
            method: 'black_scholes' or 'binomial'
        """
        original = getattr(self, parameter)
        values = np.linspace(value_range[0], value_range[1], n_points)
        option_values = []
        
        for v in values:
            setattr(self, parameter, v)
            if method == 'black_scholes':
                res = self.black_scholes(option_type, dividend_yield, verbose=False)
            else:
                res = self.binomial_tree(option_type, 50, dividend_yield, verbose=False)
            option_values.append(res['value'])
        
        setattr(self, parameter, original)  # Restore
        
        # Plot
        fig, ax = plt.subplots(figsize=(10, 5))
        ax.plot(values, option_values, color='#3B82F6', linewidth=2.5)
        ax.axhline(0, color='black', linewidth=0.5)
        ax.set_title(f'Real Option Sensitivity: {parameter} vs Option Value', 
                     fontweight='bold', fontsize=13)
        ax.set_xlabel(f'{parameter}', fontsize=12)
        ax.set_ylabel('Option Value (₹ Cr)', fontsize=12)
        ax.grid(alpha=0.3)
        plt.tight_layout()
        plt.show()
    
    def _print_result(self, result, verbose=True):
        if not verbose:
            return
        print("=" * 55)
        print(f"  {result.get('option_name', 'Real Option Valuation')}")
        print("=" * 55)
        for k, v in result.items():
            if isinstance(v, float):
                print(f"  {k:<25} {v:>20,.2f}")
            else:
                print(f"  {k:<25} {str(v):>20}")


# ========================================================
# EXAMPLE USAGE
# ========================================================
if __name__ == '__main__':
    # Option to Delay (SunPharma Patent)
    rov = RealOptionsValuer(S=2000, K=1500, T=5, sigma=0.45, r=0.07)
    bs_result = rov.black_scholes(option_type='call', dividend_yield=0.04)
    bt_result = rov.binomial_tree(option_type='call', n_steps=100, dividend_yield=0.04)
    
    # Sensitivity
    rov.sensitivity(parameter='sigma', value_range=(0.1, 0.8))
    rov.sensitivity(parameter='T', value_range=(0.5, 10))