rhydomako/chow_test.ipynb

## chow_test.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chow test for structural break evaluation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import norm, f\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = pd.date_range(start='2019-01', end='2019-11', freq='M')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_change():\n",
    "    Y = []\n",
    "    for i, date in enumerate(idx):\n",
    "        if date >= pd.to_datetime('2019-08-01'):\n",
    "            Y.append( norm(12000 - 300*(i-2), 200).rvs() )\n",
    "        else:\n",
    "            Y.append( norm(12000, 200).rvs() )\n",
    "\n",
    "    df = pd.DataFrame(Y, index=idx, columns=['Y'])\n",
    "    df['b'] = 1\n",
    "    df['X'] = [m.month for m in s.index]\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Y</th>\n",
       "      <th>b</th>\n",
       "      <th>X</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2019-01-31</th>\n",
       "      <td>12037.499493</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-02-28</th>\n",
       "      <td>11947.368249</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-03-31</th>\n",
       "      <td>12158.814023</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-04-30</th>\n",
       "      <td>11876.049085</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-05-31</th>\n",
       "      <td>12082.947888</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-06-30</th>\n",
       "      <td>11794.185150</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-07-31</th>\n",
       "      <td>12337.295125</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-08-31</th>\n",
       "      <td>10540.547065</td>\n",
       "      <td>1</td>\n",
       "      <td>8</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-09-30</th>\n",
       "      <td>10198.610762</td>\n",
       "      <td>1</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2019-10-31</th>\n",
       "      <td>9959.139408</td>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                       Y  b   X\n",
       "2019-01-31  12037.499493  1   1\n",
       "2019-02-28  11947.368249  1   2\n",
       "2019-03-31  12158.814023  1   3\n",
       "2019-04-30  11876.049085  1   4\n",
       "2019-05-31  12082.947888  1   5\n",
       "2019-06-30  11794.185150  1   6\n",
       "2019-07-31  12337.295125  1   7\n",
       "2019-08-31  10540.547065  1   8\n",
       "2019-09-30  10198.610762  1   9\n",
       "2019-10-31   9959.139408  1  10"
      ]
     },
     "execution_count": 264,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = generate_change()\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 15000)"
      ]
     },
     "execution_count": 265,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEHCAYAAAC9TnFRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8XOV97/HP19psSbbkRcbGS2yC\ngbCETWXJQhYStnILtElLkluchFuXNi1Lm9skTV8lJfttWgJpQmKWBvpKWQJpoQRCXAIhNIEgg83m\nBDs42MIGC8uWV8mW/Lt/nEfWSJbt45HGGsP3/XrpNXOe85zRb8aj+c5zznOOFRGYmZntzaiRLsDM\nzA4MDgwzM8vFgWFmZrk4MMzMLBcHhpmZ5eLAMDOzXBwYZmaWiwPDzMxycWCYmVkulSNdQLEmTZoU\ns2bNGukyzMwOKAsXLnwtIpqK2faADYxZs2bR0tIy0mWYmR1QJL1U7LbeJWVmZrk4MMzMLBcHhpmZ\n5eLAMDOzXBwYZmaWiwPDzMxycWCYmVkuDgwzM8vFgWFmZrnsNTAk3SRpjaRnB1n3SUkhaVJalqRr\nJS2T9LSkEwr6zpW0NP3MLWg/UdIzaZtrJWm4npyZmQ2fPCOM7wJnDWyUNAN4P7CioPlsYE76mQdc\nl/pOAK4ETgZOAq6UND5tc13q27vdLr/LzMxG3l4DIyIeAdoHWXU18DdAFLSdB9wSmceARklTgTOB\nBRHRHhHrgAXAWWnduIj4RUQEcAtw/tCekpmZlUJRxzAk/R7wckQsHrBqGrCyYLk1te2pvXWQ9t39\n3nmSWiS1tLW1FVO6mZkVaZ8DQ1It8Fng7wdbPUhbFNE+qIiYHxHNEdHc1FTU1XnNzKxIxYww3gzM\nBhZL+i0wHXhS0hSyEcKMgr7TgVV7aZ8+SLuZmZWZfQ6MiHgmIiZHxKyImEX2oX9CRLwC3ANclGZL\nnQJ0RMRq4AHgDEnj08HuM4AH0rqNkk5Js6MuAu4epudmZmbDKM+02luBXwCHS2qVdPEeut8HvAgs\nA64H/hwgItqBzwNPpJ+rUhvAnwE3pG1+A9xf3FMxM7NSUjY56cDT3Nwc/h/3zMz2jaSFEdFczLY+\n09vMzHJxYJiZWS4ODDMzy8WBYWZmuTgwzMwsFweGmZnl4sAwM7NcHBhmZpaLA8PMzHJxYJiZWS4O\nDDMzy8WBYWZmuTgwzMwsFweGmZnl4sAwM7NcHBhmZpaLA8PMzHJxYJjZXkUEq9Zv5YHnXmH5a5tH\nuhwbIZUjXYCZlZ+Ordt5prWDRSvXsWhlB4tb19O2sQuA6spRfOqsI/jY22YxapRGuFLbnxwYZm9w\nXd09LFm9kcUr17N45XoWta7nxba+UcQhTXW889BJHDujkSOmjOX6ny3n8/c+z09+9Spf++CxTG0Y\nM4LV2/6018CQdBNwLrAmIo5Obf8I/C9gG/Ab4GMRsT6t+wxwMdADXBoRD6T2s4BrgArghoj4Smqf\nDdwGTACeBP44IrYN55N8PdrWvYN1W7axdtM22jdvY+3mLto3997fRvumrC0CpjSMZmrDaKY0jEm3\nozm4YQxNY2uo8DfEN5QdO4LlazcXhEMHS1ZtYFvPDgAm1ddw3IxG/uCE6Rw7vZFjpjfQMKaq32Oc\nNHsCtz2xkqv+63nOvPoRvvT7x3DuWw8eiadj+5kiYs8dpNOATcAtBYFxBvCTiOiW9FWAiPiUpCOB\nW4GTgIOB/wYOSw/1AvB+oBV4AvhQRDwv6Q7gBxFxm6RvA4sj4rq9FX5ic3MsbGnZ92dcpjq39/T/\nwN/ctTMM+tq2sXZTF2s3b2NjZ/egjzNKML62mon11Uyoqwbg1Q1drFq/la7uHf36VowSk8fW7AyQ\nvmDJbqc2jGHy2BoqKw7sQ10RwYat3bRv6Xtd123pC9b2LdnrOam+Jj3v0RzcOGbnazCmumKkn0LR\n1mzsZPHKjiwgWrOQ2JDeO7XVFbx1egPHzmjkuOmNHDujkakNo5HyfYlY/tpmrrh9EYtWrueC46fx\nD+cdxbjRVXvf0EaUpIUR0VzUtnsLjPQLZgH39gbGgHUXAB+IiI+k0QUR8eW07gHgc6nr5yLizNT+\nmdT2FaANmJLC59TCfntSM3VOTP/YNYyuHMXoqgpGV1VQUzWK0ZUVjK7qaxud2mqqCtr79clua3az\nXf8+o3L/MW3d1sPagg/9nSHQ+yFVEALtm7exqWvwAKgcJcbXVTOxrjcEaphYl4XBhNQ+oWBd45iq\nQfcrRwTrt2xndUcnr2zYyuqOTlav79xleev2nn7bjRI0ja1hSsMYDi4Ik8Llg8aNpmo/hsr2nh2s\nS6/fugGv4y4/W7I+3TsGf5+PrhrFxLoa6moqWLspe6yBGmurmJqe79TGLER6w+Tgxuw1qKkc+VDZ\n3NXNMy/3hcOiFetZ1dEJZF8Ojpgytl84HDq5fsgjzO6eHfzLQ8v4xk+WMWXcaP7pD4/llEMmDsfT\nsRIZSmAMxzGMjwO3p/vTgMcK1rWmNoCVA9pPBiYC6yOie5D+u5A0D5gHMGHabC551yF0bt9B5/ae\n7La7h67e+9t72NC5vd/6ru09dHb3sL1n7yG5OzWV/YOmN3xqqiro3N6zMyAGfvD2qq4Y1fdhX1/N\nrIm1WQjUDxICdTWMG1OZO6T2RMqCZ3xdNUcePG7QPr3fxFenAHmlozMFyVZe2dDJ0jWbeOSFNjZv\n6xnw2NmujL5AGdNvlDK1YTSTx9UM+qEaEWzZ1rMzQNcVhGv75u3ptv9Ia3ejK8g+3CfUZq/fmybW\ncvzMxp2v68CfiXU1u4weOrf39D3vjux1WLU+3XZ0snDFOtZv2b7L751UX73zuWfBMqbf85/SMLyh\n2t2zg1+/upHFK7MD04tXdrB0zUZ6c3HmhFpOnDWBj09v4LgZjRx1cENJRkqVFaO4/H2H8a7Dmrji\n9kV86PrHmHfaIfzV+w8rixC14TWkwJD0WaAb+F5v0yDdgsGn78Ye+g8qIuYD8wGam5vj/555xD7V\n26tnR6QQ6aGze0ff/YJQ6RdEBW1dBX2ztr4+42urOXRyffrALxgJ1PeFQH3N8ARAKUiiobaKhtoq\njpgyeKgAbOzcnj5QO3mloy9cVnV0svy1zfz8N2sH/VDv/VBtGFPF+i3bdwbAtgG7ynpVVSh9uNcw\noa6KY8Y3MrGumvG12WvaGwwT67O28bVVQ959NrqqglmT6pg1qW63fbZs6y4YoW3dGS6r1nfy0trN\nPPbirs9fgqb6mixIxmUjlYMbxvQbsexu919E0LpuK0+l4w6LV67n2VUddG7PXrfxtVUcO6ORs46e\nwnEzstFD7+7I/eX4meO577J38oUfLuE7P32RR154jWsuPI7DDhq7X+uw0io6MCTNJTsYfnr07ddq\nBWYUdJsOrEr3B2t/DWiUVJlGGYX9S6ZilKirqaSuxpPEijF2dBVjR1ft8cNgU1d3+qbef7TySsdW\nOrZuZ8q40bxl6ridQTq+YGQ1oczDtba6kjc31fPmpvrd9tnU1c3q9VtZlUZohSOWpWs28sjSNrYM\nGKn1HlOamkYoTfU1vLR2M4tbO2hPu8pqKkdx9LQGPnzSmzhuZrZ7acaEMWXxOtVWV/KlC47hvYdP\n5lN3Pc2533jU029fZ4o6hpFmPP0z8K6IaCvodxTw7/Qd9H4QmEM2kngBOB14meyg94cj4jlJ3wfu\nKjjo/XREfGtvNTU3N0fL6+igt72xRAQbOruzEFnfyaqOrdkorWDU8uqGTqaPH7Nz1HDs9EYOnzJ2\nvx4vKtZrm7r49F1P899L1vD2Qyd6+m0ZKelBb0m3Au8GJgGvAlcCnwFqgLWp22MRcUnq/1my4xrd\nwOURcX9qPwf4Otm02psi4oup/RD6ptU+BfzviOjaW+EODLPyFhE7p99WV47iixcc7em3ZaDks6TK\nkQPD7MDg6bflZSiBUf5jWzM7oM2eVMedl5zK5e+bwz2LV3H213/G4y+u3fuGVnYcGGZWcr3Tb++8\n5FSqKsSF1z/Gl+9fQlf34NPPrTw5MMxsvzl+5nh+eOk7ufB3ZvKdn77I+d/8OS+8unGky7KcHBhm\ntl/V1VTy5d8/hhsuambNhk7O/caj3PTocnbs5mx8Kx8ODDMbEe878iB+dPlpvPPQSVx17/NcdNMv\neSVdysTKkwPDzEZM09gabpjbzJcuOIaFL63jzK8/wr1Pl/zcXSuSA8PMRpQkPnzyTH546TuYNamO\nv/j3p7ji9kVs6Nz1ml02shwYZlYWDmmq585LTuWy0z39tlw5MMysbFRVjOKK93v6bblyYJhZ2emb\nfjvD02/LiAPDzMpSNv32rVzv6bdlw4FhZmXt/Wn67Ts8/XbEOTDMrOw1ja3hxrnNfPGCoz39dgQ5\nMMzsgCCJj5z8Jk+/HUEODDM7oPROv720YPrtPYtX8cKrG9m6zbOpSsn/H4aZHbCeXLGOK25fxEtr\nt+xsaxpbw8wJtcycUMuMCbW8aUItMydmy031NW/4/y7W/4GSmb1hdW7vYcnqDaxo38LK9i2sSD8r\n27eyqmMrhR9xNZWjmJHCpN/PxFpmjK9lTHXFyD2R/WQogVE53MWYme1Po6sqOH7meI6fOX6XdV3d\nPaxa35mFyNrNO8NkRftWHn9xLZsH7MLy6GTPHBhm9rpVU1nB7El1zJ5UBzT1WxcRrNuyvWBEsoUV\na7fwUvtmfrm8nf9c9LJHJwM4MMzsDUkSE+qqmVBXzXEzGndZv617By+v37pLoKxo38Ivl7ezqau7\nX/+Bo5Pp48cweWwNTelnYl0NFQf4CGWvgSHpJuBcYE1EHJ3aJgC3A7OA3wJ/GBHrJAm4BjgH2AJ8\nNCKeTNvMBf4uPewXIuLm1H4i8F1gDHAfcFkcqAdWzOx1o7pyVMHopL/djU56w+TuRS8z8IT0UYIJ\ndVl4TKqv3hkkTfV9oTJ5bA2T6mtoGFNF9nFaXvKMML4L/AtwS0Hbp4EHI+Irkj6dlj8FnA3MST8n\nA9cBJ6eAuRJoBgJYKOmeiFiX+swDHiMLjLOA+4f+1MzMSiPP6OSVjk7aNnXRtrGr73Zj3/KLbZtp\n29jFtp4du2xfXTGqf6gUBMukgoBpGltDbfX+21G0198UEY9ImjWg+Tzg3en+zcDDZIFxHnBLGiE8\nJqlR0tTUd0FEtANIWgCcJelhYFxE/CK13wKcjwPDzA5g1ZWjsoPlE2v32C8i2NDZ3S9IBgbLy+s7\nWdzawdpNXbuMWgDqqiv6Bcik+v6jlsJdYtWVQzv1rthoOigiVgNExGpJk1P7NGBlQb/W1Lan9tZB\n2gclaR7ZaISZM2cWWbqZWXmQRMOYKhrGVHHo5Po99u3ZEbRv3rbbYGnb2MkLr27if5atpWPr4Ge/\nj6+tGlK9wz2WGWynWxTRPqiImA/Mh+w8jGIKNDM7EFWM0s7Rwt50bu9hbW+49AuWThYNoYZiA+NV\nSVPT6GIqsCa1twIzCvpNB1al9ncPaH84tU8fpL+ZmRVpdFUF0xrHMK1xzC7rvjiExy12h9Y9wNx0\nfy5wd0H7RcqcAnSkXVcPAGdIGi9pPHAG8EBat1HSKWmG1UUFj2VmZmUkz7TaW8lGB5MktZLNdvoK\ncIeki4EVwAdT9/vIptQuI5tW+zGAiGiX9HngidTvqt4D4MCf0Tet9n58wNvMrCz5WlJmZm8gQ7mW\nlC9vbmZmuTgwzMwsFweGmZnl4sAwM7NcHBhmZpaLA8PMzHJxYJiZWS4ODDMzy8WBYWZmuTgwzMws\nFweGmZnl4sAwM7NcHBhmZpaLA8PMzHJxYJiZWS4ODDMzy8WBYWZmuTgwzMwsFweGmZnl4sAwM7Nc\nhhQYkq6Q9JykZyXdKmm0pNmSHpe0VNLtkqpT35q0vCytn1XwOJ9J7b+WdObQnpKZmZVC0YEhaRpw\nKdAcEUcDFcCFwFeBqyNiDrAOuDhtcjGwLiIOBa5O/ZB0ZNruKOAs4FuSKoqty8zMSmOou6QqgTGS\nKoFaYDXwXuDOtP5m4Px0/7y0TFp/uiSl9tsioisilgPLgJOGWJeZmQ2zogMjIl4GvgasIAuKDmAh\nsD4iulO3VmBauj8NWJm27U79Jxa2D7JNP5LmSWqR1NLW1lZs6WZmVoSh7JIaTzY6mA0cDNQBZw/S\nNXo32c263bXv2hgxPyKaI6K5qalp34s2M7OiDWWX1PuA5RHRFhHbgR8AbwMa0y4qgOnAqnS/FZgB\nkNY3AO2F7YNsY2ZmZWIogbECOEVSbToWcTrwPPAQ8IHUZy5wd7p/T1omrf9JRERqvzDNopoNzAF+\nOYS6zMysBCr33mVwEfG4pDuBJ4Fu4ClgPvBD4DZJX0htN6ZNbgT+TdIyspHFhelxnpN0B1nYdAOf\niIieYusyM7PSUPYl/8DT3NwcLS0tI12GmdkBRdLCiGguZluf6W1mZrk4MMzMLBcHhpmZ5eLAMDOz\nXBwYZmaWiwPDzMxycWCYmVkuDgwzM8vFgWFmZrk4MMzMLBcHhpmZ5eLAMDOzXBwYZmaWiwPDzMxy\ncWCYmVkuDgwzM8vFgWFmZrk4MMzMLBcHhpmZ5eLAMDOzXIYUGJIaJd0p6VeSlkg6VdIESQskLU23\n41NfSbpW0jJJT0s6oeBx5qb+SyXNHeqTMjOz4TfUEcY1wI8i4gjgWGAJ8GngwYiYAzyYlgHOBuak\nn3nAdQCSJgBXAicDJwFX9oaMmZmVj6IDQ9I44DTgRoCI2BYR64HzgJtTt5uB89P984BbIvMY0Chp\nKnAmsCAi2iNiHbAAOKvYuszMrDSGMsI4BGgD/lXSU5JukFQHHBQRqwHS7eTUfxqwsmD71tS2u/Zd\nSJonqUVSS1tb2xBKNzOzfTWUwKgETgCui4jjgc307X4ajAZpiz2079oYMT8imiOiuampaV/rNTOz\nIRhKYLQCrRHxeFq+kyxAXk27mki3awr6zyjYfjqwag/tZmZWRooOjIh4BVgp6fDUdDrwPHAP0DvT\naS5wd7p/D3BRmi11CtCRdlk9AJwhaXw62H1GajMzszJSOcTt/xL4nqRq4EXgY2QhdIeki4EVwAdT\n3/uAc4BlwJbUl4hol/R54InU76qIaB9iXWZmNswUMejhgrLX3NwcLS0tI12GmdkBRdLCiGguZluf\n6W1mZrk4MMzMLBcHhpmZ5eLAMDOzXBwYZmaWiwPDzMxycWCYmVkuDgwzM8vFgWFmZrk4MMzMLBcH\nhpmZ5eLAMDOzXBwYZmaWiwPDzMxycWCYmVkuDgwzM8vFgWFmZrk4MMzMLBcHhpmZ5eLAMDOzXIYc\nGJIqJD0l6d60PFvS45KWSrpdUnVqr0nLy9L6WQWP8ZnU/mtJZw61JjMzG37DMcK4DFhSsPxV4OqI\nmAOsAy5O7RcD6yLiUODq1A9JRwIXAkcBZwHfklQxDHWZmdkwGlJgSJoO/C5wQ1oW8F7gztTlZuD8\ndP+8tExaf3rqfx5wW0R0RcRyYBlw0lDqMjOz4TfUEcbXgb8BdqTlicD6iOhOy63AtHR/GrASIK3v\nSP13tg+yTT+S5klqkdTS1tY2xNLNzGxfFB0Yks4F1kTEwsLmQbrGXtbtaZv+jRHzI6I5Ipqbmpr2\nqV4zMxuayiFs+3bg9ySdA4wGxpGNOBolVaZRxHRgVerfCswAWiVVAg1Ae0F7r8JtzMysTBQ9woiI\nz0TE9IiYRXbQ+icR8RHgIeADqdtc4O50/560TFr/k4iI1H5hmkU1G5gD/LLYuszMrDSGMsLYnU8B\nt0n6AvAUcGNqvxH4N0nLyEYWFwJExHOS7gCeB7qBT0RETwnqMjOzIVD2Jf/A09zcHC0tLSNdhpnZ\nAUXSwohoLmZbn+ltZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFh\nZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZ\nmeXiwDAzs1yKDgxJMyQ9JGmJpOckXZbaJ0haIGlpuh2f2iXpWknLJD0t6YSCx5qb+i+VNHfoT8vM\nzIbbUEYY3cBfR8RbgFOAT0g6Evg08GBEzAEeTMsAZwNz0s884DrIAga4EjgZOAm4sjdkzMysfBQd\nGBGxOiKeTPc3AkuAacB5wM2p283A+en+ecAtkXkMaJQ0FTgTWBAR7RGxDlgAnFVsXWZmVhrDcgxD\n0izgeOBx4KCIWA1ZqACTU7dpwMqCzVpT2+7aB/s98yS1SGppa2sbjtLNzCynIQeGpHrgLuDyiNiw\np66DtMUe2ndtjJgfEc0R0dzU1LTvxZqZWdGGFBiSqsjC4nsR8YPU/Gra1US6XZPaW4EZBZtPB1bt\nod3MzMrIUGZJCbgRWBIR/1yw6h6gd6bTXODugvaL0mypU4COtMvqAeAMSePTwe4zUpuZmZWRyiFs\n+3bgj4FnJC1KbX8LfAW4Q9LFwArgg2ndfcA5wDJgC/AxgIhol/R54InU76qIaB9CXWZmVgKKGPRw\nQdlrbm6OlpaWkS7DzOyAImlhRDQXs63P9DYzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFhZma5\nODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZmeXi\nwDAzs1wcGGZmlosDw8zMcnFgmJlZLmUTGJLOkvRrScskfXqk6zEzs/7KIjAkVQDfBM4GjgQ+JOnI\nka3KzMwKlUVgACcByyLixYjYBtwGnDfCNZmZWYHKkS4gmQasLFhuBU4e2EnSPGBeWtwu6en9UNu+\naAA6RrqIAVxTfuVYl2vKxzXlN6fYDcslMDRIW+zSEDEfmA8gqS0imktd2L6QND8i5u295/7jmvIr\nx7pcUz6uKT9J84vdtlx2SbUCMwqWpwOr9rLN+tKVU7T/GukCBuGa8ivHulxTPq4pv6LrUsQuX+T3\nO0mVwAvA6cDLwBPAhyPiuT1s01JuIwwzs9ezstglFRHdkv4CeACoAG7aU1gkRQ+rzMxs35XFCMPM\nzMpfuRzDMDOzMufA2AeSeiQtKviZtYe+75Z0736oKST9W8FypaS2/fG790bSBam+I0a4jrJ9jQAk\nbRrpGnZnb7VJelhSyY8llst7aSBJn5X0nKSn02fCLqcDjEBN0yXdLWmppN9IukZS9R76Xy6pNs9j\nl31glNkf09aIOK7g57cjXRCwGTha0pi0/H6yiQO5pUkHpfAh4FHgwn3ZKJ35P5yG/BrZiCvqvVRK\nkk4FzgVOiIi3Au+j//lkI1GTgB8A/xkRc4DDgHrgi3vY7HLg9REY5U5ShaR/lPRE+pbxpwWrx0n6\nD0nPS/q2pFK93vcDv5vufwi4taC+kyT9XNJT6fbw1P5RSd+X9F/Aj4e7IEn1wNuBi0l/5GnU9chg\nr4mkTZKukvQ4cOpw10Nxr9HPJB1X0O9/JL21BLXtMiKV9C+SPpru/1bSP0h6UtIz+/tb9p5q20+/\nf3fvpd29XudI+pWkRyVdW8KR5FTgtYjoAoiI1yJilaQTJf1U0kJJD0iamup6WNLX03vsWUknlaCm\n9wKdEfGvqaYe4Arg45LqJH0tvYeelvSXki4FDgYekvTQ3h78gAgMSfWSHiz4gzkvtc+StETS9WlY\n+OOCb5GlMEZ9u6P+I7VdDHRExO8AvwP8iaTZad1JwF8DxwBvBn6/RHXdBlwoaTTwVuDxgnW/Ak6L\niOOBvwe+VLDuVGBuRLy3BDWdD/woIl4A2iWdkNp395rUAc9GxMkR8WgJ6inmNboB+CiApMOAmogY\nqasLvBYRJwDXAZ8coRpGyu7eS7tI/77fAc6OiHcATSWs68fADEkvSPqWpHdJqgK+AXwgIk4EbqL/\nt/u6iHgb8Odp3XA7ClhY2BARG4AVwP8BZgPHpxHR9yLiWrJz3t4TEe/Z24MfEIEBdAIXpD+Y9wD/\nlIZekJ3m/s2IOIrsZL4/KGEdhbukLkhtZwAXSVpE9iE0kb5T73+Zro/VQ/aN9h2lKCp9iM0i++Z8\n34DVDcD3JT0LXE32huq1ICLaS1FTquW2dP+2tAy7f016gLtKVEuxr9H3gXPTh8DHge+Wqr4cfpBu\nF5I9jzeS3b2XBnME8GJELE/Lt+6h75BExCbgRLLLFbUBtwN/ChwNLEifCX9HdiJyr1vTto+Q7YFo\nHOayxCBXyUjtpwHfjojuVMM+/+2XxXkYOQj4kqTTgB1k1546KK1bHhGL0v2R+GMS8JcR8UC/Rund\n7PoPV8o5zPcAXwPeTRZavT4PPBQRFyg7SP9wwbrNpShE0kSyofHRkoLs3Jog+6De3WvSmUKklPbp\nNYqILZIWkF0I8w+BUh7c7ab/F7jRA9Z3pdse9v/f7d5qK5k9vJfu2U1Ng11mqGTSe/Zh4GFJzwCf\nAJ6LiN3tVi31Z8JzDPjSLGkc2ZU0Xhzq7ztQRhgfIRtanhgRxwGv0vcG6SroNxJ/TA8Af5a+hSLp\nMEl1ad1Jkman/fR/RHbQrlRuAq6KiGcGtDfQd4D3oyX8/YU+ANwSEW+KiFkRMQNYTjaa2J+vyUDF\nvEY3ANcCT5RwNAbwEnCkpBpJDWRXPSgXI1nb7t5L7KamXwGHqG8G4x+VqjBJh0sqvJDfccASoEnZ\nAXEkVUkqHNX/UWp/B9mu7OG+OOGDQK2ki9LvqQD+iWx0/GPgEqVJLpImpG02AmPzPPiBEhgNwJqI\n2C7pPcCbRrqgAjcAzwNPpl0a36EvtH4BfAV4luxN/h+DPsIwiIjWiLhmkFX/D/iypP8h+3a2P3yI\nXZ/rXcCH2Y+vyUDFvEYRsRDYAPxrKWpKf7xdEbESuAN4Gvge8FQpft++KJPa9vRe2qWmiNhKdnzg\nR5IeJftyWaorxtYDNyubwPE02f/l8/dkIfdVSYuBRcDbCrZZJ+nnwLfJjn8Oq8jOxL4A+KCkpWSX\nXOoE/pbss2oF8HSq7cNps/nA/XkOepf1md7pDfsqcDjZBbOqyP4B3k72ny0B3BsRR6f+nwTqI+Jz\n+79a25u0m+6TEXHuSNeSl6SDyXY5HBERO0rw+McC10dEKWbMDEk517YnkuojYlM6zvlNYGlEXF0G\ndT1M9v5vGelailXuxzCOAn6dlg/pAAACC0lEQVQTEa+x+6mWR/feiYiv7Zeq7A0hDeu/CPxVicLi\nEuBSsnnwZaWca8vhTyTNBarJRh7fGeF6XjfKdoRR+IaNiGE/T8DMzPZN2QaGmZmVlwPloLeZmY2w\nsgkMSTMkPZTO3H5O0mWpfYKkBcoupLVA0vjUfoSkX0jqSge7Cx/rsnTq/XOSDsR9sGZmZadsAoPs\n5KC/joi3AKcAn5B0JPBp4MF0Ia0H0zJAO9kxjn4HuiUdDfwJ2SUojiU7U7fo//TczMwyZRMYEbE6\nIp5M9zeSnQAzjews25tTt5vJritDRKyJiCeA7QMe6i3AYxGxJZ0C/1OyeclmZjYEZRMYhdJZmseT\nXZvpoIhYDVmoAJP3svmzwGmSJiq7xvs5ZKfFm5nZEJTdeRjKLmV8F9l02g191xjMJyKWSPoqsADY\nBCwm291lZmZDUFYjjHQ9prvILrvbe3XOV9V3PfmpwJq9PU5E3BgRJ0TEaWTHOpaWqmYzszeKsgmM\ndBr/jcCSiPjnglX3AHPT/bnA3Tkea3K6nUn2/y2U7BLHZmZvFGVz4l66euPPgGfILmEO2QWzHie7\nyNhMsgtnfTAi2iVNAVqAcan/JuDItBvrZ2SXr95OdlmHB/frkzEzex0qm8AwM7PyVja7pMzMrLw5\nMMzMLBcHhpmZ5eLAMDOzXBwYZmaWiwPDzMxycWCYmVku/x+C0P8mZYWOXAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e7954e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.Y.plot()\n",
    "plt.ylim([0,15000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Chow test\n",
    "\n",
    "https://en.wikipedia.org/wiki/Chow_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_rss(df):\n",
    "    # linear regression of the input series\n",
    "    # return the residual sum of squares\n",
    "    rss = np.linalg.lstsq(df[['X','b']].values, df.Y.values, rcond=-1)[1]\n",
    "    return rss, len(df)\n",
    "\n",
    "def chow(df, k=2, change_point='2019-07-01'):\n",
    "    # find the two sub-series\n",
    "    df_1 = df[df.index < change_point]\n",
    "    df_2 = df[df.index >= change_point]\n",
    "    \n",
    "    # linear regressions for the total series\n",
    "    rss, N = find_rss(df)\n",
    "    # and the two sub-series\n",
    "    rss_1, N_1 = find_rss(df_1)\n",
    "    rss_2, N_2 = find_rss(df_2)\n",
    "    # degrees-of-freedom for the sub-series\n",
    "    k_2 = (N_1 + N_2 - 2*k)\n",
    "    # Chow statistic\n",
    "    C = ((rss - (rss_1 + rss_2)) / (rss_1 + rss_2) * (k_2/k))[0] \n",
    "    # p-value from the F-statistic\n",
    "    p_value = f.sf(C, k, k_2)\n",
    "    \n",
    "    return C, k, k_2, p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 263,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P-value: 0.0024087096046\n"
     ]
    }
   ],
   "source": [
    "C, k_1, k_2, p = chow(df, change_point='2019-07-01')\n",
    "print(\"P-value:\", p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Show how the Chow test performs under the null hypothesis of a single generating distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_null():\n",
    "    Y = [norm(12000, 200).rvs() for _ in idx]\n",
    "\n",
    "    df_null = pd.DataFrame(Y, index=idx, columns=['Y'])\n",
    "    df_null['b'] = 1\n",
    "    df_null['X'] = [m.month for m in s.index]\n",
    "    return df_null"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df_null = generate_null()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 15000)"
      ]
     },
     "execution_count": 270,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEHCAYAAAC9TnFRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xuc3XV95/HXO3PNJDOZBCYJZIJB\njYRbUJgGrF2q0nKr22BXt2i3RMs2vbhVu/XRat1H6art6taKsluxKWChDxfESxfqajFFrLUVygQk\nIQSSCJYMIclgkpnc5v7ZP37fmTmZzOWXc+ZkzsD7+XjM4/x+39/3d+Zzzpxz3uf7u40iAjMzs6nM\nmekCzMxsdnBgmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcqme6QKK\ndfrpp8eKFStmugwzs1ll06ZNL0VESzHrztrAWLFiBe3t7TNdhpnZrCLp34pd15ukzMwsFweGmZnl\n4sAwM7NcHBhmZpaLA8PMzHJxYJiZWS4ODDMzy8WBYWZmuTgwzMwslykDQ9IdkvZJenKcZR+SFJJO\nT/OSdIuknZI2S7q4oO86STvSz7qC9kskbUnr3CJJ0/XgzMxs+uQZYfw1cPXYRknLgZ8Hni9ovgZY\nmX7WA7emvouAm4BLgTXATZIWpnVuTX2H1zvhd5mZ2cybMjAi4nvA/nEW3Qz8PhAFbWuBuyLzMNAs\n6QzgKmBjROyPiAPARuDqtKwpIn4QEQHcBVxX2kMyM7NyKGofhqRfBF6IiCfGLFoG7CqY70htk7V3\njNM+0e9dL6ldUntnZ2cxpZuZWZFOOjAkNQAfBf5ovMXjtEUR7eOKiA0R0RYRbS0tRV2d18zMilTM\nCOM1wNnAE5J+DLQCj0laSjZCWF7QtxXYPUV76zjtZmZWYU46MCJiS0QsjogVEbGC7EP/4ojYA9wP\n3JCOlroM6IqIF4EHgCslLUw7u68EHkjLDkm6LB0ddQNw3zQ9NjMzm0Z5Dqu9G/gBcI6kDkk3TtL9\nm8CzwE7gr4DfBoiI/cDHgUfTz8dSG8BvAbeldX4EfKu4h2JmZuWk7OCk2aetrS38H/fMzE6OpE0R\n0VbMuj7T28zMcnFgmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFg\nmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJm9QpT6H1arp6kO\nMzOrIMf6Btm+9xDP7DnEtj3dPP3iIZ7e013SfTowzMxmsaGhoOPAsZFQeGZvdvvcT44wPKCYW1PF\nOUsbuer8pfywhN81ZWBIugN4G7AvIi5IbX8G/HugD/gR8N6IOJiWfQS4ERgE3h8RD6T2q4HPAVXA\nbRHxydR+NnAPsAh4DPjViOgr4TGZmb0sdR3r55k92UhhWxoxbN9ziCN9gwBI8KpFDaxa2sQvvv5M\nVi1tZNXSJs5a1MCcOQLgUyX8fk21TUvS5cBh4K6CwLgS+E5EDEj6FEBE/IGk84C7gTXAmcA/AK9L\nd7Ud+HmgA3gUeFdEPCXpXuDrEXGPpC8AT0TErVMV3tbWFu3t7Sf/iM3MKtzA4BDPvXSEbXsO8fSL\n3Tydbnd39Yz0WTC3hlVLGzn3jCZWLW3knKWNvG5JI/PqJh8HSNoUEW3F1DXlCCMividpxZi2bxfM\nPgy8I02vBe6JiF7gOUk7ycIDYGdEPJsKvgdYK2kb8Fbg3anPncAfA1MGRv/gEBGBpKm6vuL09A+y\n7cVuNnd0sae7hzmCORKSqJKy+TlCqb1w+ZyCtmx+tP/xfYf7FfZNbXOOX7eq8L7niKo5Ym5NFQ21\nVcytraKhtpqGmqqRb0AvJxHBsf5Buo8N0N3Tz6Ge/pHp7mP9dPcMTx/fdqinn0M9AzTWV7O4sY4l\nTfUjty2NdSxurGdJUx2Lm+qZP8UHxMvd0FAQQNUsfP1EBJ2He7NNSQX7GnbuO0zf4BAA1XPEa1rm\n81NnL2LV0iZWndHIuUubWNJUd8o//6bjlfZrwJfT9DKyABnWkdoAdo1pvxQ4DTgYEQPj9D+BpPXA\neoDapa/lkk/8AxcsW8CFy5q4cNkCLli2gGXNc19RITIwOMT2vYfZ3HGQJzq62NxxkGf2HGJgKBs5\nVqc30VAEQ6UdIFF2ddVzaBgOkNrjA2VubRUNIyEzdnkVc2tOXGd4vqG2uugPk6Gh4EjfQPbBfuz4\nD/SRD/xj/SkMTvzwP9QzMPK3mOxxN9bX0DS3mqb6GhbMrWH5wrnMr6vmUM8Ae7t7ePz5g+zt7qF3\nYOiE9efVVrE4BcripnqWNNaxuCkLleHbJU11zK+rrvj3xsDgEF3H+jlwtI8DR/s5cKRvdPpoX5rv\n5+DRPvYf6ePg0X4OHutncCioq57DvLrs7z4/3c4ruJ1XW01DXVV2O9ynrpp56TUyr+74fg01VVRX\nTd+BpD39g+zYe/iEfQ0/OTK6BX5JUx3nLG3i3608nVVnZJuTXtMyn9rqyjigtaTAkPRRYAD40nDT\nON2C8Q/fjUn6jysiNgAbAM4+d3VcsWoxW17o4gs7X2IwvSkXNtSkEFkwEiKtC18eITI0FDz3kyNs\n7jjI5o4uNnd0sXV3Fz392YdIU301q1ubWX/5q1nd2sxFyxewtKn+uMc+NBQj4TEUQcRwmGRtUbCs\ncPngUGHf4XUL+g6Nd39p+dDx990/NERP3yBH+gY51jfA0b5BjvYNcqx/kKNp/thwW98g+w71ZH16\ns+XH+gfpHzy59KsdDqOaMSGUfuqqqzjcO3DCCOBQ7wBTHYnYUFtFU30NjfXVNM2t4fT5tby6Zd5x\nbU0FgVDY1lhfTX1NVa7HEBF09wywr7uHfYd62TvmtrO7l80dWbAMvyYKza2pykYlY4KkcLqlsZ6m\n+ukJlt6BQQ6mD/rhD/fCD/3xAqC7Z2DC+6utmsPCeTUsbKiluaGGc5Y20txQy8KGGmqrqjjaN8CR\nvgGO9g5mt32DHO4dYF9373HzfeOE7kTqquekYBkNmhPCJ93OGw6fdDswFGzfcyjbnLSnm+deOjLy\npa2+Zg7nLGnk585dwjlLG0fCYdG82lKf9rIqOjAkrSPbGX5FjO4I6QCWF3RrBXan6fHaXwKaJVWn\nUUZh/0mdNq+WP3vnRcDoJpgnX+hiywtdbHmhmw3fe3bkm91wiBQGSaWHSESwu6uHzbtGRw5bXuji\nUHpDza2p4oJlTbx7zau4aPkCVrc2s+K0hikf05w5Ys64OT279A8OFQTLQEHgHB9CYwMo6ze6/CeH\n+9jVN0DvwBDz67IP9DOb61lV35g+1KuP+/Zf+EHfNDe7rZnGb6GTkcSCudkIZOWSxgn7RQSH0gfl\n2HAZnt66u5vvdO/jaNpZWqi+Zs5omDRmm8CWjIxg6gAmDoCjfRw4kgXAkXHue1hDbRULG2pHAuCs\nRQ0sbKihuaGWRfOyQFg4Zrqhtmpa3rPDr52jfQMc6R3kSO/xQXOkd3TZcAgN9xsOnb3dPQXLBycN\nobMWNbBqaSO/sHp4J3Qjrzpt3qzchDblTm+AtA/jGwU7va8GPgP8bER0FvQ7H/g/jO70fhBYSTaS\n2A5cAbxAttP73RGxVdJXgK8V7PTeHBGfn6qmqXZ69/QP8vSeQ2x5oYsnO7Ig2b53dFNNc0MNF5x5\nfIgsXzRzIfLS4ezb4RO7RsPhpcPZULWmSqxa2sTq1gVc1NrM6uULeG3L/GkdLtsr0/CH377uXvYd\nGr3de9x8L4d7J/7m31RfzcJ5tVkApA/3bH7iAKirzjeqmi36B4cKRjZZwATw2sXzK24fUyk7vfMc\nJXU38GbgdGAvcBPwEaAO+Enq9nBE/Gbq/1Gy/RoDwAcj4lup/Vrgs2SH1d4REX+S2l/N6GG1jwP/\nKe00n1QxR0n19A/yzHCIpNFI4fb+BXNruGBZ00iIrF7WXJYQ6e7p58mOrpGRw+aOLl44eAzIDotb\nuXh+tkmpdQEXtjZz7hmNL7s3mM0uR3oHRkYnVXM0EgbNc2v8xWWWKWtgVKrpOqy2d+D4ENnckY1E\nhreRjw2RC5ct4KxFU2/6GXasb5CnXuwaGTls7uji2ZeOjCw/a1HD6MihNRvxTHVYnJlZscp6WO3L\nXV11Fatbm1nd2jzSNjZEtrzQxR3ff24kRJrqq0cCZPj2Vac1MDAUPLPnEE90HGRLGkFs33toZIf8\nkqY6Vrc280sXL2N1azMXLlvAwgrfyWVmNuwVHxjjmShEtu85nHaqZ0HyxX/+8cix0o311fQODI3s\n/GpuqGF1azNXrFqcjSCWN7OkqX5GHo+Z2XRwYORUV13Fha0LuLB1wUhb38AQ2/eOjkQaaqvSvofy\n7PswM5tJDowS1FbPGTlc18zs5c6HN5iZWS4ODDMzy8WBYWZmuTgwzMwsFweGmZnl4sAwM7NcHBhm\nZpaLA8PMzHJxYJiZWS4ODDMzy8WBYWZmuTgwzMwsFweGmZnl4sAwM7NcHBhmZpaLA8PMzHKZMjAk\n3SFpn6QnC9oWSdooaUe6XZjaJekWSTslbZZ0ccE661L/HZLWFbRfImlLWucW+d/UmZlVpDwjjL8G\nrh7T9mHgwYhYCTyY5gGuAVamn/XArZAFDHATcCmwBrhpOGRSn/UF6439XWZmVgGmDIyI+B6wf0zz\nWuDONH0ncF1B+12ReRholnQGcBWwMSL2R8QBYCNwdVrWFBE/iIgA7iq4LzMzqyDF7sNYEhEvAqTb\nxal9GbCroF9HapusvWOc9nFJWi+pXVJ7Z2dnkaWbmVkxpnun93j7H6KI9nFFxIaIaIuItpaWliJL\nNDOzYhQbGHvT5iTS7b7U3gEsL+jXCuyeor11nHYzM6swxQbG/cDwkU7rgPsK2m9IR0tdBnSlTVYP\nAFdKWph2dl8JPJCWHZJ0WTo66oaC+zIzswpSPVUHSXcDbwZOl9RBdrTTJ4F7Jd0IPA+8M3X/JnAt\nsBM4CrwXICL2S/o48Gjq97GIGN6R/ltkR2LNBb6VfszMrMIoOzhp9mlra4v29vaZLsPMbFaRtCki\n2opZ12d6m5lZLg4MMzPLxYFhZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4M\nMzPLxYFhZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFhZma5ODDM\nzCyXkgJD0u9K2irpSUl3S6qXdLakRyTtkPRlSbWpb12a35mWryi4n4+k9mckXVXaQzIzs3IoOjAk\nLQPeD7RFxAVAFXA98Cng5ohYCRwAbkyr3AgciIjXAjenfkg6L613PnA18HlJVcXWZWZm5VHqJqlq\nYK6kaqABeBF4K/DVtPxO4Lo0vTbNk5ZfIUmp/Z6I6I2I54CdwJoS6zIzs2lWdGBExAvAp4HnyYKi\nC9gEHIyIgdStA1iWppcBu9K6A6n/aYXt46xzHEnrJbVLau/s7Cy2dDMzK0Ipm6QWko0OzgbOBOYB\n14zTNYZXmWDZRO0nNkZsiIi2iGhraWk5+aLNzKxopWyS+jnguYjojIh+4OvATwPNaRMVQCuwO013\nAMsB0vIFwP7C9nHWMTOzClFKYDwPXCapIe2LuAJ4CngIeEfqsw64L03fn+ZJy78TEZHar09HUZ0N\nrAT+tYS6zMysDKqn7jK+iHhE0leBx4AB4HFgA/D/gHskfSK13Z5WuR34G0k7yUYW16f72SrpXrKw\nGQDeFxGDxdZlZmbloexL/uzT1tYW7e3tM12GmdmsImlTRLQVs67P9DYzs1wcGGZmlosDw8zMcnFg\nmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFh\nZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLiUFhqRmSV+V9LSkbZLeKGmRpI2S\ndqTbhamvJN0iaaekzZIuLrifdan/DknrSn1QZmY2/UodYXwO+PuIWAVcBGwDPgw8GBErgQfTPMA1\nwMr0sx64FUDSIuAm4FJgDXDTcMiYmVnlKDowJDUBlwO3A0REX0QcBNYCd6ZudwLXpem1wF2ReRho\nlnQGcBWwMSL2R8QBYCNwdbF1mZlZeZQywng10Al8UdLjkm6TNA9YEhEvAqTbxan/MmBXwfodqW2i\n9hNIWi+pXVJ7Z2dnCaWbmdnJKiUwqoGLgVsj4g3AEUY3P41H47TFJO0nNkZsiIi2iGhraWk52XrN\nzKwEpQRGB9AREY+k+a+SBcjetKmJdLuvoP/ygvVbgd2TtJuZWQUpOjAiYg+wS9I5qekK4CngfmD4\nSKd1wH1p+n7ghnS01GVAV9pk9QBwpaSFaWf3lanNzMwqSHWJ6/8O8CVJtcCzwHvJQuheSTcCzwPv\nTH2/CVwL7ASOpr5ExH5JHwceTf0+FhH7S6zLzMymmSLG3V1Q8dra2qK9vX2myzAzm1UkbYqItmLW\n9ZneZmaWiwPDzMxycWCYmVkuDgwzM8vFgWFmZrk4MMzMLBcHhpmZ5eLAMDOzXBwYZmaWiwPDzMxy\ncWCYmVkuDgwzM8vFgWFmZrk4MMzMLBcHhpmZ5eLAMDOzXBwYZmaWiwPDzMxycWCYmVkuDgwzM8ul\n5MCQVCXpcUnfSPNnS3pE0g5JX5ZUm9rr0vzOtHxFwX18JLU/I+mqUmsyM7PpNx0jjA8A2wrmPwXc\nHBErgQPAjan9RuBARLwWuDn1Q9J5wPXA+cDVwOclVU1DXWZmNo1KCgxJrcAvALeleQFvBb6autwJ\nXJem16Z50vIrUv+1wD0R0RsRzwE7gTWl1GVmZtOv1BHGZ4HfB4bS/GnAwYgYSPMdwLI0vQzYBZCW\nd6X+I+3jrHMcSesltUtq7+zsLLF0MzM7GUUHhqS3AfsiYlNh8zhdY4plk61zfGPEhohoi4i2lpaW\nk6rXzMxKU13Cum8CflHStUA90EQ24miWVJ1GEa3A7tS/A1gOdEiqBhYA+wvahxWuY2ZmFaLoEUZE\nfCQiWiNiBdlO6+9ExK8ADwHvSN3WAfel6fvTPGn5dyIiUvv16Siqs4GVwL8WW5eZmZVHKSOMifwB\ncI+kTwCPA7en9tuBv5G0k2xkcT1ARGyVdC/wFDAAvC8iBstQl5mZlUDZl/zZp62tLdrb22e6DDOz\nWUXSpohoK2Zdn+ltZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFh\nZma5ODDMzCwXB4aZmeXiwDAzs1wcGGZmlosDw8zMcnFgmJlZLg4MMzPLxYFhZma5ODDMzCwXB4aZ\nmeXiwDAzs1yKDgxJyyU9JGmbpK2SPpDaF0naKGlHul2Y2iXpFkk7JW2WdHHBfa1L/XdIWlf6wzIz\ns+lWyghjAPi9iDgXuAx4n6TzgA8DD0bESuDBNA9wDbAy/awHboUsYICbgEuBNcBNwyFjZmaVo+jA\niIgXI+KxNH0I2AYsA9YCd6ZudwLXpem1wF2ReRholnQGcBWwMSL2R8QBYCNwdbF1mZlZeUzLPgxJ\nK4A3AI8ASyLiRchCBVicui0DdhWs1pHaJmof7/esl9Quqb2zs3M6Sjczs5xKDgxJ84GvAR+MiO7J\nuo7TFpO0n9gYsSEi2iKiraWl5eSLNTOzopUUGJJqyMLiSxHx9dS8N21qIt3uS+0dwPKC1VuB3ZO0\nm5lZBSnlKCkBtwPbIuIzBYvuB4aPdFoH3FfQfkM6WuoyoCttsnoAuFLSwrSz+8rUZmZmFaS6hHXf\nBPwqsEXSD1PbHwKfBO6VdCPwPPDOtOybwLXATuAo8F6AiNgv6ePAo6nfxyJifwl1mZlZGShi3N0F\nFa+trS3a29tnugwzs1lF0qaIaCtmXZ/pbWZmuTgwzMwsFweGmZnl4sAwM7NcHBhmZpaLA8PMzHJx\nYJiZWS4ODDMzy8WBYWZmuTgwzMwsFweGmZnl4sAwM7NcHBhmZpaLA8PMzHJxYJiZWS4ODDMzy8WB\nYWZmuTgwzMwsFweGmZnl4sAwM7NcKiYwJF0t6RlJOyV9eKbrMTOz41VEYEiqAv4CuAY4D3iXpPNm\ntiozMytUEYEBrAF2RsSzEdEH3AOsneGazMysQPVMF5AsA3YVzHcAl47tJGk9sD7N9kvafApqOxkL\ngK6ZLmIM15RfJdblmvJxTfmtLHbFSgkMjdMWJzREbAA2AEjqjIi2chd2MiRtiIj1U/c8dVxTfpVY\nl2vKxzXlJ2lDsetWyiapDmB5wXwrsHuKdQ6Wr5yi/d1MFzAO15RfJdblmvJxTfkVXZciTvgif8pJ\nqga2A1cALwCPAu+OiK2TrNNeaSMMM7OXs4rYJBURA5L+C/AAUAXcMVlYJEUPq8zM7ORVxAjDzMwq\nX6XswzAzswrnwDgJkgYl/bDgZ8Ukfd8s6RunoKaQ9DcF89WSOk/F756KpLen+lbNcB0V+xwBSDo8\n0zVMZKraJH1XUtn3JVbKa2ksSR+VtFXS5vSZcMLpADNQU6uk+yTtkPQjSZ+TVDtJ/w9Kashz3xUf\nGBX2ZjoWEa8v+PnxTBcEHAEukDQ3zf882YEDuaWDDsrhXcD3getPZqV05v90Kvk5shlX1GupnCS9\nEXgbcHFErAZ+juPPJ5uJmgR8Hfi/EbESeB0wH/iTSVb7IPDyCIxKJ6lK0p9JejR9y/iNgsVNkv5W\n0lOSviCpXM/3t4BfSNPvAu4uqG+NpH+R9Hi6PSe1v0fSVyT9HfDt6S5I0nzgTcCNpDd5GnV9b7zn\nRNJhSR+T9Ajwxumuh+Keo3+S9PqCfv8saXUZajthRCrpf0t6T5r+saT/LukxSVtO9bfsyWo7Rb9/\notfSRM/XtZKelvR9SbeUcSR5BvBSRPQCRMRLEbFb0iWS/lHSJkkPSDoj1fVdSZ9Nr7EnJa0pQ01v\nBXoi4ouppkHgd4FfkzRP0qfTa2izpN+R9H7gTOAhSQ9NdeezIjAkzZf0YMEbZm1qXyFpm6S/SsPC\nbxd8iyyHuRrdHPW3qe1GoCsifgr4KeDXJZ2dlq0Bfg+4EHgN8Etlquse4HpJ9cBq4JGCZU8Dl0fE\nG4A/Av60YNkbgXUR8dYy1HQd8PcRsR3YL+ni1D7RczIPeDIiLo2I75ehnmKeo9uA9wBIeh1QFxEz\ndXWBlyLiYuBW4EMzVMNMmei1dIL09/1L4JqI+BmgpYx1fRtYLmm7pM9L+llJNcD/At4REZcAd3D8\nt/t5EfHTwG+nZdPtfGBTYUNEdAPPA/8ZOBt4QxoRfSkibiE75+0tEfGWqe58VgQG0AO8Pb1h3gL8\neRp6QXaa+19ExPlkJ/P9hzLWUbhJ6u2p7UrgBkk/JPsQOo3RU+//NV0fa5DsG+3PlKOo9CG2guyb\n8zfHLF4AfEXSk8DNZC+oYRsjYn85akq13JOm70nzMPFzMgh8rUy1FPscfQV4W/oQ+DXgr8tVXw5f\nT7ebyB7HK8lEr6XxrAKejYjn0vzdk/QtSUQcBi4hu1xRJ/Bl4DeAC4CN6TPhv5GdiDzs7rTu98i2\nQDRPc1linKtkpPbLgS9ExECq4aTf+xVxHkYOAv5U0uXAENm1p5akZc9FxA/T9Ey8mQT8TkQ8cFyj\n9GZO/MOV8xjm+4FPA28mC61hHwceioi3K9tJ/92CZUfKUYik08iGxhdICrJza4Lsg3qi56QnhUg5\nndRzFBFHJW0kuxDmfwTKuXN3gOO/wNWPWd6bbgc59e/bqWorm0leS/dPUNN4lxkqm/Sa/S7wXUlb\ngPcBWyNios2q5f5M2MqYL82SmsiupPFsqb9vtowwfoVsaHlJRLwe2MvoC6S3oN9MvJkeAH4rfQtF\n0uskzUvL1kg6O22n/2WynXblcgfwsYjYMqZ9AaM7eN9Txt9f6B3AXRHxqohYERHLgefIRhOn8jkZ\nq5jn6DbgFuDRMo7GAP4NOE9SnaQFZFc9qBQzWdtEryUmqOlp4NUaPYLxl8tVmKRzJBVeyO/1wDag\nRdkOcSTVSCoc1f9yav8Zsk3Z031xwgeBBkk3pN9TBfw52ej428BvKh3kImlRWucQ0JjnzmdLYCwA\n9kVEv6S3AK+a6YIK3AY8BTyWNmn8JaOh9QPgk8CTZC/yvx33HqZBRHRExOfGWfQ/gf8h6Z/Jvp2d\nCu/ixMf6NeDdnMLnZKxinqOI2AR0A18sR03pzdsbEbuAe4HNwJeAx8vx+05GhdQ22WvphJoi4hjZ\n/oG/l/R9si+X5bpi7HzgTmUHcGwm+18+f0QWcp+S9ATwQ+CnC9Y5IOlfgC+Q7f+cVpGdif124J2S\ndpBdcqkH+EOyz6rngc2ptnen1TYA38qz07uiz/ROL9i9wDlkF8yqIfsDvInsny0BfCMiLkj9PwTM\nj4g/PvXV2lTSZroPRcTbZrqWvCSdSbbJYVVEDJXh/i8C/ioiynHETEkqubbJSJofEYfTfs6/AHZE\nxM0VUNd3yV7/7TNdS7EqfR/G+cCPIuIlJj7U8oLhiYj49Cmpyl4R0rD+T4D/Wqaw+E3g/WTHwVeU\nSq4th1+XtA6oJRt5/OUM1/OyUbEjjMIXbERM+3kCZmZ2cio2MMzMrLLMlp3eZmY2wyomMCQtl/RQ\nOnN7q6QPpPZFkjYqu5DWRkkLU/sqST+Q1Jt2dhfe1wfSqfdbJc3GbbBmZhWnYgKD7OSg34uIc4HL\ngPdJOg/4MPBgupDWg2keYD/ZPo7jdnRLugD4dbJLUFxEdqZu0f/03MzMMhUTGBHxYkQ8lqYPkZ0A\ns4zsLNs7U7c7ya4rQ0Tsi4hHgf4xd3Uu8HBEHE2nwP8j2XHJZmZWgooJjELpLM03kF2baUlEvAhZ\nqACLp1j9SeBySacpu8b7tWSnxZuZWQkq7jwMZZcy/hrZ4bTdo9cYzCcitkn6FLAROAw8Qba5y8zM\nSlBRI4x0PaavkV12d/jqnHs1ej35M4B9U91PRNweERdHxOVk+zp2lKtmM7NXiooJjHQa/+3Atoj4\nTMGi+4F1aXodcF+O+1qcbs8i+38LZbvEsZnZK0XFnLiXrt74T8AWskuYQ3bBrEfILjJ2FtmFs94Z\nEfslLQXagabU/zBwXtqM9U+Q2AUyAAAAR0lEQVRkl6/uJ7usw4On9MGYmb0MVUxgmJlZZauYTVJm\nZlbZHBhmZpaLA8PMzHJxYJiZWS4ODDMzy8WBYWZmuTgwzMwsl/8P6O1Xy8vnthAAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10eb232b0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_null.Y.plot()\n",
    "plt.ylim([0,15000])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate a large number of null distributions\n",
    "\n",
    "- The p-values from distributions conforming to the null hypothesis should form a uniform distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [],
   "source": [
    "ps = []\n",
    "for _ in range(10000):\n",
    "    df_null = generate_null()\n",
    "    C, k_1, k_2 = chow(df_null)\n",
    "    ps.append(f.sf(C, k_1, k_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 491.,  524.,  509.,  500.,  511.,  505.,  479.,  564.,  484.,\n",
       "         540.,  501.,  462.,  465.,  485.,  528.,  485.,  510.,  486.,\n",
       "         491.,  480.]),\n",
       " array([  5.99848496e-05,   5.00553245e-02,   1.00050664e-01,\n",
       "          1.50046004e-01,   2.00041344e-01,   2.50036683e-01,\n",
       "          3.00032023e-01,   3.50027363e-01,   4.00022702e-01,\n",
       "          4.50018042e-01,   5.00013382e-01,   5.50008721e-01,\n",
       "          6.00004061e-01,   6.49999401e-01,   6.99994740e-01,\n",
       "          7.49990080e-01,   7.99985420e-01,   8.49980759e-01,\n",
       "          8.99976099e-01,   9.49971439e-01,   9.99966778e-01]),\n",
       " <a list of 20 Patch objects>)"
      ]
     },
     "execution_count": 235,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAADw9JREFUeJzt3H+MZWddx/H3hy4FFWShO22a3a2D\nYUloSIBmUteQKLCEtMV0+0dLSsQuzcZNsBoUoqz6B/76o2i0SEKA1RK2RKAFxW6wik1/BDVuZUqh\ntlTSpdZ2sht2oe0KaUALX/+4z+KwO9s5s3PvzM6z71cyuec855l7v8/8+Mwzzzn3pKqQJPXrOatd\ngCRpsgx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUufWrXYBABs2bKjp6enVLkOS\n1pR77733m1U1tVi/0yLop6enmZ2dXe0yJGlNSfJfQ/q5dCNJnTPoJalzBr0kdc6gl6TOGfSS1DmD\nXpI6Z9BLUucMeknqnEEvSZ07Ld4ZqzPD9O6/O+XPffT6N4+xEunM4oxekjpn0EtS5wx6SeqcQS9J\nnTPoJalzBr0kdc6gl6TOGfSS1DnfMKXu+UYtnemc0UtS5wx6SeqcQS9JnXONXtKPWM45DfC8xuno\njA96T9RJ6p1LN5LUOYNekjpn0EtS5wx6SercGX8ydi3yBLKkpRgU9EkeBb4NfB94pqpmkrwEuBmY\nBh4F3lJVTyYJ8OfAZcDTwNur6kvjL331GbjSePk7NRlLmdG/vqq+OW9/N3BHVV2fZHfbfw9wKbCl\nffwM8KH2qNOA10hLZ57lrNFvB/a27b3AFfPab6qR/cD6JOcv43UkScswdEZfwD8mKeAjVbUHOK+q\nDgFU1aEk57a+G4HH533uXGs7NP8Jk+wCdgFccMEFpz6CNWq5M2utDP8DUg+GBv1rq+pgC/Pbk/zH\ns/TNAm11QsPoj8UegJmZmROOS5LGY9DSTVUdbI+Hgc8CFwPfOLYk0x4Pt+5zwOZ5n74JODiugiVJ\nS7No0Cf5iSQvPLYNvAl4ANgH7GjddgC3tu19wDUZ2QocPbbEI0laeUOWbs4DPju6apJ1wCeq6h+S\nfBG4JclO4DHgqtb/NkaXVh5gdHnltWOveh7XuiUtV++XdS4a9FX1CPCqBdq/BWxboL2A68ZSnSQN\n5KTv5HxnrNYEf4mXZjW/Xn6vTj/e60aSOueMXkvibE36UWvhvRbO6CWpcwa9JHXOoJekzhn0ktQ5\nT8ZKpylPfGtcDHppggxrnQ5cupGkzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCX\npM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1LnBQZ/krCT3Jflc239p\nknuSPJzk5iRnt/bntf0D7fj0ZEqXJA2xlBn9O4GH5u2/D7ihqrYATwI7W/tO4MmqehlwQ+snSVol\ng4I+ySbgzcBftv0AbwA+07rsBa5o29vbPu34ttZfkrQKhs7o3w/8FvCDtn8O8FRVPdP254CNbXsj\n8DhAO3609ZckrYJFgz7JLwCHq+re+c0LdK0Bx+Y/764ks0lmjxw5MqhYSdLSDZnRvxa4PMmjwKcY\nLdm8H1ifZF3rswk42LbngM0A7fiLgCeOf9Kq2lNVM1U1MzU1taxBSJJObtGgr6rfrqpNVTUNXA3c\nWVW/CNwFXNm67QBubdv72j7t+J1VdcKMXpK0MpZzHf17gHclOcBoDf7G1n4jcE5rfxewe3klSpKW\nY93iXf5fVd0N3N22HwEuXqDPd4GrxlCbJGkMfGesJHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxB\nL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS\n1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOrdo0Cd5fpJ/S/KVJA8m\n+f3W/tIk9yR5OMnNSc5u7c9r+wfa8enJDkGS9GyGzOi/B7yhql4FvBq4JMlW4H3ADVW1BXgS2Nn6\n7wSerKqXATe0fpKkVbJo0NfId9ruc9tHAW8APtPa9wJXtO3tbZ92fFuSjK1iSdKSDFqjT3JWki8D\nh4Hbga8DT1XVM63LHLCxbW8EHgdox48C54yzaEnScIOCvqq+X1WvBjYBFwOvWKhbe1xo9l7HNyTZ\nlWQ2yeyRI0eG1itJWqIlXXVTVU8BdwNbgfVJ1rVDm4CDbXsO2AzQjr8IeGKB59pTVTNVNTM1NXVq\n1UuSFjXkqpupJOvb9o8BbwQeAu4CrmzddgC3tu19bZ92/M6qOmFGL0laGesW78L5wN4kZzH6w3BL\nVX0uyVeBTyX5I+A+4MbW/0bg40kOMJrJXz2BuiVJAy0a9FV1P/CaBdofYbRef3z7d4GrxlKdJGnZ\nfGesJHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn\n0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9\nJHXOoJekzhn0ktQ5g16SOrdo0CfZnOSuJA8leTDJO1v7S5LcnuTh9vji1p4kH0hyIMn9SS6a9CAk\nSSc3ZEb/DPDuqnoFsBW4LsmFwG7gjqraAtzR9gEuBba0j13Ah8ZetSRpsEWDvqoOVdWX2va3gYeA\njcB2YG/rthe4om1vB26qkf3A+iTnj71ySdIgS1qjTzINvAa4Bzivqg7B6I8BcG7rthF4fN6nzbW2\n459rV5LZJLNHjhxZeuWSpEEGB32SFwB/Dfx6Vf33s3VdoK1OaKjaU1UzVTUzNTU1tAxJ0hINCvok\nz2UU8n9VVX/Tmr9xbEmmPR5u7XPA5nmfvgk4OJ5yJUlLNeSqmwA3Ag9V1Z/NO7QP2NG2dwC3zmu/\npl19sxU4emyJR5K08tYN6PNa4JeAf0/y5db2O8D1wC1JdgKPAVe1Y7cBlwEHgKeBa8dasSRpSRYN\n+qr6ZxZedwfYtkD/Aq5bZl2SpDHxnbGS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJek\nzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6Seqc\nQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucWDfokH01yOMkD89pekuT2JA+3xxe39iT5\nQJIDSe5PctEki5ckLW7IjP5jwCXHte0G7qiqLcAdbR/gUmBL+9gFfGg8ZUqSTtWiQV9VXwCeOK55\nO7C3be8FrpjXflON7AfWJzl/XMVKkpbuVNfoz6uqQwDt8dzWvhF4fF6/udYmSVol4z4ZmwXaasGO\nya4ks0lmjxw5MuYyJEnHnGrQf+PYkkx7PNza54DN8/ptAg4u9ARVtaeqZqpqZmpq6hTLkCQt5lSD\nfh+wo23vAG6d135Nu/pmK3D02BKPJGl1rFusQ5JPAq8DNiSZA94LXA/ckmQn8BhwVet+G3AZcAB4\nGrh2AjVLkpZg0aCvqree5NC2BfoWcN1yi5IkjY/vjJWkzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0md\nM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmD\nXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnJhL0SS5J8rUkB5Ls\nnsRrSJKGGXvQJzkL+CBwKXAh8NYkF477dSRJw0xiRn8xcKCqHqmq/wE+BWyfwOtIkgaYRNBvBB6f\ntz/X2iRJq2DdBJ4zC7TVCZ2SXcCutvudJF87xdfbAHzzFD93rXLMZwbHfAbI+5Y15p8a0mkSQT8H\nbJ63vwk4eHynqtoD7FnuiyWZraqZ5T7PWuKYzwyO+cywEmOexNLNF4EtSV6a5GzgamDfBF5HkjTA\n2Gf0VfVMkl8FPg+cBXy0qh4c9+tIkoaZxNINVXUbcNsknnsBy17+WYMc85nBMZ8ZJj7mVJ1wnlSS\n1BFvgSBJnVszQb/YbRWSPC/Jze34PUmmV77K8Row5ncl+WqS+5PckWTQpVans6G3z0hyZZJKsuav\n0Bgy5iRvad/rB5N8YqVrHLcBP9sXJLkryX3t5/uy1ahzXJJ8NMnhJA+c5HiSfKB9Pe5PctFYC6iq\n0/6D0UndrwM/DZwNfAW48Lg+vwJ8uG1fDdy82nWvwJhfD/x4237HmTDm1u+FwBeA/cDMate9At/n\nLcB9wIvb/rmrXfcKjHkP8I62fSHw6GrXvcwx/xxwEfDASY5fBvw9o/chbQXuGefrr5UZ/ZDbKmwH\n9rbtzwDbkiz05q21YtExV9VdVfV0293P6D0La9nQ22f8IfDHwHdXsrgJGTLmXwY+WFVPAlTV4RWu\ncdyGjLmAn2zbL2KB9+KsJVX1BeCJZ+myHbipRvYD65OcP67XXytBP+S2Cj/sU1XPAEeBc1akuslY\n6q0kdjKaEaxli445yWuAzVX1uZUsbIKGfJ9fDrw8yb8k2Z/kkhWrbjKGjPn3gLclmWN0Bd+vrUxp\nq2ait46ZyOWVEzDktgqDbr2whgweT5K3ATPAz0+0osl71jEneQ5wA/D2lSpoBQz5Pq9jtHzzOkb/\ntf1TkldW1VMTrm1Shoz5rcDHqupPk/ws8PE25h9MvrxVMdH8Wisz+iG3VfhhnyTrGP2792z/Kp3u\nBt1KIskbgd8FLq+q761QbZOy2JhfCLwSuDvJo4zWMvet8ROyQ3+2b62q/62q/wS+xij416ohY94J\n3AJQVf8KPJ/RfXB6Nej3/VStlaAfcluFfcCOtn0lcGe1sxxr1KJjbssYH2EU8mt93RYWGXNVHa2q\nDVU1XVXTjM5LXF5Vs6tT7lgM+dn+W0Yn3kmygdFSziMrWuV4DRnzY8A2gCSvYBT0R1a0ypW1D7im\nXX2zFThaVYfG9eRrYummTnJbhSR/AMxW1T7gRkb/3h1gNJO/evUqXr6BY/4T4AXAp9t558eq6vJV\nK3qZBo65KwPH/HngTUm+Cnwf+M2q+tbqVb08A8f8buAvkvwGoyWMt6/liVuSTzJaetvQzju8F3gu\nQFV9mNF5iMuAA8DTwLVjff01/LWTJA2wVpZuJEmnyKCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9J\nnTPoJalz/wdTAuwEFD8KggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e4fba90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(ps, bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like we would expect"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Show how the distribution of p-values changes with a sample that has a break point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 238,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_change = []\n",
    "for _ in range(1000):\n",
    "    df = generate_change()\n",
    "    C, k_1, k_2 = chow(df)\n",
    "    p_change.append(f.sf(C, k_1, k_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 628.,  215.,   73.,   37.,   18.,   12.,    7.,    3.,    0.,\n",
       "           5.,    0.,    1.,    0.,    0.,    0.,    0.,    0.,    0.,\n",
       "           0.,    1.]),\n",
       " array([  1.12148413e-05,   2.06047951e-02,   4.11983753e-02,\n",
       "          6.17919555e-02,   8.23855358e-02,   1.02979116e-01,\n",
       "          1.23572696e-01,   1.44166276e-01,   1.64759857e-01,\n",
       "          1.85353437e-01,   2.05947017e-01,   2.26540597e-01,\n",
       "          2.47134178e-01,   2.67727758e-01,   2.88321338e-01,\n",
       "          3.08914918e-01,   3.29508499e-01,   3.50102079e-01,\n",
       "          3.70695659e-01,   3.91289239e-01,   4.11882820e-01]),\n",
       " <a list of 20 Patch objects>)"
      ]
     },
     "execution_count": 239,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAD49JREFUeJzt3W2MZmV9x/HvT1Zs69PyMBCyu3a0\nblpNU4VM6TYmprq24aFheQEJTVs2ZJNNW2xtMKnbh6TpwwusSVESQ92I7dJogdIaNkptyQIxvoA6\nKKKwGkZC2clSdiyw1hI11H9fzLU6XWb3PrPzfO33k0zOOdf53/f9vy+H35w997mPqSokSf16xWo3\nIElaXga9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMbVrsBgHPPPbfGx8dXuw1J\nWlcefvjhb1XV2Ki6NRH04+PjTE5OrnYbkrSuJPmPIXWeupGkzhn0ktQ5g16SOmfQS1LnDHpJ6pxB\nL0mdM+glqXMGvSR1zqCXpM6tiW/GLsb4ns8u6vFP3Xj5EnUiSWuTR/SS1DmDXpI6Z9BLUucMeknq\n3KCgT7IxyV1Jvp7kYJJfTHJ2knuTPNGWZ7XaJLk5yVSSR5NctLxvQZJ0MkOP6D8CfK6qfgZ4G3AQ\n2AMcqKqtwIG2DXApsLX97AZuWdKOJUkLMjLok7wOeCdwK0BVfb+qXgB2APta2T7gyra+A7itZj0I\nbExywZJ3LkkaZMgR/ZuAGeBvk3w5yceTvBo4v6qeAWjL81r9JuDQnMdPtzFJ0ioYEvQbgIuAW6rq\nQuB/+NFpmvlknrF6WVGyO8lkksmZmZlBzUqSFm5I0E8D01X1UNu+i9ngf/bYKZm2PDKnfsucx28G\nDh//pFW1t6omqmpibGzk/7etJOkUjQz6qvpP4FCSn25D24HHgf3Azja2E7i7re8Hrm1X32wDjh47\nxSNJWnlD73Xzu8Ank5wJPAlcx+wfiTuT7AKeBq5utfcAlwFTwIutVpK0SgYFfVU9AkzMs2v7PLUF\nXL/IviRJS8RvxkpS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEv\nSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLU\nOYNekjo3KOiTPJXkq0keSTLZxs5Ocm+SJ9ryrDaeJDcnmUryaJKLlvMNSJJObiFH9O+qqrdX1UTb\n3gMcqKqtwIG2DXApsLX97AZuWapmJUkLt5hTNzuAfW19H3DlnPHbataDwMYkFyzidSRJizA06Av4\ntyQPJ9ndxs6vqmcA2vK8Nr4JODTnsdNtTJK0CjYMrHtHVR1Och5wb5Kvn6Q284zVy4pm/2DsBnjD\nG94wsA1J0kINOqKvqsNteQT4NHAx8OyxUzJteaSVTwNb5jx8M3B4nufcW1UTVTUxNjZ26u9AknRS\nI4M+yauTvPbYOvArwNeA/cDOVrYTuLut7weubVffbAOOHjvFI0laeUNO3ZwPfDrJsfpPVdXnknwR\nuDPJLuBp4OpWfw9wGTAFvAhct+RdS5IGGxn0VfUk8LZ5xv8L2D7PeAHXL0l3kqRF85uxktQ5g16S\nOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalz\nBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzg0O+iRnJPlyks+0\n7TcmeSjJE0nuSHJmG39V255q+8eXp3VJ0hALOaJ/H3BwzvYHgZuqaivwPLCrje8Cnq+qNwM3tTpJ\n0ioZFPRJNgOXAx9v2wHeDdzVSvYBV7b1HW2btn97q5ckrYKhR/QfBv4A+EHbPgd4oapeatvTwKa2\nvgk4BND2H231kqRVMDLok/wqcKSqHp47PE9pDdg393l3J5lMMjkzMzOoWUnSwg05on8HcEWSp4Db\nmT1l82FgY5INrWYzcLitTwNbANr+1wPPHf+kVbW3qiaqamJsbGxRb0KSdGIjg76q/rCqNlfVOHAN\ncF9V/TpwP3BVK9sJ3N3W97dt2v77quplR/SSpJWxmOvoPwDckGSK2XPwt7bxW4Fz2vgNwJ7FtShJ\nWowNo0t+pKoeAB5o608CF89T813g6iXoTZK0BPxmrCR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6Seqc\nQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0\nktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMjgz7JjyX59yRfSfJYkj9r429M8lCSJ5LckeTMNv6q\ntj3V9o8v71uQJJ3MkCP67wHvrqq3AW8HLkmyDfggcFNVbQWeB3a1+l3A81X1ZuCmVidJWiUjg75m\nfadtvrL9FPBu4K42vg+4sq3vaNu0/duTZMk6liQtyKBz9EnOSPIIcAS4F/gm8EJVvdRKpoFNbX0T\ncAig7T8KnLOUTUuShhsU9FX1v1X1dmAzcDHwlvnK2nK+o/c6fiDJ7iSTSSZnZmaG9itJWqAFXXVT\nVS8ADwDbgI1JNrRdm4HDbX0a2ALQ9r8eeG6e59pbVRNVNTE2NnZq3UuSRhpy1c1Yko1t/ceB9wAH\ngfuBq1rZTuDutr6/bdP231dVLzuilyStjA2jS7gA2JfkDGb/MNxZVZ9J8jhwe5K/BL4M3NrqbwX+\nPskUs0fy1yxD35KkgUYGfVU9Clw4z/iTzJ6vP378u8DVS9KdJGnR/GasJHXOoJekzhn0ktQ5g16S\nOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalz\nBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ3bsNoNrLbxPZ895cc+dePlS9iJJC2PkUf0\nSbYkuT/JwSSPJXlfGz87yb1JnmjLs9p4ktycZCrJo0kuWu43IUk6sSGnbl4C3l9VbwG2AdcneSuw\nBzhQVVuBA20b4FJga/vZDdyy5F1LkgYbGfRV9UxVfamt/zdwENgE7AD2tbJ9wJVtfQdwW816ENiY\n5IIl71ySNMiCPoxNMg5cCDwEnF9Vz8DsHwPgvFa2CTg052HTbUyStAoGB32S1wD/BPx+VX37ZKXz\njNU8z7c7yWSSyZmZmaFtSJIWaFDQJ3klsyH/yar65zb87LFTMm15pI1PA1vmPHwzcPj456yqvVU1\nUVUTY2Njp9q/JGmEIVfdBLgVOFhVfz1n135gZ1vfCdw9Z/zadvXNNuDosVM8kqSVN+Q6+ncAvwl8\nNckjbeyPgBuBO5PsAp4Grm777gEuA6aAF4HrlrRjSdKCjAz6qvoC8593B9g+T30B1y+yL0nSEvEW\nCJLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEv\nSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM6NDPok\nn0hyJMnX5oydneTeJE+05VltPEluTjKV5NEkFy1n85Kk0YYc0f8dcMlxY3uAA1W1FTjQtgEuBba2\nn93ALUvTpiTpVI0M+qr6PPDcccM7gH1tfR9w5Zzx22rWg8DGJBcsVbOSpIU71XP051fVMwBteV4b\n3wQcmlM33cYkSatkqT+MzTxjNW9hsjvJZJLJmZmZJW5DknTMqQb9s8dOybTlkTY+DWyZU7cZODzf\nE1TV3qqaqKqJsbGxU2xDkjTKqQb9fmBnW98J3D1n/Np29c024OixUzySpNWxYVRBkn8Afgk4N8k0\n8KfAjcCdSXYBTwNXt/J7gMuAKeBF4Lpl6FmStAAjg76qfu0Eu7bPU1vA9YttSpK0dPxmrCR1zqCX\npM4Z9JLUOYNekjo38sNYndj4ns+e8mOfuvHyJexEkk7MI3pJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z\n9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalz3utmlXifHEkrxSN6SeqcQS9JnTPoJalzBr0kdc6g\nl6TOedXNOuQVO5IWwiN6SercshzRJ7kE+AhwBvDxqrpxOV5HC7eYfw2A/yKQ1qMlD/okZwAfBX4Z\nmAa+mGR/VT2+1K+lledpI2n9WY4j+ouBqap6EiDJ7cAOwKDXqvCPk053yxH0m4BDc7angV9YhtfR\nOnM6Bu7p+J5PN+vhdOhyBH3mGauXFSW7gd1t8ztJvnGKr3cu8K1TfOzpZF3PUz64Yi/1/+ZpBV/3\nZVbztQdY179PK2jkPC3yf+efHFK0HEE/DWyZs70ZOHx8UVXtBfYu9sWSTFbVxGKfp3fO0zDO0zDO\n0zBrZZ6W4/LKLwJbk7wxyZnANcD+ZXgdSdIAS35EX1UvJXkv8K/MXl75iap6bKlfR5I0zLJcR19V\n9wD3LMdzz2PRp39OE87TMM7TMM7TMGtinlL1ss9JJUkd8RYIktS5dRP0SS5J8o0kU0n2zLP/VUnu\naPsfSjK+8l2uvgHz9M4kX0ryUpKrVqPHtWDAPN2Q5PEkjyY5kGTQZWy9GTBPv5Xkq0keSfKFJG9d\njT5X26h5mlN3VZJKsrJX4lTVmv9h9kPdbwJvAs4EvgK89bia3wH+pq1fA9yx2n2v0XkaB34OuA24\narV7XsPz9C7gJ9r6b/v7dMJ5et2c9SuAz61232txnlrda4HPAw8CEyvZ43o5ov/hbRWq6vvAsdsq\nzLUD2NfW7wK2J5nvy1s9GzlPVfVUVT0K/GA1GlwjhszT/VX1Ytt8kNnvg5xuhszTt+dsvpp5vhx5\nGhiSTwB/AfwV8N2VbA7Wz6mb+W6rsOlENVX1EnAUOGdFuls7hsyTFj5Pu4B/WdaO1qZB85Tk+iTf\nZDbEfm+FeltLRs5TkguBLVX1mZVs7Jj1EvRDbqsw6NYLnXMOhhk8T0l+A5gAPrSsHa1Ng+apqj5a\nVT8FfAD4k2Xvau056TwleQVwE/D+FevoOOsl6IfcVuGHNUk2AK8HnluR7taOQbef0LB5SvIe4I+B\nK6rqeyvU21qy0N+n24Erl7WjtWnUPL0W+FnggSRPAduA/Sv5gex6Cfoht1XYD+xs61cB91X7BOQ0\n4u0nhhk5T+2f2h9jNuSPrEKPa8GQedo6Z/Ny4IkV7G+tOOk8VdXRqjq3qsarapzZz3yuqKrJlWpw\nXQR9O+d+7LYKB4E7q+qxJH+e5IpWditwTpIp4AbghJc49WrIPCX5+STTwNXAx5KcdrenGPj79CHg\nNcA/tksHT7s/mAPn6b1JHkvyCLP/3e08wdN1a+A8rSq/GStJnVsXR/SSpFNn0EtS5wx6SeqcQS9J\nnTPoJalzBr0kdc6gl6TOGfSS1Ln/Az17BxV43AlsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e31e3c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(p_change, bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}