clear;

%% Data

thetas = [0, 1/2, 1];
X0 = 1;
T= 1;
M = 1e4;
epsilons = 1./[20, 40, 60];
Ns = 2.^[2, 4, 6];

%% Stability

epsilon = 1/20;
sigma = sqrt(2/epsilon);
b = 1/epsilon;

for i = 1:length(thetas)
    theta = thetas(i);
    for j = 1:length(Ns)
        N = Ns(j);
        h = 1/N;
        t = linspace(0, T, N+1);
        C1 =  1/(1 + b*h*theta);
        C2 = (1 - b*(1-theta)*h);
        C3 =  sigma*sqrt(h);
        X = nan(1, N+1);
        X(1) = X0;
        for k = 1:N
            X(k+1) = C1 * (C2 * X(k) + C3 * randn);
        end
        figure;
        plot(t, X, 'Linewidth', 3);
        title(['$\theta = $', num2str(theta), ', $h = $', num2str(h)], 'interpreter', 'latex');
        set(gca, 'FontSize', 20);
        saveas(gcf, ['figures/exercise6_3_', num2str(i), num2str(j), '.png']);
    end
end

%% Approximation density function

N = 2^6;
h = 1/N;

for i = 1:length(thetas)
    theta = thetas(i);
    for j = 1:length(epsilons)
        epsilon = epsilons(j);
        sigma = sqrt(2/epsilon);
        b = 1/epsilon;
        C1 =  1/(1 + b*h*theta);
        C2 = (1 - b*(1-theta)*h);
        C3 =  sigma*sqrt(h);
        X = X0 * ones(M, 1);
        for k = 1:N-1
            X = C1 * (C2 * X + C3 * randn(M, 1));
        end
        figure;
        histogram(X, 'Normalization','pdf');
        hold on;
        x_pdf = linspace(min(X), max(X), 1000);
        plot(x_pdf, normpdf(x_pdf, 0, sqrt(sigma^2/(2*b))), 'Linewidth', 3);
        title(['$\theta = $', num2str(theta), ', $\epsilon = $', num2str(epsilon)], 'interpreter', 'latex');
        set(gca, 'FontSize', 20);
        saveas(gcf, ['figures/exercise6_4_', num2str(i), num2str(j), '.png']);
    end
end