Brave et al., A Practitioner's Guide and MATLAB Toolbox for Mixed Frequency State Space Models (2022-10)
Summary
Article by Scott A. Brave, R. Andrew Butters and David Kelley, "A Practitioner's Guide and MATLAB Toolbox for Mixed Frequency State Space Models," in the Journal of Statistical Software, October 2022, Volume 104, Issue 10. The authors show how state space methods can facilitate estimation and inference with mixed frequency data and present applications of the authors' MATLAB toolbox. Examples include quarterly gross domestic product (GDP) modeled jointly with monthly information on employment and industrial production. The article states that the opinions are those of the authors and not necessarily those of the Federal Reserve Bank of New York or the Federal Reserve System. An appendix describes how filter and smoother weights are computed and used.
Summary drafted by a model from the document's text below and checked by script against that text before publication. It is a navigation aid, not a reading of what the document proves. Where AI is used
Cite this document
“Brave et al., A Practitioner's Guide and MATLAB Toolbox for Mixed Frequency State Space Models (2022-10).”
Original at www.jstatsoft.org. Archived at Pandemic Darlings, https://pandemicdarlings.com/source-documents/brave-et-al-a-practitioner-s-guide-and-matlab-toolbox-for-mixed-freque-9a596ad59c/
, SHA-256 1dcbd2380623861d712f318a808ab3fe47bd0e834a13d75dacf6af4c46aa4122.Full text
JSS Journal of Statistical Software
October 2022, Volume 104, Issue 10. doi: 10.18637/jss.v104.i10
A Practitioner’s Guide and MATLAB Toolbox for
Mixed Frequency State Space Models
Scott A. Brave R. Andrew Butters David Kelley
Morning Consult Indiana University Federal Reserve Bank
of New York
Abstract
The use of mixed frequency data is now common in many applications, ranging from
the analysis of high frequency financial time series to large cross-sections of macroeconomic
time series. In this article, we show how state space methods can easily facilitate both
estimation and inference in these settings. After presenting a unified treatment of the state
space approach to mixed frequency data modeling, we provide a series of applications to
demonstrate how our MATLAB toolbox can make the estimation and post-processing of
these models straightforward.
Keywords: mixed frequency, state space, Kalman filter, maximum likelihood, MATLAB.
1. Introduction
The use of mixed frequency data is now common in many applications, ranging from high
frequency financial time series to large cross-sections of macroeconomic time series. For
example, a model of a security might jointly involve hourly spread data with the daily closing
price. Or, in a macroeconomic model, quarterly gross domestic product (GDP) might be
incorporated jointly with monthly information on employment and industrial production.
In this article, we show how state space methods can easily facilitate both estimation and
inference in mixed frequency settings. By preserving the frequency of the data generating
process, this approach offers unique identification and inference possibilities.
In the past, integrating data observed at differing frequencies was often considered a pre-
analysis step of data collection, either by temporally aggregating the higher frequency time
series or by narrowing the selection of time series to a single frequency. However, given the re-
sulting potential for temporal aggregation bias, the demand for methods that instead explicitly
The opinions expressed here are those of the authors and not necessarily those of the Federal Reserve Bank
of New York or the Federal Reserve System.
2 MATLAB Toolbox for Mixed Frequency State Space Models
accommodate the mixed frequency nature of a dataset has grown substantially. For example,
mixed frequency extensions to traditional dynamic factor models (Mariano and Murasawa
2003; Aruoba, Diebold, and Scotti 2009; Brave and Butters 2012), dynamic stochastic gen-
eral equilibrium models (Kim 2010; Foroni and Marcellino 2014), and vector auto-regressions
(Kuzin, Marcellino, and Schumacher 2011; Schorfheide and Song 2015; Brave, Butters, and
Justiniano 2016) have found their way into both economic and financial applications.
All of the above extensions share in common the use of state space methods in order to
accommodate the mixed frequency nature of their applications.1 Generally speaking, state
space methods are well suited to overcome the challenges posed by mixed frequency data. A
critical source of their comparative advantage centers on their flexible accommodation of han-
dling missing observations through the use of latent state variables. Somewhat fortuitously,
state space methods can treat missing observations due to the mixed frequency or staggered
release of data (another kind of missing observation) similarly — provided the appropriate
temporal restrictions have been embedded into the state space system. Paying this upfront
cost allows the researcher to model the system in the highest (or an even higher) frequency
observed in the dataset and opens up the full array of inference and estimation capabilities
of a standard state space model. Thus, it should come as no surprise that these methods are
quickly becoming a standard in the applied econometrician’s toolkit.
Despite the increased use of state space methods in mixed frequency applications, a unified,
self-contained description of these methods and the challenges faced in a mixed frequency
environment does not, to the best of our knowledge, exist in the extant literature. The same
is true for the existence of a singular software package that can readily handle setting up,
estimating, and post-estimation processing of these models with limited user cost. Instead,
practitioners are often left having to piece together many different programming languages
with far from uniform standards and default settings for key specification choices.
In this paper, we present a unified framework that enables low-frequency data to be integrated
into state space models with higher frequency data in a manner that respects the temporal
aggregation properties of each time series. To further lower the entry cost of utilizing this
approach, we also introduce and illustrate the use of our MATLAB2 (The MathWorks Inc.
2021) toolbox in applying these methods — what we call the Mixed Frequency State Space
(MFSS) toolbox. The MFSS toolbox is available from GitHub at https://github.com/
davidakelley/MFSS.
The toolbox offers similar estimation and post-processing capabilities of popular alternatives,
e.g., STAMP (Koopman, Harvey, Doornik, and Shephard 2010b); SsfPack in Ox (Koopman,
Doornik, and Shephard 2010a; Doornik 2007), KFAS (Helske 2018), statsmodels (Perktold,
Seabold, and Taylor 2017), and the Econometrics Toolbox (in MATLAB). For a more complete
overview of available alternatives, see Commandeur, Koopman, and Ooms (2011). In addition,
MFSS integrates leading computational methods, e.g., univariate treatment of multivariate
time series (Koopman and Durbin 2000), exact initialization of the state (Koopman and
1
Accommodating mixed frequency data has also been the focus of statistical procedures for temporal dis-
aggregation (Harvey and Pierse 1984; Harvey 1989; Harvey and Chung 2000; Proietti 2006). In these cases,
the object of interest is the higher frequency (latent) estimate of one of the lower frequency variables in the
model. Furthermore, within the regression framework, the mixed data sampling (MIDAS) approach provides
an alternative methodology to the state space approach in explicitly incorporating the mixed frequency nature
of time series (Ghysels, Santa-Clara, and Valkanov 2004).
2
The toolbox is currently limited to use in MATLAB. Introspection of anonymous functions in Octave is
limited, preventing necessary steps in estimating parameters.
Journal of Statistical Software 3
Durbin 2003), improved calculation of starting values for maximum likelihood estimation
(Carpenter, Gelman, Hoffman, Lee, Goodrich, Betancourt, Brubaker, Guo, Li, and Riddell
2017), a compiled implementation of the filter and smoother (providing a roughly 15× speedup
compared to native MATLAB code), and the computation of filtering and smoothing weights
(Koopman and Harvey 2003), all within an environment that can accommodate an arbitrary
mix of data frequencies — a feature unmatched by existing software packages. It should be
noted, that the breadth of these capabilities does not come without a cost, as our package is
restricted to the set of linear state space methods.3
Our goal in this paper is to outline and highlight the foundational elements of setting up and
estimating a mixed frequency state space model through a series of applications (e.g., a small
monetary policy mixed frequency vector autoregression, MF-VAR) that illustrate the utility
of the MFSS toolbox. Given that goal, we leave the reader to find most of the theoretical
derivations and evidence behind the methodological advancements of time series state space
modeling to the traditional treatments (e.g., Harvey 1989; Durbin and Koopman 2012). Each
of the applications are selected for their comparative advantage in illustrating the different
types of requirements for estimating and analyzing a mixed frequency model. In some cases,
such as the trend-cycle decomposition of GDP (Section 4.1) and the estimation of the natural
rate of interest (Section 4.5), they are the first mixed frequency estimates of their kind.
To ensure consistency throughout the article, and to avoid any confusion that could be caused
by the minor variations used in the literature, we outline the standard normal-linear state
space model in Section 2. In Section 3, we outline the key departures and extensions from
the standard state space model required to accommodate mixed frequency data, leaving some
of the underlying details not likely found in usual references to an online appendix. With
the foundational elements in place, in Section 4 we provide a series of topical applications
that benefit from the use of the mixed frequency approach (e.g., a dynamic factor model of
the business cycle). For each application, we provide code-examples and visualizations of key
results that showcase the approachable syntax of our MFSS toolbox in achieving many of
the standard tasks associated with setting up, estimating, and conducting post-estimation
analysis of a mixed frequency state space model. In Section 5, we discuss other likely uses for
the toolbox and how it could be extended in the future.
2. A normal-linear state space model
Our approach to mixed frequency data modeling centers around normal-linear state space
models of the form shown below, where data are noisy observations of an unobserved state.4
yt = Zt αt + dt + βt xt + εt εt ∼ N (0, Ht ) (1)
αt = Tt αt−1 + ct + γt wt + Rt ηt ηt ∼ N (0, Qt ) (2)
The observation equation, Equation 1, specifies the vector of data yt as a linear combination of
the unobserved (latent) states αt , a deterministic component dt , possible regressors βt xt , and
3
This restriction is made less problematic because of the many opportunities for nonlinear state space
models to be formulated and estimated through linear state space approaches (e.g., see Section 4.3), and the
full time varying parameter capabilities of our package. For a MATLAB package capable of handling nonlinear
state space models and a limited set of temporal aggregation schemes, see Villegas and Pedregal (2018).
4
For our purposes, we describe everything in the normal-linear state space framework typically used in most
empirical applications. However, much of the utility from this article applies more generally under alternative
distributional and functional forms with only modest modifications.
4 MATLAB Toolbox for Mixed Frequency State Space Models
a measurement error εt . The state, or transition, Equation 2 specifies that the unobserved
state evolves according to a univariate or vector autoregressive process with a determinis-
tic component ct , possible regressors γt wt , and shocks ηt . For an extended discussion and
examples that build intuition around these models, see Sargent and Stachurski (2018).
A few additional remarks on our specification of these models:
• Exogenous variables (βt xt and γt wt ) could be subsumed into dt and ct , but are included
to simplify specification of models with exogenous series.
• The elements of εt may be correlated with each other and the elements of ηt can be
correlated with each other, but we assume independence across the two sets of shocks
and across all time periods.
• The time indexing of system matrices, Tt , Rt , ct , γt and state equation shocks ηt are a
slight deviation from the treatment in Durbin and Koopman (2012), and aligns more
closely with the treatment in Harvey (1989).
• To allow for potentially fewer shocks than states, the Rt matrix selects which shocks
are applied to which states. However, it may have elements other than zero or one. The
Rt matrix also allows for lagged variables to be included in the state. For example, it
allows a VAR(p) model to be written in companion form in the state equation.
• Due to the recursive nature of the state equation, we must specify an initial condition,
α0 ∼ N (a0 , P0 ). These parameters may be set arbitrarily based on outside information,
implied from the steady state of the transition equation parameters, or the diagonal
elements of P0 can be arbitrarily large, in which case the values of a0 are irrelevant.
Given a set of model parameters and observations, we are interested in recovering the implied
latent state αt . The Kalman filter is a set of recursive equations that provides the minimum
mean-squared error estimate of αt at each point in time conditional on all observed data
prior to that time, i.e., at = E(αt | y1 , . . . yt−1 ), while the Kalman smoother provides the
corresponding estimate given the full sample of data, α̂t = E(αt | y1 , . . . , yT ) (see Appendix A).
Additionally, the Kalman filter recursion provides the one-step ahead forecast errors, which
can be used to compute the likelihood of a set of model parameters. While in our experience
maximum likelihood estimation for most simple models is possible with standard gradient
ascent algorithms, some models require alternative methods (see Appendix B). Furthermore,
because the Kalman filter and smoother recursions are linear, it is possible to decompose
estimates of the state according to a set of weights based on the observed data and parameters
(see Appendix C).
3. Mixed frequency state space modeling
Our approach to modeling data observed at differing frequencies can be viewed as interpolat-
ing any lower frequency time series according to a set of dynamics that are consistent across
all observed time series. In other words, we appropriately aggregate the unobserved dynamic
process captured in the transition equation of the state space model for each of these time
series before it is related to its lower frequency observations. Doing so separates the obser-
vation frequency of the data from the frequency of the model and allows the specification of
the model to be time-consistent.
Journal of Statistical Software 5
3.1. Defining accumulators and augmenting state spaces
In handling mixed frequency data, we follow the direction of Harvey (1989) and Mariano and
Murasawa (2003). The latent state process αt is determined at a base frequency that matches
the frequency of the highest-frequency observation.5 To handle data observed at differing
frequencies, we partition yt into two components. Series observed at the base frequency are
denoted ytb . Series at a lower frequency are denoted ytl . To move from the base frequency to
a lower frequency, define an aggregation as a linear function A that takes the most recent
n base-frequency values and creates a low-frequency version, ytl = A(yth , yt−1 t−n+1 ),
h , . . . , yh
where yt is an unobserved high-frequency version of yt .
h l
For the time being, we restrict the measurement equation loadings (Z and β) on time series
requiring aggregation to be time-invariant. Given the state αt and some set of parameters
{Z h , dht , β h xht , Hth }, we could create an interpolated version of the observations with yth =
Z h αt + dht + β h xht + εht . Because A is linear, we can factor this expression as
ytl = A(Z h αt + dht + β h xht + εht , . . . , Z h αt−n+1 + dht−n+1 + β h xht−n+1 + εht−n+1 )
= Z h A({αt−i }n−1
i=0 ) + A({dt−i }i=0 ) + β A({xt−i }i=0 ) + A({εt−i }i=0 )
h n−1 h h n−1 h n−1
= Z l ζt + dlt + β l xlt + εlt .
We next define an accumulator, ζt , that can be included in the state to relate the base
frequency state αt to the low-frequency observations ytl , implying that the high-frequency
loading matches the lower frequency loading, Z l = Z h . In order to incorporate a low-frequency
observation in a state space model, we need to be able to integrate it into the Kalman filter
recursions, and so require a recursive formula for an accumulator,
ζt = ft (ζt−1 , αt )
= Aζt ζt−1 + ATt αt−1 + Act + Aγt wt + AR
t ηt ,
where {Aζt , ATt , Act , Aγt , AR
t } are functions of the system matrices. Formulas are given below
for each accumulator type that we consider. Combined with a transition equation for αt and
an observation equation for the base-frequency data, a state space model that accommodates
mixed frequency observations is given by6
" # " #" # " # " #" # " #
ytb Zb 0 αt db βb 0 xbt εb
= t + tl + t + tl
ytl 0 Zl ζt dt 0 βl xtl εt
0
" # " #" # " # " # " #
αt Tt αt−1 c γ Rt
= + tc + tγ wt + ηt .
ζt ATt Aζt ζt−1 At At AR
t
This makes clear that in the mixed frequency context, examining the effect of ytb and ytl on
αt becomes more difficult. Instead of comparing the estimated elements of Z to each other,
as in a single-frequency context, we instead must examine how the data affects the estimated
states using the Kalman filter and smoother weights to find the effect of ytl on αt .
5
It is possible to specify the transition equation at a higher frequency than the data. For example, Aruoba
et al. (2009) specify a daily frequency despite the highest frequency data in their model being weekly; and in
Section 4.4 we estimate a trend-cycle decomposition of quarterly GDP at a monthly frequency.
6
To keep the state dimension as small as possible, we only need to augment the states that will be aggregated
for the time series under consideration. Specifically, for any observation requiring an accumulator, we can create
augmented states only for those where Z contains non-zero elements.
6 MATLAB Toolbox for Mixed Frequency State Space Models
Sampling type Level ∆1 ∆H≥1
Point-in-time None None Sum
Sum Sum Sum Sum
Average Average Average Triangle Average
Table 1: Sampling types and accumulator specifications.
3.2. Aggregation types and accumulator construction
In general, the matrices in the recursion for ζt are time-varying even if Tt , ct , γt , and Rt are
not. They depend on the calendar that the low-frequency data follows, a scalar function of
time t specific to each time series. The calendar definitions will also be specific to each type
of accumulator, but some common notation is helpful. Define P1 as the set of time indexes
that start a new low-frequency period. For example, with a monthly base frequency and an
observed series at a quarterly frequency, these would be the indexes of the first month of each
quarter (January, April, July, and September). Accordingly, for the observations to line up
with the calculations for ζt , the data must be placed in the last base-frequency period per
low-frequency period (March, June, September, and December).
Construction of the accumulators differs based on the nature of the aggregation of each time
series. For the purposes of this article (and the MFSS toolbox), there will exist three types.
The appropriate accumulator type for any low-frequency time series is a function of the way
in which it is sampled as well as any transformation applied to it. These dimensions are
intricately linked to a time series’ interpretation of being either a stock or flow variable.
As far as sampling is concerned, point-in-time sampling (common for “stock” variables) ap-
plies to data series which are measured at a specific sub-set of time periods constituting the
underlying base frequency (e.g., the exchange rate on business days in a daily base frequency
model). Extending the typology of Marcellino (1999), the other two types of sampling pro-
cesses include sums (e.g., total monthly job separations from weekly separations) and averages
(e.g., annual GDP from annualized quarterly GDP).
A series’ transformation, or lack thereof, also affects its temporal aggregation properties. For
example, while the (untransformed) level of payroll employment constitutes a point-in-time
sampling process, its one-period change (∆1 ) reflects a sum sampling process. Finally, if
the series’ transformation constitutes a “long-difference” of length H greater than the base
frequency (∆H≥1 ), a triangle average accumulator is required. As the accumulator types used
in the MFSS toolbox are introduced below, the manner in which they fit data sampling and
transformation types (summarized in Table 1) are discussed.7
Sum accumulators
Sum accumulators are those where the lower frequency is the sum of a time series of higher
frequency observations. For a low-frequency time series that is the sum of the previous n
7
Given that the triangle average collapses to the simple average when H = 1, we adopt the convention
(with a slight abuse of notation) of setting H = 1 for any series that is a level in the MFSS toolbox.
Journal of Statistical Software 7
observations, we would have
ytl = yth + yt−1
h
+ · · · + yt−n+1
h
= Z l (αt + αt−1 + · · · + αt−n+1 ) + dlt + βtl xlt + εlt ,
where ytl is the observed low-frequency time series and yth is an unobserved, high-frequency
version of that series. We define a sum accumulator, St , recursively as follows.
(
αt t ∈ P1
St = (3)
St−1 + αt t∈
/ P1
To set this in the state space context, we define a calendar parameter st that is equal to 0 in
the first period of a low-frequency period and equal to 1 in all other periods. We can then
rewrite Equation 3 as St = st St−1 + αt . The augmented transition equation would then be
! " # ! " # " # " #
αt T 0 αt−1 c γ Rt
= t + t + t wt + η.
St Tt st St−1 ct γt Rt t
Simple average accumulators
Simple average accumulators are those where the lower frequency is the average of a time
series of high-frequency observations. For a low-frequency time series that is the average of
the previous n observations, we would have
n−1
1X
ytl = yh .
n i=0 t−i
To arrive at a recursive formulation for the average accumulator, we write out the averages
for a 3- and 4-period case below.
1 h 1 2 1 h 1 h
(yt + yt−1
h
+ yt−2
h
) = yth + yt−1 + yt−2
3 3 3 2 2
1 h 1 3 1 2 1 h 1 h
(y + yt−1 + yt−2 + yt−3 ) = yt +
h h h h
y h
+ y + y
4 t 4 4 3 t−1 3 2 t−2 2 t−3
From these examples, we can see that the (m − 1)-period average, appropriately weighted,
can be combined with the next observation to produce an m-period average. For the general
case, we define an averaging calendar mt equal to the number of high-frequency periods
within the low-frequency period.8 For example, with a monthly base frequency, a quarterly
average accumulator’s calender would be the repetition of (1, 2, 3). Differing numbers of
high-frequency periods within a low-frequency period are allowed by changing how high the
mt counts before it resets to 1, as is required for instance with monthly observations of a
weekly process. We then recursively define the average accumulator Mt as
1 mt − 1
Mt = αt + Mt−1 .
mt mt
8
mt = t − pt + 1 where pt = max P1pt ≤t .
8 MATLAB Toolbox for Mixed Frequency State Space Models
To augment the state space with this accumulator, we can write
! " # ! " # " # " #
αt Tt 0 αt−1 ct γt Rt
= + + wt + η.
Mt Tt /mt (mt − 1)/mt Mt−1 ct /mt γt /mt Rt /mt t
Triangle average accumulators
A long-differenced transformation with an average sampling process requires an alteration
to the simple average accumulator, and thus an additional accumulator type. To illustrate,
consider a model where we observe quarterly GDP with a monthly base frequency. Denote
monthly levels of GDP as Xt and its base frequency-differences as xt = Xt − Xt−1 . We want
to relate the observed year-over-year change in quarterly GDP (i.e., the Q4/Q4 change in
GDP) to the latent monthly first-differences.9 At the end of the second year of the sample,
this is the 12-month change (a long-difference) in the 3-month average of Xt :
(3,12) 1 1
x24 = (X22 + X23 + X24 ) − (X10 + X11 + X12 )
3 3
22
1
" #
= 3X10 + 3 xt + 2x23 + x24 − (3X10 + 2x11 + x12 )
X
3 t=11
22
1 2 2 1
= x11 + x12 + xt + x23 + x24
X
3 3 t=13
3 3
We can see that the base frequency-differences are averaged using a set of weights that forms
a trapezoid ending at time t. The name “triangle” average comes from the most common
case in which the horizon of the change in the averaged quantities matches the difference in
frequencies, in which case the inner portion of the trapezoid collapses and the weights form
a triangle. This accumulator, with an averaging period S and a differencing period H (and
D = min{S, H}), is defined on the differenced time series so that the following relationship
to its level holds,
(S,H) 1 S−1 1 S−1
=
X X
xt Xt−i − Xt−i−H
S i=0 S i=0
t−H−S+D+1
1
= xt + 2xt−1 + · · · + Dxi + · · · + 2xt−H−S+3 + xt−H−S+2 .
X
S i=t−D+1
To put this into recursive form, for a given horizon the m-period average can be written in
terms of the (m − 1)-period average, similar to the earlier result, as
(m,H) 1 (1,H) m − 1 (m−1,H)
xt = x + xt−1 .
m t m
9
If this model is run on log-levels of U.S. GDP, as in Mariano and Murasawa (2010), this treatment involves
an approximation given that the standard accounting identity used for quarterly GDP is the arithmetic average
and not the geometric average of the time series. The approximation is necessary, however, if the linearity of
the state space model is to be preserved. Mitchell, Smith, Weale, Wright, and Salazar (2005) finds this to be
a good first-order approximation in the case of GDP.
Journal of Statistical Software 9
(1,H)
Since the 1-period average of an H-difference series is xt = Xt − Xt−H = H−1i=0 xt−i , we
P
have a slightly modified form of the simple average accumulator where we have to take the
sum of the past H base-frequency observations into account,
1 mt − 1
Mt = (αt + αt−1 + · · · + αt−H+1 ) + Mt−1 .
mt mt
This implies that we will need to make sure that there are at least H − 1 lags of the state we
are accumulating in αt . To augment the state with this accumulator, we can write
0 0
αt Tt αt−1 ct γt Rt
{α } H =
t−i i=1 L H−1,n 0 t−i i=2 t 0 ηt
{α } H+1
+ 0 + 0
w +
Tt +Im I(H−1)m mt −1 ct γt Rt
Mt m m mt
Mt−1 mt mt mt
t t
where LH−1,n is a lag matrix ensuring that we have a total of H − 1 lags of the n states.10
Accumulated observations
Restricting Z l to be time-invariant above allowed us to focus on creating low-frequency ver-
sions of the latent states that can be used in the observation equation of low-frequency time
series. This can provide substantial performance gains when estimating these models, since
the number of observations is commonly far larger than the number of states (e.g., in dynamic
factor models) and keeping the size of the state small has computational advantages. It is
also possible to define accumulators with a time-varying observation matrix, Ztl , but only at
the expense of expanding the state by as many observables requiring accumulation. Following
Harvey (1989), we instead bring the observation into the state equation as follows:
" # " #" # " # " #" # " #
ytb Zb 0 0 αt db βb 0 xbt εb
= t + t + t + t
ytl 0 0 1 l
yt 0 0 0 0 0
0
" # " #" # " # " # " #" #
αt Tt αt−1 ct γt Rt 0 ηt
= + + γ wt + .
l
yt Ztl ATt Zt Aζt
l l
yt−1 Zt At + dt + βt xt
l c l l l l
Z t At Z t At 1
l R εlt
For models with a time-varying Z l , this redefinition allows for the accumulators to still be
calculated in the state, since it forces the redefined Z l to be time-invariant.
4. Mixed frequency state space applications
The primary contribution of the MATLAB toolbox MFSS is to allow for easy manipulation of
the accumulator specifications described above across a variety of mixed frequency datasets
and applications. The toolbox also readily handles the computation of the state decompo-
sitions, in addition to being able to match the estimation capabilities and computational
performance of other packages. While each application to follow includes code snippets that
illustrate the functionality and appropriate syntax of key functions found in the toolbox,
Table 2 provides an overall summary of the primary functions with a short description.
To illustrate how to use the MFSS toolbox and appropriately handle a variety of different
modeling aspects of a mixed frequency application, we use it in five different empirical appli-
cations of topical interest in the rest of this section. The five applications are: (i) a dynamic
10
The triangle average collapses to the simple average when H = 1.
10 MATLAB Toolbox for Mixed Frequency State Space Models
Function name Description
Accumulator Temporal aggregation accumulator structure
.GenerateRegular() Generate accumulator for a regularly spaced intervals
.GenerateFromDates() Generate accumulator from a specified date vector
.augmentStateSpaceEstimation() Augment estimation object with required accumulators
.augmentStateSpace() Augment system with known parameters
StateSpaceEstimation Build state space structure with unknown parameters
.estimate() Estimate the parameters of the state space system
ThetaMap Controls StateSpaceEstimation parameter restrictions
.addRestrictions() Add restrictions to state space parameters
.addStructuralRestriction() Add restrictions on structural parameters
StateSpace Build state space structure with known parameters
.filter() Construct filtered estimates of the state
.smooth() Construct smoothed estimates of the state
.decompose_smoothed() Construct contributions decomposition of state
.smoothSample() Draw from the state using the simulation smoother
Table 2: Main functions for setting up state space. This table summarizes most of the main
functions/structures used in the MFSS toolbox for setting up the state space accommodating
any mixed frequency aspects of the model. For each function, a short description is provided.
factor model of the business cycle, (ii) a small monetary policy VAR, (iii) a stochastic volatil-
ity model of world trade, (iv) a trend-cycle decomposition of GDP, and (v) an estimate of
the natural rate of interest.
Table 3 summarizes the key data used for each of these applications. Almost all the data
used in the examples and included in the MFSS toolbox are retrieved from FRED, Federal
Reserve Bank of St. Louis database. There are two noticeable exceptions. First, the volume
of world trade measure, used in the third example, was taken from the CPB Netherlands
Bureau for Economic Policy Analysis. Second, for the single frequency version of the natural
rate of interest model (i.e., example 5), we use a historical release of the publicly available
data maintained by the Federal Reserve Bank of New York associated with the Laubach and
Williams (2003)’s R-star model.11
While each of these applications vary in their object of interest and empirical strategy, what
they share is a tradition of being estimated in a single frequency set-up. And, while the recent
literature has developed mixed frequency estimates of (i) and (ii), the results to follow offer
the first mixed frequency estimates of the volatility of trade, the trend-cycle decomposition
of GDP, and the natural rate of interest.
4.1. Dynamic factor model (DFM) of the business cycle
Our first application introduces the MFSS toolbox and its core functionality. For this in-
troduction, we estimate perhaps one of the very first successful mixed frequency economic
applications; a dynamic factor model of the business cycle. Specifically, we estimate a model
similar to the mixed frequency dynamic factor model of Mariano and Murasawa (2003), which
11
For the mixed frequency version of the model in example 5, we pull all the data from the FRED database.
Journal of Statistical Software 11
Variable Application Frequency Source
GDP 1, 2, 4, 5⋆ Quarterly @GDPC1
Payroll Employment 1, 2 Monthly @PAYEMS
Pers. Income Less Transfers 1 Monthly @W875RX1
Industrial Production 1 Monthly @INDPRO
Real Mnfg & Trade Sales 1 Monthly @CMRMTSPL
CPI 2 Monthly @CPIAUCSL
Commodity Price Index 2 Monthly @PPIACO
Federal Funds Rate 2, 5⋆ Monthly @DFF
Trade Weighted Dollar 3 Daily @DTWEXB
World Trade 3 Monthly ⋆⋆
PCE Price Index 5⋆ Monthly @PCE
Import Prices 5⋆ Monthly @IR
Oil Prices 5⋆ Monthly @DCOILWTICO
Table 3: Data sources. This table summarizes the data used across our set of applications, as
well as their source. Series labeled with @ represent the mnemonic associated with accessing
this series on the FRED, Federal Reserve Bank of St. Louis database. ⋆ The mixed frequency
version in example 5 uses data from FRED, using the mnemonics reported, while the single
frequency version uses data we gathered directly from the Laubach and Williams web-page
for the R-star model, which now can be found at: https://www.newyorkfed.org/research/
policy/rstar. ⋆⋆ The volume of world trade series is made available by the CPB Netherlands
Bureau for Economic Policy Analysis, and can be found here: https://www.cpb.nl/en/
worldtrademonitor.
involves using quarterly GDP together with a selection of monthly indicators of economic ac-
tivity to extract a common (monthly) factor that is linked to the business cycle.
The following code outlines the few steps involved in estimating the dynamic factor model
using MFSS. The steps towards estimating the model involves: (i) pre-allocating the system
matrices of the state space model — with nan’s for parameters to be estimated, (ii) build-
ing the state space estimation structure, (iii) setting up the lower and upper bounds for the
relevant parameter values, (iii) setting the properties of the required accumulator and aug-
menting the state space structure with this information, and (iv) estimating the model. Each
of these steps is facilitated by the MFSS functions including: StateSpaceEstimation(),
addRestrictions(), GenerateRegular(), augmentStateSpaceEstimation().
Z = [1; nan(4,1)];
d = [nan; zeros(4,1)];
H = diag(nan(5,1));
T = nan;
Q = nan;
ssE = StateSpaceEstimation(Z, H, T, Q, 'd', d);
LB = ssE.ThetaMapping.LowerBound;
UB = ssE.ThetaMapping.UpperBound;
LB.Z(2:5,1) = 0;
LB.T(1,1) = -1;
12 MATLAB Toolbox for Mixed Frequency State Space Models
Figure 1: DFM estimation results. This figure displays the time series used in the estimation
of the dynamic factor model together with the estimate of the common factor. In the top
panel: (i) payroll employment (blue), (ii) real personal income less transfers (orange), (iii)
industrial production (yellow), (iv) real manufacturing and trade sales (purple), as well as
the estimate of the factor (black) are reported in standardized units. In the bottom panel,
GDP (blue) and the factor are reported in annualized growth units.
UB.T(1,1) = 1;
ssE.ThetaMapping = ssE.ThetaMapping.addRestrictions(LB, UB);
accum = Accumulator.GenerateRegular(y, {'avg', ' ', ' ', ' ', ' '},...
[3 1 1 1 1]);
ssEA = accum.augmentStateSpaceEstimation(ssE);
ssML = ssEA.estimate(y);
alphaHat = ssML.smooth(y);
The observed time series in this model is a panel of five time series at a monthly base frequency.
The first series in y is quarterly with observations in the third month of the quarter. Each
series has been log-differenced and standardized. We require a triangle average accumulator
with a horizon of three months given the monthly base frequency of this model. After creating
the accumulator object, we augment the state space to align with the mixed frequency nature
the data, estimate the parameters via maximum likelihood, and obtain the smoothed state
given the estimated parameters.
When estimating reduced form models like the one above, we often need to restrict the values
the parameters can take to ensure identification. The simplest parameter restrictions occur
Journal of Statistical Software 13
when elements of the state space matrices need to be restricted to obey explicit bounds. This
restriction above requires that the AR(1) controlling the dynamics of the latent factor ensures
that it remains stationary throughout estimation by restricting the appropriate element of T
to be between −1 and 1.
Results for the estimation of the DFM are presented in Figure 1. A couple of points are
worth highlighting. First, the estimated factor is a clear coincident indicator of the business
cycle. Positive realizations of the factor generally coincide with business cycle expansions,
while negative realizations are typically associated with contractionary periods (gray shaded
periods), with very few false signals in either case. Given the idiosyncratic variation each
series exhibits beyond the common factor, it is clear that being able to jointly include both
key quarterly (e.g., GDP) and monthly variables (e.g., payroll employment) into the dynamic
factor model facilitates the extraction of the business cycle “signal” from the inherent noise
present in many macroeconomic time series.
4.2. Small monetary policy mixed frequency VAR (MF-VAR)
The use of VARs to identify monetary policy shocks has a long tradition. Ideally, these VARs
inform the extent of the impact of unanticipated monetary policy interventions on the ultimate
path of economic activity and inflationary pressures. In its most concise form, a monetary
policy VAR is comprised of an economic activity variable (e.g., payroll employment or GDP),
an aggregate price variable (e.g., the Consumer Price Index, CPI), and the federal funds rate.
Once estimated, one can examine impulse response functions, variance decompositions and/or
forecasts of economic activity and price levels to historical policy surprises.
Even in this admittedly small system, a critical trade-off exists in the modeling decisions
involving frequency and variable selection. While monthly variations in monetary policy and
prices are observed, some of the most critical aggregate measures of economic activity are
only observed at the quarterly frequency, most notably GDP. Alternatively, other perhaps
less desirable measures of aggregate activity are available at the higher monthly frequency;
such as monthly payroll employment. So, without the ability of estimating the VAR in a
mixed frequency setting, the econometrician is left with a less than ideal trade-off that is
likely to distort the ultimate set of inferences trying to be recovered. A perfect example of
the ambiguous nature of this trade-off comes from the different choices made by Bańbura,
Giannone, and Reichlin (2010) and Koop (2013) in their examination of the role of Bayesian
shrinkage techniques on the performance of a small monetary policy VAR.
By leveraging a mixed frequency setup and the relative ease in which a model of this size can
be estimated using the MFSS toolbox, we alleviate the potential drawbacks of their alternative
modeling choices. Specifically, we estimate a four variable monthly VAR comprised of GDP,
CPI, commodity prices, and the federal funds rate. Here, the mixed frequency aspect of the
application stems from the fact that we observe the price indexes and the federal funds rate at
a monthly frequency, and GDP at a quarterly frequency. To illustrate the different inferences
such an approach can yield relative to the most natural single frequency alternatives, we
examine the impulse response functions of GDP to a monetary policy shock for our mixed
frequency VAR(6) model in addition to a quarterly version of the model (with two lags),
as well as a monthly version that replaces GDP with monthly payroll employment as the
measure of economic activity (Bańbura et al. 2010).
14 MATLAB Toolbox for Mixed Frequency State Space Models
Figure 2: Monetary Policy VAR estimation results. This figure plots the estimated impulse
response functions to a 100 basis point increase to the federal funds rate for GDP and payroll
employment for the three different small monetary policy VARs. The impulse for payroll
employment (yellow) and GDP from the mixed frequency VAR (blue) are plotted over 48
months after the shock, while the impulse for GDP from the quarterly VAR (red crosses) are
shown for the parallel 12 quarters.
Figure 2 reports the resulting impulse response functions for economic activity from a 100
basis point increase to the federal funds rate using data from 1965–2017.12 Each of the three
models result in similar interpretations regarding an increase in the federal funds rate: a
contraction in economic activity. The differences across the single frequency quarterly and
monthly models highlight the potential ambiguity facing policymakers faced with having to
make assessments in real time given the publication delays associated with many headline, but
lower frequency, releases like GDP. Here, the quarterly VAR that used GDP as its measure of
economic activity predicts a much steeper descent as a result of a federal funds rate increase
than the monthly VAR which used payroll employment.
Interestingly, the mixed frequency model seems to strike a balance between the two models
in terms of its interpretation of the depth of the contraction in GDP following a federal
funds rate increase. Taking a closer look at the mixed frequency model’s impulse response,
it is clear that the more immediate effects (within two quarters after the shock) resemble
12
The code supplement reports the estimated VAR coefficients for each model. Alone, these coefficient
estimates can often be difficult to interpret, and is why we choose to focus on the impulse response functions.
Journal of Statistical Software 15
the short run impacts of the quarterly model (which also used GDP as its measure), while
the longer run dynamics (beyond two years) more closely resemble the monthly frequency
model. In practice, the mixed frequency model would also facilitate a richer set of possible
identification strategies to the extent that higher frequency movements in monetary policy
are able to identify critical responses of interest for policymakers.
4.3. Stochastic volatility
A common feature of financial and high frequency time series is the clustering of high and low
volatility periods. In many cases, this feature of economic and financial times series represents
an object of interest. For example, recently a lot of interest amongst both policy makers and
financial markets centers on the volatility of international trade. Conveniently, a common
modelling approach to accommodate such time series behavior, the stochastic volatility model,
translates well to a state space approach. Often times the relevant indicators that might
provide information on the evolution of volatility occur across a range of frequencies. In
the case of international trade, while exchange rate information can be readily found at the
daily frequency, the most comprehensive measures of trade typically are reported at lower
frequencies (e.g., monthly).
In the case of trying to recover the volatility in international trade, if a single frequency
approach is taken a critical tradeoff exists. One must decide between averaging over the
possible identifying variation in the higher frequency variation embedded in exchange rates,
versus using a less comprehensive or representative series that measures trade. By adopting
a mixed frequency setup, we overcome this tradeoff. Specifically, we estimate a bi-variate
common stochastic volatility model comprised of the trade weighted dollar exchange rate and
the total volume of world trade. Here, the mixed frequency aspect of the application stems
from the fact that we observe the exchange rate at the daily frequency, and volume of world
trade at the monthly frequency.
Specifically, we assume that a common stochastic volatility component (θt ) affects both the
volume of trade as well as the real exchange rate (in log first differences), summarized by the
following relationship
1
yit = µi + σi exp θt εit
2
where εit is distributed as a standard normal distribution, and µi and σi represent series
specific mean and variance components. At least in terms of the time varying log volatility
component, θt , this is a nonlinear model, and could require a nonlinear state space system.13
We adopt the approach outlined in Durbin and Koopman (2012) in estimating the stochas-
tic volatility component, by way of an approximation via data transformation. While this
approach leads to an approximation, it ensures that the ultimate state space model that is
estimated remains linear.
After assuming the common volatility component (θt ) admits an AR(1) process, a transfor-
mation of the data allows us to use the following state space model:
log yt2 = κ + θt + ξt
θt = ϕθt−1 + ηt ηt ∼ N (0, ση2 ),
13
For a discussion on nonlinear state space models, see Durbin and Koopman (2012). For a MATLAB package
that has some capabilities on nonlinear filtering, see Villegas and Pedregal (2018).
16 MATLAB Toolbox for Mixed Frequency State Space Models
where κi = log σi2 + E(log ε2it ), and ξit = log ε2it − E(log ε2it ), and −1 < ϕ < 1. Now, while the ξt
will not be normally distributed, it will be the case that we can estimate the parameters using
the same linear Kalman filter techniques capable in our toolbox through a quasi-maximum
likelihood (QML) approach (e.g., Harvey, Ruiz, and Shephard 1994).
We apply this approach to our bi-variate system of the daily trade weighted exchange rate, and
the monthly trade volume of world trade. Importantly, this mixed frequency state space model
requires a temporal aggregation scheme that must take into account the irregular number of
days that occur within each month of the year. To facilitate this temporal aggregation of
θt for the monthly volume of world trade, we assume that the series we observe log ỹit 2 =
M (t) (log yit + . . . + log yit−M (t)+1 ), which states the monthly realization of log volatility is the
1 2 2
geometric standard deviation of the M (t) underlying base frequency realizations, the number
of which–of course–will vary for each month.14
The following code outlines the few steps involved in estimating the stochastic volatility model
using MFSS. In addition to the already discussed steps of pre-allocating the system matri-
ces of the state space model and building the (augmented) state space estimation structure
with the relevant temporal aggregation schemes, the code here also illustrates the use of the
GenerateFromDates() function for building the accumulator structure using a pre-specified
date vector, in addition to taking a draw from the state (post-estimation) using the simulation
smoother using the .smooth() function.
accum = Accumulator.GenerateFromDates(dates, {'Monthly', ''}, [1 0]));
Z = ones(2,1);
d = nan(2,1);
H = diag(nan(2,1));
T = nan;
Q = nan;
mdl = StateSpaceEstimation(Z, H, T, Q, 'd', d);
mdlA = accum.augmentStateSpaceEstimation(mdl);
mdlA.P0 = diag(Inf(2,1));
theta0 = [-0.9347; -4.1606; 1.5648; 1.8696; 0.0203; -0.0004];
mdlMLE = mdlA.estimate(Y, theta0);
alpha = mdlMLE.smooth(Y);
alphaSim = mdlMLE.smoothSample(Y, [], [], [], 100);
sigmaHat = exp(0.5 .* alpha(:,1));
sigmaHatDraws = squeeze(exp(0.5 .* alphaSim(:,1,:)));
sigmaHatBands = prctile(sigmaHatDraws, [5 95], 2);
Figure 3 plots the (smoothed) estimate of the common log-volatility of international trade
at the daily frequency, or more specifically exp 12 θt . In addition to the smoothed estimate
of the log-volatility component, we also report the implied 95 percent confidence bands from
1000 draws of the state using the simulation smoother. Given the daily base frequency our
mixed frequency approach affords us, our measure of volatility balances the high frequency
information in the daily time series of exchange rate with the lower frequency information of
14
Like the temporal aggregation schemes involving the geometric average of the level of variables, this
treatment involves an approximation given that the standard accounting identity is the arithmetic average.
The approximation is necessary, however, if the linearity of the state space model is to be preserved, and is
typically found to be an innocuous approximation (e.g., Corsi 2009, see footnote 4).
Journal of Statistical Software 17
Figure 3: Stochastic volatility in trade results. This figure plots the smoothed estimates of
the common log-volatility component of the trade-weighted dollar exchange and the world
volume of trade. In addition to the smoothed estimate of the log-volatility component, we
also report the implied 95 percent confidence bands from 100 draws from the state using the
simulation smoother.
world trade. Inspection of Figure 3 indicates an increase in the volatility of trade occurred
both around the Great Recession and the Iraq War. Furthermore, an escalation in trade
volatility also appears in conjunction with the escalating trade tensions between the U.S. and
China since 2016.
As the availability of higher frequency information continues to become increasingly available
to researchers, the value of being able to handle irregular temporal aggregation schemes are
likely to increase and become more nuanced across applications. For instance, in the case of
daily and weekly time series, the temporal aggregation properties to lower frequencies will have
an irregular scheme due to calendar effects, and these effects permeate to the characterization
of the uncertainty of the state. Thus, the built-in capabilities of the MFSS package to handle
a broad set of mixed frequency schemes in addition to characterizing the uncertainty of
(latent) elements of the state–highlighted in this application–will be increasingly valuable to
researchers.
4.4. Trend-cycle decomposition of GDP
The next application takes a closer look at separating the trend and business cycle components
of GDP. Due to the inherent non-stationary behavior of time series like GDP, establishing
empirical regularities as they relate to business cycle phenomenon versus the longer run
movements in GDP often require such a decomposition. The trend-cycle decomposition, in
particular, serves both as an object of interest in its own right (e.g., Antolin-Diaz, Drechsel,
18 MATLAB Toolbox for Mixed Frequency State Space Models
and Petrella 2017), and also as a critical input to any multivariate analysis that attempts to
explain the cyclical co-movements in the economy (e.g., Hodrick and Prescott 1997).
Given these demands, it seems costly to require that such decomposition only be available at
the frequency of observation — especially for a series as important as GDP. To estimate such
a decomposition at a higher frequency than the frequency of observation naturally lends itself
to a mixed frequency set up. In particular, we attempt to decompose GDP into its trend
and cycle components as in Harvey and Jaeger (1993), but at the monthly frequency. This
application has the added benefit of illustrating the MFSS toolbox’s flexibility in being able
to handle more complicated restrictions in the estimation process.
Evaluating the likelihood of a vector of time series given a proposed parameter vector θU is, in
this instance, a 3-step process — an unrestricted “structural” parameter vector θU is domain
restricted according to element-wise restriction functions R to θ, transformed to reduced-form
parameters ψ, then element-wise transformed into the state space parameters
θi = Ri (θiU ) → ψi = Ψi (θ) → Xj,i = TX,j (ψi )
for X ∈ {Zt , dt , βt , Ht , Tt , ct , γt , Rt , Qt }. In cases where multiple parameters are related, we
define the state space matrices as functions of the underlying structural parameters.
The stochastic trend-cycle decomposition of Harvey and Jaeger (1993) is exactly a setting
where such functionality is required. For instance, a version of the stochastic trend-cycle
decomposition model is as follows:
yt = µt + ψt
µt = µt−1 + ϕt−1
ϕt = ϕt−1 + ξt
" # " #" # " #
ψt cos λ sin λ ψt−1 κ
=ρ + ∗t ,
∗
ψt − sin λ cos λ ∗
ψt−1 κt
where ξt is normally distributed and κt and κ∗t are independently normally distributed with
a common variance. The structural parameters ρ and λ are restricted so that the cyclical
component remains stationary with an expected period between 1.5 and 12 years. Casting
the system into state space form, we have
h i
yt = 1 0 1 0 α t
1 1 0 0 0
0 1 0 0 ξ
αt = α + t.
0 0 ρ cos λ ρ sin λ t−1 κt
0 0 −ρ sin λ ρ cos λ κ∗t
The 4 structural parameters in this model are stacked to form θ = [λ, ρ, log σξ2 , log σκ2 ]⊤ and
transformed to their reduced form ψ = [ρ cos(λ), ρ sin(λ), −ρ sin(λ), log σξ2 , log σκ2 ]⊤ . The user
needs only to write out the model in terms of the underlying structural parameters to be
estimated, as the following code demonstrates. The steps unique to this application worth
highlighting are the use of MATLAB’s symbolic expressions capabilities (e.g., through the syms
command) to define the state space system matrices in terms of their structural parameters
(e.g., in system matrix T), and then constructing the parameter restrictions on those structural
Journal of Statistical Software 19
Figure 4: Trend-cycle decomposition estimation results. This figure plots the components
of the three trend cycle decompositions estimated for GDP. Across all of the panels, the
mixed frequency Harvey-Jaeger decomposition is reported in blue, while the quarterly (single-
frequency) Harvey-Jaeger decomposition is reported in orange. The Hodrick-Precott filter
trend-cycle decomposition is reported in yellow. Recessions as defined by the NBER are
shaded gray for the panels involving the cycle component (right panels).
objects themselves using the .addStructuralRestriction() functionality. This application
also benefits from well chosen set of initial conditions, which the example code illustrates how
to implement.
syms lambda rho sigmaKappa sigmaXi
Z = [1, 0, 1, 0];
H = 0;
T = blkdiag([1 1; 0 1], ...
rho .* [cos(lambda), sin(lambda); -sin(lambda) cos(lambda)]);
R = [zeros(1, 3); eye(3)];
Q = diag([sigmaXi; sigmaKappa; sigmaKappa]);
ssE = StateSpaceEstimation(Z, H, T, Q, 'R', R);
accum = Accumulator.GenerateRegular(y, {'avg'}, 1);
ssEA = accum.augmentStateSpaceEstimation(ssE);
ssEA.ThetaMapping = ssEA.ThetaMapping.addStructuralRestriction(rho, 0, 1);
ssEA.ThetaMapping = ...
ssEA.ThetaMapping.addStructuralRestriction(lambda, pi/72, pi/9);
ssML = ssEA.estimate(y, [0.0943; 0.9610; log(0.00003379); ...
log(0.0000003789)]);
alpha = ssML.smooth(y);
20 MATLAB Toolbox for Mixed Frequency State Space Models
Figure 4 displays the mixed frequency version of the Harvey and Jaeger (1993) trend and
cycle decomposition of GDP (top panels) as well as their growth rates (bottom panels), us-
ing GDP data from 1947q1–2018q2. For comparison purposes, we also report the same time
series components for the Harvey and Jaeger (1993) trend-cycle decomposition estimated at
the quarterly frequency, as well as the Hodrick-Precott (HP) filtered estimate of the decom-
position (Hodrick and Prescott 1997).15 An appealing fact involving these comparisons is
the strong similarity across each of the different models estimates. Their similarities suggest
that the mixed frequency approach’s ability to provide an estimate of the underlying monthly
decomposition (something the other two models cannot provide) of the trend and cycle might
lead practitioner’s to favor the mixed frequency approach.
Though the three models do largely interpret the trend and cycle similarly, one interesting
departure of the mixed frequency version is the smaller variation in the trend for GDP (see
bottom left panel) relative to both of the quarterly frequency models. This departure while
influencing the relative amplitude of the cycle component of GDP (see right panels), largely
does not change the strong coincidence of the cycle with NBER defined business cycles (shaded
in gray). Additionally, the mixed frequency model’s parameter estimates (e.g., the variance of
the shock to trend, σξ2 ) are more natural candidates for semi-calibrated models of the economy
modeled at the monthly frequency.
4.5. Estimating the natural rate of interest
Our final application showcases the MFSS toolbox’s ability to calculate the full set of Kalman
filter and smoother weights to create useful measures of how each series contributes to the
movements in the underlying latent state. To showcase this functionality, we estimate the
natural rate of interest (e.g., Laubach and Williams 2003) — an inherently latent object of
substantial economic importance. What makes this effort particularly suitable for a mixed
frequency setup is that many of the variables incorporated in the estimation of the natural
rate of interest (e.g., GDP, inflation, and the federal funds rate) are observed at very different
frequencies, ranging from the daily frequency (e.g., the federal funds rate) to the quarterly
frequency (e.g., GDP). To date, what this has meant is that publicly available estimates of
the natural rate of interest have been restricted to a quarterly frequency. By fully leveraging
the capabilities of the MFSS toolbox, we provide the first estimate, to our knowledge, of the
natural rate of interest at the monthly frequency.
In order to estimate this model at a monthly frequency, the state space specification of
Laubach and Williams (2003) must be adapted so that monthly GDP is a latent variable.16
We initialize the structural parameters at the estimated values found by Laubach and Williams
(2003) and attempt to mimic their handling of the initial state. Further details can be found
in the accompanying User Guide. Finally, to compute the full set of smoother weights, the
code base utilizes the function decompose_smoothed() taking the exogenous variables as
additional inputs (i.e., X, W ) beyond the time series of the model (i.e., Y ).
15
For the Hodrick-Prescott (HP) filter we use a smoothing parameter of 1600 which is standard for quarterly
frequency series.
16
The mixed-frequency version of the model follows the tradition of Laubach and Williams (2003), and uses
the monthly versions of changes in the PCE price index, the monthly real interest rate constructed using
inflation expectations from an AR model of observed inflation, and measures of import prices and oil prices,
together with quarterly GDP.
Journal of Statistical Software 21
syms a1 a2 a3 b1 b2 b3 b4 b5 c sigma2Ystar sigma2IS sigma2PC
Z = [0 0 1 0 zeros(1,5);
-b3 0 b3 0 zeros(1,5)];
beta = [zeros(1,5); b1 b2 (1-b1-b2) b4 b5];
H = diag([0 sigma2PC]);
T = [1 0 0 0 1 0 0 0 0; 1 zeros(1,8);
1-a1 -a2 a1 a2 1-(12*c*a3/2) -12*c*a3/2 -a3/2 -a3/2 0;
0 0 1 0 zeros(1,5); zeros(1,4) 1 zeros(1,4); zeros(1,4) 1 zeros(1,4);
zeros(1,6) 1 zeros(1,2); zeros(1,6) 1 zeros(1,2);
0 0 0 0 12*c 0 1 0 0];
gamma = [zeros(2); a3/2 a3/2; zeros(6,2)];
R = [1 0 0 0; zeros(1,4);
1 0 1 0 ; zeros(1,4);
0 lambda_g/3 0 0; zeros(1,4);
0 0 0 lambda_z/(3*a3); zeros(1,4);
zeros(1,4)];
Q = diag([sigma2Ystar sigma2Ystar sigma2IS sigma2IS]);
ssE = StateSpaceEstimation(Z, H, T, Q, 'beta', beta, 'gamma', gamma, 'R', R);
ssE.ThetaMapping = ssE.ThetaMapping.addStructuralRestriction(a3, [], ...
-0.0025);
ssE.ThetaMapping = ssE.ThetaMapping.addStructuralRestriction(b3, 0.025, []);
accum = Accumulator.GenerateRegular(YM, {'avg', ' '}, [1 0]);
ssEA = accum.augmentStateSpaceEstimation(ssE);
interpGDP = interp1(3:3:size(gdpM,1), gdpM(3:3:end), 1:size(gdpM))';
lagGDP = 100*interpGDP(23:24)';
gInit = 0.9394933/3;
ystarInit = [809.7823659 809.7823659-gInit];
rstarInit = 12*gInit;
ssE.a0 = [ystarInit, lagGDP, repmat(gInit, [1 2]), zeros(1, 2), rstarInit]';
ssE.P0 = 0.2*eye(9);
ssE.P0(3,3) = 0; ssE.P0(4,4) = 0;
ssEA.a0 = [ssE.a0; 100*mean(lagGDP)];
ssEA.P0 = blkdiag(ssE.P0, 0.2);
[ssMOpt, ~ , thetaMOpt] = ssEA.estimate(Y, thetaR, X, W);
alphaM = ssMOpt.smooth(Y, X, W);
[yContribM, paramContribM, inflContribM, ratesContribM] = ...
ssMOpt.decompose_smoothed(Y, X, W);
Figure 5 reports state decompositions using the filter/smoother weights for each of the time
series involved in the estimation of the natural rate of interest over the full period of the esti-
mation sample (1959q1–2018q1).17 For comparison purposes, we report along with the mixed
frequency (monthly) estimate of the natural rate of interest (top panel), an estimate from the
traditional quarterly frequency model (bottom panel). In each panel, the contributions from
each series involved in the estimation of the natural rate of interest (in black) is reported
(e.g., GDP, inflation, real rate of interest, oil import prices, and import prices).
17
For an explicit formulation of the model equations involved with estimating the natural rate of interest,
see the accompanying code in addition to Laubach and Williams (2003).
22 MATLAB Toolbox for Mixed Frequency State Space Models
Figure 5: Natural rate of interest estimation results. This figure reports contributions of each
time series to an estimate of the natural rate of interest (r⋆ ). The smoothed estimate of the
natural rate of interest is reported in black, while contributions to it driven by (i) GDP are
reported in blue, (ii) inflation are reported in orange, (iii) the real rate of interest are reported
in yellow, (iv) oil import prices are reported in purple, and (v) import prices are reported in
green.
Across both panels, it is clear that the primary drivers of the fluctuations in the natural rate
of interest have been the real rate of interest and GDP. While the real rate of interest has
supported higher levels of the natural rate of interest in some periods, 1980s–2000s, GDP has
largely been bringing the natural rate of interest lower since the 1960s. While both versions
of the model are in agreement with the countervailing roles of GDP’s and the real rate of
interest’s influence on the natural rate of interest, the mixed frequency version interprets a
stronger divergence between the two series in the 1980–2000 period and less divergence in the
period since 2010 compared to the quarterly version of the model.
The capabilities of the MFSS package to incorporate the filter/smoother weight calculations
for many types of variables within a mixed frequency state space formulation are underscored
by this application. As the development of higher frequency identification strategies become
more prevalent in economic and financial applications, these capabilities of the MFSS package
will only further elucidate identification challenges and implications of modeling assumptions.
5. Conclusion
Mixed frequency modeling through the use of state space methods is a valuable tool for a
wide array of economic and financial applications. Previously, using these methods required
Journal of Statistical Software 23
substantial effort by practitioners to implement even in the simplest of models. It is our hope
that this article together with the availability of our MFSS toolbox will increase the accessi-
bility of mixed frequency state space modeling and facilitate its adoption in a wider array of
empirical applications. By outlining its use across four very different types of applications,
we were able to highlight the unique set of capabilities of the MFSS toolbox. In each case,
adopting a mixed frequency framework was shown to facilitate better variable selection, more
parsimonious modeling choices, as well as a wider variety of inference possibilities.
The User’s Guide accompanying this article contains many more examples that may be of
interest to readers interested in using these methods. In the future, it will be important
to extend the capabilities of this package to incorporate the frontier in shrinkage estimation
strategies given the growth of densely parameterized models. Furthermore, as micro-economic
applications gain more access to higher frequency panel datasets, it will also be likely that
the relevant applications for this suite of methods will extend to these additional fields within
economics and other social sciences.
Acknowledgments
This research did not receive any specific grant from funding agencies in the public, commer-
cial, or not-for-profit sectors. We would like to thank the section editor, and a referee for
thoughtful suggestions, and Bill Kluender for excellent research assistance.
References
Antolin-Diaz J, Drechsel T, Petrella I (2017). “Tracking the Slowdown in Long-Run GDP
Growth.” The Review of Economics and Statistics, 99(2), 343–356. doi:10.1162/rest_a_
00646.
Aruoba SB, Diebold FX, Scotti C (2009). “Real-Time Measurement of Business Conditions.”
Journal of Business & Economic Statistics, 27(4), 417–427. doi:10.1198/jbes.2009.
07205.
Bańbura M, Giannone D, Reichlin L (2010). “Large Bayesian Vector Auto Regressions.”
Journal of Applied Econometrics, 25(1), 71–92. doi:10.1002/jae.1137.
Brave SA, Butters RA (2012). “Diagnosing the Financial System: Financial Conditions and
Financial Stress.” International Journal of Central Banking, 8(2), 191–293. ISSN 1815-
4654. doi:10.2139/ssrn.1668786.
Brave SA, Butters RA, Justiniano A (2016). “Forecasting Economic Activity with Mixed
Frequency Bayesian VARs.” Working Paper Series WP-2016-5, Federal Reserve Bank of
Chicago.
Carpenter B, Gelman A, Hoffman M, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo
J, Li P, Riddell A (2017). “Stan: A Probabilistic Programming Language.” Journal of
Statistical Software, 76(1), 1–32. ISSN 1548-7660. doi:10.18637/jss.v076.i01.
24 MATLAB Toolbox for Mixed Frequency State Space Models
Commandeur J, Koopman S, Ooms M (2011). “Statistical Software for State Space Methods.”
Journal of Statistical Software, 41(1), 1–18. ISSN 1548-7660. doi:10.18637/jss.v041.
i01.
Corsi F (2009). “A Simple Approximate Long-Memory Model of Realized Volatility.” Jour-
nal of Financial Econometrics, 7(2), 174–196. ISSN 1479-8409. doi:10.1093/jjfinec/
nbp001.
Doornik JA (2007). Object-Oriented Matrix Programming Using Ox. 3rd edition. Timberlake
Consultants Press and Oxford, London.
Durbin J, Koopman SJ (2012). Time Series Analysis by State Space Methods. 2nd edition.
Oxford University Press. doi:10.1093/acprof:oso/9780199641178.001.0001.
Foroni C, Marcellino M (2014). “Mixed-Frequency Structural Models: Identification, Esti-
mation, and Policy Analysis.” Journal of Applied Econometrics, 29(7), 1118–1144. doi:
10.1002/jae.2396.
Ghysels E, Santa-Clara P, Valkanov R (2004). “The MIDAS Touch: Mixed Data Sampling
Regression Models.” CIRANO Working Papers 2004s-20, CIRANO.
Harvey A, Chung CH (2000). “Estimating the Underlying Change in Unemployment in
the UK.” Journal of the Royal Statistical Society A, 163(3), 303–309. doi:10.1111/
1467-985x.00171.
Harvey A, Jaeger A (1993). “Detrending, Stylized Facts and the Business Cycle.” Journal of
Applied Econometrics, 8(3), 231–47. doi:10.1002/jae.3950080302.
Harvey A, Ruiz E, Shephard N (1994). “Multivariate Stochastic Variance Models.” The
Review of Economic Studies, 61(2), 247–264. doi:10.2307/2297980.
Harvey AC (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cam-
bridge University Press.
Harvey AC, Pierse RG (1984). “Estimating Missing Observations in Economic Time Se-
ries.” Journal of the American Statistical Association, 79(385), 125–131. doi:10.1080/
01621459.1984.10477074.
Helske J (2018). KFAS: Kalman Filter and Smoother for Exponential Family State Space
Models. URL https://cran.r-project.org/package=KFAS.
Hodrick RJ, Prescott EC (1997). “Postwar U.S. Business Cycles: An Empirical Investigation.”
Journal of Money, Credit and Banking, 29(1), 1–16. doi:10.2307/2953682.
Kim TB (2010). “Bayesian Mixed Frequency Estimation of DSGE Models.” Technical report,
Duke University.
Koop GM (2013). “Forecasting with Medium and Large Bayesian VARs.” Journal of Applied
Econometrics, 28(2), 177–203. doi:10.1002/jae.1270.
Koopman SJ, Doornik JA, Shephard N (2010a). SsfPack. URL http://www.ssfpack.com/.
Journal of Statistical Software 25
Koopman SJ, Durbin J (2000). “Fast Filtering and Smoothing for Multivariate State Space
Models.” Journal of Time Series Analysis, 21(3), 281–296. ISSN 1467-9892. doi:10.1111/
1467-9892.00186.
Koopman SJ, Durbin J (2003). “Filtering and Smoothing of State Vector for Diffuse State-
Space Models.” Journal of Time Series Analysis, 24(1), 85–98. doi:10.1111/1467-9892.
00294.
Koopman SJ, Harvey A (2003). “Computing Observation Weights for Signal Extraction and
Filtering.” Journal of Economic Dynamics and Control, 27(7), 1317–1333. doi:10.1016/
s0165-1889(02)00061-1.
Koopman SJ, Harvey A, Doornik JA, Shephard N (2010b). STAMP. URL http://www.
stamp-software.com/.
Kuzin V, Marcellino M, Schumacher C (2011). “MIDAS vs. Mixed-Frequency VAR: Now-
casting GDP in the Euro Area.” International Journal of Forecasting, 27(2), 529–542.
doi:10.1016/j.ijforecast.2010.02.006.
Laubach T, Williams JC (2003). “Measuring the Natural Rate of Interest.” The Review of
Economics and Statistics, 85(4), 1063–1070. doi:10.1162/003465303772815934.
Marcellino M (1999). “Some Consequences of Temporal Aggregation in Empirical Analysis.”
Journal of Business & Economic Statistics, 17(1), 129–136. doi:10.1080/07350015.1999.
10524802.
Mariano RS, Murasawa Y (2003). “A New Coincident Index of Business Cycles Based on
Monthly and Quarterly Series.” Journal of Applied Econometrics, 18(4), 427–443. doi:
10.1002/jae.695.
Mariano RS, Murasawa Y (2010). “A Coincident Index, Common Factors, and Monthly
Real GDP.” Oxford Bulletin of Economics and Statistics, 72(1), 27–46. doi:10.1111/j.
1468-0084.2009.00567.x.
Mitchell J, Smith RJ, Weale MR, Wright S, Salazar EL (2005). “An Indicator of Monthly
GDP and an Early Estimate of Quarterly GDP Growth.” The Economic Journal, 115(501),
F108–F129. ISSN 1468-0297. doi:10.1111/j.0013-0133.2005.00974.x.
Perktold J, Seabold S, Taylor J (2017). StatsModels. URL https://www.statsmodels.org/
stable/.
Proietti T (2006). “Temporal Disaggregation by State Space Methods: Dynamic Regres-
sion Methods Revisited.” The Econometrics Journal, 9(3), 357–372. doi:10.1111/j.
1368-423x.2006.00189.x.
Sargent TJ, Stachurski J (2018). “A First Look at the Kalman Filter.” URL https://
lectures.quantecon.org/py/kalman.html.
Schorfheide F, Song D (2015). “Real-Time Forecasting with a Mixed-Frequency VAR.” Journal
of Business & Economic Statistics, 33(3), 366–380. doi:10.1080/07350015.2014.954707.
26 MATLAB Toolbox for Mixed Frequency State Space Models
Shumway RH, Stoffer DS (1982). “An Approach to Time Series Smoothing and Forecasting
Using the EM Algorithm.” Journal of Time Series Analysis, 3(4), 253–264. doi:10.1111/
j.1467-9892.1982.tb00349.x.
The MathWorks Inc (2021). MATLAB – The Language of Technical Computing, Version
R2021a. Natick. URL https://www.mathworks.com/products/matlab/.
Villegas M, Pedregal D (2018). “SSpace: A Toolbox for State Space Modeling.” Journal of
Statistical Software, 87(5), 1–26. doi:10.18637/jss.v087.i05.
Journal of Statistical Software 27
A. Estimating the latent state
We are interested in the estimate of the latent state from a state space model with observed
data {yt , xt , wt } and time-varying parameters of the form
yt = Zt αt + dt + βt xt + εt εt ∼ N (0, Ht )
αt = Tt αt−1 + ct + γt wt + Rt ηt ηt ∼ N (0, Qt )
α1 ∼ N (a1 , P1 )
A.1. Multivariate filter
As shown in Durbin and Koopman (2012, §4.3), the filtered estimates of the state, at = E(αt |
Yt ) and Pt = VAR(αt | Yt ), are given by
at+1 = Tt+1 at + ct+1 + γt+1 wt+1 + Kt vt
Ft = Zt Pt Zt⊤ + Ht
Pt+1 = Tt+1 Pt L⊤ ⊤
t + Rt+1 Qt+1 Rt+1
Kt = Tt+1 Pt Zt⊤ Ft−1
vt = yt − Zt at − dt − βt xt
Lt = Tt+1 − Kt Zt
and the smoothed estimates, α̂t = E(αt | Yn ) and Vt = VAR(αt | Yn ) are given by
α̂t = at + Pt rt
Vt = Pt − Pt Nt Pt
rt = Zt⊤ Ft−1 vt + L⊤
t rt+1
Nt = Zt⊤ Ft−1 Zt + L⊤
t Nt+1 Lt
where rn+1 = 0 and Nn+1 = 0.
A.2. Univariate filter
Substantial computational gains are available using the univariate Kalman filter by avoiding
the inversion of the Ft matrix above (as well as enabling the use of the exact initial filter, see
below). For more details, see Durbin and Koopman (2012, §6.4).
When any Ht is non-diagonal, the observation equation is transformed by taking the LDL
factorization of the Ht matrices. The transformed parameters are marked with u to denote
that they have been transformed to be suitable for the univariate filter.
ytu = Ztu αt + dut + βtu xt + εut εut ∼ N (0, Htu )
ytu = Ct−1 yt Ztu = Ct−1 Zt dut = Ct−1 dt βtu = Ct−1 βt εut = Ct−1 εt
Ht = Ct Htu Ct⊤
28 MATLAB Toolbox for Mixed Frequency State Space Models
The multivariate filter and smoother is then modified by computing several quantities via the
univariate filter on the transformed system:
at,i+1 = at,i + Kt,i vt,i
⊤
Pt,i+1 = Pt,i − Kt,i Ft,i Kt,i
vt,i = yt,i − Zt,i at,i − dt,i − βt,i xt
⊤
Ft,i = Zt,i Pt,i Zt,i + Ht,i
⊤ −1
Kt,i = Pt,i Zt,i Ft,i
at+1,1 = Tt+1 at,p+1 + ct+1 + γt+1 wt+1
⊤ ⊤
Pt+1,1 = Tt+1 Pt,p+1 Tt+1 + Rt+1 Qt+1 Rt+1
where yi,t is the ith element of ytu , Zt,i is the ith row of Ztu , dt,i is the ith element of dut ,
βt,i is the ith row of βtu , and Ht,i is the ith diagonal element of Htu . The filtered estimates
of the state at = at,1 and Pt = Pt,1 are equivalent to those computed in the multivariate
−1
filter above. Note that Ft,i is a scalar, so that Ft,i is simply scalar division instead of more
computationally expensive matrix inversion. However, this comes at the cost of there being
no convenient transformation of the quantities vt,i , Ft,i or Kt,i that recovers the multivariate
versions computed above.
Similarly, the univariate smoother provides rt = rt,0 and Nt = Nt,0 :
⊤ −1
rt,i−1 = Zt,i Ft,i vt,i + L⊤
t,i rt,i
⊤ −1
Nt,i−1 = Zt,i Ft,i Zt,i + L⊤
t,i Nt,i Lt,i
Lt,i = Im − Kt,i Zt,i
rt−1,p = Tt⊤ rt,0
⊤
Nt−1,p = Tt−1 Nt,0 Tt−1
where rt+1,p = 0 and Nt+1,p = 0.
A.3. Exact initial Kalman filter
When the state αt is stationary, the initial values a1 and P1 can be computed as the uncondi-
tional mean and variance of the state given the system parameters by inverting the transition
equation.
To handle cases where some states are non-stationary, the state is separated into states with
known variance (those that are stationary) and those that are initialized as diffuse (the non-
stationary states). Let R̃ be a selection matrix with columns from the identity such that
the initial shock η1 is applied to the states with known variance. The selection matrix A is
composed of the columns of the identity matrix associated with the diffuse states such that
taking the limit as κ → ∞ allows the diffuse states to have infinite initial variance where the
initial values are given by
a1 = a + R̃η1 + Aδ η1 ∼ N (0, Q̃) δ ∼ N (0, κI)
⊤
P1 = P∗,1 + κP∞,1 P∗,1 = R̃Q̃R̃ P∞,1 = AA⊤
Journal of Statistical Software 29
The unconditional mean, a, and variance, P∗,1 , of the state are computed by inverting the
stationary portion of the system. Any elements of a associated with non-stationary or diffuse
states are set to 0. Letting T̃ = R̃⊤ T1 R̃ and c̃ = R̃⊤ c1 , this is accomplished by
R̃⊤ a = (Im − T̃ )−1 c̃
vec(R̃Q̃R̃⊤ ) = (Im2 − T̃ ⊗ T̃ )−1 vec(R̃⊤ R1 Q1 R1⊤ R̃)
When all states are stationary, this initialization collapses down to the simple case where
α1 ∼ N (a1 , P1 ), where a1 and P1 are the unconditional mean and variance of the state,
determined by inverting the full system.
Given this initialization, the univariate filter recursions must be altered to separate the states
with finite and infinite variances (see Koopman and Durbin 2000; Durbin and Koopman 2012,
§5.2):
⊤ ⊤
F∗,t,i = Zt,i P∗,t,i Zt,i + Ht,i F∞,t,i = Zt,i P∞,t,i Zt,i
⊤ −1 ⊤ −1
K∗,t,i = P∗,t,i Zt,i F∗,t,i K∞,t,i = P∞,t,i Zt,i F∞,t,i
at,i + K∗,t,i vt,i F∞,t,i = 0
(
at,i+1 =
at,i + K∞,t,i vt,i F∞,t,i ̸= 0
F∞,t,i = 0
(
⊤ F
P∗,t,i − K∗,t,i K∗,t,i ∗,t,i
P∗,t,i+1 =
P∗,t,i − (K∗,t,i K∞,t,i + K∞,t,i K∗,t,i
⊤ ⊤ −K
∞,t,i K∞,t,i )F∗,t,i
⊤ F∞,t,i ̸= 0
F∞,t,i = 0
(
P∞,t,i
P∞,t,i+1 = ⊤
P∞,t,i − K∞,t,i K∞,t,i F∞,t,i F∞,t,i ̸= 0
at+1,1 = Tt+1 at,p+1 + ct+1 + γt+1 wt+1
⊤ ⊤ ⊤
P∞,t+1,1 = Tt+1 P∞,t,p+1 Tt+1 P∗,t+1,1 = Tt+1 P∗,t,p+1 Tt+1 + Rt+1 Qt+1 Rt+1
For any set of system parameters where the state can be identified, there exists some time d
such that F∞,d,i = 0 for all i. For time t ≥ d, the simpler Kalman filter recursion above can
be employed with Pt,i = P∗,t,i .
The smoother must be similarly altered so that beginning at t = d, the computation of rt,i is
expanded to account for the initialization:
L∞,t,i = Im − K∞,t,i Zt,i L∗,t,i = Im − K∗,t,i Zt,i
(0) −1
Lt,i = (K∞,t,i − K∗,t,i )Zt,i F∗,t,i F∞,t,i
Z ⊤ F −1 v + L⊤ r (0) F∞,t,i = 0
(0) t,i ∗,t,i t,i ∗,t,i t,i
rt,i−1 =
L⊤ r (0) F∞,t,i ̸= 0
∞,t,i t,i
r (1) F∞,t,i = 0
(1)
rt,i−1 = t,i
Z ⊤ F −1 vt,i + L(0)⊤ r (0) + L⊤ r (1) F∞,t,i ̸= 0
t,i ∞,t,i t,i t,i ∞,t,i t,i
(0) (0) (1) (1)
rt−1,p = Tt⊤ rt,0 rt−1,p = Tt⊤ rt,0
30 MATLAB Toolbox for Mixed Frequency State Space Models
(0) (1)
α̂t = at + P∗,t,1 rt,0 + P∞,t,1 rt,0
(0) (1) (0)
where rd,p = rd,p and rd,p = 0. Note that L∞,t,i and Lt,i only need to be computed when
F∞,t,i ̸= 0 and L∗,t,i only needs to be computed when F∞,t,i = 0.
For the smoothed variance of the state,
⊤
(0) (1) (1) (2)
Vt = P∗,t,1 − P∗,t,1 Nt,0 P∗,t,1 − P∞,t,1 Nt,0 P∗,t,1 − P∞,t,1 Nt,0 P∗,t,1 − P∞,t,1 Nt,0 P∞,t,1
(0) (0) (1) (1) (2) (1)
Nt−1,p = Tt⊤ Nt,0 Tt Nt−1,p = Tt⊤ Nt,0 Tt Nt−1,p = Tt⊤ Nt,0 Tt
where when F∞,t,i = 0,
(0) ⊤ −1 (0)
Nt,i−1 = Zt,i F∗,t,i Zt,i + L⊤
∗,t,i Nt,i L∗,t,i
(1) (1)
Nt,i−1 = Nt,i
(2) (2)
Nt,i−1 = Nt,i
and when F∞,t,i ̸= 0,
(0) (0)
Nt,i−1 = L⊤
∞,t,i Nt,i L∞,t,i
(1) ∗⊤ −1 (0) (0) (1)
Nt,i−1 = Zt,i F∞,t,i Zt,i + L⊤ ⊤
∞,t,i Nt,i Lt,i + L∞,t,i Nt,i L∞,t,i
(2) ⊤ −2 (0)⊤ (1) (0)
Nt,i−1 = Zt,i F∞,t,i Zt,i F∗,t,i + Lt,i Nt,i Lt,i
(1) (0) (0)⊤ (1) (2)
+ L⊤ ⊤
∞,t,i Nt,i Lt,i + Lt,i Nt,i L∞,t,i + L∞,t,i Nt,i L∞,t,i
(0) (1) (2)
where Nd,p = Nd,p and Nd,p = Nd,p = 0.
B. Parameter estimation
B.1. General maximum likelihood estimation
The likelihood of of data y1 , . . . , yn in the standard filter as shown in Durbin and Koopman
(2012, §7.2) is given by the prediction error decomposition:
n
np 1X −1
log L(Yn ) = − log 2π − log|Ft | + vt⊤ Ft,i vt
2 2 t=1
For the univariate filter, the same decomposition works in the univariate context:
p
np 1 n
log L(Yn ) = − log 2π − log Ft,i + vt,i
XX
2
/Ft,i
2 2 t=1 i=1
The likelihood for the exact initial filter allows for some simplifications in F∞,t,i :
n X p d X n Xp p
1X 1X 1X
log Ld (Yn ) = − ιt,i log 2π − wt,i − ιt,i log Ft,i + vt,i
2
/Ft,i
2 t=1 i=1 2 t=1 i=1 2 t=d i=1
Journal of Statistical Software 31
where ιt,i = 1 if F∗,t,i ̸= 0 or t > d, and
(0)2
ιt,i (log(F∗,t,i ) + vt,i /F∗,t,i ) F∞,t,i = 0
(
wt,i =
log F∞,t,i F∞,t,i ̸= 0
Since these quantities are naturally produced by the Kalman filter, this is the preferred
method to calculate the likelihood.
For parameter estimation, the elements of the parameter matrices Zt , dt , βt , Ht , Tt , ct , γt ,
Rt , and Qt are divided into those that are known and those that are to be estimated as a
function of the parameter vector θ. We are interested in the maximum likelihood estimate of
those parameters given data Yn = {y1 , . . . , yn }.
We are interested in the maximum likelihood estimate of a set of structural parameters (θ)
given data y1 , . . . , yn where elements of the state space parameters (Zt , dt , βt , Ht , Tt , ct , γt , Rt ,
and Qt ) depend on θ. For simplicity, the state space parameters will be restricted such that
each of their scalar elements must be a function of a single element of a vector of the reduced
form parameters (ψ) which depend on θ. The reduced form parameters are unrestricted in
how they may depend on θ, allowing for a rich specification of parameter restrictions. In
many estimations, there will be the same number of structural parameters and reduced form
parameters.
More formally, let θU ∈ Rnθ and ψ ∈ Rnψ where nθ ≤ nψ . Define a set of bounds θ and θ̄
¯
such that each element θi of θ ∈ Rnθ , θi ≤ θi ≤ θ̄i . Define a set of functions Ψi : [θ, θ̄] → R
¯ ¯
such that each element of ψ is a function of the structural parameters, ψi = Ψi (θ). Except in
cases of cross-parameter restrictions, nθ = nψ and Ψi (θ) = θi so that ψ = θ.
Additionally define a set of functions τXi,j : R → R for X ∈ {Zt , dt , βt , Ht , Tt , ct , γt , Rt , Qt }
such that each element of the parameter matrices to be estimated is a transformation of an
element of ψ. Common τXi,j transformations include the identity, exponential, negative ex-
ponential, and logistic transformations to allow for an unbounded ψ ∈ Rnψ while maintaining
bounds on the parameter matrices. When estimating models with mixed frequency data, all
state space parameters shared between the high-frequency and low-frequency states depend
on the same reduced form parameters simply by using different τXi,j functions. In almost all
cases, the specification of the τXi,j functions will be done automatically to account for state
space parameter bounds and accumulator definitions. User definitions of τXi,j should be rare.
With the likelihood of an unconstrained parameter vector θU defined, the optimization can be
performed using a variety of methods. For simple models, gradient ascent is usually sufficient.
In practice, it has been observed that further improvements in the likelihood are sometimes
possible using simplex based optimization methods once the gradient ascent method has
converged. Due to this, the default behavior in MFSS is to repeat these two methods until
neither is able to improve the likelihood using standard convergence criteria. Alternative
optimization routines are also supported in MFSS including simulated annealing and particle
swarm optimizations, though we have found these to be less useful than gradient ascent and
simplex optimization.
It has been our experience that initial values for the parameters to be estimated can have
meaningful influence in achieving optima in maximum likelihood estimation. In that regard,
it is important to provide reasonable starting values for the optimization, as was done in
Sections 4.1 and 4.5. These initial values are usually most easily generated from approximate
32 MATLAB Toolbox for Mixed Frequency State Space Models
models with low-frequency data that avoids the latent states that result in the mixed fre-
quency case, and being mindful of the type of accumulation being modeled (e.g., a simple or
triangle average). We should also note that the common hurdles towards maximum likelihood
estimation including badly scaled variables, and models that are not identified are likely to
lead to convergence issues. In the case of dramatically different scales amongst variables,
standardizing the variables before estimation can be advantageous. Ensuring your model is
identified is always an important first step before estimation is attempted.
B.2. VAR estimation
For VAR models, estimation via the EM algorithm of Shumway and Stoffer (1982) is compu-
tationally more convenient. To do so, we need to compute Jt = Cov(αt , αt+1 | Yn ) as shown
in Durbin and Koopman (2012, §4.7),
Jt = Pt L†⊤
t (Im − Nt+1 Pt+1 )
where L†t is defined in C.1. Using this and quantities from A, the algorithm requires iterating
between the two steps:
1. Calculate α̂ti+1 , Vti+1 , Jti+1 using parameters T i , ci , Qi and ai0 .
2. Calculate ai+1
0 = α̂0i and estimate state space parameters T i+1 , ci+1 , Qi+1 given the
estimates of the state α̂ti+1 , Vti+1 , Jti+1 according to a modified OLS estimate.
This method avoids the computation of gradients of the likelihood at the expense of computing
the variance and covariance of the state. As a result, it scales much better as the size of the
state increases. While it does require a fixed value for the variance of the initial state, if P0
needs to be estimated the optimal values from the EM algorithm can be used as a starting
point for general maximum likelihood estimation.
C. State decompositions
Since all of the Kalman filter and smoother calculations are linear, we can decompose the
estimated states by the effects of the data and parameters on the state in each period. In
general, we will be creating a four-dimensional object since each state k ∈ {1, . . . , m} at
time t is affected by data series i ∈ {1, . . . , p} from time j. Depending on the use of the
decompositions, this 4D object will be collapsed in different ways. For example, if we were
suspicious of the validity of an observation, we could inspect how much it affects a state
variable of interest, in which case we could select a particular j, k, and i and plot across t.
More commonly, we will be interested in plots of how a given state variable is determined by
the data, in which case we sum across the origin dates of the data j and plot the effect of
each data series on a given state i across time t.
C.1. Filtered state decompositions
For the filter, we are interested in
t−1
at = 1 1 1ωtjc + 1ωtjd + 1ωtjx + 1ωtjw + 1ωtj
ωta0 + +
X
c
ωtt
j=1
Journal of Statistical Software 33
where
• ωta0 is an (m × m) matrix of the effect of the initial conditions on at ,
• ωtj
c is an (m × m) matrix of the effect of c on a ,
j t
• ωtj
d is an (m × p) matrix of the effect of d on a ,
j t
• ωtj
x is an (m × k) matrix of the effect of x on a ,
j t
• ωtj
w is an (m × l) matrix of the effect of w on a , and
j t
• ωtj is an (m × p) matrix of the effect of yj on at .
In each case, the period of the state being affected is denoted by t and the period of the
observation of parameter causing the effect is denoted by j.
We begin by incorporating the observed data from a single period. Using the filtering equa-
tions from above, we examine the integration of information contained in yt on at+1 = at+1,1
given at,1 :
at,2 = at,1 + Kt,1 vt,1
= at,1 + Kt,1 (yt,1 − Zt,1 at,1 − dt,1 − βt,1 xt )
= Lt,1 at,1 + Kt,1 (yt,1 − dt,1 − βt,1 xt )
at,3 = at,2 + Kt,2 vt,2
= at,2 + Kt,2 (yt,2 − Zt,2 at,2 − dt,2 − βt,2 xt )
= Lt,2 (Lt,1 at,1 + Kt,1 (yt,1 − dt,1 − βt,1 xt )) + Kt,2 (yt,2 − dt,2 − βt,2 xt )
= Lt,2 Lt,1 at,1 + Kt,2 (yt,2 − dt,2 − βt,2 xt ) + Lt,2 Kt,1 (yt,1 − dt,1 − βt,1 xt )
..
.
1 p i+1
at,p+1 = Lt,k at,1 + Lt,k Kt,i (yt,i − dt,i − βt,i xt )
Y X Y
k=p i=1 k=p
at+1,1 = Tt+1 at,p+1 + ct+1 + γt+1 wt+1
1 p i+1
= Tt+1 Lt,k at,1 + ct+1 + γt+1 wt+1 + Tt+1
Y X Y
Lt,k Kt,i yt,i
k=p i=1 k=p
p
X i+1
Y p
X i+1
Y
− Tt+1 Lt,k Kt,i dt,i − Tt+1 Lt,k Kt,i βt,i xt
i=1 k=p i=1 k=p
= L†t at,1 + ct+1 + γt+1 wt+1 + Kt† yt − Kt† dt − Kt† βt xt
where the products over Lt,k proceed from higher to lower values of k and
1
L†t = Tt+1
Y
Lt,i
i=p
hQ i
Kt† = Tt+1 2
i=p Lt,i Kt,1 . . . Lt,p Kt,p−1 Kt,p
34 MATLAB Toolbox for Mixed Frequency State Space Models
To accommodate the exact initialization, when t ≤ d we adjust these expressions for the
condition on F∞,t,i by creating selection objects
1
L†t = Tt+1 ∗
Y
St,i
i=p
hQ i
Kt† = Tt+1 2
i=p Lt,i K̃t,1 . . . Lt,p K̃t,p−1 K̃t,p
F∞,t,i = 0 F∞,t,i = 0
( (
∗ L∗,t,i K∗,t,i
St,i = K̃t,i =
L∞,t,i F∞,t,i ̸= 0 K∞,t,i F∞,t,i ̸= 0
which allows the earlier recursion on at,1 to hold.
To find how information propagates across time, we use this new formulation of the filter,
starting with the first period:
a1 =L†0 a0 + c1 + γ1 w1
a2 =L†1 a1 + c2 + γ2 w2 + K1† y1 − K1† d1 − K1† β1 x1
=L†1 (L†0 a0 + c1 + γ1 w1 ) + c2 + γ2 w2 + K1† y1 − K1† d1 − K1† β1 x1
=L†1 L†0 a0 + L†1 c1 + L†1 γ1 w1 + c2 + γ2 w2 + K1† y1 − K1† d1 − K1† β1 x1
a3 =L†2 a2 + c3 + γ3 w3 + K2† y2 − K2† d2 − K2† β2 x2
=L†2 L†1 L†0 a0 + L†2 L†1 c1 + L†2 c2 + c3 + L†2 L†1 γ1 w1 + L†2 γ2 w2 + γ3 w3 +
L†2 K1† y1 + K2† y2 − L†2 K1† d1 − K2† d2 − L†2 K1† β1 x1 − K2† β2 x2
From which the recursion becomes clear and we can see that information propagates with the
product of the L†t matrix.
When using the univariate treatment of multivariate series, each series element of yt can
potentially affect other elements of ytu . To account for this, we could repeatedly consider the
effects of each observation where each element of yt is transformed and considered in turn.
To do so in a single operation, we can simply consider the product Ct−1 diag(yt ), which gives
us the full contribution from each element of yt across all transformed series ytu .
From the above expressions, we can infer that the weights are
j+1 j+1
ωtj = L∗k Kj∗ Cj−1 diag(yj ) = − L∗k Kj∗ Cj−1 βj diag(xj )
Y Y
x
ωtj
k=t−1 k=t−1
j+1 j
= − L∗k Kj∗ Cj−1 diag(dj ) = L∗k γj diag(wj )
Y Y
d w
ωtj ωtj
k=t−1 k=t−1
j 0
= L∗k diag(cj ) ωta0 = L∗k diag(a0 )
Y Y
c
ωtj
k=t−1 k=t−1
completing the earlier desired decomposition of the filtered state.
Journal of Statistical Software 35
C.2. Smoothed state decompositions
For the smoother, we similarly want to compute quantities {ω̂tj , ω̂tc , ω̂td , ω̂tx , ω̂tw , ω̂ta0 } such that
T
α̂t = 1ω̂ta0 + 1ω̂tj + 1ω̂tjc + 1ω̂tjd + 1ω̂tx + 1ω̂tw
X
j=1
We first begin by breaking this expression into contributions from at and rt ,
T
r a0 rc rd rx rw r
rt = 1 ω t + 1 ωtj +1 ωtj +1 ωtj +1 ωtj +1 ωtj
X
j=t
α̂t = at + Pt rt
T
r a0 rc rw
= 1ωta0 + 1Pt ω t + 1ωtjc + 1Pt ωtj +1ωtjw + 1Pt ωtj
X
j=1
rd rx r
+1 d
ωtj + 1
Pt ω tj + 1 x
ωtj + 1
Pt ω tj + 1ωtj + 1Pt ωtj
T
= 1ω̂ta0 + 1ω̂tjc + 1ω̂tjw + 1ω̂tjd + 1ω̂tjx + 1ω̂tj
X
j=1
where
r rd rc
ω̂tj = ωtj + Pt ω tj d
ω̂tj = ωtj
d
+ Pt ω tj c
ω̂tj = ωtj
c
+ Pt ω tj
rx rw r a0
x
ω̂tj = ωtj
x
+ Pt ω tj w
ω̂tj = ωtj
w
+ Pt ω tj ω̂ta0 = ωta0 + Pt ω t
Before examining the smoother recursion, we first rewrite the observation errors in matrix
form,
yt,1 − Zt,1 at,1 − dt,1 − βt,1 xt
y − Z a − d − β x
t,2 t,2 t,2 t,2 t,2 t
vtu = yt,3 − Zt,3 at,3 − dt,3 − βt,3 xt
..
.
yt,p − Zt,p at,p − dt,p − βt,p xt
Zt,1 at,1
Zt,2 [Lt,1 at,1 + Kt,1 (yt,1 − dt,1 − βt,1 xt )]
Zt,3 [Lt,2 Lt,1 at,1 + Lt,2 Kt,1 (yt,1 − dt,1 − βt,1 xt ) + Kt,2 (yt,2 − dt,2 − βt,2 xt )]
= yt − dt −
u u
..
.
hQ Pp−1 Qi+1 i
at,1 + k=p−1 Lt,k Kt,i (yt,i − dt,i − βt,i xt )
1
Zt,p k=p−1 Lt,k i=1
y u
= yt − dt − βt xt − Ãt (yt − dt − βt xt ) − Aat at
u u u u u
= Ayt Ct−1 (yt − dt − βt xt ) − Aat at
36 MATLAB Toolbox for Mixed Frequency State Space Models
where Ayt = Ip − Ãyt and
0 0 ... 0 0 0
Zt,2 Kt,1 0 ... 0 0 0
y
Zt,3 Lt,2 Kt,1 Zt,3 Kt,2 ... 0 0 0
Ãt = .. .. .. .. ..
. . . . .
Zt,p−1 ( 2k=p−2 Lt,k )Kt,1 Zt,p−1 ( 3k=p−2 Lt,k )Kt,2 . . . 0 0
Q Q
Zt,p−1 Kt,p−2
Zt,p ( k=p−1 Lt,k )Kt,1 Zt,p ( 3k=p−1 Lt,k )Kt,2 . . . Zt,p Lt,p−1 Kp−2 Zt,p Kt,p−1 0
Q2 Q
Zt,1
Zt,2 Lt,1
Aat = ..
Q .
1
Zt,p k=p−1 Lt,k
To adjust for the diffuse filter, this must be altered so that the lower-diagonal elements of Ãyt
are adjusted depending on the condition on F∞,t,i and the propagation of at,0 respects the
selection between L∗,t,i and L∞,t,i :
j
A˜yt (i,j) = Zt,i ∗
Y
St,k K̃t,j
k=i−1
Zt,1
∗
Zt,2 St,2
t =
Aa∗ ..
Q .
1 ∗
Zt,p k=p−1 St,k
Similarly to how at was handled for the filter, the incorporation of the information at time t
to form rt,0 given rt+1,0 can be written as a single operation in matrix for as
⊤
rt,p = Tt+1 rt+1,0
⊤ −1
rt,p−1 = Zt,p Ft,p vt,p + L⊤ ⊤
t,p Tt+1 rt+1,0
⊤ −1
rt,p−2 = Zt,p−1 Ft,p−1 vt,p−1 + L⊤
t,p−1 rt,p−1
⊤ −1 ⊤ −1
= Zt,p−1 Ft,p−1 vt,p−1 + L⊤ ⊤ ⊤ ⊤
t,p−1 Zt,p Ft,p vt,p + Lt,p−1 Lt,p Tt+1 rt+1,0
..
.
rt,0 = Mt† vtu + L†⊤
t rt+1,0
= Mt† Ayt Ct−1 yt − Mt† Ayt Ct−1 dt − Mt† Ayt Ct−1 βt xt − Mt† Aat at + L†⊤
t rt+1,0
where
h Q i
Mt† = Zt,1
⊤ F −1 L⊤ Z ⊤ F −1 . . .
t,1 t,1 t,2 t,2
p−1 ⊤
L
i=1 t,i Z ⊤ F −1
t,p t,p
h i
−1 −1 −1
= Aa⊤
t diag Ft,1 Ft,2 ... Ft,p
Journal of Statistical Software 37
Noting the similarity between propagation through time like at , we have that
r r
n o n o
† a
ω tj = L†⊤
t ω t+1,j − Mt At ωtj + Mt† Ayt Ct−1 diag(yt )
t>j t=j
rd †⊤ r d
n o n o
† a d † y −1
ω tj = Lt ω t+1,j − Mt At ωtj − Mt At Ct diag(dt )
t>j t=j
rx †⊤ r x
n o n o
† a x † y −1
ω tj = Lt ω t+1,j − Mt At ωtj − Mt At Ct βt diag(xt )
t>j t=j
rc †⊤ r c
n o
† a c
ω tj = Lt ω t+1,j − Mt At ωtj
t≥j
rw †⊤ r w
n o
† a w
ω tj = Lt ω t+1,j − Mt At ωtj
t≥j
r a0 †⊤ r a0 † a a0
ω t = Lt ω t+1 −Mt At ωt
To accommodate the exact initialization, we expand rt as as done in the smoother:
(0) (1)
α̂t = at + P∗,t,1 rt,0 + P∞,t,1 rt,0
where
Z ⊤ F −1 vt,i + L⊤ r (0) F∞,t,i = 0
(0) t,i ∗,t,i ∗,t,i t,i
rt,i−1 =
L⊤ r (0) F∞,t,i ̸= 0
∞,t,i t,i
(1)
(1)
r
t,i F∞,t,i = 0
rt,i−1 =
Z ⊤ F −1 v + L⊤ r (1) + L(0)⊤ r (0) F∞,t,i ̸= 0
t,i ∞,t,i t,i ∞,t,i t,i t,i t,i
To handle the selection in these equations based on the condition F∞,t,i = 0, define a set of
selection quantities:
F∞,t,i = 0 F∞,t,i = 0 0m×m F∞,t,i = 0
( (
∗ L∗,t,i ∞ Im (0)
St,i = St,i = St,i =
L∞,t,i F∞,t,i ̸= 0 L∞,t,i F∞,t,i ̸= 0 L(0)
t,i F∞,t,i ̸= 0
0m×1 F∞,t,i = 0
( ⊤ −1
F∞,t,i = 0
(
∗
Zt,1 F∗,t,1 ∞
Mt,i = Mt,i = ⊤ −1
0m×1 F∞,t,i ̸= 0 Zt,1 F∞,t,1 F∞,t,i ̸= 0
This allows the recursion of the diffuse smoother to be written as
(0) ∗ ∗⊤ (0)
rt,i−1 = Mt,i vt,i + St,i rt,i
(1) ∞ ∞⊤ (1) (0)⊤ (0)
rt,i−1 = Mt,i vt,i + St,i rt,i + St,i rt,i
(0)
Using these quantities, we can rewrite the recursion for rt,0 as we did rt,0 . To arrive at the
(1)
similar recursion for rt,0 , first make the simplifications that were done for rt,0 :
p i−1
!
(1) (1) (0)⊤ (0)
rt,0 = M̃t∞ vtu + L∞⊤ + ∞⊤
X Y
t rt+1,0 St,k St,i rt,i
i=1 k=1
38 MATLAB Toolbox for Mixed Frequency State Space Models
Also note that we can infer from the manipulation of rt,0 earlier that
p p j−1
rt,i = L⊤ ⊤ L⊤ ⊤ −1
t,k Tt+1 rt+1,0 +
Y X Y
t,k Zt,j Ft,j vt,j
k=i+1 j=i+1 k=i+1
Making the required notation adjustments, these two equations give the recursion
(1) (1)
rt,0 = M̃t∞ vtu + L∞⊤
t rt+1,0
p i−1
! p
∞⊤ (0)⊤ ∗⊤ ⊤ (0)
+
X Y Y
St,k St,i St,k Tt+1 rt+1,0
i=1 k=1 k=i+1
p i−1
! p j−1
∞⊤ (0)⊤
+ S ∗⊤ M ∗ vt,j
X Y X Y
St,k S t,i
t,k t,j
i=1 k=1 j=i+1 k=i+1
(0) (0)⊤ (0) (1)
= (M̃t∞ + Mt )vtu + Lt rt+1,0 + L∞⊤
t rt+1,0
This allows us to write the recursions required for the diffuse smoother as
(0) (0)
rt,0 = Mt∗ vtu + L∗⊤
t rt+1,0
(1) (0)⊤ (0) (1)
rt,0 = Mt∞ vtu + Lt rt+1,0 + L∞⊤
t rt+1,0
1
∗
L∗t = Tt+1
Y
St,i
i=p
1
L∞
t = Tt+1
∞
Y
St,i
i=p
⊤
p i−1
! p
(0) ∞⊤ (0)⊤ ∗⊤
L = Tt+1
X Y Y
t
St,k St,i St,k
i=1 k=1 k=i+1
h Q i
p−1 ∗⊤
Mt∗ = Mt,1
∗ ∗⊤ M ∗
St,1 t,2 . . . i=1 St,i
∗
Mt,p
(0)
Mt∞ = M̃t∞ + Mt
h Q i
p−1 ∞⊤
M̃t∞ = Mt,1
∞ S ∞⊤ M ∞ . . .
t,1 t,2 i=1 St,i
∞
Mt,p
p i−1
!
(0) ∞⊤ (0)⊤ (0)
Mt =
X Y
St,k St,i Mt,i
i=1 k=1
h Q i
(0) p−1
Mt,i = 0 ... 0 Mt,i+1
∗ ∗⊤ M ∗
St,i+1 t,i+2 . . .
∗⊤
k=i+1 St,k
∗
Mt,p
Combined with the expressions above gives the recursion
(0) (0)
rt,0 = Mt∗ Ayt Ct−1 yt − Mt∗ Ayt Ct−1 dt − Mt∗ Aat at + L∗⊤
t rt+1,0
(1) (0)⊤ (0) (1)
rt,0 = Mt∞ Ayt Ct−1 yt − Mt∞ Ayt Ct−1 dt − Mt∞ Aat at + Lt rt+1,0 + L∞⊤
t rt+1,0
Journal of Statistical Software 39
from which the weights for the diffuse smoother are similar to above
r(0) r(0)
n o
ω tj = L∗⊤
t ω t+1,j − {Mt∗ Aat ωtj }t>j + Mt∗ Ayt Ct−1 diag(yt )
t=j
d d
r(0) r(0)
n o n o
ω tj = L∗⊤
t ω t+1,j − Mt∗ Aat ωtj
d
− Mt∗ Ayt Ct−1 diag(dt )
t>j t=j
x x
r(0) r(0)
n o n o
ω tj = L∗⊤
t ω t+1,j − Mt∗ Aat ωtj
x
− Mt∗ Ayt Ct−1 βt diag(xt )
t>j t=j
c c
r(0) r(0)
n o
ω tj = L∗⊤
t ω t+1,j − Mt∗ Aat ωtj
c
t≥j
w w
r(0) r(0)
n o
ω tj = L∗⊤
t ω t+1,j − Mt∗ Aat ωtj
w
t≥j
a0 a0
r(0) r(0)
ω tj = L∗⊤
t ω t+1 −Mt∗ Aat ωta0
r(1) r(1) (0) r
(0)
n o
ω tj = L∞⊤
t ω t+1,j − {Mt∞ Aat ωtj }t>j + Mt∞ Ayt Ct−1 diag(yt ) + Lt ω t+1,j
t=j
d d (0) d
r(1) r(1) (0) r
n o n o
ω tj = L∞⊤
t ω t+1,j − Mt∞ Aat ωtj
d
− Mt∞ Ayt Ct−1 diag(dt ) + Lt ω t+1,j
t>j t=j
r(1)
x
r(1)
x (0) x
(0) r
n o n o
ω tj = L∞⊤
t ω t+1,j − Mt∞ Aat ωtj
x
− Mt∞ Ayt Ct−1 βt diag(xt ) + Lt ω t+1,j
t>j t=j
c c c
r(1) r(1) r(0)
n o
(0)
ω tj = L∞⊤
t ω t+1,j − Mt∞ Aat ωtj
c
+ Lt ω t+1,j
t≥j
w w (0) w
r(1) r(1) (0) r
n o
ω tj = L∞⊤
t ω t+1,j − Mt∞ Aat ωtj
w
+ Lt ω t+1,j
t≥j
a0 a0 a0
r(1) r(1) r(0)
ω t = L∞⊤
t ω t+1 −Mt∞ Aat ωta0 + L(0)
t ω t+1
These can then be used similarly to those above to find weights for αt for t ≤ d:
d d
r(0) r(1) r(0) r(1)
ω̂tj = ωtj + Pt∗ ω tj +Pt∞ ω tj d
ω̂tj = ωtj
d
+ Pt∗ ω tj +Pt∞ ω tj
x x c c
r(0) r(1) r(0) r(1)
x
ω̂tj = ωtj
x
+ Pt∗ ω tj +Pt∞ ω tj c
ω̂tj = ωtj
c
+ Pt∗ ω tj +Pt∞ ω tj
w w a0 a0
r(0) r(1) r(0) r(1)
w
ω̂tj = ωtj
w
+ Pt∗ ω tj +Pt∞ ω tj ω̂ta0 = ωta0 + Pt∗ ω t +Pt∞ ω t
C.3. Using the weights
Once the filter and smoother weights have been computed, we can use them in several ways
to summarize the parameters of our model. Most directly comparable to the coefficients in a
traditional state space context, we can create a (p × m) matrix w that summarizes the total
absolute contributions for each state from all of the observed data:
T X
T
w=
X
|ω̂tj |
t=1 j=1
This matrix is often more informative after it has been row normalized to sum to one so that
each element gives the proportion of a given state attributable to a time series.
40 MATLAB Toolbox for Mixed Frequency State Space Models
Additionally, plots of the contributions to the state can help diagnose what data are in-
formative about them. In particular, for dynamic factor models we may be interested in
understanding how a given state is affected by different time series so that decomposing the
state at time t by all influences from the p time series across all observation periods j is
informative. Plots of the cumulative effects on the state are also often informative, especially
when the scale of the time series differs substantially.
Affiliation:
Scott A. Brave
Morning Consult
Chicago, IL, United States of America
E-mail: sbrave@morningconsult.com
R. Andrew Butters
Business Economics and Public Policy, Kelley School of Business, Indiana University
1309 E. Tenth Street, Bloomington, IN, 47405, United States of America
E-mail: rabutter@indiana.edu
David Kelley
Federal Reserve Bank of New York
230 S. LaSalle St., Chicago, IL, 60604, United States of America
E-mail: david.kelley@ny.frb.org
Journal of Statistical Software https://www.jstatsoft.org/
published by the Foundation for Open Access Statistics https://www.foastat.org/
October 2022, Volume 104, Issue 10 Submitted: 2018-11-09
doi:10.18637/jss.v104.i10 Accepted: 2022-10-13
File and source
- File
- brave-butters-kelley-2022-jss-104-10-mixed-frequency-state-space-models.pdf
- Size
- 1,016,668 bytes
- SHA-256
- 1dcbd2380623861d712f318a808ab3fe47bd0e834a13d75dacf6af4c46aa4122
- Original
- www.jstatsoft.org