{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 分位数VAR模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 原理讲解\n",
    "\n",
    "参考文献：\n",
    "1. VAR for VaR: Measuring tail dependence using multivariate regression quantiles\n",
    "2. 分位数向量自回归分布滞后模型及脉冲响应分析"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用 statsmodels 库实现\n",
    "\n",
    "### QVAR 模型估计"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>hs300_1</th>\n",
       "      <th>hs300_2</th>\n",
       "      <th>sz_1</th>\n",
       "      <th>sz_2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>hs300</th>\n",
       "      <td>17.229937</td>\n",
       "      <td>1.390869</td>\n",
       "      <td>-0.369268</td>\n",
       "      <td>-0.589803</td>\n",
       "      <td>0.549461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sz</th>\n",
       "      <td>11.392195</td>\n",
       "      <td>0.345127</td>\n",
       "      <td>-0.324991</td>\n",
       "      <td>0.502069</td>\n",
       "      <td>0.463664</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           const   hs300_1   hs300_2      sz_1      sz_2\n",
       "hs300  17.229937  1.390869 -0.369268 -0.589803  0.549461\n",
       "sz     11.392195  0.345127 -0.324991  0.502069  0.463664"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import statsmodels.regression.quantile_regression as qr\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, Yname='', Xname='', exogenous=[]):\n",
    "    '''\n",
    "    待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q) + exogenous\n",
    "    获取自回归分布滞后模型的估计向量\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    Y : 被估计变量\n",
    "    X : 估计变量\n",
    "    p : ADL 模型 Y 的滞后阶数，标量默认为1\n",
    "    q : ADL 模型 X 的滞后阶数，标量默认为1\n",
    "    Yname : 被响应变量名\n",
    "    Xname : 响应变量名\n",
    "    exogenous : 外生变量\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    ADLy : ADL 模型被解释变量\n",
    "    ADLx : ADL 模型解释变量\n",
    "\n",
    "    '''\n",
    "    if not Yname:\n",
    "        Yname = 'y'\n",
    "    if not Xname:\n",
    "        Xname = 'x'\n",
    "        \n",
    "    ADLx = pd.DataFrame()\n",
    "    T = len(Y)\n",
    "    ADLy = pd.Series(list(Y[max(p, q):T]), name=Yname)\n",
    "    for i in range(1, p+1):\n",
    "        name = f'{Yname}_{i}'\n",
    "        ADLx[name] = list(Y[max(p, q)-i:T-i])\n",
    "    for i in range(1, q+1):\n",
    "        name = f'{Xname}_{i}'\n",
    "        ADLx[name] = list(X[max(p, q)-i:T-i])\n",
    "    \n",
    "    # 增加控制变量\n",
    "    if type(exogenous) == pd.Series:\n",
    "        ADLx[exogenous.name] = list(exogenous[:0-max(p, q)])\n",
    "    elif type(exogenous) == pd.DataFrame:\n",
    "        for name in exogenous.columns:\n",
    "            ADLx[name] = list(exogenous[name][:0-max(p, q)])\n",
    "    \n",
    "    # 增加常数项\n",
    "    ADLx = sm.add_constant(ADLx)\n",
    "    return ADLy, ADLx\n",
    "\n",
    "\n",
    "def qvar(Y, X, P, Q, Yname='', Xname='', exogenous=[]):\n",
    "    '''\n",
    "    待估计方程：y = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q) + exogenous\n",
    "               x = c + y(-1) +....+y(-p) + x(-1) + ... + x(-q) + exogenous\n",
    "    估计 QVAR 模型\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    Y : 被估计变量\n",
    "    X : 估计变量\n",
    "    P : QVAR 模型的滞后阶数\n",
    "    Q : 分位点\n",
    "    Yname : 被响应变量名\n",
    "    Xname : 响应变量名\n",
    "    exogenous : 外生变量\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    res1, res2 : 两模型估计结果\n",
    "    \n",
    "    '''\n",
    "    ADLy, ADLx = lag_list(Y, X ,P, P, Yname, Xname, exogenous)\n",
    "    mod = qr.QuantReg(ADLy, ADLx)\n",
    "    res1 = mod.fit(Q)\n",
    "    \n",
    "    ADLy, ADLx = lag_list(X, Y ,P, P, Xname, Yname, exogenous)\n",
    "    mod = qr.QuantReg(ADLy, ADLx)\n",
    "    res2 = mod.fit(Q)\n",
    "    \n",
    "    return res1, res2\n",
    "\n",
    "res1, res2 = qvar(Y, X, 2, 0.3, 'hs300', 'sz')\n",
    "beta1 = res1.params\n",
    "beta2 = res2.params\n",
    "pd.DataFrame([beta1, beta2], ['hs300', 'sz'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 脉冲响应分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "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>0</th>\n",
       "      <th>1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.224023e+06</td>\n",
       "      <td>868179.613056</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8.681796e+05</td>\n",
       "      <td>645978.438493</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              0              1\n",
       "0  1.224023e+06  868179.613056\n",
       "1  8.681796e+05  645978.438493"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def OIRF(res1, res2, p):\n",
    "    '''\n",
    "    估计脉冲响应函数\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    res1, res2 : 两模型估计结果\n",
    "    p : 滞后阶数\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    OIRF\n",
    "    \n",
    "    '''\n",
    "    pass\n",
    "p = 2\n",
    "resid = pd.DataFrame([res1.resid, res2.resid]).T   # 残差序列\n",
    "# 正交化分解，估计P2矩阵\n",
    "a = resid.T @ resid\n",
    "a=a/(len()-2*p-1);\n",
    "P2 = chol(a, 'lower');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matlab实现\n",
    "可以参考 matlab 代码：[QVAR 模型](https://github.com/lei940324/econometrics/tree/master/matlab代码/分位数回归/VAR模型)"
   ]
  }
 ],
 "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
}