\documentclass[12pt,a4paper]{article}
\usepackage{amsmath,amssymb,amsthm,natbib,graphicx,mathrsfs,tikz,amsfonts,mathdots}
\usetikzlibrary[decorations.pathreplacing]
\usetikzlibrary{patterns}

\linespread{1.6}
\addtolength{\hoffset}{-2cm}
\addtolength{\textwidth}{4cm}
\addtolength{\voffset}{-2cm}
\addtolength{\textheight}{3cm}


\begin{document}

\title{Online Supplement: The Linear Systems Approach to Linear Rational Expectations Models}
\author{Majid M.\ Al-Sadoon\\Universitat Pompeu Fabra \& Barcelona GSE}

\theoremstyle{definition}
\newtheorem{lem}{Lemma}[section]
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}{Corollary}[section]
\newtheorem{prop}{Proposition}[section]
\newtheorem{exmp}{Example}[section]
\newtheorem{defn}{Definition}[section]
\newtheorem{claim}{Claim}[section]
\newtheorem{alg}{Algorithm}[section]
\newtheorem{rmk}{Remark}[section]
\newtheorem{res}{Result}[section]


\renewcommand{\labelenumi}{(\roman{enumi})}


\maketitle

\abstract{
This document contains: (i) supplementary results for the ILWHF, (ii) additional proofs of results from the paper, and (iii) an algorithm for computing the ILWHF.

\appendix
\section{Supplementary Results on the ILWHF}
\subsection{Existence of ILWHF}

Various proofs of the existence of WHF can be adopted to prove the existence of ILWHF. The simplest approaches can be found in \cite{ggk1} and \cite{gks}. However, a much more direct proof can be constructed utilizing the concept of column properness from linear system theory. The importance of this proof is that: (i) it clarifies the structure of ILWHFs of Laurent matrix polynomials, (ii) it demonstrates a simple relationship between the Smith-McMillan form of a rational matrix and the Smith-McMillan forms of its backward component, and (iii) it forms the basis of the numerical implementation of the ILWHF presented in Section \ref{sec:comp}.

Column properness, a concept commonly attributed to \cite{wolovich}, can be explained as follows. For $M(z)\in\mathbb{R}^{n\times n}[z]$, let $\nu_i$ denotes the degree of the $i$-th column of $M(z)$. Then we may express $M(z)$ uniquely as $M_{hc}\mathrm{diag}(z^{\nu_1},\ldots,z^{\nu_n})+L(z)$, where the degree of the $i$-th column of $L(z)$ is strictly less than $\nu_i$. We say that $M(z)$ is column proper (or column reduced) if $M_{hc}$ is of full rank. Clearly, $\deg(\det(M(z)))=\sum_{i=1}^n\nu_i$ if and only if $M(z)$ is column proper. There are many instances in linear system theory when we are interested in factoring $z^{\nu_i}$ from the $i$-th column of $M(z)$ to arrive at $N(z)=M_{hc}+L(z)\mathrm{diag}(z^{-\nu_1},\ldots,z^{-\nu_n})\in\mathbb{R}^{n\times n}(z)$. Constructed in this way, $N(z)$ can be ensured to have no pole at infinity but it cannot be ensured to have no zero at infinity. $N(z)$ will have no zero at infinity if and only if $M_{hc}$ is of full rank; that is, if and only if $M(z)$ is column proper. We recall, for future reference, that every non-singular $M(z)\in\mathbb{R}^{n\times n}[z]$ can be brought to column proper form by either left or right multiplication by a unimodular matrix \cite[Theorem 2.5.14]{wolovich}. 

\begin{thm}\label{thm:whexistence}
Let $M(z)\in\mathbb{R}^{n\times n}(z)$ be non-singular and $\rho>0$, then an ILWHF exists for $M(z)$ relative to $\rho \mathbb{T}$.
\end{thm}
\begin{proof}
Let $q(z)$ be the greatest common denominator of all of the elements of $M(z)$, and let $P(z)=q(z)M(z)$. Then, an ILWHF of $M(z)=M_f(z)M_0(z)M_b(z)$ is obtained from any ILWHFs of $P(z)=P_f(z)P_0(z)P_b(z)$ and $q(z)=q_f(z)q_0(z)q_b(z)$ as $M_f(z)=P_f(z)/q_f(z)$, $M_0(z)=P_0(z)/q_0(z)$, and $M_b(z)=P_b(z)/q_b(z)$. Thus, the existence of an ILWHF for a non-singular rational matrix follows from the existence of an ILWHF for a non-singular polynomial matrix.

Let $P(z)\in\mathbb{R}[z]$. The result is trivial if $P(z)$ is a non-zero constant. Therefore assume $P(z)$ has a non-empty set of zeros, $\{\zeta_i\}$. Counting multiplicities, let $\kappa_1$ be the number of roots of $P(z)$ inside $\rho \mathbb{D}$ and let $K$ be the leading coefficient of $P(z)$.
\begin{align*}
P(z)&=\underbrace{\prod_{|\zeta_i|<\rho}(z-\zeta_i)}_{P_{\rho \mathbb{D}}(z)}\underbrace{K\prod_{|\zeta_i|\geq\rho}(z-\zeta_i)}_{P_{\rho \mathbb{D}^c}(z)}=\underbrace{\prod_{|\zeta_i|<\rho}(1-\zeta_iz^{-1})}_{P_f(z)}\underbrace{\vphantom{\prod_{0<|\zeta_i|<\rho}}z^{\kappa_1}}_{P_0(z)}\underbrace{K\prod_{|\zeta_i|\geq\rho}(z-\zeta_i)}_{P_b(z)}.
\end{align*}
It is clear that an ILWHF with respect to $\rho\mathbb{T}$ is obtained. Notice that two steps are required for the factorization, a polynomial factorization relative to $\rho \mathbb{T}$ into $P_{\rho \mathbb{D}}(z)P_{\rho \mathbb{D}^c}(z)$, followed by division of $P_{\rho \mathbb{D}}(z)$ by its degree.

The factorization of a non-singular $P(z)\in\mathbb{R}^{n\times n}[z]$ follows exactly the same logic. First, obtain the Smith form of $P(z)=U(z)\Lambda(z)V(z)$. Next, obtain the ILWHF relative to $\rho \mathbb{T}$ of the $i$-th diagonal element of $\Lambda(z)$ as $\Lambda_{ii}(z)=\Lambda_{iif}(z)\Lambda_{ii0}(z)\Lambda_{iib}(z)$ and set
\begin{align*}
P_{\rho \mathbb{D}}(z)&=U(z)\mathrm{diag}(\Lambda_{11f}(z)\Lambda_{110}(z),\ldots,\Lambda_{nnf}(z)\Lambda_{nn0}(z))\\
P_{\rho \mathbb{D}^c}(z)&=\mathrm{diag}(\Lambda_{11b}(z),\ldots,\Lambda_{nnb}(z))V(z).
\end{align*}
Thus, $P(z)=P_{\rho \mathbb{D}}(z)P_{\rho \mathbb{D}^c}(z)$, where $P_{\rho \mathbb{D}}(z)\in\mathbb{R}^{n\times n}[z]$ contains all the zeros of $P(z)$ that are in $\rho \mathbb{D}$ and $P_{\rho \mathbb{D}^c}(z)\in\mathbb{R}^{n\times n}[z]$ contains all the zeros in $\rho \mathbb{D}^c$. We may now attempt to divide each column of $P_{\rho \mathbb{D}}(z)$ by its degree to form $P_b(z)$. However, if $P_{\rho \mathbb{D}}(z)$ is not column reduced, this may result in a rational matrix that has a zero at infinity. Thus, let $W(z)\in\mathbb{R}^{n\times n}[z]$ be a unimodular matrix such that $P_{\rho \mathbb{D}}(z)W(z)$ is column proper \citep[Theorem 2.5.7]{wolovich} and let $\Pi\in\mathbb{R}^{n\times n}$ be a permutation matrix such that $P_{\rho \mathbb{D}}(z)W(z)\Pi$ has column degrees $\kappa_1\geq\cdots\geq\kappa_n$. Then $P(z)=P_f(z)P_0(z)P_b(z)$, with $P_0(z)=\mathrm{diag}(z^{\kappa_1},\ldots,z^{\kappa_n})$, $P_f(z)=P_{\rho \mathbb{D}}(z)W(z)\Pi P_0^{-1}(z)$, and $P_b(z)=\Pi^{-1}W^{-1}(z)P_{\rho \mathbb{D}^c}(z)$. To see that this is an ILWHF with respect to $\rho\mathbb{T}$, note that any finite zeros or poles of $P_f(z)$ occur only in $\rho\mathbb{D}$ and due to the column reduction step and subsequent division by the column degrees, $P_f(z)$ has no zeros or poles at infinity; on the other hand, $P_b(z)$ has all its zeros and poles in $\rho\mathbb{D}^c$.
\end{proof}

\begin{cor}\label{cor:1}
Under the assumptions of Theorem \ref{thm:whexistence}, the Smith-McMillan form of a backward component of $M(z)$ is a backward component of the Smith-McMillan form of $M(z)$.
\end{cor}
\begin{proof}
Using the notation above, $M_b(z)=\Pi^{-1}W^{-1}(z)\mathrm{diag}(\Lambda_{11b}(z),\ldots,\Lambda_{nnb}(z))V(z)/q_b(z)$, which implies that the Smith-McMillan form of $M_b(z)$ is $\mathrm{diag}(\Lambda_{11b}(z),\ldots,\Lambda_{nnb}(z))/q_b(z)$.
\end{proof}

\subsection{Uniqueness of ILWHF}

Having proven existence, the next question that ought to be answered concerns the uniqueness of the factorization.

\begin{thm}\label{thm:whuniqueness}
Let $M(z)\in\mathbb{R}^{n\times n}(z)$ be non-singular, $\rho>0$, and let $M(z)=M_f(z)M_0(z)M_b(z)$ be an ILWHF relative to $\rho\mathbb{T}$, then $\breve M_f(z)\breve M_0(z)\breve M_b(z)$ is also an ILWHF relative to $\rho\mathbb{T}$ if and only if $M_0(z)=\breve M_0(z)$, $\breve M_f(z)=M_f(z)M_0(z)U^{-1}(z)M_0^{-1}(z)$, and $\breve M_b(z)=U(z)M_b(z)$, where $U(z)\in\mathbb{R}^{n\times n}[z]$ is unimodular, $U_{ij}(z)=0$ for $\kappa_i>\kappa_j$, $U_{ij}(z)\in\mathbb{R}$ for $\kappa_i=\kappa_j$, and $\deg(U_{ij}(z))\leq\kappa_j-\kappa_i$ for $\kappa_i<\kappa_j$.
\end{thm}
\begin{proof}
By Proposition 3.1 (i), both ILWHFs can be considered WHFs relative to a contracted contour. The result then follows from the corresponding result for WHF (e.g.\ Theorems I.1.1 and I.1.2 of \cite{cg}).
\end{proof}

It follows from Theorem \ref{thm:whuniqueness} that the partial indices of a non-singular matrix rational function are unique but its forward and backward components are not. The backward component is determined only up to left multiplication by $U(z)$, a non-singular, block lower triangular matrix polynomial, with constant blocks on the diagonal, and subdiagonal blocks of bounded degrees. The structure of $M_0(z)U^{-1}(z)M_0^{-1}(z)$, which determines the equivalence class of the forward component, is similar in that it is a non-singular, block lower triangular matrix polynomial in $z^{-1}$, with constant blocks on the diagonal, and subdiagonal blocks of bounded degrees in $z^{-1}$.

An important special case of the theorem occurs when the partial indices are all zero, in which case $\breve M_f(z)\breve M_b(z)$ is an ILWHF of $M(z)$ relative to $\rho\mathbb{T}$ if and only if there exists an invertible matrix $U\in\mathbb{R}^{n\times n}$ such that $M_f(z)=\breve M_f(z)U^{-1}$ and $M_b(z)=U\breve M_b(z)$. In that case, a unique choice of ILWHF can be obtained by setting either $M_f(\infty)=I_n$ or $M_b(0)=I_n$.

\subsection{ILWHF of Generic Systems}
A natural question that arises in relation to the ILWHF is whether the partial indices are stable under small perturbations of the matrix function. To answer this question, we endow $\mathbb{R}_{pq}^{n\times n}(z)$ with the metric $d(M(z),N(z))=\sum_{i=-q}^p\|M_i-N_i\|$. \cite{gk} prove a general result that specializes in our context to the fact that an infinitesimally small perturbation to an $M(z)\in\mathbb{R}_{pq}^{n\times n}(z)$ with no zeros on $\rho\mathbb{T}$ leaves its partial indices unchanged if and only if its largest and smallest partial indices differ by no more than $1$. \cite{gks} provide a simplified proof of this result. A generalization of this result for the ILWHF is the following.

\begin{thm}\label{thm:stability}
For fixed $\rho>0$, non-negative integers $p$ and $q$, and $n\geq1$, the set of all non-singular $M(z)\in\mathbb{R}_{pq}^{n\times n}(z)$ whose partial indices in ILWHFs relative to $\rho \mathbb{T}$ satisfy $\kappa_1\leq\kappa_n+1$ contains an open and dense subset.
\end{thm}
\begin{proof}
Define
\begin{align*}
\frak{S}_{pq}&=\{M(z)\in\mathbb{R}_{pq}^{n\times n}(z):\det(M(z))\neq0, z\in \rho \mathbb{T}\},\\
\overline{\frak S}_{pq}&=\{M(z)\in\mathbb{R}_{pq}^{n\times n}(z): M(z)\text{ is non-singular}\},\\
A&=\{M(z)\in \frak S_{pq}: \kappa_1\leq\kappa_n+1\},\\
\overline A&=\{M(z)\in \overline{\frak S}_{pq}: \kappa_1\leq\kappa_n+1\}. 
\end{align*}
$A\subset\overline A$, since $\frak{S}_{pq}\subset\overline{\frak S}_{pq}$. We now claim that $A$ is open and dense in $\bar{\frak{S}}_{pq}$. The fact that it is open follows from the fact that $A$ is open in $\frak{S}_{pq}$ \citep{gk} and the fact that $\frak{S}_{pq}$ is open in $\bar{\frak{S}}_{pq}$. The latter fact follows from continuity of the set of zeros of a polynomial as a function of its own coefficients \citep[p.\ 166]{stewartsun}. The fact that $A$ is dense in $\bar{\frak{S}}_{pq}$, follows from the fact that $A$ is dense in $\frak{S}_{pq}$ \citep{gk} and the fact that $\frak{S}_{pq}$ is dense in $\overline{\frak{S}}_{pq}$. The latter fact follows from taking $r$ to be as in Proposition 3.1, then $M((r/\rho)z)\in\mathbb{R}_{pq}^{n\times n}(z)$ has no zeros on $\rho T$ and can be made arbitrarily close to $M(z)$.
\end{proof}

Theorem \ref{thm:stability} implies that a generic non-singular Laurent polynomial has partial indices that satisfy $\kappa_1\leq\kappa_n+1$ in its ILWHF relative to $\rho\mathbb{T}$, thus its partial indices are either all non-negative or all non-positive. This has important implications for the existence and uniqueness of solutions to generic LREMs.

\section{Additional Proofs}
\begin{proof}[Proof of Lemma 1 of the Paper]
Before we begin, we will need to state the following useful inequality. For $Z\in\mathcal{S}^n$, $t\in\mathbb{Z}$, and $0\leq\varphi<1$, there is a constant $C(Z,t,\varphi)>0$ such that
\begin{align}
\varphi^{t+i}E\|Z_{t+i}\|\leq C(Z,t,\varphi),\qquad i\geq0.\tag{$\ast$}\label{eq:ineq}
\end{align}

(i) Since $N(z)$ has no poles in $\mathbb{D}^c$, there exists a $\rho\in(0,1)$ such that $N(z)$ has a Laurent series expansion $\sum_{i=0}^\infty N_iz^{-i}$ in $\rho\mathbb{D}^c$. Now the Monotone Convergence Theorem \citep[Theorem 5.3]{williams} and inequality \eqref{eq:ineq} imply that $E\left(\sum_{i=0}^\infty\|N_i\|\|Y_{t+i}\|\right)=\sum_{i=0}^\infty\|N_i\|E\|Y_{t+i}\|\leq \sum_{i=0}^\infty\|N_i\|\varphi^{-t-i}C(Y,t,\varphi)=C(Y,t,\varphi)\varphi^{-t}\sum_{i=0}^\infty\|N_i\|\varphi^{-i}$ for $0<\varphi<1$. If we choose $\varphi$ so that $\rho<\varphi<1$, then it lies in the convergence region of the Laurent series for $N(z)$ and so $\sum_{i=0}^\infty\|N_i\|\varphi^{-i}<\infty$. It follows that $E\left(\sum_{i=0}^\infty\|N_i\|\|Y_{t+i}\|\right)<\infty$ for all $t\in\mathbb{Z}$. This implies that $\sum_{i=0}^\infty\|N_i\|\|Y_{t+i}\|$ converges almost surely for all $t\in\mathbb{Z}$ \cite[Result 6.5.(c)]{williams} and therefore $\sum_{i=0}^\infty N_iY_{t+i}$ converges almost surely and is in $\mathcal{L}^1$ for all $t\in\mathbb{Z}$. To prove that $N(L)Y\in\mathcal{S}^n$, let $0\leq\theta<1$ and now choose $\varphi$ such that $\max\{\rho,\theta\}<\varphi<1$, then inequality \eqref{eq:ineq} implies that $\theta^tE\|\sum_{i=0}^\infty N_iY_{t+i}\|\leq \theta^t\sum_{i=0}^\infty \|N_i\|E\|Y_{t+i}\|\leq \theta^t\sum_{i=0}^\infty\varphi^{-t-i}C(Y,0,\varphi)\|N_i\|=C(Y,0,\varphi)(\theta/\varphi)^t\sum_{i=0}^\infty\|N_i\|\varphi^{-i}$ for all $t\geq0$. Since $\sum_{i=0}^\infty\|N_i\|\varphi^{-i}<\infty$, the last term tends to zero as $t\rightarrow\infty$ so $N(L)Y\in\mathcal{S}^n$.

(ii) Just as in (i), there is a $\rho\in(0,1)$ such that both $N(z)$ and $M(z)$ have Laurent series expansions in $\rho\mathbb{D}^c$. Denote the Laurent series expansion for $M(z)$ in $\rho\mathbb{D}^c$ by $\sum_{i=0}^\infty M_iz^{-i}$. Now absolute summability implies that the order of summation of a series is irrelevant \cite[Theorem 3.55]{rudin}. A simple extension of that result implies that for $t\in\mathbb{Z}$, if $\sum_{j=0}^\infty\sum_{i=0}^\infty \|M_j\|\|N_i\|\|Y_{t+i+j}\|<\infty$ a.s., then $M(L)(N(L)Y_t)=\sum_{j=0}^\infty M_j\sum_{i=0}^\infty N_iY_{t+i+j}=\sum_{k=0}^\infty (\sum_{i+j=k}M_jN_i) Y_{t+k}=(M(L)N(L))Y_t$ a.s.\ By the Monotone Convergence Theorem and inequality \eqref{eq:ineq} again, $E\left(\sum_{j=0}^\infty\sum_{i=0}^\infty \|M_j\|\|N_i\|\|Y_{t+i+j}\|\right)=\sum_{j=0}^\infty\sum_{i=0}^\infty \|M_j\|\|N_i\|E\|Y_{t+i+j}\|\leq\sum_{j=0}^\infty\sum_{i=0}^\infty \|M_j\|\|N_i\|\varphi^{-t-i-j}C(Y,t,\varphi)\leq\varphi^{-t}C(Y,t,\varphi)\sum_{j=0}^\infty\|M_j\|\varphi^{-j}\sum_{i=0}^\infty \|N_i\|\varphi^{-i}$ for $0<\varphi<1$. This last term is finite if we choose $\rho<\varphi<1$. Thus $\sum_{j=0}^\infty\sum_{i=0}^\infty \|M_j\|\|N_i\|\|Y_{t+i+j}\|<\infty$ a.s.\ by Result 6.5.(c) of \cite{williams}. Finally, the fact that $M(L)N(L)Y\in \mathcal{S}^n$ follows from (i).

(iii) By the conditional version of Jensen's inequality and the Tower Property \cite[Properties 9.7.(h) and 9.7.(i)]{williams}, $E\|E(N(L)Y_t|\mathscr{I}_t)\|\leq E(E(\|N(L)Y_t\||\mathscr{I}_t))=E\|N(L)Y_t\|$; the fact that $\{E(N(L)Y_t|\mathscr{I}_t):t\in\mathbb{Z}\}\in\mathcal{S}^n$ then follows from (i). Next, we found in (i) that $\sum_{i=0}^\infty N_iY_{t+i}$ converges a.s.\ and is bounded by $\sum_{i=0}^\infty\|N_i\|\|Y_{t+i}\|\in \mathcal{L}^1$, thus, the conditional version of the Dominated Convergence Theorem \cite[Property 9.7.(g)]{williams} implies that $\sum_{i=0}^\infty N_iE(Y_{t+i}|\mathscr{I}_t)=E(N(L)Y_t|\mathscr{I}_t)$ a.s.
\end{proof}

\begin{proof}[Proof of Lemma 2 of the Paper]
(i) Let $N(z)=\sum_{i=0}^pN_iz^i$. Then $\det(N_0)=\det(N(0))\neq0$ because $0\in \mathbb{D}$. Therefore, the process $X$ can be obtained recursively as $X_t=-N_0^{-1}N_1X_{t-1}-\cdots-N_0^{-1}N_pX_{t-p}+N_0^{-1}Y_t$ for all $t\geq0$, with $X_t=\tilde X_t$ for $t<0$. To prove that it is in $\mathcal{S}^n$, first note that $X_t\in \mathcal{L}^1$ for all $t\geq0$. Now let $0<\theta<1$ and define $Q(z)=N(\theta z)$. Then $Q(L)(\theta^tX_t)=\theta^tN(L)X_t=\theta^tY_t$ for all $t\geq0$. Since $\det(Q(z))\neq0$ for all $|z|<\theta^{-1}$, $Q^{-1}(z)=\sum_{i=0}^\infty Q^iz^i$ converges in $\theta^{-1}\mathbb{D}$ and $\theta^tX_t=G_{-1t}\tilde X_{-1}+\cdots+G_{-pt}\tilde X_{-p}+\sum_{i=0}^t Q^i\theta^{t-i}Y_{t-i}$ for $t\geq0$, where the matrices $G_{it}$ are exponentially decaying. It follows that $E\|\theta^tX_t\|\leq \|G_{-1t}\|E\|\tilde X_{-1}\|+\cdots+\|G_{-pt}\|E\|\tilde X_{-p}\|+\sum_{i=0}^t \|Q^i\|\theta^{t-i}E\|Y_{t-i}\|$ for $t\geq0$. Since $Y\in\mathcal{S}^n$, inequality \eqref{eq:ineq} implies that $\sum_{i=0}^t \|Q^i\|\theta^{t-i}E\|Y_{t-i}\|\leq \sum_{i=0}^t \|Q^i\|(\theta/\varphi)^{i-t}C(Z,0,\varphi)=C(Z,0,\varphi)(\theta/\varphi)^{-t}\sum_{i=0}^t \|Q^i\|(\theta/\varphi)^i$ for $0<\varphi<1$ and $t\geq0$. If we further require that $\theta^2<\varphi<\theta$, then $\sum_{i=0}^\infty \|Q^i\|(\theta/\varphi)^i<\infty$ and $(\theta/\varphi)^{-t}\rightarrow0$, so $E\|\theta^tX_t\|\rightarrow0$ and $X\in\mathcal{S}^n$.

(ii) If $\hat X$ is another solution, then $N(L)(X_t-\hat X_t)=0$ a.s.\ for $t\geq0$. Since $\det(N_0)\neq0$, $X_{\bar t+1}=\hat X_{\bar t+1}$ a.s.\ whenever $X_t=\hat X_t$ a.s.\ for $t\leq\bar t$. But $X_t=\hat X_t$ a.s.\ for $t<0$, therefore $X_t=\hat X_t$ a.s.\ for all $t\in\mathbb{Z}$.
\end{proof}

\begin{proof}[Proof of Proposition 5.2 of the Paper]
Let $t\in\mathbb{Z}$ and note that
\begin{align*}
\Psi P_t=\left[\begin{array}{ccccc}
M_1 & M_2 & \cdots & M_{p-1} & M_p\\
M_2 & M_3 & \cdots & M_p & 0\\
\vdots & \vdots & \iddots & \iddots & \vdots\\
M_{p-1} & M_p & \iddots &  & \vdots\\
M_p & 0 & \cdots & \cdots & 0\\
\hline
0 & 0 & \cdots & \cdots & 0\\
\vdots & \vdots & & & \vdots
\end{array}\right]
\left[\begin{matrix}
X_t\\
X_{t-1}\\
\vdots\\
X_{t+1-p}
\end{matrix}\right].
\end{align*}
Thus $\Psi P_t$ is a.s.\ bounded over the first $p$ blocks and equal to zero in all subsequent blocks of the sequence. It follows that $\Psi P_t\in l_n^\infty$ a.s. The fact that $F_{t+1|t}$ is a solution to $\Theta F_{t+1|t}+\Psi P_t=0$ has already been demonstrated. To see that $F_{t+1|t}\in l_n^\infty$ a.s.\ note that because $X$ is generated by a VAR, we can use basic state space methods (see e.g.\ Section 2.2 of \cite{lutkepohl}) to find a triple $(A,B,C)\in\mathbb{R}^{n\times k}\times\mathbb{R}^{k\times k}\times\mathbb{R}^{k\times pn}$ such that
\begin{align*}
F_{t+1|t}
=
\left[\begin{matrix}
CAB\\
CA^2B\\
CA^3B\\
\vdots
\end{matrix}\right]
\left[\begin{matrix}
X_t\\
X_{t-1}\\
\vdots\\
X_{t+1-p}
\end{matrix}\right].
\end{align*}
Since $M_b(z)$ is stable, $A$ can be chosen to have all its eigenvalues inside the unit circle. Thus, $F_{t+1|t}$ is a.s.\ an element of $l_n^\infty$. Finally, since $M(z)$ has no zeros or poles on $\mathbb{T}$, its ILWHF relative to $\mathbb{T}$ is also a WHF and since, additionally, its partial indices are zeros, this implies that $\Theta$ is an invertible operator on $l_n^\infty$ \citep[Theorem VIII.4.2]{gf}. Thus $F_{t+1|t}$ is a.s.\ the unique solution in $l_n^\infty$ of $\Theta F_{t+1|t}+\Psi P_t=0$.
\end{proof}

\section{Computing the ILWHF}\label{sec:comp}
There is surprisingly little in the linear systems or the linear operators literature on the computation of WHFs. The methods used to prove Theorem \ref{thm:whexistence} are non-starters due to their high complexity. Results using state space methods are available but, as far as the author is aware, these tend to impose restrictive assumptions on $M(z)\in\mathbb{R}_{pq}^{n\times n}(z)$ such as partial indices that are identically zero, strict properness, and non-singularity at infinity \citep{ggk2,gks}. \cite{adukov} suggested a method for obtaining the WHF that requires computing moments of $M^{-1}(z)$. We will derive a simpler solution that connects nicely to the proof of Theorem \ref{thm:whexistence} and to the \cite{sims} method for solving LREMs. The algorithm is implemented in the Matlab program \texttt{ilwhf.m} accompanying this paper.

Let $M(z)\in\mathbb{R}_{pq}^{n\times n}(z)$ be non-singular and $\rho>0$, then clearly $M(z)=M_f(z)M_0(z)M_b(z)$ is an ILWHF relative to $\rho \mathbb{T}$ if and only if $z^{q+1}M(z)=M_f(z)(z^{q+1}M_0(z))M_b(z)$ is also an ILWHF relative to $\rho \mathbb{T}$. Therefore, we can obtain an ILWHF of $M(z)$ relative to $\rho \mathbb{T}$ from that of $z^{q+1}M(z)$. Note that although multiplying $M(z)$ by $z^q$ is sufficient to turn it into a matrix polynomial on which the polynomial operations below are applicable, the fact that $\deg(z^{q+1}M(z))\geq1$ will serve a crucial purpose later on. Define the following matrices
\begin{gather*}
\Gamma_0=\left[\begin{array}{c|ccccccc}
0 & M_{-q} & \cdots & M_{p-1}\\
\hline
0 \\
\vdots & & I_{n(p+q)}\\
0
\end{array}\right],\quad \Gamma_1=\left[\begin{array}{ccc|c}
0 & \cdots & 0 & M_p\\
\hline
& & & 0 \\
& -I_{n(p+q)} & & \vdots \\
& & & 0 
\end{array}\right],
\end{gather*}
and $\Gamma(z)=\Gamma_0+\Gamma_1z\in\mathbb{R}^{l\times l}[z]$, with $l=n(p+q+1)$. Then, $E(z)\Gamma(z)=\left[\begin{smallmatrix}
z^{q+1}M(z) & 0 \\ 0 & I_{n(p+q)}
\end{smallmatrix}\right]F(z)$, where
\begin{gather*}
E(z)=\left[\begin{array}{c|ccccccc}
I_n & E_1(z) & \cdots & E_{p+q}(z)\\
\hline
0 \\
\vdots & & I_{n(p+q)}\\
0
\end{array}\right], 
F(z)=\left[\begin{matrix}
I_n   & 0   & \cdots & 0\\
-zI_n & I_n & \cdots & 0\\
\vdots & \ddots & \ddots & \vdots\\
0 & \cdots & -zI_n & I_n
\end{matrix}\right],
\end{gather*}
and $E_{p+q+1}(z)=-M_p$ and $E_i(z)=-M_{i-q-1}+zE_{i+1}(z)$ for $i=1,\ldots,p+q$. Since $\det(E(z))=\det(F(z))=1$, it follows that $\det(z^{q+1}M(z))=\det(\Gamma(z))$ so $\Gamma(z)$ is non-singular. \cite{ggk1} refer to $E(z)\Gamma(z)F^{-1}(z)$ as a linearisation of $z^{q+1}M(z)$.

Our construction will proceed as in the proof of Theorem \ref{thm:whexistence}. We will factor $\Gamma(z)$ as $\Gamma_{\rho \mathbb{D}}(z)\Gamma_{\rho \mathbb{D}^c}(z)$, where $\Gamma_{\rho \mathbb{D}}(z)\in\mathbb{R}^{l\times l}[z]$ has zeros only in $\rho \mathbb{D}$ and $\Gamma_{\rho \mathbb{D}^c}(z)\in\mathbb{R}^{l\times l}[z]$ has zeros only in $\rho \mathbb{D}^c$. By Theorem 2.5.7 of \cite{wolovich}, there exists a unimodular matrix $W(z)\in\mathbb{R}^{l\times l}[z]$ such that $E(z)\Gamma_{\rho \mathbb{D}}(z)W(z)$ is column proper. Let $\Pi$ be a permutation matrix so that the column degrees of $E(z)\Gamma_{\rho \mathbb{D}}(z)W(z)\Pi$ are $\nu_1\geq\cdots\geq\nu_l\geq0$. Then, an ILWHF relative to $\rho\mathbb{T}$ of $N(z)=\left[\begin{smallmatrix}
z^{q+1}M(z) & 0 \\ 0 & I_{n(p+q)}
\end{smallmatrix}\right]=N_f(z)N_0(z)N_b(z)$ is given by
\begin{align*}
N_f(z)&=E(z)\Gamma_{\rho \mathbb{D}}(z)W(z)\Pi\mathrm{diag}(z^{-\nu_1},\ldots,z^{-\nu_l})\\
N_0(z)&=\mathrm{diag}(z^{\nu_1},\ldots,z^{\nu_l})\\
N_b(z)&=\Pi^{-1}W^{-1}(z)\Gamma_{\rho \mathbb{D}^c}(z)F^{-1}(z). 
\end{align*}
Theorem \ref{thm:whuniqueness} then implies that $N_0(z)=\left[\begin{smallmatrix} z^{q+1}M_0(z) & 0 \\ 0 & I_{l-n}\end{smallmatrix}\right]$. Thus, the partial indices of $M(z)$, satisfy $\kappa_i=\nu_i-q-1$ for $i=1,\ldots, n$ and $\nu_i=0$ for $i=n+1,\ldots,l$. We now claim that we may read $M_f(z)$ directly from the top left $n\times n$ block of $N_f(z)$ and likewise $M_b(z)$ from the top left $n\times n$ block of $N_b(z)$.  To see this, note that the unimodular transformation that determines the set of all ILWHFs of $N(z)$ relative to $\rho \mathbb{T}$ in Theorem \ref{thm:whuniqueness} takes the form $U(z)=\left[\begin{smallmatrix} U_{11}(z) & 0 \\ U_{21}(z) & U_{22} \end{smallmatrix}\right]$. Here, $U_{22}\in\mathbb{R}^{(l-n)\times (l-n)}$ is invertible and $U_{11}(z)\in\mathbb{R}^{n\times n}[z]$ is of the general type of unimodular transformation that $M(z)=\breve M_f(z)\breve M_0(z)\breve M_b(z)$, with
\begin{align*}
\breve M_f(z)&=M_f(z)M_0(z)U_{11}^{-1}(z)M_0^{-1}(z)\\
\breve M_0(z)&=M_0(z)\\
\breve M_b(z)&=U_{11}(z)M_b(z) 
\end{align*}
is an ILWHF relative to $\rho \mathbb{T}$. This partitioning of $U(z)$ relies crucially on the fact that the non-zero partial indices of $N(z)$ are bounded below by 1, which is made possible by the extra power of $z$ we mentioned earlier; without this extra power of $z$, the number of columns (or rows) of $U_{11}(z)$, which is equal to the number of non-zero partial indices of $N(z)$, can be of smaller size than $n$ and the algorithm cannot proceed. Since
\begin{align*}
N(z)=\left[\begin{matrix} M_f(z) & 0 \\ 0 & I_{l-n}\end{matrix}\right]\left[\begin{matrix} z^{q+1}M_0(z) & 0 \\ 0 & I_{l-n}\end{matrix}\right]\left[\begin{matrix} M_b(z) & 0 \\ 0 & I_{l-n}\end{matrix}\right]
\end{align*}
is an ILWHF relative to $\rho \mathbb{T}$. It follows that
\begin{align*}
N_f(z)&=\left[\begin{matrix} M_f(z) & 0 \\ 0 & I_{l-n}\end{matrix}\right]\left[\begin{matrix} M_0(z)U_{11}^{-1}(z)M_0^{-1}(z) & 0 \\ -U_{22}^{-1}U_{21}(z)U_{11}^{-1}(z)M_0^{-1}(z) & U_{22}^{-1}\end{matrix}\right]\\
N_b(z)&=\left[\begin{matrix} U_{11}(z) & 0 \\ U_{12}(z) & U_{22}\end{matrix}\right]\left[\begin{matrix} M_b(z) & 0 \\ 0 & I_{l-n}\end{matrix}\right]
\end{align*}
for some unimodular $U(z)$ of the aforementioned form. But now note that the top left $n\times n$ blocks of $N_f(z)$ and $N_b(z)$ have the forms $\breve M_f(z)$ and $\breve M_b(z)$ respectively.

Now all that remains is to factorize $\Gamma(z)$. This can be accomplished using the real QZ decomposition. By Theorem VI.1.9 and Exercise IV.1.3 of \cite{stewartsun}, there are orthogonal matrices $Q,Z\in\mathbb{R}^{l\times l}$ such that
\begin{align*}
Q\Gamma_0 Z&=\left[\begin{matrix} \Lambda_{11} & \Lambda_{12} \\ 0 & \Lambda_{22} \end{matrix}\right]& Q\Gamma_1 Z&=\left[\begin{matrix} \Omega_{11} & \Omega_{12} \\ 0 & \Omega_{22} \end{matrix}\right] 
\end{align*}
are partitioned conformably, $\det(\Lambda_{11}+\Omega_{11}z)$ has all its zeros in $\rho \mathbb{D}^c$, and $\det(\Lambda_{22}+\Omega_{22}z)$ has all its zeros in $\rho \mathbb{D}$ (thus, $\Omega_{22}$ is non-singular). It follows that 
\begin{align*}
\Gamma(z)&=\Gamma_0+\Gamma_1z\\
&=Q'(\Lambda+\Omega z)Z'\\
&=Q'\left(\left[\begin{matrix} I_s & 0 \\ 0 & \Lambda_{22} \end{matrix}\right]+\left[\begin{matrix} 0 & 0 \\ 0 & \Omega_{22} \end{matrix}\right]z\right)
\left(\left[\begin{matrix} \Lambda_{11} & \Lambda_{12} \\ 0 & I_{l-s} \end{matrix}\right]+\left[\begin{matrix} \Omega_{11} & \Omega_{12} \\ 0 & 0 \end{matrix}\right]z\right)Z'
\end{align*}
and so we may take
\begin{align*}
\Gamma_{\rho \mathbb{D}}(z)&=Q'\left(\left[\begin{matrix} I_s & 0 \\ 0 & \Lambda_{22} \end{matrix}\right]+\left[\begin{matrix} 0 & 0 \\ 0 & \Omega_{22} \end{matrix}\right]z\right)\\
\Gamma_{\rho \mathbb{D}^c}(z)&=\left(\left[\begin{matrix} \Lambda_{11} & \Lambda_{12} \\ 0 & I_{l-s} \end{matrix}\right]+\left[\begin{matrix} \Omega_{11} & \Omega_{12} \\ 0 & 0 \end{matrix}\right]z\right)Z'.
\end{align*}
The existence of this representation is proven under more general conditions in \cite{gk87}. The notation here highlights the connection to the \cite{sims} method, which evidently implicitly computes an ILWHF.

The implementation of the algorithm in \texttt{ilwhf.m} uses the Matlab command \texttt{qz} with option \texttt{`real'} to obtain an initial QZ decomposition, then uses the Matlab command \texttt{ordqz} to obtain the final form discussed above. The column reduction part of the algorithm utilizes an implementation of the \cite{geurts} correction to the \cite{kc} algorithm, which requires a tolerance to be specified.

\begin{exmp}\label{exmp:consc}
Consider Example 3.5 with $R=1.05$, corresponding to an interest rate of $5\%$ per period. Then we compute the ILWHF as follows:
\begin{verbatim}
>> M(:,:,1)=[1 0; 0 0]; M(:,:,2)=[-1 0; 1 1]; M(:,:,3)=[0 0; 0 -1.05]; q=1;
>> [Mf,Mb,kappa]=ilwhf(M,q);
>> Mf
Mf(:,:,1) =
   -0.0408         0
         0   -0.2300
Mf(:,:,2) =
    0.0389   -0.0110
    0.0000         0
>> Mb
Mb(:,:,1) =
   24.4898   -1.2245
   -4.3471   -4.3471
Mb(:,:,2) =
    0.0000    0.0000
    0.0000    4.5645
>> kappa
kappa =
     0     0
\end{verbatim}
Here $M_f(z)=\verb|Mf(:,:,1)|+\verb|Mf(:,:,2)|z^{-1}$ and $M_b(z)=\verb|Mb(:,:,1)|+\verb|Mb(:,:,2)|z$. These are exactly equal to the factor we obtained in Example 3.5 up to a non-singular transformation. To obtain the factors in Example 3.5 one simply computes $M_f(z)M_b(0)$ and $M_b^{-1}(0)M_b(z)$. \qed
\end{exmp}

According to Theorem \ref{thm:stability}, the set of non-singular elements of $M(z)\in\mathbb{R}^{n\times n}_{pq}(z)$ whose largest and smallest partial indices differ by more than 1 is non-generic. That is, the partial indices of such systems are sensitive to certain small perturbations and this will introduce numerical problems as we explore in the following example.

\begin{exmp}
Consider the Laurent matrix polynomial $M_\epsilon(z)=\left[\begin{smallmatrix} z & \epsilon\\ 0 & z^{-1}\end{smallmatrix}\right]$, which is discussed in Example 4.10 of the paper. The ILWHF of $M_\epsilon(z)$ relative to $\mathbb{T}$ has partial indices of $\{+1,-1\}$ for $\epsilon=0$ and $\{0,0\}$ for $\epsilon\neq0$. With tolerance set at machine epsilon, we obtain the following output:
\begin{verbatim}
>> M(:,:,1)=[0 0; 0 1]; M(:,:,2)=[0 10^-15; 0 0]; M(:,:,3)=[1 0; 0 0]; q=1;
>> [Mf,Mb,kappa]=ilwhf(M,q,eps);
>> kappa
kappa =
     0     0
>> M(:,:,1)=[0 0; 0 1]; M(:,:,2)=[0 10^-16; 0 0]; M(:,:,3)=[1 0; 0 0]; q=1;
>> [Mf,Mb,kappa]=ilwhf(M,q,eps);
>> kappa
kappa =
     1    -1
\end{verbatim}
Thus, the algorithm gives the correct partial indices for $\epsilon$ as small as $10^{-15}$ but provides incorrect partial indices when $\epsilon=10^{-16}$.\qed
\end{exmp}

Although the example above may be reassuring, the algorithm does fail in other situations where $M(z)$ is near the set of systems with $\kappa_1>\kappa_n+1$, which are exceptional by Theorem \ref{thm:stability}. In particular, the example given in \cite{so} cannot be computed with the algorithm above using a reasonable tolerance. Thus, the problem of formulating the optimal approach to computing the ILWHF must be left for future research.

It is worth emphasizing that even in the region where the algorithm is expected to work well (i.e.\ in the interior of the set of systems with $\kappa_1\leq\kappa_n+1$), it should be used only if the factors of the ILWHF are of interest. If the researcher is interested only in obtaining the representations (2) or (3) of the paper, then it is much quicker to use the \cite{sims} algorithm.

\setlength{\bibsep}{-2pt}
\bibliographystyle{newapa}
\bibliography{gire}

\end{document}