{ "cells": [ { "cell_type": "markdown", "id": "342fd358", "metadata": {}, "source": [ "# Creating MA, AR, and ARMA Graphs" ] }, { "cell_type": "code", "execution_count": 1, "id": "981b9837", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "d5c4a57e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "By using this Application, You are agreeing to be bound by the terms and conditions of the Halerium End-User License Agreement that can be downloaded here: https://erium.de/halerium-eula.txt\n" ] } ], "source": [ "import halerium.core as hal" ] }, { "cell_type": "markdown", "id": "628435ca", "metadata": {}, "source": [ "## the MA model (moving average)" ] }, { "cell_type": "markdown", "id": "a493264b", "metadata": {}, "source": [ "An MA model is defined as\n", "\n", "$$\n", "y_t = c + \\epsilon_t + \\sum_{j=1}^q b_j \\epsilon_{t-j}\\\\\n", "\\epsilon_t \\hookleftarrow N(\\epsilon_t| 0, \\sigma^2)\n", "$$\n", "\n", "we can create this explicitly in Halerium using the ``TimeIndex``." ] }, { "cell_type": "code", "execution_count": 3, "id": "0bb3ef3c", "metadata": {}, "outputs": [], "source": [ "q = 5\n", "sigma = 1.\n", "sigma_obs = 1.e-2\n", "\n", "ma_graph = hal.Graph(\"graph\")\n", "with ma_graph:\n", " hal.StaticVariable(\"b\", shape=(q,), mean=0, variance=1)\n", " hal.StaticVariable(\"c\", shape=(), mean=0, variance=1)\n", " \n", " hal.Variable(\"epsilon\", mean=0, variance=sigma**2)\n", " hal.Variable(\"y\")\n", " y.mean = sum([\n", " b[i] * epsilon[hal.TimeIndex-i-1]\n", " for i in range(q)\n", " ]) + c\n", " y.variance = (sigma/100)**2\n", " # TODO: or the next\n", " hal.Variable(\"observed\", mean=y, variance=sigma_obs**2)\n", " # this \"observed\" variable is strictly speaking not part of the model,\n", " # but it is recommended when we want to fit data to the model as it allows us\n", " # to include observation noise." ] }, { "cell_type": "markdown", "id": "fabbace8", "metadata": {}, "source": [ "The ``[TimeIndex...]`` is a shorthand convenience for the ``TimeShift`` operator with initial values set to 0.\n", "\n", "So instead of ``epsilon[TimeIndex-i-1]`` we could write ``TimeShift(epsilon, shift=-i-1, initial_values=0)``." ] }, { "cell_type": "markdown", "id": "d145cdbf", "metadata": {}, "source": [ "Note that in the Graph definition we did no specify the length of the time series. Halerium will utilize the dynamic axis of the ``Variable`` class for the time steps, which means we specify the length of the time series with n_data." ] }, { "cell_type": "code", "execution_count": 4, "id": "df3b5e68", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABSnklEQVR4nO29eZQk113n+70ZW+6VtVd1V69Sd0ut1t6W5AXbWoxk4wXGw9gY8wzM2IwHxgZjA8YDDO89HsOBMea9BzY6ZvV4AcvCCw+wJNvyLlndklpbqzepl+qqrj33JSIy7/sj4kZGZkZmRmZGrnU/5/TprsrqzFuZEd/4xff+FkIpBYfD4XCGF1+/F8DhcDiczuBCzuFwOEMOF3IOh8MZcriQczgczpDDhZzD4XCGHLEfLzo1NUX37t3bj5fmcDicoeX48ePrlNLp6u/3Rcj37t2LY8eO9eOlORwOZ2ghhFxw+j63VjgcDmfI4ULO4XA4Qw4Xcg6HwxlyuJBzOBzOkMOFnMPhcIYcLuQcDocz5HAh53A4nCGHC/k25vGXNnDqSqrfy+BwOB3ChXwb85tfegZ/9o3T/V4Gh8PpEC7k25RSiWIpnkemUOz3UjgcTod4IuSEkBgh5AFCyIuEkJOEkFd68byc7rGZVaEWS8hrXMg5nGHHq14rfwbg3yil/54QIgMIevS8nC6xHM8DABdyDmcE6FjICSFRAK8F8PMAQClVAaidPi+nuywncgCAvFbq80o4HE6neGGt7AewBuBvCCFPEUI+TQgJefC8nC6ynDAi8hyPyDmcoccLIRcB3ALgk5TSmwFkAPxW9Q8RQt5HCDlGCDm2trbmwctyOmHJisi5kHM4w44XQr4IYJFS+rj59QMwhL0CSun9lNKjlNKj09M1fdE5PeYKj8g5nJGhYyGnlF4BcIkQcsj81t0AXuj0eTndhW12FrhHzuEMPV5lrfxXAJ81M1ZeAvALHj0vp0ssJw1rRS2WUCxRCD7S5xVxOJx28UTIKaVPAzjqxXNxuk+pRHElkYckEGhFirxWREjpy9Q/DofjAbyycxuykVGhFSn2TBrJRXzDk8MZbriQb0NYDvleU8j5hieHM9xwId+GsBzy/dMsIucbnhzOMMOFfBuyHDci8n1T3FrhcEYBLuTbkOVkHrLgw45YAAAXcg5n2OFCvg1ZjucxN+ZHUBYAcI+cwxl2uJBvQ64kDCH3i4aQc4+cwxluuJBvQ64k85iL+hGQjY+fR+QcznDDhXwbksxrGAtIUKyInAs5hzPMcCHfhmQKOsJ+EQHTIy9wIedwhhou5NuMgl6EVqQIKyL8Et/s5HBGAS7k24x0XgcAQ8hF4+Pnm50cznDDhXybkSkY0XdIESEKPkgC4RE5hzPkcCHfZqQKGgAgrBi2il8U+GYnhzPkcCHfZrCIPKxIAAC/zIWcM9q8tJZGMq/1exldhQv5NiNTMDzyEIvIJR/3yDkjzX/4y8fwqUfP9XsZXYUL+TYjVShvdgJAQBKQU3lEzhldtrIqVlOFfi+jq3Ah32awiDzsN4TcLwnI61zIOaOJbo4yTHFrhTNKlK0Vm5Bzj5wzouR1wzZMmWm3owoX8m0GO6BDclnIc9wj54worGqZCzlnpMgUdARlAYKPAAACko+X6HNGloIVkXNrhTNCZFTdslUAFpFzIeeMJnkekXNGkVReR8Qu5LwgiDPCFLhHzhlFMoXKiDwg8/RDzujCghS1WBrpgIUL+TYjXdCtYiAAUCSftbPP4YwaBduxPcpRORfybUa6ULTK8wGjIEjVjVxbDmfUqBTy0d3w5EK+zcgUdKthFgCrJ3mBFwVxRhC7nVIdkX/+Rxfx9eev9HpJXYEL+TYjXe2RS3wAM2d0aWSt/NX3XsYXfnSx10vqClzItxlpc8wbwy/xAcyc0aVQEZFXWis5tYjMiGz0eybkhBCBEPIUIeSfvXpOjreoegmqXkJYrswjB/gAZs5okm8Qkee1IrLqaGyAehmRfxDASQ+fj+Mx1Q2zgLKQ8xREzihij8ire5LntCKyhdE47j0RckLIAoCfAPBpL56P0x3SVQ2zgLJHzjc7OaNIPY+cUoqcVrTOiWHHq4j8EwB+A0DdHTNCyPsIIccIIcfW1tY8ellOK6SrepED9oicb3ZyRg8WkYdkoULIC3oJlALZEbkT7VjICSFvBrBKKT3e6OcopfdTSo9SSo9OT093+rKcNsg4CrlxCHCPnDOKFPQSFNGHaECq2OwsmFlaGVUHpcNfQ+FFRP5qAG8lhJwH8AUAdxFC/pcHz8vxmEbWCs9a4Ywiea0IvyQg4hcrInJ2vFM6Gqm3HQs5pfSjlNIFSuleAO8E8E1K6bs7XhnHcxpZKzwi54wiLCKP+KWKzU574JIZgcwVnke+jWiUtcL7rXBGkYJegiL5aiNymzeeGYENT0+FnFL6KKX0zV4+J8c70maqVWUeuemRj8imD4djJ68V4RcFRPyVHnlFRD4CKYg8It9GpNmYN4deK9xa4Ywi9SJy+/E+CkVBXMi3ERlVh1/yQRTKH7sk+CD6CN/s5Iwk5Yi8vpCPQpk+F/JtRCqvV2x0MgKSMBI79xxONSwij/qliuES9sAlyz1yzjBhtLCtFXKFz+3kjCgFvQjFjMiBcnVnxWYnj8g5w0T1mDeGX/JV9KTgcEaFvFaC3/TIgXIHxAprhUfknGEiVUfIAzwi54woVkRuTsWyInKeR84ZVjIFHREnIZeFoe058dEHn8XHHzrV72VwBpS8xgqCqq0VY0/IRzASHRC5kG8j6lkrMxEFK8l8H1bUOY+cXMGDT13u9zI4A0rBKtFnEblprehFyIIPIUXkETlnuEjlnYV8RyyApXiuDyvqjLxWxFqqgMWtHFZTw3kh4nSXcol+7WanIvkQkkUekXOGB1UvYTOrYiai1Dy2IxZAMq8PXW/my7aLz5MX4v1bCGcgoZSa6YcCogEjIk/aNjsDkoCgIiDNI3LOsLCayoNSYH7MX/PYjlgAALA8ZFH54pZNyC9u9XElnEGEDZVQRJ+Vdmvf7AzIAsKKyPPIOcPDlYRhPcybom1nhynul4dMyC+bQr5jzI8nL3Ah51RiF3LBRxBWxAprJSAJCMoCzyMfFooliodfWBmJBvLtssyEvEFEvhQfLp95cSsL0Udw75E5PHM5AZV3cOTYYLURrJ9QxC9a1krO3AQNySLvtTIsfPfMGt7798fw9KV4v5fSN5YTRvQ65yDkMxEFgo8M3Ybn4lYOO2IB3LZ3AqpewvNLiX4viTNA2CNyABgLSEjkDCEvaCXTI+ebnUPDWqpQ8fd2ZDmRR0gWHPPIRcGHuagfS4lhE/IsFsYDuGXPOADgyYvx/i6IM1DkqyLysYCERLYckQdkASFZGLpNfie2hZDHzQ+P/b0duZLIYz4WACHE8fH5Mf9QRuQ7YwHMRv3YGQtwn5xTQXVEHguWI/KcmbUSUsShLYazsz2EPKcCALayap9X0j+WE3lHf5xh5JIPj0de0ItYTRWwMB4EAFy/cwwvLCf7vCrOIFHQDYFWbBE504JyHrkwEgOYt4WQb5mR+NY2jsiXEznMRRsL+ZVEHqXScBzQ7KKzMG5s1B6ai+D8Rqaiqx1ne8NaM/utiFy2IvJyHrk4EgOYt4WQx81IfCuzPSNyrVjCaqrgmHrI2BnzQy2WsJ4Zjn2Exa0sgLKQXzsfAaXA6ZVUP5fFGSCcIvK8ZvQkt6wV2Xhs2Mv0t4mQs4i8e0L+zRdX8J3Ta117/k5YSxXqFgMx5seGKwWRFQMtTBjWyqG5KADg1BUu5ByDglabtQIAiZxmROSygKA5v3bYW9luCyHf6sFm5yceOYO/ePRs156/E1gOuVPqIaOcSz4cG54sh3zWbDmweyIIv+TDi1zIOSZ5vTZrBTACmxI1vs/m1w77AObaXLQRxLJWuhiRJ3MaigPqL7Mc8h1jjayVYRPyHOZjfmv+qOAjODQbwYtX+IbnsFEqUfh8ztlUnVAdkceChpCzKmeWtQIM/wDmbRGRx3uw2Zka4KZTV1xE5NGAiJAsDI21cnkrh4VYsOJ7h+Yi3FoZMr7y9GXc/offwGYX9q+c8sgBYNls2VxhrQz5JvnICznb2BB9BPGs2rU0o1RBRzo/mEK+nMgjKAuI+uvfgBFCMD9E7WyvJPM1F6ZDc1FsZNRtXfg1bDx1MY61VAF//8Pznj93TR55QAYArJiBjV/yWdbKsDfOGnkhZ9H47okg9BJFqlA7fLVTCnoRql6ynnvQuJIwRK9eMRBjLurHypD09Y5nNYwH5YrvXTMXAQBurwwRlzaN7KO//cF5z+0NpxJ9wAgCANNa4RH5cMB88X1TIQBAPKPh/HoGN/z+13Hco0pA1lFN1UtWytMgsZTINfTHGbFguYR5kNGKJaQLuuV5MpiQc3tleLi4mcXOWADxrIZ/eOKSp8+dN+/E2T5KxC+CEFjTsPxm90OAZ60MPCwiZ0K+lVXx/FISWpF61mQpZbNUBnH3m0XkzYgFJcRzgy/k7DMdrxLyybCCqbDCM1eGBEopFrdyuO/IHI7uGcenv/uypwkDbDoQw+cjiPolx81Onkc+4LCMlX3ThpBvZlVc2MwAqBxM0Al2b3zQfPJSiWI1VWhY1cmIBWTEs+rAV3cmzDLrsSprBQAOzITx0lq610vitMF6WkVOK2L3RBA/dctOXI7nrAwrLyjoRWujkxEL2oRcFqxe5cPeAXHkhXyrKiKPZ1Vc3DB8OVYd2ClsoCtQHiU1KGxkVBRLFDPR2hFv1cSCEkoUAz/6ikXksYBU89hYQBrIuyJOLZfM82/XRABTYeP49LLWI69VRuSAcXywvayAJIAQYg6XGOxjvhkdCzkhZBch5FuEkJOEkOcJIR/0YmFewZrkWNZKRsMFU8gvexSRJ+0R+YB5bcwPnIm4sVaMCHfQffIty1qpjciDyvCflNsFttG5azyIiZDxWXpZ68HmddoZs138WbQ+CgOYvYjIdQC/Tim9FsAdAH6ZEHLYg+f1hHhWg1/yYSbih4+YEfkmi8i9EXJ7RD5o1gpLxXMVkZsH+aB3iWR2WfVmJwAEZWEk2pJuB5iQL4wHrf0OL/PJ81rRMSJnMCEfhQHMHQs5pXSZUvqk+e8UgJMAdnb6vF6xlVERC8gQfARjAQkryQKWEjkEZQEbGdWTlCd7FD5oEflqikXk7qwVYPD7trMOdmMOQh6SxYHPQPjpT/0Ad/3PR/H/PbM88PsR3eTSZg5TYQUBWbDurrw89pwicvvFPyDbI/LBPmaa4alHTgjZC+BmAI87PPY+QsgxQsixtbXeNZeK5zTrwxsPynjmcgKUArftmwDgjb1iz1rpdi755XgOR37v665T7FaTRkQ+3YqQt5i5Qim1quh6wVZWheAjjtOOgrKIgl4a2HYJBb2I4xe2sLiVwy9/7kn8yUOn+r2kvnFxM4vdE0ZaLIuUvYzIC80icvOxkDL8A5g9E3JCSBjAlwD8KqW0piKDUno/pfQopfTo9PS0Vy8LwPCBv3jMOQc1nlWtq30sKFltTl991RQAb+yVVF6DYPaK6La1cnEji3RBd92udTVVQCwoQRGFpj87FmAeeWsn01dPLOGOP/xGz3qBx7MaYgHJscDJqtQb0Fvl8+tZlCjwR2+/HjcujOHJi9t3qtGlrSx2md0rRcFnDH7w0NbL66WarBUm5LLgs/LLR2EAsydCTgiRYIj4ZymlD3rxnK3wwPFFfOSBZywbwc5WthyRT4RkK1J71dWTALzJXEnldUyEZIg+gnShu7ZETjMOOLdR80oy78pWAcoReas9aV5ezyCe1TxNHWtE3PaZVsN6Z9TzydMFHRc2Ml1bWzPOrhqpkQdmIlgYD458O4Hz6xk8cHwRL15JVtwl6cUSlhN57Bov98sZD0rY9NJacYjIWZm+Xyp/fxQGMHuRtUIA/BWAk5TSj3e+pNZhV3GWVlj5mGZlY7C/g7KAa+eikAWfRxG5johfRMQvdj0iz6lG2bHbqHk1VXCVsQIAkuBDWBFb9imZJ81Kn7tNPKdan2U15bakzp/DX377HN78f3+vp1aQnbOraRACXDUdxnREGXkh//jDp/HhL57AfZ/4Lv7dX3zf2hNYTuRRLFHsmihXHI+HZE+Hv6gOEXnUjMiZPw5gJAYwexGRvxrAzwG4ixDytPnnTR48r2vY5hdLKzx+YQv3/ul3sJrMm9YK88iNv3dPBOHzEewcD2DRgyZRqYKOiCIi7Bcr/PJuwG4B3YrtWqrgKmOFYZ9r6BbmL670SshNa8WJgMSsFWehXk7kkSroeOL8ZtfW14hza2nsjAUQkAVMhWUk83rfLiq9YHErixsXxvCBu67GicUEvn9uHQCszDFmrQDARFD2NGPKKWuF3ckFbAIflId/ALMXWSvfo5QSSukNlNKbzD//4sXi3MJEjd0yf+f0Gk6tpPDn3zoLvUStD49FcXsmjYNnYTzgmUce8UsIK1LXNzvZSe/GWqGUGkLuMiIH2uu3YkXkid5El/a7rGqskus6n0PSfN++fao/05zOrqZx9UwYQHkDej09ulH5UjyPA7MR/Jc7r0YsKOELZj8Vew45Ixb0NiKvLtEHyh65PVIPK8M/gHkkKjutiNw8OJgP+bkfXQRQFnBWdLBn0igOWhgP4LJHHnnELyKi9MBaYULuQmzjWQ1qseTaIwcMIW81KmKi2buIXG3gkTeOyNmx8p0zvRfyUonipfU0rpquFvLBzttvF61Ywkoqjx2xAPySgJ+6eSceev6KkZxwfBFhRawYPzgRkqz9mbVUAf/9q89D1dsfipzXnEv0gUohH4UBzKMl5Ka1cnolhYXxALSicYVlWSt2awUwChHW0youbWbx6/94Aufa7NFhROSGtdJtr40JVMKF/cFa0rZircSCcsvph+keCrmql5BRizUNsxjliS/OQs6qcE+vpHvee/1yPIe8VipH5GFDxEbVJ19J5kGpMdgbAN75it3QihT//lM/wPELW/jDf3e9lTkCGMdeTisirxXx8Asr+NsfnO9omHajiNxurdQbwPzJR8/hv3352bZfv5eMhJCz2+WLm1loxRJeXs/gLTfuwG17jVxxdtIvmLdx185Hza+NjZa3f/IH+NKTi/j+2fW2Xj+d1xFWJISV7gt5KxE5yyFvyVoJtGOtGGvqxWZnvEHDLKAckdcr00/mNFy/cwwAej4sm90pMiGfihi/w6gKOZs2xQZ7H5qL4ObdMVzazOE/v+4qvOXGHRU/by/TZzZpu961XixBL9GaiDwgCZAFX8Vmp9MA5hevJPEnD53CN0+utvX6vWYkhDyRM/K4NzMqnr2cgF6iODgbxgfvOYDJkGxZKUd2juE7H7kTt+4ZB1AWcrZB2c5GZbFEkVGLVkTe7c1OlqvtJmpeZeX5LVor8ZzWkl/ITgB24egmiQYNswBYgwLqVeolcxpu3TOOuai/5/aKJeSmtTIZMj6X0RVyc1ZsrJyZ8jtvPoxfufNqfOTeQzU/by/TZ3fX7fbNyZuWjD3NEDAmYUUDUmVEXjWAmVKK3/3y8yiWKLJDshE99EKuFY1b7UOzxlCBb5xcAWDk6b766ikc+2/3VFQ17p4sb64cnh/Dm66fw2f+422QBV9bnQuZJ2555N3OI2fWSra52K62Y60EjFz7Vu4s2Mm2ksx3veS8UcMsoJxW5lSpVzQnRI0FJLzu4DS+e2a9pxtcZ1fTmAzJGDcjT1n0YTwoYS09HFOZWuWyJeTlO8Jbdo/jw/cesgro7LDPdCujWftd7eZ3s/OkOiIHgLfeuAOvPThlfV2uPTCO43966jJ+dH4TO2OBockvH3ohZ/74DQvG7fIjL6xaeboAGo43C8gC/uJnb8XRvROItBlNM/GP+g1rJa+VoBW7t2nCrBW1WLL+nS4477ivJgsIK6J1oLphrI1+K5lCEbLog16i2Ggz6yCZb3xh+vR3X8Kx85sNG2YBKPeXdojkWHOzsYCEIwtjSOV1666lF5xdS+Mq01ZhTEcUrKdGc7NzKZ7DeFByffyxC9ymzVppOyLX6gv5777lMN7xit3W1+XhEsb/+dqJJeybCuGnjy5ALZagd/F89oqREfIjpu95aiWF3RPBCg/MDe0KOYtcmbUCdLdM314GH89qWEsVcOv/8TC+7eD3rqbyLUXjQNmycCvklFJkVB37TPuqnQ3PRFbDbX/wCB56YaXuz3z84dP482+dtSylekLO+ks7eavJnPG5RAMS9pp3ZufXe1PlSSnFmZWU5Y8zpsIK1kY0/XApnquwVZrBIvKzKynr82u3mRUbuegk5NVUD2BO5DTsjAUQZhvnQ2CvjIyQ7xwPYNK8oh+oOlncEPFLFe1o3cLEP+wXrQ++mxueOa1SyM+tpVHQS4758KvJQkv+OFCOitwWBWXVIigF9k+3L+Rr6TzyWsnykKthTbkef3nT8pPr5ZED9ftLs2Ml6hex17zwnO9Ruf7FzSySed3aaGX0srpzKZ7DvX/6nZ5NUFpO5K2NTjewi/NTl+LW99ptZsVSCf1ic4mrHsCcyGkYC0hWMNirHkKdMDJCPhaQrEKfA6Zf3grtRuRM/CN+CRG/ZH6ve0KeVYuQzZSteE61NpScqgNbKc9ntBqRs41OZmW1k7mSNkW3Xuc7tVhCiRq/+zdfXIXoI1bKmBP1hkskbdbK/JgfkkBw3qGtQzd4ZtGYD1sj5GFDyHvh1X//7DpOraTwb89f6fprAYZHvjPm/viTBB8ifhEnbELutplVqUTx+1973goGWMDj5s68egBzIqcjGhCb1iQMEkMv5MkKITeirIOz7UTkYkcROeu1AnQ3Is9rRWuQciKrWW14qw82SqlhrbQYkZc9cncROYti9kwGzQnlrUeX7ASqJ+T2Qo3jF7YQC8oN9z5CdUqurYg8IEEUfNg1HuxZA61nLycgCz4crAoypiMKclqxK21Uz6ykcM/Hv41V8+L6/JLRlPSxl7rfniCZ15DK6y1ZK4BhryTzOnzEEFi3Y/tWUnn8zffP49FTRrpgI4+8muoBzMm8Zma2VG6CDjJDL+QschwLSFahz4GZdiJyqcOI3G6tdC9zJacVrWq4eE7DktlxMFcVkacLOvJayVUfcjtjbUbkYwEJU2EFK4nWI3L2vtfbKC1U/W71/HFGQBYcS/TtF33AuPicX+9NRP7sYgLXzkcgV93qs1mV612wV569nMDZ1TS+c8aoj3jBFPJj5ze7uiEPAMtmDnnLQm5ae/NjAYwFJNciyo5XdgEvWyvNhdw+gDmvFaHqJUT9khWRc2ulB9itlTccnsV9183VRD1uaNtaYZudimRtdnbbWrGEPKtZ3nj1wcbKvplQuEURBQRlwXV1J7v7CCsi5qL+tqyVckTuLGbspGR5//VyyBmhOpud9ogcMFo1XNjIdN3WKJUonrucwPULYzWPsQttNzY82e977PwmSiWKF5aTmI0qyKpFPHs54fnr2VlySD10A8sl3zsVNIciuxNRJuQsoMlZEXlzibMPYLZf7Lm10kMSOQ0hWYAk+HBk5xg+9XO31kQ9boj4JaQLesuTZVJ5HaKPwC/5rIk1XbVW1CLGQzJk0dfQI2cbaFMtRuSAcXvbakQeUkTMRv1tbXay92uzTs+RvJmBcM+1swAab3QCZn/pOh65YPPX904GkVGLXe91cn4jg1RBxw07YzWPWULehYicCfkT5zdxcdMYSPK/vXIvAODxLtsr7E6x1Yh8wvxsd0+EEFLcj2BjLStyVkTu3loByhvk9sAwwIW8d7Ad5k6Jtulvsz4rhJDepB9qRQRlAbGAhHhGs4ouqq0V1lFvusWIHDAOYje9XIDy+xVSBMyNKR0J+UZGdYyO2cl5x/4JRBQR05HGQl4vIk/mdETNzwoA9kwZeyrd9slZ9NswIu+ikJ9by+C7ZvuJ1x6YxoGZMB57acPz17OzFM9B8JHWN9tNId872V5Ezi7ghRaFnA1gtupCAuX8dzbMZZAZCSGPeiDkEcsWac3fNjoflhvx+Ej3rBVVN/pHBCQBsaCElzcylu1QLVxMyKeaiJ4TRgdEtxG58bohRcRsxI+trNa0v3ZeK+I3H3jGEi8W1Rf0kqMAs+eL+CV89r234wN3H2j4/ME6A5irL/r7rBTE7vrkzywmoIg+x7TY8aAMH+mukAPAZ354HqKP4MBsGHfsn+y6T355K4e5qN+xgrMRE6Hy/kUrI9jY75ozzwfLI3dhrQDlAcwJB2vF7YZrPxl+Ic96E5G3mzpoNMwyLgKEkK42zsrZooxYQMbJpfJo1GrxXE8VQEj5VrUVQor7SfTsRAspIubN2+hmXQWfX0rgH45dwg/NqND+fjllrtj7ZtywEGuamxxSBDO/vTK6Z9kIjJ3jAQg+0vWioGcvJ3B4R7Si0x9D8BFMhpWu9CRP5jRcNR2CLPpwesXog+6XBNyxfxIZtYjnuuiTv7yesdKBW4Ftdu6eCLU0go3t6eTM4zHXqrViDmC21xrwPPIeksjVn9/YCpE2NypZL/Ly87SX/eIGdkAFZRFjwfIQi6hfrDnY1tIqJoKyo3g0wy8JKLjsA21ZK7JoVUuyPhn1YLfB7KRpKuTmSelmgDRgvD96iUKtijirI3JJ8GFhPND1oqCTS0lctyNa9/HpcHuWVDMSOQ3TEQU3mpbOYXMN7O+Xu3QBo5Ti3FqmrcK8Nx2Zx8fedC2umYu0NIKtNmulCNFHILk8/ln0z6p/xwISgk2mTQ0SIyHk3kbkLVorhUohD3excVa5yMFXkblx9UzY0SNvNWOFoYg+1+PHMgUdAUmA4CNWQzKn2al2mICzDAH7nkIjIXftd9aJpJI5DVF/5bFiZK50z1rJa0WkCnrDu4h90yGcqVPVChjTdB58chEXN7ItZdiwc+Oo2c75uh2GoHfTlweMorB0Qa9pR+CG8ZCM9752P3w+0tIINmuzUyunH7o9XoDyAObqWgNZ9CHbgUfeTiO+duBCbtJuRJ7XigjI9oi8i9aKeVAHJLFi9uD8WMBZyNvwxwHDwnAfkRetXhXTYQVBWWga4bKThRUdZVTduug45ZLnW6jSA2pLrsuvq9fsp+ydDOJ8F1MQ2e/Dem07cXg+isWtXN2T/v/55hl86B9P4LV//C38/N884fq12bnxmquNTn+37I4BMDaDA5LQNSFn1ZXVDcJaJdTCCDYr/ZBF5HrRtT8OlAcwJ3MagmYWHGAEBe1aK8cvbOLm//1ha6xdNxlqIS/oReS0okdZK+1F5Dm1iIDtgIkGpJan0Dfi3Foad//PR7GWKli75wFZsHb3d8T8CDgcbJ1E5H5RaCkiD9n2CNxEuDXWSl63/FSnXPJW+mYARgYCUNtwyfDIKzvxXT0TRiqvezK71YktF0J+7bxR9/DisvM0nPW0in1TIbzzFbvw7dNrroWBCfmrr57Cox9+PW7ebfThJ4QYPV661KyreoBGuwRl9yPYqvPI82rRtRXHXitreuT2u7ag5JwB5YYzK2kUSxTLbRTJtcpQC7l9h7lTWESebDEiz2nFiib1kyEZGx7mJb+wlMS5tQzOrKSQU40DOiAJ1u+8c9xI07JH5JRSrKfU9oVcMoTcTSSUKehWBAwAeyaal72zz83ukc9GFciCr2FE3kpOMFAZkbOKvepj5fZ9kwBgbbx6Dft9JhsKueFZv3gl6fj4VlbFzlgAP//qvQCAx19ungNe0IvIa+Xfd6+ZasmYCstdG/p8djWNqF9sK/XVjjXwwUXmCjuWsraIvJUOqCz6r77DdwqS3MLe316U+A+1kFtVWG1kZlTjN0dAtepp5bUi/LYDZiqiYCPjXROkjC3Hmh0QQVmwrJWdMT8CUuXBllGNO5VWy/MZiuhDicKaedpwfWo5awcA9kwFcWkz17CwqlrIM4UiwoqIiZDsWBRUTiVrzSO3R+RJKxuhUsgPzoYxGZLx2LnuCDm7w2gUkc9F/YgFJZxcdhbyRFbDWFDCwZkIYkEJj7u46Ng37ZzoZtfFM6tGhkyjfjhusAY+uMhcYTZduSCo1Jq1Yg5gXkkVKu7agi2kQFbDCs16sVk61ELuZUQOtF6mXypRFPRSRT+HyZAMrUhbjuzrwaLKraxak34IADvGAlaWCZvOw/p2dBKRA+Wezg3XZ/PIAWDPRAhqsYTlRH2roizkxnuUNu2ZiZDsuNmZ04qQBOI6JznoEJHbOx/aIYTgjv2T+OFLG13xydndWSMhJ4TgmrkIXqhjrWxlVYwHJfh8BLftnXAVkVe3I6imm0J+bjXdVr+jauoNRa6GTQnzEeNYYW2P3fRZqX6t5XiuJiJvV4iZdeU2lbcThlLIKaXQi6W+CzkrHbffwllNkDy6bbUi8rRasek3bhZOLEwErNdn67GKgcLtb3YC7rzJTEFH0BaRWymIDXxye9ZKqWQMpogoIibDcl1rpbUMBJY2Vv4sGwnbHVdNYjmR70r2ylZWheAjNXcC1Vw7H8WpK8maO5lSiRoptuaF+/b9k7i4mW2aq9/s3JgOG8VbXhcFbWVUbGTUjv1xANZx1SwiZr/rTMSPopl2mtNas1bYxX8tXag4Rqpty1bYsKwVHpE78slvn8Mr/uARfPuUMRXHOyFvbbhEOYvEFpGb4umVT25v8coOiKAk4PB8FH/09uvxxiPzNc19ykLerrViXhhcHMDpgo6w3SO3yt7riyK7DU7kNGQ1YzBFo4i8oLcm5NYAZntE3sBqeOX+7vnkmxnVqN5scjdx7XwUea1Uk/GTyuso0XLHx9v3GamEj7/ceK3VnR6rYRlNXu7nAMY4O6DzjU7AFpE3sVbYRidr75xXS8hrpZY2O9ldJaWV9lu9aVNu4NZKAyil+OKxRWxlNfzdDy8AaN4Nzy2tR+S1ZcBWGp1XEbkZjWzarJWALIAQgne8Yjf8kmCJHLuwrJkHUNseufn7uLNWylkrgOH3yoIPFzYyOL+ewX/9/FM16Zh2S4VFU2G/iPGgbGV52GnV7wxaU9EdInK/WPPzV02HMB1R8MMu+OQbabXhRifjMNvwrLJX2KQmNgbt2vkoon6xadOrptZKuDu55F5lrACVQ5HTBR2/8+XnHAMtlkPOOi1mNR0FrbX0Q/tc0QprRXLfuKsavtnZgDOraby8nsEH7j6Ag7NhyOZUES9odbiE06RuFpF7Z62Y03PSKvJqEYQYm5F22B0Bi6CZR97Il22E33q+xrfdhi1SRNjmkQs+gl0TRrXkH/zLSXztxBKOX9iyHqeUIpFTrbsIZhGEFRGTIRmpgl5zAWnV73SqyLM3Q6qGEIJXdskn38yolg3WiKtnwhB8pGbDk/W8YRG54CO4bd9E07uHptaK1T7X29S4s6tp+CUfdrbY9dAJK2ulUMRj5zbwmccu4EcO+wNWRB41XjNrbva3dBdnC0bsx0hIEdqa2akVS9a6etGrZeiE/OvPXQEhwLtv340H3v8qfOn9r2qrDN2JVsvrLc/adsCw3iZetUattlYCklCTDcBen0Xs6+kCxoOS6/LkatxudrLXs58EALB3MoQfnNvAw+Yw5TMrqYr/oxWpNQTELuQT5kVwK1N5MW3V7xQFHxTRV7FJlsg6Z60wbt4dw1qq4HlL282MislQ8zsjvyTgqumQg5Ab67G37n3doRlc2MjimcV43edzK+TrKW9/35fXM9g7GWpqJbnBHpEvmy0MnKwgJpgsIs+pxoCIQEtCXv5ZLzY77RZhL7oneqKAhJD7CCGnCCFnCSG/5cVz1uPrL1zBzbtimIn6EfVLjq1B26Vla8Wh4lAUfBgPStioMyShVaqtFaeDs7q5TyfFQEA54m8Wkdt7kdvZPRlEKq9jKiwjFpRw2ibk7KRjQs4KcUJmRA6g5r1rNSJnz2dPW0vmNQQkoW6vehbxtjPurxGbWdX1ndGhuShOr1ZaK4mqiBwA3nbTDgQkAZ97/GLd50pUVShWw44Pr4uC4tn26xeqKeeRF3HFzIJadzivWMMs5pHntGLr6Yc2a8VuvwUlEapeanlOgd2yGoqInBAiAPhzAG8EcBjAzxBCDnf6vE4sbmXx3OUk7r1urhtP3/JwiXod1qbCimebSGww8RaLyB0iU6sBvhWRq23744DdWml8ANp7kdthE+p/5c6rcc1cBKdXyn1EWKToGJGbkWv1hmdeK1m+vVsCtoo8vVjCo6fWGnbjiyjeD87Wzdtrt0J+cCaMS5u5Ck+VReTjtog86pfwlhvn8dUTS3UvPM1aV/glARG/6LlHbrSV9sbq9IsCCDHqAdjoOKfzKpHTQAgwGzU9crVolui3krXiHJGXEwlaOy7s1uqweOS3AThLKX2JUqoC+AKAt3nwvDU89Lxxq94tIW91uIRT1gpg+OReeeRso0UvUawk884RORNettmZ6iwi91ubnc0ictNakStP3DddP48PveEg3nX7HhycjeDsatrynpmQM1GtFHJDrGqFvLWTEmCtbI337oHjizizmsav3lO/j3m4C4Ozmb896TIN9IA5NPzcajlzxT6T1s67bt+DrFrEl59ecnwuNz2IulGmn8zrTVMt3eLzEQQlo70sK3N3SiJIZFVE/ZJ1HCZyGih1X0AG1PfI221lu2FLOBiWrJWdAC7Zvl40v1cBIeR9hJBjhJBja2trbb1QKq/j5t2xmnJjr2h1uES9iHzSw4g8U9Ahm7fHS/Gcc0Tu4JF3JOQu0w/t8zrtTEcUfODuA5BFHw7MRpAu6NaJaFkrZtTOJhxVWCtV711Bb62THWAOl1CLyKo6Pv7wadyyO9YwAGi3aVoj2AVp3GXl8QFz1mylFaUi6hdriqFuXBjD4floXXvFzcCVqbD3RUFJj5rYMdjYPjYL1qnOIG62smbRM8t8auWYUUQf2FvsHJG3JsYskNszEXQ95agTvBByp12NGm+CUno/pfQopfTo9PR0Wy/0wXsO4MH3v6qt/+uGVodLFOpMIZkKeReRZ9Qido6zgQ11InK5LOQ5tYisWmy78yFQTj9s5pHbh0rU46CZhsbEKVllrVw2PfKIX8RYQIIkkJr3rroxmRtCioCLGxl84PNPYTVVwMd+4tqGJeNla8U7j3zTRZ8VO3smgpAFX0VL23hOs4Yt2CGE4M03zuPkcrJiEhDDjaBORxQrw8kL8loRBb3kycQuhtGVsGhVCjttRsezGmK2GZublpC7P2YIIdZxXF0QBLQn5IroMyLyIansXASwy/b1AgDn+z0P6LR/QyNajcpyDlkrgBHpJPO1aXStQilFpqBb0+PVYqmhR55Ti5YIuhUPJ1hE3mz96TqbnXYOmlHmGdMnZ3nRU2EZAfO2WfARIyLyEccosVW/EzB85PMbWTz20iY+cPcB3LpnouHPR7pgrTBBmXBprYiCD/unQxVZPltZre6w6V3jxsXwikN3PVfWiscROTtvnHL12yUoi1iK56ygwslaiec0jAVl6zxg+wotb5DLYsVwbgBWi+pWM0/W08ambys91TvBi3f8CQAHCCH7AFwG8E4A7/LgeXtOq8Mlcg5ZK4BhrQDGidxsLFkj1KIxo3PXRHmTrpFHnlOL1kE84SLlrR5u88jL8zrrnzDjIRlTYdmKyBM5Y5J9WDEi8JxmNMxiF2gn37Ydj/wj9x7C227aidcdnHaVuhjuirXSej7/1TNhPLNYHsEWz6p1rZl5M0tjOZHDobnK3iZuPfJUQW/r/XWiUa5+u4QUwbpD2TcVwqXNLEolWpHemMiq2DMRtM6DzTasFcAoJBsrShXBYicR+VREqdir6SYdR+SUUh3ArwD4OoCTAP6RUvp8p8/bD1qNyK32qmK1kHtT/syEkkVegPNwBUnwQfQR5LSiq0EGzSinHzY+eOulH1ZzYCaC0+bJGDdnrBJCLKGxe+zVUaLRAKnkuhc5Y/90GPcdmXOdfy4JPvgln6cR+UaLHjlg3MFc2spam2vxrIbxOqMM2YzU6ohcKxpDrN1E5IB31Z31Okx2QlAWrX2V63ZEoZdoTYdS5pFLgg+SQKxgJiC3aMfJYs3dBLs4tJpCuJ5WMR027hKGxSMHpfRfKKUHKaVXUUr/wIvn7AflnuTuI3LZtATseNU4iwnlZFi2IoN6RQ4ByWjuw9rAdmKt+HwEsuCzmnDVXR/zyOXGQn5wNoyzKymzqlOzWiowobFH9NXWCsuc8bdQENQuYaW1XjvN2MwYG5WtFGYdmAmDUmOgCGDYBPWslZmIAkKApSohd9tMrlzd6ZGQM2vFo/RDoPLYuH6nUTNi98lLJYqk7ZgKSIJVUNaqtRJWxJqW2CxIad1aKWAypCAkG3nousfNyarx7h0fASaCMghxX5WZV50LdKbC3lR3MqEMK0YfkqyaqyvkftkYBmHlHXcg5ICx4VloutlZhCL6mraXPTAbQUYtYimRr8imGAs6ROQRBRsZ1bp9rnfX0w2iLRaENWMzo1o2m1vsmSvXzEWQyut1h4tLgg/TYcUqlmG0LOQDHpEDRmuCa8x+NBvpgtXLJVUwmooxAQ7IgnUOKK3acfcdqqkhacdaKZUoNjMqpiLlACyrFRH1qALdiaEr0e8mouDDVFjBisvRTHmt5Ciskx41zmIReVAWLLsmWCcyZbMFNzIqJIF0vOFk9DhvfPBmVb3ueuywCfLHzm9WeLfliLxSyIslap2MrQ6V6IRwF4S8VYtrz2QQkkBwZjVtCXIja2Y+FqgZJdY3Ie+GR24eXzMRBTOR2rmuL5otDdh+QVAWLY+8lRJ9ALhl9zhesbdyU7ydPPJ4TkOxRK3NTsDdcIxO4EJexVzUb+WsNiNXp8NaSBagiM5jy1qB+XL2Ypl6FgOrZNxMG5tjnWb3KKLPRfphsaJrXD1uWIhhIiTjmy+uGtZKsFLI7U3Pqm/3y2Peun+oej04m7WwbQVJ8GH/VBhnVlJW6Xm9iBwA5qP+Go+8WedDxmRIho8Aqy6P92awVsGeRuTmRX5+zG/beypfeL5yYgkBScDrDhopzWzIivHvzo8ZpwZszbAyx8JKS+PqOoELeRWzUQUrLQl5rbASYqTRdZqjW47IRasZV7CetWJ65BttRIH1nq/ZZmdOLbqKyAUfwZ2HZvDoqTVsptXaiFx2EHLzvWM+fU8icqW17pfN2Mi4a2FbzXU7ozh+Ycs6fup55IDRX6Q6Il+z/l9jQbXuQJPeReSSQDy96IYtIQ/UNKRT9RL+5dllvOHwrHVXZz8evThmRMEHWfC1JuTWhC7ZCnTanfvpFi7kVcxG/a6FPN+gK99UWMZ6pxG5WhuR13u9gFT2yL0Rcl/zEn2XQg4Ad187g0ROQ6qg12x2hv2VWStAWYzqtUHoBhG/hLRH1gqlFFsZ1XUOuZ17r5vDVlbDvz1/BUDjfvvzY36kC3rFBei7Z9YxGZKtnjeNmIkqWE15FZEbE+i9rPVgx9fcmL+mId13Tq8hntXwkzfvqPl5wLuLvzGA2f1xwe7oZyKKtZ5uj3vjQl7FXNQYgeWmmKdRq8xJDyPykCJYglDv4GQjqdrxZZ3wi24ict11et+PHZiCJBgneLRayJUGEblp77TaNKsdjIjcmxMuntWgl2hbk+Rfd3AaQVnAl44vAmjukQPlFEStWMKjp1Zx1zUzrmaczkb8Hkbkuqf+OFC+W2MeuL39xVdOLGE8KOHHDpQrxe3nh1cX/1anBL2wlIQs+rBnMtR2HnqrcCGvgnVQW3VxcDdqXj8Rkq0Nu3axV05a1kodT9rPNjvThY5SDxmK5Gsq5Fm12DT1kBHxS7h9nzFSjVkFTpudIUVEQBL6Yq1E/SLSqm4Nse4E5vG304XSLwm465oZK50v1mAwRbkoyBDyJ85vIpXXcfe1s65ey4jIvRFyN/1dWoVNe2KFdZMhGRtpFemCjodfuIKfuGG+Ir3THpFXD2Bpl4Dc2nCJZy8ncHg+Cknw2XqqcyHvKbPmieFmwzNXJ/0QgDV/spOJM1lVt8rXLWulQR55Mq8jmdc7qupkGBF5Y2slV6etbj3uvnYGQFnAndIPgcrqzkIP0w/DfhGUerMxxQKBdtsJv+n6eQDG/kKkQcHVXLRc3QkA3zi5Clnw4ccOTLl6nZmIHxuZgidDmA1rxduMZnZHs3fKKIqbCivYyBTw8AtXkNdKeNtNlf35mJA71Xe0C8sIc0OpRPH8UtLKebesFb7Z2VvYieHGJzea19cX8oJe6uhKnCkYHjQhBLfsGcdd18xYqXzVBCTB2i2fcDFarBlu0g8zLtMPGW+5cQfeeGQON++OATDea5+tjzRjOqLUWCutXDDahbVo8CJzhY1Qa1fI7zw0g4AkIBZo7DnPRv0gxIjIKaV45OQKXnX1ZNNqW/v/p9Sb0YTJvPcR+W37JvDQr70W1+0whHEyLGMjo+LLTy1hZyyAW3ePV/w8Ox9brQRuRFASXZfZv7yRQbqgWwNv2OfQ7cZZXMirmI0aJ55TI6Jq8g0GvDIrxGkqvFsyBd2KVqfCCv76519Rt9AnKAtgwb8XEblhrXiTfsiYCiv45LtvtSpfd8QCePhDr8Pd18xU/Nx0WLGEJdfD9EP2Xnvhk7ML0UybQh6QBbzxyFzDYRiAEXlOhRVcSeRxbi2NCxtZ17aKfX1urMRmJHPe9SJnEEKsxmsAMBlSEM9q+N7Zdbz1ph01UbdVAe3hhb+VcW/Pmn1yqiPyduZ+tgKv7KxiLCBBEX2ufMN6o9cAVAxJsDe9aoWMqruOrOx3Bl5sdipi84jcbfphI66arp22Ph1R8PjLxnDhXlZ2etk4ay1VgF/y1dhGrfCHb78eerG5NTdvpiD+1ffOw0eAe66dafp/GLMt3IE2w4jIuyspLJe8WKL4yZtqxh5Y56OXeyohRcDluEshv5yAIvpwwKw8ZX3OeUFQjyGEYNahyKIaSmnDocAsct7sYMMzXShWtNRshH0dbifSNMLfpERf1Y3OjJ0KuRPTEQVbWQ2qXuppZWe0xcEijVhNFTAdUTpKxVNEwdWFfC7qx5MXtvD5H13EL756X0sdN9kdqD1wKZUo3vGXP7SGZ7shrxWh6iXPI/JqWPuLa+YiNR0fgXLbWS8v/AFJdO2RP7uYwOEdUWsgPCEEIVnkHnk/cFPdqRZLDcdJscyRzQ76rWQL7iNy+51Bq9WETvgloWHTLOYZBlqwVtwybZViF6yI3KsMhEaElc488tVU3urVsZYqtJV62A47YgGkCjp2xgL4tTccbOn/ToaVmurOrFbE4y9v4vtn110/TzfK851g7S+qNzkZLLDw0oqLBSVXiQvGRmfCslUYgRY2S9uFC7kDs2P+pmXLebVxpMgi8lZTEJ84v4m3f/IHyBR0pAu6aw+6Usg7P5kU0QetSOsOomaeYTci8ilbUVBeL1pDJ7pNJ+PeFreyeM0ffQsPPmnkfq+lCpiJ+Jv8L29gg0f+z5884vrCzxDMYR72XHJWv9DorvR//OuL+MKPymPmyuX53bVWbtoVw4d//CB+9o7djo93w1qZjSrIacWmF/iX1jPIqMUaIQ8pYtdb2XIhd2A2ouBKMt/wClxvOhAj6hch+kjL/Vb+4ltncfzCFk4sxpFViwg3GNpgh1krsaBk3dZ1Qnm4hPMB2E0htxcF5VVvhh64wRrA3IaQf/mpy1D1kjUUglkrveAdr9iFz733dtx5jXtv3M5s1F9R3cmEfLlBMPO1E0t45OSq9XWvInJJ8OFX7jpQ18IJyN4LObsgN9s3e37J3OhcqBTyoCzwrJV+MDfmR14rWQUZTjCBq9e8nhCC8ZBsDYJ1w+JWFo+eNgZTP385iUwb1ooXG51AOX2rXpl+zhLy7lkra6mCmeLZm8M0LLfnkVNK8eBTlwEAZ1fTKOhFJHJaz4Q84pfwqqvc5Y07MROpjsiNz7ZRF9DqtgDdaGHbDuWI3LtjhmX2NNsQZo8vjFcmN7RaGdoOXMgdmHGxk59zkU0xGZJbisj/4YlLIDCi+eeWEi1lrbBIxIuqTqB5RM42b7pjrRhd+ZbiubbmdbaLzxxBl2oxejqxmMBLaxlE/CLOrqWtpk7tph72mpnqiNz8bFdTeceBCGyWrN1qYEHPWJezVprRbABLOzA9aNbudyOjQhZ8NQkKxtxOHpH3HFYU1MgjtIS8gZCNB5tH5F95+jKO/N7X8alvn8M/PHEJrz80g9v2TeLEpTjyWsl1CTwTOy82OoFyb5N6Qm41s+qCkCuigIOzETy9mGjYz6YbRFroSc6st396chGy6MN7XrkXa6kCzppj7XoVkXfKjDnMg1V3MmulRJ2nBxXMjKUKIR+UiLwb1krUXa79VkbFeKi2gCukdH/cG88jd8BNdWe+iUcOGNPTTy4lG77Wc5cTSBd0/I9/fREA8DO37cbzSwk8ctJI/Wo02NgOi0S8SD0Eynca9YqCuumRA8DNu8fxz88s4aZdsZYnvXRCWBFdeeQb6QLu+fi3EQvKWE3m8YbDs1bF6g/PGTnwwyLk9urO+bFAhUAvJ/I16YxM6O0XvF555M3oxmZnxOz/08xaMRrW1X7mAUnkHnk/YFfgjoU8KDfNI9/IqNgZC+Avf+5W/NJr9+POQ9M4sqO8WdI3j9x8vnpFQVmX8zrb5ebdMaTyOk4uJz0tt25GxC8iVWjukZ/fyGIra5Skz8cC+IVX7bXGj/3wpWET8sqI0+7nOvnkzEO3X/CSOR2y6OuZDVYPtmfjZbdMQoir5mKGkNdeyEJKa0232oFH5A74JQERv9hw5mauSfohYIhqPKtBL5bqZpKwtrP3XjeHe6+bAwAc2dm6kLMucZMelOcDdmulcUTerR4ot5g9NNbTKq7b0cOI3C8h4SJlNJEzfub333odbtoVA2BUG8qiD88uxgGU0ygHHZaVwQKXTFVEDgCPvbSB3RNB7IiVI3a1WEJBL0IRBaPzYZ9tFaB8PHptx81Emvdt38youH48VvP9oCzyys5+MR1WGk4Xb5Z+CJSjYzayy4mNdG3/8NmoYlWwua3snIn48Udvvx4/ebNzoUSrKMxaqRuRd9da2T8VsnKSe5W1Arj3yONZcwybzUoQfAT7p0IoUeOzl7o4bNdLZqqqO1nELQs+XEnmoeolvOevf4T7v/OS8bht4469V70oz3dDWBFx9UwYh2z9WbxgJupv6pFvZlRMONRwBGUBarHkSYfJegzHkdYHmo1qs3qA1Ek/BCr7rdTDmLReKeSEEKvbWysFHu94xW4PrRUz/bDuZqdxAnerB4rPR3CTGZX3dLPTZdbKVtZ5niazV3pV1ekFrK0w87kzqg5F9GE+ZrSqeGE5iYJeso5ju4fO7JXkgETkgo/gkQ+9Dm802wB7hRGR19cDrViq20K6F8MluJDXYSoiN2ztWR4K3DwibyTkGxnnQRBHdhrtarvlQTej7JHXt1aCstDVisubTcuil76rEZE398gTWRWElFvfMiwhHxJ/HDDeX1n0WdWZrOsm6zn09MUtAOX8erv1wiLyRE6zLgijyGzUGKlXb2Qbq+B29sjZcInubXhyIa/DVFhpmDfqZpbkeJNWtllVR14rOV7FX3tgGhG/iB2x3pR5V9M8j7zzzofNuGXPeMVaekFYkZDXmt8Gx80ItHqcGhPyYckhZ0T9EhK5slCHFNHoqpjM4alLcQDlXPEKITc3hreyqietIQYVq91vHU3Yyhjvg1Ob6fLcTh6R95zpsIJkXq+btZHXixB9pKEPyiyTekLOZg86pQzevn8Sz/zej1tNgnoNa1JVb7OzlXmd7XLTQgyEdM+HdyLiskw/ntUcp9QfmDG82WGKyAGjkIdZK2lzoMncmB8riQKeuhgHUM4VT9sEib1P8YxmjfAbRawy/TqZbGwgtJO1yTJputk4iwt5HaZYB746mSs5tdTUu2Unel0hN79frxrTy2nkrdI8/bCIoNRd22csKOGTP3sL3nW7c4OkbsCEPNnEXonnnIVr31QIV02HrJzyYSEakCyhzqqGtTIX9UMtlnBxMwtCbB56lbWiFUtIFXTP9mcGEZaiudIkInd6D3aaA7Lv/+5LdZvQdQoX8jqw1LF6PnlOKzYtVFFEARFFrCvkmw2u4v3G3ywi14pWymM3ue/IfE3vim7C7DCWlVKPRFatyFhhyKIP3/j11+O+I95utnWbqL8s5JmCjqBprTBuXIhVeOiMdEG3/OHRtlYaR+SNzuXDO6L4zfuuwddOLOF3vvJcR3N868GFvA4s/a+ekOe1Yt2GWXbGzSHMTrA8da9yv71EFHwQfKRh98NeWh69YtzcrGpWyGVE5KMjXGMBm0dudt2cMys6BR/Bq6+eRM4cHpEu6LaWv1o5FXOErZVoQIQi+urum20yj7zOe/D+11+F97/+Knzu8Yv42jPLnq+vo3tjQsgfA3gLABXAOQC/QCmNe7CuvmPvwOeE2x4gEyG5bk9yJvBeldV7jV+sP7czU9AxHuxdpNwryhF5EyHPao4R+bASDYgVm5khWbRaVVw7H7FGwqXymvnZG8PFUwXd6ifkVZ+fQYRVd9ar9t7Kqoj4xYZ7Zr9x7yEcno/iJzxOjQQ6j8gfBnCEUnoDgNMAPtr5kgaDsrVSxyNvQcjr+eybGRWK6BvYyNYv1Z/bmdNGNCK3Mo3qWyvFEkUyr2FshISLReSss2FIETEdUSCLPtyye9zKEU/mdWMEoSIiYvalsawVh9S7UWIm4q+btbKRUZt2HiWE4C031g6M9oKOhJxS+hCllBlmjwFY6HxJg4FfMvztehF5Tm3ukQOGRVPPntlIGx9+Pzc1G+GXhIYl+m4beg0T0YAEH6mNyDczKt7798ewkS4glddAKUYrIvdLKJYoMmoRGfOzFXwEn/nF2/CBuw9YVZvJnGbmmQtWFSwrjhrliBwwNjzrpx+qjqmHvcJLj/wXAfxrvQcJIe8jhBwjhBxbW1vz8GW7x1REqe+R682zVgCjkGA9XXDs67yRKfQtvdANiuirW6KfU4sIdDlrpR8IPoKxgFRjhz11cQsPv7CCx1/etHnCIyTk5kVpLVVAsUStlLnb909iKqzYInLN6pMf9otVm52jLeS7xoM4v57BN07WDqV2E5F3k6ZCTgh5hBDynMOft9l+5mMAdACfrfc8lNL7KaVHKaVHp6envVl9l5muUxRUKlGk85prIa/X15k1zBpUFElwLNGnlCKj6iNprQAwJztVWitMvC9tZq3eOaMk5KwqczmeA2D0LLHDhD6ZMwZKGNaKZFgrGRV+ydf1uoJ+80uvuwqHd0Txvs8cx5eOL1Y8tpVR+3ohaxpSUUrvafQ4IeQ9AN4M4G7ajbyaPjIVkXHqSgoA8JEvnsBqqoC9k0F87+w6zq1l8GMHml+Q7EMqqvs6b6RVXD0d9n7hHuGXfI4l+gW9BEq71/mw34wHazeo2deXtrI4lDWKfsYCg3sRbhUWcS+Z3Q6re/xUROQFHWFZhOov4dKm0c531KNxwNjv+tx778B/+rsn8NEHn8Wbrp9HQBZAKe17UNaRtUIIuQ/AbwJ4K6U0682SBoepsIL1tIrz6xl88fgizqyk8I/HFhFSRHziHTfht990bdPnmBurP6RiI1MY6IjcLwqO6Yes+Y/bzozDxniwNmWUpeZd2sxZ/x6liJx54Cwir/5srUKpnIaMfbOzoCOeVUc69dBOWBHxn16zH2qxZA1bzqhFqMVSX8/lTk3O/xeAAuBhc8PuMUrpf+54VQPCVFhBIqfhoReuAAC+8L5XYtdEoKXNyVlr2lCltcL6rAy0Ry75kE7XlqpnrXmdo+eRA0Zhy3OXK60Ve0Tu1MJ22GHWSr2IPCgbm5+JnOGRhxUBxVJ5s3OUi4GqudFs5vb0pTiO7p0op18Oq5BTSq/2aiGDCMslf+D4IvZPh7B7svW86cmQDEkguFIVkVt9VoY4Ih9Va2UiZEx2opRaF20m3oubOau1wih1+2PWyXLCjMirMpIIIYj6RVxJ5kGpIfRFaszt3MqouHZHtOdr7hfTEQU7YwE8bTYTa9Zqoxfwys4GsFzy0ytp3Hlopq3n8PkIZiL+mpFZ7NZ9kK2VWFCyUsvsdHuoRL+JBWWoeskaHgKUhVwtlnBmJYWIX6w79WkYYdbJctw5IgeMDU/74xEzZXEpkdtWETkA3LQrZgn5IETko3MkdoEpW8Vlu0IOGBNYaiJyszfDoFZ1AsaFbDNTsBr9fPXEEj74hadG3lphPaXtF7F4ToVs9p959nJipPxxwGjJEFZELLGI3OGzjfolK2IPK6KV2ZLXSttis9POTbtiWNzKYT1d4BH5oMMi8qAs4BX7xtt+nrmov4G1Mrge+VRYRomW/eFvnFzBV55estqajnJEDpQjLcCIyA/PG/bB4lYOsRHKWGFEbWPunCNy0ZrhaUTk5Z/ZLpudDOaTn7gUxw/OrkMWfH1tXcyFvAHsg3nN1VPWDMt2mI02sFYGOSKPVHaAZDMLv3jsEoDRFXIWXdpTEONZDdfZfOBRi8iBcq44UOuRA0ZEztJRQ2ZlJ8NpMs4oc2RnFIKP4LOPX8SDT13GL7xmb1/vULmQN8AvCfgvr78Kv/S6/R09z9yYHxm1WDFCbDOrQhZ8A53CZ/WbSRmCxqaIn98wMk1Hd7Ozso+8VjQ6/s1G/VZdwChtdDKYkIs+AtnB/7fP5DSslfLX2y0iD8oiDs5G8M0XVzERkvHLd/Y374MLeRN+475rcOueiY6eYy5am0uezOmIBqSB7bMC1PZkX00WKiLRUfXIY1U9yeNWLxEJuyYC5s+MnpCzi1NIER2PS5Zrzn7GHpFvN48cMHxyAPi1ew70ffA0F/IeMGtVd5ZzyZM5DWOBwRbCaZu1klV1pAo6fvrWBWtO5ahaKyw/nEXkiZyZbhiUscsccjGaHrnxe1eX51c/zn7G/nMT21DI/8PRBbzr9t34mdt6N8GqHoOtJCOCU3VnIqdVeJKDSNQvQhZ8WEsXLH/80FwUr756Ct87s2bN9Rw1RMGHqF+0OiBu2SLyhQlTyEcwImcRd70LdKWHLkIRy+0bYtvMIweAm3eP4+bd7SdBeAkX8h7A5v3ZM1eSeW2gc8gBowhkKixjPaVa7Ttnowo+9IaDeOX+yYG2hTrFGAhSaa3EAjJ2m0I+ih653Vpxwm6tBCUBJfNCLvoIInX+D6c38He/BwRlw0+sjsj3Tob6uCp3sFa+bKNzJuLHobmI5Q+OKjFb4yz2dywoYd+UIeRTA9xaoV2YdVKvz7z1uCzA5yPwgSAgCQgpwkhf1IcBLuQ9Yi7qx5WEfbNTG4qobiqs4Eoib/WKmeljrmwvmQjJ1sUrYes/vjAewP0/dytee3A4WjG3ArNOnIqBKh63Rd8RvzgUx/GoM5om5wAyN+a3InJKKZJ5veJWdVBhE45WU3nIgm8kvWEnYkHJ6kkez6kQfQRhM5vjx6+bszZ8RwkmyM02O+2Ph/3itsxYGTQGX0lGhNmoH6dXjN7mGbWIYokORSQzFVawkVGxkshjOqJsm1voiQprRUMsONipol4QNdMJg/WsFTPwsEfkd+yfHOjGb9sFLuQ9YjqiYCNtdNRj/az7nXvqhqmwgmKJ4tRK2tq03Q6Mh2Rk1SLyWhGJ7HDYYJ0yFmyy2engof9fP3V99xfGaQq3VnrEZEiGXjJEPGkK+TCIAyvTP7uawkzE3+fV9I5xW1HQVra/Y7x6RXkz01nIWU/yetYLp39wIe8R5SpJtRyRD4GQT5vr1ooUM9soIp83c/9fWksjbloro85ESMZsVMGBGefxg4QQRPxi3Yid0z/4J9IjmJBvpAtDFZFPR8qRKKtQ3Q7ctm8CsuDDt06tIp5VcXgbDE7wSwIe/+2GI3rxjqO7cGTnWI9WxHELF/IewfqOr6dVq5/3sHjkjH626ew1IUXE7fsn8M0XVxHPaSM11q0TPupiTi2n93BrpUdYEXmmgKTZ83kYIvKxgARJMLI1tlNEDhjDRM6tZZBVi32d/sLhNIMLeY8YD0ogBFhPFSyPPOwf/BsiQog1/GK7FAMx7rqmPBVqGC66nO0LF/IeIQo+TARlrGdUJHMaIn5xaIpKpkyffLsJ+d6pEPZPGW0UtsNmJ2d44ULeQybDsrXZOQz+OGMqrEASyLZIwavmTjMq346/O2d44ELeQ6bCCtbTKpL54SowuWo6jP1TYfiG5A7CS9520w5MhmTsnx78Bmec7cvgm7QjxGRYwbOLcfgIhqLPCuMj9x7CB+460O9l9IUbFmI4/jtv6PcyOJyG8Ii8h0yFZWyYBUHDFJH7JcEq3+ZwOIMHF/IeMhVWkCroWEsVhsoj53A4gw0X8h7CusRtbZMmTBwOpzdwIe8h9irJYeizwuFwhgNPhJwQ8mFCCCWETHnxfKMKK9MHeIEJh8Pxjo6FnBCyC8AbAFzsfDmjTWVEPjxZKxwOZ7DxIiL/UwC/AYB68FwjjV3IeUTO4XC8oiMhJ4S8FcBlSukJj9Yz0gRkASHZmK7ChZzD4XhF0/t7QsgjAOYcHvoYgN8G8ONuXogQ8j4A7wOA3bt3t7DE0WIyrCCzmeXphxwOxzOaCjml1LHTPCHkegD7AJwwh9IuAHiSEHIbpfSKw/PcD+B+ADh69Oi2tWGmwjIubmZ5RM7hcDyj7R03SumzAKw+n4SQ8wCOUkrXPVjXyDJp+uQ8/ZDD4XgFzyPvMVNhGbLgg18Smv8wh8PhuMCzHDhK6V6vnmuU+ZnbduPQbKTfy+BwOCMET2buMTcsxHDDQqzfy+BwOCMEt1Y4HA5nyOFCzuFwOEMOF3IOh8MZcriQczgczpDDhZzD4XCGHC7kHA6HM+RwIedwOJwhhws5h8PhDDmE0t73ryKErAG40OZ/nwIw6P1c+Bq9ga+xcwZ9fQBfYyvsoZROV3+zL0LeCYSQY5TSo/1eRyP4Gr2Br7FzBn19AF+jF3BrhcPhcIYcLuQcDocz5AyjkN/f7wW4gK/RG/gaO2fQ1wfwNXbM0HnkHA6Hw6lkGCNyDofD4djgQs7hcDhDzlAJOSHkPkLIKULIWULIbw3AenYRQr5FCDlJCHmeEPJB8/sThJCHCSFnzL/HB2CtAiHkKULIPw/iGgkhMULIA4SQF83385UDuMZfMz/n5wghnyeE+Pu9RkLIXxNCVgkhz9m+V3dNhJCPmufPKULIvX1c4x+bn/UzhJB/IoTEBm2Ntsc+TAihhJCpfq6xEUMj5IQQAcCfA3gjgMMAfoYQcri/q4IO4NcppdcCuAPAL5tr+i0A36CUHgDwDfPrfvNBACdtXw/aGv8MwL9RSq8BcCOMtQ7MGgkhOwF8AMaA8SMABADvHIA1/i2A+6q+57gm89h8J4DrzP/zF+Z51Y81PgzgCKX0BgCnAXx0ANcIQsguAG8AcNH2vX6tsS5DI+QAbgNwllL6EqVUBfAFAG/r54IopcuU0ifNf6dgiM9Oc11/Z/7Y3wH4yb4s0IQQsgDgJwB82vbtgVkjISQK4LUA/goAKKUqpTSOAVqjiQggQAgRAQQBLKHPa6SUfgfAZtW3663pbQC+QCktUEpfBnAWxnnV8zVSSh+ilOrml48BWBi0NZr8KYDfAGDPCunLGhsxTEK+E8Al29eL5vcGAkLIXgA3A3gcwCyldBkwxB7ATB+XBgCfgHEwlmzfG6Q17gewBuBvTPvn04SQ0CCtkVJ6GcCfwIjMlgEkKKUPDdIabdRb06CeQ78I4F/Nfw/MGgkhbwVwmVJ6ouqhgVkjY5iEnDh8byByJwkhYQBfAvCrlNJkv9djhxDyZgCrlNLj/V5LA0QAtwD4JKX0ZgAZ9N/qqcD0md8GYB+AHQBChJB393dVLTNw5xAh5GMwLMrPsm85/FjP10gICQL4GIDfdXrY4Xt9fR+HScgXAeyyfb0A49a2rxBCJBgi/llK6YPmt1cIIfPm4/MAVvu1PgCvBvBWQsh5GHbUXYSQ/4XBWuMigEVK6ePm1w/AEPZBWuM9AF6mlK5RSjUADwJ41YCtkVFvTQN1DhFC3gPgzQB+lpYLWgZljVfBuGifMM+dBQBPEkLmMDhrtBgmIX8CwAFCyD5CiAxjs+Gr/VwQIYTA8HVPUko/bnvoqwDeY/77PQC+0uu1MSilH6WULlBK98J4z75JKX03BmuNVwBcIoQcMr91N4AXMEBrhGGp3EEICZqf+90w9kQGaY2Memv6KoB3EkIUQsg+AAcA/KgP6wMh5D4AvwngrZTSrO2hgVgjpfRZSukMpXSvee4sArjFPFYHYo0VUEqH5g+AN8HY4T4H4GMDsJ7XwLilegbA0+afNwGYhJEtcMb8e6LfazXX+3oA/2z+e6DWCOAmAMfM9/LLAMYHcI2/D+BFAM8B+AwApd9rBPB5GJ69BkNs/mOjNcGwC84BOAXgjX1c41kYPjM7bz41aGusevw8gKl+rrHRH16iz+FwOEPOMFkrHA6Hw3GACzmHw+EMOVzIORwOZ8jhQs7hcDhDDhdyDofDGXK4kHM4HM6Qw4Wcw+Fwhpz/HzpPhpw3gVg3AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data = hal.get_data_linker(\n", " {ma_graph.b: np.ones(q),\n", " ma_graph.c: 0.},\n", " n_data=150\n", ")\n", "m = hal.get_generative_model(ma_graph, data)\n", "plt.plot(m.get_example(ma_graph.y))" ] }, { "cell_type": "markdown", "id": "4fe97dce", "metadata": {}, "source": [ "## the AR model (autoregressive)" ] }, { "cell_type": "markdown", "id": "120c1b0e", "metadata": {}, "source": [ "An AR model is defined as\n", "\n", "$$\n", "y_t = c + \\epsilon_t + \\sum_{i=1}^p a_i y_{t-i}\n", "$$\n", "\n", "which as an Halerium graph is" ] }, { "cell_type": "code", "execution_count": 5, "id": "8aab7b97", "metadata": {}, "outputs": [], "source": [ "p = 7\n", "sigma = 1.\n", "sigma_obs = 1.e-2\n", "\n", "ar_graph = hal.Graph(\"graph\")\n", "with ar_graph:\n", " hal.StaticVariable(\"a\", shape=(p,), mean=0, variance=1)\n", " hal.StaticVariable(\"c\", shape=(), mean=0, variance=1)\n", " \n", " hal.Variable(\"y\")\n", " y.mean = sum([\n", " a[i] * y[hal.TimeIndex-i-1]\n", " for i in range(p)\n", " ]) + c\n", " y.variance = sigma**2 # since there is no moving average in the AR model\n", " # time we include epsilon as the variance of y.\n", " \n", " hal.Variable(\"observed\", mean=y, variance=sigma_obs**2)\n", " # this \"observed\" variable is strictly speaking not part of the model,\n", " # but it is recommended when we want to fit data to the model as it allows us\n", " # to include observation noise." ] }, { "cell_type": "markdown", "id": "120ea63c", "metadata": {}, "source": [ "Let's look at an example for the AR process" ] }, { "cell_type": "code", "execution_count": 6, "id": "195dd345", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABhfElEQVR4nO29eZwk11Xn+7ux5VZ7Ve/daqml1i55k2XjRd7HwjAwLI8BZoAZA3oMMINnYVj8eczyHvMZYGB48wYGNOyMMB/GCzbGNl5k2eBFsiRbm9WSWnuvVV175RbbfX/cODduREZkRlVlVmVW3+/no4+6qrIyb2VGnPjF75x7DuOcQ6PRaDSji7HbC9BoNBrN9tCBXKPRaEYcHcg1Go1mxNGBXKPRaEYcHcg1Go1mxLF240Xn5ub4lVdeuRsvrdFoNCPLQw89dIlzvi/9/V0J5FdeeSUefPDB3XhpjUajGVkYYy9mfV9bKxqNRjPi6ECu0Wg0I44O5BqNRjPi6ECu0Wg0I44O5BqNRjPi6ECu0Wg0I44O5BqNRjPi6ECu0WgGxun5dfzdM5d2exl7Hh3INRrNwPjt+57FXX/6INp+sNtL2dNsO5AzxsqMsQcYY48wxp5gjP2HfixMo9GMPvW2j4Yb4P7nlnZ7KXuafijyNoC3c85fAeCVAO5kjL2+D8+r0WhGnJYXAgDuPTW/yyvZ22w7kHPBRvSlHf2n58dpNBo0PWGp3HtqHnqs5ODoi0fOGDMZY98AMA/gM5zz+zMecxdj7EHG2IMLCwv9eFmNRjPktKNA/tJSA88u1Hd5NXuXvgRyznnAOX8lgKMAbmeM3ZzxmLs557dxzm/bt6+jC6NGo9mDNL0Atx6dBADce+riLq9m79LXqhXO+QqA+wDc2c/n1Wg0o0nLC3FirobrD47jc09qn3xQ9KNqZR9jbCr6dwXAOwGc2u7zajSa0afpBSjbJl5zfBpPX1zf7eXsWfoxWOIQgD9mjJkQF4a/4Jx/vA/Pq9FoRpxWFMhNg8ELdLJzUGw7kHPOHwXwqj6sRaPR7DEokAOAG4S7vJq9y2Wxs3N+rYUf++MHsdr0duw1vSDUu9k0lzV+EMILOCq2Ccdk8IJQlyAOiMsikD/44jI+++RFnDq/tmOv+e8/9gR+/E8e2rHX02iGjZYvFHjZNuBYBjgH/FAH8kFwWQRyUuINb+cU8pnlJl5a1HWzmsuXVnS+VRwTtilCjaftlYFwWQTylUYUyNs7F8hdP0Td1daK5vKlGR3/ZUsJ5L5W5IPgsgjkpMjrrr9jr9n2A9TbO/d6Gs2wQTmismPCtkSo0QnPwXCZBHIXQKwQdoK2H6LhBgi1J6i5TGm6kUduGXBMBkBbK4PisgjkZK3spCJ3o0RPcwd9eY1mmGj5sUfukCL3dSAfBJdFICdrZacVObCzFw+NZpigZGfZ1snOQXNZBHKpyHc42bnTr6m5PHj+Uh0feujMbi+jJyScKkog1x75YLgsArksP9zhZCcAnfDU9J0//9pL+LkPPbrby+hJoo5cKnKdMxoEl1kg33lrZSdfU3N50HQD+CFHMOSJ9Jbbaa1oj3ww7PlA7gUhNiJVvJOKPLZWtCLX9Bfynoc9KFKys2zHyU7tkQ+GPR/I15T+KqSO11se7vqTB3FxrTWQ1/SDUG5F1slOTb9pRnMwhz2QJz1yUX6oPfLBsOcD+YoSyGmn5akL6/j0Ny/iGy+vDOQ11YN1J3eTai4PpCIf8qBIg5cTVStDfvEZVfZ8ICd/vGKbaEQ2B6n0Qd3mqUpJK3JNvxmVQN70AjimAdNgirUy3L7+qLL3A3lUenh4qiytlbVW/wP5y0sNnF1pAogTnYD2yDX9hyyL3bRW7v7is/jEY+e7PqblBSjZIsQ4svxQ36EOgr0fyJsUyCsy2bnWFP/vZwOff/vBR/FLf/k4gLQiFwfuBx86g3/7wUf69nqayxdKIu5kIH/4pWXMKzml3/+75/Gxb5zr+jstL0AlGipBvVZ006zBsOcD+UpD9Fk5NFmWQXU9UuT9vDVdbrhYjl5LHShBds59T83jM9/UU8Q122c3FPl7/+hr+O+fPw0ACEOOxQ235+AUdTqQTnYOln7M7BxqKNl5cLIC1w/hByHWWpEi7+NB1fZDGEwcrJTkAWJFvtLwdN8VTV+g42ungqLrh1hpeDi/KhT5atODH/Ker99UFLmjt+gPlD0fyFebHsZLFibK4k9teMFAkp0tLwCL/q0e4OSRLzdctLwQYchhGCzjGTSaYux0HTnZkwvrbQDApQ3x/7bX/fVbXohy5JHrDUGDZc9bK6sND5NVG1VHBPKmGyjJzv75dU0vkIq7naPIgWQiVKPZCjtdtUL2ZBzIyULsrcjJWtEbggbL3g/kTQ9TVRtVRxxQ9bYvk539VActNZBH3mHViUseyT/X9opmO3DO5TG0U+p2RVHknHOpyHu9flsJ5JZBHrlOdg6CbQdyxtgxxtjnGWNPMsaeYIz9TD8W1i9Wmh4mK3EgbyQUeX9OBM45Wl7Yccs7XXWw0fbR9gNZ+rjZQP7swgbe/uv3yZNnO3zwoTN4/Ozqtp9Hs3u4QQhqsbJT6na57srXXmv6sbXSM9kZSo+cMQbHNLQiHxD9UOQ+gH/NOb8BwOsB/BRj7MY+PG9fWGm4mKo40lppuAHW+5zspFtM8sDp6+majYYbSFsF2HxP9CfPr+G5hTpeuLT9Qc7/8a+ewD33v7Tt59HsHmoifccUuXL8Lmy0lEBexFqJQ4xtMr2zc0BsO5Bzzs9zzh+O/r0O4EkAR7b7vP1iteljomKjWoqsFddXkp39uc1rKSq77YcJRd5wfWmrpB9bBNriv11LhnOOhhts+vU1w4X6+e2ctRIfv/PrbSxGHnmv1295ASrRnTAgasl1+eFg6KtHzhi7EsCrANyf8bO7GGMPMsYeXFhY6OfL5sI5x2rTxVTVRo0UeTu2Vvp1UKlBtukFUqnM1BzU2wGW6m7mY4tAm5i2O93IjRp57eSUJE3/SYiGnbJWVEW+3t6UIi9ZcSDX1srg6FsgZ4yNAfgQgPdxztfSP+ec3805v41zftu+ffv69bJdaXoBvIBjSvHIl+ptqcT7dZun3u62vEB6h9NVB00vFcg3GUgbXn8Ueb+UvWZ3ae6GIm94KEVVJwvrbSxsdG58y6LthUlFbhpw9c7OgdCXQM4YsyGC+D2c8w/34zn7AXl7arKTNjUA/fPI1eDc9AJ5gs3UHADA2eVm4ufdWG14+PMHXgLn4oCXAbjHBaDpBnjdf/os7ntqPvPn/bogaHYXVTTslLpdabg4NlOFYxlCkUdliF7AEeYMtwiiDUNlVZFbWpEPin5UrTAAvw/gSc75b2x/Sf2DAvlU1UatJKyVC0q/CK9PE1ZaijJpurG1Ml21AUA20wJ6e+R/+OXn8fMffkz+DlW79Jo0dGmjjYtrbTxzcSPz51QGqT3y0Ua9oO+kIp+u2tg3VhIeeb2NaBNzrj1Jx1nFSSU7dSAfCP1Q5G8E8EMA3s4Y+0b033v68LzbhnakTVRslCwDjCExTKJ/1kqQ+HfbD2AZDBOVKJCrirxHQP7C0yJ/QJU1TY/+37uvBQA5DSkNbUzSHvlosxvJzuWGi8mKg33jJbywWEfLC3FgvAwgf3cnHa9URw4IRb7VNV9YbeHO3/wizimiSBOz7S36nPO/AzCUe84pkE9WbDDGUHMsXBiAtZIM5KJqxbEMWfJ4ZrmJim0mdn9msVx38Ug07KIux9MVC8Cxcs8O5I3U82lGk0Qg3zFrxcMtR2wYDPi705cAiLbQF9ZakU9u565TDeS2ufWqlWfm13HqwjpOXVjD4anKlp5jL7Ond3ZSMoY2JVQcUwby6ardx/LD+OCkqpWSZaAWlTyeXWni0FS547Fp/u70JbnZgxR0vWCSsikVefbjKIBra2W02ZVkZ1T5tW+8JI8jCqZ5lSt5gXyr4onOm7zj+3Jnbwfy6MOnPg81x5QBcnas1L/yw1Sys+2FKFmmLHncaPs4OFGWP8+DbBUgVuRkrfRS0s0eipwmFelk52iT2BC0A4q85QVoeSGmqsJaIY5M9wrk4vsVO11+uDXxRKJso6UHtWQx0oH8i08v4CfveUhWeKShOluqZSWrAwBma07/rBUl2dlyA7iBsFZIkQPAdM1BxTZzFTHnHF94egGvPDYFIPa6iyppCtB5E4nUFgF575dm+KHPeaxk7Ygip4KB6XQgl4o8+7iMPfJksnOra44VudfjkZcnIx3IP//UPD7x2AUEOdUndNCQIqcSRMc0MF62+uiRK3Xkvkh2lhSPHBBWTsUxc73uJ8+vY2G9jffcchCA4pEXLD9spKyYNPR8nOsOjMNOw/Xl8JM0dEGfKO9MIKddyVNVG/ujBCcAHJoUgTxvDbJqxe5P+WGrh3V4uTPSgXw+qmfNu8UktUCbGapRCeJExRJ+XZ82J6hquekGMtlJJY+AUDSU8MzioReXAAB/70YRyGXykqyVooo8x1pRLwTaJx9ufumjT+Cf/a+HM3/W8gIwBtR2SJGrgZwU+XTVlnebeaKAjrd+JTvpdbS1ks1IB/KFtWhjQk5Aloo8ampfjQ6qibK9rcRLGjq5gGSys6rsapuqOijbRm4gX48U88HJMmyTSWuFTohWT49cPD7PWqmnfHzN8HJxrYXzq9lldjQHc6c219Dw8qlKbK3MjpWkOMr1yKPvlzs88u0pcj3MPJuRDuTz66ICJe8q7/ohbJPJiTzUOGu8Ym9LHaRputHJZRpoeaFMdtqmIW2dmZqwVvICskzMmkLJp8sPSZnnr0H8fp61oiZBdQnicOMFYW51U1MJ5EWP3/ufW8QP/s+vysEQm4H6rEzXbMyNiZ3Kc2OOzDu1c0RBy83yyLd+FywVuQ7kmYx4IKetwvmBnNQ4AFlFMlG24Fj922XW8kUD/bJtiA1BUbITEEkpQCjybtZKyw/gWAYMQ9S7b7R9hCEvXkceBfo8a0UN8HpT0HDjBTzX/mq6IcqRaCiS6/jU4xfwQ3/wAL787CK+Ee1R2AzU+XC6KoL3VNXGnKLIc3d2+tke+Zatlej9WNeBPJORDeT1ti+DXF5AbvtxQAXiZOdEhayV/tWRV2xTJjPbXhD78tFrTlcdlJVAvtr08NiZeMiDUPFRmWTJRL3td2z977oGN771zKpKaSqKXnvkw40XhF0v+GXbKLRL8pvn1vCT9zyEE3M1ANjSrkhqmEUWyc/deT1++FuulOdV3s7O3DryLVetbN1a+eBDZ/D0xfUtve6oMLKBfF65TeymyNU2mlVFkW/noErT9AKUbEMqbtcPUYoOYLoLmK7a4udRwP3DLz2P7/vdr8ig2/bjsVi1koWGG08VMlhvX5seG+ZUpSQUuQ7kQ43ri2lTWRfklhtIRd4rkD9/qY6QA7/+fa+AYxlbDOQupquO/PoHbr8Ct181E1srucnOTo/cttiOJzs55/jFDz+GP9vjA1VGN5ArPVPyWmO6QVKRU6adkp398sjbXoCyZaIc1Ym3FUuHfPmpqiM88iiIXlxri6AfrUFV5GMlYa1Q6eFMzenpa6tVLVk+YsP1MVGOB1BrhhcvGueWdXy2/OLJTj+M91EcniwnmrcVZbkhZt6mkdZKTh15yw/gmAZMI+7e0Y9k52Y98o22DzcIZe+ivcroBvICirwd+c5EJWGt9M8jb0aTUMg6afshSnYclE2DYaJsJTzy1ch7pKRWS1Xkjkh2ku89U3PQ9sPclqFAsqqlkZHwrLcDzI2V5Ho1wwtZflkJzyYp8gJ+Mz2PYxo4PFXZsiLPDOR296oVsc5keLFNAyFH7r6PbsQbgjYXkGlD016vdtnTgbxbsnM7B1WalheiHFkrLS+AG6kRQHjkU1HTrrJirVBDL0riqIq8WjJRb8fWymxNBOBWl0b+qmLPOtibXoDZqOpAe+TDDVkmWZ+TONaKWSt0XlgmiwJ5q+vjsxAtbJ2O79PxnRfIVatQ/o5U8ZsXUHKLfk4OKA8a6pJXBLBXGOFArlgrXZKdJTtPkYvv90OVU21vxelU5K+9cgZ3XCsmIgnrRbweKYUsRU7WCgV9CsDd7JWmF8jb2Kx+K/W2Ly8I2lrZHb7vd7+C3/nCsz0fR8dkdiAvnuz0U4H84npr08d7nrViRbZJ3hqabnJeJwB5zm3F0qTzJFCGmxdhKdrQtNcV+bbb2O4Wak1s3sGUr8htLJjxrtC0ctgsItlpgnMuB0tQMujH3nxCPq5im2J2ZhDKQN5UFDndiopkpy8Pvtlo0lC3ANx0A8zUHCyst3M88gAzdEHQinxXePLcGkqWgZ94y9VdH0eBLssCozpyu5Aij62VI1NlcC76eh+bqRZar+uHWGm4UgCkKVlGbq+Vl5YamEopecdk0bq2rsgB0au/6Dm7TIp8j2/tH1lFvrDelkE6r4wwXX5467FJfPerjuDVV0zL7/ejcqXthTLZSUmVktX51tK0lJYfSmuFVFfLjwfVjpUseAGXj5kt4G03vdgDTyt3zjnqro+Z6MTqtUtUMxiaXoDnFuo9Hxcr8oxkZ5SPKRXwyCnZaUUeObC5EsSnL67DDzmuPzSe+XPHyq5lf+FSHQ+/tII7bzqY+P527oLV92IzPjlZK3t9I9HIBvL5tbbs8d29/DD+EyfKNn7jH74Sk1XVWtmcR35updlhXYhkp/DIV6LgmxnIIxWx3vLkgdXKUuTRLSnlAcha6abIG24gd96lD9q2H4JzofS7bUrSDA4/COGHHGczjp80dExmfd5NL5ns7OYX0/NYBpPdCs/lbP3P4tFon8MtRyYzf16yjMw68g89fAYGA77rVUcS35fWyhbEU8uP92ZsxiahXjHaIx9S5tdb8uCkQL7a9PDxR8/Jx7hBso5cZbPq4PlLdfz0nz2MN/zne/Ebn3468bNWVH5Ysc2OjosqdDt4cS22hagnharIqbkX2Ud0a9vNI295AfaRIk8d6HTg10qxj6/ZWVpK8OqmyoOQywR8OrnNOU8kOzkH/C7Jejq27YQiL57wfOzsKibKFq7IsWJKltlxVxCGHB9++CzeeM0cDk6WEz+Td8FdzrmWF+Az37zY8f22F8o7zs2UEi7rqpXhxfVDLDc8GcgpeP7VI+fw03/2dRkA0+WHKnbk1xVNvPzQ79+Pe0/No+aYiXpccXLF5YdEtrUifn5htXOGp6rIaVv/wkZSkXfrZd5wfcxFTY3qqYBPF4CKbUabknQb251G/eyeXcgekA0kg1zaAiMbg5KdQHd16wccBgNMQ1RMzdacTdWSP352FbccnQRj2ZMcszzyrz6/iLMrTXzva452PD5W5PkXn7965Bx+/E8exMtLjcT32358x7kpRR5ZK17Ac/38vcBIBnIKcDSlJH0rSreu6WSnitNDkV9YjTP8XhDizHITP/GWq3H9oQnpXQPiQhByobbVLH3WnQBZK+eVuaF0cLW8WJFT+9tL623YJsNE2Y7+ruwDkdYwWbFhGazjQKfbylrJSmxK6jfnVpr4jU8/pQdXZJAM5PmKXBUWaUWu9vguckfphSEs5fjfTC152w9w6sIabjkylfsYJ8Na+ejXz2GsZMl2zMnH9052kp2oVqUBwiOnXNFWPHJgbyc8RyqQN1wfT5xblYo7ba1QUCTrwE2VH6rIEyFDHbS8AO/49fvwgQfEtl7y2aZrDiYrNtaUpv+UhFH7UQDZ1goFcnUAtPTIE5uIxOMW1ttSRat/V5qmorjVzokEXQCqjjlQj/yzT17Ef7v3NC6sbb5eedTwgxC/97fPJY6FbqiB/LluilxR2Ok7J3UyfRFF7vkctrKz8vBUGWeXiwXypy9swAt4rj8OIDPh+vylOm4+MtFReggUszNJQS9uuInvt/1AVm9tpnEWnbvA3rZX+hLIGWN/wBibZ4w93o/ny+P9H3kcP3D3V/HV5xYBdAZyOqgpsHVT5HaX7m3za23U3QAvXBK3d/G4KxuTFTuhyKVKcszETrYsa6VM1sqaGshFwkotWaSeMAvrbVQdS54UuYFcWYM6l5SgnZ4y2TmgqpWWcgHd6/ztM5fw//z1k/jck51+bhZ0wTcN1lWRq8n39J2TesGWTau6WSs5irzIHdOjZ1cAALce7RbIzQ5FvtJ0MVXp3EAExHfB3Y4PUtCqkg5CDi/gUpFvJiAv1eM6+L2c8OyXIv8jAHf26bly+VfvuhaOZeBXPnUKQGytyH4lFMgVldvLI89SBwsbItAu1dvR/8VBNVN1MFG2ZLN9QOnyZpmJlp1FrRXaQASgwyNfb/uolmLLpplzIKqKO0uR0wFcsU2UHXNgdeR0Ul8OgfzzT80D6FSOedBxcnL/GJ5b2Mhtt6Aej+kLN10MyrbRs42seC4uVTAghE/dDbDW7B3QHj+7ismKjaPROZaFk+GRr+RsIAK6iyeCNvAsKoGc3rupqg2DFW+cxTnHSsOVf4NW5D3gnH8RwFI/nqsbx2aq+L0feS1KlgHG4rmBZI9QQGx7oi+JH/LcQN7NI5+PqkroYFqR466EtbIe9QoHlGnhjtnRezmNaq1Q/qjlBTIApj1yILZDgM5bbUIdq1UtWR2KnH4uFLkxsDpy+f7v8UDOOcd9Ty0AAC4VDuTiPbnx8ATafpibdFSDXHpoQ8JaKaBu/SCUggWArFwpkvB87Owqbu2S6AQo2Zl8/dWmh8lKdiDvte8DyFbkUuhYhtz1XIT1tg8/5Dg2Lapu9vK8zx3zyBljdzHGHmSMPbiwsLDl53nlsSnc/UO34V+8/SQcK9omHES39MqOOPr3VsoPKeESH1RCgc/UHExUbHAel0Cp08LLTtGqlRYmK3Y8UchPTlOhLo0AULVF0y3HMnKnBNEaqo6JsaiXuYpMdg7YI1f9/r3Mc5fqeCmqqqC7tl7Qe3PTYWFV5FWudFPkbSWQF0p2BiGsjEDeK+HpBSGeurAu15pHSSm3BSA7f07mKfICa6Y7nKUMRV62TYyX7cKBnPx22smqFXkf4JzfzTm/jXN+2759+7b1XHdcuw//8l3XAkDUxTBS5F7skcvRabnWSn4pFGXM6WBSB9CS2kjvzCxbJspWd0VOP3eDEFMVG6VoolArpcjFmDhxAlLwr3YZE6daK1UnI9kZKZFqyULFsQYWyCmA73VrhdT4bM0pbK3Qe37joQkA+ZUravI9vbNT5kKKJjtDDttQPXJR191rU9DZ5Sa8gOPqfbWuj0sr8hVlvmcW3exMYjnDWqHXKNkGaiWzsLVC5++xyFrZy7s7R6pqJQu1C5yqyNuRSs+1VrqUQqnWCuccy3VX+Mu22RHIpSJ30uWHWcnO+HuT0eg3oWKSihyIE56kziu2mVt+qForYyWrI6mjeuQVO/+CsF3o7+hXn/dh5b6n5nH1vhpuPDyBS/XNeeRHpiqYqtq5ipzuLoF8j7ziFAvkwlqJj6m5WgmOaSSslW+8vNLRAfTF6G7jyrnugTztkdNYuFxrpdfAZi/u+Kne6ahiKev4zoMuCke1Ih9+1Ab7dOvZiqb0AEApr2ol+j71o1Aha8X1Q9TdAMsNDzNR6RMdpFR21lYOskSyM6Opj2MaoGqwyYotB1GkFTkQJzwrNv0/3xKhMW5Vx0LVMTv6kVNvaNNgqDjGAK2Vva/IG66P+59bwtuu24/ZmlPcWiGf1zEwN1ZKJMxVXL9L1YpyrMmg2OWi6Qc8Ya0YBsOhqbLc3fnCpTr+wW99CX/x4MuJ33txUdwtHO/RXCutyOlvykt29tq7QQraYMDSRra1UitZhXd2kiV6TCc7i8EY+wCArwC4jjF2hjH2o/143iLYytQRUoKJQL6FOnK11/nShotlpbk++X+xtZKT7My4gDDG5GOmKnY0rLnTIwdiJU4zP2keaBaUBK1Eijx9C1l3fdn5sWKb8EPet6EaKlKRbyOQz6+1ZPvVYeTpixtwgxC3XzWD2bFSYWtF9bdFiWh2UEns7MwL5I4RB8Uu77UbJMsPAXFHQB75M/PiruALTyVzVi8uNlCxTewbz+56SJQsM2mtROdEniK3e6yZAvnx2Zq8GwYUa8UyMF4unuwkj3z/RBm2yXSysxec8x/gnB/inNuc86Oc89/vx/MWQR2iLOvIlZK+3DryLr2RF9ZbODgh/MTFehvLDVcqctpl2WGt2EayjjznAkL2y1Q1UuR+tiKnypWqYq3kKWnayVqJPPK2HyaCYaMd94amTUu9RsdtBemRB1t77obr463/5T585Otn+7msvkLv9XjZliP4itTlq/ZA1bEypzgBan8U1uGRqxcDp0Apnx9w2TqWUHd3vnBJKO8vPXspcby8uFjH8dlq14oVALInOgVcUuS5gVz2WsmuWiFf/Op9Y2j7oTxG6b0rkXVYNJA3XFgGw3jJkq2h9yojb63YJuusI3fDrs2rgPzbPD8IsVh3cd1B0bpzqe5iue7K3spdk51Od0UOxIF0KrJWmm6g3DrGv0PWSpWsFSffI1e3bpOSV0sQE4o8WuMgtumrnRy3wkbLR8MNMsvj/vBLz+PrLy1va339QN0ARr0/FgvYK00vgMHE8VotoMjHy3bHhTuxIahI+WEYwjKSx+HhqQouron2Ey9EFsp6y8ejZ1flY15cbOQ2ylJJ17LTOZFbR96jvxEp6JMHxgDECl2tn6+VrMLJzuWGi+maA8YYak5xJT+K7IFAbshbNVWR08GSW7WSk+y8tOGCc8gezIt1F8sND9PRwVl1TFgG61DkaWslK9kJxLXkk1UnUuShcuuolB06KWvFzu+R0nADWFGJIil5VX003CCh7IHBTAlqp5LOm4VO2CwP9Fc/9RT+4sEzW19cn1CDykzUlXKpQMKTuhYyxlAtWV365gi1OlG2Oj7vli+mQNlmsaZZXsojB4AjU2WEHLi41sILi3VcNVcDY8DfPn0JgOhe+NJSo2eiE4iPcfrcV5ouTINJEZKm18WHFPk1+8YSX5NlV7JMjJcsbLjJcW+ffOx8pkpfqruyB/9mlPwoMvKBPJHsVBpQpTfZpMnrR06lhzccFGViC+ttrLXiuYWMscQ2fbXXim0aMqDm3ZZKa6Vio2wZaHvZijxtrVS7tJ+lqTHq76kHbcMNpCKnC8MgEp7bTXbKuYypQB6GHE0vkAOrdxNVFVNXyiI+eUv9jJzOWn+CRMlEJUuRhyhHwbOIteKlqlYAJNrZvnCpgVcem8ItRybxd6eFT35xvYW2HxZT5NHfQ+faSkNsBso79hljXYeeL9fFhYAuIpRIbqcUOeexNXh2pYl/ds/D+MRj5zOez8N0TQiwWjQHd68y8oE80yN3A+nT5ilyKyofSQcdKj28cq4GxzLw/KU6OIdU5AASgbztBdFOU/F8ZdvMrZShnwOKR674+VlVK2qys1v5IV0gaCiFetDW236HRz6IQE6BeKsbguj31tvJig5a60pOpcd2CEJeuPEVEHckpLawQLLmOff3vHikYNXppsijQF62Oyyqlh9/zkXb2NoZHjkAPH9pA+dWmzg+W8WbrpnDwy+tYL3l4cVFUXp4fLZAIDc7rZWpHH+cUIsT0izWXUxXbdlXny6QqiIfK4vzgmwSsmOy3s8lJbdV28SO0FFkDwTyTo+85cdVK3leNWMMTsZBRRUr+8dLmK05OB1l9qdr8SaHiYqNNcVaUevHy7aZm+gEYmtjqmrLBGY3RR6XH1q59d/qGvIVedJaGUQt+XZ7rdD7kLZW6CRdHkAg/9DDZ/DmX/l84SoetWZ/Vgac3h55yw/kcVErCY88q3kVrWOi0rlxS2117HRJ1svnCjurVg5HbS2++twSOAeumqvhzSf3IQg5vnR6UZYeXjlbwFqJ/h5Kwq42vdxdnYQqvNIs1duYrjpytmyWR04Ch4KyLAPO6DWu5ra0tTLkqENoVUWu7gbLw8q4zSNrZW6shJmaIzduTCuDZCeVQE7TgYiKY+TaOYDikVccpfwwo2ol7ZE7BhpekHnyN1z1tj0K5EqgFh55Mtk5SEW+VY+c3oe1jkAuvl5t9N9aubDawmrTK5z8VRuc1aKNOYU8cjc+TqqOsAey7lzoGB4v2R1rantxq+NivVaSbWwB8fnP1Bx86bTwxI/P1vCa49PYN17CPfe/iBcXG7AMhkOp6T5ZdHjkjfw+K0TenE9AWCEzNafjfVXryGUgj44Ruuin716CkItqs+i8zdrxvJcY+UCuqmo6qFW7Ik+RA9nqYH69jZmaA8cyMFNz5IEyU0sG8jjZGSYUubp9Ogt6bHJDEN06dipydWcntfNM00oocrJWxLo559hoex2KvEgg/83PPo2fuufhno8j+qfIk8p7kIrcDza35qYrqk8cU9hpczWnUOMs1RZJf0Yq9PmSIlcv3G0/vhgYBoNlsB7Jzk5FDoit+nTnedWssBD/yRuuxN8+cwmfffIijk5XMn8vTXqnZhFrJesumFisi3OPMSbaH8hkZwiDCTt0LHXHKQN56n24sNZCyOMOqWMlU1srw0wy2alUrSibCPKwzc7G+PNrbeyPNkLMKsFbLamaqFiJ8kP1Ncq22fU1y3YcyEu22FDR8gKx61NRTx3WSqS0s6pNGm4glbu0ViIVe2GthZYXJ682U0f+0IvLeOCFpZ6PI+Luh1tT+/T76WQnKXLVhuoXVCVSdAh3Kxp+TDmRmbFiuztbyig/ar+Q9Rmo5YecpyYGKc8BxHXceaTb2BJkr0xVbWmF/OPXHUfNMfH0xQ0cL2CrAPEdpCsVeWxl5NE12ansoJ6pOQlFTu/5eDm5s5rujNPH3JmozQC1sK1FXUH36vSqkQ/kpKo558nuhz3qyAHAMVnHLrOF9Zbc0UblZUCGtdLy5bxOdTJQuYcin605mBsTip/U8VrL67CAbr9yBm+7bh+unBMBuJuSVq0VUiyUGHzqwjoA4NoDopxyM3XkVENf5OAPwvj9H5RHDsQnbr8gRV7YI0993rO1UsFkZ6ymZUI6o5bcC4T6pAtyyw2Tz5Fqldxt3X4YdiQ7gTjhqfrgk1Ub33/7FQCKJToB1VoJoqSxj4ktJjvJCplVAjm9ry0/SBQJAHHNOh0r6c1TZ6JJSEejFra1koUg5Hu2M+eeCOSuUosNINd37vjdjBNhfr2N/ePCH6TyMscypOIFRCAPQo6Ntp8oKwOAd91wAO+64UDua971lhP43z/xBgBxcnOl4XWs88q5Gv7wn94u1Ru9ftbuNGGtkBVj4chUBaeiAP70xWQgl+WHBRT5Ut2FH52gvVCD99bLD+MLsfq5qBU4/bZX0u0detHywsTnXbQDohqEKxmVRYQblQzSsaHO7Wz5ybs/uqN8/OwqfvZ/P5KxuY13bAgC4slaV6Vqxd/7pqtQsU3c3KN9LSGtFS+UdliRqpWsjqMrDbF/g4oK1D42bS+UfzcFchIqeclOCuTU8TFtyew1siv3RwjHErdqiVtQt5giT3vkYcixsN7G/glS5OKgmq4ma2PV3Z0tL8TcWPw2/vgdJ7qud6Jsy23+dGKvNr3ELXMW3coGG66PivL7tx6dxKNnVgAAT13YwP7xkjxBSBX28sg551IRLdXdnkks9UTacrJTWdNGy5drbip92Ff6nPD0QrJWCgZypfoEEBf7Ijs7W0qiMmvTFkHjCbM2bqkljIDwm9t+iE8+fh7/+6Ez+J7XHMXrT8zGzxVkK3LyjdPK+8hUBV/9xXdgPGdDT5qS0pZ5pUfDLLnmjDmfQNypMLZWSrJxVsuP/27a1UoX9HUZyNOKvIEDE6WOYS31doDZsUJ/3kixJxS5F4Qy0TZRtqKdnWIXnGnk94tIe+RLDaFA94+nA3nS96OgdmnDxen5jUKbJ7JIKvLuHwUp6eV6pyJtuoFU7gBwy9FJvLjYwGrDw9MX12W7AUAkyUpW7w6IDeViWHTnIrFla0X5PTUxNVBFHr1mVvO0LFpu8g5splZCywt79vFQ79yqXRS5F4SwlUHeqiJv+0lrpRR55FT7fe+p+cRzpdvYEjQx58S+zog2WbETuZpuqNbKSo/t+YSj7MRWobsaOudmxxzUo/YVah6KMYbJqi03h+VVrZxZbkpbBYjtrL2a8NwjgTz2ZyerNvyQo94OulasAJFHrgRymjBOt56zOYGcfMB7T82j6QV488mtDcqoJBR5vgUEAK84NoWJsoU/+NLzHT9L+7a3HpkCAHzjzAqemV+Xtop83S5DKgg1eBcJ5O1EwNmqIo9/T92ko6rSfu/u9CNFXrTRVysVTIvu7lStlZqTr8g9n/dQ5MlkpxeEeHmpM5CHIUfI0bFFHwBuPjKB3/rBV+POmw52/2N7ENeRh9Kz7nXnZkd30J987Dx+8H9+VeZfOhV5XEve9sNEW+ipit3bWllpJOaNposA9hp7IpC7QdwkS/YLb3YmELN+NxHIo2ZNdOuZPqgIeo2PP3oOlsHw+qtnsRXo4FxpuD0V+WTFxk+89Wrce2oeDzwfV5J4QQgv4AkP/5YjwuP8xKPn0fJCXJcO5LbZMdczjRq8lwsF8n4o8nhNasJTPfn6vbszTtAWU+TNlCIvurtT2ANR1UpGYzNCKHImj131TkfdEATEVSsvLTVQtg2cnt/AS5E696I++1mKnDGGb7v1UFfbsQjqpiSyvCZzpgMRdM594Gsv48vPLsq/Tx2nCIh9HIDIWYlEcbzWqWocyLPKD/0gxPmVVmYg14p8SCFVTVdkGjO12vR6KnLRcCs+gUmRH50St2SzUdVK+naRPO7nFup49fHp3CZBvSC/uu4GPRU5APzTN1yFAxMl/OdPPimVjDr+i5is2jg+W8VfR/0nrj2YDOSHpypy5mQeaiAvWpVB9EORq4G86QZRLxvWd2tls1UraVVcZHenF4QIQq5UrVApaYZHHtkhcgeu+r4qPjsgjl9q6vZdrzoKALj31MXoNXn0mGI2yVZQe6306nxI2KaB9ZaP+59bBNC5Q5POLQrCZ5YbHYp8suJIK2ct1YUUAC6ut+GHPGGt7PVk58gHcts0wHnsN6qJyF6Kw04lXs6uNDFWsjBRER/6RMXC3FipI7uvbkO+4+Tcltee6F9eQB1VHBPve+e1ePilFXzptDgRyCJRNyUBwK1Hp+RJcnJ/0gu99sAYnrm43rWsMKHICyQY1R2P222aBQAbSr+VuuujVrIwWXH6bq14wSaTnSkbi0pVu02mj3vWJzdl5Xnkjql45NHvUnlnYi6saeDZqIXEHSfncGJfDZ+L7BW6QGVVrfSLhEfeoxc54VgGnrtUj/cMRMfoRsuHweL8QRzIm2hnKHLa5ZulyNM15ED3TVh7gdEP5NEHTAcE+dcrzd4JxLRHfma5iSNTFVmhwhjDvf/mLfgnb7gy8XtjjiVHtm3VHwfQ0aOlCN/1qiOoOqZU2y9Et9KzKfvn1sheOTpdkbeVxDX7x7Hc8LoqbQrk42WrUHkdqenxsr2tNrb0mamKnOrkp6u2TPbe/9wiTs+vb+l1VLxNK/JkID88Wcb+8RK+9kJ+r3S5xdyJd2VWHTPbI4828aT3DbjyQpm0VsieuWK2irdftx/3P7eEphvsiCK3DAbGRBBdbYrdw1lWjkr6Lpk2f220fYyVLHnujZdtTFVtnFludLznUxUbK00PnPPsQJ6qIQdURb43OyCOfiCPDgy60iYVeffgmOWRH1Gu4oC41UtvVzYMhomKONBuPlKs5jYLVV0VUeSAOJHfdv1+fOabFxCEHH/5jbOo2CbuuDZ5QbnlqFhX2h8HhCIH4hrzLBbrLmyT4fhsteDORXGCTJStnor8wmoL3/e7X5ETaoi2H0hvNG2t1Eqm8EYjRf7PP/B1/LfPne65rl7EdeQFPfLUvgHGGF5/Yhb3P7eYuMPhnOP7fucr+OtHz8dtWJXPuOpYmR6564uSwViRx20nAHQkO4ljM1Vce2AcbhBisd6Ws2iLbLXfKowxWTmz0vB67uoE4gsLvYfU6XKt5cldm8TR6QpeXmoKayWlyBtugPW2rzTMi9/LdA05EO+m1Yp8SKGDma7s6s6vntZKqo787HJDVqz04oqZKt5x/YGu5Y29SO8ILcqdNx3EpQ0XX3l2EX/96Hm8+6YDHar75iOTsE2Gmw5PdPz+yf0iuFNnxyyW667oRFcrYamAL02KaLxs9wzkv33faTzw/BKeOLeW+H7LCzFRseFYRqJqpe6KDU+TFQcrDQ8rDRfz6+2+nJR+qgVyL9IeOQC8/sQs5tfbeF65MK21fDzwwhIeeH4x0fSJEEOye3vkpMjV9rkEHd/TVbE3QU2QUu6nl0LeLjS3c7Xp9tzVqa7n7dfvB6Ao8pbfkWs6Nl3NVOST0QWDqnUYS+ZX0jXkgHivHNPAhq5aGU5oJuFGSpG7fti1LziQ7Jy43vKw1vI7FHkef/qjr8Mvf9fNW102gM175MTbrt8PxzTw7z72OFabHv7Bq450PGasZOEjP/nGzA1KByZKGC9bPRX5TM3BTNUupMhJEY2Xra7WyvxaC3/+tZcBdLYJaEc7F8dT47yaro+aI6yVlYYnhwa3ttjTRWUz1grnvKPUEwBef2IGgGgNS1Dy89KGq7RhTQbyvKoVxzKUoEzDUjp7B5FNcUW01Z4CV9sPlKqVwVkrtJ7Vpocnzq3J/RfdoIvPt94iSh+pImmj7cte48TR6QrOLDej91xR5NE5/vKSUN7TVafDWlFtFUIMl9CBfCjJs1aA7i1sgXhXKKCUHhZU5NS9cDtsVZGPlSy8+eQcnl2oY26shDddk51wvfnIZMftKiBuiU/uH8MzF7so8oaL2TGhyLM2IaWhQDNR6a7I7/7ic4mRfCrtSO2Ol61k+WFbNAUja4XWne6vsRU2k+xsZ/jUgNjqvm+8hK9GlRhAXOmzsNFODOgm8oYBU7JTDCtR5qBmKXIK5LIhWqzI6U5jkMlOQATmjz96DudXW/hnb7265+OvPziOW45M4vYrxcVP9cjHOwJ5NWoqFybUNd11kyLfN1ZKiIJ0DTlRK+UPvR519kwgTytyoHsLW/pdGciXkzXkOwGdrAAS5VVFePfNQtF8xysOb8kHPbl/XCrbLJaktWJjo+337GhIP58o27mPXa67uOf+l/DOqBdNtiI3MV62E61sm1Evmamqg5YX4rGzK5m/vxXSLZC70c5Q1kDsk39V8cljRd5ODMcmqk726DHPF8lOxhjKltmhyLM88uOpzpZtP+5Vk7UhqJ+ULGFPfverjiTaA+TxD197Bf7qn79J2jDr7S7Wykx8LiYVeWStLItAPjcuFDnnHEHIcX6llSnIqjlTtv7s/pfwoYd2fx7sduhLIGeM3ckYe4oxdpox9vP9eM6idA3km/DISZEfLajI+wEli4DNWSsA8K03H8S333qoo6KmKCcPjGGp7ubWPy/VRSe6ogOGSa1OlC14AUcYdiYPT11YR9ML8IOvO5b4HfU5SlaWIhfWCikxqhDpRyD3Za+V3snOrJp94vUnZhI+OfUov7TezvTIa06+IqdKrIoyp1U+h9XpkZMil+WAXigDeS8xs13KtomJsoVfeM8Nm/q9kiXm28oBETmKPH58pyKnvRCUIHeDEBttH37IOzbxAZRg7nzP/+QrL+B/3f/iptY/bGy7aRZjzATwWwDeBeAMgK8xxj7GOf/mdp+7CI6V7ZGLnxXwyBVF7piGPCh2CjFcIty0TTNetvHff/DVW37dk1E1y9MXN/Atqb/ZC0Q52XTNSWyVPjSZf5GjQEMnoxuEKBtm5mOmqw4M1tmBUU1qvVSPNyzRTFJSYpSkLWKtfOVZMb6MWrSmISVexFrJqhwhXneVUKMPvrCME/vGZMnmWsuXFyX196qlbHWoNroqW0ZH1UopQ5FfMdupyP1QfA6DVuT/5t3XoWQZsp6+KIwxjJXjOZrrLa9Dkav2iPreTWZYK4AQAnRxTCf/xfey3/NLG713Vg87/Vj97QBOc86f45y7AP4cwHf24XkL0eGRK5t1itaRc85xZqWJw1Plwg2D+gWpu50+kKgEMasWmzYAzaYCeTfafnIjS1bCUypax0TFNjOsFVLktgx+nHOxIcixEgOwxeN7K/IPPPAS/sunn879OZXpFbFWuinyq+ZqMA0mVaLaEZHu9pKj/LIHMFP3Q0DUnTc7kp3xc9hmtiJvKYp80B75267bjzdcvbVNcTVHBHIvED74WCn5+VYdS+6PUK3H8ZIF02B4ObJD6SLS9kL5nladzs+oYne+52HIsVRvY2G9PdJDJ/rxKR8B8LLy9ZnoewkYY3cxxh5kjD24sLDQh5cV0MFMJ37NsWBFwbiIIudc7Jo7u9xZQ74TUODbbuJ0sxycKGO8ZOHpjIQnJTfTirwb1DO623R31Ssu22ZH1Qkp8vGypTRDChFyEfzVi/SJuVohRd72AyzV2wgyrB5gc8nOLIuEMA2GgxNlnIuCtrqJ6kzk5SaqVnIqKDylY2HZMmVr36xk52uOT+Mt1+7DwYly4mdtP5DJzkFXrWyH8bKoTqL3IV21AsSqXBU6jDFMRkl1g8U9zNt+IJOZajdQQijy5Hu+0vQQRpOYVvs8tGQn6UcgzzpSOs4azvndnPPbOOe37du39d2QaVSP3Ira1tIB7Zg9NgRFB4cXcLEZaAf9cWKrHvl2YYzh2ExVBh4VUpMzmwjk1Ke721DgZjqQe/ke+UbbRxhyqaCERx77nrccnSzkkbvRhSBv/ZsZLJHeap/m0GQZ51bF+7mw0ZafKW1QUXfyVm0LbT+UW+nj9XB5MVQ98nZGsvMt1+7DH7/3dnkXmaXIB11Hvh3GSuJzJhGW9sgB4GgqkUtQCeJ4Oa4ea3mh9MBrGYpcDGBOHjNqjojmmI4i/fiUzwA4pnx9FMC5PjxvIRzFWqETgD7YIt0PAXERWFhv48jU1vqKb4fdUuRAcog0APzW50/jVz91SirymZqDqYoNg/XugNiOSsTSA3lVyBMvOybKdrInujo2b7wspszXXV8qqKpirVQdEyfmxuCHvCMQdqwrWselnKRuuo6cc577t2YFU5XDUxWcW2kBEAGC2gdTIFd3dlLvj0bqYuSqitxWPPIMRZ5G7c9CdxqD9si3A3nk5JNnDbQgRZ7+u+nubKJiJXq+yOMl47mqjtnRqEwdnD2/dnkH8q8BOMkYu4ox5gD4fgAf68PzFkLu7Gz78gOtOOL/vcsPxUH+3IKwF9Ryp51itzxyIGoHqgTyv3niAn7nC8/i4ZdEVchMzYFhMExXxfzEhfU2Hnl5JfO5SJGnB/ImHpNS5OpEID/qn00eOSA+U+l5lkw5Heaa/WPyM2718LbdHoGcLAjaCfmZb17E7f/ps3j87GrHY3sq8qkyLqy2EIZiuhIN9Di30oRlsESZqBzA3E5ezEQdebyNnS5+MtnZ5TgpKRdRv0sb22FhLNr4tdHVWkn6/4RU5CU78Xerd3Bpao6JhhckKqrUXMb8ems7f86usu1PmXPuA/hpAH8D4EkAf8E5f2K7z1sUW9nZKW9JyVop4JEDwGPRSXtdqt3rTkDqbhgU+XLDRciBP/2KKMWigRrTNQfnV1v4od+/H+/9o69lPldakbtRXa/aT6XpialNdpQUVa0VdbMN3WKvt/xE8ooxhkNTZdx0eLKjO2Ae3RQ551wpPxSPu7jWghdw/OrfPNXx+Kx6cJUjUxW4QYiL6y2sNDzRsMwx4Ye84/OV3fgUhRiEHJxDUeRZdeT5x4llipK+lqd65MMbyMlCoz0DWe2gbzkyCdNgOJyqmCKbTSjyuJ0uXRjT3UABodI5T+4IVnMZl7u1As75Jzjn13LOr+ac/3I/nrModKCqu7+Kqlz63UfOrMIyGK5JtXvdCcq7qMgnqzZWG57M1q/UPRhM3N5PlC35/szUHNx7ah6nLqxjse5mJgbFGDIl2RkE+NLpRbzt1++Tww6abjy4OG2tqOV1dEKvtzx5q1yxxffu+bHX4ee/9fp49miPARlSka932iVq7Xg7oDI/8f8vPr2ALz97KfH43h65CDaPnxU9ZGbHSpiLKirSdkyWIpcdC6P3UCR9qbtffBHsRskSczxdWbUyvNYKVa1088hfeWwKj//7d8sSS2JS8cjlpCI/UDzybGsFQKJyZXGjDYOJmHG5Wyu7iqq66d+lwoFcHOSPnlnBNfvHOibZ7wS76ZFPVRy4UemXF4RYb/v4nlcfhWmwxIaKmUj9UDI4a0qPrFox49vccytNcA4sbIhbVnVMWrr8kJSzaq2stXwZ6EjBHp2uYrKSPHm7QT9fyFDk6gWJ5kjS4w9MlPArn3oqUZJGQT5PkR+aFNUjdIc3V3PkvoQORe50KnK66FCwnijbsnqn5YUJjz0PUvGjoMjHyqIckIZDpMsPiSx1TZuCxsuWvKir1gpNYVLJunheinoKHZgoXd7Wym6jHqilTVorFHReXGzg+l2wVYBdVuSyd7srLZabj0zi/7zjBN5+/QH5uDdfO4d33rAfP/vu6wBkD5qg7fWqtUJBiE6ulhtIb7uUCuRqad9EpMw2Wr5MBqbrgtNtXvOIFXlnIPcVRe4pitxgwE++9Ro88vIKnl2IyzOzNuWo0IXusTMrACJFHs30TAdySsap5XCkoskjp741NIC4yMWeFHncxnZ4FTndeV1YEwE0S5HnQR55suujSHZaBsvMj1UzLp6LG23M1krYP14ubK187JFz+N7/8eXCa90Jtr2zc7dR62S36pEDwA2HOtu97gS76ZGrLX/pFnyqauNHUtv+/9HrjuMfve44vnxaWA1ZpXztaCZlSQ3k0cWBSr7UXt6ij4jikXudinw9mhoDdNYFZ41Cy4KUfqYiDxVFHgV1uiDRrfya0iqg5QVgLP+iO1W1UbYNPBZZK3NjqiJP/k4t4zY/XTI4oeQK1KEb3SBF7o2AIqfAfX61BcayN/HkIT3yspVIdtbbYhcwDahQybZWRHO46ZqDb6baKufxxLlVPPjiMsKQ7/gGwjyG91MuiHqg0lWYbsWK1pEDwPW7FsiHQJE3PGmXTHcZDkAbL7LK82gwsHpSURAk1akG8opjJBW5T2rXlNULay0vd6fephV5xpQj1VpxFUVeto140n07eddQtrKDBCBq8w9PVWRiVSjyKJBbOYq83RnISYBQY6m1lpewpbpRskV/8PiiMByBJguyUi6utRLTgYoQlx/aSvtesUU/yx8H4m376l3QYt3F7FgJ+8dLmF8rZq1QhdNWJ2ENgj0VyEspv7moRw4ANxzaHWulsoseuTpNiYYadxueKzcHZVornTs7SZFTMG66QfwZWSmPXFHkNUfYKy8vNRJ15Cpl5Xa6G92qVrKslXYUMOnCoU5dF10Yu39OVF1hmwwTZUtJdqYCud15m9+pyKNA3vSiwcvFrBXhke/MFv3tQBfs8yutzBrybsQbgqw4X+IFqLtBpj8OZM9KvbTRxmzNwf7xMupuUKhfuRvQbtvhCeQjb62Y0W7OIOSxIt+kRz43Jj7I3eA7XnEYVcfsGSAGgQzkDU/uz+2myCnIZylyYa0oHnmgeuTi5Gh5gbwlLtti1yLnHIyxxBZ0xhiuOziOpy+uY7JiC88z9VlKRR793pefvQTLMHD7VTPyMZyLgcWMCTsoCHliopOqqEi5k4UxlqHeiiQcabzYbK0Exhj25Xrknbf5bmqqDw0BX2v5siqoF2XbiBT58G/RH4veg/OrrU3v4TgxN4YT+2q4+chk4i6w6QY9FXnTiyuB1lt+dP6LC+78ehtX9biokCIXx2zvqUg7wfBerjcBHawl6TdHG4IKeuTXH9wdWwUArpyr4cfefGJXXlv1yFcLKPKSZWKsZGEpY9CEsFaSW/TXmuKEIQXU8kLFWjER8tibTk/Aue7gOE5dWEfDDTK9U7Iq6Pd+7W+ewn/9TLI5Fimm/eMlBCHvSNJmJzsjRS7rvONA2/QCOUA5DypBnI0CeJ5H7kQ136oCjK2VKNmpKHJ6f3tRivqz+GEYDUce5kAu/r6mF2TWkHdjsmrj3n/9Vtx0eBKOacgh0PW2nyuKZKVQdDxSrmd2rIT9E1EgL2CvbKaH/U6xRwJ5VAmRUuRF68h3y1bZbcaiLnIrTRfLDReWwXqeUNM1O6dqJRS9VjKqVqj+WrUmZF8QP9kUivzO6w6MY73l49mFjcwGSGlrZb3ld0wcIsVN1SRpe4VOyIptKsnOMLJ3Oof1tiOPvBv0WrNRAM8rP2SMdQw6cNPWiuqRF2x1TNv6vYAPdcUKkNzJmTXJqijU119s0Q8yd3UCce6M9h7QZiCyVoBim4JcHcgHg2w8ZaeSnT0COSVMbj06NbjFDTHURW616WGl6WGqavdUcDNVp6NqxQ9CBCFH2UpZK7JqJU52puvmW252U6jrorukb7y8kqnIS6mqlY2W3zmoInrOwxTIU5uCKJDXSqZirQgfv2KbYAyJAcnp2ZFZHIqslbkon5C3IUi8rpVU5Bl15ACw1vRlorUXYhiymBBkD7E/DiR3cmZtz98M4k5ENM3K6rMCxHkWykuoSWnVWulFnE8ZnkA+8h45EB/4so9zQUV+ZKqCT/yLN+9aDfkwMFkRA41DzhNDOfKYqjodipz6najdD5NVK0odubIhCIitkSxFDgilfXy2s5lZ2Y5fBxBJyXSSSyry6TxFLlR4xTFl0G/7IcbLFgyDoWonByS3vLB3snMqaa3UHBMzNQeztc7BC1XHTExCSivysm3ANhnWWp4s7+wFKXJ/FBS5EnA3m+xMQ4q82UWRi86ohjweqZJpbszBVNWGYxqFNgXRcaMVeZ+RgTwK3EemKnBMo2vijrjx8MTQ1ILuBqTIl+teofdrptapyKn5VSkqzXMsA003boaklh+m6+ZjayWpyCertuyzXbU7T3LyRVtRE6SNti+DcXpdR3OsFarsqDlWh0cOiBJBVTE33WLWyljJwlVzot0DYwwf/ak34q47OvMgJ/eP4/FzcXMuChAkQBhjYndn5JEX2xBkRhOCwqGuIQdEYKUL+mY98jQl2xCKvO1nWnFEVRmxt6gocsYY9o2XcPriBn77vtP48MP5Mzy1Ih8QMtkZnWRvvW4fvvILb5c+pSafyYrwvL2AF+rHPl11OqpW0kG4ZBqJZkQNV9zq+yFP9FoBkNHdLw5W1x4cx4W1VmY5mTqcmG6V09v1SeHOjZXgmAYWUrfN9POKE1srrh970WMlK6nI/d7JzrJt4gs/+9bE3c2xmez2yLdfNYNPPXEB51aaODxVyewhPlGxsdbyCyc7SZG70RDnYWesbIlkZx+slVbkkXfbWFR1TFm7v1h3ZbkrICYNfe7UPD53ah7jJQt//xWHM99D9VgZFob/ky5AWpEzxnQQL8hUNfLIG27XihVipmaj7gaJoJkOwo5lJHZS1tt+3DnQSXnkXqzIDZYsl7suGkeXd2JS4y1S/p2KPLZ85sacjt2dVLUiFDlVz8QBU5z0SvlhAUUOCIVnFQiiVCr5tReWACBzE89E2RKK3C+W7FQV+bBbK0BsqWwn2QmIu5j1lhi8nDWvk6gpA5gvbbQxF6lxAPixN1+Fu+44gV98z/VYb/t46MXlzOeQyc5g+8O/+8WeCOSyWdaID1DdDcgjX264HTMxs6DdnWrjLLXhFSA+D/IfTYOh4QYdnQNja0VJMqZ2TVLCM+9WuRK1wt2QHQKTgTzuXWJi33ipY3enpyryaHZrS/GiayUrsSGo5YeyV0w/uOHQBMZKFh54XgTydNMsQCjy1aYH1y+2IahsGwi5uNMZ5s6HBCnx7XrkZduUx2Q3RV5RKoVoez7x7bcexi++5wb8wO1XwDYZPv/UfOZzSGulwKjBnWJPRL60ItcUZyoKFC0vTIxSy4M6IS7VXXzysfP47t/+kvSR5Yg9y5B+9IHxEhpugJab7ByYtlayknmU8MxX5MJaWY9e3w3CxNCApCIvdTTO8kJS5OL5/ZCjHV1Q6PtqeWARj3wzmAbDa45Px4E8tUUfEJUr9F4WS3bGO1JHwlqJAvj2rRVD5m66BXIxt5OslbYc7qwyXrbx2itncN+p7NnCeov+gEh75JriTChebhFrRe238onHL+Dhl1bwpdOLABRFbsYn1cHJMhqun6vIyaJpKQGUOHlgDJbBZBlemlJKkQOp3ZrRra9jRoE8XbUSKWAqV/OCMKHIq6X4Nlyo9d5b9DfL7VfN4Jn5DSzV3Y7yQ0Ds7qSSuGLlh+J3Ry6Q96FqJQ7k3ZOdJDwurLZxYCJ7R/fbr9+Ppy6u42zGTFutyAeEVuRbR1XhRatWANFv5dGoXesnHz8PQOkDbxtyar0I5LG1QtZEunshbShSKdsm/uRHb8cPv+F45lrEdvQgaX9k9W+xxUDndB8NavVaUzaKUD08AIwpJ72wXvrfE+e1V8Y+eda2+omyLS2XXjNoxWMiRd7yR8Ij758ijwdV13J6rQDR3E4vgOuHWKy3cXAyO5C/9br9AIDPn+q0V0gstLUi7y+ODuRbZkpV5IXqyMVjnl+o48Vo8s+pC+sAkoqcODhRQcONh+KmFXlT2RCUpTjfcPVcbh8cqlpRA7nqk8ceuRGp6+S8RlfWkcftYoE4YFZLcYUDWUP9zsPcenQSjmXga88vZVsrymeyJxV5nzxy9SLXW5EHmF9vgfN4GEiaq/fVcGymkhnI9Rb9ASG36OtAvmkmq6q10luRk2r/wtPCP3znDfEAinSPG4OJki4AWI76s6Q9cpnsjIY3bwYqtVOtFfV2N1bkpmzQpG7jj+vIxc8okNNFhiocyFYBsqfVbIeybeIVRyfx0EvLcdWKsiNzQlGqxbboqx756CjyflStEF09csdEw/VxYVVs/MmzVhhjeO3xGTx5fq3jZ2rv+mFhT0Q+x9KKfKuoKny61vtksk1hUzz8kijNet87T8qfxeWH4v/jZVsG0KVoWrksP7RS1kqOIu8GJTuTijzo+LdjGvH27IwmVVUZyL3o74irVkIudnTSnUM/k53EzUcm8eT5NbQ80ehK3aCmBrgiyU5ae8MNhrqFLTFTc2AabFPTgbJQL3LdFbmwVs5FgfzQZP7eidkxR7Z3VhnGOvI9siFIK/KtMpmwVnorckCcfOstX7YRPbGvhucW6lJRk7UyUbHkSUWlf6TIjag1LSnklr/5DnhlW2wCSXrkiiJXvGV6bnWDDykrWuNaWpErk+4HpcgB4ObDk2h5IZ66sNZhh1ArW6BYMl8NaKOgyL//9ivwymNTXWu/i6Ce+1098pIFzoEXLtUBINcjB8QdajNjzJ47hDs7txX5GGP/B2PsCcZYyBi7rV+L2izU9lMH8s1DHmzJMgoHKbJXXhE1G3vH9SIxlO46OVG25Um1WKcSuvg1ypYR9zjxwk1XHZG1ovYrSSpyxSOXLUw7FTmtkS4IcflhrOKlIt+k/VOEW45OAgC+/vJKR/Cd2KIiB4Z7zBsxVrJwW5Tw3Q7qsdNNkZONJrpqmgnrKg0l9tXeQpzzofTIt6vIHwfw3QB+tw9r2TKxItflh5ulbJso20ahhlkEHeC3RgHop952DV59xbS0ARwlkNNJRaVhaiCvOPGUoK145KVeyU5loxIpcvWxfiAGTdBxI60VuSEo7l8dB/L+H2Mn5moo2wZWGp4c1kwkkp2b8MgBFNpdulco6pFTYvu5hToOTpS7dvucloNUPGnBBCEHj/LlwxTIt/VJc86f5Jw/1a/FbBVdfrg9pipOodJDgh5L7X+nqg6+9ZZD8udSkVcseVItpqwVIJ4SBGzdI297ITZaHuh8TCty0VyLZU6t9wLhSdPxI5Od0TroItRwfVyqU6e8/rd+sExDDv/usFa2o8hHYGdnv1BtvW53IqTIn1vY6GqrAHHyX1XknjKMRCc7+0y6ja1mc0xW7EKbgYgDEyU4loGbDmdPVspS5It1F6bBEtaBOrczq468F2XbgBuEWG3GnRvVqhXXjyfPj0n7JOmRixNfrKlTkccqnibHUN/qfnPzYXF300+PfBTqyPsFvTd58zqJqpIr6RXIs6yVrPGAw0BPa4Ux9lkABzN+9H7O+UeLvhBj7C4AdwHAFVdcUXiBRXBSo940m+PH7ziR28M58/FvPoE7bz6Ye6sfJztjj3yp7kbDGpRA7phxP/KC3f1USN1f2nAxG7XXbaWqVuiiItV1yiO3TAZbqb0GYkVeU+ZqLqy34Vibs6A2wy1HRCBP31VWbBOWweCHvNiGoBHzyPuFrDTq4o8DSdvlYE7pIZE1o9YLOpPpw0DPQM45f2c/XohzfjeAuwHgtttu4z0evim0It8e3/uao5t6/HTNkVv1s1AVOSVQlxtux3CFslK10i7Y3S/x+zbZNm3ccnQSz8znK/Jalkce9eym42YttSFITXbOr7exT+mU129uOpJtrTDGMFGxsVR3N+2RX06BnP7uXgl7NZDnbQYipqW1Epcguhk5mGFgT3zStWj25G5Motd04igeOQVDsb29cwt+2xPb4t0g3LQip+eru4FsW6yqpLbSMZDuONJT620znjPasSGopAbyltzcNAhO7h+HYxry7lKFKis2s7MTwEh0P+wXsSLv/h6piv1glxpyINozUbJSHvmIKvJuMMa+C8D/B2AfgL9mjH2Dc/7uvqxsE3zvbUdx4+GJrmVHmp2DAuN42U4kNyspRRkPQaDBFFtT5EA8I1NNQLlRshMQCcWSZSTKD4Uij5OdG5FHXlb6kQPiQrGw3saVs7VNrW8zOJaBGw6NZ4qRiYrd0as9D8s0pBVzWVWt2EkLLY/NKHIAmKrZudbKMCnybUU+zvlHAHykT2vZMhNlG68/Mbvby9BEUOJpIpp9WYmqU9JBir4fD6bYfPkhQYo8uSEoWdJYU7oZAojmWqrJTrJW4np40U9dWCs0CGJQ/Mr33oowIzZMlG2UU/mFbpRtExttP1Pd71VKqbxGHupg5rzt+SozVSdlrShVK0PUNEtLWE3fia0VkSyqlcxoXmdakSfrwLvV/2ahWjXTVRumwZKKPAgTeZNayUS9nfy56pHH5Yfia8YYao6JpbqHlYaX27yrX1x/MLsKaKJibeoiV7IMbLQvzzrySg9FTneFtskye5GnSQ8bJ0VuGUzOhB0GLp9PWrNjlMw42QnEt7ud1ooI5M9H26WPb9K6UC8MY2UR7NJNsxKK3Em2svWDpLWy3vJgGiwRAGslCy8uivUNqvSwF9fsH9/Ue0Pvy+VYftjLIzcNhrJt4MBEudDQ9ZladiCvlayhGiyhFbmm79x4eALXHxzHsRmRTCKlnZXsbPkhnpnfAACc3D+2qddJBPKSLQJ5qo2t2owpba14gUh2Uvlh3Q06AkGtZMm+HPsndieQv+8dJ/Ez7zjZ+4ERWe2E9zrlgh45IC7oRfxxQJQgrtQVayUK3mMlSw+W0Oxtbj4yiU+97w65ZZ8CeVay0/VFs6iZmrPpgdnqhWGsZMnBw0S6f0vVSVor8c5OpbY9tcaaY8pOefvGBmut5GEYDOYmKlDI47+8qlaKeeSAUNnHZqqFnne66mC97cvEJu3sHNOKXHO5QWV86WQnBc3Hzq7hmk2qcSBZjjdetmQVDOEGYWKDzVjJkn2oARHIayUroVzTXrSq8HZLkW8W+hsuK4/cJo+8dyD/3R96TeH+53LYeNPF/vGyDOhjZQsX1lrdfnVHuXw+ac2uEQ+T6KxaAYCnL65v2lZJP18tU5EHqWZKVqKO3A85LIOBsViVdyjy6CLEGAolx4YBulMZhTa2/WKibGO8ZOH4TO9cwol9Y4X3BKiNs4CkRz5MvVa0ItcMHKnIM6wVQHSU21ogT1krdqdHnlTkZmJnp+uHMtFpmwa8IOjYIk+36rO10sgoXLIZLqednRXHxAPvf2ff2wzPpBpnUSAfL1lw/RCc84Ht9t0Ml88nrdk18j3y+OtrD4xv+nnV3x/Pq1pRPfKSleh+6Ic8EcjTzynWLi5Cu1WxshUomI3KhadfVJzitfZFoQ6IK1Egl9ZKND3KD/vabWTLXF6ftGZXiKtW8gP5NQc2r8hLlgHGRFKvZBlyYhDRzvDIvYDLW2IvKj8E1ECePCWoa+Ko+ONA/L5eTm1sBwWNP1yS1ooI3HSXOSy7O3Ug1wwcUrXlnGTnZMXGvi30+WZMBPCxsiX/TYqcc55omiXWEfVbiSpXaGcnoHTQtEZfkV+Oyc5BMZ1jrYyVdSDXXGaQz9xhrUQB5+T+sS3fEpdtU07/UZOdVBqmKvJ0B0Ta2ak+rlORi98ZZMOsflNWdi9qtkfZNlGxTdlvRQby6JgelsZZOpBrBk6ly85OADi5BVtFPoelBvI42dlWxrwRNTnxhxR5p7XSocjJWhnw9vx+Qn/z5ZTsHCRid6ewVui40taK5rKDdktWnNTQhOj7J/dvPtFJlG1D7t4s2aY80dysQJ4askw7O4F8j7w2gtZK+TLcEDRIpqq2THZ6ys5OYHjGvelArhk40iNPKfITczW8941X4dtuPZT1a4WoOJbc3FGyjMToOPG9ZK05EM/tpAlBAOQ2/bQipykxh6a6964eJrRH3l+mqw6WlEAu+rUMl7Wi68g1A+f4bBWWwXAkFQwt08Av/f0bt/Xc/9e33SADtFpHTorcybBW6u04kDupZGdakb/pmjn8zx++Da84Ormtde4kFGQup14rg2S65uDsShMA3cUxebEclm36OpBrBs4NhybwxH98d6HhwZvlDdfMyX+XLFNu0qBb3ixrpd4WU4lCDlhG9zpyyzTwrhsP9H3dgyRW5Npa6QeTFQurTeGR07ASEgjD0jhLX7I1O8Iggnjna0Qnlx9mK3I5Qd2XXqdtscTjNjvcYhgp6aqVvqK2P/aivQl0PA+LIh/9o1ajiZC+pRdme+TSWgnkjjy7hyIfRW46PIHrDowXmoCj6Y3oqxLCD8JoE5kRi4YhGS6hrRXNniFW5EGmIi/bBgwmPHIv+jmpVsfcO4r8psOT+Jt/ecduL2PPoM5upf48zpB55KN/1Go0Eaq1kuWRM8ZQK1nYaPvwouGYliw/jHZ27gFFrukvarVTR7JzSKpWdCDX7BlKsiQskEmojm6GjiVPSCBW4nvJWtH0F5lbafvxnFdFNAwDOpBr9gy05b/lhfKWN22V0ABmn4bodtSR61NCk4Q2tNXbQWeycy8EcsbYrzHGTjHGHmWMfYQxNtWndWk0m6aQIo/mdsqqFVlHrhW5Jpt0tVNSkQ9HsnO78uMzAG7mnN8K4GkAv7D9JWk0W6Ok1Pa2g86qFSAuJSNrRSY7tSLX5KBWO3m+8Mjpwr8nFDnn/NOcc+rU/1UAR7e/JI1maySSnVFZWNbEH7pFBqD0Wske9abR0EayhutHU6dM2CYDY3skkKd4L4BP5v2QMXYXY+xBxtiDCwsLfXxZjUZQVqyVfI+crBWhyK1UslMrck0atf2xaOsg5rw6pjE0yc6edeSMsc8COJjxo/dzzj8aPeb9AHwA9+Q9D+f8bgB3A8Btt902HPORNHuKkpLslB55qt9I1bFSijw9IUgrck0SWX7YDhJzXtW2ybtNz0DOOX9nt58zxn4EwLcDeAfnXAdoza5RSilyxzRgpFq5TlQsrDU9eQJ2Jju1ItckoT76pMjjYSTm6ATybjDG7gTwcwDewjlv9GdJGs3WSHrkYYc/Doje524Q4vT8BoC4Z/f+iRIqtilb4mo0hGkwVGxT2RAUK/Jh8ci3u0X/vwMoAfhMNKrrq5zzn9j2qjSaLaD2WnGDINPvvv6gGGLx2JkVALEi//ZbD+Nbrp6VAwM0GhWRW4nu9Cya72oMzRb9bR21nPNr+rUQjWa7qL1W8hT5NfvHYDDgsbOrAOJAbhpspMa5aXYWUe2UtlaMoWmapQ1BzZ7BMhgMFu/szFLkZdvEiX1jeO5SHYBu9aopRo2S5Klk57Aoch3INXsGxhhKlom2H6DeDnJ7oF9/cByUltcDijVFiBU5l3d6QpHrQK7R9J2ybaDhBnjoxSXceHgi8zE3HIq/rwO5pgi0/8ANVEVuakWu0QyCkmXiay8sYbnh4W3X7898DCU8AT0OTVOMmmNhpSHGvTlKW4dhqVrRgVyzpyjZBp6+uAGDAXecnMt8zPVakWs2Sa1kYqXhAkjuPdgrTbM0mqGCEpyvOT6NqaqT+ZjDk2WMl0XBlk52aopQdSystURbKWmt2MUU+VrLw0MvLmOQ+yV10axmT0G15Hm2CiCSojccnMADLyxpRa4phLq/gHrX5/Vauf+5Rfyrv3gEc2MOqo6FB19cghdw/OVPvRGvPDY1kPXpo1izpyBF/rbr8gM5AFx/aByMxTs7NZpuVEtxBZQjxwJ2KnLOOf7TJ55E2w8xUbGx3vbwnlsOAQAurDYHtj6tyDV7iopj4dBkOZHQzOLH3nQCrzg6hWhHskbTFepJDqgeudkRyD/75DweObOKX/2eW/F9rz0GADi30sRHv3EOy1GydBDoQK7ZU/zs37sOLT/oGaCvmK3iitnqDq1KM+rUFGtFDiGxk9ZKGHL8+qefwpWzVXz3q4/I709HuZrlKFk6CHQg1+wpbjk6udtL0OxBaG4nkKxacYMQnHMwxnDf0/M4dWEdv/kPXyn73ANAxTFRsgxZvjgItEeu0Wg0PUgo8ihIVx2aHCRKEJ+5KDpqvuOGzvzMdNXBcn1wilwHco1Go+lBrdSpyGdqwjJZigL0Ut2FYxmZHTSnqvZAPXIdyDUajaYHqiKnvQdzYyUAwKWNdvR/F3M1JzM/M1115IaiQaADuUaj0fQgUbUSJTtnx4QiX9wQAXqx3sbMWPYmtOmaPdBkpw7kGo1G04Msj3w2UuSLdaHIl+ouZmulzN+fqjo62anRaDS7STWjamU28sgvkSLfcKVKTzNdtbHS9Aa2TV8Hco1Go+lByTJgGnHXQ0C0gxgvWbi00QbnHJc22jK4p5muOghCLvu19BsdyDUajaYHjDFZS642Wpsdc7C44aLhBmj7obRb0lADt0ElPHUg12g0mgKQT+4om31mx0pYrLdlwjNPkU9VbAAYWAmiDuQajUZTAArkasfM2ZpQ5JeihGeuR14TgVwrco1Go9lFpLViJRX5pQ0XS1KR97JWhlCRM8b+b8bYo4yxbzDGPs0YO9yvhWk0Gs0wESvy2COfG3OwVG9jYaOHIh9w46ztKvJf45zfyjl/JYCPA/il7S9Jo9Foho9qtCnINpLWSsiBZ+c3oq+zFflkxQZjQ+qRc87XlC9rAAY3y0ij0Wh2kbGSCctgMAy1akUE7qcurqPqmKgo9eYqpsEwUbYH5pFvu40tY+yXAfwwgFUAb+vyuLsA3AUAV1xxxXZfVqPRaHaUasnqGA1I/Vaevriea6sQ0wNsnNVTkTPGPssYezzjv+8EAM75+znnxwDcA+Cn856Hc3435/w2zvlt+/bt699foNFoNDvAHSfn8G23Hkp8by4K3hfX2pjJsVWIqQE2zuqpyDnn7yz4XH8G4K8B/LttrUij0WiGkDtvPoQ7b04GcnUD0FxODTkxXbVlUrTfbLdq5aTy5XcAOLW95Wg0Gs3oMFWxQZZ5b2vFwXJ9MNbKdj3y/8wYuw5ACOBFAD+x/SVpNBrNaGAYDDO1kuizkrM9n9hVa6UbnPPv6ddCNBqNZhSZG3O6Nswipqs26m4A1w9l461+oXd2ajQazTYgS6WXtTJVG1zjLB3INRqNZhtQCWLeZiBiujq4xlk6kGs0Gs02oAA+09NaGdw2fR3INRqNZhsUtVYOTJTwiqOTsIzO4czbZds7OzUajeZy5u/fehheEOLgRLnr467ZP46P/vSbBrIGHcg1Go1mG1wxW8X73nntrq5BWysajUYz4uhArtFoNCOODuQajUYz4uhArtFoNCOODuQajUYz4uhArtFoNCOODuQajUYz4uhArtFoNCMO43zn5yUzxhYg+pdvhTkAl/q4nEGg19gf9Bq3z7CvD9Br3AzHOecdszJ3JZBvB8bYg5zz23Z7Hd3Qa+wPeo3bZ9jXB+g19gNtrWg0Gs2IowO5RqPRjDijGMjv3u0FFECvsT/oNW6fYV8foNe4bUbOI9doNBpNklFU5BqNRqNR0IFco9FoRpyRCuSMsTsZY08xxk4zxn5+CNZzjDH2ecbYk4yxJxhjPxN9f4Yx9hnG2DPR/6eHYK0mY+zrjLGPD+MaGWNTjLEPMsZORe/ntwzhGv9l9Dk/zhj7AGOsvNtrZIz9AWNsnjH2uPK93DUxxn4hOn+eYoy9exfX+GvRZ/0oY+wjjLGpYVuj8rN/wxjjjLG53VxjN0YmkDPGTAC/BeBbAdwI4AcYYzfu7qrgA/jXnPMbALwewE9Fa/p5AJ/jnJ8E8Lno693mZwA8qXw9bGv8fwF8inN+PYBXQKx1aNbIGDsC4F8AuI1zfjMAE8D3D8Ea/wjAnanvZa4pOja/H8BN0e/8dnRe7cYaPwPgZs75rQCeBvALQ7hGMMaOAXgXgJeU7+3WGnMZmUAO4HYApznnz3HOXQB/DuA7d3NBnPPznPOHo3+vQwSfI9G6/jh62B8D+Ae7ssAIxthRAN8G4PeUbw/NGhljEwDuAPD7AMA5dznnKxiiNUZYACqMMQtAFcA57PIaOedfBLCU+nbemr4TwJ9zztuc8+cBnIY4r3Z8jZzzT3PO/ejLrwI4OmxrjPivAP4tALUqZFfW2I1RCuRHALysfH0m+t5QwBi7EsCrANwP4ADn/Dwggj2A/bu4NAD4TYiDMVS+N0xrPAFgAcAfRvbP7zHGasO0Rs75WQD/BUKZnQewyjn/9DCtUSFvTcN6Dr0XwCejfw/NGhlj3wHgLOf8kdSPhmaNxCgFcpbxvaGonWSMjQH4EID3cc7Xdns9Koyxbwcwzzl/aLfX0gULwKsB/A/O+asA1LH7Vk+CyGf+TgBXATgMoMYY+8e7u6pNM3TnEGPs/RAW5T30rYyH7fgaGWNVAO8H8EtZP8743q6+j6MUyM8AOKZ8fRTi1nZXYYzZEEH8Hs75h6NvX2SMHYp+fgjA/G6tD8AbAXwHY+wFCDvq7Yyx/4XhWuMZAGc45/dHX38QIrAP0xrfCeB5zvkC59wD8GEAbxiyNRJ5axqqc4gx9iMAvh3AP+LxhpZhWePVEBftR6Jz5yiAhxljBzE8a5SMUiD/GoCTjLGrGGMORLLhY7u5IMYYg/B1n+Sc/4byo48B+JHo3z8C4KM7vTaCc/4LnPOjnPMrId6zeznn/xjDtcYLAF5mjF0XfesdAL6JIVojhKXyesZYNfrc3wGRExmmNRJ5a/oYgO9njJUYY1cBOAnggV1YHxhjdwL4OQDfwTlvKD8aijVyzh/jnO/nnF8ZnTtnALw6OlaHYo0JOOcj8x+A90BkuJ8F8P4hWM+bIG6pHgXwjei/9wCYhagWeCb6/8xurzVa71sBfDz691CtEcArATwYvZd/CWB6CNf4HwCcAvA4gD8FUNrtNQL4AIRn70EEmx/ttiYIu+BZAE8B+NZdXONpCJ+ZzpvfGbY1pn7+AoC53Vxjt//0Fn2NRqMZcUbJWtFoNBpNBjqQazQazYijA7lGo9GMODqQazQazYijA7lGo9GMODqQazQazYijA7lGo9GMOP8/K4rJ+D/LOyEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data = hal.get_data_linker(\n", " {ar_graph.a: np.ones(p)/p,\n", " ar_graph.c: 0.},\n", " n_data=150\n", ")\n", "m = hal.get_generative_model(ar_graph, data)\n", "plt.plot(m.get_example(ar_graph.y))" ] }, { "cell_type": "markdown", "id": "6337ad88", "metadata": {}, "source": [ "## the ARMA model\n", "The ARMA model combines the autoregressive and moving average concept\n", "\n", "$$\n", "y_t = c + \\epsilon_t + \\sum_{i=1}^p a_i y_{t-i} + \\sum_{j=1}^q b_j \\epsilon_{t-j}\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "id": "394a4868", "metadata": {}, "outputs": [], "source": [ "p = 7\n", "q = 5\n", "sigma = 1.\n", "sigma_obs = 1.e-2\n", "\n", "arma_graph = hal.Graph(\"graph\")\n", "with arma_graph:\n", " hal.StaticVariable(\"a\", shape=(p,), mean=0, variance=1)\n", " hal.StaticVariable(\"b\", shape=(q,), mean=0, variance=1)\n", " hal.StaticVariable(\"c\", shape=(), mean=0, variance=1)\n", " \n", " hal.Variable(\"epsilon\", mean=0, variance=sigma**2)\n", " \n", " hal.Variable(\"y\")\n", " \n", " ar_contrib = sum([\n", " a[i] * y[hal.TimeIndex-i-1]\n", " for i in range(p)\n", " ])\n", " \n", " ma_contrib = sum([\n", " b[i] * epsilon[hal.TimeIndex-i-1]\n", " for i in range(q)\n", " ])\n", " \n", " y.mean = c + ar_contrib + ma_contrib\n", " y.variance = (sigma/100)**2\n", " # OR\n", " hal.Variable(\"observed\", mean=y, variance=sigma_obs**2)\n", " # this \"observed\" variable is strictly speaking not part of the model,\n", " # but it is recommended when we want to fit data to the model as it allows us\n", " # to include observation noise." ] }, { "cell_type": "code", "execution_count": 8, "id": "c659cf29", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6MElEQVR4nO3dd3yb13X4/8/FJgFwAtxTe0uWZHmvWN7xTJw6iZu0qZM2cX/fZjZJna6MNs1o6jTTTdomTVLHTbxix7Zs2Y48ZVF7ixSHOMS9Jwjg/v4AQJESNwFinffrxZfIB8DzHIni4cV57j1Xaa0RQgiRmAzRDkAIIUTkSJIXQogEJkleCCESmCR5IYRIYJLkhRAigZmiHcB4LpdLl5WVRTsMIYSIK3v37m3XWrsneyymknxZWRkVFRXRDkMIIeKKUqpuqsekXCOEEAlMkrwQQiQwSfJCCJHAJMkLIUQCkyQvhBAJTJK8EEIkMEnyQgiRwCTJxzi/X/NYRT2tfcPRDkUIEYckyce4f3nhBH/9m0M8+k59tEMRQsQhSfIx7L/fqOHHf6gGoKZ9IMrRCCHiUUy1NRDnnGrp48vPHOOGNbn0DY9SLUleCDEPMpKPIaM+/9jn//5yFSlmI994zwaW5Tioaetnrls1fvLR/Xzu/w6GO0whRByRJB8jqlr72PSPO/jey5VUtfbzzKEmPnR5GZl2C+UuB73DXroGR2d9Pq01O0+0sr++O3JBCyFinpRrYsRv9jYy4PHxrR2n+HVFPTaTkQeuLAdgicsOQE17P1n2rFmdr7ZjkL5hL2iZlSNEMpORfAzw+zVPH2jkmhVu7tpUQH3nEH98WSnZDisAZcEkX902+7r8oYZuAPpGvAyMeMMesxAiPshIPgbsqe2kqWeYz9+yindvKODGtXlctzJn7PGizBRMBjWnGTaHGnrGPm/uHWap2xHWmIUQ8UFG8jHgyQNNpFqM3LAmF6NBcev6fFIsxrHHzUYDJVmp1HbMbSRvMQa+vc09UrIRIllJko8yj9fP7w+f5cY1uaRapn5jVe6yz7pc4/NrjjT2cvmybECSvBDJTJJ8lO2p7aRnaJTbNxZM+7wyl53ajgH8/pmnUVa19jM06uPGNXlAoFwjhEhOEa/JK6VqgT7AB3i11lsjfc14EqqdbynNnPZ55S47w6N+mnuHKchImfa5B4M3XbeVZ5FmM9EiSV6IpLVYN16v01q3L9K14srhxm6Ks1LISLVM+7zQNMra9oEZk/zhhh6cVhNLXHby0m1SrhEiiUm5JsoON/awoTBjxueVu4PTKGcxw2bfmS7WFaZjMChy02wykhciiS1GktfADqXUXqXUx85/UCn1MaVUhVKqoq2tbRHCiR3dgx7qO4dYV5g+43NznTYsRgP1nYPTPq+pe4ijTb1cvcINQH66TWryQiSxxUjyV2itNwO3AA8qpa4e/6DW+hGt9Vat9Va3270I4cSOw42Bevz6WSR5g0GRl27j7Ayll5eOtwBw49pcAPLSbLT1jeAd1xdHCJE8Ip7ktdZNwT9bgSeAbZG+ZryYS5KHwKj8bM/QtM/ZcbSFpW772OKn3HQbfg1t/SPzjnPQ4+XFYy3zfr0QInoimuSVUnallDP0OXAjcCSS14wnhxt6KMlKJT3VPKvnF2Sk0NQ99Ui+Z3CUt6s7uHFt3tixvDQbsLC58v/y3Ak++vMKTrf1z/scQojoiPRIPhd4XSl1EHgHeFZr/XyErxnT+oZH+fWeM/QNj3K4sYf1RbMbxQPkpQduovqCc+WrWvs51NDN8KgPgFdOtuL1a25ckzv2mtxgkp/vzdea9gF+ufsMAKdbJckLEW8iOoVSa10NbIzkNeLJ8KiPB35Wwe6aTr694xStfSPcf2nprF9fkG7D69d09I/gdlq590dv0jU4ikGdm0ef47SysShj7DX56QsbyX/zhROYjAqvX8+prYIQIjbIFMpF4vNrPvXrA+yu6eQzN6wY6zA50yKo8fLTA/Pjm3qGaekdoWtwlPdtLeIv37WcJW4HJqPiQ5eVYjCosddk2S1YjAaae+dekz9Q383vDzfz8WuWkWW3yBaEQsQh6UK5CI409vDQE4c52NDDl25bzQNXLeHj1y7lRHPfrKZPhuQFR+Vnu4foHw60D75rUyGXL3NN+RqlFDlp1nmVa5452ITFZOCBq8r5w6lWSfJCxCFJ8hH2ZlU79/90N1l2Kw/ft4k7NxUCYDIa5pTggbGVrmd7hmntC4zMl+XM3EI4L802rwT9VnUHW0oysVtNlLscvFEli5aFiDdSromw5440k2oxsfMz14wl+PnKTDVjNRk42zNEVWs/TqsJt9M64+u2r8nlQH03T+5vnPW1egZHOXa2l0uXBDpZlrtSae4dZtAjG5AIEU8kyUfYntpOLirJID1ldtMkp6OUIj/dRlPPMKfb+lma40ApNePrPnrVEi4uy+Rvnzwy44rZkHdqO9EaLl0S2G6w3BV4x1DbPrvXCyFigyT5COoZGuVkSx8Xl81uX9bZyE9P4Wx3YCQ/m1INgNGg+M4fbQLgS0/ObpnCW6c7sJoMbCrJAKDMlQogdXkh4owk+Qjad6YLrWFr2exn0MwkP8NGVWs/rX0jc9rSrygzlfu2FfPW6Q483plbHLxd3cGW0kyspsAOVWXZwS6YMo1SiLgiST6CKmo7MRkUFxWHL8kXpKfQG5xZM9uRfMjG4gw8Pj8nm/umfV73oIfjzefq8QB2q4ncNKuM5IWIM5LkI2hPTaDl7/j9WhcqNI0SYGmw/fBshRZJhTYVmcrumlA9PnvC8bJsuyR5IeKMJPkIGfH6ONDQzcVhLNUAFGQEkrwluLn3XBRlppBlt3Cwvnva571W2UaqxcjG4olTPJe47dRKkhcirsg8+TDy+zW/3dfArsp2slLNeLx+tobxpiucW/Va5krFZJzb72ilFBuK0se2HJyM1pqXj7dy5TLXWD0+pCzbTseAh0d2nSbLbuWeiwonrK4VQsQeSfLz1No7TFv/CGsLAqPdxu4hPvXoAd6p7STbbqFjwIPFaGDrHNoWzEaoF81cbrqOt6Eog12nKhn0eEm1XPjtP9HcR1PPMH+1ffkFj11UkolBwT/9/kQwBjsXlYT37yeECC9J8vP0ud8cYl9dF3u+tB2b2cg3nz/BkaYevvGeDbx3SxFdgx76R7xjPWrCJT3FzBK3/YJ6+WxtLErHr+FIYy/byi98l/HyiVYArluZc8Fj28qzOPqPN7PvTBcf/MluugY984pBCLF4JMnPQ2P3ELsq29AaXjnRynWrcnjxWAt3bCzgfRcXA5DtsIY9wUOg5PLyZ66d9+s3hG6+1ndPmuR3Hm9hY1E6OWm2Cx4DSLEYx95N9A7J6lchYp0k+Xn4TUUDAE6bid8dakIpGPD4ePeGgihHNjO300phRsqkM2w6+kfYX9/NJ69fMe050oKrd3uHRyMRohAijCTJz5Hfr3msop4rl7kod9n59Z56BkZ8ZNstYy0AYt3G4nT21Hbi82uM426c7jzRitZw/eoLSzXjpdkCSb5nUJK8ELFOplDO0Run22nsHuJ9W4u5Y2MBI14/fzjVxi3r8+Y82yVabt9QQEvvyFj9HQL7uD78UiUrch2sLUib9vUWk4EUs1FG8kLEgfjISjHC6/Pz7zuryEg1c+PaXDaXZFIYbP97exyUakJuWJNLfrqNn79VO3bs4Z2VNHYP8bW718+q6Vlaiklq8kLEAUnyc/D1507wTm0nf3/7GqwmIwaD4gOXlLAqzxnWJmSRZjIa+OAlJbxW2c7ptn4O1nfz09dq+KOtxbP+e6SnmOkZkpG8ELFOavKz9LuDTfzk9Rr+5PIy7r6oaOz4g9ct48HrlkUxsvm5b1sJ391ZxZ//z15q2gfIslv4wi2rZv36NJtZyjVCxAEZyc/CiNfHV589xsaidB66bXW0wwkLl8PKHZsKqGkf4L6Li/n9/7uKTLtl1q9PS5EkL0Q8kJH8LDy1v4mW3hG+8d6NmOPk5upsfPWudfz1TSunnBM/nfQUM5Wt03ezFEJEnyT5Gfj8mh/tOs3agjSuXj71htnxyGY2YjPPr0Nmmk1uvAoRDxJnWBohLx5rprptgI9fu3RWs06SRVqKmb7hUfx+He1QhBDTkCQ/g0d2VVOancot6/KjHUpMSU8x49fQLxt7CxHTIp7klVI3K6VOKqWqlFJfiPT1wuloUw/7znTzocvKJqwMFedWvfZGeRrl8KiPlt7hqMYgRCyLaJJXShmB7wO3AGuA9yul1kTiWpEoG/xq9xmsJgPv2VwY9nPHu7SUwO2caNflv/rsMW59+DUpGwkxhUiP5LcBVVrraq21B3gUuDPcFznW1MsN3/kDp1rCN9ujf8TLk/sbefeGAjJSZz+1MFmEmpRFc0HU8KiPpw400THgoa5zMGpxCBHLIp3kC4H6cV83BI+NUUp9TClVoZSqaGtrm9dFctKsdA+O8slHD+Dx+ucf7ThP7m9kwOPj/ktLwnK+RDNWrlnkufJDHh/dwT72r55spS+4qfnRpql3uxIimUU6yU9WyJ7wvlpr/YjWeqvWeqvb7Z7XRVwOK/98z3qOne3l4Z2n5nWO8z19oIlVeU42FWeE5XyJJj0lOjX5T/36ANd+61Uau4d4Yn8jLocFs1FxtKl3UeMQIl5EOsk3AMXjvi4CmiJxoRvX5nHvliJ++OppjjQufFRX3d7PpuIMmTY5hWiUa5p7htlxrJnuwVE+8Yu9vHKijds3FrAsxylJXogpRDrJ7wGWK6XKlVIW4D7g6Uhd7G9vX4PFZODRPWcWdJ6BES/t/R6Ks1LDFFnicVpNKAW9w4t34/U3e+vxa/jcTSs52NCDx+fn7osKWVuQxrGmHrSWm69CnC+iSV5r7QX+EngBOA48prU+GqnrpdnM3LAmj2cPnWXUN//afEPXEIAk+WkYDAqn1bRo5Rq/X/PrinouX5odbAq3lGtXullfmM7agjTa+z209o0sSixCxJOItzXQWv8e+H2krxNy58YCfnewidcr27lu1fQ7HE3lTHCmRokk+WmlpZgXLcm/Vd1BfecQn71xJQCfu+lcx8y1BelA4OZr7jz68AiRyBJuxevVK9ykp5h56kDjvM8RSvLFmSnhCishLWa74cf3NZKeYuamtXkXPLY63wnA0UapywtxvoRL8haTgVvX57PjWAuDwSX3Xp+fn7xWTUPX7OZS13cOYrcYyZpD691klJ5iXrTFUCdbetlUnDFpQzWnzUxZdqrcfBViEgmX5AHu3FTAoMfHv+44hd+v+dunjvLVZ4/z7R2zm15Z3zlIcVaqzKyZQVqKaVFm12itqWsfpDR76vLZ2oJ0jshceSEukJBJ/pLyLN6/rZifvF7DLQ+/xv++c4Ycp5XnjzQzMDJx5Onx+i9YEl/fNSg3XWdhsco1XYOj9I14Kc22T/mcrWWZNHQNyaIoIc6TkEleKcU/3b2eL96yilOtfdx9USHf/+BmhkZ9PH+keex5WmtufngXX3j80IRj9Z1DctN1FtIX6cZrXccAAKXTfE/uuaiIFLORn71ZG/F4hIgnCZnkIZDo//yapbzx+Xfx7Xs3srU0k5KsVB7f3zD2nNqOQarbBnisooHXK9sBaO/3MDTqk5uus5CWYmbA41vQdNXZCN0In65ck55q5u7NhTx5oInOAU9E4xEiniRskg8pyEjBYFAopbj7okLePN3B2Z7APPjd1R0AuBwWvvTkYYZHfeemT06TUERAmi0wA7cvwgui6jqCs51meHf1J5eX4fH6F7wYTohEkvBJfry7LypEa3hif2B65e6aTlwOCw/fdxG1HYN87+Uq6semT0qSn0l66uz71zy+r4Ha9oF5XaeuY5C8NNuMWxWuyHVyxbJsfv5mnYzmhQhKqiRf5rKzrSyLx/bUo7Vmd3UHl5Rnc8UyF/dsLuTHu07z8olWAIokyc8o1Ilyphk2PUOjfPqxg9z5/TfYXd2Bx+vnQH03fbO8aXumc2DW76w+fcMKugY9fOA/3qa9X1bACpFUSR7gvm3F1HYM8n97G2jqGeaSJVkAfOm2NTisJp4+2ITbaSXFMr8NrpNJRnAk3zU4/ai5uSewc9OI18f9P93Npi/v4K7vv8F3d1bO6jq1HYPT3nQdb0tpFv/5JxdT2zHA/T/ZHfH7BULEuqRL8reuzyfNZuJrzx4H4JLybACy7Bb+5tbVgLQzmC23I9BCoG2GnjHNwe35vv+Bzbx3SxH3bC6kNDuV42dn3uRl0OOlrW+EMtfU0yfPd8UyF1++cx0nmvs4JgukRJKLeO+aWGMzG7n7okJ+9lYdmalmluc4xh5775YiXq9qZ3V+WhQjjB85aVaAGRuDNQdvdK/IdXL96lwAPv3rA7wVvPE9nfn2EbpquQuAvXVdbJQ9AUQSS7qRPMAfXRzY7WlbeRaGcRt0K6V4+L6L+ItrlkYrtLhiMxtx2kwzj+R7Ao+Pbx62NMfB2Z5h+kemn5kTmlkz3fTJyeSnp1CQbmPfma45vU6IRJOUSX5NQRqfuWEFH71qSbRDiXtup3UW5ZohXA4LFtO5/25L3YF3UNVt/dO+9kwoyWfNvlwTsrk0k311kuRFckvKJA/w/12/nK1lWdEOI+7lOK209g1P+5zmnmHy0ie2AF6WE0jaVa3TJ/m6zgHSU8xj0zXnYktpJk09w2PrIoRIRkmb5EV45DhtM9bkz/YMk3den/fSbDsmg+L0NCP5F4+18PyRZpa65z6KB9hckgnAvrrueb1eiEQgSV4sSI7TSmvvyLRb77X0XjiSNxsNlGSnTjmS/8bzJ/jozytwOax89a7184ptTUEaNrOBvVKyEUks6WbXiPByO60MjfoY8PhwWC/87zQ86qNrcPSCkTzAMreD020XroL1+zW/3H2G7atz+MEHt0yo5c+F2WhgQ2GG3HwVSU1G8mJBxqZR9k5el28JHs9Lv7Dh29IcB7XtAxcsWKps7adnaJSb1+XPO8GHbC7N5GhTD8OjvgWdR4h4JUleLEiOMzBCn6oufza42nWqkbzXr8fmwoe8U9sJwLYw3BjfWJTOqE9T2TL9DV4hEpUkebEgOc7pF0SdG8lfmOSXBheinT6vLr+nppPcNCvFWQtv97wiL7D/68mWmVfXCpGIJMmLBXEHk/xUc+XHRvKTJfngrJmqcTNstNa8U9PJxWVZYdl+sTQrFYvRQKUkeZGkJMmLBUlPMWMxGaacK9/cM4zTapr0pqzTZqYwI4XDDee27GvoGqK5d5ht5eFZw2AyGlia45CRvEhakuTFgiilcDustPVOPpJv7hkmd5JRfMhVy128XtmOxxu4+fpOTaAef3EYF6qtyHVITV4kLUnyYsFy0qxT33jtHSZ/miT/rlU59I14qQjebN1T20mazcTKXGfY4luR66Sxe2jW/euFSCQRS/JKqX9QSjUqpQ4EP26N1LVEdLkd5/rXeLx+Xq9s5zsvnmLfmS5aeoYnNCY73xXLXFhMBnaeaGXI42PniVa2lWdPaBy3UCuCvzBOyWheJKFIL4b6jtb6WxG+hoiynDQre2o7OdMxyD0/fHNsR6aHg5uCTDeSt1tNXLYkm5dPtJJlt9DWN8KfXxPexnGhdwWVLX1sKc0M67mFiHVSrhELluO00TU4ypefOcbAiJcf//EWKr60nc/dtJLCjJQZG8FdvzqHmvYB/v3lSravzg1rPR6gKDOFFLNRbr6KpBTpkfxfKqU+BFQAn9FaX7C+XCn1MeBjACUlJREOR0RCaK78S8db+NT2Fdy0Ng+AB69bxoPXLZvx9detzAGO4vH6+eubV4Y9PoNBsVxuvooktaCRvFLqJaXUkUk+7gR+CCwFNgFngW9Pdg6t9SNa661a661ut3sh4YgoCc2Vz02z8tGry+f8+uKsVK5c5uLDl5eN1c/DbXmOU0byIiktaCSvtd4+m+cppf4DeGYh1xKxq9xlRyn4wi2rSLXM77/ULx64JMxRTbQyz8Fv9zXQOeAhy26J6LWEiCWRnF2TP+7Lu4EjkbqWiK4lbgcVD23n7ouKoh3KlLaUBur8O442RzkSEUlaa/7nrVqeP3I22qHEjEjeeP2GUuqwUuoQcB3wqQheS0RZtsMa7RCmtbkkg+U5Dv53T320QxERorXm2ztO8bdPHeWHr56OdjgxI2JJXmv9x1rr9VrrDVrrO7TW8qtVRI1SivdvK+FgfTfHmnqjHY6IgH97qZLvvVJFeoqZ020D025kk0xkCqVIGvdsLsRiMvDonjPRDkWE2ZmOQb73ShV3bSrgszeuoH/ES/MUexwkG9kZSiSNjFQLt67L4/F9jTR1D9E37OV7H9g8NjtIxK/vvVKJ0aD44q2rqQ7uNlbV2k/+JJvVJBsZyYuk8pEry0mxGDnV0s/umk4ONXRHOySxQGc6BvntvkY+sK2E3DQby4L7FEy1f3CykSQvksqGogz2PLSdXwanbIZaMIj4FRrFf/zapQC4HBbSbCZJ8kGS5EVSmmmzExEf/H7Nc4ebuWtTwVgjPKUUy3IckuSDJMmLpGQzG3FaTbT3e6IdiliA6vYB+ka8F/RHWpbj4HSbJHmQJC+SmNtppU3KNXHtcGM3ABuLMiYcX5bjoL3fQ/eg/BKXJC+SlmtcH3wRnw7W95BiNo7tFxwSuvkqo3lJ8iKJuZwWufEa5w439rCuMA2TcWIqW+YONLqTurwkeZHE3DKSj2ten5+jTT2sL8y44LHCzBQsJoMkeWQxlEhiLoeVvmEvw6M+bGZjtMMRc3SqpZ/hUT8bi9MveMxoUCxx2fndwbOYjAZu31DAmoK0KEQZfTKSF0krNI1SSjbxKXTTdX3hhUkeApvW5GfY+I9d1Xz6sQOLF1iMkZG8SFouRyjJeyjKTI1yNGK2PviTt8m2W/FpjdNmoizbPunzbt9YwO0bC/jas8f4+Vt1+P06rBvExwtJ8iJpyYKo+NM/4uWNqo6xry9fmj1j4i5z2Rnx+jnbO0xhRvL1spFyjUhaLinXxJ2z3UMAfOLapVy5zMU9m2feqKY8ONKvbR8ISwzDoz5++noNoz5/WM4XaZLkRdJyOQLbALaPG8lrrXl8XwN9w6PRCktMozGY5N+1KodfPHAJ790yc5IvcwWSfE2YkvyrJ1v5yjPH+MPJtrCcL9IkyYukZTUZSbOZJqx6fa2ynU8/dpCfvl4TxcjEVEJJvmAOZZe8NBtWkyFsI/najkEAjsbJ5jOS5EVSczmtE8o1v9nbAMDTB5tkZ6EY1NQ9hNGgxpqRzYbBoCjNTqW2IzxJvi54nqNNPWE5X6RJkhdJbfyCqJ6hUV442ozbaaW6bSBuRmrJpKl7mLw0G8Y5zpIpy7aPjcAXqk5G8kLEj8BIPtDE6tlDZxnx+vnmezdgMih+d7Bp7Hk+v+bO773OUwcaoxWqIFCumc8MmXKXnTMdg/j8C393Fkryjd1DdA96qO8c5IqvvxyzI3tJ8iKpjR/J/3ZfA8tzHFyzws1Vy1387mAT/mBSaOoe4mBDD0/ulyQfTU3dQxRkzL5UE1LmsuPx+WkK1vTna3jUR1PPEFtLMwE41tTL0webaOwe4mhjbI7sJcmLpOZ2Wukf8fKHU23sreviPVuKUEpxx6YCmnqG2XumCzg3M+Odms64mTqXaHx+TXPP8JxuuoaEFkwttC7f0DWI1nDr+nwgULJ55tBZgJhtWy1JXiQ1d3DV6yd+sZdyl50PXFICwPbVuSgFb50OLLypDrasHfD4ONQQm2/LE11b3whev6Ywc37lGlj4XPlQqWZTSQb56TaeOdTE8bOBEXysrreQJC+SmssZmCuvlOI/PrSFNJsZAKfNTGlWKieaAz/ANe0DWE2BH5e3TrdHJ9gkN5/pkyG5aVZsZgM17YP0DI7S2js8rxhCN2/Lsu2sLUjjYPAXfpotdncZW1CSV0rdq5Q6qpTyK6W2nvfYF5VSVUqpk0qpmxYWphCRsTIvjdw0Kw/ft4llOc4Jj63KS+PE2T4gsM3cilwnq/PTePN0x2SnEhEWSvLzufGqlKIs287vD5/l8q/v5L0/emteMdR1DOC0mchMNbOmINAYbWtpJstznXTMYSRf1dpHz+DiLLhb6Ej+CHAPsGv8QaXUGuA+YC1wM/ADpZT0chUxpzAjhd1/s53rV+de8NiqfCc1HQMMerxUtw2wxG3n8qXZ7K3rYnjUF4Vok1vopml++txvvAKsynPS3DtMtsPKmc5BugbmPvKu7RikLNuOUoq1wdbFt67Px+WY2wY09z2ym++8dGrO15+PBSV5rfVxrfXJSR66E3hUaz2ita4BqoBtC7mWEIttdX4aWsOhhh6aeoYodwWS/IjXz/4z3dEOL+k0dQ+RZjPhDJbU5uof7ljLq5+9lq/ctQ6Aky19cz7HmY4BSrMDHUuvWeHmk9uXc+/WIlwO66zLNSNeH+39I2OlwEiLVE2+EKgf93VD8NgFlFIfU0pVKKUq2trioxeESA6r8wIjteePNKN14ObdtvIsjAbF61Xyf3WxNXUPUbiAltAZqRbKXHZW5QXKciebZ5fkewZHeWTXaboGPDR0DY0leZvZyCe3r8BpM5PtsNI16ME7i5lXXQOBMs3ptvCswJ3JjEleKfWSUurIJB93TveySY5NugpBa/2I1nqr1nqr2+2ebdxCRFxRZgp2i5HnjzQDsNTtwGkzs60sixeOtkQ5uuTT2D1M4TzmyJ8vx2klPcXMiVkm+V/sruOffn+CW7/7Gl6/pnSS/vVuhwWtoXNw5tF8x0CgrNPWN7IojfBmTPJa6+1a63WTfDw1zcsagOJxXxcBTVM8V4iYZDAoVuWn0RyciRHqZnjL+jyqWvupap37230xf4GFUAvvB6+UYmWek1OzLNe8UdVOfrqNEW9glD7ZJiVjG9D0zZzkO8fdC6hehNF8pMo1TwP3KaWsSqlyYDnwToSuJUTEhN7a5zitOKyBPXZuWpsHwHOHm6MWV7LpGRylZ2iUonnMkZ/Mqjwnp5r7ZmxCN+TxUVHbxW3r83n6L6/g7969hi3B1a7jZQeTfGiUPp3xSf50W+Q3Gl/oFMq7lVINwGXAs0qpFwC01keBx4BjwPPAg1prmY4g4s7q/EBdfon73OgtN83GltJMnjsiSX6xnG4PJMOlbkdYzrci10nfiHdsWuZUKuo68fj8XLHcRVFmKh+5snzS5mhjexPMYoZNR38cjeS11k9orYu01latda7W+qZxj31Na71Ua71Sa/3cwkMVYvGtzg+M5MtdE5PLLevyOHa2lzNh6mwopne6NZDkl4QpyYfeoZ1q6WNvXRf/vrNy0ue9XtmO2ai4pDxr2vON7TI2i3JNx8AIRoOiLDuV6uAvr0i2tZYVr0JMY2VeGk6biU3F6ROOj5VsjpyNRlhJp7p9ALNRURymcs3y3ECSr6jt4i9/tY9vv3hqrIvkyeY+/uHpo3QPeni9qp3NJZmkWqbfDttpNWExGmY1ku8c8JCZamFZjoPTrYGR/Mf+Zy9ff+7EAv9Wk5MkL8Q0HFYTb37hXdy7pXjC8eKsVFbnp/HKydYoRRafWnuH+fSvD9A/4p3T60639lOabcdkDE/KSk8xU5Bu48e7qmnpHcZsVDy2JzDr+8vPHOW/36zlnh+8ydGmXq5c5prxfEqp4IKoWYzk+z1k2y0scTuo6RjgVEsfLx5rwW6JzHpRSfJCzMBpM2OYpA579XIXe+u6GJhjwkomJ5p7+ejPKxj0BP6Nnj7YxOP7G6mo7ZzTearbB1jiunBWy0KszHPi82seuGoJN6/L54n9jeyu7uCNqg7u2lQwNiq/YvnMSR4u3GVsKp0DHrLsFpa67Xi8fr727HEsRgPvDzbHCzdJ8kLM09Ur3Iz6NLtrpJfNVJ460MSLx1p4vTLQ1O3t6kByn8um2l6fn7qOAZbmhKceH3L96ly2lGbyqe0ruO/iYnqHvTz4q32kp5j56t3refwTV/D3t69hU1HGrM6XbZ9da4POAQ9ZDsvY/YU/nGrjjk0FY9Mww02SvBDztKU0E5vZwK5T0pVyKvvqAv34d1W24fdr9gRH8HNp+VvfNcSoT4d9JH//paX89uOXk2IxctmSbIqzUmjv9/Dhy8twWE0sy3Hwp1eUT/oubjIuh3XCzJmpdAwEyzXj/j5/cnnZfP8aM5IkL8Q82cxGLinP5rVKaXEwGa/PP9Z7f9epdo4399IzFFjhWT2HJB/q5R+umTWTMRgUH7q0jMxUM386z4TrclrpGBjB59f8aveZSbtSjvr89AyNkmW3jH1sK8tiXWH6JGcMj+lvGQshpnX1CjdfeeYYjd1DNHUPkZ9uo2gB/VUSyYnmPoZGfWwtzaSirov/q2gAAq1557JDU2jB0FJ3eEfy53vgqnL++LJSbOb53QDNtlsY9Wn+640avvrscYZGffzZleUTntMVbHuQbbcE9zDYSm5aZMo0ITKSF2IBrg7elHvfj97i3h+9xd88cSTKEcWO/cGtEz+5fQUAv9xdR3FWCpcvc9HYNcSId3brI6vbBsi2W8hItUQsVgjMkJlvgofAVpIA33wh0Ji3vvPCNRSh1a5Z9sBzt5RmRnxQIEleiAVYluOg3GVnxOtjQ1E6e2s7Z9WJMBnsO9ONy2HlimXZlGSlMurTXFKezRKXHb+ePAlO5nRb/4QVx7EqdON0xOsnPcVMQ9eFq2k7+0NJPrK/sMaTJC/EAiilePLBK3j98+/igauWMODxcezs4vQJj3X7z3SxuSQDpRTXrAh0mL2kPGus0dtsl/RXtw2ErZ1BJIVG8vdsLuTisiwaui78JdYRHMlnOyTJCxE30lPM2MxGtpUFlr6/UzO3OeDx5HcHm/juFC0AxuvoH6G2Y5DNwWZed2wqwO20cvUKN+XBLo6zqcvvru6gY8Az1kMoli3PcfDVu9bxd+9eQ1FmCg1dQxe0KzhXrpEkL0TcyUu3UZKVOjZNMBF9d2clP3z1NH7/9L1WQjtnbS4JJPmLy7LY89B2ctNspKeaybJbZpwrP+rz83dPHaUwI4X3bS2e9rmxQCnF/ZeWkpFqoSgzhf4R79hsopCOAQ9KQWaE7y+MJ0leiDC6uCyLitquiDacipba9gEqW/sZGvVNWm8e77kjzdgtRjYUTT41sCw7dcYk/7M3aznZ0sff376GlAgt+Y+U4qzAzdT6zsC/074zXfSPeOkcGCEjxTxpJ8tIkSQvRBhtK8+kY8CzaFu7LaaXjp/bDWu6/VG7Bz08c6iJuzcXTjlbpdzlmDbJt/eP8G8vVXLdSjc3rLlwk/VYF+p739A1SPegh3t/9BYPPXF4rKXBYpIkL0QYXRysyydiyWbn8VaKswLJa7pdlX67r5ERr58PbCud8jnlrlRaekem7Pvzo1dPM+jx8tBta1Bq8Ua94RKaFtnQNcT+M934/JqnDzax/0w32fbIzos/nyR5IcKo3GXH5bAk3M3XnsFR3qnt5PYNBRRmpEyZ5LXW/HJ3HZtLMlhTMPXN0lB//slG8809w/zP23Xcs7mIZWHuV7NY0lPMpNlM1HcNsreuC6NB4bCaONszLCN5IeKZUorLl7p4rbJ9xpuT8eTVU634/Jrta3JZmefk5BSbYL9d3Ul12wAfvGTqUTzAitxA8p5sM+3vv1KFz6/5q+uXLzzwKCrKTKWha4i9dV2syU/jL65ZCkDWIk6fBEnyQoTdtSvdtPePcLQpMebLt/eP8MiualwOC5uKMliR66S6bYDRSRZ9PX2wEYfVxG0b8qc95xK3A7vFyKGG7gnHOwc8PLrnDPduLR67eRmvijJTqG0f4EB9N1tKM/nTK8pYmeucdVfLcJEkL0SYXb3CjVLwagJsKFLZ0sed33uD0239/PM9GzAYFCvzHHiC7X/H8/s1Lx1v5ZoV7hnbAxgNirWF6RwMNjALOdLYw6hPc/vG6X9JxIPirFSq2wcYGvWxuTSwu9QLn7qa9128uNNBJckLEWYuh5UNhekJsWvUl585xtCoj8f+/LKxWS7LcwJb551s7uexPfV8/Bd78fk1hxp7aOsbmfVsmI1F6Rw/24vHe+4dQagMtCov9hc/zaRo3FaFW4KLwqJBulAKEQHXrszhuy9X0jXgIXORb7SF0/GzfWxfncOGcSWGZTkODAqePdzES8da8fj8PHWgkeq2AYwGxbUr3bM694aiDDzeGk619I212j3R3EeO07roNycjITTDJj/dRmFGePamnQ8ZyQsRAdeudKN1YLOMeNU96KG9f2Rs5B5iMxspy7bz+8PNOGwmVuQ6eHhnJS8cbebissxZd4vcGPzFcXBcXf5Ecy8r85yTvyDOhKabbo7iKB4kyQsRERuKMsiyW3jxWMvMTw6zQY+XP/vvPRxt6pn5ydOoag30cZ9sGuOK3EAi/sqd6/jrm1ZR1zFIZWs/N6zJm/X5i7NSyEg1c6g+EKfX56eytZ9VCZLkS7JSSbOZuHbF7N7ZRIqUa4SIAKNBcdemQv7rzRoeqO9mU3HGol37lRNt7DzRysbiDNYWzH/HocppkvxHrixnc2kGt23IR2vNxuIMDtZ3s311zqzPr5RifWE6hxoDSb62YxCP158Q9XiAVIuJdx7ajtUU3bG0jOSFiJBP3bCcXKeNLz5+eNLphpHy/NFm4NyOSvNV2dJPitk4aT15W3kWH7s6MO9bKcW/vGc9X7ptNaXZc+v7vrEog1MtfQx5fGM3XROlXAOB0la0V+wuKMkrpe5VSh1VSvmVUlvHHS9TSg0ppQ4EP3608FCFiC9Om5l/uGMtx8/28h+vVS/KNUe8Pl45EZjVs9AkX9XWz9Ic+6w2sl6Vl8YDVy2Z8zU2FKXj82v2n+niZHMvRoOK21WusWqh5ZojwD3Ajyd57LTWetMCzy9EXLt5XR63rMvj2ztOsa4gnasjXJ99s6qD/hEvS912TrcO4PfrWSXpyVS19HHJkuwwRzjRZUuzcTksfHPHSbLtVsqyUxe0BZ+40IJG8lrr41rrk+EKRohE9I33bmB5joNP/HIfxyO8a9TzR5pxWk3cf2kpQ6M+mnuH53We/hEvTT3DER9VO21mPn/zKvaf6eblEy2sioPNQeJNJGvy5Uqp/UqpPyilrprqSUqpjymlKpRSFW1t8TvdTIipOG1m/utPL8ZuNfLgr/bhC/a0OdzQw2MV9WG7js+vefF4C9etyhm7eTnfks3paW66htt7NhexuSQDv4ZVuYlTj48VMyZ5pdRLSqkjk3zcOc3LzgIlWuuLgE8Dv1JKTforWmv9iNZ6q9Z6q9sd3alGQkRKfnoKf/vuNVS3DfDisWZ8fs0nf72fL/z2EB39I2G5xqmWPjoHPLxrVQ5LcwI3QEPJeq6mmz4ZbgaD4it3rcNpM3Hp0siWh5LRjDV5rfX2uZ5Uaz0CjAQ/36uUOg2sACrmHKEQCeLmtXkUZ6Xw413V9A17xzYW2Xm8NSz9TEK9ZJblOHA7rDhtpnlvXlLZ2o/ZqChdpCZhawvSOfh3N877/oGYWkTKNUopt1LKGPx8CbAcWJzpBULEKJPRwANXLmH/mW6+8swx1hemU5iRwo5jzZM+3+vz09o3+5p6XccgACXZqSilWOp2UN0+35F8H+UuOybj4s2ylgQfGQudQnm3UqoBuAx4Vin1QvChq4FDSqmDwG+Av9BaJ9YuCkLMw71bi8hINdM77OUzN67ghjW57Kpsn7BD0qjPz3++XsM133yVK7/+Cg1dgxPO0Tc8yr/uOEn/ebsq1XUOkplqJs1mBmCp28Hp1vmN5Os6Bil3zW3Ou4hNC51d84TWukhrbdVa52qtbwoe/63Weq3WeqPWerPW+nfhCVeI+JZqMfGFm1fx/m3FXLPCzY1rc/F4/ew6dW7SwX+/UcuXnzlGTpoVj8/P0webJpxjx9EWvvtyFd97uWrC8bqOAUrGLUZammOnuXf4gl8GM9Fa09g9RGFGfPdzFwGy4lWIRXbfthL++Z4NKKXYVpZFZqqZF46eK9nsqmxjZa6TJz5xBZtLMnj6wMQkH9qM5D/fqKGxe2jseF3H4IQa+pLgFnvVc5xh0z04yqDHR0GGbc5/NxF7JMkLEUUmo4HrV+ey80QrHq+fUZ+fitouLl0S2BD8zk2FnGjum7Dd3tGmHkqzU1HAt14ILFPxeP00dQ9Rmn0uyS8LzrCpmuMMm9AvjvH90EX8kiQvRJTduj6PvmEvb1S1c6ihh6FR39hK01vX52M0KJ4+2AgEdl861tTLVctd/NmV5Tyxv5GTzX00dg/h10zoHVOabSfFbORgffec4gkleSnXJAZJ8kJE2ZXL3DhtJp45dJbdNR1AoAEYgNtp5YplLp4+2ITWmvquQfpGvKwtSOcjV5YD8OKx5rHpk+NH8majga1lmbxdPbc5D41dwSQvI/mEIEleiCizmAzctDaPHceaee1UO8tzHLgc1rHH79xYQH3nEBV1XWP1+LUFabgcVtYVprHrVPvY9Mnz57VfuiSbky19c1pw1dg9RIrZSGaqOQx/OxFtkuSFiAG3bcinb9jLW9UdXBKsx4fcvC4Pu8XI/1XUc7SpB6NBjW3acdVyN/vOdHGsqZcUsxG30zrhtZcGyz67a2Y/mm/sGqIgwxb1FrkiPCTJCxEDrljqIj0lMHK+9LzOj3arids25PPsobPsqelieY5jrFPj1cvdeP2aZw41UZKVekFi3lCUTqrFyNvVHbOOpbF7iMJMqccnCknyQsSAQMkmFzhXjx/vfVuLGfD4eKe2kzUF59pAbSnNxG4xMuDxUZJ9YWIO1OWzeOv07JN8U/dQVDeeFuElSV6IGPHZG1fyo/u3kOO8cH76ltJMlgRXoI7f0s9iMnBZsKlX2SRJHuDSJVlUtvbTPou6/JDHR8eAR6ZPJhBJ8kLEiJw0Gzevm3wjbKUU79lSBMC6gokNXUMbkZRMsfXeZcHyz2xKNqHpk7IQKnFIkhciTvzpFWX88z3rubhsYjnnxjV5LMtxcMkkZR6A9YXp2MwG9tV1z3gNmSOfeBa6/Z8QYpGkWky8f1vJBcfz0m289OlrpnydyWhgdX4aRxp7ZryGzJFPPDKSFyIJrC9M52hTD/7grlRTaewexGhQ5J43FVPEL0nyQiSBdYXpDHh8VLdP33q4qXuYvDTbovaRF5El30khksD6wsCMnJlKNo1dMn0y0UiSFyIJLM9xYDUZODxDkj/TOUhRliT5RCJJXogkELr5Ol2SH/L4aO4dpnyKqZgiPkmSFyJJrC9M51hT75Q3X+s6A/X6Mtn2L6FIkhciSawvTKd/xEttx+Q3X2vbA50sy2Qkn1AkyQuRJNYWBlbKTlWyCSX/UpcshEokkuSFSBIrcp2YDGrCVoLj1XUMkG23kGaTPvKJRJK8EEnCbDTgclhp7Zu8UVlt++CEnaVEYpAkL0QScTktU3ajrO0YkHp8ApIkL0QScTuskyb54VEfZ3uGZWZNApIkL0QScTmstE1SrjnTGdwjVso1CWdBSV4p9U2l1Aml1CGl1BNKqYxxj31RKVWllDqplLppwZEKIRbM5bTS0e+5YK58TbCnTbmM5BPOQkfyLwLrtNYbgFPAFwGUUmuA+4C1wM3AD5RSxgVeSwixQG6HFa9f0zM0OuF4XWj6ZJYk+USzoCSvtd6htfYGv3wbKAp+fifwqNZ6RGtdA1QB2xZyLSHEwrmCLYTPr8vXdgySmWomPVWmTyaacNbkPwI8F/y8EKgf91hD8NgFlFIfU0pVKKUq2trawhiOEOJ8LocF4IK6fG37gNx0TVAzJnml1EtKqSOTfNw57jkPAV7gl6FDk5xq0oYZWutHtNZbtdZb3W73fP4OQohZcjsCI/m2cSP5Ea+Po029LHM7ohWWiKAZt//TWm+f7nGl1IeBdwPXa61DibwBKB73tCKgab5BCiHCwz1WrvGMHXv5eCs9Q6PctiE/WmGJCFro7Jqbgc8Dd2itB8c99DRwn1LKqpQqB5YD7yzkWkKIhUtPMWM2qgk1+d/uayDHaeWq5fJOOhEtdCPv7wFW4EWlFMDbWuu/0FofVUo9BhwjUMZ5UGvtW+C1hBALpJQi235urnx7/wivnmzjz64sx2iYrMoq4t2CkrzWetk0j30N+NpCzi+ECL/xrQ2eOtCE1695z5aiGV4l4pWseBUiyYxvbfD4vgY2FKWzItcZ5ahEpEiSFyLJuBxW2vs8tPePcLSpl5vW5kU7JBFBkuSFSDIuZ2Ak/9bpDgAuW5od5YhEJEmSFyLJhFobPH+0GbvFyPrC9GiHJCJIkrwQSSbU2mDn8RYuLs/CbJQ0kMjkuytEkgm1Nhge9XPZEinVJDpJ8kIkmVBrA5B6fDKQJC9Ekgm1NnDaTKwtkHp8opMkL0SSCbU2uKQ8S1a5JoGFtjUQQsQZpRRfvGU1m0oyoh2KWASS5IVIQh+5sjzaIYhFIuUaIYRIYJLkhRAigUmSF0KIBCZJXgghEpgkeSGESGCS5IUQIoFJkhdCiAQmSV4IIRKY0lpHO4YxSqk2oG4Bp3AB7WEKJxJiPT6QGMNFYgwPiXF2SrXW7skeiKkkv1BKqQqt9dZoxzGVWI8PJMZwkRjDQ2JcOCnXCCFEApMkL4QQCSzRkvwj0Q5gBrEeH0iM4SIxhofEuEAJVZMXQggxUaKN5IUQQowjSV4IIRJYQiR5pdTNSqmTSqkqpdQXoh0PgFKqWCn1ilLquFLqqFLqr4LHs5RSLyqlKoN/ZkY5TqNSar9S6pkYjS9DKfUbpdSJ4L/lZTEY46eC3+MjSqn/VUrZoh2jUuo/lVKtSqkj445NGZNS6ovBn5+TSqmbohjjN4Pf60NKqSeUUhmxFuO4xz6rlNJKKVc0Y5xJ3Cd5pZQR+D5wC7AGeL9Sak10owLAC3xGa70auBR4MBjXF4CdWuvlwM7g19H0V8DxcV/HWnwPA89rrVcBGwnEGjMxKqUKgf8HbNVarwOMwH0xEON/Azefd2zSmIL/L+8D1gZf84Pgz1U0YnwRWKe13gCcAr4YgzGilCoGbgDOjDsWrRinFfdJHtgGVGmtq7XWHuBR4M4ox4TW+qzWel/w8z4CyamQQGw/Cz7tZ8BdUQkQUEoVAbcBPxl3OJbiSwOuBn4KoLX2aK27iaEYg0xAilLKBKQCTUQ5Rq31LqDzvMNTxXQn8KjWekRrXQNUEfi5WvQYtdY7tNbe4JdvA0WxFmPQd4C/BsbPXIlKjDNJhCRfCNSP+7oheCxmKKXKgIuA3UCu1vosBH4RADlRDO3fCPxH9Y87FkvxLQHagP8KlpR+opSyx1KMWutG4FsERnRngR6t9Y5YinGcqWKK1Z+hjwDPBT+PmRiVUncAjVrrg+c9FDMxjpcISV5Ncixm5oUqpRzAb4FPaq17ox1PiFLq3UCr1npvtGOZhgnYDPxQa30RMED0y0cTBOvadwLlQAFgV0rdH92o5izmfoaUUg8RKHn+MnRokqcteoxKqVTgIeDvJnt4kmNRz0WJkOQbgOJxXxcReLscdUopM4EE/0ut9ePBwy1Kqfzg4/lAa5TCuwK4QylVS6DE9S6l1C9iKD4IfG8btNa7g1//hkDSj6UYtwM1Wus2rfUo8DhweYzFGDJVTDH1M6SU+jDwbuCD+txCnliJcSmBX+gHgz87RcA+pVQesRPjBImQ5PcAy5VS5UopC4EbH09HOSaUUopALfm41vpfxz30NPDh4OcfBp5a7NgAtNZf1FoXaa3LCPybvay1vj9W4gPQWjcD9UqplcFD1wPHiKEYCZRpLlVKpQa/59cTuP8SSzGGTBXT08B9SimrUqocWA68E4X4UErdDHweuENrPTjuoZiIUWt9WGudo7UuC/7sNACbg/9XYyLGC2it4/4DuJXAnfjTwEPRjicY05UE3qodAg4EP24FsgnMbKgM/pkVA7FeCzwT/Dym4gM2ARXBf8cngcwYjPEfgRPAEeB/AGu0YwT+l8A9glECiejPpouJQAniNHASuCWKMVYRqGuHfmZ+FGsxnvd4LeCKZowzfUhbAyGESGCJUK4RQggxBUnyQgiRwCTJCyFEApMkL4QQCUySvBBCJDBJ8kIIkcAkyQshRAL7/wH/k5RbdKQZRAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data = hal.get_data_linker(\n", " {arma_graph.a: np.ones(p)/p,\n", " arma_graph.b: np.ones(q),\n", " arma_graph.c: 0.},\n", " n_data=150\n", ")\n", "m = hal.get_generative_model(arma_graph, data)\n", "plt.plot(m.get_example(arma_graph.y))" ] }, { "cell_type": "markdown", "id": "4b639643", "metadata": {}, "source": [ "## the ARIMA model\n", "The arima model is an integrated ARMA model. So the ARMA process is simply intergrated $d$ times.\n", "\n", "In Halerium it can be defined like this" ] }, { "cell_type": "code", "execution_count": 9, "id": "7ef900f9", "metadata": {}, "outputs": [], "source": [ "d = 3\n", "p = 7\n", "q = 5\n", "sigma = 1.\n", "sigma_obs = 1.e-2\n", "\n", "arima_graph = hal.Graph(\"graph\")\n", "with arima_graph:\n", " hal.StaticVariable(\"a\", shape=(p,), mean=0, variance=1)\n", " hal.StaticVariable(\"b\", shape=(q,), mean=0, variance=1)\n", " hal.StaticVariable(\"c\", shape=(), mean=0, variance=1)\n", " \n", " hal.Variable(\"epsilon\", mean=0, variance=sigma**2)\n", " \n", " hal.Variable(\"x\")\n", " \n", " ar_contrib = sum([\n", " a[i] * x[hal.TimeIndex-i-1]\n", " for i in range(p)\n", " ])\n", " \n", " ma_contrib = sum([\n", " b[i] * epsilon[hal.TimeIndex-i-1]\n", " for i in range(q)\n", " ])\n", " \n", " x.mean = c + ar_contrib + ma_contrib\n", " \n", " # x follows the ARMA model\n", " \n", " # to integrate we define a helper variable int_i\n", " temp = x\n", " for i in range(d):\n", " temp_integrator = hal.Variable(f\"int_{i}\")\n", " temp_integrator.mean = temp_integrator[hal.TimeIndex-1] + temp\n", " temp = temp_integrator\n", " \n", " hal.Variable(\"y\")\n", " y.mean = temp\n", " \n", " hal.Variable(\"observed\", mean=y, variance=sigma_obs**2)\n", " # this \"observed\" variable is strictly speaking not part of the model,\n", " # but it is recommended when we want to fit data to the model as it allows us\n", " # to include observation noise." ] }, { "cell_type": "code", "execution_count": 10, "id": "dec676a9", "metadata": {}, "outputs": [], "source": [ "data = hal.get_data_linker(\n", " {arma_graph.a: np.ones(p)/p,\n", " arma_graph.b: np.ones(q),\n", " arma_graph.c: 0.},\n", " n_data=150\n", ")\n", "m = hal.get_generative_model(arima_graph, data)\n", "x_ex, y_ex = m.get_example([arima_graph.x, arima_graph.y])" ] }, { "cell_type": "code", "execution_count": 11, "id": "18384029", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgE0lEQVR4nO3deXxV9Z3/8dcnNwkhYSdhMWEXF0QRiKDSWmx1xGWkm63bdBkttuNM2xlnuv6m0878Or/+Op1OtXVkrN2sVluXWuqGa+sKEhCQsIY9kEASSAImZLn3M3/cG73EhAS4N+fem/fz8biPe+85J/e8We47537vWczdERGR9JcVdAAREUkMFbqISIZQoYuIZAgVuohIhlChi4hkCBW6iEiGCLTQzeznZrbfzNb1cvlPmNl6Mys3s98kO5+ISDqxIPdDN7OLgMPAve4+vYdlpwK/Az7o7gfNbJS77++LnCIi6SDQLXR3fwk4ED/NzKaY2dNmttLMXjazM2KzPgfc6e4HYz+rMhcRiZOKY+h3A3/n7rOBfwT+Ozb9NOA0M3vVzJaZ2YLAEoqIpKDsoAPEM7NBwIXAQ2bWMXlA7D4bmArMB0qAl81survX93FMEZGUlFKFTvQTQ727n9vFvEpgmbu3AdvNbBPRgl/Rh/lERFJWSg25uHsj0bK+BsCiZsRmPwZcHJteSHQIZlsQOUVEUlHQuy0+ALwOnG5mlWZ2E3ADcJOZrQHKgYWxxZcCdWa2HngR+Cd3rwsit4hIKgp0t0UREUmclBpyERGRExfYl6KFhYU+ceLEoFYvIpKWVq5cWevuRV3NC6zQJ06cSFlZWVCrFxFJS2a2s7t5GnIREckQKnQRkQyhQhcRyRAqdBGRDKFCFxHJECp0EZEMoUIXEckQKnQRkT50+3NbeG1rbVJeW4UuItJH9h86wo+e30zZjoNJeX0VuohIH3lhw37c4dJpo5Py+ip0EZE+8tyGfRQPG8gZYwYn5fVV6CIifaC5NczLW2q5dNpo4i6xmVAqdBGRPvDylhpa2iNJG24BFbqISJ94bsM+BudlM2fSiKStQ4UuIpJk4Yjzwsb9zD99FDmh5NWuCl1EJMlW7jxI7eHWpA63gApdRCTplpZXkxvK4uLTu7zQUMKo0EVEksjdWVpezbxTRzI4Lyep61Khi4gk0fqqRioPNnPZWWOSvi4VuohIEi1dV02WwSVJHj8HFbqISFItLd9H6cQRFA4akPR1qdBFRJJke+3bbNp3qE+GW0CFLiKSNEvLqwG47KzkD7eACl1EJGmWllczvXgIJcPz+2R9KnQRkSTY13iEN3fVs6CPhlugF4VuZuPM7EUz22Bm5Wb2pS6WmW9mDWa2Onb7VnLiioikh2fW7wPos/FzgOxeLNMO3Obuq8xsMLDSzJ519/WdlnvZ3a9KfEQRkfSzdF01kwsLOHXUoD5bZ49b6O5e5e6rYo8PARuA4mQHExFJV/VNrSzbVsdl08ck7dznXTmuMXQzmwjMBJZ3MfsCM1tjZk+Z2Vnd/PwiMyszs7KamprjTysikgaeWb+P9ohz+fS+G26B4yh0MxsEPAJ82d0bO81eBUxw9xnAj4HHunoNd7/b3UvdvbSoKLknqRERCcoTa6soGT6Qs4uH9ul6e1XoZpZDtMzvd/dHO89390Z3Pxx7/CSQY2aFCU0qIpIG6ptaebWilivPGdunwy3Qu71cDPgZsMHdf9jNMmNiy2Fmc2KvW5fIoCIi6eCZ8uhwy5Vnj+3zdfdmL5d5wF8Bb5nZ6ti0bwDjAdx9MfBx4Atm1g40A9e6uyc+rohIanvirSrGjej74RboRaG7+yvAMT83uPtPgJ8kKpSISDrqGG656f2T+ny4BXSkqIhIwnQMt1x19imBrF+FLiKSIB3DLdOLhwSyfhW6iEgCvLN3y9mnBDLcAip0EZGECHLvlg4qdBGRBHj8rSrGj8gPbLgFVOgiIift4NutvFZRyxVn9/3BRPFU6CIiJ+mpddXRvVvOCW64BVToIiInbcmaPUwpKuCsU4IbbgEVuojISalqaGb59gMsPLc40OEWUKGLiJyUx9dU4Q5XzwjmYKJ4KnQRkZPwhzV7mFEylImFBUFHUaGLiJyoiv2HWbenkavPTY2LuKnQRURO0JI1ezGDvwx475YOKnQRkRPg7ixZvYcLp4xk1JC8oOMAKnQRkROytrKBHXVNLJyRGsMtoEIXETkhf1i9l9xQFpf18YWgj0WFLiJynMIR549r93LxGUUMHZgTdJx3qNBFRI7T61vrqDnUwtUpNNwCKnQRkeP26KpKBudl86EzRwUd5SgqdBGR43C4pZ2n1lVz1TmnkJcTCjrOUVToIiLH4el11TS3hfnYrNQabgEVuojIcXlkZSUTR+Yze8LwoKO8hwpdRKSXKg828fq2Oj46qyTwMyt2RYUuItJLj725B4CPzEy94RZQoYuI9Iq788iqPcydNIJxI/KDjtMlFbqISC+s2lXP9tq3+djskqCjdKvHQjezcWb2opltMLNyM/tSF8uYmd1hZhVmttbMZiUnrohIMB5dVUleThaXp9Ch/p1l92KZduA2d19lZoOBlWb2rLuvj1vmcmBq7DYXuCt2LyKS9o60hfnjmr0sOGsMg/NS51D/znrcQnf3KndfFXt8CNgAdP5GYCFwr0ctA4aZWWqcIFhE5CQ9t2EfjUfa+eis1B1ugeMcQzezicBMYHmnWcXA7rjnlby39DGzRWZWZmZlNTU1xxlVRCQYv12xm1OG5jHv1MKgoxxTrwvdzAYBjwBfdvfGzrO7+BF/zwT3u9291N1Li4qKji+piEgAKg828UpFLdeUjiOUlXr7nsfrVaGbWQ7RMr/f3R/tYpFKYFzc8xJg78nHExEJ1kNllQBcU5rawy3Qu71cDPgZsMHdf9jNYkuAT8X2djkfaHD3qgTmFBHpc+GI81DZbt53aiElw1Nz3/N4vdnLZR7wV8BbZrY6Nu0bwHgAd18MPAlcAVQATcBnE55URKSPvVJRy96GI3zzymlBR+mVHgvd3V+h6zHy+GUcuDVRoUREUsFvV+xiREEul0xLrfOed0dHioqIdKHucAvPrt/HR2YWMyA7tc573h0VuohIF37/5h7aws4nzxvX88IpQoUuItKJu/Pgit3MGj+M00YPDjpOr6nQRUQ6WbWrnor9h7n2vPFBRzkuKnQRkU7uX76TgtwQV56TXmcwUaGLiMQ5+HYrj6+t4qOzSigY0Js9u1OHCl1EJM7DKytpbY9w4/kTgo5y3FToIiIxkYhz3/KdnDdxOKePSZ8vQzuo0EVEYl6pqGVnXVNabp2DCl1E5B33LdvJyIJcFqTwVYmORYUuIgJUNTTz3IZ9fOK8cWlzZGhnKnQREeCB5btw4Po56bXveTwVuoj0e23hCA+u2M3804oYNyL1T5PbHRW6iPR7z67fx/5DLWn7ZWgHFbqI9Hu/fHUHJcMHMv/09DhNbndU6CLSr63b08AbOw7w6Qsmpvw1Q3uiQheRfu0Xr+4gPzfEJ9LoNLndUaGLSL9Vc6iFP67Zy8dmlTB0YE7QcU6aCl1E+q3fLN9FazjCZ+ZNDDpKQqjQRaRfamkPc9/yncw/vYgpRYOCjpMQKnQR6ZeeWFtFzaEWPjtvUtBREkaFLiL9jrvzi1d3MKWogIumFgYdJ2FU6CLS76zceZC39jTw2XmTMEvvXRXjqdBFpN/52SvbGZKXzUdnFQcdJaFU6CLSr+yofZuny6u58fwJ5Oem1yXmeqJCF5F+5acvbyMnKytjdlWM12Ohm9nPzWy/ma3rZv58M2sws9Wx27cSH1NE5OTVHm7h4ZWVfGx2MaMG5wUdJ+F683njl8BPgHuPsczL7n5VQhKJiCTJva/toDUc4eb3Tw46SlL0uIXu7i8BB/ogi4hI0jS1tnPvsp1ceubojDmQqLNEjaFfYGZrzOwpMzuru4XMbJGZlZlZWU1NTYJWLSLSs9+t2E19Uxu3fCAzt84hMYW+Cpjg7jOAHwOPdbegu9/t7qXuXlpUVJSAVYuI9Kw9HOGeV7ZTOmE4syeMCDpO0px0obt7o7sfjj1+Esgxs8w59EpE0t6T66qpPNjMoosyd+scElDoZjbGYodamdmc2GvWnezriogkgrtz15+2MrmogEvOHB10nKTqcS8XM3sAmA8Umlkl8C9ADoC7LwY+DnzBzNqBZuBad/ekJRYROQ7Pb9jPhqpG/vOaGWSl+RWJetJjobv7dT3M/wnR3RpFRFKKu/PjF7YwbsRArj73lKDjJJ2OFBWRjPXyllrWVDbwN/NPJSeU+XWX+X9CEemXOrbOxw7Ny7iTcHVHhS4iGWn59gOs2HGQz39gCgOyQ0HH6RMqdBHJSD9+YQtFgwfwyfPGBR2lz6jQRSTjrNx5kFcr6lj0/snk5fSPrXNQoYtIBvrxC1sYUZDLDeePDzpKn1Khi0hGWbnzAH/aVMPn3j854y5g0RMVuohklB8s3UzhoAF8+sIJQUfpcyp0EckYr1XU8vq2Om69eEq/2zoHFbqIZAh35wfPbGLs0Dyum9O/xs47qNBFJCP8aVMNq3bV83cfnNqv9myJp0IXkbQXiUS3zsePyOea0pKg4wRGhS4iaW9peTXlexv58iVT+8U5W7rTf//kIpIR2sMR/vPZzUwpKmDhuf3jnC3dUaGLSFp7aGUlFfsP80+XnU4ow8933hMVuoikrabWdn747GZmTxjOZWeNCTpO4FToIpK2fvrSdmoOtfCNK84gdiXMfk2FLiJpqeZQC//z0lYunz6G2RNGBB0nJajQRSQt/ei5zbS2R/jKgjOCjpIyVOgiknYq9h/mwRW7uWHueCYVFgQdJ2Wo0EUk7fz/pzcyMCfEFz80NegoKUWFLiJp5ZUttTy7fh9fmD+FkYMGBB0npajQRSRttIUjfOeP5Ywfkc9N75sUdJyUo0IXkbTx69d3smX/Yf75qmn99gRcx6JCF5G0UHu4hf96bjMXnVbEJWeOCjpOSlKhi0ha+MHSTTS3hvnWVdN0EFE3eix0M/u5me03s3XdzDczu8PMKsxsrZnNSnxMEenP1lbW89uy3XzmwomcOmpQ0HFSVm+20H8JLDjG/MuBqbHbIuCuk48lIhIViTjfXlLOyIJcvniJdlM8lh4L3d1fAg4cY5GFwL0etQwYZmZjExVQRPq335XtZtWuer6y4AyG5OUEHSelJWIMvRjYHfe8MjbtPcxskZmVmVlZTU1NAlYtIpms5lAL//7kBuZMGsE1s/vvlYh6KxGF3tW3E97Vgu5+t7uXuntpUVFRAlYtIpnsu0+sp7ktzL9/5Gx9EdoLiSj0SmBc3PMSYG8CXldE+rFXttTy2Oq9fOEDU/RFaC8lotCXAJ+K7e1yPtDg7lUJeF0R6aeOtIX5P4+9xcSR+fzNxacGHSdtZPe0gJk9AMwHCs2sEvgXIAfA3RcDTwJXABVAE/DZZIUVkf7hzhcr2FHXxP03z9URocehx0J39+t6mO/ArQlLJCL92qbqQyz+81Y+OrOYeacWBh0nrehIURFJGW3hCLc9tJoheTl888ozg46TdnrcQhcR6SuL/7SVdXsaueuGWTo17gnQFrqIpIQNVY3c8cIW/nLGKVx+to5NPBEqdBEJXFs4wj8+tIahA3P4ztVnBR0nbWnIRUQCd9eftlK+t5HFN85mREFu0HHSlrbQRSRQ6/Y0cMfzW7h6xiksmD4m6DhpTYUuIoFpam3niw++ychBuRpqSQANuYhIYP7t8Q1sr32b+2+ay3ANtZw0baGLSCCeXlfNA2/s4paLpnChDiBKCBW6iPS5qoZmvvboWs4pGco/XHpa0HEyhgpdRPpUOOL8w2/X0Noe4fZrZ5KbrRpKFI2hi0if+skLFby+rY7vf+wcJhUWBB0no+hXo4j0mZc21/Cj5zfz0ZnFXFOqKxAlmgpdRPrEnvpmvvTgm5w2ajDf1RWIkkKFLiJJ19oe4db7V9EWdu66cRYDc3WO82TQGLqIJN13n1jP6t31LL5xFpOLdDm5ZNEWuogk1aOrKvnV6zu5+X2TWDBdZ1FMJhW6iCTNql0H+dojb3H+5BF89fIzgo6T8VToIpIUe+qbWXTvSsYMzeOuG2aTE1LdJJvG0EUk4Zpa2/ncr8poaQvzwOd0npa+okIXkYSKxI4E3VjdyM8+cx5TRw8OOlK/oc9AIpJQ31+6iafLq/nGFWdy8emjgo7Tr6jQRSRhfvHqdhb/eSs3zB3PTe+bFHScfkeFLiIJ8fjavfzr4+v5i2mj+deF03UkaABU6CJy0l7bWss//HYNs8cP547rZhLKUpkHQYUuIidl/d5Gbrl3JRNG5nPPp0vJy9Fh/UHpVaGb2QIz22RmFWb2tS7mzzezBjNbHbt9K/FRRSTVbNl3iBt/tpxBedn86q/nMCxfuycGqcfdFs0sBNwJXApUAivMbIm7r++06MvuflUSMopICtpWc5jr71lOKMu4/+a5nDJsYNCR+r3ebKHPASrcfZu7twIPAguTG0tEUtmuuiau/+lyIhHnNzfP1Qm3UkRvCr0Y2B33vDI2rbMLzGyNmT1lZmd19UJmtsjMysysrKam5gTiikjQ9tQ3c91Pl3GkPcx9N8/VgUMppDeF3tXX1d7p+SpggrvPAH4MPNbVC7n73e5e6u6lRUVFxxVURIK3s+5tPrH4dRqPtHHfTXM5c+yQoCNJnN4UeiUwLu55CbA3fgF3b3T3w7HHTwI5ZlaYsJQiErjN+w5xzeLXaWpt5zc3n8/04qFBR5JOelPoK4CpZjbJzHKBa4El8QuY2RiLHUVgZnNir1uX6LAiEoy3Khv45P+8DsBvb7mAs0tU5qmox71c3L3dzP4WWAqEgJ+7e7mZfT42fzHwceALZtYONAPXunvnYRkRSUNvbD/ATb9cwdD8HO6/eS4TRhYEHUm6YUH1bmlpqZeVlQWybhHpnSfWVvH3v1tNyfCB3H/zXMYO1a6JQTOzle5e2tU8nT5XRN7D3bn7pW38v6c2UjphOD/9VKnOaZ4GVOgicpT2cIRv/7Gc+5bt4spzxvKf18zQ4fxpQoUuIu9oaG7jyw++yYubarjlA5P56mVnkKUTbaUNFbqIANHzsiz69Up2H2jiux+Zzg1zJwQdSY6TCl1EeHpdFbf9bg0Dc7N5YNH5nDdxRNCR5ASo0EX6sfZwhP96bjN3vriVc8cNY/GNsxkzNC/oWHKCVOgi/VTlwSa+/OBqynYe5Lo54/j21WcxIFtffqYzFbpIP/TUW1V89ZG1RBxuv/ZcFp7b1fn2JN2o0EX6kbdb2vm/T2zggTd2MaNkKHdcN1NHfmYQFbpIP/FaRS1feWQte+qbueUDk7nt0tPJzdZVKDOJCl0kwx1uaed7T23gvmW7mFRYwEO3XECp9mLJSCp0kQz27Pp9fHtJOXsbmrn5fZO47S9OZ2CuvvjMVCp0kQy0q66J7/yxnOc37ue00YN4+PMXMHuCtsoznQpdJIM0t4a5+6Vt/PefKsjOMr55xZl8Zt5EckIaK+8PVOgiGSAccR5ZVckPn9lMdeMRrjxnLP985TQdJNTPqNBF0pi78+fNNXzvqY1srD7EueOGccd1M5kzScMr/ZEKXSQNuTsvb6nl9ue3sHLnQSaMzOfO62dxxdljiF0NUvohFbpIGunYIr/9+S28uaueU4bm8W8fns4nS8dpn3JRoYukg9b2CEvW7OWel7exsfoQxcMG8t2PTOfjs0t0/hV5hwpdJIXVHW7hwRW7+dVrO9h/qIXTRg/i+x87hw/PLNYWubyHCl0kxbg7r2+r44E3drN0XTWt4Qjvn1rIf1wzg4umFmqMXLqlQhdJEXvqm/nD6j08XFbJttq3GZKXzfVzx3P93PGcNnpw0PEkDajQRQLU0NTGU+uq+P2be1i+/QAApROGc+vFp3LlOWN1cWY5Lip0kT62p76ZZ8ureWb9PpZvP0A44kwuLOC2S09j4bnFjB+ZH3RESVMqdJEkC0ec9XsbeWHjfp5ZX0353kYApo4axC0XTWbB9DGcXTxUY+Ny0lToIgnm7mytOcxrW+t4taKWZdsO0NDchhnMGj+cr19+BpdOG83kokFBR5UM06tCN7MFwO1ACLjH3b/Xab7F5l8BNAGfcfdVCc4qkpKOtIUp39vI6t31rN5dzxvb69jX2AJA8bCB/MW00cw7tZB5pxZSNHhAwGklk/VY6GYWAu4ELgUqgRVmtsTd18ctdjkwNXabC9wVuxfJKM2tYSr2H2ZjdSNrKxtYU1nPhqpG2sIOwNiheZw3cQTzTi3kwikjGT8iX0Mp0md6s4U+B6hw920AZvYgsBCIL/SFwL3u7sAyMxtmZmPdvSrRgXcfaGLZtjpiWeh4q5hFbwCGEf8ees9ysWfRx+9Op8vp3a8DgywzsrNit1BW7N7ICWURyjJysrLIDsXN73iclUVOyPRmT1H1Ta3sOtDEjromKvYdYmP1ITbvO8TOA014tLspyA1xTskwbn7/ZGaUDGPm+GGMHqKzG0pwelPoxcDuuOeVvHfru6tlioGjCt3MFgGLAMaPH3+8WQFYW9nAPz289oR+NhWFsowB2VkMyM4iLydEXk4o+jwnRF4X9x3zo8tmMSA7RF5uiIGx5wNzYo/fmRaKm5ZFbiir3/8SaQ9HqDncQnXDEfY1trCv8QhVDUfYfbCJXXVN7Kx7m8Yj7e8sn2UwqbCAaacM4cMzizl99GCmjh7MpMICQln9++9SUktvCr2r/7F+Asvg7ncDdwOUlpa+Z35vXHxGES9/5eLY63WsyOMeR7+Uevfxu1Hc3w0VfRw3vYvX6s06Iu60h532SOwWjtAWdsIRpz3S8Th63x6OHLVce8RpC0doaYvQ0h7hSFuYI+0RWmL3R9rCNDS3sb8t/O78uMeRE/gbzDLeKfq8nBADc0NH/RLIy856Z1rH/LzsEANzs46adtRrHDXt3V86yfjF4e60tEdobg3T1BamuTV2a4ve3m5pp76pjfrm1uh9U+y+uY2GpjYONLVSe7jlqH9jgOwso2T4QMaPLODcccMYPyKf8SPzmTAyn4kjC7Q/uKSF3hR6JTAu7nkJsPcElkmI/Nxs8kdo5xx3py3sHGkP09IWLfjmtnfLraP8o9MiR0+LK8CWtsg7P9cY++UR/1pH2iK0hiPHnc8M8rJDZGcZWVlGKMvIMiOUBSGLTos+N7Is+ksyHPH33twJh6P3Hb8AO5dxd/Jyshg2MJdh+TkMy89hUmEBM/OHMWpIHmOG5DF6yABGD8lj9JA8RhbkkqWtbUlzvWnGFcBUM5sE7AGuBa7vtMwS4G9j4+tzgYZkjJ/Lu8yM3GyLnqApycO24Yh3Kvlo0TfHTWtpP/oXxZHY4/ZI9JNNRzlH4oo6EnHCDpGIYxYdfgplGSGLfg/R8f1EVta797mh6CeI/I5PBbnZDMwJkZ8b/bRQMCDE8Pxchg7M0Va19Ds9Frq7t5vZ3wJLie62+HN3Lzezz8fmLwaeJLrLYgXR3RY/m7zI0tdCWUbBgGwKBuiTkUgq69U71N2fJFra8dMWxz124NbERhMRkeOhEyqLiGQIFbqISIZQoYuIZAgVuohIhlChi4hkCBW6iEiGUKGLiGQI894eR53oFZvVADtP8McLgdoExkkGZUwMZUwMZTx5qZJvgrsXdTUjsEI/GWZW5u6lQec4FmVMDGVMDGU8eameDzTkIiKSMVToIiIZIl0L/e6gA/SCMiaGMiaGMp68VM+XnmPoIiLyXum6hS4iIp2o0EVEMkTaFbqZLTCzTWZWYWZfCzoPgJmNM7MXzWyDmZWb2Zdi00eY2bNmtiV2PzzgnCEze9PMHk/RfMPM7GEz2xj7u7wgBTP+fezfeJ2ZPWBmeUFnNLOfm9l+M1sXN63bTGb29dj7Z5OZXRZgxv+I/VuvNbPfm9mwVMsYN+8fzczNrDDIjD1Jq0I3sxBwJ3A5MA24zsymBZsKgHbgNnc/EzgfuDWW62vA8+4+FXg+9jxIXwI2xD1PtXy3A0+7+xnADKJZUyajmRUDXwRK3X060St4XZsCGX8JLOg0rctMsf+X1wJnxX7mv2PvqyAyPgtMd/dzgM3A11MwI2Y2DrgU2BU3LaiMx5RWhQ7MASrcfZu7twIPAgsDzoS7V7n7qtjjQ0SLqJhotl/FFvsV8OFAAgJmVgJcCdwTNzmV8g0BLgJ+BuDure5eTwpljMkGBppZNpBP9GLogWZ095eAA50md5dpIfCgu7e4+3ail42cE0RGd3/G3dtjT5cRvbh8SmWM+S/gK0SvZd4hkIw9SbdCLwZ2xz2vjE1LGWY2EZgJLAdGd1wsO3Y/KsBoPyL6nzISNy2V8k0GaoBfxIaF7jGzglTK6O57gB8Q3VKrInox9GdSKWOc7jKl6nvor4GnYo9TJqOZXQ3scfc1nWalTMZ46Vbo1sW0lNnv0swGAY8AX3b3xqDzdDCzq4D97r4y6CzHkA3MAu5y95nA2wQ/BHSU2Dj0QmAScApQYGY3BpvquKXce8jMvkl02PL+jkldLNbnGc0sH/gm8K2uZncxLfAuSrdCrwTGxT0vIfqRN3BmlkO0zO9390djk/eZ2djY/LHA/oDizQOuNrMdRIepPmhm96VQPoj+21a6+/LY84eJFnwqZbwE2O7uNe7eBjwKXJhiGTt0lyml3kNm9mngKuAGf/egmFTJOIXoL+81sfdOCbDKzMaQOhmPkm6FvgKYamaTzCyX6JcSSwLOhJkZ0bHfDe7+w7hZS4BPxx5/GvhDX2cDcPevu3uJu08k+nf2grvfmCr5ANy9GthtZqfHJn0IWE8KZSQ61HK+meXH/s0/RPT7klTK2KG7TEuAa81sgJlNAqYCbwSQDzNbAHwVuNrdm+JmpURGd3/L3Ue5+8TYe6cSmBX7v5oSGd/D3dPqBlxB9BvxrcA3g84Ty/Q+oh+31gKrY7crgJFE9zDYErsfkQJZ5wOPxx6nVD7gXKAs9vf4GDA8BTN+B9gIrAN+DQwIOiPwANEx/TaipXPTsTIRHUbYCmwCLg8wYwXRceiO98ziVMvYaf4OoDDIjD3ddOi/iEiGSLchFxER6YYKXUQkQ6jQRUQyhApdRCRDqNBFRDKECl1EJEOo0EVEMsT/AgUWk63l0H8fAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(y_ex)" ] }, { "cell_type": "markdown", "id": "c1ced6a0", "metadata": {}, "source": [ "Note that an ARIMA model is not stationary, but can show significant growth and trends due to the integration.\n", "It is therefore often used to fit non-stationary time-series. In practice one usually differentiates the data and uses an ARMA model on the differetiated data though." ] }, { "cell_type": "code", "execution_count": null, "id": "d53db19f", "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.8.10" } }, "nbformat": 4, "nbformat_minor": 5 }