% Towards a Polyalgorithm for Land Use Change Detection,
% R. Saxena, L. T. Watson, R. H. Wynne, E. B. Brooks, V. A. Thomas,
% Y. Zhiqiang; Computers & Geosciences, 12/2017.
%
\magnification=1000
\hoffset=-.1875truein \voffset=0.0truein 
\parindent=1pc \parskip=0pt
\baselineskip=12pt plus 2pt minus 1pt \lineskiplimit=2pt
\lineskip=2pt plus 2pt \tolerance=1000 %\nopagenumbers
\abovedisplayskip=4pt plus 4pt \belowdisplayskip=4pt plus 4pt

%fonts
\font\bfXII=cmbx12 \font\caps=cmcsc10 

% allows spaces to stretch
%\spaceskip=3.3333pt plus 3pt minus 1.1111pt \hyphenpenalty=20

%\input epsf % For Postscript figures: {\epsfxsize=  \epsffile{}}

% TeX editing macros:  to underline, double-underline, wavy underline,
% horizontal strike out, slanted strike out "T", use
% \uline{T}, \uuline{T}, \uwave{T}, \sout{T}, \xout{T}, respectively.
%\input ulem.sty \input colordvi
%\def\BC#1#2{\sout{#1}{\textPeach #2\textBlack}} %TeX source markup.
%to maKE a change in this, enter
%to \BC{maKE}{make} a change in this

% Color definitions.
%\input colordvi
%\def\BEGINC{\textPeach} \def\ENDC{\textBlack}
%\let\BEGINC=\relax \let\ENDC=\relax

