notebooks/Module 2.2.2 - ARIMA.ipynb · sinomiko/QuantFinance

Create your Gitee Account

Explore and code with more than 12 million developers，Free private repositories ！：）

文件

Clone or Download

Module 2.2.2 - ARIMA.ipynb 99.73 KB

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Quantitative Finance\n",
    "\n",
    "Copyright (c) 2019 Python Charmers Pty Ltd, Australia, <https://pythoncharmers.com>. All rights reserved.\n",
    "\n",
    "<img src=\"img/python_charmers_logo.png\" width=\"300\" alt=\"Python Charmers Logo\">\n",
    "\n",
    "Published under the Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC 4.0) license. See `LICENSE.md` for details.\n",
    "\n",
    "Sponsored by Tibra Global Services, <https://tibra.com>\n",
    "\n",
    "<img src=\"img/tibra_logo.png\" width=\"300\" alt=\"Tibra Logo\">\n",
    "\n",
    "\n",
    "## Module 2.2: Modelling Techniques\n",
    "\n",
    "### 2.2.2 ARIMA\n",
    "\n",
    "In this module we will introduce ARIMA, an improvement over the ARMA model we reviewed in a past module (1.6.4). We will quickly review the key parts of ARMA as we build to the new ARIMA model:\n",
    "\n",
    "\n",
    "### Autoregressive Model (AR)\n",
    "\n",
    "The autoregressive model is used for predicting the value of a variable in a time series. We use the annotation $AR(p)$ for an autoregressive model with $p$ periods.\n",
    "\n",
    "$AR(p) X_t = c + \\sum_{i=1}^p{\\beta_i X_{t-i}} + u_t$\n",
    "\n",
    "We can simplify in the case of an AR(1) model, that is $p=1$. This simplifies further if we also assume a zero mean (which can be done by demeaning the data beforehand, giving $c=0$) and an error term $u$ that is white noise:\n",
    "\n",
    "$AR(1) = \\beta X_{t-1}$\n",
    "\n",
    "### Moving Average (MA)\n",
    "\n",
    "A Moving Average (MA) model is given as:\n",
    "\n",
    "$MA(q) X_t = \\mu + \\epsilon_{t} + \\sum_{i=1}^{q}\\theta_i\\epsilon_{t-i}$\n",
    "\n",
    "and specifically the MA(1) process as:\n",
    "\n",
    "$MA(1)X_t = \\epsilon_t + \\theta \\epsilon_{t-1}$\n",
    "\n",
    "Here, the values $\\epsilon_t$ are the error terms for a given time step $t$ and $\\mu$ is the average of the values of $X$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ARMA model\n",
    "\n",
    "The ARMA model is a combination of both the AR model and the MA. It is quite a simple combination - we just concatenate the models, but the training of the model does become more complicated. \n",
    "\n",
    "An $ARMA(p, q)$ model, where $p$ is the lag in the autoregressive model and $q$ is the lag in the moving-average model, given as:\n",
    "\n",
    "$X_t = c + \\epsilon_t + \\sum_{i=1}^{p}{\\beta_i X_{t-i}} + \\sum_{i=1}^{q}\\theta_i\\epsilon_{t-i}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extended Exercise\n",
    "\n",
    "It is often argued that cryptocurrency prices for non-Bitcoin coins \"follow\" Bitcoin, i.e. they lag. Test this hypothesis on a daily level by checking if the values for Ethereum (ETH) and Ripple (XRP) can be modelled using a ARMA model. Test different values for p and q and find the best model.\n",
    "\n",
    "You can get daily prices from Quandl, via BITFINEX: https://www.quandl.com/data/BITFINEX-Bitfinex\n",
    "\n",
    "You can review module 1.6.4 for code on how to run the ARMA model in statsmodels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\statespace\\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.\n",
      "  warn('Non-invertible starting MA parameters found.'\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\statespace\\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.\n",
      "  warn('Non-stationary starting autoregressive parameters'\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\statespace\\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.\n",
      "  warn('Non-invertible starting MA parameters found.'\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XRP:  {'bic':              0            1            2\n",
      "0  1299.589528 -1630.192257 -3575.422884\n",
      "1 -7854.080210 -7893.240446 -7892.451647\n",
      "2 -7888.674767 -7900.853077 -7886.244724\n",
      "3 -7884.022783 -7891.970767 -7898.979903\n",
      "4 -7935.348262 -7916.563303 -7927.404542, 'bic_min_order': (4, 0)}\n",
      "ETH: {'bic':               0             1             2\n",
      "0  47430.355561  43798.197288  40825.073345\n",
      "1  31475.917478  31482.202234  31489.359942\n",
      "2  31482.144164  31479.759838  31486.869773\n",
      "3  31489.322247  31486.929108  31490.332885\n",
      "4  31497.149284  31494.754230  31489.856497, 'bic_min_order': (1, 0)}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:471: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  self._init_dates(dates, freq)\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:471: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  self._init_dates(dates, freq)\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:471: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  self._init_dates(dates, freq)\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:471: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  self._init_dates(dates, freq)\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:471: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  self._init_dates(dates, freq)\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:471: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  self._init_dates(dates, freq)\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\statespace\\sarimax.py:966: UserWarning: Non-stationary starting autoregressive parameters found. Using zeros as starting parameters.\n",
      "  warn('Non-stationary starting autoregressive parameters'\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\tsa\\statespace\\sarimax.py:978: UserWarning: Non-invertible starting MA parameters found. Using zeros as starting parameters.\n",
      "  warn('Non-invertible starting MA parameters found.'\n",
      "c:\\Python311\\Lib\\site-packages\\statsmodels\\base\\model.py:604: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  warnings.warn(\"Maximum Likelihood optimization failed to \"\n"
     ]
    }
   ],
   "source": [
    "#From solutions\n",
    "%run setup.ipy\n",
    "\n",
    "import quandl\n",
    "import my_secrets\n",
    "quandl.ApiConfig.api_key = my_secrets.QUANDL_API_KEY\n",
    "\n",
    "xrp = quandl.get(\"BITFINEX/XRPUSD\")[\"Last\"]\n",
    "eth = quandl.get(\"BITFINEX/ETHUSD\")[\"Last\"]\n",
    "\n",
    "# Should we replace the index? Some dates are missing so we can't set a frquency to 'D' without actually\n",
    "# replacing the entire index? This generates a lot of warnings...\n",
    "\n",
    "\n",
    "# When we've played with values a bit, we can then use the below.\n",
    "from statsmodels.tsa import stattools\n",
    "xrp_stats = stattools.arma_order_select_ic(xrp) \n",
    "eth_stats = stattools.arma_order_select_ic(eth) \n",
    "print(\"XRP: \", xrp_stats)\n",
    "print(\"ETH:\", eth_stats)\n",
    "\n",
    "from statsmodels import api as sms\n",
    "xrp_model = sms.tsa.ARIMA(xrp, order=(4, 0, 2))\n",
    "eth_model = sms.tsa.ARIMA(eth, order=(3, 0, 2))\n",
    "\n",
    "xrp_results = xrp_model.fit()\n",
    "eth_results = eth_model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>SARIMAX Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>Last</td>       <th>  No. Observations:  </th>   <td>2387</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>            <td>ARIMA(4, 0, 2)</td>  <th>  Log Likelihood     </th> <td>3994.813</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Fri, 09 Feb 2024</td> <th>  AIC                </th> <td>-7973.627</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>00:49:18</td>     <th>  BIC                </th> <td>-7927.405</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Sample:</th>                  <td>0</td>        <th>  HQIC               </th> <td>-7956.807</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                      <td> - 2387</td>     <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>        <td>opg</td>       <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "     <td></td>       <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>  <td>    0.4986</td> <td>    0.132</td> <td>    3.766</td> <td> 0.000</td> <td>    0.239</td> <td>    0.758</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L1</th>  <td>    1.2222</td> <td>    0.044</td> <td>   27.640</td> <td> 0.000</td> <td>    1.135</td> <td>    1.309</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L2</th>  <td>   -0.6241</td> <td>    0.066</td> <td>   -9.429</td> <td> 0.000</td> <td>   -0.754</td> <td>   -0.494</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L3</th>  <td>    0.5904</td> <td>    0.047</td> <td>   12.681</td> <td> 0.000</td> <td>    0.499</td> <td>    0.682</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L4</th>  <td>   -0.2025</td> <td>    0.006</td> <td>  -35.731</td> <td> 0.000</td> <td>   -0.214</td> <td>   -0.191</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L1</th>  <td>   -0.0934</td> <td>    0.048</td> <td>   -1.961</td> <td> 0.050</td> <td>   -0.187</td> <td>-5.03e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L2</th>  <td>    0.3363</td> <td>    0.046</td> <td>    7.313</td> <td> 0.000</td> <td>    0.246</td> <td>    0.426</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sigma2</th> <td>    0.0021</td> <td> 1.43e-05</td> <td>  143.853</td> <td> 0.000</td> <td>    0.002</td> <td>    0.002</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Ljung-Box (L1) (Q):</th>     <td>0.00</td> <th>  Jarque-Bera (JB):  </th> <td>418606.36</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Q):</th>                <td>0.99</td> <th>  Prob(JB):          </th>   <td>0.00</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Heteroskedasticity (H):</th> <td>0.21</td> <th>  Skew:              </th>   <td>2.83</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(H) (two-sided):</th>    <td>0.00</td> <th>  Kurtosis:          </th>   <td>67.63</td>  \n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using the outer product of gradients (complex-step)."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                               SARIMAX Results                                \n",
       "==============================================================================\n",
       "Dep. Variable:                   Last   No. Observations:                 2387\n",
       "Model:                 ARIMA(4, 0, 2)   Log Likelihood                3994.813\n",
       "Date:                Fri, 09 Feb 2024   AIC                          -7973.627\n",
       "Time:                        00:49:18   BIC                          -7927.405\n",
       "Sample:                             0   HQIC                         -7956.807\n",
       "                               - 2387                                         \n",
       "Covariance Type:                  opg                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const          0.4986      0.132      3.766      0.000       0.239       0.758\n",
       "ar.L1          1.2222      0.044     27.640      0.000       1.135       1.309\n",
       "ar.L2         -0.6241      0.066     -9.429      0.000      -0.754      -0.494\n",
       "ar.L3          0.5904      0.047     12.681      0.000       0.499       0.682\n",
       "ar.L4         -0.2025      0.006    -35.731      0.000      -0.214      -0.191\n",
       "ma.L1         -0.0934      0.048     -1.961      0.050      -0.187   -5.03e-05\n",
       "ma.L2          0.3363      0.046      7.313      0.000       0.246       0.426\n",
       "sigma2         0.0021   1.43e-05    143.853      0.000       0.002       0.002\n",
       "===================================================================================\n",
       "Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):            418606.36\n",
       "Prob(Q):                              0.99   Prob(JB):                         0.00\n",
       "Heteroskedasticity (H):               0.21   Skew:                             2.83\n",
       "Prob(H) (two-sided):                  0.00   Kurtosis:                        67.63\n",
       "===================================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n",
       "\"\"\""
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xrp_results.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>SARIMAX Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>Last</td>       <th>  No. Observations:  </th>    <td>2814</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>            <td>ARIMA(3, 0, 2)</td>  <th>  Log Likelihood     </th> <td>-15717.368</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Fri, 09 Feb 2024</td> <th>  AIC                </th>  <td>31448.736</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>00:49:48</td>     <th>  BIC                </th>  <td>31490.333</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Sample:</th>                  <td>0</td>        <th>  HQIC               </th>  <td>31463.747</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                      <td> - 2814</td>     <th>                     </th>      <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>        <td>opg</td>       <th>                     </th>      <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "     <td></td>       <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>  <td> 1009.0201</td> <td> 1315.347</td> <td>    0.767</td> <td> 0.443</td> <td>-1569.012</td> <td> 3587.052</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L1</th>  <td>   -0.4165</td> <td>    0.011</td> <td>  -38.570</td> <td> 0.000</td> <td>   -0.438</td> <td>   -0.395</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L2</th>  <td>    0.4609</td> <td>    0.008</td> <td>   57.773</td> <td> 0.000</td> <td>    0.445</td> <td>    0.477</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L3</th>  <td>    0.9512</td> <td>    0.010</td> <td>   91.210</td> <td> 0.000</td> <td>    0.931</td> <td>    0.972</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L1</th>  <td>    1.4376</td> <td>    0.009</td> <td>  154.759</td> <td> 0.000</td> <td>    1.419</td> <td>    1.456</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L2</th>  <td>    0.9717</td> <td>    0.009</td> <td>  107.671</td> <td> 0.000</td> <td>    0.954</td> <td>    0.989</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sigma2</th> <td> 4109.7234</td> <td>   36.911</td> <td>  111.343</td> <td> 0.000</td> <td> 4037.380</td> <td> 4182.067</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Ljung-Box (L1) (Q):</th>     <td>5.08</td>  <th>  Jarque-Bera (JB):  </th> <td>34153.62</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Q):</th>                <td>0.02</td>  <th>  Prob(JB):          </th>   <td>0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Heteroskedasticity (H):</th> <td>11.55</td> <th>  Skew:              </th>   <td>-0.67</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(H) (two-sided):</th>    <td>0.00</td>  <th>  Kurtosis:          </th>   <td>20.02</td> \n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using the outer product of gradients (complex-step)."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                               SARIMAX Results                                \n",
       "==============================================================================\n",
       "Dep. Variable:                   Last   No. Observations:                 2814\n",
       "Model:                 ARIMA(3, 0, 2)   Log Likelihood              -15717.368\n",
       "Date:                Fri, 09 Feb 2024   AIC                          31448.736\n",
       "Time:                        00:49:48   BIC                          31490.333\n",
       "Sample:                             0   HQIC                         31463.747\n",
       "                               - 2814                                         \n",
       "Covariance Type:                  opg                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const       1009.0201   1315.347      0.767      0.443   -1569.012    3587.052\n",
       "ar.L1         -0.4165      0.011    -38.570      0.000      -0.438      -0.395\n",
       "ar.L2          0.4609      0.008     57.773      0.000       0.445       0.477\n",
       "ar.L3          0.9512      0.010     91.210      0.000       0.931       0.972\n",
       "ma.L1          1.4376      0.009    154.759      0.000       1.419       1.456\n",
       "ma.L2          0.9717      0.009    107.671      0.000       0.954       0.989\n",
       "sigma2      4109.7234     36.911    111.343      0.000    4037.380    4182.067\n",
       "===================================================================================\n",
       "Ljung-Box (L1) (Q):                   5.08   Jarque-Bera (JB):             34153.62\n",
       "Prob(Q):                              0.02   Prob(JB):                         0.00\n",
       "Heteroskedasticity (H):              11.55   Skew:                            -0.67\n",
       "Prob(H) (two-sided):                  0.00   Kurtosis:                        20.02\n",
       "===================================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n",
       "\"\"\""
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eth_results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/arma_cryptocurrency.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The ARIMA model\n",
    "\n",
    "One of the issues with the ARMA model is that it requires data to be stationary before the algorithm begins. A normal step for turning non-stationary data into stationary data is to difference the data, either once or twice. Therefore, a normal process is to difference the data and then run ARMA. A problem with this is that we have a two step methodology to describe. ARIMA captures this idea in a single model.\n",
    "\n",
    "The ARIMA model is a more general form of the ARMA model. Specifically, ARMA(p, q) is an ARIMA(p, 0, q) model. The $p$ and $q$ values are the same as the ARMA model. The middle parameter is $d$, for differencing. The ARIMA model will apply differencing to turn a non-stationary dataset into a stationary one, allowing the AR and MA processes to model them better. That's the same as our two step process above, just captured nicely. And a lot more useful.\n",
    "\n",
    "That gives the full ARIMA model as:\n",
    "\n",
    "$ARIMA(p, d, q)$\n",
    "\n",
    "where $p$ is the lag in the autoregressive model and $q$ is the lag in the moving-average model, and $d$ is the order of differencing applied (i.e. how many times to difference the data to get a stationary series).\n",
    "\n",
    "The ARIMA model can also be identified by having components \"turned off\", i.e. set to zero. We saw above that ARIMA(p, 0, q) is simply the ARMA(p, q) model. Further, ARIMA(0, 0, q) is simply MA(q). ARIMA(0, 1, 1) is IMA(1, 1), although that is more obscure than other variants you'll see."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "1. What can the ARIMA(1, 0, 0) model also be known as?\n",
    "2. What does an ARIMA(0, 0, 0) model actually model? Hint: use the ARMA equation, and set p and q to zero.\n",
    "3. What does an ARIMA(0, 1, 0) model represent - i.e. what other term have we seen that refers to the same idea?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. An AR(1) model\n",
    "2. Constant mean plus white noise\n",
    "3. Difference model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/arima_types.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extended Exercise\n",
    "\n",
    "The ARIMA model is implemented in statsmodels under `statsmodels.tsa.arima_model.ARIMA` with a similar use case to the `ARMA` model previously used. Perform an ARIMA modelling on the cryptocurrency data from the previous exercise.\n",
    "\n",
    "Normally, the value for $d$ is determined before running the model, but performing a test of stationarity. See Module 1.6.2 for information on performing these tests. Simply difference the data, check for stationarity, and if it isn't, difference it again. Values more than 3 are abnormal - if you still aren't getting stationary data at that point, check your assumptions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\stattools.py:1685: FutureWarning: The behavior of using lags=None will change in the next release. Currently lags=None is the same as lags='legacy', and so a sample-size lag length is used. After the next release, the default will change to be the same as lags='auto' which uses an automatic lag length selection method. To silence this warning, either use 'auto' or 'legacy'\n",
      "  warn(msg, FutureWarning)\n",
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\stattools.py:1708: InterpolationWarning: p-value is smaller than the indicated p-value\n",
      "  warn(\"p-value is smaller than the indicated p-value\", InterpolationWarning)\n",
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\stattools.py:1708: InterpolationWarning: p-value is smaller than the indicated p-value\n",
      "  warn(\"p-value is smaller than the indicated p-value\", InterpolationWarning)\n",
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:219: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  ' ignored when e.g. forecasting.', ValueWarning)\n",
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:219: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  ' ignored when e.g. forecasting.', ValueWarning)\n",
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:219: ValueWarning: A date index has been provided, but it has no associated frequency information and so will be ignored when e.g. forecasting.\n",
      "  ' ignored when e.g. forecasting.', ValueWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XRP p_value: 0.01\n",
      "ETH p_value: 0.01\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ARMA Model Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>Last</td>       <th>  No. Observations:  </th>    <td>2126</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>            <td>ARMA(3, 1)</td>    <th>  Log Likelihood     </th> <td>-11896.681</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>css-mle</td>     <th>  S.D. of innovations</th>   <td>65.071</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Sat, 19 Feb 2022</td> <th>  AIC                </th>  <td>23805.363</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>14:48:00</td>     <th>  BIC                </th>  <td>23839.335</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Sample:</th>                <td>0</td>        <th>  HQIC               </th>  <td>23817.798</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                       <td> </td>        <th>                     </th>      <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>      <td>  980.3539</td> <td>  770.587</td> <td>    1.272</td> <td> 0.203</td> <td> -529.970</td> <td> 2490.678</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L1.Last</th> <td>    0.1663</td> <td>    0.050</td> <td>    3.296</td> <td> 0.001</td> <td>    0.067</td> <td>    0.265</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L2.Last</th> <td>    0.8629</td> <td>    0.040</td> <td>   21.468</td> <td> 0.000</td> <td>    0.784</td> <td>    0.942</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L3.Last</th> <td>   -0.0318</td> <td>    0.023</td> <td>   -1.395</td> <td> 0.163</td> <td>   -0.076</td> <td>    0.013</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L1.Last</th> <td>    0.8140</td> <td>    0.046</td> <td>   17.794</td> <td> 0.000</td> <td>    0.724</td> <td>    0.904</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>Roots</caption>\n",
       "<tr>\n",
       "    <td></td>   <th>            Real</th>  <th>         Imaginary</th> <th>         Modulus</th>  <th>        Frequency</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AR.1</th> <td>           1.0014</td> <td>          +0.0000j</td> <td>           1.0014</td> <td>           0.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AR.2</th> <td>          -1.1509</td> <td>          +0.0000j</td> <td>           1.1509</td> <td>           0.5000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AR.3</th> <td>          27.2954</td> <td>          +0.0000j</td> <td>          27.2954</td> <td>           0.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>MA.1</th> <td>          -1.2285</td> <td>          +0.0000j</td> <td>           1.2285</td> <td>           0.5000</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                              ARMA Model Results                              \n",
       "==============================================================================\n",
       "Dep. Variable:                   Last   No. Observations:                 2126\n",
       "Model:                     ARMA(3, 1)   Log Likelihood              -11896.681\n",
       "Method:                       css-mle   S.D. of innovations             65.071\n",
       "Date:                Sat, 19 Feb 2022   AIC                          23805.363\n",
       "Time:                        14:48:00   BIC                          23839.335\n",
       "Sample:                             0   HQIC                         23817.798\n",
       "                                                                              \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const        980.3539    770.587      1.272      0.203    -529.970    2490.678\n",
       "ar.L1.Last     0.1663      0.050      3.296      0.001       0.067       0.265\n",
       "ar.L2.Last     0.8629      0.040     21.468      0.000       0.784       0.942\n",
       "ar.L3.Last    -0.0318      0.023     -1.395      0.163      -0.076       0.013\n",
       "ma.L1.Last     0.8140      0.046     17.794      0.000       0.724       0.904\n",
       "                                    Roots                                    \n",
       "=============================================================================\n",
       "                  Real          Imaginary           Modulus         Frequency\n",
       "-----------------------------------------------------------------------------\n",
       "AR.1            1.0014           +0.0000j            1.0014            0.0000\n",
       "AR.2           -1.1509           +0.0000j            1.1509            0.5000\n",
       "AR.3           27.2954           +0.0000j           27.2954            0.0000\n",
       "MA.1           -1.2285           +0.0000j            1.2285            0.5000\n",
       "-----------------------------------------------------------------------------\n",
       "\"\"\""
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#From solutions\n",
    "# Use KPSS to test for stationarity.\n",
    "from statsmodels.tsa.stattools import kpss\n",
    "\n",
    "kpss_stat, p_value, lags, crit = kpss(xrp)\n",
    "print(\"XRP p_value:\", p_value)\n",
    "\n",
    "kpss_stat, p_value, lags, crit = kpss(eth)\n",
    "print(\"ETH p_value:\", p_value)\n",
    "\n",
    "\n",
    "from statsmodels.tsa.arima_model import ARIMA as arima\n",
    "\n",
    "xrp_model_arima = arima(xrp, order=(4,1, 2))  #Remember if we use d=0 we're back to using the ARMA model!\n",
    "eth_model_arima = sms.tsa.ARMA(eth, order=(3, 1, 2))\n",
    "\n",
    "xrp_results_arima = xrp_model_arima.fit()\n",
    "eth_results_arima = eth_model_arima.fit()\n",
    "\n",
    "xrp_results_arima.summary()\n",
    "eth_results_arima.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/arima_cryptocurrency.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choosing parameters\n",
    "\n",
    "When choosing parameters for the ARIMA model, the normal rule of thumb is \"keep them small\". Robert Nau summarised this as:\n",
    "\n",
    "<i>In most cases either p is zero or q is zero, and p+q is less than or equal to 3, so there aren’t very many terms on the right-hand-side of this equation</i>\n",
    "\n",
    "<small>See https://people.duke.edu/~rnau/Notes_on_nonseasonal_ARIMA_models--Robert_Nau.pdf</small>\n",
    "\n",
    "Galit Shmueli has some further insight into this problem in her video on YouTube: https://www.youtube.com/watch?v=0xHf-SJ9Z9U\n",
    "\n",
    "When performing a grid search (in other words, \"try all combinations of these values for the parameters\") you can keep your search space small by trying 0, 1, 2 and 3 as the only options (and just 0, 1, 2 for $d$). If that model isn't sufficient, **and** you have a good theoretical reason for a different value. For seasonal change, for instance, if you suspect a yearly trend (this June's data will be dependent on last June's), then seasonal ARIMA is needed, which we will cover later in this module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import quandl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels import api as sms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "employment = quandl.get(\"FRED/NROUST\").diff().dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "changes = employment.diff().dropna()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.arima_model import ARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:165: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.\n",
      "  % freq, ValueWarning)\n",
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:165: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.\n",
      "  % freq, ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "model = ARIMA(changes, order=(1, 1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>ARIMA Model Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>      <td>D.Value</td>     <th>  No. Observations:  </th>    <td>329</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>          <td>ARIMA(1, 1, 1)</td>  <th>  Log Likelihood     </th> <td>1018.750</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>css-mle</td>     <th>  S.D. of innovations</th>   <td>0.011</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>          <td>Sat, 19 Feb 2022</td> <th>  AIC                </th> <td>-2029.500</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>              <td>14:49:49</td>     <th>  BIC                </th> <td>-2014.316</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Sample:</th>           <td>10-01-1949</td>    <th>  HQIC               </th> <td>-2023.443</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                 <td>- 10-01-2031</td>   <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "        <td></td>           <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th>         <td> 3.741e-07</td> <td> 5.24e-06</td> <td>    0.071</td> <td> 0.943</td> <td>-9.89e-06</td> <td> 1.06e-05</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L1.D.Value</th> <td>   -0.1977</td> <td>    0.054</td> <td>   -3.661</td> <td> 0.000</td> <td>   -0.304</td> <td>   -0.092</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L1.D.Value</th> <td>   -1.0000</td> <td>    0.008</td> <td> -125.172</td> <td> 0.000</td> <td>   -1.016</td> <td>   -0.984</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<caption>Roots</caption>\n",
       "<tr>\n",
       "    <td></td>   <th>            Real</th>  <th>         Imaginary</th> <th>         Modulus</th>  <th>        Frequency</th>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>AR.1</th> <td>          -5.0574</td> <td>          +0.0000j</td> <td>           5.0574</td> <td>           0.5000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>MA.1</th> <td>           1.0000</td> <td>          +0.0000j</td> <td>           1.0000</td> <td>           0.0000</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                             ARIMA Model Results                              \n",
       "==============================================================================\n",
       "Dep. Variable:                D.Value   No. Observations:                  329\n",
       "Model:                 ARIMA(1, 1, 1)   Log Likelihood                1018.750\n",
       "Method:                       css-mle   S.D. of innovations              0.011\n",
       "Date:                Sat, 19 Feb 2022   AIC                          -2029.500\n",
       "Time:                        14:49:49   BIC                          -2014.316\n",
       "Sample:                    10-01-1949   HQIC                         -2023.443\n",
       "                         - 10-01-2031                                         \n",
       "=================================================================================\n",
       "                    coef    std err          z      P>|z|      [0.025      0.975]\n",
       "---------------------------------------------------------------------------------\n",
       "const          3.741e-07   5.24e-06      0.071      0.943   -9.89e-06    1.06e-05\n",
       "ar.L1.D.Value    -0.1977      0.054     -3.661      0.000      -0.304      -0.092\n",
       "ma.L1.D.Value    -1.0000      0.008   -125.172      0.000      -1.016      -0.984\n",
       "                                    Roots                                    \n",
       "=============================================================================\n",
       "                  Real          Imaginary           Modulus         Frequency\n",
       "-----------------------------------------------------------------------------\n",
       "AR.1           -5.0574           +0.0000j            5.0574            0.5000\n",
       "MA.1            1.0000           +0.0000j            1.0000            0.0000\n",
       "-----------------------------------------------------------------------------\n",
       "\"\"\""
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "Play around with the values for order and test out several types. How does the summary results change? What extra variables are added? How are they added, for instance, when you set $p=2$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: There are no solutions for this exercise - just try different values and examine the result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing your ARIMA model\n",
    "\n",
    "With your ARIMA model, you can test the residuals of the model to confirm they are white noise, or whether there is an additional correlation in the residuals that needs to be modelled. The Ljung-Box tests the following hypotheses:\n",
    "\n",
    "$H_0$: The data are independently distributed\n",
    "\n",
    "$H_A$: The data is not independently distributed, that is that they have a serial correlation.\n",
    "\n",
    "Such a serial correlation indicates that the ARIMA model hasn't done its job well, indicative of a bad choice of parameters for $p$ and $q$, mainly $q$.\n",
    "\n",
    "The test itself is:\n",
    "\n",
    "$Q = n(n + 2)\\sum_{k=1}^{h}\\frac{p^2}{n-k}$\n",
    "\n",
    "Where $n$ is the size of the dataset, $h$ is the number of lags being tested"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.stats.diagnostic import acorr_ljungbox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_statistics, p_values = acorr_ljungbox(results.resid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True, False, False, False, False, False, False, False,\n",
       "       False, False, False, False, False, False, False, False, False,\n",
       "       False, False, False, False, False, False, False, False, False,\n",
       "       False, False, False, False, False, False, False, False, False,\n",
       "       False, False, False, False])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_values > 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "Review carefully these results, compared against the hypothesis stated for the Ljung-Box test. What are the results saying?\n",
    "\n",
    "Hint: for documentation on the function itself, see https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.acorr_ljungbox.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/ljungbox.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prediction with ARIMA\n",
    "\n",
    "ARIMA models can be used for predicting future values. \n",
    "The prediction confidence interval from an ARIMA model will be wider for data with a higher volatility.\n",
    "\n",
    "An ARIMA model is often harder to analyse than a simple Linear Regression model. For this reason, analysis of the fitted parameters (and meta-parameters, internal states of the model) are not often performed to understand \"why\" the model fit the way it did. Compare this to Linear Regression models, where we can interpret the $\\beta$ values to understand why the model fit the way it did."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Value</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>2030-10-01</td>\n",
       "      <td>0.000026</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2031-01-01</td>\n",
       "      <td>0.000039</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2031-04-01</td>\n",
       "      <td>0.000048</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2031-07-01</td>\n",
       "      <td>0.000051</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2031-10-01</td>\n",
       "      <td>0.000042</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Value\n",
       "Date                \n",
       "2030-10-01  0.000026\n",
       "2031-01-01  0.000039\n",
       "2031-04-01  0.000048\n",
       "2031-07-01  0.000051\n",
       "2031-10-01  0.000042"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "changes.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Value</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>2000-01-01</td>\n",
       "      <td>-0.000088</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2000-04-01</td>\n",
       "      <td>-0.000089</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2000-07-01</td>\n",
       "      <td>-0.000085</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2000-10-01</td>\n",
       "      <td>-0.000069</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Value\n",
       "Date                \n",
       "2000-01-01 -0.000088\n",
       "2000-04-01 -0.000089\n",
       "2000-07-01 -0.000085\n",
       "2000-10-01 -0.000069"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "changes[\"2000\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Predict values within the sample\n",
    "from datetime import datetime\n",
    "\n",
    "start_date = datetime(2000, 1, 1)\n",
    "end_date = datetime(2001, 10, 1)\n",
    "\n",
    "y_pred = results.predict(start_date, end_date)\n",
    "y_pred.name = \"Predictions\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run setup.ipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEUCAYAAADdvgZNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhU9d3+8fcnewi7YQlLCAjIEjYFQVxKad1R0artYxe1PrX+ulqXLlqtG7VatbaPttZu2Npa64K7VqugKIgGlR1kC/sS9oSErJ/fHzOhQYFMQiZnlvt1XXMlOTNncgeP95z5zjnna+6OiIjEp5SgA4iISPOpxEVE4phKXEQkjqnERUTimEpcRCSOqcRFROKYSlwSgplNNbM7wt+fbGbLmvk8D5nZTS2bTiR6VOLSqsys2MwqzKzMzLaY2V/MrG1L/g53n+nux0SQ5TIze/sT617l7re3ZB6RaFKJSxDOcfe2wLHAGOCnDe80s7RAUonEIZW4BMbdNwAvA4Vm5mb2bTNbDiwHMLNJZvaRme0ys1lmNrx+XTMbZWYfmFmpmT0OZDW4b4KZrW/wc28ze9rMSsxsu5k9YGaDgYeAE8LvCnaFH7t/WCb88zfMbIWZ7TCz58ysR4P73MyuMrPlZrbTzB40Mwvf19/M3jSz3Wa2LZxRpMWpxCUwZtYbOAv4MLxoMjAWGGJmxwJ/Br4JHAX8HnjOzDLNLAN4Bvgb0Bl4AvjCIX5HKvACsAYoAHoC/3T3JcBVwGx3b+vuHQ+y7kTgTuBiIC/8HP/8xMMmEXo3MSL8uNPDy28HXgU6Ab2A/4v030WkKaJW4mb2ZzPbamYLW+j5asN7ZR+Z2XMt8ZwSmGfCe75vA28CPw8vv9Pdd7h7BfAN4PfuPsfda939EaASGBe+pQP3u3u1uz8JvH+I33U80AO43t33uvs+d3/7EI/9pC8Df3b3D9y9EvgJoT33ggaP+YW773L3tcB0YGR4eTXQB+jRxN8p0iTR3BOfCpzRgs9X4e4jw7dzW/B5pfVNdveO7t7H3b8VLm2AdQ0e0we4NjyUsitc+r0JFXIPYIMfePW2NYf4Xb2BNe5e04ycPRo+r7uXAdsJ7c3X29zg+3Kg/kPaHwIGvGdmi8zs6834/SKNilqJu/tbwI6Gy8zsaDN7xczmmtlMMxsUrd8vcalhKa8DpoTLvv7Wxt0fAzYBPevHn8PyD/Gc64D8Q3xY2tglPDcSejEBwMxyCA3tbGj0D3Hf7O7fcPcehIaEfmtm/RtbT6SpWntM/GHgu+5+HHAd8NsmrJtlZkVm9q6ZTY5OPIkhfwCuMrOxFpJjZmebWTtgNlADfM/M0szsAkLDJgfzHqHS/0X4ObLM7MTwfVuAXuEx9oP5B3C5mY00s0xCwz5z3L24sfBmdpGZ9Qr/uJPQC0Zt43+2SNO02qFc4WOBxwNPNNiBygzfdwFw20FW2+Du9R8U5bv7RjPrB7xhZgvcfWW0c0sw3L3IzL4BPAAMACoIjaG/5e5V4W3mD8AdwEvA04d4nlozOwf4DbCWUJn+A3gHeANYBGw2szp3z/3Euq+HT/x5itAHlLOAL0X4J4wB7jezDoReLL7v7qsj/gcQiZBFc1KI8AdAL7h7oZm1B5a5e14LPO/U8PM+eaTPJSISz1ptOMXd9wCrzewigPBb5BGRrGtmncJvZzGzXOBEYHHUwoqIxIloHmL4GKGxy2PMbL2ZXUHokK0rzGweobex50X4dIOBovB60wkd1qUSF5GkF9XhFBERiS6dsSkiEsdU4iIicSwqhxjm5uZ6QUFBNJ5aRCQhzZ07d5u7d2nqelEp8YKCAoqKiqLx1CIiCcnMDnXpiMPScIqISBxTiYuIxDGVuIhIHFOJi4jEMZW4iEgcU4mLiMQxlXgCqa6to7yqORPYiEi8arXriUvL2rOvmqWbSlm8cTeLN+1hyaZSlm0pBYdJw/O4dHwBI3p/au5fEUkwKvEY5+5s3L2PxRv3sGTTHhZv3MPiTXtYu6N8/2M652QwtEd7Lh9fQHlVLU9/sJ6nP9zAqPyOXDa+gDML88hI05sukUQUlasYjh492nXGZtNV1dSxsqRsf1HXf91dUQ2AGfQ9KofBee0Z0qM9Q8Jfu7bLpOF0k6X7qnly7nr+OnsNq7ftpUu7TL48Np9LxubTtV1WUH+eiByGmc1199FNXk8lHozdFdUH7Fkv3riH5VtLqa4N/ffITEthUN5/i3pIXnsGdW9HTmbkb57q6pw3l5fwyKxiZiwrIT3VOHtYaKhlVH6naP1pItIMKvEY5e6s31mxv6iXbAqV9vqdFfsfk9s2gyE9OjQo7HYUHJVDWmrLDYGsKinjr7PX8OTc9ZRV1jCid0cuH1/AWcM01CISC1TiMaCyppblW8r2F3X9XnbpvtARI2bQNzfngL3r0HBI6w1xlFXW8NTc9Twyu5hVJXvJbZvJJWPz+crYfLq211CLSFBiqsSPPe44L3q/iJQUa/zBcWpXedUBRb144x5WbC2jpi7075mdnsqgvHYHFPYx3dvRJiM2Pkuuq3NmrtjGI7OKmb5sK6lmnDUsj8tOLGBU744HjLGLSPTFVIln5g3wvEvvJzMtheyMVNqkp5KVkUp2eviWceDXrPRU2jRY1vDn+vXaNHhsdkbo56y01Ki/UNTV1Q+H7D6gsDfu3rf/MV3bZe4v6voPHQuOyiE1Tl7Eirft5a+z1/BE0TpKK2sY3qsDl55QwKQReWSmpQYdTyQpxFSJ9zlmmP/gwaeoqK6loip8q65lX3Ut5eHvK6oO/Hlfde3+D/WaItIXiuxP3HeoF4raOufjLaXh8etSlmzaQ2llaDgkxeDoLm0PODpkcF57urTLbOl/wkCUVdYw7YP1TJ1VzMqSveS2zeCS4/P58rg+dNNQi0hUxVSJN3dMvLq2LlTo9UUfLvlP/Rx+Edj/ff3Pn1i3ourTX+uHOxrTJiM1VNYNCntgt3ZkZyT+nqm783Z4qOX1paGhljOH5XHZ+D4cm99JQy0iUdDcEo+NAdqw9NQU0lNTaJ+VHrXf0fCForxB4df/XOfOgG7t6NO5TUKP6R+OmXHygC6cPKALa7bv5W+z1/B40Tqen7eRwp7tuWx8XyYNzyMrPfFf0ERiXUztiUvs2ltZw7QPN/DIrGKWby2jc079UEs+eR2yg44nEvcSYjhFYp+7M2vldv7yTjGvL91CihlnFHbnsvEFjO6joRaR5kqI4RSJfWbGif1zObF/Lmu3l/O3d4t5/P11vDh/E0N7tOfS8QWcO6KHhlpEWklEe+JmVgyUArVATWOvFtoTTy7lVTU88+FGps5azcdbyujUJp3/OT6fr4zrQ4+OGmoRiURUh1PCJT7a3bdF8qQq8eTk7sxetZ2p7xTznyVbMDNOH9qNS08o4Pi+nTXUInIYGk6RwJkZ44/OZfzRuazbUc6jc9bwz/fW8dKCzQzOa89l4/tw3sieGmoRaUGR7omvBnYCDvze3R8+3OO1Jy71KqpqefajDUydVczSzaV0bJPOl8bk89UT+tBTQy0i+0V7OKWHu280s67Aa8B33f2tTzzmSuBKgPz8/OPWrFnT1CySwNydOat3MPWdYl5dvBmA04Z059LxBYzrp6EWkVY7xNDMbgHK3P2eQz1Ge+JyOOt3lvPou2v55/tr2VVezaDu7bh0fAGTR/ZMijNiRQ4maiVuZjlAiruXhr9/DbjN3V851DoqcYnEvur6oZY1LNm0hw7Z6Vxz6kAuHV8QdDSRVhfNDza7AdPCb3fTgH8crsBFIpWVnsoXx+Rz8ejevF+8k/v/8zG3PL+IUfkdGd5LkzyLRKLRKV3cfZW7jwjfhrr7lNYIJsnDzDi+b2ce+upxdGmbyQ3TFlAb4YXKRJKd5uWSmNE+K52bzxnCwg17+Nvs4qDjiMQFlbjElLOH5XHKwC7c8+rHbNmzr/EVRJKcSlxiiplx+3lDqaqt47YXFgcdRyTmqcQl5vQ5KofvfrY/L87fxIxlW4OOIxLTVOISk678TD/6dcnh5mcXsa+6Nug4IjFLJS4xKTMtlTsmF7J2RzkPvLEi6DgiMUslLjFr/NG5XDCqJ79/ayUrtpYGHUckJqnEJabdcPZgstNTuXHaQqIxC5VIvFOJS0zLbZvJj88czJzVO3j6gw1BxxGJOSpxiXlfGtObY/M7MuWlJewqrwo6jkhMUYlLzEtJMaacP4zdFdXc9crSoOOIxBSVuMSFwXntueKkvjz23jrmrtkRdByRmKESl7jx/c8NoEeHLG6ctpDq2rqg44jEBJW4xI2czDRuOXcoSzeX8pd3VgcdRyQmqMQlrpw2tDufH9yNX722nA27KoKOIxI4lbjEnVvOHRL6+tyigJOIBE8lLnGnV6c2XP35Aby2eAuvLtocdByRQKnEJS59/aS+HNOtHbc8t4i9lTVBxxEJjEpc4lJ6agpTzi9k4+59/Ob15UHHEQmMSlzi1uiCznxpTG/++PZqlm7eE3QckUCoxCWu/eiMQXTITueGpxdQp8mVJQmpxCWudcrJ4IazBvPB2l08XrQu6DgirU4lLnHvC8f2ZGzfzvzi5aVsK6sMOo5Iq1KJS9wzM6acX0h5VQ0/f2lJ0HFEWpVKXBJC/67t+OYpR/P0BxuYtXJb0HFEWo1KXBLGdyb2J79zG376zEIqazS5siQHlbgkjKz0VG47byirSvbyh7dWBR1HpFWoxCWhTDimK2cPy+P/3ljBmu17g44jEnUqcUk4N00aQnpqCjc/u0iTK0vCU4lLwuneIYtrTxvImx+X8NICXSBLEptKXBLSV8f1obBne259fhGl+6qDjiMSNSpxSUhpqSlMmTyMkrJK7n3146DjiERNxCVuZqlm9qGZvRDNQCItZUTvjnx1XB/+OruYBet3Bx1HJCqasif+fUCnw0lcue70YziqbSY3PrOAWl0gSxJQRCVuZr2As4E/RjeOSMtqn5XOTZOGMH/9bv4+Z03QcURaXKR74vcDPwTqophFJCrOGZ7HyQNy+eUry9iyZ1/QcURaVKMlbmaTgK3uPreRx11pZkVmVlRSUtJiAUWOlJlx23mFVNbWcfsLi4OOI9KiItkTPxE418yKgX8CE83s0U8+yN0fdvfR7j66S5cuLRxT5Mj0zc3h2xP688L8Tbz5sXYyJHE0WuLu/hN37+XuBcCXgDfc/StRTybSwq6a0I9+uTnc/OxC9lXrAlmSGHScuCSNzLRU7phcyJrt5fx2+oqg44i0iCaVuLvPcPdJ0QojEm3j++dy/qie/O7NlazYWhZ0HJEjpj1xSTo3nDWY7PRUbnpmoS6QJXFPJS5Jp0u7TH505iBmr9rOMx9tCDqOyBFRiUtS+p8x+YzK78gdLyxhd7kukCXxSyUuSSklxZgyeRi7Kqq5699Lg44j0mwqcUlaQ3q05/LxBfxjzlrmrtkZdByRZlGJS1K7+tSB5HXI4sZpC6ip1VUlJP6oxCWptc1M42fnDGXp5lKmzioOOo5Ik6nEJemdPrQbnxvUlfte+5iNuyqCjiPSJCpxSXpmxi3nDqXOnVufXxR0HJEmUYmLAL07t+H7nxvIvxdt4T+LtwQdRyRiKnGRsP89uS8Du7XlZ88toryqJug4IhFRiYuEpaemMOX8YWzYVcGvX18edByRiKjERRoYU9CZL47uzZ9mrmbp5j1BxxFplEpc5BN+fOYg2mWl8dNpC6nT5MoS41TiIp/QKSeDG84aTNGanTwxd13QcUQOSyUuchAXHteL4/t25s6Xl7K9rDLoOCKHpBIXOQgzY8rkQsr21XDny7pAlsQulbjIIQzo1o4rT+nHk3PX8+6q7UHHETkolbjIYXx34gB6dcrmp88spKpGF8iS2KMSFzmM7IxUbj+vkBVby/jDzFVBxxH5FJW4SCM+O6grZxZ25zevL2ft9vKg44gcQCUuEoGbzxlCWopx83OaXFlii0pcJAJ5HbK55rRjmLGshFcWbg46jsh+KnGRCF16Qh+G5LXnlucXUVapC2RJbFCJi0QoLTWFKecXsrW0kvte/TjoOCKASlykSUbld+LLY/OZOms1CzfsDjqOiEpcpKmuP30QnXMyuXHaAmp1gSwJmEpcpIk6ZKdz06TBzFu/m3/MWRN0HElyKnGRZjh3RA9O6p/L3a8sY2vpvqDjSBJTiYs0g5lx++RCKmvruOOFJUHHkSSmEhdppr65OXxrwtE8N28jM5eXBB1HkpRKXOQIXPWZo+mbm8NNzyxkX3Vt0HEkCTVa4maWZWbvmdk8M1tkZre2RjCReJCVHrpAVvH2cn43Y2XQcSQJRbInXglMdPcRwEjgDDMbF91YIvHjpAG5nDeyB7+bsZJVJWVBx5Ek02iJe0j9lpkevungWJEGbjx7MJnpKfz46QWaXFlaVURj4maWamYfAVuB19x9TnRjicSXru2yuHnSEN5bvYM/v7M66DiSRCIqcXevdfeRQC/geDMr/ORjzOxKMysys6KSEn1SL8nnwuN68fnB3bj738tYvqU06DiSJJp0dIq77wJmAGcc5L6H3X20u4/u0qVLC8UTiR9mxp0XDKNtZhrX/Gse1bWazk2iL5KjU7qYWcfw99nA5wFN/y1yEF3aZTJlciELNuzmwekrgo4jSSCSPfE8YLqZzQfeJzQm/kJ0Y4nErzOH5TF5ZA8eeGMFC9brSocSXZEcnTLf3Ue5+3B3L3T321ojmEg8u/XcQnLbZnLNvz7SSUASVTpjUyQKOrRJ564Lh7N8axn3vros6DiSwFTiIlHymYFd+PLYfP749mrmrNoedBxJUCpxkSi64azB9O7UhuuenKd5OSUqVOIiUZSTmca9F49g/c4KpryoS9ZKy1OJi0TZmILOXHlyPx57by3Tl20NOo4kGJW4SCv4wakDGditLT96cj67yquCjiMJRCUu0gqy0lO57+KR7Nhbxc3PLgo6jiQQlbhIKyns2YHvfW4Az83byIvzNwUdRxKESlykFX1rwtGM6NWBnz6zQBMsS4tQiYu0orTUFO69eCTlVbX85KkFuOva43JkVOIirax/17b88IxBvL50K08UrQ86jsQ5lbhIAC4fX8C4fp257YXFrNtRHnQciWMqcZEApKQYv7xwBO7O9U/O05Ru0mwqcZGA9O7chpsmDeHdVTuYOqs46DgSp1TiIgH64pjeTBzUlbteWcqKrWWNryDyCSpxkQCZGb+4YBjZGalc+8Q8ajSlmzSRSlwkYF3bZ3H7eYXMW7eL381YGXQciTMqcZEYcM6IHkwansevX1/Owg2a0k0ipxIXiRG3n1dIp5wMrv3XPCprNKWbREYlLhIjOuVkcPcXhrNsSyn3vfZx0HEkTqjERWLIZwd15UtjevPwW6soKt4RdByJAypxkRjz00lD6Nkxm2ufmMdeTekmjVCJi8SYtplp3HPRCNbuKOfOlzWlmxyeSlwkBo3rdxRfP7Evj767lrc+Lgk6jsQwlbhIjLr+9GNCVzx8cj67y6uDjiMxSiUuEqNCU7qNoKSsklue15RucnAqcZEYNrxXR7792f5M+3ADryzUlG7yaSpxkRj33Yn9KezZnhumLaSktDLoOBJjVOIiMS49NYX7Lh5JWWUNN0zTlG5yIJW4SBwY2K0d1502kNcWb+GpDzYEHUdiiEpcJE5ccVI/ji/ozK3PLWLDroqg40iMUImLxInUFOOei0ZQ684PNaWbhDVa4mbW28ymm9kSM1tkZt9vjWAi8mn5R7XhxrMH886K7fzt3TVBx5EYEMmeeA1wrbsPBsYB3zazIdGNJSKHcsnx+XxmYBfufHkJq0o0pVuya7TE3X2Tu38Q/r4UWAL0jHYwETk4M+OuLwwnIzVFU7pJ08bEzawAGAXMiUYYEYlM9w5Z3D65kA/X7uL3b60KOo4EKOISN7O2wFPA1e6+5yD3X2lmRWZWVFKiC/aIRNu5I3pw1rDu3P+fj1m88VP/S0qSiKjEzSydUIH/3d2fPthj3P1hdx/t7qO7dOnSkhlF5CDMjDsmD6NDdgbX/OsjTemWpCI5OsWAPwFL3P2+6EcSkUh1zsngFxcMY+nmUn79n+VBx5EARLInfiLwVWCimX0Uvp0V5VwiEqHPD+nGRcf14qE3VzJ3zc6g40gri+TolLfd3dx9uLuPDN9eao1wIhKZm88ZQl6HbK57Yh7lVZrSLZnojE2RBNAuK51fXjSc1dv2ctfLS4OOI61IJS6SIMYfnctl4wt4ZPYa3lmxLeg40kpU4iIJ5EdnDKJfbg7XPzGPPfs0pVsyUImLJJDsjFTuvXgEm/fs49bnFgcdR1qBSlwkwYzK78S3JvTnqQ/W8+qizUHHkShTiYskoO99bgBD8tpzw7QFbC/TlG6JTCUukoAy0lK474sj2FNRw43TFmpKtwSmEhdJUIO6t+cHpw7klUWbeeYjTemWqFTiIgnsylP6cVyfTtz87CI27daUbolIJS6SwFJTjHsvGkFNrfPDJ+drWCUBqcRFElxBbg43nDWImcu38eictUHHkRamEhdJAl8Z14eTB+Ty8xeXULxtb9BxpAWpxEWSgJlx94XDSUs1rntiHrV1GlZJFCpxkSSR1yGbW88dStGanfxhpqZ0SxQqcZEkcv6onpw+tBv3vfoxSzdrSrdEoBIXSSJmxs/PH0a7rDSueXweVTV1QUeSI6QSF0kyR7XN5OcXDGPxpj383xua0i3eqcRFktDpQ7tzwbE9+e2MlXy0blfQceQIqMRFktTPzhlK13aZXPOvj9hXXRt0HGkmlbhIkuqQnc4vLxzBqpK93PWKpnSLVypxkSR20oBcvnZCH/7yTjGzVmpKt3ikEhdJcj8+cxAFR7Xh+ifma0q3OKQSF0lybTLSuPfikWzaXcHZv5nJKws360JZcUQlLiIc16cTj14xlqy0VK56dC5f+dMclm0uDTqWREAlLiIAjO+fy8vfP5lbzhnCgvW7OfPXb3HzswvZVV4VdDQ5DJW4iOyXlprCZSf2Zcb1n+WSsfk8+u4aJtwzg7/OLqamVmd3xiKVuIh8SuecDO6YPIwXv3cyg7q34+ZnF3H2b95m1godwRJrVOIickiD89rz2DfG8dBXjmVvVQ2X/HEOV/1tLut2lAcdTcLSgg4gIrHNzDijMI8Jx3TljzNX8eD0lbyxbCvfOLkv35rQn5xM1UiQtCcuIhHJSk/lOxMH8MZ1n+Gswu48OH0lE++dwbQP1+uQxACpxEWkSfI6ZHP/l0bx1P87gW7ts/jB4/P4wu9mMU8X0gqESlxEmuW4Pp155lsncveFw1m7o4LzHnyH65+Yx9bSfUFHSyqNlriZ/dnMtprZwtYIJCLxIyXFuHh0b6Zf9xm+eUo/nvloAxPveZOH3lxJZY2ujNgaItkTnwqcEeUcIhLH2mWl85OzBvPqDz7D2L6d+cXLSzn9V2/xn8VbNF4eZY2WuLu/BexohSwiEuf65ubwp8vGMPXyMaSmGP/71yIu/cv7rNiqU/ijRWPiItLiJhzTlVeuPoWbJg3hw7U7Of3+mdz6/CJ2V+gqiS2txUrczK40syIzKyopKWmppxWROJWemsIVJ/VlxnUTuHh0b6bOKuaz98zg73PWUFunIZaWYpGMV5lZAfCCuxdG8qSjR4/2oqKiI0smIgll4Ybd3Pb8Yt4r3sHgvPbccs4QxvY7KuhYMcPM5rr76Kaup+EUEWkVhT078Pg3x/HAJaPYXV7FFx9+l2//4wPW79Qp/EcikkMMHwNmA8eY2XozuyL6sUQkEZkZk4b34PVrJ3D15wfw+pItfO7eN7nvtY+pqNIhic0R0XBKU2k4RUQisWFXBXe+tIQX5m+iR4csfnzWYM4ZnoeZBR2t1Wk4RUTiTs+O2TxwybH865sn0Ckng+899iEX/342CzfsDjpa3FCJi0jgju/bmee+cxJ3XjCMlSV7OeeBt/nJ0/PZVlYZdLSYpxIXkZiQmmL8z/H5TL9uAl8/sS9PFK3ns7+cwR9nrqKqRrMKHYpKXERiSofsdG6aNIRXrj6FY/t04o4Xl3DGr99i+rKtQUeLSSpxEYlJ/bu25ZGvH8+fLxuNO1z+l/f5+tT3WVVSFnS0mKISF5GYNnFQN/599SnceNZg3lu9g9Pvf4spLy5mzz6dwg8qcRGJAxlpKXzjlH5Mv24C54/qyR/fXs3Ee2bw+PtrqUvyU/hV4iISN7q0y+TuC0fw7LdPpM9ROfzoqQWc9+A7FBUn74VWdbKPiMQld+e5eRu586WlbN6zj7OH51HYowNtMlLJzkilTUYqORlp+78PLU+jTXro/sy0lJg6qai5J/tommoRiUtmxnkje3LqkG48NGMlf5i5mhfnb4p4/RSDNg1KPju9vuzTDiz9/d+n0iY99cB1Gjw+Oz2VnMzQ9635AqESF5G41iYjjWtOO4YfnDqQfdV1lFfVUF5VS0V1LeVVtZRX1VBRFfq+oqqWvfX31y+rrgk/rn5ZDdvKKv+7fmUN5dW1NGXQwozwHv8nXgQyUslO//SyNhnNr2KVuIgkBDMjO1yMLX2BW3ensqbuUy8KB38RqN3/QhJaduCLyo69Ff9dVlVLeXXtEV1fXSUuItIIMyMrPZWs9FQ652S06HO7O1W1dWTd1bz1dXSKiEiAzIzMtNRmr68SFxGJYypxEZE4phIXEYljKnERkTimEhcRiWMqcRGROKYSFxGJY1G5AJaZlQLLWvyJW0cHIJ5naVX+YOUC24IOcQTi/d8/nvMPcPcOTV0pWmdsLmvO1bhigZk97O5XBp2juZQ/WGZWFK/bPiTEv3/c5jezh5uznoZTPu35oAMcIeWXIxHv//7xnL9Z2aM1nBLXeyMizaVtX1pbtPbEm/W2QCQBaNuXVhWVPXEREWkdST0mbma9zOxZM1tuZqvM7AEzyzSzU81srpktCH+dGHTWgzlM/uPN7KPwbZ6ZnR901oM5VP4G9+ebWZmZXRdkzkSkbT9YLbntJ22JW2jupKeBZ9x9ADAAyAbuJnSI2DnuPgy4FPhbYEEPoZH8C4HR7j4SOAP4vZnF1LXjG8lf71fAywHES2ja9oPV4tu+uzf7BvQCngWWA6uAB4BM4ChgOlAGPHAkvyNaN+BzwFufWNYe2Am0bekml8oAAAQYSURBVLDMgO1AZtCZm5m/L7AFSAs6c1PyA5OBXwK3ANcFnfcg+bXtx37+pNj2m70n3siryT7gJiCW3wYPBeY2XODue4BioH+DxV8APnT3ytaLFpHD5jezsWa2CFgAXOXuNa0f8bAOl38E8CPg1taP1Tht+4HTtt/AkQynTAT2uftfwiFqgR8AXyP0genbhDboWGXAwT7V3T9FtZkNBe4CvtlaoZrgsPndfY67DwXGAD8xs6zWDBeBw+W/FfiVu5e1bqSIadsPlrb9Bo6kxCN9NY9Vi4ADjuc1s/ZAN2CZmfUCpgFfc/eVAeRrzGHz1y9z9yXAXqCwVdM17nD5OwB3m1kxcDVwg5l9p9UTHpq2/WBp22/gSEq80VfzGPc60MbMvgZgZqnAvfx3bPNF4Cfu/k5wEQ/rcPm713+YY2Z9gGMIFUwsOWR+dx/j7gXuXgDcD/zc3R8ILuqnaNsPlrb9Bo6kxCN6NYxVHvo04XzgQjNbTugDnDp3nwJ8h9Ae1U0NDlfqGmDcT2kk/0nAPDP7iNAe1bfcPaYuytRI/linbT9A2vY//YTN/YTVgCJCb7kAUoE/ADc2eMxlxOgn9Af5e8YDa4Djgs6i/MHnaSSrtv0YuiV7/iM6Y9PMegMPAoOBLsDj7v7N8H3FhA6byQB2Aae5++Jm/zKRGKJtX2JFi512b2bjgceAC9x9bmOPF0kU2vYlSLp2iohIHEva0+5FRBJBRCVuZr3NbLqZLTGzRWb2/fDyzmb2WvgiLq+ZWafwcjOz35jZCjObb2bHNniuS8OPX25ml0bnzxJpGS287b9iZrvM7IWg/h5JPJHuidcA17r7YGAc8G0zGwL8GHjdQ6cevx7+GeBMQqciDwCuBH4HoQ0f+BkwFjge+Fn9xi8So1pk2w/7JfDV1gouySGiEnf3Te7+Qfj7UmAJ0BM4D3gk/LBHCF24hfDyv3rIu0BHM8sDTgdec/cd7r4TeI3QlcZEYlILbvu4++tAaWvml8TX5DFxMysARgFzgG7uvglCGztQf1JAT2Bdg9XWh5cdarlIzDvCbV8kKppU4mbWFngKuNpD14o45EMPsswPs1wkprXAti8SFRGXuJmlE9qI/+7uT4cXb6l/qxj+ujW8fD3Qu8HqvYCNh1kuErNaaNsXiYpIj04x4E/AEne/r8FdzxGa/YPw12cbLP9a+JP6ccDu8FvOfwOnmVmn8Aeap4WXicSkFtz2RaIiopN9zOwkYCahi6zXhRffQGhs8F9APrAWuMjdd4Q3/AcIfWhZDlzu7kXh5/p6eF2AKR6+JrNILGrhbX8mMIjQ7C3bgSvcXTsxckR0xqaISBzTGZsiInFMJS4iEsdU4iIicUwlLiISx1TiIiJxTCUuIhLHVOIiInFMJS4iEsf+P1zjoQ+RQm3HAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEiCAYAAAAcSqIJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xUdf798dc7hSSEhEDogRA6CFKDvS3Y14q6Lu6qrAX72tb1p67rusX1q+5aFhuWtayI2MDu2gtrCyR0kF4SSgglISF1Pr8/ZtCogUBmkjvlPB+PPJjcmXvnjN45ufO5d+415xwiIhKd4rwOICIizUclLyISxVTyIiJRTCUvIhLFVPIiIlFMJS8iEsVU8iIhZGZ/MrP/eJ1DZBeVvEQdM/vYzLaaWdJePHaCmX3eErlEvKCSl6hiZjnA4YADTvE0jEgYUMlLtDkP+BJ4Cjh/10Qz62Fmr5hZsZmVmNkkMxsEPAIcbGY7zGxb4LEfm9lF9eb9wda+md1vZmvNrNTMZpnZ4S314kT2lUpeos15wHOBn+PMrLOZxQNvAKuBHCALmOqcWwRcCnzhnGvjnMvYy+f4BhgOtAemAC+aWXJoX4ZIaHhW8mb2pJltMrP5IVpenZkVBH5eC8UyJbKY2WFAT2Cac24WsBw4BzgA6Abc4Jwrd85VOueaPA7vnPuPc67EOVfrnPsHkAQMCMFLEAk5L7fknwKOD+Hydjrnhgd+NBYbm84H/uuc2xz4fUpgWg9gtXOuNhRPYmbXm9kiM9seGOJpC3QIxbJFQi3Bqyd2zn0a2En2HTPrAzwIdAQqgIudc4tbPp1EGjNLAX4BxJvZhsDkJCAD2Ahkm1lCA0Xf0GlYy4HW9X7vUu95DgduBMYCC5xzPjPbClhoXolIaIXbmPxk4Crn3Cjgd8BD+zBvspnlmdmXZnZa88STMHYaUAfsh3+8fDgwCPgscN964E4zSzWzZDM7NDDfRqC7mbWqt6wCYJyZtTazvsCF9e5LA2qBYiDBzP4IpDfj6xIJimdb8j9mZm2AQ/DvxNo1OSlw3zjgzw3MVuicOy5wO9s5V2RmvYEPzWyec255c+eWsHE+8G/n3Jr6E81sEvAA/nH5B4A1+LfepwAzgQ+BBcAGM/M55zoA9wKj8f8BmIt/J+7RgUW+C7wNfIt/i/9eYG2zvjKRIJiXFw0JDNe84ZwbYmbpwBLnXNcQLPepwHJfCnZZIiKRLGyGa5xzpcBKMzsLwPyG7c28ZtZu17cbzawDcCiwsNnCiohECC8PoXwe+AIYYGbrzOxC4FfAhWY2B/9H6FP3cnGDgLzAfB8BdzrnVPIiEvM8Ha4REZHmFTbDNSIiEnoqeRGRKObJIZQdOnRwOTk5Xjy1iEjEmjVr1mbnXMd9mceTks/JySEvL8+LpxYRiVhmtnpf59FwjYhIFFPJi4hEMZW8iEgUC5tz19TU1LBu3ToqKyu9jtLikpOT6d69O4mJiV5HEZEoEzYlv27dOtLS0sjJyaHeCcqinnOOkpIS1q1bR69evbyOIyJRJmyGayorK8nMzIypggcwMzIzM2PyE4yINL+wKXkg5gp+l1h93SKy91ZtLm/SfCEpeTP7nZm5wBkgI9JRRx3Fu++++4Np9913H5dffvluH69j/UWkuTnnePbL1Zxw/2dNmj/okjezHsAx+C/GELHGjx/P1KlTfzBt6tSpjB8/3qNEIhLrNmyv5Px/f8Ot0+eTm9OuScsIxZb8vcDvafhamRHjzDPP5I033qCqqgqAVatWUVRUxJQpU8jNzWXw4MHcdtttDc7bpk2b726/9NJLTJgwAYDi4mLOOOMMRo8ezejRo5k5c2azvw4RiQ4zCgo59t5P+HplCX85dTDPXHBAk5YT1NE1ZnYK/kvwzWlsXNnMJgITAbKzs/f42NtfX8DCotJgov3Eft3Sue3kwbu9PzMzkwMOOIB33nmHU089lalTp3L22Wdz00030b59e+rq6hg7dixz585l6NChe/WcV199Nddeey2HHXYYa9as4bjjjmPRokWhekkiEoW2lldz64z5vDF3PSOyM/jnL4bTq0Nqk5fXaMmb2fvUu1p9PbcANwPH7s0TOecm479QN7m5uWG51b9ryGZXyT/55JNMmzaNyZMnU1tby/r161m4cOFel/z777/PwoXfX7uktLSUsrIy0tLSmusliEgE+2jJJm58aS5byqu54bgBXHJEbxLigxtwabTknXNHNzTdzPYHegG7tuK7A7PN7ADn3IZgQu1pi7s5nXbaaVx33XXMnj2bnTt30q5dO+655x6++eYb2rVrx4QJExo81LH+p5j69/t8Pr744gtSUlJaJL+IRKbyqlr+9tYipny1hv6d2/DkhNEMyWobkmU3+U+Ec26ec66Tcy7HOZcDrANGBlvwXmrTpg1HHXUUF1xwAePHj6e0tJTU1FTatm3Lxo0befvttxucr3PnzixatAifz8err7763fRjjz2WSZMmffd7QUFBs78GEYkseau2cML9n/H812uYeERvXrvysJAVPITRN17Dxfjx4xk3bhxTp05l4MCBjBgxgsGDB9O7d28OPfTQBue58847Oemkk+jRowdDhgxhx44dADzwwANcccUVDB06lNraWo444ggeeeSRlnw5IhKmqmrruO/9pTz6yXK6ZaQw9eKDOLB3Zsifx5NrvObm5rofH2O+aNEiBg0a1OJZwkWsv36RWLJofSnXvlDA4g1l/HJ0D/5w0n60SWp8m9vMZjnncvflubQlLyLSQup8jsc+W8E///st6SmJPHF+LmMHdW7W51TJi4i0gNUl5Vw/bQ55q7dywpAu/O30/Wmf2qrZn1clLyLSjJxzPP/1Wv765kLi44x7zx7GacOzWuycVWFV8s65mDxZlxf7RUSk+W0qreTGl+fy0ZJiDu2byd1nDqNbRsseUh02JZ+cnExJSUnMnW541/nkk5OTvY4iIiH05tz13DJ9HpU1ddx+ymDOPagncXEt321hU/Ldu3dn3bp1FBcXex2lxe26MpSIRL7tFTX88bX5zCgoYlj3tvzz7OH06dim8RmbSdiUfGJioq6MJCIR7dNvi/n9S3PZvKOKa4/uzxU/6xP0aQmCFTYlLyISqSqqa/n7W4t59svV9O3UhsnnjWJo9wyvYwEqeRGRoMxes5Xrp81h5eZyLjysFzccN4DkxHivY31HJS8i0gTVtT4e+GApD328jK5tU5hy8YEc0if8Lo6nkhcR2UdLNpRx7QsFLFxfypmjuvPHk/cjPTnR61gNUsmLiOylOp/jic9XcM+735KWnMCj547iuMENXW4jfKjkRUT2wtotFVz/4hy+XrmFY/brzN/H7U+HNklex2qUSl5EZA+cc0zLW8ufX1+ImXHPWcM4Y2TLnZYgWCp5EZHd2FRWyU0vz+ODxZs4uHcmd581lO7tWnsda5+o5EVEGvD2vPXc/Oo8yqvruPWk/fjNITmenJYgWCp5EZF6tu+s4fbXFvBKfiH7Z7Xln78YRr/OaV7HajKVvIhIwOdLN3PDS3PYVFbF1WP7ceWYviR6fFqCYKnkRSTm7ayu4//eWcxT/1tF746pvHLZIQzrER6nJQiWSl5EYlrB2m1cN62AFcXlTDgkhxuPH0hKq/A5LUGwVPIiEpNq6nz868NlPPjRMjqlJfHcRQdyaN/wOy1BsFTyIhJzlm4s49ppBcwvLGXcyCxuO3kwbVPC87QEwVLJi0jM8PkcT85cyV3vLqFNUgKP/Hokxw/p6nWsZqWSF5Go5ZyjssbH9p01bCqr5I63FvHlii0cPagTd4zbn05p0X/ZTZW8iIS1mjofpTtrKK2sZfvOmsDtmsDt2nq3v39M2c7AtMoaaurcd8tKbRXPXWcM5azc7hFzWoJgqeRFpFn5fI4d1bWU7vy+mHcVcGm9ci6tV8z1H1NRXbfH5SfEGW1TEmmbkkha4N8e7VJID9xOT04kPSWB9OREcnPa0bVtSgu98vCgkheRvVJVW8eakoofbEXvccs6UNZllTX43O6XawZpSQmkBwq5bUoiOR1af3f7u7IOFPWuabtuJyfGxcxWeVOo5EVkj3ZW1zHl6zU88slyisuqGnxMSmI86SkJ3205d05Ppn/nNNKTE360Rf3Tsk5LSojIc8JECpW8iDSoorqW575cw6OfrmDzjioO6t2em08cSPvUJNKTE36wRd0qIbK/+h/NVPIi8gMV1bU8+8VqHvtsBZt3VHNo30weHDOCA3tneh1NmiCokjezPwEXA8WBSTc7594KNpSItLzyqlqeCZT7lvJqDu/XgavH9iM3p73X0SQIodiSv9c5d08IliMiHiirrOGZL1bz+Gcr2FpRwxH9O3L12H6M6tnO62gSAhquEYlRpZU1PD1zFY9/vpLtO2v42YCO/HZsP0Zkq9yjSShK/kozOw/IA653zm0NwTJFpJls31nDv2eu5MnPV1JaWcvYgZ347dh+UXNqXfmhRkvezN4HujRw1y3Aw8BfABf49x/ABbtZzkRgIkB2dnYT44pIU22vqOGJmSv598yVlFXWcsx+nbl6bD+GZLX1Opo0I3NuD99S2JcFmeUAbzjnhjT22NzcXJeXlxeS5xWRPdtWUc0Tn6/kqZmrKKuq5bjBnfnt2H4M7qZyjzRmNss5l7sv8wR7dE1X59z6wK+nA/ODWZ6IhM6W8moe/2wFT/9vFeXVdZy4fxeuGtOPQV3TvY4mLSjYMfm7zGw4/uGaVcAlQScSkaCU7Kjisc9W8swXq9hZU8eJ+3flt2P6MaBL5F6MWpouqJJ3zp0bqiAiEpzNO6qY/OkKnv1iNZW1dZw0tBu/HdOXfp1V7rFMh1CKRLhNZZVM/mQF//lqNdW1Pk4Z1o0rx/SlbyeVu6jkRSLWptJKHvlkBc99tZqaOh+nDc/iijF96dOxjdfRJIyo5EUizIbtlTzyyXKmfL2GOp/j9BFZXPGzvvTqkOp1NAlDKnmRCLF++04e/ng5U79ei885xo30l3vPTJW77J5KXiTMFW7bycMfL2PaN+vwOcdZud25/Ki+9Gjf2utoEgFU8iJhau2WCh76eDkvzVoLwFm5Pbj8qD50b6dyl72nkhcJM2tKKnjo42W8NGsdcWb8cnQ2lx7Vh6yM2Lo2qYSGSl4kTKwuKWfSh8t4Jb+Q+DjjVwf6yz3WLjwtoaWSF/HYys3+cp9eUEhCnHHuQT257Kg+dE5P9jqaRAGVvIhHlhfvYNKHy5hRUEirhDgmHJLDJUf0ppPKXUJIJS/SwpZtKuNfHy7j9TlFtEqI48LDenHxEb3plKZyl9BTyYu0kPXbd3LHW4t5Y24RyQnxXHx4by4+ojcd2iR5HU2imEpepAWUVdZw/pNfs3bLTi49sg8XHdaLTJW7tACVvEgzq/M5rplawPLicp7+zQEc1q+D15EkhsR5HUAk2t31zmI+WLyJP528nwpeWpxKXqQZvZi3lkc/XcG5B/Xk3INzvI4jMUglL9JM8lZt4ZZX53No30z+ePJ+XseRGKWSF2kG67ZWcMmzs+iWkcyD54wkMV5vNfGG1jyRECuvquWip/OorvPx+PmjyWjdyutIEsN0dI1ICPl8jmteKGDpph38e8Jo+nbSVZrEW9qSFwmhe/67hPcWbuTWnw/iiP4dvY4jopIXCZVX89fx0MfLGX9ANucfkuN1HBFAJS8SErPXbOXGl+dxUO/2/PnUwZiZ15FEAJW8SNAKt+1k4jOz6JKezMO/GqUjaSSsaMerSBAqqmu5+Ok8qmrqeP7iA2mXqiNpJLyo5EWayOdzXPfCHBZvKOWJCaPp1znN60giP6HPlSJNdO/73/LOgg3cfOIgfjagk9dxRBqkkhdpghkFhfzrw2WcnduDCw/r5XUckd1SyYvso4K127jhpbkckNOev5w2REfSSFhTyYvsgw3bK5n4TB6d0pJ4+NcjaZWgt5CEN62hIntpZ3UdFz+TR3lVLU+cP1pXdpKIEHTJm9lVZrbEzBaY2V2hCCUSbnw+x+9enMP8ou08MH4EA7roSBqJDEEdQmlmPwNOBYY656rMTIcYSFR64MOlvDlvPTedMJCxgzp7HUdkrwW7JX8ZcKdzrgrAObcp+Egi4eXNueu57/2lnDGyOxOP6O11HJF9EmzJ9wcON7OvzOwTMxu9uwea2UQzyzOzvOLi4iCfVqRlzFu3netfLGBUz3bcMU5H0kjkaXS4xszeB7o0cNctgfnbAQcBo4FpZtbbOed+/GDn3GRgMkBubu5P7hcJNxtLK7nomW/ITE3i0XNHkZQQ73UkkX3WaMk7547e3X1mdhnwSqDUvzYzH9AB0Ka6RLTKmjomPpNHWWUtL192CB10JI1EqGCHa6YDYwDMrD/QCtgcbCgRLznnuOGlucwt3M59Zw9nUNd0ryOJNFmwJyh7EnjSzOYD1cD5DQ3ViESSSR8u4/U5Rfz++AEcO7ihkUqRyBFUyTvnqoFfhyiLiOfemb+ef7z3LaePyOKyI/t4HUckaPrGq0jA/MLtXPvCHEZkZ/D3cfvrSBqJCip5EWBTWSUXP5NHRutEHj13FMmJOpJGooMuGiIxz38kzSy2VdTw4qUH0ykt2etIIiGjkpeY5pzjplfmUbB2G4/8eiRDstp6HUkkpDRcIzHt4U+W82p+Idcf05/jh3T1Oo5IyKnkJWb9d8EG7n53CScP68aVY/p6HUekWajkJSYtWl/KNS8UMDSrLXefOVRH0kjUUslLzNm8o4qLns4jPTmRyefl6kgaiWra8Soxpaq2jkufnUVJeRUvXnIIndN1JI1EN5W8xAznHDe/Mp+81VuZdM4I9u+uI2kk+mm4RmLGY5+t4OXZ67h6bD9OGtrN6zgiLUIlLzHhg0Ub+fvbi/n5/l25emw/r+OItBiVvES9JRvK+O3z+Qzuls49Zw0jLk5H0kjsUMlLVCvZUcWFT39DalICj52XS0orHUkjsUU7XiVqVdf6uOw/sykuq+KFSw6ma9sUryOJtDiVvEQl5xy3Tp/P16u2cP8vhzO8R4bXkUQ8oeEaiUpPfL6SF/LWctWYvpw6PMvrOCKeUclL1PloySbueGsRxw3uzLVH9/c6joinVPISVZZuLOO3U/IZ2CWde88eriNpJOap5CVqbC2v5qJn8khKjOex83Np3Uq7nERU8hIVqmt9XPbcLNZvq+TRc0eRlaEjaURAR9dIFHDOcdtrC/hyxRbuPXsYo3q28zqSSNjQlrxEvKf/t4rnv17DZUf14fQR3b2OIxJWVPIS0T79tpg/v7GQY/brzA3HDvA6jkjYUclLxFq2aQdXTJlN/85pOpJGZDdU8hKRtlVUc9HT39AqPo7Hz8+lTZJ2L4k0RO8MiTg1dT6umDKbom2VTLn4QLq3a+11JJGwpZKXiPPn1xcyc1kJd585lNyc9l7HEQlrGq6RiPLsF6t49svVTDyiN2fl9vA6jkjYU8lLxPhyRQl/en0hYwd24sbjB3odRyQiBDVcY2YvALuOW8sAtjnnhgedSuRHisuquOr5fHq2b819vxxOvI6kEdkrQZW8c+7sXbfN7B/A9qATifxInc9x9dR8SnfW8OyFB5CWnOh1JJGIEZIdr2ZmwC+AMaFYnkh9D3ywlP8tL+GuM4YysEu613FEIkqoxuQPBzY655aGaHkiAHy2tJgHPlzKGSO7c1auTlkgsq8a3ZI3s/eBLg3cdYtzbkbg9njg+UaWMxGYCJCdnb2PMSUWbSyt5JqpBfTt2Ia/nDYY/wdGEdkXjZa8c+7oPd1vZgnAOGBUI8uZDEwGyM3NdfuQUWJQbZ2Pq57Pp6K6jqkTR+rc8CJNFIp3ztHAYufcuhAsSwSAf773LV+v9J86uF/nNK/jiESsUIzJ/5JGhmpE9sVHSzbx0MfL+eXoHjp1sEiQgt6Sd85NCEEOEQCKtu3k2hcKGNgljT+dMtjrOCIRT994lbBRU+fjyimzqan18dCvRpKcGO91JJGIp71ZEjbufncJs9ds41/jR9C7Yxuv44hEBW3JS1h4b+FGJn+6gnMP6snJw7p5HUckaqjkxXNrt1Rw/bQChmSl84eTBnkdRySqqOTFU9W1/nF45+DBc0aSlKBxeJFQ0pi8eOqOtxYxZ912Hvn1SHpmpnodRyTqaEtePPP2vPU89b9V/ObQHI4f0tXrOCJRSSUvnlhdUs7vX5rLsB4Z3HSCxuFFmotKXlpcZU0dlz83m7g4Y9L4EbRK0Goo0lw0Ji8t7q9vLmRBUSmPn5dLj/atvY4jEtW0CSUtakZBIf/5cg2XHNGbo/fr7HUckainkpcWs7x4Bze/Mo9RPdvxu+MGND6DiARNJS8torKmjiuem02rhDgmnTOCxHiteiItQWPy0iJum7GAxRvKeOo3o+naNsXrOCIxQ5tT0uxenrWOF/LWcsXP+nDUgE5exxGJKSp5aVZLN5bxh+nzObBXe649ur/XcURijkpemk1FdS2XPzeb1KR4/jV+BAkahxdpcRqTl2bhnOMPr85nWfEO/nPhgXRKT/Y6kkhM0qaVNItpeWt5Jb+Qq8f249C+HbyOIxKzVPIScovWl/LHGQs4rG8HrhrTz+s4IjFNJS8htaOqliuem016SiL3nj2c+DjzOpJITFPJS8g457jplXmsKinnX+NH0DEtyetIIjFPJS8h89xXa3h9ThHXHzuAg3pneh1HRFDJS4jML9zOn19fyJH9O3LZkX28jiMiASp5CVppZQ2XPzebzDatuPfs4cRpHF4kbOg4eQmKc44bX5pL4badvDDxINqntvI6kojUoy15CcpT/1vF2/M3cOPxA8jNae91HBH5EZW8NFnB2m3c8dYijh7UiYsP7+11HBFpgEpemmR7RQ1XPDebTmnJ3HPWMMw0Di8SjjQmL/vMOcf1L85hU1kl0y45mIzWGocXCVfakpd99vhnK3l/0UZuOmEQI7LbeR1HRPYgqJI3s+Fm9qWZFZhZnpkdEKpgEp5mrd7Cne8s5vjBXfjNoTlexxGRRgS7JX8XcLtzbjjwx8DvEqW2lFdz5ZR8sjJSuOusoRqHF4kAwY7JOyA9cLstUBTk8iRM+XyO66YVULKjmlcuP4T05ESvI4nIXgi25K8B3jWze/B/Kjhkdw80s4nARIDs7Owgn1Za2sOfLOfjJcX85bQhDMlq63UcEdlLjZa8mb0PdGngrluAscC1zrmXzewXwBPA0Q0txzk3GZgMkJub65qcWFrcVytK+Md/l3DS0K78+kD9gRaJJI2WvHOuwdIGMLNngKsDv74IPB6iXBImNu+o4qrn8+mZmcrfx+2vcXiRCBPsjtci4MjA7THA0iCXJ2Gkzue4ZmoB23fW8OA5I0nTOLxIxAl2TP5i4H4zSwAqCYy5S3SY9OEyPl+2mTvH7c9+3dIbn0FEwk5QJe+c+xwYFaIsEkZmLtvMfR98y7gRWZw9uofXcUSkifSNV/mJTaWVXD01nz4d2/DX04doHF4kguncNfIDtXU+rno+n/KqOqZcPJLWrbSKiEQyvYPlB+57fylfrdzCPWcNo3/nNK/jiEiQNFwj3/nk22Ie/HgZv8jtzpmjunsdR0RCQCUvAKzfvpNrXyigf6c0bj9liNdxRCREVPJCTZ2Pq6bkU1VTx0O/HklKq3ivI4lIiGhMXrjnv0vIW72V+385nD4d23gdR0RCSFvyMe6DRRt59JMV/OrAbE4dnuV1HBEJMZV8DFu3tYLrps1hcLd0bj1pP6/jiEgzUMnHqOpaH1dOyafO53jwnJEkJ2ocXiQaaUw+Rt359mIK1m7joV+NJKdDqtdxRKSZaEs+Br0zfwNPzlzJhENyOHH/rl7HEZFmpJKPMWtKKrjhpTkM696Wm04c6HUcEWlmGq6JEc45Zq/Zyq3TF2DApHNGkpSgcXiRaKeSj3LLi3cwI7+Q6QVFrNlSQUpiPJPOGUGP9q29jiYiLUAlH4WKy6p4fU4R0wsKmbtuO3EGh/btwNVj+3HckC60SdL/dpFYoXd7lCivquW/Czfwan4RM5dtps7nGJKVzh9+PohThnWjU3qy1xFFxAMq+QhWW+fjs2WbmZFfyLsLNrKzpo6sjBQuPbI3pw3Pop9OFSwS81TyEcY5x5x125meX8gbc4vYvKOatimJnD4yi9NHZDEqux1xcbqSk4j4qeQjxOqScqbn+8fZV24up1VCHEcP6sRpw7M4ckBHHSkjIg1SyYexkh1VvDlvPa/mF5K/ZhtmcFCvTC49sjfHD+lK25REryOKSJhTyYeZndV1vLdoI9PzC/n022JqfY6BXdL4fycM5JRh3eiWkeJ1RBGJICr5MFDnc3yxvIRX8wt5Z/56yqvr6JKezIWH9+K04VkM6prudUQRiVAqeY8451hQVMr0/EJem1PEprIq0pISOGloN04d0Y2DemVqB6qIBM2Tkl9QVMrx931KdvvW5HRIJbt9a3pmtqZn+1S6ZSSTEB+9p9RZu6WC1+YU8Wp+Ics27SAx3jhqQCdOH5HFmIGddMpfEQkpT0q+XWoiWRkprNhczsffFlNd6/s+UJzRvV0K2Zmp9NxV/pmp9MxsTXb71hFZgtsqqnlz3npm5Bfx9aotAIzOacffTh/Cz/fvSkbrVh4nFJFo5UnJd2ubwhMTRgPg8zk2lFayuqSCNVvKWV1S4f/ZUk7+mq2UVdb+YN4u6clkZ7amZwOfAtq2Dp+jTSpr6vho8SZezS/koyWbqKlz9OmYyg3HDeCUYd107hgRaRGej8nHxRndMlLolpHCwX0yf3Cfc45tFTWsKilnzZaK7/8AlPg/ARTPWveDx2e0TqRn+9YNfgrolJaEWfOOcft8jq9WbmF6fiFvzV9PWWUtHdOSOO/gHE4fkcXgbunNnkFEpD7PS35PzIx2qa1ol9qKEdntfnJ/RXUta7ZUsGrzDz8FFKzdyptzi/C57x+bnBhHz/ap330KqP8HICsjJaj9AIs3lPJqfiGvFRSxfnslqa3iOW5IF04fkcUhfToQrx2oIuKRsC75xrRulcDALukM7PLTQwxr6nwUbt35k08BqzaX8+m3xVTV2w8Qv2s/QL2hn11/BLLbtyal1U/3A6zfvpMZBUVMzy9k8YYy4uOMI/t35KYTB3HMoM4Nzpsq4MsAAAhjSURBVCMi0tKCKnkzGwY8ArQBVgG/cs6VhiBX0BLj48jpkNrg9Ut9PsemsipWl5R/N/7v3ydQwWsFRZT+aD9Ap7QkcjL9nwK6ZaTwzcotfLmyBOdgRHYGt58ymJOGdiWzTVJLvTwRkb1izrnGH7W7mc2+AX7nnPvEzC4Aejnnbm1svtzcXJeXl9fk521u2yqqA+VfwerN5azeUsGawB+DjaVV9OqQyqnDu3Ha8CxdBFtEWoyZzXLO5e7LPMEO1wwAPg3cfg94F2i05MNdRutWZLRuxbAeGT+5r6q2jlbxcdqBKiIRIdhvHc0HTgncPgvoEeTywl5SQrwKXkQiRqNb8mb2PtClgbtuAS4AHjCzPwKvAdV7WM5EYCJAdnZ2k8KKiMi+CWpM/gcLMusP/Mc5d0Bjjw33MXkRkXDUlDH5oIZrzKxT4N844A/4j7QREZEwEeyY/Hgz+xZYDBQB/w4+koiIhEpQR9c45+4H7g9RFhERCbHoPaeviIio5EVEolnIjq7Zpyc1KwOWtPgTh05bYLvXIYIQyfkjOTtAB2Cz1yGCEOn//SM9fz/nXNt9mcGrE5Qt2dfDgMKJmU12zk30OkdTRXL+SM4OYGZ5Wve9Ew3593UeDdc0zeteBwhSJOeP5OzRINL/+8dcfq+GayJ6a0akqbTuS0vzakt+nz9yiEQJrfvSojzZkhcRkZahMfk9MLPuZjbDzJaa2Qozm2RmSWZ2jJnNMrN5gX/HeJ21IXvIf4CZFQR+5pjZ6V5nbcju8te7P9vMdpjZ77zMGa0ief3Xuv89lfxumP98wq8A051z/YB+QApwF/5D4E52zu0PnA8861nQ3Wgk/3wg1zk3HDgeeNTMwupSkI3k3+Ve4G0P4kW9SF7/te7/iHOu2X6A7sAMYCmwApgEJAGZwEfADmBSc2YIIvtY4NMfTUsHtgJt6k0zoARI8jpzE/P3AjYCCV5n3pf8wGnA3cCf8F+dzPPMDbwGrf/hnT0m1v1m25Jv5K9RJf4rSIXzx+zBwKz6E5z/+rWrgL71Jp8B5Dvnqlou2l7ZY34zO9DMFgDzgEudc7U/XYSn9pR/GHAjcHvLx9o7Wv89pXW/nuYcrhkDVDrn/g3gnKsDrgXOw7/D93P8K3u4MqChvdLfXRbKzAYD/wdc0lKh9sEe8zvnvnLODQZGAzeZWXJLhtsLe8p/O3Cvc25Hy0baJ1r/vaN1v57mLPm93RIIVwuAHxzPbGbpQGdgiZl1B14FznPOLfcgX2P2mH/XNOfcIqAcGNKi6Rq3p/xtgbvMbBVwDXCzmV3Z4gn3TOu/d7Tu19OcJd/olkCY+wBobWbnAZhZPPAPvh9XfRO4yTk307uIe7Sn/F127Wwys574L8i+yqOcu7Pb/M650c65HOdcDnAfcIdzbpJ3URuk9d87Wvfrac6S36u/puHK+fd2nA6caWZL8e9c8jnn/gZciX9r7NZ6h2N18jDuTzSS/zBgjpkV4N8au9w5F1YnzWokfyTQ+u8Rrfs/XWBz7SE2IA//xzmAeOAx4JZ6j5lAmB5d0MDrOQRYDYzyOovye59nL/Jq/Q+Tn0jOHor8zfqNVzPrATwIDAI6Ai845y4J3LcK/2FBrYBtwLHOuYXNFkakhWn9l3DQYqc1MLNDgOeBcc65WY09XiSaaP0Xr+jcNSIiUUynNRARiWIhKXkz62FmH5nZIjNbYGZXB6a3N7P3AifZec/M2gWmm5k9YGbLzGyumY2st6zzA49fambnhyKfSHMJ8br/jpltM7M3vHo9En1CtSVfC1zvnBsEHARcYWb7Af8P+MD5v9b9QeB3gBPwf827HzAReBj8bwzgNuBA4ADgtl1vDpEwFZJ1P+Bu4NyWCi6xISQl75xb75ybHbhdBiwCsoBTgacDD3sa/4l1CEx/xvl9CWSYWVfgOOA959wW59xW4D38Z4oTCUshXPdxzn0AlLVkfol+IR+TN7McYATwFdDZObce/G8GYNcXJrKAtfVmWxeYtrvpImEvyHVfpFmEtOTNrA3wMnCN85+nY7cPbWCa28N0kbAWgnVfpFmErOTNLBH/Sv6cc+6VwOSNuz6KBv7dFJi+DuhRb/buQNEepouErRCt+yLNIlRH1xjwBLDIOffPene9hv/KMQT+nVFv+nmBIw0OArYHPtK+CxxrZu0CO1yPDUwTCUshXPdFmkVIvgxlZocBn+E/Cb8vMPlm/GOT04BsYA1wlnNuS+CNMQn/TtUK4DfOubzAsi4IzAvwNxc4H7dIOArxuv8ZMBD/1X9KgAudc9rIkaDoG68iIlFM33gVEYliKnkRkSimkhcRiWIqeRGRKKaSFxGJYip5iXpmVhe4DukCM5tjZteZ2R7XfTPLMbNzWiqjSHNRyUss2OmcG+6cGwwcA5yI/2yne5IDqOQl4uk4eYl6ZrbDOdem3u+9gW+ADkBP4FkgNXD3lc65/5nZl/ivzboS/1kkHwDuBI4CkoAHnXOPttiLEGkilbxEvR+XfGDaVvzfLi0DfM65SjPrBzzvnMs1s6OA3znnTgo8fiLQyTn3VzNLAmbi/xbryhZ9MSL7KMHrACIe2XU2yERgkpkNB+qA/rt5/LHAUDM7M/B7W/wX/lDJS1hTyUvMCQzX1OE/M+RtwEZgGP59VJW7mw24SueSkUijHa8SU8ysI/AIMMn5xyrbAuudcz78l96LDzy0DEirN+u7wGWB0wpjZv3NLBWRMKcteYkFKWZWgH9ophb/jtZdpwV+CHjZzM4CPgLKA9PnArVmNgd4Crgf/xE3swNnkizm+0v6iYQt7XgVEYliGq4REYliKnkRkSimkhcRiWIqeRGRKKaSFxGJYip5EZEoppIXEYliKnkRkSj2/wHx/d6yDaWHWwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_pred.plot(title=\"Predictions\")\n",
    "changes[start_date:end_date].plot(title=\"Actual\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform out-of-sample prediction of data past the end of the data we have\n",
    "y_new = results.forecast(steps=10)  # You can also use the predict method with a start/end date in the future:\n",
    "# y_new = results.predict(len(changes), end=len(changes)+10, dynamic=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.forecast?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2153f372f88>]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3TU93nn8fejGwIkBEJCkgVCSNyRbMAyNsIxNvgCJI0T2m7TJt7GTZY2TbNx66Rpvbs9u83paXLS49PsaVOXjdvUrU+6m5ikaYrwPb5IGBswtgTCtsQdRkgjgZDERZd59o8Z25gINAKJuX1e53AyzHznNw8T+Pin7+/7e77m7oiISOJLi3UBIiIyNhToIiJJQoEuIpIkFOgiIklCgS4ikiQU6CIiSSKmgW5m/2Bm7WbWNEbHGzKzPZFfPxuLY4qIJAqL5Tp0M7sD6AWecPeqMTher7vnXHtlIiKJJ6Zn6O7+MtB18XNmVmlm28xsl5m9YmYLY1SeiEhCicc59M3AV9z9ZuBrwPdG8d5sM9tpZq+Z2afGpzwRkfiUEesCLmZmOUAt8CMze//pCZHXNgJ/Pszbjrv7fZHHZe5+wswqgBfMrNHdW8e7bhGReBBXgU74J4bT7r700hfcfQuw5UpvdvcTkf89YGa/AJYBCnQRSQlxNeXi7meAg2b26wAWdlM07zWzaWb2/tl8AbAK2DduxYqIxJlYL1v8IbAdWGBmx8zsC8BngS+Y2VvAXuD+KA+3CNgZed+LwLfcXYEuIikjpssWRURk7Ix4hm5ms8zsRTNrNrO9ZvbVYcZMM7OfmNnbZva6mV3zmnIRERmdEc/QzawEKHH33WaWC+wCPnXxdIaZfQfodff/FVk3/rfuvvZKxy0oKPDy8vJr/gOIiKSSXbt2Bd29cLjXRlzl4u4BIBB53GNmzUApH73guBj4y8iY/WZWbmZF7n7ycsctLy9n586do/hjiIiImR2+3GujuihqZuWElwLuuOSlt4CNkTErgNnAzGHevyly48/Ojo6O0Xy0iIiMIOpAj9z08xTwUGR54cW+BUwzsz3AV4A3gcFLj+Hum929xt1rCguH/YlBRESuUlQ3FplZJuEwfzJyg89HRAL+wchYAw5GfomIyHUSzSoXAx4Hmt390cuMmWpmWZHffhF4eZizeBERGUfRnKGvAh4AGiNTKgCPAGUA7v4Y4Zt6njCzIcIXS78wDrWKiMgVRLPK5VXARhizHZg3VkWJiMjoxVUvFxERuXoKdBGR6+RwZx+PvdRKQ2twXI4fb+1zRUSSSkt7D3WNbWxtaqM5EF4r8qU7K6mtLBjzz1Kgi4iMIXenOdBDXVOAuqY2Wtp7AVheNpX/tmER66qKmZU/aVw+W4EuInKN3J23jnVT1xRgW1MbhzvPkmZwS3k+D3xyCfctKaY4L3vc61Cgi4hchVDI2XXkFFsbAzzd1MaJ7vNkpBkrK6fzu3dUcu+SIgpyJlzXmhToIiJRGhwKseNgF3VNAZ7ee5KOngtkpafxsXkF/NG9C7h70QymTsoa+UDjRIEuInIF/YMh6luD1DUGeHbfSU6dHSA7M427FsxgXVUxaxbOIDc7M9ZlAgp0EZFfcn5giJfe7WBbUxvPNZ+k5/wgORMyWLNwBhuqi1k9fwYTs9JjXeYvUaCLiAB9FwZ58Z126hrbePGdds72D5E3MZP7lhSzvqqYVXMLyM6MvxC/mAJdRFJW97kBnm8+SV1TGy+/28GFwRAFOVncv7SUDdXF3FYxncz0xLn/UoEuIimlq6+fZ/e1sbWxjYbWIANDTvGUbH5zRRnrqoq5pTyf9LQrtq+KWwp0EUl67WfO8/TeNuqa2thxsIuhkDNz2kQ+X1vO+uoSls6cSlqChvjFFOgikpSOnz7HtqY26hoD7DpyCneoKJjM762uYH1VCUtumEJ4u4fkoUAXkaRxKNhHXVMb25oCvHWsG4CFxbl8de08NlSXMG9GTtKF+MUU6CKS0N472cPWxjbqmgLsb+sB4MaZefzxugWsryphTsHkGFd4/SjQRSShuDt7T5wJT6c0BWjt6APg5tnT+O8fX8R9S8av+VW8U6CLSNwLhZy3jp2OTKe0caQr3Pzq1jnT+e3acu5bUkzRlPFvfhXvFOgiEpeGQs7OQ13UNbXx9N42ApHmV7VzC/jSnZXcu7iI6de5+VW8U6CLSNwYGAqx40AXW5sCPLP3JMHeC2RlpHHHvEK+du8C7l5URN6k+OibEo9GDHQzmwU8ARQDIWCzu3/3kjF5wL8AZZFj/pW7/+PYlysiyebC4BD1LUHqGtt4tvkkp88OMDEznbsWFrKuqoQ1C2eQM0HnntGI5lsaBB52991mlgvsMrNn3X3fRWO+DOxz918xs0LgHTN70t37x6NoEUls5/rDza/qmgK80NxOz4VBcidksHbRDNZVlbB6fmFcNr+KdyMGursHgEDkcY+ZNQOlwMWB7kCuhRd45gBdhP9DICICQO+FQV7Y3862pgAv7u/g3MAQUydlsr66mPVVJdTOnc6EDIX4tRjVzzFmVg4sA3Zc8tLfAD8DTgC5wG+4e2gM6hORBNZ9doDnmk9S1xTg5feC9A+GKMiZwMblpayvKuHWivyEan4V76IOdDPLAZ4CHnL3M5e8fB+wB1gDVALPmtkrl44zs03AJoCysrJrqVtE4lRn7wWe2RfuYNjQEmQw5JTkZfNbK8rYUF3CzbOnJWzzq3gXVaCbWSbhMH/S3bcMM+RB4Fvu7kCLmR0EFgKvXzzI3TcDmwFqamr8WgoXkfhx8v3mV41t7DjYScihLH8SX7h9DuuqirkpSZpfxbtoVrkY8DjQ7O6PXmbYEWAt8IqZFQELgANjVqWIxJ1jp85G7tZsY3ek+VVl4WR+/865rK8uZnFJ8jW/infRnKGvAh4AGs1sT+S5RwgvUcTdHwO+CfzAzBoBA77h7sFxqFdEYuhgsI+6pgDbmtp4O9L8alHJFP7w7vmsrypmXlFujCtMbdGscnmVcEhfacwJ4N6xKkpE4oO78157L3WXNL+6aWYe31i3kPVVxZSnUPOreKfV+iLyERc3v9raFOBARx9mcHNZuPnVuqpiZk5LzeZX8U6BLiK4O3uOnv5gTvz95le3VUznwUjzqxlqfhX3FOgiKWoo5Ow6fIqtjYEPml9lphu1lQX8/p2V3KPmVwlHgS6SQgaHQrx2oIu6pgBPX9L86uv3LWDtoiLyJqr5VaJSoIskuQuDQzS0dFLXFOCZfR82v1qzcAbrqoq5S82vkob+XxRJQucHws2vtjW18dy+k2p+lSIU6CJJ4nLNr9ZVFbOhWs2vUoECXSSBdZ8d4Nnmk2y7qPlVYe4EfvXmUtYtUfOrVKNAF0kwHT0XeGZfeG/N7a2dDIac0qkT+dyts1lfXczyMjW/SlUKdJEEEOg+98Ea8Z2Hugg5lE+fxBc/VsH6qmJunJmnvimiQBeJV4c7+z4I8T1HTwOwoCiXr6yZx7qqYhYW5yrE5SMU6CJx5L2TPdRFQrw5EN5OoLo0j6/ft4D1VcVUFObEuEKJZwp0kRi6uG9KXVOA1o4+AG6eHe6bct+SYmblq2+KREeBLnKdhULOm0dPhzeEaApwtOvcB31TPl9bzr1LiilS3xS5Cgp0ketgKOS8frCLbZFb7tvOhPumrJpbwB/cNZd7FheTPzkr1mVKglOgi4yT/sEQ2w90sq0pwDN7T9LZ18+EjDRWzy/kG9ULWLNQfVNkbCnQRcbQ+YEhXn0vyNamAM/tO8mZ84NMzkpnzaIi1lcVs3p+IZPVN0XGif5miVyj9/um1DUGeK65nd4Lg0zJzuCexcWsryrm9nkFZGfqlnsZfwp0katwfmCIX7zTzn80tvFC80n6+sN9Uz5eXcL66mJqKwvIytAt93J9KdBFonS2f5AX93ewtSnAi/vbOds/RP7kLD659AY2VJdwW8V09U2RmFKgi1xBX6SD4dbGAC++0875gRAFOVl8elkpG6pLuHVOPhkKcYkTIwa6mc0CngCKgRCw2d2/e8mYrwOfveiYi4BCd+8a23JFxl/vhUGebz7J1sYAv3ingwuDIQpyJvDrN89ifXUxt86ZruZXEpeiOUMfBB52991mlgvsMrNn3X3f+wPc/TvAdwDM7FeAP1SYSyI5c34gEuJtvPRuB/2DIWbkTuAzt8xiQ3UJNeX5CnGJeyMGursHgEDkcY+ZNQOlwL7LvOU3gR+OWYUi46T73ADP7Qufib/yXpD+oRDFU7L57K1lbKgu4eayaaQpxCWBjGoO3czKgWXAjsu8PglYB/zBZV7fBGwCKCsrG81Hi4yJ02f7eWbfSeoaA7zaEmRgyLkhL5sHVs5mQ3UJy2ZNVYhLwoo60M0sB3gKeMjdz1xm2K8A9ZebbnH3zcBmgJqaGh9lrSJX5VRfP8/sa2NrYxv1LcEPNoT4fG05G6pLWDprqtrQSlKIKtDNLJNwmD/p7luuMPQzaLpF4kBn7wWe3nuSuqYADa2dDIWcWfkT+cLH5rChqkQbQkhSimaViwGPA83u/ugVxuUBq4HPjV15ItELdJ/juebwJsmvHehiKOTMnj6JTXdU8PHqEpbcMEUhLkktmjP0VcADQKOZ7Yk89whQBuDuj0We+zTwjLv3jXmVIsMIhZy3j3fzfPNJnm9uZ19kQ4g5BZP5vdUVbKguYXGJQlxSRzSrXF4FRvwX4e4/AH5w7SWJXF7vhUFefa+D55vbefGddoK9/aRZeEOIP1m/kDULZzBvRo5CXFKS7hSVuHe062z4LHx/OzsOdNE/FCI3O4M7F8xg7cIZrJ5fyDT1EhdRoEv8GRwK8ebR0zzf3M4L+0/y7sleACoKJ/PbtbNZs7CImvJp6psicgkFusSF7nMDvPxuBy/sD0+lnD47QEaasWJOPv+pZhZrFxUxp2ByrMsUiWsKdImZAx29vLC/neeb23njUBeDIWfapEzWLJjB2kVFfGx+AVOytaOPSLQU6HLdDAyFeONQFy80t/P8/nYOBsMLohYW57LpjgrWLprB0lnT1DNF5Cop0GVcnerr5xfvtvNcczsvv9NBz4VBstLTWFk5nQdXlXPXghnMyp8U6zJFkoICXcaUu/Neey/PNZ/kheZ2dh85RcihMHcCG6pLWLNoBrfPLdC+miLjQP+q5Jq5OzsOdlHXGOD5/e0cO3UOgKrSKXxlzTzWLppB1Q15anolMs4U6HLVBodC1DW18fcvt9J0/AzZmWncPreAL981l7sWzKA4LzvWJYqkFAW6jNq5/iF+vOso/+eVgxzpOktFwWT+cmM1n1paysQs7W4vEisK9Ks0FHJ6LwySNzF1ltWd6uvnie2H+afth+jq62fprKk8smER9ywu0soUkTigQL9Kj73Uyl898w4rK6azcflM1lUVk5OkF/qOdp3l8VcP8n/fOMq5gSHWLpzB766u5JbyaeqZIhJHkjOBroPnmk9SPCWb46fP8bUfvcX/+GkT66qK2bi8lNrKgqQ4Y917opvNLx/g528HMOD+paVsuqOCBcW5sS5NRIahQL8KPecHePtYN19aXcnD985n95FTPLX7OD9/6wQ/efM4RVMm8KllpWxcNjPhws/d2d7ayd+91Mor7wWZnJXO76wq58FVc7hh6sRYlyciV6BAvwo7Ipsn1M6djplx8+x8bp6dz599YjEv7G9ny+5jPP7KQf7+pQNUlU5h47KZfHLpDRTkTIh16Zc1OBRi2942/v6lAzQe76YgZwJfv28Bn7ttdkpdJxBJZAr0q1DfGmRCRhrLy6Z95PnszHQ2VJewobqEYO8F/j1yxv7nP9/HX2xtZvX8QjYuL+XuRUVkZ8bHapDzA0P8aOeHK1bmRFasfHpZadzUKCLRUaBfhYaWTm4pz79i4BXkTODBVXN4cNUc3jvZw5Y3j/OT3cd5YX87udkZfOLGEjYun0nN7NhcWDzV188/v3aYf2o4RKdWrIgkBQX6KHX0XOCdkz3cv+yGqN8zryiXb6xbyNfuXcBrBzp5avcx/m3PCX74+lHK8ifx6WWlbFxeyuzp498e9tip8IqVf309vGJlzcIZ/O4dFayYk68VKyIJToE+Sg2tQQBWVRaM+r3pacaquQWsmlvAN+8f5Om9bWzZfZz//cJ7fPf596iZPY2Ny2fy8RtLxnzeet+JM2x+uZV/14oVkaSlQB+lhpZOpmRnUFWad03HmTwhg43LZ7Jx+UwC3ef46ZsneGr3MR75SSP/89/3cs+iIjYuL+WO+YVXvTPP+ytWHnv5AC+/28HkrHQerC3nd27XihWRZDRioJvZLOAJoBgIAZvd/bvDjLsT+GsgEwi6++qxLTU+1LcGua1i+pjOM5fkTeRLd1bye6sraDp+hqd2H+Nnb53gPxoDTJ+cxSeX3sCvLp/Jkhui28F+KORsi/RYefvYRStWbp1N3iStWBFJVtGcoQ8CD7v7bjPLBXaZ2bPuvu/9AWY2FfgesM7dj5jZjHGqN6aOdJ7l2Klz/JePVYzL8c2M6pl5VM/M4799fBEvvdPBljeP8eRrR/jH+kPML8ph4/KZfGpp6bCNr84PDPGjXcf4/isHONypFSsiqWbEQHf3ABCIPO4xs2agFNh30bDfAra4+5HIuPZxqDXm6t+fP587fdw/KzM9jbsXF3H34iK6zw7w88YTbNl9nG/V7efb2/azqrKAjctLuW9JMQNDIf55+2F+EFmxctOsqfzp+oXcs7hYK1ZEUsio5tDNrBxYBuy45KX5QKaZ/QLIBb7r7k8M8/5NwCaAsrKy0VcbY/UtQWbkTqCyMOe6fm7epEw+e+tsPnvrbA4F+9jy5nG27D7GH/2/t5iU1QTA2X6tWBFJdVEHupnlAE8BD7n7mWGOczOwFpgIbDez19z93YsHuftmYDNATU2NX0vh11soFL7AeMf8wpiGZXnBZP7onvk8tHYeOw+f4idvHgec364tZ2HxlJjVJSKxF1Wgm1km4TB/0t23DDPkGOELoX1An5m9DNwEvDvM2IT0bnsPnX391FaO/3RLNNLSjBVz8lkxJz/WpYhInBhxPZyFT0cfB5rd/dHLDPs34GNmlmFmk4BbgeaxKzP26ls6AVg1d/Trz0VErodoztBXAQ8AjWa2J/LcI0AZgLs/5u7NZrYNeJvw0sbvu3vTeBQcKw0tQeYUTNb6bRGJW9GscnkVGHHS2N2/A3xnLIqKN4NDIXYc7OKTS6O/3V9E5Hq7ulsQU8xbx7rpvTB4Vbf7i4hcLwr0KDS0hNefr4yTC6IiIsNRoEehvjXI4pIp5E/OinUpIiKXpUAfwbn+IXYfPn1d7g4VEbkWCvQR7DzcRf9QiFotVxSROKdAH0F9SycZacaKct3AIyLxTYE+gobWIMvKpjJ5glrHi0h8U6BfQffZARqPd1Or5YoikgAU6Few/UAn7rrdX0QSgwL9Chpag0zMTGfprKmxLkVEZEQK9CuobwmyYk4+WRn6mkQk/impLqOt+zytHX1afy4iCUOBfhkNke3mdEFURBKFAv0y6ls6mTYpk8Ul2gVIRBKDAn0Y7s721iArK6eTpk2WRSRBKNCHcajzLCe6z2u6RUQSigJ9GPWRdrlafy4iiUSBPoyG1iA35GVTPn1SrEsREYmaAv0SoZCzvbWT2rkFhPfHFhFJDAr0S+wLnOHU2QFqtTuRiCSYEQPdzGaZ2Ytm1mxme83sq8OMudPMus1sT+TXn41PuePv/fXnmj8XkUQTTU/YQeBhd99tZrnALjN71t33XTLuFXf/xNiXeH3Vt3RSWTiZoinZsS5FRGRURjxDd/eAu++OPO4BmoHS8S4sFvoHQ7x+sEtn5yKSkEY1h25m5cAyYMcwL680s7fMrM7MloxBbdfdnqOnOTcwpPXnIpKQot6Gx8xygKeAh9z9zCUv7wZmu3uvmW0AfgrMG+YYm4BNAGVlZVdd9HipbwmSZrCyQhdERSTxRHWGbmaZhMP8SXffcunr7n7G3Xsjj7cCmWb2S6e57r7Z3WvcvaawsPAaSx97Da1BqkrzyJuUGetSRERGLZpVLgY8DjS7+6OXGVMcGYeZrYgct3MsCx1vfRcGefPIaU23iEjCimbKZRXwANBoZnsizz0ClAG4+2PArwFfMrNB4BzwGXf3cah33Lx+qIvBkKv/uYgkrBED3d1fBa54y6S7/w3wN2NVVCw0tATJSk+jZnZ+rEsREbkqulM0or6lk+WzpzIxKz3WpYiIXBUFOtDV18++wBlWaf5cRBKYAh3Y3hq+flurG4pEJIEp0IH61iA5EzK4aWZerEsREblqCnTCZ+i3zsknI11fh4gkrpRPsBOnz3Ew2KfpFhFJeCkf6B9uN6f15yKS2FI+0BtaOynIyWJBUW6sSxERuSYpHejuTn1LkJWV2m5ORBJfSgd6a0cv7T0XWKXt5kQkCaR0oNe3RNaf64YiEUkCKR7oQWZOm0jZ9EmxLkVE5JqlbKAPhZzXDnTqdn8RSRopG+hNx7s5c36QWi1XFJEkkbKBXt8aXn+u+XMRSRYpG+gNLZ0sKMqlMHdCrEsRERkTKRno5weGeONQl6ZbRCSppGSg7z5yiguDIV0QFZGkkpKB3tDSSXqacWuFtpsTkeSRkoFe3xrkxpl55GZnxroUEZExk3KB3nN+gLePdWu6RUSSzoiBbmazzOxFM2s2s71m9tUrjL3FzIbM7NfGtsyxs+NAF0Mh1wVREUk6GVGMGQQedvfdZpYL7DKzZ91938WDzCwd+Dbw9DjUOWYaWjuZkJHG8rJpsS5FRGRMjXiG7u4Bd98dedwDNAOlwwz9CvAU0D6mFY6xhtYgt5Tnk52ZHutSRETG1Kjm0M2sHFgG7Ljk+VLg08BjI7x/k5ntNLOdHR0do6t0DAR7L7C/rUfTLSKSlKIOdDPLIXwG/pC7n7nk5b8GvuHuQ1c6hrtvdvcad68pLCwcfbXXqKE13C5XF0RFJBlFM4eOmWUSDvMn3X3LMENqgH+N7PpTAGwws0F3/+mYVToGGlqCTMnOoKo0L9aliIiMuRED3cIp/TjQ7O6PDjfG3edcNP4HwM/jLcwhvP78torppKdpuzkRST7RnKGvAh4AGs1sT+S5R4AyAHe/4rx5vDjadZajXef44u0VsS5FRGRcjBjo7v4qEPUprbt//loKGi/1LeF2uat0QVREklTK3Cla39rJjNwJVBbmxLoUEZFxkRKB7u5sbw1SWzmdyIVbEZGkkxKB/s7JHoK9/dTO1XJFEUleKRHo9S2R9ecKdBFJYikR6A0tQcqnT6J06sRYlyIiMm6SPtAHh0LsONil6RYRSXpJH+hvHeum98KgbvcXkaSX9IHeEFl/vrJS689FJLklfaDXtwZZXDKF/MlZsS5FRGRcJXWgn+sfYvfh07o7VERSQlIH+q7Dp+gfCumCqIikhKQO9PrWIBlpxory/FiXIiIy7pI60Btagiwrm8rkCVG1fRcRSWhJG+jd5wZoPN5NrZYrikiKSNpAf+1AJyHX7f4ikjqSNtAbWoJMzExn6aypsS5FROS6SNpAr2/tZMWcfLIykvaPKCLyEUmZdifPnKelvVfrz0UkpSRloDe0hm/31wVREUklSRno9S2dTJ2UyeKSKbEuRUTkuhkx0M1slpm9aGbNZrbXzL46zJj7zextM9tjZjvN7PbxKXdk7k5DS5CVFdNJS9N2cyKSOqI5Qx8EHnb3RcBtwJfNbPElY54HbnL3pcDvAN8f2zKjd6jzLCe6z+t2fxFJOSMGursH3H135HEP0AyUXjKm19098tvJgBMj9ZF2uavULldEUsyo5tDNrBxYBuwY5rVPm9l+4D8In6UP9/5NkSmZnR0dHaOvNgoNrUFK8rKZUzB5XI4vIhKvog50M8sBngIecvczl77u7j9x94XAp4BvDncMd9/s7jXuXlNYWHi1NV9WKORsb+2ktrIAM82fi0hqiSrQzSyTcJg/6e5brjTW3V8GKs3suk9i7wuc4dTZAa0/F5GUFM0qFwMeB5rd/dHLjJkbGYeZLQeygM6xLDQa768/V/8WEUlF0fSVXQU8ADSa2Z7Ic48AZQDu/hjwq8B/NrMB4BzwGxddJL1u6ls6qSycTNGU7Ov90SIiMTdioLv7q8AVJ6Td/dvAt8eqqKvRPxjijUNd/NrNM2NZhohIzCTNnaJvHTvN2f4h3e4vIikraQK9viVImsHKCl0QFZHUlDSB3tDSSVVpHnmTMmNdiohITCRFoJ/tH+TNo6c03SIiKS0pAv31g10MDLnWn4tISkuKQG9o7SQrPY2a2fmxLkVEJGaSItDrW4Isnz2ViVnpsS5FRCRmEj7QT/X1sy9whlWaPxeRFJfwgb79QCfuqP+5iKS8hA/0+pYgORMyuGlmXqxLERGJqYQP9IbWTlbMyScjPeH/KCIi1yShU/DE6XMcDPZRq92JREQSO9A/2G5O8+ciIokd6A2tnUyfnMWCotxYlyIiEnMJG+juTn1LkJWV00lL03ZzIiIJG+itHb2091zQdIuISETCBnp9S3iHO91QJCISlrCB3tAaZOa0iZRNnxTrUkRE4kJCBvpQyNne2qmzcxGRiyRkoO890c2Z84PUql2uiMgHEjLQ358/14YWIiIfGjHQzWyWmb1oZs1mttfMvjrMmM+a2duRXw1mdtP4lBvW0BpkQVEuhbkTxvNjREQSSjRn6IPAw+6+CLgN+LKZLb5kzEFgtbvfCHwT2Dy2ZX7owuAQbxzq0nSLiMglMkYa4O4BIBB53GNmzUApsO+iMQ0XveU1YOYY1/mB3YdPc34gpAuiIiKXGNUcupmVA8uAHVcY9gWg7jLv32RmO81sZ0dHx2g++gOZ6cadCwq5tULbzYmIXMzcPbqBZjnAS8BfuPuWy4y5C/gecLu7d17peDU1Nb5z585RlisiktrMbJe71wz32ohTLpEDZAJPAU9eIcxvBL4PrB8pzEVEZOxFs8rFgMeBZnd/9DJjyoAtwAPu/u7YligiItGI5gx9FfAA0GhmeyLPPQKUAbj7Y8CfAdOB74Xzn8HL/UggIiLjI5pVLq8CV+xP6+5fBL44VkWJiMjoJeSdoiIi8ssU6CIiSUKBLiKSJBToIiJJIuobi8b8g806gKTV0IQAAALASURBVMNX+fYCIDiG5SQ6fR8fpe/jQ/ouPioZvo/Z7l443AsxC/RrYWY7tSzyQ/o+Pkrfx4f0XXxUsn8fmnIREUkSCnQRkSSRqIE+bv3WE5S+j4/S9/EhfRcfldTfR0LOoYuIyC9L1DN0ERG5hAJdRCRJJFygm9k6M3vHzFrM7E9iXU8sRbOBd6oxs3Qze9PMfh7rWmLNzKaa2Y/NbH/k78jKWNcUK2b2h5F/I01m9kMzy451TeMhoQLdzNKBvwXWA4uB3xxmw+pUEs0G3qnmq0BzrIuIE98Ftrn7QuAmUvR7MbNS4L8CNe5eBaQDn4ltVeMjoQIdWAG0uPsBd+8H/hW4P8Y1xYy7B9x9d+RxD+F/sKWxrSp2zGwm8HHCO2elNDObAtxBeHMa3L3f3U/HtqqYygAmmlkGMAk4EeN6xkWiBXopcPSi3x8jhQPsYlFu4J3s/hr4YyAU60LiQAXQAfxjZArq+2Y2OdZFxYK7Hwf+CjgCBIBud38mtlWNj0QL9OE22kj5dZeRDbyfAh5y9zOxricWzOwTQLu774p1LXEiA1gO/J27LwP6gJS85mRm0wj/JD8HuAGYbGafi21V4yPRAv0YMOui388kSX90ilY0G3iniFXAJ83sEOGpuDVm9i+xLSmmjgHH3P39n9h+TDjgU9HdwEF373D3AcL7H9fGuKZxkWiB/gYwz8zmmFkW4QsbP4txTTETzQbeqcLd/9TdZ7p7OeG/Fy+4e1KehUXD3duAo2a2IPLUWmBfDEuKpSPAbWY2KfJvZi1JeoE4mk2i44a7D5rZHwBPE75S/Q/uvjfGZcXSsBt4u/vWGNYk8eMrwJORk58DwIMxricm3H2Hmf0Y2E14ZdibJGkLAN36LyKSJBJtykVERC5DgS4ikiQU6CIiSUKBLiKSJBToIiJJQoEuIpIkFOgiIkni/wOVmvhpzDzQiQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(y_new[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "Compute the SSE value on the predicted values from 2015 onwards. Remember to retrain your model after doing a train/test split before you evaluate!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions, see `solutions/arima_sse.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seasonal ARIMA\n",
    "\n",
    "Seasonal variations are not included by default in an ARIMA model. You can approximate this with lag periods of your season, for instance, setting $d=12$. A Seasonal ARIMA model has parameters for both the standard ARIMA model, and additional ones for the seasonality. This gives the full model as: \n",
    "\n",
    "$ARIMA(p, d, q)(P, D, Q)$\n",
    "\n",
    "Alternatively you may see this model written as $ARIMA(p, d, q)x(P, D, Q)$.\n",
    "\n",
    "Where $p$, $d$ and $q$ are as with ARIMA, and $P$, $D$, and $Q$ are the same as their lowercase version, except with a seasonal lag added in. For example,  D is seasonal lag, setting $y_t = y_t - y_{t-M}$ where $M$ is the seasonal lag period (if you have monthly data and want yearly differencing, $M=12$).\n",
    "\n",
    "As with ARIMA, and all the model we have seen so far, any of these values could be zero, effectively turning that part of the model off. Unlike non-seasonal ARIMA, we normally still refer to the model by it's full name, i.e. ARIMA(0, 0, 1)(1, 1, 0), (as opposed to, say, MASAR(1, 1) or some weird combination).\n",
    "\n",
    "In statsmodels, the name of this model is SARIMAX - Seasonal ARIMA, with eXogenous regressors (additional independent variables)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.statespace.sarimax import SARIMAX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\marin\\Anaconda3\\lib\\site-packages\\statsmodels\\tsa\\base\\tsa_model.py:165: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.\n",
      "  % freq, ValueWarning)\n"
     ]
    }
   ],
   "source": [
    "seasonal_model = SARIMAX(changes, order=(1, 1, 1), seasonal_order=(0, 1, 0, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = seasonal_model.fit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Statespace Model Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>               <td>Value</td>             <th>  No. Observations:  </th>    <td>322</td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>           <td>SARIMAX(1, 1, 1)x(0, 1, 0, 4)</td> <th>  Log Likelihood     </th>  <td>891.773</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>                  <td>Mon, 13 May 2019</td>        <th>  AIC                </th> <td>-1777.546</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                      <td>11:14:20</td>            <th>  BIC                </th> <td>-1766.269</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Sample:</th>                   <td>07-01-1949</td>           <th>  HQIC               </th> <td>-1773.041</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th></th>                         <td>- 10-01-2029</td>          <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>              <td>opg</td>              <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "     <td></td>       <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ar.L1</th>  <td>   -0.3695</td> <td>    0.020</td> <td>  -18.912</td> <td> 0.000</td> <td>   -0.408</td> <td>   -0.331</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>ma.L1</th>  <td>   -0.9999</td> <td>    1.169</td> <td>   -0.856</td> <td> 0.392</td> <td>   -3.290</td> <td>    1.291</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>sigma2</th> <td>    0.0002</td> <td>    0.000</td> <td>    0.850</td> <td> 0.395</td> <td>   -0.000</td> <td>    0.001</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Ljung-Box (Q):</th>          <td>196.94</td> <th>  Jarque-Bera (JB):  </th> <td>7969.37</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Q):</th>                 <td>0.00</td>  <th>  Prob(JB):          </th>  <td>0.00</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Heteroskedasticity (H):</th> <td>713.28</td> <th>  Skew:              </th>  <td>0.66</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(H) (two-sided):</th>     <td>0.00</td>  <th>  Kurtosis:          </th>  <td>27.53</td> \n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Covariance matrix calculated using the outer product of gradients (complex-step)."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                                 Statespace Model Results                                \n",
       "=========================================================================================\n",
       "Dep. Variable:                             Value   No. Observations:                  322\n",
       "Model:             SARIMAX(1, 1, 1)x(0, 1, 0, 4)   Log Likelihood                 891.773\n",
       "Date:                           Mon, 13 May 2019   AIC                          -1777.546\n",
       "Time:                                   11:14:20   BIC                          -1766.269\n",
       "Sample:                               07-01-1949   HQIC                         -1773.041\n",
       "                                    - 10-01-2029                                         \n",
       "Covariance Type:                             opg                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "ar.L1         -0.3695      0.020    -18.912      0.000      -0.408      -0.331\n",
       "ma.L1         -0.9999      1.169     -0.856      0.392      -3.290       1.291\n",
       "sigma2         0.0002      0.000      0.850      0.395      -0.000       0.001\n",
       "===================================================================================\n",
       "Ljung-Box (Q):                      196.94   Jarque-Bera (JB):              7969.37\n",
       "Prob(Q):                              0.00   Prob(JB):                         0.00\n",
       "Heteroskedasticity (H):             713.28   Skew:                             0.66\n",
       "Prob(H) (two-sided):                  0.00   Kurtosis:                        27.53\n",
       "===================================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n",
       "\"\"\""
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "1. Check the documentation for SARIMAX on the statsmodels website. These values were set in the previous code example?\n",
    "2. Choose a seasonal commodity from Quandl, such as Wheat, and apply a Seasonal ARIMA to the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For solutions see `solutions/arima_seasonal.py`*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extended Exercise Automating parameter selection\n",
    "\n",
    "As we saw above, choosing parameters for ARIMA is a fairly straight-forward process, using AIC to choose the best from a subset of features. This process is easily automatable, and has been automated in a few libraries. In the R programming language, the ARIMA implementation already includes it, while it is a separate module for statsmodels in Python. You can get the code from https://github.com/tgsmith61591/pmdarima\n",
    "\n",
    "Install the package on your system and run on the data you received for the exercises in this module (cryptocurrency and seasonal commodity). What parameters does it choose, and how effective was the algorithm?\n",
    "\n",
    "As a warning, automated parameter selection is basically a brute force selection. Generally it composes the following steps:\n",
    "\n",
    "1. Try all parameter combinations from those given\n",
    "2. Evaluate all of them using some metric\n",
    "3. Choose the best one.\n",
    "\n",
    "Some algorithms exist separately to optimise and improve this process by being a bit more clever about their search, but these options are not available in `pmdarima`.\n",
    "\n",
    "Further, the documenation for this library also contains additional information on choosing parameter bounds that is worth reading for a more complete picture of ARIMA."
   ]
  }
 ],
 "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

Copy Edit Raw Blame History

authored 2024-02-12 21:11 . dharshan

马建仓 AI 助手

尝试更多

代码解读

代码找茬

代码优化

https://gitee.com/sinomiko/QuantFinance.git

git@gitee.com:sinomiko/QuantFinance.git

sinomiko

QuantFinance

master

sinomiko/QuantFinance

About

Releases

Contributors

Activities