% Exercise sheet 5
% X_t = X_0 + B_t, X_0 ~ mu_0
% mu_0 uniform distribution on [-1,1]
% estimate P(X_T in A) for A=[2,3] and T=1

%% Monte Carlo
clear
close all

T = 1;
A = [2,3];
M = 10000; % number of realisations
X_T = zeros(1,M);
B_T = sqrt(T)*randn(1,M);
c = zeros(1,M);
x0 = zeros(1,M);

for m = 1:M
    r = rand;
    x0(m) = -1+2*r; % initial value 
    %B_T = normrnd(0,T);
    X_T(m) = x0(m) + B_T(m);
    if A(1) <= X_T(m) && X_T(m) <= A(2)
        c(m) = 1;
    else
        c(m) = 0;
    end
end

hist(x0)
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w')
%xlabel('x_0')
%ylabel('')
title('x_0')
figure
hist(X_T)
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w')
%xlabel('x_0')
%ylabel('')
title('X_T')

P = sum(c)/M % approx 0.04

%% Kolmogorov forward equation

% discretisation of the state space
clear L X

T = 1;
A = [2,3];
dx = 0.1;
dt = 0.005;
Time = 0:dt:T;
X = -5:dx:5;
s = length(X);
t = length(Time);
rho = zeros(t, s);

% create discretization of the generator (self adjoint)

for i = 2:s-1
    L(i,i) = -2;
    L(i,i+1) = 1;
    L(i,i-1) = 1;
end
L(1,1) = -1;
L(1,2) = 1;
L(s,s) = -1;
L(s,s-1) = 1;

L = 1/(2*dx^2)*L; 

mu_0 = -1:dx:1;
s_0 = length(mu_0);

for i = 1:s
    if mu_0(1) <= X(i) && X(i) <= mu_0(s_0)
        rho(1,i) = 1/(s_0*dx);
    else
        rho(1,i) = 0;
    end
end

for i=1:(length(Time)-1)
    rho(i+1,1:s) = rho(i,1:s)+dt*rho(i,1:s)*L;
end;
rho_T = rho(t,1:s);% distribution at time T
figure
plot(X,rho_T)
%h = findobj(gca,'Type','patch');
%set(h,'FaceColor','b','EdgeColor','w')
%xlabel('x_0')
%ylabel('')
title('rho_T, p=')

p = 0;

for i = 1:s
    if A(1) <= X(i) && X(i) <= A(2)
        p = p + dx*rho_T(i);
    else
    end
end

[Xg, Tg] = meshgrid(X, Time);
figure;
mesh(Xg, Tg, rho);
p % für dx=0.1: 0.0435, für dx=0.01: 0.0378 
        
