MC Simulation – Python Code

This is a simple Python code for pricing European call and put options using Monte Carlo simulation. According to the law of large numbers, as the number of simulations increase , the actual ratio will converge to the expected ratio of outcomes. In this case, for a large number of simulations the concluded prices will converge to the values generated by the Black-Scholes closed-form solution.

## European call / put option

import numpy as np
#import pandas as pd

# inputs
S0 = 10 # stock price
K = 10 # strike price
sigma = 0.2 # volatility
t = 0.5 # number of years
r = 0.02 # risk-free rate
n = 1000 # number of simulations
Dt = 1/252 # frequency; 1/252 is daily
h = int(t / Dt) # number of simulated periods

# simulate stock prices using MC
rnd = np.array([[0.0 for x in range(h + 1)] for y in range(n)])
for j in range(0, h + 1, 1):
rnd[:, j] = np.random.normal(0, 1, n)

S = np.array([[0.0 for x in range(h + 1)] for y in range(n)])
S[:, 0] = S0
for i in range(0, n, 1):
for j in range(1, h + 1, 1):
S[i, j] = S[i, j - 1] * np.exp((r - 0.5 * sigma ** 2) * Dt + sigma * rnd[i, j - 1] * np.sqrt(Dt))

# calculate call price
c = np.mean(np.maximum(0, S[:, h - 1] - K))

# calculate put price
p = np.mean(np.maximum(0, K - S[:, h - 1]))

Black-Scholes – Python Code

This is a simple Python code for calculating the price of a European call or put option using the Black-Scholes-Merton closed-form solution. The inputs can be easily modified in the corresponding block (initial stock price, strike price, annual volatility, time, risk-free rate and dividend yield).

## European call / put option

import numpy as np
from scipy.stats import norm

# inputs
S = 10 # stock price
K = 10 # strike price
sigma = 0.2 # volatility
t = 1.5 # number of years
r = 0.02 # risk-free rate
d = 0 # expected dividend yield

# calculate d1 & d2
d1 = (np.log(S/K) + (r - d + 0.5*sigma ** 2)*t) / (sigma * np.sqrt(t))
d2 = d1 - sigma * np.sqrt(t)

# calculate call price
c = S * norm.cdf(d1) * np.exp(-d * t) - K * norm.cdf(d2) * np.exp(-r * t)

# calculate put price
p = K * norm.cdf(-d2) * np.exp(-r * t) - S * norm.cdf(-d1) * np.exp(-d * t)

Correlated Random Numbers – R Code

This is simple R code for generating correlated random numbers using the package MASS and its built-in function mvrnorm. In the example below 3 correlated vectors of values are being generated with correlations of 0.5, 0.6 and 0.7, respectively.

## Generate Correlated Random Numbers

# clear environment and console
rm(list = ls())
cat("\014")

# install and load required packages
install.packages("MASS")
require(MASS)

# inputs
n = 1000 # number of random numbers
correl = matrix(c(1, 0.5, 0.6, 0.5, 1, 0.7, 0.6, 0.7, 1), ncol = 3)
mu = c(0, 0, 0) # Gaussian vectors of values

num = mvrnorm(n, mu, correl, empirical = TRUE)

Multivariable Optimization – R code

This is simple R code on optimizing a function with two variables and a lower and upper bound for each one of them.

## Optimization - Minimum of function

# clear environment and console
rm(list = ls())
cat("\014")

# inputs
x = rnorm(1:100)
y = rnorm(1:100) * 100

fun = function(z, x, y)
{
sum(x * z[1] - (y ^ 2) * z[2])
}

# estimate optimal values for z[1] and z[2]
opt = optim(par = c(1, 1), fun, x = x, y = y, lower = c(-1, -1), upper = c(10, 10), method="L-BFGS-B")
opt$par
opt$value

Bayesian Regression Analysis – R Code

This is an R code on regression analysis: a Bayesian regression model is being implemented. The time series are randomly created so no meaningful results can be expected. Moreover, non-informative priors are utilized. Using proper priors the estimates would differ (the estimated slope is the mean of the OLS estimate and the information derived by the priors).

