← back to Watches

analytics/statistical_insights.py

544 lines

#!/usr/bin/env python3
"""
Statistical Insights and Time Series Analysis
Advanced statistical methods for watch price analysis
"""

import json
import numpy as np
import pandas as pd
from datetime import datetime
from scipy import stats
from collections import defaultdict
import warnings
warnings.filterwarnings('ignore')


class StatisticalAnalyzer:
    """Advanced statistical analysis for watch prices"""

    def __init__(self, data_path):
        self.data_path = data_path
        self.watches = None

    def load_data(self):
        """Load watch data"""
        with open(self.data_path, 'r') as f:
            data = json.load(f)
        self.watches = data['watches']
        return self.watches

    def calculate_moving_averages(self, window_years=[5, 10, 20]):
        """Calculate moving averages for price trends"""
        ma_results = []

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 5:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])

            # Convert to time series
            years = [p['year'] for p in prices_sorted]
            prices_arr = [p['price'] for p in prices_sorted]

            watch_ma = {
                'watch_id': watch['id'],
                'model': watch['model'],
                'series': watch['series'],
                'moving_averages': {}
            }

            # Calculate moving averages for each window
            for window in window_years:
                if len(prices_arr) >= window:
                    ma_values = []
                    for i in range(window - 1, len(prices_arr)):
                        window_prices = prices_arr[i - window + 1:i + 1]
                        ma_values.append({
                            'year': years[i],
                            'ma': round(np.mean(window_prices), 2)
                        })

                    watch_ma['moving_averages'][f'{window}_year'] = ma_values

            if watch_ma['moving_averages']:
                ma_results.append(watch_ma)

        return ma_results

    def analyze_best_worst_periods(self):
        """Identify best and worst performing decades for each watch"""
        period_analysis = []

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 3:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])

            # Calculate decade-by-decade performance
            decade_performance = defaultdict(list)

            for i in range(1, len(prices_sorted)):
                prev_price = prices_sorted[i - 1]['price']
                curr_price = prices_sorted[i]['price']
                prev_year = prices_sorted[i - 1]['year']
                curr_year = prices_sorted[i]['year']

                if prev_price > 0:
                    years_diff = curr_year - prev_year
                    if years_diff > 0:
                        annualized_return = ((curr_price / prev_price) ** (1 / years_diff) - 1) * 100
                        decade = (prev_year // 10) * 10

                        decade_performance[decade].append({
                            'period': f"{prev_year}-{curr_year}",
                            'return': annualized_return,
                            'start_price': prev_price,
                            'end_price': curr_price
                        })

            # Find best and worst decades
            decade_summary = []
            for decade, periods in decade_performance.items():
                avg_return = np.mean([p['return'] for p in periods])
                decade_summary.append({
                    'decade': f"{decade}s",
                    'avg_annual_return': round(avg_return, 2),
                    'periods': periods
                })

            if decade_summary:
                decade_summary.sort(key=lambda x: x['avg_annual_return'], reverse=True)

                best_decade = decade_summary[0]
                worst_decade = decade_summary[-1]

                period_analysis.append({
                    'watch_id': watch['id'],
                    'model': watch['model'],
                    'series': watch['series'],
                    'best_decade': best_decade,
                    'worst_decade': worst_decade,
                    'all_decades': decade_summary
                })

        return period_analysis

    def calculate_advanced_statistics(self):
        """Calculate comprehensive statistical metrics"""
        stats_results = []

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 3:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])
            price_values = [p['price'] for p in prices_sorted]

            # Basic statistics
            mean_price = np.mean(price_values)
            median_price = np.median(price_values)
            std_price = np.std(price_values)
            var_price = np.var(price_values)

            # Skewness and Kurtosis
            skewness = stats.skew(price_values)
            kurtosis = stats.kurtosis(price_values)

            # Coefficient of variation
            cv = (std_price / mean_price) * 100 if mean_price > 0 else 0

            # Quartiles and IQR
            q1, q3 = np.percentile(price_values, [25, 75])
            iqr = q3 - q1

            # Calculate returns distribution
            returns = []
            for i in range(1, len(prices_sorted)):
                prev_price = prices_sorted[i - 1]['price']
                curr_price = prices_sorted[i]['price']
                if prev_price > 0:
                    ret = ((curr_price - prev_price) / prev_price) * 100
                    returns.append(ret)

            if returns:
                mean_return = np.mean(returns)
                std_return = np.std(returns)
                sharpe_like_ratio = mean_return / std_return if std_return > 0 else 0
            else:
                mean_return = std_return = sharpe_like_ratio = 0

            stats_results.append({
                'watch_id': watch['id'],
                'model': watch['model'],
                'series': watch['series'],
                'price_statistics': {
                    'mean': round(mean_price, 2),
                    'median': round(median_price, 2),
                    'std_dev': round(std_price, 2),
                    'variance': round(var_price, 2),
                    'coefficient_of_variation': round(cv, 2),
                    'min': round(min(price_values), 2),
                    'max': round(max(price_values), 2),
                    'range': round(max(price_values) - min(price_values), 2),
                    'q1': round(q1, 2),
                    'q3': round(q3, 2),
                    'iqr': round(iqr, 2)
                },
                'distribution_metrics': {
                    'skewness': round(skewness, 4),
                    'kurtosis': round(kurtosis, 4),
                    'interpretation': {
                        'skewness': 'right-skewed' if skewness > 0.5 else ('left-skewed' if skewness < -0.5 else 'symmetric'),
                        'kurtosis': 'heavy-tailed' if kurtosis > 1 else ('light-tailed' if kurtosis < -1 else 'normal-tailed')
                    }
                },
                'return_statistics': {
                    'mean_return': round(mean_return, 2),
                    'std_return': round(std_return, 2),
                    'sharpe_ratio': round(sharpe_like_ratio, 4)
                }
            })

        return stats_results

    def calculate_percentile_rankings(self):
        """Calculate percentile rankings for key metrics"""
        # Compile all watches with key metrics
        watch_metrics = []

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 2:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])
            first_price = prices_sorted[0]['price']
            last_price = prices_sorted[-1]['price']
            years_span = prices_sorted[-1]['year'] - prices_sorted[0]['year']

            if years_span == 0 or first_price == 0:
                continue

            cagr = ((last_price / first_price) ** (1 / years_span) - 1) * 100
            total_appreciation = ((last_price - first_price) / first_price) * 100

            watch_metrics.append({
                'watch_id': watch['id'],
                'model': watch['model'],
                'series': watch['series'],
                'current_price': last_price,
                'cagr': cagr,
                'total_appreciation': total_appreciation
            })

        df = pd.DataFrame(watch_metrics)

        # Calculate percentiles
        rankings = []
        for _, watch in df.iterrows():
            price_percentile = stats.percentileofscore(df['current_price'], watch['current_price'])
            cagr_percentile = stats.percentileofscore(df['cagr'], watch['cagr'])
            appreciation_percentile = stats.percentileofscore(df['total_appreciation'], watch['total_appreciation'])

            rankings.append({
                'watch_id': watch['watch_id'],
                'model': watch['model'],
                'series': watch['series'],
                'percentile_rankings': {
                    'price': round(price_percentile, 1),
                    'cagr': round(cagr_percentile, 1),
                    'total_appreciation': round(appreciation_percentile, 1)
                },
                'interpretation': {
                    'price': self._interpret_percentile(price_percentile),
                    'cagr': self._interpret_percentile(cagr_percentile),
                    'appreciation': self._interpret_percentile(appreciation_percentile)
                }
            })

        return rankings

    def _interpret_percentile(self, percentile):
        """Interpret percentile ranking"""
        if percentile >= 90:
            return 'Top 10%'
        elif percentile >= 75:
            return 'Top 25%'
        elif percentile >= 50:
            return 'Above Average'
        elif percentile >= 25:
            return 'Below Average'
        else:
            return 'Bottom 25%'

    def hypothesis_testing(self):
        """Perform statistical hypothesis tests"""
        # Compile data by series
        series_data = defaultdict(list)

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 2:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])
            first_price = prices_sorted[0]['price']
            last_price = prices_sorted[-1]['price']
            years_span = prices_sorted[-1]['year'] - prices_sorted[0]['year']

            if years_span > 0 and first_price > 0:
                cagr = ((last_price / first_price) ** (1 / years_span) - 1) * 100
                series_data[watch['series']].append(cagr)

        # Test: Are Speedmaster returns different from Seamaster?
        test_results = []

        series_list = list(series_data.keys())
        for i in range(len(series_list)):
            for j in range(i + 1, len(series_list)):
                series1 = series_list[i]
                series2 = series_list[j]

                data1 = series_data[series1]
                data2 = series_data[series2]

                if len(data1) >= 2 and len(data2) >= 2:
                    # T-test
                    t_stat, p_value = stats.ttest_ind(data1, data2)

                    # Effect size (Cohen's d)
                    pooled_std = np.sqrt((np.var(data1) + np.var(data2)) / 2)
                    cohens_d = (np.mean(data1) - np.mean(data2)) / pooled_std if pooled_std > 0 else 0

                    test_results.append({
                        'comparison': f"{series1} vs {series2}",
                        'series1': {
                            'name': series1,
                            'mean_cagr': round(np.mean(data1), 2),
                            'sample_size': len(data1)
                        },
                        'series2': {
                            'name': series2,
                            'mean_cagr': round(np.mean(data2), 2),
                            'sample_size': len(data2)
                        },
                        't_statistic': round(t_stat, 4),
                        'p_value': round(p_value, 4),
                        'cohens_d': round(cohens_d, 4),
                        'significant': p_value < 0.05,
                        'interpretation': self._interpret_hypothesis_test(p_value, cohens_d)
                    })

        return test_results

    def _interpret_hypothesis_test(self, p_value, effect_size):
        """Interpret hypothesis test results"""
        if p_value < 0.05:
            if abs(effect_size) > 0.8:
                return "Statistically significant with large effect size"
            elif abs(effect_size) > 0.5:
                return "Statistically significant with medium effect size"
            else:
                return "Statistically significant with small effect size"
        else:
            return "No statistically significant difference"

    def calculate_risk_metrics(self):
        """Calculate risk-adjusted performance metrics"""
        risk_metrics = []

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 3:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])
            years = [p['year'] for p in prices_sorted]
            price_values = [p['price'] for p in prices_sorted]

            # Calculate returns
            returns = []
            for i in range(1, len(price_values)):
                if price_values[i - 1] > 0:
                    ret = (price_values[i] - price_values[i - 1]) / price_values[i - 1]
                    returns.append(ret)

            if not returns:
                continue

            mean_return = np.mean(returns)
            std_return = np.std(returns)

            # Downside deviation (only negative returns)
            negative_returns = [r for r in returns if r < 0]
            downside_deviation = np.std(negative_returns) if negative_returns else 0

            # Maximum drawdown
            cumulative_returns = np.cumprod([1 + r for r in returns])
            running_max = np.maximum.accumulate(cumulative_returns)
            drawdowns = (cumulative_returns - running_max) / running_max
            max_drawdown = np.min(drawdowns) if len(drawdowns) > 0 else 0

            # Sharpe ratio (assuming 0 risk-free rate)
            sharpe_ratio = mean_return / std_return if std_return > 0 else 0

            # Sortino ratio (using downside deviation)
            sortino_ratio = mean_return / downside_deviation if downside_deviation > 0 else 0

            # Calmar ratio (return / max drawdown)
            calmar_ratio = mean_return / abs(max_drawdown) if max_drawdown != 0 else 0

            risk_metrics.append({
                'watch_id': watch['id'],
                'model': watch['model'],
                'series': watch['series'],
                'risk_metrics': {
                    'mean_return': round(mean_return * 100, 2),
                    'volatility': round(std_return * 100, 2),
                    'downside_deviation': round(downside_deviation * 100, 2),
                    'max_drawdown': round(max_drawdown * 100, 2),
                    'sharpe_ratio': round(sharpe_ratio, 4),
                    'sortino_ratio': round(sortino_ratio, 4),
                    'calmar_ratio': round(calmar_ratio, 4)
                },
                'risk_adjusted_ranking': round(sharpe_ratio, 4)
            })

        # Sort by Sharpe ratio
        risk_metrics.sort(key=lambda x: x['risk_adjusted_ranking'], reverse=True)

        return risk_metrics

    def generate_summary_statistics(self):
        """Generate market-wide summary statistics"""
        all_cagrs = []
        all_prices = []
        all_appreciations = []

        for watch in self.watches:
            prices = watch.get('priceHistory', [])
            if len(prices) < 2:
                continue

            prices_sorted = sorted(prices, key=lambda x: x['year'])
            first_price = prices_sorted[0]['price']
            last_price = prices_sorted[-1]['price']
            years_span = prices_sorted[-1]['year'] - prices_sorted[0]['year']

            if years_span > 0 and first_price > 0:
                cagr = ((last_price / first_price) ** (1 / years_span) - 1) * 100
                appreciation = ((last_price - first_price) / first_price) * 100

                all_cagrs.append(cagr)
                all_prices.append(last_price)
                all_appreciations.append(appreciation)

        return {
            'cagr': {
                'mean': round(np.mean(all_cagrs), 2),
                'median': round(np.median(all_cagrs), 2),
                'std': round(np.std(all_cagrs), 2),
                'min': round(min(all_cagrs), 2),
                'max': round(max(all_cagrs), 2),
                'percentiles': {
                    '25th': round(np.percentile(all_cagrs, 25), 2),
                    '75th': round(np.percentile(all_cagrs, 75), 2),
                    '90th': round(np.percentile(all_cagrs, 90), 2)
                }
            },
            'current_prices': {
                'mean': round(np.mean(all_prices), 2),
                'median': round(np.median(all_prices), 2),
                'std': round(np.std(all_prices), 2),
                'min': round(min(all_prices), 2),
                'max': round(max(all_prices), 2)
            },
            'total_appreciation': {
                'mean': round(np.mean(all_appreciations), 2),
                'median': round(np.median(all_appreciations), 2),
                'std': round(np.std(all_appreciations), 2)
            }
        }


