From fad52f0628f7ad3f9d7943914d19d6438f67f39e Mon Sep 17 00:00:00 2001 From: JulienPascal Date: Sat, 12 Jan 2019 12:30:11 +0100 Subject: [PATCH] adding tests for the inference part --- .../exampleLinearModel-checkpoint.ipynb | 92 +++++--- examples/exampleLinearModel.ipynb | 92 +++++--- src/econometrics.jl | 4 +- test/runtests.jl | 202 ++++++++++++++++++ 4 files changed, 321 insertions(+), 69 deletions(-) diff --git a/examples/.ipynb_checkpoints/exampleLinearModel-checkpoint.ipynb b/examples/.ipynb_checkpoints/exampleLinearModel-checkpoint.ipynb index 7dc15f2..9311033 100644 --- a/examples/.ipynb_checkpoints/exampleLinearModel-checkpoint.ipynb +++ b/examples/.ipynb_checkpoints/exampleLinearModel-checkpoint.ipynb @@ -316,9 +316,7 @@ { "cell_type": "code", "execution_count": 9, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# x[1] corresponds to the intercept\n", @@ -372,7 +370,7 @@ " output[\"mean_x1y\"] = mean(simX[startT:nbDraws,1] .* y[startT:nbDraws])\n", " output[\"mean_x2y\"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws])\n", " output[\"mean_x1y^2\"] = mean((simX[startT:nbDraws,1] .* y[startT:nbDraws]).^2)\n", - " output[\"mean_x2y^2\"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)\n", + " output[\"mean_x2y^2\"] = mean((simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)\n", "\n", " return output\n", "end" @@ -463,7 +461,7 @@ "name": "stdout", "output_type": "stream", "text": [ - " 67.708742 seconds (236.31 M allocations: 50.348 GiB, 52.94% gc time)\n" + " 87.049598 seconds (236.31 M allocations: 50.348 GiB, 54.97% gc time)\n" ] }, { @@ -650,9 +648,9 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 16, "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [ { @@ -664,7 +662,7 @@ " -4.75213 -1.2148 3.25637" ] }, - "execution_count": 27, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -679,32 +677,41 @@ "metadata": {}, "source": [ "Once the asymptotic variance has been calculated, a summary table can be obtained using the\n", - "function `summary_table`" + "function `summary_table`. This function has four inputs:\n", + "* a SMMProblem\n", + "* the minimizer of the objective function\n", + "* the length of the empirical sample\n", + "* the confidence level associated to the test H0: $\\theta_i = 0$, H1: $\\theta_i != 0$" ] }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 47, "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [ { - "ename": "LoadError", - "evalue": "\u001b[91mArgumentError: Length of iterable does not match DataFrame column count.\u001b[39m", - "output_type": "error", - "traceback": [ - "\u001b[91mArgumentError: Length of iterable does not match DataFrame column count.\u001b[39m", - "", - "Stacktrace:", - " [1] \u001b[1mpush!\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::DataFrames.DataFrame, ::Array{Float64,1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/julien/.julia/v0.6/DataFrames/src/dataframe/dataframe.jl:1028\u001b[22m\u001b[22m", - " [2] \u001b[1msummary_table\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::SMM.SMMProblem, ::Array{Float64,1}, ::Int64, ::Float64\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/julien/.julia/v0.6/SMM/src/econometrics.jl:154\u001b[22m\u001b[22m", - " [3] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:522\u001b[22m\u001b[22m" - ] + "data": { + "text/html": [ + "
EstimateStdErrortValuepValueConfIntervalLowerConfIntervalUpper
10.5876030.018627431.54510.00.5569640.618243
20.581410.0061311294.82930.00.5713250.591494
30.7626680.00570646133.650.00.7532810.772054
" + ], + "text/plain": [ + "3×6 DataFrames.DataFrame. Omitted printing of 1 columns\n", + "│ Row │ Estimate │ StdError │ tValue │ pValue │ ConfIntervalLower │\n", + "├─────┼──────────┼────────────┼─────────┼────────┼───────────────────┤\n", + "│ 1 │ 0.587603 │ 0.0186274 │ 31.5451 │ 0.0 │ 0.556964 │\n", + "│ 2 │ 0.58141 │ 0.00613112 │ 94.8293 │ 0.0 │ 0.571325 │\n", + "│ 3 │ 0.762668 │ 0.00570646 │ 133.65 │ 0.0 │ 0.753281 │" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "summary_table(myProblem, minimizer, T, alpha)" + "df = summary_table(myProblem, minimizer, T, 0.05)" ] }, { @@ -723,7 +730,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 18, "metadata": { "collapsed": true }, @@ -735,11 +742,28 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "StatsModels.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}\n", + "\n", + "Formula: y ~ 1 + x1 + x2\n", + "\n", + "Coefficients:\n", + " Estimate Std.Error t value Pr(>|t|)\n", + "(Intercept) 0.571592 0.00836568 68.3259 <1e-99\n", + "x1 0.588239 0.00218787 268.863 <1e-99\n", + "x2 0.767633 0.00219593 349.572 <1e-99\n" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ols = lm(@formula(y ~ x1 + x2), data)" ] @@ -762,7 +786,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 20, "metadata": { "scrolled": true }, @@ -816,7 +840,7 @@ "name": "stdout", "output_type": "stream", "text": [ - " 41.372656 seconds (72.39 M allocations: 13.920 GiB, 30.75% gc time)\n" + " 55.815914 seconds (72.39 M allocations: 13.916 GiB, 29.16% gc time)\n" ] }, { @@ -828,7 +852,7 @@ " Plot{Plots.PyPlotBackend() n=1}" ] }, - "execution_count": 25, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -847,7 +871,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 21, "metadata": { "scrolled": false }, @@ -856,7 +880,7 @@ "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVWX+wPHPc85VBBGU3BgRTBZNMwFFy13LMsslGcvKJq1QK9sbbXPStnGammmb0tKxxZ+2qJVZlpZLOGopaplmLEqA4RYqIgrec57fH4yMCCI7d/m+X6/z0nvvOfc+Xy98/Z7nPOd5lNZaI4QQQgghaoxR3w0QQgghhPA0UmCJOrFmzRqUUkyfPr1a7zN9+nSUUqxZs6ZG2iWE8Ew1lXOEqCopsIQQQohKateuHe3ataux93v11VcZP348l1xyCQ6HQ04kPYCjvhsghBBCeLt7770XgODgYFq0aMG+ffvquUWiuqQHSwghhKhny5YtIzs7m99++40RI0bUd3NEDZACS1RZYWEhr776KldddRVt27bFx8eHli1bMmrUKLZu3Vqh9zjdzX7kyBEmTpxI69atadSoETExMSxcuLDcYxcsWEB0dDS+vr4EBwdz3333ceLEiRpvoxDCva1bt44BAwbQpEkTmjZtSnx8PKmpqaX2O3DgAA888AARERH4+PjQvHlz4uPj+emnn4r3SU9PRynFr7/+yq+//opSqng7Pd6rKnnnmmuuoXXr1rUSv6gfcolQVFlOTg73338/ffv2ZejQoTRr1ozdu3ezdOlSli9fzrfffktcXNx536ewsJArrriCvLw8brnlFo4fP86HH37ITTfdxKFDh7jnnntKHfPaa6/x5ZdfMmLECAYNGsSXX37JK6+8wqFDh/i///u/Gm+jEMI9bdy4kb/+9a8MGTKEe+65hx07dvDxxx+TmJjIxo0bad++PQBpaWkMGDCArKwsrrzySkaOHMmBAwdYvHgxX331Fd988w09e/akadOmPPnkk7z00ksA3H///cWfNWDAAEDyjvgvLUQVnTx5UmdlZZV6/qefftL+/v76iiuuKH5u9erVGtBPPvlkiX3DwsI0oPv166cLCgqKn8/MzNTNmzfXPj4+JT7jySef1IAODAzUu3btKn4+Pz9fR0VFacMw9N69e6vURiGE5zidcwA9a9asEq/NmjVLA/raa68tfq5Xr17aNE395Zdfltj3l19+0U2aNNFdunQp8XxYWJgOCwsr87Orm3cmTpyoAb169epy9xOuTS4Riirz8fGhTZs2pZ7v3LkzAwcO5Ntvv+XUqVMVeq/nnnuOhg0bFj8OCQnhvvvuo6CggPfff7/U/vfddx8dOnQofuzr68uNN96IbdskJSXVShuFEO4nKiqKhISEEs8lJCQQGRnJ559/zsGDB9m6dSvr16/n1ltv5aqrrirz+O3bt5e4VFgeyTsC5BKhqKZt27bx/PPPs27dOvbt21cqaRw6dIjg4OBy38PhcHDZZZeVer5v374AZY5Z6NatW6nnQkJCADhy5EiNt1EI4Z569+6NYZTsSzAMg969e5OSksIPP/xASkoKAPv37y9z3qxdu3YV/3nxxRdX6HMl7wgpsESVrV+/nkGDBgFw5ZVXEhkZib+/P0opPvnkE3744QcKCgrO+z7NmzcvlQABWrVqBcDRo0dLvRYQEFDqOYej6MfZsqwab6MQwj2dziPnev7o0aPk5OQA8Pnnn/P555+f872OHz9eoc+UvCNACixRDc8++ywFBQUkJibSp0+fEq9t3LiRH374oULvc+jQIWzbLlVk7d+/H4DAwMB6b6MQwj2dziPnej4wMLD4hO3VV19l8uTJ1f5MyTsCZJoGUQ1paWkEBQWVSiD5+fls2bKlwu/jdDrZsGFDqecTExMBiImJqfc2CiHc03/+8x9s2y7xnG3brF+/HqUUXbt2pWfPngBl5qFzMU2zRG/5mSTvCJACS1RDWFgYhw8fZseOHcXPWZbFww8/zMGDByv1Xo899hiFhYXFj7Oysnj55Zfx8fFhzJgxLtFGIYT7SU5O5q233irx3FtvvUVycjLXXHMNLVq0oEePHvTs2ZOFCxfywQcflHoP27ZZu3ZtieeCgoI4dOgQJ0+eLLW/5B0BcolQVMM999zDihUr6NOnD9dffz2NGjVizZo17N27lwEDBlR4Ha3g4GCOHz/OJZdcwrBhw4rnwfr999955ZVXyrwbp67bKIRwT1dddRX33nsvX3zxBZ07d2bHjh189tlnNG/enJdffrl4v4ULFzJw4EDGjBnDSy+9RGxsLL6+vmRkZLBhwwYOHjxYopgaNGgQmzdv5uqrr6Zv3740bNiQfv360a9fvyrlnZkzZxYPpj/dkzZz5kzefvttAEaOHMnIkSNr7x9K1Lz6nidCuLdFixbp2NhY7efnp5s3b66vv/56nZaWpm+99VYN6D179mity58HKywsTOfk5OgJEyboVq1aaR8fH921a1e9YMGCUp93eh6ssuaHmTdvngb0vHnzqtRGIYTnODPnJCYm6v79++vGjRvrgIAAfd111+mUlJRSx+Tk5OgnnnhCX3zxxdrX11f7+/vryMhIfdNNN+klS5aU2PfYsWM6ISFBBwcHa9M0S+W3yuad/v37F8/bVdZ2du4Urk9prXU91HVCABSvRp+enl6v7RBCCCFqkozBEkIIIYSoYVJgCSGEEELUMCmwhBBCCCFqmIzBEkIIIYSoYdKDJYQQQghRw6TAEkIIIYSoYS5bYJ1eUiA/P7++myKEECVIfhJCnI/LFli7du2iW7duxTPbnu3o0aN13KL6540xg8TtTdwlZslPZfPGuL0xZpC4K8JlC6zzOdcim57MG2MGidubeErMnhJHZXlj3N4YM0jcFeG2BZYQQgghhKuSAksIIYQQooa5ZYH1yPcWf93RsL6bIYQQpWTmaQZ87cfuXJliUAhv5pYFVotG8NIvDUk9KglMCOFaLmgEvxcqHtzonWNUhBBF3LLAuquTQXMfzZNJksCEEK7Fz6F4+pICPv1VszzTru/mCCHqiVsWWL4OxZ8vKmRhmubH36UXSwjhWka0cTIwWHHfeosCS3KUEN7ILQssgJvbnSI8AB7fLL1YQgjXohS82stk9zF4abv0Ygnhjdy2wGpgwFPdTJZlaNbvlwQmhHAtnYMU93Q2eHqrzd7j0oslhLdx2wIL4IZwxSVB8NgmG60lgQkhXMv0bgaNG8Cfv5OediG8jVsXWIZSPNvdZG22ZuVeKbCEEK4lsKHibz1MFqZp1mZLT7sQ3sStCyyAa0IVvVop6cUSQrikP0UqLm2puGe9hdOWHCWEt3D7AkspxXNxBkmHNIv3SPISwlulpKTQq1cvoqKiiIuLY8eOHaX2WbVqFT169KBTp0507tyZKVOmYNv/61latmwZHTt2JDIyklGjRpGbm1vtdhlK8Vovk59y4I2d0oslhLdw+wILoH+wwVUhimlJcoYohLeaOHEiEyZMIDk5malTpzJu3LhS+zRr1oz333+fnTt3kpSUxPr163n33XcByMvL4/bbb+eTTz4hJSWFP/zhDzz99NM10rZuLRQJHQ2mJdkcOCE5Sghv4BEFFsCz3U12HYH3UiR5CeFtDhw4wObNmxk7diwA8fHxZGZmkpqaWmK/mJgY2rdvD0CjRo2Ijo4mPT0dgOXLlxMTE0PHjh0BuOuuu1i4cGGNtfHZOANDwWObZMC7EN7AYwqsbi0Uf7xQMX2LTOwnhLfJzMwkODgYh8MBFA0dCA0NJSMj45zH7Nu3j0WLFnHttdcCkJGRQVhYWPHr7dq1Izs7G6fTec73yMvLIzc3t3grKCg4577NGyme6W4w9xfN9wfkUqEQns5R3w2oSU93N+m8yMnsn23uvdis7+YIIVxUbm4uw4YNY8qUKXTv3r3K79O/f/8Sj6dMmcLUqVM5fPhwmfuPbgWzAv2Y9K3FioH5GKrKH+2SzhW3J/PGmEHiLktQUFCJxx5VYHVsqhgXqXhmq81tHQz8G3hY9hJClKlt27bFvU0OhwOtNRkZGYSGhpba99ixYwwZMoQRI0bw4IMPFj8fGhrKypUrix+np6eX6BUry9q1a4mOji5+7OPjg4+PD1A62Z72Rj+bPp9ZfHqwKbd39JiLCMXOFbcn88aYwbvitrXm6S02o1opwoOaVeiYUr/dFbkTB2Du3LlERkYSHh5OQkICp06dAsq/Syc9PR3TNImOji7e0tLSqhpvmZ7sZnK0UJanEMKbtGzZktjYWObPnw/A4sWLCQkJISIiosR+eXl5DBkyhCFDhvDEE0+UeG3IkCFs2bKFXbt2AfD6668zZsyYcj/X39+fgICA4u10cVWe3q0NxkYoHtlkcbhAhjMI4Q4S92mmb7H57UTFO25KFVgVuRNnz549TJs2jcTERFJTU9m/fz9vvvkmUP5dOgBNmjRh27ZtxVt4eHgVQj23UH/FnZ0M/v6jTc5JSV5CeIvZs2cze/ZsoqKimDlzJvPmzQPgjjvuYOnSpQC8/PLLfP/99yxZsqT4JO/ZZ58FinLTnDlzGDlyJBEREWRlZTFt2rRaaevzPU1OWvBkkpwICuEO3k62CQ+A7kEV/50t0fd9+k6cFStWAEV34kyePJnU1NQSZ4KLFi1i+PDhtG7dGoBJkybx3HPPcffddxMTE1O839l36dSVx6IN5uyy+dsPNn/rKWOxhPAGHTp0YMOGDaWenzNnTvHfH3/8cR5//PFzvsfw4cMZPnx4rbTvTMF+iidjDaZ+b3NHB4NLLpDhDEK4qrxTmo92ax7paqAq8ataogeronfilHW3TVl365x9lw7A8ePHiYuLIzY2lqeeegrLKv+W5crcpXNaS1/FA10MXtlh85sssiqEcEH3djaICoR71luyCoUQLmzRHk2+E/4UVbkxk7U2yL2su3SCg4PZu3cvLVu2JCcnhxtuuIEXX3yRKVOmnPN9KnuXzmm3t4V/7fDniY35vBBz/qLMHchdG97FG+M+X8yeNKi2oal45TKTK5dbvJ+muTFCerGEcEVvJ9tc3kYR6q/Iyan4cSUKrIreiRMaGlpicHp6enqJfc51l46Pjw8tW7YEihLlbbfdxoIFC8otsKpylw5AEPBotMVjmxSPx/kRHuAZycuT/oOpDInbe3hTzINDDEa1s3n4O4thYUrufBbCxezO1azN1swfWPnhRiX6uyp6J058fDxLly5l3759aK2ZNWtW8d025d2lc+DAgeK7DQsKCliyZEmJMVtlqcpdOqfd3dmghS88mSQzJwshXNM/LjXJKYBntsqAdyFczTspNgEN4Lp2lT/5KXVBsSJ34rRv354ZM2bQu3dvIiIiaNGiBRMnTgTKv0tn3bp1xMTE0LVrV2JjY2ndunW5A06ry8+h+EuMwYJUzfYcGeMghHA9YU0Uj0Yb/GO7zc+HJU8J4SpsrXkn2eaGcIWfo/IFltIuOrpyy5YtdOvWjaSkJGJjY0u9npOTU6FLCYWW5qKPnFwcpPj0SveeV7WiMXsaidt7uEvMNZWfTjvp1HRZ7KRNY8Xqa0xUZW5VciHu8v3VJG+MGbwj7tW/2Qz63OI/w016tSrqj6pM3J43jfBZGpqKp7qbLP1Vs2G/dMELIVxPI4fijT4ma7M178iC9UK4hHm/2EQGwGUtq3bC4/EFFsCN4YouQfDn72y5HVoI4ZKuaGNwc4Ti4Y0Wh2SSZCHqVW6hZtEezbgoo8o9yl5RYBlK8Y9LTf6zX/PhbklcQgjX9GJPE0vDn7+TG3OEqE+L9mhOWvCnyKqXSV5RYEHR2eHwMMWU7y1OOKXIEkK4nlZ+iud7mrydrFmbLUMahKgv836xGdxGEeJf9fGQXlNgAbzQ0yQ7H174URKXEMI13d5B0buVYtI6iwJLTgaFqGupRzXr9hddHqwOryqwIgMV911sMPMHm72yhI4QwgUZSjGrj0nqUXj+BzkZFKKuvZNiE9gQRlZh7qszeVWBBfBEjEFjBzzyvYxxEEK4pouDFA9dYvDsNpuUo3IyKERdseyiua/GtDfwrcLcV2fyugIrsKHi2TiT+amajTJtgxDCRf0l1iDYF+76jywGLURdWZ2tyTwO46KqPxed1xVYALdFKboGwf0bbWxJXEIIF+TnULzex+TrvZoFaZKnhKgLbyfbdAiEnlWc++pMXllgmYbi5V4m3x3QLEiVxCWEcE1XtzUYfaHiwY0WhwskVwlRm44Wahbv0YyvxtxXZ/LKAgugf7BB/IWKRzZZHD8liUsI4ZpeuszkpBMe+V6GNAhRmz7crSm0YWw15r46k9cWWAB/72Fy6CT8Te7UEUK4qD80VjwXZ/DmLpv1Mm5UiFrzdrLNlW0UbRrXzFqgXl1gXRigeLCLwd9/tPn1mPRiCSFc06SLDOJaKCYmWpyyJVcJUdOSj2jW18DcV2fy6gIL4NGuBk0bwlSZtkEI4aJMQzG7j8nOI/APmShZiBr3dopN04YwIqxmeq9ACiyaNFT8tYfJB7s16/ZJ4hJCuKaY5or7OhvM2GKzJ1d6sYSoKZateTfF5sZwg0bVnPvqTF5fYAH8KVLRvbni/g0ybYMQwnU91d2geSO4e73MjSVETfnmN83eGpr76kxSYFG0NMVLlxkkHdK8kyxJSwjhmvwbKF7tZbI8U7Noj+QqIWrCvGSbTk0hrkUtF1gpKSn06tWLqKgo4uLi2LFjR5kHzp07l8jISMLDw0lISODUqVMArFq1ih49etCpUyc6d+7MlClTsO3/XXpbtmwZHTt2JDIyklGjRpGbm1ujAVVV79YGN4YrHttkcaxQEpcQwjWNaGcwIkxx3waLo5KrhKiWIwWaj9OLBrfXxNxXZypVYE2cOJEJEyaQnJzM1KlTGTduXKmD9uzZw7Rp00hMTCQ1NZX9+/fz5ptvAtCsWTPef/99du7cSVJSEuvXr+fdd98FIC8vj9tvv51PPvmElJQU/vCHP/D000/XaEDV8bceJkcL4bltMhZLCOG6Xu1lklsIj2+SXCVEdXyw2+ZUDc59daYS73jgwAE2b97M2LFjAYiPjyczM5PU1NQSBy1atIjhw4fTunVrlFJMmjSJhQsXAhATE0P79u0BaNSoEdHR0aSnpwOwfPlyYmJi6NixIwB33XVX8XGuoK2/YkpXg39st9ktg0iFcCsV6X1PT09nwIABBAYGEh0dXeK1NWvW4OvrS3R0dPF24sSJump+pbT1Vzzd3eD1nTI3lhDVMS9ZMyREEexXs71XcFaBlZmZSXBwMA6HAwClFKGhoWRkZJQ4KCMjg7CwsOLH7dq1K7UPwL59+1i0aBHXXnvtOY/Lzs7G6XSes4F5eXnk5uYWbwUFBVUIs+KmdDVo6Qt//k6mbRDCnVSk9z0gIIBnnnmGBQsWlPkeHTp0YNu2bcWbr69vLbe66u7tbNCjpeK2tRYnnHJCKERl/XxY892BoqVxaoOjVt4VyM3NZdiwYUyZMoXu3btX+X369+9f4vGUKVOYOnUqhw8frm4Tz2laJwcTN/my9Jej9GnhOoVWbcbsyiRu73G+mIOCgsp8/nTv+4oVK4Ci3vfJkyeTmppKREREieP79OnDmjVraqzN9cU0FP/uZxKzxMn0JJu/9TTru0lCuJV3Umya+cCwGpz76kwlCqy2bdsW9yg5HA601mRkZBAaGlrioNDQUNLS0oofp6enl9jn2LFjDBkyhBEjRvDggw+WOG7lypUljjuzx6wsa9euLdGV7+Pjg4+PD3DuZFtdCc00b/9q8Zef/Ei6zoFp1M4/flXUVsyuTuL2HlWJubze9zMLrPNJS0sjNjYW0zQZP348d911V7n7n+5hP+3M/FQXOjVTTO9m8MRmm1EXKnq2lBvDhagI53/nvrop3MDHrIMCq2XLlsTGxjJ//nzGjRvH4sWLCQkJKZWg4uPj6dOnD9OnT6dVq1bMmjWLMWPGAEUJZ8iQIQwZMoQnnniixHFDhgzh7rvvZteuXXTs2JHXX3+9+Lhz8ff3JyAgoCZirTClFC/3MujxicVbu2wmdZIzQyE8XWxsLFlZWQQGBpKVlcXQoUNp3rw5119//TmPqY8e9rPdHgIfpPhx6yqL1Zfn41OP6Up6Xb2Hu8f99T6T7Hw/rmt9jJycio9jLC/us08OS3UdzZ49m3HjxvHcc88REBDAvHnzALjjjjsYPnw4w4cPp3379syYMYPevXsDMGDAACZOnAjAyy+/zPfff8/x48dZsmQJAKNHj+bxxx+nSZMmzJkzh5EjR+J0Orn44ot55513KhxYXYprYXBblM0jm2xGtDNqZQCcEKJmVLT3vTxnnsiFhIRw4403kpiYWG6BVR897GWZf7km9mMnr6UH8mxc/Z4QSq+r93DnuBdvdXJxM82g9oGVnp6honGXKrA6dOjAhg0bSu04Z86cEo8TEhJISEgotd/jjz/O448/fs4PPF2kuYPne5p8luHk/g0WH1xea8PVhBDVVNHe9/JkZ2fTqlUrDMPg2LFjLFu2jNtvv73cY+qjh70sFwcppsUULaMzqp1BtxqeMFEIT7I/X/NJuua5uJqf++pMcsG+HBc0Urx8mcmHuzXLfpVboYVwZbNnz2b27NlERUUxc+bMEr3vS5cuBSA/P5+QkBBGjx7Nzp07CQkJ4dFHHwVg8eLFdOnSha5du3LppZcyePBgxo8fX2/xVNYj0QZdgmD8t04KLbmrUIhzeW2nTQMDbutQuyWQdMucx5hwxbspirv+Y9E/WNGkoZwZCuGKKtL77ufnR1ZWVpnHT548mcmTJ9da+2pbA0Pxdn8H3T928uw2mxndZOyoEGfLd2pe32lzR0eDZj61+/+59GCdh1KK13ub/F4A05KkF0sI4bq6XqB4LMbgua02236XXiwhzvZ2ss2RQrivc+2XP1JgVcCFAYqnuhm88pPN9wekyBJCuK7How0uagbj1zo5ZUuRJcRplq3553abP16ouDCg9q9GSYFVQfddbBDTHCYkWpK0hBAuq6FZdKlwew7MlHVVhSi29FdNai481KVuSh8psCrIYSje6utg+2H453ZJWkII1xXbXDG1q8HTW22258gJoRAAL2y36dta0aOOJuSVAqsSYpsrHrjY4MkkmzRZDFoI4cL+EmsQGQDj11o4pdddeLkN+23W79c8fEndlT1SYFXSjG4GrXxh0joLrSVpCSFck4+pmNffZOvvmr//KL3uwru9uN0mKhCuDa27mQCkwKqkxg0Us/qYfL1XMz9VCiwhhOvq0dLg4S4G05Nsdh6WfCW8U1qu5uN0zYNdDIxanFj0bFJgVcGQtgY3hSse2GBx6KQkLSGE65rRzeDCJnKpUHivl7bbBPnAnyLrtuSRAquK/nmZiQ08uMGq76YIIcQ5NXIUXSrcdFDLDTrC6+Sc1Pw72ebuTga+jrqdKFwKrCpq6at4safJe6malVmStIQQruuyVgYPdDGYlmTzyxHpxRLeY9bPNraGuzrVfbkjBVY1jItSDAxWTFpnke+UpCWEcF1Pdzdo2xhu+9bCkkuFwgsUWJpXd9j8KdKgpW/dL3MnBVY1KKWY3ddkbz48tUV6sYQQrsvPoXi7v8nGA5q/ygSkwgssSNXsOwEP1tHEomeTAquaIgMVf4kxeOFHWftLiOpyyv/7tap3a4PHog2mb7H5Tpb9Eh5Ma82L2y2GhSo6NK373iuQAqtGPHyJwUVNIUG63oWosoc2Wty43re+m+Hx/hJr0L2F4qZVFscKJV8Jz/RVlmbHYep0YtGzSYFVAxqairf6miQd0ry2U84KhaisDftt/rnd5orWzvpuisdrYCj+b6DJgZNwz3q5C1p4phd+tIlroejbun56r6CMAislJYVevXoRFRVFXFwcO3bsKPPAuXPnEhkZSXh4OAkJCZw6dQqA9PR0BgwYQGBgINHR0SWOWbNmDb6+vkRHRxdvJ06cqIWw6t6lrQzu7mTw+CabjDw5KxSiogotzYREi+4tFHeEn6rv5niF8ADFa71M3knRfJAmJ4XCs2z7XfPNb5qHuhioOpxY9GylCqyJEycyYcIEkpOTmTp1KuPGjSt10J49e5g2bRqJiYmkpqayf/9+3nzzTQACAgJ45plnWLBgQZkf2KFDB7Zt21a8+fp6ziWBZ+MMmvnAbWstbFlGR4gKeeFHm5+PwFt9Tcz6y4Ve50+RihvaKyaus+SkUHiUF3+0CPOH+AvrN6GUKLAOHDjA5s2bGTt2LADx8fFkZmaSmppa4qBFixYxfPhwWrdujVKKSZMmsXDhQgCCgoLo06cPjRs3rqMQXEdAw6K7dL75TSb0E6IiUo5qntpq81AXg64XSHVVl5RSvNHHJKABjF0t40eFZ8jK07yfprn/YgOH4UIFVmZmJsHBwTgcDqDoFzA0NJSMjIwSB2VkZBAWFlb8uF27dqX2OZe0tDRiY2OJi4vj9ddfP+/+eXl55ObmFm8FBQUV+pz6cnkbg4e6GDy2Se4qFKI8WmsmrbP4gx882U2Gg9aHZj6K+QNN/rNf87cf5KRQuL9Xdtg0bgC3d6j/nOKoyw+LjY0lKyuLwMBAsrKyGDp0KM2bN+f6668/5zH9+/cv8XjKlClMnTqVw4cP13Zzq+yhcPgqw48bVlqsujwfX7Nm3teVY65NErdnWvirg1W/+bKoTz4ncy1Ocv6Yg4KC6qZxXqRfsMGjXTVPJtlc0UbRo2X9/8ckRFXkFmpm/2xzZyeDJg3rv0e8RIHVtm1bsrOzcTqdOBwOtNZkZGQQGhpa4qDQ0FDS0tKKH6enp5fapywBAQHFfw8JCeHGG28kMTGx3AJr7dq1JQbL+/j44OPjA7h2sv1gsKbbx05mJgfyau8aqrBw7Zhrk8TtWQ6e0Pxlu5ObIxTxFwWWeM1TY3ZlT3YzWLlXc9Nqi63XKZf4z0mIypr7i02+E+7p7BonCSVa0bJlS2JjY5k/fz4AixcvJiQkhIiIiBIHxcfHs3TpUvbt24fWmlmzZjFmzJjzflh2dja2XdQNfezYMZYtW0ZMTEy5x/j7+xMQEFC8nS6uXF2nZooXehq8ttPmiwzpehfiTA9tLJoe4B+X1tzJh6i601M37D8B98kC9sINnbI1L/1kc1OD5mgPAAAgAElEQVSEok1j1zhBKFXmzZ49m9mzZxMVFcXMmTOZN28eAHfccQdLly4FoH379syYMYPevXsTERFBixYtmDhxIgD5+fmEhIQwevRodu7cSUhICI8++ihQVLB16dKFrl27cumllzJ48GDGjx9fV7HWubs6GQxtqxj/rcWBEzIeSwiAlVk276VqXuhp1sv6YKJsEYGKV3uZzEvWfLRbTgqFe1m0W5ORBw92cZ2TNqW1a84nsGXLFrp160ZSUhKxsbGlXs/JyXGLSwn78zVdFjvp2VKx9EqzWnNyuEvMNU3i9hz5Tk2XRU5C/RWrrin9++AuMXtKfjqb1pobvrFYuVfzQ7yDUP/K5St3jbs6vDFmcK24tdZ0/8TJBT6KFUNrd2h5ZeJ2jQuVHqyVn+Lf/UyWZWhm/SxnhcK7Pb3FZm8+zO5TvZMNUTtOL2DfpAH8aY1M3SDcw9pszZZD9bssTllcqzUe6towgzsvMnhoo83PhyVhCe/04++av/9o80S0QVQtLL5akVUoyltpAs69QoU3aeajeG+gybfZmud/lJNC4fpe3G7TJQgGt3GtkzYpsOrIC5cahDWBm1c7KbSkyBLexbI1E9ZZdAiEKV1rJ+1UZBWK8laaKG+FCm/TP9jgkWiDv2y22XRQiizhujbst1mWoZlyiev1ikuBVUf8HIoFAx38dBimbZaEJbzLGz/bfHdA82Zfk4a1sB5ORVehKG+lifJWqPBGM7oZxDRX3LTKIu+UnBQK16O15qGNNtEXwI3hrlVcgRRYdSqmueLZ7gZ//9Fm9W9SZAnvkJWneWyTzcSOBr1b107KqegqFOWpygoV7rbSRGWcnrohOx/uXS9TNwjXs3iPZsOBojuSzXpeFqcsdTqTu4CHLjFYnqm5ZY3Fj6MUQY1c74dCiJp0z3qLxg6Y2cPzzufccaWJyrgA+GtXB/cm+RLTJJeb2znL3d9T4q4Mb4wZ6j/uQhv+vLExg1vbxPgeIyenbj63vLjPvrtQCqw6ZijFOwNMui5xMmmdxQeXu951YyFqysd7bD75VfPR5SZNfWrv57yiq1CUpyorVLjrShOVcU8Q/JDn5M/bfOkd6iC2efnfo6fEXRneGDPUb9wvbbfIyLf5/OqGBDXzrdPPlmkaXFhbf8XsPiYf7dG8myJjG4Rnyi3UTF5vcW2oIv7C2j2JqOgqFOWpygoV7rrSRGW91sukSzPFqJVOfj8pOUvUr8MFmqe22iR0MOjUzHU7KKTAqiej2xuMi1JMXm+RlisJS3iexzbZHC2Ef/Wum17aiqxCUd5KE+WtUOHtGjkUi64wyTsFN6+W+bFE/Xpmq80pG6Z3c+0SRi4R1qNXLjP5NtvJ2NUWicNMHC44SE+Iqtiw3+b1nTb/uNSo9GzgVdWhQwc2bNhQ6vk5c+YU/93Pz4+srKxzvkdCQgIJCQm10j53F9ZEsXCQyZAvLaZvsXm6u+ssSSK8x+5czas7bP4SY9Daz7X/z3Tt8s/DNWmomD/QZNNBLVM3CI9xtFBz82qLuBbKZVa1FzVjcIjBM90Nntlq89mvkrNE3Xt0k0WLRvCgi83aXhbXb6GHu6yVwcweBjN/sPkwTRKWcG9aayYkWvx+EhYOcs1bp0X1PNLVYGSY4pY1FqlH5VKhqDsb9tt8uFvzbJyJn8P1c4sUWC7goS4GN4Yrxn9r8cPvkrCE+3prl+bD3Zo5/UzaB7h+AhSVp5Ti7QEmrXzhupVOjsskpKIOaK15+LuiSUVviXCP3CIFlgtQSjGnn0mHQBi5Qu7SEe5pe47mvg0WEzsajG4vqcWTBTZULLnCwZ5jMCHRQmvJWaJ2LUnXrN/vupOKlkWyoIvwcyg+HuwgzwnXf2PhlLt0hBs5fkpzwzdOIgPhn5dJWvEGnYMU/+5nsiCtaNCxELWl0NJM/d5iaFvF5W3cJ7+4T0u9QFgTxUeXm6zN1kz5ThKWcB/3rrf4NQ8+GOTA1w3GRoiacX24wYNdDB7aaLNun+QsUTte32mz5xg838O97lyVAsvFDPiDwUuXGfzzJ5t3kyVhCde3INXm38ma13qZXOTCk/6J2jGzh0GvVorRX1vsOyHfv6hZpycVvaODQecg9/r5kgLLBd3dyWB8lGLCOovNB6XIEq4r5ahm4jqLmyMU46LcK/mJmtHAUHx4uYmh4PbvGnFKhjeIGvTsfycVneHik4qWpVSLU1JS6NWrF1FRUcTFxbFjx44yD5w7dy6RkZGEh4eTkJDAqVOngKL1uwYMGEBgYGCJNbrOd5z4H6UUr/c2ib5Acd1Ki/35krCE6ymwNGNWOQn2hTfqaLZ24Zpa+RXN9J6UY/LwRjkpFDXj9KSiUy9x/UlFy1KqwJo4cSITJkwgOTmZqVOnMm7cuFIH7dmzh2nTppGYmEhqair79+/nzTffBCAgIIBnnnmGBQsWVOo4UVIjh2LxFSZOG/74tUWhJUWWcC1Tv7f5KQc+uNxBk4bul/xEzbqslcEzXQt4ZYfNglQpskT1PbbJormbTCpalhKtPnDgAJs3b2bs2LFA0eKnmZmZpKamljho0aJFDB8+nNatW6OUYtKkSSxcuBAoWmW6T58+NG7cuNSHlXecKK1NY8XiwSbfHdTcv0ESlnAdS3+1efknmxd6GsQ0l+JKFLm9/SluiVAkJFpsz5GTQlF1G/fbfOBGk4qWpUSBlZmZSXBwMA5H0RKFSilCQ0PJyMgocVBGRgZhYWHFj9u1a1dqn7JU5bi8vDxyc3OLt4KCgvNH5UF6tTJ4vbfJGz/bvLunQX03Rwgy8zTj11qMDFNMlqVwxBmUgll9TSIDiub0O3hCiixReVprHvrOpmuQ+0wqWhaXX+y5f//+JR5PmTKFqVOncvjw4XpqUd0b1RLWt/dhylYfOgYcoccF3tWb5U3f9ZlcMW6nDaO/9cXPNHjhkuPUdBPPF3NQUFDNfqCocX4OxSdXOrj0UyfDVlisusZ9eyBE/Tg9qejKoe4zqWhZShRYbdu2JTs7G6fTicPhQGtNRkYGoaGhJQ4KDQ0lLS2t+HF6enqpfcpSlePWrl1bYrC8j48PPj4+gHcl21kDNLtyTzL+u8Zsvs5Bm8bu+0NXFd70XZ/J1eJ+YpPF5hybb4eZhLeqnba5Wsyi8to1USy7ymTAMosbV1ksucK9/6MUdefMSUWvcKNJRctSovUtW7YkNjaW+fPnA7B48WJCQkKIiIgocVB8fDxLly5l3759aK2ZNWsWY8aMOe+HVeU4f39/AgICirfTxZW3aWgq/t3zJA4DRq20OOmUrndRt77ea/PcNpunuxv0auXeiU/Uvu4tDD683OTzDM09621ZTkdUyBs/u+ekomUplSVnz57N7NmziYqKYubMmcybNw+AO+64g6VLlwLQvn17ZsyYQe/evYmIiKBFixZMnDgRgPz8fEJCQhg9ejQ7d+4kJCSERx999LzHifNr2UjzyWAHP+Zo7vyPrP8l6s7+fM3Y1RZXtFFM7SrFlaiYoaEGs/oUjSF9/gfvGtogKi8zTzM9yeb2DsrtJhUtS6kxWB06dGDDhg2ldpwzZ06JxwkJCSQkJJTaz8/Pj6ysrHN+4LmOExXTrYXirb4mt6yxuKipzZSu7l/lC9dma82f1lgAvDfAxJD5rkQl3NHRIPO45pFNNiH+ipsjpEAXpVl2UZ7xbwB/84DeK3CDQe6itLGRBr8c1Uz93qZpQ8WEiyRhidozPclm5V7NiqEmrdxwsj9R/6bHGmT89+7TYF8Y5OZja0TNe3G7zdpszTfXmDTz8Yw8IwWWm3qqm8HRQpi0ziKgIYwJl4Qlat6cXTZPb7X5a5zh9gNORf1RSvFmX5Pfjltct9Ji3XBFFw+4BCRqxpZDmic22/z5EoOBf/CcPOM5kXgZpRQvXWZwS6TiltUWn2fI+AZRs77IsJm0zuKuToaMuxLV1sAoWk6nfQBc/aWTrDwZQyog36m5ebWTi5vB0909K894VjRexlCKuf1Mrg1V/PFrizW/SZElasamgzajv7G4NlTxymWGrDMoakSThoovhjgwVVGRdbRQiixv9+fvbNKPwf8NdNDQ9Kw8IwWWm3MYivcvN+nTWjFshcWmg1JkiepJy9Vc86XFJUGKBYNk/iJRs4L9FF8OcZB1HK5bIeuserNlv9q8vtPmxUsNLmrmeXlGCiwP4GMqPhls0qWZYshyix2yBpioooMnNFcvd9K0IXx2lczALWrHRc0US680+c9+zW3fypQz3mh/ftF3f01bxZ0eeqOWZ0blhRo3UHw+xKRtYxi83MnuXElYonLynZphKyyOnoIvr3bQvJEUV6L29A02eG+Ayf+lah7fLD3v3kTrouIKYG4/02OHIEiB5UGa+Si+utqBvwOu+MLJ3uNSZImKsWzNjasstudoPr/KpH2AZyY84VquDzd4safBX7fZvLHTqu/miDryxs82X2Rq/t3Ps6d+kQLLw7TyU3x9jQOnDVd+4eTQSSmyRPm0LlrK5PMMzUeXm3RvIWlB1J0Huhjcd7HB5PU2S3+VnixP9/NhzUMbbe68yODaMM/ONZ4dnZcK9Vd8PdTBwZMwZLlFrtypI8ox8webN362md3HZGio+6aElJQUevXqRVRUFHFxcezYsaPM/ebOnUtkZCTh4eEkJCRw6tQpANasWYOvry/R0dHF24kTJ+oyBK+klOLFngbXtVNc/43FFzLljMcqtIqmZGjXBF641H1zTUV5foReKqqpYsVQB6m5mmFfWZyQxaFFGd5LsXlsk82TsQa3d3TvdDBx4kQmTJhAcnIyU6dOZdy4caX22bNnD9OmTSMxMZHU1FT279/Pm2++Wfx6hw4d2LZtW/Hm6+tbhxF4L9NQ/N9AkyEhipErLT5NlyLLE03bbPPTYVgw0OEVN9C4d0YV5Yq+QPHFEJPNhzR//FpuhxYlfb3X5ra1FrdFKZ6Mde9UcODAATZv3szYsWMBiI+PJzMzk9TU1BL7LVq0iOHDh9O6dWuUUkyaNImFCxfWR5PFWXxMxUdXmIwIK5rX76PdUmR5kjW/2fz9R5tnuhvENPf84gqkwPJ4vVoZfDLY5Ou9RQtpWrYUWQK2/a4ZtdLiijaKWX3d/y6ezMxMgoODcTiKVv9SShEaGkpGRkaJ/TIyMggLCyt+3K5duxL7pKWlERsbS1xcHK+//vp5PzcvL4/c3NziraCgoIYi8k4NDMXCQSbXt1eMWWWxIFWKLE9wuEBzyxqL/sGKh7p4T9khaxF6gcEhBgsHwehvLGxt8e4Ak0Ze0D0rypaRpxn6pZPIQPjoCpMGMpEoALGxsWRlZREYGEhWVhZDhw6lefPmXH/99ec8pn///iUeT5kyhalTp3L48OHabq5Lqqm4X+oK2tmIsas1Obknuamds0betzbId10+reGO7xtxrNDBK9F5HD3i3if55cUdFBRU4rEUWF5i1IUGiy6Hm1ZbXLXc4pMrPWfFclFxOSeLJhL1MeDzqxz4N/CMn4G2bduSnZ2N0+nE4XCgtSYjI4PQ0NAS+4WGhpKWllb8OD09vXifgICA4udDQkK48cYbSUxMLLfAWrt2LdHR0cWPfXx88PHxAUonW29RU3HPH6wJWGdzT5IvDX1NJrjwZJTyXZ/bu8k2n2RZvD/IpEtIszpoVe2r6Pftuj+xosZdd6HBN9eY/HRY03upkwxZbNWr7D2u6b/Myf4TRROJtvag+WdatmxJbGws8+fPB2Dx4sWEhIQQERFRYr/4+HiWLl3Kvn370Foza9YsxowZA0B2dja2XXRJ6tixYyxbtoyYmJhyP9ff35+AgIDi7XRxJarPUIo3+hjc09lg4jqL13bIPFnuZk+uZvJ6i1siFDeEe1+54X0Re7lerQzWD3dwwgmXfurkh9+lyPIGvxwpKqqPFsK64Q46NPWc4uq02bNnM3v2bKKiopg5cybz5s0D4I477mDp0qUAtG/fnhkzZtC7d28iIiJo0aIFEydOBIqKsi5dutC1a1cuvfRSBg8ezPjx4+stHlE0lu7lywwevsTgnvU2L/4oRZa7OH5KM2aVxQU+8Fpvs76bUy9KXSJMSUnh1ltv5dChQwQGBvL222/TuXPnUgfOnTuXmTNnYts2gwYN4vXXX6dBgwblvrZmzRquvvpqOnToUPw+GzZskFuh61iHpooNIxxc85WTvp85WTLY5Io2Umt7qu8P2Az90qKVL3x1tYMQf88rrqBoioUNGzaUen7OnDklHickJJCQkFBqv8mTJzN58uRaa5+oGqUUz/cw8DHg4e9sCix4LMY7/8N2FyedmpErLXYe0ay6xiSgoWfmnPMp9b9qdeeSkXlm3ENrP8Waaxz0aqW4ernF/BS5W8cTfZVpM+hziw5NFYnDPLe4Ep5NKcUzcSZPdTN4fLPNk0myQLSrOmVrblhlsW6fZtlVJnFevDJEichrYi4ZmWfGfTRpqPjsKpNbIhW3rLGYuU2SlidZkGpz7VcWA/+gWDnUJEgWbxZublqsycw4g6e22Dy6yZZ85WIsW3PrGovlmZolg036B3tvcQVnXSIsby6ZMweLljeXTEXnmTFNk/Hjx3PXXXeV28DT88ycduZdOqL6GhiKuf1M2voXJayMPHi1l4Ept+67tZd/srh/g824KMWbfWUqBuE5pkabNDThwY1Flwv/canh9vO4eQKtNZPWWXywW/Ph5SZXt/Xu4grqeJoGmWememoz5vsuhGY04OGtPqQfKeDNHifxc5FJPLzxu4aqxa01PL2jIS//4sO9UQX85eJCjh2phcbVkvPF7K23w4uSHuhi4mPC3f+xyXfCa70NOYmoR1oXLeA85xfN2/1N4i+U4grOKrBqYi4ZmWemdtVmzPcHQWQLm+u/UYze0JDPrjJp7iKXlbzxu4bKxe20NRMSLeYla/5xqcEDXfxrsWW1x1u/a1E5d3Uy8TUVExItdh0p6jVp5UFTj7iT6Vts/vmTzb96G9waJcXVaSX+JWpiLhmZZ8a9XRNqsOZak7RcTa9PnezOlTEO7iDfWbT0zXspmvkDTR7oIndZCc83voPB6mtNko9qun3i5PsDcrNOXXvhR4unttjMjDO4q5PknTOVKjWrO5eMzDPj/uJaGGwY4UADly118m22JC1XdrhAc+UXFt/8pvnsKpObI+QMUniPPq0Nkq5zENpY0fczi7m7JF/VlVk7Lf78nc3j0QZTo6W4OpvSLnobxpYtW+jWrRtJSUnExsaWej0nJ8frLiXUdcwHT2j++LVF4j7N1K4GM7oZNDTrvgveG79rqFjce49rrlruZF8+fD7EpGdL9y6u3OW7lvxUtvqMu8DS3L/BZtbPNhM7Grzcy8CnDvKVt37Xs7blctcmX+7pbPDSZd5zo0Flvm/3zsaiVrXwVay6xuS5OIMXfrS59FMnPx92yXrcK237XdNrqZNjp4pmZ3f34kqI6vAxFW/0MZnT12Ress2AZRZ7j0u+qg0f77GZvLkR46MU//Si4qqyJCOLcpmG4pFok+9GOjhhQezHTv61Q+bLqk+WrfnrNosenzgJ8oH1wx109MClb4Soits7GiQOM8k8run2sZN1++SSYU36KtPmhlUWw9s4ebOviSHF1TlJgSUqJLa5Iuk6B7d3MJi83uaaryz25UuRVdfScjX9llk8vsnmwS4GG0c4aNNYEpwQZ+rR0iBppIMOgYqBy4oWipaTwur7NtvmupUWV4Uo3og7KfMlnocUWKLC/ByK13qbfH6VyZZDmi6LnXyaLmeHdUFrzVu7bLoudpKdr/l2mMnMHmadjDERwh218lN8fY3J3Z2LFooet9bihFOKrKradLBoZYjLWik+utykgVQP5yX/RKLShoYabI930LuVYuRKi4RvneSdksRVW/bla4avsJiQaHFjuOKHUQ76tJZfXSHOp4GheOkyk/cGmHy0W9PnMye/HpNcVRlaa/79S9GYts7NFJ9eadLIISd2FSFZWlRJC1/Fx4NN3uprsiBNE73Eycb90ptV0z7eY9NlsZPvD2o+vdLkrX4OmnjpyvRCVNXYSIP1IxzknIRuHztZkSW5qiIOF2jGrLK4/duik7uVQ038G0j+qSgpsESVKaW4o6PBtlEOLvBR9PnMYkaShdOWM8TqOlqouXtzI0Z9bdGnleKneAfDw+TXVYiqir5Asfk6B7HNFVctt7hplZPf5C7Dc0rMtole4mRFVtEs+XP6OaS4qiTJ2KLaIgMV64abPBFj8NRWm76fWWw9JImrqtb8ZnPJYifL9jqY199kyWCTFr6S2ISorgsaKb682mRef5NvftN0+MjJ33+wKLQkX53mtDVPJlkM+NwizF/xQ7yD0e2lVKgK+VcTNaKBoZjezWTdMJPDBZrYj5388WsnO3IkcVXUSafm4Y0Wgz63aOevSLziOOOiZI4ZIWqSoRTjogx+Ge3gtiiDRzbZdF3i5Ou9ctlwT66m32cWz261mR5rsPoak1B/yT9VJQWWqFGXtTL46Y8O/t3PZPPBojsNx652knpUCq1z0VqzIssm7hMnr+6web6nwaprTEIby7+ZELWlqY/i5V4mW69z0KKRYvAXFqO/dpKZ552/dwtTiy4Jnr5LeVqsKdMwVJMUWKLGOQzF+A4Gydc7+Fdvg9W/aTp+5CThWycZXpq8ymJrzSfpNj0+sbhquYWfo2iMyMOXSGIToq5ccoFi7bUm8wearNtXlKue22pR4CWXDY8Vam5d4+Sm1RbXhCq2xTvo1UpKg5og/4qi1jQ0FXd2Mkm9wcHfexp8+qsm8gMn9/zHItuLJyl12poFqUXjrK5badG4Aay42mTjCJMuQVJYCVHXlFLcHGHwy/UO7rzI4Mkkm4sXOVme6dmXDb8/YBPzsZMl6Zp3B5j830CTQLlLucZIgSVqna9D8UAXk91jHEzvZjA/1Sb8fSdTvrM4dNJ7Cq1CSzNnl03Hj5zcvNoi1F+xbpjJmmsdDA6RsVZC1LeAhooXLjX5Id5BqL9i6JcWI1c42ZPrWXnq9HJbvZdaXOCj2DbKwS2RkoNqmhRYos74N1A8Gm2yZ4yDhy8xeONnmwvfd/KXzRZHCjwrgZ0p36l55SeL8A+cTEi0iL6gaNmhL4Y46C0Thgrhcjo1U3w91OTDy02SDmk6LSo6IUw+4t556oRT826yTc9Pi5bbmtLVYN1wk/AAKaxqg6O+GyC8T1MfxVPdTe692OD5H2xe+NHm1R02o9sr4tsZDGqjaOABY5ByCzX/2mnzz+02OQVwc4Tika4mFzVz/9iE8HRKKUa3V1zdVvHXbTav77T5+482vVopxkcZXN9eEeAml9NSjmpm/WzzdnJRLhrcRrHmWpN+wXKCV5ukwBL1pnkjxfM9TR7sYvDyTzYf7bF5a5dF04YwIkwRf6HB4DbukcBO01qz8zC8v7uoaDzhhNs6GEy5xOBCOUsUwu34N1A8G2cyLaZoHOm8ZJsJiRb3rof4CxXxwSbDm2kMF7u85rQ1S3/VvPGzzdd7NUE+MD7KYOJFBpGBrtVWTyUFlqh3rf0Uf+1h8lycwY85sGiPzeI9Nu+kWPg3gCtbNeKmjjZDQhSNXXAm4QMnNF/v1azca7MiS/NbPvg5YNJFBg91MfhDY9drsxCicho5FDeEK24IN8jK07yXajPvF5v5qX6EbXVya5TBrZEG7ev5RGrv8aKF4d/aZfNbPlzWUvHuAJPRFypZQ7COleofTElJoVevXkRFRREXF8eOHTvKPHDu3LlERkYSHh5OQkICp06dqvZrwrsppeh6geLp7iY7Rzdgxx8dTLnEIPmYwR+/tmjxnpP4lU4WpNrkFtbfWIiTTs03e20e+d4idskpWs0vGrSedEhzY7jBV1ebHLrFwYuXmlJc1aHazl1CnBbiXzSe9JfrHSwfcJzBbRT/3G4T/oGTAcucvJNsk3eq7nKU/d+59K5b4SRsoZMXfrQZFlq0jNn6EUUD2KW4qgf6LAMHDtTz5s3TWmv90Ucf6e7du5+9i969e7cODg7W2dnZ2rZtPWzYMP3aa69V67WzJSUlaUAnJSWV+frvv/9e5vOezBtj1roo7pQjtp651anjPj6lebNQN5xTqK9efkpP2ejUs3Y69VeZlk4+YusCp13jn2/btt7+u63/8aNTD/nilPadW6h5s1C3eq9Qj111Sr+bbOnfjtf853rj912dmGszd51N8lPZvDHu0zHnFdr6vWRLD1pWlKMa/7tQD/vylH5wg1P/a4dTf5lh6ZQjti60qpcr8k/Z+ufDtl6eYek3djj1nzc6dcT7RTnp4o8K9b92OPXRgprPR2fzxu9a68rFXeIS4YEDB9i8eTMrVqwAID4+nsmTJ5OamkpERETxfosWLWL48OG0bt0agEmTJvHcc89x9913V/k1IcoTEaiYGm0yNRp+PaZZkm6zPFPz4W6bzONwek5AQ0EbP2gfoLiwCVzYRNG+yf/+3toPCizIKYDDBZBToMkpOP247L9nHdfsPwE+JvRrrXiqu8GVbQy6BCG3NbuI2s5dQpxP4waKsZGKsZEG6cc076bYrN+vWfqrTfoxcP43R5kKwvwhPED9d4PwJv/7u6Hg1zxIP6aLttN//++f+0/87zNNBW0bQ+/Wirf7G/RqpSQnuZASBVZmZibBwcE4HEVPK6UIDQ0lIyOjRJLKyMggLCys+HG7du3IyMio1mvnkpeXR25ubvFjHx8fAP72t7/x1FNPFT/2dAUFBV4XM5Qdd1iTonm1HuhStM8pW5OZB3uOafYcg93HNHuOaX4+Ap9n2Bw8+b/3MxTY5+i5b9IAgnyKtmY+iiCfoqLs6raKfq0VfVorfOuom90bv+/qxFzbuetcJD/9j/zM/i/mdk0Uf4k1ix87/5uj0o5pUo9q0o5BWq5m/X6b91LguLPs9zcVhPpDO39Fp6ZwdYhBuyaKdk2KnmvTuGjljLrmjd81VD5ulx/k3r9//xKPp0yZwp133pUd0VYAACAASURBVMnzzz/PnXfeSUBAQD21rG7l5uZ6XcxQ8bibAjG+RRstS76W54TM4wbpxxX7Tho0dmiaNdA0bVi0NWsIgQ00Dc5zx/KJXDhR/i41xhu/74rEHBQUVMetKp/kp/+Rn9nyYw4EYn2LNlr/73mt4WCBYk+ewZ7jCg2E+WlCG9u0bqRxnCsvnYLcIzUUSCV543cN54/77PxUosBq27Yt2dnZOJ1OHA4HWmsyMjIIDQ0tcVBoaChpaWnFj9PT04v3qepr57J27Vqio6OLH/v4+FBQUFAcjLd8uafPzL0pZqiZuIOAUKB3zTWr1nnj912dmGs7d52L5Kf/kZ/Zqsd8AdCxhtpUF7zxu4bKx12iNm7ZsiWxsbHMnz8fgMWLFxMSElKiix2KxjcsXbqUffv2obVm1qxZjBkzplqvnYu/vz8BAQHFmzd1RwohKqa2c9e5SH4SQpxLqUuEs2fPZty4cTz33HMEBAQwb948AO644w6GDx/O8OHDad++PTNmzKB376I+gQEDBjBx4kSAKr92thMnii7G/Pzzz6Vey8vLA2Dbtm34+/tXPXo34o0xg8TtTXFXNOaOHTvi5+dX6vnazF1nk/xUmjfG7Y0xg8RdXtwl8lNt3MZYE+bPn68B2WSTTbYS27mmRpD8JJtsstX3dmZ+UlprjQs6dOgQX331Fe3atcPX17e+myOEcBHn6sGqS5KfhBBlOTM/uWyBJYQQQgjhrmQpbSGEEEKIGiYFlhBCCCFEDXOpAstbF2utSNyrVq2iR48edOrUic6dOzNlyhRs2waK5usxTZPo6Oji7cz5fFxRRWJes2YNvr6+JeI6ffcWeO53PW/evBIxN2/enFGjRgHu+V3fe++9tGvXDqUU27ZtO+d+7vB77Y05yhvzE3hnjvLG/AS1mKPq6SacMtXlYq2upCJxb9myRf9/e/ceF1W57w/886w1OF5wUEKMBCRULDBBFDVRQUslb2Sal45t2abSxbMvnn6w2zuPu7Z5zP3r/OrYaYflce+ynSZoorRTSyRNzbsZxuYSApqIhgooF2et5/fH5BwRxREGmMvn/XrNq2bNmjXP19U8flrrmecpKCiQUkpZXV0to6Ojre8pLCyUXl5erdVcu7Cl5szMTBkeHn7L97vyub5ZWFiYTE1NlVI657nOysqSJSUlsmfPnvLo0aO33MdZvtfu2Ee5Y/8kpXv2Ue7YP0nZcn2UwwSsc+fOyc6dO8tr165JKaXUdV12795d5uXl1dtvxYoVMjEx0fo8IyNDRkdH3/E1R2Vr3Td74YUX5JIlS6SUzvcfta01N9Z5ucu53r9/v+zWrZusq6uTUjrfub5RY52XM3yv3bGPcsf+SUr37KPcvX+S0v59lMPcImxssdYb2Xux1rZma903Ki0tRWpqKiZOnGjdduXKFURFRSEyMhKvvvoqNE1r8bY31d3UXFBQgMjISERFReGdd96xbneXc7169Wo8/fTT8PDwsG5zpnNtK2f4XrtjH+WO/RPgnn0U+6fGNeV77fCLPVN9FRUVmDRpEpKSkjBo0CAAgJ+fH86cOQNfX1+Ul5djxowZeOONN5CUlNTGrW2eyMhInD59Gl5eXjh9+jTGjx8PHx8fTJ8+va2b1iquXLmCdevWYf/+/dZtrnquyTW4U/8EuHcfxf7pzhzmCtaNi7UCaHSx1qKiIuvzmxdrvd1rjsrWugGgsrIScXFxiI+Px6JFi6zbjUYjfH19AVgWoZw7dy52797dOgU0ga01m0wmeHl5AQD8/f0xa9Ysa12ufq4BYMOGDQgLC0NoaKh1m7Oda1s5w/faHfsod+yfAPfso9g/Na5J3+uWuZPZNDExMfUG2A0cOLDBPgUFBQ0Gk61cufKOrzkyW+qurKyUw4YNk6+88kqD186dO2e9B15TUyOnTZsmFy9e3KJtbi5bav7xxx+lpmlSSikrKirksGHD5OrVq6WUrn2urxs+fLh877336m1zxnN9XWPjG5zle+2OfZQ79k9Sumcf5c79k5T276McKmDl5OTIoUOHyj59+siBAwfKb7/9Vkop5TPPPCM3b95s3W/VqlUyODhYBgcHy7lz51pP6J1ec1S21L106VJpMBhkeHi49bF06VIppZRpaWkyLCxM9u/fX4aGhsqFCxfKmpqaNqvHFrbUvHLlShkaGmqta8mSJVLXdesxXPVcX9/P09NTVlRU1Hu/M57rBQsWyB49ekhVVaWvr6/s1auXlNI5v9fu2Ee5Y/8kpXv2Ue7YP0nZcn0Ul8ohIiIisjOHGYNFRERE5CoYsIiIiIjsjAGLiIiIyM4YsIiIiIjsjAGLiIiIyM4YsKjJdu3aBSEE/vjHP7Z1U4iIiBwKAxY5rKCgIAQFBdnlWHl5eVi2bBlGjhyJ++67D+3atUNAQAB+8YtfICcnxy6fQUREdB3XIiS3sHjxYqxfvx79+vVDfHw8TCYTTpw4gQ8//BCpqan4/PPPMXLkyLZuJhERuQgGLHILcXFxSE5OxoABA+ptX7duHWbNmoXnnnsO2dnZbdQ6IiJyNbxFSHaxZ88exMbGonPnzujSpQumTp2K/Pz8BvuVlZXht7/9LXr37g2j0QgfHx9MnToV3333nXWfU6dOQQiBoqIiFBUVQQhhfVwf71VXV4eVK1di3LhxCAgIsC4y+sQTT+Do0aMNPjchIaFBuAKAmTNnIiQkBCdPnsSFCxfs9wdCRERujUvlUJPt2rULo0aNwrhx45CZmYm4uDiEhYUhOzsbW7ZsgY+PD/bv34/g4GAAQEFBAWJjY3H69GmMHTsWDz30EMrKypCWlgYhBL788ksMGTIEly5dwptvvok333wTAPCb3/zG+pmxsbGIjY1FaWkpevTogREjRqBv377o2rUrfvjhB6Snp0MIga+++gpRUVE21dGvXz9kZ2fj4sWL6NKli/3/oIiIyO0wYFGTXQ9YAPDuu+8iMTHR+lpKSgqeffZZTJw4EVu2bAEAREdH45tvvkFGRgbGjRtn3Tc3NxeDBg1CUFAQvv32W+v26wPcT5061eCza2trceHCBfTo0aPe9uzsbAwdOhRDhw7Fjh077ljDgQMHMGTIEERFReHAgQM2105ERNQY3iKkZgsJCcH8+fPrbZs/fz769OmDjIwMnD9/HkePHsXevXsxZ86ceuHqxvefOHGi3q3CxhiNxgbhCgDCwsIwatQofPXVV7h27Vqjx7h8+TLmzJkDRVGwYsUKmz6XiIjIFhzkTs0WHR0NRamf1RVFQXR0NPLy8nD8+HHk5eUBAM6dO3fLebOuT5WQk5ODfv362fS5x44dw4oVK7Bnzx6UlpY2CFQXLlyAn5/fLd9bXV2NKVOmICcnB6+99hpiY2Nt+kwiIiJbMGBRs3Xv3r3R7ZcvX0Z5eTkAICMjAxkZGbc91pUrV2z6zL1792L06NEAgLFjx6JPnz7w9PSEEAKffvopjh8/jtra2lu+t6amBvHx8cjMzMRLL72E3//+9zZ9JhERka0YsKjZzp071+h2Ly8vmEwmAMDKlSuxcOHCZn/ma6+9htraWuzevRvDhw+v99r+/ftx/PjxW76vuroa8fHx2LFjB5KSkrBs2bJmt4WIiOhmHINFzfb1119D1/V623Rdx969eyGEQHh4OIYMGQIA2Ldvn83HVVUVmqbd8rWCggJ4e3s3CFdXr17FkSNHbvmeG8PViy++iNdff93mthAREd0NBixqttzcXLz33nv1tr333nvIzc3FhAkT0K1bNwwePBhDhgzBxx9/jPXr1zc4hq7ryMrKqrfN29sbFy5cQE1NTYP9e/bsiYsXL9abHFTTNLz44os4f/58g/2v3xbcsWMHFi1ahD//+c9NLZeIiOiOOE0DNdmd5sG655578M0331jnwSosLMSoUaNQVFSEoUOHIjIyEh06dEBxcTH27duH8+fP1wtTycnJWLFiBWJjYzFixAi0a9cOI0eOxMiRI7F161ZMmjQJXbp0wfTp09G+fXvs2rULZ86cwUMPPYRdu3ahsLDQOtVDQkIC/va3v+Hee++tN53EjRISEuy29iEREbk5SdREmZmZEoBcsmSJ3L17t4yJiZGdOnWSJpNJTpkyRebl5TV4T3l5uXz55Zdlv379ZIcOHaSnp6fs06ePfOqpp+TGjRvr7VtZWSnnz58v/fz8pKqq1s+6LjU1VUZGRsqOHTtKHx8fOX36dFlQUCDnzJkjAcjCwkLrvjExMRJAo4/MzMwW+pMiIiJ3wytYRERERHbGMVhEREREdsaARURERGRnDFhEREREdsaARURERGRnDFhEREREdsaARURERGRnDhuwri95cvXq1bZuChEREdFdcdiAlZOTg4EDByInJ+eWr1++fLmVW9T23LFmgHW7E3esmYhck8MGrDu53SLArswdawZYtztxx5qJyDU5bcAiIiIiclQMWERERER25pQB63cHNCQfM4LLKBIREZEjcsqAFXmPwPsF7fCfJ/S2bgoRERFRA04ZsKb3UvDrvrVIOqBjWwlDFhERETkWpwxYAPCHsDo85i8wY6eG3Eu8VUhERESOw2kDliqAj0ar8OsITN5uxuU6hiwiIiJyDE4bsADAq53A5jEGlFYDT+3UoOkMWURERNT2nDpgAUBIF4H1o1V8flriD4c4HouIiIjantMHLAAYF6BgxWAFrx/X8XE+QxYRERG1LUNbN8BeFj2k4PhPEnO/0hDiJTCwm2jrJhEREZGbcokrWAAghMCqESr6ews8vsOM0qscj0VERERtw2UCFgC0NwhsGqNCk8DULzTUagxZRERE1PpcKmABwH2dLCHr0HmJF77WuJwOERERtTqXC1gAMMRXwaoRKlb/U+LtbA56JyIiotblMoPcbzYnxDLo/bf7dYR1FRjdwyWzJBERETkgl04dK4YoGH2fwJNfavihgrcKiYiIqHW4dMAyKALrH1HhbQTit5tRyeV0iIiIqBXYLWDl5eVh2LBhCAkJQVRUFLKzsxvss2/fPkRERCAiIgJhYWFITExEbW2tvZpwS12NAuljDSiqAp7epUHnoHciIiJqYXYLWImJiViwYAFyc3ORnJyMhISEBvuEh4fj4MGDOHbsGE6cOIGysjK888479mrCbT3YVeDj0SrSiyT+cJCD3omIiKhl2SVglZWV4dChQ5g9ezYAYOrUqSgpKUF+fn69/Tp27AgPDw8AQF1dHaqrqyFE68y4PiFQwf8domD5cR0f5DJkERERUcuxS8AqKSmBn58fDAbLjxKFEAgMDERxcXGDfU+dOoXw8HD4+PjAy8sLzz//fKPHrqqqQkVFhfXRnFuKv31IwTN9Bebv1vB1KUMWERERtYxWn6YhKCgIx48fR1VVFWbPno2NGzdi5syZt90/Jiam3vOkpCQkJyfj4sWLTfr8paHA9z91wOPbdewYdRWBnZxnTFZTa3Z2rNt93Klmb2/vVmoJEVHz2CVgBQQE4OzZszCbzTAYDJBSori4GIGBgbd9j6enJ2bOnImPPvqo0YCVlZWFiIgI63Oj0Qij0Qig6Z3t5sckhnxqxtPfeGLvZAM6t3OehaHd9S8Y1u0+3LFmInI9drlF6Ovri8jISKxduxYAkJaWBn9/f/Tu3bvefvn5+bh27RoAyxisTZs2oX///o0e29PTEyaTyfq4Hq6aw6e9wJZxBhRXAU9latB057mKRURERI7Pbr8iTElJQUpKCkJCQrB8+XKsWbMGADBv3jykp6cDAHbu3IkBAwYgPDwcAwYMQPfu3bF48WJ7NeGuhHYVWPeIis9KJF7iLwuJiIjIjoR00NWQjxw5goEDB+Lw4cOIjIxs8Hp5ebldbiW89Z2G3+zT8T8jVfyyr2PPu2qvmp0N63Yf7lgzEbkml12L0Fa/ClNw8iKQuEdDbxMwws+xQxYRERE5PrdPE0IIvB2tILq7wJQdXLOQiIiIms/tAxYAeCgCqY+q6GIEJm0zo4JrFhIREVEzMGD97J72AlvHGXDmKjBzJ39ZSERERE3HgHWDB7oIfPKIiu2nJf7PN/xlIRERETUNA9ZNxvorePNhBf/vOx3v5TBkERER0d1z+18R3soLoQqyLwLP79HQxwTE3sccSkRERLZjcrgFIQT+a5iCGD+BJ77Q8M9LHI9FREREtmPAuo3rvyy8twMw/nMzzlczZBEREZFtGLAa0cUo8FmcAVVmIH67hhozQxYRERHdGQPWHQR1FtgyVsXRnyTmZGnQHXNlISIiInIgDFg2GOyrYO0oFRt+kHj5EH9ZSERERI1jwLLR1PsVrBii4D+O6VjN6RuIiIioEZym4S7820MKCiqAZ/do6NkZeLQH8ykRERE1xIRwF4QQWDlMwaM9BKbu0JBdzvFYRERE1BAD1l0yKALrH1HRszMwYZsZpVcZsoiIiKg+BqwmMLUTyBhnQJ0OTN6u4SqnbyAiIqIbMGA1UYCnwNZxBmRflJidqUHTGbKIiIjIggGrGSJ9BNaNVrG5SCL5AH9ZSERERBYMWM00qaeCN4cqeOOEjr+c1Nq6OUREROQAOE2DHfxrPxX5FcDCvTp6egqMD2RuJSIicmdMAnbyn0MVTAgQmLFTw7GfOB6LiIjInTFg2YmqCPx9tIoQL2D852YUVTJkERERuSsGLDvy9LBM39BeBeI+N+OnGoYsIiIid8SAZWf3dhTY9pgBF2osc2RVc44sIiIit8OA1QL6eAlkjFNx7CeJWTs1mDlHFhERkVthwGohg30VfPKIiq3FEgu/1iElQxYREZG7YMBqQRMCFawaoSIlR8drRzkRKRERkbvgPFgtbG5fBWeuSCw+rOO+TgJz+zLTEhERuToGrFbw8gAFZ64AC3Zr6N7BcmWLiIiIXBf/pm8FQgj8d7SCiYEC07/UcKCMtwuJiIhcGQNWK7k+EWm4t8CEbRpyL3HQOxERkatiwGpFHQ0CW8ap8GlvmYj03FWGLCIiIlfEgNXK7mkv8HmcATUaMH6bGZV1DFlERESuxm4BKy8vD8OGDUNISAiioqKQnZ3dYJ+dO3di8ODBCA0NRVhYGJKSkqDr7jceqWdngX/EGZB/GZj2hYY6jSGLiIjIldgtYCUmJmLBggXIzc1FcnIyEhISGuzTtWtXrFu3DidPnsThw4exd+9efPDBB/ZqglMJv0dg0xgVmWcl5n2lcSJSIiIiF2KXgFVWVoZDhw5h9uzZAICpU6eipKQE+fn59fYbMGAAgoODAQDt27dHREQETp06ZY8mOKXRPRR8EKviw3yJpAOc7Z2IiMhV2GUerJKSEvj5+cFgsBxOCIHAwEAUFxejd+/et3xPaWkpUlNTsXXr1kaPXVVVhYqKCutzo9EIo9Foj2Y7hJm9FJRVS/x6n46u7YDfD1DbuklERETUTG0y0WhFRQUmTZqEpKQkDBo0qNF9Y2Ji6j1PSkpCcnIyLl682JJNbFWz7wPOhrbDHw4ZYTBfxbxe1265nyvVfDdYt/u4U83e3t6t1BIiouaxS8AKCAjA2bNnYTabYTAYIKVEcXExAgMDG+xbWVmJuLg4xMfHY9GiRXc8dlZWFiIiIqzPb7yC5Uqd7bJhEjWKjuRj7XFfl06Y3efWd29dqea7wbrdhzvWTESuxy5jsHx9fREZGYm1a9cCANLS0uDv79/g9mBVVRXi4uIQFxeHl19+2aZje3p6wmQyWR+udHvwRkIIvDFUwS9DBBKyNKQXud+vK4mIiFyF3X5FmJKSgpSUFISEhGD58uVYs2YNAGDevHlIT08HALz11ls4cOAANm7ciIiICEREROC1116zVxOcniIEVo1Q8XiQZUmdzB8ZsoiIiJyRkA7607UjR45g4MCBOHz4MCIjIxu8Xl5e7rK3Emo1icnbNOwtk/hyvIrBvpYc7Mo1N4Z1uw93rJmIXBNncndARlVg4xgV/b0FHvtcw3flDpmBiYiI6DYYsBxUJw+BjHEqAjoBY/9hxg8VDFlERETOggHLgXUxCmx7zABPA/DoZ2acrRZt3SQiIiKyAQOWg+veUeCLCQaYdWDq7g74qYZXsoiIiBwdA5YTCPQU2DHegJ/qLGOyKusYsoiIiBwZA5aT6NtFYMPwavzzksTk7RqqzQxZREREjooBy4n076IjI07FN2USM77UcE1nyCIiInJEDFhOZvi9CjaOUfH5aYnZmRrMDFlEREQOhwHLCcUFKPjkERUbCyWe3sWQRURE5GgYsJzU40EK1j+iYsMPEnN2adAYsoiIiBwGA5YTe+J+BetGq1j/g0RCFkMWERGRozC0dQOoeaYFK/gYwKydGgQ0rIlRoSqckJSIiKgtMWC5gCeDFegSeCpTgyI0rB7JkEVERNSWGLBcxIxelpA1e5clZL0/UoUiGLKIiIjaAgOWC5nVW4EE8PTPIWvVCIYsIiKitsCA5WKe6m25kvWLn0PWu8MZsoiIiFobA5YLmt1HgSaBX2ZpUKDjneEKQxYREVErYsByUXNCLLcL52ZpUATw39EKBEMWERFRq2DAcmEJIZbbhc98ZQlZK4cxZBEREbUGBiwXN7evJWTN320JWW89zJBFRETU0hiw3MC8BywhK3GPBrMOvB3NMVlEREQtiQHLTSx4UIFBsVzJqrgmsSZGhQcnIyUiImoRDFhuZG5fBZ09gH/J1FB5TcP60SraGxiyiIiI7I2LPbuZJ4MVpI9VseO0xIRtGqqucYFoIiIie2PAckNxAQq2Pabi4HmJRz/TUF7DkEVERGRPDFhuaoSfgswJBhRUSMRmmFF6lSGLiIjIXhiw3NjAbgJfTTTgpxpgxBYziioZsoiIiOyBAcvNPdhVYM8kAySA4VvMyLnEkEVERNRcDFiE+00CuycZ4NUOGLnFjKMXGLKIiIiagwGLAAB+HQWyJhoQ1FkgdqsZe0r1tm4SERGR02LAIqt72gt8OV5FpI/A2M80bCthyCIiImoKBiyqp3M7gc/iVDzSQ2DSdg1phQxZREREd4sBixroYBDYOEbFtPsFpn+pIeV7ra2bRERE5FS4VA7dkoci8GGsCp/2Op7do6OgAlg+mItEExER2YIBi25LVQT+a5iKXibgt/t0/FAp8WGsig5cv5CIiKhRdrtFmJeXh2HDhiEkJARRUVHIzs5usM+pU6cQGxsLLy8vRERE2OujqYX9up+KTWNU/KNEYlSGhnOc9Z2IiKhRdgtYiYmJWLBgAXJzc5GcnIyEhIQG+5hMJixduhR///vf7fWx1ErigxRkTVRRVCkxdLMZ319kyCIiIroduwSssrIyHDp0CLNnzwYATJ06FSUlJcjPz6+3n7e3N4YPH45OnTrZ42OplQ3qpmB/vAGeHsCwdDMyf+QvDImIiG7FLgGrpKQEfn5+MBgsQ7qEEAgMDERxcXGzj11VVYWKigrro7a2ttnHpKbr2Vlgz2QDorpZ5sr6Wy5DFhER0c0cfpB7TExMvedJSUlITk7GxYsX26hFbceRav5wMJB01IiErHb47txV/C60Di31A0NHqrs1uWPdd6rZ29u7lVpCRNQ8dglYAQEBOHv2LMxmMwwGA6SUKC4uRmBgYLOPnZWVVW9AvNFohNFoBOCena0j1fzXRyVCj+v43UEjfrzWHv8To8KotkzKcqS6W5M71u2ONROR67HLLUJfX19ERkZi7dq1AIC0tDT4+/ujd+/ezT62p6cnTCaT9XE9XFHbE0IgOULFJ4+oSDslMeYzDT/VcPA7ERGR3X5FmJKSgpSUFISEhGD58uVYs2YNAGDevHlIT08HAFy9ehX+/v548skncfLkSfj7++Oll16yVxOojTwZrCBzgorvL0k8vNmM/MsMWURE5N7sNgarb9++2LdvX4Pt77//vvXfO3bsiNOnT9vrI8mBPNxdwf54gQmfmzF0sxkbHlUx6j6uxERERO6JfwOS3fQyCeydbEDEPQKPfqbh9WMapOTVLCIicj8MWGRX3u0Ftj2m4qVwBb87qGPKDg2XahmyiIjIvTBgkd2pisDSKBVbxqrIOisx6FMzjv3EkEVERO6DAYtazMSeCg5PMcDkATy82Yy/clJSIiJyEwxY1KKCTQJfTzbgX3oL/DJLw4LdZtSYeTWLiIhcGwMWtbgOBoH3RxqweqSKD/MkoreYUVjBkEVERK6LAYtazdy+CvZONuBSLTDwUzM+K+YtQyIick0MWNSqBvgIHJpiwPDuAhO2afj3Qxo0nVeziIjItTBgUavrahT4dKyKZVEKXjum47HPNVzgEjtERORCGLCoTShC4KUIFdsfU3HsJ4kBG83Y9SNvGRIRkWtgwKI29UgPBUemGBDcWWBUhoZf7dVw5RqvZhERkXNjwKI25+8pkDlRxZsPK3g/R0fERjO+LuXVLCIicl4MWOQQFCHw634qjj1hQLf2AiO2aHhxv4ZqzpnlFqSU+MtJDX/+vl1bN4WIyC4YsMihhHQR2D1JxYohCt4+qSNykxmHyvmfqSu7apZ4epeG57/W4cFTTUQugt0ZORxVEXixv4ojUwzo7CHwWGZHvHRAQ63Gq1muprBCIjrdjE2nJD4ereI3fevauklERHbBgEUOK7SrwN7JKl4Kq8MbJ3QM3GTG4fMMWa5i+2kdAz81o6IO2DfZgJm92B0Rketgj0YOzaAILHqgDoenGNBOAYZsNmPJYQ11vJrltKSU+I9jGuL+oWGor2Xi2f73iLZuFhGRXTFgkVN4yFvgm8cNWDxAwbKjOoZsNuP4TwxZzqayTmLaFxp+f1DH7yMUbBmroquR4YqIXA8DFjkND0VgyUAV3zxugCaBqE/N+Lf9Gi7WMmg5g39ekhiy2YwdZyQ2jVGxNEqFqjBcEZFrYsAipxPpI3DwcQOWRCpI+V5H7/VmrPxOwzWuaeiw0ot0DP7UDAngwOMGPB7EroeIXBt7OXJKRlXgDwNU5M0w4IkggV/v0/FQqhlbi3RIyaDlKHQpseSwhvjtGkbfJ/BNvAEPdOFVKyJyfQxYuCnzKgAACytJREFU5NT8Ogq8N9KAo08Y0KOTwKTtGsZ8pnF8lgO4VCsxaZuGPx3RsXSQgrQxKkztGK6IyD0wYJFLCL9H4IvxKraMVXH6imXx6HlfmXH2KoNWWzh4XkfUp2bsPSeREafiDwNUKILhiojcBwMWuQwhBCb2VHBimgFvPaxg0ymJPuvNWHqES+60lh+vSCTsMmPwpxo8PYBDUwx4LIDdDBG5H/Z85HI8FIF/7acif4YBiQ8qePWojr6fmPFRvg6d47NaRLVZ4rWjGkI+MSOjROIv0QoOPm5ALxOvWhGRe2LAIpfV1SjwxlAVJ6cZENVNYHamhqGbNWwq1KHxF4d2IaXEJwU6Htxgxh8P60h8UEHedAOeDVVh4BQMROTGGLDI5fX2EkgbY0DWRBXtFOCJLyxXWt76TkNlHYNWUx25IBGzVcOMnRr6ewtkTzPgjaEqunDiUCIiBixyHyP9FOyZbMA38SqG+Aq8uF+H/98tk5WeqmTQslXpVYm5WWYM2mRGea3EtsdUpI8zIITTLxARWTFgkdsZ7Kvg76MNKJxpwHOhCtbk6ui13oxpX5jxdSnn0bqdGrNlDcE+n5iRXiTxdrSCY08YMNaf3QgR0c3YM5Lb8vcUWD5YRcksA94epuBEucTwLRqGbNbwcb7OmeF/JqVEWqGO0FQz/v2Qjmf6KsibYcDzHGdFRHRbDFjk9jp5CDwXquL7Jw3YOk6FyQN4KlPD/evMWH5MQ3mNewat4iqJZUc1PLjBjGlfaHiwi8CJaQa8+TAXaCYiuhNDWzeAyFEoQmBCoMCEQMvVrDdPaPjjER2vHNExpodAfE8FEwMFund03XBRWSeRdkrig1wdmWclOhqAqUEC7w5XEHsf/3+MiMhWDFhEt/CQt8DqGAP+Y7DEB3k60oskFuzRICUwxFcgvqfA5J4KHuximeDUmWm6xM4fLXVuPCVRbQZG3Sfw1xgVU+8X8PRw7vqIiNoCAxZRI3w7CLzYX8WL/YHz1RKflUikF+lYelTHSwd19DIB8T0VTO4pEN1dONWYpOxyS6ham6/jx6tAXy/g5QEK/qW3gkBP56mDiMgR2S1g5eXlYc6cObhw4QK8vLzw17/+FWFhYQ32W716NZYvXw5d1zF69Gi888478PDwsFcziFpMtw4Cc0IE5oQoqDFbrvpsLpL4uEDHf54AvI3AhADLla1x/gKdHWxh4zpNIucSsOusjg/yJA5fkPA2ArN6KfhFH4GobsLpr8YRETkKuwWsxMRELFiwAAkJCUhNTUVCQgIOHjxYb5/CwkIsXrwYR44cQffu3REfH49Vq1bhhRdesFcziFpFe4PA+ECB8YHAX6SCwxck0oskNhfp+DBfgyqA4M7AA10EHuwifv6n5XlLT8QppcSZK8C35RLflkuc+PmfOZcAswQ8FEsQfHmAivEBAu1UhioiInuzS8AqKyvDoUOHsH37dgDA1KlTsXDhQuTn56N3797W/VJTUzF58mTce++9AIBnn30Wy5YtY8Aip6YIy9WfqG7AnwapKKyQ+OJHie8vSnx/SeKTH3Scqvrf/e/tgAah68EuAj063f14rso6ie8uSpwoR70wdanO8npnD6C/t8DwexU8H2oZW9bfW8DkYFfXiIhcjV0CVklJCfz8/GAwWA4nhEBgYCCKi4vrBazi4mL07NnT+jwoKAjFxcWNHruqqgoVFRXW50ajEQDw+uuv49VXX7U+d3W1tbVuVzPgnHXfbxKYf9Mix1fNErmXge8vSuRctoSvPed0/E8uUKtZ9vFQAEUAlnlOJcxaByhKHSQsx5KwvHarSSNUAYR4WcLUWH8F/X8OUoGezjMI3xnPNRHR7Tj8IPeYmJh6z5OSkvDcc89hxYoVeO6552AymdqoZa2roqLC7WoGXKvuQAEEegPjvAHcb9mmSaD4ikBupYKSqwqkBIQAampqsGTxYrzypz+hffv2EAK4HpMEYH1uVCUeMOkI6ayjvXrTB14DLl5stfKazZZz7e3t3cqtIiJqGrsErICAAJw9exZmsxkGgwFSShQXFyMwMLDefoGBgSgoKLA+P3XqVIN9bpaVlYWIiAjrc6PRiNraWgCWztbZ/9K11fWrg+5UM+AedXe7Bxh407aKCoklWSn4VfgKmEyd26Rdrc0dzjURuQ+7zBzo6+uLyMhIrF27FgCQlpYGf3//ercHAcvYrPT0dJSWlkJKiXfffRczZ85s9Nienp4wmUzWB28dEBERkaOz2y3ClJQUJCQkYNmyZTCZTFizZg0AYN68eZg8eTImT56M4OBgvPLKK4iOjgYAxMbGIjEx8ZbHq66uBgB8//33DV6rqrKMGD527Bg8PT3tVYJDc8eaAdbtTnXbWvMDDzyAjh07tlaziIiaREgpHXKhtY8++gizZ89u62YQkYM5fPgwIiMj27oZRESNctiAdeHCBWzbtg1BQUHo0KFDWzeHiBwEr2ARkTNw2IBFRERE5KzsMsidiIiIiP4XAxYRERGRnTlUwMrLy8OwYcMQEhKCqKgoZGdn33K/1atXo0+fPujVqxfmz5+Pa9eu2fSao7Kl7p07d2Lw4MEIDQ1FWFgYkpKSoOs6AMt8YqqqIiIiwvq4cb4xR2RLzbt27UKHDh3q1XX916WA657rNWvW1KvZx8cHTzzxBADnPNe/+tWvEBQUBCEEjh07dtv9XO17TURuTjqQUaNGyTVr1kgppdywYYMcNGhQg31++OEH6efnJ8+ePSt1XZeTJk2Sb7/99h1fc2S21H3kyBFZUFAgpZSyurpaRkdHW99TWFgovby8Wqu5dmFLzZmZmTI8PPyW73flc32zsLAwmZqaKqV0znOdlZUlS0pKZM+ePeXRo0dvuY8rfq+JyL05TMA6d+6c7Ny5s7x27ZqUUkpd12X37t1lXl5evf1WrFghExMTrc8zMjJkdHT0HV9zVLbWfbMXXnhBLlmyRErpfH/p2lpzYwHLXc71/v37Zbdu3WRdXZ2U0vnO9Y0aC1iu9r0mInKYW4SNLRh9o8YWjG7KYtJtzda6b1RaWorU1FRMnDjRuu3KlSuIiopCZGQkXn31VWia1uJtb6q7qbmgoACRkZGIiorCO++8Y93uLud69erVePrpp+Hh4WHd5kzn2lau9r0mInL4xZ6pvoqKCkyaNAlJSUkYNGgQAMDPzw9nzpyBr68vysvLMWPGDLzxxhtISkpq49Y2T2RkJE6fPg0vLy+cPn0a48ePh4+PD6ZPn97WTWsVV65cwbp167B//37rNlc910RErsZhrmDduGA0gEYXjC4qKrI+v3HB6MZec1S21g0AlZWViIuLQ3x8PBYtWmTdbjQa4evrC8CyUO7cuXOxe/fu1imgCWyt2WQywcvLCwDg7++PWbNmWety9XMNABs2bEBYWBhCQ0Ot25ztXNvK1b7XREQOE7DssWB0UxaTbmu21l1VVYW4uDjExcXh5ZdfrvdaWVmZ9VdVtbW12LhxIwYMGNA6BTSBrTWfPXvW+kvJyspKbN261VqXK5/r61avXo1nnnmm3jZnO9e2crXvNRGRwwxyl1LKnJwcOXToUNmnTx85cOBA+e2330oppXzmmWfk5s2brfutWrVKBgcHy+DgYDl37lzrAOA7veaobKl76dKl0mAwyPDwcOtj6dKlUkop09LSZFhYmOzfv78MDQ2VCxculDU1NW1Wjy1sqXnlypUyNDTUWteSJUukruvWY7jqub6+n6enp6yoqKj3fmc81wsWLJA9evSQqqpKX19f2atXLyml63+vici9cakcIiIiIjtzmFuERERERK6CAYuIiIjIzhiwiIiIiOzs/wN+oHRKUAqYRwAAAABJRU5ErkJggg==" }, - "execution_count": 26, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } diff --git a/examples/exampleLinearModel.ipynb b/examples/exampleLinearModel.ipynb index 7dc15f2..9311033 100644 --- a/examples/exampleLinearModel.ipynb +++ b/examples/exampleLinearModel.ipynb @@ -316,9 +316,7 @@ { "cell_type": "code", "execution_count": 9, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# x[1] corresponds to the intercept\n", @@ -372,7 +370,7 @@ " output[\"mean_x1y\"] = mean(simX[startT:nbDraws,1] .* y[startT:nbDraws])\n", " output[\"mean_x2y\"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws])\n", " output[\"mean_x1y^2\"] = mean((simX[startT:nbDraws,1] .* y[startT:nbDraws]).^2)\n", - " output[\"mean_x2y^2\"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)\n", + " output[\"mean_x2y^2\"] = mean((simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)\n", "\n", " return output\n", "end" @@ -463,7 +461,7 @@ "name": "stdout", "output_type": "stream", "text": [ - " 67.708742 seconds (236.31 M allocations: 50.348 GiB, 52.94% gc time)\n" + " 87.049598 seconds (236.31 M allocations: 50.348 GiB, 54.97% gc time)\n" ] }, { @@ -650,9 +648,9 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 16, "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [ { @@ -664,7 +662,7 @@ " -4.75213 -1.2148 3.25637" ] }, - "execution_count": 27, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -679,32 +677,41 @@ "metadata": {}, "source": [ "Once the asymptotic variance has been calculated, a summary table can be obtained using the\n", - "function `summary_table`" + "function `summary_table`. This function has four inputs:\n", + "* a SMMProblem\n", + "* the minimizer of the objective function\n", + "* the length of the empirical sample\n", + "* the confidence level associated to the test H0: $\\theta_i = 0$, H1: $\\theta_i != 0$" ] }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 47, "metadata": { - "scrolled": true + "scrolled": false }, "outputs": [ { - "ename": "LoadError", - "evalue": "\u001b[91mArgumentError: Length of iterable does not match DataFrame column count.\u001b[39m", - "output_type": "error", - "traceback": [ - "\u001b[91mArgumentError: Length of iterable does not match DataFrame column count.\u001b[39m", - "", - "Stacktrace:", - " [1] \u001b[1mpush!\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::DataFrames.DataFrame, ::Array{Float64,1}\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/julien/.julia/v0.6/DataFrames/src/dataframe/dataframe.jl:1028\u001b[22m\u001b[22m", - " [2] \u001b[1msummary_table\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::SMM.SMMProblem, ::Array{Float64,1}, ::Int64, ::Float64\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m/home/julien/.julia/v0.6/SMM/src/econometrics.jl:154\u001b[22m\u001b[22m", - " [3] \u001b[1minclude_string\u001b[22m\u001b[22m\u001b[1m(\u001b[22m\u001b[22m::String, ::String\u001b[1m)\u001b[22m\u001b[22m at \u001b[1m./loading.jl:522\u001b[22m\u001b[22m" - ] + "data": { + "text/html": [ + "
EstimateStdErrortValuepValueConfIntervalLowerConfIntervalUpper
10.5876030.018627431.54510.00.5569640.618243
20.581410.0061311294.82930.00.5713250.591494
30.7626680.00570646133.650.00.7532810.772054
" + ], + "text/plain": [ + "3×6 DataFrames.DataFrame. Omitted printing of 1 columns\n", + "│ Row │ Estimate │ StdError │ tValue │ pValue │ ConfIntervalLower │\n", + "├─────┼──────────┼────────────┼─────────┼────────┼───────────────────┤\n", + "│ 1 │ 0.587603 │ 0.0186274 │ 31.5451 │ 0.0 │ 0.556964 │\n", + "│ 2 │ 0.58141 │ 0.00613112 │ 94.8293 │ 0.0 │ 0.571325 │\n", + "│ 3 │ 0.762668 │ 0.00570646 │ 133.65 │ 0.0 │ 0.753281 │" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ - "summary_table(myProblem, minimizer, T, alpha)" + "df = summary_table(myProblem, minimizer, T, 0.05)" ] }, { @@ -723,7 +730,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 18, "metadata": { "collapsed": true }, @@ -735,11 +742,28 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "StatsModels.DataFrameRegressionModel{GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,Base.LinAlg.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}\n", + "\n", + "Formula: y ~ 1 + x1 + x2\n", + "\n", + "Coefficients:\n", + " Estimate Std.Error t value Pr(>|t|)\n", + "(Intercept) 0.571592 0.00836568 68.3259 <1e-99\n", + "x1 0.588239 0.00218787 268.863 <1e-99\n", + "x2 0.767633 0.00219593 349.572 <1e-99\n" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ols = lm(@formula(y ~ x1 + x2), data)" ] @@ -762,7 +786,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 20, "metadata": { "scrolled": true }, @@ -816,7 +840,7 @@ "name": "stdout", "output_type": "stream", "text": [ - " 41.372656 seconds (72.39 M allocations: 13.920 GiB, 30.75% gc time)\n" + " 55.815914 seconds (72.39 M allocations: 13.916 GiB, 29.16% gc time)\n" ] }, { @@ -828,7 +852,7 @@ " Plot{Plots.PyPlotBackend() n=1}" ] }, - "execution_count": 25, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -847,7 +871,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 21, "metadata": { "scrolled": false }, @@ -856,7 +880,7 @@ "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAYAAAByNR6YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVWX+wPHPc85VBBGU3BgRTBZNMwFFy13LMsslGcvKJq1QK9sbbXPStnGammmb0tKxxZ+2qJVZlpZLOGopaplmLEqA4RYqIgrec57fH4yMCCI7d/m+X6/z0nvvOfc+Xy98/Z7nPOd5lNZaI4QQQgghaoxR3w0QQgghhPA0UmCJOrFmzRqUUkyfPr1a7zN9+nSUUqxZs6ZG2iWE8Ew1lXOEqCopsIQQQohKateuHe3ataux93v11VcZP348l1xyCQ6HQ04kPYCjvhsghBBCeLt7770XgODgYFq0aMG+ffvquUWiuqQHSwghhKhny5YtIzs7m99++40RI0bUd3NEDZACS1RZYWEhr776KldddRVt27bFx8eHli1bMmrUKLZu3Vqh9zjdzX7kyBEmTpxI69atadSoETExMSxcuLDcYxcsWEB0dDS+vr4EBwdz3333ceLEiRpvoxDCva1bt44BAwbQpEkTmjZtSnx8PKmpqaX2O3DgAA888AARERH4+PjQvHlz4uPj+emnn4r3SU9PRynFr7/+yq+//opSqng7Pd6rKnnnmmuuoXXr1rUSv6gfcolQVFlOTg73338/ffv2ZejQoTRr1ozdu3ezdOlSli9fzrfffktcXNx536ewsJArrriCvLw8brnlFo4fP86HH37ITTfdxKFDh7jnnntKHfPaa6/x5ZdfMmLECAYNGsSXX37JK6+8wqFDh/i///u/Gm+jEMI9bdy4kb/+9a8MGTKEe+65hx07dvDxxx+TmJjIxo0bad++PQBpaWkMGDCArKwsrrzySkaOHMmBAwdYvHgxX331Fd988w09e/akadOmPPnkk7z00ksA3H///cWfNWDAAEDyjvgvLUQVnTx5UmdlZZV6/qefftL+/v76iiuuKH5u9erVGtBPPvlkiX3DwsI0oPv166cLCgqKn8/MzNTNmzfXPj4+JT7jySef1IAODAzUu3btKn4+Pz9fR0VFacMw9N69e6vURiGE5zidcwA9a9asEq/NmjVLA/raa68tfq5Xr17aNE395Zdfltj3l19+0U2aNNFdunQp8XxYWJgOCwsr87Orm3cmTpyoAb169epy9xOuTS4Riirz8fGhTZs2pZ7v3LkzAwcO5Ntvv+XUqVMVeq/nnnuOhg0bFj8OCQnhvvvuo6CggPfff7/U/vfddx8dOnQofuzr68uNN96IbdskJSXVShuFEO4nKiqKhISEEs8lJCQQGRnJ559/zsGDB9m6dSvr16/n1ltv5aqrrirz+O3bt5e4VFgeyTsC5BKhqKZt27bx/PPPs27dOvbt21cqaRw6dIjg4OBy38PhcHDZZZeVer5v374AZY5Z6NatW6nnQkJCADhy5EiNt1EI4Z569+6NYZTsSzAMg969e5OSksIPP/xASkoKAPv37y9z3qxdu3YV/3nxxRdX6HMl7wgpsESVrV+/nkGDBgFw5ZVXEhkZib+/P0opPvnkE3744QcKCgrO+z7NmzcvlQABWrVqBcDRo0dLvRYQEFDqOYej6MfZsqwab6MQwj2dziPnev7o0aPk5OQA8Pnnn/P555+f872OHz9eoc+UvCNACixRDc8++ywFBQUkJibSp0+fEq9t3LiRH374oULvc+jQIWzbLlVk7d+/H4DAwMB6b6MQwj2dziPnej4wMLD4hO3VV19l8uTJ1f5MyTsCZJoGUQ1paWkEBQWVSiD5+fls2bKlwu/jdDrZsGFDqecTExMBiImJqfc2CiHc03/+8x9s2y7xnG3brF+/HqUUXbt2pWfPngBl5qFzMU2zRG/5mSTvCJACS1RDWFgYhw8fZseOHcXPWZbFww8/zMGDByv1Xo899hiFhYXFj7Oysnj55Zfx8fFhzJgxLtFGIYT7SU5O5q233irx3FtvvUVycjLXXHMNLVq0oEePHvTs2ZOFCxfywQcflHoP27ZZu3ZtieeCgoI4dOgQJ0+eLLW/5B0BcolQVMM999zDihUr6NOnD9dffz2NGjVizZo17N27lwEDBlR4Ha3g4GCOHz/OJZdcwrBhw4rnwfr999955ZVXyrwbp67bKIRwT1dddRX33nsvX3zxBZ07d2bHjh189tlnNG/enJdffrl4v4ULFzJw4EDGjBnDSy+9RGxsLL6+vmRkZLBhwwYOHjxYopgaNGgQmzdv5uqrr6Zv3740bNiQfv360a9fvyrlnZkzZxYPpj/dkzZz5kzefvttAEaOHMnIkSNr7x9K1Lz6nidCuLdFixbp2NhY7efnp5s3b66vv/56nZaWpm+99VYN6D179mity58HKywsTOfk5OgJEyboVq1aaR8fH921a1e9YMGCUp93eh6ssuaHmTdvngb0vHnzqtRGIYTnODPnJCYm6v79++vGjRvrgIAAfd111+mUlJRSx+Tk5OgnnnhCX3zxxdrX11f7+/vryMhIfdNNN+klS5aU2PfYsWM6ISFBBwcHa9M0S+W3yuad/v37F8/bVdZ2du4Urk9prXU91HVCABSvRp+enl6v7RBCCCFqkozBEkIIIYSoYVJgCSGEEELUMCmwhBBCCCFqmIzBEkIIIYSoYdKDJYQQQghRw6TAEkIIIYSoYS5bYJ1eUiA/P7++myKEECVIfhJCnI/LFli7du2iW7duxTPbnu3o0aN13KL6540xg8TtTdwlZslPZfPGuL0xZpC4K8JlC6zzOdcim57MG2MGidubeErMnhJHZXlj3N4YM0jcFeG2BZYQQgghhKuSAksIIYQQooa5ZYH1yPcWf93RsL6bIYQQpWTmaQZ87cfuXJliUAhv5pYFVotG8NIvDUk9KglMCOFaLmgEvxcqHtzonWNUhBBF3LLAuquTQXMfzZNJksCEEK7Fz6F4+pICPv1VszzTru/mCCHqiVsWWL4OxZ8vKmRhmubH36UXSwjhWka0cTIwWHHfeosCS3KUEN7ILQssgJvbnSI8AB7fLL1YQgjXohS82stk9zF4abv0Ygnhjdy2wGpgwFPdTJZlaNbvlwQmhHAtnYMU93Q2eHqrzd7j0oslhLdx2wIL4IZwxSVB8NgmG60lgQkhXMv0bgaNG8Cfv5OediG8jVsXWIZSPNvdZG22ZuVeKbCEEK4lsKHibz1MFqZp1mZLT7sQ3sStCyyAa0IVvVop6cUSQrikP0UqLm2puGe9hdOWHCWEt3D7AkspxXNxBkmHNIv3SPISwlulpKTQq1cvoqKiiIuLY8eOHaX2WbVqFT169KBTp0507tyZKVOmYNv/61latmwZHTt2JDIyklGjRpGbm1vtdhlK8Vovk59y4I2d0oslhLdw+wILoH+wwVUhimlJcoYohLeaOHEiEyZMIDk5malTpzJu3LhS+zRr1oz333+fnTt3kpSUxPr163n33XcByMvL4/bbb+eTTz4hJSWFP/zhDzz99NM10rZuLRQJHQ2mJdkcOCE5Sghv4BEFFsCz3U12HYH3UiR5CeFtDhw4wObNmxk7diwA8fHxZGZmkpqaWmK/mJgY2rdvD0CjRo2Ijo4mPT0dgOXLlxMTE0PHjh0BuOuuu1i4cGGNtfHZOANDwWObZMC7EN7AYwqsbi0Uf7xQMX2LTOwnhLfJzMwkODgYh8MBFA0dCA0NJSMj45zH7Nu3j0WLFnHttdcCkJGRQVhYWPHr7dq1Izs7G6fTec73yMvLIzc3t3grKCg4577NGyme6W4w9xfN9wfkUqEQns5R3w2oSU93N+m8yMnsn23uvdis7+YIIVxUbm4uw4YNY8qUKXTv3r3K79O/f/8Sj6dMmcLUqVM5fPhwmfuPbgWzAv2Y9K3FioH5GKrKH+2SzhW3J/PGmEHiLktQUFCJxx5VYHVsqhgXqXhmq81tHQz8G3hY9hJClKlt27bFvU0OhwOtNRkZGYSGhpba99ixYwwZMoQRI0bw4IMPFj8fGhrKypUrix+np6eX6BUry9q1a4mOji5+7OPjg4+PD1A62Z72Rj+bPp9ZfHqwKbd39JiLCMXOFbcn88aYwbvitrXm6S02o1opwoOaVeiYUr/dFbkTB2Du3LlERkYSHh5OQkICp06dAsq/Syc9PR3TNImOji7e0tLSqhpvmZ7sZnK0UJanEMKbtGzZktjYWObPnw/A4sWLCQkJISIiosR+eXl5DBkyhCFDhvDEE0+UeG3IkCFs2bKFXbt2AfD6668zZsyYcj/X39+fgICA4u10cVWe3q0NxkYoHtlkcbhAhjMI4Q4S92mmb7H57UTFO25KFVgVuRNnz549TJs2jcTERFJTU9m/fz9vvvkmUP5dOgBNmjRh27ZtxVt4eHgVQj23UH/FnZ0M/v6jTc5JSV5CeIvZs2cze/ZsoqKimDlzJvPmzQPgjjvuYOnSpQC8/PLLfP/99yxZsqT4JO/ZZ58FinLTnDlzGDlyJBEREWRlZTFt2rRaaevzPU1OWvBkkpwICuEO3k62CQ+A7kEV/50t0fd9+k6cFStWAEV34kyePJnU1NQSZ4KLFi1i+PDhtG7dGoBJkybx3HPPcffddxMTE1O839l36dSVx6IN5uyy+dsPNn/rKWOxhPAGHTp0YMOGDaWenzNnTvHfH3/8cR5//PFzvsfw4cMZPnx4rbTvTMF+iidjDaZ+b3NHB4NLLpDhDEK4qrxTmo92ax7paqAq8ataogeronfilHW3TVl365x9lw7A8ePHiYuLIzY2lqeeegrLKv+W5crcpXNaS1/FA10MXtlh85sssiqEcEH3djaICoR71luyCoUQLmzRHk2+E/4UVbkxk7U2yL2su3SCg4PZu3cvLVu2JCcnhxtuuIEXX3yRKVOmnPN9KnuXzmm3t4V/7fDniY35vBBz/qLMHchdG97FG+M+X8yeNKi2oal45TKTK5dbvJ+muTFCerGEcEVvJ9tc3kYR6q/Iyan4cSUKrIreiRMaGlpicHp6enqJfc51l46Pjw8tW7YEihLlbbfdxoIFC8otsKpylw5AEPBotMVjmxSPx/kRHuAZycuT/oOpDInbe3hTzINDDEa1s3n4O4thYUrufBbCxezO1azN1swfWPnhRiX6uyp6J058fDxLly5l3759aK2ZNWtW8d025d2lc+DAgeK7DQsKCliyZEmJMVtlqcpdOqfd3dmghS88mSQzJwshXNM/LjXJKYBntsqAdyFczTspNgEN4Lp2lT/5KXVBsSJ34rRv354ZM2bQu3dvIiIiaNGiBRMnTgTKv0tn3bp1xMTE0LVrV2JjY2ndunW5A06ry8+h+EuMwYJUzfYcGeMghHA9YU0Uj0Yb/GO7zc+HJU8J4SpsrXkn2eaGcIWfo/IFltIuOrpyy5YtdOvWjaSkJGJjY0u9npOTU6FLCYWW5qKPnFwcpPj0SveeV7WiMXsaidt7uEvMNZWfTjvp1HRZ7KRNY8Xqa0xUZW5VciHu8v3VJG+MGbwj7tW/2Qz63OI/w016tSrqj6pM3J43jfBZGpqKp7qbLP1Vs2G/dMELIVxPI4fijT4ma7M178iC9UK4hHm/2EQGwGUtq3bC4/EFFsCN4YouQfDn72y5HVoI4ZKuaGNwc4Ti4Y0Wh2SSZCHqVW6hZtEezbgoo8o9yl5RYBlK8Y9LTf6zX/PhbklcQgjX9GJPE0vDn7+TG3OEqE+L9mhOWvCnyKqXSV5RYEHR2eHwMMWU7y1OOKXIEkK4nlZ+iud7mrydrFmbLUMahKgv836xGdxGEeJf9fGQXlNgAbzQ0yQ7H174URKXEMI13d5B0buVYtI6iwJLTgaFqGupRzXr9hddHqwOryqwIgMV911sMPMHm72yhI4QwgUZSjGrj0nqUXj+BzkZFKKuvZNiE9gQRlZh7qszeVWBBfBEjEFjBzzyvYxxEEK4pouDFA9dYvDsNpuUo3IyKERdseyiua/GtDfwrcLcV2fyugIrsKHi2TiT+amajTJtgxDCRf0l1iDYF+76jywGLURdWZ2tyTwO46KqPxed1xVYALdFKboGwf0bbWxJXEIIF+TnULzex+TrvZoFaZKnhKgLbyfbdAiEnlWc++pMXllgmYbi5V4m3x3QLEiVxCWEcE1XtzUYfaHiwY0WhwskVwlRm44Wahbv0YyvxtxXZ/LKAgugf7BB/IWKRzZZHD8liUsI4ZpeuszkpBMe+V6GNAhRmz7crSm0YWw15r46k9cWWAB/72Fy6CT8Te7UEUK4qD80VjwXZ/DmLpv1Mm5UiFrzdrLNlW0UbRrXzFqgXl1gXRigeLCLwd9/tPn1mPRiCSFc06SLDOJaKCYmWpyyJVcJUdOSj2jW18DcV2fy6gIL4NGuBk0bwlSZtkEI4aJMQzG7j8nOI/APmShZiBr3dopN04YwIqxmeq9ACiyaNFT8tYfJB7s16/ZJ4hJCuKaY5or7OhvM2GKzJ1d6sYSoKZateTfF5sZwg0bVnPvqTF5fYAH8KVLRvbni/g0ybYMQwnU91d2geSO4e73MjSVETfnmN83eGpr76kxSYFG0NMVLlxkkHdK8kyxJSwjhmvwbKF7tZbI8U7Noj+QqIWrCvGSbTk0hrkUtF1gpKSn06tWLqKgo4uLi2LFjR5kHzp07l8jISMLDw0lISODUqVMArFq1ih49etCpUyc6d+7MlClTsO3/XXpbtmwZHTt2JDIyklGjRpGbm1ujAVVV79YGN4YrHttkcaxQEpcQwjWNaGcwIkxx3waLo5KrhKiWIwWaj9OLBrfXxNxXZypVYE2cOJEJEyaQnJzM1KlTGTduXKmD9uzZw7Rp00hMTCQ1NZX9+/fz5ptvAtCsWTPef/99du7cSVJSEuvXr+fdd98FIC8vj9tvv51PPvmElJQU/vCHP/D000/XaEDV8bceJkcL4bltMhZLCOG6Xu1lklsIj2+SXCVEdXyw2+ZUDc59daYS73jgwAE2b97M2LFjAYiPjyczM5PU1NQSBy1atIjhw4fTunVrlFJMmjSJhQsXAhATE0P79u0BaNSoEdHR0aSnpwOwfPlyYmJi6NixIwB33XVX8XGuoK2/YkpXg39st9ktg0iFcCsV6X1PT09nwIABBAYGEh0dXeK1NWvW4OvrS3R0dPF24sSJump+pbT1Vzzd3eD1nTI3lhDVMS9ZMyREEexXs71XcFaBlZmZSXBwMA6HAwClFKGhoWRkZJQ4KCMjg7CwsOLH7dq1K7UPwL59+1i0aBHXXnvtOY/Lzs7G6XSes4F5eXnk5uYWbwUFBVUIs+KmdDVo6Qt//k6mbRDCnVSk9z0gIIBnnnmGBQsWlPkeHTp0YNu2bcWbr69vLbe66u7tbNCjpeK2tRYnnHJCKERl/XxY892BoqVxaoOjVt4VyM3NZdiwYUyZMoXu3btX+X369+9f4vGUKVOYOnUqhw8frm4Tz2laJwcTN/my9Jej9GnhOoVWbcbsyiRu73G+mIOCgsp8/nTv+4oVK4Ci3vfJkyeTmppKREREieP79OnDmjVraqzN9cU0FP/uZxKzxMn0JJu/9TTru0lCuJV3Umya+cCwGpz76kwlCqy2bdsW9yg5HA601mRkZBAaGlrioNDQUNLS0oofp6enl9jn2LFjDBkyhBEjRvDggw+WOG7lypUljjuzx6wsa9euLdGV7+Pjg4+PD3DuZFtdCc00b/9q8Zef/Ei6zoFp1M4/flXUVsyuTuL2HlWJubze9zMLrPNJS0sjNjYW0zQZP348d911V7n7n+5hP+3M/FQXOjVTTO9m8MRmm1EXKnq2lBvDhagI53/nvrop3MDHrIMCq2XLlsTGxjJ//nzGjRvH4sWLCQkJKZWg4uPj6dOnD9OnT6dVq1bMmjWLMWPGAEUJZ8iQIQwZMoQnnniixHFDhgzh7rvvZteuXXTs2JHXX3+9+Lhz8ff3JyAgoCZirTClFC/3MujxicVbu2wmdZIzQyE8XWxsLFlZWQQGBpKVlcXQoUNp3rw5119//TmPqY8e9rPdHgIfpPhx6yqL1Zfn41OP6Up6Xb2Hu8f99T6T7Hw/rmt9jJycio9jLC/us08OS3UdzZ49m3HjxvHcc88REBDAvHnzALjjjjsYPnw4w4cPp3379syYMYPevXsDMGDAACZOnAjAyy+/zPfff8/x48dZsmQJAKNHj+bxxx+nSZMmzJkzh5EjR+J0Orn44ot55513KhxYXYprYXBblM0jm2xGtDNqZQCcEKJmVLT3vTxnnsiFhIRw4403kpiYWG6BVR897GWZf7km9mMnr6UH8mxc/Z4QSq+r93DnuBdvdXJxM82g9oGVnp6honGXKrA6dOjAhg0bSu04Z86cEo8TEhJISEgotd/jjz/O448/fs4PPF2kuYPne5p8luHk/g0WH1xea8PVhBDVVNHe9/JkZ2fTqlUrDMPg2LFjLFu2jNtvv73cY+qjh70sFwcppsUULaMzqp1BtxqeMFEIT7I/X/NJuua5uJqf++pMcsG+HBc0Urx8mcmHuzXLfpVboYVwZbNnz2b27NlERUUxc+bMEr3vS5cuBSA/P5+QkBBGjx7Nzp07CQkJ4dFHHwVg8eLFdOnSha5du3LppZcyePBgxo8fX2/xVNYj0QZdgmD8t04KLbmrUIhzeW2nTQMDbutQuyWQdMucx5hwxbspirv+Y9E/WNGkoZwZCuGKKtL77ufnR1ZWVpnHT548mcmTJ9da+2pbA0Pxdn8H3T928uw2mxndZOyoEGfLd2pe32lzR0eDZj61+/+59GCdh1KK13ub/F4A05KkF0sI4bq6XqB4LMbgua02236XXiwhzvZ2ss2RQrivc+2XP1JgVcCFAYqnuhm88pPN9wekyBJCuK7How0uagbj1zo5ZUuRJcRplq3553abP16ouDCg9q9GSYFVQfddbBDTHCYkWpK0hBAuq6FZdKlwew7MlHVVhSi29FdNai481KVuSh8psCrIYSje6utg+2H453ZJWkII1xXbXDG1q8HTW22258gJoRAAL2y36dta0aOOJuSVAqsSYpsrHrjY4MkkmzRZDFoI4cL+EmsQGQDj11o4pdddeLkN+23W79c8fEndlT1SYFXSjG4GrXxh0joLrSVpCSFck4+pmNffZOvvmr//KL3uwru9uN0mKhCuDa27mQCkwKqkxg0Us/qYfL1XMz9VCiwhhOvq0dLg4S4G05Nsdh6WfCW8U1qu5uN0zYNdDIxanFj0bFJgVcGQtgY3hSse2GBx6KQkLSGE65rRzeDCJnKpUHivl7bbBPnAnyLrtuSRAquK/nmZiQ08uMGq76YIIcQ5NXIUXSrcdFDLDTrC6+Sc1Pw72ebuTga+jrqdKFwKrCpq6at4safJe6malVmStIQQruuyVgYPdDGYlmTzyxHpxRLeY9bPNraGuzrVfbkjBVY1jItSDAxWTFpnke+UpCWEcF1Pdzdo2xhu+9bCkkuFwgsUWJpXd9j8KdKgpW/dL3MnBVY1KKWY3ddkbz48tUV6sYQQrsvPoXi7v8nGA5q/ygSkwgssSNXsOwEP1tHEomeTAquaIgMVf4kxeOFHWftLiOpyyv/7tap3a4PHog2mb7H5Tpb9Eh5Ma82L2y2GhSo6NK373iuQAqtGPHyJwUVNIUG63oWosoc2Wty43re+m+Hx/hJr0L2F4qZVFscKJV8Jz/RVlmbHYep0YtGzSYFVAxqairf6miQd0ry2U84KhaisDftt/rnd5orWzvpuisdrYCj+b6DJgZNwz3q5C1p4phd+tIlroejbun56r6CMAislJYVevXoRFRVFXFwcO3bsKPPAuXPnEhkZSXh4OAkJCZw6dQqA9PR0BgwYQGBgINHR0SWOWbNmDb6+vkRHRxdvJ06cqIWw6t6lrQzu7mTw+CabjDw5KxSiogotzYREi+4tFHeEn6rv5niF8ADFa71M3knRfJAmJ4XCs2z7XfPNb5qHuhioOpxY9GylCqyJEycyYcIEkpOTmTp1KuPGjSt10J49e5g2bRqJiYmkpqayf/9+3nzzTQACAgJ45plnWLBgQZkf2KFDB7Zt21a8+fp6ziWBZ+MMmvnAbWstbFlGR4gKeeFHm5+PwFt9Tcz6y4Ve50+RihvaKyaus+SkUHiUF3+0CPOH+AvrN6GUKLAOHDjA5s2bGTt2LADx8fFkZmaSmppa4qBFixYxfPhwWrdujVKKSZMmsXDhQgCCgoLo06cPjRs3rqMQXEdAw6K7dL75TSb0E6IiUo5qntpq81AXg64XSHVVl5RSvNHHJKABjF0t40eFZ8jK07yfprn/YgOH4UIFVmZmJsHBwTgcDqDoFzA0NJSMjIwSB2VkZBAWFlb8uF27dqX2OZe0tDRiY2OJi4vj9ddfP+/+eXl55ObmFm8FBQUV+pz6cnkbg4e6GDy2Se4qFKI8WmsmrbP4gx882U2Gg9aHZj6K+QNN/rNf87cf5KRQuL9Xdtg0bgC3d6j/nOKoyw+LjY0lKyuLwMBAsrKyGDp0KM2bN+f6668/5zH9+/cv8XjKlClMnTqVw4cP13Zzq+yhcPgqw48bVlqsujwfX7Nm3teVY65NErdnWvirg1W/+bKoTz4ncy1Ocv6Yg4KC6qZxXqRfsMGjXTVPJtlc0UbRo2X9/8ckRFXkFmpm/2xzZyeDJg3rv0e8RIHVtm1bsrOzcTqdOBwOtNZkZGQQGhpa4qDQ0FDS0tKKH6enp5fapywBAQHFfw8JCeHGG28kMTGx3AJr7dq1JQbL+/j44OPjA7h2sv1gsKbbx05mJgfyau8aqrBw7Zhrk8TtWQ6e0Pxlu5ObIxTxFwWWeM1TY3ZlT3YzWLlXc9Nqi63XKZf4z0mIypr7i02+E+7p7BonCSVa0bJlS2JjY5k/fz4AixcvJiQkhIiIiBIHxcfHs3TpUvbt24fWmlmzZjFmzJjzflh2dja2XdQNfezYMZYtW0ZMTEy5x/j7+xMQEFC8nS6uXF2nZooXehq8ttPmiwzpehfiTA9tLJoe4B+X1tzJh6i601M37D8B98kC9sINnbI1L/1kc1OD5mgPAAAgAElEQVSEok1j1zhBKFXmzZ49m9mzZxMVFcXMmTOZN28eAHfccQdLly4FoH379syYMYPevXsTERFBixYtmDhxIgD5+fmEhIQwevRodu7cSUhICI8++ihQVLB16dKFrl27cumllzJ48GDGjx9fV7HWubs6GQxtqxj/rcWBEzIeSwiAlVk276VqXuhp1sv6YKJsEYGKV3uZzEvWfLRbTgqFe1m0W5ORBw92cZ2TNqW1a84nsGXLFrp160ZSUhKxsbGlXs/JyXGLSwn78zVdFjvp2VKx9EqzWnNyuEvMNU3i9hz5Tk2XRU5C/RWrrin9++AuMXtKfjqb1pobvrFYuVfzQ7yDUP/K5St3jbs6vDFmcK24tdZ0/8TJBT6KFUNrd2h5ZeJ2jQuVHqyVn+Lf/UyWZWhm/SxnhcK7Pb3FZm8+zO5TvZMNUTtOL2DfpAH8aY1M3SDcw9pszZZD9bssTllcqzUe6towgzsvMnhoo83PhyVhCe/04++av/9o80S0QVQtLL5akVUoyltpAs69QoU3aeajeG+gybfZmud/lJNC4fpe3G7TJQgGt3GtkzYpsOrIC5cahDWBm1c7KbSkyBLexbI1E9ZZdAiEKV1rJ+1UZBWK8laaKG+FCm/TP9jgkWiDv2y22XRQiizhujbst1mWoZlyiev1ikuBVUf8HIoFAx38dBimbZaEJbzLGz/bfHdA82Zfk4a1sB5ORVehKG+lifJWqPBGM7oZxDRX3LTKIu+UnBQK16O15qGNNtEXwI3hrlVcgRRYdSqmueLZ7gZ//9Fm9W9SZAnvkJWneWyTzcSOBr1b107KqegqFOWpygoV7rbSRGWcnrohOx/uXS9TNwjXs3iPZsOBojuSzXpeFqcsdTqTu4CHLjFYnqm5ZY3Fj6MUQY1c74dCiJp0z3qLxg6Y2cPzzufccaWJyrgA+GtXB/cm+RLTJJeb2znL3d9T4q4Mb4wZ6j/uQhv+vLExg1vbxPgeIyenbj63vLjPvrtQCqw6ZijFOwNMui5xMmmdxQeXu951YyFqysd7bD75VfPR5SZNfWrv57yiq1CUpyorVLjrShOVcU8Q/JDn5M/bfOkd6iC2efnfo6fEXRneGDPUb9wvbbfIyLf5/OqGBDXzrdPPlmkaXFhbf8XsPiYf7dG8myJjG4Rnyi3UTF5vcW2oIv7C2j2JqOgqFOWpygoV7rrSRGW91sukSzPFqJVOfj8pOUvUr8MFmqe22iR0MOjUzHU7KKTAqiej2xuMi1JMXm+RlisJS3iexzbZHC2Ef/Wum17aiqxCUd5KE+WtUOHtGjkUi64wyTsFN6+W+bFE/Xpmq80pG6Z3c+0SRi4R1qNXLjP5NtvJ2NUWicNMHC44SE+Iqtiw3+b1nTb/uNSo9GzgVdWhQwc2bNhQ6vk5c+YU/93Pz4+srKxzvkdCQgIJCQm10j53F9ZEsXCQyZAvLaZvsXm6u+ssSSK8x+5czas7bP4SY9Daz7X/z3Tt8s/DNWmomD/QZNNBLVM3CI9xtFBz82qLuBbKZVa1FzVjcIjBM90Nntlq89mvkrNE3Xt0k0WLRvCgi83aXhbXb6GHu6yVwcweBjN/sPkwTRKWcG9aayYkWvx+EhYOcs1bp0X1PNLVYGSY4pY1FqlH5VKhqDsb9tt8uFvzbJyJn8P1c4sUWC7goS4GN4Yrxn9r8cPvkrCE+3prl+bD3Zo5/UzaB7h+AhSVp5Ti7QEmrXzhupVOjsskpKIOaK15+LuiSUVviXCP3CIFlgtQSjGnn0mHQBi5Qu7SEe5pe47mvg0WEzsajG4vqcWTBTZULLnCwZ5jMCHRQmvJWaJ2LUnXrN/vupOKlkWyoIvwcyg+HuwgzwnXf2PhlLt0hBs5fkpzwzdOIgPhn5dJWvEGnYMU/+5nsiCtaNCxELWl0NJM/d5iaFvF5W3cJ7+4T0u9QFgTxUeXm6zN1kz5ThKWcB/3rrf4NQ8+GOTA1w3GRoiacX24wYNdDB7aaLNun+QsUTte32mz5xg838O97lyVAsvFDPiDwUuXGfzzJ5t3kyVhCde3INXm38ma13qZXOTCk/6J2jGzh0GvVorRX1vsOyHfv6hZpycVvaODQecg9/r5kgLLBd3dyWB8lGLCOovNB6XIEq4r5ahm4jqLmyMU46LcK/mJmtHAUHx4uYmh4PbvGnFKhjeIGvTsfycVneHik4qWpVSLU1JS6NWrF1FRUcTFxbFjx44yD5w7dy6RkZGEh4eTkJDAqVOngKL1uwYMGEBgYGCJNbrOd5z4H6UUr/c2ib5Acd1Ki/35krCE6ymwNGNWOQn2hTfqaLZ24Zpa+RXN9J6UY/LwRjkpFDXj9KSiUy9x/UlFy1KqwJo4cSITJkwgOTmZqVOnMm7cuFIH7dmzh2nTppGYmEhqair79+/nzTffBCAgIIBnnnmGBQsWVOo4UVIjh2LxFSZOG/74tUWhJUWWcC1Tv7f5KQc+uNxBk4bul/xEzbqslcEzXQt4ZYfNglQpskT1PbbJormbTCpalhKtPnDgAJs3b2bs2LFA0eKnmZmZpKamljho0aJFDB8+nNatW6OUYtKkSSxcuBAoWmW6T58+NG7cuNSHlXecKK1NY8XiwSbfHdTcv0ESlnAdS3+1efknmxd6GsQ0l+JKFLm9/SluiVAkJFpsz5GTQlF1G/fbfOBGk4qWpUSBlZmZSXBwMA5H0RKFSilCQ0PJyMgocVBGRgZhYWHFj9u1a1dqn7JU5bi8vDxyc3OLt4KCgvNH5UF6tTJ4vbfJGz/bvLunQX03Rwgy8zTj11qMDFNMlqVwxBmUgll9TSIDiub0O3hCiixReVprHvrOpmuQ+0wqWhaXX+y5f//+JR5PmTKFqVOncvjw4XpqUd0b1RLWt/dhylYfOgYcoccF3tWb5U3f9ZlcMW6nDaO/9cXPNHjhkuPUdBPPF3NQUFDNfqCocX4OxSdXOrj0UyfDVlisusZ9eyBE/Tg9qejKoe4zqWhZShRYbdu2JTs7G6fTicPhQGtNRkYGoaGhJQ4KDQ0lLS2t+HF6enqpfcpSlePWrl1bYrC8j48PPj4+gHcl21kDNLtyTzL+u8Zsvs5Bm8bu+0NXFd70XZ/J1eJ+YpPF5hybb4eZhLeqnba5Wsyi8to1USy7ymTAMosbV1ksucK9/6MUdefMSUWvcKNJRctSovUtW7YkNjaW+fPnA7B48WJCQkKIiIgocVB8fDxLly5l3759aK2ZNWsWY8aMOe+HVeU4f39/AgICirfTxZW3aWgq/t3zJA4DRq20OOmUrndRt77ea/PcNpunuxv0auXeiU/Uvu4tDD683OTzDM09621ZTkdUyBs/u+ekomUplSVnz57N7NmziYqKYubMmcybNw+AO+64g6VLlwLQvn17ZsyYQe/evYmIiKBFixZMnDgRgPz8fEJCQhg9ejQ7d+4kJCSERx999LzHifNr2UjzyWAHP+Zo7vyPrP8l6s7+fM3Y1RZXtFFM7SrFlaiYoaEGs/oUjSF9/gfvGtogKi8zTzM9yeb2DsrtJhUtS6kxWB06dGDDhg2ldpwzZ06JxwkJCSQkJJTaz8/Pj6ysrHN+4LmOExXTrYXirb4mt6yxuKipzZSu7l/lC9dma82f1lgAvDfAxJD5rkQl3NHRIPO45pFNNiH+ipsjpEAXpVl2UZ7xbwB/84DeK3CDQe6itLGRBr8c1Uz93qZpQ8WEiyRhidozPclm5V7NiqEmrdxwsj9R/6bHGmT89+7TYF8Y5OZja0TNe3G7zdpszTfXmDTz8Yw8IwWWm3qqm8HRQpi0ziKgIYwJl4Qlat6cXTZPb7X5a5zh9gNORf1RSvFmX5Pfjltct9Ji3XBFFw+4BCRqxpZDmic22/z5EoOBf/CcPOM5kXgZpRQvXWZwS6TiltUWn2fI+AZRs77IsJm0zuKuToaMuxLV1sAoWk6nfQBc/aWTrDwZQyog36m5ebWTi5vB0909K894VjRexlCKuf1Mrg1V/PFrizW/SZElasamgzajv7G4NlTxymWGrDMoakSThoovhjgwVVGRdbRQiixv9+fvbNKPwf8NdNDQ9Kw8IwWWm3MYivcvN+nTWjFshcWmg1JkiepJy9Vc86XFJUGKBYNk/iJRs4L9FF8OcZB1HK5bIeuserNlv9q8vtPmxUsNLmrmeXlGCiwP4GMqPhls0qWZYshyix2yBpioooMnNFcvd9K0IXx2lczALWrHRc0US680+c9+zW3fypQz3mh/ftF3f01bxZ0eeqOWZ0blhRo3UHw+xKRtYxi83MnuXElYonLynZphKyyOnoIvr3bQvJEUV6L29A02eG+Ayf+lah7fLD3v3kTrouIKYG4/02OHIEiB5UGa+Si+utqBvwOu+MLJ3uNSZImKsWzNjasstudoPr/KpH2AZyY84VquDzd4safBX7fZvLHTqu/miDryxs82X2Rq/t3Ps6d+kQLLw7TyU3x9jQOnDVd+4eTQSSmyRPm0LlrK5PMMzUeXm3RvIWlB1J0Huhjcd7HB5PU2S3+VnixP9/NhzUMbbe68yODaMM/ONZ4dnZcK9Vd8PdTBwZMwZLlFrtypI8ox8webN362md3HZGio+6aElJQUevXqRVRUFHFxcezYsaPM/ebOnUtkZCTh4eEkJCRw6tQpANasWYOvry/R0dHF24kTJ+oyBK+klOLFngbXtVNc/43FFzLljMcqtIqmZGjXBF641H1zTUV5foReKqqpYsVQB6m5mmFfWZyQxaFFGd5LsXlsk82TsQa3d3TvdDBx4kQmTJhAcnIyU6dOZdy4caX22bNnD9OmTSMxMZHU1FT279/Pm2++Wfx6hw4d2LZtW/Hm6+tbhxF4L9NQ/N9AkyEhipErLT5NlyLLE03bbPPTYVgw0OEVN9C4d0YV5Yq+QPHFEJPNhzR//FpuhxYlfb3X5ra1FrdFKZ6Mde9UcODAATZv3szYsWMBiI+PJzMzk9TU1BL7LVq0iOHDh9O6dWuUUkyaNImFCxfWR5PFWXxMxUdXmIwIK5rX76PdUmR5kjW/2fz9R5tnuhvENPf84gqkwPJ4vVoZfDLY5Ou9RQtpWrYUWQK2/a4ZtdLiijaKWX3d/y6ezMxMgoODcTiKVv9SShEaGkpGRkaJ/TIyMggLCyt+3K5duxL7pKWlERsbS1xcHK+//vp5PzcvL4/c3NziraCgoIYi8k4NDMXCQSbXt1eMWWWxIFWKLE9wuEBzyxqL/sGKh7p4T9khaxF6gcEhBgsHwehvLGxt8e4Ak0Ze0D0rypaRpxn6pZPIQPjoCpMGMpEoALGxsWRlZREYGEhWVhZDhw6lefPmXH/99ec8pn///iUeT5kyhalTp3L48OHabq5Lqqm4X+oK2tmIsas1Obknuamds0betzbId10+reGO7xtxrNDBK9F5HD3i3if55cUdFBRU4rEUWF5i1IUGiy6Hm1ZbXLXc4pMrPWfFclFxOSeLJhL1MeDzqxz4N/CMn4G2bduSnZ2N0+nE4XCgtSYjI4PQ0NAS+4WGhpKWllb8OD09vXifgICA4udDQkK48cYbSUxMLLfAWrt2LdHR0cWPfXx88PHxAUonW29RU3HPH6wJWGdzT5IvDX1NJrjwZJTyXZ/bu8k2n2RZvD/IpEtIszpoVe2r6Pftuj+xosZdd6HBN9eY/HRY03upkwxZbNWr7D2u6b/Myf4TRROJtvag+WdatmxJbGws8+fPB2Dx4sWEhIQQERFRYr/4+HiWLl3Kvn370Foza9YsxowZA0B2dja2XXRJ6tixYyxbtoyYmJhyP9ff35+AgIDi7XRxJarPUIo3+hjc09lg4jqL13bIPFnuZk+uZvJ6i1siFDeEe1+54X0Re7lerQzWD3dwwgmXfurkh9+lyPIGvxwpKqqPFsK64Q46NPWc4uq02bNnM3v2bKKiopg5cybz5s0D4I477mDp0qUAtG/fnhkzZtC7d28iIiJo0aIFEydOBIqKsi5dutC1a1cuvfRSBg8ezPjx4+stHlE0lu7lywwevsTgnvU2L/4oRZa7OH5KM2aVxQU+8Fpvs76bUy9KXSJMSUnh1ltv5dChQwQGBvL222/TuXPnUgfOnTuXmTNnYts2gwYN4vXXX6dBgwblvrZmzRquvvpqOnToUPw+GzZskFuh61iHpooNIxxc85WTvp85WTLY5Io2Umt7qu8P2Az90qKVL3x1tYMQf88rrqBoioUNGzaUen7OnDklHickJJCQkFBqv8mTJzN58uRaa5+oGqUUz/cw8DHg4e9sCix4LMY7/8N2FyedmpErLXYe0ay6xiSgoWfmnPMp9b9qdeeSkXlm3ENrP8Waaxz0aqW4ernF/BS5W8cTfZVpM+hziw5NFYnDPLe4Ep5NKcUzcSZPdTN4fLPNk0myQLSrOmVrblhlsW6fZtlVJnFevDJEichrYi4ZmWfGfTRpqPjsKpNbIhW3rLGYuU2SlidZkGpz7VcWA/+gWDnUJEgWbxZublqsycw4g6e22Dy6yZZ85WIsW3PrGovlmZolg036B3tvcQVnXSIsby6ZMweLljeXTEXnmTFNk/Hjx3PXXXeV28DT88ycduZdOqL6GhiKuf1M2voXJayMPHi1l4Ept+67tZd/srh/g824KMWbfWUqBuE5pkabNDThwY1Flwv/canh9vO4eQKtNZPWWXywW/Ph5SZXt/Xu4grqeJoGmWememoz5vsuhGY04OGtPqQfKeDNHifxc5FJPLzxu4aqxa01PL2jIS//4sO9UQX85eJCjh2phcbVkvPF7K23w4uSHuhi4mPC3f+xyXfCa70NOYmoR1oXLeA85xfN2/1N4i+U4grOKrBqYi4ZmWemdtVmzPcHQWQLm+u/UYze0JDPrjJp7iKXlbzxu4bKxe20NRMSLeYla/5xqcEDXfxrsWW1x1u/a1E5d3Uy8TUVExItdh0p6jVp5UFTj7iT6Vts/vmTzb96G9waJcXVaSX+JWpiLhmZZ8a9XRNqsOZak7RcTa9PnezOlTEO7iDfWbT0zXspmvkDTR7oIndZCc83voPB6mtNko9qun3i5PsDcrNOXXvhR4unttjMjDO4q5PknTOVKjWrO5eMzDPj/uJaGGwY4UADly118m22JC1XdrhAc+UXFt/8pvnsKpObI+QMUniPPq0Nkq5zENpY0fczi7m7JF/VlVk7Lf78nc3j0QZTo6W4OpvSLnobxpYtW+jWrRtJSUnExsaWej0nJ8frLiXUdcwHT2j++LVF4j7N1K4GM7oZNDTrvgveG79rqFjce49rrlruZF8+fD7EpGdL9y6u3OW7lvxUtvqMu8DS3L/BZtbPNhM7Grzcy8CnDvKVt37Xs7blctcmX+7pbPDSZd5zo0Flvm/3zsaiVrXwVay6xuS5OIMXfrS59FMnPx92yXrcK237XdNrqZNjp4pmZ3f34kqI6vAxFW/0MZnT12Ress2AZRZ7j0u+qg0f77GZvLkR46MU//Si4qqyJCOLcpmG4pFok+9GOjhhQezHTv61Q+bLqk+WrfnrNosenzgJ8oH1wx109MClb4Soits7GiQOM8k8run2sZN1++SSYU36KtPmhlUWw9s4ebOviSHF1TlJgSUqJLa5Iuk6B7d3MJi83uaaryz25UuRVdfScjX9llk8vsnmwS4GG0c4aNNYEpwQZ+rR0iBppIMOgYqBy4oWipaTwur7NtvmupUWV4Uo3og7KfMlnocUWKLC/ByK13qbfH6VyZZDmi6LnXyaLmeHdUFrzVu7bLoudpKdr/l2mMnMHmadjDERwh218lN8fY3J3Z2LFooet9bihFOKrKradLBoZYjLWik+utykgVQP5yX/RKLShoYabI930LuVYuRKi4RvneSdksRVW/bla4avsJiQaHFjuOKHUQ76tJZfXSHOp4GheOkyk/cGmHy0W9PnMye/HpNcVRlaa/79S9GYts7NFJ9eadLIISd2FSFZWlRJC1/Fx4NN3uprsiBNE73Eycb90ptV0z7eY9NlsZPvD2o+vdLkrX4OmnjpyvRCVNXYSIP1IxzknIRuHztZkSW5qiIOF2jGrLK4/duik7uVQ038G0j+qSgpsESVKaW4o6PBtlEOLvBR9PnMYkaShdOWM8TqOlqouXtzI0Z9bdGnleKneAfDw+TXVYiqir5Asfk6B7HNFVctt7hplZPf5C7Dc0rMtole4mRFVtEs+XP6OaS4qiTJ2KLaIgMV64abPBFj8NRWm76fWWw9JImrqtb8ZnPJYifL9jqY199kyWCTFr6S2ISorgsaKb682mRef5NvftN0+MjJ33+wKLQkX53mtDVPJlkM+NwizF/xQ7yD0e2lVKgK+VcTNaKBoZjezWTdMJPDBZrYj5388WsnO3IkcVXUSafm4Y0Wgz63aOevSLziOOOiZI4ZIWqSoRTjogx+Ge3gtiiDRzbZdF3i5Ou9ctlwT66m32cWz261mR5rsPoak1B/yT9VJQWWqFGXtTL46Y8O/t3PZPPBojsNx652knpUCq1z0VqzIssm7hMnr+6web6nwaprTEIby7+ZELWlqY/i5V4mW69z0KKRYvAXFqO/dpKZ552/dwtTiy4Jnr5LeVqsKdMwVJMUWKLGOQzF+A4Gydc7+Fdvg9W/aTp+5CThWycZXpq8ymJrzSfpNj0+sbhquYWfo2iMyMOXSGIToq5ccoFi7bUm8wearNtXlKue22pR4CWXDY8Vam5d4+Sm1RbXhCq2xTvo1UpKg5og/4qi1jQ0FXd2Mkm9wcHfexp8+qsm8gMn9/zHItuLJyl12poFqUXjrK5badG4Aay42mTjCJMuQVJYCVHXlFLcHGHwy/UO7rzI4Mkkm4sXOVme6dmXDb8/YBPzsZMl6Zp3B5j830CTQLlLucZIgSVqna9D8UAXk91jHEzvZjA/1Sb8fSdTvrM4dNJ7Cq1CSzNnl03Hj5zcvNoi1F+xbpjJmmsdDA6RsVZC1LeAhooXLjX5Id5BqL9i6JcWI1c42ZPrWXnq9HJbvZdaXOCj2DbKwS2RkoNqmhRYos74N1A8Gm2yZ4yDhy8xeONnmwvfd/KXzRZHCjwrgZ0p36l55SeL8A+cTEi0iL6gaNmhL4Y46C0Thgrhcjo1U3w91OTDy02SDmk6LSo6IUw+4t556oRT826yTc9Pi5bbmtLVYN1wk/AAKaxqg6O+GyC8T1MfxVPdTe692OD5H2xe+NHm1R02o9sr4tsZDGqjaOABY5ByCzX/2mnzz+02OQVwc4Tika4mFzVz/9iE8HRKKUa3V1zdVvHXbTav77T5+482vVopxkcZXN9eEeAml9NSjmpm/WzzdnJRLhrcRrHmWpN+wXKCV5ukwBL1pnkjxfM9TR7sYvDyTzYf7bF5a5dF04YwIkwRf6HB4DbukcBO01qz8zC8v7uoaDzhhNs6GEy5xOBCOUsUwu34N1A8G2cyLaZoHOm8ZJsJiRb3rof4CxXxwSbDm2kMF7u85rQ1S3/VvPGzzdd7NUE+MD7KYOJFBpGBrtVWTyUFlqh3rf0Uf+1h8lycwY85sGiPzeI9Nu+kWPg3gCtbNeKmjjZDQhSNXXAm4QMnNF/v1azca7MiS/NbPvg5YNJFBg91MfhDY9drsxCicho5FDeEK24IN8jK07yXajPvF5v5qX6EbXVya5TBrZEG7ev5RGrv8aKF4d/aZfNbPlzWUvHuAJPRFypZQ7COleofTElJoVevXkRFRREXF8eOHTvKPHDu3LlERkYSHh5OQkICp06dqvZrwrsppeh6geLp7iY7Rzdgxx8dTLnEIPmYwR+/tmjxnpP4lU4WpNrkFtbfWIiTTs03e20e+d4idskpWs0vGrSedEhzY7jBV1ebHLrFwYuXmlJc1aHazl1CnBbiXzSe9JfrHSwfcJzBbRT/3G4T/oGTAcucvJNsk3eq7nKU/d+59K5b4SRsoZMXfrQZFlq0jNn6EUUD2KW4qgf6LAMHDtTz5s3TWmv90Ucf6e7du5+9i969e7cODg7W2dnZ2rZtPWzYMP3aa69V67WzJSUlaUAnJSWV+frvv/9e5vOezBtj1roo7pQjtp651anjPj6lebNQN5xTqK9efkpP2ejUs3Y69VeZlk4+YusCp13jn2/btt7+u63/8aNTD/nilPadW6h5s1C3eq9Qj111Sr+bbOnfjtf853rj912dmGszd51N8lPZvDHu0zHnFdr6vWRLD1pWlKMa/7tQD/vylH5wg1P/a4dTf5lh6ZQjti60qpcr8k/Z+ufDtl6eYek3djj1nzc6dcT7RTnp4o8K9b92OPXRgprPR2fzxu9a68rFXeIS4YEDB9i8eTMrVqwAID4+nsmTJ5OamkpERETxfosWLWL48OG0bt0agEmTJvHcc89x9913V/k1IcoTEaiYGm0yNRp+PaZZkm6zPFPz4W6bzONwek5AQ0EbP2gfoLiwCVzYRNG+yf/+3toPCizIKYDDBZBToMkpOP247L9nHdfsPwE+JvRrrXiqu8GVbQy6BCG3NbuI2s5dQpxP4waKsZGKsZEG6cc076bYrN+vWfqrTfoxcP43R5kKwvwhPED9d4PwJv/7u6Hg1zxIP6aLttN//++f+0/87zNNBW0bQ+/Wirf7G/RqpSQnuZASBVZmZibBwcE4HEVPK6UIDQ0lIyOjRJLKyMggLCys+HG7du3IyMio1mvnkpeXR25ubvFjHx8fAP72t7/x1FNPFT/2dAUFBV4XM5Qdd1iTonm1HuhStM8pW5OZB3uOafYcg93HNHuOaX4+Ap9n2Bw8+b/3MxTY5+i5b9IAgnyKtmY+iiCfoqLs6raKfq0VfVorfOuom90bv+/qxFzbuetcJD/9j/zM/i/mdk0Uf4k1ix87/5uj0o5pUo9q0o5BWq5m/X6b91LguLPs9zcVhPpDO39Fp6ZwdYhBuyaKdk2KnmvTuGjljLrmjd81VD5ulx/k3r9//xKPp0yZwp133pUd0VYAACAASURBVMnzzz/PnXfeSUBAQD21rG7l5uZ6XcxQ8bibAjG+RRstS76W54TM4wbpxxX7Tho0dmiaNdA0bVi0NWsIgQ00Dc5zx/KJXDhR/i41xhu/74rEHBQUVMetKp/kp/+Rn9nyYw4EYn2LNlr/73mt4WCBYk+ewZ7jCg2E+WlCG9u0bqRxnCsvnYLcIzUUSCV543cN54/77PxUosBq27Yt2dnZOJ1OHA4HWmsyMjIIDQ0tcVBoaChpaWnFj9PT04v3qepr57J27Vqio6OLH/v4+FBQUFAcjLd8uafPzL0pZqiZuIOAUKB3zTWr1nnj912dmGs7d52L5Kf/kZ/Zqsd8AdCxhtpUF7zxu4bKx12iNm7ZsiWxsbHMnz8fgMWLFxMSElKiix2KxjcsXbqUffv2obVm1qxZjBkzplqvnYu/vz8BAQHFmzd1RwohKqa2c9e5SH4SQpxLqUuEs2fPZty4cTz33HMEBAQwb948AO644w6GDx/O8OHDad++PTNmzKB376I+gQEDBjBx4kSAKr92thMnii7G/Pzzz6Vey8vLA2Dbtm34+/tXPXo34o0xg8TtTXFXNOaOHTvi5+dX6vnazF1nk/xUmjfG7Y0xg8RdXtwl8lNt3MZYE+bPn68B2WSTTbYS27mmRpD8JJtsstX3dmZ+UlprjQs6dOgQX331Fe3atcPX17e+myOEcBHn6sGqS5KfhBBlOTM/uWyBJYQQQgjhrmQpbSGEEEKIGiYFlhBCCCFEDXOpAstbF2utSNyrVq2iR48edOrUic6dOzNlyhRs2waK5usxTZPo6Oji7cz5fFxRRWJes2YNvr6+JeI6ffcWeO53PW/evBIxN2/enFGjRgHu+V3fe++9tGvXDqUU27ZtO+d+7vB77Y05yhvzE3hnjvLG/AS1mKPq6SacMtXlYq2upCJxb9myRf9/e/ceF1W57w/886w1OF5wUEKMBCRULDBBFDVRQUslb2Sal45t2abSxbMvnn6w2zuPu7Z5zP3r/OrYaYflce+ynSZoorRTSyRNzbsZxuYSApqIhgooF2et5/fH5BwRxREGmMvn/XrNq2bNmjXP19U8flrrmecpKCiQUkpZXV0to6Ojre8pLCyUXl5erdVcu7Cl5szMTBkeHn7L97vyub5ZWFiYTE1NlVI657nOysqSJSUlsmfPnvLo0aO33MdZvtfu2Ee5Y/8kpXv2Ue7YP0nZcn2UwwSsc+fOyc6dO8tr165JKaXUdV12795d5uXl1dtvxYoVMjEx0fo8IyNDRkdH3/E1R2Vr3Td74YUX5JIlS6SUzvcfta01N9Z5ucu53r9/v+zWrZusq6uTUjrfub5RY52XM3yv3bGPcsf+SUr37KPcvX+S0v59lMPcImxssdYb2Xux1rZma903Ki0tRWpqKiZOnGjdduXKFURFRSEyMhKvvvoqNE1r8bY31d3UXFBQgMjISERFReGdd96xbneXc7169Wo8/fTT8PDwsG5zpnNtK2f4XrtjH+WO/RPgnn0U+6fGNeV77fCLPVN9FRUVmDRpEpKSkjBo0CAAgJ+fH86cOQNfX1+Ul5djxowZeOONN5CUlNTGrW2eyMhInD59Gl5eXjh9+jTGjx8PHx8fTJ8+va2b1iquXLmCdevWYf/+/dZtrnquyTW4U/8EuHcfxf7pzhzmCtaNi7UCaHSx1qKiIuvzmxdrvd1rjsrWugGgsrIScXFxiI+Px6JFi6zbjUYjfH19AVgWoZw7dy52797dOgU0ga01m0wmeHl5AQD8/f0xa9Ysa12ufq4BYMOGDQgLC0NoaKh1m7Oda1s5w/faHfsod+yfAPfso9g/Na5J3+uWuZPZNDExMfUG2A0cOLDBPgUFBQ0Gk61cufKOrzkyW+qurKyUw4YNk6+88kqD186dO2e9B15TUyOnTZsmFy9e3KJtbi5bav7xxx+lpmlSSikrKirksGHD5OrVq6WUrn2urxs+fLh877336m1zxnN9XWPjG5zle+2OfZQ79k9Sumcf5c79k5T276McKmDl5OTIoUOHyj59+siBAwfKb7/9Vkop5TPPPCM3b95s3W/VqlUyODhYBgcHy7lz51pP6J1ec1S21L106VJpMBhkeHi49bF06VIppZRpaWkyLCxM9u/fX4aGhsqFCxfKmpqaNqvHFrbUvHLlShkaGmqta8mSJVLXdesxXPVcX9/P09NTVlRU1Hu/M57rBQsWyB49ekhVVaWvr6/s1auXlNI5v9fu2Ee5Y/8kpXv2Ue7YP0nZcn0Ul8ohIiIisjOHGYNFRERE5CoYsIiIiIjsjAGLiIiIyM4YsIiIiIjsjAGLiIiIyM4YsKjJdu3aBSEE/vjHP7Z1U4iIiBwKAxY5rKCgIAQFBdnlWHl5eVi2bBlGjhyJ++67D+3atUNAQAB+8YtfICcnxy6fQUREdB3XIiS3sHjxYqxfvx79+vVDfHw8TCYTTpw4gQ8//BCpqan4/PPPMXLkyLZuJhERuQgGLHILcXFxSE5OxoABA+ptX7duHWbNmoXnnnsO2dnZbdQ6IiJyNbxFSHaxZ88exMbGonPnzujSpQumTp2K/Pz8BvuVlZXht7/9LXr37g2j0QgfHx9MnToV3333nXWfU6dOQQiBoqIiFBUVQQhhfVwf71VXV4eVK1di3LhxCAgIsC4y+sQTT+Do0aMNPjchIaFBuAKAmTNnIiQkBCdPnsSFCxfs9wdCRERujUvlUJPt2rULo0aNwrhx45CZmYm4uDiEhYUhOzsbW7ZsgY+PD/bv34/g4GAAQEFBAWJjY3H69GmMHTsWDz30EMrKypCWlgYhBL788ksMGTIEly5dwptvvok333wTAPCb3/zG+pmxsbGIjY1FaWkpevTogREjRqBv377o2rUrfvjhB6Snp0MIga+++gpRUVE21dGvXz9kZ2fj4sWL6NKli/3/oIiIyO0wYFGTXQ9YAPDuu+8iMTHR+lpKSgqeffZZTJw4EVu2bAEAREdH45tvvkFGRgbGjRtn3Tc3NxeDBg1CUFAQvv32W+v26wPcT5061eCza2trceHCBfTo0aPe9uzsbAwdOhRDhw7Fjh077ljDgQMHMGTIEERFReHAgQM2105ERNQY3iKkZgsJCcH8+fPrbZs/fz769OmDjIwMnD9/HkePHsXevXsxZ86ceuHqxvefOHGi3q3CxhiNxgbhCgDCwsIwatQofPXVV7h27Vqjx7h8+TLmzJkDRVGwYsUKmz6XiIjIFhzkTs0WHR0NRamf1RVFQXR0NPLy8nD8+HHk5eUBAM6dO3fLebOuT5WQk5ODfv362fS5x44dw4oVK7Bnzx6UlpY2CFQXLlyAn5/fLd9bXV2NKVOmICcnB6+99hpiY2Nt+kwiIiJbMGBRs3Xv3r3R7ZcvX0Z5eTkAICMjAxkZGbc91pUrV2z6zL1792L06NEAgLFjx6JPnz7w9PSEEAKffvopjh8/jtra2lu+t6amBvHx8cjMzMRLL72E3//+9zZ9JhERka0YsKjZzp071+h2Ly8vmEwmAMDKlSuxcOHCZn/ma6+9htraWuzevRvDhw+v99r+/ftx/PjxW76vuroa8fHx2LFjB5KSkrBs2bJmt4WIiOhmHINFzfb1119D1/V623Rdx969eyGEQHh4OIYMGQIA2Ldvn83HVVUVmqbd8rWCggJ4e3s3CFdXr17FkSNHbvmeG8PViy++iNdff93mthAREd0NBixqttzcXLz33nv1tr333nvIzc3FhAkT0K1bNwwePBhDhgzBxx9/jPXr1zc4hq7ryMrKqrfN29sbFy5cQE1NTYP9e/bsiYsXL9abHFTTNLz44os4f/58g/2v3xbcsWMHFi1ahD//+c9NLZeIiOiOOE0DNdmd5sG655578M0331jnwSosLMSoUaNQVFSEoUOHIjIyEh06dEBxcTH27duH8+fP1wtTycnJWLFiBWJjYzFixAi0a9cOI0eOxMiRI7F161ZMmjQJXbp0wfTp09G+fXvs2rULZ86cwUMPPYRdu3ahsLDQOtVDQkIC/va3v+Hee++tN53EjRISEuy29iEREbk5SdREmZmZEoBcsmSJ3L17t4yJiZGdOnWSJpNJTpkyRebl5TV4T3l5uXz55Zdlv379ZIcOHaSnp6fs06ePfOqpp+TGjRvr7VtZWSnnz58v/fz8pKqq1s+6LjU1VUZGRsqOHTtKHx8fOX36dFlQUCDnzJkjAcjCwkLrvjExMRJAo4/MzMwW+pMiIiJ3wytYRERERHbGMVhEREREdsaARURERGRnDFhEREREdsaARURERGRnDFhEREREdsaARURERGRnDhuwri95cvXq1bZuChEREdFdcdiAlZOTg4EDByInJ+eWr1++fLmVW9T23LFmgHW7E3esmYhck8MGrDu53SLArswdawZYtztxx5qJyDU5bcAiIiIiclQMWERERER25pQB63cHNCQfM4LLKBIREZEjcsqAFXmPwPsF7fCfJ/S2bgoRERFRA04ZsKb3UvDrvrVIOqBjWwlDFhERETkWpwxYAPCHsDo85i8wY6eG3Eu8VUhERESOw2kDliqAj0ar8OsITN5uxuU6hiwiIiJyDE4bsADAq53A5jEGlFYDT+3UoOkMWURERNT2nDpgAUBIF4H1o1V8flriD4c4HouIiIjantMHLAAYF6BgxWAFrx/X8XE+QxYRERG1LUNbN8BeFj2k4PhPEnO/0hDiJTCwm2jrJhEREZGbcokrWAAghMCqESr6ews8vsOM0qscj0VERERtw2UCFgC0NwhsGqNCk8DULzTUagxZRERE1PpcKmABwH2dLCHr0HmJF77WuJwOERERtTqXC1gAMMRXwaoRKlb/U+LtbA56JyIiotblMoPcbzYnxDLo/bf7dYR1FRjdwyWzJBERETkgl04dK4YoGH2fwJNfavihgrcKiYiIqHW4dMAyKALrH1HhbQTit5tRyeV0iIiIqBXYLWDl5eVh2LBhCAkJQVRUFLKzsxvss2/fPkRERCAiIgJhYWFITExEbW2tvZpwS12NAuljDSiqAp7epUHnoHciIiJqYXYLWImJiViwYAFyc3ORnJyMhISEBvuEh4fj4MGDOHbsGE6cOIGysjK888479mrCbT3YVeDj0SrSiyT+cJCD3omIiKhl2SVglZWV4dChQ5g9ezYAYOrUqSgpKUF+fn69/Tp27AgPDw8AQF1dHaqrqyFE68y4PiFQwf8domD5cR0f5DJkERERUcuxS8AqKSmBn58fDAbLjxKFEAgMDERxcXGDfU+dOoXw8HD4+PjAy8sLzz//fKPHrqqqQkVFhfXRnFuKv31IwTN9Bebv1vB1KUMWERERtYxWn6YhKCgIx48fR1VVFWbPno2NGzdi5syZt90/Jiam3vOkpCQkJyfj4sWLTfr8paHA9z91wOPbdewYdRWBnZxnTFZTa3Z2rNt93Klmb2/vVmoJEVHz2CVgBQQE4OzZszCbzTAYDJBSori4GIGBgbd9j6enJ2bOnImPPvqo0YCVlZWFiIgI63Oj0Qij0Qig6Z3t5sckhnxqxtPfeGLvZAM6t3OehaHd9S8Y1u0+3LFmInI9drlF6Ovri8jISKxduxYAkJaWBn9/f/Tu3bvefvn5+bh27RoAyxisTZs2oX///o0e29PTEyaTyfq4Hq6aw6e9wJZxBhRXAU9latB057mKRURERI7Pbr8iTElJQUpKCkJCQrB8+XKsWbMGADBv3jykp6cDAHbu3IkBAwYgPDwcAwYMQPfu3bF48WJ7NeGuhHYVWPeIis9KJF7iLwuJiIjIjoR00NWQjxw5goEDB+Lw4cOIjIxs8Hp5ebldbiW89Z2G3+zT8T8jVfyyr2PPu2qvmp0N63Yf7lgzEbkml12L0Fa/ClNw8iKQuEdDbxMwws+xQxYRERE5PrdPE0IIvB2tILq7wJQdXLOQiIiIms/tAxYAeCgCqY+q6GIEJm0zo4JrFhIREVEzMGD97J72AlvHGXDmKjBzJ39ZSERERE3HgHWDB7oIfPKIiu2nJf7PN/xlIRERETUNA9ZNxvorePNhBf/vOx3v5TBkERER0d1z+18R3soLoQqyLwLP79HQxwTE3sccSkRERLZjcrgFIQT+a5iCGD+BJ77Q8M9LHI9FREREtmPAuo3rvyy8twMw/nMzzlczZBEREZFtGLAa0cUo8FmcAVVmIH67hhozQxYRERHdGQPWHQR1FtgyVsXRnyTmZGnQHXNlISIiInIgDFg2GOyrYO0oFRt+kHj5EH9ZSERERI1jwLLR1PsVrBii4D+O6VjN6RuIiIioEZym4S7820MKCiqAZ/do6NkZeLQH8ykRERE1xIRwF4QQWDlMwaM9BKbu0JBdzvFYRERE1BAD1l0yKALrH1HRszMwYZsZpVcZsoiIiKg+BqwmMLUTyBhnQJ0OTN6u4SqnbyAiIqIbMGA1UYCnwNZxBmRflJidqUHTGbKIiIjIggGrGSJ9BNaNVrG5SCL5AH9ZSERERBYMWM00qaeCN4cqeOOEjr+c1Nq6OUREROQAOE2DHfxrPxX5FcDCvTp6egqMD2RuJSIicmdMAnbyn0MVTAgQmLFTw7GfOB6LiIjInTFg2YmqCPx9tIoQL2D852YUVTJkERERuSsGLDvy9LBM39BeBeI+N+OnGoYsIiIid8SAZWf3dhTY9pgBF2osc2RVc44sIiIit8OA1QL6eAlkjFNx7CeJWTs1mDlHFhERkVthwGohg30VfPKIiq3FEgu/1iElQxYREZG7YMBqQRMCFawaoSIlR8drRzkRKRERkbvgPFgtbG5fBWeuSCw+rOO+TgJz+zLTEhERuToGrFbw8gAFZ64AC3Zr6N7BcmWLiIiIXBf/pm8FQgj8d7SCiYEC07/UcKCMtwuJiIhcGQNWK7k+EWm4t8CEbRpyL3HQOxERkatiwGpFHQ0CW8ap8GlvmYj03FWGLCIiIlfEgNXK7mkv8HmcATUaMH6bGZV1DFlERESuxm4BKy8vD8OGDUNISAiioqKQnZ3dYJ+dO3di8ODBCA0NRVhYGJKSkqDr7jceqWdngX/EGZB/GZj2hYY6jSGLiIjIldgtYCUmJmLBggXIzc1FcnIyEhISGuzTtWtXrFu3DidPnsThw4exd+9efPDBB/ZqglMJv0dg0xgVmWcl5n2lcSJSIiIiF2KXgFVWVoZDhw5h9uzZAICpU6eipKQE+fn59fYbMGAAgoODAQDt27dHREQETp06ZY8mOKXRPRR8EKviw3yJpAOc7Z2IiMhV2GUerJKSEvj5+cFgsBxOCIHAwEAUFxejd+/et3xPaWkpUlNTsXXr1kaPXVVVhYqKCutzo9EIo9Foj2Y7hJm9FJRVS/x6n46u7YDfD1DbuklERETUTG0y0WhFRQUmTZqEpKQkDBo0qNF9Y2Ji6j1PSkpCcnIyLl682JJNbFWz7wPOhrbDHw4ZYTBfxbxe1265nyvVfDdYt/u4U83e3t6t1BIiouaxS8AKCAjA2bNnYTabYTAYIKVEcXExAgMDG+xbWVmJuLg4xMfHY9GiRXc8dlZWFiIiIqzPb7yC5Uqd7bJhEjWKjuRj7XFfl06Y3efWd29dqea7wbrdhzvWTESuxy5jsHx9fREZGYm1a9cCANLS0uDv79/g9mBVVRXi4uIQFxeHl19+2aZje3p6wmQyWR+udHvwRkIIvDFUwS9DBBKyNKQXud+vK4mIiFyF3X5FmJKSgpSUFISEhGD58uVYs2YNAGDevHlIT08HALz11ls4cOAANm7ciIiICEREROC1116zVxOcniIEVo1Q8XiQZUmdzB8ZsoiIiJyRkA7607UjR45g4MCBOHz4MCIjIxu8Xl5e7rK3Emo1icnbNOwtk/hyvIrBvpYc7Mo1N4Z1uw93rJmIXBNncndARlVg4xgV/b0FHvtcw3flDpmBiYiI6DYYsBxUJw+BjHEqAjoBY/9hxg8VDFlERETOggHLgXUxCmx7zABPA/DoZ2acrRZt3SQiIiKyAQOWg+veUeCLCQaYdWDq7g74qYZXsoiIiBwdA5YTCPQU2DHegJ/qLGOyKusYsoiIiBwZA5aT6NtFYMPwavzzksTk7RqqzQxZREREjooBy4n076IjI07FN2USM77UcE1nyCIiInJEDFhOZvi9CjaOUfH5aYnZmRrMDFlEREQOhwHLCcUFKPjkERUbCyWe3sWQRURE5GgYsJzU40EK1j+iYsMPEnN2adAYsoiIiBwGA5YTe+J+BetGq1j/g0RCFkMWERGRozC0dQOoeaYFK/gYwKydGgQ0rIlRoSqckJSIiKgtMWC5gCeDFegSeCpTgyI0rB7JkEVERNSWGLBcxIxelpA1e5clZL0/UoUiGLKIiIjaAgOWC5nVW4EE8PTPIWvVCIYsIiKitsCA5WKe6m25kvWLn0PWu8MZsoiIiFobA5YLmt1HgSaBX2ZpUKDjneEKQxYREVErYsByUXNCLLcL52ZpUATw39EKBEMWERFRq2DAcmEJIZbbhc98ZQlZK4cxZBEREbUGBiwXN7evJWTN320JWW89zJBFRETU0hiw3MC8BywhK3GPBrMOvB3NMVlEREQtiQHLTSx4UIFBsVzJqrgmsSZGhQcnIyUiImoRDFhuZG5fBZ09gH/J1FB5TcP60SraGxiyiIiI7I2LPbuZJ4MVpI9VseO0xIRtGqqucYFoIiIie2PAckNxAQq2Pabi4HmJRz/TUF7DkEVERGRPDFhuaoSfgswJBhRUSMRmmFF6lSGLiIjIXhiw3NjAbgJfTTTgpxpgxBYziioZsoiIiOyBAcvNPdhVYM8kAySA4VvMyLnEkEVERNRcDFiE+00CuycZ4NUOGLnFjKMXGLKIiIiagwGLAAB+HQWyJhoQ1FkgdqsZe0r1tm4SERGR02LAIqt72gt8OV5FpI/A2M80bCthyCIiImoKBiyqp3M7gc/iVDzSQ2DSdg1phQxZREREd4sBixroYBDYOEbFtPsFpn+pIeV7ra2bRERE5FS4VA7dkoci8GGsCp/2Op7do6OgAlg+mItEExER2YIBi25LVQT+a5iKXibgt/t0/FAp8WGsig5cv5CIiKhRdrtFmJeXh2HDhiEkJARRUVHIzs5usM+pU6cQGxsLLy8vRERE2OujqYX9up+KTWNU/KNEYlSGhnOc9Z2IiKhRdgtYiYmJWLBgAXJzc5GcnIyEhIQG+5hMJixduhR///vf7fWx1ErigxRkTVRRVCkxdLMZ319kyCIiIroduwSssrIyHDp0CLNnzwYATJ06FSUlJcjPz6+3n7e3N4YPH45OnTrZ42OplQ3qpmB/vAGeHsCwdDMyf+QvDImIiG7FLgGrpKQEfn5+MBgsQ7qEEAgMDERxcXGzj11VVYWKigrro7a2ttnHpKbr2Vlgz2QDorpZ5sr6Wy5DFhER0c0cfpB7TExMvedJSUlITk7GxYsX26hFbceRav5wMJB01IiErHb47txV/C60Di31A0NHqrs1uWPdd6rZ29u7lVpCRNQ8dglYAQEBOHv2LMxmMwwGA6SUKC4uRmBgYLOPnZWVVW9AvNFohNFoBOCena0j1fzXRyVCj+v43UEjfrzWHv8To8KotkzKcqS6W5M71u2ONROR67HLLUJfX19ERkZi7dq1AIC0tDT4+/ujd+/ezT62p6cnTCaT9XE9XFHbE0IgOULFJ4+oSDslMeYzDT/VcPA7ERGR3X5FmJKSgpSUFISEhGD58uVYs2YNAGDevHlIT08HAFy9ehX+/v548skncfLkSfj7++Oll16yVxOojTwZrCBzgorvL0k8vNmM/MsMWURE5N7sNgarb9++2LdvX4Pt77//vvXfO3bsiNOnT9vrI8mBPNxdwf54gQmfmzF0sxkbHlUx6j6uxERERO6JfwOS3fQyCeydbEDEPQKPfqbh9WMapOTVLCIicj8MWGRX3u0Ftj2m4qVwBb87qGPKDg2XahmyiIjIvTBgkd2pisDSKBVbxqrIOisx6FMzjv3EkEVERO6DAYtazMSeCg5PMcDkATy82Yy/clJSIiJyEwxY1KKCTQJfTzbgX3oL/DJLw4LdZtSYeTWLiIhcGwMWtbgOBoH3RxqweqSKD/MkoreYUVjBkEVERK6LAYtazdy+CvZONuBSLTDwUzM+K+YtQyIick0MWNSqBvgIHJpiwPDuAhO2afj3Qxo0nVeziIjItTBgUavrahT4dKyKZVEKXjum47HPNVzgEjtERORCGLCoTShC4KUIFdsfU3HsJ4kBG83Y9SNvGRIRkWtgwKI29UgPBUemGBDcWWBUhoZf7dVw5RqvZhERkXNjwKI25+8pkDlRxZsPK3g/R0fERjO+LuXVLCIicl4MWOQQFCHw634qjj1hQLf2AiO2aHhxv4ZqzpnlFqSU+MtJDX/+vl1bN4WIyC4YsMihhHQR2D1JxYohCt4+qSNykxmHyvmfqSu7apZ4epeG57/W4cFTTUQugt0ZORxVEXixv4ojUwzo7CHwWGZHvHRAQ63Gq1muprBCIjrdjE2nJD4ereI3fevauklERHbBgEUOK7SrwN7JKl4Kq8MbJ3QM3GTG4fMMWa5i+2kdAz81o6IO2DfZgJm92B0Rketgj0YOzaAILHqgDoenGNBOAYZsNmPJYQ11vJrltKSU+I9jGuL+oWGor2Xi2f73iLZuFhGRXTFgkVN4yFvgm8cNWDxAwbKjOoZsNuP4TwxZzqayTmLaFxp+f1DH7yMUbBmroquR4YqIXA8DFjkND0VgyUAV3zxugCaBqE/N+Lf9Gi7WMmg5g39ekhiy2YwdZyQ2jVGxNEqFqjBcEZFrYsAipxPpI3DwcQOWRCpI+V5H7/VmrPxOwzWuaeiw0ot0DP7UDAngwOMGPB7EroeIXBt7OXJKRlXgDwNU5M0w4IkggV/v0/FQqhlbi3RIyaDlKHQpseSwhvjtGkbfJ/BNvAEPdOFVKyJyfQxYuCnzKgAACytJREFU5NT8Ogq8N9KAo08Y0KOTwKTtGsZ8pnF8lgO4VCsxaZuGPx3RsXSQgrQxKkztGK6IyD0wYJFLCL9H4IvxKraMVXH6imXx6HlfmXH2KoNWWzh4XkfUp2bsPSeREafiDwNUKILhiojcBwMWuQwhBCb2VHBimgFvPaxg0ymJPuvNWHqES+60lh+vSCTsMmPwpxo8PYBDUwx4LIDdDBG5H/Z85HI8FIF/7acif4YBiQ8qePWojr6fmPFRvg6d47NaRLVZ4rWjGkI+MSOjROIv0QoOPm5ALxOvWhGRe2LAIpfV1SjwxlAVJ6cZENVNYHamhqGbNWwq1KHxF4d2IaXEJwU6Htxgxh8P60h8UEHedAOeDVVh4BQMROTGGLDI5fX2EkgbY0DWRBXtFOCJLyxXWt76TkNlHYNWUx25IBGzVcOMnRr6ewtkTzPgjaEqunDiUCIiBixyHyP9FOyZbMA38SqG+Aq8uF+H/98tk5WeqmTQslXpVYm5WWYM2mRGea3EtsdUpI8zIITTLxARWTFgkdsZ7Kvg76MNKJxpwHOhCtbk6ui13oxpX5jxdSnn0bqdGrNlDcE+n5iRXiTxdrSCY08YMNaf3QgR0c3YM5Lb8vcUWD5YRcksA94epuBEucTwLRqGbNbwcb7OmeF/JqVEWqGO0FQz/v2Qjmf6KsibYcDzHGdFRHRbDFjk9jp5CDwXquL7Jw3YOk6FyQN4KlPD/evMWH5MQ3mNewat4iqJZUc1PLjBjGlfaHiwi8CJaQa8+TAXaCYiuhNDWzeAyFEoQmBCoMCEQMvVrDdPaPjjER2vHNExpodAfE8FEwMFund03XBRWSeRdkrig1wdmWclOhqAqUEC7w5XEHsf/3+MiMhWDFhEt/CQt8DqGAP+Y7DEB3k60oskFuzRICUwxFcgvqfA5J4KHuximeDUmWm6xM4fLXVuPCVRbQZG3Sfw1xgVU+8X8PRw7vqIiNoCAxZRI3w7CLzYX8WL/YHz1RKflUikF+lYelTHSwd19DIB8T0VTO4pEN1dONWYpOxyS6ham6/jx6tAXy/g5QEK/qW3gkBP56mDiMgR2S1g5eXlYc6cObhw4QK8vLzw17/+FWFhYQ32W716NZYvXw5d1zF69Gi888478PDwsFcziFpMtw4Cc0IE5oQoqDFbrvpsLpL4uEDHf54AvI3AhADLla1x/gKdHWxh4zpNIucSsOusjg/yJA5fkPA2ArN6KfhFH4GobsLpr8YRETkKuwWsxMRELFiwAAkJCUhNTUVCQgIOHjxYb5/CwkIsXrwYR44cQffu3REfH49Vq1bhhRdesFcziFpFe4PA+ECB8YHAX6SCwxck0oskNhfp+DBfgyqA4M7AA10EHuwifv6n5XlLT8QppcSZK8C35RLflkuc+PmfOZcAswQ8FEsQfHmAivEBAu1UhioiInuzS8AqKyvDoUOHsH37dgDA1KlTsXDhQuTn56N3797W/VJTUzF58mTce++9AIBnn30Wy5YtY8Aip6YIy9WfqG7AnwapKKyQ+OJHie8vSnx/SeKTH3Scqvrf/e/tgAah68EuAj063f14rso6ie8uSpwoR70wdanO8npnD6C/t8DwexU8H2oZW9bfW8DkYFfXiIhcjV0CVklJCfz8/GAwWA4nhEBgYCCKi4vrBazi4mL07NnT+jwoKAjFxcWNHruqqgoVFRXW50ajEQDw+uuv49VXX7U+d3W1tbVuVzPgnHXfbxKYf9Mix1fNErmXge8vSuRctoSvPed0/E8uUKtZ9vFQAEUAlnlOJcxaByhKHSQsx5KwvHarSSNUAYR4WcLUWH8F/X8OUoGezjMI3xnPNRHR7Tj8IPeYmJh6z5OSkvDcc89hxYoVeO6552AymdqoZa2roqLC7WoGXKvuQAEEegPjvAHcb9mmSaD4ikBupYKSqwqkBIQAampqsGTxYrzypz+hffv2EAK4HpMEYH1uVCUeMOkI6ayjvXrTB14DLl5stfKazZZz7e3t3cqtIiJqGrsErICAAJw9exZmsxkGgwFSShQXFyMwMLDefoGBgSgoKLA+P3XqVIN9bpaVlYWIiAjrc6PRiNraWgCWztbZ/9K11fWrg+5UM+AedXe7Bxh407aKCoklWSn4VfgKmEyd26Rdrc0dzjURuQ+7zBzo6+uLyMhIrF27FgCQlpYGf3//ercHAcvYrPT0dJSWlkJKiXfffRczZ85s9Nienp4wmUzWB28dEBERkaOz2y3ClJQUJCQkYNmyZTCZTFizZg0AYN68eZg8eTImT56M4OBgvPLKK4iOjgYAxMbGIjEx8ZbHq66uBgB8//33DV6rqrKMGD527Bg8PT3tVYJDc8eaAdbtTnXbWvMDDzyAjh07tlaziIiaREgpHXKhtY8++gizZ89u62YQkYM5fPgwIiMj27oZRESNctiAdeHCBWzbtg1BQUHo0KFDWzeHiBwEr2ARkTNw2IBFRERE5KzsMsidiIiIiP4XAxYRERGRnTlUwMrLy8OwYcMQEhKCqKgoZGdn33K/1atXo0+fPujVqxfmz5+Pa9eu2fSao7Kl7p07d2Lw4MEIDQ1FWFgYkpKSoOs6AMt8YqqqIiIiwvq4cb4xR2RLzbt27UKHDh3q1XX916WA657rNWvW1KvZx8cHTzzxBADnPNe/+tWvEBQUBCEEjh07dtv9XO17TURuTjqQUaNGyTVr1kgppdywYYMcNGhQg31++OEH6efnJ8+ePSt1XZeTJk2Sb7/99h1fc2S21H3kyBFZUFAgpZSyurpaRkdHW99TWFgovby8Wqu5dmFLzZmZmTI8PPyW73flc32zsLAwmZqaKqV0znOdlZUlS0pKZM+ePeXRo0dvuY8rfq+JyL05TMA6d+6c7Ny5s7x27ZqUUkpd12X37t1lXl5evf1WrFghExMTrc8zMjJkdHT0HV9zVLbWfbMXXnhBLlmyRErpfH/p2lpzYwHLXc71/v37Zbdu3WRdXZ2U0vnO9Y0aC1iu9r0mInKYW4SNLRh9o8YWjG7KYtJtzda6b1RaWorU1FRMnDjRuu3KlSuIiopCZGQkXn31VWia1uJtb6q7qbmgoACRkZGIiorCO++8Y93uLud69erVePrpp+Hh4WHd5kzn2lau9r0mInL4xZ6pvoqKCkyaNAlJSUkYNGgQAMDPzw9nzpyBr68vysvLMWPGDLzxxhtISkpq49Y2T2RkJE6fPg0vLy+cPn0a48ePh4+PD6ZPn97WTWsVV65cwbp167B//37rNlc910RErsZhrmDduGA0gEYXjC4qKrI+v3HB6MZec1S21g0AlZWViIuLQ3x8PBYtWmTdbjQa4evrC8CyUO7cuXOxe/fu1imgCWyt2WQywcvLCwDg7++PWbNmWety9XMNABs2bEBYWBhCQ0Ot25ztXNvK1b7XREQOE7DssWB0UxaTbmu21l1VVYW4uDjExcXh5ZdfrvdaWVmZ9VdVtbW12LhxIwYMGNA6BTSBrTWfPXvW+kvJyspKbN261VqXK5/r61avXo1nnnmm3jZnO9e2crXvNRGRwwxyl1LKnJwcOXToUNmnTx85cOBA+e2330oppXzmmWfk5s2brfutWrVKBgcHy+DgYDl37lzrAOA7veaobKl76dKl0mAwyPDwcOtj6dKlUkop09LSZFhYmOzfv78MDQ2VCxculDU1NW1Wjy1sqXnlypUyNDTUWteSJUukruvWY7jqub6+n6enp6yoqKj3fmc81wsWLJA9evSQqqpKX19f2atXLyml63+vici9cakcIiIiIjtzmFuERERERK6CAYuIiIjIzhiwiIiIiOzs/wN+oHRKUAqYRwAAAABJRU5ErkJggg==" }, - "execution_count": 26, + "execution_count": 21, "metadata": {}, "output_type": "execute_result" } diff --git a/src/econometrics.jl b/src/econometrics.jl index f0487ef..7da5c51 100644 --- a/src/econometrics.jl +++ b/src/econometrics.jl @@ -44,7 +44,9 @@ function calculate_Avar!(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData # sandwich formula, where tau capture the extra noise introduced by simulation, # compared to the usual GMM estimator - sMMProblem.Avar = (1+tau)*inv_Sigma1*Sigma2*inv_Sigma1 + # I use the function Symmetric, because otherwise a test issymmetric() + # on the asymptotic variance will return false because of numerical approximations + sMMProblem.Avar = Symmetric((1+tau)*inv_Sigma1*Sigma2*inv_Sigma1) end diff --git a/test/runtests.jl b/test/runtests.jl index 0543971..584a13a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -837,3 +837,205 @@ end end + + +@testset "Testing Inference" begin + + tolLinear = 0.05 + + # Inference in the linear model + #------------------------------ + srand(1234) #for replicability reasons + T = 100000 #number of periods + P = 2 #number of dependent variables + beta0 = rand(P) #choose true coefficients by drawing from a uniform distribution on [0,1] + alpha0 = rand(1)[] #intercept + theta0 = 0.0 #coefficient to create serial correlation in the error terms + println("True intercept = $(alpha0)") + println("True coefficient beta0 = $(beta0)") + println("Serial correlation coefficient theta0 = $(theta0)") + + # Simulation of a sample: + # Generation of error terms + #-------------------------- + # row = individual dimension + # column = time dimension + U = zeros(T) + d = Normal() + U[1] = rand(d, 1)[] #first error term + # loop over time periods + for t = 2:T + U[t] = rand(d, 1)[] + theta0*U[t-1] + end + # Let's simulate x_t + #------------------- + x = zeros(T, P) + + d = Uniform(0, 5) + for p = 1:P + x[:,p] = rand(d, T) + end + + # Let's calculate the resulting y_t + y = zeros(T) + + for t=1:T + y[t] = alpha0 + x[t,1]*beta0[1] + x[t,2]*beta0[2] + U[t] + end + + myProblem = SMMProblem(options = SMMOptions(maxFuncEvals=500, saveSteps = 500, globalOptimizer = :dxnes, localOptimizer = :LBFGS, minBox = false, showDistance = false)); + + # Empirical moments + #------------------ + dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}() + dictEmpiricalMoments["mean"] = [mean(y); mean(y)] #informative on the intercept + dictEmpiricalMoments["mean_x1y"] = [mean(x[:,1] .* y); mean(x[:,1] .* y)] #informative on betas + dictEmpiricalMoments["mean_x2y"] = [mean(x[:,2] .* y); mean(x[:,2] .* y)] #informative on betas + dictEmpiricalMoments["mean_x1y^2"] = [mean((x[:,1] .* y).^2); mean((x[:,1] .* y).^2)] #informative on betas + dictEmpiricalMoments["mean_x2y^2"] = [mean((x[:,2] .* y).^2); mean((x[:,2] .* y).^2)] #informative on betas + + set_empirical_moments!(myProblem, dictEmpiricalMoments) + + dictPriors = OrderedDict{String,Array{Float64,1}}() + dictPriors["alpha"] = [0.5, 0.001, 1.0] + dictPriors["beta1"] = [0.5, 0.001, 1.0] + dictPriors["beta2"] = [0.5, 0.001, 1.0] + + set_priors!(myProblem, dictPriors) + + # x[1] corresponds to the intercept + # x[1] corresponds to beta1 + # x[3] corresponds to beta2 + @everywhere function functionLinearModel(x; nbDraws::Int64 = 1000000, burnInPerc::Int64 = 10) + + # Structural Model + #----------------- + srand(1234) #for replicability reasons + T = nbDraws + P = 2 #number of dependent variables + + alpha = x[1] + beta = x[2:end] + theta = 0.0 #coefficient to create serial correlation in the error terms + + + # Creation of error terms + # row = individual dimension + # column = time dimension + U = zeros(T) + d = Normal() + U[1] = rand(d, 1)[] #first error term + # loop over time periods + for t = 2:T + U[t] = rand(d, 1)[] + theta*U[t-1] + end + + simX = zeros(T, P) + d = Uniform(0, 5) + for p = 1:P + simX[:,p] = rand(d, T) + end + + # Let's calculate the resulting y_t + y = zeros(T) + + for t=1:T + y[t] = alpha + simX[t,1]*beta[1] + simX[t,2]*beta[2] + U[t] + end + + # Get rid of the burn-in phase: + #------------------------------ + startT = div(nbDraws, burnInPerc) + + # Moments: + #--------- + output = OrderedDict{String,Float64}() + output["mean"] = mean(y[startT:nbDraws]) + output["mean_x1y"] = mean(simX[startT:nbDraws,1] .* y[startT:nbDraws]) + output["mean_x2y"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws]) + output["mean_x1y^2"] = mean((simX[startT:nbDraws,1] .* y[startT:nbDraws]).^2) + output["mean_x2y^2"] = mean((simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2) + + return output + end + + set_simulate_empirical_moments!(myProblem, functionLinearModel) + + # Construct the objective function using: + #* the function: parameter -> simulated moments + #* emprical moments values + #* emprical moments weights + construct_objective_function!(myProblem) + + # Run the optimization in parallel using n different starting values + # where n is equal to the number of available workers + #-------------------------------------------------------------------- + @time listOptimResults = local_to_global!(myProblem, verbose = true) + + minimizer = smm_local_minimizer(myProblem) + + # The minimizer should not be too far from the true values + #--------------------------------------------------------- + @test minimizer[1] ≈ alpha0[1] atol = tolLinear + + @test minimizer[2] ≈ beta0[1] atol = tolLinear + + @test minimizer[3] ≈ beta0[2] atol = tolLinear + + # Empirical Distance matrix + #-------------------------- + X = zeros(T, 5) + + X[:,1] = y + X[:,2] = (x[:,1] .* y) + X[:,3] = (x[:,2] .* y) + X[:,4] = (x[:,1] .* y).^2 + X[:,5] = (x[:,2] .* y).^2 + + Sigma0 = cov(X) + + set_Sigma0!(myProblem, Sigma0) + + @test myProblem.Sigma0 == Sigma0 + + nbDraws = 1000000 #number of draws in the simulated data + calculate_Avar!(myProblem, minimizer, T, nbDraws) + + # The asymptotic variance should be + # * symmetric + # * positive semi-definite + #----------------------------------- + @test issymmetric(myProblem.Avar) == true + + # test for positive semi-definiteness + # all eigenvalues should be non-negative + eigv = eigvals(myProblem.Avar) + + counterNegativeEigVals = 0 + for i=1:length(eigv) + if eigv[i] < 0 + counterNegativeEigVals += 1 + end + end + + @test counterNegativeEigVals == 0 + + # summary table: + #--------------- + df = summary_table(myProblem, minimizer, T, 0.05) + + # first column : point estimates + @test df[:Estimate] == minimizer + + # 2nd column : std error + for i =1:size(df,1) + @test df[:StdError][i] > 0. + end + + # confidence interval + for i =1:size(df,1) + @test df[:ConfIntervalLower][i] <= df[:ConfIntervalUpper][i] + end + + +end