%% Asset Model naive attempt

RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', 42));

T = 1;
n = 1000;
ts = linspace(0,T,n+1);

S0=100;
mu = 0.1;
sigma = 0.15;

Y = randn(1,n);
Y = [0, Y];

S = S0 * exp( (mu-sigma^2/2) * ts + sigma * sqrt(ts) .* Y);

figure;
plot(ts, S);
xlabel('t');
ylabel('S_t');

%% Brownian motion via partial sums: one path

RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', 42));

T = 10;
S0=100;
mu = 0.08;
sigma = 0.1;

n = 1000;
h = T/n;
ts = linspace(0,T,n+1);
Bt = [0, cumsum(randn(1,n) * sqrt(h))];

figure;
subplot(2,1,1);
plot(ts, Bt);
ylabel('W_t');

subplot(2,1,2);
plot(ts, S0 * exp( (mu-sigma^2/2) * ts + sigma * Bt));
xlabel('t');
ylabel('S_t');

%% Simulation of Brownian motion via Paley-Wiener representation

M = 1000;
n = 10000;
RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', 42));

% Fourier coefficients
coef0 = randn(1,1) / sqrt(2*pi);
coef = randn(1, M) * 2 / sqrt(pi) ./ (1:M);

ts = linspace(0, 2*pi, n);

figure;
hold on;
cc = {'r', 'b', 'k'};

Mt = [10, 50, 1000];
thick = [3,2,1];

for j = 1:length(Mt)
    mts = (1:Mt(j))' * ts/2;
    y = coef0 * ts + coef(1:Mt(j)) * sin(mts);
    
    plot(ts, y,cc{j}, 'linewidth',thick(j));
end
legend(sprintf('%d coefficients', Mt(1)),...
    sprintf('%d coefficients', Mt(2)),...
    sprintf('%d coefficients', Mt(3)));
xlim([0, 2*pi]);
xlabel('t');
ylabel('W_t');

%% Timescale invariance: same path sampled at different frequencies

RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', 42));

figure;
n = 10000;
T = 1;
h = T/n;
ts = linspace(0,T,n+1);

Bt = [0, cumsum(randn(1,n) * sqrt(h))];

for j = 1:3
    subplot(3, 1, j);
    ev = 1+(0:99)*10^(3-j);
    plot(ts(ev), Bt(ev));
    axis([0, 10^(1-j), -sqrt(10^(1-j))*2, sqrt(10^(1-j))]);
end

%% Brownian motion: many paths
figure;
RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', 42));

T = 1;
n = 10000;
h = T/n;
ts = linspace(0,T,n+1);

paths = 50;

Bt = [zeros(paths,1), cumsum(randn(paths,n) * sqrt(h), 2)];

% Plot the paths
subplot(1,3,1);
plot(ts, Bt, 'b');

% Sample mean and std
subplot(1,3,2);
plot(ts, mean(Bt), 'k', ts, std(Bt), 'r');
hold on;
plot(ts, zeros(1,n+1), 'k--', ts, sqrt(ts), 'r--');

% Sum of difference-squares
subplot(1,3,3);
plot(ts(2:end), cumsum(diff(Bt,1,2).^2, 2), 'k');


%% Geometric Brownian motion
RandStream.setGlobalStream(RandStream('mt19937ar', 'Seed', 42));

T = 20;
n = 5000;
h = T/n;
ts = linspace(0,T,n+1);

Bt = [0, cumsum(randn(1,n) * sqrt(h))];

figure;

mus = [0.1, 0.1, 0.1, 0.05, 0.05, 0.05];
sigmas = [0.1, 0.2, 0.4, 0.1, 0.2, 0.4];

for j = 1:6
    subplot(2,3,j);
    mu = mus(j);
    sigma = sigmas(j);

    Xt = exp((mu-sigma^2/2)*ts + sigma*Bt);
    plot(ts, Xt);
    ylim([0,6]);
    title(['Drift ', num2str(mu), '  Vol ', num2str(sigma)]);
end
