Residual-Based Monitoring of Autocorrelated Data
In this section, we demonstrate how the ResidualStatistic
and AbstractPhase2
interfaces can be extended to accommodate custom data types. Specifically, we focus on monitoring the residuals of an autoregressive $\text{AR}(1)$ model, $ y_t = \phi y_{t-1} + \varepsilon_{t}, \quad \varepsilon_{t} \sim N(0,1), $ using an EWMA control chart.
First, the required packages are imported,
using StatisticalProcessMonitoring, Distributions, Random, NLopt, Plots, Parameters, LaTeXStrings
We then define a specialization of the ResidualStatistic
type, which will be used to apply the control chart to the $\text{AR}(1)$ residuals.
julia> import StatisticalProcessMonitoring.residual!,
StatisticalProcessMonitoring.new_data!
julia> mutable struct AR1Statistic{S} <: ResidualStatistic
stat::S
phi::Float64
ym1::Float64
end
julia> function residual!(x, S::AR1Statistic)
yhat = x - S.phi * S.ym1
S.ym1 = x
return yhat
end
To sample observations from an $\text{AR}(1)$ process, we define a new type called Phase2AR1
.
julia> @with_kw mutable struct Phase2AR1 <:
StatisticalProcessMonitoring.AbstractPhase2
phi::Float64
y::Float64 = 0.0
init::Bool = false
end
The initial observation $y_0$ required to sample $y_1$ is initialized using the stationary distribution $N(0, 1 / (1 + \phi^2 ))$.
julia> function new_data!(PH2::Phase2AR1)
if !PH2.init
PH2.y = randn() / sqrt(1 - PH2.phi^2)
PH2.init = true
end
yhat = PH2.phi * PH2.y + randn()
PH2.y = yhat
return yhat
end
We consider an $\text{AR}(1)$ model with $\phi = 0.5$.
julia> seed = 4398354798
julia> Random.seed!(seed)
julia> phi = 0.5
julia> PH2 = Phase2AR1(phi = phi)
We define an EWMA control chart applied to the residuals of the $\text{AR}(1)$ model using the AR1Statistic
object.
julia> STAT = AR1Statistic(EWMA(λ = 0.1), phi, 0.0)
We use a two-sided control limit for the control chart.
julia> LIM = TwoSidedFixedLimit(1.0)
We set the ARLIC to 500 and create the ControlChart
object.
julia> NOM = ARL(500)
julia> CH = ControlChart(STAT, LIM, NOM, PH2)
The smoothing constant of the EWMA control chart is optimized against an anticipated persistent mean shift of $\delta = 2$.
julia> delta = 2.0
julia> rlsim_oc = x -> run_sim_oc(x, shift = delta)
julia> settings = OptSettings(verbose = false, minpar = [0.001],
maxpar = [0.99])
julia> opt = optimize_design!(CH, rlsim_oc, settings, optimizer = :LN_BOBYQA)
julia> print(opt)
1-element Vector{Float64}:
0.13830568163003895
We consider 100 Phase II observations from the $\text{AR}(1)$ process with $\phi = 0.5$.
julia> n = 100
julia> tau = 50
julia> DGP = Phase2AR1(phi = phi)
julia> y = [new_data!(DGP) for _ in 1:n]
julia> y[(tau+1):n] .+= delta
The control chart is then applied to the data and the results are plotted.
julia> proc = apply_chart(CH, y)
julia> plt = plot_series(proc, dpi = 300, label = "", xlab = L"t",
ylab = L"C_t")
julia> vline!(plt, [tau], label = "tau", linestyle = :dot, colour = "black")