A General Framework for Observation Driven
Time-Varying Parameter Models∗
Drew Creala , Siem Jan Koopmana,c , André Lucasb,c(a) Department of Econometrics, VU University Amsterdam(b) Department of Finance, VU University Amsterdam, and Duisenberg school of finance(c) Tinbergen Institute, Amsterdam

October 23, 2008

Abstract
We propose a new class of observation driven time series models that we refer to as Generalized Autoregressive Score (GAS) models. The driving mechanism of the GAS modelis the scaled likelihood score. This provides a unified and consistent framework for introducing time-varying parameters in a wide class of non-linear models. The GAS model
encompasses other well-known models such as the generalized autoregressive conditional
heteroskedasticity, autoregressive conditional duration, autoregressive conditional intensity and single source of error models. In addition, the GAS specification gives rise toa wide range of new observation driven models. Examples include non-linear regression
models with time-varying parameters, observation driven analogues of unobserved components time series models, multivariate point process models with time-varying parametersand pooling restrictions, new models for time-varying copula functions and models fortime-varying higher order moments. We study the properties of GAS models and provide
several non-trivial examples of their application.
Keywords: dynamic models, time-varying parameters, non-linearity, exponential family,marked point processes, copulas.
JEL classification codes: C10, C22, C32, C51.
∗

Corresponding author: André Lucas, VU University Amsterdam, FEWEB, De Boelelaan 1105, 1081 HV

Amsterdam, Netherlands, phone: +31 20 598 6039, fax: +31 20 598 6020, email: alucas@feweb.vu.nl.

1

Introduction

Time series models with time-varying parameters have been categorized by Cox (1981) into twoclasses: observation driven models and parameter driven models. In this paper we develop a
new, general framework for building observation driven time series models. In the observation
driven approach, time variation of the parameters is introduced by making them dependenton (functions of) their own lagged values, past observations, and exogenous variables. Although the parameters are stochastic, they are perfectly predictable given past information.
This greatly simplifies likelihood evaluation and explains why the models have become popularin the applied econometrics and statistics literature. Typical examples of observation driven
models are the generalized autoregressive conditional heteroskedasticity (GARCH) model ofEngle (1982), Bollerslev (1986) and Engle and Bollerslev (1986), the autoregressive conditional
duration and intensity (ACD and ACI, respectively) models of Engle and Russell (1998) and
Russell (2001), the dynamic conditional correlation (DCC) model of Engle (2002a) and the
Poisson count model discussed by Davis, Dunsmuir, and Streett (2003). Our approach encompasses many of the existing observation driven models as mentioned above. In addition, it gives
rise to a wide range of new models.
The alternative to observation driven models are parameter driven models. In parameter
driven models, the parameters are stochastic processes which are subject to their own sourceof error. Given past and concurrent observations, the parameters are not perfectly predictable.
Typical examples include the stochastic volatility (SV) model, see Shephard (2005), and thestochastic intensity models of Bauwens and Hautsch (2006) and Koopman, Lucas, and Monteiro(2008). Estimation is usually more involved for these models because the associated likelihood
functions are not available in closed-form. Exceptions include linear Gaussian state space models and discrete-state hidden Markov models, see Harvey (1989) and Hamilton (1989). In mostother cases, computing the likelihood function requires evaluating a high-dimensional integral
usually by simulation methods such as importance sampling and Markov chain Monte Carlo;see, e.g., Shephard and Pitt (1997). However, parameter driven models offer a conceptuallystraightforward way of introducing time-varying parameters in a wide class of non-linear andnon-Gaussian models.
It is widely accepted that in many settings of empirical interest, time variation in a selectionof parameters of the model is important for capturing the dynamic behavior of (multivariate)
time series processes. Obvious examples include the models for volatility clustering as often
found in financial return data, see, e.g., the surveys of Bollerslev, Engle, and Nelson (1994),
Ghysels, Harvey, and Renault (1996), Engle (2001), Shephard (2005) and Bauwens, Laurent,2and Rombouts (2006). There are also other examples. In the case of univariate and multivariate (marked) point processes, the intensities can be made time-varying, see, e.g., Engle
and Russell (1998), Russell (2001), Bauwens and Hautsch (2006) and Koopman, Lucas, andMonteiro (2008). In the case of copula models, the dependence parameters can be considered as time-varying, see Dias and Embrechts (2004), van den Goorbergh, Genest, and Werker(2005), and, in particular, Patton (2006). The higher order moments of the distribution canalso be made to vary over time, see Hansen (1994), Harvey and Siddique (1999), Jondeau andRockinger (2003) and Brooks, Burke, Heravi, and Persand (2005). Correlations in multivariate
density functions can be allowed to vary over time, see Engle (2002a). Finally, researchers haveconsidered settings where the density itself or its quantiles are time-varying, see Park and Qian
(2006) and Engle and Manganelli (2004).
The main contribution of this paper is the proposal of a common framework for developing
observation driven models with time-varying parameters. The main difficulty in formulatinga unified framework lies in the choice of a function that links the past observations to future
parameter values. Such a function should be applicable to a wide class of non-linear and nonGaussian models. In this paper, we argue that the scaled score of the model density at time t isan effective choice for the driving mechanism of the time-varying parameters. By choosing thescaling appropriately, standard observation driven models such as the GARCH, ACD, and ACI
models are recovered. The scaled score is equally applicable to non-standard models and leads
to the formulation of new observation driven models that mimic the behavior of parameter
driven models for which no observation driven equivalents have been specified.
We will refer to our observation driven model with a scaled score function as the generalized autoregressive score (GAS) model. The GAS model has similar advantages as the GARCH
model. Likelihood evaluation is straightforward. Extensions to asymmetric, long memory, andother more complicated dynamics can be considered without introducing further complexities.
Other general frameworks for observation driven models have been suggested in the literatureincluding the generalized linear autoregressive (GLAR) models of Shephard (1995), the generalized autoregressive moving average (GARMA) models of Benjamin, Rigby, and Stanispoulos(2003) and the vector multiplicative error models (MEM) of Cipollini, Engle, and Gallo (2006).
In contrast to these proposals, GAS models exploit the complete density structure rather than
only means and higher moments.
To illustrate the applicability of GAS models, we study a number of interesting, non-trivial
settings for observation driven models. We consider linear and nonlinear regression models with
time-varying coefficients a typical class of models that we can treat within the GAS framework.
An example is the Nelson and Siegel (1987) model for analyzing the term structure of inter3

est rates which emphasizes that GAS can also treat multivariate models, see Subsection 4.1.
Multivariate non-Gaussian models for pooled marked point-processes with a GAS specification
for latent factors driving the log-intensities is a new model specification that can be used forthe modeling of credit rating transitions, see Subsection 4.2. The GAS representations of unobserved components time series models are discussed in Subsections 4.3 and 4.4. A new classof time-varying copulas models based on the GAS framework is formulated in Subsection 4.5.
A challenging direction of research is to model higher-order moments of financial returns astime-varying processes. In Subsection 4.6 we show that GAS provides a generic tool to develop
models that have time-varying higher-order moments. A particular case is to consider linear
and GARCH models with Student t distributions where the degrees of freedom parameter ismade time-varying. Observing trade by trade transaction prices on a discrete grid leads to someinteresting research directions. For example, price changes can be viewed as realizations of amultinomial distribution which need to be subject to time-varying processes. In Subsection 4.7,
we discuss a GAS treatment for the modeling of discrete price changes. Finally, in Subsection
4.8 we discuss a methodology for dynamic mixture distributions based on the GAS framework.
The remainder of the paper is organised as follows. In Section 2 we provide the basic GAS
specification together with a set of motivating examples. Section 3 includes a discussion of
the statistical properties of GAS models. Section 4 contains a range of non-trivial examples of
GAS models, where we develop several new observation driven models. In Section 5 we provide
simulation evidence for the statistical properties of the estimators. Finally, Section 6 concludes
and provides directions for future research.

2

The model

2.1

Basic model specification

Denote yt as the dependent variable of interest, ft as the time-varying parameter vector, xt as avector of exogenous variables (covariates), all at time t, and θ as a vector of static parameters.
The available information set at time t consists of Y1t−1 = {y1 , . . . , yt−1 }, F1t−1 = {f1 , . . . , ft−1 }and X1t = {x1 , . . . , xt }. We assume that observation yt is generated by the observation densityp(yt |ft−1 , Y1t−1 , X1t , F1t−2 ; θ).

(1)Furthermore, we assume that the mechanism for updating the time-varying parameter ft isgiven by the familiar autoregressive updating equation
ft = ω +p−1
X

Ai st−i +

q
X
j=1

i=0

4

Bj ft−j ,(2)where ω is a vector of constants, coefficient matrices Ai and Bj have appropriate dimensions
for i = 0, . . . , p − 1 and j = 1, . . . , q while st is the scaled score function and depends on pastobservations Y1t−1 , the time-varying parameters in F1t−1 , and the static parameter θ. Furthermore, the unknown coefficients in (2) are functions of θ as well, that is ω = ω(θ), Ai = Ai (θ)
and Bj = Bj (θ). Our main contribution is the particular choice for the driving mechanism stthat is applicable uniformly over a wide class of densities and, possibly non-linear, models.
Equation (1) differs from a typical parameter driven model specification. In a parameter
driven model, the parameter ft would evolve subject to its own independent source of error,say ηt . In particular, st would be replaced by ηt in (2). Estimation of ft is then based on theconditional (filtered) density p(ft |Y1t , X1t , F1t−1 ; θ). For linear Gaussian state space models, thisdensity can be computed in closed form by the Kalman filter. In non-linear and non-Gaussianmodels, conditional densities are generally evaluated via simulation methods; see, e.g., Durbin
and Koopman (2001) and Doucet, de Freitas, and Gordon (2001). The simulations are often
most effective when they are based on second order expansions of the log observation density(1). We use the same intuition to update the time-varying parameter from ft−1 to ft via (2)withst = S(t, Y1t−1 , F1t−1 ) · ∇t = St−1 · ∇t ,

(3)

where
∇t = ∂ ln p(yt |ft−1 , Y1t−1 , X1t , F1t−2 ; θ) / ∂ft−1 ,St−1 = S(t, Y1t−1 , F1t−1 ; θ),(4)with time dependent scaling matrix S(·). Given the reliance of the driving mechanism in(2) on the scaled score vector (3), we let the equations (1) – (3) constitute the generalized
autoregressive score model with orders p and q which we abbreviate by GAS(p, q).
There are several intuitive choices for the scaling matrix that we investigate here. Our first
choice is to set St−1 equal to the (pseudo)-inverse information matrix based on the density (1),that is−1−1
St−1 = It−1
= Et−1 [∇t ∇′t ] = −Et−1

−1
 2∂ ln p(yt |ft−1 , Y1t−1 , X1t , F1t−2 ; θ).
′∂ft−1 ∂ft−1(5)

The updating mechanism (2) for ft reduces closely to a Gauss-Newton updating step for every
new observation yt that becomes available through time. Using this particular choice for scaling
the score vector, the GAS model encompasses the well-known observation driven GARCH,
ACD, and ACI models as well as most of the Poisson count models considered by Davis et al.
(2003). When the scaling matrix is the identity matrix, that is St−1 = I in (3), the recursion
captures models such as the autoregressive conditional multinomial (ACM) model of Russell
and Engle (2005). In addition, it gives rise to a number of useful observation driven models that5

have not been investigated before. We first give some introductory examples of GAS models.
In Section 4, we provide a more systematic review of new, non-trivial models within the GAS
family that provide interesting directions for future research.

2.2

Some examples

Example 1 (GARCH models): Consider the model yt = σt−1 εt where the Gaussian disturbance εt has mean zero and variance one while σt is a time-varying standard deviation. It is a−1
basic exercise to show that the GAS(1, 1) model with St−1 = It−1
for ft = σt2 reduces toft = ω + A0 yt2 − ft−1 + B1 ft−1 ,(6)which is equivalent to the standard GARCH(1,1) specification of Bollerslev (1986). However,
if we assume that εt follows a Student t distribution scaled to have variance one and with ν
degrees of freedom, that is εt ∼ tν , the GAS(1,1) specification for the conditional variance leads
to the updating equation
ft = ω + A0 1 + 3ν

−1



·



(1 + ν −1 )y 2 − ft−1
(1 − 2ν −1 )(1 + ν −1 yt2 /(1 − 2ν −1 ) ft−1 ) t+ B1 ft−1 .

(7)This collapses to (6) in case of the Gaussian distribution when ν −1 = 0. The recursion in(7), however, has an important difference with the standard t-GARCH(1,1) model of Bollerslev(1987) which has the Student t density in (1) and the updating equation (6). The denominatorof the second term in the right-hand side of (7) causes a more moderate increase in the variancefor a large realization of |yt | as long as ν is finite. The intuition is clear: if the errors are modeledby a fat-tailed distribution, a large realization in yt does not necessitate a substantial increasein the variance. The GAS updating mechanism for the model with Student t errors therefore
is substantially different from its familiar GARCH counterpart. We return to this example inmore detail in Section 4.
Example 2 (MEM, ACI, and ACD models): Consider the model yt = µt−1 εt where εt
has a gamma distribution with shape and rate parameters both expressed as a parameter γ.
Given standard properties of the gamma distribution, µt−1 is the conditional mean of yt . Let
ft = µt . The GAS(1,1) model in this case is
ft = ω + A0 /2 (yt − ft−1 ) + B1 ft−1 ,(8)which is equivalent to the multiplicative error model (MEM) proposed by Engle (2002b) andextended in Engle and Gallo (2006). The exponential distribution is a special case of the gamma
distribution making the ACD and ACI models a special case of MEM and consequently GAS.
6To see this, suppose yt is exponentially distributed with time-varying mean λt = exp(ft−1 ) anddensity function p(yt |ft−1 ) = λt exp(−λt yt ). The GAS(1,1) model has the updating equation
ft = ω + A0 (1 − yt exp(ft−1 )) + B1 ft−1 ,(9)which is equivalent to the standard ACI(1, 1) model of Russell (2001). If we parameterize the−1
expected duration rather than the intensity, so that p(yt |ft−1 ) = λ−1t exp(−λt yt ), and define

