Python-计算自相关系数.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cmaps\n",
    "import numpy as np\n",
    "import scipy.signal as signal\n",
    "import matplotlib.pyplot as plt\n",
    "import netCDF4 as nc\n",
    "import os\n",
    "import xarray as xr\n",
    "import glob\n",
    "import sys\n",
    "# sys.path.append('/DatadiskExt/xpji/CCKW_code/')\n",
    "# import ma as mp\n",
    "import pickle\n",
    "import datetime\n",
    "import pandas as pd\n",
    "import calendar\n",
    "import cartopy.crs as ccrs\n",
    "import matplotlib.ticker as mticker\n",
    "import matplotlib.ticker as ticker\n",
    "import cartopy.feature as cfeature\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
    "from statsmodels.tsa.stattools import acf\n",
    "\n",
    "import time  # 导入 time 模块"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time taken by auto2: 0.001860 seconds\n",
      "Time taken by acf_by_hand: 0.000839 seconds\n",
      "Time taken by acf: 0.000565 seconds\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1280x960 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/UAAAK2CAYAAAAPALyPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABbOUlEQVR4nO3deVhV5f7//9dmBpGtokwKijNGNjikmKmVmFPapKaJmlZmHo9iw9FOTpnYoMesHCoNK1Mrq6PfTPNYljklJjZI5hia4IAKjiCwfn/0c3/aAQYJbG98Pq5rX6d9r3vd670WXHVe3Pday2ZZliUAAAAAAGAcN1cXAAAAAAAA/h5CPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AOCKNXPmTNlsNkVHR1/2WCtWrNCECRMuvygD7N+/XzabTYmJiSXe99ChQ5owYYKSk5MLbJswYYJsNtvlF/g3ZWdn69VXX9XNN9+sqlWrysvLSzVr1lSvXr301Vdfuayu8nI5P1cAQMVFqAcAXLHmz58vSfrpp5+0efPmyxprxYoVmjhxYmmUVaEdOnRIEydOLDTUDxkyRBs3biz/oiQdO3ZMbdq0UXx8vKKjo5WYmKg1a9Zo2rRpcnd312233abt27e7pLbyEhoaqo0bN6pr166uLgUAcAXxcHUBAAAUJikpSdu3b1fXrl316aefat68ebrppptcXVaZOnv2rPz8/Aq0W5al8+fPy9fX1wVV/Z9atWqpVq1aLjl2XFyctm/frlWrVunWW2912tanTx/Fx8eratWqLqmtrOXl5Sk3N1fe3t5q1aqVq8sBAFxhmKkHAFyR5s2bJ0maOnWqYmJitHjxYp09e9apz9q1a2Wz2bR27Vqn9j8vUx44cKBee+01SZLNZnN89u/fL0k6f/68xowZo8jISMeS7scee0wnT54sUNd7772n1q1by9/fX/7+/rr++usdtV40f/58XXfddfLx8VG1atV01113KSUlxanPwIED5e/vrx9++EGxsbGqXLmybrvtNkeNw4cP15w5cxQVFSVvb28tWLBAkrRr1y717dtXQUFB8vb2VlRUlOPcLmX37t0aNGiQGjRoID8/P9WsWVPdu3fXDz/84HQ9W7RoIUkaNGiQ4zpdvG2hsOX3+fn5euGFF9S4cWN5e3srKChIcXFxOnjwoFO/9u3bKzo6Wlu2bFHbtm3l5+enunXraurUqcrPz79k7Vu3btVnn32mwYMHFwj0F7Vo0UIRERGO7z/++KN69OihqlWrysfHR9dff73jGv7xfG02m9577z099dRTCg0Nlb+/v7p3767Dhw/r1KlTevjhh1W9enVVr15dgwYN0unTp53GuPizmjt3rho2bChvb281adJEixcvdup39OhRDRs2TE2aNJG/v7+CgoJ06623at26dU79Lv7uvvDCC5o8ebIiIyPl7e2tL7/8stDl90ePHtXDDz+s8PBweXt7q0aNGmrTpo3+97//OY1bkt/J3bt3q0uXLvL391d4eLhGjx6t7OzsS/6MAACuw0w9AOCKc+7cOS1atEgtWrRQdHS0HnzwQQ0ZMkQffPCBBgwYUOLxnnnmGZ05c0Yffvih0/Lx0NBQWZalnj17as2aNRozZozatm2r77//XuPHj9fGjRu1ceNGeXt7S5LGjRunZ599VnfffbdGjx4tu92uH3/8Ub/++qtjzISEBI0dO1b333+/EhISlJGRoQkTJqh169basmWLGjRo4Oibk5OjO++8U4888oj+9a9/KTc317Htk08+0bp16zRu3DiFhIQoKChIO3bsUExMjCIiIjRt2jSFhIRo1apVGjFihI4dO6bx48cXeQ0OHTqkwMBATZ06VTVq1NDx48e1YMEC3XTTTdq2bZsaNWqkG2+8UW+99ZYGDRqkf//7345l3peanX/00Uf1+uuva/jw4erWrZv279+vZ555RmvXrtV3332n6tWrO/qmp6erX79+Gj16tMaPH6+PP/5YY8aMUVhYmOLi4oo8xueffy5J6tmzZ5F9/mjnzp2KiYlRUFCQZs6cqcDAQL377rsaOHCgDh8+rCeffNKp/9ixY9WhQwclJiZq//79evzxx3X//ffLw8ND1113nRYtWqRt27Zp7Nixqly5smbOnOm0/7Jly/Tll19q0qRJqlSpkmbNmuXY/95775UkHT9+XJI0fvx4hYSE6PTp0/r444/Vvn17rVmzRu3bt3cac+bMmWrYsKFeeuklBQQEOP3e/FH//v313Xff6bnnnlPDhg118uRJfffdd8rIyHD0Kcnv5IULF3TnnXdq8ODBGj16tL7++ms9++yzstvtGjduXLGuPwCgnFkAAFxh3n77bUuSNWfOHMuyLOvUqVOWv7+/1bZtW6d+X375pSXJ+vLLL53a9+3bZ0my3nrrLUfbY489ZhX2n72VK1dakqwXXnjBqX3JkiWWJOv111+3LMuy9u7da7m7u1v9+vUrsu4TJ05Yvr6+VpcuXZzaU1NTLW9vb6tv376OtgEDBliSrPnz5xcYR5Jlt9ut48ePO7V36tTJqlWrlpWZmenUPnz4cMvHx8fRv7Dz/7Pc3FwrJyfHatCggTVq1ChH+5YtW4rcd/z48U7XMCUlxZJkDRs2zKnf5s2bLUnW2LFjHW3t2rWzJFmbN2926tukSROrU6dORdZpWZY1dOhQS5L1888/X7LfRX369LG8vb2t1NRUp/bOnTtbfn5+1smTJy3L+r/fn+7duzv1GzlypCXJGjFihFN7z549rWrVqjm1SbJ8fX2t9PR0R1tubq7VuHFjq379+kXWmJuba124cMG67bbbrLvuusvRfvFnV69ePSsnJ8dpn8J+rv7+/tbIkSOLPM7f+Z18//33nfp26dLFatSoUZHHAAC4FsvvAQBXnHnz5snX11d9+vSRJPn7++u+++7TunXrtGvXrlI91hdffCHp96XHf3TfffepUqVKWrNmjSRp9erVysvL02OPPVbkWBs3btS5c+cKjBUeHq5bb73VMdYf3XPPPYWOdeuttzrdI37+/HmtWbNGd911l/z8/JSbm+v4dOnSRefPn9emTZuKrC03N1dTpkxRkyZN5OXlJQ8PD3l5eWnXrl0FlmEX15dffimp4LVr2bKloqKiCpxvSEiIWrZs6dTWtGlTp5UOpeGLL77QbbfdpvDwcKf2gQMH6uzZswUe9tetWzen71FRUZJU4IF0UVFROn78eIEl+LfddpuCg4Md393d3dW7d2/t3r3b6TaEOXPm6MYbb5SPj488PDzk6empNWvWFHr977zzTnl6ev7lubZs2VKJiYmaPHmyNm3apAsXLjhtL+nvpM1mU/fu3Z3ayuJnBAAoPYR6AMAVZffu3fr666/VtWtXWZalkydP6uTJk45lzBefiF9aMjIy5OHhoRo1aji122w2hYSEOJYxHz16VNKll6Jf7BsaGlpgW1hYmNOSaEny8/NTQEBAoWP9eYyMjAzl5ubqlVdekaenp9OnS5cukn5/QnxR4uPj9cwzz6hnz55avny5Nm/erC1btui6667TuXPnitzvUkp6voGBgQX6eXt7/+XxL94rv2/fvmLXVVRNf6z7omrVqjl99/LyumT7+fPnndpDQkIKHOti28VjTZ8+XY8++qhuuukmLV26VJs2bdKWLVt0xx13FHr+hdVfmCVLlmjAgAF688031bp1a1WrVk1xcXFKT093On5Jfid9fHyc2ry9vQucMwDgysE99QCAK8r8+fNlWZY+/PBDffjhhwW2L1iwQJMnT5a7u7sjfPz5IV6XCrd/FhgYqNzcXB09etQp2FuWpfT0dMeD4y5uO3jwYIEZ4D+OJUlpaWkFth06dMjp/nJJl3zn+5+3Va1aVe7u7urfv3+RqwUiIyOLHO/dd99VXFycpkyZ4tR+7NgxValSpcj9LuWP5/vnP3YUdr5/V6dOnTR27Fh98sknuuOOO4pVV1E/A0mlVtdFFwN0YW0Xr9G7776r9u3ba/bs2U79Tp06VeiYl/rd+KPq1atrxowZmjFjhlJTU7Vs2TL961//0pEjR7Ry5coS/04CAMzDTD0A4IqRl5enBQsWqF69evryyy8LfEaPHq20tDR99tlnkqQ6depIkr7//nuncZYtW1Zg7IsPu/vzrOjFJ86/++67Tu1Lly7VmTNnHNtjY2Pl7u5eIJT9UevWreXr61tgrIMHDzqWhP9dfn5+6tChg7Zt26amTZuqefPmBT6FzYRfZLPZHNfgok8//VS//fabU1tR16kwF59E/+fz3bJli1JSUi7rfP/oxhtvVOfOnTVv3jzH7RJ/lpSUpNTUVEm//0y/+OILR4i/6O2335afn1+pvxZuzZo1Onz4sON7Xl6elixZonr16jn+2FHY9f/+++8L3ApwOSIiIjR8+HB17NhR3333naSy/Z0EAFwZmKkHAFwxPvvsMx06dEjPP/98gaeBS1J0dLReffVVzZs3T926dVNISIhuv/12JSQkqGrVqqpdu7bWrFmjjz76qMC+1157rSTp+eefV+fOneXu7q6mTZuqY8eO6tSpk5566illZWWpTZs2jqff33DDDerfv7+k3/+AMHbsWD377LM6d+6c7r//ftntdu3YsUPHjh3TxIkTVaVKFT3zzDMaO3as4uLidP/99ysjI0MTJ06Uj4/PJZ9OXxwvv/yybr75ZrVt21aPPvqo6tSpo1OnTmn37t1avnx5kYFX+v2+8cTERDVu3FhNmzbV1q1b9eKLLxaYYa9Xr558fX21cOFCRUVFyd/fX2FhYY6l63/UqFEjPfzww3rllVfk5uamzp07O55+Hx4erlGjRl3W+f7R22+/rTvuuEOdO3fWgw8+qM6dO6tq1apKS0vT8uXLtWjRIm3dulUREREaP368/t//+3/q0KGDxo0bp2rVqmnhwoX69NNP9cILL8hut5daXdLvs+W33nqrnnnmGcfT73/++Wen19p169ZNzz77rMaPH6927dpp586dmjRpkiIjI53eelASmZmZ6tChg/r27avGjRurcuXK2rJli1auXKm7775bksr8dxIAcAVw8YP6AABw6Nmzp+Xl5WUdOXKkyD59+vSxPDw8HE8bT0tLs+69916rWrVqlt1utx544AErKSmpwFPCs7OzrSFDhlg1atSwbDabJcnat2+fZVmWde7cOeupp56yateubXl6elqhoaHWo48+ap04caLA8d9++22rRYsWlo+Pj+Xv72/dcMMNBZ4U/+abb1pNmza1vLy8LLvdbvXo0cP66aefnPoMGDDAqlSpUqHnKMl67LHHCt22b98+68EHH7Rq1qxpeXp6WjVq1LBiYmKsyZMnO/X58/mfOHHCGjx4sBUUFGT5+flZN998s7Vu3TqrXbt2Vrt27ZyOsWjRIqtx48aWp6enJckaP368ZVkFn35vWZaVl5dnPf/881bDhg0tT09Pq3r16tYDDzxgHThwwKlfu3btrGuuuabA+QwYMMCqXbt2oef6Z+fOnbNmzpxptW7d2goICLA8PDyssLAw6+6777Y+/fRTp74//PCD1b17d8tut1teXl7WddddV+DndPHp9x988IFT+1tvvWVJsrZs2eLUfvH8jx496mi7+LOaNWuWVa9ePcvT09Nq3LixtXDhQqd9s7Ozrccff9yqWbOm5ePjY914443WJ598UuD8L/7sXnzxxQLn/+ef6/nz562hQ4daTZs2tQICAixfX1+rUaNG1vjx460zZ8447Xs5v5OF/dwBAFcOm2VZlkv+mgAAAGA4m82mxx57TK+++qqrSwEAXKW4px4AAAAAAEMR6gEAAAAAMBQPygMAAPibuIsRAOBqRs3Uf/311+revbvCwsJks9n0ySefXLL/Rx99pI4dO6pGjRoKCAhQ69attWrVqvIpFgAAAACAMmZUqD9z5oyuu+66Yj+M5uuvv1bHjh21YsUKbd26VR06dFD37t21bdu2Mq4UAAAAAICyZ+zT7202mz7++GP17NmzRPtdc8016t27t8aNG1c2hQEAAAAAUE6uqnvq8/PzderUKVWrVq3IPtnZ2crOznba5/jx4woMDJTNZiuPMgEAAAAAVzHLsnTq1CmFhYXJze3SC+yvqlA/bdo0nTlzRr169SqyT0JCgiZOnFiOVQEAAAAAUNCBAwdUq1atS/a5apbfL1q0SEOGDNF///tf3X777UX2+/NMfWZmpiIiInTgwAEFBARcbtkAAAAAAFxSVlaWwsPDdfLkSdnt9kv2vSpm6pcsWaLBgwfrgw8+uGSglyRvb295e3sXaA8ICCDUAwAAAADKTXFuATfq6fd/x6JFizRw4EC999576tq1q6vLAQAAAACg1Bg1U3/69Gnt3r3b8X3fvn1KTk5WtWrVFBERoTFjxui3337T22+/Len3QB8XF6eXX35ZrVq1Unp6uiTJ19f3L5cwAAAAAABwpTNqpj4pKUk33HCDbrjhBklSfHy8brjhBsfr6dLS0pSamuroP3fuXOXm5uqxxx5TaGio4/PPf/7TJfUDAAAAAFCajH1QXnnJysqS3W5XZmYm99QDAAAAKJH8/Hzl5OS4ugxcYTw9PeXu7l7k9pLkUKOW3wMAAACAKXJycrRv3z7l5+e7uhRcgapUqaKQkJBiPQzvUgj1AAAAAFDKLMtSWlqa3N3dFR4eLjc3o+58RhmyLEtnz57VkSNHJEmhoaGXNR6hHgAAAABKWW5urs6ePauwsDD5+fm5uhxcYXx9fSVJR44cUVBQ0CWX4v8V/lwEAAAAAKUsLy9PkuTl5eXiSnCluvjHngsXLlzWOIR6AAAAACgjl3u/NCqu0vrdINQDAAAAAGAoQj0AAAAAAIbiQXkAAAAAUE5so5eX6/Gsad3L9XiSlJiYqJEjR+rkyZMl2m/t2rX6z3/+o2+//VZZWVlq0KCBnnjiCfXr169sCq0gmKkHAAAAALjchg0b1LRpUy1dulTff/+9HnzwQcXFxWn58vL9Q4hpCPUAAAAAAIeVK1fq5ptvVpUqVRQYGKhu3bppz549kn6fTbfZbE6z8MnJybLZbNq/f7/Wrl2rQYMGKTMzUzabTTabTRMmTJAknThxQnFxcapatar8/PzUuXNn7dq1yzHO2LFj9eyzzyomJkb16tXTiBEjdMcdd+jjjz8uz9M3DqEeAAAAAOBw5swZxcfHa8uWLVqzZo3c3Nx01113KT8//y/3jYmJ0YwZMxQQEKC0tDSlpaXp8ccflyQNHDhQSUlJWrZsmTZu3CjLstSlS5dLvtItMzNT1apVK7Vzq4i4px4AAAAA4HDPPfc4fZ83b56CgoK0Y8eOv9zXy8tLdrtdNptNISEhjvZdu3Zp2bJlWr9+vWJiYiRJCxcuVHh4uD755BPdd999Bcb68MMPtWXLFs2dO/cyz6hiY6YeAAAAAOCwZ88e9e3bV3Xr1lVAQIAiIyMlSampqX97zJSUFHl4eOimm25ytAUGBqpRo0ZKSUkp0H/t2rUaOHCg3njjDV1zzTV/+7hXA2bqAQAAAAAO3bt3V3h4uN544w2FhYUpPz9f0dHRysnJkb+/vyTJsixH/0stn7/oj/3/3G6z2ZzavvrqK3Xv3l3Tp09XXFzcZZzJ1YGZegAAAACAJCkjI0MpKSn697//rdtuu01RUVE6ceKEY3uNGjUkSWlpaY625ORkpzG8vLyUl5fn1NakSRPl5uZq8+bNTsf65ZdfFBUV5Whbu3atunbtqqlTp+rhhx8uzVOrsAj1AAAAAABJUtWqVRUYGKjXX39du3fv1hdffKH4+HjH9vr16ys8PFwTJkzQL7/8ok8//VTTpk1zGqNOnTo6ffq01qxZo2PHjuns2bNq0KCBevTooYceekjffPONtm/frgceeEA1a9ZUjx49JP1foB8xYoTuuecepaenKz09XcePHy/Xa2Aalt8DAAAAQDmxpnV3dQmX5ObmpsWLF2vEiBGKjo5Wo0aNNHPmTLVv316S5OnpqUWLFunRRx/VddddpxYtWmjy5MlOD7qLiYnR0KFD1bt3b2VkZGj8+PGaMGGC3nrrLf3zn/9Ut27dlJOTo1tuuUUrVqyQp6enJCkxMVFnz55VQkKCEhISHOO1a9dOa9euLc/LYBSbVdTNDZAkZWVlyW63KzMzUwEBAa4uBwAAAIABzp8/r3379ikyMlI+Pj6uLgdXoEv9jpQkh7L8HgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAAKXm9ddfV3h4uNzc3DRjxoxL9p0wYYKuv/76MqslMTFRVapUKbPxL2Xt2rWy2Ww6efJkmR7Ho0xHBwAAAAA42N56vFyPZw16qVyPl5WVpeHDh2v69Om65557ZLfby/X4VyNCPQAAAACgVKSmpurChQvq2rWrQkNDXV3OVYHl9wAAAAAAh5UrV+rmm29WlSpVFBgYqG7dumnPnj2O7QcPHlSfPn1UrVo1VapUSc2bN9fmzZuVmJioa6+9VpJUt25d2Ww27d+/v1jHnDt3rsLDw+Xn56f77rvPsWT966+/lqenp9LT0536jx49Wrfcckuxz2nVqlWKioqSv7+/7rjjDqWlpTm2bdmyRR07dlT16tVlt9vVrl07fffdd07722w2vfnmm7rrrrvk5+enBg0aaNmyZU59VqxYoYYNG8rX11cdOnQo9rlfLkI9AAAAAMDhzJkzio+P15YtW7RmzRq5ubnprrvuUn5+vk6fPq127drp0KFDWrZsmbZv364nn3xS+fn56t27t/73v/9Jkr799lulpaUpPDz8L4+3e/duvf/++1q+fLlWrlyp5ORkPfbYY5KkW265RXXr1tU777zj6J+bm6t3331XgwYNKtb5nD17Vi+99JLeeecdff3110pNTdXjj//fbRCnTp3SgAEDtG7dOm3atEkNGjRQly5ddOrUKadxJk6cqF69eun7779Xly5d1K9fPx0/flySdODAAd19993q0qWLkpOTNWTIEP3rX/8qVn2Xi+X3AAAAAACHe+65x+n7vHnzFBQUpB07dmjDhg06evSotmzZomrVqkmS6tev7+gbGBgoSapRo4ZCQkKKdbzz589rwYIFqlWrliTplVdeUdeuXTVt2jSFhIRo8ODBeuutt/TEE09Ikj799FOdPXtWvXr1Ktb4Fy5c0Jw5c1SvXj1J0vDhwzVp0iTH9ltvvdWp/9y5c1W1alV99dVX6tatm6N94MCBuv/++yVJU6ZM0SuvvKJvv/1Wd9xxh2bPnq26devqP//5j2w2mxo1aqQffvhBzz//fLFqvBzM1AMAAAAAHPbs2aO+ffuqbt26CggIUGRkpKTf75dPTk7WDTfc4Aj0pSEiIsIR6CWpdevWys/P186dOyX9HqZ3796tTZs2SZLmz5+vXr16qVKlSsUa38/PzxHoJSk0NFRHjhxxfD9y5IiGDh2qhg0bym63y2636/Tp00pNTXUap2nTpo5/rlSpkipXruwYJyUlRa1atZLNZnM6j/LATD0AAAAAwKF79+4KDw/XG2+8obCwMOXn5ys6Olo5OTny9fUt8+NfDMYX/zcoKEjdu3fXW2+9pbp162rFihVau3Ztscfz9PQsML5lWY7vAwcO1NGjRzVjxgzVrl1b3t7eat26tXJycv5ynPz8fElyGq+8EeoBAAAAAJKkjIwMpaSkaO7cuWrbtq0k6ZtvvnFsb9q0qd58800dP3681GbrU1NTdejQIYWFhUmSNm7cKDc3NzVs2NDRZ8iQIerTp49q1aqlevXqqU2bNqVybElat26dZs2apS5dukj6/f74Y8eOlWiMJk2a6JNPPnFqu7iyoKyx/B4AAAAAIEmqWrWqAgMD9frrr2v37t364osvFB8f79h+//33KyQkRD179tT69eu1d+9eLV26VBs3bvzbx/Tx8dGAAQO0fft2rVu3TiNGjFCvXr2c7snv1KmT7Ha7Jk+eXOwH5BVX/fr19c477yglJUWbN29Wv379SrwiYejQodqzZ4/i4+O1c+dOvffee0pMTCzVOovCTD0AAAAAlBNr0EuuLuGS3NzctHjxYo0YMULR0dFq1KiRZs6cqfbt20uSvLy89Pnnn2v06NHq0qWLcnNz1aRJE7322mt/+5j169d3PDn++PHj6tKli2bNmlWgroEDB2rKlCmKi4u7nFMsYP78+Xr44Yd1ww03KCIiQlOmTHF6On5xREREaOnSpRo1apRmzZqlli1basqUKXrwwQdLtdbC2CxXLv43QFZWlux2uzIzMxUQEODqcgAAAAAY4Pz589q3b58iIyPl4+Pj6nIqhIceekiHDx8u8H54U13qd6QkOZSZegAAAADAFSszM1NbtmzRwoUL9d///tfV5VxxuKceAAAAAFAmrrnmGvn7+xf6WbhwYbHG6NGjh+6880498sgj6tixo9O2zp07Fzn+lClTyuKUrjjM1AMAAAAAysSKFSt04cKFQrcFBwcXa4xLvb7uzTff1Llz5wrdVlpP57/SEeoBAAAAAGWidu3aZTp+zZo1y3R8E7D8HgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAApeb1119XeHi43NzcNGPGDFeXU+HxSjsAAAAAKCc7BtjK9XhNFljlerysrCwNHz5c06dP1z333CO73V6ux78aEeoBAAAAAKUiNTVVFy5cUNeuXRUaGurqcq4KLL8HAAAAADisXLlSN998s6pUqaLAwEB169ZNe/bscWw/ePCg+vTpo2rVqqlSpUpq3ry5Nm/erMTERF177bWSpLp168pms2n//v0uOourB6EeAAAAAOBw5swZxcfHa8uWLVqzZo3c3Nx01113KT8/X6dPn1a7du106NAhLVu2TNu3b9eTTz6p/Px89e7dW//73/8kSd9++63S0tIUHh7u4rOp+Fh+DwAAAABwuOeee5y+z5s3T0FBQdqxY4c2bNigo0ePasuWLapWrZokqX79+o6+gYGBkqQaNWooJCSk/Iq+ijFTDwAAAABw2LNnj/r27au6desqICBAkZGRkn6/Xz45OVk33HCDI9DD9ZipBwAAAAA4dO/eXeHh4XrjjTcUFham/Px8RUdHKycnR76+vq4uD3/CTD0AAAAAQJKUkZGhlJQU/fvf/9Ztt92mqKgonThxwrG9adOmSk5O1vHjx11YJf6IUA8AAAAAkCRVrVpVgYGBev3117V792598cUXio+Pd2y///77FRISop49e2r9+vXau3evli5dqo0bN7qw6qsby+8BAAAAoJw0WWC5uoRLcnNz0+LFizVixAhFR0erUaNGmjlzptq3by9J8vLy0ueff67Ro0erS5cuys3NVZMmTfTaa6+5tvCrGKEeAAAAAOBw++23a8eOHU5tlvV/f4yoXbu2Pvzww0L3vf766536ouyx/B4AAAAAAEMZFeq//vprde/eXWFhYbLZbPrkk0/+cp+vvvpKzZo1k4+Pj+rWras5c+aUfaEAAAAAAJQDo0L9mTNndN111+nVV18tVv99+/apS5cuatu2rbZt26axY8dqxIgRWrp0aRlXCgAAAABA2TPqnvrOnTurc+fOxe4/Z84cRUREaMaMGZKkqKgoJSUl6aWXXtI999xTRlUCAAAAAFA+jJqpL6mNGzcqNjbWqa1Tp05KSkrShQsXCt0nOztbWVlZTh8AAAAAAK5ERs3Ul1R6erqCg4Od2oKDg5Wbm6tjx44pNDS0wD4JCQmaOHFieZUo2+jlJdsh+qsSdf9p7bRi9z1Xc1GJxm42pU+J+gMAAAAASleFnqmXJJvN5vT94usV/tx+0ZgxY5SZmen4HDhwoMxrBAAAAADg76jQM/UhISFKT093ajty5Ig8PDwUGBhY6D7e3t7y9vYuj/IAAAAAALgsFXqmvnXr1lq9erVT2+eff67mzZvL09PTRVUBAAAAAFA6jAr1p0+fVnJyspKTkyX9/sq65ORkpaamSvp96XxcXJyj/9ChQ/Xrr78qPj5eKSkpmj9/vubNm6fHH3/cFeUDAAAAAFCqjFp+n5SUpA4dOji+x8fHS5IGDBigxMREpaWlOQK+JEVGRmrFihUaNWqUXnvtNYWFhWnmzJm8zg4AAACAS2wdu7hcj+eKh1snJiZq5MiROnnyZLkf+2pkVKhv376940F3hUlMTCzQ1q5dO3333XdlWBUAAAAA4Ep04cKFCn/rtVHL7wEAAAAAZWvlypW6+eabVaVKFQUGBqpbt27as2ePJGnt2rWy2WxOs/DJycmy2Wzav3+/1q5dq0GDBikzM1M2m002m00TJkyQJJ04cUJxcXGqWrWq/Pz81LlzZ+3atUuSlJmZKV9fX61cudKplo8++kiVKlXS6dOnJUlPPfWUGjZsKD8/P9WtW1fPPPOMLly44Og/YcIEXX/99Zo/f77q1q0rb2/vS04MVwRGzdQDpcE2ennJdoj+qkTdrUEvlWx8AAAA4Apy5swZxcfH69prr9WZM2c0btw43XXXXY5nm11KTEyMZsyYoXHjxmnnzp2SJH9/f0nSwIEDtWvXLi1btkwBAQF66qmn1KVLF+3YsUN2u11du3bVwoULdccddzjGe++999SjRw/HGJUrV1ZiYqLCwsL0ww8/6KGHHlLlypX15JNPOvbZvXu33n//fS1dulTu7u6leGWuTIR6AAAAAIDDn59BNm/ePAUFBWnHjh1/ua+Xl5fsdrtsNptCQkIc7RfD/Pr16xUTEyNJWrhwocLDw/XJJ5/ovvvuU79+/RQXF6ezZ8/Kz89PWVlZ+vTTT7V06VLHOP/+978d/1ynTh2NHj1aS5YscQr1OTk5euedd1SjRo2/fQ1MwvJ7AAAAAIDDnj171LdvX9WtW1cBAQGKjIyUJKeHkpdUSkqKPDw8dNNNNznaAgMD1ahRI6WkpEiSunbtKg8PDy1btkyStHTpUlWuXFmxsbGOfT788EPdfPPNCgkJkb+/v5555pkCddWuXfuqCfQSM/VAqdsxwFbsvudqLirR2K54eikAAACuLt27d1d4eLjeeOMNhYWFKT8/X9HR0crJyXEsg//jfep/vKe9KEXd125Zlmy23///s5eXl+69916999576tOnj9577z317t1bHh6/x9ZNmzapT58+mjhxojp16iS73a7Fixdr2rRpTmNWqlTpb523qQj1AAAAAABJUkZGhlJSUjR37ly1bdtWkvTNN984tl+cAU9LS1PVqlUlqcC99l5eXsrLy3Nqa9KkiXJzc7V582bH8vuMjAz98ssvioqKcvTr16+fYmNj9dNPP+nLL7/Us88+69i2fv161a5dW08//bSj7ddffy2FszYby+8BAAAAAJKkqlWrKjAwUK+//rp2796tL774QvHx8Y7t9evXV3h4uCZMmKBffvlFn376aYGZ8jp16uj06dNas2aNjh07prNnz6pBgwbq0aOHHnroIX3zzTfavn27HnjgAdWsWVM9evRw7NuuXTsFBwerX79+qlOnjlq1auV07NTUVC1evFh79uzRzJkz9fHHH5f9RbnCMVMPAAAAAOXkSr+d0s3NTYsXL9aIESMUHR2tRo0aaebMmWrfvr0kydPTU4sWLdKjjz6q6667Ti1atNDkyZN13333OcaIiYnR0KFD1bt3b2VkZGj8+PGaMGGC3nrrLf3zn/9Ut27dlJOTo1tuuUUrVqxweo+8zWbT/fffrxdffFHjxo1zqq1Hjx4aNWqUhg8fruzsbHXt2lXPPPOM45V5VyubVdFf2neZsrKyZLfblZmZqYCAgFIfv6xfr/bT2ml/3en/d7Xc3801L3+8RhAAAFxtzp8/r3379ikyMlI+Pj6uLgdXoEv9jpQkh7L8HgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAygjPJUdRSut3g1APAAAAAKXM3d1dkpSTk+PiSnClOnv2rCQ5vdLv7+A99QCMt2OArdh9eY1gEXiNIAAApcrDw0N+fn46evSoPD095ebGfCp+Z1mWzp49qyNHjqhKlSqOPwD9XYR6AAAAAChlNptNoaGh2rdvn3799VdXl4MrUJUqVRQSEnLZ4xDqAQAAAKAMeHl5qUGDBizBRwGenp6XPUN/EaEeAAAAAMqIm5ubfHx8XF0GKjBu7AAAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMxYPyAAAoBbbRy0u2Q/RXJepuDXqpZOMDAICrAjP1AAAAAAAYilAPAAAAAIChWH4PACixHQNsxe57ruaiEo3dbEqfkpYDAABw1WKmHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUD8oDAABGso1eXrIdor8qUXdr0EslGx8AABdgph4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFA/KAwDAADsG2Ird91zNRSUau9mUPiUtBwAAXCGYqQcAAAAAwFCEegAAAAAADMXyewAAABSLbfTyku0Q/VWJuluDXirZ+AAAZuoBAAAAADAVoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAU76kHAAAoxI4BtmL3PVdzUYnGbjalT0nLwVXKNnp5yXaI/qpE3a1BL5VsfABXHGbqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFO+pBwAAAID/n2308pLtEP1Vibpbg14q2fjAX2CmHgAAAAAAQxHqAQAAAAAwFMvvAQAAcEXYMcBW7L7nai4q0djNpvQpaTkAYARm6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwlHGhftasWYqMjJSPj4+aNWumdevWXbL/woULdd1118nPz0+hoaEaNGiQMjIyyqlaAAAAAADKjlGhfsmSJRo5cqSefvppbdu2TW3btlXnzp2VmppaaP9vvvlGcXFxGjx4sH766Sd98MEH2rJli4YMGVLOlQMAAAAAUPqMevr99OnTNXjwYEconzFjhlatWqXZs2crISGhQP9NmzapTp06GjFihCQpMjJSjzzyiF544YVyrRsAAAAAUDjb6OUl2yH6qxJ1twa9VLLxDWPMTH1OTo62bt2q2NhYp/bY2Fht2LCh0H1iYmJ08OBBrVixQpZl6fDhw/rwww/VtWvXIo+TnZ2trKwspw8AAAAAAFciY0L9sWPHlJeXp+DgYKf24OBgpaenF7pPTEyMFi5cqN69e8vLy0shISGqUqWKXnnllSKPk5CQILvd7viEh4eX6nkAAAAAAFBajAn1F9lsNqfvlmUVaLtox44dGjFihMaNG6etW7dq5cqV2rdvn4YOHVrk+GPGjFFmZqbjc+DAgVKtHwAAAACA0mLMPfXVq1eXu7t7gVn5I0eOFJi9vyghIUFt2rTRE088IUlq2rSpKlWqpLZt22ry5MkKDQ0tsI+3t7e8vb1L/wQAAACAK8yOAYVPjhXmXM1FJRq72ZQ+JS0HwN9gzEy9l5eXmjVrptWrVzu1r169WjExMYXuc/bsWbm5OZ+iu7u7pN9n+AEAAAAAMJkxoV6S4uPj9eabb2r+/PlKSUnRqFGjlJqa6lhOP2bMGMXFxTn6d+/eXR999JFmz56tvXv3av369RoxYoRatmypsLAwV50GAAAAAAClwpjl95LUu3dvZWRkaNKkSUpLS1N0dLRWrFih2rVrS5LS0tKc3lk/cOBAnTp1Sq+++qpGjx6tKlWq6NZbb9Xzzz/vqlMAAAAAAKDUGBXqJWnYsGEaNmxYodsSExMLtP3jH//QP/7xjzKuCgAAAACA8mfU8nsAAAAAAPB/CPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYyrhX2gEAAACAqXYMsBW777mai0o0drMpfUpaDioAZuoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQ3m4ugAAAAAAAMrKjgG2Yvc9V3NRicZuNqVPScspdczUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKONC/axZsxQZGSkfHx81a9ZM69atu2T/7OxsPf3006pdu7a8vb1Vr149zZ8/v5yqBQAAAACg7Hi4uoCSWLJkiUaOHKlZs2apTZs2mjt3rjp37qwdO3YoIiKi0H169eqlw4cPa968eapfv76OHDmi3Nzccq4cAAAAAIDSZ1Sonz59ugYPHqwhQ4ZIkmbMmKFVq1Zp9uzZSkhIKNB/5cqV+uqrr7R3715Vq1ZNklSnTp3yLBkAAAAAgDJjzPL7nJwcbd26VbGxsU7tsbGx2rBhQ6H7LFu2TM2bN9cLL7ygmjVrqmHDhnr88cd17ty5Io+TnZ2trKwspw8AAAAAAFciY2bqjx07pry8PAUHBzu1BwcHKz09vdB99u7dq2+++UY+Pj76+OOPdezYMQ0bNkzHjx8v8r76hIQETZw4sdTrBwAAAACgtBkzU3+RzWZz+m5ZVoG2i/Lz82Wz2bRw4UK1bNlSXbp00fTp05WYmFjkbP2YMWOUmZnp+Bw4cKDUzwEAAAAAgNJgzEx99erV5e7uXmBW/siRIwVm7y8KDQ1VzZo1ZbfbHW1RUVGyLEsHDx5UgwYNCuzj7e0tb2/v0i0eAAAAAIAyYMxMvZeXl5o1a6bVq1c7ta9evVoxMTGF7tOmTRsdOnRIp0+fdrT98ssvcnNzU61atcq0XgAAAAAAypoxoV6S4uPj9eabb2r+/PlKSUnRqFGjlJqaqqFDh0r6fel8XFyco3/fvn0VGBioQYMGaceOHfr666/1xBNP6MEHH5Svr6+rTgMAAAAAgFJhzPJ7Serdu7cyMjI0adIkpaWlKTo6WitWrFDt2rUlSWlpaUpNTXX09/f31+rVq/WPf/xDzZs3V2BgoHr16qXJkye76hQAAAAAACg1RoV6SRo2bJiGDRtW6LbExMQCbY0bNy6wZB8AAAAAgIrAqOX3AAAAAADg/xDqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwFKEeAAAAAABDEeoBAAAAADAUoR4AAAAAAEMR6gEAAAAAMBShHgAAAAAAQxHqAQAAAAAwlHGhftasWYqMjJSPj4+aNWumdevWFWu/9evXy8PDQ9dff33ZFggAAAAAQDkxKtQvWbJEI0eO1NNPP61t27apbdu26ty5s1JTUy+5X2ZmpuLi4nTbbbeVU6UAAAAAAJQ9o0L99OnTNXjwYA0ZMkRRUVGaMWOGwsPDNXv27Evu98gjj6hv375q3bp1OVUKAAAAAEDZMybU5+TkaOvWrYqNjXVqj42N1YYNG4rc76233tKePXs0fvz4Yh0nOztbWVlZTh8AAAAAAK5ExoT6Y8eOKS8vT8HBwU7twcHBSk9PL3SfXbt26V//+pcWLlwoDw+PYh0nISFBdrvd8QkPD7/s2gEAAAAAKAt/O9Tv3r1bq1at0rlz5yRJlmWVWlGXYrPZnL5bllWgTZLy8vLUt29fTZw4UQ0bNiz2+GPGjFFmZqbjc+DAgcuuGQAAAACAslC86es/yMjIUO/evfXFF1/IZrNp165dqlu3roYMGaIqVapo2rRpZVGnqlevLnd39wKz8keOHCkwey9Jp06dUlJSkrZt26bhw4dLkvLz82VZljw8PPT555/r1ltvLbCft7e3vL29y+QcAAAAAAAoTSWeqR81apQ8PDyUmpoqPz8/R3vv3r21cuXKUi3uj7y8vNSsWTOtXr3aqX316tWKiYkp0D8gIEA//PCDkpOTHZ+hQ4eqUaNGSk5O1k033VRmtQIAAAAAUB5KPFP/+eefa9WqVapVq5ZTe4MGDfTrr7+WWmGFiY+PV//+/dW8eXO1bt1ar7/+ulJTUzV06FBJvy+d/+233/T222/Lzc1N0dHRTvsHBQXJx8enQDsAAAAAACYqcag/c+aM0wz9RceOHSvzZeu9e/dWRkaGJk2apLS0NEVHR2vFihWqXbu2JCktLe0v31kPAAAAAEBFUeLl97fccovefvttx3ebzab8/Hy9+OKL6tChQ6kWV5hhw4Zp//79ys7O1tatW3XLLbc4tiUmJmrt2rVF7jthwgQlJyeXeY0AAAAAAJSHEs/Uv/jii2rfvr2SkpKUk5OjJ598Uj/99JOOHz+u9evXl0WNAAAAAACgECWeqW/SpIm+//57tWzZUh07dtSZM2d09913a9u2bapXr15Z1AgAAAAAAApR4pl6SQoJCdHEiRNLuxYAAAAAAFACJQ71X3/99SW3//EedwAAAAAAUHZKHOrbt29foM1mszn+OS8v77IKAgAAAAAAxVPie+pPnDjh9Dly5IhWrlypFi1a6PPPPy+LGgEAAAAAQCFKPFNvt9sLtHXs2FHe3t4aNWqUtm7dWiqFAQAAAACASyvxTH1RatSooZ07d5bWcAAAAAAA4C+UeKb++++/d/puWZbS0tI0depUXXfddaVWGAAAAAAAuLQSh/rrr79eNptNlmU5tbdq1Urz588vtcIAAAAAAMCllTjU79u3z+m7m5ubatSoIR8fn1IrCgAAAAAA/LUSh/ratWuXRR0AAAAAAKCEihXqZ86cWewBR4wY8beLAQAAAAAAxVesUP+f//ynWIPZbDZCPQAAAAAA5aRYof7P99EDAAAAAADXK7X31AMAAAAAgPJV4gflSdLBgwe1bNkypaamKicnx2nb9OnTS6UwAAAAAABwaSUO9WvWrNGdd96pyMhI7dy5U9HR0dq/f78sy9KNN95YFjUCAAAAAIBClHj5/ZgxYzR69Gj9+OOP8vHx0dKlS3XgwAG1a9dO9913X1nUCAAAAAAAClHiUJ+SkqIBAwZIkjw8PHTu3Dn5+/tr0qRJev7550u9QAAAAAAAULgSh/pKlSopOztbkhQWFqY9e/Y4th07dqz0KgMAAAAAAJdU4nvqW7VqpfXr16tJkybq2rWrRo8erR9++EEfffSRWrVqVRY1AgAAAACAQpQ41E+fPl2nT5+WJE2YMEGnT5/WkiVLVL9+ff3nP/8p9QIBAAAAAEDhShzqn332WT3wwAOyLEt+fn6aNWtWWdQFAAAAAAD+Qonvqc/IyFDXrl1Vq1YtjR49WsnJyWVQFgAAAAAA+CslDvXLli1Tenq6xo8fr61bt6pZs2Zq0qSJpkyZov3795dBiQAAAAAAoDAlDvWSVKVKFT388MNau3atfv31Vw0aNEjvvPOO6tevX9r1AQAAAACAIvytUH/RhQsXlJSUpM2bN2v//v0KDg4urboAAAAAAMBf+Fuh/ssvv9RDDz2k4OBgDRgwQJUrV9by5ct14MCB0q4PAAAAAAAUocRPv69Vq5YyMjLUqVMnzZ07V927d5ePj09Z1AYAAAAAAC6hxKF+3Lhxuu+++1S1atWyqAcAAAAAABRTiUP9ww8/XBZ1AAAAAACAErqsB+UBAAAAAADXIdQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwAAAABgKEI9AAAAAACGItQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhjIu1M+aNUuRkZHy8fFRs2bNtG7duiL7fvTRR+rYsaNq1KihgIAAtW7dWqtWrSrHagEAAAAAKDtGhfolS5Zo5MiRevrpp7Vt2za1bdtWnTt3VmpqaqH9v/76a3Xs2FErVqzQ1q1b1aFDB3Xv3l3btm0r58oBAAAAACh9RoX66dOna/DgwRoyZIiioqI0Y8YMhYeHa/bs2YX2nzFjhp588km1aNFCDRo00JQpU9SgQQMtX768nCsHAAAAAKD0GRPqc3JytHXrVsXGxjq1x8bGasOGDcUaIz8/X6dOnVK1atWK7JOdna2srCynDwAAAAAAVyJjQv2xY8eUl5en4OBgp/bg4GClp6cXa4xp06bpzJkz6tWrV5F9EhISZLfbHZ/w8PDLqhsAAAAAgLJiTKi/yGazOX23LKtAW2EWLVqkCRMmaMmSJQoKCiqy35gxY5SZmen4HDhw4LJrBgAAAACgLHi4uoDiql69utzd3QvMyh85cqTA7P2fLVmyRIMHD9YHH3yg22+//ZJ9vb295e3tfdn1AgAAAABQ1oyZqffy8lKzZs20evVqp/bVq1crJiamyP0WLVqkgQMH6r333lPXrl3LukwAAAAAAMqNMTP1khQfH6/+/furefPmat26tV5//XWlpqZq6NChkn5fOv/bb7/p7bfflvR7oI+Li9PLL7+sVq1aOWb5fX19ZbfbXXYeAAAAAACUBqNCfe/evZWRkaFJkyYpLS1N0dHRWrFihWrXri1JSktLc3pn/dy5c5Wbm6vHHntMjz32mKN9wIABSkxMLO/yAQAAAAAoVUaFekkaNmyYhg0bVui2Pwf1tWvXln1BAAAAAAC4iDH31AMAAAAAAGeEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADGVcqJ81a5YiIyPl4+OjZs2aad26dZfs/9VXX6lZs2by8fFR3bp1NWfOnHKqFAAAAACAsmVUqF+yZIlGjhypp59+Wtu2bVPbtm3VuXNnpaamFtp/37596tKli9q2batt27Zp7NixGjFihJYuXVrOlQMAAAAAUPqMCvXTp0/X4MGDNWTIEEVFRWnGjBkKDw/X7NmzC+0/Z84cRUREaMaMGYqKitKQIUP04IMP6qWXXirnygEAAAAAKH3GhPqcnBxt3bpVsbGxTu2xsbHasGFDofts3LixQP9OnTopKSlJFy5cKHSf7OxsZWVlOX0AAAAAALgS2SzLslxdRHEcOnRINWvW1Pr16xUTE+NonzJlihYsWKCdO3cW2Kdhw4YaOHCgxo4d62jbsGGD2rRpo0OHDik0NLTAPhMmTNDEiRMLtGdmZiogIKCUzgYAAAAAgMJlZWXJbrcXK4caM1N/kc1mc/puWVaBtr/qX1j7RWPGjFFmZqbjc+DAgcusGAAAAACAsuHh6gKKq3r16nJ3d1d6erpT+5EjRxQcHFzoPiEhIYX29/DwUGBgYKH7eHt7y9vbu3SKBgAAAACgDBkzU+/l5aVmzZpp9erVTu2rV692Wo7/R61bty7Q//PPP1fz5s3l6elZZrUCAAAAAFAejAn1khQfH68333xT8+fPV0pKikaNGqXU1FQNHTpU0u9L5+Pi4hz9hw4dql9//VXx8fFKSUnR/PnzNW/ePD3++OOuOgUAAAAAAEqNMcvvJal3797KyMjQpEmTlJaWpujoaK1YsUK1a9eWJKWlpTm9sz4yMlIrVqzQqFGj9NprryksLEwzZ87UPffc46pTAAAAAACg1Bjz9HtXKclTBwEAAAAAuFwV+un3AAAAAADgd4R6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUMaE+hMnTqh///6y2+2y2+3q37+/Tp48WWT/Cxcu6KmnntK1116rSpUqKSwsTHFxcTp06FD5FQ0AAAAAQBkyJtT37dtXycnJWrlypVauXKnk5GT179+/yP5nz57Vd999p2eeeUbfffedPvroI/3yyy+68847y7FqAAAAAADKjs2yLMvVRfyVlJQUNWnSRJs2bdJNN90kSdq0aZNat26tn3/+WY0aNSrWOFu2bFHLli3166+/KiIiolj7ZGVlyW63KzMzUwEBAX/7HAAAAAAAKI6S5FAjZuo3btwou93uCPSS1KpVK9ntdm3YsKHY42RmZspms6lKlSpF9snOzlZWVpbTBwAAAACAK5ERoT49PV1BQUEF2oOCgpSenl6sMc6fP69//etf6tu37yX/0pGQkOC4b99utys8PPxv1w0AAAAAQFlyaaifMGGCbDbbJT9JSUmSJJvNVmB/y7IKbf+zCxcuqE+fPsrPz9esWbMu2XfMmDHKzMx0fA4cOPD3Tg4AAAAAgDLm4cqDDx8+XH369Llknzp16uj777/X4cOHC2w7evSogoODL7n/hQsX1KtXL+3bt09ffPHFX96P4O3tLW9v778uHgAAAAAAF3NpqK9evbqqV6/+l/1at26tzMxMffvtt2rZsqUkafPmzcrMzFRMTEyR+10M9Lt27dKXX36pwMDAUqsdAAAAAABXM+Ke+qioKN1xxx166KGHtGnTJm3atEkPPfSQunXr5vTk+8aNG+vjjz+WJOXm5uree+9VUlKSFi5cqLy8PKWnpys9PV05OTmuOhUAAAAAAEqNEaFekhYuXKhrr71WsbGxio2NVdOmTfXOO+849dm5c6cyMzMlSQcPHtSyZct08OBBXX/99QoNDXV8SvLEfAAAAAAArlRGvKfelXhPPQAAAACgPFW499QDAAAAAICCCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABjKmFB/4sQJ9e/fX3a7XXa7Xf3799fJkyeLvf8jjzwim82mGTNmlFmNAAAAAACUJ2NCfd++fZWcnKyVK1dq5cqVSk5OVv/+/Yu17yeffKLNmzcrLCysjKsEAAAAAKD8eLi6gOJISUnRypUrtWnTJt10002SpDfeeEOtW7fWzp071ahRoyL3/e233zR8+HCtWrVKXbt2La+SAQAAAAAoc0bM1G/cuFF2u90R6CWpVatWstvt2rBhQ5H75efnq3///nriiSd0zTXXFOtY2dnZysrKcvoAAAAAAHAlMiLUp6enKygoqEB7UFCQ0tPTi9zv+eefl4eHh0aMGFHsYyUkJDju27fb7QoPD/9bNQMAAAAAUNZcGuonTJggm812yU9SUpIkyWazFdjfsqxC2yVp69atevnll5WYmFhkn8KMGTNGmZmZjs+BAwf+3skBAAAAAFDGXHpP/fDhw9WnT59L9qlTp46+//57HT58uMC2o0ePKjg4uND91q1bpyNHjigiIsLRlpeXp9GjR2vGjBnav39/oft5e3vL29u7+CcBAAAAAICLuDTUV69eXdWrV//Lfq1bt1ZmZqa+/fZbtWzZUpK0efNmZWZmKiYmptB9+vfvr9tvv92prVOnTurfv78GDRp0+cUDAAAAAOBiRjz9PioqSnfccYceeughzZ07V5L08MMPq1u3bk5Pvm/cuLESEhJ01113KTAwUIGBgU7jeHp6KiQk5JJPywcAAAAAwBRGPChPkhYuXKhrr71WsbGxio2NVdOmTfXOO+849dm5c6cyMzNdVCEAAAAAAOXLZlmW5eoirmRZWVmy2+3KzMxUQECAq8sBAAAAAFRwJcmhxszUAwAAAAAAZ4R6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAAAAwFCEegAAAAAADEWoBwAAAADAUIR6AAAAAAAMRagHAAAAAMBQhHoAAAAAAAzl4eoCrnSWZUmSsrKyXFwJAAAAAOBqcDF/Xsyjl0Ko/wunTp2SJIWHh7u4EgAAAADA1eTUqVOy2+2X7GOzihP9r2L5+fk6dOiQKleuLJvN5upySiQrK0vh4eE6cOCAAgICXF3OVYFrXv645uWPa17+uOblj2te/rjm5Y9rXv645uXP1GtuWZZOnTqlsLAwubld+q55Zur/gpubm2rVquXqMi5LQECAUb/AFQHXvPxxzcsf17z8cc3LH9e8/HHNyx/XvPxxzcufidf8r2boL+JBeQAAAAAAGIpQDwAAAACAoQj1FZi3t7fGjx8vb29vV5dy1eCalz+uefnjmpc/rnn545qXP655+eOalz+uefm7Gq45D8oDAAAAAMBQzNQDAAAAAGAoQj0AAAAAAIYi1AMAAAAAYChCPQAAAAAAhiLUAwBwheOZtgAAoCgeri4ApePgwYOaPXu2NmzYoPT0dNlsNgUHBysmJkZDhw5VeHi4q0sEAPxN3t7e2r59u6KiolxdCgAAV6y0tDTNnj1b33zzjdLS0uTu7q7IyEj17NlTAwcOlLu7u6tLLBO80q4C+Oabb9S5c2eFh4crNjZWwcHBsixLR44c0erVq3XgwAF99tlnatOmjatLvaocOHBA48eP1/z5811dSoVy7tw5bd26VdWqVVOTJk2ctp0/f17vv/++4uLiXFRdxZSSkqJNmzapdevWaty4sX7++We9/PLLys7O1gMPPKBbb73V1SVWGPHx8YW2v/zyy3rggQcUGBgoSZo+fXp5lnXVOXHihBYsWKBdu3YpNDRUAwYM4I/jpWzbtm2qUqWKIiMjJUnvvvuuZs+erdTUVNWuXVvDhw9Xnz59XFxlxfKPf/xDvXr1Utu2bV1dylXllVdeUVJSkrp27apevXrpnXfeUUJCgvLz83X33Xdr0qRJ8vBgnrU0JCUl6fbbb1dkZKR8fX21efNm9evXTzk5OVq1apWioqK0atUqVa5c2dWllj4LxmvevLk1cuTIIrePHDnSat68eTlWBMuyrOTkZMvNzc3VZVQoO3futGrXrm3ZbDbLzc3NateunXXo0CHH9vT0dK55Kfvss88sLy8vq1q1apaPj4/12WefWTVq1LBuv/1267bbbrM8PDysNWvWuLrMCsNms1nXX3+91b59e6ePzWazWrRoYbVv397q0KGDq8uscEJDQ61jx45ZlmVZe/futUJCQqyQkBCrY8eOVq1atSy73W6lpKS4uMqK5YYbbrC++OILy7Is64033rB8fX2tESNGWLNnz7ZGjhxp+fv7W/PmzXNxlRXLxf92NmjQwJo6daqVlpbm6pIqvEmTJlmVK1e27rnnHiskJMSaOnWqFRgYaE2ePNmaMmWKVaNGDWvcuHGuLrPCaNOmjTVhwgTH93feece66aabLMuyrOPHj1vXX3+9NWLECFeVV6aYqa8AfH19lZycrEaNGhW6/eeff9YNN9ygc+fOlXNlFduyZcsuuX3v3r0aPXq08vLyyqmiiu+uu+5Sbm6u3nrrLZ08eVLx8fH68ccftXbtWkVEROjw4cMKCwvjmpeimJgY3XrrrZo8ebIWL16sYcOG6dFHH9Vzzz0nSXr66ae1ZcsWff755y6utGJISEjQG2+8oTfffNNpBYSnp6e2b99eYHUKSoebm5vS09MVFBSk+++/X+np6fr000/l5+en7Oxs3XvvvfLx8dEHH3zg6lIrjEqVKiklJUURERG68cYbNXToUD388MOO7e+9956ee+45/fTTTy6ssmJxc3PT6tWrtXz5ci1cuFCZmZnq3LmzHnroIXXp0kVubjxqq7TVq1dPL774ou6++25t375dzZo104IFC9SvXz9J0scff6wnn3xSu3btcnGlFYOfn59+/PFH1a1bV5KUn58vHx8fHThwQMHBwVq9erUGDhyo3377zcWVlgFX/1UBly8yMtKaP39+kdvnz59vRUZGlmNFV4eLf/G22WxFfpg1Ll1BQUHW999/79Q2bNgwKyIiwtqzZw8z9WUgICDA2rVrl2VZlpWXl2d5eHhYW7dudWz/4YcfrODgYFeVVyF9++23VsOGDa3Ro0dbOTk5lmVZloeHh/XTTz+5uLKKy2azWYcPH7Ys6/f/pv559cmmTZusWrVquaK0CiswMNBKSkqyLOv3f7cnJyc7bd+9e7fl6+vritIqrD/+nufk5FhLliyxOnXqZLm7u1thYWHW2LFjHf++R+nw9fW1fv31V8d3T09P68cff3R8379/v+Xn5+eK0iqk2rVrW998843j+6FDhyybzWadPXvWsizL2rdvn+Xj4+Oq8soUf5KrAB5//HENHTpUw4cP13//+19t2rRJmzdv1n//+18NHz5cjz76qJ588klXl1nhhIaGaunSpcrPzy/0891337m6xArn3LlzBe47e+2113TnnXeqXbt2+uWXX1xU2dXBzc1NPj4+qlKliqOtcuXKyszMdF1RFVCLFi20detWHT16VM2bN9cPP/wgm83m6rIqvIvXODs7W8HBwU7bgoODdfToUVeUVWF17txZs2fPliS1a9dOH374odP2999/X/Xr13dFaVcFT09P9erVSytXrtTevXv10EMPaeHChUWu+sTfExISoh07dkiSdu3apby8PMd3Sfrpp58UFBTkqvIqnJ49e2ro0KFauXKlvvzyS/Xr10/t2rWTr6+vJGnnzp2qWbOmi6ssGzyVoQIYNmyYAgMD9Z///Edz5851LD12d3dXs2bN9Pbbb6tXr14urrLiadasmb777jv17Nmz0O02m43XUJWyxo0bKykpqcATwF955RVZlqU777zTRZVVXHXq1NHu3bsd/+d648aNioiIcGw/cOCAQkNDXVVeheXv768FCxZo8eLF6tixI7eUlIPbbrtNHh4eysrK0i+//KJrrrnGsS01NVXVq1d3YXUVz/PPP682bdqoXbt2at68uaZNm6a1a9cqKipKO3fu1KZNm/Txxx+7usyrQkREhCZMmKDx48frf//7n6vLqVD69u2ruLg49ejRQ2vWrNFTTz2lxx9/XBkZGbLZbHruued07733urrMCmPy5MlKS0tT9+7dlZeXp9atW+vdd991bLfZbEpISHBhhWWHe+ormAsXLujYsWOSpOrVq8vT09PFFVVc69at05kzZ3THHXcUuv3MmTNKSkpSu3btyrmyiishIUHr1q3TihUrCt0+bNgwzZkzR/n5+eVcWcU1Z84chYeHq2vXroVuf/rpp3X48GG9+eab5VzZ1ePgwYPaunWrbr/9dlWqVMnV5VRIEydOdPreqlUrderUyfH9iSee0MGDB7Vo0aLyLq1CO3nypKZOnarly5dr7969ys/PV2hoqNq0aaNRo0apefPmri6xQomMjFRSUpLjLRooe3l5eZo6dao2bdqkm2++WU899ZQWL16sJ598UmfPnlX37t316quv8u/2Unb+/Hnl5ubK39/f1aWUG0I9AAAAAACG4p56AAAAAAAMRagHAAAAAMBQhHoAAAAAAAxFqAcAAH9L+/btNXLkSFeXAQDAVY1QDwAAAACAoQj1AAAAAAAYilAPAAAu27vvvqvmzZurcuXKCgkJUd++fXXkyBGnPsuWLVODBg3k6+urDh06aMGCBbLZbDp58qRrigYAoAIg1AMAgMuWk5OjZ599Vtu3b9cnn3yiffv2aeDAgY7t+/fv17333quePXsqOTlZjzzyiJ5++mnXFQwAQAXh4eoCAACA+R588EHHP9etW1czZ85Uy5Ytdfr0afn7+2vOnDlq1KiRXnzxRUlSo0aN9OOPP+q5555zVckAAFQIzNQDAIDLtm3bNvXo0UO1a9dW5cqV1b59e0lSamqqJGnnzp1q0aKF0z4tW7Ys7zIBAKhwCPUAAOCynDlzRrGxsfL399e7776rLVu26OOPP5b0+7J8SbIsSzabzWk/y7LKvVYAACoalt8DAIDL8vPPP+vYsWOaOnWqwsPDJUlJSUlOfRo3bqwVK1Y4tf25DwAAKDlm6gEAwGWJiIiQl5eXXnnlFe3du1fLli3Ts88+69TnkUce0c8//6ynnnpKv/zyi95//30lJiZKUoEZfAAAUHyEegAAcFlq1KihxMREffDBB2rSpImmTp2ql156yalPZGSkPvzwQ3300Udq2rSpZs+e7Xj6vbe3tyvKBgCgQrBZ3NAGAABc4LnnntOcOXN04MABV5cCAICxuKceAACUi1mzZqlFixYKDAzU+vXr9eKLL2r48OGuLgsAAKMR6gEAQLnYtWuXJk+erOPHjysiIkKjR4/WmDFjXF0WAABGY/k9AAAAAACG4kF5AAAAAAAYilAPAAAAAIChCPUAAAAAABiKUA8AAAAAgKEI9QAAAAAAGIpQDwAAAACAoQj1AAAAAAAYilAPAAAAAIChCPUAAAAAABjq/wPCpx/lA5Ie8QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def autocovariance(x, lag):\n",
    "    \"\"\"\n",
    "    计算给定滞后期下的自协方差\n",
    "    \"\"\"\n",
    "    n = len(x)\n",
    "    mean_x = np.mean(x)\n",
    "    # 计算自协方差\n",
    "    cov = np.sum((x[:n-lag] - mean_x) * (x[lag:] - mean_x)) / n\n",
    "    return cov\n",
    "\n",
    "def autocorrelation_from_cov(x, lag):\n",
    "    \"\"\"\n",
    "    通过自协方差计算自相关系数\n",
    "    \"\"\"\n",
    "    var_x = np.var(x)  # 计算方差\n",
    "    cov_lag = autocovariance(x, lag)  # 计算自协方差\n",
    "    return cov_lag / var_x  # 标准化自协方差，得到自相关系数\n",
    "\n",
    "\n",
    "# 自相关系数计算方法\n",
    "def autocorrelation2(x, max_lag):\n",
    "    n = len(x)\n",
    "    mean_x = np.mean(x)  # 均值\n",
    "    s = np.std(x, ddof=0)  # 标准差\n",
    "\n",
    "    autocorr = np.zeros(max_lag + 1)  \n",
    "\n",
    "    for j in range(max_lag + 1):\n",
    "        autocorr[j] = np.sum((x[:n-j] - mean_x) * (x[j:] - mean_x)) / ((n - j) * s**2)\n",
    "\n",
    "    return autocorr\n",
    "\n",
    "\n",
    "\n",
    "def acf_by_hand(x, lag):\n",
    "    y1 = x[:(len(x)-lag)]\n",
    "    y2 = x[lag:]\n",
    "    sum_product = np.sum((y1-np.mean(x))*(y2-np.mean(x)))\n",
    "    return sum_product / ((len(x) - lag) * np.var(x))\n",
    "\n",
    "x = np.linspace(0, 10, 101)\n",
    "\n",
    "results = {}\n",
    "lags = 10\n",
    "\n",
    "# 记录每个方法的计算时间\n",
    "times = {}\n",
    "\n",
    "# 比较每个方法的自相关系数计算时间\n",
    "start_time = time.perf_counter()\n",
    "results[\"auto2\"] = [autocorrelation2(x, lag)[lag] for lag in range(lags)]\n",
    "times[\"auto2\"] = time.perf_counter() - start_time\n",
    "\n",
    "start_time = time.perf_counter()\n",
    "results[\"acf_by_hand\"] = [acf_by_hand(x, lag) for lag in range(lags)]\n",
    "times[\"acf_by_hand\"] = time.perf_counter() - start_time\n",
    "\n",
    "\n",
    "start_time = time.perf_counter()\n",
    "results[\"acf\"] = acf(x, nlags=9)\n",
    "times[\"acf\"] = time.perf_counter() - start_time\n",
    "\n",
    "results[\"autovar\"] = [autocorrelation_from_cov(x, lag) for lag in range(lags)]\n",
    "\n",
    "# 打印每个方法的计算时间\n",
    "for method, elapsed_time in times.items():\n",
    "    print(f\"Time taken by {method}: {elapsed_time:.6f} seconds\")\n",
    "\n",
    "# 对比不同方法的计算结果\n",
    "plt.figure(dpi=200)\n",
    "pd.DataFrame(results).plot(kind=\"bar\", figsize=(12,8), )\n",
    "plt.xlabel(\"lag\")\n",
    "plt.ylim([-.5, 1.2])\n",
    "plt.ylabel(\"value\")\n",
    "plt.title(\"Autocorrelation Comparison\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'auto2': [1.0,\n",
       "  0.9998000000000002,\n",
       "  0.9995999600079986,\n",
       "  0.9993998800239955,\n",
       "  0.9991997600479907,\n",
       "  0.9989996000799842,\n",
       "  0.998799400119976,\n",
       "  0.9985991601679666,\n",
       "  0.9983988802239555,\n",
       "  0.9981985602879425],\n",
       " 'acf_by_hand': [1.0,\n",
       "  0.9998000000000002,\n",
       "  0.9995999600079986,\n",
       "  0.9993998800239955,\n",
       "  0.9991997600479907,\n",
       "  0.9989996000799842,\n",
       "  0.998799400119976,\n",
       "  0.9985991601679666,\n",
       "  0.9983988802239555,\n",
       "  0.9981985602879425],\n",
       " 'acf': array([1.        , 0.99970003, 0.99940006, 0.99910009, 0.99880012,\n",
       "        0.99850015, 0.99820018, 0.99790021, 0.99760024, 0.99730027]),\n",
       " 'autovar': [1.0,\n",
       "  0.9997000299970004,\n",
       "  0.9994000600059972,\n",
       "  0.9991000900389868,\n",
       "  0.9988001201079655,\n",
       "  0.9985001502249298,\n",
       "  0.9982001804018759,\n",
       "  0.9979002106508007,\n",
       "  0.9976002409837004,\n",
       "  0.9973002714125709]}"
      ]
     },
     "execution_count": 186,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geocat",
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}