1+ # Black-Scholes Option Pricing Algorithm in R
2+ # Implements the Black-Scholes-Merton model for European option pricing
3+ # Features: Call/Put pricing, Greeks calculation, and implied volatility estimation
4+
5+ library(R6 )
6+
7+ # ' BlackScholesCalculator Class
8+ # ' @description R6 class for option pricing using Black-Scholes model
9+ # ' @details Calculates option prices and Greeks for European options
10+ # ' Assumptions:
11+ # ' - No dividend payments
12+ # ' - European-style options (can only be exercised at expiration)
13+ # ' - Log-normal distribution of stock prices
14+ # ' - Constant risk-free rate and volatility
15+ # ' - No transaction costs or taxes
16+ # ' - Perfectly divisible securities
17+ BlackScholesCalculator <- R6Class(
18+ " BlackScholesCalculator" ,
19+
20+ public = list (
21+ # ' @description Initialize calculator with market parameters
22+ # ' @param r Risk-free interest rate (annualized)
23+ # ' @param include_checks Whether to perform parameter validation
24+ initialize = function (r = 0.05 , include_checks = TRUE ) {
25+ private $ risk_free_rate <- r
26+ private $ validate_params <- include_checks
27+ invisible (self )
28+ },
29+
30+ # ' @description Calculate call option price
31+ # ' @param S Current stock price
32+ # ' @param K Strike price
33+ # ' @param T Time to expiration (in years)
34+ # ' @param sigma Volatility (annualized)
35+ calculate_call_price = function (S , K , T , sigma ) {
36+ if (private $ validate_params ) {
37+ private $ validate_inputs(S , K , T , sigma )
38+ }
39+
40+ d1 <- private $ calculate_d1(S , K , T , sigma )
41+ d2 <- private $ calculate_d2(d1 , sigma , T )
42+
43+ call_price <- S * stats :: pnorm(d1 ) - K * exp(- private $ risk_free_rate * T ) * stats :: pnorm(d2 )
44+ return (call_price )
45+ },
46+
47+ # ' @description Calculate put option price
48+ # ' @param S Current stock price
49+ # ' @param K Strike price
50+ # ' @param T Time to expiration (in years)
51+ # ' @param sigma Volatility (annualized)
52+ calculate_put_price = function (S , K , T , sigma ) {
53+ if (private $ validate_params ) {
54+ private $ validate_inputs(S , K , T , sigma )
55+ }
56+
57+ d1 <- private $ calculate_d1(S , K , T , sigma )
58+ d2 <- private $ calculate_d2(d1 , sigma , T )
59+
60+ put_price <- K * exp(- private $ risk_free_rate * T ) * stats :: pnorm(- d2 ) - S * stats :: pnorm(- d1 )
61+ return (put_price )
62+ },
63+
64+ # ' @description Calculate all Greeks for a call option
65+ # ' @param S Current stock price
66+ # ' @param K Strike price
67+ # ' @param T Time to expiration (in years)
68+ # ' @param sigma Volatility (annualized)
69+ calculate_call_greeks = function (S , K , T , sigma ) {
70+ if (private $ validate_params ) {
71+ private $ validate_inputs(S , K , T , sigma )
72+ }
73+
74+ d1 <- private $ calculate_d1(S , K , T , sigma )
75+ d2 <- private $ calculate_d2(d1 , sigma , T )
76+
77+ # Calculate Greeks
78+ delta <- stats :: pnorm(d1 )
79+ gamma <- stats :: dnorm(d1 ) / (S * sigma * sqrt(T ))
80+ theta <- (- S * stats :: dnorm(d1 ) * sigma / (2 * sqrt(T )) -
81+ private $ risk_free_rate * K * exp(- private $ risk_free_rate * T ) * stats :: pnorm(d2 ))
82+ vega <- S * sqrt(T ) * stats :: dnorm(d1 )
83+ rho <- K * T * exp(- private $ risk_free_rate * T ) * stats :: pnorm(d2 )
84+
85+ return (list (
86+ delta = delta ,
87+ gamma = gamma ,
88+ theta = theta ,
89+ vega = vega ,
90+ rho = rho
91+ ))
92+ },
93+
94+ # ' @description Calculate all Greeks for a put option
95+ # ' @param S Current stock price
96+ # ' @param K Strike price
97+ # ' @param T Time to expiration (in years)
98+ # ' @param sigma Volatility (annualized)
99+ calculate_put_greeks = function (S , K , T , sigma ) {
100+ if (private $ validate_params ) {
101+ private $ validate_inputs(S , K , T , sigma )
102+ }
103+
104+ d1 <- private $ calculate_d1(S , K , T , sigma )
105+ d2 <- private $ calculate_d2(d1 , sigma , T )
106+
107+ # Calculate Greeks
108+ delta <- stats :: pnorm(d1 ) - 1
109+ gamma <- stats :: dnorm(d1 ) / (S * sigma * sqrt(T ))
110+ theta <- (- S * stats :: dnorm(d1 ) * sigma / (2 * sqrt(T )) +
111+ private $ risk_free_rate * K * exp(- private $ risk_free_rate * T ) * stats :: pnorm(- d2 ))
112+ vega <- S * sqrt(T ) * stats :: dnorm(d1 )
113+ rho <- - K * T * exp(- private $ risk_free_rate * T ) * stats :: pnorm(- d2 )
114+
115+ return (list (
116+ delta = delta ,
117+ gamma = gamma ,
118+ theta = theta ,
119+ vega = vega ,
120+ rho = rho
121+ ))
122+ },
123+
124+ # ' @description Estimate implied volatility using Newton-Raphson method
125+ # ' @param market_price Observed market price of the option
126+ # ' @param S Current stock price
127+ # ' @param K Strike price
128+ # ' @param T Time to expiration (in years)
129+ # ' @param is_call Whether the option is a call (TRUE) or put (FALSE)
130+ # ' @param tolerance Convergence tolerance
131+ # ' @param max_iter Maximum iterations
132+ estimate_implied_volatility = function (market_price , S , K , T ,
133+ is_call = TRUE , tolerance = 1e-5 , max_iter = 100 ) {
134+ if (private $ validate_params ) {
135+ if (market_price < = 0 ) stop(" Market price must be positive" )
136+ private $ validate_inputs(S , K , T , 0.5 ) # Initial volatility check
137+ }
138+
139+ # Initial guess for volatility
140+ sigma <- sqrt(2 * abs(log(S / K ) + private $ risk_free_rate * T ) / T )
141+ sigma <- min(max(0.1 , sigma ), 5 ) # Bound initial guess
142+
143+ for (i in 1 : max_iter ) {
144+ # Calculate price and vega
145+ if (is_call ) {
146+ price <- self $ calculate_call_price(S , K , T , sigma )
147+ greeks <- self $ calculate_call_greeks(S , K , T , sigma )
148+ } else {
149+ price <- self $ calculate_put_price(S , K , T , sigma )
150+ greeks <- self $ calculate_put_greeks(S , K , T , sigma )
151+ }
152+
153+ diff <- price - market_price
154+
155+ if (abs(diff ) < tolerance ) {
156+ return (sigma )
157+ }
158+
159+ # Update volatility estimate using Newton-Raphson
160+ sigma <- sigma - diff / greeks $ vega
161+
162+ # Bound the volatility
163+ sigma <- min(max(0.001 , sigma ), 5 )
164+ }
165+
166+ warning(" Implied volatility did not converge" )
167+ return (sigma )
168+ }
169+ ),
170+
171+ private = list (
172+ risk_free_rate = NULL ,
173+ validate_params = NULL ,
174+
175+ calculate_d1 = function (S , K , T , sigma ) {
176+ (log(S / K ) + (private $ risk_free_rate + sigma ^ 2 / 2 ) * T ) / (sigma * sqrt(T ))
177+ },
178+
179+ calculate_d2 = function (d1 , sigma , T ) {
180+ d1 - sigma * sqrt(T )
181+ },
182+
183+ validate_inputs = function (S , K , T , sigma ) {
184+ if (S < = 0 ) stop(" Stock price must be positive" )
185+ if (K < = 0 ) stop(" Strike price must be positive" )
186+ if (T < = 0 ) stop(" Time to expiration must be positive" )
187+ if (sigma < = 0 ) stop(" Volatility must be positive" )
188+ }
189+ )
190+ )
191+
192+ # Demonstration
193+ demonstrate_black_scholes <- function () {
194+ cat(" === Black-Scholes Option Pricing Demo ===\n\n " )
195+
196+ # Initialize calculator
197+ bs <- BlackScholesCalculator $ new(r = 0.05 )
198+
199+ # Example parameters
200+ S <- 100 # Current stock price
201+ K <- 100 # Strike price
202+ T <- 1 # One year to expiration
203+ sigma <- 0.2 # 20% volatility
204+
205+ # Calculate option prices
206+ call_price <- bs $ calculate_call_price(S , K , T , sigma )
207+ put_price <- bs $ calculate_put_price(S , K , T , sigma )
208+
209+ cat(" Parameters:\n " )
210+ cat(sprintf(" Stock Price: $%.2f\n " , S ))
211+ cat(sprintf(" Strike Price: $%.2f\n " , K ))
212+ cat(sprintf(" Time to Expiration: %.1f years\n " , T ))
213+ cat(sprintf(" Volatility: %.1f%%\n " , sigma * 100 ))
214+ cat(sprintf(" Risk-free Rate: %.1f%%\n\n " , bs $ risk_free_rate * 100 ))
215+
216+ cat(" Option Prices:\n " )
217+ cat(sprintf(" Call Option: $%.2f\n " , call_price ))
218+ cat(sprintf(" Put Option: $%.2f\n\n " , put_price ))
219+
220+ # Calculate and display Greeks
221+ call_greeks <- bs $ calculate_call_greeks(S , K , T , sigma )
222+ put_greeks <- bs $ calculate_put_greeks(S , K , T , sigma )
223+
224+ cat(" Call Option Greeks:\n " )
225+ cat(sprintf(" Delta: %.4f\n " , call_greeks $ delta ))
226+ cat(sprintf(" Gamma: %.4f\n " , call_greeks $ gamma ))
227+ cat(sprintf(" Theta: %.4f\n " , call_greeks $ theta ))
228+ cat(sprintf(" Vega: %.4f\n " , call_greeks $ vega ))
229+ cat(sprintf(" Rho: %.4f\n\n " , call_greeks $ rho ))
230+
231+ cat(" Put Option Greeks:\n " )
232+ cat(sprintf(" Delta: %.4f\n " , put_greeks $ delta ))
233+ cat(sprintf(" Gamma: %.4f\n " , put_greeks $ gamma ))
234+ cat(sprintf(" Theta: %.4f\n " , put_greeks $ theta ))
235+ cat(sprintf(" Vega: %.4f\n " , put_greeks $ vega ))
236+ cat(sprintf(" Rho: %.4f\n\n " , put_greeks $ rho ))
237+
238+ # Demonstrate implied volatility calculation
239+ test_market_price <- call_price * 1.1 # Use 10% higher price for demonstration
240+ implied_vol <- bs $ estimate_implied_volatility(test_market_price , S , K , T , is_call = TRUE )
241+ cat(" Implied Volatility Estimation:\n " )
242+ cat(sprintf(" Market Price: $%.2f\n " , test_market_price ))
243+ cat(sprintf(" Implied Volatility: %.1f%%\n " , implied_vol * 100 ))
244+
245+ cat(" \n === Demo Complete ===\n " )
246+ }
247+
248+ # Run demonstration if not in interactive mode
249+ if (! interactive()) {
250+ demonstrate_black_scholes()
251+ }
0 commit comments