Skip to content

Instantly share code, notes, and snippets.

@MuslemRahimi
Created September 11, 2025 12:18
Show Gist options
  • Select an option

  • Save MuslemRahimi/3f07b7322e65749fceec5f37ab1aabe0 to your computer and use it in GitHub Desktop.

Select an option

Save MuslemRahimi/3f07b7322e65749fceec5f37ab1aabe0 to your computer and use it in GitHub Desktop.
import os
from datetime import datetime, timedelta
import orjson
from typing import List, Tuple, Optional
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import yfinance as yf
try:
from tqdm import tqdm
except Exception:
tqdm = lambda x: x
def read_price_history(symbol: str) -> Optional[pd.DataFrame]:
try:
ticker = yf.Ticker(symbol)
df = ticker.history(period="2y")
if df.empty:
return None
df = df.reset_index()
df = df.rename(columns={"Date": "date", "Close": "adjClose"})
df = df[["date", "adjClose"]]
df = df.sort_values("date").reset_index(drop=True)
return df
except Exception as e:
print("Failed reading price history:", e)
return None
def calculate_rolling_var_history(df: pd.DataFrame, window_days: int = 252, step_days: int = 21,
num_simulations: int = 10000, trading_days_per_month: int = 21,
confidence_level: float = 0.95) -> List[dict]:
history = []
for i in range(window_days, len(df) + 1, step_days):
window_df = df.iloc[i-window_days:i].copy()
var_pct, _ = compute_var_monte_carlo(window_df, num_simulations, trading_days_per_month, confidence_level)
history.append({
"date": window_df.iloc[-1]["date"].strftime("%Y-%m-%d"),
"var": var_pct
})
return history
def compute_var_monte_carlo(df: pd.DataFrame, num_simulations: int = 10000, trading_days_per_month: int = 21,
confidence_level: float = 0.95) -> Tuple[float, np.ndarray]:
df = df.copy()
if "adjClose" not in df.columns:
raise ValueError("DataFrame must contain adjClose column")
df["Returns"] = df["adjClose"].pct_change()
df = df.dropna()
if len(df) < 2:
return 0.0, np.array([])
daily_returns = df["Returns"].values
mean_return = np.mean(daily_returns)
std_return = np.std(daily_returns, ddof=0)
rng = np.random.default_rng()
monthly_returns = rng.normal(loc=mean_return, scale=std_return, size=(num_simulations, trading_days_per_month))
monthly_compounded = np.prod(1 + monthly_returns, axis=1) - 1
var_percentile = (1 - confidence_level) * 100
monthly_var = np.percentile(monthly_compounded, var_percentile)
var_percentage = round(abs(monthly_var) * 100, 2)
if var_percentage > 99:
var_percentage = 99.0
return var_percentage, monthly_compounded
def plot_monte_carlo_distribution(monthly_returns: np.ndarray, confidence_level: float = 0.95, symbol: str = "SYMB"):
if monthly_returns.size == 0:
print("No simulated returns to plot")
return
var_pctile = (1 - confidence_level) * 100
var_value = np.percentile(monthly_returns, var_pctile)
plt.figure(figsize=(10, 5))
plt.hist(monthly_returns, bins=60, density=True, alpha=0.75)
plt.axvline(var_value, linestyle="--", linewidth=2, label=f"VaR {int(confidence_level*100)}% = {var_value:.2%}")
plt.title(f"Monte Carlo Simulated Monthly Returns - Distribution ({symbol})")
plt.xlabel("Monthly Return")
plt.ylabel("Density")
plt.legend()
plt.grid(alpha=0.25)
plt.tight_layout()
plt.savefig(f"{symbol}_monte_carlo_distribution.png")
plt.close()
def plot_monte_carlo_paths(df: pd.DataFrame, num_simulations: int = 100, horizon_days: int = 21, symbol: str = "SYMB"):
df = df.copy()
df["Returns"] = df["adjClose"].pct_change()
df = df.dropna()
if df.empty:
print("Not enough data for paths plot")
return
mean_return = np.mean(df["Returns"].values)
std_return = np.std(df["Returns"].values, ddof=0)
last_price = float(df["adjClose"].iloc[-1])
rng = np.random.default_rng()
plt.figure(figsize=(10, 5))
x = np.arange(horizon_days + 1)
for i in range(num_simulations):
simulated = rng.normal(loc=mean_return, scale=std_return, size=horizon_days)
prices = [last_price]
for r in simulated:
prices.append(prices[-1] * (1 + r))
plt.plot(x, prices, alpha=0.08)
plt.title(f"Monte Carlo Price Paths ({symbol}) - next {horizon_days} days")
plt.xlabel("Days ahead")
plt.ylabel("Price")
plt.grid(alpha=0.25)
plt.tight_layout()
plt.savefig(f"{symbol}_monte_carlo_paths.png")
plt.close()
def plot_var_history(history: List[dict], symbol: str = "SYMB"):
if not history:
print("Empty VaR history")
return
df = pd.DataFrame(history)
df["date"] = pd.to_datetime(df["date"])
df = df.sort_values("date")
plt.figure(figsize=(10, 4))
plt.plot(df["date"], df["var"], marker="o")
plt.title(f"VaR (95%) over time - {symbol}")
plt.xlabel("Date")
plt.ylabel("VaR (%)")
plt.grid(alpha=0.25)
plt.tight_layout()
plt.savefig(f"{symbol}_var_history.png")
plt.close()
def plot_combined(df: pd.DataFrame, monthly_returns: np.ndarray, history: List[dict], symbol: str = "SYMB"):
import matplotlib.dates as mdates
df_local = df.copy()
df_local["Returns"] = df_local["adjClose"].pct_change()
df_local = df_local.dropna()
last_price = float(df_local["adjClose"].iloc[-1]) if not df_local.empty else np.nan
fig, axes = plt.subplots(2, 2, figsize=(14, 9))
ax = axes[0, 0]
if monthly_returns.size > 0:
ax.hist(monthly_returns, bins=60, density=True, alpha=0.75)
var_value = np.percentile(monthly_returns, 5)
ax.axvline(var_value, linestyle="--", linewidth=2, label=f"VaR 95% = {var_value:.2%}")
ax.set_title("Simulated Monthly Returns")
ax.set_xlabel("Monthly return")
ax.set_ylabel("Density")
ax.legend()
ax.grid(alpha=0.25)
ax = axes[0, 1]
mean_return = np.mean(df_local["Returns"].values) if not df_local.empty else 0
std_return = np.std(df_local["Returns"].values) if not df_local.empty else 0
rng = np.random.default_rng()
horizon = 21
x = np.arange(horizon + 1)
for i in range(100):
simulated = rng.normal(loc=mean_return, scale=std_return, size=horizon)
prices = [last_price]
for r in simulated:
prices.append(prices[-1] * (1 + r))
ax.plot(x, prices, alpha=0.06)
ax.set_title(f"Simulated Price Paths (next {horizon} days)")
ax.set_xlabel("Days ahead")
ax.set_ylabel("Price")
ax.grid(alpha=0.25)
ax = axes[1, 0]
if not df_local.empty:
rolling_vol = df_local["Returns"].rolling(window=21).std() * np.sqrt(21)
ax.plot(df_local["date"].iloc[-len(rolling_vol.dropna()):], rolling_vol.dropna())
ax.set_title("21-day Rolling Volatility (annualized approx)")
ax.set_xlabel("Date")
ax.set_ylabel("Volatility (approx)")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
fig.autofmt_xdate(rotation=30)
else:
ax.text(0.5, 0.5, "Not enough data", ha="center")
ax = axes[1, 1]
if history:
hdf = pd.DataFrame(history)
hdf["date"] = pd.to_datetime(hdf["date"])
hdf = hdf.sort_values("date")
ax.plot(hdf["date"], hdf["var"], marker="o")
ax.set_title("VaR History")
ax.set_xlabel("Date")
ax.set_ylabel("VaR (%)")
ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
fig.autofmt_xdate(rotation=30)
else:
ax.text(0.5, 0.5, "No VaR history", ha="center")
plt.tight_layout()
plt.savefig(f"{symbol}_combined_analysis.png")
plt.close()
print(f"Plots saved as {symbol}_*.png")
def build_var_history_entry(var_pct: float) -> dict:
return {"date": datetime.today().strftime("%Y-%m-%d"), "var": var_pct}
def process_symbol(symbol: str, recent_days: int = 252,
num_simulations: int = 10000, trading_days_per_month: int = 21,
horizon_days: int = 21) -> None:
df = read_price_history(symbol)
if df is None or df.empty:
print(f"No data available for {symbol}")
return
if recent_days and len(df) > recent_days:
df_recent = df.tail(recent_days).reset_index(drop=True)
else:
df_recent = df.copy()
var_pct, monthly_returns = compute_var_monte_carlo(df_recent, num_simulations=num_simulations,
trading_days_per_month=trading_days_per_month)
print(f"{symbol} - Current Monte Carlo monthly VaR (95%): {var_pct}%")
history = calculate_rolling_var_history(df, window_days=252, step_days=30,
num_simulations=1000, trading_days_per_month=trading_days_per_month)
if history:
print(f"Historical VaR calculated for {len(history)} periods")
print(f"VaR range: {min(h['var'] for h in history):.2f}% - {max(h['var'] for h in history):.2f}%")
plot_monte_carlo_distribution(monthly_returns, symbol=symbol)
plot_monte_carlo_paths(df_recent, num_simulations=200, horizon_days=horizon_days, symbol=symbol)
plot_var_history(history, symbol=symbol)
plot_combined(df_recent, monthly_returns, history, symbol=symbol)
def main():
simulations = 100_000
recent_days = 252
symbols = ["AVGO"]
horizon = 30
for symbol in tqdm(symbols):
process_symbol(symbol, recent_days=recent_days, num_simulations=simulations, horizon_days=horizon)
if __name__ == "__main__":
main()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment