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