the factor as λt = exp(−ft−1 ), we obtain
ft = ω + A0 yt / exp(ft−1 ) + B1 ft−1 ,(10)which is equivalent to the standard ACD(1, 1) model of Engle and Russell (1998), see also
Russell and Engle (2005).
Example 3 (regression models): We next consider the linear regression model yt = xt βt−1 +εt with normally distributed disturbances εt ∼ N(0, σ 2 ). Let ft = βt , then the scaled score−1
function with St−1 = It−1isst = xt (yt − xt ft−1 ) / (x′t xt ),(11)

and the GAS(1, 1) specification becomes
ft = ω + A0 xt

yt − xt ft−1
+ B1 ft−1 .
x′t xt(12)

For the time-varying intercept βt with xt ≡ 1, the updating equation (12) reduces to the
exponentially weighted moving average (EWMA) recursion by setting ω = 0 and B1 = 1, thatis
ft = ft−1 + A0 (yt − ft−1 ).

(13)

In this case, we obtain the observation driven analogue of the local level (parameter driven)model,
yt = µt−1 + εt ,µt = µt−1 + ηt ,where the unobserved level component µt is modeled by a random walk process and the disturbances εt and ηt are mutually and serially independent, and normally distributed, see Durbin
and Koopman (2001). A direct link between the parameter and observation driven models isestablished when we set ηt = α(yt − µt−1 ) = αεt while in (13) we set α ≡ A0 and consider
ft−1 as the (filtered) estimate of µt−1 . The local level model illustration shows that the GAS
representation is closely related to the single source of error (SSOE) framework as advocatedby Ord, Koehler, and Snyder (1997). However, the GAS framework allows straightforward
extensions for this class of models. For example, the EWMA scheme in (13) can be extended

7

by including σ 2 as a time-varying factor and recomputing the scaled score function in (11) for2
the new time-varying parameter vector ft−1 = (βt−1 , σt−1)′ .

We note that for the linear regression model, the updating function (12), with ω = 0 and
B1 = I, does not reduce to the recursive least squares update which is given by!−1tXyt − xt ft−1
+ ft−1 ,with Pt = σ 2
x′j xj.
ft = Pt−1 xt ′xt Pt−1 xt + σ 2
j=1
The GAS updating function (12) however reveals that if x′t xt is close to zero, the GAS driving
mechanism might become unstable. In this case, the recursive least squares update remains
stable, at least when σ 2 is sufficiently large. As a remedy for such instabilities, we provide
an information smoothed variant of the GAS driving mechanism which we discuss in the next
section. Alternatively, we may want to consider using the identity matrix to scale the scorewith St−1 = I and st = xt (yt − xt ft−1 ).
Example 4 (Dynamic exponential family models): The exponential family of distributions has a long history in statistics due to many convenient mathematical properties. Parameter driven models for some members of this family, like the normal (GARCH), the exponential(ACI and ACD), and the multinomial (ACM), have been developed separately in the literature. Shephard (1995) and Benjamin et al. (2003) proposed parameter-driven models for other
distributions in this family. The GAS framework introduces time-varying parameters for othermembers of the exponential family in a natural way. As mentioned earlier, each GAS model
also has the advantage that it exploits the full density structure to update the time-varying
parameters.
The main obstacle for using GAS models is often computing the information matrix given theparameterization used. To facilitate this task, we provide Table 1 which contains the gradientsand elements of the information matrix for a variety of models in the exponential family. In
addition to the GARCH and MEM classes of models, the GAS framework also encompasses
the time-varying binomial models of Cox (1958) and Rydberg and Shephard (2003), the ACM
model of Russell and Engle (2005), and some of the Poisson models in Davis et al. (2003).
The latter three models can be obtained by scaling the relevant score vector from Table 1 with−1
either an identity scaling matrix St−1 = I or the matrix square root of St−1 = It−1.

2.3

Other GAS specifications

An important advantage of the GAS(p, q) specification is that its applicability is not restricted
to one specific model or choice of model parameterization. In contrast, the recursion scheme

8

Table 1: Dynamics for some Exponential Family Distributions
Distribution
Normal (1)2exp(−0.5(y−µ) )(2πσ 2 )1/2

ft

∇t

µt

0.5(yt − µt )/σt2−0.5σt−2 + 0.5σt−4 (yt − µt )2

σt2

ItIt,11 = 0.5σt−2
It,22 = 0.5σt−4It,12 = 0µt0.5(yt − µt )/σt2

ln(σt2 )−0.5 + 0.5σt−2 (yt − µt )2

Normal (2)2exp(−0.5(y−µ) )(2πσ 2 )1/2

It,11 = 0.5σt−2
It,22 = 0.5
It,12 = 0ln(λt )

1 − λt yt

Gamma

ln(αt )

αt (ln(yt ) − ln(βt ) − Ψ(αt , 1))y α−1 exp(−y/β)β α Γ(α)

ln(βt )

yt /βt − αtln(αit )Pαit (Ψ ( αjt , 1) − Ψ(αit , 1))Exponential

It = 1

λ exp(−λy)

It,11 = α2t Ψ(αt , 2)
It,22 = αtIt,12 = αtln(µt )

yt − µt

It,ii = αit (1 + Ψ(αit , 1) + αit Ψ(αit , 2)−PPΨ ( αjt , 1) − αit Ψ ( αjt , 2))PIt,ij = αit αjt Ψ ( αjt , 2)

Negativeln(rt )rt (ln(pt ) + Ψ(yt + rt , 1) − Ψ(rt , 1))

It,11 = rt2 (Ψ(rt , 2) − E[Ψ(rt + yt , 2)])

Binomialln(pt /(1 − pt ))rt (1 − pt ) − yt pt

Dirichlet+αit ln(yit )PoissonI = µt

−µ y

e

µ
y!

y+r−1 r
p (1 − p)ykMultinomial

Q
n!

yj
J
j=1 pj

P
y =n−
P
p =1−
y1 !···yJ !

·Jj<J yjJ

j<J pjIt,22 = rt (1 − pt )It,12 = −pt

ln



1−

Pp

it
J−1
j=1 pjtyit − npit

j = 1, . . . , J − 1

It,ii = npit (1 − pit )It,ij = −npit pjt

We have defined ∇t in (4) and It in (5). The (i, j) element of It is denoted by It,ij . We further note that
Ψ(x, k) = ∂ k ln Γ(x)/∂xk .

is applicable to a wide class of models that are characterized by a parametric likelihood specification. This turns out to be particularly relevant for the applications in Section 4, where
we generalize some well-known models with time-varying parameters outside their usual areaof application. For example, if the time-varying parameter is common across different observations, the specification in (3) gives an automatic and model consistent way to weight the
information provided by different observations.
A possible difficulty with using the inverse information matrix as the scaling matrix St−1
is that it may be difficult to compute for specific models. As an alternative one could use aunit scaling matrix St−1 ≡ I. In this case, the updating mechanism only uses the unscaled
gradient making it close to a steepest-descent optimization step of the likelihood at time t. Our

9experience, however, is that this type of updating mechanism is often less stable. Therefore,
we prefer to stick to the information scaled score step whenever possible.
A further useful feature of the GAS model in (3) is that under the correct model specification
st is a martingale difference, Et−1 [st ] = 0. This follows directly from the properties of the score′
vector. Due to scaling by the matrix St−1 , we also obtain Et−1 [st s′t ] = St−1 · It · St−1
. This−11
simplifies to It−1
and It−1 for St−1 = It−1
and St−1 = I, respectively. The first two moments ofthe driving mechanism st thus easily link up with the theoretical properties of the postulated
model density (1).
An important ingredient of the GAS model is the scaling of the score. As mentioned in
Example 3 above, in some cases scaling by the inverse information matrix may cause numericalinstability. The simplest example is given by a time-varying AR(1) model without intercept,
yt = φt−1 yt−1 + εt ,with εt standard normally distributed. The information scaled score step in this case reduces
2
to st = (yt − φt−1 yt−1 )/yt−1. The GAS updating scheme becomes numerically unstable if the2
process yt−1 is close to zero. In this case, the information matrix I = yt−1
is close to zero andst can jump to extreme values. To accommodate this numerical problem, we introduce a formcof information smoothing over the most recent stretch of observations, that is St−1 = (It−1
)−1

where
c
c
It−1
= αIt−2+ (1 − α)It−1 .

(14)for some 0 ≤ α ≤ 1. This is an EWMA smoothing scheme. Other weighting schemes forsmoothing the information are also possible. The smoothing parameter α determines the number of observations that St−1 takes into account. For α → 0, we recover the standard GAS
model with information scaling. For α → 1, the model tends to average the information over all
past observations. The optimal smoothing parameter could be fixed a priori, or be determinedfrom the data itself by treating α as part of the static parameter vector θ in the likelihood.
The basic dynamics of (2) may be further extended in obvious directions. For example, it
may be interesting to include exogenous variables in (2), or to generalize the evolution of ft toinclude other non-linearities such as regime-switching. In addition, it may be useful in some
applications to consider long-memory versions of (2), for example
ft = ω +∞X(i + d − 1)!
i=0

i!(d − 1)!

st−1−i ,for fractional integration parameter d < 1/2, such that we obtain a fractionally integrated
GAS or FIGAS model specification, similar to the ARFIMA and FIGARCH literature, see theseminal paper of Hosking (1981). We leave such extensions for future research.
10

3

Statistical properties

3.1

Estimation and inference

An important advantage of observation driven models is that parameter estimation is relatively
straightforward using maximum likelihood (ML). This also holds for the GAS model. For an
observed time series y1 , . . . , yn and adopting the standard prediction error decomposition, we
can express the maximization problem as
max
θ

n
X

ℓ(θ; yt , ft , Y1t−1 , X1t , F1t−1 ),

(15)

t=1

where ℓ(θ; yt , ft , Y1t−1 , X1t , F1t−1 ) = p(yt |ft , Y1t−1 , X1t , F1t−1 ; θ) for an observed value yt . Similar
to the GARCH model, the GAS model defines the filter for the time-varying parameters. This
makes likelihood evaluation particularly simple. It only requires the implementation of the GAS
updating function (2) and the evaluation of p(yt |ft , Y1t−1 , X1t , F1t−1 ; θ∗ ) for a particular value θ∗of θ.
It is possible to formulate recursions for computing the gradient of the likelihood with respect
to the static parameters θ. Gradient recursions for the GARCH model have been developed by
Fiorentini, Calzolari, and Panattoni (1996). For the GAS(1,1) specification, we obtain∂ ln pt ∂ ln pt ∂ft−1
dℓt
=
+
,
′
′
dθ
∂θ′
∂ft−1
∂θ′

(16)

~0
~1
 ∂B
∂ft
∂ω
∂st
∂ft−1
∂A
′
′
(17)
=
+
A
+
B
+
(s
⊗
I)
+
f
⊗
I
,
0
1
t
t−1
∂θ′
∂θ′
∂θ′
∂θ′
∂θ′
∂θ′
~t−1
∂st
∂∇t
∂S
′
(18)
=
S
+
(∇
⊗
I)
,
t−1
t
∂θ′
∂θ′
∂θ′
~ = vec(A) denotes thewhere ℓt = ℓ(θ; yt , ft , Y t−1 , X t , F t−1 ), pt = p(yt |ft , Y t−1 , X t , F t−1 ; θ), A
1

1

1

1

1

1

vector with the stacked columns of the matrix A and ⊗ is the Kronecker matrix product. Higher
order GAS specifications can be dealt with similarly by formulating the GAS model updating
equation in companion form. The log-likelihood derivatives can thus be computed simultaneously with the time-varying parameters ft . However, computing the analytic derivatives, inparticularly for (18), may be cumbersome. In practice, we therefore often turn to likelihood
maximization based on numerical derivatives.
The easiest way to conduct inference for GAS models is to apply a standard limiting result
and use the inverse information matrix at the optimum to compute standard errors and t-values
for the estimated parameters. In particular, if θ gathers the static parameters of the model, we
conjecture that under standard regularity conditions, the maximum likelihood estimator θ̂ of θis consistent and satisfies
T 1/2 (θ̂ − θ) → N(0, H −1),11with H = E[∂ 2 ℓ/∂θ∂θ′ ].
It is not clear, however, that standard statistical results apply directly. As an example,
even though {st } forms a martingale difference sequence, it is not directly evident that theGAS(1,1) model will be stable as long as |B| < 1. This follows because the volatility of st
also changes over time in a stochastic way. For example, the GAS specification per se does
not ex clude the variance of st from becoming unbounded if st becomes large or small. If themodel density is such that the inverse information matrix with respect to ft−1 is uniformly
bounded, standard stability results apply for |B| < 1. This holds for a number of examples
we discuss in Section 4. It is clear that given the generality of the GAS specification, theconditions for standard asymptotic theory to apply need to be checked on a case by case basis.
We leave this to future research and concentrate here on the conceptual issues. To provide someindications of the statistical convergence, we complement our empirical examples of Section 4
by a Monte-Carlo simulation experiment for a selected set of examples in Section 5.

3.2

Parameterization and identification issues

The GAS specification also allows freedom of choice regarding the model’s parameterization.
In the GARCH example from Section 2, the time-varying parameter was ft = σt2 . If one would
like to account for the positivity of σt2 , an obvious alternative would be to parameterize themodel in terms of ft = ln(σt2 ). After some manipulations, the GAS(1,1) specification for this‘new’ model is
ft = ω + A1 ·




yt2
− 1 + B1 ft−1 .
2
σt−1

(19)

The GAS dynamics automatically adapt to the choice of parameterization. In general, assume
that one prefers a different parameterization f˜t = h(ft ) for some invertable mapping h(·). Letḣt = ∂h(ft )/∂ft′ . Note that ḣt is deterministic given all information up to and including time t.
Let sft denote the GAS information scaled score step for parameterization ft . For well behaved
densities, the information matrix equals both the expected outer product of scores and the
expected second derivative of the log density. This allows the information scaled score step tobe written as
˜
sft =



−1

Et−1 [ḣt−1 ∇t ∇′t ḣ′t−1 ]ḣt−1 ∇t = ḣ′t−1 sft .

(20)

Reparameterizing the model thus reduces to rescaling the information scaled score step by thegradient of the mapping h(·) in each period.
Another important issue concerns parameter identification. Consider a model density of thefollowing form,p(yt ; Φft−1 ),
12

(21)

where ft follows a GAS(1,1) specification and Φ is a matrix of constants. For example, Φft
can be a vector of volatilities of a vector time series, while ft is the common factor driving
all of them. Then it is not possible to define both Φ and all GAS parameters ω, A, and B,
simultaneously. Take the model in (21) and introduce an invertible matrix K. Define f˜t = Kft ,s̃t = Kst , Φ̃ = ΦK −1 , ω̃ = Kω, Ã = KAK −1 , B̃ = KBK −1 . The likelihoods for p(yt ; Φft−1 )and p(yt ; Φ̃f˜t−1 ) are obviously identical. Pre-multiplying the GAS(1,1) transition equation forthe original parameterization by K, we obtain
Kft = Kω + KBK −1 Kst + KAK −1 Kft−1 ⇔f˜t = ω̃ + Ãf˜t−1 + B̃s̃t−1 .

(22)
(23)

Introducing the notation ht = Φft = Φ̃f˜t , however, we note that for the parameterizationp(yt ; Φ̃f˜t−1 ) the GAS driver equalsE[∇f˜t−1 ∇′f˜t−1 ]−1 ∇f˜t−1 = E[Φ̃′ ∇ht−1 ∇′ht−1 Φ̃]−1 Φ̃′ ∇ht−1= KE[Φ′ ∇ht−1 ∇′ht−1 Φ]−1 Φ′ ∇ht−1= s̃t .
The GAS equation for the model p(yt ; Φ̃f˜t−1 ) is thus also given by (23). To solve this identification issue, restrictions must be imposed on the matrix Φ. For example, one can equate
specific rows of Φ to corresponding rows from the identity matrix. Note that the issue cannot
be solved by only imposing restrictions on the matrix A in the GAS equation, for example
A = I. To see this, put A = I, but let K be an arbitrary invertable matrix. The normalization
constraint then also holds for the re-parameterized model as Ã = KAK −1 = KK −1 = I. The
other parameter values, however, are clearly different (e.g., Φ versus Φ̃ = ΦK −1 ), while thelikelihood value remains unchanged. The fact that the identification issue can be solved via
restrictions on Φ but not on A in this case is a direct effect of the equivariance of the score andinformation matrix as a basis for the recursions in the GAS model. Care is therefore needed innormalizing the parameter spaces of models with a factor structure, such as (21).

4Applications and new models

In this section we apply the GAS model to a number of empirical examples. The examples
provide more insight into the potential of the GAS specification for various models. We also
provide several new models and non-trivial extensions of existing models.

13

4.1

Time-varying nonlinear regression models

The term structure of interest rates plays a central role in both macroeconomics and finance as
it describes the linkage between monetary policymakers’ impact on the short term interest rate
and firms’ investment decisions at longer horizons. We develop a parameter driven analogue
for the popular term structure model of Nelson and Siegel (1987) given by the partial nonlinear
regression modelyt,τ = xτ (λ) β + εt,τwhere yt,τ is the interest rate at time t for an investment that matures after τ months. The
1 × 3 covariate vector xτ (λ) is defined asxτ (λ) = 1 , (λτ )−1 (1 − exp(−λτ )) , (λτ )−1 (1 − exp(−λτ )) − exp(−λτ ) ,where the coefficients λ and β are unknown and with independent disturbance εt,τ ∼ N(0, σ 2) fora given time t. For the particular choice of xτ (λ), the three coefficients in β can be interpreted
as the level, slope, and curvature of the term structure, respectively. The slope and curvaturefactors depend on the parameter λ that is nonlinear in yt,τ . At time t, interest rates for mmaturities can be observed such that τ = τj for j = 1, . . . , m. Based on these m observations(for a given t), parameters β, λ and σ 2 can be estimated by nonlinear least squares methods.
This estimation procedure can be repeated at each time t, resulting in a time series of parameter
estimates, see Diebold and Li (2006) for further details.
Recently, Diebold, Rudebusch, and Aruoba (2006) proposed analyzing the time series dimension simultaneously with the maturity dimension by considering the Nelson-Siegel model asa multivariate state space model. For this purpose, they treat β as a 3 × 1 vector of unobservables, βt , that are modeled as a vector autoregressive (VAR) process. Furthermore, the interest
rates for all m maturities at time t are put into the observation vector yt = (yt,τ1 , . . . , yt,τm )′and modeled by yt = Xt βt + εt where Xt = [xτ1 (λ)′ , . . . , xτm (λ)′ ]′ , with serially independent
m × 1 disturbance vector εt ∼ N(0, D) and m × m positive diagonal matrix D. The unknown
parameters in the VAR model for βt , as well as D and λ are then estimated by ML using the
Kalman filter in this parameter driven approach. In an observation driven approach, we set
ft−1 = βt and consider the Gaussian density for yt − Xt ft−1 ∼ N(0, D) with the GAS(1, 1)
updating equation ft = ω + A0 st + B1 ft−1 wherest = Xt′ D −1 Xt

−1Xt′ D −1 (yt − Xt ft−1 ) ,assuming that m > 3. The static parameter vector is given by
′

′
′
′
′
~
~
θ = diag(D) , ω , A0 , B1 , λ .
14

To illustrate our GAS approach to the dynamic Nelson-Siegel model, we consider the dataset from Diebold et al. (2006) which consists of monthly U.S. Treasury yields with maturities
3, 6, 9, 12, 15, 18, 21, 24, 30, 36, 48, 60, 72, 84, 96, 108, and 120 months over the period
from January 1972 to December 2000. For comparison purposes, we have estimated both theparameter driven (DNS) and observation driven (GAS) models by ML. The DNS estimates areclose to those reported by Diebold et al. (2006) while the GAS estimates are different thanthose obtained for the parameter driven DNS model. For example, the estimate of λ is 0.0778
in the DNS model while it is 0.0948 in the GAS model. The estimates of other coefficientsare also different which emphasizes that the interpretation of “comparable” coefficients in bothmodels are different. Nevertheless, the estimates of the three factors in βt are similar as we
observe from the graphs in Figure 1. The three estimated factors in both models accurately
track the changes of the yield curve over time.
15

10

GAS factor
DNS factor
Data approx

Level

5

5

1975

1980

1985

1990

1995

2000

1975

1980

1985

1990

1995

2000

1980

1985

1990

1995

2000

Slope

0
−5

5

Curvature

0
−5
1975

Figure 1: Three factor dynamic Nelson-Siegel model: level, slope and curvature factors. Each panel contains
the estimated factor from the GAS(1,1) observation driven model, the estimated (filtered) factor from the DNS
parameter driven model, and the data proxy which for the level factor is the 120 month yield, for the slope factor
is the spread from the 3 month yield minus the 120 month yield and for the curvature factor is the 24 monthyield minus the 3 and 120 month yield.

15

4.2

Pooled marked point-process models

Models with time-varying intensities have been given much attention in the finance and microeconometric literature. The principal areas of application in economics include intraday trade
data (market microstructure), defaults of firms, credit rating transitions and (un)employment
spells over time. To illustrate the GAS model in this setting, we consider an application fromthe credit risk literature in which pooled marked point-processes are playing an important role.
We develop a new and useful modeling framework that is based on the GAS specification.
Recently, a number of promising models with stochastically evolving intensities have been
proposed, see Bauwens and Hautsch (2006), Koopman, Lucas, and Monteiro (2008), Duffie,
Eckner, Horel, and Saita (2006), and Koopman, Lucas, and Schwaab (2008). The econometric handling of these parameter driven models is intricate while parameter estimation can becomputationally demanding. In particular, likelihood evaluation for these models requires computing a (very) high-dimensional integral. Importance sampling techniques or Markov chain
Monte Carlo algorithms have been adopted for likelihood-based analyses of these models. The
use of such simulation-based techniques, however, obstructs the widespread application of thesemodels in practice. A computationally less-demanding alternative can be based on developing
observation driven analogues of these models.
The first step would then be to consider multivariate generalizations of Russell (2001).
However, this is not straightforward. Most of the models of Russell (2001) are developed in thecontext of high frequency data and in particular for stock trades. The structure of data sets oftrades is substantially different from the data sets that are used in credit risk. Whereas in highfrequency data one typically observes many spells for a limited number of stocks, in modelingcredit data one typically works with many different companies that only have very few spells
each. This requires the pooling of data over different companies in the sample. Consequently,
different events might carry information on where the dynamic parameter is heading at any
point in time. The GAS model provides an easy and consistent methodology to address thisissue.
Let yk (t)′ = (y1k (t), . . . , ynk (t)) be a vector of marks of n competing risk processes for firm
k = 1, . . . , N. We have yjk (t) = 1 if event type j materializes for firm k at time t, and zerootherwise. By following the application in Koopman, Lucas, and Monteiro (2008), we model
the log intensities of these processes by
λjk (t) = ηj + αj′ ft−1 ,(24)where ηj is the baseline intensity and αj is the vector of loadings for ft which is the vector ofdynamic factors and specified by the GAS(1,1) updating equation (2) with ω = 0. Since ft is an16

unobserved process, we may impose sign restrictions for αj to obtain economic interpretations
for the factors. This GAS specification states that the intensities of all firms are driven by
the same vector of time-varying systematic factors ft . Model (24) nests the model of Russell(2001) when we set the dimension of ft equal to the number of firms N. In a credit risk
context, we typically have dim(ft ) << N. Further parameter restrictions are needed for model
identification, see the discussion in Section 3. In the illustration below we set one of the αj ’sequal to unity.
The log-likelihood specification using (24) isℓt =Xyjk (t)λjk (t) − Rjk (t) · (t − t∗ ) · exp(λjk (t)),

(25)

j,k

where t∗ is the last event time before t, and Rjk (t) is a zero-one variable indicating whether
company k is potentially subject to risk j at time t. Based on the first and second derivativeof ℓt , we obtainst+1 ="
Xj,kwjk (t)αj αj′#−1where wjk (t) = Rjk (t) · exp(λjk (t)) /X

yjk (t)αj − Rjk (t) · (t − t∗ ) · exp(λjk (t))αj

j,k

!

,

(26)

P

j,k Rjk (t) · exp(λjk (t)) = P[yjk (t) = 1] is the probabilityof the next event being of type j for company k. Combining all these elements into a GAS
specification, we have obtained a new observation driven model for credit rating transitions.
As an illustration, we adopt the model described above for the CreditPro 7.0 data set whichcontains the Standard and Poor’s (S&P) rating histories of all US corporates over the period
1981–2005. We distinguish two complementary credit rating classes: the investment grade (IG)and the sub-investment grade (SIG). The GAS(1,1) model has a univariate (single) factor ft
and the updating equation has the scaled score function (26). The resulting model is estimatedunder the restrictions α3 < 0 and α4 = 1. The estimation results are presented in Table 2.
The GAS parameter A0 is estimated close to unity which implies a persistent dynamic process
for ft . Given the estimates of αj , the downgrades appear to be most sensitive to the common
factor ft . In particular, the baseline downgrade from investment grade to default is small withan estimate of -7.4 while it is strongly sensitive to the common factor ft with a loading estimateof 1.2. Interestingly, the estimated pattern of the systematic intensity factor ft is close to theestimated pattern of the parameter driven model of Koopman, Lucas, and Monteiro (2008).
However, in a GAS framework we do not require computationally intensive methods such asimportance sampling for parameter and factor estimation.
It is straightforward in our GAS framework to generalize the model to a three-factor model.
In this case, A0 and B1 in the GAS updating equation become 3 × 3 matrices. To obtain17

Table 2: Estimation results for the parameters in the one-factor GAS(1,1) intensity model (24) in a twograde system, with α4 = 1, with the scaled scoring function (26) and based on the S&P ratings of all US
corporates between 1981 and 2005. The estimates are reported with asymptotic standard errors in parantheses
below the estimates. Parameter A0 is subject to a logistic transformation during estimation and α3 is subject toa identifying restriction α3 < 0.

