Trends-totPrecip.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6533af15-0b31-489e-80f6-94d8ef65b9a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/global/homes/d/duan0000/.conda/envs/duan/lib/python3.8/site-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.5\n",
      "  warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion}\"\n"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "import xarray as xa\n",
    "import cartopy\n",
    "from cartopy import crs as ccrs\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import cm, colors\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib as mpl\n",
    "import scipy.stats\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "47be826f-5710-459a-a0d3-87b7514c3219",
   "metadata": {},
   "outputs": [],
   "source": [
    "def start_end_time(sub_id):\n",
    "    start = []\n",
    "    end = []\n",
    "    ext_times = ext_days[sub_id]\n",
    "    length = len(ext_times)\n",
    "    i = 0\n",
    "    while i <= length - 1:\n",
    "        j = 1\n",
    "        while (j + i <= length - 1) and ((ext_times[i + j] - ext_times[i + j - 1]) / np.timedelta64(1, \"D\") == 1):\n",
    "            j += 1\n",
    "        start.append(ext_times[i] - np.timedelta64(12, \"h\"))\n",
    "        end.append(ext_times[i + j - 1] + np.timedelta64(12, \"h\"))\n",
    "        i += j\n",
    "    return start, end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6f4a4be9-b3f0-49fc-9397-bf60f479618c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sub_id:  1\n",
      "39 21.363590636973544\n",
      "sub_id:  2\n",
      "38 7.98353434112839\n",
      "sub_id:  3\n",
      "37 7.6309751208430425\n",
      "sub_id:  4\n",
      "77 6.2968948209451305\n",
      "sub_id:  5\n",
      "118 6.469523059079158\n",
      "sub_id:  6\n",
      "124 10.647874758983226\n",
      "sub_id:  7\n",
      "153 17.252742222601025\n"
     ]
    }
   ],
   "source": [
    "def is_monsoon_precip(month):\n",
    "    return (month <= 9) & (month >= 6)\n",
    "data = xa.open_dataarray('cpc-global/NAM_sub_precip', engine='netcdf4')\n",
    "monsoon_precip = data.sel(time=(is_monsoon_precip(data.time.dt.month)))\n",
    "monsoon_precip = monsoon_precip.sel(time=(monsoon_precip.time.dt.year < 2019))\n",
    "ext_days = {}\n",
    "for sub_id in range(1, 8):\n",
    "    print(\"sub_id: \", sub_id)\n",
    "    precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "    precip_data = precip.data\n",
    "    q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "    q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "    ext_time = precip.where(precip > q2, drop=True).time.data\n",
    "    print(len(ext_time), q2)\n",
    "    ext_days[sub_id] = ext_time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f03761b-4cff-424c-b9fe-d2e3aab0062b",
   "metadata": {},
   "source": [
    "# Number of Event days\n",
    "* number of epe days, instead of number of epes. (consecutive days are counted multiple times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dd2b9a82-bbef-4224-a747-4f2a9fa644e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mk_score(x, n):\n",
    "    s = 0\n",
    "    demo = np.ones(n)\n",
    "    for k in range(n-1):\n",
    "        s = s+np.sum(demo[k+1:n][x[k+1:n] > x[k]])-np.sum(demo[k+1:n][x[k+1:n] < x[k]])\n",
    "    return s\n",
    "\n",
    "def variance_s(x, n):\n",
    "    unique_x = np.unique(x)\n",
    "    g = len(unique_x)\n",
    "    if n==g:\n",
    "        var_s = (n*(n-1)*(2*n+5))/18\n",
    "    else:\n",
    "        tp = np.zeros(unique_x.shape)\n",
    "        demo = np.ones(n)\n",
    "        for i in range(g):\n",
    "            tp[i] = np.sum(demo[x==unique_x[i]])\n",
    "        var_s = (n*(n-1)*(2*n+5) - np.sum(tp*(tp-1)*(2*tp+5)))/18\n",
    "    return var_s\n",
    "\n",
    "def z_score(s, var_s):\n",
    "    if s > 0:\n",
    "        z = (s - 1)/np.sqrt(var_s)\n",
    "    elif s == 0:\n",
    "        z = 0\n",
    "    elif s < 0:\n",
    "        z = (s + 1)/np.sqrt(var_s)\n",
    "    return z\n",
    "\n",
    "def p_value(z, alpha):\n",
    "    # two tail test\n",
    "    p = 2*(1-norm.cdf(abs(z)))  \n",
    "    h = abs(z) > norm.ppf(1-alpha/2)\n",
    "\n",
    "    if (z < 0) and h:\n",
    "        trend = 'decreasing'\n",
    "    elif (z > 0) and h:\n",
    "        trend = 'increasing'\n",
    "    else:\n",
    "        trend = 'no trend'\n",
    "    \n",
    "    return p, h, trend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "033c7c63-327a-4d82-8af1-3aa4c5c8e6d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sens_slope(x):\n",
    "    \"\"\"\n",
    "    This method proposed by Theil (1950) and Sen (1968) to estimate the magnitude of the monotonic trend. Intercept calculated using Conover, W.J. (1980) method.\n",
    "    Input:\n",
    "        x:   a one dimensional vector (list, numpy array or pandas series) data\n",
    "    Output:\n",
    "        slope: Theil-Sen estimator/slope\n",
    "        intercept: intercept of Kendall-Theil Robust Line\n",
    "    Examples\n",
    "    --------\n",
    "      >>> import numpy as np\n",
    "\t  >>> import pymannkendall as mk\n",
    "      >>> x = np.random.rand(120)\n",
    "      >>> slope,intercept = mk.sens_slope(x)\n",
    "    \"\"\"\n",
    "    # res = namedtuple('Sens_Slope_Test', ['slope','intercept'])\n",
    "    \n",
    "#     x, n = __missing_values_analysis(x, method = 'skip')\n",
    "    n = len(x)\n",
    "    slope = np.nanmedian(sens_estimator(x))\n",
    "    intercept = np.nanmedian(x) - np.median(np.arange(n)[~np.isnan(x.flatten())]) * slope  # or median(x) - (n-1)/2 *slope\n",
    "    \n",
    "    return slope, intercept\n",
    "def sens_estimator(x):\n",
    "    idx = 0\n",
    "    n = len(x)\n",
    "    print(n)\n",
    "    d = np.ones(int(n*(n-1)/2))\n",
    "\n",
    "    for i in range(n-1):\n",
    "        j = np.arange(i+1,n) # all the combination of two points selected from original series. \n",
    "        d[idx : idx + len(j)] = (x[j] - x[i]) / (j - i)\n",
    "        idx = idx + len(j)\n",
    "        # print(i, idx)\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2778439-7fbd-4c2f-b3a7-49b913f95b0e",
   "metadata": {},
   "source": [
    "# Sub1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b2d142c4-823d-4fa0-a7fd-ff910659396c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40\n",
      "4593.61250524041\n",
      "40\n",
      "0.5369017953048782 False no trend\n",
      "0.6000506909101648 93.53823762690985\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f2f4bb665b0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABHtklEQVR4nO29eXRb93Xv+92YZ44gxUmiRFHzYDuqa8eJWzt249hNnQ5p7TapkyZNB2e9pL139Ta97W1ze93V2/ua3LRvJa9JM7hpUl+/Jk0Uj7FlZ3A8RZYlkRopa+IIzgRAEvPv/XHOAQ9BDOcAB8S0P2tpkQQB6qcjcGPju/f+bhJCgGEYhqkvTJU+AMMwDGM8HNwZhmHqEA7uDMMwdQgHd4ZhmDqEgzvDMEwdYqn0AQCgvb1d9Pf3V/oYDMMwNcUbb7wxK4TwZ/teVQT3/v5+HD9+vNLHYBiGqSmI6Fqu77EswzAMU4dwcGcYhqlDOLgzDMPUIRzcGYZh6hAO7gzDMHUIB3eGYZg6hIM7wzBMHcLBnWGYhuaJ0xOYX45V+hiGw8GdYZiGZXElho9/8008fny00kcxHA7uDMM0LEurcQBAIBip8EmMh4M7wzANSyiSAABMh6IVPonxcHBnGKZhUYL7DAd3hmGY+iEcbeDgTkQOInqdiE4R0Rki+rR8eysRPUdEI/LHFtVjPkVEl4joAhG9u5z/AIZhmGIJRSTNvSGDO4AogDuFEIcB3ADgHiK6BcCfAjgmhBgEcEz+GkS0D8ADAPYDuAfA54nIXIazMwzDlIQiy4SjCazEEhU+jbEUDO5CIix/aZX/CAD3A3hUvv1RAO+TP78fwGNCiKgQ4gqASwBuNvLQDMMwRqDIMgAwHayv7F2T5k5EZiI6CWAawHNCiNcAdAohJgFA/tgh370HgLppdEy+LfNnfoyIjhPR8ZmZmRL+CQzDMMURlGUZAJgJN2BwF0IkhRA3AOgFcDMRHchzd8r2I7L8zC8KIY4IIY74/Vm3RDEMw5SVcKTBM3cFIcQigB9A0tIDRNQFAPLHafluYwD6VA/rBTBR6kEZhmGMJhRJwOuQto3OhOprkElLt4yfiJrlz50A7gJwHsBRAA/Jd3sIwHflz48CeICI7ES0HcAggNcNPjfDMEzJhCJxbG11wWyiuhtk0rIguwvAo3LHiwnA40KIJ4joFQCPE9FHAFwH8H4AEEKcIaLHAZwFkADwsBAiWZ7jMwzDFE84moDPYUW7x1Z37ZAFg7sQ4jSAG7PcPgfgXTke8wiAR0o+HcMwTBkJRRLY2upCh9fRkJk7wzBMXRKKJOBxWOBP2TG11GCaO8MwTL0SisThc1jR4bVz5s4wDFMPCCEQjkrdMl6HBfPLUSRTAmZTtm7u2oMzd4ZhGpKVWBIpAXjsFnR47UgJYK6OBpk4uDMM05AovjJehxV+rwNAffm6c3BnGKYhCUcl6wGPwwK/1w6gvtwhObgzDNOQBNOZuyTLAMB0HU2pcnBnGKYhUXxlvHbO3BmGYeoGtebusJrhc1hYc2cYhql1FM1dMQ7ze+2cuTMMw9Q6SubukYN7vVkQcHBnGKYhUQqqHhtn7gzDMHVDOJKAx26BSZ5IlSwIIhBiw26hmoSDO8MwDUkoEk/r7YCUuUfiqXV7VWsZDu4MwzQk4aiUuSt0+JRe9/qQZji4MwzTkKhX7AFSQRWon153Du4MwzQkoWgCXoc1/bXfy5k7wzBMzROKxNNtkADSFgScuTMMw9QwoUgCPlVwb3JaYTOb6sZfhoM7wzANidIKqUBEddXrzsGdYZiGI55MYTWeXKe5A0A7B3eGYZjaZTm6ZverpsNrx3SQgzvDMExNkvaVsa8P7n6vHTN1smqPgzvDMA2H2u5XTYfXjvnlGGKJVCWOZSgc3BmGaThCkfV2vwpKr/vccu1n7wWDOxH1EdGLRHSOiM4Q0Sfk2/+KiMaJ6KT8517VYz5FRJeI6AIRvbuc/wCGYRi9hCK5NHd5UXYd6O5aMvcEgP8khNgL4BYADxPRPvl7nxVC3CD/eQoA5O89AGA/gHsAfJ6IzGU4e1UzG47iwS++imtzy5U+CsMwGSjmYJmaez0NMhUM7kKISSHECfnzEIBzAHryPOR+AI8JIaJCiCsALgG42YjD1hJHT07glctzODm6WOmjMHXC3zx1Ds8MT1X6GHXBmiyzXnOvJwsCXZo7EfUDuBHAa/JNHyei00T0FSJqkW/rATCqetgYsrwYENHHiOg4ER2fmZnRf/Iq5+nhSQBAcDVe4ZMw9cI3X7uOZ+TnFVMaoRytkO2eBsrcFYjIA+BbAD4phAgC+AKAAQA3AJgE8PfKXbM8fIP7vRDii0KII0KII36/X++5q5pAMILj1xYAAEsc3BkDSKYEwtEE5pZjlT5KXRCKJGA1E+yW9SHQZjGhxWWtCwsCTcGdiKyQAvs3hBDfBgAhREAIkRRCpAB8CWvSyxiAPtXDewFMGHfk6ufpoUkoy1yUVV4MUwqKRjzPwd0QpEUdVhBtzEU7vI7GyNxJ+td/GcA5IcRnVLd3qe72ywCG5c+PAniAiOxEtB3AIIDXjTty9fPU8BR2dXrQ4bWzLMMYgqIRz4U5uBtBpq+MGr/X3jCa+20APgjgzoy2x78joiEiOg3gDgB/BABCiDMAHgdwFsAzAB4WQiTLc/zqYzoYwU+vzuPeg13wOa0IRji4M6WjtO7NL8fqZsdnJclc1KGmo078ZbL/61QIIV5Cdh39qTyPeQTAIyWcq2Z55swUhADuO9iFH12cQXCVZRmmdJTgHktKOz4zuzwYfYSi+TP3mVAUQoissk2twBOqBvPU0CR2dngw2OmFz2nlgipjCCHVO0DW3UtHytyzv0D6vXbEkqmaT8w4uBvITCiK169Ikgwgmf+zLMMYQUhVmOeOmdIJR+PrFnWoWet1r+2OGQ7uBvLMmSmkZEkGAHwOKxdUGUNYl7lzUbVkQpHEuhV7atIWBDWuu3NwN5Cnhyaxw+/Grk4PAMDntCAYSXABjCkZZegGYFmmVIQQeQuq/jqxIODgbhCz4ShevTyH+w52pYswPocVyZTASqxhmoWYMhGKJKDU9liWKY1IPIVkSsBjz665d/hYlmFUPCtLMoreDgA+p/Tk4aIqUyqhSBzNTiucVjPm68COtpLksvtV8NotcFhNnLkzEk8NTWJ7uxt7tnjTt/nkajwXVZlSUbo7Wt02ztxLJJevjIKyKJs1dwZz4ShevTyPew9uWdcX2yRn7rXeUsVUHkUjbvPYWHMvkVxe7mrqwYKAg7sBfP9sAMmUwHsOdK273eeUnjzcMcOUiuSFYkGrm4N7qYRzrNhT4/dw5s5AkmS2tbmwv9u37naWZRijWCfLcCtkSSiae64JVUAqqnLm3uAsLMfw8ltzuFfVJaPgS8syHNyZ0kjLMm5bXez3rCRaZBm/x46l1Tgi8drtdOPgXiLfPzuFZErg3gxJBlh78iyx5s6USDASh89hRavbjkg8hZUYP6eKJV1QzdEKCay1Q86Ga/eFlIN7iTw5NIW+VicO9Pg2fM9qNsFtM7Msw5RESl7UoWTuAFv/lkJalsmXudfBuj0O7iWwuBLDy5dms0oyCj4nWxAwpbEcS0AIpAuqAE+plkI4koDLZobZlNvxUbEgqGXdnYN7CXz/bACJHJKMgs/B5mFMaYRU3R2tHg7upZLPekCBM/cG56mhSfS2OHGotynnfXxOC/e5MyWhBHePXSXLcHAvGi1++G1uG4iAmWDtWhBwcC+SSDyJn1yaxbv3b8lr6O9zsKc7Uxrqcfk1WaZ2M0qjGVtYwcVASPP9g5F43jZIALCYTWhz2zHDBdXG4/xUCPGkwJFtLXnvx57uTKmsjctb4bFbYDObOHNX8bdPn8fD3zih+f5aZBlA3qUa5ODecAyNLwEADvTklmQALqgypaPIMj6HBUQkTalyt0yauXAMowsrmq21lc6jQnR4OXNvSM6ML6HZZUVvizPv/XwOC0LRBFIp9nRnimNNlpF0YvaXWU8oGkckntIsf4Yi8bw97gqcuTcoQ+NLONDdVHCBrs9phRBAmIdOmCLJnKhkZ8j1KNdnSmPxM6xRlunw2jEbjtZsYsbBvQiiiSQuBkIFJRlgzV9maYWlGaY4QpE4zCaCy2YGALYgyEAJ7pNLhYN7MiWwHEvmHWBS8HvtSKQEFlZq84WUg3sRXJwKI54UOKgluDvZPIwpjVAkAY/dkn6X2Oq2s+YuI63Mk363pjQEdy2OkArpQaYa1d05uBfBWjF1o+VAJmu2vyzLMMWR2d3R5rFhOZasaVMro4gmUognJdlES+Yeisr1iwKtkIBqkKlGdfeCwZ2I+ojoRSI6R0RniOgT8u2tRPQcEY3IH1tUj/kUEV0iogtE9O5y/gMqwdD4EnwOC7a2ugrel21/mVKRvNzXMk22IFhD/XsV0BLcNThCKnTU+KJsLZl7AsB/EkLsBXALgIeJaB+APwVwTAgxCOCY/DXk7z0AYD+AewB8nojM5Th8pTgzsYQDPYWLqYB6GxMHd6Y4ghmZOwf3NZRgDQCTGgqqYXlmQKvmDtSuBUHB4C6EmBRCnJA/DwE4B6AHwP0AHpXv9iiA98mf3w/gMSFEVAhxBcAlADcbfO6KEUukcH5SWzEVUBVUObgzRRKKJOBTyzJsQZBGSZpsFhOmllYL3j+zrTQfbrsFbpu5rjP3NETUD+BGAK8B6BRCTALSCwCADvluPQBGVQ8bk2+rCy4GQoglU5qDu9dhAZGUfTFMMeSWZWoz6BiJkrkP+D3aNHcdsgwAdPgcmA7Vpr+M5uBORB4A3wLwSSFEMN9ds9y2oVGUiD5GRMeJ6PjMzIzWY1ScMxNSMVVLpwwAmEwEj93CsgxTNBsKqm5JLmBP97VgvavTg1AkgeVo/iQqHdw1FFQBaSNTXWfuRGSFFNi/IYT4tnxzgIi65O93AZiWbx8D0Kd6eC+AicyfKYT4ohDiiBDiiN/vL/b8m87Q+BI8dgu2aSimKrDtL1MsQogN4/I+pwUWE7HmjjWZZVenF0DhQaaQjlZIAPDX8C5VLd0yBODLAM4JIT6j+tZRAA/Jnz8E4Luq2x8gIjsRbQcwCOB1445cWYbGg9jf7YMpj9F/JpK/DMsyjH5WYkkkUwIe1bg8EaHFzRYEgDpzl4N7AWkmHJUGwhxWbaKF32Ov34IqgNsAfBDAnUR0Uv5zL4C/BXA3EY0AuFv+GkKIMwAeB3AWwDMAHhZC1EVDbjyZwrnJoGZJRsHnYFmGKY5wNLtG3MYWBACkzJ0IGPC7ARTudVckLi2dboC0SzUcTdTkztqCwpMQ4iVk19EB4F05HvMIgEdKOFdVcmk6jFgihYN5lnNko8lpxfX5lTKdiqln1F7ualo5cwcgNSp4bBZ0N0sGfoECskxYnvbVit+z1uu+rU3746oBnlDVgTKZur9bZ+bOtr9MkQTTdr/rNeI2jx1zNToWbyRKJu6wmtHismKyQDukNDOgTW8HpG4ZAAjU4JQqB3cdDI8vwW0zY0e7W9fjpIJq7b2tYypPrtY9lmUk1G2iW5qcmjR3rW2QANK/65dnwsUfskJwcNfB8PgS9nc36SqmAlJ3QziaQCKZKtPJmHol19BNq9uGUCSBWKKxn1PqNtEtPrs2zV2HLNPT7ITDasLFAAf3uiWRTOHsZBD7NZiFZaK8pQ5x9s7oJFfmrgwy1aodrVGEVJn4liZnQc1d64o9BZOJsLPDg5Fp7TtaqwUO7hp5a2YZkXhKd6cMoPKX4V53Rie5CqppC4IGH2QKqTT0riYHZsMxRBO5m/PC0YQmXxk1uzq8GOHMvX4ZHtc3maom7enOve6MTkKRBIgAty175t7oHTPrZJkmqfiZy6JX8X7XU1AFgJ2dHkwFIzWXnHFw18jQ+BKcVjN2+D26H6uYPtXak4OpPMqijsw6T5tHMQ+rvS4Oo8gM1lvkzpZcurvi/a6nFRKQMndAaoWuJTi4a2R4fAn7un0w6yymAurMnYM7o49gJL6hDRKQtjEBjZ25R+JSsFYy9y45c89lQRBKt5XqC+6DnVJCNxKoLd2dg7sGkimBs0VMpirwqj2mWHIVAJudVpiosYO7Uo/wZcgyuax/16Z99ckyvS0uOKymmtPdObhr4MpsGCuxJPZ36++UAdYKquzpzuhFkh02BneTidDiauxe92CGCZjXYYXHbskpyygvBnplGbOJMOD34CLLMvWHMpmq13ZAwW0zw0RcUGX0E8ozUdnqtjX0ouxsnUSdPnvOQSa9Xu5qdnV6cYllmfpjaCwIh9WEnUUUUwHJxc/nZNtfRj+hPF4oje4vk9bQnWsvfl1NzoKau95WSADY2eHBxFIk/YJSC3Bw18DwxBL2dvlgMRd/uXwO9pdh9JPp5a6mzWPDbAN3y2TLxLc0OfJk7opGr09zB4DBDimxq6WOGQ7uBUilBM5OBHFAp1lYJj6nhf1lGF0U6svmzH2jNUNXkwPToWhWq4/0cmydmjuw5hdfS0VVDu4FuDK3jHA0UXSnjEKT08oFVUYXSl92zszdbcfiSrxhPYuyZe6dPgeSKYHZLLWIUmSZvlYX7BZTTdkQcHAvgDKZqnUhdi5YlmH0Esxo9ctEGWRaWGnM55WyqMOjmt7N1+seisThtJphLUJeVTpmRliWqR+Gx5dgs5jSgwzFwntUq59kSuDhb57A905tWPlbEQrt+2x0CwJlUYd6ejdfr3sxvjJqBjs9LMvUE0PjS9i7xVvUq70an9PCrZBVzjdfu4YnT0/ixQvThe+8CRRq3VOCe6NaEGQb8OpqkjYyZet1D+p0hMxkV6cX44urae2+2uHgnodUSuDMeLBkSQaQMvfVeLLh/berlblwFP/r2Qvy59WRCefycldoa3ALgmzF5haXFTaLKWvHTFinl3smO2usY4aDex6uz68gFE0YEtybXGxBUM38r2cvYCWWxIDfXTWZsNbMvXGD+8ZMnIiwxefIqbnrtR5Qo7RD1orHDAf3PAyVYPObidJby0XV6uPk6CL+z/FRfPi2fty0taUKM/fswb1FThiq5bybTSjHyrwtTY6sskw4qm85diZbW12wWUw1U1Tl4J6H4YklWM2U7nEtBZ9Tsf2tDb2uUUilBP7yu8No99jxf71rUF48HYMQotJHK1hQtZhNaHZZGzZzD65mt2boyjHIpHcLUyYWswk72t2cudcDw+NL2L3FC5ul9MvEmXt18v+9MYpTY0v4s3v3wOuwot1jQyyZqooXYeUM+bLNRh5kymWqpsgymS/Q+Xx6tLKr08uZe60jhMDpsSVDJBmgtmx/j1+dxy1/cwxz4erQnsvF0koc//OZC/iZ/ha874YeAKolGFXwbw9F4vDYLXl3CLS5bVVTI9hMpOnd7MF6S5MDsURqXf9/KiVKboUEJN19bGEVyzXQMcPBPQfX5lYQiiRwqLfZkJ+nZO61MKX65vVFTAUjOHF9sdJHKSufee4CFldi+PQvHQCRFEDbPVIHSrYJx81Gi4zQ6rY1pOYeiaeQSGWf3lUGmSZVve7hWHGLOjIZlCXat2aqP3svGNyJ6CtENE1Ew6rb/oqIxonopPznXtX3PkVEl4joAhG9u1wHLzenxhYBAIeKtPnNpKmG9qgqxSiloFyPnJlYwtdfvYYP3rIN+1Q+/Up7YTVk7mFNwd3ekLJM2gTMmS1zl3rd1bp7WIPEpQVlmPFiDQwzacncvwbgniy3f1YIcYP85ykAIKJ9AB4AsF9+zOeJyGzUYTeTobEl2C0mQ4qpAOCwmmA1U03IMgG5jexMmYL782cDGJ1fKcvP1oIQAn/53TNodtnwx3fvXve9dlmWma2CgBmKxgsGoza3DQsrMaRSlS8AbybBPCvzsu1SLVSc1sq2Vhds5trwmCkY3IUQPwIwr/Hn3Q/gMSFEVAhxBcAlADeXcL6KcVremVrqZKoCEdWMv4zydrYcmftMKIqPff04/ulHbxn+s7XynZPjOH5tAf/lnt3p+QOF9NRnGTL3K7PLePz4qOb7aykAtrptSAlgsQaeV0aSr03U77XDbKJ0kgIA4Wj+tlKtWMwm7PC7a8KGoJTI9XEiOi3LNi3ybT0A1M/eMfm2DRDRx4joOBEdn5mZKeEYxpNMCQyPL+GQQcVUBWlhR/XLMoFgFCYCpkNRTOdYfFAsz5yZQkoAo/PZ91yWm1Akjr956jwO9zbh/W/r2/B9i9mEFpcVs2UI7l/7yRX8yb+fRiSe1HR/LZq7UgCeb7Ciar5M3GwidHjt6zL3YAmOkJns7PDUR+aegy8AGABwA4BJAH8v356trJ/1/aIQ4otCiCNCiCN+v7/IY5SHyzPSztSDBhVTFXwOS9UXVFMpgUAwgrdtk16vz0wEDf35T56WTLlGFyojy/w/L17CbDiK/37/gXWGU2qUXnejuTIn/Ztz7fjMRMtE5VqNoPIy0mZSaHo3c2lHKI+Mo5ddnV6MLaxiJVY4UfvGa9fw10+cLfnvLIaigrsQIiCESAohUgC+hDXpZQyAOh3qBVAdFns6ODUmyRGHDSqmKvic1S/LzC5HkUgJ3LGnA0TGSjPToQheuzIPu8WE8YXVigwKPTs8hTt3d+BwX3PO+7SVqQPl6uwyAGBiUdu7lmAkUTAYNaoFQSHfnS0+x/pumXRBtTTNHZDaIYUA3ppeznu/lVgCf/fMBXz5pSu4MLX5mX5RwZ2IulRf/jIApZPmKIAHiMhORNsBDAJ4vbQjbj5DY4tw2czYUeTO1FzUwh5VJdvZ6fdge7s77WdvBM8OT0EI4P1HehFNpDAT2lwpYX45hqtzKzjS35r3fu1eu+Hr62KJFMbkdyvjGoJ7NCGZzGmVZeYaLrgXztwDwajq/sZo7sBaO2QhaebbJ8axtBqH2UT42stXSv579aKlFfLfALwCYDcRjRHRRwD8HRENEdFpAHcA+CMAEEKcAfA4gLMAngHwsBBCm8BYRZweX8KB7qa8wyPFIBVUq1tzV4L7liYHDnQ3GRrcnzg9icEOD+7c0wEAGF3YXN391OgiAOCGPFk7ALS7bZg1+IVndGEFSkOLlsxda3dHi6txM/fMRR1qupocCEcT6aAejiZgIsBlK715b1ubC1Yz5W2HTKUEvvqTKzjY04RfP9KLb58Yx8Im/x9p6ZZ5UAjRJYSwCiF6hRBfFkJ8UAhxUAhxSAjxS0KISdX9HxFCDAghdgshni7v8Y0nnkzh7ETQsP52NdIe1SrP3INrwf1gTxMmliKGdI5MByN4/eo87jvUhd4WFwCkM9nN4s3RRZio8OxCm8eOYCRhqD2zIskAeoN7/kzTZjHB67A0XHDPtqhDTWaveygimYYpw2qlYDWbsL3djUt5MvcfjszgrZllfOQd2/Hh27Yjmkjhm69fL/nv1gNPqGZwMRBCNJHCwXIEd4cVsURKc7dEJZhaisBiIrS77djfIw33DBtQVH1almTuO9iF3hbpF29skzP3N68vYFenF+5CveMe47Phq3IxdVubCxOLhQuqhTRlNZIFQWMF90KdRJm97kb4yqgZ7PTmzdy/8tIVdHjtuPdgF3Z1evGOne34+ivXEN/Efbcc3DMYkoupRtkOqFmbUq3e7H1qKYJOnwMmE2F/t/QCZ4Q08+TQJHZ1ejDY6YXLZkGb27apmXsqJXBqdBE3bm0ueN81CwLjpJmrs8vwOiw40N1kaOYOKOZhjdUKGSzQSZS5SzWXyVixDHZ4MLqwgtXYxkTtYiCEH4/M4rdv3ZY2Hfzwbf2YCkbwzPCUYWcoBAf3DE6NLcHrsKC/zWX4z64F87CpYCS9h7LJacW2NlfJwT0QjOCnV+dx38Hu9G29ra5N7XW/PLuMYCSBG/taCt43PaVqZHCfW8b2djd6WpwYXyzcKaSnANjqLk/rZjVTKFh3+KQX6Kl1mbtxwX1Xp1fqmMniMfPVn1yF3WLCgzdvTd92x+4O9Le58NWfbF5hlYN7BkPjizjU22SINpeJ0ta2VMVF1amlSPotLQAc6GnC8ERpwf3poUlJkjm0JX1bb4tzUzP3k0oxVUPmXo7e8Suzy+hvc6O7yYFoIlVQ8lkbr2dZJhuFgrXdYka7x5aWZUpd1JFJeitThu6+sBzDt0+M4Zdv7EGb/A4QAEwmwofe3o8T1xfTz8Vyw8FdRSSexIWpUFkkGaD6M3chxLrMHQAOdDdhdH4ViyvFB48nhyaxZ4sXOzvWfHr6WlwYX1xFcpM8UU6OLsBjt2BAQ3vrWnuhMZl7NJHExOIq+tvd6G6W6g2FdHddsozHhoXl6lgwsllo0dA7fQ5Myb3upa7Yy6S/3Q2LiTbYEHzz9euIJlL48G3bNzzm1470wWu3bFr2zsFdxfmpEOJJYbjtgEK1L+wIRhJYiSUzMnepqFrspOrUUgTHry3gvoNd627vbXEinhSYDhlrb5CLN68v4nCftvZWj90Cu8VkWOY+Oi+1QW5vd6WDe6Fedz0uhm1uGxIpUfVttkYSisTT281y0aVat2eEl7sapWNGXVSNJ1P4l1eu4h0727F7y0bDQY/dgvcf6cOTpyfX+d6UCw7uKoZkm99ydMoA1V9QDajaIBUOyEXVYidVnx6WJJl7D20M7sDmdMysxpI4PxUq2N+uQERo99gxY5DmfmVW6ZRxoyeduef/d4cicbhsZlg0GNelzc4apKiab1GHGmmQSXpOBw3W3AFJd1e3Qz41NIlAMIqPvGNj1q7wobf3IykE/vXVa4aeJRsc3FWcGltCm9uW/gU0GuXJVa3mYZNLG4N7i3w9ii2qPnlakmQy5ZC+VqlgvRnWv8MTS0imhKZiqkKbxzgLAqXHfXubG80uK5xWs4bgrl0jbjQLgnyLOtR0NTmxsBJHMBJHLJHSVL/Qw84OD67NryAST0IIga+8dAU72t34uV25vbK2trlw195OfOO162VviebgrmJobAkHy1RMBQCH1Qy7xVS1mbuiT6plGQA42FPcpOrk0iqOX1vAL2Zk7QDSL6Cbkbm/eX0BgLZiqkK7x25YJnx1bhlNTita3DYQEbqbHZhYKhDco9pb95TWzUYpqmqdAVCex4oubmRBFZAWdygdMyeuL+LU2BI+fFt/zsEqhQ/f1o/55RiOniyv7RYHd5mVWAIj06Gy6e0K1ewvM7UkBbPOjOB+oMeHq3Mrus/99JDU03vvwY3B3WE1o8Nr35TM/eToInpbnOkgqIU2tw2zIYMy97ll9Le70193NzsxrqGgqrUA2GiZe75FHWqUd6CKdFIOWUb6+WF85aUr8Dks+JWbegs+7tYdbdizxYuv/ORKWYvgHNxlzkwEkRLlGV5S43NYqrbwNRVcRbvHlh68UDggv+CdGddXVH1yaBJ7u3w5DdikdsjNyNwXceNW7ZIMINv+LkcN+eW7OruC7aq5iZ5mZ0FZRo9G3GjBXesMgBLcy5W597dJHTM/uDCDp4cn8eDNWwtOPwNSTefDt/Xj/FQIr1yeM/RMaji4y5xOT6aWN3Nvclqr1tNdmU7NJB3cdfS7Tyyu4o0ckoxCX6ur7L7ugWAEk0sRzcVUhXaPDfGkKLk+EoknMbG0uiFznwlFEU3k1lxDkbhmjdhhNcNtMzfMIJNWUzVFlrk4HdZ0f73YLCb0t7vxH2+Og4jw22/v1/zY+2/oQavbhq/+5KqhZ1LDwV1maGwRW3wOdGQJbkZSzbLM5FIkPbatpt1jR1eTQ1fHzFNDkpdcNklGobfFicmlCBJl9Nt48/oigMJOkJmkdewSO2auz69ACGB7RnAH1i9wzkTvRGWrp3otCJIpgc89P2KYBKd1BsBtt8DnsOBSoDyyDLA2zHTPgS26GjEcVjN+8+ateP5cANfnypPgcHCXOS0XU8tNNe9RDQSzZ+4AsF+n/e9TQ5PY3+1bF9Qy6WtxIZkSaf+PcvDm6AKsZsL+bp+ux7WlLQhKy4avyJ0y/W3q4C5d43y97nq9UFrd9qotqJ64voDPPn8R//HmuCE/T4+p2pYmBybkF9FyBHdFd/+dLENLhfjgrdtgJsKjr1w1+FQSNR3cV2IJvHFtoWSZIxiJ4/LsctmLqYBi+1t9mnsknsTCSjxr5g5IHTOXZ5cRjhY++/jiKk5cX8ybtQNIW/+W02Pm5PVF7OvywWHV5+O9ZkFQWjastEGqZZmeAlOq8WQKkXhKl4xQru1RRvD82QCA7D4sxaBnelex/pXub6wsAwAfuGUbPvfADbhJRyeWQqfPgY+8Y3veBKgUajq4n5sM4Ve/8DKOX50v6ecMK3q7zrfuxaBk7tU2Kq4Me+TK3A/0+CAEcG6ycFH1aVmSyZxKzWRtkKk8b0sTyRSGxpd0SzKAyjysxGz46twyWt229AAbsFboy1VU1RO8FCRnyOoM7s+dMza4Bwss6lDTpXo+G11QBQC/1477b+gpun36U/fuxQdu2WbwqSRqOrhvkzsQrpeo5Z2W5YaDm5C5NzmtSKQEVrJYhVYSZYCpqym7bqgUVRVL5FzE5KUEB3p867LVbHQ3O0FUvo1MFwPSonO9nTKAauqz5Mx9ZYPDqN1iht9rzxPctcsOCm1ycK+2pOGtmTAuzyzD57Dg8syyIecLFVjUoUZ5IbVbTBu6wOqdmv7XtrltcNnMuFZiQWJobAm9Lc70L3Q5qVbzsDXrgey94J0+B/xee0GHyC+/dAWXZ5bxx3fvKvh32iwmbPE5ypa5n9S4Vi8bFrMJLS5ryba/V+eW1+ntClKvu7GZeyyZ0iSbbSbH5Kz9A7dsw0osaUh9JaijHqEE93Lo7dVOTQd3IsLWVlfJVfhTY4s4XOb+doU187Dq+iVcsx7IXfE/0O3L2+s+trCCfzg2gnfv78Sdezo1/b19LS6MlUlzf/P6Alpc1vQ7PL20eUrzSV+NJTG5FMn6Dqan2ZEzcw8Wscy5Wnvdnz87jX1dPrxzUBrJf2t6ucAjChOKJNJJUiHWgrvxenu1U9PBHZCkmWslBPf55RjGFlY3pVMGQNrJrtoy96mlCDx2S15d8mBPE0amQ1m3zwDAp793FgDw3967X/PfW05f95Oji7ihr7loPbTUIuW1+Y3FVIXuJicmFiNZZYqwDi93hTWb4uoJ7vPLMRy/No+79nVioEO6Bkbo7no6iZQGgXLo7dVOzQd3JXNPFekLrvRul3t4SaFabX+nltb7uGdjf08TUgI4N7Uxe3/+bADPnQ3gk3cN6ur37W11YSoYMXQZNSC9eF6aCeMGHWZhmbR77ZgtoXdcbRiWSXezE6vxJBZXNj4PipNlJDltvoo6Zl48P42UAO7e2wm/xw6vw2JQcNduzdDlk56LLMvUIFvb3IgmUpgOFfdLeFrWZQ9sQjEVWLP9rbYp1algZINhWCZKwTmz3301lsRfHj2DXZ0e/E4eu9Ns9LY4kRKSyZiRnB5dghDQtDM1F+1uG2aLfF4Ba1a//e0bZaF8vu5KQVVPttnmNn41YKk8fy6ATp8dB3p8ICLs8HsMDO7aro3PaYHDauLgXotsbS2tY+b0+BJ2tLsNtwPNha9KPd21ZO5dTQ60um0bgvs/vjCC8cVV/I/3HYRVg/+4mj65191oj5mTo5IT5OES2lvbPHYEI4mi31VcnV1Gu8eWNcvM5+uudbxeTVeTA1Yz4WqZph31Eokn8cOLM7hrb2daFhvwuw3S3LXLMkSEm7a2YPcWfUNs9UDNB/dtcnC/Nlfck+b02OKmSTJAdXq6J5IpzISjBTN3IsKBniYMqYqql6ZD+NKPL+NXb+rFzdtbdf/dSq+70e6Qb15fxA6/e11/uV4UC4Jii5RXcnTKAGtTqlmDezShu3XPYjahv82NS9PG9JKXyquX57ASS+KufWuF9QG/B1PBSEkdPVoXdaj55u/eoql7q96o+eDe3eyEiYrL3APBCALBKA5uUqcMIK3nctnMVZW5z4ZjSKZEwcwdkDpmRgKh9IKCP//OMFw2C/7s3j1F/d1dTQ6YTWRo5i6EwMnRRV3LObKxZkFQnNRxLcPqV02r2wa7xZRTlimmu2NnhzGyhxE8fy4Al82MW3e0pW9TFrZcLuGMWhd1MBqCOxF9hYimiWhYdVsrET1HRCPyxxbV9z5FRJeI6AIRvbtcB1ewWUzobnYWFdyHNskJMhOfo7rMw5Te40KZOyDp7omUwIWpEL57cgKvXp7Hn9yze92mdz1YzCZ0NRnb6z62sIq55Ziu5RzZaC8huK/EEggEozlHy4lItv7d2PcdjCQKepVnY2eHB9fmlvO6TeZieHwJ//OZ84YMGQkh8PzZadw+6F9n+7DTgI6ZYga8GhUtmfvXANyTcdufAjgmhBgEcEz+GkS0D8ADAPbLj/k8Eekz9SiCra2uogaZTo8twkTQbSpVKk1Oa1X1uac3MGnJ3OWi6stvzeF/PHkOh/ua8eDPbC3p7+9tcRo6pXpC3rx0Y4l2Emv+MvplmatKMTWHLAPkHmTS6wipsLPDg5RY+7v18NhPr+MLP3hr3cLnYjkzEcRUMLJOkgGAra1umE1Uku6udVEHoyG4CyF+BCDTvOV+AI/Knz8K4H2q2x8TQkSFEFcAXAJwszFHzc22tuIGmc5OhrC93Q2XBo8KI/E5LVXVLTOVth4oHNx7W5xoclrx2ecvYn45ikfed0DTGHg++lpchmbuJ0cX4bCasCfLBno9tHuV9XX6M/ercg0o3wBVd45BpmJlGUX2KEZ3vzAl2eIeOx/Q/dhMnjsbgImAO3av3yVqs5iwrdVlUObOwb0QxWrunUKISQCQP3bIt/cAGFXdb0y+bQNE9DEiOk5Ex2dmZoo8hkRfqwtzyzHdhZqR6RB2lxgAiqHaZJnJYAQ2s0mT/YJUVPUhlkjht2/tN6SFtLfFhUAwatjC4JOjizjY0wSLzs6dTNw2aedtMba/V7K4QWbS3ezEdJalHcVm7gN+D4j0B3chBM7Lwf2Fc9O6/95Mnj8XwNu2tWSV6kpthwwW0UnUqBhdUM2WwmUV8YQQXxRCHBFCHPH7c28L18K2VukXSI/p/WosievzKxjsqEBwr7KFHYGlCDqb7JonOW8f9KOv1Yk//gVjOhD6WnO3BeolmkjizHiwKLOwTIgI7R57UZr71dll+L32vL3qSq97YGn9z9fr5a7gtJnR0+zEJZ3Bc3IpglAkgU6fHSeuL5RkYTCxuIozE0HctTe7/cRAhxtXZ1eQLHLokDN37RQb3ANE1AUA8kfl5X4MQJ/qfr0AyrviG+ped+1a3lszYQixZra/mVTbHtXJpcIDTGp+7+cG8MP/fIdhswFpX3cDdPdzkyHEkqmizMKy0eYpzoLg6txy1slUNT05Bpn0tvqp2dnh0Z25K5LM775zB1IC+OHF4rN3xSgsU29XGPB7EEumipbhipkBaFSKDe5HATwkf/4QgO+qbn+AiOxEtB3AIIDXSztiYbYWYf17UV69tasz+/LmctLktCIUiRdtmQAAn3jsTfzFd4YL31EDgWAkr2FYNkrV2dUombsRuvuJa1Ix1ajg3i4vytbLldmVrJOparqzDDIlkimsxJJFZ6Y7/R5cngnryowVSebX3taLdo8dL5wvXiZ97tw0drS70/p/JgP+0jpmOHPXjpZWyH8D8AqA3UQ0RkQfAfC3AO4mohEAd8tfQwhxBsDjAM4CeAbAw0KIshuXNzmtaHZZdXXMXAyEYTVTQc/xcuBzWpESQDhWXPY+HYrge6cm8MyZqZJb14QQcuZeXCujEXR4penKUnvdUymBx4+PYmeHR1NxWAttbhtmQ/oy93A0gdlwtOBzqyvL0g6lblRK5h5NpDCu41pemAqiq8mBZpcNd+z244cXphEvYq9tKBLHK2/N5szaAWBHuxT0i+2YCUUSmhd1NDpaumUeFEJ0CSGsQoheIcSXhRBzQoh3CSEG5Y/zqvs/IoQYEELsFkI8Xd7jr7G11aUrcx8JhLCj3aN7XN4ISjUPe+r0JFICmAlFS/bHXlyJI5pI6c7cjcRsInQ3O0ueUn3uXADnp0J4+I6Bop0gM2mTM3c9L6L5DMPUOKxmtHtsmFD56hRjGqZmp7yw+dJMSPNjzk+tNRa8a28HghFpfaVefjwyi3hS5NTbAaDFbUOb21ZC5p6Ax65tUUejU/MTqgp6g/vF6RAGKyDJACrb3yJ196OnJtKFulOj2pdWZ0PPAFM5kdohi8/chRD4h2Mj6G9z4b2Hug07V7vHhnhS6LKLUNogtbwrlHrd116gQyX2caeDu0bdPZ5M4a2ZcDq4v2PQD6uZ8MJ5/br782cDaHFZC+4THSihYyYYiW+aD1StU1fBfXxhFQkNbyeXowmMzq9WpJgKqDL3IjpmRudXcOL6Ij7yju2wmAinxxZLOstUeklHZYN7qb7uL16YxpmJIB6+Y2fJLZBqFH8ZPev20kuxC2TugOLrrs7cS5vAbHbZ0O6xaQ7uV2aXEU+K9EyAx27BLTva0oVRrSSSKbxwYRp37OkoeP0HOtx4a6Z4WYb1dm3UTXDf1uZCIiXSG4XyoTzxK1FMBUpzhnzitLR8+tfe1os9XV6cLrDTtBDpzL3Cwb2v1YXZcCznIpB8CCHwuWOX0NvixPtuzDpWUTRr/jLadfcrsyvY4nPAaSs8nN3dLAV3RfZRMvdSlksM+LV3zCjF1N2da1Pad+7pwFszy7rM+N64toDFlTjuziPJqM83vxwrquWy2DbRRqRugnufDutfpVNmsEKZeyme7kdPTeDGrc3oa3XhUG8zTo8tltR1M7kUARHQ4a1cQRVYc4csJnv/0cgsTo0u4uE7dhpeQ1mzINCRuc8ta17t193swEosmX4uhKKld4MMyO2QWuoEF6dCMJsovSkJkII7AF3SzPPnArCZTXjnrsIzK6UYiJXSJtpo1E1w3ya/BdbSMTMyHYbNbErbBW82a7KMPs390nQI5yaD+KXDkqZ8uLcJwUgirfEWQ2ApgnaPvSKFZTW9Rfq6K1p7d5MDv3pTr+HnavfKmbuOLPPq7HJOw7BMMnvdjejj3un3IBhJYEbDC9L5qRB2tLtht6y9y9jW5saA3605uMcSKTw1NIVbBto0veNQgnsxujvLMtqpm+C+xSe102nN3Hf43YZqs3rwKJ7uOjP3o6cmYSLgvoNdAIBDslVxKdLMZDBiWNtgKfQVmbm/8tYc3ri2gD/4+QFd/udaaXXJu0k1Zu7BSBxzyzHNLbZrve6SPFZqtwygr6h6IRDMasHxrr2dePXynCZLj3955SrGF1fx4bf3azpfT4sTNoupKN2dZRnt1E1wN5sIfS0uTVOqI4FwxYqpgHRWr92iq6AqhMD3Tk3glh1t6JA7WwY7PHBYTThVQlE1sBRBZ4U7ZQDA77XDbjHpnlL93LERdPrseP+RvsJ3LgKL2YQWl1WzBcE1DW6QajIHmYKROGxm0zqrXL0owf2tAsE9LDcWZDNYu3NPB+JJgZdG8g80zS/H8LljI/i5XX7csacj730VzCbC9jZ3wfNlUsyijkamboI7IOnuhTL3cDSB8cXVihiGqfHptP09MxHEldllvPfwWpufxWzCge6m0jL3pdWqyNyJCD06O2ZeuzyH167M4/duHygpGBaizWPXbEFwRZbItMoybW4bbBZTOrgbITt0NTngtpkLZu5rU9obfxfetq0FXocFxwoYiX32uYtYiSXx5/ft1XXGgQ43Ls/qy9x5UYc+6iq4b2uTfN3zFZJGlGJqR2U6ZRR8TquugurRUxOwmAjvObBl3e2HeptxZmJJUwtoJiuxBIKRRFVk7oCku4/Oa8/c//GFS2j32PHgzaX5yReiXYe/jNIGqbWgajIRupsc6zT3UoMXEUlF1QKatuIpsyfLflGr2YSf2+XHixdmchbsL0yF8I3XruEDP7tVd3PCgN+D6/MruhaLBHlRhy7qKrhvbXUhFElgcSV30BwJKG2QFc7cHdplmVRKkmRu3+VHs2u9Le/hviZE4qmilizo8XHfDPp0ZO5vXFvAS5dm8Xu379DUclgKbTqcIa/OLqO7yaHrnYTSDgkU7+WeyU4N7ZAXpkJw2czpTqVM3rW3A7PhKIbGN74zFELgr584C6/Dik/epd8ddMDvQTIldDm5KjMAvKhDG3UX3IH87ZAXAyHYLaZ062SlkGQZbcH9jesLmFyKpLtk1ByWi6rF6O7VMp2q0NviwsJKXFMR7x9fGEGr24bfuqW8WTsAtLttmoP7lTx7U3PRrVq3Z1Q3yECHB4FgNG8CcX4qiF2d3pyj/D+3qwMmAo5l6Zp54fw0Xro0i0/eNYgWDXsANpyviI6ZoAHF5kairoJ7uh0yT3C/EAhhZ4cH5gp7U/gc1nRnRCGOnpyA3WLKasi0rc2FJqe1qEnVaplOVdDqDnlqdBE/uDCDj75z+6Zs0Wrz2BGMJBBLFJa+rs4WF9wDoQjiyZRh3SBKUfVyjo4UIaQ9uPm2VbW6bbhpawteyNjOFEuk8MiT57DD78YHbtlW1Pl2pN0htevubPerj7oK7kpwyGdAVelOGQWf06Ipc08kU3hqaBJ37e3M2kNMRDjU21SUx0y1TKcqpHvdC+ju//jCJTS7rPjtW/s34VRrFgSFJiqXVuJYWIkXNAzLpKfZASGkF1ujukEKtUPOhKJYWIkXbCy4Y08HhseDCKgM6r7+6jVcnl3GX9y3r+j5CLfdgq4mh66OmTVZhoO7FuoquLtsFvi99pxj00urcUwFIxUzDFPT5LQiFE0U9N1++a05zC3H1nXJZHKotwkXAiHda+qmliLwOSybvkM2F0qv+2iezP3lS7N4/lwAv3Pb9pJG9PWwZkGQX5q5omFvajbU7ZBhg2SZba0uWM2UM7inbQcKBPd37ZXaG1+UpZn55Rg+9/xF3L7Lj5/fXXgaNR96DcSMmAFoJOoquAP53SEvTcutXxVYrZeJkn0srOTPBo+emoDXbsn7i3SotxnJlMCZiaCuM0wtRaomawckGcBpNeecUn3i9AQ+9LWfYoffjQ/d1r9p52rXGNyHZGlM7ztDJbiPLawiHDMmc7eYTehvc+cM7vk6ZdTs7vSip9mZ1t0/+9xFLMeS+Iv79pZsqzzglwzEtNop86IOfdRdcN/W6spZgVc6Sird4w4Ae7qkM/z6//sKXrs8l/U+0UQSzw5P4Rf2b8nbfaFsHdKru08VsYGpnBAR+lo3+roLIfClH13Gx7/5Jg71NOFbv//2TX1rvuYvk/+F+JkzUxjwu/Vr7vL/wcXpEIQwrhtkZ0fuzPj8VAh+r73gUnQiwp17OvDSyCyGxpaKbn3MxkCHB+FoAtMhbcVqZVGHu0reaVY7dRfc+1pdmAxGsvbPXgyE4LSa034eleTtA+34+kduRjyVwm988VX82X8Mbehs+MGFGYSiCbz3cFfen9Xpc6DTZ8ep0UVdZ5iq8AambPRm+LonUwKf/t5ZPPLUOdx3sAv/+tGfLao7oxTaZVO1fOv2FpZjePXyPO7JmEPQgtNmRqvbls6mjcpMd3Z4cG1uOevvwoVAMG8xVc2dezuwGk/io//yU3jslqJaH7OxtpVJmzTDizr0UXfBfVubC0JkN6AaCYQx2OmpmifHOwf9ePaTt+N337kdj71+HXd/5od49sxU+vtHT02g1W3DbTvbC/4sySFSe1E1nkxhJhytqswdkNwhFc19NZbEH/zrG/jay1fx0Xdsxz8+eGNZJ1Fz4baZYbeY8tr+PncugGRK4J79+V+Ic9Hd7EgHd4/dmHclOzs8SAng6uz6d0LJlMBIIIzdGrPvW3e0wWE1IRCM4pN37TLsxVVxotSqu/OiDn3UXXDP1+t+MRDCYBXo7WpcNgv+63378J2Hb0Or247f+/ob+MNvvIGrs8s4di6A9xzYoqkj4XBvEy7PLmueep0JRSFE9QwwKfS1SINoV2aX8Zv//CqeOxfAX753H/78F/dV7EWZiNBeYJDp2eEp9DQ7caAnv4adi+4mZ3oXgVGZu9JLnqm7X51bRjSR0ixPOqxm3LW3E7s7vfjgrcW1PmZji88Bl82suR2SHSH1UXdXaqvcqZCpuy+txDEdilZsQUchDvU24+jHb8OXfnwZ//v5EXz/TACJlMg6uJTr8QAwPL6kKdNXAkm1DDApKNOSv/L5n2AllsQXfusm3HOguGzYSNryWBCEown8eGQWH7x1W9FFxm6VVGhkcCfaGNwvaiymqvnMr9+AZEoYag1NRLo6ZtgRUh91l7n7PXY4reYNmfvF6dwmSdWC1WzCH/78TjzziXfi5u2t2Nvlw8/0t2p67KHeJgDaJ1UDVdbjrqCeHP7m7/5sVQR2QOp1z6W5v3B+GrFkqii9XaFnXXA3Rnpw2qT6UqbHzPmpEEwEXS3BNoupLDYPA353zkGrTNgRUh919zJIRNja6tqwtEPRM6uhx70QO/wefPN3b4EQQnMm2Oyyob/NhdMah5mqNXPf2+XDJ+8axP039Gh2VtwM2tw2nM3Ravrs8BT8XjvetrWl6J+vztyN9E7Z2bHRY+bCVAj9be6K1C8yGfB78J2TE1iJJQrOW4QiCezsqLuQVTbqLnMHpOwvs51uJBCC21YdnTJa0fsW/1Bvs67M3W4xodlVXZmQ2UT45F27qiqwA7Lt73J0Q092JJ7Eixem8Qv7OkuqCXQ3r73IGpmd7vR7cHkmvG5Y7kIgVBXtwIDUDgnktklQE2RZRhd1Gdy3tUmDTOpfxIuBMHZ2eksevKhmDvU2YXIpgulQ4SXhk/IAUz1fDyNp99gQT4oNqxF/dHEGK7FkSZIMsCbLWEwEh9W4X8udHR5EEymMy91jq7Ekrs4tV09w12ggxos69FPSs4iIrhLREBGdJKLj8m2tRPQcEY3IH4t/r1ok29pcWI0n1+2QHJkOYXcNSDKlcFgZZtIgzQSWIlUnyVQzir9MZsfMM2em0OS04pYdbSX/fKuZ4HVYDH3BVTLjSzOSLDkiD0pp7XEvN9vaXDBRYQOx1XgSSV7UoQsjUoQ7hBA3CCGOyF//KYBjQohBAMfkrzcVpSindMzML8cwG45VdTHVCPZ3+2AibZOqk8HVqiumVjOKv4y6YyaeTOH5swHctbez5C4Sk4nQ1eQ0PDPdmdEOqXjKVMvvgsNqRm+Lq2Dmzo6Q+imHLHM/gEflzx8F8L4y/B152ZbR666sEzNiZLqacdks2NXpxakCw0xCCASWohzcdbBmQbCWub96eQ7BSKJkSUaht8VpeA2kxW1Dm9uWDu4XpkJwWE1pe+xqYMBfeJ8qL+rQT6lXSgD4PhEJAP8khPgigE4hxCQACCEmiSjr1lwi+hiAjwHA1q3GLlzoaXGCCOmOmZH0rsj6lmUAaXnH989O5e20mV+OIZZMsSyjg3avbB6msv19engKLpsZ7xwsPFeghb/4xX2IavCM18uAqmPmwpQ0yFfpfQZqBvwevPzWHJIpkfNcSq2DJ1S1U2rmfpsQ4iYA7wHwMBHdrvWBQogvCiGOCCGO+P2lWYdmYreY0d3kVGXuYXjtloYIZof6mrCwEs+5i3Q1lsSf/PtpAPqGWBqdVnm94axscpVMCXz/TAB37OkwrKVwb5cvbQJnJEo7pBAC56eqp1NGYSCj6JsNtvvVT0nBXQgxIX+cBvAfAG4GECCiLgCQP+Zfn14m+lqd62SZwU5PQ3SG5Fu7N78cw2/+86t44cI0/vp9B3DrQGlFwEbCYjahxWVNDzKduL6A2XAU9+w3RpIpJzv9HgQjCVwIhDAbjlZNMVVBGcD74chMzvuEeDm2booO7kTkJiKv8jmAXwAwDOAogIfkuz0E4LulHrIYtrW6cW1Oaoe8GAhVTQGp3Oze4oXNYtpQVB2dX8GvfeFlnJkI4gu/9TZ8sMj1aI1Mm8eeLqg+PTQFm8WEO/ZkVR2rCmUr01OnJwFUh+W1mn1dPuzr8uHfXrue09udM3f9lJK5dwJ4iYhOAXgdwJNCiGcA/C2Au4loBMDd8tebztY2F2bDUYzOr2JhJV73xVQFq9mEfV2+dUXVMxNL+JUvvIzZcBTf+OjPGlYAbDTaZX8ZIQSePTOF2wfbN20bVCkowf2JoeoM7kSEB2/uw9nJIIbGszcD8KIO/RQd3IUQl4UQh+U/+4UQj8i3zwkh3iWEGJQ/zht3XO0o7pDH5OW+Wu1N64HDvU0YHl9CMiXwk0uz+I1/ehVWE+Fbf/B2zV41zEbaZGfI4fEgxhdX8e4akGQAyfnTbTPj8swyWt02+D3V5eEPAPff2AOH1YR/e3006/d5UYd+6nJCFVgL7s+fk4J7I3TKKBzua8ZKLInPPX8RH/rq6+hpduJbf/j2hnn3Ui7a3TbMhqN4engSZhPhrr2dlT6SJogoPcy0u0qntH0OK+472I2jJ8exHE1s+D4v6tBP3QZ3ZUnxa5fn0eS0wu+tvmylXCj2v//wwiXcuLUFj//+reiqsqUctUibx45gJIEnTk/i1h1tm74RqhSUYaZqk2TUPHhzH5ZjSTxxemLD93hRh37qNrg3Oa3wOixIpAR2NUinjMKOdjcGOzz4pcPd+JffuRlNTv6lMALFguD6/AreXWN1CyVzr7ZOGTVv29aCnR2erNIML+rQT90GdyJKZ++NJkeYTITv/9Ht+IcKraWrVxQLAiLg3ftqQ5JRONgjtRsq7+qqEamwuhUnRxdxfmq9vXJwlR0h9VK3wR2Q2iEBYFdH4+jtCo30TmWzaJeD+01bW9BRYwNx7xxsx4v/+eexr7u6B9d+5cYe2MwmPJaRvbMjpH7qOrgrBmKN0uPOlBdlmfh7akySAaQX+2rzyM9Gi9uGew5swbdPjCEST6ZvD0U5c9dLXQf3w71NcNvM2NtV3dkKUxv0NDvx779/Kx56e3+lj1LXPHBzH4KRBJ6S+/IB1tyLoa6D+z0HtuCnf35XTXU1MNXNkf5WQ5dEMxu5dUcb+ttcaWlGWdTB3TL6qOtnKREV3MvIMEx1QUT4jZ/ZitevzuPSdFi1qIODux7qOrgzDFOb/NrbemExEf7PT6+zr0yRcHBnGKbq8HvtuGtvJ751Yjxt1sbBXR8c3BmGqUoe/NmtmF+O4VsnxgDwog69cHBnGKYqeefOdvQ0O/HY69cBcOauFw7uDMNUJSYT4Td+RvKbAXhRh144uDMMU7W8/0gvFCNIztz1wcGdYZiqpavJiTt2S9uuOLjrg68WwzBVzZ/csweHeptrYutVNcFXi2GYqmb3Fm9V+9BXKyzLMAzD1CEc3BmGYeoQDu4MwzB1CAd3hmGYOoSDO8MwTB3CwZ1hGKYO4eDOMAxTh3BwZxiGqUNICFHpM4CIZgBcK+FHtAOYNeg4RsNnKw4+W3Hw2YqjVs+2TQjhz/aNqgjupUJEx4UQRyp9jmzw2YqDz1YcfLbiqMezsSzDMAxTh3BwZxiGqUPqJbh/sdIHyAOfrTj4bMXBZyuOujtbXWjuDMMwzHrqJXNnGIZhVHBwZxiGqUNqOrgT0T1EdIGILhHRn1b6PGqI6CoRDRHRSSI6XuGzfIWIpoloWHVbKxE9R0Qj8seWKjrbXxHRuHztThLRvRU6Wx8RvUhE54joDBF9Qr694tcuz9kqfu2IyEFErxPRKflsn5Zvr4brlutsFb9uqjOaiehNInpC/rqo61azmjsRmQFcBHA3gDEAPwXwoBDibEUPJkNEVwEcEUJUfDCCiG4HEAbwL0KIA/JtfwdgXgjxt/ILY4sQ4r9Uydn+CkBYCPF/b/Z5Ms7WBaBLCHGCiLwA3gDwPgAfQoWvXZ6z/ToqfO2IiAC4hRBhIrICeAnAJwD8Cip/3XKd7R5UwXMOAIjojwEcAeATQvxisb+rtZy53wzgkhDishAiBuAxAPdX+ExViRDiRwDmM26+H8Cj8uePQgoMm06Os1UFQohJIcQJ+fMQgHMAelAF1y7P2SqOkAjLX1rlPwLVcd1yna0qIKJeAPcB+GfVzUVdt1oO7j0ARlVfj6FKntwyAsD3iegNIvpYpQ+ThU4hxCQgBQoAHRU+TyYfJ6LTsmxTEclIDRH1A7gRwGuosmuXcTagCq6dLC2cBDAN4DkhRNVctxxnA6rgugH43wD+BEBKdVtR162Wgztlua1qXoEB3CaEuAnAewA8LMsPjDa+AGAAwA0AJgH8fSUPQ0QeAN8C8EkhRLCSZ8kky9mq4toJIZJCiBsA9AK4mYgOVOIc2chxtopfNyL6RQDTQog3jPh5tRzcxwD0qb7uBTBRobNsQAgxIX+cBvAfkGSkaiIg67aKfjtd4fOkEUIE5F/AFIAvoYLXTtZlvwXgG0KIb8s3V8W1y3a2arp28nkWAfwAkqZdFddNQX22KrlutwH4Jble9xiAO4noX1Hkdavl4P5TAINEtJ2IbAAeAHC0wmcCABCRWy5ygYjcAH4BwHD+R206RwE8JH/+EIDvVvAs61CeyDK/jApdO7n49mUA54QQn1F9q+LXLtfZquHaEZGfiJrlz50A7gJwHtVx3bKerRqumxDiU0KIXiFEP6R49oIQ4gMo9roJIWr2D4B7IXXMvAXgv1b6PKpz7QBwSv5zptJnA/BvkN5qxiG94/kIgDYAxwCMyB9bq+hsXwcwBOC0/MTuqtDZ3gFJ6jsN4KT8595quHZ5zlbxawfgEIA35TMMA/hv8u3VcN1yna3i1y3jnD8P4IlSrlvNtkIyDMMwuallWYZhGIbJAQd3hmGYOoSDO8MwTB3CwZ1hGKYO4eDOMAxTh3BwZxiGqUM4uDMMw9Qh/z/57lAd0C+AbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sub_id = 1\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    year_precip = year_precip[year_precip>1]\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "    \n",
    "n = len(n_epe)\n",
    "print(n)\n",
    "print(np.sum(n_epe))\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)\n",
    "plt.plot(n_epe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0e4ec5d2-9fc1-4f6f-a7bb-838ff2d8278d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.015\n",
      "Model:                            OLS   Adj. R-squared:                 -0.010\n",
      "Method:                 Least Squares   F-statistic:                    0.5974\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.444\n",
      "Time:                        15:46:54   Log-Likelihood:                -224.55\n",
      "No. Observations:                  40   AIC:                             453.1\n",
      "Df Residuals:                      38   BIC:                             456.5\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        100.7867     21.129      4.770      0.000      58.013     143.560\n",
      "x1             0.7207      0.932      0.773      0.444      -1.167       2.608\n",
      "==============================================================================\n",
      "Omnibus:                        6.419   Durbin-Watson:                   1.789\n",
      "Prob(Omnibus):                  0.040   Jarque-Bera (JB):                5.386\n",
      "Skew:                           0.880   Prob(JB):                       0.0677\n",
      "Kurtosis:                       3.371   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "from scipy import stats\n",
    "\n",
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9bedd9ee-4603-4577-977e-5dc343ef1b24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.72069731]\n"
     ]
    }
   ],
   "source": [
    "model = LinearRegression().fit(np.arange(n).reshape(-1, 1), n_epe)\n",
    "print(model.coef_)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a2f9b3b-5ad3-419d-8c8e-d37bb199fcbb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "cecfee69-7ef2-4362-9eae-b5df7761cd24",
   "metadata": {},
   "source": [
    "# Sub2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0afb5be9-1aff-4392-8ba1-93ab3e07d4b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2272.0989094232436\n",
      "40\n",
      "0.8248046168337546 False no trend\n",
      "-0.12251368985787539 50.329824135776526\n"
     ]
    }
   ],
   "source": [
    "sub_id = 2\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    year_precip = year_precip[year_precip>1]\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "print(np.sum(n_epe))\n",
    "n = len(n_epe)\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "174da059-d264-4d88-85f2-ee816b1650f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.026\n",
      "Model:                            OLS   Adj. R-squared:                  0.000\n",
      "Method:                 Least Squares   F-statistic:                     1.013\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.321\n",
      "Time:                        15:46:54   Log-Likelihood:                -194.33\n",
      "No. Observations:                  40   AIC:                             392.7\n",
      "Df Residuals:                      38   BIC:                             396.0\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         65.3966      9.925      6.589      0.000      45.305      85.488\n",
      "x1            -0.4407      0.438     -1.006      0.321      -1.327       0.446\n",
      "==============================================================================\n",
      "Omnibus:                        9.232   Durbin-Watson:                   1.651\n",
      "Prob(Omnibus):                  0.010   Jarque-Bera (JB):                8.514\n",
      "Skew:                           1.096   Prob(JB):                       0.0142\n",
      "Kurtosis:                       3.551   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5de6e9a5-bc8d-4790-a0cd-b9551e6a4d56",
   "metadata": {},
   "source": [
    "# Sub3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f397ce69-71d1-4c94-a107-e66329ad7562",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2131.012547078346   37\n",
      "40\n",
      "0.33352610308834585 False no trend\n",
      "-0.28104336005880537 47.46489045479774\n"
     ]
    }
   ],
   "source": [
    "sub_id = 3\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    year_precip = year_precip[year_precip>1]\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "print(np.sum(n_epe), ' ', len(ext_days[sub_id]))\n",
    "n = len(n_epe)\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a8550e5d-3ddf-4867-9379-0bfdd2fc4b84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.023\n",
      "Model:                            OLS   Adj. R-squared:                 -0.003\n",
      "Method:                 Least Squares   F-statistic:                    0.8785\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.355\n",
      "Time:                        15:46:55   Log-Likelihood:                -191.19\n",
      "No. Observations:                  40   AIC:                             386.4\n",
      "Df Residuals:                      38   BIC:                             389.8\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const         60.6755      9.175      6.613      0.000      42.102      79.249\n",
      "x1            -0.3795      0.405     -0.937      0.355      -1.199       0.440\n",
      "==============================================================================\n",
      "Omnibus:                        9.743   Durbin-Watson:                   1.152\n",
      "Prob(Omnibus):                  0.008   Jarque-Bera (JB):                9.313\n",
      "Skew:                           1.158   Prob(JB):                      0.00950\n",
      "Kurtosis:                       3.469   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edbef3a9-1cd4-486a-b7f4-db16c7dc1938",
   "metadata": {},
   "source": [
    "# Sub4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1743472c-6b6f-4aca-9307-8efc07c8cdb8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4129.093102444584   77\n",
      "40\n",
      "0.8067106222830902 False no trend\n",
      "-0.09392418216157161 102.45366076193893\n"
     ]
    }
   ],
   "source": [
    "sub_id = 4\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    year_precip = year_precip[year_precip>1]\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "print(np.sum(n_epe), ' ', len(ext_days[sub_id]))\n",
    "n = len(n_epe)\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "445acfeb-c9c5-4033-859e-526a6fedae3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.001\n",
      "Model:                            OLS   Adj. R-squared:                 -0.026\n",
      "Method:                 Least Squares   F-statistic:                   0.02606\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.873\n",
      "Time:                        15:46:55   Log-Likelihood:                -194.15\n",
      "No. Observations:                  40   AIC:                             392.3\n",
      "Df Residuals:                      38   BIC:                             395.7\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        101.8548      9.881     10.308      0.000      81.852     121.858\n",
      "x1             0.0704      0.436      0.161      0.873      -0.812       0.953\n",
      "==============================================================================\n",
      "Omnibus:                        0.597   Durbin-Watson:                   1.508\n",
      "Prob(Omnibus):                  0.742   Jarque-Bera (JB):                0.700\n",
      "Skew:                           0.246   Prob(JB):                        0.705\n",
      "Kurtosis:                       2.578   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21af63bb-acba-4117-986d-523519379f61",
   "metadata": {},
   "source": [
    "# Sub5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "17316867-edba-41f4-8d1d-5c94af401a42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6747.659112155144   118\n",
      "40\n",
      "0.7887194899032053 False no trend\n",
      "0.15573527436594942 157.04395857235403\n"
     ]
    }
   ],
   "source": [
    "sub_id = 5\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    year_precip = year_precip[year_precip>1]\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "print(np.sum(n_epe), ' ', len(ext_days[sub_id]))\n",
    "n = len(n_epe)\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0826d740-cb00-48bb-8cd1-cfe97e0f321c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.003\n",
      "Model:                            OLS   Adj. R-squared:                 -0.024\n",
      "Method:                 Least Squares   F-statistic:                    0.1029\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.750\n",
      "Time:                        15:46:55   Log-Likelihood:                -204.72\n",
      "No. Observations:                  40   AIC:                             413.4\n",
      "Df Residuals:                      38   BIC:                             416.8\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        165.1388     12.867     12.834      0.000     139.090     191.187\n",
      "x1             0.1822      0.568      0.321      0.750      -0.967       1.332\n",
      "==============================================================================\n",
      "Omnibus:                        3.645   Durbin-Watson:                   1.839\n",
      "Prob(Omnibus):                  0.162   Jarque-Bera (JB):                2.145\n",
      "Skew:                           0.332   Prob(JB):                        0.342\n",
      "Kurtosis:                       2.080   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a81ed5fa-dbee-4691-b391-2e6961d423a7",
   "metadata": {},
   "source": [
    "# Sub6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d6a7444b-a4a7-4ce3-9deb-cdecad390d57",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10107.429640396474   124\n",
      "40\n",
      "0.06395197805243646 False no trend\n",
      "1.7033664485416855 215.42671109966136\n"
     ]
    }
   ],
   "source": [
    "sub_id = 6\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    year_precip = year_precip[year_precip>1]\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "print(np.sum(n_epe), ' ', len(ext_days[sub_id]))\n",
    "n = len(n_epe)\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7f6e7289-abe5-4efc-a6a2-6943c50fd404",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.026\n",
      "Model:                            OLS   Adj. R-squared:                  0.000\n",
      "Method:                 Least Squares   F-statistic:                     1.019\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.319\n",
      "Time:                        15:46:55   Log-Likelihood:                -226.17\n",
      "No. Observations:                  40   AIC:                             456.3\n",
      "Df Residuals:                      38   BIC:                             459.7\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        233.5764     21.998     10.618      0.000     189.044     278.108\n",
      "x1             0.9800      0.971      1.009      0.319      -0.985       2.945\n",
      "==============================================================================\n",
      "Omnibus:                        8.173   Durbin-Watson:                   1.893\n",
      "Prob(Omnibus):                  0.017   Jarque-Bera (JB):                7.001\n",
      "Skew:                           0.862   Prob(JB):                       0.0302\n",
      "Kurtosis:                       4.107   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a78947ff-cd92-4b55-8392-0d397089fcc2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8b08e6ce-83ab-4bd7-8008-d5a3e64d360b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.97996876]\n"
     ]
    }
   ],
   "source": [
    "model = LinearRegression().fit(np.arange(n).reshape(-1, 1), n_epe)\n",
    "print(model.coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1491c720-7062-406f-9cc8-4bebfa4a027c",
   "metadata": {},
   "source": [
    "## category"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e3bcbdd1-8192-42c0-92d7-05527fc8d00b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_id = 6\n",
    "surge_induced_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_surgeWINDOW_induced_flag.npy\")\n",
    "# tutt_induced_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_TUTT_induced_flag.npy\")\n",
    "tutt_induced_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_TUTT_induced-10degree.npy\")\n",
    "# tc_induced_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_tc_induced_flag.npy\")\n",
    "tc_induced_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_tc_corrected_flag.npy\")\n",
    "MCS_induced_flag = np.load(\"ERA5/PAU/MCS_\" + str(sub_id) + \".npy\")\n",
    "MCS_induced_flag[MCS_induced_flag < 0] = 0\n",
    "midtro_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_midTropo.npy\")\n",
    "front_flag = np.load(\"ERA5/PAU/\" + str(sub_id) + \"_front.npy\")\n",
    "start, end = start_end_time(sub_id)\n",
    "all_time = ext_days[sub_id]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "085dda5b-06e1-4f6d-bb1b-1a103f737a2c",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'find_categories_events' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_100533/901669615.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# TC\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mn_epe\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfind_categories_events\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtc_induced_flag\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstart\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_epe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_epe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'find_categories_events' is not defined"
     ]
    }
   ],
   "source": [
    "# TC\n",
    "n_epe = find_categories_events(tc_induced_flag, end, start)\n",
    "plt.plot(n_epe)\n",
    "plt.show()\n",
    "n = len(n_epe)\n",
    "print(n)\n",
    "print(np.sum(n_epe), ' Sum n_epe')\n",
    "print(np.sum(tc_induced_flag), ' Sum Category')\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05/6)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0904828d-8c4f-45a4-aa0b-0a02dbdf6281",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37fb69ea-f361-43c3-b249-aa5ce3a6cbfa",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4243563b-8c0a-4898-af64-4285720c3e84",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f425c1-51c7-4c13-a1af-8026361a0e9e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "afddbc7f-f23e-4ad4-bded-1e795061c50d",
   "metadata": {},
   "source": [
    "# Sub7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fab56b05-d4c6-41af-a2d7-9c597314ae89",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20041.438409143604   153\n",
      "40\n",
      "0.5680669608513922 False no trend\n",
      "-1.0839621797581709 525.462273121139\n"
     ]
    }
   ],
   "source": [
    "sub_id = 7\n",
    "n_epe = []\n",
    "precip = monsoon_precip.sel(sub_id=sub_id)\n",
    "precip_data = precip.data\n",
    "q1 = np.quantile(precip_data[precip_data > 1], 0.05)\n",
    "q2 = np.quantile(precip_data[precip_data > 1], 0.95)\n",
    "\n",
    "for year in range(1979, 2019):\n",
    "    year_precip = precip.sel(time=precip.time.dt.year==year)\n",
    "    total_precip = year_precip.sum()\n",
    "    # print(len(ext_time), q2)\n",
    "    n_epe.append(total_precip.data)\n",
    "print(np.sum(n_epe), ' ', len(ext_days[sub_id]))\n",
    "n = len(n_epe)\n",
    "n_epe = np.array(n_epe)\n",
    "s = mk_score(n_epe, n)\n",
    "var_s = variance_s(n_epe, n)\n",
    "Tau = s/(.5*n*(n-1))\n",
    "\n",
    "z = z_score(s, var_s)\n",
    "p, h, trend = p_value(z, alpha=0.05)\n",
    "slope, intercept = sens_slope(n_epe)\n",
    "print(p, h, trend)\n",
    "print(slope, intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "29edee06-0ecc-4ceb-bb8d-0f78f33f24b5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.019\n",
      "Model:                            OLS   Adj. R-squared:                 -0.007\n",
      "Method:                 Least Squares   F-statistic:                    0.7336\n",
      "Date:                Fri, 02 Jun 2023   Prob (F-statistic):              0.397\n",
      "Time:                        15:47:26   Log-Likelihood:                -242.85\n",
      "No. Observations:                  40   AIC:                             489.7\n",
      "Df Residuals:                      38   BIC:                             493.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const        525.6423     33.386     15.745      0.000     458.057     593.228\n",
      "x1            -1.2619      1.473     -0.856      0.397      -4.244       1.721\n",
      "==============================================================================\n",
      "Omnibus:                        3.944   Durbin-Watson:                   1.430\n",
      "Prob(Omnibus):                  0.139   Jarque-Bera (JB):                2.677\n",
      "Skew:                           0.486   Prob(JB):                        0.262\n",
      "Kurtosis:                       3.813   Cond. No.                         44.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "X = np.arange(n).reshape(-1, 1)\n",
    "y = n_epe\n",
    "\n",
    "X2 = sm.add_constant(X)\n",
    "est = sm.OLS(y, X2)\n",
    "est2 = est.fit()\n",
    "print(est2.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "debec065-af57-4583-9234-ff7c53c0696d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "duan",
   "language": "python",
   "name": "duan"
  },
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}