{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 自回归分布滞后模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 原理讲解\n",
    "\n",
    "在自回归模型中,也可以引入其他变量来构成“自回归分布滞后模型”(Autoregressive Distributed Lag Model,简记 $ADL(p,q)$:\n",
    "\n",
    "$$y_{t}=\\beta_{0}+\\beta_{1} y_{t-1}+\\cdots+\\beta_{p} y_{t-p}+\\gamma_{1} x_{t-1}+\\cdots+\\gamma_{q} x_{t-q}+\\varepsilon_{t}\\tag{1}$$\n",
    "\n",
    "如果自回归分布滞后模型满足以下假定则万事大吉,可以用OLS来估计它。\n",
    "1. $\\mathrm{E}\\left(\\varepsilon_{t} | y_{t-1}, y_{t-2}, \\cdots, x_{1, t-1}, x_{1, t-2}, \\cdots, x_{K, t-1}, x_{K, t-2}, \\cdots\\right)=0$。这个假定类似于严格外生性假设,它意味着扰动项 $\\varepsilon_{t}$ 与所有解释变量的整个历史全部无关。这保证了对滞后期数 $\\left(p, q_{1}, \\cdots, q_{K}\\right)$ 的设定是正确的。如果滞后期数的设定不正确,比如,真实模型还应该包括 $y_{t-(p+1)}$,但该项 $\\beta_{p+1} y_{t-(p+1)}$ 却被纳入扰动项中,则扰动项 $\\varepsilon_{t}$ 便与解释变量相关,导致OLS不致。\n",
    "2. $\\left\\{y_{t}, x_{1 t}, \\cdots, x_{K t}\\right\\}$ 为渐近独立的平稳序列。\n",
    "3. $\\left\\{y_{t}, x_{1 t}, \\cdots, x_{K t}\\right\\}$ 有非零的有限四阶矩。\n",
    "4. 解释变量无完全多重共线性。\n",
    "\n",
    "对滞后期数的选择可以使用信息准则(最小化 AIC 或 BIC),或使用 $t,F$ 检验来检验最后期系数的显著性。更一般地,可以在 $ARMA$ 模型中引入其他变量,称为“$ARMAX$”模型。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.980</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.980</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   5540.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 29 May 2020</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>19:37:29</td>     <th>  Log-Likelihood:    </th> <td> -2414.2</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   458</td>      <th>  AIC:               </th> <td>   4838.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   453</td>      <th>  BIC:               </th> <td>   4859.</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     4</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td>   53.2565</td> <td>   26.908</td> <td>    1.979</td> <td> 0.048</td> <td>    0.376</td> <td>  106.137</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>y_1</th>   <td>    1.7194</td> <td>    0.202</td> <td>    8.511</td> <td> 0.000</td> <td>    1.322</td> <td>    2.116</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>y_2</th>   <td>   -0.7140</td> <td>    0.203</td> <td>   -3.517</td> <td> 0.000</td> <td>   -1.113</td> <td>   -0.315</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x_1</th>   <td>   -1.0193</td> <td>    0.274</td> <td>   -3.727</td> <td> 0.000</td> <td>   -1.557</td> <td>   -0.482</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x_2</th>   <td>    0.9939</td> <td>    0.275</td> <td>    3.620</td> <td> 0.000</td> <td>    0.454</td> <td>    1.534</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>28.071</td> <th>  Durbin-Watson:     </th> <td>   2.008</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  85.486</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.167</td> <th>  Prob(JB):          </th> <td>2.73e-19</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 5.090</td> <th>  Cond. No.          </th> <td>8.10e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 8.1e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.980\n",
       "Model:                            OLS   Adj. R-squared:                  0.980\n",
       "Method:                 Least Squares   F-statistic:                     5540.\n",
       "Date:                Fri, 29 May 2020   Prob (F-statistic):               0.00\n",
       "Time:                        19:37:29   Log-Likelihood:                -2414.2\n",
       "No. Observations:                 458   AIC:                             4838.\n",
       "Df Residuals:                     453   BIC:                             4859.\n",
       "Df Model:                           4                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const         53.2565     26.908      1.979      0.048       0.376     106.137\n",
       "y_1            1.7194      0.202      8.511      0.000       1.322       2.116\n",
       "y_2           -0.7140      0.203     -3.517      0.000      -1.113      -0.315\n",
       "x_1           -1.0193      0.274     -3.727      0.000      -1.557      -0.482\n",
       "x_2            0.9939      0.275      3.620      0.000       0.454       1.534\n",
       "==============================================================================\n",
       "Omnibus:                       28.071   Durbin-Watson:                   2.008\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               85.486\n",
       "Skew:                          -0.167   Prob(JB):                     2.73e-19\n",
       "Kurtosis:                       5.090   Cond. No.                     8.10e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 8.1e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "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",
    "ADLy, ADLx = lag_list(Y, X ,p=2, q=2)\n",
    "ADLx = sm.add_constant(ADLx)\n",
    "mod = sm.OLS(ADLy, ADLx)\n",
    "res = mod.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## statsmodels 库实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "未找到,官方介绍基本模型只有 AR、VAR、ARMA\n",
    "\n",
    "statsmodels.tsa contains model classes and functions that are useful for time series analysis. Basic models include univariate autoregressive models (AR), vector autoregressive models (VAR) and univariate autoregressive moving average models (ARMA). Non-linear models include Markov switching dynamic regression and autoregression. It also includes descriptive statistics for time series, for example autocorrelation, partial autocorrelation function and periodogram, as well as the corresponding theoretical properties of ARMA or related processes. It also includes methods to work with autoregressive and moving average lag-polynomials. Additionally, related statistical tests and some useful helper functions are available."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matlab 实现\n",
    "\n",
    "ADLxx.m\n",
    "```matlab\n",
    "function [ADLy,ADLx] = ADLxx(Y,X,p,q)\n",
    "\n",
    "% 待估计方程:y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q)\n",
    "% 获取自回归分布滞后模型的估计向量\n",
    "%  X   - 解释变量,列向量\n",
    "%  Y   - 被解释变量,列向量\n",
    "%  p   - ADL模型Y的滞后阶数,标量\n",
    "%  q   - ADL模型X的滞后阶数,标量\n",
    "\n",
    "T = length(Y);\n",
    "N = max(p,q)+1;\n",
    "\n",
    "ADLy = Y(N:T);\n",
    "c = ones(length(ADLy),1);  %常数项\n",
    "ADLx = zeros(length(ADLy),1+p+q);\n",
    "ADLx(:,1) = c;\n",
    "for i=1:p, ADLx(:,1+i) = Y(N-i:T-i);end\n",
    "for i=1:q, ADLx(:,1+p+i) = X(N-i:T-i);end\n",
    "\n",
    "% % 1.载入数据\n",
    "% data = xlsread('../数据/上证指数与沪深300.xlsx');\n",
    "% f1 = data(:,1); f2 = data(:,2)\n",
    "% \n",
    "% % 2.估计ADL模型f1 = c + f1(-1) + f1(-2) + e(-1)\n",
    "% [ADLy,ADLx] = ADLxx(f1,e,2,1);\n",
    "% beta = regress(ADLy,ADLx);\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
}