{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# VAR模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 原理讲解\n",
    "\n",
    "可以参考公众号：[Stata: VAR (向量自回归) 模型](https://mp.weixin.qq.com/s/OxG-sk1MJUB8HV9V6y6Byw)\n",
    "\n",
    "我们常常同时关心几个经济变量的预测，比如 GDP 增长率与失业率。一种方法是用单变量时间序列的方法对每个变量分别作预测。另一种方法则是将这些变量放在一起，作为一个系统来预测，以使得预测相互自洽（ mutually consistent），这称为\"多变量时间序列\"（multivariate time series）。由 Sims（1980）所提倡的\"向量自回归\"（ Vector Autoregression，简记 VAR）正是这样一种方法。\n",
    "假设有两个时间序列变量 {$y_{11}, y_{2 t}$}, 分别作为两个回归方程的被解释变量, 而解释变量为这两个变量的 $p$ 阶滞后值, 构成一个二元 ( bivariate) 的 $\\operatorname{VAR}(p)$ 系统：\n",
    "\n",
    "$$\\left\\{\\begin{array}{l}\n",
    "y_{1 t}=\\beta_{10}+\\beta_{11} y_{1, t-1}+\\cdots+\\beta_{1 p} y_{1, t-p}+\\gamma_{11} y_{2, t-1}+\\cdots+\\gamma_{1 p} y_{2, t-p}+\\varepsilon_{1 t} \\\\\n",
    "y_{2 t}=\\beta_{20}+\\beta_{21} y_{1, t-1}+\\cdots+\\beta_{2 p} y_{1, t-p}+\\gamma_{21} y_{2, t-1}+\\cdots+\\gamma_{2 p} y_{2, t-p}+\\varepsilon_{2 t}\n",
    "\\end{array}\\right.\\tag{1}$$\n",
    "\n",
    "其中 {$\\varepsilon_{1t}$} 与 {$\\varepsilon_{2t}$} 均为白噪声过程，但允许两个方程的扰动项之间存在同期相关性\n",
    "\n",
    "$$\\operatorname{Cov}\\left(\\varepsilon_{1 t}, \\varepsilon_{2 s}\\right)=\\left\\{\\begin{array}{c}\n",
    "\\sigma_{12}, \\text { 若 } t=s \\\\0, \\text { ,其他 }\\end{array}\\right.\\tag{2}$$\n",
    "\n",
    "注意到上面两个方程的解释变量完全一样.将两个方程写在一起:\n",
    "\n",
    "$$\\left(\\begin{array}{l}\n",
    "y_{1 t} \\\\\n",
    "y_{2 t}\n",
    "\\end{array}\\right)=\\left(\\begin{array}{l}\n",
    "\\beta_{10} \\\\\n",
    "\\beta_{20}\n",
    "\\end{array}\\right)+\\left(\\begin{array}{ll}\n",
    "\\beta_{11} & \\gamma_{11} \\\\\n",
    "\\beta_{21} & \\gamma_{21}\n",
    "\\end{array}\\right)\\left(\\begin{array}{l}\n",
    "y_{1, t-1} \\\\\n",
    "y_{2, i-1}\n",
    "\\end{array}\\right)+\\cdots+\\left(\\begin{array}{ll}\n",
    "\\beta_{1 p} & \\gamma_{1 p} \\\\\n",
    "\\beta_{2 p} & \\gamma_{2 p}\n",
    "\\end{array}\\right)\\left(\\begin{array}{l}\n",
    "y_{1, t-p} \\\\\n",
    "y_{2, t-p}\n",
    "\\end{array}\\right)+\\left(\\begin{array}{c}\n",
    "\\varepsilon_{1 t} \\\\\n",
    "\\varepsilon_{2 t}\n",
    "\\end{array}\\right)\\tag{3}$$\n",
    "\n",
    "记 $\\boldsymbol{y}_{t} \\equiv\\left(\\begin{array}{l}y_{1 t} \\\\ y_{2 t}\\end{array}\\right), \\boldsymbol{\\varepsilon}_{t} \\equiv\\left(\\begin{array}{c}\\boldsymbol{\\varepsilon}_{1 t} \\\\ \\boldsymbol{\\varepsilon}_{2 t}\\end{array}\\right)$ ，则有：\n",
    "\n",
    "$$\\boldsymbol{y}_{t}=\\underbrace{\\left(\\begin{array}{l}\n",
    "\\beta_{10} \\\\\n",
    "\\beta_{20}\n",
    "\\end{array}\\right)}_{\\Gamma_{0}}+\\underbrace{\\left(\\begin{array}{ll}\n",
    "\\beta_{11} & \\gamma_{11} \\\\\n",
    "\\beta_{21} & \\gamma_{21}\n",
    "\\end{array}\\right)}_{\\Gamma_{1}} \\boldsymbol{y}_{t-1}+\\cdots+\\underbrace{\\left(\\begin{array}{ll}\n",
    "\\beta_{1 p} & \\gamma_{1 p} \\\\\n",
    "\\beta_{2 p} & \\gamma_{2 p}\n",
    "\\end{array}\\right)}_{\\Gamma_{p}} \\boldsymbol{y}_{t-p}+\\boldsymbol{\\varepsilon}_{t}\\tag{4}$$\n",
    "\n",
    "这个形式与 AR($p$) 很相似，故名“ VAR（p）”。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 滞后阶数选择\n",
    "在进行VAR建模时，需要确定变量的滞后阶数，以及VAR系统中包含几个变量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 信息准则\n",
    "根据残差 $\\hat{\\varepsilon}_{t}$ 可以估计协方差矩阵 $\\boldsymbol{\\Sigma},$ 记为 $\\hat{\\boldsymbol{\\Sigma}}_{\\circ}$ 矩阵 $\\hat{\\boldsymbol{\\Sigma}}$ 的 $(i, j)$ 元素为 $\\hat{\\mathbf{\\Sigma}}_{i j} \\equiv \\frac{1}{T} \\sum_{i=1}^{T} \\hat{\\varepsilon}_{i i} \\hat{\\varepsilon}_{j t},$ 其中 $T$ 为样本容量。则 $\\mathrm{VAR}$ 模型的 AIC 与 BIC 分别为：\n",
    "\n",
    "$$\\operatorname{AIC}(p) \\equiv \\ln |\\hat{\\mathbf{\\Sigma}}|+n(n p+1) \\frac{2}{T}\\tag{5}$$\n",
    "\n",
    "$$\\operatorname{BIC}(p) \\equiv \\ln |\\hat{\\mathbf{\\Sigma}}|+n(n p+1) \\frac{\\ln T}{T}\\tag{6}$$\n",
    "\n",
    "其中, $n$ 为 VAR 系统中变量的个数，p 为滞后阶数 $,|\\hat{\\mathbf{\\Sigma}}|$ 为 $\\hat{\\mathbf{\\Sigma}}$ 的行列式 $,$ 而 $n(n p+1)$ 为 $\\operatorname{VAR}$ 模型中待估系数之总数 (每个方程共有 ( $n p+1$ ) 个系数, 共有 $n$ 个方程)。将 | $\\hat{\\mathbf{\\Sigma}} |$ 视为多维情形下的残差平方和,则以上表达式为单一方程的信息准则向多方程情形的推广。比如，当 $n=1$ 时,$\\hat{\\mathbf{\\Sigma}}$ 为 1 阶矩阵,故 $|\\hat{\\mathbf{\\Sigma}}|=\\hat{\\mathbf{\\Sigma}}=\\frac{1}{T} \\sum_{t=1}^{T} \\hat{\\varepsilon}_{t}^{2}=\\operatorname{SSR} / T,$ 故 $\\operatorname{AIC}(p)=\\ln (\\operatorname{SSR} / T)+(p+1) \\frac{2}{T},$ 这正是单一方程的 AIC 表达式(其中,p + 1 为解释变量个数,含常数项)。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 检验最后一阶系数的显著性\n",
    "类似于由大到小的序贯 t 规则\n",
    "\n",
    "假设要确定使用 VAR( $p$ ) 还是 VAR $(p-1),$ 则可以检验原假设“ $H_{0}: \\beta_{1 p}=\\beta_{2 p}=\\gamma_{1p}=\\gamma_{2 p}=0 \"$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 检验VAR模型的残差是否为白噪声\n",
    "方法之三是检验 VAR 模型的残差是否为白噪声，即是否存在自相关。如果真实模型为 VAR( $p$ ) ,但被错误地设置为 $\\operatorname{Var}(p-1),$ 则解释变量的最后一阶滞后 $\\boldsymbol{y}_{t-p}$ 被纳入扰动项 $\\boldsymbol{\\varepsilon}_{t},$，导致扰动项出现自相关。更糟糕的是，由于 $y_{t-p}$ 的相关性, 包含 $\\boldsymbol{y}_{t-p}$ 的扰动项 $\\boldsymbol{\\varepsilon}_{t}$ 将与解释变量 {$ \\boldsymbol{y}_{t-1}, \\cdots\\boldsymbol{y}_{t-(p-1)} $} 相关, 导致 OLS 估计不一致。为此,需要检验 $\\mathrm{VAR}$ 模型的残差是否存在自相关。如果存在自相关,则可能意味着应该在解释变量中加入更高阶的滞后变量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### VAR 变量个数的选择\n",
    "在设定 VAR 模型是，主要应根据经济理论来确定哪些变量应在 VAR 模型中。比如，经济理论告诉我们，通货膨胀率、失业率、短期利息率互相关联，可以构成一个三变量的 VAR 模型。如果 VAR 模型包含不相关的变量，则会增大估计量方差，降低预测能力。另外，也可以在 VAR 系统中引入其他外生变量，比如 {$ z_{1:}, z_{2 t}, \\cdots, z_{k t},$} 与扰动项不相关。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 脉冲响应函数\n",
    "由于 VAR 模型包含许多参数，而这些参数的经济意义很难解释，故将注意力集中于脉冲响应函数。为了研究 VAR 的脉冲响应函数，首先定义向量移动平均过程。将 MA（∞）向多维推广，可以定义 n 维“无穷阶向量移动平均过程“（Vector Moving Average，简记 VMA（∞）：\n",
    "\n",
    "$$\\boldsymbol{y}_{t}=\\boldsymbol{\\alpha}+\\boldsymbol{\\psi}_{0} \\boldsymbol{\\varepsilon}_{t}+\\boldsymbol{\\psi}_{1} \\boldsymbol{\\varepsilon}_{t-1}+\\boldsymbol{\\psi}_{2} \\boldsymbol{\\varepsilon}_{t-2}+\\cdots=\\boldsymbol{\\alpha}+\\sum_{j=0}^{\\infty} \\boldsymbol{\\psi}_{j} \\boldsymbol{\\varepsilon}_{t-j}\\tag{7}$$\n",
    "\n",
    "其中 $\\boldsymbol{\\psi}_{0}=\\boldsymbol{I}_{n}, \\boldsymbol{\\psi}_{j}$ 皆为 $n$ 维方阵。定义“多维滤波”( multivariate filter)为(无穷项)“满后矩阵多项式\" (lag matrix polynomial)：\n",
    "\n",
    "$$\\boldsymbol{\\psi}(L) \\equiv \\boldsymbol{\\psi}_{0}+\\boldsymbol{\\psi}_{1} L+\\boldsymbol{\\psi}_{2} L^{2}+\\cdots\\tag{8}$$\n",
    "\n",
    "因此,可以将 VMA (∞) 简洁地写为 $y_{t}=\\boldsymbol{\\alpha}+\\boldsymbol{\\psi}(L) \\boldsymbol{\\varepsilon}_{t \\circ}$ 对于两个多维滤波,可以类似于一维滤波那样地定义其乘积，以及多维滤波 $\\boldsymbol{\\psi}(L)$ 的“逆” $(\\text { inverse }) \\psi(L)^{-1}$ \n",
    "\n",
    "使用滞后算子,可以把上述 VAR( $p$ )系统“ $y_{t}=\\Gamma_{0}+\\Gamma_{1} y_{t-1}+\\cdots+\\Gamma_{p} y_{t-p}+\\varepsilon_{t}$ \"写成 $\\operatorname{VMA}(\\infty)$ 的形式：\n",
    "\n",
    "$$\\left(I-\\Gamma_{1} L-\\cdots-\\Gamma_{p} L^{p}\\right) y_{t}=\\Gamma_{0}+\\varepsilon_{t}\\tag{9} $$\n",
    "\n",
    "$$\\Gamma(L) y_{t}=\\Gamma_{0}+\\varepsilon_{t}\\tag{10}$$\n",
    "\n",
    "其中, $\\boldsymbol{\\Gamma}(L) \\equiv \\boldsymbol{I}-\\boldsymbol{\\Gamma}_{1} L-\\cdots-\\boldsymbol{\\Gamma}_{p} L^{\\prime \\prime} \\circ$ 在方程 (10) 两边同时左乘 $\\boldsymbol{\\Gamma}(L)^{-1}$ 可得\n",
    "\n",
    "$$\\boldsymbol{y}_{t}=\\boldsymbol{\\Gamma}(L)^{-1} \\boldsymbol{\\Gamma}_{0}+\\boldsymbol{\\Gamma}(L)^{-1} \\boldsymbol{\\varepsilon}_{t}\\tag{11}$$\n",
    "\n",
    "记 $\\Gamma(L)^{-1} \\equiv \\boldsymbol{\\psi}(L)=\\boldsymbol{\\psi}_{0}+\\boldsymbol{\\psi}_{1} L+\\boldsymbol{\\psi}_{2} L^{2}+\\cdots, \\boldsymbol{\\Gamma}(L)^{-1} \\boldsymbol{\\Gamma}_{0} \\equiv \\boldsymbol{\\alpha},$ 则得到 $\\mathrm{VAR}$ 模型的 $\\mathrm{VMA}$ 表示法(Vector Moving Average Representation):\n",
    "\n",
    "$$y_{i}=\\alpha+\\psi_{0} \\varepsilon_{t}+\\psi_{1} \\varepsilon_{t-1}+\\psi_{2} \\varepsilon_{t-2}+\\cdots=\\alpha+\\sum_{i=0}^{\\infty} \\psi_{i} \\varepsilon_{t-i}\\tag{12}$$\n",
    "\n",
    "其中, $\\boldsymbol{\\psi}_{0} \\equiv \\boldsymbol{I}_{n},$ 而其余的 $\\boldsymbol{\\psi}_{i}$ 可通过递推公式 $\\boldsymbol{\\psi}_{i}=\\sum_{j=1}^{i} \\boldsymbol{\\psi}_{i-j} \\boldsymbol{\\Gamma}_{j}$ 来确定。\n",
    "\n",
    "推导如下：\n",
    "\n",
    "递推公式 $\\psi_{i}=\\sum_{j=1}^{i} \\psi_{i-j} \\Gamma_{j}$ 证明:\n",
    "\n",
    "利用 $\\left(\\psi_{0}+\\psi_{1} L+\\psi_{2} L^{2}+\\cdots\\right) \\Gamma(L)=I_{n}$   \n",
    "\n",
    "$\\Rightarrow\\left(\\psi_{0}+\\psi_{1} L+\\psi_{2} L^{2}+\\cdots\\right)\\left(I-\\Gamma_{1} L-\\cdots-\\Gamma_{p} L^{p}\\right)=I_{n}$  \n",
    "\n",
    "$\\Rightarrow \\psi_{0}-\\psi_{0} \\Gamma_{1} L-\\psi_{0} \\Gamma_{2} L^{2} \\cdots+\\psi_{1} L-\\psi_{1} \\Gamma_{1} L^{2} \\cdots=I_{n}$  \n",
    "\n",
    "$\\Rightarrow \\psi_{0}+\\left(\\psi_{1}-\\psi_{0} \\Gamma_{1}\\right) L \\cdots=I_{n}$  \n",
    "\n",
    "$\\Rightarrow \\psi_{1}-\\psi_{0} \\Gamma_{1}=0$\n",
    "\n",
    "根据向量微分法则可知：\n",
    "\n",
    "$$\\frac{\\partial \\boldsymbol{y}_{t+s}}{\\partial \\boldsymbol{\\varepsilon}_{t}^{\\prime}}=\\boldsymbol{\\psi}_{s}\\tag{13}$$\n",
    "\n",
    "其中, ( $\\left.\\partial \\boldsymbol{y}_{t+s} / \\boldsymbol{\\partial} \\boldsymbol{\\varepsilon}_{t}^{\\prime}\\right)$ 为 $n$ 维列向量 $\\boldsymbol{y}_{t+s}$ 对 $n$ 维行向量$\\boldsymbol{\\varepsilon}_{t}$ 求偏导数,故得到 $n \\times n$ 矩阵 $\\boldsymbol{\\psi}$ 。矩阵 $\\boldsymbol{\\psi},$ 是一维情形下相隔 s 期的动态乘子向多维的推广,其第 $i$ 行、第 j 列元素等于 $\\left(\\partial y_{i, t+s} / \\partial \\varepsilon_{i t}\\right)$ 。它表示的是，当第 j 个变量在第 $t$ 期的扰动项 $\\varepsilon_{jt}$ 增加 1 单位时(而其他变量与其他期的扰动项均不变) 对第 i 个变量在第 $(t+s)$ 期的取值 $y_{i, t+s}$ 的影响。将 $\\left(\\partial y_{i, t+s} / \\partial \\varepsilon_{j t}\\right)$ 视为时间间隔 $s$ 的函数, 就是脉冲响应函数” (IRF)。\n",
    "\n",
    "脉冲响应函数的缺点是,它假定在计算 ( $\\left.\\partial y_{i, t+s} / \\partial \\varepsilon_{j t}\\right)$ 时,只让 $\\varepsilon_{jt}$ 变动，而所有其他同期扰动项均不变。此假定只有当扰动项的协方差矩阵 $\\boldsymbol{\\Sigma} \\equiv \\mathrm{E}\\left(\\boldsymbol{\\varepsilon}_{t} \\boldsymbol{\\varepsilon}_{t}^{\\prime}\\right)$ 为对角矩阵时才成立（即同期扰动项之间正交）;否则，当 $\\varepsilon_{jt}$ 变动时，必然伴随着其他方程的同期扰动项发生相应的变动。为此,需要 考虑“正交化的脉冲响应函数”( Orthogonalized Impulse Response Function,简记 OIRF)，即从扰动项 $\\boldsymbol{\\varepsilon}_{t}$ 中分离出相互正交的部分,记为 $\\boldsymbol{u}_{t},$ 然后计算当 $\\boldsymbol{u}_{t}$ 中的某个分量变动时,对各变量在不时期的影响。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 原理实现\n",
    "#### VAR 模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>const</th>\n",
       "      <th>y_1</th>\n",
       "      <th>y_2</th>\n",
       "      <th>x_1</th>\n",
       "      <th>x_2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>53.256500</td>\n",
       "      <td>1.719442</td>\n",
       "      <td>-0.714049</td>\n",
       "      <td>-1.019274</td>\n",
       "      <td>0.993932</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>39.431266</td>\n",
       "      <td>0.241814</td>\n",
       "      <td>0.738772</td>\n",
       "      <td>0.564473</td>\n",
       "      <td>-0.560130</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       const       y_1       y_2       x_1       x_2\n",
       "0  53.256500  1.719442 -0.714049 -1.019274  0.993932\n",
       "1  39.431266  0.241814  0.738772  0.564473 -0.560130"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "\n",
    "data = pd.read_excel('../数据/上证指数与沪深300.xlsx')\n",
    "Y = data['hs300']\n",
    "X = data['sz']\n",
    "\n",
    "def lag_list(Y, X, p=1, q=1):\n",
    "    '''\n",
    "    待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q)\n",
    "    获取自回归分布滞后模型的估计向量\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    Y : 被估计变量\n",
    "    X : 估计变量\n",
    "    p : ADL 模型 Y 的滞后阶数，标量默认为1\n",
    "    q : ADL 模型 X 的滞后阶数，标量默认为1\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ADLy : ADL 模型被解释变量\n",
    "    ADLx : ADL 模型解释变量\n",
    "\n",
    "    '''\n",
    "    ADLx = pd.DataFrame()\n",
    "    T = len(Y)\n",
    "    ADLy = list(Y[max(p, q):T])\n",
    "    for i in range(1, p+1):\n",
    "        name = f'y_{i}'\n",
    "        ADLx[name] = list(Y[max(p, q)-i:T-i])\n",
    "    for i in range(1, q+1):\n",
    "        name = f'x_{i}'\n",
    "        ADLx[name] = list(X[max(p, q)-i:T-i])\n",
    "    return ADLy, ADLx\n",
    "\n",
    "def VAR(Y, X, lag):\n",
    "    ADLy, ADLx = lag_list(Y, X ,p=lag, q=lag)\n",
    "    ADLx = sm.add_constant(ADLx)\n",
    "    mod = sm.OLS(ADLy, ADLx)\n",
    "    a = mod.fit().params\n",
    "    ADLy, ADLx = lag_list(X, Y ,p=lag, q=lag)\n",
    "    ADLx = sm.add_constant(ADLx)\n",
    "    mod = sm.OLS(ADLy, ADLx)\n",
    "    b = mod.fit().params\n",
    "    return pd.DataFrame([a,b])\n",
    "\n",
    "VAR(Y, X, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## statsmodels 库实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  Summary of Regression Results   \n",
       "==================================\n",
       "Model:                         VAR\n",
       "Method:                        OLS\n",
       "Date:           Sun, 31, May, 2020\n",
       "Time:                     19:38:46\n",
       "--------------------------------------------------------------------\n",
       "No. of Equations:         2.00000    BIC:                    11.9914\n",
       "Nobs:                     458.000    HQIC:                   11.9368\n",
       "Log likelihood:          -4015.14    FPE:                    147453.\n",
       "AIC:                      11.9013    Det(Omega_mle):         144286.\n",
       "--------------------------------------------------------------------\n",
       "Results for equation sz\n",
       "===========================================================================\n",
       "              coefficient       std. error           t-stat            prob\n",
       "---------------------------------------------------------------------------\n",
       "const           39.431266        19.842955            1.987           0.047\n",
       "L1.sz            0.241814         0.201697            1.199           0.231\n",
       "L1.hs300         0.564473         0.148984            3.789           0.000\n",
       "L2.sz            0.738772         0.202499            3.648           0.000\n",
       "L2.hs300        -0.560130         0.149698           -3.742           0.000\n",
       "===========================================================================\n",
       "\n",
       "Results for equation hs300\n",
       "===========================================================================\n",
       "              coefficient       std. error           t-stat            prob\n",
       "---------------------------------------------------------------------------\n",
       "const           53.256500        26.908303            1.979           0.048\n",
       "L1.sz           -1.019274         0.273514           -3.727           0.000\n",
       "L1.hs300         1.719442         0.202032            8.511           0.000\n",
       "L2.sz            0.993932         0.274602            3.620           0.000\n",
       "L2.hs300        -0.714049         0.203001           -3.517           0.000\n",
       "===========================================================================\n",
       "\n",
       "Correlation matrix of residuals\n",
       "               sz     hs300\n",
       "sz       1.000000  0.973265\n",
       "hs300    0.973265  1.000000\n",
       "\n"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from statsmodels.tsa.vector_ar.var_model import VAR\n",
    "import pandas as pd\n",
    "\n",
    "data = pd.read_excel('../数据/上证指数与沪深300.xlsx')\n",
    "estimate_data = data[['sz', 'hs300']]\n",
    "res = VAR(estimate_data).fit(maxlags=5, method='ols', ic='bic')\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "对上证指数的脉冲示意图\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqoAAAK8CAYAAAAnCp3sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3yV5f3/8dcne0AWI0SGQNgExIKgqKgFFaWKq9ZRtzjAUVtb+7V1fPVr66q/atVapcVRrYh7VmVogCoURCSALJkSZgIEspPr98d9zu05GRAgkAN5Px+P+3Fy39c9rnNyNG+u+7qu25xziIiIiIhEmqimroCIiIiISF0UVEVEREQkIimoioiIiEhEUlAVERERkYikoCoiIiIiEUlBVUREREQikoKqiNRiZs7MLmjqesgPzCzdzDaaWfYBOv+9ZpZ3IM7dFMysc+B7PKiu9QN0zUGBa3QOrP/EzL42M/2tFdlH+o9HJAKY2fNm9n5T10Mi2p3Ah865Fft7omb6D5G1QBbw9cG6oHPufaAKuPRgXVPkcKOgKiJNzszimroOkczMkoBrgb/v53ma7efsnKtyzm1wzlUe5EtPAG45yNcUOWwoqIpEoGALq5ndYWYbzGy7mT1oZlGBW7SbAtvvqHGcM7ObzOwDMys2s9Vm9vOQ8jpvf+6phc3M7g6cqyxw3RdDyszMfmNmK8ysxMwWhF6zAe9vHbAusD3OzB4ys3VmtsvM/mtmp4ccF2tmT5jZ+kBd1prZgyHlqwKfzz/NbGegrrfXuHYnM3vLzIoCy5tm1iGk/F4zyzOziwLvqcjM3jaz1iH79DOzKWa2I1A+38xOCSnvE/gdFAV+V/8ys3YNPb4OZwLVwMwa72WYmc0ys9JAt4D/FxpGzewzM/urmT1qZpuBmWa2KlA8KfB7X1XjnLt731Fmdlfgcy8L/K5H1zh+iJl9FajTPDM7M3Cdk/ey3k+b2R/MbEvgM3zUQm6hm9nPA9+P4Gc8ycza1/cB1vzuB67h6lhODpTv9rsY2GekmX0beB/TgR51XPpdYJCZdauvbiJSPwVVkcg1DOgCnAzcAPwG+BCIB04A7gUeNLOBNY77X7w/jgOAZ4EXbT/65ZnZ+cDtwFigO/ATYHbILv8HXAOMA/oAfwT+Zmaj9nDqk4D+wEhgeGDbhMD2S4B+wAvAe2Z2VKD8FuBc4KJAXX4GLKlx3l8Ci4EfAfcAfzCz8wLvxYC3gUzgx8ApwBHA24GyoM6Bc58LnAYcDTwQUv4KkA8MDpTdC5QGrpEF5AJ5gfIRQAvg3ZCgVe/x9TgRmOtCnnkdCGUfAfMC57gGuBjv8w/1c8AC57gcOCawfQzerfBjQvbd0/u+Ffg1cAfe7+ct4E0zGxCoUwvgfeBbYCDed/aR0MrsRb0vBSqBocBNwC8CdQuKw/v9HoX3nWwN/IuGOy/w/oPLM8DGQN1hD99FM+uI9136FO+/tb8AD9e8iHNuTeC8J+1F3UQkyDmnRYuWJl6A54H3a6yvBaJDts0Bvqlx3Crg9pB1BzxXY5/JwD8DP3cO7DOoxj4OuKCudbzgtwSIraPeyUAJcGKN7X/G60+5u/e7GYgP2ZaN12rYqca+bwNPB35+ApgCWD3nXQV8WmPbeGBG4OdT8foMdg4p7xq47ojA+r14oTE1ZJ/fActD1ncAV9RTh/uAKTW2pQc+08F7Or6ec74NvFBj2wPAciAqZNuVQBmQFFj/rOZ3pq7f91687++Bu2sc91nI9+t6oABIDCm/JHC9k/ey3l/UuM6nwPjdfEa9AtfpUNd3veZ6jWN/hvc9PnYvvot/AJaGfheB3weu0bnGcV8B9zf0961Fi5YfFrWoikSuRc65qpD1jcCCGvtsBNrW2PZFHet99qMek4AEYKWZ/d3Mfmpm8YGyPoGyfwdute80s53AjXh/7HcnzzlXFrL+I7yWv0U1zjUq5FzP47VeLTWzp8xslNUeUb27998bWO+cWxUsdM59B6wn/DNa7ZzbHrK+nvDP+TFgvJlNNbPfmVmvkLKBwLAa72FtoCy7AcfXJZHaLa698cJcdci2GXgtjaG3mefu4dyh6n3fZpaC1/o8s8YxM/jhs+uF93stCSmftY/1/qbGcWG/AzP7kZm9Y16XlCK8f8gBdNrN+6slcLfhH8A1zrkvA5sb8l3sDXzpnHMhp6v53QsqwfsdisheimnqCohIvSpqrLt6tu3NPziD4cC/zW1msbs7wDm31sx64t2eHwH8CbjHzIaEXPssYE2NQ2vWtaZdNdaj8N7PMXUcWxKoy1fmTf0zEu/W/QvAfDM7tUbwqY8FrlGX0O27/Zydc/ea2cvAGcDpeJ/HDc65fwT2+wCvu0RNGxtwfF224LXK7st7qfk5705Dvl91XTO4bXd1ogH7NOh3YGbJwMd4dwsuAzbh3fqfjhd4G8TMjsBrJX3MOfdKSNEev4uE/DfUABl4dxBEZC+pRVXk8HNsHeuLAz8H/1hmhZQP2NMJnXOlzrkPnHO34f3x7gscDyzCu2V7pHNueY1l9V7Wex7eH/92dZzr+5C6FDnnJjnnbsRr4fox4S1xu3v/i4D2gbALgJl1xWspXLQ3lXXOLXPOPeGcG4U3Gv/aQNFXeJ/P6jreR1EDjq/LPGq3ii8CjqvRonwCUA7saQqrCiB6T+8xlHNuB16r5gk1ik7gh89uMdDPzEJbDwc3Yr2DeuEF0zudc7nOuW+pfWdht8wsAS+kfgncXaO4Id/FRcCQGn2ba373gtfJxvteiMheUlAVOfycZ2ZjzKy7mf0PXkvonwECt2S/BO4ws75mNhR4dHcnM7Mrzexa80aqdwGuwgs6ywLB61HgUTO72sy6mdkAM7vBzK7bm0o755YCLwPPm9kFZtbVvAnUbw8ZDPVLM7vYzHoHRlFfgtffc13IqY41s/8JvP8xeAOI/l+gbDIwH3jZzAYGbvu+jBcipjaknmaWGOh2cHJgJPkQwsPaU0AqMNG8EfBdzWyEmT1rZi0bcHxdPgZ6m1mrkG1P4wXspwOfxyjgQeBJ51zxHt7GKmC4mbUzs5ottbvzCHB74HfQw8zuwxuk9adA+ct4fYCfM2/mgxF487/CD62l+1PvoDV4/0C6KfD5jgLu34v3AfA3IA1vwFdm4LNoZ2ZxDfku4g2+6gz82cx6mjdrxg11XOfYQF1rdpkQkQZQUBU5/NwLnI/Xx+9G4Crn3H9Dyq8OvP4X74/17/dwvm14I7On441kPx84zzm3MlB+V+CatwML8Qa9nA+srHWmPbsKb7T1w3ijr9/Hm/0g2DpbhDfqfDZeuBwAnFEj4DyGN5vAPLwZCe52zr0OEOhPeA5ey/JnwDRgA3BOjb6Gu1OFdxv+BbxBZm/h9U38ZeAa6/Fam6uBf+N9Jk/hhZWyPR1fF+fcgsB7vihk2/d4XQeOxpvE/h94o97vrOscNfwKb8aDtXifU0M9gRdWH8b7LpwLnO+c+zpQp5143UD6Bs77CN53AwJ9bPez3gTOsRm4Au93uQhv9H+9n189TsKbOWIF3gwMwWVooHy330XnjeY/D68bynzgNuC3dVznYuDlvQjhIhLCGv7/ZhGJdGbmgJ8Gg1lzY96coE8653bbSnwoMrORwONAnxqD7CKaefOsvgW0dc5taer6HExm1gYv5A4K+YediOwFDaYSETkEOOf+bWZPAR34oYU54pjZFcB3eK21OXjdTt5rbiE1oAswViFVZN8pqIqIHCKcc080dR0aIBPvoRNZeN0qPsB7QECz45ybTfjDMURkL+nWv4iIiIhEJA2mEhEREZGIpKAqIiIiIhFJQVVEREREIpKCqoiIiIhEJAVVEREREYlICqoiIiIiEpEUVEVEREQkIimoioiIiEhEUlAVERERkYikoCoiIiIiEUlBVUREREQikoKqiIiIiEQkBVURkQYyszvM7HszKzKzJWY23Mx+ZmY7Q5YyM/usgeezg3UtEZFDkYKqiEgDmFlP4CbgGOdcS+B0YJVzbqJzroVzrgVwBPAd8K8GnO9E4EMzSzjQ1xIROVQpqIqINEwVEA/0MbNY59wq59yKYKGZRQGvAJ855/7WgPPNBDYB79YRVhv7WiIihyQFVRGRBnDOLQd+AdwLbDKzV83siJBdHgBaAreEHmdmI83M1VzwwujlwKnAjY1xLRGRw40555q6DiIihxQzSwH+BlQ65y4zs4uAB/Fu1W9u4DmigAlAFjDaOVdyoK4lInKoimnqCoiIHAoC/Ubb492yLwVKgCgzOxr4C3DqXgbH44FM6gipB+BaIiKHJN36FxFpmHi8lswtwAagLXAnMBpIB2aEjMb/aE8nc85NB86opyW1Ua8lInKo0q1/EREREYlIalEVERERkYikoCoiIiIiEUlBVUREREQiUoODqplFm9k8M3s/sJ5hZp+a2bLAa/qBq6aIiIiINDcNHkxlZr8EBgEpzrmfmNnDQIFz7kEz+y2Q7py7Y3fnaN26tevcufP+1llEREREDhNz587d4pxrU1dZg+ZRNbMOwCi8p6H8MrB5NHBy4OcXgM+A3QbVzp07M2fOnIZcUkRERESaATNbXV9ZQ2/9/xn4DVAdsi3TOZcPEHhtW8/FrzOzOWY2Z/NmzU8tIiIiIg2zx6BqZj8BNjnn5u7LBZxzzzrnBjnnBrVpU2erroiIiIhILQ259X88cLaZnQkkAClm9k9go5llOefyzSwL2HQgKyoiIiIizcseW1Sdc//jnOvgnOsMXARMdc79HHgXuCKw2xXAOwesliIiIiLS7OzPPKoPAqea2TLg1MC6iIiIiEijaNCo/yDn3Gd4o/txzm0Fhjd+lURERERE9GQqEREREYlQCqoiIiIiEpEUVEVEREQkIimoioiIiEhE2qvBVNL0ysvLiYmJISoqitzcXN5++22WL1/OihUr2LZtGz179uSNN94gPT2dDRs2kJiYSGpqalNXW0RERGSvKahGsCVLlvDOO++wYsUKf1mzZg1Lly4lOzubuXPn8re//Y3s7Gx69OhBamoqK1eu9IPpXXfdxfjx4+nYsSM5OTnk5OTQv39/fv7znzfxOxMRERHZM3POHbSLDRo0yM2ZM+egXS/Sbd68mX//+99hQXTFihW89tprnHTSSUyaNIkLL7yQVq1akZ2d7S833ngjRxxxBBUVFcTExGBmdZ5/5syZTJ8+nby8PPLy8li8eDFHHnkkS5cuBeCGG25g8+bNfojt27cv3bt3JzY29mB+DCIiItKMmdlc59ygusrUonoAlZSU8J///McPoMFb9L/73e/46U9/yooVK7j88ssxMzp06EC3bt04++yzSUtLA2DUqFFs27at3lv3ewqUxx9/PMcff7y/XllZycaNG/316OhoFi5cyNtvv011dTUAP/7xj5kyZQoAzz33HO3atSMnJ4cjjzySqCh1aRYREZGDR0F1P1RXV/Ptt9+GhdAVK1Zw7rnnct1111FQUMCIESMAiIuLo2vXrmRnZ9OiRQsABgwYwOLFi+ncuTMJCQm1zp+UlNSo9Y2JiaF9+/b++lNPPQVAaWkp3377LXl5ebRs2RLwQu1NN91EeXk5AMnJyfTt25drr72WMWPG4Jxj48aNZGZm1tuiKyIiIrI/FFT3oLCwMOy2/PLly+nXrx+33XYb1dXVHHXUUVRWVgKQkpJCdnY2VVVVAGRlZTF16lSys7Np37490dHRYedOSEigV69eB/091ZSQkMCAAQMYMGCAvy0mJobNmzezcOFCv+tAXl6eH1w3btxIVlYWrVq18rsO5OTkcNppp9G1a9emeisiIiJyGGn2QbW6upr8/PywFtHU1FR+85vfAF6r55o1a/z927Vr59+Kj4mJ4fXXX6ddu3ZkZ2fTqlWrsNbFqKgoTjnllIP7hhpRSkoKxx13HMcdd1ytsvj4eJ544gk/wL700kvs2LGDF198ka5du7JgwQJ+/etfh4XY3r17k5yc3ATvRERERA5FzWIwVXl5OatXr/aDaHFxMb/+9a8Br0/mtGnT/H2jo6P58Y9/zCeffALAxIkTiY+PJzs7m65duypo1cM5x7p160hJSSE1NZUZM2Zw6623smjRIkpLSwEwM3JzcznhhBNYvHgx8+fPJycnhx49ehAXF9fE70BERESawu4GUx32QfXmm2/m6aef9gcLAWRmZpKfn4+Z8a9//YvCwkKys7Pp1q0bnTp10qj3RlRVVcV3333nt7yOHTuWVq1a8cADD/D73/8e8Fqme/bsSU5ODn/9619JT09n586dJCYm1uouISIiIoeXZh1UX3vtNRYsWEC3bt386Z3atWunAUBNrKysjKVLl4b1f12yZAl5eXnExMQwbtw4JkyYQJ8+ffyuA/369eP0009v6qqLiIhII2rWQVUOTR999BGffvqpH2Lz8/Pp1KkTq1evBuCmm25i2bJlpKWl+UuPHj245pprAPjyyy8BwsrrmllBREREmpbmUZVDzhlnnMEZZ5zhr2/dupX8/Hx/3czYvn07q1evZtu2bRQWFjJ48GA/qF5zzTUsWrSo1jk//PBDAM477zyKi4vDguwxxxzD+eefD8C0adNISkoKK4+Pjz/Qb1tERERCKKjKIaFVq1a0atXKX//LX/5Sa5+Kigr/5xdffJFNmzaxbds2f+nYsaNfHh0dTWFhIStXrvSD7sUXX+wH1TPPPNMfBBY0btw4nnzySaqrqznhhBNITU0NC7IjRoxg+PDhVFRUMG3atLCytLQ0DRgTERHZSwqqctgIHQQ3cODA3e47adKksHXnnD//rXOOKVOmhIXcbdu2+fPMlpaW0qJFC7Zu3cqKFSv8oJuQkMDw4cPZunVrnX1pH3nkEW6//XbWrVvHhRdeGBZi09PTOe+88xg4cCDbt29n9uzZfllKSgpJSUkkJyfr6WAiItKsKKiK4HUliImJ8X8eOnRovfsmJSX505cFOef8mSXS09OZMWNGraAbfJxtRUUFycnJbN68mWXLlvnlPXv2ZODAgSxcuJDTTjut1nVff/11zj//fKZNm8Zll11GUlJS2PLQQw8xcOBAvvrqKyZMmOBvT05OJikpifPPP5/MzEzWrVvH0qVLw45NTk6mTZs2/mcgIiISCfRXSaQRmJk/lVZ8fLwfSuvSpUsXPv3007BtoUE3JyeH6dOn+wF2x44dFBcX079/f8DrBnH66adTXFwctgStWrWKV155heLi4rDuC0OGDCEzM5MPP/yQ66+/vla9Fi9eTK9evXjiiSe46667agXd9957jzZt2vDGG2/w3nvvhZUlJSXxi1/8gvj4eObPn8+qVatqlXfv3h0zo7KykujoaM28ISIie6RR/yKHsaqqKkpKSiguLiY9PZ3Y2Fg2bNjAkiVLagXdiy66iNTUVKZNm8a7777Lrl27wsonTpxIamoqjz32GI8//nitkFxaWkp8fDy33HJLrT7E0dHRVFRUYGZcc801PP/882EtupmZmf5MDX/84x+ZPXs2SUlJJCYmkpSURLt27bjzzjsBeP/999m4cWNYeUZGBoMGeQNGN2zYgJn55WolFhGJbJqeSkQOGOccpaWlJCQkYGbk5+eTn58fFmTLysr42c9+BsA777zDnDlzwspjY2MZP348ALfffjuffPIJxcXFfsjOysryZ3EYPnw4U6dODavDUUcdxddffw14LcezZ8/2y2JjYzn55JP97hqjR49m/fr1fshNSkpiyJAh3HHHHQA8+OCDlJWV+eWJiYn06NGDE088EfCmPouNjQ07vmXLliQmJh7AT1lE5PCloCoih42CggKKioooKSnxg2xcXBzHHHMMAG+//Tbr16/3g25JSQkdOnRg7NixgDd7w8qVK/1jS0pKGDp0KM888wwAnTp1Yu3atWHXvOSSS3j55ZcBSE5ODutqATBmzBieffZZnHOkp6eTkJDgh9jExESuuOIKbrrpJkpKSrj22mv97QkJCSQmJjJ8+HCGDRvGrl27mDRpUlhZQkIC3bt3Jysri4qKCjZt2hRWpgF2InKo0zyqInLYyMjIICMjo97yc845Z7fHP/XUU7stX7NmDVVVVZSWloYF4aBgt4jQoNy7d2/A62px5ZVX+tuDQTjY2lpaWsrs2bP940tLS/3W6GHDhrFp0yauuuqqWnV6/PHHueWWW1iyZAn9+vULK4uLi2P8+PFcdtllzJs3j0suuaRW0P3d737H0KFDWbx4MU8//XRYWWJiIueddx5HHnkk33//PXPmzKl1fLdu3UhMTKS8vJzKykoFZBE5aBRURURqiI6OJjk5meTk5Fplw4cPr/e4mJgY/vznP9dbnp6ezrJly8K2hQ6k69ChA999950fYoOv3bt3ByArK4u//e1vYWUlJSX06dMHgISEBPr37++H6NLSUrZv305ZWRkA69at4+WXX/aPCzrqqKM48sgjmT59OhdffHGtes+aNYvBgwfz4osvMmbMGMALyMEwm5ubS48ePXjllVf4y1/+UivoPvHEE7Rq1YrJkyczdepU4uPjSUhI8F+vuuoq4uPjWbhwIatXrw4ri4+Pp2/fvpgZO3fuxDlHfHw8sbGxGpAn0gzsMaiaWQKQC8QH9n/dOXePmd0LjAE2B3a90zn34YGqqIjI4Sh0xojY2Fi6dOlS776tWrXiuuuuq7e8d+/eTJw4sd7yU089lYKCAsALyOXl5ZSWlpKUlATA6aefzty5c8Nae0tKSvygPHjwYB566CF/e3C/1NRUwJvxomXLlpSWlrJ582Z/v2AQ/+KLL3j00UfDHs4BcNlllxEfH8/48ePrDPrB42+77Ta/L7OZER8fT+vWrf2uGrfffjuTJ08mPj7eD7pZWVm88MILADz55JN8++23YUE4KyuLa6+9FoCPP/6Ybdu2hYXkjIwMjjrqKAC+//57gLDjNVhP5MDaYx9V8/7Jmuyc22lmscAM4FZgJLDTOfdoQy+mPqoiIlJVVUVZWRllZWWUlpbSrl07zIy1a9eyfv16f3tpaSnl5eVccMEFAEyePJn58+dTWlrq7xMdHc0f//hHAP70pz8xffp0v6ysrIz09HQ++OADwOtr/Mknn/hllZWV9O3bl7y8PACGDh3KF198EVbXIUOG+DNS9O/fnwULFoSVn3rqqf5AvaFDh5Kfn098fDxxcXHEx8czYsQIv35XXHGFPztGsHzo0KFceumlADz88MNER0f7QTs+Pp4+ffowaNAgqqurmTZtmn9ccGnTpg2tWrWiurqaXbt2qbVZDkn71UfVeUl2Z2A1NrAcvBFY+yEvDyZNgpQUaNnSew39OfRVj3EXETk4oqOj/cFmoTp27Bj2qOOaRowYwYgRI+ot/9WvfsWvfvWrestfeeWVsPWqqqqw1t3XXnuNoqKisCAcWsf77ruPTZs2hQXpI4880i8/4YQT2LBhA2VlZZSXl1NWVhZ2/IoVK9iyZYsf0svLywH8oPo///M/futx0C233MKgQYMoKyur873feeedPPDAA2zZsoXMzEx/ezDQ3nffffziF79g7dq1nHbaaX7ADZbfeuutnHXWWaxZs4a77747LAjHxcVx4YUXMmDAANavX8/bb78ddmxcXBxDhgwhKyuLwsJCli5dSlxcXFh5ZmYm8fHxVFZW4pwjJiZGIVr2SoPuWZhZNDAX6AY85ZybZWZnADeZ2eXAHOBXzrnCOo69DrgOvNG0B9OCBXDffQ3bNy6u/hBbX8Ctb/+QJ3mKiEiEio6O9rtdgNdHeHf2NFDv4Ycf3m35jBkzdlsenMotdGnZsiXgBc/PP/88LASXlZX5A/mSkpJ45JFHapUHB9/FxMTQr1+/sJBcUlJCZWUlANu3b2fatGlhx5aXl5OTk8OAAQNYunQp48aNq1Xnd955h7PPPpuZM2dy1lln1SqfOnUqp5xyCpMmTeKSSy7x30twmTx5MkcffTSvvfYa9913X1hZXFwcEyZMoGPHjnzwwQe8+uqrtYLw7373O1JTU5k5cyazZs3ytweXCy64gNjYWJYuXcr69etrnb9nz55+/+eqqip/e+j3QprWXk1PZWZpwFvAzXh9U7fgta7eD2Q5567e3fFNceu/uhp27YIdO6CoyHsN/Xlvtu3a1bBrJiTsfeitr0z/rYiISFNxzmFmVFRUUFBQUCvIdunShbS0NDZt2sScOXMoLy8PW8444wyysrLIy8vjnXfeqVX+m9/8ho4dO/LJJ5/w7LPP+tuD5584cSIdOnTgueee449//GOt41etWkXbtm256667+L//+79a9d+1axdJSUncdttttfo/mxlVVVWYGddeey1///vf/bLo6GhatWrFxo0bARg7diwfffSRH2RjY2Np37497733HgD33HMP8+bNCwvBHTt25P777wfgmWeeYd26dcTGxvrl7du356KLLgLggw8+YOfOnf654+LiaNu2LQMGDAC8Jwc658Kun5ycTEpKCuD14z6UZ+Jo1HlUzeweYFdo31Qz6wy875zL2d2xh3of1aoq2Llz/8Ju8DVkwO1uJSU1PNi2agV9+kCvXl5YFhERaQ6CLcQ1g2z37t2JiopixYoVrF27NiwIV1ZW+g8imTx5MgsWLPADcnl5OdHR0fzv//4vAE8//TSzZs0KO3d6ejovvvgi4HXRmD59OhUVFX55165d/YeTnHLKKUyfPp2qqiq/zscee6zfJ7pfv35+X+mgESNG+I/b7tKlC6tWrQorP+ecc3jrrbcAaN26Ndu2bQsLwhdffDFPPPEEAD/60Y/8oBvc57zzzuOmm26iqqqKSy+91N+2pzsHB8J+9VE1szZAhXNum5klAiOAh8wsyzmXH9jtXCCv3pMcJqKjITXVW/ZXRYUXWPcl7K5eHV4W6OYUVs/u3SEnB/r1815zciA7Wy20IiJy+AmGs/pkZ2eTnZ1db/me+j+PHTvWf2hIXYKBsD7Tpk0DfugXXVFREdYf+f3336e4uDgsCAe7fYAXlHfs2BEWhEP7R//617+mqKiI8vJyf5+BAwf65d26dfNDeLA82D+7oqKCefPmUVFREXZMpGjIqP/+wAtANBAFvOacu8/MXgIG4N36XwVcHxJc63Sot6hGqrIyL7hu3AgLF3qDyBYs8F5XrIDgrzghwWtxDetQYdQAACAASURBVAbXYIht3x7Ut11ERESagh6h2owVF8OiRV5oDQ2w69f/sE9a2g/hNTTA7ubhPyIiIiKNQo9QbcaSkmDQIG8JVVDwQ3gNBthXX4Vt237YJyurdveBPn2gjof1iIiIiDQ6BdVmKiMDhg3zliDnvJbW0JbXvDx4+mkoLfX2MYOuXWsH2B49NC2XiIiINC4FVfGZef1V27eH00//YXtVFXz3Xe0A+/77Xhl4IbVXr9rdB448Eg7hGTNERESkCamPquyz0lJYsqR2gF29+od9WrSAvn1rB9i2bTWAS0RERDSYSg6yHTtqzz6wYAFs2fLDPq1b1+4+kJPjzQcrIiIizYcGU8lBlZICxx3nLaE2bqw9+8A//hH+xK9OnWoHWD3AQEREpHlSUJWDJjPTW4YP/2FbdbXXVaDmDASffuo9FAHCH2AQDLHduv3wVK4WLSA+vmnek4iIiBw4CqrSpKKioEsXbznrrB+2V1TAsmXhfV/nzYM33vjhAQahYmO9wBoMrqE/78trcrIGgYmIiDS1wz6o7tq1i6SkJEwjdw4psbHenK19+kDgUcyA101g8WJYuRJ27vSeyFXXa/DnzZvDy4LTbDVEUtL+hd2agTk+XgPIRERE9sZhH1TvuOMOJk6cyLBhwxg2bBgnnXQS/fr1I1oPvT8kJSfX/QCDhqqsDA+ye/O6c6cXfFeuDC8LeVzzbsXE7HvIDd0WusQc9v8Fi4hIc3bY/5k79dRT2blzJ7m5ubz55psA9O3bl7y8PAC+++47OnbsSKxmq28WYmK8R8ampTXO+ZzzWmn3NvSGht8tW8K3lZQ0/PoJCbXDa32htqHbFX5FRCRSHPZ/kkaPHs3o0aMBWLt2Lbm5uZSEJIGTTz6ZrVu3MnToUL/FdfDgwSRomLk0gBkkJnpL27aNc87KSq+LQ12htuZSX1l+fvj63nR5iI9vnMAb2t9X/w4UEZF90aznUa2urubNN98kNzeX3NxcvvnmG5xz3HjjjTz99NNUVVUxdepUjjvuOFq0aNHU1RXZZ8Hwu7twu7eBeG9afoPht6Fht2XL3S+JiervKyJyuNCE/w1UUFDAzJkz6dChA0cffTRff/01Rx99NNHR0QwcOJCTTjrJ7+uaopnppZmrqqo/2O5rIC4ubti1o6J+CK0NCbahS137a6CbiEjTUVDdR8XFxcyYMYPc3Fw+//xzZs+eTXl5OR999BEjR45k6dKlLFiwgBNPPJG2jXXfV6QZC4bfoqI9Lw3Zr6ysYdeNidm7YLunIBwXd2A/JxGRw4meTLWPkpKSOO200zjttNMAKCkpYfbs2QwcOBCAiRMncvfddwPQq1cvv8X1/PPPJ14z0IvstehoSE31lsZQUbH34TZ0yc8PXw8+hGJP4uL2rXW3vnX18RWR5kotqvuhvLycuXPn+n1cZ8yYQUVFBdu2bSMuLo6JEydSUlLCsGHD6NKli+ZyFTnElZXtfytv6NLQqc1CB7jtbcitq0wtviISSXTr/yCpqqriu+++o3v37gAMHz6cqVOnAtChQweGDRvGqFGjuOSSS5qymiISAUKnNqsZcmsG3oauV1Y27NqxsQ0PtQ1ZVx9fEdkfuvV/kERHR/shFeDTTz9l8eLFfP755+Tm5jJ16lRKS0v9oDpu3Dh69uzpP4QgSs/sFGk2GntqM+e8Ft99Dbn709Uh2Me3Ia26dc3uUHNbUpKCr4h41KJ6EDnn2LVrFy1atKCoqIj+/fuzatUqANLS0jjhhBMYN24cI0eObNqKiogA5eX1h9q9CcDBnxs6uM2s/qezNTTs1nzVgyxEIpdaVCOEmfnzsbZs2ZKVK1eyZs0av49rbm4umzZtAmDJkiXcfPPN/kMIjjnmGD2EQEQOqrg4aNXKWxpDRUXdjzBu6JPc8vNh2bLwKc4a2tYSH994obdFC83lK3KwqEU1wjjnMDNmzJjBuHHj+OabbwCIj49nyJAhPPvss/Ts2bOJayki0vSqq70HT+xN2N3Ta3l5w64dFdWwQNuQn4Ovmt1Bmiu1qB5CgjMDnHDCCcyfP5+CggJ/Ltfc3Fxat24NwGOPPcZrr73mT4l1/PHHk9ZYD7AXETkEREV5j+hNToZ27RrnnMHuDvsadtetq/1Qi4YKTmvWkFDbkJ+Tk70p30QOZWpRPUS9+OKLPPvss8yePZuKigrMjGOPPZaZM2diZixcuJDU1FTat2+vabFERJpIdbX3xLWaT2Wr6+c9lQd/Li1t+PWTkvY+4O4uGKvLgxwImp7qMFZcXMysWbPIzc1l69atPPHEEwAMHTqUL774grS0NHJycujXrx8nnngiF198cRPXWERE9kdlZf2PJd6XMLw3U5sFB7qFBti61uvbVtd2hV/Zr6BqZglALhCP11XgdefcPWaWAUwEOgOrgAudc4W7O5eC6sHz5Zdf8tVXX5GXl8eCBQtYsGABJ598Mm+//TYAxx13HKmpqfTr188Psr179yYxMbGJay4iIgdbcGqzhgbcoiLYtat2Weg+DX2gRbC/774GXYXfQ9/+BlUDkp1zO80sFpgB3AqcBxQ45x40s98C6c65O3Z3LgXVphM6NVZVVRVXX301CxYsYNGiRZQF5owZO3YsTz31FJWVlTzwwAPk5OSQk5NDt27diFZHJxERaaDgAy3qa/nd1237G373pvVX4ffgabRb/2aWhBdUbwReBE52zuWbWRbwmXNut8PRFVQjT2VlJcuXLycvL48jjzySY445hhUrVtC9e3eC342EhAT69OnDvffey1lnnUVpaSlbt27liCOOUP9XERE5KGqG34aE3Ibs09DwGzq/756W+lqDay7JyV6obu72e9S/mUUDc4FuwFPOuVlmlumcywcIhNU6n61iZtcB1wF06tRpX+ovB1BMTAy9evWiV69e/rbs7Gx27tzJokWLwroOBLsFzJw5kxEjRpCenh7WdeDcc88lMzOzqd6KiIgcxkKf5tamTeOcM/RRxvUF2tBuDnUtmzfDypXhx1VVNbwOwQFv+7uEzvZwOE11trctqmnAW8DNwAznXFpIWaFzLn13x6tF9fCwdu1a3nnnHT/E5uXlsWPHDubMmcPAgQN5++23efbZZ2v1f42Pj2/qqouIiBxQznnTnNUXbPd2CQbmhj7ZDbwHXOxr4D3hBAjMhHnQNNo8qs65bWb2GTAS2GhmWSG3/jftf1XlUNCxY0duuukmf905x9q1a8nKygKgtLSU9evXM2XKFMoDs2dHR0eTn59PmzZtmDlzJhs3bqRfv3507dpV/V9FROSwYeYFxfj4xnuqG3hPdttdy25Dli1bwtd37ap9nc8/h2HDGq/e+2uPQdXM2gAVgZCaCIwAHgLeBa4AHgy8vnMgKyqRy8zCunVcdNFFXHTRRVRUVLB8+XIWLFjA0qVL/YcVPP3007zyyisAJCYm0rdvXwYOHMgzzzwDwK5du0hKSlL/VxERkYDYWEhL85bGEjrPb3Dp1q3xzt8YGjLqvz/wAhANRAGvOefuM7NWwGtAJ2AN8FPnXMHuzqVb/wJeEF24cGFY/1eAyZMnA3DKKaewYMGCsK4DgwYN4kc/+lFTVltEREQOAE34L4eUCRMm8MUXX/j9X3fu3MkZZ5zBhx9+CMCYMWNISkoiOTmZxMREEhMTOfrooxk+fDgA77zzDvHx8X5ZYmIimZmZtG3bFuccJSUlJCQkEKWhliIiIk1OQVUOWdXV1axZs4aSkhJ69+5NVVUVgwYNYsWKFZSUlFAZeJzKDTfcwF//+leqqqqIiando+X222/nkUceYfv27aQF7puEhtk77riDW2+9lc2bN3PeeeeFhdzExEQuvfRSRowYwZYtW3jmmWdITEwkISHBLz/22GPp0qULO3fuZPHixbWOb9myZZ31EhERae4abTCVyMEWFRVF586d/fXo6GjmzZvnr1dWVlJSUuL3ZzUzvvrqK0pLSykpKfGXboFON7GxsTz44INhZaWlpWRnZwP4QbeoqIhNmzb5+5x44okArF+/nrvuuqtWPSdMmECXLl345ptvOP7442uVT5o0iQsuuICpU6dywQUX1Aq6Tz75JEOGDOHLL7/k8ccf97cnJCSQnJzMmDFj6NixI2vWrGHhwoWkpqaSkpIStqiFWEREDjcKqnJIi4mJoWXLlv56VFQURx99dL37JyUlcccd9T9ArV27dkybNq3e8v79+1NeXl4r6Abnj+3VqxfvvfdeWHlJSQkDBgzwz3/ppZfWOj44dVdBQQFz584NKy8uLmb06NF07NiRf//731x//fW16rVo0SJ69+7N+PHjefTRR2sF2ccee4z09HRmzZrF/Pnza4Xcvn37Eh0djXNOg9hERCRi6Na/SISrDjw2JSoqii1btrB8+XJ27NjBjh072L59Ozt27ODqq68mNTWV9957j5deeskvDy7ffPMNaWlp/Pa3v+Whhx6qdY1gWL7lllv4+9//TkpKih9209LS+PjjjzEzXn31VRYuXBhWnpGRwWmnnQbAli1biI6OVlcHERFpMPVRFREASkpKKCgoqBVkzz//fMAbiDZ9+vSwIFxVVcUnn3wCwJVXXslLL73kh2eArKws1q9fD8BZZ53F+++/D3it1ykpKeTk5PDpp58CcM8997B69Wq/JTc1NZUuXbpwwQUXAPDNN98QFRVFeno6aWlpmqZMRKQZUFAVkUbjnGPXrl1+mC0rK+Ooo44C4KOPPmLp0qV+S++OHTtIT0/3W3EvuugivvjiC7+surqa448/nhkzZgDQt29fFi1a5F8rNjaWs88+m9dffx3wgnJJSYkfZNPT0znqqKMYOXIkAPPmzSM5OZn09HRSU1OJi4s7mB+NiIjsAwVVEYk4zjmKi4spKysjIyMDgOnTp7NhwwYKCwvZtm0bhYWFdOnSheuuuw6A008/ndWrV1NYWEhhYSEVFRVcdtllvPjii4DXiltSUuJfIykpibFjx/LII49QXV3NOeecQ2pqqh9009LSGDp0KMceeyyVlZUsWLDAL9MANRGRg0Oj/kUk4pgZycnJJCcn+9uCsyvU5+OPP/Z/Ds6JG5yizDnHpEmT/IAbfB08eDDgdXtYt24deXl5bNu2jW3btuGc45577uHYY49ly5YtYQ+VMDNSU1N54IEHGDt2LPn5+YwbNy6sNTctLY0RI0bQq1cviouLWblypb89MTFR3RZERPaTgqqIHJLMjKSkpLD1UaNG1bt/cnIyX331lb9eXV1NUVGR32qakpLCm2++GRZ0t23bRq9evQAoKipi6dKlfnlxcTEAzz//PL169eLrr78Om5osLi6OtLQ0JkyYwJlnnsnXX3/Ngw8+WCvo/uQnP+GII45gx44dbN26lbS0NFJTU9WaKyKCgqqINFNRUVGkpqb660lJSZx77rn17t+jRw/y8vL89fLycrZv3+6H5e7duzNx4sSw1txt27bRvn17AAoLC5k3b57fbSHYEpybm8sRRxzBe++9x89//nPgh9bc9PR03nvvPfr27cvUqVN59dVXSU9PD1vOPPNMWrRowY4dO6ioqCAtLY3o6OhG/7xERJqCgqqIyD6Ii4ujTZs2/nqbNm248MIL693/lFNOYcmSJcAP/XMLCwtp3bo1AMcddxwTJkzwg2xwCT5JbeXKlbz77rsUFhZSXl7un3ft2rW0aNGCxx9/nLvvvhvwWoeDQfbzzz8nJSWFt956i//85z+1gu6pp55KVFQUxcXFxMXFaVoxEYko+j+SiMhBVlf/3K5du9K1a9d6j7nmmmu45ppr/L65wVbbdu3aAXDGGWeQkpJSK+gGW3xnzZrFk08+SWlpqX/OmJgYP/SOGzeO559/nhYtWvghtmPHjv50Y//85z/9PrjBJTMzk4EDBwJeVwp1VxCRxqZR/yIizUhpaakfYouKihgyZAgAH3zwAXPnzg0LuXFxcUyaNAmA0aNH8+6774adq2fPnnz77bcAnHTSScydOzcsyA4cOJDHHnsMgH/84x/+1GIZGRmkp6eTlZVFp06dDuK7F5FIpOmpRERkv5WXl4f1wa2urmbo0KEAPPfcc3z77bdhQbdbt26MHz8e8Pr4Llu2LOx8I0eO5KOPPgKgX79+7Ny5MyzI/vjHP+bGG28E4KWXXiIpKckPwRkZGbRu3TqsVVpEDk2ankpERPZbXFwcmZmZZGZm1iobM2bMbo8NTgtWWFhIQUEBBQUFYYPZRo0axfr16ykoKKCwsJDvv/+ejh07Al6f3quvvtofgBY0duxYnnrqKSorK+nbt29YyM3IyGDUqFGMHDmS8vJyPv74Y397cJ/4+PhG+FRE5EBSUBURkQMuLi6Otm3b0rZt2zrLH3zwwd0e/9133/khN/jas2dPAMrKyhgwYACFhYVs2rSJJUuWUFhYSFZWFiNHjmTjxo2cffbZtc752GOPcdttt7F69Wouv/zysJCbnp7O6NGjycnJoaioiG+//dbfnpqaqpkVRA4SBVUREYloZkbHjh39FtaakpOTmThxYr3Ht23blv/+979hIbewsNCf97a8vBwzY8WKFX55cXExXbp0IScnh6+++oqTTz45rD6pqan861//YuTIkcydO5eHH37Yb60NBtozzzyTdu3aUVRURFFRERkZGSQkJDTqZyNyuFNQFRGRw1p8fDyDBtXZ/Q3w5sD97LPPwraVlZX5TxbLycnxpwYLDbvBgWAFBQV8/fXXfllVVRUA//nPf2jXrh1vvPEGV111FQCJiYl+mH3jjTfo3r07n3/+OR988EFY0M3IyGDo0KEkJCRQWVlJdHS0nnQmzZKCqoiISA2h/VdbtWrFWWedVe++p556atgcuUVFRRQWFvp9eY877jieeeYZv29ucGnRogUA8+fP54knnqCsrCzsvBs2bCAhIYH777+fP/zhD7WC7KRJk0hISGDKlCl+14TQJTs7u7E/FpGDTqP+RUREIkBJSUlYkD3++OOJiYlh8uTJTJkypVbXhf/+979ERUUxZswYf3aFoKSkJHbt2gV4g84++eSTsBDbqVMnv1/w559/7ndNCO26EBsbe9A/A2meND2ViIjIYaqioiJsNoWCggJKSkr46U9/CsAzzzzD9OnTw8pTU1MJ/j0eMWIEU6ZMCTtn3759/UcGX3311axZsyYsxPbp04fLLrsMgLlz5xIXF0erVq3UD1f2iYKqiIiI1Gnt2rVs2LAhLMgmJydz5ZVXAt5Ty4J9cIPLKaecwieffAJ4T1VbuXKlf77ExEQuvPBCnn/+eQCuvfZaoqKiwlps+/fvz+DBgwHIz88nPT1dAbcZ0zyqIiIiUqfdzagA8NRTT4WtO+f8R+8CPP/882zcuDEsyPbq1csvnzdvHuvXr2fr1q1UVFQAcP311zN48GCqqqpo3749zrmwgWZjxozh5ptvpqysjLvvvtvfHmy17dGjB0cccQTOOQ0yO8wpqIqIiEiDmVnYYLNhw4btdv+5c+cCXsAtLi5m69atfv/X6upq/vrXv/oBd+vWrRQUFNCyZUsACgsL+fOf/xwWjMGbd/eOO+5g5cqV9O3bNyzEZmRkMHbsWEaMGMHmzZt56623wsoyMjLIzMwkLi6uMT8WOUAUVEVEROSAMzOSk5PDHnsbGxvL9ddfX+8x7dq1o7S0lJKSEj/EFhQU0LlzZ8AbNHbzzTeHlS1btozt27cDsGTJkjrPP3HiRC688EJmzpzJDTfcEBZiW7VqxZgxY+jevTv5+fksXrw4LAgnJSWpFfcg2mNQNbOOwItAO6AaeNY597iZ3QuMATYHdr3TOffhgaqoiIiIND9mRlJSEklJSbW6KLRr146HH3643mOHDBnC2rVrw7olbN261e8fGx8fT/fu3SkoKGDFihXMnj2brVu3cvbZZ9O9e3cmT57M5ZdfHnbO+Ph4vvzySwYMGMAHH3zA+PHj/RAbfP3Zz35GSkqKP7AtIyODxMTExv9wmoGGtKhWAr9yzn1lZi2BuWb2aaDs/znnHj1w1RMRERHZN7GxsXTo0IEOHTrUWT5o0CDefPPNWtuDA81HjhzJtGnTwrolFBQUcMQRRwCwfft2li9f7gfc4Fy4o0aNIiUlhSeffJJ77rkHCH/Yw8yZM2nZsiVvvfUWs2bNCgu5rVq14sQTT8TMqK6uJioq6kB8NIeMPQZV51w+kB/4ucjMFgPtD3TFRERERJpC8NZ+mzZtwh6fW9Mll1zCJZdcAnjhNjgXbvBhD2effTbt2rULC7lbt24lKSkJ8J5e9vjjj/uDzADi4uIoLS0FvKnBJk2aFBZiO3bs6M+o8P7777Np06awPritW7f2r3842KvpqcysM5AL5AC/BK4EdgBz8FpdC+s45jrgOoBOnToNXL169f7WWUREROSw4Jxj165dfogtKiryB6hNmjSJL7/8MqxFNzExkU8/9W5sn3baaf7PQb169WLx4sUAXHjhhSxbtiysxbZfv36MGzcOgNzcXKKiosKCblM86KFR5lE1sxbA58ADzrk3zSwT2AI44H4gyzl39e7OoXlURURERBrHzp072bJlS1gf3NjYWM4991wAfv/73/PNN9/4IXfr1q0cffTRfPzxxwB0796d5cuX++f7/e9/z/3333/Q38d+B1UziwXeBz52zj1WR3ln4H3nXM7uzqOgKiIiIhIZ5s+fHzYH7sCBAxkyZMhBr8d+TfhvXkeNvwOLQ0OqmWUF+q8CnAvkNUZlRUREROTAO+qoo5q6CnvUkFH/xwOXAQvM7OvAtjuBi81sAN6t/1VA/ROhiYiIiIjspYaM+p8B1DWzreZMFREREZEDpnlPziUiIiIiEUtBVUREREQikoKqiIiIiEQkBVURERERiUgKqiIiIiISkRRURURERCQiKaiKiIiISERSUBURERGRiKSgKiIiIiIRSUFVRERERCKSgqqIiIiIRCQFVRERERGJSAqqIiIiIhKRFFRFREREJCIpqIqIiIhIRFJQFREREZGIpKAqIiIiIhFJQVVEREREIpKCqoiIiIhEJAVVEREREYlICqoiIiIiEpEUVEVEREQkIimoioiIiEhEUlAVERERkYi0x6BqZh3NbJqZLTazhWZ2a2B7hpl9ambLAq/pB766IiIiItJcNKRFtRL4lXOuN3AsMM7M+gC/BaY457oDUwLrIiIiIiKNYo9B1TmX75z7KvBzEbAYaA+MBl4I7PYCcM6BqqSIiIiIND971UfVzDoDRwOzgEznXD54YRZoW88x15nZHDObs3nz5v2rrYiIiIg0Gw0OqmbWAngD+IVzbkdDj3POPeucG+ScG9SmTZt9qaOIiIiINEMNCqpmFosXUl92zr0Z2LzRzLIC5VnApgNTRRERERFpjhoy6t+AvwOLnXOPhRS9C1wR+PkK4J3Gr56IiIiINFcxDdjneOAyYIGZfR3YdifwIPCamV0DrAF+emCqKCIiIiLN0R6DqnNuBmD1FA9v3OqIiIiIiHj0ZCoRERERiUgKqiIiIiISkRRURURERCQiKaiKiIiISERSUBURERGRiKSgKiIiIiIRSUFVRERERCKSgqqIiIiIRCQFVRERERGJSAqqIiIiIhKRFFRFREREJCIpqIqIiIhIRFJQFREREZGIpKAqIiIiIhFJQVVEREREIpKCqoiIiIhEJAVVEREREYlICqoiIiIiEpEUVEVEREQkIimoioiIiEhEUlAVERERkYikoCoiIiIiEUlBVUREREQikoKqiIiIiEQkBVURERERiUh7DKpm9g8z22RmeSHb7jWz783s68By5oGtpoiIiIg0Nw1pUX0eGFnH9v/nnBsQWD5s3GqJiIiISHO3x6DqnMsFCg5CXUREREREfPvTR/UmM/sm0DUgvb6dzOw6M5tjZnM2b968H5cTERERkeZkX4PqX4FsYACQD/ypvh2dc8865wY55wa1adNmHy8nIiIiIs3NPgVV59xG51yVc64aeA4Y3LjVEhEREZHmbp+CqpllhayeC+TVt6+IiIiIyL6I2dMOZvYv4GSgtZmtA+4BTjazAYADVgHXH8A6ioiIiEgztMeg6py7uI7Nfz8AdRERERER8enJVCIiIiISkRRURURERCQiKaiKiIiISERSUBURERGRiKSgKiIiIiIRSUFVRERERCKSgqqIiIiIRCQFVRERERGJSAqqIiIiIhKRFFRFREREJCIpqIqIiIhIRFJQFREREZGIpKAqIiIiIhFJQVVEREREIpKCqoiIiIhEJAVVEREREYlICqoiIiIiEpEUVEVEREQkIimoioiIiEhEUlAVERERkYikoCoiIiIiEUlBVUREREQikoKqiIiIiEQkBVURERERiUgKqiIiIiISkfYYVM3sH2a2yczyQrZlmNmnZrYs8Jp+YKspIiIiIs1NQ1pUnwdG1tj2W2CKc647MCWwLiIiIiLSaPYYVJ1zuUBBjc2jgRcCP78AnNPI9RIRERGRZm5f+6hmOufyAQKvbevb0cyuM7M5ZjZn8+bN+3g5EREREWluDvhgKufcs865Qc65QW3atDnQlxMRERGRw8S+BtWNZpYFEHjd1HhVEhERERHZ96D6LnBF4OcrgHcapzoiIiIiIp6GTE/1L+ALoKeZrTOza4AHgVPNbBlwamBdRERERKTRxOxpB+fcxfUUDW/kuoiIiIiI+PRkKhERERGJSAqqIiIiIhKRFFRFREREJCIpqIqIiIhIRFJQFREREZGIpKAqIiIiIhFJQVVEREREIpKCqoiIiIhEJAVVEREREYlICqoiIiIiEpEUVEVEREQkIimoioiIiEhEUlAVERERkYikoCoiIiIiEUlBVUREREQikoKqiIiIiEQkBVURERERiUgKqiIiIiISkRRURURERCQiKaiKiIiISERSUBURERGRiKSgKiIiIiIRSUFVRERE5P+zd+fRVVX3+8efneRmnicEAgQEZAgCEhlERVQULRVrHdCqVKUoitr2W6tinavWZl+1AgAAIABJREFUijOyFFuK1CIOrQoqg+KIMv9ESAQUZRASICQxCWRO9u+Pm3vIJQQCCeSQvF9rnXXvPZ8z7OQu4uM+++wDVyKoAgAAwJUIqgAAAHCloMbsbIzZLKlIUpWkSmttelM0CgAAAGhUUK0x3Fq7uwmOAwAAADi49A8AAABXamxQtZIWGmNWGWPGH2gDY8x4Y8xKY8zKnJycRp4OAAAArUVjg+pQa+0pki6QdIsx5sz9N7DWTrPWpltr05OSkhp5OgAAALQWjQqq1tqsmtddkt6WNLApGgUAAAAccVA1xkQYY6J87yWdJymjqRoGAACA1q0xd/23kfS2McZ3nFnW2vlN0ioAAAC0ekccVK21P0rq24RtAQAAABxMTwUAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXalRQNcaMNMZsMMZsNMbc1VSNAgAAAI44qBpjAiW9IOkCSb0kXWmM6dVUDQMAAEDrFtSIfQdK2mit/VGSjDGzJY2W9G19O2zYsEFnnXVWI04JAACA1qIxl/7bS/qp1udtNev8GGPGG2NWGmNWVlRUNOJ0AAAAaE0a06NqDrDO1llh7TRJ0yQpPT3dfvrpp404JQAAAFoSYw4UKb0a06O6TVKHWp9TJGU14ngAAACAozFBdYWkbsaYzsaYYEljJM1pmmYBAACgtTviS//W2kpjzERJCyQFSppurc1sspYBAACgVWvMGFVZaz+Q9EETtQUAAABw8GQqAAAAuBJBFQAAAK5EUAUAAIArEVQBAADgSgRVAAAAuBJBFQAAAK5EUAUAAIArGWvtsTuZMTmSthyzE3olStp9jM+Jo4/vtWXie22Z+F5bJr7Xlqe5vtNO1tqkAxWOaVBtDsaYldba9OZuB5oW32vLxPfaMvG9tkx8ry2PG79TLv0DAADAlQiqAAAAcKXWEFSnNXcDcFTwvbZMfK8tE99ry8T32vK47jtt8WNUAQAAcHxqDT2qAAAAOA4RVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBYBjwBiz2RhzbnO3AwCOJwRVAHAxY8yrxphsY0yhMeY7Y8y4WrV4Y8zbxpi9xpgtxpir9tv3oHUAcLug5m4AAOCgHpN0g7W2zBjTQ9KnxpivrbWrJL0gqVxSG0n9JL1vjPnGWptZs++h6gDgavSoAmi1jDF3GmO2G2OKjDEbjDHnGGOuMMbsqbWUGWM+beDxzCE26WeMWWOMKTDGvG6MCa2vHb4drLWZ1toy38ea5URjTISkX0u611q7x1q7WNIcSdfUHPOgdQA4HhBUAbRKxpiTJE2UdKq1NkrS+ZI2W2tft9ZGWmsjJbWT9KOk1xpwvDMkfeALn/W4XNJISZ0lnSzpt/W1Y79jTzXGFEtaLylb0geSukuqstZ+V2vTbyT1rnl/qDoAuB5BFUBrVSUpRFIvY4zHWrvZWvuDr2iMCZA0S9Kn1tqXGnC8LyXtkjTnIGH1OWttlrU2T9JceS/HH7QdkmStvVlSlKQzJP1PUpmkSEkF+x2/oGY7NaAOAK5HUAXQKllrN0r6vaQHJO0yxsw2xrSrtckj8oa622rvZ4wZaYyx+y/yBs5rJY2QNKGe0+6o9b5YUmQD2uFrb1XN5fuUmuPvkRS932bRkopq3h+qDgCuR1AF0GpZa2dZa0+X1EnesZ+PS5IxZoykKyVdaq2t2G+f+dZas/8iKVDSTEkfSnqxKdpRjyBJJ0r6TlKQMaZbrVpfSb4bpQ5VBwDXI6gCaJWMMScZY842xoRIKpVUIqnKGNNf0vOSLrbW5hzGIYfKe3f9aGttSWPbUVNLNsaMMcZEGmMCjTHnyxugP7bW7pV3GMBDxpgIY8xQSaMl/VuSDlUHgOMBQRVAaxUi6W+Sdst7ST5Z0iR5w1ycpMW17vyfd6iDWWu/kHTB4YTUQ7RD8vauTpC0TVK+pMmSfm+tfbemfrOkMHnHxr4macJ+U08dqg4Armastc3dBgAAAKAOelQBAADgSgRVAAAAuBJBFQAAAK5EUAUAAIArBR3LkyUmJtrU1NRjeUoAAAC42KpVq3Zba5MOVDumQTU1NVUrV648lqcEAACAixljttRX49I/AAAAXImgCgAAAFciqAIAAMCVCKoAAABwpQbdTGWM2SypSFKVpEprbboxJl7S65JSJW2WdLm1Nv/oNBMAAACtzeH0qA631vaz1qbXfL5L0iJrbTdJi2o+AwAAAE2iMZf+R0t6peb9K5IubnxzAAAAAK+GBlUraaExZpUxZnzNujbW2mxJqnlNPtCOxpjxxpiVxpiVOTk5jW8xAAAAWoWGTvg/1FqbZYxJlvShMWZ9Q09grZ0maZokpaen2yNoIwAAAFqhBvWoWmuzal53SXpb0kBJO40xbSWp5nXX0WokAAAAWp9D9qgaYyIkBVhri2renyfpIUlzJI2V9Lea13ePZkOP1Nq1a7V9+3YFBgYqKChIQUFBCg4O1qBBgyRJP/30k4qKiurU27ZtK0nau3evrLUKCgpytjHGNOePBAAA0Co05NJ/G0lv14SzIEmzrLXzjTErJL1hjLlB0lZJlx29Zh65KVOmaNq0aX7rIiIitGfPHknSnXfeqddee82vfsIJJyg7O1uSNGbMGL333nt+9e7du2vDhg2SpFGjRumrr75yQm5gYKD69u3r7DNmzBitX7/eqQcFBemUU07Rc889J0m68cYblZ2d7ReEBwwYoD/96U+SpLvuukt79uzxC9L9+/fXmDFjJElPPvmkqqqqFBQUpIiICPXs2VN9+vRRXFxcU/4aAQAAjrlDBlVr7Y+S+h5gfa6kc45Go5rSXXfdpeuuu06VlZWqqqpSZWWlX/3222/XxRdfrMrKSmeb4OBgp37DDTdo2LBhTr2ystIvBJ5//vnq0qWLX71Dhw5OPSUlRaWlpX7HDwra92vftWuXfvrpJ6dtVVVVio6Odurvv/++srOz/fa//PLLnaB63333qbi42O9nGj9+vF566SVVV1frz3/+s3r16qW0tDT16tVLkZGRTfOLBQAAOMqMtcfu/qb09HS7cuXKY3a+1qCsrMwJuT///LO+/fZbJScn65RTTtG2bdvUvXt3lZSUONunpqbq0Ucf1ZVXXqni4mJ9//336tGjh0JCQprxpwAAAK2VMWZVrXn6/TT0rn+4VO2AGR0drY4dOzqfU1JSVFRUpM2bNysjI8NZkpO9M4mtXLlSw4YNU2BgoLp166a0tDSlpaXp2muvVefOnY/5zwIAAFAbPaqt2O7du/Xxxx/7hdiNGzdqyZIlGjRokN566y09/PDDToD1LZ06dVJAQGOeFQEAAOBFjyoOKDExUZdffrkuv/xyZ11xcbEzRjcyMlIpKSlavHixZs2a5Wzz008/KSUlRQsWLFBmZqYTYNu2bcuMCAAAoMkQVOEnPDzceT9y5EiNHDlSklRQUKBvv/1W3377rdq3by9JmjNnjqZOnepsHxcXp379+mnRokUyxujHH39UTEyMEhISju0PAQAAWgQu/aNRdu3apczMTGVkZCgzM1N79uzRq6++KkkaMWKEPvroI7Vt21a9e/dWWlqaBg8erCuuuKKZWw0AANziYJf+Cao4aj777DOtXLnSGf+amZmp008/XQsXLpTkndorODjYGTrQu3dv9ejRQ6Ghoc3ccgAAcKwwRhXNYtiwYRo2bJjzubq6WgUFBZIka63atGmj1atXa8GCBaqoqJAkjR07VjNmzJC1Vn/961/Vo0cPpaWlqWvXrvJ4PM3ycwAAgOZBjyqaXUVFhb7//ntlZGSoffv2Gjp0qLKzs5WSkqLq6mpJUnBwsHr06KFJkybpiiuuUGFhoZYtW6aYmBjFxMQoOjpaMTExCgsL44YuAACOI/SowtU8Ho969eqlXr16Oevatm2rPXv2aP369X7TZ/lu9vr222913nnn1TnWa6+9pjFjxmjZsmWaMGFCnSB78803q2fPntq6dau++uorv1pMTIxOOOEEem4BAHAJgipcKywsTP3791f//v3r1Hr37q3PP/9cBQUFKigoUGFhoQoKCpxtg4KC1L59exUUFGjTpk1O/bLLLlPPnj315Zdf6qqrrqpz3KVLl2rQoEH6z3/+o0mTJtUJsk888YTat2+vVatWaenSpXXqvXv3lsfjkbWWnl0AABqJoIrjUlRUlM4444x66wMGDNDcuXPrrY8aNUqZmZl1gm7Xrl0lSe3atdPw4cOd2s6dO/Xdd9/JN1Rm/vz5+stf/lLnuLt27VJSUpLuvfdePfPMM3WC7Ny5cxUcHKw5c+Zo9erVfrXY2FidffbZkqQ9e/bI4/HwaFsAQKvGGFXgCJSXlys/P98v5BYUFOiiiy5SUFCQ5s+frw8//NCvXlRUpC+//FLGGE2YMEEvvvii3zEjIyNVVFQkSbrqqqv02muvKTg4WNHR0QoLC1OnTp30xRdfSJLuv/9+ffPNNwoPD3eWDh066I477pAkvfvuu8rLy/OrJyYmqm/fvpK8gTooKEjh4eEKCQmh9xcA0GyYngpwoYqKChUWFjpBtqSkREOGDJEkvf/++1qzZo0TcIuLixUZGannn39eknTbbbfp888/V3FxsbN06dJFvn9fQ4YM0dKlS/3ON3jwYC1ZskSS1KdPH2VkZEiSjDEKDw/XyJEj9dZbb0mSRo8erfz8fL+gO2TIEN1yyy2SpKeeekrWWr96165dnaEXa9euVVhYmF/d4/EQiAEAdbTaoLpypfTEE9Lgwd6lf3+JKTrRGuTl5TkB17eEhobq1FNPlSS9/vrr2rFjh1+9a9eumjBhgiTp6quv1vbt21VSUuLUzzvvPOdJZFFRUdqzZ4/fOceNG6eXX35Z1loFBgZq/78tv//97/X000+rpKREvXv3Vnh4uF+Yveaaa3TVVVepsLBQ9957r1P3LWeccYb69u2rvXv3avHixX618PBwJScnKzIy0jkvoRgAjg+t9q7/rCxp6VLpjTe8nz0eb1j1BdfBg6XUVIn/nqGliY+PV3x8fL31Qz0dzPd0sfrk5+f7hdji4mJFRUVJ8s6R++abbzrrfdulp3v/BlVXV+v000/32zc3N9cJvoWFhXrllVdUUlKi8vJy55xPP/20+vbtqy1btjiP9q3t5Zdf1rhx47RixQoNGTLEL8iGhYXpmWee0QUXXKDVq1frL3/5i1/IDQsL0/jx49WzZ0/98MMPWrBgQZ36qaeeqtjYWBUWFiovL8+vHhTUov+UAkCzadF/XS+6yLtkZUnLlnlD69Kl0ssvS889590mOdk/uKanSzX/vQVQj6CgIEVFRTnhtLaAgAD9+te/rnffiIgIzZw5s956SkqKfv75Z0lSVVWVSkpKVFJSorCwMElSamqqvvzySycA++qnnXaaJOmEE07Q3Xff7az3bRcXFydJKikpUXZ2dp39R40apZ49e2rVqlXOEIfavvrqKw0ZMkT//e9/df3119f5fXz99ddKS0vTzJkz9dhjj9XpEf7HP/6h5ORkLViwQAsXLlRYWJhCQ0MVGhqqsLAw3XDDDQoNDVVmZqY2b97sVwsNDVXPnj1ljFFxcbEkKTQ0VAEBAQf7mgDguNeiL/3Xp7JSWrt2X3BdulT67jtvLSBA6t3bP7z26OFdD6Dl890ot3/Q7du3r6KiorRx40YtXry4TtD9/e9/r6SkJL3//vuaMWNGnf0XLlyo5ORkPf7443r44YdVUlLiPNBC8vYkR0VF6U9/+pOefPLJOu2qqqpSQECAbrrpJr300kuSvA/CCA0NVXx8vDZt2iRJuueee7Ro0SK/oNumTRvn5r1//vOf2rhxo18ITk5O1pgxYyRJS5Ys0Z49e5z9Q0NDFRMTo44dO0ryBn2Px0MvMoAm02rHqB6OvDxp+fJ9wXXZMqmmU0cxMdLAgfuC66BBUkJC87YXwPGvoqJCpaWlKi0tVWJioowx2r59u7KyslRSUuLUSktLdfnll0uSFi5cqK+//lqlpaXONgEBAXrqqackSY8//rg++eQTv/3j4uL06aefSpJ+9atf6f3333ceWyxJvXr1UmZmpiTp9NNP15dffunXzlNPPVXLly+XJPXr10/ffPONAgMDnaB79tln6/XXX5ckXXLJJcrJyfHrMR4yZIj+8Ic/SJIeeOABlZeXKyQkRKGhoQoJCVFaWprOPfdcSdI777yjoKAghYSEOEu7du3UoUMHWWuVnZ3tVwsKCmI8MnCcI6gegepqby9r7V7XtWu96yWpWzf/Xtc+fbxjYAHgeFBVVaWysjKVlJSoqqpKycnJkqR169YpLy/PL+hGRUXp/PPPlyRNnz5d27dvd2olJSXq1q2bE0SvvfZabdu2zamVlJTo3HPP1ZQpUyR5h3bs3LlTlZWVTlt++9vf6l//+pck75PqatckaeLEiXr++edVVlam0P3uiA0ICNA999yjhx56SPn5+erXr59fCA4JCdGECRP0m9/8Rjk5Obrtttv8gm5ISIguueQSDRkyRLt379Zrr73mrPcdIz09XSkpKSoqKtKGDRv8jh0SEqK4uDgFBwcfnS8KaAVa7c1UjREQ4L3k36OH9Nvfetft2eOdScAXXBcskP79b28tLMw7vrV2eG3XrtmaDwAHFRgY6My4UFvPnj0Put/+43P3d7Dxx5K0bds2SfuCcllZmQIDA536qlWrnPW+pUOHDpK8ofTFF1/0q5WWlur000936sOHD3fW+7bxjeUtLi4+4PG7deumIUOGaMuWLbrtttvqtPnf//63rr76an399dcaNmxYnfrbb7+tiy++WB988IFGjx6t4OBghYSEKDg4WMHBwZo9e7ZOO+00ffjhh7rvvvuc9b5t/v73v6tr165avHixM3+yrxYSEqIbb7xRiYmJWrt2rVasWOF37ODgYA0fPlyhoaHKyspSTk5Onf0TEhIUEBDAE/NwXCKoHobISOmss7yLJFkrbdni3+v6zDOS74pahw7+wfWUU5geCwCk+oPyySefXO8+Ho9HN954Y731mJgYzZgxo956p06d9J3vhoQD6Nu3r3bv3u0XgsvKypzxub1799acOXPqBN1+/fpJkjp37qw///nPKisrU3l5ufPqm4EjKChIMTExKisrU3FxsfLz81VeXu7MbvHDDz/ojTfe8Ntf8s7SkZiYqPnz5+vPf/5znXZnZWWpbdu2evHFF/Xwww/XqRcVFSkyMlL/93//p+eee65OkN66dauMMfrrX/+quXPn+tViYmL02muvSZKmTZum1atX++0fHx/v9Ka/99572r59u1+IjouLc4Z1rFmzRnv37lVwcLA8Ho+Cg4MVGRmplJQUp50BAQEKDg5mSAccXPpvYqWl0urV/mNdN2/21jweqV8///DauTPTYwEA6rLWqrKyUoGBgQoICFBRUZHy8vKccOsLtAMGDJDH49H69eu1bt26OkH5xhtvVFBQkObNm6cvv/zSWV9eXq7Kykq9/PLLkqRnn31W8+bN8zt+aGio80S88ePH6+233/Y7dkpKirZu3SpJGjlypBYsWOD3M9Qe/zx06FB99dVXfvWBAwdq2bJlkrz/o7BmzRqn5vF4dN555+m9996T5B0/nZWV5YRgj8ejc889V48//rgk7xP9iouL/YLy6aef7lwFuPfee2WM8QvK/fv311lnnaXq6mrNnj3br+bxeNSlSxd16dJFlZWV+vbbb/1qHo9HMTExCg8Pl7VWVVVVCgwMJGAfAcaoNrMdO/ynx1q+XKqZYUZJSftu0Bo8WDr1VCk6unnbCwBAQ/jCmST9/PPPKi4uVkVFhRN2AwMD1aNHD0nSypUrlZub69QqKioUGxvrzIs8c+ZM7dy506mXl5erS5cu+t3vfidJuuOOO7Rjxw6/4w8ePFh/+ctfJElnn322cnNz/eq//vWv9fTTT0vyTo3nm97N59Zbb9Vzzz2n0tJSZwq82u6++249+uij2r17t5KSkurUH3nkEU2aNEmbNm1Sly5dJMkJsR6PR5MnT9a4ceO0YcMG/fKXv/SreTwe3Xvvvbrgggu0bt063XXXXc56Xxi++eabdcopp2jDhg3617/+VScoX3rpperUqZM2b96szz//vM7xTzvtNMXGxionJ0dbtmzxqwUHB6tdu3byeDzO9+G7QfFYY4xqMzvhBGn0aO8ieafHysz0HzIwd663Zkzd6bF69mR6LACA+9QeXxwbG6vY2Nh6t/U99KM+11577UHrTzzxxEHrH3/88UHre/fudXo+fUHWF8qCg4O1fv16J+T6Xtu3by9JioyM1FtvvaWKigq/bQYOHCjJO+zkoYcecuq+5aSTTpLknfd4wIABdeq+85eUlGjLli1+tfLycl166aWSpE2bNunpp59WRUWF31P/+vXrp06dOmnp0qUaO3ZsnZ952bJlGjhwoN59910n8NeWmZmpXr16aerUqfrDH/6g+++/Xw888MBBf4/HGj2qLpGf7z891tKl+6bHio6uOz1WYmLzthcAABxbvqDtC7NhYWHyeDzau3ev0xtdO+z27t1bkZGR2rp1q7755ps6QfmSSy5RTEyMVq5cqY8//lhDhw7V0KFDj/nPxaX/41B1tfT99/7Bdc2afdNjde3q3+t68slMjwUAAI4/BNUWYu9eadWqfcF1yRLv+FfJO5vAKad453Pt00dKS/MuPJgAAAC4GUG1hbJW+umnfcF1xQopI2PfkAHJOz7WF1p9S+/e3qm2AAAAmhs3U7VQxkgdO3qXmqcrylopO9v7FK2MjH3LSy9JJSX79u3cuW6APekkKSSkeX4WAACA/RFUWxhjvE/EatdOqnnioSSpqso7n6svuPqC7Lx53lkIJCkwUOre3Rtaaw8f6NLFWwMAADiWWnxQfeyxx7R+/Xrdeuuth5waoyULDJROPNG7+KbJkqTycum77/x7X//f/5PeesvbOyt5x7/26uXf+9qnj9S+PQ8rAAAAR0+LD6qlpaX673//q5kzZ2rQoEGaOHGiLrvsMoVwjVuSFBy8L3zWtnevtG6d/xCCjz6Saj/GOyam7vCBtDSmzgIAAE2jVdxMVVBQoJkzZ2rKlCn67rvvdO211+qVV1455u1oCfLyvA8rqN0Du3atdx5YnzZt6g4f6NVLiopqvnYDAAB34q7/GtXV1Vq0aJESExPVv39/ff/995o0aZImTpyoM888k+fzHiHfDVy1w2tGhjfQ1n5aXWpq3d7XHj24gQsAgNaMoFqPuXPnauzYscrPz1daWpomTpyoq6++WhEREc3dtBahutr/Bi5f7+v69Qe+gav2cuKJ3MAFAEBrQFA9iOLiYs2ePVvPP/+8Vq9erRNOOEGbN29mDOtRVF7uferW/j2wP/zgfwNXz551hxCkpHADFwAALQlBtQGstVqyZIkyMjI0fvx4SdKf/vQnnXPOOTr//PMVEBDQzC1s+YqL697AlZEhbd++b5voaG9g7drV+9CCyEgpIqLuUt/6iAgpqMXfQggAwPGDoHoEdu3apb59+2rHjh068cQTdcstt+i6665TbGxsczet1cnPr3sD16ZN0p493tkJysoO73ghIQ0PtgcLvAeqMVwBAIDDQ1A9QuXl5Xr77bc1ZcoULV68WOHh4Zo3b57OPPPM5m4aaqms9PbG7t3rXXwBdv/lSNaXlx9eW0JCDj/0Hmp9ZKQUHi7RqQ8AaIl4hOoRCg4O1hVXXKErrrhCq1ev1rRp03TKKadIkubMmaPy8nKNHj1aHo+nmVvaugUFeYcEREc3/bErKxsWbA8VgnfsqLu+ouLw2uILrbWXqKi66w62vnYtIoLwCwBwN3pUj9CFF16oefPmqX379rrppps0fvx4JScnN3ezcBypqGh4CK5vKSqq+9k3o0JDhIcfXrg91HqGPwAADheX/o+CqqoqzZs3T1OmTNGCBQsUHBys+++/X5MmTWrupqGVKy+vG2APFm4Ptf5wh0CEhTU++EZF7VsiIpjpAQBaMi79HwWBgYEaNWqURo0ape+++05Tp05Vz549JUk5OTn64IMPdMUVVyg0NLSZW4rWJjhYSkjwLk2lvPzww23tWkGBd/aG2usbGn6NqRteay8Hqx1oCQ0l+ALA8aLBParGmEBJKyVtt9aOMsbES3pdUqqkzZIut9bm13+EltWjejBTp07VLbfcooSEBP3ud7/TTTfdpE6dOjV3swBXKS+vO7ShqKjhy/7bl5Y27LyBgYcXbg+1LVMuA0DjNMmlf2PMHyWlS4quCap/l5Rnrf2bMeYuSXHW2jsPdozWElSttfr000/1/PPP691335UkXXzxxXr99dcVxCSewFFRUXHwsHu4QbihN7t5PEfew1t7W9/74OCj+3sCALdp9KV/Y0yKpF9IekTSH2tWj5Z0Vs37VyR9KumgQbW1MMZo+PDhGj58uLZu3aoXX3xR27dvd0Lq+++/rzPPPFNRUVHN3FKg5fB4pLg479IUysqOvJe3sNA71KF2raqqYecNDj5wgD1YuD3Y+5BvDvnbAAAf4klEQVQQhjoAOH41qEfVGPOWpMckRUn6U02P6s/W2tha2+Rba+v8J8IYM17SeEnq2LHjgC1btjRZ449H27dvV8eOHRUREaHf/va3uuWWW3TSSSc1d7MAHEXWeocm1NfDe7jvD2d2h6CgA9+odqRhODyc4AugaTXq0r8xZpSkC621NxtjztJhBtXaWsul/4Ox1mr58uWaMmWKXn/9dVVUVGjEiBF69tlnnZuxAOBQysqaJvT63jf0CW8BAf6zNBxuL+/+r0xpBqCxl/6HSrrIGHOhpFBJ0caYVyXtNMa0tdZmG2PaStrVdE1uuYwxGjRokAYNGqTJkyfr5Zdf1vTp0xVXc71yw4YNSkxMVEJT3rINoMUJCfEuiYlNc7wDjfE9nKC7dav/+pKShp/bN5/vwQLt4bwy3AFoOQ5rHtX9elSfkJRb62aqeGvtnw+2Pz2qB2atlan5qzp8+HAtXbpUV111lSZOnKj+/fs3c+sA4PBVVflPR9YUr9XVDTv3gYY7NOY1MpKnuAFHU5NN+L9fUE2Q9IakjpK2SrrMWpt3sP0Jqoe2du1avfDCC/r3v/+t4uJiDR06VPfcc48uuOCC5m4aADSb2uN8myr4NnRKM8nb61s7uB5OyD3Qe+bzBfbhyVTHofz8fM2YMUMvvPCCbr31Vt1+++0qKSlRfn6+2rVr19zNA4DjXmXlgR9ScSx6fX3z+db3xLYjec/UZjheEVSPY9XV1aqsrFRwcLCmT5+uG2+8UZdeeqkmTpyo0047zRkyAABoXr5e3/2fzNaY93v2NPz8Hk/jA+/+n5n6G8cCQbWF+PHHH/XCCy/on//8pwoKCtSvXz9NnDhR1113nQIYQAUALU51tVRc3DSB1/f5cG50Cw1tWLitb93+6yMiCL+oi6Dawuzdu1f/+c9/NGXKFAUHB2vFihUyxignJ0eJiYn0sgIA6uW70a0xPb5FRf6PQG7o9GZS/eG3IaG3vhrh9/hGUG2hrLXKzc1VYmKi8vPzlZKSIo/Ho7S0NGc599xz1aNHj+ZuKgCgBfNNb3agZf+e3UOtP9y5fSXvlGSHG24PFYg9nqP3+4K/Rj9CFe5kjFFirUkUJ0+erIyMDGVkZOiNN97QSy+9pGeeeUY9evTQ1q1bNW7cOL8Q27t3b0VERDTjTwAAaAma+hHGkjf81u61bUi43X/dzp3+6w9npoeQkPpD7JEuzPF7+AiqLURcXJwmTJjgfLbWKjs7WyEhIZKkvLw85eXl6cUXX1RJrQFKc+fO1ahRo7Rp0yYtW7ZMffr0Uffu3eXhfyUBAM3I45FiY71LU6ms9A+/Dent3X+YQ06O/7bFxQ0//4FmeziSJSLC/31Lvk2FoNpCGWP8prHq16+fVq5cqaqqKm3atElr165VRkaG+vbtK0maP3++br75ZkmSx+PRSSedpD59+mjy5Mlq166dSkpKFBISwk1bAIDjVlCQFBPjXZpKVdW+G96OdMnO9r76AvHhTHUm7Xu6W1Msbdq4a6ozxqhCklRWVqb169c7QwcyMjK0du1arVmzRtHR0br33nv19NNPq3fv3kpLS1OfPn2Ulpam4cOHK5AHdQMA0GSs9Y7RbUz4PdDSkHG/H38sDR9+9H/G2hijikMKCQlR3759nR7W/Z155pkqLCxURkaG5s6dq+nTpysqKkoFBQWSvONjt27d6jcGNjo6+lj+CAAAtAjGeGdHCA2Vat2K0mgHGve7/9KzZ9OdrykQVNEgI0aM0IgRI5zPu3bt0pYtW5ypsNatW6c33nhDe2rNTn322Wdr0aJFkqSPPvpIiYmJ6tGjh0JDQ49t4wEAwFEZ93u0cekfTaa6ulpbt251xr+Gh4fr9ttvlyS1b99eWVlZCgwMVLdu3ZSWlqaLLrpI11xzjbMv418BAGh9uPSPYyIgIECpqalKTU3VL3/5S7/ahx9+6Ix7zcjI0OrVq5WSkqJrrrlG5eXlSkpKUteuXf3Gvw4YMEBJSUnN9NMAAIDmRo8qmo21VsYYFRQU6KGHHnJu4srKypIk/f3vf9cdd9yhnTt36sEHH/Qb/xofH9/MrQcAAE2BHlW4km98a0xMjJ588klnfW5urjIzM9WxY0dJ0qZNmzRr1iznxi1Jio2N1Ztvvqlzzz1Xn3zyif74xz8qNDRUISEhCg0NVWhoqB599FH16tVLy5cv16uvvurUfK/XXnutkpOT9f333+vrr7/2q4WEhKhfv34KDQ1VYWGh9u7d61djpgMAAI4+gipcJyEhQWeeeabzefDgwcrPz9f27dudXtfNmzc788SGhISoQ4cOKi0tVVlZmfLz81VWVqaKigpJ0g8//KBXX31VpaWlKi0tle8qwsiRI5WcnKx58+Y5Y2lr27Rpk1JTU/XCCy9o0qRJfrWgoCBlZ2crMTFRjz/+uKZNm+YXkkNCQrRgwQIFBwdrxowZ+uSTT/yCcEREhO6//35J0scff6wff/zRLwhHRUVpeM38INu2bVNVVZWio6MVHR1NSAYAtBpc+kerYq1VZWWlSktLFR4ersDAQOXn5ysrK0tlZWVOmC0rK9NZZ52lsLAwff3111q+fLmz3rfNfffdp9DQUL3xxhuaM2eOU/O9fvrppwoICNB9992nmTNn+tWCgoJUVFQkSfrNb36jWbNm+bWzTZs22rFjhyRp9OjRmjNnjlOLiIhQnz59tGTJEknS3XffrY0bNzpBNjo6WieeeKKuvfZaSdLy5ctVXV2t6OhoxcTEKDo6WhEREdy8BgBwhYNd+ieoAs3ANz5XkgoKClRUVOQXkq21Sk/3/pv97LPPtHHjRhUWFqqwsFAFBQWKjo7WAw88IEkaN26clixZooKCAhUWFqqoqEiDBw92guzJJ5+stWvX+p3/nHPO0UcffSRJuvDCC51j+sLs4MGDdf3110uSZs+eLY/H44TcmJgYJSUlMU4YANAkCKpAK1JVVaXS0lJFRERIklatWqWcnBwnyBYUFKhdu3a66qqrJEnXX3+9tm7d6heEL7zwQv3zn/+UJEVGRmrv3r1+5xg3bpxefvllWWvVpk0bRURE+PXYXnbZZRo7dqzKy8s1efJkv1pMTIy6du2qlJQUVVdXq6KiQiEhIcf2lwQAcA1upgJakcDAQCekStKAAQMOuv306dMPWl+7dq0TYH1hNjU1VZJUWVmpyy+/3C/kZmdnKy8vT5K3t/iee+6pc8xHHnlEkyZN0k8//aTU1FQFBwcrJiZGMTExiouL05133qlf//rX2rFjh5555hnFxcUpNjbWee3bt6/atGmjqqoqWWsVFMSfMgBoifjrDuCgOnfuXG/N4/FoypQp9daTkpJUWlrqF2QLCwvVqVMnSd7e2kceecSpFRQUKD8/3+lh3b59u5566innxjifWbNm6corr9QXX3yh4cOHKyoqyi/IPv744xo8eLDWrVunN9980y/oxsXF6eSTT1ZUVJSqq6tljHGGYQAA3IWgCuCoCgkJUVJS0gEf3pCQkFBnRoXaBgwYoLKyMpWUlOjnn39Wfn6+8vPz1b17d0lShw4d9OCDDyo/P9+v7uthXbt2rTO7Qm1Lly7VoEGDNGPGDE2YMKFOkH3xxRfVsWNHLVu2TF988YWz3rdNWlqaPB5PE/2GAAD1YYwqgBatsrLSCbG+18GDBys6OlrLly/X//73Pyfg+urvvPOO2rdvr0cfffSAQxd27dqlpKQk3X///Xr22WfrBNlZs2YpJCREixYt0oYNG/x6e+Pi4tSjR49m+E0AgDtxMxUAHAFrrYqKiuoE3VGjRikoKEjvv/++5s+f79ebW1hYqDVr1sgYo3Hjxjk3pflERkY6U5ONGzdO8+fPd4JufHy8OnXqpGeffVaSNG/ePOXn5/vVExISlJCQcMx/FwBwtBBUAaAZlJeXOwHWF3RLS0v1q1/9SpL0j3/8Q1999ZXfNrGxsfrss88kSWeffbY++eQTv2P26dNHa9askSRdcskl+uGHHxQfH++E2b59++q2226TJM2fP1/GGKfm69Xl5jMAbkJQBYDjUG5urnbv3q38/Hzl5eUpPz9fYWFhuuSSSyRJkyZNUmZmpl/QTU9P1zvvvCNJSk1N1ZYtW/yOOXr0aKd+7rnnqrq62umtjYuL09ChQzV69GhJ3qemRUdHOyE3JiaGJ6MBaHJMTwUAx6FDXeZ/9NFHD7r/woUL6wRd34wLkhQTE6OdO3dq3bp1TtAtKirS6NGjVVVVpXPOOcfveMYY3XnnnXrsscdUUlKiX/7yl369tfHx8Ro+fLgGDhyo8vJyZWRkOOujo6OZXQHAYSOoAkAL1b17d2eGhAP573//W2ddVVWV8/6LL75wAq5vGTJkiCSppKRExcXF2r59uxOEKyoq9Pe//10DBw7U1q1b/ebwDQwMVGxsrJ588kmNHTtWW7Zs0T333OP05MbHxys+Pl7Dhg1Tx44dnWnN4uLimGEBaMUIqgAAh+/SfmBgoE4//fR6t4uPj9dXX33lfLbWqri42Ok1bdOmjd5++20n4Obl5SkvL08nnniiJCkvL09LlixRXl6eCgoK5BuG9uabb6pjx4764osvdN5550mSoqKinDA7depUDRkyRGvWrNGsWbOcgOurp6enM0cu0IIQVAEAjWaM8XsiWlRUlC6++OJ6t+/fv79++OEHSd5e3IKCAuXl5Sk5OVmSdNJJJ2nKlClOwPWF3cjISEnSunXrDvgwiFWrVumUU07RP/7xD916661+vbVxcXF66aWX1K5dOy1fvlwrVqxwar4lNTWVcbiAi3AzFQDguOTrxa0dZtPT0xUZGamlS5fqnXfeqRN058+frzZt2uiBBx7Qgw8+WOeYvpkXHn74Yc2YMaNOkH3uuecUGBio5cuXKzs7u07d91Q1AA3HXf8AANRSUVHhNyTBt1x99dUKCAjQ7Nmz9d577/nViouLtW3bNknS2LFjNXPmTL9jxsbGKj8/X5J05513Oj22CQkJzhy5N910kyQpMzNT1lon4IaGhh7bXwDgIgRVAACa0I4dO5SVleUXZKuqqnTLLbdIku677z598sknys3NderdunVTZmamJOmMM87Q4sWLneOFh4frzDPP1Lx58yR5g25+fr5fb223bt00bNgwSVJOTo6ioqIIuGgRmJ4KAIAmdMIJJ+iEE06ot/7QQw/5fbbWqrS01Pk8efJkbd26VXl5eU6YbdOmjVP/5ptv9M033yg3N9cZh/uLX/zCCar9+vVTVlaWwsLCnF7biy++2BnO8OCDD/rV4uPj1blzZ3Xs2LHJfgfAsUBQBQDgKDPGKCwszPk8aNAgDRo0qN7t58+fL2nfONzc3Fy/+sMPP6ydO3f6BV3fjWbV1dX629/+5heMJenWW2/Vc889p7Kysjpja+Pj43XllVfq0ksvVWlpqV599VUn4PqWpKQkBQcHN9WvBGgQgioAAC7lm02h9owKknT99dfXu09AQICKi4udgOsbeuDrAa6qqtLNN9/sV9uwYYN27twpyTus4Xe/+12d4z711FP6wx/+oB9++EG/+tWvnCDrex0zZoz69eungoICrVmzxi/oEnBxpAiqAAC0MLUD7v6X+8PDw/XEE0/Uu29KSoq2bNniN/42NzdXp512miRvL2+XLl2Ul5en9evXO/X+/furX79+WrVqVZ2nmkVFRenNN9/U+eefrxUrVmjy5MnOk9d8YXfkyJFKTk7W3r17VVpaqtjYWKYKA0EVAADsExQUpI4dO9Y7nrVr16565513/NZZa52HNvTv318LFy50Aqyv5zY1NVWSlJubq9WrVys3N1f5+fmqrq6WJC1btkzJycmaPXu2xo0bJ2OMYmNjnTD7+uuvKzU1VV988YU+/vhjv5CbkJCgvn378hSzFoigCgAAGqX2U8Di4uI0YsSIercdOXKkNmzYIMk7nragoEC5ublKSUmRJA0ePFjPPvusX8jNzc11xvh++eWXeuCBB+ocNycnR4mJiXrggQc0derUOkMTpk2bpuDgYC1ZskTbtm3zC7nx8fF1hlfAHZieCgAAHFcqKyuVn5/vF2YvvPBCBQYGas6cOZo3b55fyC0sLNTGjRtljNF1112nGTNm+B0vKipKhYWFkqS77rpLS5Ys8QuxHTt2dKYeW7dunaqrq50a428bj3lUAQAAJO3evVs7duzwC7nl5eW6+eabJXlnVPjoo4+cem5urjp37uz0Ap911ln67LPPnONFRUVp6NChzhy49913n/Ly8vx6dE888UQNGTJEklRUVKSIiAgFBAQc45/cvQiqAAAAR8Baq7KyMufhCkuXLtWWLVv8gmxycrImTZokSTrvvPO0atUq5efnO+N2R44c6QTZjh07KisrS3FxcU6v7YUXXqh77rlHkvTMM8/4zYGbkJCg9u3bKzExsRl++mODCf8BAACOgDHG7wlggwcP1uDBg+vdfuHChZK804D9/PPPys3N9es9vfvuu5WVleUXdKuqqiR5Q/Edd9yhyspKv2NOmDBBU6dOVWVlpTp37uw39jYhIUEXXXSRfvGLX6i8vFzz5893Am5CQoLi4uIUFHT8xr3jt+UAAAAuFRgY6ITF2iZMmFDvPsYYFRUV+c2YkJubq06dOkmSysvLNWLECGfIwrfffusMTfjFL36hHTt2aPTo0XWO++yzz+q2227T1q1bdeONN9aZMeG8885T9+7dm/YX0EQIqgAAAC4RGhqqdu3aqV27dnVq4eHhmj59er37tmnTRitWrPALuXl5eRo4cKAkqbi4WLt379Z3332n3NxcFRQUSJL+85//HL9B1RgTKulzSSE1279lrb3fGBMv6XVJqZI2S7rcWpt/9JoKAACA+oSEhCg9/YBDPSVJPXr00IoVK5zPFRUVys/Pd/XUXA255axM0tnW2r6S+kkaaYwZLOkuSYustd0kLar5DAAAgOOAx+NRcnLy8R1Urdeemo+emsVKGi3plZr1r0i6+Ki0EAAAAK1SgybxMsYEGmNWS9ol6UNr7TJJbay12ZJU85pcz77jjTErjTErc3JymqrdAAAAaOEaFFSttVXW2n6SUiQNNMakNfQE1tpp1tp0a216UlLSkbYTAAAArcxhPRbBWvuzpE8ljZS00xjTVpJqXnc1eesAAADQah0yqBpjkowxsTXvwySdK2m9pDmSxtZsNlbSu0erkQAAAGh9GjKPaltJrxhjAuUNtm9Ya98zxiyR9IYx5gZJWyVddhTbCQAAgFbmkEHVWrtGUv8DrM+VdM7RaBQAAABwWGNUAQAAgGOFoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcCWCKgAAAFyJoAoAAABXIqgCAADAlQiqAAAAcKVDBlVjTAdjzCfGmHXGmExjzO016+ONMR8aY76veY07+s0FAABAa9GQHtVKSf9nre0pabCkW4wxvSTdJWmRtbabpEU1nwEAAIAmccigaq3Nttb+v5r3RZLWSWovabSkV2o2e0XSxUerkQAAAGh9DmuMqjEmVVJ/ScsktbHWZkveMCspuZ59xhtjVhpjVubk5DSutQAAAGg1GhxUjTGRkv4r6ffW2sKG7metnWatTbfWpiclJR1JGwEAANAKNSioGmM88obU/1hr/1ezeqcxpm1Nva2kXUeniQAAAGiNGnLXv5H0T0nrrLVP1SrNkTS25v1YSe82ffMAAADQWgU1YJuhkq6RtNYYs7pm3SRJf5P0hjHmBklbJV12dJoIAACA1uiQQdVau1iSqad8TtM2BwAAAPDiyVQAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVDhlUjTHTjTG7jDEZtdbFG2M+NMZ8X/Mad3SbCQAAgNamIT2qMySN3G/dXZIWWWu7SVpU8xkAAABoMocMqtbazyXl7bd6tKRXat6/IuniJm4XAAAAWrkjHaPaxlqbLUk1r8n1bWiMGW+MWWmMWZmTk3OEpwMAAEBrc9RvprLWTrPWpltr05OSko726QAAANBCHGlQ3WmMaStJNa+7mq5JAAAAwJEH1TmSxta8Hyvp3aZpDgAAAODVkOmpXpO0RNJJxphtxpgbJP1N0ghjzPeSRtR8BgAAAJpM0KE2sNZeWU/pnCZuCwAAAODgyVQAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVCKoAAABwJYIqAAAAXImgCgAAAFciqAIAAMCVGhVUjTEjjTEbjDEbjTF3NVWjAAAAgCMOqsaYQEkvSLpAUi9JVxpjejVVwwAAANC6BTVi34GSNlprf5QkY8xsSaMlfVvfDhs2bNBZZ53ViFMCAACgtWjMpf/2kn6q9XlbzTo/xpjxxpiVxpiVFRUVjTgdAAAAWpPG9KiaA6yzdVZYO03SNElKT0+3n376aSNOCQAAgJbEmANFSq/G9Khuk9Sh1ucUSVmNOB4AAADgaExQXSGpmzGmszEmWNIYSXOaplkAAABo7Y740r+1ttIYM1HSAkmBkqZbazObrGUAAABo1RozRlXW2g8kfdBEbQEAAAAcPJkKAAAArkRQBQAAgCsRVAEAAOBKBFUAAAC4EkEVAAAArkRQBQAAgCsRVAEAAOBKxlp77E5mTI6kLcfshF6JknYf43Pi6ON7bZn4XlsmvteWie+15Wmu77STtTbpQIVjGlSbgzFmpbU2vbnbgabF99oy8b22THyvLRPfa8vjxu+US/8AAABwJYIqAAAAXKk1BNVpzd0AHBV8ry0T32vLxPfaMvG9tjyu+05b/BhVAAAAHJ9aQ48qAAAAjkMEVQAAALhSiw6qxpiRxpgNxpiNxpi7mrs9aDxjTAdjzCfGmHXGmExjzO3N3SY0DWNMoDHma2PMe83dFjQNY0ysMeYtY8z6mn+zQ5q7TWg8Y8wfav7+ZhhjXjPGhDZ3m3D4jDHTjTG7jDEZtdbFG2M+NMZ8X/Ma15xtlFpwUDXGBEp6QdIFknpJutIY06t5W4UmUCnp/+z/b+9+Qq2oAyiOfw+9glRaSZHvCRqICkEYEQ+FkHQhKD43gQtFxKWWrULduG0hois3agqJIiroQlSwhTsJLYh0Iyp67fkHxAo3Jh4XM4HIA+HN4G/ecD6bO/NbHRjm3nPn95sZeyEwCmzOce2NrcD10iGiVXuBc7YXAJ+R4zvlSRoGvgO+sP0p8A6wtmyqmKRDwIrXxrYBF23PAy7W+0X1tqgCXwI3bN+0/Qw4BowVzhQN2R63fbXe/pfqh2+4bKpoStIIsBLYXzpLtEPSB8BXwAEA289sPymbKloyBLwvaQiYBvxVOE9Mgu1LwOPXhseAw/X2YWDNWw01gT4X1WHg7iv7A1JoekXSHGARcLlskmjBHuAH4EXpINGaT4BHwE/1ko79kqaXDhXN2L4H7ALuAOPA37YvlE0VLfrI9jhUF4aADwvn6XVR1QRjeRZXT0iaAZwEvrf9T+k8MXmSVgEPbV8pnSVaNQR8DuyzvQh4SgemEaOZes3iGDAXmAVMl7SubKrosz4X1QEw+5X9ETI90QuS3qUqqUdsnyqdJxpbAqyWdJtqic7Xkn4uGylaMAAGtv+f8ThBVVxjalsO3LL9yPZ/wClgceFM0Z4Hkj4GqD8fFs7T66L6KzBP0lxJ71Et9j5TOFM0JElUa96u295dOk80Z3u77RHbc6jO019s5wrNFGf7PnBX0vx6aBlwrWCkaMcdYFTStPr7eBm5Sa5PzgAb6u0NwOmCWYBqaqaXbD+XtAU4T3VX4kHbfxaOFc0tAdYDf0j6vR7bYftswUwRMbFvgSP1xYKbwMbCeaIh25clnQCuUj2F5Tc6+NrNeDNJR4GlwExJA2An8CNwXNImqj8l35RLWMkrVCMiIiKik/o89R8RERERU1iKakRERER0UopqRERERHRSimpEREREdFKKakRERER0UopqRERERHRSimpEREREdNJLw4HI6BhnA2sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 脉冲响应分析\n",
    "irf = res.irf()\n",
    "\n",
    "fig = irf.plot(orth=True, impulse='sz')\n",
    "print('对上证指数的脉冲示意图')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matlab 实现\n",
    "\n",
    "### VAR 模型\n",
    "\n",
    "VAR.m\n",
    "```matlab\n",
    "function [beta,resid,cov_mat,AIC] = VAR(Y,X,p)\n",
    "% 注意:只是实现了两变量的VAR模型，思考如何实现任意个变量的VAR\n",
    "% p - 滞后阶数\n",
    "% AIC - 信息准则\n",
    "\n",
    "% 1.估计系数\n",
    "[ADLy,ADLx] = ADLxx(Y,X,p,p);\n",
    "[beta1,~,resid1] = regress(ADLy,ADLx);\n",
    "\n",
    "[ADLy,ADLx] = ADLxx(X,Y,p,p);\n",
    "[beta2,~,resid2] = regress(ADLy,ADLx);\n",
    "\n",
    "beta = [beta1';beta2'];\n",
    "resid = [resid1,resid2];\n",
    "\n",
    "% 2.估计AIC值(eviews调整自由度之后的)\n",
    "T = length(ADLy);\n",
    "cov_mat = resid'*resid/(T-3);\n",
    "AIC = log(det(cov_mat))+2*(2*p+1)*2/T;\n",
    "\n",
    "% % 1.载入数据\n",
    "% data = xlsread('C:\\Users\\Administrator\\Desktop\\hourse.xlsx');\n",
    "% f1 = data(:,2); f2 = data(:,3); e = data(:,6);\n",
    "% \n",
    "% % 2.估计VAR模型\n",
    "% [beta,resid,AIC] = VAR(f1,e,2);\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 脉冲响应函数\n",
    "\n",
    "OIRF1.m\n",
    "```matlab\n",
    "function [beta,IR] = OIRF1(Y,X,num,IMP)\n",
    "\n",
    "% num - 脉冲期数\n",
    "\n",
    "% 1.估计VAR模型参数\n",
    "[beta,~,cov_mat] = VAR(Y,X,1);\n",
    "\n",
    "% 2.脉冲响应函数\n",
    "% 2.1 正交化分解，估计P矩阵\n",
    "P = chol(cov_mat, 'lower');\n",
    "% 2.2 估计IR,s=1期为ADt,s=k，则为A^(k)Dt\n",
    "b = beta(:,2:3);\n",
    "SHOCK = zeros(2,1);\n",
    "if IMP == 1\n",
    "    SHOCK(1,1) = 1; \n",
    "elseif IMP == 2\n",
    "    SHOCK(2,1) = 1; \n",
    "end\n",
    "IR = zeros(num,2);\n",
    "IR(1,:) = b*(P*SHOCK);\n",
    "for s=2:num, IR(s,:) = (b*IR(s-1,:)')'; end\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}