1\documentclass[10pt,a4paper]{article}
2\usepackage{a4wide}
3\usepackage{graphicx,color,alltt}
4\usepackage[authoryear,round,longnamesfirst]{natbib}
5\usepackage{hyperref}
6
7\definecolor{Red}{rgb}{0.7,0,0}
8\definecolor{Blue}{rgb}{0,0,0.8}
9\definecolor{hellgrau}{rgb}{0.55,0.55,0.55}
10
11\newcommand{\E}{\mbox{$\mathsf{E}$}}
12\newcommand{\VAR}{\mbox{$\mathsf{VAR}$}}
13\newcommand{\COV}{\mbox{$\mathsf{COV}$}}
14\newcommand{\p}{\mbox{$\mathsf{P}$}}
15\newcommand{\email}[1]{\href{mailto:#1}{\normalfont\texttt{#1}}}
16\newenvironment{smallexample}{\begin{alltt}\small}{\end{alltt}}
17
18\setlength{\parskip}{0.5ex plus0.1ex minus0.1ex}
19\setlength{\parindent}{0em}
20
21\bibpunct{(}{)}{;}{a}{}{,}
22
23\newcommand{\ui}{\underline{i}}
24\newcommand{\oi}{\overline{\imath}}
25
26\RequirePackage{color}
27\definecolor{Red}{rgb}{0.5,0,0}
28\definecolor{Blue}{rgb}{0,0,0.5}
29\definecolor{hellgrau}{rgb}{0.55,0.55,0.55}
30
31\hypersetup{%
32hyperindex,%
33colorlinks,%
34linktocpage,%
35plainpages=false,%
36linkcolor=Blue,%
37citecolor=Blue,%
38urlcolor=Red,%
39pdfstartview=Fit,%
40pdfview={XYZ null null null}%
41}
42
43
44\begin{document}
45
46\SweaveOpts{engine=R,eps=FALSE}
47%\VignetteIndexEntry{strucchange: An R Package for Testing for Structural Change in Linear Regression Models}
48%\VignetteDepends{strucchange}
49%\VignetteKeywords{structural change, CUSUM, MOSUM, recursive estimates, moving estimates, monitoring, R, S}
50%\VignettePackage{strucchange}
51
52\title{\texttt{strucchange}: An \textsf{R}  Package for Testing for Structural
53Change in Linear Regression Models}
54\author{Achim Zeileis \hspace{0.4cm} Friedrich Leisch \hspace{0.4cm} Kurt Hornik \hspace{0.4cm} Christian Kleiber}
55\date{}
56\maketitle
57
58\begin{abstract}
59This introduction to the \textsf{R} package \texttt{strucchange} is a (slightly)
60modified version of
61\cite{Z-papers:Zeileis+Leisch+Hornik:2002}, which
62reviews tests for structural change in linear regression models from
63the generalized fluctuation test framework as well as from the $F$ test (Chow
64test) framework. Since \cite{Z-papers:Zeileis+Leisch+Hornik:2002} various extensions
65were added to the package, in particular related to breakpoint estimation
66\citep[also know as ``dating'', discussed in][]{Z-papers:Zeileis+Kleiber+Kraemer:2003}
67and to structural change tests in other parametric models \citep{Z-papers:Zeileis:2006}.
68A more unifying view of the underlying theory is presented in \cite{Z-papers:Zeileis:2005}
69and \cite{Z-papers:Zeileis+Shah+Patnaik:2010}.
70
71Here, we focus on the linear regression model and
72introduce a unified approach for implementing tests from the fluctuation test
73and $F$~test framework for this model, illustrating how this approach
74has been realized in {\tt strucchange}.
75Enhancing the standard significance test approach the package contains
76methods to fit, plot and test empirical fluctuation processes (like CUSUM, MOSUM
77and estimates-based processes) and to compute, plot and test
78sequences of $F$ statistics with the sup$F$, ave$F$ and exp$F$ test.
79Thus, it makes powerful tools available to display information about
80structural changes in regression relationships and to assess their significance.
81Furthermore, it is described how incoming data can be monitored.
82\end{abstract}
83
84\noindent
85{\bf Keywords:} structural change, CUSUM, MOSUM, recursive estimates, moving
86estimates, monitoring, \textsf{R}, \textsf{S}.\\
87
88\section{Introduction} \label{sec:intro}
89
90The problem of detecting structural changes in linear regression relationships
91has been an important topic in statistical and econometric research.
92The most important classes of tests on structural change are the tests from the
93generalized fluctuation test framework \citep{Z:Kuan+Hornik:1995} on the
94one hand and tests based on $F$ statistics
95\citep{Z:Hansen:1992,Z:Andrews:1993,Z:Andrews+Ploberger:1994} on the other.
96The first class includes in particular the CUSUM and MOSUM tests and the
97fluctuation test, while the Chow and the sup$F$ test belong to the latter.
98A topic that gained more interest rather recently is to monitor structural
99change, i.e., to start after a history phase (without structural changes) to
100analyze new observations and to be able to detect a structural change as soon
101after its occurrence as possible.\\
102
103This paper concerns ideas and methods for implementing generalized fluctuation
104tests as well as $F$ tests in a comprehensive and flexible way, that reflects
105the common features of the testing procedures. It also offers facilities to
106display the results in various ways.\\
107
108This paper is organized as follows: In Section~\ref{sec:model} the standard
109linear regression model upon which all tests are based will be described and the
110testing problem will be specified. Section~\ref{sec:data} introduces a data set
111which is also available in the package and which is used for the
112examples in this paper. The following sections \ref{sec:fluctests},
113\ref{sec:Ftests} and \ref{sec:monitor} will then explain the tests, how they are
114implemented in {\tt strucchange} and give examples for each. Section
115\ref{sec:fluctests} is concerned with computing empirical fluctuation processes,
116with plotting them and the corresponding boundaries and finally with testing for
117structural change based on these processes. Analogously,
118Section~\ref{sec:Ftests} introduces the $F$ statistics and their plotting and
119testing methods before Section~\ref{sec:monitor} extends the tools from
120Section~\ref{sec:fluctests} for the monitoring case.
121
122
123\section{The model} \label{sec:model}
124Consider the standard linear regression model
125\begin{equation} \label{model1} y_i = x_i^\top \beta_i + u_i
126\qquad (i = 1, \dots, n), \end{equation}
127where at time $i$, $y_i$ is the observation of the dependent
128variable, $x_i = (1, x_{i2}, \dots, x_{ik})^\top$ is a $k \times 1$ vector of
129observations of the independent variables, with the first component equal to
130unity, $u_i$ are iid(0, $\sigma^2$), and $\beta_i$ is the $k \times 1$ vector of
131regression coefficients. Tests on structural change are concerned with testing
132the null hypothesis of ``no structural change''
133\begin{equation} \label{null-hypothesis}
134H_0: \quad \beta_i = \beta_0 \qquad (i = 1, \dots, n)
135\end{equation}
136against
137the alternative that the coefficient vector varies over time, with
138certain tests being more or less suitable (i.e., having good or poor power) for
139certain patterns of deviation from the null hypothesis.\\
140
141It is assumed
142that the regressors are nonstochastic with $||x_i|| = O(1)$ and that
143\begin{equation} \label{A5} \frac{1}{n} \sum_{i=1}^n x_i x_i^\top \quad
144\longrightarrow \quad Q\end{equation}
145for some finite regular matrix $Q$. These
146are strict regularity conditions excluding trends in the data which are assumed
147for simplicity. For some tests these assumptions can be
148extended to dynamic models without changing the main properties of the tests;
149but as these details are not part of the focus of this work they are omitted
150here.\\
151
152In what follows $\hat \beta^{(i, j)}$ is the ordinary least
153squares (OLS) estimate of the regression coefficients based on the
154observations $i+1, \dots, i+j$, and $\hat \beta^{(i)} = \hat \beta^{(0,i)}$ is
155the OLS estimate based on all observations up to $i$. Hence $\hat \beta^{(n)}$
156is the common OLS estimate in the linear regression model.
157Similarly $X^{(i)}$ is the regressor matrix based on all observations up to $i$.
158The OLS residuals
159are denoted as $\hat u_i = y_i - x_i^\top \hat \beta^{(n)}$ with the
160variance estimate $\hat{\sigma}^2 = \frac{1}{n-k}
161\sum_{i=1}^n \hat u_i^2$. Another type of residuals that are often used in
162tests on structural change are the recursive residuals \begin{equation}
163\label{rr} \tilde u_i \; = \; \frac{y_i - x_i^\top \hat
164\beta^{(i-1)}}{\sqrt{1+x_i^\top \left(X^{(i-1)^\top} X^{(i-1)} \right)^{-1}x_i}}
165\qquad (i = k+1, \dots, n),\end{equation} which have zero mean and variance
166$\sigma^2$ under the null hypothesis. The corresponding variance estimate is
167$\tilde{\sigma}^2 = \frac{1}{n-k} \sum_{i=k+1}^n (\tilde u_i - \bar{ \tilde
168u})^2$.
169
170\section{The data} \label{sec:data}
171The data used for examples throughout this paper are macroeconomic time
172series from the USA. The data set contains the aggregated monthly personal
173income and
174personal consumption expenditures (in billion US dollars) between January 1959
175and February 2001, which are seasonally adjusted at annual rates.
176It was originally taken from \url{http://www.economagic.com/}, a web site for
177economic times series. Both time series are depicted in Figure
178\ref{fig:USIncExp}.
179
180\setkeys{Gin}{width=0.6\textwidth}
181\begin{figure}[htbp]
182\begin{center}
183<<data,fig=TRUE,height=4,width=6,echo=FALSE,results=hide>>=
184library("strucchange")
185data("USIncExp")
186plot(USIncExp, plot.type = "single", col = 1:2, ylab = "billion US$")
187legend(1960, max(USIncExp), c("income", "expenditures"),
188       lty = c(1,1), col = 1:2, bty = "n")
189@
190\caption{\label{fig:USIncExp} Personal income and personal consumption
191expenditures in the US}
192\end{center}
193\end{figure}
194
195The data is
196available in the {\tt strucchange} package: it can be loaded and a
197suitable subset chosen by
198<<subset>>=
199library("strucchange")
200data("USIncExp")
201USIncExp2 <- window(USIncExp, start = c(1985,12))
202@
203
204We use a simple error correction model (ECM) for the consumption function
205similar to \cite{sc:Hansen:1992a}:
206\begin{eqnarray} \label{ecm.model}
207\Delta c_t & = & \beta_1 + \beta_2 \; e_{t-1} + \beta_3 \; \Delta i_t + u_t, \\
208\label{coint.model} e_t & = & c_t - \alpha_1 - \alpha_2 \; i_t,
209\end{eqnarray}
210where $c_t$ is the consumption expenditure and $i_t$ the income. We estimate
211the cointegration equation (\ref{coint.model}) by OLS and use the residuals
212$\hat e_t$ as regressors in equation (\ref{ecm.model}), in which we will test
213for structural change. Thus, the dependent variable is the increase in
214expenditure and the regressors are the cointegration residuals and the
215increments of income (and a constant).
216To compute the cointegration residuals
217and set up the model equation we need the following steps in \textsf{R}:
218<<ecm-setup>>=
219coint.res <- residuals(lm(expenditure ~ income, data = USIncExp2))
220coint.res <- lag(ts(coint.res, start = c(1985,12), freq = 12), k = -1)
221USIncExp2 <- cbind(USIncExp2, diff(USIncExp2), coint.res)
222USIncExp2 <- window(USIncExp2, start = c(1986,1), end = c(2001,2))
223colnames(USIncExp2) <- c("income", "expenditure", "diff.income",
224                         "diff.expenditure", "coint.res")
225ecm.model <- diff.expenditure ~ coint.res + diff.income
226@
227
228Figure~\ref{fig:ts} shows the transformed time series necessary for estimation
229of equation (\ref{ecm.model}).
230
231\begin{figure}[htbp]
232\begin{center}
233<<ts-used,fig=TRUE,echo=FALSE>>=
234plot(USIncExp2[,3:5], main = "")
235@
236\caption{\label{fig:ts} Time series used -- first differences and
237cointegration residuals} \end{center}
238\end{figure}
239
240
241In the following sections we will apply the methods introduced to test for
242structural change in this model.
243
244
245\section{Generalized fluctuation tests} \label{sec:fluctests}
246
247The generalized fluctuation tests fit a model to the given data and derive an
248empirical process, that captures the fluctuation either in residuals or in
249estimates. For these empirical processes the limiting processes are known, so
250that boundaries can be computed, whose
251crossing probability under the null hypothesis is $\alpha$. If the empirical
252process path crosses these boundaries, the fluctuation is improbably large and
253hence the null hypothesis should be rejected (at significance level $\alpha$).
254
255\subsection{Empirical fluctuation processes: function \texttt{efp}}
256
257Given a formula that describes a linear regression model to be tested the
258function {\tt efp} creates an object of class {\tt "efp"} which contains
259a fitted empirical fluctuation process of a specified type. The types
260available will be described in detail in this section.\\
261
262{\bf CUSUM processes}:
263The first type of processes that can be computed are CUSUM processes, which
264contain cumulative sums of standardized residuals.
265\cite{Z:Brown+Durbin+Evans:1975} suggested to consider cumulative sums of
266recursive residuals:
267\begin{equation}\label{Rec-CUSUM} W_n(t) \; = \; \frac{1}{\tilde \sigma
268\sqrt{\eta}}\sum_{i=k+1}^{k +\lfloor t\eta \rfloor} \tilde u_i \qquad (0
269\le t \le 1),\end{equation}
270where $\eta = n-k$ is the number of recursive residuals and $\lfloor t\eta
271\rfloor$ is the integer part of $t\eta$.\\
272
273Under the null hypothesis the limiting process for the empirical fluctuation
274process $W_n(t)$ is the Standard Brownian Motion (or Wiener Process) $W(t)$.
275More precisely the following functional central limit theorem (FCLT) holds:
276\begin{equation}\label{lim(Rec-CUSUM)}  W_n \Longrightarrow W, \end{equation}
277as $n \rightarrow \infty$, where $\Rightarrow$ denotes weak convergence
278of the associated probability measures.\\
279
280Under the alternative, if there is just a single structural change point $t_0$,
281the recursive residuals will only have zero mean up to $t_0$. Hence the path of
282the process should be close to 0 up to $t_0$ and leave its mean afterwards.
283\cite{Z:Kraemer+Ploberger+Alt:1988}
284show that the main properties of the CUSUM quantity remain the same even under
285weaker assumptions, in particular in dynamic models. Therefore {\tt efp} has the
286logical argument {\tt dynamic}; if set to {\tt TRUE} the lagged observations
287$y_{t-1}$ will be included as regressors.\\
288
289\cite{Z:Ploberger+Kraemer:1992} suggested to
290base a structural change test on cumulative sums of the common OLS residuals.
291Thus, the OLS-CUSUM type empirical fluctuation process is defined by:
292\begin{equation} \label{OLS-CUSUM} W_n^0(t)
293\quad = \quad \frac{1}{\hat \sigma \sqrt{n}} \sum_{i=1}^{\lfloor nt \rfloor}
294\hat u_i \qquad (0 \le t \le 1). \end{equation}
295The limiting process for $W_n^0(t)$ is the standard Brownian bridge $W^0(t) =
296W(t) - t W(1)$. It starts in 0 at $t = 0$ and it also returns to 0 for $t =
2971$. Under a single structural shift alternative the path should have a peak
298around $t_0$.\\
299
300These processes are available in the function {\tt efp} by specifying the
301argument {\tt type} to be either {\tt "Rec-CUSUM"} or {\tt "OLS-CUSUM"},
302respectively.\\
303
304{\bf MOSUM processes}:
305Another possibility to detect a
306structural change is to analyze moving sums of residuals (instead of using
307cumulative sums of the same residuals).
308The resulting empirical fluctuation process does then not contain the sum of all
309residuals up to a certain time~$t$ but the sum of a fixed number of residuals in
310a data window whose size is determined by the bandwidth parameter $h \in (0,1)$
311and which is moved over the whole sample period. Hence the Recursive MOSUM
312process is defined by
313\begin{eqnarray} \label{Rec-MOSUM}
314M_n(t | h) & = & \frac{1}{\tilde \sigma \sqrt{\eta}} \sum_{i =
315k+\lfloor N_\eta t \rfloor +1}^{k+\lfloor N_\eta t \rfloor +\lfloor
316\eta h \rfloor} \tilde u_i \qquad (0 \le t \le 1-h) \\
317\label{Rec-MOSUM2}
318& = &  W_n \left( \frac{ \lfloor N_\eta t \rfloor + \lfloor \eta
319h\rfloor}{\eta} \right) - W_n \left( \frac{ \lfloor N_\eta t
320\rfloor}{\eta} \right), \end{eqnarray}
321where $N_\eta = (\eta - \lfloor \eta h \rfloor)/(1-h)$. Similarly the OLS-based
322MOSUM process is defined by
323\begin{eqnarray} \label{OLS-MOSUM}
324M_n^0 (t | h) & = & \frac{1}{\hat \sigma \sqrt{n}} \left( \sum_{i =
325\lfloor N_n t \rfloor +1}^{\lfloor N_n t \rfloor +\lfloor nh \rfloor} \hat u_i
326\right) \qquad (0 \le t \le 1 - h) \\
327\label{OLS-MOSUM2}
328& = & W_n^0 \left( \frac{ \lfloor N_n t \rfloor + \lfloor n
329h\rfloor}{n} \right) - W_n^0 \left( \frac{ \lfloor N_n t \rfloor}{n}
330\right),
331\end{eqnarray}
332where $N_n = (n - \lfloor n h \rfloor)/(1-h)$. As the representations
333(\ref{Rec-MOSUM2}) and (\ref{OLS-MOSUM2}) suggest, the limiting process
334for the empirical MOSUM processes are the increments of a Brownian motion or a
335Brownian bridge respectively. This is shown in detail in
336 \cite{Z:Chu+Hornik+Kuan:1995}.\\
337
338If again a single structural shift is assumed at $t_0$, then both MOSUM paths
339should also have a strong shift around $t_0$.\\
340
341The MOSUM processes will be computed if {\tt type} is set to {\tt "Rec-MOSUM"}
342or {\tt "OLS-MOSUM"}, respectively.\\
343
344{\bf Estimates-based processes}:
345Instead of defining fluctuation processes on the basis of residuals they
346can be equally well based on estimates of the unknown regression coefficients.
347With the same ideas as for the residual-based CUSUM- and MOSUM-type processes
348the $k \times 1$-vector $\beta$ is either estimated recursively with a growing
349number of observations or with a moving data window of constant bandwidth $h$
350and then compared to the estimates based on the whole sample.
351The former idea leads to the fluctuation process  in the spirit of
352\cite{Z:Ploberger+Kraemer+Kontrus:1989} which is defined by
353\begin{equation}
354\label{fluctuation} Y_n \left(t \right) = \frac{\sqrt{i}}{\hat \sigma
355\sqrt{n}}  \left({X^{(i)}}^\top X^{(i)} \right)^{\frac{1}{2}} \left( \hat
356\beta^{(i)} - \hat \beta^{(n)} \right), \end{equation}
357where $i = \lfloor k + t(n-k) \rfloor$ with $t \in [0,1]$. And the latter gives
358the moving estimates (ME) process introduced by \cite{Z:Chu+Hornik+Kuan:1995a}:
359\begin{equation} \label{ME}
360Z_n \left( \left. t \right| h \right) = \frac{\sqrt{\lfloor nh
361\rfloor}}{\hat \sigma \sqrt{n}} \left({X^{(\lfloor nt \rfloor, \lfloor nh
362 \rfloor)}}^\top X^{(\lfloor nt \rfloor,
363\lfloor nh \rfloor)} \right)^{\frac{1}{2}} \left( \hat \beta^{(\lfloor nt
364 \rfloor,\lfloor nh
365\rfloor)} - \hat \beta^{(n)} \right), \end{equation}
366where $0 \le t \le 1-h$.
367Both are $k$-dimensional empirical processes. Thus, the limiting processes are a
368$k$-dimensional Brownian Bridge or the increments thereof respectively. Instead
369of rescaling the processes for each $i$ they can also be standardized by
370$\left( {X^{(n)}}^\top X^{(n)} \right)^{\frac{1}{2}}$. This has the
371advantage that it has to be calculated only once, but \cite{Z:Kuan+Chen:1994}
372showed that if there are dependencies between the regressors the rescaling
373improves the empirical size of the resulting test. Heuristically the rescaled
374empirical fluctuation process ``looks'' more like its theoretic counterpart.\\
375
376Under a single shift alternative the recursive estimates processes should
377have a peak and the moving estimates process
378should again have a shift close to the shift point $t_0$.\\
379
380For {\tt type="fluctuation"} the function {\tt efp} returns the recursive
381estimates process, whereas for {\tt "ME"} the moving estimates process is
382returned.\\
383
384All six processes may be fitted using the function {\tt efp}.
385For our example we want to fit an OLS-based CUSUM process, and a moving
386estimates (ME) process with bandwidth $h = 0.2$. The commands are simply
387<<efp>>=
388ocus <- efp(ecm.model, type="OLS-CUSUM", data=USIncExp2)
389me <- efp(ecm.model, type="ME", data=USIncExp2, h=0.2)
390@
391These return objects of class {\tt "efp"} which contain mainly the
392empirical fluctuation processes and a few further attributes like the process
393type. The process itself is of class {\tt "ts"} (the basic time series class in
394\textsf{R}), which either preserves the time
395properties of the dependent variable if this is a time series (like in our
396example), or which is standardized to the interval $[0,1]$ (or a
397subinterval). For the MOSUM and ME processes the centered interval $[h/2,
3981-h/2]$ is chosen rather than $[0,1-h]$ as in (\ref{Rec-MOSUM}) and
399(\ref{OLS-MOSUM}).\\
400
401Any other process type introduced in this section can be fitted by setting the
402{\tt type} argument. The fitted process can then be printed, plotted or tested
403with the corresponding test on structural change. For the latter appropriate
404boundaries are needed; the concept of boundaries for fluctuation processes is
405explained in the next section.
406
407\subsection{Boundaries and plotting}
408
409The idea that is common to all generalized fluctuation tests is that the null
410hypothesis of ``no structural change'' should be rejected when the fluctuation
411of the empirical process $\mathit{efp}(t)$ gets improbably large compared to the
412fluctuation of the limiting process. For the one-dimensional residual-based
413processes this comparison is performed by some appropriate boundary $b(t)$, that
414the limiting process just crosses with a given probability $\alpha$. Thus, if
415$\mathit{efp}(t)$ crosses either $b(t)$ or $-b(t)$ for any $t$ then it has to be
416concluded that the fluctuation is improbably large and the null hypothesis can
417be rejected at confidence level $\alpha$. The procedure for the $k$-dimensional
418estimates-based processes is similar, but instead of a boundary for the process
419itself a boundary for $||\mathit{efp}_i(t)||$ is used,
420where $||\cdot||$ is an appropriate functional which is applied component-wise.
421We have implemented the functionals `max' and `range'. The null hypothesis is
422rejected if $||\mathit{efp}_i(t)||$ gets larger than a constant $\lambda$, which
423depends on the confidence level $\alpha$, for any $i = 1, \dots, k$.\\
424
425The boundaries for the MOSUM processes are also constants, i.e., of form $b(t) =
426\lambda$, which seems natural as the limiting processes are stationary. The
427situation for the CUSUM processes is different though. Both limiting processes,
428the Brownian motion and the Brownian bridge, respectively, are not stationary.
429It would seem natural to use boundaries that are proportional to the standard
430deviation function of the corresponding theoretic process, i.e.,
431\begin{eqnarray}
432\label{bound:Rec-CUS1} b(t) &=& \lambda \cdot \sqrt{t}\\
433\label{bound:OLS-CUS1} b(t) &=& \lambda \cdot \sqrt{t(1-t)}
434\end{eqnarray}
435for the Recursive CUSUM and the OLS-based CUSUM path respectively,
436where $\lambda$ determines the confidence level. But the
437boundaries that are commonly used are linear, because a closed form solution for
438the crossing probability is known. So the standard boundaries for the
439two proccess are of type
440\begin{eqnarray}
441\label{bound:Rec-CUS} b(t) &=& \lambda \cdot (1+2t)\\
442\label{bound:OLS-CUS} b(t) &=& \lambda.
443\end{eqnarray}
444They were chosen because they are tangential to the boundaries
445(\ref{bound:Rec-CUS1}) and (\ref{bound:OLS-CUS1}) respectively in $t = 0.5$.
446However, \cite{Zo:Zeileis:2000a} examined the properties of the alternative
447boundaries (\ref{bound:Rec-CUS1}) and (\ref{bound:OLS-CUS1}) and showed that
448the resulting OLS-based CUSUM test has better power
449for structural changes early and late in the sample period.\\
450
451Given a fitted empirical fluctuation process the boundaries can be
452computed very easily using the function {\tt boundary}, which returns a time
453series object with the same time properties as the given fluctuation process:
454<<efp-boundary>>=
455bound.ocus <- boundary(ocus, alpha=0.05)
456@
457It is also rather convenient to plot the process with its boundaries for some
458confidence level $\alpha$ (by default 0.05) to see whether the path exceeds the
459boundaries or not: This is demonstrated in Figure~\ref{fig:ocus}.
460
461\begin{figure}[hbtp]
462\begin{center}
463<<OLS-CUSUM,fig=TRUE,height=4,width=6>>=
464plot(ocus)
465@
466\caption{\label{fig:ocus} OLS-based CUSUM process}
467\end{center}
468\end{figure}
469
470It can be seen that the OLS-based CUSUM process exceeds its boundary; hence
471there is evidence for a structural change. Furthermore the process seems to
472indicate two changes: one in the first half of the 1990s and another one at the
473end of 1998.\\
474
475It is also possible to suppress the boundaries and add them afterwards, e.g. in
476another color
477<<efp-boundary2,eval=FALSE>>=
478plot(ocus, boundary = FALSE)
479lines(bound.ocus, col = 4)
480lines(-bound.ocus, col = 4)
481@
482For estimates-based processes it is only sensible to use time series plots if
483the functional `max' is used because it is equivalent to rejecting the null
484hypothesis when $\max_{i=1, \dots, k} ||\mathit{efp}(t)||$ gets large or when
485$\max_t \max_{i=1, \dots, k} \mathit{efp}_i(t)$ gets large. This again is
486equivalent to any one of the (one-dimensinal) processes $\mathit{efp}_i(t)$ for
487$i = 1, \dots, k$ exceeding the boundary. The $k$-dimensional process can also
488be plotted by specifying the parameter {\tt functional} (which defaults to {\tt
489"max"}) as {\tt NULL}:
490
491\begin{figure}[h]
492\begin{center}
493<<ME-null,fig=TRUE>>=
494plot(me, functional = NULL)
495@
496\caption{\label{fig:me-null} 3-dimensional moving estimates process}
497\end{center}
498\end{figure}
499
500The output from \textsf{R} can be seen in Figure~\ref{fig:me-null}, where the
501three parts of the plot
502show the processes that correspond to the estimate of the
503regression coefficients of the intercept, the cointegration residuals and the
504increments of income, respectively. All three paths show two shifts: the first
505shift starts at the beginning of the sample period and ends in about 1991 and
506the second shift occurs at the very end of the sample period. The shift that
507causes the significance seems to be the strong first shift in the process for
508the intercept and the cointegration residuals, because these cross their
509boundaries. Thus, the ME test leads to similar results as the
510OLS-based CUSUM test, but provides a little more information about the nature of
511the structural change.
512
513\subsection{Significance testing with empirical fluctuation processes}
514
515Although calculating and plotting the empiricial fluctuation process with its
516boundaries provides and visualizes most of the information, it might still be
517necessary or desirable to carry out a traditional significance test. This can be
518done easily with the function {\tt sctest} (\underline{s}tructural
519\underline{c}hange \underline{test}) which returns an object of class {\tt
520"htest"} (\textsf{R}'s standard class for statistical test results) containing
521in particular the test statistic and the corresponding $p$ value. The test
522statistics reflect what was described by the crossing of boundaries in the
523previous section. Hence the test statistic is $S_r$ from (\ref{statr}) for the
524residual-based processes and $S_e$ from (\ref{state}) for the estimates-based
525processes: \begin{eqnarray} \label{statr} S_r & = & \max_t
526\frac{\mathit{efp}(t)}{f(t)},\\ \label{state} S_e & = & \max
527||\mathit{efp}(t)||, \end{eqnarray}
528where $f(t)$ depends on the shape of the boundary, i.e.,
529$b(t) = \lambda \cdot f(t)$. For most boundaries is $f(t) \equiv 1$, but the
530linear boundary for the Recursive CUSUM test for example has shape $f(t) = 1 +
5312t$.\\
532
533It is either possible to supply {\tt sctest} with a fitted
534empirical fluctuation process or with a formula describing the model that should
535be tested. Thus, the commands
536<<efp-sctest,eval=FALSE>>=
537sctest(ocus)
538@
539and
540<<efp-sctest2>>=
541sctest(ecm.model, type="OLS-CUSUM", data=USIncExp2)
542@
543lead to equivalent results.
544{\tt sctest} is a generic function which has methods not only for
545fluctuation tests, but all
546structural change tests (on historic data) introduced in this paper including
547the $F$ tests described in the next section.
548
549
550\section{$F$ tests} \label{sec:Ftests}
551
552A rather different approach to investigate whether the null hypothesis
553of ``no structural change'' holds, is to use $F$ test statistics. An important
554difference is that the alternative is specified: whereas the generalized
555fluctuation tests are suitable for various patterns of structural changes, the
556$F$ tests are designed to test against a single shift alternative. Thus, the
557alternative can be formulated on the basis of the model (\ref{model1})
558\begin{equation} \label{H1}
559\beta_i = \left\{ \begin{array}{l} \beta_A \quad (1 \le i \le i_0) \\
560\beta_B \quad (i_0 < i \le n) \end{array} \right.,
561\end{equation}
562where $i_0$ is some change point in the interval $(k, n-k)$. \cite{Z:Chow:1960}
563was the first to suggest such a test on structural change for the case where the
564(potential) change point $i_0$ is known.
565He proposed to fit two separate regressions for the two subsamples defined by
566$i_0$ and to reject whenever
567\begin{equation} \label{Fstat}
568F_{i_0} \; = \; \frac{{\hat u}^\top {\hat u} - {\hat e}^\top {\hat
569e}}{{\hat e}^\top {\hat e}/(n-2k)}.
570\end{equation}
571is too large, where $\hat e = (\hat u_A, \hat u_B)^\top$ are the residuals from
572the full model, where the coefficients in the subsamples are estimated
573separately, and $\hat u$ are the residuals from the restricted model, where the
574parameters are just fitted once for all observations. The test statistic
575$F_{i_0}$ has an asymptotic $\chi^2$ distribution with $k$ degrees of freedom
576and (under the assumption of normality) $F_{i_0}/k$ has an exact $F$
577distribution with $k$ and $n-2k$ degrees of freedom. The major drawback of this
578``Chow test'' is that the change point has to be known in advance, but there are
579tests based upon $F$ statistics (Chow statistics), that do not require a
580specification of a particular change point and which will be introduced in the
581following sections.
582
583\subsection{$F$ statistics: function \texttt{Fstats}}
584
585A natural idea to extend the ideas from the Chow test is to calculate the $F$
586statistics for all potential change points or for all potential change points in
587an interval $[\ui, \oi]$ and to reject if any of those statistics get too large.
588Therefore the first step is to compute the $F$ statistics $F_i$ for $k < \ui \le
589i \le \oi < n-k$, which can be easily done using the function {\tt Fstats}.
590Again the model to be tested is specified by a formula interface and the
591parameters $\ui$ and $\oi$ are respresented by {\tt from} and {\tt to},
592respectively. Alternatively to indices of observations these two parameters can
593also be specified by fractions of the sample; the default is to take {\tt from =
5940.15} and implicitly {\tt to = 0.85}. To compute the $F$ test statistics for all
595potential change points between January 1990 and June 1999 the appropriate
596command would be:
597<<Fstats>>=
598fs <- Fstats(ecm.model, from = c(1990, 1), to = c(1999,6), data = USIncExp2)
599@
600This returns an object of class {\tt "Fstats"} which mainly contains a time
601series of $F$ statistics. Analogously to the empiricial fluctuation processes
602these objects can be printed, plotted and tested.
603
604\subsection{Boundaries and plotting}
605
606The computation of boundaries and plotting of $F$ statistics is rather similar
607to that of empirical fluctuation processes introduced in the previous section.
608Under the null hypthesis of no structural change boundaries can be computed
609such that the asymptotic probability that the supremum (or the mean) of the
610statistics $F_i$ (for $\ui \le i \le \oi$) exceeds this boundary is $\alpha$.
611So the following command
612\begin{figure}[htbp]
613\begin{center}
614<<Fstats-plot,fig=TRUE,height=4,width=6>>=
615plot(fs)
616@
617\caption{\label{fig:Fstats} $F$ statistics}
618\end{center}
619\end{figure}
620plots the process of $F$ statistics with its boundary;
621the output can be seen in
622Figure~\ref{fig:Fstats}. As the $F$ statistics cross their boundary, there is
623evidence for a structural change (at the level $\alpha = 0.05$). The process
624has a clear peak in 1998, which mirrors
625the results from the analysis by empirical fluctuation processes and tests,
626respectively, that also indicated a break in the late 1990s.\\
627
628It is also possible to plot the $p$ values instead of
629the $F$ statistics themselves by
630<<pval-plot,eval=FALSE>>=
631plot(fs, pval=TRUE)
632@
633which leads to equivalent results. Furthermore it is also possible to set up
634the boundaries for the average instead of the supremum by:
635<<aveF-plot,eval=FALSE>>=
636plot(fs, aveF=TRUE)
637@
638In this case another dashed line for the observed mean of the $F$ statistics
639will be drawn.
640
641\subsection{Significance testing with $F$ statistics}
642
643As already indicated in the previous section, there is more than one
644possibility to aggregate the series of $F$ statistics into a test statistic.
645\cite{Z:Andrews:1993} and \cite{Z:Andrews+Ploberger:1994} respectively suggested
646three different test statistics and examined their asymptotic distribution:
647\begin{eqnarray}
648\label{supF} \mbox{\textnormal{sup}}F & = & \sup_{\ui\le i \le \oi} F_i, \\
649\label{aveF} \mbox{\textnormal{ave}}F & = & \frac{1}{\oi - \ui+ 1}
650\sum_{i = \ui}^{\oi} F_i, \\
651\label{expF} \mbox{\textnormal{exp}}F & = & \log \left( \frac{1}{\oi - \ui+
6521} \sum_{i = \ui}^{\oi} \exp ( 0.5 \cdot F_i) \right).
653\end{eqnarray}
654The sup$F$ statistic in (\ref{supF}) and the ave$F$ statistic from
655(\ref{aveF}) respectively reflect the testing procedures that have been
656described above. Either the null hypothesis is rejected when the maximal or the
657mean $F$ statistic gets too large. A third possibility is to reject when the
658exp$F$ statistic from (\ref{expF}) gets too large. The ave$F$ and exp$F$ test
659have certain optimality properties \citep{Z:Andrews+Ploberger:1994}.
660The tests can be carried out
661in the same way as the fluctuation tests: either by supplying the fitted {\tt
662Fstats} object or by a formula that describes the model to be tested. Hence the
663commands
664<<Fstats-sctest,eval=FALSE>>=
665sctest(fs, type="expF")
666@
667and
668<<Fstats-sctest2>>=
669sctest(ecm.model, type = "expF", from = 49, to = 162, data = USIncExp2)
670@
671lead to equivalent output.
672
673The $p$ values are computed based on \cite{Z:Hansen:1997}.\footnote{The authors
674thank Bruce Hansen, who wrote the original code for computing $p$ values for $F$
675statistics in \textsf{GAUSS}, for putting his code at disposal for porting to
676\textsf{R}.}
677
678
679\section{Monitoring with the generalized fluctuation test}
680\label{sec:monitor}
681
682In the previous sections we were concerned with the retrospective detection
683of structural changes in \emph{given} data sets. Over the last years
684several structural change tests have been extended to monitoring
685of linear regression models where new data arrive over time
686\citep{Z:Chu+Stinchcombe+White:1996,Z:Leisch+Hornik+Kuan:2000}. Such
687forward looking tests are closely related to sequential tests. When
688new observations arrive, estimates are computed sequentially from all
689available data (historical sample plus newly arrived data) and compared
690to the estimate based only on the historical sample. As in the
691retrospective case, the hypothesis of no structural change is
692rejected if the difference between these two estimates gets too large.
693
694The standard linear regression model~(\ref{model1}) is generalized to
695\begin{equation}
696  y_i = x_i^\top \beta_i + u_i
697  \qquad (i = 1, \dots, n, n+1, \ldots),
698\end{equation}
699i.e., we expect new observations to arrive after time $n$ (when the
700monitoring begins). The sample $\{(x_1,y_1),\ldots,(x_n,y_n)\}$ will
701be called the \emph{historic sample}, the corresponding time period
702$1,\ldots,n$ the \emph{history period}.
703
704Currently monitoring has only been developed for recursive
705\citep{Z:Chu+Stinchcombe+White:1996} and moving
706\citep{Z:Leisch+Hornik+Kuan:2000} estimates tests. The respective
707limiting processes are---as in the retrospective case---the Brownian
708Bridge and increments of the Brownian Bridge. The empirical
709processes are rescaled to map the history period to the interval [0,1]
710of the Brownian Bridge. For recursive estimates there exists a closed
711form solution for boundary functions, such that the limiting Brownian
712Bridge stays within the boundaries on the interval $(1,\infty)$ with
713probability $1-\alpha$. Note that the monitoring period consisting of
714all data arriving after the history period corresponds to the Brownian
715Bridge after time 1. For moving estimates, only the growth rate of the
716boundaries can be derived analytically and critical values have to be
717simulated.
718
719Consider that we want to monitor our ECM during the
7201990s for structural change, using years 1986--1989 as the history
721period. First we cut the historic sample from the complete data set
722and create an object of class \texttt{"mefp"}:
723<<mefp>>=
724USIncExp3 <- window(USIncExp2, start = c(1986, 1), end = c(1989,12))
725me.mefp <- mefp(ecm.model, type = "ME", data = USIncExp3, alpha = 0.05)
726@
727Because monitoring is a sequential test procedure, the
728significance level has to be specified \emph{in advance}, i.e., when
729the object of class \texttt{"mefp"} is created. The \texttt{"mefp"}
730object can now be monitored repeatedly for structural changes.
731
732Let us assume we get new observations for the year 1990. Calling function
733\texttt{monitor} on \texttt{me.mefp} automatically updates our
734monitoring object for the new observations and runs a sequential test
735for structural change on each new observation (no structural break is
736detected in 1990):
737<<monitor1>>=
738USIncExp3 <- window(USIncExp2, start = c(1986, 1), end = c(1990,12))
739me.mefp <- monitor(me.mefp)
740@
741Then new data for the years 1991--2001 arrive and we repeat the
742monitoring:
743<<monitor2>>=
744USIncExp3 <- window(USIncExp2, start = c(1986, 1))
745me.mefp <- monitor(me.mefp)
746me.mefp
747@
748The software informs us that a structural break has been detected at
749observation \#72, which corresponds to December 1991. Boundary and
750plotting methods for \texttt{"mefp"} objects work (almost) exactly as
751their \texttt{"efp"} counterparts, only the significance level
752\texttt{alpha} cannot be specified, because it is specified when the
753\texttt{"mefp"} object is created. The output of
754\texttt{plot(me.mefp)} can be seen in Figure~\ref{fig:monitor}.
755
756\begin{figure}[htbp]
757\begin{center}
758<<monitor-plot,fig=TRUE,echo=FALSE,height=4,width=6>>=
759plot(me.mefp)
760@
761\caption{\label{fig:monitor} Monitoring structural change with bandwidth $h=1$}
762\end{center}
763\end{figure}
764
765Instead of creating an {\tt "mefp"} object using the formula interface like
766above, it could also be done re-using an existing \texttt{"efp"} object, e.g.:
767<<mefp2>>=
768USIncExp3 <- window(USIncExp2, start = c(1986, 1), end = c(1989,12))
769me.efp <- efp(ecm.model, type = "ME", data = USIncExp3, h = 0.5)
770me.mefp <- mefp(me.efp, alpha=0.05)
771@
772If now again the new observations up to February 2001 arrive, we can monitor
773the data
774<<monitor3>>=
775USIncExp3 <- window(USIncExp2, start = c(1986, 1))
776me.mefp <- monitor(me.mefp)
777@
778and discover the structural change even two observations earlier as we used
779the bandwidth {\tt h=0.5} instead of {\tt h=1}. Due to this we have not one
780history estimate that is being compared with the new moving estimates, but we
781have a history process, which can be seen on the left in
782Figure~\ref{fig:monitor2}. This plot can simply be generated by
783\texttt{plot(me.mefp)}.
784
785\begin{figure}[htbp]
786\begin{center}
787<<monitor-plot2,fig=TRUE,height=4,width=6,echo=FALSE>>=
788plot(me.mefp)
789@
790\caption{\label{fig:monitor2} Monitoring structural change with
791bandwidth $h=0.5$} \end{center}
792\end{figure}
793
794The results of the monitoring emphasize the results of the historic tests: the
795moving estimates process has two strong shifts, the first around 1992 and the
796second around 1998.
797
798\section{Conclusions} \label{sec:conclusions}
799
800In this paper, we have described the {\tt strucchange} package that implements
801methods for testing for structural change in linear regression relationships.
802It provides a unified framework for displaying information
803about structural changes flexibly and for assessing their significance
804according to various tests.\\
805
806Containing tests from the generalized fluctuation test framework as well as
807tests based on $F$ statistics (Chow test ststistics) the package extends
808standard significance testing procedures: There are methods for fitting
809empirical fluctuation processes (CUSUM, MOSUM and estimates-based
810processes), computing an appropriate boundary, plotting these results and
811finally carrying out a formal significance test. Analogously a sequence of $F$
812statistics with the corresponding boundary can be computed, plotted and tested.
813Finally the methods for estimates-based fluctuation processes have extensions to
814monitor incoming data.
815
816In addition to these methods for the linear regression model, the \texttt{strucchange}
817package contains infrastructure for testing, monitoring, and dating structural
818changes in other parametric models, e.g., estimated by maximum likelihood. Details
819about the underlying theory can be found in \cite{Z-papers:Zeileis:2005},
820\cite{Z-papers:Zeileis+Hornik:2007}, and \cite{Z-papers:Zeileis+Shah+Patnaik:2010}.
821The corresponding functions in \texttt{strucchange} are presented in
822\cite{Z-papers:Zeileis+Kleiber+Kraemer:2003} and \cite{Z-papers:Zeileis:2006}.
823
824
825\section*{Acknowledgments}
826
827The research of Achim Zeileis, Friedrich Leisch and Kurt Hornik was supported
828by the Austrian Science Foundation (FWF) under grant SFB\#010 (`Adaptive
829Information Systems and Modeling in Economics and Management Science').\\
830The work of Christian Kleiber was supported by the
831Deutsche Forschungsgemeinschaft, Sonderforschungsbereich 475.
832
833\bibliography{strucchange}
834\bibliographystyle{abbrvnat}
835
836\newpage
837\begin{appendix}
838\section{Implementation details for $p$ values}
839
840An important and useful tool concerning significance tests are $p$ values,
841especially for application in a software package. Their implementation is
842therefore crucial and in this section we will give more detail about the
843implementation in the {\tt strucchange} package.\\
844
845For the CUSUM tests with linear boundaries there are rather good
846approximations to the asymptotic $p$ value functions given in
847\cite{Zo:Zeileis:2000a}. For the recursive estimates fluctuation test
848there is a series expansion, which is evaluated for the first hundred terms. For
849all other tests from the generalized fluctuation test framework the $p$ values
850are computed by linear interpolation from tabulated critical values. For the
851Recursive CUSUM test with alternative boundaries $p$ values from the interval
852$[0.001, 1]$ and $[0.001, 0.999]$ for the OLS-based version respectively are
853approximated from tables given in \cite{Zo:Zeileis:2000}. The critical values
854for the Recursive MOSUM test for levels in $[0.01, 0.2]$ are taken from
855\cite{Z:Chu+Hornik+Kuan:1995}, while the critical values for the levels
856in $[0.01, 0.1]$ for the OLS-based MOSUM and the ME test are given in
857\cite{Z:Chu+Hornik+Kuan:1995a}; the parameter $h$ is in both cases interpolated
858for values in $[0.05, 0.5]$.\\
859
860The $p$ values for the sup$F$, ave$F$ and exp$F$ test are approximated based on
861\cite{Z:Hansen:1997}, who also wrote the original code in \textsf{GAUSS}, which
862we merely ported to \textsf{R}. The computation uses tabulated simulated
863regression coefficients.
864
865\end{appendix}
866
867\end{document}
868