1# Core imports
 2import numpy as np
 3import pandas as pd
 4import matplotlib.pyplot as plt
 5import seaborn as sns
 6import yfinance as yf
 7from datetime import datetime, timedelta
 8import warnings
 9warnings.filterwarnings('ignore')
10
11# Statistical tools
12from scipy import stats
13from statsmodels.tsa.stattools import adfuller, kpss, acf, pacf
14from statsmodels.stats.diagnostic import het_arch
15
16# Plotting style
17plt.style.use('seaborn-v0_8-whitegrid')
18sns.set_palette(sns.color_palette())
19
20## 1. The Two Competing Approaches <a name="approaches"></a>
21
22When transforming price data $P_t$ into returns, we have two options:
23
24### Simple (Percentage) Returns
25
26$$r_t = \frac{P_t - P_{t-1}}{P_{t-1}} = \frac{P_t}{P_{t-1}} - 1$$
27
28### Log Returns
29
30$$R_t = \log\left(\frac{P_t}{P_{t-1}}\right) = \log(P_t) - \log(P_{t-1})$$
31
32### Why This Matters
33
34Statistical models (including HMMs) require **stationary** observations. The wrong choice breaks fundamental assumptions.
35
36We'll prove that log returns satisfy **5 critical mathematical properties** that percentage returns do not.
37
38## 2. Data Collection <a name="data"></a>
39
40We'll use 4 assets with different characteristics:
41
42- **SPY**: S&P 500 ETF (market benchmark, ~\$500)
43- **AAPL**: Apple (~\$180)
44- **TSLA**: Tesla (~\$250, high volatility)
45- **KO**: Coca-Cola (~\$60, defensive)
46
47The different price levels ranging from \$60 to \\$500 will demonstrate scale invariance.
48
49
50```python
51# Download 2 years of data
52tickers = ['SPY', 'AAPL', 'TSLA', 'KO']
53end_date = datetime.now()
54start_date = end_date - timedelta(days=730)
55
56print(f"Downloading data from {start_date.date()} to {end_date.date()}\n")
57
58data = {}
59for ticker in tickers:
60    print(f"Downloading {ticker}...", end=' ')
61    df = yf.download(ticker, start=start_date, end=end_date, progress=False, auto_adjust=True)
62    data[ticker] = df
63    price_delta = df['Close'].max()-df['Close'].min()
64    print(f"({len(df)} days, price range: ${price_delta[ticker]:.2f})")
65
66print(f"\nDownloaded {len(tickers)} assets")

 1# Visualize raw prices to show scale differences
 2fig, axes = plt.subplots(2, 2, figsize=(15, 10))
 3axes = axes.ravel()
 4
 5for idx, ticker in enumerate(tickers):
 6    ax = axes[idx]
 7    prices = data[ticker]['Close'].squeeze()
 8    prices.plot(ax=ax, linewidth=2, color=f'C{idx}')
 9    
10    ax.set_title(f'{ticker} - Raw Price Data', fontsize=14, fontweight='bold')
11    ax.set_xlabel('Date')
12    ax.set_ylabel('Price ($)')
13    ax.grid(True, alpha=0.3)
14    
15    # Add price statistics
16    mean_price = prices.mean()
17    min_price = prices.min()
18    max_price = prices.max()
19    
20    ax.text(0.02, 0.98, 
21            f'Mean: ${mean_price:.2f}\nRange: ${min_price:.2f}-${max_price:.2f}',
22            transform=ax.transAxes, fontsize=10, verticalalignment='top',
23            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
24
25plt.tight_layout()
26plt.suptitle('Raw Price Data - Notice Vastly Different Scales', 
27             y=1.02, fontsize=16, fontweight='bold')
28plt.show()

 1# Demonstration: Correct vs Incorrect Display of Log Returns
 2ticker = 'SPY'
 3prices = data[ticker]['Close'].squeeze()
 4log_returns = np.log(prices / prices.shift(1)).dropna()
 5simple_returns = prices.pct_change().dropna()
 6
 7print("HOW TO DISPLAY LOG RETURNS")
 8print("=" * 80)
 9
10# Calculate statistics
11log_mean = log_returns.mean()
12log_std = log_returns.std()
13simple_mean = simple_returns.mean()
14simple_std = simple_returns.std()
15
16print("\n1. MEAN (Expected Daily Return)")
17print("-" * 80)
18print(f"   Log return mean (log space):        {log_mean:.6f}")
19print(f"   Simple return mean:                 {simple_mean:.6f}")
20print(f"\n   CORRECT:   (e^{log_mean:.6f} - 1) * 100 = {(np.exp(log_mean) - 1) * 100:.4f}%")
21print(f"   WRONG:     {log_mean:.6f} * 100 = {log_mean * 100:.4f}%")
22print(f"   Difference: {abs((np.exp(log_mean) - 1) * 100 - log_mean * 100):.4f}%")
23
24print("\n2. STANDARD DEVIATION (Volatility)")
25print("-" * 80)
26print(f"   Log return std (log space):         {log_std:.6f}")
27print(f"   Simple return std:                  {simple_std:.6f}")
28print(f"\n   Approximation: {log_std:.6f} * 100 ≈ {log_std * 100:.4f}%")
29print(f"   Simple std:    {simple_std:.6f} * 100 = {simple_std * 100:.4f}%")
30print(f"   Difference: {abs(log_std * 100 - simple_std * 100):.4f}% (very small!)")
31print(f"\n   Note: For daily returns (~1%), the approximation std*100 is acceptable")
32print(f"      But ALWAYS use (e^mean - 1)*100 for the mean!")
33
34print("\n3. EXAMPLE: Large Return")
35print("-" * 80)
36large_log_return = 0.10  # 10% in log space
37print(f"   Log return: {large_log_return:.4f}")
38print(f"   CORRECT:  (e^{large_log_return:.4f} - 1) * 100 = {(np.exp(large_log_return) - 1) * 100:.2f}%")
39print(f"   WRONG:    {large_log_return:.4f} * 100 = {large_log_return * 100:.2f}%")
40print(f"   Error: {abs((np.exp(large_log_return) - 1) * 100 - large_log_return * 100):.2f}%")
41
42print("\n" + "=" * 80)
43print("KEY TAKEAWAYS:")
44print("=" * 80)
45print("Mean: ALWAYS convert with (e^mean - 1) * 100")
46print("Std: Multiplying by 100 is acceptable for small daily returns")
47print("Individual returns: ALWAYS convert with (e^r - 1) * 100 for display")
48print("NEVER just multiply log returns by 100 without conversion!")
49print("=" * 80)

1log_return = 0.01  # This is in log space
2percentage = (np.exp(log_return) - 1) * 100  # ≈ 1.005%

1log_return = 0.01
2percentage = log_return * 100  # = 1.0% (INCORRECT!)

 1# Calculate both return types for SPY
 2ticker = 'SPY'
 3prices = data[ticker]['Close'].squeeze()
 4
 5# Simple returns: (P_t - P_{t-1}) / P_{t-1}
 6simple_returns = prices.pct_change().dropna()
 7
 8# Log returns: log(P_t / P_{t-1})
 9log_returns = np.log(prices / prices.shift(1)).dropna()
10
11print(f"{ticker} Returns Summary:")
12print("=" * 60)
13print(f"{'Metric':<30} {'Simple':<15} {'Log':<15}")
14print("=" * 60)
15# Convert log returns to percentage: (e^r - 1) * 100
16print(f"{'Mean (daily %)':<30} {simple_returns.mean()*100:>14.4f} {(np.exp(log_returns.mean()) - 1)*100:>14.4f}")
17print(f"{'Std Dev (daily %, approx.)':<30} {simple_returns.std()*100:>14.4f} {log_returns.std()*100:>14.4f}")
18print(f"{'Min (daily %)':<30} {simple_returns.min()*100:>14.2f} {(np.exp(log_returns.min()) - 1)*100:>14.2f}")
19print(f"{'Max (daily %)':<30} {simple_returns.max()*100:>14.2f} {(np.exp(log_returns.max()) - 1)*100:>14.2f}")
20print(f"{'Correlation':<30} {np.corrcoef(simple_returns, log_returns)[0,1]:>14.6f}")
21print("=" * 60)

 1# Visual comparison
 2fig, axes = plt.subplots(2, 3, figsize=(18, 10))
 3
 4# 1. Time series comparison
 5ax = axes[0, 0]
 6simple_returns.plot(ax=ax, label='Simple Returns', alpha=0.7, linewidth=1)
 7log_returns.plot(ax=ax, label='Log Returns', alpha=0.7, linewidth=1)
 8ax.set_title('Returns Over Time', fontweight='bold')
 9ax.set_xlabel('Date')
10ax.set_ylabel('Return')
11ax.legend()
12ax.grid(True, alpha=0.3)
13
14# 2. Scatter plot showing relationship
15# Convert log returns to percentage for comparison
16log_returns_pct = (np.exp(log_returns) - 1) * 100
17ax = axes[0, 1]
18ax.scatter(simple_returns * 100, log_returns_pct, alpha=0.3, s=10)
19ax.plot([-5, 5], [-5, 5], 'r--', label='y=x line', linewidth=2)
20ax.set_xlabel('Simple Returns (%)')
21ax.set_ylabel('Log Returns (converted to %)')
22ax.set_title('Log vs Simple Returns', fontweight='bold')
23ax.legend()
24ax.grid(True, alpha=0.3)
25ax.text(0.05, 0.95, f'Correlation: {np.corrcoef(simple_returns, log_returns)[0,1]:.6f}',
26        transform=ax.transAxes, verticalalignment='top',
27        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
28
29# 3. Distribution comparison
30ax = axes[0, 2]
31ax.hist(simple_returns * 100, bins=50, alpha=0.5, label='Simple', density=True, color='blue')
32ax.hist(log_returns_pct, bins=50, alpha=0.5, label='Log', density=True, color='orange')
33ax.set_xlabel('Return (%)')
34ax.set_ylabel('Density')
35ax.set_title('Distribution Comparison', fontweight='bold')
36ax.legend()
37ax.grid(True, alpha=0.3)
38
39# 4. Q-Q plot: Simple returns
40ax = axes[1, 0]
41stats.probplot(simple_returns, dist="norm", plot=ax)
42ax.set_title('Q-Q Plot: Simple Returns', fontweight='bold')
43ax.grid(True, alpha=0.3)
44
45# 5. Q-Q plot: Log returns
46ax = axes[1, 1]
47stats.probplot(log_returns, dist="norm", plot=ax)
48ax.set_title('Q-Q Plot: Log Returns', fontweight='bold')
49ax.grid(True, alpha=0.3)
50
51# 6. Cumulative returns (testing additivity)
52ax = axes[1, 2]
53cum_simple = (1 + simple_returns).cumprod() - 1
54cum_log = np.exp(log_returns.cumsum()) - 1
55cum_simple.plot(ax=ax, label='From Simple Returns', linewidth=2)
56cum_log.plot(ax=ax, label='From Log Returns', linewidth=2, linestyle='--')
57ax.set_title('Cumulative Returns', fontweight='bold')
58ax.set_xlabel('Date')
59ax.set_ylabel('Cumulative Return')
60ax.legend()
61ax.grid(True, alpha=0.3)
62
63plt.tight_layout()
64plt.show()

 1# Comprehensive stationarity testing
 2from statsmodels.tsa.stattools import adfuller, kpss
 3from statsmodels.tools.sm_exceptions import InterpolationWarning
 4
 5def test_stationarity(series, name):    
 6    """Run multiple stationarity tests with graceful handling of edge cases"""    
 7    print(f"\n{name}:")    
 8    print("-" * 70)        
 9    # ADF Test    
10    adf_result = adfuller(series.dropna(), autolag='AIC')    
11    print(f"ADF Test:")    
12    print(f"  Statistic: {adf_result[0]:.4f}")    
13    print(f"  p-value: {adf_result[1]:.4f}")    
14    print(f"  Result: {'STATIONARY' if adf_result[1] < 0.05 else 'NON-STATIONARY'}")        
15    # KPSS Test with edge case handling    
16    print(f"\nKPSS Test:")        
17    # Use conservative nlags to avoid edge cases while remaining statistically valid    
18    n_obs = len(series.dropna())    
19    nlags = min(int(12 * (n_obs / 100) ** 0.25), n_obs // 4)        
20    # Catch interpolation warnings (they actually reinforce our conclusions!)    
21    with warnings.catch_warnings(record=True) as w:        
22        warnings.simplefilter("always", InterpolationWarning)        
23        kpss_result = kpss(series.dropna(), regression='c', nlags=nlags)                
24        interp_warnings = [warn for warn in w if issubclass(warn.category, InterpolationWarning)]        
25    
26    print(f"  Statistic: {kpss_result[0]:.4f}")    
27    print(f"  p-value: {kpss_result[1]:.4f}")    
28    print(f"  Critical value (5%): {kpss_result[3]['5%']:.4f}")        
29    # Explain if test statistic was extreme    
30    if interp_warnings:        
31        if kpss_result[0] > kpss_result[3]['10%']:            
32            print(f"  [WARNING] Extremely high statistic (very strong non-stationarity)")        
33        else:            
34            print(f"  [INFO] Extremely low statistic (very strong stationarity)")       
35    
36    print(f"  Result: {'STATIONARY' if kpss_result[1] > 0.05 else 'NON-STATIONARY'}")        
37    
38    return {        
39        'adf_stat': adf_result[0],        
40        'adf_pvalue': adf_result[1],        
41        'kpss_stat': kpss_result[0],        
42        'kpss_pvalue': kpss_result[1],        
43        'is_stationary': (adf_result[1] < 0.05) and (kpss_result[1] > 0.05)    
44    }
45            
46# Test both prices and log returns for SPY
47print("=" * 70)
48print("STATIONARITY TESTS (SPY)")
49print("=" * 70)
50prices_result = test_stationarity(prices, "Raw Prices")
51log_result = test_stationarity(log_returns, "Log Returns")
52print("=" * 70)

 1# Comprehensive stationarity testingfrom statsmodels.tsa.stattools 
 2from statsmodels.tools.sm_exceptions import InterpolationWarning
 3
 4def test_stationarity(series, name):    
 5    """Run multiple stationarity tests with graceful handling of edge cases"""    
 6    print(f"\n{name}:")    
 7    print("-" * 70)        
 8    # ADF Test    
 9    adf_result = adfuller(series.dropna(), autolag='AIC')    
10    print(f"ADF Test:")    
11    print(f"  Statistic: {adf_result[0]:.4f}")    
12    print(f"  p-value: {adf_result[1]:.4f}")    
13    print(f"  Result: {'STATIONARY' if adf_result[1] < 0.05 else 'NON-STATIONARY'}")        
14    # KPSS Test with edge case handling    
15    print(f"\nKPSS Test:")        
16    # Use conservative nlags to avoid edge cases while remaining statistically valid    
17    n_obs = len(series.dropna())    
18    nlags = min(int(12 * (n_obs / 100) ** 0.25), n_obs // 4)        
19    # Catch interpolation warnings (they actually reinforce our conclusions!)    
20    with warnings.catch_warnings(record=True) as w:        
21        warnings.simplefilter("always", InterpolationWarning)        
22        kpss_result = kpss(series.dropna(), regression='c', nlags=nlags)                
23        interp_warnings = [warn for warn in w if issubclass(warn.category, InterpolationWarning)]        
24        print(f"  Statistic: {kpss_result[0]:.4f}")    
25        print(f"  p-value: {kpss_result[1]:.4f}")    
26        print(f"  Critical value (5%): {kpss_result[3]['5%']:.4f}")        
27        # Explain if test statistic was extreme    
28        if interp_warnings:        
29            if kpss_result[0] > kpss_result[3]['10%']:            
30                print(f"  [WARNING] Extremely high statistic (very strong non-stationarity)")        
31            else:            
32                print(f"  [INFO] Extremely low statistic (very strong stationarity)")       
33                print(f"  Result: {'STATIONARY' if kpss_result[1] > 0.05 else 'NON-STATIONARY'}")        
34                return {        
35                    'adf_stat': adf_result[0],        
36                    'adf_pvalue': adf_result[1],        
37                    'kpss_stat': kpss_result[0],        
38                    'kpss_pvalue': kpss_result[1],        
39                    'is_stationary': (adf_result[1] < 0.05) and (kpss_result[1] > 0.05)    
40            }
41            
42# Test all three for SPY
43print("=" * 70)
44print("STATIONARITY TESTS (SPY)")
45print("=" * 70)
46prices_result = test_stationarity(prices, "Raw Prices")
47log_result = test_stationarity(log_returns, "Log Returns")
48print("\n" + "=" * 70)
49print("CONCLUSION:")
50print("=" * 70)
51if log_result['is_stationary']:    
52    print("Log returns are STATIONARY (both tests pass)")    
53    print("  -> Safe to use in statistical models")
54else:    
55    print("Log returns show some non-stationarity")    
56    if not prices_result['is_stationary']:    
57        print("Prices are NON-STATIONARY (as expected)")    
58        print("  -> Cannot use directly in statistical models")

 1# Rolling statistics to visualize stationarity
 2window = 60  # 60-day rolling window
 3
 4fig, axes = plt.subplots(2, 2, figsize=(16, 10))
 5
 6# Top row: Prices
 7# 1. Price level (shows trend)
 8ax = axes[0, 0]
 9prices.plot(ax=ax, linewidth=1.5, label='Price', color='blue')
10prices.rolling(window).mean().plot(ax=ax, linewidth=2, label=f'{window}-day MA', color='red')
11ax.set_title('Prices Show Trend (Non-Stationary)', fontweight='bold')
12ax.set_ylabel('Price ($)')
13ax.legend()
14ax.grid(True, alpha=0.3)
15
16# 2. Price rolling std (shows changing variance)
17ax = axes[0, 1]
18prices.rolling(window).std().plot(ax=ax, linewidth=2, color='purple')
19ax.set_title('Price Volatility Changes Over Time', fontweight='bold')
20ax.set_ylabel('Rolling Std Dev ($)')
21ax.grid(True, alpha=0.3)
22
23# Bottom row: Log Returns
24# 3. Log returns (no trend)
25ax = axes[1, 0]
26log_returns.plot(ax=ax, linewidth=0.8, alpha=0.7, label='Log Return', color='green')
27log_returns.rolling(window).mean().plot(ax=ax, linewidth=2, label=f'{window}-day MA', color='red')
28ax.axhline(y=0, color='black', linestyle='--', linewidth=1)
29ax.set_title('Log Returns Mean-Reverting (Stationary)', fontweight='bold')
30ax.set_ylabel('Log Return')
31ax.legend()
32ax.grid(True, alpha=0.3)
33
34# 4. Log return rolling std (relatively stable)
35ax = axes[1, 1]
36log_returns.rolling(window).std().plot(ax=ax, linewidth=2, color='orange')
37ax.set_title('Log Return Volatility More Stable', fontweight='bold')
38ax.set_ylabel('Rolling Std Dev')
39ax.grid(True, alpha=0.3)
40
41plt.tight_layout()
42plt.show()

 1# Rolling statistics to visualize stationarity
 2window = 60  # 60-day rolling window
 3
 4fig, axes = plt.subplots(2, 2, figsize=(16, 10))
 5
 6# Top row: Prices
 7# 1. Price level (shows trend)
 8ax = axes[0, 0]
 9prices.plot(ax=ax, linewidth=1.5, label='Price', color='blue')
10prices.rolling(window).mean().plot(ax=ax, linewidth=2, label=f'{window}-day MA', color='red')
11ax.set_title('Prices Show Trend (Non-Stationary)', fontweight='bold')
12ax.set_ylabel('Price ($)')
13ax.legend()
14ax.grid(True, alpha=0.3)
15
16# 2. Price rolling std (shows changing variance)
17ax = axes[0, 1]
18prices.rolling(window).std().plot(ax=ax, linewidth=2, color='purple')
19ax.set_title('Price Volatility Changes Over Time', fontweight='bold')
20ax.set_ylabel('Rolling Std Dev ($)')
21ax.grid(True, alpha=0.3)
22
23# Bottom row: Log Returns
24# 3. Log returns (no trend)
25ax = axes[1, 0]
26log_returns.plot(ax=ax, linewidth=0.8, alpha=0.7, label='Log Return', color='green')
27log_returns.rolling(window).mean().plot(ax=ax, linewidth=2, label=f'{window}-day MA', color='red')
28ax.axhline(y=0, color='black', linestyle='--', linewidth=1)
29ax.set_title('Log Returns Mean-Reverting (Stationary)', fontweight='bold')
30ax.set_ylabel('Log Return')
31ax.legend()
32ax.grid(True, alpha=0.3)
33
34# 4. Log return rolling std (relatively stable)
35ax = axes[1, 1]
36log_returns.rolling(window).std().plot(ax=ax, linewidth=2, color='orange')
37ax.set_title('Log Return Volatility More Stable', fontweight='bold')
38ax.set_ylabel('Rolling Std Dev')
39ax.grid(True, alpha=0.3)
40
41plt.tight_layout()
42plt.show()
43
44print("\nPrices: Mean and variance change over time -> NON-STATIONARY")
45print("Log Returns: Mean ~0, variance relatively stable -> STATIONARY")

 1# Demonstrate time additivity
 2print("TIME ADDITIVITY DEMONSTRATION")
 3print("=" * 70)
 4
 5# Pick a random 5-day period
 6start_idx = 100
 7period = 5
 8period_dates = prices.index[start_idx:start_idx+period+1]
 9period_prices = prices.iloc[start_idx:start_idx+period+1]
10
11print(f"\n5-Day Period: {period_dates[0].date()} to {period_dates[-1].date()}")
12print("-" * 70)
13
14# Calculate daily returns
15daily_simple = period_prices.pct_change().dropna()
16daily_log = np.log(period_prices / period_prices.shift(1)).dropna()
17
18print("\nDaily Returns:")
19for i, (date, simple_ret, log_ret) in enumerate(zip(daily_simple.index, daily_simple, daily_log)):
20    # Convert log return to percentage for display
21    log_ret_pct = (np.exp(log_ret) - 1) * 100
22    print(f"  Day {i+1} ({date.date()}): Simple={simple_ret*100:>6.2f}%, Log={log_ret_pct:>6.2f}%")
23
24# Total return (direct calculation)
25total_return_direct = (period_prices.iloc[-1] / period_prices.iloc[0]) - 1
26
27# Total return from simple returns (compounding)
28total_simple_compound = (1 + daily_simple).prod() - 1
29
30# Total return from log returns (summing)
31total_log_sum = daily_log.sum()
32total_log_to_simple = np.exp(total_log_sum) - 1
33
34print("\n" + "=" * 70)
35print("MULTI-PERIOD RETURN CALCULATION:")
36print("=" * 70)
37print(f"\nDirect calculation:      {total_return_direct*100:>6.2f}%")
38print(f"Simple returns (compound): {total_simple_compound*100:>6.2f}%  <- Requires multiplication")
39print(f"Log returns (sum):         {total_log_to_simple*100:>6.2f}%  <- Just add them up!")
40
41print("\nLog returns: R_total = R_1 + R_2 + R_3 + R_4 + R_5")
42print("Simple returns: r_total = (1+r_1)(1+r_2)(1+r_3)(1+r_4)(1+r_5) - 1")
43print("\nNote: Log returns are linearly additive across time!")

 1# Demonstrate scale invariance
 2print("SCALE INVARIANCE DEMONSTRATION")
 3print("=" * 70)
 4
 5# Compare all 4 assets
 6print("\nAsset Price Levels:")
 7for ticker in tickers:
 8    price = data[ticker]['Close'].squeeze().iloc[-1]
 9    print(f"  {ticker}: ${price:.2f}")
10
11print("\n" + "-" * 70)
12print("Scenario: Each asset moves $5 in one day")
13print("-" * 70)
14
15dollar_move = 5.0
16
17for ticker in tickers:
18    current_price = data[ticker]['Close'].squeeze().iloc[-1]
19    new_price = current_price + dollar_move
20    
21    # Simple return
22    simple_ret = dollar_move / current_price
23    
24    # Log return
25    log_ret = np.log(new_price / current_price)
26    
27    # Convert log return to percentage for display
28    log_ret_pct = (np.exp(log_ret) - 1) * 100
29    
30    print(f"\n{ticker} (${current_price:.2f} → ${new_price:.2f}):")
31    print(f"  Simple return: {simple_ret*100:>6.2f}%")
32    print(f"  Log return:    {log_ret_pct:>6.2f}%")
33
34print("\n" + "=" * 70)
35print("OBSERVATION:")
36print("=" * 70)
37print("Both simple and log returns correctly show:")
38print("  - KO ($60): $5 move = larger percentage (8.3%)")
39print("  - SPY ($500): $5 move = smaller percentage (1.0%)")
40print("\nBoth adjust for scale differences")
41print("But log returns have better mathematical properties (additivity, etc.)")

 1# Demonstrate symmetric treatment
 2print("SYMMETRIC TREATMENT OF GAINS AND LOSSES")
 3print("=" * 70)
 4
 5# Example: $100 stock
 6price_0 = 100
 7moves = [1.10, 1.20, 1.50, 2.00]  # 10%, 20%, 50%, 100% gains
 8
 9print("\nScenario: Up then down by same factor\n")
10print(f"{'Factor':<10} {'Up (Simple)':<15} {'Down (Simple)':<15} {'Up (Log)':<15} {'Down (Log)':<15}")
11print("-" * 70)
12
13for factor in moves:
14    # Up move
15    price_up = price_0 * factor
16    simple_up = (price_up / price_0) - 1
17    log_up = np.log(price_up / price_0)
18    
19    # Down move (reverse)
20    price_down = price_up / factor  # Should return to price_0
21    simple_down = (price_down / price_up) - 1
22    log_down = np.log(price_down / price_up)
23    
24    # Convert log returns to percentage for display
25    log_up_pct = (np.exp(log_up) - 1) * 100
26    log_down_pct = (np.exp(log_down) - 1) * 100
27    
28    print(f"{factor:<10.2f} {simple_up*100:>13.1f}% {simple_down*100:>13.1f}% {log_up_pct:>13.1f}% {log_down_pct:>13.1f}%")
29
30print("\n" + "=" * 70)
31print("OBSERVATION:")
32print("=" * 70)
33print("Simple returns: Gains and losses have different magnitudes")
34print("  Example: +100% gain requires -50% loss to break even (asymmetric)")
35print("\nLog returns: Gains and losses have equal magnitudes (symmetric)")
36print("  Example: +100% log gain requires -100% log loss to break even")
37print("\nNote: This symmetry is crucial for statistical modeling!")

 1# Calculate log returns for all assets
 2all_log_returns = {}
 3for ticker in tickers:
 4    prices = data[ticker]['Close'].squeeze()
 5    all_log_returns[ticker] = np.log(prices / prices.shift(1)).dropna()
 6
 7# Create comparison table
 8print("CROSS-ASSET LOG RETURN STATISTICS")
 9print("=" * 80)
10print(f"{'Asset':<8} {'Mean (daily %)':<15} {'Std (daily %)':<15} {'Ann. Return %':<15} {'Ann. Vol %':<15}")
11print("=" * 80)
12
13for ticker in tickers:
14    lr = all_log_returns[ticker]
15    # Convert mean log return to percentage properly: (e^mean - 1) * 100
16    mean_daily = (np.exp(lr.mean()) - 1) * 100
17    # Std can be multiplied by 100 (volatility approximation is valid for small returns)
18    std_daily = lr.std() * 100
19    # Annualized return: convert from log space
20    ann_return = (np.exp(lr.mean() * 252) - 1) * 100
21    ann_vol = lr.std() * np.sqrt(252) * 100
22    
23    print(f"{ticker:<8} {mean_daily:>14.3f} {std_daily:>14.2f} {ann_return:>14.1f} {ann_vol:>14.1f}")
24
25print("=" * 80)
26print("\nNote: All assets expressed in same units (log returns)")
27print("   Can now meaningfully compare volatilities and distributions")

 1# Visual comparison of distributions
 2fig, axes = plt.subplots(2, 2, figsize=(15, 10))
 3
 4# 1. Distribution overlay
 5ax = axes[0, 0]
 6for ticker in tickers:
 7    all_log_returns[ticker].hist(bins=50, alpha=0.5, label=ticker, ax=ax, density=True)
 8ax.set_xlabel('Log Return')
 9ax.set_ylabel('Density')
10ax.set_title('Log Return Distributions (Comparable!)', fontweight='bold')
11ax.legend()
12ax.grid(True, alpha=0.3)
13
14# 2. Volatility comparison
15ax = axes[0, 1]
16vols = [all_log_returns[t].std() * np.sqrt(252) * 100 for t in tickers]
17colors = ['C0', 'C1', 'C2', 'C3']
18ax.bar(tickers, vols, color=colors, alpha=0.7)
19ax.set_ylabel('Annualized Volatility (%)')
20ax.set_title('Volatility Comparison (Directly Comparable)', fontweight='bold')
21ax.grid(True, alpha=0.3, axis='y')
22
23# 3. Correlation matrix
24ax = axes[1, 0]
25returns_df = pd.DataFrame(all_log_returns)
26corr = returns_df.corr()
27im = ax.imshow(corr, cmap='RdYlGn', aspect='auto', vmin=-1, vmax=1)
28ax.set_xticks(range(len(tickers)))
29ax.set_yticks(range(len(tickers)))
30ax.set_xticklabels(tickers)
31ax.set_yticklabels(tickers)
32for i in range(len(tickers)):
33    for j in range(len(tickers)):
34        text = ax.text(j, i, f'{corr.iloc[i, j]:.2f}',
35                       ha="center", va="center", color="black", fontweight='bold')
36plt.colorbar(im, ax=ax)
37ax.set_title('Log Return Correlation Matrix', fontweight='bold')
38
39# 4. Rolling correlation (SPY vs others)
40ax = axes[1, 1]
41window = 60
42for ticker in ['AAPL', 'TSLA', 'KO']:
43    rolling_corr = returns_df['SPY'].rolling(window).corr(returns_df[ticker])
44    rolling_corr.plot(ax=ax, label=f'SPY-{ticker}', linewidth=2)
45ax.set_ylabel('Rolling Correlation')
46ax.set_title(f'{window}-Day Rolling Correlation with SPY', fontweight='bold')
47ax.legend()
48ax.grid(True, alpha=0.3)
49ax.axhline(y=0, color='black', linestyle='--', linewidth=1)
50
51plt.tight_layout()
52plt.show()
53
54print("\nLog returns enable meaningful cross-asset comparison")
55print("Can identify common patterns and regime correlations")
56print("Foundation for multi-asset regime detection")

 1# Comprehensive distribution tests
 2from scipy.stats import shapiro, jarque_bera, normaltest
 3
 4print("DISTRIBUTION ANALYSIS")
 5print("=" * 80)
 6
 7for ticker in tickers:
 8    lr = all_log_returns[ticker]
 9    
10    print(f"\n{ticker}:")
11    print("-" * 80)
12    
13    # Moments
14    print(f"  Mean:     {lr.mean():.6f}")
15    print(f"  Std Dev:  {lr.std():.6f}")
16    print(f"  Skewness: {lr.skew():.4f} (0 = symmetric)")
17    print(f"  Kurtosis: {lr.kurtosis():.4f} (0 = normal, >0 = fat tails)")
18    
19    # Normality tests
20    jb_statistic, jb_pval = jarque_bera(lr)
21    print(f"\n  Jarque-Bera test p-value: {jb_pval:.4f} and statistic: {jb_statistic:.4f}")
22    print(f"    {'Reject normality (p < 0.05)' if jb_pval < 0.05 else 'Cannot reject normality'}")
23
24    shapiro_statistic, shapiro_pvalue = shapiro(lr)
25    print(f"\n  Shapiro-Wilk test p-value: {shapiro_pvalue:.4f} and statistic: {shapiro_statistic:.4f}")
26    print(f"    {'Reject normality (p < 0.05)' if shapiro_pvalue < 0.05 else 'Cannot reject normality'}")