IG → SIG IG → DEF

SIG → IG

SIG → DEF

j

1

2

3

4

η

-3.920-7.360-3.360-3.330(0.118)(0.353)(0.109)(0.217)

0.520

1.190

-0.4701.000(0.076)(0.330)(0.086)—–logit(A0 )

A0

B1

6.415

0.998

0.024

α

(0.537)(0.003)IG → SIG

IG → DEF

−3.0

−6

−3.5

−7

−4.0

−8

−4.5

−9
1985
SIG → IG

1990

1995

2000

2005

1985
1990
SIG → DEF

1995

2000

2005

1995

2000

2005

−2
−3.5
−3

−4.0

−4

−5
1985

1990

1995

2000

2005

1985

1990

Figure 2: Marked point-process model: the estimated intensities in a two-grade system of the GAS(1,1) modelwith three credit risk factors, based on the scaled score function (26) and using the S&P rating histories of US
corporates for the period 1981–2005.

18

interpretable dynamic factors and to acquire the necessary restrictions, we set the loading
vector αj equal to the jth column of a 4 × 4 diagonal matrix for j = 1, 2, 3 while α4 = α3 . This
parsimonious specification implies that upgrades and downgrades between IG and SIG havedifferent factors while transitions to default also have a distinct factor. The static parameter
vector θ contains the elements of A0 and B0 together with the baseline intensities ηj and thethree unknown coefficients α1 , . . . , α3 . These parameters can be estimated in the usual way by
ML. After parameter estimation, we obtain similar estimated patterns for the three factors in ft
as for the more involved parameter driven model of Koopman, Lucas, and Monteiro (2008). In
particular, we corroborate the finding that the dynamics of upgrades are substantially different
from those of downgrades and defaults as can be clearly viewed in Figure 2 where the estimatedintensities ηj + αj′ ft−1 are displayed.

4.3

Unobserved component models with a single source of error

Unobserved components or structural time series models are a popular class of parameter driven
models where the unobserved components (UC) have a direct interpretation, see Harvey (1989).
In this section, we describe observation-driven analogues to UC models. For a univariate timeseries y1 , . . . , yn , a univariate signal ψt can be extracted. The dynamic properties of ψt canbe broken into a vector of factors ft−1 that are specified by the updating equation (2). For
example, we can specify the signal as the sum of r factors, that isψt = f1,t−1 + . . . + fr,t−1(27)with ft = (f1,t , . . . , fr,t )′ . In the case r = 2, we can specify the first factor as a time-varying trend
component (random walk plus drift) and the second factor as a second-order autoregressive
process with possibly cyclical dynamics. For this decomposition we obtain the GAS(1,2) modelwith observation model yt = ψt + εt = f1,t−1 + f2,t−1 + εt , observation density p(yt |ψt ; θ) =N(f1,t−1 + f2,t−1 , σ 2 ) and updating equation
  





ω
a1
1 0
0 0
 st +  ft−1 +  ft−2 .
ft =   + 
0
a2
0 φ1
0 φ2

(28)

The constant ω is the drift of the random walk trend factor f1,t and the autoregressive coefficients φ1 and φ2 impose stationarity for the second factor f2,t . The scaled score function is
given byst = yt − ψt = yt − f1,t−1 − f2,t−1 = εt ,(29)

and can be interpreted as the single source of error. The static parameter vector θ, consistingof coeffients ω, a1 , a2 , φ1 , φ2 and σ, can be estimated straightforwardly by ML. The estimates
19

Table 3: Estimation results for the parameters in the trend-cycle GAS(1,2) decomposition model (27) with theupdating equation (28) and the scaled scoring function (29) based on quarterly log U.S. real GDP from 1947(1)to 2008(2). The estimates are obtained by ML and reported with asymptotic standard errors in parantheses
below the estimates. Furthermore, the ML estimates of parameters in the parameter driven trend-cycle UC
model (30)–(31) are reported which are based on the same data set.

ω

a1

a2

φ1

φ2

σ

log-like

GAS 0.825 0.723 0.563 1.328 -0.424 0.905 -324.51(0.043)(0.206)(0.202)(0.130)(0.142)(0.041)

UC 0.825 0.604 0.621 1.501 -0.573(0.040)(0.098)(0.112)(0.102)(0.106)–-324.06–of ft result in a decomposition of yt into trend, cycle, and noise. This GAS decomposition canbe regarded as the observation driven equivalent of the UC models of Watson (1986) and Clark(1989), who also aim to decompose macroeconomic time series into trend and cycle factors.
The UC trend-cycle decomposition model is then given by yt = f1,t + f2,t with
f1,t = ω + f1,t−1 + a1 ξ1,t ,ξ1,t ∼ N(0, 1),(30)

f2,t = φ1 f2,t−1 + φ2 f2,t−2 + a2 ξ2,t ,ξ2,t ∼ N(0, 1),(31)

where the disturbances ξ1,t and ξ2,t are mutually and serially independent.
To illustrate the GAS trend-cycle decomposition model, we consider the time series ofquarterly log U.S. real GDP from 1947(1) to 2008(2) obtained from the Federal Reserve Bankof St. Louis. The vector of static coefficients θ is estimated by ML and the results are reported
in Table 3. The estimated autoregressive polynomial for factor f2,t has roots in the complex
range and therefore factor f2,t has cyclical properties. We may interpret f2,t−1 as a real-timebusiness cycle indicator for time t which is displayed in Figure 3. To compare this indicatorwith the indicator produced by the Watson (1986) model, we also report the ML estimates of
the corresponding coefficients in an UC trend-cycle model. These estimates are obtained byusing the Kalman filter for likelihood evaluation. Parameter estimates for the UC model arereported in Table 3 and the one-step ahead predicted estimate of f2,t is plotted in Figure 3. We
find that the parameter estimates from each model correspond closely. The second factor fromeach model exhibits cyclical behavior and the growth rate of the trend is estimated to be thesame. Estimates of the GAS and UC cycle factors in Figure 3 are almost indistinguishable.
The GAS framework is sufficiently general to provide an observation driven alternative forthe decomposition of univariate and multivariate time series based on UC models including
models with trend, seasonal, cycle and irregular components. For example, the GAS updating
equation can also be designed to incorporate the trend and cycle dynamics as formulated by
205(i) GAS cycle

UC cycle

4
3
2
1
0
−1
−2
−3
−4
1950

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

2005

Figure 3: Trend-cycle example: estimated cycles from the GAS and UC trend-cycle models based on quarterly
log of U.S. real gdp from 1947(1) through 2008(1). NBER recession dates are indicated by the shaded regions.

Harvey and Jaeger (1993). Regression and intervention effects can also be incorporated in the
GAS specification, see the discussion in Subsection 2.1. Since the resulting GAS models areequivalent to single source of error models, we refer to Ord et al. (1997) for a more detailed
discussion on this class of models.

4.4

State space models with time-varying GAS parameters

In this subsection, we argue that the GAS framework can also be applied to the static parameters of linear Gaussian state space models. Similar ideas have been considered for combining
state space models with ARCH disturbances; see, e.g. Harvey, Ruiz, and Sentana (1992). To
illustrate its relevance, we consider the local level model as specified by
y t = µt + εt ,µt = µt−1 + ξt ,εt ∼ N(0, σε2 ),ξt ∼ N(0, σξ2 ),t = 1, . . . , n.

(32)

We treat the variances of the irregular and level disturbances, σε2 and σξ2 respectively, as GAS
factors in order to obtain a time-varying UC model. The resulting model is similar in spirit to
the model considered by Stock and Watson (2007) for forecasting quarterly U.S. inflation. The
details of our time-varying local level model are given below.
The variances σε2 and σξ2 of the local level model (32) can be replaced by two GAS factors to
obtain a model with time-varying variances. Since variances must remain positive, we specify
the two GAS factors as log-variances. Replacing the constant variances in (32) with the GAS

21

factors, new disturbances for the local level model are given byεt ∼ N {0, exp(f1,t−1 )} ,ξt ∼ N {0, exp(f2,t−1 )} ,t = 1, . . . , n,(33)

with ft = (f1,t , f2,t )′ . Conditional on ft−1 , the unobserved level µt in (32) remains linear in theobservations y1 , . . . , yt and therefore the Kalman filter can be adopted to produce the optimal
estimate of µt which is given by at = E(µt |ft−1 , Y1t−1 , F1t−2 ; θ) with its mean square error pt .
The log-likelihood function for the local level model (32) is given byℓ(θ) =nXℓt (θ; yt , Y1t−1 ),

t=1

1
1
1
ℓt = ℓt (θ; yt , Y1t−1 ) = − ln 2π − ln dt − e2t / dt ,
2
2
2

(34)

where the prediction error et and its variance dt are evaluated by the Kalman filter as given by
et = yt − at ,at+1 = at + kt vt ,dt = pt + exp(f1,t−1 ),k t = pt / d t ,pt+1 = (1 − kt )pt + exp(f2,t−1 ),

(35)
(36)

for t = 1, . . . , n with diffuse initializations a1 = 0 and p1 = κ while κ → ∞, see Durbin andKoopman (2001). The Kalman filter update can be carried out simultaneously with the GAS
updating equation (2) for ft and with the scaled score function defined by (3). In this case, we
have∇i,t =∂ℓt∂ exp(fi,t−1 )∂ℓt×= exp(fi,t−1 )∂ exp(fi,t−1 )∂fi,t−1∂ exp(fi,t−1 )(37)

for i = 1, 2. Given the Kalman filter equations (35)–(36), the latter term of (37) is evaluated
by

1
1
∂ℓt
= − d˙it + 2et ėit / dt + d˙it (et / dt )2 ,(38)
∂ exp(fi,t−1 )22where ėit = ∂et / ∂ exp(fi,t−1 ) and d˙it = ∂dt / ∂ exp(fi,t−1 ) are evaluated by the additional
recursionsėit = −ȧit ,d˙it = ṗit + 1(i = 1),ȧi,t+1 = (1 − kt )ȧit + k̇it vt ,k̇it = (ṗit − kt d˙it ) / dt ,ṗi,t+1 = −k̇it pt + (1 − kt )ṗit + 1(i = 2),

(39)
(40)

where 1(i = j) equals one if i = j and zero otherwise, with initializations ȧi1 = 0 and ṗi1 = 0, fort = 1, . . . , n and i = 1, 2. For the local level model, we therefore obtain two additional recursions
to evaluate the score and they can also be carried out simultaneously with the Kalman filter.
By following Harvey (1989, p 140–2), we approximate the information matrix by




1˙ ˙
∂ 2 ℓt2≈ exp(fi,t−1 ) exp(fj,t−1 )× dit djt / dt + ėit ėjt / dt , (41)It−1 (i, j) = −Et−1′∂fi,t−1 ∂fj,t−12where It−1 (i, j) is the (i, j) element of It−1 for i, j = 1, 2. The computation of this approximation is feasible given the additional recursions (39)–(40). However, in practice, the information22

matrix may become singular or close to singular. We therefore have adopted the EWMA
smoothing scheme for the information matrix to obtain the scaled score function, see the discussion in Subsection 2.3.
The Kalman filter, the additional recursions for the score, and the GAS updating equation
for ft can be carried out simultaneously. The parameter vector θ consists of the GAS updating
coefficients in ω, A0 , . . . , Ap−1, B1 , . . . , Bq and, possibly, the smoothing coefficient α of the
EWMA smoothing recursion for the information matrix. The estimation of θ is done by ML
via the maximization of ℓ(θ) in (34) with respect to θ. Given a value of θ, the recursions in realtime provide estimates of µt via the Kalman filter and estimates of ft via the score recursions
and the GAS updating equation simultaneously.
To illustrate the new GAS model, we consider a time series of quarterly U.S. Consumer
Price Index inflation from 1959(1) to 2007(2) obtained from the FRED database. The local level
model with a GAS(1,1) updating equation for the log-variances is adopted and the methodologyof estimation as discussed above is implemented. The EWMA smoothing scheme is applied tothe information matrix depends on α which is estimated as part of θ. The GAS coefficientmatrices A0 and B1 are chosen to be diagonal so that we need to estimate a total of sevencoefficients. The estimation results are given by




0.122
0.426
0
,
,
ω̂ = 
Â0 = 
0.003
0
0.081



B̂1 = 

0.569

0

0

0.916



,

α̂ = 0.674,

with the maximized loglikelihood value given by −372.97. The loglikelihood value for standard

local level model is given by −394.19 indicating a substantial improvement in fit with the GAS
model specification.
Panels (i) and (ii) of Figure 4 present the estimated factors exp(f1,t /2) and exp(f2,t /2),
which are the standard deviations for εt and ξt , respectively. The standard deviation of theobservation disturbance is moderate until the end of 2006 onwards at which time it has become
relatively high. The standard deviation of the level disturbance increased in the 1970’s during
the periods of higher inflation and then decreased steadily over the remaining sample. The
signal to noise ratio in our GAS framework is defined by the ratio qt = exp(f2,t ) / exp(f1,t ).
When it is low, the estimate of µt is based on a long range of past observations. When it ishigh, µt is estimated using only a small set of recent observations. The third graph in Figure 4
displays qt based on the estimate of ft . As the properties of the model suggested, the estimateof µt is only based on recent observations during the years of the oil-crisis, 1974–1976. From1980 onwards, the level of inflation is more stable and a longer stretch of past observations
are used in the estimator of the trend. Finally, the (filtered) estimate of µt is displayed ingraph (iv) of Figure 4. The GAS framework captures the overall development of U.S. inflation234(i) volatility irregular component(ii) volatility level component

