Options Lab räknar allt i webbläsaren, i JavaScript. Den här sidan är samma matematik skriven i Python, ben för ben, så att du kan läsa den, köra den och hitta felen. Python-versionen är testad mot verktygets egen kod: normalfördelningen, grekerna, binomialträdet och den implicita volatiliteten ger identiska tal ner till sista decimalen. Monte Carlo-delen är vektoriserad med numpy och drar därför andra slumptal, men följer samma scheman.
Koden hänger ihop uppifrån och ner. Klistra in blocken i ordning i en Colab-cell eller en Fabric-notebook så fungerar allt. Det enda som krävs utöver standardbiblioteket är numpy.
1. Normalfördelningen
Allt nedan vilar på den kumulativa normalfördelningen N(x). Verktyget använder Wests algoritm, en rationell approximation med ungefär femton siffrors precision, eftersom en webbläsare inte har scipy. Python-koden använder samma algoritm av ett enda skäl: att talen ska gå att jämföra exakt.
import math
import numpy as np
SQRT2PI = math.sqrt(2.0 * math.pi)
def norm_pdf(x: float) -> float:
# Klockkurvan. Allt annat på den här sidan är i praktiken fotnoter till den här raden.
return math.exp(-0.5 * x * x) / SQRT2PI
def norm_cdf(x: float) -> float:
# Wests algoritm (Hart 1968 i botten). Ja, scipy har en färdig, men verktyget på sajten
# kör i en webbläsare utan scipy, och jag vill att Python och JavaScript ger samma
# decimaler. Femton siffrors precision, vilket är fjorton fler än indatan förtjänar.
ax = abs(x)
if ax > 37.0:
c = 0.0
else:
e = math.exp(-ax * ax / 2.0)
if ax < 7.07106781186547:
b = 0.0352624965998911 * ax + 0.700383064443688
b = b * ax + 6.37396220353165
b = b * ax + 33.912866078383
b = b * ax + 112.079291497871
b = b * ax + 221.213596169931
b = b * ax + 220.206867912376
d = 0.0883883476483184 * ax + 1.75566716318264
d = d * ax + 16.064177579207
d = d * ax + 86.7807322029461
d = d * ax + 296.564248779674
d = d * ax + 637.333633378831
d = d * ax + 793.826512519948
d = d * ax + 440.413735824752
c = e * b / d
else:
# Så här långt ut i svansen räcker en kedjebråksutveckling. Här bor bara
# händelser som "aldrig kan inträffa" och som inträffar vart tionde år.
b = ax + 0.65
b = ax + 4.0 / b
b = ax + 3.0 / b
b = ax + 2.0 / b
b = ax + 1.0 / b
c = e / b / 2.506628274631
return 1.0 - c if x > 0 else c2. Prissättning: en formel, tre namn
Europeiska optioner prissätts med Black-Scholes-Merton med kontinuerlig direktavkastning q. Valutaoptioner (Garman-Kohlhagen) och optioner på terminer och räntor (Black-76) är samma formel med olika q. Det är inte en förenkling, det är så modellerna är definierade.
d1 = [ ln(S/K) + (r - q + sigma^2/2) T ] / ( sigma sqrt(T) ) d2 = d1 - sigma sqrt(T) Call = S e^(-qT) N(d1) - K e^(-rT) N(d2) Put = K e^(-rT) N(-d2) - S e^(-qT) N(-d1) Aktie/index: q = direktavkastning Valuta: q = utländsk ränta Råvara/ränta: q = r (S tolkas som terminspris)
def carry_yield(instrument: str, r: float, dividend_yield: float) -> float:
# Tre "olika" modeller, en formel. Skillnaden är vad man stoppar in som q:
# aktie/index -> q = direktavkastning (Black-Scholes-Merton)
# valuta -> q = utländsk ränta (Garman-Kohlhagen, samma formel med ny namnskylt)
# råvara/ränta -> q = r, så att S beter sig som ett terminspris (Black-76)
# Finansbranschen har alltid varit bra på att sälja samma sak tre gånger.
if instrument in ("commodity", "rate"):
return r
return dividend_yield
def d1_d2(S, K, T, r, q, sigma):
if T <= 1e-10 or sigma <= 1e-10:
return 0.0, 0.0
sqrt_t = math.sqrt(T)
d1 = (math.log(S / K) + (r - q + 0.5 * sigma * sigma) * T) / (sigma * sqrt_t)
return d1, d1 - sigma * sqrt_t
def bs_price(kind, S, K, T, r, q, sigma) -> float:
# På lösendagen finns ingen modell kvar att gömma sig bakom, bara realvärdet.
if T <= 1e-10:
return max(S - K, 0.0) if kind == "call" else max(K - S, 0.0)
d1, d2 = d1_d2(S, K, T, r, q, sigma)
if kind == "call":
return S * math.exp(-q * T) * norm_cdf(d1) - K * math.exp(-r * T) * norm_cdf(d2)
return K * math.exp(-r * T) * norm_cdf(-d2) - S * math.exp(-q * T) * norm_cdf(-d1)
def warrant_dilution(shares_outstanding: float, new_shares: float) -> float:
# Warranter ger nya aktier vid lösen, så kakan delas på fler. Faktorn n/(n+m) multipliceras
# på positionens storlek. Grovt, men bättre än att låtsas att utspädning inte finns.
if shares_outstanding > 0 and new_shares > 0:
return shares_outstanding / (shares_outstanding + new_shares)
return 1.03. Grekerna
Analytiska derivator av formeln ovan. Enheterna är valda för att gå att använda: theta per kalenderdag, vega och rho per en procentenhet. Positionens greker är summan över benen, viktad med antal och tecken.
Delta(call) = e^(-qT) N(d1) Delta(put) = e^(-qT) (N(d1) - 1) Gamma = e^(-qT) n(d1) / ( S sigma sqrt(T) ) Vega = S e^(-qT) n(d1) sqrt(T) / 100 Theta = [ -S e^(-qT) n(d1) sigma / (2 sqrt(T)) -/+ r K e^(-rT) N(+/-d2) +/- q S e^(-qT) N(+/-d1) ] / 365 Rho = +/- K T e^(-rT) N(+/-d2) / 100
def bs_greeks(kind, S, K, T, r, q, sigma) -> dict:
# Enheter, eftersom det är här folk gör bort sig:
# theta = per kalenderdag, vega = per 1 procentenhet IV, rho = per 1 procentenhet ränta.
# Läroböckerna anger per år och per 1,00. Ingen handlar en hel volatilitetsenhet åt gången.
if T <= 1e-10 or sigma <= 1e-10:
intrinsic = max(S - K, 0.0) if kind == "call" else max(K - S, 0.0)
if kind == "call":
delta = 1.0 if S > K else 0.0
else:
delta = -1.0 if S < K else 0.0
return {"price": intrinsic, "delta": delta, "gamma": 0.0, "theta": 0.0, "vega": 0.0, "rho": 0.0}
d1, d2 = d1_d2(S, K, T, r, q, sigma)
pdf = norm_pdf(d1)
disc_s, disc_k, sqrt_t = math.exp(-q * T), math.exp(-r * T), math.sqrt(T)
if kind == "call":
price = S * disc_s * norm_cdf(d1) - K * disc_k * norm_cdf(d2)
delta = disc_s * norm_cdf(d1)
rho = K * T * disc_k * norm_cdf(d2) / 100.0
else:
price = K * disc_k * norm_cdf(-d2) - S * disc_s * norm_cdf(-d1)
delta = disc_s * (norm_cdf(d1) - 1.0)
rho = -K * T * disc_k * norm_cdf(-d2) / 100.0
gamma = disc_s * pdf / (S * sigma * sqrt_t)
vega = S * disc_s * pdf * sqrt_t / 100.0
decay = -(S * disc_s * pdf * sigma) / (2.0 * sqrt_t)
if kind == "call":
theta = (decay - r * K * disc_k * norm_cdf(d2) + q * S * disc_s * norm_cdf(d1)) / 365.0
else:
theta = (decay + r * K * disc_k * norm_cdf(-d2) - q * S * disc_s * norm_cdf(-d1)) / 365.0
return {"price": price, "delta": delta, "gamma": gamma, "theta": theta, "vega": vega, "rho": rho}4. Amerikanska optioner: binomialträd
Amerikanska optioner kan lösas in i förtid och har ingen sluten formel. Verktyget använder ett Cox-Ross-Rubinstein-träd med 80 steg. I varje nod jämförs värdet av att vänta med värdet av att lösa in nu. Delta, gamma och theta läses direkt ur trädets första noder. Vega och rho fås genom att knuffa på indatan och prissätta om.
u = e^(sigma sqrt(dt)), d = 1/u, p = ( e^((r-q)dt) - d ) / ( u - d ) V(i,j) = max( lösenvärde, e^(-r dt) [ p V(i+1,j+1) + (1-p) V(i+1,j) ] )
CRR_STEPS = 80 # 80 steg: felet hamnar på någon tiondels procent och webbläsaren överlever
def crr_core(kind, S, K, T, r, q, sigma, n_steps=CRR_STEPS, american=True) -> dict:
dt = T / n_steps
u = math.exp(sigma * math.sqrt(dt))
d = 1.0 / u
disc = math.exp(-r * dt)
# Sannolikheten kan rymma utanför (0,1) när driften är enorm mot sigma*sqrt(dt).
# Att klämma fast den är fult. NaN i en prislapp är fulare.
pu = min(max((math.exp((r - q) * dt) - d) / (u - d), 0.0), 1.0)
pd = 1.0 - pu
dpu, dpd = disc * pu, disc * pd
is_call = kind == "call"
# Prisstege: ladder[k] = S * u^(k - N). Nod (i, j) ligger då på ladder[2j - i + N],
# och jag slipper ropa på pow() några tusen gånger i den inre loopen.
ladder = np.empty(2 * n_steps + 1)
ladder[n_steps] = S
for k in range(1, n_steps + 1):
ladder[n_steps + k] = ladder[n_steps + k - 1] * u
ladder[n_steps - k] = ladder[n_steps - k + 1] * d
terminal = ladder[0::2] # noderna på slutdagen: j = 0..N
values = np.maximum(terminal - K, 0.0) if is_call else np.maximum(K - terminal, 0.0)
v20 = v21 = v22 = v10 = v11 = 0.0
for i in range(n_steps - 1, -1, -1):
cont = dpu * values[1:i + 2] + dpd * values[0:i + 1]
if american:
spots = ladder[(n_steps - i)::2][: i + 1]
exercise = (spots - K) if is_call else (K - spots)
# Hela poängen med amerikanskt: i varje nod frågar man sig om det är dags att gå hem.
cont = np.maximum(cont, exercise)
values = cont
if i == 2:
v20, v21, v22 = values[0], values[1], values[2]
if i == 1:
v10, v11 = values[0], values[1]
price = float(values[0])
if n_steps < 2:
if is_call:
delta = 1.0 if S > K else 0.0
else:
delta = -1.0 if S < K else 0.0
return {"price": price, "delta": delta, "gamma": 0.0, "theta": 0.0}
su, sd, suu, sdd = S * u, S * d, S * u * u, S * d * d
delta = (v11 - v10) / (su - sd)
gamma = ((v22 - v21) / (suu - S) - (v21 - v20) / (S - sdd)) / (0.5 * (suu - sdd))
# Nod (2,1) har samma spotpris som idag, fast två steg senare. Skillnaden är rent tidsförfall.
theta = (v21 - price) / (2.0 * dt) / 365.0
return {"price": price, "delta": float(delta), "gamma": float(gamma), "theta": float(theta)}
def crr_price(kind, S, K, T, r, q, sigma, n_steps=CRR_STEPS) -> float:
if T <= 1e-10 or sigma <= 1e-10:
return bs_price(kind, S, K, T, r, q, sigma)
return crr_core(kind, S, K, T, r, q, sigma, n_steps, True)["price"]
def crr_greeks(kind, S, K, T, r, q, sigma, n_steps=CRR_STEPS) -> dict:
if T <= 1e-10 or sigma <= 1e-10:
return bs_greeks(kind, S, K, T, r, q, sigma)
base = crr_core(kind, S, K, T, r, q, sigma, n_steps, True)
# Ett träd har ingen analytisk vega. Man knuffar på indatan och ser vad som ramlar ut.
vega = crr_core(kind, S, K, T, r, q, sigma + 0.01, n_steps, True)["price"] - base["price"]
rho = (crr_core(kind, S, K, T, r + 0.0001, q, sigma, n_steps, True)["price"] - base["price"]) * 100.0
return {**base, "vega": vega, "rho": rho}5. Implicit volatilitet
Givet ett marknadspris söks det sigma som återskapar priset. För europeiska optioner används Newton-Raphson med vega som derivata, inlåst i ett intervall som halveras när Newton hoppar utanför. För amerikanska optioner används ren intervallhalvering mot trädet. Priser utanför arbitragegränserna ger inget svar alls, vilket är rätt svar.
def implied_vol(kind, S, K, T, r, q, market_price, guess=0.25, style="eu"):
if not (T > 1e-8 and market_price > 0 and S > 0 and K > 0):
return None
disc_s, disc_k = S * math.exp(-q * T), K * math.exp(-r * T)
lower = max(disc_s - disc_k, 0.0) if kind == "call" else max(disc_k - disc_s, 0.0)
upper = disc_s if kind == "call" else disc_k
tol = 1e-9 * max(1.0, market_price)
if style == "am":
# Trädets vega är för hackig för Newton, så här blir det ren intervallhalvering.
# Långsamt, tråkigt, konvergerar alltid. Som en bra revisor.
intrinsic = max(S - K, 0.0) if kind == "call" else max(K - S, 0.0)
if market_price < intrinsic - 1e-10 or market_price > upper + 1e-10:
return None
lo, hi = 1e-4, 5.0
for _ in range(60):
mid = 0.5 * (lo + hi)
p = crr_price(kind, S, K, T, r, q, mid, 60)
if abs(p - market_price) < tol:
return mid
if p > market_price:
hi = mid
else:
lo = mid
return 0.5 * (lo + hi)
# Ett pris utanför arbitragegränserna har ingen volatilitet, bara en felskrivning.
if market_price < lower - 1e-10 or market_price > upper + 1e-10:
return None
lo, hi = 1e-4, 5.0
sigma = min(max(guess or 0.25, lo), hi)
for _ in range(100):
g = bs_greeks(kind, S, K, T, r, q, sigma)
diff = g["price"] - market_price
if abs(diff) < tol:
return sigma
if diff > 0:
hi = sigma
else:
lo = sigma
vega = g["vega"] * 100.0 # bs_greeks ger vega per 1 %, Newton vill ha per 1,00
nxt = sigma - diff / vega if vega > 1e-12 else float("nan")
if not (lo < nxt < hi):
nxt = 0.5 * (lo + hi) # Newton sprang ut i skogen, halvera intervallet i stället
if abs(nxt - sigma) < 1e-12:
return nxt
sigma = nxt
return None6. Volatilitetsytan
Ytan i verktyget är en parametrisk leksak, inte en kalibrerad marknadsyta. Tre reglage flyttar volatiliteten som funktion av lösenpris och löptid. Syftet är att se hur skev och leende påverkar en position, inte att prissätta mot marknaden.
x = ln(K/S) / 0,10 sigma(K,T) = sigma0 + term (T - 30/365) + skew x + curv x^2
def surface_vol(base_sigma, S, K, T, skew=0.0, curv=0.0, term=0.0) -> float:
# Tre rattar: skew (lutning), curv (leende) och term (löptid), mätt per 10 % log-moneyness
# och runt en 30-dagars ankarpunkt. Ingen kalibrerad SVI-yta, utan en leksak för att se
# vad skev och leende gör med en position. Säljs inte som något annat.
if (skew == 0.0 and curv == 0.0 and term == 0.0) or K <= 0 or S <= 0:
return base_sigma
x = math.log(K / S) / 0.10
return max(base_sigma + term * (T - 30.0 / 365.0) + skew * x + curv * x * x, 1e-4)7. Positionen: ben, resultat och breakeven
En strategi är en lista av ben. Resultatet på slutdagen är summan av benens realvärde minus betald eller plus erhållen premie. Breakeven hittas numeriskt genom att leta teckenbyten längs prisaxeln. Om vinst eller förlust är obegränsad avgörs av nettopositionen i köpoptioner, inte av kurvans lutning i kanten av diagrammet.
def leg_payoff(leg: dict, ST: float) -> float:
# Ett ben: {"kind": "call"|"put"|"stock", "side": "long"|"short", "strike", "qty", "premium"}
# premium är det du betalade (lång) eller fick (kort) per enhet. För aktier: köpkursen.
sign = 1.0 if leg["side"] == "long" else -1.0
if leg["kind"] == "stock":
value = ST
elif leg["kind"] == "call":
value = max(ST - leg["strike"], 0.0)
else:
value = max(leg["strike"] - ST, 0.0)
return sign * leg["qty"] * (value - leg["premium"])
def strategy_payoff(legs: list, ST: float) -> float:
return sum(leg_payoff(leg, ST) for leg in legs)
def find_breakevens(payoff, low: float, high: float, n: int = 900) -> list:
# Går längs prisaxeln och letar teckenbyten, med linjär interpolation i sista steget.
step = (high - low) / n
hits, prev = [], payoff(low)
for i in range(1, n + 1):
x = low + i * step
curr = payoff(x)
if prev * curr < 0:
hits.append(round(x - step * curr / (curr - prev + 1e-15), 2))
prev = curr
return sorted(set(hits))
def analyze_extremes(legs: list, payoff, low: float, high: float) -> dict:
# Om svansarna skenar avgörs av boken, inte av att kisa på en lutning: netto långa calls
# betyder obegränsad uppsida, netto korta calls obegränsad nedsida. Aktier räknas som en
# call med lösenpris noll. Puttar tar slut vid noll, så dem fångar avsökningen.
net_calls = 0.0
for leg in legs:
if leg["kind"] in ("call", "stock"):
net_calls += (1.0 if leg["side"] == "long" else -1.0) * leg["qty"]
unlimited_profit, unlimited_loss = net_calls > 1e-9, net_calls < -1e-9
probes = [1e-6] + [low + i * (high - low) / 700 for i in range(701)] + [high * 3.0]
values = [payoff(x) for x in probes]
return {
"max_profit": math.inf if unlimited_profit else max(values),
"max_loss": -math.inf if unlimited_loss else min(values),
"unlimited_profit": unlimited_profit,
"unlimited_loss": unlimited_loss,
}8. Sannolikhet för vinst och väntevärde
Sannolikheten för vinst (POP) och väntevärdet (EV) beräknas under det riskneutrala måttet, där underliggande i snitt växer med r – q. Integralen över slutprisets lognormalfördelning tas numeriskt med trapetsregeln. Kom ihåg vad det betyder: det är marknadens prissatta sannolikheter, inte en prognos. Under det måttet tjänar en korrekt prissatt option i väntevärde bara räntan på premien, varken mer eller mindre.
ln S_T ~ N( ln S + (r - q - sigma^2/2) T , sigma^2 T ) EV = e^(-rT) * integral av resultat(S_T) * täthet POP = sannolikhetsmassan där resultat(S_T) > 0
def pop_and_ev(payoff, S, T, r, q, sigma, n: int = 1400) -> dict:
# Riskneutralt mått: driften är r - q, inte vad du hoppas på. Integralen tas över z från
# -7 till 7 med trapetsregeln. Utanför det ligger händelser med sannolikhet 1e-12, och
# har du en position som hänger på dem har du andra problem än numerisk integration.
mu = math.log(S) + (r - q - 0.5 * sigma * sigma) * T
std = sigma * math.sqrt(T)
z_lo, z_hi = -7.0, 7.0
dz = (z_hi - z_lo) / n
z_prev = z_lo
p_prev = payoff(math.exp(mu + std * z_prev))
c_prev = norm_cdf(z_prev)
pop = 0.0
ev = p_prev * norm_pdf(z_prev) * dz * 0.5
for i in range(1, n + 1):
z = z_lo + i * dz
p = payoff(math.exp(mu + std * z))
cdf = norm_cdf(z)
ev += p * norm_pdf(z) * dz * (0.5 if i == n else 1.0)
if p_prev > 0 and p > 0:
pop += cdf - c_prev
elif (p_prev > 0) != (p > 0):
# Teckenbyte mitt i ett steg: interpolera fram nollstället i stället för att
# ge hela steget till den ena sidan. Breakeven ligger sällan snällt på en gridpunkt.
cx = norm_cdf(z_prev + dz * p_prev / (p_prev - p))
pop += (cx - c_prev) if p_prev > 0 else (cdf - cx)
z_prev, p_prev, c_prev = z, p, cdf
return {"pop_pct": pop * 100.0, "ev": ev * math.exp(-r * T)}9. Monte Carlo: tre modeller för underliggande
Simuleringen drar slutpriser och värderar strategin i vart och ett. Tre processer finns. Geometrisk brownsk rörelse är Black-Scholes egen värld. Mertons hoppdiffusion lägger till plötsliga hopp med Poisson-fördelad frekvens. Heston låter volatiliteten själv vara slumpmässig och korrelerad med kursen, vilket ger skev och feta svansar. Alla tre använder antitetiska slumptal. Ur utfallen beräknas medel, median, POP, VaR 95 %, CVaR 95 % och skevhet.
GBM: dS/S = (r - q) dt + sigma dW
Merton: dS/S = (r - q - lambda kappa) dt + sigma dW + (J - 1) dN, kappa = E[J - 1]
Heston: dS/S = (r - q) dt + sqrt(v) dW1
dv = kap (theta - v) dt + xi sqrt(v) dW2, corr(dW1, dW2) = rho
def merton_params(jump_mean=-0.05, jump_vol=0.10, lam=1.0) -> dict:
# kappa är det förväntade relativa hoppet. Det dras av från driften (kompensatorn) så att
# processen fortfarande är en martingal. Annars har man uppfunnit gratispengar, igen.
jump_vol = max(jump_vol, 1e-4)
return {"lam": max(lam, 0.0), "mu_j": jump_mean, "sig_j": jump_vol,
"kappa": math.exp(jump_mean + 0.5 * jump_vol * jump_vol) - 1.0}
def heston_params(sigma, kappa=2.0, theta_vol=None, xi=0.50, rho=-0.60) -> dict:
theta_vol = sigma if theta_vol is None else theta_vol
return {"kap": max(kappa, 0.0), "theta": max(theta_vol * theta_vol, 1e-6),
"xi": max(xi, 0.0), "rho": min(max(rho, -0.99), 0.99)}
def mc_terminal_prices(model, S, T, r, q, sigma, n_sims=3000, mp=None, seed=None) -> np.ndarray:
# Antitetiska par överallt: varje slumptal z får en tvilling -z. Halverar bruset gratis,
# vilket är den enda gratislunch jag känner till i det här yrket.
rng = np.random.default_rng(seed)
half = (n_sims + 1) // 2
vol = sigma * math.sqrt(T)
if model == "gbm":
# Slutpriset under GBM har sluten form. Ingen anledning att traska 180 dagssteg
# för att landa i samma lognormalfördelning.
drift = (r - q - 0.5 * sigma * sigma) * T
z = rng.standard_normal(half)
st = np.concatenate([S * np.exp(drift + vol * z), S * np.exp(drift - vol * z)])
return st[:n_sims]
if model == "merton":
mp = mp or merton_params()
drift = (r - q - mp["lam"] * mp["kappa"] - 0.5 * sigma * sigma) * T
z = rng.standard_normal(half)
z_pair = np.concatenate([z, -z])[:n_sims]
# Summan av n normalfördelade hopp är en enda normal med medel n*mu_j och
# standardavvikelse sig_j*sqrt(n). Ett slumptal räcker, oavsett hur många hopp det blev.
n_jumps = rng.poisson(mp["lam"] * T, n_sims)
jumps = n_jumps * mp["mu_j"] + mp["sig_j"] * np.sqrt(n_jumps) * rng.standard_normal(n_sims)
return S * np.exp(drift + vol * z_pair + jumps)
if model == "heston":
mp = mp or heston_params(sigma)
n_steps = max(20, min(250, round(T * 365)))
dt = T / n_steps
sq = math.sqrt(1.0 - mp["rho"] ** 2)
z1 = rng.standard_normal((n_steps, half))
z2 = rng.standard_normal((n_steps, half))
z1 = np.concatenate([z1, -z1], axis=1)[:, :n_sims]
z2 = np.concatenate([z2, -z2], axis=1)[:, :n_sims]
x = np.full(n_sims, math.log(S))
v = np.full(n_sims, sigma * sigma)
for t in range(n_steps):
# Euler med "full truncation": variansen får bli negativ i bokföringen men klipps
# till noll där den används. Eulerschemat respekterar inte att varians är positiv,
# så någon måste vara vuxen i rummet.
vp = np.maximum(v, 0.0)
sv = np.sqrt(vp * dt)
x = x + (r - q - 0.5 * vp) * dt + sv * z1[t]
v = v + mp["kap"] * (mp["theta"] - vp) * dt + mp["xi"] * sv * (mp["rho"] * z1[t] + sq * z2[t])
return np.exp(x)
raise ValueError(f"okänd modell: {model}")
def mc_summary(payoff, terminal_prices: np.ndarray) -> dict:
pnl = np.sort(np.array([payoff(float(st)) for st in terminal_prices]))
n = len(pnl)
tail = max(1, n // 20) # de sämsta fem procenten
mean, std = float(pnl.mean()), float(pnl.std())
return {
"mean": mean,
"median": float(pnl[n // 2]),
"pop_pct": float((pnl > 0).mean() * 100.0),
"var95": float(pnl[tail - 1]), # dörren
"cvar95": float(pnl[:tail].mean()), # hur långt ner det är till golvet bakom dörren
"skew": float((((pnl - mean) / std) ** 3).mean()) if std > 0 else 0.0,
}10. Attribution: varför tjänade eller förlorade positionen?
Resultatförändringen delas upp på två sätt. Taylor-uppdelningen multiplicerar varje grek med sin förändring och lämnar en rest, som består av korstermer som vanna, charm och volga. Den sekventiella uppdelningen prissätter om positionen exakt med en faktor i taget: spot, volatilitet, tid, ränta. Delarna summerar då exakt till helheten, men fördelningen beror på ordningen. Att de två metoderna inte ger samma svar är informationen, inte ett fel.
dV ~ Delta dS + 1/2 Gamma dS^2 + Theta dt + Vega dsigma + Rho dr + rest
def attribution(value_fn, greeks: dict, S, d_spot, d_iv, d_days, d_rate) -> dict:
# value_fn(spot, dagar_framåt, iv_skift, ränte_skift) ger positionens teoretiska värde.
# Två sätt att förklara samma resultat, och de ska inte stämma överens:
base = value_fn(S, 0.0, 0.0, 0.0)
full = value_fn(S + d_spot, d_days, d_iv, d_rate) - base
# (a) Taylor: grekerna gånger sina knuffar, gamma till andra ordningen. Det som blir över
# är korstermer (vanna, charm, volga) och allt annat Taylor aldrig fick se.
taylor = {
"delta": greeks["delta"] * d_spot,
"gamma": 0.5 * greeks["gamma"] * d_spot * d_spot,
"theta": greeks["theta"] * d_days,
"vega": greeks["vega"] * d_iv * 100.0,
"rho": greeks["rho"] * d_rate * 100.0,
}
taylor["residual"] = full - sum(taylor.values())
# (b) Sekventiell exakt omprisning: en faktor i taget, spot -> vol -> tid -> ränta.
# Ordningen spelar roll (det är poängen) och delarna summerar till helheten per konstruktion.
s1 = value_fn(S + d_spot, 0.0, 0.0, 0.0)
s2 = value_fn(S + d_spot, 0.0, d_iv, 0.0)
s3 = value_fn(S + d_spot, d_days, d_iv, 0.0)
sequential = {"spot": s1 - base, "vega": s2 - s1, "theta": s3 - s2, "rho": full - (s3 - base)}
return {"total": full, "taylor": taylor, "sequential": sequential}11. Positionsstorlek: egen drift, Kelly och säkerhetskrav
Här lämnas det riskneutrala måttet. Du anger en egen förväntad drift, och fördelningen räknas om under den. Kapital i risk är maxförlusten när den är definierad, annars CVaR 95 % som ersättning. Kelly-andelen beräknas över hela utfallsfördelningen, inte med formeln för ett enkelt vad. Säkerhetskravet är en schablon i Reg T-stil och motsvarar inte någon specifik mäklares regler.
f* = argmax över f av E[ ln(1 + f x) ], x = resultat / kapital i risk Löses ur E[ x / (1 + f x) ] = 0 med intervallhalvering
def real_world_stats(payoff, S, T, mu, q, sigma, n: int = 1200) -> dict:
# Här byts r mot din egen driftgissning mu. Det är P-måttet: vad du tror ska hända, inte
# vad marknaden prisar. Skräp in, skräp ut, fast med fler decimaler.
m = math.log(S) + (mu - q - 0.5 * sigma * sigma) * T
std = sigma * math.sqrt(T)
z_lo, z_hi = -7.0, 7.0
dz = (z_hi - z_lo) / n
pts, ev, pop = [], 0.0, 0.0
for i in range(n + 1):
z = z_lo + i * dz
w = norm_pdf(z) * dz * (0.5 if i in (0, n) else 1.0)
p = payoff(math.exp(m + std * z))
pts.append((p, w))
ev += p * w
if p > 0:
pop += w
# VaR och CVaR 95 %: sortera utfallen och ät sannolikhetsmassa underifrån tills 5 % är slut.
acc, c_sum, c_w = 0.0, 0.0, 0.0
ordered = sorted(pts, key=lambda o: o[0])
var95 = ordered[0][0]
for p, w in ordered:
take = min(w, 0.05 - acc)
if take <= 0:
break
c_sum += p * take
c_w += take
acc += take
var95 = p
return {"ev": ev, "pop_pct": pop * 100.0, "var95": var95,
"cvar95": c_sum / c_w if c_w > 0 else var95, "pts": pts}
def kelly_fraction(pts: list, stake: float) -> float:
# Generaliserad Kelly över hela utfallsfördelningen: maximera E[log(1 + f*x)] där x är
# avkastning per satsad krona. Derivatan g(f) = E[x / (1 + f*x)] är avtagande, så
# intervallhalvering hittar nollstället. f kapas precis innan värsta utfallet ger konkurs,
# eftersom log(0) är matematikens sätt att säga "gör inte så".
if not stake > 0:
return 0.0
xs = [(p / stake, w) for p, w in pts]
worst = min(x for x, _ in xs)
f_max = 0.999 / (-worst) if worst < 0 else 10.0
def g(f):
return sum(w * x / (1.0 + f * x) for x, w in xs)
if g(0.0) <= 0:
return 0.0 # ingen positiv förväntan: rätt insats är noll, hur kul strategin än ser ut
if g(f_max) > 0:
return f_max
lo, hi = 0.0, f_max
for _ in range(60):
mid = 0.5 * (lo + hi)
if g(mid) > 0:
lo = mid
else:
hi = mid
return 0.5 * (lo + hi)
def estimate_margin(legs: list, S: float, max_loss: float, unlimited_loss: bool) -> float:
# Tumregler i Reg T-stil. Din mäklare räknar annorlunda och har alltid rätt, framför allt
# klockan 03 en natt när volatiliteten just har dubblats.
if not unlimited_loss:
return max(-max_loss, 0.0) # definierad risk: kravet är maxförlusten
total = 0.0
for leg in legs:
qty = abs(leg["qty"])
if leg["kind"] == "stock":
total += 0.5 * S * qty
elif leg["side"] == "long":
total += leg["premium"] * qty
elif leg["kind"] == "call":
otm = max(leg["strike"] - S, 0.0)
total += max(0.2 * S - otm + leg["premium"], 0.1 * S + leg["premium"]) * qty
else:
otm = max(S - leg["strike"], 0.0)
total += max(0.2 * S - otm + leg["premium"], 0.1 * leg["strike"] + leg["premium"]) * qty
return total
def position_sizing(payoff, legs, S, T, mu, q, sigma, account, risk_budget_pct, contract_mult=100.0) -> dict:
ext = analyze_extremes(legs, payoff, 0.3 * S, 1.7 * S)
rw = real_world_stats(payoff, S, T, mu, q, sigma)
# Kapital i risk per enhet: maxförlusten om den finns, annars CVaR 95 % som nödlösning.
if ext["unlimited_loss"]:
stake_unit = max(-rw["cvar95"], 0.0)
else:
stake_unit = max(-ext["max_loss"], 0.0)
stake_pos = stake_unit * contract_mult
f_star = kelly_fraction(rw["pts"], stake_unit)
budget = account * risk_budget_pct / 100.0
return {
"stake_per_position": stake_pos,
"n_by_budget": math.floor(budget / stake_pos) if stake_pos > 0 else None,
"kelly_fraction": f_star,
"n_by_kelly": math.floor(f_star * account / stake_pos) if stake_pos > 0 else None,
# Halv-Kelly: tre fjärdedelar av tillväxten, halva magsåret. Hel Kelly förutsätter att
# du kan dina sannolikheter exakt. Det kan du inte.
"n_by_half_kelly": math.floor(0.5 * f_star * account / stake_pos) if stake_pos > 0 else None,
"margin_per_unit": estimate_margin(legs, S, ext["max_loss"], ext["unlimited_loss"]),
"ev_real_world": rw["ev"],
"pop_real_world_pct": rw["pop_pct"],
}Exempel: prova själv
En köpt call, prissatt, analyserad och simulerad med koden ovan.
S, K, T, r, q, sigma = 100.0, 100.0, 1.0, 0.05, 0.0, 0.20
premie = bs_price("call", S, K, T, r, q, sigma) # 10,4506, facit i varenda lärobok
ben = [{"kind": "call", "side": "long", "strike": K, "qty": 1, "premium": premie}]
resultat = lambda ST: strategy_payoff(ben, ST)
print(bs_greeks("call", S, K, T, r, q, sigma))
print(crr_price("put", S, K, T, r, q, sigma)) # amerikansk put, cirka 6,08
print(find_breakevens(resultat, 30, 170)) # [110.45]
print(pop_and_ev(resultat, S, T, r, q, sigma)) # POP cirka 36 %, EV cirka 0,5: räntan på premien, inget mer
for modell in ("gbm", "merton", "heston"):
slutpriser = mc_terminal_prices(modell, S, T, r, q, sigma, n_sims=20000, seed=7)
print(modell, mc_summary(resultat, slutpriser))
# Och så den obekväma frågan: hur mycket ska man satsa om man tror på 8 % drift?
print(position_sizing(resultat, ben, S, T, mu=0.08, q=q, sigma=sigma,
account=1_000_000, risk_budget_pct=2.0))Svagheter i modellerna
- Konstant volatilitet. Black-Scholes och trädet antar en enda sigma. Marknaden har en hel yta, och den rör sig. Verktygets yta är parametrisk och okalibrerad.
- Diskreta utdelningar saknas. Direktavkastningen är kontinuerlig. För amerikanska köpoptioner runt en utdelningsdag ger det fel tidpunkt för förtida lösen.
- Trädet har 80 steg. Priset svänger något med antalet steg, och gamma och theta ur trädet är hackigare än de analytiska.
- Heston med Euler. Schemat har diskretiseringsfel, särskilt när Feller-villkoret bryts och variansen ofta slår i noll. Parametrarna sätts för hand, de är inte kalibrerade mot optionspriser.
- Hopp och stokastisk volatilitet gäller bara simuleringen. Premierna prissätts fortfarande med Black-Scholes eller trädet. Simulerar du med Heston eller Merton mot BSM-premier är modellerna inte inbördes konsistenta, och väntevärdet speglar den skillnaden.
- Riskneutralt är inte verkligt. POP och EV under Q-måttet säger vad marknaden prisar, inte vad som kommer att hända. Under din egen drift säger de vad du tror, vilket är ännu mindre värt.
- Kelly förutsätter känd fördelning. Med skattade sannolikheter överskattar hel Kelly systematiskt rätt insats. Halv Kelly finns med av det skälet.
- Säkerhetskravet är en schablon. Verkliga krav beror på mäklare, portföljmarginal och stressnivå, och de höjs när du minst vill det.
- Inga transaktionskostnader. Ingen spread, inget courtage, ingen likviditetsbrist, ingen tilldelningsrisk.
Källor
Black och Scholes (1973), The Pricing of Options and Corporate Liabilities. Merton (1973), Theory of Rational Option Pricing. Black (1976), The Pricing of Commodity Contracts. Garman och Kohlhagen (1983), Foreign Currency Option Values. Cox, Ross och Rubinstein (1979), Option Pricing: A Simplified Approach. Merton (1976), Option Pricing When Underlying Stock Returns Are Discontinuous. Heston (1993), A Closed-Form Solution for Options with Stochastic Volatility. Lord, Koekkoek och van Dijk (2010), A Comparison of Biased Simulation Schemes for Stochastic Volatility Models (full truncation). Kelly (1956), A New Interpretation of Information Rate. West (2005), Better Approximations to Cumulative Normal Functions. Implementationen är egen.
Ansvarsfriskrivning
Koden och verktyget är skrivna i utbildningssyfte för att illustrera optionsteori. De är inte investeringsrådgivning, inte en rekommendation och inte ett beslutsunderlag för handel. Modellerna är förenklingar, koden kan innehålla fel och ingenting garanteras vara korrekt, fullständigt eller lämpligt för något ändamål. Koden tillhandahålls i befintligt skick utan garantier av något slag. Handel med optioner och andra derivat innebär hög risk, och utfärdade optioner kan ge förluster som överstiger insatsen. Alla beslut du fattar är dina egna, och skribenten tar inget ansvar för förluster eller andra konsekvenser av att någon använder koden, verktyget eller texten. Prata med en licensierad rådgivare innan du handlar.