퀀트 도전- 월간 퀀트스터디, LPPL) 파이썬 LPPL 코드로 과거 TSLA 주가 분석해보기(코드 수정됐으니 받아가신 분들 필독!!)

*****경량화 해서 트레이딩뷰에서 돌리려고 숫자 줄여봤던 부분이 남아있어서 full LPPL로 코드 수정해서 올립니다.

글 내용도 수정됨*****

2019-01-01 ~ 2021-01-01 까지 TSLA 수정 종가를 바탕으로 한 분석 결과입니다.

LPPL Confidence 23.08로 시장 과열에 가까운 경고단계

Critical time은 2021-01-30로 나왔습니다

이후 진짜로 그 쯤에서 버블이 꺼졌네요

자산군에 적용하는 것이 좋지만 만든김에 검증용으로 과거 테슬라 차트에 분석해봤습니다.

추가로 이미 퀀트팀에서 해주셨지만 SPY에 적용해볼까요?

오늘 날짜로 SPY는 LPPL confidence 6.38%로 주의필요 초기단계 쯤 되네요

앞으로 종종 모델 구현해서 공유해드리겠습니다.

아래는 만들어 본 코드입니다.(수정됨)

valley fellow 분들만 재미로 사용해보셨으면 합니다!

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.linalg import lstsq
from datetime import datetime, timedelta
import warnings
import os
import time
import random

warnings.filterwarnings('ignore')


