{
"cells": [
{
"cell_type": "markdown",
"id": "4c4ce752",
"metadata": {},
"source": [
"$$\n",
"\\newcommand{\\argmax}{arg\\,max}\n",
"\\newcommand{\\argmin}{arg\\,min}\n",
"\\newcommand{\\col}{col}\n",
"\\newcommand{\\Span}{span}\n",
"\\newcommand{\\epsilon}{\\varepsilon}\n",
"\\newcommand{\\EE}{\\mathbb{E}}\n",
"\\newcommand{\\PP}{\\mathbb{P}}\n",
"\\newcommand{\\RR}{\\mathbb{R}}\n",
"\\newcommand{\\NN}{\\mathbb{N}}\n",
"\\newcommand{\\ZZ}{\\mathbb{Z}}\n",
"\\newcommand{\\aA}{\\mathcal{A}}\n",
"\\newcommand{\\bB}{\\mathcal{B}}\n",
"\\newcommand{\\cC}{\\mathcal{C}}\n",
"\\newcommand{\\dD}{\\mathcal{D}}\n",
"\\newcommand{\\eE}{\\mathcal{E}}\n",
"\\newcommand{\\fF}{\\mathcal{F}}\n",
"\\newcommand{\\gG}{\\mathcal{G}}\n",
"\\newcommand{\\hH}{\\mathcal{H}}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "a3222374",
"metadata": {},
"source": [
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "396169cf",
"metadata": {},
"source": [
"# Additive and Multiplicative Functionals\n",
"\n",
"\n",
"\n",
"In addition to what’s in Anaconda, this lecture will need the following libraries:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d4ce2c16",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"!pip install --upgrade quantecon"
]
},
{
"cell_type": "markdown",
"id": "99ff2a73",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"Many economic time series display persistent growth that prevents them from being asymptotically stationary and ergodic.\n",
"\n",
"For example, outputs, prices, and dividends typically display irregular but persistent growth.\n",
"\n",
"Asymptotic stationarity and ergodicity are key assumptions needed to make it possible to learn by applying statistical methods.\n",
"\n",
"But there are good ways to model time series that have persistent growth that still enable statistical learning based on a law of large numbers for an asymptotically stationary and ergodic process.\n",
"\n",
"Thus, [[Hansen, 2012](https://python-advanced.quantecon.org/zreferences.html#id246)] described two classes of time series models that accommodate growth.\n",
"\n",
"They are\n",
"\n",
"1. **additive functionals** that display random “arithmetic growth” \n",
"1. **multiplicative functionals** that display random “geometric growth” \n",
"\n",
"\n",
"These two classes of processes are closely connected.\n",
"\n",
"If a process $ \\{y_t\\} $ is an additive functional and $ \\phi_t = \\exp(y_t) $, then $ \\{\\phi_t\\} $ is a multiplicative functional.\n",
"\n",
"In this lecture, we describe both additive functionals and multiplicative functionals.\n",
"\n",
"We also describe and compute decompositions of additive and multiplicative processes into four components:\n",
"\n",
"1. a **constant** \n",
"1. a **trend** component \n",
"1. an asymptotically **stationary** component \n",
"1. a **martingale** \n",
"\n",
"\n",
"We describe how to construct, simulate, and interpret these components.\n",
"\n",
"More details about these concepts and algorithms can be found in Hansen [[Hansen, 2012](https://python-advanced.quantecon.org/zreferences.html#id246)] and Hansen and Sargent [[Hansen and Sargent, 2024](https://python-advanced.quantecon.org/zreferences.html#id247)].\n",
"\n",
"Let’s start with some imports:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f34ad80c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.linalg as la\n",
"import quantecon as qe\n",
"import matplotlib.pyplot as plt\n",
"from scipy.stats import norm, lognorm"
]
},
{
"cell_type": "markdown",
"id": "812c799c",
"metadata": {},
"source": [
"## A particular additive functional\n",
"\n",
"[[Hansen, 2012](https://python-advanced.quantecon.org/zreferences.html#id246)] describes a general class of additive functionals.\n",
"\n",
"This lecture focuses on a subclass of these: a scalar process $ \\{y_t\\}_{t=0}^\\infty $ whose increments are driven by a Gaussian vector autoregression.\n",
"\n",
"Our special additive functional displays interesting time series behavior while also being easy to construct, simulate, and analyze\n",
"by using linear state-space tools.\n",
"\n",
"We construct our additive functional from two pieces, the first of which is a **first-order vector autoregression** (VAR)\n",
"\n",
"\n",
"\n",
"$$\n",
"x_{t+1} = A x_t + B z_{t+1} \\tag{32.1}\n",
"$$\n",
"\n",
"Here\n",
"\n",
"- $ x_t $ is an $ n \\times 1 $ vector, \n",
"- $ A $ is an $ n \\times n $ stable matrix (all eigenvalues lie within the open unit circle), \n",
"- $ z_{t+1} \\sim {\\cal N}(0,I) $ is an $ m \\times 1 $ IID shock, \n",
"- $ B $ is an $ n \\times m $ matrix, and \n",
"- $ x_0 \\sim {\\cal N}(\\mu_0, \\Sigma_0) $ is a random initial condition for $ x $ \n",
"\n",
"\n",
"The second piece is an equation that expresses increments\n",
"of $ \\{y_t\\}_{t=0}^\\infty $ as linear functions of\n",
"\n",
"- a scalar constant $ \\nu $, \n",
"- the vector $ x_t $, and \n",
"- the same Gaussian vector $ z_{t+1} $ that appears in the VAR [(32.1)](#equation-old1-additive-functionals) \n",
"\n",
"\n",
"In particular,\n",
"\n",
"\n",
"\n",
"$$\n",
"y_{t+1} - y_{t} = \\nu + D x_{t} + F z_{t+1} \\tag{32.2}\n",
"$$\n",
"\n",
"Here $ y_0 \\sim {\\cal N}(\\mu_{y0}, \\Sigma_{y0}) $ is a random\n",
"initial condition for $ y $.\n",
"\n",
"The nonstationary random process $ \\{y_t\\}_{t=0}^\\infty $ displays\n",
"systematic but random *arithmetic growth*."
]
},
{
"cell_type": "markdown",
"id": "e0757182",
"metadata": {},
"source": [
"### Linear state-space representation\n",
"\n",
"A convenient way to represent our additive functional is to use a [linear state space system](https://python-intro.quantecon.org/linear_models.html).\n",
"\n",
"To do this, we set up state and observation vectors\n",
"\n",
"$$\n",
"\\hat{x}_t = \\begin{bmatrix} 1 \\\\ x_t \\\\ y_t \\end{bmatrix}\n",
"\\quad \\text{and} \\quad\n",
"\\hat{y}_t = \\begin{bmatrix} x_t \\\\ y_t \\end{bmatrix}\n",
"$$\n",
"\n",
"Next we construct a linear system\n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
" 1 \\\\\n",
" x_{t+1} \\\\\n",
" y_{t+1}\n",
" \\end{bmatrix} =\n",
" \\begin{bmatrix}\n",
" 1 & 0 & 0 \\\\\n",
" 0 & A & 0 \\\\\n",
" \\nu & D & 1\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
" 1 \\\\\n",
" x_t \\\\\n",
" y_t\n",
"\\end{bmatrix} +\n",
"\\begin{bmatrix}\n",
" 0 \\\\ B \\\\ F\n",
"\\end{bmatrix}\n",
"z_{t+1}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
" x_t \\\\\n",
" y_t\n",
"\\end{bmatrix}\n",
"= \\begin{bmatrix}\n",
" 0 & I & 0 \\\\\n",
" 0 & 0 & 1\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
" 1 \\\\ x_t \\\\ y_t\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"This can be written as\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\hat{x}_{t+1} &= \\hat{A} \\hat{x}_t + \\hat{B} z_{t+1} \\\\\n",
" \\hat{y}_{t} &= \\hat{D} \\hat{x}_t\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"which is a standard linear state space system.\n",
"\n",
"To study it, we could map it into an instance of [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) from [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
"\n",
"But here we will use a different set of code for simulation, for reasons described below."
]
},
{
"cell_type": "markdown",
"id": "fb70e327",
"metadata": {},
"source": [
"## Dynamics\n",
"\n",
"Let’s run some simulations to build intuition.\n",
"\n",
"\n",
"\n",
"In doing so we’ll assume that $ z_{t+1} $ is scalar and that $ \\tilde x_t $ follows a 4th-order scalar autoregression.\n",
"\n",
"\n",
"\n",
"$$\n",
"\\tilde x_{t+1} = \\phi_1 \\tilde x_{t} + \\phi_2 \\tilde x_{t-1} +\n",
"\\phi_3 \\tilde x_{t-2} +\n",
"\\phi_4 \\tilde x_{t-3} + \\sigma z_{t+1} \\tag{32.3}\n",
"$$\n",
"\n",
"in which the zeros $ z $ of the polynomial\n",
"\n",
"$$\n",
"\\phi(z) = ( 1 - \\phi_1 z - \\phi_2 z^2 - \\phi_3 z^3 - \\phi_4 z^4 )\n",
"$$\n",
"\n",
"are strictly greater than unity in absolute value.\n",
"\n",
"(Being a zero of $ \\phi(z) $ means that $ \\phi(z) = 0 $)\n",
"\n",
"Let the increment in $ \\{y_t\\} $ obey\n",
"\n",
"$$\n",
"y_{t+1} - y_t = \\nu + \\tilde x_t + \\sigma z_{t+1}\n",
"$$\n",
"\n",
"with an initial condition for $ y_0 $.\n",
"\n",
"While [(32.3)](#equation-ftaf) is not a first order system like [(32.1)](#equation-old1-additive-functionals), we know that it can be mapped into a first order system.\n",
"\n",
"- For an example of such a mapping, see [this example](https://python.quantecon.org/linear_models.html#second-order-difference-equation). \n",
"\n",
"\n",
"In fact, this whole model can be mapped into the additive functional system definition in [(32.1)](#equation-old1-additive-functionals) – [(32.2)](#equation-old2-additive-functionals) by appropriate selection of the matrices $ A, B, D, F $.\n",
"\n",
"You can try writing these matrices down now as an exercise — correct expressions appear in the code below."
]
},
{
"cell_type": "markdown",
"id": "299c77d3",
"metadata": {},
"source": [
"### Simulation\n",
"\n",
"When simulating we embed our variables into a bigger system.\n",
"\n",
"This system also constructs the components of the decompositions of $ y_t $ and of $ \\exp(y_t) $ proposed by Hansen [[Hansen, 2012](https://python-advanced.quantecon.org/zreferences.html#id246)].\n",
"\n",
"All of these objects are computed using the code below\n",
"\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3371f11c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"class AMF_LSS_VAR:\n",
" \"\"\"\n",
" This class transforms an additive (multiplicative)\n",
" functional into a QuantEcon linear state space system.\n",
" Handles both matrix and scalar inputs.\n",
" \"\"\"\n",
"\n",
" def __init__(self, A, B, D, F=None, ν=None):\n",
" # Handle scalar inputs by converting to arrays\n",
" self.scalar_case = np.isscalar(A)\n",
" if self.scalar_case:\n",
" A = np.atleast_2d(A)\n",
" B = np.atleast_2d(B)\n",
" D = np.atleast_1d(D)\n",
" if F is not None:\n",
" F = np.atleast_2d(F)\n",
" if ν is not None and np.isscalar(ν):\n",
" ν = float(ν)\n",
"\n",
" # Unpack required elements\n",
" self.nx, self.nk = B.shape\n",
" self.A, self.B = A, B\n",
"\n",
" # Checking the dimension of D (extended from the scalar case)\n",
" if len(D.shape) > 1 and D.shape[0] != 1:\n",
" self.nm = D.shape[0]\n",
" self.D = D\n",
" elif len(D.shape) > 1 and D.shape[0] == 1:\n",
" self.nm = 1\n",
" self.D = D\n",
" else:\n",
" self.nm = 1\n",
" self.D = np.expand_dims(D, 0)\n",
"\n",
" # Create space for additive decomposition\n",
" self.add_decomp = None\n",
" self.mult_decomp = None\n",
"\n",
" # Set F\n",
" if not np.any(F):\n",
" self.F = np.zeros((self.nk, 1))\n",
" else:\n",
" self.F = F\n",
"\n",
" # Set ν\n",
" if not np.any(ν):\n",
" self.ν = np.zeros((self.nm, 1))\n",
" elif type(ν) == float:\n",
" self.ν = np.asarray([[ν]])\n",
" elif len(ν.shape) == 1:\n",
" self.ν = np.expand_dims(ν, 1)\n",
" else:\n",
" self.ν = ν\n",
"\n",
" if self.ν.shape[0] != self.D.shape[0]:\n",
" raise ValueError(\"The dimension of ν is inconsistent with D!\")\n",
"\n",
" # Construct BIG state space representation\n",
" self.lss = self.construct_ss()\n",
"\n",
" def construct_ss(self):\n",
" \"\"\"\n",
" This creates the state space representation that can be passed\n",
" into the quantecon LSS class.\n",
" \"\"\"\n",
" # Pull out useful info\n",
" nx, nk, nm = self.nx, self.nk, self.nm\n",
" A, B, D, F, ν = self.A, self.B, self.D, self.F, self.ν\n",
" if self.add_decomp:\n",
" ν, H, g = self.add_decomp\n",
" else:\n",
" ν, H, g = self.additive_decomp()\n",
"\n",
" # Auxiliary blocks with 0's and 1's to fill out the lss matrices\n",
" nx0c = np.zeros((nx, 1))\n",
" nx0r = np.zeros(nx)\n",
" nx1 = np.ones(nx)\n",
" nk0 = np.zeros(nk)\n",
" ny0c = np.zeros((nm, 1))\n",
" ny0r = np.zeros(nm)\n",
" ny1m = np.eye(nm)\n",
" ny0m = np.zeros((nm, nm))\n",
" nyx0m = np.zeros_like(D)\n",
"\n",
" # Build A matrix for LSS\n",
" # Order of states is: [1, t, xt, yt, mt]\n",
" A1 = np.hstack([1, 0, nx0r, ny0r, ny0r]) # Transition for 1\n",
" A2 = np.hstack([1, 1, nx0r, ny0r, ny0r]) # Transition for t\n",
" # Transition for x_{t+1}\n",
" A3 = np.hstack([nx0c, nx0c, A, nyx0m.T, nyx0m.T])\n",
" # Transition for y_{t+1}\n",
" A4 = np.hstack([ν, ny0c, D, ny1m, ny0m])\n",
" # Transition for m_{t+1}\n",
" A5 = np.hstack([ny0c, ny0c, nyx0m, ny0m, ny1m])\n",
" Abar = np.vstack([A1, A2, A3, A4, A5])\n",
"\n",
" # Build B matrix for LSS\n",
" Bbar = np.vstack([nk0, nk0, B, F, H])\n",
"\n",
" # Build G matrix for LSS\n",
" # Order of observation is: [xt, yt, mt, st, tt]\n",
" # Selector for x_{t}\n",
" G1 = np.hstack([nx0c, nx0c, np.eye(nx), nyx0m.T, nyx0m.T])\n",
" G2 = np.hstack([ny0c, ny0c, nyx0m, ny1m, ny0m]) # Selector for y_{t}\n",
" # Selector for martingale\n",
" G3 = np.hstack([ny0c, ny0c, nyx0m, ny0m, ny1m])\n",
" G4 = np.hstack([ny0c, ny0c, -g, ny0m, ny0m]) # Selector for stationary\n",
" G5 = np.hstack([ny0c, ν, nyx0m, ny0m, ny0m]) # Selector for trend\n",
" Gbar = np.vstack([G1, G2, G3, G4, G5])\n",
"\n",
" # Build H matrix for LSS\n",
" Hbar = np.zeros((Gbar.shape[0], nk))\n",
"\n",
" # Build LSS type\n",
" x0 = np.hstack([1, 0, nx0r, ny0r, ny0r])\n",
" S0 = np.zeros((len(x0), len(x0)))\n",
" lss = qe.LinearStateSpace(Abar, Bbar, Gbar, Hbar, mu_0=x0, Sigma_0=S0)\n",
"\n",
" return lss\n",
"\n",
" def additive_decomp(self):\n",
" \"\"\"\n",
" Return values for the martingale decomposition\n",
" - ν : unconditional mean difference in Y\n",
" - H : coefficient for the (linear) martingale component (κ_a)\n",
" - g : coefficient for the stationary component g(x)\n",
" - Y_0 : it should be the function of X_0 (for now set it to 0.0)\n",
" \"\"\"\n",
" I = np.identity(self.nx)\n",
" A_res = la.solve(I - self.A, I)\n",
" g = self.D @ A_res\n",
" H = self.F + self.D @ A_res @ self.B\n",
"\n",
" return self.ν, H, g\n",
"\n",
" def multiplicative_decomp(self):\n",
" \"\"\"\n",
" Return values for the multiplicative decomposition (Example 5.4.4.)\n",
" - ν_tilde : eigenvalue\n",
" - H : vector for the Jensen term\n",
" \"\"\"\n",
" ν, H, g = self.additive_decomp()\n",
" ν_tilde = ν + (.5)*np.expand_dims(np.diag(H @ H.T), 1)\n",
"\n",
" return ν_tilde, H, g\n",
"\n",
" def loglikelihood_path(self, x, y):\n",
" A, B, D, F = self.A, self.B, self.D, self.F\n",
" k, T = y.shape\n",
" FF = F @ F.T\n",
" FFinv = la.inv(FF)\n",
" temp = y[:, 1:] - y[:, :-1] - D @ x[:, :-1]\n",
" obs = temp * FFinv * temp\n",
" obssum = np.cumsum(obs)\n",
" scalar = (np.log(la.det(FF)) + k*np.log(2*np.pi))*np.arange(1, T)\n",
"\n",
" return -(.5)*(obssum + scalar)\n",
"\n",
" def loglikelihood(self, x, y):\n",
" llh = self.loglikelihood_path(x, y)\n",
"\n",
" return llh[-1]"
]
},
{
"cell_type": "markdown",
"id": "a766e8fa",
"metadata": {},
"source": [
"#### Plotting\n",
"\n",
"The code below adds some functions that generate plots for instances of the `AMF_LSS_VAR` [class](#amf-lss)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1420460f",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def plot_given_paths(amf, T, ypath, mpath, spath, tpath,\n",
" mbounds, sbounds, horline=0, show_trend=True):\n",
"\n",
" # Allocate space\n",
" trange = np.arange(T)\n",
"\n",
" # Create figure\n",
" fig, ax = plt.subplots(2, 2, sharey=True, figsize=(15, 8))\n",
"\n",
" # Plot all paths together\n",
" ax[0, 0].plot(trange, ypath[0, :], label=\"$y_t$\", color=\"k\")\n",
" ax[0, 0].plot(trange, mpath[0, :], label=\"$m_t$\", color=\"m\")\n",
" ax[0, 0].plot(trange, spath[0, :], label=\"$s_t$\", color=\"g\")\n",
" if show_trend:\n",
" ax[0, 0].plot(trange, tpath[0, :], label=\"$t_t$\", color=\"r\")\n",
" ax[0, 0].axhline(horline, color=\"k\", linestyle=\"-.\")\n",
" ax[0, 0].set_title(\"One Path of All Variables\")\n",
" ax[0, 0].legend(loc=\"upper left\")\n",
"\n",
" # Plot Martingale Component\n",
" ax[0, 1].plot(trange, mpath[0, :], \"m\")\n",
" ax[0, 1].plot(trange, mpath.T, alpha=0.45, color=\"m\")\n",
" ub = mbounds[1, :]\n",
" lb = mbounds[0, :]\n",
"\n",
" ax[0, 1].fill_between(trange, lb, ub, alpha=0.25, color=\"m\")\n",
" ax[0, 1].set_title(\"Martingale Components for Many Paths\")\n",
" ax[0, 1].axhline(horline, color=\"k\", linestyle=\"-.\")\n",
"\n",
" # Plot Stationary Component\n",
" ax[1, 0].plot(spath[0, :], color=\"g\")\n",
" ax[1, 0].plot(spath.T, alpha=0.25, color=\"g\")\n",
" ub = sbounds[1, :]\n",
" lb = sbounds[0, :]\n",
" ax[1, 0].fill_between(trange, lb, ub, alpha=0.25, color=\"g\")\n",
" ax[1, 0].axhline(horline, color=\"k\", linestyle=\"-.\")\n",
" ax[1, 0].set_title(\"Stationary Components for Many Paths\")\n",
"\n",
" # Plot Trend Component\n",
" if show_trend:\n",
" ax[1, 1].plot(tpath.T, color=\"r\")\n",
" ax[1, 1].set_title(\"Trend Components for Many Paths\")\n",
" ax[1, 1].axhline(horline, color=\"k\", linestyle=\"-.\")\n",
"\n",
" return fig\n",
"\n",
"def plot_additive(amf, T, npaths=25, show_trend=True):\n",
" \"\"\"\n",
" Plots for the additive decomposition.\n",
" Acts on an instance amf of the AMF_LSS_VAR class\n",
"\n",
" \"\"\"\n",
" # Pull out right sizes so we know how to increment\n",
" nx, nk, nm = amf.nx, amf.nk, amf.nm\n",
"\n",
" # Allocate space (nm is the number of additive functionals -\n",
" # we want npaths for each)\n",
" mpath = np.empty((nm*npaths, T))\n",
" mbounds = np.empty((nm*2, T))\n",
" spath = np.empty((nm*npaths, T))\n",
" sbounds = np.empty((nm*2, T))\n",
" tpath = np.empty((nm*npaths, T))\n",
" ypath = np.empty((nm*npaths, T))\n",
"\n",
" # Simulate for as long as we wanted\n",
" moment_generator = amf.lss.moment_sequence()\n",
" # Pull out population moments\n",
" for t in range (T):\n",
" tmoms = next(moment_generator)\n",
" ymeans = tmoms[1]\n",
" yvar = tmoms[3]\n",
"\n",
" # Lower and upper bounds - for each additive functional\n",
" for ii in range(nm):\n",
" li, ui = ii*2, (ii+1)*2\n",
" mscale = np.sqrt(yvar[nx+nm+ii, nx+nm+ii])\n",
" sscale = np.sqrt(yvar[nx+2*nm+ii, nx+2*nm+ii])\n",
" if mscale == 0.0:\n",
" mscale = 1e-12 # avoids a RuntimeWarning from calculating ppf\n",
" if sscale == 0.0: # of normal distribution with std dev = 0.\n",
" sscale = 1e-12 # sets std dev to small value instead\n",
"\n",
" madd_dist = norm(ymeans[nx+nm+ii], mscale)\n",
" sadd_dist = norm(ymeans[nx+2*nm+ii], sscale)\n",
"\n",
" mbounds[li:ui, t] = madd_dist.ppf([0.01, .99])\n",
" sbounds[li:ui, t] = sadd_dist.ppf([0.01, .99])\n",
"\n",
" # Pull out paths\n",
" for n in range(npaths):\n",
" x, y = amf.lss.simulate(T)\n",
" for ii in range(nm):\n",
" ypath[npaths*ii+n, :] = y[nx+ii, :]\n",
" mpath[npaths*ii+n, :] = y[nx+nm + ii, :]\n",
" spath[npaths*ii+n, :] = y[nx+2*nm + ii, :]\n",
" tpath[npaths*ii+n, :] = y[nx+3*nm + ii, :]\n",
"\n",
" add_figs = []\n",
"\n",
" for ii in range(nm):\n",
" li, ui = npaths*(ii), npaths*(ii+1)\n",
" LI, UI = 2*(ii), 2*(ii+1)\n",
" add_figs.append(plot_given_paths(amf, T,\n",
" ypath[li:ui,:],\n",
" mpath[li:ui,:],\n",
" spath[li:ui,:],\n",
" tpath[li:ui,:],\n",
" mbounds[LI:UI,:],\n",
" sbounds[LI:UI,:],\n",
" show_trend=show_trend))\n",
"\n",
" add_figs[ii].suptitle(f'Additive decomposition of $y_{ii+1}$',\n",
" fontsize=14)\n",
"\n",
" return add_figs\n",
"\n",
"\n",
"def plot_multiplicative(amf, T, npaths=25, show_trend=True):\n",
" \"\"\"\n",
" Plots for the multiplicative decomposition\n",
"\n",
" \"\"\"\n",
" # Pull out right sizes so we know how to increment\n",
" nx, nk, nm = amf.nx, amf.nk, amf.nm\n",
" # Matrices for the multiplicative decomposition\n",
" ν_tilde, H, g = amf.multiplicative_decomp()\n",
"\n",
" # Allocate space (nm is the number of functionals -\n",
" # we want npaths for each)\n",
" mpath_mult = np.empty((nm*npaths, T))\n",
" mbounds_mult = np.empty((nm*2, T))\n",
" spath_mult = np.empty((nm*npaths, T))\n",
" sbounds_mult = np.empty((nm*2, T))\n",
" tpath_mult = np.empty((nm*npaths, T))\n",
" ypath_mult = np.empty((nm*npaths, T))\n",
"\n",
" # Simulate for as long as we wanted\n",
" moment_generator = amf.lss.moment_sequence()\n",
" # Pull out population moments\n",
" for t in range(T):\n",
" tmoms = next(moment_generator)\n",
" ymeans = tmoms[1]\n",
" yvar = tmoms[3]\n",
"\n",
" # Lower and upper bounds - for each multiplicative functional\n",
" for ii in range(nm):\n",
" li, ui = ii*2, (ii+1)*2\n",
" Mdist = lognorm(np.sqrt(yvar[nx+nm+ii, nx+nm+ii]).item(),\n",
" scale=np.exp(ymeans[nx+nm+ii] \\\n",
" - t * (.5)\n",
" * np.expand_dims(\n",
" np.diag(H @ H.T),\n",
" 1\n",
" )[ii]\n",
" ).item()\n",
" )\n",
" Sdist = lognorm(np.sqrt(yvar[nx+2*nm+ii, nx+2*nm+ii]).item(),\n",
" scale = np.exp(-ymeans[nx+2*nm+ii]).item())\n",
" mbounds_mult[li:ui, t] = Mdist.ppf([.01, .99])\n",
" sbounds_mult[li:ui, t] = Sdist.ppf([.01, .99])\n",
"\n",
" # Pull out paths\n",
" for n in range(npaths):\n",
" x, y = amf.lss.simulate(T)\n",
" for ii in range(nm):\n",
" ypath_mult[npaths*ii+n, :] = np.exp(y[nx+ii, :])\n",
" mpath_mult[npaths*ii+n, :] = np.exp(y[nx+nm + ii, :] \\\n",
" - np.arange(T)*(.5)\n",
" * np.expand_dims(np.diag(H\n",
" @ H.T),\n",
" 1)[ii]\n",
" )\n",
" spath_mult[npaths*ii+n, :] = 1/np.exp(-y[nx+2*nm + ii, :])\n",
" tpath_mult[npaths*ii+n, :] = np.exp(y[nx+3*nm + ii, :]\n",
" + np.arange(T)*(.5)\n",
" * np.expand_dims(np.diag(H\n",
" @ H.T),\n",
" 1)[ii]\n",
" )\n",
"\n",
" mult_figs = []\n",
"\n",
" for ii in range(nm):\n",
" li, ui = npaths*(ii), npaths*(ii+1)\n",
" LI, UI = 2*(ii), 2*(ii+1)\n",
"\n",
" mult_figs.append(plot_given_paths(amf,T,\n",
" ypath_mult[li:ui,:],\n",
" mpath_mult[li:ui,:],\n",
" spath_mult[li:ui,:],\n",
" tpath_mult[li:ui,:],\n",
" mbounds_mult[LI:UI,:],\n",
" sbounds_mult[LI:UI,:],\n",
" 1,\n",
" show_trend=show_trend))\n",
" mult_figs[ii].suptitle(\n",
" f'Multiplicative decomposition of $y_{ii+1}$', fontsize=14)\n",
"\n",
" return mult_figs\n",
"\n",
"def plot_martingale_paths(amf, T, mpath, mbounds, \n",
" horline=1, show_trend=False):\n",
"\n",
" # Allocate space\n",
" trange = np.arange(T)\n",
"\n",
" # Create figure\n",
" fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n",
"\n",
" # Plot Martingale Component\n",
" ub = mbounds[1, :]\n",
" lb = mbounds[0, :]\n",
" ax.fill_between(trange, lb, ub, color=\"#ffccff\")\n",
" ax.axhline(horline, color=\"k\", linestyle=\"-.\")\n",
" ax.plot(trange, mpath.T, linewidth=0.25, color=\"#4c4c4c\")\n",
"\n",
" return fig\n",
"\n",
"def plot_martingales(amf, T, npaths=25):\n",
"\n",
" # Pull out right sizes so we know how to increment\n",
" nx, nk, nm = amf.nx, amf.nk, amf.nm\n",
" # Matrices for the multiplicative decomposition\n",
" ν_tilde, H, g = amf.multiplicative_decomp()\n",
"\n",
" # Allocate space (nm is the number of functionals -\n",
" # we want npaths for each)\n",
" mpath_mult = np.empty((nm*npaths, T))\n",
" mbounds_mult = np.empty((nm*2, T))\n",
"\n",
" # Simulate for as long as we wanted\n",
" moment_generator = amf.lss.moment_sequence()\n",
" # Pull out population moments\n",
" for t in range (T):\n",
" tmoms = next(moment_generator)\n",
" ymeans = tmoms[1]\n",
" yvar = tmoms[3]\n",
"\n",
" # Lower and upper bounds - for each functional\n",
" for ii in range(nm):\n",
" li, ui = ii*2, (ii+1)*2\n",
" Mdist = lognorm(np.sqrt(yvar[nx+nm+ii, nx+nm+ii]).item(),\n",
" scale= np.exp(ymeans[nx+nm+ii] \\\n",
" - t * (.5)\n",
" * np.expand_dims(\n",
" np.diag(H @ H.T),\n",
" 1)[ii]\n",
"\n",
" ).item()\n",
" )\n",
" mbounds_mult[li:ui, t] = Mdist.ppf([.01, .99])\n",
"\n",
" # Pull out paths\n",
" for n in range(npaths):\n",
" x, y = amf.lss.simulate(T)\n",
" for ii in range(nm):\n",
" mpath_mult[npaths*ii+n, :] = np.exp(y[nx+nm + ii, :] \\\n",
" - np.arange(T) * (.5)\n",
" * np.expand_dims(np.diag(H\n",
" @ H.T),\n",
" 1)[ii]\n",
" )\n",
"\n",
" mart_figs = []\n",
"\n",
" for ii in range(nm):\n",
" li, ui = npaths*(ii), npaths*(ii+1)\n",
" LI, UI = 2*(ii), 2*(ii+1)\n",
" mart_figs.append(plot_martingale_paths(amf, T, mpath_mult[li:ui, :],\n",
" mbounds_mult[LI:UI, :],\n",
" horline=1))\n",
" mart_figs[ii].suptitle(f'Martingale components for many paths of \\\n",
" $y_{ii+1}$', fontsize=14)\n",
"\n",
" return mart_figs"
]
},
{
"cell_type": "markdown",
"id": "fe23cccd",
"metadata": {},
"source": [
"For now, we just plot $ y_t $ and $ x_t $, postponing until later a description of exactly how we compute them.\n",
"\n",
"\n",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "59f031a3",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"ϕ_1, ϕ_2, ϕ_3, ϕ_4 = 0.5, -0.2, 0, 0.5\n",
"σ = 0.01\n",
"ν = 0.01 # Growth rate\n",
"\n",
"# A matrix should be n x n\n",
"A = np.array([[ϕ_1, ϕ_2, ϕ_3, ϕ_4],\n",
" [ 1, 0, 0, 0],\n",
" [ 0, 1, 0, 0],\n",
" [ 0, 0, 1, 0]])\n",
"\n",
"# B matrix should be n x k\n",
"B = np.array([[σ, 0, 0, 0]]).T\n",
"\n",
"D = np.array([1, 0, 0, 0]) @ A\n",
"F = np.array([1, 0, 0, 0]) @ B\n",
"\n",
"amf = AMF_LSS_VAR(A, B, D, F, ν=ν)\n",
"\n",
"T = 150\n",
"x, y = amf.lss.simulate(T)\n",
"\n",
"fig, ax = plt.subplots(2, 1, figsize=(10, 9))\n",
"\n",
"ax[0].plot(np.arange(T), y[amf.nx, :], color='k')\n",
"ax[0].set_title('Path of $y_t$')\n",
"ax[1].plot(np.arange(T), y[0, :], color='g')\n",
"ax[1].axhline(0, color='k', linestyle='-.')\n",
"ax[1].set_title('Associated path of $x_t$')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "07e0de50",
"metadata": {},
"source": [
"Notice the irregular but persistent growth in $ y_t $."
]
},
{
"cell_type": "markdown",
"id": "585e3391",
"metadata": {},
"source": [
"### Decomposition\n",
"\n",
"Hansen and Sargent [[Hansen and Sargent, 2024](https://python-advanced.quantecon.org/zreferences.html#id247)] describe how to construct a decomposition of\n",
"an additive functional into four parts:\n",
"\n",
"- a constant inherited from initial values $ x_0 $ and $ y_0 $ \n",
"- a linear trend \n",
"- a martingale \n",
"- an (asymptotically) stationary component \n",
"\n",
"\n",
"To attain this decomposition for the particular class of additive\n",
"functionals defined by [(32.1)](#equation-old1-additive-functionals) and [(32.2)](#equation-old2-additive-functionals), we first construct the matrices\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" H & := F + D (I - A)^{-1} B\n",
" \\\\\n",
" g & := D (I - A)^{-1}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"Then the Hansen [[Hansen, 2012](https://python-advanced.quantecon.org/zreferences.html#id246)], [[Hansen and Sargent, 2024](https://python-advanced.quantecon.org/zreferences.html#id247)] decomposition is\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" y_t\n",
" &= \\underbrace{t \\nu}_{\\text{trend component}} +\n",
" \\overbrace{\\sum_{j=1}^t H z_j}^{\\text{Martingale component}} -\n",
" \\underbrace{g x_t}_{\\text{stationary component}} +\n",
" \\overbrace{g x_0 + y_0}^{\\text{initial conditions}}\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"At this stage, you should pause and verify that $ y_{t+1} - y_t $ satisfies [(32.2)](#equation-old2-additive-functionals).\n",
"\n",
"It is convenient for us to introduce the following notation:\n",
"\n",
"- $ \\tau_t = \\nu t $ , a linear, deterministic trend \n",
"- $ m_t = \\sum_{j=1}^t H z_j $, a martingale with time $ t+1 $ increment $ H z_{t+1} $ \n",
"- $ s_t = g x_t $, an (asymptotically) stationary component \n",
"\n",
"\n",
"We want to characterize and simulate components $ \\tau_t, m_t, s_t $ of the decomposition.\n",
"\n",
"A convenient way to do this is to construct an appropriate instance of a [linear state space system](https://python-intro.quantecon.org/linear_models.html) by using [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) from [QuantEcon.py](http://quantecon.org/quantecon-py).\n",
"\n",
"This will allow us to use the routines in [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) to study dynamics.\n",
"\n",
"To start, observe that, under the dynamics in [(32.1)](#equation-old1-additive-functionals) and [(32.2)](#equation-old2-additive-functionals) and with the\n",
"definitions just given,\n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
" 1 \\\\\n",
" t+1 \\\\\n",
" x_{t+1} \\\\\n",
" y_{t+1} \\\\\n",
" m_{t+1}\n",
"\\end{bmatrix} =\n",
"\\begin{bmatrix}\n",
" 1 & 0 & 0 & 0 & 0 \\\\\n",
" 1 & 1 & 0 & 0 & 0 \\\\\n",
" 0 & 0 & A & 0 & 0 \\\\\n",
" \\nu & 0 & D & 1 & 0 \\\\\n",
" 0 & 0 & 0 & 0 & 1\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
" 1 \\\\\n",
" t \\\\\n",
" x_t \\\\\n",
" y_t \\\\\n",
" m_t\n",
"\\end{bmatrix} +\n",
"\\begin{bmatrix}\n",
" 0 \\\\\n",
" 0 \\\\\n",
" B \\\\\n",
" F \\\\\n",
" H\n",
"\\end{bmatrix}\n",
"z_{t+1}\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\begin{bmatrix}\n",
" x_t \\\\\n",
" y_t \\\\\n",
" \\tau_t \\\\\n",
" m_t \\\\\n",
" s_t\n",
"\\end{bmatrix} =\n",
"\\begin{bmatrix}\n",
" 0 & 0 & I & 0 & 0 \\\\\n",
" 0 & 0 & 0 & 1 & 0 \\\\\n",
" 0 & \\nu & 0 & 0 & 0 \\\\\n",
" 0 & 0 & 0 & 0 & 1 \\\\\n",
" 0 & 0 & -g & 0 & 0\n",
"\\end{bmatrix}\n",
"\\begin{bmatrix}\n",
" 1 \\\\\n",
" t \\\\\n",
" x_t \\\\\n",
" y_t \\\\\n",
" m_t\n",
"\\end{bmatrix}\n",
"$$\n",
"\n",
"With\n",
"\n",
"$$\n",
"\\tilde{x} := \\begin{bmatrix} 1 \\\\ t \\\\ x_t \\\\ y_t \\\\ m_t \\end{bmatrix}\n",
"\\quad \\text{and} \\quad\n",
"\\tilde{y} := \\begin{bmatrix} x_t \\\\ y_t \\\\ \\tau_t \\\\ m_t \\\\ s_t \\end{bmatrix}\n",
"$$\n",
"\n",
"we can write this as the linear state space system\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
" \\tilde{x}_{t+1} &= \\tilde{A} \\tilde{x}_t + \\tilde{B} z_{t+1} \\\\\n",
" \\tilde{y}_{t} &= \\tilde{D} \\tilde{x}_t\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"By picking out components of $ \\tilde y_t $, we can track all variables of\n",
"interest."
]
},
{
"cell_type": "markdown",
"id": "e8a1483f",
"metadata": {},
"source": [
"## Code\n",
"\n",
"The class `AMF_LSS_VAR` mentioned [above](#amf-lss) does all that we want to study our additive functional.\n",
"\n",
"In fact, `AMF_LSS_VAR` does more\n",
"because it allows us to study an associated multiplicative functional as well.\n",
"\n",
"(A hint that it does more is the name of the class – here AMF stands for\n",
"“additive and multiplicative functional” – the code computes and displays objects associated with\n",
"multiplicative functionals too.)\n",
"\n",
"Let’s use this code (embedded above) to explore the [example process described above](#addfunc-eg1).\n",
"\n",
"If you run [the code that first simulated that example](#addfunc-egcode) again and then the method call\n",
"you will generate (modulo randomness) the plot"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8cb35166",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_additive(amf, T)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "c41e38c8",
"metadata": {},
"source": [
"When we plot multiple realizations of a component in the 2nd, 3rd, and 4th panels, we also plot the population 95% probability coverage sets computed using the LinearStateSpace class.\n",
"\n",
"We have chosen to simulate many paths, all starting from the *same* non-random initial conditions $ x_0, y_0 $ (you can tell this from the shape of the 95% probability coverage shaded areas).\n",
"\n",
"Notice tell-tale signs of these probability coverage shaded areas\n",
"\n",
"- the purple one for the martingale component $ m_t $ grows with\n",
" $ \\sqrt{t} $ \n",
"- the green one for the stationary component $ s_t $ converges to a\n",
" constant band "
]
},
{
"cell_type": "markdown",
"id": "541dab99",
"metadata": {},
"source": [
"### Associated multiplicative functional\n",
"\n",
"Where $ \\{y_t\\} $ is our additive functional, let $ M_t = \\exp(y_t) $.\n",
"\n",
"As mentioned above, the process $ \\{M_t\\} $ is called a **multiplicative functional**.\n",
"\n",
"Corresponding to the additive decomposition described above we have a multiplicative decomposition of $ M_t $\n",
"\n",
"$$\n",
"\\frac{M_t}{M_0}\n",
"= \\exp (t \\nu) \\exp \\Bigl(\\sum_{j=1}^t H \\cdot Z_j \\Bigr) \\exp \\biggl( D(I-A)^{-1} x_0 - D(I-A)^{-1} x_t \\biggr)\n",
"$$\n",
"\n",
"or\n",
"\n",
"$$\n",
"\\frac{M_t}{M_0} = \\exp\\left( \\tilde \\nu t \\right) \\Biggl( \\frac{\\widetilde M_t}{\\widetilde M_0}\\Biggr) \\left( \\frac{\\tilde e (X_0)} {\\tilde e(x_t)} \\right)\n",
"$$\n",
"\n",
"where\n",
"\n",
"$$\n",
"\\tilde \\nu = \\nu + \\frac{H \\cdot H}{2} ,\n",
"\\quad\n",
"\\widetilde M_t = \\exp \\biggl( \\sum_{j=1}^t \\biggl(H \\cdot z_j -\\frac{ H \\cdot H }{2} \\biggr) \\biggr), \\quad \\widetilde M_0 =1\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\tilde e(x) = \\exp[g(x)] = \\exp \\bigl[ D (I - A)^{-1} x \\bigr]\n",
"$$\n",
"\n",
"An instance of class `AMF_LSS_VAR` ([above](#amf-lss)) includes this associated multiplicative functional as an attribute.\n",
"\n",
"Let’s plot this multiplicative functional for our example.\n",
"\n",
"If you run [the code that first simulated that example](#addfunc-egcode) again and then the method call in the cell below you’ll\n",
"obtain the graph in the next cell."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7b7d4dd8",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"plot_multiplicative(amf, T)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "3433aed1",
"metadata": {},
"source": [
"As before, when we plotted multiple realizations of a component in the 2nd, 3rd, and 4th panels, we also plotted population 95% confidence bands computed using the LinearStateSpace class.\n",
"\n",
"Comparing this figure and the last also helps show how geometric growth differs from\n",
"arithmetic growth.\n",
"\n",
"The top right panel of the above graph shows a panel of martingales associated with the panel of $ M_t = \\exp(y_t) $ that we have generated\n",
"for a limited horizon $ T $.\n",
"\n",
"It is interesting to how the martingale behaves as $ T \\rightarrow +\\infty $.\n",
"\n",
"Let’s see what happens when we set $ T = 12000 $ instead of $ 150 $."
]
},
{
"cell_type": "markdown",
"id": "8d5c0195",
"metadata": {},
"source": [
"### Peculiar large sample property\n",
"\n",
"Hansen and Sargent [[Hansen and Sargent, 2024](https://python-advanced.quantecon.org/zreferences.html#id247)] (ch. 8) describe the following two properties of the martingale component\n",
"$ \\widetilde M_t $ of the multiplicative decomposition\n",
"\n",
"- while $ E_0 \\widetilde M_t = 1 $ for all $ t \\geq 0 $,\n",
" nevertheless $ \\ldots $ \n",
"- as $ t \\rightarrow +\\infty $, $ \\widetilde M_t $ converges to\n",
" zero almost surely \n",
"\n",
"\n",
"The first property follows from the fact that $ \\widetilde M_t $ is a multiplicative martingale with initial condition\n",
"$ \\widetilde M_0 = 1 $.\n",
"\n",
"The second is a **peculiar property** noted and proved by Hansen and Sargent [[Hansen and Sargent, 2024](https://python-advanced.quantecon.org/zreferences.html#id247)].\n",
"\n",
"The following simulation of many paths of $ \\widetilde M_t $ illustrates both properties"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "260ec23e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"np.random.seed(10021987)\n",
"plot_martingales(amf, 12000)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "a786073b",
"metadata": {},
"source": [
"The dotted line in the above graph is the mean $ E \\tilde M_t = 1 $ of the martingale.\n",
"\n",
"It remains constant at unity, illustrating the first property.\n",
"\n",
"The purple 95 percent frequency coverage interval collapses around zero, illustrating the second property."
]
},
{
"cell_type": "markdown",
"id": "633675a2",
"metadata": {},
"source": [
"## More about the multiplicative martingale\n",
"\n",
"Let’s drill down and study probability distribution of the multiplicative martingale $ \\{\\widetilde M_t\\}_{t=0}^\\infty $ in\n",
"more detail.\n",
"\n",
"As we have seen, it has representation\n",
"\n",
"$$\n",
"\\widetilde M_t = \\exp \\biggl( \\sum_{j=1}^t \\biggl(H \\cdot z_j -\\frac{ H \\cdot H }{2} \\biggr) \\biggr), \\quad \\widetilde M_0 =1\n",
"$$\n",
"\n",
"where $ H = [F + D(I-A)^{-1} B] $.\n",
"\n",
"It follows that $ \\log {\\widetilde M}_t \\sim {\\mathcal N} ( -\\frac{t H \\cdot H}{2}, t H \\cdot H ) $ and that consequently $ {\\widetilde M}_t $ is log normal.\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"id": "40e831d8",
"metadata": {},
"source": [
"### Moments, skewness, and kurtosis\n",
"\n",
"Because $ \\log \\widetilde M_t $ is Gaussian, the moments of $ \\widetilde M_t $ have simple closed forms.\n",
"\n",
"Let\n",
"\n",
"$$\n",
"a_t := t (H \\cdot H) = \\mathrm{Var}(\\log \\widetilde M_t).\n",
"$$\n",
"\n",
"Then for any $ k \\geq 1 $,\n",
"\n",
"$$\n",
"E\\left[\\widetilde M_t^k\\right] = \\exp\\left(\\frac{1}{2} k(k-1) a_t\\right).\n",
"$$\n",
"\n",
"In particular, $ E[\\widetilde M_t]=1 $ for all $ t $, while higher-order moments grow rapidly in $ t $ whenever $ H \\neq 0 $.\n",
"\n",
"Standard summary statistics can be written in terms of $ a_t $ as\n",
"\n",
"$$\n",
"\\mathrm{Var}(\\widetilde M_t) = e^{a_t} - 1,\n",
"$$\n",
"\n",
"$$\n",
"\\mathrm{Skew}(\\widetilde M_t) = (e^{a_t} + 2)\\sqrt{e^{a_t}-1},\n",
"$$\n",
"\n",
"$$\n",
"\\mathrm{Kurt}(\\widetilde M_t) = e^{4a_t} + 2e^{3a_t} + 3e^{2a_t} - 3.\n",
"$$\n",
"\n",
"(Here, [excess kurtosis](https://en.wikipedia.org/wiki/Kurtosis) is obtained by subtracting 3 from kurtosis)\n",
"\n",
"These formulas show that everything depends on $ a_t = t(H \\cdot H) $.\n",
"\n",
"To illustrate these formulas, let’s plot the growth of a few moments and summary statistics over time.\n",
"\n",
"The following functions compute $ \\log_{10} E[\\widetilde M_t^k] $ and $ \\log_{10} \\mathrm{Var}(\\widetilde M_t) $, $ \\log_{10} \\mathrm{Skew}(\\widetilde M_t) $, and $ \\log_{10} \\mathrm{Kurt}(\\widetilde M_t) $"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2098e30c",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def squared_norm(H):\n",
" H = np.atleast_1d(H).astype(float)\n",
" return H @ H\n",
"\n",
"def log10_raw_moment(H, t, k):\n",
" a = t * squared_norm(H)\n",
" moment = np.exp(0.5 * k * (k - 1) * a)\n",
" return np.log10(moment)\n",
"\n",
"def log10_var_skew_kurt(H, t):\n",
" a = t * squared_norm(H)\n",
"\n",
" var = np.exp(a) - 1\n",
" skew = (np.exp(a) + 2) * np.sqrt(var)\n",
" kurt = np.exp(4*a) + 2*np.exp(3*a) + 3*np.exp(2*a) - 3\n",
"\n",
" return np.log10(var), np.log10(skew), np.log10(kurt)"
]
},
{
"cell_type": "markdown",
"id": "729fe493",
"metadata": {},
"source": [
"Let’s illustrate the growth of $ k $th-order raw moments on a log scale for $ H = 2.0 $"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e1820026",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"H_example = 2.0\n",
"t_grid = np.arange(1, 9)\n",
"ks = [1, 2, 3, 4]\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 5))\n",
"for k in ks:\n",
" ax.plot(t_grid,\n",
" [log10_raw_moment(H_example, t, k) for t in t_grid],\n",
" marker=\"o\",\n",
" label=f\"$k={k}$\")\n",
"\n",
"ax.set_xlabel(\"$t$\")\n",
"ax.set_ylabel(r\"$\\log_{10} E[\\widetilde M_t^k]$\")\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "bc25a662",
"metadata": {},
"source": [
"Next we plot the growth of variance, skewness, and kurtosis"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "58dd816b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"log10_var = []\n",
"log10_skew = []\n",
"log10_kurt = []\n",
"\n",
"for t in t_grid:\n",
" lv, ls, lk = log10_var_skew_kurt(H_example, t)\n",
" log10_var.append(lv)\n",
" log10_skew.append(ls)\n",
" log10_kurt.append(lk)\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 5))\n",
"ax.plot(t_grid, log10_var, marker=\"o\", \n",
" label=r\"$\\log_{10}\\,\\mathrm{Var}$\")\n",
"ax.plot(t_grid, log10_skew, marker=\"o\", \n",
" label=r\"$\\log_{10}\\,\\mathrm{Skew}$\")\n",
"ax.plot(t_grid, log10_kurt, marker=\"o\", \n",
" label=r\"$\\log_{10}\\,\\mathrm{Kurt}$\")\n",
"ax.set_xlabel(\"$t$\")\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "55595ae0",
"metadata": {},
"source": [
"These plots help explain the *peculiar property* of the multiplicative martingale.\n",
"\n",
"The rapid growth of variance, skewness, and kurtosis reveals that the distribution of $ \\widetilde M_t $ becomes increasingly right-skewed over time while $ E[\\widetilde M_t] = 1 $ for all $ t $.\n",
"\n",
"This means that most probability density concentrates near zero, while a long right tail preserves the unit mean."
]
},
{
"cell_type": "markdown",
"id": "3576b622",
"metadata": {},
"source": [
"### Simulating a multiplicative martingale again\n",
"\n",
"Next, we want a program to simulate the likelihood ratio process $ \\{ \\tilde{M}_t \\}_{t=0}^\\infty $.\n",
"\n",
"In particular, we want to simulate 5000 sample paths of length $ T $ for the case in which $ x $ is a scalar and\n",
"$ [A, B, D, F] = [0.8, 0.001, 1.0, 0.01] $ and $ \\nu = 0.005 $.\n",
"\n",
"After accomplishing this, we want to display and study histograms of $ \\tilde{M}_T^i $ for various values of $ T $.\n",
"\n",
"Let’s first write a program to simulate sample paths of $ \\{ x_t, y_{t} \\}_{t=0}^{\\infty} $.\n",
"\n",
"We’ll do this by formulating the additive functional as a linear state space model and putting the [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to work via our `AMF_LSS_VAR` class [defined above](#amf-lss).\n",
"\n",
"The following code adds some simple functions that make it straightforward to generate sample paths from an instance of `AMF_LSS_VAR`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "78a9fa8e",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def simulate_xy(amf, T):\n",
" \"Simulate individual paths.\"\n",
" foo, bar = amf.lss.simulate(T)\n",
" x = bar[0, :]\n",
" y = bar[1, :]\n",
"\n",
" return x, y\n",
"\n",
"def simulate_paths(amf, T=150, I=5000):\n",
" \"Simulate multiple independent paths.\"\n",
"\n",
" # Allocate space\n",
" storeX = np.empty((I, T))\n",
" storeY = np.empty((I, T))\n",
"\n",
" for i in range(I):\n",
" # Do specific simulation\n",
" x, y = simulate_xy(amf, T)\n",
"\n",
" # Fill in our storage matrices\n",
" storeX[i, :] = x\n",
" storeY[i, :] = y\n",
"\n",
" return storeX, storeY\n",
"\n",
"def population_means(amf, T=150):\n",
" # Allocate Space\n",
" xmean = np.empty(T)\n",
" ymean = np.empty(T)\n",
"\n",
" # Pull out moment generator\n",
" moment_generator = amf.lss.moment_sequence()\n",
"\n",
" for tt in range (T):\n",
" tmoms = next(moment_generator)\n",
" ymeans = tmoms[1]\n",
" xmean[tt] = ymeans[0]\n",
" ymean[tt] = ymeans[1]\n",
"\n",
" return xmean, ymean"
]
},
{
"cell_type": "markdown",
"id": "fbdb222e",
"metadata": {},
"source": [
"Now that we have these functions in our toolkit, let’s apply them to run some\n",
"simulations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "718fe664",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def simulate_martingale_components(amf, T=1000, I=5000):\n",
" # Get the multiplicative decomposition\n",
" ν, H, g = amf.multiplicative_decomp()\n",
"\n",
" # Allocate space\n",
" add_mart_comp = np.empty((I, T))\n",
"\n",
" # Simulate and pull out additive martingale component\n",
" for i in range(I):\n",
" foo, bar = amf.lss.simulate(T)\n",
"\n",
" # Martingale component is third component\n",
" add_mart_comp[i, :] = bar[2, :]\n",
"\n",
" mul_mart_comp = np.exp(add_mart_comp - (np.arange(T) * H**2)/2)\n",
"\n",
" return add_mart_comp, mul_mart_comp\n",
"\n",
"\n",
"# Build model\n",
"amf_2 = AMF_LSS_VAR(0.8, 0.001, 1.0, 0.01,.005)\n",
"\n",
"amc, mmc = simulate_martingale_components(amf_2, 1000, 5000)\n",
"\n",
"amcT = amc[:, -1]\n",
"mmcT = mmc[:, -1]\n",
"\n",
"print(\"The (min, mean, max) of additive Martingale component in period T is\")\n",
"print(f\"\\t ({np.min(amcT)}, {np.mean(amcT)}, {np.max(amcT)})\")\n",
"\n",
"print(\"The (min, mean, max) of multiplicative Martingale component \\\n",
"in period T is\")\n",
"print(f\"\\t ({np.min(mmcT)}, {np.mean(mmcT)}, {np.max(mmcT)})\")"
]
},
{
"cell_type": "markdown",
"id": "f8c9030c",
"metadata": {},
"source": [
"Let’s plot the probability density functions for $ \\log {\\widetilde M}_t $ for\n",
"$ t=100, 500, 1000, 10000, 100000 $.\n",
"\n",
"Then let’s use the plots to investigate how these densities evolve through time.\n",
"\n",
"We will plot the densities of $ \\log {\\widetilde M}_t $ for different values of $ t $.\n",
"\n",
">**Note**\n",
">\n",
">`scipy.stats.lognorm` expects you to pass the standard deviation\n",
"first $ (tH \\cdot H) $ and then the exponent of the mean as a\n",
"keyword argument `scale` (`scale=np.exp(-t * H2 / 2)`).\n",
"\n",
"- See the documentation [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html#scipy.stats.lognorm). \n",
"\n",
"\n",
"This is peculiar, so make sure you are careful in working with the log normal distribution.\n",
"\n",
"Here is some code that tackles these tasks"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc422261",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"def Mtilde_t_density(amf, t, xmin=1e-8, xmax=5.0, npts=5000):\n",
"\n",
" # Pull out the multiplicative decomposition\n",
" νtilde, H, g = amf.multiplicative_decomp()\n",
" H2 = (H @ H.T).item() \n",
"\n",
" # The distribution\n",
" mdist = lognorm(np.sqrt(t*H2), scale=np.exp(-t*H2/2))\n",
" x = np.linspace(xmin, xmax, npts)\n",
" pdf = mdist.pdf(x)\n",
"\n",
" return x, pdf\n",
"\n",
"\n",
"def logMtilde_t_density(amf, t, xmin=-15.0, xmax=15.0, npts=5000):\n",
"\n",
" # Pull out the multiplicative decomposition\n",
" νtilde, H, g = amf.multiplicative_decomp()\n",
" H2 = (H @ H.T).item() \n",
"\n",
" # The distribution\n",
" lmdist = norm(-t*H2/2, np.sqrt(t*H2))\n",
" x = np.linspace(xmin, xmax, npts)\n",
" pdf = lmdist.pdf(x)\n",
"\n",
" return x, pdf\n",
"\n",
"\n",
"times_to_plot = [10, 100, 500, 1000, 2500, 5000]\n",
"dens_to_plot = map(lambda t: Mtilde_t_density(amf_2, t, xmin=1e-8, xmax=6.0),\n",
" times_to_plot)\n",
"ldens_to_plot = map(lambda t: logMtilde_t_density(amf_2, t, xmin=-10.0,\n",
" xmax=10.0), times_to_plot)\n",
"\n",
"fig, ax = plt.subplots(3, 2, figsize=(14, 14))\n",
"ax = ax.flatten()\n",
"\n",
"fig.suptitle(r\"Densities of $\\tilde{M}_t$\", fontsize=18, y=1.02)\n",
"for (it, dens_t) in enumerate(dens_to_plot):\n",
" x, pdf = dens_t\n",
" ax[it].set_title(f\"Density for time {times_to_plot[it]}\")\n",
" ax[it].fill_between(x, np.zeros_like(pdf), pdf)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "d12f0f18",
"metadata": {},
"source": [
"These probability density functions again help us understand mechanics underlying the **peculiar property** of our multiplicative martingale\n",
"\n",
"- As $ T $ grows, most of the probability mass shifts leftward toward zero. \n",
"- For example, note that most mass is near $ 1 $ for $ T =10 $ or $ T = 100 $ but\n",
" most of it is near $ 0 $ for $ T = 5000 $. \n",
"- As $ T $ grows, the tail of the density of $ \\widetilde M_T $ lengthens toward the right. \n",
"- Enough mass moves toward the right tail to keep $ E \\widetilde M_T = 1 $\n",
" even as most mass in the distribution of $ \\widetilde M_T $ collapses around $ 0 $. "
]
},
{
"cell_type": "markdown",
"id": "8dbe43a2",
"metadata": {},
"source": [
"### Multiplicative martingale as likelihood ratio process\n",
"\n",
"[This lecture](https://python.quantecon.org/likelihood_ratio_process.html) studies **likelihood processes**\n",
"and **likelihood ratio processes**.\n",
"\n",
"A **likelihood ratio process** is a multiplicative martingale with mean unity.\n",
"\n",
"Likelihood ratio processes exhibit the peculiar property that naturally also appears\n",
"[here](https://python.quantecon.org/likelihood_ratio_process.html)."
]
},
{
"cell_type": "markdown",
"id": "cc49ef23",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"id": "00e258a2",
"metadata": {},
"source": [
"## Exercise 32.1\n",
"\n",
"Fix $ H = 0.5 $ and consider the log-normal martingale $ \\widetilde M_t $.\n",
"\n",
"1. Using the formula for $ E[\\widetilde M_t^k] $, plot $ \\log_{10} E[\\widetilde M_t^k] $ against $ t $ for $ k = 1,2,3,4 $ and $ t = 1,\\ldots,8 $. \n",
"1. Verify from your plot (and the formula) that $ \\log_{10} E[\\widetilde M_t^k] $ is linear in $ t $, and compute its slope as a function of $ k $ and $ H $. "
]
},
{
"cell_type": "markdown",
"id": "400a1069",
"metadata": {},
"source": [
"## Solution\n",
"\n",
"Here is one solution using the `log10_raw_moment` function"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6d2e77b",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"H = 0.5\n",
"ks = [1, 2, 3, 4]\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 5))\n",
"for k in ks:\n",
" ax.plot(t_grid,\n",
" [log10_raw_moment(H, t, k) for t in t_grid],\n",
" marker=\"o\", label=f\"$k={k}$\")\n",
"ax.legend()\n",
"ax.set_xlabel(\"$t$\")\n",
"ax.set_ylabel(r\"$\\log_{10} E[\\widetilde M_t^k]$\")\n",
"plt.show()\n",
"\n",
"H2 = H**2\n",
"for k in ks:\n",
" slope = np.log10(np.exp(0.5 * k * (k - 1) * H2))\n",
" print(f\"k={k}: slope = {slope:.4f}\")"
]
},
{
"cell_type": "markdown",
"id": "205a67fe",
"metadata": {},
"source": [
"When $ k=1 $, the slope is zero, reflecting $ E[\\widetilde M_t]=1 $ for all $ t $."
]
},
{
"cell_type": "markdown",
"id": "16f88750",
"metadata": {},
"source": [
"## Exercise 32.2\n",
"\n",
"For the same martingale $ \\widetilde M_t $ (take $ H=0.5 $ as in the previous exercise),\n",
"\n",
"1. Plot $ \\log_{10}\\mathrm{Var}(\\widetilde M_t) $, $ \\log_{10}\\mathrm{Skew}(\\widetilde M_t) $, and $ \\log_{10}\\mathrm{Kurt}(\\widetilde M_t) $ for $ t = 1,\\ldots,8 $. \n",
"1. Plot $ \\log_{10}(\\mathrm{Var}(\\widetilde M_t)+1) $ and explain why it is exactly linear in $ t $. "
]
},
{
"cell_type": "markdown",
"id": "657d1d42",
"metadata": {},
"source": [
"## Solution\n",
"\n",
"Here is one solution using the `log10_var_skew_kurt` function"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1b58c703",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"H = 0.5\n",
"\n",
"log10_var = []\n",
"log10_skew = []\n",
"log10_kurt = []\n",
"log10_var_plus_1 = []\n",
"\n",
"for t in t_grid:\n",
" lv, ls, lk = log10_var_skew_kurt(H, t)\n",
" log10_var.append(lv)\n",
" log10_skew.append(ls)\n",
" log10_kurt.append(lk)\n",
" log10_var_plus_1.append(np.log10(np.exp(t * H**2)))\n",
"\n",
"fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
"\n",
"ax[0].plot(t_grid, log10_var, marker=\"o\",\n",
" label=r\"$\\log_{10}\\,\\mathrm{Var}$\")\n",
"ax[0].plot(t_grid, log10_skew, marker=\"o\",\n",
" label=r\"$\\log_{10}\\,\\mathrm{Skew}$\")\n",
"ax[0].plot(t_grid, log10_kurt, marker=\"o\",\n",
" label=r\"$\\log_{10}\\,\\mathrm{Kurt}$\")\n",
"ax[0].set_xlabel(\"$t$\")\n",
"ax[0].legend()\n",
"\n",
"ax[1].plot(t_grid, log10_var_plus_1, marker=\"o\",\n",
" label=r\"$\\log_{10}(\\mathrm{Var}+1)$\")\n",
"ax[1].set_xlabel(\"$t$\")\n",
"ax[1].legend()\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "8f639ff2",
"metadata": {},
"source": [
"Since $ \\mathrm{Var}(\\widetilde M_t) + 1 = e^{a_t} $ and $ a_t = t(H \\cdot H) $, we have\n",
"$ \\log_{10}(\\mathrm{Var}(\\widetilde M_t)+1) = a_t \\log_{10}(e) $, which is linear in $ t $."
]
},
{
"cell_type": "markdown",
"id": "210cfd05",
"metadata": {},
"source": [
"## Exercise 32.3\n",
"\n",
"The formulas in [Moments, skewness, and kurtosis](#mult-mart-moment) depend on $ H $ only through $ H \\cdot H $.\n",
"\n",
"1. Choose two different vectors $ H_1 $ and $ H_2 $ with the same value of $ H \\cdot H $ (for example, $ H_1 = 1 $ and $ H_2=(0.6,0.8) $). \n",
"\n",
"\n",
"Show that $ \\log_{10} E[\\widetilde M_t^k] $ agrees for all $ t $ and $ k $.\n",
"\n",
"1. Compare your results to a larger value, such as $ H_3=(1,2,1.5) $, and comment on how quickly the higher-order moments grow. "
]
},
{
"cell_type": "markdown",
"id": "dd6a93fb",
"metadata": {},
"source": [
"## Solution\n",
"\n",
"Here is one solution using the `squared_norm` and `log10_raw_moment` functions"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "02ef34da",
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"H1 = 1.0\n",
"H2 = np.array([0.6, 0.8])\n",
"H3 = np.array([1.0, 2.0, 1.5])\n",
"\n",
"print(\"H1.H1 =\", squared_norm(H1))\n",
"print(\"H2.H2 =\", squared_norm(H2))\n",
"print(\"H3.H3 =\", squared_norm(H3))\n",
"\n",
"k = 4\n",
"\n",
"fig, ax = plt.subplots(figsize=(9, 5))\n",
"ax.plot(t_grid, [log10_raw_moment(H1, t, k) for t in t_grid],\n",
" marker=\"o\", label=r\"$H_1=1$\",\n",
" alpha=0.7, markersize=9)\n",
"ax.plot(t_grid, [log10_raw_moment(H2, t, k) for t in t_grid],\n",
" marker=\"v\", label=r\"$H_2=(0.6,0.8)$\",\n",
" alpha=0.7)\n",
"ax.plot(t_grid, [log10_raw_moment(H3, t, k) for t in t_grid],\n",
" marker=\"o\", label=r\"$H_3=(1,2,1.5)$\",\n",
" alpha=0.7)\n",
"ax.set_xlabel(\"$t$\")\n",
"ax.set_ylabel(rf\"$\\log_{{10}} E[\\widetilde M_t^{{{k}}}]$\")\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "f1726967",
"metadata": {},
"source": [
"The curves for $ H_1 $ and $ H_2 $ coincide because they have the same squared norm $ H \\cdot H $.\n",
"\n",
"This exercise shows that only the squared norm $ H \\cdot H $ and $ t $ matters for the moments of $ \\widetilde M_t $."
]
}
],
"metadata": {
"date": 1769597367.3400552,
"filename": "additive_functionals.md",
"kernelspec": {
"display_name": "Python",
"language": "python3",
"name": "python3"
},
"title": "Additive and Multiplicative Functionals"
},
"nbformat": 4,
"nbformat_minor": 5
}