Document Text (Pages 21-30) Back to Document

Discontinuities due to survey redesigns : a structural time series approach

by Goesjenov, Ibragim, MS

Page 21

an irregular term ɛt, where t = (1, . . . , t) is the observation period. In this paper for the victimization
series there are several observations at one time point, so the subscript k is used to refer to series k at
time point t.
In this paper the measurement error model as described in van den Brakel and Renssen (2005) is implemented.
They consider the observed population value obtained via complete enumeration under two or
more survey designs as the sum of the true population value and a bias induced by the survey design.
Hence, Yt,k = Y true
t,k + t,k, where Yt,k is the population value observed under a survey at time t, Y true
is the true unobserved population value at time t and t,k is the measurement bias introduced by the
survey at time t to measure the true population value Y true
t,k . More specific, using the decomposition as
proposed by Durbin and Koopman (2001), let

Y true
t,k = µt,k + ιt,k + γt,k + ɛt,k, with t = (1, . . . , t) and k = (1, . . . , k),

where ɛt,k N (0, σ2)

ɛ and k refers to the number of series observed at a time point. If a survey is

redesigned at time τ the value for the new survey Y ∗ ∗


becomes Yτ,k


= Yτ,k +
τ,k. Then, the difference
τ,k Yτ,k =
τ,k τ,k is the systematic difference between two survey approaches caused by the
introduction of a new survey as discussed in section 3.3.
As proposed by Durbin and Koopman (2001), it is possible to include a regression variable δt with its
time independent coefficient βk, such that explanatory variables can be allowed in the structural model
framework. Following van den Brakel and Roels (2010), in this paper the regression variable δt is used as
an intervention variable to quantify discontinuities. The intervention coefficient βk models the difference


Hence, the model boils down to

Yt,k = µt,k + ιt,k + γt,k + βkδt + ɛt,k, with t = (1, . . . , t) and k = (1, . . . , k). (4.1)

Note that (4.1) does not account for any sampling error. Following Scott and Smith (1974), Yt,k can
be considered as the true population value which can be properly described with the structural time
series model (4.1). Then, let ˆy(greg)
t,k = Yt,k + et,k. Here, ˆy(greg)
t,k is the parameter value estimated by
the generalized regression estimator as discussed in section 2.4, Yt,k the unobserved population value
described by (4.1) and et,k is the sampling error since not the entire population but a sample is observed.
Summarized, the estimated values for the target parameters become

t,k = µt,k + ιt,k + γt,k + βkδt + ɛt,k + et,k. (4.2)


Page 22

To improve readability, the superscript in ˆy(greg)
t,k will be left out in this section unless it is necessary.
Also, note that in this entire section the components containing the time index t should also contain the
label l to denote the survey type considered at time t. Since in the dataset considered there is only one
overlapping survey, the SMIV , the time index t is sufficient to indicate the survey type. In case SMIV is
used, it will be explicitly stated to to avoid ambiguities.

4.2.2 Application to victimization data

In this paper for the victimization series a relatively short time span of 19 years is considered. It is
therefore unlikely that the cyclic effects of (4.2) will have significant effects, and the same holds for the
seasonal effects since we have annual observations. Therefore, following the approach in van den Brakel
and Roels (2010), it has been decided to disregard the seasonal component ιt,k and the cyclic component

γt,k. For the dataset considered in this paper, it is possible to check for this assumption by looking at

the Durbin-Watson statistics or the ACF plots.
In (4.2), the sampling errors et,k can be correlated through time. As proposed by Binder and Dick (1990),
using the variances of ˆyt,k the following form can be obtained which allows for nonhomogeneous variance
in the sampling errors :


et,k =

V ar( )

ˆyt,k × ˜et,k, (4.3)

V ar( )

ˆyt,k is the standard error of ˆyt,k defined in section 2.4, and ˜et,k is an ARMA process

which models the serial correlation between the sampling errors. Since the standard error of ˆyt,k is not
available for all surveys considered in this paper, this approach cannot be implemented here. To account
for nonhomogeneity between the sampling errors it is assumed that the sampling errors are inversely
proportional to the sample size of the survey.
Since in this paper only cross sectional series are considered, it is not possible to separate et,k from ɛt,k,
which would be possible if a panel dataset was available. Therefore, it has been decided to combine the
sampling error etk and the irregular term ɛt,k into one irregular term ϕt,k = et,k + ɛt,k. The variance of

