qfinpy 0.0.3

A powerful, easy-to-use library for Quantitative Finance

Introduction

QFinPy is a powerful, easy-to-use Python library designed for quantitative finance research, analysis, and modeling. It provides a set of tools for creating options payoff diagrams, pricing derivatives, Monte Carlo Simulations, time series analysis, and constructing portfolios.

Installation

pip install qfinpy

Usage and Examples

import qfinpy as qf
import numpy as np
import matplotlib.pyplot as plt

Interest Rate

rate_to_continuous(r, t=1, compounded=1)
rate_to_compounded(r, t)
future_value(M, r, T, a=0.0)
present_value(M, r, T)

A nominal interest rate of 8% per year compounded quarterly --> continuously compounded rate --> effective annual rate (EAR)

r = qf.rate_to_continuous(0.08, t=1, compounded=1/4)
print('continuous rate = ', r)
ear = qf.math.exp(r) - 1
print('effective annual rate = ', ear)

continuous rate =  0.07921050918471892
effective annual rate =  0.08243215999999998

Continuously compounded rate --> quarterly rate --> nominal yearly rate

r_quarterly = qf.rate_to_compounded(r, 1/4)
print('3 months rate = ', r_quarterly)
# Nominal annual rate
r_quarterly_nominal = 4 * r_quarterly
print('Nominal annual rate = ', r_quarterly_nominal)

3 months rate =  0.020000000000000018
Nominal annual rate =  0.08000000000000007

Present and Future Values: 1000 units invested for a period of 10 years at a continuously compounded rate of 7%. (future_value function also takes an optional argument 'a', which is the continuous deposit rate)

fv = qf.future_value(1000, 0.07, 10)
print('Future value = ', fv)

Future value =  2013.7527074704767

pv = qf.present_value(fv, 0.07, 10)
print('Present value = ', pv)

Present value =  1000.0

Options

Options Payoff Diagrams

C (European Calls), P (European Puts), BC (Binary Calls), and BP (Binary Puts) classes needs to be instantiated using strike price. The S (stock) class doesn't require strike price. The function options_payoff_diagram takes an optional argument 'u_price' for the current underlying price. u_price is required when S (stock) is included.

European Call with strike 100

qf.options_payoff_diagram(qf.C(100))

European Call with strike 60 and European Put with strike 40, and current underlying price at 50.

folio = qf.C(60) + qf.P(40)
qf.options_payoff_diagram(folio, 50)

4 P(50) short, 4 P(70) long, 6 C(90) long, 2 C(110) short, 4 C(120) short and 1 S() short.

folio = -4*qf.P(50) + 4*qf.P(70) + 6*qf.C(90) - 2*qf.C(110) - 4*qf.C(120) - qf.S()
qf.options_payoff_diagram(folio, u_price=80)

Black Scholes Oprion Pricing and Greeks

black_scholes_value(type, strike, t, u_price, vol, rf_rate, u_yield=0)
  type: 'C' (European Calls), 'P' (European Puts), 'BC' (Binary Calls), and 'BP' (Binary Puts)
  strike: strike/exersise price
  t: time to exersise
  u_price: underlying asset price
  vol: underlying asset volatility
  rf_rate: riskfree rate of interest
  u_yield: continuous dividend yield on the asset

black_scholes_delta(type, strike, t, u_price, vol, rf_rate, u_yield=0)
black_scholes_gamma(type, strike, t, u_price, vol, rf_rate, u_yield=0)
black_scholes_theta(type, strike, t, u_price, vol, rf_rate, u_yield=0)
black_scholes_speed(type, strike, t, u_price, vol, rf_rate, u_yield=0)
black_scholes_vega(type, strike, t, u_price, vol, rf_rate, u_yield=0)
black_scholes_rho(type, strike, t, u_price, vol, rf_rate, u_yield=0)
black_scholes_yield_sensitivity(type, strike, t, u_price, vol, rf_rate, u_yield=0)

