correlation_.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import pearsonr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "year_start = 1998\n",
    "year_end = 2018\n",
    "years_iter = range(year_start,year_end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read cleaned data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## MSA data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = {}\n",
    "for year in years_iter:\n",
    "    df[year] = pd.read_csv('data/cleaned/msa_naicssector_cleanedv2_totEmploy_'+str(year)+'.csv')\n",
    "    \n",
    "# used later to rank the MSA's by employment\n",
    "df_2017 = df[2017].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## State-country data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_state = {}\n",
    "for year in years_iter:\n",
    "    df_state[year] = pd.read_csv('data/cleaned/state_naicssector_cleanedv2_totEmploy_'+str(year)+'.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merge data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## MSA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for year in years_iter:\n",
    "    df[year] = df[year].set_index(['MSA','NAICS Code'],drop=True)\n",
    "    \n",
    "tot_empl_all = df[year_start]\n",
    "for year in range(year_start+1,year_end):\n",
    "    tot_empl_all = pd.merge(tot_empl_all, df[year],how='inner' ,on=['MSA','NAICS Code'],copy=False )\n",
    "    #df_all.to_csv('data/draft/df_all_merge_'+str(year)+'.csv')\n",
    "    \n",
    "tot_empl_all = tot_empl_all.T\n",
    "tot_empl_all.index = list(years_iter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## state-country"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "for year in years_iter:\n",
    "    df_state[year] = df_state[year].set_index(['State','NAICS Code'],drop=True)\n",
    "\n",
    "tot_empl_state_all = df_state[year_start]\n",
    "for year in range(year_start+1,year_end):\n",
    "    tot_empl_state_all = pd.merge(tot_empl_state_all, df_state[year],how='inner' ,on=['State','NAICS Code'],copy=False )\n",
    "    \n",
    "tot_empl_state_all = tot_empl_state_all.T\n",
    "tot_empl_state_all.index = list(years_iter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate PPC: point percent change"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## MSA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_size = 2\n",
    "\n",
    "tot_empl_ppc = tot_empl_all.pct_change().drop( tot_empl_all.index[0] )*100\n",
    "tot_empl_ppc.loc[2003] = (tot_empl_ppc.loc[2002]+tot_empl_ppc.loc[2004])/2\n",
    "#tot_empl_ppc = tot_empl_ppc.rolling(window=window_size,center=False).mean().drop( tot_empl_ppc.index[:window_size-1] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## State-Country"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "tot_empl_state_ppc = tot_empl_state_all.pct_change().drop( tot_empl_state_all.index[0] )*100\n",
    "tot_empl_state_ppc.loc[2003] = (tot_empl_state_ppc.loc[2002]+tot_empl_state_ppc.loc[2004])/2\n",
    "#tot_empl_state_ppc = tot_empl_state_ppc.rolling(window=window_size,center=False).mean().drop( tot_empl_state_ppc.index[:window_size-1] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Utilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## MSA to State "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "MSA_df = pd.read_csv('data/MSA_wikipedia.csv',delimiter='\\t')\n",
    "\n",
    "MSAs = MSA_df['Metropolitan statistical area ']\n",
    "MSA_shorts = MSAs.map(lambda x: x.split(',')[0].split('-')[0].lower() )\n",
    "MSA_shorts_dict = {MSA_shorts.iloc[i]:\n",
    "                   {'MSA': MSAs.iloc[i][:-1],\n",
    "                    'State': MSAs.iloc[i].split(', ')[1].split(' ')[0] } for i in range(len(MSAs))  }\n",
    "\n",
    "\n",
    "state_abbr_df = pd.read_csv('data/state_abbr.csv',delimiter='\\t',names=['State','Abbr.'])\n",
    "state_to_abbr_dict = { state_abbr_df['State'].iloc[i].title()[0:-1] :state_abbr_df['Abbr.'].iloc[i] for i in range(len(state_abbr_df)) }\n",
    "state_to_abbr_dict['United States'] = 'US'\n",
    "\n",
    "abbr_to_state_dict ={ v:k for k,v in state_to_abbr_dict.items()}\n",
    "\n",
    "MSA_to_state_dict = { MSAs.iloc[i][:-1]: MSAs.iloc[i].split(', ')[1].split(' ')[0]  for i in range(len(MSAs))  }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Compare MSA data to state and country data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ppc_relative( MSA,NAICS_code ):\n",
    "    return (tot_empl_ppc[ (MSA,NAICS_code) ]-tot_empl_state_ppc[ ( abbr_to_state_dict[MSA_to_state_dict[MSA]] ,NAICS_code) ],\n",
    "            tot_empl_ppc[ (MSA,NAICS_code) ]-tot_empl_state_ppc[ ( 'United States' ,NAICS_code) ] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## NAICS Codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-11-48819ea87482>:1: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.\n",
      "  NAICS_df = pd.read_csv('data/NAICS_codes.csv',delimiter='    ')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'11': 'Agriculture, Forestry, Fishing and Hunting',\n",
       " '21': 'Mining, Quarrying, and Oil and Gas Extraction',\n",
       " '22': 'Utilities',\n",
       " '23': 'Construction',\n",
       " '31-33': 'Manufacturing',\n",
       " '42': 'Wholesale Trade',\n",
       " '44-45': 'Retail Trade',\n",
       " '48-49': 'Transportation and Warehousing',\n",
       " '51': 'Information',\n",
       " '52': 'Finance and Insurance',\n",
       " '53': 'Real Estate and Rental and Leasing',\n",
       " '54': 'Professional, Scientific, and Technical Services',\n",
       " '55': 'Management of Companies and Enterprises',\n",
       " '56': 'Administrative and Support and Waste Management and Remediation Services',\n",
       " '61': 'Educational Services',\n",
       " '62': 'Health Care and Social Assistance',\n",
       " '71': 'Arts, Entertainment, and Recreation',\n",
       " '72': 'Accommodation and Food Services',\n",
       " '81': 'Other Services (except Public Administration)',\n",
       " '92': 'Public Administration',\n",
       " '99': 'Unclassified',\n",
       " '--': 'Total'}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NAICS_df = pd.read_csv('data/NAICS_codes.csv',delimiter='    ')\n",
    "NAICS_dict = { NAICS_df['NAICS Code'].iloc[i]:NAICS_df['Description'].iloc[i] for i in NAICS_df.index }\n",
    "NAICS_dict['--'] = 'Total'\n",
    "NAICS_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## MSA rank by total Employment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "msa_emp_rank_df =  df_2017[ df_2017['NAICS Code']=='--' ].sort_values(by='Total Employment - 2017',ascending=False).reset_index(drop=True)\n",
    "del msa_emp_rank_df['NAICS Code'], msa_emp_rank_df['Total Employment - 2017']\n",
    "msa_emp_rank_df_dict = {i:msa_emp_rank_df['MSA'].iloc[i] for i in range(len(msa_emp_rank_df))  }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Print Info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# House data\n",
    "- ## Read, clean and Calculate PPC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "house_df = pd.read_csv('data/other-data/Metro_zhvi_uc_sfrcondo_tier_0.33_0.67_sm_sa_mon.csv')\n",
    "del house_df['RegionID'], house_df['SizeRank'], house_df['RegionType'], house_df['StateName']\n",
    "house_df.set_index('RegionName',drop=True,inplace=True)\n",
    "house_df = house_df.T\n",
    "house_df.index = pd.to_datetime(house_df.index)\n",
    "\n",
    "house_df = house_df.rolling(12,center=True).mean()\n",
    "\n",
    "house_df['Year'] = house_df.index.year\n",
    "house_df = house_df.groupby(house_df['Year']).mean()\n",
    "house_df = house_df.loc[ list( range(year_start,2020) ) ]\n",
    "\n",
    "# calculate PPC\n",
    "house_ppc_df = house_df.diff()/house_df.shift()*100\n",
    "house_ppc_df = house_ppc_df.drop( house_ppc_df.index[0] )\n",
    "#window_size = 2\n",
    "#house_ppc_df = house_ppc_df.rolling(window=window_size,center=False).mean().drop( house_ppc_df.index[:window_size-1] )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Some more utilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def house_ppc_relative(MSA_rank):\n",
    "    return house_ppc_df[house_ppc_df.columns[MSA_rank]]-house_ppc_df[house_ppc_df.columns[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "houseRank_to_empRank = {1:None, 2:0, 3:1, 4:2, 5:4, 6:5, 7:3, 8:8, 9:7,\n",
    "                        10:6, 11:9, 12:None, 13:17, 14:12, 15:11, 16:10, 17:14, 18:15, 19:18 ,\n",
    "                        20:16, 21:13, 22:19, 23:None, 24:22, 25:31, 26:26, 27:20, 28:24, 29:None, 30:23 }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def available_industries(MSA):\n",
    "    code_list = list(tot_empl_all[MSA].columns)\n",
    "    return pd.Series( [ NAICS_dict[ key ] for key in code_list ] , index=code_list )\n",
    "\n",
    "def print_info(house_rank):\n",
    "    print( \"MSA: \", house_ppc_df.columns[ house_rank ] )\n",
    "    if houseRank_to_empRank[house_rank] != None: \n",
    "        print( \"Available Industires:\\n\", available_industries( msa_emp_rank_df_dict[ houseRank_to_empRank[house_rank] ] ) )\n",
    "    else:\n",
    "        print(\"Employment data no available at the moment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Compare PPC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Compare Relative PPC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSA:  Los Angeles-Long Beach-Anaheim, CA\n",
      "Available Industires:\n",
      " --                                                   Total\n",
      "11              Agriculture, Forestry, Fishing and Hunting\n",
      "21           Mining, Quarrying, and Oil and Gas Extraction\n",
      "22                                               Utilities\n",
      "23                                            Construction\n",
      "31-33                                        Manufacturing\n",
      "42                                         Wholesale Trade\n",
      "44-45                                         Retail Trade\n",
      "48-49                       Transportation and Warehousing\n",
      "51                                             Information\n",
      "52                                   Finance and Insurance\n",
      "53                      Real Estate and Rental and Leasing\n",
      "54        Professional, Scientific, and Technical Services\n",
      "55                 Management of Companies and Enterprises\n",
      "56       Administrative and Support and Waste Managemen...\n",
      "61                                    Educational Services\n",
      "62                       Health Care and Social Assistance\n",
      "71                     Arts, Entertainment, and Recreation\n",
      "72                         Accommodation and Food Services\n",
      "81           Other Services (except Public Administration)\n",
      "99                                            Unclassified\n",
      "dtype: object\n"
     ]
    }
   ],
   "source": [
    "MSA_rank = 2\n",
    "print_info(MSA_rank)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poverty and income data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Read and calculate PPC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "poverty_df = pd.read_csv( 'data/saipe/poverty_percent_county_years_.csv',header=[0,1] ,index_col=0)\n",
    "income_df = pd.read_csv('data/saipe/income_county_years_.csv',header=[0,1],index_col=0)\n",
    "\n",
    "window_size = 2\n",
    "\n",
    "poverty_df.dropna(axis=1,inplace=True)\n",
    "poverty_ppc_df = poverty_df.pct_change()*100\n",
    "poverty_ppc_df = poverty_ppc_df.drop( poverty_ppc_df.index[0] )\n",
    "#poverty_ppc_df = poverty_ppc_df.rolling(window=window_size,center=False).mean().drop( poverty_ppc_df.index[:window_size-1] )\n",
    "\n",
    "income_df.dropna(axis=1,inplace=True)\n",
    "income_ppc_df = income_df.pct_change()*100\n",
    "income_ppc_df = income_ppc_df.drop( income_ppc_df.index[0] )\n",
    "#income_ppc_df = income_ppc_df.rolling(window=window_size,center=False).mean().drop( income_ppc_df.index[:window_size-1] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relative_ppc_poverty(MSA_rank):\n",
    "    return -poverty_ppc_df[ houseRank_to_county[MSA_rank] ] + poverty_ppc_df['US','United States']\n",
    "\n",
    "def relative_ppc_income(MSA_rank):\n",
    "    return income_ppc_df[ houseRank_to_county[MSA_rank] ] - income_ppc_df['US','United States']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "houseRank_to_county = {1:('US','United States'), 2:('CA','Los Angeles County'), 3:('IL','Cook County'), 4:('TX','Dallas County'), 5:('PA',''), \n",
    "                       6:('TX','Harris County'), 7:('DC','District of Columbia'), 8:('FL','Miami-Dade County'), 9:('GA','Fulton County'), 10:('MA','Suffolk County'), \n",
    "                       11:('CA','San Francisco County'), 12:('MI','Wayne County'), 13:('CA','Riverside County'), 14:('AZ','Maricopa County'), 15:('WA','King County'), \n",
    "                       16:('',''), 17:('CA','San Diego County'), 18:('',''), 19:('FL','Hillsborough County') ,  20:('MD','Baltimore city'), \n",
    "                       21:('CO','Denver County'), 22:('PA','Allegheny County'), 23:'NaN', 24:22, 25:('CA','Sacramento County'),\n",
    "                       26:('TX','Bexar County'), 27:('FL','Orange County')\n",
    "                        }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ## Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSA:  Los Angeles-Long Beach-Anaheim, CA\n",
      "Available Industires:\n",
      " --                                                   Total\n",
      "11              Agriculture, Forestry, Fishing and Hunting\n",
      "21           Mining, Quarrying, and Oil and Gas Extraction\n",
      "22                                               Utilities\n",
      "23                                            Construction\n",
      "31-33                                        Manufacturing\n",
      "42                                         Wholesale Trade\n",
      "44-45                                         Retail Trade\n",
      "48-49                       Transportation and Warehousing\n",
      "51                                             Information\n",
      "52                                   Finance and Insurance\n",
      "53                      Real Estate and Rental and Leasing\n",
      "54        Professional, Scientific, and Technical Services\n",
      "55                 Management of Companies and Enterprises\n",
      "56       Administrative and Support and Waste Managemen...\n",
      "61                                    Educational Services\n",
      "62                       Health Care and Social Assistance\n",
      "71                     Arts, Entertainment, and Recreation\n",
      "72                         Accommodation and Food Services\n",
      "81           Other Services (except Public Administration)\n",
      "99                                            Unclassified\n",
      "dtype: object\n"
     ]
    }
   ],
   "source": [
    "available_MSAs = [2,4,6,8,9,11,13,14,15,17,19,20,21,22,25,26,27]\n",
    "NAICS_codes_sig =  [ '--','11', '23','31-33','42','44-45','48-49','51','52','53','54','62','71','72' ]\n",
    "MSA_rank = 2\n",
    "MSA = msa_emp_rank_df_dict[ houseRank_to_empRank[MSA_rank] ]\n",
    "print_info(MSA_rank)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "NAICS_code = '--'\n",
    "ppc_relative_state, ppc_relative_country =  ppc_relative(MSA,NAICS_code)\n",
    "house_ppc_relative(MSA_rank)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.15228322,   2.2616642 ,   3.65815492,   5.7433837 ,\n",
       "        10.26725076,  16.63117116,  15.50893525,   5.44100332,\n",
       "        -3.28812701, -16.83779024, -17.69525293,  -0.61410168,\n",
       "         2.42790778,   1.71912517,  12.21029542,  11.56142221,\n",
       "         1.11534061,   0.92927852,   0.66086375])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#ppc_relative_country\n",
    "house_ppc_relative(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'available_MSAs' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-432123c00352>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0myear_range\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0myear_end\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0myear_start\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mlag\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0;32mfor\u001b[0m \u001b[0mMSA_rank\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mavailable_MSAs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m     \u001b[0mMSA\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmsa_emp_rank_df_dict\u001b[0m\u001b[0;34m[\u001b[0m \u001b[0mhouseRank_to_empRank\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mMSA_rank\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'available_MSAs' is not defined"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 864x432 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "NAICS_code = '--'\n",
    "\n",
    "fig = plt.figure()\n",
    "fsize = 6\n",
    "plt.gcf().set_size_inches(fsize*2,fsize)\n",
    "\n",
    "year_range = year_end-year_start-1\n",
    "lag = 0\n",
    "for MSA_rank in available_MSAs:\n",
    "    MSA = msa_emp_rank_df_dict[ houseRank_to_empRank[MSA_rank] ]\n",
    "    \n",
    "    ppc_relative_state, ppc_relative_country =  ppc_relative(MSA,NAICS_code)\n",
    "    x=ppc_relative_country.values[:year_range-lag]\n",
    "    \n",
    "    plt.subplot(121)\n",
    "    y=house_ppc_relative(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "    plt.scatter( x,y )\n",
    "    \n",
    "    plt.subplot(122)\n",
    "    y=relative_ppc_poverty(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "    plt.scatter( x,y )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Professional, Scientific, and Technical Services\n",
      "Lag = 0  (0.14424118617320927, 0.014804450158424352) (0.12372391949398245, 0.03683864374104193)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-32-2b54fcb76460>:23: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  plt.subplot(121)\n",
      "<ipython-input-32-2b54fcb76460>:28: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  plt.subplot(122)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lag = 1  (0.13711995235979796, 0.024239548305978995) (0.008257091511642169, 0.8925718044662494)\n",
      "Lag = 2  (0.1124532745930971, 0.07303652946133038) (0.08059023368830191, 0.19960160039815986)\n",
      "Lag = 3  (0.02978446299291954, 0.6461538879122724) (-0.025256515588397338, 0.697060455374101)\n",
      "Lag = 4  (-0.13920170105544827, 0.03692747152941683) (0.032871635323396875, 0.6238078786197612)\n",
      "Lag = 5  (-0.16191954829280553, 0.018877205834938886) (-0.15335526737438818, 0.0262680401507252)\n",
      "Lag = 6  (-0.09853104187504694, 0.17056143769630475) (-0.047401528509764826, 0.5105119992775395)\n",
      "Lag = 7  (-0.10249510919238307, 0.17095405944417122) (0.04175061853469658, 0.5778802439884522)\n",
      "Lag = 8  (-0.0537494577108009, 0.49292037416198814) (-0.18010967078522375, 0.020616954894588875)\n",
      "Lag = 9  (-0.017721446365892977, 0.8295772850140861) (-0.04735769105096309, 0.5649654480450154)\n",
      "Lag = 10  (0.05425242701626425, 0.532000525976768) (-0.13066162141226073, 0.13091153011974183)\n",
      "Lag = 11  (0.011591197082245331, 0.9000086856796956) (0.08372795146599063, 0.3632491456179836)\n",
      "Lag = 12  (-0.0020985434534231194, 0.9830491941469842) (-0.21951217206340065, 0.024453934956646283)\n",
      "Lag = 13  (0.05683486364990093, 0.5946731501375951) (0.18670289411791943, 0.07807122425618711)\n",
      "Lag = 14  (0.03958002175418302, 0.7360053326833913) (-0.09326768287503798, 0.42609458470207867)\n",
      "Lag = 15  (-0.033489272800769625, 0.7994792311647838) (0.0038445925807236614, 0.9767419871858236)\n",
      "Lag = 16  (0.10496970464962324, 0.4925601296818251) (0.2946201778905774, 0.04945991748896638)\n",
      "Lag = 17  (0.2743415579954408, 0.14234895149261637) (-0.26577855361474134, 0.15573799385851883)\n",
      "Lag = 18  (0.12964635910194067, 0.6451506836584735) (0.30452190890060155, 0.26978252672994346)\n"
     ]
    }
   ],
   "source": [
    "NAICS_code = '54'\n",
    "print( NAICS_dict[NAICS_code] )\n",
    "\n",
    "year_range = year_end-year_start-1\n",
    "for lag in range( int(year_range) ):\n",
    "    fig = plt.figure()\n",
    "    fsize = 6\n",
    "    plt.gcf().set_size_inches(fsize*2,fsize)\n",
    "\n",
    "    emp_relppc_s = []\n",
    "    house_relppc_s = []\n",
    "    poverty_relppc_s = []\n",
    "    \n",
    "    for MSA_rank in available_MSAs:\n",
    "        MSA = msa_emp_rank_df_dict[ houseRank_to_empRank[MSA_rank] ]\n",
    "        \n",
    "        \n",
    "        if NAICS_code in list(tot_empl_all[MSA].columns):\n",
    "            ppc_relative_state, ppc_relative_country =  ppc_relative(MSA,NAICS_code)\n",
    "            x=ppc_relative_country.values[:year_range-lag]\n",
    "            emp_relppc_s += list(x)\n",
    "\n",
    "            plt.subplot(121)\n",
    "            y1=house_ppc_relative(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "            plt.scatter( x,y1 )\n",
    "            house_relppc_s += list(y1)\n",
    "\n",
    "            plt.subplot(122)\n",
    "            y2=relative_ppc_poverty(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "            plt.scatter( x,y2 )\n",
    "            poverty_relppc_s += list(y2)\n",
    "            \n",
    "    print('Lag = {} '.format(lag), pearsonr(emp_relppc_s,house_relppc_s),pearsonr(emp_relppc_s,poverty_relppc_s) )\n",
    "    \n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Professional, Scientific, and Technical Services\n"
     ]
    }
   ],
   "source": [
    "NAICS_code = '54'\n",
    "print( NAICS_dict[NAICS_code] )\n",
    "\n",
    "max_lag = 12# int(year_range)\n",
    "\n",
    "#corr_auto = np.zeros(max_lag)\n",
    "#corr_house = np.zeros(max_lag)\n",
    "#corr_poverty = np.zeros(max_lag)\n",
    "#lags = list(range( max_lag ) ) # []\n",
    "lags = []\n",
    "corr_ = []\n",
    "corr_auto = []\n",
    "corr_house = []\n",
    "corr_poverty = []\n",
    "\n",
    "\n",
    "year_range = year_end-year_start-1\n",
    "for lag in range( max_lag ):\n",
    "    emp_relppc_s = []\n",
    "    emp_relppc_s_shifted = []\n",
    "    house_relppc_s = []\n",
    "    poverty_relppc_s = []\n",
    "    emp_us_ppc_s = []\n",
    "    \n",
    "    for MSA_rank in available_MSAs:\n",
    "        MSA = msa_emp_rank_df_dict[ houseRank_to_empRank[MSA_rank] ]\n",
    "        \n",
    "        \n",
    "        if NAICS_code in list(tot_empl_all[MSA].columns):\n",
    "            ppc_relative_state, ppc_relative_country =  ppc_relative(MSA,NAICS_code)\n",
    "            x=ppc_relative_country.values[:year_range-lag]\n",
    "            #x=ppc_relative_state.values[:year_range-lag]\n",
    "            emp_relppc_s += list(x)\n",
    "            \n",
    "            x_shifted = ppc_relative_country.values[lag:year_range]\n",
    "            #x_shifted = ppc_relative_state.values[lag:year_range]\n",
    "            emp_relppc_s_shifted += list(x_shifted)\n",
    "\n",
    "            y1=house_ppc_relative(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "            house_relppc_s += list(y1)\n",
    "\n",
    "            y2=relative_ppc_poverty(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "            poverty_relppc_s += list(y2)\n",
    "            \n",
    "            y3 = tot_empl_state_ppc[ ( 'United States' ,NAICS_code) ].values[lag:year_range]\n",
    "            emp_us_ppc_s += list(y3)\n",
    "            \n",
    "    #pcorr = pearsonr(emp_relppc_s,emp_relppc_s_shifted)\n",
    "    pcorr = pearsonr(emp_relppc_s,house_relppc_s)\n",
    "    #pcorr = pearsonr(emp_relppc_s,emp_us_ppc_s)\n",
    "    if len( emp_relppc_s )>40:\n",
    "        lags += [lag]\n",
    "        corr_ += [pcorr[0] ]\n",
    "    #corr_house += [pearsonr(emp_relppc_s,house_relppc_s)[0] ]\n",
    "    #corr_poverty[lag] = pearsonr(emp_relppc_s,poverty_relppc_s)[0]\n",
    "    #print('Lag = {} '.format(lag), corr_auto[lag], corr_house[lag], corr_poverty[lag])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Correlation')"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjwAAADUCAYAAACYh/lcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeyUlEQVR4nO3df3hdVZ3v8feHUDSCGJACbQDLaI3UW6EYQMRHQYUA3msrjg6OI+A4VK/gb+Nt0VHu+CgdqyM4DwIF+TUiyPWW0qvVyA8ZRwehxTi0gJFaoCQpNFgCitGW9nv/2Dt4cjhJzknOyd7n5PN6nvOcvddea5/vOUD7Za2111JEYGZmZtbIdss6ADMzM7Nac8JjZmZmDc8Jj5mZmTU8JzxmZmbW8JzwmJmZWcNzwmNmZmYNb/esA8jSfvvtF3PmzMk6DDMzM6uCe+6554mImFnqWq4SHkknAxcBTcAVEbGs6PqrgKuAI4HPRsRXC649DPwe2Ak8GxHt433enDlzWLduXfW+gJmZmWVG0iOjXctNwiOpCbgYOBHoBdZKWh0R9xdU2wZ8FFg0ym1OiIgnahqomZmZ1Z08zeE5GtgYEZsiYjtwA7CwsEJEbI2ItcCOLAI0MzOz+pSbHh6gFXi04LwXOKaC9gH8WFIAl0XEilKVJC0GFgMccsghEwx1dKu6+1je1UP/4BCzW5rp7Ghj0YLWqn+OmZmZlS9PCY9KlFWy0ddxEdEvaX/gFkm/joifPu+GSSK0AqC9vb2qG4mt6u5j6cr1DO3YCUDf4BBLV64HcNJjZmaWoTwNafUCBxecHwT0l9s4IvrT963ATSRDZFNqeVfPc8nOsKEdO1ne1TPVoZiZmVmBPCU8a4G5kg6VtAdwOrC6nIaS9pT04uFj4CRgQ80iHUX/4FBF5WZmZjY1cjOkFRHPSjoX6CJ5LP3KiLhP0ofS65dKOhBYB+wN7JL0cWAesB9wkyRIvtN3IuJHU/0dZrc001ciuZnd0jzVoYzJ84zMzGy6yU3CAxARa4A1RWWXFhw/RjLUVexp4PDaRje+zo62EXN4AJpnNNHZ0ZZhVCN5npGZmU1HeRrSqnuLFrRywWnz2aMp+VlbW5q54LT5uUokPM/IzMymo1z18DSCRQtauf7uzQB894PHZhzN83mekZmZTUfu4ZlmRptPlLd5RmZmZtXkhGea6exoo3lG04iyvM0zMjMzqzYPaU0zw/OJPvO9e9m+cxetfkrLzMymASc801De5xmZmZlVm4e0zMzMrOE54TEzM7OG54THzMzMGp4THjMzM2t4nrRsueU9v8zMrFqc8Fguec8vMzOrJg9pWS55zy8zM6smJzyWS97zy8zMqskJj+WS9/wyM7NqylXCI+lkST2SNkpaUuL6qyTdKenPkj5dSVurL97zy8zMqik3k5YlNQEXAycCvcBaSasj4v6CatuAjwKLJtDW6oj3/DIzs2rKTcIDHA1sjIhNAJJuABYCzyUtEbEV2CrpbZW2tfrjPb/MzKxa8jSk1Qo8WnDem5bVuq2ZmZk1uDwlPCpRFtVuK2mxpHWS1g0MDJQdnJmZmdWvPCU8vcDBBecHAf3VbhsRKyKiPSLaZ86cOaFAzczMrL7kKeFZC8yVdKikPYDTgdVT0NbMzMwaXG4mLUfEs5LOBbqAJuDKiLhP0ofS65dKOhBYB+wN7JL0cWBeRDxdqm0mX8TMzMxyJzcJD0BErAHWFJVdWnD8GMlwVVltzczMzCBfQ1pmZmZmNeGEx8zMzBqeEx4zMzNreE54zMzMrOE54TEzM7OG54THzMzMGp4THjMzM2t4TnjMzMys4TnhMTMzs4bnhMfMzMwanhMeMzMza3hlJzySXirpKkl9kp6VtLPwVcsgzczMzCajks1DvwUcAXwD6AOiFgGZmZmZVVslCc8JQEdE/KJWwZiZmZnVQiUJz5PA07UKxMzMzCZmVXcfy7t66B8cYnZLM50dbSxa0Jp1WLlSyaTlLwGflVRJklQRSSdL6pG0UdKSEtcl6Rvp9XslHVlw7WFJ6yX9StK6WsVoZmaWJ6u6+1i6cj19g0ME0Dc4xNKV61nV3Zd1aLlSSfLybuAooE/SA8D2wosRcdJkApHUBFwMnAj0AmslrY6I+wuqnQLMTV/HAJek78NOiIgnJhOHmZlZPVne1cPQjpHPDg3t2Mnyrh738hSoJOHpTV+1cjSwMSI2AUi6AVgIFCY8C4FrIyKAX0hqkTQrIrbUMC4zM7Pc6h8cqqh8uio74YmI99cyEKAVeLTgvJeRvTej1WkFtpA8NfZjSQFcFhErahirmZlZLsxuaaavRHIzu6U5g2jyq+KFByUdLKlD0kmSDqpiLCpRVvzo+1h1jouII0mGvc6R9MaSHyItlrRO0rqBgYGJR2tmZpYDnR1tNM9oGlHWPKOJzo62jCLKp0oWHnyRpGuAh4EfAj8CHpZ0taRqpJG9wMEF5wcB/eXWiYjh963ATSRDZM8TESsioj0i2mfOnFmFsM3MzLKzaEErF5w2nz2akr/SW1uaueC0+Z6/U6SSHp7lwPHAO4B90tc7SdbnWV6FWNYCcyUdKmkP4HRgdVGd1cAZ6dNarwOeiogtkvaU9GIASXsCJwEbqhCTmZlZ7i1a0MqCQ1o45tB9+fmSNzvZKaGSSct/DZwREV0FZTdL+jNwDXDuZAKJiGclnQt0AU3AlRFxn6QPpdcvBdYApwIbgT8Cw/OKDgBukgTJd/pORPxoMvGYmZlZ46gk4XkJyXBWsYeAvasRTESsIUlqCssuLTgO4JwS7TYBh1cjBjMzM2s8lQxpbQDOLlG+GA8fmZmZWY5V0sPzeZIhrDcAPyV5OupNwJHA22sQm5mZmVlVVLIOzxpJrwU+DbwlLb4PODsi1tciOLO88/41Zmb1oaJ9sSLiXuCMGsViVleG968ZXtJ9eP8awEmPmVnOVLzwoJklxtq/xszM8mXMHh5J24HWiBiQtIPnr3z8nIjYo9rBmeWZ968xM6sf4w1pnQ08XXA8asJjNt14/xozs/oxZsITEdcUHF9d82jM6khnR9uIOTzg/WvMzPKqkr20Nkl6aYnyFkmbqhuWWf55/xozs/pRyVNac0i2fCj2AsB/wtu0tGhBK9ffvRmA737w2IyjMTOz0Yyb8Eh6Y8HpsZKeLDhvItmos7fagZmZmVnjyHrdsnJ6eO4gmawcwE0lrv8B+HAVYzIzM7MGkod1y8pJeA4GBGwm2UZioODaduCJdFNPMzMzs+cZa92y3CQ8EdGXHnqRQjMzM6tYHtYtq2hrCUm7A0cBLwNGLDQYEddWMS4zMzNrEHlYt6ySx9LnAvcDPweuA64CrgauBK6oRjCSTpbUI2mjpCUlrkvSN9Lr90o6sty2ZmZmlo3OjjaaZ4x80Huq1y2rZJjq6yQJz37AH4F5wBuAX/KX3dMnTFITcDFwSnrv90iaV1TtFGBu+loMXFJBWzMzM8tAHtYtq2RI6xjgrRGxTVIARMR/SloKXAi8dpKxHA1sjIhNAJJuABaSJFnDFgLXppOkf5EuejiLZI2g8dqamZlZRrJet6yShGcGf9lX6wngQKAHeAg4rAqxtAKPFpz3kiRZ49VpLbPtlDn5ju9w4MCjPPKzvbMKYVxnbUn+UeY5RqiPOOshRjNrfPXwZ9FZW57msZkHQ84Tnl+TJDYPAb8CPiLpMeCjjEw2Jkolyoofdx+tTjltkxtIi0mGwzjkkEMqia9sp8yfxZ8feKom966WebPy+x9EoXqIsx5iBLg//cMwz/HWQ4xQH3HWQ4xQH3HWQ4yQ//ggiXHBYbMy+exKEp6LgJnp8T8BXSRDRn8G/q4KsfSSrPkz7CCgv8w6e5TRFoCIWAGsAGhvb6/J+kEHnndeLW5rNimfuexOIN9bYNRDjFAfcdZDjFAfcdZDjDa+shOeiLi+4PhXkuaQ9Pg8EhG/q0Isa4G5kg4F+oDTgb8tqrMaODedo3MM8FREbJE0UEZbMzMzm6YqWoenUEQMkTyhVRUR8aykc0l6jpqAKyPiPkkfSq9fCqwBTgU2kjwp9v6x2lYrNjMzM6tvYyY8ksoem4mIL082mIhYQ5LUFJZdWnAcwDnltjUzMzOD8Xt4zi7zPgFMOuExMzMzq4UxE56IOHSqAjEzM8ubVd19dG8eZPvOXRy37HY6O9qmdLE8qx5vCGpmZlbCqu4+lq5cz/aduwDoGxxi6cr1rOruG6el5VFFCY+k90vqlvR0+kQUkj4j6Z21Cc/MzCwby7t6GNqxc0TZ0I6dLO/qySgim4xKNg9dDHwNWEmy6vLwYn8DwLnVD83MzCw7/SV29x6r3PKtkh6ejwAfjIgvAs8WlN8DvLqqUZmZmWVsdktzReWWb5UkPK8A7i5R/gyQ//WszczMKtDZ0UbzjKYRZc0zmujsaMsoIpuMShYe3EKS9DxSVH4ssKlqEZmZmeXA8NNYy7t66B8cYnZLs5/SqmOVJDzXAl+T9G6SdXeaJZ0K/DPw9VoEZ2ZmlqVFC1qd4DSIShKeLwJzgAdIJizfm5ZfRTKZ2czMzCyXykp4JO0OnAR8GjgfaCeZ/3NPRPy2ZtGZmZmZVUFZCU+6OedK4FUR8TDwcC2DMjMzM6umSp7Sup9kSMvMzMysrlSS8HwS+GdJr5e0R60CMjMzM6u2SiYt30KSIP0HgKQR621HhJMgMzMzy6VKEp5/qFkUZmZmZjVU7lNaM4BTgc9GxMZqByFpX+C7JHOEHgbeHRFPlqh3MnAR0ARcERHL0vLzgbNJ9vUCOC8i1lQ7TjMzM6tPZc3hiYgdQAewc7y6E7QEuC0i5gK3pecjSGoCLgZOAeYB75E0r6DK1yPiiPTlZMfMzMyeU8mk5R+Q9PLUwkLgmvT4GmBRiTpHAxsjYlNEbAduSNuZmVkdWtXdR/fmQe56aBvHLbudVd19WYdkDaySOTy/AM6XdDiwlmTT0OdExHcmEccBEbElvc8WSfuXqNMKPFpw3gscU3B+rqQzgHXAp0oNiQFIWgwsBjjkkEMmEbKZmU3Uqu4+lq5cz/aduwDoGxxi6cr1AN7KwWqikh6ei4CXkkxevgz4dsHr38ZrLOlWSRtKvMrtpVGJskjfLwFeDhxBssnpqFtdRMSKiGiPiPaZM2eW+dFmZlZNy7t6GNoxcpbE0I6dLO/qySgia3Rl9/BERCXJUan2bx3tmqTHJc1Ke3dmAVtLVOsFDi44PwjoT+/9eMG9Lge+P5lYzcystvoHhyoqN5usSSUxVbQaODM9PhO4uUSdtcBcSYemCx+enrYjTZKGvQPYUMNYzcxskma3NFdUbjZZFSU8kk5Ih6a2SOqXdIuk46sQxzLgREkPAiem50iaLWkNJPt5AecCXSQ7tt8YEfel7b8iab2ke4ETgE9UISYzM6uRzo42mmc0jShrntFEZ0dbRhFZoyt7SEvSe0jm69xMkpAIeBNwq6T3RsR3JxpERPwOeEuJ8n4KngxLHzd/3iPnEfG+iX62mZlNveGJycu7eugfHGJ2SzOdHW2esGw1U8lTWp8DPhcRFxSUXSjpPOAfSRYONDMzK8uiBa1OcGzKVDKk9Qrg/5QovzG9ZmZmZpZLlSQ8A8BrSpQfwV+2dDAzMzPLnUqGtL4NXCZpJsmO6UEyh+eLwOU1iM3MzMysKiqdw9NEsgDhDJJJy38GvgF8vvqhmZmZmVXHuENaknaT9BpgRkR0AvuQDGMdTrLy8nXArloGaWZmZjYZ5czheS9wLbAdICKGImJ9RKxPy64F/qZ2IZqZmZlNTjkJzweAr0XEzuIL6WKAXyXdjNPMzMwsj8pJeF4F/OcY1+8EDqtOOGZmZmbVV07C8xKSScqj2QPYuzrhmJmZmVVfOQnPIySTlEdzBLC5GsGYmZmZ1UI5Cc9q4IuS9iq+IGlv4H8D/6/agZmZmZlVSznr8CwjeQrrN5L+lWSn8gBeTbJ7+Y60jpmZmVkujZvwRMQ2Sa8HLiVZVXm4V2gX8EPgw+lu52ZmZma5VNZKyxHRB/wPSfuQbBQq4MGIeLKWwZmZmZlVQyWbhxIRT0bE2oi4u5rJjqR9Jd0i6cH0fZ9R6l0paaukDRNpb2Y2Wau6++jePMhdD23juGW3s6q7L+uQzKwMFSU8NbQEuC0i5gK3peelXA2cPIn2ZmYTtqq7j6Ur17N9Z7KbTt/gEEtXrnfSY1YH8pLwLASuSY+vARaVqhQRPwW2TbS9mdlkLO/qYWjHyEXnh3bsZHlXT0YRleZeKLPny0vCc0BEbAFI3/evVXtJiyWtk7RuYGBgwgGb2fTTPzhUUXkW3AtlVtqUJTySbpW0ocRr4VTFABARKyKiPSLaZ86cOZUfbWZ1bnZLc0XlWaiXXiizqVbWU1rVEBFvHe2apMclzYqILZJmAVsrvP1k25uZjauzo42lK9ePSCiaZzTR2dGWYVQj1UMvlFkW8jKktRo4Mz0+E7h5itubmY1r0YJWLjhtPq0tzQhobWnmgtPms2hBa9ahPaceeqHMsjBlPTzjWAbcKOkDJPtyvQtA0mzgiog4NT2/Hjge2E9SL/CFiPjWaO3NzKpt0YLWXCU4xeqhF8osC7lIeNKVmt9SorwfOLXg/D2VtDczm26Gk7HlXT30Dw4xu6WZzo62XCdpZlMhFwmPmZlVT957ocyykJc5PGZmZmY144THzMzMGp4THjMzM2t4TnjMzMys4TnhMTMzs4bnhMfMzMwanhMeMzMza3hOeMzMzKzhOeExMzOzhueEx2waWNXdR/fmQe56aBvHLbudVd19WYdkZjalnPCYNbhV3X0sXbme7Tt3AdA3OMTSleud9JjZtOKEx6zBLe/qGbFzNsDQjp0s7+rJKCIzs6nnhMeswfUPDlVUbmbWiJzwmDW42S3NFZWbmTWiXCQ8kvaVdIukB9P3fUapd6WkrZI2FJWfL6lP0q/S16lTE7lZ/nV2tNE8o2lEWfOMJjo72jKKyMxs6uUi4QGWALdFxFzgtvS8lKuBk0e59vWIOCJ9ralBjGZ1adGCVi44bT6tLc0IaG1p5oLT5rNoQWvWoZmZTZndsw4gtRA4Pj2+BrgD+F/FlSLip5LmTFlUZg1i0YJWJzhmNq3lpYfngIjYApC+7z+Be5wr6d502KvkkJiZmZlNT1OW8Ei6VdKGEq+FVbj9JcDLgSOALcDXxohjsaR1ktYNDAxU4aPNzMws76ZsSCsi3jraNUmPS5oVEVskzQK2VnjvxwvudTnw/THqrgBWALS3t0cln2NmZmb1KS9DWquBM9PjM4GbK2mcJknD3gFsGK2umeWTt78ws1rKS8KzDDhR0oPAiek5kmZLeu6JK0nXA3cCbZJ6JX0gvfQVSesl3QucAHxiasM3s8nw9hdmVmuKmL6jOu3t7bFu3bqswzCb9o5bdjt9JVZ+bm1p5udL3pxBRGZWjyTdExHtpa7lpYfHzKYxb39hZrXmhMfMMuftL8ys1pzwmFnmvP2FmdVaXlZaNrNpbHgV6OVdPfQPDjG7pZnOjjavDm1mVeOEx8xywdtfmFkteUjLzMzMGp4THjMzM2t4TnjMzMys4U3rhQclDQCP1Oj2+wFP1Oje041/y+rxb1k9/i2rx79l9Uz33/JlETGz1IVpnfDUkqR1o632aJXxb1k9/i2rx79l9fi3rB7/lqPzkJaZmZk1PCc8ZmZm1vCc8NTOiqwDaCD+LavHv2X1+LesHv+W1ePfchSew2NmZmYNzz08ZmZm1vCc8FSZpJMl9UjaKGlJ1vHUK0kHS/qJpAck3SfpY1nHVO8kNUnqlvT9rGOpZ5JaJH1P0q/Tfz+PzTqmeiXpE+l/3xskXS/phVnHVE8kXSlpq6QNBWX7SrpF0oPp+z5ZxpgnTniqSFITcDFwCjAPeI+kedlGVbeeBT4VEYcBrwPO8W85aR8DHsg6iAZwEfCjiHgVcDj+TSdEUivwUaA9Iv4b0AScnm1Udedq4OSisiXAbRExF7gtPTec8FTb0cDGiNgUEduBG4CFGcdUlyJiS0T8Mj3+PclfKt5ZcoIkHQS8Dbgi61jqmaS9gTcC3wKIiO0RMZhpUPVtd6BZ0u7Ai4D+jOOpKxHxU2BbUfFC4Jr0+Bpg0VTGlGdOeKqrFXi04LwX/yU9aZLmAAuAuzIOpZ5dCHwG2JVxHPXur4AB4Kp0ePAKSXtmHVQ9iog+4KvAZmAL8FRE/DjbqBrCARGxBZL/cQT2zzie3HDCU10qUebH4CZB0l7A/wU+HhFPZx1PPZL034GtEXFP1rE0gN2BI4FLImIB8AweMpiQdG7JQuBQYDawp6S/yzYqa2ROeKqrFzi44Pwg3EU7YZJmkCQ710XEyqzjqWPHAW+X9DDJMOubJX0725DqVi/QGxHDvY3fI0mArHJvBR6KiIGI2AGsBF6fcUyN4HFJswDS960Zx5MbTniqay0wV9KhkvYgmYC3OuOY6pIkkcyTeCAi/iXreOpZRCyNiIMiYg7Jv5O3R4T/T3oCIuIx4FFJbWnRW4D7Mwypnm0GXifpRel/72/BE8CrYTVwZnp8JnBzhrHkyu5ZB9BIIuJZSecCXSRPHFwZEfdlHFa9Og54H7Be0q/SsvMiYk12IZkB8BHguvR/ajYB7884nroUEXdJ+h7wS5KnMrvxKsEVkXQ9cDywn6Re4AvAMuBGSR8gSSrflV2E+eKVls3MzKzheUjLzMzMGp4THjMzM2t4TnjMzMys4TnhMTMzs4bnhMfMzMwanhMeMxuVpOMlRboXVxaff4ekKdv/S9Kn6mU3eUlNku5PV9I2s3F4HR6zaUrSeGtSPAK8EpgWq7VK2hf4R+CErGMpR0TslHQ+8DVJayLC+6SZjcE9PGbT16yC18K07OiCsqPS3cAfmyZ/mX4A+E1EdGcdiKTdJDWVUfUm4KXAqTUOyazuOeExm6bSROaxdLuEbWnxQEH5QPGQVsH5qZLulDQk6R5Jr05fP5P0R0l3S5pX+HmSXivpx5L+IGlA0kpJL6skZkknpsNc2yQ9JenfJR1dVOfQ9HP+JGmzpHPKHBp7L7Cq4D4nSNopqXB/PCSdKen3kl6cnh8g6er0O/1e0s8lvbGgviRdLum36e+1SdKXJb2goM75kjZK+htJvwa2A4elv2mXpEFJz0h6QNL7htule1D9APBWIWbjcMJjZhPxJeCzwGtJ/nK+HriEZGn74bKrhiunyc+/A3cC7cCbgZ3ALZJeWMHn7gVcDLyOZKPJB4EfSXpp+jki6fV4CfBG4O3A24AFY9003bn7NcDdw2UR8ZP0/n9fVP0fgBsi4veSmoGfAC8GTkk/Z036vQ4bvj3wOPC3wGHAx0m2oziv6L6zgQ8DZwHzSIYUrwd+l37X+cAngSeL2t1FnQzDmWUqIvzyy69p/gLeAAQwp6j8+LT8oKLzRQV13pWWvbOg7B1p2V7p+dUkSULhvV8A/LHwXiXiugO4Yozru5EkAO9Nz09MP/cVBXX2TT9nrPsckbY7rKj8kySJx27peVta76j0/CySHdR3L2p3O3DhGJ/3CeDBgvPzgV3AIUX1ngLOGuef3dvTmPbM+t8jv/zK88uTls1sIv6r4Pix9P3eEmX7A38AjgJeIekPRfd5ITC33A+VdCjwT8Cx6b13A14EDA+NzQOeiIiNw20iYpuknnFu3Zy+/6mo/GqS3qwO4IfA2cB/RcTa9PpRwIHAYNK59JwXAEMFcZ9N0jM0B9iT5IGR4h72xyNic1HZV4ErJJ1FkvytjohfFtUZjrkZeGb0r2g2vTnhMbOJ2FFwHGOU7Vbw/m8kOzkX+10Fn/t94AngHOBRkqGznwF7lPjsSgyk7/sCDz13oyRZ+h5wtqRbgTNIemOG7QY8QNKjVeyPAJLeRTIMt4RkWO9pkl6xLxXVf16yEhFflHQdcDLJMOB5kr4SEZ8rqLYvyfDgtuL2ZvYXTnjMbCqsI5kj89uImEhCQjpPZx5wakR0pWUHkfT0DLsfmCnpFcO9POn8nFcC94xx+03AIPDqEvUuI5mn8yGS3pnrir7XGcDTETHao/tvBLoj4l8KvsucMWIZISI2Ad8EvilpCdAJFCY889P7T4cn6cwmzJOWzWwqfJlkwu63JR2dPkl1gqSLJP1Vmfd4kqQn5mxJr5R0LMmk3qGCOreSDLddK+koSYeT9Cw9yxg9P2my0AW8qcS1nwE9JMNLN0bEUwWXryPpEfqBpJMkzZF0jKSlkhaldXqA+ZIWSnq5pI8Bp433ZSXtJeliSW9Of68FJD099xdVPZ7kSS0zG4MTHjOruYh4gORJo71IEov7gctJ5p0MlnmPXSRDQS8nmS90NXAhsKWgTpAMLz0D/AfJENgPSZKO4vk5xS4B/jp98qrY5STDZiuKYvoTSZK0juSptN8AK0nWM3okrXYZSdJ1FdANHMPIYbHRPAvsA3yLZNisi7887QVAmiwendYxszFogr3LZmZ1IV0vpxf4XET86zh1bwW+HxEXFpV/BTglIubXLNAJkPRNkj/H/2fWsZjlnefwmFlDkfR2kt6RB0jm93yBZDjrxjKaf5jk0fbhe72EZI7M2SSPkueGpN1IErkV49U1M/fwmFmDkXQ68HmSR8CfIZmE/OmI2DCBe91BMgT1XeDvPTHYrH454TEzM7OG50nLZmZm1vCc8JiZmVnDc8JjZmZmDc8Jj5mZmTU8JzxmZmbW8JzwmJmZWcP7/1p0IHWzIA8CAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 648x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "fsize = 3\n",
    "plt.gcf().set_size_inches(fsize*3,fsize)\n",
    "\n",
    "plt.stem( lags, corr_ )\n",
    "#plt.stem( lags, corr_house )\n",
    "plt.xlabel( 'Time lag (years)' ,fontsize='x-large')\n",
    "plt.ylabel( 'Correlation' ,fontsize='x-large')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAHgCAYAAABtixjUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABi5klEQVR4nO3de7xUdb3/8ddbQCOtyPCCKGJFKmaF4S0qLTPUPGF2s4upp7KbZafSH1qn7Conq1OdTCMzrcxLHkNOWmSamlYqiFeQJES5qaihoggCn98f3+/oYpw9e2/27Fl7Zt7Px2Mee9Z3rVnrs9asmb0+870sRQRmZmZmZmbtZJOyAzAzMzMzM2s0JzpmZmZmZtZ2nOiYmZmZmVnbcaJjZmZmZmZtx4mOmZmZmZm1HSc6ZmZmZmbWdpzomLUgSftLCkn7N3i9o/N6j27kevuLpIMlzZa0Ksc9us6yIemU5kXXNUlHdxfvQNcO+9AqJF0t6eoeLlv5bvhIH7e3eGNf3wyt9l1Vj6RTJPleH2b9wImOdTxJwyR9JV8wPybpKUnzJf1U0riy4+sPkj7X6hcIkl4CXAQI+DRwJLC81KCsY0jaMl+g7l/CtjfJ2z6s2ds2M2slg8sOwKxMknYDfg+MAH4D/Ax4ChgDvBv4sKRRETGgf93cCJ8D5gPnVJXfCwwFnm52QBthT2AL4BsRcXHZwVjH2RL4Sn5+dT9v661V05vkbZ8LTOvnbQ9UrfRd1Z1vAFPKDsKsHTnRsY4laQtgOrA5sHdE3Fw1/4vACaQag0Zs7/kR8WRv5zVTRAQp0WsFW+e/K8oMwqy/RcSasmMYaFrsu6qmyvd+RKwF1pYdj1k7ctM162THAi8FvlCd5ABExNqIODUiFlXKJI2UdI6kByStljRH0n9I2iAZqrRxl7SzpMslPQZc3t28PH8/SX+U9Gjue/J3SYd2tzOSdpf0M0l3S3pS0gpJf5C0Z9VyAYwE9stt3EPSwjyvZrv3jdjvUZKmSXpc0iOSzpS0WXf7kNcxTNIP83pWS/qnpK8VX5/jPTdPXpFjvron6+/ttgrL7ijpbElL8rKLJJ0naWRhmc9JukbSg3mZ+ZK+KmlIb2Or2vZHJd2cz4cVki6VNLZqmVPycXiVpB/lGB6TdJGkF0sapNREc5FS88wrJe3UxTrGSZoq6SFJT0j6naSX9jDWffL5+1h+7fWSDi7MH5SP9x+7eH3lHNqkanqnHMfj+Tz8upKXSPplPs8ey5+B59VY76sl/VbSw3n/b5N0TNUylfP/G5LeKen2vOw/JL2nsNz+wN158it69nN0Tp3jcqGkpVVlX8qvO6Oq/O+S/lJ1TK6uxMiztRhHFbZ9dY1tfjyfg6sl3SLpTV3F10XMPfocK/WVuz6/349JmiFp76pluuxXKGlh9bGT9DFJt0paqfRdeLukrxbmP+e7Ss/2G5uYj23lXL9e0qtrbPcVkn6f416e9++V1eutc3wWSrpO0t55G08qfT98XdLgqmXr/U+o2UdH6Tv9IqXP8lOSFih9Ll9QtVxPvh+2lvQTSffm8+EBSVdJenN3+2nWylyjY53sHcBq4Nc9WVipT8hfgW2B04EFwKHA94CXAcdVveT5wJ+AGcAXgHXdzZP0TuBC4G/AV0m/8r0PmC7pfRFxYZ0QJwKvzq+/F9gG+DBwjaQ9IuKuvNyRwA+BB4Bv5rKVDdzvocCVwDWkGrF9gY+R+s/8Z534yRdRVwLjSM0IZwNvzK8bB/xbXvSzwMGkZHUKcGfenx7rxbaQtDNwPel9+ykwh1Sj9Dbg5cCSvOgJpPd0OvAk8HrgS8AoYIOL6l7E+X1SH6QLgLOAFwGfAv4qaXxEzK96yc9Jx+KrwCuBjwObAYvyfp0GbAd8HjgPeF2Nzf4ceAz4Gikp/jRwraRXRcQjdWKdQDqmD+XtrCLt92WSjoiIiyJinaRfAP9P0vbFZqGSdiS9B/8VEesLqx5K+rz8ETiR9Nn9Eum8fS8wF/gisB/w73n/Ty6s93XAFcA/gW8Dj5Pe37MlDY+I06p25aAc95mkGsNjgfMl3RIR/8jb+wLwHeC3wCX5df/s6tiQmre9R9LOETEvl+0HrAeeSUAkbQ68FvivLtazHDiKlOj/BZiay6vP/48BLyCdr0+TPjOXStoxIv5VJ86KHn2OcwJ4ATAPOAXYlHTOXSPpgIi4vgfb2oBSAnomcGn+C7Az6Xj1xDeAIH1HDSW9V9Mkjcm1J0jaGriWdIx+CCwFDgd+0ctwtyM1f/51fhxMOjeHA5+oWrbe/4QNSHpDXm4V6T1cAGyfY3wJ6RzuzffDb0if/9NJzZa3BPYG9gCu6uU+m7WOiPDDj458AA8Dt/Zi+W+T/nm+s1Am0kVOALsXyq/OZZNrrKfmPNI/wYeAi6vKBwE3AfcBymX753XsX1hu8xrbGg48CJxZVb4YuLrG8qPzeo/u435/tmq904EHe3CMP5lf//mq8v/O5YcWyo7OZW/p4fsXwCkbua0/AWuA19RYr7p5D75CuqAZWSP20d3EvHde7riq8u1JicivCmWn5GUvqVr2N7n878CgQvlpNd6/yjquBwYXyv8tl/9XvX0AbgSeAEYVyl5ESryXAUNy2Rhqfwa+lMt3rnFOHV8oG0K6MF0PfLdqHTcBD1edq3eSfjwYXLXsJTneF1Wd/yuB7QvLjSD9KPLtQtnLq8+pbt7LXfPyxxb2YSVwfi7fNpe/NU8fWHUMri5MD87LnFNjO/vneYuALQrle+TyT/Yg1soxr/s5znEszdsaVnV+Pg7MrBHX/jW2t7C4L6Tk8c5uYqy8V8Xvqso5eRMbnr+H5/JDCmXfzWVvLZQNKuz70fW2X4g7gE9UlV+cy8fWOKa1/iecQm6Nl6c3Af4BPALsUGP5yv+BHn0/kD6DAZzQk3PVDz/a6eGma9bJXkj6Z9BTbwfmR8T/VgoiIkgXjFCoASg4o0ZZV/PeQvql7leShlcewItJvxjuALyiq5VFxBOV55Ken2tiIF187ln7VT3S2/1eD/ykquzPwFbVTS662NYTpF8di75dmN8oPdpWfg/eDFwYEbdUryQfi8rzJ/JrBik1ixtO2vdNSBeavfVeUoL126pz4inSxdwBNV5zZtV05Rf1n0XEuhrlL6+xjh9F/uUbICL+j9RUq9Y5DoCkbUjn2fkRcV/htY/mmLbN84mIu4HrSDUTRR8C/h7P1nhUbHBORcTTpPNa1N7fLSW9OE/vDowl/do+rOo4Xk76gWGfqnVMj0JNU0QsA+4i1WBulIiYS/rRYf9cNJ7UP/BbpBqXSvl+pJrcv27strJfRsQzNbWRmuc+Rs/3oSef4/GkJPAnEbGisK3FpOP9WknbbUTsjwIjJVW/Lz310+L5S4obNtz3twF3RcQzTSjz5+NHvdzW46Qa4aLv5b+1mhzX+59Q8RrSjwGnR6HpdEXhO6en3w+r8nL7F/4vmHUEJzrWyR4jJTs9NZp0sVNtTv67U1X5I/kir5Za83bJf39Lah5SfFSaimxNFyS9UKm/yTLSBfxD+bVvIyVLG2s0vdvvByNiVVXZv/LfLXuwrXsiYoNOxvlCc0WNbfVFT7f1MtIF9W3drVDSQZL+Srqw+Bfp+F+TZ2/Me7ALqSnQYp57TryZ2ufDvVXTK/Lf+7oor/WeVCcalbJ6x78yb26NebXOlZ8Du0jaC0DSvqSLu3NqvP7B6veJnu9X5XP1Q557DH+a51Ufx+pjCOn97O787c41PNv8aj9gSUTcTroo3b9QPrP4w8VG6us+9ORzPDr/7el73lNTSN/Pf8v9Us6W9HZpwz6BdWyw7/FsU73ivo/m2X5WRf/oZawL47mDRVQ+P735n1A0Jv/t7junR98POb4TSbWF90v6q1JfxF1qrtWsjbiPjnWyOcCekp5X4yKqK9GLedUXCd3Nq/zw8AlSG+pa7qizzvNJv+D9N3Az6VfR9cBJ9OGX6Kw3+91lu3N6NoJdV9tSN3FsjN5sq+6286/PlwE3kNrMLyY1dxpJunjfmB+WNiH19ZnUi9d0dfy7Kq/1nvTlONd6rWrMu4iUfBxFqp35EOmX6Fr90OqdU93tV+W4n8KztVjV7uzlOjfWNcC7JY0hJTRX5/KrgXdKGkqq9fpezVf3Tl/3oTef45685/XOqUEbrCzirtwvbiLp4nwiqc/UHyUdUlUzWUtP971e3D3Vm89Kvf8JG7PuHn8/RMQPJP2WVFN9APAfwEmSjo2In/cyLrOW4UTHOtk0Umfx9wNn92D5haR29tV2Lczvi0py80hE/Kk3L5Q0DDgE+GpEnFI17+s1XtKbf84L6d/9rt7WhOrkU9K2pHbmZWxrPul4vaqb9VWakRxQ/CVc0sQ+xDifdJF3e0T0arCFPtqFNDhD0c7APXVeszD/HVtj3i5VyxARKyX9BjhC0mTS8ZtWbALVIJXP1arefq66sTHJ4NX57wHABNKAEJASoJNJfUk25dlawEZuuz8szH/H8uyADBXV73mlVmWDmk2lEfJGVK84f4amkQYREHAq8P9Ix67miH29tJDaTYHH1CirZydJm1bV6uyc/9b7vNRTqWl6Fam/T1d69f2Qm5T+CPhRbtr5N9KANE50rG256Zp1sqmkf3anSRpXPVPSYEn/T9L2uej/gJdLekdhGZFGz6nM74sZpIuBk/Mvu9XxdNlsjfQLZlD1mVYaynXvGsuvpOdNqfp7v6u3tQXPHa3ohLK2FREPk0YlOkLSc5KdQnOaynswqDBvUGF9G+P8/PdrtWZK2qoP667nOBWGx5X0b6QLwN919YKIuB+YSTpO2xde+wLSaF33k5poFZ1Nak40lXQ+ntOg+ItuJjUlOj73X9hAH45hpf9Lj5skRsSdpGZFx5NG+6okNNeT+uV8iXQe1R2pLNdoPNWbbfeTmaRBJo6V9Ewz4Nwv5wPAzRFRGVJ7IWkfq4cz/gxVNTrV/Uhyn5Rb8mRfmw9WXE5qOvnMzVjz5/VTvVzPC0ijWxZ9Lv/t8vPSjVtIyc6nVBi+vqLwndOj74fcZ3OD/ym5Od9Cyj+HzPqVa3SsY0XE45LeTurof0P+dfmvpAuIlwPvIt1n51f5JVOA95CGma0Ms/w20nCip+e29n2N5yOkpjt3SvolaTSjEaTO0rvQRRO0/NorgRPzP7R/kDphH0NqllM9CMBNwAclfSUvuzJ3OK+lX/e7ylmk4YG/m9uP3wK8gTTE9u/Y+AuHvm7r06SLz79JqgwvPZxUi/ZF0gXrpaQLnKsknUsa1va99OEHpYi4XtL3gM8p3Rfj/0h9UHYkHf/bSCNNNdrmpP34Danp3WdIF7Tfrvuq1BzmT8Dfle4N8xTpHNwROCIPIvCMiPiLpLuBI0hDdF/R0L1I21ifhyv+IzBH0s9I5/Bw0nC7k0jDb/d2vQ9Iuo+U2P2DNIrjPRFxQzcvvRZ4J7As0lDVRMQTkmaSPuc3RcTjPQjhJuAtkr5Aaib5YEQ0dZjgiFgr6bOkoY3/Lulsnh1eeggpoass+5ik84BP5gv1O0hDVr+B1J+w6ApJy0mfuSWkgVg+RUoSG3WO/BcpGfutpOLw0pWErae1ZvcA35S0G+l74WDSIAQ/jYg5dV/ZhXzOfhT4A3Br/s5ZQBrK+nDSObuwF98PrwD+LOniHONK0jDuE+lZawazluVExzpaRNwuaXfSPSYm5ccQUoJxJWlI5SV52YeV7sfxLVJ/gheS/vl8ntQvphHxXKJ0L5KTSP/YX0C6P8YtpAvqej5AGjL1KNKF6q3AYaT75uxftexJpBHePp+3cS9d1JY0Y78L21ot6QDg66R/6MeQLuK+AXyjOMJZM7cVEXOVbrx6CikRGkZ6X64hNzPJF+3vI71P3yENDfsbUof3en2ruovz8/ki+FOkQSk2IV2UXcdzR8VqlGNIw2+fAjyP1OTqMxFRfUFaHet1uRbx68Bk0i/1s0lDdV/excvOITWf+WVseO+chomIv+X370ukC7+XkC6a55CSs411JKk/zXdJydK5pD5a9VxNSnSqm6ddQ0p0umu2VvEJ4Mek+yU9P7+u6fdDiYiLJK0knfdfI9VI/Q14T0T8vWrxz5KuO44kncdXke4hVL3PZ5A+Z8eRmpE+QOr/9vVcw9qIuB+QtB/wfVJC9iSpmdg5Of6e9ttcmmP9Hqlm51+k78qv9DG+a/L37pdJ93HanGd/DHiosFxPvh8WAb8k1aZVfny5h1Tb/MO+xGk20KmB1w1mZtbCJJ1CukAbE8+9EWl/bfMEUk3RLjWGlTZrqtxE9xJgQkTUHeJb0kJgcUS8vhmxmVnvuY+OmZmVQtImwEeBvzrJsWar7reS++gcTxqx8uZSgjKzhnLTNTMza6o8sMZbgINIgxz0ZcAGs431F0mzSM18tyD1y9wT+EIvbjlgZgOYEx0zM2u2scB5pA78X4uIS0uOxzrTdODdpFsMDCbdGPkjEfGzUqMys4ZxHx0zMzMzM2s77qNjZmZmZmZtx4mOmZmZmZm1HSc6ZmZmZmbWdpzomJmZmZlZ23GiY2ZmZmZmbceJjpmZmZmZtR0nOmZmZmZm1nac6JiZmZmZWdtxomNmZmZmZm3HiY6ZmZmZmbUdJzpmZmZmZtZ2nOiYmZmZmVnbcaJjZmZmZmZtx4mOmZmZmZm1HSc6ZmZmZmbWdpzomJlZW5M0WlJIen3ZsZiZWfM40TEzs6bIyUa9x8JuXn+WpKubE62ZmbW6wWUHYGZmHWNE4flewKX576Jctq7pEZmZWdtyjY6ZmTVFRNxfeQCP5OLlhbI9Jc2StFrSg5J+LGlzAEmnAB8G9ivUAB2d5x0v6RZJKyXdL+kCSSNqhGBmZh3EiY6ZmZVO0quA6cC1wGuAo4BDgTPzIt8Bfg38jVQzNAK4sLCKLwC7A+8ARgEXNCNuMzMbuNx0zczMBoITgJsj4j/y9FxJnwZ+K+lLEXGvpFXAmlz784yI+EFh8h5JnwJuljQyIpY0J3wzMxtoXKNjZmYDwW6k2pyiawABY+u9UNL+kmZIWiTpceC6PGvHxodpZmatwomOmZkNFNHLciSNAi4HFgJHAOOBt+fZmzYyODMzay1uumZmZgPBncB+VWX7kZKcOXl6DTCoapk9gaHAZyNiFYCk1/ZjnGZm1iJco2NmZgPBacAekr4naRdJBwH/A5wXEfflZe4BdpG0m6ThkjYD7iYlQ5+XtJOkw4Avl7EDZmY2sDjRMTOz0kXEbaQmZ/sBtwK/BC4DPl5Y7GfATcBfgeXA+/LrPg18jFTz8wXgs00L3MzMBixFdNn02czMzMzMrCW5RsfMzMzMzNqOEx0zMzMzM2s7TnTMzMzMzKztONExMzMzM7O2M6DvozN8+PAYPXp02WGYmZmZmdkANWvWrIciYqvq8gGd6IwePZqZM2eWHYaZmZmZmQ1Qku6tVd6QpmuSDpI0T9J8SZNrzN9f0qOSbskP38zNzMzMzMz6TZ9rdCQNAk4HDgQWAzdJmh4Rc6oW/UtEHNrX7ZmZmZmZmXWnETU6ewHzI2JBRKwBLgAmNWC9ZmZmZmZmG6URic5IYFFhenEuq7avpFsl/V7Sbl2tTNKxkmZKmrl8+fIGhGdmZmZmZp2mEYMRqEZZVE3fDOwYESslHQJMA8bUWllETAWmAowfP756PU03bfYSTpsxj6UrVrHdsKGcMHFnDhtXK48zMzMzM7OBohE1OouBHQrT2wNLiwtExGMRsTI/vxwYIml4A7bdr6bNXsJJl9zOkhWrCGDJilWcdMntTJu9pOzQzMzMzMysjkYkOjcBYyTtJGlT4AhgenEBSdtKUn6+V97uww3Ydr86bcY8Vj29boOyVU+v47QZ80qKyMzMzMzMeqLPTdciYq2k44AZwCDg7Ii4U9LH8/wzgXcBn5C0FlgFHBERpTdL687SFat6VW5mZmZmZgNDQ24YmpujXV5Vdmbh+Y+AHzViW8203bChLKmR1Gw3bGgJ0ZTH/ZTMzMzMrNU05Iah7eqEiTszdMigDcqGDhnECRN3Limi5nM/JTMzMzNrRU506jhs3EhOPXx3Nh2UDtPIYUM59fDdO6o2w/2UzMzMzKwVNaTpWjs7bNxIzr/xPgAu/Ni+JUfTfO6nZGZmZmatyDU6VldX/ZE6rZ+SmZmZmbUWJzpWl/spmZmZmVkrctM1q6vSH+nEi29jzbr1jPSoa2ZmZmbWApzoWLc6vZ+SmZmZmbUeN10zMzMzM7O24xods274hqlmZmZmrceJjlkdlRumVu4lVLlhKuBkx8zMzGwAc9M1szp8w1QzMzOz1uREx6wO3zDVzMzMrDU50TGrwzdMNTMzM2tNTnTM6vANU83MzMxakwcjMKvDN0w1M7MKj8Jp1lqc6Jh1wzdMNTMzj8Jp1nrcdM3MzMysGx6F06z1NCTRkXSQpHmS5kuaXGO+JP0wz79N0h6N2K6ZmZlZM3gUzlSrNWHKVew0+TImTLmKabOXlB2SWV19bromaRBwOnAgsBi4SdL0iJhTWOxgYEx+7A2ckf+amZlZC+j0/inbDRvKkhpJTaeMwumme9aKGtFHZy9gfkQsAJB0ATAJKCY6k4BfREQAf5c0TNKIiFjWgO33u4Ou/jXbLl/Evde9sOxQSnP0sscAOvYYdPr+m1lne2jlap5+6Ak+sz6eKXv6MjFr+OYM32KzEiNrnh+uXM2Ch55gfeEYbLKJeOnwzbn3yHPKC6xZ7lvBKWvXPbf82kHcO2pY08Oxcm226y5se/LJZYfRrUYkOiOBRYXpxTy3tqbWMiOB5yQ6ko4FjgUYNWpUA8Lru4N3H8HquY+WHUapxo7o7Av8Tt//h1au5p/LnyAi2GzwIHbYcmjHXNxUzMnJbiefC51+DDp5/xc9smqDC3yA9euDRY+s6pjvgsp+LnpkFavXruu478LVtZKcOuXtrJO/CyDt//1rl3FM2YH0QCMSHdUoi41YJhVGTAWmAowfP77mMs3WChmrWX+pbq4A6V5Cpx6+e0c1VzjxJ38DOnvkvU4/Bp28//tPvqzmP20B90x5W7PDKc2OwGvLDqIk759yVc2meyOHDeW9k99cQkTlmDZ7ScffcqLyXdgKiU4jBiNYDOxQmN4eWLoRy5jZAOSRhsysq34ondI/xXwDbXj2h78169YDz/ZT8qAMA1cjEp2bgDGSdpK0KXAEML1qmenAh/Loa/sAj7ZK/xyzTueRhszMF7l22LiRnHr47owcNhSRanI6rWbfP/ylZG/2fSu44Z5HWmLkvT43XYuItZKOA2YAg4CzI+JOSR/P888ELgcOAeYDT9IatV1mhkcaMrNnR9Xq5FHXLJ0Hnfyed/oPf13VaMHAHXmvEX10iIjLSclMsezMwvMAPtWIbZlZc50wceeafXT8S65ZZ+n0i1yzTv/hr16N1kD9bmjIDUPNrH25uYJZ6zXXMLPG6/QmnK1Yo9WQGh0za2/+Jdc6WSs21zCzxuv0JpytWKPlRMfMzKyOVmyuYWb9o5N/+GvFpuxOdMzMzOpoxeYaZmaN1oo1Wk50zMy6UemfsWbdeiZMuWrAf7FbY7Vicw0zs/7QajVaHozAzKwO3yDOOr0DsplZq3KiY2ZWh28Ql3TyqGMeedDMrDW56ZqZWR3un+FRx6D1mmuYmZlrdMzM6uqqH0Yn9c9wrZaZmbUiJzpmZnW4f4ZrtczMrDU50TEzq8P9M1yrZWZmrcl9dMzMutHp/TNa8SZxZmZmTnTMzKyuVrxJnJmZmRMdMzPrVqfXapmZWetxHx0zMzMzM2s7TnTMzMzMzKztONExMzMzM7O206c+OpK2BC4ERgMLgfdExL9qLLcQeBxYB6yNiPF92a6ZmZmZmVk9fa3RmQxcGRFjgCvzdFfeFBGvcZJjZmZmZmb9ra+JziTg3Pz8XOCwPq7PzMzMzMysz/qa6GwTEcsA8t+tu1gugD9KmiXp2HorlHSspJmSZi5fvryP4ZmZmZmZWSfqto+OpD8B29aY9cVebGdCRCyVtDVwhaS7IuLaWgtGxFRgKsD48eOjF9swMzMzMzMDepDoRMRbupon6QFJIyJimaQRwINdrGNp/vugpN8CewE1Ex0zMzMzM7O+6mvTtenAUfn5UcCl1QtI2lzSCyrPgbcCd/Rxu2ZmZmZmZl3qa6IzBThQ0t3AgXkaSdtJujwvsw1wnaRbgRuByyLiD33crpmZmZmZWZf6dB+diHgYOKBG+VLgkPx8AfDqvmzHzMzMzMysNxQxcPv7S1oO3Ft2HNlw4KGyg7BS+RwwnwPmc8DA54H5HBhodoyIraoLB3SiM5BImumbnXY2nwPmc8B8Dhj4PDCfA62ir310zMzMzMzMBhwnOmZmZmZm1nac6PTc1LIDsNL5HDCfA+ZzwMDngfkcaAnuo2NmZmZmZm3HNTpmZmZmZtZ2nOj0gKSDJM2TNF/S5LLjseaStIOkP0uaK+lOSceXHZOVQ9IgSbMl/a7sWKz5JA2TdLGku/L3wb5lx2TNJek/8v+BOySdL+l5Zcdk/UvS2ZIelHRHoWxLSVdIujv/fXGZMVrXnOh0Q9Ig4HTgYGAs8D5JY8uNyppsLfD5iNgV2Af4lM+BjnU8MLfsIKw0PwD+EBG7kG6E7XOhg0gaCXwGGB8RrwQGAUeUG5U1wTnAQVVlk4ErI2IMcGWetgHIiU739gLmR8SCiFgDXABMKjkma6KIWBYRN+fnj5MubkaWG5U1m6TtgbcBZ5UdizWfpBcCbwR+BhARayJiRalBWRkGA0MlDQaeDywtOR7rZxFxLfBIVfEk4Nz8/FzgsGbGZD3nRKd7I4FFhenF+CK3Y0kaDYwDbig5FGu+7wMnAutLjsPK8VJgOfDz3HzxLEmblx2UNU9ELAG+A9wHLAMejYg/lhuVlWSbiFgG6cdQYOuS47EuONHpnmqUeai6DiRpC+B/gc9GxGNlx2PNI+lQ4MGImFV2LFaawcAewBkRMQ54AjdX6Si5H8YkYCdgO2BzSR8sNyozq8eJTvcWAzsUprfHVdUdR9IQUpJzXkRcUnY81nQTgLdLWkhqvvpmSb8qNyRrssXA4oio1OZeTEp8rHO8BbgnIpZHxNPAJcDrSo7JyvGApBEA+e+DJcdjXXCi072bgDGSdpK0Kanj4fSSY7ImkiRSu/y5EfG9suOx5ouIkyJi+4gYTfoOuCoi/EtuB4mI+4FFknbORQcAc0oMyZrvPmAfSc/P/xcOwANSdKrpwFH5+VHApSXGYnUMLjuAgS4i1ko6DphBGmHl7Ii4s+SwrLkmAEcCt0u6JZedHBGXlxeSmZXg08B5+UevBcAxJcdjTRQRN0i6GLiZNBrnbGBquVFZf5N0PrA/MFzSYuArwBTgIkkfJiXA7y4vQqtHEe5uYmZmZmZm7cVN18zMzMzMrO040TEzMzMzs7bjRMfMzMzMzNqOEx0zMzMzM2s7TnTMzMzMzKztONExMzMzM7O240THzMzMzMzajhMdMzMzMzNrO050zMzMzMys7TjRMTMzMzOztuNEx8zMzMzM2o4THTMzswJJIemDZcdhZmZ940THzKzDSXqJpG9LmifpKUkPSrpW0ockDW5SDK/PCcboZmwvb/NPks6pMWsEcHGz4jAzs/7RlH9gZmY2MEnaHrgeWAt8GZgNPA28DvgCcBtwS1nxVZM0BFgbEdFf24iI+/tr3WZm1jyu0TEz62xnAJsBe0TEeRExJyLujohzgdcCd0saImmKpCWS1kiaI+n9xZXk2phPSvqlpMclLZJ0YtUykyTNlvSkpBWSbpQ0Ltfi/CUvdk9e19X5NefkmpdPS1oIrAY2l3S1pLOq1v+lvEyx7L2SZuWaqocl/V7Si3NNzgHAUXl7IWn/wr58sLCOEZIuyDGvytseX5i/f37Ngbkm7Ml8jCZu1DtiZmYN4UTHzKxDSdoSOAT4UUQ8Wj0/Ip6OiCeAbwEfBT4LvBL4FfArSQdUveQrwLXAa4DTgP+S9Ka8rW2B3wDnA7sB+wLfJ9UkLQIm5XXsRWo6dnhhvXsBbwYOA14NPNXD/TsmxzoN2AN4E/AHYBBwPCm5uihvbwTw1xrrUH79LsChOZYHgCskDa9a/DukY/VqYCZwoaRhPYnVzMwaz03XzMw618tJP3jN6WoBSc8HPgP8R0T8Jhd/S9KewBeBKwuLXxgRP83Pfyjpk8BbgT+TEokhwEURsTAvM7ewnUfy0+U1mo6tB46MiJWF5Xuyf18FfhIRXy+U3VZYxxpgVTdN1d5MSm52i4g5+XUfAhYCnwS+VtxeRPwhL3MicCSwNzCjJ8GamVljuUbHzKxzVbKFev1dXg5sSqqpKbqGVDNTdEvV9BJgm/z8NtIF/x2SfivpeEk79DDOucUkpyckbQ3sAPyxN6+rYTfg4UqSAxARq4EbqLP/OXlax7P7b2ZmTeZEx8ysc91Nqi2pvmCvpToZUo2yNTVeswlARKwDDibVkNwEvBP4h6RDe7DtJ2qUrefZRK1iSI3lGjFoQa119GT/wf9nzcxK4y9gM7MOFRGPAL8HjpP0our5eYSzf5IGANivavYbgTt7ub2IiBsj4lsR8UZSrdAxeXYlSRjUw9U9CGxXVbZHYVsPAouBegMCrOnB9u4EhksaWymQtBmpOVuv9t/MzJrLiY6ZWWf7JGk46VmS3i9prKSX51HHZgJjgB8CX5f0bkljJJ1MGjzgWz3diKTXSfpPSXtLGpUHMngVz/YPupdUS3OIpK1rJV5V/gS8RdJ7cryTgTdULfNV4GN5u7tK2k3ScYVBBO4BXivpZZKG58Su2lXAjcCvJU2Q9ErgF8DzSCPWmZnZAOVEx8ysg0XEfaSakEuBU4CbSaOPfZQ0ctodpEEHfkoaJe1O4IPAByPiyueusUuPkkZau5TUZO5s4Dzg6zmOB4CTgMnAsrxcPecCpwM/IiVkO5ASsuK+nQUcDbyL1H/mWlLzubV5ke8CDwG3AsuBCdUbyffrOQy4C7iM1OxuW+DAiHioJztuZmblUD/ec83MzMzMzKwUrtExMzMzM7O240THzMzMzMzajhMdMzMzMzNrO050zMzMzMys7QwuO4B6hg8fHqNHjy47DDMzMzMzG6BmzZr1UERsVV0+oBOd0aNHM3PmzLLDMDMzMzOzAUrSvbXKG9J0TdJBkuZJmp9v2lY9f39Jj0q6JT++3IjtNsO02UuYMOUqdpp8GROmXMW02UvKDsnMzMzMzLrR5xodSYNIN207EFgM3CRpekTMqVr0LxFxaF+310zTZi/hpEtuZ9XT6wBYsmIVJ11yOwCHjRtZZmhmZmZmZlZHI2p09gLmR8SCiFgDXABMasB6S3fajHnPJDkVq55ex2kz5pUUkZmZmZmZ9UQjEp2RwKLC9OJcVm1fSbdK+r2k3Rqw3X63dMWqXpWbmZmZmdnA0IjBCFSjLKqmbwZ2jIiVkg4BpgFjaq5MOhY4FmDUqFENCG/jbTdsKEtqJDXbDRtaQjTlmTZ7CafNmMfSFavYbthQTpi4s5vumZmZmdmA1ogancXADoXp7YGlxQUi4rGIWJmfXw4MkTS81soiYmpEjI+I8Vtt9ZxR4prqhIk7M3TIoA3Khg4ZxAkTdy4pouar9FNasmIVwbP9lDwog5mZmZkNZI1IdG4CxkjaSdKmwBHA9OICkraVpPx8r7zdhxuw7X512LiRnHr47mw6KB2mkcOGcurhu3dUbYb7KZmZmZlZK+pz07WIWCvpOGAGMAg4OyLulPTxPP9M4F3AJyStBVYBR0REdfO2AemwcSM5/8b7ALjwY/uWHE3zuZ+SmZmZmbWihtwwNDdHu7yq7MzC8x8BP2rEtqy53E/JzMzMzFpRQ24Yau3L/ZTMzMzMrBU1pEbH2lelP9KJF9/GmnXrGelR18zMzMysBTjRsW51ej8lMzMzM2s9brpmZmZmZmZtx4mOmZmZmZm1HSc6ZmZmZmbWdpzomJmZmZlZ2/FgBGbdmDZ7CafNmMfSFavYzqPOmZmZmbUEJzpmdUybvYSTLrmdVU+vA2DJilWcdMntAE52zMzMzAYwN10zq+O0GfOeSXIqVj29jtNmzCspIjMzMzPrCSc6ZnUsXbGqV+VmZmZmNjA40TGrY7thQ3tVbmZmZmYDgxMdszpOmLgzQ4cM2qBs6JBBnDBx55IiMjMzM7Oe8GAEZnVUBhw48eLbWLNuPSM96pqZmZlZS3CiY9aNw8aN5Pwb7wPgwo/tW3I0ZmZmZtYTbrpmZmZmZmZtx4mOmZmZmZm1nYY0XZN0EPADYBBwVkRMqZqvPP8Q4Eng6Ii4uRHbNjMzM2uGabOXcNqMeSxdsYrtOrDPZqfvv7WePic6kgYBpwMHAouBmyRNj4g5hcUOBsbkx97AGfmvmbUA/3Mzs043bfYSTrrk9mduIr1kxSpOuuR2gI74Puz0/bfW1Iima3sB8yNiQUSsAS4AJlUtMwn4RSR/B4ZJGtGAbZtZP6v8c1uyYhXBs//cps1eUnZoZmZNc9qMec9c5Fesenodp82YV1JEzdXp+2+tqRFN10YCiwrTi3lubU2tZUYCyxqw/X530NW/Ztvli7j3uheWHUppjl72GEDHHoOO3v/7VnDK2nXPLb92EPeOGtb0cMzMyvCZBQ93Oe/eO89pXiAl6fT9r3ho5WoWPbKK1WvXsdngQeyw5VCGb7FZ2WE13Wa77sK2J59cdhjdakSioxplsRHLpAWlY4FjAUaNGtW3yBrk4N1HsHruo2WHUaqxIzrwAr+gk/d/da0kp055u5qTk91OPhc6/Rh0+v5DZx+DzQYPqvm9t9ngQTWWbj+dvv+QkpwFDz3B+vXpEnb12nUseOgJgI5KduYse4z71y7jmLID6YFGJDqLgR0K09sDSzdiGQAiYiowFWD8+PE1k6Fma4WM1ay/vH/KVSxZseo55SOHDeW9k99cQkTlOPEnfwM6+15KnX4MOn3/obOPwezZSzil0EcFYOiQQZx6+O7s2AF9VDp9/6H+/8PrO+T/4bTZS565ifpZU64a8H12G9FH5yZgjKSdJG0KHAFMr1pmOvAhJfsAj0ZESzRbM+t0J0zcmaFDNvzFbuiQQZwwceeSIjKzMkybvYTZ963ghnseYcKUqzqun95h40Zy6uG7M3LYUES6uD318N0H9EVeI3X6/gMsrZHk1CtvN5U+u2vWrQdao89un2t0ImKtpOOAGaThpc+OiDslfTzPPxO4nDS09HzS8NKtUNtlZjw7mo5HXTPrXF1d4EBnjbh12LiRHbW/1Tp9/7cbNrRmjc52w4aWEE3z1RuQYqCeFw25j05EXE5KZoplZxaeB/CpRmzLzJqv0/+5mXW6VrzAMWu0EybuvMEQ29BZLRxasUarEU3XzMzMrI214gWOWaN1evO9rmquBnKNVkNqdMzM2lmlb8KadeuZ0AKdL80ardOb7JhVdHILh1as0XKNjplZHa3Y+dKs0TwoiZm1Yo2Wa3TMzOpw3wQD1+p5UBIzg9ar0XKiY2ZWh/smmEccS1rtAsfMzE3XzMzqaMXOl9ZY9Wr1zMxs4HKiY2ZWh/smmGv1zMxakxMdM7M6WrHzZX+o9FG54Z5HmDDlqo4ajMG1emZmrcl9dMzMutHpfRM6vY9KKw6pamZmrtExM7NudHofFdfqmZm1JtfomJlZXe6j4lo9M7NW5BodMzOry31UzMysFTnRMTOzujzynJmZtSI3XTMzs7oqTbZOmzGPpStWsd2woZwwcWc35TIzswHNiY6ZmXXLfVTMzKzVuOmamZmZmZm1HSc6ZmZmZmbWdvrUdE3SlsCFwGhgIfCeiPhXjeUWAo8D64C1ETG+L9s1MzMzMzOrp681OpOBKyNiDHBlnu7KmyLiNU5yzMzMzMysv/U10ZkEnJufnwsc1sf1mZmZmZmZ9VlfE51tImIZQP67dRfLBfBHSbMkHdvHbZqZmZmZmdXVbR8dSX8Ctq0x64u92M6EiFgqaWvgCkl3RcS1XWzvWOBYgFGjRvViE2ZmZmZmZkm3iU5EvKWreZIekDQiIpZJGgE82MU6lua/D0r6LbAXUDPRiYipwFSA8ePHR/e7YGZmZmZmtqG+Nl2bDhyVnx8FXFq9gKTNJb2g8hx4K3BHH7drZmZmZmbWpb4mOlOAAyXdDRyYp5G0naTL8zLbANdJuhW4EbgsIv7Qx+2amZmZmZl1qU/30YmIh4EDapQvBQ7JzxcAr+7LdszMzMzMzHpDEQO3G4yk5cC9ZceRDQceKjsIK5XPAfM5YD4HDHwemM+BgWbHiNiqunBAJzoDiaSZvtlpZ/M5YD4HzOeAgc8D8znQKvraR8fMzMzMzGzAcaJjZmZmZmZtx4lOz00tOwArnc8B8zlgPgcMfB6Yz4GW4D46ZmZmZmbWdlyjY2ZmZmZmbceJTg9IOkjSPEnzJU0uOx5rLkk7SPqzpLmS7pR0fNkxWTkkDZI0W9Lvyo7Fmk/SMEkXS7orfx/sW3ZM1lyS/iP/H7hD0vmSnld2TNa/JJ0t6UFJdxTKtpR0haS7898Xlxmjdc2JTjckDQJOBw4GxgLvkzS23KisydYCn4+IXYF9gE/5HOhYxwNzyw7CSvMD4A8RsQvpRtg+FzqIpJHAZ4DxEfFKYBBwRLlRWROcAxxUVTYZuDIixgBX5mkbgJzodG8vYH5ELIiINcAFwKSSY7ImiohlEXFzfv446eJmZLlRWbNJ2h54G3BW2bFY80l6IfBG4GcAEbEmIlaUGpSVYTAwVNJg4PnA0pLjsX4WEdcCj1QVTwLOzc/PBQ5rZkzWc050ujcSWFSYXowvcjuWpNHAOOCGkkOx5vs+cCKwvuQ4rBwvBZYDP8/NF8+StHnZQVnzRMQS4DvAfcAy4NGI+GO5UVlJtomIZZB+DAW2Ljke64ITne6pRpmHqutAkrYA/hf4bEQ8VnY81jySDgUejIhZZcdipRkM7AGcERHjgCdwc5WOkvthTAJ2ArYDNpf0wXKjMrN6nOh0bzGwQ2F6e1xV3XEkDSElOedFxCVlx2NNNwF4u6SFpOarb5b0q3JDsiZbDCyOiEpt7sWkxMc6x1uAeyJieUQ8DVwCvK7kmKwcD0gaAZD/PlhyPNYFJzrduwkYI2knSZuSOh5OLzkmayJJIrXLnxsR3ys7Hmu+iDgpIraPiNGk74CrIsK/5HaQiLgfWCRp51x0ADCnxJCs+e4D9pH0/Px/4QA8IEWnmg4clZ8fBVxaYixWx+CyAxjoImKtpOOAGaQRVs6OiDtLDsuaawJwJHC7pFty2ckRcXl5IZlZCT4NnJd/9FoAHFNyPNZEEXGDpIuBm0mjcc4GppYblfU3SecD+wPDJS0GvgJMAS6S9GFSAvzu8iK0ehTh7iZmZmZmZtZe3HTNzMzMzMzajhMdMzMzMzNrO050zMzMzMys7TjRMTMzMzOztuNEx8zMzMzM2o4THTMzMzMzaztOdMzMzMzMrO040TEzMzMzs7bjRMfMzMzMzNqOEx0zMzMzM2s7TnTMzMzMzKztONExM+sQkhZK+lID1nOKpPmNiKlMjToe/UFSSPpg2XHU0woxmllnc6JjZi1P0raSnpJ0v6QhZcfTLiS9Pl/Mjq6a9R1gnxJCaipJ++f9r378oezYACQdI2mWpMckPS5prqSfNjGEEcDFTdyemVmvDC47ADOzBvh34DJgF2ASvvjqVxGxElhZdhxNtAewrDC9uqxAKiQdDZwJfB74Yy7elXT+92W9mwCKiHXdLRsR9/dlW2Zm/c01OmbW0vKF2UeBc/Pj2BrLbC3p55IeyDU/8yT9e2H+yyT9RtIjkp6UdJukQwvzD8m/nK+W9KCkH0vavDD/HEl/kvRpSYslrZR0lqQhkj4u6V5J/5I0VdKmhdddLelnkr6R17tC0jclbSLpyzne5ZK+WbU/L5D0kzzvKUkzJb21aplXS/prnv8PSe+pcVyOl3RLjvd+SRdIGpHnjQb+khe9J9dkXJ3nPafpmqSjJM3Jx2hx3qfBhflX52Pyn3lbj+Tjtjl11Isxz6/Uuhwo6dr8/s2RNLG3x6OO5RFxf+Hxr7zOIZKmSFoiaU3e7vurtjsix7xC0qp8HMZXLfOmfM49lf++qQcxHQb8X0T8KCL+kR+XRsS/FxeS9FpJf8zHb7mkSyTtWJh/iqT5kt4r6S5gDfBJSesk7VC1rvfmGIfl6Q2arknaQtL3JS3K58FCSScX5m+T3/PlSjVQ10t6Y2H+EEnfy+fPaknLJF3Qg2NhZlaTEx0za3VvBTYHLgd+Cewv6aWVmZKGAtcArwY+AIwFPg08medvC/wVeDHwdmB34D+B9Xn+q4DpwLXAa4CjgENJv6YX7QmMBw4E3g98ELgUeB1wMHBkfny46nXvAoYArwc+B5wM/A7YAngD8AXgZEkHF15zNjAxb2MccD3wO0m7FPb5cmAFsHeO+QRg6xrH7wt5n98BjAIqF5aLeLZ2YC9SM6XDa7weSW/LMf0yr+vzwKeAr9TY1y2B/fMxOgw4sdY6exhj0XeAb5He55nAhYUL8t4cj974FinJ/izwSuBXwK8kHZC3K2AaqabxUNJxfAC4QtLwvMx2pPd7Fqnm6PPAD3qw7WXAeEmv6GoBSWNJ5/7fSOfmm4F1efvPKyy6HfBJ4GjS5+MXef3V/W+OBC6NiBU1tqW8H28nfb52BT4ELM/zhwJ/Bl5A+jyMI70nV0jaNa/m08B78nbH5HX9vbsDYWbWpYjwww8//GjZB/Bb4L8L05cD3ypMfxh4Cti+i9d/Hbgf2LyL+b8Ebqwqm0RKhHbM0+cADwKbFpa5DHgI2KxQdilwcWH6auCWqnXfCdxeVXYr8J38/OVAAIdULXMzcHZ+/hFS07IXF+a/Mr/uS3WO5bi8zMg8/fo8PbpquVOA+YXpvwAXVS1zPLCqckzyvt5WtcyZwN96+X5Xx7h/nj68sMy2uWxiH49HZd1P5NdXHm8Bnk9qwvbJGufjVfn5Afn1YwvzNyMlEV/O098A7gUGF5Y5NL/ug3Vi2xa4Li+3ELiQVJv5/MIy5wAXVL1uM1KSf1jhvVwPjKpabgowpzC9NfA08LZC2TMxFvZ1fBfxHg0sLu5nLr8K+H5+/oM8rd6cE3744YcfXT1co2NmLSs3YTqU1GSt4hzgmEKzqdeSLtgWd7Ga1wJ/jYgnupi/G6k2p+gaQKRfvyvmRsSawvT9wLyIWF1VVl2LcGvV9P3AbTXKKq+rbLM6pmtzrJVl5kZuYgUQEXcAjxZfkJt9zchNjR4nXTgD7EjvdHWMnge8rFB2S9UyS4Bt6q24FzE+s+5IfUfWFdbdo+NRx0RSbV7l8VdSwrkptfe78j7sBjwcEXMK210N3MCG79WNEbG2sI7r6EakJnSvz68/lZSMfRu4U1LlXNkTeEdutrZS0krgYdL7Mqawugci4r6qTZwL7Cppzzz9vvzaGV2E9FrgXxExs4v5e5KSsxVV8byhEMvPSTV38yWdKemdKjT1NDPrLQ9GYGat7MOk77GZqeXMMwaRmr1ckqejm/Vs7Pxi+dM15tUqq/6BaWNfV02FeIrPay8sjeLZ5n5fI9U+bQ/8iXQB31vV21ON8jVVy9Tdr17GWL1uCuvu9nh0Y2F1olw432rtd7Gs1na7e696HGtEzAXmAj+R9HXgH8AngK+S9v+XpNqZag8Xnj8nyY+IuZJmkpqf3ZT//roqIXvOy+rM2yTH+Y4a857M27xF0k6k5p9vItXwfF3SPhHxWJ11m5nV5BodM2tJSoMQfITUT+I1VY9f8eygBLOA3SRt38WqZgET6nSKvxPYr6psP9JF3ZznLt7v7sx/31hV/obCvDuBsZU+KgCSdgNeVFh+T2Ao8NmIuD4i5vHc2pVK8jCoBzFVH6M3kpquLejmtfX0JMae6Mnx6K35pKZrtfa7+D4Mz31lKtvdjNRXp7jM3pKKx/j1GxnTQlLSUKnRmQm8CvhnRMyvevyrq5UU/AI4QtKrSf2Hzq2z7Cxgy+qBFgpmAi8FHqsRy9LKQhGxMiJ+GxGfIfUr2pXnHmMzsx5xomNmreogUsf0n0TEHcUHqQnMgUojh51P6gMxXdJbJO0k6QBJ783r+THpu/BSSRPy/EMLnf9PA/bIo0HtIukg4H+A82o09+l3EfFP4DfAjyVNzDH9gNTn5LS82K+Bx0kd418taR/SYAGrCqu6m5SsfT7v82HAl6s2dy+p/8YhSiPXdZUYnAq8U9JkSa9QGtHsFOC7Vc35eqsnMfZET45Hr0TEk8APSTUO75Y0Jo8wNomUfEPqb3Ij8Ot8br2SlDw8DzgjL3MGsBUwVdKueSCDDUbZq0XSGZK+IukNknaU9FpSIvJC0gAI5Dh2zfu9Vz6Gb5L0g+KAHXWcT0oGzyH1r6puZll0Famv1oWSJuVtTZD0kTz/POAe4DJJb5U0WtLekk7K7yuSTpD0AUm75Zqdfyc1QfxHD2I1M3sOJzpm1qo+BtzQRbJxDWm0p4/kC9L9gDtIo3XNBU4n1RQQEctIv6A/TmomdSfpQlN5/m2kZnD7kfrT/JI00MDH+2vHeuAjpL4Sv8oxTQAOjYi74JmL8EOAl5AutM8D/ps0YAJ5mdtIo1x9jFQz9QXS6GEUlnkAOAmYTOpAf2mtYCLictJF6VGk4/zfpATyq33ZyZ7E2MP1dHs8NtIXgZ8C3yedNx8kdc6/Mm83SCPL3UU6Z24i9VM5MCIeysssAf6NVMtzC6m51ud6sO0rSP1iziclApeTRsY7JCKuyOueSxr1bwvS+TInxzuUNAJdXTnGy0i1pL/oZtkA3pbjOBOYRzo/h+f5T5E+QzNJP0T8g9S0dC9SQg3wGGnf/wbcTmrm9s5ck2dm1mtK301mZmZmZmbtwzU6ZmZmZmbWdpzomJmZmZlZ23GiY2ZmZmZmbceJjpmZmZmZtZ0BfcPQ4cOHx+jRo8sOw8zMzMzMBqhZs2Y9FBFbVZcP6ERn9OjRzJw5s+wwzMzMzMxsgJJ0b63yhjRdk3SQpHmS5kuaXGP+/pIelXRLfmzMDd/MzMzMzMx6pM81OpIGkW6+dyCwGLhJ0vSImFO16F8i4tC+bs/MzMzMzKw7jajR2QuYHxELImIN6c7jkxqwXjMzMzMzs43SiERnJLCoML04l1XbV9Ktkn4vabcGbNfMzMzMzKymRgxGoBplUTV9M7BjRKyUdAgwDRhTc2XSscCxAKNGjWpAeGZmZmZm1mkaUaOzGNihML09sLS4QEQ8FhEr8/PLgSGShtdaWURMjYjxETF+q62eM0qcmZmZmZlZtxqR6NwEjJG0k6RNgSOA6cUFJG0rSfn5Xnm7Dzdg22ZmZmZmZs/R56ZrEbFW0nHADGAQcHZE3Cnp43n+mcC7gE9IWgusAo6IiOrmbWZmZmZmZg2hgZxvjB8/PnzDUDMzMzMz64qkWRExvrq8ITcMNTMzMzMzG0ic6JiZmZmZWdtxomNmZmZmZm2nEffRsTY3bfYSTpsxj6UrVrHdsKGcMHFnDhtX656wZmZmZmYDgxMdq2va7CWcdMntrHp6HQBLVqzipEtuB3CyY2ZmZmYDlpuuWV2nzZj3TJJTserpdZw2Y15JEZmZmZmZdc+JjtW1dMWqXpWbmZmZmQ0EbrrWjU7vn7LdsKEsqZHUbDdsaAnRlKPTzwEzMzOzVuQanToq/VOWrFhF8Gz/lGmzl5QdWtOcMHFnhg4ZtEHZ0CGDOGHiziVF1Fw+B8zMzMxakxOdOtw/JQ04cOrhu7PpoHSqjBw2lFMP371jajR8DpiZmZm1Jjddq8P9U5LDxo3k/BvvA+DCj+1bcjTN5XPAzMzMrDW5RqeOrvqhdFL/lE7nc8DMzMysNTnRqaPT+6eYzwEzM7OKabOXMGHKVew0+TImTLnK/VVtwHPTtToq/VBOvPg21qxbz0iPuNVxfA6YmZn5BuLWmpzodKOT+6dY4nPAzMw6Xb3BeZzo2EDlRMfMzMysBzr5vmoenCfp5HOgFTnRMTMzM+tGpzfd8g3EfQ60ooYMRiDpIEnzJM2XNLnGfEn6YZ5/m6Q9GrFdMzMzs2bo9PuqeXAenwOtqM+JjqRBwOnAwcBY4H2SxlYtdjAwJj+OBc7o63bNzMzMmqXTm251+g3EwedAK2pE07W9gPkRsQBA0gXAJGBOYZlJwC8iIoC/SxomaURELGvA9s3MzMz6lZtueXAenwOtpxGJzkhgUWF6MbB3D5YZCbREonPQ1b9m2+WLuPe6F5YdSmmOXvYYQMceg07ffzOzTvfDlatZ8NATrF8fz5Rtsol46fDNuffIc8oLrMk6+f+hz4FnbbbrLmx78sllh9GtRiQ6qlEWG7FMWlA6ltS8jVGjRvUtsgY5ePcRrJ77aNlhlGrsiM77Qivq9P0HmJP/uXXqsej0/YfOPgYPrVzNP5c/QUSw2eBB7LDlUIZvsVnZYTVdJ58Dlfd70SOrWL12XceeB5343lf4HEjmLHuM+9cu45iyA+mBRiQ6i4EdCtPbA0s3YhkAImIqMBVg/PjxNZOhZmuFjNWsP02bvaTjb5p64k/+BnRmc42KTj0G1SMtQeqE3Wn9E/w9ADsCry07CCuVz4Fn/xe0QqLTiFHXbgLGSNpJ0qbAEcD0qmWmAx/Ko6/tAzzq/jlmraFykbdm3Xrg2eE0p81eUnJk1kzTZi9h9n0ruOGeR5gw5aqOev890pK/B8ysNfU50YmItcBxwAxgLnBRRNwp6eOSPp4XuxxYAMwHfgp8sq/bNbPm8EWedfpFrkda8veAmSWt9qNXQ24YGhGXk5KZYtmZhecBfKoR2zKz5vJFntW7yO2EpkseacnfA2bW9Y9eMHBvmNqQG4aaWfvq6mKuky7yOl2nX+T6Ron+HjCz1qzZdaJjZnX5Is86/SK3cqPEkcOGIjrzRon+HjCzVvzRqyFN18ysfVUu5k6bMY+lK1axXQeOtlRpk7xm3XomTLmq4/b/hIk71xx1rJMucg8bN7Kj3vNq/h4ws1ZsxqvUfWZgGj9+fMycObPsMMysg3lo4WTa7CW+yDUz62AD+f+hpFkRMb663DU6ZmZ1dHpH/IpOr9EwM+t0rViz60THzKyOVmyTbGZm1h9a7UcvD0ZgZlZHp3fENzMza1VOdMzM6vBoU2ZmZq3JTdfMzOpoxTbJZmZm5kTHzKxbrdYm2czMzNx0zczMzMzM2pATHTMzMzMzaztOdMzMzMzMrO040TEzMzMzs7bjRMfMzMzMzNqOEx0zMzMzM2s7TnTMzMzMzKzt9Ok+OpK2BC4ERgMLgfdExL9qLLcQeBxYB6yNiPF92a6ZmZmZmVk9fa3RmQxcGRFjgCvzdFfeFBGvcZJjZmZmZmb9ra+JziTg3Pz8XOCwPq7PzMzMzMysz/qa6GwTEcsA8t+tu1gugD9KmiXp2D5u08zMzMzMrK5u++hI+hOwbY1ZX+zFdiZExFJJWwNXSLorIq7tYnvHAscCjBo1qhebMDMzMzMzS7pNdCLiLV3Nk/SApBERsUzSCODBLtaxNP99UNJvgb2AmolOREwFpgKMHz8+ut8FMzMzMzOzDfW16dp04Kj8/Cjg0uoFJG0u6QWV58BbgTv6uF0zMzMzM7Mu9TXRmQIcKOlu4MA8jaTtJF2el9kGuE7SrcCNwGUR8Yc+btfMzMzMzKxLfbqPTkQ8DBxQo3wpcEh+vgB4dV+2Y2ZmZmZm1huKGLjdYCQtB+4tO45sOPBQ2UFYqXwOmM8B8zlg4PPAfA4MNDtGxFbVhQM60RlIJM30zU47m88B8zlgPgcMfB6Yz4FW0dc+OmZmZmZmZgOOEx0zMzMzM2s7TnR6bmrZAVjpfA6YzwHzOWDg88B8DrQE99ExMzMzM7O24xodMzMzMzNrO050ekDSQZLmSZovaXLZ8VhzSdpB0p8lzZV0p6Tjy47JyiFpkKTZkn5XdizWfJKGSbpY0l35+2DfsmOy5pL0H/n/wB2Szpf0vLJjsv4l6WxJD0q6o1C2paQrJN2d/764zBita050uiFpEHA6cDAwFnifpLHlRmVNthb4fETsCuwDfMrnQMc6HphbdhBWmh8Af4iIXUg3wva50EEkjQQ+A4yPiFcCg4Ajyo3KmuAc4KCqssnAlRExBrgyT9sA5ESne3sB8yNiQUSsAS4AJpUckzVRRCyLiJvz88dJFzcjy43Kmk3S9sDbgLPKjsWaT9ILgTcCPwOIiDURsaLUoKwMg4GhkgYDzweWlhyP9bOIuBZ4pKp4EnBufn4ucFgzY7Kec6LTvZHAosL0YnyR27EkjQbGATeUHIo13/eBE4H1Jcdh5XgpsBz4eW6+eJakzcsOyponIpYA3wHuA5YBj0bEH8uNykqyTUQsg/RjKLB1yfFYF5zodE81yjxUXQeStAXwv8BnI+KxsuOx5pF0KPBgRMwqOxYrzWBgD+CMiBgHPIGbq3SU3A9jErATsB2wuaQPlhuVmdXjRKd7i4EdCtPb46rqjiNpCCnJOS8iLik7Hmu6CcDbJS0kNV99s6RflRuSNdliYHFEVGpzLyYlPtY53gLcExHLI+Jp4BLgdSXHZOV4QNIIgPz3wZLjsS440eneTcAYSTtJ2pTU8XB6yTFZE0kSqV3+3Ij4XtnxWPNFxEkRsX1EjCZ9B1wVEf4lt4NExP3AIkk756IDgDklhmTNdx+wj6Tn5/8LB+ABKTrVdOCo/Pwo4NISY7E6BpcdwEAXEWslHQfMII2wcnZE3FlyWNZcE4Ajgdsl3ZLLTo6Iy8sLycxK8GngvPyj1wLgmJLjsSaKiBskXQzcTBqNczYwtdyorL9JOh/YHxguaTHwFWAKcJGkD5MS4HeXF6HVowh3NzEzMzMzs/bipmtmZmZmZtZ2nOiYmZmZmVnbcaJjZmZmZmZtx4mOmZmZmZm1HSc6ZmZmZmbWdpzomJmZmZlZ23GiY2ZmZmZmbceJjpmZmZmZtR0nOmZmZmZm1nac6JiZmZmZWdtxomNmZmZmZm3HiY6ZmXVL0kJJXyo7jr6QtL+kkLR92bH01ECKWdIpkuaXHYeZWU850TGztpUvEOs9FpYdY3+RtFbS0RvxurMkXV1j1p7Af/c1roFM0oclPS3pBVXlt9UpP7e5UZbqO8A+ZQdhZtZTTnTMrJ2NKDwm5bK9CmV7FheWtGlTo+sH/bUPEbE8Ip7oj3UPIH8CBgP7VQokDQd2A5bVKH8lcOXGbEjSEEnqU7RNFhErI+KhsuMwM+spJzpm1rYi4v7KA3gkFy8vlD0o6TOSfi3pUeA8AEnflDRX0pOSFkk6U9KLKuuVdHSuMZkg6ea83E2SXltYZoik70laLGm1pGWSLijMP0fSnyR9TtKSvI7/zRfQlWUk6QuSFkhaI+mfkj5b3MfcpOwbkn4s6WHg+lxTNQj4eaX2Ki/7Ykm/knSfpFWS5kn6fOWCW9IpwIeB/Qq1XkcXtvOlwnZfIOknkpZLekrSTElvLcwfnV//Hkn/l/dvgaQj671n3cVYdeyOlXSvpMckXSppq6p1fTof/yclzQBG1dt2RNwL/BM4oFD8ZuAO4NIa5QKu7GXMn87vz2pgc0nb5HnLJT0u6XpJb6wR3q6Srs37MkfSxKp93VnSZZJW5sf/SXp5Yf7RktZWvWb7/B7tn6e7O2c3aLpWmZY0SdJdkp6Q9GdJL6vazvvyufuUpL9KOjRv9/V13g4zsz5zomNmne4rwN+APYAv5rJVwLHAWOBoYH/gh1Wv2wQ4FTg+v/ZfwEWSBuf5nwbeA3wQGAO8Hfh71Tr2yus+CDgEeBVwdmH+J4GvA1NItQqnAVMkfbhqPZ8BHgT2BY4i1VStAz7Ls7VXAJsBtwOH5X37OvDVvI+Qmib9Oh+PyusupLazgYl5/8YB1wO/k7RL1XJTgF/mfbuIlHyN6WKdPYmxYk/gTcDbSMfvNTl+ACRNIjW1+16edxHp+HXnSp6b0FyVH9Xld0XEkl7EvFd+3WHAq0mJ0p+BFwAHk47j5cAVknateu13gG/l180ELpQ0LO/rUOCPwPNItU77AVsAf1Dvavh6cs5WGwF8AvgA8DpgGIVzWCn5Pw84P8f+beD7vYjJzGzjRYQffvjhR9s/gNcDAYwulAXwsx689h2kX+A3ydNH59fuUVhmn1y2c57+AeniWF2s8xxgJfCiQtlb8zrG5OlFwLerXvffwILC9ELgyhrrXwsc3YN9+wFwRWH6LODqGsstBL6Un788x3lI1TI3A2fn56PzMp8rzB+c9/ljvXzvqmM8B1gObFYomwwsK0xfB5xXtZ7v5Ji2r7Ot9wDrga3z9D9IF/xb5mNaLP+fXsa8AtiiUHY0sBgYXPXaq4Dv5+f755gPL8zfNpdNzNMfBp4EhheW2YaUsH+osK21VdvZPq9n/x6es6cA86um1wJbFcqOyMfveXn6POAvVev5eN7u63tzHvjhhx9+9PbhGh0z63Q3VhdIOjw3E1oqaSXpYm1T0gVmRQC3FqaX5L/b5L8/B3YH5is1fXtnjV/X50TEo4Xp6/PfXSW9kHQhem3Va64BRkt6fr19qEXSJpImS7pF0kN53z4O7NiT1xeMzX+rY7uWVPNUdEvlSUSsBR7g2WPUlxjnRsTqwvSSqvWOBf5a9ZrrutpuwVX57wFKI529FLgmIh4BbiuUjyH16eltzCsL03uSzqkVhSZnK4E35PUX3VJ5EqnZ5brC/u5GOpceKizzADCP574f9fTknK22NCKWF6aXkGqqts7TY3lurdDfehGTmdlGc6JjZp1ugw72kvYGfkO6aH8HqVnax/Ps4kXf+ohYV5iO/HcTgIi4BdgJ+AKwhvRr+S05gemNqJqu1YG9p4MEfB44Cfgf4EBSk66z2HC/+kI8N941VdNB/f89PY2x1nqrj011LN3KycKtpGZqBwA3F5LRPxfK1wFX9zLm6vdpE2BuXr742BX4aNWy1ftbef0zodeYX3w/1teYP6Q4sZHnbK33oSexmZn1Oyc6ZmYbej3wUER8KSJuiIh/kGpWei3SKFW/jYjPAONJF7D7FRbZteoi8nX579yIeIzUrKm4PMAbgXsi4sluNr+GNCBB9Wv/EBE/i4jZETGf59Yc1HpdtTsL6yt6Q2HexupJjD0xB5hQVVY93ZVKP51K/5yKPxfKZxYSoI2NeSapxuixiJhf9Vjaw1ghHfPdtOFAFtsAr+DZ9+NBYFAur9ijekU9OGd7aw6p71iRh6g2s6ZwomNmtqF5wFZK91R5qaQPkQYF6BVJJ0j6gKTdJO0E/DupFuAfhcUC+IWkV+aRtk4HLouIu/P8U4FPS/qopDGSPkbq+P2tHoRwD/AmSdsVLoDnAftLepOkV0j6BrB3jdftkuMeLmmz6hVHxD9JtV4/ljRR0i6SfkAabrknHf7r6UmMPfFd4L2Sjs/H7hig7ohvBVeS+hi9gw0TnWuBHXJ5cVjpjY35PNLxvkzSW5VGqttb0kmSDuthrJAGkFhOGqBgjzwAwAWkZmSVwSRuBB4nDWYxRtJBwJeLK+nhOdtb3wMmSPpaPjZvJ9WAgWt6zKyfOdExMyuIiN8B3yQlE7eTOlefsBGregz4HKk/wu2ki+N3RsS8wjI3kvqNXAHMIP36fkxh/hmki9GTSb+M/z9gckT8rAfb/zzwWtKFdKUPxddJfXwuzXG9mOeOJvcz4CZS/5blwPu6WP9Hcsy/IjX1mgAcGhF39SC2enoSY7ci4rekY3AiqW/NB0jHryeuBZ4mjWL2TL+eXMs2izRK2p/6GnNEPEWqLZlJ6h/zD+AS0uhs9/YwViJiFWkgi9U59mtIzeQOiog1eZlHSO/lPqTj8Z+kY1PUk3O2VyJiFunYfyCv8ySgMkz5Uxu7XjOznlCEf1AxM2s2SeeQRv96S9mxmDVTriX9OfCSiFhRcjhm1sYGd7+ImZmZ2caR9AVS/6ZHSCPN/RfwGyc5ZtbfnOiYmZlZf3oVqRnhlqR7Q/2KdKNeM7N+5aZrZmZmZmbWdjwYgZmZmZmZtZ0B3XRt+PDhMXr06LLDMDMzMzOzAWrWrFkPRcRW1eUDOtEZPXo0M2fOLDsMMzMzMzMboCTVHJK/IU3XJB0kaZ6k+ZIm15i/v6RHJd2SH1+utR4zMzMzM7NG6HONjqRBpLt5HwgsBm6SND0i5lQt+peIOLSv2zMzMzMzM+tOI2p09gLmR8SCfAfmC4BJDVivmZmZmZnZRmlEojOSNC5+xeJcVm1fSbdK+r2k3RqwXTMzMzMzs5oaMRiBapRV35znZmDHiFgp6RBgGjCm5sqkY4FjAUaNGtWA8MzMzMzMrNM0okZnMbBDYXp7YGlxgYh4LCJW5ueXA0MkDa+1soiYGhHjI2L8Vls9Z5Q4MzMzMzOzbjUi0bkJGCNpJ0mbAkcA04sLSNpWkvLzvfJ2H27Ats3MzMzMzJ6jz03XImKtpOOAGcAg4OyIuFPSx/P8M4F3AZ+QtBZYBRwREdXN28zMzMzMzBpCAznfGD9+fPiGoWZmZmZm1hVJsyJifHV5Q24YamZmZmZmNpA40TEzMzMzs7bjRMfMzMzMzNqOEx0zMzMzM2s7TnTMzMzMzKztONExMzMzM7O240THzMzMzMzajhMdMzMzMzNrO050zMzMzMys7TjRMTMzMzOztuNEx8zMzMzM2s7gsgMwMzMzM7OBb9rsJZw2Yx5LV6xiu2FDOWHizhw2bmTZYXXJiU43Wu0NNTMzMzNrtGmzl3DSJbez6ul1ACxZsYqTLrkdYMBeGzvRqaMV31Cz/uCE38zMrLOdNmPeM9fEFaueXsdpM+YN2GsC99Gpo94batYpKgn/khWrCJ5N+KfNXlJ2aGZmZtYkS1es6lX5QOAanTpa8Q01a7RW/AXHzBrPNbtmnf052G7YUJbUuAbebtjQEqLpGSc6dbTiG9ofOvlDDd5/J/xm5qbcZv4cnDBx5w32H2DokEGcMHHnEqOqryFN1yQdJGmepPmSJteYL0k/zPNvk7RHI7bb306YuDNDhwzaoGygv6GN1unNljp9/6HrxL7TEn7rbNNmL2HClKvYafJlTJhyVUd9B4CbcpuBPweHjRvJqYfvzqaDUvowcthQTj189wGd5PU50ZE0CDgdOBgYC7xP0tiqxQ4GxuTHscAZfd1uM7TiG9ponf6h7vT9Byf8Zv7BwzW7ZuDPAaRr43GjhrH3Tlty/eQ3D/hr4kbU6OwFzI+IBRGxBrgAmFS1zCTgF5H8HRgmaUQDtt3vWu0NbbRO/1B3+v6DE34z/+Dhml0z8OegFTWij85IYFFhejGwdw+WGQksa8D2+91BV/+abZcv4t7rXlh2KE333/etYPXadc8p32zwIO498sISImquTt//inHAD5Y9BsDYES+EO+HeckNqqodWrmbRI6tYvXYdmw0exA5bDmX4FpuVHZY1yWcWPNzlvHvvPKd5gZTohytXs+ChJ1i/Pp4p22QT8dLhm3PvkeeUF1iT+bugs/lzkBy97DHu32oH+Ni+ZYfSrUYkOqpRFhuxTFpQOpbUvI1Ro0b1LbIGOXj3Eaye+2jZYZRihy2H1vxQ77BlZ/x60en7XzR2ROcl+pAubIrnwOq161jw0BMAHXeBM6eY7HaQzQYP6vIHj05ROdc7+SLf3wXpGPxz+RNEREeeA/4cJGNHvJBxu7ZEw6yGJDqLgR0K09sDSzdiGQAiYiowFWD8+PE1k6Fm2/bkk8sOoTQ7AotqjDr22g5pttTp+2/w/ilX1Rx9ceSwoVw/+c0lRFSeE3/yNwAubIFf8Rpp9uwlnFJjpKFTD9+dHTvou2BH4LVlB1GiTv8uqB5xDJ79HHRSU+ZO/xy0mkYkOjcBYyTtBCwBjgDeX7XMdOA4SReQmrU9GhEt0WzNUh+NTvoSq9bp+9/p3E/LKp//Th5m3vxd4HuqWSvqc6ITEWslHQfMAAYBZ0fEnZI+nuefCVwOHALMB54Ejunrds3MmsH300qmzV7C7PtWsGbdeiZMuarjLvT9g4d1+ndBpyd61poach+diLg8Il4RES+LiG/msjNzkkMebe1Tef7uETGzEds1M+tvHl772SYra9atBzpzeGWzTv8u8Ihj1ooakuiYmbWryvDaI4cNRXTm8NoeXtnM3wWdnuhZa2pEHx0zs7bW6c2W3GTFLOnk7wL3VbNW5ETHzMzq6vS+CWaWdHKiZ63JTdfMzKwuN1kxM7NW5BodMzOry01WzMysFTnRMTOzbrnJipmZtRo3XTMzMzMzs7bjRMfMzMzMzNqOEx0zMzMzM2s7TnTMzMzMzKztONExMzMzM7O240THzMzMzMzajhMdMzMzMzNrO050zMzMzMys7TjRMTMzMzOztuNEx8zMzMzM2o4THTMzMzMzazuD+/JiSVsCFwKjgYXAeyLiXzWWWwg8DqwD1kbE+L5s18zMzMzMrJ6+1uhMBq6MiDHAlXm6K2+KiNc4yTEzMzMzs/7W10RnEnBufn4ucFgf12dmZmZmZtZnfU10tomIZQD579ZdLBfAHyXNknRsH7dpZmZmZmZWV7d9dCT9Cdi2xqwv9mI7EyJiqaStgSsk3RUR13axvWOBYwFGjRrVi02YmZmZmZkl3SY6EfGWruZJekDSiIhYJmkE8GAX61ia/z4o6bfAXkDNRCcipgJTAcaPHx/d74KZmZmZmdmG+tp0bTpwVH5+FHBp9QKSNpf0gspz4K3AHX3crpmZmZmZWZf6muhMAQ6UdDdwYJ5G0naSLs/LbANcJ+lW4Ebgsoj4Qx+3a2ZmZmZm1qU+3UcnIh4GDqhRvhQ4JD9fALy6L9sxMzMzMzPrDUUM3G4wkpYD95YdRzYceKjsIKxUPgfM54D5HDDweWA+BwaaHSNiq+rCAZ3oDCSSZvpmp53N54D5HDCfAwY+D8znQKvoax8dMzMzMzOzAceJjpmZmZmZtR0nOj03tewArHQ+B8zngPkcMPB5YD4HWoL76JiZmZmZWdtxjY6ZmZmZmbUdJzo9IOkgSfMkzZc0uex4rLkk7SDpz5LmSrpT0vFlx2TlkDRI0mxJvys7Fms+ScMkXSzprvx9sG/ZMVlzSfqP/H/gDknnS3pe2TFZ/5J0tqQHJd1RKNtS0hWS7s5/X1xmjNY1JzrdkDQIOB04GBgLvE/S2HKjsiZbC3w+InYF9gE+5XOgYx0PzC07CCvND4A/RMQupBth+1zoIJJGAp8BxkfEK4FBwBHlRmVNcA5wUFXZZODKiBgDXJmnbQByotO9vYD5EbEgItYAFwCTSo7JmigilkXEzfn546SLm5HlRmXNJml74G3AWWXHYs0n6YXAG4GfAUTEmohYUWpQVobBwFBJg4HnA0tLjsf6WURcCzxSVTwJODc/Pxc4rJkxWc850eneSGBRYXoxvsjtWJJGA+OAG0oOxZrv+8CJwPqS47ByvBRYDvw8N188S9LmZQdlzRMRS4DvAPcBy4BHI+KP5UZlJdkmIpZB+jEU2LrkeKwLTnS6pxplHqquA0naAvhf4LMR8VjZ8VjzSDoUeDAiZpUdi5VmMLAHcEZEjAOewM1VOkruhzEJ2AnYDthc0gfLjcrM6nGi073FwA6F6e1xVXXHkTSElOScFxGXlB2PNd0E4O2SFpKar75Z0q/KDcmabDGwOCIqtbkXkxIf6xxvAe6JiOUR8TRwCfC6kmOycjwgaQRA/vtgyfFYF5zodO8mYIyknSRtSup4OL3kmKyJJInULn9uRHyv7His+SLipIjYPiJGk74DrooI/5LbQSLifmCRpJ1z0QHAnBJDsua7D9hH0vPz/4UD8IAUnWo6cFR+fhRwaYmxWB2Dyw5goIuItZKOA2aQRlg5OyLuLDksa64JwJHA7ZJuyWUnR8Tl5YVkZiX4NHBe/tFrAXBMyfFYE0XEDZIuBm4mjcY5G5hablTW3ySdD+wPDJe0GPgKMAW4SNKHSQnwu8uL0OpRhLubmJmZmZlZe3HTNTMzMzMzaztOdMzMzMzMrO040TEzMzMzs7bjRMfMzMzMzNqOEx0zMzMzM2s7TnTMzDqApP0lhaTtS9r+1ZLOauL2Pi/pd83aXl9IGiRpjqRDy47FzKyd+D46ZmYtTlJ39wm4F3gF0BF38Ja0JfCfwJvKjqUnImKdpFOA70q6PCLWlx2TmVk7cI2OmVnrG1F4TMplexXK9oyINRFxf4dcRH8Y+EdEzC47EEmbSBrUg0V/C7wEOKSfQzIz6xhOdMzMWlxOYO6PiPuBR3Lx8kL58uqma4XpQyT9TdIqSbMk7ZYf10l6UtKNksYWtyfptZL+KGmlpOWSLpG0Y29ilnRgbs72iKRHJV0jaa+qZXbK23lK0n2SPtXDJnAfAKYV1vMmSesk7VC1/qMkPS7pBXl6G0nn5H16XNL1kt5YWF6Sfirpn/l4LZD0LUmbFZY5RdJ8Se+VdBewBtg1H9MZklZIekLSXElHVl4XEU8DlwEf7M1xNDOzrjnRMTPrbN8Evgi8lnRRfj5wBvCVQtnPKwvnpOca4G/AeODNwDrgCknP68V2twBOB/YBXgfcDfxB0kvydkSq5XgR8Ebg7cDbgHH1VirpxcCrgBsrZRHx57z+f69a/CPABRHxuKShwJ+BFwAH5+1cnvdr18rqgQeA9wO7Ap8FjgFOrlrvdsAngaOBsaSmg+cDD+d93R34HPCvqtfdQIs0tzMzawXuo2Nm1tm+GhFXAUj6HnAR8K6IuDKXfRe4RNIWEbESOBH4XUR8pbICSR8kXbQfRKEmpZ6I+G1xWtKxwDvzOs4D3gK8GhgTEfML21nczap3JCUkS6rKpwLHS/p6RKyXtDPwelLCAfBe4IXAeyNibS77pqQDgI8Bn83N/r5UWOdCSS8jJTVfKZQ/DzgyIu4r7N+OwPciYk4uWlAj9sXA1pI2j4gnutlPMzPrhhMdM7POdmvh+f357201yrYGVgJ7Ai+XtLJqPc8DxvR0o5J2Ar4G7JvXvQnwfFKiAqkm5KFKkgMQEY9ImtfNqofmv09VlZ9Dqr2aCPwe+Chwa0TclOfvCWwLrEiVSc/YDFhViPujpJqg0cDmpP+j1a0jHigmOdl3gLMkHQ1cDUyPiJurlqnEPBRwomNm1kdOdMzMOtvThedRp2yTwt9fAlNqrOvhXmz3d8BDwKeARaQmctcBm9bYdm8sz3+3BO55ZkUpSboY+KikPwEfAk4pvG4TYC7wjhrrfBJA0rtJze0mk5rvPQa8m5RAFT0nSYmIr0s6j1Rj9WbgZEnfjohiDdGWpGaAj1S/3szMes+JjpmZ9cZMUh+Yf0bExiQi5H44Y4FDImJGLtueVLNTMQfYStLLC03XXkwaJntWndUvAFYAu9VY7iekfjgfJ9XGnFe1Xx8CHouIrobgfiMwOyK+V9iX0XVi2UBELAB+DPxY0mTgBDZsCrd7Xn8njIxnZtbvPBiBmZn1xrdIHfF/JWmvPDLamyT9QNJLe7iOf5FqXj4q6RWS9iV11l9VWOZPpGZ1v5C0p6RXk2qS1lKnpicnCTOA/WrMuw6YR2pGdlFEPFqYfR6pBugySW+VNFrS3pJOknRYXmYesLukSZJeJul44PDudlbSFpJOl/TmfLzGkWp25lQtuj9p5DUzM2sAJzpmZtZjETGXNHLYFqSEYg7wU1K/khU9XMd6UpOvl5H6A50DfB9YVlgmSM3IngD+Qmrq9ntSslHd/6baGcC78khq1X5Kah43tSqmp0jJ0UzSKHP/AC4h3Y/o3rzYT0jJ1s+B2cDebNj8rStrgRcDPyM1j5vBs6O3AZCTxL3yMmZm1gDayJYHZmZmTZXvd7MY+FJE/E83y/6JNDrc96vKvw0cHBG791ugG0HSj0n/kz9RdixmZu3CfXTMzGxAkvR2Um3IXFL/na+Qmq1d1IOXfxI4sLCuF5H6wHwU+I+GB9sHkjYhJXBTu1vWzMx6zjU6ZmY2IEk6AvgyaSjnJ0iDC3whIu7YiHVdTWpqdiHw7+7wb2bW/pzomJmZmZlZ2/FgBGZmZmZm1nac6JiZmZmZWdtxomNmZmZmZm3HiY6ZmZmZmbUdJzpmZmZmZtZ2nOiYmZmZmVnb+f8DPyIDhF235AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x504 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "year_range = year_end-year_start-1\n",
    "max_lag = 12\n",
    "#corr_ = {}\n",
    "\n",
    "fig = plt.figure()\n",
    "fsize = 7\n",
    "plt.gcf().set_size_inches(fsize*2,fsize)\n",
    "plt.suptitle( 'Correlation of local employment with housing prices' ,fontsize='xx-large')\n",
    "\n",
    "#NAICS_codes_sig =  [ '--', '23','44-45','48-49','52','54','62','72' ]\n",
    "#NAICS_codes_sig =  [ '--', '23','72','44-45','48-49','52','54' ]\n",
    "NAICS_codes_sig =  [ '--', '23','72','48-49' ]\n",
    "\n",
    "for i,NAICS_code in enumerate(NAICS_codes_sig):\n",
    "    lags = []\n",
    "    corr_ = []\n",
    "    corr_auto = []\n",
    "    corr_house = []\n",
    "    corr_poverty = []\n",
    "\n",
    "\n",
    "\n",
    "    for lag in range( max_lag ):\n",
    "        emp_relppc_s = []\n",
    "        emp_relppc_s_shifted = []\n",
    "        house_relppc_s = []\n",
    "        poverty_relppc_s = []\n",
    "        emp_us_ppc_s = []\n",
    "\n",
    "        for MSA_rank in available_MSAs:\n",
    "            MSA = msa_emp_rank_df_dict[ houseRank_to_empRank[MSA_rank] ]\n",
    "\n",
    "\n",
    "            if NAICS_code in list(tot_empl_all[MSA].columns):\n",
    "                ppc_relative_state, ppc_relative_country =  ppc_relative(MSA,NAICS_code)\n",
    "                x=ppc_relative_country.values[:year_range-lag]\n",
    "                #x=ppc_relative_state.values[:year_range-lag]\n",
    "                emp_relppc_s += list(x)\n",
    "\n",
    "                x_shifted = ppc_relative_country.values[lag:year_range]\n",
    "                #x_shifted = ppc_relative_state.values[lag:year_range]\n",
    "                emp_relppc_s_shifted += list(x_shifted)\n",
    "\n",
    "                y1=house_ppc_relative(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "                house_relppc_s += list(y1)\n",
    "\n",
    "                y2=relative_ppc_poverty(MSA_rank).loc[ list( range(year_start+1,year_end) ) ].values[lag:year_range]\n",
    "                poverty_relppc_s += list(y2)\n",
    "\n",
    "                #y3 = tot_empl_state_ppc[ ( 'United States' ,NAICS_code) ].values[lag:year_range]\n",
    "                #emp_us_ppc_s += list(y3)\n",
    "\n",
    "        #pcorr = pearsonr(emp_relppc_s,emp_relppc_s_shifted)\n",
    "        pcorr = pearsonr(emp_relppc_s,house_relppc_s)\n",
    "        #pcorr = pearsonr(emp_relppc_s,emp_us_ppc_s)\n",
    "        if len( emp_relppc_s )>40:\n",
    "            lags += [lag]\n",
    "            corr_ += [pcorr[0] ]\n",
    "            \n",
    "    if len(corr_)>10:\n",
    "        plt.subplot(  len(NAICS_codes_sig ) ,1,i+1  )\n",
    "        plt.stem( lags, corr_ )\n",
    "        #plt.stem( lags, corr_house )\n",
    "        #plt.xlabel( 'Time lag (years)' ,fontsize='x-large')\n",
    "        #plt.ylabel( 'Correlation' ,fontsize='x-large')\n",
    "        plt.title( NAICS_dict[NAICS_code] ,fontsize='x-large')\n",
    "        plt.ylim( [-0.6,0.6] )\n",
    "        plt.subplots_adjust(hspace = 0.7)\n",
    "\n",
    "plt.xlabel( 'Time lag (years)' ,fontsize='x-large')\n",
    "\n",
    "\n",
    "plt.savefig('data/plots/corr_emp_house_2.png')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}