% ----------------------------------------------------------------
% AMS-LaTeX Paper ************************************************
% **** -----------------------------------------------------------
\documentclass[12pt,reqno,oneside]{amsart}
\usepackage{srcltx} % SRC Specials: DVI [Inverse] Search
\usepackage[left=1.25in,top=1.2in]{geometry}
%\usepackage{harvard}
\usepackage{natbib}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{mathptmx}
\usepackage[scaled]{helvet}
\usepackage[T1]{fontenc}
%\linespread{1.05} % Palatino needs more leading
\usepackage[scaled]{helvet}
\usepackage{courier}
\normalfont
\usepackage{hyperref}
%\usepackage{calc}
\renewcommand{\baselinestretch}{1.1}
% ----------------------------------------------------------------
\vfuzz2pt % Don't report over-full v-boxes if over-edge is small
\hfuzz2pt % Don't report over-full h-boxes if over-edge is small
% THEOREMS -------------------------------------------------------
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}
\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}
%\numberwithin{equation}{section}
% MATH -----------------------------------------------------------
% \divby command defined below produces nice-looking, right-sized "a/b" template.
\newcommand{\divby}[2]{#1 \mathord{\left/ \vphantom{#1 #2} \right.}
\kern-\nulldelimiterspace #2}
\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\Real}{\mathbb R}
\newcommand{\eps}{\varepsilon}
\newcommand{\To}{\xrightarrow}
\newcommand{\Cstar}{{C^\ast}}
\newcommand{\tcol}{\text{:}}
%e.g.: \[X_t \To[t\to \infty]{q.m.} Z\]
%\newcommand{\BX}{\mathbf{B}(X)}
\newcommand{\A}{\mathcal{A}}
\DeclareMathOperator{\Span}{span}
\DeclareMathOperator{\diag}{diag}
% ----------------------------------------------------------------
\begin{document}
%\setlength{\parskip}{.5\baselineskip plus 20pt - .5\baselineskip}
\title[Linear RE Models]{Solving Linear Rational Expectations Models}%
\author{Christopher A. Sims}%
\address{Department of Economics, Yale University}%
\email{sims@econ.yale.edu}%
\date{August 1995. Revised February 1996, September 1997, January 2000.}
\thanks{This paper has benefited from so many comments from so many people over its
overlong gestation that I cannot list everyone. Particularly important were
Jinill Kim, Gary Anderson, and an anonymous referee. \copyright 1995, 1996,
1997, 2000 by Christopher A. Sims. This material may be reproduced for
educational and research purposes so long as the copies are not sold, even to
recover costs, the document is not altered, and this copyright notice is
included in the copies.} \maketitle
\section{General form of the models}
The models we are interested in can be cast in the form
\begin{equation}\label{eq:mainEq}
\Gamma _0 y(t) = \Gamma _1 y(t - 1) + C + \Psi z(t) + \Pi \eta(t)
\end{equation}
$t = 1,\ldots,T$, where $C$ is a vector of constants, $z(t)$
is an exogenously evolving, possibly serially correlated, random
disturbance, and $\eta(t)$ is an expectational error, satisfying $E_t \eta (t + 1) =
0$, all $t$. The $\eta(t)$ terms are not given exogenously, but instead
are treated as determined as part of the model solution. Models
with more lags, or with lagged expectations, or with expectations
of more distant future values, can be accommodated in this
framework by expanding the y vector. This paper's analysis is
similar to that of \citet{BlanchKahn} with four important
differences:
\begin{enumerate}\renewcommand{\theenumi}{\roman{enumi}}
\item They assume regularity conditions as they proceed that leave
some models encountered in practice outside the range of their
analysis, while this paper covers all linear models with
expectational error terms.
\item They require that the analyst specify which elements of the
y vector are predetermined, while this paper recognizes that the
structure of the $ \Gamma _0$, $ \Gamma _1$,
$\Psi$, and $\Pi$ matrices fixes the list of predetermined
variables. Our approach therefore handles automatically
situations where linear combinations of variables, not individual
variables, are predetermined.
\item This paper makes an explicit extension to continuous time,
which raises some distinct analytic difficulties.
\item They assume that boundary conditions at infinity are given
in the form of a maximal rate of growth for any element of the $y$
vector, whereas this paper recognizes that in general only
certain linear combinations of variables are required to grow at
bounded rates and that different linear combinations may have
different growth rate restrictions.
\end{enumerate}
There are other more recent papers dealing with models like these
\citep{KingWatsonAlgo,KingWatsonTheory,AndersonCT,KleinPSchur} that, like this one,
expand the range of models covered beyond what is covered in Blanchard and
Kahn's paper, particularly to include singular $\Gamma_0$ cases. All these
papers, though, follow Blanchard and Kahn in requiring the specification of
``jump" and ``predetermined" \emph{variables}, rather than recognizing that in
equilibrium models expectational residuals more naturally are attached to
equations. Also, only \citet{AndersonCT} has discussed the continuous
time case.
Less fundamentally, this paper uses a notation in which time
arguments or subscripts relate consistently to the information
structure: variables dated $t$ are always known at $t$. Blanchard
and Kahn's use of a different convention often leads to confusion
in the application of their method to complex models.
An instructive example to illustrate how we get a model into
the form \eqref{eq:mainEq} is a model of overlapping contracts in wage setting
along the lines laid out by Taylor.
\begin{equation}\label{eq:Taylor}
\begin{aligned}
w(t) &= \tfrac{1}
{3}E_t \left[ {W(t) + W(t + 1) + W(t + 2)} \right] - \alpha
(u(t) - u_n ) + \nu (t) \\
W(t) &= \tfrac{1}
{3}\left( {w(t) + w(t - 1) + w(t - 2)} \right) \\
u(t) &= \theta u(t - 1) + \gamma W(t) + \mu + \varepsilon
(t)\text{ ,}
\end{aligned}
\end{equation}
where $E_t \nu (t + 1) = E_t \varepsilon (t + 1) = 0$. To cast
\eqref{eq:Taylor} into the form \eqref{eq:mainEq} requires using the expanded state vector
\begin{equation}
y(t) = \begin{bmatrix}
w(t) \\
w(t - 1) \\
W(t) \\
u(t) \\
E_t W(t + 1) \\
\end{bmatrix}
\end{equation}
With this definition of $y$, \eqref{eq:Taylor} can be written in the matrix
notation of \eqref{eq:mainEq}, with definitional equations added, by defining
\begin{equation}\label{eq:TaylorMats}
\begin{gathered}
\Gamma _0 = {\begin{bmatrix}
0 & 0 & \frac{1}{3} & 0 & \frac{1}{3} \\
-\frac{1}{3} & -\frac{1}{3} & 1 & 0 & 0 \\
0 & 0 & -\gamma & 1 & 0 \\
0 & 1 & 0 & 0 & 0 \\
0 & 0 & 1 & 0 & 0 \\
\end{bmatrix} }\:,\quad \Gamma _1 =
{\begin{bmatrix}
1 & 0 & -\frac{1}{3} & \alpha & 0 \\
0 & \frac{1}{3}& 0 & 0 & 0 \\
0 & 0 & 0 & 0 & \theta & 0 \\
1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1
\end{bmatrix} } \:,\\
C =
{\begin{bmatrix}
\alpha\cdot u_n \\
0 \\
\mu \\
0 \\
0 \\
\end{bmatrix} }\:,\quad
\Psi = {\begin{bmatrix}
1 & 0 \\
0 & 0 \\
0 & 1 \\
0 & 0 \\
0 & 0 \\
\end{bmatrix} }\:,\quad\Pi =
{\begin{bmatrix}
0 & 1 \\
0 & 0 \\
0 & 0 \\
0 & 0 \\
1 & 0 \\
\end{bmatrix} }\:,\\
\quad z(t) =
{\begin{bmatrix}
{\varepsilon (t)} \\
{\nu (t-1)} \\
\end{bmatrix} }\:\,. \\
\end{gathered}
\end{equation}
This example illustrates the principle that we can always get a linear model
into the form \eqref{eq:mainEq} by replacing terms of the form $ E_t x\left( {t + 1}
\right)$ with $ y(t) = E_t x(t + 1)$ and adding to the system an equation
reading $ x(t) = y(t - 1) + \eta (t)$. When terms of the form $ E_t x(t +
s)$ appear, we simply make a sequence of such variable and equation
creations. It is often possible to reach the form \eqref{eq:mainEq} with fewer new
variables by replacing expectations of variables in
equations with actual values of the same variables, while adding $ \eta
(t)$ disturbance terms to the equation. While this approach always
produces valid equations, it may lose information contained in the
original system and thereby imply spurious indeterminacy. The danger
arises only in situations where a single equation involves some variables of
the form $E_tX(t+1)$ and other variables of the form $Y(t+1)$. In this
case, dropping the $E_t$ operators from the equation and adding an $\eta(t)$
error term loses the information that some variables are entering as actual
values and others as expectations.
The formulation that most writers on this subject have used, following
Blanchard and Kahn, is
\begin{equation}\label{eq:BKForm}
\Gamma_0E_ty(t+1)=\Gamma_1y(t)+C+\Psi z(t)\:.
\end{equation}
One then adds conditions that certain variables in the $y$ vector are
``predetermined", meaning that for them $E_ty(t+1)=y(t+1)$. The formulation
adopted here embodies the useful notational convention, that all variables
dated $t$ are observable at $t$; thus no separate list of what is
predetermined is needed to augment the information that can be read off from
the equations themselves. It also allows straightforward handling of the
commonly occurring systems in which a variable $y_i(t)$ and its expectation
$E_{t-1}y_i(t)$ both appear, with the same $t$ argument, in the same or
different equations. To get such systems into the Blanchard-Kahn form
requires introducing an extended state vector. For example, in our simple
overlapping contract model \eqref{eq:Taylor} above, the fact that $W$ appears
as an expected value in the first equation and/or in the dummy equation
defining the artificial state variable $\zeta(t)=E_tW(t+1)$ (the bottom row
of \eqref{eq:TaylorMats}), yet also in two other dynamic equations without any
expectation, seems to require that, to cast the model into standard
Blanchard-Quah form, we add $w(t-2)$ and $u(t-1)$ to the list of variables,
with two corresponding additional dummy equations.
Standard difference equations of the form \eqref{eq:mainEq} have a single, exogenous
disturbance vector. That is, they have $\Gamma_0=I$ and $\Pi=0$. They can therefore be
interpreted as determining $y(t)$ for $ t > t_0$ from given initial conditions
$ y(t_0 )$ and random draws for $ z(t)$. In \eqref{eq:mainEq}, however, the disturbance
vector $ \Psi z(t) + \Pi \eta (t)$ defined in \eqref{eq:mainEq} is not exogenously given
the way $z(t)$ itself is. Instead, it depends on $y(t)$ and its expectation,
both of which are generally unknown before we solve the model. Because we
need to determine $ \eta (t)$ from z(t) as we solve the model, we generally
need to find a number of additional equations or restrictions equal to the
rank of the matrix in order to obtain a solution.
\section{Canonical Forms and Matrix Decompositions}
Solving a system like \eqref{eq:mainEq} subject to restrictions on the rates
of growth of components of its solutions requires breaking its solutions
into components with distinct rates of growth. This is best done with
some version of an eigenvalue-eigenvector decomposition. We are already
imposing a somewhat stringent canonical form on the equation system by
insisting that it involve just one lag and just one-step-ahead expectations.
Of course as we made clear in the example above, systems with more lags or
with multi-step and lagged expectations can be transformed into systems of
the form given here, but there may be some computational work in making the
transformation. Further tradeoffs between system simplicity and simplicity
of the solution process are possible.
To illustrate this point, we begin by assuming a very stringent
canonical form for the system: $ \Gamma _0 = I$, $ E_t z\left( {t + 1}
\right) = 0$, all t. Systems derived from even moderately large rational
expectation equilibrium models often have singular $ \Gamma _0$ matrices, so
that simply ``multiplying through by $ \Gamma _0^{ - 1}$" to achieve this
canonical form is not possible. In most economic models, there is little
guidance available from theory in specifying properties for $z$. The
requirement that $ E_t z\left( {t + 1} \right) = 0$ is therefore extremely
restrictive.
On the other hand, it is usually possible, by solving for some variables in
terms of others and thereby reducing system size, to manipulate the system
into a form with non-singular $ \Gamma _0$. And it is common practice
to make a model tractable by assuming a simple flexible parametric form
for the process generating exogenous variables, so that the exogenous
variables themselves are incorporated into the $y$ vector, while the serially
uncorrelated $z$ vector in the canonical form is just the disturbance
vector in the process generating the exogenous variables. So this
initial very simple canonical form is in some sense not restrictive.
Nonetheless, getting the system into a form with non-singular
$ \Gamma _0$
may involve very substantial computational work if it is done
ad hoc, on an even moderately large system. And it is often
useful to display the dependence of $y$ on current, past, and
expected future exogenous variables directly, rather than to
make precise assumptions on how expected future $z$'s depend on
current and past $z$'s. For these reasons, we will below display
solution methods that work on more general canonical forms.
While these more general solution methods are themselves harder
to understand, they shift the burden of analysis from the
individual economist/model-solver toward the computer, and are
therefore useful.
\section{Using the Jordan decomposition with serially uncorrelated
shocks} \label{sect:Jordan}
\subsection{Discrete time}
In this section we consider the special case of
\begin{equation}
y(t) = \Gamma _1 y(t - 1) + C + \Psi z(t) + \Pi \eta (t)
\:,\end{equation}
with $E_t z\left( {t + 1} \right) = 0$. The system matrix $\Gamma _1$
has a Jordan decomposition
\begin{equation}\label{eq:G1Jord}
\Gamma _1 = P\Lambda P^{ - 1}\:,
\end{equation}
where $P$ is the matrix of right-eigenvectors of $\Gamma _1$, $P^{ - 1}$ is
the matrix of left-eigenvectors, and $\Lambda$ has the eigenvalues of $\Gamma
_1$ on its main diagonal and 0's everywhere else, except that it may have
$\lambda _{i,i + 1} = 1$ in positions where the corresponding diagonal
elements satisfy $\lambda _{ii} = \lambda _{i + 1,i + 1}$. Multiplying the
system on the left by $P^{ - 1}$, and defining $w = P^{ - 1} y$, we arrive at
\begin{equation}\label{eq:wForm}
w(t) = \Lambda w(t-1) + P^{ - 1} C + P^{ - 1}\cdot\left( {\Psi
z(t) + \Pi \eta (t)} \right) \:.\end{equation} In this setting, we can easily
consider non-homogeneous growth rate restrictions. That is, suppose that we
believe that a set of linear combinations of variables, $ \phi _i y(t)$,
$i=1,\ldots,m$, have bounded growth rates, with possibly different bounding growth
rates for each $i$. That is we believe that a solution must satisfy
\begin{equation}
E_s \left[ {\phi _i y(t)\xi _i^{ - t} } \right]\mathop \to
\limits_{t \to \infty } 0
\end{equation}
for each $i$ and $s$, with $\xi _i > 1$
for every $i$. Equation \eqref{eq:wForm} has isolated components of the system
that grow at distinct exponential rates. The matrix has a block-
diagonal structure, so that the system breaks into unrelated
components, with a typical block having the form
\begin{equation}\label{eq:wj}
w_j \left( t \right) = \begin{bmatrix}
{\lambda _j } & 1 & 0 & \cdots & 0 \\
0 & {\lambda _j } & 1 & \ddots & \vdots \\
0 & 0 & \ddots & \ddots & 0 \\
\vdots & \ddots & \ddots & {\lambda _j } & 1 \\
0 & \cdots & 0 & 0 & {\lambda _j } \\
\end{bmatrix}w_j \left( {t - 1} \right) + P^{j \cdot }
C + P^{j \cdot } \cdot \left( {\Psi z(t) + \Pi \eta (t)} \right)
\:,\end{equation} where $P^{j \cdot }$ is the block of rows of $P^{ - 1}$
corresponding to the $j$'th diagonal block of \eqref{eq:wForm}. If
the disturbance term (including the combined effects of $z$ and $\eta$ on the
equation) is zero and $\lambda _j \ne 1$, this block has the deterministic
steady-state solution
\begin{equation}\label{eq:wjSol}
w_j (t) = \left[ {I - \Lambda _j } \right]^{ - 1} P^{j \cdot }C
\:,\end{equation}
where $\Lambda _j$ is the Jordan block displayed in
\eqref{eq:wj}. If $\left| {\lambda _j } \right| > 1$, then $E_s \left[ {w_j
\left( {t + s} \right)} \right]$ grows in absolute value at the rate of
$\left| {\lambda _j } \right|^t$ as $t \to \infty$, for any solution other
than that given in \eqref{eq:wjSol}. Now consider the $k$'th restriction on growth,
\begin{multline}\label{eq:kthGrowth}
E_s \left[ {\phi _k y( t)} \right]\xi _k^{ - t} = \phi _k
PE_s \left[ {w(t)} \right]\xi _k^{ - t} \\
= \phi _k P\bigl[\Lambda^{t - s} ( {w(s) - \left( {I - \Lambda } \right)^{ - 1} P^{
- 1} C} )+(I-\Lambda)^{-1}P^{-1}C\bigr]\xi _k^{ - t} \to 0
\:.\end{multline}
In order for this condition to hold, every one of the vectors $w_j$
corresponding to a $\left| {\lambda _j } \right| > \xi _k$
and to a $\phi _k P_{ \cdot j} \ne 0$
must satisfy \eqref{eq:wjSol}. This is obvious for cases where $\Lambda _j$
is scalar. When $\Lambda _j$
is a matrix with ones on the first diagonal above the main
diagonal, a somewhat more elaborate argument is required. We
need to observe that we can expand terms of the form on the right
of \eqref{eq:kthGrowth} as
\begin{equation}\label{eq:JordSeries}
x\Lambda _j^s q = \lambda _j^s x \cdot \left( {q + s\lambda _j^{
- 1} q_{ - 1} + \frac{{s \cdot \left( {s - 1} \right)}}
{2}\lambda _j^{ - 2} q_{ - 2} + \ldots + c_n (s)\lambda _j^{ -
n + 1} q_{ - n + 1} } \right)
\:,\end{equation}
where $c_k (s)$
is the $(k, s)$'th binomial coefficient, i.e. $
{{s!} \mathord{\left/
{\vphantom {{s!} {\left( {\left( {s - k} \right)! \cdot k!}
\right)}}} \right.
\kern-\nulldelimiterspace} {\left( {\left( {s - k} \right)!
\cdot k!} \right)}}
$,
for $s \leqslant k$, and is 0 otherwise, and $q_{ - k}$
is the vector $q$ shifted up by $k$ with the bottom elements filled
out with zeros, i.e.
\begin{equation}
q_{ - k} = \begin{bmatrix}
{\mathop 0\limits_{(n - 1) \times 1} } &\vline & I \\
\hline
{\mathop 0\limits_{1 \times 1} } &\vline & {\mathop
0\limits_{1 \times (n - 1)} } \\
\end{bmatrix} \cdot q_{ - k + 1}
\:.\end{equation}
Using \eqref{eq:JordSeries}, it is straightforward, though still some work, to show
that indeed, for \eqref{eq:kthGrowth} to hold, every one of the vectors $w_j$
corresponding to a $\left| {\lambda _j } \right| > \xi _k$
and to a $\phi _k P_{ \cdot j} \ne 0$
must satisfy \eqref{eq:wjSol}.
Of course a problem must have special structure in order for it to turn out
that there is a $k, j$ pair such that $ \phi _k P_{ \cdot j} = 0$.
This is the justification for the (potentially misleading) common practice of
assuming that if any linear combination of $y$'s is constrained to grow slower
than $ \xi ^t$, then all roots exceeding $\xi_i$ in absolute value
must be suppressed in the solution. If \eqref{eq:wjSol} does hold for
all $t$ then we can see from \eqref{eq:wj} that this entails
\begin{equation}\label{eq:PjShock0}
P^{j \cdot } \cdot \left( {\Psi z + \Pi \eta } \right) = 0
\:.\end{equation} Collecting all the rows of $P^{ - 1}$ corresponding to $j$'s
for which \eqref{eq:PjShock0} holds into a single matrix $P^{U \cdot }$
(where the U stands for "unstable"), we can write
\begin{equation}\label{eq:PUShock0}
P^{U \cdot } \left( {\Psi z + \Pi \eta } \right) = 0
\:.\end{equation}
Existence problems arise if the endogenous shocks $\eta$ cannot adjust
to offset the exogenous shocks $z$ in \eqref{eq:PUShock0}. We might expect this to
happen if $P^{U \cdot }$
has more rows than has columns. This accounts for the usual
notion that there are existence problems if the number of
unstable roots exceeds the number of "jump variables". However,
the precise condition is that columns of $P^{U \cdot } \Pi$
span the space spanned by the columns of $P^{U \cdot } \Psi$, i.e.
\begin{equation}\label{eq:JordExist}
\Span\left( {P^{U \cdot } \Psi } \right) \subset
\Span\left( {P^{U \cdot } \Pi } \right) \:.\end{equation} In order
for the solution to be unique, it must be that \eqref{eq:PUShock0} pins down
not only the value of $ P^{U \cdot } \Pi \eta$, but also all the
other error terms in the system that are influenced by $\eta$. That is, from
knowledge of $ P^{U \cdot } \Pi \eta$ we must be able to determine $
P^{S \cdot } \Pi \eta$, where $ P^{S \cdot }$ is made up of all
the rows of $ P^{ - 1}$ not included in $ P^{U \cdot }$. Formally,
the solution is unique if and only if
\begin{equation}\label{eq:JordUnique}
\Span\left( {\Pi '\left( {P^{S \cdot } } \right)^\prime
} \right) \subset \Span\left( {\Pi '\left( {P^{U \cdot }
} \right)^\prime } \right)
\:.\end{equation}
In this case we will have
\begin{equation}
P^{S \cdot } \Pi \eta = \Phi P^{U \cdot } \Pi \eta
\end{equation}
for some matrix $\Phi$.
Usually we aim at writing the system in a form that can be simulated from
arbitrary initial conditions, delivering a solution path that does not
violate the stability conditions. We can construct such a system by
assembling the equations of the form delivered by the stability conditions
\eqref{eq:wjSol}, together with the lines of \eqref{eq:wj} that determine
$w_S$, the components of $w$ not determined by the stability
conditions, and use \eqref{eq:PUShock0} to eliminate dependence on $\eta$.
Specifically, we can use the system
\begin{equation}\label{eq:wSol}
\begin{bmatrix}
{w_S \left( t \right)} \\
{w_U \left( t \right)} \\
\end{bmatrix} = \begin{bmatrix}
{\Lambda _S } \\
0 \\
\end{bmatrix}w_S (t - 1) + \begin{bmatrix}
{P^{S \cdot } C} \\
{\left( {I - \Lambda _U } \right)^{ - 1} P^{U \cdot } C} \\
\end{bmatrix} + \begin{bmatrix}
I & { - \Phi } \\
0 & 0 \\
\end{bmatrix}P^{ - 1} \Psi z
\:.\end{equation}
To arrive at an equation in $y$, we use $y = Pw$
to transform \eqref{eq:wSol} into
\begin{multline}\label{eq:ySol}
y(t) = P_{ \cdot S} \Lambda _S P^{S \cdot } y(t - 1) + \left(
{P_{ \cdot S} P^{S \cdot }+ P_{ \cdot U} \left( {I -
\Lambda _U } \right)^{ - 1} P^{U \cdot } } \right)C \\ + \left(
{P_{ \cdot S} P^{S \cdot } - P_{ \cdot S} \Phi P^{U
\cdot } } \right)\Psi z\:.
\end{multline}
Labeling the three matrix coefficients in \eqref{eq:ySol}, we can give it the
form
\begin{equation}
y(t) = \Theta _1 y(t - 1) + \Theta _c C + \Theta _z z(t)
\:,\end{equation}
which can in turn be used to characterize the impulse responses
of $y$, according to
\begin{equation}\label{eq:yImpulses}
y(t + s) - E_t y(t + s) = \sum\limits_{v = 0}^{s - 1} {\Theta
_1^v \Theta _z z(t + s - v)}
\:.\end{equation}
Of course to completely characterize the mapping from initial
conditions and $z$ realizations to $y$, we need in addition to \eqref{eq:yImpulses} a
formula for the predictable part of $y$, i.e.
\begin{equation}\label{eq:yDet}
E_t y(t + s) = \Theta _1^s y(t) + (I - \Theta _1^{s + 1} ) \cdot \left( {I -
\Theta _1 } \right)^{ - 1} \Theta _C C \:.\end{equation} However all the
information needed to compute both \eqref{eq:yImpulses} and \eqref{eq:yDet}
is contained in a report of the coefficient matrices for \eqref{eq:ySol}.
Note that \eqref{eq:ySol}, while it insures that the second row of \eqref{eq:wSol},
\begin{equation}\label{eq:wUsol}
w_U (t) = P^{U \cdot } y(t) = \left( {I - \Lambda _U }
\right)^{ - 1} P^{U \cdot } C
\:,\end{equation}
holds for all $t$ after the initial date $t=0$, it does not in itself
impose \eqref{eq:wUsol} at $t=0$, which in fact is required by the solution.
\subsection{Continuous time}
In this type of canonical form, the analysis for continuous
time is nearly identical to that for discrete time. The equation
we start with is
\begin{equation}\label{eq:ctJordEq}
\dot y = \Gamma _1 y + C + \Psi z + \Pi \eta
\:,\end{equation}
where $z$ and $\eta$ are both assumed to be derivatives of martingale
processes, i.e. white noise. Just as in discrete time, we form a
Jordan decomposition of $\Gamma _1$
of the form \eqref{eq:G1Jord}. Again we use it to change variables to $w = P^{ - 1} y$
and split w into components $w_S$
and $w_U$
that need not be suppressed, and need to be suppressed,
respectively, to satisfy the stability conditions. The stability
conditions in continuous time are naturally written, analogously
to \eqref{eq:kthGrowth},
\begin{equation}
E_s \left[ {\phi _k y(t)} \right]e^{ - \xi _k t} = \phi _k PE_s \left[
{w(t)} \right]e^{ - \xi _k t}
\\= \phi _k Pe^{\Lambda \cdot \left( {t -
s} \right)} \left[ {w(s) - \Lambda ^{ - 1} P^{ - 1} C} \right]e^{ - \xi _k t}
\to 0 \:.
\end{equation}
The restricted components of $w$ are then those that
correspond both to a non-zero $\phi _k P^{ \cdot j}$ and to a $\lambda
_j$ with real part exceeding the corresponding $\xi _k$. Once we have
partitioned $w$, $P$, and $P^{ - 1}$ into $S$ and $U$ components according
to this criterion, the analysis proceeds as in the discrete case, with
the conditions for existence and uniqueness applying in unchanged form ---
\eqref{eq:JordExist} and \eqref{eq:JordUnique} . The final form of the
equation usable for simulation is
\begin{equation}
\dot y = P_{ \cdot S} \Lambda _S P^{S \cdot } y + P_{ \cdot
S} P^{S \cdot } C + P_{ \cdot S} \left( {P^{S \cdot } -
\Phi P^{U \cdot } } \right)\Psi z
\:.\end{equation}
This equation is usable to compute impulse responses, according
to
\begin{equation}
y\left( {t + s} \right) - E_t y(t + s) = \int\limits_0^s
{e^{\Theta _1 v} \Theta _z z\left( {t + s - v} \right)ds}
\end{equation}
and the deterministic part of $y$, according to
\begin{equation}
E_t y\left( {t + s} \right) = e^{\Theta _1 s} y(t) - \left( {I -
e^{\Theta _1 s} } \right)\Theta _1^{ - 1} \Theta _c C
\:.\end{equation}
In this case, in contrast to the discrete time case, the
preceding three equations contain no information about the
stability conditions restricting the levels of $y$ at all. We need
to specify separately the condition
\begin{equation}
w_U (t) = P^{U \cdot } y(t) = - \Lambda _U^{ - 1} P^{U \cdot
} C
\end{equation}
\section{Discrete time, solving forward}
In this section we consider the generic canonical form \eqref{eq:mainEq},
allowing for possibly singular $ \Gamma _0$
and not requiring $z$ to be serially uncorrelated. At first,
though, we consider only the case where there is a single bound
$ \bar \xi$
on the maximal growth rate of any component of $y$. We find
conditions that prevent such explosive growth as follows. First
we compute a QZ decomposition
\begin{equation}\label{eq:qz}
\begin{gathered}
Q'\Lambda Z' = \Gamma _0 \hfill \\
Q'\Omega Z' = \Gamma _1 \text{ }\text{.} \hfill \\
\end{gathered}
\end{equation}
In this decomposition, $Q'Q = Z'Z = I$, where $Q$ and $Z$ are both possibly
complex and the $\prime$ symbol indicates transposition and complex
conjugation. Also $\Omega$ and $\Lambda$ are possibly complex and are upper
triangular. The QZ decomposition always exists. Letting $w(t) = Z'y(t)$,
we can multiply \eqref{eq:mainEq} by $Q$ to obtain
\begin{equation}\label{eq:ctWform}
\Lambda w(t) = \Omega w(t - 1) + QC + Q\Pi \eta (t) + Q\Psi z(t)\:.
\end{equation}
Though the QZ decomposition is not unique, the collection of values for the
ratios of diagonal elements of $\Omega$ and $\Lambda$, $ \set{ \omega _{ii}
/\lambda _{ii} }$ (called the set of generalized eigenvalues), is usually
unique (if we include $\infty$ as a possible value). The generalized
eigenvalues are indeterminate only when $ \Gamma _0$ and $ \Gamma _1$ have
zero eigenvalues corresponding to the same eigenvector.\footnote{This would
imply that a linear combination of the equations contains no y's, i.e. that
there is effectively one equation fewer for determining the y's than would
appear from the order of the system.} We can always arrange to have the
largest of the generalized eigenvalues in absolute value appear at the
lower right. In particular, let us suppose that we have partitioned \eqref{eq:wForm} so
that $ \left| {{{\omega _{ii} } \mathord{\left/ {\vphantom {{\omega _{ii} }
{\lambda _{ii} }}} \right. \kern-\nulldelimiterspace} {\lambda _{ii} }}}
\right| \geqslant \bar \xi$
for all i>k and $\left| {{{\omega _{ii} } \mathord{\left/
{\vphantom {{\omega _{ii} } {\lambda _{ii} }}} \right.
\kern-\nulldelimiterspace} {\lambda _{ii} }}} \right| < \bar \xi$
for all $i\le k$. Then \eqref{eq:wForm} can be expanded as
\begin{equation}\label{eq:wFormPart}
\begin{gathered}
\begin{bmatrix}
{\Lambda _{11} } & {\Lambda _{12} } \\
0 & {\Lambda _{22} } \\
\end{bmatrix} \cdot \begin{bmatrix}
{w_1 (t)} \\
{w_2 (t)} \\
\end{bmatrix} = \begin{bmatrix}
{\Omega _{11} } & {\Omega _{12} } \\
0 & {\Omega _{22} } \\
\end{bmatrix} \cdot \begin{bmatrix}
{w_1 (t - 1)} \\
{w_2 (t - 1)} \\
\end{bmatrix} \hfill \\
\text{ } +
\begin{bmatrix}
{Q_{1 \cdot } } \\
{Q_{2 \cdot } } \\
\end{bmatrix}\left( {C + \Psi z(t) + \Pi \eta (t)}
\right) \hfill \\
\end{gathered}
\end{equation}
Note that some diagonal elements of $ \Lambda _{22}$, but not of $ \Lambda
_{11}$, may be zero. Also note that zeros at the same position $i$ on the
diagonals of both $\Lambda$ and $\Omega$ cannot occur unless the equation system is
incomplete, meaning that some equation is exactly a linear combination of
the others.
Because of the way we have grouped the generalized eigenvalues,
the lower block of equations in \eqref{eq:ctWform} is purely
explosive. It has a solution that does not explode any faster than the
disturbances $z$ so long as we solve it "forward" to make $ w_2$ a function
of future z's. That is, if we label the last additive term in
\eqref{eq:wFormPart} $x(t)$ and set $M=\Omega _{22}^{ - 1}\cdot\Lambda_{22}$,
\begin{multline}\label{eq:forward}
Z'_{ \cdot 2} y(t) = w_2 (t) = Mw_2 (t + 1) - \Omega
_{22}^{ - 1} x_2 (t + 1) \\
= M^2 \cdot w_2 (t + 2) - M \cdot \Omega _{22}^{ - 1}
\cdot x_2 (t + 2) - \Omega _{22}^{ - 1} \cdot x_2 (t + 1) \\
= - \sum\limits_{s = 1}^\infty {M^{s - 1} \cdot \Omega
_{22}^{ - 1} \cdot x_2 (t + s)}\:.
\end{multline}
The last equality in \eqref{eq:forward} follows on the assumption that $M^t
w_2 (t) \to 0$ as $t\to\infty$. Note that in the special case of $\lambda
_{ii} = 0$ there are equations in \eqref{eq:wFormPart} containing no
current values of $w$. While these cases do not imply explosive
growth, the corresponding components of \eqref{eq:forward} are still valid.
For example, if the lower right element of $\Lambda$ is zero, the last equation of
\eqref{eq:wFormPart} has the form
\begin{equation}
0 \cdot w_n (t) = \omega _{nn} \cdot w_n (t - 1) + x_n (t).
\end{equation}
Solving for $w_n (t - 1)$ produces the corresponding component of
\eqref{eq:forward}. Since $\lambda _{ii} = 0$ corresponds to a singularity
in $\Gamma _0$, the method we are describing handles such singularities
transparently.
Note that \eqref{eq:forward} asserts the equality of something
on the left that is known at time $t$ to something on the right that
is a combination of variables dated $t+1$ and later. Since taking
expectations as of date $t$ leaves the left-hand side unchanged, we can write
\begin{multline}\label{eq:decision}
Z'_{ \cdot 2} y(t) = w_2 (t) = - E_t \left[ {\sum\limits_{s
= 1}^\infty {M^{s - 1} \cdot \Omega _{22}^{ - 1} \cdot x_2 (t
+ s)} } \right]\\
= - \sum\limits_{s = 1}^\infty {M^{s -
1} \cdot \Omega _{22}^{ - 1} \cdot x_2 (t + s)} \,\text{.}
\end{multline}
If the original system \eqref{eq:mainEq} consisted of first order
conditions from an optimization problem, \eqref{eq:decision} will be what is
usually called the decision rule for the problem. When the original
system described a dynamic market equilibrium, \eqref{eq:decision} will be composed of
decision rules of the various types of agents in the economy, together
with pricing functions that map the economy's state into a vector of prices.
The last equality in \eqref{eq:decision} imposes conditions on $x_2$
that may or may not be consistent with the economic
interpretation of the model. Recall that $ x_2$
is made up of constants, terms in $z$, and terms in $\eta$. But $z$ is
an exogenously given stochastic processes whose properties cannot
be taken to be related to $ \Gamma _0$
and $ \Gamma _1$. Thus if it should turn out that $ x_2$
contains no $\eta$ component, the assertion in \eqref{eq:decision} will be impossible
--- it requires that exogenously evolving events that occur in the
future be known precisely now. If $ x_2$
does contain an $\eta$ component, then \eqref{eq:decision} asserts that this
endogenously determined component of randomness must fluctuate as
a function of future $z$'s so as exactly to prevent any deviation
of the right-hand side of \eqref{eq:decision} from its expected value.
Replacing $x$'s in \eqref{eq:decision} with their definitions, that equation
becomes
\begin{multline}
Z'_{ \cdot 2} y(t) = \left( {\Lambda _{22} - \Omega _{22}
} \right)^{ - 1} Q_{2 \cdot } C - E_t \left[ {\sum\limits_{s =
1}^\infty {M^{s - 1} \Omega _{22}^{ - 1} Q_{2 \cdot } \Psi z(t
+ s)} } \right] \\
= \left( {\Lambda _{22} - \Omega _{22} } \right)^{ - 1}
Q_{2 \cdot } C - \sum\limits_{s = 1}^\infty {M^{s - 1} \Omega
_{22}^{ - 1} Q_{2 \cdot } \bigl( {\Psi z(t + s) + \Pi \eta (t +
s)} \bigr)}
\:.\end{multline}
The latter equality is satisfied if and only if
\begin{equation}\label{eq:etaz}
Q_{2 \cdot } \Pi \eta (t + 1) = \sum\limits_{s = 1}^\infty
{\Omega _{22} M^{s - 1} \Omega _{22}^{ - 1} Q_{2 \cdot } \Psi
\cdot \left( {E_{t + 1} z(t + s) - E_t z(t + s)} \right)} \:.
\end{equation}
A leading special case is that of serially uncorrelated $z$'s, i.e. $E_t z(t +
s) = 0$, all $s>1$. In this case the right-hand-side of \eqref{eq:etaz} is
just $Q_{2\cdot}\Psi z(t+1)$, so a necessary condition for the existence of
a solution\footnote{Note that here it is important to our analysis that
there are no hidden restrictions on variation in $z$ that cannot be
deduced from the structure of the equation system. In an equation
system in which there are two exogenous variables with $z_1 (t) = 2z_2 (t -
2)$, for example, our analysis requires that this restriction connecting
the two exogenous variables be included as an equation in the system and the
number of exogenous variables be reduced to one.} satisfying \eqref{eq:etaz}
is that the column space of $Q_{2 \cdot } \Psi$ be contained in that of $Q_{2
\cdot } \Pi$, i.e.
\begin{equation}\label{eq:qzNscExist}
\Span\left( {Q_{2 \cdot } \Psi } \right) \subset \Span\left( {Q_{2
\cdot } \Pi } \right) \:.
\end{equation}
This condition is necessary and sufficient when $E_t z(t + 1) = 0$, all $t$.
A necessary and sufficient condition for a solution to exist regardless of the
pattern of changes in the expected future time path of $z$'s is\footnote{It
may appear that $s$ should range from $1$ to infinity, rather than $1$ to
$n-k$, in this expression. But $\Omega _{22}^{ - 1} Q_{2 \cdot } \Psi$ is
in some invariant subspace of $M$, if only the trivial full $n-k$-dimensional
space Euclidean space. The invariant space containing it, say of
dimension $j$, will be spanned by $j$ elements of the $M^{s - 1} \Omega _{22}^{ -
1} Q_{2 \cdot } \Psi$ sequence.}
\begin{equation}\label{eq:qzScExist}
\Span\left( {\left\{ {\Omega _{22} M^{s - 1} \Omega
_{22}^{ - 1} Q_{2 \cdot } \Psi } \right\}_{s = 1}^{n - k} }
\right) \subset \Span\left( {Q_{2 \cdot } \Pi } \right)
\:.
\end{equation}
In most economic models this latter necessary and sufficient
condition is the more meaningful one, even if it is always true
that $E_t z(t + 1) = 0$, because ordinarily our theory places no reliable restrictions
on the serial dependence of the $z$ process, even if we have made
some assumption on it for the purpose at hand.
Assuming a solution exists, we can combine \eqref{eq:etaz}, or its equivalent
in terms of $w$ \eqref{eq:forward}, with some linear combination of equations
in \eqref{eq:wFormPart} to obtain a new complete system in $w$ that is stable.
However, the resulting system will not be directly useful for
generating simulations or distributions of $y$ from specified processes for
$z$ unless we can free it from references to the endogenous error term $\eta$. From
\eqref{eq:etaz}, we have an expression that will determine $ Q_{2 \cdot } \Pi \eta
(t)$ from information available at $t$ and a known stochastic process for
$z$. However the system also involves a different linear transformation of
$\eta$, $ Q_{1 \cdot } \Pi \eta (t)$. It is possible that knowing $ Q_{2
\cdot } \Pi \eta (t)$ is not enough to tell us the value of $ Q_{1
\cdot } \Pi \eta (t)$, in which case the solution to the model is not
unique. In order that the solution be unique it is necessary and sufficient
that the row space of $ Q_{1 \cdot } \Pi$ be contained in that of $
Q_{2 \cdot } \Pi$. In that case we can write
\begin{equation}\label{eq:etaLink}
Q_{1 \cdot } \Pi = \Phi Q_{2 \cdot } \Pi
\end{equation}
for some matrix $\Phi$. Premultiplying \eqref{eq:wFormPart} by $[I \quad - \Phi]$
gives us a new set of equations, free of references to $\eta$, that
can be combined with \eqref{eq:forward} to give us
\begin{multline}\label{eq:wSolution}
\begin{bmatrix}
{\Lambda _{11} } & {\Lambda _{12} - \Phi \Lambda _{22} }
\\
0 & I \\
\end{bmatrix} \cdot \begin{bmatrix}
{w_1 (t)} \\
{w_2 (t)} \\
\end{bmatrix} \\= \begin{bmatrix}
{\Omega _{11} } & {\Omega _{12} - \Phi \Omega _{22} } \\
0 & 0 \\
\end{bmatrix} \cdot \begin{bmatrix}
{w_1 (t - 1)} \\
{w_2 (t - 1)} \\
\end{bmatrix}
+ \begin{bmatrix}
{Q_{1 \cdot } - \Phi Q_{2 \cdot } } \\
{(\Omega _{22} - \Lambda _{22} )^{ - 1} Q_{2 \cdot } }
\end{bmatrix}C \\
+ \begin{bmatrix}
{Q_{1 \cdot } - \Phi Q_{2 \cdot } } \\
0 \\
\end{bmatrix}\Psi z(t) - E_t \begin{bmatrix}
0 \\
{\sum\nolimits_{s = 1}^\infty {M^{s - 1} \Omega _{22}^{ -
1} Q_{2 \cdot } \Psi z(t + s)} } \\
\end{bmatrix} \:.
\end{multline}
This can be translated into a system in $y$ of the form
\begin{equation}\label{eq:ySolution}
y(t) = \Theta _1 y(t - 1) + \Theta _c + \Theta _0 z(t) +
\Theta _y \sum\limits_{s = 1}^\infty {\Theta _f^{s - 1} \Theta
_z E_t z(t + s)}
\:.\end{equation}
The details of the translation are
\begin{equation}\label{eq:translate}
\begin{gathered}
H = Z\begin{bmatrix}
{\Lambda _{11}^{ - 1} } & { - \Lambda _{11}^{ - 1} \left(
{\Lambda _{12} - \Phi \Lambda _{22} } \right)} \\
0 & I \\
\end{bmatrix}\:;\quad\Theta _1 = Z_{ \cdot 1}
\Lambda _{11}^{ - 1} \begin{bmatrix}
{\Omega _{11} } & {(\Omega _{12} - \Phi \Omega _{22} )}
\\
\end{bmatrix}Z\:; \\
\Theta _c = H \cdot \begin{bmatrix}
{Q_{1 \cdot } - \Phi Q_{2 \cdot } } \\
{\left( {\Omega _{22} - \Lambda _{22} } \right)^{ - 1}
Q_{2 \cdot } }
\end{bmatrix}C\:;\quad \Theta _0 = H \cdot
\begin{bmatrix}
{Q_{1 \cdot } - \Phi Q_{2 \cdot } } \\
0 \\
\end{bmatrix} \cdot \Psi\:; \\
\Theta _y = - H_{ \cdot 2} \:;\quad\Theta _f =
M\:;\quad\Theta _z = \Omega _{22}^{ - 1} Q_{2 \cdot } \Psi \:.
\end{gathered}
\end{equation}
The system defined by \eqref{eq:ySolution} and \eqref{eq:translate} can
always be computed and is always a complete equation system for $y$
satisfying the condition that its solution grow slower than $ \bar \xi ^t$,
even if there is no solution for $\eta$ in \eqref{eq:etaz} or the solution for
$Q_{1\cdot}\Pi$ in \eqref{eq:etaLink} is not unique. If there is no solution to
\eqref{eq:etaz}, then \eqref{eq:ySolution}-\eqref{eq:translate} implicitly restricts the
way $z$ enters the system, achieving stability by contradicting the original specification.
If the solution to \eqref{eq:etaLink} is not unique, then the absence
of $\eta$ from \eqref{eq:ySolution}-\eqref{eq:translate} contradicts the
original specification. If the solution is not unique, but $\Phi$ is
computed to satisfy \eqref{eq:etaLink}, the
\eqref{eq:ySolution}-\eqref{eq:translate} system generates one of the multiple
solutions to \eqref{eq:mainEq} that grows slower that $ \bar \xi ^t$. If
one is interested in generating the full set of non-unique solutions, one has
to add back in, as additional ``disturbances", the components of
$Q_{1\cdot}\Pi\eta$ left undetermined by \eqref{eq:etaLink}.
To summarize, we have the following necessary and sufficient
conditions for existence and uniqueness of solutions satisyfing
the bounded rate of growth condition:
\begin{enumerate}\renewcommand{\theenumi}{\Alph{enumi}}
\item \label{item:existGen} A necessary and sufficient condition that \eqref{eq:mainEq} has a
solution meeting the growth condition for arbitrary
assumptions on the expectations of future $z$'s is that the
column space spanned by
\[
\left\{ {\Omega _{22} M^{s - 1} \Omega _{22}^{ - 1} Q_{2
\cdot } \Psi } \right\}_{s = 1}^{n - k}
\]
be contained in that of $ Q_{2 \cdot } \Pi$.
\item \label{item:uniqueGen} A necessary and sufficient condition that any solution to
\eqref{eq:mainEq} be unique is that the row space of $ Q_{1 \cdot } \Pi$
be contained in that of $ Q_{2 \cdot } \Pi$.
\end{enumerate}
Condition \ref{item:existGen} takes a simpler form if the system has fully
specified the dynamics of exogenous variables:
\begin{list}{\ref{item:existGen}$'$.}{}
\item A necessary and sufficient condition that \eqref{eq:mainEq} has a solution
meeting the growth condition for arbitrary assumptions on the
covariance matrix of serially uncorrelated $z$'s is \eqref{eq:qzNscExist}, i.e. that
the column space spanned by $ Q_{2 \cdot } \Psi$
be contained in that of $ Q_{2 \cdot } \Pi$.
\end{list}
When condition A is met, a solution is defined by
\eqref{eq:ySolution}-\eqref{eq:translate}. In the special case of $ E_t
z\left( {t + 1} \right) = 0$, the last term of \eqref{eq:ySolution} (involving $ \Theta _y$, $
\Theta _f$ and $ \Theta _z$ ) drops out.
\section{Computational details}
If $A$ has full column rank, we can check to see whether the column
space of a matrix $A$ includes that of a matrix $B$ by ``regressing" $B$
on $A$ to see if the residuals are zero, i.e. by checking
\begin{equation}
\left( {I - A(A'A)^{ - 1} A'} \right)B = 0
\:.\end{equation}
If $A$ has full row rank, then its column space automatically
includes any other space of the same dimension. If A has neither
full row nor full column rank, other methods are required. The
singular value decomposition (svd) of a matrix $A$ is a
representation
\begin{equation}
A = UDV'
\end{equation}
in which $U$ and $V$ satisfy $U'U=I=V'V$ (but are not in general
square) and $D$ is square and diagonal.\footnote{This is not actually the
standard version of the svd. Matlab returns with $U$ and $V$ square, $D$ of the
same order as $A$. From such an svd, the form discussed in the text can be
obtained by replacing $D$ with a square non-singular matrix that has the non-
zero diagonal elements of the original $D$ on the diagonal and by forming
$U$ and $V$ from the columns of the original $U$ and $V$ corresponding to
non-zeros on the diagonal of the original $D$.} If $B$'s svd is $TCW$, $A$'s column
space includes $B$'s if and only if
\begin{equation}\label{eq:spanCondition}
\left( {I - UU'} \right)T = 0
\:.\end{equation}
If \eqref{eq:spanCondition} holds, then
\begin{equation}\label{eq:spanCalc}
B = A \cdot VD^{ - 1} U'B \:.\end{equation} Equation \eqref{eq:spanCondition}
gives us a computationally stable way to check the column and row space
spanning conditions summarized in \ref{item:existGen} and \ref{item:uniqueGen}
at the end of the previous section, and \eqref{eq:spanCalc} is a
computationally stable way to compute the matrix transforming $A$ to $B$, and thus
to compute the $\Phi$ matrix in \eqref{eq:etaLink}.\footnote{Of course since $\Phi$ is applied
to rows rather than columns, corresponding adjustments have to be made in the formula.}
Though the QZ decomposition is available in Matlab, it emerges
with the generalized eigenvalues not sorted along the diagonals
of $\Lambda$ and $\Omega$. Since our application of it depends on getting the
unstable roots into the lower right corner, an auxiliary routine
is needed to sort the roots around a $ \bar \xi$
level.\footnote{The routine \texttt{qzdiv} does this and is available, with other
Matlab routines that implement the methods of this paper, at
\texttt{http://www.princeton.edu/$\sim$sims/\#gensys}.}
If $ \Gamma _0$ in \eqref{eq:mainEq} is invertible we can multiply through
by its inverse to obtain a system which, like that in section \ref{sect:Jordan}, has $ \Gamma
_0 = I$. In such a system the QZ decomposition delivers $ Q = Z'$, and the
decomposition $ \Gamma _1 = Q'\Omega Q$ is what is known as the (complex)
Schur decomposition of the matrix $ \Gamma _1$. Since the Schur
decomposition is somewhat more widely available than the QZ, this may be
useful to know. Also in such a system we may be able to use an
ordinary eigenvalue decomposition of $ \Gamma _1$ in place of either the
Schur or the QZ. Most such routines return a matrix $V$ whose columns are
eigenvectors of $ \Gamma _1$, together with a vector $\lambda$ of eigenvalues. If
$V$ turns out to be non-singular, as it will if $ \Gamma _1$ has no repeated
roots, then $ \Gamma _1 = V\diag(\lambda) V^{ - 1}$, and this decomposition can
be used to check existence and uniqueness and to find a system of the form \eqref{eq:ySolution}
that generates the stable solutions.
With $V$, $V^{ - 1}$ and $\lambda$ partitioned to put the excessively
explosive roots in the lower right, we can write the conditions for
existence and uniqueness just as in \ref{item:existGen} and
\ref{item:uniqueGen} of the previous section but setting the matrices in those
conditions to the following:
\begin{equation}\label{eq:om2V}
\Omega _{22} = \diag(\lambda _{22})\:;\quad M = \diag(\lambda _{22}^{ - 1})\:;\quad
Q_{2 \cdot } = V^{2 \cdot }\:;\quad
Q_{1 \cdot } = V^{1 \cdot }
\:.
\end{equation}
Here $V^{i \cdot }$ refers to the i'th block of rows of $V^{ - 1}$. The
calculations in \eqref{eq:translate} can be carried out based on
\eqref{eq:om2V} also, with $\Lambda=I$. The advantage of this approach to
the problem is that, even though \eqref{eq:om2V} is written with $V$, $V^{ - 1}$ and
$\lambda$ ordered in a particular way, if the roots to be re-ordered are
distinct, they can be re-ordered simply by permuting the elements of $\lambda$, the
columns of $V$, and the rows of $V^{ - 1}$. There is no need here, as there
is with QZ, to recompute elements of the matrices when the decomposition is
re-ordered.
An intermediate form of simplification is available when $ \Gamma _0$
is singular, but the generalized eigenvalues are all
distinct. In that case it is possible to diagonalize $\Lambda$ and $\Omega$
in \eqref{eq:qz} to arrive at
\begin{equation}\label{eq:GenJord}
\begin{gathered}
P\Lambda _0 R = \Gamma _0 \hfill \\
P\Lambda _1 R = \Gamma _1 \hfill \\
\end{gathered}
\end{equation}
in which the $\Lambda _i$
are diagonal and $P$ and $R$ are both nonsingular. Here the
conditions, which we omit in detail, are almost exactly as in the
QZ case, but as with \eqref{eq:om2V}, any re-ordering of the decomposition that
may be required can be done by permutations rather than requiring
new calculations.
\section{More general growth rate conditions}
Properly formulated dynamic models with multiple state variables generally do
not put a uniform growth rate restriction on all the equations of the system. Transversality
conditions, for example, usually restrict the ratio of a wealth variable to
marginal utility not to explode. In a model with several assets, there is no
constraint that individual assets not exponentially grow or shrink, so long as
they do so in such a way as to leave total wealth satisfying the
transversality condition. Ignoring this point can lead to missing
sources of indeterminacy.
No program that simply calculates roots and
counts how many fall in various partitions will be able to distinguish such
cases. To handle conditions like these formally, we proceed as follows.
We suppose that boundary conditions at infinity are given in the form
\begin{equation}\label{eq:GenGrowthLim}
\xi _i^{ - t} H_i y(t)\mathop \to \limits_{t \to \infty } 0 \:,
\end{equation}
$i=1,\ldots,p$. Here $H_i y$ is the set of linear combinations of y that are constrained to
grow slower than $\xi _i^t$. Suppose we have constructed the QZ
decomposition \eqref{eq:qz} for our system and the j'th generalized eigenvalue $
{{\omega _{jj} } \mathord{\left/ {\vphantom {{\omega _{jj} } {\lambda _{jj}
}}} \right. \kern-\nulldelimiterspace} {\lambda _{jj} }}
$
exceeds $\xi _i$
in absolute value. To see whether this root needs to be put in
the forward part of the solution or instead belongs in the
backward part -- i.e. to see whether the boundary condition
generates a restriction -- we must re-order the QZ decomposition
so that the j'th root appears at the upper left. Then we can
observe that
\begin{equation}
H_i y = H_i ZZ'y = H_i Zw \:,
\end{equation}
where $w = Z'y$ as in \eqref{eq:wForm}. Assuming that no component of $w$
other than the first produces explosive growth in $H_i y$ faster than $\xi
_i^t$, the first component produces such growth if and only if the first
column of $H_i Z$ is non-zero. If this first column is non-zero, the j'th
root generates a restriction, otherwise it does not. A test of this sort,
moving the root in question to the upper left of the QZ decomposition,
needs to be done for each root that exceeds in absolute value any of the
$\xi _i$'s. Roots can be tested this way in groups, with a block of
roots moved to the upper left together, and if it turns out that the block
generates no restrictions, each one of the roots generates no
restrictions. However if the block of roots generates restrictions,
each root must then be tested separately to see whether it generates
restrictions by itself. The only exception is complex pairs of roots,
which in a real system should generate restrictions or not, jointly.
When we have this type of boundary condition, there may be a particular
appeal to trying to achieve the decomposition \eqref{eq:GenJord}, as with
that decomposition the required re-orderings become trivial. Note that all
that is essential to avoiding the work of re-ordering is that the
generalized eigenvalues that are candidates for violating boundary conditions
all be distinct and distinct from the eigenvalues that are not
candidates for violating boundary conditions. In that case we can
achieve a block diagonal version of \eqref{eq:GenJord}, with all the roots
not candidates for violating boundary conditions grouped in one block and
the other roots entering separately on the diagonal.
The Matlab programs \texttt{gensys} and \texttt{gensysct} that accompany this
paper do not automatically check the more general boundary conditions
discussed in this section of the paper. Each has a component, however,
called respectively \texttt{rawsys} and \texttt{rawsysct}, that computes the forward and
backward parts of the solution from a QZ decomposition that has been sorted
with all roots that generate restrictions grouped at the lower right. There
is also a routine available, called qzswitch, that allows the user to re-order
the QZ in any desired way. With these tools, it is possible to check
existence and uniqueness and find the form of the solution when boundary
conditions take the form \eqref{eq:GenGrowthLim}.
\section{Extension to a general continuous time system}
Consider a system in continuous time\footnote{In this section we assume the
reader is familiar with the definition of a continuous time white noise and
understands how integrals of them over time can be used to represent
continuous time stochastic processes.} of the form
\begin{equation}\label{eq:contGenEq}
\Gamma _0 \dot y = \Gamma _1 y + C + \Psi z + \Pi \eta \:,
\end{equation}
with an endogenous white noise error $\eta$ and $z$ an exogenous process that may
include a white noise component with arbitrary covariance matrix. By a
"white noise" here we mean the time-derivative of a martingale. The
martingales corresponding to $z$ and could have jump discontinuities in
their paths without raising problems for this paper's analysis. A fully
general analysis of continuous-time models is messier than for discrete-time
models, because when $z$ is an arbitrary exogenous process, the solution may
in general contain both a ``forward component" like that in the discrete
time case and a ``differential component" that relates $y$ to the
non-white-noise components of first and higher- order derivatives of $z$. We
therefore first work through the case of white-noise $z$, which is only a
slight variation on the analysis for the discrete-time model with
serially uncorrelated $z$.
An example of such a system is the continuous time analogue of
\eqref{eq:Taylor},
\begin{equation}
\begin{aligned}
w(t) &= .3E_t \int\limits_{s = 0}^\infty {e^{ - .3s} W(t +
s)ds} - \alpha \cdot (u(t) - u_n ) + \nu (t) \\
W(t) &= .3\int\limits_{s = 0}^\infty {e^{ - .3s} w(t - s)ds}
\\
\dot u(t) &= - \theta \cdot u(t) + \gamma \cdot W(t) + \mu
+ \varepsilon (t)\text{ ,} \\
\end{aligned}
\end{equation}
which can be rewritten as
\begin{equation}\label{eq:ctTaylor}
\begin{aligned}
\dot w &= .3w - .3W - \alpha \dot u + .3\alpha \cdot (u -
u_n ) + z_1 - .3\nu + \eta _1 \\
\dot \nu &= z_1 \\
\dot W &= - .3W + .3w \\
\dot u &= - \theta u + \gamma W + \mu + z_2 \:. \\
\end{aligned}
\end{equation}
Note that to fit the framework with all exogenous disturbances
white noise, we have made $\nu$ a martingale. We could make other
assumptions instead on the process generating $\nu$, but to stay in
this simple framework we need an explicit form for the process.
In the notation of \eqref{eq:contGenEq}, \eqref{eq:ctTaylor} has
\begin{gather*}
\Gamma _0 = \begin{bmatrix}
1 & 0 & 0 & \alpha \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{bmatrix}\text{ , }\Gamma _1 = \begin{bmatrix}
{.3} & { - .3} & { - .3} & {.3\alpha } \\
0 & 0 & 0 & 0 \\
{.3} & 0 & { - .3} & 0 \\
0 & 0 & \gamma & { - \theta } \\
\end{bmatrix}\text{ ,}\\
C = \begin{bmatrix}
{ - .3\alpha u_n } \\
0 \\
0 \\
\mu \\
\end{bmatrix}\text{ , }\Psi = \begin{bmatrix}
1 & 0 \\
1 & 0 \\
0 & 0 \\
0 & 1 \\
\end{bmatrix}\text{ , }\Pi \text{ = }\begin{bmatrix}
\text{1} \\
\text{0} \\
\text{0} \\
\text{0} \\
\end{bmatrix}
\:.\end{gather*}
We can proceed as before to do a QZ decomposition of $ \Gamma _0$
and $ \Gamma _1$
using the notation of \eqref{eq:qz}, arriving at the analogue of \eqref{eq:wForm}
\begin{equation}\label{eq:ctWformGen}
\Lambda \dot w = \Omega w + QC + Q\Pi \eta + Q\Psi z \:.
\end{equation}
Again we want to partition the system to arrive at an analogue to
\eqref{eq:wFormPart}, but now instead of putting the roots largest in
absolute value in the lower right corner of $\Omega$, we want to put there
the roots with algebraically largest (most positive) real parts.
Cases where $\Lambda$ has zeros on the diagonal present somewhat
different problems here than in the discrete case. In the
discrete case, we can think of the ratio of a non-zero $ \Omega _{ii}$
to a zero $ \Lambda _{ii}$
as infinite, which is surely very large in absolute value,
and then treat these pairs as corresponding to explosive roots.
Here, on the other hand, we are not concerned with the absolute
values of roots, only with their real parts. The ratio of a non-
zero value to zero cannot be treated as having a well-defined
sign on its real part. Nonetheless, it turns out that we want to
treat zeros on the diagonal of $\Lambda$ as producing ``unstable" roots.
If $ \lambda _{ii}$
is close to zero and $
{{\omega _{ii} } \mathord{\left/
{\vphantom {{\omega _{ii} } {\lambda _{ii} }}} \right.
\kern-\nulldelimiterspace} {\lambda _{ii} }}
$
has positive real part, there is no doubt that we classify
the corresponding generalized eigenvalue as unstable. If $
{{\omega _{ii} } \mathord{\left/
{\vphantom {{\omega _{ii} } {\lambda _{ii} }}} \right.
\kern-\nulldelimiterspace} {\lambda _{ii} }}
$
has negative, extremely large real part, or a negative real part of any size
combined with an extremely large imaginary part, we may still be justified
in treating it as an unstable root. This is true as a matter of numerical
analysis because extremely small $ \lambda _{ii}$ values may differ from
zero only by accumulated rounding error. But also as a matter of
economic interpretation, extremely large generalized eigenvalues correspond
to components of $y$ that, though they are technically defined as
random variables at each date, behave so erratically that they approach
white-noise like behavior.
We now proceed with our plan to partition \eqref{eq:ctWformGen} to obtain
an analogue to \eqref{eq:wFormPart}. In the lower right we place all cases
of $ \lambda _{ii}$ extremely close to zero, then just above that all cases
of $ {{\omega _{ii} } \mathord{\left/ {\vphantom {{\omega _{ii} } {\lambda
_{ii} }}} \right. \kern-\nulldelimiterspace} {\lambda _{ii} }}$ positive and
exceeding some boundary level $ \bar \xi \geqslant 0$. The resulting
system is
\begin{equation}\label{eq:ctWformPartGen}
\begin{bmatrix}
{\Lambda _{11} } & {\Lambda _{12} } & {\Lambda _{13} } \\
0 & {\Lambda _{22} } & {\Lambda _{23} } \\
0 & 0 & {\Lambda _{33} } \\
\end{bmatrix} \cdot \begin{bmatrix}
{\dot w_1 } \\
{\dot w_2 } \\
{\dot w_3 } \\
\end{bmatrix} = \begin{bmatrix}
{\Omega _{11} } & {\Omega _{12} } & {\Omega _{13} } \\
0 & {\Omega _{22} } & {\Omega _{23} } \\
0 & 0 & {\Omega _{33} } \\
\end{bmatrix} \cdot \begin{bmatrix}
{w_1 } \\
{w_2 } \\
{w_3 } \\
\end{bmatrix} + \begin{bmatrix}
{Q_{1 \cdot } } \\
{Q_{2 \cdot } } \\
{Q_{3 \cdot } } \\
\end{bmatrix}\left( {C + \Psi z + \Pi \eta } \right)
\:.\end{equation}
The last block of equations in \eqref{eq:ctWformPartGen} can be written as
\begin{equation}\label{eq:w3ct}
w_3 = \Omega _{33}^{ - 1}\bigl( \Lambda _{33} \dot w_3 - Q_{3 \cdot } (
{C + \Psi z + \Pi \eta } )\bigr) \:.
\end{equation}
Because we require that $w_3$ be a
random variable observable at $t$, $E_t w_3 \left( t \right) = w_3 \left( t
\right)$. For white-noise $z$ and $\eta$, $E_t z(t) = E_t \eta \left( t \right) = 0$.
Since $\Lambda _{33}$ is upper triangular with zeros on the diagonal, its
bottom row is zero. Thus we can see from \eqref{eq:w3ct} that the bottom element in the $w$
vector, $w_{3n}$, satisfies $w_{3n} = - \omega_{nn}^{-1}Q_{3n \cdot } C$, where
$\omega_{nn}$ is the lower right diagonal element of $\Omega$. But now
we can proceed recursively up to the second-to-last row of \eqref{eq:w3ct}, etc. to conclude
that in fact
\begin{equation}
\begin{gathered}
w_3 = - \Omega_{33}^{ - 1} Q_{3 \cdot } C \hfill \\
Q_{3 \cdot } \left( {\Psi z + \Pi \eta } \right) = 0\:.
\\
\end{gathered}
\end{equation}
Proceeding to the second block of equations in \eqref{eq:ctWformPartGen}, we see it is purely
explosive at a growth rate exceeding $\bar \xi$, so that as in the case
of the explosive component in the discrete time model, we must insist that
it follows its unique non-explosive solution, which is simply a constant.
Thus the full set of stability conditions is
\begin{equation}
\begin{bmatrix}
{w_2 } \\
{w_3 } \\
\end{bmatrix} = - \begin{bmatrix}
{\Omega_{22} } & {\Omega_{12} } \\
0 & {\Omega_{33} } \\
\end{bmatrix}^{ - 1} \begin{bmatrix}
{Q_{2 \cdot } C} \\
{Q_{3 \cdot } C} \\
\end{bmatrix}
\end{equation}
\begin{equation}\label{eq:ctGenzeta}
\begin{bmatrix}
{Q_{2 \cdot } } \\
{Q_{3 \cdot } } \\
\end{bmatrix}\left( {\Psi z + \Pi \eta } \right) = 0\:.
\end{equation}
Just as in the discrete case, we can now translate the equations
and stability conditions in terms of $w$ back in to a stable
equation in $y$, here taking the form
\begin{equation}\label{eq:ctSolGenTheta}
\dot y = \Theta _1 y + \Theta _c + \Theta _0 z\:.
\end{equation}
As in the discrete case also, \eqref{eq:ctGenzeta} will allow us, when conditions for
existence are met, to write
\begin{equation}
Q_{1 \cdot } \Pi \eta = \Phi \cdot
\begin{bmatrix}
{Q_{2 \cdot } \Pi } \\
{Q_{3 \cdot } \Pi } \\
\end{bmatrix}
\end{equation}
for some $\Phi$. Letting the subscript $u$ (for ``unstable") refer to the
range of indexes in blocks 2 and 3 of \eqref{eq:ctWformPartGen}, and then letting the
$u$ subscript play the role of the 2 subscript in the formulas of \eqref{eq:translate},
we get the correct answers for the coefficients of \eqref{eq:ctSolGenTheta}.
The conditions for existence and uniqueness are exactly the same as in the
discrete case without serial correlation in $z$, i.e. A$'$ and B.
\section{Continuous time, unrestricted $z$}
In discrete time, a first-order polynomial of the form $\Lambda+\Omega L$, in
which $\Lambda$ and $\Omega$ are upper triangular and for each $i$
either $\abs{\omega_{ii}/\lambda_{ii}}>1$ or
$\lambda_{ii}=0$, $\omega_{ii}\neq 0$, always
has a convergent ``forward" inverse, in the form of the sum that appears in
the right-hand side of \eqref{eq:ySolution}. In continuous time, the
operation analogous to inverting such an ``unstable" $\Lambda+\Omega L$ is
inverting a polynomial in the differential operator $D$ of the form $\Lambda D
+\Omega$ in which for every $i$, either the real part of
$\omega_{ii}/\lambda_{ii}$ is negative or $\lambda_{ii}=0$.
If there are 0's on the diagonal of $\Lambda$ in this case, the resulting
inverse is not simply a ``forward" operator in continuous time, with integrals
replacing sums, but is a convolution of such operators with finite-order
polynomials in the differential operator $D$. To be specific, if $\Lambda$ is
upper triangular and has zeros on its main diagonal, then
\begin{equation}\label{eq:lemma}
(I-\Lambda D)^{-1}=\sum_{s=0}^{n-1}\Lambda^s D^s\:,
\end{equation}
where $n$ is the order of the $\Lambda$ matrix.
This can be checked easily, and follows from the fact that for a $\Lambda$ of
this form it is guaranteed that $\Lambda^n=0$. Unless the upper triangle of
$\Lambda$ is dense, it is likely that terms in \eqref{eq:lemma} for larger $s$
turn out to be zero.
If as before we group blocks 2 and 3 of \eqref{eq:wFormPart} into a $u$
block, still using 1 to label the other block, we can write that equation, using the
differential operator $D$, as
\begin{equation}\label{eq:Dform}
\left(\begin{bmatrix}
\Lambda_{11} & \Lambda_{1u}\\
0 & \Lambda_{uu}
\end{bmatrix}D -\begin{bmatrix}
\Omega_{11} \Omega_{1u}\\
0 \Omega_{uu}
\end{bmatrix}\right)\begin{bmatrix}
w_1\\w_u
\end{bmatrix}=\begin{bmatrix}
Q_{1\cdot}\\Q_{2\cdot}
\end{bmatrix}\left(C+\Psi z + \Pi \eta\right)\:.
\end{equation}
Using the fact we observed in \eqref{eq:lemma}, we can see that the
differential operator on the left-hand side of this equation, in the lower
block, $\Lambda_uD-\Omega_{uu}$, has a stable inverse that is a convolution of
a finite-order matrix polynomial in $D$ with exponentially weighted averages
of future values. While it is not worthwhile to write out the whole inverse
explicitly, we can observe that the lower block can be solved recursively as
\begin{gather}\label{eq:Dsolve}
w_3=-(I-\Omega_{33}^{-1}\Lambda_{33}D)^{-1}\Omega_{33}^{-1}Q_{3\cdot}(C+\Psi
z(t))\\
\begin{split}w_2(t) = -\int_0^{\infty}e^{-\Lambda_{22}^{-1}\Omega_{22}s}\lambda_{22}^{-1}\cdot& \\
\bigl(Q_{2\cdot}\left(C+\Psi z(t+s)+\right.&\left.\Pi
\eta(t+s)\right)-\Lambda_{23}Dw_3(t+s)+\Omega_{23}w_3(t+s)\bigr)\, ds
\:.\label{eq:Fsolve}
\end{split}
\end{gather}
The inverted polynomial in \eqref{eq:Dsolve} becomes a finite order polynomial
in $D$ according to \eqref{eq:lemma}, and the $D$ operators appearing on right
of \eqref{eq:Fsolve} are interpreted as applying to $E_tz(t+s)$, considered as
a function of $s$ and (at $s=0$) taken as right-derivatives.
Existence and uniqueness questions are determined by exactly the same
conditions as in conditions A and B for the discrete case, with the $u$
subscript playing the role of the 2 subscript in those conditions. To get a
complete solution, we combine the recursive solution for $w_u$ in
\eqref{eq:Fsolve} and \eqref{eq:Dsolve} with the stable equation in $w_1$,
free of occurrences of $\eta$, that is obtained by multiplying
\eqref{eq:ctWformPartGen} by $[I\quad -\Phi]$, where $\Phi$ has been computed
to solve \eqref{eq:etaLink}.
The continuous time version of the matlab software, \texttt{gensysct}, does
not attempt to provide an explicit complete solution for this general case. It
checks existence and uniqueness and provides the matrices needed to write out
an explicit solution for the case of pure white noise $z$ in the form
\eqref{eq:ctSolGenTheta}, and it provides the QZ decomposition ordered as in
\eqref{eq:ctWformPartGen}.
Note that if one wants the two-sided (on future and past) projection of $y$ on $z$, and
if one is confident (perhaps because of having run \texttt{gensys}) that there
are no existence or uniqueness problems, one can find the projection by
Fourier methods directly. That is, one can form $\tilde \Xi(\omega)=i\omega
\Gamma_0-\Gamma_1$ and form the projection as the inverse Fourier transform of
$\tilde\Xi^{-1}\Psi$, where the inverse in this expression is
frequency-by-frequency matrix inversion. Of course this only produces stable
results if the projection does not involve delta functions or derivative
operators.
%------------------------------------------
\bibliographystyle{econometrica}
\bibliography{cas}
\end{document}
TRASH PILE
\footnote{%
To put the analysis of this section on a mathematically more
general and rigorous foundation, we would require that $z$ be the
time-derivative of a vector of semimartingales with absolutely
continuous non-martingale components and $\eta$ the time-derivative of
a martingale. Equation \eqref{eq:contGenEq} is then interpreted as defining a
relation among stochastic integrals. The requirement that $z$ be
free to have white noise components then is formalized as the
notion that $y(t)$ itself be measurable with respect to the $\sigma$-field
on which the model's processes are defined, while $z$ and $\eta$ are
allowed to produce components in $y$'s time-derivative that are not
random variables on this field.}