function [Tr] = ccpoisson_process(mu, sigma, N, lambda, T, ... 
                                 graph, n_instances)
% Simulation of the Compensated Compound Poisson process
%
%[Tr]
%   Tr  Trajectory of the process
%
%CCPOISSON_PROCESS(N,LAMBDA,T,GRAPH,N_INSTANCES)
%   MU, SIGMA   are parameters for normal distribution of the jump's size
%   N   granularity of the [0,T] interval
%   All the other input parameters correspond to 
%   the POISSON_PROCESS function.

%% Input parameters handling
if nargin < 7; n_instances = 1; end
if nargin < 6; graph = false; end
if nargin < 5; T = 10; end
if nargin < 4; lambda = 1; end
if nargin < 3; N = 1000; end
if nargin < 2; sigma = 1; end
if nargin < 1; mu = 0; end

%% Computation task
dt = T/N;
t = (0:dt:T);
% -lambda*mu*t part will be the same for all of the elements
Tr = repmat(-lambda*mu*t',1,n_instances);  

for i = 1:n_instances
    p = poisson_process(0, lambda, T, false, 1, false);
    jumps = mu + sigma*randn(length(p),1);  % jumps are from N(mu,sigma)
    jumps(1) = 0;
    for j = 2:length(p)
        % adds the jump size from t to the end
        Tr(t > p(j), i) = Tr(t > p(j), i) + jumps(j);
    end
end

st = length(t);
%% Plotting
if graph
    figure
    % Statistics
    maxTr = max(Tr, [], 1); minTr = min(Tr, [], 1);
    % Keep the indices for plotting
    % To plot the run with maximum and minimum value
    [maxTr, imax] = max(maxTr); [minTr, imin] = min(minTr);
    lastTr = Tr(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,Tr(1:st,[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 compensated compound Poisson process")
    xlabel("\it t")
    ylabel("\it C*_t")
    grid on;
    message = "Minimal value of C*: " + minTr + newline + ...
              "Maximal value of C*: " + maxTr + newline + ...
              "Mean C*(T): " + mean(lastTr) + newline + ...
              "St. dev. C*(T): " + std(lastTr);
    annotation('textbox',[.1 .19 0 0],'String',message,'FitBoxToText','on')
    message2 = "\lambda: " + lambda + newline + ...
               "\mu: " + mu + newline + ...
               "\sigma: " + sigma;    
    annotation('textbox',[.5 .19 0 0], ...
               'String',message2, ...
               'FitBoxToText','on')
end
