LPPL 코드 가져가셨던 분들 꼭 봐주세요!! (2025-03-03 추가 업데이트)

*****경량화 해서 트레이딩뷰에서 돌리려고 숫자 줄여봤던 부분이 남아있어서 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, ...

모든 시점에 모든 자산에서 0%가 나오는 듯... ㅠㅜ 어쩌면 야후 파이낸스 데이터 문제일지도 일단 도전해봤다는 것에 의의를 둡니다 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 warnings.filterwarnings('ignore') class LPPL: def __init__(self, observations): """ LPPL 모델 초기화. Parameters: observations (pd.DataFrame): 'Close' 가격 데이터와 DatetimeIndex를 포함하는 DataFrame. """ self.observations = observations self.t = self.get_days_since_start() # 시작일로부터의 일수 계산 self.price = self.observations['Close'].values self.log_price = np.log(self.price) self.ts = self.t / self.t[-1] # 최적화를 위해 t를 [0,1] 범위로 스케일 조정 self.bounds = [ (0.01, 0.99), # m의 범위 (6, 13), # omega의 범위 (0, 2 * np.pi), # phi의 범위 # tc의 범위는 fit()에서 추가 ] def get_days_since_start(self): """ DatetimeIndex를 시작일로부터의 일수 배열로 변환. Returns: np.array: 시작일로부터의 일수를 나타내는 정수 배열. """ start_date = self.observations.index[0] return (self.observations.index - start_date).days.values def lppl_matrix(self, t, tc, m, omega, phi): """ 주어진 파라미터에 대해 LPPL 행렬을 계산. 이 버전은 다음과 같은 3항식 형태를 사용합니다: log(p(t)) = A + B*(tc-t)^m + C*(tc-t)^m * cos(omega * log(tc-t) - phi) Args: t (np.ndarray): 스케일 조정된 시간 배열 (0~1). tc (float): 임계 시간. m (float): 버블 성장 지수. omega (float): 진동 주파수. phi (float): 위상 이동. Returns: np.ndarray: (n, 3) 모양의 LPPL 행렬. """ t = np.array(t) dt = (tc - t).clip(0) # 음수가 되지 않도록 함 dt_pow_m = dt**m oscillatory = dt_pow_m * np.cos(omega * np.log(dt + 1e-8) - phi) return np.vstack((np.ones_like(t), dt_pow_m, oscillatory)).T def _func_(self, x, *args): """ 최적화에서 최소화할 목적 함수 (비선형 파라미터에 대해). Args: x (np.ndarray): 비선형 파라미터 [m, omega, phi, tc] 배열. *args: 추가 인자 (t, log_price). Returns: float: 제곱 오차 합. """ m, omega, phi, tc = x t, log_price = args lppl_m = self.lppl_matrix(t, tc, m, omega, phi) try: coefs = lstsq(lppl_m, log_price)[0] except np.linalg.LinAlgError: return 1e10 # lstsq 실패 시 큰 값을 반환 return np.sum((log_price - lppl_m @ coefs) ** 2) def fit(self): """ L-BFGS-B 최적화 기법을 사용하여 가격 데이터에 LPPL 모델을 피팅. Returns: dict: 피팅된 파라미터와 기타 지표를 포함하는 딕셔너리. """ t_end = self.ts[-1] # tc의 경계: 마지막 데이터 이후로 설정 bounds = self.bounds + [(t_end + 0.01, t_end + 5)] init_guess = np.array([np.random.uniform(low, high) for low, high in bounds]) opt = minimize( self._func_, x0=init_guess, args=(self.ts, self.log_price), method='L-BFGS-B', bounds=bounds, ) m, omega, phi, tc = opt.x lppl_m = self.lppl_matrix(self.ts, tc, m, omega, phi) coefs = lstsq(lppl_m, self.log_price)[0] A, B, C = coefs[0], coefs[1], coefs[2] model = lppl_m @ coefs residuals = self.log_price - model mse = np.mean(residuals**2) damping_factor = abs(B) / abs(C) if C != 0 else np.inf num_oscillations = omega / (2 * np.pi) * np.log((tc - self.ts[0]) / (tc - self.ts[-1])) rel_error = np.std(residuals) / np.std(self.log_price) params = { 'A': A, 'B': B, 'C': C, 'm': m, 'omega': omega, 'phi': phi, 'tc': tc, 'tc_date': self.observations.index[0] + pd.Timedelta(days=int(tc * self.t[-1])), 'mse': mse, 'damping_factor': damping_factor, 'num_oscillations': num_oscillations, 'rel_error': rel_error, 'fit_success': opt.success, 'residuals': residuals } return params def check_bubble_conditions(self, params): """ 피팅된 파라미터가 버블 조건을 만족하는지 확인. Parameters: params (dict): 피팅된 LPPL 파라미터 딕셔너리. ...