gibbs-parameters_init_v4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter initialization using Gibbs Sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import truncnorm, invgamma, norm\n",
    "import pickle\n",
    "import time\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "from elec_forecast.gibbs_sampling_model import GibbsSamplingModel\n",
    "from elec_forecast.rm_estimation import RmEstimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Required data loading\n",
    "\n",
    "### Electricity data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mypath=\"data/\"\n",
    "df=pd.read_csv(mypath + \"cdc_conso_daytypes.csv\",sep=\",\",encoding='latin-1')\n",
    "df.index=range(0,len(df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=df.sort_values(by=['Date'])\n",
    "df.index=range(0,len(df))\n",
    "#Valeurs intermédiaires à partir de 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#On se concentre d'abord sur une année afin d'initialiser le filtre\n",
    "df=df[(df.Date >'2016-01-01') & (df.Date <'2017-01-01')]\n",
    "df.index=range(0,len(df))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Temperatures data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_temp=pd.read_csv(mypath + 'temp_2016.csv',sep=',')\n",
    "df_temp=df_temp.drop(df_temp.columns[0],axis=1)\n",
    "df_temp=df_temp.sort_values(by='date')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Common parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_days = 15\n",
    "nb_particles = 10**3\n",
    "N_daytype = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs Sampling For 3AM (03:00)\n",
    "\n",
    "### Parameters required"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#temperatures every day at 3am\n",
    "three_am = \"03:00\"\n",
    "temp_day_3=df_temp[df_temp.hour==three_am]\n",
    "temp_day_3.head()\n",
    "T_h_3=list(temp_day_3.t)\n",
    "\n",
    "temperatures_3 = T_h_3[:nb_days] #france mean temp in jan 2016\n",
    "daytypes_3 = list(df[df['Heure']==three_am]['daytype'][:nb_days])\n",
    "consumptions_3 = list(df[df['Heure']==three_am]['Consommation (MW)'][:nb_days])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#k_day_3 = np.random.dirichlet(np.ones(N_daytype),1)\n",
    "#kappa_3 = k_day_3[0].tolist()\n",
    "kappa_3 = np.array([1/N_daytype]*N_daytype)\n",
    "\n",
    "#u_heat_3 = np.random.normal(14,1)  #fixed parameter, invariant to time\n",
    "u_heat_3 = 13\n",
    "\n",
    "#sigma2_3 = 27000000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Execute Gibbs sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gibbs sampling for particle 1/1000\n",
      "Gibbs sampling for particle 100/1000\n",
      "Gibbs sampling for particle 200/1000\n",
      "Gibbs sampling for particle 300/1000\n",
      "Gibbs sampling for particle 400/1000\n",
      "Gibbs sampling for particle 500/1000\n",
      "Gibbs sampling for particle 600/1000\n",
      "Gibbs sampling for particle 700/1000\n",
      "Gibbs sampling for particle 800/1000\n",
      "Gibbs sampling for particle 900/1000\n",
      "--- 3.549649477005005 seconds ---\n",
      "Gibbs sampling for particle 1/1000\n",
      "Gibbs sampling for particle 100/1000\n",
      "Gibbs sampling for particle 200/1000\n",
      "Gibbs sampling for particle 300/1000\n",
      "Gibbs sampling for particle 400/1000\n",
      "Gibbs sampling for particle 500/1000\n",
      "Gibbs sampling for particle 600/1000\n",
      "Gibbs sampling for particle 700/1000\n",
      "Gibbs sampling for particle 800/1000\n",
      "Gibbs sampling for particle 900/1000\n",
      "--- 3.598416328430176 seconds ---\n",
      "w_init_sum = 696.138838634\n",
      "i = 2\n",
      "x_init_median = 65905.015211\n"
     ]
    }
   ],
   "source": [
    "w_init_sum = 1000\n",
    "x_init_median = 0\n",
    "i=0\n",
    "while (w_init_sum > 0.9*nb_particles or x_init_median > 10**5 or x_init_median < 10**4 ):\n",
    "    sigma2_3 = invgamma.rvs(a=10**(-2), scale=10**(2), size = 1)[0]\n",
    "    GSM_3 = GibbsSamplingModel(temperatures_3,daytypes_3,consumptions_3,nb_days,nb_particles,sigma2_3,kappa_3,u_heat_3)\n",
    "\n",
    "    start_time = time.time()\n",
    "    GSM_3.gibbs_init()\n",
    "    print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "\n",
    "    w_init_sum = np.sum(GSM_3.w)\n",
    "    x_init_median = np.median(GSM_3.x)\n",
    "    i+=1\n",
    "\n",
    "print(\"w_init_sum =\", w_init_sum)\n",
    "print(\"i =\", i)\n",
    "print(\"x_init_median =\", x_init_median)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write initial parameters from Gibbs sampling in a file\n",
    "#timestr = time.strftime(\"%Y%m%d-%H%M%S\")\n",
    "\n",
    "output_file_3 = mypath + 'parameters_init_3am.pkl'\n",
    "parameters = {\n",
    "    's_init':GSM_3.s[nb_days-1,],\n",
    "    'g_heat_init':GSM_3.g_heat[nb_days-1,:],\n",
    "    'sigma_s_init':np.sqrt(GSM_3.sigma_s_star_2[0,nb_particles-1]),\n",
    "    'sigma_g_init':np.sqrt(GSM_3.sigma_g_star_2[0,nb_particles-1]),\n",
    "    'x_init':GSM_3.x,\n",
    "    'w_init':GSM_3.w,\n",
    "    'sigma2': sigma2_3,\n",
    "    'u_h':u_heat_3,\n",
    "    'kappa': kappa_3\n",
    "}\n",
    "output_3 = open(output_file_3, 'wb')\n",
    "pickle.dump(parameters, output_3)\n",
    "output_3.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gibbs Sampling For 3PM (15:00)\n",
    "\n",
    "### Parameters required"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#temperatures every day at 3am\n",
    "three_pm = \"15:00\"\n",
    "temp_day_15=df_temp[df_temp.hour==three_pm]\n",
    "temp_day_15.head()\n",
    "T_h_15=list(temp_day_15.t)\n",
    "\n",
    "temperatures_15 = T_h_15[:nb_days] #france mean temp in jan 2016\n",
    "daytypes_15 = list(df[df['Heure']==three_pm]['daytype'][:nb_days])\n",
    "consumptions_15 = list(df[df['Heure']==three_pm]['Consommation (MW)'][:nb_days])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#sigma2_15 = invgamma.rvs(a=10**(-2), scale=10**(2), size = 1)[0]\n",
    "sigma2_15 = 27000000\n",
    "\n",
    "#k_day_15 = np.random.dirichlet(np.ones(N_daytype),1)\n",
    "#kappa_15 = k_day[0].tolist()\n",
    "kappa_15 = np.array([1/N_daytype]*N_daytype)\n",
    "\n",
    "#u_heat_15 = np.random.normal(14,1)  #fixed parameter, invariant to time\n",
    "u_heat_15 = 13"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Execute Gibbs sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "GSM_15 = GibbsSamplingModel(temperatures_15,daytypes_15,consumptions_15,nb_days,nb_particles,sigma2_15,kappa_15,u_heat_15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gibbs sampling for particle 1/1000\n",
      "Gibbs sampling for particle 100/1000\n",
      "Gibbs sampling for particle 200/1000\n",
      "Gibbs sampling for particle 300/1000\n",
      "Gibbs sampling for particle 400/1000\n",
      "Gibbs sampling for particle 500/1000\n",
      "Gibbs sampling for particle 600/1000\n",
      "Gibbs sampling for particle 700/1000\n",
      "Gibbs sampling for particle 800/1000\n",
      "Gibbs sampling for particle 900/1000\n",
      "--- 4.150750160217285 seconds ---\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "GSM_15.gibbs_init()\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  6.27587845e+37,   6.95676015e+04,   6.53183041e+04,\n",
       "         7.21392914e+04,   6.98536625e+04,   7.27918585e+04,\n",
       "         6.79722981e+04,   7.24051805e+04,   6.87381481e+04,\n",
       "         6.15175092e+04,   6.44169932e+04,   6.76497114e+04,\n",
       "         7.02544870e+04,   7.58951525e+04,   5.88115695e+04,\n",
       "         6.55264875e+04,   6.66222436e+04,   6.60591903e+04,\n",
       "         6.74993977e+04,   5.41190487e+04,   6.42969987e+04,\n",
       "         6.52033008e+04,   6.63350290e+04,   6.75999175e+04,\n",
       "         7.20752981e+04,   6.23426257e+04,   6.46051750e+04,\n",
       "         6.27691066e+04,   5.86754204e+04,   5.81531569e+04,\n",
       "         6.46640137e+04,   7.36408485e+04,   5.84822266e+04,\n",
       "         6.17458532e+04,   6.53182598e+04,   6.29600286e+04,\n",
       "         6.99762416e+04,   7.29280087e+04,   6.76724497e+04,\n",
       "         7.49992052e+04,   6.95162393e+04,   6.47369673e+04,\n",
       "         6.41707950e+04,   6.42332185e+04,   6.55969958e+04,\n",
       "         5.17568253e+04,   6.63900082e+04,   6.65490388e+04,\n",
       "         6.13469962e+04,   7.37166463e+04,   7.27566984e+04,\n",
       "         5.79639464e+04,   6.69083458e+04,   7.65727950e+04,\n",
       "         6.77192770e+04,   6.94785704e+04,   6.41415426e+04,\n",
       "         6.78734904e+04,   6.89703957e+04,   7.42560123e+04,\n",
       "         5.98136331e+04,   5.50412165e+04,   7.21544929e+04,\n",
       "         6.33624037e+04,   5.73715453e+04,   7.46959035e+04,\n",
       "         6.62368640e+04,   6.54832888e+04,   6.94604025e+04,\n",
       "         6.22969302e+04,   6.66139206e+04,   6.32286754e+04,\n",
       "         6.90824968e+04,   7.14124945e+04,   7.20750014e+04,\n",
       "         6.67892554e+04,   6.81020146e+04,   5.28022985e+04,\n",
       "         6.41080874e+04,   7.00611520e+04,   6.67161457e+04,\n",
       "         6.39047022e+04,   6.83345732e+04,   7.60822397e+04,\n",
       "         6.85768049e+04,   6.14623979e+04,   6.03066600e+04,\n",
       "         6.58409899e+04,   6.14207944e+04,   5.86243058e+04,\n",
       "         5.25867098e+04,   6.16555557e+04,   7.12314068e+04,\n",
       "         7.61911881e+04,   6.45080600e+04,   6.61983518e+04,\n",
       "         6.70225498e+04,   5.97084554e+04,   6.85061163e+04,\n",
       "         6.41826050e+04,   6.91838246e+04,   6.65387030e+04,\n",
       "         7.82972489e+04,   7.50097494e+04,   6.72932396e+04,\n",
       "         6.03955450e+04,   6.26097176e+04,   7.30815237e+04,\n",
       "         6.36382891e+04,   7.06299890e+04,   5.42441860e+04,\n",
       "         6.96945433e+04,   6.78604254e+04,   6.04278906e+04,\n",
       "         6.81347584e+04,   6.27269523e+04,   6.80297434e+04,\n",
       "         6.48349368e+04,   7.01153548e+04,   6.00114184e+04,\n",
       "         7.61272967e+04,   6.85399084e+04,   7.07954174e+04,\n",
       "         6.22253148e+04,   6.60541296e+04,   6.96508574e+04,\n",
       "         6.32057867e+04,   6.22751443e+04,   6.78763613e+04,\n",
       "         6.55479963e+04,   7.50840666e+04,   6.91056716e+04,\n",
       "         7.56141766e+04,   6.59553734e+04,   7.14255261e+04,\n",
       "         7.52914174e+04,   6.76028231e+04,   6.14692649e+04,\n",
       "         5.29533241e+04,   6.48051817e+04,   6.16812446e+04,\n",
       "         6.76338879e+04,   6.45236357e+04,   6.91629202e+04,\n",
       "         6.09633751e+04,   6.50373636e+04,   6.48886937e+04,\n",
       "         6.68523836e+04,   7.06908171e+04,   6.32440289e+04,\n",
       "         6.44493961e+04,   5.11687359e+04,   6.57062692e+04,\n",
       "         6.42773880e+04,   6.33571569e+04,   7.09639303e+04,\n",
       "         6.15012231e+04,   7.40011777e+04,   5.53767846e+04,\n",
       "         6.09516481e+04,   6.58719879e+04,   6.26389792e+04,\n",
       "         5.84044438e+04,   8.03359331e+04,   6.24248287e+04,\n",
       "         6.21960497e+04,   6.97900594e+04,   7.13409452e+04,\n",
       "         5.92964874e+04,   7.21995084e+04,   7.33701926e+04,\n",
       "         7.30916205e+04,   6.10187460e+04,   6.20134388e+04,\n",
       "         7.06024180e+04,   6.12556627e+04,   5.11117507e+04,\n",
       "         5.54058357e+04,   6.40233965e+04,   7.65502584e+04,\n",
       "         6.12237122e+04,   6.24122695e+04,   6.40442192e+04,\n",
       "         6.03493092e+04,   7.51568595e+04,   6.55778230e+04,\n",
       "         7.06885631e+04,   7.07374250e+04,   6.68482987e+04,\n",
       "         5.19470826e+04,   6.93352079e+04,   6.76214886e+04,\n",
       "         6.81517742e+04,   6.74496285e+04,   6.17372914e+04,\n",
       "         6.32275820e+04,   6.00533047e+04,   6.20404127e+04,\n",
       "         6.34164984e+04,   6.09807215e+04,   6.64526786e+04,\n",
       "         6.25195102e+04,   5.69729093e+04,   6.94215684e+04,\n",
       "         7.19432170e+04,   6.76633640e+04,   6.64059441e+04,\n",
       "         6.77603108e+04,   5.11178693e+04,   6.20470711e+04,\n",
       "         6.76629429e+04,   7.03115059e+04,   6.60994155e+04,\n",
       "         6.63874541e+04,   7.35313708e+04,   5.64892012e+04,\n",
       "         6.37486127e+04,   7.05021578e+04,   6.57083631e+04,\n",
       "         6.43042310e+04,   6.51545718e+04,   6.62939700e+04,\n",
       "         6.55272822e+04,   6.31000808e+04,   7.02505563e+04,\n",
       "         5.92164530e+04,   6.12789660e+04,   6.42326562e+04,\n",
       "         6.91099692e+04,   6.11676256e+04,   6.81607827e+04,\n",
       "         5.74745553e+04,   5.98256455e+04,   7.81649046e+04,\n",
       "         6.20054588e+04,   6.15755873e+04,   7.41236274e+04,\n",
       "         5.76318719e+04,   6.51610110e+04,   7.18077523e+04,\n",
       "         6.58437629e+04,   6.49647173e+04,   6.14417273e+04,\n",
       "         6.48539663e+04,   6.65450258e+04,   6.41259932e+04,\n",
       "         6.85477107e+04,   7.27064350e+04,   6.79478863e+04,\n",
       "         6.41034030e+04,   6.18814183e+04,   6.53504635e+04,\n",
       "         6.75875187e+04,   6.48308637e+04,   6.47733502e+04,\n",
       "         6.12562144e+04,   6.44078177e+04,   5.77790890e+04,\n",
       "         6.70315652e+04,   6.05280916e+04,   6.30712986e+04,\n",
       "         7.00198726e+04,   6.71786309e+04,   6.44403915e+04,\n",
       "         7.05884386e+04,   7.47556113e+04,   6.38782657e+04,\n",
       "         6.41107165e+04,   6.95449097e+04,   7.49337734e+04,\n",
       "         6.67901753e+04,   7.24716420e+04,   7.25401063e+04,\n",
       "         6.31964559e+04,   6.49954826e+04,   6.39064694e+04,\n",
       "         6.60897999e+04,   6.58526984e+04,   6.49344884e+04,\n",
       "         5.39231975e+04,   6.77073023e+04,   6.81349545e+04,\n",
       "         5.99679430e+04,   6.46901606e+04,   7.44973141e+04,\n",
       "         6.18218118e+04,   6.39673140e+04,   6.85558084e+04,\n",
       "         6.74254204e+04,   6.40620242e+04,   6.52632760e+04,\n",
       "         5.87515693e+04,   6.98685514e+04,   7.59245233e+04,\n",
       "         6.36934417e+04,   6.89778541e+04,   6.42292598e+04,\n",
       "         5.82890376e+04,   6.86630059e+04,   6.94021089e+04,\n",
       "         6.31412990e+04,   7.27071032e+04,   6.83686009e+04,\n",
       "         6.84648000e+04,   7.15306310e+04,   6.70607823e+04,\n",
       "         6.29535862e+04,   5.89143188e+04,   6.78532683e+04,\n",
       "         5.62578150e+04,   6.22922484e+04,   6.37409712e+04,\n",
       "         7.22283020e+04,   6.79940579e+04,   6.79974753e+04,\n",
       "         6.11015341e+04,   6.30590711e+04,   5.75137327e+04,\n",
       "         8.63463892e+04,   7.52879310e+04,   6.08921096e+04,\n",
       "         5.96029130e+04,   7.31612711e+04,   5.79359263e+04,\n",
       "         6.35259931e+04,   6.32488236e+04,   7.16042704e+04,\n",
       "         6.39231518e+04,   6.14512621e+04,   6.87843343e+04,\n",
       "         7.21200856e+04,   6.60860655e+04,   6.57036427e+04,\n",
       "         6.18283817e+04,   6.58869628e+04,   6.34476113e+04,\n",
       "         6.34887537e+04,   7.37066100e+04,   7.46162411e+04,\n",
       "         5.97896258e+04,   6.57397527e+04,   6.56968282e+04,\n",
       "         7.29964495e+04,   7.08137719e+04,   6.44738338e+04,\n",
       "         7.11770882e+04,   6.64894160e+04,   7.25510788e+04,\n",
       "         5.84275745e+04,   5.84409931e+04,   5.72881594e+04,\n",
       "         6.79278416e+04,   6.41738718e+04,   6.67431725e+04,\n",
       "         6.24569597e+04,   7.25373601e+04,   6.86034882e+04,\n",
       "         7.06002026e+04,   6.32937743e+04,   7.21002461e+04,\n",
       "         7.08577919e+04,   6.44829752e+04,   7.49554495e+04,\n",
       "         7.01881223e+04,   6.56906889e+04,   6.39012553e+04,\n",
       "         6.66545336e+04,   5.65731894e+04,   6.12434713e+04,\n",
       "         7.41623927e+04,   6.97512453e+04,   6.56056609e+04,\n",
       "         6.13348659e+04,   6.81730520e+04,   6.34738133e+04,\n",
       "         6.64560262e+04,   7.15772560e+04,   6.86972681e+04,\n",
       "         6.90470113e+04,   5.29912162e+04,   7.19737983e+04,\n",
       "         6.75199657e+04,   6.80721202e+04,   6.35518398e+04,\n",
       "         6.94581179e+04,   6.23443566e+04,   5.62620578e+04,\n",
       "         6.93629289e+04,   6.05956229e+04,   6.19755805e+04,\n",
       "         6.97351968e+04,   6.41072943e+04,   6.36684053e+04,\n",
       "         5.83198492e+04,   6.78767938e+04,   6.62599949e+04,\n",
       "         5.83064249e+04,   6.06111365e+04,   7.24729228e+04,\n",
       "         6.60758994e+04,   6.15715739e+04,   6.29171271e+04,\n",
       "         6.35555242e+04,   6.46345755e+04,   7.28551497e+04,\n",
       "         7.66250720e+04,   6.89410470e+04,   6.99474070e+04,\n",
       "         6.80402415e+04,   6.75812583e+04,   6.76044147e+04,\n",
       "         6.33296807e+04,   6.89204297e+04,   6.72273556e+04,\n",
       "         6.48223032e+04,   6.17816140e+04,   6.34289635e+04,\n",
       "         6.43736316e+04,   7.09967944e+04,   6.49038590e+04,\n",
       "         6.04185919e+04,   6.68615420e+04,   6.37444411e+04,\n",
       "         7.46225970e+04,   6.32436133e+04,   6.21725495e+04,\n",
       "         6.08223961e+04,   6.73142828e+04,   6.45984911e+04,\n",
       "         6.95460801e+04,   6.71646821e+04,   6.80626831e+04,\n",
       "         5.62837333e+04,   6.14579351e+04,   6.59180070e+04,\n",
       "         6.37550740e+04,   6.61965119e+04,   7.16314950e+04,\n",
       "         6.58984455e+04,   6.72174341e+04,   5.77417082e+04,\n",
       "         6.52424102e+04,   6.42899007e+04,   6.81391742e+04,\n",
       "         6.34572388e+04,   7.32906756e+04,   7.18972544e+04,\n",
       "         5.75283784e+04,   6.08069338e+04,   7.36030220e+04,\n",
       "         7.32442972e+04,   7.07096229e+04,   6.76050941e+04,\n",
       "         6.21752550e+04,   6.32351026e+04,   6.15838692e+04,\n",
       "         6.25820016e+04,   5.62222443e+04,   6.26247049e+04,\n",
       "         8.03025877e+04,   6.33385456e+04,   6.07340534e+04,\n",
       "         6.96181620e+04,   6.65554694e+04,   6.91588962e+04,\n",
       "         7.00059880e+04,   6.53417657e+04,   6.02841599e+04,\n",
       "         6.43416421e+04,   6.75271783e+04,   6.81982610e+04,\n",
       "         7.30441916e+04,   6.53143328e+04,   7.29374223e+04,\n",
       "         6.36840360e+04,   6.87079631e+04,   6.88030622e+04,\n",
       "         6.39145272e+04,   6.98456361e+04,   6.32462391e+04,\n",
       "         6.98054139e+04,   7.27317575e+04,   6.97632889e+04,\n",
       "         7.79235246e+04,   6.81666482e+04,   6.70244753e+04,\n",
       "         5.89751939e+04,   7.00059451e+04,   6.79632558e+04,\n",
       "         6.87561308e+04,   6.21946254e+04,   6.59631049e+04,\n",
       "         5.89390341e+04,   6.89839460e+04,   6.70195837e+04,\n",
       "         6.21329806e+04,   7.27241015e+04,   7.45468648e+04,\n",
       "         6.39386796e+04,   6.16913908e+04,   6.63523955e+04,\n",
       "         5.71473352e+04,   6.59072812e+04,   6.74371220e+04,\n",
       "         7.82010330e+04,   6.46624441e+04,   6.84736316e+04,\n",
       "         6.51817476e+04,   6.81118838e+04,   6.14835690e+04,\n",
       "         6.44388174e+04,   6.53604859e+04,   6.75182881e+04,\n",
       "         7.01585703e+04,   6.39427444e+04,   6.31001508e+04,\n",
       "         6.91204196e+04,   7.58881451e+04,   6.64041855e+04,\n",
       "         7.08039197e+04,   6.02465349e+04,   6.06421687e+04,\n",
       "         6.23736436e+04,   6.52051063e+04,   6.83858014e+04,\n",
       "         6.10719903e+04,   6.77413797e+04,   6.81657805e+04,\n",
       "         6.90003033e+04,   6.81983580e+04,   6.47824090e+04,\n",
       "         6.61847055e+04,   6.93895001e+04,   6.08419645e+04,\n",
       "         7.03160363e+04,   6.11992331e+04,   7.09687033e+04,\n",
       "         7.65550548e+04,   6.17355606e+04,   5.99452218e+04,\n",
       "         8.09218466e+04,   5.14210413e+04,   6.81508552e+04,\n",
       "         7.42107413e+04,   5.78682011e+04,   6.20624044e+04,\n",
       "         7.37888744e+04,   6.64260276e+04,   6.96065981e+04,\n",
       "         6.24910353e+04,   6.82893220e+04,   6.00737079e+04,\n",
       "         6.59557768e+04,   7.19603738e+04,   6.89918040e+04,\n",
       "         6.84117818e+04,   6.84362412e+04,   6.14987626e+04,\n",
       "         6.72466851e+04,   7.19971459e+04,   6.52558446e+04,\n",
       "         7.18450661e+04,   6.45075046e+04,   6.19861751e+04,\n",
       "         6.30402369e+04,   6.95934380e+04,   7.31247755e+04,\n",
       "         7.01242678e+04,   6.30835845e+04,   6.42171683e+04,\n",
       "         6.69951731e+04,   6.60457424e+04,   6.32516049e+04,\n",
       "         6.66652443e+04,   7.20302225e+04,   6.84877007e+04,\n",
       "         5.99232048e+04,   6.46793682e+04,   7.24128097e+04,\n",
       "         6.14062569e+04,   5.95452275e+04,   6.70836219e+04,\n",
       "         7.03364206e+04,   6.36134345e+04,   6.10322727e+04,\n",
       "         6.61791503e+04,   6.32087792e+04,   7.49523118e+04,\n",
       "         6.52159842e+04,   6.06716328e+04,   6.98592777e+04,\n",
       "         6.22980672e+04,   6.12665550e+04,   5.87152134e+04,\n",
       "         7.45629928e+04,   6.80257096e+04,   6.47921974e+04,\n",
       "         6.95691491e+04,   5.96559313e+04,   6.33089956e+04,\n",
       "         6.49788702e+04,   6.74099110e+04,   6.60218019e+04,\n",
       "         6.93767576e+04,   6.13857850e+04,   7.57073584e+04,\n",
       "         6.24631560e+04,   6.41016417e+04,   6.94013964e+04,\n",
       "         7.23543634e+04,   6.40735486e+04,   6.87331007e+04,\n",
       "         6.50305753e+04,   6.93816018e+04,   6.44002033e+04,\n",
       "         6.76614880e+04,   5.97698711e+04,   6.13809832e+04,\n",
       "         6.88382926e+04,   6.73694226e+04,   5.97126488e+04,\n",
       "         6.84276365e+04,   6.06257278e+04,   6.43848300e+04,\n",
       "         5.98003312e+04,   6.44693494e+04,   5.80894218e+04,\n",
       "         6.41419053e+04,   6.41899877e+04,   6.26876462e+04,\n",
       "         7.67257040e+04,   5.42297001e+04,   5.76462277e+04,\n",
       "         6.80857180e+04,   7.43852675e+04,   6.44888663e+04,\n",
       "         6.25067939e+04,   7.05903149e+04,   6.93968735e+04,\n",
       "         5.65299814e+04,   6.14543004e+04,   6.81510980e+04,\n",
       "         7.03549946e+04,   6.86369005e+04,   6.26596246e+04,\n",
       "         7.32049856e+04,   5.38600606e+04,   6.12797691e+04,\n",
       "         7.31973388e+04,   6.91570417e+04,   7.35936447e+04,\n",
       "         6.50190743e+04,   7.61839280e+04,   5.85720284e+04,\n",
       "         6.54281338e+04,   6.80317711e+04,   7.67386235e+04,\n",
       "         6.24037705e+04,   6.28424751e+04,   7.27819051e+04,\n",
       "         6.24877130e+04,   7.46382874e+04,   7.21229321e+04,\n",
       "         6.12793935e+04,   7.03446649e+04,   6.21752099e+04,\n",
       "         6.51940352e+04,   6.37853235e+04,   7.03072296e+04,\n",
       "         5.72481467e+04,   6.95677089e+04,   5.91043126e+04,\n",
       "         6.76103466e+04,   6.60447616e+04,   6.96968224e+04,\n",
       "         5.59491974e+04,   6.62976927e+04,   6.30146214e+04,\n",
       "         7.20594027e+04,   5.89838571e+04,   5.91074689e+04,\n",
       "         6.19497011e+04,   5.81894919e+04,   6.16218242e+04,\n",
       "         6.51982464e+04,   6.69818917e+04,   7.48894653e+04,\n",
       "         6.61631008e+04,   6.60447784e+04,   5.90922760e+04,\n",
       "         7.29735882e+04,   6.20845651e+04,   6.70432200e+04,\n",
       "         6.21920922e+04,   6.79851815e+04,   6.34791617e+04,\n",
       "         6.59302875e+04,   6.19517795e+04,   6.30497636e+04,\n",
       "         7.40337310e+04,   6.42816670e+04,   6.48859100e+04,\n",
       "         6.27205118e+04,   6.61350259e+04,   7.12300127e+04,\n",
       "         7.20522556e+04,   7.67045795e+04,   6.78837056e+04,\n",
       "         6.61849858e+04,   6.61833255e+04,   6.34881788e+04,\n",
       "         6.60284389e+04,   6.88762046e+04,   7.00536371e+04,\n",
       "         6.21767003e+04,   6.81025346e+04,   7.87230467e+04,\n",
       "         6.19020833e+04,   6.89411675e+04,   6.99424203e+04,\n",
       "         6.19527692e+04,   6.22079333e+04,   6.97931127e+04,\n",
       "         6.23800214e+04,   6.78820209e+04,   6.94094887e+04,\n",
       "         7.31854785e+04,   7.20418142e+04,   6.32178217e+04,\n",
       "         5.49836766e+04,   6.24596617e+04,   6.57113878e+04,\n",
       "         5.94734001e+04,   6.73366495e+04,   6.81007066e+04,\n",
       "         6.64187204e+04,   5.55731591e+04,   6.71798700e+04,\n",
       "         5.58577176e+04,   6.44276158e+04,   7.65644720e+04,\n",
       "         6.73937479e+04,   6.86787639e+04,   7.43232733e+04,\n",
       "         7.96688825e+04,   7.01041517e+04,   6.27133125e+04,\n",
       "         6.88485244e+04,   6.36306699e+04,   5.78928526e+04,\n",
       "         6.37923443e+04,   6.34573932e+04,   6.26111575e+04,\n",
       "         6.36644403e+04,   6.85798394e+04,   7.10876833e+04,\n",
       "         6.53050982e+04,   7.77079808e+04,   6.88792035e+04,\n",
       "         7.07205497e+04,   6.68163446e+04,   5.81659477e+04,\n",
       "         6.22026085e+04,   7.26882162e+04,   6.50631949e+04,\n",
       "         6.79778712e+04,   7.21838530e+04,   6.49019943e+04,\n",
       "         5.20207050e+04,   5.82119714e+04,   7.52488435e+04,\n",
       "         6.97680551e+04,   6.53559307e+04,   6.85113724e+04,\n",
       "         6.79945405e+04,   7.09423150e+04,   7.34438594e+04,\n",
       "         6.53313519e+04,   6.02428211e+04,   7.17347428e+04,\n",
       "         6.78808959e+04,   6.45028739e+04,   6.88665932e+04,\n",
       "         6.63963667e+04,   7.81703500e+04,   6.75052341e+04,\n",
       "         6.32929106e+04,   6.39814375e+04,   5.93678192e+04,\n",
       "         7.26621184e+04,   6.18717126e+04,   7.00110910e+04,\n",
       "         6.43667113e+04,   7.08942672e+04,   7.36535260e+04,\n",
       "         7.30134101e+04,   5.80228279e+04,   6.72801617e+04,\n",
       "         7.08616570e+04,   7.49307358e+04,   6.58543012e+04,\n",
       "         7.08342981e+04,   6.61899911e+04,   6.29807167e+04,\n",
       "         6.92212442e+04,   6.30667345e+04,   8.05450966e+04,\n",
       "         6.34884706e+04,   7.09496615e+04,   6.50318677e+04,\n",
       "         6.24532456e+04,   6.96385061e+04,   6.37153940e+04,\n",
       "         5.90777650e+04,   6.54528807e+04,   7.07830092e+04,\n",
       "         6.08887914e+04,   6.28078785e+04,   7.25242037e+04,\n",
       "         6.62047151e+04,   6.64287486e+04,   5.72900664e+04,\n",
       "         6.48488657e+04,   7.24792155e+04,   6.91273381e+04,\n",
       "         6.84047570e+04,   6.13436316e+04,   5.86600625e+04,\n",
       "         6.84308933e+04,   6.94569536e+04,   6.26852187e+04,\n",
       "         5.92548387e+04,   7.26831584e+04,   6.16265130e+04,\n",
       "         6.57597492e+04,   5.83990670e+04,   7.10883408e+04,\n",
       "         6.64854338e+04,   6.09392132e+04,   6.83778940e+04,\n",
       "         5.91203537e+04,   7.50557252e+04,   6.50860238e+04,\n",
       "         6.22763587e+04,   6.60728244e+04,   6.33025992e+04,\n",
       "         6.42740680e+04,   6.85208112e+04,   7.25336303e+04,\n",
       "         5.96921027e+04,   7.66079937e+04,   7.18533330e+04,\n",
       "         6.71224514e+04,   6.50619426e+04,   6.37627584e+04,\n",
       "         6.60913300e+04,   5.82554486e+04,   6.41312828e+04,\n",
       "         6.25875306e+04,   7.01715068e+04,   6.15233775e+04,\n",
       "         6.41155476e+04,   7.17858219e+04,   6.38953499e+04,\n",
       "         6.19827187e+04,   6.99116383e+04,   6.44703372e+04,\n",
       "         5.68763379e+04,   5.99850986e+04,   6.74488208e+04,\n",
       "         6.71915209e+04,   7.22199113e+04,   7.57004305e+04,\n",
       "         6.36135741e+04,   7.10730086e+04,   7.43776883e+04,\n",
       "         6.94088184e+04,   6.45732393e+04,   6.41245421e+04,\n",
       "         7.00886690e+04,   6.48611948e+04,   6.74276955e+04,\n",
       "         6.63359711e+04,   6.24457258e+04,   6.82291185e+04,\n",
       "         6.32657983e+04,   6.10417465e+04,   6.06643502e+04,\n",
       "         7.04301171e+04,   7.23567700e+04,   6.34282783e+04,\n",
       "         7.74488563e+04,   6.48320293e+04,   6.54180694e+04,\n",
       "         6.93132515e+04,   6.83916087e+04,   6.25670261e+04,\n",
       "         6.44216814e+04,   6.78933014e+04,   6.34519399e+04,\n",
       "         6.77466449e+04,   5.97016380e+04,   6.28655192e+04,\n",
       "         6.93287071e+04,   7.26306246e+04,   6.38178561e+04,\n",
       "         5.50645729e+04,   6.81126096e+04,   6.89377180e+04,\n",
       "         7.18068376e+04,   6.60712399e+04,   5.79348400e+04,\n",
       "         6.51805755e+04,   6.38722181e+04,   5.92339482e+04,\n",
       "         7.35089588e+04,   6.79741164e+04,   6.29852778e+04,\n",
       "         6.15508343e+04,   6.74771875e+04,   6.21878769e+04,\n",
       "         5.30368075e+04,   6.39095739e+04,   6.40727213e+04,\n",
       "         5.80415832e+04,   5.96564720e+04,   6.62095307e+04,\n",
       "         6.70191905e+04,   6.61930678e+04,   6.63641651e+04,\n",
       "         6.19772567e+04,   6.53677474e+04,   6.85527819e+04,\n",
       "         7.38583044e+04,   5.75249807e+04,   6.25593060e+04,\n",
       "         6.50436652e+04,   6.62483612e+04,   6.99773561e+04,\n",
       "         7.49504419e+04,   6.15188769e+04,   6.37149620e+04,\n",
       "         5.10631327e+04,   7.10398899e+04,   6.37827707e+04,\n",
       "         6.23636327e+04,   6.80457326e+04,   6.76107641e+04,\n",
       "         6.39637983e+04,   5.99196487e+04,   7.06559098e+04,\n",
       "         6.33189493e+04,   5.69906953e+04,   6.39008916e+04,\n",
       "         6.21420560e+04,   6.86472602e+04,   5.90980615e+04,\n",
       "         7.57456691e+04,   6.53855896e+04,   7.22312892e+04,\n",
       "         6.38004120e+04,   6.48274317e+04,   6.29099575e+04,\n",
       "         5.44082882e+04,   6.21807265e+04,   7.20309150e+04,\n",
       "         6.75420687e+04,   7.29757770e+04,   7.07939920e+04,\n",
       "         6.92587530e+04,   6.82445932e+04,   6.59159422e+04,\n",
       "         6.49749894e+04,   7.44632911e+04,   6.42737131e+04,\n",
       "         5.61951532e+04,   6.31452920e+04,   6.74417951e+04,\n",
       "         5.38705366e+04,   5.43084458e+04,   6.46255576e+04,\n",
       "         5.77052490e+04,   6.31469556e+04,   7.00293466e+04,\n",
       "         6.94156814e+04,   6.80864212e+04,   6.19660823e+04,\n",
       "         7.29581233e+04,   6.71107700e+04,   6.18885525e+04,\n",
       "         5.80533024e+04,   7.19602684e+04,   7.56234242e+04,\n",
       "         6.32077686e+04,   6.49055351e+04,   6.44415253e+04,\n",
       "         7.07986814e+04,   6.25217281e+04,   6.93594005e+04,\n",
       "         6.43770857e+04])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.median(BFM_15.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Write initial parameters from Gibbs sampling in a file\n",
    "#timestr = time.strftime(\"%Y%m%d-%H%M%S\")\n",
    "\n",
    "output_file_15 = mypath + 'parameters_init_3pm.pkl'\n",
    "parameters = {\n",
    "    's_init':GSM_15.s[nb_days-1,],\n",
    "    'g_heat_init':GSM_15.g_heat[nb_days-1,:],\n",
    "    'sigma_s_init':np.sqrt(GSM_15.sigma_s_star_2[0,nb_particles-1]),\n",
    "    'sigma_g_init':np.sqrt(GSM_15.sigma_g_star_2[0,nb_particles-1]),\n",
    "    'x_init':GSM_15.x,\n",
    "    'w_init':GSM_15.w,\n",
    "    'sigma2': sigma2_15,\n",
    "    'u_h':u_heat_15,\n",
    "    'kappa': kappa_15\n",
    "}\n",
    "output_15 = open(output_file_15, 'wb')\n",
    "pickle.dump(parameters, output_15)\n",
    "output_15.close()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}