FitFullQuantiles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tssimul\n",
    "import tspred\n",
    "import tspred_qtl\n",
    "import rnn_analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import (mean_squared_error,\n",
    "                             log_loss,\n",
    "                             confusion_matrix)\n",
    "from sklearn.linear_model import (LinearRegression,\n",
    "                                  LogisticRegression,\n",
    "                                  SGDClassifier,\n",
    "                                  SGDRegressor)\n",
    "from sklearn.decomposition import FastICA, PCA\n",
    "from sklearn.ensemble import (RandomForestRegressor,\n",
    "                              RandomForestClassifier)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "burn_in = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AR(1) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "ar = [0.0868, 0.3667]\n",
    "ma = [-0.1150, -0.4068]\n",
    "sigma = .0112"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation = tssimul.simulate_arima(ar, ma, sigma, (5000, 1000), 'n')\n",
    "simulation_test = tssimul.simulate_arima(ar, ma, sigma, (5000, 1000), 'n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 3, 1000), (5000, 3, 1000))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simulation.shape, simulation_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_t = tssimul.simulate_arima(ar, ma, sigma, (5000, 1000), 't', 4.5)\n",
    "simulation_t_test = tssimul.simulate_arima(ar, ma, sigma, (5000, 1000), 't', 4.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_exp = tssimul.simulate_arima(ar, ma, sigma, (5000, 1000), 'exp')\n",
    "simulation_exp_test = tssimul.simulate_arima(ar, ma, sigma, (5000, 1000), 'exp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.savez_compressed('simulation', data=simulation)\n",
    "# np.savez_compressed('simulation_test', data=simulation_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = simulation_exp\n",
    "data_test = simulation_exp_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Volatility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv_beta = .98\n",
    "sv_sigma = .168\n",
    "sv_intercept = -3.926"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv, svvol = tssimul.simulate_sv(sv_beta, sv_sigma,\n",
    "                                sv_intercept, (5000, 1000))\n",
    "sv_test, svvol_test = tssimul.simulate_sv(sv_beta, sv_sigma,\n",
    "                                          sv_intercept, (500, 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DescribeResult(nobs=5000000, minmax=(-0.16251504124838245, 0.16889895057080634), mean=6.3666536513315898e-06, variance=0.00010975681166093846, skewness=0.008468034416795452, kurtosis=3.1227177163589674)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.describe(sv.ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DescribeResult(nobs=5000000, minmax=(0.018045324639038317, 1.1020478390646546), mean=0.15329274226257966, variance=0.0044917155658833405, skewness=1.4174405203109106, kurtosis=3.7752985793496885)"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.describe(svvol.ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DescribeResult(nobs=5000000, minmax=(0.00032563374132828335, 1.214509439587075), mean=0.027990379497921915, variance=0.00080035720864075953, skewness=4.054450387095615, kurtosis=36.16007596911728)"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.describe(svvol.ravel() ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DescribeResult(nobs=500000, minmax=(0.00035299657692129102, 0.89837317286388985), mean=0.027837523641173794, variance=0.00081043106975812117, skewness=4.255823225893204, kurtosis=39.43058636697476)"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.describe(svvol_test.ravel() ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = tssimul.simulate_arima_given_innov(ar, [], sv)\n",
    "data_test = tssimul.simulate_arima_given_innov(ar, [], sv_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.savez_compressed(<file>, data=data, svvol=svvol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regime-Switching Process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs_p0 = .8\n",
    "rs_p00 = .9\n",
    "rs_p10 = .2\n",
    "rs_mu = [0, -.005]\n",
    "rs_sigma = [.08, .3]\n",
    "rs, rs_prob0, rs_regime = tssimul.simulate_rs(rs_p0, rs_p00,\n",
    "                                              rs_p10, rs_mu,\n",
    "                                              rs_sigma, (5000, 1000))\n",
    "rs_test, rs_prob0_test, rs_regime_test = tssimul.simulate_rs(rs_p0,\n",
    "                                                             rs_p00,\n",
    "                                                             rs_p10,\n",
    "                                                             rs_mu,\n",
    "                                                             rs_sigma,\n",
    "                                                             (500, 1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 1000), (5000, 1000), (5000, 1000))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rs.shape, rs_prob0.shape, rs_regime.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DescribeResult(nobs=5000000, minmax=(-0.10027182433879982, 0.1031418152339469), mean=-0.0016676499245161524, variance=0.00013972994882308483, skewness=-0.6644921381648347, kurtosis=3.8469821872161205)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats.describe(rs.ravel())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.011789826122551595"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(0.000139)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12b1d9c18>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXe8FcX5/z/P7XDpvQkXREEQC17BhhWRiBF7NIkhRr8mUX8xMYlirFGjxlgSE01C1K9+1VhijyKIKDZEAUU6cmkC0gSkXS63nPn9cXbO2d0z23fPnvK8Xy+4Z3dnd2d3ZueZp8wMCSHAMAzDMJKSuDPAMAzD5BYsGBiGYRgDLBgYhmEYAywYGIZhGAMsGBiGYRgDLBgYhmEYAywYGIZhGAMsGBiGYRgDLBgYhmEYA2VxZ8APXbp0ETU1NXFng2EYJq+YO3fuN0KIrk7p8lIw1NTUYM6cOXFng2EYJq8gojVu0rEpiWEYhjHAgoFhGIYxwIKBYRiGMcCCgWEYhjHAgoFhGIYxwIKBYRiGMcCCgWEYhjHAgoHJaxqaWtDUkog7GwyTYsfeJrz2xddxZyMQLBiYvGbwTVNw+l8+iDsbDJPi18/Pwy+e+Rwrt+yOOyu+YcHA5D3LN+fvB8gUHuu27wUA7GvOX02WBQPDMAxjgAUDwzBMiAiR/EsUbz6CwIKBYRgmRASSkoGQv5KBBQPDMEwEsMbAMAzDANCZkuLNRiBYMDAMw4SIJhdYY2AYhmEKBxYMDMMwISJESmeINR9BYMGQJyz6eoeuwjEMk6uwKUmDiMYS0TIiqiOiiYrjlUT0nHb8EyKqMR3vS0S7ieg3YeSn0Hhn6SaMe/BDPD9nbdxZYRjGiQLovwUWDERUCuAhAN8BMATARUQ0xJTsUgDbhRADATwA4I+m4/cDeDNoXgqVlVv2AACWbeSpHxgm18l/Q1I4GsMIAHVCiJVCiEYAzwIYb0ozHsAT2u8XAJxClFS0iOgsAKsALAohLwzDMLEiTb6Ux7akMARDbwB6G8c6bZ8yjRCiGcAOAJ2JqA2A6wD8PoR8MAzD5Az5Kxbidz7fCuABIYSjjYSILieiOUQ0Z8uWLdHnjGFCoKklgbcWbeTAgSKiEEo6DMGwHsB+uu0+2j5lGiIqA9AewFYAIwHcQ0SrAfwSwO+I6CrVTYQQk4QQtUKI2q5du4aQbYaJngenL8flT87FjC+5M1MsFMIkemUhXGM2gAOIqD+SAuBCAN83pXkNwAQAHwM4D8A7ItmFGiUTENGtAHYLIf4WQp4YJieQc/Nv39MYc06YbJPPk+gFFgxCiGatlz8VQCmAx4QQi4joNgBzhBCvAXgUwJNEVAdgG5LCg2EKHjYhFR+iAIxJYWgMEEJMBjDZtO9m3e8GAOc7XOPWMPLCMLlIPpsVGG/IvkA+C4i4nc+MB7hxyV9YcSgeCqGsWTAwTITkcyw7E4x8FhAsGBiGYSIgj+UCC4ZcZfOuBnyx9tu4s8EwjEcKIeCABUOOMuaB9zH+oY/izgbDMB6RYiGfBQQLBpcs37QLh9/2FjbtbMjK/b6tb8rKfZhoyefGgfHOll37sGFHdtqIKGHB4JLHZ67G9vomvLV4U9xZYfIQ9kEXB8ff827qdz53CVgweMWhB/ifOWvxvX9+nKXMMMXGuu316H/9G1j89c64s8Io2NvUEncWQoEFg0vc9vh++8J8fLJqW7SZYYqWaYs3QQjgudlfxZ0VxoF8tiKyYPBIHpc14xIhBN5csAHNLYm4s2IJj4/IB/K3tWDB4JJ8nhArWzw3+yu8Os88sW7+8ebCjfj505/hn++vjDsrDBMLocyVxGSHXFdNr3txAQBg/GHmdZryi6279wEANuzYG/haYRdZrtcBJk0+lxVrDExO0tSSwN7GeB15ufhh52CWmAKEBUOWeHXe+oIws2SLCY99ioNunhLPzUO030dlgGQXQ+6Tz0KcBYNH/PYir352Hq5+dl6gexdTYzBzxda4s5DXH3YYrNm6Bze/uhCJhPObqNu8G8s37cpCrphswD4GlxRTo1zs5HJRZ3Mk9f975nPMX7cD2/Y04tzhfXDS4G6WaUff/x4AYPXd47KVvZwnF02RbmGNgWEsyOUPOxtRcgntBbw+fwMueXx25PcrNHihHoYpIEpS6mHuftiswTJRwoLBIzwpWuEjG90wipprS/GSz00FCwaX5H4fkgmbMD9sHiDJ5BMsGFzCUxAUD7lc0lJY+c3jATdMxsMz6kLLD2ONl47F2m31mLPaOMdaIiFis1CwYGBcs3FHA/Y1F8bskW7IReehzJPffkpTi8A9U5aFmCMmDEbd8y7O+0d6VuaGphYM+N1kPDDty1jyw4IhB9i4owFPzVqjPJYrPo2WhMBRd03HNc99EXdWIoeVQyYMgnQs9uxrBgA8adEuRA2PY/BIFO30JY/PxpINOzFmaHd0a1sV+f380JxIzjQ6rYgWKsqVd89kh29270O7qnJUlIXTXw5Sf6TpOq4qyBpDDrBtT3LSthbFCNNcaZtSlbwIetPSURzmuw/LLJXyMbBaEzq1d7yNXzzzedzZAACUaMXrZtR5JPeP5a5FgBcTkF3ESi6ZkoB0hS1ocvgZg8jnXKlLuYh8N1MWbYw5J0mi6Jx4gQVDDrJm657U71z5lFuEFAw53GqGRCo0ucDCVWPqfOYFOSczY46PD0UwENFYIlpGRHVENFFxvJKIntOOf0JENdr+U4loLhEt0P6eHEZ+cgE/FU2ec8KfZhj25UKlTSTiEQxx9HLT9t0cePEhwhqDNS0RvJtAlxSGP1knsGAgolIADwH4DoAhAC4ioiGmZJcC2C6EGAjgAQB/1PZ/A+C7QohhACYAeDJofqIiNRrWZXovBWrX1uobJ5msuSWBZz79SumTiIKVW3bjxc/WG/KQLWJty3KwHQ3yPnLwcXKGRI4Jzbg7JWFoDCMA1AkhVgohGgE8C2C8Kc14AE9ov18AcAoRkRDicyHE19r+RQBaEVFlCHkKnZTNz2UFmrVSPW30gnU7PN1XdbvHZ67G9S8twL8/yU4o29i/fIDbX18MIPuhnHF8sH4esW7zblz/0vzIhXWqwfCRSa/vMsfaykiJ4lmDNO6yGuXzALfeANbqttdp+5RphBDNAHYA6GxKcy6Az4QQ+0LIU+z84JFPlPu/+7cPLc9RVQFVvfi2vsnwN2oamxOp3yVZ9j7ni8Jw5dOf4ZlP12L55uysSeDHZ1FMDb1XpEAPs+MTSLvTTo6ryHJiHAMRDUXSvDTGJs3lAC4HgL59+2YpZ9kh7exUhataV404Kk22fQx2vdyoelPpSfTcX9+qnMLOYrDGxlt6t0X9yuf5vzKhrGdeandDUwuIgMqy0tDzI4sqLmEehsawHsB+uu0+2j5lGiIqA9AewFZtuw+AlwH8SAixwuomQohJQohaIURt165dQ8h2fpBrvbxc8jFE9W78yD6p+lsJzijk6cL1O9DQ5H6Kkqjs1r98LtjKhGFQt3lXoBXkZPl5GR8y+KYpOP6edy2PB/I9p5zP+WtKmg3gACLqT0QVAC4E8JopzWtIOpcB4DwA7wghBBF1APAGgIlCiI9CyEvBoaoWYU4L7ZVsD6yyFQxR39tLWh89ziBs3tmAM/76IW58ZaHrc3KtkxEmo+9/H6c+8L7v8xM+x+ls2pm0fAsh8MgHK33f34wUCA1NCYeU0RBYMGg+g6sATAWwBMDzQohFRHQbEZ2pJXsUQGciqgNwDQAZ0noVgIEAbiaiedo/6/UDYyTK9jAVHqn4cFXmjDgj4rM9wC0e57MP+708N0vv59u9Sf/SwvXugxmCvMtCH9iYNiX5e9A5a7bjjjeWGPYFMnXqTl3/7V7/1/FJKD4GIcRkAJNN+27W/W4AcL7ivDsA3BFGHgoB1YdrV7XiUDOLwceQvr6XxPJHtO9HPnOz1sMt9dBiB3lbVveJa8qGsEk9hs/iqyjN7GOHMIwBQDI8PdvwyOccQq0xKBLaaBhRk/1wVetjUT2+1zErqnPNhF1WskEu8yIYAuTBSjBEMTAsDmQHxK9m1LoiXAd03K+VBYOJ3fuabWPRoywwZe9YsSvO0fJZtyjE4HxOX99LVFKSHXubcN0L81HfmJw2OS1kwsmszJKc7dZL+HAQDausRN1UZGuQZdQENSW5DTV3fz39wNbs2/FYMOgQQuDgW6biuhfnx3J/1TcWpblIzvnuhWw7n21NSRG9Gz9THstG94FpX+K5OWvx4tx12n4Y/oaFJheypjFY3aa5YARD8q/f6h2WLyw1fkF3uTimJ2PBoEP2fl6Yuw6ffbXdcCzdS1dXgJl13+C/X3ytPOYehY/Bpr6t21bv+04L1+/A0FumYvKCDZ7O45HPamQupbA1z+kf1mPIy0iNwYuPIci7LFPY0IH4NIbJCzZg086GjP2bFfvsEELgqVlrsH1PIwD/GrH61Xp/N/I6+rKKw6zEgkGHvvdzzsMzDcecGsTvP/IJ/p/PudzltdUaQ1oYfVj3jSH9SwEGFi3Qolne/3KLp/OyPome3bEcileVealvTI4rKCspwcYd6UYq7Ky2aBcst2iwVUThfI5DMDQ2J3DF05/hwkmzMo6NuHO6p87OvLXf4sZXFuKMvyZnJPCrEYfVgUkoNIY4OkcsGHSErRa/Om893l262TFdWjCoNIb0vqUbd2HGss2h2Bz91rWiCFf18YxSeEvBcOfkJTjqrukpDSLsCKoWHxpDIOezxUuRmks2kXXCKoxz7prtyv0qZHlJwtQYfM2wrNgXh4OfBYOOlhbnAtCX0bTFm/DNbuupna5+dh4ueXy260ZB9Y2Zz1ywbgceeDueBcIBYPXWenz+lfsPLyhWr25XQxP+o9nxI7u3jz62bGi2aqYJuR3Wp51yPrf4iUrynwuvGoMQAm/M34B9ze5HZjshhEBjcyJwBJGejI6H7ppPzVrjetCay7gR1/kxaAwxaGUsGHQ02fR+zCrmjr1N+J//m4PL/2+O43Xto5wEdu7VepUufAyzVqlnbfWLn97xnZOXOCfyiRACz376lW5bne7+aV/iJg+jfq3Qm3skUiNbuWVPxjErZD73Nhod+qlgntB8DMkLyQbEi8awp9F/I211n2ZTZ0oKn5krtuLKf3+Ge6Ys831PM/dP+xIH3vgmdmtamJXm7EX+mb9Nvan0xlcWZgxasyI857P2V1dh4nDjsGDQ4cVeuqshOfJU1bCYsbvsk7PWYIc2ilXd6zDuDEtzDxLR07F1RTiZULDo652Y+NKC1HaUpqQX567DUXdNx9w125THl250P/eOzKa58U0v0ejuOYQQePLj1djZYD9zbnNqHIP7T/jsh73NOqN/9Vaaibl8pKa0vT7518334ZZnPk1O4rxT+17CcHdlKAw+rxlWuKoqii0OPw4LBh1NHkYYShNB60rnweN2jdu0xZtSv9dtV9hMhXkzukpSt3k3PtIc3HZ0bpMWDGu27sGr88KbXdNserB6dz3aVTlea8++ZlvzyRxNICzbuNuw3887trJ3e53Xas6a7bjp1UW44WV7bajFx8hnp2naV27Zjd37mlG3eRdqJr6BRV/vTB2zGi9h9svV3vE2gOAjie2Qn2kYlzY3un6vGbrzOYJre4EFgw4vklk6FatdCIbBN02xPKY3Uf3sqbkZx805MmdxR30TfvDILGzY4Xc+lfT9R9//nuU6Enra6J55/EMf4epn54XoXDV+mlaXdWoQ126rx9BbpuLpT76yTRcGbmY4dft2ZIfj2/pGrP92L1ZsMQkts4+hNLyW9+T73sOExz7FO4qACSuNQfXNCCEinVRQ3jOMCDlzo+v3mqp6esE/P/adH/33xBpDzLiJSpIpUhpDebCh8E7VsLE5ga/04xVMWXzp83X4qG4r/jHDcsZyJWG147IXuq/Zn43r4xVbMfz2aSnTnPm7tMqnU1ihbFDfWrwJW20CBIBMDcHru3HTo3N7TX2DcM7DH+GU+94zaLLClM4qWsgr0sE5d812ZeNoFsT7mluwt7Elw8cAGDsv4Q6INPpXwpA65rLzbUoK8EHt0pkNhekvwFFJsXL764vxr/etIxDM9SWtMQQUDA4V8cZXFuKpWeler7kip74R04X2NrZkRDMsWLcDay0GxUk/hxv0WSjXeqx+RlEDwH1vLcO2PY34qG4rhBAZ79mvGi1Pe//LLTjijrdTS6r+v2c+x9RFG7VUmXNOvT7/a7zscXxIk5toNrc+BpkzotSUzvMVy8E2aWW7bU+jp7KzolEnfFSNuVkwnHzvezjo5inKMQPNiQRCbLsBJHvN6elAkj92NTTj5c8zI9O8mALNfUEhgL/PWOH5nbrt1Nc3NuPw297CjGVJrWzttnoMu/Wt9P0T6Xykrs0aQ3ao27wrY7riRz9chWdnr7U4IxPZQw66etOMZfYDzN4zDUAzV5GMbSGwt7EFB908JSN66Lt/+xCjLBYW+Z4HtVd/zyrt+c3x4GYamlqw6pvMKB/Z8P/sqbl45tO1GY2SSjDs2dfs+LGYGwe55OZ/v/gaP30yabKTt6rbvBurtbxd9e/PlaYUO9zMreX629Y1qIN7tAUALNmwMyNZg/a+py/djJPuneHy4kYSCZGK+d/XlB4XobIamQWD9Kn87d06xXXT73/qoo14fb73GQGaWxIpE92SDTux/+8mpxzbeg3qV8994fnaesxlt3VPI/44ZSkemOY+JDwptNwV8Mote7C9vgl/mpqM1jJ31NL1naOSss7o+99PjXT0iiy3lJ3TRzD15l3+IzXMFTBly9Wy8cvn5uGgm5M+jefnuBd0XiJw9FSWJ6vQnkZ7jeHX//kCJ907Aw1NLajbvAs3vLwAiYSx+X7/yy0ZPUzzN7Ho6+RUHv+1aGw272rA/HXfZphuSksoY/piea/HZ67GiRaN65Cbp2DKQvuRtHaDvFJP6LLhkOmTS0Ym362sa+Me/AAPTl8OAKjX+TW2aQ2mV/7+3gqc+/eZ+HTVNuxrSV6vvJTUpiQPNha9xrCvOYGr/m0/I8A9U5bifz9aZdh30b9mpXxz5rDhJgezpRcl00ojtRsfsnX3PkMnpzmRcN14J0zfq1ng/uLZz7V06X3sY8gyX27a5T6ixjRLZnPKAeb9vmMCrDSV6YyWTr5kRl6dl24wo5rw7tEPV+Gqf38GIK0x7dlnrzF8oGk+extb8LOnPsPTn3yFVVv3ZEwWluljMD7xovXJ3rPKvAIAY//8Ac7820cZH2p5aYnBXOKW+sYW/HHKMjS1JFAz8Q08NWtN6tj/fbwaNRPfQEOjjWDQNZBuMJhgtJchGwZ9lJCXxkIIgYcUPfvF2vU272pIaQwVpSXKzo6X6KemFoEtu+z9OnoenrECv//vYsO+2avTgyjLTQ52cznOW/ut7fVnrdyKmRbRdlaCoWvbSsvrnXjvDIOW1twiHM1Xi7/eibXb6jOWgDW/1w+WJ/PJU2JkmR26kL0xD7yPq5/1tl6tNGGkpiTw0fhKh60f26G5jkgBpcqG3Pfel1uUtlir8yST3l+BU+9/T3ns9fnJXnSVpjHUO2gM8gO4Z+pS1G1OOob10StA8mMxv5KMV+TwumXv2SxQykoo1filLuWh6GTs/H1vpQds3a+ZG77da91jl9mQA6WEELj7zaWo26zW0N7Q2exlm+GlYZj0/gp8uDzdCM5f9y3eX/5NynShR5pkykooJbgqykpcmZLsePjdOtz15lLlsZaEMGhuu3W+KflO9urMkuu212dos2aH91kP2Y/PuHDSLHxfi7ZLJAT+8MZirNm6R9tWn9PKZn2FXQ3Gut7cIhy1lNMf/ACj7nlXp+ET3pi/Af+08GvqBQ1rDFng3re8jcTcsmsfaia+gde/2KCd/yU27mjQNciEmSucY/9VqKIN7KbYADIbPPmRqD5buW/CY596ssX+VTNX3Dl5KZZv3m2btk1VOQBgu0OMvGxY5CAlIBmPbngayvwIbGYsMDD2z+8bhL75zZYQ4R/vJSO3pEBwO+eUECLVS9XPmirzatco6Bv1Jz9ejSUbduEf762wNK+89Nn6VH5TM/p6aBfunLwUP3w0HXJ85t8+woTHPlWm1Q+Sa9QEQ3lpidKU5CWMc9qSTZbHzvjrhxh4w5sAkpFNB98yNXVs9P1JTVqf3+P++G5KAEu8jDcys2LLbvzrg1UpP5NVxI+Xd540Jbk7ITW9N4Ar//2ZYRyT1f05KikLeFXLpNNSP4Bp6cadhsXDv/8v59h/FaqewGkOZibzGfIjKSkhzFmtHsHr5jp6/vpOptnBig6tkoLBKSRU1bC8Mf9rQ3kQMsvnrUUb8W19I37y+Gx8s3ufZQO1dOMu3PZ62hxhLuY3F25M9c7KPYwWBpK9WhmlIgXDwvU7Uj1Hu3epP3bTq4twy2vJgWtVDmHOm3Y14LOvkiaSqBqGlMZQSqmBhRVlJbbics7qbVi+yd4fpcquHDind6Rbmdc+dajHTiZB2XlaunEnaia+YTym/W1OCXX1u/XSTjQnnDUGcw7s5GxyYGZ6e29jC/73o1VZ1RxCWfM5n/AyTTGgVjX37GvxteaunlfnrcfJg7tl7N/q4Eg0q9EpzQXAef8wRhZZ+RjqNu/GUu0DVaXw8lG0rUpWoa277fOtek8PvlOHob3apbaJKKPy3zftSyzZuBPvLN2MMx78ENeOHWR5jxc/S5vLzB+8nKIBSNusza9n/jq1rfqb3Y0Y++cPAKTX9tUHL9iuyWt6lWu3JTsYnartpxVZuD7dgEZlY9aPnt7noDFI04a5jrnlsidmY9bKdIN/4aSPcd3YwRnp3Kxv7CY8GAAmL9hoeUzWD6vb6avh24s3YfSQ7pbX2rJrH/7lcbI9Ow1sqE6LAoCH3q3Doq93ok1lGc6v3c/VfYJSdIKhssybYFD11q7UHK+Afwfv1c/Ow9EDOns+zxwB02wzP4BVzkbr/AZPf/IVbht/sOG4lx6qTLmt3l4wWH0IBucz1H4X6djeuLMBc1xOqWzuVepHJ1stOvPzpz5T7tejCk+2a6jMTskGrWfeXtO03OCmOMw945aEMPhDVMhOxlfb6lP5aVVequzNflS31XCPztUVjp0YPXqhILfPVwiZdx3Ct5P5dmdKkv4vO6zqeovuO7vs/+Zg9d3jLK/x6IerUk5jJ9LOZ1fJAaRDwbc4aOVhUnSmJJXGcObfrENXnXprQab+/Xil95lSrTWGzIy4lVnmqA63cuHxj1alpipvNJkFmloShn1W1huDKYnUH6o+dLDK5bgRs5lCLxhSGoPpHKv5jvSYV2YD7G3e5seRz+LFZu/HhPDu0s142GE0vJxN+IaXF+I1LZptSM92rvJm55xVUaH47lQzDeze5zywzK2PwVxXbnplYcb3Y54NV+LFjeFlQKSftaWl2dHvIFI/FJ3GoBIMVqGPgHPkkF9Tkl/MU4OnnM8BsuG3wt3638WpeZP0QkAIgaPvegd79jVjye1jAVhHb5k1BlUjqA+ffMwU727FPtP8RXpBIcvMj7anKm/bAW6mbTkbqv7Wjc0J3Dl5CX5xygHKa/gxJblpPPWN5GLNtFhSQpZCXI/dq1OFbpaXEtzM+t3Q5JzvRpemJLPwenLWGny6yqi5WI3fabEIV7Kaidct6QGP7stUaj5OIeFhUnyCocx9Y3DdC/Nx1P6dbNNke6lLaaOWSNOSOhfu8uY0atkOObBN3xBt3dOYEV1lNRBQ34gQkfKD8RMSvNckGPQaQ5D5hfbr2Cpjn10jbH6elFDS7ZuyaCMen7kaO/c2oX+X6owR4gkBR4evGTfNjj5vm1PjDoSrOm3X41W1qeVlJXAjGb508ZxOpiQB4PnZa5Wjl5dp15dPbjUtuFmbEUKAiHDu3735WGbWfWNYWEv6NtyaRIG0tuUUEh4mRWdKUqm0Vjw3Zy0e+cC+h5pluZCBncbgNm+qWV3vtohDNyPbFn3jeNwf30n9flsLx1O990Hd2xqcfC9/vl4ppPxoZeae5wrd6FkiwpsLNhji5d0yb+23GbOpbrIZzKUagQ0Yl5+UT7evOaGMkkkkhHIcgh1etQyp8SUS7jQpuyQqk5zb7+5/P1rtmMZqjISea1+crxN21lhpe39+e7lhW77Oti5mU9bz/Uc+MQzW8xNYJIVU2EsP21FUGsN5f5/pSVIDziNWYwgxNmDrYwhwXRnz7xb9e9I3yk9/sgajh3RXNjZNLYmMBmbO6szyqXMYS6HCbn6e9d/uxc+fdnY0q1i9tR7XvTjfsM/LSnJSQ1r5zR48P2ctavt1TPk8GppasHpr5iSHCSE8C0fzYD5lXhR1V0BEsq6312jAIHj5Jt2mbWxJYN039ejarhK7tvjvufsJPZYdmGxOpldUgqHBx/qzQQbTZAPp+FJ9zNnUZqze07vLtmDR1zuU6n9jSyLDmfv4zNUZ6ZZ5NKMAwJebvAsTt+inHXHC/CnrfTHXvpAUMI/8qBZAckI8FQnhfU4usynNLZt37XMMPQa8dzpUTvtcwG1D/eD05Xh4xorAPkUr34UduzSHvEvXSiiEUlpENJaIlhFRHRFNVByvJKLntOOfEFGN7tj12v5lRHRaGPmxonWFdzmomm9ej6ohiwVljzx7NcnuXuMe/FCpsm/d3eiqZ5vPmE1DKg3UabGdhGI6cif8mMmA5Gy/t7y2yDGdSrOxwzzfUZS4mnZbS+I24muBNhtz0EFm89ZaB7pYsSPAFDp+CSwYiKgUwEMAvgNgCICLiGiIKdmlALYLIQYCeADAH7VzhwC4EMBQAGMBPKxdLxLaeLQPArmvMdgRxjz9bpnrYKIzR1MByV6tmxDRfMbcIVUJQiebfotuemy3BAkoiIJsagy7G5xNPfuaE1i7rd5yUKOZVgEX5JLI2XG9sFN7nu31jfj6271ZERBhlNYIAHVCiJVCiEYAzwIYb0ozHsAT2u8XAJxCya9hPIBnhRD7hBCrANRp14sEN8twmskXwaBqWloSwnFR+WyxaYfREWg3rXG+0bayzPJ5zL1X85rWgHNP8OXP12ODRfSMFfpIGBWff7Xdl2nVL9n0MbziwtRX39iMUfe869oZ7KftCJuZK7bimLvf8b1aohfCKK3eAPQT/6+lorLvAAAgAElEQVTT9inTCCGaAewA0NnluQAAIrqciOYQ0ZwtW5xHR6po42O1tWyaY4Jg1em8Z4p1BMfAbm0iyk0m5pHIXkeg5zK79jXj6P3Vo9jNipKqIbrk8dm21/e73oIdZz88M2OdgyjxOj9V1HhdYMvrgL6oIHI3ojsouVVaNgghJgkhaoUQtV27dvV1DV8+Bh/OojiwiivXLwtqpjrGyl4egmA4aoD9GJNsIueMMpMf3YrocfKjZJv9OmWOR7HDqiPzo6P7Zey7bfxQX3lyQ3LKkujfZRiCYT0A/cxOfbR9yjREVAagPYCtLs8NDT/qoJuRmLmAn7riNMNnlMgRwKcNtZ6czIkxQ3qElR0AsBx57IayHOsR5xrZniHACbejpyVWgyJVnc0ozWZh+TqcCOMJZgM4gIj6E1EFks7k10xpXgMwQft9HoB3RDJc4zUAF2pRS/0BHABAPXl8CPgxJeUL5mrbpY397J0AUBmjYJBRKkFGjlv10v0SxLxldW5uNYfx4cen1LuDt169F5yWBzVj5YtQRVtFWebZMmkFFgyaz+AqAFMBLAHwvBBiERHdRkRnaskeBdCZiOoAXANgonbuIgDPA1gMYAqAK4UQkXnEcsGBFBXmOPf9uzr7D7yMAg8b2YP005NsU1mGBbeOCb0XGkQwHNGvo3J/QojQ3vO/tPEO+UipD43qjEN6RpCTJF6DSqxMyipNMcppclrni2AAACHEZCHEgUKI/YUQf9D23SyEeE373SCEOF8IMVAIMUIIsVJ37h+08wYJId4MIz9WVPvwMeQL5tDUAV2rHc9RNYTdbNa6DRNzD9KLv6NNZRnaVpWHLhjsQiqdzEzfO1I9T74Q4YRq/viYGow6oEvg68SFG42hj3keqojaVyI4rkxoxkqQKH0nEaoM+WRKyhvsNIaxQ8O1V2ebSaa1YztVV+DYgfbrPZSXEqb/+gTDvkP6tMdPTxgQev7MyEa9fatydGxdjnvOO9T1uTIqI2wnnJ1t2C6ktLKsxDIvLUKEEoFVVkKoKi/F4tvCHQN6UM/0QklRDkIrdXHt358ZndNWj9up2/X894sNGftm3zBaKfDYlJRnVNv4GMKIkvGD3QIgQVGp7/06t0793tecMJicLj6qH+674DB0qY5ea5AqeEVZCT6/eQzGeTAbyIV2gsySqsLO5GM3gZmd5pIQwlJj8CKA5TM7Obl7tq9yfU0AOGy/9qnfQZ2mJw1KRwua36UbjcF8fy9rFnjBT7jnbsXU9F3bVioXfYoyaihb40GKSjDYqWGFNOAKSH5UqmfSV6w92uhYOSL8J8f1R/tW5Y4RTi/+/GjTNb2/uyDfTnqxG//XUGHXObCb40YKqMd+XItJFx9hOCaEte/i8P06ZOzr36VaOYOn1XKkZq46eaB9AhN6QRPU5KXviPzp/EMs7wMAHVtnrmA3oGs1vj+yL8Yf1sv1PUf074SLRvT1lM8wzTGquh9lU5Kt6K6iEgz6RnH2DaMNx7K9rkLUEKkrkb4nJxssmcytyeOIfsbxA356MUF6VfJ+eod7lzYV+MPZB1ud4gq7zoGtxqA1DicP7o6R/Y3mu4QQlu9HpdGdX9tHOQGe2wbB60Cydq3SQihob1T//sYf1tvy2G9PG4RJCkd6RVkJ7jx7GA7s3tb1PXt3aOU55DnMMG1VuUTRlMi10bPVThWtYCBKOvQkcQToeC3j+853b4cH1A2d7BUP7tEWd549DEC6gZXvx2ujrb/PkJ7tlGGkFx9lHAjktXpfceL+6ftpDbHelHT8gV3Rrsr9OsoqzB/5r0YfiB+MTPZGm1uMJiG9T8rwnk0PlhDWC9+oyqeESCmEUmVj/wieFqICgH6d00EKQaOn7Aax6YV4+1blymf3c/+vttU7Cs1bv2ucus1LmPZ3Drb3PaoEcRQmMKnVZ8uwUWSCIf1WCcYoCD/hdEHx6pT06nhSTdV8VP9kb//ucw9BVy0CqWf75HtQrS5mhT5CRt9gXnpcf7x9zQkZ6W8/y9ibl6/b7fT0lx7XP/Vb9p6G9Eo7Tmv7dXLsTfVoV4WzD1fOuAIg08RwyXE1qXs0JxKYe2Nay9Q3Rvrf5leeXPhGfT9Vg1ZCwI3jDsrYLxtSJ6HtdaBdv05pn1NQU5K9xpEu6IN6tlWWVYaPwUVFrCwrsfQ1nTy4GwCgY7VxTI8XH0OH1vbjgVTCMIpOvRRm2Rj1DBSdYNBrDGRolORcN4f2aW8+zZHHLzkS5x/Rx/L4c5cfhZ8en+lo9Dpfi5feAkHdI/3pCfvj7WtOwGE6+/bjlxyJP513CDppH5C+7n1/pNp+++SlI/HCz5K+Bv17LSlx92F4VYn1jei1YwcDAHp1aIVVd52Ol644BheN2M9R65v1u1PwwPcOszzerlV56pmA5DuUPcLmFoG2VeV45cpj8bMT9jdIT33DZP5wrVZSO3lwt1T5tCovRV+tgSYQLhs1AK9ceawhvVuhbdU4/1yncemp6ZLWGIJGJdn13PWv4Yh+nZRpM53PzlSUlViuVTHp4iOw5LaxGd+ZW3Pbgd3b4HiHEGG1KSn8xlt2IsMOuLCiqASDvkdUQsaZL4f1bo9Vd52O4w/0Pg/TsQO74Pc286NYVRSvGoPXCqeqtG2ryjImz+vergrn16bj8PVn7dexNayQFg+DYCBCZWnyQ5RjIn572qCMc70+i/z4q8pLDCYjIsLwvh1BRIE/yFIi1Nak/SdElHqHch7+w/brgInfGWw8z6SJ6rEyJT324yNTz1RaQqmR6jKp+Qy73nh/h8Z99d3j8D+j1BFQ3dpWooPmCHbyMRClNbfLFR0dOx/Ns7PXGrZVRSXzrlre1IoSImtTXWkJWlWUotKkIbid/+yMQ3o5dnJU78x8yi9H+59qRZLyB2apxS4qwaCvuGY7ICHZEOj3dqquwD9+ONyQ7p5zjdEWQLJy2vWWrCqXucI64amXTZlRSbeNH+rKwUiGHrB1Otkb1qvTRIT2rcvx0PeHY/LVo7D67nG48qRkpIx+JKvXJtxNT8nq/ZxzeG/M+M2JzuebXg0h/WxNZru/blOfN3MeEsL6HZaVZJ4n/5rPSWkMimu9efWoVJSPVe9ZtXvOjaNBRHjp58fg92c61427zxmGG04/CPece4hS2HuJmDGnLaF0vZPzD7mdqcCpStf264gBXavRSwvldZoxWYZOC+G8RrNyHINpl1vLwK9PPdAyYqqKTUnRYQhHJKN6m/4g0y/+iH4dMfZgY3y9qrKWkH2jTan/jFhVmBJS98i8Op7MeXIrWPTJ7M6QgqGi1KiJAcmPq0sb43iIBy88HEfWdDSkc+KfFx+BRyfUusq71TW7t68ymEysMDdW+sguc7iqXtvUn2fOphDC0ccgdFpFWmMwnpQOV828WFV5Kdq1SgoGKwGqcojK8hnQtQ0mHFPj2DsmIpSUEC44cj+Ul5bgmlMPNBx3U0YyiMCc9gSdpv7Do/rht6cNwmWj+mec/8G1Jxm2hY1zX9K2qhzv/PrElDZoF2EGAPtrdUVAWJoCJWqNwZgfK5/GL0yhxXbZkhpDdsRCkQmGClNUkr4cZN3SV7JDemf6G04c1BXHDeyCc4annZhko87K46oPdlCPZFje/RccigO7p807AmozkxeNIakBGfe57dHpU9ndUliYkqwoKaF0qCm5MxucNrQHTjmoe6o3f2ifzNh//fVVOH3cEnMZEQijDuiKA7q1yZgSQ39JvcPX/PgJYR2lIs8TyKx/5us4OZVlcsvX76LonZKY3495Cgs3S2pK/5BdVawoK8GVJw00dJzm3Dgan/zuFOzXKdO06TqUt1T6ixxMSam66awxqAJCzNkZ3jfZGTL7635oitSzW39avotsTeNeVILB4HyG8eNOf5jJv4f37YArTsocLFRdWYanLhuJ+y8wOjHt6iYRDCMkfzX6QDw6oTZlljpneB+cfXjaeS1EMn/m3pH37oLxBNeOK70pyeam8n121s3k6vSNpnrEinQvX3GM5XmVZaV46Ypj8K8J1hPJWQkltyZrs2AhSoZWTrvmBAzu0c5wzKBt2pgoEwqNQQ5gS2sMmfXPjNN6BlKTsHoHYQQEONm3PbgGXHdyHrzocAzu0RadqyvQvV3mqG7h4VoTjkk2xMc5OJTl1QScOy4dFAP1zNnp3q4Kq+8eh1O0KKn0uRUoofRsBHb3klqHF/9LEIpKMBhVfjL0cMjUUxvet6Mnm6md7Y9gtEUKCJxyUHd7GypRKkpI4kVj6NWhKqOCWvWoM25tzEYKc5z5kTUdUzbndHqHBky7ukrgONlih/ftaDtWQT6eOQtnHWYdoqrHS8SH0ZSUmQdJIiFASNrn7zpnGH572iC88YtR2nmkO8+Yea8aQ8Z1XO7X45TEfI0Ms5njHayvZXXumYf2wpRfHm9br1Tf6alDMge9HdKnA1bfPS41nfcxFqvu/eTY/jh9WA/85NgaR2HXUQtnNc56auqQWcwkXFFWgpV3jUtFNNpptvLbyNa6YYU73agD5rqUblTSvbiwICIXPT7jdlNLAtWVZfho4sk49u53DHl0om1lGS6o3Q/z1+0w7Hd7vsHHoNuYe9NoQ8UkIvzP8QNSETvJe7gUPopkQSMuZMNeSoRmIdCmsgwLf+9+0jmVj8EKg7apawgyw1WT+y5UTNugGguRqoemxsWp/vTqUIVV3+yxtGe7KRUnoZ4hGAJYvMOa2kEIdZ3760WHO55rVVfbty7Hwz9ITm3iZIZsr/l26hvTo9UzhXpyh95i0VnX6ZPv3c7CJQNV3JpFg1JUGoMeAik/7pTtO0RrHsEYO+2mbL/cuAtAuuLp8+bESYO7KT9y1422vqHT7W9bVY72CtVZ/417FT4GP09A15p8ZnMj6xa3GhVg3cPN0BhsnM+ywdDXNWsfg33e/nbRcPz5e4elxkOYcVP0TkmcBKeXNivIuYbzoC5nu2kvUhFgLsrbycdQVV6KLm0qDct5mr8zGc6sf38zrz85fVwXhGB5nzIpGByzHArFKxgsNYbMtFaDvLzcS9/j81K2VmYdNzip/pb3JPVv6/Q6k4jDx6Zy8qfODdiJNJ+vmvnSy/l2gkrln0r+NmsMwvIqeh9DSiAorgk4P0vH6gqcdXhvy564u6guJ43B6QrJl+I0jYTqXkHaO6/fhUwv1wC5aMR+mDnxZGVaNz30OTeOxo+Orklf33RcCnX5t6KsxGA2le+1xabVT0dUZkcyFK0pKYNUg5V56M6zh+Hfn3wV4NLkuZFKnWvojbv7AmTV8auuk8VvN7jNoypd4AFqJluu59HVngSpzj9lk6qhKWH5XCnBgHS9SwsI4znlLsvSi9ZjxruPIfMEt9PIh2dKEp7rjUzdqrwUX9wyBm0qyyzz48fZaxUNaHWPXprPo7d5oSIA140djAXrv00JFdYYIobIWOjmAUbh+hhMpgAPF9c3EEGmkUjmw6UpSS+MPH7AfvwYXs+1vmbyAnLE9Xk205QAwEemXmJGVJLNuZ7MJhb7VQPcrKK23DakQaKSrNJIZ21GfXKVI2/38oPX76JEJ5Dbt7JfCdBPQ2zOjwwcsBpAOG5YTzx+yZGYoNM6JCcP7oaHf3BEai63bPkYilZjUI18jhInG7H5qKr4vX5M5nu6brQtfAxucO98zhTAQTUGGZ/epU0lJl89ynGlLvNi85kaQzi1wuoy5ig5w19TWrcap/4ZHvr+cPRoX6ldz78pSTZGjlFJHtos87sOEobptZS8pO/QyseMvaYbyGK2EkBEhBMHdbM4lvzLGkOWKCF1oxTFfOdExt6CJx+DwZTk7b6ZPTzXkkGdATenukw+Upvl9ZSD0h9EUI2hXlvHoHVlGVpXlPnQdjxoDIaETkJffTwVgir0piT1Jd0uJKWP7Bp3SM/U2hlB3m1KMDjUpyjGMbjB86U8nDDWwl/yjmlJXMPlM26X3JMqQ0/vKfnXjYM6TIpWMJh7g+YP0lwAL/78aNx1zjB/94IxXFVVtqNNcdeqCuC2ByvPVc1H4way+O0Gtx/8kJ7tsOqu0w09paCNRf2+pGCo9rsurun29uGq1h/ovy8biQd14ZKWGoOcNA7pyKX0OzCZI1zOfGptSrLv7dudKx/Vqf54ieQLUtb/vmxkarJLq3BV+3u7T0tEGK11XuSyqdeOHYQBXdvYnqPCn5/RKFSyJBeKWDBAHSopi9T8/o/o530JwdS1yXmA0v5d2yjn4ffjfJZ4ics3ptObODzd0jGPqWuTu8bKC3IE6UkWarkT5vvbCWJj3TFyzMAuKT+H6roS1eR7Vj4Gt6urWc+VZL+tuqdEmi+cnPOeNIYALc8xA7vgR7rpJDxHJXns7siyGdqrPd6+5gT87Hj1FOZOuNH6zMJVPltpypTEPoZIIbMpqUTuj8aU1KtDeji/Vc9K36jKFEbns7f7miuia+ezRZ7c4EcrSe0L+O4P7t0es28YnVqAyCtOPgk9Tt+nUetSP5d+sSInx7dr57OlHduUThUVZnFNqR1lCHJXOVKTEa7qsb1LW+EyJ9EzT05nxms1k+++vJQypqz3cx0vmpXMqtQYeYBbxGRMiWE6HmpUEgiH9OmAkwal1V8VqsClIBqDlygbQzpdQu/mW5fmLsW+MCIY/QqFA7q18bSCmWoCRhj2OWtdctDjj4+psbRLm9O6xRyRlWk6zcyUk/PZcYCbh/yZr+V1QKlxtLnx2M8sFiWSpEz9pg/Ran0U+V4Cr4ntwhyYERSj3TsdlRQoC64pWo0BsJp2WzsW4kASWRFrazrh3WVbLK+salQNPU+XDae8fmZUkkuNgdS/3eA4iZ7tueFra27xWtpOTkA377CkhLDiztNRQsDVz87TrqudY0rr1scAAEtuG+so5FRyxiqf8kkzy9Z/eQUt6rQvMLlk66XH9cejH67ScuVgzrQ4/ukNo5Uzr8oOltsykKlqOrfG69q8WID3ZVf11ypj53M86EzfkV3biSaH6YD92kZT2y5L2xiuGuyelijqd4xyAbePT69JfeO4g/w7sDXcmJKAZM85OW17cjuRMtsY07mNSgKSU0E7mZ68DDBMJCyikoL4GAIWtlEjI9x0xhDLfGWerN7dvlU5OrfJ1DilPPCqtZWVlqCNbqLMUh8O5Ewfg6cs+CaQxkBEnQA8B6AGwGoAFwghtivSTQBwo7Z5hxDiCSJqDeA/APYH0ALgv0KIiUHy4xX9OzYv1OM0Zbs33NlT9UPiUz4GvYPSoxjP8DG4bOTJ2Kp5wqlBsvtos6UxfPK7UzLK4GjdTJuXjRqAyyyWwpQY59nKxKvWZR5YmTmJXrh9OG8+BvU5GeZXD3pX5jgG16dqefF3XvJcf50dL1obkPl+/KypLetBtp3PQWvbRADThRAHAJiubRvQhMctAEYCGAHgFiLqqB2+VwgxGMDhAI4lou8EzI83FC9ZFp55xa4gWEWamFGtLKU/xWuFHtHfOK2wn3aXkFxI6OnLRrpL71phsH/WKOnergo92mfO7e8F1ZTtRsjhuCk1GT/8DI3BR6Nifz93+/R5CnMAYGBTkk1tcbq211uXpJzP7ppLyyg0H040cxhzvgxwGw/gCe33EwDOUqQ5DcA0IcQ2TZuYBmCsEKJeCPEuAAghGgF8BsB+DoMIkS9eOnmclv/zgtueVbNuLVqhaCCc6tXpw3rIGwAAhvRqhytPSjvi3E+JYUx3zvA+OHag/eIm6Tz6/+Lj9DF4xTEqSa8xuLieUw/Yy1oRbvByNfkpqNbE9n3/sDQGxbfkHDLt7V7y3bs156kCR5LnJ1/gSYO9h1Onbp0n4ardhRAbtN8bAWSujgH0BrBWt71O25eCiDoA+C6AvwTMjydUkSWy8O1mOvSK/AgmHF2DlVv24IoT1OF0zQotxWxLldR0bo3VW+sNaVV1xs/0Fvp0Xt+C40cZ4NxcwlO4qovHkmmypTF4GRVuNSVGBgE+Gc/BHjZZcXQxpISKO6RA9GrOM2s1pSWED649yVPkXFwag6NgIKK3AajGhd+g3xBCCCLynG0iKgPwDIAHhRArbdJdDuByAOjbN9g02CrSGkPyr773bsW5w/tgV0OTYzpZPaory3Dv+YdapmvS3VNlS5R57N2hFZ74yQic8KcZhuOqxspo8nDMamY6nz05J5RCrCR9jWx9AH6xC3UGjI2oO43B+OHrOwEzfnOi4+p2XlFpINYzjKqPBwlXDQtlPQq5g6Efx+AFVTZUa1bbYZ7cM2cGuAkhRlsdI6JNRNRTCLGBiHoC2KxIth7AibrtPgBm6LYnAVguhPizQz4maWlRW1sbyttRORBl4TtFCAHAfRdYN/J63NZTfaicSmPRT9uh6r2pKo0qJNcJfU/Hc3y5oxpvYxvW/paWEBIuBHOcXDSiL2at3GZ53Oh8dn7vJSbnon7upJou1f4zapk/hfPZIp+yDmSsV2EWDNmarwEOPgbX13BHyvnsMvojTLlk9k/mi4/hNQATtN8TALyqSDMVwBgi6qg5ncdo+0BEdwBoD+CXAfPhC1VvWhZ+qKYkl1VQ79dQ3d4qOiR1PPU3M7opmQ93BKnYQQapmXtHucz4w3rjsR/XAlC/L6eoJTNpE7JmSorYFa8qJys/RsKi3oWZxyAruJkJEDGtRF4ubHOeu3sbv4l8GcdwN4BTiWg5gNHaNoiologeAQAhxDYAtwOYrf27TQixjYj6IGmOGgLgMyKaR0SXBcyPb1IjDKXGEKqPwV06vY9BVQH0vUmnxkidD7cag/trmnG9mJDNdfNBMAD2DaPRf+VCY8gwJbnLw9wbLRV6V/fTY+VctfQxhGBKGtyjra9zU1nxYUryK9D8hJsGJeVj0E2fkg0COZ+FEFsBnKLYPwfAZbrtxwA8ZkqzDtmLUFSieslyEItqBGTUNLc4aAza3xIii8bT3pQ0tHc7V/kIpjH4dz7LrFZXlmFvU4tNyhxBqvmKp9Kb9dyNYzCe56YI7j3/UOWALDeoZICVQ9rSx2CRzi2r7x6Hj1dsxUX/muXtRMW9s4GfkctBkc+ZbR9DUY98VjViKedzDBqDvtDV/oL0fAmqax5Zk5x7/5Jj+6fP0Zrb68YORrsqt4uO6H0M3ggiVNpUluG6sYPx3E+P8n+RLGIr5PSCwc21TAPc3JzktEKdm/sZ9jme43yNbBPm1DXW90ji1UwayvvR+ZqSf7PzzotaMPz0hMzRrVJdDFNjcFtBrh07OPXbyceg6qV2qq7A6rvHpQRE8iSZB9fZ9dW4e528zupz/vmJ+2N/m7nucwnbabn1PgZXGoOxRxilj+GiEX3xyITajP0tDr1RJ43BT2/WziRkf178QsmJMHIo68GALm3wP6P6458XHxHCVZ0pasHQVtGDlrHK4Tqf3dFFZxZIKO6fSCsMyt6LOlzVO358DOcOT/ZerWaoTF07979n11DGjzT64nO3rGbyr9XgqDC565xhOKhnplnRqaydlvZ0E+JtJi0X/H1vWQyE8iyMQlEYUj4Gwg3jhkQSoaaiqGdXVSEdcE0hhkv6Uf9Uva8ubSoAAGcc0suiMbLOs9866vaDvfa0QfjhUX3RrZ27qSayGdoYNap369XHkDHALYR8ecWpM+RUj92EeJvx2/MPo9H1WgXjqLNx9aNYMJiIYkEMP5X4D2dnLiPauU0lFtw6BtUVZdi6pzHjuCrHQSuz29NLSgh9OroZvFM4KoNduXqtP7IH31/rEcZhKvFuSjJuhxnJ5xZfGnHErzZMM2BcJjMWDCYO7NYWPzm2Py4+ul9o1/RTtIft10G5X5q/VKYkuwFuXupX4fTlo8W2AdC9RDdC4uzDe2NQj7YY2qu9du3s49SJcHK+BvHLee2/qLLy5tWj8PGKrVm5l6vzQijEMBav8gMLBhMlJYSbvzvEOaEXIihcVU/i1IMyp6pKB7m4z0Q2NGanW1xybA0O79vRIVW82GsM6d9u2ksiSgkF/bWzKaQdTUkOq675MyV5PsWYB11lPahnO6XvJJ+JeqCjFUXtfA6buTeOxmXH9c/YH0Xhqq5oZ9/3+wGG3TC5zcct3x2KMw/tFfLdw0U+iuqZnEKPna8dgynJoV13Wqc5iF8um2HR2bLOhFKGMWkMLBhCpHObSnSsrsjYH0VFDGOEsc1ZQS9Q9Ojfmq8ItxgaBJUA01czp6m/VbMDO3Fg9+TI5ysc1mm2otBrZ1xRfGxKygKRlK3Li/oJA9S3D1F9eIUkb5xGPgeK7w+Rpy8biV4dWlkeVwoG6FcTNB4zp/YTrtq+VTlW3z3O83lhfFVRVcHBPdtiRP9OuOH0gwJfK65QDdYYskAUkQXykk6rQqWdz7kRDSTXUo5j3pmwsWtYOus0Rz8aQxRv59iBXVJRTyq6aYMUraqKua6ZndWNMUwjk40OhlwJcWjv9g4pk1SWleD5nx6NQy0CSLzAUUkFTBRFK69ZVkKuGh4veYjyW7v1zKHo17kaJw3yvopVrqL6dg/p0wHnHN4bL32+3qfGkP0G4ZbvDsWI/p1x+rAe2LCjAafc914yHxaT6IWhMfglxNkmHBl3SE+MHDDaMADV9rohlh1rDAVMFN+4lAUVIS8SD5hMSSF/6x1aV+BXpx7oaQWxfOUYbTnUPHExoLqyDOcd0QetK8rQV7GgTEaRmZ4rjokng1RPL+e6FQpAuGUXl6JfdIKhprO3FZTCIIoIEykQzjzMPnJHtXa0FwpphHLYOL0aOYpeNb2JE3Fb/uTtq3RTnGSs0wyzKSmLGkPW7hQvcYWrFp0p6Y1fjEJ9Y5andI6gbFtVlOLzm05F26oyjB7SHf0slgz0MFmn7hyh+81YId+Tk03eaUSxCrsG4bLj+mPKoo2er+mFstISXDt2EEYf1B1jHnhfmUb/WFXlJfj9mUMjzZNjJnKEKFZwyzZFJxiqK8tQXRndY6t62FEVrgyNtbPX+3E+R2lKKkSsGvGoNIYbzxiCG88IeRCmgitOHPORe7UAAAzHSURBVGh7XF83Zl1/Cjq0zgzVjorUNOVZu6N74urlh0nRmZLiIM5q0rF1cgqNtlVF1weInSAaA8MA8a1myK1FFohzqcorTx6Irm0rcdZhvV2fw81YOMgJGf34ZOP2MbhBX0+y3UtOTdedg5W1EExJrDGEzI+P7Y/zTStrxfmRV5aV4uKjazxFAenNYTn43eUMbtcv8GVKygNzhMFsmuXshjPtdu7Xbg5XLRDaVJbhT+cfatiXDx+5Ffnw8cSF1YhgiVwjOFdGPvuljYVPLs6aIYWun4GSuTLYU0Ubk8mXB7gVMDlcD5kICeJjkA1fLsjladccj7Xb9mbs97p8aZgM7dUOV500EBeN7Ov53Kg7O0HexfM/PRrTl2zCHW8sSV4rpDx5hTUGJoNcaIwKgYqy5GcdZMnLXKBn+1YY0b+T4kh8FYWI8JvTBqG3zdxPcRHEQtC/SzUuG5Vei559DAUMawyFiVPPU4Zvbq/PXG3PCVln2kYYWh2UGF0MgchlU5IZNiUVMPnmY9APcOvXOTuLj+czVh9vl+rkNAq7Gpp9XfOW7w7BqAO6BspblBiikvKosY2aQngVRS8YBnSpxspv9kR6Dy8VZfqvT8C+puzPOaNH9gT371qNU4dkrgpXSJw4qCvWbK2P5NrtWgX7vC45NnPRp1wi302OUWW/AOQCC4bJV4/yt5CKB7xUlP27toksH145pE/waYNznccvGeH7XKdaQ0QYf1gvjBnSw/c98oVCaAyZNEUvGKrKSyO/R76p2fneE8wacroRmyR/ufDwrGQlDvQmRzdVfNb1p0SYG+9E9VXm2/euIpDzmYg6EdE0Ilqu/VWu3k5EE7Q0y4loguL4a0S0MEhecpl8qyZ+Jt4rZgqgHfCF2w7EHWcdjJvPGIIe7a3XJI8DNiVZEzQqaSKA6UKIAwBM17YNEFEnALcAGAlgBIBb9AKEiM4BsDtgPnKavG048jXfTFZwOyXGD4/qh58clzv+Eq7WzgQVDOMBPKH9fgLAWYo0pwGYJoTYJoTYDmAagLEAQERtAFwD4I6A+chZpv/6hLxTLXm0M+MGridq8uxzVxJUMHQXQmzQfm8EoAph6Q1grW57nbYPAG4HcB+AaMJCcoB8rCNpU1I+5j57pNZjiDkfcXG8LpQ2nxrDEwd1xfC+HXDNqQdGcv186wiqcHQ+E9HbAFRhFTfoN4QQgohcdyGI6DAA+wshfkVENS7SXw7gcgDo29f7MPi4yOdKksdZzwp+1rooJOR6IPlG26pyvHTFsXFnI6dxFAxCiNFWx4hoExH1FEJsIKKeADYrkq0HcKJuuw+AGQCOBlBLRKu1fHQjohlCiBOhQAgxCcAkAKitrWUdNkJOHNQVrStK8eNjauLOCsMwMRDUlPQaABllNAHAq4o0UwGMIaKOmtN5DICpQoi/CyF6CSFqABwH4EsroZDP5GNfslvbKiy+bSwO7t0+7qzkNGxiT1OkSlPBElQw3A3gVCJaDmC0tg0iqiWiRwBACLENSV/CbO3fbdq+ooA/GIZh8o1AA9yEEFsBZIxaEULMAXCZbvsxAI/ZXGc1gIOD5CVXYQdu4cMlzPW80ODZVRnGJ2xJSsOacWHBgiFi+IMpXGQcP5cxU2iwYGCYwLBk4DdQWLBgYBif9NJWDxvG0VtFO5ajUCn62VWjhr+XwuXg3u3xxi+Ow+Ae7eLOCsOECguGiOGeVGEztBdrCwCbkgoNNiVFRIn2pfAHwzBMvsGCISLKSvnVMsUDK8aFBbdeEVGmqQz8wTDFAJtMgatOGhh3FkKDfQwRkRIMIRiTHv7BcHTO05ksGaZY+M1pg/Cb0wbFnY1QYMEQEeWaKSmMjtTpw3oGvwjDZJnnLj8Ks1YWzbRoBQULhogoK01KhARPwckUKSMHdMbIAZ3jzgbjA/YxRERZSfLVNrewYGAYJr9gwRAR5ZrG0NSSiDknDMMw3mDBEBEyXLU5wRoDwzD5BQuGiHjggsNw8uBu6N+lOu6sMAzDeIKdzxExrE97PPbjI+POBsMwjGdYY2AYhmEMsGBgGIZhDLBgYBiGYQywYGAYhmEMsPOZYRgmx7j+O4OxdOOu2O7PgoFhGCbH+OkJ+8d6fzYlMQzDMAZYMDAMwzAGWDAwDMMwBlgwMAzDMAYCCQYi6kRE04houfa3o0W6CVqa5UQ0Qbe/gogmEdGXRLSUiM4Nkh+GYRgmOEE1hokApgshDgAwXds2QESdANwCYCSAEQBu0QmQGwBsFkIcCGAIgPcC5odhGIYJSFDBMB7AE9rvJwCcpUhzGoBpQohtQojtAKYBGKsd+wmAuwBACJEQQnwTMD8MwzBMQIIKhu5CiA3a740AuivS9AawVre9DkBvIuqgbd9ORJ8R0X+ISHU+AICILieiOUQ0Z8uWLQGzzTAMw1jhKBiI6G0iWqj4N16fTgghAHhZlaYMQB8AM4UQwwF8DOBeq8RCiElCiFohRG3Xrl093IZhGIbxguPIZyHEaKtjRLSJiHoKITYQUU8AmxXJ1gM4UbfdB8AMAFsB1AN4Sdv/HwCXuss2wzAMExVBTUmvAZBRRhMAvKpIMxXAGCLqqDmdxwCYqmkY/0VaaJwCYHHA/DAMwzABCSoY7gZwKhEtBzBa2wYR1RLRIwAghNgG4HYAs7V/t2n7AOA6ALcS0XwAFwP4dcD8MAzDMAGhZMc9v6itrRVz5syJOxsMU/QsWLcD89Zux8VH18SdFcYFRDRXCFHrlI5nV2UYxjfD+rTHsD7t484GEzI8JQbDMAxjgAUDwzAMY4AFA8MwDGOABQPDMAxjgAUDwzAMY4AFA8MwDGOABQPDMAxjgAUDwzAMYyAvRz4T0RYAa3ye3gVAsa37wM9cHPAzFwdBnrmfEMJxeuq8FAxBIKI5boaEFxL8zMUBP3NxkI1nZlMSwzAMY4AFA8MwDGOgGAXDpLgzEAP8zMUBP3NxEPkzF52PgWEYhrGnGDUGhmEYxoaiEQxENJaIlhFRHRFNjDs/YUFE+xHRu0S0mIgWEdHV2v5ORDSNiJZrfztq+4mIHtTew3wiGh7vE/iHiEqJ6HMiel3b7k9En2jP9hwRVWj7K7XtOu14TZz59gsRdSCiF4hoKREtIaKjC72ciehXWr1eSETPEFFVoZUzET1GRJuJaKFun+dyJaIJWvrlRDRBdS+3FIVgIKJSAA8B+A6AIQAuIqIh8eYqNJoB/FoIMQTAUQCu1J5tIoDpQogDAEzXtoHkOzhA+3c5gL9nP8uhcTWAJbrtPwJ4QAgxEMB2AJdq+y8FsF3b/4CWLh/5C4ApQojBAA5F8tkLtpyJqDeAXwCoFUIcDKAUwIUovHJ+HMBY0z5P5UpEnQDcAmAkgBEAbpHCxBdCiIL/B+BoAFN129cDuD7ufEX0rK8COBXAMgA9tX09ASzTfv8TwEW69Kl0+fQPQB/tgzkZwOsACMlBP2XmMgcwFcDR2u8yLR3F/Qwen7c9gFXmfBdyOQPoDWAtgE5aub0O4LRCLGcANQAW+i1XABcB+KduvyGd139FoTEgXcEk67R9BYWmOh8O4BMA3YUQG7RDGwF0134Xyrv4M4BrASS07c4AvhVCNGvb+udKPbN2fIeWPp/oD2ALgP/VzGePEFE1CrichRDrAdwL4CsAG5Ast7ko7HKWeC3XUMu7WARDwUNEbQC8COCXQoid+mMi2YUomPAzIjoDwGYhxNy485JFygAMB/B3IcThAPYgbV4AUJDl3BHAeCSFYi8A1cg0uRQ8cZRrsQiG9QD202330fYVBERUjqRQeFoI8ZK2exMR9dSO9wSwWdtfCO/iWABnEtFqAM8iaU76C4AORFSmpdE/V+qZtePtAWzNZoZDYB2AdUKIT7TtF5AUFIVczqMBrBJCbBFCNAF4CcmyL+Rylngt11DLu1gEw2wAB2jRDBVIOrBeizlPoUBEBOBRAEuEEPfrDr0GQEYmTEDS9yD3/0iLbjgKwA6dypoXCCGuF0L0EULUIFmW7wghfgDgXQDnacnMzyzfxXla+rzqWQshNgJYS0SDtF2nAFiMAi5nJE1IRxFRa62ey2cu2HLW4bVcpwIYQ0QdNU1rjLbPH3E7XbLo3DkdwJcAVgC4Ie78hPhcxyGpZs4HME/7dzqSttXpAJYDeBtAJy09IRmhtQLAAiQjPmJ/jgDPfyKA17XfAwB8CqAOwH8AVGr7q7TtOu34gLjz7fNZDwMwRyvrVwB0LPRyBvB7AEsBLATwJIDKQitnAM8g6UNpQlIzvNRPuQL4ifbsdQAuCZInHvnMMAzDGCgWUxLDMAzjEhYMDMMwjAEWDAzDMIwBFgwMwzCMARYMDMMwjAEWDAzDMIwBFgwMwzCMARYMDMMwjIH/D2q/YNXXGGaZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11dd89eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(rs[75, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = tssimul.simulate_arima_given_innov(ar, [], rs)\n",
    "data_test = tssimul.simulate_arima_given_innov(ar, [], rs_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantile Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess = tf.Session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_features = lambda x: x[:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "features_func = simple_features\n",
    "targets_func = lambda x: x[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "symbols_qtl = {}\n",
    "cell_qtl = {}\n",
    "quantiles = [.05, .25, .5, .75, .95]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Either we use `features_func` and `targets_func` to compute `x` and `y`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Iter 0 Train loss 0.0179 Dev loss 0.0179\n",
      "Iter 50 Train loss 0.0174 Dev loss 0.0174\n",
      "Iter 100 Train loss 0.0174 Dev loss 0.0173\n",
      "Iter 150 Train loss 0.0173 Dev loss 0.0172\n",
      "Iter 200 Train loss 0.0171 Dev loss 0.0172\n",
      "Iter 250 Train loss 0.0171 Dev loss 0.0172\n",
      "Iter 0 Train loss 0.3545 Dev loss 0.3537\n",
      "Iter 50 Train loss 0.2873 Dev loss 0.2860\n",
      "Iter 100 Train loss 0.2180 Dev loss 0.2160\n",
      "Iter 150 Train loss 0.1463 Dev loss 0.1451\n",
      "Iter 200 Train loss 0.0864 Dev loss 0.0855\n",
      "Iter 250 Train loss 0.0532 Dev loss 0.0528\n",
      "Iter 0 Train loss 0.4299 Dev loss 0.4274\n",
      "Iter 50 Train loss 0.2282 Dev loss 0.2242\n",
      "Iter 100 Train loss 0.0729 Dev loss 0.0710\n",
      "Iter 150 Train loss 0.0450 Dev loss 0.0441\n",
      "Iter 200 Train loss 0.0450 Dev loss 0.0440\n",
      "Iter 250 Train loss 0.0441 Dev loss 0.0440\n",
      "Iter 0 Train loss 0.5563 Dev loss 0.5548\n",
      "Iter 50 Train loss 0.2646 Dev loss 0.2594\n",
      "Iter 100 Train loss 0.0581 Dev loss 0.0574\n",
      "Iter 150 Train loss 0.0404 Dev loss 0.0403\n",
      "Iter 200 Train loss 0.0399 Dev loss 0.0400\n",
      "Iter 250 Train loss 0.0398 Dev loss 0.0398\n",
      "Iter 0 Train loss 0.0422 Dev loss 0.0422\n",
      "Iter 50 Train loss 0.0397 Dev loss 0.0397\n",
      "Iter 100 Train loss 0.0374 Dev loss 0.0373\n",
      "Iter 150 Train loss 0.0351 Dev loss 0.0350\n",
      "Iter 200 Train loss 0.0328 Dev loss 0.0328\n",
      "Iter 250 Train loss 0.0306 Dev loss 0.0306\n"
     ]
    }
   ],
   "source": [
    "for i, qtl in enumerate(quantiles):\n",
    "    symbols_qtl[qtl], cell_qtl[qtl] = (\n",
    "        tspred_qtl.train_lstm(sess, data[:, 0, :] * 10,\n",
    "                              features_func=features_func,\n",
    "                              targets_func=targets_func,\n",
    "                              quantiles=[qtl], burn_in=burn_in,\n",
    "                              batch_size=50,\n",
    "                              lr0=3e-3, lr_decay=(50, .99), n_iter=300,\n",
    "                              num_units=10, variable_scope=f'rs_run{i}',\n",
    "                              print_every=50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we hand make them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 999, 2), (5000, 999))"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.stack([data[:, 0, :-1] * 10,\n",
    "              rs_regime[:, :-1]], axis=2)\n",
    "y = data[:, 0, 1:] * 10\n",
    "x.shape, y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Iter 0 Train loss 1.5250 Dev loss 1.5176\n",
      "Iter 50 Train loss 0.9165 Dev loss 0.9051\n",
      "Iter 100 Train loss 0.3849 Dev loss 0.3781\n",
      "Iter 150 Train loss 0.0844 Dev loss 0.0793\n",
      "Iter 200 Train loss 0.0222 Dev loss 0.0223\n",
      "Iter 250 Train loss 0.0211 Dev loss 0.0209\n",
      "Iter 0 Train loss 0.7959 Dev loss 0.7284\n",
      "Iter 50 Train loss 0.0792 Dev loss 0.0778\n",
      "Iter 100 Train loss 0.0654 Dev loss 0.0658\n",
      "Iter 150 Train loss 0.0636 Dev loss 0.0640\n",
      "Iter 200 Train loss 0.0624 Dev loss 0.0623\n",
      "Iter 250 Train loss 0.0605 Dev loss 0.0607\n",
      "Iter 0 Train loss 0.5997 Dev loss 0.5966\n",
      "Iter 50 Train loss 0.4097 Dev loss 0.4062\n",
      "Iter 100 Train loss 0.1910 Dev loss 0.1858\n",
      "Iter 150 Train loss 0.0683 Dev loss 0.0682\n",
      "Iter 200 Train loss 0.0632 Dev loss 0.0641\n",
      "Iter 250 Train loss 0.0601 Dev loss 0.0610\n",
      "Iter 0 Train loss 0.4671 Dev loss 0.4590\n",
      "Iter 50 Train loss 0.0623 Dev loss 0.0626\n",
      "Iter 100 Train loss 0.0560 Dev loss 0.0559\n",
      "Iter 150 Train loss 0.0535 Dev loss 0.0523\n",
      "Iter 200 Train loss 0.0503 Dev loss 0.0498\n",
      "Iter 250 Train loss 0.0473 Dev loss 0.0478\n",
      "Iter 0 Train loss 0.0866 Dev loss 0.0863\n",
      "Iter 50 Train loss 0.0835 Dev loss 0.0836\n",
      "Iter 100 Train loss 0.0807 Dev loss 0.0810\n",
      "Iter 150 Train loss 0.0784 Dev loss 0.0784\n",
      "Iter 200 Train loss 0.0757 Dev loss 0.0758\n",
      "Iter 250 Train loss 0.0732 Dev loss 0.0732\n"
     ]
    }
   ],
   "source": [
    "for i, qtl in enumerate(quantiles):\n",
    "    symbols_qtl[qtl], cell_qtl[qtl] = (\n",
    "        tspred_qtl.train_lstm(sess, x, y=y,\n",
    "                              quantiles=[qtl], \n",
    "                              burn_in=burn_in,\n",
    "                              batch_size=50,\n",
    "                              lr0=3e-3, lr_decay=(50, .99), n_iter=300,\n",
    "                              num_units=10, variable_scope=f'rs_run{i + 10}',\n",
    "                              print_every=50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RNN Internal Contents "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Either"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "lstm_input = np.apply_along_axis(features_func, 1, data[:, 0, :])\n",
    "lstm_output = np.apply_along_axis(targets_func, 1, data[:, 0, :])\n",
    "lstm_input_test = np.apply_along_axis(features_func, 1, data_test[:, 0, :])\n",
    "lstm_output_test = np.apply_along_axis(targets_func, 1, data_test[:, 0, :])\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Or"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "lstm_input = x\n",
    "lstm_output = y\n",
    "\n",
    "lstm_input_test = np.stack([data_test[:, 0, :-1] * 10,\n",
    "                            rs_regime_test[:, :-1]], axis=2)\n",
    "lstm_output_test = data_test[:, 0, 1:] * 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 999, 2), (5000, 999, 1))"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "if lstm_input.ndim == 2:\n",
    "    lstm_input = lstm_input[:, :, None]\n",
    "    lstm_input_test = lstm_input_test[:, :, None]\n",
    "if lstm_output.ndim == 2:\n",
    "    lstm_output = lstm_output[:, :, None]\n",
    "    lstm_output_test = lstm_output_test[:, :, None]\n",
    "lstm_input.shape, lstm_output.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def concat_lstm_states(sess, cell, lstm_input):\n",
    "    ''' Return concatenated LSTM states '''\n",
    "    c_states, h_states = rnn_analysis.lstm_states(sess, cell, lstm_input)\n",
    "    return np.concatenate([c_states, h_states], axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "states = {}\n",
    "states_test = {}\n",
    "\n",
    "for qtl in quantiles:\n",
    "    states[qtl] = concat_lstm_states(sess, cell_qtl[qtl], lstm_input)\n",
    "    states_test[qtl] = concat_lstm_states(sess, cell_qtl[qtl],\n",
    "                                          lstm_input_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((5000, 999, 20), (500, 999, 20))"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "states[.5].shape, states_test[.5].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Baseline Explanation of Latent Variable "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Continuous Latent States"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the latent state is continuous, we regress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "n1 = states[.5].shape[0]\n",
    "n2 = states_test[.5].shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [],
   "source": [
    "def reg_latent_cont(regressor,\n",
    "                    states, states_test, quantiles_subset,\n",
    "                    latent, latent_test, burn_in,\n",
    "                    id_traj=random.sample(range(n1), 50),\n",
    "                    fname=None):\n",
    "    ''' Fit and predict states '''\n",
    "    data = np.concatenate([states[q][id_traj]\n",
    "                           for q in quantiles_subset], axis=2)\n",
    "    data_test = np.concatenate([states_test[q]\n",
    "                                for q in quantiles_subset], axis=2)\n",
    "    p = data.shape[-1]\n",
    "    \n",
    "    print(f'data shape: {data.shape}')\n",
    "    print(f'data_test shape: {data_test.shape}')\n",
    "    \n",
    "    # train\n",
    "    if regressor.lower().startswith('l'):\n",
    "        reg = LinearRegression()\n",
    "    else:\n",
    "        reg = RandomForestRegressor(n_jobs=3, n_estimators=100)\n",
    "        \n",
    "    reg.fit(data[:, burn_in:, :].reshape((-1, p)),\n",
    "            latent[id_traj, (burn_in + 1):].ravel())\n",
    "    \n",
    "    # test\n",
    "    correl_test = [stats.spearmanr(reg.predict(s[burn_in:, :]),\n",
    "                                   lat[(burn_in + 1):]).correlation\n",
    "                   for s, lat in zip(data_test, latent_test)]\n",
    "    \n",
    "    plt.hist(correl_test)\n",
    "    plt.grid()\n",
    "    plt.xlabel('Spearman Correlation')\n",
    "    plt.ylabel('Count')\n",
    "    plt.title('Spearman Correlation on Test Data')\n",
    "    # plt.savefig('doc/hist-spearman-full-qtl.png')\n",
    "    if fname is not None:\n",
    "        plt.savefig(fname)\n",
    "    \n",
    "    return reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data shape: (50, 999, 100)\n",
      "data_test shape: (500, 999, 100)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XuYHVWZ7/Hvj4Rwa0iAYA8kkQ7KIEi80QIKM3ZEPSgjoMMojDMmCuZ4ZgYYjZcoKjhnGAFllPEekRN0GEJEkCheQKBBkYCJoiEENGKQhEi4JRBAIPCeP9ZqU9lU997d6b2r0/v3eZ79dN3Xu2pX1Vurqrq2IgIzM7Na21QdgJmZjUxOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCCsqSR1SQpJY4c4/0clnT/ccVlrSZon6UNVx2GDM2oThKTDJf1M0npJD0m6UdIrq46r1SSNk3SGpN9KekzSSkkXSOqqOrZaknokrSoOi4j/iIiTqoqpHkk/kLQhf56W9FSh/ytbsNyz6iVGSX+U9LikRyU9LOknkk6SpAbLeJGkjUONsWR5R0paUTYuImZGxDnDVdaWKKy3DZLWSDpf0g5VxzUSjcoEIWkX4HvA54HdgEnAJ4EnWxzHkM6ah9mlwNHA3wPjgZcCS4AjBrugsvqMkDpWJiLeGBEdEdEBXASc09cfEe9tQQhviIidganAZ4FPAF9qQbkjnqRtJPV3jHtD/s66gVcDH2hSDGOasdxBlL9l+2dEjLoP6UtfN8D4mcCNwBeA9cAdwBGF8eOBrwNrgNXAvwNj8rgXANcCDwIPkA4KEwrzrgQ+DPyalJDG5mEfzMMey8vuBH4APAr8GNi1sIxvAX/Msd0AvLgwbh7wReDKPO/NwAv6qefrgCeAKQOsi72AhcBDwArgPYVxZ5ASzH8DjwAn9TNsG2AO8Lu8XhYAu+VldAEBjM397wKW59jvAv53Hr5TjvVZYEP+7JXL++9CTEcDy4B1QC+wf826/0Bez+uBS4Dt+6n3NsDHgLuBtcA3gPE1Mc8A/pC/59Ma2O7mAf9eMvwtOaZ1wE+AAwrjPk7azh7J6+WvgGOBp4Cn83q4pZ/y/ggcXjPsr/I63LdQ9q/y8v8AfLQw7dpcz771/XLgRXm9PgTcD1wI7NzgfncksKKfcfOBjxWnAz6ay1gNvKMw7Q7A54B7ch0/D2yXx+1B2m/uzzFeAexZmHcR8G+k/eJPwOR66w34L+DbjZSfx38MuA9YBczK63ByoZ7/BVxF2tcPr1OfvwB+mLeNB4FrB9o2CvF9MY9bBXwa2LZm3X48x/i1LTqWbsnMI/UD7JJX9oXAGykcfPP4mcBG4H3AtsDbSQeUvoPa5cBXSQet5wG3sOlA9kLg9cB2eWO9AfhcYdkrgVuBKcAOhWGLSElhEmnH/AVph9yelHBOLyzj3cDOuYzPAbcWxs3LdTuYlHwuAub3sx7OAq6vs65uIJ1xbg+8jLTjvTaPO4N0kDqWdEDdoZ9hp+b6Tc4xfxW4OC+ji80TxFGkJCvgNcDjwCvyuB5gVU18Z5ATBPCXpJ3u9fl7+1DeGcYV1vMtpMSyG2mnem8/9X53nncfoAO4DPhmTcxfy/V7KSnZ719nXc6jJkEAh5J25IOAMaQDym/yd/dSUpLszOtjH2Bq4bs7v055z0kQefha4F25+wjgxfm7egXpoHpkHvciYGPNvC8CXguMIx28FgFnNbjfDSZBPA2clr/Ht5BOGDry+C+TTkImkE7WfkTeP/K6OiZ/L+NJCWJ+oZxFeZ3ul5c9dqD1BjyfdIJ4dmH8QOUfSzoo70c6PizguQniIeCQvM63q7O8zwLn5e1hHPDXefhA28Y5pBONiXn8z8knMHndbiQlyXHkY9CQj6VbMvNI/gD7k3bYVXmFLQQ687iZwL2ACtPfAvxjXuFPFlcscAJwXT/lHAv8stC/Enh3zTQr2fwM6dvAlwv9JwPf6Wf5E/IG2Hd2O4/CgQN4E3BHP/N+jX6SRx4/BXiGwhki8ClgXu4+A7ihZp6yYcvZvAW2J+kAMJaaBFESw3eAU3N3DwMniI8DCwrjtiGdffYU1vM/FMafA3yln3KvAf6p0L9fScyTC+NvAY6vs83N47kJ4v9R0/ogtVoOIR241wDTa9cPW5YgbgVm9zPPV4BP5e7nJIiS6Y8HbmpwnxtMglgPbFMY/wjpBGUsqfU0qTBuOrC8n+UeCqwp9C+i0EoaYL09mj9BOoPfOY8bsHzgf9j8ZO5Anpsg5hbG11veOaQrBvvUxDjQtrGafBKX+48hHwPyun2M3KLY0s+ovAcBEBHLI90Ym0z6EvcinY33WR15jWZ352n2Jp15rJG0TtI60hnx8wAkdUqaL2m1pEdIl1om1hR/T0lI9xW6nyjp78jLH5NvUP4uL39lnqZYxh8L3Y/3zVviQdLBuj97AQ9FxKOFYXeTWjl9yupSO2xv4PLC+lpOSjydtTNKeqOkRfnBgXWkBFe7/gaK9+6+noh4NsdSjLfRdbPZsnL32JqYG13WQPYGPtq3bnKd9yAdMJaRLs2dCayVdJGk56yzIZhEOotF0mGSrpd0v6T1pJOjfte3pL0kfauwfZ8/0PRb4P78/fXpW797kfa/ZYX19R027X8754cs/pDju6okvrJtttYbI927eQPp+LBbHj5g+Xl8cfn19o96yzuTdLJ6naQVkt4P0N+2kR9A+Aueu+1utg9ExNMNrIO6Rm2CKIqIO0hndwcWBk+qedrj+aQv6h5SC2JiREzIn10i4sV5uv8gnTFMi4hdgH8gNQE3K3ILwv170hnB60jN0a48vKEnU2r8GDhY0uR+xt8L7CZp58Kw55POUPqU1aV22D2kHW5C4bN9RBSXg6TtSK2nz5BacxOA77OpbvXW272kA27f8kRqBa3ud44Gl0Wq90Y2T9zD4R7gEzXrZseIuAwgIi6MiFeTLiFsT7rfBUPchiQdDuwO/DQPWkC6FzMlIsaT9oOB1venSWegB+bt+ySGtu0N1RrS9/CCwvoaHxG75/FzSJcyX5nje0NJfA2vu4i4mrR+zm6w/DW5/D5TyhbbaH0iYn1EnBoRewN/C3xM0mF53HO2jXxS+0eeu+3W22eHZFQmiPz43uy+A6OkKaTLRIsKkz0POEXStpL+jnRJ6vsRsYZ0VnKupF3ykxAvkPSaPN/OpBt66yVNIt18Hk47kxLUg8COpIQ0JBHxY+Bq0tn9QZLG5jOw90p6d0TcA/wM+JSk7SW9BDiR1CoajK8AZ0raG0DSHpKOKZluHOma7P3ARklvJO3gfe4Ddpc0vp9yFgBHSTpC0rbAbNK6+tkg4wW4GHifpKmSOkjr+ZKIGLbHPrO5wMmSupV0SDpa0o6SDpD0mpw4n2DTTXpI62LqIB5ZHS/pWNJ3d35E/DbP2wE8GBF/kvRq4O8Ks60Fxkh6fmFY3/b9SB7+/ppy5mvgx3eVt6U/fxqJv08+870AOE/SxLzOpkh6fSG+x4F1kiaSbhhvqXOBYyTt30D5C4CTJO0raad65ddbXt4W9snf1XpSy/vZOtvGxcDpknaX9DzSvZzB7rMNGZUJgnRt8RDgZkmPkRLDbaQDSp+bgX1JT6icCRwXEQ/mce8kHcxuBx4m3WDqu1TzSdLNvvWkJ4kuG+bYv0FqMq7O5S8aePK6jiOdpV9Civk20lNeP87jTyC1Uu4l3Zw/PSeWwTiPdI/nKkmP5pgPqZ0oX8o6hbSTPUxqLS0sjL+DtPHflZvje9XMfyepxfZ50vf2ZuDNEfHUIOOFtNN+k3ST/vekJ15OHsJyBhQRN5Lq/FXSkyq/IdU7SDdazyXVZQ3pYP7xPOt80gnCQ5IGSoBXSdpA2mY+SLqH9N5cduTuz+Tv5UOk6919sT1Muga+JK/vl5Eekz2ctK1cTmrxFU0hPQHYn33YdEB7AnhigBZsf/6VtD0uznH8kPRwCKTW50TSCdRPSdv2FomIe8n3SOqVHxGXk55CvJH0Xf4kzzPQI/QD1Wd/4DrSMesG4DMRcRMDbxufIB0blpHuN91I+h6HnTa/DN8eJM0EToqIw6uOxWxrIWlH0hMzL4mIZ6qOZySQ9HLgJtJDLaPuYDpaWxBmNswi4vGIeHG7JwdJb1V6Q8FEUovtO6MxOYAThJnZYJ1MuvRzJ+nS0CnVhtM8bXmJyczM6nMLwszMSm3VL1qbOHFidHV1DTjNY489xk477dSagEaQdq03tG/d27Xe0L51H2q9lyxZ8kBE7FFvuq06QXR1dbF48eIBp+nt7aWnp6c1AY0g7VpvaN+6t2u9oX3rPtR6S7q7/lS+xGRmZv1wgjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZXaqv+T2sxGjq45V1ZW9rwj2+81G63gFoSZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr1bQEIekCSWsl3VYY9mlJd0j6taTLJU0ojPuIpBWS7pT0v5oVl5mZNaaZLYh5wJE1w64GDoyIlwC/AT4CIOkA4HjgxXmeL0ka08TYzMysjqYliIi4AXioZthVEbEx9y4CJufuY4D5EfFkRPweWAEc3KzYzMysvirvQbwb+EHungTcUxi3Kg8zM7OKVPKTo5JOAzYCFw1h3lnALIDOzk56e3sHnH7Dhg11pxmN2rXe0L51r7res6dtrD9Rk1Rd96o0u94tTxCSZgJ/AxwREZEHrwamFCabnIc9R0TMBeYCdHd3R09Pz4Dl9fb2Um+a0ahd6w3tW/eq6z2z4t+k9nc+/Fp6iUnSkcCHgKMj4vHCqIXA8ZK2kzQV2Be4pZWxmZnZ5prWgpB0MdADTJS0Cjid9NTSdsDVkgAWRcR7I2KZpAXA7aRLT/8cEc80KzYzM6uvaQkiIk4oGfz1AaY/EzizWfGYmdng+D+pzcyslBOEmZmVcoIwM7NSThBmZlaqkn+UM7Pm6arw/xFsdHELwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZqaYlCEkXSFor6bbCsN0kXS3pt/nvrnm4JP2XpBWSfi3pFc2Ky8zMGtPMFsQ84MiaYXOAayJiX+Ca3A/wRmDf/JkFfLmJcZmZWQOaliAi4gbgoZrBxwAX5u4LgWMLw78RySJggqQ9mxWbmZnVp4ho3sKlLuB7EXFg7l8XERNyt4CHI2KCpO8BZ0XET/O4a4APR8TikmXOIrUy6OzsPGj+/PkDxrBhwwY6OjqGr1JbiXatN7Rv3fvqvXT1+qpDabmp48e09Xc+WNOnT18SEd31phs7pKiGQUSEpEFnp4iYC8wF6O7ujp6engGn7+3tpd40o1G71hvat+599Z4558qqQ2m52dM2cu5PH2t5uSvPOqrlZRY1e1tv9VNM9/VdOsp/1+bhq4Ephekm52FmZlaRVieIhcCM3D0DuKIw/J35aaZDgfURsabFsZmZWUHTLjFJuhjoASZKWgWcDpwFLJB0InA38LY8+feBNwErgMeBdzUrLjMza0zTEkREnNDPqCNKpg3gn5sVi5mZDZ7/k9rMzEo5QZiZWSknCDMzK+UEYWZmpZwgzMyslBOEmZmVcoIwM7NSThBmZlbKCcLMzEo5QZiZWSknCDMzK+UEYWZmpZwgzMyslBOEmZmVcoIwM7NSThBmZlbKCcLMzEo5QZiZWSknCDMzK+UEYWZmpZwgzMyslBOEmZmVqiRBSHqfpGWSbpN0saTtJU2VdLOkFZIukTSuitjMzCxpeYKQNAk4BeiOiAOBMcDxwNnAZyPihcDDwImtjs3MzDap6hLTWGAHSWOBHYE1wGuBS/P4C4FjK4rNzMwARUTrC5VOBc4EngCuAk4FFuXWA5KmAD/ILYzaeWcBswA6OzsPmj9//oBlbdiwgY6OjuGtwFagXesN7Vv3vnovXb2+6lBarnMHuO+J1pc7bdL41hdaMNRtffr06UsiorvedGOHFNUWkLQrcAwwFVgHfAs4stH5I2IuMBegu7s7enp6Bpy+t7eXetOMRu1ab2jfuvfVe+acK6sOpeVmT9vIuUtbfjhj5Tt6Wl5mUbO39SouMb0O+H1E3B8RTwOXAYcBE/IlJ4DJwOoKYjMzs6z1KRf+ABwqaUfSJaYjgMXAdcBxwHxgBnBFBbGZDZuuFp/Jz562sS1bD9Y8LW9BRMTNpJvRvwCW5hjmAh8G3i9pBbA78PVWx2ZmZps01IKQdFhE3FhvWKMi4nTg9JrBdwEHD2V5ZmY2/BptQXy+wWFmZjZKDNiCkPQq4NXAHpLeXxi1C+kf3MzMbJSqd4lpHNCRp9u5MPwR0g1lMzMbpQZMEBFxPXC9pHkRcXeLYjIzsxGg0cdct5M0F+gqzhMRr21GUGZmVr1GE8S3gK8A5wPPNC8cMzMbKRpNEBsj4stNjcTMzEaURh9z/a6kf5K0p6Td+j5NjczMzCrVaAtiRv77wcKwAPYZ3nDMzGykaChBRMTUZgdiZmYjS6Ov2nhn2fCI+MbwhmNmZiNFo5eYXlno3p70BtZfAE4QZmajVKOXmE4u9kuaQHott5mZjVJDfd33Y6RfhDMzs1Gq0XsQ3yU9tQTpJX37AwuaFZSZmVWv0XsQnyl0bwTujohVTYjHzMxGiIYuMeWX9t1BeqPrrsBTzQzKzMyq11CCkPQ24Bbg74C3ATdL8uu+zcxGsUYvMZ0GvDIi1gJI2gP4Mem3pc3MbBRq9CmmbfqSQ/bgIOY1M7OtUKMtiB9K+hFwce5/O/D95oRkZmYjQb3fpH4h0BkRH5T0VuDwPOom4KJmB2dmZtWp14L4HPARgIi4DLgMQNK0PO7NTY3OzMwqU+8+QmdELK0dmId1DbVQSRMkXSrpDknLJb0q/8bE1ZJ+m//uOtTlm5nZlquXICYMMG6HLSj3POCHEfEi4KXAcmAOcE1E7Atck/vNzKwi9RLEYknvqR0o6SRgyVAKlDQe+Gvg6wAR8VRErAOOAS7Mk10IHDuU5ZuZ2fBQRPQ/UuoELif953RfQugGxgFviYg/DrpA6WXAXOB2UuthCXAqsDoiJuRpBDzc118z/yxgFkBnZ+dB8+cP/FLZDRs20NHRMdgwt3rtWm8YOXVfunp9S8vr3AHue6KlRY4YVdV92qTxrS+0YKjb+vTp05dERHe96QZMEH+eSJoOHJh7l0XEtYOOaNOyuoFFwGERcbOk84BHgJOLCUHSwxEx4H2I7u7uWLx48YDl9fb20tPTM9Rwt1rtWm8YOXXvmnNlS8ubPW0j5y5t9Mn10aWquq8866iWl1k01G1dUkMJotHfg7gOuG7QUZRbBayKiJtz/6Wk+w33SdozItZI2hNY2+8SzMys6Vr+39D5stQ9kvbLg44gXW5aCMzIw2YAV7Q6NjMz26Sq9ujJwEWSxgF3Ae8iJasFkk4E7ia9FNDMzCpSSYKIiFtJN7trHdHqWMzMrJxfuGdmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKVZYgJI2R9EtJ38v9UyXdLGmFpEskjasqNjMzq7YFcSqwvNB/NvDZiHgh8DBwYiVRmZkZUFGCkDQZOAo4P/cLeC1waZ7kQuDYKmIzM7NEEdH6QqVLgU8BOwMfAGYCi3LrAUlTgB9ExIEl884CZgF0dnYeNH/+/AHL2rBhAx0dHcMa/9agXesNI6fuS1evb2l5nTvAfU+0tMgRo6q6T5s0vvWFFgx1W58+ffqSiOiuN93YIUW1BST9DbA2IpZI6hns/BExF5gL0N3dHT09Ay+it7eXetOMRu1abxg5dZ8558qWljd72kbOXdryXXpEqKruK9/R0/Iyi5q9rVexNR0GHC3pTcD2wC7AecAESWMjYiMwGVhdQWxmZpa1/B5ERHwkIiZHRBdwPHBtRLwDuA44Lk82A7ii1bGZmdkmI+n/ID4MvF/SCmB34OsVx2Nm1tYqvWAZEb1Ab+6+Czi4ynjMzGyTkdSCMDOzEcQJwszMSjlBmJlZKScIMzMr5QRhZmal2vPfLq1tdLX4v5nNRhO3IMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCMDOzUk4QZmZWygnCzMxKOUGYmVkpJwgzMyvlBGFmZqVaniAkTZF0naTbJS2TdGoevpukqyX9Nv/dtdWxmZnZJlW0IDYCsyPiAOBQ4J8lHQDMAa6JiH2Ba3K/mZlVpOU/ORoRa4A1uftRScuBScAxQE+e7EKgF/hwq+MzM2tUlT9pu/Kso5pehiKi6YX0W7jUBdwAHAj8ISIm5OECHu7rr5lnFjALoLOz86D58+cPWMaGDRvo6OgY3sC3Au1ab9i87ktXr684mtbp3AHue6LqKKrRjnWfNmn8kPfz6dOnL4mI7nrTVZYgJHUA1wNnRsRlktYVE4KkhyNiwPsQ3d3dsXjx4gHL6e3tpaenZzhC3qq0a71h87pXeYbXarOnbeTcpS2/KDAitGPdV5511JD3c0kNJYhKnmKStC3wbeCiiLgsD75P0p55/J7A2ipiMzOzpIqnmAR8HVgeEf9ZGLUQmJG7ZwBXtDo2MzPbpIo22WHAPwJLJd2ah30UOAtYIOlE4G7gbRXEZk3Syks9s6dtZGYbXVoya5YqnmL6KaB+Rh/RyljMzKx//k9qMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmalnCDMzKyUE4SZmZVygjAzs1JOEGZmVsoJwszMSo2tOoBako4EzgPGAOdHxFkVhzTsuuZc2fQyZk/byMwWlGNmo9eIShCSxgBfBF4PrAJ+LmlhRNw+3GW14iBtZrY1G2mXmA4GVkTEXRHxFDAfOKbimMzM2pIiouoY/kzSccCREXFS7v9H4JCI+JfCNLOAWbl3P+DOOoudCDzQhHBHunatN7Rv3du13tC+dR9qvfeOiD3qTTSiLjE1IiLmAnMbnV7S4ojobmJII1K71hvat+7tWm9o37o3u94j7RLTamBKoX9yHmZmZi020hLEz4F9JU2VNA44HlhYcUxmZm1pRF1iioiNkv4F+BHpMdcLImLZFi624ctRo0y71hvat+7tWm9o37o3td4j6ia1mZmNHCPtEpOZmY0QThBmZlZqVCQISUdKulPSCklzBpjubyWFpFHzOFy9ukuaKel+Sbfmz0lVxDncGvnOJb1N0u2Slkn6n1bH2CwNfOefLXzfv5G0roo4m6GBuj9f0nWSfinp15LeVEWcw62Beu8t6Zpc515Jk4el4IjYqj+km9m/A/YBxgG/Ag4omW5n4AZgEdBdddytqjswE/hC1bFWUO99gV8Cu+b+51Udd6vqXjP9yaSHPSqPvUXf+1zg/+TuA4CVVcfdonp/C5iRu18LfHM4yh4NLYhGX8/xf4GzgT+1Mrgma9dXkzRS7/cAX4yIhwEiYm2LY2yWwX7nJwAXtySy5muk7gHskrvHA/e2ML5maaTeBwDX5u7rSsYPyWhIEJOAewr9q/KwP5P0CmBKRIy2N/TVrXv2t7npeamkKSXjtzaN1Psvgb+UdKOkRfktwaNBo985kvYGprLpwLG1a6TuZwD/IGkV8H1SC2pr10i9fwW8NXe/BdhZ0u5bWvBoSBADkrQN8J/A7Kpjqch3ga6IeAlwNXBhxfG0yljSZaYe0ln01yRNqDSi1jseuDQinqk6kBY6AZgXEZOBNwHfzMeA0e4DwGsk/RJ4DekNFFv8vY+GFVfv9Rw7AwcCvZJWAocCC0fJjeq6ryaJiAcj4sncez5wUItia6ZGXsmyClgYEU9HxO+B35ASxtZuMK+jOZ7Rc3kJGqv7icACgIi4Cdie9EK7rVkj+/m9EfHWiHg5cFoetsUPJ4yGBDHg6zkiYn1ETIyIrojoIt2kPjoiFlcT7rCq+2oSSXsWeo8GlrcwvmZp5JUs3yG1HpA0kXTJ6a5WBtkkDb2ORtKLgF2Bm1ocXzM1Uvc/AEcASNqflCDub2mUw6+R/XxioaX0EeCC4Sh4q08QEbER6Hs9x3JgQUQsk/Rvko6uNrrmarDup+THPH8FnEJ6qmmr1mC9fwQ8KOl20k27D0bEg9VEPHwGsb0fD8yP/FjLaNBg3WcD78nb+8XAzK19HTRY7x7gTkm/ATqBM4ejbL9qw8zMSm31LQgzM2sOJwgzMyvlBGFmZqWcIMzMrJQThJmZlXKCsMpJOi0/ivvr/AbSQ6qOaThI6pD0VUm/k7Qkv2WzKXWTNE/ScXWmmSlpr0L/+ZIOaEY8NjqMqJ8ctfYj6VXA3wCviIgn8z+1jWtieWPzc+WtcD7we2DfiHhW0lTSS9XqkiTSY+jPFoaN2cLXZswEbiO/wC4iRsWr36153IKwqu0JPND3OpCIeCAi7gWQtFLSOZKWSrpF0gvz8D0kfVvSz/PnsDz8YEk35d8C+Jmk/fLwmZIWSroWuEZSj6TrJV0h6S5JZ0l6Ry5jqaQX5PneLOnmvLwfS+rMw8+QdEFuEdwl6ZTaSuVlHAJ8rO8gHxG/73thpKT3S7otf/41D+tSeuf/N0gH8imSNkg6N//j16skHZRjXyLpRzX/Kd9X9ifyerlN0lwlxwHdwEW5lbZDjr87z3NCrvttks4uLGuDpDMl/UrppYedW/yN29aj6ned+9PeH6ADuJX0rqQvAa8pjFsJnJa73wl8L3f/D3B47n4+sDx37wKMzd2vA76du2eS3s20W+7vAdaRktN2pPfafDKPOxW2bq5pAAAC4ElEQVT4XO7elU3/THoScG7uPgP4WZ53IvAgsG1NvY4GLu+nzgcBS4Gdcv2XAS8HuoBngUML0wbwtty9bS53j9z/dvJvPQDzgONy926F+b8JvDl391L4LZS+fmAv0isq9iBdVbgWOLZQft/855ASXuXbjT+t+fgSk1UqIjZIOgj4K2A6cImkORExL09yceHvZ3P364AD0lUYAHaR1EF6//+FkvYlHdi2LRR1dUQ8VOj/eUSsAZD0O+CqPHxpjgPSS9EuyWfp40iXi/pcGanV86SktaTXG6xqsNqHk5LHY7n8y3L9FwJ3R8SiwrTPAN/O3fuRXjx5da77GGBNyfKnS/oQsCOwGykBfXeAeF4J9EbE/Tmei4C/Jr3P6inge3m6JcDrG6yjjQJOEFa5SNfVe0lv3F0KzCCdEUM60FPTvQ3pLHuzH3+S9AXguoh4i6SuvMw+j9UU+2Sh+9lC/7Ns2i8+D/xnRCyU1ENqOZTN/wzP3ZeWAS8dwn2D2jj/VJhfwLKIeFV/M0vantQS646IeySdQXph3VA9HRF9672snjaK+R6EVUrSfvmMv8/LgLsL/W8v/O17M+lVFH4IRtLLcud4Nr0GeeYwhFdc3ozBzBgRvwMWA5/MN5z77jEcBfwEOFbSjpJ2Iv3Ay08aWOydwB75xj6StpX04ppp+pLBA7lVVXyy6VHS6+9r3UL6LYGJksaQflPh+oYqaqOaE4RVrYN0Weh2Sb8mPeVzRmH8rnn4qcD78rBTgG6lx2JvB96bh58DfErpR1OG40z3DOBbkpYADwxh/pNIl55WSLqN1CpaGxG/yN23ADcD50fEL+stLNLPTR4HnJ1vWt8KvLpmmnXA10g3uX9EelV0n3nAV/puUhfmWQPMIb319lfAkoi4Ygj1tVHGb3O1EUvpB566I2IoB2cz20JuQZiZWSm3IMzMrJRbEGZmVsoJwszMSjlBmJlZKScIMzMr5QRhZmal/j9IZC8BoWR0KAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1278103c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "reg = reg_latent_cont('l',\n",
    "                      states, states_test,\n",
    "                      [.05, .25, .5, .75, .95],\n",
    "                      svvol, svvol_test, burn_in,\n",
    "                      id_traj=random.sample(range(n1), 50))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discrete Latent States"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "n1 = states[.5].shape[0]\n",
    "n2 = states_test[.5].shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clf_latent_disc(classifier,\n",
    "                    states, states_test, quantiles_subset,\n",
    "                    latent, latent_test, burn_in,\n",
    "                    id_traj=random.sample(range(n1), 50),\n",
    "                    fname=None):\n",
    "    ''' Fit and predict states '''\n",
    "    data = np.concatenate([states[q][id_traj]\n",
    "                           for q in quantiles_subset], axis=2)\n",
    "    data_test = np.concatenate([states_test[q]\n",
    "                                for q in quantiles_subset], axis=2)\n",
    "    p = data.shape[-1]\n",
    "    \n",
    "    print(f'data shape: {data.shape}')\n",
    "    print(f'data_test shape: {data_test.shape}')\n",
    "    \n",
    "    # train\n",
    "    if classifier.lower().startswith('l'):\n",
    "        clf = LogisticRegression(C=1e8)\n",
    "    else:\n",
    "        clf = RandomForestClassifier(n_jobs=3, n_estimators=100,\n",
    "                                     criterion='entropy')\n",
    "        \n",
    "    clf.fit(data[:, burn_in:, :].reshape((-1, p)),\n",
    "            latent[id_traj, (burn_in + 1):].ravel())\n",
    "    \n",
    "    # test\n",
    "    cross_entropy = [log_loss(g, clf.predict_proba(s))\n",
    "                     for g, s in zip(latent_test[:, (burn_in + 1):],\n",
    "                                     data_test[:, burn_in:, :])]\n",
    "    \n",
    "    plt.hist(cross_entropy)\n",
    "    plt.grid()\n",
    "    plt.xlabel('Cross Entropy')\n",
    "    plt.ylabel('Count')\n",
    "    plt.title('Cross Entropy on Test Data')\n",
    "    # plt.savefig('doc/hist-ce-full-qtl.png')\n",
    "    if fname is not None:\n",
    "        plt.savefig(fname)\n",
    "    \n",
    "    return clf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data shape: (50, 999, 100)\n",
      "data_test shape: (500, 999, 100)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHiFJREFUeJzt3XuYHGWd9vHvDeEQMkg4vSMkaFAQhQRcMrK6qNcMrCsKC9FFAVGCi5vXVVRMVGD1EnRxN6Io6nrYvICEBRkOgqDgETOCroAJCCEIGiFIAiYcgwPhMPB7/6hnTNM8PdMzPd1dmdyf6+prqp463d0z1b95qqqrFRGYmZlV26TdAczMrJxcIMzMLMsFwszMslwgzMwsywXCzMyyXCDMzCzLBcLMzLJcICxL0rskLZbUL+l+ST+U9Po25jlX0tMpz+DjljqXPVXS+c3O2EqS3lDxOjwuKapem5eMcr1bpnVNHWKe90saqNjWXZLOkvTyEWynV9KnRpPRWscFwl5A0lzgTOA/gE7gJcA3gMNqzD+hRdFOj4iOisc+Y7FSFTaofSEirht8HYC9UvPkitfmT02O0Je2vQ3w5tS2RNIeTd6utVJE+OHHXx8UO3w/8I4h5jkVuBQ4H3gMeB+wBUVRuS89zgS2SPPvAPwAeBR4GLgO2CRNOxFYBfwFuBM4sMY2zwVOqzFtGhDAbOBPwIPAJ9O0g4CngWfS87oltfcBnwN+BawDdgN2Bq5MGZcD/5J5zhelrDcB+6RpHwe+W5Xpq8BXauR9Vdr+o8Ay4NCq5/l14Kq0nRuAlw/zOxt8/hOq2rcDzgP+DNwLnFLxur8S+CWwFngAOC+135jW9Xh6vWZltvd+4GeZ9p8B56fhCcB3gdXpeS4C9kjTPpx+H0+lbVyS2j8N3J2e923Awe3eHzb2R9sD+FGuR3pDHah+s6ma59S0g8+i6IVOBD4LXA/8H2BH4H+Bf0/z/yfwLWCz9HgDIGCP9Ma1c5pvWq03Q+orEP8vZdknvfm8qiLv+VXL9FEUk73Sm9lmwLUUPaUtgVenN84Dqp7z4Wnej6U3s82AndIb6uQ07wRgDTAzk3UziuLzb8DmwAHpDXGPiuf5ELBfWs8FQO8wv7PB519dIH4IfA3YKmW8GZidpl2enoPSa7Z/at8yrWvqENurVSA+ANxT8RocA3SkdX4TuL5i3l7gU1XLH5FybgK8J70uO7R7n9iYHxtUt9paYnvgwYgYGGa+X0fE9yLiuYhYBxwNfDYi1kTEA8BnKHZyKN5YdwJeGhHPRHF4JIBnKXoee0raLCJWRMQfh9jmxyQ9WvFYWDX9MxGxLiJuAW6hKBRDOTcilqXn+mJgf+DEiHgyIn4LnEXxJjdoSURcGhHPAF+ieON7bUTcT1Fc3pHmO4jiNVyS2eZrKd4050fE0xHxc4re1VEV81weETemXBdQFKsRkfRS4I3A3Ih4ImX8KnBkmuUZisLy4vSa/Wqk28i4j6LXQkQMRMR5EdEfEU9S/D3sJ2nLWgtHxEURcX/6m/ofip7lzDHIZaPkAmHVHgJ2qOO8wr1V4zsD91SM35PaAL5A8V/zT9IJzZMAImI5cALFf+dr0onLnantixExueIxu2r6nyuGn6B4I673OewMPBwRf6l6DlNy80fEc8DKiue4EHh3Gn438D81trkzcG9avtZ2Rvo8cl5KUcAeGCyowFcozikBfJSiZ3GzpFslvbvGekZiCsXhOSRNkHRG+n0/BtxB0VvZvtbCko5LWQbz7kZxeNLaxAXCqv2a4vDMrGHmq74N8H0Ub0qDXpLaiIi/RMS8iHgZcCgwV9KBadp3IuL1adkAPt/4Uxg2a679PmA7SVtXtL2E4r/YQbsMDqST2lPTcgDfA/aWNB04hOI//5z7gF2qTopXb2cs3EtxfH/bioL6oojYFyAiVkXEP1P07D4MnJOufGrk9s6zKM4vAbwXeBPQQ3Fe65WpXenn87Yj6RUUh8PmANtFxGSKfyqEtY0LhD1PRKylOFn4dUmzJG0laTNJb5F0+hCLXgh8StKOknZI6zgfQNIhknaTJIqTos8Cz0naQ9IBkrYAnqQ4WfxcjfU3YjUwbagrlSLiXorzJv+ZLvXcGzhu8DkkMyW9PfWuTqAopNen5Z+kOIn9HeDGqH0V0Q0UvYJPpNe1G/hHimPyYyYi7k7ZTpe0taRNJO0+eKmypCMk7ZwO9T2aFns2Ip6i+B29rJ7tSNpU0ssl/TfFeZPT0qStKX6nDwGTKtoHra7aRgfF7/4BYBNJ76foQVgbuUDYC0TEGcBc4FMUO+y9wPEU/yXXchqwGLgVWEpxlc/gm8LuFFe49FP0UL4REYsozj/Mp7jq6M8UJ7hPHmIbn6i61v/BOp/SJennQ5JuGmK+oyiOy99HcRL3lIj4WcX0KyhOpD5CcX7l7el8xKCFwAxqH14iIp6mKAhvoXje3wCOiYg76nwuI3EUMJni8M7DFFdgDR5ieh3FZan9FK/PnIgY7MV8GrgkHeo5tMa6u9OyjwHXUJxw76p4HmdT/O38meLv4ZdVyy8AXpO20RsRN1FcyLAYuB/YNQ1bG6n4B8LMhiLpVGC3iKh5rD4dormD4sTvY63KZtYs7kGYjYF0+GouxSWpLg42LrTqE7Bm45akSRTH1O+huMTVbFzwISYzM8vyISYzM8vaoA8x7bDDDjFt2jQAHn/8cSZNmtTeQEMoez4of0bna1zZMzpfY+rNt2TJkgcjYsdhZ2z3vT4aecycOTMGLVq0KMqs7Pkiyp/R+RpX9ozO15h68wGLw/diMjOz0XKBMDOzLBcIMzPLcoEwM7MsFwgzM8tygTAzsywXCDMzy3KBMDOzLBcIMzPL2qBvtWEbjmknXdXwOubNGODYEa5nxfyDG96u2cbKPQgzM8tygTAzsywXCDMzy3KBMDOzLBcIMzPLcoEwM7MsFwgzM8tygTAzsywXCDMzy3KBMDOzrKYVCEnnSFoj6baKti9IukPSrZIulzS5YtrJkpZLulPSm5uVy8zM6tPMHsS5wEFVbT8FpkfE3sDvgZMBJO0JHAnslZb5hqRNm5jNzMyG0bQCERHXAg9Xtf0kIgbS6PXA1DR8GNAbEU9FxN3AcmC/ZmUzM7PhKSKat3JpGvCDiJiemfZ94KKIOF/SfwHXR8T5adrZwA8j4tLMcnOAOQCdnZ0ze3t7Aejv76ejo6NZT6VhZc8Hzc24dNXahtfRORFWrxvZMjOmbNPwduu1sf+Ox4LzNabefD09PUsiomu4+dpyu29JnwQGgAtGumxELAAWAHR1dUV3dzcAfX19DA6XUdnzQXMzjvQ23TnzZgxwxtKR/cmuOLq74e3Wa2P/HY8F52vMWOdreYGQdCxwCHBgrO++rAJ2qZhtamozM7M2aellrpIOAj4BHBoRT1RMuhI4UtIWknYFdgdubGU2MzN7vqb1ICRdCHQDO0haCZxCcdXSFsBPJUFx3uH9EbFM0sXA7RSHnj4YEc82K5uZmQ2vaQUiIo7KNJ89xPyfAz7XrDxmZjYy/iS1mZlluUCYmVlWWy5zNWuVaWNweW295s0YeN7lvCvmH9yybZs1g3sQZmaW5QJhZmZZLhBmZpblAmFmZlkuEGZmluUCYWZmWS4QZmaW5QJhZmZZLhBmZpblAmFmZlkuEGZmluUCYWZmWS4QZmaW5QJhZmZZLhBmZpblAmFmZlkuEGZmluUCYWZmWS4QZmaW5QJhZmZZTSsQks6RtEbSbRVt20n6qaQ/pJ/bpnZJ+qqk5ZJulbRvs3KZmVl9mtmDOBc4qKrtJOCaiNgduCaNA7wF2D095gDfbGIuMzOrQ9MKRERcCzxc1XwYsDANLwRmVbSfF4XrgcmSdmpWNjMzG54ionkrl6YBP4iI6Wn80YiYnIYFPBIRkyX9AJgfEb9M064BToyIxZl1zqHoZdDZ2Tmzt7cXgP7+fjo6Opr2XBpV9nzQ3IxLV61teB2dE2H1ujEI0yTV+WZM2aZ9YWoo+9+h8zWm3nw9PT1LIqJruPkmjEmqUYiIkDTi6hQRC4AFAF1dXdHd3Q1AX18fg8NlVPZ80NyMx550VcPrmDdjgDOWtu1PdljV+VYc3d2+MDWU/e/Q+Roz1vlafRXT6sFDR+nnmtS+CtilYr6pqc3MzNqk1QXiSmB2Gp4NXFHRfky6mum1wNqIuL/F2czMrELT+uuSLgS6gR0krQROAeYDF0s6DrgHeGea/WrgrcBy4Angvc3KZWZm9WlagYiIo2pMOjAzbwAfbFYWMzMbOX+S2szMslwgzMwsywXCzMyyXCDMzCzLBcLMzLJcIMzMLMsFwszMslwgzMwsywXCzMyyXCDMzCzLBcLMzLJcIMzMLMsFwszMslwgzMwsywXCzMyyXCDMzCzLBcLMzLJcIMzMLMsFwszMslwgzMwsywXCzMyyXCDMzCzLBcLMzLLaUiAkfVTSMkm3SbpQ0paSdpV0g6Tlki6StHk7spmZWaHlBULSFODDQFdETAc2BY4EPg98OSJ2Ax4Bjmt1NjMzW69dh5gmABMlTQC2Au4HDgAuTdMXArPalM3MzABFROs3Kn0E+BywDvgJ8BHg+tR7QNIuwA9TD6N62TnAHIDOzs6Zvb29APT399PR0dGaJzAKZc8Hzc24dNXahtfRORFWrxuDME1SnW/GlG3aF6aGsv8dOl9j6s3X09OzJCK6hptvwpikGgFJ2wKHAbsCjwKXAAfVu3xELAAWAHR1dUV3dzcAfX19DA6XUdnzQXMzHnvSVQ2vY96MAc5Y2vI/2bpV51txdHf7wtRQ9r9D52vMWOdrx97298DdEfEAgKTLgP2ByZImRMQAMBVY1YZsZmNm2hgUxdFYMf/gtmzXxp92nIP4E/BaSVtJEnAgcDuwCDg8zTMbuKIN2czMLGl5gYiIGyhORt8ELE0ZFgAnAnMlLQe2B85udTYzM1uvLQd0I+IU4JSq5ruA/doQx8zMMvxJajMzy3KBMDOzrLoKhKT962kzM7Pxo94exNfqbDMzs3FiyJPUkl4H/B2wo6S5FZNeRHEPJTMzG6eGu4ppc6Ajzbd1RftjrP/MgpmZjUNDFoiI+AXwC0nnRsQ9LcpkTTLcJ3vnzRgYk1timNn4UO/nILaQtACYVrlMRBzQjFBmZtZ+9RaIS4BvAWcBzzYvjpmZlUW9BWIgIr7Z1CRmZlYq9V7m+n1JH5C0k6TtBh9NTWZmZm1Vbw9idvr58Yq2AF42tnHMzKws6ioQEbFrs4OYmVm51FUgJB2Ta4+I88Y2jpmZlUW9h5heUzG8JcWX/NwEuECYmY1T9R5i+lDluKTJQG9TEpmZWSmM9nbfjwM+L2FmNo7Vew7i+xRXLUFxk75XARc3K5SZmbVfvecgvlgxPADcExErm5DHzMxKoq5DTOmmfXdQ3NF1W+DpZoYyM7P2q/cb5d4J3Ai8A3gncIMk3+7bzGwcq/cQ0yeB10TEGgBJOwI/Ay5tVjAzM2uveq9i2mSwOCQPjWBZMzPbANXbg/iRpB8DF6bxI4CrR7vR9DmKs4DpFFdH/TNwJ3ARxXdOrADeGRGPjHYbZmbWmCF7AZJ2k7R/RHwc+G9g7/T4NbCgge1+BfhRRLwS2Af4HXAScE1E7A5ck8bNzKxNhjtMdCbF908TEZdFxNyImAtcnqaNmKRtgDcCZ6f1Ph0RjwKHAQvTbAuBWaNZv5mZjQ1FRO2J0m8i4jU1pi2NiBkj3qD0aorex+0UvYclwEeAVRExOc0j4JHB8arl5wBzADo7O2f29hZ3/Ojv76ejo2OkcVqmDPmWrlo75PTOibB6XYvCjILz1WfGlG1qTivD3+FQnK8x9ebr6elZEhFdw8033DmIF7xBV5g4bIra29wX+FBE3CDpK1QdToqIkJStXBGxgHR4q6urK7q7uwHo6+tjcLiMypDv2JOuGnL6vBkDnLG03tNSred89VlxdHfNaWX4OxyK8zVmrPMNd4hpsaR/qW6U9D6K//xHYyWwMiJuSOOXUhSM1ZJ2SuvfCVhTY3kzM2uB4f7dOQG4XNLRrC8IXcDmwNtGs8GI+LOkeyXtERF3Utw6/Pb0mA3MTz+vGM36zcxsbAxZICJiNfB3knooLkkFuCoift7gdj8EXCBpc+Au4L0UvZmLJR0H3EPxiW0zM2uTer8PYhGwaKw2GhG/peiJVDtwrLZhZmaN8aehzcwsywXCzMyyXCDMzCzLBcLMzLJcIMzMLMsFwszMslwgzMwsywXCzMyyXCDMzCzLBcLMzLJcIMzMLMsFwszMslwgzMwsywXCzMyyXCDMzCzLBcLMzLJcIMzMLMsFwszMslwgzMwsq67vpDazDce0k66qOW3ejAGOHWJ6I1bMP7gp67X2cQ/CzMyyXCDMzCzLBcLMzLLaViAkbSrpZkk/SOO7SrpB0nJJF0navF3ZzMysvT2IjwC/qxj/PPDliNgNeAQ4ri2pzMwMaFOBkDQVOBg4K40LOAC4NM2yEJjVjmxmZlZoVw/iTOATwHNpfHvg0YgYSOMrgSntCGZmZgVFRGs3KB0CvDUiPiCpG/gYcCxwfTq8hKRdgB9GxPTM8nOAOQCdnZ0ze3t7Aejv76ejo6Mlz2E0ypBv6aq1Q07vnAir17UozCg4X+OamXHGlG0aXkcZ9pOhjJd8PT09SyKia7j52vFBuf2BQyW9FdgSeBHwFWCypAmpFzEVWJVbOCIWAAsAurq6oru7G4C+vj4Gh8uoDPmG+4DUvBkDnLG0vJ+ddL7GNTPjiqO7G15HGfaToWxs+Vp+iCkiTo6IqRExDTgS+HlEHA0sAg5Ps80Grmh1NjMzW69Mn4M4EZgraTnFOYmz25zHzGyj1tb+cET0AX1p+C5gv3bmMTOz9crUgzAzsxJxgTAzsywXCDMzy3KBMDOzLBcIMzPLcoEwM7MsFwgzM8tygTAzsywXCDMzy3KBMDOzLBcIMzPLKve9icepacPcdtvMrAzcgzAzsywXCDMzy3KBMDOzLBcIMzPLcoEwM7MsFwgzM8tygTAzsywXCDMzy3KBMDOzLBcIMzPLcoEwM7MsFwgzM8tqeYGQtIukRZJul7RM0kdS+3aSfirpD+nntq3OZmZm67WjBzEAzIuIPYHXAh+UtCdwEnBNROwOXJPGzcysTVpeICLi/oi4KQ3/BfgdMAU4DFiYZlsIzGp1NjMzW08R0b6NS9OAa4HpwJ8iYnJqF/DI4HjVMnOAOQCdnZ0ze3t7Aejv76ejo6M1wUehMt/SVWvbnCavcyKsXtfuFLU5X+OamXHGlG0aXseGtB+XUb35enp6lkRE13Dzta1ASOoAfgF8LiIuk/RoZUGQ9EhEDHkeoqurKxYvXgxAX18f3d3dzYzckMp8Zf3CoHkzBjhjaXm/Q8r5GtfMjCvmH9zwOjak/biM6s0nqa4C0ZarmCRtBnwXuCAiLkvNqyXtlKbvBKxpRzYzMyu04yomAWcDv4uIL1VMuhKYnYZnA1e0OpuZma3Xjv7w/sB7gKWSfpva/g2YD1ws6TjgHuCdbchmZmZJywtERPwSUI3JB7Yyi5mZ1eZPUpuZWZYLhJmZZblAmJlZVrkv2jazDcZYfL5n3owBjh3FesbiMxj2Qu5BmJlZlguEmZlluUCYmVmWC4SZmWW5QJiZWZYLhJmZZblAmJlZlguEmZlluUCYmVnWRvtJ6lZ/q9toPyFqZtYu7kGYmVmWC4SZmWVttIeYzGz8aNUh4+pDxeP9JoHuQZiZWZYLhJmZZblAmJlZlguEmZlluUCYmVmWC4SZmWWVrkBIOkjSnZKWSzqp3XnMzDZWpfochKRNga8DbwJWAr+RdGVE3N7eZGZmL9TqW/ZUasVnMMrWg9gPWB4Rd0XE00AvcFibM5mZbZQUEe3O8FeSDgcOioj3pfH3AH8bEcdXzDMHmJNG9wDuTMM7AA+2MO5IlT0flD+j8zWu7BmdrzH15ntpROw43EylOsRUj4hYACyobpe0OCK62hCpLmXPB+XP6HyNK3tG52vMWOcr2yGmVcAuFeNTU5uZmbVY2QrEb4DdJe0qaXPgSODKNmcyM9soleoQU0QMSDoe+DGwKXBORCyrc/EXHHYqmbLng/JndL7GlT2j8zVmTPOV6iS1mZmVR9kOMZmZWUm4QJiZWdYGUSCGu/2GpPdLWirpt5J+KWnP1P4mSUvStCWSDihTvorpL5HUL+ljZcsnaW9Jv5a0LM2zZZkyStpM0sI07XeSTm5Hvor5/klSSOqqaDs5LXenpDeXKV9Z9pFa+Sram7qPNJqxFftJA7/j0e8jEVHqB8XJ6j8CLwM2B24B9qya50UVw4cCP0rDfwPsnIanA6vKlK+i7VLgEuBjZcpHcRHDrcA+aXx7YNOSZXwX0JuGtwJWANNanS/NtzVwLXA90JXa9kzzbwHsmtYzpq9hg/lKsY/UylcxrWn7yBi8hk3fTxrMN+p9ZEPoQQx7+42IeKxidBIQqf3miLgvtS8DJkraoiz5ACTNAu5O+ZqhkXz/ANwaEbek+R6KiGdLljGASZImABOBp4HKeVuSL/l34PPAkxVth1HsnE9FxN3A8rS+UuQryz5SKx+0ZB9pNGMr9pNG8o16H9kQCsQU4N6K8ZWp7XkkfVDSH4HTgQ9n1vNPwE0R8VRZ8knqAE4EPjPGmcYkH/AKICT9WNJNkj5RwoyXAo8D9wN/Ar4YEQ+3Op+kfYFdIqL67m11Pbc25qvUtn2kVr4W7SMNZaQ1+0kj+Ua9j2wIBaIuEfH1iHg5xR/TpyqnSdqLoqr+33Zkg5r5TgW+HBH97co1qEa+CcDrgaPTz7dJOrBNEWtl3A94FtiZ4hDOPEkva2UuSZsAXwLmtXK79aonXzv3kWHynUoJ9pFhMrZ9Pxkm36j3kVJ9UK6Gkd5+oxf45uCIpKnA5cAxEfHHkuX7W+BwSacDk4HnJD0ZEf9VknwrgWsj4kEASVcD+wLXjGG+RjO+i+J8xDPAGkm/ArqAu1qYb2uK4/d9kgBeDFwp6dA6lm1rvohYXIJ9ZKjXrxX7SKMZW7GfNJJv9PtIM074jPHJmQnpiezK+pMze1XNs3vF8D8Ci9Pw5DT/28uYr2qeU2nOSepGXr9tgZsoTmxNAH4GHFyyjCcC307Dk4Dbgb1bna9q/j7WnyDci+efpL6LsT+B2Ui+UuwjtfJVtTdlHxmD17Dp+0mD+Ua9j5S+BxE1br8h6bMUbxJXAsdL+nvgGeARYHZa/HhgN+DTkj6d2v4hItaUJF/TNZIvIh6R9CWKe2QFcHUMfQy75RkpvmDq25KWAaLYEW5tQ75ayy6TdDHFTjkAfDDG+ARmI/kozz7SVg3+jpu+nzT4Go56H/GtNszMLGvcnKQ2M7Ox5QJhZmZZLhBmZpblAmFmZlkuEGZmluUCYeOOpBdL6pX0x3SH0qslvaIF2z1W0gMq7jg7+NhzmGVOkLRVs7OZjYYvc7VxRcXHSP8XWBgR30pt+1DcDfa6ivkmRMTAGG/7WIoPJx0/gmVWpGUezEzbdKw/M2E2Eu5B2HjTAzwzWBwAIuKWiLhOUrek6yRdSfHBNSTNlXRbepyQ2iZJukrSLan9iNQ+X9Ltkm6V9MV6A6Xt9km6VNIdki5Q4cMU98dZJGlRmrdf0hmSbgFeJ+lASTeruJf/OUp3WpW0QtLpqf1GSbtJ2lrS3ZI2S/O8qHLcbKRK/0lqsxGaDiwZYvq+wPSIuFvSTOC9FPf7EXCDpF9Q3HP/vog4GEDSNpK2B94GvDIiQtLkGus/QtLrK8Zfl37+DcVtN+4DfgXsHxFflTQX6KnoQUwCboiIeSq+dOYPwIER8XtJ5wH/CpyZ5l0bETMkHQOcGRGHSOoDDga+BxwJXBbFPXjMRsw9CNvY3BjF9zJAcefNyyPi8SjuFnoZ8AZgKfAmSZ+X9IaIWAuspbjH/tmS3g48UWP9F0XEqyse6yq2uzIingN+C0yrsfyzwHfT8B7A3RHx+zS+EHhjxbwXVvwcLERnURQ90s9v134pzIbmAmHjzTJg5hDTHx9uBekNeV+KQnGapE+n8xX7Udxb/xDgRyPMVfkdC89Su/f+5AjOO0T1cET8CpgmqZvipoC3jTCn2V+5QNh483NgC0lzBhtUfF/wGzLzXgfMkrSVpEkUh5Cuk7Qz8EREnA98AdhXxRfXbBMRVwMfBfYZo7x/obhVc86dFG/2u6Xx9wC/qJh+RMXPX1e0nwd8B/cerEE+B2HjSjo/8DbgTEknUhwWWgGcQNU3cEXETZLOBW5MTWdFxM2S3gx8QdJzFHeP/VeKN/Er0nkBAXNrRKg+B/GBYSIvAH4k6b6I6KnK96Sk9wKXqPi6yN8A36qYZVtJt1L0To6qaL8AOI31h6DMRsWXuZptgIa5PPZw4LCIeE/Lg9m44h6E2Tgi6WvAW4C3tjuLbfjcgzAzsyyfpDYzsywXCDMzy3KBMDOzLBcIMzPLcoEwM7Os/w/YCV39Ft5K6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x148f576a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "clf = clf_latent_disc('l',\n",
    "                      states, states_test,\n",
    "                      [.05, .25, .5, .75, .95],\n",
    "                      rs_regime, rs_regime_test, burn_in,\n",
    "                      id_traj=random.sample(range(n1), 50),\n",
    "                      fname='doc/hist-ce-full-qtl.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If no regime-related information is provided, the high-vol state can only be identified 10% of time. It is possibly due to the fact that fitting high-variance variables take precedence over low-variance variables in regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[622,  58],\n",
       "       [ 58, 211]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "__states_test = np.concatenate([states_test[q]\n",
    "                                for q in [.05, .25, .5, .75, .95]], axis=2)\n",
    "confusion_matrix(rs_regime_test[227, (burn_in + 1):],\n",
    "                 clf.predict(__states_test[227, burn_in:, :]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}