%Assignemnt 1 Question 2 Part 2

%Part 2

clc; close all;

k = (2/3);
h3 = 0.01;
mystoeul(1) = 10;
sig = 0.2;
timesteps3 = 5/h3;

for t = [2:1:timesteps3+1];
    mystoeul(t) = mystoeul(t-1) - h3*k*mystoeul(t-1) + sig*sqrt(h3)*randn;
end

timevect3 = [0:h3:5];
figure;
plot(timevect3,mystoeul)
title('Figure 7: h = 0.01 Stochastic Implementation of Euler')

%Part 2b

storevals = zeros(15,5/h3);
mystoeul1(1) = 10;
storevals(:,1) = mystoeul1(1);

figure;
for n = 1:1:15
    for t = [2:1:timesteps3+1];
        mystoeul1(t) = mystoeul1(t-1) - h3*k*mystoeul1(t-1) + sig*sqrt(h3)*randn;
        storevals(n,t) = mystoeul1(t);
    end
end

for n = 1:1:15
    hold on;
    plot(timevect3,storevals(n,:))
end
title('Figure 8: h = 0.01 15 Repetitions of SDE with Average')

for i = 1:1:timesteps3+1
    sumvals(i) = sum(storevals(:,i)); %sums values
end

avgvals = sumvals./15; %calculates the average values
plot(timevect3,avgvals,'red','LineWidth',2)
hold off;

%The average of an infinite number of stochastic runs will result in the
%deterministic solution as calculated in part a.

%Part 2c - for time values greater than 5

k = (2/3);
step = 0.01;
mysto(1) = 0;
sig = 0.2;
timeval = 10;
tstep = timeval/step;

for t = [2:1:tstep+1];
    mysto(t) = mysto(t-1) - step*k*mysto(t-1) + sig*sqrt(step)*randn;
end

tvect = [0:step:timeval];
figure;
subplot(1,2,1), plot(tvect,mysto)
title('Figure 9: sig=0.2 w/ x0=0, Noise')
subplot(1,2,2), hist(mysto,40)
title('Figure 10: Histogram of sig = 0.2')

mean02 = mean(mysto)
std02 = std(mysto)

sig2 = 5;
mysto1(1) = 0;
for t = [2:1:tstep+1];
    mysto1(t) = mysto1(t-1) - step*k*mysto1(t-1) + sig2*sqrt(step)*randn;
end

tvect = [0:step:timeval];
figure;
subplot(1,2,1), plot(tvect,mysto1)
title('Figure 11: sig=5 w/ x0=0, Noise')
subplot(1,2,2), hist(mysto1,40)
title('Figure 12: Histogram of sig = 5')

mean5 = mean(mysto1)
std5 = std(mysto1)

%The width of the two histograms are dependent upon the value of sigma
%(standard deviation) that is plugged in.  On average, the width (std) of
%the histogram will equal to the standard deviation that is used to
%calculate the gaussian term.  The mean is always around zero, but has some
%error due to the stochastic nature of the function.