% Illustration av autokorrelerade data
% 20061126
close all hidden

%% Utgår från exempel med tanken kapitel 9.4

V=3;            %(m^3)
f=0.03;         %(m^3/s)
delta_t=25;      % Samplingstid
T=V/f
e=exp(1);
a=1-e^(-delta_t/T);
rho=1-a;

x=66;
N=1000;

for k=2:N
    x(k)=(1-a)*x(k-1)+a*normrnd(66,1,1,1);
end
figure
plot(x,'.-b')
ylabel('Utkoncentrationen x_t','FontSize',12)
xlabel('Prov nummer','FontSize',12)
title('Blandningstank. Korrelation \rho = 0.25','FontSize',18)
%%

figure

plot(x(1:N-1),x(2:N),'o')
axis equal
xlabel('x_{t-1}','FontSize',12)
ylabel('x_{t}','FontSize',12)
title('Plot av x_t mot x_{t-1}, \rho=0.25','FontSize',18)

figure
r=[];
korrtid=25;
for k=0:korrtid;
    KORR=corrcoef(x(1:N-k),x(1+k:N));
    r=[r;KORR(1,2) ];
end
k=0:korrtid;
stem(k,r,'Filled')
set(gca,'Ylim',[-1,1])

xlabel('Förskjutning k','FontSize',12)
ylabel('Autokorrelationen r(k)','FontSize',12)
title('Autokorrelationsfunktionen','FontSize',18)

%% Moving range
MRM= mean(abs(diff(x)));

UCL=mean(x)+3*MRM/1.128;
LCL=mean(x)-3*MRM/1.128;

figure
subplot(2,1,1)
plot(x,'.-')
hold
plot(x*0+UCL,'--r','LineWidth',2)
plot(0*x+LCL,'--r','LineWidth',2)

ylabel('x','FontSize',12)
title('MR-diagram för korrelerade data, \rho=0.25','FontSize',18)

subplot(2,1,2)
MR=abs(diff(x));
plot(MR,'.-')
hold

UCL=3.267*MRM;
LCL=0;
plot(x*0+UCL,'--r','LineWidth',2)
plot(0*x+LCL,'--r','LineWidth',2)

ylabel('MR','FontSize',12)

%% EWMA
L=0.1;  % Parametern Lambda

Mu0=mean(x);
z=Mu0;  % Initiera

for i=1: length(x)
    UCL(i)=Mu0+ 2.7*sqrt(L/(2-L)*(1-(1-L).^(2*i)));
    LCL(i)=Mu0- 2.7*sqrt(L/(2-L)*(1-(1-L).^(2*i)));
    CL(i)=Mu0;
end
z(1)=[];

figure
hold

plot(x,'o-k','MarkerFaceColor','b')

k=0.5:1:length(x);
stairs(k,UCL,'r','LineWidth',2)
stairs(k,LCL,'r','LineWidth',2)
plot(CL,'r--')

ylabel('x','FontSize',12)
title('EWMA för korrelerade data, \rho=0.25','FontSize',18)



%% Tag bort värden
% Jag mäter var K:te mätpunkt. K= tiden tills att autokorrelationen är
% noll.

K=10;   % 

x_arkiv=x;

x=x(1:K:N);

%% Moving range

MRM= mean(abs(diff(x)));

UCL=mean(x)+3*MRM/1.128;
LCL=mean(x)-3*MRM/1.128;

figure
subplot(2,1,1)
plot(x,'.-')
hold
plot(x*0+UCL,'--r','LineWidth',2)
plot(0*x+LCL,'--r','LineWidth',2)

ylabel('x','FontSize',12)
title('MR-diagram för var 10:e prov, \rho=0.25','FontSize',18)

subplot(2,1,2)
MR=abs(diff(x));
plot(MR,'.-')
hold

UCL=3.267*MRM;
LCL=0;
plot(x*0+UCL,'--r','LineWidth',2)
plot(0*x+LCL,'--r','LineWidth',2)

ylabel('MR','FontSize',12)

%% EWMA
L=0.1;  % Parametern Lambda

Mu0=mean(x);
z=Mu0;  % Initiera
UCL=[];
LCL=[];
CL=[];

for i=1: length(x)
    UCL(i)=Mu0+ 2.7*sqrt(L/(2-L)*(1-(1-L).^(2*i)));
    LCL(i)=Mu0- 2.7*sqrt(L/(2-L)*(1-(1-L).^(2*i)));
    CL(i)=Mu0;
end
z(1)=[];

figure
hold

plot(x,'o-k','MarkerFaceColor','b')

k=0.5:1:length(x);
stairs(k,UCL,'r','LineWidth',2)
stairs(k,LCL,'r','LineWidth',2)
plot(CL,'r--')

ylabel('x','FontSize',12)
title('EWMA för var 10:prov, \rho=0.25','FontSize',18)

%%

figure
N=length(x);
plot(x(1:N-1),x(2:N),'o')
axis equal
xlabel('x_{t-1}','FontSize',12)
ylabel('x_{t}','FontSize',12)
title('Plot av x_t mot x_{t-1}, \rho=0.25','FontSize',18)
%% Autokorrelationen 
figure
r=[];
korrtid=25;
for k=0:korrtid;
    KORR=corrcoef(x(1:N-k),x(1+k:N));
    r=[r;KORR(1,2) ];
end
k=0:korrtid;
stem(k,r,'Filled')
set(gca,'Ylim',[-1,1])

xlabel('Förskjutning k','FontSize',12)
ylabel('Autokorrelationen r(k)','FontSize',12)
title('Autokorrelationsfunktionen','FontSize',18)