% Exercise sheet 4
% SDE dX_t = a X_t dt + b dB_t 
% Simple Euler Method



%% nur eine Realisierung
clear
close all

x0 = 0;
t(1) = 0;
dt = 0.01; % timestep

N = 100000;

a = -1;
b = 1;

mu = exp(a*N*dt)*x0
sigma = sqrt(b^2/(2*a)*(exp(2*a*N*dt)-1))


X(1) = x0; % initial value
for n = 2:N
    r = rand; 
    xi = -1+2*round(r); % P(xi=1)=P(xi=-1)=0.5
    X(n) = X(n-1) + a*X(n-1)*dt + b*sqrt(dt)*xi; % simple euler method for SDE
    t(n) = t(n-1) + dt;
    Y(n) = normrnd(mu,sigma);
end


histfit(X)

figure
histfit(Y)

figure
plot(t,X)
hold on
xlabel('t')
ylabel('X_t')

%% mehrere Realisierungen

clear
close all

x0 = 0;
t(1) = 0;
dt = 0.01; % timestep

N = 100000;
M = 100;

a = -1;
b = 1;

mu = exp(a*N*dt)*x0;
sigma = b^2/(2*a)*(exp(2*a*N*dt)-1);

for m = 1:M
    X(1) = x0; % initial value
    for n = 2:N
        r = rand; 
        xi = -1+2*round(r); % P(xi=1)=P(xi=-1)=0.5
        X(n) = X(n-1) + a*X(n-1)*dt + b*sqrt(dt)*xi; % simple euler method for SDE
        t(n) = t(n-1) + dt;
    end
    XN(m) = X(N);
    Y(m) = normrnd(mu,sigma);
end

histfit(XN)

figure
histfit(Y)

figure
plot(t,X)
hold on
xlabel('t')
ylabel('X')
