x=seq(-3,3,by=0.1)
plot(x,pnorm(x),type="l")
x = collect(linspace(-3, 3, 61))
using Plots;
using Distributions;
return plot(x, cdf.(Normal(0,1), x))
set.seed(12) # set seed
S=10
runif(S)
rnorm(S)
rt(S,4)
srand(12); # set seed
S = 10;
println(rand(Uniform(0,1), S)) # alternatively, rand(S)
println(rand(Normal(0,1), S)) # alternatively, randn(S)
println(rand(TDist(4), S))
yield=c(5.00, 5.69, 6.09, 6.38, 6.61,
6.79, 6.94, 7.07, 7.19, 7.30) # yield curve
r=0.07 # initial yield rate
Par=10 # par value
coupon=r*Par # coupon payments
cc=1:10*0+coupon # vector of cash flows
cc[10]=cc[10]+Par # add par to cash flows
P=sum(cc/((1+yield/100)^(1:length(yield)))) # calculate price
print(P)
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 = length(yield_c)
r = 0.07 # initial yield rate
Par = 10 # par value
coupon = r * Par # coupon payments
cc = repeat([coupon], outer = 10) # vector of cash flows
cc[10] += Par # add par to cash flows
P = sum(cc./((1+yield_c/100).^(1:T))) # calc price
set.seed(12) # set seed
sigma = 1.5 # daily yield volatiltiy
S = 8 # number of simulations
r = rnorm(S,0,sigma) # generate random numbers
ysim = matrix(nrow=length(yield),ncol=S)
for (i in 1:S) ysim[,i]=yield+r[i]
matplot(ysim,type='l')
srand(12) # set seed
sigma = 1.5 # daily yield volatility
S = 8 # number of simulations
r = rand(Normal(0,1), S) # generate random numbers
ysim = fill!(Array{Float64}(T,S), NaN)
for i in range(1, S)
ysim[:,i] = yield_c + r[i]
end
using Plots;
plot(ysim)
SP = vector(length=S)
for (i in 1:S){ # S simulations
SP[i] = sum(cc/((1+ysim[,i]/100)^(length(yield))))
}
SP = SP-(mean(SP) - P) # correct for mean
par(mfrow=c(1,2), pty="s")
barplot(SP)
hist(SP,probability=TRUE)
x=seq(6,16,length=100)
lines(x, dnorm(x, mean = mean(SP), sd = sd(SP)))
SP = fill!(Array{Float64}(S,1), NaN)
for i in range(1,S) # S simulations
SP[i] = sum(cc./(1+ysim[:,i]./100).^T)
end
SP -= (mean(SP) - P) # correct for mean
using Plots;
bar(SP)
S = 50000
r = randn(S) * sigma
ysim = fill!(Array{Float64}(T,S), NaN)
for i in range(1, S)
ysim[:,i] = yield_c + r[i]
end
SP = fill!(Array{Float64}(S,1), NaN)
for i in range(1,S)
SP[i] = sum(cc./(1+ysim[:,i]./100).^T)
end
SP -= (mean(SP) - P)
using Plots;
histogram(SP,nbins=100,normed=true,xlims=(7,13))
res = fit_mle(Normal, SP)
plot!(Normal(res.μ, res.σ), linewidth = 4)
P0 = 50 # initial spot price
sigma = 0.2 # annual volatility
r = 0.05 # annual interest
TT = 0.5 # time to expiration
X = 40 # strike price
f = bs(X,P0,r,sigma,length(yield)) # analytical call price
## this calculation uses the Black-Scholes pricing function (Listing 6.1/6.2)
print(f)
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)
set.seed(12) # set seed
S = 1e6 # number of simulations
F = P0*exp(r*length(yield)) # futures price
ysim = rnorm(S,-0.5*sigma^2*length(yield),sigma*sqrt(length(yield))) # sim returns, lognorm corrected
F=F*exp(ysim) # sim futures price
SP = F-X # payoff
SP[SP<0] = 0 # set negative outcomes to zero
fsim = SP*exp(-r*length(yield)) # discount
call_sim = mean(fsim) # simulated price
print(call_sim)
srand(12) # set seed
S = 10^6 # number of simulations
F = P0 * exp(r * T) # futures price
ysim = randn(S)*sigma*sqrt(T)-0.5*sigma^2*T # sim returns, lognorm corrected
F = F * exp.(ysim) # sim futures price
SP = F - X # payoff
SP[SP.<0] = 0 # set negative outcomes to zero
fsim = SP * exp(-r * T) # discount
call_sim = mean(fsim) # simulated price
par(mfrow=c(1,2), pty="s")
hist(F,probability=TRUE,ylim=c(0,0.06))
x=seq(min(F),max(F),length=100)
lines(x, dnorm(x, mean = mean(F), sd = sd(SP)))
hist(fsim,nclass=100,probability=TRUE)
using Plots;
histogram(F, bins = 100, normed = true, xlims = (20,80))
res = fit_mle(Normal, F)
plot!(Normal(res.μ, res.σ), linewidth = 4)
vline!([X], linewidth = 4, color = "black")
histogram(fsim, bins = 110, normed = true, xlims = (0,35))
vline!([f["Call"]], linewidth = 4, color = "black")
set.seed(1) # 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 = rnorm(S,r/365-0.5*s2,sqrt(s2)) # sim returns
Psim = P*exp(ysim) # sim future prices
q = sort(Psim-P) # simulated P/L
VaR1 = -q[p*S]
print(VaR1)
srand(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 = randn(S) * 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[convert(Int, p*S)]
TT = 0.25; # time to expiration
X = 100; # strike price
sigma = sqrt(s2*250); # annual volatility
f = bs(X,P,r,sigma,length(yield)) # analytical call price
fsim = bs(X,Psim,r,sigma,length(yield)-(1/365)) # sim option prices
q = sort(fsim$Call-f$Call) # simulated P/L
VaR2 = -q[p*S]
print(VaR2)
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[convert(Int, p*S)]
X1 = 100
X2 = 110
f1 = bs(X1,P,r,sigma,TT)
f2 = bs(X2,P,r,sigma,TT)
f2sim = bs(X2,Psim,r,sigma,TT-(1/365))
f1sim = bs(X1,Psim,r,sigma,TT-(1/365))
q = sort(f1sim$Call+f2sim$Put+Psim-f1$Call-f2$Put-P);
VaR3 = -q[p*S]
print(VaR3)
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 = sort(f1sim["Call"] + f2sim["Put"] + Psim - f1["Call"] - f2["Put"] - P)
VaR2 = -q[convert(Int, p*S)]
library (MASS)
set.seed(12) # set seed
mu = c(r/365,r/365) # return mean
Sigma = matrix(c(0.01,0.0005,0.0005,0.02),ncol=2) # covariance matrix
y = mvrnorm(S,mu,Sigma) # simulated returns
using Distributions;
srand(12); # set seed
mu = Vector([r/365, r/365]); # return mean
Sigma = [0.01 0.0005; 0.0005 0.02]; # covariance matrix
y = rand(MvNormal(mu,Sigma), S); # simulated returns
K=2
P = c(100,50) # prices
x = c(1,1) # number of assets
Port = P %*% x # portfolio at t
Psim = matrix(t(matrix(P,K,S)),ncol=K)*exp(y) # simulated prices
PortSim = Psim %*% x # simulated portfolio value
q = sort(PortSim-Port[1,1]) # simulated P/L
VaR4 = -q[S*p]
print(VaR4)
K = 2
P = [100 50] # prices
x = [1 1] # number of assets
Port = reshape(P * x', 1)[1] # portfolio at t
Psim = repmat(P, S, 1).*exp.(y)' # simulated prices
PortSim = reshape(Psim * x', S) # simulated portfolio value
q = sort(PortSim - Port) # simulated P/L
VaR4 = -q[convert(Int, p * S)]
f = bs(P[2],P[2],r,sigma,TT)
fsim = bs(P[2],Psim[,2],r,sigma,TT-(1/365))
q = sort(fsim$Call+Psim[,1]-f$Call-P[1]);
VaR5 = -q[p*S]
print(VaR5)
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[convert(Int, p * S)]