function []=blatt1()

%Parameter
r=0.0001;  
T=30;         %Tage
K=0.5;        %Strike
Q=1;          %Payoff
mu=0.01;      %Erwartungswert
sigma=0.1;    %Standardabweichung
S0=0.43;         %Startwert
n=10;         %Anzahl normalverteilter ZVe

%Wert: V = exp(-rT)*Q*p
%Aktienwert S_T = S0 exp(sT), s normalverteilt bzgl  (mu- sigma^2/2 , sigma^2/T)

n=[10,100,1000,10000];
S0=(0:0.001:1)';

for i=1:4
     p(:,i) = wkeit(r,T,K,Q,mu,sigma,S0,n(i));
end

%approximativer Wert
Vapprox = exp(-r*T )*Q*p;

%exakter Wert
d = ( log(S0/K) + T*(mu - sigma^2/2) )/(sigma*sqrt(T));
phid = normcdf(d);             %normal cumulative distribution function
V = exp( -r*T ) *Q*phid ;

%Werte
figure
hold on

for i=1:4
  plot(S0,Vapprox(:,i))
end

plot(S0,V,'r')

title('Riskioabschaetzung einer binaerer Option')
xlabel('S0')
ylabel('V')
legend('n=10','n=100','n=1000','n=10000','exakt')

%Finite Differenzen
figure
hold on

for i=1:4
   Vapproxdiff(:,i) = ( Vapprox(2:end,i) - Vapprox(1:end-1,i))/0.001;
   plot(S0(1:end-1),Vapproxdiff(:,i))
end

Vdiff = (V(2:end) - V(1:end-1))/0.001;
plot(S0(1:end-1),Vdiff,'r','Linewidth',1.5)
title('Risikoabschaetzung - finite Differenzen')
xlabel('S0')
ylabel('dV/dS0')
legend('n=10','n=100','n=1000','n=10000','exakt')
%=> Differentiation schlecht gestellt

end

function p=wkeit(r,T,K,Q,mu,sigma,S0,n)

mu = mu - sigma^2/2;
sigma = sqrt( sigma^2/T );

%Erzeugen der ZVe
R = normrnd(mu,sigma,n,1);
R =  exp(R*T);
S = R*S0';

for i=1:length(S0)
   p(i) = sum( S(:,i) >= K )/n;    % R >= K ist ein Vektor, der in der i-ten Komponente eine 1 enthält, wenn R(i) >= K ist und ansonsten eine 0
end

end