diff --git a/notebooks/Uncertainty-and-the-Saving-Rate.ipynb b/notebooks/Uncertainty-and-the-Saving-Rate.ipynb deleted file mode 100644 index 48589728..00000000 --- a/notebooks/Uncertainty-and-the-Saving-Rate.ipynb +++ /dev/null @@ -1,758 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Uncertainty and Saving in Partial Equilibrium\n", - "\n", - "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/econ-ark/DemARK/master?filepath=notebooks%2FUncertainty-and-the-Saving-Rate.ipynb)\n", - "\n", - "Saving rates vary widely across countries, but there is no consensus about the main causes of those differences.\n", - "\n", - "One commonly mentioned factor is differences across countries in the degree of uncertainty that individuals face, which should induce different amounts of precautionary saving.\n", - "\n", - "Uncertainty might differ for \"fundamental\" reasons, having to do with, say, the volatility of demand for the goods and services supplied by the country, or might differ as a result of economic policies, such as the strucutre of the social insurance system.\n", - "\n", - "A challenge in evaluating the importance of precautionary motives for cross-country saving differences has been a lack of consensus about what measures of uncertainty ought, in principle, to be the right ones to look at in any attempt to measure a relationship between uncertainty and saving.\n", - "\n", - "This notebook uses [a standard model](https://econ.jhu.edu/people/ccarroll/papers/cstwMPC) to construct a theoretical benchmark for the relationship of saving to two kinds of uncertainty: Permanent shocks and transitory shocks to income. \n", - "\n", - "Conclusions:\n", - "1. The model implies a close to linear relationship between the variance of either kind of shock (transitory or permanent) and the saving rate\n", - "2. The _slope_ of that relationship is much steeper for permanent than for transitory shocks\n", - " * Over ranges of values calibrated to be representative of microeconomically plausible magnitudes\n", - "\n", - "Thus, the quantitative theory of precautionary saving says that the principal determinant of precautionary saving should be the magnitude of permanent (or highly persistent) shocks to income.\n", - "\n", - "(Because the result was obtained in a partial equilibrium model, the conclusion applies also to attempts to measure the magnitude of precautionary saving across groups of people who face different degrees of uncertainty within a country).\n", - "\n", - "@authors: Derin Aksit, Tongli Zhang, Christopher Carroll" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "code_folding": [ - 0, - 11 - ] - }, - "outputs": [], - "source": [ - "# Boring non-HARK setup stuff\n", - "\n", - "Generator = True # This notebook can be used as a source for generating derivative notebooks\n", - "nb_name = 'Uncertainty-and-the-Saving-Rate'\n", - "\n", - "# This is a jupytext paired notebook that autogenerates BufferStockTheory.py\n", - "# which can be executed from a terminal command line via \"ipython BufferStockTheory.py\"\n", - "# But a terminal does not permit inline figures, so we need to test jupyter vs terminal\n", - "# Google \"how can I check if code is executed in the ipython notebook\"\n", - "\n", - "from IPython import get_ipython # In case it was run from python instead of ipython\n", - "def in_ipynb():\n", - " try:\n", - " if str(type(get_ipython())) == \"\":\n", - " return True\n", - " else:\n", - " return False\n", - " except NameError:\n", - " return False\n", - "\n", - "# Determine whether to make the figures inline (for spyder or jupyter)\n", - "# vs whatever is the automatic setting that will apply if run from the terminal\n", - "if in_ipynb():\n", - " # %matplotlib inline generates a syntax error when run from the shell\n", - " # so do this instead\n", - " get_ipython().run_line_magic('matplotlib', 'inline')\n", - "else:\n", - " get_ipython().run_line_magic('matplotlib', 'auto')\n", - " print('You appear to be running from a terminal')\n", - " print('By default, figures will appear one by one')\n", - " print('Close the visible figure in order to see the next one')\n", - "\n", - "# Import the plot-figure library matplotlib\n", - "\n", - "import matplotlib.pyplot as plt\n", - "\n", - "# In order to use LaTeX to manage all text layout in our figures, we import rc settings from matplotlib.\n", - "from matplotlib import rc\n", - "plt.rc('font', family='serif')\n", - "\n", - "# LaTeX is huge and takes forever to install on mybinder\n", - "# so if it is not installed then do not use it \n", - "from distutils.spawn import find_executable\n", - "iflatexExists=False\n", - "if find_executable('latex'):\n", - " iflatexExists=True\n", - " \n", - "plt.rc('font', family='serif')\n", - "plt.rc('text', usetex=iflatexExists)\n", - "\n", - "# The warnings package allows us to ignore some harmless but alarming warning messages\n", - "import warnings\n", - "warnings.filterwarnings(\"ignore\")\n", - "\n", - "import os\n", - "from copy import copy, deepcopy\n", - "\n", - "# Define (and create, if necessary) the figures directory \"Figures\"\n", - "if Generator:\n", - " nb_file_path = os.path.dirname(os.path.abspath(nb_name+\".ipynb\")) # Find pathname to this file:\n", - " FigDir = os.path.join(nb_file_path,\"Figures/\") # LaTeX document assumes figures will be here\n", - "# FigDir = os.path.join(nb_file_path,\"/tmp/Figures/\") # Uncomment to make figures outside of git path\n", - " if not os.path.exists(FigDir):\n", - " os.makedirs(FigDir)\n", - "\n", - "from copy import deepcopy\n", - "from scipy.optimize import golden, brentq\n", - "from time import time\n", - "import numpy as np\n", - "import scipy as sp" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Import HARK tools and cstwMPC parameter values\n", - "from HARK.utilities import plotFuncsDer, plotFuncs\n", - "from HARK.ConsumptionSaving.ConsIndShockModel import PerfForesightConsumerType\n", - "import HARK.cstwMPC.cstwMPC as cstwMPC\n", - "import HARK.cstwMPC.SetupParamsCSTW as Params\n", - "\n", - "# Double the default value of variance\n", - "# Params.init_infinite['PermShkStd'] = [i*2 for i in Params.init_infinite['PermShkStd']]" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": { - "code_folding": [ - 0 - ], - "scrolled": false - }, - "outputs": [], - "source": [ - "# Setup stuff for general equilibrium version\n", - "\n", - "# Set targets for K/Y and the Lorenz curve\n", - "lorenz_target = cstwMPC.getLorenzShares(Params.SCF_wealth,weights=\n", - " Params.SCF_weights,percentiles=\n", - " Params.percentiles_to_match)\n", - "\n", - "lorenz_long_data = np.hstack((np.array(0.0),\\\n", - " cstwMPC.getLorenzShares(Params.SCF_wealth,weights=\\\n", - " Params.SCF_weights,percentiles=\\\n", - " np.arange(0.01,1.0,0.01).tolist()),np.array(1.0)))\n", - "KY_target = 10.26" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Setup and calibration of the agent types\n", - "\n", - "# The parameter values below are taken from\n", - "# http://econ.jhu.edu/people/ccarroll/papers/cjSOE/#calibration\n", - "\n", - "Params.init_cjSOE = Params.init_infinite # Get default values of all parameters\n", - "# Now change some of the parameters for the individual's problem to those of cjSOE\n", - "Params.init_cjSOE['CRRA'] = 2\n", - "Params.init_cjSOE['Rfree'] = 1.04**0.25\n", - "Params.init_cjSOE['PermGroFac'] = [1.01**0.25] # Indiviual-specific income growth (from experience, e.g.)\n", - "Params.init_cjSOE['PermGroFacAgg'] = 1.04**0.25 # Aggregate productivity growth \n", - "Params.init_cjSOE['LivPrb'] = [0.95**0.25] # Matches a short working life \n", - "\n", - "PopGroFac_cjSOE = [1.01**0.25] # Irrelevant to the individual's choice; attach later to \"market\" economy object\n", - "\n", - "# Instantiate the baseline agent type with the parameters defined above\n", - "BaselineType = cstwMPC.cstwMPCagent(**Params.init_cjSOE)\n", - "BaselineType.AgeDstn = np.array(1.0) # Fix the age distribution of agents\n", - "\n", - "# Make desired number of agent types (to capture ex-ante heterogeneity)\n", - "EstimationAgentList = []\n", - "for n in range(Params.pref_type_count):\n", - " EstimationAgentList.append(deepcopy(BaselineType))\n", - " EstimationAgentList[n].seed = n # Give every instance a different seed" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Make an economy for the consumers to live in\n", - "\n", - "EstimationEconomy = cstwMPC.cstwMPCmarket(**Params.init_market)\n", - "EstimationEconomy.print_parallel_error_once = True # Avoids a bug in the code\n", - "\n", - "EstimationEconomy.agents = EstimationAgentList\n", - "EstimationEconomy.act_T = Params.T_sim_PY # How many periods of history are good enough for \"steady state\"" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Uninteresting parameters that also need to be set \n", - "EstimationEconomy.KYratioTarget = KY_target\n", - "EstimationEconomy.LorenzTarget = lorenz_target\n", - "EstimationEconomy.LorenzData = lorenz_long_data\n", - "EstimationEconomy.PopGroFac = PopGroFac_cjSOE # Population growth characterizes the entire economy\n", - "EstimationEconomy.ignore_periods = Params.ignore_periods_PY # Presample periods\n", - "\n", - "#Display statistics about the estimated model (or not)\n", - "EstimationEconomy.LorenzBool = False\n", - "EstimationEconomy.ManyStatsBool = False\n", - "EstimationEconomy.TypeWeight = [1.0]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "code_folding": [ - 0 - ], - "lines_to_next_cell": 2 - }, - "outputs": [], - "source": [ - "# construct spread_estimate and center_estimate if true, otherwise use the default values\n", - "Params.do_param_dist=True # Whether to use a distribution of ex-ante heterogeneity\n", - "\n", - "# Discount factors assumed to be uniformly distributed around center_pre for spread_pre on either side\n", - "\n", - "spread_pre=0.0019501105739768 #result under the default calibration of cjSOE\n", - "center_pre=1.0065863855906343 #result under the default calibration of cjSOE\n", - "\n", - "do_optimizing=False # Set to True to reestimate the distribution of time preference rates\n", - "\n", - "if do_optimizing: # If you want to rerun the cstwMPC estimation, change do_optimizing to True\n", - " # Finite value requires discount factor from combined pure and mortality-induced\n", - " # discounting to be less than one, so maximum DiscFac is 1/LivPrb\n", - " DiscFacMax = 1/Params.init_cjSOE['LivPrb'][0] # \n", - " param_range = [0.995,-0.0001+DiscFacMax] \n", - " spread_range = [0.00195,0.0205] # \n", - "\n", - " if Params.do_param_dist: # If configured to estimate the distribution\n", - " LorenzBool = True\n", - " # Run the param-dist estimation\n", - " paramDistObjective = lambda spread : cstwMPC.findLorenzDistanceAtTargetKY(\n", - " Economy = EstimationEconomy,\n", - " param_name = Params.param_name,\n", - " param_count = Params.pref_type_count,\n", - " center_range = param_range,\n", - " spread = spread,\n", - " dist_type = Params.dist_type) # Distribution of DiscFac\n", - " t_start = time()\n", - " \n", - " spread_estimate = golden(paramDistObjective \n", - " ,brack=spread_range\n", - " ,tol=1e-4) \n", - " center_estimate = EstimationEconomy.center_save\n", - " t_end = time()\n", - " else: # Run the param-point estimation only\n", - " paramPointObjective = lambda center : cstwMPC.getKYratioDifference(Economy = EstimationEconomy,\n", - " param_name = Params.param_name,\n", - " param_count = Params.pref_type_count,\n", - " center = center,\n", - " spread = 0.0,\n", - " dist_type = Params.dist_type)\n", - " t_start = time()\n", - " center_estimate = brentq(paramPointObjective # Find best point estimate \n", - " ,param_range[0]\n", - " ,param_range[1],xtol=1e-6)\n", - " spread_estimate = 0.0\n", - " t_end = time()\n", - " \n", - " print(spread_estimate)\n", - " print('****************')\n", - " print(center_estimate)\n", - " print('****************')\n", - "else: # Just use the hard-wired numbers from cstwMPC\n", - " center_estimate=center_pre\n", - " spread_estimate=spread_pre" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Construct the economy at date 0\n", - "EstimationEconomy.distributeParams( # Construct consumer types whose heterogeneity is in the given parameter\n", - " 'DiscFac',\n", - " Params.pref_type_count,# How many different types of consumer are there \n", - " center_estimate, # Increase patience slightly vs cstwMPC so that maximum saving rate is higher\n", - " spread_estimate, # How much difference is there across consumers\n", - " Params.dist_type) # Default is for a uniform distribution" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Function to calculate the saving rate of a cstw economy\n", - "def calcSavRte(Economy,ParamToChange,NewVals):\n", - " '''\n", - " Calculates the saving rate as income minus consumption divided by income.\n", - " \n", - " Parameters\n", - " ----------\n", - " Economy : [cstwMPCmarket] \n", - " A fully-parameterized instance of a cstwMPCmarket economy\n", - " ParamToChange : string\n", - " Name of the parameter that should be varied from the original value in Economy\n", - " NewVals : [float] or [list]\n", - " The alternative value (or list of values) that the parameter should take\n", - "\n", - " Returns\n", - " -------\n", - " savRte : [float]\n", - " The aggregate saving rate in the last year of the generated history\n", - " '''\n", - " for NewVal in NewVals:\n", - " if ParamToChange in [\"PermShkStd\",\"TranShkStd\"]:\n", - " ThisVal = [NewVal]\n", - " else:\n", - " ThisVal = NewVal # If they asked to change something else, assume it's a scalar\n", - " \n", - " for j in range(len(Economy.agents)): # For each agent, set the new parameter value\n", - " setattr(Economy.agents[j],ParamToChange,ThisVal)\n", - " cstwMPC.cstwMPCagent.updateIncomeProcess(Economy.agents[j]) \n", - " \n", - " Economy.solve()\n", - " \n", - " C_NrmNow=[]\n", - " A_NrmNow=[]\n", - " M_NrmNow=[]\n", - " for j in range (len(Economy.agents)): # Combine the results across all the agents\n", - " C_NrmNow=np.hstack((C_NrmNow,Economy.agents[j].cNrmNow))\n", - " A_NrmNow=np.hstack((A_NrmNow,Economy.agents[j].aNrmNow))\n", - " M_NrmNow=np.hstack((M_NrmNow,Economy.agents[j].mNrmNow))\n", - " CAgg=np.sum(np.hstack(Economy.pLvlNow)*C_NrmNow) # cNrm times pLvl = level of c; sum these for CAgg\n", - " AAgg=np.sum(np.hstack(Economy.pLvlNow)*A_NrmNow) # Aggregate Assets\n", - " MAgg=np.sum(np.hstack(Economy.pLvlNow)*M_NrmNow) # Aggregate Market Resources\n", - " YAgg=np.sum(np.hstack(Economy.pLvlNow)*np.hstack(Economy.TranShkNow)) # Aggregate Labor Income\n", - " BAgg=MAgg-YAgg # Aggregate \"Bank Balances\" (at beginning of period; before consumption decision)\n", - " IncAgg=(BaselineType.Rfree-1)*BAgg+YAgg # Interest income plus noninterest income\n", - " savRte=(IncAgg-CAgg)/IncAgg # Unspent income divided by the level of income\n", - " return savRte" - ] - }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Function to plot relationship between x and y; x is the parameter varied and y is saving rate\n", - "def plotReg(x,y,xMin,xMax,yMin,yMax,xLbl,yLbl,Title,fileName):\n", - " # Result_data_path = os.path.join(Folder_path,'SavingVSPermShr_Youth_MPC_15.png')\n", - " plt.ylabel(yLbl)\n", - " plt.xlabel(xLbl)\n", - " plt.title(Title)\n", - " plt.xlim(xMin,xMax)\n", - " plt.ylim(yMin,yMax)\n", - " plt.scatter(x,y)\n", - " # Draw the linear fitted line\n", - " m, b = np.polyfit(x, y, 1)\n", - "# plt.plot(x, m*np.asarray(x) + b, '-')\n", - " if Generator:\n", - " plt.savefig(FigDir + nb_name + '-' + fileName + '.png')\n", - " plt.savefig(FigDir + nb_name + '-' + fileName + '.svg')\n", - " plt.savefig(FigDir + nb_name + '-' + fileName + '.pdf')\n", - " slope, intercept, r_value, p_value, std_err = sp.stats.linregress(x,y)\n", - " print('Slope=' + str(slope) + ', intercept=' + str(intercept) + ', r_value=' + str(r_value) + ', p_value=' + str(p_value)+', std=' + str(std_err))" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [], - "source": [ - "# Proportion of base value for uncertainty parameter to take (up to 1 = 100 percent)\n", - "# Do not go above one to avoid having to worry about whether the most patient consumer violates the \n", - "# Growth Impatience Condition (https://econ.jhu.edu/people/ccarroll/papers/BufferStockTheory/#GIC)\n", - "bottom=0.5\n", - "points=np.arange(bottom,1.+0.025,0.025)" - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": { - "code_folding": [ - 0 - ], - "scrolled": true - }, - "outputs": [], - "source": [ - "# Calculate variance of permanent shock vs saving measures\n", - "savRteList = []\n", - "KtoYList = []\n", - "pVarList = []\n", - "pVarBase = BaselineType.PermShkStd[0] ** 2\n", - "for pVar in points * pVarBase:\n", - " pVarList.append(pVar) # Variance is square of standard deviation\n", - " pStd = pVar ** 0.5\n", - "# print(pStd)\n", - " savRteList.append(calcSavRte(EstimationEconomy,\"PermShkStd\",[pStd]))\n", - " KtoYList.append(0.25*np.mean(np.array(EstimationEconomy.KtoYnow_hist)[EstimationEconomy.ignore_periods:]))" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Halving the magnitude of the permanent variance causes target wealth to fall to 0.646\n", - "of its original value.\n" - ] - } - ], - "source": [ - "# Calculate how much net worth shrinks when permanent variance is halved \n", - "ShrinksBy = KtoYList[1]/KtoYList[-1]\n", - "print('Halving the magnitude of the permanent variance causes target wealth to fall to %1.3f' % ShrinksBy)\n", - "print('of its original value.')" - ] - }, - { - "cell_type": "code", - "execution_count": 14, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Slope=43.87569286985961, intercept=0.07612904106465898, r_value=0.9957594503438054, p_value=3.715857043912896e-21, std=0.9299464143538145\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEYCAYAAABY7FHWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHIpJREFUeJzt3c+SG9d1x/Hfie0FK1USRIpZmFUKBVZKK1XEIfQCEmhVdip6KOYBzGGylzhUXoAeKpusbI78AtSQjDdZ2Bz6AaLh0FXexEkGihfyhuYQ0kYLl+tk0adnusEG0PjT0w3g+6maItDoBg5aUJ/ue++5be4uAABSf1V3AACAZiExAABySAwAgBwSAwAgh8QAAMghMQAAckgMAIAcEgMAIIfEAADI+X7dAWC1mFlb0paktqTr7r5vZi1Jn0s6lLTp7v06Y0yZ2YG7X6g7jiwz24iHh/HvaXffnuJ9Gvfd0BxcMeBEuXtP0iNJe+6+H8v6km5L2qo6KWQOrGONO3BO8l7zYGZrklruvu3u9yXtS5rq4E5SwCgkBqyMuDK5UXLdtpmtz+O95ui0pHfTJ5Fk7036JuO+G0BiQOOYWdfMDuLfrpndjQNx+vpGLF+Ls2iZ2c1YthEHvux77MT2bUktM1tPt8t83pqZbUVTl5Q01WyZWasoHknvZN8rlr/IxHPXzG4O+W4vrRfvkb7XVtF+cffd2OYgtuumV11F32NETCO/28C+3or3TPc5CWUF0MeAxnH3XTPrSeq5ey8ObB1Ju3FguuDu23EQ3zSzJ+l2kmRmO+5+Nd7j0N2vxlvvm1k/mmGybsT6krQZz/ux/bB4vi8p915mtq3krF6SdtJ4Cr7bS+tFMvgy+lxG7Zur8b27ku6a2VamjyH3Pdz9xrCYxny3dF+3lTRd7ZtZNx4P7jssIa4YUIeejg9WqdM67lBNZZ+nZ7HvSvpSSppS3P2GpEuS+pkriC/T98yeUY+wGQmnM2a9oniybsd7tZR8x2Gy6x3Gd9mU1I4kV9hElZ75x/fejn6C7LpF3+Olzyr73dL+oPjcfZLC6iAxoA57Spp1slolO56/VKadPQ54jyTJ3fcjEYwapXMY23Uz/27GQS89mx6MrdR7Rfw9SRtxUC0U6/Vjvf14jw13v+PuaZIriqGdflZs01KyL4d+j6LPmlTs15eufrC8aErCiYtmmuvR3p0eQI8OPHGG2pb0kZl9IemypJ6Z7br7/UznaV9JU9H9aAM/WhbNKe044GYTxY4lo4nSzz3MfKaUXLm0M30SH5nZngrikfQfA+8lSXf1ctIrclv5M/tWJv6DEYklXU/xOZujvkfElvusgf1b+N1iX/cl3TCza/H+T6YZGovFY9zBDZhdnJ33okN4Kc6uI+ntxvdqKenXmKpuAouFKwZgPtbNbFej+xYWzZ6ktewoJWWu7LC8uGIAAOTQ+QwAyKkkMUShTreowCde34i/lwp5hm0DADgZc+9jyIy13o3RI2uD1Zk67tDayXbWxWuXJd0Z9Rmvv/66nz9/ft6hA8BSe/LkyZ/c/ey49arofL6mGFeupCOuq2Syr1Q7/rbj9bJjxo+cP39ee3t7M4YJAKvFzP5QZr0qmpIGKyzPZF+Mis10uNuajgt01pZlmB8ALLLaOp8zZfbp1cTgFAkAgBpUkRj6Oj7ItyQ9H7JeN+aHKXW1EJ3Ve2a29+zZs/lFCwDIqSIx3NNxv0Fbx/O25KZNdvc78birZAqC9ai0PJ0p6z8STVAdd++cPTu27wQAMKW5J4bMpGBdJdMSp01FjzPLt2IO+Bexzf3MzI1Fs1YCAE7IQlY+dzodZ1QSAEzGzJ64+7jp5al8BgDkkRgAADkkBgBADokBAJBDYgAA5JAYAAA5JAYAQA6JAQCQQ2IAAOSQGAAAOSQGAEAOiQEAkFPFrT0BABm/fPq1PvvV7/XH/nf6YeuUPvngLX148VzdYQ1FYgCACv3y6df69OHv9N2f/yJJ+rr/nT59+DtJamxyoCkJACr02a9+f5QUUt/9+S/67Fe/rymi8UgMAFChP/a/m2h5E5AYAKBCP2ydmmh5E5AYAKBCn3zwlk794Hu5Zad+8D198sFbNUU0Hp3PAFChtIOZUUkAsISmHXb64cVzjU4Eg0gMAFDCIg47nRZ9DABQwiIOO50WiQEASljEYafTIjEAQAmLOOx0WiQGAChhEYedTovOZwArZ5rRRYs47HRaJAYAK2WW0UWLNux0WjQlAVgpqzS6aFokBgArZZVGF02LxABgpazS6KJpkRgArJRVGl00LTqfAayUVRpdNC0SA4CFNMt9lFdldNG0SAwAFs4qTWhXB/oYACwchpxWi8QAYOEw5LRaJAYAC4chp9WqJDGY2bqZdc3s5pDXN+Jva9QyACjCkNNqzT0xmNmaJLn7rqR++jzzelfSrrtvS2pHAnlp2bzjArA8Prx4TrevvK1zrVMySedap3T7ytt0PM9JFaOSrkl6FI97krqS9jOvt+NvO15vZ5YPLgOAQgw5rU4ViaEl6TDz/Ez2xbgqSK1Juufu+4PLKogLQAPNUo+AatRWxxBNTPvZpFC0DMDyoh6hmarofO5LOh2PW5KeD1mv6+6bJZZJOuqc3jOzvWfPns0pVAB1oh6hmapIDPeU7zfYlSQza6UrmNmGu9+Jx91hy7LcfdvdO+7eOXv2bAVhAzhp1CM009wTQ9oMFAf3fqZZ6HFm+ZaZHZjZi2HLACw/6hGaqZI+hoEO5nTZpfh3V9JrBZsVLQOwxD754K1cH4NEPUITMIkegNowBXYzkRgAzMW0w06pR2geEgOAmTHsdLkwiR6AmTHsdLmQGADMjGGny4XEAGBmDDtdLqUSg5ldN7OfmtnHZvaKmb1XdWAAFgfTYC+Xsp3PB+7+uZlddPdvzazSoAAsFoadLpeyieFSJIOWmbmkS5J+U1lUABYOw06XR9nEsC3pUyVzH/2nu39WXUgA6sIU2JBKJgZ3/0bSLUkys4tm9oq7f1tpZABOFLUISJXtfL6SPnb3p0ruygZgiVCLgNTIKwYz+7Gky5I6ZnYjFveV3H7zYcWxAThB1CIgNTIxuPsDM9uV1HH3xycUE4Aa/LB1Sl8XJAFqEVbP2KYkd/8mmxTM7Hy2aQnAcqAWAalSnc9mdl3SDSW36TRJT0RTErBUqEVAqvTsqu7eMbP33f2xmb1fZVAA6kEtAqTyieHQzD6W1DOzn0hqKW7VCaB5qEfALEoNV3X3B5Ieu/tDSX8t6b8rjQrA1NJ6hK/738l1XI/wy6df1x0aFkTp2VWjfkHu/m+S/qayiADMhHoEzGpoYogK50Mze25mf29mb5rZz83sV0pqGwA0EPUImNWoPob33f20mbUl/VzSF5IeSepT0wA0F/UImNWopqRvJMnde5K23P0X7v6ApAA0G/UImNWoK4a2mb0Tj9/MPJaka+7+aYVxAZgS9QiYlbl78Qtm/6ukkK3orjxvuvu7VQY2SqfT8b29vbo+HgAWkpk9cffOuPVGXTHcGNZsZGZvTh0ZgNKoR0AdhiaGUX0J7v5VNeEASHF/BNSldB0DgJNFPQLqUvZGPT8ZeP6qmf3UzN6rJiwA1COgLmWvGF4zs3tmdj6e35J0V9JrVQQFYHjdAfUIqFrZxHDg7tckteN5O/oZ+tWEBYB6BNSl7OyqF9LmJDP7Kp6/IunVyiIDVhz1CKhL2cSwLakbt/q8GPdm+ERJnQOAinB/BNShVGJw928kPYjHT83sFXf/rNLIAAC1KHtrz4uSrklyJZXQFyV9UGFcwNKgSA2LpmxTUlfJKKTUegWxAEuHIjUsorKjkp64+1fpn5LptwGMQZEaFlHZK4ZbZrYl6VBJU9Kbkv5u2Mpmtq5kKOuau98peH0jHl5w980y2wCLiCI1LKKyiWErO3eSmb0/bEUzW5Mkd981s7aZrbn7fub1rqRdd++Z2U48Pxy1DbCouGkOFlGppqSCCfUORqx+TceFbz0l/RNZ7cyyXjwftw2wkChSwyIaesVgZvckXZd0QdKWpBfpS0pGJQ1rSmoprgDCmeyL7r6debom6Z6kS6O2ARYVRWpYRKOakm65+7dm1ldyb4ajqbZj+OpMoslp3933zYruBQQsB4rUsGhG3Y/hq/TfuHq4lnnt6Yj37Es6HY9bkp4PWa+bdjyX2SY6rDck6Y033hjx8UA1qEfAqig7XPULM7toZu8N3Pu5yD1lJtuTtCtJZtZKVzCzjXTkUXQ+F26T5e7b7t5x987Zs2dLhg3MR1qP8HX/O7mO6xF++fTrukMD5q5sYngUVwmvSfoXM/vZsBXT0URxwO9nRhc9zizfMrMDM3sxZhugEahHwCopO1x138wOJO1Iuh5zJw010MGcLrsU/+6q4D4ORdsATUE9AlZJ2cSw6e4PKo0EaDDqEbBKytYxHCUFMztvZleqCwloHuoRsErKzq56XdINJaOFTMl9GB5WGBfQKNQjYJWUbUpS3JznfXd/PGpKDGBZUY+AVVF2VNKhmX0s6dW4xefMBW4AgGaapI/hsbs/VNKUNKrADQCwwCZpSnoa/35eXThA9ahgBkYbesVgZp+Y2c/SSud4/Gszu1ei+hloJCqYgfFGNSX13P2f3f23ZvaJpNPu/iN3vyapc0LxAXNFBTMw3qjE4JnHHyl/z+cXAhYQFczAeKP6GC6Y2XuSfiTJ3P03UlLgpnzSABYGFczAeEOvGNz9MyUjkL50944kmdmbki6fUGzA3FHBDIw3clTS4C094x4NjErCwqKCGRiv9HBVYFlQwQyMVrbyGQCwIrhiwMKiUA2oxiSzq16Q9CdJ25I66SgloA5poVpak5AWqkkiOQAzKtuUdODut5TMl/RtlQEBZVCoBlSnbFPSJTOTpJaZuaRLkrhiQG0oVAOqU/aKYVtJods/SupGjQNQm2EFaRSqAbMrO+32N+5+y90/kvTYzF6pOC5gJArVgOqUSgzZezzH9NvdyiICSvjw4jndvvK2zrVOySSda53S7Stv0/EMzMHIPgYz+7GSKTA6ZnYjFvcl9cQ9n1EzCtWAaoybEuOBme0qGZ76eNS6AIDlMLYpKfoXjpKCmZ3PNi0BAJbLJAVuNyQ9VzLj6hPRlIQ5oHoZaJ5J7vncMbP33f2xmb1fZVBYDVQvA81Uto7h0Mw+lvSqmf1E0sUKY8KKoHoZaKaydQwPlEyH8VBJU9LTSqPCSqB6GWimsnUMP4n6Bbn755L2zOyncetPYCpULwPNVLYp6TUzuxf3e5akW5LuSnqtiqCwGqheBpppktlVr0lqx/N23OazX01YWAVULwPNVHZU0oXodJaZfRXPX5H0amWRYSVQvQw0zySzqx66+y8ktdy9o6SugSsGAFgyZa8YPlJyldCW9LmZvcfU2wCwnMomhgN3/9zMLrr7N3HTHuAIFczA8qjkDm5mtq6kmWnN3e8MWWfN3fcLtmm7+3bJuNAAVDADy2Xud3AzszVJcvddSf30+cA6XUk7A9v0Ypte0TZoLiqYgeVSNjFcl3Tb3T9y938ds+41HXdK91RwU580AQws3op/29krCTQfFczAcimbGHru/k36xMzeGbFuS9Jh5vmZcW8eiaBnZi8GtsUCoIIZWC5lE8M/mdn/RPXzF8o0A82DmbWUXGXcVjLqqT1mEzQIFczAcinb+Xw3JtKTJI2Zdrsv6XQ8bim5h8M4G0qaqvpm1pO0LinXaW1mG7Ge3njjjZJh4ySkHcyMSgKWQ9nEMDgn0p6Z/VTSr919cHTSPUmdeNyWtCslVwXuPrYgzt3vRxIYXL6tpBNcnU7HS8aNE0IFM7A85j6JXtpxHCOP+pmO5OztQdcldeJfxZDWDTNbN7MNhqsCQH0mKXD7LKbZ/j/FJHrD+gKKDuzufinz+L6k+wOvF9Y7AABO1qyT6LUqiwwAUIuyiWFbSWHbg5gWoxO3+qTeYMkwtQWAsonhasysKnd/GlcLr1cXFurA1BYApOk7nz9V0vlMU9ISYWoLANLsd3D7ZsQ2WDBMbQFAKp8Y0s7ntpm9Ke7gtpSY2gKANNnsqi+4g9tyY2oLAFLJzueYQO9BPH5qZh9zB7flw9QWAKTyo5IUxW1XJV1WUvE8bvptLCCmtgAwNDFEH0JXyf0VLkvak3Sg5O5tzH4KAEtq1BXDvyiZ5XTH3U9LkpldiWalpycRHADg5A1NDO5+S9ItM7sYI5JMyZXCQzN7x91/e1JBYjJULwOYxdg+Bnd/qrhCMLM3zezHSmZXfbfi2DAFqpcBzKrscFVJkrt/FTfsuVVRPJgR1csAZjVRYki5++Pxa6EOVC8DmNVUiQHNRfUygFmRGJYM1csAZlW6wA2LgeplALMiMSwhqpcBzIKmJABADokBAJBDYgAA5JAYAAA5dD43GHMeAagDiaGhmPMIQF1oSmoo5jwCUBcSQ0Mx5xGAupAYGoo5jwDUhcTQUMx5BKAudD43FHMeAagLiaHBmPMIQB1oSgIA5JAYAAA5JAYAQA59DCeAqS0ALBISQ8WY2gLAoqEpqWJMbQFg0VRyxWBm65L6ktbc/c6QddbcfT/7XFJbktz9fhVx1YGpLQAsmrlfMcQBXu6+K6mfPh9YpytpZ2Dxp5EQ2kXbLCqmtgCwaKpoSrqm5GpBknqSuoMrRNLopc/jCuPLeO1O9kpi0TG1BYBFU0ViaEk6zDw/U2KbdyWdMbM1M7tZQUy1+fDiOd2+8rbOtU7JJJ1rndLtK2/T8QygsZo0Kum5u++bWdfM1pepn4GpLQAskiquGPqSTsfjlqTnJbZ5ruOmpb6SK4gcM9swsz0z23v27NlcAgUAvKyKxHBPMboo/t2VJDNrjdjmfmablqK/Icvdt9294+6ds2fPzjFcAEDW3BND2nEcI4/6mY7kx+k60dnciX/l7j0lI5jWJZ1ZpmYkAFg05u51xzCxTqfje3t7dYcBAAvFzJ64e2fcelQ+AwBymjQqqdGYCA/AqiAxlMBEeABWCU1JJTARHoBVQmIogYnwAKwSEkMJTIQHYJWQGEpgIjwAq4TO5xLSDmZGJQFYBSSGkpgID8CqoCkJAJBDYgAA5JAYAAA5JAYAQM7KdT4z5xEAjLZSiYE5jwBgvJVqSmLOIwAYb6USA3MeAcB4K5UYmPMIAMZbqcTAnEcAMN5KdT4z5xEAjLdSiUFiziMAGGelmpIAAOORGAAAOSQGAEAOiQEAkENiAADkkBgAADnm7nXHMDEzeybpDwUvvS7pTyccziSaHB+xTa/J8RHb9Joc37Sx/a27nx230kImhmHMbM/dO3XHMUyT4yO26TU5PmKbXpPjqzo2mpIAADkkBgBAzrIlhu26AxijyfER2/SaHB+xTa/J8VUaW+P6GMxsXVJf0pq73ynz+rBtzGzN3ffTx5KeSOrFy7vufsPMttx908w23L3JPwQAOBGNumKIg7fcfVdSP30+6vVh25hZV9JOZvPT7m7ufkHSVUlbsXzDzA50nDBGxbduZl0zu1n29SHLNuJva9Jta4yvaNlW+lrNsb0UxyT7rqrY4vfpZnYQf3eHxXtC8XXjb26/u4pja8pvrii2mX5zVcY3j99doxKDpGtKzvyl5EDdLfF64TaRKI4O9vE81XH39LXr7n5h4PWX2JySliUJazeuTtrxH7R0wqspvpeWxduWSqpVxlYUxyT7ruLY5nEyMs/4rsaytRHr1bXvBmNr0m8uF1tRHDX//zoY38y/u6Ylhpakw8zzMyVeH7dNTuzILzKL0oPfuCw/r6TVzmzbi+elE15N8RUtk0om1YpjK4pjkn1XWWyznozMMz5333X3G7GsHU2ss/7uqoytEb+5IbEVxVHL/69F8c3jd9e0xHASLrt7unPl7ndiR53JnJUUmUvScvftTF/GmqS9stuOiK3S+IbELJVPqlXuu6I4Jtl3Vcc2y8nI3OLLxHJTUnogmfV3V1lsTfnNFcU2JI5a/n8dEd9Mv7umJYa+pNPxuCXpeYnXx20z6OiSzZL2y/V4+lzHZyWVi0u+/cwZSKMUxTe4bIKkWmlsdcVRJrYw7cnI3MXgjBtm1jrJzy2jKLam/LcejK1Jv7k0Hr3833Xq313TEsM9HR+c25J2JSnzZYteL9ymiJkNHvj3MutfUOYsr8C8k1bX3Ten3Pak43tp2YRJtbLYhsQxyb47if02y8nIXOLLtlEraYbYmOK7nWRsqVp/c0WxzeE3V2l8mfeY+nfXqMTgx0NLu5L6mbOux8NeH7ZN7IROZmeksh3S+5I+inUOxpy9zy1pWTI0Nh1q251k25riK1o2SVKtMraiOCbZd1Xvt1lORuYZX1f5g0pvku9WQ2xN+c0VxTbrb67q+Gb+3TUqMUhHbYvpaIR02aUxrxctu+/ur7n7/cyyXqajJrvtfS+omRhYby5JKx5vWTKM7MUk29YV34iYSyXVE9h3uTgm2XdVxpYx7cnIPE+W0tE9G7Hd/Vl/d1XG1pTf3Ij9NvVvrur4Mh8z9e+ucQVuAIB6Ne6KAQBQLxIDACCHxAAAyPl+3QEAWA0xUqalZCTNfT+uyEXDcMWAUiypmHxh+UnDbprZozLFUmbWNrOdcetVwZKJx7LjuNOJxg7MbCte32pCodK04r/FevzdjGVdM3s0xXutV/Tfak3JSJldSYPDyNEgJAaUEhWTXyhfkr+vZAKvfvFWue17kq5XFN4412Io89HY8hiuty/pXgxB3FR+Nt7GKajJSZd3lVQH34/hihek48nWJv2cgSGPcxPx9RVXDFV8BuaDxIBJ7CiZyCvVKpMUpKNmhBObcqRIiVgPCwqDGiGuyi4PeXnwtbvVRzQdO541lWakBqOPAaW5++6wJoY4mz0d623H8xtKDlJtJWeIW4oDWMH6XUmbsc6akoNHWtBzU8nZfVvJVctG+twHbq4UTV3p7KbbMV1A28zWR50Jx4G37+69eHz0GUqukrLfZX8g1p6SM/OrkjbdvV/2+xV8Vm/IfugoKvkHv0cUhH1qZk+UXAHlijXjswf36eB+OtrH2X0a73lXA1M1+/hp6ruxPz5SUmV7N/bRpqSemT2q6soEc+Du/PFX+k/J/+DrSpoDWrFsTdLNePwks+5B/JuutzNm/Sfxb1vSVjzekLQej28qOWB24/nWQGw3M691M5+xM+S77KTbxOekcb70GQXfJRvrTibW9CBc9vsVfdZL68XzR2P+26Q3p7o75jMH99O/Z/dxZt90lSSKwTjbY+JopTHEb2W97t8tf5P90ZSESaXNSUfNSJ6chfbiLHGwD0I+0IQzYv3s49QlxdmqJ2fCbUmtuBIYnHjsXR2f2fbi+Tj7HtOpZOIs+ozB75KNNXs23Zrw+xV9VtF6Q6XNX/E9rirfZFf0XoP76R3l93Ea1w0lB3lJui3pclxBjBtssKHjm8Oc1hT9HKgXiQET8aQJITd6J5olTsdrRRN4aYb1DxQHumh2+VJSz4/nicnK3tClHetOY9RnjDXh95vks/rxfoN3B+sOfMa4A/Hgfvov5fdxus51HR/gu+6+6cm8Zd2BdQed8eM+hEte0OxU8B3QIPQxYBq5ET5KDiLpbQb343E7/l3zpC19LbN82PrtWK8by1rufieGkkpKzmhjaGY6o2R2pNFmvCZJa7Fu+rnd7AEq7XuQdNXM9nxg3vrsZ6Tvl/ku3YFY0wPzZSVn6DsTfr/Bzypar6/kKmRdL8/ceZjGp+QM/Xa812Cc6XsN7qd/yO5jM+sraQ5LY7qr5FaS6efdj6TwRDECasDdSI6HanBHOIZjEj0AUzGztk8wuigS1LVITGtKrkJGzmqMetCUBGBaE90Fzo9HRK35cR0JGoimJABT8TFz+g/xXMmwW2n8TYpQE5qSAJyozBUDGorEAADIoY8BAJBDYgAA5JAYAAA5JAYAQA6JAQCQQ2IAAOSQGAAAOf8P/yep0CHD21MAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Slope=530.7061985365614, intercept=0.6766240453921681, r_value=0.9969394089293704, p_value=1.6852816164419637e-22, std=9.547588506159919\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEeCAYAAACHXhKxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAHENJREFUeJzt3U9vG3d+x/HPN2gD7GXN2HHgHjb1MuktKDYU/QQSehe9pEBWjp9ATeXUU2zHRYDNIUBWzt53JfcJ2NK2QG5dy/sAqj8pihzag+QWuQTx2lYWBQJsi3x7+P1GGo2G5JDSaGY47xdAiJwZcn4zpOY7v//m7gIAtNcLVScAAFAtAgEAtByBAABajkAAAC1HIACAliMQAEDLEQgAoOUIBADQcgSCmjCzrpktm5mb2baZ3Uqt65nZmpk9N7PlKtM5KzNbMbOVqtNxEqnv6LmZbeesK+U7MrNFM9uNv41lM+vE5cP08kxatmNahqecll5Mw1Sfm/f7jo/l9DGdQvoa/zurhLvzqNFD0q6kWyPW5S6vII3DGd4zkNQrez9ndPy3JD3P+z6m/Y6KHqOkniSX1MksXwz/xse270oalPFdx+N/OOPnHft9SxrG5Z1TSNvUvzMeTo4AM1mY9g3uvuHuO2Xv54zsS7ohadnMujnrplHoGOO525f0Xt7+zGyQWd5z940p05InL33Tfo+TPFAIXNljm+RY2mb8nbUegQCFmVknZrvPz/C+gZn1ytzPWXL3dUnrkmYqhpjxGB9IupZZ1pG0kbP8RM7yO3D3JHg+K7L9qLRN+zvDoT+rOgGYTbwDXJa0p8OL0VVJu+6+mtl2WdKmwt1jJ17EFMtl78R1VxSy+xupz96S9FDS9fj3mcI/Xy/WYewn+4r/fMnd8VVJK6k7s/OSlhQuWlcLpH2Q3U/cd1IOfs3dd+LnPJR0191vjzhPiznvW5S0JmnV3ZfiPnZi+q5Iul/wrvKGpMdmNsye89T+c89x3jGO+oyUNUkPzayTungmy5cVznGyz4P18fVQ4XxLUtfd78Z1U3/Xqc9NciG5v7ui4u/zbvK7jMvG/Z5Gnbsjv7NJx46UqsumeBx9aIo6AoXy4V2FH7cU/gE8s822UmWmCmXbvdS+upl9d+Lzg3JbhfLpXmqfazlp25a0GJ93JT3PrB8oVa48Ke15+4mfsZtZNrGMfdz74nEuZrYdW8ac3md8//P0ecv5Pked49xzOWHfz1NpH8Tz1lGoPxgkn5vz3XRSrwcKF9b0MUzzXQ9iOkb+7ib8vlfiZw9j2vLK+if9nsal7WHmc0YeO4/woGio2fYlyd334t+kvDhpVdJT+GdN390u+OFd8cF7ox2FfxQp3BHuufu+u+/45Dvktz3e0SWfmVN+XjjteTyWeSdZ/3hH+mBCuia9b0+hrH9oZl2fsozZw53olqR72XUFzvEs0sVDnfj97GtE8VBy1+6pHEQ8H8PUuZ72u5akZ9N8dxm77r4ez93bCud/MbPNtL+nYwoeO0QdQR3tSbowxfbjKif7OswSSzpyUepK2o9lqoP4T5MUHxX57CPcPfmsYc4/9SxpH2VZoahFihfCgu9bUSw6UQiOSSDaiMuvStqNTRunvUgsSVrMqbAtco6ntSZpkJPGNUnvZYuFFO7w88re9xV+H+nX0zjJMRyI38MDHX6nB8tn+D1lFT321iMQ1M+ODstGD8S72WlbQ+zlfVZqneIdcPK461O0NIltw5O77IcKF9hVT5X1nob0fhQuGkmF4N6Yt2WtKlwou+n3mdkgHvs1dzeFO+up2sjH4Hpb4WKcVvgcZ45x3L42FC5k92JaEw8UimiWM5+/p/wK346mOH9F0zejfWV+p9P8nsak7VSOvQ0IBDXjodKzl5MNvj7NRTp+1oakvfSdavJP44cVxt3Uuk6Bf/Z0cOnGYqaeQvl0urIwKZ6a9W7u2H7iMSV3kMtTFuEkxScrmfPYy9zJ3y/wccdyDB4qIPckvZZaNukc5x5jAQ+UytXEfSXHd+TCF9PQyaRhUdJ6psgqz6zpm9ZuTGPym3lHk39PE9N2wmNvFVoN1dOCQrnprg5b+hxpFRMvJknQuKVwx5tkr5fNbDn+2JMy2K5CNnk/dSF8W9IdM9tMPtfd1+NnL0nqJ5+dKkrZMbOtVGueZNlqXLYXly/F9K2l0tq30CN1q0Daj+0nZUWzZe0/zXnfvqRu6gIzslVJPIe3FcqYL2S/E4Uy+mNl3co5x/HvuGMcZ02hEjRved4FbiGVhvMKv6dr8Zim+q6n+N0dES/y7ylcvK+a2U7yO3T3VTNL0rirEOhG/p4KpK2fask18thxyNyZsxjNYmaLp138BLQZRUNoBAvj0pyktQ2AEcgRoBFiEOhKB002AZwSAgEAtBxFQwDQcgQCAGi5RjQfffnll/3y5ctVJwMAGmV7e/sP7n5x0naNCASXL1/W1tZW1ckAgEYxs/8ush1FQwDQcgQCAGg5AgEAtByBAABajkAAAC1HIACAliMQAEDLldaPII47L0mv5YzbfjCfrnQ4PjsA4OyVkiOII0VuxFEiuyOGD74TA0CZU+ABACYoq2ioKym5+B+bNzfOBrUphSn+SpwCDwAwQSlFQ5nx4ns6Pg/sFemgeGgwampAAED5Sq0sjhf6nRF3/E+T5XkTnJvZMM5JuvXkyZMykwkArVZ2q6FBXkWxpKc6nGR7XzGHkObuq+7ed/f+xYsTB88DAMyotEBgZsOkyCepLDazTly9rsN6g45ifQEA4OyV2Wpo2cx2zex5atUjSXL3PUn7sUjoAs1HAaA6ZVUWb0h6KWf5Qup5UqFMEAAwVzY/X9GPdj7TK/5E39hFfdW7qSvvLFWdrJEaMTENADTF5ucremP7I/3A/iSZdElPdG77I21KtQ0GDDEBAKfoRzufhSCQ8gP7k36081lFKZqMQAAAp+gVz2/u/or/4YxTUhyBAABO0TeW39z9G3v5jFNSHIEAAE7RV72b+s5fPLLsO39RX/VuVpSiyQgEAHCKrryzpC8XPtHXuqjv3fS1LurLhU9qW1EsSebuVadhon6/71tbW1UnAwAaxcy23b0/aTuajwLACE3rDzArAgEA5Ghif4BZUUcAADma2B9gVgQCAMjRxP4AsyIQAECOJvYHmBWBAAByNLE/wKwIBACQo4n9AWZFPwIAc68tzUCz6EcAAGpXM9BZUTQEYK61qRnorAgEAOZam5qBzopAAGCutakZ6KwIBADmWpuagc6KQABgrrWpGeisaD4KAHOK5qMA5kpb+wKcBQIBgNqjL0C5qCMAUHv0BSgXgQBA7dEXoFwEAgC1R1+AchEIANQefQHKRSAAUHv0BSgX/QgAYE7RjwBALdEfoH5KKxoys2F8LE/Y7lZZaQBQL0l/gEt6ohdif4A3tj/S5ucrVSet1UoJBGY2kLTh7quSuvH1qO2ulpEGAPVDf4B6KitH0JWUXPz34msALUd/gHoqJRC4+2rMDUhST9Kxml4z67n7Rhn7B1BP9Aeop1Kbj5pZT9KOu+/krD5f5r4B1A/9Aeqp7H4EA3e/nV1YJDcQK5q3zGzryZP87CSAZqE/QD2V1o/AzIZJ8ZCZDdx9w8w67r5vZotxs/OSliTdGJFrkEQ/AgCYRdF+BGW2Glo2s10ze55a9UiS3H3d3dfjsk4ZaQAAFFNKh7JY7PNSzvKFzOtVSavZ7QDUHx3D5gc9iwFMjYli5guDzgGYGh3D5guBAMDU6Bg2XwgEAKZGx7D5QiAAMDU6hs0XAgGAqdExbL4wMQ0AzCkmpgEwEX0BIBEIgNaiLwASheoIzOyGmf3SzD4wsx+a2VtlJwxAuegLgETRHMGuu98zszfd/Y9mVmqiAJTvFX8i5fwr0xegfYq2GlqIuYAfm9lPJC1MegOAeqMvABJFA8GqpJ9Kuq4wxwB5R6Dh6AuARNFA4JIeSnog6bGZ/bq8JAE4C/QFQKJoHcFdhXmHkxLFC+UkB8BZuvLOkhQv/JfiA+1TNBCsufuj5IWZPSwpPQCAM1Y0EHTM7L6kPYVcwduSrpSWKgBToWMYTqJoIOhK+jD1+mkJaQEwAzqG4aSKVhZvu/vj5KFQcQygBugYhpMqmiP40MyWJT1TKBr6saS/Ki1VAAqjYxhOqmggWM5UFr9dUnoATOkbu6hLOj5j2Df2Mq2AUEihoiF3f2RmN83svpl9kA4KAKpFxzCcVNFB5/5O0o5ChfEXZvZBqakCUBgdw3BSRYuGHqdyAY8ZdA6oFzqG4SQKNx81M1foR9CV9KYkiocAYA4UCgRxCOqbkt5XGJL6TrnJAtqJjmGoQqFAYGY3JMnd3zOzc2b2rrv/U7lJA9qFjmGoStEOZbvJ0NPu/m2J6QFai45hqErROoIFM+vosI7giiRyBMApomMYqlK0H8FnCkNPvy+pSx0BcPqYMQxVKVo0JHe/5+7vu/uvzOyHZSYKaCM6hqEqRSuL31SYptIVMq9vSvpZiekCWufKO0valGKroT/oG3tZXy3QagjlK1pHMJC0knq9WEJagNajYxiqUDQQbMfhpyUVm6HMzIbx6Wvufnva9QCAs1HKMNRmNpC04e57ZrZmZgN33yi6HmgyOoWhacoahrobH6s6bHI6zXqgkegUhiYa2Woo3TIoZ9jpzXEf6u6r7r4aX/YkbU2zHmgqOoWhicblCO7ECevzXJc0sS+BmfUk7bj7zrTrYx3CUJJeffXVSbsCaoFOYWiicYHgmkKRTd6Y02+qQCCQNJhQETxyfcwxrEpSv9/3AvsCKsdsYWiisYHA3b/IWxH7FYxlZkN3vxufD9x9w8w67r4/av0M6Qdq5aveTZ1L6gii7/xFfbVwk0CA2hrXs/ilUStGBYhEbBW0bGa7ZvY8terRhPVAozFbGJrI3PNLXczsnKT3FHoTP6ty2Ol+v+9bW9QnA8A0zGzb3fuTthtZNBSHm74XP+xcnJPAFfoSbLj7H08rsQCA6hSdoexIUJB01cxeEkEBABqvaIeyAzEo/FY6CAp9Sb8/5XQBlaOHMNqi8DDUaUlnM3f/1t0JApg7SQ/hS3qiF2IP4Te2P9Lm5yuT3ww0DMNQAznG9hAmV4A5wzDUQA56CKNNihYNbbv74+QhaeIw1ECTMW0k2mRkjsDMfifp+eFLuyvpqQoMQw00HT2E0SbjioaWc0YdlVRoGGqg0Zg2Em0ysmfxyDeY/UTS3ln2HaBnMQBMr2jP4kJ1BGb2bvLc3f9NofIYADAHxrYaMrOfS7oqqW9mSZ54X2FWscrGHgIAnJ6xgcDdf2tmG5L6o+oLgLqjhzAw3sSioTikxJvpqSuBpqCHMDBZ0X4ERyqHY4UxUHvMIQxMVrRn8ftmtixpR4dDTNCPALVHD2FgsqKBYMXdf5u8oB8BmoI5hIHJChUNxUrjm2Z238w+oOIYTfFV76a+8xePLPvOX9RXvZsVpQion6L9CP5OoVjoQ0lfmNkHpaYKOCXMIQxMVrRo6HEqF/DYLKfQFaipK+8sHQwdfSk+ABwqGgi6ZuYKHcm6CpXFFA8BwBwYWTRkZr9Onrv7PUkLku5KuuruvzqDtAEAzsC4HIGZ2VsKje9c0qq70/galaGHMFCOkYHA3d9Pvzazy2a2IKknaSMOPgeciaSH8A/sT1LsIXxu+yNtSgQD4ITGFQ1dTv7GVkKr8fG6pPNnkTggQQ9hoDzjioY2zOwlhQri+5KW4jSVwJmjhzBQnnH9CDYUWgh9KOkLST9m4DlUhTmEgfKMDATu/r67f+vuj2Ifgi1JS2a2SYcynDV6CAPlGTd5/bsKvYmvSbou6ZxC34EPGWICZ405hIHyjJyz2My+l/RQoYhovcr6AeYsBoDpFZ2zeFxl8bX0iKMAgPk0rh/BiYKAmQ3j09fc/XbO+kWF+Y977n73JPsCAMyu6AxlUzGzgUKns1WFcYoGmfU9SXL3DUn7yWsAwNkrJRAoNDtNLv7JQHVp1xVyA8n6gdAKm5+v6OuPX9f3vzinrz9+nbmDgRooOvroVGJOINFT6JCW1pH0LPX6QhnpQL0wTARQT1PnCMzszaIdy2KRz46770ydMswdhokA6qnoDGXvJs/d/QsVL8oZ5FUUKxQLJeMVdSQ9zdnn0My2zGzryZPjc86ieV7x/O+RYSKAao0NBGb2czP7jaR/MLN/MbPfmdl9SVcmfbCZDZPWQEllsZl14ur7Oqw36Cr0VTjC3Vfdve/u/YsX84cXQLMwTARQT2MDQWxCelvSbXf/mbv/1N2vu/udce+LF/5lM9s1s+epVY/i5+6kttun6KgdGCYCqKeJlcXu/q2ZJU1AnyoMRd1399+Pec+GpJdyli+knq9m12O+MUwEUE9FWw3tuvs9M3vT3f/I5PWYFRPJA/VTNBAsxIt/J05ivyBpZI4AANAcRQPBqqQ7ChW7/8rcxQAwP8YNQ/1rhYpiKUxe/2ny3Mx+6O5/LDtxqC8mkgfmx7hWQxcU2vonf1+StCzpsaT3yk8a6irpIXxJT/RC7CH8xvZHDBcBNNS4QHDD3f8rzkPwkqR1hZxB193/8UxSh1qihzAwX8YNQ/2tmV2WdFdhdrJrTF4PiYnkgXkzMkdgZr9UmKHsN7Ez2ePUurfOInGoJ3oIA/NlXNHQQNL7kszM3jKzt5O/CnUFaCl6CAPzZVzz0RtxgLljzOzDktKDBqCHMDBfRk5eXydMXg8A0ys6eX1ZM5QBABqCQAAALUcgAICWK2XOYjQHQ0UAIBC0GJPJA5AoGmo1hooAIBEIWo3J5AFIBIJWY6gIABKBoNUYKgKARCBotSvvLOnLhU/0tS7qezd9rYv6cuETKoqBlmGICQCYUwwxAQAohEAAAC1HIACAliMQAEDLMcTEHGC8IAAnQSBoOMYLAnBSFA01HOMFATgpAkHDMV4QgJMiEDQc4wUBOKlSA4GZ9casWzSzgZkNy0zDvGO8IAAnVVogMLOBpLUR63qS9tx9Q9LeuICB8RgvCMBJldZqyN03zGxvzCbLkq5K6saAgBldeWdJihf+S/EBAEVVUkfg7jsKOYHnkp5VkQYAQFBJIDCzjqR9SZ9Kumdm3SrSAQCorkPZUNKn7r4fi48WJd1NbxArkYeS9Oqrr559CgGgJc40RxBzAke4+7pC7iC7fNXd++7ev3gxv4kkAODkSssRmNmipL6ZLcaLvSQ9krTg7nfN7FbMDZx399Wy0tEkjBkEoAplthpal7SeWbaQen732JtajDGDAFSFnsU1wZhBAKpCIKgJxgwCUBUCQU0wZhCAqhAIaoIxgwBUhUBQE4wZBKAq5u5Vp2Gifr/vW1tbVScDABrFzLbdvT9pO3IEANByBAIAaDkCAQC0HIEAAFquqtFH5xpjBgFoEgLBKWPMIABNQ9HQKWPMIABNQyA4ZYwZBKBpCASnjDGDADQNgeCUMWYQgKYhEJwyxgwC0DSMNQQAc4qxhgAAhRAIAKDlCAQA0HIEAgBoOYaYGIHxggC0BYEgB+MFAWgTioZyMF4QgDYhEORgvCAAbUIgyMF4QQDahECQg/GCALQJgSAH4wUBaBPGGgKAOcVYQwCAQggEANBypQYCM+uNW2dmi2a2WGYaAADjlRYIzGwgaW3MJnfcfV1Sd1zAAACUq7QhJtx9w8z28tbFXMBm3O5uWWkAAExWVR3BFUkXYvHQrTJ3tPn5ir7++HV9/4tz+vrj17X5+UqZuwOAxqmysvipu+9IBzmEU5cMHndJT/RCHDzuje2PCAYAkFLV6KNPJSXFRvsKOYT19AZmNpQ0jC//x8z+c9qd/PVfvLjw5/q/Y8v/997f69//9v3tgh/zsiQGGQIwi6qvH39ZZKMzDQRm1nH3fYWLfpIL6CjWF6S5+6qk1TNMXi4z2yrSIQMAsppy/Siz1dCipH6m2OeRJLn7nqT9uO5CbD0EAKhAI4aYqFJTIjqA+mnK9YOexZNVXjwFoLEacf0gRwAALUeOIMPMukn/BjPrjtvuLNMFoJmKXlOqRCA4rqfQtHVDhy2bjohDYnTi82F8LJ9dEgE0yMRrStWq6kdQOTPrSBqkFm24+37SgimOlTSqNVPf3VfjNhvuvmdma2Y2cPeNkpMOoIZOeE2pVGsDQao/wzHpC/yEj+nGx6pCxK9ltg9A+U7pmlKJVlcWxy/nmqT3JG1JWlHo6Xxb4cL+MNvHIRkpNRkeI7X8oaTb2eUA2mOWa0odtDZHELNx19x9KV7ElfqCxhXv9GOv5/Rn9STtEASA9jrBNaVyba4sHkpKKnjPK0TtWQ3c/fbJkwSgwU7zmnKm2hwILqTK6xbyKnmzE+bEpl97mWXDZE6FmC0E0E5TX1Pqos2BYCU2+1xUKMcr4kiroHjhXzazXTN7XkoqATTFLNeUWmh1ZXGeGLGvu/vt+HyQuuMfZusHAGCccdeUumhzjiBXarKcXnyevO5IelZl2gA0z6hrSp20ttXQBE8VhtCWQhMwKXQUqXXNP4Dayrum1AZFQyOkojcAnFidrykEAgBoOeoIAKDlCAQA0HIEAgBoOVoNAShF7ImfDM28XteRN0GOACOY2cDMnpvZMLXslpk9jH0qJr2/a2Zr5aZy5L4XUz08k2W92AN8Oa5fbvKQIPG7WIyPW3HZIBnsbMrPWizpu6r9hCwICATIFYfSeKCjneh2FEZXnDiYVrz7u1FS8ia5HnuAH/T7SHXkue/u63GQwEoCVVHpQJZZPlAY7XY9jm75mnTwnU090FlZwyLH9O0r5gjK2AdOB4EA46xJup563SkSBKSDYoFKJ+opkNZndZ1DNua6ro5YnV1X23Ft6j4hCwLqCDCSu2+MKjKId6vn43ar8fWSwkWpq3AHuKx4wcrZfqAwWceyQhHCRqor/i2Fu/euQq5kmLzOmQtiqDg7XPzcnqSumS2Ou9ONF9r9OM1oJ70PhVxQ+lh2MmndU7jzvqYwGdF+0ePL2dfeiPPQV+iJeuw43H3dzO6Y2bZCDufIuDVx39lzmj1PB+c4fU7jZ64oM8rupClYJ03IYma1nJAFkbvz4DHyofAPvaiQve/EZT1Jt+Lz7dS2u/Fvst3ahO2349+upOX4fChpMT6/pXCBHMTXy5m03UqtG6T2sTbiWNaS98T9JOk8to+cY0mndS2V1uSiW/T48vZ1bLv4+uGE72YQj2llwj6z5+mf0+c4dW4GCoEhm87uhHR0kjTE38pi1b9bHtM9KBrCJEnx0EGxkIe7zL14F5itQ5BnimTGbJ83iN+C4t2ohzvdrqROvNN/mtn2ig7vXPfi60l23H3D3VdT6czbR/ZY0mlN3y13pjy+vH1NNZhhUpwVj+OajhbB5X1W9jz9REfPcZKuJYWLuiR9KulqzCFMahzQ2AlZEBAIMJaHIoEjrWtiMcP5uO7gwjTKlNvvKl7YYjHKpqS9eLHNDgG+p8OLYDduO4tx+5hoyuObZl/78fOyk5kMMvuYdOHNnqf/0NFznGxzQ4cX9IG733b3BcXvf0xrscZOyIKAOgIUcaQFjsJFo5e0XonPu/Fvz0NZeC+1fNT23WR89ris4+53Y9NOSeGONTaVPB/3nW4JdDuuk6Re3DbZb3YSoZ7Cxe+amW2lcy3ZfSSflzqWQSatyYX4qsId+NqUx5fdV952+wq5jEUdH/X2WZI+hTvwT+NnZdOZfFb2PP1N+hyb2b5C8VaSphVJ+8l6SesxCGwrtlDKWInB8JlqXHGN0Rh0DkAhZtb1KVr/WAMmZEFA0RCAoiZ2JEzzBkzIgoCiIQCF+Gxj6dd6QhYEFA0BKJXVeEIWBAQCAGg56ggAoOUIBADQcgQCAGg5AgEAtByBAABajkAAAC33/9GGLTViVgNbAAAAAElFTkSuQmCC\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "# Plot pVar vs saving measures\n", - "plotReg(pVarList,savRteList,\n", - " xMin=pVarList[1]-0.0002,xMax=pVarList[-1]+0.0002,yMin=savRteList[1]-0.01,yMax=savRteList[-1]+0.01,\n", - " xLbl=r'Variance of Permanent Shocks, $\\sigma^{2}_{\\psi}$',\n", - " yLbl='Aggregate Saving Rate',\n", - " Title='Uncertainty vs Saving',\n", - " fileName='savRtevsPermShkVar'\n", - " )\n", - "plt.show(block=False)\n", - "plotReg(pVarList,KtoYList,\n", - " xMin=pVarList[1]-0.0002,xMax=pVarList[-1]+0.0002,yMin=1.7,yMax=KtoYList[-1]+0.1,\n", - " xLbl=r'Variance of Permanent Shocks, $\\sigma^{2}_{\\psi}$',\n", - " yLbl='Net Worth/Income',\n", - " Title='Uncertainty vs Net Worth Ratio',\n", - " fileName='BvsPermShkVar'\n", - " )\n", - "plt.ylabel('Net Worth/Income')\n", - "plt.xlabel(r'Variance of Permanent Shocks, $\\sigma^{2}_{\\psi}$')\n", - "plt.title('Uncertainty vs Net Worth Ratio',fontsize=16)\n", - "plt.xlim(pVarList[1]-0.0002,pVarList[-1]+0.0002)\n", - "plt.ylim(1.6,KtoYList[-1]+0.1)\n", - "plt.scatter(pVarList,KtoYList)\n", - "plt.xticks([pVarList[1],pVarList[-1]],[r'$\\bar{\\sigma}^{2}_{\\psi}/2$',r'$\\bar{\\sigma}^{2}_{\\psi}$'])\n", - "fileName='BvsPermShkVar'\n", - "if Generator:\n", - " plt.savefig(FigDir + nb_name + '-' + fileName + '.png')\n", - " plt.savefig(FigDir + nb_name + '-' + fileName + '.svg')\n", - " plt.savefig(FigDir + nb_name + '-' + fileName + '.pdf')\n", - "plt.show(block=False) " - ] - }, - { - "cell_type": "code", - "execution_count": 15, - "metadata": { - "code_folding": [ - 0 - ], - "scrolled": false - }, - "outputs": [], - "source": [ - "# Calculate variance of transitory shock vs saving measures\n", - "# Restore benchmark solution\n", - "EstimationEconomy.distributeParams( # Construct consumer types whose heterogeneity is in the given parameter\n", - " 'DiscFac',\n", - " Params.pref_type_count,# How many different types of consumer are there \n", - " center_estimate, # Increase patience slightly vs cstwMPC so that maximum saving rate is higher\n", - " spread_estimate, # How much difference is there across consumers\n", - " Params.dist_type) # Default is for a uniform distribution\n", - "EstimationEconomy.solve()\n", - "\n", - "savRteList_Tran = []\n", - "KtoYList_Tran = []\n", - "tVarList = []\n", - "tVarBase = BaselineType.TranShkStd[0] ** 2\n", - "for tVar in points * tVarBase:\n", - " tVarList.append(tVar) # Variance is std squared\n", - " savRteList_Tran.append(calcSavRte(EstimationEconomy,\"TranShkStd\",[tVar ** 0.5]))\n", - " KtoYList_Tran.append(0.25*np.mean(np.array(EstimationEconomy.KtoYnow_hist)[EstimationEconomy.ignore_periods:]))" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": { - "code_folding": [ - 0 - ] - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Slope=0.04425458972691994, intercept=0.22993889217509514, r_value=0.9997560905067558, p_value=6.235938075080433e-33, std=0.00022427988559214007\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAGnpJREFUeJzt3c+PG+d9x/HPt3ZaCEVtemUVSF24ErdtggBBrRV9ak42lbSHAoGztnqPteqtPSSS3UPhAAWUldtD0UOidf8BWbKbSw+pVgF66CXeXRXIpW67dHtwe5C1otwCLpCk3x7mGe3zjIbkkMvhkNz3C1gsOTPP8MtZ7nz5/JhnzN0FAEDuF5oOAAAwX0gMAIAEiQEAkCAxAAASJAYAQILEAABIkBgAAAkSAwAgQWIAACSebDoAHC9m1pa0Kakt6aK775lZS9I7kg4kXXH3fpMx5sxs391Xm44jZmYb4eFB+L3i7lsT7Gfu3hvmBzUGzJS79yTdlrTj7nthWV/SVUmbdSeF6MQ60qgT5zj7mgYzW5PUcvctd78laU/SRCd3kgKGITHg2Ag1k0sVt22b2fo09jVFK5JezJ+EJHtj3J2Mem8AiQFzx8y6ZrYffnfN7Ho4EefrN8LytfAtWmZ2OSzbCCe+eB83Q/m2pJaZreflotdbM7PN0NQlZU01m2bWKotH0gvxvsLyB1E8183s8oD39th2YR/5vjbLjou7b4cy+6FcN691lb2PITENfW+FY70Z9pkfcxLKMUAfA+aOu2+bWU9Sz9174cTWkbQdTkyr7r4VTuJXzGw3LydJZnbT3V8N+zhw91fDrvfMrB+aYWKXwvaSdCU874fyg+J5UlKyLzPbUvatXpJu5vGUvLfHtgvJ4IPQ5zLs2Lwa3ndX0nUz24z6GJL34e6XBsU04r3lx7qtrOlqz8y64XHx2GEJUWNAE3o6PFnlVnTYoZqLn+ffYl+U9IGUNaW4+yVJ5yT1oxrEB/k+42/UQ1wJCaczYruyeGJXw75ayt7jIPF2B+G9XJHUDkmutIkq/+Yf3vdW6CeIty17H4+9VtX3lvcHhdfdIykcHyQGNGFHWbNOrFWx4/kDRe3s4YR3W5LcfS8kgmGjdA5CuW70+0o46eXfpouxVdpXiL8naSOcVEuF7fphu72wjw13v+bueZIri6Gdv1Yo01J2LAe+j7LXGlc4ro/VfrC8aErCzIVmmouhvTs/gT468YRvqG1Jr5nZu5LOS+qZ2ba734o6T/vKmopuhTbwR8tCc0o7nHDjRHHTstFE+eseRK8pZTWXdtQn8ZqZ7agkHkl/V9iXJF3X40mvzFWl3+xbUfz7QxJLvp3C61wZ9j5CbMlrFY5v6XsLx7ov6ZKZXQj7351kaCwWj3EHN+DowrfzXugQXopv1yHpbYf31VLWrzHRdRNYLNQYgOlYN7NtDe9bWDQ7ktbiUUqKanZYXtQYAAAJOp8BAIlampKiTrQ1d79Wsj6fSmA1DNOL110uKxN79tln/fTp09MKFwCOhd3d3U/c/dSo7aaeGKKx1tth9Mha8epMHXZo3Yw768K685KGJobTp09rZ2dn2qEDwFIzs/+osl0dTUkXlNUWpKwjrltY346W9VRtaB8AYEbqSAzFKyxPxivDFZv5cLc1HV6gs7Ysw/wAYJE11vkcXWafNzMVp0gobr9hZjtmtnPv3r36AwSAY6qOxNDX4Um+Jen+gO26ecdzldpCqGl03L1z6tTIvhMAwITqGJV0Q4eX37d1OG/Lo7lw8nlhwuOussv828oSykqxwxoAMDtTrzFEk4J1lU1LnJ/g70TLN8Mc8A9CmVvRzI1ls1YCAGZkIa987nQ6znBVABiPme26+6jp5bnyGQCQIjEAABIkBgBAgsQAAEiQGAAACRIDACBBYgAAJEgMAIAEiQEAkCAxAAASJAYAQKKWez4DAA794O7HevuHH+o/+5/p11on9O2vfUFfP/vczMqPi8SAY6Hpf0zKH9/yP7j7sd58/yf67Kc/lyR93P9Mb77/E0mqtI+jlp/EQs6u+kuf/y3v/PH1hfpwUP7o/1zT+seUpBOfe0JXX/nyRP+YlKf8OOV/97s/0sf9zx5b/lzrhP7xjZdqLx9b+tlV86z5g7sfV9o+/+N+3P9MTvmFKn/U1377hx8m/9SS9NlPf663f/gh5Slfe/n/LDmpD1s+7fKTWNjEIC3Wh4Pyk5df9H9Myh/v8r/WOjHW8mmXn8RCJwZpcT4clJ+8/KL/Y1L+eJf/9te+oBOfeyJZduJzT+jbX/vCTMpPYuETw6J8OCg/eflF/8ek/PEu//Wzz+nqK1/Wc60TMmV9A1X7J6ZRfhJPvPXWW7XtvC5//pd//davvPB7OvG5J/Rnf/AlffHzT40sc/KXf1H/8C/39LP/O+xsp/xilD/qa3/x80/p1585oZ98/FD/878/03OtE/qzP/hS5X8sylP+KOXzfXzzK2f0J93f1je/cqbS53aa5XPf+c53/uutt97aGrUdo5IovxDlZz2OG1hGVUclLWRi6HQ6vrOz03QYALBQln64KgCgHiQGAECCxAAASJAYAAAJEgMAIEFiAAAkSAwAgASJAQCQIDEAABIkBgBAgsQAAEiQGAAACRIDACBBYgAAJEgMAIAEiQEAkHiyjp2a2bqkvqQ1d79Wsn4jPFx19yuDlgEAZm/qNQYzW5Mkd9+W1M+fR+u7krbdfUtS28y6ZcumHRcAoJo6mpIuKKstSFJPUvEk346W9cLzsmUAgAbU0ZTUknQQPT8Zrwy1gtyapBvuvldcVkNcAIAKGut8Dk1Me3FSKFsWrdswsx0z27l3794sQwWAY6WOxNCXtBIetyTdH7Bdt6STuWyZpKym4e4dd++cOnVqSqECAIrqSAw3dNhH0Ja0LUlm1so3MLONfLRS3tFctgwAMHtTTwx5M1A4ufejZqE70fJNM9s3sweDlgEAmlHLdQyFDuZ82bnwe1vSMyXFypYBAGaMK58BAAkSAwAgQWIAACRIDACABIkBAJAgMQAAEpWGq5rZRUmrkj6RtCWp4+4/qjMwAEAzql7HsO/u75jZWXf/1MxqDQoA0JyqieFcSAYtM3NJ5yRRYwCAJVQ1MWxJelPZ3Ec/dve36wsJANCkSonB3R9KekOSzOysmT3l7p/WGhkAoBGVRiWZ2Sv5Y3e/q8fvygYAWBJDawxm9g1J5yV1zOxSWNxXdvvN92uODQDQgKGJwd3fM7NtZcNT78woJgBAg0Y2Jbn7wzgpmNnpuGkJALBcxrnA7ZKy23SapF3RlAQAS6nyjXrcvWNmL7v7HTN7uc6gAADNqTpX0oGZfUvS02b2uqSzNcYEAGhQpcTg7u9JuuPu70v6ZUn/UmtUAIDGVJ5dNVy/IHf/K0m/WltEAIBGDUwM4QrnAzO7b2a/Y2ZnzOz7ZvZDZdc2AACW0LDO55fdfcXM2pK+L+ldSbcl9bmmAQCW17DE8FCS3L1nZpskAwA4HoYlhraZvRAen4keS9IFd3+zxrgAAA0ZlhheVTbNdn5Xnq9G684om4YbALBkhiWGS4Oaj8zsTE3xAAAaNnBU0rA+BXf/qJ5wAABNq3wdAwDgeKh6o57XC8+fNrPvmtlL9YQFAGhK1RrDM2Z2w8xOh+dvSLou6Zk6ggIANKdqYth39wvKRilJUjv0M/TrCQsA0JSq026v5s1JZvZReP6UpKdriwwA0IiqNYYtSQ/c/W8ktdy9o+zGPdQYAGDJVKoxuPtDSe+Fx3fN7Cl3f7vWyAAAjah6a8+zki5IcmVXQp+V9LUa4wIANKRqH0NX2Sik3HoNsQAA5kDVxLAbX+1sZrdrigcA0LCqieENM9uUdKCsKemMpN8atLGZrSvrmF5z92sl6zfCw1V3v1KlDABgNqomhuR+DGb28qANzWxNktx928zaZrbm7nvR+q6k7XCfh5vh+cGwMgCA2ak0XLVkQr39IZtf0OEw1p6y/olYO1rWC89HlQEAzMjAGoOZ3ZB0UdKqpE1JD/JVykYlDWpKainUAIKT8Up334qerkm6IencsDIAgNkZ1pT0hrt/amZ9ZfdmiDufzx71hUOT056775lZle03JG1I0vPPP3/UlwcADDDsfgwfRb+/W1h3d8g++5JWwuOWpPsDtuvmHc9Vyrj7lrt33L1z6tSpIS8PADiKqlNivGtmZ83spcK9n8vcUDTZnqRtSTKzVr6BmW3kI49C53NpGQDA7FVNDLdDLeEZSX9qZt8btGE+miic8PvR6KI70fJNM9s3swcjygAAZqzqcNU9M9uXdFPSxTB30kCFDuZ82bnwe1sl93EoKwMAmL2qieGKu79XayQAgLlQ9TqGR0nBzE6b2Sv1hQQAaFLV2VUvKrv/wn1l1zHsSnq/xrgAAA2p2pQkd++Y2cvufmfYlBgAgMVWdVTSgZl9S9LT4RafR77ADQAwn8bpY7jj7u8ra0oadoEbAGCBjdOUdDf8fqe+cAAATRtYYzCzb5vZ9/IrncPjvzezGxWufgYALKhhNYaeu78tZUlC0oq7fzU8f13SP80gPgDAjA3rY/Do8WtK7/n8QACApTSsxrBqZi9J+qokc/cfSdkFbkqTBgBgiQybdvttZSOQPnD3jiSZ2RlJ52cUGwCgAUNHJRVv6RnuzcCoJABYYlUvcAMAHBMkBgBAgsQAAEiMM7vqqqRPJG1J6uSjlAAAy6XqlBj77v6OmZ1190/NrNagAADNqZoYzoVk0DIzl3ROEjUGAFhCVRPDlqQ3JbUl/TifKgMAsHwqJQZ3fyjpDUkys7Nm9pS7f1prZACARlQalRTf4zlMv92tLSIAQKOG1hjM7BvKpsDomNmlsLgvqSfu+QwAS2nUlBjvmdm2suGpd4ZtCwBYDiObktz9YZwUzOx03LQEAFgu41zgdknSfWUzru6KpiQAWErj3PO5Y2Yvu/sdM3u5zqAAAM2pOlfSgZl9S9LT4baeZ2uMCQDQoEqJwd3fk3TH3d9X1pR0t9aoAACNqXodw+vh+gW5+zuSdszsu+HWnwCAJVK1KekZM7sR7vcsZVdBX5f0TB1BAQCaUzUx7Lv7BWVzJUlSO9zms19PWACAplQdlbQaOp1lZh+F509Jerq2yAAAjahaY9iSdODufyOp5e4dZdc1UGMAgCVTtcbwmrJaQlvSO2b2ElNvA8BymuQObg+5gxsALK9a7uBmZuvKmpnW3P3agG3W3H2vpEzb3bcqxgUAmLJx+hi+KukPJXWHNSOZ2Zokufu2pH7+vLBNV9LNQpleKNMrKwMAmI2qieGipKvu/pq7/8WIbS/osFO6p5Kb+uQJoLB4M/xuxzUJAMBsVU0MvXB7T0mSmb0wZNuWpIPo+clROw+JoGdmDwplAQAzVjUx/JGZ/Wu4+vldRc1A02BmLWW1jKvKRj21RxQBANSkaufz9TCRniRpxLTbfUkr4XFL2T0cRtlQ1lTVN7OepHVJSae1mW2E7fT8889XDBsAMK7KcyUVng+bRO+GoqkzJG1Lj2oFI7n7LZVcOOfuW+7ecffOqVOnKoYNABjX1CfRyzuOw8ijftSRHN8edF1SJ/xWGNK6YWbrZrbBcFUAaM44F7i9HWoI/64wid6gvoCyE7u7n4se35J0q7C+9HoHAMBsHXUSvUrNQwCAxTHOBW4PCpPobUh6UFtkAIBGVE0Mr+ajktz9bqgtPFtfWACApkza+fymss5nmpIAYMkc9Q5uD4eUAQAsIO7gBgBIHKXzmTu4AcASqlRjCBPoxZ3P3+IObgCwnKrWGGRmL5nZ98zs35R1PgMAltDAGkPoQ+gqu7/CeUk7kvaV3b2N2U8BYEkNa0r6U2WznN509xVJMrNXQrPS3VkEBwCYvYFNSe7+hrv/pqR3zex1M7so6UVp5I16AAALbGTns7vfVaghmNkZM/uGstlVX6w5NgBAA6pexyBJChe1fWRmDFMFgCVVeVRSzN3vjN4KALCIJkoMAIDlRWIAACRIDACABIkBAJAgMQAAEiQGAECCxAAASJAYAAAJEgMAIEFiAAAkSAwAgASJAQCQIDEAABIkBgBAgsQAAEiQGAAACRIDACBBYgAAJEgMAIAEiQEAkCAxAAASJAYAQOLJOnZqZuuS+pLW3P3agG3W3H0vfi6pLUnufquOuAAAo029xhBO8HL3bUn9/Hlhm66km4XFb4aE0C4rAwCYjTqaki4oqy1IUk9St7hBSBq9/HmoYXwQ1l2LaxIAgNmqIzG0JB1Ez09WKPOipJNmtmZml2uICQBQ0Tx1Pt/PawqhBpEwsw0z2zGznXv37s0+OgA4JupIDH1JK+FxS9L9CmXu67Bpqa+sBpFw9y1377h759SpU1MJFADwuDoSww2F0UXh97YkmVlrSJlbUZmWQn8DAGD2pp4YouagrqR+1JF8J98mNBV18iYjd+8pG8G0Lukkw1UBoDnm7k3HMLZOp+M7OztNhwEAC8XMdt29M2q7eep8BgDMARIDACBBYgAAJEgMAIAEiQEAkCAxAAASJAYAQILEAABIkBgAAAkSAwAgQWIAACRIDACABIkBAJAgMQAAEiQGAECCxAAASJAYAAAJEgMAIEFiAAAkSAwAgASJAQCQIDEAABIkBgBAgsQAAEiQGAAACRIDACBh7t50DGMzs/+W9GHTcQzxrKRPmg5iiHmOb55jk4jvqIjvaI4a32+4+6lRGz15hBdo0ofu3mk6iEHMbIf4JjPPsUnEd1TEdzSzio+mJABAgsQAAEgsamLYajqAEYhvcvMcm0R8R0V8RzOT+Bay8xkAUJ9FrTEAwEIzs3Uz65rZ5XHXx8vKthu171EaTwyTHBwz2wg/myO2O9LBmXJ8Zcs283VzEN9jsczL8TOzNTNzM9sPP9cHxTyj+Lrhp9bP3xRjm6fPXll88/TZS+Kr67NnZmuS5O7bkvr58yrrzawr6fyg7Ubtu4pGE8MkBycclG1335LUDn/EWg7OFON7bFnYxYaZ7UvqjRvbNOMri2Wejp+kFXc3d1+V9Kqk/KTS1PF7NSxbG/RZO+rxm2Js8/bZS+Iri6Xhz14xvlo+e5IuSOqHxz1J3THXD9uuatmBmq4xTHJw2tF2vfC8loMzxfjKlknSRXdfDR/ESUwrvrJY5ub4FY5Px93zf8aZHz9333b3S2FZ2933BuznqMdvWrHNzWdvQHxlsTTy2SuLr8bPXkvSQfT8ZJX1ZrZWeM2y7Ubte6SmL3Ab++C4+7Xo+ZqkG5LOleznyAdnWvFF/wBxzNLhN7i1QrmZxjcglrk5fvmTEN+70fqZH78olsuSLg3Z7qjHbyqxhZpCrtHPXll8A2Jp5LM3JL46PnuTWpnFizRdY5hYqObtFU66c6MsvuIyd78Wsv/JqIrfSHxNxlJmwN/3vLvn3/IajTmcDC6ZWWuWr1tFWWzz9Pcuxjdvn70Bf9tpf/b6OjzJtyTdH7W+pLYwaD+j9j1S04lh7IMTreu6+5Uh2x354EwxvseWWdYZuB6W39dhFX/m8Q2IZR6PX9wB18jxi/sPlDVBbAzYz1GP37RiyzX+2SuLb54+eyOO37Q/ezeicm1J22HfrSHr25Z1mG9IWgmxlm1Xuu9xNJ0YJjk4MrONvPoWsnUtB2eK8ZUt24liWg3Pm4qvLJZ5O37Ff76mjl9X6QmlN2C7ox6/acU2T5+9svjm6bM36PhN/bOX19zC36Mf1YzvDFrv7rfc/VYU36DtBu27skYTwyQHJzzetGzo2INB203j4EwrviExvxa+eew3GV9ZLPN0/CKPRoA0dfyUXXnaDt/aFP5Zp/75m1Zs8/TZK4tvnj57ZfFFu5zqZy/sZyt0eG9Fy84NWx8tX42aBcv2U1q2Kq58BgAkmm5KAgDMGRIDACDR9HUMAI6B0IHbUtbBe8sPLxTDHKLGgMosm97jgaVz2lw2s9tVxvObWdvMbtYb5cDXXi8MM8xj3wzrPKy/bGE+nCm+dq3vO8S8Hn4uh2VdM7s94f7Wa4h3TVkH7rak9RHbomEkBlQWLq55V+kVo3vK5pfpl5dKyvckXawpvFEuhBEa8dDHvrtfCaNPemEkxzVJu9N84eL7jpPTUYVRNfFQxtXwmts6nPJhLIXROFMR4usr1BimvX9MF4kB47qpbJ6ZXKtKUpAeNSdMcjHQ1BRiHTT+fJJx/QPF7zvUrM5PcffF/U21tjNNdjihH81Ic47EgLGEb6KlUwBEzTUb0fPbcROHDmenLNu+G7bvhuaR+GrTy2H5hpm14uclcWzE68J+2sVv6oPGn4frKR6LvWq8djh1cz6zbvy+O5I6hSatYrzJa4ff++F9J8084dt918x2zexy8T0NOJbF1xt2LHej7YvvqVTY5npodszfR1fSFWVTTdCUNO/cnR9+xvpR9q10XVmCaIVla5Iuh8e70bb74Xe+3c0R2++G321Jm+HxhqT18PiyspNsNzzfLMR2OVrXjV7j5oj3tD9ombJv5ePEG8fXLr6+pNsV4h143AbE31VWm7s+Irbi6/1tfGyjsjfD+vag9zQgjlYeQ/iMrDf9eeVn/B9qDJhE3pz0qBnJs2+qvfDNsNgHIS80Nw3ZPn6cO6dw5alnfQBtSa3wLbg4B86LOrxKtReeT+pR7GPGe1XSeTPbVZi6YIhB8RaP243wTTuZXTM0U8mzq1xfVdpUVxZb8fVeUHpsc21lM4zm8Vd9Txs6rB2taMJ+DjSLxICxeUlzUmiGWAnryuaX0RG231faRv+Bss7ifBqDWHzPgXbY9sjGjLfrWaf2OZU3u/XDPvKROiPj9azJ6JIevzFMtxDLqBNx8fX+Wemxjbe7qMOT/GPvacBItJN+2IdwzqPZQEMT0/qopig0j+sYMKniCJ+esrtedZV9210LJ6w1y6YL3gsnwnz5oO3bYbtuWNZy92uWDSuVlH2zDe3i+bfnR3G4+5WwTgpz5Uev2/XCtMXh5PZaeN3LCmPsQ1yPYh8nXkkv5rFKuhW/73DS7IVv/9sD4i2+du62P95xe5Bvq+wb+tXwvroDjmXx9X4/PrYh3q6yprOVsK/ryu5yFm/TUjZ6a7UQz/WQRA8UdYSH7c+H19/UZBPjYUaYKwmYc1FifSyxNSlKdFW2TTq2fcLJ3TAbNCUB8++Cma3PU1IIxrlJ0WpIBqtK74SGOUSNAUDtouGyK3OY4FBAYgAAJGhKAgAkSAwAgASJAQCQIDEAABIkBgBAgsQAAEiQGAAAif8HV5mSJcS5+AwAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Slope=0.24992735556376564, intercept=2.554920083325562, r_value=0.9999887157190909, p_value=1.3035255504050127e-45, std=0.0002723909104974136\n" - ] - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEYCAYAAABRB/GsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAGPJJREFUeJzt3c9vHOd9x/HPt3EOvjgryTw0KVJliSI52sv1qbdmWaO3wqGkfyBa+Q+ISeumAAUcyumph5R0/gFFcppLD6koA71WJF0gl7oAV+kh6UERtU6AumgafHOYZ8jhcGZ39sfsr+f9Agju/P7us7vznXlm5nnM3QUAiNefzDsAAMB8kQgAIHIkAgCIHIkAACJHIgCAyJEIACByJAIAiByJAAAiRyIAgMiRCCJkZk0ze2hmR2bWCuMaYdyemTXmHWPKzE7mHUNWKLu9bFyh7PYmKTsz65jZSzPbzWznsZmdmFkzjGuFz2xrwvfQqlKume/JYzPbMrOumW2n35kxt71QnycSRhMTcTKzrqQNd7+TGdeS1Hf3Xt3bdvf9RVvXCNvsSGpJWk/LL91ZDyq7YbGGJPDC3e+H4Zakh+6+nt22ux+MEfOFbZvZY3ffrLDclqTNzPtsSHrm7ldG3SYWF2cEmKmwI7kzdEadHZGWHv2Osq5pCzvrdkgKQ1WM9YGkW5nhq5L62bM2SaejxjrlcroqaehZT9E2h32emB8SAQqFqoqT8L+Tr/YI1QSdUM2Q7qi2w7hu+NFn1/EwLN+U1AhVDa3c9lpmtpseXSvZ6e2GqpdL8Uh6I7uuTPVKGs+emW2XvLdL84V1pOvarVBMNyTtlZTfhbIoe99Z7n4c5mlmRmeTQzvMk26jmymPrcx7G6XML322Ba6GMulK2pW0kXuvRZ9d0TbPPs+y+DEfr8w7ACwmdz8ws56knrv3wo+5Lekg/GjX3X0//PB3zOwoXU6SzOyhu98I6zh19xth1cdm1nf3R7lN3gnzS9JOGO6H5cvieUVJVdbZusxsX8lRq5RUq1yqRgnrujRf2Pk/dffjEMewMuqFneiekh1kGkO3pCyK3nfeI0lbZvZIUi/8HYUyOZPuODPb2DWzXqacqpT5pc+2JKbTkIDScunnphd9dpe2mf08B8R/LMwcZwTx6ul8R5i6qstVD9nh9KjxLUlPpWRnGOqPNxSqMcKO5Wm6zoo/7p2wc2gPma8onqwPwroaSt5jmex8p+G97EhqhqRWqSolrSKSlK0iKiuLKtIzgGYo256kniVVUNn3s5kbfqHzsqta5sPKssxOfrjiZ5c1KH7MGIkgXodKTt+zGu6eP9or8lRJMpB0Vh/8WEqqN8JOaNBFwtOwXCfzfyccPaZHiPnYKq0rxN+T1B104TbM1w/zHYd1dN39vrunO/KqMdxQ5oxA5WVxIdaSuI6V7JSzO+YHknZz7+dIFz+/dSWfaZmh267oNN1uWmWn8s9u0DZHjR81IhFEKuwIb4e67K1wRHdWNRCOZJuSboYd/aakTTNrhB/9i7BcR8nR66Ow3Nm4dB1pVUnGw9y4dIfRUrIDvJpdPsRQGI+kfy5Y/54Gnw2kPtDFao5GJv6TokSSqQo7u/YQ5tvJDF8qi5L3XWZPF6tp9pUkgzPhbpz02smWpKNQHTO0zAd9trn3mk5r2/ldUY+UJMn0fRV+dkXvN7vdsvgrlA1qwO2jWBlm1gx13mPdYgnEijMCrJL0DpVan4MAVg1nBAAQOc4IACByJAIAiNxSPFD2+uuv+/Xr1+cdBgAslaOjo9+4+9qw+ZYiEVy/fl2Hh9xiDACjMLP/qjIfVUMAEDkSAQBEjkQAAJEjEQBA5JbiYvEvfvW5/vIHn+i9t7+pv33za5WX+9mnv9KHP/9Mv+5/oa82XmX5GS6/zLGzPMsv+/Kj+tK9e/dqW/m0/N3f/8M9/1ZH//qfz/VnV17Vt/70taHL/OzTX+nuT3+h0//5P0nS7/73/1l+Rssvc+wsz/LLvnzW97///f++d+/e0O5Cl6pq6Ivf/0Ef/vyzSvN++PPP9MXv/8Dyc1h+mWNneZZf9uXHsVSJQJJ+3f9iovlYvv7llzl2lmf5ZV9+HLUlgtAfaddK+n5Nx1dsn/3MVxuvTjQfy9e//DLHzvIsv+zLj6OWRBA6rTgInU80S3oo6prZiUZoMvjVL39J7739zUrzvvf2N/Xql7/E8nNYfpljZ3mWX/blx1HXXUPN8LevZEdf1OXf7QodeZ/52ohXztP5xr3yzvLjL7/MsbM8yy/78uOovT8CM3uspE/T49z4bUnHklqhA/BS7XbbaWsIAEZjZkfu3h42X60Xi0NvUcdFfZGGTsIPJF0rqjoK1xcOzezw+fPndYYJAFGr+66hjrvv5EeGnfxWGHyhgqojd99397a7t9fWhraiCgAYU613DaVVPukRv5k1wuRDSWnn4uthGAAwB3XeNbRrZidm9jIz6Ykkhaqim+Gs4KSo6ggAMBu13DUU6v6vFIzfyLwe+tgzAKB+S/dkMQBgukgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJEjEQBA5EgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJEjEQBA5EgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAELlX6lqxmXXDy3V33ymYviWpL6nl7vfrigMAMFgtZwRm1pF04O77kpphODu9JUnufiCpnw4DAGavrqqhpqR0598Lw1m3lJwNpNM7AgDMRS1VQ+FMINWS9CA3S0PSaWb4Wh1xAACGq/VicajyOXb34zGW7ZrZoZkdPn/+vIboAABS/XcNdYouFCupFroaXjckvcjP4O777t529/ba2lqdMQJA1GpLBGbWTe8GSi8Wm1kjTH6g8+sGTUkHdcUBABiszruGds3sxMxeZiY9kaS0qijM1x+n6ggAMB11XSw+kHSlYPxG5vV+fjoAYPZ4shgAIkciAIDIkQgAIHIkAgCIHIkAACJX6a4hM7staV3SbyTtS2q7+yd1BgYAmI2qt4+euPtHZvamu//WzGoNCgAwO1UTwUbY+TfMzCVtSOKMAABWQNVEsC/prqRvSHrq7h/WFxIAYJaqXix2SY8l/UTSMzP7UX0hAQBmqeoZwX1Jh5LSiwP0HwAAK6JqInjo7k/SATN7XFM8AIAZq5oIGmb2QEm3kibp25Leqi0qAMDMVE0ETUnvZ4YvdSQDAFhOVRPBkbs/SweoGgKA1VE1EbxvZrtKOpw3JbeR/kVtUQEAZqZqItjNXSz+dk3xAABmrFIicPcnZvaepLaSB8p+WG9YAIBZqfRAmZl9V9KxkgvGn5rZ92qNCgAwM1Wrhp5lqoae0egcAKyOyrePhsbmekpuJX1T0pPBiwAAlkGlqiF3/0hJi6P3JW1yjQAAVscoHdPI3W+a2VfM7B13/2m9oQEAZmGUjmk+kSR3/5xrBACwOkbpmKah82sEb0nijAAAVkDVawQfKml6+l1JTXe/W2tUAICZqXpGkF4w/kiSzOw1d/9tbVEBAGam6sXiNyXdUtJTmSm5ffTtGuMCAMxI1TOCjqS9zPBWDbEAAOaAZqgBIHI0Qw0AkaMZagCIXOnto2b2Wvo6mwSCp7VFBACYqUHPEdw1szeK/iRVeo7AzFoDpu2G/93RQgYATNOgqqEbSp4iLmpP4k0NSQZmlt5ptF4yS9fMtiTdqRAnAKAmAxOBu39aNCE8VzCQux+YWW/ALLfd/dGw9QAA6jWoauhK2YSyBDGippl1zGx7CusCAIxp0BnBUWh+2iWdTrvZaXe/L0lmtmlmHXc/mOb6AQDVlCYCd/9c520LfSWbFCQdTNLWULhAfBqqhl4ouRYBAJiDSs8R5JOCpE0zu6IRk4KZNdy9L+lQSZPWUnIxea9g3q6kriR9/etfr7J6AMAYKjVDneXun7v7x+7+YyX9FreL5gt3BLXD/9STsI5jSTfDtJMwnN/Ovru33b29trY2apgAgIoqN0OdlTZDHc4UPimaJ1T7PMqN28i83h9n2wCA6aIZagCIHM1QA0DkaIYaACJXmgjM7F8kvTwftPtKbvWkGWoAWCGDzgh2C1odlUQz1ACwSkpvHx2QBN4QzVADwMqo9ByBmb2Tvnb3f1dy8RgAsAIGXiw2s+9I2lTyYFjaXHRfyVPBU217CAAwHwMTgbt/bGYHktplVUUAgOU2tGooPD38ZrbrSgDA6qja1lAv27BcuGAMAFgBVR8oezf0MXys8yYmeI4AAFZA1USw5+4fpwM8RwAAq6NqfwQfm9l7SpqcfuruP6w3LADArFR9juC7SqqF3pf0qZl9r9aoAAAzU7Vq6Fnm9tFnZlZXPACAGauaCJpm5koeJGsquVjMcwUAsAJKq4bM7Efpa3f/SNKGpPuSNrlGAACrY9AZgZnZXym5XdQl7bv7h7MJCwAwK6WJwN3fzQ6b2XUz25DUknQQGp8DACy5QR3TXHf3X5rZdSVdU/61kusDB0ruIAIArIBBVUMHZnZFyQXiB5LuZLurBACshoGJQNKOkofIJOkbZvYi2+YQAGD5VblG8ESSQuujd8zspqQH3DkEAKth0DWCd5RcC7gh6ZakryhJCu/TNwEArI5BVUOPJD1WUkV0g+sDALCaBiWCG9kWRwEAq6n0yWKSAADEoWoPZQCAFUUiAIDIkQgAIHIjJwIzezM8UwAAWAFVeyh7J33t7p9K6tQWEQBgpgZ2TGNm35G0KaltZneUNEn9Ukn7Qz+tPzwAQN0GJoLQaf2BpPY4TxObWcvdC1sqNbMtSX1JLXe/P+q6AQDTMbRqyN0/V9JV5Qdm9j0zey10WDOQmXUkPSyZ1grrPpDUT4cBALNX9WLxibvflfSkauujYSffK5l8S8nZgMI8XHMAgDmp2nn9hplJUiN0Yr8h6ZMJttuQdJoZvjbBugAAE6iaCPYl3VXSQ9m/0XcxAKyOQc1Q/0hJxzRS0nn9B+lrM3ttwg5q+pKuhtcNSS8mWBcAYAKDrhFcU7KzTv9fkbQr6Zmkm+NszMwa4eUDJWcX0nk/yPl5u2Z2aGaHz58/H2dzAIAKBiWC2+7+y9APwRUl/RO4pKa7/3jYisPtoe3wP/VEktJbSsOdRf2iW0zdfd/d2+7eXltbq/6OAAAjGdRV5edmdl3SfSW9k43UOY27P1KSPLLjNjKv90cNFgAwfaVnBGb2AyU9lP2ju7+dTQJVniMAACyHQXcNdSS9K8nCjt+UVA2ZpB9Ieqv+8AAAdRuUCG6HBuYuMbP3a4oHADBjg7qqLEwCYdrI7Q4BABYTHdMAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJEjEQBA5EgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJEjEQBA5EgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJGrLRGY2ZaZdcxsu2T6bvjfrSsGAMBwtSQCM2tJkrsfSOqnwzldMzuR1KsjBgBANXWdEdyS1A+ve5I6BfPcdvf1kCwAAHNSVyJoSDrNDF8rmKc5qOoIADAbc7tY7O73w9nANTO7dMZgZl0zOzSzw+fPn88hQgCIQ12JoC/panjdkPQiOzHs5LfC4AtJzfwK3H3f3dvu3l5bW6spTABAXYnggc537k1JB5JkZo0w7jAdJ2k9DAMA5qCWRODux5IUqnz66bCkJ5npN8NZwUlmOgBgxl6pa8Xuvl8wbmPQdADA7PFkMQBEjkQAAJEjEQBA5EgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJEjEQBA5EgEABA5EgEARI5EAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAEDkSAQBEjkQAAJEzd593DEOZ2e8kfTbvOAZ4XdJv5h3EAIsc3yLHJhHfpIhvMpPG9+fuvjZsplcm2MAsfebu7XkHUcbMDolvPIscm0R8kyK+ycwqPqqGACByJAIAiNyyJIL9eQcwBPGNb5Fjk4hvUsQ3mZnEtxQXiwFg2ZnZlqS+pJa73x9lupltp+OK5hu27mFmfkZgZltm1jGz7arTzawb/naHzDdw3TOOr2jcbjptAeK7FMuilJ+ZtczMzewk/O2VxTyD2Drhb1G/e0XxLdJ3ryi+RfruXYivxu9eS5Lc/UBSPx2uMt3MOpI2y+Ybtu4qZpoIximMUAgH7r4vqRk+tFoKY4rxXRoXVtE1sxNJvVFjm2Z8RbEsUvlJuuru5u7rkm5ISnciY5ffBLHdCONaZd+zOZddPr5F++5diK8olkUqP9Xw3QtuKTliV1hHZ8Tpg+arumypWZ8RjFMYzcx8vTBcS2FMMb6icZJ0293XwxdvHNOKryiWhSm/XPm03T398U1SfiPH5u4H7n4njGu6+3HJeuZSdiXxLcx3ryS+olgWpvxq+u5JUkPSaWb4WpXpZtbKbbNovmHrHmrWzxGMXBi5+q6WpAeSNgrWM3FhTCu+zBc+G7N0foQ2Vj3etOIriWVhyi8dCPH9JDN9kvIb64cY4tiWdGfAfHMpu6L4wplAaq7fvaL4SmJZmPLLjJvmd28SV2exkWW5ayg9tTvO7WQXRlF8+XHufj9k92uZU/a5xDfPWIqUfL6b7p4exc0t5vDDv2NmjVltcxRF8S3S552Pb9G+eyWf77S/e32d79Qbkl4Mm15wNlC2nmHrHmrWiWDkwshM67j7zoD5Ji6MKcZ3aZwlF++2wvgXOj9ln3l8JbEsYvllL5hNWn5j/RAzdc09Sd2S9cyl7EriS839u1cU3yJ994aU3zS/e1JyZpYu15R0ENbdGDC9ackF7q6kqyHWovkK1z2KWSeCcQpDZtbN3CbVKZlv4sKYYnxF4w4zMa2H4XnFVxTLopVf/sc2afmNE1tHF3cevZL55lV2RfEt0nevKL5F+u6Vld+0v3tKz8zC59HPnPk+KZvu7o/c/VEmvrL5ytZd2UwTwTiFEV7vWnIr18uy+aZRGNOKb0DMN8ORxck84yuKZZHKL+PsDo1Jy2+c2JQ8zNMMR2QKP8yF+e4VxbdI372i+Bbpu1cUX2aVU/vuZdazHy5Q72fGbQyanhm/nqnmK1pP4bJV8UAZAERuaS4WAwDqQSIAgMiRCAAgcsvSMQ2AJRPuvmkouTvnkZ8/pYsFwxkBClnSZtJLu9gw2LaZPa7yYJWZNc3sYb1Rlm57K3fvd9qY2ImZ7Ybpu4vwMNO4wmexFf62w7iOmT0eY11bNX1WLSV33xxI2hoyL+aIRIBC4YnGn+jiY/nHShrp6hcvdWH5nqTbNYU3zK1wG93Z/ejh1rtjJU2APAoPW80lUVWVTWS58R0lTw2n95mvS+cNqo26ndxtk1MT4usrnBHUsQ1MB4kAgzxU0lhXqlElCUhn1QLjPIE5NRViPS14eGghhLOuzZLJ+Wl79Uc0HjtvDZVqoQXGNQKUcveDsiqDcLR6Ncy3H4bvKNkpNZUcAe7qvB31/PwdSTthnpaSnUX60M+2kqP3ppKzkm46nH9gJlRdpa2W7ofH8JtmtjXoSDfsaPvu3guvz7ah5Cwo+16Oc7H2lBx535C04+79qu+vYFu9knJoS2oXvY/w4NhdMztScoaT78SkU1Cm+XI6K+NsmYZ17inX3LIPaXUzbPOGpJtKnrzdC2W0I6lnZo/rOvPAFLg7f/yV/in5QW8pOb1vhHEtSdvh9VFm3pPwP53v4ZD5j8L/pqTd8LoraSu83layg+yE4d1cbNuZaZ3MNh6WvJeH6TJhO2mcl7ZR8F6ysT7MxJrudKu+v6JtXZovDD8e8tl0wnvaG7LNfDn9U7aMM2XTUZIY8nE2h8TRSGMI35WteX9v+Rvtj6ohDJNWD51VC3lylNkLR4H5awjyXJXMgPmzr1MbCkejnhzpNiU1wpF+viGxt3R+5NoLw8Mce3gUPxNn0Tby7yUba/ZouTHi+yvaVtF8pdLqrPA+buhiFVzRuvLl9IYulnEa1x2FNm0kfSBpM5whDLs5oKvzDlyuaozrFJgvEgEG8qRK4MLdNaGa4WqYVtRIlyaY/0RhxxaqUZ5K6vl52zBZ2Y5XmmHecQzaxlAjvr9RttUP68v32NXJbWPYjjdfTv+hi2WcznNb5zv0jrvveNIWTic3b941P78GsOEF1UgF7wELhGsEqOLCHThKdhppV3/H4XUz/G95Uhfeyowvm78Z5uuEcQ13vx9u7ZSUHLGGWyXTViKzdwLthGlS6DAks91OdoeUXjuQdMPMDj3X1nx2G+n6Mu+lk4s13RFvKjkCfzji+8tvq2i+vpKzjC1dbo3zNI1PyRH4B2Fd+TjTdeXL6W+yZWxmfSXVW2lMe0q6c0y39ygkgSOFO5Ry9kIyPNUCX7hGORqdA1CJmTV9hLt/QkK6FRJRS8lZxix790JFVA0BqGqkHtr8/I6llp8/x4EFRNUQgEp8vHb4Xyi5DVYar0MczABVQwBqlTkjwIIiEQBA5LhGAACRIxEAQORIBAAQORIBAESORAAAkSMRAEDkSAQAELk/AmQJviVVIxUtAAAAAElFTkSuQmCC\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "# Plot transitory variance versus saving measures\n", - "plotReg(tVarList,savRteList_Tran,\n", - " xMin=tVarList[1]-0.001,xMax=tVarList[-1]+0.001,yMin=savRteList[1]-0.01,yMax=savRteList[-1]+0.01,\n", - " xLbl=r'Variance of Transitory Shocks, $\\sigma^{2}_{\\theta}$',\n", - " yLbl='Aggregate Saving Rate',\n", - " Title='Uncertainty vs Saving',\n", - " fileName='savRteVSTranShkVar'\n", - " )\n", - "plt.show(block=False)\n", - "plotReg(tVarList,KtoYList_Tran,\n", - " xMin=tVarList[1]-0.001,xMax=tVarList[-1]+0.001,yMin=savRteList[1]-0.01,yMax=KtoYList[-1]+0.1,\n", - " xLbl=r'Variance of Permanent Shocks, $\\sigma^{2}_{\\psi}$',\n", - " yLbl='Net Worth/Income',\n", - " Title='Uncertainty vs Net Worth Ratio',\n", - " fileName='BvsTranShkVar'\n", - " )\n", - "plt.show(block=False) " - ] - } - ], - "metadata": { - "cite2c": { - "citations": { - "6202365/7MR8GUVS": { - "DOI": "10.3982/QE694", - "URL": "https://onlinelibrary.wiley.com/doi/abs/10.3982/QE694", - "abstract": "In a model calibrated to match micro- and macroeconomic evidence on household income dynamics, we show that a modest degree of heterogeneity in household preferences or beliefs is sufficient to match empirical measures of wealth inequality in the United States. The heterogeneity-augmented model's predictions are consistent with microeconomic evidence that suggests that the annual marginal propensity to consume (MPC) is much larger than the roughly 0.04 implied by commonly used macroeconomic models (even ones including some heterogeneity). The high MPC arises because many consumers hold little wealth despite having a strong precautionary motive. Our model also plausibly predicts that the aggregate MPC can differ greatly depending on how the shock is distributed across households (depending, e.g., on their wealth, or employment status).", - "accessed": { - "day": 5, - "month": 2, - "year": 2019 - }, - "author": [ - { - "family": "Carroll", - "given": "Christopher" - }, - { - "family": "Slacalek", - "given": "Jiri" - }, - { - "family": "Tokuoka", - "given": "Kiichi" - }, - { - "family": "White", - "given": "Matthew N." - } - ], - "container-title": "Quantitative Economics", - "id": "6202365/7MR8GUVS", - "issue": "3", - "issued": { - "year": 2017 - }, - "language": "en", - "note": "Citation Key: carrollDistributionWealthMarginal2017", - "page": "977-1020", - "page-first": "977", - "title": "The distribution of wealth and the marginal propensity to consume", - "type": "article-journal", - "volume": "8" - } - } - }, - "jupytext": { - "formats": "ipynb,py:percent" - }, - "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.9" - }, - "varInspector": { - "cols": { - "lenName": 16, - "lenType": 16, - "lenVar": 40 - }, - "kernels_config": { - "python": { - "delete_cmd_postfix": "", - "delete_cmd_prefix": "del ", - "library": "var_list.py", - "varRefreshCmd": "print(var_dic_list())" - }, - "r": { - "delete_cmd_postfix": ") ", - "delete_cmd_prefix": "rm(", - "library": "var_list.r", - "varRefreshCmd": "cat(var_dic_list()) " - } - }, - "types_to_exclude": [ - "module", - "function", - "builtin_function_or_method", - "instance", - "_Feature" - ], - "window_display": false - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/notebooks/Uncertainty-and-the-Saving-Rate.py b/notebooks/Uncertainty-and-the-Saving-Rate.py deleted file mode 100644 index 931cb426..00000000 --- a/notebooks/Uncertainty-and-the-Saving-Rate.py +++ /dev/null @@ -1,418 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.2' -# jupytext_version: 1.2.3 -# kernelspec: -# display_name: Python 3 -# language: python -# name: python3 -# --- - -# %% [markdown] -# # Uncertainty and Saving in Partial Equilibrium -# -# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/econ-ark/DemARK/master?filepath=notebooks%2FUncertainty-and-the-Saving-Rate.ipynb) -# -# Saving rates vary widely across countries, but there is no consensus about the main causes of those differences. -# -# One commonly mentioned factor is differences across countries in the degree of uncertainty that individuals face, which should induce different amounts of precautionary saving. -# -# Uncertainty might differ for "fundamental" reasons, having to do with, say, the volatility of demand for the goods and services supplied by the country, or might differ as a result of economic policies, such as the strucutre of the social insurance system. -# -# A challenge in evaluating the importance of precautionary motives for cross-country saving differences has been a lack of consensus about what measures of uncertainty ought, in principle, to be the right ones to look at in any attempt to measure a relationship between uncertainty and saving. -# -# This notebook uses [a standard model](https://econ.jhu.edu/people/ccarroll/papers/cstwMPC) to construct a theoretical benchmark for the relationship of saving to two kinds of uncertainty: Permanent shocks and transitory shocks to income. -# -# Conclusions: -# 1. The model implies a close to linear relationship between the variance of either kind of shock (transitory or permanent) and the saving rate -# 2. The _slope_ of that relationship is much steeper for permanent than for transitory shocks -# * Over ranges of values calibrated to be representative of microeconomically plausible magnitudes -# -# Thus, the quantitative theory of precautionary saving says that the principal determinant of precautionary saving should be the magnitude of permanent (or highly persistent) shocks to income. -# -# (Because the result was obtained in a partial equilibrium model, the conclusion applies also to attempts to measure the magnitude of precautionary saving across groups of people who face different degrees of uncertainty within a country). -# -# @authors: Derin Aksit, Tongli Zhang, Christopher Carroll - -# %% {"code_folding": [0, 11]} -# Boring non-HARK setup stuff - -Generator = True # This notebook can be used as a source for generating derivative notebooks -nb_name = 'Uncertainty-and-the-Saving-Rate' - -# This is a jupytext paired notebook that autogenerates BufferStockTheory.py -# which can be executed from a terminal command line via "ipython BufferStockTheory.py" -# But a terminal does not permit inline figures, so we need to test jupyter vs terminal -# Google "how can I check if code is executed in the ipython notebook" - -from IPython import get_ipython # In case it was run from python instead of ipython -def in_ipynb(): - try: - if str(type(get_ipython())) == "": - return True - else: - return False - except NameError: - return False - -# Determine whether to make the figures inline (for spyder or jupyter) -# vs whatever is the automatic setting that will apply if run from the terminal -if in_ipynb(): - # %matplotlib inline generates a syntax error when run from the shell - # so do this instead - get_ipython().run_line_magic('matplotlib', 'inline') -else: - get_ipython().run_line_magic('matplotlib', 'auto') - print('You appear to be running from a terminal') - print('By default, figures will appear one by one') - print('Close the visible figure in order to see the next one') - -# Import the plot-figure library matplotlib - -import matplotlib.pyplot as plt - -# In order to use LaTeX to manage all text layout in our figures, we import rc settings from matplotlib. -from matplotlib import rc -plt.rc('font', family='serif') - -# LaTeX is huge and takes forever to install on mybinder -# so if it is not installed then do not use it -from distutils.spawn import find_executable -iflatexExists=False -if find_executable('latex'): - iflatexExists=True - -plt.rc('font', family='serif') -plt.rc('text', usetex=iflatexExists) - -# The warnings package allows us to ignore some harmless but alarming warning messages -import warnings -warnings.filterwarnings("ignore") - -import os -from copy import copy, deepcopy - -# Define (and create, if necessary) the figures directory "Figures" -if Generator: - nb_file_path = os.path.dirname(os.path.abspath(nb_name+".ipynb")) # Find pathname to this file: - FigDir = os.path.join(nb_file_path,"Figures/") # LaTeX document assumes figures will be here -# FigDir = os.path.join(nb_file_path,"/tmp/Figures/") # Uncomment to make figures outside of git path - if not os.path.exists(FigDir): - os.makedirs(FigDir) - -from copy import deepcopy -from scipy.optimize import golden, brentq -from time import time -import numpy as np -import scipy as sp - -# %% {"code_folding": [0]} -# Import HARK tools and cstwMPC parameter values -from HARK.utilities import plotFuncsDer, plotFuncs -from HARK.ConsumptionSaving.ConsIndShockModel import PerfForesightConsumerType -import HARK.cstwMPC.cstwMPC as cstwMPC -import HARK.cstwMPC.SetupParamsCSTW as Params - -# Double the default value of variance -# Params.init_infinite['PermShkStd'] = [i*2 for i in Params.init_infinite['PermShkStd']] - -# %% {"code_folding": [0]} -# Setup stuff for general equilibrium version - -# Set targets for K/Y and the Lorenz curve -lorenz_target = cstwMPC.getLorenzShares(Params.SCF_wealth,weights= - Params.SCF_weights,percentiles= - Params.percentiles_to_match) - -lorenz_long_data = np.hstack((np.array(0.0),\ - cstwMPC.getLorenzShares(Params.SCF_wealth,weights=\ - Params.SCF_weights,percentiles=\ - np.arange(0.01,1.0,0.01).tolist()),np.array(1.0))) -KY_target = 10.26 - -# %% {"code_folding": [0]} -# Setup and calibration of the agent types - -# The parameter values below are taken from -# http://econ.jhu.edu/people/ccarroll/papers/cjSOE/#calibration - -Params.init_cjSOE = Params.init_infinite # Get default values of all parameters -# Now change some of the parameters for the individual's problem to those of cjSOE -Params.init_cjSOE['CRRA'] = 2 -Params.init_cjSOE['Rfree'] = 1.04**0.25 -Params.init_cjSOE['PermGroFac'] = [1.01**0.25] # Indiviual-specific income growth (from experience, e.g.) -Params.init_cjSOE['PermGroFacAgg'] = 1.04**0.25 # Aggregate productivity growth -Params.init_cjSOE['LivPrb'] = [0.95**0.25] # Matches a short working life - -PopGroFac_cjSOE = [1.01**0.25] # Irrelevant to the individual's choice; attach later to "market" economy object - -# Instantiate the baseline agent type with the parameters defined above -BaselineType = cstwMPC.cstwMPCagent(**Params.init_cjSOE) -BaselineType.AgeDstn = np.array(1.0) # Fix the age distribution of agents - -# Make desired number of agent types (to capture ex-ante heterogeneity) -EstimationAgentList = [] -for n in range(Params.pref_type_count): - EstimationAgentList.append(deepcopy(BaselineType)) - EstimationAgentList[n].seed = n # Give every instance a different seed - -# %% {"code_folding": [0]} -# Make an economy for the consumers to live in - -EstimationEconomy = cstwMPC.cstwMPCmarket(**Params.init_market) -EstimationEconomy.print_parallel_error_once = True # Avoids a bug in the code - -EstimationEconomy.agents = EstimationAgentList -EstimationEconomy.act_T = Params.T_sim_PY # How many periods of history are good enough for "steady state" - -# %% {"code_folding": [0]} -# Uninteresting parameters that also need to be set -EstimationEconomy.KYratioTarget = KY_target -EstimationEconomy.LorenzTarget = lorenz_target -EstimationEconomy.LorenzData = lorenz_long_data -EstimationEconomy.PopGroFac = PopGroFac_cjSOE # Population growth characterizes the entire economy -EstimationEconomy.ignore_periods = Params.ignore_periods_PY # Presample periods - -#Display statistics about the estimated model (or not) -EstimationEconomy.LorenzBool = False -EstimationEconomy.ManyStatsBool = False -EstimationEconomy.TypeWeight = [1.0] - -# %% {"code_folding": [0]} -# construct spread_estimate and center_estimate if true, otherwise use the default values -Params.do_param_dist=True # Whether to use a distribution of ex-ante heterogeneity - -# Discount factors assumed to be uniformly distributed around center_pre for spread_pre on either side - -spread_pre=0.0019501105739768 #result under the default calibration of cjSOE -center_pre=1.0065863855906343 #result under the default calibration of cjSOE - -do_optimizing=False # Set to True to reestimate the distribution of time preference rates - -if do_optimizing: # If you want to rerun the cstwMPC estimation, change do_optimizing to True - # Finite value requires discount factor from combined pure and mortality-induced - # discounting to be less than one, so maximum DiscFac is 1/LivPrb - DiscFacMax = 1/Params.init_cjSOE['LivPrb'][0] # - param_range = [0.995,-0.0001+DiscFacMax] - spread_range = [0.00195,0.0205] # - - if Params.do_param_dist: # If configured to estimate the distribution - LorenzBool = True - # Run the param-dist estimation - paramDistObjective = lambda spread : cstwMPC.findLorenzDistanceAtTargetKY( - Economy = EstimationEconomy, - param_name = Params.param_name, - param_count = Params.pref_type_count, - center_range = param_range, - spread = spread, - dist_type = Params.dist_type) # Distribution of DiscFac - t_start = time() - - spread_estimate = golden(paramDistObjective - ,brack=spread_range - ,tol=1e-4) - center_estimate = EstimationEconomy.center_save - t_end = time() - else: # Run the param-point estimation only - paramPointObjective = lambda center : cstwMPC.getKYratioDifference(Economy = EstimationEconomy, - param_name = Params.param_name, - param_count = Params.pref_type_count, - center = center, - spread = 0.0, - dist_type = Params.dist_type) - t_start = time() - center_estimate = brentq(paramPointObjective # Find best point estimate - ,param_range[0] - ,param_range[1],xtol=1e-6) - spread_estimate = 0.0 - t_end = time() - - print(spread_estimate) - print('****************') - print(center_estimate) - print('****************') -else: # Just use the hard-wired numbers from cstwMPC - center_estimate=center_pre - spread_estimate=spread_pre - - -# %% {"code_folding": [0]} -# Construct the economy at date 0 -EstimationEconomy.distributeParams( # Construct consumer types whose heterogeneity is in the given parameter - 'DiscFac', - Params.pref_type_count,# How many different types of consumer are there - center_estimate, # Increase patience slightly vs cstwMPC so that maximum saving rate is higher - spread_estimate, # How much difference is there across consumers - Params.dist_type) # Default is for a uniform distribution - - -# %% {"code_folding": [0]} -# Function to calculate the saving rate of a cstw economy -def calcSavRte(Economy,ParamToChange,NewVals): - ''' - Calculates the saving rate as income minus consumption divided by income. - - Parameters - ---------- - Economy : [cstwMPCmarket] - A fully-parameterized instance of a cstwMPCmarket economy - ParamToChange : string - Name of the parameter that should be varied from the original value in Economy - NewVals : [float] or [list] - The alternative value (or list of values) that the parameter should take - - Returns - ------- - savRte : [float] - The aggregate saving rate in the last year of the generated history - ''' - for NewVal in NewVals: - if ParamToChange in ["PermShkStd","TranShkStd"]: - ThisVal = [NewVal] - else: - ThisVal = NewVal # If they asked to change something else, assume it's a scalar - - for j in range(len(Economy.agents)): # For each agent, set the new parameter value - setattr(Economy.agents[j],ParamToChange,ThisVal) - cstwMPC.cstwMPCagent.updateIncomeProcess(Economy.agents[j]) - - Economy.solve() - - C_NrmNow=[] - A_NrmNow=[] - M_NrmNow=[] - for j in range (len(Economy.agents)): # Combine the results across all the agents - C_NrmNow=np.hstack((C_NrmNow,Economy.agents[j].cNrmNow)) - A_NrmNow=np.hstack((A_NrmNow,Economy.agents[j].aNrmNow)) - M_NrmNow=np.hstack((M_NrmNow,Economy.agents[j].mNrmNow)) - CAgg=np.sum(np.hstack(Economy.pLvlNow)*C_NrmNow) # cNrm times pLvl = level of c; sum these for CAgg - AAgg=np.sum(np.hstack(Economy.pLvlNow)*A_NrmNow) # Aggregate Assets - MAgg=np.sum(np.hstack(Economy.pLvlNow)*M_NrmNow) # Aggregate Market Resources - YAgg=np.sum(np.hstack(Economy.pLvlNow)*np.hstack(Economy.TranShkNow)) # Aggregate Labor Income - BAgg=MAgg-YAgg # Aggregate "Bank Balances" (at beginning of period; before consumption decision) - IncAgg=(BaselineType.Rfree-1)*BAgg+YAgg # Interest income plus noninterest income - savRte=(IncAgg-CAgg)/IncAgg # Unspent income divided by the level of income - return savRte - - -# %% {"code_folding": [0]} -# Function to plot relationship between x and y; x is the parameter varied and y is saving rate -def plotReg(x,y,xMin,xMax,yMin,yMax,xLbl,yLbl,Title,fileName): - # Result_data_path = os.path.join(Folder_path,'SavingVSPermShr_Youth_MPC_15.png') - plt.ylabel(yLbl) - plt.xlabel(xLbl) - plt.title(Title) - plt.xlim(xMin,xMax) - plt.ylim(yMin,yMax) - plt.scatter(x,y) - # Draw the linear fitted line - m, b = np.polyfit(x, y, 1) -# plt.plot(x, m*np.asarray(x) + b, '-') - if Generator: - plt.savefig(FigDir + nb_name + '-' + fileName + '.png') - plt.savefig(FigDir + nb_name + '-' + fileName + '.svg') - plt.savefig(FigDir + nb_name + '-' + fileName + '.pdf') - slope, intercept, r_value, p_value, std_err = sp.stats.linregress(x,y) - print('Slope=' + str(slope) + ', intercept=' + str(intercept) + ', r_value=' + str(r_value) + ', p_value=' + str(p_value)+', std=' + str(std_err)) - - -# %% {"code_folding": [0]} -# Proportion of base value for uncertainty parameter to take (up to 1 = 100 percent) -# Do not go above one to avoid having to worry about whether the most patient consumer violates the -# Growth Impatience Condition (https://econ.jhu.edu/people/ccarroll/papers/BufferStockTheory/#GIC) -bottom=0.5 -points=np.arange(bottom,1.+0.025,0.025) - -# %% {"code_folding": [0]} -# Calculate variance of permanent shock vs saving measures -savRteList = [] -KtoYList = [] -pVarList = [] -pVarBase = BaselineType.PermShkStd[0] ** 2 -for pVar in points * pVarBase: - pVarList.append(pVar) # Variance is square of standard deviation - pStd = pVar ** 0.5 -# print(pStd) - savRteList.append(calcSavRte(EstimationEconomy,"PermShkStd",[pStd])) - KtoYList.append(0.25*np.mean(np.array(EstimationEconomy.KtoYnow_hist)[EstimationEconomy.ignore_periods:])) - -# %% {"code_folding": [0]} -# Calculate how much net worth shrinks when permanent variance is halved -ShrinksBy = KtoYList[1]/KtoYList[-1] -print('Halving the magnitude of the permanent variance causes target wealth to fall to %1.3f' % ShrinksBy) -print('of its original value.') - -# %% {"code_folding": [0]} -# Plot pVar vs saving measures -plotReg(pVarList,savRteList, - xMin=pVarList[1]-0.0002,xMax=pVarList[-1]+0.0002,yMin=savRteList[1]-0.01,yMax=savRteList[-1]+0.01, - xLbl=r'Variance of Permanent Shocks, $\sigma^{2}_{\psi}$', - yLbl='Aggregate Saving Rate', - Title='Uncertainty vs Saving', - fileName='savRtevsPermShkVar' - ) -plt.show(block=False) -plotReg(pVarList,KtoYList, - xMin=pVarList[1]-0.0002,xMax=pVarList[-1]+0.0002,yMin=1.7,yMax=KtoYList[-1]+0.1, - xLbl=r'Variance of Permanent Shocks, $\sigma^{2}_{\psi}$', - yLbl='Net Worth/Income', - Title='Uncertainty vs Net Worth Ratio', - fileName='BvsPermShkVar' - ) -plt.ylabel('Net Worth/Income') -plt.xlabel(r'Variance of Permanent Shocks, $\sigma^{2}_{\psi}$') -plt.title('Uncertainty vs Net Worth Ratio',fontsize=16) -plt.xlim(pVarList[1]-0.0002,pVarList[-1]+0.0002) -plt.ylim(1.6,KtoYList[-1]+0.1) -plt.scatter(pVarList,KtoYList) -plt.xticks([pVarList[1],pVarList[-1]],[r'$\bar{\sigma}^{2}_{\psi}/2$',r'$\bar{\sigma}^{2}_{\psi}$']) -fileName='BvsPermShkVar' -if Generator: - plt.savefig(FigDir + nb_name + '-' + fileName + '.png') - plt.savefig(FigDir + nb_name + '-' + fileName + '.svg') - plt.savefig(FigDir + nb_name + '-' + fileName + '.pdf') -plt.show(block=False) - -# %% {"code_folding": [0]} -# Calculate variance of transitory shock vs saving measures -# Restore benchmark solution -EstimationEconomy.distributeParams( # Construct consumer types whose heterogeneity is in the given parameter - 'DiscFac', - Params.pref_type_count,# How many different types of consumer are there - center_estimate, # Increase patience slightly vs cstwMPC so that maximum saving rate is higher - spread_estimate, # How much difference is there across consumers - Params.dist_type) # Default is for a uniform distribution -EstimationEconomy.solve() - -savRteList_Tran = [] -KtoYList_Tran = [] -tVarList = [] -tVarBase = BaselineType.TranShkStd[0] ** 2 -for tVar in points * tVarBase: - tVarList.append(tVar) # Variance is std squared - savRteList_Tran.append(calcSavRte(EstimationEconomy,"TranShkStd",[tVar ** 0.5])) - KtoYList_Tran.append(0.25*np.mean(np.array(EstimationEconomy.KtoYnow_hist)[EstimationEconomy.ignore_periods:])) - -# %% {"code_folding": [0]} -# Plot transitory variance versus saving measures -plotReg(tVarList,savRteList_Tran, - xMin=tVarList[1]-0.001,xMax=tVarList[-1]+0.001,yMin=savRteList[1]-0.01,yMax=savRteList[-1]+0.01, - xLbl=r'Variance of Transitory Shocks, $\sigma^{2}_{\theta}$', - yLbl='Aggregate Saving Rate', - Title='Uncertainty vs Saving', - fileName='savRteVSTranShkVar' - ) -plt.show(block=False) -plotReg(tVarList,KtoYList_Tran, - xMin=tVarList[1]-0.001,xMax=tVarList[-1]+0.001,yMin=savRteList[1]-0.01,yMax=KtoYList[-1]+0.1, - xLbl=r'Variance of Permanent Shocks, $\sigma^{2}_{\psi}$', - yLbl='Net Worth/Income', - Title='Uncertainty vs Net Worth Ratio', - fileName='BvsTranShkVar' - ) -plt.show(block=False)