1.5
3
1.0
2
0.5
1
1960

1970

1980

1990

2000

1960
15

(iii) signal−to−noise−ratio

1970

1980

1990

2000

(iv) U.S. inflation and level component

6
10
4
5
2
0
1960

1970

1980

1990

2000

1960

1970

1980

1990

2000

Figure 4: Estimation results for the time-varying local level model applied to quarterly U.S. inflation from1959(1) to 2007(2). (i) estimated volatility exp(f1,t /2) of the irregular component; (ii) estimated volatility
exp(f2,t /2) of the trend component; (iii) estimated signal-to-noise ratio qt ; (iv) U.S. inflation and estimatedtrend µt .

effectively. The estimated patterns of the time-varying standard deviations are very similar to
those obtained by Stock and Watson (2007) who used Markov chain Monte Carlo to estimate
the time-varying variances in their parameter driven model.

4.5

Dynamic copula models

Copulas have become popular over the last decade in the literature on financial risk management. A copula is a multivariate distribution function over a hypercube with uniform marginals.
The copula can be used to link marginal distributions into a multivariate distribution using
Sklar’s theorem. In this subsection, we demonstrate that the GAS framework provides newmodel specifications for simple copulas such as the bivariate Gaussian copula. We then illustrate some of the numerical extensions of the GAS specification to mixture copulas that allow
for asymmetric tail behavior.

24

We start with simple Gaussian copulas where the GAS model suggests an alternative dynamic structure compared to earlier suggestions in the literature. Patton (2006) introduced thenotion of time-varying copulas, see also Dias and Embrechts (2004) and van den Goorbergh,
Genest, and Werker (2005). Patton (2006) models1 the driving mechanism for the dynamic
bivariate Gaussian copula as
ft = ω + A

m−1X

Φ−1 (u1,t−i )Φ−1 (u2,t−i) + Bft−1 ,

(42)

i=0

where Φ−1 is the inverse normal distribution function, u1t and u2t are the probability integral
transforms using the univariate marginals, m is a smoothing parameter, and the (Gaussian)correlation parameter ρt is a transform of the factor ft , e.g., ρt = (1−exp(−ft ))/(1+exp(−ft )).
Equation (42) is intuitively appealing and builds on our understanding of covariances: if thetransformed marginals covary positively, the correlation should increase. The reverse holds if
the transformed marginals are of opposite sign.
By using the density of the Gaussian copula, we can derive the GAS specification for thetime-varying correlation parameter. Note that the score with respect to the correlation parameter is the same for the Gaussian copula and the bivariate normal. Our results thus also
directly apply to the Dynamic Conditional Correlation (DCC) framework of (Engle 2002a).
Define xt = Φ−1 (u1t )2 + Φ−1 (u2t )2 and yt = Φ−1 (u1t )Φ−1 (u2t ). For m = 1, Patton’s model(42) then reduces to
ft = ω + A · yt + B · ft−1 .

(43)

Deriving the score and information matrix of the bivariate normal for the transformed correlation parameter, one can easily show that the GAS(1,1) specification for ft equalsxt − 2A+ Bft−1 .
· yt − ρt−1 − ρt−1
ft = ω +(1 − ρt−1 )(1 + ρt−1 )1 + ρ2t−1(44)

The similarities and differences between (43) and (44) are clear. Both models are driven by yt
as positively clustered transformed marginals should increase the correlation parameter. The
additional scaling factor (1−ρt−1 )(1+ρt−1 ) in (44) is a consequence of modeling the transformed
correlation parameter ft rather than ρt directly. The additional ρt−1 term in parenthesis is just atranslation parameter such that yt −ρt−1 becomes a martingale difference. The most interesting
difference between the two models lies in the final term in parenthesis involving xt . Clearly,
xt − 2 is a martingale difference. Note that xt is large if we observe an extreme observationin u1t , u2t , or particularly in both. The effect of this depends on the current estimate of the1

We adapt Patton’s notation here slightly to correspond with the timing convention used in the current

paper, i.e., using ft−1 in the copula at time t rather than ft .

25

correlation parameter ρt−1 . If the correlation is positive, then xt−1 has a negative effect. In
this case the xt term offsets part of the effect of yt if the latter is positive, i.e., if yt correspondswith the current positive estimate of ρt−1 . If yt is negative, however, the xt term reinforces themagnitude of the GAS step triggered by yt .
The effects are most clearly visualized by plotting the GAS step for different values of(u1t , u2t ) and three different values of the correlation parameter, ρt−1 = −0.5, 0.2, 0.9. Thisis done in Figure 5. If we consider the plot for ρ = 0.9, we see two clear differences. First,
the GAS step results in a smaller increase in the correlation parameter along the u1t = u2t
axis. Particularly if uit are both large or small, the step based on yt alone results in a more
pronounced increase of the transformed correlation parameter. The same holds for the smaller
positive correlation parameter of ρ = 0.2.
ρ= 0.2; GAS

ρ= 0.9; GAS

1.0
0.5

0.5

1.0

−50

−5

0

5

0

−25

10

0

5

ρ=−0.5; GAS

1.0
0.5

0.5

1.0

0
−25

0

−50

−5
0.5

1.0

0.5

ρ= 0.9; Patton

10
5
0

1.0
0.5

1.0

ρ= 0.2; Patton5ρ=−0.5; Patton

0.5

1.0

1.0
0.5

0.5

1.0

1.0
0.5

0.5

1.0

Figure 5: Comparisons between the GAS and Patton drivers for the bivariate Gaussian copula as a functionof the uniforms (u1t , u2t ). The top panels contain the graphs for the GAS step in (44) for ρt−1 = −0.5, 0.2, 0.9(left, middle, right). The lower graphs contain the (re-centered) steps yt − ρt−1 of the Patton model, (43). The
vertical axes have the same scale for each column of graphs.

The much more striking feature, however, is the increased sensitivity along the anti-diagonal
for positive ρ. If the current estimate of ρ is positive, but one observes a combination of (u1t , u2t )26

that signals negative rather than positive dependence, the GAS specification is more sensitive
to this and is more inclined to rapidly adjust the current estimate of ρ downward. Analogously,
for negative values of ρ the left panels show that the effects are reversed. The GAS specification
becomes more sensitive to observations along the diagonal than the specification based on yt
alone.
To illustrate this empirically, we extend the example from Patton (2006). We use the daily
exchange rates between the German Mark (later Euro) and the US dollar, and between theJapanese Yen and the US dollar. The sample period is January 1986 through August 2008. We
model the log returns on the exchange rates by AR-GARCH models and obtain the transformed
sample of u1t and u2t as input for the Gaussian copula model. Apart from (43) and (44), we
also estimate an ad-hoc implementation of the DCC framework of Engle (2002a). In particular,
we model the correlation parameter directly using
ρt = ω + A · yt + B · ρt−1 .

(45)

To enforce stationarity, we estimate the logit transform of B. The results are presented inTable 4 and Figure 6.
1.0

German Mark (Euro) and British Pound versus US Dollar

German Mark (Euro) and Japanese Yen versus US Dollar

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2
GAS
DCC

0.0
−0.2

Patton

0.0
GAS
DCC

Patton1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008

1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 2006 2008

Figure 6: Comparisons of the correlation parameter for the GAS, Patton, and DCC drivers in (43)–(45).
The data are the marginal AR-GARCH transforms of log exchange rates for the German Mark-US dollar andJapanese Yen-US dollar (left panel) and for the German Mark-US dollar and British Pound-US dollar (rightpanel). The sample period is January 1986–August 2008.

Table 4 shows that the GAS specification increases the log-likelihood value 25 to 125 points
for the same number of parameters. The figures show the empirical estimates of the timevarying correlation. Judging on the value of the parameter B, the GAS specification is themost persistent one, followed by the DCC and the Patton specification. However, the increased
sensitivity of the score mechanism to correlation breakdowns, reveals an opposite pattern in the27

Table 4: Parameter estimates for the GAS, Patton, and DCC drivers in (43)–(45). The data are the marginal
AR-GARCH transforms of log exchange rates for the German Mark-US dollar and Japanese Yen-US dollar (leftpanel) and for the German Mark-US dollar and British Pound-US dollar (right panel). The sample period is
January 1986–August 2008. Confidence interval in parentheses for B, otherwise standard errors in parentheses.
103 ω

