arma.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='arma'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Covariance Stationary Processes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Covariance Stationary Processes](#Covariance-Stationary-Processes)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Introduction](#Introduction)  \n",
    "  - [Spectral Analysis](#Spectral-Analysis)  \n",
    "  - [Implementation](#Implementation)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture we study covariance stationary linear stochastic processes, a\n",
    "class of models routinely used to study economic and financial time series.\n",
    "\n",
    "This class has the advantage of being\n",
    "\n",
    "1. simple enough to be described by an elegant and comprehensive theory  \n",
    "1. relatively broad in terms of the kinds of dynamics it can represent  \n",
    "\n",
    "\n",
    "We consider these models in both the time and frequency domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARMA Processes\n",
    "\n",
    "We will focus much of our attention on linear covariance stationary models with a finite number of parameters.\n",
    "\n",
    "In particular, we will study stationary ARMA processes, which form a cornerstone of the standard theory of time series analysis.\n",
    "\n",
    "Every ARMA process can be represented in [linear state space](https://lectures.quantecon.org/linear_models.html) form.\n",
    "\n",
    "However, ARMA processes have some important structure that makes it valuable to study them separately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral Analysis\n",
    "\n",
    "Analysis in the frequency domain is also called spectral analysis.\n",
    "\n",
    "In essence, spectral analysis provides an alternative representation of the\n",
    "autocovariance function of a covariance stationary process.\n",
    "\n",
    "Having a second representation of this important object\n",
    "\n",
    "- shines a light on the dynamics of the process in question  \n",
    "- allows for a simpler, more tractable representation in some important cases  \n",
    "\n",
    "\n",
    "The famous *Fourier transform* and its inverse are used to map between the two representations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Reading\n",
    "\n",
    "For supplementary reading, see\n",
    "\n",
    "- [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012), chapter 2  \n",
    "- [[Sar87]](https://lectures.quantecon.org/zreferences.html#sargent1987), chapter 11  \n",
    "- John Cochrane’s <a href=/_static/lecture_specific/arma/time_series_book.pdf download>notes on time series analysis</a>, chapter 8  \n",
    "- [[Shi95]](https://lectures.quantecon.org/zreferences.html#shiryaev1995), chapter 6  \n",
    "- [[CC08]](https://lectures.quantecon.org/zreferences.html#cryerchan2008), all  \n",
    "\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import quantecon as qe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Consider a sequence of random variables $ \\{ X_t \\} $ indexed by $ t \\in \\mathbb Z $ and taking values in $ \\mathbb R $.\n",
    "\n",
    "Thus, $ \\{ X_t \\} $ begins in the infinite past and extends to the infinite future — a convenient and standard assumption.\n",
    "\n",
    "As in other fields, successful economic modeling typically assumes the existence of features that are constant over time.\n",
    "\n",
    "If these assumptions are correct, then each new observation $ X_t, X_{t+1},\\ldots $ can provide additional information about the time-invariant features, allowing us to  learn from as data arrive.\n",
    "\n",
    "For this reason, we will focus in what follows on processes that are *stationary* — or become so after a transformation\n",
    "(see for example [this lecture](https://lectures.quantecon.org/additive_functionals.html)).\n",
    "\n",
    "\n",
    "<a id='arma-defs'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definitions\n",
    "\n",
    "\n",
    "<a id='index-3'></a>\n",
    "A real-valued stochastic process $ \\{ X_t \\} $ is called *covariance stationary* if\n",
    "\n",
    "1. Its mean $ \\mu := \\mathbb E X_t $ does not depend on $ t $.  \n",
    "1. For all $ k $ in $ \\mathbb Z $, the $ k $-th autocovariance $ \\gamma(k) := \\mathbb E (X_t - \\mu)(X_{t + k} - \\mu) $ is finite and depends only on $ k $.  \n",
    "\n",
    "\n",
    "The function $ \\gamma \\colon \\mathbb Z \\to \\mathbb R $ is called the *autocovariance function* of the process.\n",
    "\n",
    "Throughout this lecture, we will work exclusively with zero-mean (i.e., $ \\mu = 0 $) covariance stationary processes.\n",
    "\n",
    "The zero-mean assumption costs nothing in terms of generality since working with non-zero-mean processes involves no more than adding a constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: White Noise\n",
    "\n",
    "Perhaps the simplest class of covariance stationary processes is the white noise processes.\n",
    "\n",
    "A process $ \\{ \\epsilon_t \\} $ is called a *white noise process* if\n",
    "\n",
    "1. $ \\mathbb E \\epsilon_t = 0 $  \n",
    "1. $ \\gamma(k) = \\sigma^2 \\mathbf 1\\{k = 0\\} $ for some $ \\sigma > 0 $  \n",
    "\n",
    "\n",
    "(Here $ \\mathbf 1\\{k = 0\\} $ is defined to be 1 if $ k = 0 $ and zero otherwise)\n",
    "\n",
    "White noise processes play the role of **building blocks** for processes with more complicated dynamics.\n",
    "\n",
    "\n",
    "<a id='generalized-lps'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: General Linear Processes\n",
    "\n",
    "From the simple building block provided by white noise, we can construct a very flexible family of covariance stationary processes — the *general linear processes*\n",
    "\n",
    "\n",
    "<a id='equation-ma-inf'></a>\n",
    "$$\n",
    "X_t = \\sum_{j=0}^{\\infty} \\psi_j \\epsilon_{t-j},\n",
    "\\qquad t \\in \\mathbb Z \\tag{1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ \\{\\epsilon_t\\} $ is white noise  \n",
    "- $ \\{\\psi_t\\} $ is a square summable sequence in $ \\mathbb R $ (that is, $ \\sum_{t=0}^{\\infty} \\psi_t^2 < \\infty $)  \n",
    "\n",
    "\n",
    "The sequence $ \\{\\psi_t\\} $ is often called a *linear filter*.\n",
    "\n",
    "Equation [(1)](#equation-ma-inf) is said to present  a **moving average** process or a moving average representation.\n",
    "\n",
    "With some manipulations, it is possible to confirm that the autocovariance function for [(1)](#equation-ma-inf) is\n",
    "\n",
    "\n",
    "<a id='equation-ma-inf-ac'></a>\n",
    "$$\n",
    "\\gamma(k) = \\sigma^2 \\sum_{j=0}^{\\infty} \\psi_j \\psi_{j+k} \\tag{2}\n",
    "$$\n",
    "\n",
    "By the [Cauchy-Schwartz inequality](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Schwarz_inequality), one can show that $ \\gamma(k) $ satisfies equation [(2)](#equation-ma-inf-ac).\n",
    "\n",
    "Evidently, $ \\gamma(k) $ does not depend on $ t $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wold Representation\n",
    "\n",
    "Remarkably, the class of general linear processes goes a long way towards\n",
    "describing the entire class of zero-mean covariance stationary processes.\n",
    "\n",
    "In particular, [Wold’s decomposition theorem](https://en.wikipedia.org/wiki/Wold%27s_theorem) states that every\n",
    "zero-mean covariance stationary process $ \\{X_t\\} $ can be written as\n",
    "\n",
    "$$\n",
    "X_t = \\sum_{j=0}^{\\infty} \\psi_j \\epsilon_{t-j} + \\eta_t\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ \\{\\epsilon_t\\} $ is white noise  \n",
    "- $ \\{\\psi_t\\} $ is square summable  \n",
    "- $ \\psi_0 \\epsilon_t $ is the one-step ahead prediction error in forecasting $ X_t $ as a linear least-squares function of the infinite history $ X_{t-1}, X_{t-2}, \\ldots $  \n",
    "- $ \\eta_t $ can be expressed as a linear function of $ X_{t-1}, X_{t-2},\\ldots $ and is perfectly predictable over arbitrarily long horizons  \n",
    "\n",
    "\n",
    "For the method of constructing a Wold representation, intuition, and further discussion, see [[Sar87]](https://lectures.quantecon.org/zreferences.html#sargent1987), p. 286."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AR and MA\n",
    "\n",
    "\n",
    "<a id='index-8'></a>\n",
    "General linear processes are a very broad class of processes.\n",
    "\n",
    "It often pays to specialize to those for which there exists a representation having only finitely many parameters.\n",
    "\n",
    "(Experience and theory combine to indicate that models with a relatively small number of parameters typically perform better than larger models, especially for forecasting)\n",
    "\n",
    "One very simple example of such a model is the first-order autoregressive or AR(1) process\n",
    "\n",
    "\n",
    "<a id='equation-ar1-rep'></a>\n",
    "$$\n",
    "X_t = \\phi X_{t-1} + \\epsilon_t\n",
    "\\quad \\text{where} \\quad\n",
    "| \\phi | < 1\n",
    "\\quad \\text{and } \\{ \\epsilon_t \\} \\text{ is white noise} \\tag{3}\n",
    "$$\n",
    "\n",
    "By direct substitution, it is easy to verify that $ X_t = \\sum_{j=0}^{\\infty} \\phi^j \\epsilon_{t-j} $.\n",
    "\n",
    "Hence $ \\{X_t\\} $ is a general linear process.\n",
    "\n",
    "Applying [(2)](#equation-ma-inf-ac) to the previous expression for $ X_t $, we get the AR(1) autocovariance function\n",
    "\n",
    "\n",
    "<a id='equation-ar-acov'></a>\n",
    "$$\n",
    "\\gamma(k) = \\phi^k \\frac{\\sigma^2}{1 - \\phi^2},\n",
    "\\qquad k = 0, 1, \\ldots \\tag{4}\n",
    "$$\n",
    "\n",
    "The next figure plots an example of this function for $ \\phi = 0.8 $ and $ \\phi = -0.8 $ with $ \\sigma = 1 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "num_rows, num_cols = 2, 1\n",
    "fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 8))\n",
    "plt.subplots_adjust(hspace=0.4)\n",
    "\n",
    "for i, ϕ in enumerate((0.8, -0.8)):\n",
    "    ax = axes[i]\n",
    "    times = list(range(16))\n",
    "    acov = [ϕ**k / (1 - ϕ**2) for k in times]\n",
    "    ax.plot(times, acov, 'bo-', alpha=0.6,\n",
    "            label=f'autocovariance, $\\phi = {ϕ:.2}$')\n",
    "    ax.legend(loc='upper right')\n",
    "    ax.set(xlabel='time', xlim=(0, 15))\n",
    "    ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another very simple process is the MA(1) process (here MA means “moving average”)\n",
    "\n",
    "$$\n",
    "X_t = \\epsilon_t + \\theta \\epsilon_{t-1}\n",
    "$$\n",
    "\n",
    "You will be able to verify that\n",
    "\n",
    "$$\n",
    "\\gamma(0) = \\sigma^2 (1 + \\theta^2),\n",
    "\\quad\n",
    "\\gamma(1) = \\sigma^2 \\theta,\n",
    "\\quad \\text{and} \\quad\n",
    "\\gamma(k) = 0 \\quad \\forall \\, k > 1\n",
    "$$\n",
    "\n",
    "The AR(1) can be generalized to an AR($ p $) and likewise for the MA(1).\n",
    "\n",
    "Putting all of this together, we get the"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARMA Processes\n",
    "\n",
    "A stochastic process $ \\{X_t\\} $ is called an *autoregressive moving\n",
    "average process*, or ARMA($ p,q $), if it can be written as\n",
    "\n",
    "\n",
    "<a id='equation-arma'></a>\n",
    "$$\n",
    "X_t = \\phi_1 X_{t-1} + \\cdots + \\phi_p X_{t-p} +\n",
    "    \\epsilon_t + \\theta_1 \\epsilon_{t-1} + \\cdots + \\theta_q \\epsilon_{t-q} \\tag{5}\n",
    "$$\n",
    "\n",
    "where $ \\{ \\epsilon_t \\} $ is white noise.\n",
    "\n",
    "An alternative notation for ARMA processes uses the *lag operator* $ L $.\n",
    "\n",
    "**Def.** Given arbitrary variable $ Y_t $, let $ L^k Y_t := Y_{t-k} $.\n",
    "\n",
    "It turns out that\n",
    "\n",
    "- lag operators facilitate  succinct representations for linear stochastic processes  \n",
    "- algebraic manipulations that treat the lag operator as an ordinary scalar  are legitimate  \n",
    "\n",
    "\n",
    "Using $ L $, we can rewrite [(5)](#equation-arma) as\n",
    "\n",
    "\n",
    "<a id='equation-arma-lag'></a>\n",
    "$$\n",
    "L^0 X_t - \\phi_1 L^1 X_t - \\cdots - \\phi_p L^p X_t\n",
    "= L^0 \\epsilon_t + \\theta_1 L^1 \\epsilon_t + \\cdots + \\theta_q L^q \\epsilon_t \\tag{6}\n",
    "$$\n",
    "\n",
    "If we let $ \\phi(z) $ and $ \\theta(z) $ be the polynomials\n",
    "\n",
    "\n",
    "<a id='equation-arma-poly'></a>\n",
    "$$\n",
    "\\phi(z) := 1 - \\phi_1 z - \\cdots - \\phi_p z^p\n",
    "\\quad \\text{and} \\quad\n",
    "\\theta(z) := 1 + \\theta_1 z + \\cdots + \\theta_q z^q \\tag{7}\n",
    "$$\n",
    "\n",
    "then [(6)](#equation-arma-lag)  becomes\n",
    "\n",
    "\n",
    "<a id='equation-arma-lag1'></a>\n",
    "$$\n",
    "\\phi(L) X_t = \\theta(L) \\epsilon_t \\tag{8}\n",
    "$$\n",
    "\n",
    "In what follows we **always assume** that the roots of the polynomial $ \\phi(z) $ lie outside the unit circle in the complex plane.\n",
    "\n",
    "This condition is sufficient to guarantee that the ARMA($ p,q $) process is covariance stationary.\n",
    "\n",
    "In fact, it implies that the process falls within the class of general linear processes [described above](#generalized-lps).\n",
    "\n",
    "That is, given an ARMA($ p,q $) process $ \\{ X_t \\} $ satisfying the unit circle condition, there exists a square summable sequence $ \\{\\psi_t\\} $ with $ X_t = \\sum_{j=0}^{\\infty} \\psi_j \\epsilon_{t-j} $ for all $ t $.\n",
    "\n",
    "The sequence $ \\{\\psi_t\\} $ can be obtained by a recursive procedure outlined on page 79 of [[CC08]](https://lectures.quantecon.org/zreferences.html#cryerchan2008).\n",
    "\n",
    "The function $ t \\mapsto \\psi_t $ is often called the *impulse response function*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spectral Analysis\n",
    "\n",
    "Autocovariance functions provide a great deal of information about covariance stationary processes.\n",
    "\n",
    "In fact, for zero-mean Gaussian processes, the autocovariance function characterizes the entire joint distribution.\n",
    "\n",
    "Even for non-Gaussian processes, it provides a significant amount of information.\n",
    "\n",
    "It turns out that there is an alternative representation of the autocovariance function of a covariance stationary process, called the *spectral density*.\n",
    "\n",
    "At times, the spectral density is easier to derive, easier to manipulate, and provides additional intuition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complex Numbers\n",
    "\n",
    "Before discussing the spectral density, we invite you to recall the main properties of complex numbers (or [skip to the next section](#arma-specd)).\n",
    "\n",
    "It can be helpful to remember that, in a formal sense, complex numbers are just points $ (x, y) \\in \\mathbb R^2 $ endowed with a specific notion of multiplication.\n",
    "\n",
    "When $ (x, y) $ is regarded as a complex number, $ x $ is called the *real part* and $ y $ is called the *imaginary part*.\n",
    "\n",
    "The *modulus* or *absolute value* of a complex number $ z = (x, y) $ is just its Euclidean norm in $ \\mathbb R^2 $, but is usually written as $ |z| $ instead of $ \\|z\\| $.\n",
    "\n",
    "The product of two complex numbers $ (x, y) $ and $ (u, v) $ is defined to be $ (xu - vy, xv + yu) $, while addition is standard pointwise vector addition.\n",
    "\n",
    "When endowed with these notions of multiplication and addition, the set of complex numbers forms a [field](https://en.wikipedia.org/wiki/Field_%28mathematics%29) — addition and multiplication play well together, just as they do in $ \\mathbb R $.\n",
    "\n",
    "The complex number $ (x, y) $ is often written as $ x + i y $, where $ i $ is called the *imaginary unit* and is understood to obey $ i^2 = -1 $.\n",
    "\n",
    "The $ x + i y $ notation provides an easy way to remember the definition of multiplication given above, because, proceeding naively,\n",
    "\n",
    "$$\n",
    "(x + i y) (u + i v) = xu - yv + i (xv + yu)\n",
    "$$\n",
    "\n",
    "Converted back to our first notation, this becomes $ (xu - vy, xv + yu) $ as promised.\n",
    "\n",
    "Complex numbers can be represented in  the polar form $ r e^{i \\omega} $ where\n",
    "\n",
    "$$\n",
    "r e^{i \\omega} := r (\\cos(\\omega) + i \\sin(\\omega)) = x + i y\n",
    "$$\n",
    "\n",
    "where $ x = r \\cos(\\omega), y = r \\sin(\\omega) $, and $ \\omega = \\arctan(y/z) $ or $ \\tan(\\omega) = y/x $.\n",
    "\n",
    "\n",
    "<a id='arma-specd'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spectral Densities\n",
    "\n",
    "Let $ \\{ X_t \\} $ be a covariance stationary process with autocovariance function $ \\gamma $  satisfying $ \\sum_{k} \\gamma(k)^2 < \\infty $.\n",
    "\n",
    "The *spectral density* $ f $ of $ \\{ X_t \\} $ is defined as the [discrete time Fourier transform](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform) of its autocovariance function $ \\gamma $.\n",
    "\n",
    "$$\n",
    "f(\\omega) := \\sum_{k \\in \\mathbb Z} \\gamma(k) e^{-i \\omega k},\n",
    "\\qquad \\omega \\in \\mathbb R\n",
    "$$\n",
    "\n",
    "(Some authors normalize the expression on the right by constants such as $ 1/\\pi $ — the convention chosen  makes little difference provided you are consistent).\n",
    "\n",
    "Using the fact that $ \\gamma $ is *even*, in the sense that $ \\gamma(t) = \\gamma(-t) $ for all $ t $, we can show that\n",
    "\n",
    "\n",
    "<a id='equation-arma-sd-cos'></a>\n",
    "$$\n",
    "f(\\omega) = \\gamma(0) + 2 \\sum_{k \\geq 1} \\gamma(k) \\cos(\\omega k) \\tag{9}\n",
    "$$\n",
    "\n",
    "It is not difficult to confirm that $ f $ is\n",
    "\n",
    "- real-valued  \n",
    "- even ($ f(\\omega) = f(-\\omega) $   ),  and  \n",
    "- $ 2\\pi $-periodic, in the sense that $ f(2\\pi + \\omega) = f(\\omega) $ for all $ \\omega $  \n",
    "\n",
    "\n",
    "It follows that the values of $ f $ on $ [0, \\pi] $ determine the values of $ f $ on\n",
    "all of $ \\mathbb R $ — the proof is an exercise.\n",
    "\n",
    "For this reason, it is standard to plot the spectral density only on the interval $ [0, \\pi] $.\n",
    "\n",
    "\n",
    "<a id='arma-wnsd'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 1: White Noise\n",
    "\n",
    "Consider a white noise process $ \\{\\epsilon_t\\} $ with standard deviation $ \\sigma $.\n",
    "\n",
    "It is easy to check that in  this case $ f(\\omega) = \\sigma^2 $.  So $ f $ is a constant function.\n",
    "\n",
    "As we will see, this can be interpreted as meaning that “all frequencies are equally present”.\n",
    "\n",
    "(White light has this property when frequency refers to the visible spectrum, a connection that provides the origins of the term “white noise”)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: AR and MA and ARMA\n",
    "\n",
    "It is an exercise to show that the MA(1) process $ X_t = \\theta \\epsilon_{t-1} + \\epsilon_t $ has a spectral density\n",
    "\n",
    "\n",
    "<a id='equation-ma1-sd-ed'></a>\n",
    "$$\n",
    "f(\\omega)\n",
    "= \\sigma^2 ( 1 + 2 \\theta \\cos(\\omega) + \\theta^2 ) \\tag{10}\n",
    "$$\n",
    "\n",
    "With a bit more effort, it’s possible to show (see, e.g., p. 261 of [[Sar87]](https://lectures.quantecon.org/zreferences.html#sargent1987)) that the spectral density of the AR(1) process $ X_t = \\phi X_{t-1} + \\epsilon_t $ is\n",
    "\n",
    "\n",
    "<a id='equation-ar1-sd-ed'></a>\n",
    "$$\n",
    "f(\\omega)\n",
    "= \\frac{\\sigma^2}{ 1 - 2 \\phi \\cos(\\omega) + \\phi^2 } \\tag{11}\n",
    "$$\n",
    "\n",
    "More generally, it can be shown that the spectral density of the ARMA process [(5)](#equation-arma) is\n",
    "\n",
    "\n",
    "<a id='arma-spec-den'></a>\n",
    "\n",
    "<a id='equation-arma-sd'></a>\n",
    "$$\n",
    "f(\\omega) = \\left| \\frac{\\theta(e^{i\\omega})}{\\phi(e^{i\\omega})} \\right|^2 \\sigma^2 \\tag{12}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ \\sigma $ is the standard deviation of the white noise process $ \\{\\epsilon_t\\} $.  \n",
    "- the polynomials $ \\phi(\\cdot) $ and $ \\theta(\\cdot) $ are as defined in [(7)](#equation-arma-poly).  \n",
    "\n",
    "\n",
    "The derivation of [(12)](#equation-arma-sd) uses the fact that convolutions become products under Fourier transformations.\n",
    "\n",
    "The proof is elegant and can be found in many places — see, for example, [[Sar87]](https://lectures.quantecon.org/zreferences.html#sargent1987), chapter 11, section 4.\n",
    "\n",
    "It’s a nice exercise to verify that [(10)](#equation-ma1-sd-ed) and [(11)](#equation-ar1-sd-ed) are indeed special cases of [(12)](#equation-arma-sd)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting the Spectral Density\n",
    "\n",
    "\n",
    "<a id='index-18'></a>\n",
    "Plotting [(11)](#equation-ar1-sd-ed) reveals the shape of the spectral density for the AR(1) model when $ \\phi $ takes the values 0.8 and -0.8 respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def ar1_sd(ϕ, ω):\n",
    "    return 1 / (1 - 2 * ϕ * np.cos(ω) + ϕ**2)\n",
    "\n",
    "ωs = np.linspace(0, np.pi, 180)\n",
    "num_rows, num_cols = 2, 1\n",
    "fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 8))\n",
    "plt.subplots_adjust(hspace=0.4)\n",
    "\n",
    "# Autocovariance when phi = 0.8\n",
    "for i, ϕ in enumerate((0.8, -0.8)):\n",
    "    ax = axes[i]\n",
    "    sd = ar1_sd(ϕ, ωs)\n",
    "    ax.plot(ωs, sd, 'b-', alpha=0.6, lw=2,\n",
    "            label='spectral density, $\\phi = {ϕ:.2}$')\n",
    "    ax.legend(loc='upper center')\n",
    "    ax.set(xlabel='frequency', xlim=(0, np.pi))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These spectral densities correspond to the autocovariance functions for the\n",
    "AR(1) process shown above.\n",
    "\n",
    "Informally, we think of the spectral density as being large at those $ \\omega \\in [0, \\pi] $ at which\n",
    "the autocovariance function seems approximately to exhibit big damped cycles.\n",
    "\n",
    "To see the idea, let’s consider why, in the lower panel of the preceding figure, the spectral density for the case $ \\phi = -0.8 $ is large at $ \\omega = \\pi $.\n",
    "\n",
    "Recall that the spectral density can be expressed as\n",
    "\n",
    "\n",
    "<a id='equation-sumpr'></a>\n",
    "$$\n",
    "f(\\omega)\n",
    "= \\gamma(0) + 2 \\sum_{k \\geq 1} \\gamma(k) \\cos(\\omega k)\n",
    "= \\gamma(0) + 2 \\sum_{k \\geq 1} (-0.8)^k \\cos(\\omega k) \\tag{13}\n",
    "$$\n",
    "\n",
    "When we evaluate this at $ \\omega = \\pi $, we get a large number because\n",
    "$ \\cos(\\pi k) $ is large and positive when $ (-0.8)^k $ is\n",
    "positive, and large in absolute value and negative when $ (-0.8)^k $ is negative.\n",
    "\n",
    "Hence the product is always large and positive, and hence the sum of the\n",
    "products on the right-hand side of [(13)](#equation-sumpr) is large.\n",
    "\n",
    "These ideas are illustrated in the next figure, which has $ k $ on the horizontal axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = -0.8\n",
    "times = list(range(16))\n",
    "y1 = [ϕ**k / (1 - ϕ**2) for k in times]\n",
    "y2 = [np.cos(np.pi * k) for k in times]\n",
    "y3 = [a * b for a, b in zip(y1, y2)]\n",
    "\n",
    "num_rows, num_cols = 3, 1\n",
    "fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 8))\n",
    "plt.subplots_adjust(hspace=0.25)\n",
    "\n",
    "# Autocovariance when ϕ = -0.8\n",
    "ax = axes[0]\n",
    "ax.plot(times, y1, 'bo-', alpha=0.6, label='$\\gamma(k)$')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set(xlim=(0, 15), yticks=(-2, 0, 2))\n",
    "ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "\n",
    "# Cycles at frequency π\n",
    "ax = axes[1]\n",
    "ax.plot(times, y2, 'bo-', alpha=0.6, label='$\\cos(\\pi k)$')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set(xlim=(0, 15), yticks=(-1, 0, 1))\n",
    "ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "\n",
    "# Product\n",
    "ax = axes[2]\n",
    "ax.stem(times, y3, label='$\\gamma(k) \\cos(\\pi k)$')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set(xlim=(0, 15), ylim=(-3, 3), yticks=(-1, 0, 1, 2, 3))\n",
    "ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "ax.set_xlabel(\"k\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, if we evaluate $ f(\\omega) $ at $ \\omega = \\pi / 3 $, then the cycles are\n",
    "not matched, the sequence $ \\gamma(k) \\cos(\\omega k) $ contains\n",
    "both positive and negative terms, and hence the sum of these terms is much smaller."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = -0.8\n",
    "times = list(range(16))\n",
    "y1 = [ϕ**k / (1 - ϕ**2) for k in times]\n",
    "y2 = [np.cos(np.pi * k/3) for k in times]\n",
    "y3 = [a * b for a, b in zip(y1, y2)]\n",
    "\n",
    "num_rows, num_cols = 3, 1\n",
    "fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 8))\n",
    "plt.subplots_adjust(hspace=0.25)\n",
    "\n",
    "# Autocovariance when phi = -0.8\n",
    "ax = axes[0]\n",
    "ax.plot(times, y1, 'bo-', alpha=0.6, label='$\\gamma(k)$')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set(xlim=(0, 15), yticks=(-2, 0, 2))\n",
    "ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "\n",
    "# Cycles at frequency π\n",
    "ax = axes[1]\n",
    "ax.plot(times, y2, 'bo-', alpha=0.6, label='$\\cos(\\pi k/3)$')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set(xlim=(0, 15), yticks=(-1, 0, 1))\n",
    "ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "\n",
    "# Product\n",
    "ax = axes[2]\n",
    "ax.stem(times, y3, label='$\\gamma(k) \\cos(\\pi k/3)$')\n",
    "ax.legend(loc='upper right')\n",
    "ax.set(xlim=(0, 15), ylim=(-3, 3), yticks=(-1, 0, 1, 2, 3))\n",
    "ax.hlines(0, 0, 15, linestyle='--', alpha=0.5)\n",
    "ax.set_xlabel(\"$k$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, the spectral density is large at frequencies $ \\omega $ where the autocovariance function exhibits damped cycles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inverting the Transformation\n",
    "\n",
    "\n",
    "<a id='index-19'></a>\n",
    "We have just seen that the spectral density is useful in the sense that it provides a frequency-based perspective on the autocovariance structure of a covariance stationary process.\n",
    "\n",
    "Another reason that the spectral density is useful is that it can be “inverted” to recover the autocovariance function via the *inverse Fourier transform*.\n",
    "\n",
    "In particular, for all $ k \\in \\mathbb Z $, we have\n",
    "\n",
    "\n",
    "<a id='equation-ift'></a>\n",
    "$$\n",
    "\\gamma(k) = \\frac{1}{2 \\pi} \\int_{-\\pi}^{\\pi} f(\\omega) e^{i \\omega k} d\\omega \\tag{14}\n",
    "$$\n",
    "\n",
    "This is convenient in situations where the spectral density is easier to calculate and manipulate than the autocovariance function.\n",
    "\n",
    "(For example, the expression [(12)](#equation-arma-sd) for the ARMA spectral density is much easier to work with than the expression for the ARMA autocovariance)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mathematical Theory\n",
    "\n",
    "\n",
    "<a id='index-20'></a>\n",
    "This section is loosely based on [[Sar87]](https://lectures.quantecon.org/zreferences.html#sargent1987), p. 249-253, and included for those who\n",
    "\n",
    "- would like a bit more insight into spectral densities  \n",
    "- and have at least some background in [Hilbert space](https://en.wikipedia.org/wiki/Hilbert_space) theory  \n",
    "\n",
    "\n",
    "Others should feel free to skip to the [next section](#arma-imp) — none of this material is necessary to progress to computation.\n",
    "\n",
    "Recall that every [separable](https://en.wikipedia.org/wiki/Separable_space) Hilbert space $ H $ has a countable orthonormal basis $ \\{ h_k \\} $.\n",
    "\n",
    "The nice thing about such a basis is that every $ f \\in H $ satisfies\n",
    "\n",
    "\n",
    "<a id='equation-arma-fc'></a>\n",
    "$$\n",
    "f = \\sum_k \\alpha_k h_k\n",
    "\\quad \\text{where} \\quad\n",
    "\\alpha_k := \\langle f, h_k \\rangle \\tag{15}\n",
    "$$\n",
    "\n",
    "where $ \\langle \\cdot, \\cdot \\rangle $ denotes the inner product in $ H $.\n",
    "\n",
    "Thus, $ f $ can be represented to any degree of precision by linearly combining basis vectors.\n",
    "\n",
    "The scalar sequence $ \\alpha = \\{\\alpha_k\\} $ is called the *Fourier coefficients* of $ f $, and satisfies $ \\sum_k |\\alpha_k|^2 < \\infty $.\n",
    "\n",
    "In other words, $ \\alpha $ is in $ \\ell_2 $, the set of square summable sequences.\n",
    "\n",
    "Consider an operator $ T $ that maps $ \\alpha \\in \\ell_2 $ into its expansion $ \\sum_k \\alpha_k h_k \\in H $.\n",
    "\n",
    "The Fourier coefficients of $ T\\alpha $ are just $ \\alpha = \\{ \\alpha_k \\} $, as you can verify by confirming that $ \\langle T \\alpha, h_k \\rangle = \\alpha_k $.\n",
    "\n",
    "Using elementary results from Hilbert space theory, it can be shown that\n",
    "\n",
    "- $ T $ is one-to-one — if $ \\alpha $ and $ \\beta $ are distinct in $ \\ell_2 $, then so are their expansions in $ H $.  \n",
    "- $ T $ is onto — if $ f \\in H $ then its preimage in $ \\ell_2 $ is the sequence $ \\alpha $ given by $ \\alpha_k = \\langle f, h_k \\rangle $.  \n",
    "- $ T $ is a linear isometry — in particular, $ \\langle \\alpha, \\beta \\rangle = \\langle T\\alpha, T\\beta \\rangle $.  \n",
    "\n",
    "\n",
    "Summarizing these results, we say that any separable Hilbert space is isometrically isomorphic to $ \\ell_2 $.\n",
    "\n",
    "In essence, this says that each separable Hilbert space we consider is just a different way of looking at the fundamental space $ \\ell_2 $.\n",
    "\n",
    "With this in mind, let’s specialize to a setting where\n",
    "\n",
    "- $ \\gamma \\in \\ell_2 $ is the autocovariance function of a covariance stationary process, and $ f $ is the spectral density.  \n",
    "- $ H = L_2 $, where $ L_2 $ is the set of square summable functions on the interval $ [-\\pi, \\pi] $, with inner product $ \\langle g, h \\rangle = \\int_{-\\pi}^{\\pi} g(\\omega) h(\\omega) d \\omega $.  \n",
    "- $ \\{h_k\\} = $ the orthonormal basis for $ L_2 $ given by the set of trigonometric functions.  \n",
    "\n",
    "\n",
    "$$\n",
    "h_k(\\omega) = \\frac{e^{i \\omega k}}{\\sqrt{2 \\pi}},\n",
    "\\quad k \\in \\mathbb Z,\n",
    "\\quad \\omega \\in [-\\pi, \\pi]\n",
    "$$\n",
    "\n",
    "Using the definition of $ T $ from above and the fact that $ f $ is even, we now have\n",
    "\n",
    "\n",
    "<a id='equation-arma-it'></a>\n",
    "$$\n",
    "T \\gamma\n",
    "= \\sum_{k \\in \\mathbb Z}\n",
    "\\gamma(k) \\frac{e^{i \\omega k}}{\\sqrt{2 \\pi}} = \\frac{1}{\\sqrt{2 \\pi}} f(\\omega) \\tag{16}\n",
    "$$\n",
    "\n",
    "In other words, apart from a scalar multiple, the spectral density is just a transformation of $ \\gamma \\in \\ell_2 $ under a certain linear isometry — a different way to view $ \\gamma $.\n",
    "\n",
    "In particular, it is an expansion of the autocovariance function with respect to the trigonometric basis functions in $ L_2 $.\n",
    "\n",
    "As discussed above, the Fourier coefficients of $ T \\gamma $ are given by the sequence $ \\gamma $, and,\n",
    "in particular, $ \\gamma(k) = \\langle T \\gamma, h_k \\rangle $.\n",
    "\n",
    "Transforming this inner product into its integral expression and using [(16)](#equation-arma-it) gives\n",
    "[(14)](#equation-ift), justifying our earlier expression for the inverse transform.\n",
    "\n",
    "\n",
    "<a id='arma-imp'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Most code for working with covariance stationary models deals with ARMA models.\n",
    "\n",
    "Python code for studying ARMA models can be found in the `tsa` submodule of [statsmodels](http://statsmodels.sourceforge.net/).\n",
    "\n",
    "Since this code doesn’t quite cover our needs — particularly vis-a-vis spectral analysis — we’ve put together the module [arma.py](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/arma.py), which is part of [QuantEcon.py](http://quantecon.org/python_index.html) package.\n",
    "\n",
    "The module provides functions for mapping ARMA($ p,q $) models into their\n",
    "\n",
    "1. impulse response function  \n",
    "1. simulated time series  \n",
    "1. autocovariance function  \n",
    "1. spectral density  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Application\n",
    "\n",
    "Let’s use this code to replicate the plots on pages 68–69 of [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012).\n",
    "\n",
    "Here are some functions to generate the plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_impulse_response(arma, ax=None):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    yi = arma.impulse_response()\n",
    "    ax.stem(list(range(len(yi))), yi)\n",
    "    ax.set(xlim=(-0.5), ylim=(min(yi)-0.1, max(yi)+0.1),\n",
    "                 title='Impulse response', xlabel='time', ylabel='response')\n",
    "    return ax\n",
    "\n",
    "def plot_spectral_density(arma, ax=None):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    w, spect = arma.spectral_density(two_pi=False)\n",
    "    ax.semilogy(w, spect)\n",
    "    ax.set(xlim=(0, np.pi), ylim=(0, np.max(spect)),\n",
    "           title='Spectral density', xlabel='frequency', ylabel='spectrum')\n",
    "    return ax\n",
    "\n",
    "def plot_autocovariance(arma, ax=None):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    acov = arma.autocovariance()\n",
    "    ax.stem(list(range(len(acov))), acov)\n",
    "    ax.set(xlim=(-0.5, len(acov) - 0.5), title='Autocovariance',\n",
    "           xlabel='time', ylabel='autocovariance')\n",
    "    return ax\n",
    "\n",
    "def plot_simulation(arma, ax=None):\n",
    "    if ax is None:\n",
    "        ax = plt.gca()\n",
    "    x_out = arma.simulation()\n",
    "    ax.plot(x_out)\n",
    "    ax.set(title='Sample path', xlabel='time', ylabel='state space')\n",
    "    return ax\n",
    "\n",
    "def quad_plot(arma):\n",
    "    \"\"\"\n",
    "    Plots the impulse response, spectral_density, autocovariance,\n",
    "    and one realization of the process.\n",
    "\n",
    "    \"\"\"\n",
    "    num_rows, num_cols = 2, 2\n",
    "    fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 8))\n",
    "    plot_functions = [plot_impulse_response,\n",
    "                      plot_spectral_density,\n",
    "                      plot_autocovariance,\n",
    "                      plot_simulation]\n",
    "    for plot_func, ax in zip(plot_functions, axes.flatten()):\n",
    "        plot_func(arma, ax)\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let’s call these functions to generate plots.\n",
    "\n",
    "As a warmup, let’s make sure things look right when we for the pure white noise model $ X_t = \\epsilon_t $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = 0.0\n",
    "θ = 0.0\n",
    "arma = qe.ARMA(ϕ, θ)\n",
    "quad_plot(arma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look carefully, things look good: the spectrum is the flat line at $ 10^0 $ at the very top of the spectrum graphs,\n",
    "which is at it should be.\n",
    "\n",
    "Also\n",
    "\n",
    "> - the variance  equals $ 1 = \\frac{1}{2 \\pi} \\int_{-\\pi}^\\pi 1 d \\omega $ as it should.  \n",
    "- the covariogram and impulse response look as they should.  \n",
    "- it is actually challenging to visualize a time series realization of white noise – a sequence of surprises – but this too looks pretty good.  \n",
    "\n",
    "\n",
    "\n",
    "To get some more examples, as our laboratory\n",
    "we’ll replicate quartets of graphs that [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012) use to teach “how to read spectral densities”.\n",
    "\n",
    "Ljunqvist and Sargent’s first  model is  $ X_t = 1.3 X_{t-1} - .7 X_{t-2} + \\epsilon_t $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = 1.3, -.7\n",
    "θ = 0.0\n",
    "arma = qe.ARMA(ϕ, θ)\n",
    "quad_plot(arma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ljungqvist and Sargent’s second model is $ X_t = .9 X_{t-1} + \\epsilon_t $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = 0.9\n",
    "θ = -0.0\n",
    "arma = qe.ARMA(ϕ, θ)\n",
    "quad_plot(arma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ljungqvist and Sargent’s third  model is  $ X_t = .8 X_{t-4} + \\epsilon_t $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = 0., 0., 0., .8\n",
    "θ = -0.0\n",
    "arma = qe.ARMA(ϕ, θ)\n",
    "quad_plot(arma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ljungqvist and Sargent’s fourth  model is  $ X_t = .98 X_{t-1}  + \\epsilon_t -.7 \\epsilon_{t-1} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ϕ = .98\n",
    "θ = -0.7\n",
    "arma = qe.ARMA(ϕ, θ)\n",
    "quad_plot(arma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explanation\n",
    "\n",
    "The call\n",
    "\n",
    "> `arma = ARMA(ϕ, θ, σ)`\n",
    "\n",
    "\n",
    "creates an instance `arma` that represents the ARMA($ p, q $) model\n",
    "\n",
    "$$\n",
    "X_t = \\phi_1 X_{t-1} + ... + \\phi_p X_{t-p} +\n",
    "    \\epsilon_t + \\theta_1 \\epsilon_{t-1} + ... + \\theta_q \\epsilon_{t-q}\n",
    "$$\n",
    "\n",
    "If `ϕ` and `θ` are arrays or sequences, then the interpretation will\n",
    "be\n",
    "\n",
    "- `ϕ` holds the vector of parameters $ (\\phi_1, \\phi_2,..., \\phi_p) $.  \n",
    "- `θ` holds the vector of parameters $ (\\theta_1, \\theta_2,..., \\theta_q) $.  \n",
    "\n",
    "\n",
    "The parameter `σ` is always a scalar, the standard deviation of the white noise.\n",
    "\n",
    "We also permit `ϕ` and `θ` to be scalars, in which case the model will be interpreted as\n",
    "\n",
    "$$\n",
    "X_t = \\phi X_{t-1} + \\epsilon_t + \\theta \\epsilon_{t-1}\n",
    "$$\n",
    "\n",
    "The two numerical packages most useful for working with ARMA models are `scipy.signal` and `numpy.fft`.\n",
    "\n",
    "The package `scipy.signal` expects the parameters to be passed into its functions in a manner consistent with the alternative ARMA notation [(8)](#equation-arma-lag1).\n",
    "\n",
    "For example, the impulse response sequence $ \\{\\psi_t\\} $ discussed above can be obtained using `scipy.signal.dimpulse`, and the function call should be of the form\n",
    "\n",
    "> `times, ψ = dimpulse((ma_poly, ar_poly, 1), n=impulse_length)`\n",
    "\n",
    "\n",
    "where `ma_poly` and `ar_poly` correspond to the polynomials in [(7)](#equation-arma-poly) — that is,\n",
    "\n",
    "- `ma_poly` is the vector $ (1, \\theta_1, \\theta_2, \\ldots, \\theta_q) $  \n",
    "- `ar_poly` is the vector $ (1, -\\phi_1, -\\phi_2, \\ldots, - \\phi_p) $  \n",
    "\n",
    "\n",
    "To this end, we also maintain the arrays `ma_poly` and `ar_poly` as instance data, with their values computed automatically from the values of `phi` and `theta` supplied by the user.\n",
    "\n",
    "If the user decides to change the value of either `theta` or `phi` ex-post by assignments\n",
    "such as `arma.phi = (0.5, 0.2)` or `arma.theta = (0, -0.1)`.\n",
    "\n",
    "then `ma_poly` and `ar_poly` should update automatically to reflect these new parameters.\n",
    "\n",
    "This is achieved in our implementation by using [descriptors](https://lectures.quantecon.org/python_advanced_features.html#descriptors)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing the Autocovariance Function\n",
    "\n",
    "As discussed above, for ARMA processes the spectral density has a [simple representation](#arma-spec-den) that is relatively easy to calculate.\n",
    "\n",
    "Given this fact, the easiest way to obtain the autocovariance function is to recover it from the spectral\n",
    "density via the inverse Fourier transform.\n",
    "\n",
    "Here we use NumPy’s Fourier transform package np.fft, which wraps a standard Fortran-based package called FFTPACK.\n",
    "\n",
    "A look at [the np.fft documentation](http://docs.scipy.org/doc/numpy/reference/routines.fft.html) shows that the inverse transform np.fft.ifft takes a given sequence $ A_0, A_1, \\ldots, A_{n-1} $ and\n",
    "returns the sequence $ a_0, a_1, \\ldots, a_{n-1} $ defined by\n",
    "\n",
    "$$\n",
    "a_k = \\frac{1}{n} \\sum_{t=0}^{n-1} A_t e^{ik 2\\pi t / n}\n",
    "$$\n",
    "\n",
    "Thus, if we set $ A_t = f(\\omega_t) $, where $ f $ is the spectral density and\n",
    "$ \\omega_t := 2 \\pi t / n $, then\n",
    "\n",
    "$$\n",
    "a_k\n",
    "= \\frac{1}{n} \\sum_{t=0}^{n-1} f(\\omega_t) e^{i \\omega_t k}\n",
    "= \\frac{1}{2\\pi} \\frac{2 \\pi}{n} \\sum_{t=0}^{n-1} f(\\omega_t) e^{i \\omega_t k},\n",
    "\\qquad\n",
    "\\omega_t := 2 \\pi t / n\n",
    "$$\n",
    "\n",
    "For $ n $ sufficiently large, we then have\n",
    "\n",
    "$$\n",
    "a_k\n",
    "\\approx \\frac{1}{2\\pi} \\int_0^{2 \\pi} f(\\omega) e^{i \\omega k} d \\omega\n",
    "= \\frac{1}{2\\pi} \\int_{-\\pi}^{\\pi} f(\\omega) e^{i \\omega k} d \\omega\n",
    "$$\n",
    "\n",
    "(You can check the last equality)\n",
    "\n",
    "In view of [(14)](#equation-ift), we have now shown that, for $ n $ sufficiently large, $ a_k \\approx \\gamma(k) $ — which is exactly what we want to compute."
   ]
  }
 ],
 "metadata": {
  "filename": "arma.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Covariance Stationary Processes"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}