% The Exponentially-Weighted Control Charts for Geometrically Moving Standard Deviations
%% Authors: Saeed Maghsoodloo and Daniel F. Silva

% Inputs
format long;
w=0.25; w1=1-w; w2=w1^2; w3=w1^3;

sig0=1;
L1=3;
L2=3;

k=1.1;
tmax = 150;

nvec=[4,5,6,7,8,9,10,11,12,15,20];
nnum=length(nvec);

Vvec = zeros(tmax,nnum);
SEvec = zeros(tmax,nnum);
LCLvec = zeros(tmax,nnum);
UCLvec = zeros(tmax,nnum);
alpha3vec = zeros(tmax,nnum);
ARLvec = zeros(tmax,nnum);

%% Section 2.1
c4=@(n)sqrt((n-1)/2)*gamma(n/2)/gamma((n+1)/2);
VG=@(n,t)w*(sig0.^2)*(1-c4(n).^2)*(1-w1.^(2*t))/(2-w);
SEG=@(n,t)sqrt(VG(n,t));
LCL=@(n,t)sig0-L1*SEG(n,t);
UCL=@(n,t)sig0+L2*SEG(n,t);
Rho=@(t)w1*sqrt((1-w2^(t-1))/(1-w2^t));

%% Section 3
g=@(n,t)sig0*sqrt(w*(1-c4(n).^2)*(1-w1.^(2*t))/(2-w));
ZL=@(n,t)(1-L1*g(n,t)-k*(c4(n)*(1-w1.^t)+w1.^t))/(k*g(n,t));
ZU=@(n,t)(1+L2*g(n,t)-k*(c4(n)*(1-w1.^t)+w1.^t))/(k*g(n,t));
Beta=@(n,t)normcdf(ZU(n,t))-normcdf(ZL(n,t));
ARL=@(n,t)1/(1-Beta(n,t));

%% Check integral
fun=@(y,n,t)1/(k*g(n,t)*sqrt(2*pi)).*exp(-0.5*((y-k*(c4(n)*(1-w1.^t)+w1.^t))/(k*g(n,t))).^2);
Int=@(n,t)integral(@(y) fun(y,n,t),1-L1*g(n,t),1+L2*UCL(n,t));

%% Appendix
alpha3=@(n,t)w.^(1.5)*c4(n)*((3-2*n)/(n-1)+2*c4(n).^2)*(1-w3.^t)/((1-w1.^3)*((1-c4(n).^2)*(1-w1.^(2*t))/(1+w1)).^1.5);

%% Main code
for j = 1:nnum
for t = 1:tmax
Vvec(t,j)=VG(nvec(j),t);
SEvec(t,j)=SEG(nvec(j),t);
LCLvec(t,j)=LCL(nvec(j),t);
UCLvec(t,j)=UCL(nvec(j),t);
alpha3vec(t,j)=alpha3(nvec(j),t);
ARLvec(t,j)=ARL(nvec(j),t);
end
end