The full text on this page is automatically extracted from the file linked above and may contain errors and inconsistencies.
FEDERAL RESERVE BANK o f ATLANTA
WORKING PAPER SERIES
Perturbation Methods for Markov-Switching DSGE Models
Andrew Foerster, Juan Rubio-Ramírez, Daniel F. Waggoner, and Tao Zha
Working Paper 2013-1
March 2013
Abstract: This paper develops a general perturbation methodology for constructing high-order
approximations to the solutions of Markov-switching DSGE models. We introduce an important and
practical idea of partitioning the Markov-switching parameter space so that a steady state is well
defined. With this definition, we show that the problem of finding an approximation of any order can
be reduced to solving a system of quadratic equations. We propose using the theory of Grobner bases
in searching all the solutions to the quadratic system. This approach allows us to obtain all the
approximations and ascertain how many of them are stable. Our methodology is applied to three
models to illustrate its feasibility and practicality.
JEL classification: C6, E1
Key words: Markov-switching parameters, partition, higher order approximations, no certainty equivalence,
quadratic system, Grobner bases
The authors thank Leonardo Melosi, Seonghoon Cho, Rhys Bidder, seminar participants at Duke University, the Federal
Reserve Bank of St. Louis, the 2010 Society of Economic Dynamics meetings, the 2011 Federal Reserve System Committee
on Business and Financial Analysis Conference, the 2012 Annual Meeting of the American Economic Association, the 8th
Dynare Conference, and the 2012 National Bureau of Economic Research (NBER) Workshop on Methods and Applications
for DSGE Models for helpful comments. This research is supported in part by the National Science Foundation Grant Nos.
SES-1127665 and SES-1227397. The views expressed here are the authors’ and not necessarily those of the Federal Reserve
Banks of Atlanta or Kansas City, the Federal Reserve System, or the NBER. Any remaining errors are the authors’
responsibility.
Please address questions regarding content to Andrew Foerster, Federal Reserve Bank of Kansas City, 1 Memorial Drive,
Kansas City, MO 64198, 816-881-2751, andrew.foerster@kc.frb.org; Juan Rubio-Ramírez, Duke University, Federal Reserve
Bank of Atlanta, CEPR, FEDEA, and BBVA Research, P.O. Box 90097, Durham, NC 27708, 919-660-1865, juan.rubioramirez@duke.edu; Daniel F. Waggoner, Federal Reserve Bank of Atlanta, 1000 Peachtree Street, N.E., Atlanta, GA 303094470, dwaggoner@frbatlanta.org; or Tao Zha, Federal Reserve Bank of Atlanta, Emory University, and NBER, 1000
Peachtree Street, N.E., Atlanta, GA 30309-4470, zmail@tzha.net.
Federal Reserve Bank of Atlanta working papers, including revised versions, are available on the Atlanta Fed’s website at
frbatlanta.org/pubs/WP/. Use the WebScriber Service at frbatlanta.org to receive e-mail notifications about new papers.
Perturbation Methods
for Markov-Switching DSGE Models
1
Introduction
In this paper we show how to use perturbation methods as described in Judd (1998) and
Schmitt-Grohe and Uribe (2004) to solve Markov-switching dynamic stochastic general equilibrium (MSDSGE) models. Our contribution advances the current literature in two significant
respects. First, we develop a general methodology for approximating the solution of a larger
class of Markov-switching models than currently possible and, second, we show the feasibility
and practicality of implementing our methodology when we consider high-order approximations
to the model solution. Current methods only allow for first-order approximations.
The literature on Markov-switching linear rational expectations (MSLRE) models has been
an active field in empirical macroeconomics (Leeper and Zha (2003), Blake and Zampolli (2006),
Svensson and Williams (2007), Davig and Leeper (2007), and Farmer et al. (2009)). Building
on standard linear rational expectations models, the MSLRE approach allows model parameters to change over time according to discrete Markov chain processes. This nonlinearity has
proven to be important in explaining changes in monetary policy and macroeconomic time series (Schorfheide (2005), Davig and Doh (2008), Liu et al. (2011), and Bianchi (2010)) and in
modeling the expected effects of future fiscal policy changes (Davig et al. (2010), Davig et al.
(2011), Bi and Traum (2012)). In particular, Markov-switching models provide a tractable way
to study how agents form expectations over possible discrete changes in the economy, such as
those in technology and policy.
There are, however, two main shortcomings with the MSLRE approach. First, that approach
begins with a system of standard linear rational expectations equations that have been obtained
from linearizing equilibrium conditions as though the parameters were constant over time. Discrete Markov chain processes are then annexed to certain parameters. As a consequence, the
resultant MSLRE model may be incompatible with the optimizing behavior of agents in an
1
original economic model with Markov-switching parameters. Second, because it builds on linear
rational expectations models, the MSLRE approach does not take into account higher-order
coefficients in the approximation. Higher-order approximations improve the approximation accuracy and may be potentially important for certain economic questions.
This paper develops a general perturbation methodology for constructing first-order and
second-order approximations to the solutions of MSDSGE models in which certain parameters
vary according to discrete Markov chain processes. Our method can be easily expanded to
higher-order approximations. The key is to find the approximations using the equilibrium conditions implied by the original economic model when Markov-switching parameters are present.
Thus, our method overcomes the aforementioned shortcomings. By working with the original
MSDSGE model directly rather than taking a system of linear rational expectations equations
with constant parameters as a short-cut, we maintain congruity between the original economic
model with Markov-switching parameters and the resultant approximations to the model solution. Such congruity is necessary for researchers to derive both first-order and second-order
approximations.
Unlike the case with a standard model with constant parameters, one conceptual difficulty
while working with MSDSGE models is that certain Markov-switching parameters, such as the
mean growth rate of technology, complicate the steady-state definition. The definition of the
steady-state for MSDSGE models should be independent of the realization of the discrete Markov
chain process for any changing parameter. To this end, we introduce a new concept of steadystate. We partition the parameter space such that the subset of Markov-switching parameters
that would influence the steady-state in the constant parameter case is now a function of the
perturbation parameter, while the complementary subset is not. This new concept renders a
key to deriving first-order and second-order approximations.
Several important results are as follows. First, we find that first-order approximations to
the solutions of MSDSGE models are, in general, not certainty equivalent when a subset of
Markov-switching parameters needs to be perturbed. Second, we identify the task of finding all
the solutions to a system of quadratic equations as the only bottleneck in obtaining first-order
2
and second-order approximations to the solutions of MSDSGE models.
Third, we propose to remove the bottleneck by applying Gröbner bases, a key insight of
our approach. Gröbner bases have not only a sound mathematical theory but also many computational applications. For example, Gröbner bases have been successfully applied by Kubler
and Schmedders (2010a) and Kubler and Schmedders (2010b) to find multiple equilibria in general equilibrium models; Datta (2010) uses such bases to find all the Nash equilibria. For our
purpose, Gröbner bases provide a computationally feasible way to obtain all the solutions to a
system of polynomial equations. Once the bottleneck of solving a system of quadratic equations
is resolved, the remaining task to obtain first-order and second-order approximations involves
solving only systems of linear equations even for higher-order approximations.
Fourth, one may use a numerical algorithm to search for a solution to the quadratic system,
but there is no guarantee that such an algorithm is capable of finding all solutions. By employing
Gröbner bases to solve the quadratic system, we can first obtain all its solutions and then
determine how many approximations are stable, where the stability concept follows the earlier
work of Costa et al. (2005), Farmer et al. (2009), Farmer et al. (2011), and Cho (2011). This
procedure enables researchers to ascertain both the existence and the uniqueness of a stable
approximation.
The rest of the paper is organized as follows. Section 2 presents a general class of MSDSGE
models, outlines our methodology, and introduces a new concept of steady-state. Section 3
derives first-order approximations and discusses the feature of no certainty equivalence. Section
4 shows how to reduce the problem of finding first-order approximations to that of finding all
the solutions to a system of quadratic equations. Gröbner bases are proposed to tackle this
problem and the concept of stability is introduced in this section. Section 5 derives a secondorder approximation and shows that it involves solving only linear systems once a first-order
approximation is obtained. Section 6 illustrates how to apply our methodology to three different
MSDSGE models. Concluding remarks are offered in Section 7.
3
2
The General Framework
Our general framework is laid out as follows. In Section 2.1 we discuss a general class of
MSDSGE models with a simple example to guide the reader through our new notation. In
Section 2.2 we introduce a new concept of steady-state and discuss a new idea of partitioning
the Markov-switching parameter space. In Section 2.3 we present the general form of model
solutions and state the goal of this paper.
2.1
The Model
We study a general class of MSDSGE models in which some of the parameters follow a discrete
Markov chain process indexed by st with the transition matrix P = (ps,s′ ). The element ps,s′
represents the probability that st+1 = s′ given st = s for s, s′ ∈ {1, . . . , ns }, where ns is the
number of regimes. When st = s, the model is said to be in regime s at time t. We denote the
vector of all Markov-switching parameters by θt , which has dimension nθ × 1.1
Given (xt−1 , εt , θt ), the equilibrium conditions for MSDSGE models have the general form
Et f (yt+1 , yt , xt , xt−1 , χεt+1 , εt , θt+1 , θt ) = 0nx +ny ,
(1)
where Et denotes the mathematical expectation operator conditional on information available
at time t, 0nx +ny is an (nx + ny ) × 1 vector of zeros, xt is an nx × 1 vector of (endogenous
and exogenous) predetermined variables, yt is an ny × 1 vector of non-predetermined (control)
variables, εt is an nε × 1 vector of i.i.d. innovations to the exogenous predetermined variables
with Et−1 εt = 0nε , and χ ∈ R is the perturbation parameter. Since the total number of equations
in (1) is ny + nx , the function f maps R2(ny +nx +nε +nθ ) into Rny +nx .
Finding first-order and second-order approximations to the solutions of MSDSGE models
using perturbation methods requires us to introduce new notation and lengthy algebraic work.
To aid the reader through our new notation and derivations in the paper, we use a simple real
business cycle (RBC) model as a clear demonstration of how our proposed methodology works.
1
Note that the vector θt does not include other parameters that are constant over time. We call those
parameters “constant parameters.”
4
The RBC model is an economy with the representative household whose preferences over a
stochastic sequence of consumption goods, ct , are represented by the utility function
max E0
∞
X
β t log ct
t=0
where β denotes the discount factor. The resource constraint is
α
ct + kt = zt kt−1
+ (1 − δ) kt−1 ,
where kt is a stock of physical capital and zt represents a technological change that is a random
walk in log with a Markov-switching drift as
log zt = µt + log zt−1 + σεt ,
where the drift µt takes two discrete values dictated by the Markov chain process represented
by st ∈ {1, 2}, and εt ∼ N (0, 1).
The three equations characterizing the equilibrium are
1
1
= βEt
αzt+1 ktα−1 + 1 − δ ,
ct
ct+1
α
ct + kt = zt kt−1
+ (1 − δ) kt−1 ,
and log zt = µt + log zt−1 + σεt .
Because log zt has a unit root, the economy is non-stationary. To derive a stationary equi1
1−α
librium, define ω t = zt−1
and let c̃t =
ct
,
ωt
k̃t−1 =
kt−1
,
ωt
z̃t =
zt
.
zt−1
The transformed (re-scaled)
equilibrium conditions become
1
z̃ α−1
1
αz̃t+1 k̃tα−1 + 1 − δ ,
= βEt t
c̃t
c̃t+1
1
α
+ (1 − δ) k̃t−1 ,
c̃t + k̃t z̃t1−a = z̃t k̃t−1
and log z̃t = µt + σεt .
Substituting out z̃t leads to the following two equilibrium conditions
1
1
µt + σεt
= βEt
exp
α exp µt+1 + σεt+1 k̃tα−1 + 1 − δ ,
c̃t
c̃t+1
α−1
5
µt + σεt
α
= exp (µt + σεt ) k̃t−1
+ (1 − δ) k̃t−1 .
and c̃t + k̃t exp
1−α
Using the notation in this section, we have ny = 1, nx = 1, nε = 1, nθ = 1, yt = c̃t ,
xt−1 = k̃t−1 , and θt = µt . The equilibrium condition (1) can be specifically expressed as
Et
2.2
1
c̃t
Et f (yt+1 , yt , xt , xt−1 , χεt+1 , εt , θt+1 , θt ) =
+σεt
1
α exp µt+1 + χσεt+1 k̃tα−1 + 1 − δ
exp µtα−1
− β c̃t+1
.
µt +σεt
α
c̃t + k̃t exp 1−α − exp (µt + σ t εt ) k̃t−1 − (1 − δ) k̃t−1
(2)
Steady-State and Partition of Markov-Switching Parameters
The definition of the steady-state for MSDSGE models is more complicated than that for standard DSGE models with constant parameters. To be consistent with the definition of the
steady-state in standard DSGE models with constant parameters, the definition of the steadystate for MSDSGE models should be independent of the realization of the discrete Markov chain
process.
To achieve this objective, our key idea is to partition the vector θt of Markov-switching
parameters into two subvectors. The first part concerns a subvector of parameters that would
influence the steady-state in the constant parameter case and thus is a function of the perturbation parameter χ. The second part contains the subvector of all remaining parameters that
would not affect the steady-state in the constant parameter case. We denote the first subvector
by θ 1t and the second subvector by θ2t . Rewrite the vector θt as
θt = θ (χ, st ) ,
where θ maps R × {1, . . . , ns } into Rnθ . We have
i⊺
h
i⊺ h
⊺
⊺
⊺
⊺
θt = θ1t θ2t
,
= θ 1 (χ, st ) θ2 (χ, st )
(3)
(4)
where ⊺ indicates transpose. The vector θ1t and θ2t have the dimensions nθ1 × 1 and nθ2 × 1
respectively, and functional forms
θ1 (χ, st ) = θ1 + χb
θ1 (st ) ,
and θ2 (χ, st ) = b
θ2 (st )
6
(5)
(6)
for all st .
The specific functional form (5) is chosen for tractable derivations in the rest of the paper.2
We will discuss, below, how to choose the value of θ 1 . We apply the same partition to the vector
θt+1 , where we simply replace the subscript t with t + 1.
Two important features stand out from (5) and (6). First, b
θ1 (st ) is a deviation of θ1t from
θ1 in regime st . Second, θ2t is not a function of the perturbation parameter χ. Thus, the
perturbation parameter, χ, affects only a subset of Markov-switching parameters, θ1t . Since
the steady-state depends on θ̄1 only, a natural choice for this point is the mean of the ergodic
distribution across θ1t .
Now we are ready to define the steady-state of the MSDSGE model.
Definition 1 The nx × 1 vector xss and the ny × 1 vector yss are the steady-state variables if
b
b
f yss , yss , xss , xss , 0nε , 0nε , θ1 , θ2 (st+1 ) , θ1 , θ2 (st ) = 0nx +ny
holds for all st+1 and st .
In principle, more than one partition of Markov-switching parameters can satisfy the condition, as stated in Definition 1, such that neither θ 2 (0, st+1 ) = b
θ2 (st+1 ) nor θ2 (0, st ) = b
θ2 (st )
enters in the calculation of the steady-state for any st+1 or st . We recommend choosing θ1t as
the smallest set of Markov-switching parameters such that θ̄1 influences the steady-state.
In our simple RBC model, there is no θ 2t , but θ1t = µt takes the functional form
µ (st ) .
µt = µ (χ, st ) = µ + χb
We set both εt = 0 and χ = 0. Thus, µt = µt+1 = µ and the steady-state of the RBC model
consists of c̃ss and k̃ss such that
µ̄
1
1
α−1
− β css exp α−1 α exp (µ̄) k̃ss + 1 − δ
css
= 02 ,
µ̄
α
c̃ss + k̃ss exp 1−α − exp (µ̄) k̃ss − (1 − δ) k̃ss
2
Any other functional form, so long as θ1 (0, st ) = θ1 holds for all st , will be valid.
7
which produces the steady-state values
1
−1+δ
,
µ̄
β exp α−1
µ̄
α
.
= exp (µ̄) k̃ss + (1 − δ) k̃ss − k̃ss exp
1−α
k̃ss =
and c̃ss
2.3
1
!! α−1
1
α exp (µ̄)
Model Solution and Approximation
Once the concept of steady-state and the partition of Markov-switching parameters are clearly
defined, we can now proceed to approximate the solution of a MSDSGE model. Consider model
solutions of the form3
yt = g (xt−1 , εt , χ, st ) ,
(7)
yt+1 = g (xt , χεt+1 , χ, st+1 ) ,
(8)
and xt = h (xt−1 , εt , χ, st )
(9)
where g maps Rnx +nε +1 × {1, . . . , ns } into Rny and h maps Rnx +nε +1 × {1, . . . , ns } into Rnx . In
the RBC model, we have
c̃t = g k̃t−1 , εt , χ, st ,
c̃t+1 = g k̃t , χεt+1 , χ, st+1 ,
and k̃t = h k̃t−1 , εt , χ, st .
In general, we do not know the explicit functional forms of g and h; thus, we need to
approximate them by Taylor expansions around the steady-state.4 In this paper, we present
the results up to second-order Taylor expansions, but our procedure can be easily expanded
to higher orders. For the rest of the paper, first-order Taylor expansions are referred to as
first-order approximations and second-order Taylor expansions as second-order approximations.
3
In theory, one could allow g and h to depend on the entire history of the regimes. In practice, this is not
tractable. The assumption that g and h depend only on the current regime corresponds to the notion of a
minimal state variable (MSV) solution in Farmer et al. (2011).
4
We assume that the solutions g and h are unique.
8
Given the steady-state defined in Definition 1, we have
yss = g (xss , 0nε , 0, st ) and xss = h (xss , 0nε , 0, st )
for all st and
yss = g (xss , 0nε , 0, st+1 ) and xss = h (xss , 0nε , 0, st+1 )
for all st+1 .
Using Equations (1), (7), (8), and (9), we write the function f as
g (h (xt−1 , εt , χ, st ) , χεt+1 , χ, st+1 ) , g (xt−1 , εt , χ, st ) ,
F (xt−1 , χεt+1 , εt , χ, st+1 , st ) = f
h (xt−1 , εt , χ, st ) , xt−1 , χεt+1 , εt , θ (χ, st+1 ) , θ (χ, st )
for all xt−1 , εt+1 , εt , st+1 , and st . The function F maps Rnx +2nε +1 × {1, . . . , ns } × {1, . . . , ns } into
Rny +nx . With the assumption that innovations to the exogenous predetermined variables, εt ,
are independent of the discrete Markov chain process, st , we write (1) as
Et f (yt+1 , yt , xt , xt−1 , χεt+1 , εt , θt+1 , θt ) =
Z
ns
X
′
pst ,s
F (xt−1 , χε′ , εt , χ, s′, st ) µ (ε′ ) dε′ = 0ny +nx
G (xt−1 , εt , χ, st ) =
(10)
s′ =1
for all xt−1 , εt , and st , where µ is the density function of the innovations. The function G maps
Rnx +nε +1 × {1, . . . , ns } into Rny +nx .
This paper has two goals. First, we show how to derive first-order and second-order approximations to the functions g and h around the steady-state. Second, we show that obtaining these
approximations requires solving a quadratic system. We propose a new methodology to find all
the solutions to such a system. Each of the solutions to the system corresponds to a different set
of first-order and second-order approximations. Finding all the solutions is a difficult task but
is crucial to determine how many approximations are stable.5 Sections 3-5, below, are devoted
to accomplishing these two objectives.
5
The stability definition is discussed in Section 4.2.3.
9
3
First-Order Approximations
This section gives a detailed description of how to derive first-order approximations to the model
solution represented by (7), (8), and (9). We proceed in several steps. We first lay out the
notation in Section 3.1 and then provide detailed derivations needed to find the approximations
in Section 3.2. With the notation and derivations in hand, we show in Section 3.3 how to obtain
first-order approximations to the model solution. In a final subsection, Section 3.4, we discuss
one of the most important features of MSDSGE models: the result of no certainty equivalence
of first-order approximations.
3.1
3.1.1
Notation
Partial Derivatives of f
We begin with the notation for the derivatives of
f (yt+1 , yt , xt , xt−1 , χεt+1 , εt , θt+1 , θt )
evaluated at the steady-state. Let
Dfss (st+1 , st ) =
h
i
Dj f i yss , yss , xss , xss , 0nε , 0nε , θ 1 , b
θ 2 (st+1 ) , θ1 , b
θ2 (st )
1≤i≤ny +nx ,1≤j≤2(ny +nx +nε +nθ )
denote the (ny + nx )×(2 (ny + nx + nε + nθ )) matrix of first partial derivatives of f with respect
to all its variables evaluated at
yss , yss, xss , xss , 0nε , 0nε , θ1 , b
θ2 (st+1 ) , θ1 , b
θ2 (st )
for all st+1 and st . To simplify notation, define
Dn,m fss (st+1 , st ) = [Dn fss (st+1 , st ) · · · Dm fss (st+1 , st )]
for all st+1 and st and 1 ≤ n < m ≤ 2 (ny + nx + nε + nθ ).6
6
When n = m, Dn,n fss (st+1 , st ) = Dn fss (st+1 , st ) for all st+1 and st .
10
Appendix A shows that our derivations depend on the following matrices:
D1,ny fss (st+1 , st ) , Dny +1,2ny fss (st+1 , st ) ,
D2ny +1,2ny +nx fss (st+1 , st ) , D2ny +nx +1,2(ny +nx ) fss (st+1 , st ) ,
D2(ny +nx )+1,2(ny +nx )+nε fss (st+1 , st ) , D2(ny +nx )+nε +1,2(ny +nx +nε ) fss (st+1 , st ) ,
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (st+1 , st ) , and D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (st+1 , st )
for all st+1 and st .
Given the importance of these derivatives, we return to the RBC example for concrete
illustration. It follows from (2) that
D1,1 fss (st+1 , st ) =
1
c2ss
0
, D2,2 fss (st+1 , st ) =
− c12
ss
1
,
kα−2
ss
0
(1 − α) αβ exp
css
, D4,4 fss (st+1 , st ) = exp( µ̄ ) ,
D3,3 fss (st+1 , st ) =
µ̄
1−α
exp 1−α
−
β
α−1
σ
αµ̄
kss
−αβ exp α−1
σ
(1−α)css
css
, D6,6 fss (st+1 , st ) =
µ̄
D5,5 fss (st+1 , st ) =
exp( 1−α
)kss
α
− exp (µ̄) kss σ
0
1−α
kα−1
µ̄α
1
ss
−αβ exp α−1 css
css
, and D8,8 fss (st+1 , st ) =
D7,7 fss (st+1 , st ) =
kss
µ̄
α
0
exp 1−α 1−α − exp (µ̄) kss
αµ̄
α−1
for all st+1 and st .
3.1.2
Partial Derivatives of G
We now introduce the notation for the derivatives of G (xt−1 , εt , χ, st ). Let
DG (xt−1 , εt , χ, st ) = Dj Gi (xt−1 , εt , χ, st ) 1≤i≤ny +nx ,1≤j≤nx +nε +1
refer to the (ny + nx ) × (nx + nε + 1) matrix of first partial derivatives of G with respect to
(xt−1 , εt , χ) for all xt−1 , εt , χ, and st . Note that there are no derivatives with respect to st since
it is a discrete variable. Similarly, we let
DG (xss , 0nε , 0, st ) = Dj Gi (xss , 0nε , 0, st ) 1≤i≤ny +nx ,1≤j≤nx +nε +1
11
refer to the (ny + nx ) × (nx + nε + 1) matrix of first partial derivatives of G with respect to
(xt−1 , εt , χ) evaluated at (xss , 0nε , 0, st) for all st .
To simplify notation, we define
DGss (st ) = DG (xss , 0nε , 0, st ) and Dj Giss (st ) = Dj Gi (xss , 0nε , 0, st )
for all st and 1 ≤ i ≤ ny + nx and 1 ≤ j ≤ nx + nε + 1. Thus,
DGss (st ) = Dj Giss (st ) 1≤i≤ny +nx ,1≤j≤nx +nε +1
for all st . Let Dj Gss (st ) denote the j th column vector of DGss (st ) for all st and 1 ≤ j ≤
nx + nε + 1. It follows that the first partial derivatives of G with respect to xt−1 evaluated at
(xss , 0nε , 0, st) for all st can be expressed as
[D1 Gss (st ) · · · Dnx Gss (st )] ,
the first partial derivatives of G with respect to εt evaluated at (xss , 0nε , 0, st) for all st can be
expressed as
[Dnx +1 Gss (st ) · · · Dnx +nε Gss (st )] ,
and the first partial derivative of G with respect to χ evaluated at (xss , 0nε , 0, st ) for all st is
Dnx +nε +1 Gss (st ) .
In contrast to one single first derivative in the constant parameter case, there are ns first derivatives of G, one for each possible value of st .
To simply notation further, we define
Dn,mGss (st ) = [Dn Gss (st ) · · · Dm Gss (st )]
for all st and 1 ≤ n < m ≤ nx + nε + 1.7 Therefore, D1,nx Gss (st ) represents the first partial
derivatives of G with respect to xt−1 evaluated at (xss , 0nε , 0, st ) for all st , Dnx +1,nx +nε Gss (st )
represents the first partial derivatives of G with respect to εt evaluated at (xss , 0nε , 0, st ) for
all st , and Dnx +nε +1 Gss (st ) is the first partial derivative of G with respect to χ evaluated at
(xss , 0nε , 0, st) for all st .
7
When n = m, then Dn,n Gss (st ) = Dn Gss (st ) for all st .
12
3.1.3
Partial Derivatives of g and h
We are now ready to introduce the notation for the derivatives of g and h evaluated at the
steady-state. Denote
Dg (xss , 0nε , 0, st ) = Dj g i (xss , 0nε , 0, st ) 1≤i≤ny ,1≤j≤nx+nε +1
as the ny × (nx + nε + 1) matrix of first partial derivatives of g with respect to (xt−1 , εt , χ)
evaluated at (xss , 0nε , 0, st ) for all st and
Dh (xss , 0nε , 0, st) = Dj hi (xss , 0nε , 0, st ) 1≤i≤nx ,1≤j≤nx +nε +1
as the nx × (nx + nε + 1) matrix of first partial derivatives of h with respect to (xt−1 , εt , χ)
evaluated at (xss , 0nε , 0, st ) for all st .
Simplifying the notation further, we define
i
(st ) = Dj g i (xss , 0nε , 0, st )
Dgss (st ) = Dg (xss , 0nε , 0, st ) and Dj gss
for all st , 1 ≤ i ≤ ny , 1 ≤ j ≤ nx + nε + 1, and
Dhss (st ) = Dh (xss , 0nε , 0, st ) and Dj hiss (st ) = Dj hi (xss , 0nε , 0, st)
for all st , 1 ≤ i ≤ nx , and 1 ≤ j ≤ nx + nε + 1. Thus,
i
Dgss (st ) = Dj gss
(st ) 1≤i≤ny ,1≤j≤nx +nε +1
and Dhss (st ) = Dj hiss (st ) 1≤i≤nx ,1≤j≤nx +nε +1
for all st . Let Dj gss (st ) be the j th column vector of Dgss (st ) and Dj hss (st ) be the j th column
vector of Dhss (st ) for all st and 1 ≤ j ≤ nx + nε + 1. We then define
Dn,m gss (st ) = [Dn gss (st ) · · · Dm gss (st )] and Dn,m hss (st ) = [Dn hss (st ) · · · Dm hss (st )]
for all st and 1 ≤ n < m ≤ nx + nε + 1.8
8
When n = m, then Dn,n gss (st ) = Dn gss (st ) and Dn,n hss (st ) = Dn hss (st ) for all st .
13
As will be discussed in Section 3.3, the following matrices of first partial derivatives are our
ultimate computational objects to obtain first-order approximations to the model solution
{D1,nx gss (st ) , D1,nx hss (st )}nsts=1 , {Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1 ,
and {Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1 ,
where D1,nx gss (st ) represents the first partial derivatives of g with respect to xt−1 evaluated
at (xss , 0nε , 0, st ) for all st , D1,nx hss (st ) represents the first partial derivatives of h with respect to xt−1 evaluated at (xss , 0nε , 0, st ) for all st , Dnx +1,nx +nε gss (st ) represents the first partial
derivatives of g with respect to εt evaluated at (xss , 0nε , 0, st ) for all st , Dnx +1,nx +nε hss (st ) represents the first partial derivatives of h with respect to εt evaluated at (xss , 0nε , 0, st ) for all st ,
Dnx +nε +1 gss (st ) is the first partial derivative of g with respect to χ evaluated at (xss , 0nε , 0, st )
for all st , and Dnx +nε +1 hss (st ) is the first partial derivative of h with respect to χ evaluated at
(xss , 0nε , 0, st) for all st .
Note that
D1,nx gss (st ) ∈ Cny ×nx , D1,nx hss (st ) ∈ Cnx ×nx ,
Dnx +1,nx +nε gss (st ) ∈ Cny ×nε , Dnx +1,nx +nε hss (st ) ∈ Cnx ×nε ,
Dnx +nε +1 gss (st ) ∈ Cny ×1 , and Dnx +nε +1 hss (st ) ∈ Cnx ×1
for all st . The symbol Cm1 ×m2 denotes an m1 × m2 matrix over the complex space. These computational objects can be found by solving systems of equations obtained through the chain rule
using the first partial derivatives of G (xt−1 , εt , χ, st ). The next subsection shows the derivations
needed to obtain such systems of equations.
3.2
Derivations
We now show how to use the chain rule using the first partial derivatives of G (xt−1 , εt , χ, st ) to
obtain systems of equations that are needed to solve for the matrices
{D1,nx gss (st ) , D1,nx hss (st )}nsts=1 , {Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1 ,
and {Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1 .
14
Denote an m1 × m2 matrix of zeros by 0m1 ×m2 . Since
G (xt−1 , εt , χ, st ) = 0ny +nx
for all xt−1 , εt , χ, and st , it follows that
DG (xt−1 , εt , χ, st ) = 0(ny +nx )×(nx +nε +1)
for all xt−1 , εt , χ, and st . In particular,
DGss (st ) = 0(ny +nx )×(nx +nε +1)
for all st so that
D1,nx Gss (st ) = 0(ny +nx )×nx ,
(11)
Dnx +1,nx +nε Gss (st ) = 0(ny +nx )×nε ,
and Dnx +nε +1 Gss (st ) = 0ny +nx
for all st . Note that in contrast to one single derivative in the constant parameter case, there
are ns first partial derivatives of G, one for each possible value of st .
The three expressions in (11) imply three systems of equations that will be used to find the
needed matrices
{D1,nx gss (st ) , D1,nx hss (st )}nsts=1 , {Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1 ,
and {Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1
respectively. In what follows, we describe how to use these three systems of equations in (11)
to obtain the needed matrices. We use the following steps to highlight the dependence between
these systems.
Step 1 The condition
D1,nx Gss (st ) = 0(ny +nx )×nx
implies
a
system
of
quadratic
{D1,nx gss (st ) , D1,nx hss (st )}nsts=1 .
15
equations
that
determines
Step 2 The conditions
Dnx +1,nx +nε Gss (st ) = 0(ny +nx )×nε
and Dnx +nε +1 Gss (st ) = 0ny +nx
imply two linear systems of equations that determine
{Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1
and
{Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1
as functions of {D1,nx gss (st ) , D1,nx hss (st )}nsts=1 .
Appendix A provides a detailed description of the derivations in Steps 1 and 2. In particular,
it derives the first partial derivatives of G needed to build the three systems. It shows that
{D1,nx gss (st ) , D1,nx hss (st )}nsts=1 can be determined by solving ns systems of quadratic equations
of the form
In x
D1,nx gss (1)
A (st )
..
.
D1,nx gss (ns )
In x
D1,nx hss (st ) = B (st )
D1,nx gss (st )
(12)
for all st , where Im denotes the identity matrix of size m×m, and A (st ) and B (st ) are functions
of the first partial derivatives of f evaluated at the steady-state as
A (st ) =
and
h P
ns
s′ =1
′
pst ,s′ D2ny +1,2ny +nx fss (s , st ) pst ,1 D1,ny fss (1, st ) · · · pst ,ns D1,ny fss (ns , st )
B (st ) = −
ns
X
s′ =1
for all st .
pst ,s′
h
′
′
D2ny +nx +1,2(ny +nx ) fss (s , st ) Dny +1,2ny fss (s , st )
i
i
The fact that we need to solve a quadratic system to determine {D1,nx gss (st ) , D1,nx hss (st )}nsts=1
is one of the key discoveries in the paper. The quadratic system has, in general, many solutions.
Each solution corresponds to a different first-order approximation. Finding all the solutions to
16
this quadratic system is a difficult task but is crucial to ascertain how many of them imply
stable approximations. In Section 4 we propose a new method to solve this quadratic system
for all its solutions.
Appendix A shows how, after finding a solution to (12), Step 2 implies that
{Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1 and {Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1 can be
obtained by simply solving the following two systems of linear equations
Dnx +nε +1 gss (1)
Dnx +1,nx +nε gss (1)
..
..
.
.
h
i−1
i
−1
Dnx +nε +1 gss (ns ) h
Dnx +1,nx +nε gss (ns )
= Θ Φ
Ψχ ,
Ψε and
ε
ε
= Θχ Φχ
Dnx +nε +1 hss (1)
Dnx +1,nx +nε hss (1)
..
..
.
.
Dnx +nε +1 hss (ns )
Dnx +1,nx +nε hss (ns )
(13)
where
ns
X
Θε =
′
s =1
′
p1,s′ Dny +1,2ny fss (s , 1) · · ·
..
..
.
.
0(nx +ny )×ny
..
.
· · · pns ,s′ Dny +1,2ny fss (s′ , ns )
0(nx +ny )×ny
,
0(nx +ny )×nx
p ′ D1,ny fss (s′ , 1) D1,nx gss (s′ ) · · ·
ns 1,s
X
..
..
..
Φε =
.
.
.
s′ =1
· · · pns ,s′ D1,ny fss (s′ , ns ) D1,nx gss (s′ )
0(nx +ny )×nx
′
′
f
(s
,
1)
·
·
·
0
D
p
2ny +1,2ny +nx ss
(nx +ny )×nx
ns 1,s
X
..
..
.
.
+
.
,
.
.
′
s =1
· · · pns ,s′ D2ny +1,2ny +nx fss (s′ , ns )
0(nx +ny )×nx
ns
X
Ψε = −
′
s =1
′
p1,s′ D2(ny +nx )+nε +1,2(ny +nx +nε ) fss (s , 1)
..
,
.
′
pns ,s′ D2(ny +nx )+nε +1,2(ny +nx +nε ) fss (s , ns )
17
Θχ
ns
X
=
s′ =1
+
Φχ
p
1,s′
′
Dny +1,2ny fss (s , 1) · · ·
..
..
.
.
0(nx +ny )×ny
..
.
· · · pns ,s′ Dny +1,2ny fss (s′ , ns )
p1,1 D1,ny fss (1, 1) · · · p1,ns D1,ny fss (ns , 1)
..
..
..
.
,
.
.
pns ,1 D1,ny fss (1, ns ) · · · pns ,ns D1,ny fss (ns , ns )
0(nx +ny )×ny
0(nx +ny )×nx
p ′ D1,ny fss (s′ , 1) D1,nx gss (s′ ) · · ·
ns 1,s
X
..
..
..
=
.
.
.
s′ =1
· · · pns ,s′ D1,ny fss (s′ , ns ) D1,nx gss (s′ )
0(nx +ny )×nx
′
′
0(nx +ny )×nx
p D2ny +1,2ny +nx fss (s , 1) · · ·
ns 1,s
X
..
..
.
.
+
.
,
.
.
s′ =1
· · · pns ,s′ D2ny +1,2ny +nx fss (s′ , ns )
0(nx +ny )×nx
′
′
f (s , 1) Dθss (s ) + . . .
D
p1,s′ 2(ny +nx +nε )+1,2(ny +nx +nε )+nθ ss
′
D
f
(s
,
1)
Dθ
(1)
2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) ss
ss
ns
X
..
and Ψχ = −
.
s′ =1
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s′ , ns ) Dθss (s′ ) + . . .
pns ,s′
′
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s , ns ) Dθss (ns )
,
where Dθss (st ) denotes the derivative of θ (χ, st ) with respect to χ evaluated at χ = 0 for all st .
That is,
for all st .
Dθss (st ) = Dθ (0, st ) = Dj θi (0, st ) 1≤i≤n
θ ,j=1
Since Φε and Φχ depend on {D1,nx gss (st ) , D1,nx hss (st )}nsts=1 , we must solve (12) prior to
solving (13). Since (13) is just a linear problem, we have that, first, there is a different solution to
(13) for each solution to (12) and, second, the bottleneck to finding first-order approximations is
obtaining {D1,nx gss (st ) , D1,nx hss (st )}nsts=1 by solving the quadratic system (12). The bottleneck
will be fully discussed in Section 4.
18
3.3
Approximations
Given a solution to (12) and (13), it is straightforward to build a first-order approximation to
the solutions of the MSDSGE model represented by Equation (1). Let a g first and hfirst be
first-order approximations to g and h around the point (xss , 0nε , 0, st ). It follows that
g first (xt−1 , εt , χ, st ) − yss = [D1 gss (st ) · · · Dnx gss (st )] (xt−1 − xss )
+ [Dnx +1 gss (st ) · · · Dnx +nε gss (st )] εt + Dnx +nε +1 gss (st ) χ
and
hfirst (xt−1 , εt , χ, st ) − xss = [D1 hss (st ) · · · Dnx hss (st )] (xt−1 − xss )
+ [Dnx +1 hss (st ) · · · Dnx +nε hss (st )] εt + Dnx +nε +1 hss (st ) χ
for all xt−1 , εt , and st . Hence, first-order approximations can be rewritten as
g first (xt−1 , εt , χ, st ) − yss = D1,nx gss (st ) (xt−1 − xss ) + Dnx +1,nx +nε gss (st ) εt + Dnx +nε +1 gss (st ) χ
and
hfirst (xt−1 , εt , χ, st ) − xss = D1,nx hss (st ) (xt−1 − xss ) + Dnx +1,nx +nε hss (st ) εt + Dnx +nε +1 hss (st ) χ
for all xt−1 , εt , and st . To express them in compact form, we have
g first (xt−1 , εt , χ, st ) − yss = Dgss (st ) St ,
where St =
h
(xt−1 − xss )
⊺
ε⊺t
χ
i⊺
and
hfirst (xt−1 , εt , χ, st ) − xss = Dhss (st ) St
for all xt−1 , εt , χ, and st .
We summarize what we have developed and what our next task is. For first-order approximations, Section 3.2 lays out a procedure for obtaining the matrices Dgss (st ) and Dhss (st )
for all st . The bottleneck is to obtain {D1,nx gss (st ) , D1,nx hss (st )}nsts=1 by solving the quadratic
system (12). Once the bottleneck is removed, the task of obtaining first-order approximations
19
only involves solving linear systems. As mentioned, the quadratic system has, in general, many
solutions. Each solution corresponds to a different first-order approximation. Finding all the
solutions to a quadratic system is central to ascertaining how many approximations are stable.
Section 4 is devoted to dealing with this bottleneck by finding all the solutions to the quadratic
system (12).
3.4
Feature of No Certainty Equivalence
As pointed out by Schmitt-Grohe and Uribe (2004), the certainty equivalence of first-order
approximations is a main result for a constant parameter model. This result implies that firstorder approximations to constant parameter models are inadequate for analyzing interesting
behavior such as economic agents’ responses to risk. For example, van Binsbergen et al. (2008)
and Rudebusch and Swanson (2008) argue that, when using constant parameter models, at least
second-order approximations are needed to analyze the effects of volatility on agents’ decisions.
While second-order approximations nullify the result of certainty equivalence, they result in a
substantially higher degree of computational difficulty in performing likelihood-based estimation,
as documented by Fernández-Villaverde and Rubio-Ramirez (2007). We show in this section,
however, that first-order approximations to the solutions of MSDSGE models are not necessarily
certainty equivalent. This salient feature opens the door to analyzing risk-related behaviors using
first-order approximations.
To see how certainty equivalence arises in first-order approximations in constant parameter
models, consider Equation (13) with only one regime so that ns = 1. In this degenerate case,
we have
h
where
h
Θχ Φχ
i
Dnx +nε +1 gss (1)
Dnx +nε +1 hss (1)
h
Θχ Φχ
i
= Ψχ ,
=
Dny +1,2ny fss (1, 1) + D1,ny fss (1, 1) D1,ny fss (1, 1) D1,nx gss (1) + D2ny +1,2ny +nx fss (1, 1)
20
(14)
i
and Ψχ = −
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (1, 1) Dθss (1) + . . .
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (1, 1) Dθss (1)
.
Because θ (χ, 1) = θ̄ in any constant parameter model, we have Dss θ (1) = 0nθ , implying that
Ψχ = 0nx +ny . Therefore the linear system (14) is homogeneous. If a unique solution exists, it is
given by
Dnx +nε +1 gss (1) = 0ny and Dnx +nε +1 hss (1) = 0nx .
(15)
For the constant parameter case, first-order approximations to policy rules are
g first (xt−1 , εt , χ, 1) − yss = D1,nx gss (1) (xt−1 − xss ) + Dnx +1,nx +nε gss (1) εt + Dnx +nε +1 gss (1) χ
and
hfirst (xt−1 , εt , χ, 1) − xss = D1,nx hss (1) (xt−1 − xss ) + Dnx +1,nx +nε hss (1) εt + Dnx +nε +1 hss (1) χ
for all xt−1 and εt . Using (15) and these policy rules evaluated at xt−1 = xss and εt = 0nε , we
have
g first (xss , 0nε , 1, 1) − yss = 0ny ,
hfirst (xss , 0nε , 1, 1) − xss = 0nx .
That is, first-order approximations to the solutions of the constant parameter model are certainty
equivalent.
With this insight we turn to the MSDSGE case. It is clear from Equation (13) that a
necessary condition for first-order approximations not to display certainty equivalence is Ψχ 6=
0ns (nx +ny ) . Let us analyze the circumstance under which the condition Ψχ 6= 0ns (nx +ny ) is true.
For visual convenience, we rewrite the expression for Ψχ
′
′
f (s , 1) Dθss (s ) + . . .
D
p1,s′ 2(ny +nx +nε )+1,2(ny +nx +nε )+nθ ss
′
D
f
(s
,
1)
Dθ
(1)
2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) ss
ss
ns
X
..
Ψχ = −
.
s′ =1
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s′ , ns ) Dθss (s′ ) + . . .
pns ,s′
′
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s , ns ) Dθss (ns )
21
.
Clearly, if Dθss (st ) = 0nθ for all st , then Ψχ = 0ns (nx +ny ) . Thus, a necessary condition for
Ψχ 6= 0ns (nx +ny ) to hold is Dθss (st ) 6= 0nθ for some st . Given the partition of θt described in (5)
and (6), we have
Dθss (st ) =
h
b
θ1 (st )⊺ 0⊺nθ
2
i⊺
.
It follows that Dθss (st ) 6= 0nθ for some st if and only if b
θ1 (st ) 6= 0nθ1 for some st . In sum, a
necessary condition for Ψχ 6= 0ns (nx +ny ) to hold is b
θ1 (st ) 6= 0nθ1 for some st .
The condition b
θ1 (st ) 6= 0nθ1 for some st , however, is insufficient for Ψχ 6= 0ns (nx +ny ) to be
true. A sufficient condition is
ns
X
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s′ , st ) Dθss (s′ ) +
=
pst ,s′
6 0nx +ny
′
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s , st ) Dθss (st )
s′ =1
for some st . This additional condition holds if θ1t does not enter the equilibrium conditions
multiplicatively with a variable to which the first partial derivative of f is zero when evaluated
at the steady-state. The following proposition summarizes our findings.
Proposition 2 First-order approximations to the solution of an MSDSGE model are not certainty equivalent if and only if both b
θ1 (st ) 6= 0nθ1 and
ns
X
s′ =1
for some st .
pst ,s′
′
′
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s , st ) Dθss (s ) +
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s′ , st ) Dθss (st )
6= 0nx +ny
Proof. The “if” part of the proof has been provided in the analysis prior to the stated proposition.
For the “only if” part, we prove it by contradiction. Note that the only way for
Dnx +nε +1 gss (st ) = 0ny and Dnx +nε +1 hss (st ) = 0nx
is for Ψχ = 0ns (nx +ny ) . Suppose b
θ1 (st ) = 0nθ1 for all st . Then Ψχ = 0ns (nx +ny ) and
Dnx +nε +1 gss (st ) = 0ny and Dnx +nε +1 hss (st ) = 0nx
22
for all st . It follows that first-order approximations are certainty equivalent.
Now suppose b
θ 1 (st ) 6= 0nθ1 for some st but
ns
X
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s′ , st ) Dθss (s′ ) +
= 0nx +ny
pst ,s′
′
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s , st ) Dθss (st )
s′ =1
for all st . In this case, it is straightforward to see that Ψχ = 0ns (nx +ny ) and
Dnx +nε +1 gss (st ) = 0ny and Dnx +nε +1 hss (st ) = 0nx
for all st . It follows that first-order approximations are certainty equivalent.
These contradictions establish the proof of the “only if” portion.
Proposition 2 implies that if first-order approximations are not certainty equivalent, approximated policy rules evaluated at xt−1 = xss and εt = 0nε are either
g first (xss , 0nε , 1, st ) − yss = Dnx +nε +1 gss (st ) 6= 0ny
or
hfirst (xss , 0nε , 1, st ) − xss = Dnx +nε +1 hss (st ) 6= 0nx
for some st .
4
Generalized Quadratic System and New Solution Method
One main discovery in this paper is that we identify solving a system of quadratic equations as
the sole bottleneck in obtaining first-order approximations to the solutions of MSDSGE models.
Once this bottleneck is removed, we need solve only two systems of linear equations to obtain
first-order approximations. In this section we show how to find all the solutions to this quadratic
system by introducing a new method. Finding all the solutions is essential to determining how
many first-order approximations are stable. Section 4.1 characterizes the quadratic system.
Section 4.2 proposes a new solution method by applying the theory of Gröbner bases. We show
that the number of solutions to the quadratic system is finite in most cases and define both
the concept and the condition of stability to be used in the paper. Section 4.3 uses our RBC
example to demonstrate how to apply our methodology.
23
4.1
Generalized Quadratic System
To understand the difficulty of solving a system of quadratic equations represented by (12), we
begin with the constant parameter case in which ns = 1. In this case, it is well understood that
one can map this special system of quadratic equations to the generalized Schur decomposition
problem related to the matrices A(1) and B(1). Thus, solving this special system of quadratic
equations is equivalent to using the standard matrix algebra procedure to find all the solutions.
In the Markov-switching case when ns > 1, there are ns systems of quadratic equations
represented by (12). Stacking these systems all together expands to a generalized system
s
. Because
of (ny + nx ) ns nx quadratic equations with unknowns {D1,nx gss (s) , D1,nx hss (s)}ns=1
s
appear in every system of the ns systems, the generalized quadratic system
{D1,nx gss (s)}ns=1
can no longer be mapped to the generalized Schur decomposition problem.9
The quadratic system represented by (12) has multiple solutions in general.10 Each solution
implies a different first-order approximation. It is crucial that we find all the solutions and then
determine how many imply stable approximations. For this purpose we develop, in the next
section, a new method by applying Gröbner bases to our problem.
4.2
4.2.1
New Solution Method
Gröbner Bases
The quadratic system discussed in Section 4.1 is simply a system of quadratic polynomials.
Gröbner bases provide a computationally practical means to obtain all the solutions to a system
of polynomial equations. A more detailed description of Gröbner bases is provided in Appendix
B. In this section we provide an intuitive explanation of how to apply Gröbner bases to solving
a system of multivariate polynomials.
Suppose one wishes to find all the solutions to a system of n polynomial equations in n
9
The intuition is that the generalized Schur decomposition cannot deal with the 2ns matrices A(1), . . . , A(ns )
and B(1), . . . , B(ns ) when ns > 1.
10
For the remainder of the paper, we refer to the generalized quadratic system (12) as the quadratic system.
24
unknowns
f1 (x1 , . . . , xn ) = 0, . . . , fn (x1 , . . . , xn ) = 0.
There exist a number of computationally efficient routines that can transform the original system
of n polynomial equations to an alternative system of n polynomial equations with the same
set of solutions. The following proposition, known as the Shape Lemma, describes a useful
transformation into a reduced Gröbner basis.
Proposition 3 Let
f1 (x1 , . . . , xn ) = 0, . . . , fn (x1 , . . . , xn ) = 0
be a system of n polynomials in n unknowns. Under certain regularity conditions, there exists
the following set of n polynomials in n unknowns with the same set of roots:
x1 − q1 (xn ) = 0, . . . , xn−1 − qn−1 (xn ) = 0, qn (xn ) = 0,
where each qi (xn ) is a univariate polynomial.
Proof. See Appendix B for details.
There are several important aspects to Proposition 3. First, the alternative set of polynomials
in Proposition 3 is known as a Shape basis, a special kind of reduced Gröbner basis. Second, the
regularity conditions, referred to in Proposition 3 and detailed in Appendix B, are satisfied by
the quadratic system in most economic problems. Third, the roots of the univariate polynomial
qn (xn ) can be found by any standard root finding algorithm. Fourth, once a root of qn (xn ) is
found, one can easily compute qi (xn ) to obtain x1 , . . . , xn−1 .
A large strand of the literature deals with the computation of reduced Gröbner bases. Buchberger (1998)’s algorithm is the original technique. Subsequently, many more efficient variants
have been proposed. We refer the interested reader to Cox et al. (1997). In this paper we use
Mathematica to find a Shape basis.
To illustrate how powerful Proposition 3 is, consider the following example featuring a system
25
of quadratic polynomials in four unknown variables x1 , . . . , x4 :
x1 x2 + x3 x4 + 2 = 0,
x1 x2 + x2 x3 + 3 = 0,
x1 x3 + x4 x1 + x4 x2 + 6 = 0,
x1 x3 + 2x1 x2 + 3 = 0.
A Shape basis is
x1 −
1
(9x54 + 6x34 − 15x4 ) = 0,
28
1
(−9x54 − 6x34 + 99x4 ) = 0,
28
1
x3 − (−3x54 − 9x34 − 2x4 ) = 0,
14
3x64 + 9x44 − 19x24 − 49 = 0.
x2 −
The last polynomial is univariate of degree six in x4 . The six roots of this polynomial are
{1.55461, −1.55461, 1.39592i, −1.39592i, 1.86232i, −1.86232i} .
Each of these roots can be substituted into the first three equations to obtain the following six
solutions:
{x1 = 2.89104, x2 = 1.7728, x3 = −4.58328, x4 = 1.55461},
{x1 = −2.89104, x2 = −1.7728, x3 = 4.58328, x4 = −1.55461},
{x1 = 0.372997i, x2 = 3.81477i, x3 = 0.41342i, x4 = 1.39592i},
{x1 = −0.372997i, x2 = −3.81477i, x3 = −0.41342i, x4 = −1.39592i},
{x1 = 4.81861i, x2 = 0.768342i, x3 = −0.914097i, x4 = 1.86232i},
{x1 == −4.81861i, x2 = −0.768342i, x3 = 0.914097i, x4 = −1.86232i}.
It is straightforward to show that these roots solve the original system of quadratic equations.
26
4.2.2
A Finite Number of Solutions
One of the regularity conditions for finding reduced Gröbner bases is that the quadratic system
has finitely many solutions. This condition is met in most economic problems. To see why this
result is true, let us first consider the constant parameter case when ns = 1. If the solution of
the model has a unique stable first-order approximation, the usual practice, as in Schmitt-Grohe
and Uribe (2004), involves constructing this stable approximation by ordering the generalized
eigenvalues of the two matrices A (1) and B (1) in a particular way.
In general, however,
the quadratic system has multiple solutions. Each solution implies a different first-order approximation. Some of these solutions may imply unstable approximations. For most economic
problems, the full set of approximations (stable and unstable) can be found by changing the
order of generalized eigenvalues. Hence, the number of approximations is related to the number
of possible orderings of eigenvalues. Therefore, the number of solutions to the quadratic system,
and therefore the number of approximations (stable and unstable), is bounded by
rank (A (1)) − nexo
(rank (A (1)) − nexo )!
=
,
(nx − nexo )! (rank (A (1)) − nx )!
nx − nexo
(16)
where nexo is the number of exogenous predetermined variables so that 0 ≤ nexo ≤ nx , rank (A (1))
stands for the rank of the matrix A (1).11 This result is familiar to most readers.
Now consider the Markov-switching case when ns > 1. We have the quadratic system or
ns quadratic systems of the form (12). If st were fixed for all t, the number of solutions to the
quadratic system, and therefore the number of approximations (stable and unstable), would be
bounded as in (16). Since we have ns regimes, the total number of solutions to the quadratic
11
Note that the matrix A (1) may not be of full rank when there are redundant variables that are linearly
dependent upon others and consequently can be eliminated from the quadratic system. A simple example is a
leisureless RBC model with three variables (capital, consumption, and output) and three equations (the Euler
condition, the resource constraint, and the output definition). If we eliminate the output definition equation and
the output variable, the newly formed matrix A(1) will be of full rank.
27
system for most parameter configurations is bounded by
ns
ns
rank (A (st )) − nexo
(rank
(A
(s
))
−
n
)!
t
exo
= max
max
st
st
(nx − nexo )! (rank (A (st )) − nx )!
nx − nexo
where rank (A (st )) stands for the rank of the matrix A (st ) for all st .12
Once we know that the number of solutions to the quadratic system is finite, we can use
Mathematica to obtain reduced Gröbner bases for finding all solutions to the quadratic system
and, hence, all first-order approximations. The next question is whether any of the approximations are stable and if so, how many. We address this question in the following section.
4.2.3
Mean Square Stability
In the constant parameter case, whether a first-order approximation is stable or not can be
determined by verifying whether its largest absolute generalized eigenvalue is greater than or
equal to one, a condition that holds for most concepts of stability. In the Markov-switching case,
the problem is subtle and complicated, and there are alternative concepts of stability. Given a
first-order approximation, we use the concept of mean square stability (MSS) as defined in Costa
et al. (2005). Farmer et al. (2009) discuss several advantages of using the MSS concept over
alternative ones such as the bounded stability. First, under MSS, even if the largest generalized
eigenvalues associated with one particular regime is greater than one, the system, as a whole, can
nonetheless be stable. Second, MSS permits applications in which the system has unbounded
errors and hence unbounded state variables. This unbounded feature holds in our RBC example
with the normally distributed shocks to TFP. Third, in the case of Markov-switching, necessary
and sufficient conditions for the MSS are easily verifiable, whereas other stability concepts do
12
For rare configurations of parameters, there may exist infinitely many solutions or the algorithm used
by Mathematica to find reduced Gröbner bases ceases to converge. In such a case, numerical procedures
can be used to solve the generalized quadratic system (12) to find one or possibly more solutions, although
none of those proecdures is guaranteed to find all solutions. Appendix C describes one such numerical procedure. Alternative numerical methods, such as Newton’s algorithm used in the RISE toolbox of Maih (2013)
(https://github.com/jmaih/RISE toolbox), may be used to find solutions, but again with no guarantee of finding
all solutions.
28
not have such conditions. Specifically, the MSS requires verifying whether the following matrix
has all its eigenvalues inside the unit circle
where
Υ=
T = P ⊺ ⊗ In2x Υ,
(17)
D1,nx hss (1) ⊗ D1,nx hss (1) · · ·
0n2x ×n2x
0n2x ×n2x
..
.
···
..
.
0n2x ×n2x
..
.
0n2x ×n2x
· · · D1,nx hss (ns ) ⊗ D1,nx hss (ns )
.
In the Markov-switching case, after we use a reduced Gröbner basis to obtain all the solutions
for {D1,nx gss (st ) , D1,nx hss (st )}nsts=1 , we verify, for each solution, whether the matrix Υ has all
its eigenvalues less than one. If only one solution satisfies the MSS criterion, the model has
a unique stable first-order approximation. If there is more than one solution that satisfies the
MSS criterion, the model has multiple stable first-order approximations. If none of the solutions
satisfies the MSS criterion or if there is no solution to the reduced Gröbner basis, the model
does not have any stable first-order approximations. At this point, it is worth emphasizing
that the stability of a second-order approximation (to be defined below) depends only on the
stability of its first-order approximation component. In other words, stability depends only on
the eigenvalues of the matrix Υ. The same is true for higher-order approximations.
4.3
The RBC model
At this junction we view it as instructive to use the previous RBC example to illustrate our
methodology developed thus far. Consider the following parameterization:
α
β
δ
σ
µ̄
µ
b (1)
µ
b (2)
p1,1
p2,2
0.3300 0.9976 0.0250 0.0002 0.00333 0.00167 −0.00163 0.90 0.90
The growth rates µ̄+b
µ (1) and µ̄+b
µ (2) correspond to regimes where annual output growth rates
are 3 percent and 1 percent respectively, β corresponds to a risk free annual rate of 3 percent
in the steady-state, and σ is set to match the total volatility of TFP growth as estimated in
29
Fernández-Villaverde and Rubio-Ramirez (2007). The standard calibration of α implies a capital
share of one third, and the value of δ implies an annual depreciation rate of approximately 10
percent, both in the steady-state. Note that regimes are symmetric in the sense that p1,1 = p2,2 .
Given this parameterization, the steady-state values of capital and consumption are kss =
32.0986 and css = 2.18946. We compute the first partial derivatives of f with respect to all the
variables evaluated at the steady-state:
0.20861
−0.2086
0.00031
, D2,2 fss (st+1 , st ) =
, D3,3 fss (st+1 , st ) =
D1,1 fss (st+1 , st ) =
0
1
1.00499
D4,4 fss (st+1 , st ) =
0
−1.0074
, D5,5 fss (st+1 , st ) =
D7,7 fss (st+1 , st ) =
−0.0147
0
for all st+1 and st .
0
0
, D6,6 fss (st+1 , st ) =
, and D8,8 fss (st+1 , st ) =
0.68169
44.9953
0.00014
0.00900
Using these results, one can build the quadratic system as in (12) and solve for
{D1,1 gss (st ) , D1,1 hss (st )}nsts=1 . Remember that in this example, nx = 1, ny = 1, and ns = 2.
The new method leads to the following four solutions
D1,1 hss (1)
D1,1 gss (1)
D1,1 hss (2)
D1,1 gss (2)
(1)
0.96364
0.03896
0.96364
0.03896
(2)
1.04023
−0.0380
1.04023
−0.0380
(3)
1.11326 + 0.116871i
−0.1114 − 0.11745i 1.11326 + 0.11687i
(4)
1.11326 − 0.116871i −0.1114 + 0.11745i
−0.1114 − 0.11745i
1.11326 − 0.11687i −0.1114 + 0.11745i
Mathematica finds the four solutions in less than a hundredth of a second. Using the results
in Section 4.2.3, one can verify that only solution (1) implies a stable first-order approximation
under the MSS criterion. Normally, when you have a unique stable first-order approximation,
you call it the first-order approximation. Thus, if we let ĉt = ct − css , k̂t−1 = kt−1 − kss , and
30
St =
h
k̂t−1 εt χ
i⊺
if st = 1, and
, the first-order approximation to the solution of the model is
ĉt
k̂t
ĉt
k̂t
=
0.96364 −0.0092 −0.0843
0.03896 0.00028 −0.00972
=
0.03896 0.00028
0.96364 −0.0092
0.00972
0.0843
St
St
if st = 2, where the rest of the derivatives used to form the first-order approximation can be
obtained by solving the linear system defined in (13).
The approximation highlights Proposition 2. First-order approximations are, in general, not
certainty equivalent. Since µ
b (st ) 6= 0 for all st and
ns
X
pst ,s′ (D7,7 fss (s′ , st ) Dθss (s′ ) + D8,8 fss (s′ , st ) Dθss (st )) 6= 02
s′ =1
for all st , first partial derivatives of g and h with respect to χ will be nonzero. That is,
D3 gss (st ) 6= 0 and D3 hss (st ) 6= 0.
5
Second-Order Approximations
Having shown how to construct first-order approximations to the solutions of the MSDSGE
model, we now show how to construct second-order approximations. In Section 5.1, we introduce
additional notation. In Section 5.2, we use the notation to derive second-order approximations to
the model solution. We show that given first-order approximations, second-order approximations
can be obtained by simply solving linear systems. This fact emphasizes that the bottleneck to
find both first-order and second-order approximations is solving the quadratic system defined
by (12). Section 5.3 returns to the RBC model for illustration.
31
5.1
5.1.1
Additional Notation
Second Partial Derivatives of f
We begin with additional notation in regard to second partial derivatives of f i . Denote
i
Hfss
(st+1 , st ) =
i
h
θ2 (st+1 ) , θ1 , b
θ2 (st )
Dk Dj f i yss, yss , xss , xss , 0nε , 0nε , θ1 , b
1≤j,k≤2(ny +nx +nε +nθ )
for all st+1 and st as the 2 (ny + nx + nε + nθ ) × 2 (ny + nx + nε + nθ ) matrix of second partial
derivatives of f i for 1 ≤ i ≤ ny + nx with respect to all its variables evaluated at
yss , yss , xss , xss , 0nε , 0nε , θ1 , b
θ2 (st+1 ) , θ1 , b
θ2 (st ) .
To conserve space, we do not represent the second partial derivatives of f for our RBC example
(they are available upon request).
5.1.2
Second Partial Derivatives of G
Let
HGi (xt−1 , εt , χ, st ) = Dk Dj Gi (xt−1 , εt , χ, st ) 1≤j,k≤nx +nε +1
be the (nx + nε + 1) × (nx + nε + 1) matrix of second partial derivatives of Gi with respect to
(xt−1 , εt , χ) for all xt−1 , εt , and st and 1 ≤ i ≤ ny + nx . It follows that
HGi (xss , 0nε , 0, st ) = Dk Dj Gi (xss , 0nε , 0, st ) 1≤j,k≤nx +nε +1
is the (nx + nε + 1) × (nx + nε + 1) matrix of second partial derivatives of Gi with respect to
(xt−1 , εt , χ) evaluated at (xss , 0nε , 0, st) for all st and 1 ≤ i ≤ ny + nx .
To simplify notation we define
HGiss (st ) = HGi (xss , 0nε , 0, st ) and Dk Dj Giss (st ) = Dk Dj Gi (xss , 0nε , 0, st )
for all st and 1 ≤ i ≤ ny + nx and 1 ≤ j, k ≤ nx + nε + 1. Thus,
HGiss (st ) = Dk Dj Giss (st ) 1≤j,k≤nx+nε +1
for all st and 1 ≤ i ≤ ny + nx .
32
5.1.3
Second Partial Derivatives of g and h
We now introduce the second partial derivatives of g and h. Let
Hg i (xt−1 , εt , χ, st ) = Dk Dj g i (xt−1 , εt , χ, st ) 1≤j,k≤nx+nε +1
be the (nx + nε + 1) × (nx + nε + 1) matrix of second partial derivatives of g i with respect to
(xt−1 , εt , χ) for all xt−1 , εt , and st and 1 ≤ i ≤ ny . It follows that
Hg i (xss , 0nε , 0, st ) = Dk Dj g i (xss , 0nε , 0, st ) 1≤j,k≤nx+nε +1
is the (nx + nε + 1) × (nx + nε + 1) matrix of second partial derivatives of g i with respect to
(xt−1 , εt , χ) evaluated at (xss , 0nε , 0, st) for all st and 1 ≤ i ≤ ny .
To put in compact notation we define
i
i
i
(st ) = Dk Dj g i (xss , 0nε , 0, st )
Hgss
(st ) = Hgss
(xss , 0nε , 0, st ) and Dk Dj gss
for all st and 1 ≤ i ≤ ny and 1 ≤ j, k ≤ nx + nε + 1. Hence,
i
Hgss
(st ) = Dk Dj g i (xss , 0nε , 0, st ) 1≤j,k≤nx+nε +1
for all st and 1 ≤ i ≤ ny .
Let
Hhi (xt−1 , εt , χ, st ) = Dk Dj hi (xt−1 , εt , χ, st ) 1≤j,k≤nx+nε +1
be the (nx + nε + 1) × (nx + nε + 1) matrix of second partial derivatives of hi with respect to
(xt−1 , εt , χ) for all xt−1 , εt , and st and 1 ≤ i ≤ nx . It follows that
Hhi (xss , 0nε , 0, st ) = Dk Dj hi (xss , 0nε , 0, st ) 1≤j,k≤nx+nε +1
is the (nx + nε + 1) × (nx + nε + 1) matrix of second partial derivatives of hi with respect to
(xt−1 , εt , χ) evaluated at (xss , 0nε , 0, st) for all st and 1 ≤ i ≤ nx .
If we define
Hhiss (st ) = Hhiss (xss , 0nε , 0, st) and Dk Dj hiss (st ) = Dk Dj hi (xss , 0nε , 0, st )
33
for all st , 1 ≤ i ≤ nx , and 1 ≤ j, k ≤ nx + nε + 1, then we have
Hhiss (st ) = Dk Dj hi (xss , 0nε , 0, st ) 1≤j,k≤nx+nε +1
for all st and 1 ≤ i ≤ nx .
With these definitions, the second partial derivatives of g and h evaluated at the steady-state
can be expressed as
Hgss (st ) =
and
Hhss (st ) =
for all st .
5.2
1
vec (Hgss
(st ))
..
.
n
vec Hgssy (st )
vec (Hh1ss
⊺
⊺
(st ))
⊺
..
.
vec (Hhnssx (st ))⊺
Approximations
Let g second and hsecond be a second-order approximation to g and h around the point (xss , 0nε , 0, st ).
With the additional notation introduced in Section 5.1, we have
1
g second (xt−1 , εt , χ, st ) − yss = Dgss (st ) St + Hgss (st ) (St ⊗ St ) ,
2
h
i⊺
where St = (xt−1 − xss )⊺ ε⊺t χ
is an (nx + nε + 1) × 1 vector. Similarly, we have
1
hsecond (xt−1 , εt , χ, st ) − xss = Dhss (st ) St + Hhss (st ) (St ⊗ St ) .
2
Similar to our analysis in Section 3.2, given the first-order approximations, the remaining
task is to derive the matrices
i
Hgss (st )
ny
i=1
, Hhiss (st )
nx ns
i=1 st =1
(18)
i
for all st , where Hgss
(st ) ∈ C(nx +nε +1)×(nx +nε +1) for all st and 1 ≤ i ≤ ny and Hhiss (st ) ∈
C(nx +nε +1)×(nx +nε +1) for all st and and 1 ≤ i ≤ nx .
34
As in the case of first-order approximations, we use the chain rule and the second partial
derivatives of G (xt−1 , εt , χ, st ) to obtain the systems of equations that can be used to solve for
the Hessian matrices expressed in (18). Since
G (xt−1 , εt , χ, st ) = 0ny +nx
for all xt−1 , εt , χ, and st , it must be the case that
HGi (xt−1 , εt , χ, st ) = 0(nx +nε +1)×(nx +nε +1)
for all xt−1 , εt , χ, st and 1 ≤ i ≤ ny + nx , and in particular
HGiss (st ) = 0(nx +nε +1)×(nx +nε +1)
for all st and 1 ≤ i ≤ ny + nx .13 It follows that
H1,nx ;1,nx Giss (st ) = 0nx ×nx , H1,nx ;nx +1,nx +nε Giss (st ) = 0nx ×nε ,
H1,nx ;nx +nε +1 Giss (st ) = 0nx , Hnx +1,nx +nε ;nx +1,nx +nε Giss (st ) = 0nε ×nε ,
(19)
Hnx +1,nx +nε ;nx +nε +1 Giss (st ) = 0nε , and Hnx +nε +1;nx +nε +1 Giss (st ) = 0
for all st and 1 ≤ i ≤ ny + nx , where we have used the following definition
D D Gi (s ) . . . Dn1 Dm2 Giss (st )
n1 m1 ss t
..
..
..
Hn1 ,n2 ;m1 ,m2 Giss (st ) =
.
.
.
Dn1 Dm1 Giss (st ) . . . Dn2 Dm2 Giss (st )
for all st , 1 ≤ n1 , m1 < n2 , m2 ≤ nx + nε + 1, and 1 ≤ i ≤ ny + nx .14
13
By Young’s Theorem, all the Hessian matrices in (18) and HGiss (st )
nx +ny
i=1
for all st are symmetric. In
the following, we exploit this fact whenever it is applicable and focus on only the relevant portion of the Hessian
matrices.
14
When m1 = m2 we have
Hn1 ,n2 ;m1 Giss (st ) =
35
Dn1 Dm1 Giss (st )
..
.
Dn1 Dm1 Giss (st )
All the equations in the systems represented in (19) depend on
{D1,nx gss (st ) , D1,nx hss (st )}nsts=1 , {Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1 ,
and {Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1 .
Thus, one must obtain first-order approximations prior to obtaining second-order approximations. The following steps describe how to use the systems in (19) to obtain the Hessian matrices
expressed in (18) and highlight the dependence between these systems.
Step 1 The condition H1,nx ;1,nx Giss (st ) = 0nx ×nx for 1 ≤ i ≤ ny + nx implies a linear system of
ny
nx ns
i
.
(st )}i=1
equations that determines the solution for {H1,nx ;1,nx gss
, {H1,nx ;1,nx hiss (st )}i=1
st =1
Step 2 Given the solution from Step 1, the condition H1,nx ;nx +1,nx +nε Giss (st ) = 0nx ×nε for 1 ≤ i ≤
ny + nx implies a linear system of equations that determines the solution for
ny
nx ns
i
.
(st )}i=1
{H1,nx ;nx +1,nx +nε gss
, {H1,nx ;nx +1,nx +nε hiss (st )}i=1
st =1
Step 3 Given the solution from Step 2, the condition H1,nx ;nx +nε +1 Giss (st ) = 0nx for 1 ≤ i ≤ ny +nx
implies a linear system of equations that determines
ny
nx ns
i
.
{H1,nx ;nx +nε +1 gss
(st )}i=1
, {H1,nx ;nx +nε +1 hiss (st )}i=1
st =1
the
solution
for
Step 4 Given the solution from Step 1, the condition Hnx +1,nx +nε ;nx +1,nx +nε Giss (st ) = 0nε ×nε for
1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution for
nx ns
ny
i
, {Hnx +1,nx +nε ;nx +1,nx +nε hiss (st )}i=1
(st )}i=1
{Hnx +1,nx +nε ;nx +1,nx +nε gss
.
st =1
Step 5 Given the solutions from Steps 1 and 3, the condition Hnx +1,nx +nε ,nx +nε ;nx +nε +1 Giss (st ) =
0nε for 1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution
ny
nx ns
i
for {Hnx +1,nx +nε ;nx +nε +1 gss
(st )}i=1
, {Hnx +1,nx +nε ;nx +nε +1 hiss (st )}i=1
.
st =1
for all st , 1 ≤ n1 , m1 < n2 ≤ nx + nε + 1, and 1 ≤ i ≤ ny + nx . When n1 = n2 we have
i
h
Hn1 ;m1 ,m2 Giss (st ) = Dn1 Dm1 Giss (st ) . . . Dn1 Dm2 Giss (st )
for all st , 1 ≤ n1 , m1 < m2 ≤ nx + nε + 1, and 1 ≤ i ≤ ny + nx . When m1 = m2 and n1 = n2 we have
Hn1 ;m1 Giss (st ) = Dn1 Dm1 Giss (st )
for all st , 1 ≤ n1 , m1 ≤ nx + nε + 1, and 1 ≤ i ≤ ny + nx .
36
Step 6 Given the solutions from Steps 1, 3 and 4, the condition Hnx +nε +1,nx +nε +1 Giss (st ) = 0 for
1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution for
ny
nx ns
i
.
{Hnx +nε +1,nx +nε +1 gss
(st )}i=1
, {Hnx +nε +1,nx +nε +1 hiss (st )}i=1
st =1
In the steps described above, we have used the following definitions:
i
(st )
D D g i (s ) . . . Dn1 Dm2 gss
n1 m1 ss t
.
.
.
i
..
..
..
Hn1 ,n2 ;m1 ,m2 gss
(st ) =
i
i
Dn1 Dm1 gss
(st ) . . . Dn2 Dm2 gss
(st )
for all st , 1 ≤ n1 , m1 < n2 , m2 ≤ nx + nε + 1, and 1 ≤ i ≤ ny ;
D D hi (s ) . . . Dn1 Dm2 hiss (st )
n1 m1 ss t
..
..
..
Hn1 ,n2 ;m1 ,m2 hiss (st ) =
.
.
.
Dn1 Dm1 hiss (st ) . . . Dn2 Dm2 hiss (st )
for all st , 1 ≤ n1 , m1 < n2 , m2 ≤ nx + nε + 1, and 1 ≤ i ≤ nx .15
15
When m1 = m2 , we have
i
(st ) =
Hn1 ,n2 ;m1 gss
for all st , 1 ≤ n1 , m1 < n2 ≤ nx + nε + 1, and 1 ≤ i ≤ ny ;
Hn1 ,n2 ;m1 hiss (st ) =
i
(st )
Dn1 Dm1 gss
..
.
i
(st )
Dn1 Dm1 gss
Dn1 Dm1 hiss (st )
..
.
Dn1 Dm1 hiss (st )
for all st , 1 ≤ n1 , m1 < n2 ≤ nx + nε + 1, and 1 ≤ i ≤ nx . When n1 = n2 , we have
i
(st ) =
Hn1 ;m1 ,m2 gss
h
i
i
(st )
(st ) . . . Dn1 Dm2 gss
Dn1 Dm1 gss
i
for all st , 1 ≤ n1 , m1 < m2 ≤ nx + nε + 1, and 1 ≤ i ≤ ny ;
Hn1 ;m1 ,m2 hiss (st ) =
h
Dn1 Dm1 hiss (st ) . . . Dn1 Dm2 hiss (st )
i
for all st , 1 ≤ n1 , m1 < m2 ≤ nx + nε + 1, and 1 ≤ i ≤ nx . When m1 = m2 and n1 = n2 , we have
i
i
(st )
(st ) = Dn1 Dm1 gss
Hn1 ;m1 gss
37
Appendix D provides a detailed description of the derivations for Steps 1 to 6. In particular,
it derives the second partial derivatives of G needed to construct the six systems. It shows
that once we obtain first-order approximations to the model solution, each of the equations in
Steps 1 to 6 becomes linear. As a result, for each first-order approximation, a second-order
approximation can be easily derived from the set of linear systems described above.
In addition, Appendix D shows that if first-order approximations are certainty equivalent,
the systems of equations represented by H1,nx ;nx +nε +1 Giss (st ) (cross partial derivatives between
xt−1 and χ) and Hnx +1,nx +nε ;nx +1,nx +nε Giss (st ) (cross partial derivatives between εt and χ) for
all st and 1 ≤ i ≤ ny + nx are homogeneous in the unknown variables
i
H1,nx ;nx +nε +1 gss
(st ) , H1,nx ;nx +nε +1 hjss (st ) ,
i
Hnx +1,nx +nε ;nx +nε +1 gss
(st ) , and Hnx +1,nx +nε ;nx +nε +1 hjss (st )
for all st , 1 ≤ i ≤ ny , and 1 ≤ j ≤ nx . This property is formally stated in the following
proposition.
Proposition 4 If first-order approximations are certainty equivalent, then the systems of equations
H1,nx ;nx +nε +1 Giss (st ) = 0nx and Hnx +1,nx +nε ;nx +nε +1 Giss (st ) = 0nx
for all st and 1 ≤ i ≤ ny + nx are homogeneous in their unknown variables and hence
i
H1,nx ;nx +nε +1 gss
(st ) = 0nx , H1,nx ;nx +nε +1 hjss (st ) = 0nx ,
i
Hnx +1,nx +nε ;nx +nε +1 gss
(st ) = 0nε and Hnx +1,nx +nε ;nx +nε +1 hjss (st ) = 0nε
for all st , 1 ≤ i ≤ ny , and 1 ≤ j ≤ nx .
Proof. If first-order approximations are certainty equivalent, one can see from Proposition 2
that
Dnx +nε +1 gss (st ) = 0ny and Dnx +nε +1 hss (st ) = 0nx
for all st , 1 ≤ n1 , m1 ≤ nx + nε + 1, and 1 ≤ i ≤ ny ;
Hn1 ;m1 hiss (st ) = Dn1 Dm1 hiss (st )
for all st , 1 ≤ n1 , m1 ≤ nx + nε + 1, and 1 ≤ i ≤ nx .
38
for all st . The proof follows directly from the above results.
In the constant parameter case, the cross partial derivatives with respect to χ and other state
variables are always zero because b
θ1 (1) = 0nθ1 . In the Markov-switching case, Proposition 4
i
implies that if b
θ1 (st ) 6= 0nθ1 for some st , then it may be the case that H1,nx ;nx +nε +1 gss
(st ) 6= 0nx ,
i
(st ) 6= 0nε or Hnx +1,nx +nε ;nx +nε +1 hjss (st ) 6= 0nε
H1,nx ;nx +nε +1 hjss (st ) 6= 0nx , Hnx +1,nx +nε ;nx +nε +1 gss
for some st , some 1 ≤ i ≤ ny , and some 1 ≤ j ≤ nx . These non-trivial extra terms may result
in more accurate second-order approximations.
5.3
Returning to the RBC Model
Let us now return to the RBC example. As shown above, given our parameterization, there is
only one stable first-order approximation. Thus, there is only one stable second-order approximation. When that is the case, we call it the second-order approximation. We can find the
second-order approximation by solving a set of linear systems described in Steps 1 to 6. We
have
1
+
2
ĉt
k̂t
=
0.03896
0.00028
0.00972
0.96364
−0.0092
−0.0843
St
−0.0004
4 × 10−6
0.00016
4 × 10−6
4 × 10−8
1 × 10−6
0.00016
1 × 10−6
−0.0003
−0.0002
−0.0003
−0.0025
−0.0003
2.7 × 10−6
2 × 10−5
−0.0025
0.00002
0.00057
(St ⊗ St )
if st = 1 or
1
+
2
ĉt
k̂t
=
0.03896
0.00028
−0.00972
0.96364
−0.0092
0.0843
St
−0.0004
4 × 10−6
−0.0002
4 × 10−6
4 × 10−8
−1 × 10−6
−0.0002
−1 × 10−6
−0.0003
−0.0002
−0.0002
0.00251
−0.0003
3 × 10−6
−2 × 10−5
0.00251
−2 × 10−5
0.00057
(St ⊗ St )
if st = 2.
Since the first-order approximation is not certainty equivalent in this example, the secondorder approximation highlights the result in Proposition 4. The cross derivatives of g and h
with respect to χ and other state variables evaluated at the steady-state are
Hi,3 gss (st ) = H3,i gss (st ) 6= 0, Hi,3 hss (st ) = H3,i hss (st ) 6= 0
39
for all st = 1, 2 and 1 ≤ i ≤ 2. Moreover, since regimes are symmetric, the following terms are
symmetric across regimes
D3 gss (1) = −D3 gss (2) , D3 hss (1) = D3 hss (2) ,
Hi,3 gss (1) = −Hi,3 gss (2) , Hi,3 gss (1) = −Hi,3 gss (2) ,
Hi,3 hss (1) = −Hi,3 hss (2) , and Hi,3 hss (1) = −Hi,3 hss (2)
for 1 ≤ i ≤ 2.
6
Applications
Since our theoretical results are new, the purpose of this section is to guide the reader by illustrating how to apply our methodology in practice. For this purpose we first continue the previous RBC model with new parameterization and then study two versions of the New-Keynesian
model.
6.1
Continuing the RBC Model
Stochastic neoclassical growth models, such as the one discussed in previous sections, are the
foundation of modern macroeconomics. Understanding how to solve MSDSGE models of this
prototype would enable us to work on richer models such as those commonly used for policy
analysis. For pedagogical reasons we consider the shock standard deviation, σ, to be constant
across regimes in all the examples studied in the paper. But our approach can be easily extended,
without much computational burden, to cases allowing for Markov-switching volatilities.
The previous RBC example restricts regimes to be symmetric such that p1,1 = p2,2 . As a
result, the cross partial derivatives of g and h evaluated at the steady-state are symmetric across
regimes. In this section, we consider an asymmetric-regime case in which p1,1 6= p2,2 . As will
be shown, the cross partial derivatives of g and h evaluated at the steady-state are asymmetric
across regimes in this case.
40
Consider the same parameter configuration as in Section 4.3, except p1,1 = 0.5, µ̄ = 0.00222,
µ (1) = 0.00278, and µ (2) = 0.00052. In this case, Regime 1 has a shorter expected duration, and
Regime 2 occurs more frequently in the ergodic distribution. With this alternative parameter
configuration, the steady-state is represented by kss = 34.6774, and css = 2.24769 and the first
partial derivatives of f evaluated at this steady-state are different from those in the symmetricregime case.16 Using these results we can solve the quadratic system (12).
There are four solutions to the quadratic system under this alternative parameter configuration
D1,1 hss (1)
D1,1 gss (1)
D1,1 hss (2)
D1,1 gss (2)
(1)
0.96545
0.0370821
0.96545
0.0370821
(2)
1.03828
−0.035996
1.03828
−0.035996
(3)
2.00373 − 0.7042i −1.00465 + 0.70654i
1.11318 + 0.39122i
−0.111145 − 0.39252i
(4)
2.00373 + 0.7042i
−1.00465 − 0.70654i 1.11318 − 0.39122i −0.111145 + 0.39252i
As before, Mathematica finds the four solutions in less than a hundredth of a second. Solution
(1) is the only one associated with a stable first-order approximation. Remember that the rest
of the coefficients necessary to construct the second-order approximation can be obtained by
solving linear systems and that, given that the first-order approximation is stable, the secondorder approximation is also stable. The second-order approximation to the model solution is
0.03708 0.00029
0.00637
ĉt
=
St
k̂t
0.96545 −0.0100 −0.1412
1 −0.0004 4 × 10−6 0.00009 4 × 10−6 5 × 10−8 1 × 10−6 0.00009 1 × 10−6 −5 × 10−6
(St ⊗ St )
+
0.00065
2 −0.0002 −0.0003 −0.0040 −0.0003 3 × 10−6 0.00004 −0.0040 0.00004
if st = 1 and
1
+
2
ĉt
k̂t
=
0.03708
0.00029
−0.0013
0.96545
−0.0100
0.02823
St
−0.0004
4 × 10−6
−1 × 10−5
4 × 10−6
5 × 10−8
−2 × 10−7
−1 × 10−5
−2 × 10−7
−6 × 10−5
−0.0002
−0.0003
0.00080
−0.0003
3 × 10−6
−8 × 10−6
0.00080
−8 × 10−6
0.00009
(St ⊗ St )
if st = 2.
16
In the interest of space, we omit many matrices from this section. They are available upon request.
41
Clearly, the first-order approximation can be constructed by only considering the first matrix
of the right hand side the above expressions.
As in the symmetric regime case, the first-order approximation is not certainty equivalent
and the cross derivatives of g and h with respect to χ and other state variables evaluated at
the steady-state are non-zero. These results follow directly from Propositions 2 and 4. Because
µ
b (st ) 6= 0 for all st , we have that
ns
X
pst ,s′ (D7,7 fss (s′ , st ) Dθss (s′ ) + D8,8 fss (s′ , st ) Dθss (st )) 6= 02
s′ =1
for all st . Consequently, D3 gss (st ) 6= 0, D3 hss (st ) 6= 0, and
Hi,3 gss (st ) = H3,i gss (st ) 6= 0, Hi,3 hss (st ) = H3,i hss (st ) 6= 0
for all st = 1, 2 and 1 ≤ i ≤ 2.
The results show that µ
b (st ) plays a crucial role in determining whether the first-order ap-
proximations are certainty equivalent and whether the cross partial derivatives of g and h with
respect to χ and other state variables evaluated at the steady-state are zero. Unlike the sym-
metric regime case, however, the cross partial derivatives of g and f at the steady-state are not
the same across regimes. Specifically,
D3 gss (1) 6= −D3 gss (2) , D3 hss (1) 6= D3 hss (2) ,
Hi,3 gss (1) 6= −Hi,3 gss (2) , H3,i gss (1) 6= −H3,i gss (2) ,
Hi,3 hss (1) 6= −Hi,3 hss (2) , and H3,i hss (1) 6= −H3,i hss (2)
for 1 ≤ i ≤ 3.
Given the first-order and second-order approximations, we can compute and compare their
accuracy using Euler equation errors as suggested in Judd (1998) and Aruoba et al. (2006). The
Euler equation error at point k̃t−1 , εt ; st for the approximation order ∈ {first, second} is
EE order k̃t−1 , εt , 1; st =
c̃order (k̃t−1 ,εt ,1;st )
µ(st )+σεt
Z X
2
exp
α−1
c̃order (k̃ order (k̃t−1 ,εt ,1;st ),ε′ ;1)
′
′
α−1
1−β
pst ,s′
µ (ε ) dε
× α exp (µ (s′ ) + σε′ ) k̃ order k̃t−1 , εt , 1; st
+1−δ
s′ =1
42
where c̃order k̃t−1 , εt , 1; st and k̃ order k̃t−1 , εt , 1; st are the appropriate approximations to the
policy functions. The unconditional absolute Euler equation error is
XZ
order
EE order k̃t−1 , εt , 1; st µorder k̃t−1 , εt , 1; st dk̃t−1 dεt
EE
=
st
where µorder is the unconditional distribution of the variables k̃ and ε.
To approximate this integral numerically, we simulate a long path for the economy, saving k̃t ,
εt , and st along the path. The simulation length is 10,000 periods, with the first 1,000 periods
discarded as a burn-in. Then, for each state k̃t−1 , εt , st , we draw 10,000 normally distributed
εt+1 ’s to compute the expectation over εt+1 , and use the transition values pst ,st+1 to compute
the expectation over st+1 . This entire process takes approximately an hour.
The following table shows the base-10 logarithms of absolute values of Euler equation errors
for the first-order and second-order approximations in both symmetric and asymmetric cases.
Both the first-order and second-order approximations produce a high degree of accuracy: a value
of −5 implies an error of $1 for each $100,000 of consumption. As can be seen, the second-order
approximation achieves an even higher degree of accuracy.
Unconditional Absolute Euler Equation Errors (log10 )
6.2
p1,1 = 0.9
p1,1 = 0.5
EE first
−4.6099
−5.3929
EE second
−5.4158
−6.1983
A New-Keynesian Model
New-Keynesian models, such as those in Woodford (2003) and Christiano et al. (2005), are often
used in policy analysis. This section discusses a version of the New-Keynesian model, illustrates
how to partition the vector of Markov-switching parameters θt , and discusses the conditions
under which multiple stable approximations exist.
43
6.2.1
The Model
We consider a New-Keynesian model with quadratic price adjustment costs. The monetary
authority follows a Taylor Rule. There are two Markov-switching parameters: the technology
drift parameter and the coefficient on inflation in the Taylor rule. Davig and Leeper (2007),
Farmer et al. (2011), and Bianchi (2010), among others, argue that Markov-switching in the
Taylor rule coefficient on inflation captures a switch in U.S. policy regime since the mid 1980s.
The model features the representative consumer maximizing the expected lifetime utility
over consumption Ct and hours worked Ht
E0
∞
X
β t (log Ct − Ht )
t=0
subject to the budget constraint
Ct +
Bt
Bt−1
= Wt Ht + Rt−1
+ Tt + Dt ,
Pt
Pt
where Bt is nominal bonds, Wt is the real wage, Rt−1 is the nominal return on bonds, Tt is
lump-sum transfers, and Dt is profits from firms. The first-order conditions are
1 = βEt
Ct Rt
Ct+1 Πt+1
and Ct = Wt .
The competitive final goods producer combines a continuum of intermediate goods Yj,t into
a final good Yt according to the constant elasticity of substitution (CES) aggregation technology
Yt =
Z
0
1
η−1
η
Yj,t dj
η
η−1
.
Intermediate goods are produced by firms taking the wage and the demand function
−η
Pj,t
Yt
Yj,t =
Pt
as given. The price Pj,t is set and hours Hj,t are demanded according to
Yj,t = At Hj,t ,
44
where At is a technology shock following the law of motion
log At = µt + log At−1 ,
where, similar to the RBC model, the drift µt takes two discrete values dictated by the Markov
chain process represented by st ∈ {1, 2}.
These intermediate-goods firms face quadratic price adjustment costs according to
2
κ
Pj,t
ACj,t =
− 1 Yt .
2 Pj,t−1
The firms maximize the expected discounted profits subject to the demand function, the production function, and adjustment costs. In the symmetric equilibrium, Pj,t = Pt , Yj,t = Yt , and
Hj,t = Ht for all j, and the optimality conditions are
Wt = mct At
and κ (Πt − 1) Πt = (1 − η) + ηmct + βκEt (Πt+1 − 1) Πt+1
Ct Yt+1
,
Ct+1 Yt
where mct is the marginal cost faced by the firms.
The monetary authority sets the interest rate Rt by the following Taylor rule
ρ
Rt−1
Rt
(1−ρ)ψt
exp (σεt )
Πt
=
Rss
Rss
where the coefficient on inflation, ψ t , takes two discrete values dictated by the same Markov
chain process represented by st .
The market clearing condition is
Y t = Ct +
κ
(Πt − 1)2 Yt .
2
Substituting out mct , Wt , and Ct leads to the equilibrium conditions
1 − κ2 (Πt − 1)2
Yt Rt
,
1 = βEt
2
κ
1 − 2 (Πt+1 − 1) Yt+1 Πt+1
1 − κ2 (Πt − 1)2
κ
2 Yt
,
κ (Πt − 1) Πt = (1 − η) + η 1 − (Πt − 1)
+ βEt κ (Πt+1 − 1) Πt+1
2
At
1 − κ2 (Πt+1 − 1)2
45
Rt
=
and
Rss
Rt−1
Rss
ρ
(1−ρ)ψt
Πt
exp (σεt ) .
Since the technology shock has a unit root, the model is non-stationary. To obtain a stationary equilibrium, we define Ỹt = Yt /At , which leads to the stationary equilibrium conditions
as
1 − κ2 (Πt − 1)2 Ỹt
R
1
t ,
1 = βEt
2
κ
1 − 2 (Πt+1 − 1) Ỹt+1 exp µt+1 Πt+1
2
κ
1
−
(Π
−
1)
κ
t
2
κ (Πt − 1) Πt = (1 − η) + η 1 − (Πt − 1)2 Ỹt + βEt κ
(Πt+1 − 1) Πt+1 ,
2
1 − κ2 (Πt+1 − 1)2
ρ
Rt
Rt−1
(1−ρ)ψt
and
=
exp (σεt ) .
Πt
Rss
Rss
We partition θt into θ1t = µt and θ2t = ψ t . We have
µ (st )
µt = µ (χ, st ) = µ + χb
b (st ) .
and ψ t = ψ (χ, st ) = ψ
Hence, the drift parameter µt depends on the perturbation parameter χ, while the Taylor-rule
coefficient on inflation ψ t does not. We choose this partition because µt would enter the definition
of steady-state in a constant parameter model, while ψ t would not.
h
i⊺
Using the notation in Section 2, yt = Πt , Ỹt and xt−1 = Rt−1 , we express the stationary
equilibrium condition as
f (yt+1 , yt , xt , xt−1 , χεt+1 , εt , θt+1 , θt ) =
(1− κ (Πt −1)2 )Ỹt
Rt
1 − β 1− κ (Π2 −1)2 Ỹ exp µ1
) t+1 ( t+1 ) Πt+1
( 2 t+1
(1− κ (Πt −1)2 )
(1 − η) + η 1 − κ2 (Πt − 1)2 Ỹt + βκ 1− κ 2(Π −1)2 (Πt+1 − 1) Πt+1 − κ (Πt − 1) Πt
)
( 2 t+1
Rt−1
Rss
ρ
(1−ρ)ψt
Πt
exp (σεt ) −
Rt
Rss
The model solution takes the form of
h
i⊺
= g (Rt−1 , εt , χ, st ) ,
Πt Ỹt
h
i⊺
= g (Rt , χεt+1 , χ, st+1 ) ,
Πt+1 Ỹt+1
and Rt = h (Rt−1 , εt , χ, st ) .
46
.
6.2.2
Solving the Model
The following subsections show the sequential steps of obtaining approximations to the solution
of the model. First, we find the steady-state. Second, we derive the first partial derivatives of
f with respect to all its variables evaluated at steady-state. Third, we construct second-order
approximations to the solution. We present the results for two parameter configurations to show
how to use reduced Gröbner bases to obtain approximations and the MSS criterion to ascertain
whether we have a unique stable approximation.
Steady-State To obtain the steady-state, we set χ = 0 and εt = 0, so that Πt = Πt+1 = Πss ,
Ỹt = Ỹt+1 = Ỹss , Rt = Rt−1 = Rss , and µt+1 = µt = µ̄. Thus the equilibrium conditions in the
steady-state become
(1− κ (Πss −1)2 )Ỹss 1 Rss
1 − β 1− κ2 (Π −1)2 Ỹ exp(µ̄)
Πss
) ss
( 2 ss
2
κ
(1− (Πss −1) )
(1 − η) + η 1 − κ2 (Πss − 1)2 Ỹss + βκ 1− κ2 (Π −1)2 (Πss − 1) Πss − κ (Πss − 1) Πss
)
ρ ( 2 ss
(1−ρ)ψ
Rss
Rss
t
Πss
−R
Rss
ss
Assuming Πss = 1, we have the steady-state values as
Rss =
= 03×1 .
exp (µ̄)
η−1
, Ỹss =
.
β
η
This result confirms our partition such that µ̄ affects the steady-state, while ψ t does not. In
principle, we could let θ1t = (µt , ψ t ) so that perturbation would apply to ψ t as well. To obtain
a numerical solution that is as accurate as possible, however, one should keep the number
of perturbed parameters at the minimum. The converse is not true. That is, one cannot
let θ2t = (µt , ψt ) because the steady-state in the constant parameter case would depend on
technology drift.
47
Partial Derivatives of f
In this example, ny = 2, nx = 1, nε = 1, and nθ = 2. Thus, the
first partial derivatives of f with respect to all its variables evaluated at the steady-state are
η
η
1
0
η−1
1−η
D1,2 fss (st+1 , st ) = 0 βκ , D3,4 fss (st+1 , st ) = η
−κ
b
0
0
0 (1 − ρ) ψ (s)
−µ̄
0
0
−βe
D5,5 fss (st+1 , st ) =
, D6,6 fss (st+1 , st ) = 0 , D7,7 fss (st+1 , st ) = 0 ,
0
−µ̄
−µ̄
0
ρβe
−βe
0 0
1 0
0
D8,8 fss (st+1 , st ) = 0 , D9,10 fss (st+1 , st ) = 0 0 , and D11,12 fss (st+1 , st ) = 0 0 ,
0 0
0 0
σ
b (st ) enters the matrices because we do not perturb ψ .
for all st+1 and st . Note that ψ
t
Unique Stable Approximation Consider the following parameterization
β
κ
η
ρ
σ
µ̄
µ̂ (1)
µ̂ (2)
0.9976 161 10 0.8 0.0025 0.005 0.0025 −0.0025
b (1) ψ
b (2)
ψ
3.1
0.9
p1,1
p2,2
0.90 0.90
The growth rates µ̄ + µ̂ (1) and µ̄ + µ̂ (2) correspond to regimes where the annual growth rate
are 3 percent and 1 percent respectively, β corresponds to an annual risk-free rate of 3 percent,
η has the steady-state markup of 11 percent, and ρ and σ match the estimates in FernandezVillaverde et al. (2009).
b (1) and ψ
b (2) are such that
The two monetary policy parameters ψ
b (1) would imply a unique stable approximation if ψ = ψ
b (1) for all t and ψ
b (2) would lead to
ψ
t
b (2) for all t.
multiple stable approximations if ψ t = ψ
Given this parameter configuration, we can calculate the steady-state values of the nominal
rate and output as Rss = 1.0074 and Ỹss = 0.90. The resulting first partial derivatives of f with
respect to all its variables evaluated at the steady-state are
−1.11111
0
1.11111
1
D1,2 fss (st+1 , st ) =
10
−161.
0
160.614 , D3,4 fss (st+1 , st ) =
b (st )
0
0.2ψ
0
0
48
,
D8,8 fss (st+1 , st ) =
for all st+1 and st .
0
, D7,7 fss (st+1 , st ) = 0 ,
0
0
0
0.77941
−0.9926
0 0
1 0
0
, D9,10 fss (st+1 , st ) = 0 0 , and D11,12 fss (st+1 , st ) = 0 0
0
0 0
0 0
0.0025
D5,5 fss (st+1 , st ) =
−0.9926
0
, D6,6 fss (st+1 , st ) =
In this example, we obtain the following nine solutions for {D1,1 gss (st ) , D1,1 hss (st )}nsts=1
D1,1 gss (1)⊺
D1,1 hss (1)
D1,1 gss (2)⊺
D1,1 hss (2)
(1)
0.59517
−1.92815
−0.327932
0.699414
−2.9541
−0.554689
(2)
0.77508
−3.64018
−0.0398952
1.3018
−7.43725
2.76721
(3)
0.79559
−1.82393
−0.00706061
1.05423
1.21892
1.40196
(4)
1.0939 − 0.4363i
−0.8264 + 4.2641i
0.4706 − 0.6986i
1.3311 + 0.0574i
−10.008 − 1.9739i
2.9287 + 0.3165i
(5)
1.0939 + 0.4363i
−0.8264 − 4.2641i
0.4706 + 0.6986i
1.3311 − 0.0574i
−10.008 + 1.9739i
2.9287 − 0.3165i
(6)
1.0952 − 0.2105i
−0.9833 + 1.9595i
0.4727 − 0.3370i
1.0240 − 0.0200i
0.8689 + 0.7833i
1.2351 − 0.1103i
(7)
1.0952 + 0.2105i
−0.9833 − 1.9595i
0.4727 + 0.3370i
1.0240 + 0.0200i
0.8689 − 0.7833i
1.2351 + 0.1103i
(8)
1.2360 − 0.2511i
0.7554 + 3.0821i
0.6980 − 0.4020i
0.7507 + 0.0047i
−2.2696 + 0.6345i
−0.2718 + 0.0260i
(9)
1.2360 + 0.2511i
0.7554 − 3.0821i
0.6980 + 0.4020i
0.7507 − 0.0047i
−2.2696 − 0.6345i
−0.2718 − 0.0260i
Mathematica finds the nine solutions in less than three hundredths of a second. The only
solution that produces a stable first-order approximation is (1) and hence there is a unique stable
approximation. Given solution (1), we can solve the linear systems that allow us to obtain the
second-order approximation. Let Yb̃ t = Ỹt − Yss , Π̂t = Πt − Πss , R̂t = Rt − Rss , and define
h
i⊺
St = R̂t−1 εt χ . The second-order approximation is
1
+
2
Yb̃ t
Π̂t
R̂t
=
−1.9282
−0.0062
0.00481
−0.3279
−0.0011
0.00014
0.59517
0.00191
0.00008
St
21.3771
0.06247
−0.0188
0.06247
0.00020
−0.0001
−0.0188
0.00020
−0.0004
0.49793
0.00056
−0.0008
0.00056
2 × 10−6
−2 × 10−6
−0.0008
2 × 10−6
−0.0003
−0.1986
0.00124
−0.0004
0.00124
4 × 10−6
−1 × 10−6
−0.0004
4 × 10−6
−0.0002
49
(St ⊗ St )
if st = 1 and
1
+
2
Yb̃ t
=
Π̂t
R̂t
−2.9541
−0.0090
−0.0094
−0.5547
−0.0017
−0.0026
0.69941
0.00214
−0.0006
St
56.9733
0.16487
0.23174
0.16487
0.00050
0.00071
0.23174
0.00050
−0.0016
0.99333
0.00136
0.00033
0.00136
4 × 10−6
1 × 10−6
0.00033
4 × 10−6
−0.0010
−0.1842
0.00160
−0.0005
0.00160
5 × 10−6
−2 × 10−6
−0.0005
5 × 10−6
−0.0002
if st = 2.
(St ⊗ St )
Again, the first-order approximation is easily obtained from the above expressions. As in the
RBC model, there is no certainty equivalence and the cross derivatives of g and h with respect
to χ and other state variables evaluated at the steady-state are non-zero. Note that the first
partial derivatives of g and h with respect to Rt−1 and εt are different across regimes. This result
occurs because ψ t affects the first derivatives of f . Perturbing ψ t would make these derivatives
equal and the approximation less accurate.
Accuracy: Euler Equation Errors The Euler equation error at the point (Rt−1 , εt ; st ) is
that for order ∈ {first, second},
EE
order
(Rt−1 , εt , 1; st ) = 1 − β
Z X
2
s′ =1
1− κ2 (Πorder (Rt−1 ,εt ,1;st )−1)
2
pst ,s′ ×
×
Ỹ order (Rt−1 ,εt ,1;st )
µ (ε′ ) dε′
×
order
order
(Rt−1 ,εt ,1;st ),εt+1 ,1;st+1 )
Ỹ
(R
Rorder (Rt−1 ,εt ,1;st )
1
exp(µ′ ) Πorder (Rorder (Rt−1 ,εt ,1;st ),εt+1 ,1;st+1 )
1− κ2 (Πorder (Rorder (Rt−1 ,εt ,1;st ),εt+1 ,1;st+1 )−1)
2
where Ỹ order (Rt−1 , εt , 1; st ), Πorder (Rt−1 , εt , 1; st ), and Rorder (Rt−1 , εt , 1; st) are the approximations to the policy functions. The unconditional absolute Euler equation error is
XZ
order
EE order (Rt−1 , εt , 1; st) µorder (Rt−1 , εt , 1; st ) dRt−1 dεt
EE
=
st
where µorder is the unconditional distribution of the variables R and ε.
Numerical approximation of this integral simulates a 10, 000-period path, with the first 1, 000
periods discarded as a burn-in. For each state (Rt−1 , εt , st ) along the simulated path, we draw
50
10,000 normally distributed εt+1 ’s to compute the expectation over εt+1 , and use the transition
values pst ,st+1 to compute the expectation over st+1 . This entire process takes about an hour.
The following table shows the base-10 logarithms of absolute Euler equation errors for the
first-order and second-order approximations. Both the first-order and the second-order approximations produce a high degree of accuracy: a value of −4 implies an error of $1 for each $10,000
of consumption.
Unconditional Absolute Euler Equation Errors (log10 )
EE first
−3.7395
EE second
−4.7485
Multiple Stable Approximations As an alternative parameter configuration, we consider
b (2) = 0.7. That is, the second regime
the same parameters as in the previous section except ψ
now has a slightly lower response by the monetary authority to inflation. In this case, the
steady-state is still the same. As before, Regime 1 would imply a unique stable approximation
if considered in isolation, whereas Regime 2 would imply multiple stable approximations.
The first partial derivatives of f with respect to all its variables evaluated at the steady-state
are the same as those in the previous section except
−1.11111
0
D3,4 fss (st+1 , st ) =
10
−161.
b (st )
0
0.2ψ
b (st ). Using these results, we obtain the following 9
for all st+1 and st , which depends on ψ
51
solutions for {D1,1 gss (st ) , D1,1 hss (st )}nsts=1 .
D1,1 gss (1)⊺
D1,1 hss (1)
D1,1 gss (2)⊺
D1,1 hss (2)
(1)
0.59067
−1.9452
−0.3351
0.71244
−3.2185
(2)
0.79733
−4.7813
−0.0043
1.32443
−11.313
3.71833
(3)
0.85231
−1.7727
0.08374
1.01525
2.03718
1.52618
−0.6209
(4)
1.0444 − 0.4989i
−0.7000 + 4.6327i
0.3912 − 0.7987i
1.3523 + 0.0442i
−14.210 − 1.8733i
3.9161 + 0.3132i
(5)
1.0444 + 0.4989i
−0.7000 − 4.6327i
0.3912 + 0.7987i
1.3523 − 0.0442i
−14.210 + 1.8733i
3.9161 − 0.3132i
(6)
1.0670 − 0.1894i
−1.2845 + 1.5032i
0.4274 − 0.3033i
0.9995 − 0.0089i
1.6635 + 0.5703i
1.4141 − 0.0629i
(7)
1.0670 + 0.1894i
−1.2845 − 1.5032i
0.4274 + 0.3033i
0.9995 + 0.0089i
1.6635 − 0.5703i
1.4141 + 0.0629i
(8)
1.2374 − 0.2527i
0.8018 + 3.0980i
0.7004 − 0.4046i
0.7582 + 0.0045i
−2.3764 + 0.6699i
−0.2963 + 0.0317i
(9)
1.2374 + 0.2527i
0.8018 − 3.0980i
0.7004 + 0.4046i
0.7582 − 0.0045i
−2.3764 − 0.6699i
−0.2963 − 0.0317i
Now we have two solutions – (1) and (3) – that produce stable approximations. Thus the solution
to the model has multiple stable approximations. This example illustrates how we can use
reduced Gröbner bases and the MSS criterion to determine the number of stable approximations.
6.3
A New-Keynesian Model with Habit
Now consider a slight variant of the previously discussed New-Keynesian model, but with the
representative household having habit formation. We consider this application because deriving
first-order and second-order approximations is more involved when habit formation is present
and the reader would be interested in knowing how to accomplish this task. The household in
this economy maximizes
E0
∞
X
β t (log (Ct − ϕCt−1 ) − Ht )
t=0
where ϕ denotes habit persistence.
Letting λt denote the Lagrange multiplier on the budget
constraint, the first-order conditions are
λt =
1
ϕ
− βEt
,
Ct − ϕCt−1
Ct+1 − ϕCt
Rt
,
λt = βEt λt+1
Πt+1
and 1 = λt Wt .
Final good and intermediate goods producers face problems identical to the previous example
and
log At = µt + log At−1 .
52
Thus, the firm’s optimality conditions are
Wt = mct At
and κ (Πt − 1) Πt = (1 − η) + ηmct + βκEt (Πt+1 − 1) Πt+1
λt+1 Yt+1
.
λt Y t
The resource constraint is
Y t = Ct +
κ
(Πt − 1)2 Yt .
2
For clarity we assume that the monetary authority does not smooth interest rates and follows
the rule
ψ
Rt = Rss Πt t exp (σεt ) .
With habits, consumption appears at different dates, as Ct−1 , Ct , and Ct+1 , in the equilibrium
conditions. Substituting out mct , Wt , Yt , and Rt , we have the equilibrium conditions as
λt =
ϕ
1
− βEt
,
Ct − ϕCt−1
Ct+1 − ϕCt
ψ
Rss Πt t
λt = βEt λt+1
,
Πt+1
2
λt+1 Ct+1 1 − κ2 (Πt − 1)
η
+ βκEt (Πt+1 − 1) Πt+1
.
and κ (Πt − 1) Πt = (1 − η) +
At λt
λt Ct 1 − κ2 (Πt+1 − 1)2
The economy is nonstationary because of the presence of the unit root. If we define λ̃t = λt At
and C̃t = Ct /At , we have the stationary equilibrium conditions as
λ̃t =
ϕ
1
,
− βEt
C̃t − ϕ exp (−µt ) C̃t−1
C̃t+1 exp µt+1 − ϕC̃t
ψ
λ̃t+1
R Π t exp (σεt )
ss t
λ̃t = βEt
,
Πt+1
exp µt+1
2
λ̃t+1 C̃t+1 1 − κ2 (Πt − 1)
η
+ βκEt (Πt+1 − 1) Πt+1
and κ (Πt − 1) Πt = (1 − η) +
2.
λ̃t
λ̃t C̃t 1 − κ2 (Πt+1 − 1)
To solve the model, we define the two auxiliary variables X̃t = C̃t , and X̃t+1 = C̃t+1 . Using
h
i⊺
the notation in Section 2, we have yt = Πt X̃t λ̃t , xt−1 = C̃t−1 , θ2t = ψ t , and the
53
equilibrium conditions are
f (yt+1 , yt , xt , xt−1 , χεt+1 , εt , θt+1 , θt ) =
(1 − η) +
6.3.1
1
C̃t −ϕ exp(−µt )C̃t−1
− β X̃
ϕ
t+1 exp(µt+1 )−ϕC̃t
ψ
Rss Πt t exp(σεt )
β expλ̃t+1
Πt+1
(µt+1 )
η
λ̃t
− λ̃t
1− κ (Πt −1)2
+ βκ (Πt+1 − 1) Πt+1 λ̃λ̃t+1 X̃C̃t+1 1− κ 2(Π
t
t
2
− λ̃t
t+1 −1)
2
X̃t − C̃t
− κ (Πt − 1) Πt
Solving the Model
Similar to the previous examples, the following subsections show how to solve the model: find the
steady-state, define the matrices of first partial derivatives of f with respect to all its variables
evaluated at the steady-state, and find the second-order approximation to the policy functions.
As in the basic New-Keynesian model, we consider two parameter configurations to show how
to use reduced Gröbner bases to find all the approximations and how to use the MSS criterion
to determine how many of them are stable.
Steady-State Calculating the steady-state involves setting χ = 0 and εt = 0. Therefore, we
have Πt = Πt+1 = Πss , X̃t = X̃t+1 = X̃ss , C̃t = C̃t−1 = C̃ss , and µt+1 = µt = µ̄. The equilibrium
conditions in the steady-state are
(1 − η) +
1
C̃ss −ϕ exp(−µ̄)C̃ss
η
λ̃ss
−
λ̃ss
β exp(µ̄)
ϕ
β X̃ exp(µ̄)−ϕ
C̃ss
ss
ψ
Rss Πsst
Πss
− λ̃ss
− λ̃ss
1− κ (Π −1)2
ss
ss
2
+ βκ (Πss − 1) Πss λ̃λ̃ss X̃
− κ (Πss − 1) Πss
C̃ 1− κ (Π −1)2
ss
ss
2
ss
X̃ss − C̃ss
Assuming Πss = 1, the steady-state satisfies
exp (µ̄)
,
β
η
,
λ̃ss =
η−1
exp (µ̄) − βϕ η − 1
and C̃ss = X̃ss =
.
exp (µ̄) − ϕ η
Rss =
54
= 04×1 .
Partial Derivatives of f
In this example, ny = 3, nx = 1, nε = 1, and nθ = 2. Thus, the
first partial derivatives of f with respect to
ϕβ exp(µ̄)
0
2 (exp(µ̄)−ϕ)2
C̄ss
−λss
0
D1,3 fss (st+1 , st ) =
βκ
0
0
0
D7,7 fss (st+1 , st ) =
2
− C̃exp(2µ̄)+βϕ
2 (exp(µ̄)−ϕ)2
ss
0
0
−1
all its variables evaluated at the steady-state are
0
0 −1
0
λss ψ (st ) 0 −1
1
, D4,6 fss (st+1 , st ) =
η
−κ
0 − 2
0
λ̃ss
0
0
1
0
, D8,8 fss (st+1 , st ) =
0
λss σ
D10,10 fss (st+1 , st ) =
0
0
, D11,12 fss (st+1 , st ) =
and D13,14 fss =
for all st+1 and st .
exp(µ̄)ϕ
2 (exp(µ̄)−ϕ)2
C̃ss
− C̃
exp(µ̄)ϕ
2
ss (exp(µ̄)−ϕ)
0
0
0
0
κ
η
ϕ
σ
µ̄
µ̂ (1)
µ̂ (2)
0.9976 161 10 0.7 0.0025 0.005 0.0025 −0.0025
0
0
0
b (1) ψ
b (2)
ψ
3.1
0.9
p1,1
p2,2
0.90 0.90
All parameters other than ϕ (the degree of habit formation) are similar to those in the previous
New-Keynesian example. The steady-state values are X̃ss = C̃ss = 0.904957 and λ̃ss = 1.11111.
55
0
0
0
, D9,9 fss (st+1 , st ) =
0
0
0
0
exp(µ̄)βϕ
0
C̃ss (exp(µ̄)−ϕ)2
−λss
0
,
0
0
0
0
Unique Stable Approximation Consider the first parameter configuration
β
Consequently the numerical values of the matrices of first partial derivatives of f with respect
to all its variables evaluated at the steady-state are
0
0
9.21159 0
1.11111ψ (st )
−1.1111
0
1
, D4,6 fss (st+1 , st ) =
D1,3 fss (st+1 , st ) =
160.614
−161.
0
0
0
0
0
0
0
−1
0 −1
,
0 −8.1
1
0
0
0
0
0
, D9,9 fss (st+1 , st ) = ,
0
0
0
0
0
−1
8.33609 0
0
−1.11111 0
0.00277778
,
D10,10 fss (st+1 , st ) =
, D11,12 fss (st+1 , st ) =
0
0
0
0
0
0
−8.35615 0
0
0
and D13,14 fss (st+1 , st ) =
0
0
0
0
D7,7 fss (st+1 , st ) =
−19.6731
, D8,8 fss (st+1 , st ) =
9.23375
for all st+1 and st .
Using these matrices, we can solve the quadratic system for {D1,1 gss (st ) , D1,1 hss (st )}nsts=1 .
56
In this example there are a total of 16 solutions. To conserve space, we report D1,1 hss (st ) only:
D1,1 hss (1)
D1,1 hss (2)
(1)
0.69651
0.69651
(2)
1.43919
1.43919
(3)
X0.79309
1.5799
(4)
1.5799
0.79309
(5)
0.7613 − 0.0895i
1.1350 − 0.1129i
(6)
0.7613 + 0.0895i
1.1350 + 0.1129i
(7)
1.0207 − 0.5488i
1.0926 − 0.0218i
(8)
1.0207 + 0.5488i
1.0926 + 0.0218i
(9)
1.1188 − 0.4326i
1.0051 + 0.0997i
(10)
1.1188 + 0.4326i
1.0051 − 0.0997i
(11)
1.5474 − 0.04355i 1.1735 − 0.1404i
(12)
1.5474 + 0.04355i
1.1735 + 0.1404i
(13)
1.1021 + 0.3496i
0.8230 + 0.0927i
(14)
1.1021 − 0.3496i
0.8230 − 0.0927i
(15)
1.1786 + 0.3971i
1.6102 − 0.0786i
(16)
1.1786 − 0.3971i
1.6102 + 0.0786i
Mathematica finds 16 solutions in approximately 5 seconds. The only solution that produces
a stable approximation is (1). Given Solution (1), we calculate the rest of the matrices needed
b̃ =
to obtain the second-order approximation through solving systems of linear equations. Let C
t
h
i⊺
b̃
b̃ = X̃ − X̃ , λ
b̃
C̃t − C̃ss , Π̂t = Πt − Πss , X
. The
C
t
t
ss
t = λ̃t − λss , and define St =
t−1 εt χ
57
second-order approximation is
b̃
0.69651
−0.0002
0.00045
C
t
Π̂t
0.
−0.0001 0.00021
St
=
b̃
X t 0.69651 −0.0002 0.00045
b̃
0.
0.00239 −0.0069
λt
0 0 −0.0020 0 −2 × 10−7 −1 × 10−6 −0.0020 −1 × 10−6 5 × 10−5
0 0 0.00020 0 2 × 10−7 −4 × 10−7 0.00020 −4 × 10−7 −6 × 10−5
+
0 0 −0.0020 0 −2 × 10−7 −1 × 10−6 −0.0020 −1 × 10−6 5 × 10−5
0 0 0.00143 0 6 × 10−6 −1 × 10−5 0.00143 −1 × 10−5 −0.0003
for st = 1 and
b̃
C
t
0.69651 −0.0002
Π̂t
0
−0.0001
=
b̃
X t 0.69651 −0.0002
b̃
0
0.00261
λt
−0.0005
(St ⊗ St )
−0.0033
St
−0.0005
0.00707
0 0 0.00151 0 −2 × 10−7 1 × 10−6 0.00151 1 × 10−6 −6 × 10−5
0 0 0.00165 0 3 × 10−7 6 × 10−6 0.00165 6 × 10−6 −0.0002
+
0 0 0.00151 0 −2 × 10−7 1 × 10−6 0.00151 1 × 10−6 −6 × 10−5
0 0 0.00156 0 7 × 10−6 2 × 10−5 0.00156 2 × 10−5 0.00036
(St ⊗ St ) .
for st = 2. Again, the first-order approximation can be easily obtained from the above expressions.
Accuracy: Euler Equation Errors As in the previous examples, we compare the accuracies of the first-order and second-order approximations. The Euler equation error at point
58
C̃t−1 , εt ; st for an approximation order ∈ {first,second} is
EE
order
Z X
2
C̃t−1 , εt , 1; st = 1 − β
pst ,s′ ×
s′ =1
(C̃t−1 ,εt ,1;st ),χεt+1 ,1;st+1)
(
×
β
Πorder (C̃ order (C̃t−1 ,εt ,1;st ),χεt+1 ,1;st+1)
′
′
ψ
µ (ε ) dε
Rss Πorder (C̃t−1 ,εt ,1;st ) t exp(σεt )
order
exp(µ(s′ ))λ̃
(C̃t−1 ,εt ,1;st )
order
where C̃ order C̃t−1 , εt , χ; st , Πorder C̃t−1 , εt , χ; st , and λ̃
C̃t−1 , εt , χ; st are the approxiorder
λ̃
C̃ order
mations to the policy functions. The unconditional absolute Euler equation error is
Z
order
order
order
EE
=
C̃t−1 , εt , χ; st µ
C̃t−1 , εt , χ; st dC̃t−1 dεt dst
EE
where µorder is the unconditional distribution of the variables implied by the approximated
solution.
Numerical approximation of this integral simulates a 10, 000-period path, with the first 1, 000
periods discarded as a burn-in. For each state C̃t−1 , εt , st along the simulated path, we draw
10,000 normally distributed εt+1 ’s to compute the expectation over εt+1 , and use the transition
values pst ,st+1 to compute the expectation over st+1 . The entire process takes approximately an
hour.
The table below reports the base-10 logarithms of absolute Euler equation errors for the
first-order and the second-order approximations.
Both the first-order and the second-order
approximations produce a high degree of accuracy: a value of −3 implies an error of $1 for each
$1,000.
Unconditional Absolute Euler Equation Errors (log10 )
EE first
−2.9261
EE second
−2.9527
Multiple Stable Approximations As an alternative parameterization, consider the same
b (2) = 0.6.
parameters as above except that ψ
The second regime now has a slightly lower
response by the monetary authority to inflation.
59
The numerical values of the matrices of
partial derivatives of f with respect to all its variables evaluated at steady-state are unchanged
except
0
1.11111ψ (st )
D4,6 fss (st+1 , st ) =
−161.
0
0
−1
0 −1
0 −8.1
1
0
b (st ). Using this new matrix as well as the other matrices,
for all st+1 and st , which depends on ψ
we solve the quadratic system for {D1,1 gss (st ) , D1,1 hss (st )}nsts=1 . There are 16 solutions to the
quadratic system, and to conserve space we present D1,1 hss (st ) only:
D1,1 hss (1)
D1,1 hss (2)
(1)
0.69651
0.69651
(2)
1.43919
1.43919
(3)
0.79309
1.57990
(4)
1.57990
0.79309
(5)
0.65550
1.03904
(6)
1.67928
1.10504
(7)
0.8703 − 0.1497i 1.1985 + 0.0323i
(8)
0.8703 + 0.1497i
(9)
1.1236 − 0.4088i 0.9377 + 0.1268i
(10)
1.1236 + 0.4088i
(11)
1.4751 − 0.1292i 1.2161 − 0.0448i
(12)
1.4751 + 0.1292i
1.2161 + 0.0448i
(13)
1.1062 + 0.3393i
0.7984 + 0.1078i
(14)
1.1062 − 0.3393i 0.7984 − 0.1078i
(15)
1.1817 + 0.3950i
(16)
1.1817 − 0.3950i 1.6177 + 0.0855i
1.1985 − 0.0323i
0.9377 − 0.1268i
1.6177 − 0.0855i
There are two solutions – (1) and (5) – that produce stable approximations. The examples
in the previous New-Keynesian model and the current variant with habit might suggest that the
60
only parameter affecting the uniqueness of a stable approximation is ψ̂ (s), which governs the
monetary authority’s response to inflation. This conclusion is incorrect.
Consider a higher degree of habit persistence, ϕ = 0.9, while keeping ψ̂ (1) = 3.1 and
ψ̂ (2) = 0.6. Given these parameter values, the steady-state is C̃ss = X̃ss = 0.918512. There are
16 solutions to the quadratic system. The following table presents D1,1 hss (st ) only:
D1,1 hss (1)
D1,1 hss (2)
(1)
0.89551
0.895511
(2)
1.11937
1.11937
(3)
0.82810
1.05334
(4)
1.47489
1.16828
(5)
1.1194 − 0.0017i 1.1194 + 0.0017i
(6)
1.1194 + 0.0017i
1.1194 − 0.0017i
(7)
1.1104 + 0.3027i
0.9173 + 0.1984i
(8)
1.1104 − 0.3027i 0.9173 − 0.1984i
(9)
1.1445 − 0.3980i 0.9790 + 0.1850i
(10)
1.1445 + 0.3980i
(11)
1.0486 − 0.0839i 1.1523 + 0.0536i
(12)
1.0486 + 0.0839i
1.1523 − 0.0536i
(13)
1.1767 + 0.0742i
1.1435 − 0.0715i
(14)
1.1767 − 0.0742i 1.1435 + 0.0715i
(15)
1.1584 + 0.4344i
(16)
1.1584 − 0.4344i 1.4032 + 0.1185i
0.9790 − 0.1850i
1.4032 − 0.1185i
There is only one solution, which is (1), that produces stable approximations. Remember that
there are two stable approximations when we set ϕ = 0.7 as in the previous example.
61
7
Conclusion
Markov switching has been introduced as an essential ingredient to a large class of models usable
for analyses of structural breaks and regime changes in policy, ranging from backward-looking
models (Hamilton (1989) and Sims and Zha (2006)) to forward-looking rational expectations
models (Clarida et al. (2000), Lubik and Schorfheide (2004), Davig and Leeper (2007), Farmer
et al. (2011)). This paper expands this literature by developing a general methodology for
constructing first-order and second-order approximations to the solutions of MSDSGE models.
We resolve conceptual issues related to the steady-state and certainty equivalence of first-order
approximations; we reduce the potentially very complicated problem to a task of solving a
system of quadratic equations; we propose using Gröbner bases to solve such a system; and we
apply the MSS criterion to verify the existence and uniqueness of a stable approximation to the
solution.
The contribution of this paper is not only theoretical but also practical. We show that after
the quadratic system is successfully dealt with, obtaining up to second-order approximations
is straightforward, as the remaining task involves finding a solution to only a system of linear
equations. Our application to three MSDSGE models illustrates the practical value of our
methodology. It is our hope that the advance made in this paper enables applied researchers to
estimate MSDSGE models by focusing on improving the efficiency of our methods.
References
Aruoba, S., J. Fernandez-Villaverde, and J. Rubio-Ramirez (2006). Comparing Solution Methods
for Dynamic Equilibrium Economies. Journal of Economic Dynamics and Control 30 (12),
2477–2508.
Becker, E., M. Marianari, T. Mora, and C. Treverso (1993). The Shape of the Shape Lemma. In
International Symposium on Symbolic and Algebraic Computation, pp. 129–133. ACM Press.
62
Becker, T., V. Weispfenning, and H. Kredel (1998). Gröbner Bases: A Computational Approach
to Commutative Algebra. Springer-Verlag.
Bi, H. and N. Traum (2012). Estimating Sovereign Default Risk. European Economic Review 102 (3), 161–66.
Bianchi, F. (2010). Regime Switches, Agents’ Beliefs, and Post-World War II US Macroeconomic
Dynamics. Working Paper 12-04, Duke University.
Blake, A. and F. Zampolli (2006). Optimal Monetary Policy in Markov-Switching Models with
Rational Expectations Agents. Working Paper 298, Bank of England.
Buchberger, B. (1998). An Algorithmic Criterion for the Solvability of Algebraic Systems of
Equations. Gröbner Bases and Applications 251, 535–545.
Cho, S. (2011). Characterizing Markov-Switching Rational Expectation Models. Working Paper.
Christiano, L., M. Eichenbaum, and C. Evans (2005). Nominal Rigidities and the Dynamic
Effects of a Shock to Monetary Policy. Journal of Political Economy 113 (1), 1–45.
Clarida, R., J. Gali, and M. Gertler (2000). Monetary Policy Rules and Macroeconomic Stability:
Evidence and Some Theory. Quarterly Journal of Economics 115 (1), 147–180.
Costa, O., M. Fragoso, and R. Marques (2005). Discrete-Time Markov Jump Linear Systems.
Springer.
Cox, D., J. Little, and D. O’Shea (1997). Ideals, Varieties and Algorithms: An Introduction to
Computational Algebraic Geometry and Commutative Algebra (2 ed.). Undergraduate Texts
in Mathematics, Springer-Verlag.
Datta, R. (2010). Finding all Nash Equilibria of a Finite Game using Polynomial Algebra.
Economic Theory 42 (1), 55–96.
Davig, T. and T. Doh (2008). Monetary Policy Regime Shifts and Inflation Persistence. Research
Working Paper 08-16, Federal Reserve Bank of Kansas City.
63
Davig, T. and E. Leeper (2007). Generalizing the Taylor Principle. American Economic Review 97 (3), 607–635.
Davig, T., E. M. Leeper, and T. B. Walker (2010). Unfunded Liabilities and Uncertain Fiscal
Financing. Journal of Monetary Economics 57 (5), 600–619.
Davig, T., E. M. Leeper, and T. B. Walker (2011). Inflation and the Fiscal Limit. European
Economic Review 55 (1), 31–47.
Farmer, R., D. Waggoner, and T. Zha (2009). Understanding Markov-Switching Rational Expectations Models. Journal of Economic Theory 144 (5), 1849–1867.
Farmer, R., D. Waggoner, and T. Zha (2011). Minimal State Variable Solutions to MarkovSwitching Rational Expectations Models. Journal of Economic Dynamics and Control 35 (12),
2150–2166.
Fernandez-Villaverde, J., P. Guerron-Quintana, and J. Rubio-Ramirez (2009). Fortune or Virtue:
Time-Variant Volatilites versus Parameter Drifting in U.S. Data. Working Paper.
Fernández-Villaverde, J. and J. Rubio-Ramirez (2007). Estimating Macroeconomic Models: A
Likelihood Approach. Review of Economic Studies 74 (4), 1059–1087.
Hamilton, J. (1989). A New Approach to the Economic Analysis of Nonstationary Time Series
and the Business Cycle. Econometrica 57 (2), 357–384.
Judd, K. (1998). Numerical Methods in Economics. MIT Press.
Kubler, F. and K. Schmedders (2010a). Competitive Equilibria in Semi-Algebraic Economies.
Journal of Economic Theory 145 (1), 301–330.
Kubler, F. and K. Schmedders (2010b). Tackling Multiplicity of Equilibria with Grbner Bases.
Operations Research 58 (4), 1037–1050.
Leeper, E. and T. Zha (2003).
Modest Policy Interventions.
nomics 50 (8), 1673–1700.
64
Journal of Monetary Eco-
Liu, Z., D. F. Waggoner, and T. Zha (2011). Sources of Macroeconomic Fluctuations: A RegimeSwitching DSGE Approach. Quantitative Economics 2 (2), 251–301.
Lubik, T. and F. Schorfheide (2004). Testing for Indeterminacy: An Application to US Monetary
Policy. American Economic Review 94 (1), 190–217.
Maih, J. (2013). Rationality in a Switching Environment. mimeo.
Rudebusch, G. and E. Swanson (2008). Examining the Bond Premium Puzzle with a DSGE
Model. Journal of Monetary Economics 55 (S1), S111–S126.
Schmitt-Grohe, S. and M. Uribe (2004). Solving Dynamic General Equilibrium Models Using
a Second-Order Approximation to the Policy Function. Journal of Economic Dynamics and
Control 28 (4), 755–775.
Schorfheide, F. (2005). Learning and Monetary Policy Shifts. Review of Economic Dynamics 8,
392–419.
Sims, C. and T. Zha (2006). Were There Regime Switches in US Monetary Policy? American
Economic Review 96 (1), 54–81.
Svensson, L. and N. Williams (2007). Monetary Policy with Model Uncertainty: Distribution
Forecast Targeting. Discussion Paper 6331, CEPR.
van Binsbergen, J., J. Fernandez-Villaverde, R. Koijen, and J. Rubio-Ramirez (2008). Likelihood
Estimation of DSGE Models with Epstein-Zin Preferences. Working Paper.
Woodford, M. (2003). Interest and Prices: Foundations of a Theory of Monetary Policy. Princeton University Press.
8
Appendix A:
This appendix derives, in detail, the two steps described in Section 3.2.
65
8.1
Step 1: Obtaining the derivatives of xt−1
Taking first partial derivatives of G with respect to xt−1 in (10) produces the expression for
D1,nx Gss (st )
D1,ny fss (s′ , st ) D1,nx gss (s′ ) D1,nx hss (st )
Z
+Dny +1,2ny fss (s′ , st ) D1,nx gss (st )
pst ,s′
D1,nx Gss (st ) =
+D2ny +1,2ny +nx fss (s′ , st ) D1,nx hss (st )
s′ =1
+D2ny +nx +1,2(ny +nx ) fss (s′ , st )
ns
X
R
µ (ε′ ) dε′
µ (ε′ ) dε′ = 1 into account, one can simplify the above expression to
′
′
D f (s , st ) D1,nx gss (s ) D1,nx hss (st )
1,ny ss
n
s
′
X
+Dny +1,2ny fss (s , st ) D1,nx gss (st )
pst ,s′
D1,nx Gss (st ) =
′
+D2ny +1,2ny +nx fss (s , st ) D1,nx hss (st )
s′ =1
′
+D2ny +nx +1,2(ny +nx ) fss (s , st )
for all st . Taking
for all st . Rearranging the above expression for each st leads to
ns
X
s′ =1
pst ,s′
D1,nx Gss (st ) =
D1,ny fss (s′ , st ) D1,nx gss (s′ ) + D2ny +1,2ny +nx fss (s′ , st ) D1,nx hss (st )
+Dny +1,2ny fss (s′ , st ) D1,nx gss (st ) + D2ny +nx +1,2(ny +nx ) fss (s′ , st )
(20)
.
Putting together the ns versions of (20), one for each value of st , and equating them to zero as
implied by (11) yields a quadratic system of (ny + nx ) nx ns equations. The number of equations
is exactly the same as the number of unknowns {D1,nx gss (st ) , D1,nx hss (st )}nsts=1 . Writing (20)
in matrix form leads to a system of quadratic equations expressed in (12).
66
8.2
8.2.1
Step 2: Obtaining the derivatives of εt and χ
Obtaining the derivatives of εt
By taking first partial derivatives with respect to εt in (10) we obtain the expression for
Dnx +1,nx +nε Gss (st )
Z
ns
X
pst ,s′
s′ =1
for all st . With
pst ,s′
s′ =1
ns
X
R
Dnx +1,nx +nε Gss (st ) =
′
′
D1,ny fss (s , st ) D1,nx gss (s ) Dnx +1,nx +nε hss (st ) +
Dny +1,2ny fss (s′ , st ) Dnx +1,nx +nε gss (st ) +
D2ny +1,2ny +nx fss (s′ , st ) Dnx +1,nx +nε hss (st ) +
D2(ny +nx )+nε +1,2(ny +nx +nε ) fss (s′ , st )
µ (ε′ ) dε′ = 1, this expression simplifies to
µ (ε′ ) dε′
Dnx +1,nx +nε Gss (st ) =
′
′
′
D1,ny fss (s , st ) D1,nx gss (s ) + D2ny +1,2ny +nx fss (s , st ) Dnx +1,nx +nε hss (st )
Dny +1,2ny fss (s′ , st ) Dnx +1,nx +nε gss (st )
D2(ny +nx )+nε +1,2(ny +nx +nε ) fss (s′ , st )
for all st .
(21)
Putting together the ns versions of (21), one for each value of st , and equating them to zero
as implied by (11) yields a system of (ny + nx ) nε ns equations. The number of equations is the
same as that of unknowns {Dnx +1,nx +nε gss (st ) , Dnx +1,nx +nε hss (st )}nsts=1 . The system is linear
and can be written in matrix form as
h
Θε
Dnx +1,nx +nε gss (1)
..
.
i
Dnx +1,nx +nε gss (ns )
= Ψε .
Φε
Dnx +1,nx +nε hss (1)
..
.
Dnx +1,nx +nε hss (ns )
The solution to this system is given in (13).
67
(22)
8.2.2
Obtaining the derivatives of χ
We obtain the expression for Dnx +nε +1 Gss (st ) by taking first partial derivatives with respect to
χ in (10)
Dnx +nε +1 Gss (st ) =
D1,nx gss (s′ ) Dnx +nε +1 hss (st )
D1,ny fss (s′ , st )
+
′
′
′
+Dnx +1,nx +nε gss (s ) ε + Dnx +nε +1 gss (s )
′
Dny +1,2ny fss (s , st ) Dnx +nε +1 gss (st ) +
Z
n
s
X
′
pst ,s′
µ (ε′ ) dε′
D2ny +1,2ny +nx fss (s , st ) Dnx +nε +1 hss (st ) +
s′ =1
′
′
D2ny +nx +1,2ny +nx +nε fss (s , st ) ε +
D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s′ , st ) Dθss (s′ ) +
′
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s , st ) Dθss (st )
for all st . Note that Dθss (st ) is the derivative of θ (χ, st ) with respect to χ evaluated at χ = 0.
That is,
Dθss (st ) = Dθ (0, st ) = Dj θi (0, st ) 1≤i≤n
for all st .
Using the two equalities
for Dnx +nε +1 Gss (st ) as
R
µ (ε′ ) dε′ = 1 and
R
θ ,j=1
ε′ µ (ε′ ) dε′ = 0, one can simplify the expression
Dnx +nε +1 Gss (st ) =
D1,ny fss (s′ , st ) {D1,nx gss (s′ ) Dnx +nε +1 hss (st ) + Dnx +nε +1 gss (s′ )} +
ns
X
Dny +1,2ny fss (s′ , st ) Dnx +nε +1 gss (st ) + D2ny +1,2ny +nx fss (s′ , st ) Dnx +nε +1 hss (st )
pst ,s′
+D2(ny +nx +nε )+1,2(ny +nx +nε )+nθ fss (s′ , st ) Dθss (s′ ) +
s′ =1
D2(ny +nx +nε )+nθ +1,2(ny +nx +nε +nθ ) fss (s′ , st ) Dθss (st )
for all st .
(23)
Putting together the ns versions of (23), one for each value of st , and equating them to
zero as implied by (11), yields a system of linear (ny + nx ) ns equations. The number of these
equations is the same as that of unknowns {Dnx +nε +1 gss (st ) , Dnx +nε +1 hss (st )}nsts=1 . The linear
68
system can be written in matrix form as
h
Θχ
D
g (1)
nx +nε +1 ss
..
.
i Dn +n +1 gss (ns )
x
ε
Φχ
Dnx +nε +1 hss (1)
..
.
Dnx +nε +1 hss (ns )
The solution to (24) is given in (13).
9
= Ψχ .
(24)
Appendix B: Gröbner bases
In this appendix we give an overview of Gröbner bases and describe how they can be applied to
our problem. See Becker et al. (1998) for a more detailed description and other applications.
We wish to find all the solutions of a system of n polynomials in n variables. Let the
polynomial system under study be
f1 (x1 , . . . , xn ) = 0,
..
.
fn (x1 , . . . , xn ) = 0.
Each equation in this system defines a manifold of dimension (n − 1) in Rn and the set of
solutions of the system is the intersection of these manifolds.
When all the fi ’s are linear, the solution set consists of a linear subspace of Rn . It is well
known that there are three possible outcomes: a unique solution, no solutions, or infinitely
many solutions. More importantly, the set of linear systems with no solution or infinitely many
solutions is of measure zero in the set of all linear systems. When there is a unique solution, it
can be easily found.
When the fi ’s are higher-order polynomials, the solution set is more complicated, but the
intuition from the linear case still holds. For almost all polynomial systems of n equations in n
variables, there are only finitely many solutions and these solutions can be easily found.
69
To describe how solutions are computed for polynomial systems, we need to develop the
concept of an ideal and its Gröbner basis. Given a set of polynomials in n variables, {f1 , . . . , fm },
the ideal generated by {f1 , . . . , fm } is the set of all polynomials of the form
g1 (x1 , . . . , xn )f1 (x1 , . . . , xn ) + · · · + gm (x1 , . . . , xn )fm (x1 , . . . , xn )
(25)
where g1 , . . . , gm vary over all polynomials in n variables. We denote this ideal by hf1 , . . . , fm i.
For our purpose we focus on one important feature of an ideal. The point (a1 , . . . , an ) ∈ Rn is
a zero of the polynomials f1 , . . . , fm if and only if it is a zero of every polynomial in the ideal
hf1 , . . . , fm i. This feature implies that if two different sets of polynomials generate the same
ideal, then they have the same zeros. The goal is to find a generating set for which it is easy to
compute zeros.
Before giving the definition of a Gröbner basis, we must first define what we mean by
the leading term of a polynomial. A polynomial in x1 , . . . , xn is a sum of terms of the form
cxk11 xk22 · · · xknn , where ki is a non-negative integer and c is a non-zero real number. The product
xk11 xk22 · · · xknn is called a monomial in x1 · · · xn . The degree of a term is the sum of its exponents,
k1 + · · · + kn . For polynomials in a single variable, the leading term is defined to be the one
of highest degree. For polynomials of several variables, there may be many terms of the same
degree. Thus, one defines the leading term relative to a monomial ordering. For instance, the
mn
1
lexicographical ordering of monomials implies that xk11 · · · xknn < xm
if and only if there
1 · · · xn
is an i such that ki < mi and kj = mj for j < i. In general, a monomial order must satisfy
1. The monomial x01 x02 · · · x0n = 1 is the smallest monomial.
2. If X, Y , and Z are monomials with X < Y , then XZ < Y Z.
The leading term of a polynomial is the largest term with respect to the monomial ordering.
With these notions in hand, we are ready to define a Gröbner basis. The set {h1 , . . . , hk } is a
Gröbner basis for the ideal hf1 , . . . , fm i if
1. hf1 , . . . , fm i=hh1 , . . . , hk i
70
2. The leading term of any polynomial in hf1 , . . . , fm i is divisible by the leading term of hi
for some i.
Consider the following example. Consider the ideal generated by
{2x1 x2 − x1 , x2 − x3 , x23 − 1}.
Note that the first term of each polynomial is the leading term with respect to the lexicographical
order. Is this generating set a Gröbner basis? The answer is negative because the leading term
in 21 (2x1 x2 − x1 ) − x1 (x2 − x3 ) = x1 x3 − 12 x1 is not divisible by any of the leading terms in the
generating set.
As an illustration, we show how to use Buchberger’s Algorithm to construct a Gröbner basis.
There are many other algorithms, most of which are based on Buchberger’s Algorithm, that can
also be used to construct Gröbner bases. The algorithm begins with constructing S-polynomials.
The polynomial 21 (2x1 x2 −x1 ) −x1 (x2 −x3 ) = x1 x3 − 12 x1 is called the S-polynomial of 2x1 x2 −x1
and x2 − x3 because factors
1
2
and x1 were chosen so that the leading terms would cancel. After
the S-polynomial has been formed, it must be reduced using the elements of the generating set.
The reduction step is illustrated by the following example.
Consider the S-polynomial of 2x1 x2 −x1 and x23 −1, which is 12 x23 (2x1 x2 −x1 )−x1 x2 (x23 −1) =
x1 x2 − 12 x1 x23 . This polynomial is reduced by using the leading terms in the generating set to
eliminate terms in the S-polynomial. In this case the reduction proceeds as follows:
1
1
1
1
1
x1 x2 − x1 x23 ⇒ (x1 x2 − x1 x23 ) − (2x1 x2 − x1 ) = − x1 x23 + x1
2
2
2
2
2
1
1
1
1
1
2
2
2
− x1 x3 + x1 ⇒ (− x1 x3 + x1 ) + x1 (x3 − 1) = 0
2
2
2
2
2
So the reduction of the S-polynomial of 2x1 x2 −x1 and x23 −1 gives the zero polynomial. Readers
should convince themselves that the S-polynomial of 2x1 x2 − x1 and x2 − x3 given above cannot
be reduced further and that the S-polynomial of x2 − x3 and x23 − 1 can be reduced to zero.
Note that the reduction is in general not unique. It can depend on the order in which the terms
are eliminated and on particular elements of the generating set that are used to eliminate the
terms. One can always devise an algorithm to reduce any polynomial in finite steps.
71
Buchberger’s Algorithm proceeds as follows. Successively form the S-polynomials from pairs
of polynomials in the generating set and reduce them. If a reduced non-zero S-polynomial is
obtained, add it to the generating set. Continue until all S-polynomials formed from pairs of
the enlarged generating set can be reduced to zero. This algorithm is guaranteed to terminate
in a Gröbner basis. See Buchberger (1998) or Becker et al. (1998) for details.
Continuing with our example, the reduced S-polynomial of 2x1 x2 −x1 and x2 −x3 is x1 x3 − 12 x1 .
We add it to our generating set to obtain
1
{2x1 x2 − x1 , x2 − x3 , x23 − 1, x1 x3 − x1 }.
2
As discussed above, the S-polynomials of both the pair 2x1 x2 −x1 and x23 −1 and the pair x2 −x3
and x23 −1 are zero. Note also that the S-polynomial of 2x1 x2 −x1 and x1 x3 − 12 x1 reduces to zero,
but the S-polynomial of x2 − x3 and x1 x3 − 21 x1 is x1 x3 (x2 − x3 ) − x2 (x1 x3 − 12 x1 ) = 21 x1 x2 − x1 x23
and reduces to
1
1
1
1
x1 x2 − x1 x23 ⇒ ( x1 x2 − x1 x23 ) − (2x1 x2 − x1 ) = −x1 x23 + x1
2
2
4
4
1
3
1
2
2
2
−x1 x3 + x1 ⇒ (−x1 x3 + x1 ) + x1 (x3 − 1) = − x1 .
4
4
4
We add this non-zero polynomial to our generating set to obtain
3
{2x1 x2 − x1 , x2 − x3 , x23 − 1, 2x1 x3 − x1 , − x1 }.
4
The reader should verify that all S-polynomials of pairs from this generating set will reduce to
zero. Thus we have obtained a Gröbner basis.
Gröbner bases are not unique, because adding any element from the ideal generated by a
Gröbner basis will result in another Gröbner basis. To obtain uniqueness, with respect to the
monomial ordering, we work with a reduced Gröbner basis. A Gröbner basis is said to be reduced
if
1. The coefficient of the leading term of each polynomial in the basis is one.
2. Each polynomial in the basis cannot be further reduced with respect to the other polynomials in the basis.
72
Any Gröbner basis can be easily transformed to a reduced Gröbner basis by first reducing
each polynomial in the basis with respect to the other polynomials in the basis and then dividing
the resultant leading term by the leading coefficient. For instance, the Gröbner basis obtained
above is not reduced because both 2x1 x2 − x1 and 2x1 x3 − x1 can be reduced to zero. Thus
these polynomials must be eliminated to obtain the reduced Gröbner basis
{x1 , x2 − x3 , x23 − 1}.
The reduced basis above is called a Shape basis because it is of the form
{x1 − q1 (xn ), . . . , xn−1 − qn−1 (xn ), qn (xn )},
where q1 , . . . , qn are polynomials in a single variable with the degree of qi strictly less than the
degree of qn for 1 ≤ i ≤ n − 1. Shape bases are particularly useful because it is straightforward
to find all the zeros from this representation. One first finds the values of xn that are zeros of
qn (xn ) and then substitutes each of these values into q1 through qn−1 to obtain the values of x1
through xn−1 .
Not all reduced Gröbner bases are Shape bases. The Shape lemma, below, gives the conditions under which the reduced Gröbner basis is a Shape basis.
Lemma 5 Let f1 , . . . , fn be polynomials in x1 , . . . , xn . The reduced Gröbner basis with respect
to the lexicographical ordering of the ideal hf1 , . . . , fn i is a Shape basis if and only if the following
conditions hold.
1. The system f1 , . . . , fn has only finitely many zeros.
2. If (a1 , . . . , an ) and (b1 , . . . , bn ) are two distinct zeros, then an 6= bn .
3. Each zero is either a simple point or a multiple point of local dimension one.
4. If a zero is a multiple point, then the tangent line at the zero does not contain the hyperplane
xn = 0.
73
The meaning of Conditions 1 and 2 is clear, but Conditions 3 and 4 need further explanation. The point (a1 , . . . , an ) ∈ Rn is a zero of the polynomial system f1 , . . . , fn if and only if
hf1 , . . . , fn i ⊆ hx1 − a1 , . . . , xn − an i. If there exists an i such that
hf1 , . . . , fn i ⊆ hx1 − a1 , . . . , (xi − ai )2 , . . . , xn − an i ⊂ hx1 − a1 , . . . , xn − an i,
then we say the zero is a multiple point; otherwise, the zero is simple. One can verify that the
zero (a1 , . . . , an ) is a multiple point if and only if there exists an i such that ∂fj /∂xi |(a1 ,...,an ) = 0
for all 1 ≤ j ≤ n. The tangent space at the zero (a1 , . . . , an ) is the set of all (x1 , . . . , xn ) ∈ Rn
such that
∂f1
|
∂x1 (a1 ,...,an )
..
.
···
..
.
∂fn
|
∂x1 (a1 ,...,an )
···
∂f1
|
x
∂xn (a1 ,...,an ) 1
..
. = 0.
∂fn
|
xn
∂xn (a1 ,...,an )
..
.
Note that this matrix of partial derivatives is the Jacobian. If the zero is simple, then this
definition of the tangent space corresponds to our usual geometric notion of a tangent space,
but the correspondence breaks down if the zero is a multiple point. The local dimension of a
zero is the dimension of the tangent space. Note that the local dimension is zero if and only if
the Jacobian is of full rank. Thus, if the Jacobian is of full rank, then the zero will be simple.
The converse, however, is not necessarily true.
One can verify that if the reduced Gröbner basis is a Shape basis, then Conditions 1-4 will
hold. The converse is also true, but the verification requires much more work. See Becker et al.
(1993) for details. If the Jacobian at each zero is of full rank, then each zero is isolated and
there can only be finitely many zeros. Thus, Conditions 1-4 hold in this case. Since the set
of polynomial systems of n equations in n unknowns whose Jacobian is not of full rank is of
measure zero in the set of all such systems, Conditions 1-4 hold almost surely.
In summary, for almost all polynomial systems, there are only finitely many zeros and Buchberger’s Algorithm can be used to find them. While Buchberger’s Algorithm is instructive, it
can be inefficient for certain problems. Active research continues to develop variants of this
algorithm that improve efficiency.
74
10
Appendix C: Iterative algorithm
This section describes an iterative procedure to find solutions of the quadratic system (12).
In general, this method will not provide us with all the solutions to the system.
quadratic system to be solved is
In x
D1,nx gss (1)
A (st )
..
.
D1,nx gss (ns )
for all st , where
A (st ) =
and
Recall the
In x
D1,nx hss (st ) = B (st )
D1,nx gss (st )
h P
ns
′
s′ =1 pst ,s′ D2ny +1,2ny +nx fss (s , st ) pst ,1 D1,ny fss (1, st ) · · · pst ,ns D1,ny fss (ns , st )
B (st ) = −
ns
X
pst ,s′
s′ =1
for all st .
h
D2ny +nx +1,2(ny +nx ) fss (s′ , st ) Dny +1,2ny fss (s′ , st )
i
i
The idea of the algorithm is to guess a set of policy functions and solve for each regime’s policy
functions as in the constant parameter model case using the generalized Schur decomposition.
When the solutions are close to the guesses, then a solution has been found.
Algorithm 6 Let
n
(j)
D1,nx gss
(j)
(st ) , D1,nx hss (st )
o
denote solution at iteration j.
(0)
(0)
1. Set j = 1 and initialize D1,nx gss (st ) = 0ny ×nx and D1,nx hss (st ) = 0nx ×nx .
2. For each st , construct the transformed system
(j−1)
(s′ )
A st , D1,nx gss
ns
s′ =1,s′ 6=st
= B (st )
In x
(j)
D1,nx gss
In x
(j)
D1,nx gss
75
(st )
(st )
D1,nx h(j)
ss (st )
(26)
where
+
Pns
Pns
s′ =1
(j)
(s′ )
A st , D1,nx gss
ns
s′ =1,s′ 6=st
′
pst ,s′ D2ny +1,2ny +nx fss (s , st )
(0)
′
′
s′ =1,s′ 6=st pst ,s′ D1,ny fss (s , st ) D1,nx gss (s )
=
pst ,st D1,ny fss (st , st )
3. The set of systems (26) is identical to a constant parameter model case.
When find-
(j)
ing D1,nx hss (st ), use the ”most stable” generalized eigenvalues, those with the smallest
modulus.
4. Check max
n
(j)
(j−1)
D1,nx hss (st ) − D1,nx hss
(j)
(j−1)
(st ) , D1,nx gss (st ) − D1,nx gss
(st )
o
< crit.
If yes, then stop and check for MSS. If no, then set j =⇒ j + 1 and return to step 2.
11
Appendix D: Second partial derivatives
This appendix derives, in detail, the six steps described in Section 5.2. Let us first define
i
Dk Dj fss
(st+1 , st ) =
Dk Dj f i yss , yss , xss , xss , 0nε , 0nε , θ1 , b
θ2 (st+1 ) , θ1 , b
θ2 (st )
for all st+1 and st and 1 ≤ i ≤ ny + nx and 1 ≤ j, k ≤ 2 (ny + nx + nε + nθ ) and
i
(st+1 , st )
D D f i (s , s ) . . . Dn1 Dm2 fss
n1 m1 ss t+1 t
.
.
.
i
..
..
..
Hn1 ,n2 ;m1 ,m2 fss
(st+1 , st ) =
i
i
Dn2 Dm1 fss
(st+1 , st ) . . . Dn2 Dm2 fss
(st+1 , st )
for all st+1 and st and 1 ≤ n1 , m1 < n2 , m2 ≤ 2 (ny + nx + nε + nθ ) and 1 ≤ i ≤ ny + nx .
76
11.1
Step 1: Obtaining the derivatives of xt−1 twice
The first equation is the derivative with respect to xt−1 twice.
H1,nx ;1,nx Giss
(st ) =
ns
X
pst ,s′ ×
s′ =1
i
(s′ , st ) D1,nx gss (s′ )
H1,ny ;1,ny fss
D1,nx hss (st )
i
′
+2H1,ny ;2ny +1,2ny +nx fss (s , st )
(D1,nx gss (s′ ) D1,nx hss (st ))⊺
i
(s′ , st ) D1,nx gss (st )
+2H1,ny ;ny +1,2ny fss
i
+2H1,ny ;2ny +nx +1,2(ny +nx ) fss
(s′ , st )
i
(s′ , st ) D1,nx gss (st )
Hny +1,2ny ;ny +1,2ny fss
i
′
+D1,nx gss (st )⊺ +2Hny +1,2ny ;2ny +1,2ny +nx fss
(s , st ) D1,nx hss (st )
i
′
+2Hny +1,2ny ;2ny +nx +1,2(ny +nx ) fss (s , st )
i
(s′ , st )
H2ny +1,2ny +nx ;2ny +1,2ny +nx fss
D1,nx hss (st )
⊺
′
i
′
h
(s
)
+D
+D1,ny fss (s , st ) H1,nx ;1,nx gss (s )
t
1,nx ss
i
+2H2ny +1,2ny +nx ;2ny +nx +1,2(ny +nx ) fss
(s′ , st )
′
′
i
D1,ny fss (s , st ) D1,nx gss (s )
H1,nx ;1,nx hss (st )
+
i
′
+D2ny +1,2ny +nx fss (s , st )
i
(s′ , st ) H1,nx ;1,nx gss (st )
+Dny +1,2ny fss
i
+H2ny +nx +1,2(ny +nx );2ny +nx +1,2(ny +nx ) fss
(s′ , st )
The condition H1,nx ;1,nx Giss (st ) = 0nx ×nx for 1 ≤ i ≤ ny + nx implies a linear system of equations
that
determines
n
nx
y
i
(st )}i=1
{H1,nx ;1,nx gss
, {H1,nx ;1,nx hiss (st )}i=1
the
ns
.
st =1
77
solution
for
11.2
Step 2: Obtaining the derivatives of xt−1 and εt
The second equation is the derivative with respect to xt−1 and εt .
H1,nx ;nx +1,nx +nε Giss
(st ) =
ns
X
pst ,s′ ×
s′ =1
i
(s′ , st ) D1,nx gss (s′ )
H1,ny ;1,ny fss
Dnx +1,nx +nε hss (st )
i
′
+2H1,ny ;2ny +1,2ny +nx fss (s , st )
(D1,nx gss (s′ ) D1,nx hss (st ))⊺
i
(s′ , st ) Dnx +1,nx +nε gss (st )
+H1,ny ;ny +1,2ny fss
i
+H1,ny ;2(ny +nx )+nε +1,2(nx +ny +nε ) fss
(s′ , st )
i
(s′ , st ) D1,nx gss (s′ )
Hny +1,2ny ;1,ny fss
Dnx +1,nx +nε hss (st )
i
′
+Hny +1,2ny ;2ny +1,2ny +nx fss (s , st )
⊺
g
(s
)
+D
ss
t
1,n
x
i
′
+Hny +1,2ny ;ny +1,2ny fss (s , st ) Dnx +1,nx +nε gss (st )
i
+Hny +1,2ny ;2(ny +nx )+nε +1,2(nx +ny +nε ) fss
(s′ , st )
i
′
H2ny +1,2ny +nx ;ny +1,2ny fss (s , st ) Dnx +1,nx +nε gss (st )
i
′
+H2ny +1,2ny +nx ;2(ny +nx )+nε +1,2(nx +ny +nε ) fss (s , st )
⊺
h
(s
)
+D
ss
t
1,n
x
′
i
f (s , st )
H
+ 2ny +1,2ny +nx ;2ny +1,2ny +nx ss
Dnx +1,nx +nε hss (st )
i
(s′ , st ) H1,nx ;1,nx gss (s′ )
+D1,ny fss
i
(s′ , st ) D1,nx gss (s′ )
H2ny +nx +1,2(nx +ny );1,ny fss
Dnx +1,nx +nε hss (st )
+
i
′
+H2ny +nx +1,2(ny +nx );2ny +1,2ny +nx fss (s , st )
′
i
′
g
(s
)
f
(s
,
s
)
D
D
t
1,nx ss
1,ny ss
H1,nx ;nx +1,nx +nε hss (st )
+
i
′
+D2ny +1,2ny +nx fss (s , st )
i
(s′ , st ) Dnx +1,nx +nε gss (st )
+H2ny +nx +1,2(ny +nx );ny +1,2ny fss
i
(s′ , st ) H1,nx ;nx +1,nx +nε gss (st )
+Dny +1,2ny fss
i
+H2ny +nx +1,2(nx +ny );2(nx +ny )+nε +1,2(nx +ny +nε ) fss
(s′ , st )
Conditional on the solution from Step 1, the condition H1,nx ;nx +1,nx +nε Giss (st ) = 0nx ×nε
for 1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution
ny
nx ns
i
.
(st )}i=1
{H1,nx ;nx +1,nx +nε gss
, {H1,nx ;nx +1,nx +nε hiss (st )}i=1
st =1
78
11.3
Step 3: Obtaining the derivatives of xt−1 and χ
The third equation is the derivative with respect to xt−1 and χ.
H1,nx ;nx +nε +1 Giss (st ) =
ns
X
pst ,s′ ×
s′ =1
Dnx +nε +1 gss (s′ )
i
H1,ny ;1,ny fss
(s′ , st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
′
+H1,ny ;ny +1,2ny fss (s , st ) Dnx +nε +1 gss (st )
(D1,nx gss (s′ ) D1,nx hss (st ))⊺
i
(s′ , st ) Dnx +nε +1 hss (st )
+H1,ny ;2ny +1,2ny +nx fss
i
+H1,ny ;2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss
(s′ , st ) Dθss (s′ )
i
+H1,ny ;2(nx +ny +nε +nθ )+1,2(nx +ny +nε +nθ ) fss
(s′ , st ) Dθss (st )
′
D
g
(s
)
nx +nε +1 ss
i
H2ny +1,2ny +nx ;1,ny fss
(s′ , st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
(s′ , st ) Dnx +nε +1 hss (st )
+H2ny +1,2ny +nx ;2ny +1,2ny +nx fss
i
′
f
(s
,
s
)
D
g
(s
)
+H
t
n
+n
+1
ss
t
2n
+1,2n
+n
;n
+1,2n
x
ε
y
y
x
y
y ss
+D1,nx hss (st )⊺
′
′
i
+H2ny +1,2ny +nx ;2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss (s , st ) Dθ ss (s )
i
+H2ny +1,2ny +nx ;2(nx +ny +nε )+nθ +1,2(nx +ny +nε +nθ ) fss
(s′ , st ) Dθ ss (st )
′
H
g
(s
)
1,nx ;nx +nε +1 ss
i
′
f
(s
,
s
)
D
t
1,ny ss
+H1,nx ;1,nx gss (s′ ) Dnx +nε +1 hss (st )
′
D
g
(s
)
nx +nε +1 ss
i
Hny +1,2ny ;1,ny fss
(s′ , st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
′
f
(s
,
s
)
D
g
(s
)
+H
t
nx +nε +1 ss
t
ny +1,2ny ;ny +1,2ny ss
⊺
g
(s
)
+D
t
1,nx ss
i
(s′ , st ) Dnx +nε +1 hss (st )
+Hny +1,2ny ;2ny +1,2ny +nx fss
i
′
′
+Hny +1,2ny ;2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss (s , st ) Dθ ss (s )
i
(s′ , st ) Dθ ss (st )
+Hny +1,2ny ;2(nx +ny +nε +nθ )+1,2(nx +ny +nε +nθ ) fss
′
D
g
(s
)
nx +nε +1 ss
i
(s′ , st )
+H2ny +nx +1,2(nx +ny );1,ny fss
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
(s′ , st ) Dnx +nε +1 gss (st )
+H2ny +nx +1,2(ny +nx );ny +1,2ny fss
i
(s′ , st ) Dnx +nε +1 hss (st )
+H2ny +nx +1,2(ny +nx );2ny +1,2ny +nx fss
i
(s′ , st ) Dθss (s′ )
+H2ny +nx +1,2(nx +ny );2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss
i
+H2ny +nx +1,2(nx +ny );2(nx +ny +nε )+nθ +1,2(nx +ny +nε +nθ ) fss
(s′ , st ) θχ (s)
f i (s′ , st ) D1,nx gss (s′ )
D
1,ny ss
H1,nx ;nx +nε +1 hss (st )
+
i
(s′ , st )
+D2ny +1,2ny +nx fss
i
(s′ , st ) H1,nx ;nx +nε +1 gss (st )
+Dny +1,2ny fss
79
Given the solution from Step 2, the condition Hnx +1,nx +nε ;nx +nε +1 Giss (st )
=
0nx
for 1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution
ny
nx ns
i
.
{Hnx +1,nx +nε ;nx +nε +1 gss
(st )}i=1
, {Hnx +1,nx +nε ;nx +nε +1 hiss (st )}i=1
st =1
11.4
Step 4: Obtaining the derivatives of εt−1 twice
The fourth equation is the derivative with respect to εt−1 twice.
Hnx +1,nx +nε ;nx +1,nx +nε Giss (st )
=
ns
X
pst ,s′ ×
s′ =1
i (s′ , s ) D
′
H1,ny ;1,ny fss
t
1,nx gss (s )
Dnx +1,nx +nε hss (st )
′
i
+2H1,ny ;2ny +1,2ny +nx fss (s , st )
[D1,nx gss (s′ ) Dnx +1,nx +nε hss (st )]⊺
i (s′ , s ) D
+2H1,ny ;ny +1,2ny fss
t
nx +1,nx +nε gss (st )
i (s′ , s )
+2H1,ny ;2(ny +nx )+nε +1,2(nx +ny +nε ) fss
t
′
i
(s
,
s
)
H
f
t
2ny +1,2ny +nx ;2ny +1,2ny +nx ss
Dnx +1,nx +nε hss (st )
⊺
i
′
′
+D
h
(s
)
+D
f
(s
,
s
)
H
g
(s
)
n
+1,n
+n
ss
t
x
x
ε
1,ny ss
t
1,nx ;1,nx ss
i (s′ , s )
+2H
f
t
2ny +1,2ny +nx ;2(nx +ny )+nε +1,2(nx +ny +nε ) ss
′
i
(s
,
s
)
D
H
g
(s
)
f
t
nx +1,nx +nε ss t
ny +1,2ny ;ny +1,2ny ss
⊺
i (s′ , s ) D
+D
g
(s
)
+2H
f
h
(s
)
nx +1,nx +nε ss t
ny +1,2ny ;2ny +1,2ny +nx ss
t
nx +1,nx +nε ss t
i
′
+2Hny +1,2ny ;2(ny +nx )+nε +1,2(nx +ny +nε ) fss (s , st )
′
′
i
D1,ny fss (s , st ) D1,nx gss (s )
Hnx +1,nx +nε ;nx +1,nx +nε hss (st )
+
i
′
+D2ny +1,2ny +nx fss (s , st )
i (s′ , s ) H
+Dny +1,2ny fss
t
nx +1,nx +nε ;nx +1,nx +nε gss (st )
i (s′ , s )
+H2(ny +nx )+nε +1,2(nx +ny +nε );2(ny +nx )+nε +1,2(nx +ny +nε ) fss
t
Given the solution from Step 1, the condition Hnx +1,nx +nε ;nx +1,nx +nε Giss (st ) = 0nε ×nε
for 1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution
ny
nx ns
i
.
(st )}i=1
{Hnx +1,nx +nε ;nx +1,nx +nε gss
, {Hnx +1,nx +nε ;nx +1,nx +nε hiss (st )}i=1
st =1
80
11.5
Step 5: Obtaining the derivatives of εt−1 and χ
The fifth equation is the derivative with respect to εt−1 and χ.
Hnx +1,nx +nε ;nx +nε +1 Giss (st ) =
ns
X
pst ,s′ ×
s′ =1
Dnx +nε +1 gss (s′ )
i
H1,ny ;1,ny fss
(s′ , st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
′
+H1,ny ;ny +1,2ny fss (s , st ) Dnx +nε +1 gss (st )
(D1,nx gss (s′ ) Dnx +1,nx +nε hss (st ))⊺
i
(s′ , st ) Dnx +nε +1 hss (st )
+H1,ny ;2ny +1,2ny +nx fss
i
+H1,ny ;2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss
(s′ , st ) Dθ ss (s′ )
i
+H1,ny ;2(nx +ny +nε +nθ )+1,2(nx +ny +nε +nθ ) fss
(s′ , st ) Dθ ss (st )
′
D
g
(s
)
nx +nε +1 ss
i
H2ny +1,2ny +nx ;1,ny fss
(s′ , st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
(s′ , st ) Dnx +nε +1 gss (st )
+H2ny +1,2ny +nx ;ny +1,2ny fss
i
′
f
(s
,
s
)
D
h
(s
)
+H
t
n
+n
+1
ss
t
2n
+1,2n
+n
;2n
+1,2n
+n
x
ε
y
y
x
y
y
x ss
+Dnx +1,nx +nε hss (st )⊺
′
′
i
+H2ny +1,2ny +nx ;2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss (s , st ) Dθ ss (s )
i
+H2ny +1,2ny +nx ;2(nx +ny +nε )+nθ +1,2(nx +ny +nε +nθ ) fss
(s′ , st ) Dθss (st )
′
H
g
(s
)
1,nx ;nx +nε +1 ss
i
′
f
(s
,
s
)
D
t
1,ny ss
+H1,nx ;1,nx gss (s′ ) Dnx +nε +1 hss (st )
′
D
g
(s
)
nx +nε +1 ss
i
Hny +1,2ny ;1,ny fss
(s′ , st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
′
f
(s
,
s
)
D
g
(s
)
+H
t
nx +nε +1 ss
t
ny +1,2ny ;ny +1,2ny ss
⊺
g
(s
)
+D
t
nx +1,nx +nε ss
i
(s′ , st ) Dnx +nε +1 hss (st )
+Hny +1,2ny ;2ny +1,2ny +nx fss
i
′
′
+Hny +1,2ny ;2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss (s , st ) Dθss (s )
i
(s′ , st ) Dθss (st )
+Hny +1,2ny ;2(nx +ny +nε +nθ )+1,2(nx +ny +nε +nθ ) fss
′
D
g
(s
)
nx +nε +1 ss
i
(s′ , st )
+H2(nx +ny )+nε +1,2(nx +ny +nε );1,ny fss
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i
(s′ , st ) Dnx +nε +1 gss (st )
+H2(nx +ny )+nε +1,2(nx +ny +nε );ny +1,2ny fss
i
(s′ , st ) Dnx +nε +1 hss (st )
+H2(nx +ny )+nε +1,2(nx +ny +nε );2ny +1,2ny +nx fss
i
(s′ , st ) Dθss (s′ )
+H2(ny +nx )+nε +1,2(nx +ny +nε );2(nx +ny +nε )+1,2(nx +ny +nε )+nθ fss
i
+H2(ny +nx )+nε +1,2(nx +ny +nε );2(nx +ny +nε )+nθ +1,2(nx +ny +nε +nθ ) fss
(s′ , st ) Dθss (st )
f i (s′ , st ) D1,nx gss (s′ )
D
1,ny ss
Hnx +1,nx +nε ;nx +nε +1 hss (st )
+
i
(s′ , st )
+D2ny +1,2ny +nx fss
i
(s′ , st ) Hnx +1,nx +nε ;nx +nε +1 gss (st )
+Dny +1,2ny fss
81
Given the solution from Steps 1 and 3, the condition Hnx +1,nx +nε ,nx +nε ;nx +nε +1 Giss (st ) = 0nε
for 1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution
ny
nx ns
i
.
{Hnx +1,nx +nε ;nx +nε +1 gss
(st )}i=1
, {Hnx +1,nx +nε ;nx +nε +1 hiss (st )}i=1
st =1
11.6
Step 6: Obtaining the derivatives of χ twice
The sixth equation is the derivative with respect to χ twice.
Hnx +nε +1;nx +nε +1 Giss (st ) =
Dnx +nε +1 gss (s′ )
⊺
ns
X
pst ,s′ ×
s′ =1
Dnx +nε +1 gss (s′ )
t)
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
+D1,nx gss (s′ ) Dnx +nε +1 hss (st )
i (s′ , s ) D
H1,ny ;ny +1,2ny fss
t
nx +nε +1 gss (st )
⊺
i (s′ , s ) D
+H1,ny ;2ny +1,2ny +nx fss
Dnx +nε +1 gss (s′ )
t
nx +nε +1 hss (st )
+2
i
′
′
′
+H1,ny ;2(nx +ny +nε )+1,2(nx +ny +nε )+n fss (s , st ) Dθ ss (s )
+D1,nx gss (s ) Dnx +nε +1 hss (st )
θ
i
′
+H1,ny ;2(nx +ny +nε +n )+1,2(nx +ny +nε +n ) fss (s , st ) Dθss (st )
θ
θ
i (s′ , s ) D
Hny +1,2ny ;ny +1,2ny fss
t
nx +nε +1 gss (st )
i (s′ , s ) D
+2Hny +1,2ny ;2ny +1,2ny +nx fss
t
nx +nε +1 hss (st )
⊺
+Dnx +nε +1 gss (st )
i (s′ , s ) Dθ
′)
f
(s
+2H
t
ss
ss
n
+1,2n
;2
n
+n
+n
+1,2
n
+n
+n
+n
(
)
(
)
y
y
x
y
ε
x
y
ε
θ
i (s′ , s ) Dθ
+2Hny +1,2ny ;2(nx +ny +nε +n )+1,2(nx +ny +nε +n ) fss
t
ss (st )
θ
θ
i (s′ , s ) D
H2ny +1,2ny +nx ;2ny +1,2ny +nx fss
t
nx +nε +1 hss (st )
i (s′ , s ) Dθ
′)
+2H
f
(s
t
ss
ss
2n
+1,2n
+n
;2
n
+n
+n
+1,2
n
+n
+n
+n
(
)
(
)
y
y
x
x
y
ε
x
y
ε
θ
⊺
i (s′ , s ) Dθ
+D
h
(s
)
+2H
f
(s
)
ss
t
n
+n
+1
t
ss
t
x
ε
2ny +1,2ny +nx ;2(nx +ny +nε )+nθ +1,2(nx +ny +nε +nθ ) ss
i (s′ , s ) H
′) D
+
D
f
g
(s
h
(s
)
t
ss
ss
t
1,n
1,n
;1,n
n
+n
+1
y ss
x
x
x
ε
i
′
′
+2 D1,ny fss (s , st ) H1,nx ;nx +nε +1 gss (s )
i (s′ , s )
H2(nx +ny )+1,2(nx +ny )+nε ;2(nx +ny )+1,2(nx +ny )+nε fss
t
i (s′ , s ) H
′
′
+D1,ny fss
t
nx +1,nx +nε ;nx +1,nx +nε gss (s )
⊺
′
ε
+ (ε )
′
′
i
′
′
+Dnx +1,nx +nε gss (s ) H1,ny ;1,ny fss (s , st ) Dnx +1,nx +nε gss (s )
i (s′ , s )
+2Dnx +1,nx +nε gss (s′ )′ H1,ny ;2(ny +nx )+1,2(ny +nx )+nε fss
t
i (s′ , s ) Dθ
′
H2(nx +ny +nε )+1,2(nx +ny +nε )+n ;2(nx +ny +nε )+1,2(nx +ny +nε )+n fss
t
ss (s )
θ
θ
+Dθss (s′ )⊺
i (s′ , s ) Dθ
+2H2(nx +ny +nε )+1,2(nx +ny +nε )+n ;2(nx +ny +nε )+n +1,2(nx +ny +nε +n ) fss
t
ss (st )
θ
θ
θ
i (s′ , s ) Dθ
+Dθ ss (st )⊺ H2(nx +ny +nε )+n +1,2(nx +ny +nε +n );2(nx +ny +nε )+n +1,2(nx +ny +nε +n ) fss
t
ss (st )
θ
θ
θ
θ
i (s′ , s ) D
′)
D
f
g
(s
t
1,ny ss
1,nx ss
Hnx +nε +1;nx +nε +1 hss (st )
+
i (s′ , s )
+D
fss
t
2n
+1,2n
+n
y
y
x
i (s′ , s ) H
′
+D1,ny fss
t
nx +nε +1;nx +nε +1 gss (s )
′
i
+Dny +1,2ny fss (s , st ) Hnx +nε +1;nx +nε +1 gss (st )
′
i (s′ , s ) Hθ
+D2(nx +ny +nε )+1,2(nx +ny +nε )+n fss
t
ss (s )
θ
i (s′ , s ) Hθ
+D2(nx +ny +nε )+n +1,2(nx +ny +nε +n ) fss
t
ss (st )
i
H1,ny ;1,ny fss
θ
(s′ , s
θ
Given the solutions from Steps 1, 3 and 4, the condition Hnx +nε +1;nx +nε +1 Giss (st ) = 0
82
for 1 ≤ i ≤ ny + nx implies a linear system of equations that determines the solution
ny
nx ns
i
.
{Hnx +nε +1;nx +nε +1 gss
(st )}i=1
, {Hnx +nε +1;nx +nε +1 hiss (st )}i=1
st =1
83
@FedFRASER