x=-3:0.1:3;
plot(x,normcdf(x))
title("CDF of Normal Distribution")
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(-3,3.1, step = 0.1) # Python's arange excludes the last value
plt.plot(x, stats.norm.cdf(x))
plt.title("CDF of Normal distribution")
plt.show()
plt.close()
rng default; % set seed
S=10;
rand(S,1)
randn(S,1)
trnd(4,S,1)
np.random.seed(12) # set seed
S = 10 # simulation size
print (np.random.uniform(size = S))
print (np.random.normal(size = S))
print (np.random.standard_t(df = 4,size = S))
yield = [5.00 5.69 6.09 6.38 6.61...
6.79 6.94 7.07 7.19 7.30]; % yield curve
T = length(yield);
r=0.07; % initial yield rate
Par=10; % par value
coupon=r*Par; % coupon payments
cc=zeros(1,T)+coupon; % vector of cash flows
cc(T)=cc(T)+Par; % add par to cash flows
P=sum(cc./((1+yield./100).^(1:T))) % calculate price
yield_c = [5.00, 5.69, 6.09, 6.38, 6.61,
6.79, 6.94, 7.07, 7.19, 7.30] # yield curve
T = len(yield_c)
r = 0.07 # initial yield rate
Par = 10 # par value
coupon = r * Par # coupon payments
cc = [coupon] * 10 # vector of cash flows
cc[T-1] += Par # add par to cash flows
P=np.sum(cc/(np.power((1+np.divide(yield_c,100)),
list(range(1,T+1))))) # calc price
print(P)
randn('state',123); % set the seed
sigma = 1.5; % daily yield volatility
S = 8; % number of simulations
r = randn(1,S)*sigma; % generate random numbers
ysim=nan(T,S);
for i=1:S
ysim(:,i)=yield+r(i);
end
ysim=repmat(yield',1,S)+repmat(r,T,1);
plot(ysim)
title("Simulated yield curves")
np.random.seed(12) # set seed
sigma = 1.5 # daily yield volatility
S = 8 # number of simulations
r = np.random.normal(0,sigma,size=S) # generate random numbers
ysim = np.zeros([T,S])
for i in range(S):
ysim[:,i] = yield_c + r[i]
plt.plot(ysim)
plt.title("Simulated yield curves")
plt.show()
plt.close()
SP = nan(S,1);
for s = 1:S % S simulations
SP(s) = sum(cc./(1+ysim(:,s)'./100).^((1:T)));
end
SP = SP-(mean(SP) - P); % correct for mean
bar(SP)
S = 50000;
rng("default")
r = randn(S,1) * sigma;
ysim = nan(T,S);
for i = 1:S
ysim(:,i) = yield' + r(i);
end
SP = nan(S,1);
for i = 1:S
SP(i) = sum(cc./(1+ysim(:,i)'./100).^((1:T)));
end
SP = SP - (mean(SP)-P);
histfit(SP)
title("Histogram of simulated bond prices with fitted normal")
S = 8 # simulation size
SP = np.zeros([S]) # initializing vector with zeros
for i in range(S): # S simulations
SP[i] = np.sum(cc/(np.power((1+ysim[:,i]/100), list(range(1,T+1)))))
SP -= (np.mean(SP) - P) # correct for mean
plt.bar(range(1,S+1), SP)
plt.title("Simulated bond prices, S = 8")
plt.show()
plt.close()
S = 50000
np.random.seed(12)
r = np.random.normal(0, sigma, size = S)
ysim = np.zeros([T,S])
for i in range(S):
ysim[:,i] = yield_c + r[i]
SP = np.zeros([S])
for i in range(S):
SP[i] = np.sum(cc/(np.power((1+ysim[:,i]/100), list(range(1,T+1)))))
SP -= (np.mean(SP) - P)
plt.hist(SP, bins = 30, range = (7, 13), density = True)
fitted_norm=stats.norm.pdf(np.linspace(7,13,30),
np.mean(SP),np.std(SP,ddof=1))
plt.plot(np.linspace(7,13,30), fitted_norm)
plt.title("Histogram for simulated bond prices, S = 50000")
plt.show()
plt.close()
P0 = 50; % initial spot price
sigma = 0.2; % annual volatility
r = 0.05; % annual interest
T = 0.5; % time to expiration
X = 40; % strike price
f = bs(X,P0,r,sigma,T) % analytical call price
%% This calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
## This calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
P0 = 50 # initial spot price
sigma = 0.2 # annual volatility
r = 0.05 # annual interest
T = 0.5 # time to expiration
X = 40 # strike price
f = bs(X, P0, r, sigma, T) # analytical call price
print(f)
randn('state',0); % set seed
S = 1e6; % number of simulations
ysim = randn(S,1)*sigma*sqrt(T)-0.5*T*sigma^2; % sim returns, lognorm corrected
F = P0*exp(r*T)*exp(ysim); % sim future prices
SP = F-X; % payoff
SP(find(SP < 0)) = 0; % set negative outcomes to zero
fsim = SP * exp(-r*T) ; % discount
mean(fsim) % simulated price
np.random.seed(12) # set seed
S = 10**6 # number of simulations
ysim=np.random.normal(-0.5*sigma**2*T,
sigma*np.sqrt(T),size=S) # sim returns, lognorm corrected
F = P0 * np.exp(r * T) * np.exp(ysim) # sim future prices
SP = F - X # payoff
SP[SP < 0] = 0 # set negative outcomes to zero
fsim = SP * np.exp(-r * T) # discount
call_sim = np.mean(fsim) # simulated price
print(call_sim)
subplot(1,2,1)
histfit(F);
title("Simulated prices");
xline(X, 'LineWidth', 1, 'label', 'Strike');
subplot(1,2,2)
hist(fsim,100);
title("Option price density");
xline(mean(fsim), 'LineWidth', 1, 'label', 'Call');
plt.hist(F, bins = 60, range = (20,80), density = True)
fitted_norm=stats.norm.pdf(np.linspace(20,80,60),np.mean(F),np.std(F,ddof=1))
plt.plot(np.linspace(20,80,60), fitted_norm)
plt.axvline(x=X, color='k')
plt.title("Density of future prices")
plt.show()
plt.close()
plt.hist(fsim, bins = 60, range = (0, 35), density = True)
plt.axvline(x=f['Call'], color='k')
plt.title("Density of Call option prices")
plt.show()
plt.close()
randn('state',0); % set seed
S = 1e7; % number of simulations
s2 = 0.01^2; % daily variance
p = 0.01; % probability
r = 0.05; % annual riskfree rate
P = 100; % price today
ysim = randn(S,1)*sqrt(s2)+r/365-0.5*s2; % sim returns
Psim = P*exp(ysim); % sim future prices
q = sort(Psim-P); % simulated P/L
VaR1 = -q(S*p)
from math import ceil
np.random.seed(1) # set seed
S = 10**7 # number of simulations
s2 = 0.01**2 # daily variance
p = 0.01 # probability
r = 0.05 # annual riskfree rate
P = 100 # price today
ysim=np.random.normal(r/365-0.5*s2,np.sqrt(s2),size=S) # sim returns
Psim = P * np.exp(ysim) # sim future prices
q = np.sort(Psim - P) # simulated P/L
VaR1 = -q[ceil(p*S) - 1]
print(VaR1)
T = 0.25; % time to expiration
X = 100; % strike price
sigma = sqrt(s2*250); % annual volatility
f = bs(X,P,r,sigma,T); % analytical call price
fsim=bs(X,Psim,r,sigma,T-(1/365)); % sim option prices
q = sort(fsim.Call-f.Call); % simulated P/L
VaR2 = -q(p*S)
T = 0.25 # time to expiration
X = 100 # strike price
sigma = np.sqrt(s2 * 250) # annual volatility
f = bs(X, P, r, sigma, T) # analytical call price
fsim = bs(X, Psim, r, sigma, T-(1/365)) # sim option prices
q = np.sort(fsim['Call']-f['Call']) # simulated P/L
VaR2 = -q[ceil(p*S) - 1]
print(VaR2)
X1 = 100;
X2 = 110;
f1 = bs(X1,P,r,sigma,T);
f2 = bs(X2,P,r,sigma,T);
f1sim=bs(X1,Psim,r,sigma,T-(1/365));
f2sim=bs(X2,Psim,r,sigma,T-(1/365));
q = sort(f1sim.Call+f2sim.Put+Psim-f1.Call-f2.Put-P);
VaR3 = -q(p*S)
X1 = 100
X2 = 110
f1 = bs(X1, P, r, sigma, T)
f2 = bs(X2, P, r, sigma, T)
f2sim = bs(X2, Psim, r, sigma, T-(1/365))
f1sim = bs(X1, Psim, r, sigma, T-(1/365))
q = np.sort(f1sim['Call'] + f2sim['Put'] + Psim - f1['Call'] - f2['Put'] - P)
VaR3 = -q[ceil(p*S) - 1]
print(VaR3)
randn('state',12) % set seed
mu = [r/365 r/365]'; % return mean
Sigma=[0.01 0.0005; 0.0005 0.02]; % covariance matrix
y = mvnrnd(mu,Sigma,S); % simulated returns
np.random.seed(12) # set seed
mu = np.transpose([r/365, r/365]) # return mean
Sigma = np.matrix([[0.01, 0.0005],[0.0005, 0.02]]) # covariance matrix
y = np.random.multivariate_normal(mu, Sigma, size = S) # simulated returns
K = 2;
P = [100 50]'; % prices
x = [1 1]'; % number of assets
Port = P'*x; % portfolio at t
Psim = repmat(P,1,S)' .*exp(y); % simulated prices
PortSim=Psim * x; % simulated portfolio value
q = sort(PortSim-Port); % simulated P/L
VaR4 = -q(S*p)
import numpy.matlib
P = np.asarray([100, 50]) # prices
x = np.asarray([1, 1]) # number of assets
Port = np.matmul(P, x) # portfolio at t
Psim = np.matlib.repmat(P,S,1)*np.exp(y) # simulated prices
PortSim = np.matmul(Psim, x) # simulated portfolio value
q = np.sort(PortSim - Port) # simulated P/L
VaR4 = -q[ceil(p*S) - 1]
print(VaR4)
f = bs(P(2),P(2),r,sigma,T);
fsim=bs(P(2),Psim(:,2),r,sigma,T-(1/365));
q = sort(fsim.Call+Psim(:,1)-f.Call-P(1));
VaR5 = -q(p*S)
f = bs(X = P[1], P = P[1], r = r, sigma = sigma, T = T)
fsim = bs(X = P[1], P = Psim[:,1], r = r, sigma = sigma, T = T-(1/365))
q = np.sort(fsim['Call'] + Psim[:,0] - f['Call'] - P[0])
VaR5 = -q[ceil(p*S) - 1]
print(VaR5)