ϕt,k is inversely proportional to the sample size nt, such that ϕt,k (

= N 0, σ2


. This also implies that
the sampling error et,k dominates the irregular term ɛt,k (van den Brakel and Roels, 2010). Note that
since no panel surveys are considered in this analysis, it can be assumed that Cov(ϕt,k, ϕt,k) = 0 t ̸= t.
These simplifications result in the following linear Gaussian model for the structural time series analysis :

ˆyt,k = µt,k + βkδt + ϕt,k. (4.4)

Stochastic trend
The first term in equation (4.4), µt,k, is the component which models the stochastic trend in the model.


Page 23

More specific, the smooth trend model is assumed, which is defined as

µt,k = µt1,k + rt1,k, where rt,k = rt1,k + ηt,r,k.

Here µt,k is called the level component, rt,k the stochastic slope component and ηt,r,k
= N
0, σ2


the irregular component of the trend. The restriction that the irregular terms of the stochastic trend
are independent, Cov(ηt,r,k, ηt,r,k) = 0, t ̸= tand k ̸= k, holds in the first part of the analysis, and
is released when explanatory variables are introduced. Durbin and Koopman (2001) also considers the
local linear trend model, which also includes an irregular term in µt,k. It has been decided to disregard
this irregular term because it is not expected that the population is changing rapidly.
Intervention variable
As aforementioned, the terms βkδt are the core elements of this analysis, they model the systematic
difference between two different survey approaches, as discussed in section 3.3. To quantify the effects
of a survey redesign, the regression variable δt is used as an intervention variable with its regression
coefficient βk (van den Brakel and Roels, 2010). A first time application of an intervention analysis is the
paper by Harvey and Durbin (1986) on the implication of the British seat belt law. Depending on the
type of intervention, δt can be defined in several ways. In case of a level intervention, it is assumed that
the magnitude of the discontinuity is constant over time after a survey redesign occurred at time τ :

⎪⎨ 1 t τ

δt =

⎪⎩ 0 t < τ
. (4.5)

This approach models only one survey redesign, which is easily extended to more interventions as it is
the case in this paper :
t =

1 τi t < τi+1
0 otherwise
, (4.6)

where i refers to the ith survey redesign, τi is the time point of that survey redesign and τi+1 is the
time point of the next survey redesign. Other examples of intervention variables are the pulse and slope
intervention variables (Durbin and Koopman, 2001) :

⎪⎨ 1 t = τ

δt =

⎪⎩ 0 t ̸= τ

⎪⎨ 1 + t τ t τ
, δt = ⎪⎩ 0 t < τ

The regression coefficient of the intervention variable has no irregular term. There are examples when
a survey redesign also leads to an increase or decrease in ϕt,k as well. Model (4.3) can be used to incorporate
the estimated variances of the survey estimates as prior information. Another possibility to
account for that would be introducing separate variances of the measurement equations every time a


Page 24

new survey is introduced. More specific, if survey l1 has been redesigned at time τ to survey l2, then
V ar(ϕt,k) = σϕ,k,l1 if t < τ and V ar(ϕt,k) = σϕ,k,l2 if t τ. To test the equivalence of these variances,
a sufficient amount of observations per survey is needed. Due to the rather short time span of the victimization
series this is not possible, and as proposed by van den Brakel and Roels (2010) it has been
decided to disregard this effect.
There is one exceptional case which should be investigated more closely. Suppose a new survey is introduced
at time point τ. On top of that, suppose that at τ there is also a change in the trend of the
parameter of interest, say due to an economic crisis. If one would apply the intervention analysis, the
intervention coefficient will wrongly assign the effect of a crisis to the new survey. Therefore, a crucial assumption
is that there is no structural change in neither the stochastic trend nor in the possible seasonal
components, such that all discontinuities can be only due to the introduction of a new survey. To check
for this assumption, one could incorporate an explanatory variable which does contain any structural
changes and which trend is correlated with the trend of the parameter of interest. In this paper, this
explanatory variable is the police series as described in section 3.2.2.

4.2.3 Transition into a state space model

Up till now, the series as provided by Statistics Netherlands have been modeled in a structural
time series model. Note that the discontinuities have not been estimated yet, the components are only
classified. Using the Kalman filter the discontinuities can be estimated, as described in the next section. To
estimate the parameters of a structural time series model, it has to be put in a state space representation.
The transformation of model (4.4) to a state space model as it had been addressed by Durbin and
Koopman (2001) is described in this subsection. In detail, a linear Gaussian state space model represents
a set of k equations coupled in a time-varying system :

ˆyt = ct + Ztαt + ϕt,

αt+1 = dt + Ttαt + ηt.


Variable Dimension Variable Dimension
ˆyt, ct, ϕt k × 1 Zt k × m
αt, dt, ηt m × 1 Tt m × m
k =number of variables, m = dimension of the state variable.
Table 4 – Dimensions of the variables as described in the linear Gaussian state space form (4.7)

Here, k is the number of variables per time point t. For the victimization series, k = 6 except for the
subpopulations, where k = 5. The dimensions of the matrices and vectors are given in table 4.
Equation (4.7a) is called the measurement equation, and equation (4.7b) is referred to as the transition


Page 25

equation. Note that the constants ct and dt are set equal to zero in this paper, but they can be used to
implement known constant effects of developments.
The measurement equation describes how the observed series ˆyt depends on the underlying unobserved
state variable αt via the matrix Zt. The state variable αt contains the parameters of the level and the
slope and the regression coefficients of the intervention variables discussed in section 4.2.2. The transition
equation describes how the state variable αt develops over time via the matrix Tt. The matrices Zt and

Tt are referred to as system matrices. In case any of the matrices is time invariant, the time subscript is

To get from (4.4) to (4.7), the matrices get the following values :

αt = (µt,1, rt,1, . . . , µt,k, rt,k, β1, . . . , βk),


Zt = (Ik (1, 0)|δtIk), T = Blockdiag(Ttr, Ik), Ttr = Ik ( 1 1
0 1

) , (4.8b)

ϕt = (ϕt,1, . . . , ϕt,k), ηt = (0, ηt,r,1, . . . , 0, ηt,r,k, 0
k), E(ϕt) = 0k, E(ηt) = 03k,

Ht = Cov(ϕt) = 1 Diag(σ2
ϕ,1, . . . , σ2
Qt = Cov(ηt) = Diag(0, σ2
r,1, . . . , 0, σ2
r,k, 0


Here, Ik is a k×k identity matrix, denotes the Kronecker product, δt is defined in (4.5), Blockdiag(x, y)
creates a diagonal matrix with x and y on the diagonals and 0k is a column vector of order k. The
variances of the disturbances are summarized in covariance matrices Ht and Qt.
Note that (4.8) is only defined for one survey redesign. If a series has more than one redesign, the model
can be easily extended. In the case of the victimization series there are three survey redesigns, such that

αt = (µt,1, rt,1, . . . , µt,k, rt,k, βD1
1 , . . . , βD1

, βD2

1 , . . . , βD2

, βD3

1 , . . . , βD3
Zt = (Ik (1, 0)|(1, 1, 1) δi
tIk), T = Blockdiag(Ttr, I3k), Ttr = Ik ( 1 1
0 1


, (4.9a)
) , (4.9b)

ϕt = (ϕt,1, . . . , ϕt,k), ηt = (0, ηt,r,1, . . . , 0, ηt,r,k, 0
3k), E(ϕt) = 0k, E(ηt) = 05k,

Ht = Cov(ϕt) = 1 Diag(σ2
ϕ,1, . . . , σ2
Qt = Cov(ηt) = Diag(0, σ2
r,1, . . . , 0, σ2
r,k, 0


Recall that δi
t is defined in (4.6) for the 3 interventions denoted by D1, D2 and D3, so there are 3 × k
regression coefficients in αt.
Now, model (4.9) can be used to implement the Kalman filter equations as discussed in section 4.3. These
provide an estimate ˆαt for αt, which also contains the estimated coefficient for the intervention variables
βk, denoted by ˆ
βk. Then, the time series k after the moment of survey redesign can be adjusted for those


Page 26

estimated discontinuities as proposed by van den Brakel et al. (2008) with

˜yt,k = ˆyt,k ˆ
βk. (4.10)

Here ˜yt,k is the series where the quantified discontinuities are eliminated.

4.2.4 Extensions

As aforementioned, model (4.9) can be extended for several reasons. One could add explanatory
variables to better explain the trend or level of the series, or add data to investigate subpopulations.
Including an explanatory variable in the intervention variable
As discussed in section 3.1, with the introduction of the ISM in 2008, the SMIV was conducted parallel
to it in 2008-2010. The difference ˆyISM,t ˆySMIV ,t, which is the yearly difference between the ISM and
the SMIV surveys, can be seen as the discontinuity between the SMIV and the new ISM. It is possible
to use this difference as explanatory information in the form of a time varying intervention variable δi
More specific, in (4.6) δi
t is initially a constant after a survey is introduced, equal to 1. When including
the explanatory variable ˆydif,t = ˆyISM,t ˆySMIV ,t in the intervention variable, (4.6) should be adjusted
t =

ˆydif,t τi t < τi+1

0 otherwise
. (4.11)

Now (4.11) can be used as the intervention variable after the introduction of the ISM survey in 2008.

Including an explanatory variable as a new variable
Another way of implementing an explanatory variable is to model it simultaneously with the target
variables and let the errors correlate with the corresponding parameter of interest. In contrast with the
the previous approach, the explanatory series could also contain a discontinuity, which in turn also can
be quantified with the approach discussed in this section. It is not necessary that the series is available
at the point a survey is redesigned, such that any correlated series can be used.
The model as described in (4.9) is applied to the original series and the explanatory series simultaneously,
such that the number of variables in (4.9) becomes kexpl = k + 1. The explanatory variable used in this
paper is the police series as described in section 3.2.2. An additional constraint is that the irregular
terms of the stochastic trend are correlated, meaning that in (4.9e) Qt is not a diagonal matrix anymore
but contains off-diagonal values for the correlations between the explanatory variable and its correlated
target variable. The covariance matrix of the transition equation is implemented as Qt = CtDtC
t, where

Dt is diagonal and Ct is the lower triangular with ones on the diagonal and the correlations on the

other entries. This decomposition is known as the Cholesky decomposition and ensures that Qt is always
positive semidefinite. Implementing the police series as an explanatory variable results in the following


Page 27

state space model :

αt = (µt,1, rt,1, . . . , µt,kexpl, rt,kexpl, βD1
1 , . . . , βD1

Zt =

, βD2
1 , . . . , βD2, β

1 , . . . , βD3,


Ikexpl (1, 0)|Blockdiag((1, 1, 1) δi
tIk, δi
T = Blockdiag(Ttr, I3k+1), Ttr = I

kexpl (

1 1
0 1


, (4.12a)
)), (4.12b)

) , (4.12c)

ϕt = (ϕt,1, . . . , ϕt,k, ɛt,expl), ηt = (0, ηt,r,1, . . . , 0, ηt,r,kexpl, 0

E(ϕt) = 0kexpl, E(ηt) = 02kexpl+(3k+1), (4.12e)


Ht = Cov(ϕt) = Diag (σ2
ϕ,1, . . . , σ2
ϕ,k), σ2 )
ɛ,expl , (4.12f)
Dt = Cov(ηt) = Diag(0, d2
r,1, . . . , 0, d2
r,kexpl, 03k+1). (4.12g)

Here, δi
t is defined in (4.6), such that using (4.11) also the difference ˆyISM,tˆySMIV ,t can be implemented.
Recall from section 3.2.2 that the police series does not contain any sampling error, such that its irregular
term of the measurement equation is equal to ɛt,expl. Similarly, the variance of the measurement equation
of the police series is σ2
ɛ,expl. βexpl refers to the coefficient of the intervention variable for the police series.
Note that in case the series would not be correlated, Ct would be the identity matrix and Qt = Dt. In case
of the police series, Ct is an identity matrix with a coefficient c on the off-diagonal entry (total, police) of
Ct, where total and police refer to the series which are correlated. The correlation between the irregular
terms of the trends of the police series and the total victimization series is then computed as

Corr(ηtotal, ηpolice) =




+ dr,totaldr,police

where ηtotal and ηpolice refer to the irregular terms in the transition equation of both series. d2
r,total and
r,police are the elements in Dt for the total victimization and police series, respectively. Note that if
the trends of total victimization and the police series are cointegrated, then dr,police 0, such that

Corr(ηpolice, ηtotal) 1.

Benchmarking with series
If there is data available for separate subpopulations as well as for the national level, it might be interesting
to investigate effects on the different subpopulations. Suppose there are H subpopulations. The
victimization series considered in this paper are reported as a percentage. Then, the following relationship
between the series ˆyt at the national level and the parameters estimates ˆyt,h at subpopulation level
must hold :
ˆyt = ˆyt,h, (4.13)



where h refers to a subpopulation in H and Nt,h and Nt = H

h=1 Nt,h is the population size at time t at


Page 28

subpopulation and national level, respectively.
Applying the models as described in section 4.2.3 separately to the H subpopulations and the series
at national level might result in inconsistent series after adjustment for discontinuities as described in
(4.10). This is because the right proportion of each subpopulation at national level is not incorporated
in the regression coefficient for the intervention variables βk.
As discussed in van den Brakel and Roels (2010), there exist several methods to implement constraint
(4.13). One method is to consider the series at national and subpopulation level simultaneously and
implement (4.13) in the regression coefficients of the intervention variables. More specific, let β = Gβ
denote the transition equation for the regression coefficients of the total population and the H subpopulations,
with β = (β
National,k, β
1,k, . . . , β

H ,k) for every intervention. Here, βNational,k is the regression

coefficient for breakdown k of the intervention variable at the national level and β1,k, . . . , βH ,k are the
regression coefficients for breakdown k of the intervention variable at the subpopulation level. The matrix
G is defined as follows :

G = Ok×k f
t,H I

, (4.14)

1H Ok×k IH Ik

where Ok×k is a k × k zero matrix, ft,H =


, . . . , Nt,H



is a H × 1 vector with the subpopulation
proportions, 1H is a column vector with ones of order H and IH an identity matrix of order H .
Since the series at national level and subpopulation level are considered simultaneously in this approach,
for the gender subpopulations of the victimization series the number of variables k per time point t
becomes ksub = 3 × k, since H = 2. Here, first there are k variables on the national level, and then 2k
variables for the male and female subpopulations, respectively. Then, the state space model (4.9) with
the implemented constraint (4.13) for the gender subpopulations is represented below.

αt = (µt,1, rt,1, . . . , µt,ksub, rt,ksub, βD1
1 , . . . , βD1

sub, βD2
1 , . . . , βD2
ksub, βD3
1 , . . . , βD3

Zt = (Iksub (1, 0)|(1, 1, 1) δi
tIksub), T = Blockdiag(Ttr, I3 G), Ttr = Iksub ( 1 1
0 1


, (4.15a)
) , (4.15b)

ϕt = (ϕt,1, . . . , ϕt,ksub), ηt = (0, ηt,r,1, . . . , 0, ηt,r,ksub, 0

E(ϕt) = 0

ksub, E(ηt) = 05ksub,



Ht = Cov(ϕt) = Diag (σ
ϕ,1, . . . , σ2
ϕ,k), 1 (σ

ϕ,k+1, . . . , σ2
, (4.15e)

Qt = Cov(ηt) = Diag(0, σ2
r,1, . . . , 0, σ2
r,ksub, 0



Here, G is defined in (4.14). Note the changed structure of Ht, which is due to the assumption that the
irregular terms of the measurement equation are inversely proportional to the sample size. To meet this
assumption for the series at subpopulation level, the irregular terms for the subpopulations are divided
by the sample size of the respective subpopulations, denoted by nt,h. In this paper, nt,h is approximated


Page 29

by Nt,h
Nt × nt. It follows that once matrix (4.14) is implemented in (4.15b) the consistency requirement is
met (van den Brakel and Roels, 2010).

4.2.5 Likelihood Ratio Tests
In subsections 4.2.1 - 4.2.4 it was assumed that the irregular terms of both the measurement equation
and the transition equation are distributed independently. This means that in Ht and Qt the variances
are independent among the breakdowns. To obtain a more parsimonious model, one could assume that
the variances of the irregular terms of either the measurement equation or the transition equation are
equal among the breakdowns. Hence,

ϕ,k = σ2
ϕ,kor σ2
r,k = σ2
r,kk and k’ k, (4.16)

where σ2
ϕ,k and σ2
r,k are defined in the model descriptions in sections 4.2.1 - 4.2.4. It is possible to test
for this assumption using the classical Likelihood Ratio Test (LRT).

4.3 Kalman filter and parameter estimation

The linear Gaussian state space model (4.7) as implied by (4.9) is restated here with the necessary
assumptions for the irregular terms :

ˆyt = Ztαt + ϕt, ϕt
= N (0, Ht) ,

αt+1 = Tαt + ηt, ηt
= N (0, Qt) , (4.17)

where t = 1, . . . , t.

Note that both the covariance matrices are Normally i.i.d. distributed. In this section it is initially
assumed that the initial conditions for αt are α1 N(a1, P1), with a1 and P1 known.

4.3.1 Filtering

Using the state space model (4.17), it is possible to apply the Kalman filter to the victimization series
as invented by the Hungarian-American electrical engineer Rudolf Emil Kálmán, known for its use in
GPS and flight control systems. ”Filtering is aimed at updating our knowledge of the system as each
observation ˆyt comes in”(Durbin and Koopman, 2001).
Let the information set ˆYt denote the set of observations ˆy1, . . . , ˆyt. Then, the objective of the Kalman
filter is to determine the conditional distribution of αt+1 given ˆYt for t = 1, . . . , t. Because all distributions
are normal, conditional distributions of subsets of variables given other subsets are also normal and the
required distribution is therefore determined by at+1 = E(αt+1| ˆYt) and Pt+1 = V ar(αt+1| ˆYt). Note that


Page 30

from (4.17) it follows that αt N(at, Pt), given information set ˆYt. Then, box (4.18) shows the filtering
equations for the Kalman filter, which are the optimal estimators of at+1 and Pt+1 given the information
set ˆYt. The state vector estimator E(αt| ˆYt) and its error variance matrix are denoted by at|t and Pt|t,
respectively. For details of the derivations the reader is referred to chapter 4.2 of Durbin and Koopman

Ft = ZtPtZ
t + Ht,

at|t = at + PtZ
t vt Pt|t = Pt PtZ
tF1[ ]
t PtZt

vt = ˆyt Ztat,

at+1 = Tat|t, Pt+1 = TPt|tT+ Qt,


where a1 and P1 are known.
The recursions in box (4.18) provide optimal estimates for the conditional distribution of αt, starting from

a1 and P1 and then recursively updating the knowledge of the system every time a new observation ˆyt

comes in. vt are the so-called one-step forecast errors of ˆyt given the information set ˆYt, and Ft = V ar(vt)
is the variance matrix of vt. Also, vt’s are sometimes called innovations since they represent the new
part of ˆyt which is not known up to point t.

4.3.2 State smoothing

While filtering focuses on updating the system every time a new observation comes in and in such
only considers the information set ˆYt, the smoothing process aims at estimating the unobserved state
variables α1, . . . , αn given the entire information set ˆYt. Hence, the smoothing process also takes into
account information that becomes available after time point t. The state vector αt is estimated by its
conditional mean ˆαt and its error variance matrix Vt = V ar(αt ˆαt) is also provided. Using the Kalman
recursion in box (4.18), the recursive formulas for the fixed interval smoother are given in box (4.19) :

Lt = T TPtZ
t Zt, wt1 = Z
t vt + L

Nt1 = Z
t Zt + L
tNtLt, such that

ˆαt = at + Ptwt1,

Vt = Pt PtNt1Pt,


where t = t, . . . , 1, wt = 0 and Nt = 0.

Note that in (4.19) we proceed backwards from t to 1 to obtain ˆαt, whereas in box (4.18) we proceeded
forwards updating the knowledge of the system every time a new observation comes in. For the details
of the derivation of (4.19) and alternative approaches to smoothing, the reader is referred to chapter 4.3
of Durbin and Koopman (2001).


© 2009 All Rights Reserved.