Example: European Call option with
strike = 95
t = 3/12
u_price = 100
vol = 0.50
rf_rate = 0.01
u_yield=0

value = qf.black_scholes_value('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
delta = qf.black_scholes_delta('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
gamma = qf.black_scholes_gamma('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
theta = qf.black_scholes_theta('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
speed = qf.black_scholes_speed('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
vega = qf.black_scholes_vega('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
rho = qf.black_scholes_rho('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
yield_sensitivity = qf.black_scholes_yield_sensitivity('C', 95, 3/12, 100, 0.50, 0.01, u_yield=0)
print('option value = ', value)
print('delta = ', delta)
print('gamma = ', gamma)
print('theta = ', theta)
print('speed = ', speed)
print('vega = ', vega)
print('rho = ', rho)
print('yield_sensitivity = ', yield_sensitivity)

option value =  12.527923392521458
delta =  0.633136941899257
gamma =  0.015060599447748629
theta =  -19.33360701765983
speed =  -0.000355534473275545
vega =  18.825749309685786
rho =  12.69644269935106
yield_sensitivity =  -15.828423547481425

Implied Volatility

implied_volatility(type, deriv_price, strike, t, u_price, rf_rate, u_yield=0, x0=0.1, tol=1.48e-08)
  type: 'C' (European Calls), 'P' (European Puts), 'BC' (Binary Calls), and 'BP' (Binary Puts)
  deriv_price: Current market price of the option
  strike: strike/exersise price
  t: time to exersise
  u_price: underlying asset price
  rf_rate: risk-free rate of interest
  u_yield: continuous dividend yield on the asset
  x0: initial guess of implied volatility in the Newton's method
  tol: tol in the Newton's method

Example: A European call option, with a strike price of 50 expires in 32 days. The risk-free interest rate is 0.05 and the stock is currently trading at 51.25 and the current market price of the option is $2.00. Using a standard Black–Scholes pricing model, the volatility implied by the market price is 0.187

imp_vol = qf.implied_volatility('C', 2.00, 50, 32/365, 51.25, 0.05)
print('Implied volatility = ', imp_vol)

Implied volatility =  0.1869228434755648

Portfolio Optimization

portfolio_optim(m, c, expected_return=None, shortable=None, rf_rate=None, allow_borrow=False, max_leverage=1.0e3)
  m: Returns of assets
  c: covariance of assets
  expected_return: expected returns for the portfolio
  shortable: list of 0's and 1's were 1's represent the assets that can be shorted
  rf_rate: risk-free rate
  allow_borrow: Allow risk-free asset to be borrowed
  max_leverage: Maximum leverage

m = np.array([0.0890833, 0.213667, 0.234583])
c = np.array([[0.01080754, 0.01240721, 0.01307513],
     [0.01240721, 0.05839170, 0.05542639],
     [0.01307513, 0.05542639, 0.09422681]])

result = qf.portfolio_optim(m, c, expected_return=0.15)
result.x

array([0.53009593, 0.35639213, 0.113512  ])

Monte Carlo Simulation

Random series generation

Normal random numbers

normal(*n, mu=np.array([0.0]), sigma=np.array([1.0]), bs=None, dtype=np.float64)
  n: Number of samples
  mu: mean - array of shape () or (n,) or bs+(1,) or bs+(n,)
  sigma: standard deviation - array of shape () or (1) or (n,) or bs+(1,) or bs+(n,)
  bs: batch size/shape - int or tuple of ints
  dtype: dtype

x = qf.normal(10)
print(x)

[-0.29769324 -0.54351292 -0.34507348  0.48251412  0.22500892  0.87064283
 -0.50450492  0.30324753  1.15314968 -0.67326597]

plt.plot(qf.normal(100, mu=10.0, sigma=3.0))

[<matplotlib.lines.Line2D at 0x7f10ff9951c0>]

plt.plot(qf.normal(mu=np.linspace(0,10,100), sigma=1.0))

[<matplotlib.lines.Line2D at 0x7f1155fdf3b0>]

plt.plot(qf.normal(mu=0, sigma=np.linspace(0,10,100)))

[<matplotlib.lines.Line2D at 0x7f1100347650>]

plt.plot(qf.normal(mu=np.linspace(0,10,100), sigma=np.linspace(1,5,100)))

[<matplotlib.lines.Line2D at 0x7f11003605f0>]

x = qf.normal(10, bs=2)
print(x.shape)
print(x)

(2, 10)
[[ 0.84284109  0.08477504 -2.43723961  0.61270538  0.9801093  -0.52259732
   1.38707595 -0.73983371 -0.66849634  1.46785926]
 [ 0.13333635 -0.90268586  0.56882097  1.1391172   0.2177926   0.41764538
   0.1300695  -0.21647265 -0.32227474 -1.62438587]]

x = qf.normal(20, bs=(2,3))
print(x.shape)

(2, 3, 20)

Multivariate normal random variables

normal_multivariate(*n, mu=np.zeros(2), cov=np.eye(2), bs=None, dtype=np.float64)
  n: Number of samples
  mu: mean - array of shape (k) or (k,n) or bs+(k,1) or bs+(k,n)
  cov: covariance - array of shape (k,k) or (k,k,n) or bs+(k,k,1) or bs+(k,k,n)
  bs: batch size/shape - int or tuple of ints
  dtype: dtype

x = qf.normal_multivariate(10, mu=np.ones(2), cov=np.eye(2))
print(x.shape)
print(x)

(2, 10)
[[ 0.77333346  0.17601791 -0.01866713  0.46400412  1.51313251 -0.90469188
   0.13323023  1.54159995  0.21370799  1.8932842 ]
 [-0.35885229  0.88045584 -0.19873425  1.01257887  0.45947712 -0.20288216
   0.84080326  2.95304337  0.45257102  1.14213169]]

x = qf.normal_multivariate(mu=np.stack((np.linspace(0,10,100), np.linspace(0,5,100))), cov=[[1,0.5],[0.5,2]])
plt.plot(x[0])
plt.plot(x[1])
plt.show()

x = qf.normal_multivariate(10, mu=np.ones(2), cov=np.eye(2), bs=(5,7))
print(x.shape)

(5, 7, 2, 10)

Lognormal (=exp(normal))

x = qf.lognormal(10)
print(x)

[0.50643209 0.23344932 4.29740853 0.77359027 0.72411853 0.21328177
 3.03007789 0.69455528 3.15378758 0.43217616]

x = qf.lognormal_multivariate(10, mu=np.ones(2), cov=np.eye(2), bs=(5,7))
print(x.shape)

(5, 7, 2, 10)

Student's t dist

x = qf.students_t(100, df=4, mu=2, sigma=4)
plt.plot(x)

[<matplotlib.lines.Line2D at 0x7669ed65be30>]

Random Walk

Additive Random Walk

random_walk(series, x0=0.0)
  series: random step size
  x0: starting point

ret = qf.normal(100)
x = qf.random_walk(ret, x0=0.0)
plt.plot(x)

[<matplotlib.lines.Line2D at 0x7f11551a0ef0>]

Geometric Random Walk

random_walk_geometric(series, x0=1.0)   series: random returns
  x0: starting point

ret = qf.normal(100, mu=0.005, sigma=0.05)
x = qf.random_walk_geometric(ret, x0=1.0)
plt.plot(x)

[<matplotlib.lines.Line2D at 0x7f115516fec0>]

Monte Carlo Option Pricing

strike = 95
t = 3/12
u_price = 100
vol = 0.50
rf_rate = 0.01

n=1000
bs=10000
mu = rf_rate*t/n
sigma = vol*(t/n)**0.5
ret = qf.normal(n, mu=mu, sigma=sigma, bs=bs)
sim = qf.random_walk_geometric(ret, x0=u_price)
print(sim.shape)

(10000, 1001)

plt.plot(sim[0])

[<matplotlib.lines.Line2D at 0x7f11550620c0>]

value = qf.present_value(np.maximum(0, sim[:,-1] - strike).mean(), rf_rate, t)
print(value)

12.453812729122065

Since intermediate values of geometric random walk are not required for option value calculation, we can use lognormal distribution to save time and memory. In this case standard deviations of both log returns and normal returns are same, but (the mean of log returns) = (mean of normal returns) - (variance of returns)/2

n=1
bs=1000000
sigma = vol*(t**0.5)
mu = rf_rate*t - (sigma**2)/2
sim = u_price * qf.lognormal(n, mu=mu, sigma=sigma, bs=bs)
print(sim.shape)

(1000000, 1)

value = qf.present_value(np.maximum(0, sim - strike).mean(), rf_rate, t)
print('option value = ', value)

option value =  12.527805286619587

Time Series Analysis

import qfinpy as qf
import numpy as np
import matplotlib.pyplot as plt

from qfinpy import tsa

Moving Average MA(q)

ma(series, theta, mu=0.0, e0=None)

e = qf.normal(150, bs=5)
theta = [0.6, 0.2, 0.1]
w = tsa.ma(e, theta, mu=0.0)
print(e.shape)
print(w.shape)

(5, 150)
(5, 150)

plt.plot(e[0])

[<matplotlib.lines.Line2D at 0x763c1f1878c0>]

plt.plot(w[0])

[<matplotlib.lines.Line2D at 0x763c1f1b7140>]

e_m = qf.normal_multivariate(150, mu=np.zeros(2), cov=np.eye(2))
theta_m = [[[0.3,0.1],[0.1,0.3]], [[0.2,0.1],[0.1,0.2]], [[0.1,0.1],[0.1,0.1]]]
w_m = tsa.ma(e_m, theta_m)
print(e_m.shape)
print(w_m.shape)

(2, 150)
(2, 150)

Autoregrassive AR(p)

ar(series, phi, mu=0.0, x0=None)

phi = [0.4, 0.2]
x = tsa.ar(w, phi, mu=0.0)
print(x.shape)

(5, 150)

plt.plot(x[0])

[<matplotlib.lines.Line2D at 0x763c1dd10bf0>]

phi_m = [[[0.4, 0.1],[0.1,0.4]], [[0.2, 0.1],[0.1,0.2]]]
x_m = tsa.ar(w_m, phi_m)
print(x_m.shape)

(2, 150)

GARCH(r,s)

gh(series, w, alpha, beta, x0=None, e0=None, mu=0, initial_var=None):
(returns a tuple containing the series and the volatility)

g = tsa.gh(x, 0.01, 0.3, 0.6, mu=0)
g[0].shape

(5, 150)

plt.plot(g[0][0])

[<matplotlib.lines.Line2D at 0x763c1542e630>]

Inverse GARCH, MA and AR

ma_inverse(series, theta, mu=0.0)
ar_inverse(series, phi, mu=0.0)
gh_inverse(series, w, alpha, beta, mu=0, initial_var=None)

inverse_g = tsa.gh_inverse(g[0], 0.01, 0.3, 0.6, mu=0)
print(inverse_g[0].shape)

(5, 149)

plt.plot(x[0,1:])
plt.plot(inverse_g[0][0])

[<matplotlib.lines.Line2D at 0x763c153880e0>]

inverse_x = tsa.ar_inverse(x, phi)
print(inverse_x.shape)

(5, 148)

plt.plot(w[0][2:])
plt.plot(inverse_x[0])
plt.show()

inverse_w = tsa.ma_inverse(w, theta)
print(inverse_w.shape)

(5, 147)

plt.plot(e[0][-147:])
plt.plot(inverse_w[0])
plt.show()

Example: ARMA(1,1) + GARCH(1,1) fit using the inverse functions

# simulated data
e = tsa.gh(qf.normal(2000), w=0.1, alpha=0.3, beta=0.6)[0]
w = tsa.ma(e, theta=0.8)
data = tsa.ar(w, phi=0.5, mu=2.0)

plt.plot(data)

[<matplotlib.lines.Line2D at 0x7669f85c7b90>]

from scipy.optimize import minimize

def arma_1_1_garch_1_1_log_likelihood(params, data):
    # Parameters
    mu, phi, theta, omega, alpha, beta = params

    N = len(data)
    w = tsa.ar_inverse(data-mu, phi=phi)
    e = tsa.ma_inverse(w, theta=theta)
    z, sigma2 = tsa.gh_inverse(e, w=omega, alpha=alpha, beta=beta)

    log_likelihood = -0.5 * N * np.log(2 * np.pi)
    log_likelihood -= 0.5 * np.sum(np.log(sigma2) + z**2)
    return -log_likelihood  # Negative log-likelihood for minimization

# Initial parameter guess 
initial_params = [0.0, 0.1, 0.1, 0.1, 0.1, 0.1]

bounds = [
    (None, None),  # mu
    (-1, 1),       # phi1
    (None, None),  # theta12
    (1e-8, None),  # omega > 0
    (1e-8, 1),     # 0 <= alpha <= 1
    (1e-8, 1)      # 0 <= beta <= 1
]

y = data
# Optimize log-likelihood
result = minimize(arma_1_1_garch_1_1_log_likelihood, initial_params, args=(y,), bounds=bounds)
fitted_params = result.x

print("Fitted Parameters:", fitted_params)
print("fitted_log_likelihood = ", - result.fun)
print("AIC = ", -2 * -result.fun + 2 * len(result.x))
print("BIC = ", -2 * -result.fun + len(result.x) * np.log(y.size).item())

Fitted Parameters: [2.00207328 0.51945168 0.79517484 0.08202248 0.33269634 0.60103715]
fitted_log_likelihood =  -2575.529401410069
AIC =  5163.058802820138
BIC =  5196.66421757739

Autovariance and Autocorrelation

acovf(series, nlags)
acf(series, nlags, qstat=False, dof_offset=0)

tsa.acf(data, nlags=5, qstat=True)

(array([1.        , 0.74846902, 0.36215007, 0.15446444, 0.06303216,
        0.03298029]),
 array([1122.09319475, 1384.92367546, 1432.76167471, 1440.73166649,
        1442.91469901]),
 array([0., 0., 0., 0., 0.]))

data_acf = tsa.acf(data, nlags=15)
plt.stem(data_acf, markerfmt='')

<StemContainer object of 3 artists>

Other Tools

sliding_window(a, window)
normal_cdf(x, mu=0.0, sigma=1.0)
normal_pdf(x, mu=0.0, sigma=1.0)
chi2_cdf(x, df)
chi2_pdf(x, df)
trimmed_mean(series, trim=None, axis=-1)
mad(series, calib=1.4826, axis=-1): Mean absolute deviation
pct_change(series, axis=-1)
log_change(series, axis=-1) \

x = qf.normal(10)
x

array([-1.81506319, -0.06177757,  0.81406117, -0.13047454, -0.08597478,
       -0.18929597,  0.84590606,  0.92444569,  0.88397027, -0.39034777])

qf.sliding_window(x, 4)

array([[-1.81506319, -0.06177757,  0.81406117, -0.13047454],
       [-0.06177757,  0.81406117, -0.13047454, -0.08597478],
       [ 0.81406117, -0.13047454, -0.08597478, -0.18929597],
       [-0.13047454, -0.08597478, -0.18929597,  0.84590606],
       [-0.08597478, -0.18929597,  0.84590606,  0.92444569],
       [-0.18929597,  0.84590606,  0.92444569,  0.88397027],
       [ 0.84590606,  0.92444569,  0.88397027, -0.39034777]])

qf.normal_cdf(0, mu=0.0, sigma=1.0)

0.5

qf.normal_pdf(0, mu=0.0, sigma=1.0)

0.3989422804014327

qf.chi2_cdf(5, 4)

np.float64(0.7127025048163542)