NAG CL Interface
g13nec (cp_binary_user)
1
Purpose
g13nec detects change points in a univariate time series, that is, the time points at which some feature of the data, for example the mean, changes. Change points are detected using binary segmentation for a usersupplied cost function.
2
Specification
The function may be called by the names: g13nec or nag_tsa_cp_binary_user.
3
Description
Let ${y}_{1:n}=\left\{{y}_{j}:j=1,2,\dots ,n\right\}$ denote a series of data and $\tau =\left\{{\tau}_{i}:i=1,2,\dots ,m\right\}$ denote a set of $m$ ordered (strictly monotonic increasing) indices known as change points with $1\le {\tau}_{i}\le n$ and ${\tau}_{m}=n$. For ease of notation we also define ${\tau}_{0}=0$. The $m$ change points, $\tau $, split the data into $m$ segments, with the $i$th segment being of length ${n}_{i}$ and containing ${y}_{{\tau}_{i1}+1:{\tau}_{i}}$.
Given a cost function,
$C\left({y}_{{\tau}_{i1}+1:{\tau}_{i}}\right)$,
g13nec gives an approximate solution to
where
$\beta $ is a penalty term used to control the number of change points. The solution is obtained in an iterative manner as follows:

1.Set $u=1$, $w=n$ and $k=0$

2.Set $k=k+1$. If $k>K$, where $K$ is a usersupplied control parameter, then terminate the process for this segment.

3.Find $v$ that minimizes

4.Test

5.If inequality (1) is false then the process is terminated for this segment.

6.If inequality (1) is true, then $v$ is added to the set of change points, and the segment is split into two subsegments, ${y}_{u:v}$ and ${y}_{v+1:w}$. The whole process is repeated from step 2 independently on each subsegment, with the relevant changes to the definition of $u$ and $w$ (i.e., $w$ is set to $v$ when processing the lefthand subsegment and $u$ is set to $v+1$ when processing the righthand subsegment.
The change points are ordered to give $\tau $.
4
References
Chen J and Gupta A K (2010) Parametric Statistical Change Point Analysis With Applications to Genetics Medicine and Finance Second Edition Birkhäuser
5
Arguments

1:
$\mathbf{n}$ – Integer
Input

On entry: $n$, the length of the time series.
Constraint:
${\mathbf{n}}\ge 2$.

2:
$\mathbf{beta}$ – double
Input

On entry:
$\beta $, the penalty term.
There are a number of standard ways of setting
$\beta $, including:
 SIC or BIC
 $\beta =p\times \mathrm{log}\left(n\right)$.
 AIC
 $\beta =2p$.
 HannanQuinn
 $\beta =2p\times \mathrm{log}\left(\mathrm{log}\left(n\right)\right)$.
where
$p$ is the number of parameters being treated as estimated in each segment. The value of
$p$ will depend on the cost function being used.
If no penalty is required then set $\beta =0$. Generally, the smaller the value of $\beta $ the larger the number of suggested change points.

3:
$\mathbf{minss}$ – Integer
Input

On entry: the minimum distance between two change points, that is ${\tau}_{i}{\tau}_{i1}\ge {\mathbf{minss}}$.
Constraint:
${\mathbf{minss}}\ge 2$.

4:
$\mathbf{mdepth}$ – Integer
Input

On entry:
$K$, the maximum depth for the iterative process, which in turn puts an upper limit on the number of change points with
$m\le {2}^{K}$.
If
$K\le 0$ then no limit is put on the depth of the iterative process and no upper limit is put on the number of change points, other than that inherent in the length of the series and the value of
minss.

5:
$\mathbf{chgpfn}$ – function, supplied by the user
External Function

chgpfn must calculate a proposed change point, and the associated costs, within a specified segment.
The specification of
chgpfn is:
void 
chgpfn (Nag_TS_SegSide side,
Integer u,
Integer w,
Integer minss,
Integer *v,
double cost[],
Nag_Comm *comm, Integer *info)



1:
$\mathbf{side}$ – Nag_TS_SegSide
Input

On entry: flag indicating what
chgpfn must calculate and at which point of the Binary Segmentation it has been called.
 ${\mathbf{side}}=\mathrm{Nag\_FirstSegCall}$
 only $C\left({y}_{u:w}\right)$ need be calculated and returned in ${\mathbf{cost}}\left[0\right]$, neither v nor the other elements of cost need be set. In this case, $u=1$ and $w=\mathrm{n}$.
 ${\mathbf{side}}=\mathrm{Nag\_SecondSegCall}$
 all elements of cost and v must be set. In this case, $u=1$ and $w=\mathrm{n}$.
 ${\mathbf{side}}=\mathrm{Nag\_LeftSubSeg}$
 the segment, ${y}_{u:w}$, is a lefthand side subsegment from a previous iteration of the Binary Segmentation algorithm. All elements of cost and v must be set.
 ${\mathbf{side}}=\mathrm{Nag\_RightSubSeg}$
 the segment, ${y}_{u:w}$, is a righthand side subsegment from a previous iteration of the Binary Segmentation algorithm. All elements of cost and v must be set.
The distinction between
${\mathbf{side}}=\mathrm{Nag\_LeftSubSeg}$ and
$\mathrm{Nag\_RightSubSeg}$ may allow for
chgpfn to be implemented in a more efficient manner. See
Section 10 for one such example.
The first call to
chgpfn will always have
${\mathbf{side}}=\mathrm{Nag\_FirstSegCall}$ and the second call will always have
${\mathbf{side}}=\mathrm{Nag\_SecondSegCall}$. All subsequent calls will be made with
${\mathbf{side}}=\mathrm{Nag\_LeftSubSeg}$ or
$\mathrm{Nag\_RightSubSeg}$.

2:
$\mathbf{u}$ – Integer
Input

On entry: $u$, the start of the segment of interest.

3:
$\mathbf{w}$ – Integer
Input

On entry: $w$, the end of the segment of interest.

4:
$\mathbf{minss}$ – Integer
Input

On entry: the minimum distance between two change points, as passed to g13nec.

5:
$\mathbf{v}$ – Integer *
Output

On exit: if
${\mathbf{side}}=\mathrm{Nag\_FirstSegCall}$ then
v need not be set.
if
${\mathbf{side}}\ne \mathrm{Nag\_FirstSegCall}$ then
$v$, the proposed change point. That is, the value which minimizes
for
$v=u+{\mathbf{minss}}1$ to
$w{\mathbf{minss}}$.

6:
$\mathbf{cost}\left[3\right]$ – double
Output

On exit: costs associated with the proposed change point,
$v$.
If
${\mathbf{side}}=\mathrm{Nag\_FirstSegCall}$ then
${\mathbf{cost}}\left[0\right]=C\left({y}_{u:w}\right)$ and the remaining two elements of
cost need not be set.
If
${\mathbf{side}}\ne \mathrm{Nag\_FirstSegCall}$ then
 ${\mathbf{cost}}\left[0\right]=C\left({y}_{u:v}\right)+C\left({y}_{v+1:w}\right)$.
 ${\mathbf{cost}}\left[1\right]=C\left({y}_{u:v}\right)$.
 ${\mathbf{cost}}\left[2\right]=C\left({y}_{v+1:w}\right)$.

7:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
chgpfn.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
g13nec you may allocate memory and initialize these pointers with various quantities for use by
chgpfn when called from
g13nec (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).

8:
$\mathbf{info}$ – Integer *
Input/Output

On entry: ${\mathbf{info}}=0$.
On exit: in most circumstances
info should remain unchanged.
If
info is set to a strictly positive value then
g13nec terminates with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_USER_STOP.
If
info is set to a strictly negative value the current segment is skipped (i.e., no change points are considered in this segment) and
g13nec continues as normal. If
info was set to a strictly negative value at any point and no other errors occur then
g13nec will terminate with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NW_POTENTIAL_PROBLEM.
Note: chgpfn should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
g13nec. If your code inadvertently
does return any NaNs or infinities,
g13nec is likely to produce unexpected results.

6:
$\mathbf{ntau}$ – Integer *
Output

On exit: $m$, the number of change points detected.

7:
$\mathbf{tau}\left[\mathit{dim}\right]$ – Integer
Output
Note: the dimension,
dim, of the array
tau
must be at least
 $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left(\lceil \frac{{\mathbf{n}}}{{\mathbf{minss}}}\rceil ,{2}^{{\mathbf{mdepth}}}\right)$, when ${\mathbf{mdepth}}>0$;
 $\lceil \frac{{\mathbf{n}}}{{\mathbf{minss}}}\rceil $, otherwise.
On exit: the first
$m$ elements of
tau hold the location of the change points. The
$i$th segment is defined by
${y}_{\left({\tau}_{i1}+1\right)}$ to
${y}_{{\tau}_{i}}$, where
${\tau}_{0}=0$ and
${\tau}_{i}={\mathbf{tau}}\left[i1\right],1\le i\le m$.
The remainder of
tau is used as workspace.

8:
$\mathbf{comm}$ – Nag_Comm *

The NAG communication argument (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).

9:
$\mathbf{fail}$ – NagError *
Input/Output

The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_INT

On entry, ${\mathbf{minss}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{minss}}\ge 2$.
On entry, ${\mathbf{n}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{n}}\ge 2$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
 NE_USER_STOP

User requested termination by setting ${\mathbf{info}}=\u2329\mathit{\text{value}}\u232a$.
 NW_POTENTIAL_PROBLEM

User requested a segment to be skipped by setting ${\mathbf{info}}=\u2329\mathit{\text{value}}\u232a$.
7
Accuracy
Not applicable.
8
Parallelism and Performance
g13nec is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
g13ndc performs the same calculations for a cost function selected from a provided set of cost functions. If the required cost function belongs to this provided set then
g13ndc can be used without the need to provide a cost function routine.
10
Example
This example identifies changes in the scale parameter, under the assumption that the data has a gamma distribution, for a simulated dataset with $100$ observations. A penalty, $\beta $ of $3.6$ is used and the minimum segment size is set to $3$. The shape parameter is fixed at $2.1$ across the whole input series.
The cost function used is
where
$a$ is a shape parameter that is fixed for all segments and
${n}_{i}={\tau}_{i}{\tau}_{i1}+1$.
10.1
Program Text
10.2
Program Data
10.3
Program Results
This example plot shows the original data series and the estimated change points.