Commit f63934ea authored by Enrico Schumann's avatar Enrico Schumann

Initial commit

parents
*~
.gitignore
function x1 = BFGS(f,x1,tol,ls)
% BFGS.m -- version 2010-09-14
if nargin == 2, tol = 1e-4; ls = 0.1; end
g1 = numG(f,x1); B1 = eye(numel(x1));
x1 = x1(:); x0 = -x1; k = 1;
while ~converged(x0,x1,tol)
x0 = x1; g0 = g1; B0 = B1;
p0 = -B0\g0;
as = lineSearch(f,x0,p0,ls);
s0 = as * p0;
z0 = B0 * s0;
x1 = x0 + s0;
g1 = numG(f,x1);
y0 = g1 - g0;
B1 = B0 + (y0*y0')/(y0'*s0) - (z0*z0')/(s0'*z0);
k = k + 1; if k > 100, error('Maxit in BFGS'); end
end
function c = Bisection(f,a,b,tol)
% Bisection.m -- version 2005-05-25
% Zero finding with bisection method
if nargin == 3, tol = 1e-8; end
fa = feval(f,a); fb = feval(f,b);
if sign(fa) == sign(fb)
error('sign f(a) is not opposite of sign f(b)');
end
done = 0;
while abs(b-a) > 2*tol & ~done
c = a + (b - a) / 2;
fc = feval(f,c);
if sign(fa) ~= sign(fc)
b = c;
fb = fc;
elseif sign(fc) ~= sign(fb)
a = c;
fa = fc;
else % center and zero coincide
done = 1;
end
end
function I = Bracketing(f,xL,xU,n)
% Bracketing.m version 2010-03-27
delta = (xU - xL)/n;
a = xL; k = 0; I = [];
fa = feval(f,a);
for i = 1:n
b = a + delta;
fb = feval(f,b);
if sign(fa)~=sign(fb)
k = k + 1;
I(k,:) = [a b];
end
a = b; fa = fb;
end
\ No newline at end of file
function y1 = Broyden(f,y1,B,tol,itmax)
% Broyden.m -- version 2010-08-10
n = numel(y1);
if nargin < 3, B = eye(n,n); tol = 1e-2; itmax = 10; end
y1 = y1(:); y0 = -y1; k = 0;
F1 = feval(f,y1);
while ~converged(y0,y1,tol)
y0 = y1;
F0 = F1;
s = B \ -F0;
y1 = y0 + s;
F1 = feval(f,y1);
dF = F1 - F0;
B = B + ((dF - B*s)*s')/(s'*s);
k = k + 1;
if k > itmax, error('Iteration limit reached'); end
end
function x1 = FPI(f,x1,tol)
% FPI.m -- version 2008-05-08
% FPI(f,x0) Fixed point iteration
if nargin == 2, tol = 1e-6; end
it = 0; itmax = 200; x0 = realmax;
while ~converged(x0,x1,tol)
x0 = x1;
x1 = feval(f,x0);
it = it + 1;
if it > itmax, error('Maxit in FPI'); end
end
% FPInD.m -- version 2010-12-23
G = str2mat('-y0(2) + 3','sqrt(-y0(1)^2 + 9)');
y1 = [2 0]; y0 = -y1; tol = 1e-2; k = 1; itmax = 10;
Y = NaN(itmax,2); Y(1,:) = y1';
while ~converged(y0,y1,tol)
y0 = y1;
y1(1) = eval(G(1,:));
y1(2) = eval(G(2,:));
k = k + 1;
Y(k,:) = y1';
if k > itmax, error('Iteration limit reached'); end
end
y1
function [f,g] = FnG(x,func)
% FnG.m -- version 2010-12-21
f = func(x);
g = numG(func,x);
\ No newline at end of file
function c = GSS(f,a,b,tol)
% GSS.m -- version 2006-05-26
% Golden Section Search
if nargin == 3, tol = 1e-8; end
tau = (sqrt(5) - 1)/2;
x1 = a + (1-tau)*(b-a); f1 = feval(f,x1);
x2 = a + tau *(b-a); f2 = feval(f,x2);
while (b-a) > tol
if f1 < f2
b = x2;
x2 = x1; f2 = f1;
x1 = a + (1-tau)*(b-a); f1 = feval(f,x1);
else
a = x1;
x1 = x2; f1 = f2;
x2 = a + tau *(b-a); f2 = feval(f,x2);
end
end
c = a + (b-a)/2;
function x1 = Newton0(f,x1,tol)
% Newton0.m -- version 2010-03-30
% Newton method with numerical derivation
if nargin == 2, tol = 1e-8; end
it = 0; itmax = 10; x0 = realmax; h = 1e-8;
while ~converged(x0,x1,tol)
x0 = x1;
f0 = feval(f,x0);
f1 = feval(f,x0 + h);
df = (f1 - f0) / h;
x1 = x0 - f0/df;
it = it + 1;
if it > itmax, error('Maxit in Newton0'); end
end
\ No newline at end of file
function x1 = NewtonnD(f,x1,tol,itmax)
% NewtonnD.m -- version 2010-09-12
if nargin < 3, tol = 1e-2; itmax = 10; end
x1 = x1(:); x0 = -x1; k = 0;
while ~converged(x0,x1,tol)
x0 = x1;
F = feval(f,x0);
b = -F; J = numJ(f,x0);
s = J \ b;
x1 = x0 + s;
k = k + 1;
if k > itmax, error('Maxit in NewtonnD'); end
end
function F = OF(theta,y,X,b,c)
% OF.m -- version 2010-11-11
r = y - X*theta;
S0 = median(abs(r)) / b;
F = FPIS(S0,r,c);
function F = Rho(s)
% Rho.m -- version 2010-11-12
k = 1.58;
F = ones(length(s),1);
I = (abs(s) <= k);
temp1 = s(I)/k;
temp2 = temp1 .* temp1;
temp4 = temp2 .* temp2;
temp6 = temp4 .* temp2;
F(I) = 3*temp2 - 3*temp4 + temp6;
function F = RhoN(s)
% RhoN.m -- version 2010-11-12
d = exp(-0.5 * s.^2) ./ sqrt(2*pi);
F = Rho(s) .* d';
\ No newline at end of file
function x1 = SteepestD(f,x1,tol,ls)
% SteepestD.m -- version2010-08-31
if nargin == 2, tol = 1e-4; ls = 0.1; end
x1 = x1(:); x0 = -x1; k = 1;
while ~converged(x0,x1,tol)
x0 = x1;
g = numG(f,x1);
as = lineSearch(f,x0,-g,ls);
x1 = x0 - as * g;
k = k + 1; if k > 300, error('Maxit in SteepestD'); end
end
function F = ex2NLeqs(x)
% ex2NLeqs.m -- version 2010-12-23
F = zeros(numel(x),1);
F(1) = x(1) + x(2) - 3;
F(2) = x(1)^2 + x(2)^2 - 9;
\ No newline at end of file
function a1 = lineSearch(f,x,s,d)
% lineSearch.m -- version 2010-09-13
if nargin == 3, d = 0.1; end
done = 0; a1 = 0; k = 0;
f1 = feval(f,x);
while ~done
a0 = a1; a1 = a0 + d;
f0 = f1; f1 = feval(f,x+a1*s);
if f1 > f0, done = 1; end
k = k + 1;
if k > 100, fprintf('Early stop in lineSearch'); break, end
end
function g = numG(f,x)
% numG.m -- version 2010-08-31
n = numel(x); g = zeros(n,1); x = x(:);
Delta = diag(max(sqrt(eps) * abs(x),sqrt(eps)));
F0 = feval(f,x);
for i = 1:n
F1 = feval(f,x + Delta(:,i));
g(i) = (F1 - F0) / Delta(i,i); % forward difference
end
function J = numJ(f,x)
% numJ.m -- version 2010-09-16
n = numel(x); J = zeros(n,n); x = x(:);
Delta = diag(max(sqrt(eps) * abs(x),sqrt(eps)));
F = feval(f,x);
for j = 1:n
Fd = feval(f,(x + Delta(:,j)));
J(:,j) = (Fd - F) / Delta(j,j);
end
function C0 = AmericanCallDiv(S0,X,r,T,sigma,D,TD,M)
% AmericanCallDiv.m -- version 2010-12-28
% compute constants
f7 = 1; dt = T/M; v = exp(-r * dt);
u = exp(sigma*sqrt(dt)); d = 1 /u;
p = (exp(r * dt) - d)/(u - d);
% adjust spot for dividend
S0 = S0 - D * exp(-r * TD);
% initialize asset prices at maturity (period M)
S = zeros(M + 1,1);
S(f7+0) = S0 * d^M;
for j = 1:M
S(f7+j) = S(f7+j - 1) * u / d;
end
% initialize option values at maturity (period M)
C = max(S - X, 0);
% step back through the tree
for i = M-1:-1:0
% compute present value of dividend (dD)
t = T * i / M;
dD = D * exp(-r * (TD-t));
for j = 0:i
C(f7+j) = v * ( p * C(f7+j + 1) + (1-p) * C(f7+j));
S(f7+j) = S(f7+j) / d;
if t > TD
C(f7+j) = max(C(f7 + j), S(f7+j) - X);
else
C(f7+j) = max(C(f7 + j), S(f7+j) + dD - X);
end
end
end
C0 = C(f7+0);
\ No newline at end of file
function P0 = AmericanPut(S0,X,r,T,sigma,M)
% AmericanPut.m -- version 2010-12-28
f7 = 1; dt = T / M; v = exp(-r * dt);
u = exp(sigma * sqrt(dt)); d = 1 / u;
p = (exp(r * dt) - d) / (u - d);
% initialize asset prices at maturity (period M)
S = zeros(M + 1,1);
S(f7+0) = S0 * d^M;
for j = 1:M
S(f7+j) = S(f7+j - 1) * u / d;
end
% initialize option values at maturity (period M)
P = max(X - S, 0);
% step back through the tree
for i = M-1:-1:0
for j = 0:i
P(f7+j) = v * (p * P(f7+j + 1) + (1-p) * P(f7+j));
S(f7+j) = S(f7+j) / d;
P(f7+j) = max(P(f7 + j),X - S(f7+j));
end
end
P0 = P(f7+0);
\ No newline at end of file
function C0 = EuropeanCall(S0,X,r,T,sigma,M)
% EuropeanCall.m -- version 2010-12-28
% compute constants
f7 = 1; dt = T / M; v = exp(-r * dt);
u = exp(sigma*sqrt(dt)); d = 1 /u;
p = (exp(r * dt) - d) / (u - d);
% initialize asset prices at maturity (period M)
S = zeros(M + 1,1);
S(f7+0) = S0 * d^M;
for j = 1:M
S(f7+j) = S(f7+j - 1) * u / d;
end
% initialize option values at maturity (period M)
C = max(S - X, 0);
% step back through the tree
for i = M-1:-1:0
for j = 0:i
C(f7+j) = v * (p * C(f7+j + 1) + (1-p) * C(f7+j));
end
end
C0 = C(f7+0);
\ No newline at end of file
function C0 = EuropeanCallBE(S0,X,r,T,sigma,M)
% EuropeanCallBE.m -- version: 2011-03-14
% compute constants
dt = T / M;
u = exp(sigma*sqrt(dt)); d = 1 /u;
p = (exp(r * dt) - d) / (u - d);
% initialise asset prices at maturity (period M)
C = max(S0*d.^((M:-1:0)').*u.^((0:M)') - X,0);
% log/cumsum version
csl = cumsum(log([1;[1:M]']));
tmp = csl(M+1) - csl - csl(M+1:-1:1) + log(p)*((0:M)') + ...
log(1-p)*((M:-1:0)');
C0 = exp(-r*T)*sum(exp(tmp).*C);
\ No newline at end of file
function [C0,deltaE,gammaE,thetaE] = EuropeanCallGreeks(S0,X,r,T,sigma,M)
% EuropeanCallGreeks.m -- version 2010-12-28
% compute constants
f7 = 1; dt = T/M; v = exp(-r * dt);
u = exp(sigma*sqrt(dt)); d = 1 /u;
p = (exp(r * dt) - d)/(u - d);
% initialize asset prices at maturity (period M)
S = zeros(M + 1,1);
S(f7+0) = S0 * d^M;
for j = 1:M
S(f7+j) = S(f7+j - 1) * u / d;
end
% initialize option values at maturity (period M)
C = max(S - X, 0);
% step back through the tree
for i = M-1:-1:0
for j = 0:i
C(f7+j) = v * ( p * C(f7+j + 1) + (1-p) * C(f7+j));
S(f7+j) = S(f7+j) / d;
end
if i==2
%gamma
gammaE = ((C(2+f7) - C(1+f7)) / (S(2+f7) - S(1+f7)) - ...
(C(1+f7) - C(0+f7)) / (S(1+f7) - S(0+f7)))/ ...
(0.5 * (S(2+f7) - S(0+f7)));
%theta (aux)
thetaE = C(1+f7);
end
if i==1
%delta
deltaE = (C(1+f7) - C(0+f7)) / (S(1+f7) - S(0+f7));
end
if i==0
%theta (final)
thetaE = (thetaE - C(0+f7)) / (2 * dt);
end
end
C0 = C(f7+0);
\ No newline at end of file
# EuropeanCall.R -- version 2010-12-28
EuropeanCall <- function(S0,X,r,T,sigma,M) {
# compute constants
f7 <- 1; dt <- T/M; v <- exp(-r * dt)
u <- exp(sigma * sqrt(dt)); d <- 1/u
p <- (exp(r * dt) - d) / (u - d)
# initialize asset prices at maturity (period M)
S <- numeric(M + 1)
S[f7+0] <- S0 * d^M
for (j in 1:M) S[f7+j] <- S[f7+j - 1] * u / d
# initialise option values at maturity (period M)
C <- pmax(S - X, 0)
# step back through the tree
for (i in seq(M-1, 0, by = -1)){
C <- v * (p * C[(1+f7):(i+1+f7)] +
(1 - p)* C[(0+f7):(i+f7)])
}
C
}
# EuropeanCallBE.R -- version 2010-12-28
EuropeanCallBE <- function(S0, X, r, T, sigma, M) {
# compute constants
dt <- T/M
u <- exp(sigma * sqrt(dt))
d <- 1/u
p <- (exp(r * dt) - d)/(u - d)
# initialize asset prices at maturity (period M)
C <- pmax(S0 * d^(M:0) * u^(0:M) - X, 0)
# log/cumsum version
csl <- cumsum(log(c(1,1:M)))
tmp <- csl[M+1] - csl - csl[(M+1):1] +
log(p) * (0:M) + log(1-p) * (M:0)
C0 <- exp(-r * T) * sum(exp(tmp) * C)
C0
}
function [V C B F E] = CPPIgap(S, m, G, r_c, gap)
% CPPIgap.m version -- 2010-12-22
% S .. stock price series t = 0 .. T
% G .. guarantueed payback amount
% m .. multiplier
% r_c .. cumulated save return over entire horizon
% gap .. readjustment frequency; if blank: 1 = always
if nargin < 5, gap = 1; end;
% --initial setting
T = length(S)-1;
t = 0:T;
V = zeros(T+1,1);
V(1) = G;
F = G*exp(-r_c*((T-t)/T))';
C = zeros(T+1,1);
B = zeros(T+1,1);
n = zeros(T+1,1); % number of risky assets
% --development over time
for tau=1:T % tau = t+1
C(tau) = V(tau)-F(tau);
if mod(tau-1,gap)==0 % re-adjust now
E(tau) = min(m * C(tau), V(tau));
n(tau) = E(tau) / S(tau);
B(tau) = V(tau) - E(tau);
else
n(tau) = n(tau-1);
E(tau) = V(tau) - B(tau);
end;
B(tau+1) = B(tau)*exp(r_c/T);
V(tau+1) = n(tau)*S(tau+1) + B(tau+1);
end
\ No newline at end of file
function vv = VDC(k,b)
% VDC.m -- version 2010-12-08
nD = 1 + floor(log(max(k))/log(b)); % required digits
nN = length(k); % number of VDC numbers
vv = zeros(nN,nD);
for i = nD:-1:1
vv(:,i) = mod(k,b);