llpls.py

import numpy as np
import matplotlib.pyplot as plt
import yfinance as yf
import pandas as pd
from scipy.optimize import curve_fit
from fastdtw import fastdtw
from scipy.spatial.distance import euclidean
from scipy.signal import find_peaks, peak_widths
from scipy.interpolate import interp1d

def lppls(t, A, B, tc, m, C, w, phi):
    with np.errstate(invalid='ignore', divide='ignore'):
        return A + B * (tc - t) ** m * (1 + C * np.cos(w * np.log(tc - t) + phi))

def lppls_oscillator(t, A, B, tc, m, C, w, phi):
    return -B * (tc - t) ** (m - 1) * (1 + C * np.cos(w * np.log(tc - t) + phi))

# Download S&P 500 data from 2000 to 2022
sp500 = yf.download('^GSPC', start='2000-01-01', end='2022-12-31')
price = sp500['Adj Close']

# Rescale price data to different time scales
daily_price = price.resample('D').mean().dropna()
weekly_price = price.resample('W').mean().dropna()
monthly_price = price.resample('M').mean().dropna()

# Function to estimate initial parameters for LPPLS model
def estimate_initial_params(price_data):
    peaks, _ = find_peaks(price_data)
    widths = peak_widths(price_data, peaks, rel_height=0.5)
    mean_peak_width = np.mean(widths[0])
    frequency = 2 * np.pi / mean_peak_width
    phase = 0
    return [0.5, -0.3, 1.1, 0.3, 0.5, frequency, phase]

# Function to fit the LPPLS model and calculate the oscillator
def fit_lppls_and_calculate_oscillator(price_data):
    t = np.linspace(0, 1, len(price_data))
    initial_params = estimate_initial_params(price_data)
    with np.errstate(invalid='ignore', divide='ignore'):
        popt, _ = curve_fit(lppls, t, price_data, p0=initial_params, maxfev=10000)
    oscillator = lppls_oscillator(t, *popt)
    return oscillator

# Fit LPPLS model and calculate oscillator for different time scales
daily_oscillator = fit_lppls_and_calculate_oscillator(daily_price)
weekly_oscillator = fit_lppls_and_calculate_oscillator(weekly_price)
monthly_oscillator = fit_lppls_and_calculate_oscillator(monthly_price)

# Interpolate oscillators to have the same length
max_len = max(len(daily_oscillator), len(weekly_oscillator), len(monthly_oscillator))
t_daily = np.linspace(0, 1, len(daily_oscillator))
t_weekly = np.linspace(0, 1, len(weekly_oscillator))
t_monthly = np.linspace(0, 1, len(monthly_oscillator))
t_interp = np.linspace(0, 1, max_len)

interp_daily = interp1d(t_daily, daily_oscillator, kind='linear')(t_interp)
interp_weekly = interp1d(t_weekly, weekly_oscillator, kind='linear')(t_interp)
interp_monthly = interp1d(t_monthly, monthly_oscillator, kind='linear')(t_interp)

daily_flat = interp_daily.ravel().astype(float)
weekly_flat = interp_weekly.ravel().astype(float)
monthly_flat = interp_monthly.ravel().astype(float)

print(daily_flat.ndim)
print(weekly_flat.ndim)
print(monthly_flat.ndim)

# Calculate DTW distances between different time scales
distance_daily_weekly, _ = fastdtw(daily_flat, weekly_flat, dist=euclidean)
distance_daily_monthly, _ = fastdtw(daily_flat, monthly_flat, dist=euclidean)
distance_weekly_monthly, _ = fastdtw(weekly_flat, monthly_flat, dist=euclidean)

print(f"DTW distance between daily and weekly oscillators: {distance_daily_weekly}")
print(f"DTW distance between daily and monthly oscillators: {distance_daily_monthly}")
print(f"DTW distance between weekly and monthly oscillators: {distance_weekly_monthly}")

# Plot the results
plt.figure(figsize=(12, 6))

plt.subplot(311)
plt.plot(daily_price.index, daily_oscillator, label='Daily Oscillator')
plt.legend()
plt.title('LPPLS Oscillator - Daily Time Scale')

plt.subplot(312)
plt.plot(weekly_price.index, weekly_oscillator, label='Weekly Oscillator')
plt.legend()
plt.subplot(312)
plt.plot(weekly_price.index, weekly_oscillator, label='Weekly Oscillator')
plt.legend()
plt.title('LPPLS Oscillator - Weekly Time Scale')

plt.subplot(313)
plt.plot(monthly_price.index, monthly_oscillator, label='Monthly Oscillator')
plt.legend()
plt.title('LPPLS Oscillator - Monthly Time Scale')

plt.tight_layout()
plt.show()