ARMA init AR has same result as statsmodels

This commit is contained in:
ritchie46 committed 2018-09-21 16:20:45 +02:00
1 parent 01b769d5b6
commit 4f61176e52
1 file changed
+548
+548
View File
@@ -0,0 +1,548 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy import stats\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import statsmodels.api as sm\n",
"from sklearn import linear_model\n",
"from sklearn import metrics\n",
"from statsmodels.graphics.api import qqplot"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"::\n",
"\n",
" Number of Observations - 309 (Annual 1700 - 2008)\n",
" Number of Variables - 1\n",
" Variable name definitions::\n",
"\n",
" SUNACTIVITY - Number of sunspots for each year\n",
"\n",
" The data file contains a 'YEAR' variable that is not returned by load.\n",
"\n"
]
}
],
"source": [
"print(sm.datasets.sunspots.NOTE)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"data = sm.datasets.sunspots.load_pandas().data"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"data.index = pd.Index(sm.tsa.datetools.dates_from_range('1700', '2008'))\n",
"del data['YEAR']"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7fc8131a0c18>"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAHVCAYAAADl4K3UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvX20LWdd5/l96mXvffZ5uW+5QCCBxAi0BOEyRrtxwAk20rSLRjFrNHGcgR5dDKMMY7drOfZSl7imQ+u4GGcUB4UGgziERhAGHWxMHAwgKCQkxAAJEEjIhdzk5t6b+3Ledr0880fVU7t27ara9TxVdU7tOt/PWlk5d59z6tQ5e++nfvV9vr/vT0gpQQghhBBCCMnH2u8TIIQQQgghpMuwYCaEEEIIIaQEFsyEEEIIIYSUwIKZEEIIIYSQElgwE0IIIYQQUgILZkIIIYQQQkpgwUwIIYQQQkgJLJgJIYQQQggpgQUzIYQQQgghJTj7fQIAcNlll8mrrrpqv0+DEEIIIYT0nLvuuusJKeVxne/pRMF81VVX4c4779zv0yCEEEIIIT1HCPGw7vfQkkEIIYQQQkgJLJgJIYQQQggpgQUzIYQQQgghJXTCw0wIIYQQssx4noeTJ09iZ2dnv0+FxIxGI1xxxRVwXbf2sVgwE0IIIYTU5OTJk1hfX8dVV10FIcR+n86BR0qJM2fO4OTJk7j66qtrH4+WDEIIIYSQmuzs7ODYsWMsljuCEALHjh1rTPFnwUwIIYQQ0gAslrtFk88HC2ZCCCGEEEJKYMFMCCGEENITbr75Zlx77bV4wQtegBMnTuAf/uEfcNVVV+GJJ55IvuZv//Zv8apXvQoAcMstt8CyLNx7773J55///OfjoYceSv599913QwiBj3/84zM/69SpU7jxxhtxzTXX4HnPex5+9Ed/FF/84hdx4sQJnDhxAkePHsXVV1+NEydO4OUvfzkeeughPP/5z8fm5iaOHTuG8+fPzxzvx3/8x/GBD3wAt9xyC974xjfi5ptvTo5l23by8e/93u/hxS9+MaSUAIAgCHDixAl85jOfafrPmcCCmRBCCCGkB3z2s5/FX/7lX+ILX/gC7r33Xtx+++248sorF37fFVdcgZtvvrnw87feeite8pKX4NZbb00ek1LiNa95Da6//no8+OCD+PKXv4y3vOUtuHDhAu655x7cc889ePWrX43f+Z3fwT333IPbb789+d7V1VW84hWvwEc+8pHksfPnz+PTn/50UsgDwK/+6q8mx1pZWUk+ftOb3oRnPetZeNe73gUA+P3f/318//d/P37wB39Q6++lA1MyCCGEEEIa5Df/4kv48ncuNHrM5z19A7/xr64t/ZpHH30Ul112GYbDIQDgsssuq3TsV73qVfjkJz+JBx54AM997nNnPielxAc/+EHcdttteOlLX4qdnR2MRiN84hOfgOu6eMMb3pB87YkTJyr/PjfddBPe/va347WvfS0A4MMf/jBe+cpXYjweV/r+3/3d38VLXvISvPjFL8bb3vY2fO5zn6v8s02gwkwIIYQQ0gNe8YpX4JFHHsFznvMc/PzP/zzuuOOOSt9nWRZ++Zd/GW95y1vmPvd3f/d3uPrqq3HNNdfg+uuvx8c+9jEAwH333Yfv+77vMz7XV77ylbjrrrtw5swZAMD73/9+3HTTTZW///LLL8cv/uIv4sUvfjF+7dd+DUePHjU+lypQYSaEEEIIaZBFSnBbrK2t4a677sKnPvUpfOITn8BP/dRP4bd+67dy0yKyj/30T/80br75Znzzm9+cefzWW2/FjTfeCAC48cYb8d73vhc/8RM/UftcB4MBXv3qV+ODH/wgbrjhBtxzzz14xSteoXWMX/iFX8Cv/Mqv4HWve13t81kEC2ZCCCGEkJ5g2zauv/56XH/99fje7/1evOc978GxY8dw7ty5xKJx9uzZObuG4zj4pV/6Jfz2b/928lgQBPjQhz6Ej370o7j55puTYSAXL17Etddeiw9+8IO1zvWmm27Cv//3/x5SSvzYj/2Y9kQ+y7L2LMqPlgxCCCGEkB7wwAMP4Gtf+1ry73vuuQfPetazcP311+O9730vgKgI/tM//VO87GUvm/v+173udbj99ttx+vRpAMDtt9+OF77whXjkkUfw0EMP4eGHH8YNN9yAj3zkI/jhH/5h7O7u4p3vfGfy/Z///Ocr20AA4GUvexm+9rWv4Q/+4A+07Bj7AQtmQgghhJAecOnSJbz2ta/F8573PLzgBS/Al7/8Zbz5zW/Gr//6r+PrX/86XvjCF+JFL3oRvvu7vxs/8zM/M/f9g8EAb3rTm/D4448DiOwYr3nNa2a+5oYbbsD73vc+CCHw4Q9/GLfddhuuueYaXHvttXjzm9+Mpz/96ZXP17Is3HDDDThz5gx+6Id+qN4v3zJCZdjtJ9ddd52888479/s0CCGEEEKM+MpXvoLv+Z7v2e/TIBnynhchxF1Syut0jkOFmRBCCCG9ZdcP8E/fcjs+9o+P7vepkCWGBTMhhBBCesvmboDHLuzi019/YvEXE1IAC2ZCCCGE9BY/CAEAD5y62PrP6oLNlUxp8vlgwUwIIYSQ3uKFUdH0wKmLrRa0o9EIZ86cYdHcEVQE3mg0auR4zGEmhBBCSG9RCvOlXR8nz23jyqPVRi/rcsUVV+DkyZNJJBvZf0ajEa644opGjsWCmRBCCCG9xQumiu8Dpy62VjC7rourr766lWOT/YeWDEIIIYT0Fj8Mk48feKx9HzPpJyyYCSGEENJb/JTCfP8eNP6RfsKCmRBCCCG9xYs9zAPHwgOnLuzz2ZBlhQUzIYQQQnqLH6dkfM/T1vGN05uY+OGC7yBkHhbMhBBCCOktSmG+9hmH4IcSD56+tM9nRJYRFsyEEEII6S3Kw/z8px8CsDcDTEj/YMFMCCGEkN6iUjKuPLoCAHji0u5+ng5ZUlgwE0IIIaS3qBzmtWE0emKXHmZiAAtmQgghhPQWZckYDxxYAtjxgn0+I7KMsGAmhBBCSG9RlgzHFhi5NrYnLJiJPiyYCSGEENJblCXDtSysuDZ2fBbMRB8WzIQQQgjpLX6QVZjpYSb6sGAmhBBCSG/x4sElUcFsUWEmRrBgJoQQQkhvUQqza1kYuTZ26GEmBrBgJoQQQkhv8TKWDCrMxAQWzIQQQgjpLUnTnx01/TElg5jAgpkQQgghvUXlMDtW7GH22PRH9GHBTAghhJDeonKYbYuWDGIOC2ZCCCGE9BYvkHBtASEEm/6IMSyYCSGEENJb/CCEY0XlTjS4hJYMog8LZkIIIYT0Fj+UcGwBABi5Fpv+iBEsmAkhhBDSW7wghGtH5Y7yMEsp9/msyLLBgpkQQgghvcUPJBxLKcw2pAR2acsgmiwsmIUQ7xZCPC6EuC/12H8SQtwT//eQEOKe+PGrhBDbqc/9YZsnTwghhBBShhfOKswAsMtoOaKJU+FrbgHwNgB/oh6QUv6U+lgI8VYA51Nf/6CU8kRTJ0gIIYQQYoofTD3MK3HBvOMHOAR3P0+LLBkLC2Yp5SeFEFflfU4IIQD8JIAfbva0CCGEEELq44dhypIRKc1s/CO61PUwvxTAY1LKr6Ueu1oIcbcQ4g4hxEuLvlEI8XohxJ1CiDtPnz5d8zQIIYQQQuaJcpinsXIAOLyEaFO3YL4JwK2pfz8K4JlSyhcB+LcA3ieE2Mj7RinlO6SU10kprzt+/HjN0yCEEEIImccPwlSsXFQwU2EmuhgXzEIIB8BPAPhP6jEp5a6U8kz88V0AHgTwnLonSQghhBBigh/KZHDJMLZk7LDpj2hSR2F+OYD7pZQn1QNCiONCCDv++LsAPBvAN+qdIiGEEEKIGVEOc6bpz6PCTPSoEit3K4DPAniuEOKkEOJn40/diFk7BgD8EIB7hRBfBPBBAG+QUp5t8oQJIYQQQqoS5TDPxsqxYCa6VEnJuKng8dflPPYhAB+qf1qEEEIIIfXxQomVwXysHCE6cNIfIYQQQnqLH8wPLtme0MNM9GDBTAghhJDekh6NTQ8zMYUFMyGEEEJ6S3o0tkrJ2GbBTDRhwUwIIYSQ3pIejT10LAgB7LJgJpqwYCaEEEJIb/GDMEnJEEJg5NhUmIk2LJgJIYQQ0lu8UCY5zAAwci0OLiHasGAmhBBCSG9Jj8YGosY/KsxEFxbMhBBCCOkt6cElQBQtx5QMogsLZkIIIYT0liglI23JsGnJINqwYCaEEEJIb/EDmcTKAcrDTIWZ6MGCmRBCCCG9REoJP5RwbFoySD1YMBNCCCGkl/ihBAC4Fpv+SD1YMBNCCCGkl/hBVDBTYSZ1YcFMCCGEkF7ihVFzH5v+SF1YMBNCCCGklyQKs5UdXEKFmejBgpkQQgghvcQPIiU5bclYoSWDGMCCmRBCCCG9xFNNfxlLxrYXQEq5X6dFlhAWzIQQQgjpJYnCbM3mMIcS8AIWzKQ6LJgJIYQQ0ku8JCVjVmEGwGg5ogULZkIIIYT0Ej9JyZiNlQOAXRbMRAMWzIQQQgjpJXkpGStUmIkBLJgJIYQQ0ku8oFhhZhYz0YEFMyGEEEJ6iRqNnfYwrwyi0ofRckQHFsyEEEII6SVeXkqGQ0sG0YcFMyGEEEJ6ifIwp3OYh4klgwUzqQ4LZkIIIYT0EpWSkZ30B7BgJnqwYCaEEEJIL/FyUjJGrvIws+mPVIcFMyGEEEJ6ydSSMZ+SQQ8z0YEFMyGEEEJ6ydSSMVWYVfGsxmYTUgUWzIQQQgjpJcqS4aZSMlQDoPpc3zl1fgdnNyf7fRpLDwtmQgghhPQSpSKnFWbVAKjU577zP936Bbz5o1/a79NYepz9PgFCCCGEkDbwcgaXHDSF+cktD0KIxV9ISqHCTAghhJBeohTmGUtG/LF3QDzMQShxYdvb79NYeqgwE0IIIaSXqJSMtMJsWQKWmH6u73hhiJ1tJoLUhQozIYQQQnqJF/uLine truncated
"text/plain": [
"<Figure size 864x576 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"data.plot(figsize=(12, 8))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtUAAAHiCAYAAAAnJDDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3X2cnfVd5//XZ85kICGkw02SEkgILRgDbhswQrvVNUqroG3pz/UGVErdVOwqVWtde2N/ba3a7e7+arXKumLBIiq9cV2lSsVKjbq1IAFSbcCUNCVMCJAQMgTIkJM55/P741wnOTOZmczknDPnnJnX8/GYxznX3ZnPXLky876+1/f6XpGZSJIkSTpxfZ0uQJIkSep1hmpJkiSpSYZqSZIkqUmGakmSJKlJhmpJkiSpSYZqSZIkqUmGaknShCLi0Yh47Qlu+x0Rsa3VNUlStzJUS9IJiohNEbE/Ik6awTYZEee3s65OGP9zZeY/ZuaaTtYkSbPJUC1JJyAiVgPfASTwxo4WcxwR0T+deZKkE2eolqQT82bgHuCTwHX1mUXr9Vsbpt8SEf+3eP8PxeyvRMTzEfEjxfyfjIjtEfFMRNwRESsatr8oIr5QLHsqIt5bzD8pIn4zInYXX79ZbzGPiA0RsSsi3hURTwJ/MNG8Yt3XR8SWiBiOiH+KiFdM9MNGxKUR8eVivSci4nciYmCyn6v+/Rq2X1vsm+GI2BoRb2xY9smIuDEi/ioinouIeyPi5Sf2zyJJnWGolqQT82bgj4uv742I5cfbIDP/Q/H2lZm5ODM/HRHfDfxX4IeBs4CdwKcAIuJU4G+BvwZWAOcDdxef8cvAq4B1wCuBS4H3NXy7lwKnA+cC1080LyIuAW4Bfgo4A/g94I5JurNUgHcAZwKvBi4Hfnqyn6txw4hYAHwO+BtgGfB24I8jorF7yDXArwCnAduBX59wJ0pSlzJUS9IMRcS3Uwumn8nM+4GvAz96gh/3Y8AtmflAZh4C3gO8uuhe8nrgycz8aGa+mJnPZea9Ddt9KDP3ZOZeaoH02obPrQIfyMxDmTkyybyfBH4vM+/NzEpm3gocohbWx8jM+zPznswczcxHqQXw75zmz/gqYDHwkcwsZ+YXgb+kFqTr/iwz/zkzR6mdqKyb5mdLUlcwVEvSzF0H/E1mPl1M/wkNXUBmaAW11mkAMvN5YB9wNrCSWmA/7nbF+xUN03sz88Vx24yfdy7wzqJLxnBEDBffc8W47YiIb4qIv4yIJyPiAPBhaq3W07ECGMrM6rh6z26YfrLh/UFqIVySeoY3qkjSDETEQmpdNUpF32SAk4DBiHgl8AKwqGGTlx7nI3dTC7f1zz+FWleMx4EhxrbmTrTd1mJ6VTGvLifYZvy8IeDXM3M6XS1+F3gQuCYzn4uInwd+cBrb1WtdGRF9DcF6FfC1aW4vSV3PlmpJmpk3UetffCG1LgrrgLXAP1LrZ70F+IGIWFQMMbdx3PZPAS9rmP4T4CciYl3Rl/nDwL1FF4u/BF4aET9f3Jh4akRcVmx3O/C+iFgaEWcC7wf+aIY/y+8Db4uIy6LmlIj4/qIv93inAgeA5yPim4H/fJyfq9G91E42fikiFkTEBuANFH3HJWkuMFRL0sxcB/xBZj6WmU/Wv4DfodbP+WNAmVrIvJVa/+BGHwRuLbpb/HBm3g38v8D/Bp4AXg5cDZCZzwGvoxZAnwQeAb6r+JxfAzYD/wL8K/BAMW/aMnMztX7VvwPsp3aD4FsmWf0XqfUbf45aGP/0uOVjfq5x36dMbdjBK4Gngf8JvDkz/20m9UpSN4vMia4QSpIkSZouW6olSZKkJhmqJUmSpCYZqiVJkqQmGaolSZKkJhmqJUmSpCb15MNfzjzzzFy9enWny5AkSdIcdv/99z+dmUuns25PhurVq1ezefPmTpchSZKkOSwidk53Xbt/SJIkSU0yVEuSJElNMlRLkiRJTWprqI6IWyJiT0R8dZLlEREfj4jtEfEvEXFJO+uRJEmS2qHdLdWfBK6YYvmVwAXF1/XA77a5nhNSqSZ3P/wUH7/7Ee5++Ckq1ex0SZIkSeoibR39IzP/ISJWT7HKVcAfZmYC90TEYESclZlPtLOumahUk2tvvpctQ8OMlCssHCixbuUgt228jFJfdLo8SZIkdYFO96k+GxhqmN5VzOsam7btYcvQMAfLFRI4WK6wZWiYTdv2dLo0SZIkdYlOh+qJmnon7FsREddHxOaI2Lx37942l3XU1t0HGClXxswbKVd4aPeBWatBkiRJ3a3ToXoXsLJh+hxg90QrZuZNmbk+M9cvXTqtB9u0xEUrlrBwoDRm3sKBEheuWDJrNUiSJKm7dTpU3wG8uRgF5FXAs93Unxpgw5plrFs5SFTKkFUWFX2qN6xZ1unSJEmS1CXaeqNiRNwObADOjIhdwAeABQCZ+b+AO4HvA7YDB4GfaGc9J6LUF9y28TJe/QMbKZ+yjI++7x1sWLPMmxQlSZJ0RLtH/7jmOMsT+Jl21tAKpb5g0fAOFg3v4PK1yztdjiRJkrpMp7t/SJIkST3PUC1JkiQ1yVAtSZIkNclQLUmSJDXJUC1JkiQ1yVAtSZIkNclQLUmSJDXJUC1JkiQ1yVAtSZIkNclQLUmSJDXJUC1JkiQ1yVAtSZIkNclQLUmSJDXJUC1JkiQ1yVAtSZIkNam/0wXMdZVqsmnbHrbuPsBFK5awYc0ySn3R6bIkSZLUQm0P1RFxBfBbQAn4RGZ+ZNzyVcCtwGCxzrsz88521zUbKtXk2pvvZcvQMCPlCgsHSqxbOchtGy8zWEuSJM0hbe3+EREl4EbgSuBC4JqIuHDcau8DPpOZFwNXA/+znTXNpk3b9rBlaJiD5QoJHCxX2DI0zKZtezpdmiRJklqo3X2qLwW2Z+aOzCwDnwKuGrdOAkuK9y8Bdre5plmzdfcBRsqVMfNGyhUe2n2gQxVJkiSpHdodqs8GhhqmdxXzGn0Q+PGI2AXcCby9zTXNmotWLGHhQGnMvIUDJS5csWSSLSRJktSL2h2qJ+o4nOOmrwE+mZnnAN8H3BYRx9QVEddHxOaI2Lx37942lNp6G9YsY93KQaJShqyyqOhTvWHNsk6XJkmSpBZqd6jeBaxsmD6HY7t3bAQ+A5CZXwZOBs4c/0GZeVNmrs/M9UuXLm1Tua1V6gtu23gZSx/5HIO7vsRvX3OxNylKkiTNQe0O1fcBF0TEeRExQO1GxDvGrfMYcDlARKylFqp7oyl6Gkp9waLhHQw+fg+Xr11uoJYkSZqD2hqqM3MUuAG4C3iY2igfWyPiQxHxxmK1dwI/GRFfAW4H3pKZ47uISJIkSV2r7eNUF2NO3zlu3vsb3j8EvKbddUiSJEnt4mPKJUmSpCYZqiVJkqQmGaolSZKkJhmqJUmSpCYZqiVJkqQmGaolSZKkJhmqJUmSpCYZqiVJkqQmGaolSZKkJhmqJUmSpCa1/THlkqZWqSabtu1h6+4DXLRiCRvWLKPUF50uS5IkzYChWuqgSjW59uZ72TI0zEi5wsKBEutWDnLbxsumHawN5ZIkdZ6hWuqgTdv2sGVomIPlCgAHyxW2DA2zadseLl+7/LjbtyKUS5Kk5tmnWuqgrbsPMFIE6rqRcoWHdh+Y1vaNoTwZG8olSdLsMVRLHXTRiiUsHCiNmbdwoMSFK5ZMa/tmQ7kkSWoNQ7XUQRvWLGPdykGiUoassqjovrFhzbJpbd9sKJckSa1hqJY6qNQX3LbxMpY+8jkGd32J377m4hn1h242lEuSpNZoe6iOiCsiYltEbI+Id0+yzg9HxEMRsTUi/qTdNUndpNQXLBreweDj93D52uUzusGw2VAuSZJao62jf0RECbgReB2wC7gvIu7IzIca1rkAeA/wmszcHxE2sUkzUA/li4Z3TGvEEEmS1Hrtbqm+FNiemTsyswx8Crhq3Do/CdyYmfsBMtNhCyRJktRT2h2qzwaGGqZ3FfMafRPwTRHxpYi4JyKumOiDIuL6iNgcEZv37t3bpnIlSZKkmWt3qJ6oY2eOm+4HLgA2ANcAn4iIwWM2yrwpM9dn5vqlS5e2vFBJkiTpRLU7VO8CVjZMnwPsnmCdv8jMw5n5DWAbtZAtSZIk9YR2h+r7gAsi4ryIGACuBu4Yt86fA98FEBFnUusOsqPNdUmSJEkt09ZQnZmjwA3AXcDDwGcyc2tEfCgi3lisdhewLyIeAv4O+C+Zua+ddUmSJEmt1NYh9QAy807gznHz3t/wPoFfKL4kSZKknuMTFSVJkqQmGaolSZKkJhmqJUmSpCYZqiVJkqQmGaolSZKkJhmqJUmSpCYZqiVJkqQmtX2camk+qFSTTdv2sHX3AS5asYQNa5ZR6otOlyVJkmaJoVpqUqWaXHvzvWwZGmakXGHhQIl1Kwe5beNlBmtJkuYJu39ITdq0bQ9bhoY5WK6QwMFyhS1Dw2zatqfTpR1XpZrc/fBTfPzuR7j74aeoVLPTJUmS1JNsqe5ydivoflt3H2CkXBkzb6Rc4aHdB7h87fIOVXV8trBLktQ6huouZujpDRetWMLCgRIHG4L1woESF65Y0sGqjq+xhR3GtrB388mAJEndyO4fXayXuxXMJxvWLGPdykGiUoassqg4+dmwZlmnS5vSVC3skiRpZgzVXczQ0xtKfcFtGy9j6SOfY3DXl/jtay7uiasJ9Rb2Rr3Qwi5JUjcyVHcxQ0/vKPUFi4Z3MPj4PVy+dnnXB2ro3RZ2zS5vZpWk6bFPdRerh54vf+0Jsq+fRSctMPSoZeot7K/+gY2UT1nGR9/3Dm+E1RjN3tfhjdaS5pO2h+qIuAL4LaAEfCIzPzLJej8IfBb4tszc3O66eoGhR+1Wb2FfNLzDmxN1jGZuZvVGa0nzTVu7f0RECbgRuBK4ELgmIi6cYL1TgZ8F7m1nPb2oF7sVSJobmrmvwxutJc037e5TfSmwPTNLine truncated
"text/plain": [
"<Figure size 864x576 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(12, 8))\n",
"ax = fig.add_subplot(211)\n",
"fig = sm.graphics.tsa.plot_acf(data.values.squeeze(), lags=40, ax=ax)\n",
"ax = fig.add_subplot(212)\n",
"fig = sm.graphics.tsa.plot_pacf(data.values.squeeze(), lags=40, ax=ax)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Autocorrelation\n",
"\n",
"Let's look how autocorrelation is defined. First we make a function that computes the Pearson correlation between two sequences that are assumed to come from a Normal distribution."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"x = data.values.squeeze()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.999687011702813"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def pearson_correlation(x, y):\n",
" return np.mean((x - x.mean()) * (y - y.mean())) / (x.std() * y.std())\n",
" \n",
"pearson_correlation(x, x + np.random.randn(x.shape[0]))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The autocorrelation is the correlation between the sequence with the current last value truncated at the end `x[:-1]` and the sequence with a certain amount of lag `x[:i]`."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAD8CAYAAABnwSKjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAEzNJREFUeJzt3XGMZdddH/DvL7sYKpKQBG+C5fWybruo2SLqlKmbKi0dErtdp8hLpYTaAupICdsKTKlIKU6DUmqEFKho8gdu1W1I7YaCMWkhK7qVExy7rSqSet2YgO263pgYL7a8i0loESKum1//mOd0GM/sHufdmTcz/nyk0dxz35l7zp49mved8857t7o7AADA+b1k0R0AAICdQHAGAIABgjMAAAwQnAEAYIDgDAAAAwRnAAAYIDgDAMAAwRkAAAYIzgAAMGDvojuwkYsvvrgPHjy46G4AALDL3Xfffb/b3fsuVG/bBueDBw/m1KlTi+4GAAC7XFU9NlLPVg0AABggOAMAwADBGQAABgjOAAAwQHAGAIABgjMAAAwQnAEAYIDgDAAAAyYJzlV1pKoerqrTVXXTOo8fqKq7q+pTVfXpqnrzFO1uhuXl5SwvLy+6GwAAbDNzB+eq2pPkliTXJDmc5PqqOrym2o8kuaO7X5fkuiT/fN52AQBgK02x4nxlktPd/Wh3P5Pk9iRH19TpJC+fHX9NkicmaBcAALbM3gmucWmSx1eVzyT5i2vq/GiSj1bV9yf56iRXTdAuAABsmSlWnGudc72mfH2SW7t7f5I3J/lQVT2v7ao6VlWnqurUuXPnJugaAABMY4rgfCbJZavK+/P8rRhvT3JHknT3ryX5qiQXr71Qdx/v7qXuXtq3b98EXQMAgGlMEZzvTXKoqi6vqouy8ua/E2vq/HaSNyVJVb02K8HZkjIAADvG3MG5u59NcmOSO5M8lJVPz3igqm6uqmtn1d6Z5Huq6teT/HySt3X32u0cAACwbU3x5sB098kkJ9ece8+q4weTvGGKtgAAYBHcORAAAAYIzgAAMEBwBgCAAYIzAAAMEJwBAGCA4AwAAAMEZwAAGCA4AwDAAMEZAAAGCM4AADBAcAYAgAGCMwAADBCcAQBggOAMAAADBGeY0PLycpaXlxfdDQBgEwjOsI0I3gCwfQnOAAAwQHAGAIABgjMAAAwQnAEAYIDgDAAAAwRnAAAYIDgDAMCASYJzVR2pqoer6nRV3bRBne+oqger6oGq+rkp2gUAgK2yd94LVNWeJLckuTrJmST3VtWJ7n5wVZ1DSd6V5A3d/bmqevW87QIAwFaaYsX5yiSnu/vR7n4mye1Jjq6p8z1JbunuzyVJd5+doF0AANgyUwTnS5M8vqp8ZnZutW9I8g1V9V+r6hNVdWS9C1XVsao6VVWnzp07N0HXAABgGlME51rnXK8p701yKMlykuuTfKCqXvG8H+o+3t1L3b20b9++CboGAADTmCI4n0ly2ary/iRPrFPnI939f7r7t5I8nJUgDQAAO8IUwfneJIeq6vKquijJdUlOrKnzy0m+NUmq6uKsbN14dIK2AQBgS8wdnLv72SQ3JrkzyUNJ7ujuB6rq5qq6dlbtziRPV9WDSe5O8kPd/fS8bQMAwFaZ++PokqS7TyY5uebce1Ydd5IfnH0BAMCO486BAAAwQHAGAIABgjMAAAwQnAEAYIDgDAAAAwRnAAAYIDgDAMAAwRkAAAYIzgAAMEBwBgCAAYIzAAAMEJwBAGCA4AwAAAMEZ1hleXk5y8vLi+4GALANCc6wiwj+ALB5BGcAABggOAMAwADBeWJeKgd2Kr+/AM5PcAZgEoI3sNsJzgAAMEBwBgCAAYIzAAAMmCQ4V9WRqnq4qk5X1U3nqfeWquqqWpqiXQAA2CpzB+eq2pPkliTXJDmc5PqqOrxOvZcl+XtJPjlvmwAAsNWmWHG+Msnp7n60u59JcnuSo+vU+7EkP5nkjyZoEwAAttQUwfnSJI+vKp+ZnfuSqnpdksu6+1cmaA8AALbcFMG51jnXX3qw6iVJ3pfknRe8UNWxqjpVVafOnTs3QdcAAGAaUwTnM0kuW1Xen+SJVeWXJfnGJPdU1WeTvD7JifXeINjdx7t7qbuX9u3bN0HXAABgGlME53uTHKqqy6vqoiTXJTnx3IPd/fvdfXF3H+zug0k+keTa7j41QdsAALAl5g7O3f1skhuT3JnkoSR3dPcDVXVzVV077/UBAC7ELd/ZCnunuEh3n0xycs2592xQd3mKNgEAYCu5cyAAAAwQnAEAYIDgDAAAAwRnAAAYIDgDAMAAwRkAAAYIzgDAi57PgWaE4AywS3jiB9hcgjMAAAwQnAEAYIDgDAAAAwRnALYFe7SB7U5wBgAWzh9O7ASCMwAADBCcAQBggOAMAAADBGcAABggOAMAwADBGQAABgjOAAAwQHAGAIABgjMAAAwQnAEAYMAkwbmqjlTVw1V1uqpuWufxH6yqB6vq01V1V1V9/RTtAgDAVpk7OFfVniS3JLkmyeEk11fV4TXVPpVkqbu/KcmHk/zkvO0CAMBWmmLF+cokp7v70e5+JsntSY6urtDdd3f3H86Kn0iyf4J2AQBgy0wRnC9N8viq8pnZuY28Pcl/nKBdAADYMnsnuEatc67XrVj1XUmWkvzVDR4/luRYkhw4cGCCrgEAwDSmWHE+k+SyVeX9SZ5YW6mqrkry7iTXdvcX1rtQdx/v7qXuXtq3b98EXQMAgGlMEZzvTXKoqi6vqouSXJfkxOoKVfW6JP8yK6H57ARtAgDAlpo7OHf3s0luTHJnkoeS3NHdD1TVzVV17azaP03y0iS/WFX3V9WJDS4HAADb0hR7nNPdJ5OcXHPuPauOr5qiHQAAWBR3DgS+ZHl5OcvLy4vuBsCO5Hfo7ic4s6v4pQUAbBbBGQAABgjOAAAwQHAGAIABgjMAAAwQnAEAYIDgDAAAAwRnAAAYIDgDAMAAwRlgm3ADH4DtTXAGAObmDz9eDARnAAAYIDgDAMAAwRkAAAYIzgAAMEBwBgCAAYIzAAAMEJwB2BV8HBqw2QRnAAAYIDhvM1ZMAAC2J8EZAAAGTBKcq+pIVT1cVaer6qZ1Hv/KqvqF2eOfrKqDU7QLAABbZe7gXFV7ktyS5Jokh5NcX1WH11R7e5LPdfefTvK+JD8xb7sAALCVplhxvjLJ6e5+tLufSXJ7kqNr6hxNctvs+MNJ3lRVNUHbAACwJaYIzpcmeXxV+czs3Lp1uvvZJL+f5GsnaBsAALZEdfd8F6h6a5K/3t3vmJW/O8mV3f39q+o8MKtzZlb+zKzO02uudSzJsSQ5cODANz/22GNz9e3F6LlP5LjnnntelD/PYi36/3+n/zyLtej//53+88xv0f+HO/3n51FV93X30oXqTbHifCbJZavK+5M8sVGdqtqb5GuS/N7aC3X38e5e6u6lffv2TdA1AACYxhTB+d4kh6rq8qq6KMl1SU6sqXMiyQ2z47ck+XjPu9QNAABbaO+8F+juZ6vqxiR3JtmT5IPd/UBV3ZzkVHefSPIzST5UVaezstJ83bztAgDAVpo7OCdJd59McnLNufesOv6jJG+doi0AAFgEdw4EAIABgjMAAAwQnAEAYIDgDAAAAwRnAAAYMMmnagAAMB93bdz+rDgDAMAAwRkAAAYIzgAAMEBwBgCAAYIzAAAMEJwBAGCA4AwAAAMEZwAAGCA4AwDAAHcOBCbjrlcA7GaCMwDgD18YYKsGAAAMsOIMMBErdgC7mxVnAAAYYMWZbcWKHQCwXVlxBgCAAYIzAAAMmCs4V9WrqupjVfXI7Psr16lzRVX9WlU9UFWfrqq/NU+bAACwCPOuON+U5K7uPpTkrll5rT9M8re7+88mOZLk/VX1ijnbBQCALTXvmwOPJlmeHd+W5J4kP7y6Qnf/z1XHT1TV2ST7knx+zrYBYDLenAxcyLwrzq/p7ieTZPb91eerXFVXJrkoyWfmbBcAALbUBVecq+pXk3zdOg+9+4U0VFWXJPlQkhu6+4sb1DmW5FiSHDhw4IVcHgAANtUFg3N3X7XRY1X1VFVd0t1PzoLx2Q3qvTzJf0jyI939ifO0dTzJ8SRZWlrqC/UNAIDdYSdsl5p3q8aJJDfMjm9I8pG1FarqoiS/lOTfdPcvztkeAAAsxLzB+b1Jrq6qR5JcPSunqpaq6gOzOt+R5FuSvK2q7p99XTFnuwAAsKXm+lSN7n46yZvWOX8qyTtmxz+b5GfnaQcAABbNnQMBAGCA4AwAAAMEZwAAGCA4AwDAAMEZAAAGzPWpGgAAbA874QYiO50VZwAAGGDFeZfx1yYAwOaw4gwAAAMEZwAAGCA4AwDAAMEZAAAGCM4AADDAp2oAwC7gU5Vg8wnOADOCBwDnIzgD24bgCsB2Zo8zAAAMsOIMAMDcXgyvGlpxBgCAAYIzAAAMsFUDACbwYniZmt3NHL4wK84AADBgruBcVa+qqo9V1SOz7688T92XV9XvVNVPz9MmAAAswrwrzjcluau7DyW5a1beyI8l+U9ztsc2d88993ipBwDYleYNzkeT3DY7vi3Jt69Xqaq+Oclrknx0zvYAAGAh5n1z4Gu6+8kk6e4nq+rVaytU1UuS/FSS707ypjnbA9iQVzsA2EwXDM5V9atJvm6dh9492Mb3JjnZ3Y9X1YXaOpbkWJIcOHBg8PIAALD5Lhicu/uqjR6rqqeq6pLZavMlSc6uU+0vJfkrVfW9SV6a5KKq+oPuft5+6O4+nuR4kiwtLfXoPwIAADbbvFs1TiS5Icl7Z98/srZCd3/nc8dV9bYkS+uFZgB4MbPVCLa/ed8c+N4kV1fVI0munpVTVUtV9YF5OwcAANvFXCvO3f101nnDX3efSvKOdc7fmuTWedoEAIBFcMtt/hgvFQIArM8ttwEAYIDgDAAAAwRnAAAYIDgDAMAAwRkAAAYIzgAAMEBwBgCAAYIzAAALine truncated
"text/plain": [
"<Figure size 864x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def plot(x, lag=40):\n",
" plt.figure(figsize=(12, 4))\n",
" corr = [pearson_correlation(x[:-i], x[i:]) for i in range(1, lag)]\n",
" plt.vlines(np.arange(1, lag), 0, corr)\n",
" plt.show()\n",
"plot(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Moveing average model MA(q)\n",
"\n",
"$$ X_t = \\mu + \\epsilon_t + \\theta_1 \\epsilon_{t-1} ... + \\theta_q \\epsilon_{q} $$\n",
"\n",
"where $\\epsilon$ is white noise:\n",
"\n",
"$$\\epsilon \\sim N(0, \\alpha)$$\n",
"\n",
"$\\alpha$ is standard deviation of the sequence?\n"
]
},
{
"cell_type": "code",
"execution_count": 249,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"54.39517901679693\n",
"4833.691033366038\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAHVCAYAAAAHPLatAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXuwLdld3/dd/dh7n9d9zr2j0cxIM0JTBEk8EiQLlxOCC2MsqmyRCqYgZcAOFaiEVB7lJIa4EuMkuCg7wYEKJoiHEQ4hgHEEiQVBFshCb80ISTMajTR33vfOzH2fe1770Y+VP7pXd+/e3b17re6zd++9v58q1dWce26fPnv3Xv3r7/r+vj8hpQQhhBBCCCFkFmvZJ0AIIYQQQkhXYbFMCCGEEEJICSyWCSGEEEIIKYHFMiGEEEIIISWwWCaEEEIIIaQEFsuEEEIIIYSUwGKZEEIIIYSQElgsE0IIIYQQUgKLZUIIIYQQQkpw6n6jEOJhAL8O4A0AQgDvk1L+rBDiAoDfAvAIgBcBfK+U8q4QQgD4WQDfBeAEwN+UUn6u6mfcd9998pFHHjH4NQghhBBCCKnPE088cUtKeWne99UulgH4AP62lPJzQog9AE8IIT4E4G8C+LCU8qeFED8O4McB/B0A7wHwWPy/dwP4hfjPUh555BE8/vjjGqdECCGEEEKIPkKIl+p8X20bhpTyNaUMSykPAXwZwIMA3gvg/fG3vR/Ad8f//70Afl1GfArAOSHEA3V/HiGEEEIIIcvGyLMshHgEwL8J4NMA7pdSvgZEBTWAy/G3PQjglcw/uxp/LX+sHxFCPC6EePzmzZsmp0MIIYQQQsipoF0sCyF2AfwugP9CSnlQ9a0FX5MzX5DyfVLKd0op33np0lzbCCGEEEIIIQtDq1gWQriICuXfkFL+i/jL15W9Iv7zRvz1qwAezvzzhwC82ux0CSGEEEIIWRy1i+U43eJXAHxZSvkzmb/6fQA/FP//HwLwe5mv/6CI+BYA95RdgxBCCCGEkFVAJw3jLwD4AQBPCiE+H3/tvwXw0wB+WwjxwwBeBvDX47/7IKLYuCuIouP+VitnTAghhBBCyIKoXSxLKT+GYh8yAHx7wfdLAD9meF6EEEIIIYQsHU7wI4QQQgghpAQWy4QQQgghhJTAYpkQQgghhJASWCwTQgghhBBSAotlQgghhBBCSmCxTAghhBBCSAkslgkhhBBCCCmBxTIhhBBCCCElsFgmhBBCCCGkBBbLhBBCCCGElMBimRBCCCEbxw/+6mfwMx/66rJPg6wAzrJPgBBCCCFk0Tx34wg9Wyz7NMgKQGWZEEIIIRuHH4a4N/SWfRpkBWCxTAghhJCNIwgli2VSCxbLhBBCCNk4/FBi/4TFMpkPPcuEEEII2TiCUOLehMUymQ+VZUIIIYRsHEEoMfZDjLxg2adCOg6LZUIIIYRsHH4oAYC+ZTIXFsuEEEII2TgCFsukJiyWCSGEELJRSClZLJPasFgmhBBCyEahCmUATMQgc2GxTAghhJCNws8Uy1SWyTxYLBNCCCFkowhYLBMNWCwTQgghZKOYUpZPJks8E7IKsFgmhBBCyEZBZZnowGKZEEIIIRuFH4bJ/2exTObBYpkQQgghG8VUGgaLZTIHFsuEEEII2Sj8YNaG8bP/6ln8q6evL+uUSIdhsUwIIYSQjSKU08WyH4T4+T+5gj946vUlnhXpKiyWCSGEELJRqDSMvb6DeyceXr5zgkkQIsh4mQlRsFgmhBBCyEahPMsXdnu4N/Tw1etH0ddl1b8imwqLZUIIIYRsFMqzfHGnBz+U+OLVfQCgskwKYbFMCCGEkI0iUZZ3+gCAx1+8C2C68Y8QBYtlQgghhGwUKmf54k4PAPCFRFlmsUxmYbFMCCGEkI1CFcUXd6NieexHxbPPYpkUwGKZkAUShhJScjEmhJBl4ic2jN7U16kskyJYLBOyQP7y//pR/OrHX1z2aRBCyEaTV5YBoO9YU2OwCVGwWCZkgbxy5wQv3jpe9mkQQshGo5Tli3GDHwA8dv8ulWVSCItlQhZIEEoMvWDZp0EIIRuNiog7s+XCtgTuP9PH+e0ePcukEBbLhCwIKSV8FsuEELJ0VEScYwmc23Lx2OU9OJagskwKYbFMyIJQi/BowmKZEEKWiVqPbUvg+//cm/C973oYtmWxWCaFOMs+AUI2BbW9N/JZLBNCyDIJZKos/1ff+bUAgH/5xVdZLJNCqCwTsiDUIjykskwIIUslqywrHMuiZ5kUwmKZkAWhFuGhx2giQghZJqlnOS2DbHqWSQkslglZEIlnmQ1+hBCyVBJl2c4qy4I5y6QQFsuELAi1CNOGQQghy0Xt9DkZG4ZtCQQBlWUyC4tlQhZEwAY/QgjpBCpnecqzbAt6lkkhLJYJWRDKI0dlmRBClkupssximRTAYpmQBaEW57EfIuSCTAghS4NpGEQHFsuELIgg0zhCKwYhhCyPVFmeTsOgkEGKYLFMyILIKha0YhBCyPIoUpZti55lUgyLZUIWhJ/psh75jCcihJBlkeYs07NM5sNimZAFEVBZJoSQThCEIYQALIs5y2Q+LJYJWRDZ7T0OJiGEkOXhhxK2EFNfsy2BUIK+ZTIDi2VCFsSUssximRBClkYg5ZRfGUgtGYFksUymYbFMyILIbu/RhkEIIcsjCOSUXxkA7DgZg75lkofFMiELYqrBj8oyIYQsDT8sV5aZiEHysFgmZEHQhkEIId0gCCUce7oEUsVzELBYJtOwWCZkQbDBjxBCukGhsmwrZZmJGGQaFsuELIiAnmVCCOkEQRgWeJbZ4EeKYbFMyIKYUpY5lIQQQpZGkbKsouTY4EfysFgmZEFwKAkhhHSDICxKw4htGPQskxwslglZEEzDIISQblDlWaayTPKwWCZkQUzlLLNYJoSQpRHlLOfTMKL/ZnQcycNimZAFoRZg2xK0YRBCyBKpylmmskzysFgmZEGoBXi377DBjxBClkgQhontQpF4lhkdR3KwWCZkQSjP8m7fobJMCCFLxA8lLEFlmdSDxTIhC2JKWaZnmRBClkYoK9IwWCyTHCyWCVkQagHeHThs8COEkCXiB0We5agk6qqyfDDyIDkwZSmwWCZkQagJfrRhEELIcglCOeNZVuEYXSyWD0ce3v1TH8YfPPX6sk9lI2GxTMiC8Kca/FgsE0LIsojSMKZLoC4rywcjH0MvwNOvHiz7VDYSFsuELAg/kLAEsNWzMaKyTAghS6Nygl8Hi+UgbhB/dX+45DPZTFgsE7Ig/DAKwd9ybXqWCSFkiVTnLHcvOs6Lz+nVeyyWlwGLZUIWhMr13OqxWCaEkGUShGG5shx0UFkOlbI8WvKZbCYslglZEErJGLg2Rl7IrmZCCFkShcqy3d2cZS+IlOXX7g0Rtnh+Ukp84sqtTv7OXYLFMiELQnnktlwbADDmFD9CCFkKRZ5lp8ue5ficvEDi1tG4teP+0dPX8R/88qfxqedvt3bMdYTFMiELQnVfD9zoY8f4OEIIWQ5RzvJ0CWR3OA3Dy1hDrrXY5Pc7j18FANw5nrR2zHWkdrEshPhVIcQNIcRTma/9pBDimhDi8/H/vivzdz8hhLgihPiKEOI72z5xQlaNIJhWlulbJoSQ5bCqyjLQnm/51tEYH/nKDQDAycRv5Zjrio6y/GsA/krB1/+xlPKb4v99EACEEG8D8H0A3h7/m38ihLCbniwhq4zyyG31WCwTQsgy8UMJ2y5u8OtiGoYfpOfUVnzc733+1eTB4HjM+1EVtYtlKeVHAdyp+e3vBfB/SSnHUsoXAFwB8OcMzo+QtcGP0zAGsbI8YrFMCCFLIQhD2KKsWF7GGVWTVbvbio/7509cxdc9cAYAleV5tOFZ/k+FEF+MbRrn4689COCVzPdcjb82gxDiR4QQjwshHr9582YLp0NIN0mUZRbLhBCyVIrSMDqtLIftKsvP3zzCl187wPe962G4tsAxe2gqaVos/wKArwHwTQBeA/C/xF8XBd9baAKSUr5PSvlOKeU7L1261PB0COkuQSDhWlaiLA8n3VuQCSFkEwhXzLOssp/v2+214lm+exI19D163w62ew5OxlSWq2hULEspr0spAyllCOCXkFotrgJ4OPOtDwF4tcnPImTVySvL9CwTQshyqPYsd7BYjs/pTRe2W1GWVXSpa1vY6dlUlufQqFgWQjyQ+c9/D4BKyvh9AN8nhOgLIR4F8BiAzzT5WYSsOukEvzg6jsUyIYQsheI0jGht7qSynCmWbx9PGtv4VBRdz7Gw3XfoWZ6DU/cbhRC/CeDbANwnhLgK4O8B+DYhxDchsli8COBHAUBK+SUhxG8DeBqAD+DHpJSsDMhGk53gB9CzTAghy0BKmeTeZ+m0shx3Hb7p4g6AyLf8lku7xsebxMpyTynLTMOopHaxLKX8/oIv/0rF9/8UgJ8yOSlC1pH8BD8Wy4QQsnhULVzqWQ46WCxnlGUgylpupVh2rMizTGW5Ek7wI2RB5JVlTvAjhJDFo5Il8mkYliUgREfTMOIC/s0Xo2L52v5Jo+N5QVos7/SpLM+DxTIhC8IPQjjZNAwqy4QQsnCUzSKvLKuvddGzrAr4h89vw7EEXrrdrFieJA1+Ajv0LM+FxTIhCyKIlWXbEug5FotlQghZAqoYzivLAGAJ0UnPsmrIG7gWHr6wjRdvHzc63iSYtmEwDaMaFsuELAg/lHDjqKKebcHzu7cgE32evHoPv/XZl5d9GoSQmgRBtbLcxWI5yBT4b764jRdvpcry4y/ewdjXK3ZnG/yoLFfBYpmQBRFkJkbZluikL47o839+5iX8gw8+s+zTWDrf/fMfxz/9+AvLPg1C5pIoy/ZsCWR31Iahzsm1LTxycQcv3j6GlBJXbhzhe/73T+JDT1/XOt6Ustx3cDIJEHbw9+4KLJYJWRB+KJMcT9cW8LgwrQVH4yCJddpknr1Line truncated
"text/plain": [
"<Figure size 864x576 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def init_eps(order, y, sigma=None):\n",
" sigma = y.std() / 3 if sigma is None else sigma\n",
" return np.random.normal(0, sigma, size=(y.shape[0]))\n",
"\n",
"\n",
"def create_ma_features(eps, y, order):\n",
" eps = np.array([eps[-(i + order):][:order][::-1] for i in range(y.shape[0])])\n",
" y = y[:-order]\n",
" eps = eps[:-order]\n",
" eps = np.stack(eps)\n",
" return eps, y\n",
"\n",
" \n",
"\n",
"def ma_model(y, order=2, plot=False, intercept=False, eps=None):\n",
" \"\"\"\n",
" When combined with AR model intercept should be false as that model provides one.\n",
" \"\"\"\n",
" y = y.copy()\n",
" mu = 0 #y.mean()\n",
" \n",
" # create a sequence of random error values.\n",
" eps = init_eps(order, y) if eps is None else eps\n",
" \n",
" # create features by shifting the window of `order` size by one step.\n",
" eps, y = create_ma_features(eps, y, order)\n",
"\n",
" # move mu and epsilon_1 to the left side of the equation.\n",
"\n",
" # This should be added to the final regression output\n",
" add = eps[:, 0] + mu\n",
"\n",
" y -= add\n",
" eps = np.delete(eps, 0, axis=1)\n",
" m = linear_model.LinearRegression(fit_intercept=intercept)\n",
" m.fit(eps, y)\n",
" pred = m.predict(eps) + add\n",
" \n",
" if plot:\n",
" print(metrics.mean_absolute_error(y, pred))\n",
" print(metrics.mean_squared_error(y, pred))\n",
"\n",
" plt.figure(figsize=(12, 8))\n",
" plt.plot(np.arange(0, y.shape[0]), y)\n",
" plt.plot(np.arange(0, y.shape[0]), pred)\n",
"\n",
" plt.show()\n",
" return pred, m\n",
" \n",
"_ = ma_model(x, 2, True, intercept=False)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Autoregressive model\n",
"\n",
"$$ X_t = c + \\epsilon_t + \\theta_1 X_{t-1} ... + \\theta_p X_{q} $$\n",
"\n",
"where $\\epsilon$ is white noise:\n",
"\n",
"$$\\epsilon \\sim N(0, \\alpha)$$\n",
"\n",
"$\\alpha$ is standard deviation of the sequence?\n",
"\n",
"$$X_t = c + \\epsilon_t + \\sum_{i=1}^{p}{\\phi_i X_{t - i}} + \\sum_{i = 1}^q{\\theta_i \\epsilon_{t-i}} $$"
]
},
{
"cell_type": "code",
"execution_count": 355,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"16.94434518547299 [-0.13025039 -0.50994881 1.30172139]\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fc7f8e37940>]"
]
},
"execution_count": 355,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD9CAYAAABOQep/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXmQJOl53vf78q6jq4+Z2QO7i8W1AEiBJEiuKIRkkrJhyhCtECU7JJEhS7TFIEiJCpMhh63DEZYiZIauoGghRJMGCYqEREJU8BAZNHiAMEgQIEBgASwWi3MPYLGzO7tz9VVVeefnP77MOnqquyuzqqdrat5fxMRMZ1fnZFVXffnkk8/7vkprjSAIgiAIgiAIY6zzPgBBEARBEARBWDVEJAuCIAiCIAjCEUQkC4IgCIIgCMIRRCQLgiAIgiAIwhFEJAuCIAiCIAjCEUQkC4IgCIIgCMIRThXJSqmHlFIfUEp9Xin1WaXUD5fbd5RS71NKPVX+vV1uV0qpdyilnlZKPaGU+qazfhKCIAiCIAiCsEzmcZIz4H/RWn8N8Bbgh5RSXwv8A+D9WutHgPeXXwP8eeCR8s/bgZ9c+lELgiAIgiAIwhlyqkjWWl/RWn+y/Pch8HngAeC7gJ8vH/bzwF8q//1dwLu14aPAllLq/qUfuSAIgiAIgiCcEbUyyUqpVwHfCPwxcK/W+goYIQ3cUz7sAeD5iR+7XG4TBEEQBEEQhDsCZ94HKqW6wK8AP6K1PlBKHfvQGdtumX2tlHo7Jo5Bp9P55je+8Y3zHoogCIIgCIIgNOITn/jEda31pdMeN5dIVkq5GIH8C1rrXy03v6yUul9rfaWMU1wtt18GHpr48QeBF4/uU2v9TuCdAI8++qh+7LHH5jkUQRAEQRAEQWiMUuq5eR43T3cLBbwL+LzW+l9PfOs3gO8t//29wK9PbP+bZZeLtwD7VSxDEARBEARBEO4E5nGS/wzwN4DPKKUeL7f9I+CfA/9JKfV9wFeBv1J+773AdwJPA0Pgf1rqEQuCIAiCIAjCGXOqSNZaf4jZOWOAt854vAZ+aMHjEgRBEARBEIRzQybuCYIgCIIgCMIRRCQLgiAIgiAIwhFEJAuCIAiCIAjCEUQkC4IgCIIgCMIRRCQLgiAIgrDW/PrjL3CjH5/3YQh3GCKSBUEQBEFYW/bDlB/+j4/zjvc/dd6HItxhiEgWBEEQBGFtGSYZAL/5xBWyvDjnoxHuJEQkC4IgCIKwtoRJDsCNQczT7/sZGNw45yMS7hREJAuCIAiCsLaEqRHJ36Ce4Y0f/V/hs796zkck3CmISBYEQRAEYW2pnOS/0v4UAGl4cJ6HI9xBiEgWBEEQBGFtMU6y5i+4HwfguSvXzveAhDsGEcmCIAiCIKwtYZLzNeqrbEWXAej3D8/5iIQ7Bee8D0AQBEEQBOGsCNOct9kfQyuLqHAo4sF5H5JwhyBOsiAIgiAIa0uU5rzN+jjJK/4Uu2oLnYhIFuZDRLIgCIIgCGtLNtznDdZlitf8VyRWgEqH531Iwh2CiGRBEARBENaWJDbjqJ32JpndwsrCcz4i4U5BRLIgCIIgCGtLnCYAOI5LZrewcxHJwnycKpKVUj+rlLqqlHpyYtsvKaUeL/98RSn1eLn9VUqpcOJ7P3WWBy8IgiAIgnASaWJEsrIdcFu4IpKFOZmnu8XPAf8WeHe1QWv916p/K6V+DNifePwzWus3L+sABUEQBEEQmpKUIhnLQbsdfB2htUYpdb4HJqw8pzrJWusPAjdnfU+Zd9hfBd6z5OMSBEEQBEFYmDQdi2TLbxMQ04+z8z0o4Y5g0UzytwIva62fmtj2aqXUp5RSf6CU+tYF9y8IgiAIgtCYJE3NPywby+/SJmZvmJ7vQQl3BIuK5O9h2kW+ArxSa/2NwN8DflEp1Zv1g0qptyulHlNKPXbtmoyIFARBEARh+WQTTrLjd2gRsztMzveghDuCxiJZKeUA/x3wS9U2rXWstb5R/vsTwDPA62f9vNb6nVrrR7XWj166dKnpYQiCIAiCIBxLOnKSHbxWF19l7PaleE84nUWc5P8a+ILW+nK1QSl1SSlll/9+DfAI8OxihygIgiAIgtCMLBuLZL/dBaB/eHCORyTcKczTAu49wEeANyilLiulvq/81ndza8HetwFPKKU+Dfwy8INa65lFf4IgCIIgCGdNNpFJDtobAAz6IpKF0zm1BZzW+nuO2f4/ztj2K8CvLH5YgiAIgiAIi5NlVSbZpdURkSzMj0zcEwRBEARhbckn4ha23wEgGhye4xEJdwoikgXhdvHkr8LzHzvvoxAEQbiryLKyJ7LlgFuK5LB/jkc0J598N/zq28/7KO5qRCQLwu3iff8Y/lgmtQuCINxOigknGa8NQHIniOSvfBi+9DvL29+Xfhde+OTy9ncXICJZEG4XSR9y6c0pCIJwu8jyAorKSbbBNSI5C++AuEU6gGSwvP399t+HD/+b5e3vLkBEsiDcLpIBZCKSBUEQbhdRVmCTmy8sZyyS4yWKz7MiDaFIl3feiA7MPoW5EZEsCLeDPIM8FidZEAThNhImOQ6F+WIiblEkw3M8qjmpBG26JEGf9CGLlrOvuwQRyYJwO6gWuTw93+MQBEG4i4jSfKaTbKVDE8VYZaqoxTIiF3lqBLKI5FqISBaE20HlWoiTLAiCMD+/9rfhw+9o/ONhmuOMRLINnulu0SJmL1xx06JykpfhesdlBltEci1OHSYiCMISqJyAPD7f4xAEQbiT+PIHIWleZBcmObaaiFvYHhqLtorZGyZc7PpLOtAzIC3FcbKEThzVPlIRyXUQJ1kQbgfVAiVxC0EQhPlJ+pA1NxemnWQHlCJ3WrRI2B2u+Ho8EslLiFvE5TlInORaiEgWhNtBKnELQRCE2qTDhToyhEmOPVm4B2i3TYuI3cGKr8fJMkVyFbeQu5l1EJEsCLeDUdxixRdlQRCEVSHPzJq5LCfZLhOmXsfELVY5k1wUkC2xu0UVWcmkBVwdRCQLwu2giltIn2RBEIT5qMThAsLOOMkTcQtAuW3axAzjbNEjPDsmYxELOsm//eQVvvjVK+V+xUmug4hkQbgdSHcLQRCEelTr5sJO8nTcwvLbtIgJ0xVuAZdOdLRYUCT/+Pue4vefeNZ8kUWg9UL7u5sQkSwIt4NE+iQLa0q0D0V+3kchrCPVurlAsVl0tHAPUF6HlkoIk9V0kveGCZ997qXxhgVFcpjmxIP98QZxk+dGRLIgnDEv7UcMBwfmC3GShXUiT+H/+np4/BfP+0iEdaSKWyzQtsxM3CtFsrLNX16HjooJ09W8uPuJDzzND//7D483LCiS4ywnCw/GG6TDxdyISBaEM+YH/8Mn+NDnnjNf5Inc6hLWh2QA0R4cvHjeRyKsI0uKW7hWAcoCq5Q8bouOihkmqymS94YpO+7Y5dYL9kmO0oI2E8JYRPLcnCqSlVI/q5S6qpR6cmLbP1FKvaCUerz8850T3/uHSqmnlVJfVEr9N2d14IJwp3BzkJCFVTN8DcVq3uIThNpU4kWG5AhnwTIK99KcwNKjqAUAZeHeqjrJcVZwf3ucl07CxURynOV0mcg4i0iem3mc5J8D3jZj+49rrd9c/nkvgFLqa4HvBv5E+TP/t1Ll/Q1BuEuJsxw1WYQhkQthXahOtvKeFpbM01cPubG7a77IE9MSrQFRmuPbxbRI9joEJIQr6iQnWUHXGtevhP2DEx59MlprorSgoyadZLmonZdTRbLW+oPAzTn3913Af9Rax1rrLwNPA9+ywPEJwh1PnBXYmYhkYQ2pTrbS2lBYMj/8Hx/nNz7+9HhDQ/czTHL8GU5yQLTCTnI+EsmxdoiHzUVykpuLi+5k3GKB4Sx3G4tkkv+uUuqJMo6xXW57AHh+4jGXy22CcNcSpwVuPrEoSYcLYV2oYhZy4ScsmYMopX842ZGhoUhOK5E8cVPbbeOQk8SrGTuIs4KuZT5bN+iRRc0L9+KsEskT5yBxkuemqUj+SeC1wJuBK8CPldvVjMfOrFJSSr1dKfWYUuqxa9euNTwMQVh9krzAK2SBEtaQTESycDbEaUE8PBxvaLhuhmmBZx2NW7QB0MsY93wGJFlB2zKfqX21tVDhXlS65R0VEqqW2ShT9+amkUjWWr+stc611gXw04wjFZeBhyYe+iAws+xZa/1OrfWjWutHL1261OQwBGHlyfKCvNDTeTARFMK6ULl7cuEnLJkozXHyyWKzZsIuSnL8oyLZrUTyaorFOCvoKPOZCr1t1AJjqeNyYEpXRdykZzbK53VuGolkpdT9E1/+ZaDqfPEbwHcrpXyl1KuBR4CPLXaIgnDnUt3qajGxKEncQlgXpHBPOCPirJheNxs7yTmepcFyxxu9jvl7hZ3kFuYzlQUXcRZwfuPMOMmbVsTL+YbZKN0t5sY57QFKqfcAfxa4qJS6DPxj4M8qpd6MiVJ8BfgBAK31Z5VS/wn4HJABP6S1Xs1kvCDcBiqR3FERubuBnR6KoBDWB4lbCGeA1po4K2g7E8K4YbFZmOZ4qjiSSTaxA2uyoHqFiLOctp+A7eF0tgj2Q5KswHPq+5pR6SR3CLmhjZO8f3CIFaVsBO5JPyowh0jWWn/PjM3vOuHxPwr86CIHJQjrQnUV3yIm9TZFJAvrhTjJwhlQdWRoqyU4yUmO6+Uz4xZqRbs8xFlBy4/BbdHu9GgTc3l3yKsvdRvsK0dREOhoJJL/ze98hvjKN/Cjf/nrln3oa4dM3BOEMyQeXcVHRG6ZBxNBcWczuAH/7jth//J5H8n5U7V+kxZwwhKZGVNrmklOc1x1a59kACsP0Ss4ATXJCgJicNt0NjZxVc7z1/ca7StOCzpl+7cbZSY5T0Ju9OUzOw8ikgXhDImzAouCQKWE9hJEcpHDL38fvPDJ5RygUJ+rn4PnPgzPfeS8j+T8GTnJUggkLI9RRwYi8kqmLJBLine truncated
"text/plain": [
"<Figure size 864x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def create_ra_features(x, order):\n",
" y = x.copy()\n",
" # create features by shifting the window of `order` size by one step.\n",
" x = np.array([y[-(i + order):][:order] for i in range(y.shape[0])])\n",
" \n",
" # make sure that we end with even array so that the stack operation succeeds\n",
" x = np.stack(x)[::-1][order - 1: -1]\n",
" y = y[order:]\n",
" \n",
" return x, y\n",
"\n",
"def create_ma_features(eps, y, order):\n",
" eps = np.array([eps[-(i + order):][:order] for i in range(y.shape[0])])\n",
" eps = np.stack(eps)[::-1][order - 1:]\n",
" # eps_t is known when predicting x_t so the labels are in the same phase /lag\n",
" # as the features eps.\n",
" y = y[order - 1:] \n",
" return eps, y\n",
"\n",
"\n",
"# https://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model \n",
"def ra_model(y, eps, order=2):\n",
" x, y = create_ra_features(y) \n",
" \n",
" # Move eps_t to other side of equation\n",
" add = eps[:y.shape[0]]\n",
" y -= add\n",
" \n",
" m = linear_model.LinearRegression(fit_intercept=True)\n",
" m.fit(x, y)\n",
" \n",
" return m.predict(x) + add, m\n",
"\n",
"\n",
"class ARMA:\n",
" def __init__(self, p, q, sigma=None):\n",
" self.p = p\n",
" self.q = q\n",
" self.MA = None\n",
" self.RA = None\n",
" self.m = None\n",
" self.sigma = sigma\n",
" \n",
" def prepare_features(self, x):\n",
" eps = init_eps(self.q, x)\n",
" \n",
" # create features by shifting the window of `order` size by one step.\n",
" x_ = None\n",
" ma_features = ra_features = np.zeros((x.shape[0], 1))\n",
" if self.q > 1:\n",
" ma_features, x_ = create_ma_features(eps, x, self.q + 1)\n",
" if self.p > 0:\n",
" ra_features, x = create_ra_features(x, self.p)\n",
" x = x if x_ is None or len(x) < len(x_) else x_\n",
" del x_\n",
" \n",
" eps = eps[:x.shape[0]]\n",
" ma_features = ma_features[:x.shape[0]]\n",
" ra_features = ra_features[:x.shape[0]]\n",
" eps_t = ma_features[:, :1]\n",
" ma_features = ma_features[:, 1:]\n",
" return np.hstack((ma_features, ra_features)), eps_t.ravel(), x\n",
" \n",
" def fit_predict(self, x): \n",
" # move eps_t to left side of the equation\n",
" \n",
" features, eps_t, x = self.prepare_features(x)\n",
"\n",
" x = x - eps_t\n",
" model = linear_model.LinearRegression(fit_intercept=True)\n",
" model.fit(features, x)\n",
" self.m = model\n",
" \n",
" return model.predict(features) + eps_t\n",
" \n",
" \n",
"q = 0\n",
"p = 3\n",
"\n",
"m = ARMA(p, q, sigma=1)\n",
"pred = m.fit_predict(x)\n",
"print(m.m.intercept_, m.m.coef_)\n",
"# print(m.RA.intercept_)\n",
"# m.RA.coef_\n",
"\n",
"plt.figure(figsize=(12,4))\n",
"plt.plot(data)\n",
"plt.plot(data.index[data.shape[0] - pred.shape[0]:], pred)"
]
},
{
"cell_type": "code",
"execution_count": 354,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/opt/miniconda3/lib/python3.6/site-packages/statsmodels/tsa/base/tsa_model.py:171: ValueWarning: No frequency information was provided, so inferred frequency A-DEC will be used.\n",
" % freq, ValueWarning)\n",
"/opt/miniconda3/lib/python3.6/site-packages/scipy/signal/signaltools.py:1341: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" out_full[ind] += zi\n",
"/opt/miniconda3/lib/python3.6/site-packages/scipy/signal/signaltools.py:1344: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" out = out_full[ind]\n",
"/opt/miniconda3/lib/python3.6/site-packages/scipy/signal/signaltools.py:1350: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" zf = out_full[ind]\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fc7f8e42860>]"
]
},
"execution_count": 354,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD9CAYAAABOQep/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXmQbOl51vn7zpJ5cq2qu/SibrUWqyVbWLYw7WVmsCEwZgzDYCCAwTEBmhkHGg0eBoKJGMMQMUMEKIYZxoawDTjEyCF7bGTACxKMbRCWLSG0uWW1Wq3NvUjq7XbfrbbMPNt3zjd/fOfkUjerKs/JrFtZed9fxI2ue6rq9Km6md95zvM97/sqYwyCIAiCIAiCIExwzvsCBEEQBEEQBGHdEJEsCIIgCIIgCEcQkSwIgiAIgiAIRxCRLAiCIAiCIAhHEJEsCIIgCIIgCEcQkSwIgiAIgiAIRzhVJCulXquU+k2l1JeUUl9QSv3V4vglpdSHlFJPF//dKY4rpdSPK6WeUUo9qZT6trP+IQRBEARBEARhlSziJGvgfzbGfBPwXcAPK6XeCvwN4DeMMY8Cv1H8HeCPAo8Wf94J/JOVX7UgCIIgCIIgnCGnimRjzDVjzO8UHx8CXwIeAn4A+Jniy34G+JPFxz8A/KyxfBLYVko9uPIrFwRBEARBEIQzolImWSn1euD3Ap8C7jfGXAMrpIH7ii97CHhh6tteLI4JgiAIgiAIwoXAW/QLlVJd4JeAv2aMOVBKHfulc47dMftaKfVObByDTqfz+77xG79x0UsRBEEQBEEQhFp85jOfuWmMuXra1y0kkpVSPlYg/7wx5peLw68qpR40xlwr4hTXi+MvAq+d+vaHgZePntMY8x7gPQCPPfaYefzxxxe5FEEQBEEQBEGojVLq64t83SLdLRTwXuBLxpgfm/rUB4F3FB+/A/jA1PG/WHS5+C5gv4xlCIIgCIIgCMJFYBEn+T8D/gLweaXUE8Wx/xX4e8C/UEr9EPA88GeLz/0q8MeAZ4AR8N+u9IoFQRAEQRAE4Yw5VSQbYz7G/JwxwPfO+XoD/PCS1yUIgiAIgiAI54ZM3BMEQRAEQRCEI4hIFgRBEARBEIQjiEgWBEEQBEEQhCOISBYEQRAEQRCEI4hIFgRBEARho/nAEy9xaxCf92UIFwwRyYIgCIIgbCz7Ycpf/YUn+PHfePq8L0W4YIhIFgRBEARhYxklGoB/8+Q1dJaf89UIFwkRyYIgCIIgbCxhkgFwa5jwpY/+IkQH53xFwkVBRLIgCIIgCBtLmFqR/E3q67ztI38JPv8vz/mKhIuCiGRBEARBEDaW0kn+M+0nAEhH++d5OcIFQkSyIAiCIAgbS+kk/4ngdwD4+is3zvNyhAuEiGRBEARBEDaWMMl4rXqVq0Pb3WJwKJlkYTFEJAuCIAiCsLGEacYfcR4HIDEueTw85ysSLgoikgVBEARB2FiiNOOPuJ8hvfJN3FCXMenovC9JuCCISBYEQRAEYWNJwgGPqa+g3/T9pE4AiYhkYTFEJAuCIAiCsLGkcYirDF7/AbTbwtUikoXFEJEsCIIgCMLGkiQJAJ7vW5Gched8RcJF4VSRrJT6aaXUdaXUU1PH/rlS6oniz9eUUk8Ux1+vlAqnPvdTZ3nxgiAIgiAIJ5GmViQrx8P4bbwsOucrEi4K3gJf8z7gJ4GfLQ8YY/6r8mOl1I8C0525nzXGvH1VFygIgiAIglCXJLYiGccDv03DRBhjUEqd74UJa8+pTrIx5qPA7XmfU/YV9ueA96/4ugRBEARBEJYmSSci2Wm2aRExiPX5XpRwIVg2k/zdwKvGmKenjr1BKfVZpdRHlFLfveT5BUEQBEEQapPq1H7guDjNLm1i9kbp+V6UcCFYViT/ILMu8jXgEWPM7wX+OvDPlFL9ed+olHqnUupxpdTjN27IiEhBEARBEFaPTiZOshd0aBGzO0rO96KEC0FtkayU8oA/Dfzz8pgxJjbG3Co+/gzwLPDmed9vjHmPMeYxY8xjV69erXsZgiAIgiAIx5KmpZPs0Wh1aSrN7kA6XAins4yT/IeBLxtjXiwPKKWuKqXc4uM3Ao8Czy13iYIgCIIgCPVI9cRJbra6AAwOD87xioSLwiIt4N4PfAJ4i1LqRaXUDxWf+vPcWbD3PcCTSqnPAb8IvMsYM7foTxAEQRAE4azRU05y0O4BMByISBZO59QWcMaYHzzm+H8z59gvAb+0/GUJgiAIgiAsT1YW7rkerY6IZGFxZOKeIAiCIAgbi9YTJ9lt2rhFNDw8xysSLgoikgXhbvH8p+DmM+d9FYIgCPcU0yIZvw1AFA7O8YoW5Cu/Bh/+u+d9Ffc0IpIF4W7xr/4H+OjfP++rEARBuKfIdTE4xPGgYUVychFE8hc/AL/93tWd7+XPwu2vru589wAikgXhbhHtQzo676sQBEG4Z9BZDvlkmEjpJOvwAsQt0hEkw9Wd71feBb/57tWd7x5ARLIg3C2SIWTSwF4QBOFuEekcj9z+xfGg0QFAxysUn2dFGkIWQ7ai6YCj2xBfgIeDNUJEsiDcDfIMdCgiWRAE4S4SJhkumf3LVCbZrNKhPSuSYudxVdcaH4KOVnOuewQRyYJwNyhjFqtyBARBEIRTidJs1kn2WwCoNLRRjHWmvG+sIqaXaWvUpCKSqyAiWRDuBqUTIE6yIAjC4nz43fC5X6j97WE67SS747hFi5i9cM1Ni3SFTnJSFCqKk1wJEcmCcDcoFzkdn+91CIIgXCQ+93748v9X+9vDJMNTU3ELt0GuXNoqZm+05qbFWCSvoBOHiORaiEgWhLvB2Elec+dCEARhnUiGS5kL1kmeilsoRe61aBOzO1rz9XiVmeSyYE9EciVEJAvC3WCcSV5z50IQBGGdSEc2S1uTMMnwpgv3AOO1aBGzO1zz9Tgtfu5kBZnkuHSSZTezCiKSBeFuUG51iUgWBEFYjDyzzucqnWSARoeWWvNMsjErjlsUTnJa/4HjXkREsiDcDaRwTxAEoRrlurmEsJvnJKtGhzYxo1gve4Vnh44AYz9eMm7x609d48vPXyvOK05yFUQkC8LdIJG4hSAIQiVKJ3VpJ/moSG7TIiZM17gF3PSDwZIt4P7Bh57mtz73rP2LDq1LLSyEiGRBuBuM4xZrvL0nCIKwToy7AtV3km2f5FmR7DTatFVMmKynk7w3SvjC169NDiwZtwjTjGi4PzkgZs3CiEgWhDPmlf2I0bCsLJatLmGDyDT8o++EL37wvK9E2ERW4SQnR4aJYOMWHRUTptmyV3gm/KPffIb/6Wf/4+TAknGLWGdk4dQ4aulwsTAikgXhjHnXz32G//CFr9m/ZIlsdQmbQzqEG1+2fwRh1Yyd5PqiLkwzfCcHFDiF5PHbtFXCKFlPkbw3StnxJy63iZdzkqM0p8V0fENE8qKcKpKVUj+tlLqulHpq6tjfVkq9pJR6ovjzx6Y+9zeVUs8opb6ilPrPz+rCBeGicHuYoMdP8cZWbAvCJlA6fLJDIpwF48K95URy4OSTzhYAjTZt1tdJjnXOg+1JXjoJlxPJsc7oTotkcZIXZhEn+X3A9885/g+MMW8v/vwqgFLqrcCfB35P8T3/WCnlrupiBeEiEusMNV14IXkwYVMoxbG8poUV88z1Q27t7dm/ZHHtHbgozWi6ZlYk+x1buLemTnKic7rOpH4lHB7UPpcxhijN6SgRyXU4VSQbYz4K3F7wfD8A/IIxJjbGfBV4BviOJa5PEC48sc5x9bRIFtdN2BBEJAtnxF/9hSf4wKefnhyoKezCJKPpHBXJLQKitS3ci3VGz7HvqdS4xEuI5CSzjnSXqd+fiOSFWSaT/D8qpZ4s4hg7xbGHgBemvubF4pgg3LPEaY6XTT3FS4cLYVMob7YStxBWzEGUMjic6shQVySP4xZTm9qNNi45SbKer9tY53QKJ/k2PXRUP24RayuSOzNxi/X8udeRuiL5nwDfALwduAb8aHFczfnauXskSql3KqUeV0o9fuPGjZqXIQjrT5LlNPJpkSyum7AhjJ1kefATVkuc5sSjqY4MNXPJYZrTuMNJ7gBgluwacVYkOqdTOMl7apt8ieuMitx1V4XENO1Bmbq3MLVEsjHmVWNMZozJgX/KJFLxIvDaqS99GHj5mHO8xxjzmDHmsatXr9a5DEFYe3SWk+WGjppa4EUkC5tC6e5JhEhYMVGa4erlc7RRktGYU7gHTIY8rRmxzuko+54aNS7hpPVFclwMTOmqiNv07UFxkhemlkhWSj049dc/BZSdLz4I/HmlVFMp9QbgUeDTy12iIFxcyq2uFlOLkhaRLGwIY5Esr2lhtcQ6p62m1836cYuGk4PrTw766y2SE53TUvY9pYPLeLr+dcbaOslbTsz1vGcPLjGc5V7DO+0LlFLvB/4gcEUp9SLwvwN/UCn1dmyU4mvAfw9gjPmCUupfAF8ENPDDxpj1LB8VhLvAJA8Wkfkd3HQogkLYHMrXsjz4CSvEGEOsc1re8sVmYZrRUEcyyYVIdpYQn2dJrDPabgzKxetcItgPSXROw6vua0alk0zILWM3+vcPD3GilF7gn/StAguIZGPMD845/N4Tvv7dwLuXuShB2BTKp/iWikn9fiGSJb8pbAgStxDOgLIjQ3tmB67eayxMMhqNY+IWa5rNjXVOqxlDo0Or06dFzEt7IW+40qlxrgwwBCbklrFxi3/4658neflbefefetuKr3zzkIl7gnCGlHmwDjGxv2UPipN8sUlG8Gs/AtH+6V+76UjhnnAGlDtw7elajpqCNkozfJXNLdxzsxFmDSegJjonIAa/Rae3RUNlvHBjr9a54jSnTYzCcAt7D8qSkFsDuQ8tgohkQThDYp2jsNm6kVsUTSzjuhkDH3433Hp2NRcoVOflz8Knfgqe/c3zvpLzR1rACWdA2ZGhTUxeNs2q6ySnGb6a7yS3iMeLine truncated
"text/plain": [
"<Figure size 864x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"p = 3\n",
"q = 0\n",
"m = sm.tsa.ARMA(data.copy(), (p, q)).fit()\n",
"plt.figure(figsize=(12,4))\n",
"plt.plot(data)\n",
"plt.plot(data.index, m.model.predict(m.params, '1700', '2008'))"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/opt/miniconda3/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
" return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAHjCAYAAADhZpMBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl4nVW5///3SucAZShlJi0CiiAUMaKIMk8yFZC5bXLQrxV6VPQ406Ny1IoIeET8gfQc9ezd7k5QZiiTgIjKUAYZlEGkLQhIaaHQlo5Zvz9WQ9OQYWcnO0+y835dF9fOfvJkP3ctyKeLe90rxBiRJEmS1HFVWRcgSZIk9VaGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRP2zLqAjttxyyzhy5Misy5AkSVKFe/jhh1+PMQ5v775eFaZHjhzJ3Llzsy5DkiRJFS6EML+Y+2zzkCRJkkpkmJYkSZJKZJiWJEmSSmSYliRJkkqU6QbEEMI84G1gLbAmxlibZT2SJElSR/SEaR4Hxxhfz7oISZIkqaNs85AkSZJKlHWYjsDtIYSHQwjjW7ohhDA+hDA3hDB34cKF3VyeJEmS1Lqsw/T+McZ9gE8D/x5COKD5DTHGyTHG2hhj7fDh7R5CI0mSJHWbTMN0jPHlda+vAdcC+2ZZjyRJktQRmYXpEMJGIYRNGr8GjgCezKoeSZIkqaOynOaxNXBtCKGxjmkxxlszrEeSJEnqkMzCdIzxH8CorJ4vSZIkdVbWGxAlSZKkXsswLUmSJJXIMC1JkqSKVCjAyJFQVZVeC4Wuf0ZPOE5ckiRJ6lKFAowfD8uXp/fz56f3AGPGdN1zXJmWJElSxZk4cX2QbrR8ebrelQzTkiRJqjgLFnTseqkM05IkSao4NTUdu14qw7QkSZIqzqRJUF294bXq6nS9KxmmJUmSVHHGjIHJk2HECAghvU6e3LWbD8FpHpIkSapQY8Z0fXhuzpVpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSGaUmSJKlEhmlJkiSpRIZpSZIkqUSZh+kQQr8QwqMhhJuyrkWSJEnqiMzDNHAu8Lesi5AkSZI6KtMwHULYATgG+N8s65AkSZJKkfXK9M+BbwINrd0QQhgfQpgbQpi7cOHC7qtMkiRJakdmYTqEcCzwWozx4bbuizFOjjHWxhhrhw8f3k3VSZIkSe3LcmV6f+D4EMI8YAZwSAhhaob1SJIkSR2SWZiOMX4nxrhDjHEkcDpwV4xxbFb1SJIkSR2Vdc+0JEmS1Gv1z7oAgBjjPcA9GZchSZIkdYgr05IkSVKJDNOSJElSiQzTkiRJUokM05IkSVKJDNOSJElSiQzTkiRJUokM05IkSVKJDNOSJElSiQzTkiRJUokM05IkSVKJDNOSJElSiQzTkiRJUokM05IkSVKJDNOSJElSiQzTkiRJUokM05IkSVKJDNOSJElSiQzTkiRJUokM05IkSVKJDNOSJEnqFoUCjBwJVVXptVDIuqLO6591AZIkSap8hQKMHw/Ll6f38+en9wBjxmRXV2e5Mi1JkqSymzhxfZButHx5ut6bGaYlSZJUdgsWdOx6b2GYliRJUtnV1HTsem9hmJYkSVLZTZoE1dUbXquuTtd7M8O0JEmSym7MGJg8GUaMgBDS6+TJvXvzITjNQ5IkSd1kzJjeH56bc2VakiRJKpFhWpIkSSqRYVqSJEltqsSTC7uKPdOSJElqVaWeXNhVXJmWJElSqyr15MKuYpiWJElSqyr15MKuYpiWJElSqyr15MKuYpiWJElSqyr15MKuYpiWJElSqyr15MKu4jQPSZIktakSTy7sKq5MS5IkSSUyTEuSJEklMkxLkiRJJTJMS5IkSSUyTEuSJEklMkxLkiRJJcosTIcQBocQHgwh/CWE8FQI4b+yqkWSJEkqRZZzplcCh8QYl4YQBgD3hRDmxBjvz7AmSZIkqWiZhekYYwSWrns7YN1fMat6JEmSpI7KtGc6hNAvhPAY8BpwR4zxgSzrkSRJkjoi0zAdY1wbY9wb2AHYN4Twoeb3hBDGhxDmhhDmLly4sPuLlCRJklrRI6Z5xBjfBO4Bjmrhe5NjjLUxxtrhw4d3e22SJEmVpFCAkSOhqiq9FgpZV9S7ZTnNY3gIYbN1Xw8BDgOezqoeSZKkSlcowPjxMH8+xJhex483UHdGlivT2wJ3hxAeBx4i9UzflGE9kiRJFW3iRFi+fMNry5en6ypNltM8Hgc+nNXzJUmS+poFCzp2Xe3rET3TkiRJKr+amo5dV/sM05IkSX3EpElQXb3hterqdF2lMUxLkiT1EWPGwOTJMGIEhJBeJ09O11Uaw7QkSVIFKHbk3ZgxMG8eNDSkV4N052S2AVGSJEldo3HkXeOkjsaRd2BYLjdXpiVJkno5R95lxzAtSZLUyznyLjuGaUmSpF7OkXfZMUxLkiT1co68y45hWpIkqZdz5F12nOYhSZJUAcaMMTxnwZVpSZIkqUSGaUmSJKlEhmlJkqQertjTDdX97JmWJEnqwTzdsGdzZVqSJKkH83TDns0wLUmS1IN5umHPZpiWJEnqwTzdsGczTEuSJPVgnm7YsxmmJUmSejBPN+zZnOYhSZLUw3m6Yc/lyrQkSZJUIsO0JEmSVCLDtCRJklQiw7QkSVLGPC6893IDoiRJUoY8Lrx3c2VakiQpQx4X3rsZpiVJkrpJS+0cHhfeuxmmJUmSulhLobmxnWP+fIhxfTvHFlu0/BkeF9472DMtSZLUhVrrgR4ypOV2jiFD0vHgTb/nceG9hyvTkiRJXai1HuhFi1q+f/FijwvvzVyZliRJ6kId7XWuqfG48N7MlWlJkqQu1Fqv87BhqX2jKds5ej/DtCRJUheaNKnl0HzppbZzVCLDtCRJUolamtoxZkzroXnMGJg3Dxoa0qtBuvezZ1qSJKkE7Z1caFDuG1yZliRJKoEnF/YSy5bBihVl+3jDtCRJUgk8ubAHa2iAe+6Bs86CbbaB6dPL9ijbPCRJkkpQU5NaO1q6row89xzk8zBlSvrN2WQTOO002Hvvsj3SMC1JklSCSZM27JkGR91l4o03YNYsyOXgz39Ou0EPPxwuuABGj37vaJUuZpiWJEkqQeMGw4kTU2tHTU0K0m487AarV8Ntt6VV6BtugJUrYY894Kc/Tb8B223XbaUYpiVJkkrk1I5u9thjaQV62jR47TXYckv4whegvh4+/OE0i7CbuQFRkiSpCC3NlFY3ePVVuOQSGDUqBebLL4dPfQquvx5efjmdhrPPPpkEaXBlWpIkqV3tzZRWF3vnndS+kculdo6GBvjYx1KQPu002GKLrCt8V4gxZvPgEHYE8sA2QAMwOcZ4aVs/U1tbG+fOndsd5UmSJL1r5MiWJ3eMGJFOMlQXiBH++MfUBz1rFixZAjvuCOPGpb92261bywkhPBxjrG3vvixXptcAX4sxPhJC2AR4OIRwR4zxrxnWJEmS9B7OlC6jF15Io+zyeXj+edhoI/jMZ1If9EEHpb6aHiyz6mKMr8QYH1n39dvA34Dts6pHkiSpUfP+6Na6CpwpXaK33oJf/xoOPBDe9z44//y0zJ/LpR7pXA4OOaTHB2noIT3TIYSRwIeBB1r43nhgPECNf8dKkqQyKBTWj7jbYgt4+21YtSp9b/58GDAABg5cfw2cKd1ha9fCnXemoHzttemI7/e/P/2POHZsr/2TSeZhOoSwMTAb+EqM8a3m348xTgYmQ+qZ7ubyJElShWu+uXDRovfes3o1DBsGG2/sTOkOe+qpFKCnToVXXoHNN4fPfhbq6mDffTObwtFVMg3TIYQBpCBdiDFek2UtkiSpb5o4ccNTDFuzeDG8/nr566kICxfC9OkpRD/yCPTvD0cfnQL0scfCoEFZV9hlMgvTIYQA/Br4W4zxZ1nVIUmS+rZiNxH20i6E7rNyJdx0U9pIeMstsGZNmv986aVw+umw1VZZV1gWWa5M7w+MA54IITy27tp5McZbMqxJkiT1MTU1LY+9a8r+6FbECA8+mFagZ8yAN96AbbeFr341rUJ/6ENZV1h2mYXpGON9QO9ukpEkSb1S8w2HzTcXDhgAQ4em1g77o1vw4ovrx9k98wwMHgwnnpjG2R16aGrr6CP6zq9UkiSJljccDhiQNhgantuwdClcc01ahb777rQqfcAB8I1vwCmnpD999EGGaUmS1Ke0tOFw9eo0qcMNhs00NKTgnM/D7NmwbNn6udBjx6av+7h2w3QIYWfgpRjjyhDCQcBeQD7G+Ga5i5MkSeoKTds6YiuDdj3NsIlnnkkBesqU1NIxdCiceWbqg95//14/zq4rFbMyPRuoDSHsQpq+cQMwDTi6nIVJkiSVqq1DWFrT56d1LF6cNhHm8/DAA+n0wSOPhIsuguOPhyFDsq6wRyomTDfEGNeEEE4Efh5jvCyE8Gi5C5MkSSpFMYewNNdnp3WsXg1z5qQ+6JtuSn/i2HNPuPjitBK97bZZV9jjFROmV4cQzgDqgePWXRtQvpIkSZI6rnE1ur0xd02F0Ac3HMYIjz6aAvS0aalRfPhwmDAhTeMYNco2jg4oJkyfBZwNTIoxvhBC2AmYWt6yJEmS2lZKK0dTI0bAvHllK6/nefnl9D9aLpeO+B44EEaPTn3QRx6ZRpqow9oN0zHLine truncated
"text/plain": [
"<Figure size 864x576 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(12,8))\n",
"ax = fig.add_subplot(111)\n",
"fig = qqplot(m.resid, line='q', ax=ax, fit=True)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m.resid"
]
}
],
"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.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}