Questioning the Lognormal Assumption in BSM¶

import numpy as np
import pandas as pd
import yfinance as yf

ticker = "TSLA"
start_date = "2020-01-01"
end_date = "2025-01-01"

stock = yf.Ticker(ticker)

# 가격 데이터 다운로드
price_df = yf.download(
    ticker,
    start=start_date,
    end=end_date,
    auto_adjust=False,
    progress=False
)

# 배당 데이터 확인
dividends = stock.dividends

print("Price data head:")
print(price_df.head())
print("\nDividends:")
print(dividends.tail())

# 종가 기준 로그수익률 계산
price_df["Log Return"] = np.log(price_df["Adj Close"] / price_df["Adj Close"].shift(1))

print("\nLog return head:")
print(price_df[["Adj Close", "Log Return"]].dropna().head())

total_dividends = dividends.loc[start_date:end_date].sum()

print(f"Total dividends during sample period: {total_dividends}")

if total_dividends == 0:
    print(f"{ticker} can be treated as a non-dividend-paying stock in this sample period.")
else:
    print(f"{ticker} paid dividends during this sample period. Consider choosing another stock.")

import numpy as np
import pandas as pd
import yfinance as yf
import scipy.stats as st
import matplotlib.pyplot as plt

ticker = "TSLA"
start_date = "2020-01-01"
end_date = "2025-01-01"

# 1. 가격 데이터 다운로드
raw_df = yf.download(
    ticker,
    start=start_date,
    end=end_date,
    auto_adjust=False,
    progress=False
)
# 2. TSLA만 따로 데이터프레임으로 정리
# yfinance 버전에 따라 컬럼 구조가 달라질 수 있어서 분기 처리
if isinstance(raw_df.columns, pd.MultiIndex):
    TSLA = raw_df.xs(ticker, axis=1, level=-1).copy()
else:
    TSLA = raw_df.copy()
# 컬럼 이름 확인용
print("TSLA columns:")
print(TSLA.columns)
# 3. 배당 확인
stock = yf.Ticker(ticker)
dividends = stock.dividends.loc[start_date:end_date]
total_dividends = dividends.sum()
print(f"\nTotal dividends during sample period: {total_dividends}")
if total_dividends == 0:
    print(f"{ticker} can be treated as a non-dividend-paying stock in this sample period.")
else:
    print(f"{ticker} paid dividends during this sample period.")
# 4. 로그수익률 계산
TSLA["Log Return"] = np.log(TSLA["Adj Close"] / TSLA["Adj Close"].shift(1))
# 5. 표본을 절반으로 분할 (시간순)
split_idx = len(TSLA) // 2
TSLA_first = TSLA.iloc[:split_idx].copy()
TSLA_second = TSLA.iloc[split_idx:].copy()
print("\nSample split:")
print("Full sample size:", len(TSLA))
print("First half size:", len(TSLA_first))
print("Second half size:", len(TSLA_second))
print("First half period:", TSLA_first.index.min().date(), "to", TSLA_first.index.max().date())
print("Second half period:", TSLA_second.index.min().date(), "to", TSLA_second.index.max().date())
# 6. 앞 절반의 로그수익률 추출
r = TSLA_first["Log Return"].dropna().values

print("\nFirst-half log return summary:")
print(pd.Series(r).describe())



def fit_normal(x):
    params = st.norm.fit(x)   # (mu, sigma)
    loglik = np.sum(st.norm.logpdf(x, *params))
    k = len(params)
    return {
        "dist": "Normal",
        "params": params,
        "loglik": loglik,
        "aic": 2 * k - 2 * loglik,
        "bic": k * np.log(len(x)) - 2 * loglik
    }

def fit_t(x):
    params = st.t.fit(x)      # (df, loc, scale)
    loglik = np.sum(st.t.logpdf(x, *params))
    k = len(params)
    return {
        "dist": "Student-t",
        "params": params,
        "loglik": loglik,
        "aic": 2 * k - 2 * loglik,
        "bic": k * np.log(len(x)) - 2 * loglik
    }

def fit_laplace(x):
    params = st.laplace.fit(x)   # (loc, scale)
    loglik = np.sum(st.laplace.logpdf(x, *params))
    k = len(params)
    return {
        "dist": "Laplace",
        "params": params,
        "loglik": loglik,
        "aic": 2 * k - 2 * loglik,
        "bic": k * np.log(len(x)) - 2 * loglik
    }

fit_results = [
    fit_normal(r),
    fit_t(r),
    fit_laplace(r)
]

fit_table = pd.DataFrame(fit_results).sort_values("aic").reset_index(drop=True)

print("\nDistribution fit comparison (first-half log returns):")
print(fit_table[["dist", "loglik", "aic", "bic"]])

best_dist_name = fit_table.loc[0, "dist"]
print(f"\nBest fit by AIC: {best_dist_name}")

# -------------------------------------------------------
# KS test
# -------------------------------------------------------
ks_results = []

for res in fit_results:
    if res["dist"] == "Normal":
        ks_stat, ks_p = st.kstest(r, "norm", args=res["params"])
    elif res["dist"] == "Student-t":
        ks_stat, ks_p = st.kstest(r, "t", args=res["params"])
    elif res["dist"] == "Laplace":
        ks_stat, ks_p = st.kstest(r, "laplace", args=res["params"])
    else:
        ks_stat, ks_p = np.nan, np.nan

    ks_results.append({
        "dist": res["dist"],
        "KS stat": ks_stat,
        "KS p-value": ks_p
    })

ks_table = pd.DataFrame(ks_results).sort_values("KS stat").reset_index(drop=True)

print("\nKS test comparison:")
print(ks_table)

# -------------------------------------------------------
# 정규성 자체 확인
# -------------------------------------------------------
jb_stat, jb_p = st.jarque_bera(r)
sw_stat, sw_p = st.shapiro(r)

print("\nNormality tests for first-half log returns:")
print(f"Jarque-Bera stat = {jb_stat:.4f}, p-value = {jb_p:.6f}")
print(f"Shapiro-Wilk stat = {sw_stat:.4f}, p-value = {sw_p:.6f}")

if jb_p < 0.05:
    print("Jarque-Bera rejects normality at the 5% level.")
else:
    print("Jarque-Bera does not reject normality at the 5% level.")

if sw_p < 0.05:
    print("Shapiro-Wilk rejects normality at the 5% level.")
else:
    print("Shapiro-Wilk does not reject normality at the 5% level.")

# -------------------------------------------------------
# 히스토그램 + 적합 분포
# -------------------------------------------------------
x_grid = np.linspace(r.min(), r.max(), 500)

plt.figure(figsize=(10, 6))
plt.hist(r, bins=40, density=True, alpha=0.6, edgecolor="black", label="First-half log returns")

# Normal
norm_params = [res["params"] for res in fit_results if res["dist"] == "Normal"][0]
plt.plot(x_grid, st.norm.pdf(x_grid, *norm_params), linewidth=2, label="Normal fit")

# Student-t
t_params = [res["params"] for res in fit_results if res["dist"] == "Student-t"][0]
plt.plot(x_grid, st.t.pdf(x_grid, *t_params), linewidth=2, label="Student-t fit")

# Laplace
lap_params = [res["params"] for res in fit_results if res["dist"] == "Laplace"][0]
plt.plot(x_grid, st.laplace.pdf(x_grid, *lap_params), linewidth=2, label="Laplace fit")

plt.title("First-half TSLA Log Returns: Histogram with Fitted Distributions")
plt.xlabel("Log Return")
plt.ylabel("Density")
plt.legend()
plt.show()

# -------------------------------------------------------
# QQ plot: Normal
# -------------------------------------------------------
plt.figure(figsize=(6, 6))
st.probplot(r, dist="norm", sparams=norm_params, plot=plt)
plt.title("QQ Plot vs Normal")
plt.show()

# -------------------------------------------------------
# QQ plot: Student-t
# -------------------------------------------------------
plt.figure(figsize=(6, 6))
st.probplot(r, dist=st.t, sparams=t_params, plot=plt)
plt.title("QQ Plot vs Student-t")
plt.show()

# -------------------------------------------------------
# QQ plot: Laplace
# -------------------------------------------------------

plt.figure(figsize=(6,6))

st.probplot(
    r,
    dist=st.laplace,
    sparams=lap_params,
    plot=plt
)

plt.title("QQ Plot vs Laplace")
plt.show()
# -------------------------------------------------------
# 근거 자동 요약
# -------------------------------------------------------
best_aic = fit_table.loc[0, "dist"]
best_ks = ks_table.loc[0, "dist"]

print("\nInterpretation:")
print(f"1. By AIC, the best-fitting distribution is: {best_aic}")
print(f"2. By KS statistic, the closest distribution is: {best_ks}")

if jb_p < 0.05 and sw_p < 0.05:
    print("3. Both Jarque-Bera and Shapiro-Wilk reject normality, suggesting that the normal assumption is weak.")
else:
    print("3. Normality is not strongly rejected by all tests, so the evidence against normality is limited.")

if best_aic != "Normal" or best_ks != "Normal":
    print("4. Since at least one goodness-of-fit criterion prefers a non-normal distribution, the BSM lognormal assumption may be questionable in this subsample.")
else:
    print("4. Since normal distribution remains competitive, the BSM lognormal assumption is not strongly contradicted in this subsample.")

r2 = TSLA_second["Log Return"].dropna().values

print(pd.Series(r2).describe())

fit_results2 = [
    fit_normal(r2),
    fit_t(r2),
    fit_laplace(r2)
]

fit_table2 = pd.DataFrame(fit_results2).sort_values("aic").reset_index(drop=True)

print("\nDistribution fit comparison (second-half log returns):")
print(fit_table2[["dist", "loglik", "aic", "bic"]])

ks_results2 = []

for res in fit_results2:
    if res["dist"] == "Normal":
        ks_stat, ks_p = st.kstest(r2, "norm", args=res["params"])
    elif res["dist"] == "Student-t":
        ks_stat, ks_p = st.kstest(r2, "t", args=res["params"])
    elif res["dist"] == "Laplace":
        ks_stat, ks_p = st.kstest(r2, "laplace", args=res["params"])

    ks_results2.append({
        "dist": res["dist"],
        "KS stat": ks_stat,
        "KS p-value": ks_p
    })

ks_table2 = pd.DataFrame(ks_results2).sort_values("KS stat").reset_index(drop=True)

print("\nKS test comparison (second half):")
print(ks_table2)

jb_stat2, jb_p2 = st.jarque_bera(r2)
sw_stat2, sw_p2 = st.shapiro(r2)

print("\nNormality tests (second-half log returns):")
print(f"Jarque-Bera stat = {jb_stat2:.4f}, p-value = {jb_p2:.6f}")
print(f"Shapiro-Wilk stat = {sw_stat2:.4f}, p-value = {sw_p2:.6f}")

if jb_p2 < 0.05:
    print("Jarque-Bera rejects normality at the 5% level.")
else:
    print("Jarque-Bera does not reject normality at the 5% level.")

if sw_p2 < 0.05:
    print("Shapiro-Wilk rejects normality at the 5% level.")
else:
    print("Shapiro-Wilk does not reject normality at the 5% level.")

x_grid = np.linspace(r2.min(), r2.max(), 500)

plt.figure(figsize=(10,6))
plt.hist(r2, bins=40, density=True, alpha=0.6, edgecolor="black")

norm_params2 = [res["params"] for res in fit_results2 if res["dist"] == "Normal"][0]
plt.plot(x_grid, st.norm.pdf(x_grid, *norm_params2), linewidth=2, label="Normal")

t_params2 = [res["params"] for res in fit_results2 if res["dist"] == "Student-t"][0]
plt.plot(x_grid, st.t.pdf(x_grid, *t_params2), linewidth=2, label="Student-t")

lap_params2 = [res["params"] for res in fit_results2 if res["dist"] == "Laplace"][0]
plt.plot(x_grid, st.laplace.pdf(x_grid, *lap_params2), linewidth=2, label="Laplace")

plt.legend()
plt.title("Second-half TSLA Log Returns Distribution Fit")
plt.show()

plt.figure(figsize=(6,6))
st.probplot(r2, dist="norm", sparams=norm_params2, plot=plt)
plt.title("QQ Plot vs Normal")
plt.show()

plt.figure(figsize=(6,6))
st.probplot(r2, dist=st.laplace, sparams=lap_params2, plot=plt)
plt.title("QQ Plot vs Laplace")
plt.show()