%% A further note on the exponentially weighted moving average control
%% charts for monitoring gradual shifts in a process mean
% Authors: Saeed Maghsoodloo, Daniel F Silva
% Quality and Reliability Engineering International 38.8 (2022): 4182-4195

%% Input parameters
% w - omega: the weight of the present observation
% sig - sigma: the specified/target variance
% Tmax: number of periods of observations
% kvec: values of k, number of standard deviations in alternate hypothesis
% Lvec: values of L, parameter for width of the control limit
format longE
w=0.30;
w1=1-w;
sig=1;
Tmax = 100;
kvec = [0.00, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.5, 3, 3.25];
Lvec = [2.50 , 2.60 , 2.70 , 2.80 , 2.90 , 3.00 , 3.10 , 3.20];

knum = length(kvec);
Lnum = length(Lvec);

Vvec = zeros(1,Tmax);
SEvec = zeros(1,Tmax);
CovGvec = zeros(1,Tmax);
CorGvec = zeros(1,Tmax);
ARLmat= zeros(knum,Lnum);

%% Section 2 - Variance and Std. Error when sigma is specified
% VG: Variance of the G statistic at time t. Eq (4)
% SEG: Standard error of the G statistic at time t. Eqs (5, 6)
VG=@(t)(sig.^2)*w*(1-w1.^(2*t))/(1+w1);
SEG=@(t)sqrt(VG(t));

%% Section 3 - Covariance and correlation when sigma is specified
% CovG: Covariance of the G statistic between times t-1 and t. Eq (10)
% CorG: Corelation of the G statistic between times t-1 and t. Eqs (11, 12)
CovG=@(t) w1.*w.*sig.^2.*((1-w1.^(2.*t-2))/(1+w1));
CorG=@(t) w1.*sqrt((1-w1.^(2.*t-2))/(1-w1.^(2.*t)));

%% Section 4 - Control limits when sigma is unknown
% VG4: Variance of the G statistic at time t. Eq (15)
% SE4: Process standard error of the G statistic at time t. Eq (16)
% se44: Sample standard error of the G statistic at time t. Eq (17)
% LCT: Lower control limit at time t. Eq (18)
% UCT: Upper control limit at time t. Eq (18)
VG4=@(t)(sig^2*w)*(1-(1-w^(2*t)))/(1-w1)+sig^2*w1*(2-w1)/m
SE4=@(t) sig*sqrt(w(1-w1^(2*t))/(1+w1)+w1*(2-w1)/m)
se44=@(t) MR/d2*sqrt(w(1-w1^(2*t))/(1+w1)+w1*(2-w1)/m)
LCT=@(t) xbar -L*se44(t)
UCT=@(t) xbar +L*se44(t)

%% Section 6 - The average run-length when sigma is specified
% BetaLimit: Type II error probability for large t. Eq (25)
% g, fun : intermediate calculations to obtain BetaLimit
g=sqrt(w/(1+w1));
fun = @(Z1,Z2) 1/(2.*pi.*sqrt(1-w1.^2)).*exp(-(Z1.^2-2.*w1.*Z1.*Z2+Z2.^2)/(2.*(1-w1.^2))); %integrand
BetaLimit =@(k,L) integral2(fun,-L ,L ,-L-k/g, L-k/g )/0.9973; %eq.6.6.

%% Main code
% Calculates a table of average run lengths for specified sigma for
% multiple values of k and L provided by the user. (Table 1).
for k = 1:knum
for L = 1:Lnum
ARLmat(k,L)=1./(1-BetaLimit(kvec(k),Lvec(L)));
end
end

ARLmat;

for t = 1:Tmax
Vvec(1,t)=VG(t);
SEvec(1,t)=SEG(t);
CovGvec(1,t)=CovG(t);
CorGvec(1,t)=CorG(t);
end