% other definitions.
%
\def\lheading#1#2{\par\vskip-\lastskip\medbreak % > 1 line heading 
  {\bfXII\setbox11=\hbox{#1 }\hangindent=\wd11 \hangafter=1 
  \noindent #1 #2\par}\nobreak\smallskip} 

% For automatic section numbering:
\newcount\secNUM \secNUM=0 \newcount\subsecNUM \subsecNUM=0
\newcount\subsubsecNUM \subsubsecNUM=0
\def\cursection{\number\secNUM.}
\def\cursubsection{\number\secNUM.\number\subsecNUM.}
\def\cursubsubsection{\number\secNUM.\number\subsecNUM.\number\subsubsecNUM.}
%
% produce a numbered (sub)(sub)section heading: \[sub][sub]section{name}
%
\def\lsection#1{\advance\secNUM by 1\subsecNUM=0 % > 1 line heading
 \lheading{\cursection}{#1}}
% end of numbering macros.

\def\sgn{{\rm sgn}\,}
\def\:{\mathrel{:\mathord=}}

\lsection{Exponentially Weighted Moving Average Change Detection (EWMACD)}
EWMACD, proposed in Brooks et al. (2014), is a kernel regression approach
modeling the time series as a linear combination of
trigonometric polynomials. The model is trained over data collected in the
initial two (or more) years. When the observations in the subsequent years
deviate from the values forecast by the model for a `substantial' length
of time (persistence), a change is declared (recovery or disturbance). The
training period as well as the persistence are parameters
of the algorithm. A positive flag is raised in instances of sustained
growth while sustained losses are indicated by negative flags.
EWMACD is able to capture seasonal changes, but is also
very sensitive to several algorithm parameters.

\bigskip
{\baselineskip=12pt plus 2pt minus 1pt \lineskiplimit=2pt
\leftline{\bf Algorithm EWMACD.}
\medskip

{\parindent=0pt

\leftskip 0pt {\bf for} row $r\:1$ {\bf step} 1 {\bf until} $R$

\leftskip 10pt {\bf for} column $c\:1$ {\bf step} 1 {\bf until} $C$ {\bf do}

\leftskip 20pt

{\bf Step 1:}
Write the time series data in the column $(D_{rc})_{\cdot b}$ as
$$(D_{rc})_{\cdot b} = \pmatrix{u \cr v\cr},$$
where the $M$-dimensional vector $u$ is deemed training
data and the $(S - M)$-dimensional vector $v$ as the test data. Let
$$\leftskip 0pt
X = \left[\matrix{ 1 & \sin t_1 & \cos t_1 & \cdots & \sin Kt_1 &
\cos Kt_1 \cr
 \vdots &\vdots &\vdots &\vdots &\vdots &\vdots \cr
 1 & \sin t_M & \cos t_M & \cdots & \sin Kt_M & \cos Kt_M \cr}\right]$$
be the Gram matrix for the time points $t_1$, $\ldots$,
$t_M$, using $K$ harmonics, where $M > 2K + 1$.  The least squares fit to
the training data $u$ is then written as
$$u(t) = \alpha_0 + \sum_{i=1}^K(\alpha_{2i-1} \sin it + \alpha_{2i}
\cos it)$$
with coefficients
$$\alpha = \left(X^tX\right)^{-1} X^tu$$
and residual
$$E(\alpha) = u - X\alpha.$$

\smallskip
{\caps Remark 1}. In practice $\alpha$ is computed via a QR factorization
of $X$, not by computing $(X^tX)^{-1}$ explicitly.
\smallskip

Next let
$I = \{ i \mid |E(\alpha)_i| < \gamma_1\},$
where $\gamma_1$ is a user defined threshold and $|I| > 2K + 1$. Calculate
the coefficients for an improved fit to the underlying signal as
$$\alpha^* = \bigl((X_{I\cdot})^tX_{I\cdot}\bigr)^{-1} (X_{I\cdot})^t
\,u_I.$$
With the refined coefficients $\alpha^*$, calculate the residuals for

{\parindent=20pt
\item{(i)}the complete time series $(D_{rc})_{\cdot b}$ as
$$E^*(\alpha^*) = (D_{rc})_{\cdot b} - \bar X \alpha^*,$$
where $\bar X_{s\cdot} = (1$, $\sin t_s$, $\cos t_s$, $\ldots$,
$\sin Kt_s$, $\cos Kt_s)$, for $s = 1$, $\ldots$, $S$.

\item{(ii)}the outlier-free time series as
$\bigl(E^*(\alpha^*)\bigr)_{\bar I}$, where
$\bar I = \{ s \mid |E^*(\alpha^*)_s| < \gamma_2\}$, $\gamma_2$ is a user
defined threshold, and

\item{(iii)} the outlier-free training set $\hat I
= \bar I \cap \{1,\ldots,M\}$ as
$$\bigl(E^*(\alpha^*)\bigr)_{\hat I}
= u_{\hat I} - X_{\hat I\cdot} \alpha^*,$$
where $|\hat I| > 2K + 1$.

}

\smallskip
{\caps Remark 2}. In the present implementation,
$$ \gamma_2 =  \cases{1.5 \eta, & $i \in [1, M]$,\cr
20 \eta, & $i \in (M,  S]$,\cr}$$
where $\eta$ is the standard deviation of the
first $M$ elements of the residual vector $E^*(\alpha^*)$.
\smallskip

{\bf Step 2:} Define the control limit vector $\tau$ by
$$\tau_i = \mu + \sigma L \sqrt{{\lambda \over 2-\lambda}
\bigl(1 - (1 - \lambda)^{2i}\bigr)},$$
for $i=1$, $2$, $\ldots$, $|\bar I|$,
where $\mu=0$ is used here, $\sigma$ is the standard deviation of the
outlier-free training data errors $\bigl(E^*(\alpha^*)\bigr)_{\hat I}$,
$L$ is the multiple of this standard deviation $\sigma$, and $\lambda\in
(0,1]$ is the weight given to the most recent residual in the exponentially
weighted moving average (EWMA) defined next.  $L$ is typically set to $3$
or slightly smaller depending on the value of $\lambda$.

{\bf Step 3:} Let $\bar I = \{j_1$, $j_2$, $\ldots$, $j_{|\bar I|}\}$,
$j_1 < j_2< \cdots < j_{|\bar I|}$. Define the vector $z$ by
$$\eqalign{z_1 &= \bigl(E^*(\alpha^*)\bigr)_{j_1},\cr
z_i &= (1-\lambda)z_{i-1} + \lambda\bigl(E^*(\alpha^*)\bigr)_{j_i}, \quad
i = 2, \ldots, |\bar I|.\cr}$$
This is the exponentially weighted moving average (EWMA) of the residual
$\bigl(E^*(\alpha^*)\bigr)_{\bar I}$.

{\bf Step 4:} Define the flag history $S$-vector $f$ by
$$f_s = \cases{\sgn(z_i) \big\lfloor|z_i/\tau_i|\big\rfloor,
&$s = j_i \in \bar I$,\cr
0, &otherwise.\cr}$$

If there is a run of $+1$ or $-1$ in the values $\sgn(\Delta f_s) =
\sgn(f_{s+1} - f_s)$ of length $\varpi$, called the `persistence', signal a
change at the index $s$ beginning the (nonzero) run.

\medskip
{\caps Remark 3}. Missing data is automatically handled by not assuming
that the time points $t_i$ are equally spaced. Alternatively, missing data
for time point $t_k$ can be handled by including $t_k$ in the sequence
$(t_1$, $t_2$, $\ldots$, $t_S)$, but excluding $t_k$ from the training
sequence $(t_1$, $t_2$, $\ldots$, $t_M)$ and $k$ from the sets $I$, $\bar
I$, and $\hat I$, which is equivalent to treating $(D_{rc})_{kb}$ as an
outlier and to setting the flag $f_k=0$.

\leftskip 10pt {\bf end}

\leftskip 0pt {\bf end}

}
}

\bigskip
The complexity of this approach is ${\cal O}\bigl((2K+1)^2S\bigr)$.
EWMACD outcomes are expected to be zeros when the time series trajectory
matches the trajectory in the training period. Nonzeros are expected
to appear when the trajectory deviates from the training period for a
substantial length of time (determined by control chart parameters and
the persistence). The first timepoint when the outcome
trajectory deviates from the prior constant (or zero) value is the 
estimated date of disturbance, the subsequent nonzero sequence of values
provides recovery (or, loss) period information. A sharp loss (e.g., harvest)
is expected to manifest as a sudden drop in the outcome as well. A gentler loss 
(e.g., mountain pine beetle) is expected to appear in the form of the
outcome trajectory gradually drifting away from zero (with an overall
negative slope). Similarly, recovery will appear in the form of
the trajectory following an overall positive slope, revealing
recovery period information. In Brooks et al.(2017), Edyn, an improved
version of EWMACD that retrains data after a disturbance is
sensed, is proposed (this development not implemented here).
In any case, since the training data in the postshock
period is unstable, the training period of Edyn after the first break is
expected to be unstable.
For all the experiments presented in this work, the first two years of data
is used for training EWMACD. The rest of the parameter values for EWMACD
are as follows: $K=2$, $L = 0.5$, $\lambda=0.3$,  $\varpi = 7$.


\vfill\eject\end
