!pip -q install yfinance scipy matplotlib pandas numpy

import yfinance as yf
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar
from datetime import datetime

GENESIS = pd.Timestamp("2009-01-03")
START_DATE = "2014-09-17"
END_DATE = "2032-01-01"
FIT_YEARS = 5

LAMBDA_1 = 2.081
LAMBDA_2 = 1.305

OMEGA_1 = 2 * np.pi / np.log(LAMBDA_1)
OMEGA_2 = 2 * np.pi / np.log(LAMBDA_2)

DACOINMINSTER_COEFFICIENT = 4.42e-17
DACOINMINSTER_EXPONENT = 5.6

PERRENOD_COEFFICIENT = 0.00688
PERRENOD_EXPONENT = 5.848

btc = yf.download(
    "BTC-USD",
    start=START_DATE,
    auto_adjust=True,
    progress=False
)

if isinstance(btc.columns, pd.MultiIndex):
    btc.columns = btc.columns.get_level_values(0)

btc = btc[["Close"]].dropna()
btc.index = pd.to_datetime(btc.index).tz_localize(None)
btc = btc.rename(columns={"Close": "BTC"})

forecast_dates = pd.date_range(
    start=btc.index.min(),
    end=END_DATE,
    freq="D"
)

data = pd.DataFrame(index=forecast_dates)
data = data.join(btc, how="left")

data["Days"] = (data.index - GENESIS).days
data["Age"] = data["Days"] / 365.2425

data = data[data["Days"] > 0].copy()

data["Dacoinminster"] = (
    DACOINMINSTER_COEFFICIENT
    * data["Days"] ** DACOINMINSTER_EXPONENT
)

data["Perrenod"] = (
    PERRENOD_COEFFICIENT
    * data["Age"] ** PERRENOD_EXPONENT
)

fit_end = btc.index.max()
fit_start = fit_end - pd.DateOffset(years=FIT_YEARS)

fit_data = data.loc[
    (data.index >= fit_start)
    & (data.index <= fit_end)
].dropna(subset=["BTC"]).copy()

t = fit_data["Age"].values
actual = fit_data["BTC"].values
baseline = fit_data["Perrenod"].values

y = np.log10(actual / baseline)

X = np.column_stack([
    np.ones(len(t)),
    np.cos(OMEGA_1 * np.log(t)),
    np.sin(OMEGA_1 * np.log(t)),
    np.cos(OMEGA_2 * np.log(t)),
    np.sin(OMEGA_2 * np.log(t))
])

coefficients, residuals, rank, singular_values = np.linalg.lstsq(
    X,
    y,
    rcond=None
)

c0, c1, c2, c3, c4 = coefficients

age = data["Age"].values
log_age = np.log(age)

oscillation_1 = (
    c1 * np.cos(OMEGA_1 * log_age)
    + c2 * np.sin(OMEGA_1 * log_age)
)

oscillation_2 = (
    c3 * np.cos(OMEGA_2 * log_age)
    + c4 * np.sin(OMEGA_2 * log_age)
)

data["Oscillation_1"] = (
    data["Perrenod"]
    * 10 ** (c0 + oscillation_1)
)

data["Oscillation_2"] = (
    data["Perrenod"]
    * 10 ** (c0 + oscillation_2)
)

data["Combined"] = (
    data["Perrenod"]
    * 10 ** (c0 + oscillation_1 + oscillation_2)
)

data["Average"] = (
    data["Oscillation_1"]
    + data["Oscillation_2"]
    + data["Combined"]
) / 3

forecast_start = pd.Timestamp("2029-01-01")
forecast_end = pd.Timestamp("2030-12-31")

forecast_window = data.loc[
    (data.index >= forecast_start)
    & (data.index <= forecast_end)
]

peak_date = forecast_window["Combined"].idxmax()
peak_price = forecast_window.loc[peak_date, "Combined"]

plt.figure(figsize=(18, 10))

plt.plot(
    data.index,
    data["BTC"],
    color="orange",
    linewidth=1.5,
    label="Bitcoin Actual Price"
)

plt.plot(
    data.index,
    data["Dacoinminster"],
    color="blue",
    linewidth=1.8,
    label="Dacoinminster Power Law"
)

plt.plot(
    data.index,
    data["Perrenod"],
    color="purple",
    linewidth=1.8,
    label="Perrenod Global Power Law"
)

plt.plot(
    data.index,
    data["Oscillation_1"],
    color="cyan",
    linewidth=1.5,
    label="Oscillation 1"
)

plt.plot(
    data.index,
    data["Oscillation_2"],
    color="gold",
    linewidth=1.5,
    label="Oscillation 2"
)

plt.plot(
    data.index,
    data["Combined"],
    color="limegreen",
    linewidth=2.5,
    label="Combined Model"
)

plt.plot(
    data.index,
    data["Average"],
    color="darkgreen",
    linewidth=1.8,
    linestyle="--",
    label="Average of Three Models"
)

plt.axvline(
    fit_start,
    color="gray",
    linestyle=":",
    linewidth=1
)

plt.axvline(
    fit_end,
    color="gray",
    linestyle=":",
    linewidth=1
)

plt.scatter(
    peak_date,
    peak_price,
    color="red",
    s=80,
    zorder=5
)

plt.annotate(
    f"${peak_price:,.0f}\n{peak_date.strftime('%d %B %Y')}",
    xy=(peak_date, peak_price),
    xytext=(15, 15),
    textcoords="offset points",
    fontsize=11,
    color="red",
    arrowprops=dict(arrowstyle="->", color="red")
)

plt.yscale("log")
plt.xlim(data.index.min(), data.index.max())
plt.xlabel("Date")
plt.ylabel("Bitcoin Price (USD, logarithmic scale)")
plt.title(
    "Bitcoin: Power Laws and Log-Periodic Oscillations",
    fontsize=17
)

plt.grid(True, which="both", alpha=0.25)
plt.legend(loc="upper left")
plt.tight_layout()
plt.show()

print("FIT PERIOD")
print(f"{fit_start:%d %B %Y} to {fit_end:%d %B %Y}")

print("\nOSCILLATION PARAMETERS")
print(f"Lambda 1: {LAMBDA_1}")
print(f"Omega 1: {OMEGA_1:.6f}")
print(f"Lambda 2: {LAMBDA_2}")
print(f"Omega 2: {OMEGA_2:.6f}")

print("\nFITTED COEFFICIENTS")
print(f"c0 = {c0:.6f}")
print(f"c1 = {c1:.6f}")
print(f"c2 = {c2:.6f}")
print(f"c3 = {c3:.6f}")
print(f"c4 = {c4:.6f}")

print("\nCOMBINED MODEL PEAK")
print(f"Date: {peak_date:%d %B %Y}")
print(f"Price: ${peak_price:,.0f}")