## Bayesian linear regression

# clear environment and console
rm(list = ls())
cat("\014")

#install and load required packages
install.packages("BAS")
library("BAS")

# create random series for the independent variable (var1) and dependent
# variables (var2 var3)
v1 = data.frame(var1 = runif(100, 20, 120))
v2 = data.frame(var2 = runif(100, 20, 120))
v3 = data.frame(var3 = runif(100, 20, 120))

# concatenate the time series
vv = c(v1, v2, v3)

# run a Bayesian regression
reg1 = bas.lm(var1 ~ ., data = vv,
prior = 'BIC', modelprior = uniform())
reg1

# run a linear regression
reg2 = lm(var1 ~ ., data = vv)
summary(reg2)

Principal Component Analysis – R Code

This is an implementation of the principal component analysis (“PCA”) method on stock prices log returns in R.

# Principal Component Analysis

# clear environment and console
rm(list = ls())
cat("\014")

# install and load required packages
# install.packages("stats")
library("stats")

data(EuStockMarkets)

# calculate log returns
log_prices <- log(EuStockMarkets[, 1:4])
ret = data.frame(diff(as.matrix(log_prices,lag=1,difference=1)))

# implement PCA on log returns of index levels
pca <- princomp(ret)

# plot method & print summary
plot(prcomp(ret))
biplot(prcomp(ret, scale = TRUE))
summary(prcomp(ret, scale = TRUE))

# view the scores for each date and the loadings for the components
pca$scores
pca$loadings

Pricing Structured Products – Matlab Code

In this Matlab code we implement a prepayment modeling method and eventually price Asset Backed Securities (ABS).


function price = MBS(k, theta)

%1a
P = 100000;
WACC = 0.08;
r = WACC / 12;
n = 360;
SY = [0.94 0.76 0.74 0.95 0.98 0.92 0.98 1.10 1.18 1.22 1.23 0.98];
SY = repmat(SY, 1, n);
WACC = repmat(WACC, n, 1)';
r10 = cir_rates(360, 120, 0.078, k, theta, 0.09);

RI = zeros(n + 1, 1); BU = zeros(n + 1, 1); SG = zeros(n + 1, 1);
CPR = zeros(n + 1, 1); w = zeros(n, 1); PV = zeros(n, 1);
c = zeros(n, 1); TPP = zeros(n + 1, 1);
PV(1) = P;

for i = 1 : n
w(i) = - 8.57 + 430 * (WACC(i) - r10(i));
RI(i + 1) = 0.28 + 0.14 * atan(w(i));
BU(i + 1) = 0.3 + 0.7 * PV(i) / PV(1);
SG(i + 1) = min(1, (i + 1) / 30);
CPR(i + 1)= RI(i + 1) * BU(i + 1) * SG(i + 1) * SY(i + 1);
TPP(i + 1) = PV(i) * r * (1 / (1 - (1 + r) ^ (- n + (i - 1))) - 1) + ...
(PV(i) - PV(i) * r * (1 / (1 - (1 + r) ^ (- n + (i - 1))) - 1)) * ...
(1 - (1 - CPR(i + 1)) ^ (1 / 12));
PV(i + 1) = PV(i) - TPP(i + 1);
c(i + 1) = PV(i) * r / (1 - (1 + r) ^ (- n + (i - 1))) + (PV(i) - ...
PV(i) * r * (1 / (1 - (1 + r) ^ (- n + (i - 1))) - 1)) * ...
(1 - (1 - CPR(i + 1)) ^ (1 / 12));
end

l = c(2 : 361).* cir_exp(360, 12, 0.078, k, theta, 0.09)';
price = mean(sum(l));

end

function r2 = cir_rates(N, S, r0, k, theta, sigma)

dt = 1 / S;
rep = r0 * ones(N, S);
z = sqrt(dt) * randn(N, S);

