function [S, t] = geom_wiener_process(n, S0, alpha, sigma, T, graph, n_instances)
% Simulation of the geometrical Wiener process for option pricing.
% 
%[S, t]
%   S   Vector with length of 'n+1' representing one scenario of the geometric 
%       Wiener process
%   t   Time series
%
%GEOM_WIENER_PROCESS(N,S0,ALPHA,SIGMA,T,GRAPH,N_INSTANCES)
%   N   Number of subintervals of [0,T] (granularity)
%       Default value: 1000
%   S0  Initial price of the option
%       Defualt value: 20
%   ALPHA Expected rate of return of time period 'T'
%       Default value: 0.1
%   SIGMA Volatility in the option price
%       Default value: 0.3
%   T   Time interval [0,T]
%       Default value: 1
%   GRAPH   'true' or 'false' whether to plot the result or not
%           Default value: false
%   N_INSTANCES Number of instances of the process to be plotted
%               Default value: 1
%               Note: If the number of instances to plot is bigger than 2,
%               then the returned vector W is highlighted with RED colored
%               DASH-DOT line

%% Input parameters handling
if nargin < 7; n_instances = 1; end
if nargin < 6; graph = false; end
if nargin < 5; T = 1; end
if nargin < 4; sigma = 0.3; end
if nargin < 3; alpha = 0.1; end
if nargin < 2; S0 = 100; end
if nargin < 1; n = 1000; end

%% Computation task
S = zeros(n+1, n_instances);
S(1,:) = S0;
dt = T/n;  % Time step
t = (0:dt:T);

% Simulation via dS
for i = 1:n
    dS = alpha*S(i,:)*dt + sigma*S(i,:)*sqrt(dt).*randn(1,n_instances);
    S(i+1,:) = S(i,:) + dS;
end

% Simulation via Wiener process
% t = t(1:end-1);
% W = wiener_process(n,T,false,n_instances);
% S = S0*exp((alpha-0.5*sigma^2)*repmat(t',1,n_instances) + sigma*W);


%% Plotting
% TODO: can be parametrized for both geom. and std. wiener process
if graph
    figure
    % Statistics
    maxW = max(S, [], 1); minW = min(S, [], 1);
    % Keep the indices for plotting
    % To plot the run with maximum and minimum value
    [maxW, imax] = max(maxW); [minW, imin] = min(minW);
    lastS = S(end,:);
    % Plot itself
    % Plot only reasonable number of instances (1000.)
    % If the max and min runs are in those first 1000, then don't plot
    % them again.
    n_to_print = min(1e3, n_instances);
    if imax <= n_to_print; imax = []; end
    if imin <= n_to_print; imin = []; end
    plot(t,S(:,[1:n_to_print imax imin])); hold on; 
    % Plot editing
    a = gca;  % Get the current axis
    a.Position([1 2 4]) = [.11 .3 .6];  % [1, 2] position, [3, 4] w, h
    title("Simulation of the geometric Wiener process")
    xlabel("\it t")
    ylabel("\it S_t")
    grid on;
    message = "Minimal value of S: " + minW + newline + ...
              "Maximal value of S: " + maxW + newline + ...
              "Mean of values of S at time T: " + mean(lastS) + newline + ...
              "St. dev. of values of S at time T: " + std(lastS);
    annotation('textbox',[.1 .19 0 0],'String',message,'FitBoxToText','on')
end