{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Poisson-Zero Hidden Markov Model with Time-Varying Transitions\n", "\n", "## Summary\n", "\n", "The following exposition uses [`pymc3-hmm`](https://github.com/AmpersandTV/pymc3-hmm) to simulate and estimate a hidden Markov model (HMM) with a time-varying transition matrix and emissions consisting of a Poisson random variable and the point at zero." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "\n", "import patsy\n", "\n", "import theano.tensor as tt\n", "\n", "import pymc3 as pm\n", "\n", "from pymc3_hmm.utils import multilogit_inv\n", "from pymc3_hmm.distributions import SwitchingProcess, DiscreteMarkovChain\n", "from pymc3_hmm.step_methods import FFBSStep" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "Our observation model can be described as follows:\n", "\n", "\\begin{align}\n", " \\label{eq:pois-zero-model}\n", " Y_t &= S_t \\epsilon_t,\\quad\n", " \\epsilon_t \\sim \\operatorname{Pois}\\left( \\mu \\right)\n", " \\\\\n", " S_t &\\sim \\operatorname{Bern}\\left( \\pi_t \\right)\n", " \\;,\n", "\\end{align}\n", "where $y_t \\sim Y_t$ are the observed values sampled from the observation distribution, $Y_t$, spanning $t \\in \\left\\{0, \\dots, T \\right\\}$.\n", "\n", "The \"hidden\" state sequence, $\\{S_t\\}$, is driven by the following Markov relationship:\n", "\n", "\\begin{equation*}\n", " \\operatorname{P}\\left(S_t \\mid S_{t-1}\\right) \\triangleq \\Gamma_{t,t-1} \\in \\mathbb{R}^{2 \\times 2}_{[0, 1]}\n", "\\end{equation*}\n", "\n", "The marginal state probability, $\\pi_t$, is then given by\n", "\n", "\\begin{equation*}\n", " \\begin{aligned}\n", " \\operatorname{P}\\left( S_t \\right)\n", " &= \\int_{S_{t-1}} \\operatorname{P}\\left(S_t \\mid S_{t-1}\\right)\n", " \\operatorname{dP}\\left(S_{t-1}\\right)\n", " \\\\\n", " &=\n", " \\begin{pmatrix}\n", " \\Gamma^{(0, 0)}_{t,t-1} & \\Gamma^{(0, 1)}_{t,t-1}\n", " \\\\\n", " \\Gamma^{(1, 0)}_{t,t-1} & \\Gamma^{(1, 1)}_{t,t-1}\n", " \\end{pmatrix}^\\top\n", " \\begin{pmatrix}\n", " \\pi_{t-1}\n", " \\\\\n", " 1 - \\pi_{t-1}\n", " \\end{pmatrix}\n", " \\\\\n", " &=\n", " \\begin{pmatrix}\n", " \\pi_{t}\n", " \\\\\n", " 1 - \\pi_{t}\n", " \\end{pmatrix}\n", " \\;.\n", " \\end{aligned}\n", "\\end{equation*}\n", "\n", "In this example, the rows of our transition matrix, $\\Gamma^{(r)}_{t, t-1}$ for $r \\in \\{0, 1\\}$, are driven by a logistic regression:\n", "\n", "\\begin{gather*}\n", " \\Gamma^{(r)}_{t, t-1} = \\operatorname{logit^{-1}}\\left( X_t \\xi_r \\right)\n", "\\end{gather*}\n", "\n", "where $X_t \\in \\mathbb{R}^{T \\times N}$ is a covariate matrix and $\\xi_r \\in \\mathbb{R}^N$ is the regression parameter vector for row $r$.\n", "\n", "In the remainder of this exposition, we will assume normal priors for each $\\xi_r$, the conjugate Gamma prior for $\\mu$, and a Dirichlet prior for $\\pi_0$.\n", "\n", "## Simulation\n", "\n", "For these simulations, we will generate a time series and make the $\\xi_r$ regression consist of fixed-effects for a seasonal component based on weekdays. In other words, the transition probabilities will vary based on the day of week." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def create_poisson_zero_hmm_tv(mu, xis, pi_0, observed):\n", " \n", " z_tt = tt.tensordot(X_df.values, xis, axes=((1,), (0,)))\n", "\n", " P_tt = multilogit_inv(z_tt)\n", "\n", " P_rv = pm.Deterministic(\"P_t\", P_tt)\n", "\n", " S_rv = DiscreteMarkovChain(\"S_t\", P_rv, pi_0, shape=np.shape(observed)[-1])\n", "\n", " Y_rv = SwitchingProcess(\n", " \"Y_t\", [pm.Constant.dist(0), pm.Poisson.dist(mu)], S_rv, observed=observed\n", " )\n", " return Y_rv\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "np.random.seed(2032)\n", "\n", "start_date = pd.Timestamp(\"2020-01-01 00:00:00\")\n", "time_index = pd.date_range(\n", " start=start_date, end=start_date + pd.Timedelta(\"30W\"), closed=\"left\", freq=\"1h\"\n", ")\n", "X_ind_df = pd.DataFrame(\n", " {\n", " \"weekday\": time_index.weekday,\n", " },\n", " index=time_index,\n", ")\n", "\n", "formula_str = \"~ 1 + C(weekday)\"\n", "X_df = patsy.dmatrix(formula_str, X_ind_df, return_type=\"dataframe\")\n", "\n", "xi_0_true = pd.Series(\n", " # The coefficients used to compute the state zero-to-zero transition probabilities\n", " np.array([2.0, -5.0, -3.0, 0.0, 0.0, -5.0, -5.0]),\n", " index=X_df.columns,\n", ")\n", "\n", "xi_1_true = pd.Series(\n", " # The coefficients for the state one-to-zero transition probabilities\n", " np.array([-2.0, -1.0, 3.0, 0.0, 0.0, 5.0, 5.0]),\n", " index=X_df.columns,\n", ")\n", "\n", "xis_true = tt.as_tensor(np.stack([xi_0_true, xi_1_true], axis=1)[..., None], name=\"xis\")\n", "\n", "mu_true = 50\n", "p_0_true = tt.as_tensor(np.r_[0.0, 1.0])\n", "\n", "with pm.Model(theano_config={\"compute_test_value\": \"ignore\"}) as sim_model:\n", " _ = create_poisson_zero_hmm_tv(mu_true, xis_true, p_0_true, np.zeros(X_df.shape[0]))\n", "\n", "sim_point = pm.sample_prior_predictive(samples=1, model=sim_model)\n", "sim_point[\"Y_t\"] = sim_point[\"Y_t\"].squeeze()\n", "\n", "y_t = sim_point[\"Y_t\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimation\n", "\n", "We will use the \"true\" data-generating observation model to estimate the parameters $\\mu$ and $\\Gamma_{t, t-1}$ (the latter as rows denoted by `p_0` and `p_1`). For demonstration purposes, we choose hyper-parameters for the $\\mu$ prior that are \"far\" from the true $\\mu$ value.\n", "\n", "The sampling steps for $S_t$ are performed using forward-filtering backward-sampling (FFBS), while sampling for $\\mu$, $\\pi_0$, and the $\\xi_r$ use NUTS.\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Sequential sampling (1 chains in 1 job)\n", "CompoundStep\n", ">FFBSStep: [S_t]\n", ">NUTS: [p_0, xis, mu]\n" ] }, { "data": { "text/html": [ "\n", "