class LPPL:
    def __init__(self, observations):
        """
        Initialize the LPPL model with historical price observations.
        
        Args:
            observations (pd.DataFrame): DataFrame with a DatetimeIndex and 'Adj Close' column.
        """
        self.observations = observations.sort_index()
        self.t = self.get_days_since_start()
        self.price = self.observations['Adj Close'].values
        self.log_price = np.log(self.price)
        # Scale time array to [0, 1] for optimization purposes
        self.ts = self.t / self.t[-1]

        # Parameter bounds based on typical LPPL references
        self.bounds = {
            'm': (0.1, 0.9),
            'omega': (2, 25),
            'phi': (0, 2 * np.pi)
            # tc bounds will be set dynamically during fitting
        }

        # Genetic algorithm settings
        self.population_size = 200
        self.generations = 100
        self.crossover_rate = 0.7
        self.mutation_rate = 0.1
        self.elitism_rate = 0.1

    def get_days_since_start(self):
        """
        Convert the DatetimeIndex into an array of days since the start date.
        
        Returns:
            np.ndarray: Array of days elapsed since the first date.
        """
        start_date = self.observations.index[0]
        return (self.observations.index - start_date).days.values

    def lppl_matrix(self, t, tc, m, omega, phi):
        """
        Construct the LPPL matrix for the given parameters.
        
        Args:
            t (np.ndarray): Scaled time array [0, 1].
            tc (float): Critical time.
            m (float): Bubble growth exponent (equivalent to beta in literature).
            omega (float): Angular frequency of oscillations.
            phi (float): Phase shift.
        
        Returns:
            np.ndarray: LPPL matrix of shape (n, 3).
        """
        t = np.array(t)
        dt = (tc - t).clip(0.00001)  # Avoid numerical instability near zero
        dt_pow_m = dt ** m
        oscillatory = dt_pow_m * np.cos(omega * np.log(dt) - phi)
        return np.vstack((np.ones_like(t), dt_pow_m, oscillatory)).T

    def _func_(self, params, t, log_price):
        """
        Objective function for non-linear parameter optimization.
        
        Args:
            params (np.ndarray): Parameter array [m, omega, phi, tc].
            t (np.ndarray): Scaled time array.
            log_price (np.ndarray): Log of the price series.
        
        Returns:
            float: Sum of squared errors.
        """
        m, omega, phi, tc = params
        try:
            lppl_m = self.lppl_matrix(t, tc, m, omega, phi)
            coefs = lstsq(lppl_m, log_price)[0]
            return np.sum((log_price - lppl_m @ coefs) ** 2)
        except (np.linalg.LinAlgError, ValueError):
            return 1e10  # Large penalty on error

    def simple_model_fit(self):
        """
        Phase 1: Fit a simple model (C=0, beta=1) to obtain an initial seed for tc.
        
        Returns:
            tuple: (tc, A, B) as initial guesses.
        """
        t_end = self.ts[-1]

        def simple_obj(tc_val, t, log_price):
            dt = tc_val - t
            X = np.vstack((np.ones_like(t), dt)).T
            try:
                coefs = lstsq(X, log_price)[0]
                return np.sum((log_price - X @ coefs) ** 2)
            except np.linalg.LinAlgError:
                return 1e10

        # Set tc candidates from t_end + 0.01 to t_end + 0.1 (per strict criteria)
        tc_candidates = np.linspace(t_end + 0.01, t_end + 0.1, 100)
        errors = [simple_obj(tc_val, self.ts, self.log_price) for tc_val in tc_candidates]
        best_tc = tc_candidates[np.argmin(errors)]

        # Solve for A, B with best_tc
        dt = best_tc - self.ts
        X = np.vstack((np.ones_like(self.ts), dt)).T
        try:
            A, B = lstsq(X, self.log_price)[0]
        except np.linalg.LinAlgError:
            A = np.mean(self.log_price)
            B = -0.01

        return best_tc, float(A), float(B)

    def initialize_population(self, tc_seed):
        """
        Phase 2: Initialize the population for the genetic algorithm.
        Modified: Set tc range to [t_end + 0.01, t_end + 0.1].
        
        Args:
            tc_seed (float): Initial tc seed from Phase 1.
        
        Returns:
            list: A list of tuples (params, fitness).
        """
        population = []
        t_end = self.ts[-1]
        tc_low, tc_high = t_end + 0.01, t_end + 0.1

        # First individual using tc_seed
        initial_individual = np.array([
            np.random.uniform(*self.bounds['m']),
            np.random.uniform(*self.bounds['omega']),
            np.random.uniform(*self.bounds['phi']),
            tc_seed
        ])
        fitness = self._func_(initial_individual, self.ts, self.log_price)
        population.append((initial_individual, fitness))

        # Remaining individuals
        for _ in range(self.population_size - 1):
            individual = np.array([
                np.random.uniform(*self.bounds['m']),
                np.random.uniform(*self.bounds['omega']),
                np.random.uniform(*self.bounds['phi']),
                np.random.uniform(tc_low, tc_high)
            ])
            fitness = self._func_(individual, self.ts, self.log_price)
            population.append((individual, fitness))

        return population

    def tournament_selection(self, population, tournament_size=3):
        """Tournament selection for the genetic algorithm."""
        tournament = random.sample(population, min(tournament_size, len(population)))
        return min(tournament, key=lambda x: x[1])

    def crossover(self, parent1, parent2):
        """Crossover operation."""
        alpha = np.random.random()
        return alpha * parent1 + (1 - alpha) * parent2

    def mutate(self, individual):
        """Mutation operation."""
        m_strength = (self.bounds['m'][1] - self.bounds['m'][0]) * 0.1
        omega_strength = (self.bounds['omega'][1] - self.bounds['omega'][0]) * 0.1
        phi_strength = (self.bounds['phi'][1] - self.bounds['phi'][0]) * 0.1
        tc_strength = 0.005  # Mutation strength for tc

        result = individual.copy()

        if np.random.random() < 0.3:
            result[0] = np.clip(result[0] + np.random.normal(0, m_strength), *self.bounds['m'])
        if np.random.random() < 0.3:
            result[1] = np.clip(result[1] + np.random.normal(0, omega_strength), *self.bounds['omega'])
        if np.random.random() < 0.3:
            result[2] = np.clip(result[2] + np.random.normal(0, phi_strength), *self.bounds['phi'])
        if np.random.random() < 0.3:
            t_end = self.ts[-1]
            tc_low, tc_high = t_end + 0.01, t_end + 0.1
            result[3] = np.clip(result[3] + np.random.normal(0, tc_strength), tc_low, tc_high)

        return result

    def genetic_algorithm_fit(self, tc_seed):
        """
        Phase 2-3: Optimize LPPL parameters using the genetic algorithm.
        
        Args:
            tc_seed (float): Initial tc seed from Phase 1.
        
        Returns:
            tuple: (best_params, best_fitness)
        """
        population = self.initialize_population(tc_seed)
        best_params, best_fitness = min(population, key=lambda x: x[1])

        for generation in range(self.generations):
            population.sort(key=lambda x: x[1])
            elite_count = max(1, int(self.population_size * self.elitism_rate))
            new_population = population[:elite_count].copy()

            while len(new_population) < self.population_size:
                if np.random.random() < self.crossover_rate:
                    parent1 = self.tournament_selection(population)[0]
                    parent2 = self.tournament_selection(population)[0]
                    child = self.crossover(parent1, parent2)
                    if np.random.random() < self.mutation_rate:
                        child = self.mutate(child)
                    fitness = self._func_(child, self.ts, self.log_price)
                    new_population.append((child, fitness))
                else:
                    new_population.append(self.tournament_selection(population))

            population = new_population[:self.population_size]
            current_best, current_fitness = min(population, key=lambda x: x[1])
            if current_fitness < best_fitness:
                best_params, best_fitness = current_best, current_fitness

            # Early stopping if no improvement after 10 generations
            if generation > 10 and abs(current_fitness - best_fitness) < 1e-6:
                break

        return best_params, best_fitness

    def multi_run_genetic_algorithm(self, tc_seed, runs=3):
        """
        Run the genetic algorithm multiple times to find the best result.
        
        Args:
            tc_seed (float): Initial tc seed from Phase 1.
            runs (int): Number of runs.
        
        Returns:
            tuple: (best_params, best_fitness)
        """
        best_overall_params = None
        best_overall_fitness = float('inf')
        original_mutation_rate = self.mutation_rate
        original_population_size = self.population_size

        for run in range(runs):
            print(f"Genetic Algorithm run {run+1}/{runs}...")
            self.mutation_rate = original_mutation_rate + (run * 0.05)
            self.population_size = original_population_size + (run * 50)
            params, fitness = self.genetic_algorithm_fit(tc_seed)
            if fitness < best_overall_fitness:
                best_overall_fitness = fitness
                best_overall_params = params

        self.mutation_rate = original_mutation_rate
        self.population_size = original_population_size
        return best_overall_params, best_overall_fitness

    def fine_tune_parameters(self, params):
        """
        Phase 3: Fine-tune parameters using local optimization.
        Modified: Set tc bounds to [t_last + 0.01, t_last + 0.1].
        
        Args:
            params (np.ndarray): Initial parameters [m, omega, phi, tc].
        
        Returns:
            np.ndarray: Refined parameters.
        """
        tuned_params = params.copy()
        param_bounds = [
            self.bounds['m'],
            self.bounds['omega'],
            self.bounds['phi'],
            (self.ts[-1] + 0.01, self.ts[-1] + 0.1)
        ]
        for idx in range(4):
            def objective(x):
                params_copy = tuned_params.copy()
                params_copy[idx] = x[0]
                return self._func_(params_copy, self.ts, self.log_price)

            result = minimize(objective, x0=[tuned_params[idx]],
                            ...