{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the last few posts we looked at how government spending and taxation interact with the circular flow of money to increase aggregate incomes. This process was predicated on an injection of money that is gradually depleted as it is taxed out of circulation during repeated spending. In the interim period an amount of aggregate income is generated which is greater than the injected amount. \n",
      "\n",
      "Just for the sake of it, we'll revisit a simpler case, without government money, and see if we can understand the circulation of money in terms of a similar multiplier description and accounting framework.\n",
      "\n",
      "When we first looked at the circular flow of money, we envisaged a fixed amount of money already in existence simply being passed around. We used the simple equations:\n",
      "\n",
      "$$Y = C$$\n",
      "$$C = Y$$\n",
      "\n",
      "The first equation simply describes the identity which pairs spending with incomes - every act of spending generates an equivalent amount of income for someone else. The second equation is a statement about the behaviour of folk in our economy. It states that folk simply spend all of their income: spending is equal to income. We call this type of equation a *consumption function*.\n",
      "\n",
      "This is certainly a highly simplified model of an economy with ridiculously implausible assumptions and ludicrous ommissions. but even such daft models can have plenty of heuristic value. In this case, there are some issues that we can improve on and that might help to nail down some of the concepts that might be useful in more complicated models.\n",
      "\n",
      "Money stock/wealth not reference\n",
      "Only reference to flows suggests infinite speed\n",
      "have to specify or infer velocity, but this always implies that folk are holding money at some stage\n",
      "Where did money start?\n",
      "\n",
      "had to know how many times the total money stock was spent in a particular time period in order to know the total cumulative spending and income. This was described in terms of the *velocity of money*. Then when we initially looked at the fiscal multiplier, we found that  \n",
      "\n",
      "\n",
      "\n",
      "$$H_t = H_0$$\n",
      "\n",
      "$$Y_t = C_t$$\n",
      "\n",
      "$$C_t = \\alpha_Y Y_t + \\alpha_H H_{t-1}$$\n",
      "\n",
      "\n",
      "...\n",
      "\n",
      "\n",
      "$$Y_t = \\alpha_Y Y_t + \\alpha_H H_{t-1}$$\n",
      "\n",
      "$$Y_t (1 - \\alpha_Y) = \\alpha_H H_{t-1}$$\n",
      "\n",
      "$$Y_t = \\frac {\\alpha_H H_{t-1}}{1 - \\alpha_Y}$$\n",
      "\n",
      "$$Y_t = \\frac {\\alpha_H }{1 - \\alpha_Y} H_{t-1}$$\n",
      "\n",
      "$$Y_t = \\frac {\\alpha_H }{1 - \\alpha_Y} H_{0}$$\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "alpha_h/(1-alpha_y)*100"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "50.000000000000014"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "%matplotlib inline  \n",
      "\n",
      "N = 200\n",
      "\n",
      "alpha_y = 0.9\n",
      "alpha_h = 0.5\n",
      "\n",
      "C   = np.zeros(N) # consumption\n",
      "Y   = np.zeros(N) # income\n",
      "H_h = np.zeros(N) # private wealth\n",
      "\n",
      "C[0] = 0\n",
      "Y[0] = 0\n",
      "H_h[0] = 100\n",
      "\n",
      "for t in range(1, N):\n",
      "    \n",
      "    # calculate consumer spending\n",
      "    C[t] = alpha_y*Y[t-1] + alpha_h*H_h[t-1] \n",
      "        \n",
      "    # calculate total income (consumer spending plus constant government spending)\n",
      "    Y[t]   = alpha_h/(1-alpha_y)*H_h[t-1] \n",
      "    \n",
      "    H_h[t] = H_h[t-1] +(1-alpha_y)*Y[t] - alpha_h*H_h[t-1] \n",
      "   "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(figsize=(12, 8))\n",
      "\n",
      "consumption_plot = fig.add_subplot(131, xlim=(0, N), ylim=(0, np.max(C)*1.1))\n",
      "consumption_plot.plot(range(N), C, lw=3)\n",
      "consumption_plot.grid()\n",
      "# label axes\n",
      "plt.xlabel('time')\n",
      "plt.ylabel('consumption')\n",
      "\n",
      "income_plot = fig.add_subplot(132, xlim=(0, N), ylim=(0, np.max(Y)*1.1))\n",
      "income_plot.plot(range(N), Y, lw=3)\n",
      "income_plot.grid()\n",
      "# label axes\n",
      "plt.xlabel('time')\n",
      "plt.ylabel('income')\n",
      "\n",
      "\n",
      "debt_plot = fig.add_subplot(133, xlim=(0, N), ylim=(0,np.max(H_h)*1.5))\n",
      "debt_plot.plot(range(N), H_h, lw=3)\n",
      "debt_plot.grid()\n",
      "# label axes\n",
      "plt.xlabel('time')\n",
      "plt.ylabel('government debt')\n",
      "\n",
      "# space subplots neatly\n",
      "plt.tight_layout()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAI5CAYAAAC8QJvjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADx0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wcmMx\nLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvyHfiKQAAIABJREFUeJzs3Xt01NW9///XJIEUES+gITFR\nhAYjAREBEVtLY2O0FhtUagTxGESbfg+rS1vaKj059fJbPSa6jlbr4adNF2oq1ojWRbwGFZlSqR5E\nrHxLvmC+FE5JiBFKKIRLIJfvH9TpjKBOYPb+zN7zfKzlWjCSfPYQPu/Zr733Z+9Qb29vrwAAAAAA\nxywt6AYAAAAAgC8IWAAAAACQIAQsAAAAAEgQAhYAAAAAJAgBCwAAAAAShIAFAAAAAAlCwAIAAACA\nBCFgAQAAAECCELAAAAAAIEEygm7AsTjppJOUn58fdDOM2bNnjwYOHBh0M4zh/blt8+bN2r59e9DN\nsIqa4zben7uoN/7x+d+rxPtz3bHWHKcD1tChQ7V69eqgm2FMOBxWUVFR0M0whvfntokTJwbdBOuo\nOW7j/bmLeuMfn/+9Srw/1x1rzWGJIAAAAAAkCAELAAAAABKEgAUAAAAACULAAgAAAIAEIWABAAAA\nQIIQsAAAAAAgQQhYAAAAAJAgBCwAAAAASBACFgAAAAAkCAELAAAAABKEgAUAAAAACULAAgAAAIAE\nIWABAAAAQIIQsAAAAAAgQQhYAAAAAJAgBCwAAAAASBACFgAAAAAkCAELAAAAABKEgAUAAAAACULA\nAgAAAIAEIWABAAAAQIIQsAAAAAAgQQhYAAAAAJAgBCwAAAAASBACFgAAAAAkCAELAAAAABKEgAUA\nAAAACULAAgAAAIAEIWABAAAAQIIQsAAAAPpgzpw5ysrK0pgxYw77f/fff79CoZC2b98eea2qqkr5\n+fkqKCjQ0qVLbTYVQAAIWAAAAH0we/ZsNTQ0HPb6li1b9Nprr+mMM86IvNbY2Ki6ujqtW7dODQ0N\nmjt3rrq7u202F4BlBCwAAIA+mDJligYPHnzY6z/84Q913333KRQKRV6rr6/XjBkzlJmZqeHDhys/\nP1+rVq2y2VwAlhGwAAAAjlF9fb1yc3N17rnnxrze0tKi008/PfL7vLw8tbS02G4eAIsygm4AAACA\ny/bu3at77rlHr7322jF9n5qaGtXU1EiS2tvbFQ6HE9C65NTR0cH7c5jv7+9YEbAAAACOwcaNG7Vp\n06bI7FVzc7PGjx+vVatWKTc3V1u2bIn82ebmZuXm5h7x+1RUVKiiokKSVFBQoKKiIuNtD0o4HOb9\nOcz393esWCIIAABwDM455xx9/PHH2rx5szZv3qy8vDytWbNG2dnZKi0tVV1dnTo7O7Vp0yY1NTVp\n0qRJQTcZgEEELAAAgD6YOXOmLrzwQm3YsEF5eXlauHDhZ/7Z0aNHq6ysTIWFhfrmN7+pBQsWKD09\n3WJrAdjGEkEAAIA+ePrppz/3/2/evDnm95WVlaqsrDTYIgDJhBksAAAAAEgQAhYAAAAAJAgBCwAA\nAAAShIAFAAAAAAliNGCdeeaZOuecczRu3DhNnDhRkrRjxw6VlJRo5MiRKikpUXt7e+TPV1VVKT8/\nXwUFBVq6dKnJpgHwDPUGAAAkA+O7CC5fvlynnHJK5PfV1dUqLi7W/PnzVV1drerqat17771qbGxU\nXV2d1q1bp61bt+qSSy7Rhx9+GNdWpns6u7T/YLfJtxGIXQd69beOzqCbYcQJA/oF3QR4yHS96ZW8\nvSclf2tOZr90HZ/JprkAADusf+LU19crHA5LksrLy1VUVKR7771X9fX1mjFjhjIzMzV8+HDl5+dr\n1apVuvDCCz/3+z36+426/7UNOtjda6H1AXjzjaBbYMSpgzL1/XNCQTcDnkt0vdmyu0cTfu7nPRnh\nYc3JSAtp7sX5Gs+4DgDAAqMBKxQK6ZJLLlF6erq+973vqaKiQm1tbcrJyZEkZWdnq62tTZLU0tKi\nyZMnR742Ly9PLS0th33Pmpoa1dTUSJLa29v1/y/zOFx5bNvuTr25qVfD/tH59VFHR0ekcw/zTNQb\nKbbm9FBqnNTV06tHw0164MJer+9Jag4AJAejAeutt95Sbm6uPv74Y5WUlOjss8+O+f+hUEihUN9m\nMSoqKlRRUSFJKigoUE8oXVKXJOmk4/oprY/fL5kdPHBA/fr3D7oZCdV5sFt7DhxaznlAGSoqKgq2\nQQaFw2Gv31+yMVFvpNiak5kzUpI0sH+6Mvt98fJl1/hYc3bsOSBJOtAtHX/88V7fk9QcAEgORgNW\nbm6uJCkrK0tXXXWVVq1apaFDh6q1tVU5OTlqbW1VVlZW5M9u2bIl8rXNzc2Rr4/XH267WIO+5M8a\nEB8/LBv+3Kr/tWiNJGm/f4/NIUA268293xmrK8aeltg3kAR8rDlnzn858uveXqYgAQDmGdtFcM+e\nPdq9e3fk16+99prGjBmj0tJS1dbWSpJqa2s1bdo0SVJpaanq6urU2dmpTZs2qampSZMmTfrC6/CB\n6ZaBUQ+ad7K0Ewliq97APdGTllQcAIANxmaw2tradNVVV0mSurq6dN111+mb3/ymzj//fJWVlWnh\nwoUaNmyYFi9eLEkaPXq0ysrKVFhYqIyMDC1YsCCuHQSjPzCPZvkP7Dqu/z//ye3vCrAh8IqtevOJ\nkKg1rkgLhdT9j4E4xuMAADYYC1gjRozQBx98cNjrQ4YM0bJly474NZWVlaqsrDzqa9LlSX4DM//Z\niWUGC4kSRL2BG9JC0ierkak4AAAbjB40bAMjkm4ZyAwWPMBkuTuiVzawCyQAwAb3A1bUmCSdnuR3\nXH9msOA+So070vhhAQAscz5gReO5iOQXvckFuwgCMC36c4EVDwAAG5wPWHxguiUzIy0yotzVIx3s\n7gm2QcBRYLbcHdEzWFQbAIAN7gesqF/T6Ul+oVAo5jmsvZ1MYwEwJ/rweQbkAAA2OB+w4J7jonYS\n3HOAnS7gIkZzXME5WAAA29wPWHxiOif6Oay9BCw4iNlyd6SlMYMFALDL+YDFLoLuiV4iuIclggAM\niv5YIF8BAGxwPmBFYxdBN0Rv1c4SQbiISuMOnsECANjmfMDiA9M9MUsEmcECYFDMQcPMYQEALHA/\nYEX9miWCbmAGC64LUWycEb1NOwNyAAAbnA9Y0ejyuCFmm/YDzGDBPdQad8QsEQywHQCA1OF8wOpl\nSNI5Mdu0dzKDBcAcZrAAALa5H7Cifs2yHTcwgwXXUWrcEWIGCwBgmfMBKxp9HjcwgwXAlhAzWAAA\ny5wPWHxguuf4qF0E2eQCLmIGyx08gwUAsM35gBWNTo8bjuvPNu1wG2fuuSP6GaweEhYAwALPAhad\nHhcMZJt2AJZw0DAAwDavAhbckNnvn//sOrt6AmwJcJQYy3FH9DNYwbUCAJBCCFiwLnqmkSU7AEzi\nGSwAgG3eBCxWB7ojPWbJDl0euIdy4w7OwQIA2OZNwII70mJmsOjxwD087+kOZrAAALZ5E7Do7rgj\nZlcvHsECYFCIAR0AgGX+BCxGlJ1Bhweuo9q4I40fFgDAMqcDFl1zN/FMBFzHeI47QtQbAIBlTges\naPR33JGWxgwWADtinvkMsB0AgNThT8AiYTkjegarm4AFB4UY0nFGiIOGAQCWeROw4A7OwQJgC0uS\nAQC2eROwGFF2B+dgwXXMmLuDbdoBALZ5E7DIV+7gHCy4jnLjjuifFdUGAGCDPwELzghxDhYAS9J4\nBgsAYJk3AYsRZXcwgwXnUXCcETOgQ7kBAFjgT8Ciw+OMtKh/deQrACbxDBYAwDa3Axaflk6K7vCw\nTTtcxKY67mBABwBgm9sBKwodHnekxSzZoccD9zBj7g5msAAAtvkTsOjwOIOHzgEEgWMhAAA2OB2w\n+Kh0E5tcwHWM57iDGSwAgG1OB6xodHjcQcACYEsauwgCACzzJ2CxRtAZnIMF11Fv3MEMFgDANm8C\nFtyRlhb9DBZdHriHfOWOEM98AgAs8yZg0d9xB0t2ANgSHYYpNwAAG7wJWCQsd3AOFlxHuXFH9IAO\n5QYAYIM/AQvOiN2mnR4PAHN4BgsAYJs3AYsRZXewRBCu4xksd3DuHgDANn8CFj0eZ7BNO9xHvXFF\nzK6lwTUDAJBCvAlYcEdMwGIKC4BBIXa5AABY5k3AYgLLHaGof3VMYMFF1Bt3pDGDBQCwzOmAFd03\np7/jDpYIArCFTXUAALY5HbDgpugRZbZph4sY0HEHKwRhwpw5c5SVlaUxY8ZEXvvJT36is88+W2PH\njtVVV12lnTt3Rv5fVVWV8vPzVVBQoKVLlwbRZAAWeROw2OTCHbEzWAE2BDhK1Bt3sIsgTJg9e7Ya\nGhpiXispKdGf//xnrV27VmeddZaqqqokSY2Njaqrq9O6devU0NCguXPnqru7O4hmA7DEn4AVdAMQ\nN5bsALCFZ7BgwpQpUzR48OCY1y699FJlZGRIkiZPnqzm5mZJUn19vWbMmKHMzEwNHz5c+fn5WrVq\nlfU2A7AnI+gGIPVwDhZcx4COO0JijSDse+yxx3TttddKklpaWjR58uTI/8vLy1NLS8sRv66mpkY1\nNTWSpPb2doXDYeNtDUpHRwfvz2G+v79j5U3AYsWOO9jkAoAtadG7lgbXDKSQ//iP/1BGRoZmzZrV\n56+tqKhQRUWFJKmgoEBFRUUJbl3yCIfDvD+H+f7+jpU3AYsxZXfEPHTee2iZIM+0wCX8c3VHiGc+\nYdETTzyhl156ScuWLYv828vNzdWWLVsif6a5uVm5ublBNRGABd48gwV3hEKhw0IW4JIQAzrOSGOF\nICxpaGjQfffdpxdeeEHHHXdc5PXS0lLV1dWps7NTmzZtUlNTkyZNmhRgSwGY5s0MFiPKbkkLhSJb\ntPf09iqNDisAA9hFECbMnDlT4XBY27dvV15enu6++25VVVWps7NTJSUlkg5tdPHoo49q9OjRKisr\nU2FhoTIyMrRgwQKlp6cH/A4AmORPwAq6AeiT9FBI3f8YT+7u7fXnHyJSAgM67oj+UZGvkChPP/30\nYa/ddNNNn/nnKysrVVlZabJJAJKIN0sE6fC4hSWCAGwIMYMFALDMm4AFt7CTIAAbYpYIBtgOAEDq\ncDpgRffLeejcLZyFBZcxY+4Oag0AwDanA1Y0OjxuYQYLgA1padEzWNQaAIB53gQsuCXmGaye4NoB\nHA1mzN0RYpcLAIBl3gQsujtuiR5VZgYLrmHG3B3RYZixHACADf4ELHo8TkmP+nl1E7AAGJLGjqUA\nAMu8CVhwS4hnsOAwxnPcwS6CAADbCFgIBKPKAGyg1gAAbPMmYDGi7BZ2EYTL2OTCHbGz5QE2BACQ\nMrwJWHALZ9PAZQzouIOfFQDANm8CFh+ibokZVSZhATAkZrY8wHYAAFKHPwGLJTtOSYv6l8cKQbiG\nauMOnsECANjmTcCCW9J5BguABdGz5ZQaAIANTges6M9Klgi6JY1zsOAw6o072KYdAGCb0wErGv0d\nt4Rilu3Q7YFrqDiuiFkiGFwzAAApxJuABbeksXUyAAsYzAEA2OZNwAqxZscpnIMFl1Fu3JHGM1gA\nAMv8CVhBNwB9Et1B7WHvZACGhNimHQBgmTcBC25hBgsuY0DHHWzTDgCwzZ+ARY/HKZyDBZexJNkd\n7CIIALDNm4BFd8ct6WzTDsCCELsIAgAs8yZgwS0hlgjCYQzouIODhgEAtnkTsFiy45Y0tk4GYAHP\nYAEAbPMnYAXdAPQJ52DBZYznuCONXQQBAJZ5E7DglphODwkLjgkxpOMMZrAAALY5HbCiPysZUXZL\nzDlYdHoAGMLycQCAbU4HrGiMKLslZutkhpXhGPrs7oj+UTGYAwCwwZuABbekp/EMFgDzOAcLAGCb\n8YDV3d2t8847T1dccYUkaceOHSopKdHIkSNVUlKi9vb2yJ+tqqpSfn6+CgoKtHTp0i/+5lGflowo\nuyX658U5WEgUo/UGToo91JxaAwAwz3jAeuihhzRq1KjI76urq1VcXKympiYVFxerurpaktTY2Ki6\nujqtW7dODQ0Nmjt3rrq7u003DwFJ4xwsGGCr3jCg4w5msAAAthkNWM3NzXr55Zd18803R16rr69X\neXm5JKm8vFxLliyJvD5jxgxlZmZq+PDhys/P16pVq0w2DwHiHCwkGvUGRxLiSAgAgGUZJr/5D37w\nA913333avXt35LW2tjbl5ORIkrKzs9XW1iZJamlp0eTJkyN/Li8vTy0tLYd9z5qaGtXU1EiSdu3e\npUH/eH3Pnj0Kh8Nm3khAOjo6vHtPn9ixY3/k1x+s/d9K++j/BNgaM3z++SUjE/VGiq05n/jvd97R\n/x3g3yOsPv6b/T+tXZFfd3V1eff+ovn48wMAFxkLWC+99JKysrI0YcKEzyz4oVCoz1voVlRUqKKi\nQpI0PL8gsuRj0PHHq6joa8fQ4uQTDodVVFQUdDOMWPQ/q6Vthzq7o0ePUdHo7IBblHg+//ySjal6\nI8XWnMyckZKkyRdeqNyTBhx1e5OVj/9m96xtlT5YI0lKS8/w7v1F8/HnBwAuMhawVq5cqRdeeEGv\nvPKK9u/fr127dun666/X0KFD1draqpycHLW2tiorK0uSlJubqy1btkS+vrm5Wbm5uaaah4ClcQ4W\nEsh2veERLHfELEcOrhkAgBRibI1LVVWVmpubtXnzZtXV1ekb3/iGFi1apNLSUtXW1kqSamtrNW3a\nNElSaWmp6urq1NnZqU2bNqmpqUmTJk2K+3o8dO6W6G3aeQYLx4p6g88SijlzL8CGAABShtFnsI5k\n/vz5Kisr08KFCzVs2DAtXrxYkjR69GiVlZWpsLBQGRkZWrBggdLT0+P+vnR43BK9sxfbtMMUU/UG\n7mAGCwBgm5WAVVRUFFkXPmTIEC1btuyIf66yslKVlZVxf18+LN0VYokgDDFVb6KFWCTojDRmsAAA\nlnmzDRYdHrfEdnro9cAtzJi7I8QMFgDAMn8CFh0ep8RuckG3B4AZzGABAGzzJmDBLdGdnp6eABsC\nHAXGc9wRsxw5uGYAAFKINwGLDo9bonf2YgYLgCksRwYA2OZNwGKNoFvSo/7l0eeBcyg3zogJWAG2\nAwCQOvwJWHBKGjNYcBib6rgjZpMLSg0AwAJvAhbdHbeEOAcLgAXsIggAsM2fgEXCckoa52DBYdQb\nd7CLIADANm8CFtzCg+cAbOAZLACAbd4ELAaU3RIzg8UUFhxDvXEHs+UAANucDljRn5Uh1uw4JXab\n9gAbAhwF6o07+FkBAGxzOmDBXelp7CIIwLwQM1gAAMu8CViMUbolja2T4TDqjTt4BgsAYJs/AYse\nj1PS2KYdgAUM5gAAbPMmYMEtIQ4ahsMY0HEHM1gAANu8CVghFu04hVFluIx64w6ewQIA2OZ2wIrZ\nRjCwVuAoRI8qs007AFMIwwAA29wOWHAWZ9PAafTZnZEW9SnHoeYAABu8CVj0d9ySxjbtACyImS0P\nsB0AgNThT8AiYTkl5sFzAhYcQ71xB897AgBs8yZgwS0sEYTLyFfuCMUM5gTYEABAyvAmYPEgs1tC\nnIMFwAK2aQcA2OZ0wIrZRJB85ZQ0zsGCw0IUHGdE/6SoNEiUOXPmKCsrS2PGjIm8tmPHDpWUlGjk\nyJEqKSlRe3t75P9VVVUpPz9fBQUFWrp0aRBNBmCR0wEL7uK5CAA2pLFEEAbMnj1bDQ0NMa9VV1er\nuLhYTU1NKi4uVnV1tSSpsbFRdXV1WrdunRoaGjR37lx1d3cH0WwAlngTsBhQdgvnYMFllBt3cNAw\nTJgyZYoGDx4c81p9fb3Ky8slSeXl5VqyZEnk9RkzZigzM1PDhw9Xfn6+Vq1aZb3NAOzJCLoBx+af\nn5Y8g+WW2G3aA2wIcBQY0HFHdK2h1MCktrY25eTkSJKys7PV1tYmSWppadHkyZMjfy4vL08tLS1f\n+P027+rRmfNfNtPYZNHA+3Oax+/vlGP8escDFlwVu4sg3R4AZrAcGUEIhUJH9axmTU2NampqDLQI\ngE0sEUQgOAcLLmPG3B3RPysqDUwaOnSoWltbJUmtra3KysqSJOXm5mrLli2RP9fc3Kzc3Nwjfo+K\nigqtXr1aq1evNt9gAMYwg4VARI8qs007AFNiZrCCawZSQGlpqWprazV//nzV1tZq2rRpkdevu+46\nzZs3T1u3blVTU5MmTZr0hd/vzBPStKF6qulmByYcDquoqCjoZhjD+3PbxIl3HtPXE7AQiFCIZ7Dg\nLmbM3RFithwGzJw5U+FwWNu3b1deXp7uvvtuzZ8/X2VlZVq4cKGGDRumxYsXS5JGjx6tsrIyFRYW\nKiMjQwsWLFB6enrA7wCASd4ELM6lcQtLBAHYwDNYMOHpp58+4uvLli074uuVlZWqrKw02SQAScSf\nZ7CCbgD6JGaTi57g2gHAbzFHQgTYDgBA6vAmYMEtsdu0M6wMtzBh7g5+VgAA25wOWNHdcj5E3ZLG\nM1hwGLsIuoPnPQEAtjkdsKLR3XFL7HMR9HoAmMEzWAAA2/wJWExhOSV2BoteD9xCuXFHzIY6AbYD\nAJA6vAlYcEso5hys4NoBwG+xO5YG2BAAQMrwJmAxoOwWZrDgMuqNO6IHc9hFEABggz8Bix6PUzgH\nCy5jSbI7Yn5UlBoAgAXeBCy4hXOwANjAM1gAANs8CliMKLuEc7DgMqqNOzgSAgBgmzcBixU7bqHT\nA8CGmG3ag2sGACCFeBOw4BbOwYLLGNBxx6efl6PeAABMcztgRX1O0t9xS/QMVjcdHjiGTS7cEuKw\nYQCARW4HrCj0d9wSs3UyHR4ABnEsBADAJqcDFh+T7mKbdgC2pDGgAwCwyOmAFS3EIkGnMKIMwJYQ\n9QYAYJE/AYt85ZS0qH95nIMFl1Br3JPGM1gAAIu8CVhwCzNYcBX5yj2xhw1TbwAAZnkTsBhVdkvs\nM1gBNgSA96I/HngGCwBgmj8Bi3Flp8Q+dE6PB+5gi3b3MGMOALDJm4AFt4Q4BwuAJTHnYPHMJwDA\nMH8CFoPKTmHbZLiKUuOetDRmsAAA9ngTsOj0uIVzsOAqVgi6hyWCAACbvAlYcEs6I8oALInOxFQb\nAIBpGUE3IFF48Nwt0T+uv+87qOXrPw6uMYas3dalXt/eV0gaf/rJOvG4fkG3JDBsqOOe6M+HPzRt\n00kD+gfYGnN8rDk5J31JZ2efEHQzAKBP/AlYQTcAfRK9ZGfLjn268Yl3A2yNQe/5976O65+uFbdd\nHHQzgLhFP/P5w2c+CK4hNnhYc352RWHQTQCAPnF6iSBLPdw19IQvKYNU7KS9B7r1x41/C7oZweHf\nrXPOGHxc0E3AMXi98aOgmwAAfeLPDBadHqcMHthf3zs3Uxs6T1JXj5/7Ju/42w4NHjI46GYkTOPW\nXfp4d6ckqSeFt36k1Ljn/5s2Rr/+w1+0cctHXt2Tn+ZTzdm596D+tGWnJMnTjwgAHvMnYAXdAPTZ\n+dkZ+knRxKCbYUw4HFZR0aSgm5EwP6h7X0v+tFWS1Mv8MRxSeNoJ+sW147y7Jz/Np/f333/5m66t\neUcS9QaAe5xeIgjAnlDM1voBNiRgzJYD5lFvALjMm4DFLoKAWdF3WAqvEARgQSjmMHoKDgC3+BOw\ngm4A4LuomyyVD4dmm3bAPM4uA+AybwIWALOig0Uqd3iYLAfMC8UM6ATXDgA4Gv4ELDo9gFEhhpQB\nWMOADgB3eROwWLYDmBWbr1K3y0OlAcyLHdBJ3XoDwE3eBCwAZrFk5xA21AHMY8IcgMu8CVj0eQCz\neAYLgC1s0w7AZf4ErKAbAHiOGaxDqDWAeWnR9YYhHQCOcTpgUXIBe2JGlLn7ABgUM2NOuQHgGKcD\nVnQfjyWCgFnMYP0DtQYwzla9efbZZ+N6DQD6wu2AFYVdBAGzeOj8ECoNYJfJelNVVRXXawDQFxlB\nNwCAG9g2GYAtsTNYia83r776ql555RW1tLTolltuiby+a9cuZWTQNQJwbLypIiwRBMxiF8FD2KYd\nMM/0qpTTTjtNEydO1AsvvKAJEyZEXh80aJB+8YtfGL02AP8RsADEhWewANhiut6ce+65Ovfcc3Xd\nddept7dX69evVygUUkFBgfr375/4CwJIKd4ELABmRY9h9KRwwmIwBzAvZGmb9tdff13f+9739OUv\nf1m9vb3atGmTfvWrX+nyyy83dk0A/vMoYNHrAUzi4M9DqDSAeba2aZ83b56WL1+u/Px8SdLGjRs1\ndepUAhaAY+LPLoL0egBrUjhfAbAg+jPd5Iz5oEGDIuFKkkaMGKFBgwYZux6A1ODRDBYAk0zv6uUK\nNrkAzDN9LMTzzz8vSZo4caK+9a1vqaysTKFQSM8++6zOP/98A1cEkEq8CVh0eQCzOGsOgC0hwwnr\nxRdfjPx66NCh+v3vfy9JOvXUU7V///7EXxBASnE6YEXXXAaVAbPYRfAQSg1gg9ljIR5//HED3xUA\nDvHmGSwAZsUOKKduwmIwBzDP1pLkDz/8UMXFxRozZowkae3atfr5z39u7HoAUoM3AYvlS4BZzGAB\nsMX0M1if+O53v6uqqir169dPkjR27FjV1dUZvCKAVBDXEsFt27bp17/+tTZv3qyurq7I64899pix\nhvUVo8qAWTHbtAfYjuBRbADT0iwdC7F3715NmjQp5rWMDKefngCQBOKqItOmTdPXvvY1XXLJJUpP\nTzfdJgBJiBksALbYOmj4lFNO0caNGyMDSM8995xycnKMXQ9AaogrYO3du1f33ntvn77x/v37NWXK\nFHV2durAgQOaNm2aqqurtWPHDl177bXavHmzzjzzTC1evFgnn3yyJKmqqkoLFy5Uenq6fvnLX+qy\nyy6L+3qMKQNmxRz8mWRzWDbrDbPlgHm2DhpesGCBKioqtH79euXm5mr48OFatGiRuQsCSAlxPYN1\nxRVX6JVXXunTN87MzNSbb76pDz74QGvXrtXy5cv1hz/8QdXV1SouLlZTU5OKi4tVXV0tSWpsbFRd\nXZ3WrVunhoYGzZ07V93d3XFfj7NpALOSeQbLZr2h0gDm2ao3I0aM0BtvvKFt27Zp/fr1euutt3Tm\nmWeauyCAlBDXDNZDDz2ke+65R/22/jKvAAAgAElEQVT79488CBoKhbRr167P/JpQKKTjjz9eknTw\n4EF1d3fr5JNPVn19vcLhsCSpvLxcRUVFuvfee1VfX68ZM2YoMzNTw4cPV35+vlatWqULL7zwGN8i\ngERI5mBBvQHQFw888MDn/v958+ZZagkAH8U1g7V792719PRo//792r17t3bv3v254eoT3d3dGjdu\nnLKyslRUVKQxY8aora0tsr45OztbbW1tkqSWlhadfvrpka/Ny8tTS0vL0bwnAAbY2jb5aNmqN0yW\nA+aZrjef9GVWr16tRx55RC0tLWppadGjjz6qNWvWJPx6AFJL3FvlvPDCC1qxYoUkqaioSFdcccUX\nfk16err+9Kc/aefOnbrsssu0fPnymP8fCoX6vLSvpqZGNTU1kqS9e/cp8x+vt7Q0Kxze1qfvlew6\nOjoio+8+4v255a//cyDy679s2hRgS47MRL2RYmuOJHV2HvDq5xrNt3+zn8b7c8ff9vVEfr2/szPh\nZ8rceeedkqQpU6ZozZo1GjRokCTprrvu0tSpUxN8NQCpJq6ANX/+fL377ruaNWuWpENLBleuXKmq\nqqq4LnLSSSdp6tSpWr16tYYOHarW1lbl5OSotbVVWVlZkqTc3Fxt2bIl8jXNzc3Kzc097HtVVFSo\noqJCkpQz/KzI63l5eSoqGh1Xe1wRDodVVFQUdDOM4f255b0DG6S//F9J0rBhwwNuzWdLZL2RYmtO\nZs5IfSkz06ufazTf/s1+Gu/PHVt37pN+/6YkqX//THV9wZ8/Wm1tberfv3/k9/3794/MdAPA0Ypr\nUOiVV17R66+/rjlz5mjOnDlqaGjQyy+//Llfs23bNu3cuVOStG/fPr3++usaN26cSktLVVtbK0mq\nra3VtGnTJEmlpaWqq6tTZ2enNm3apKampsPOpjhM1KoBDhoGzIo9+DO5lghaqTf/wBJBwDxb27Tf\ncMMNmjRpku666y7ddddduuCCCzR79mxj1wOQGuJeIrhz504NHjxYkvT3v//9C/98a2urysvL1dPT\no56eHl1//fUqKSnR+PHjVVZWpoULF2rYsGFavHixJGn06NEqKytTYWGhMjIytGDBgi88cyu65NLp\nAQyzdPDn0bBRbwDYY2ub9srKSl1++eX6wx/+IEl6/PHHdd555x3T96yqqtKTTz6ptLQ0nXPOOXr8\n8ce1d+/ezzwyAoB/4gpYP/3pT3Xeeefp4osvVm9vr1asWBHZ7vizjB07Vu+///5hrw8ZMkTLli07\n4tdUVlaqsrIyniYBsCx2Biu52Kw3jOUA5sXOYJk1fvx4jR8/PiHfa/PmzaqpqVFjY6MGDBigsrIy\n1dXVqbGxUcXFxZo/f76qq6tVXV3d5/NFAbgjriWCM2fO1DvvvKOrr75a06dP19tvv61rr73WdNv6\nhE4PYFbMLHGyTWEB8EpsuXGn3pxwwgnq16+f9u3bp66uLu3du1ennXaa6uvrVV5eLunQkRFLliwJ\nuKUATPrcGaz169fr7LPPjmxZmpeXJ0naunWrtm7dmrARn0RgiSBgVsySnQDbETQONQcs+NRBw67c\ndYMHD9aPf/xjnXHGGRowYIAuvfRSXXrppZ95ZMSnRe9a2t7e7s2ukEfi066XR8L7S22fG7AeeOAB\n1dTU6Ec/+tFh/y8UCunNN9801jAAySX0qQ4PAJjy6QEdUwHr9ttvP2yp3pFei9fGjRv1i1/8Qps2\nbdJJJ52ka665RosWLYr5M593ZET0rqUFBQXe7Ap5JD7tenkkvL/U9rlLBD8ZRXn11Ve1fPnymP9e\neeUVKw2MF6PKgFk2dxH88MMPVVxcrDFjxkiS1q5dq5///OdGrwkgeaRZOtj89ddfP+y1V1999ai/\n3+rVq/WVr3xFp556qvr166err75af/zjHyNHRkiKOTICgJ/iegbrK1/5SlyvBYl4BZhlcwbru9/9\nrqqqqtSvXz9JhzaxqKurM3vRODGWA5gXPWhqotw88sgjOuecc7RhwwaNHTs28t/w4cM1duzYo/6+\nBQUFeuedd7R371719vZq2bJlGjVq1GceGQHAT5+7RPCjjz5SS0uL9u3bp/fffz8yirRr1y7t3bvX\nSgMBJAfTHZ5oe/fuPexcqoyMuE+VMIqABZhnek+d6667Tpdffrl++tOfxuyKPGjQoMiRNEdj3Lhx\nuuGGGzRx4kSlpaXpvPPOU0VFhTo6Oo54ZAQAP31uj2Xp0qV64okn1NzcrHnz5kVeHzRokO655x7j\njesTOj2AUTZnsE455RRt3LgxEuqee+65yAPiAPwXMrxE8MQTT9SJJ56op59+Wt3d3Wpra1NXV5c6\nOjrU0dGhM84446i/9+23367bb7895rXMzMzPPDICgH8+N2CVl5ervLxcv/vd7zR9+nRbbToqIRIW\nYFTsQ+dmE9aCBQtUUVGh9evXKzc3V8OHDz/sQfGgUGsA82ztWvpf//VfuuuuuzR06FClpR16aiIU\nCmnt2rUGrwrAd3GtuSkqKtItt9yit956S6FQSBdddJHuuOMODRkyxHT7ACSJ2HOwzF5rxIgReuON\nN7Rnzx719PRo0KBBZi8IILlYqjcPPvigNmzYQH8GQELFFbBmzJihKVOm6He/+50k6amnntK1116r\nN954w2jj+oLnIgCzLOYr7dy5U7/5zW+0efNmdXV1RV7/5S9/afjKX4xaA5gXs0TQ4HVOP/10nXji\niQavACAVxRWwWltb9bOf/Szy+3//93/XM888Y6xR8YouuvR5ALOiOzw9PWYj1re+9S1NnjxZ55xz\nTmTZTrKg1gDmxW5yYa7ejBgxQkVFRZo6daoyMzMjr0c/dw4AfRVXwLr00ktVV1ensrIySYceOL/s\nssuMNgxAcrH1TIQk7d+/Xw888IDhqwBIVrZ2LT3jjDN0xhln6MCBAzpw4IDBKwFIJXEFrF//+td6\n8MEHdf3110uSenp6NHDgQP3qV79SKBTSrl27jDYyHizbAcyyuYvgrFmz9Otf/1pXXHFFzKjysWyf\nnCgcag6YZ3qb9k/ceeedkg4dDXHccceZuxCAlBLX2pvdu3erp6dHXV1d6urqUk9Pj3bv3q3du3cn\nRbiS2NkLsMn0LoKZmZn68Y9/rAsvvFATJkzQhAkTNHHiRKPXBJA8Yp/BMldv3n77bRUWFurss8+W\nJH3wwQeaO3eusesBSA1xn9y5du3awx44v/rqq4006mgwqAyYFbNkx/AM1gMPPKCNGzfqlFNOMXuh\no0CpAcyLWZJssN784Ac/0NKlS1VaWipJOvfcc7VixQpzFwSQEuIKWHPmzNHatWs1evTomHMikilg\nATDLZrDIz89P3uU6JCzAOFu7CEqHdhKMlp6ebviKAHwXV8B655131NjYaLotx4Q+D2BW7DNYZrs8\nAwcO1Lhx43TxxRfHPIOVDNu0A7DLZL05/fTT9cc//lGhUEgHDx7UQw89pFGjRhm7HoDUEFfAmjRp\nkhobG1VYWGi6PUePNYKAUTbPwbryyit15ZVXGr7K0aHSAOalWVqS/Oijj+rWW29VS0uLcnNzdeml\nl2rBggXmLgggJcQVsGbPnq3JkycrJydHmZmZ6u3tVSgU0tq1a023D0CSsPkMVnl5uQ4cOKAPP/xQ\nklRQUKB+/fqZvSiApGFrieApp5yip556yuAVAKSiuALWzTffrEWLFiXloZ+fYFQZMMvWrl6SFA6H\nVV5erjPPPFO9vb3asmWLamtrNWXKFKPXjQfbtAPm2TpoeNOmTXr44YcP28TrhRdeMHZNAP6LK2Cd\neuqpkR12khV9HsAsmzNYP/rRj/Taa6+poKBAkvThhx9q5syZeu+998xeOA6UGsA8WwcNX3nllbrp\nppv07W9/O2kHkAG4J66Add555+m6667Tt7/97ZgHztlFEEgdNp/BOnjwYCRcSdJZZ52lgwcPGr4q\ngGRh66DhzMxM3XLLLeYuACAlxRWw9u3bp8zMTL322muR15Jhm/bomstBw4BZsbsImr3WxIkTdfPN\nN+v666+XJD311FNJc9Aws+WAebbus1tuuUV33XWXLrvsspgB5PHjx9tpAAAvxRWwHn/8cdPtOGZ0\negCzQhbnsB555BEtWLAgsi371772Nc2dO9foNQEkD1vPOv75z3/Wk08+qeXLl8ec8/nmm29auT4A\nP8UVsG688cYjFrvHHnss4Q3qE9PrlABE2JzB6urq0q233qp58+ZJkrq7u9XZ2Wn2onFithzwx3PP\nPadNmzapf//+QTcFgEfieqLziiuu0NSpUzV16lQVFxdr165dOv744023rU/o8gBmRd9jPYYTVnFx\nsfbt2xf5/b59+3TJJZcYvWa8mC0H7LBxr40ZM0Y7d+40fyEAKSWuGazp06fH/H7mzJm66KKLjDTo\naNHpAcyyOYO1f//+mEGc448/Xnv37jV7UQBJJSTzC1V27typs88+W+eff37MM1hs0w7gWMQVsD6t\nqalJH3/8caLbAiCJRS+NM93pGThwoNasWRN50Py9997TgAEDDF8VQDIJhULGR3Puvvtuo98fQGqK\nK2ANGjQo5hms7Oxs3XvvvcYadTQ4/BMwzOIM1oMPPqhrrrlGp512mnp7e/XRRx/pmWeeMXtRAEnF\nxqf617/+dQtXAZBq4gpYu3fvNt0OAEkudg9Bswnr/PPP1/r167VhwwZJUkFBgfr162f0mvFiMAew\nw8at9vzzz+v222/Xxx9/rN7eXvX29ioUCmnXrl3mLw7AW3EFrJUrV2rcuHEaOHCgFi1apDVr1ujW\nW2/VsGHDTLcPQJKICRYWdvB89913tXnzZnV1dWnNmjWSpBtuuMH8hb8A8QqwI2ThKazbbrtNL774\nokaNGmX0OgBSS1y7CP7rv/6rjjvuOH3wwQe6//779eUvfzkpOjrRGFQGzLJ3Cpb0L//yL/rxj3+s\nt956S++++67effddrV692vBVASQVC5/rQ4cOJVwBSLi4ZrAyMjIUCoVUX1+v73//+7rpppu0cOFC\n023rE86mAcyK3UXQbMRavXq1Ghsbk3I5XhI2CfBSmoV7beLEibr22mt15ZVXxuwiePXVV5u/OABv\nxb3JRVVVlRYtWqQVK1aop6dHBw8eNN02AEnE5grBMWPG6KOPPlJOTo7hKwFIVjYGTnft2qXjjjtO\nr7322j+vGwoRsAAck7gC1jPPPKPf/va3WrhwobKzs/XXv/5VP/nJT0y37QtFd/IYVQbMitmm3XDC\n2r59uwoLCzVp0qSkO5uGWgPYYfpe6+7u1tixY/XDH/7Q7IUApJy4AlZ2drbmzZsX+f0ZZ5yRfM9g\nBd0AwHM2Z7Duuusuw1c4eixHBuwwfaelp6fr6aefJmABSLi4AhbbmAKIfh7K9DNYnE0DwMYzmF/9\n6lf1/e9/X9dee60GDhwYef2TQ84B4GjEFbBc2MaUZTuAWTZ2Ebzooov01ltvHXa4eTIN6lBrADts\n3Gp/+tOfJEl33HHHP68bCunNN9+0cHUAvoorYLmwjSnLdgCzQhYS1ltvvSWJw80ByErCWr58ufmL\nAEg5cQUstjEFED2I0WN6l4skxlAOYIeNe62trU3/9m//pq1bt+rVV19VY2Oj3n77bd10000Wrg7A\nV3EdNBy9jemLL76oF198US+99JLptvUJy3YAs2LPwQquHYGj2ABW2HgGa/bs2brsssu0detWSdJZ\nZ52lBx980Ph1Afgtrhmsxx9/3HQ7ACS52BWCqZywANhgYyxj+/btKisrU1VVlSQpIyND6enp5i8M\nwGtxzWA1NzfrqquuUlZWlrKysjR9+nQ1NzebbhuAJMIM1iHMXwF22LjXBg4cqL/97W+R2bJ33nlH\nJ554ooUrA/BZXAHrxhtvVGlpqbZu3aqtW7fq29/+tm688UbTbesTG0sJgNQWtatfgK0IGqUGsMPG\n5/r999+v0tJSbdy4UV/96ld1ww036OGHHzZ+XQB+i2uJ4LZt22IC1ezZs5NujTJ9HsAsZrAA2GTj\nc33ChAn6/e9/rw0bNqi3t1cFBQXq16+fhSsD8FlcM1hDhgzRokWL1N3dre7ubi1atEhDhgwx3TYA\nSSS2s5O6CYvBHMAOG7PFY8eO1X333acvfelLGjNmDOEKQELEFbAee+wxLV68WNnZ2crJydFzzz2n\nJ554wnDTvlh0F49lO4BZsQf/BtgQACnBxhLBF198URkZGSorK9P555+v//zP/9Rf//pX49cF4Le4\nAtYdd9yh2tpabdu2TR9//LEee+wx3Xnnnabb1ifkK8AsC+cMO4HnPQE7bNxpw4YN02233ab33ntP\nv/3tb7V27VoNHz7cwpUB+CyuZ7DWrl2rk08+OfL7wYMH6/333zfWKADJJ/YZrNSNWMQrwA5bYxn/\n8z//o2eeeUbPPPOM0tPTdd9999m5MABvxRWwenp61N7eHglZO3bsUFdXl9GG9RWjyoBZMQEruGYA\nSBEhC8MZF1xwgQ4ePKhrrrlGzz77rEaMGGH8mgD8F1fA+tGPfqQLL7xQ11xzjSTp2WefVWVlpdGG\nxSWql0e+AsyK7uyk8AQWtQawxMa99pvf/EYFBQXmLwQgpcQVsG644QZNnDhRb775piTp+eefV2Fh\nodGGAUguzGABsMnGWEZ2drbmzZunFStWSJK+/vWv64477uCwYQDHJK6AJUmFhYVJHaoYVAbMit1F\nMHUjlo1lSwDsLP2fM2eOxowZo8WLF0uSnnzySd144416/vnnjV8bgL/iDlhJj3U7gFHcYf/AXwTg\njY0bN+p3v/td5Pd33nmnxo0bF2CLAPggrm3aASB6DKMnhWewANhhY9x0wIABeuuttyK/X7lypQYM\nGGD+wgC85s0MFoPKgFlscnEItQaww0bAeuSRR1ReXq6///3v6u3t1eDBg/XEE0+YvzAAr/kTsOj1\nAEbFnoMVXDsApAYbzzuOGzdOH3zwgXbt2iVJOuGEE4xfE4D/vAlYAMyK7ur0pvA+ggzmAHbYuNce\neOCBw1478cQTNWHCBJ7FAnDUnH4GK7qLx85egGHMYEmi1gC22LjTVq9erUcffVQtLS1qaWnRr371\nKzU0NOi73/2u7rvvPgstAOAjb2awGFUGzIp5BivAdgBIDTa2aW9ubtaaNWt0/PHHS5LuvvtuTZ06\nVStWrNCECRN02223GW8DAP84PYMFwJ5Q7BrBlMVgDmCHjXvt448/VmZmZuT3/fr1U1tbmwYMGBDz\nOgD0hT8zWEE3APAcz2ABsMnG5/qsWbN0wQUXaNq0aZKkF198Udddd5327NmjwsJCCy0A4CN/AhYJ\nCzAqerlOSj+DRa0BrLCxRPBnP/uZLr/8cq1cuVKS9Oijj2rixImSpKeeesr49QH4yZuABcCsmG3a\ng2tG4NjkArDD1p02ceLESKgCgETw5hksOj2AWTFLBFN5CguAFcwWA3CVNwGLfAWYxQzWIXT6ADsY\nOAXgKn8CFgCjeAYLgE0MZgBwlTcBizoMmMUu7QDwxXbu3KnvfOc7OvvsszVq1Ci9/fbb2rFjh0pK\nSjRy5EiVlJSovb096GYCMMifgMVQF2BU7AxW6kYsag1gh6v32q233qpvfvObWr9+vT744AONGjVK\n1dXVKi4uVlNTk4qLi1VdXR10MwEY5E3AAmBW7CYXgTUDQIpwMV79/e9/14oVK3TTTTdJkvr376+T\nTjpJ9fX1Ki8vlySVl5dryZIlQTYTgGFOB6zoPp6LhRhwSewmF6mbsKg1gB0uTmBt2rRJp556qm68\n8Uadd955uvnmm7Vnzx61tbUpJydHkpSdna22traAWwrAJG/OwXKxEAMuid7RixksAKa5+Lne1dWl\nNWvW6OGHH9YFF1ygW2+99bDlgKFQ6DOXP9bU1KimpkaS1N7ernA4bLrJgeno6OD9Ocz393esvAlY\nAMyKmcFK4YDlYqcPcJGL27Tn5eUpLy9PF1xwgSTpO9/5jqqrqzV06FC1trYqJydHra2tysrKOuLX\nV1RUqKKiQpJUUFCgoqIiW023LhwO8/4c5vv7O1ZOLxGMRqcHsCeF85WDXT7ATS5+rmdnZ+v000/X\nhg0bJEnLli1TYWGhSktLVVtbK0mqra3VtGnTgmwmAMO8mcFycaQLcEnsDFYqRywANrj6qf7www9r\n1qxZOnDggEaMGKHHH39cPT09Kisr08KFCzVs2DAtXrw46GYCMMibgAXALAYxDnF162jANa7ea+PG\njdPq1asPe33ZsmUBtAZAENxeIhg1iO5oHQacwTNYh1BqADv4XAfgKrcDFgBr2KYdgE3kKwCuImAB\niAvbtB/CqDpgh6tLBAHAm4BFIQbMip3BAgCz+FQH4Cp/AlbQDQA8F32PpfYuglQbwAbGTQG4ypuA\nBcCs6FnilI5XdPoAK9i5FICrnA5Y0Z08Oj2AWTH3WConLAB28LkOwFFOB6xojHQBZkXfYT0pvESQ\nSgPYwb0GwFXeBCwAZrFEEIBNrEwB4CpvAhaFGDArdpOLwJoROGoNYAcrUwC4yljA2rJliy6++GIV\nFhZq9OjReuihhyRJO3bsUElJiUaOHKmSkhK1t7dHvqaqqkr5+fkqKCjQ0qVL+3Q9yjBgVjIfNGyz\n3tDpA+xgMAOAq4wFrIyMDN1///1qbGzUO++8owULFqixsVHV1dUqLi5WU1OTiouLVV1dLUlqbGxU\nXV2d1q1bp4aGBs2dO1fd3d1xX49CDJiVzAcN2643AMzjcx2Aq4wFrJycHI0fP16SNGjQII0aNUot\nLS2qr69XeXm5JKm8vFxLliyRJNXX12vGjBnKzMzU8OHDlZ+fr1WrVplqHoA+ipnBSrKAZbPe0OkD\n7GC2GICrMmxcZPPmzXr//fd1wQUXqK2tTTk5OZKk7OxstbW1SZJaWlo0efLkyNfk5eWppaXlsO9V\nU1OjmpoaSdKBgwcir69bt05f2r7B5NuwrqOjQ+FwOOhmGMP7c8u2vT2RX+/fvz9puz6JrDdSbM2R\npG3btnn1c43m27/ZT+P9uWVn+/6gmwAAR8V4wOro6ND06dP14IMP6oQTToj5f6FQKGZnsnhUVFSo\noqJCknTy6WdFXh8zZoyKRmcfe4OTSDgcVlFRUdDNMIb355bm9r3SiuWSpMzMTB34gj8fhETXGym2\n5mTmjFRW1qkqKpqQkPYmG9/+zX4a788tj/1llfS3bUE3AwD6zOguggcPHtT06dM1a9YsXX311ZKk\noUOHqrW1VZLU2tqqrKwsSVJubq62bNkS+drm5mbl5uaabB6APkj2bdpt1RuWLQF2cKcBcJWxgNXb\n26ubbrpJo0aN0rx58yKvl5aWqra2VpJUW1uradOmRV6vq6tTZ2enNm3apKamJk2aNOnzrxH1awox\nYFYyb9Nuo94AsIvnHQG4ytgSwZUrV+rJJ5/UOeeco3HjxkmS7rnnHs2fP19lZWVauHChhg0bpsWL\nF0uSRo8erbKyMhUWFiojI0MLFixQenp63Nc7mqU/AOKXzNu0W603lBrACm41AK4yFrAuuugi9X7G\nMPeyZcuO+HplZaUqKytNNQnAMUjmbdqpN4B/GDgF4Cqjz2DZRBkGzIqdwUpd1BrADu41AK7yJ2BR\niQGjkvkcLJsYVQfs4FYD4CpvAhYAs2KXCKZwwgJgCQkLgJu8CViMdAFmsUTwEEoNYAef6wBc5U/A\notsDGBW7TXsqRywANvCpDsBV3gQsAGYl+0HDtjCqDtjBvQbAVf4ELAoxYFQyHzRsE6UGsIOVKQBc\n5XbAiurkUYYBs2J3EUzhhAXACmawALjK7YAFwJqYXQQDbEfQ2KYdsCONew2Ao5wOWL1R3Tw6PYBh\nMWsEA2sFgFTBxzoARzkdsKJRhwGz2Kb9EGoNYAf3GgBXeROwAJjFNu3/QK8PsIKVKQBc5U3Aog4D\nZrFNOwCb+FgH4Cp/AhalGDCKbdoPodYAdjBwCsBV3gQsAGbFPoOVugmLTh9gB7caAFd5E7Do9ABm\nxWzTnrr5CoAlPIMFwFX+BKygGwB4jl0ED6HWAHZwrwFwlTcBC4BZMQGLKSwAppGwADjKn4BFIQaM\nYongIaxaAuxgQxkArvImYFGIAbNYIngItQawg8EMAK7yJmABMIuDhgHYRL4C4CqnA1Z0F4+RLsAs\nDho+hFoD2MG9BsBVTgesaNRhwCwOGgZgE8txAbjKm4AFwCxGkw/h7wGwI40eCgBHeVO+OJAQMIt7\n7BP8PQB2cK8BcJNHASvoFgAAgEThcx2Aq7wJWADMo8PD3wFgC7caAFd5E7AoxIB53GcAbGEwA4Cr\n/AlYFGLAOJ7DImQCtrCLIABXeROwAJhHd4fBHMAW7jUArvIoYFGJAdPo8ACwhXIDwFXeBCw6foB5\nLBFk2RJgC/UGgKu8CVgAzKO7AwAA8PmcDli9Ub+m4weYx4AyfweALdxrAFzldMCKTlgsJQDMY3kc\ngzmALdQbAK5yO2ABsIpxDAC2UG8AuMqbgEUdBszjPmO2HLCFOw2Aq/wJWFRiwDjCBQBb0tKoNwDc\n5E3AAmAe3R0AtlBvALjKm4DFw7CABdxmzJYDtnCvAXCUPwGLQgwYx20GwBYGTgG4ypuABcA8nsGi\n0wfYQrkB4CoCFoC40eEBYIvL5aa7u1vnnXeerrjiCknSjh07VFJSopEjR6qkpETt7e0BtxCASd4E\nLDp+gHncZtQawBaX77WHHnpIo0aNivy+urpaxcXFampqUnFxsaqrqwNsHQDTvAlYAMxjiSAhE7DF\n1eW4zc3Nevnll3XzzTdHXquvr1d5ebkkqby8XEuWLAmqeQAsyAi6AceiN+rXrhZiwCXcZQBscXU8\n5wc/+IHuu+8+7d69O/JaW1ubcnJyJEnZ2dlqa2s74tfW1NSopqZGktTe3q5wOGy8vUHp6Ojg/TnM\n9/d3rJwOWNFcLcSAS5jBotYAtrh4q7300kvKysrShAkTPrPzGQqFPrOWVlRUqKKiQpJUUFCgoqIi\nQy0NXjgc5v05zPf3d6y8CRhkk0gAABVESURBVFgAzCNcEDIBaxy811auXKkXXnhBr7zyivbv369d\nu3bp+uuv19ChQ9Xa2qqcnBy1trYqKysr6KYCMMibZ7AcrMOAc7jNANjiYr2pqqpSc3OzNm/erLq6\nOn3jG9/QokWLVFpaqtraWklSbW2tpk2bFnBLAZjkT8ByshQDbmEgw81OH+Ain+rN/Pnz9frrr2vk\nyJF64403NH/+/KCbBMAglggCiBsDGQBscb3eFBUVRZ5RGTJkiJYtWxZsgwBY488Mltt1GHAC95mY\nwgIsSeNeA+AofwJW0A0AUgD3mfuj6oArGNAB4CpvAhYA89hBD4At1BsArvImYFGHAdhArQEAAJ/H\nm4DF4iXAPMIFAFuoNwBc5VHAAmAaHR6GcgBbeN4RgKucDli9Ub+m4weYR4eHWgPYwr0GwFVOB6xo\n1GHAPDo8AGyh3ABwlTcBC4B5dHiYxQNsYUAHgKu8CVhs5wqYx30GwBYGMwC4yu2AFfUQFmUYMI98\nxd8BYAv3GgBXuR2wolCIAfO4zfg7AAAAn8+bgAXAPJYIArCFegPAVd4ELNZqA+Zxl4npcsAS7jQA\nrvInYFGJAeO4zwDYkka9AeAobwIWAPOYKWZUHbCFJYIAXEXAAhA3+jv8HQC2cK8BcJXTAStql3YK\nMQAAHuFjHYCrnA5YAOxiyQ7LJAFrqDcAHOVNwKLjB5jHXQbAFuoNAFf5E7CCbgCQAhjH4O8AsIV7\nDYCrvAlYAMyjw8NgDmALy3EBuMqbgEXHDzCPDg8AW/hcB+AqfwIWHT/AODo8/B0AtnCrAXCVNwEL\ngHl0eADYwmAGAFd5E7AoxIAF3GjsWApYwsoUAK7yJ2AF3QAgBaRxowGwhXoDwFHeBCwA5tHfAWAL\n9QaAq/wJWFRiwDiWx7FKErCFegPAVU4HrN6oX7NWGzCPuwyALSxJBuAqpwMWALsYUGYwB7CFegPA\nVd4ELAoxYB7hgloD2EK9AeAqfwJW0A0AUgE3GgBLGMwA4CpvAhYA8+jv8HcAAAA+nzcBi92GAPO4\nzfg7AGzhcx2Aq4wFrDlz5igrK0tjxoyJvLZjxw6VlJRo5MiRKikpUXt7e+T/VVVVKT8/XwUFBVq6\ndGmfr0cZBsxL5mcibNccAGYlb7UBgM9nLGDNnj1bDQ0NMa9VV1eruLhYTU1NKi4uVnV1tSSpsbFR\ndXV1WrdunRoaGjR37lx1d3ebahqAo5TMA8q2ak4yh0zAJ8lcbwDg8xgLWFOmTNHgwYNjXquvr1d5\nebkkqby8XEuWLIm8PmPGDGVmZmr48OHKz8/XqlWrvvgiUQdhUYgB85L5PrNScwBYw2AGAFdl2LxY\nW1ubcnJyJEnZ2dlqa2uTJLW0tGjy5MmRP5eXl6eWlpYjfo+amhrV1NRIknp7/5mwVq5cqYH9/CrG\nHR0dCofDQTfDGN6fe3a27wu6CX2S6JojSRs3blS4568GWx0cH//NRuP9uaXxo66gmwAAR8VqwIoW\nCoWO6gHWiooKVVRUSJIG5IyMvH7RVy/Sicf1S1j7kkE4HFZRUVHQzTCG9+eehRv/W/rb9qCbcVQS\nUXMyc0YqP//LKvraiEQ3Lyn4+G82Gu/PLfv+d6v0pzVBNwMA+szqLoJDhw5Va2urJKm1tVVZWVmS\npNzcXG3ZsiXy55qbm5Wbm9u3b+7X5BWABDBacwAYlcxLkgHg81gNWKWlpaqtrZUk1dbWatq0aZHX\n6+rq1NnZqU2bNqmpqUmTJk36wu8X9QgWhRiwIM2xGy3RNUdi62jAHu41AG4ytkRw5syZCofD2r59\nu/Ly8nT33Xdr/vz5Kisr08KFCzVs2DAtXrxYkjR69GiVlZWpsLBQGRkZWrBggdLT0001DcBRSuZs\nQc0B/JLM9QYAPo+xgPX0008f8fVly5Yd8fXKykpVVlYe9fWow4B5yXyf2ao5yfx3APjEtRlzAPiE\n1SWCJrFsBzCP+4xRdcAWbjUArvImYAEwjw4PAFsYzADgKm8CFnUYMI8OD7UGsIV6A8BV/gQsCjFg\nATcaADtC1BsAjvImYAEwj4EMnkMDrOFWA+AobwIWI12AedxlhEzAFm41AK7yJ2BRiQHjuM8A2MJs\nMQBXeROwAJjHTDGj6oAt3GsAXEXAAhA3BpQB2EK9AeAqpwNWb9SvKcSAedxn4i8BsIQZcwCucjpg\nAbCLDg/LlgBbGMsA4CpvAhYdP8ACbjMAllBuALjKn4BFJQaMS+NGo9YAtnCvAXCUNwELgHn0dwDY\nwoAOAFd5E7Aow4B59HdYjgzYwp0GwFX+BCx6foBx3GWETMAWPtcBuMqbgAXAPDo8AGxxsdxs2bJF\nF198sQoLCzV69Gg99NBDkqQdO3aopKREI0eOVElJidrb2wNuKQCTvAlYDtZhwDncZ/wdALa4eK9l\nZGTo/vvvV2Njo9555x0tWLBAjY2Nqq6uVnFxsZqamlRcXKzq6uqgmwrAIH8ClouVGHAN9xkAS1z8\nXM/JydH48eMlSYMGDdKoUaPU0tKi+vp6lZeXS5LKy8u1ZMmSIJsJwLCMoBsAwB1s8OBmpw9wk9s3\n2+bNm/X+++/rggsuUFtbm3JyciRJ2dnZamtrO+LX1NTUqKamRpLU3t6ucDhsq7nWdXR08P4c5vv7\nO1beBCyeDQHM4zYjZAK2uFxvOjo6NH36dD344IM64YQTYv5fKBT6zD5LRUWFKioqJEkFBQUqKioy\n3dTAhMNh3p/DfH9/x8qbJYIAzHO4vwPAMa7Wm4MHD2r69OmaNWuWrr76aknS0P/X3v2FtlnvcRz/\nJE1bGHTsTGxXsjEXU2vXJg1u0wtl0OmOHHAZOpWCwqDDHWSIsjLZnYMz7MALFf/BdGiRadELl4t1\nU/DCv0hxWGWdO5auky4ts+pkf8DOzt+5OJhtusY0TX9Pfr+8X1dr2sDv26f57Pk8z5MnDQ2amJiQ\nJE1MTKi+vj7IJQKYZxQsAAVz+YhyyfA7AKxw8coUY4y2bNmilpYWbd++Pfd4Op1Wb2+vJKm3t1cb\nN24MaokALPDiEkEHMxhwEpfH0a8AW1x8rX322Wd68803lUgklEqlJElPP/20du7cqQcffFD79u3T\n8uXL9c477wS8UgDzyY+CFfQCgArBwQwAtriYN3fccYeMMdf83ocffmh5NQCCwiWCAArm4g5Pqbl4\n2RLgIs6YA3CVFwWLHR7AFl5rAOzgv3YArvKjYAW9AKBChHmxkTcAACAvLwoWADs4oszvALAlzIsN\ngKO8KFhkMGAH74kAYAv/twNwlR8Fi50+wAp2ePgdALbwWgPgKi8KFgA72N8BYAsHTwG4yo+CRQYD\nVnDHTnb6AFuIGwCu8qJgkcEAbGGnD7CDlxoAV3lRsADYQbkAYAt5A8BVXhQsQhiwg8vjANhD3gBw\nkx8FixAGrOBgBgBbyBsArvKiYAGwg/0dbvQB2MIrDYCrvChY7O8AdvBaY6cPsIWDGQBc5UfBCnoB\nQIVghweALaQNAFf5UbDY6QOs4JXGWTzAFl5rAFzlR8EKegFApeDFBsASbmAFwFWRoBdQCtURL3oi\nUPZqqsKqjYQr+sgyO32AHeGwKj5vALjJi2ZSXUX6AjZ0/7NZ/939Lx3/z7+CXkpg2NkD7Fj6jwUV\nnzcA3ORJwfJiDAAAAACO86KZ1FCwAFjCCSwAAJCPF82EM1gAAAAAyoEXzaQ6wjFlAHbwHiwAAJCP\nHwWLM1gArKFhAQCAmXnRTChYAAAAAMqBF82Em1wAsIVLBAEAQD5eNBM+BwsAAABAOfCkYHkxBgAH\ncDgHAADk40UzqY54MQYAB4S4RhAAAOThRTPhPVgAAAAAyoEXzYT3YAGwhbQBAAD5eFKwvBgDgAO4\nQhAAAOTjRTOp4T1YAAAAAMqAF82E92ABsIUzWAAAIB8vmgmXCAIAAAAoB140EwoWAFtC3OYCAADk\n4UUzqY6wwwPAEuIGAADk4UXB4j1YAAAAAMqBF82ESwQB2MIJLAAAkI8XzYSCBQAAAKAceNFMqqs4\npgzAjhD3aQcAAHl4UbD4oGEAtlCvAABAPl40Ey4RBAAAAFAOvGgmFCwAtnCFIAAAyMeLZsJ7sAAA\nAACUAy8KFp+DBcCWEO/CAgAAeXjRTKq5yQUAS7hEEAAA5ONFM+E9WAAAAADKgRfNhEsEAdjCCSwA\nAJCPF82kJsIuDwAAAIDgeVGwuEQQgDUczwEAAHl40UwoWABs4S6CAAAgHy+aCQULAAAAQDnwoplw\nkwsAtnCbdgAAkI8XzaSam1wAAAAAKAN+FCzOYAGwhMM5AAAgHy+aSSTMLg8AO0JcIwgAAPJwvmDV\nVIXZ4QEAAABQFpwvWNVVlCsA9nA8BwAA5ON+wYo4PwIAAAAAT5RdOzl8+LCam5sVj8e1Z8+ev/15\nbnABoFizzRuJm1wAKE4xeQPATWXVTi5duqRt27bp0KFDOnbsmN5++20dO3Ys73P4DCwAxSgmbyQu\nEQQwe8XmDQA3lVU7GRgYUDweVywWU01NjTo7O5XJZPI+h/dgAShGMXkDAMUgb4DKUlYFK5vNatmy\nZbmvly5dqmw2m/c5XCIIoBjF5M3/cVAHwOwUnzcAXBQJegGztXfvXu3du1eSFDozpjNvdWv1WwEv\nap5MTk7q+uuvD3oZ84b53Hb8+PGgl2DFnzPn35vWB7yi+eP73yzzuasS8+b777/X6tWrA17R/PH5\n71ViPtfNNXPKqmBFo1GNjY3lvj516pSi0ehVP7N161Zt3bpVkrR69Wp9+eWXVtdoE/O5rRLmc1kh\neSOROT5hPneRN/5hPrdVwnxzUVbX161Zs0bDw8MaHR3VxYsX1dfXp3Q6HfSyAHiIvAFgC3kDVJay\nOoMViUT04osv6u6779alS5fU1dWl1tbWoJcFwEPkDQBbyBugslTt2rVrV9CLuFJTU5Mee+wxPf74\n41q7du3f/vyqVassrCo4zOc25itvs80byf2Z/w7zuc3n+Vyfjbz5K+ZzG/PNLGSMMSVcCwAAAABU\nrLJ6DxYAAAAAuMzZgnX48GE1NzcrHo9rz549QS+nJG644QYlEgmlUqnc3Ut+/vlnrV+/Xk1NTVq/\nfr3OnDkT8CoL19XVpfr6erW1teUeyzdPT0+P4vG4mpub9f777wex5Fm51ny7du1SNBpVKpVSKpVS\nf39/7nsuzTc2NqaOjg6tXLlSra2tev755yX5tf1mg7wpf+SNu3kjkTl/RuaUPzLH3cyxkjfGQdPT\n0yYWi5mRkREzNTVlksmkGRoaCnpZc7Z8+XIzOTl51WM7duwwPT09xhhjenp6zJNPPhnE0ory0Ucf\nmSNHjpjW1tbcYzPNMzQ0ZJLJpPn111/NiRMnTCwWM9PT04Gsu1DXmu+pp54yzzzzzF9+1rX5xsfH\nzZEjR4wxxpw9e9Y0NTWZoaEhr7ZfocgbN5A3l7k4H5lzGZnjBjLnMtfms5E3Tp7BGhgYUDweVywW\nU01NjTo7O5XJZIJe1rzIZDLavHmzJGnz5s06cOBAwCsq3Nq1a7V48eKrHptpnkwmo87OTtXW1mrF\nihWKx+MaGBiwvubZuNZ8M3FtvsbGRt1yyy2SpLq6OrW0tCibzXq1/QpF3riBvLnMxfnInMvIHDeQ\nOZe5Np+NvHGyYGWzWS1btiz39dKlS5XNZgNcUWmEQiHdddddWrVqVe6T3E+fPq3GxkZJ0pIlS3T6\n9OkglzhnM83j0zZ94YUXlEwm1dXVlTu97PJ8J0+e1FdffaXbbrutIrbfn/k6G3njxzb1LW8kMsfX\n2cgcP7apb5kzX3njZMHy1aeffqrBwUEdOnRIL730kj7++OOrvh8KhRQKhQJaXen5No8kPfroozpx\n4oQGBwfV2Nio7u7uoJc0J+fPn9emTZv03HPPaeHChVd9z8ftV0nIG/f5ljcSmeMzMsd9vmXOfOaN\nkwUrGo1qbGws9/WpU6cUjUYDXFFp/DFDfX297r33Xg0MDKihoUETExOSpImJCdXX1we5xDmbaR5f\ntmlDQ4OqqqoUDof1yCOP5E4huzjfb7/9pk2bNumhhx7SfffdJ8n/7Xctvs5G3ri/TX3KG4nM+YOv\ns5E57m9TnzJnvvPGyYK1Zs0aDQ8Pa3R0VBcvXlRfX5/S6XTQy5qTCxcu6Ny5c7l/f/DBB2pra1M6\nnVZvb68kqbe3Vxs3bgxymXM20zzpdFp9fX2amprS6OiohoeHdeuttwa51KL88cKUpPfeey939x3X\n5jPGaMuWLWppadH27dtzj/u+/a6FvHGX73+vvuSNROZcicxxl+9/r75kjpW8mcebdMyrgwcPmqam\nJhOLxczu3buDXs6cjYyMmGQyaZLJpFm5cmVuph9//NGsW7fOxONxc+edd5qffvop4JUWrrOz0yxZ\nssREIhETjUbNa6+9lnee3bt3m1gsZm666SbT398f4MoLc635Hn74YdPW1mYSiYTZsGGDGR8fz/28\nS/N98sknRpJJJBKmvb3dtLe3m4MHD3q1/WaDvCl/5I27eWMMmfNnZE75I3PczRwbeRMyxpj56YcA\nAAAAUFmcvEQQAAAAAMoRBQsAAAAASoSCBQAAAAAlQsECAAAAgBKhYAEAAABAiVCwUFK//PKLXn75\nZUnS+Pi47r///oBXBMBnZA4AW8gbFIrbtKOkTp48qXvuuUdHjx4NeikAKgCZA8AW8gaFigS9APhl\n586dGhkZUSqVUlNTk7799lsdPXpUb7zxhg4cOKALFy5oeHhY3d3dmpqa0v79+1VbW6v+/n4tXrxY\nIyMj2rZtmyYnJ7VgwQK9+uqruvnmm4MeC0CZInMA2ELeoGDz9SnJqEyjo6OmtbX1L/9+/fXXzY03\n3mjOnj1rfvjhB1NXV2deeeUVY4wxTzzxhHn22WeNMcasW7fOfPfdd8YYY7744gvT0dERwBQAXEHm\nALCFvEGhOIMFazo6OlRXV6e6ujotWrRIGzZskCQlEgl98803On/+vD7//HM98MADuedMTU0FtVwA\njiNzANhC3uBKFCxYU1tbm/t3OBzOfR0OhzU9Pa3ff/9dixYt0uDgYFBLBOARMgeALeQNrsRdBFFS\ndXV1OnfuXFHPXbhwoVasWKF3331XkmSM0ddff13K5QHwDJkDwBbyBoWiYKGkrrvuOt1+++1qa2vT\njh07Zv38/fv3a9++fWpvb1dra6symcw8rBKAL8gcALaQNygUt2kHAAAAgBLhDBYAAAAAlAgFCwAA\nAABKhIIFAAAAACVCwQIAAACAEqFgAQAAAECJULAAAAAAoEQoWAAAAABQIhQsAAAAACiR/wHibr2d\nx2QGIQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f6656d0ed10>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "N = 200\n",
      "\n",
      "alpha_y = np.zeros(N)\n",
      "alpha_y[0:(N/2)] =0.9\n",
      "alpha_y[(N/2):-1] =0.85\n",
      "alpha_h = 0.5\n",
      "\n",
      "C   = np.zeros(N) # consumption\n",
      "Y   = np.zeros(N) # income\n",
      "H_h = np.zeros(N) # private wealth\n",
      "\n",
      "C[0] = 0\n",
      "Y[0] = 0\n",
      "H_h[0] = 100\n",
      "\n",
      "for t in range(1, N):\n",
      "    \n",
      "    # calculate consumer spending\n",
      "    C[t] = alpha_y[t]*Y[t-1] + alpha_h*H_h[t-1] \n",
      "        \n",
      "    # calculate total income (consumer spending plus constant government spending)\n",
      "    Y[t]   = alpha_h/(1-alpha_y[t])*H_h[t-1] \n",
      "    \n",
      "    H_h[t] = H_h[t-1] +(1-alpha_y[t])*Y[t] - alpha_h*H_h[t-1] "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "alpha_y = np.zeros(N)\n",
      "alpha_y[0:(N/2)] =0.9\n",
      "alpha_y[(N/2):-1] =0.8 \n",
      "alpha_y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "array([ 0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,  0.9,\n",
        "        0.9,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,  0.8,\n",
        "        0.8,  0. ])"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "N = 200\n",
      "\n",
      "alpha_y = 0.9\n",
      "alpha_h = 0.5\n",
      "\n",
      "C   = np.zeros(N) # nominal consumption\n",
      "c   = np.zeros(N) # real consumption\n",
      "Y   = np.zeros(N) # nominal income\n",
      "H_h = np.zeros(N) # private wealth\n",
      "\n",
      "sigma_T = 0.5\n",
      "chi     = 0.5\n",
      "beta    = 0.5\n",
      "\n",
      "y     = np.zeros(N) # real production\n",
      "s_e   = np.zeros(N) # expected sales\n",
      "s     = np.zeros(N) # realised real sales\n",
      "inv   = np.zeros(N) # real value of inventories\n",
      "inv_e = np.zeros(N) # expected inventories\n",
      "inv_T = np.zeros(N) # targeted inventories\n",
      "\n",
      "pr = 1           # productivity\n",
      "E = np.zeros(N)  # employment\n",
      "WB = np.zeros(N) # wage bill\n",
      "UC = np.zeros(N) # unit cost\n",
      "\n",
      "INV   = np.zeros(N) # nominal value of inventories\n",
      "S     = np.zeros(N) # realised nominal sales\n",
      "p     = np.zeros(N) # price level\n",
      "NHUC  = np.zeros(N) # normal historical unit cost\n",
      "\n",
      "C[0] = 0\n",
      "Y[0] = 0\n",
      "H_h[0] = 100\n",
      "\n",
      "for t in range(1, N):\n",
      "    \n",
      "    # calculate consumer spending\n",
      "    C[t] = alpha_y*Y[t-1] + alpha_h*H_h[t-1] \n",
      "        \n",
      "    # calculate total income (consumer spending plus constant government spending)\n",
      "    Y[t]   = alpha_h/(1-alpha_y)*H_h[t-1] \n",
      "    \n",
      "    H_h[t] = H_h[t-1] +(1-alpha_y)*Y[t] - alpha_h*H_h[t-1] "
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}