def main():
    """Run comprehensive statistical analysis"""
    print("=" * 80)
    print("OMEGA WATCH STATISTICAL INSIGHTS")
    print("=" * 80)

    analyzer = StatisticalAnalyzer('/root/Projects/watches/data/watches.json')

    # Load data
    print("\n[1/8] Loading watch data...")
    watches = analyzer.load_data()
    print(f"   Loaded {len(watches)} watches")

    # Moving averages
    print("\n[2/8] Calculating moving averages...")
    ma_results = analyzer.calculate_moving_averages()
    print(f"   Calculated MA for {len(ma_results)} watches")

    # Best/worst periods
    print("\n[3/8] Analyzing best and worst performing periods...")
    period_analysis = analyzer.analyze_best_worst_periods()
    print(f"   Analyzed {len(period_analysis)} watches")

    # Advanced statistics
    print("\n[4/8] Computing advanced statistics...")
    advanced_stats = analyzer.calculate_advanced_statistics()
    print(f"   Generated statistics for {len(advanced_stats)} watches")

    # Percentile rankings
    print("\n[5/8] Calculating percentile rankings...")
    rankings = analyzer.calculate_percentile_rankings()
    print(f"   Ranked {len(rankings)} watches")

    # Hypothesis testing
    print("\n[6/8] Performing hypothesis tests...")
    hypothesis_tests = analyzer.hypothesis_testing()
    print(f"   Performed {len(hypothesis_tests)} statistical tests")

    # Risk metrics
    print("\n[7/8] Calculating risk-adjusted metrics...")
    risk_metrics = analyzer.calculate_risk_metrics()
    print(f"   Calculated risk metrics for {len(risk_metrics)} watches")

    # Summary statistics
    print("\n[8/8] Generating market summary...")
    summary = analyzer.generate_summary_statistics()
    print(f"   Market median CAGR: {summary['cagr']['median']}%")
    print(f"   Market median price: ${summary['current_prices']['median']:,.0f}")

    # Compile results
    output = {
        'generated_at': datetime.now().isoformat(),
        'summary_statistics': summary,
        'moving_averages': ma_results,
        'period_analysis': period_analysis,
        'advanced_statistics': advanced_stats,
        'percentile_rankings': rankings,
        'hypothesis_tests': hypothesis_tests,
        'risk_metrics': risk_metrics
    }

    output_path = '/root/Projects/watches/analytics/statistical_insights.json'
    with open(output_path, 'w') as f:
        json.dump(output, f, indent=2)

    print(f"\n" + "=" * 80)
    print(f"Statistical analysis complete! Results saved to: {output_path}")
    print("=" * 80)

    return output


if __name__ == '__main__':
    results = main()