for i = 1 : S
for j = 1 : N
rep(j, i + 1) = rep(j, i) + dt * k * (theta - rep(j, i)) + sigma * sqrt(rep(j, i))
* z(j, i) ;
end
end

w = zeros(N, 1);
for j = 1 : N
w(j) = mean(rep(j, :));
end

r2 = w;
end

function r2 = cir_exp(N, S, r0, k, theta, sigma)

dt = 1 / S;
rep = r0 * ones(N, S);
z = sqrt(dt) * randn(N, S);

for i = 1 : S
for j = 1 : N
rep(j, i + 1) = rep(j, i) + dt * k * (theta - rep(j, i)) + sigma * ...
sqrt(rep(j, i)) * z(j, i) ;
end
end

Rf = zeros(1, N);
for i = 1 : N
Rf(i) = mean(rep(i, :)) * i * dt;
end

r2 = exp(- Rf);

end

CIR Model – Matlab

This is the implementation of the CIR model in Matlab (interest rate modeling).

function r2 = cir3(M, S)
% dr_t = k (rt - r_t) dt + sigma dW_t

M = 10000; S = 180;
N = 365;
t = N / 365;
FV = 1000;
K = 950;
k = 0.92;
dt = 1 / S;
theta = 0.055;
sigma = 0.12 / 365;
r0 = 0.05;
T = S / 365 ;

rep = r0 * ones(M + 1, S+1);
z = sqrt(dt) * randn(M + 1, S+1);

for i = 1 : S - N
for j = 1 : M + 1
rep(j, i + 1) = rep(j, i) + dt * k * (theta - rep(j, i)) + sigma * sqrt(rep(j, i)) * z(j, i) ;
end
end

Rf = zeros(1, M + 1);
for i = 1 : M + 1
Rf(i) = dt * sum(rep(i, :));
Rf(i) = exp(-Rf(i) * T);
end

r = mean(rep);
V = FV * mean(Rf);
r2 = mean(max(0, exp(r * t) * (V - K)));
end




Vasicek Model – Matlab Code

This is an implementation of the Vasicek model for interest rate modeling in Matlab.

function r1 = Vasicek3(M, N, r0)
% dr_t = k (rt - r_t) dt + sigma dW_t

k = 0.82;
dt = 1 / N;
theta = 0.05;
sigma = 0.10 / 365;
T = N / 365 ;

rep = r0 * ones(M + 1, N + 1);
z = sqrt(dt) * randn(M + 1, N + 1);

for i = 1 : N
rep(:, i + 1) = rep(:, i) + dt * k * (theta - rep(:, i)) + sigma * z(:, i) ;
end

Rf = zeros(1, M + 1);

for i = 1 : M + 1
Rf(i) = dt * sum(rep(i, :));
Rf(i) = exp(-Rf(i) * T);
end

r1 = mean(Rf);

end

Finite Differences Method – Matlab Code

In this Matlab code we implement the following finite differences methods: implicit, explicit and Crank-Nicolson to estimate option values. We also plot our concluded values vs. the ones generated by the Black-Scholes closed-form solution.


function FiniteDifferences()

S0 = 10; sigma = 0.20;r = 0.04; dT = 0.002;X0 = log(S0);T =0.5
delX = [sigma*sqrt(dT),sigma*sqrt(3*dT),sigma*sqrt(4*dT)];
figure;
for i =1:3
dX = delX(i);
Path = ceil(log(18/2) / dX) ;

Time = (T/dT)+1;

X = zeros(Path,1);
P = zeros(Path,Time);

X(1) = log(2);

P(1,Time) = 10 - exp(X(1));
for j =2:Path
X(j) = X(j-1) + dX;
P(j,Time) = max(10 - exp(X(j)),0);
end
% Explicit different method
[p_d, p_m, p_u] = compute_explicit_prob(r,sigma,dX,dT);
A = GetExplicitA(p_u, p_m, p_d,Path);

PutValue = zeros(Path,3);

% Main Block<span></span>
for j = Time-1:-1:1
P(2:Path-1,j) = A*P(:,j+1);
P(1,j) = P(2,j) -(exp(X(1)) - exp(X(2)));
P(Path,j) = P(Path-1,j) ;
end
Y = zeros(size(X,1),1);
for m =1:size(X,1)
[~, Y(m)] = blsprice(exp(X(m)),S0,0.04,0.5,0.20,0);
end
subplot(3,3,(i-1)*3+1);
plot(exp(X),P(:,1),exp(X),Y);
legend('Explicit Difference Price','BS Price');

P(:,1:Time-1) = zeros(Path,Time-1);
[p_d, p_m, p_u] = compute_implicit_prob(r,sigma,dX,dT);
A_I = GetImplicitA(p_u, p_m, p_d,Path);
for j = Time-1:-1:1
b = [exp(X(2)) - exp(X(1)); P(2:Path-1,j+1);0 ];
P(:,j) = A_I\b;
end

subplot(3,3,(i-1)*3+2);
plot(exp(X),P(:,1),exp(X),Y);
legend('Implicit Difference Price','BS Price');

P(:,1:Time-1) = zeros(Path,Time-1);
[p_d, p_m, p_u] = compute_crank_nicol_prob(r,sigma,dX,dT);
Z = zeros(Path,1);
A_I = GetImplicitA(p_u, p_m, p_d,Path);
for j = Time-1:-1:1

Z(2:Path-1) = GetCrankZ(p_u, p_m, p_d,Path,P(:,j+1));
Z(1) = -(exp(X(1)) - exp(X(2)));
Z(Path) =Z(Path-1) ;
P(:,j) = A_I\Z;
end

subplot(3,3,(i-1)*3+3);
plot(exp(X),P(:,1),exp(X),Y);
legend('Crank Nicolson Price','BS Price');

end

end

function [p_u, p_m, p_d] = compute_explicit_prob(r,sigma,dX,dT)
t1 = (sigma^2)/(2*(dX^2));
t2 = (r - ((sigma^2)/2))/(2*dX);
t3 = dT * ((sigma^2)/(dX^2));
t4 = r*dT;
p_u = dT*(t1 +t2);
p_m = 1 - t3 -t4;
p_d = dT*(t1 -t2);
end

function [p_u, p_m, p_d] = compute_implicit_prob(r,sigma,dX,dT)
v = r - ((sigma^2)/2);
t1 = (sigma^2)/((dX^2));
t2 = (v)/(dX);
t3 = dT * ((sigma^2)/(dX^2));
t4 = r*dT;
p_u = -0.5* dT*(t1 +t2);
p_m = 1 + t3 + t4;
p_d = -0.5* dT*(t1 -t2);
end

function [p_u, p_m, p_d] = compute_crank_nicol_prob(r,sigma,dX,dT)
v = r - ((sigma^2)/2);
t1 = (sigma^2)/((dX^2));
t2 = (v)/(dX);
t3 = dT * ((sigma^2)/(2*dX^2));
t4 = r*dT*0.5;
p_u = -0.25* dT*(t1 +t2);
p_m = 1 + t3 + t4;
p_d = -0.25* dT*(t1 -t2);
end

function [A] = GetExplicitA(p_u, p_m ,p_d, N)
A = zeros(N-2,N);
A(1,1:3)= [p_u p_m p_d];
for i = 2:N-2
A(i,:) = circshift(A(i-1,:),[0,1]);
end

end
function [A] = GetImplicitA(p_u, p_m ,p_d, N)
A = zeros(N,N);
A(1,1:2) = [ 1 -1];
A(2,1:3)= [p_u p_m p_d];
for i = 3:N-1
A(i,:) = circshift(A(i-1,:),[0,1]);
end
A(N,N-1:N) = [ 1 -1];
end

function [A] = GetCrankZ(p_u, p_m ,p_d, N,P)
A = zeros(N-2,N);
A(1,1:3)= [-p_u -(p_m-2) -p_d];
for i = 2:N-2
A(i,:) = circshift(A(i-1,:),[0,1]);
end
A = A*P;
end