% ----------------------------------------------------------------
% 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.}