B
)
ln( 1−BA

B

log-like

German Mark (Euro)–US $, Japanese Yen–US $GAS
Patton
DCC

6.11

0.058

5.30

0.995

(2.48)(0.009)(0.37)(0.990,0.998)

-1.60

0.036

4.27

0.986

(0.85)(0.003)(0.10)(0.983,0.989)

1.03

0.008

4.65

0.991

(0.29)(0.001)(0.09)(0.989,0.992)

1218.161191.511184.13

German Mark (Euro)–US $, British Pound–US $GAS
Patton
DCC

12.55

0.082

4.97

0.993

(3.55)(0.008)(0.26)(0.988,0.996)

-0.97

0.025

4.71

0.991

(0.84)(0.002)(0.11)(0.989,0.993)

2.64

0.004

4.84

0.992

(0.39)(0.000)(0.11)(0.990,0.994)

2218.822090.422060.43

figures. Due to the sharpe decline at the edges as visualized in Figure 5, the GAS specification
reacts much more fiercely to exchange rate returns of opposite sign if the current correlation
estimate is positive. This is most clearly seen for the Mark-Pound example, but also the MarkYen example shows similar features at the end of 1993 and 2003. The DCC dynamics, and toa lesser extent the Patton dynamics, are much smoother in this sense. The difference between
the dynamics for the different specifications may be highly relevant for risk managers, wherechanges in correlations and in particular correlation breakdowns are a major concern.
Stepping away from the Gaussian copula, it is also interesting to see how the GAS specification generalizes to non-Gaussian copulas. We consider here a mixture of Clayton-type copulas.
Patton proposes a generally applicable driving mechanism for copula parameters,−1

ft = ω − m AmX|u1,t−i − u2,t−i | + Bft−1 ,

(46)

i=1

where ft captures the dependence between the coordinates. The intuition for (46) is clear.
If the most recent u1t and u2t are close together, this is a signal of strong dependence and,
therefore, ft is increased. Similarly ft is decreased if u1t and u2t are far apart.
28

Though the driving mechanism in (46) is intuitively straightforward, two issues are lessclear. First, (46) uses no information contained in the particular choice of the copula. As with
the Gaussian copula, such information may be helpful in specifying the dynamics. Second,although (46) provides an easy updating scheme for the bivariate case, the extension to themultivariate case is less obvious. In particular, if one has an Archimedian copula characterized
by a single dependence parameter, there are many different ways in which one could use thedifferences |uit − ujt | for i 6= j to update the dependence parameter. Equation (46) provides
little guidance as to how these different and possibly conflicting signals should be weighed.
The Clayton copula for our example is a member of the Archimedian family. Its specificationin dimension d is given by
C(u1 , . . . , ud) =

1−d+

d
X

u−α
i

i=1

!−1/α

.

(47)

The Clayton copula is characterized by the dependence parameter α. Low values of α indicate
high levels of dependence. This is also captured by the tail dependence coefficient, whichmeasures the probability of joint extreme exceedances. For the Clayton, extreme joint crashes
receive positive probability, while joint extreme upward shocks have zero probability.
P
We specify α = ft−1 and define S(α) = di=1 u−α
i . The Clayton copula has pdf
c(u1 , . . . , ud) = (1 − d + S(α))−1/α−d ·

d−1
Y
i=0


(1 − i · α)u−α−1
.
i

(48)

We obtain the score vector∇t = −

d−1 
X
i=0


1
i
− ln(ui) + 2 ln (1 − d + S(α)) +
1−i·α
α

(49)

 Pd
−α
1
i=1 ui ln(ui )+d.
α
1 − d + S(α)
The principal difficulty for some GAS-based dynamic copula models is deriving a closed-form

expression for the information matrix. Even for simple copula models, this may quickly become
unmanageable analytically. This certainly holds for mixtures of copulas that we consider next.
To solve this analytical issue, we compute the information matrix numerically. In our current
example, the information matrix can be written as
It−1 = Et−1 (∇t )2 ≡ h(ft−1 ),(50)with the score vector ∇t as defined in (49). Note that the function h(·) in (50) does not dependon time or on any parameter other than ft−1 . We can therefore construct a grid of values
f (0) < . . . < f (n) and compute the function value h(f (j) ) at each of the grid points. Values29

at intermediate points can be obtained by cubic spline interpolation or non-parametric kernelsmoothing to ensure continuity of first and second derivatives of the likelihood function. The
numerical procedure is then as follows. First, choose starting values of the parameter θ andset the starting value f0 . Using interpolation, compute h(f0 ) and use it to scale the score steps1 = ∇1 /h(f0 ). Compute the new parameter value f1 through the GAS recursion, and again
use interpolation to obtain h(f1 ). This process is repeated for the complete sample. Finally,
the likelihood can be computed.
The disadvantage of the Clayton copula is that it only has lower tail dependence and no
upper tail dependence. Therefore, it is useful to use a symmetrized version of the Clayton
copula that allows for non-zero, but different upper and lower tail dependence. The symmetrized
Clayton copula is a mixture of the Clayton and the survival Clayton copula. Consider a general
mixture of copulas,C(u1 , . . . , ud) =mXpi Ci (u1 , . . . , ud).

(51)

i=1

Define wi = pi Ci /Pm

j=1 pj Cj as the weight of copula i.

some parameter vector θ,andIt is straightforward to derive that form∂ ln ci∂ ln c X
=
wi ·
,
∂θ
∂θ
i=1


 2
m
∂ 2 ln c X∂ ln ci ∂ ln ci ∂ ln ci
−
=
wi ·
+
∂θ∂θ′
∂θ∂θ′
∂θ ∂θ′
i=1m
X∂ ln ci
wi ·
∂θ
i=1

(52)

!

m
X

∂ ln ci
wi ·
∂θ
i=1

!′

,

(53)

and thus,
! m
!′ #
" m

 2
X
X
∂ ln ci∂ ln ci∂ ln cwi ·,= −Et−1
wi ·
Et−1
∂θ∂θ′
∂θ
∂θ
i=1
i=1

such that the scores and information matrices of the individual copulas can be used directly tobuild the driving mechanism of the mixture copula.
The GAS mechanism for the mixture of copulas has an intuitive interpretation. A given
observation may have a contribution to the evolution of either αL or αU , i.e., to either theupper or lower tail dependence. The contributions are measured in terms of the likelihood ofeach mixture vis-a-vis the total likelihood. As a result, observations that cluster in the upper
tail automatically will contribute to the evolution of αU , and similarly in the lower tail for αL .
By contrast, Patton’s methodology for the symmetrized copula cannot make automatic use ofsuch features, as its driving mechanism is given by averages of |uit − ujt| for both upper andlower tail dependence.
To illustrate the differences between these two, we construct a simulated example. We
generate data from the symmetrized Clayton copula. The lower tail dependence coefficient
30

12

GAS
Patton(10)Patton(1)
true

12

10

10

8

8

6

6

4

4

2

2

0

500

1000

1500

2000

2500

0

GAS
Patton(10)500Patton(1)
true

1000

1500

2000

2500

Figure 7: Comparisons between estimates fron the GAS framework and the Patton model for the symmetrized
Clayton copula based on a simulated data set.

follows a sinusoid pattern. The pattern of the upper tail dependence is also sinusoid, but witha period that is half as long. This makes it difficult for a model with a uniform observation
driving mechanism to capture both upper and lower tail dependence within a single model. We
plot the results in Figure 7 for smoothing parameter values m = 1 and m = 10.
As expected, the driving mechanism based only on averages of |uit − ujt | does capture someof the variation in the dependence coefficients. However, as the same mechanism underlies
both types of dependence, it has difficulty in capturing the upper and lower tail dependence
dynamics simultaneously. The GAS specification on the other hand is more successful in picking
up both dynamics. The estimate is noisier than that of the Patton model, but follows the true
dependence pattern much more closely, resulting in a significant increase in the likelihood.

4.6

Time-varying higher order moments

Following the success of the GARCH(1, 1) model, many authors have suggested further generalizations, in particular to the model with Student t errors. Hansen (1994) proposed making the31degrees of freedom time-varying, while Harvey and Siddique (1999), Jondeau and Rockinger(2003) and Brooks et al. (2005) consider models with time-varying skewness and kurtosis. We
develop a t-GAS(1, 1) model for yt = σt−1 εt where εt ∼ tνt . The error term is scaled to have2
unit variance such that σt−1
is the conditional variance while νt is a time-varying degrees ofon−1) where the latter factor is
freedom. Define the vector of factors as ft = (σt2 , − ln νb−at −a

the inverse of the logit transformation which is used to keep νt in the interval [a, b]. In ourempirical work, we select the interval [2.01, 30] to ensure that the conditional variance exists,
i.e. νt > 2. We note that it is possible to select the conditional kurtosis as a factor instead ofνt but for some time series the conditional kurtosis may not exist.
Taking derivatives of the observation density w.r.t. σt2 and νt , we obtain the score vector
−1



yt2
yt2
(νt +1)1,− 2σ2 + 2
1 + (νt −2)σ2(νt −2)σt4
t
,


t
−1
∇t =  
2
yt2
yt
yt2
(νt +1)1
1
1′ νt +1′ νt
1 + (νt −2)σ4Γ ( 2 ) − Γ ( 2 ) − 2νt − 2 ln 1 + (νt −2)σ2 + 2.
2(ν −2)2 σ4
t

t

t

t

and with some additional derivations the elements of the information matrix are


3
νt
−
−
2σt4 (νt +3)
2σt2 (νt +1)(νt +3)(νt −2) .
Et−1 [∇t ∇′t ] = 
 ′′ νt +1 (νt +4)(νt −3)
1
3
ν
′′
t
− 2σ2 (νt +1)(ν+3)(νt −2) 4 Γ
−
Γ
+
2
2
2
2(ν −2) (ν +1)(ν +3)
t

t

t

t

The functions Γ′ and Γ′′ are the digamma and trigamma functions, which are standard in any
matrix programming software. Given the results above and the derivatives of the logit transformation, it is straightforward to construct a GAS(1,1) recursion using the reparameterizationargument from (20). We label this model the tv-t-GAS(1,1) model.
We consider daily returns on the S&P 500 from February 1989 through April 2008 as anillustration. We compare the tv-t-GAS(1,1) model described above to a t-GAS(1,1) model withconstant ν, that is equation (7), and a standard t-GARCH(1,1) model with constant ν as inBollerslev (1987). Parameter estimates from each of these models are reported in Table 4.6 andestimates of the conditional variance are plotted in panel (i) of Figure 8. Focusing momentarilyon the t-GAS(1,1) model versus the t-GARCH(1,1) model, we can see that the log-likelihood
values are close to one another with the GARCH(1,1) slightly greater. Both the persistence
parameter b11 and degrees of freedom are estimated to be larger for the t-GAS(1,1) model thanfor the t-GARCH(1,1) model. Estimates of the conditional variance in panel (i) are hard to
distinguish from one another with the exception of those periods when there are outliers. To
see this more clearly, we also plot the differences between the estimates from the two GAS
models minus the GARCH model in panel (ii) of Figure 8. In the first half of the sample
before 1998, the level of volatility is lower and there are several outliers in the series. The
estimated conditional variance from the t-GARCH(1,1) model is larger than from both GAS
models. These are the large negative values in panel (ii). The difference in estimated degrees
32(i) tv−t−GAS
t−GARCH

6

1.0

t−GAS

(ii) t.v. GAS minus GARCH

constant GAS minus GARCH

0.5
4

0.0

2

−0.5

1990
30

1995

2000

2005

1990

1995

2000

2005

(iii) tv−t−GAS d.o.f.

7
20

6

5

10

(iv) BBHP(2005) d.o.f,

1990

1995

2000

2005

1990

1995

2000

2005

Figure 8: Time-varying degrees of freedom example: (i) estimated conditional variances from the t-GAS(1,1),t-GARCH(1,1), and tv-t-GAS(1,1) models; (ii) differences between the two GAS(1,1) models and the tGARCH(1,1) model; (iii) estimated time-varying degrees of freedom from the tv-t-GAS(1,1) model; (iv) estimated time-varying degrees of freedom from the GARCH model of Brooks et al. (2005).

of freedom is due to the fact that the t-GAS model does not treat outliers like a standard
t-GARCH model. From 1998-2003, volatility increases and, relative to this level, large returns
are not outliers. Estimates of the conditional variance from the GAS and GARCH models arestill significantly different and economically meaningful during this period.
Turning our attention to the tv-t-GAS(1,1) model, the estimated time-varying degrees of
freedom from this model is plotted in panel (iii) of Figure 8 and these estimates demonstratesignificant variability. The log-likelihood for our new time-varying GAS model increases appreciably relative to the t-GAS(1,1) model. Estimates of the conditional variance in panels (i) and(ii) are reasonably similar to the t-GAS(1,1) model with some differences in 1998-2004 when
the time-varying degrees of freedom increases. We compare this model with the time-varying
higher-order GARCH model of Brooks et al. (2005), which we label as the tv-t-GARCH(1,1)
model. In their model, the conditional kurtosis evolves independently from σt2 according to its
own GARCH(1,1) recursion. The implied estimates of νt can be calculated straightforwardly.

33

Table 5: Estimates from the t-GARCH(1,1), t-GAS(1,1), and tv-t-GAS(1,1) models applied todaily returns of the S&P500 from Feb. 1989 - April 2008. The tv-t-GARCH(1,1) model isfrom Brooks et. al. (2005). The full sample results are on the left. Split sample results for the
t-GAS(1,1) model are on the right.
tv-t-GAS t-GARCH

t-GAS tv-t-GARCH

t-GAS

t-GAS

pre-1998 post-1998
ω1
ω2

0.006

0.003

0.004

0.003

0.002

0.007

(0.005)(0.001)(0.001)(0.001)(0.001)(0.003)-2.373

-

-

-

-

0.057

0.047

0.044

0.049

0.026

0.061

(0.007)(0.007)(0.006)(0.007)(0.006)(0.009)-0.128

-

-

-

-

-

-

-

-

-

-

-

-

0.005

-

-

(0.310)a11
a12(0.043)a21-0.219(0.033)a22-1.498(0.002)b11
b12(0.006)

0.994

0.951

0.997

0.949

0.997

0.995

(0.003)(0.007)(0.002)(0.007)(0.003)(0.004)

0.000

-

-

-

-

-

-

-

-

-

-

-

-

0.965

-

-

5.367

10.96

(0.610)(2.074)-2359.55-3778.63(0.000)b21

0.982(0.154)b22

0.026(0.121)ν
log-like-6138.18(0.024)

6.699

7.032

(0.622)(0.677)-6153.02 -6156.46-6153.44

It is notable result that the estimates of νt from our model shown in panel (iii) are significantly different than the implied estimates of νt from the tv-t-GARCH(1,1) model of Brooks
et al. (2005). In the literature on time-varying higher-order moments, researchers typicallyforce the factors to evolve independently by imposing zero restrictions on b12 and b21 . The
34

estimated autoregressive coefficients b21 and b22 reported in Table 4.6 for the GAS model imply
that both σt2 and νt are driven by the same factor because b22 is close to zero. Accordingly, theestimates of νt in panel (iii) exhibit a similar pattern with the conditional variance in panel (i).
Estimates of νt from the tv-t-GARCH(1,1) model, which imposes these restrictions, result in adifferent behavior for the time-varying degrees of freedom. The parameter b22 is estimated tobe significant and persistent in this model.
To investigate this result further, we split the sample in half before and after 1998 and estimated ν using the t-GAS(1,1) model with constant degrees of freedom on the two sub-samples.
Estimates from this model on the two sub-samples are reported in the right-hand columns ofTable 4.6. The degrees of freedom parameter and its standard error clearly increase in the second half of the sample. Estimates of ν on the two sub-samples from the t-GARCH(1,1) model(not reported) are similar. Although this result may seem counterintuitive initially, the reason
is that large returns during this period are no longer extreme outliers because the conditional
volatility σt2 is higher. This provides some support for estimates of νt from our model and someevidence that modelling higher-order moments independently of the conditional variance may
be inappropriate. The models described in this section might be improved further by linking
the time-varying behavior of the degrees of freedom with a time-varying level parameter ωt in
the variance. We leave this extension to future research.

4.7

Time-varying multinomial

Trade by trade financial transaction prices lie on a discrete grid with most price changes taking
only a small number of values. Russell and Engle (2005) proposed modelling this behavior using
a conditional multinomial distribution with time-varying probabilities in conjunction with their
ACD model. Here, we construct a GAS version of their model. Consider the case where theobserved series yt , for t = 1, . . . , n, has a J-dimensional multinomial distribution with vector
of probabilities πt and let πj,t be the j − th elementof this vector. The vector of factors ft
will have dimension J-1 with elements fjt = ln

1−

πjt
PJ −1

j=1 πjt

where the final probability πJ,t isdetermined by the constraint that they sum to one. Next, denote ỹt and π̃t as the corresponding
J-1 dimensional vectors with the J-th elemented omitted. The score with respect to fj,t−1 is∇jt = ỹjt − π̃j,t−1 ,(54)while the diagonal and off-diagonal elements of the information matrix are
Iii,t−1 = π̃i,t−1 (1 − π̃i,t−1 ),(55)

Iij,t−1 = −π̃i,t−1 π̃j,t−1 .

(56)

35

Combining these results, a GAS(p, q) model for the multinomial distribution reduces to
ft = ω +q−1
X

Ai St−i−1 (ỹt−i − π̃t−i−1 ) +p
XBi ft−j ,

(57)

j=1

i=0

−1
where the scale matrix St−1 = It−1
can be constructed from (55) and (56). The ACM model of

Russell and Engle (2005) can be obtained as a special case of the GAS model (57) by selecting
the scale matrix St−1 to be the identity matrix. These authors also add the expected durations
from an ACD model as explanatory variables in (57).
As an empirical illustration, we use transaction data from the NYSE TAQ database on
Royal Dutch Shell A (RDSA) for the month of November 2007. After retaining trades between
9:30 and 4:00, there are 61,690 trades remaining. Panels (i)-(ii) of Figure 9 contain the observed
price changes and observed durations for the first 23,500 trades, while panel (iii) is a histogramof all the trades. The observed durations give evidence of diurnal patterns that are typicalof transactions data. In addition, the observed price changes indicate that the probabilities
should contain a similar diurnal pattern, as trades with large tick sizes are less likely during
openning and closing of the market when volume is higher.
In our sample, 98% of the price changes fall within a ± 5 tick range of zero (see panel(iii)), where a tick is now 1 cent after decimilization of the market in 2001. Decimilization
unfortunately causes an increase in the required dimension of the factor ft and a corresponding
increase in the number of parameters to estimate. For this example, ft will have a minimumof 10 dimensions meaning that the A0 matrix in an ACM(1,1) model will have 100 parameters.
Our solution to this problem is to define new factors f˜t as ft = d + Z f˜t where f˜t has dim(f˜t ) <<dim(ft ). The GAS(1,1) model reduces to−1
f˜t = A0 Z ′ It−1
ZZ ′ (ỹt − π̃t−1 ) + B1 f˜t−1 ,(58)

where the matrix Z must be restricted to identify the model. For illustration purposes, weselected dim(f˜t ) = 3 and set the upper 3 × 3 elements of Z equal to the identity matrix for
identification. Following Russell and Engle (2005), we include expected durations in (58) andjointly estimate the ACD model. We also restrict the matrices Bj to be diagonal. Specifying
a multinomial-GAS(1,2)-ACD(1,2) model for this series, some of the estimated time-varying
probabilities for the first third of the data set are shown in panels (v) and (vi) of Figure 9.
Panel (v) is a plot of the probability of a price increase of 5 ticks or more while panel (vi)plots the probability of no price movement. The model picks up the diurnal dynamics of theprice changes reasonably well with the probability of an increase of 5 ticks or more changing
considerably throughout the day. An alternative observation driven model for trade-by-trade
data has been proposed by Rydberg and Shephard (2003) using the GLAR methodology of
36

0.1

300

(i) observed price changes(ii) observed durations

200
0.0
100

0

5000

10000

15000

20000

0
40

100 (iii) histogram of price changes

5000

10000

15000

20000

(iv) expected durations − ACD model

30
20

50
37

10
−0.10

−0.05

0.00

0.05

0.10

(v) probability up 5 ticks or more

0
0.5

0.10

5000

10000

15000

20000

15000

20000

(vi) probability no price change

0.4
0.05

0

0.3
5000

10000

15000

20000

0

5000

10000

Figure 9: Time-varying multinomial GAS(1,2)-ACD(1,2) model. (i) observed price changes; (ii) observed durations; (iii) histogram of price changes; (iv) estimatedexpected duration from the ACD model; (v) estimated probability of an increase of 5 ticks or more; (vi) estimated probability of a trade with no change in price.

Shephard (1995). We note that a GAS version of their model will be slightly different but close
to their specification.

4.8

Dynamic mixtures of models

The GAS specification can provide a mixture framework for probabilities of several competing,possibly, time-varying models. Assume we have a mixture model with J components where
each component or sub-model has a likelihood Ljt. Define the vector of GAS factors as thetime-varying mixture probabilities πjt , which defines a new mixture model
Lt =

J
X

πjt Ljt.

(59)

j=1

We parameterize the πjt ’s using the logit transformation to ensure that the probabilities remainin the zero-one interval. The GAS factors are
!
J−1
X
efi
πkt ,
⇔ fjt = ln(πjt ) − ln 1 −
πjt =
P
fk
1 + J−1
k=1 e
k=1

(60)

for j = 1, . . . , J − 1 with the probability of the last component determined by the constraintP
πJt = 1 − J−1
k=1 πkt . Taking the derivative of the log-likelihood with respect to fj,t−1 , we obtain
the elements of the score vector

πj,t−1 Ljt∂Lt− πj,t−1 ,= PJ∂fj,t−1k=1 πk,t−1 Lkt(61)

for j = 1, . . . , J − 1. The interpretation of (61) is intuitive. The probability of model j isincreased if the relative likelihood of model j is above its expectation πj,t−1 . Otherwise, it isdecreased. The information matrix for this GAS model is not easy to compute analytically.
In our empirical example below, we use a mixture of two normal densities φj (y) for j = 1, 2
implying an information matrix of the form
Et−1 [∇t ∇′t ] = π1,t (1 − π1,t )Et−1"

φ1 (y) − φ2 (y)π1,t φ1 (y) + (1 − π1,t )φ2 (y)

2 #

,

where the expectation is taken with respect to the mixture distribution. We use numericalintegration to compute the information matrix, which is feasible when the mixture model (59)contains say J = 5 components or fewer.
To illustrate the methodology, we consider a time series of quarterly log U.S. real GDP
growth rates from 1947(2) to 2008(2) obtained from the Federal Reserve Bank of St. Louis.
The GAS model is a mixture of two normals with different means µi for i = 1, 2 and a common
variance σ 2 . The GAS factor is the probability that the data comes from the normal distribution
38

with low mean indicating the probability of a recession. The initial conditions of the factors fk,0are estimated as part of the static parameter vector θ, while the GAS(1,1) updating equationis adopted with an information smoothed scaling matrix St as in (14) with α = 0.05. ThisGAS model provides an observation driven alternative to a hidden Markov model (HMM). We
compare it to a simplied version of Hamilton’s (89) model with no autoregressive dynamics,
that isy t = µt + εt ,εt ∼ N (0, σ 2), µif St = 0
1
µt =
 µ2
if St = 1

pij = P (St = j|St−1 = i),i = 0, 1

j = 0, 1

In this model, the latent variable St is a regime-switching variable indicating whether the
economy is in a recession or expansion. We base our comparison on the one-step ahead predicted
estimates produced by the hidden Markov model because the GAS factor is effectively a one-stepahead predictor.
Table 6: Estimates from the GAS(1,1) mixture and hidden Markov models applied to U.S. log
real gdp growth rates from 1947(2) to 2008(2). Standard errors are in parenthesis.
µ1
GAS

µ2

σ

ω

A

B

log-like

0.203 1.129 0.868 0.364 2.415 0.673 -329.45(0.008)(0.005)(0.003)(0.018)(0.121)(0.006)

µ1

µ2

σ

p11

p22

-

HMM -0.090 1.106 0.830 0.741 0.918(0.019)(0.007)(0.003)(0.007)-333.17(0.003)

Estimates of the parameters of both models are reported in Table 4.8. The estimated values
for each mean are reasonably close with the recession parameter µ1 for the HMM model beingslightly smaller and negative. Panel (i) of Figure 10 presents the growth rate of log U.S. real
GDP along with the estimated conditional mean πt µ1 + (1 − πt )µ2 from the GAS and HMM
models. The GAS and HMM estimates nicely follow the changes in the mean of the series.
The estimated probabilities of a recession from each model are plotted in panel (ii) of Figure
10. The estimated probabilities from the GAS model reflect the possibility of the model torapidly adapt to new signals concerning the current behavior of the time series. As a result, we
obtain a clear division of regimes (switches) over time as depicted in the graph. In contrast, theone-step ahead predicted probabilities produced by the hidden Markov model do not change39

4

hidden Markov model mean
rgdp growth rate

GAS model mean

2

0

−2
1950

1955

1960

1965

1970

1975

1980

1985

1990

1995

2000

2005

1975

1980

1985

1990

1995

2000

2005

1.00
hidden Markov modelGAS model

0.75
0.50
0.25

1950

1955

1960

1965

1970

Figure 10: Mixture model example: (i) growth rate of log U.S. real GDP from 1947(2)-2008(2) and theestimated conditional mean from the GAS(1,1) model and the hidden Markov model; (ii) one-step ahead predictedprobability of a recession from each model. NBER recession dates are represented by the shaded regions.

as rapidly and are not as clear. The GAS model offers a convenient method for forecasting
economic downturns. A multivariate model incorporating leading economic variables would be
an interesting extension of the GAS model presented here.

5

Simulation experiments

In this section, we provide simulation evidence on the statistical properties of the GAS ML
estimators for a selected set of three examples. We concentrate on the marked point processmodel, the linear state space model with time varying variances, and the Gaussian copula with
time-varying correlation.

5.1

The marked point process model

To investigate the statistical properties of the GAS model for the marked point processes ofSection 4.2, we consider a simplified version of this model. We consider a cross-section of firms
40with two possible ratings, A and B, and possible transitions between them. Neither of thestates is absorbing, such that there is no attrition of the panel of firms over time. We consider
a panel of N = 250 and N = 2, 500 firms. As the results are comparable, we only present thegraphs for N = 2, 500.
The relevant equations are
λAt = η1 + ft ,λBt = η2 + αft ,ft = Ast + Bft−1 ,with η1 = −3.5, η2 = −4.0, α = −1, A = 0.025, B = 0.95, and where λAt and λBt denote theintensities of an A firm becoming a B firm, and of a B firm becoming an A firm, respectively.
The parameter values are roughly in line with the empirical estimates concerning the levels ofintensities and the magnitude of the systematic factor.
We consider three sample sizes for the time series dimension, namely T = 20, 50, 100. We
use 1,000 simulations of the true dgp. For each simulation, we compute the ML estimate aswell as its t-value based on the numerical second derivative of the likelihood at the optimum.
As in the empirical application, we enforce stationarity by parameterizing and estimating thelogit transform of B in the GAS equation.
The results are presented in Figure 11 and are very encouraging. The densities of theparameter estimates reveal that for increasing sample sizes T , the estimates peak more at theirtrue values. There is some skewness in the densities of α and A, particularly for smaller sample
sizes. If we consider the t-values, however, it appears that the approximation by the normal
distribution for purposes of inference is reasonable, even for sample sizes as small as T = 20.

5.2

State space models with time-varying variances

The finite sample properties are investigated for the UC model with time-varying variances,
see Section 4.4 where the model is considered as an illustration. In particular, the local level
model (32) is adopted with the log-variances of the irregular and level disturbances, ln σε2 andln σξ2 respectively, treated as GAS factors. The model for the two log-variance factors is given
by ft = ω + Ast−1 + Bft−1 where ft is a 2 × 1 vector and the 2 × 2 matrices A and B are both
diagonal. The true parameter values in the Monte Carlo study below are chosen to be close to
those obtained in the illustration of Section 4.4. To enforce a stable behaviour of the factors
over time, the parameters are subject to transformations. The true values for the resulting
coefficients are given byln ω1 = −4.5,

ln ω2 = 5.0,

logit A11 = 0.1,
41logit A22 = 1.0,

η1

25
20

η2

25

T=20
T=50
T=100

20

15

α

T=20
T=50
T=100

2.0

0.4

1.0

0.4

T=20
T=50
T=100
normal−5Density
0.4

0.5

25

0.5

−4.1 −4 −3.9
Density

T=20
T=50
T=100
normal00.025 0.075
Density
0.4

T=20
T=50
T=100
normal

0.0 2.5 5.0
Density

T=20
T=50
T=100
normal

0.4

0.3

0.3

0.2

0.2

0.2

0.2

0.2

0.1

0.1

0.1

0.1

0.1

−5

0

5

T=20
T=50
T=100

50

5
−3.6 −3.5 −3.4
Density

B

1.0

10

5

75

1.5

T=20
T=50
T=100

1.5

15

10

A
T=20
T=50
T=100

−5

0.3

0.3

0

5

−5

0

5

−5

T=20
T=50
T=100
normal

0.3

0

5

−5

0

5

Figure 11: Simulation densities over 1,000 simulations of a marked point process model. The top panel
contains the densities of the parameter estimates, the bottom panels contain the densities of t-values computed
using the inverted second derivative of the Hessian at the optimum.

logit B11 = 3.5,

logit B22 = 3.2,

logit α = −1.0,

where α is the parameter of the EWMA smoothing scheme for the information matrix, see thediscussion in Section 2.1. The subscript refers to the element of the associated vector or matrix.
Given an observed time series, all coefficients can be estimated using the methods describedearlier, see Section 4.4 for further the details.
The Monte Carlo design is similar to the one presented in the previous section. The numberof simulations equals 1, 000 while the time series dimension is set to T = 100, 250 and 1000.
Smaller values for T are less interesting for this model. The results of the Monte Carlo study
are given in Figure 12.

42

0.3

ω1

0.3

T=100
T=250
T=1000

A1

0.3

T=100
T=250
T=1000

B1

0.4

T=100
T=250
T=1000

ω2
T=100
T=250
T=1000

0.3
0.2

0.2

0.2
0.2

0.1

0.1

0.1
0.1

−5

0

A2
0.3

T=100
T=250
T=1000−5.0 −2.5
B2
0.75

0.0

2.5

T=100
T=250
T=1000

0
α
0.5
0.4

0.50

5

10

−10.0

−7.5

−5.0

T=100
T=250
T=1000

0.3

0.2

0.2
0.25

0.1

0.1
−5

0

5

0

2

4

−2.5

0.0

2.5

Figure 12: Simulation densities over 1,000 simulations of a local level model with time-varying log variances.
The top panel contains the densities of the parameter estimates, the bottom panels contain the densities oft-values computed using the inverted second derivative of the Hessian at the optimum.

5.3

Time-varying Gaussian copula model

In this section, we inspect the finite sample properties of the time-varying Gaussian copula
model described in Section 4.5. The model is specified as
xt − 2A+ Bft−1 .
· yt − ρt−1 − ρt−1
ft = ω +(1 − ρt−1 )(1 + ρt−1 )1 + ρ2t−1
where xt = Φ−1 (u1t )2 + Φ−1 (u2t )2 and yt = Φ−1 (u1t )Φ−1 (u2t ). The parameter settings forthe true dgp were ω = 0.02, A = 0.15, and B = 0.96, while the simulation sample sizes
were T = 200, 400, 600 observations. To ensure stationarity of the factor ft and for numericalstability, we restricted A and B during estimation using a logit transformation.
The results from the experiment over 1,000 simulations are shown in Figure 13. The densityof the parameter estimates are converging toward their true values as T increases. The rateof convergence appears to be slower for this model than the marked point process model.
The densities of the t-values appear slightly biased for the ω and B parameters with the bias
diminishing as the sample size increases.
43

ω

A
0.75

T = 200
T = 400
T = 600

20
15

B
0.6

T = 200
T = 400
T = 600

0.50

0.4

0.25

0.2

T = 200
T = 400
T = 60010
5−0.1
Density
0.5

0.0

0.1

−5.0
DensityT = 200
T = 400
T = 600
normal

0.4

−2.5

0.0

0.0

2.5
Density

0.4

T = 200
T = 400
T = 600
normal

0.4

5.0

7.5

T = 200
T = 400
T = 600
normal

0.3

0.3
0.2
0.2

0.2
0.1

0.1
−5

0

5

−5

0

5

−5

0

5

Figure 13: Simulation densities over 1,000 simulations of a time-varying Gaussian copula model. The top
panel contains the densities of the parameter estimates, the bottom panels contain the densities of t-values
computed using the inverted second derivative of the Hessian at the optimum.

6

Conclusions

In this paper, we have introduced Generalized Autoregressive Score (GAS) models. A GASmodel is a uniformly applicable observation driven model specification to capture time variationin parameters. We have shown how GAS models encompass other well-known models, such as
generalized autoregressive conditional heteroskedasticty models and autoregressive conditional
duration and intensity models as well as multiplicative error models and single source of error
models. The advantage of the GAS model is that it exploits the full likelihood information.
By making a scaled (local density) score step, the time-varying parameter automatically tries
to reduce its prediction error at the current observation with respect to current parametervalues. Though based on a completely different paradigm, the GAS model provides a powerfuland highly competitive alternative to other observation driven models as well as parameter
driven models. We have illustrated this extensively by describing a number of non-trivial
empirical and simulated examples. Some of these examples are interesting in their own rightand provide interesting extensions or alternative specifications for parameter driven models with44

time-varying parameters, in particular for state space models with stochastically time-varying
parameters, for multivariate marked point processes, and for time-varying copula models.
There are many interesting future research directions. The issues of identification, consistency, stationarity, and asymptotic distribution theory require more work than presented here.
Due to its generality and applicability for a wide class of models, however, it appears difficult tocome up with an uniform set of conditions for stationarity and consistency that is applicable toall situations of interest. A more promising route may be to formulate conditions for particularsub-sets of models with a GAS specification.
To investigate the finite-sample properties of GAS models in more detail is a second direction
for further research. Although we have provided a number of interesting empirical and simulated
examples, a more systematic study into the statistical properties of parameter estimates for
GAS models may be appropriate.
A third direction for future research concerns the development of misspecification tests for
GAS models. On the one hand, we require goodness-of-fit tests and model selection criteria for
GAS models. Many of such tests and diagnostics are already developed for the class of GARCH
models. On the other hand, the GAS model itself might provide a powerful basis for dynamic
misspecification tests. A similar approach to test for the presence of possible GARCH effects
is already widely applied in empirical studies. Lagrange multiplier based tests for the possible
presence of GAS effects are straightforward extensions of these. Such tests might provide a
useful empirical tool for testing for possible time variation in parameters in the context of alarge class of nonlinear and non-Gaussian models.
A fourth direction of research is the application of the GAS specification to new models.
In this paper, we have tried to review a number of interesting directions of new models with
time-varying parameters. However, the GAS framework is not restricted to these, and othernew and empirically relevant models with time-varying parameters would provide additional
support for the usefulness of GAS as an empirical modeling tool.
Long-memory versions of the GAS model would be another interesting research project.
However, the long-memory specification for GAS models is not trivial and therefore more theoretical and empirical research is needed. A related issue is that GAS models may be interpreted
as discrete time approximation of their parameter driven counterparts. An interesting research
project may be to bridge the gap between GAS models and parameter driven models in asimilar continuous time limiting sense as obtained by Nelson (1996) who has bridged the gap
between GARCH models and stochastic volatility model specifications.
A sixth direction of future research is to provide a systematic comparison of the advantages
and disadvantages of parameter driven versus observation driven models in a wider setting45

than GARCH and ACI. Given numerical advances for non-linear and non-Gaussian state spacemodels, and given the general applicability of the current GAS specification, such comparisons
have become feasible.
Finally, there are various computational details that need to be studied in further detail.
Three issues of particular interest are: finding starting values, finding the required degree of
information smoothing for the GAS updating step in particular models, and finding better
numerical approximations to the scaling matrix if it cannot be computed analytically. With
respect to the first issue, our findings so far are mixed. In relatively straightforward models, theproblem of finding appropriate starting values does not exist. In particular, if the information
matrix is clearly non-singular for all sample observations, the maximum likelihood maximization
algorithm converges quickly and robustly. Introducing information smoothing as well as finding
reasonable starting values become more relevant when an observation contains limited or no
information on the parameter of interest. This is particularly relevant if there are regions witha degenerate information matrix. In our experience, some degree of information smoothing is
indispensable in such cases. In addition, automatic smoothing by estimating the smoothing parameter directly from the data has increased the likelihood value in several cases. In our current
implementations, however, the information smoothing is rather rigid. One could consider more
involved specifications, where the degree of smoothing also depends on the current position in
the sample and the parameter space. The third issue concerns further progress that is needed
for models where the information matrix cannot be computed analytically. In the illustrationof time-varying copulas in Subsection 4.5, we provided some suggestions based on numericalinterpolation techniques using kernel smoothing in low-dimensional parameter spaces. Further
extensions are needed to develop computationally feasible estimation methods for GAS models
with large parameter spaces and (perhaps even) more complicated specifications.

References
Bauwens, L. and N. Hautsch (2006). Stochastic conditional intensity processes. Journal of Financial Econometrics 4 (3), 450.
Bauwens, L., S. Laurent, and J. Rombouts (2006). Multivariate GARCH models: a survey. Journal of Applied
Econometrics 21 (1), 79–109.
Benjamin, M. A., R. A. Rigby, and M. Stanispoulos (2003). Generalized autoregressive moving averagemodels. Journal of the American Statistical Association 98 (461), 214–223.
Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31,
307–327.
Bollerslev, T. (1987). A conditionally heteroskedastic time series model for speculative prices and rates ofreturn. The Review of Economics and Statistics 69 (3), 542–547.

46

Bollerslev, T., R. F. Engle, and D. B. Nelson (1994). ARCH models. Handbook of Econometrics 4, 2959–3038.
Brooks, C., S. P. Burke, S. Heravi, and G. Persand (2005). Autoregressive conditional kurtosis. Journal of
Financial Econometrics 3 (3), 399–421.
Cipollini, F., R. F. Engle, and G. M. Gallo (2006). Vector multiplicative error models: representation andinference. NBER Working Paper .
Clark, P. K. (1989). Trend reversion in real output and unemployment. Journal of Econometrics 40, 15–32.
Cox, D. R. (1958). The regression analysis of binary sequences (with discussion). Journal of the Royal Statistical Society, Series B 20 (2), 215–242.
Cox, D. R. (1981). Statistical analysis of time series: some recent developments. Scandinavian Journal of
Statistics 8, 93–115.
Davis, R. A., W. T. M. Dunsmuir, and S. Streett (2003). Observation driven models for Poisson counts.
Biometrika 90, 777–790.
Dias, A. and P. Embrechts (2004). Dynamic copula models for multivariate high-frequency data in finance.
Manuscript, ETH Zurich.
Diebold, F. X. and C. Li (2006). Forecasting the term structure of government bond yields. Journal of
Econometrics 130, 337–364.
Diebold, F. X., G. Rudebusch, and S. Aruoba (2006). The macroeconomy and the yield curve. Journal of
Econometrics 131, 309–338.
Doucet, A., N. de Freitas, and N. Gordon (2001). Sequential Monte Carlo Methods in Practice. New York,
NY: Springer-Verlag.
Duffie, D., A. Eckner, G. Horel, and L. Saita (2006). Frailty correlated default. manuscript Stanford University.
Durbin, J. and S. J. Koopman (2001). Time Series Analysis by State Space Methods. Oxford: Oxford University Press.
Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United
Kingdom inflation. Econometrica 50 (4), 987–1007.
Engle, R. F. (2001). GARCH 101: The use of ARCH/GARCH models in applied econometrics. Journal of
Economic Perspectives 15 (4), 157–168.
Engle, R. F. (2002a). Dynamic conditional correlation: a simple class of multivariate generalized autoregressive conditional heteroskedasticity models. Journal of Business & Economic Statistics 20 (3), 339–350.
Engle, R. F. (2002b). New frontiers for ARCH models. Journal of Applied Econometrics 17, 425–446.
Engle, R. F. and T. Bollerslev (1986). Modelling the persistence of conditional variances. Econometric Reviews 5 (1), 1–50.
Engle, R. F. and G. M. Gallo (2006). A multiple indicators model for volatility using intra-daily data. Journalof Econometrics 131, 3–27.
Engle, R. F. and S. Manganelli (2004). CAViaR: conditional autoregressive value at risk by regression quantiles. Journal of Business & Economic Statistics 22 (4), 367–382.
Engle, R. F. and J. R. Russell (1998). Autoregressive conditional duration: a new model for irregularly spacedtransaction data. Econometrica 66 (5), 1127–1162.

47

Fiorentini, G., G. Calzolari, and L. Panattoni (1996). Analytic derivatives and the computation of GARCH
estimates. Journal of Applied Econometrics 11, 399–417.
Ghysels, E., A. C. Harvey, and E. Renault (1996). Stochastic volatility. Handbook of Statistics 14, 119–191.
Hamilton, J. (1989). A new approach to the economic analysis of nonstationary time series and the businesscycle. Econometrica 57 (2), 357–384.
Hansen, B. E. (1994). Autoregressive conditional density estimation. International Economic Review 35 (3),705–730.
Harvey, A. C. (1989). Forecasting, structural time series models and the Kalman filter. Cambridge, UK:
Cambridge University Press.
Harvey, A. C. and A. Jaeger (1993). Detrending, stylised facts and the business cycle. Journal of Applied
Econometrics 8, 231–247.
Harvey, A. C., E. Ruiz, and E. Sentana (1992). Unobserved component time series models with ARCH
disturbances. Journal of Econometrics 52 (1-2), 129–157.
Harvey, C. R. and A. Siddique (1999). Autoregressive conditional skewness. Journal of Financial and Quantitative Analysis 34, 465–487.
Hosking, J. R. M. (1981). Fractional differencing. Biometrika 68 (1), 165–176.
Jondeau, E. and M. Rockinger (2003). Conditional volatility, skewness, and kurtosis: existence, persistence,and comovements. Journal of Economic Dynamics and Control 27 (10), 1699–1737.
Koopman, S. J., A. Lucas, and A. Monteiro (2008). The multi-state latent factor intensity model for credit
rating transitions. Journal of Econometrics 142, 399–424.
Koopman, S. J., A. Lucas, and B. Schwaab (2008). Forecasting cross-sections of frailty-correlated default.
08-029/4.
Nelson, C. R. and A. Siegel (1987). Parsimonious modelling of yield curves. Journal of Business 60 (4),473–489.
Nelson, D. B. (1996). ARCH Models as diffusion approximations. Modelling Stock Market Volatility: Bridging
the Gap to Continuous Time.
Ord, J. K., A. B. Koehler, and R. D. Snyder (1997). Estimation and prediction for a class of dynamic nonlinear
statistical models. Journal of the American Statistical Association 92 (440).
Park, J. Y. and J. Qian (2006). Autoregressive modeling of time-varying densities in functional space. Working
paper Texas A&M University.
Patton, A. J. (2006). Modelling asymmetric exchange rate dependence. International Economic Review 47 (2),527–556.
Russell, J. R. (2001). Econometric modeling of multivariate irregularly-spaced high-frequency data. Unpublished manuscript, University of Chicago, Graduate School of Business.
Russell, J. R. and R. F. Engle (2005). A discrete-state continuous-time model of financial transactions prices
and times: the autoregressive conditional multinomial-autoregressive conditional duration model. Journalof Business & Economic Statistics 23 (2), 166–180.
Rydberg, T. H. and N. Shephard (2003). Dynamics of trade-by-trade price movements: decomposition andmodels. Journal of Financial Econometrics 1, 2–25.

48

Shephard, N. (1995). Generalized linear autoregressions. Unpublished manuscript, Nuffield College, Universityof Oxford .
Shephard, N. (2005). Stochastic Volatility: Selected Readings. Oxford: Oxford University Press.
Shephard, N. and M. K. Pitt (1997). Likelihood analysis of non-Gaussian measurement time series. Biometrika 84, 653–67.
Stock, J. H. and M. W. Watson (2007). Why Has U.S. Inflation Become Harder to Forecast? Journal of
Money, Credit, and Banking 39, 3–33.
van den Goorbergh, R. W. J., C. Genest, and B. J. M. Werker (2005). Bivariate option pricing using dynamic
copula models. Insurance Mathematics and Economics 37 (1), 101–114.
Watson, M. W. (1986). Univariate detrending methods and stochastic trends. Journal of Monetary Economics 18, 49–75.

49