classical_filtering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9ffb6ea4",
   "metadata": {},
   "source": [
    "\n",
    "<a id='classical-filtering'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63bedda2",
   "metadata": {},
   "source": [
    "# Classical Prediction and Filtering With Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4da0d458",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This is a sequel to the earlier lecture [Classical Control with Linear Algebra](https://tools-techniques.quantecon.org/lu_tricks.html).\n",
    "\n",
    "That lecture used linear algebra – in particular,  the [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition)  – to formulate and solve a class of linear-quadratic optimal control problems.\n",
    "\n",
    "In this lecture, we’ll be using a closely related decomposition,\n",
    "the [Cholesky decomposition](https://en.wikipedia.org/wiki/Cholesky_decomposition), to solve linear prediction and filtering problems.\n",
    "\n",
    "We exploit the useful fact that there is an intimate connection between two superficially different classes of problems:\n",
    "\n",
    "- deterministic linear-quadratic (LQ) optimal control problems  \n",
    "- linear least squares prediction and filtering problems  \n",
    "\n",
    "\n",
    "The first class of problems involves no randomness, while the second is all about randomness.\n",
    "\n",
    "Nevertheless,  essentially the same mathematics  solves both types of problem.\n",
    "\n",
    "This connection, which is often termed “duality,” is present whether one uses “classical” or “recursive” solution procedures.\n",
    "\n",
    "In fact, we saw duality at work earlier when we formulated control and prediction problems recursively in lectures [LQ dynamic programming problems](https://python-intro.quantecon.org/lqcontrol.html), [A first look at the Kalman filter](https://python-intro.quantecon.org/kalman.html), and [The permanent income model](https://python-intro.quantecon.org/perm_income.html).\n",
    "\n",
    "A useful consequence of duality is that\n",
    "\n",
    "- With every LQ control problem, there is implicitly affiliated a linear least squares prediction or filtering problem.  \n",
    "- With every linear least squares prediction or filtering problem there is implicitly affiliated a LQ control problem.  \n",
    "\n",
    "\n",
    "An understanding of these connections has repeatedly proved useful in cracking interesting applied problems.\n",
    "\n",
    "For example, Sargent [[Sargent, 1987](https://tools-techniques.quantecon.org/zreferences.html#id197)] [chs. IX, XIV] and Hansen and Sargent [[Hansen and Sargent, 1980](https://tools-techniques.quantecon.org/zreferences.html#id107)] formulated\n",
    "and solved control and filtering problems using $ z $-transform methods.\n",
    "\n",
    "In this lecture, we begin to investigate these ideas by using mostly elementary linear algebra.\n",
    "\n",
    "This is the main purpose and focus of the lecture.\n",
    "\n",
    "However, after showing matrix algebra formulas, we’ll summarize classic infinite-horizon formulas built on $ z $-transform  and lag\n",
    "operator methods.\n",
    "\n",
    "And we’ll occasionally refer to some of these formulas from the infinite dimensional problems as we present the finite time\n",
    "formulas and associated linear algebra.\n",
    "\n",
    "We’ll start with the following standard import:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ae35b75",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d028611d",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "Useful references include [[Whittle, 1963](https://tools-techniques.quantecon.org/zreferences.html#id106)], [[Hansen and Sargent, 1980](https://tools-techniques.quantecon.org/zreferences.html#id107)], [[Orfanidis, 1988](https://tools-techniques.quantecon.org/zreferences.html#id108)], [[Athanasios and Pillai, 1991](https://tools-techniques.quantecon.org/zreferences.html#id109)], and [[Muth, 1960](https://tools-techniques.quantecon.org/zreferences.html#id110)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df08ea47",
   "metadata": {},
   "source": [
    "## Finite Dimensional Prediction\n",
    "\n",
    "Let $ (x_1, x_2, \\ldots, x_T)^\\prime = x $ be a $ T \\times 1 $ vector of random variables with mean $ \\mathbb{E} x = 0 $ and covariance matrix $ \\mathbb{E} xx^\\prime = V $.\n",
    "\n",
    "Here $ V $ is a $ T \\times T $ positive definite matrix.\n",
    "\n",
    "The $ i,j $ component $ E x_i x_j $ of $ V $ is the **inner product**   between $ x_i $ and $ x_j $.\n",
    "\n",
    "We regard the random variables as being\n",
    "ordered in time so that $ x_t $ is thought of as the value of some\n",
    "economic variable at time $ t $.\n",
    "\n",
    "For example, $ x_t $ could be generated by the random process described  by the Wold representation presented in equation [(6.16)](#equation-eq-31)\n",
    "in the section below on infinite dimensional prediction and filtering.\n",
    "\n",
    "In that case, $ V_{ij} $ is given by the coefficient on $ z^{\\mid i-j \\mid} $ in the expansion of $ g_x (z) = d(z) \\, d(z^{-1}) + h $, which equals\n",
    "$ h+\\sum^\\infty_{k=0} d_k d_{k+\\mid i-j \\mid} $.\n",
    "\n",
    "We want to  construct $ j $ step ahead linear least squares predictors of the form\n",
    "\n",
    "$$\n",
    "\\mathbb{\\hat E}\n",
    "\\left[\n",
    "    x_T\\vert x_{T-j}, x_{T-j + 1}, \\ldots, x_1\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "where $ \\mathbb{\\hat E} $ is the linear least squares projection operator.\n",
    "\n",
    "(Sometimes $ \\mathbb{\\hat E} $ is called the wide-sense expectations operator)\n",
    "\n",
    "To find linear least squares predictors it is helpful  first to construct a $ T \\times 1 $ vector $ \\varepsilon $\n",
    "of random variables that form an orthonormal basis   for the vector of random variables $ x $.\n",
    "\n",
    "The key insight here comes from noting that because the covariance matrix $ V $ of $ x $ is a positive definite and symmetric,\n",
    "there exists a (Cholesky) decomposition of $ V $ such that\n",
    "\n",
    "$$\n",
    "V = L^{-1} (L^{-1})^\\prime\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "L \\, V \\, L^\\prime = I\n",
    "$$\n",
    "\n",
    "where $ L $  and $ L^{-1} $ are both lower triangular.\n",
    "\n",
    "Form the $ T \\times 1 $ random vector $ \\varepsilon = Lx $.\n",
    "\n",
    "The random vector $ \\varepsilon $ is an orthonormal basis for $ x $ because\n",
    "\n",
    "- $ L $ is nonsingular  \n",
    "- $ \\mathbb{E} \\, \\varepsilon \\, \\varepsilon^\\prime = L \\mathbb{E} xx^\\prime L^\\prime = I $  \n",
    "- $ x = L^{-1} \\varepsilon $  \n",
    "\n",
    "\n",
    "It is enlightening  to write out and interpret the equations $ Lx = \\varepsilon $ and $ L^{-1} \\varepsilon = x $.\n",
    "\n",
    "First, we’ll write $ Lx = \\varepsilon $\n",
    "\n",
    "\n",
    "<a id='equation-eq-53'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    L_{11} x_1 &= \\varepsilon_1 \\\\\n",
    "    L_{21}x_1 + L_{22} x_2 &= \\varepsilon_2 \\\\ \\, \\vdots \\\\\n",
    "    L_{T1} \\, x_1 \\, \\ldots \\, + L_{TT} x_T &= \\varepsilon_T\n",
    "\\end{aligned} \\tag{6.1}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-eq-54'></a>\n",
    "$$\n",
    "\\sum^{t-1}_{j=0} L_{t,t-j}\\, x_{t-j} = \\varepsilon_t, \\quad t = 1, \\, 2, \\ldots T \\tag{6.2}\n",
    "$$\n",
    "\n",
    "Next, we write $ L^{-1} \\varepsilon = x $\n",
    "\n",
    "\n",
    "<a id='equation-eq-55a'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "     x_1 &= L_{11}^{-1} \\varepsilon_1 \\\\\n",
    "     x_2 &= L_{22}^{-1} \\varepsilon_2 + L_{21}^{-1} \\varepsilon_1  \\\\ \\, \\vdots \\\\\n",
    "     x_T &=  L_{TT}^{-1} \\varepsilon_T + L_{T,T-1}^{-1} \\varepsilon_{T-1} \\, \\ldots \\, + L_{T,1}^{-1} \\varepsilon_1\n",
    "\\end{aligned}, \\tag{6.3}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-eq-55'></a>\n",
    "$$\n",
    "x_t = \\sum^{t-1}_{j=0} L^{-1}_{t,t-j}\\, \\varepsilon_{t-j} \\tag{6.4}\n",
    "$$\n",
    "\n",
    "where $ L^{-1}_{i,j} $ denotes the $ i,j $ element of $ L^{-1} $.\n",
    "\n",
    "From [(6.2)](#equation-eq-54), it follows that $ \\varepsilon_t $ is in the linear subspace spanned by $ x_t,\\, x_{t-1}, \\ldots,\\, x_1 $.\n",
    "\n",
    "From [(6.4)](#equation-eq-55) it follows that  that $ x_t $ is in the linear subspace spanned by\n",
    "$ \\varepsilon_t, \\, \\varepsilon_{t-1}, \\ldots, \\varepsilon_1 $.\n",
    "\n",
    "Equation [(6.2)](#equation-eq-54) forms  a sequence of **autoregressions**  that for $ t = 1, \\ldots, T $ express\n",
    "$ x_t $  as linear functions of $ x_s, s = 1, \\ldots, t-1 $ and a random variable $ (L_{t,t})^{-1} \\varepsilon_t $\n",
    "that is orthogonal to each componenent of $ x_s, s = 1, \\ldots, t-1 $.\n",
    "\n",
    "(Here $ (L_{t,t})^{-1} $ denotes the reciprocal of $ L_{t,t} $ while $ L_{t,t}^{-1} $ denotes the $ t,t $ element\n",
    "of $ L^{-1} $).\n",
    "\n",
    "The equivalence of the subspaces spanned by $ \\varepsilon_t, \\ldots, \\varepsilon_1 $ and $ x_t, \\ldots, x_1 $ means that\n",
    "for $ t-1\\geq m \\geq 1 $\n",
    "\n",
    "\n",
    "<a id='equation-eq-56'></a>\n",
    "$$\n",
    "\\mathbb{\\hat E}\n",
    "[ x_t \\mid x_{t-m},\\, x_{t-m-1}, \\ldots, x_1 ] =\n",
    "\\mathbb{\\hat E}\n",
    "[x_t \\mid \\varepsilon_{t-m}, \\varepsilon_{t-m-1},\\ldots, \\varepsilon_1] \\tag{6.5}\n",
    "$$\n",
    "\n",
    "To proceed, it is useful to drill down and note that for $ t-1 \\geq m \\geq 1 $ we can  rewrite [(6.4)](#equation-eq-55)  in the form of the\n",
    "**moving average representation**\n",
    "\n",
    "\n",
    "<a id='equation-eq-57'></a>\n",
    "$$\n",
    "x_t = \\sum^{m-1}_{j=0} L_{t,t-j}^{-1}\\, \\varepsilon_{t-j} + \\sum^{t-1}_{j=m}\n",
    "L^{-1}_{t, t-j}\\, \\varepsilon_{t-j} \\tag{6.6}\n",
    "$$\n",
    "\n",
    "Representation [(6.6)](#equation-eq-57)  is an orthogonal decomposition of $ x_t $ into a part $ \\sum^{t-1}_{j=m} L_{t, t-j}^{-1}\\, \\varepsilon_{t-j} $\n",
    "that lies in the space spanned by $ [x_{t-m},\\, x_{t-m+1},\\, \\ldots, x_1] $ and an orthogonal component\n",
    "$ \\sum^{t-1}_{j=m} L^{-1}_{t, t-j}\\, \\varepsilon_{t-j} $ that does not lie in that space but instead in a linear space knowns as its **orthogonal complement**.\n",
    "\n",
    "It follows that\n",
    "\n",
    "$$\n",
    "\\mathbb{\\hat E}  [ x_t \\mid x_{t-m},\\, x_{t-m-1}, \\ldots, x_1 ] =  \\sum^{m-1}_{j=0} L_{t,t-j}^{-1}\\, \\varepsilon_{t-j}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97f260c6",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "\n",
    "Here’s the code that computes solutions to LQ control and filtering problems using the methods described here and in [Classical Control with Linear Algebra](https://tools-techniques.quantecon.org/lu_tricks.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a2c7fcd",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as spst\n",
    "import scipy.linalg as la\n",
    "\n",
    "class LQFilter:\n",
    "\n",
    "    def __init__(self, d, h, y_m, r=None, h_eps=None, β=None):\n",
    "        \"\"\"\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "            d : list or numpy.array (1-D or a 2-D column vector)\n",
    "                    The order of the coefficients: [d_0, d_1, ..., d_m]\n",
    "            h : scalar\n",
    "                    Parameter of the objective function (corresponding to the\n",
    "                    quadratic term)\n",
    "            y_m : list or numpy.array (1-D or a 2-D column vector)\n",
    "                    Initial conditions for y\n",
    "            r : list or numpy.array (1-D or a 2-D column vector)\n",
    "                    The order of the coefficients: [r_0, r_1, ..., r_k]\n",
    "                    (optional, if not defined -> deterministic problem)\n",
    "            β : scalar\n",
    "                    Discount factor (optional, default value is one)\n",
    "        \"\"\"\n",
    "\n",
    "        self.h = h\n",
    "        self.d = np.asarray(d)\n",
    "        self.m = self.d.shape[0] - 1\n",
    "\n",
    "        self.y_m = np.asarray(y_m)\n",
    "\n",
    "        if self.m == self.y_m.shape[0]:\n",
    "            self.y_m = self.y_m.reshape(self.m, 1)\n",
    "        else:\n",
    "            raise ValueError(\"y_m must be of length m = {self.m:d}\")\n",
    "\n",
    "        #---------------------------------------------\n",
    "        # Define the coefficients of ϕ upfront\n",
    "        #---------------------------------------------\n",
    "        ϕ = np.zeros(2 * self.m + 1)\n",
    "        for i in range(- self.m, self.m + 1):\n",
    "            ϕ[self.m - i] = np.sum(np.diag(self.d.reshape(self.m + 1, 1) \\\n",
    "                                           @ self.d.reshape(1, self.m + 1),\n",
    "                                           k=-i\n",
    "                                           )\n",
    "                                    )\n",
    "        ϕ[self.m] = ϕ[self.m] + self.h\n",
    "        self.ϕ = ϕ\n",
    "\n",
    "        #-----------------------------------------------------\n",
    "        # If r is given calculate the vector ϕ_r\n",
    "        #-----------------------------------------------------\n",
    "        if r is None:\n",
    "            pass\n",
    "        else:\n",
    "            self.r = np.asarray(r)\n",
    "            self.k = self.r.shape[0] - 1\n",
    "            ϕ_r = np.zeros(2 * self.k + 1)\n",
    "            for i in range(- self.k, self.k + 1):\n",
    "                ϕ_r[self.k - i] = np.sum(np.diag(self.r.reshape(self.k + 1, 1) \\\n",
    "                                                 @ self.r.reshape(1, self.k + 1),\n",
    "                                                 k=-i\n",
    "                                                 )\n",
    "                                        )\n",
    "            if h_eps is None:\n",
    "                self.ϕ_r = ϕ_r\n",
    "            else:\n",
    "                ϕ_r[self.k] = ϕ_r[self.k] + h_eps\n",
    "                self.ϕ_r = ϕ_r\n",
    "\n",
    "        #-----------------------------------------------------\n",
    "        # If β is given, define the transformed variables\n",
    "        #-----------------------------------------------------\n",
    "        if β is None:\n",
    "            self.β = 1\n",
    "        else:\n",
    "            self.β = β\n",
    "            self.d = self.β**(np.arange(self.m + 1)/2) * self.d\n",
    "            self.y_m = self.y_m * (self.β**(- np.arange(1, self.m + 1)/2)) \\\n",
    "                                   .reshape(self.m, 1)\n",
    "\n",
    "    def construct_W_and_Wm(self, N):\n",
    "        \"\"\"\n",
    "        This constructs the matrices W and W_m for a given number of periods N\n",
    "        \"\"\"\n",
    "\n",
    "        m = self.m\n",
    "        d = self.d\n",
    "\n",
    "        W = np.zeros((N + 1, N + 1))\n",
    "        W_m = np.zeros((N + 1, m))\n",
    "\n",
    "        #---------------------------------------\n",
    "        # Terminal conditions\n",
    "        #---------------------------------------\n",
    "\n",
    "        D_m1 = np.zeros((m + 1, m + 1))\n",
    "        M = np.zeros((m + 1, m))\n",
    "\n",
    "        # (1) Constuct the D_{m+1} matrix using the formula\n",
    "\n",
    "        for j in range(m + 1):\n",
    "            for k in range(j, m + 1):\n",
    "                D_m1[j, k] = d[:j + 1] @ d[k - j: k + 1]\n",
    "\n",
    "        # Make the matrix symmetric\n",
    "        D_m1 = D_m1 + D_m1.T - np.diag(np.diag(D_m1))\n",
    "\n",
    "        # (2) Construct the M matrix using the entries of D_m1\n",
    "\n",
    "        for j in range(m):\n",
    "            for i in range(j + 1, m + 1):\n",
    "                M[i, j] = D_m1[i - j - 1, m]\n",
    "\n",
    "        #----------------------------------------------\n",
    "        # Euler equations for t = 0, 1, ..., N-(m+1)\n",
    "        #----------------------------------------------\n",
    "        ϕ = self.ϕ\n",
    "\n",
    "        W[:(m + 1), :(m + 1)] = D_m1 + self.h * np.eye(m + 1)\n",
    "        W[:(m + 1), (m + 1):(2 * m + 1)] = M\n",
    "\n",
    "        for i, row in enumerate(np.arange(m + 1, N + 1 - m)):\n",
    "            W[row, (i + 1):(2 * m + 2 + i)] = ϕ\n",
    "\n",
    "        for i in range(1, m + 1):\n",
    "            W[N - m + i, -(2 * m + 1 - i):] = ϕ[:-i]\n",
    "\n",
    "        for i in range(m):\n",
    "            W_m[N - i, :(m - i)] = ϕ[(m + 1 + i):]\n",
    "\n",
    "        return W, W_m\n",
    "\n",
    "    def roots_of_characteristic(self):\n",
    "        \"\"\"\n",
    "        This function calculates z_0 and the 2m roots of the characteristic\n",
    "        equation associated with the Euler equation (1.7)\n",
    "\n",
    "        Note:\n",
    "        ------\n",
    "        numpy.poly1d(roots, True) defines a polynomial using its roots that can\n",
    "        be evaluated at any point. If x_1, x_2, ... , x_m are the roots then\n",
    "            p(x) = (x - x_1)(x - x_2)...(x - x_m)\n",
    "        \"\"\"\n",
    "        m = self.m\n",
    "        ϕ = self.ϕ\n",
    "\n",
    "        # Calculate the roots of the 2m-polynomial\n",
    "        roots = np.roots(ϕ)\n",
    "        # Sort the roots according to their length (in descending order)\n",
    "        roots_sorted = roots[np.argsort(abs(roots))[::-1]]\n",
    "\n",
    "        z_0 = ϕ.sum() / np.poly1d(roots, True)(1)\n",
    "        z_1_to_m = roots_sorted[:m]     # We need only those outside the unit circle\n",
    "\n",
    "        λ = 1 / z_1_to_m\n",
    "\n",
    "        return z_1_to_m, z_0, λ\n",
    "\n",
    "    def coeffs_of_c(self):\n",
    "        '''\n",
    "        This function computes the coefficients {c_j, j = 0, 1, ..., m} for\n",
    "                c(z) = sum_{j = 0}^{m} c_j z^j\n",
    "\n",
    "        Based on the expression (1.9). The order is\n",
    "            c_coeffs = [c_0, c_1, ..., c_{m-1}, c_m]\n",
    "        '''\n",
    "        z_1_to_m, z_0 = self.roots_of_characteristic()[:2]\n",
    "\n",
    "        c_0 = (z_0 * np.prod(z_1_to_m).real * (- 1)**self.m)**(.5)\n",
    "        c_coeffs = np.poly1d(z_1_to_m, True).c * z_0 / c_0\n",
    "\n",
    "        return c_coeffs[::-1]\n",
    "\n",
    "    def solution(self):\n",
    "        \"\"\"\n",
    "        This function calculates {λ_j, j=1,...,m} and {A_j, j=1,...,m}\n",
    "        of the expression (1.15)\n",
    "        \"\"\"\n",
    "        λ = self.roots_of_characteristic()[2]\n",
    "        c_0 = self.coeffs_of_c()[-1]\n",
    "\n",
    "        A = np.zeros(self.m, dtype=complex)\n",
    "        for j in range(self.m):\n",
    "            denom = 1 - λ/λ[j]\n",
    "            A[j] = c_0**(-2) / np.prod(denom[np.arange(self.m) != j])\n",
    "\n",
    "        return λ, A\n",
    "\n",
    "    def construct_V(self, N):\n",
    "        '''\n",
    "        This function constructs the covariance matrix for x^N (see section 6)\n",
    "        for a given period N\n",
    "        '''\n",
    "        V = np.zeros((N, N))\n",
    "        ϕ_r = self.ϕ_r\n",
    "\n",
    "        for i in range(N):\n",
    "            for j in range(N):\n",
    "                if abs(i-j) <= self.k:\n",
    "                    V[i, j] = ϕ_r[self.k + abs(i-j)]\n",
    "\n",
    "        return V\n",
    "\n",
    "    def simulate_a(self, N):\n",
    "        \"\"\"\n",
    "        Assuming that the u's are normal, this method draws a random path\n",
    "        for x^N\n",
    "        \"\"\"\n",
    "        V = self.construct_V(N + 1)\n",
    "        d = spst.multivariate_normal(np.zeros(N + 1), V)\n",
    "\n",
    "        return d.rvs()\n",
    "\n",
    "    def predict(self, a_hist, t):\n",
    "        \"\"\"\n",
    "        This function implements the prediction formula discussed in section 6 (1.59)\n",
    "        It takes a realization for a^N, and the period in which the prediction is \n",
    "        formed\n",
    "\n",
    "        Output:  E[abar | a_t, a_{t-1}, ..., a_1, a_0]\n",
    "        \"\"\"\n",
    "\n",
    "        N = np.asarray(a_hist).shape[0] - 1\n",
    "        a_hist = np.asarray(a_hist).reshape(N + 1, 1)\n",
    "        V = self.construct_V(N + 1)\n",
    "\n",
    "        aux_matrix = np.zeros((N + 1, N + 1))\n",
    "        aux_matrix[:(t + 1), :(t + 1)] = np.eye(t + 1)\n",
    "        L = la.cholesky(V).T\n",
    "        Ea_hist = la.inv(L) @ aux_matrix @ L @ a_hist\n",
    "\n",
    "        return Ea_hist\n",
    "\n",
    "    def optimal_y(self, a_hist, t=None):\n",
    "        \"\"\"\n",
    "        - if t is NOT given it takes a_hist (list or numpy.array) as a\n",
    "          deterministic a_t\n",
    "        - if t is given, it solves the combined control prediction problem \n",
    "          (section 7)(by default, t == None -> deterministic)\n",
    "\n",
    "        for a given sequence of a_t (either deterministic or a particular\n",
    "        realization), it calculates the optimal y_t sequence using the method\n",
    "        of the lecture\n",
    "\n",
    "        Note:\n",
    "        ------\n",
    "        scipy.linalg.lu normalizes L, U so that L has unit diagonal elements\n",
    "        To make things consistent with the lecture, we need an auxiliary\n",
    "        diagonal matrix D which renormalizes L and U\n",
    "        \"\"\"\n",
    "\n",
    "        N = np.asarray(a_hist).shape[0] - 1\n",
    "        W, W_m = self.construct_W_and_Wm(N)\n",
    "\n",
    "        L, U = la.lu(W, permute_l=True)\n",
    "        D = np.diag(1 / np.diag(U))\n",
    "        U = D @ U\n",
    "        L = L @ np.diag(1 / np.diag(D))\n",
    "\n",
    "        J = np.fliplr(np.eye(N + 1))\n",
    "\n",
    "        if t is None:   # If the problem is deterministic\n",
    "\n",
    "            a_hist = J @ np.asarray(a_hist).reshape(N + 1, 1)\n",
    "\n",
    "            #--------------------------------------------\n",
    "            # Transform the 'a' sequence if β is given\n",
    "            #--------------------------------------------\n",
    "            if self.β != 1:\n",
    "                a_hist =  a_hist * (self.β**(np.arange(N + 1) / 2))[::-1] \\\n",
    "                                    .reshape(N + 1, 1)\n",
    "\n",
    "            a_bar = a_hist - W_m @ self.y_m            # a_bar from the lecture\n",
    "            Uy = np.linalg.solve(L, a_bar)             # U @ y_bar = L^{-1}\n",
    "            y_bar = np.linalg.solve(U, Uy)             # y_bar = U^{-1}L^{-1}\n",
    "\n",
    "            # Reverse the order of y_bar with the matrix J\n",
    "            J = np.fliplr(np.eye(N + self.m + 1))\n",
    "            # y_hist : concatenated y_m and y_bar\n",
    "            y_hist = J @ np.vstack([y_bar, self.y_m])\n",
    "\n",
    "            #--------------------------------------------\n",
    "            # Transform the optimal sequence back if β is given\n",
    "            #--------------------------------------------\n",
    "            if self.β != 1:\n",
    "                y_hist = y_hist * (self.β**(- np.arange(-self.m, N + 1)/2)) \\\n",
    "                                    .reshape(N + 1 + self.m, 1)\n",
    "\n",
    "            return y_hist, L, U, y_bar\n",
    "\n",
    "        else:           # If the problem is stochastic and we look at it\n",
    "\n",
    "            Ea_hist = self.predict(a_hist, t).reshape(N + 1, 1)\n",
    "            Ea_hist = J @ Ea_hist\n",
    "\n",
    "            a_bar = Ea_hist - W_m @ self.y_m           # a_bar from the lecture\n",
    "            Uy = np.linalg.solve(L, a_bar)             # U @ y_bar = L^{-1}\n",
    "            y_bar = np.linalg.solve(U, Uy)             # y_bar = U^{-1}L^{-1}\n",
    "\n",
    "            # Reverse the order of y_bar with the matrix J\n",
    "            J = np.fliplr(np.eye(N + self.m + 1))\n",
    "            # y_hist : concatenated y_m and y_bar\n",
    "            y_hist = J @ np.vstack([y_bar, self.y_m])\n",
    "\n",
    "            return y_hist, L, U, y_bar"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf714cfc",
   "metadata": {},
   "source": [
    "Let’s use this code to tackle two interesting examples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a75c84ad",
   "metadata": {},
   "source": [
    "### Example 1\n",
    "\n",
    "Consider a stochastic process with moving average representation\n",
    "\n",
    "$$\n",
    "x_t = (1 - 2 L) \\varepsilon_t\n",
    "$$\n",
    "\n",
    "where $ \\varepsilon_t $ is a serially uncorrelated random process with mean zero and variance unity.\n",
    "\n",
    "If we were to use the tools associated with infinite dimensional prediction and filtering to be described below,\n",
    "we would use the Wiener-Kolmogorov formula [(6.21)](#equation-eq-36) to compute the linear least squares forecasts $ \\mathbb{E} [x_{t+j} \\mid x_t, x_{t-1}, \\ldots] $, for $ j = 1,\\, 2 $.\n",
    "\n",
    "But we can do everything we want by instead using our finite dimensional tools and\n",
    "setting $ d = r $, generating an instance of LQFilter, then invoking pertinent methods of LQFilter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45e13209",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "m = 1\n",
    "y_m = np.asarray([.0]).reshape(m, 1)\n",
    "d = np.asarray([1, -2])\n",
    "r = np.asarray([1, -2])\n",
    "h = 0.0\n",
    "example = LQFilter(d, h, y_m, r=d)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4ab3b6e",
   "metadata": {},
   "source": [
    "The Wold representation is computed by `example.coeffs_of_c()`.\n",
    "\n",
    "Let’s check that it “flips roots” as required"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afeee280",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "example.coeffs_of_c()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "883d7066",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "example.roots_of_characteristic()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "151f980e",
   "metadata": {},
   "source": [
    "Now let’s form the covariance matrix of a time series vector of length $ N $\n",
    "and put it in $ V $.\n",
    "\n",
    "Then we’ll take a Cholesky decomposition of $ V = L^{-1} L^{-1} $ and use it to form the vector of\n",
    "“moving average representations” $ x = L^{-1} \\varepsilon $ and the vector of “autoregressive representations” $ L x = \\varepsilon $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8a9d67c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "V = example.construct_V(N=5)\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c80b46d6",
   "metadata": {},
   "source": [
    "Notice how the lower rows of the “moving average representations” are converging to the appropriate infinite history Wold representation\n",
    "to be described below when we study infinite horizon-prediction and filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4560d24",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Li = np.linalg.cholesky(V)\n",
    "print(Li)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d5bee4e",
   "metadata": {},
   "source": [
    "Notice how the lower rows of the “autoregressive representations” are converging to the appropriate infinite-history\n",
    "autoregressive representation to be described below when we study infinite horizon-prediction and filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0de1eb2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "L = np.linalg.inv(Li)\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5d6f43",
   "metadata": {},
   "source": [
    "### Example 2\n",
    "\n",
    "Consider a stochastic process $ X_t $ with moving average\n",
    "representation\n",
    "\n",
    "$$\n",
    "X_t = (1 - \\sqrt 2 L^2) \\varepsilon_t\n",
    "$$\n",
    "\n",
    "where $ \\varepsilon_t $ is a serially uncorrelated random process\n",
    "with mean zero and variance unity.\n",
    "\n",
    "Let’s find a Wold moving average representation for $ x_t $ that will prevail in the infinite-history context to be studied in\n",
    "detail below.\n",
    "\n",
    "To do this, we’ll  use the Wiener-Kolomogorov formula [(6.21)](#equation-eq-36) presented below to compute the linear least squares forecasts\n",
    "$ \\mathbb{\\hat E}\\left[X_{t+j} \\mid X_{t-1}, \\ldots\\right] \\hbox { for } j = 1,\\, 2,\\, 3 $.\n",
    "\n",
    "We proceed in the same way as in example 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48cab849",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "m = 2\n",
    "y_m = np.asarray([.0, .0]).reshape(m, 1)\n",
    "d = np.asarray([1, 0, -np.sqrt(2)])\n",
    "r = np.asarray([1, 0, -np.sqrt(2)])\n",
    "h = 0.0\n",
    "example = LQFilter(d, h, y_m, r=d)\n",
    "example.coeffs_of_c()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73c846eb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "example.roots_of_characteristic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d8c497c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "V = example.construct_V(N=8)\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31448785",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Li = np.linalg.cholesky(V)\n",
    "print(Li[-3:, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbfb538d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "L = np.linalg.inv(Li)\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7fb5a7a",
   "metadata": {},
   "source": [
    "### Prediction\n",
    "\n",
    "It immediately follows from the “orthogonality principle” of least squares (see [[Athanasios and Pillai, 1991](https://tools-techniques.quantecon.org/zreferences.html#id109)] or [[Sargent, 1987](https://tools-techniques.quantecon.org/zreferences.html#id197)] [ch. X]) that\n",
    "\n",
    "\n",
    "<a id='equation-eq-58'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\mathbb{\\hat E} & [x_t \\mid x_{t-m},\\, x_{t-m+1}, \\ldots x_1]\n",
    "                    = \\sum^{t-1}_{j=m} L^{-1}_{t,t-j}\\, \\varepsilon_{t-j} \\\\\n",
    "               & = [L_{t, 1}^{-1}\\, L^{-1}_{t,2},\\, \\ldots, L^{-1}_{t,t-m}\\ 0 \\ 0 \\ldots 0] L \\, x\n",
    "\\end{aligned} \\tag{6.7}\n",
    "$$\n",
    "\n",
    "This can be interpreted as a finite-dimensional version of the Wiener-Kolmogorov $ m $-step ahead prediction formula.\n",
    "\n",
    "We can use [(6.7)](#equation-eq-58)  to represent the linear least squares projection of\n",
    "the vector $ x $ conditioned on the first $ s $ observations\n",
    "$ [x_s, x_{s-1} \\ldots, x_1] $.\n",
    "\n",
    "We have\n",
    "\n",
    "\n",
    "<a id='equation-eq-59'></a>\n",
    "$$\n",
    "\\mathbb{\\hat E}[x \\mid x_s, x_{s-1}, \\ldots, x_1]\n",
    "= L^{-1}\n",
    "\\left[\n",
    "    \\begin{matrix}\n",
    "        I_s & 0 \\\\\n",
    "        0 & 0_{(t-s)}\n",
    "    \\end{matrix}\n",
    "\\right] L x \\tag{6.8}\n",
    "$$\n",
    "\n",
    "This formula will be convenient in representing the solution of control problems under uncertainty.\n",
    "\n",
    "Equation [(6.4)](#equation-eq-55)  can be recognized as a finite dimensional version of a moving average representation.\n",
    "\n",
    "Equation  [(6.2)](#equation-eq-54) can be viewed as a finite dimension version of an autoregressive representation.\n",
    "\n",
    "Notice that even\n",
    "if the $ x_t $ process is covariance stationary, so that $ V $\n",
    "is such that $ V_{ij} $ depends only on $ \\vert i-j\\vert $, the\n",
    "coefficients in the moving average representation are time-dependent,\n",
    "there being a different moving average for each $ t $.\n",
    "\n",
    "If\n",
    "$ x_t $ is a covariance stationary process, the last row of\n",
    "$ L^{-1} $ converges to the coefficients in the Wold moving average\n",
    "representation for $ \\{ x_t\\} $ as $ T \\rightarrow \\infty $.\n",
    "\n",
    "Further, if $ x_t $ is covariance stationary, for fixed $ k $\n",
    "and $ j > 0, \\, L^{-1}_{T,T-j} $ converges to\n",
    "$ L^{-1}_{T-k, T-k-j} $ as $ T \\rightarrow \\infty $.\n",
    "\n",
    "That is,\n",
    "the “bottom” rows of $ L^{-1} $ converge to each other and to the\n",
    "Wold moving average coefficients as $ T \\rightarrow \\infty $.\n",
    "\n",
    "This last observation gives one simple and widely-used practical way of\n",
    "forming a finite $ T $ approximation to a Wold moving average\n",
    "representation.\n",
    "\n",
    "First, form the covariance matrix\n",
    "$ \\mathbb{E}xx^\\prime = V $, then obtain the Cholesky decomposition\n",
    "$ L^{-1} L^{-1^\\prime} $ of $ V $, which can be accomplished\n",
    "quickly on a computer.\n",
    "\n",
    "The last row of $ L^{-1} $ gives the approximate Wold moving average coefficients.\n",
    "\n",
    "This method can readily be generalized to multivariate systems.\n",
    "\n",
    "\n",
    "<a id='fdcp'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94af7e11",
   "metadata": {},
   "source": [
    "## Combined Finite Dimensional Control and Prediction\n",
    "\n",
    "Consider the finite-dimensional control problem, maximize\n",
    "\n",
    "$$\n",
    "\\mathbb{E} \\, \\sum^N_{t=0} \\,\n",
    "\\left\\{\n",
    "     a_t y_t - {1 \\over 2} h y^2_t - {1 \\over 2} [d(L) y_t ]^2\n",
    "\\right\\},\\  \\quad h > 0\n",
    "$$\n",
    "\n",
    "where $ d(L) = d_0 + d_1 L+ \\ldots + d_m L^m $, $ L $ is the\n",
    "lag operator, $ \\bar a = [ a_N, a_{N-1} \\ldots, a_1, a_0]^\\prime $ a\n",
    "random vector with mean zero and $ \\mathbb{E}\\,\\bar a \\bar a^\\prime = V $.\n",
    "\n",
    "The variables $ y_{-1}, \\ldots, y_{-m} $ are given.\n",
    "\n",
    "Maximization is over choices of $ y_0,\n",
    "y_1 \\ldots, y_N $, where $ y_t $ is required to be a linear function\n",
    "of $ \\{y_{t-s-1}, t+m-1\\geq 0;\\ a_{t-s}, t\\geq s\\geq 0\\} $.\n",
    "\n",
    "We saw in the lecture [Classical Control with Linear Algebra](https://tools-techniques.quantecon.org/lu_tricks.html)  that the solution of this problem under certainty could be represented in the feedback-feedforward form\n",
    "\n",
    "$$\n",
    "U \\bar y\n",
    "   = L^{-1}\\bar a + K\n",
    "   \\left[\n",
    "     \\begin{matrix}\n",
    "         y_{-1}\\\\\n",
    "         \\vdots\\\\\n",
    "         y_{-m}\n",
    "     \\end{matrix}\n",
    "   \\right]\n",
    "$$\n",
    "\n",
    "for some $ (N+1)\\times m $ matrix $ K $.\n",
    "\n",
    "Using a version of formula [(6.7)](#equation-eq-58), we can express $ \\mathbb{\\hat E}[\\bar a \\mid a_s,\\, a_{s-1}, \\ldots, a_0 ] $ as\n",
    "\n",
    "$$\n",
    "\\mathbb{\\hat E}\n",
    "[ \\bar a \\mid a_s,\\, a_{s-1}, \\ldots, a_0]\n",
    "= \\tilde U^{-1}\n",
    "\\left[\n",
    "    \\begin{matrix}\n",
    "        0 & 0 \\\\\n",
    "        0 & I_{(s+1)}\n",
    "    \\end{matrix}\n",
    "\\right]\n",
    "\\tilde U \\bar a\n",
    "$$\n",
    "\n",
    "where $ I_{(s + 1)} $ is the $ (s+1) \\times (s+1) $ identity\n",
    "matrix, and $ V = \\tilde U^{-1} \\tilde U^{-1^{\\prime}} $, where\n",
    "$ \\tilde U $ is the *upper* triangular Cholesky factor of the\n",
    "covariance matrix $ V $.\n",
    "\n",
    "(We have reversed the time axis in dating the $ a $’s relative to earlier)\n",
    "\n",
    "The time axis can be reversed in representation [(6.8)](#equation-eq-59) by replacing $ L $ with $ L^T $.\n",
    "\n",
    "The optimal decision rule to use at time $ 0 \\leq t \\leq N $ is then\n",
    "given by the $ (N-t +1)^{\\rm th} $ row of\n",
    "\n",
    "$$\n",
    "U \\bar y = L^{-1} \\tilde U^{-1}\n",
    "    \\left[\n",
    "        \\begin{matrix}\n",
    "            0 & 0 \\\\\n",
    "            0 & I_{(t+1)}\n",
    "        \\end{matrix}\n",
    "    \\right]\n",
    "    \\tilde U \\bar a + K\n",
    "    \\left[\n",
    "    \\begin{matrix}\n",
    "        y_{-1}\\\\\n",
    "        \\vdots\\\\\n",
    "        y_{-m}\n",
    "    \\end{matrix}\n",
    "    \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e791356",
   "metadata": {},
   "source": [
    "## Infinite Horizon Prediction and Filtering Problems\n",
    "\n",
    "It is instructive to compare the finite-horizon formulas based on linear algebra decompositions of finite-dimensional covariance matrices\n",
    "with classic formulas for infinite horizon and infinite history prediction and control problems.\n",
    "\n",
    "These classic infinite horizon formulas used the mathematics of $ z $-transforms and lag operators.\n",
    "\n",
    "We’ll meet interesting lag operator and $ z $-transform  counterparts to our finite horizon matrix formulas.\n",
    "\n",
    "We pose two related prediction and filtering problems.\n",
    "\n",
    "We let $ Y_t $ be a univariate $ m^{\\rm th} $ order moving average, covariance stationary stochastic process,\n",
    "\n",
    "\n",
    "<a id='equation-eq-24'></a>\n",
    "$$\n",
    "Y_t = d(L) u_t \\tag{6.9}\n",
    "$$\n",
    "\n",
    "where $ d(L) = \\sum^m_{j=0} d_j L^j $, and $ u_t $ is a serially uncorrelated stationary random process satisfying\n",
    "\n",
    "\n",
    "<a id='equation-eq-25'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\mathbb{E} u_t &= 0\\\\\n",
    "    \\mathbb{E} u_t u_s &=\n",
    "    \\begin{cases}\n",
    "        1 & \\text{ if } t = s \\\\\n",
    "        0 & \\text{ otherwise}\n",
    "    \\end{cases}\n",
    "\\end{aligned} \\tag{6.10}\n",
    "$$\n",
    "\n",
    "We impose no conditions on the zeros of $ d(z) $.\n",
    "\n",
    "A second covariance stationary process is $ X_t $ given by\n",
    "\n",
    "\n",
    "<a id='equation-eq-26'></a>\n",
    "$$\n",
    "X_t = Y_t + \\varepsilon_t \\tag{6.11}\n",
    "$$\n",
    "\n",
    "where $ \\varepsilon_t $ is a serially uncorrelated stationary\n",
    "random process with $ \\mathbb{E} \\varepsilon_t = 0 $ and $ \\mathbb{E} \\varepsilon_t \\varepsilon_s $ = $ 0 $ for all distinct $ t $ and $ s $.\n",
    "\n",
    "We also assume that $ \\mathbb{E} \\varepsilon_t u_s = 0 $ for all $ t $ and $ s $.\n",
    "\n",
    "The **linear least squares prediction problem** is to find the $ L_2 $\n",
    "random variable $ \\hat X_{t+j} $ among linear combinations of\n",
    "$ \\{ X_t,\\  X_{t-1},\n",
    "\\ldots \\} $ that minimizes $ \\mathbb{E}(\\hat X_{t+j} - X_{t+j})^2 $.\n",
    "\n",
    "That is, the problem is to find a $ \\gamma_j (L) = \\sum^\\infty_{k=0} \\gamma_{jk}\\, L^k $ such that $ \\sum^\\infty_{k=0} \\vert \\gamma_{jk} \\vert^2 < \\infty $ and $ \\mathbb{E} [\\gamma_j \\, (L) X_t -X_{t+j}]^2 $ is minimized.\n",
    "\n",
    "The **linear least squares filtering problem** is to find a $ b\\,(L) = \\sum^\\infty_{j=0} b_j\\, L^j $ such that $ \\sum^\\infty_{j=0}\\vert b_j \\vert^2 < \\infty $ and $ \\mathbb{E} [b\\, (L) X_t -Y_t ]^2 $ is minimized.\n",
    "\n",
    "Interesting versions of these problems related to the permanent income theory were studied by [[Muth, 1960](https://tools-techniques.quantecon.org/zreferences.html#id110)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "722e0b5c",
   "metadata": {},
   "source": [
    "### Problem Formulation\n",
    "\n",
    "These problems are solved as follows.\n",
    "\n",
    "The covariograms of $ Y $ and $ X $ and their cross covariogram are, respectively,\n",
    "\n",
    "\n",
    "<a id='equation-eq-27'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    C_X (\\tau) &= \\mathbb{E}X_t X_{t-\\tau} \\\\\n",
    "    C_Y (\\tau) &= \\mathbb{E}Y_t Y_{t-\\tau}  \\qquad \\tau = 0, \\pm 1, \\pm 2, \\ldots \\\\\n",
    "    C_{Y,X} (\\tau) &= \\mathbb{E}Y_t X_{t-\\tau}\n",
    "\\end{aligned} \\tag{6.12}\n",
    "$$\n",
    "\n",
    "The covariance and cross-covariance generating functions are defined as\n",
    "\n",
    "\n",
    "<a id='equation-eq-28'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    g_X(z) &= \\sum^\\infty_{\\tau = - \\infty} C_X (\\tau) z^\\tau \\\\\n",
    "    g_Y(z) &= \\sum^\\infty_{\\tau = - \\infty} C_Y (\\tau) z^\\tau \\\\\n",
    "    g_{YX} (z) &= \\sum^\\infty_{\\tau = - \\infty} C_{YX} (\\tau) z^\\tau\n",
    "\\end{aligned} \\tag{6.13}\n",
    "$$\n",
    "\n",
    "The generating functions can be computed by using the following facts.\n",
    "\n",
    "Let $ v_{1t} $ and $ v_{2t} $ be two mutually and serially uncorrelated white noises with unit variances.\n",
    "\n",
    "That is, $ \\mathbb{E}v^2_{1t} = \\mathbb{E}v^2_{2t} = 1, \\mathbb{E}v_{1t} = \\mathbb{E}v_{2t} = 0, \\mathbb{E}v_{1t} v_{2s} = 0 $ for all $ t $ and $ s $, $ \\mathbb{E}v_{1t} v_{1t-j} = \\mathbb{E}v_{2t} v_{2t-j} = 0 $ for all $ j \\not = 0 $.\n",
    "\n",
    "Let $ x_t $ and $ y_t $ be two random processes given by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    y_t &= A(L) v_{1t} + B(L) v_{2t} \\\\\n",
    "    x_t &= C(L) v_{1t} + D(L) v_{2t}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Then, as shown for example in [[Sargent, 1987](https://tools-techniques.quantecon.org/zreferences.html#id197)] [ch. XI], it is true that\n",
    "\n",
    "\n",
    "<a id='equation-eq-29'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    g_y(z) &= A(z) A(z^{-1}) + B (z) B(z^{-1}) \\\\\n",
    "    g_x (z) &= C(z) C(z^{-1}) + D(z) D(z^{-1}) \\\\\n",
    "    g_{yx} (z) &= A(z) C(z^{-1}) + B(z) D(z^{-1})\n",
    "\\end{aligned} \\tag{6.14}\n",
    "$$\n",
    "\n",
    "Applying these formulas to [(6.9)](#equation-eq-24) – [(6.12)](#equation-eq-27), we have\n",
    "\n",
    "\n",
    "<a id='equation-eq-30'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    g_Y(z) &= d(z)d(z^{-1}) \\\\\n",
    "    g_X(z) &= d(z)d(z^{-1}) + h\\\\\n",
    "    g_{YX} (z) &= d(z) d(z^{-1})\n",
    "\\end{aligned} \\tag{6.15}\n",
    "$$\n",
    "\n",
    "The key step in obtaining solutions to our problems is to factor the covariance generating function  $ g_X(z) $ of $ X $.\n",
    "\n",
    "The solutions of our problems are given by formulas due to Wiener and Kolmogorov.\n",
    "\n",
    "These formulas utilize the Wold moving average representation of the $ X_t $ process,\n",
    "\n",
    "\n",
    "<a id='equation-eq-31'></a>\n",
    "$$\n",
    "X_t = c\\,(L)\\,\\eta_t \\tag{6.16}\n",
    "$$\n",
    "\n",
    "where $ c(L) = \\sum^m_{j=0} c_j\\, L^j $, with\n",
    "\n",
    "\n",
    "<a id='equation-eq-32'></a>\n",
    "$$\n",
    "c_0 \\eta_t\n",
    "= X_t - \\mathbb{\\hat E} [X_t \\vert X_{t-1}, X_{t-2}, \\ldots] \\tag{6.17}\n",
    "$$\n",
    "\n",
    "Here $ \\mathbb{\\hat E} $ is the linear least squares projection operator.\n",
    "\n",
    "Equation [(6.17)](#equation-eq-32)  is the condition that $ c_0 \\eta_t $ can be the one-step-ahead error in predicting $ X_t $ from its own past values.\n",
    "\n",
    "Condition [(6.17)](#equation-eq-32)  requires that $ \\eta_t $ lie in the closed\n",
    "linear space spanned by $ [X_t,\\  X_{t-1}, \\ldots] $.\n",
    "\n",
    "This will be true if and only if the zeros of $ c(z) $ do not lie inside the unit circle.\n",
    "\n",
    "It is an implication of [(6.17)](#equation-eq-32) that $ \\eta_t $ is a serially\n",
    "uncorrelated random process and that normalization can be imposed so\n",
    "that $ \\mathbb{E}\\eta_t^2 = 1 $.\n",
    "\n",
    "Consequently, an implication of [(6.16)](#equation-eq-31)  is\n",
    "that the covariance generating function of $ X_t $ can be expressed\n",
    "as\n",
    "\n",
    "\n",
    "<a id='equation-eq-33'></a>\n",
    "$$\n",
    "g_X(z) = c\\,(z)\\,c\\,(z^{-1}) \\tag{6.18}\n",
    "$$\n",
    "\n",
    "It remains to discuss how $ c(L) $ is to be computed.\n",
    "\n",
    "Combining [(6.14)](#equation-eq-29)  and [(6.18)](#equation-eq-33)  gives\n",
    "\n",
    "\n",
    "<a id='equation-eq-34'></a>\n",
    "$$\n",
    "d(z) \\,d(z^{-1}) + h = c \\, (z) \\,c\\,(z^{-1}) \\tag{6.19}\n",
    "$$\n",
    "\n",
    "Therefore, we have already shown constructively how to factor the covariance generating function $ g_X(z) = d(z)\\,d\\,(z^{-1}) + h $.\n",
    "\n",
    "We now introduce the **annihilation operator**:\n",
    "\n",
    "\n",
    "<a id='equation-eq-35'></a>\n",
    "$$\n",
    "\\left[\n",
    "    \\sum^\\infty_{j= - \\infty} f_j\\, L^j\n",
    "\\right]_+\n",
    "\\equiv \\sum^\\infty_{j=0} f_j\\,L^j \\tag{6.20}\n",
    "$$\n",
    "\n",
    "In words, $ [\\phantom{00}]_+ $ means “ignore negative powers of $ L $”.\n",
    "\n",
    "We have defined the solution of the prediction problem as $ \\mathbb{\\hat E} [X_{t+j} \\vert X_t,\\, X_{t-1}, \\ldots] = \\gamma_j\\, (L) X_t $.\n",
    "\n",
    "Assuming that the roots of $ c(z) = 0 $ all lie outside the unit circle, the Wiener-Kolmogorov formula for $ \\gamma_j (L) $ holds:\n",
    "\n",
    "\n",
    "<a id='equation-eq-36'></a>\n",
    "$$\n",
    "\\gamma_j\\, (L) =\n",
    "\\left[\n",
    "    {c (L) \\over L^j}\n",
    "\\right]_+ c\\,(L)^{-1} \\tag{6.21}\n",
    "$$\n",
    "\n",
    "We have defined the solution of the filtering problem as $ \\mathbb{\\hat E}[Y_t \\mid X_t, X_{t-1}, \\ldots] = b (L)X_t $.\n",
    "\n",
    "The Wiener-Kolomogorov formula for $ b(L) $ is\n",
    "\n",
    "$$\n",
    "b(L) = \\left[{g_{YX} (L) \\over c(L^{-1})}\\right]_+ c(L)^{-1}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-eq-37'></a>\n",
    "$$\n",
    "b(L) = \\left[ {d(L)d(L^{-1}) \\over c(L^{-1})} \\right]_+ c(L)^{-1} \\tag{6.22}\n",
    "$$\n",
    "\n",
    "Formulas [(6.21)](#equation-eq-36) and [(6.22)](#equation-eq-37)  are discussed in detail in  [[Whittle, 1983](https://tools-techniques.quantecon.org/zreferences.html#id85)] and [[Sargent, 1987](https://tools-techniques.quantecon.org/zreferences.html#id197)].\n",
    "\n",
    "The interested reader can there find several examples of the use of these formulas in economics\n",
    "Some classic examples using these formulas are due to [[Muth, 1960](https://tools-techniques.quantecon.org/zreferences.html#id110)].\n",
    "\n",
    "As an example of the usefulness of formula [(6.22)](#equation-eq-37), we let $ X_t $ be a stochastic process with Wold moving average representation\n",
    "\n",
    "$$\n",
    "X_t = c(L) \\eta_t\n",
    "$$\n",
    "\n",
    "where $ \\mathbb{E}\\eta^2_t = 1, \\hbox { and } c_0 \\eta_t = X_t - \\mathbb{\\hat E} [X_t \\vert X_{t-1}, \\ldots], c (L) = \\sum^m_{j=0} c_j L $.\n",
    "\n",
    "Suppose that at time $ t $, we wish to predict a geometric sum of future $ X $’s, namely\n",
    "\n",
    "$$\n",
    "y_t \\equiv \\sum^\\infty_{j=0} \\delta^j X_{t+j} = {1 \\over 1 - \\delta L^{-1}}\n",
    "X_t\n",
    "$$\n",
    "\n",
    "given knowledge of $ X_t, X_{t-1}, \\ldots $.\n",
    "\n",
    "We shall use [(6.22)](#equation-eq-37)  to obtain the answer.\n",
    "\n",
    "Using the standard formulas  [(6.14)](#equation-eq-29), we have that\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    g_{yx}(z) &= (1-\\delta z^{-1})c(z) c (z^{-1}) \\\\\n",
    "    g_x (z) &= c(z) c (z^{-1})\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Then [(6.22)](#equation-eq-37)  becomes\n",
    "\n",
    "\n",
    "<a id='equation-eq-38'></a>\n",
    "$$\n",
    "b(L)=\\left[{c(L)\\over 1-\\delta L^{-1}}\\right]_+ c(L)^{-1} \\tag{6.23}\n",
    "$$\n",
    "\n",
    "In order to evaluate the term in the annihilation operator, we use the following result from [[Hansen and Sargent, 1980](https://tools-techniques.quantecon.org/zreferences.html#id107)].\n",
    "\n",
    "**Proposition** Let\n",
    "\n",
    "- $ g(z) = \\sum^\\infty_{j=0} g_j \\, z^j $ where $ \\sum^\\infty_{j=0} \\vert g_j \\vert^2 < + \\infty $.  \n",
    "- $ h\\,(z^{-1}) = $ $ (1- \\delta_1 z^{-1}) \\ldots (1-\\delta_n z^{-1}) $, where $ \\vert \\delta_j \\vert < 1 $, for $ j = 1, \\ldots, n $.  \n",
    "\n",
    "\n",
    "Then\n",
    "\n",
    "\n",
    "<a id='equation-eq-39'></a>\n",
    "$$\n",
    "\\left[{g(z)\\over h(z^{-1})}\\right]_+ = {g(z)\\over h(z^{-1})} - \\sum^n_{j=1}\n",
    "\\ {\\delta_j g (\\delta_j) \\over \\prod^n_{k=1 \\atop k \\not = j} (\\delta_j -\n",
    "\\delta_k)} \\ \\left({1 \\over z- \\delta_j}\\right) \\tag{6.24}\n",
    "$$\n",
    "\n",
    "and, alternatively,\n",
    "\n",
    "\n",
    "<a id='equation-eq-40'></a>\n",
    "$$\n",
    "\\left[\n",
    "    {g(z)\\over h(z^{-1})}\n",
    "\\right]_+\n",
    "=\\sum^n_{j=1} B_j\n",
    "\\left(\n",
    "    {zg(z)-\\delta_j g (\\delta_j) \\over z- \\delta_j}\n",
    "\\right) \\tag{6.25}\n",
    "$$\n",
    "\n",
    "where $ B_j = 1 / \\prod^n_{k=1\\atop k+j} (1 - \\delta_k / \\delta_j) $.\n",
    "\n",
    "Applying formula [(6.25)](#equation-eq-40)  of the proposition to evaluating  [(6.23)](#equation-eq-38)  with $ g(z) = c(z) $ and $ h(z^{-1}) = 1 - \\delta z^{-1} $ gives\n",
    "\n",
    "$$\n",
    "b(L)=\\left[{Lc(L)-\\delta c(\\delta)\\over L-\\delta}\\right] c(L)^{-1}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "b(L) =\n",
    "\\left[\n",
    "    {1-\\delta c (\\delta) L^{-1} c (L)^{-1}\\over 1-\\delta L^{-1}}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "Thus, we have\n",
    "\n",
    "\n",
    "<a id='equation-eq-41'></a>\n",
    "$$\n",
    "\\mathbb{\\hat E}\n",
    "\\left[\n",
    "    \\sum^\\infty_{j=0} \\delta^j X_{t+j}\\vert X_t,\\, x_{t-1},\n",
    "    \\ldots\n",
    "\\right] =\n",
    "\\left[\n",
    "    {1-\\delta c (\\delta) L^{-1} c(L)^{-1} \\over 1 - \\delta L^{-1}}\n",
    "\\right]\n",
    "\\, X_t \\tag{6.26}\n",
    "$$\n",
    "\n",
    "This formula is useful in solving stochastic versions of problem 1 of lecture [Classical Control with Linear Algebra](https://tools-techniques.quantecon.org/lu_tricks.html) in which the randomness emerges because $ \\{a_t\\} $ is a stochastic\n",
    "process.\n",
    "\n",
    "The problem is to maximize\n",
    "\n",
    "\n",
    "<a id='equation-eq-42'></a>\n",
    "$$\n",
    "\\mathbb{E}_0\n",
    "\\lim_{N \\rightarrow \\infty}\\\n",
    "\\sum^N_{t-0} \\beta^t\n",
    "\\left[\n",
    "    a_t\\, y_t - {1 \\over 2}\\ hy^2_t-{1 \\over 2}\\ [d(L)y_t]^2\n",
    "\\right] \\tag{6.27}\n",
    "$$\n",
    "\n",
    "where $ \\mathbb{E}_t $ is mathematical expectation conditioned on information\n",
    "known at $ t $, and where $ \\{ a_t\\} $ is a covariance\n",
    "stationary stochastic process with Wold moving average representation\n",
    "\n",
    "$$\n",
    "a_t = c(L)\\, \\eta_t\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "c(L) = \\sum^{\\tilde n}_{j=0} c_j L^j\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\eta_t =\n",
    "a_t - \\mathbb{\\hat E} [a_t \\vert a_{t-1}, \\ldots]\n",
    "$$\n",
    "\n",
    "The problem is to maximize [(6.27)](#equation-eq-42)  with respect to a contingency plan\n",
    "expressing $ y_t $ as a function of information known at $ t $,\n",
    "which is assumed to be\n",
    "$ (y_{t-1},\\  y_{t-2}, \\ldots, a_t, \\ a_{t-1}, \\ldots) $.\n",
    "\n",
    "The solution of this problem can be achieved in two steps.\n",
    "\n",
    "First, ignoring the uncertainty, we can solve the problem assuming that $ \\{a_t\\} $ is a known sequence.\n",
    "\n",
    "The solution is, from above,\n",
    "\n",
    "$$\n",
    "c(L) y_t = c(\\beta L^{-1})^{-1} a_t\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-eq-43'></a>\n",
    "$$\n",
    "(1-\\lambda_1 L) \\ldots (1 - \\lambda_m L) y_t\n",
    "= \\sum^m_{j=1} A_j\n",
    "\\sum^\\infty_{k=0} (\\lambda_j \\beta)^k\\, a_{t+k} \\tag{6.28}\n",
    "$$\n",
    "\n",
    "Second, the solution of the problem under uncertainty is obtained by\n",
    "replacing the terms on the right-hand side of the above expressions with\n",
    "their linear least squares predictors.\n",
    "\n",
    "Using [(6.26)](#equation-eq-41) and [(6.28)](#equation-eq-43), we have\n",
    "the following solution\n",
    "\n",
    "$$\n",
    "(1-\\lambda_1 L) \\ldots (1-\\lambda_m L) y_t =\n",
    "\\sum^m_{j=1} A_j\n",
    " \\left[\n",
    "     \\frac{1-\\beta \\lambda_j \\, c (\\beta \\lambda_j) L^{-1} c(L)^{-1} }\n",
    "     { 1-\\beta \\lambda_j L^{-1} }\n",
    " \\right] a_t\n",
    "$$\n",
    "\n",
    "**Blaschke factors**\n",
    "\n",
    "The following is a useful piece of mathematics underlying “root flipping”.\n",
    "\n",
    "Let $ \\pi (z) = \\sum^m_{j=0} \\pi_j z^j $ and let $ z_1, \\ldots,\n",
    "z_k $ be the zeros of $ \\pi (z) $ that are inside the unit circle, $ k < m $.\n",
    "\n",
    "Then define\n",
    "\n",
    "$$\n",
    "\\theta (z) = \\pi (z) \\Biggl( {(z_1 z-1) \\over (z-z_1)} \\Biggr)\n",
    "\\Biggl( { (z_2 z-1) \\over (z-z_2) } \\Biggr ) \\ldots \\Biggl({(z_kz-1) \\over\n",
    "(z-z_k) }\\Biggr)\n",
    "$$\n",
    "\n",
    "The term multiplying $ \\pi (z) $ is termed a “Blaschke factor”.\n",
    "\n",
    "Then it can be proved directly that\n",
    "\n",
    "$$\n",
    "\\theta (z^{-1}) \\theta (z) = \\pi (z^{-1}) \\pi (z)\n",
    "$$\n",
    "\n",
    "and that the zeros of $ \\theta (z) $ are not inside the unit circle."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18ca08a9",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b0afa81",
   "metadata": {},
   "source": [
    "## Exercise 6.1\n",
    "\n",
    "Let $ Y_t = (1 - 2 L ) u_t $ where $ u_t $ is a mean zero\n",
    "white noise with $ \\mathbb{E} u^2_t = 1 $. Let\n",
    "\n",
    "$$\n",
    "X_t = Y_t + \\varepsilon_t\n",
    "$$\n",
    "\n",
    "where $ \\varepsilon_t $ is a serially uncorrelated white noise with\n",
    "$ \\mathbb{E} \\varepsilon^2_t = 9 $, and $ \\mathbb{E} \\varepsilon_t u_s = 0 $ for all\n",
    "$ t $ and $ s $.\n",
    "\n",
    "Find the Wold moving average representation for $ X_t $.\n",
    "\n",
    "Find a formula for the $ A_{1j} $’s in\n",
    "\n",
    "$$\n",
    "\\mathbb{E} \\widehat X_{t+1} \\mid X_t, X_{t-1}, \\ldots = \\sum^\\infty_{j=0} A_{1j}\n",
    "X_{t-j}\n",
    "$$\n",
    "\n",
    "Find a formula for the $ A_{2j} $’s in\n",
    "\n",
    "$$\n",
    "\\mathbb{\\hat E} X_{t+2} \\mid X_t, X_{t-1}, \\ldots = \\sum^\\infty_{j=0} A_{2j}\n",
    "X_{t-j}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a04d2cc",
   "metadata": {},
   "source": [
    "## Exercise 6.2\n",
    "\n",
    "**Multivariable Prediction:** Let $ Y_t $ be an $ (n\\times 1) $\n",
    "vector stochastic process with moving average representation\n",
    "\n",
    "$$\n",
    "Y_t = D(L) U_t\n",
    "$$\n",
    "\n",
    "where $ D(L) = \\sum^m_{j=0} D_j L^J, D_j $ an $ n \\times n $\n",
    "matrix, $ U_t $ an $ (n \\times 1) $ vector white noise with\n",
    "$ \\mathbb{E} U_t =0 $ for all $ t $, $ \\mathbb{E} U_t U_s' = 0 $ for all $ s \\neq t $,\n",
    "and $ \\mathbb{E} U_t U_t' = I $ for all $ t $.\n",
    "\n",
    "Let $ \\varepsilon_t $ be an $ n \\times 1 $ vector white noise with mean $ 0 $ and contemporaneous covariance matrix $ H $, where $ H $ is a positive definite matrix.\n",
    "\n",
    "Let $ X_t = Y_t +\\varepsilon_t $.\n",
    "\n",
    "Define the covariograms as $ C_X\n",
    "(\\tau) = \\mathbb{E} X_t X^\\prime_{t-\\tau}, C_Y (\\tau) = \\mathbb{E} Y_t Y^\\prime_{t-\\tau},\n",
    "C_{YX} (\\tau) = \\mathbb{E} Y_t X^\\prime_{t-\\tau} $.\n",
    "\n",
    "Then define the matrix\n",
    "covariance generating function, as in [(5.21)](https://tools-techniques.quantecon.org/lu_tricks.html#equation-onetwenty), only interpret all the\n",
    "objects in [(5.21)](https://tools-techniques.quantecon.org/lu_tricks.html#equation-onetwenty) as matrices.\n",
    "\n",
    "Show that the covariance generating functions are given by\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    g_y (z) &= D (z) D (z^{-1})^\\prime \\\\\n",
    "    g_X (z) &= D (z) D (z^{-1})^\\prime + H \\\\\n",
    "    g_{YX} (z) &= D (z) D (z^{-1})^\\prime\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "A factorization of $ g_X (z) $ can be found (see [[Rozanov, 1967](https://tools-techniques.quantecon.org/zreferences.html#id83)] or [[Whittle, 1983](https://tools-techniques.quantecon.org/zreferences.html#id85)]) of the form\n",
    "\n",
    "$$\n",
    "D (z) D (z^{-1})^\\prime + H = C (z) C (z^{-1})^\\prime, \\quad C (z) =\n",
    "\\sum^m_{j=0} C_j z^j\n",
    "$$\n",
    "\n",
    "where the zeros of $ \\vert C(z)\\vert $ do not lie inside the unit\n",
    "circle.\n",
    "\n",
    "A vector Wold moving average representation of $ X_t $ is then\n",
    "\n",
    "$$\n",
    "X_t = C(L) \\eta_t\n",
    "$$\n",
    "\n",
    "where $ \\eta_t $ is an $ (n\\times 1) $ vector white noise that\n",
    "is “fundamental” for $ X_t $.\n",
    "\n",
    "That is, $ X_t - \\mathbb{\\hat E}\\left[X_t \\mid X_{t-1}, X_{t-2}\n",
    "\\ldots\\right] = C_0 \\, \\eta_t $.\n",
    "\n",
    "The optimum predictor of $ X_{t+j} $ is\n",
    "\n",
    "$$\n",
    "\\mathbb{\\hat E} \\left[X_{t+j} \\mid X_t, X_{t-1}, \\ldots\\right]\n",
    " = \\left[{C(L) \\over L^j} \\right]_+ \\eta_t\n",
    "$$\n",
    "\n",
    "If $ C(L) $ is invertible, i.e., if the zeros of $ \\det $\n",
    "$ C(z) $ lie strictly outside the unit circle, then this formula can\n",
    "be written\n",
    "\n",
    "$$\n",
    "\\mathbb{\\hat E} \\left[X_{t+j} \\mid X_t, X_{t-1}, \\ldots\\right]\n",
    "    = \\left[{C(L) \\over L^J} \\right]_+ C(L)^{-1}\\, X_t\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "date": 1718249790.4545488,
  "filename": "classical_filtering.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Classical Prediction and Filtering With Linear Algebra"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}