%% sketch of the proposed approach with changing pdfs in time
clear all; close all; clc
%parameters
L0=1/4; %mean frequency
G0=10; %mean jump
tau=5; %decay timescale

N=1E4; %Number of events
lags = 10; %for autocorrelation function

Color=[0 0.45 0.74;
       1 0.41 0.16];
%% plot processes (first 100 events)
F1=figure(1); F1.WindowState='maximized'; pause(1)
Fsize=24;

dt = exprnd(1/L0,[1 N]); %intertime values
dX = exprnd(G0,[1,N+1]); %jump values
t = [0 cumsum(dt)]; %time values at event occurrence

%plot rainfall Poisson Process
subplot(2,1,1); box on; hold on
for i=1:N
    if i<=1E2 %plot only the first 100 events (graphical purposes)
        plot([t(i) t(i)],[0 dX(i)],'-b','LineWidth',2)
    end
end
set(gca,'FontSize',Fsize,'TickLabelInterpreter','latex','LineWidth',1,'Position',[0.13 0.61 0.775 0.365]);
xlabel('$t$','Interpreter','latex');
ylabel('$h$','Interpreter','latex');

%plot Flow rate Compound Poisson Process
subplot(2,1,2,'replace'); box on; hold on
X0 = dX(1);
Xvec=[];
Xvec2=[];
for i=2:N+1
    Xvec=[Xvec X0*exp(-(t(i)-t(i-1))/tau)]; %vector for histogram statistics
    if i<=1E2  %plot only the first 100 events (graphical purposes)
        plot([t(i-1) t(i-1)],[X0-dX(i-1) X0],'-b','LineWidth',2); 
        fplot(@(x) X0 * exp(-(x-t(i-1))/tau),[t(i-1) t(i)],'-b','LineWidth',2);
    elseif isempty(Xvec2)
        Xvec2=Xvec; %vector for histogram statistics for the first 10^2 data
    end
    X0 = Xvec(end) + dX(i);
end
set(gca,'FontSize',Fsize,'TickLabelInterpreter','latex','LineWidth',1,'Position',[0.13 0.11 0.775 0.365]);
xlabel('$t$','Interpreter','latex');
ylabel('$q$','Interpreter','latex');

%% plot statistics and autocorrelation function
F2=figure(2); F2.WindowState='maximized'; pause(1)

%Poisson Process histogram from data (N=10^2, N=10^4 and pdf)
subplot(2,2,1,'replace'); hold on; box on
%pdf from data with N=10^4
    histogram(dX,'NumBins',100,'Normalization','pdf','DisplayName','$p_h(h,N=10^4)$','FaceColor',Color(1,:));
%pdf from data with N=10^2
    histogram(dX(1:1E2),'Normalization','pdf','BinWidth',F2.Children(1).Children(1).BinWidth,'DisplayName','$p_h(h,N=10^2)$','FaceColor',Color(2,:));
%theoretical pdf 
    fplot(@(x) 1/G0 * exp(-x/G0),[0 max(dX)],'-r','LineWidth',2,'DisplayName','$p_h(h)=\frac{1}{\alpha} e^{-\frac{h}{\alpha}}$')

xlim([0 60])
set(gca,'FontSize',Fsize,'TickLabelInterpreter','latex','LineWidth',1);
xlabel('$h$','Interpreter','latex');
ylabel('$p_h(h)$','Interpreter','latex');
legend('Interpreter','latex');

%Poisson Process autocorrelation
subplot(2,2,2,'replace'); hold on; box on
% autocorrelation function for N=10^4 data
    [ac,Lvec]=xcorr(dX-mean(dX),lags,'normalized');
    stem(Lvec,ac,'Color',Color(1,:),'LineWidth',4,'MarkerFaceColor',Color(1,:),'DisplayName','$\rho(h, N=10^4)$');
% autocorrelation function for N=10^2 data
    [ac,Lvec]=xcorr(dX(1:1E2)-mean(dX(1:1E2)),lags,'normalized');
    stem(Lvec,ac,'Color',Color(2,:),'LineWidth',2,'MarkerFaceColor',Color(2,:),'DisplayName','$\rho(h, N=10^2)$');

axis([0 lags -0.2 1]);
xticks(0:1:lags);
yticks(-0.2:0.1:1);
set(gca,'FontSize',Fsize,'TickLabelInterpreter','latex','LineWidth',1);
legend('Interpreter','latex');

%Compound Poisson Process histogram from data (N=10^2, N=10^4 and pdf)
subplot(2,2,3,'replace'); hold on; box on
%pdf from data with N=10^4
    histogram(Xvec,'BinWidth',0.5,'Normalization','pdf','DisplayName','$p_q(q,N=10^4)$');
%pdf from data with N=10^2
    histogram(Xvec2,'Normalization','pdf','BinWidth',0.5,'DisplayName','$p_q(q,N=10^2)$');
%theoretical pdf 
    fplot(@(x) exp(-x/G0) * (x/G0)^(L0*tau) / x / gamma(L0*tau),[0 max(dX)],'-r','LineWidth',2,'DisplayName','$p_q(q)=\frac{e^{-\frac{q}{\alpha}} \, \left(\frac{q}{\alpha}\right)^{\lambda \tau}}{q \, \Gamma[\lambda \tau]}$')

xlim([0 60])
set(gca,'FontSize',Fsize,'TickLabelInterpreter','latex','LineWidth',1);
xlabel('$q$','Interpreter','latex');
ylabel('$p_q(q)$','Interpreter','latex');
legend('Interpreter','latex');

%Poisson Process autocorrelation
subplot(2,2,4,'replace'); hold on; box on
% autocorrelation function for N=10^4 data
    [ac,Lvec]=xcorr(Xvec-mean(Xvec),lags,'normalized');
    stem(Lvec,ac,'Color',Color(1,:),'LineWidth',4,'MarkerFaceColor',Color(1,:),'DisplayName','$\rho(h, N=10^4)$');
% autocorrelation function for N=10^2 data
    [ac,Lvec]=xcorr(Xvec(1:1E2)-mean(Xvec(1:1E2)),lags,'normalized');
    stem(Lvec,ac,'Color',Color(2,:),'LineWidth',2,'MarkerFaceColor',Color(2,:),'DisplayName','$\rho(h, N=10^2)$');

axis([0 lags -0.2 1]);
xticks(0:1:lags)
yticks(-0.2:0.1:1)
set(gca,'FontSize',Fsize,'TickLabelInterpreter','latex','LineWidth',1);
legend('Interpreter','latex');

%% export figures
% exportgraphics(F1,'PP_and_CPP.png','Resolution',600)
% exportgraphics(F2,'pdf_and_AutoCorr.png','Resolution',600)
