clear;

%% Data

rng(2022);

T = 1;
nSteps_vec = 2.^[4, 6, 8];
nTrajectories_vec = [20, 200, 2000];

lambda = 2;
mu = 1;
X0 = 1;

%% Brownian motion

for i = 1:length(nSteps_vec)
    nSteps = nSteps_vec(i);
    dt = T/nSteps;
    tVec = 0:dt:T;
    for j = 1:length(nTrajectories_vec)       
        nTrajectories = nTrajectories_vec(j); 
        W = BM(T, nSteps, nTrajectories);
        figure;
        plot(tVec, W, 'color', 0.75*ones(1, 3));
        xlabel('t', 'Fontsize', 20);
        title(['{\Delta}t = ', num2str(dt), ', ', num2str(nTrajectories), ' sample paths'], 'Fontsize', 20);
        hold on;
        plot(tVec, mean(W), 'r', 'Linewidth', 2);
        plot(tVec, zeros(size(tVec)), 'g', 'Linewidth', 2);
        set(gca, 'FontSize', 20);
        saveas(gcf, ['figures/exercise1W_', num2str(i), num2str(j), '.png']);
    end
end
   
%% Exponential stochastic process

for i = 1:length(nSteps_vec)
    nSteps = nSteps_vec(i);
    dt = T/nSteps;
    tVec = 0:dt:T;
    for j = 1:length(nTrajectories_vec)       
        nTrajectories = nTrajectories_vec(j); 
        W = BM(T, nSteps, nTrajectories);
        X = X0*exp((lambda - (mu^2)/2)*tVec + mu*W);
        figure;
        plot(tVec, X, 'color', 0.75*ones(1, 3));
        xlabel('t', 'Fontsize', 20);
        title(['{\Delta}t = ', num2str(dt), ', ', num2str(nTrajectories), ' sample paths'], 'Fontsize', 20);
        hold on;
        plot(tVec, mean(X), 'r', 'Linewidth', 2);
        plot(tVec, X0*exp(lambda*tVec), 'g', 'Linewidth', 2);
        set(gca, 'FontSize', 20);
        %saveas(gcf, ['figures/exercise1Exp_', num2str(i), num2str(j), '.png']);
    end
end